Local MQRBF-FD elastic wave propagation simulation method based on multi-scale modeling and parallel optimization

Abstract

The invention discloses a local MQRBF-FD elastic wave propagation simulation method based on multi-scale modeling and parallel optimization. The method comprises the following steps: 1, based on sine function combination characteristics of elastic wave Fourier expansion, optimizing shape parameters of an MQRBF by adopting a parallel improved random walk algorithm, and training an Adam-BP neural network model; 2, dividing a computational domain into a plurality of local sub-domains, and completing high-precision spatial discretization of a wave field in each sub-domain; 3, discretizing time by adopting a finite difference method, and finishing time integration in combination with a central difference method; and 4, converting wave field evolution into a solving problem of a sparse linear equation set, and updating a global wave field value through efficient iteration to complete high-precision simulation of elastic wave propagation. According to the method, the precision and calculation efficiency of elastic wave propagation simulation in the complex heterogeneous medium are remarkably improved, and the problems that a traditional finite difference method is large in numerical error and low in calculation efficiency under the complex boundary condition and the heterogeneous medium are effectively solved.

Description

Technical Field

[0001] The present invention relates to the technical field of acoustic wave propagation simulation, and in particular to a local MQRBF-FD (Multi-Quadratic Radial Basis Function-Finite Difference Method) elastic wave propagation simulation method based on multi-scale modeling and parallel optimization. Background Art

[0002] When traditional numerical methods simulate elastic wave propagation problems, they have significant limitations when facing complex boundary conditions and heterogeneous media. For example, the Finite Difference Method (FD) is widely used because of its simple algorithm and high computational efficiency. However, when dealing with complex geometric structures or heterogeneous media, it often requires a significant increase in grid resolution to ensure accuracy, which leads to a sharp increase in computational costs. However, such methods face severe challenges in scenarios such as deep-earth resource exploration and shale gas development in China: the global grid encryption strategy results in an exponential growth in the computational scale, making it difficult to meet the real-time requirements of large-scale three-dimensional simulations; at the same time, in areas where the wave velocity gradient changes violently (such as fault zones or fracture-developed areas), the FD method is prone to numerical oscillations caused by grid distortion, severely restricting the high-precision imaging ability of seismic inversion software. In addition, the rigid constraint of the FD method on the time step limits its application in cross-scale simulations in kilometer-level deep target areas, with significant cumulative errors and difficulty in supporting the refined requirements of deep resource exploration.

[0003] In recent years, the Radial Basis Function (RBF) has shown unique advantages in elastic wave propagation problems due to its grid-free characteristics. The Finite Difference Method based on Radial Basis Function (MQRBF-FD) combines the flexibility of the radial basis function in irregular geometries with the high efficiency of the finite difference method in time discretization, providing a new idea for solving elastic wave propagation problems in complex media. However, there are still two core bottlenecks in the existing technology:

[0004] First, the optimization of the shape parameter depends on empirical trial and error or global algorithms, and it is difficult to dynamically adapt to local medium characteristics, resulting in a sharp increase in interpolation errors at heterogeneous interfaces and unable to meet the high-fidelity reconstruction requirements of the wave field in complex structural areas;

[0005] Second, the traditional parallel framework lacks a collaborative mechanism for multi-scale modeling and sparse solution, and its efficiency is insufficient in the GPU / CPU heterogeneous computing power environment, restricting the potential release of high-performance computing platforms in large-scale grid simulations.

[0006] The above problems have seriously restricted the practical application of elastic wave numerical simulation technology in complex medium simulation. This invention breaks through the traditional technical framework and pioneers a local MQRBF-FD method based on the fusion of multi-scale modeling and intelligent optimization. Through the collaborative innovation of shape parameter adaptive learning model, heterogeneous parallel acceleration architecture and multi-resolution dynamic coupling mechanism, it significantly improves the accuracy and efficiency of complex medium simulation, providing core technical support for national strategic needs such as deep earth resource exploration and major engineering safety monitoring. Summary of the invention

[0007] In order to overcome the above technical problems, the present invention proposes a local MQRBF-FD elastic wave propagation simulation method based on multi-scale modeling and parallel optimization. This method dynamically optimizes shape parameters through an improved random walk algorithm (IRW) and an Adam-BP neural network model, and combines adaptive mesh refinement (AMR) and an efficient time integration strategy. It significantly improves the accuracy and computational efficiency of elastic wave propagation simulation in complex heterogeneous media, and effectively solves the problems of large numerical errors and low computational efficiency of traditional finite difference methods under complex boundary conditions and heterogeneous media.

[0008] In order to achieve the above object, the technical solution adopted by the present invention is:

[0009] The local MQRBF-FD elastic wave propagation simulation method based on multi-scale modeling and parallel optimization includes the following steps:

[0010] Step 1: Based on the sine function combination characteristics of elastic wave Fourier expansion, a parallel improved random walk algorithm is used to optimize the shape parameters of MQRBF, and the "LCSF-MQSP900K" data set containing the sine function combination and its optimal shape parameters is generated. The Adam-BP neural network model is trained based on this data set.

