Underground imaging method based on reformed Marchenko method

Abstract

The invention belongs to the technical field of earth exploration, and relates to an underground imaging method based on a reformed Marchenko method, which comprises the following steps: removing direct waves from seismic data, and carrying out seismic deconvolution to obtain preprocessed seismic data; selecting an underground imaging point, setting the underground imaging point as an underground focus point, and estimating a direct wave Green function from the underground focus point to a ground surface receiving point according to a background speed model; the direct wave Green function is reversed on a time axis to serve as an initial downlink focusing function, and an initial uplink focusing function is set to be zero; and performing iterative calculation according to the initial downlink focusing function and the initial uplink focusing function to obtain a downlink focusing function and an uplink focusing function after iteration, calculating a downlink green function and an uplink green function according to the uplink focusing function and the downlink focusing function, and imaging the underground imaging point. According to the method, iteration absolute convergence can be ensured by mediating convergence parameters, and finally an underground imaging result which is not influenced by multiple waves is obtained.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of exploring the earth, and particularly relates to a method for imaging the underground based on a renormalized Marchenko method. BACKGROUND

[0002] In the field of exploration seismology, the use of the Marchenko method to calculate the underground Green function was proposed in 2014. The Marchenko method is a method for calculating the Green function of a virtual point source in the underground that is excited and received at the ground surface using surface artificial seismic excitation and reception data, which can effectively process the multiple waves between the underground layers. The calculated underground Green function contains correct multiple wave information compared with the traditional calculation method proposed in the 1980s and 1990s. The Marchenko method is derived from the reciprocity theorem of seismology to obtain the Marchenko equation, and the focusing function is calculated by solving the Marchenko equation through the Neumann series iteration. The underground Green function is calculated from the focusing function according to the Marchenko method. The Neumann series iteration needs to meet certain convergence conditions. For seismic data that does not meet the convergence conditions, the Green function cannot be normally iterated and converged. The seismic imaging technology based on the Marchenko method can effectively suppress the multiple wave artifacts in the underground imaging results. In particular, when facing seismic data of the underground strong scattering medium containing surface-related multiple waves, the conventional Marchenko method almost cannot converge, and cannot suppress the multiple wave artifacts. SUMMARY

[0003] The embodiment of the application provides a method for imaging the underground based on a renormalized Marchenko method, which solves the problem of non-convergence of the conventional Marchenko method when facing seismic data.

[0004] According to the method for seismic imaging based on the renormalized Marchenko method provided in the embodiment of the application, the method comprises the following steps: The direct wave is removed from the seismic data, and seismic deconvolution is performed to obtain preprocessed seismic data; An underground imaging point is selected, and the underground imaging point is set as an underground focusing point. The direct wave Green function from the underground focusing point to the ground surface receiving point is estimated according to a background velocity model; The direct wave Green function is inverted on the time axis as an initial downward focusing function, and the initial upward focusing function is set to zero. Iterative calculation is performed according to the initial downward focusing function and the initial upward focusing function to obtain the iterated downward focusing function and the upward focusing function. The iterative calculation comprises the following steps: The downward focusing function of the last iteration is calculated, and the downward focusing function of the last iteration is subtracted from the initial downward focusing function to obtain the first iteration downward focusing function; times the first difference of the uplink focusing function; calculating the time convolution of the reflection data and the first difference to obtain the convolution result; performing spatial integration operation on the convolution result, multiplying it with the time window function and the convergence parameter, and then multiplying it with the uplink focusing function of the previous iteration. Add them up and get the new iterative upward focusing function. is the convergence parameter, is the surface reflection coefficient; Calculate the upward focusing function of the last iteration and The second difference of the downlink focusing function is times the second difference, and the time convolution of the reflected data inverted on the time axis and the second difference is calculated to obtain the convolution result; the convolution result is spatially integrated and multiplied with the time window function and the convergence parameter, and then multiplied with the downlink focusing function of the previous iteration. times and add them together with the initial downward focusing function Multiple superposition is performed to obtain the downward focusing function of a new iteration; Calculating the downward Green's function and the upward Green's function according to the upward focusing function and the downward focusing function; The underground imaging points are imaged according to the upward focusing function, the downward focusing function, the downward Green's function and the upward Green's function.

