Inversion Methods for Seismic Stratigraphic Imaging

There is an integral relationship between seismic observations and Earth parameters. For example, the seismic wave velocity and the attenuation properties of the earth are related to the measured travel time and the seismic wave amplitude, and this relationship can be obtained by line integration along the ray path. Similarly, the transmission loss of high-frequency electromagnetic waves is related to the electrical nature of the earth. Thus, the geophysical technique of stratigraphy is to invert these integral relationships to obtain an estimate of the wave velocity field v (x, y) or the electromagnetic depletion coefficient field α (x, y) in some region of space through which the rays travel (Figure 4-18).

Figure 4-18 Block diagram of geophysical stratigraphic inversion techniques

This section introduces the principles of matrix inversion and Fourier transform methods, and other seismic stratigraphic imaging inversion methods can be found in the relevant literature.

(I) Matrix inversion

Multiple methods such as generalized matrix inversion, damped least squares, or linear programming are easily applied to this problem. The area to be studied is divided into a grid array of many cells, and it is assumed that v(x, y) or α(x, y) is constant within each cell. Taking the seismic traveltime equation as an example, its line integral can be approximated as

Introduction to Solid-State Geophysics

Where: ?sj is the distance of the ray through the jth cell; vj is the seismic wave velocity in the jth cell, and the summation is implemented for the kth ray that actually passes through the cell. In specific applications, the traveltime equations are usually built in accordance with the amount of velocity perturbation for solving the initial model, i.e.

Introduction to Solid State Geophysics

Here

Introduction to Solid State Geophysics

that is, the amount of perturbation for the slowness within the jth cell. This is treated in this way because the amount of specific information is not sufficient to uniquely determine all parameter values within the cell mesh array. In other words, the number of unknown quantities is greater than the number of linearly independent equations. In this case additional constraints must be imposed on the solution of the system of equations, and the usual solution procedure is actually to search for a solution that not only satisfies the data, but is also "close" to the preset initial model (e.g., for all values of cell j).

The matrix equation is of the form

Y=AX

Where: Y is the m-dimensional column vector of walk-time perturbations ?tk (k=1, 2, ..., m); X is the n-dimensional column vector of slowness perturbations ?Pj (j=1, 2, ..., n); and A is the matrix of ?s (kmax×jmax), with a matrix of ?s (kmax×jmax). jmax) matrix, where kmax is the total number of rays that penetrate the study area and jmax is the total number of cells. A is a relatively sparse matrix (containing a large number of zero elements) because any one ray can only penetrate a small number of cells in the entire study area.

To form the A matrix it is necessary to know the path of the ray from the source to the reception point, which in turn depends on the wave velocity field, so the answer to this problem is exactly the solution to be found. In the actual inversion, the ray paths will be traced through the initial (envisioned) model, and the A matrix will be constructed from this. The equation Y = AX can be solved by some matrix inversion procedure.The amount of velocity perturbation in X, which is kept small by appropriate damping in the inversion calculations, and the resulting deviation of the ray paths, is assumed to be very small and negligible, so that the seismic rays can be traced through the new, revised velocity model once it is available and the new A-matrix can be constructed. This process may need to be repeated several times to greatly reduce the residuals between the observed and calculated travel times.

For a successful stratigraphic inversion, the angular coverage of the seismic rays must be as wide as possible, and in most geophysical applications, the distribution of source-receiver locations is severely limited. This inevitably results in confusion and ambiguity in the resulting image. A rather extreme example is given in Figure 4-19. With the simple velocity model shown in Figure 4-19(a), 140 rays were traced from the locations of 10 blast points to the locations of 14 receiving points.

The study area is 8 cells wide and 12 cells long, A is a 140 x 96 matrix, and Figure 4-19(b) shows the image obtained by the generalized inverse matrix inversion method. Figure 4-19(c) shows the image obtained by inversion of the walk time using only one explosion point. For the resulting image the appearance of indistinct contours in the ray path direction is a natural consequence of the extremely narrow angular coverage in this case.

Figure 4-19 Velocity field obtained by generalized inverse matrix inversion method of tomography

