A method for ultra-relaxed iterative prestack migration imaging based on thin-bed approximation

Abstract

The application relates to a kind of ultra relaxation iteration prestack migration imaging methods based on thin plate approximation, first, the velocity model is divided into thin plates;In the first thin plate, the first-order ultra relaxation iteration solution G D (x,z,w) and the normal derivative B (x,z,w) of the first-order ultra relaxation iteration solution of reverse migration Green function are respectively solved by using source wavelet and prestack seismic record;The calculation result is used as the boundary condition of the next thin plate until the model bottom;G D (x,z,w) and B (x,z,w) are cross-correlated to obtain the final migration result.The application applies the method for solving scattering wave field by using ultra relaxation iteration solution to wave field continuation of prestack migration, overcomes the problem that the traditional Born scattering series is limited by weak scattering assumption and is only suitable for short-range propagation;The model is divided by referring to the thin plate approximation idea, the global relaxation factor is converted into local relaxation factor, the storage space is reduced, the calculation efficiency is improved without affecting the precision;The migration precision is improved.

Description

Technical Field

[0001] This invention belongs to the field of seismic scattered wave migration imaging technology, specifically involving the use of the ultra-relaxation iterative method to obtain the Green's function solution and its normal derivative of the Lippmann–Schwinger integral equation within a thin plate, and particularly relates to an ultra-relaxation iterative pre-stack migration imaging method based on thin plate approximation. Background Technology

[0002] As seismic work areas expand from simple to complex regions, and especially with the development of reflection seismic technology and the urgent need to address more complex geological problems, extremely complex seismic wave fields with interference from multiple wave groups are formed when the target geological body has a complex structure, well-developed faults, steep strata dips, abrupt lateral changes in lithology, or coexistence of heterogeneous geological bodies of different scales. Reflected and diffracted waves alone cannot provide detailed imaging of complex areas. However, scattered waves also carry geometric and physical information related to complex structures and lithologies. Therefore, it is crucial to achieve migration imaging of complex structures using both conventional reflected waves and scattered waves.

[0003] To achieve offset imaging of scattered waves, the problem of extending the scattered wave field is transformed into solving the Green's function. Using perturbation theory, the Helmholtz equations describing the scattering problem can be transformed into second-kind Fredholm-type integral equations, namely the Lippmann-Schwinger equations (LS equations). Therefore, the numerical simulation problem of scattered waves becomes a problem of solving the LS equations. Compared to other methods such as the transfer (propagation) matrix method, perturbation method, finite element method, and geometric ray method, the integral equation method is a powerful tool for scattered wave analysis due to its semi-analytical characteristics and the simplicity and ease of theoretical research and formula derivation.

[0004] Existing technology discloses a generalized over-relaxation iterative solution method for the Lippmann-Schwinger integral equation and its convergence characteristics. Applying this generalized over-relaxation iterative method to forward modeling of seismic scattering fields yields numerical simulation results comparable to finite-difference methods, even with large grid spacing and time sampling intervals. Existing technology also discloses a method utilizing… A parallel computational method for scattered seismic waves using discrete-time and FFT fast convolution is proposed, which rewrites the original LS equations and constructs equivalent LS equations. The method discretizes the matrix and, based on the Toeplitz property of the Green's function matrix, accelerates matrix-vector multiplication through Fast Fourier Transform, thereby reducing storage and computational complexity. However, the computational efficiency of this method remains low, and improving computational efficiency is a pressing issue. To date, the over-relaxation iterative solution method for solving the scattered wavefield has not been applied to wavefield extrapolation in pre-stack migration imaging. Summary of the Invention

[0005] The purpose of this invention is to provide a pre-stack migration imaging method based on thin plate approximation with super-relaxation iteration, in order to solve the problem of large storage space and low computational efficiency caused by the full rank of the discrete coefficient matrix when solving the Lippermann-Schwinger (LS) equation in the migration imaging process of the generalized super-relaxation iteration method.

[0006] The objective of this invention is achieved through the following technical solution:

[0007] First, the velocity model is divided into thin plates. Second, within the first thin plate, the first-order over-relaxation iterative solution G of the forward migration Green's function is obtained using the source wavelet and pre-stack seismic records. D The normal derivative B(x,z,w) of the first-order over-relaxation iterative solution of the inverse offset Green's function is calculated; then, the result is used as the boundary condition for the next thin plate until the bottom of the model; finally, G... D Cross-correlation imaging of (x,z,w) and B(x,z,w) is performed to obtain the final migration imaging result.