[0005] Furthermore, the time window function is obtained by determining the travel time of the seismic wave, which is expressed as: ,in, is half the time length of the seismic wavelet, is the travel time of the seismic wave, is the absolute value of time, represents the surface receiving point, represents the underground focal point, is the time window function.

[0006] Further, the downlink Green function is calculated based on the uplink focusing function and the downlink focusing function, specifically including: calculating the downlink focusing function and The third difference of the uplink focusing function is times the third difference, and the third difference is temporally convolved with the reflection data to obtain the convolution result; the convolution result is spatially integrated and added to the negative of the uplink focusing function to obtain the downlink Green's function.

[0007] Furthermore, the uplink Green's function is calculated based on the uplink focusing function and the downlink focusing function, specifically including: reversing the reflection data on the time axis, calculating the uplink focusing function and the downlink Green's function. The fourth difference of the downlink focusing function is times the fourth difference, and the fourth difference is time-convolved with the reflection data after inversion on the time axis to obtain the convolution result; the convolution result is spatially integrated and added to the negative of the downlink focusing function to obtain the uplink Green's function.

[0008] Further, the first imaging result is obtained by double focusing imaging using the uplink focusing function, the downlink focusing function and the downlink Green function, specifically including: calculating the fifth difference value of the downlink focusing function and the uplink focusing function with a factor of 2; performing time convolution on the fifth difference value and the downlink Green function to obtain a convolution result; performing spatial integral operation on the convolution result to obtain the underground imaging result.

[0009] Further, the second imaging result is obtained by cross-correlation imaging of the uplink Green function and the downlink Green function.

[0010] Further, the background velocity model is the longitudinal wave velocity of the underground medium.

[0011] Compared with the prior art, the application has the beneficial effects that the application breaks through the convergence condition of the conventional Marchenko algorithm, the conventional Marchenko algorithm cannot iteratively converge the seismic data containing surface-related multiple waves, and the method of the application can ensure absolute convergence of iteration by adjusting the convergence parameter, and finally obtain the underground imaging result not affected by multiple waves. BRIEF DESCRIPTION OF DRAWINGS

[0012] Figure 1 A method flowchart of an underground imaging method based on a renormalization Marchenko method provided by an embodiment of the application; Figure 2 A schematic diagram of a velocity model and a density model provided by an embodiment of the application; Figure 3 A schematic diagram of reflection data provided by an embodiment of the application; Figure 4 A schematic diagram of a background velocity model provided by an embodiment of the application; Figure 5 A drawing provided by an embodiment of the application, in which (a) is an initial focusing function, (b) is a downlink focusing function and (c) is an uplink focusing function corresponding to a focusing imaging point; Figure 6 A drawing provided by an embodiment of the application, in which (a) is a Green function, (b) is an uplink Green function and (c) is a downlink Green function; Figure 7 An underground imaging result diagram provided by an embodiment of the application. DETAILED DESCRIPTION

[0013] In order to make the purpose, technical scheme and advantages of the application clearer, the application will be further described in detail below with reference to the embodiments. It should be understood that the specific embodiments described herein are only used to explain the application and not to limit the application. ​

[0014] See also Figure 1 A seismic imaging method based on the renormalized Marchenko method is shown, comprising: S1 removes direct waves from seismic data and performs seismic deconvolution to obtain preprocessed seismic data; Seismic data contains free surface multiple waves, which are data collected by blasting on the surface to stimulate surface reception. Direct waves are waves that propagate directly from the earthquake source to the surface receiving point without reflection or refraction at any interface. They propagate along the shortest path, usually a straight line. Removing direct waves from seismic data usually refers to removing the direct waves with strong energy near the earthquake source. The preprocessing process is the denoising process. After removing the direct waves with strong energy near the earthquake source, the remaining seismic data is subjected to seismic deconvolution. Seismic deconvolution is a conventional processing method that eliminates the influence of seismic wavelets from seismic data and restores the original form of the underground surface reflection coefficient, thereby improving the resolution and signal-to-noise ratio of the seismic data. The preprocessed seismic data is recorded as: , Represents the meaning of reflection data. The following reflection data refers to the pre-processed seismic data. represents the free boundary condition of the surface. represents the surface receiving point, represents the coordinates of the excitation point, Represents time. Includes reflections and multiples (interlayer multiples and surface-related multiples).

