A Subsurface Imaging Method Based on the Renormalized Marchenko Method

Abstract

This application belongs to the field of exploration earth technology and is a subsurface imaging method based on the renormalized Marchenko method. The method includes: removing direct waves from seismic data and performing seismic deconvolution to obtain preprocessed seismic data; selecting subsurface imaging points and setting them as subsurface focal points; estimating the direct wave Green's function from the subsurface focal point to the surface receiver point based on a background velocity model; reversing the direct wave Green's function on the time axis to obtain the initial downlink focusing function, and setting the initial uplink focusing function to zero; iteratively calculating the initial downlink and uplink focusing functions to obtain the iterated downlink and uplink focusing functions; calculating the downlink and uplink Green's functions based on the uplink and downlink focusing functions; and imaging the subsurface imaging points. This method ensures absolute convergence of the iteration by adjusting the convergence parameters, ultimately obtaining subsurface imaging results unaffected by multiple waves.

Description

Technical Field

[0001] This application belongs to the field of exploration earth technology, and in particular relates to a subsurface imaging method based on the renormalized Marchenko method. Background Technology

[0002] In the field of exploration seismology, the Marchenko method for calculating subsurface Green's function was proposed in 2014. The Marchenko method is a method that uses surface-generated artificial seismic excitation and reception data to calculate the Green's function of subsurface virtual point sources excited and received at the surface. It can effectively handle subsurface interlayer multiples, and the calculated subsurface Green's function contains accurate multiple information compared to traditional calculation methods proposed in the 1980s and 1990s. The Marchenko method derives the Marchenko equation from the seismological reciprocity theorem, calculates the focusing function using the Newman series iterative solution, and then calculates the subsurface Green's function based on the focusing function using the Marchenko method. The Newman series iterative solution requires specific convergence conditions; for seismic data that do not meet these conditions, the Green's function cannot be calculated through normal iterative convergence. Seismic imaging technology based on the Marchenko method can effectively suppress multiple artifacts in subsurface imaging results. In particular, when dealing with seismic data in subsurface strongly scattering media containing surface-correlated multiples, the conventional Marchenko method almost fails to converge and cannot suppress multiple artifacts. Summary of the Invention

[0003] This application provides a subsurface imaging method based on the renormalized Marchenko method, which solves the problem of non-convergence of conventional Marchenko methods with seismic data.

[0004] A seismic imaging method based on the renormalized Marchenko method, according to an embodiment of this application, includes:

[0005] Direct waves are removed from the seismic data, and seismic deconvolution is performed to obtain preprocessed seismic data.

[0006] Select underground imaging points and set them as underground focal points. Estimate the direct Poglin function from the underground focal point to the surface receiving point based on the background velocity model.

[0007] The direct-to-Borglin function is inverted on the time axis to serve as the initial downlink focusing function, and the initial uplink focusing function is set to zero. Iterative calculations are then performed based on the initial downlink and uplink focusing functions to obtain the iteratively derived downlink and uplink focusing functions. These iterative calculations include:

[0008] Calculate the downlink focusing function of the previous iteration and The first difference of the uplink focusing function is calculated; the temporal convolution of the reflected data and the first difference is calculated to obtain the convolution result; the convolution result is spatially integrated, multiplied by the temporal window function and the convergence parameter, and then compared with the uplink focusing function of the previous iteration. By doubling the sum, we obtain the upward focusing function for the next iteration. For convergence parameters, The surface reflectance coefficient;

[0009] Calculate the uplink focusing function of the previous iteration and The second difference of the downlink focusing function is calculated, and the temporal convolution of the reflected data inverted on the time axis with the second difference is obtained. The convolution result is then spatially integrated, multiplied by the temporal window function and the convergence parameter, and finally combined with the downlink focusing function from the previous iteration. Double the sum and add it to the initial downlink focusing function. By multiplying by a factor of two, we obtain the downlink focusing function for the next iteration;

[0010] Calculate the downlink Green's function and the uplink Green's function based on the uplink focus function and the downlink focus function;

[0011] Underground imaging points are imaged using the upward focusing function, downward focusing function, downward Green's function, and upward Green's function.