[0011] Step 2: Based on the optimized Adam-BP neural network model, a multi-scale space partitioning strategy is adopted to optimize the shape parameters in parallel according to the local characteristics of each subdomain, and the MQRBF interpolation method is used to complete the high-precision spatial discretization of elastic waves;

[0012] Step 3: Based on the spatial discretization result of step 2, the finite difference method (FD) is used to discretize the time, and the central difference method (CDM) is combined to complete the time integration;

[0013] Step 4: Based on the time discretization of step 3 and the space discretization of step 2, the evolution of the elastic wave field is transformed into the problem of solving a set of sparse linear equations; finally, the global wave field value is updated through efficient iteration to complete the high-precision simulation of elastic wave propagation.

[0014] The step 1 is specifically as follows:

[0015] The propagation behavior of elastic wave fields is described by the two-dimensional elastic wave equation:

[0016]

[0017] In the local subdomain, the wave velocity is set to be approximately uniform (v(x,y)≈v0), and the solution of equation (1) can be expanded into a Fourier series:

[0018]

[0019] Among them A mn and B mn are coefficients related to the frequency components of the fluctuations, corresponding to the amplitudes of the cosine and sine terms, ω mn is the angular frequency of the wave, which determines the vibration frequency of the wave field. x and L y is the length of the computational domain in the x and y directions, affecting the spatial scale of the wave mode. m and n are the indices of the spatial mode, controlling the number of modes of the wave field in the x and y directions, respectively. Formula (2) is derived from the homogeneous medium solution of equation (1) and characterizes the vibration characteristics of the local wave field. The formula represents the solution of the wave field as a sine function of multiple different frequency and mode combinations by Fourier expansion.

[0020] In complex media, due to the drastic spatial variation of the wave velocity v(x,y), the global Fourier expansion cannot effectively represent the wave field. Therefore, the local Fourier expansion is used to divide the complex domain into multiple small subdomains. In each subdomain, the wave velocity is assumed to be approximately a constant v0. In order to screen the dominant frequency ω mn , define the energy density:

[0021]

[0022] Formula (3) is based on the sum of the squares of the Fourier coefficients of formula (2), the dominant frequency ω dom corresponds to the maximum energy density.

[0023] For different frequency components ω mn The improved random walk (IRW) algorithm is used to optimize the corresponding MQRBF shape parameter c, with the goal of minimizing the maximum error MaxError(c):

[0024] MaxError(c)=max x |f(x)-s(x,c)| (4).

[0025] Formula (4) compares the actual wave field f(x) with the MQRBF interpolation result s(x,c).

[0026] The optimization process of the IRW algorithm for the shape parameter c is as follows:

[0027] Step (1) Initialization: Set the initial shape parameter c0, the number of iterations i and k, and define the step size accuracy θ and the error tolerance ∈;

[0028] Step (2) Generate random vectors: In each iteration, generate N random vectors u ik , whose distribution interval is (-q,q), is normalized and used to update the shape parameter c;

[0029] Step (3) Error comparison: Calculate each shape parameter c i The corresponding maximum error MaxError(c i ), and the error MaxError(c i-1 ) for comparison; if MaxError(c i ) <MaxError(c i-1 ), then update the shape parameters; if the error is larger and the step size λ exceeds the threshold θ, then halve the step size and start again;

[0030] Step (4) Termination condition: When MaxError(c) is less than ∈ or reaches the maximum number of iterations, stop the optimization and output c opt .

[0031] Through the above optimization process, the IRW algorithm can minimize the maximum error and optimize the accuracy of MQRBF interpolation. The IRW algorithm generates the "LCSF-MQSP900K" data set containing 900,000 sets of optimal parameters. Based on this data set, the Adam-BP neural network model is trained to output the dynamically optimized optimal shape parameters c opt , and successfully built an adaptive optimization model for shape parameters. Finally, the model was applied in various subdomains to dynamically adjust the shape parameters, thereby further improving the numerical stability and simulation accuracy in complex media.

[0032] The second step dynamically adjusts the grid density through an adaptive grid refinement (AMR) strategy: the grid is refined in areas with large velocity gradients, and the grid is coarsened in areas with smooth velocity gradients to optimize computing resources;

[0033] In each subdomain, the elastic wave equation solution is expanded into a combination of sine functions optimized in step 1. The MQRBF interpolation is performed in parallel with the dynamic shape parameters optimized by the model. The spatial variation of the wave field is captured through the neighborhood nodes, and a high-precision spatial derivative discrete matrix is ​​generated to complete the global spatial discretization, laying the foundation for time integration.

[0034] The step 2 is specifically as follows:

[0035] The computational domain is divided into multiple local subdomains, and the mesh density is dynamically adjusted through adaptive mesh refinement (AMR):

[0036]

