- Home
- Documents
*Iterative asymptotic inversion in the acoustic construction of the asymptotic inverse operators by...*

prev

next

out of 17

View

0Download

0

Embed Size (px)

GEOPHYSICS, VOL. 57. NO. 9 (SEPTEMBER 1992): P. 1138-l 154. I3 FIGS

Iterative asymptotic inversion in the acoustic approximation

Gilles Lambar&*, Jean Virieuxz, Raul Madariaga”, and Side Jin*

ABSTRACT

We propose an iterative method for the linearized prestack inversion of seismic profiles based on the asymptotic theory of wave propagation. For this pur- pose, we designed a very efficient technique for the downward continuation of an acoustic wavefield by ray methods. The different ray quantities required for the computation of the asymptotic inverse operator are estimated at each diffracting point where we want to recover the earth image. In the linearized inversion, we use the background velocity model obtained by velocity analysis. We determine the short wavelength components of the impedance distribution by linear- ized inversion of the seismograms observed at the surface of the model. Because the inverse operator is not exact, and because the source and station distri- bution is limited, the first iteration of our asymptotic inversion technique is not exact. We improve the images by an iterative procedure. Since the back-

ground velocity does not change between iterations, there is no need to retrace rays, and the same ray quantities are used in the iterations. For this reason our method is very fast and efficient. The results of the inversion demonstrate that iterations improve the spatial resolution of the model images since they mainly con- tribute to the increase in the short wavelength contents of the final image. A synthetic example with one-dimen- sional (I-D) velocity background illustrates the main features of the inversion method. An example with two-dimensional (2-D) heterogeneous background dem- onstrates our ability to handle multiple arrivals and a nearly perfect reconstruction of a flat horizon once the perturbations above it are known. Finally, we consider a seismic section taken from the Oseberg oil field in the North Sea off Norway. We show that the iterative asymptotic inversion is a reasonable and accurate alter- native to methods based on finite differences. We also demonstrate that we are able to handle an important amount of data with presently available computers.

INTRODUCTION the integral equation that relates the model to the observed

The construction of images of the subsurface of the earth by the inversion of seismic reflection profiles has been investigated by many authors. In a pioneering work, Cohen and Bleistein (1977) obtained an approximate solution of the inverse problem for almost vertical incidence using asymp- totic methods. In a major contribution to inverse theory, Beylkin (1985) showed how to use asymptotic ray theory to construct an inverse operator for the case of a single source point and a continuous distribution of receivers on the surface of the model. Several improvements of Beylkin’s original method were proposed in the literature (e.g., Miller et al., 1987; Bleistein, 1987b: Beylkin and Burridge, 1990). A common feature of these techniques is that the inverse operator was constructed by mathematical manipulation of

seismograms. These methods are very fast because they are based on asymptotic theory. The main drawback for the construction of the asymptotic inverse operators by Bey- Ikin’s method is the need to establish a one-to-one relation- ship between the observed seismograms and the earth model parameters.

A very different approach to inversion, based on the theory of optimization, was proposed in a series of papers by Tarantola and coworkers (see, e.g., Tarantola, 1984; Pica et al., 1990: Crase et al., 1990). In the optimization approach, the hope of building an explicit inverse operator is aban- doned and instead, an iterative method is developed to find the model that best fits observations within a certain error criterion. The optimization approach leads to algorithms that

Manuscript received by the Editor July 15, 1991; revised manuscript received January 14. 1992. *Laboratorie de Sismologie, university Paris 7, institutede Physique du Globe. 75252 Paris Cedex 05, France. SInsmut de Geodynamique, Universite de Nice-CNRS, Av. Albert Einstein, 06560 Valbonne. France. 0 1992 Society of Exploration Geophysicists. All rights reserved.

1138

Iterative Asymptotic Inversion 1139

are much slower than those of Beylkin, but are much more robust and can handle incomplete and redundant data sets without problems.

