Seismic tomography inversion method

There is an inseparable relationship between seismic observation data and earth parameters. For example, the velocity of seismic waves and the attenuation properties of the earth are related to the measured travel time and amplitude of seismic waves, which can be obtained by line integration along the ray path. Similarly, the transmission loss of high-frequency electromagnetic waves is also related to the electrical properties of the earth. Therefore, geophysical tomography technology is to invert these integral relations, so as to obtain the estimation of wave velocity field v(x, y) or electromagnetic loss coefficient field α(x, y) in some spatial areas where the ray passes (Figure 4- 18).

Fig. 4- 18 block diagram of geophysical tomography inversion technology and method

This section introduces the principles of matrix inversion and Fourier transform, and other seismic tomography inversion methods can refer to related literature.

Matrix inversion

Many methods, such as generalized matrix inversion, damped least square method or linear programming, can be 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 in each cell. Taking the travel time equation of earthquake as an example, its line integral can be approximately expressed as

Introduction to solid geophysics

Where:? Sj is the distance that the ray passes through the j-th unit; Vj is the seismic wave velocity in the j-th cell, and the k-th ray actually passing through the cell is summed. In practical application, the travel time equation is usually established according to the velocity disturbance of the initial model, that is,

Introduction to solid geophysics

here

Introduction to solid geophysics

That is, the disturbance of the j th unit slowness. This is because the amount of specific data is not enough to uniquely determine all the parameter values in the cell network array. In other words, there are more unknowns than linear independent equations. In this case, when solving the equations, additional constraints must be given. In fact, the usual solution process is to find such a solution, that is, not only to meet the data, but also to "approach" the preset initial model (such as the values of all j cells).

The form of matrix equation is

Y=AX

Where: y is the propagation time disturbance value. Tk (k = 1，2，...，m)； X is the slowness disturbance value? N-dimensional column vector composed of pj (j = 1, 2, …, n); A yes? The matrix of s (kmax×jmax), where kmax is the total number of rays passing through the study area and jmax is the total number of pixels. A is a sparse matrix (containing a large number of zero elements), because any ray can only pass through a few cells in the whole research area.

To form a matrix, we need to know the path of the ray from the source to the receiving point, and the path depends on the wave velocity field, so the answer to this question is the solution of the required solution. In the actual inversion, the ray path will be traced through the initial (hypothetical) model, and the A matrix will be constructed from it. The equation Y=AX can be solved by some matrix inversion programs. In the inversion calculation, the velocity disturbance on the X axis is kept small by proper damping treatment, so the ray path deviation caused by it is assumed to be very small and can be ignored. Once a new modified velocity model is obtained, it can be used to track seismic rays and construct a new A matrix. This process may need to be repeated many times to greatly reduce the remaining time for observation and calculation of travel.

For successful tomographic inversion, the angular coverage of seismic rays must be as wide as possible, and in most geophysical applications, the location distribution of source and receiver points is seriously limited. This inevitably leads to confusion and blurring of the obtained image. Figure 4- 19 gives a rather extreme example. Through the simple velocity model shown in Figure 4- 19(a), 140 rays were tracked from the explosion point of 10 to the receiving point of14.

The study area is 8 units wide and 12 units long. A is the matrix of 140×96, and Figure 4- 19(b) is the image obtained by the generalized inverse matrix inversion method. Fig. 4- 19(c) is the travel time inversion image with only one explosion point. In this example, the unclear outline of the image obtained in the direction of the ray path is the inevitable result of extremely narrow angular coverage.

Fig. 4- 19 velocity field tomography obtained by generalized inverse matrix inversion method

(a) There is a square abnormal area in the velocity field, and the rest are homogeneous areas, and the velocity difference between them is10%; (b) inversion results; (c) Inversion results obtained by using only one explosion point to emit 14 rays.

The biggest disadvantage of matrix inversion method in tomography is that the calculation workload is too large. If the number of cells in the study area is n, the computation of matrix inversion is n3. All 10× 10 units will need 106 time unit, while 100× 100 unit will need10/2 time unit. If 1μs is needed for each operation, then the reciprocal of 100× 100 takes 277h, and the reciprocal of 200×200 units takes two years. In many practical problems, the model of 100× 100 unit is not unrealistic, so it is very important to find other calculation methods with faster operation speed.

(2) Fourier transform method

Fourier transform method is developed on the basis of projection slice theorem, which shows: "When the angle is θ, the one-dimensional Fourier transform of projection is the tangent of the two-dimensional Fourier transform of its original object along the same angle (Mersereau et al., 1974)." Figure 4-20(a) is the velocity field distribution. There is a high-speed area of a square in the uniform velocity field (the four sides of this square are slightly smoother in the figure for calculation reasons, but this does not hinder the observation of its main characteristics). The projection of this velocity field is the travel time of many parallel rays passing through the field at a fixed angle. When the ray is parallel to the X-axis or Y-axis, its propagation time is abnormally rectangular, while when the ray is incident at an angle of 45 with the X-axis or Y-axis, its projection propagation time is abnormally triangular. Figure 4-20(b) is the two-dimensional Fourier transform of Figure 4-20(a). According to the projection section theorem, the plane passing through Figure 4-20(b) is the one-dimensional Fourier transform of Figure 4-20(a) projected at the same angle. People familiar with one-dimensional Fourier transform can easily find that the sections along line A-A and line B-B in Figure 4-20(b) will be rectangular and triangular Fourier transforms respectively.

Figure 4-20 Schematic Diagram of Projection Slice Theorem

(b) is the two-dimensional Fourier transform of (a), which passes through the plane of (b).

For example, A-A or B-B is a one-dimensional Fourier transform projected at the same angle by warp (A).

Figure 4-2 1 Fourier transform parallel beam projection along the x direction

Let the parallel ray beam follow the X 1 direction, and its projection value is (Figure 4-2 1).

Introduction to solid geophysics

Where: p stands for propagation time or amplitude; F(X 1, X2) is the attribute of the medium to be inverted. The Fourier transform of PX 1(X2) is

Introduction to solid geophysics

But the two-dimensional Fourier transform of f(X 1, X2) is

Introduction to solid geophysics

Substitute formula (4-49) into formula (4-50) and compare with formula (4-5 1).

PX 1(k2)=F(k 1，k2)k 1

The above proves the slice projection theorem. In fact, it can be extended to the projection of any angle θ, as shown in Figure 4-22.

Figure 4-22 Projection Diagram of Slice Projection Theorem Extending to Any Angle

The above discussion shows that the fast back projection method can be established in principle. The basic idea is that for the one-dimensional Fourier transform of projection functions with different angles, a two-dimensional Fourier transform function can be constructed, and then the two-dimensional inverse Fourier transform is carried out on this function to obtain the velocity field distribution.

The advantage of this method is fast operation speed, and the data volume can be halved by using the symmetrical characteristics of Fourier spectrum. The disadvantage of this method is that the high frequency loss is serious in frequency domain interpolation, which will affect the imaging resolution and the interpolation calculation is time-consuming. In view of the incomplete projection of seismology, we can use the method of limited angle projection reconstruction or the method of fan filtering to avoid the missing area of observation point (Menke, 1985), which will reduce the resolution of the image without losing the basic shape. In addition, the linearity of the ray and the structural form of the source-receiver system should be strictly limited.

Other inversion methods of seismic tomography mainly include convolution method or filtered back projection method (Andersen et al.,1984; Liu Futian et al., 1989), algebraic reconstruction method (Gordon et al.,1970; Liu Futian, 1988) and joint inversion of waveform and travel time (Dziewonski et al., 1993). And will not be described further.