[0037] Equation (5) serves the local Fourier expansion requirement of Equation (2). The greater the wave speed gradient , the denser the grid. Within each subdomain, MQRBF interpolation achieves high-precision calculation of spatial derivatives through neighborhood nodes:

[0038]

[0039] where k represents the number of neighborhood nodes, (x i , y i ) are the coordinates of neighborhood nodes, and w i (t) is the time-dependent interpolation weight. The numerical stability is enhanced through the polynomial extension matrix (Equation (12)). The interpolation calculations are independently completed for each subdomain, significantly improving the parallel efficiency. Finally, the spatial discretization of the wave field is completed, providing an accurate spatial discretization result for subsequent time integration and wave field evolution.

[0040] The specific content of Step 3 is as follows:

[0041] After completing the spatial discretization, the finite difference method (FD) is used for the discretization in the time direction. The second-order time derivative of the wave field is discretized using the second-order central difference method (CDM), and its expression is:

[0042]

[0043] where and represent the wave field values at the current, previous, and next steps respectively, △t is the time step, and Equation (7) substitutes the spatial derivative (Equation 6) into the time domain and dynamically adjusts the step size in combination with the adaptive time step strategy:

[0044]

[0045] where ∈ actual is the current error and ∈ tol is the preset tolerance. Equation (8) ensures that the step size is reduced to improve the accuracy when the wave field changes violently. The time integration is independently completed for each subdomain, optimizing the calculation efficiency.

[0046] The specific content of Step 4 is as follows:

[0047] Based on the spatio-temporal discretization results of Step 3, the wave field evolution is transformed into a sparse linear equation system:

[0048] Ax n+1 =b n (9)

[0049] The matrix A in Equation (9) is jointly constructed by the spatial derivative matrix in Equation (6) and the time discretization term in Equation (7).

[0050] The preconditioned conjugate gradient method (PCG) combined with incomplete LU decomposition (ILU) is used for parallel solution:

[0051] Preconditioning acceleration: The ILU decomposition improves the matrix condition number;

[0052] Boundary synchronization: Based on the local wave field values obtained during the time integration process, by directly sharing the wave field values of adjacent subdomain boundary nodes, the continuity of the wave field between subdomains is ensured.

[0053] Efficient iteration: PCG minimizes the residual ||Ax n+1 -b n ||, and each subdomain is solved independently.

[0054] Through the efficient solution of the sparse linear equations, the solution of the current time step of the global wave field is successfully updated, and the propagation behavior of elastic waves in complex media is fully simulated, finally achieving high-precision and high-efficiency elastic wave propagation simulation.

[0055] Advantages of the present invention:

[0056] Through an innovative synchronization mechanism, the present invention realizes the seamless connection and global consistency of the wave field between subdomains. Based on the multi-level parallel architecture design, the calculation time is significantly shortened, making the efficient evolution of the wave field possible. Combining the cooperative optimization strategy of the preconditioned conjugate gradient method (PCG) and incomplete LU decomposition (ILU), and the efficient parallel solution framework, the present invention demonstrates excellent computational efficiency and numerical accuracy in large-scale elastic wave propagation simulations, providing a stable and reliable solution for elastic wave propagation in complex media.

[0057] Through the above steps, the present invention constructs a complete elastic wave propagation simulation system. Its core innovation lies in: dynamically optimizing the shape parameters of multi-quadratic radial basis functions (MQRBF) through an Adam-BP neural network, combining multi-scale modeling, parallel computing, and adaptive time step strategies to comprehensively improve the computational accuracy, numerical stability, and overall simulation efficiency. Specifically, through the parallel improved random walk (IRW) algorithm, it quickly converges to the global optimal shape parameters, effectively creating a large-scale dataset of shape parameters, breaking through the efficiency bottleneck of traditional empirical parameter tuning; based on the wave velocity gradient-driven adaptive mesh refinement (AMR) technology, the mesh is accurately refined in areas with drastic changes in wave velocity (such as strong reflection interfaces) and coarsened in stable areas to achieve the optimal allocation of computing resources; through the full-process parallel acceleration framework, from mesh generation to the solution of sparse equations, the tasks of each subdomain are independently assigned to heterogeneous computing units, significantly reducing the calculation time of large-scale three-dimensional simulations.

[0058] In the present invention, the multi-scale modeling technology is realized through an adaptive mesh refinement strategy, which can accurately capture the wave field characteristics in regions with drastic wave speed changes, while reducing computational redundancy in stable regions, thereby optimizing the overall computational resource allocation. Through parallel processing, the mesh generation, shape parameter optimization, local interpolation, and solution of sparse linear equations for each sub-domain can be carried out independently, and the "ghost node" technology is combined to ensure data synchronization between sub-domains and the continuity of the wave field. This framework combining parallel and multi-scale modeling enables efficient and accurate elastic wave propagation simulation under complex boundary conditions, heterogeneous media, and irregular meshes.