[0015] S2 selects an underground imaging point and sets it as an underground focusing point, and estimates the direct wave Green's function from the underground focusing point to the surface receiving point based on the background velocity model; For the underground area to be imaged, a background velocity model needs to be prepared and recorded as , Represents depth, represents the lateral distance, and the background velocity model is the longitudinal wave velocity of the underground medium. The accuracy requirements are consistent with those of conventional imaging methods, that is, a smooth background velocity model is sufficient, which can be obtained through travel-time tomography or full waveform inversion technology in conventional seismic data processing procedures.

[0016] Select an underground imaging point and set it as the underground focal point. The coordinates are expressed as According to a background velocity model Estimating underground focal points To the surface receiving point The direct wave Green function , subscript Represents the direct wave component. Direct wave Green's function Green function of direct wave propagating from a subsurface focus point to a surface receiver.

[0017] The direct wave Green function from a subsurface focus point to a surface receiver is estimated according to a background velocity model, which can be obtained by solving wave equation directly or calculating travel time by solving eikonal equation The travel time curve is obtained by convolving the reflection data with the wavelet.

[0018] S3 reverses the direct wave Green function on the time axis as the initial down-going focusing function, and sets the initial up-going focusing function as zero; S3 reverses the direct wave Green function on the time axis as the initial down-going focusing function, and sets the initial up-going focusing function as zero; Herein represents the focusing function, and the superscript represents the down-going wave field, and the subscript represents the iteration number, and the subscript represents the iteration number is zero, i.e. the initial down-going focusing function, and represents the direct part of the focusing function, and the direct part of the down-going focusing function is the same as the initial down-going focusing function . The initial up-going focusing function is directly set as zero, i.e. , and the superscript represents the up-going wave field; S4 performs iterative calculation according to the initial down-going focusing function and the initial up-going focusing function to obtain the down-going focusing function and the up-going focusing function after iteration, and the iterative calculation comprises: calculating a first difference value of the down-going focusing function after the last iteration and times the up-going focusing function; calculating a convolution result by time-convolving the reflection data and the first difference value; multiplying the convolution result after spatial integration operation with a time window function and a convergence parameter, and adding the product to times the up-going focusing function after the last iteration to obtain the up-going focusing function of the new iteration, is the convergence parameter; calculating a second difference value of the up-going focusing function after the last iteration and times the down-going focusing function; calculating a convolution result by time-convolving the reflection data reversed on the time axis and the second difference value; multiplying the convolution result after spatial integration operation with a time window function and a convergence parameter, and adding the product to times the down-going focusing function after the last iteration, and adding the product to times the initial down-going focusing function to obtain the down-going focusing function of the new iteration; S5 calculates the down-going Green function and the up-going Green function according to the up-going focusing function and the down-going focusing function; S6 images underground imaging points based on the up-focusing function, down-focusing function, down-Green's function, and up-Green's function. In one embodiment, an iterative calculation is performed based on the initial downlink focusing function and the initial uplink focusing function, and the calculation formulas of the iterative downlink focusing function and the uplink focusing function are expressed as follows: , in, The number of iterations is The upward focusing function of The number of iterations is -1 upward focusing function, is the convergence parameter, is the time window function, Represents the surface, and the integration process is along the surface Integrate the excitation samples, The time variable is The reflection data at Indicates the number of iterations , time is The downward focusing function, is the surface reflection coefficient, which reflects the reflection effect of the free surface and is generally set to -1. The number of iterations is The upward focusing function of The number of iterations is The downward focusing function, The number of iterations is The downward focusing function, The time variable is The reflection data at .

[0019] In this embodiment of the present application, the convergence parameters , which changes the calculation steps of each iteration.

[0020] In one embodiment, the time window function is obtained by determining the travel time of the seismic wave, which is expressed as: ,in, is half the time length of the seismic wavelet, is the travel time of the seismic wave, is the absolute value of time, represents the surface receiving point, represents the underground focal point, is the time window function, and Refers to the same meaning.

[0021] In order to improve the operation speed, the convolution calculation can use fast Fourier transform to convert the two to the frequency domain, perform multiplication operation, and then perform inverse Fourier transform to convert back to the time domain, so that the calculation is equivalent to the direct time convolution calculation.