Jin et al. (1992) proposed an asymptotic inversion method where the main advantages of both approaches to inversion were exploited to construct a fast and robust inversion method. In this approach, the inverse problem was formu- lated using classical optimization theory [see Tarantola (1987) for a recent review] but the gradient was computed using the methods of asymptotic ray theory. In their final result, Jin et al. (1992) found an iterative algorithm based on Newton’s optimization method. This algorithm allows for the use of redundant and incomplete data sets, as well as for the limited frequency band of the sources used in exploration geophysics.

The inverse problem requires a very efficient algorithm for the solution of the forward problem. Most authors have adopted finite difference methods for forward modeling (see, e.g., Kolb et al., 1984; Mora, 1986; Pica et al., 1990, etc), but the computer time required by finite difference is simply too long. An alternative to finite difference is to use ray theory that is probably just as accurate as the full numerical solution in the frequency band used in vertical seismic profiling. It should be pointed out, however, that ray theory is more economical than finite differences only if ray tracing is done very efficiently. In this article, we propose an efficient downward ray-tracing strategy to compute the different parameters-traveltime, slowness vector, and geometrical spreading-required at each diffracting point by the asymp- totic inverse method. We discuss the iterative asymptotic inversion described by Jin et al. (1992) in the acoustic approximation for a single parameter. It turns out that, because in linearized inversion the reference velocity model does not change, ray tracing can be done once for all the iterations. The iterations are thus run at a minimum compu- tational cost.

In the following, we present first a synthetic example with one-dimensional (1-D) velocity background that is used to demonstrate the practical importance of iterations. A two- dimensional (2-D) background is adopted for the second synthetic example to demonstrate the use of our method in the presence of a more realistic velocity background. Fi- nally, the acoustic algorithm will be applied to a seismic profile recorded on the North Sea, off-shore from Norway. With this example we intend to demonstrate the efficiency of our two-dimensional ray-tracing method and the importance of iterations for improving the image in a practical example.

DEPTH CONTINUATION OF ACOUSTIC FIELDS BY RAY THEORY

The high-frequency approximation of the acoustic field is required for both the forward and inverse problem. Two time scales are involved: the time of propagation and the time of the source wavelet. Using ray theory, we are able to distinguish the separate effects of these two time scales in the forward and inverse problems.

When a source has been fired at the point r, with the temporal signature S(t), the pressure P recorded at the point r outside the source area satisfies the Helmholtz equation in the frequency domain:

w* - P(r, rJ, 0) + V*P(r, r,Y, w) = 0, c2(r)

(1)

where c is the velocity and V2 denotes the laplacian. In a 2-D medium, the high-frequency approximation to the solution of this equation is:

P(r, rJ. w) = S(o)As(r, r,5)e’we”~“’ &-

(2)

where 8 satisfies the eikonal equation (On)2 = r-* and A0 the transport equation (Cerveny et al., 1977). The typical tail associated with 2-D propagation arises from the term l/v-i”. The slowness vector p = Vt3, the gradient of 0, is perpendicular to the phase fronts 8 = constant. The rays, which are the orthogonal trajectories to the phase fronts, are consequently tangent to the slowness vector. Selecting a sampling parameter for phase fronts, as well as for rays, will precisely define the eikonal. We choose the sampling param- eter 7, related to the traveltime by dn = u2d7, where u* is the square of slowness (cerveny, 1987; Virieux et al., 1988; Farra et al., 1989). We always perform kinematic and dynamic or paraxial ray tracing simultaneously to be able to compute the geometrical spreading, as well as control the ray sampling inside the medium.

Starting from one point at the free surface, we sample the medium down to a specified horizon. If caustics are de- tected, the different ray branches (cerveny et al., 1977; Hanyga, 1988) are located on this horizon and an interpola- tion is constructed inside each branch. After this first step, ray tracing is continued to the next horizon for each branch. With the help of interpolation, the ray may leave the horizon from any position suitable for a good discretization of the following horizon. We avoid the redundant procedure of computing branches on this new horizon starting again from the initial point. In this way, we guarantee a uniform sampling of the whole medium and avoid oversampling around the initial point and undersampling far away from the starting point. On every horizon,, small ray branches may be suppressed if they are poorly sampled. A minimum of five rays are required for example to accept a ray branc