[0008] A super-relaxed iterative pre-stack migration imaging method based on thin-plate approximation includes the following steps:

[0009] A. Read the smoothed 2D velocity model. The number of horizontal grid points is x, the vertical length is z, and the velocity at each position in the model is v. xz The velocity model is divided into m thin plates along the z-direction.

[0010] B. Set the source parameters. Here, a Delta pulse is used with a dominant frequency of f Hz. The position of the first shot is (x a ,z a ), where a represents the number of shots, a = 1, 2, ..., n, channel spacing hm, a total of n shots, r detectors, channel spacing sm, sampling interval Δt, sampling length N, grid spacing dx and dz respectively, and the thickness of each thin plate is dz;

[0011] C. Within the first thin plate, a Delta pulse is used to excite the source location, serving as the initial value G(x,z) for the Green's function. j (x, z)), where t is the propagation time, (x, z)j＝1 ) represents spatial location, z j Let z be the depth of the j-th thin plate (j = 1, 2, ..., m). j The initial value is the shot point position (x0, z0);

[0012] D. Using Fourier transform to transform the Green's function G(x,z) j ,t) transforms from the time-space domain to the frequency-wavenumber domain The specific formula is as follows:

[0013]

[0014] Where i is the imaginary unit, k x is the transverse wavenumber, and w is the angular frequency, which is numerically equal to 2πf;

[0015] E. Solving for the Green's function value within the thin plate, i.e., the forward-scattered Green's function field, based on the thin plate approximation theory and using the first-order over-relaxation iterative method:

[0016]

[0017] in, This is the inverse Fourier transform. k0(z) j ) = w / c j c j For z j Minimum speed within a thin plate. For the velocity disturbance of the medium, Let L be the relaxation factor, where L = IA, I is the identity matrix, and A = G0(x, z) j +dz;x,z j ;w)O(x,z)dz, where Equation (2) is to perform calculations alternately in the frequency wavenumber domain and the frequency space domain, with phase shift processing in the frequency wavenumber domain and time shift correction in the frequency space domain;

[0018] F. Calculate the result G(x,z) from step E. j +dz,w) is used as the boundary condition for the next thin plate. Step E is repeated to obtain the forward scattering Green's function values ​​in different thin plates.

[0019] G. When the calculation reaches the bottom of the model, the loop stops, and the Green's function values ​​at different frequencies for each depth are saved (i.e., the forward propagation extension field G). D (x,z,w));

[0020] H. Perform a Fourier transform on the seismic record d(x,z,t) obtained at the epicenter after removing the direct wave, transforming it from the time-space domain to the frequency-wavenumber domain. This serves as the initial value for the first thin plate;

[0021] I. Solving for the reverse Green's function value within the thin plate using the first-order over-relaxation iterative method:

[0022]

[0023] Here, because the solution is for the upward wave, the sign of i is positive. Since the reverse recursion is performed, time is recursively reversed, so the sign is negative. Therefore, the reverse recursion formula and the forward recursion formula are the same in form, but have different physical meanings. Equation (3) is also calculated alternately in the frequency wavenumber domain and the frequency space domain.

[0024] J. Calculate the normal derivative of the Green's function within the thin plate:

[0025]

[0026] in, Representing the normal derivative, sinθ=k x / k0;

[0027] K. Calculate the result D(x,z) from step I. j +dz,w), as the boundary condition for the next thin plate, repeat steps nine and ten to obtain the backscattering Green's function values ​​D(x,z) within different thin plates. j +dz,w) and its normal derivative value B(x,z) j ,w);

[0028] L. When the calculation reaches the bottom of the model, the loop stops and the derivative values ​​B(x,z,w) of the Green's function at different frequencies for each depth are saved.

[0029] M. Perform cross-correlation imaging to obtain the final migration imaging result I1(x,z), calculated as follows:

[0030] I1(x,z)=∫dωG D (x,z,w)B(x,z,w) (5)

[0031] N. Repeat steps B through M for the remaining n-1 shots to obtain the offset imaging results for each single shot: I2(x,z), I3(x,z), ..., I n (x,z);

