High-efficiency forward algorithm of discrete fourier tdm viscoelastic wave equation

Abstract

The application provides a high-efficiency forward algorithm for discrete Fourier TDM viscoelastic wave equation, comprising the following steps: step 1, frequency domain modeling of an initial model; step 2, time modeling (TDM) by using an integration scheme with an explicit format; step 3, frequency domain modeling (DSM) based on a direct solver; step 4, frequency domain modeling (ISM) based on an iterative solver; step 5, frequency domain modeling (HSM) based on a hybrid solver; and step 6, comparison of the three methods. The high-efficiency forward algorithm for discrete Fourier TDM viscoelastic wave equation is a forward algorithm for time modeling based on an explicit integration scheme, can accurately propagate in any heterogeneous layer medium without instability, in a three-dimensional case, frequency response is extracted by combining a discrete Fourier transform, and three-dimensional viscoelastic wave equation is solved based on frequency domain modeling of a direct and iterative or hybrid solver.

Description

Technical Field

[0001] This invention relates to the field of seismic wave imaging modeling technology, and specifically to an efficient forward modeling algorithm based on the discrete Fourier Transform (TDM) viscoelastic wave equation. Background Technology

[0002] Modern oil and gas exploration began in the mid-19th century, during the Second Industrial Revolution. For over a century since then, oil and gas resources have been one of the most important energy sources for humankind. They are not only a major source of daily transportation and a primary raw material for chemical production, but also a vital strategic resource for nations. As a major consumer of oil and gas resources, China's demand is increasing year by year. The oil and gas we use includes both conventional and unconventional sources. China is rich in unconventional oil and gas resources, but the quality of these reservoirs is declining, leading to increased technological requirements. Common methods in actual oil and gas exploration include seismic exploration and sonic logging. These methods primarily utilize wave information propagating in the reservoir to conduct qualitative and quantitative studies of reservoir parameters and structure.

[0003] Classical hyperbolic acoustic and elastic wave equations cannot accurately characterize the geological models of this type of reservoir, and traditional numerical simulation algorithms struggle to handle complex coupled characterization equations. Therefore, the development of new, stable, and efficient algorithms is crucial to meet the demands of unconventional oil and gas reservoir exploration and development.

[0004] Chinese patent application CN202110609064.1 discloses a method, system, device, and medium for seismic inversion based on wave equations. The method includes the following steps: obtaining and determining the calculation parameters of a single-frequency sensitivity kernel function based on seismic observation parameters and initial velocity model parameters; generating a velocity model with random boundaries; performing forward modeling of the wave equation based on the calculation parameters and velocity model to obtain a synthetic seismic record and a frequency domain incident wave field; determining the travel-time associated source based on the synthetic seismic record, and determining the single-frequency gradient of the travel-time target functional based on the travel-time associated source; extracting the low wavenumber portion from the single-frequency gradient, and performing seismic velocity inversion based on the low wavenumber. Compared with wave equation travel-time tomography based on a wavefield reconstruction framework, this invention has significant advantages in terms of computational efficiency, memory usage, and algorithm complexity, and can be widely applied in the field of seismic velocity modeling technology.

[0005] Chinese patent application CN202110177528.6 discloses a method, apparatus, device, and storage medium for forward modeling of frequency-domain viscoelastic waves. The method includes: obtaining a one-dimensional viscoelastic wave equation in the frequency-wavenumber domain; converting it into an equivalent weak integral form and then discretizing it into a finite element control equation; applying a three-dimensional point force source and absorbing boundary conditions of the target area to the finite element control equation; solving the finite element control equation with the applied three-dimensional point force source and absorbing boundary conditions to obtain the frequency-wavenumber domain wavefield of the target area; determining the wavenumber set corresponding to the frequency-wavenumber domain wavefield according to an adaptive wavenumber sampling strategy; and performing a non-equally spaced two-dimensional inverse Fourier transform on the frequency-wavenumber domain wavefield corresponding to each wavenumber in the wavenumber set as input to obtain the corresponding three-dimensional elastic wavefield in the frequency domain. This invention can improve the accuracy of forward modeling of frequency-domain viscoelastic waves.

