Parallel forward method of multi-source anisotropic elastic wave based on one-dimensional grouped convolution

Abstract

The application provides an anisotropic multi-source elastic wave parallel forward method based on one-dimensional transverse grouping convolution, belongs to the technical field of geoscience and seismic detection service, and comprises the following steps: 1, anisotropic elastic wave model preparation; 2, multi-source same kind physical quantity storage reorganization; 3, parallel computing operator based on one-dimensional transverse grouping convolution; 4, parallel dispersion suppression filter design; and 5, anisotropic elastic wave forward numerical simulation precision and parallel computing efficiency test. In the multi-shot forward numerical simulation process, the same kind physical quantity matrix is reorganized into a three-dimensional tensor, parallel computing is carried out by combining a one-dimensional transverse grouping convolution operator, then a filter is used for synchronous processing of the wave field of multiple sources at a specific moment, so that the synchronous calculation of the wave field of multiple sources is realized in one sampling time. The method is a parallel scheme which takes into account the calculation efficiency, simulation precision and application cost.

Description

Technical Field

[0001] This invention belongs to the field of Earth Science and Earthquake Detection Service Technology, specifically relating to a parallel forward modeling method for multi-source anisotropic elastic waves based on one-dimensional grouped convolution. Background Technology

[0002] Elastic wave forward modeling is a crucial tool for studying seismic wave propagation, reverse time migration, and full waveform inversion, playing a vital role in research fields such as petroleum seismic exploration. Among various numerical solution methods, the finite difference method is widely used for elastic wave forward modeling due to its simplicity and high computational efficiency. However, in actual seismic imaging and inversion processes, to meet the wavefield data requirements for imaging condition construction and gradient calculation, it is usually necessary to perform forward modeling calculations on a large number of sources using the same model. This significantly increases the computational cost of forward modeling and becomes one of the main computational bottlenecks in seismic imaging and inversion.

[0003] In traditional finite-difference forward modeling, multi-source calculations typically employ a shot-by-shot approach. For each source, independent calculations of spatial derivatives and temporal step updates are required. Due to the highly similar numerical computational structures across different sources, this shot-by-shot approach results in repetitive execution in the temporal dimension. Furthermore, each simulation necessitates rescheduling and reactivating GPU computing resources, incurring additional hardware startup and scheduling overhead. This hinders operator fusion and data reuse within a unified computational graph, limiting the full utilization of GPU parallel computing capabilities. Therefore, researching accelerated forward modeling methods for elastic waves from multiple sources is crucial for improving the efficiency of elastic wave numerical simulations and supporting seismic exploration, microseismic monitoring, and other related seismic detection services.

[0004] Currently, methods for accelerating multi-shot forward modeling of elastic waves can be mainly divided into three categories: The first category is the multi-GPU parallel computing strategy. This method utilizes multiple GPUs or GPU clusters to distribute forward modeling tasks for different shot points across different computing devices, thereby achieving parallel computation of multi-shot forward modeling. By increasing the number of computing devices, this strategy can significantly improve overall computational efficiency. However, this method usually requires the design of an efficient multi-device communication and task scheduling mechanism to coordinate the allocation of computational tasks and result aggregation among different GPUs, thus ensuring the efficiency of parallel computation, which to some extent increases the complexity of system implementation. At the same time, while multi-source acceleration strategies based on multi-GPU collaboration can achieve high computational efficiency, they usually place higher demands on hardware resources and deployment costs. Therefore, research on accelerating multi-source elastic wave forward modeling using a single GPU still has significant practical implications. The second category is the source coding parallel computing strategy. This method can mathematically encode multiple source signals on a single GPU, combining multiple sources into a coded source, and simultaneously performing calculations for multiple shot points in a single forward modeling simulation, thereby significantly reducing the number of forward modeling iterations and improving computational efficiency. However, due to the superposition and propagation of multiple seismic sources in the same wavefield, crosstalk noise can easily arise between different shot points. Therefore, this type of method usually requires the design of reasonable source encoding or decoding strategies to reduce or suppress the impact of crosstalk. The third category is wavefield prediction methods based on deep learning. In recent years, deep learning has emerged in many fields and has also provided a disruptive technical path for realizing fast calculation of elastic wavefields on a single GPU. Many studies have proposed methods such as autoencoders and PINN-based neural networks, which use neural networks as function approximators to learn and simulate the physical processes of elastic wave propagation to achieve fast wavefield prediction. Compared with traditional direct numerical solutions, these methods can achieve significant speedups in specific tasks. However, their accuracy and generalization ability are still limited, especially when applied to complex elastic media or geological models outside the training distribution, and they often cannot completely replace high-fidelity finite difference forward modeling.

[0005] In summary, while existing methods can reduce the simulation time of multi-source seismic waves at different levels, achieving a balanced approach to elastic wave multi-shot forward modeling remains a valuable research area, especially for finite-difference-based physics-driven methods. Due to discretization errors, finite-difference forward modeling is inevitably affected by numerical dispersion; although increasing the difference order or refining the mesh can mitigate this problem, it usually incurs additional computational overhead. Therefore, how to improve the efficiency of multi-shot forward modeling while suppressing numerical dispersion at low cost remains a crucial issue.

[0006] Based on the above analysis, the main problems of current multi-source elastic wave forward modeling are: how to achieve efficient parallel computation of multiple sources, avoid coupling interference and application costs between multiple sources, and improve the accuracy of numerical simulation.

