Pre-stack reverse-time migration imaging method based on GPU parallel technology and boundary storage strategy

Abstract

The present application belongs to the field of seismic exploration data processing, and relates to a pre-stack reverse-time migration imaging method based on GPU parallel technology and boundary storage strategy, comprising: setting basic model parameters and convolution complete matching layer absorption boundary condition parameters, and determining an absorption boundary area; performing forward propagation of a seismic source wave field and saving the wave field absorption boundary area and the seismic source wave field data at the final time; synchronously advancing reverse reconstruction of the seismic source wave field and reverse propagation of a detector wave field, performing cross-correlation calculation on two sound pressure fields at the same time, and obtaining reverse-time migration imaging values after stacking the values at each time; parallelizing wave field propagation calculation based on GPU; obtaining reverse-time migration imaging values corresponding to each gather and performing stacking, obtaining a stacked reverse-time migration imaging value, and performing deep gain and multi-stage filtering processing to obtain a high-precision migration imaging profile. The method has the beneficial effect of being able to maintain high calculation efficiency while obtaining a high-precision migration imaging profile, thereby providing a basis for fine interpretation of complex underground structures.

Description

Technical Field

[0001] This invention relates to the field of seismic exploration data processing technology, and in particular to a pre-stack reverse time migration imaging method based on GPU parallel technology and boundary storage strategy. Background Technology

[0002] Seismic exploration is one of the core methods for detecting mineral resources and understanding deep structures. Migration imaging, as a crucial component, aims to reorient surface-received seismic wavefield data to their actual subsurface reflection locations to obtain high-precision structural images. Reverse-time migration imaging, based on the two-way wave equation, offers advantages such as no dip limitations and adaptability to strong lateral velocity variations, making it one of the most accurate techniques for imaging complex structures.

[0003] Reverse-time migration (RTM) imaging typically involves reconstructing subsurface structures through steps such as forward propagation of the source wavefield, backward propagation of the receiver wavefield, imaging, and post-processing. Due to the limited model region in numerical calculations, seismic waves reach the model boundaries during wavefield propagation, generating reflected waves that affect the RTM imaging results. Absorbing boundary conditions are commonly used to mitigate this impact. Traditional sponge boundaries and conventional perfectly matched layers have limited absorption capacity for reflections, while convolutional perfectly matched layers (CPMLs) can more effectively absorb boundary reflections through complex frequency stretching functions, improving the accuracy of RTM imaging. Furthermore, RTM imaging is susceptible to low-frequency noise and various artifacts due to the two-way wave equation, severely affecting image quality. Therefore, filtering methods are commonly used for post-processing. Laplace filtering is a commonly used denoising method; while it effectively suppresses low-frequency noise, it inevitably damages the effective signal, leading to fuzziness of steep-dipping structures and distortion of amplitude information, affecting the reliability of subsequent interpretation. Nonlinear methods such as median filtering and morphological filtering, although possessing edge-preserving characteristics, easily obscure fault information. Overall, existing filtering methods face a trade-off between noise suppression and signal fidelity.

[0004] In reverse-time migration (RTM) imaging calculations, the conventional process requires storing source wavefield data for all time periods. For example, a 1000×1000 model with 4000 time sampling points requires approximately 29.80 GB of memory to store the data of a single component of the source wavefield in double-precision floating-point format. High-precision RTM imaging involves even more precise model settings and a greater number of time sampling points, leading to significant memory consumption and a heavy computational burden. To alleviate memory pressure, the boundary storage method is widely adopted, which involves saving the wavefield boundary values ​​during the forward propagation of the source wavefield and simultaneously reconstructing them during the reverse propagation of the wavefield at the receiving point. However, while this method reduces memory requirements, it significantly reduces the efficiency of RTM imaging because it requires repeatedly executing the forward propagation process of the source wavefield.

[0005] In summary, there is an urgent need to develop a high-performance reverse time migration imaging method that can synergistically optimize imaging accuracy and computational efficiency, so as to achieve more accurate and efficient imaging of complex geological targets and provide more reliable technical support for subsurface structural analysis and resource exploration. Summary of the Invention

[0006] Technical problems to be solved

[0007] In view of the above-mentioned shortcomings and deficiencies of the prior art, the present invention provides a pre-stack reverse time migration imaging method based on GPU parallel technology and boundary storage strategy, which solves the technical problems of high computational resource consumption, low efficiency, and low imaging accuracy caused by noise interference and weak deep signal energy in the existing reverse time migration imaging methods.

[0008] Technical solution

[0009] To achieve the above objectives, the main technical solutions adopted by the present invention include:

[0010] This invention provides a pre-stack reverse time migration imaging method based on GPU parallel technology and boundary storage strategy, comprising the following steps:

[0011] Step 1: Based on the data of the study area, set the basic model parameters of reverse time migration imaging and the absorption boundary condition parameters of the convolution perfectly matched layer, and determine the absorption boundary region;

[0012] Step 2: Solve the first-order velocity-sound pressure acoustic wave equation using the high-order staggered grid finite difference method, perform forward propagation calculation of the source wavefield, and apply the convolution fully matched layer CPML absorbing boundary condition in the absorbing boundary region to absorb boundary reflections. At the same time, save the wavefield data in the absorbing boundary region at each time step, as well as the source wavefield data at the final time.

[0013] Step 3: Reconstruct the source wavefield in reverse based on the saved wavefield data of the absorbing boundary region and the source wavefield data at the final moment;

[0014] Step 4: Load the common shot gather according to the detector location, and perform reverse propagation of the detector wavefield using the same calculation method as the forward propagation of the source wavefield;

