Rapid Marchenko imaging method for undulating surface conditions

Abstract

The invention provides a rapid Marchenko imaging method for undulating surface conditions, and the method comprises the following steps: S1, inputting an earthquake shot record of a fixed surface and a depth domain undulating migration surface, and correcting the earthquake shot record to the migration surface through a static correction method; s2, inputting a depth domain seismic velocity model, and constructing an initial downlink focusing function from an underground target area imaging point to an undulating surface migration surface through a ray tracing algorithm; and S3, taking the shot record on the fluctuating offset surface and the initial downlink focusing function as input, and calculating uplink and downlink Green functions from the underground imaging point to the ground through an iterative method. And S4, applying an energy normalization cross-correlation imaging condition to the reconstructed uplink and downlink Green functions to obtain a Marchenko imaging result. According to the Marchenko imaging method, the problems that conventional Marchenko imaging is low in calculation efficiency and undulating surface conditions cannot be processed are solved, and the Marchenko imaging quality under the undulating surface conditions is ensured while the calculation efficiency is improved.

Description

Technical Field

[0001] The invention provides a fast Marchenko imaging method for undulating surface conditions, belonging to the technical field of imaging methods. Background Art

[0002] Traditional migration imaging methods, such as Kirchhoff migration and reverse time migration imaging, are based on the single scattering assumption. During the migration process, they can only image a single reflection wave, but cannot correctly distinguish and process multiple reflection waves between layers. They will mistakenly image the multiple waves between layers as single reflection waves, resulting in the appearance of migration artifacts related to the multiple waves between layers in the imaging results, thus affecting the imaging quality. In recent years, the Marchenko imaging method has been developed based on the inverse scattering theory. In the case of a known macroscopic velocity model, the Marchenko imaging method reconstructs the uplink and downlink Green's functions from the seismic data by solving the Marchenko equation, and then uses the reconstructed Green's functions for structural imaging. The Marchenko imaging method can effectively suppress the migration artifacts caused by multiple waves between layers in the seismic data in the imaging domain.

[0003] With the continuous improvement of the exploration and development of oil and gas resources in eastern my country and the increasing shortage of oil and gas resources, simple structural exploration can no longer meet industrial needs. The seismic exploration area has gradually shifted from the east to the west where exploration is more difficult. The surface conditions have changed from nearly horizontal surface conditions to complex undulating surface conditions, such as mountain piedmont belts, hills, loess plateaus, gravel Gobi and other exploration areas. However, these areas have brought great difficulties and challenges to seismic data processing due to the violent surface fluctuations, complex near-surface structures and large lateral velocity changes.

[0004] At present, the research and application of Marchenko imaging methods are mainly based on the assumption of a horizontal surface. In actual data applications, especially for seismic data on the undulating surface in the western part of the land, it is necessary to correct the shot points and detection points of the seismic data to the same reference plane before Green's function reconstruction and imaging processing. When the surface elevation and near-surface lateral velocity change greatly, the use of static correction to correct the shot points of the seismic data to the same reference plane will inevitably destroy the hyperbolic shape of the reflected wave and distort the seismic wave field, resulting in a decrease in the quality of Marchenko imaging after elevation static correction, and the imaging accuracy cannot meet the needs of geological interpretation.

[0005] In addition, the Marchenko imaging method uses a single underground imaging point as a unit, and reconstructs and images the Green's function of each imaging point in the underground target area. A necessary step in this process is to construct an initial downward focusing function from each imaging point underground to the ground. The conventional approach is to use a finite difference forward modeling method or a ray tracing method to calculate the first arrival wave from the imaging point to the ground through a background velocity model. However, this process is very time-consuming, especially when processing three-dimensional data, which results in extremely low calculation efficiency of Marchenko imaging. Summary of the invention

[0006] In view of the shortcomings of the above-mentioned Marchenko imaging method, the present invention proposes a fast Marchenko imaging technology for undulating surface conditions, which aims to calculate the travel time volume from the undulating surface to the underground target area by using a fast ray tracing algorithm in a step-by-step manner, and then extract the initial arrival travel time from the underground imaging point to the undulating surface, and then use the travel time and seismic wavelet to construct the initial downward focusing function from the imaging point to the undulating surface, and finally perform Green's function reconstruction and Marchenko imaging. This technology not only improves the Marchenko imaging accuracy under undulating surface conditions, but also greatly improves the computational efficiency of Marchenko imaging.

