A fast Marchenko imaging method for undulating terrain conditions

Abstract

The present invention provides a fast Marchenko imaging method for undulating surface conditions, comprising the following steps: S1: inputting seismic shot records of a fixed surface and an undulating offset surface in the depth domain, and correcting the seismic shot records to the offset surface by a static correction method. S2: inputting a seismic velocity model in the depth domain, and constructing an initial downward focusing function from the imaging point of the underground target area to the offset surface of the undulating surface by a ray tracing algorithm. S3: using the shot records and the initial downward focusing function on the undulating offset surface as input, and calculating the upward and downward Green's functions from the underground imaging point to the ground by an iterative method. S4: applying the energy normalized cross-correlation imaging condition to the reconstructed upward and downward Green's functions to obtain the Marchenko imaging result. The present invention solves the problem that conventional Marchenko imaging has low computational efficiency and cannot handle undulating surface conditions, and ensures the quality of Marchenko imaging under undulating surface conditions while improving computational efficiency.

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, 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 interlayer multiple reflection waves. Consequently, they mistakenly image interlayer multiple waves as single reflection waves, resulting in migration artifacts related to interlayer multiple waves in the imaging results, thus affecting imaging quality. In recent years, the Marchenko imaging method has been developed based on inverse scattering theory. Given a known macroscopic velocity model, the Marchenko imaging method reconstructs uplink and downlink Green's functions from seismic data by solving the Marchenko equation. The reconstructed Green's functions are then used for structural imaging. The Marchenko imaging method can effectively suppress migration artifacts caused by interlayer multiple waves in seismic data within the imaging domain.

[0003] With the increasing exploration and development of oil and gas resources in eastern my country and their increasing scarcity, simple structural exploration can no longer meet industrial needs. Seismic exploration has gradually shifted from the east to the more challenging west, where surface conditions have changed from nearly level to complex, undulating terrain, such as mountainous piedmont zones, hills, loess plateaus, and gravelly Gobi deserts. However, the dramatic surface fluctuations, complex near-surface structures, and large lateral velocity variations in these areas pose significant challenges to seismic data processing.

[0004] Current research and application of the Marchenko imaging method are primarily based on the assumption of a horizontal surface. In practical data applications, particularly for seismic data from rugged western onshore surfaces, Green's function reconstruction and imaging processing require calibrating the shot and receiver points of the seismic data to the same datum. When surface elevation and near-surface lateral velocity vary significantly, using static correction to align the shot and receiver points of the seismic data to the same datum inevitably destroys the hyperbolic shape of the reflected waves, distorting the seismic wavefield. This results in a decrease in the quality of the Marchenko imaging after elevation static correction, and the imaging accuracy cannot meet the requirements of geological interpretation.

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

[0006] To address the shortcomings of the aforementioned Marchenko imaging method, the present invention proposes a fast Marchenko imaging technique for rugged terrain. This technique, in a step-by-step approach, first uses a fast ray tracing algorithm to calculate the traveltime volume from the rugged surface to the underground target area. It then extracts the initial arrival traveltimes from the underground imaging point to the rugged surface. These traveltimes and seismic wavelets are then used to construct the initial descending focusing function from the imaging point to the rugged surface. Finally, Green's function reconstruction and Marchenko imaging are performed. This technique not only improves the accuracy of Marchenko imaging on rugged terrain but also significantly enhances 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 by using static correction to correct the shot check points of seismic data to the same reference plane