[0015] Step 5: Simultaneously advance the reverse reconstruction of the source wavefield and the reverse propagation of the detector wavefield, perform cross-correlation calculations on the two sound pressure fields at the same time, and superimpose the calculation results at each time along the time direction to obtain the reverse time migration imaging value of the current common shot point gather;

[0016] Step 6: Parallelize wave field propagation calculations based on GPU architecture;

[0017] Step 7: Perform reverse time migration imaging on the seismic records of multiple common shot point gathers, obtain the reverse time migration imaging values ​​corresponding to each gather, and stack them to obtain the final stacked reverse time migration imaging values, thus completing the pre-stack reverse time migration imaging.

[0018] Step 8: Perform depth gain and multi-level filtering on the superimposed reverse time migration imaging values ​​to obtain the final high-precision migration imaging profile.

[0019] As a further improvement to the method of the present invention, in step 1, the basic model parameters of the study area are set according to the study area data, including at least the geometric model, density model, P-wave velocity model, spatial position parameters of the shot point and the detector, and the time-frequency parameters.

[0020] A perfectly matched convolutional layer is introduced as an absorption boundary condition. The number of CPML layers (npml), absorption coefficient (R), and maximum complex frequency shift parameter (α) are set according to the basic model parameters and time-frequency parameters. max Maximum tensile parameter κ max This allows us to determine the absorbing boundary region, set a high-order staggered mesh with a difference order of 2N, and calculate the difference coefficients a.

[0021] As a further improvement to the method of the present invention, in step 2, the first-order velocity-sound pressure acoustic wave equation is solved using a high-order staggered grid finite difference method, and the source wavefield is propagated forward by applying convolutionally perfectly matched layer absorbing boundary condition parameters within the absorbing boundary region. The wavefield data of the absorbing boundary region under the staggered grid and the source wavefield data at the final moment are then saved. Specifically:

[0022] Wavefield simulation was performed based on the first-order velocity-sound pressure acoustic wave equation. Numerical discretization was performed using a high-order staggered grid finite difference scheme. A seismic source was loaded at the shot point, and the source wavefield velocity and sound pressure fields were calculated using basic model parameters. Furthermore, CPML absorbing boundary conditions were applied to the defined absorbing boundary region, calculated according to the following formula:

[0023] ;

[0024] In the formula, Let represent the sound pressure field at index (i, j) at time t. Let be the horizontal component of the velocity field at the position indexed (i, j+0.5) at time t-0.5. Let the vertical component of the velocity field at the index (i+0.5, j) at time t-0.5 be the velocity field component. This represents the model density at index (i,j+0.5). , Let represent the wave velocity of the model P-wave at index (i,j), Δt represent the time sampling interval, and Δx and Δz represent the spatial sampling intervals. Let be the CPML convolution term applied to the horizontal component v of the velocity field at time t. Let w be the CPML convolution term applied to the vertical component w of the velocity field at time t. Let be the CPML convolution term applied to the horizontal component p of the sound pressure field at time t. κ is the CPML convolution term applied to the vertical component of the sound pressure field p at time t. x κ is the horizontal tensile parameter. z For vertical stretching parameters;

[0025] A boundary storage strategy is adopted to save the wavefield data of the source wavefield, sound pressure field, and velocity field at each time step within the absorbing boundary region, as well as the source wavefield data at the final moment.

[0026] As a further improvement to the method of the present invention, in step 2, the storage formula for the wavefield boundary under the staggered grid finite difference scheme is as follows:

[0027] ;

[0028] In the formula, v bd w bd px represents the part where boundary conditions are applied to the velocity field. bd This represents the part of the sound pressure field with boundary conditions applied in the horizontal direction, pz bd This indicates the part where boundary conditions are applied in the vertical direction of the sound pressure field. nx and nz are the total number of sampling points in each direction, and npml is the number of CPML layers set.

[0029] As a further improvement to the method of the present invention, in step 3, the source wavefield is reconstructed in reverse based on the saved wavefield data of the absorbing boundary region and the source wavefield data at the final moment. Specifically:

[0030] Using the source wavefield at the final moment saved during forward propagation as the initial condition for reverse reconstruction of the source wavefield, the source wavefield is reconstructed based on the first-order velocity-sound pressure acoustic wave equation. The absorption boundary region wavefield data saved during forward propagation of the source wavefield is used to replace the absorption boundary region of the reconstructed source wavefield in order to restore the complete source wavefield.

[0031] The initial condition is the source wavefield at the final moment preserved during the forward propagation of the source wavefield, calculated according to the following formula:

[0032] ;

[0033] In the formula, v fd nt wfd nt p fd nt Represents the source wave field at the final moment, v re nt w re nt p re nt This indicates the final moment of reconstruction of the source wave field.

[0034] As a further improvement to the method of the present invention, in step 4, based on the common shot gather loaded at the detector location, the detector wavefield is propagated in reverse using the same calculation method as the forward propagation of the source wavefield. Specifically:

[0035] Using the seismic records of each shot point gather loaded at the detector location as the source, the reverse propagation of the detector wavefield is calculated using the forward propagation calculation method of the source wavefield and the application of CPML absorbing boundary conditions, specifically according to the following formula:

[0036] ;

[0037] In the formula, (zrec, xrec) represents the detector location, Rec(zrec, xrec, t) represents the seismic record received by the detector at location (zrec, xrec) at time t, ρ(zrec, xrec) represents the model density at location (zrec, xrec), and v t+0.5 (zrec, xrec), w t+0.5 (zrec, xrec) represent the horizontal and vertical components of the velocity field at position (zrec, xrec) at time t+0.5, respectively.

