One-way wave pre-stack depth migration imaging method

Abstract

The invention belongs to the technical field of seismic exploration, and particularly relates to a one-way wave pre-stack depth migration imaging method. The one-way wave pre-stack depth migration imaging method provided by the invention comprises the following steps: uniformly dividing a certain section of an underground half space into N areas with included angles not greater than 60 degrees, taking an included angle bisector in each area as a paraxial direction, and performing approximate one-way wave recursion along the paraxial direction by adopting a split step Fourier (SSF) operator to obtain a forward transmission wave field; performing back propagation wave field recursion on all depths to obtain a back propagation wave field; and applying cross-correlation imaging conditions to the forward wave fields and the backward wave fields of all positions to obtain a single-shot imaging result of a certain seismic source, and superposing the single-shot imaging results of all seismic sources to obtain a whole-region seismic imaging result. By adopting the SSF operator to carry out stepping recursion along the paraxial direction and adopting a partition mode of multiple paraxial approximation, flexible approximation of a wave field along any angle is realized, a more accurate large-angle seismic wave field is obtained, and the imaging precision of a steep dip angle structure is improved.

Description

Technical Field

[0001] This invention relates to the field of seismic exploration technology, specifically to a one-way wave pre-stack depth migration imaging method. Background Technology

[0002] In recent years, as oil exploration targets have shifted from simple structures such as anticlines and basins to complex structures with steep dips, such as folds, faults, and fracture zones, the problem of imaging steep dips in complex media has gradually become one of the core issues of concern in the academic and industrial communities of seismic data processing. Especially with the rapid development of shale oil and gas resources in recent years, the demand for high-precision imaging of complex steep dip structures in the oil and gas sector is increasing to avoid drilling through layers due to inaccurate seismic dip angles in shale gas horizontal wells. Among these methods, pre-stack depth migration imaging based on one-way waves offers the flexibility of step-by-step recursion along spatial directions and avoids the inherent problems of strong low-frequency noise and high computational cost in shallow layers inherent in conventional wave equation migration imaging methods. It is a mature, reliable, cost-effective seismic imaging method particularly suitable for industrial production applications.

[0003] However, the one-way wave migration method originates from the parabolic approximation theory, which can only solve the first 2-3 terms of the full wave equation using a limited one-way wave equation expansion. This results in the recursive operator of the one-way wave equation having high approximation accuracy only within a certain angular range of the equation's approximate direction. As the propagation angle increases, the wave field recursion error also gradually increases (specifically manifested as amplitude and phase errors in the wave field), leading to insufficient imaging accuracy of one-way wave imaging technology in steeply dipping structures. This reduces the quality of seismic data, increases the risk to subsequent geological interpretation and well location design, and severely limits its industrial application in complex steeply dipping structures.

[0004] To address the insufficient imaging accuracy of pre-stack depth migration imaging with single-pass waves at large angles and steep dips, the focus is on how to accurately solve for the large-angle information of the single-pass wavefield with a finite order, and then use this accurate large-angle wavefield information to improve the imaging accuracy of steep dip structures. To this end, academia and industry have conducted extensive research and technological breakthroughs in areas such as increasing the expansion order of the single-pass wavefield, optimizing coefficients, improving solution methods, and using multi-paraxial approximations.

[0005] Regarding increasing the order of equation expansion: In December 1994, *Geophysics* published Ristow and Ruehl's "Fourier finite-difference migration," proposing the Fourier finite difference (FFD) method, which offers greater angular adaptability. Compared to the conventional split-step Fourier method, this method adds higher-order difference compensation terms in the spatial domain, resulting in a wider angular adaptability range. However, this method requires transforming the equations to multiple domains for separate solutions (from the frequency-wavenumber domain to the frequency-space domain and then to the time-space domain), and the higher-order difference compensation terms also face a series of problems, such as solution stability.

[0006] Regarding coefficient optimization: The IEEE Transactions on Propagation, Vol. 8, 2020, published Tang et al.'s "Optimized Pseudo-Pade Fourier Migrator in Terms of Propagation Angles." The core idea of ​​this method is to simultaneously solve for the propagation angle of the wavefield during wavefield extrapolation and to use different difference coefficients for wavefield components at different propagation angles. This avoids the problem that a single difference coefficient cannot adapt to solving wavefields at all angles, thus improving the ability to image steep dip angles for one-way wave equations. However, the approximate accuracy of large-angle wavefields using this method is still limited, generally not exceeding 60°.

[0007] Regarding improvements in equation-solving methods: The April 1983 issue of the *Chinese Journal of Geophysics* published Ma Zaitian's "Split Algorithm for Migration of Higher-Order Equations." This article proposed a finite-difference higher-order split algorithm for the one-way wave equation, decomposing the higher-order equation into multiple lower-order equations, and then equivalently solving the higher-order equation through these lower-order equations. This solved the problem of solving higher-order equations and enabled accurate solutions for large-angle one-way wave information. However, this method requires multiple splits of the equations and still cannot achieve accurate solutions for the entire angular wave field.