[0006] Chinese patent application CN201911315597.8 discloses a multi-wave joint pre-stack waveform inversion method. Pre-stack seismic inversion methods mainly fall into two categories: those based on ray tracing forward operators and those based on wave equation forward operators. The former is the most widely used pre-stack AVO inversion, while the latter is called pre-stack waveform inversion. Unlike pre-stack AVO inversion, pre-stack waveform inversion methods can better utilize the complex wavefield responses in the input data, reduce the processing difficulty of the input data, and improve the sensitivity of seismic records to the target inversion parameters. Currently, pre-stack waveform inversion mostly employs nonlinear inversion strategies, resulting in low computational efficiency. This invention proposes a linear strategy pre-stack inversion method to reduce computation time. Simultaneously, to improve inversion accuracy and stability, this invention adopts a joint inversion strategy using PP and PS multi-wave data, which effectively improves the inversion accuracy of the three parameters, especially for shear wave velocity and density estimation. Well logging model testing verifies the effectiveness of the method.

[0007] The existing technologies described above are significantly different from the present invention and have failed to solve the technical problem we want to address. Therefore, we have invented a new efficient forward modeling algorithm based on the discrete Fourier TDM viscoelastic wave equation. Summary of the Invention

[0008] The purpose of this invention is to provide an efficient forward modeling algorithm based on the Discrete Fourier Transform (TDM) viscoelastic wave equation that can accurately propagate in any heterogeneous layered medium without instability.

[0009] The objective of this invention can be achieved through the following technical measures: an efficient forward modeling algorithm based on the Discrete Fourier Transform (DFT) TDM viscoelastic wave equation, which includes:

[0010] Step 1: Perform frequency domain modeling on the initial model;

[0011] Step 2: Time modeling (TDM) using an integrated scheme with an explicit format;

[0012] Step 3: Calculate the frequency domain modeling DSM based on the direct solver;

[0013] Step 4: Calculate the frequency domain modeling ISM based on the iterative solver;

[0014] Step 5: Calculate the frequency domain modeling HSM based on the hybrid solver;

[0015] Step 6: Compare these three methods.

[0016] The objective of this invention can also be achieved through the following technical measures:

[0017] In step 1, for two-dimensional frequency domain FWI, the forward modeling problem is realized in the frequency domain by a direct solver. In the three-dimensional case, the optimal strategy for forward modeling is not obvious. It can be improved by using time-domain modeling (TDM) combined with discrete Fourier transform to extract the frequency response, or by using frequency domain modeling based on direct and hybrid or iterative solvers.

[0018] In step 1, the three-dimensional viscous acoustic wave equation is considered in the frequency domain:

[0019]

[0020] Where density is represented by ρ(x), bulk modulus by κ(x), and angular frequency by ω; the monochromatic pressure wave field and source are represented by p(x, ω) and s(x, ω), respectively; in the expression for bulk modulus, the inherent attenuation can be easily realized in the frequency domain using the complex wave velocity;

[0021] Equation (1) can be reproduced in matrix form as follows:

[0022] Ap = s (2)

[0023] The complex impedance matrix A depends on the angular frequency ω and the parameters κ and ρ.

[0024] In step 2, the partial differential operator is transformed into an algebraic operation using existing iterative numerical methods, which can be represented by matrices. The wave equation attempts to estimate the vector f in time through an explicit system, where the vector f is pressure, solid particle velocity, or fluid / solid particle velocity.

[0025] In step 2, the matrix expression is:

[0026]

[0027] Where xyz are the three-dimensional spatial coordinate parameters, and t is the time parameter;

[0028] The quality matrix M is a diagonal matrix;

[0029] The stiffness matrix A should be operated on back and forth in the spatial domain or the spectral domain to obtain the desired solution;

[0030] In controlled-source seismology, the source term S is usually a local point source, and the corresponding formula in the frequency domain is a generalization of the Helmholtz equation.

[0031]

[0032] Where xyz are the three-dimensional spatial coordinate parameters, and w is the angular velocity;

[0033] When considering finite discretization, the impedance matrix B is complex and has a symmetric mode.