[0038] As a further improvement to the method of the present invention, in step 5, the reverse reconstruction of the source wavefield and the reverse propagation of the detector wavefield are carried out simultaneously. Cross-correlation calculations are performed on the two sound pressure fields at the same time, and the calculation results at each time are superimposed along the time direction to obtain the reverse-time migration imaging value of the current common shot point gather. Specifically:

[0039] Based on the principle of time consistency, that is, taking advantage of the simultaneous existence of the source wave field as the incident wave and the detector wave field as the reflected wave at the reflection point, the reverse reconstruction of the source wave field and the reverse propagation of the detector wave field are synchronously promoted from the maximum time. The cross-correlation calculation of the two sound pressure fields at the same time is performed to obtain the dimensionless instantaneous imaging value until the time returns to zero.

[0040] During the acquisition of imaging values, the instantaneous imaging values ​​obtained at each moment are saved. The instantaneous imaging values ​​at all moments are superimposed along the time direction to form the reverse-time migration imaging values ​​of the current common shot point gather seismic record, calculated according to the following formula:

[0041] ;

[0042] In the formula, i_shot(t) represents the instantaneous imaging value at time t, and I_shot represents the reverse time migration imaging value of the current common shot point gather seismic record.

[0043] As a further improvement to the method of this invention, the wave field propagation calculation is parallelized based on a GPU architecture, specifically:

[0044] The GPU-based CUDA thread architecture parallelizes wave field propagation calculations by dividing a two-dimensional spatial grid into multiple thread blocks. Threads within a thread block collaborate to read data and perform parallel local difference calculations through shared memory.

[0045] As a further improvement to the method of the present invention, in step 7, reverse time migration imaging is performed on the seismic records of multiple common shot point gathers, and the reverse time migration imaging values ​​corresponding to each gather are obtained and superimposed to obtain the final superimposed reverse time migration imaging values, thus completing the pre-stack reverse time migration imaging. Specifically:

[0046] Reverse-time migration imaging was performed on multiple common shot point gather seismic records to obtain the reverse-time migration imaging value corresponding to each common shot point gather seismic record. The spatially consistent point values ​​in the reverse-time migration imaging values ​​corresponding to each common shot point gather seismic record were algebraically summed to obtain the final stacked reverse-time migration imaging value, thus completing the pre-stack reverse-time migration imaging.

[0047] As a further improvement to the method of the present invention, in step 8, the superimposed reverse-time migration imaging values ​​are subjected to depth gain and multi-level filtering to obtain the final high-precision migration imaging profile, specifically:

[0048] First, depth gain processing is applied to the superimposed reverse time migration imaging values. Then, the imaging values ​​after depth gain are sequentially subjected to a first Laplace filter to suppress high-frequency noise, a morphological filter to enhance the structural boundaries, and a second Laplace filter to suppress secondary noise introduced by morphological operations, thereby obtaining a high-precision migration imaging profile with low noise levels and clear structural boundaries.

[0049] Beneficial effects

[0050] The beneficial effects of this invention are:

[0051] The first-order velocity-sonic pressure wave equations were discretized using a high-order staggered grid, and CPML absorbing boundary conditions were employed to effectively suppress boundary reflection interference and improve the numerical accuracy of wavefield simulation. By using boundary storage and wavefield reconstruction strategies, the computational memory requirements were reduced while maintaining wavefield propagation accuracy. A GPU parallel computing architecture was used to achieve efficient execution of wavefield propagation and imaging, significantly improving the computational speed of reverse-time migration imaging. A combination of Laplace filtering and morphological filtering was employed to enhance the continuity and recognizability of structural boundaries, resulting in clearer and more reliable final imaging results. Furthermore, high-precision migration imaging profiles were obtained while maintaining high computational efficiency, providing a foundation for the detailed interpretation of complex subsurface structures and significantly contributing to improving the efficiency of seismic data processing and the reliability of geological interpretation. Attached Figure Description

[0052] Figure 1 A flowchart of a pre-stack reverse time migration imaging method based on GPU parallel technology and boundary storage strategy provided in an embodiment of the present invention;

[0053] Figure 2 This is the P-wave velocity and density model used in the embodiments of the present invention; wherein, Figure 2 Part (a) in the figure represents the P-wave velocity model. Figure 2 Part (b) in the diagram represents the density model;

[0054] Figure 3 This embodiment of the invention presents the forward-propagating source wavefield, the reverse-reconstructed source wavefield, and the residual between the two at t=1s during the reverse-time migration imaging of the common shot point gather for shot 51; wherein, Figure 3 Part (a) in the figure represents the forward propagating source wave field at t=1s; Figure 3 Part (b) in the figure represents the source wave field reconstructed in reverse at t=1s; Figure 3 Part (c) in the figure represents the residual between the forward-propagating source wavefield and the reverse-reconstructed source wavefield;

[0055] Figure 4 This refers to a portion of the common shot point gathers used in the embodiments of the present invention; wherein, Figure 4 Part (a) in the text indicates the 51st shot; Figure 4 Part (b) in the text represents the 100th shot;

[0056] Figure 5 These are the reverse-time migration imaging values ​​of some common shot point gathers in embodiments of the present invention; wherein, Figure 5 Part (a) in the text indicates the 51st shot; Figure 5 Part (b) in the text represents the 100th shot;

[0057] Figure 6These are the superimposed reverse time-shifted imaging values ​​obtained by performing reverse time shifting on the 100 common shot point gathers in this embodiment of the invention.