[0032] O. Superimpose all single-shot migration imaging results to obtain the final migration imaging result I(x,z):

[0033] I(x,z)=I1(x,z)+I2(x,z)+…+In (x,z) (6).

[0034] Further, step E specifically includes:

[0035] E1. Divide the velocity model into m thin plates along the z-axis;

[0036] E2. When within the first thin plate, the detla function is given as the initial value. Perform a Fourier transform on it to convert it from the time domain to the frequency domain;

[0037] E3. Rewrite the Lippmann–Schwinger integral equation as follows:

[0038] G(x,z j ,w)=G0(x a ,z a ,w)+AG(x,z j ,w) (7)

[0039] in, in The Green function for the dielectric background within the first thin plate is expressed as follows: G0(x a ,z a ,w) is the initial value of the detla function, and the above expression is written in the following form:

[0040] (IA)G(x,z j ,w)=G0(x a ,z a ,w) (8)

[0041] E4, the first-order iterative solution of equation (8) is in the form of:

[0042]

[0043] Where β is the hyperrelaxation factor; IL = A;

[0044] E5. To obtain an exact solution, we need to find the remaining vector r1 = G(x,z). a ,w)-βLG(x,z a Find the minimum value of ,w), and take β such that ||r1||=||G(x,z) a ,w)-βLG(x,z a If ,w)|| is minimal in the Hilbert space, then we have:

[0045]

[0046] Where G0=G(x,z) a ,w), G0 is the complex conjugate of G0. A necessary condition for satisfying the above equation is:

[0047]

[0048] get:

[0049]

[0050] E6. Substituting L=IA into equation (9), we get:

[0051]

[0052] E7. Using the thin-plate approximation theory, pass the result to the next thin plate using the following formula:

[0053]

[0054] in, This is the inverse Fourier transform. k0(z) j ) = w / c j c j For z j Minimum speed within a thin plate.

[0055]

[0056] For the velocity disturbance of the medium, Let L be the relaxation factor, where L = IA, I is the identity matrix, and A = G0(x, z) j +dz;x,z j ;w)O(x,z)dz, where Equation (2) involves alternating calculations in the frequency wavenumber domain and the frequency space domain, with phase shifting in the frequency wavenumber domain and time shift correction in the frequency space domain.

[0057] Compared with the prior art, the beneficial effects of the present invention are:

[0058] 1. This invention applies the over-relaxation iterative solution method for solving the scattered wave field to the wave field extension of pre-stack migration imaging, overcoming the problem that the traditional Born scattering series is limited by the weak scattering assumption and is only applicable to short-range propagation.

[0059] 2. Drawing on the idea of ​​thin plate approximation, the model is divided into parts, and the global relaxation factor is transformed into a local relaxation factor, reducing storage space and thus greatly reducing the amount of computation, improving computational efficiency without affecting accuracy.

[0060] 3. The normal derivative of the inverse propagation field based on seismic records is obtained (similar to high-pass filtering of the signal), which improves the accuracy of migration imaging. Attached Figure Description

[0061] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0062] Figure 1 This is a flowchart of the algorithm for the super-relaxed iterative pre-stack migration imaging method based on thin-plate approximation of the present invention;

[0063] Figure 2 It is a numerical model;

[0064] Figure 3 It is a single-frequency forward propagation seismic wave field;

[0065] Figure 4a This is the seismic record of the 11th blast;

[0066] Figure 4b Single-shot seismic records after removing direct waves;

[0067] Figure 5 This is the result of offset imaging.

[0068] In the diagram: 1. Earthquake source 2. Detector 3. First wave impedance interface 4. Second wave impedance interface 5. Third wave impedance interface 6. Fourth wave impedance interface. Detailed Implementation

[0069] The present invention will be further described below with reference to embodiments:

[0070] The present invention will now be described in further detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and not intended to limit it. Furthermore, it should be noted that, for ease of description, the accompanying drawings show only the parts relevant to the present invention, and not all of the structures.

[0071] It should be noted that similar reference numerals and letters in the following figures indicate similar items; therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures. Furthermore, in the description of this invention, terms such as "first," "second," etc., are used only to distinguish descriptions and should not be construed as indicating or implying relative importance.