[0007] Therefore, it is essential to design and develop a method and system that can accelerate the forward modeling speed of multi-source anisotropic elastic waves at a lower application cost, while simultaneously improving the accuracy of finite difference numerical simulations with lower overhead. Summary of the Invention

[0008] The purpose of this invention is to provide a parallel forward modeling method for multi-source anisotropic elastic waves. The main contributions include: (1) Based on the commonalities in the calculation of wave field evolution of different sources by solving the first-order velocity-stress equations for each shot, the physical quantity matrices of the same type corresponding to different sources are organized along the source dimension to construct a unified three-dimensional tensor structure, thereby realizing centralized management of multi-shot data during the calculation process and improving the efficiency of data access and supply. (2) For high-dimensional tensor data structures, a parallel computing strategy based on grouped convolution is designed. In this structure, each two-dimensional sub-tensor in the three-dimensional tensor corresponds to an independent one-dimensional transverse convolution operator, which is used to solve the first-order partial derivative of space, thereby realizing the parallel independent calculation of multiple physical quantity matrices of the same type and improving the utilization rate of GPU computing units. (3) In order to suppress the artifacts caused by the finite difference decomposition, a plug-and-play filter is designed by combining the Nyquist sampling theorem and the Butterworth response. This filter can effectively suppress numerical artifacts in the wave field at a specific time while minimizing additional computational overhead and reducing the impact on the effective signal.

[0009] This invention provides the following technical solution:

[0010] A parallel forward modeling method for anisotropic multi-source elastic waves based on one-dimensional transverse grouping convolution includes the following steps:

[0011] Step (1): Establish an anisotropic elastic wave model:

[0012] The complex anisotropic elastic wave model is constructed using a two-dimensional parameter field containing various irregular formation media. For the stiffness coefficients C11, C13, C33, P-wave Vp, S-wave Vs, density rho, and other parameter fields in the model, the same set of coordinate mapping relationships is reused. Combined with interpolation, scaling, and clipping, the multi-parameter field is uniformly resampled from the original grid size to different target grid sizes. This processing method can effectively ensure that the spatial correspondence between the physical parameter fields remains unchanged after resampling, which facilitates the parallel computing efficiency test and numerical simulation effect verification under different scale velocity models.

[0013] Step (2): Storage and recombination of similar physical quantities from multiple earthquake sources:

[0014] Based on the first-order velocity-stress equation of anisotropic elastic waves, the same physical quantities corresponding to different earthquake sources in the same formula are uniformly organized to construct a high-dimensional tensor along the earthquake source dimension. Independent and parallel operators are used on each layer of tensor to realize the synchronous update of the same velocity field or stress field of multiple earthquake sources within the same time step. This data organization method reorganizes the matrix of the same physical quantities into a continuous data structure, thereby providing a stable data input to the GPU computing unit and forming a continuous data computation flow in the computation process.

[0015] Step (3): Parallel computation operator based on one-dimensional horizontal grouping convolution: In this structure, each two-dimensional sub-tensor in the three-dimensional tensor corresponds to an independent one-dimensional horizontal convolution operator, which is used to solve the first-order partial derivative in space, thereby realizing the parallel independent computation of multiple physical quantity matrices of the same kind and improving the utilization of GPU computing units.

[0016] Step (4): Design a parallel dispersion suppression filter:

[0017] To suppress numerical dispersion errors introduced during finite difference calculations, a plug-and-play parallel wavenumber domain low-pass filter was designed. This filter performs a two-dimensional Fourier transform on the high-dimensional tensor composed of multi-shot wavefields at a specific time, mapping the spatial domain wavefields to the wavenumber domain in batches. Based on the Nyquist sampling theorem, it filters out wavelengths smaller than [a certain value]. The components; the nth-order Butterworth response is used to minimize the Gibbs effect in the spatial domain, and finally the processed multi-source wavefield is synchronously restored to the spatial domain; since the wavefields of different sources are independent of each other;

[0018] Step (5): Testing the accuracy and parallel computing efficiency of forward numerical simulation of anisotropic elastic waves:

[0019] A PML layer of thickness c is wrapped around velocity models of different sizes, and in (X) 0+i Y 0+i Simultaneously inject a Ricker wavelet of h Hz at the location to excite the wave field; compare the numerical accuracy and computational efficiency of the proposed method and the traditional shot-by-shot convolution method before and after using the filter at the co-source location.

[0020] The specific steps in step (3) above are as follows:

[0021] Step (3.1) Convolution kernel setting:

[0022] Divide the finite difference coefficients by the spatial step size The values ​​are used as kernel weights for the one-dimensional horizontal convolution. The convolution size is set to 1×R according to the finite difference order R. Based on the number of earthquake sources being computed in parallel, multiple one-dimensional horizontal convolutions are stacked into one-dimensional horizontal grouped convolutions to ensure that each earthquake source has a one-dimensional horizontal convolution operator on its corresponding physical quantity matrix. The formula for the one-dimensional horizontal grouped convolution and the definition of the convolution kernel are as follows:

[0023]

[0024] In the formula: , These represent the row and column indices of the input and output matrices, respectively. Where is the convolution kernel length, i.e., the finite difference order, and dim represents the number of shots that need to be computed in parallel. This represents the position index of the one-dimensional horizontally grouped convolutional kernel and the input matrix in the row direction; Represents the output matrix; The output matrix represents the output matrix in The value at the position; Represents a one-dimensional horizontal convolution kernel; Represents the input matrix. These are the difference coefficients;

[0025] Step (3.2) Calculation of partial derivatives in the x-direction:

[0026] One-dimensional lateral grouping convolution is used to achieve parallel computation of the partial derivatives of a three-dimensional input matrix in the x-direction. A one-dimensional lateral convolution is applied to each submatrix of the input three-dimensional matrix. The convolution starts from the top left corner of each sub-tensor and slides row by row and column by column with a set stride. The computation is performed on a local region of the same size as the convolution kernel on each sub-tensor. For each sliding position, the input matrix in the group is multiplied element-wise with the corresponding convolution kernel and summed to obtain the partial derivative value at that position. Since there is no data interaction between the sub-tensors, this process can complete the calculation of the partial derivatives of multiple submatrices simultaneously in a single convolution operation, achieving efficient parallel computation.

[0027] Step (3.3) z-direction partial derivative principle:

[0028] To address the slow calculation of z-direction partial derivatives due to the low efficiency of GPUs reading non-contiguous data column-wise, this method uses matrix transposition to transform the originally non-contiguous data in the column direction into contiguous data in the row direction. Then, one-dimensional horizontal grouped convolution is used for calculation, and finally, the result is transposed back to complete the z-direction partial derivative calculation. This method does not change the actual storage location of the data, but improves data reading efficiency by changing the data reading index. The matrix transposition is defined as:

[0029]

[0030] In the formula: This represents the value in the i-th row and j-th column of the matrix, where i represents the row index and j represents the column index; m represents the number of rows in the matrix and n represents the number of columns in the matrix.

[0031] In step (4) above, the filter can achieve complete data parallelism at the source dimension during the filtering process, thereby ensuring overall processing efficiency; the formula is as follows:

[0032]

[0033] in, Represents Fourier transform, Represents the inverse Fourier transform. This represents the filter input wave field at a specific moment. This represents the output wave field of the filter at a specific moment. and These represent the horizontal and vertical directions, respectively. and The spatial distances from the walk in the x and z directions are respectively. and These represent the wavenumbers along the horizontal and vertical directions, respectively. Represents the Butterworth response order. This is the wave field safety factor.

[0034] The beneficial effects of this invention are as follows: This application provides a parallel forward numerical simulation method for multi-source anisotropic elastic waves that combines similar data matrix recombination, one-dimensional lateral grouping convolution, GPU data reading optimization, and filters, aiming to make full use of the computing units and cache resources of GPUs to reduce application costs.

[0035] First, by utilizing the commonalities in the calculation of wavefield evolution processes from different sources for each shot, the matrices of similar physical quantities corresponding to different sources are organized along the source dimension to construct a unified three-dimensional tensor structure. This enables centralized management of multi-shot data during the calculation process, thereby improving data access and supply efficiency.

[0036] Secondly, an independent one-dimensional transverse convolution operator is used on each two-dimensional subtensor in the three-dimensional tensor to solve the first-order partial derivative in space, thereby enabling parallel and independent computation of multiple physical quantity matrices of the same type and improving the utilization of GPU computing units.

[0037] Then, to suppress artifacts caused by finite difference decomposition, a plug-and-play parallel filter was designed. This filter can simultaneously and effectively suppress numerical artifacts in the wavefield at multiple specific times while minimizing additional computational overhead and reducing the impact on the effective signal.

[0038] Finally, the parallel computing efficiency and numerical simulation accuracy were verified on Hess models of different sizes. Compared with the current mainstream technologies, this invention improves the accuracy of finite difference numerical simulation while exhibiting high parallel computing efficiency and low application cost, providing reliable and efficient forward modeling technology support for subsequent full waveform inversion and reverse time migration. Attached Figure Description

[0039] Figure 1 This is a flowchart of Example 1;

[0040] Figure 2 These are model diagrams of different scales during the preparation of the anisotropic elastic wave velocity model in Example 1;

[0041] In the figure: a is the original size Vp model of the Hess model, b is the 400×400 size Vp model, and c is the 1000×1000 size Vp model;

[0042] Figure 3 This is a schematic diagram of the principle of storing and recombining similar physical quantities from multiple earthquake sources, as shown in Example 1.

[0043] In the figure: a) shows the physical quantity data corresponding to the shot-by-shot calculation method in a scattered layout in the GPU, and b) shows the physical quantity data corresponding to the method of this application in a continuous layout in the GPU.

[0044] Figure 4 This is a diagram illustrating the one-dimensional horizontal grouping convolution calculation principle of Example 1;

[0045] Figure 5 This is the complex anisotropic elastic wave parameter model structure of Example 2;

[0046] In the figure: a is the trimmed stiffness coefficient matrix C11, b is the trimmed stiffness coefficient matrix C13, c is the trimmed stiffness coefficient matrix C33, d is the trimmed stiffness coefficient matrix C44, e is the trimmed P-wave velocity model, f is the trimmed density model, g is the trimmed P-wave anisotropy intensity, h is the transition parameter between the trimmed P-wave and S-wave, i is the original stiffness coefficient matrix C11, j is the original stiffness coefficient matrix C13, k is the original stiffness coefficient matrix C33, l is the original stiffness coefficient matrix C44, m is the original P-wave velocity model, n is the original density model, o is the original P-wave anisotropy intensity, and p is the original P-wave and S-wave transition parameter.

[0047] Figure 6 This is a comparison chart of the numerical simulation accuracy of Example 2 and the shot-by-shot finite difference forward modeling method;

[0048] In the figure: a is the original Vz wave field snapshot of this method, b is the Vz wave field snapshot of this method + filter, c is the Vz wave field snapshot of the shot-by-shot method, d is the Vz wave field snapshot of the 20th order finite difference, e is the difference between a and c, f is the difference between b and d, g is the difference between c and d, and h is the difference between b and c.