[0008] Regarding increasing the number of paraxial axes: The journal *Geophysics*, 2022, published Meng XY et al.'s work, "Calculating multiple scattering from a time-varying undulating sea surface using a full-wavefield multistage algorithm." The core idea of ​​this research is to use Ristow and Ruehl's FFD algorithm for solving large-angle wavefields as the core, approximating and recursively applying the one-way wavefield along the horizontal and depth directions respectively. Then, by weighted superposition of wavefields in different directions, the weighted one-way wavefield simultaneously possesses high approximation accuracy in all directions. However, this method only expands the paraxial axes along the horizontal x and depth z directions, and cannot flexibly expand the paraxial axes along multiple axes (e.g., 10°, 30°, etc.) in specific spatial directions. Therefore, to compensate for the shortcomings caused by too few approximation axes and unreasonable approximation axis directions, this method can only achieve wavefield solutions at all angles based on the FFD operator with large-angle approximation, and suffers from the inherent problem of differential instability in higher-order correction terms of the FFD operator.

[0009] Therefore, the above-mentioned method of increasing the number of paraxial elements still suffers from insufficient imaging accuracy at large-angle and steep-tilt structures. Summary of the Invention

[0010] This invention provides a single-pass wave pre-stack depth migration imaging method, which solves the problem of insufficient imaging accuracy at large-angle and steep-tilt structures in the prior art.

[0011] To solve the above-mentioned technical problems, the technical solution of the single-pass pre-stack depth migration imaging method of the present invention is as follows:

[0012] A single-pass pre-stack depth migration imaging method includes the following steps:

[0013] (1) Load the migration imaging parameters, velocity model and seismic data, take a certain source as the origin, and divide a certain section of the underground half space into N regions with an angle of no more than 60°. Take the angle bisector of each region as the paraxial axis.

[0014] (2) Using the source function as the boundary condition, the one-way wave recursion is performed along a certain paraxial direction using the split-step Fourier operator to obtain the one-way wave field along a certain paraxial direction.

[0015] Repeat the above one-way wave recursion method to obtain the one-way wave field along all paraxial directions one by one. The one-way wave field is cut off to obtain the one-way wave field of the region where all paraxial directions are located. The one-way wave fields of the regions where all paraxial directions are located are added together to obtain the propagating wave field.

[0016] (3) Using seismic data as boundary conditions, the back propagation wavefield at all depths is recursively derived using an approximate one-way operator along the depth direction to obtain the back propagation wavefield.

[0017] (4) Apply cross-correlation imaging conditions to the forward propagation wave field of step (2) and the reverse propagation wave field of step (3) to obtain the single-shot imaging result of a certain source.

[0018] (5) Repeat steps (1)-(4) to obtain single-shot imaging results of all seismic sources in the entire area;

[0019] The single-shot imaging results of all earthquake sources are superimposed to obtain the seismic imaging results for the entire region.

[0020] This invention improves upon existing technologies by providing a single-pass pre-stack depth migration imaging method. It employs a SSF operator with an approximate angle of ±30°, it iteratively iteratively along the paraxial direction, using an unconditionally stable split-step Fourier (SSF) operator as its core. This eliminates the need for solving higher-order finite difference terms in the existing Fourier finite difference (FFD) operator, thus theoretically avoiding the instability problem associated with solving higher-order difference terms in the FFD operator. Furthermore, it extends the existing paraxial approach, which approximates only the horizontal and depth directions, to multiple paraxial directions within a half-space. The approximate partitioning mode enables flexible approximation of the wavefield along any angle, thereby achieving high-precision solution of the wavefield across all angles and obtaining a more accurate large-angle seismic wavefield. The partitioning method is more reasonable and can avoid the waste of accurate small-angle wavefield near the paraxial region, enabling better solution of large-angle one-way wavefield information. Introducing the large-angle one-way wavefield information of multi-paraxial approximation into seismic imaging improves the imaging accuracy of steep-dip structures. A one-way wave migration imaging method for steep-dip structures is proposed, which makes up for the lack of imaging accuracy of steep-dip structures in traditional one-way wave imaging.

[0021] This invention realizes a more flexible and reasonable multi-paraxial approximation single-pass wavefield migration imaging method, which has important theoretical and practical significance for high-precision numerical simulation of large-angle single-pass wavefields and complex steep-dip structural seismic imaging.

[0022] To further improve imaging accuracy and reduce computational load, preferably, the N regions are 3, 4, 5, or 6. When N is 3, the included angle is 60°, and the paraxial region is the bisector of the included angle. In this case, the angle range of the paraxial region is within ±30°. The split-step Fourier (SSF) operator used has high accuracy within ±30°. The other regions are similar and will not be elaborated here. By limiting the N regions to 3, 4, 5, or 6, the accuracy of wavefield recursion is improved while the workload is greatly reduced.