[0072] This invention draws on the idea of ​​thin-plate approximation. By dividing the velocity model into thin plates, the global relaxation factor is transformed into a local relaxation factor, significantly reducing the computational load. Furthermore, compared to the global relaxation factor, the velocity perturbation within the thin plate is relatively smaller, resulting in higher accuracy. Therefore, firstly, the relaxation factor within the first thin plate is obtained. Then, the forward scattering field within the first thin plate is solved using the over-relaxation iterative method. This result serves as the initial condition for the next thin plate, and the calculation is repeated sequentially until the bottom of the model is reached, obtaining the first-order over-relaxation iterative solution of the forward migration Green's function. Next, pre-stack seismic records are used as the initial condition for the first thin plate for directional recursion. Unlike the forward propagation wavefield, the normal derivative B(x,z,w) of the first-order over-relaxation iterative solution of the reverse migration Green's function is used until the bottom of the model is reached. Finally, the derivatives of the forward and reverse propagation continuation fields are cross-correlated for imaging, ultimately obtaining the pre-stack migration imaging result.

[0073] This invention relates to an ultra-relaxed iterative pre-stack migration imaging method based on thin-plate approximation, comprising the following steps:

[0074] Step 1: Read the smoothed 2D velocity model. The number of horizontal grid points is x, the vertical length is z, and the velocity at each position in the model is v. xz The velocity model is divided into m thin plates along the z-direction.

[0075] Step 2: Set the source parameters. Here, a Delta pulse is used with a dominant frequency of f Hz. The position of the first shot is (x... a ,z a ), where a represents the number of shots, a = 1, 2, ..., n. The channel spacing is hm, there are a total of n shots, r detectors, channel spacing sm, sampling interval Δt, sampling length N, grid spacings dx and dz, and the thickness of each plate is dz.

[0076] Step 3: Within the first thin plate, use a Delta pulse to excite at the source location, which serves as the initial value G(x,z) for the Green's function. j (x, z)), where t is the propagation time, (x, z) j＝1 ) represents spatial location, z j Let z be the depth of the j-th thin plate (j = 1, 2, ..., m). j The initial value is the shot point position (x). a ,z a );

[0077] Step 4: Use Fourier transform to transform the Green's function G(x,z) j ,t) transforms from the time-space domain to the frequency-wavenumber domain The specific formula is as follows:

[0078]

[0079] Where i is the imaginary unit, k x denoted as transverse wavenumber, and w as angular frequency, which is numerically equal to 2πf.

[0080] Step 5: Solve for the Green's function value within the thin plate, i.e., the forward-scattered Green's function field, based on the thin plate approximation theory and using the first-order over-relaxation iterative method:

[0081]

[0082] in, This is the inverse Fourier transform. k0(z) j ) = w / c j c j For z j Minimum speed within a thin plate. For the velocity disturbance of the medium, Let L be the relaxation factor, where L = IA, I is the identity matrix, and A = G0(x, z) j +dz;x,z j ;w)O(x,z)dz, where Equation (2) involves alternating calculations in the frequency wavenumber domain and the frequency space domain, with phase shifting in the frequency wavenumber domain and time shift correction in the frequency space domain.

[0083] Step Six: Calculate the result G(x,z) from Step Five. j +dz,w) is used as the boundary condition for the next thin plate. Step 5 is repeated to obtain the forward scattering Green function values ​​in different thin plates.

[0084] Step 7: When the calculation reaches the bottom of the model, stop the loop and save the Green's function values ​​(i.e., the forward extension field G) at different frequencies for each depth. D (x,z,w)).

[0085] Step 8: Perform a Fourier transform on the seismic record d(x,z,t) obtained at the epicenter after removing the direct wave, transforming it from the time-space domain to the frequency-wavenumber domain. This serves as the initial value for the first thin plate.

[0086] Step 9: Solve for the inverse Green's function value within the thin plate using the first-order over-relaxation iterative method:

[0087]

[0088] Here, the solution makes the upward wave i positive, and since the reverse recursion is performed, time is recursively reversed, so the reverse recursion formula is negative. Therefore, the reverse recursion formula and the forward recursion formula are the same in form, but have different physical meanings. Equation (3) is also calculated alternately in the frequency wavenumber domain and the frequency space domain.