[0049] Figure 7 This is a graph showing the computational efficiency of Example 2 and the shot-by-shot method in models with different numbers of shots and different sizes.

[0050] In the figure: a represents the average computation time of the method of this application and the shot-by-shot method under different numbers of shots and different velocity model sizes; b represents the average time saving and speedup ratio of the method of this application and the shot-by-shot method under different numbers of shots and different velocity model sizes. Detailed Implementation

[0051] To make the technical problems to be solved, the technical solutions, and the beneficial effects of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the present invention and are not intended to limit the present invention.

[0052] Example 1

[0053] As attached Figure 1 As shown in the flowchart, a parallel forward modeling method for multi-source anisotropic elastic waves based on one-dimensional grouped convolution includes the following steps:

[0054] Step (1): Preparation of anisotropic elastic wave velocity model:

[0055] The Hess model is used for the velocity model of complex anisotropic elastic waves. For the stiffness coefficients C11, C13, C33, P-wave Vp, S-wave Vs, and density rho, the same set of coordinate mapping relationships is used, combined with bilinear interpolation, scaling, and clipping to reconstruct the original model in multiple sizes. Starting with a 200×200 grid size, the original model is clipped by gradually increasing the number of grid points in the row and column directions, successively constructing sub-models of different sizes such as 200×200 and 400×400. When the row dimension increases to the maximum range of the original model, the model size is further expanded only along the column direction by increasing the number of grid points by 400 each time. The constructed velocity models of different scales of Vp are shown below. Figure 2 As shown, the model prepared by this method is used to verify the parallel computing efficiency and numerical simulation effect of the method of the present invention under different shot counts and different model sizes.

[0056] Step (2): Storage and recombination of similar physical quantities from multiple earthquake sources:

[0057] Based on the first-order velocity-stress equation for anisotropic elastic waves, similar physical quantities corresponding to different earthquake sources in the same formula are uniformly organized to construct a high-dimensional tensor along the source dimension. Independent and parallel operators are then used on each layer of the tensor to achieve synchronous updates of similar velocity or stress fields from multiple earthquake sources within the same time step. This data organization method reorganizes similar physical quantity matrices into a continuous data structure, thereby continuously providing stable data input to the GPU computing unit and enabling the computation process to form a continuous data computation flow.

[0058] Taking the stress fields Txx and Txz as examples, Figure 3 The data organization of the shot-by-shot calculation method and the method of this invention in multi-source forward modeling is demonstrated. Figure 3 In the image, 'a' represents the distribution of the Txx and Txz stress fields corresponding to different earthquake sources in the GPU memory during the shot-by-shot calculation process. Since each shot forward modeling task is executed sequentially, the stress field data of the same type are relatively scattered in the memory, with poor continuity, which is not conducive to subsequent unified call of operators for parallel calculation. Figure 3 In this invention, 'b' represents the data organization method. The Txx and Txz stress fields corresponding to multiple shots are reorganized into three-dimensional tensors of the form [Ns, Nx, Nz] according to the shot number dimension, where Ns represents the shot number. Each two-dimensional tensor of size [Nx, Nz] represents the stress field corresponding to a source location. Different two-dimensional tensors are physically independent of each other, but they are stored in memory in a centralized manner according to the same type of physical quantity, resulting in a more continuous data distribution. Subsequently, the same operator can be used on each two-dimensional tensor for synchronous and independent calculation, thereby realizing the simultaneous calculation of the partial derivatives of the Txx and Txz stress fields of multiple sources.

[0059] Step (3): Parallel computation operator based on one-dimensional horizontal grouping convolution:

[0060] (3.1) Convolution kernel setting:

[0061] Divide the finite difference coefficients by the spatial step size The values ​​are used as kernel weights for the one-dimensional horizontal convolution. The convolution size is set to 1×6 with a finite difference order R=6. Based on the number of earthquake sources being computed in parallel, multiple one-dimensional horizontal convolutions are stacked into one-dimensional horizontal grouped convolutions, ensuring that each earthquake source's corresponding physical quantity matrix has a one-dimensional horizontal convolution operator. The formula for one-dimensional horizontal grouped convolution and the definition of the convolution kernel are as follows:

[0062]

[0063] In the formula: , These represent the row and column indices of the input and output matrices, respectively. Where is the convolution kernel length (finite difference order), and dim represents the number of shots that need to be computed in parallel. This represents the position index of the one-dimensional horizontally grouped convolutional kernel and the input matrix in the row direction; Represents the output matrix; The output matrix represents the output matrix in The value at the position; Represents a one-dimensional horizontal convolution kernel; Represents the input matrix. These are the difference coefficients.

[0064] (3.2) Calculation of partial derivatives in the x-direction:

[0065] The parallel computation process of the partial derivatives of the 3D input matrix in the x-direction using one-dimensional lateral grouping convolution is as follows: Figure 4 As shown, a one-dimensional horizontal convolution is applied to each submatrix of the input 3D matrix. The convolution starts from the top-left corner of each sub-tensor and slides row-by-row and column-by-column with a set stride, performing calculations on a local region of the same size as the convolution kernel on each sub-tensor. For each sliding position, the input matrix is ​​element-wise multiplied with the corresponding convolution kernel weights and summed. The result is then saved to the corresponding output matrix, thus obtaining the partial derivative value at that position. Since there is no data interaction between the sub-tensors, this process can simultaneously calculate the partial derivatives of multiple sub-matrices in a single convolution operation, achieving efficient parallel computation.