[0059] In traditional methods, the selection of shape parameters relies on experience or trial and error, making it difficult to dynamically adapt to changes in local medium characteristics, which easily leads to numerical instability or reduced accuracy. The shape parameter adaptive optimization model constructed by the Adam-BP neural network in the present invention can dynamically adjust the shape parameters according to local wave field characteristics without manual intervention, ensuring high accuracy and stability in numerical simulation. Based on the optimized shape parameters, local MQRBF interpolation avoids the high computational cost of global mesh refinement and can efficiently capture the detailed changes in the wave field, especially in regions with drastic wave speed changes or significant medium heterogeneity. Combined with the adaptive time step strategy, the time step is dynamically adjusted, reducing the step size during drastic wave field changes to ensure accuracy and increasing the step size in stable regions to optimize efficiency, significantly improving the computational performance of long-term simulation.

[0060] Parallel computing technology runs through the entire process of the present invention, including mesh generation, shape parameter optimization, local interpolation, and solution of sparse linear equations. The computational tasks of each sub-domain can be completed independently, and boundary data is synchronized through the "ghost node" technology to ensure the continuity and consistency of the global wave field. In the solution of sparse linear equations, the preconditioned conjugate gradient method (PCG) and incomplete LU decomposition (ILU) are combined to effectively reduce the solution cost and significantly accelerate the numerical simulation of large-scale three-dimensional elastic wave propagation.

[0061] The present invention has successfully constructed an efficient, stable, and complex medium-adaptive elastic wave propagation simulation framework, providing a high-precision numerical tool for fields such as seismic exploration, ultrasonic imaging, and engineering simulation. Brief Description of the Drawings

[0062] Figure 1 It is a schematic flow chart of the local MQRBF-FD elastic wave propagation simulation method based on multi-scale modeling and parallel optimization.

[0063] Figure 2 It is the velocity model diagram of the embodiment.

[0064] Figure 3 It is the propagation behavior of the distribution field of the embodiment at multiple receiving points.

[0065] Figure 4 It is the comparison of wavefront propagation under different methods of the embodiments.

[0066] Figure 5 It is the comparison result between the method of the present invention and the finite difference (FD) method at five receiving points in the embodiments. Detailed implementation manners

[0067] The present invention will be further described in detail below with reference to the accompanying drawings.

[0068] The local MQRBF-FD elastic wave propagation simulation method of the present invention based on multi-scale modeling and parallel optimization realizes an efficient and high-precision simulation scheme of elastic waves based on an adaptive optimization shape parameter model, local interpolation technology and parallel computing.

[0069] The present invention will be introduced in detail in the following steps, and the specific operation process is as Figure 1 shown:

[0070] Step 1: Adaptive optimization of MQRBF shape parameters

[0071] In the Fourier expansion, the solution of the elastic wave equation can be expressed as a combination of frequency components of sine functions by formula (1). Since it is necessary to optimize the shape parameters corresponding to a large number of different sine function combinations, the optimization tasks can be assigned to multiple computing units by parallelizing the IRW algorithm introduced above, and the shape parameter search and optimization can be performed on multiple sine function combinations simultaneously, greatly shortening the optimization time.

[0072] The main steps of the parallel IRW algorithm are as follows:

[0073] Parallel step 1: In step 2 of the IRW algorithm, generating N random vectors u ik and the process of normalizing them can be executed simultaneously on multiple processors or threads, thus significantly reducing the time required for generating random vectors;

[0074] Parallel step 2: The optimization of the shape parameters involves calculating the maximum error MaxError(c) of each candidate shape parameter c, and this calculation can be completed in parallel by multiple computing units. Each processor can evaluate the errors of different shape parameters respectively, thus accelerating the convergence of the algorithm;

[0075] Parallel step 3: The two-part strategy for updating the shape parameters (i.e., step size adjustment and error comparison) can also be parallelized. When evaluating multiple candidate solutions, decisions can be made on different candidate solutions simultaneously, improving the overall optimization efficiency.

[0076] Through the above parallelization strategy, the IRW algorithm optimized the shape parameters of 900,000 sine function combinations and generated the "LCSF-MQSP900K" dataset. Each data instance encapsulates the frequency combination ω1, ω2, …, ω n and its corresponding optimal shape parameter c opt . Using this dataset, the backpropagation (BP) algorithm was trained and optimized with the adaptive moment estimation (Adam) optimizer to establish an adaptive optimization model for the shape parameters. By inputting the dominant frequency of the local subdomain, this model dynamically outputs the optimal shape parameter, providing adaptive parameter support for high-precision interpolation in subsequent steps.

[0077] The network structure of the Adam-BP model includes an input layer, two hidden layers, and an output layer:

[0078] - Input layer: Input the dominant frequency ω mn ;

[0079] - Hidden layer: Each layer contains 100 neurons, and the activation function is ReLU;

[0080] - Output layer: Output the corresponding optimal shape parameter c opt .

[0081] The key parameter settings for the training process are as follows:

[0082] α = 0.003, batch size = 128, number of training epochs = 50