[0012] Furthermore, the time window function obtained by determining the travel time of seismic waves is expressed as: ,in, It is half the time length of the seismic wavelet. For the travel time of seismic waves, The absolute value of time. Indicates the surface receiving point. Indicates an underground focal point. This is a time window function.

[0013] Furthermore, the downlink Green's function is calculated based on the uplink and downlink focusing functions, specifically including: calculating the downlink focusing function and... The third difference of the uplink focusing function is obtained by temporally convolving the third difference with the reflection data; the convolution result is obtained by spatial integration of the convolution result and added to the negative of the uplink focusing function to obtain the downlink Green's function.

[0014] Furthermore, the upward Green's function is calculated based on the upward and downward focusing functions. Specifically, this includes reversing the reflected data on the time axis and calculating the upward focusing function and... The fourth difference of the downlink focusing function is obtained by temporal convolution of the fourth difference with the inverted reflection data on the time axis. The convolution result is obtained by spatial integration of the convolution result and added to the negative of the downlink focusing function to obtain the uplink Green's function.

[0015] Furthermore, a first imaging result is obtained by using an upward focusing function, a downward focusing function, and a downward Green's function for dual-focusing imaging, specifically including:

[0016] Calculate the downlink focusing function and The fifth difference of the upward focusing function;

[0017] The fifth difference is convolved with the descending Green's function over time to obtain the convolution result.

[0018] The underground imaging results are obtained by performing spatial integration on the convolution results.

[0019] Furthermore, the second imaging result is obtained by cross-correlating the upward Green's function and the downward Green's function.

[0020] Furthermore, the background velocity model is the longitudinal wave velocity of the underground medium.

[0021] Compared with the prior art, the beneficial effects of this application are as follows: This application breaks through the convergence condition of the conventional Marchenko algorithm. The conventional Marchenko algorithm cannot iteratively converge seismic data containing surface-correlated multiples. The method of this application can ensure absolute convergence of the iteration by adjusting the convergence parameters, and finally obtain subsurface imaging results that are not affected by multiples. Attached Figure Description

[0022] Figure 1 A flowchart illustrating a subsurface imaging method based on the renormalized Marchenko method, provided in this application embodiment;

[0023] Figure 2 A schematic diagram of the velocity model and density model provided in the embodiments of this application;

[0024] Figure 3 A schematic diagram of reflection data provided in an embodiment of this application;

[0025] Figure 4 A schematic diagram of the background velocity model provided in the embodiments of this application;

[0026] Figure 5 The diagram shows the (a) initial focusing function, (b) downlink focusing function, and (c) uplink focusing function corresponding to the focusing imaging points provided in the embodiments of this application.

[0027] Figure 6A plot of (a) Green's function, (b) upward Green's function, and (c) downward Green's function calculated for the embodiments of this application;

[0028] Figure 7 The image shows the underground imaging results provided in the embodiments of this application. Detailed Implementation

[0029] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0030] See Figure 1 The seismic imaging method shown includes:

[0031] S1 removes direct waves from the seismic data and performs seismic deconvolution to obtain preprocessed seismic data.

[0032] Seismic data includes free-surface multiples, which are acquired by surface receivers after blasting at the surface. Direct waves are waves that propagate directly from the epicenter to the receiving point without reflection or refraction at any interface. They propagate along the shortest path, typically a straight path. Removing direct waves from seismic data usually refers to removing high-energy direct waves near the epicenter. Preprocessing is a noise reduction process. After removing high-energy direct waves near the epicenter, the remaining seismic data is subjected to seismic deconvolution. Seismic deconvolution eliminates the influence of seismic wavelets from seismic data, restoring the original form of the subsurface surface reflection coefficients, thereby improving the resolution and signal-to-noise ratio of the seismic data. This is a conventional processing method. The preprocessed seismic data is recorded as follows: , This represents reflection data. All reflection data below refers to preprocessed seismic data. (Subscript) The Earth's surface is represented by free boundary conditions. Represents the surface receiving point. Represents the coordinates of the excitation point. Represents time. Includes reflected waves and multiples (interlayer multiples and surface-related multiples).

[0033] S2 selects underground imaging points and sets them as underground focal points. Based on the background velocity model, the direct Poglin function from the underground focal point to the surface receiving point is estimated.