[0034] In step 2, TDM is typically performed using an explicit time planning algorithm; at each time step, the solution for each spatial grid point is estimated from the solution of the previous time step.

[0035] In step 2, for frequency domain FWI, the core memory storage of the full time series is useless because the frequency domain wave field is obtained by discrete Fourier summation; by using the standard domain decomposition method, the computational domain is divided into finite-dimensional subdomains, which can effectively parallelize the time domain algorithm; in the framework of multi-source simulation, the low memory requirements of TDM also allow for coarse-grained parallelism of the sources, and if the number of processors is significantly greater than the number of sources, this parallelism can be combined with domain decomposition parallelism.

[0036] In step 2, the dimension of the 3D N3 computational grid is represented by n; if only parallelism on the source is implemented, the memory and time complexity of TDM are given; real 3D surveys require a large amount of memory to store N distributed across the processor. rhs Wave field, O(N) 3 N rhs )=O(N 5 ), where N rhs The number of sources is on one side of the grid; the TDM algorithm can accurately achieve propagation in any heterogeneous medium without instability, provided that a sufficiently fine time discretization method is used; despite the significant increase in memory requirements, the expansion of elasticity, anisotropy and decay modeling is possible.

[0037] In step 2, the dimension of the matrix is ​​the number of unknowns in the computation grid. The numerical bandwidth and the number of non-zero coefficient matrices depend on the template of the numerical discretization. Differential operators are embedded in the frequency domain. Wave modeling reduces the solution of large sparse linear equation systems and multiple right-hand side terms (RHS). Each RHS corresponds to a source, resulting in a linear system.

[0038] In step 2, time-domain modeling (TDM) of the explicit integration scheme is performed; within the framework of multi-source simulation, the low memory requirements of TDM also allow for coarse-grained parallelism on the sources, and if the number of processors is significantly greater than the number of sources, domain decomposition parallelism can be combined.

[0039] In step 3, by performing LU decomposition on the impedance matrix B for each frequency, the finite number of frequencies required for the frequency domain FWI can be effectively modeled for a large number of sources; the parallelism in the direct solution method DSM is achieved by using a large-scale parallel direct solver; and practical experience shows that speedups exceeding 15 are difficult to achieve regardless of the number of processors used in 2D and 3D applications; the memory complexity and time complexity of the direct solver for the two-dimensional finite difference problem are O(N) and O(N) respectively. 2 log2N) and O(N) 3 In three-dimensional space, the values ​​increase dramatically to O(N). 4 ) and O(N 6 ).

[0040] In step 4, another approach to frequency domain modeling is based on the iterative solver ISM. Compared to DSM, its main advantage lies in its smaller memory requirements, typically O(N) for 3D. 3 ).

[0041] In step 5, the hybrid frequency domain modeling method, HSM, offers a good trade-off between DSM and ISM in terms of memory requirements and multi-rhs simulation efficiency; HSM is based on the domain decomposition method and uses a hybrid direct / iterative solver.

[0042] In step 6, a suitable solver is selected based on computer resources and the number of sources and receivers in the seismic experiment.

[0043] The efficient forward modeling algorithm based on the Discrete Fourier Transform (TDM) viscoelastic wave equation in this invention obtains the frequency domain wave field based on time-based discrete Fourier summation, and then selects an appropriate solver according to computer resources and the number of sources / receivers in the seismic experiment.