(a) the velocity field in a square anomaly, the rest of the homogeneous area, the two velocity difference of 10%; (b) inversion results; (c) using only one explosion point to send out 14 rays from the inversion results

the matrix inversion method used in tomography, the biggest disadvantage is that the computational effort is too great. If the number of cells in the study area is n, then the operation of inverting the matrix is on the order of n3. All 10 × 10 cells would require 106 time units, while 100 × 100 cells would require 1012 time units. If each step of the operation takes 1 μs, then the inverse of 100 × 100 units will take 277 h, while 200 × 200 units will take two years. In many real-world problems, the 100 × 100-unit model is not unrealistic, so it is crucial to seek other computational methods that operate faster.

(ii) Fourier transform method

The Fourier transform method was developed based on the projection slicing theorem, which states that "the one-dimensional Fourier transform of a projection at an angle θ is the tangent value of the two-dimensional Fourier transform of its original object along the same angle (Mersereau et al., 1974). " Figure 4-20(a) shows the distribution of a velocity field. There is a square region of high velocity in a uniform velocity field (for computational reasons, the four sides of this square are slightly smoothed in the figure, but this does not prevent observation of its main features). The projection of this velocity field is the travel time of many parallel rays passing through the field at a fixed angle. When the rays are parallel to the X- or Y-axis, the timing anomaly is rectangular, whereas when the rays are incident at an angle of 45° to the X- or Y-axis, the projected timing anomaly is triangular. Figure 4-20(b) is the two-dimensional Fuchsian transformation of Figure 4-20(a). According to the projective tangent theorem, the plane through Figure 4-20(b) is the one-dimensional Fuchs transform of Figure 4-20(a) projected at the same angle. Those familiar with the one-dimensional Fuchsian transform will easily realize that the tangent planes along the A-A and B-B lines in Figure 4-20(b) will be the Fuchsian transforms of rectangles and triangles, respectively.

Figure 4-20 Illustration of the projected slice theorem

(b) is the two-dimensional Fourier transform of (a) across the plane of (b).

For example, A-A or B-B, is the one-dimensional Fuchsian transform through (a) projected at the same angle

Figure 4-21 Illustration of Fuchsian transform of a parallel beam of rays projected along the X-direction

Let the parallel beam of rays be along the direction of X1, and the value of its projection be (Figure 4-21)

Introduction to Solid State Geophysics

Where: P denotes the travel time or amplitude; f(X1, X2) is the medium property to be inverted. the Fuchsian transform of PX1 (X2) is

Introduction to Solid State Geophysics

But the two-dimensional Fuchsian transform of f(X1, X2)

Introduction to Solid State Geophysics

Substituting Equation (4-49 ) into Eq. (4-50) and comparing it with Eq. (4-51) yields

PX1(k2) = F(k1, k2)k1

The above proves the sliced projection theorem, and in fact it can be generalized to projections at any angle θ, as shown in Figure 4-22.

Figure 4-22 Schematic diagram of the generalization of the slice projection theorem to projections at any angle

The above discussion illustrates that, in principle, it is possible to set up a fast inverse projection method. The basic idea is as follows: for the one-dimensional Fourier transform of the projection function of each different angle, a two-dimensional Fourier transform function can be constructed, and then the two-dimensional inverse Fourier transform of this function can be used to obtain the velocity field distribution.

The advantage of this method is that the operation speed is fast, but also can use the symmetry of the Fourier spectrum to halve the amount of data. The disadvantage of this method is that when interpolation is done in the frequency domain, the high-frequency loss is more serious, which will affect the imaging resolution, and, the interpolation calculation is also more time-consuming. For the case of seismological incomplete projection, the approach of finite angle projection reconstruction can be used or the use of sector filters to avoid the missing area of the observation point (Menke, 1985), which in turn will reduce the resolution of the image, but without losing the basic form. In addition, the straightness of the rays and the structural morphology of the source-receiver system should be strictly limited.

Other inversion methods for seismic stratigraphic imaging mainly include plethysmography or filtered inverse projection (Andersen et al., 1984; Liu Futian et al., 1989), algebraic reconstruction (Gordon et al., 1970; Liu Futian, 1988), and combined waveform and traveltime inversion (Dziewonski et al., 1993). etc., will not be described further.