[0034] For the underground area that needs to be imaged, a background velocity model needs to be prepared and recorded as follows: , Represents depth, The background velocity model represents the P-wave velocity of the subsurface medium, and its accuracy requirements are consistent with those of conventional imaging methods, i.e., a smooth background velocity model is sufficient. In conventional seismic data processing, it can be obtained through travel-time tomography or full waveform inversion techniques.

[0035] Select an underground imaging point and set it as the underground focal point. The coordinates are represented as follows: Based on a background velocity model Estimate underground focal point To the surface receiving point direct Bogle function subscript Represents the direct wave component. Direct wave Green's function. The Green's function refers to the direct wave that propagates directly from the underground focal point to the surface receiving point.

[0036] Estimate the direct Poglin function from the underground focal point to the surface receiver point based on the background velocity model. This step can be directly obtained by solving the wave equation using the finite difference method, or by solving the equation to calculate the travel time. The travel time curve is obtained by convolving it with the wavelet.

[0037] S3 reverses the direct Bogle function on the time axis as the initial downlink focusing function, and sets the initial uplink focusing function to zero;

[0038] The direct-to-Borglin function is inverted on the time axis and used as the initial downlink focusing function, i.e. .here Represents the focusing function, superscript Represents the downlink wave field, in the subscript Represents the number of iterations, indicated by the subscript. This represents an iteration count of zero, i.e., the initial downlink focusing function. This represents the direct portion of the focusing function; here, it refers to the direct portion of the downlink focusing function. With the initial downlink focusing function They are the same. The initial uplink focus function is directly set to zero. superscript Indicates an upward wave field;

[0039] S4 performs iterative calculations based on the initial downlink focusing function and the initial uplink focusing function to obtain the iteratively calculated downlink focusing function and uplink focusing function. The iterative calculations include:

[0040] Calculate the downlink focusing function of the previous iteration and The first difference of the uplink focusing function is calculated; the reflection data is temporally convolved with the first difference to obtain the convolution result; the convolution result is spatially integrated, multiplied by the temporal window function and the convergence parameter, and then compared with the uplink focusing function of the previous iteration. By doubling the sum, we obtain the upward focusing function for the next iteration. These are the convergence parameters;

[0041] Calculate the uplink focusing function of the previous iteration and The second difference of the downlink focusing function is calculated, and the reflection data inverted on the time axis is convolved with the second difference to obtain the convolution result. The convolution result is then spatially integrated, multiplied by the time window function and the convergence parameter, and finally combined with the downlink focusing function from the previous iteration. Double the sum and add it to the initial downlink focusing function. By multiplying by a factor of two, we obtain the downlink focusing function for the next iteration;

[0042] S5 calculates the downlink Green's function and the uplink Green's function based on the uplink focus function and the downlink focus function;

[0043] S6 images the underground imaging points based on the upward focusing function, downward focusing function, downward Green's function, and upward Green's function.

[0044] In one embodiment, the calculation formulas for the iterative downlink focusing function and the uplink focusing function are obtained by iterative calculation based on the initial downlink focusing function and the initial uplink focusing function.

[0045] ,

[0046] in, For the number of iterations The upward focusing function, For the number of iterations -1 is the upward focusing function. For convergence parameters, For time window functions, Representing the Earth's surface, the integration process proceeds along the Earth's surface. Integrating the excitation sampling, The time variable is Reflection data at that time Indicates the number of iterations. The time is The downlink focusing function, This is the surface reflectance, which reflects the reflection effect of a free surface; it is generally set to -1. Indicates the number of iterations. The upward focusing function, For the number of iterations The downlink focusing function, For the number of iterations The downlink focusing function, The time variable is Reflection data at that time.

[0047] This application embodiment uses convergence parameters. This changed the calculation steps for each iteration.

[0048] In one embodiment, the time window function is obtained by determining the travel time of the seismic wave, and is expressed as: ,in, It is half the time length of the seismic wavelet. For the travel time of seismic waves, The absolute value of time. Indicates the surface receiving point. Indicates an underground focal point. For time window function, and They refer to the same meaning.

[0049] In the entire iterative process described above, in order to improve the computation speed, the convolution calculation can also use the Fast Fourier Transform to convert the two to the frequency domain for multiplication and then perform an inverse Fourier transform to convert them back to the time domain. This calculation is equivalent to directly calculating the time convolution.

