1. 1.5 degree finite difference method wave equation migration
Wave equation migration with 15-degree finite difference method is a migration method, which takes the horizontal superimposed time profile obtained on the ground as the boundary condition, uses difference instead of differential to solve the approximate wave equation containing only upward waves, obtains the wave field values of underground points, and then obtains the real image of underground interface. Migration process is also a continuation and imaging process.
Derivation of 1) continuation equation
Starting with the following two-dimensional wave equation
Seismic wave field and seismic exploration
According to the explosion reflection surface model, the speed is halved, that is, v/2 is used instead of V, and the following results can be obtained:
Seismic wave field and seismic exploration
This equation has two solutions, corresponding to the upward wave and the downward wave respectively. Seismic record is a simple upgoing wave record, which can not be extended by this equation, but must be converted into a simple upgoing wave equation before it can be used. The usual method is approximate after coordinate transformation. The first step is coordinate transformation, so that
Seismic wave field and seismic exploration
There is no change in the first transformation of the above formula; The second transformation only changed the space depth z into the time depth τ, and there was no substantial change. The key is the third transformation, that is to say, instead of using the traditional old clock to measure time, we use a new clock with the same speed as the old clock, but with a different depth at the beginning. When the new clock is used, the upper and lower waves show the difference.
Because coordinate transformation will not change the actual wave field, the wave field u(x, z, t) in the original coordinate system is exactly the same as that in the new coordinate system, that is
Seismic wave field and seismic exploration
Through the differential method of composite function, we get:
Seismic wave field and seismic exploration
Substitute the above second-order partial differential result into equation (4-4-3), and after completion, we get:
Seismic wave field and seismic exploration
For writing convenience, if u, x and t are replaced respectively, the formula (4-4-5) can be written as follows.
Seismic wave field and seismic exploration
Where uxx, uπτ and uτt respectively represent the second derivative of U. attention. This equation still contains upward waves and downward waves, which cannot be used for continuation, so there is a second step.
After coordinate transformation, although the wave field remains unchanged, the upgoing wave and downgoing wave show differences in the new coordinate system, mainly in the different sizes of uπτ: when the included angle between the propagation direction of upgoing wave and the vertical direction is small (less than 15), uπτ can be ignored; For downgoing waves, u▽ can't be ignored. Ignoring the uπτ term, an approximate equation containing only upward waves is obtained.
Seismic wave field and seismic exploration
This is an approximate equation of 15 upward wave (because it is only applicable to upward waves whose included angle between the running direction and the vertical direction is less than 15, or only the upward reflection wave formed by the interface whose dip angle is less than 15 can be satisfied), which is a common continuation equation.
In order to solve this equation, we must also give the definite solution conditions. Due to the finite source intensity, the following explicit solution conditions can be given:
A the wave field outside the two ends of the survey line is zero, that is
U(x, τ, t)≡0 when x > xmax or x < xmin.
B the wave field outside the maximum recording time is zero, i.e.
When t > tmax, U(x, τ, t)≡0.
C. self-excited and self-received records (horizontal superimposed profile) with given boundary conditions, that is, the wavefield value u(x, 0, t) when the time depth τ = 0 is known.
Figure 4-4-5 12 Point Difference Scheme
With these definite conditions, equation (4-4-7) can be solved to obtain the wave field value u(x, τ, t) at any depth underground, which is a continuation process. Then, according to the imaging principle mentioned above, the zero-time wavefield value of the traditional old clock, that is, the wavefield value u(x, τ, τ) when the new clock time t = τ, constitutes the shifted output profile.
2) Establishment of difference equation
In order to solve the differential equation (4-4-7), approximate differential by difference, and adopt the 12 point difference scheme as shown in Figure 4-4-5, we can get:
Seismic wave field and seismic exploration
Substituting (4-4-8) and (4-4-9) into (4-4-7) gives:
Seismic wave field and seismic exploration
Seismic wave field and seismic exploration
Define vectors i, t:
I=[0， 1，0] T=[- 1，2，- 1]
Let the vector u (x, j, l) be
u(x，j，l)=[u(i- 1，j，l)，u(i，j，l)，u(i+ 1，j，l)]
Then (4-4- 10) can be abbreviated as
Seismic wave field and seismic exploration
And make
Seismic wave field and seismic exploration
Then (4-4- 1 1) can be written as follows:
[I-(α+β)T]u(x，j+ 1，l+ 1)-[I+(α-β)T]u(x，j，l+ 1)+[I-(α+β)T]u(x，j，l)=[I+(α-β)T]u(x，j+ 1，l)
Therefore, there are:
Seismic wave field and seismic exploration
This is a difference equation suitable for computer calculation.
3) calculation steps and migration results
The difference equation (4-4- 12) is an implicit equation in form, that is, the wave field value at time depth τ = (j+ 1) τ cannot be obtained by combining the wave field values at time depth τ = j τ, and there is another term τ = (j+ 1) τ on the right side of the equation. As shown in figure 4-4-6, in order to obtain a row of data u(x, j+ 1, l), three rows of data u(x, j+ 1, L+ 1) and u(x, j, L+ 1) must be used. Generally speaking, the implicit equation must be solved by solving simultaneous equations, which is troublesome, but favorable conditions for definite solution can be used here without complicated simultaneous operations.
Using the definite solution condition B, when calculating the wave field value at the new depth τ = (j+ 1) τ, the value line of t = tmax is calculated from the maximum time. Because u(x, j+ 1, tmax+δ t) ≡ 0 and u(x, j, tmax+δ t) ≡ 0, there are:
Seismic wave field and seismic exploration
It is very easy to calculate u(x, j+ 1, tmax) only by using the known values of u(x, j, tmax). Then using the formula (4-4- 12), it is not difficult to recursively calculate the wave field value at any time at the depth of τ = (j+1).
In concrete calculation, the wave field value at the depth Δ τ is calculated by extending downward from the ground: firstly, the wave field at this depth is calculated at t = tmax, and then it is calculated in the direction of decreasing t until all the wave field values at this depth are calculated. After calculating the wave field value of depth, the calculation is continued by extending downward by a step Δ τ. By analogy, the wave field values of all points under the ground at different times can be obtained.
As mentioned above, the wavefield value at the new clock t = τ is the expected "image". Therefore, every time the wave field value of a certain depth τ is recursively calculated, the calculation from t = tmax to t = τ can be completed, and u (x, τ, τ) is the "image" of this depth. The "images" with different depths constitute the migrated output profile.
Figure 4-4-6 One-step finite difference migration solution
①u(x，j，l+ 1)，②u(x，j，l)，③u(x，j+ 1，l+ 1)，④u(x，j+ 1，l)
Figure 4-4-7 Location Map of Migration Results
Figure 4-4-7 shows the calculation relationship and the location of the results when the migration occurs. A represents the superimposed profile observed on the ground, and the wave field value b of the next depth Δ is calculated by A. When calculating B, the value of 1' line is calculated first (only the value of 1 line in A is used), and then the value of 2' line is calculated (using the values of 1 and 2 lines in A and 1' line in B).
When the continuation calculation step Δ τ is the same as the sampling interval Δ t of seismic records, it can be seen from the geometric relationship in Figure 4-4-7 that the migration profile is the value on the diagonal of 45 in the figure. In practical work, Δ τ is not necessarily equal to Δ t, but should be determined according to the inclination angle of the interface. When the inclination angle is large, it should be smaller, and when the inclination angle is small, it can be larger to reduce the calculation workload. The intermediate value is obtained by interpolation.
Compared with other wave equation migration methods, the finite difference method has the advantages of adapting to the lateral velocity change, low migration noise and working well under the condition of low profile signal-to-noise ratio, but the 15 finite difference method can not get good migration effect when the interface dip angle is too large. Therefore, the finite difference migration methods of 45, 60 or even 90 have been developed, and interested readers can refer to relevant literature.
2. Wave equation migration in frequency wave number domain
The finite difference migration method is used to calculate in time domain and space domain. Fourier transform can also be used to realize frequency wave number domain migration.
Just like the idea of finite difference migration, it is considered that the horizontal stack profile is the wave field value u(x, 0, t) of the upward waves emitted by countless sources at the same time on the interface, and it is a continuation process to use it to invert the wave field value u(x, z, t) at any point underground. Then, according to the imaging principle, take the value u(x, z, 0) when t = 0 to form the offset output profile.
Starting from the wave equation (4-4-3) after the speed is halved, two-dimensional Fourier transform is carried out on both sides of the equation on X and T, and an ordinary differential equation is obtained:
Seismic wave field and seismic exploration
Where: u = u (kx, z, ω) is the two-dimensional Fourier transform of the wavefield function u(x, z, t), ω = 2π f is the angular frequency, and kx is the spatial wavenumber in the X direction.
Equation (4-4- 13) is an ordinary differential equation, which is easy to solve. There are two solutions, corresponding to upward wave and downward wave respectively. Migration studies the downward continuation of upward wave, so only the upward wave solution is considered.
Seismic wave field and seismic exploration
Where U(kx, 0, ω) is the initial value of the solution, that is, the Fourier transform of the ground (z = 0) uplink wave record. Therefore, equation (4-4- 14) represents the process of obtaining the Fourier transform of the wave field at any depth underground from the Fourier transform of the wave field at z = 0, which is the continuation of the wave field in the frequency wavenumber domain.
The wave field value of any depth underground can be obtained from U(kx, z, ω) by inverse Fourier transform:
Seismic wave field and seismic exploration
According to the imaging principle, the migration result should be the wave field value at this depth t = 0:
Seismic wave field and seismic exploration
This is the mathematical model of frequency wave number domain migration. Specific implementation steps will not be described in detail.
When solving the ordinary differential equation (4-4- 13), the initial value does not take the Fourier transform of the wave field value at z = 0, but takes the Fourier transform value of the wave field at any shallow place, so that:
Seismic wave field and seismic exploration
Thus, the mathematical model of phase shift offset is obtained:
Seismic wave field and seismic exploration
Using this formula, imaging can be expanded downward step by step, and the speed used for each expansion can be changed, so the phase shift method can adapt to the longitudinal change of speed.
Because of the application of fast Fourier transform, the frequency wavenumber domain migration method has high efficiency and less running time. This is the most economical method in wave equation migration algorithm, which is suitable for large dip areas. Because the calculation is carried out in the frequency wavenumber domain, we should pay attention to the problem of false frequency, and this method is not suitable for areas where the lateral velocity changes.
3. Shikhov integral migration.
G-Shikhov integral migration is a migration method based on G-Shikhov integral solution of wave equation.
The Shikhov integral solution of the three-dimensional longitudinal wave equation (see chapter 1) is as follows
Seismic wave field and seismic exploration
Where q is a closed surface around the point (x, y, z), n is the outer normal of q, r is the distance from the point (x, y, z) to every point on the q plane, and [] represents the delay bit.
The essence of this solution is to calculate the wave field value of any point on the plane from the known wave field value of each point on the closed surface Q, which is the strict mathematical form of Huygens principle.
The selected closed surface Q consists of an infinite flat surface Q0 and an infinite hemispherical surface Q 1. Q 1 The contribution of the area of the wave field value of each point on the plane to the wave field function of one point on the plane is zero, so the wave field value of each point underground is only calculated by the wave field value of each point on the flat Q0. Under this condition, the wave field value of any point underground is
Seismic wave field and seismic exploration
At this point, the term in the original formula disappears, and the negative sign before the integral sign becomes positive because the positive direction of the z axis is opposite to N.
Knowing the source function, it is a positive problem to find the wave field value of the wave propagating to a certain point. The above is the formula for calculating the K-Shikhov integral of the forward problem. Migration processing is an inverse problem. The wave field value received on the ground is regarded as the secondary source, and the underground wave field value is "backward" in time, and the reflection interface is determined by the wave field value when t = 0. The inverse problem can also be solved by the above formula, the only difference is that [] is no longer a delay bit but an advance bit. According to this understanding, the continuation formula of K-Shikhov integral should be
Seismic wave field and seismic exploration
According to the imaging principle, the wave field value at t = 0 is the migration result. Only two-dimensional offset is considered, and y coordinate is ignored. The space depth z is converted into time depth t0 = 2z/v, and the Kr-Shikhov integral migration formula is obtained.
Seismic wave field and seismic exploration
Where:; Xl is the abscissa of the ground trajectory; X is the abscissa of the migration profile (see Figure 4-4-8).
Figure 4-4-8 Schematic diagram of each quantity in Shikhov shift formula.
From, from:
Seismic wave field and seismic exploration
It can be seen that K-Shikhov integral migration is very similar to diffraction scanning superposition, which is based on the value of diffraction hyperbola and placed at the vertex of hyperbola. The difference is:
A. Not only the amplitude of each trace, but also the derivative value of the amplitude of each trace to time should be taken to participate in superposition;
B. When stacking the amplitudes corresponding to each channel, it is not simple addition, but weighted superposition according to Formula (4-4-22).
Because of this, K-Shikhov integral migration method is similar to diffraction scanning superposition migration in form, but it is different in essence. The former is based on wave equation, which can preserve the dynamic characteristics of waves; The latter belongs to the category of geometric seismology, which only retains the kinematic characteristics of waves.
Compared with other wave equation migration methods, K-Shikhov integration method is easy to understand and can adapt to large dip strata. It is difficult to use in areas with large lateral speed change and large offset noise.
Seismic wave actually propagates in three-dimensional space, so it must be three-dimensional migration to achieve complete migration. At present, the three-dimensional migration method has made great progress, from two-step method to split method, and now it has realized one-step complete three-dimensional migration without approximation.
The above post-stack migration has a basic assumption, that is, the horizontal stack profile is self-excited and self-collected, which is not true when the underground interface is complex. In order to realize real stacking and migration, a prestack migration method is developed, which is carried out simultaneously to ensure the real stacking of * * * reflection points. However, completing migration and stacking tasks at one time is not conducive to obtaining velocity parameters. Therefore, a prestack partial migration method is developed, and only partial migration is carried out, so that the * * * central trace is integrated into a real * * * reflection trace set, and the superposition of * * * reflection points is ensured. Pre-stack partial migration method plus post-stack migration is equal to pre-stack migration. With * * * reflection point gathers, velocity analysis is easy.
At present, most migration is still time migration, that is, time profile is obtained after migration. Time migration does not consider the bending phenomenon when seismic waves pass through the interface. If this phenomenon is considered, the depth profile is obtained after migration, which is called depth migration. 3D prestack depth migration is the most advanced migration method at present.
Seismic waves are actually elastic waves. At present, all migration methods use acoustic wave equation, which is only a kind of acoustic wave migration. To realize real complete migration, elastic wave migration is needed. Due to the need of multi-wave seismic data, elastic wave migration is rarely used at present.