[0007] The technical problems to be solved by the present invention are:

[0008] (1) Solve the problem of introducing errors when Marchenko imaging faces undulating surface conditions and using static correction to correct the shot check points of seismic data to the same reference surface

[0009] (2) Solve the huge computational complexity problem caused by the need to calculate the travel time from the imaging point to the ground when constructing the initial downlink focusing function in Marchenko imaging.

[0010] The technical solution adopted by the present invention is a fast Marchenko imaging method for undulating surface conditions, comprising the following steps:

[0011] S1: Input the seismic shot records of the fixed surface and the depth domain fluctuation offset surface, and correct the seismic shot records to the offset surface through the static correction method.

[0012] S2: Input the depth domain seismic velocity model and construct the initial downward focusing function from the imaging point of the underground target area to the offset surface of the undulating surface through the ray tracing algorithm.

[0013] S3: Taking the shot records on the undulating offset surface and the initial down-going focusing function as input, the up-going and down-going Green's functions from the underground imaging point to the ground are calculated by an iterative method.

[0014] S4: The energy-normalized cross-correlation imaging condition is applied to the reconstructed uplink and downlink Green's functions to obtain the Marchenko imaging results.

[0015] Preferably, S1 includes the following sub-steps:

[0016] S11: Input data, including time domain seismic shot records and depth domain ups and downs offset surface. The time domain seismic shot records include the horizontal coordinates and elevations of the shot points and receiver points, and the depth domain ups and downs offset surface includes the horizontal coordinates and elevations of each sample point.

[0017] S12: Calculate the static correction value T of each seismic data according to the replacement speed from the offset surface to the fixed surface, the coordinates of the shot point detection points in the shot record, and the elevation of the undulating offset surface. statics , correction amount T statics The calculation formula is:

[0018]

[0019] Where sy and sx are the coordinates of the shot point; gy and gx are the coordinates of the detector point; E(sy, sx) is the elevation of the projected point of the shot point on the offset surface; E(gy, gx) is the elevation of the projected point of the detector point on the offset surface; v rep is the replacement speed between the undulating offset surface and the fixed surface.

[0020] S13: Calculate the number of sample points Nt of each seismic data upward translation according to the static correction value statics , and then perform static correction on the seismic data according to the number of sample points. The calculation formula for the number of translation sample points is:

[0021]

[0022] Where dt is the temporal sampling rate of the seismic data.

[0023] Preferably, S2 includes the following sub-steps:

[0024] S21: Count the shot point coordinates and receiver point coordinates of all the trace data in the shot record and determine the range of the coordinates. Then, according to the coordinate range, determine the coordinates of the source point of the ray tracing on the undulating offset surface. Finally, calculate the travel time from each source point to the imaging point in the underground target area. The travel time calculation formula is:

[0025]

[0026] Where s is the ray path; s0 is the initial position of the ray path; s i is the end position of the ray path, τ(s i ,s0) is from point s0 to point s i ’s travel time; v is the speed of the ray.

[0027] S22: The travel time τ(s0,s i ).

[0028] S23: Generate synthetic seismic wavelet, the calculation formula is:

[0029]

[0030] Among them, r(t) is the seismic wavelet expression; f p is the main frequency of the seismic wavelet; t is time; t0 is the time corresponding to the wavelet peak position; π is pi; exp is the natural exponential function;

[0031] S24: Using the travel time τ(s0,s i ) replaces t0 in the seismic wavelet r(t) to construct the initial downward focusing function

[0032] Preferably, S3 includes the following sub-steps:

[0033] S31: Calculate the complete downlink focusing function using the constructed initial downlink focusing function, the calculation formula is as follows:

[0034]

[0035] Among them, f1 + is the downward focusing function including the coda wave; R is the shot record corrected to the offset plane; K is the number of iterations; * indicates time reversal, and the superscript “+” indicates the downward direction; Θ a and θ b is the time window function, which is defined as follows:

[0036] Θ a (t) = θ(t d -ε-t) (6)

[0037] Θ b (t) = θ(t + t d -ε) (7)

[0038] Where, t is time; t d is the travel time of the first arrival wave, from τ(s0,s i ) is determined; ε is a positive time constant, usually half the duration of the sub-wave; θ(t) is the unit step function.

[0039] S32: Calculate the uplink focusing function according to the downlink focusing function calculated in the previous step, and the calculation formula is:

[0040] f1 - =Θa Rf1 + (8)