[0050] By adjusting the convergence parameters This ensures absolute convergence of the iterative process, and the convergence parameters... The value of is between 0 and 1. The smaller the value, the slower the convergence speed. Reducing the convergence parameter... It is possible to improve convergence performance while reducing convergence speed. The settings can be adjusted according to the data to ensure convergence of results at a higher convergence speed.

[0051] In one embodiment, calculating the downlink Green's function based on the uplink focusing function and the downlink focusing function specifically includes: calculating the downlink focusing function and... The third difference of the uplink focusing function is obtained by temporally convolving the third difference with the reflection data; the convolution result is obtained by spatial integration of the convolution result and added to the negative of the uplink focusing function to obtain the downlink Green's function.

[0052] The upward Green's function is calculated based on the upward and downward focusing functions, specifically including: reversing the reflected data on the time axis, calculating the upward focusing function and... The fourth difference of the downlink focusing function is obtained by temporal convolution of the fourth difference with the inverted reflection data on the time axis. The convolution result is obtained by spatial integration of the convolution result and added to the negative of the downlink focusing function to obtain the uplink Green's function.

[0053] After several iterations and stabilization, the upward Green's function of the surface reception excited by the underground focal point is obtained using the obtained downlink and uplink focusing functions. With the downlink Green's function Green's function subscript This indicates that it includes free surface multiple waves:

[0054] ,

[0055] ,

[0056] The above calculation process does not require iteration and does not have a time window function. It is the expression for the upward focusing function obtained after iteration when performing integration. This is the expression for the downward focusing function obtained after iteration when performing integration. It is the expression for the downlink focusing function obtained after iteration. It is the expression for the downlink focusing function obtained after iteration. Represents a time variable.

[0057] In one embodiment, a first imaging result is obtained by performing dual-focusing imaging using an up-focusing function, a down-focusing function, and a down-Green's function, including:

[0058] Calculate the downlink focusing function and The fifth difference of the upward focusing function;

[0059] The fifth difference is convolved with the descending Green's function over time to obtain the convolution result.

[0060] The underground imaging results are obtained by performing spatial integration on the convolution results.

[0061] The dual-focusing imaging formula used is:

[0062] ,

[0063] This application achieves its effect by adding a free surface correction term, namely... The underground focal point can be obtained. By analyzing the imaging values ​​at a given location and imaging all underground points, the underground imaging results can be obtained.

[0064] In another embodiment, a second imaging result is obtained by cross-correlating the upward Green's function and the downward Green's function. The first imaging result and the second imaging result can be cross-referenced and verified.

[0065] Taking model testing as an example, the following methods are used: Figure 2The velocity and density models shown are used as test models to illustrate the implementation process of this application, wherein... Figure 2 (a) in the figure represents the velocity model. Figure 2 (b) in the figure represents the density model, such as Figure 3 The diagram shown illustrates the reflection data, with the horizontal axis representing the surface receiving point. The interval between receiving points is The coordinate range is from arrive ,total The channel receives data, and the vertical axis represents time. The sampling interval is Time from arrive ,total One sampling point. Figure 3 The data displayed is the reflection data of a single shot, namely... The reflection data at that time, in addition to the data from arrive interval One shot total Artillery data;

[0066] Preparation Figure 4 The background velocity model shown;

[0067] Select an underground focal point, such as the horizontal distance. ,depth The point at that location, i.e. Based on the background velocity model, the travel time is obtained by solving the equation of motion. By performing a convolution operation between the travel time curve and the Ricker wavelet, and then reversing the time axis, the initial downlink focusing function is obtained. The time range of the initial downlink focusing function is set to... arrive The sampling interval is The coordinates of the surface receiving point are from arrive , interval is ,like Figure 5 As shown in (a), and using time travel Set the time window function;

[0068] Time window function Similarly, set the range from arrive The sampling interval is The coordinates of the surface receiving point from arrive , interval is ;

[0069] Iteratively solve for the convergent downlink and uplink focusing functions, and set the convergence parameters. This ensures that the results converge and the iteration continues. The subsequent result is as follows Figure 5 The graph of the downlink focusing function shown in Figure (b) and... Figure 5 (c) shows a plot of the upward focusing function. It can be seen that the focusing function result is stable.