[0044] When considering elasticity, anisotropy, and other extensions, DSM can efficiently implement two-dimensional frequency domain FWI. Three-dimensional frequency domain FWI requires more fine-grained distinctions. From previous numerical experiments, we can conclude that 3D DSM is still suitable for handling small models (i.e., if a large amount of resources needs to be considered, a suitable distributed memory platform consisting of a finite number of processors and a large number of shared memory processors can be used). Due to the low memory requirements of TDM, coarse-grained parallelization of sources can be seen on distributed memory platforms with a large number of processors, usually in the same order as the number of sources. If the number of processors significantly exceeds the number of sources, additional levels of parallelism can be easily achieved through classical domain decomposition. TDM also allows us to extract any number of frequency components via discrete Fourier transform, which would require additional computation if multiple frequencies need to be inverted simultaneously in frequency domain FWI. TDM should also be the most robust to the complexity of the medium and best suited for development towards elastic and anisotropic wave models. As Pratt and Sirgue stated in EAGE 2008... As noted at the FWI workshop, a key feature of time-division multiplexing (TDM) is its ability to flexibly select specific arrival points via time windows during modeling, while leaving the inversion in the frequency domain. The frequency domain remains the domain of choice for realizing attenuation effects of arbitrary complexity. In frequency-domain FWI, HSM adds additional computational overhead when inversions of multiple frequencies need to be considered simultaneously. TDM is also arguably the most robust approach to considering the complexity of the medium and is best suited for the development of elastic and anisotropic wave models. As Pratt and Sirgue pointed out at the 2008 EAGE FWI workshop, a key feature of time-division multiplexing is its ability to flexibly select specific arrival points via time windows during modeling, while leaving the inversion in the frequency domain. The frequency domain remains the domain of choice for realizing attenuation effects of arbitrary complexity.

[0045] To date, the HSM method has shown low efficiency, although it reduces the memory requirements of the DSM for large-scale problems. However, the conclusions achievable by iterative solver-based methods are highly dependent on the efficiency of the pre-tuning factor, which remains an active area of ​​research. Recent advancements in iterative solvers have yielded very promising results for the three-dimensional Helmholtz equations with optimal time complexity. Extending the three-dimensional elastic wave equations remains an ongoing problem.

[0046] This invention can accurately propagate in any heterogeneous layered medium without instability. In the three-dimensional case, it combines discrete Fourier transform to extract the frequency response and solves the three-dimensional viscoelastic wave equation based on frequency domain modeling using direct and iterative or hybrid solvers. Attached Figure Description

[0047] Figure 1This is a flowchart of a specific embodiment of the efficient forward modeling algorithm based on the discrete Fourier TDM viscoelastic wave equation of the present invention;

[0048] Figure 2 This is a schematic diagram illustrating the memory and time complexity of the efficient forward modeling algorithm TDM based on the discrete Fourier TDM viscoelastic wave equation of the present invention.

[0049] Figure 3 This is a schematic diagram illustrating the running time calculation of the efficient forward modeling algorithms DSM, HSM, and TDM based on the Discrete Fourier Transform (TDM) viscoelastic wave equation of the present invention.

[0050] Figure 4 This diagram illustrates the time consumed by the efficient forward modeling algorithm based on the Discrete Fourier TDM viscoelastic wave equation of the present invention as a function of the total number of DSM, HSM, and TDM sources. Detailed Implementation

[0051] To make the above and other objects, features, and advantages of the present invention more apparent and understandable, preferred embodiments are described below in detail with reference to the accompanying drawings. The descriptive information provided in the specific embodiments and applications is for illustrative purposes only and should not be construed as limiting the invention. The general principles defined in this invention may be applied in other embodiments without departing from the scope of the invention.

[0052] Full-waveform inversion FWI can be a good tool for seismic imaging. Frequency domain FWI has a significant advantage over time domain FWI because it introduces a natural multi-scale method from low frequency to high order and can manage and process compact data volumes. However, inversion is based on forward modeling.

[0053] Various forward modeling tools have been designed for 2D and 3D time-domain FWI and frequency-domain FWI. For 2D frequency-domain FWI, the forward modeling problem is naturally solved in the frequency domain using a direct solver. In the 3D case, the optimal strategy for forward modeling is not obvious. It can be approached from time-domain modeling (TDM) combined with discrete Fourier transform to extract the frequency response, to frequency-domain modeling based on direct and hybrid or iterative solvers.

[0054] like Figure 1 As shown, Figure 1This is a flowchart of the efficient forward modeling algorithm for the viscoelastic wave equation based on Discrete Fourier Transform (TDM) of the present invention. The flowchart includes: first, determining the initial model; then, performing time-domain and frequency-domain modeling on the initial model; then, employing a time-domain modeling (TDM) using an integrated scheme with an explicit format; next, calculating frequency-domain modeling (DSM) based on a direct solver; then calculating frequency-domain modeling (ISM) based on an iterative solver; and finally calculating frequency-domain modeling (HSM) based on a hybrid solver. Finally, these three methods are compared to determine the appropriate solver selection based on different computer resources and the number of sources / receivers in the seismic experiment.

