Marchenko imaging method based on optimal integration range constraint

Abstract

The application discloses a Marchenko imaging method based on optimal integral range constraint and belongs to the technical field of Marchenko imaging, and comprises the following steps: S1, acquiring an initial imaging profile of an underground target region; S2, determining an underground structure dip angle field according to the initial imaging profile of the underground target region; S3, determining an optimal integral range of each imaging point according to the underground structure dip angle field; S4, determining a Marchenko imaging result of the imaging point with the optimal integral range as a constraint; and S5, obtaining an imaging profile result of the underground target region according to the Marchenko imaging result of all the imaging points with the optimal integral range as a constraint.

Description

Technical Field

[0001] This invention belongs to the field of Marchenko imaging technology, specifically relating to a Marchenko imaging method based on optimal integration range constraints. Background Technology

[0002] Given a macroscopic velocity model, the Marchenko imaging method reconstructs the uplink and downlink Green's functions from seismic data by solving the Marchenko equations, and then uses the reconstructed Green's functions to construct images. The Marchenko imaging method can effectively suppress migration artifacts caused by interlayer multiples in seismic data in the imaging domain.

[0003] Marchenko imaging primarily utilizes the direct wave components of the ascending and descending Green's functions to image subsurface structures. When the imaging point is on a reflecting interface, the travel times of the ascending and descending Green's function direct waves are equal, and their cross-correlation phase axis lies at the zero-time axis, contributing energy to the imaging point. When the imaging point is not on a reflecting interface, there is a time difference between the ascending and descending Green's function direct waves. Therefore, the cross-correlation phase axis of the direct waves deviates from the zero-time axis, theoretically contributing no energy to the imaging point. However, due to the limited frequency band of the Green's function direct waves, the frequency further decreases after correlation. Thus, even before the imaging point reaches the reflecting interface, the sidelobes of the cross-correlation phase axis contribute energy to the imaging point. Furthermore, for the same imaging point, the time difference between the ascending and descending Green's function direct waves at the far offset position is smaller than that at the near offset position. The wavelet at the far offset of the cross-correlation phase axis is closer to the zero-time axis, causing imaging points on non-reflecting interfaces to also be imaged, resulting in a decrease in the final imaging resolution. Based on the implementation process of Marchenko imaging conditions, the simplest way to solve the resolution reduction problem is to directly remove the far offset response in the cross-correlation gather of the direct wave, thereby reducing the contribution of the far offset of the cross-correlation gather to the energy of the imaging point and thus improving the resolution of the imaging result.

[0004] Limiting the integration range can effectively improve the shallow resolution of Marchenko imaging, but a small integration range cannot accurately image steep structures. Inappropriate selection of the integration range will inevitably lead to suboptimal imaging results, thus affecting subsequent geological interpretation.

[0005] To address the shortcomings of the aforementioned Marchenko imaging method, this invention proposes a Marchenko imaging method based on optimal spatial integration range constraints. Summary of the Invention

[0006] To address the issues of low resolution in shallow imaging when the spatial integration range is too large and poor imaging quality of steep structures when the spatial integration range is too small in the Marchenko imaging method, this invention proposes a Marchenko imaging method based on optimal integration range constraints.

[0007] The technical solution of this invention is: a Marchenko imaging method based on optimal integration range constraints, comprising the following steps:

[0008] S1. Obtain the initial imaging profile of the underground target area;

[0009] S2. Determine the dip field of the underground structure based on the initial imaging profile of the underground target area;

[0010] S3. Determine the optimal integration range for each imaging point based on the dip field of the underground structure;

[0011] S4. Determine the Marchenko imaging results with the imaging points constrained by the optimal integration range;

[0012] S5. Based on the Marchenko imaging results of all imaging points constrained by the optimal integration range, the imaging profile results of the underground target area are obtained.

[0013] Furthermore, S1 includes the following sub-steps:

[0014] S11. Collect seismic data of the underground target area and preprocess the seismic data;

[0015] S12. Based on the preprocessed seismic data, construct a depth domain velocity model using a grid tomography modeling algorithm;

[0016] S13. Based on the depth domain velocity model, determine the direct wave from each imaging point in the underground target area to the surface;

[0017] S14. Initialize the downlink focusing function using the direct wave;

[0018] S15. Based on the preprocessed seismic data and the initial downlink focusing function, the downlink focusing function and the uplink focusing function are solved iteratively.

[0019] S16. Calculate the downlink Green's function and the uplink Green's function based on the downlink focusing function and the uplink focusing function obtained from the iterative solution;

[0020] S17. Determine the initial imaging profile of the underground target area based on the downlink Green's function and the uplink Green's function.