[0066] (3.3) Partial derivative in the z-direction:

[0067] To address the slow speed of z-direction partial derivative calculation due to the low efficiency of GPUs reading non-contiguous data column-wise, this method uses matrix transposition to transform the originally non-contiguous data in the column direction into contiguous data in the row direction. Then, one-dimensional horizontal grouped convolution is used for calculation, and finally, the result is transposed back to complete the z-direction partial derivative calculation. This method does not change the actual storage location of the data, but improves data reading efficiency by changing the data reading index. Matrix transposition is defined as:

[0068]

[0069] In the formula: This represents the value in the i-th row and j-th column of the matrix, where i represents the row index and j represents the column index. m represents the number of rows in the matrix, and n represents the number of columns in the matrix.

[0070] Step (4): Design a parallel dispersion suppression filter:

[0071] To suppress numerical dispersion errors introduced during finite difference calculations, a plug-and-play parallel wavenumber domain low-pass filter was designed. This filter performs a two-dimensional Fourier transform on the high-dimensional tensor composed of multi-shot wavefields at a specific time, mapping the spatial domain wavefields to the wavenumber domain in batches. Based on the Nyquist sampling theorem, it filters out wavelengths smaller than... The filter employs an nth-order Butterworth response to minimize the Gibbs effect in the spatial domain, and then synchronously restores the processed multi-source wavefield to the spatial domain. Since the wavefields of different sources are independent of each other, this filtering process can achieve complete data parallelism along the source dimension, thus ensuring overall processing efficiency. The principle of this filter is as follows:

[0072]

[0073] in, Represents Fourier transform, Represents the inverse Fourier transform. This represents the filter input wave field at a specific moment. This represents the output wave field of the filter at a specific moment. and These represent the horizontal and vertical directions, respectively. and The spatial distances from the walk in the x and z directions are respectively. and These represent the wavenumbers along the horizontal and vertical directions, respectively. Represents the Butterworth response order. This is the wave field safety factor.

[0074] Step (5): Testing the accuracy and parallel computing efficiency of forward numerical simulation of anisotropic elastic waves:

[0075] Hess models of varying sizes were surrounded by 100-layer PML layers, and several shot sources were sequentially injected at the model's center, left, and right sides, with the wave field excited using a 40 Hz Ricker wavelet. The numerical accuracy and computational efficiency of the proposed method at the same source location and the shot-by-shot method were compared before and after using filters.

[0076] Example 2:

[0077] This embodiment uses the method of Embodiment 1 to verify the parallel computing efficiency and numerical simulation accuracy on Hess models of different sizes.

[0078] The experimental platform configuration for this embodiment is as follows: AMD EPYC 7773X, NVIDIA A100, Ubuntu 24.04 operating system, CUDA 11.7 acceleration, Python 3.9 experimental environment, and PyTorch 2.0.1 deep learning framework. The specific process of this embodiment is as follows:

[0079] 2.1 Preparation of Hess models of different sizes:

[0080] The system reads .mat files containing parameters such as stiffness coefficients C11, C13, C33, longitudinal wave Vp, transverse wave Vs, and density rho into the system and performs a consistency check on each parameter field to ensure that their mesh size and spatial correspondence remain consistent. Based on this, the same set of coordinate mapping relationships is reused, and bilinear interpolation, scaling, and clipping are used to reconstruct the original Hess model in multiple dimensions. Starting with a 200×200 mesh size, the original model is clipped by gradually increasing the number of mesh points in the row and column directions, successively constructing sub-models of different sizes such as 200×200 and 400×400. When the row dimension reaches the maximum range of the original model, the model size is further expanded only in the column direction by increasing the number of mesh points by 400 each time, thus obtaining a series of complex anisotropic elastic wave velocity models of different scales. Because the parameter fields reuse a unified coordinate mapping relationship during processing, it can be guaranteed that the parameter fields C11, C13, C33, Vp, Vs, rho, etc., after reconstruction and cropping, will still maintain a strict spatial correspondence. Finally, the processed parameter fields are saved, and the results are as follows. Figure 5 As shown, this is used for subsequent parallel computing efficiency testing and numerical simulation effect verification under different shot counts and model sizes. Figure 5 The ah sub-image in the image represents the processed model, while ip represents the original size model. From... Figure 5 As can be seen, the processed parameter fields are consistent with their corresponding original parameter fields in terms of spatial distribution, structural morphology, and stratigraphic variation.

[0081] The steps for preparing Hess models of different sizes are as follows, with the input being a Hess model of the original size and the output being a Hess model of a specified size:

[0082] Step 1: Read the multiple parameter fields of the original Hess model.

[0083] fields <- Load(C11, C13, C33, C44, Vp, Vs, rho).

[0084] Step 2: Check that all field sizes are consistent. Check_Size_Consistency(fields).

[0085] Step 3: Construct a unified coordinate mapping based on the target size.

[0086] coords <- Build_Coordinates(original_size, target_size).

[0087] Step 4: Resample all parameter fields using bilinear interpolation and downsampling methods.

[0088] for each field in fields do. resampled_field <- Interpolate(field,coords, method). end for.

[0089] Step 5: Check if the correspondence between the multiple parameters is consistent. Check_Size_Consistency(resampled_field).

[0090] Step 6: Output the new Hess model file. Save_To_Output(resampled_field).

[0091] 2.2 Storage and recombination of similar physical quantities from multiple earthquake sources:

[0092] Based on the first-order velocity-stress equation for anisotropic elastic waves, similar physical quantities corresponding to different earthquake sources in the same formula are uniformly organized to construct a high-dimensional tensor along the source dimension. Independent and parallel operators are then used on each layer of the tensor to achieve synchronous updates of similar velocity or stress fields from multiple earthquake sources within the same time step. This data organization method reorganizes similar physical quantity matrices into a continuous data structure, thereby continuously providing stable data input to the GPU computing unit and enabling the computation process to form a continuous data computation flow.

[0093] For example, the Txx stress fields corresponding to six earthquake sources, each with a 200x200 grid size, are stacked into a three-dimensional tensor of the form [6, 200, 200]. Each two-dimensional tensor of size [200, 200] represents the stress field corresponding to a specific earthquake source location. Subsequently, the same operator is used to perform synchronous and independent calculations on each two-dimensional tensor, thereby enabling the simultaneous calculation of the partial derivatives of the Txx stress fields from multiple earthquake sources.

[0094] 2.3 One-dimensional horizontal grouping convolution setting and calculation:

[0095] 2.3.1 Setting of one-dimensional horizontal grouping convolution kernel:

[0096] Divide the 6th-order finite difference coefficients by the space step size. The values ​​after that are used as the kernel weights for the one-dimensional horizontally grouped convolution. The convolution kernel weights are:

[0097]

[0098] 2.3.2 Calculation of partial derivatives in the x-direction:

[0099] One-dimensional horizontal grouping convolution is used to realize the x-square of a three-dimensional input matrix.

[0100] Parallel computation of partial derivatives. A one-dimensional transverse convolution is applied to each submatrix of the input 3D matrix. The convolution starts from the top left corner of each subtensor and slides row by row and column by column with a set stride, performing computation on a local region of the same size as the convolution kernel on each subtensor. For each sliding position, the input matrix within the group is element-wise multiplied with the corresponding convolution kernel and summed to obtain the partial derivative value at that position. Since there is no data interaction between the subtensors, this process can simultaneously complete the calculation of partial derivatives of multiple submatrices in a single convolution operation, achieving efficient parallel computation.

[0101] 2.3.3 Partial derivative in the z-direction:

[0102] To address the slow speed of z-direction partial derivative calculation due to the low efficiency of GPUs reading non-contiguous data column-wise, this method uses matrix transposition to transform the originally non-contiguous data in the column direction into contiguous data in the row direction. Then, one-dimensional horizontal grouped convolution is used for calculation, and finally, the result is transposed back to complete the z-direction partial derivative calculation. This method does not change the actual storage location of the data, but improves data reading efficiency by changing the data reading index. Matrix transposition is defined as:

[0103]

[0104] In the formula: This represents the value in the i-th row and j-th column of the matrix, where i represents the row index and j represents the column index. m represents the number of rows in the matrix, and n represents the number of columns in the matrix.

[0105] 2.4 Parallel Dispersion Suppression Filter:

[0106] To suppress numerical dispersion errors introduced during finite difference calculations, a plug-and-play parallel wavenumber domain low-pass filter was designed. This filter performs a two-dimensional Fourier transform on the high-dimensional tensor composed of multi-shot wavefields at a specific time, mapping the spatial domain wavefields to the wavenumber domain in batches. Furthermore, it filters out wavelengths smaller than [a certain value] based on the Nyquist sampling theorem. The component (α=4) is used. The n=4th order Butterworth response is employed to minimize the Gibbs effect in the spatial domain, and finally, the processed multi-source wavefield is synchronously restored to the spatial domain. Since the wavefields of different sources are independent of each other, this filtering process can achieve complete data parallelism in the source dimension, thus ensuring overall processing efficiency. The filter is as follows:

[0107]

[0108] in, Represents Fourier transform, Represents the inverse Fourier transform. This represents the filter input wave field at a specific moment. This represents the output wave field of the filter at a specific moment. and These represent the horizontal and vertical directions, respectively. and These represent the wavenumbers along the horizontal and vertical directions, respectively. Represents the Butterworth response order. This is the wave field safety factor.

[0109] The calculation steps for the parallel dispersion suppression filter are as follows: its input is the original high-dimensional tensor of the multi-shot wavefield, and its output is the filtered high-dimensional tensor of the multi-shot wavefield.

[0110] Step 1: Generate wavenumber domain coordinates based on the input multi-shot wavefield

[0111] dx,dz,nx,nz <- reshape(input_fields)

[0112] kz <- FFT_Frequency(nz, dz)kx <- RFFT_Frequency(nx, dx)

[0113] Step 2: Calculate the two-dimensional wavenumber modulus based on the wavenumber domain coordinates.

[0114] kz_grid, kx_grid <- Meshgrid(kz, kx)k_mag <- Sqrt(kz_grid^2 + kx_grid^2)

[0115] Step 3: Calculate the truncated wavenumber range based on the two-dimensional wavenumber modulus and generate the mask code.

[0116] k_cutoff <- 1 / (4 ×Max(dx, dz))

[0117] mask <- 1 / (1 + (k_mag / k_cutoff)^8)

[0118] Step 4: Apply the mask code to the high-dimensional tensor of the input wavefield to suppress dispersion artifacts.

[0119] output_fields <- mask ×input_fields

[0120] 2.5 Accuracy and Parallel Computation Efficiency Tests of Anisotropic Elastic Wave Forward Modeling Numerical Simulation:

[0121] Hess models of varying sizes were surrounded by 100-layer PML layers, and several shot sources were sequentially injected at the model's center, left, and right sides, sampling the Ricker wavelet excitation wavefield at 40 Hz. The numerical accuracy and computational efficiency of the proposed method and the shot-by-shot convolution method at the same source location were compared before and after using filters.

[0122] 2.6 Implementation Results:

[0123] use Figure 5 The ah model shown is surrounded by a 100-grid-thick PML layer. The seismic source is positioned at the center, and a single 40 Hz Ricker wavelet is injected with a sampling time of 5 seconds. The accuracy of single-run numerical simulations is compared between the conventional shot-by-shot method and the parallel convolution method of this invention. The results are as follows... Figure 6 As shown.

[0124] from Figure 6 The difference in Vz wavefield snapshots between the original method and the shot-by-shot method shows that the two methods perform similarly in terms of numerical simulation accuracy. Figure 6 The 20th-order finite-difference wavefield snapshot shown in d is used as a reference. Even without a filter, noticeable artifacts still exist in the wavefield (see [reference]). Figure 6 (Points 1 and 2 marked in c); after introducing the filter of this invention, the artifacts at the corresponding positions are significantly suppressed (see...). Figure 6 (at points 1 and 2 in b). Figure 6 f in Figure 6 The expression `g` in the figure gives the difference distributions relative to the 20th-order finite difference results for the cases with and without filters, respectively. Compared to... Figure 6 g in Figure 6 The differences at points 1 and 2 marked in f are significantly smaller, indicating that the filter effectively reduces artifact errors. Furthermore, Figure 6 The value of 'h' in the figure represents the difference in wavefield before and after filtering. It can be seen that the filter suppresses artifacts while having minimal impact on the effective signal. In summary, the filter proposed in this invention maintains numerical accuracy while exhibiting good numerical artifact suppression performance.

[0125] The seismic source position is then set at the center, and 6, 10, and 14 shot sources are sequentially excited on the left and right sides of the center, using a 40Hz Ricker wavelet and a sampling time of 5 seconds. The computational efficiency of the traditional shot-by-shot method and the parallel convolution method of this invention is compared. Figure 7 As shown.

[0126] Figure 7 Figure 'a' shows the comparison between the shot-by-shot method and the method of this invention on an A100 GPU with different numbers of shots, at speed 401. 2The computation time varies across models ranging from mesh size to full mesh size. The time taken for the 6-shot sequential method ranges from 29.43s to 150.87s, while the time taken for parallel computation of 6 shots using the method of this invention ranges from 6.23s to 113.85s. The time taken for the 10-shot sequential method ranges from 50.61s to 259.35s, while the time taken for parallel computation of 10 shots using the method of this invention ranges from 9.02s to 186.09s. The time taken for the 14-shot sequential method ranges from 75.64s to 366.65s, while the time taken for parallel computation of 14 shots using the method of this invention ranges from 11.84s to 257.90s. As can be seen from the figure, when the model size is less than 1201... 2 When meshing, the computation time of the shot-by-shot method remains largely consistent with the document, while for model sizes exceeding 1201... 2 After meshing, the overall computation time gradually increases with the increase in the computational scale of a single shot. In contrast, although the computation time of the method in this invention continues to increase with the increase in model size, it is always lower than that of the shot-by-shot method.

[0127] exist Figure 7 Figure b further illustrates the relationship between the computation time difference and speedup ratio between the two methods as model size changes. As the model size increases, the computation time difference represented by the bar chart generally shows a trend of first decreasing and then increasing. The computation time differences for shots 6, 10, and 14 decrease from the initial 23.19s, 41.58s, and 63.79s to 1201s. 2 The lowest values ​​for the grid size were 9.08s, 13.91s, and 24.12s, which then gradually increased to 37.02s, 73.25s, and 108.75s. Meanwhile, the speedup ratio, represented by the line graph, gradually decreased from 4.72, 5.6, and 6.3, reaching its lowest value of 1.20, 1.18, and 1.21 at the 1701×1901 grid size, before gradually increasing to 1.32, 1.39, and 1.45.

[0128] The main reason for this phenomenon is related to the utilization of GPU computing resources and changes in memory access pressure. When the model mesh size is less than 1200², the A100 GPU still has a large margin in computing units and cache resources. The time consumption of the shot-by-shot serial computation method remains basically unchanged, while the time consumption of the method of this invention continues to increase. Therefore, the time difference between the two methods gradually narrows. As the model size increases to the range of 1200² to 1701×1901, the computational scale of the shot-by-shot method gradually increases, and its computation time begins to rise slowly. However, because the parallel method processes multiple data simultaneously, the memory access pressure increases more significantly, and the time consumption growth trend is more obvious, which further narrows the performance gap between the two methods, resulting in a continuous decrease in the speedup ratio. When the model size exceeds 1701×1901, the shot-by-shot method needs to repeat large-scale calculations shot by shot, and its computation time growth rate is significantly higher than that of the method of this invention. The parallel method can make fuller use of GPU computing resources by processing multiple shot data simultaneously, so the overall computational advantage is revealed again, and the speedup ratio shows an upward trend again.

[0129] In summary, the method of this invention can more fully leverage the advantages of parallel computing and more effectively utilize the GPU's cache and computing resources, thereby exhibiting more stable and superior overall performance in large-scale model computing tasks.

Claims