[0058] Figure 7 In this embodiment of the invention, the superimposed reverse-time offset imaging values ​​after depth gain are used.

[0059] Figure 8 This is an image obtained by superimposing reverse-time offset imaging values ​​and performing multi-level filtering in an embodiment of the present invention; wherein, Figure 8 Part (a) in the image represents the image after the first Laplace filter; Figure 8 Part (b) in the image represents the image after morphological filtering based on the first Laplace filter; Figure 8 Part (c) in the image represents the image after Laplace filtering based on morphological filtering. Detailed Implementation

[0060] To better understand the above technical solutions, exemplary embodiments of the present invention will be described in more detail below with reference to the accompanying drawings. Although exemplary embodiments of the present invention are shown in the drawings, it should be understood that the present invention can be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that the present invention can be understood more clearly and thoroughly, and that the scope of the present invention can be fully conveyed to those skilled in the art.

[0061] In this embodiment, the computing device CPU includes an Intel(R) Core(TM) i7-14650HX and an NVIDIA GeForce RTX 4060 Laptop GPU.

[0062] like Figure 1 As shown, this embodiment of the invention provides a pre-stack reverse time migration imaging method based on GPU parallel technology and boundary storage strategy, including the following steps:

[0063] Step 1: Based on the data of the study area, set the basic model parameters of reverse time migration imaging and the absorption boundary condition parameters of the convolution perfectly matched layer, and determine the absorption boundary region;

[0064] Specifically, based on the data from the study area, basic model parameters for the study area are set, including parameters such as the geometric model, density model, P-wave velocity model, spatial positions of the shot point and detector, and time-frequency parameters. A perfectly matched convolutional layer is introduced as an absorbing boundary condition. Through a complex frequency stretching function, the wave field energy entering from the boundary is efficiently absorbed, effectively suppressing false signals caused by boundary reflection during wave field propagation and reducing the impact of boundary reflection on imaging quality.

[0065] Based on the basic model parameters and time-frequency parameters, set the CPML layer number npml, absorption coefficient R, and maximum complex frequency shift parameter α. max Maximum tensile parameter κ max This allows us to determine the absorbing boundary region, set a high-order staggered mesh with a difference order of 2N, and calculate the difference coefficients a.

[0066] In this embodiment, the Marmousi acoustic model is used to resample the basic model of the study area. The spatial sampling of the basic model is changed from 4m×4m to 10m×10m. The portion from 3000m to 9200m in the horizontal direction and from 0m to 3000m in the depth is selected as the basic model for reverse time migration imaging.

[0067] like Figure 2 As shown in the following description, the basic model ranges from 0m to 6200m laterally and from 0m to 3000m in depth, with 100 shot points at a depth of 300m, uniformly distributed laterally between 776.25m and 5433.75m. To account for geometric sampling limitations, the lateral positions of the shot points are moved to the nearest sampling point. A total of 601 geophones are set, at a depth of 100m, uniformly distributed laterally between 100m and 6100m, with an interval of 10m. The seismic source is a Riker wavelet with a dominant frequency of 20Hz, a time step of 0.001s, and a total of 5000 sampling points. The top surface is a free boundary, and the other three surfaces are absorbing boundaries. The npml is 15 layers, and R is set to 10. -6 α is calculated based on time-frequency parameters. max、 κ max The difference order is second in time and sixth in space, i.e., N is 3. Calculate the difference coefficients a.

[0068] Step 2: Solve the first-order velocity-sound pressure acoustic wave equation using the high-order staggered grid finite difference method, calculate the forward propagation of the source wave field in the entire model space, and apply the convolution fully matched layer CPML absorbing boundary condition in the absorbing boundary region to absorb boundary reflection. At the same time, save the wave field data in the absorbing boundary region at each time step, as well as the source wave field data at the final time.

[0069] Specifically, wave field simulation is performed based on the first-order velocity-sound pressure acoustic wave equation. A high-order staggered grid finite difference scheme is used for numerical discretization to improve computational accuracy and stability. A seismic source is loaded at the shot point, and parameters such as density and P-wave velocity are used to calculate the source wave field velocity field and sound pressure field. Furthermore, CPML absorbing boundary conditions are applied in the set absorbing boundary region, which are calculated according to equation (1):

[0070] (1)

[0071] In the formula, Let represent the sound pressure field at index (i, j) at time t. Let be the horizontal component of the velocity field at the position indexed (i, j+0.5) at time t-0.5. Let the vertical component of the velocity field at the index (i+0.5, j) at time t-0.5 be the velocity field component. This represents the model density at index (i,j+0.5). , Let represent the wave velocity of the model P-wave at index (i,j), Δt represent the time sampling interval, and Δx and Δz represent the spatial sampling intervals. Let be the CPML convolution term applied to the horizontal component v of the velocity field at time t. Let w be the CPML convolution term applied to the vertical component w of the velocity field at time t. Let be the CPML convolution term applied to the horizontal component p of the sound pressure field at time t. κ is the CPML convolution term applied to the vertical component of the sound pressure field p at time t. x κ is the horizontal tensile parameter. z For vertical stretching parameters;

[0072] A boundary storage strategy is adopted to save the wavefield data of the source wavefield, sound pressure field, and velocity field at each time step within the absorbing boundary region, as well as the source wavefield data at the final moment.

[0073] The storage formula for the wavefield boundary under the staggered grid finite difference scheme is shown in equation (2):

[0074] (2)

[0075] In the formula, v bd w bd px represents the part where boundary conditions are applied to the velocity field. bd This represents the part of the sound pressure field with boundary conditions applied in the horizontal direction, pz bd This indicates the part where boundary conditions are applied in the vertical direction of the sound pressure field. nx and nz are the total number of sampling points in each direction, and npml is the number of CPML layers set.