[0083] The Adam-BP model is optimized by minimizing the mean squared error (MSE) loss function, and the formula is:

[0084]

[0085] The training results of the final model show that the predicted mean squared error is 0.163512, and the prediction accuracy reaches 97.98%. Thus, we established an adaptive optimization model for the shape parameters. By applying this model in parallel, the shape parameters can be dynamically adjusted within each subdomain, further improving the numerical stability and simulation accuracy in complex media.

[0086] Step 2: Local subdomain division and high-precision spatial discretization.

[0087] To accurately simulate the elastic wave propagation behavior in complex media, this step will use the optimization model constructed in Step 1 to further complete the high-precision spatial discretization of the wave field within each subdomain.

[0088] First, the computational domain is divided into multiple local subdomains. The division of local subdomains is based on the wave velocity distribution and spatial variation characteristics of the medium. The adaptive mesh refinement (AMR) technology is used to significantly improve the resolution in areas with drastic wave velocity changes (such as strong reflection interfaces or high gradient medium areas), while a coarser grid is used in areas with gentle wave velocity changes to optimize computing resources. Grid refinement dynamically adjusts the grid node spacing through the wave velocity gradient. The specific calculation is formula (5). In areas with drastic wave velocity changes, is larger, △x and △y are reduced, thus achieving mesh encryption; in the area where the wave velocity changes slowly, The smaller the value, the larger the △x and △y are, thus reducing the amount of calculation. The spatial distribution of the grid nodes satisfies the following relationship:

[0089] x i =x0+i△x,y j =y0+j△y (11)

[0090] Among them, i, j are grid indices, and △x, △y are position-dependent grid spacings.

[0091] After completing the local subdomain division and mesh refinement, this step uses the shape parameters optimized in step 1 to further complete the high-precision spatial discretization of the wave field in each subdomain. Specifically, combined with the parallelization strategy, the MQRBF interpolation method independently completes the interpolation calculation in each local subdomain, and each computing unit is responsible for the node interpolation of its corresponding subdomain, thereby significantly improving the computational efficiency. In the calculation of spatial derivatives, the derivative form of the radial basis function is used for the first-order and second-order derivatives.

[0092] In order to further improve the interpolation accuracy in the area where the wave field changes dramatically, the present invention introduces an extended matrix representation combining polynomial basis functions and radial basis functions:

[0093]

[0094] Among them, Φ is the radial basis function matrix, P is the polynomial basis function term, and λ is the regularization parameter used to improve numerical stability.

[0095] Through the above process, the MQRBF interpolation method finally completes the spatial discretization of the wave field, provides accurate discrete results for subsequent calculations in the time direction, and ensures high-precision simulation of elastic wave propagation in complex media.

[0096] Step 3: Time discretization and optimization of adaptive time step strategy. Through the adaptive time step strategy, the step size is dynamically adjusted according to the local error: when the wave field changes drastically, the step size is reduced to improve accuracy, and when it changes gently, the step size is increased to optimize efficiency. Each subdomain completes interpolation and time integration independently to generate the local wave field solution of the current time step, ensuring computational efficiency and numerical stability.

[0097] After the spatial discretization is completed in Step 2, the time discretization adopts the finite difference method (FD) to discretize the second-order time derivative term in the elastic wave equation (see Equation (7)). By combining Equation (7) with the spatial derivative discretization result (Equation (6)) in Step 2, a spatio-temporal coupled discrete equation is constructed, providing an explicit recurrence basis for wave field evolution.

[0098] To improve the computational efficiency and take into account numerical stability, the present invention introduces an error-driven adaptive time step strategy.

[0099] The specific process is as follows:

[0100] Local error estimation: Calculate the local error ∈ at the current time step within each subdomain actual ,

[0101] Dynamic step size adjustment: According to the ratio of the error to the preset tolerance ∈ tol , use Equation (8) to update the next time step size. When the wave field changes violently (∈ actual > ∈ tol ): Reduce the step size to improve the accuracy; when the wave field is stable (∈ actual < ∈ tol ): Increase the step size to optimize the efficiency.

[0102] The time integration tasks of each subdomain are independently completed through distributed computing. In addition, the maximum time step size is restricted by the CFL condition (Equation 13) to avoid numerical instability caused by too large a step size:

[0103] where v

[0104]

[0105] is the maximum wave speed in the medium, and △x and △y are determined by the adaptive mesh refinement in Step 2 (Equation 5). max

[0106] Through parallel local interpolation and time integration calculations, each subdomain independently completes its internal numerical simulation tasks, effectively improving the overall computational efficiency while ensuring computational accuracy.

[0107]

[0108]