[0089] Step 10: Calculate the normal derivative of the Green's function within the thin plate:

[0090]

[0091] in, Representing the normal derivative, sinθ=k x / k0.

[0092] Step 11: Calculate the result D(x,z) from Step 9. j +dz,w), as the boundary condition for the next thin plate, repeat steps nine and ten to obtain the backscattering Green's function values ​​D(x,z) within different thin plates. j +dz,w) and its normal derivative value B(x,z) j ,w).

[0093] Step 12: When the calculation reaches the bottom of the model, stop the loop and save the derivative values ​​B(x,z,w) of the Green's function at different frequencies for each depth.

[0094] Step 13: Perform cross-correlation imaging to obtain the final migration imaging result I1(x,z), calculated as follows:

[0095] I1(x,z)=∫dωG D (x,z,w)B(x,z,w) (5)

[0096] Step Fourteen: Repeat Steps Two through Thirteen for the remaining n-1 shots to obtain the offset imaging results I2(x,z), I3(x,z), ..., I... n (x,z).

[0097] Step 15: Stack all single-shot migration imaging results to obtain the final migration imaging result I(x,z):

[0098] I(x,z)=I1(x,z)+I2(x,z)+…+I n (x,z) (6).

[0099] To better illustrate the implementation effect of the above specific algorithm process, a specific example is given below.

[0100] Specifically, step five includes the following steps:

[0101] 1. Divide the velocity model into m thin plates along the z-axis;

[0102] 2. When within the first thin plate, use the `detla` function as the initial value. Perform a Fourier transform on it to convert it from the time domain to the frequency domain;

[0103] 3. Rewrite the Lippmann–Schwinger integral equation as follows:

[0104] G(x,z j ,w)=G0(x a ,z a ,w)+AG(x,z j ,w) (7)

[0105] in in The Green function for the dielectric background within the first thin plate is expressed as follows: G0(x a ,z a ,w) is the initial value of the detla function, and the above expression is written in the following form:

[0106] (IA)G(x,z j ,w)=G0(x a ,z a ,w) (8)

[0107] 4. The first-order iterative solution of equation (8) is in the form of:

[0108]

[0109] Where β is the hyperrelaxation factor; IL = A;

[0110] 5. To obtain an accurate solution, we need to find the remaining vector r1 = G(x,z). a ,w)-βLG(x,z a Find the minimum value of ,w), and take β such that ||r1||=||G(x,z) a ,w)-βLG(x,z a If ,w)|| is minimal in the Hilbert space, then we have:

[0111]

[0112] Where G0=G(x,z) a ,w), G0 is the complex conjugate of G0. A necessary condition for satisfying the above equation is:

[0113]

[0114] get:

[0115]

[0116] 6. Substituting L=IA into equation (9), we get:

[0117]

[0118] 7. Using the thin-plate approximation theory, the result can be passed to the next thin plate using the following formula:

[0119]

[0120] in This is the inverse Fourier transform. k0(z) j ) = w / c j c j For z j Minimum speed within a thin plate.

[0121] For the velocity disturbance of the medium, Let L be the relaxation factor, where L = IA, I is the identity matrix, and A = G0(x, z) j +dz;x,z j ;w)O(x,z)dz, where Equation (2) involves alternating calculations in the frequency wavenumber domain and the frequency space domain, with phase shifting in the frequency wavenumber domain and time shift correction in the frequency space domain.

[0122] Example 1

[0123] like Figure 1 As shown, the present invention provides an ultra-relaxed iterative pre-stack migration imaging method based on thin-plate approximation, comprising the following steps:

[0124] 1. Read the smoothed two-dimensional velocity model (e.g.) Figure 2 As shown, the horizontal grid has 401 points and the vertical length is 301. The velocity model is divided into thin plates along the z-direction, totaling 300 plates.

[0125] 2. Set the source parameters. Here, Delta pulse is used with a main frequency of 20Hz. The first source position is (0,0). The channel spacing is 105m, with a total of 21 shots. There are 401 detectors with a channel spacing of 5m, a sampling interval of 0.001, a sampling length of 1750, and grid spacings of 5m and 5m respectively. The thickness of each plate is 5m.