[0021] Furthermore, in S15, the downlink focusing function f1 +The expression is:

[0022]

[0023] In the formula, Θ a Represents the first time window function, Θ b This represents the second time window function, K represents the iteration number, k represents the iteration number, and R represents the preprocessed seismic data. * This indicates a reversal in time. Indicates the initial downlink focusing function;

[0024] In S15, the upward focusing function f1 - The expression is:

[0025] f1 - =Θ a Rf1 + ;

[0026] In S16, the downlink Green's function G + (x R ,x A The expression for ,t) is:

[0027]

[0028] In the formula, Ψ a This represents the difference between 1 and the first time window function. Let x represent the range of integration. A Indicates the underground imaging point, x R Indicates the surface receiving point, x S t represents the surface artillery position, t represents the first time point, and t′ represents the second time point;

[0029] In S16, the ascending Green's function G - (x R ,x A The expression for (, -t) is:

[0030]

[0031] In the formula, Ψ b This represents the difference between 1 and the second time window function;

[0032] In S17, the initial imaging profile I(x) of the underground target area A The expression for ) is:

[0033] I(x A ) = D d (x A ,x A (t=0);

[0034]

[0035] In the formula, D d (·) represents the cross-correlation function for energy compensation, and ε represents the positive time constant. G represents the direct wave component of the descending Green's function. + (x R ,x A ,tt′) represents the downlink Green's function corresponding to the second time.

[0036] Furthermore, S2 includes the following sub-steps:

[0037] S21. Extract the offset imaging profile construction information based on the initial imaging profile;

[0038] S22. The offset imaging profile construction information is used as the reflection coefficient model. The reflection coefficient model is transformed to the time domain using the depth-time conversion method and then convolved with the Ricker wavelet.

[0039] S23. The convolution results are converted to the depth domain using the time-depth conversion method to obtain the synthetic seismic record;

[0040] S24. Based on synthetic seismic records, the similarity coefficient of the dip angle is calculated using the discrete similarity dip angle scanning algorithm, and the dip angle corresponding to the maximum value of the similarity coefficient is taken as the tectonic dip angle.

[0041] S25. Based on the structural dip angle, interpolation is performed to obtain the dip angle information of the entire profile, and the dip angle information of the entire profile is smoothed to obtain the underground structural dip angle field.

[0042] Furthermore, S3 includes the following sub-steps:

[0043] S31. Calculate the integral radius of each imaging point in the initial imaging profile of the underground target area;

[0044] S32. Based on the structural dip angle of the underground structural dip field where the imaging point is located, determine the intersection point of the imaging point normal and the ground surface, and take the intersection point as the center point of the optimal integration range.

[0045] S33. Determine the optimal integration range of the imaging point based on the center point of the optimal integration range.

[0046] Furthermore, in S31, the expression for the integral radius L of the imaging point is:

[0047]

[0048] In the formula, s represents the stretching coefficient of the wavelet, λ represents the seismic wavelength, and h represents the depth of the randomly selected imaging point;

[0049] In S33, the optimal integration range of the imaging point The expression is:

[0050]

[0051] In the formula, x R Indicates the surface receiving point, x O This represents the center point of the optimal integration range.

[0052] Furthermore, S4 includes the following sub-steps:

[0053] S41. Based on the optimal integration range of the imaging points, select the seismic data required to reconstruct the Green's function;

[0054] S42. Based on the seismic data required to reconstruct the Green's function, determine the downlink Green's function and the uplink Green's function constrained by the optimal integration range;

[0055] S43. Based on the downlink Green function and uplink Green function constrained by the optimal integration range, determine the Marchenko imaging result of the imaging point constrained by the optimal integration range.

[0056] The beneficial effects of this invention are:

[0057] (1) The Marchenko imaging method with optimal integration range constraint proposed in this invention can adaptively optimize the integration radius of the imaging point based on the travel time difference relationship of the uplink and downlink Green's function and the depth of the imaging point, thereby improving computational efficiency and reducing storage space requirements.

[0058] (2) The Marchenko imaging method with optimal integration range constraint proposed in this invention can effectively improve the shallow imaging resolution, while taking into account the wave field propagation direction at the steeply dipping structure location, thereby ensuring the imaging quality at the steeply dipping structure location. Attached Figure Description

[0059] Figure 1 The flowchart shows the Marchenko imaging method based on the optimal integration range constraint.

[0060] Figure 2 This is a schematic diagram of the velocity model;

[0061] Figure 3 A schematic diagram of Marchenko imaging results under the condition of fixed integration range;

[0062] Figure 4 This is a schematic diagram of the imaging structure;

[0063] Figure 5 A schematic diagram of the tilt angle for constructing the imaging position;

[0064] Figure 6 This is a schematic diagram of the dip angle field of the entire cross-section;

[0065] Figure 7 This is a schematic diagram of the Marchenko imaging profile under the constraint of optimal integration range. Detailed Implementation

[0066] The embodiments of the present invention will be further described below with reference to the accompanying drawings.

[0067] like Figure 1 As shown, this invention provides a Marchenko imaging method based on optimal integration range constraints, comprising the following steps:

[0068] S1. Obtain the initial imaging profile of the underground target area;

[0069] S2. Determine the dip field of the underground structure based on the initial imaging profile of the underground target area;

[0070] S3. Determine the optimal integration range for each imaging point based on the dip field of the underground structure;

[0071] S4. Determine the Marchenko imaging results with the imaging points constrained by the optimal integration range;

[0072] S5. Based on the Marchenko imaging results of all imaging points constrained by the optimal integration range, the imaging profile results of the underground target area are obtained.

[0073] First, the integration radius of the Marchenko image at each imaging point is determined based on the travel time relationship of the upper and lower Green's functions. Then, the center point of integration is determined based on the structural dip information of the imaging point. Finally, the final integration range of the Marchenko image is determined based on the integration center point and the integration radius. By using this method to provide the optimal integration range for imaging, the resolution of shallow imaging can be improved, while also compensating for the imaging energy of steep structures.

[0074] In this embodiment of the invention, S1 includes the following sub-steps:

[0075] S11. Collect seismic data of the underground target area and preprocess the seismic data;

[0076] S12. Based on the preprocessed seismic data, construct a depth domain velocity model using a grid tomography modeling algorithm;

[0077] S13. Based on the depth domain velocity model, determine the direct wave from each imaging point in the underground target area to the surface;

[0078] S14. Initialize the downlink focusing function using the direct wave;

[0079] S15. Based on the preprocessed seismic data and the initial downlink focusing function, the downlink focusing function and the uplink focusing function are solved iteratively.

[0080] S16. Calculate the downlink Green's function and the uplink Green's function based on the downlink focusing function and the uplink focusing function obtained from the iterative solution;

[0081] S17. Determine the initial imaging profile of the underground target area based on the downlink Green's function and the uplink Green's function.

[0082] Mesh tomography is a seismic velocity modeling method based on ray travel time or waveform data. By dividing the subsurface medium into regular meshes and iteratively optimizing the velocity values ​​of each mesh cell, it reconstructs the subsurface velocity structure and achieves accurate characterization of complex structures.

[0083] In S11, preprocessing includes denoising, amplitude compensation, and static correction.

[0084] In this embodiment of the invention, in S15, the downlink focusing function f1 + The expression is:

[0085]

[0086] In the formula, Θ a Represents the first time window function, Θ b This represents the second time window function, K represents the iteration number, k represents the iteration number, and R represents the preprocessed seismic data. * This indicates a reversal in time. This represents the initial downlink focusing function; the number of iterations is typically set to 8.

[0087] The expression for the first time window function is:

[0088] Θ a (x R ,x A ,t)=θ(t d -ε-t);

[0089] In the formula, θ(·) represents the unit step function, ε represents a positive time constant, usually half the duration of the wavelet, and t d Indicates the initial arrival of the wave;

[0090] The expression for the second time window function is:

[0091] Θ b (x R ,x A ,t)=θ(t+t d -ε).

[0092] In S15, the upward focusing function f1 - The expression is:

[0093] f1 - =Θ a Rf1 + ;

[0094] In S16, the downlink Green's function G + (x R ,x A The expression for ,t) is:

[0095]

[0096] In the formula, Ψ a This represents the difference between 1 and the first time window function. Let x represent the range of integration. A Indicates the underground imaging point, x R Indicates the surface receiving point, x S t represents the surface artillery position, t represents the first time point, and t′ represents the second time point;

[0097] In S16, the ascending Green's function G - (x R ,x A The expression for (, -t) is:

[0098]

[0099] In the formula, Ψ b This represents the difference between 1 and the second time window function;

[0100] In S17, the initial imaging profile I(x) of the underground target area A The expression for ) is:

[0101] I(x A ) = D d (x A ,x A (t=0);

[0102]

[0103] In the formula, D d (·) represents the cross-correlation function for energy compensation, and ε represents the positive time constant. G represents the direct wave component of the descending Green's function. + (x R ,x A ,tt′) represents the downlink Green's function corresponding to the second time.

[0104] In this embodiment of the invention, S2 includes the following sub-steps:

[0105] S21. Extract the offset imaging profile construction information based on the initial imaging profile;

[0106] S22. The offset imaging profile construction information is used as the reflection coefficient model. The reflection coefficient model is transformed to the time domain using the depth-time conversion method and then convolved with the Ricker wavelet.

[0107] S23. The convolution results are converted to the depth domain using the time-depth conversion method to obtain the synthetic seismic record;

[0108] S24. Based on synthetic seismic records, the similarity coefficient of the dip angle is calculated using the discrete similarity dip angle scanning algorithm, and the dip angle corresponding to the maximum value of the similarity coefficient is taken as the tectonic dip angle.

[0109] S25. Based on the structural dip angle, interpolation is performed to obtain the dip angle information of the entire profile, and the dip angle information of the entire profile is smoothed to obtain the underground structural dip angle field.

[0110] The discrete similarity tilt angle scanning algorithm determines the optimal tilt angle by scanning a series of discrete tilt angle directions within a local region and evaluating the similarity measure of the data in each direction.

[0111] In S21, the offset imaging profile construction information is extracted by constructing and interpreting the initial imaging results. The value at the imaging construction location is set to 0.5, and the value at the non-imaging construction location is set to 0.0.

[0112] In S24, in order to reduce the amount of calculation and reduce the calculation error, the dip angle is only calculated for the structural positions, and the dip angle is not calculated for non-structural positions.

[0113] In this embodiment of the invention, S3 includes the following sub-steps:

[0114] S31. Calculate the integral radius of each imaging point in the initial imaging profile of the underground target area;

[0115] S32. Based on the structural dip angle of the underground structural dip field where the imaging point is located, determine the intersection point of the imaging point normal and the ground surface, and take the intersection point as the center point of the optimal integration range.

[0116] S33. Determine the optimal integration range of the imaging point based on the center point of the optimal integration range.

[0117] In this embodiment of the invention, in S31, the expression for the integral radius L of the imaging point is:

[0118]

[0119] In the formula, s represents the stretching factor of the wavelet, λ represents the seismic wavelength, and h represents the depth of the randomly selected imaging point; the stretching factor of the wavelet is a given positive constant, usually given as 1.2.

[0120] In S33, the optimal integration range of the imaging point The expression is:

[0121]

[0122] In the formula, x R Indicates the surface receiving point, x O This represents the center point of the optimal integration range.

[0123] In this embodiment of the invention, S4 includes the following sub-steps:

[0124] S41. Based on the optimal integration range of the imaging points, select the seismic data required to reconstruct the Green's function;

[0125] S42. Based on the seismic data required to reconstruct the Green's function, determine the downlink Green's function and the uplink Green's function constrained by the optimal integration range;

[0126] S43. Based on the downlink Green function and uplink Green function constrained by the optimal integration range, determine the Marchenko imaging result of the imaging point constrained by the optimal integration range.

[0127] The present invention will now be described with reference to specific embodiments.

[0128] Taking the test model as an example, the first step is to input the collected seismic data and preprocess it, including denoising, amplitude compensation, and static correction. A depth-domain velocity model is then established using a mesh tomography modeling algorithm, such as... Figure 2 As shown.

[0129] The second step involves using the Marchenko imaging method with a fixed spatial integration range to obtain an initial imaging profile of the underground target area, such as... Figure 3 As shown.

[0130] The third step is to extract structural information based on the initial imaging profile, such as... Figure 4 As shown.

[0131] Step four, as Figure 5 As shown, the dip angle of the underground structure is obtained based on the structural information, and the structural dip field of the entire profile is calculated by interpolation, as follows. Figure 6 As shown.

[0132] The fifth step is to select an imaging point, determine the integration radius of the imaging point based on the time difference function, determine the integration center point of the imaging point based on the constructed tilt field, and determine the optimal integration range by combining the integration radius and the integration center point.

[0133] The fifth step is to calculate the upward and downward Green's functions under the optimal integration range constraint.

[0134] The sixth step is to calculate the Marchenko imaging values ​​using the reconstructed uplink and downlink Green's functions.

[0135] Step 7: Determine if this is the last imaging point. If not, proceed to step 5 to image the next imaging point. If it is, output the Marchenko imaging profile, as shown below. Figure 7 As shown.

[0136] Those skilled in the art will recognize that the embodiments described herein are intended to help the reader understand the principles of the invention, and should be understood that the scope of protection of the invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific modifications and combinations based on the technical teachings disclosed in this invention without departing from the spirit of the invention, and these modifications and combinations are still within the scope of protection of this invention.

Claims

1. A Marchenko imaging method based on optimal integration range constraints, characterized in that, Includes the following steps: S1. Obtain the initial imaging profile of the underground target area; S2. Determine the dip field of the underground structure based on the initial imaging profile of the underground target area; S3. Determine the optimal integration range for each imaging point based on the dip field of the underground structure; S4. Determine the Marchenko imaging results with the imaging points constrained by the optimal integration range; S5. Based on the Marchenko imaging results of all imaging points constrained by the optimal integration range, the imaging profile results of the underground target area are obtained. S2 includes the following sub-steps: S21. Extract the offset imaging profile construction information based on the initial imaging profile; S22. The offset imaging profile construction information is used as the reflection coefficient model. The reflection coefficient model is transformed to the time domain using the depth-time conversion method and then convolved with the Ricker wavelet. S23. The convolution results are converted to the depth domain using the time-depth conversion method to obtain the synthetic seismic record; S24. Based on synthetic seismic records, the similarity coefficient of the dip angle is calculated using the discrete similarity dip angle scanning algorithm, and the dip angle corresponding to the maximum value of the similarity coefficient is taken as the tectonic dip angle. S25. Based on the structural dip angle, interpolation is performed to obtain the dip angle information of the entire profile, and the dip angle information of the entire profile is smoothed to obtain the underground structural dip angle field. S3 includes the following sub-steps: S31. Calculate the integral radius of each imaging point in the initial imaging profile of the underground target area; S32. Based on the structural dip angle of the underground structural dip field where the imaging point is located, determine the intersection point of the imaging point normal and the ground surface, and take the intersection point as the center point of the optimal integration range. S33. Determine the optimal integration range of the imaging point based on the center point of the optimal integration range; In S31, the expression for the integral radius L of the imaging point is: ； In the formula, s represents the stretching coefficient of the wavelet, λ represents the seismic wavelength, and h represents the depth of the randomly selected imaging point; In S33, the optimal integration range of the imaging point The expression is: ； In the formula, x R Indicates the surface receiving point, x O This represents the center point of the optimal integration range.

2. The Marchenko imaging method based on optimal integration range constraints according to claim 1, characterized in that, S1 includes the following sub-steps: S11. Collect seismic data of the underground target area and preprocess the seismic data; S12. Based on the preprocessed seismic data, construct a depth domain velocity model using a grid tomography modeling algorithm; S13. Based on the depth domain velocity model, determine the direct wave from each imaging point in the underground target area to the surface; S14. Initialize the downlink focusing function using the direct wave; S15. Based on the preprocessed seismic data and the initial downlink focusing function, the downlink focusing function and the uplink focusing function are solved iteratively. S16. Calculate the downlink Green's function and the uplink Green's function based on the downlink focusing function and the uplink focusing function obtained from the iterative solution; S17. Determine the initial imaging profile of the underground target area based on the downlink Green's function and the uplink Green's function.

3. The Marchenko imaging method based on optimal integration range constraints according to claim 2, characterized in that, In S15, the downlink focusing function The expression is: ； In the formula, Θ a Represents the first time window function, Θ b This represents the second time window function, K represents the iteration number, k represents the iteration number, and R represents the preprocessed seismic data. * This indicates a reversal in time. Indicates the initial downlink focusing function; In S15, the upward focusing function The expression is: ； In S16, the downlink Green's function G + (x R ,x A The expression for ,t) is: ； In the formula, Ψ a This represents the difference between 1 and the first time window function. Let x represent the range of integration. A Indicates the underground imaging point, x R Indicates the surface receiving point, x S This represents a surface artillery point, and t represents the first time. Indicates the second time; In S16, the upward Green's function G - (x R ,x A The expression for (, -t) is: ； In the formula, Ψ b This represents the difference between 1 and the second time window function; In S17, the initial imaging profile I(x) of the underground target area A The expression for ) is: ； ； In the formula, D d (·) represents the cross-correlation function for energy compensation, and ε represents the positive time constant. This represents the direct wave portion of the descending Green's function. This represents the downlink Green's function corresponding to the second time step.

4. The Marchenko imaging method based on optimal integration range constraints according to claim 1, characterized in that, S4 includes the following sub-steps: S41. Based on the optimal integration range of the imaging points, select the seismic data required to reconstruct the Green's function; S42. Based on the seismic data required to reconstruct the Green's function, determine the downlink Green's function and the uplink Green's function constrained by the optimal integration range; S43. Based on the downlink Green function and uplink Green function constrained by the optimal integration range, determine the Marchenko imaging result of the imaging point constrained by the optimal integration range.

Images