[0055] The following are several specific embodiments of the application of the present invention.

[0056] Example 1

[0057] In a specific embodiment 1 of the present invention, the efficient forward modeling algorithm based on the discrete Fourier TDM viscoelastic wave equation includes the following steps:

[0058] In step 1, the initial model is determined.

[0059] In step 2, the initial model is modeled in both the time and frequency domains. Numerical methods transform partial differential operators into algebraic operations that can be represented by matrices. The wave equation attempts to estimate vectors, including pressure, solid particle velocity, or fluid / solid particle velocity, through an explicit system over time.

[0060]

[0061] Where x, y, and z are three-dimensional spatial coordinate parameters, and t is a time parameter.

[0062] The quality matrix M is a diagonal matrix;

[0063] The stiffness matrix A should be operated on back and forth in the spatial domain or the spectral domain to obtain the desired solution;

[0064] In controlled-source seismology, the source term S is usually a local point source, and the corresponding formula in the frequency domain is a generalization of the Helmholtz equation.

[0065]

[0066] Where x, y, and z are the three-dimensional spatial coordinate parameters, and w is the angular velocity.

[0067] When considering finite discretization, the impedance matrix B is complex and has a symmetric mode.

[0068] The dimension of the matrix is ​​the number of unknowns in the computational grid. The numerical bandwidth and the number of non-zero coefficient matrices depend on the template of the numerical discretization, which is embedded in the frequency domain using differential operators. Wave modeling reduces the difficulty of solving large sparse linear systems of equations and multiple right-hand sides (RHS), with each RHS corresponding to a source. A linear system is obtained.

[0069] In step 3, time modeling (TDM) with an explicit integrated scheme is employed. TDM with an explicit integrated scheme is performed; within the framework of multi-source simulation, the low memory requirements of TDM also allow for coarse-grained parallelism on the sources, which can be combined with domain decomposition parallelism if the number of processors is significantly greater than the number of sources.

[0070] In step 4, the frequency domain model (DSM) based on the direct solver is computed. Parallelism in the direct solver method DSM is achieved by using a large-scale parallel direct solver; this approach looks quite attractive for low frequencies, although it involves significant memory requirements. To date, 3D DSM has been limited to acoustic formulations.

[0071] In step 5, frequency domain modeling (ISM) based on iterative solver is computed; another method for frequency domain modeling is based on iterative solver ISM, whose main advantage over DSM is its smaller memory requirement, typically O(N) for 3D. 3 The drawback is that the impedance matrix is ​​indeterminate (the real eigenvalues ​​of B with varying signs), and therefore ill-conditioned.

[0072] Designing effective preconditions for equation (2) is currently a hot research topic. A recently developed multilevel Krylov method is employed, which generates numerous iterations that are almost independent of the problem size or frequency.

[0073] In step 6, frequency domain modeling (HSM) based on a hybrid solver is computed. HSM, a hybrid frequency domain modeling method, offers a good trade-off between DSM and ISM in terms of memory requirements and multi-RHS simulation efficiency. The HSM is based on a domain decomposition approach using a hybrid direct / iterative solver. The management idea is to divide the unknowns into two subsets: unknowns related to the subdomain interface and unknowns related to the subdomain interior. One processor is allocated to each subdomain.

[0074] In step 7, the results of the three solvers are compared. The appropriate solver selection is determined based on different computer resources and the number of sources / receivers in the seismic experiment.

[0075] Example 2

[0076] In a specific embodiment 2 of the present invention, the efficient forward modeling algorithm based on the discrete Fourier TDM viscoelastic wave equation of the present invention includes the following steps:

[0077] In step 1, the initial model is determined.