[0126] 3. Within the first thin plate, a Delta pulse is used to excite the source location, serving as the initial value G(x,z) for the Green's function. j (x, z)), where t is the propagation time, (x, z) j＝1 ) represents spatial location, z j Let z be the depth of the j-th thin plate (j = 1, 2, ..., m).j The initial value is the shot point position (x0, z0).

[0127] Using Fourier transform to transform the Green's function G(x,z) j ,t) transforms from the time-space domain to the frequency-wavenumber domain The specific formula is as follows:

[0128]

[0129] Where i is the imaginary unit, k x is the transverse wavenumber, and w is the angular frequency, which is numerically equal to 2πf;

[0130] Obtain the initial frequency domain values ​​of the first thin plate.

[0131] 4. Using the following steps, alternate calculations are performed in the frequency wavenumber domain and frequency space domain within the first thin plate to obtain the forward Green function value based on the first-order over-relaxation iteration of the thin plate approximation theory.

[0132] Specifically, the Green function values ​​within the thin plate, i.e. the forward-scattered Green function field, are solved based on the thin plate approximation theory and using the first-order over-relaxation iterative method:

[0133]

[0134] in, This is the inverse Fourier transform. k0(z) j ) = w / c j c j For z j Minimum speed within a thin plate. For the velocity disturbance of the medium, Let L be the relaxation factor, where L = IA, I is the identity matrix, and A = G0(x, z) j +dz;x,z j ;w)O(x,z)dz, where Equation (2) involves alternating calculations in the frequency wavenumber domain and the frequency space domain, with phase shifting in the frequency wavenumber domain and time shift correction in the frequency space domain.

[0135] 5. Using the result of step 4 above as the initial value for the next thin plate, the forward Green function value of the second thin plate is obtained by alternating calculations in the frequency wavenumber domain and frequency space domain using the above method.

[0136] 6. Repeat steps 4 and 5 to obtain the Green's function values ​​at different frequencies for each depth, and the forward propagation extension field G. D (x,z,w) (e.g.) Figure 3 (as shown);

[0137] 7. Perform a Fourier transform on the seismic record after removing direct waves (as shown in Figure 4) to obtain the initial boundary values ​​of the first thin plate.

[0138] 8. Using the method described below, alternately calculate in the frequency wavenumber domain and frequency space domain to obtain the backscattering Green's function value D(x,z) within the first thin plate. j +dz,w) and its normal derivative value B(x,z) j ,w).

[0139] The reverse Green's function value within the thin plate is obtained using a first-order over-relaxation iterative method:

[0140]

[0141] Find the normal derivative of the Green's function within the thin plate:

[0142]

[0143] in Representing the normal derivative, sinθ=k x / k0.

[0144] 9. Using the result of step 4 as the initial value for the next thin plate, repeat step 8, alternating between calculations in the frequency wavenumber domain and frequency space domain, to obtain the backscattering Green's function value D(x,z) within the second thin plate. j +dz,w) and its normal derivative value B(x,z) j ,w).

[0145] 10. Repeat steps 4 and 5 to obtain the derivative values ​​B(x,z,w) of the Green's function at different frequencies for each depth;

[0146] 11. Obtain single-shot offset imaging results using the following steps;

[0147] Cross-correlation imaging is performed to obtain the final migration imaging result I1(x,z), calculated as follows:

[0148] I1(x,z)=∫dωG D (x,z,w)B(x,z,w) (5)