[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: Apply the energy-normalized cross-correlation imaging condition to the reconstructed uplink and downlink Green's functions to obtain Marchenko imaging results.

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

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

[0017] S12: Calculate the static correction value T for each seismic data according to the replacement speed from the offset surface to the fixed surface, the coordinates of the shot point receiver 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 receiver 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 receiver 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 for each seismic data upward translation based on the static correction value statics , and then perform static correction on the seismic data based on 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 seismic data.

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

[0024] S21: Count the shot point coordinates and receiver point coordinates of all trace data in the shot record and determine the coordinate range. Then, determine the coordinates of the ray-tracing source point on the undulating offset surface based on the coordinate range. 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 travel time; v is the speed of the ray.

[0027] S22: The travel time τ(s0,s1) of the underground imaging point to the undulating ground surface extracted from the travel time volume of the common earthquake source point calculated in the previous step. 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 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 wavelet; θ(t) is the unit step function.

[0039] S32: Calculate the uplink focusing function based on the downlink focusing function calculated in the previous step. 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. The calculation formula is:

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

[0044]

[0045] Among them, G + and G - are the downward and upward Green 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: Calculate the cross-correlation function using the calculated first arrival waves of the downgoing Green's function and the upgoing Green's 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 is 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 present invention solves the problem that conventional Marchenko imaging has low computational efficiency and cannot process undulating surface conditions, and ensures the quality of Marchenko imaging under undulating surface conditions while improving computational 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 Schematic diagram of a ray tracing path from a source point on an undulating surface to the underground in an embodiment;

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

[0061] Figure 5 A seismic wavelet synthesized for the 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] Figure 9 This is the Marchenko imaging result of the target area under the embodiment. DETAILED DESCRIPTION

[0066] The specific technical solutions of the present invention are described with reference to the embodiments.

[0067] The present invention will be further described below with reference to 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, undulating offset surface, and velocity model, such as Figure 2 As shown in Figure 2, the seismic shot records are corrected to the offset surface using the 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 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] Step 6: 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 and image the next imaging point. If it is, output the Marchenko imaging section, such as Figure 9 shown.

Claims

1. A fast Marchenko imaging method for undulating terrain conditions, characterized by: 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 fluctuation offset surface through the 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 undulating offset 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; The following sub-steps are included: S31: Calculate the complete downlink focusing function using the constructed initial downlink focusing function. The calculation formula is as follows: in, is the downward focusing function including the coda wave; R is the shot record corrected to the undulating offset surface; 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 the duration of the wavelet; θ(t) is a unit step function; S32: Calculate the uplink focusing function based on the downlink focusing function calculated in the previous step. The calculation formula is: in, is the upward focusing function; the superscript "-" indicates the upward direction; S33: Calculate the uplink and downlink Green functions according to the calculated uplink and downlink focusing functions. The calculation formula is: Among them, G + and G - are the downward and upward Green functions respectively; Ψ a and Ψ b is the time window function, which is defined as follows: P a,b (t)=1-Θ a,b (t) (11) S4: Apply the energy-normalized cross-correlation imaging condition to the reconstructed uplink and downlink Green's functions to obtain Marchenko imaging results; The following sub-steps are included: S41: Calculate the cross-correlation function using the calculated first arrival waves of the downgoing Green's function and the upgoing Green's 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.

2. The fast Marchenko imaging method for undulating terrain 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 offset surfaces. Time-domain seismic shot records include the horizontal coordinates and elevations of shot points and receiver points, while depth-domain offset surfaces include the horizontal coordinates and elevations of each sample point. S12: Calculate the static correction value T for each seismic data according to the replacement speed of the undulating offset surface to the fixed surface, the coordinates of the shot point receiver in the shot record, and the elevation of the undulating offset surface. statics , static correction value 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 receiver point; E(sy, sx) is the elevation of the projected point of the shot point on the undulating offset surface; E(gy, gx) is the elevation of the projected point of the receiver point on the undulating 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 for each seismic data upward translation based on 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 seismic data.

3. The fast Marchenko imaging method for undulating terrain 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 trace data in the shot record and determine the coordinate range. Then, determine the coordinates of the ray-tracing source point on the undulating offset surface based on the coordinate range. 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 from point s0 to point s i The travel time of ; v is the speed of the ray; S22: The travel time τ(s0,s1) of the underground imaging point to the undulating ground surface extracted from the travel time volume of the common earthquake source point calculated in the previous step. 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 travel time τ(s0,s i ) replaces t0 in the seismic wavelet r(t) to construct the initial downward focusing function

Images