[0022] By adjusting the convergence parameter The iteration process can be ensured to converge absolutely, and the convergence parameter is between 0 and 1, and the smaller the value is, the slower the convergence speed is, and reducing the convergence parameter can improve the convergence performance while reducing the convergence speed, and is set according to the data, to ensure the convergence of the result at a higher convergence speed.

[0023] In an embodiment, the downgoing Green function is calculated according to the upgoing focusing function and the downgoing focusing function, specifically including: calculating a third difference value of the downgoing focusing function and the upgoing focusing function times the upgoing focusing function, performing time convolution on the third difference value and the reflection data to obtain a convolution result, and adding the convolution result after spatial integration operation and the negative of the upgoing focusing function to obtain the downgoing Green function.

[0024] The upgoing Green function is calculated according to the upgoing focusing function and the downgoing focusing function, specifically including: inverting the reflection data on the time axis, calculating a fourth difference value of the upgoing focusing function and the downgoing focusing function times the downgoing focusing function, performing time convolution on the fourth difference value and the reflection data inverted on the time axis to obtain a convolution result, and adding the convolution result after spatial integration operation and the negative of the downgoing focusing function to obtain the upgoing Green function.

[0025] After several iterations are stable, the upgoing Green function and the downgoing Green function are calculated according to the obtained downgoing focusing function and the upgoing focusing function, the subscript of the Green function indicates that the free surface multiple waves are included: , , The above calculation process does not need iteration and does not have a time window function, is an expression of the upgoing focusing function obtained after iteration when integration is performed, is an expression of the downgoing focusing function obtained after iteration when integration is performed, is an expression of the downgoing focusing function obtained after iteration, is an expression of the downgoing focusing function obtained after iteration, indicates a time variable.

[0026] In an embodiment, a first imaging result is obtained by using the upgoing focusing function, the downgoing focusing function and the downgoing Green's function for double focusing imaging, including: calculating a fifth difference value of the downgoing focusing function and times the upgoing focusing function; performing time convolution on the fifth difference value and the downgoing Green's function to obtain a convolution result; performing spatial integral operation on the convolution result to obtain an underground imaging result.

[0027] The double focusing imaging formula used is: , By adding the free surface correction term, i.e. , the imaging value at the underground focusing point can be obtained, and all imaging points underground are imaged, so that the underground imaging result can be obtained.

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

[0029] Taking a model test as an example, the velocity model and the density model shown in Figure 2 are used as test models to illustrate the implementation process of the present application, wherein Figure 2 (a) is the velocity model, Figure 2 (b) is the density model, and Figure 3 is a schematic diagram of reflection data, wherein the horizontal coordinate represents the surface receiving point , the receiving point interval is , the coordinate range is from to , a total of channel receiving data, the vertical coordinate represents time , the sampling interval is , the time is from to , a total of sampling points. Figure 3 The reflection data shown is one shot reflection data, i.e. the reflection data at , in addition, it also contains data from to at intervals , a total of shots of data; Prepare the background velocity model as shown in Figure 4 ; Select the underground focusing point, for example, the point at a horizontal distance , a depth , i.e. According to the background velocity model, the travel time is obtained by solving the Eikonal equation , the travel time curve is convolved with the Ricker wavelet, and the time axis is reversed to obtain the initial downward focusing function. The time range of the initial downward focusing function is set to arrive , the sampling interval is , the surface receiving point coordinates are from arrive , the interval is ,like Figure 5 As shown in (a), and using the travel time Set the time window function; Time of the time window function Also set the range from arrive , the sampling interval is , and the surface receiving point coordinates from arrive , the interval is ; Iteratively solve the converged downward focusing function and upward focusing function, and set the convergence parameters , which can ensure the result convergence, iteration The results after the second Figure 5 (b) shows the plot of the downlink focusing function, and Figure 5 (c) in Figure 3 shows the plot of the upward focusing function. It can be seen that the focusing function results are stable.

[0030] Find the upgoing Green's function and the downgoing Green's function as Figure 6 (b) and Figure 6 As shown in (c) in the figure, the Green function can be obtained by summing the uplink Green function and the downlink Green function. The result is also in Figure 6 As shown in (a), it can be seen that the Green function result is stable.