[0149] 12. Repeat steps 2-11 to obtain the migration imaging results for all individual shots, and then superimpose the results according to the following steps to obtain the final migration imaging result (e.g., Figure 5 (As shown).

[0150] The final migration imaging result I(x,z) is obtained by superimposing all the single-shot migration imaging results:

[0151] I(x,z)=I1(x,z)+I2(x,z)+…+I n (x,z) (6).

[0152] Note that the above description is merely a preferred embodiment of the present invention and the technical principles employed. Those skilled in the art will understand that the present invention is not limited to the specific embodiments described herein, and various obvious changes, readjustments, and substitutions can be made without departing from the scope of protection of the present invention. Therefore, although the present invention has been described in detail through the above embodiments, the present invention is not limited to the above embodiments, and may include many other equivalent embodiments without departing from the concept of the present invention, the scope of which is determined by the scope of the appended claims.

Claims

1. A super-relaxed iterative pre-stack migration imaging method based on thin-plate approximation, characterized in that, Includes the following steps: A. Read the smoothed 2D velocity model; the number of horizontal grid points is... x The longitudinal length is z The velocity at each position of the model is The velocity model is divided into m thin plates along the z-direction. B. Set the source parameters, using a Delta pulse with a dominant frequency of [value missing]. f Hz, the position of the first shot is ( ), where 'a' represents the number of guns, lane spacing h m, n guns in total, r detectors, channel spacing s m, sampling interval is The sampling length is N, and the grid spacing is respectively , Each sheet is 100mm thick. ; C. Within the first thin plate, a Delta pulse is used to excite the source location, serving as the initial value for the Green's function. t is the propagation time, ( ) represents the location in space. z j Let m be the depth of the j-th thin plate, where j = 1, 2, ..., m. z j To determine the initial shot location ; D. Using Fourier transform to transform the Green's function From the time-space domain to the frequency-wavenumber domain The specific formula is as follows: （1） in, i The imaginary unit, For transverse wavenumber, w It is the angular frequency, numerically equal to ; E. Solving for the Green's function value within the thin plate, i.e., the forward-scattered Green's function field, based on the thin plate approximation theory and using the first-order over-relaxation iterative method: （2） in, This is the inverse Fourier transform; c j for Minimum speed within a thin plate. , For the velocity disturbance of the medium, Let be the relaxation factor, where I is the identity matrix. ,in Equation (2) is to perform calculations alternately in the frequency wavenumber domain and the frequency space domain, with phase shift processing in the frequency wavenumber domain and time shift correction in the frequency space domain; F. Calculate the results of step E. As the boundary condition for the next thin plate, repeat step E to obtain the forward scattering Green's function values ​​in different thin plates; G. When the calculation reaches the bottom of the model, the loop stops, and the Green's function values ​​at different frequencies for each depth are saved, which are the forward propagation extension fields. ; H. Seismic records obtained at the epicenter, after removing direct waves. Perform a Fourier transform to transform from the time-space domain to the frequency-wavenumber domain. , as the initial value for the first thin plate; I. Solving for the reverse Green's function value within the thin plate using the first-order over-relaxation iterative method: （3） J. Calculate the normal derivative of the Green's function within the thin plate: （4） in, Representing the normal derivative, , , ; K. The calculation results of step I Using this as the boundary condition for the next thin plate, repeat step IJ to obtain the backscattered Green's function values ​​within different thin plates. and its normal guide value ; L. When the calculation reaches the bottom of the model, the loop stops, and the derivative values ​​of the Green's function at different frequencies for each depth are saved. ; M. Perform cross-correlation imaging to obtain the final migration imaging result. The calculation formula is as follows: （5） N. Repeat steps B through M for the remaining n-1 shots to obtain the offset imaging result for each single shot. ; O. Superimpose all single-shot offset imaging results to obtain the final offset imaging result: （6）。 2. The ultra-relaxed iterative pre-stack migration imaging method based on thin-plate approximation according to claim 1, characterized in that, Step E specifically includes: E1. Divide the velocity model into m thin plates along the z-axis; E2. When within the first thin plate, the detla function is given as the initial value. Then perform a Fourier transform on it to convert it from the time domain to the frequency domain; E3. Rewrite the Lippmann–Schwinger integral equation as follows: （7） in, ;in The Green function for the dielectric background within the first thin plate is expressed as follows: , Let be the initial value of the detla function, and write the above expression in the following form: （8） E4, the first-order iterative solution of equation (8) is in the form of: （9） in, It is a relaxation factor; E5. To obtain an exact solution, we need to find the remaining vector. The minimum value is taken. make If it is minimal in Hilbert space, then: （10） in, , for Complex conjugate; a necessary condition for satisfying the above equation is: （11） get: （12） E6, will Substituting into equation (9), we get: （13） E7. Using the thin-plate approximation theory, pass the result to the next thin plate using the following formula: （2） in, For inverse Fourier transform, c j The minimum speed within a thin plate. , , For the velocity disturbance of the medium, Let be the relaxation factor, where I is the identity matrix. ,in Equation (2) is to perform calculations in the frequency wavenumber domain and the frequency space domain alternately. The frequency wavenumber domain is phase-shifted and the frequency space domain is time-shifted.

Images