[0108] Step 4: Parallel sparse linear equation system solution and wave field evolution. This equation system is jointly constructed by the MQRBF spatial derivative matrix and the time discretization term. The right-hand side contains the current wave field values and boundary data. The equation system is solved in parallel by the preconditioned conjugate gradient method (PCG) combined with incomplete LU decomposition (ILU): First, the sparse matrix is preconditioned to accelerate convergence, and then based on the local wave field solutions in Step 3, the boundary data of adjacent subdomains are directly shared to ensure wave field continuity.After the space-time discretization in Steps 2 and 3, the wave field evolution is transformed into a problem of solving a sparse linear equation system, whose mathematical form is Formula (9). To efficiently solve this equation system, the preconditioned conjugate gradient method (PCG) combined with incomplete LU decomposition (ILU) is adopted as the core solver. The preconditioning matrix M -1 is constructed by ILU decomposition:

[0109] M = L·U (14)

[0110] where L and U are lower triangular and upper triangular matrices that preserve the sparsity of the original matrix, significantly improving the matrix condition number and reducing the number of iterations. The linear equation systems of each subdomain are allocated to heterogeneous computing units (CPU / GPU) for independent solution. The CPU is responsible for logical control and task scheduling, and the GPU accelerates intensive computing tasks (such as matrix decomposition and iterative operations), and hybrid parallel optimization is achieved by combining CUDA or OpenMP.

[0111] In addition, the boundary data between subdomains are synchronously updated through the "ghost node" technology. Ghost nodes are virtual nodes introduced at the boundaries of subdomains. These nodes do not directly participate in the calculation but are used to store the boundary data shared with neighboring subdomains. During the numerical solution process, ghost nodes ensure the continuity and consistency of the wave field between different subdomains through data synchronization with adjacent subdomains. Specifically, the wave field values of adjacent subdomains are stored in virtual nodes and dynamically updated:

[0112]

[0113] where represents the boundary data stored in the virtual node of subdomain i, is the actual boundary wave field value of the adjacent subdomain j. This formula ensures the continuity and consistency of the wave field between different subdomains by regularly updating the boundary data, avoiding phase errors caused by boundary discontinuities. The synchronization frequency is strictly matched with the time step, ensuring the time continuity of the global wave field evolution. The solution of the sparse linear equation system and the wave field update are efficiently completed within the global scope, and finally, a high-precision simulation of elastic wave propagation is achieved.

[0114] Through the above steps, a numerical method for simulating elastic wave propagation based on the local MQRBF-FD and Adam-BP optimization model is finally constructed, which can accurately describe the behavior of elastic wave propagation in complex media. The key MQRBF shape parameter is adaptively optimized by the Adam-BP model, and combined with multi-scale modeling, heterogeneous parallel computing, and ghost node technology, the present invention demonstrates excellent computational efficiency and numerical accuracy in large-scale three-dimensional elastic wave propagation simulations, providing reliable technical support for seismic exploration and engineering detection.

[0115] In the present invention, the multi-scale modeling technique is realized through an adaptive mesh refinement (AMR) strategy, which can accurately capture the characteristics of the wave field in areas with drastic changes in wave velocity, while reducing computational redundancy in stable areas, thereby optimizing the overall allocation of computational resources. Through parallel processing, the mesh generation, shape parameter optimization, local interpolation, and solution of sparse linear equations for each sub-domain can be carried out independently, and the "ghost node" technology is combined to ensure data synchronization between sub-domains and the continuity of the wave field. This framework combining parallel and multi-scale modeling significantly reduces the computational time of large-scale three-dimensional simulations while maintaining high accuracy and numerical stability.

[0116] In summary, the local MQRBF-FD elastic wave propagation simulation method based on multi-scale modeling and parallel optimization proposed in the present invention constructs an efficient and stable elastic wave propagation simulation framework through a training optimization model combined with efficient numerical calculation techniques. This framework can achieve dynamic optimization of shape parameters and high-precision interpolation under complex wave fields, accurately simulate wave propagation behaviors in complex media and irregular grids, and provide a high-precision numerical tool for fields such as seismic exploration, ultrasonic imaging, and non-destructive testing.

[0117] In Step 1, the shape parameters of MQRBF with a large number of sine function combinations are optimized through a parallel improved random walk (IRW) algorithm, and a "LCSF-MQSP900K" dataset containing sine function combinations and corresponding optimal shape parameters is established. Based on this dataset, an Adam-BP neural network is trained to construct a shape parameter adaptive optimization model to realize the dynamic adjustment of shape parameters in the local medium. In Step 2, the computational domain is divided into multiple local sub-domains. Using the adaptive mesh refinement (AMR) technique, the mesh is refined in areas with drastic changes in wave velocity (such as strong reflection interfaces) to improve the resolution, and the solution of the elastic wave equation is expressed as a combination of sine functions through Fourier expansion. Combining the optimization model in Step 1 to dynamically adjust the shape parameters, and finally completing the spatial discretization through parallel MQRBF interpolation. In Step 3, based on the spatial discretization results in Step 2, the second-order time derivative is discretized using the finite difference method (FD), and the time integration calculation is completed in combination with the central difference method (CDM). The time step is dynamically adjusted through an adaptive time step strategy, and the interpolation and time integration of each sub-domain are independently completed to obtain the local wave field solution at the current time step. Based on the time discretization in Step 3 and the spatial discretization in Step 2, the wave field evolution is transformed into a problem of solving a sparse linear equation system. Using the preconditioned conjugate gradient method (PCG) combined with incomplete LU decomposition (ILU), the sparse linear equation system is efficiently solved, and finally, a high-precision and high-efficiency simulation of elastic wave propagation in complex media is realized.