[0023] To further improve imaging accuracy, preferably, the migration imaging parameters include the spatial sampling interval and number of sampling points of the model in the horizontal and depth directions, the time sampling interval and number of sampling points for wavefield recursion, the frequency sampling interval and number of frequencies, and the dominant frequency of the source function; the velocity model includes the seismic wavefield propagation velocity at all spatial sampling points; and the seismic data includes the source location and depth, and the recorded wavefield values ​​of all sampling points on the Earth's surface at different times.

[0024] To further improve the flexibility of paraxial rotation and the flexibility of approximating the wave field along any angle, thereby improving the high-precision solution across all angles, preferably, the specific recursive method for the one-way wave field along a certain paraxial direction is as follows:

[0025] S1: Define the z-axis angle A along the depth direction as 0°, with A > 0° on the right and A < 0° on the left. i Let the angle be a certain paraxial angle, centered on the source s(sx, sz) and passing through the mapping function R1(A i Rotate the velocity field v(x, z) clockwise by A i The velocity field v(x1, z1) in the mapped coordinates is obtained.

[0026] S2: Using the source function as the boundary condition, an approximate one-way wave recursion is performed along a certain paraxial direction using the split-step Fourier operator:

[0027]

[0028] in, For the positive wave field in the frequency domain under mapped coordinates, E z1 Let ω be the one-way recursive operator for the wavefield along the depth z1 direction in the mapped coordinates, where ω is the frequency and c is the background velocity.

[0029] S3: After completing the wave field recursion along the z1 direction at all positions, obtain the one-way wave field along a certain paraxial direction through counterclockwise rotation and inverse Fourier transform.

[0030] Understandably, one can first use the inverse Fourier transform to obtain the wave field in the time-space domain, and then perform a counterclockwise rotation of A. i To obtain the forward propagation wave field, you can first rotate A counterclockwise. i The forward-scattered wave field in the frequency domain is obtained, and then the forward propagation wave field in the time-space domain is obtained by inverse Fourier transform.

[0031] To further improve the flexibility of paraxial rotation and the flexibility of wavefield approximation along any angle, thereby improving the high-precision solution across all angles and making the large-angle wavefield more accurate, preferably, step S3 completes the wavefield at all positions z1. After recursion, the inverse Fourier transform is used to transform the forward scattered wave field in the frequency domain. Transform to the time-space domain Then, taking the epicenter s(sx, sz) as the center, Through the mapping function R -1 (A i v) Rotate A counterclockwise i Degrees, to obtain the paraxial direction A at different times t. i One-way wave field of angle

[0032] To further reduce wavefield error, preferably, in step (2), the one-way wavefield cut-off involves cutting off the wavefield outside the region containing all paraxial directions, resulting in the one-way wavefield of the region containing all paraxial directions being called paraxial A. i Nearby (A) i The wave field within the range of ±180 / 2N° is specifically as follows:

[0033]

[0034] Where α is the angle between the imaging point (x, z) and the seismic source s (sx, sz).

[0035] To further improve the accuracy of the wavefield along all directions, preferably, the wavefield summation involves summing the one-way wavefields of all paraxial directions, specifically:

[0036]

[0037] Among them, A1, A2…A N Para-axis pointing to A i Forward propagation wave field at an angle.

[0038] To further improve paraxial rotation flexibility, preferably, the mapping function R1(A) i v) is a clockwise rotation matrix, specifically:

[0039] x1=(x-sx)cos(A i )+(z-zs)sin(A i )

[0040] z1=-(x-sx)sin(A i )+(z-zs)cos(A i )

[0041] Where (x, z) represents the true coordinates and (x1, z1) represents the mapped coordinates.

[0042] It is understandable that the mapping function R -1 (A iv) is a counterclockwise rotation matrix, specifically:

[0043] x=(x1-sx)cos(-A i )+(z1-zs)sin(-A i )

[0044] z = -(x1-sx)sin(-A) i )+(z1-zs)cos(-A i )

[0045] Where (x, z) represents the true coordinates and (x1, z1) represents the mapped coordinates.

[0046] To further improve the accuracy of the reverse propagation wavefield recursion, preferably, the recursion method for the reverse propagation wavefield recursion is as follows:

[0047] Step 1: Using seismic data as boundary conditions, perform backpropagation wavefield recursion using an approximate one-way operator along the depth direction:

[0048]

[0049] Among them, U b (ω,x,z) represents the inverse wave field in the frequency domain under the mapped coordinates, E z1 Let v(x, z) be the one-way recursive operator for the wave field along the depth z direction in the mapped coordinates, where ω is the frequency, c is the background velocity, and v(x, z) is the velocity field at different coordinate positions.