[0078] In step 2, the initial model is modeled in both the time and frequency domains. The dimension of the matrix is ​​the number of unknowns in the computational grid. The numerical bandwidth and the number of non-zero coefficient matrices depend on the template of the numerical discretization, which is embedded in the frequency domain using differential operators. Wave modeling reduces the difficulty of solving large sparse linear equation systems and multiple right-hand sides (RHS), with each RHS corresponding to a source. A linear system is obtained.

[0079] For two-dimensional frequency domain FWI, the forward modeling problem is realized in the frequency domain by a direct solver. In the three-dimensional case, the optimal strategy for forward modeling is not obvious. It can be improved by using time-domain modeling (TDM) combined with discrete Fourier transform to extract the frequency response, or by using frequency domain modeling based on direct and hybrid or iterative solvers.

[0080] Considering the three-dimensional viscous acoustic wave equation in the frequency domain

[0081]

[0082] Here, density is represented by ρ(x), bulk modulus by κ(x), and angular frequency by ω. The monochromatic pressure wave field and source are represented by p(x, ω) and s(x, ω), respectively. In the expression for bulk modulus, the inherent attenuation can be easily realized in the frequency domain using complex-valued wave velocities.

[0083] Equation (1) can be reproduced in matrix form as follows:

[0084] Ap = s (2)

[0085] The complex impedance matrix A depends on the angular frequency ω and the parameters κ and ρ.

[0086] The discretization of the wave equation was performed using the compact finite difference template of Opera et al. (2007), which was originally designed for direct methods but is suitable for substructured methods because its local support allows minimizing the interface size (width) between adjacent subdomains, unlike higher-order finite difference methods.

[0087] A discretization criterion of four grid points per wavelength is used below in this study. For the three-dimensional problem, this template involves 27 coefficients spanning two grid intervals across three Cartesian directions, resulting in a numerical band accuracy of O(N) for A. 2 The absorbing boundary conditions are implemented by a fully matched layer.

[0088] In step 3, time modeling (TDM) with an explicit format is employed. Partial differential operators are transformed into algebraic operations using existing iterative numerical methods, which can be represented by matrices. The wave equation attempts to estimate the vector f over time through the explicit system, where f is pressure, solid particle velocity, or fluid / solid particle velocity.

[0089] The matrix expression is:

[0090]

[0091] Where xyz are the three-dimensional spatial coordinate parameters, and t is the time parameter;

[0092] The quality matrix M is a diagonal matrix;

[0093] The stiffness matrix A should be operated on back and forth in the spatial domain or the spectral domain to obtain the desired solution;

[0094] In controlled-source seismology, the source term S is usually a local point source, and the corresponding formula in the frequency domain is a generalization of the Helmholtz equation.

[0095]

[0096] Where xyz are the three-dimensional spatial coordinate parameters, and w is the angular velocity;

[0097] When considering finite discretization, the impedance matrix B is complex and has a symmetric mode.

[0098] TDM is typically performed using an explicit time planning algorithm; at each time step, the solution for each spatial grid point is estimated from the solution of the previous time step.

[0099] For frequency-domain FWI, core memory storage of the full time series is useless because the frequency-domain wavefield is obtained through discrete Fourier summation over time. At each time step, the solution for each spatial grid point is estimated from the solution of the previous time step. By using standard domain decomposition methods to divide the computational domain into finite-dimensional subdomains, time-domain algorithms can be effectively parallelized. The efficiency of these algorithms is typically close to 1.0. If parallelism is only achieved on the source, then... Figure 2 The memory and time complexity of TDM are given.

[0100] By employing standard domain decomposition methods, the computational domain is divided into finite-dimensional subdomains, which effectively parallelizes time-domain algorithms. The efficiency of these algorithms is generally close to 1.0. In the framework of multi-source simulation, the low memory requirements of TDM also allow for coarse-grained parallelism on the sources. If the number of processors is significantly greater than the number of sources, this parallelism can be combined with domain decomposition parallelism. In the following discussion, the dimension of the three-dimensional N3 computational grid is denoted as n. If parallelism is only implemented on the sources, Figure 1 The memory and time complexity of TDM are given. Real-world 3D surveys require a large amount of memory to store N distributed across processors. rhs Wave field, O(N) 3 N rhs )=O(N 5 ), where N rhs The number of sources is on one side of the grid. With sufficiently fine time discretization methods, the TDM algorithm can accurately achieve propagation in any unstable heterogeneous layer medium. Despite a significant increase in memory requirements, extensions to resilience, anisotropy, and attenuation modeling are possible.