[0041] Among them, f1 - is the upward focusing function; the superscript “-” indicates the direction is upward.

[0042] S33: Calculate the uplink and downlink Green functions according to the calculated uplink and downlink focusing functions, and the calculation formula is:

[0043] G + =Ψ a Rf1 + (9)

[0044]

[0045] Among them, G + and G - are the downward and upward Green's functions respectively; Ψ a and b is the time window function, which is defined as follows:

[0046] Ψ a,b (t) = 1 - Θ a,b (t) (11)

[0047] Preferably, S4 includes the following sub-steps:

[0048] S41: Using the calculated downlink Green's function and uplink Green's function first arrival wave to obtain the cross-correlation function, the calculation formula is:

[0049]

[0050] Where, C(t) is the cross-correlation function; is the first arrival wave of the upward Green's function; t and t' are time; N represents the number of channels.

[0051] S42: Use the upgoing Green's function first arrival wave to calculate the amplitude compensation factor. The calculation formula is:

[0052]

[0053] S42: Perform amplitude compensation on the cross-correlation function, and take the value at time t=0 as the imaging result. The imaging result expression is:

[0054]

[0055] Where I is the Marchenko imaging result of the underground imaging point.

[0056] The invention solves the problem that conventional Marchenko imaging has low calculation efficiency and cannot process undulating surface conditions, and ensures the quality of Marchenko imaging under undulating surface conditions while improving the calculation efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] Figure 1 It is a flowchart of the present invention;

[0058] Figure 2 The model used for the embodiment is a velocity model of surface undulation;

[0059] Figure 3 It is a schematic diagram of a ray tracing path from a source point on an undulating surface to underground in an embodiment;

[0060] Figure 4 A schematic diagram of a ray tracing path from an underground imaging point extracted from a travel-time body to an undulating surface in an embodiment;

[0061] Figure 5 A seismic wavelet synthesized for an embodiment;

[0062] Figure 6 An initial downward focusing function from the underground imaging point to the ground constructed for the embodiment;

[0063] Figure 7 The upgoing Green's function from the underground imaging point to the undulating surface calculated for the embodiment;

[0064] Figure 8 The downward Green's function from the underground imaging point to the undulating surface calculated for the embodiment;

[0065] Fig. 9 This is the Marchenko imaging result of the target area in the embodiment. DETAILED DESCRIPTION

[0066] The specific technical solution of the present invention is explained in conjunction with embodiments.

[0067] The present invention will be further described below in conjunction with the accompanying drawings and specific embodiments:

[0068] Taking the test model as an example, a fast Marchenko imaging method for undulating surface conditions, such as Figure 1 As shown, the following steps are included:

[0069] The first step is to input fixed surface shot records, heaving offset surfaces, and velocity models, such as Figure 2 As shown, the seismic shot records are corrected to the offset plane by static correction method.

[0070] In the second step, the ray tracing algorithm is used to calculate the travel time volume from the offset surface to the imaging point of the underground target area, such as Figure 3As shown, the travel time table from the underground imaging point to the undulating offset surface is extracted from the travel time body, as shown in Figure 4 .

[0071] The third step is to synthesize seismic wavelets, such as Figure 5 shown.

[0072] The fourth step is to select an imaging point, extract the travel time from the point to the undulating offset surface, and then construct the initial downward focusing function from the underground imaging point to the undulating offset surface, such as Figure 6 shown.

[0073] The fifth step is to calculate the uplink Green's function and downlink focusing function from the underground imaging point to the ground, such as Figure 7 and Figure 8 .

[0074] The sixth step is to calculate the Marchenko imaging value using the reconstructed uplink Green's function and downlink focusing function.

[0075] Step 7: Determine whether it is the last imaging point. If it is not the last imaging point, proceed to step 4 to image the next imaging point. If it is, output the Marchenko imaging profile, such as Fig. 9 shown.

Claims

1. A fast Marchenko imaging method for undulating surface conditions, characterized in that: The steps include: S1: Input the seismic shot record of the fixed surface and the depth domain fluctuation offset surface, and correct the seismic shot record to the offset surface by static correction method; S2: Input the depth domain seismic velocity model and construct the initial downward focusing function from the imaging point of the underground target area to the offset surface of the undulating surface through the ray tracing algorithm; S3: Taking the shot records on the undulating offset surface and the initial up-going focusing function as input, the up-going and down-going Green's functions from the underground imaging point to the ground are calculated by an iterative method; S4: The energy-normalized cross-correlation imaging condition is applied to the reconstructed uplink and downlink Green's functions to obtain the Marchenko imaging results.