[0076] Using the boundary storage method described above, if the total number of time sampling points is nt, v bd w bd Then nt, v bd w bd The wavefield data sizes are nz×(npml×2)×(nt-1) and (npml×2)×nx×(nt-1), respectively; px bd pz bdThe wavefield data sizes are nz×(npml×2+1)×(nt-1) and (npml×2+1)×nx×(nt-1), respectively. Therefore, the wavefield data size that needs to be saved is (nz+nx)×(npml×4+1)×(nt-1). Compared to the wavefield data size of nz×nx×(nt-1) that needs to be saved without using the boundary storage strategy, the memory required with the boundary storage strategy is (1 / nz+1 / nx)×(npml×4+1). As the total number of sampling points nx and nz in each direction of the wavefield increases, the memory saving ratio of the boundary storage strategy is greater, and its advantage is more significant.

[0077] In this embodiment, based on the geometric model, density, P-wave velocity, etc., the forward propagation of the source wavefield in reverse time migration imaging of a single common shot point gather is performed. The propagation of seismic waves caused by the Riker wavelet source in the model is calculated, and the wavefield data of the absorbing boundary region in the sound pressure field and velocity field at all time points are saved. Through the boundary storage strategy, this embodiment requires approximately 2.09 GB of memory when the data type is double-precision floating point, compared to approximately 6.96 GB of memory required to store the entire sound pressure field in the traditional method, effectively reducing memory requirements.

[0078] Step 3: Reconstruct the source wavefield in reverse based on the saved wavefield data of the absorbing boundary region and the source wavefield data at the final moment;

[0079] Specifically, the source wavefield at the final moment saved during the forward propagation of the source wavefield is used as the initial condition for the reverse reconstruction of the source wavefield. The reverse reconstruction of the source wavefield is carried out based on the first-order velocity sound pressure acoustic wave equation. The absorption boundary region wavefield data saved during the forward propagation of the source wavefield is used to replace the absorption boundary region of the reconstructed source wavefield in order to restore the complete source wavefield.

[0080] Among them, the source wave field at the final moment preserved during the forward propagation of the source wave field is used as the initial condition, and the calculation is performed according to equation (3):

[0081] (3)

[0082] In the formula, v fd nt w fd nt p fd nt Represents the source wave field at the final moment, v re nt w re nt p re nt This indicates the final moment of reconstruction of the source wave field.

[0083] like Figure 3 As shown, the source wavefield is reconstructed inversely based on the wavefield at the final moment of forward propagation and the preserved wavefield data of the absorbing boundary region. The accuracy of the inverse reconstruction method is evaluated by comparing the forward-propagated source wavefield with the reconstructed source wavefield. Figure 3 It can be seen that the reconstructed wavefield is numerically similar to the original wavefield, and the residual is more than five orders of magnitude different from the wavefield value, which can be considered to have little impact on the migration results.

[0084] Step 4: Load the common shot gather according to the detector location, and perform reverse propagation of the detector wavefield using the same calculation method as the forward propagation of the source wavefield;

[0085] like Figure 4 As shown, the seismic records of each shot point gather loaded at the detector location are used as the source. The reverse propagation of the detector wavefield is calculated using the forward propagation calculation method of the source wavefield and the application of CPML absorbing boundary conditions. Specifically, it is calculated according to equation (4):

[0086] (4)

[0087] In the formula, (zrec, xrec) represents the detector location, Rec(zrec, xrec, t) represents the seismic record received by the detector at location (zrec, xrec) at time t, ρ(zrec, xrec) represents the model density at location (zrec, xrec), and v t+0.5 (zrec, xrec), w t+0.5 (zrec, xrec) represent the horizontal and vertical components of the velocity field at position (zrec, xrec) at time t+0.5, respectively.

[0088] In this embodiment, the detector wavefield in reverse time migration imaging of a single common shot point gather is propagated in reverse based on geometric model, density, P-wave velocity, seismic records, etc., and the seismic records received by 601 detectors are used as the source.

[0089] Step 5: Simultaneously advance the reverse reconstruction of the source wavefield and the reverse propagation of the detector wavefield, perform cross-correlation calculations on the two sound pressure fields at the same time, and superimpose the calculation results at each time along the time direction to obtain the reverse time migration imaging value of the current common shot point gather;

[0090] Specifically, based on the principle of time consistency, the source wave field as the incident wave and the detector wave field as the reflected wave exist simultaneously at the reflection point. Starting from the maximum time, the reverse reconstruction of the source wave field and the reverse propagation of the detector wave field are synchronously promoted. The two sound pressure fields at the same time are cross-correlated and cross-correlated to obtain the dimensionless instantaneous imaging value until the time returns to zero.

[0091] In the process of acquiring the imaging values ​​described above, the instantaneous imaging values ​​obtained at each moment are saved, and the instantaneous imaging values ​​at all moments are superimposed along the time direction to form the reverse time migration imaging values ​​of the current common shot point gather seismic record, which are calculated according to equation (5):

[0092] (5)

[0093] In the formula, i_shot(t) represents the instantaneous imaging value at time t, and I_shot represents the reverse time migration imaging value of the current common shot point gather seismic record.

[0094] In this embodiment, as Figure 5 As shown, this is the reverse time migration imaging value of the current common shot point gather seismic record obtained through cross-correlation imaging.

[0095] Step 6: Parallelize wave field propagation calculations based on GPU architecture;