[0031] Calculate the imaging value of underground imaging points; Repeat the operation until all the required imaging points in the underground are calculated. Here, the horizontal grid spacing is used. Spacing from vertical grid , set the imaging grid to ,total imaging points, and obtain the final underground imaging result based on the renormalized Marchenko method, such as Figure 7 As shown in the figure, the imaging process is completely stable without any non-convergence, and the multiple wave artifacts in the final imaging results are suppressed.

[0032] The above only describes preferred embodiments of the present application and is not used to limit the present application, and any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application should be included in the protection scope of the present application.

Claims

1. A subsurface imaging method based on the renormalized Marchenko method, characterized in that: include: Remove direct waves from seismic data and perform seismic deconvolution to obtain pre-processed seismic data; Selecting an underground imaging point and setting the underground imaging point as an underground focusing point, estimating a direct wave Green's function from the underground focusing point to a surface receiving point based on a background velocity model; The direct wave Green's function is inverted on the time axis as the initial downward focusing function, and the initial upward focusing function is set to zero; An iterative calculation is performed based on the initial downlink focusing function and the initial uplink focusing function to obtain an iterated downlink focusing function and an iterative uplink focusing function, wherein the iterative calculation includes: Calculate the downward focusing function of the last iteration and times the first difference of the uplink focusing function; calculating the time convolution of the reflection data and the first difference to obtain the convolution result; performing spatial integration operation on the convolution result, multiplying it with the time window function and the convergence parameter, and then multiplying it with the uplink focusing function of the previous iteration. Add them up and get the new iterative upward focusing function. is the convergence parameter, is the surface reflection coefficient; Calculate the upward focusing function of the last iteration and The second difference of the downlink focusing function is times the second difference, and the time convolution of the reflected data inverted on the time axis and the second difference is calculated to obtain the convolution result; the convolution result is spatially integrated and multiplied with the time window function and the convergence parameter, and then multiplied with the downlink focusing function of the previous iteration. times and add them together with the initial downward focusing function Multiple superposition is performed to obtain the downward focusing function of a new iteration; Calculating the downward Green's function and the upward Green's function according to the upward focusing function and the downward focusing function; The underground imaging points are imaged according to the upward focusing function, the downward focusing function, the downward Green's function and the upward Green's function.

2. The underground imaging method based on the renormalized Marchenko method according to claim 1, characterized in that: The time window function is obtained by judging the travel time of the seismic wave, which is expressed as: ,in, is half the time length of the seismic wavelet, is the travel time of the seismic wave, is the absolute value of time, represents the surface receiving point, represents the underground focal point, is the time window function.

3. The underground imaging method based on the renormalized Marchenko method according to claim 1, characterized in that: The downlink Green function is calculated based on the uplink focusing function and the downlink focusing function, specifically including: calculating the downlink focusing function and The third difference of the uplink focusing function is times the third difference, and the third difference is temporally convolved with the reflection data to obtain the convolution result; the convolution result is spatially integrated and added to the negative of the uplink focusing function to obtain the downlink Green's function.

4. The underground imaging method based on the renormalized Marchenko method according to claim 1, characterized in that: The uplink Green's function is calculated based on the uplink focusing function and the downlink focusing function, specifically including: reversing the reflection data on the time axis, calculating the uplink focusing function and the downlink Green's function. The fourth difference of the downlink focusing function is times the fourth difference, and the fourth difference is time-convolved with the reflection data after inversion on the time axis to obtain the convolution result; the convolution result is spatially integrated and added to the negative of the downlink focusing function to obtain the uplink Green's function.

5. The underground imaging method based on the renormalized Marchenko method according to claim 1, characterized in that: The first imaging result is obtained by performing dual focusing imaging using an uplink focusing function, a downlink focusing function, and a downlink Green's function, specifically including: Calculate the downlink focusing function and times the fifth difference of the upward focusing function; Performing time convolution on the fifth difference value and the descending Green's function to obtain a convolution result; The convolution result is spatially integrated to obtain the underground imaging result.

6. The underground imaging method based on the renormalized Marchenko method according to claim 1, characterized in that: The uplink Green's function and the downlink Green's function are cross-correlated and imaged to obtain a second imaging result.

7. The underground 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