[0118] Example:

[0119] The locally MQRBF-FD elastic wave propagation simulation method based on multi-scale modeling and parallel optimization proposed by the present invention has been verified for elastic wave propagation in the complex Marmousi model. This model simulates the geological environment of a complex heterogeneous medium, demonstrating the ability of the method in dealing with complex wave fields and capturing medium heterogeneity.

[0120] The range of the simulation calculation domain is 2 km to 11 km in the horizontal direction and 0 km to 3 km in the vertical direction. The calculation grid uses 200 horizontal units and 150 vertical units, with grid spacings of 45 m and 20 m respectively. This resolution fully captures the geological heterogeneity while ensuring the calculation efficiency.

[0121] The elastic parameters of the medium are set as follows: The P-wave velocity (v p ) varies between 1500 m / s and 4500 m / s, corresponding to sedimentary layers and consolidated rock layers. The shear wave velocity (v s ) is calculated according to , and the density (ρ) is calculated according to the Gardner formula. These parameters construct an elastic medium model with significant heterogeneity.

[0122] The source is set at a horizontal position of 7 km and a depth of 40 m, and a Ricker wave with a dominant frequency of 5 Hz is used as the source waveform. The receivers are evenly distributed along the top surface, with a total of 100 receivers at a depth of 10 m, used to capture the time evolution of the wave field. The boundary condition adopts the perfectly matched layer (PML) technique to effectively absorb the boundary reflected waves and avoid artificial boundary interference.

[0123] The time step is determined according to the CFL condition as:

[0124]

[0125] where v max = 4500 m / s is the maximum P-wave velocity and CFL = 0.5.

[0126] The simulation results include the velocity model (as shown in Figure 2 ), the stress distribution (as shown in Figure 3 ), and the receiver records (as shown in Figure 4 ). The stress distribution and receiver waveforms of the complex wave field reveal the severe wave field interference and significant medium effects in the Marmousi model, verifying the accuracy and stability of the method of the present invention in capturing the complex behavior of the wave field.

[0127] To further evaluate the performance of the method, five receiving points are selected to compare the performance of the present invention and the finite difference (FD) method in terms of maximum amplitude, total energy, mean square error, and calculation time (see Figure 5)。The results show that this method is superior to the FD method in terms of maintaining the maximum amplitude and total energy, and can better retain waveform features and energy transfer. At the same time, the mean square error of this method is significantly lower than that of the FD method, with higher numerical accuracy, but the calculation time increases slightly. However, this computational overhead is acceptable in the context of a significant improvement in accuracy and energy conservation.

[0128] This embodiment verifies the excellent performance of the method of the present invention in complex heterogeneous media: through the deep integration of multi-scale modeling and intelligent parameter optimization, the nonlinear propagation behaviors such as multi-path scattering, interface refraction, and diffraction of the wave field in the Marmousi model are accurately analyzed. Compared with the conventional finite difference method, this method, through dynamic parameter adaptation and global parallel architecture, while maintaining the energy conservation characteristics, significantly improves the numerical stability and phase accuracy, providing a high-performance simulation tool that is autonomous and controllable for deep geological resource exploration and major engineering safety assessment.

Claims

1. A local MQRBF-FD elastic wave propagation simulation method based on multi-scale modeling and parallel optimization, characterized in that: The steps include: Step 1: Based on the sine function combination characteristics of elastic wave Fourier expansion, a parallel improved random walk algorithm is used to optimize the shape parameters of MQRBF, a data set containing the sine function combination and its optimal shape parameters is generated, and the Adam-BP neural network model is trained based on the data set; Step 2: Based on the optimized Adam-BP neural network model, a multi-scale space partitioning strategy is adopted to optimize the shape parameters in parallel according to the local characteristics of each subdomain, and the MQRBF interpolation method is used to complete the high-precision spatial discretization of elastic waves; Step 3: Based on the spatial discretization results of step 2, the time evolution of elastic waves is discretized using the finite difference method, and the time integration is completed in combination with the central difference method; Step 4: Based on the time discretization in step 3 and the space discretization in step 2, the elastic wave field evolution is converted into the solution of a sparse linear equation system. By efficiently iteratively updating the global wave field value, a high-precision simulation of elastic wave propagation is completed.