[0096] The most time-consuming part of the reverse time migration imaging algorithm is the wave field propagation calculation, which requires differential calculation of the huge spatial grid at each time step to obtain the wave field derivative. The computational load increases exponentially with the model size, time sampling, and differential order.

[0097] Based on this, this invention parallelizes the wave field propagation calculation algorithm using GPUs, especially CUDA thread architecture, and maps thread coordinates to global cells through thread block indices and thread indices. During computation, the entire two-dimensional spatial mesh is divided into multiple thread blocks responsible for processing sub-regions within the mesh. Each thread within a block corresponds to one or more specific mesh points, directly mapped to mesh coordinates through thread indices. Before computation, thread blocks collaboratively load the required data from global memory into shared memory, and the threads then access the data from shared memory, reducing the number of repeated accesses to global memory. During computation, each thread reads the wave field values ​​of its corresponding mesh point and neighboring points, and then independently performs local difference calculations, directly updating the physical quantities it is responsible for. This transforms the spatial derivative calculations, which traditionally require multiple nested loops on a CPU, into a parallel operation executed simultaneously by tens of thousands of threads.

[0098] The pseudocode for the parallel algorithm for calculating the horizontal derivative of the 2Nth-order forward differential wavefield in a spatial grid is shown in Table 1. Here, threadIdx.x and threadIdx.y are the indices of the current thread, blockIdx.x and blockIdx.y are the indices of the current thread block, blockDim.x and blockDim.y are the dimensions of the thread block, and BLOCK_SIZE_Y and BLOCK_SIZE_X are the set thread block sizes. The vertical derivative of the wavefield is calculated in a similar manner and applied to the forward propagation of the source wavefield, the inverse reconstruction of the source wavefield, and the inverse propagation of the detector wavefield.

[0099] Table 1. Pseudocode for the parallel algorithm for calculating the horizontal derivative of the 2Nth-order forward differential wavefield in a spatial grid.

[0100]

[0101] The computation time of reverse time migration imaging of single common shot gather seismic records under parallelization was statistically analyzed under different parallel configurations. The computation time is the average time of 20 shot data. The results are shown in Table 2, which shows that efficient computation of reverse time migration imaging was achieved by GPU parallel technology.

[0102] Table 2. Statistical Table of Parallelization Algorithm Performance under Different Parallel Configurations

[0103]

[0104] Step 7: Perform reverse time migration imaging on the seismic records of multiple common shot point gathers, obtain the reverse time migration imaging values ​​corresponding to each gather, and stack them to obtain the final stacked reverse time migration imaging values, thus completing the pre-stack reverse time migration imaging.

[0105] Step 5 above describes the process of performing reverse time migration imaging on a single common shot point gather seismic record. The same process is used to perform reverse time migration imaging on multiple common shot point gather seismic records to obtain the reverse time migration imaging value corresponding to each common shot point gather seismic record. Since the calculation method is the same as that for reverse time migration imaging of a single common shot point gather seismic record, it will not be described again here.

[0106] Based on this, in order to suppress random noise, highlight the true structural response, and improve the signal-to-noise ratio and lateral continuity of the imaging results, the spatially consistent point values ​​in the reverse time migration imaging values ​​corresponding to the seismic records of each common shot point gather are algebraically summed to obtain the final stacked reverse time migration imaging values, thus completing the pre-stack reverse time migration imaging.

[0107] In this embodiment, reverse time migration (RTM) imaging was performed on the seismic records of 100 common shot gathers, generating corresponding RTM imaging values ​​for each shot. The RTM imaging values ​​of the 100 shots were then superimposed to form a stacked RTM imaging value with high signal-to-noise ratio and clear structural features, as shown below. Figure 6 As shown.

[0108] Step 8: Perform multi-stage processing on the superimposed reverse time migration imaging values, including depth gain and Laplace filtering-morphological filtering-Laplace filtering, to obtain the final high-precision migration imaging profile.

[0109] To address the issue of weak deep signal energy in the migration imaging profile, a deep gain processing is applied to the superimposed reverse-time migration imaging values ​​to enhance the imaging intensity and continuity of deep reflection structures, effectively compensating for amplitude attenuation caused by wavefront propagation and medium absorption.

[0110] Deep gain processing is calculated according to equation (6):

[0111] (6)

[0112] In the formula, I(z, x) represents the stacked reverse time migration image value at spatial location (z, x), I g (z, x) represents the stacked reverse-time migration image value at spatial location (z, x) after deep gain; q is the gain factor, typically greater than 1; z max This indicates the maximum vertical position of the superimposed reverse time-shifted imaging values.

[0113] In this embodiment, a depth gain processing is applied to the superimposed reverse-time migration imaging values ​​to enhance the imaging intensity and continuity of deep structures, such as... Figure 7 As shown.

[0114] Subsequently, Laplace filtering was applied to the superimposed reverse-time migration (RTM) images of the deep gain to suppress high-frequency noise generated by the numerical calculation of the two-way wave equation during RTM imaging. This significantly reduced high-frequency interference components while preserving the main structural features. Building upon this, morphological processing methods were introduced, employing top-hat and bottom-hat filtering to enhance the Laplace-filtered results. The top-hat filter extracts local bright structures from the migrated imaging profile, while the bottom-hat filter enhances the contrast between dark structures and the background. The synergistic effect of both effectively separates structural information at different scales and gray levels, significantly improving the clarity of structural boundaries and the ability to identify geological features.

[0115] To further optimize the morphological processing results, Laplace filtering is applied to the output results again to suppress secondary noise that may be introduced during morphological operations, thereby further enhancing the clarity of the reflective interface and the overall quality of the imaging results.