[0070] Find the ascending and descending Green's functions as follows: Figure 6 (b) and Figure 6 As shown in (c), the Green's function can be obtained by summing the ascending Green's function and the descending Green's function, and the result is also in Figure 6 As shown in (a), it can be seen that the Green's function result is stable.

[0071] Calculate the imaging values ​​of underground imaging points;

[0072] Repeat this process until all required imaging points underground are calculated, based on the horizontal grid spacing. Vertical grid spacing Set the imaging grid to ,total By analyzing several imaging points, the final subsurface imaging results based on the renormalized Marchenko method are obtained, as shown below. Figure 7 As shown, the imaging process is completely stable with no non-convergence issues, and the multi-wave artifacts in the final imaging result are suppressed.

[0073] The above description is merely a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. A subsurface imaging method based on the renormalized Marchenko method, characterized in that, include: Direct waves are removed from the seismic data, and seismic deconvolution is performed to obtain preprocessed seismic data. Select underground imaging points and set them as underground focal points. Estimate the direct Poglin function from the underground focal point to the surface receiving point based on the background velocity model. The direct Bogle function is inverted on the time axis and used as the initial downlink focusing function, while the initial uplink focusing function is set to zero; Iterative calculations are performed based on the initial downlink focusing function and the initial uplink focusing function to obtain the iterative downlink focusing function and the uplink focusing function. The iterative calculations include: Calculate the downlink focusing function of the previous iteration and The first difference of the uplink focusing function is calculated; the temporal convolution of the reflected data and the first difference is calculated to obtain the convolution result; the convolution result is spatially integrated, multiplied by the temporal window function and the convergence parameter, and then compared with the uplink focusing function of the previous iteration. By doubling the sum, we obtain the upward focusing function for the next iteration. For convergence parameters, The surface reflectance coefficient; Calculate the uplink focusing function of the previous iteration and The second difference of the downlink focusing function is calculated, and the reflection data inverted on the time axis is convolved with the second difference to obtain the convolution result. The convolution result is then spatially integrated, multiplied by the time window function and the convergence parameter, and finally combined with the downlink focusing function from the previous iteration. Double the sum and add it to the initial downlink focusing function. By multiplying by a factor of two, we obtain the downlink focusing function for the next iteration; Calculate the downlink Green's function and the uplink Green's function based on the uplink focus function and the downlink focus function; Imaging of underground imaging points is performed based on the upward focusing function, downward focusing function, downward Green's function, and upward Green's function; The first imaging result is obtained by using dual-focusing imaging with an upward focusing function, a downward focusing function, and a downward Green's function, specifically including: Calculate the downlink focusing function and The fifth difference of the upward focusing function; The fifth difference is convolved with the descending Green's function over time to obtain the convolution result. The underground imaging results are obtained by performing spatial integration on the convolution results. The second imaging result is obtained by cross-correlating the upward Green's function and the downward Green's function.

2. The subsurface imaging method based on the renormalized Marchenko method according to claim 1, characterized in that, The time window function is obtained by determining the travel time of seismic waves, and is expressed as: ,in, It is half the time length of the seismic wavelet. For the travel time of seismic waves, The absolute value of time. Indicates the surface receiving point. Indicates an underground focal point. This is a time window function.

3. The subsurface imaging method based on the renormalized Marchenko method according to claim 1, characterized in that, Calculating the downlink Green's function based on the uplink and downlink focusing functions specifically includes: calculating the downlink focusing function and... The third difference of the uplink focusing function is obtained by temporally convolving the third difference with the reflection data; the convolution result is obtained by spatial integration of the convolution result and added to the negative of the uplink focusing function to obtain the downlink Green's function.

4. The subsurface imaging method based on the renormalized Marchenko method according to claim 1, characterized in that, The upward Green's function is calculated based on the upward and downward focusing functions, specifically including: reversing the reflected data on the time axis, calculating the upward focusing function and... The fourth difference of the downlink focusing function is obtained by temporal convolution of the fourth difference with the inverted reflection data on the time axis. The convolution result is obtained by spatial integration of the convolution result and added to the negative of the downlink focusing function to obtain the uplink Green's function.

5. The subsurface imaging method based on the renormalized Marchenko method according to claim 1, characterized in that, The background velocity model is the longitudinal wave velocity of the underground medium.

Images