[0101] In step 4, the frequency domain modeling (DSM) based on the direct solver is computed. The finite number of frequencies required for the frequency domain FWI can efficiently model a large number of sources, simply by performing LU decomposition on the impedance matrix B at each frequency. Parallelism in the direct solver method DSM is achieved using a large-scale parallel direct solver. Furthermore, our practical experience shows that speedups exceeding 15 are difficult to achieve regardless of the number of processors used in 2D and 3D applications. The memory complexity and time complexity of the direct solver for the two-dimensional finite difference problem are O(N^2). 2 log2N) and O(N) 3 In three-dimensional space, the values ​​increase dramatically to O(N). 4 ) and O(N 6 The parentheses following the O here contain a function that indicates the relationship between the time / space consumption of an algorithm and the amount of data growth. N represents the amount of input data. The complexity order is O(1). <O(log2N)<O(N)<O(N 2 ).

[0102] In step 5, frequency domain modeling (ISM) based on the iterative solver is computed. Its main advantage over DSM is its smaller memory requirement, typically O(N) for 3D. 3 The drawback is that the impedance matrix is ​​indeterminate (the real eigenvalues ​​of B with varying signs), and therefore ill-conditioned. A recent, novel multi-level Krylov method is employed, which generates numerous iterations that are almost independent of the problem size or frequency.

[0103] In this ideal scenario, for 3D problems, the time complexity of the iterative solver is O(N). 3 The previous version of the method showed a time complexity of O(N), while the previous version showed a time complexity of O(N). 4 (Note that the number of iterations increases linearly with the increase of N.)

[0104] For multi-source problems and coarse-grained parallelism on the sources, the time complexity of the iterative solver is O(N). 4 In theory, iterative solvers should provide the most efficient numerical format. However, attenuation effects can easily be introduced.

[0105] In step 6, frequency domain modeling (HSM) based on a hybrid solver is computed. First, a direct solver is used on each processor to perform sequential logical unit decomposition of the local matrices assembled on each subdomain. Second, an iterative solver such as GMRES is used to solve a simplified system, the so-called Schur's complement system, which is better positioned than the full system conditions handled by the ISM. The solution is the interface unknown, and the RHS is derived from the previous decomposition steps. Once the interface unknowns are determined, the internal unknowns can be efficiently computed by performing forward / backward substitutions on each processor. The parameter ε (residual divided by the norm of the RHS) represents the stopping criterion for the GMRES iteration.

[0106] In step 7, the results of the three solvers are compared. The appropriate solver selection is determined based on different computer resources and the number of sources / receivers in the seismic experiment.

[0107] N p = Total number of processors.

[0108] N rhs = Total number of sources.

[0109] For DSM:

[0110] N DSM = The number of processors dedicated to a single LU decomposition.

[0111] T LU =Time consumed by LU decomposition

[0112] T rhs = The time taken to complete an RHS solution.

[0113] For HSM:

[0114] N HSM = The number of processors in a domain decomposition,

[0115] T LU+M =Running time of RHS standalone task

[0116] T GMRES =GMRES runtime

[0117] For TDM:

[0118] T seq = 1 - rhs sequential simulation runtime.

[0119] In our implementation, we have Figure 4 The TDM was observed to have a slower slope than the HSM, which makes the former superior, although any improvement in iterative preprocessing could challenge this conclusion. Due to the efficiency of the solution step, the DSM performs better with a large number of sources (over 500 sources).

[0120] Example 3

[0121] In a specific embodiment 3 of the present invention, the invention is applied. Figure 2 This is a schematic diagram illustrating the memory and time complexity of the efficient forward modeling algorithm TDM based on the discrete Fourier TDM viscoelastic wave equation of the present invention.

[0122] Figure 3 This is a schematic diagram illustrating the running time calculation of the efficient forward modeling algorithms DSM, HSM, and TDM based on the Discrete Fourier Transform (TDM) viscoelastic wave equation of the present invention.

[0123] Figure 4 This diagram illustrates the time consumed by the efficient forward modeling algorithm based on the Discrete Fourier Transform (DFT) TDM viscoelastic wave equation of this invention as a function of the total number of DSM, HSM, and TDM sources. The vertical axis represents the estimated time (in hours), and the horizontal axis represents the lens number.

[0124] Elapsed time is a function of the total number of sources for DSM, HSM, and TDM. This curve calculates the number of three available processors: 192,500 and 2000 (represented in boxes). 192 processors are grouped for parallel logic unit decomposition of the DSM and domain decomposition of the HSM. The initial LU factorization is negligible for the HSM. For the TDM, parallelization is performed via lenses.

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

[0126] Except for the technical features described in the specification, all other technologies are known to those skilled in the art.

Claims

1. An efficient forward modeling algorithm based on the Discrete Fourier Transform (TDM) viscoelastic wave equation, characterized in that, This efficient forward modeling algorithm based on the Discrete Fourier Transform (TDM) viscoelastic wave equation includes: Step 1: Perform frequency domain modeling on the initial model; Step 2: Employ time-mapping (TDM) using an integrated scheme with an explicit format. Specifically, the partial differential operators in the viscoelastic wave equation are transformed into algebraic operations, and the vector f is estimated using an explicit time-progression scheme. The vector f represents pressure, solid particle velocity, or fluid / solid particle velocity. Its time-domain matrix expression is as follows: （1） Where xyz are the three-dimensional spatial coordinate parameters, t is the time parameter, M is the mass matrix, A is the stiffness matrix, and S is the source term; The source term S corresponds to a generalization of the Helmholtz equation in the frequency domain, and its frequency domain matrix expression is as follows: （2） Where w is the angular frequency and B is the impedance matrix; when finite discretization is performed, the impedance matrix B is a complex-valued matrix and has a symmetric mode. Step 3: Calculate the frequency domain modeling DSM based on the direct solver; Step 4: Calculate the frequency domain modeling ISM based on the iterative solver; Step 5: Calculate the frequency domain modeling HSM based on the hybrid solver; Step 6: Compare these three methods.

2. The efficient forward modeling algorithm based on the Discrete Fourier Transform (TDM) viscoelastic wave equation according to claim 1, characterized in that, In step 1, for two-dimensional frequency domain FWI, the forward modeling problem is implemented in the frequency domain by a direct solver. In the three-dimensional case, the optimal strategy for forward modeling is not obvious. It can be improved from time domain modeling (TDM) combined with discrete Fourier transform to extract the frequency response, to frequency domain modeling based on direct and hybrid or iterative solvers.

3. The efficient forward modeling algorithm based on the discrete Fourier TDM viscoelastic wave equation as described in claim 2, characterized in that, In step 1, the three-dimensional viscous acoustic wave equation is considered in the frequency domain: (1) Density is represented by ρ(x), bulk modulus by κ(x), and angular frequency by ω; the monochromatic pressure wave field and source are represented by p(x, ω) and s(x, ω), respectively; in the expression for bulk modulus, the inherent attenuation can be easily achieved in the frequency domain using complex wave velocity; Equation (1) can be reproduced in matrix form as follows: (2) The complex impedance matrix A depends on the angular frequency ω and the parameters κ and ρ, p is the spatially discretized monochromatic pressure wave field vector, and s is the spatially discretized source term vector.

4. The efficient forward modeling algorithm based on the discrete Fourier TDM viscoelastic wave equation according to claim 1, characterized in that, In step 2, TDM is typically performed using an explicit time planning algorithm; at each time step, the solution for each spatial grid point is estimated from the solution of the previous time step.

5. The efficient forward modeling algorithm based on the Discrete Fourier Transform (TDM) viscoelastic wave equation according to claim 1, characterized in that, In step 6, a suitable solver is selected based on computer resources and the number of sources and receivers in the seismic experiment.

Images