[0116] The above process is calculated according to formula (7):

[0117] (7)

[0118] In the formula, I f These are the superimposed reverse-time migration image values ​​after multi-stage filtering. Represents the Laplace operator. • denotes the morphological opening operation, and • denotes the morphological closing operation.

[0119] In this embodiment, as Figure 8 The image shown is the result of multi-level filtering after superimposing reverse-time migration imaging values; where, Figure 8 Part (a) in the image represents the image after the first Laplace filter; Figure 8 Part (b) in the image represents the image after morphological filtering based on the first Laplace filter; Figure 8 Part (c) in the image represents the image after Laplace filtering based on morphological filtering.

[0120] Through the aforementioned deep gain and the cascaded application of multi-level filtering, a high-precision migration imaging profile with low noise level, clear structural boundaries, and rich geological information is finally obtained, laying a reliable foundation for subsequent structural interpretation.

[0121] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0122] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, then this invention should also include these modifications and variations.

[0123] Although embodiments of the present invention have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those skilled in the art can make modifications, alterations, substitutions and variations to the above embodiments within the scope of the present invention.

Claims

1. A pre-stack reverse time migration imaging method based on GPU parallel technology and boundary storage strategy, characterized in that, Includes the following steps: Step 1: Based on the data of the study area, set the basic model parameters of reverse time migration imaging and the absorption boundary condition parameters of the convolution perfectly matched layer, and determine the absorption boundary region; Step 2: Solve the first-order velocity-sound pressure acoustic wave equation using the high-order staggered grid finite difference method, perform forward propagation calculation of the source wavefield, and apply the convolution fully matched layer CPML absorbing boundary condition in the absorbing boundary region to absorb boundary reflections. At the same time, save the wavefield data in the absorbing boundary region at each time step, as well as the source wavefield data at the final time. Step 3: Reconstruct the source wavefield in reverse based on the saved wavefield data of the absorbing boundary region and the source wavefield data at the final moment; Step 4: Load the common shot gather according to the detector location, and perform reverse propagation of the detector wavefield using the same calculation method as the forward propagation of the source wavefield; Step 5: Simultaneously advance the reverse reconstruction of the source wavefield and the reverse propagation of the detector wavefield, perform cross-correlation calculations on the two sound pressure fields at the same time, and superimpose the calculation results at each time along the time direction to obtain the reverse time migration imaging value of the current common shot point gather; Step 6: Parallelize wave field propagation calculations based on GPU architecture; Step 7: Perform reverse time migration imaging on the seismic records of multiple common shot point gathers, obtain the reverse time migration imaging values ​​corresponding to each gather, and stack them to obtain the final stacked reverse time migration imaging values, thus completing the pre-stack reverse time migration imaging. Step 8: Perform depth gain and multi-level filtering on the superimposed reverse time migration imaging values ​​to obtain the final high-precision migration imaging profile.

2. The pre-stack reverse time migration imaging method based on GPU parallel technology and boundary storage strategy according to claim 1, characterized in that, In step 1, the basic model parameters of the study area are set based on the study area data, including at least the geometric model, density model, P-wave velocity model, spatial location parameters of the shot point and the detector, and time-frequency parameters. A perfectly matched convolutional layer is introduced as an absorption boundary condition. The number of CPML layers (npml), absorption coefficient (R), and maximum complex frequency shift parameter (α) are set according to the basic model parameters and time-frequency parameters. max Maximum tensile parameter κ max This allows us to determine the absorbing boundary region, set a high-order staggered mesh with a difference order of 2N, and calculate the difference coefficients a.

3. The pre-stack reverse time migration imaging method based on GPU parallel technology and boundary storage strategy according to claim 2, characterized in that, In step 2, the first-order velocity-sound pressure acoustic wave equation is solved using the high-order staggered grid finite difference method to calculate the forward propagation of the source wavefield. A fully matched convolutional layer (CPML) absorbing boundary condition is applied within the absorbing boundary region to absorb boundary reflections. Simultaneously, the wavefield data at each time step within the absorbing boundary region, as well as the source wavefield data at the final time step, are saved. Specifically: Wavefield simulation was performed based on the first-order velocity-sound pressure acoustic wave equation. Numerical discretization was performed using a high-order staggered grid finite difference scheme. A seismic source was loaded at the shot point location. The source wavefield velocity and sound pressure fields were calculated using basic model parameters. Furthermore, a convolutionally perfectly matched layer absorbing boundary condition was applied to the defined absorbing boundary region, calculated according to the following formula: ； ； ； In the formula, Let represent the sound pressure field at index (i, j) at time t. Let be the horizontal component of the velocity field at the position indexed (i, j+0.5) at time t-0.

5. Let the vertical component of the velocity field at the index (i+0.5, j) at time t-0.5 be the velocity field component. This represents the model density at index (i,j+0.5). , Let represent the wave velocity of the model P-wave at index (i,j), Δt represent the time sampling interval, and Δx and Δz represent the spatial sampling intervals. Let be the CPML convolution term applied to the horizontal component v of the velocity field at time t. Let w be the CPML convolution term applied to the vertical component w of the velocity field at time t. Let be the CPML convolution term applied to the horizontal component p of the sound pressure field at time t. κ is the CPML convolution term applied to the vertical component of the sound pressure field p at time t. x κ is the horizontal tensile parameter. z For vertical stretching parameters; A boundary storage strategy is adopted to save the wavefield data of the source wavefield, sound pressure field, and velocity field at each time step within the absorbing boundary region, as well as the source wavefield data at the final moment.