1. A parallel forward modeling method for anisotropic multi-source elastic waves based on one-dimensional transverse grouping convolution, characterized in that: Includes the following steps: Step (1): Establish an anisotropic elastic wave model: The complex anisotropic elastic wave model is constructed using a two-dimensional parameter field containing various irregular formation media. The same set of coordinate mapping relationships is reused for the stiffness coefficients C11, C13, C33, P-wave Vp, S-wave Vs, and density rho parameter fields in the model. Combined with interpolation, scaling, and clipping, the multi-parameter field is uniformly resampled from the original grid size to different target grid sizes. Step (2): Storage and recombination of similar physical quantities from multiple earthquake sources: Based on the first-order velocity-stress equation of anisotropic elastic waves, the same physical quantities corresponding to different sources in the same formula are uniformly organized to construct a high-dimensional tensor along the source dimension. Independent and parallel operators are used on each layer of tensor to realize the synchronous update of the same velocity field or stress field of multiple sources in the same time step. Step (3): Parallel computation operator based on one-dimensional horizontal grouping convolution: In this structure, each two-dimensional sub-tensor in the three-dimensional tensor corresponds to an independent one-dimensional horizontal convolution operator, which is used to solve the first-order partial derivative in space, thereby realizing the parallel independent computation of multiple physical quantity matrices of the same kind and improving the utilization of GPU computing units. Step (4): Design a parallel dispersion suppression filter: To suppress numerical dispersion errors introduced during finite difference calculations, a plug-and-play parallel wavenumber domain low-pass filter was designed. This filter performs a two-dimensional Fourier transform on the high-dimensional tensor composed of multi-shot wavefields at a specific time, mapping the spatial domain wavefields to the wavenumber domain in batches. Based on the Nyquist sampling theorem, it filters out wavelengths smaller than [a certain value]. The components; the nth-order Butterworth response is used to minimize the Gibbs effect in the spatial domain, and finally the processed multi-source wavefield is synchronously restored to the spatial domain; since the wavefields of different sources are independent of each other; Step (5): Testing the accuracy and parallel computing efficiency of forward numerical simulation of anisotropic elastic waves: A PML layer of thickness c is wrapped around velocity models of different sizes, and in (X) 0+i Y 0+i Simultaneously inject a Ricker wavelet of h Hz at the location to excite the wave field; compare the numerical accuracy and computational efficiency of the proposed method and the traditional shot-by-shot convolution method before and after using the filter at the co-source location.

2. The parallel forward modeling method for anisotropic multi-source elastic waves based on one-dimensional transverse grouping convolution according to claim 1, characterized in that: The specific steps in step (3) are as follows: Step (3.1) Convolution kernel setting: Divide the finite difference coefficients by the spatial step size The values ​​are used as kernel weights for the one-dimensional horizontal convolution. The convolution size is set to 1×R according to the finite difference order R. Based on the number of earthquake sources being computed in parallel, multiple one-dimensional horizontal convolutions are stacked into one-dimensional horizontal grouped convolutions to ensure that each earthquake source has a one-dimensional horizontal convolution operator on its corresponding physical quantity matrix. The formula for the one-dimensional horizontal grouped convolution and the definition of the convolution kernel are as follows: In the formula: , These represent the row and column indices of the input and output matrices, respectively. Where is the convolution kernel length, i.e., the finite difference order, and dim represents the number of shots that need to be computed in parallel. This represents the position index of the one-dimensional horizontally grouped convolutional kernel and the input matrix in the row direction; Represents the output matrix; The output matrix represents the output matrix in The value at the position; Represents one-dimensional horizontal Convolution kernel; represents the input matrix, These are the difference coefficients; Step (3.2) Calculation of partial derivatives in the x-direction: One-dimensional lateral grouping convolution is used to achieve parallel computation of the partial derivatives of a three-dimensional input matrix in the x-direction. A one-dimensional lateral convolution is applied to each submatrix of the input three-dimensional matrix. The convolution starts from the top left corner of each sub-tensor and slides row by row and column by column with a set stride. The computation is performed on a local region of the same size as the convolution kernel on each sub-tensor. For each sliding position, the input matrix in the group is multiplied element-wise with the corresponding convolution kernel and summed to obtain the partial derivative value at that position. Since there is no data interaction between the sub-tensors, this process can complete the calculation of the partial derivatives of multiple submatrices simultaneously in a single convolution operation, achieving efficient parallel computation. Step (3.3) z-direction partial derivative: To address the slow calculation of z-direction partial derivatives due to the low efficiency of GPUs reading non-contiguous data column-wise, this method uses matrix transposition to transform the originally non-contiguous data in the column direction into contiguous data in the row direction. Then, one-dimensional horizontal grouped convolution is used for calculation, and finally, the result is transposed back to complete the z-direction partial derivative calculation. This method does not change the actual storage location of the data, but improves data reading efficiency by changing the data reading index. The matrix transposition is defined as: In the formula: This represents the value in the i-th row and j-th column of the matrix, where i represents the row index and j represents the column index. m represents the number of rows in the matrix, and n represents the number of columns in the matrix.

3. The parallel forward modeling method for anisotropic multi-source elastic waves based on one-dimensional transverse grouping convolution according to claim 1, characterized in that: In step (4), the filter can achieve complete data parallelism at the source dimension during the filtering process, thereby ensuring overall processing efficiency; the formula is as follows: in, Represents Fourier transform, Represents the inverse Fourier transform. This represents the filter input wave field at a specific moment. This represents the output wave field of the filter at a specific moment. and These represent the horizontal and vertical directions, respectively. and The spatial distances from the walk in the x and z directions are respectively. and These represent the wavenumbers along the horizontal and vertical directions, respectively. Represents the Butterworth response order. This is the wave field safety factor.

Images