[0050] Step 2: After completing the wave field recursion in the z-direction at all positions, obtain the backpropagating wave field U by counterclockwise rotation and inverse Fourier transform. b (x, z, t).

[0051] To further improve the accuracy of single-shot imaging results, preferably, the single-shot imaging results are calculated as follows:

[0052]

[0053] Among them: U f (x, z, t) represents the forward propagation wave field at spatial location (x, z) at time t, recursively derived from the earthquake source; U b (x, z, t) represents the backpropagating wavefield at spatial location (x, z) at time t, derived from seismic data; I(x, z) represents the imaging result at spatial location (x, z), and nt represents the number of time sampling points for wavefield recursion. Attached Figure Description

[0054] Figure 1 A comparison diagram of the multi-directional paraxial approximation provided by the present invention and the prior art;

[0055] Figure 2 A flowchart of Embodiment 1 of the single-pass pre-stack depth migration imaging method provided by the present invention;

[0056] Figure 3 This is a schematic diagram of the velocity field loaded in Embodiment 1 of the present invention;

[0057] Figure 4 The first shot reflection seismic data loaded in Embodiment 1 of the present invention;

[0058] Figure 5 This is an approximate one-way wavefield diagram along paraxial 1 (A = 67.5°) of Embodiment 1 of the present invention;

[0059] Figure 6 This is a wavefield diagram of the area near paraxial A1 retained after large-angle wavefield cutoff in Embodiment 1 of the present invention;

[0060] Figure 7 This is a propagating wavefield diagram approximated along all paraxial directions in Embodiment 1 of the present invention;

[0061] Figure 8 This is a full-area seismic imaging result map of Embodiment 1 of the present invention;

[0062] Figure 9 This is the forward propagation wavefield diagram obtained by approximating along all paraxial directions using the traditional single paraxial wide-angle approximation approach.

[0063] Figure 10 This is an image of the multi-shot data obtained using the traditional single paraxial wide-angle approximation method. Detailed Implementation

[0064] The technical concept of the single-pass wave pre-stack depth migration imaging method of the present invention is as follows:

[0065] A single-pass pre-stack depth migration imaging method includes the following steps:

[0066] (1) Load the migration imaging parameters, velocity model and seismic data, take a certain source as the origin, and divide a certain section of the underground half space into N regions with an angle of no more than 60°. Take the angle bisector of each region as the paraxial axis.

[0067] (2) Using the source function as the boundary condition, the one-way wave recursion is performed along a certain paraxial direction using the split-step Fourier operator to obtain the one-way wave field along a certain paraxial direction.

[0068] Repeat the above one-way wave recursion method to obtain the one-way wave field along all paraxial directions one by one. The one-way wave field is cut off to obtain the one-way wave field of the region where all paraxial directions are located. The one-way wave fields of the regions where all paraxial directions are located are added together to obtain the propagating wave field.

[0069] (3) Using seismic data as boundary conditions, the back propagation wavefield at all depths is recursively derived using an approximate one-way operator along the depth direction to obtain the back propagation wavefield.

[0070] (4) Apply cross-correlation imaging conditions to the forward propagation wave field of step (2) and the reverse propagation wave field of step (3) to obtain the single-shot imaging result of a certain source.

[0071] (5) Repeat steps (1)-(4) to obtain single-shot imaging results of all seismic sources in the entire area;

[0072] The single-shot imaging results of all earthquake sources are superimposed to obtain the seismic imaging results for the entire region.