4. The pre-stack reverse time migration imaging method based on GPU parallel technology and boundary storage strategy according to claim 3, characterized in that, In step 2, the storage formula for the wavefield boundary under the staggered grid finite difference scheme is as follows: ； In the formula, v bd w bd px represents the part where boundary conditions are applied to the velocity field. bd This represents the part of the sound pressure field with boundary conditions applied in the horizontal direction, pz bd This indicates the part where boundary conditions are applied in the vertical direction of the sound pressure field. nx and nz are the total number of sampling points in each direction, and npml is the number of CPML layers set.

5. The pre-stack reverse time migration imaging method based on GPU parallel technology and boundary storage strategy according to claim 4, characterized in that, In step 3, the source wavefield is reconstructed in reverse based on the saved absorbing boundary region wavefield data and the source wavefield data at the final moment. Specifically: Using the source wavefield at the final moment saved during forward propagation as the initial condition for reverse reconstruction of the source wavefield, the source wavefield is reconstructed based on the first-order velocity-sound pressure acoustic wave equation. The absorption boundary region wavefield data saved during forward propagation of the source wavefield is used to replace the absorption boundary region of the reconstructed source wavefield in order to restore the complete source wavefield. The initial condition is the source wavefield at the final moment preserved during the forward propagation of the source wavefield, calculated according to the following formula: ； In the formula, v fd nt w fd nt p fd nt Represents the source wave field at the final moment, v re nt w re nt p re nt This indicates the final moment of reconstruction of the source wave field.

6. The pre-stack reverse time migration imaging method based on GPU parallel technology and boundary storage strategy according to claim 5, characterized in that, In step 4, the common shot gather is loaded according to the detector location, and the detector wavefield is propagated backward using the same calculation method as the forward propagation of the source wavefield. Specifically: Using the seismic records of each shot point gather loaded at the detector location as the source, the reverse propagation of the detector wavefield is calculated using the forward propagation calculation method of the source wavefield and the application of CPML absorbing boundary conditions, specifically according to the following formula: ； In the formula, (zrec, xrec) represents the detector location, Rec(zrec, xrec, t) represents the seismic record received by the detector at location (zrec, xrec) at time t, ρ(zrec, xrec) represents the model density at location (zrec, xrec), and v t+0.5 (zrec, xrec), w t+0.5 (zrec, xrec) represent the horizontal and vertical components of the velocity field at position (zrec, xrec) at time t+0.5, respectively.

7. The pre-stack reverse time migration imaging method based on GPU parallel technology and boundary storage strategy according to claim 6, characterized in that, In step 5, the reverse reconstruction of the source wavefield and the reverse propagation of the detector wavefield are carried out simultaneously. Cross-correlation calculations are performed on the two sound pressure fields at the same time, and the calculation results at each time are superimposed along the time direction to obtain the reverse-time migration imaging value of the current common shot point gather. Specifically: Based on the principle of time consistency, that is, taking advantage of the simultaneous existence of the source wave field as the incident wave and the detector wave field as the reflected wave at the reflection point, the reverse reconstruction of the source wave field and the reverse propagation of the detector wave field are synchronously promoted from the maximum time. The cross-correlation calculation of the two sound pressure fields at the same time is performed to obtain the dimensionless instantaneous imaging value until the time returns to zero. During the acquisition of imaging values, the instantaneous imaging values ​​obtained at each moment are saved. The instantaneous imaging values ​​at all moments are superimposed along the time direction to form the reverse-time migration imaging values ​​of the current common shot point gather seismic record, calculated according to the following formula: ； In the formula, i_shot(t) represents the instantaneous imaging value at time t, and I_shot represents the reverse time migration imaging value of the current common shot point gather seismic record.

8. The pre-stack reverse time migration imaging method based on GPU parallel technology and boundary storage strategy according to claim 7, characterized in that, In step 6, the wave field propagation calculation is parallelized based on the GPU architecture, specifically as follows: The GPU-based CUDA thread architecture parallelizes wave field propagation calculations by dividing a two-dimensional spatial grid into multiple thread blocks. Threads within a thread block collaborate to read data and perform parallel local difference calculations through shared memory.

9. The pre-stack reverse time migration imaging method based on GPU parallel technology and boundary storage strategy according to claim 8, characterized in that, In step 7, reverse time migration imaging is performed on the seismic records of multiple common shot point gathers, and the reverse time migration imaging values ​​corresponding to each gather are obtained and stacked to obtain the final stacked reverse time migration imaging values, thus completing the pre-stack reverse time migration imaging. Specifically: Reverse-time migration imaging was performed on multiple common shot point gather seismic records to obtain the reverse-time migration imaging value corresponding to each common shot point gather seismic record. The spatially consistent point values ​​in the reverse-time migration imaging values ​​corresponding to each common shot point gather seismic record were algebraically summed to obtain the final stacked reverse-time migration imaging value, thus completing the pre-stack reverse-time migration imaging.

10. The pre-stack reverse time migration imaging method based on GPU parallel technology and boundary storage strategy according to claim 9, characterized in that, In step 8, the superimposed reverse-time migration imaging values ​​are subjected to depth gain and multi-level filtering to obtain the final high-precision migration imaging profile, specifically: First, depth gain processing is applied to the superimposed reverse time migration imaging values. Then, the imaging values ​​after depth gain are sequentially subjected to a first Laplace filter to suppress high-frequency noise, a morphological filter to enhance the structural boundaries, and a second Laplace filter to suppress secondary noise introduced by morphological operations, thereby obtaining a high-precision migration imaging profile with low noise levels and clear structural boundaries.

Images