2. The local MQRBF-FD elastic wave propagation simulation method based on multi-scale modeling and parallel optimization according to claim 1 is characterized in that: The step 1 is specifically as follows: The propagation behavior of elastic wave fields is described by the two-dimensional elastic wave equation: In the local subdomain, the wave velocity is set to be approximately uniform (v(x,y)≈v0), and the solution of equation (1) can be expanded into a Fourier series: Among them A mn and B mn are coefficients related to the frequency components of the fluctuations, corresponding to the amplitudes of the cosine and sine terms, ω mn is the angular frequency of the wave, which determines the vibration frequency of the wave field, L x and L y are the lengths of the computational domain in the x and y directions, and m and n are the indices of the spatial pattern; The complex domain is divided into multiple small subdomains. In each subdomain, the wave velocity is assumed to be approximately a constant v0. To screen the dominant frequency ω mn , define the energy density: Formula (3) is based on the sum of the squares of the Fourier coefficients of formula (2), the dominant frequency ω dom corresponds to the maximum energy density; The improved random walk algorithm is used to optimize different frequency components ω mn The sine function combination corresponds to the MQRBF shape parameter c, and its goal is to minimize the maximum error MaxError(c): MaxError(c)=max x |f(x)-s(x,c)| (4); Formula (4) compares the actual wave field f(x) with the MQRBF interpolation result s(x,c).

3. The local MQRBF-FD elastic wave propagation simulation method based on multi-scale modeling and parallel optimization according to claim 2 is characterized in that: The optimization process of the IRW algorithm for the shape parameter c is as follows: Step (1) Initialization: Set the initial shape parameter c0, the number of iterations i and k, and define the step size accuracy θ and the error tolerance ∈; Step (2) Generate random vectors: In each iteration, generate N random vectors u ik , whose distribution interval is (-q,q), is normalized and used to update the shape parameter c; Step (3) Error comparison: Calculate each shape parameter c i The corresponding maximum error MaxError(c i ), and the error MaxError(c i-1 ) for comparison; if MaxError(c i ) <MaxError(c i-1 ), then update the shape parameters; if the error is larger and the step size λ exceeds the threshold θ, then halve the step size and start again; Step (4) Termination condition: When MaxError(c) is less than ∈ or reaches the maximum number of iterations, stop the optimization and output c opt .

4. The local MQRBF-FD elastic wave propagation simulation method based on multi-scale modeling and parallel optimization according to claim 1 is characterized in that: The second step dynamically adjusts the grid density through an adaptive grid refinement strategy: the grid is refined in areas with large wave velocity gradients, and the grid is coarsened in areas with smooth wave velocity gradients to optimize computing resources; In each subdomain, the elastic wave equation solution is expanded into a combination of sine functions optimized in step 1. The MQRBF interpolation is performed in parallel with the dynamic shape parameters optimized by the model. The spatial variation of the wave field is captured through the neighborhood nodes, and a high-precision spatial derivative discrete matrix is ​​generated to complete the global spatial discretization.

5. The local MQRBF-FD elastic wave propagation simulation method based on multi-scale modeling and parallel optimization according to claim 4 is characterized in that: The step 2 is specifically as follows: Divide the computational domain into multiple local subdomains and dynamically adjust the mesh density through adaptive mesh refinement: Wave velocity gradient The larger the value, the denser the grid. In each subdomain, MQRBF interpolation achieves high-precision calculation of spatial derivatives through neighborhood nodes: Among them, k represents the number of neighboring nodes, (x i ,y i ) are the coordinates of the neighborhood nodes, w i (t) is the time-dependent interpolation weight.

6. The local MQRBF-FD elastic wave propagation simulation method based on multi-scale modeling and parallel optimization according to claim 1 is characterized in that: The step three is specifically as follows: After completing the spatial discretization, the time direction is discretized using the finite difference method, and the time second-order derivative of the wave field The discretization is performed using the second-order central difference method, which is expressed as: in, and They represent the wave field values ​​of the current, previous step and next step respectively, △t is the time step, and formula (7) substitutes the spatial derivative (formula 6) into the time domain, and dynamically adjusts the step size in combination with the adaptive time step strategy: Among them, ∈ actual is the current error, ∈ tol The preset tolerance.

7. The local MQRBF-FD elastic wave propagation simulation method based on multi-scale modeling and parallel optimization according to claim 1 is characterized in that: The step 4 is specifically as follows: Based on the time-space discretization results of step 3, the wave field evolution is transformed into a sparse linear equation system: Ax n+1 =b n (9)。 8. The local MQRBF-FD elastic wave propagation simulation method based on multi-scale modeling and parallel optimization according to claim 7 is characterized in that: The preconditioned conjugate gradient method is combined with incomplete LU decomposition to solve in parallel: Preconditioning acceleration: ILU decomposition improves matrix condition number; Boundary synchronization: Based on the local wave field values ​​obtained during the time integration process, the wave field values ​​of the adjacent sub-domain boundary nodes are directly shared to ensure the continuity of the wave field between sub-domains; Efficient iteration: PCG minimizes residual || Ax n+1 -b n ||, each subdomain is solved independently; By efficiently solving the sparse linear equations, the current time step solution of the global wave field is updated, the propagation behavior of elastic waves in complex media is simulated, and the elastic wave propagation simulation is realized.

Images