[0073] Unlike the traditional approximation along a single paraxial direction ( Figure 1 (a) shown) and the approximations of multi-paraxial coordinates along the horizontal and depth directions in recent years ( Figure 1 (b) As shown, this invention is based on a more flexible partitioning method ( Figure 1 (c) presents a method and process for calculating large-angle wavefields using multiple paraxial approaches, thereby improving the imaging accuracy of steep-tilt structures. The core idea is as follows: First, the FFD single-pass approximation operator, which can approximate angles up to ±45° and is widely used in previous research, is replaced with the SSF operator, which can approximate angles up to ±30°. Then, the paraxial approximation in previous research, which only approximates along two directions (horizontal and depth), is extended to a partitioned pattern of multiple paraxial approximations within half-space.

[0074] The single-pass pre-stack depth migration imaging method provided by this invention employs a paraxial recursive stepwise operation of the SSF operator with an approximate angle of ±30°. Using an unconditionally stable split-step Fourier (SSF) operator as its core, it eliminates the need for solving higher-order finite difference terms in the existing Fourier finite difference (FFD) operator, thus theoretically avoiding the instability problem associated with solving higher-order difference terms in the FFD operator. Furthermore, it extends the existing paraxial approximation method, which only approximates along two directions (horizontal and depth), to multiple paraxial approximation partitions within a half-space. This method enables flexible approximation of the wavefield along any angle, thereby achieving high-precision solution of the wavefield across all angles. This yields a more accurate large-angle seismic wavefield, with a more reasonable partitioning method. It avoids wasting the accurate small-angle wavefield near the paraxial region and can better solve for large-angle one-way wavefield information. By introducing the large-angle one-way wavefield information of multi-paraxial approximation into seismic imaging, the imaging accuracy of steep-dip structures is improved. A one-way wave migration imaging method for steep-dip structures is proposed, which makes up for the lack of imaging accuracy of steep-dip structures in traditional one-way wave imaging.

[0075] Based on previous research on multi-paraxial approximation, we have developed a large-angle single-pass wavefield calculation method that allows for flexible selection of paraxial direction and flexible partitioning, and provides accurate approximation in all directions. This method enables single-pass pre-stack depth migration imaging for steep dip angles, which has extremely important theoretical and practical significance.

[0076] The embodiments of the present invention will be further described below with reference to specific examples.

[0077] I. Specific Embodiments of the Single-Wave Pre-Stack Depth Migration Imaging Method of the Present Invention

[0078] Example 1

[0079] The single-pass pre-stack depth migration imaging method of this embodiment is illustrated in the flowchart below. Figure 2 As shown, it includes the following steps:

[0080] (1) Load migration imaging parameters, velocity model, and seismic data. The migration imaging parameters include: the spatial sampling intervals dx and dz and the number of sampling points nx and nz in the horizontal and depth directions of the model; the time sampling interval dt and the number of sampling points nt for wavefield recursion; the frequency sampling interval dω and the number of frequencies nω; and the dominant frequency fm of the source function (generally the same as the dominant frequency of the seismic data). The velocity model includes: the seismic wavefield propagation velocity v(x, z) at all spatial sampling points; and the seismic data includes: the source location and depth s(sx, sz), and the wavefield values ​​recorded at all sampling points on the surface at different times.

[0081] The above parameters constitute the known conditions of this embodiment. Specifically, in this embodiment, the offset imaging parameters are: spatial sampling interval dx = dz = 10m; number of spatial sampling points nx = 401, nz = 301; time sampling interval dt = 1ms, number of sampling points nt = 4000 for wavefield recursion; frequency sampling interval dω = 1Hz and number of frequencies nω = 155; velocity model: the applied velocity field v(x, z) is as follows... Figure 3 As shown, the spatial sampling interval dx = dz = 10m, the number of spatial sampling points nx = 401, and nz = 301; to verify the steep-dipping imaging capability of the imaging algorithm, this embodiment selected a spherical anomaly velocity model, with a seismic wave propagation velocity of 3000m / s and a background velocity of 1500m / s; seismic data: in terms of the observation system, there are a total of 10 shots, with initial shot point positions of (50m, 10m) and a shot spacing of 350m; Figure 4 For single-shot seismic data with shot points at (50 m, 10 m), the time sampling interval dt = 1 ms, the number of sampling points nt = 4000, the total data duration is 4 seconds, and the total length of the data observation array is 4 kilometers.

[0082] (2) Taking a certain earthquake source as the origin, a certain section of the underground half space is uniformly divided into N regions (the included angle of each region is (180 / N)°); it is stipulated that the angle A along the depth direction is 0°, the angle A on the right side of the z-axis is 0°, and the angle A on the left side is 0°;

[0083] The principle for dividing N regions is to ensure that the included angle ((180 / N)°) of each region is no greater than 60°. This is because: the present invention uses an SSF operator with an approximate angle of ±30°, which can achieve accurate solution of the one-way wave field in the range of ±30° near the paraxial axis.

[0084] It is understood that the present invention can be realized by dividing a certain section of the underground half-space into 3, 5, 6 or more regions. In order to further reduce the calculation difficulty, it is preferable that the certain section of the underground half-space is divided into 3, 4, 5 or 6 regions. More preferably, the certain section of the underground half-space is divided into 4 regions, which can improve the accuracy of subsequent single-way wave recursion without increasing the workload of calculation.

[0085] In this embodiment, taking a certain seismic source as the origin, a certain cross-section of the underground half-space is uniformly divided into 4 regions. The included angle ranges of the 4 regions are 90°-45°, 45°-0°, 0°--45°, and -45°--90°, respectively. Paraxial lines 1, 2, 3, and 4 are the angle bisectors of each region. Specifically, the directions of paraxial lines 1, 2, 3, and 4 are: A1 = 67.5°, A2 = 22.5°, A3 = -22.5°, and A4 = -67.5°. The included angle of each region is (180 / N)° = 45°. The approximate schematic diagram of the multi-directional paraxial lines of the 4 regions is shown below. Figure 1 As shown in (c).

[0086] (3) Centered on paraxial 1 with angle A = A1 = 67.5°, the split-step Fourier operator is used to perform approximate one-way wave recursion along the A1 direction. The specific recursion method is as follows:

[0087] S1: Coordinate transformation is performed using the mapping function R1(A1, v), rotating the velocity field by A1 degrees with the source s(sx, sz) as the center. The specific expression is as follows:

[0088] v(x1,z1)=R1(A1,v(x,z))

[0089] The physical meaning of the mapping function R1(A1, v) is a clockwise rotation of the matrix v, and its specific expression is:

[0090] x1=(x-sx)cos(A1)+(z-zs)sin(A1)

[0091] z1=-(x-sx)sin(A1)+(z-zs)cos(A1)

[0092] Where (x, z) represents the true coordinates and (x1, z1) represents the mapped coordinates.

[0093] Following this method, the velocity field v(x, z) loaded in step (1) can be rotated from the real coordinates to the corresponding paraxial direction to obtain the velocity field v(x1, z1) in the mapped coordinates; starting from the first shot, the velocity field obtained by rotating around the source position s(sx, sz) = s(50 m, 10 m);

[0094] S2: Using the source function as the boundary condition, an approximate one-way wave recursion is performed along the A1 direction using the split-step Fourier operator. The specific expression is as follows:

[0095]

[0096] in, For the positive wave field in the frequency domain under mapped coordinates, E z1 Let v(x1, z1) be the one-way recursive operator for the wave field along the depth z1 direction in the mapped coordinates, where ω is the frequency, c is the background velocity, and v(x1, z1) is the velocity field obtained in the mapped coordinates as described above.

[0097] S3: Complete the forward propagation wave field in the z1 direction at all positions. After recursion, the inverse Fourier transform is used to transform the forward scattered wave field in the frequency domain. Transform to the time-space domain Will Through the mapping function R -1 The coordinate transformation of (A1, v) is performed by rotating it counterclockwise by A1 degrees around the source s(sx, sz). Its expression is similar to the mapping function R1(A1, v) mentioned above, specifically:

[0098]

[0099] Among them, R -1 for:

[0100] x=(x1-sx)cos(-A1)+(z1-zs)sin(-A1)

[0101] z=-(x1-sx)sin(-A1)+(z1-zs)cos(-A1)

[0102] Where (x, z) represents the true coordinates and (x1, z1) represents the mapped coordinates;

[0103] Thus, the approximate one-way wave recursion along the paraxial A1 direction of the propagating wavefield has been completed, yielding... That is, the paraxial propagating wave field pointing to A1 degrees at different times t in space (x, z), which can be considered accurate within a certain range (±30°) near the paraxial A1 degree.

[0104] The approximate one-way wavefield along paraxial 1 (A = 67.5°) with the source s (50 m, 10 m) as the origin is as follows: Figure 5 As shown, from Figure 5 It can be seen that the approximate accuracy of the wave field is highest along the paraxial direction; while in the direction away from the paraxial direction, the error of the wave field gradually increases with the increase of the propagation angle; specifically, it manifests as: wavefront shape error (no longer circular) and attenuation of wave field energy.

[0105] (4) Determine the angles of all imaging points (x, z) in space relative to the source position, and perform wavefield cutoff based on these angles, retaining only the wavefield within the range of (A1±180 / 2N)° near paraxial A1. The specific expression is as follows:

[0106]

[0107]

[0108] Where α is the angle between the imaging point (x, z) and the earthquake source s (sx, sz);

[0109] right Figure 5 Wavefield clipping is performed on the wavefield, preserving the wavefield within the range of (A1±180 / 2N)° (±22.5°) around A1=67.5°. Figure 6 As shown, the wave field has the highest approximate accuracy within this angular range, while outside this angular range, the wave field exhibits a certain degree of error.

[0110] (5) Repeat steps (3)-(4) until all paraxial axes A = A1, A2...A are completed. N The approximate recursion of the wavefield in the directions (N=4; A1=67.5°, A2=22.5°, A3=-22.5° and A4=-67.5°) yields the one-way wavefield along all paraxial directions. The one-way wavefields in the regions containing all paraxial directions (i.e., all areas near the paraxial axis (Ai±180 / 2N)° (±22.5°)) are then added together to form a wavefield that accurately approximates the propagation wavefield within the range of -90° to +90° in space. The specific expression of this wavefield is as follows:

[0111]

[0112] Among them, A1, A2…A N Para-axis pointing to A i Forward propagation wave field at an angle.

[0113] In this embodiment:

[0114]

[0115] The forward propagation wavefield (one-way wavefield) approximated along all paraxial directions with the source s (50 m, 10 m) as the origin is as follows: Figure 7 As shown, from Figure 7 It can be seen that when approximations and weighted combinations are made along all paraxial axes, the resulting wave field is accurate in all directions, and the wave phenomenon of the spherical scatterer is clearer and more reasonable.

[0116] (6) The backpropagation wavefield at all depths is recursively derived using an approximate one-way operator along the depth direction. The specific recursion method is as follows:

[0117] Step 1: Using seismic data as boundary conditions, perform backpropagation wavefield recursion using an approximate one-way operator along the depth direction:

[0118]

[0119] The physical meaning of this formula is the same as that of the recursion and single-pass recursion operators in step (3). Since it is a reverse recursion, the sign changes from +i to -i, and the boundary conditions of the recursion, which used the earthquake source as the boundary condition, are replaced by the earthquake data as the boundary condition; U b (ω,x,z) represents the inverse wave field in the frequency domain under the mapped coordinates, E z1 Let v(x, z) be the one-way recursive operator for the wave field along the depth z direction in the mapped coordinates, where ω is the frequency, c is the background velocity, and v(x, z) is the velocity field at different coordinate positions.

[0120] Step 2: Complete the reverse propagation wave field U in the z-direction at all positions. b After recursion, the inverse Fourier transform is used to transform the frequency domain inverse propagated wave field U. b (x, z, ω) is transformed into the time-space domain antipropagating wave field U. b (x, z, t).

[0121] (7) Apply the cross-correlation imaging condition to the forward propagation wavefield from step (5) and the reverse propagation wavefield from step (6) to obtain the single-shot imaging result of a certain source (source number A = A1), the specific expression of which is:

[0122]

[0123] Among them: U f (x, z, t) represents the forward propagation wave field at spatial location (x, z) at time t, recursively derived from the earthquake source; U b (x, z, t) represents the backpropagating wave field recursively derived from seismic data at spatial location (x, z) at time t; I A＝A1(x, z) represents the imaging result at the lower spatial location (x, z), and nt represents the number of time sampling points for wavefield recursion.

[0124] (8) Based on the number of shots in the seismic data for the entire region, repeat steps (1)-(7) to obtain the single-shot imaging results of all sources in the region. Then, superimpose the single-shot imaging results of all sources A1-AN to obtain the seismic imaging results for the entire region. The superposition process is as follows:

[0125]

[0126] The seismic imaging results for the entire region are as follows: Figure 8 As shown, from Figure 8 It can be seen that the imaging accuracy of the steeply tilted spherical target structure obtained by the single-pass pre-stack depth migration imaging method provided by this invention is significantly improved. This is because this invention adopts a more rationally partitioned multi-directional parabolic approximation method, and the recursive single-pass wavefield is more accurate near the spherical target (e.g., ...). Figure 7 As shown in the figure, this can provide a higher quality and more accurate forward propagation wave field for subsequent cross-correlation imaging steps, so the imaging results at steep tilt angles are more accurate.

[0127] II. Comparative Example

[0128] This comparative example uses a single-path wave pre-stack depth migration imaging method, employing the traditional single paraxial wide-angle approximation approach (such as...). Figure 1 (a) Perform one-way pre-stack depth migration imaging, with the source s (50 m, 10 m) as the origin, and approximate the forward propagation wavefield (one-way wavefield) along all paraxial directions as follows: Figure 9 As shown, from Figure 9 It can be seen that, due to the angular limitation of parabolic unfolding in traditional methods, the wave field error is large at large angles far from the paraxial direction. Specifically, the wave field attenuates severely near the spherical target, and the wavefront shape does not match the theoretical wave phenomenon.

[0129] Imaging results from multi-shot data are as follows Figure 10 As shown, from Figure 10 It can be seen that, due to the lack of large-angle information in the wavefield derived by traditional methods, the wavefield near the spherical target body is severely attenuated and inaccurate (see...). Figure 9 This leads to inaccurate cross-correlation imaging results in subsequent step 10, reducing the imaging quality of steep angles near the spherical target.

[0130] Finally, it should be noted that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A single-pass wave pre-stack depth migration imaging method, characterized in that, Includes the following steps: (1) Load the migration imaging parameters, velocity model and seismic data, take a certain source as the origin, and divide a certain section of the underground half space into N regions with an angle of no more than 60°. Take the angle bisector of each region as the paraxial axis. (2) Using the source function as the boundary condition, the one-way wave recursion is performed along a certain paraxial direction using the split-step Fourier operator to obtain the one-way wave field along a certain paraxial direction. Repeat the above one-way wave recursion method to obtain the one-way wave field along all paraxial directions one by one. The one-way wave field is cut off to obtain the one-way wave field of the region where all paraxial directions are located. The one-way wave fields of the regions where all paraxial directions are located are added together to obtain the propagating wave field. (3) Using seismic data as boundary conditions, the back propagation wavefield at all depths is recursively derived using an approximate one-way operator along the depth direction to obtain the back propagation wavefield. (4) Apply cross-correlation imaging conditions to the forward propagation wave field of step (2) and the reverse propagation wave field of step (3) to obtain the single-shot imaging result of a certain source. (5) Repeat steps (1)-(4) to obtain single-shot imaging results of all seismic sources in the entire area; The single-shot imaging results of all earthquake sources are superimposed to obtain the seismic imaging results for the entire region.

2. The single-pass wave pre-stack depth migration imaging method as described in claim 1, characterized in that, The N regions can be 3, 4, 5, or 6.

3. The single-pass wave pre-stack depth migration imaging method as described in claim 1 or 2, characterized in that, The migration imaging parameters include the spatial sampling interval and number of sampling points in the horizontal and depth directions, the time sampling interval and number of sampling points for wavefield recursion, the frequency sampling interval and number of frequencies, and the dominant frequency of the source function; the velocity model includes the seismic wavefield propagation velocity at all spatial sampling points; the seismic data includes the source location and depth, and the recorded wavefield values ​​at all surface sampling points at different times.

4. The single-pass wave pre-stack depth migration imaging method as described in claim 1 or 2, characterized in that, The specific recursive method for the one-way wavefield along a certain paraxial direction is as follows: S1: Define the z-axis angle A along the depth direction as 0°, with A > 0° on the right and A < 0° on the left. i Let the angle be a certain paraxial angle, centered on the source s(sx, sz) and passing through the mapping function R1(A i Rotate the velocity field v(x, z) clockwise by A i The velocity field v(x1, z1) in the mapped coordinates is obtained. S2: Using the source function as the boundary condition, an approximate one-way wave recursion is performed along a certain paraxial direction using the split-step Fourier operator: in, For the positive wave field in the frequency domain under mapped coordinates, E z1 Let ω be the one-way recursive operator for the wavefield along the depth z1 direction in the mapped coordinates, where ω is the frequency and c is the background velocity. S3: After completing the wave field recursion along the z1 direction at all positions, obtain the one-way wave field along a certain paraxial direction through counterclockwise rotation and inverse Fourier transform.

5. The single-pass pre-stack depth migration imaging method as described in claim 4, characterized in that, In step S3, the wave field at all positions z1 is completed. After recursion, the inverse Fourier transform is used to transform the forward scattered wave field in the frequency domain. Transform to the time-space domain Then, taking the epicenter s(sx, sz) as the center, Through the mapping function R -1 (A i v) Rotate A counterclockwise i Degrees, to obtain the paraxial direction A at different times t. i One-way wave field of angle 6. The single-pass wave pre-stack depth migration imaging method as described in claim 1 or 2, characterized in that, In step (2), the one-way wavefield removal involves removing the wavefield outside the region containing all paraxial directions, resulting in the one-way wavefield of the region containing all paraxial directions, which is called paraxial A. i Nearby (A) i The wave field within the range of ±180 / 2N° is specifically as follows: Where α is the angle between the imaging point (x, z) and the seismic source s (sx, sz).

7. The single-pass pre-stack depth migration imaging method as described in claim 6, characterized in that, The wavefield summation involves adding the one-way wavefields of all regions containing the paraxial directions, specifically: Among them, A1, A2…A N Para-axis pointing to A i Forward propagation wave field at an angle.

8. The single-pass wave pre-stack depth migration imaging method as described in claim 4, characterized in that, The mapping function R1(A) i v) is a clockwise rotation matrix, specifically: x1=(x-sx)cos(A i )+(z-zs)sin(A i ) z1=-(x-sx)sin(A i )+(z-zs)cos(A i ) Where (x, z) represents the true coordinates and (x1, z1) represents the mapped coordinates.

9. The single-pass wave pre-stack depth migration imaging method as described in claim 1 or 2, characterized in that, The recursive method for the reverse propagation wavefield is as follows: Step 1: Using seismic data as boundary conditions, perform backpropagation wavefield recursion using an approximate one-way operator along the depth direction: Among them, U b (ω,x,z) represents the inverse wave field in the frequency domain under the mapped coordinates, E z1 Let v(x, z) be the one-way recursive operator for the wave field along the depth z direction in the mapped coordinates, where ω is the frequency, c is the background velocity, and v(x, z) is the velocity field at different coordinate positions. Step 2: After completing the wave field recursion in the z-direction at all positions, obtain the backpropagating wave field U by counterclockwise rotation and inverse Fourier transform. b (x, z, t).

10. The single-pass wave pre-stack depth migration imaging method as described in claim 1 or 2, characterized in that, The single-shot imaging results are calculated as follows: Among them: U f (x, z, t) represents the forward propagation wave field at spatial location (x, z) at time t, recursively derived from the earthquake source; U b (x, z, t) represents the backpropagating wavefield at spatial location (x, z) at time t, derived from seismic data; I(x, z) represents the imaging result at spatial location (x, z), and nt represents the number of time sampling points for wavefield recursion.

Images