2. The fast Marchenko imaging method for undulating surface conditions according to claim 1, characterized in that: S1 includes the following sub-steps: S11: Input data, including time domain seismic shot records and depth domain ups and downs offset surface; the time domain seismic shot records include the horizontal coordinates and elevations of the shot points and receiver points, and the depth domain ups and downs offset surface includes the horizontal coordinates and elevations of each sample point; S12: Calculate the static correction value T of each seismic data according to the replacement speed from the offset surface to the fixed surface, the coordinates of the shot point detection points in the shot record, and the elevation of the undulating offset surface. statics , correction amount T statics The calculation formula is: Where sy and sx are the coordinates of the shot point; gy and gx are the coordinates of the detector point; E(sy, sx) is the elevation of the projected point of the shot point on the offset surface; E(gy, gx) is the elevation of the projected point of the detector point on the offset surface; v rep is the replacement speed between the undulating offset surface and the fixed surface; S13: Calculate the number of sample points Nt of each seismic data upward translation according to the static correction value statics , and then perform static correction on the seismic data according to the number of sample points; the calculation formula for the number of translation sample points is: Where dt is the temporal sampling rate of the seismic data.

3. The fast Marchenko imaging method for undulating surface conditions according to claim 1, characterized in that: S2 includes the following sub-steps: S21: Count the shot point coordinates and receiver point coordinates of all the trace data in the shot record and determine the range of the coordinates. Then, according to the coordinate range, determine the coordinates of the source point of the ray tracing on the undulating offset surface. Finally, calculate the travel time from each source point to the imaging point in the underground target area. The travel time calculation formula is: Where s is the ray path; s0 is the starting position of the ray path; s i is the end position of the ray path, τ(s i ,s0) is the distance from point s0 to point s i The travel time of ; v is the speed of the ray; S22: The travel time τ(s0,s i ); S23: Generate synthetic seismic wavelet, the calculation formula is: Among them, r(t) is the seismic wavelet expression; f p is the main frequency of the seismic wavelet; t is time; t0 is the time corresponding to the wavelet peak position; π is pi; exp is the natural exponential function; S24: Using the travel time τ(s0,s i ) replaces t0 in the seismic wavelet r(t) to construct the initial downward focusing function 4. The fast Marchenko imaging method for undulating terrain conditions according to claim 1, characterized in that: S3 includes the following sub-steps: S31: Calculate the complete downlink focusing function using the constructed initial downlink focusing function, and the calculation formula is as follows: Among them, f1 + is the downward focusing function including the coda wave; R is the shot record corrected to the offset plane; K is the number of iterations; * indicates time reversal; the superscript "+" indicates the downward direction; Θ a and θ b is the time window function, which is defined as follows: I a (t)=θ(t d -e-t) (6) I b (t)=θ(t+t d -e) (7) Where, t is time; t d is the travel time of the first arrival wave, from τ(s0,s i ) is determined; ε is a positive time constant, usually half of the duration of the sub-wave; θ(t) is a unit step function; S32: Calculate the uplink focusing function according to the downlink focusing function calculated in the previous step, and the calculation formula is: f1 - =Θ a Rf1 + (8) Among them, f1 - is the upward focusing function; the superscript "-" indicates the direction is upward; S33: Calculate the uplink and downlink Green functions according to the calculated uplink and downlink focusing functions, and the calculation formula is: G + =Ψ a Rf1 + (9) Among them, G + and G - are the downward and upward Green's functions respectively; Ψ a and b is the time window function, which is defined as follows: P a,b (t)=1-Θ a,b (t) (11).

5. The fast Marchenko imaging method for undulating terrain conditions according to claim 1, characterized in that: S4 includes the following sub-steps: S41: Using the calculated downlink Green's function and uplink Green's function first arrival wave to obtain the cross-correlation function, the calculation formula is: Where, C(t) is the cross-correlation function; is the first arrival wave of the upward Green's function; t and t' are time; N is the number of channels; S42: Use the upgoing Green's function first arrival wave to calculate the amplitude compensation factor. The calculation formula is: S42: Perform amplitude compensation on the cross-correlation function, and take the value at time t=0 as the imaging result. The imaging result expression is: Where I is the Marchenko imaging result of the underground imaging point.

Images