CN117150703A - Efficient forward algorithm based on discrete Fourier TDM viscoelastic wave equation - Google Patents
Efficient forward algorithm based on discrete Fourier TDM viscoelastic wave equation Download PDFInfo
- Publication number
- CN117150703A CN117150703A CN202210566067.6A CN202210566067A CN117150703A CN 117150703 A CN117150703 A CN 117150703A CN 202210566067 A CN202210566067 A CN 202210566067A CN 117150703 A CN117150703 A CN 117150703A
- Authority
- CN
- China
- Prior art keywords
- tdm
- frequency domain
- discrete fourier
- modeling
- wave equation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 51
- 238000000034 method Methods 0.000 claims abstract description 43
- 230000010354 integration Effects 0.000 claims abstract description 10
- 230000004044 response Effects 0.000 claims abstract description 7
- 239000011159 matrix material Substances 0.000 claims description 38
- 238000000354 decomposition reaction Methods 0.000 claims description 20
- 238000013459 approach Methods 0.000 claims description 11
- 238000004364 calculation method Methods 0.000 claims description 10
- 239000002245 particle Substances 0.000 claims description 8
- 238000004088 simulation Methods 0.000 claims description 8
- 239000007787 solid Substances 0.000 claims description 8
- 238000002474 experimental method Methods 0.000 claims description 7
- 230000014509 gene expression Effects 0.000 claims description 6
- 239000013598 vector Substances 0.000 claims description 6
- 239000012530 fluid Substances 0.000 claims description 4
- 230000036961 partial effect Effects 0.000 claims description 4
- 230000003595 spectral effect Effects 0.000 claims description 4
- 230000001133 acceleration Effects 0.000 claims description 3
- 230000005055 memory storage Effects 0.000 claims description 3
- 208000030514 Leukocyte adhesion deficiency type II Diseases 0.000 claims 1
- 239000000284 extract Substances 0.000 abstract description 2
- 230000008901 benefit Effects 0.000 description 6
- 238000010586 diagram Methods 0.000 description 6
- 230000006870 function Effects 0.000 description 5
- 230000008569 process Effects 0.000 description 5
- 238000010521 absorption reaction Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 230000018109 developmental process Effects 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 230000008859 change Effects 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000035945 sensitivity Effects 0.000 description 2
- 238000012512 characterization method Methods 0.000 description 1
- 238000012824 chemical production Methods 0.000 description 1
- 238000004587 chromatography analysis Methods 0.000 description 1
- 230000004154 complement system Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000009472 formulation Methods 0.000 description 1
- 230000000670 limiting effect Effects 0.000 description 1
- 238000007726 management method Methods 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- -1 pressure Substances 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000001902 propagating effect Effects 0.000 description 1
- 239000002994 raw material Substances 0.000 description 1
- 230000002829 reductive effect Effects 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 230000000630 rising effect Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000003860 storage Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
- G06F17/141—Discrete Fourier transforms
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Abstract
The invention provides a TDM viscoelastic wave equation high-efficiency positive algorithm based on discrete Fourier, comprising the following steps: step 1: carrying out frequency domain modeling on the initial model; step 2: time modeling (TDM) employing an integration scheme with explicit format; step 3: computing a direct solver-based frequency domain modeling (DSM); step 4: computing an iterative solver-based frequency domain modeling (ISM); step 5: computing a hybrid solver-based frequency domain modeling (HSM); step 6: a comparison was made for these three methods. The efficient forward algorithm based on the discrete Fourier TDM viscoelastic wave equation is a forward algorithm based on an explicit integration scheme for time modeling, can accurately propagate in any heterogeneous layer medium without instability, extracts frequency response by combining discrete Fourier transform under the three-dimensional condition, and solves the three-dimensional viscoelastic wave equation based on frequency domain modeling of a direct and iterative or hybrid solver.
Description
Technical Field
The invention relates to the technical field of modeling of seismic wave imaging, in particular to a high-efficiency positive algorithm based on a discrete Fourier TDM viscoelastic wave equation.
Background
The exploration of modern oil and gas resources originates in the middle of the 19 th century, the second industrial revolution. Oil and gas resources have become one of the most important energy sources for human society for over a century. The oil gas resource is not only a main source of daily traffic and transportation in the human society and a main raw material for chemical production, but also a national important strategic resource. China is a large country of consumption of oil and gas resources, and the demand for the oil and gas resources is rising year by year. The oil gas used by us is conventional oil gas and unconventional oil gas. The Chinese unconventional oil gas resources are rich, the quality of the oil gas stock is reduced, and the technical requirements are increased. In practical oil and gas exploration, common means are seismic exploration, acoustic logging and the like. The qualitative and quantitative study of reservoir parameters and structures is mainly performed by using wave information propagating in the reservoir.
Classical hyperbolic acoustic and elastic wave equations cannot accurately describe geologic models of such reservoirs, and traditional numerical simulation algorithms are also difficult to process complex coupling characterization equations. In order to meet the requirements of unconventional oil and gas reservoir exploration and development, the construction work of a new stable and efficient algorithm is very important.
In application number: in CN202110609064.1, the method, system, device and medium for seismic inversion based on wave equation are related, and the method comprises the following steps: acquiring and determining calculation parameters of a single-frequency sensitivity kernel function according to the seismic observation parameters and the initial speed model parameters; generating a velocity model of the random boundary; forward wave equation modeling is carried out according to the calculation parameters and the speed model, so that a synthetic seismic record and a frequency domain incident wave field are obtained; determining a travel time accompanying seismic source according to the synthetic seismic records, and determining a single-frequency gradient of a travel time target functional according to the travel time accompanying seismic source; extracting a low wave number part from the single-frequency gradient, and carrying out seismic velocity inversion according to the low wave number; compared with wave equation travel time chromatography technology based on wave field reconstruction frame, the method has remarkable advantages in the aspects of calculation efficiency, memory occupation and algorithm complexity, and can be widely applied to the technical field of seismic velocity modeling.
In application number: in the chinese patent application CN202110177528.6, a method, an apparatus, a device and a storage medium for forward modeling of viscoelastic waves in a frequency domain are related, where the method includes: acquiring a frequency wave number domain one-dimensional viscoelasticity wave equation, and discretizing the wave equation into a finite element control equation after converting the wave equation into an equivalent integral weak form; loading a three-dimensional point force source and an absorption boundary condition of a target work area to a finite element control equation; solving a finite element control equation loaded with a three-dimensional point force source and an absorption boundary condition to obtain a frequency wave number domain wave field of a target work area; determining a wave number set corresponding to the wave field of the frequency wave number domain according to the self-adaptive wave number sampling strategy; and performing unequal interval two-dimensional Fourier inverse transformation on the frequency wave number domain wave field corresponding to each wave number in the wave number set when the wave number is used as input to obtain the corresponding frequency domain three-dimensional elastic wave field. The invention can improve the accuracy of forward simulation of the frequency domain viscoelastic wave.
In application number: in CN201911315597.8, a method for inversion of multi-wave combined prestack waveforms is disclosed. The pre-stack seismic inversion method mainly comprises two main categories, namely a positive operator based on ray tracing and a positive operator based on wave equation. The former is the most widely used pre-stack AVO inversion, and the latter is called pre-stack waveform inversion. Different from the pre-stack AVO inversion, the pre-stack waveform inversion method can better apply complex wave field response in input data, reduce the processing difficulty of the input data and improve the sensitivity of the seismic record to target inversion parameters. At present, nonlinear inversion strategies are mostly adopted for pre-stack waveform inversion, and the calculation efficiency is too low. The invention provides a pre-stack inversion method of a linear strategy to reduce calculation time. Meanwhile, in order to improve inversion precision and stability, the invention adopts a joint inversion strategy of PP and PS multi-wave data, and the method can effectively improve the inversion precision of three parameters, especially for transverse wave speed and density estimation. Logging model testing verifies the effectiveness of the method.
The prior art is greatly different from the invention, the technical problem which is needed to be solved is not solved, and a novel efficient forward algorithm based on the discrete Fourier TDM viscoelastic wave equation is invented.
Disclosure of Invention
The invention aims to provide a high-efficiency forward algorithm based on a discrete Fourier TDM viscoelastic wave equation, which can accurately propagate in any heterogeneous layer medium without instability.
The aim of the invention can be achieved by the following technical measures: the efficient forward algorithm based on the discrete Fourier TDM viscoelastic wave equation comprises the following steps:
step 1, carrying out frequency domain modeling on an initial model;
step 2: time modeling TDM with an integration scheme having an explicit format;
step 3: calculating a direct solver-based frequency domain modeling, DSM;
step 4: calculating a frequency domain modeling ISM based on an iteration solver;
step 5: calculating a frequency domain modeling HSM based on a hybrid solver;
step 6: a comparison was made for these three methods.
The aim of the invention can be achieved by the following technical measures:
in step 1, for two-dimensional frequency domain FWI, the forward problem is implemented in the frequency domain by a direct solver, and in the three-dimensional case, the optimal strategy for forward modeling is less obvious, and the frequency response can be extracted from the time domain modeling TDM in combination with the discrete fourier transform to the frequency domain modeling based on the direct and hybrid or iterative solver.
In step 1, consider the three-dimensional viscous acoustic wave equation in the frequency domain:
wherein density is represented by ρ (x), bulk modulus is represented by κ (x), and angular frequency is represented by ω; the monochromatic pressure wavefield and source are denoted by p (x, ω) and s (x, ω), respectively; in bulk modulus expressions, inherent attenuation can be easily achieved in the frequency domain using complex-valued wave velocities;
formula (1) can be reproduced in matrix form
Ap＝s (2)
Wherein the complex-valued impedance matrix a depends on the angular frequency ω and the median parameters κ and ρ.
In step 2, the partial differential operator is converted to an algebraic operation using existing iterative numerical methods, which can be represented by a matrix, the wave equation attempting to estimate the vector f over time by an explicit system, the vector f being the pressure, the velocity of the solid particles, or the velocity of the fluid/solid particles.
In step 2, the matrix expression is:
where xyz is a three-dimensional spatial coordinate parameter and t is a time parameter;
the quality matrix M is a diagonal matrix;
the stiffness matrix a should operate on the desired solution back and forth in the spatial domain or in the spectral domain;
the source term S is usually a local point source in vibroseis seismology, and the corresponding formula in the frequency domain is popularization by adopting a Helmholtz equation;
where xyz is a three-dimensional spatial coordinate parameter and w is angular velocity;
wherein the impedance matrix B is complex valued and has a symmetrical pattern when considering the finite discretization.
In step 2, TDM is typically performed according to an explicit time planning algorithm; in each time step, the solution for each spatial grid point is estimated from the solution for the previous time step.
In step 2, for frequency domain FWI, the core memory storage of the full time sequence is useless, since the frequency domain wavefield is obtained by discrete fourier summation; the calculation domain is divided into subdomains with limited dimensions by adopting a standard domain decomposition method, so that a time domain algorithm can be effectively parallelized; in the framework of multi-source emulation, the low memory requirements of TDM also allow coarse-grained parallelism to the sources, which can be combined with domain decomposition parallelism if the number of processors is significantly greater than the number of sources.
In step 2, the dimension of the three-dimensional N3 computational grid is denoted as N; if only parallelism on the source is realized, the memory and time complexity of TDM are given; real 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 ) Wherein N is rhs The number of sources is on one side of the grid; the TDM algorithm can accurately achieve propagation in any heterogeneous layer medium without instability as long as a sufficiently fine time discretization method is used; although memory requirements increase significantly, extensions to elasticity, anisotropy, and decay modeling are possible.
In step 2, the dimensions of the matrix are the number of unknowns in the computational grid, the numerical bandwidth and the number of non-zero coefficient matrices are embedded in the frequency domain by using differential operators depending on templates of numerical discretization, wave modeling reduces solving a large sparse linear equation set and multiple right terms RHS, each RHS corresponds to one seismic source, and a linear system is obtained.
In step 2, time domain modeling TDM of the explicit integration scheme is performed; in the framework of multi-source simulation, the low memory requirements of TDM also allow coarse-grained parallelism on the sources, which can be combined if the number of processors is significantly greater than the number of sources.
In step 3, a large number of sources can be effectively modeled as long as LU decomposition is performed on the impedance matrix B of each frequency, the finite frequency required for the frequency domain FWI; parallelism in the direct solution method DSM is achieved by using a massive parallel direct solver; moreover, practical experience shows that no matter how many processors are used by 2D and 3D applications, acceleration of over 15 is difficult to achieve; the memory complexity and time complexity of the two-dimensional finite difference problem direct solver are O (N) 2 log 2 N) and O (N) 3 ) In three dimensions, the expansion is greatly exaggerated to O (N 4 ) And O (N) 6 )。
In step 4, another approach to frequency domain modeling is based on an iterative solver ISM, which has the major advantage of smaller memory requirements compared to DSM, typically O (N for 3D 3 )。
In step 5, the hybrid frequency domain modeling method, HSM, provides a good compromise between DSM and ISM in terms of memory requirements and multiple rhs simulation efficiency; the HSM is based on a region decomposition approach, using a hybrid direct/iterative solver.
In step 6, a suitable solver is selected based on the computer resources and the number of sources and receivers for the seismic experiment.
The efficient forward algorithm based on the discrete Fourier TDM viscoelastic wave equation is based on the time discrete Fourier summation to obtain a frequency domain wave field, and then a proper solver is selected according to the computer resources and the number of sources/receivers of the seismic experiment.
The use of DSM can effectively achieve two-dimensional frequency domain FWI when considering elasticity, anisotropy, and other expansions. Three-dimensional frequency domain FWI requires finer differentiation. From previous numerical experiments we can conclude that 3D DSM is still a small model to handle (i.e. if a large number of resources need to be considered, if a suitable distributed memory platform consists of a limited number of processors and a large number of shared memory handlers, this platform can be used, because of the low demands of TDM on memory, coarse-grained parallelization of sources can be seen on a distributed memory platform with a large number of processors, typically 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 by classical domain decomposition TDM also allows us to extract any number of frequency components by discrete fourier transform, if simultaneous inversion of multiple frequencies needs to be considered in the frequency domain FWI, additional computation is needed TDM should also be most robust to take the complexity of the medium into account, as noted by Pratt and Sirgue at the EAGE FWI seminar in 2008, a good feature of time division multiplexing is that specific points of arrival are flexibly selected through a time window during modeling while leaving inversion in the frequency domain, which is still the selected domain for achieving any complexity attenuation effect. While the inversion is retained in the frequency domain. The frequency domain is still the chosen domain for achieving any complexity attenuation effect.
So far, the HSM approach has been shown to be less efficient, although it reduces the memory requirements of DSM for large scale problems. However, the conclusions that can be drawn by iterative solver-based methods are highly dependent on the efficiency of the preconditioning factor, which remains an active area of research. Recent developments in iterative solvers provide very promising results for three-dimensional helmholtz equations, with optimal time complexity. Expansion of the three-dimensional elastic wave equation remains a continuing problem.
The invention can accurately propagate in any heterogeneous layer medium without instability, extracts frequency response by combining discrete Fourier transform under the three-dimensional condition, and solves the three-dimensional viscoelastic wave equation based on frequency domain modeling of a direct and iterative or hybrid solver.
Drawings
FIG. 1 is a flow chart of an embodiment of the efficient forward algorithm of the present invention based on a discrete Fourier TDM viscoelastic wave equation;
FIG. 2 is a schematic diagram of the memory and time complexity of the efficient forward algorithm TDM based on the discrete Fourier TDM viscoelastic wave equation of the present invention;
FIG. 3 is a schematic diagram of the run-time calculation of the discrete Fourier TDM viscoelastic wave equation based efficient forward algorithm DSM, HSM, TDM of the present invention;
FIG. 4 is a schematic diagram of the time spent by the discrete Fourier TDM based viscoelastic wave equation high efficiency forward algorithm of the present invention as a function of the total number of DSM, HSM and TDM sources.
Detailed Description
The foregoing and other objects, features and advantages of the invention will be apparent from the following more particular description of preferred embodiments, as illustrated in the accompanying drawings. The description of the embodiments and applications is provided only for illustrating the invention and should not be taken as limiting the invention. The general principles defined herein may be applied to other embodiments without departing from the scope of the invention.
Full waveform inversion FWI can be a good tool for seismic imaging, where frequency domain FWI has significant advantages over time domain FWI because it introduces a natural multi-scale approach from low frequency to high order and is able to manage and process compact data volumes. But the inversion is based on forward modeling.
Various forward modeling tools were designed for two-dimensional and three-dimensional time-domain FWI and frequency-domain FWI. For two-dimensional frequency domain FWI, the forward problem is naturally implemented in the frequency domain by a direct solver. In three dimensions, the optimal strategy for forward modeling is less obvious, and the frequency response can be extracted from Time Domain Modeling (TDM) in combination with discrete fourier transforms to frequency domain modeling based on direct and hybrid or iterative solvers.
As shown in fig. 1, fig. 1 is a flowchart of the efficient forward algorithm of the present invention based on the discrete fourier TDM viscoelastic wave equation. The flow chart based on the discrete Fourier TDM viscoelastic wave equation efficient forward algorithm comprises the following steps: firstly, determining an initial model; then modeling the time domain and the frequency domain of the initial model; then time modeling (TDM) with an integration scheme in explicit format is employed; then calculating a direct solver-based frequency domain modeling (DSM); and computing an iterative solver-based frequency domain modeling (ISM); and computing a hybrid solver-based frequency domain modeling (HSM); finally, the three methods are compared. The number of sources/receivers based on different computer resources and seismic experiments was derived, with each appropriate solver selection.
The following are several embodiments of the invention
Example 1
In one embodiment 1 to which the present invention is applied, the discrete Fourier TDM based efficient forward algorithm for viscoelastic wave equation includes the steps of:
in step 1, an initial model is determined.
In step 2, time-domain and frequency-domain modeling is performed on the initial model. Numerical methods convert partial differential operators into algebraic operations, which can be represented by matrices, wave equations attempt to estimate vectors, including pressure, solid particle velocity, or fluid/solid particle velocity, through an explicit system of times.
Wherein x, y and z are three-dimensional space coordinate parameters, and t is time parameter
The quality matrix M is a diagonal matrix;
the stiffness matrix a should operate on the desired solution back and forth in the spatial domain or in the spectral domain;
the source term S is typically a local point source in vibroseis seismology, and the corresponding formula in the frequency domain is a generalization using the helmholtz equation.
Where xyz is a three-dimensional space coordinate parameter and w is angular velocity
When finite discretization is considered, the impedance matrix B is complex valued and has a symmetrical pattern.
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 number discretized templates embedded in the frequency domain using differential operators, wave modeling reduces resolution of large sparse linear equations and multiple right terms (RHS), one for each source. A linear system is obtained.
In step 3, time modeling (TDM) of an integration scheme with explicit format is employed. Performing time domain modeling TDM of an explicit integration scheme; in the framework of multi-source simulation, the low memory requirements of TDM also allow coarse-grained parallelism on the sources, which can be combined if the number of processors is significantly greater than the number of sources.
In step 4, a direct solver based frequency domain modeling (DSM) is calculated. Parallelism in the direct solution method DSM is achieved by using a massive parallel direct solver; for low frequencies, this approach appears to be quite attractive, despite the large number of memory requests, so far 3D DSM has been limited to acoustic formulations.
In step 5, computing an iterative solver-based frequency domain modeling (ISM); another approach to frequency domain modeling is based on an iterative solver ISM, which has the major advantage of smaller memory requirements compared to DSM, typically O (N for 3D 3 ). The disadvantage is that the impedance matrix is indefinite (real eigenvalue of the B change sign) and therefore ill-conditioned.
The design of effective pretreatment conditions for equation (2) is a hotspot of current research. With the recent development of a multi-layer Krylov process, this method produces many iterations that are almost independent of the problem size or frequency.
In step 6, a hybrid solver based frequency domain modeling (HSM) is calculated. A hybrid frequency domain modeling approach, HSM, provides a good compromise between DSM and ISM in terms of memory requirements and multiple rhs simulation efficiency. The HSM is based on a region decomposition approach, using a hybrid direct/iterative solver. The management idea is to divide the unknowns into two subsets, an unknown related to the subfield interfaces and an unknown related to the inside of the subfields. One processor is assigned to each sub-domain.
In step 7, the results of the three solvers are compared. The number of sources/receivers based on different computer resources and seismic experiments was derived, with each appropriate solver selection.
Example 2
In a specific embodiment 2 to which the present invention is applied, the discrete fourier TDM viscoelastic wave equation-based efficient positive algorithm of the present invention includes the steps of:
in step 1, an initial model is determined.
In step 2, time-domain and frequency-domain modeling is performed on the initial model. 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 number discretized templates embedded in the frequency domain using differential operators, wave modeling reduces resolution of large sparse linear equations and multiple right terms (RHS), one for each source. A linear system is obtained.
For two-dimensional frequency domain FWI, the forward problem is implemented in the frequency domain by a direct solver, and in the three-dimensional case, the optimal strategy for forward modeling is less obvious, and the frequency response can be extracted from the time domain modeling TDM in combination with the discrete fourier transform to the frequency domain modeling based on the direct and hybrid or iterative solver.
Considering three-dimensional viscous acoustic wave equation in frequency domain
Where density is denoted by ρ (x), bulk modulus by κ (x), and angular frequency by ω. The monochromatic pressure wavefield and source are denoted by p (x, ω) and s (x, ω), respectively. In bulk modulus expressions, inherent attenuation in the frequency domain can be easily achieved with complex-valued wave velocities.
Formula (1) can be reproduced in matrix form
Ap＝s (2)
Wherein the complex-valued impedance matrix a depends on the angular frequency ω and the median parameters κ and ρ.
Discretization of wave equations was performed with a compact finite difference template of Operto et al (2007), which was originally designed for a straightforward method, suitable for the substructure method, because its local support allows minimizing the interface size (width) between adjacent subfields, unlike the higher-order finite difference method.
A discrete criterion of four grid points per wavelength is used below the study. For three-dimensional problems, the template involves two grid intervals of 27 coefficients spanning three Cartesian directions, which results in the accuracy of the value band of A being O (N 2 ) The absorption boundary conditions are realized by completely matched layers
In step 3, time modeling (TDM) of an integration scheme with explicit format is employed. The partial differential operator is converted to an algebraic operation using existing iterative numerical methods, which can be represented by a matrix, and the wave equation attempts to estimate the vector f, over time, as pressure, solid particle velocity, or fluid/solid particle velocity by an explicit system.
The matrix expression is:
where xyz is a three-dimensional spatial coordinate parameter and t is a time parameter;
the quality matrix M is a diagonal matrix;
the stiffness matrix a should operate on the desired solution back and forth in the spatial domain or in the spectral domain;
the source term S is usually a local point source in vibroseis seismology, and the corresponding formula in the frequency domain is popularization by adopting a Helmholtz equation;
where xyz is a three-dimensional spatial coordinate parameter and w is angular velocity;
wherein the impedance matrix B is complex valued and has a symmetrical pattern when considering the finite discretization.
TDM is typically performed in accordance with an explicit time planning algorithm; in each time step, the solution for each spatial grid point is estimated from the solution for the previous time step.
For frequency domain FWI, the core memory storage of the full time sequence is useless, since the frequency domain wavefield is obtained by discrete fourier summation over time. At each time step, the solution for each spatial grid point is estimated from the solution for the previous time step. The computational domain is divided into subdomains of finite dimensions by a standard domain decomposition method, so that the time domain algorithm can be effectively parallelized. The efficiency of these algorithms is typically close to 1.0. If only parallelism on the source is implemented, then figure two gives the memory and time complexity of TDM.
The calculation domain is divided into subdomains with limited dimensions by adopting a standard domain decomposition method, so that the time domain algorithm can be effectively parallelized. The efficiency of these algorithms is typically close to 1.0. In the framework of multi-source emulation, the low memory requirements of TDM also allow coarse-grained parallelism to the sources, which can be combined with domain decomposition parallelism if the number of processors is significantly greater than the number of sources. In the discussion that follows, the dimension of the three-dimensional N3 computational grid is denoted as N. If only parallelism on the source is implemented, FIG. 1 gives the memory and time complexity of TDM. Real 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 ) Wherein N is rhs The number of sources is on one side of the grid. The TDM algorithm can accurately achieve propagation in any heterogeneous layer medium without instability as long as a sufficiently fine time discretization method is used. Although memory requests increase significantly, elasticity, anisotropy, andan extension of the attenuation modeling is possible.
In step 4, a direct solver based frequency domain modeling (DSM) is calculated. The finite frequencies required for frequency domain FWI can effectively model a large number of sources, as long as the impedance matrix B for each frequency is LU decomposed. Parallelism in the direct solution method DSM is achieved by using a massive parallel direct solver. Furthermore, it has been shown by our practical experience that acceleration of over 15 is difficult to achieve, no matter how many processors are used by 2D and 3D applications. The memory complexity and time complexity of the two-dimensional finite difference problem direct solver are O (N) 2 log 2 N) and O (N) 3 ) In three dimensions, the expansion is greatly exaggerated to O (N 4 ) And O (N) 6 ). Here, a function is included in brackets after O to indicate the relationship between time/space consumed and the amount of data increase for an algorithm. Where N represents the amount of input data. And the magnitude relation of the complexity is: o (1)<O(log 2 N)<O(N)<O(N 2 )。
In step 5, an iterative solver based frequency domain modeling (ISM) is calculated. The main advantage over DSM is the smaller memory requirement, typically O (N for 3D 3 ). The disadvantage is that the impedance matrix is indefinite (real eigenvalue of the B change sign) and therefore ill-conditioned. With a recent new multi-layer Krylov process, the process produces many iterations that are almost independent of the problem size or frequency.
In this ideal case, for a 3D problem, the time complexity of the iterative solver is O (N 3 ) Whereas previous versions of the method show a temporal complexity of O (N 4 ) (illustrating the number of iterations increases linearly with increasing N)
For multi-source problems and coarse-grained parallelism on sources, the time complexity of the iterative solver is O (N 4 ). Theoretically, an iterative solver should provide the most efficient numerical format. But attenuation effects are easily introduced.
In step 6, a hybrid solver based frequency domain modeling (HSM) is calculated. First, a direct solver is used on each processor to perform sequential logical element decomposition on the local matrix assembled on each sub-domain. Secondly, a simplified system, the so-called Schur complement system, is solved by using an iterative solver such as GMRES, which is better than the complete system condition handled by ISM. The solution is the interface unknowns and rhs is derived from the previous decomposition step. Once the interface unknowns are determined, the internal unknowns may be effectively computed by performing a forward/reverse substitution on each processor. The parameter epsilon (residual divided by the norm of RHS) represents the stopping criterion of GMRES iterations.
In step 7, the results of the three solvers are compared. The number of sources/receivers based on different computer resources and seismic experiments was derived, with each appropriate solver selection.
N p Total number of processors.
N rhs Total source number.
For DSM:
N DSM =number of processors dedicated to one LU decomposition,
T LU time-consuming LU decomposition,
T rhs time-consuming of one rhs solution.
For HSM:
N HSM =number of processors of one domain decomposition,
T LU+M run time of =rhs independent task
T GMRES Runtime of GMRES
For TDM:
T seq run time of =1-rhs sequential simulation.
In our implementation we observe from figure four that TDM is slower than HSM in slope, which makes the former approach superior, although any improvement in iterative preprocessing may challenge this conclusion. The DSM performs better for a large number of sources (over 500 sources) due to the efficiency of the solution steps.
Example 3
In a specific embodiment 3 to which the present invention is applied, fig. 2 is a schematic diagram of the memory and time complexity of the high-efficiency forward algorithm TDM based on the discrete fourier TDM viscoelastic wave equation of the present invention;
FIG. 3 is a schematic diagram of the run-time calculation of the discrete Fourier TDM viscoelastic wave equation based efficient forward algorithm DSM, HSM, TDM of the present invention;
FIG. 4 is a schematic diagram of the time spent by the discrete Fourier TDM based viscoelastic wave equation high efficiency forward algorithm of the present invention as a function of the total number of DSM, HSM and TDM sources. The vertical axis coordinates are the estimated time (hours), and the horizontal axis coordinates are the lens numbers.
Elapsed time is a function of the total number of DSM, HSM and TDM sources. This curve computes the number of three available processors, 192,500 and 2000 (represented in the box). 192 processor packets are used for parallel logical unit decomposition of the DSM and domain decomposition of the HSM. The initial LU factorization is negligible for HSM. For TDM, parallelization is performed by lens.
Finally, it should be noted that: the foregoing description is only a preferred embodiment of the present invention, and is not intended to limit the present invention, but although the present invention has been described in detail with reference to the foregoing embodiment, it will be apparent to those skilled in the art that modifications may be made to the technical solution described in the foregoing embodiment, or equivalents may be substituted for some of the technical features thereof. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.
Other than the technical features described in the specification, all are known to those skilled in the art.
Claims (14)
1. A discrete fourier TDM viscoelastic wave equation based efficient forward algorithm, comprising:
step 1, carrying out frequency domain modeling on an initial model;
step 2: time modeling TDM with an integration scheme having an explicit format;
step 3: calculating a direct solver-based frequency domain modeling, DSM;
step 4: calculating a frequency domain modeling ISM based on an iteration solver;
step 5: calculating a frequency domain modeling HSM based on a hybrid solver;
step 6: a comparison was made for these three methods.
2. The efficient forward algorithm of discrete fourier TDM viscoelastic wave equation based according to claim 1, where in step 1, for two-dimensional frequency domain FWI, the forward problem is implemented in the frequency domain by a direct solver, and in three dimensions, the optimal strategy for forward modeling is less obvious, and frequency response can be extracted from time domain modeling TDM in combination with discrete fourier transform to frequency domain modeling based on direct and hybrid or iterative solvers.
3. The discrete fourier TDM viscoelastic wave equation based efficient positive algorithm according to claim 2 wherein in step 1, the three-dimensional viscous acoustic wave equation is considered in the frequency domain:
wherein density is represented by ρ (x), bulk modulus is represented by κ (x), and angular frequency is represented by ω; the monochromatic pressure wavefield and source are denoted by p (x, ω) and s (x, ω), respectively; in bulk modulus expressions, inherent attenuation can be easily achieved in the frequency domain using complex-valued wave velocities;
formula (1) can be reproduced in matrix form
Ap＝s (2)
Wherein the complex-valued impedance matrix a depends on the angular frequency ω and the median parameters κ and ρ.
4. The efficient forward algorithm of discrete fourier TDM viscoelastic wave equation based on claim 1, where in step 2, the partial differential operator is converted to algebraic operation using existing iterative numerical methods, which can be represented by a matrix, the wave equation attempting to estimate the vector f over time by an explicit system, where the vector f is pressure, velocity of solid particles, or velocity of fluid/solid particles.
5. The discrete fourier TDM viscoelastic wave equation based efficient positive algorithm according to claim 1, wherein in step 2, the matrix expression is:
where xyz is a three-dimensional spatial coordinate parameter and t is a time parameter;
the quality matrix M is a diagonal matrix;
the stiffness matrix a should operate on the desired solution back and forth in the spatial domain or in the spectral domain;
the source term S is usually a local point source in vibroseis seismology, and the corresponding formula in the frequency domain is popularization by adopting a Helmholtz equation;
where xyz is a three-dimensional spatial coordinate parameter and w is angular velocity;
wherein the impedance matrix B is complex valued and has a symmetrical pattern when considering the finite discretization.
6. The discrete fourier TDM viscoelastic wave equation based efficient forward algorithm according to claim 5 wherein in step 2 TDM is performed generally according to an explicit time planning algorithm; in each time step, the solution for each spatial grid point is estimated from the solution for the previous time step.
7. The discrete fourier TDM viscoelastic wave equation based efficient positive algorithm according to claim 6 wherein in step 2, the core memory storage of the full time sequence is useless for the frequency domain FWI because the frequency domain wavefield is derived by discrete fourier summation; the calculation domain is divided into subdomains with limited dimensions by adopting a standard domain decomposition method, so that a time domain algorithm can be effectively parallelized; in the framework of multi-source emulation, the low memory requirements of TDM also allow coarse-grained parallelism to the sources, which can be combined with domain decomposition parallelism if the number of processors is significantly greater than the number of sources.
8. The discrete fourier TDM viscoelastic wave equation based efficient positive algorithm according to claim 7, wherein in step 2, the dimension of the three-dimensional N3 computational grid is denoted as N; if only parallelism on the source is realized, the memory and time complexity of TDM are given; real 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 ) Wherein N is rhs The number of sources is on one side of the grid; the TDM algorithm can accurately achieve propagation in any heterogeneous layer medium without instability as long as a sufficiently fine time discretization method is used; although memory requirements increase significantly, extensions to elasticity, anisotropy, and decay modeling are possible.
9. The efficient forward algorithm of discrete fourier TDM viscoelastic wave equation based on claim 8, wherein in step 2, the dimension of the matrix is the number of unknowns in the computational grid, the number of numerical bandwidths and non-zero coefficient matrices are embedded in the frequency domain using differential operators depending on the numerical discretized templates, wave modeling reduces solving large sparse linear equations and multiple right term RHSs, one for each source, resulting in a linear system.
10. The efficient forward algorithm of discrete fourier TDM viscoelastic wave equation based on claim 9, wherein in step 2, time domain modeling TDM of explicit integration scheme is performed; in the framework of multi-source simulation, the low memory requirements of TDM also allow coarse-grained parallelism on the sources, which can be combined if the number of processors is significantly greater than the number of sources.
11. The efficient forward algorithm of discrete fourier TDM viscoelastic wave equation based on claim 1, wherein in step 3, a large number of sources can be modeled effectively as long as LU decomposition is performed on the impedance matrix B for each frequency, the finite frequency required for frequency domain FWI; parallelism in the direct solution method DSM is achieved by using a massive parallel direct solver; moreover, practical experience shows that no matter how many processors are used by 2D and 3D applications, acceleration of over 15 is difficult to achieve; the memory complexity and time complexity of the two-dimensional finite difference problem direct solver are O (N) 2 log 2 N) and O (N) 3 ) In three dimensions, the expansion is greatly exaggerated to O (N 4 ) And O (N) 6 )。
12. The efficient forward algorithm of discrete fourier TDM viscoelastic wave equation based on claim 1, wherein in step 4, another approach to frequency domain modeling is based on an iterative solver ISM, which is mainly advantageous in terms of smaller memory requirements compared to DSM, typically O (N 3 )。
13. The efficient forward algorithm of discrete fourier TDM viscoelastic wave equation based on claim 1, wherein in step 5, the hybrid frequency domain modeling method, HSM, provides a good tradeoff between DSM and ISM in terms of memory requirements and multiple rhs simulation efficiency; the HSM is based on a region decomposition approach, using a hybrid direct/iterative solver.
14. The efficient forward algorithm of discrete fourier TDM viscoelastic wave equation based on claim 1, wherein in step 6, a suitable solver is selected based on the computer resources and the number of sources and receivers for the seismic experiment.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210566067.6A CN117150703A (en) | 2022-05-23 | 2022-05-23 | Efficient forward algorithm based on discrete Fourier TDM viscoelastic wave equation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210566067.6A CN117150703A (en) | 2022-05-23 | 2022-05-23 | Efficient forward algorithm based on discrete Fourier TDM viscoelastic wave equation |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN117150703A true CN117150703A (en) | 2023-12-01 |
Family
ID=88910657
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210566067.6A Pending CN117150703A (en) | 2022-05-23 | 2022-05-23 | Efficient forward algorithm based on discrete Fourier TDM viscoelastic wave equation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117150703A (en) |
-
2022
- 2022-05-23 CN CN202210566067.6A patent/CN117150703A/en active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Bohlen | Parallel 3-D viscoelastic finite difference seismic modelling | |
| US10459117B2 (en) | Extended subspace method for cross-talk mitigation in multi-parameter inversion | |
| Li et al. | 2D and 3D frequency-domain elastic wave modeling in complex media with a parallel iterative solver | |
| US9195783B2 (en) | Reducing the dimensionality of the joint inversion problem | |
| Larsen et al. | Elastic modeling initiative, Part III: 3-D computational modeling | |
| Klin et al. | Numerical simulation of seismic wave propagation in realistic 3-D geo-models with a Fourier pseudo-spectral method | |
| Liu et al. | Source wavefield reconstruction using a linear combination of the boundary wavefield in reverse time migration | |
| Golubev et al. | Simulation of seismic wave propagation in a multicomponent oil deposit model | |
| MX2011003850A (en) | Image domain signal to noise estimate. | |
| Vamaraju et al. | Enriched Galerkin finite element approximation for elastic wave propagation in fractured media | |
| Wu et al. | Seismic wave propagation and scattering in heterogeneous crustal waveguides using screen propagators: I SH waves | |
| Belonosov et al. | 3D numerical simulation of elastic waves with a frequency-domain iterative solver | |
| Angus | The one-way wave equation: a full-waveform tool for modeling seismic body wave phenomena | |
| Xu et al. | Adaptive 27-point finite-difference frequency-domain method for wave simulation of 3D acoustic wave equation | |
| Ravasi et al. | An open-source framework for the implementation of large-scale integral operators with flexible, modern HPC solutions-enabling 3D Marchenko imaging by least squares inversion | |
| Wu et al. | A new full waveform inversion method based on shifted correlation of the envelope and its implementation based on OPENCL | |
| Meng et al. | Numerical dispersion analysis of discontinuous Galerkin method with different basis functions for acoustic and elastic wave equations | |
| JP2001519932A (en) | Frequency Domain Seismic Data Processing Method Using Massively Parallel Computer | |
| Operto et al. | Is 3D frequency-domain FWI of full-azimuth/long-offset OBN data feasible? The Gorgon data FWI case study | |
| US20160034612A1 (en) | Re-ordered Interpolation and Convolution for Faster Staggered-Grid Processing | |
| Xu et al. | A simplified calculation for adaptive coefficients of finite-difference frequency-domain method | |
| Raknes et al. | Challenges and solutions for performing 3D time-domain elastic full-waveform inversion | |
| CN117150703A (en) | Efficient forward algorithm based on discrete Fourier TDM viscoelastic wave equation | |
| Bassi et al. | High order ADER-DG schemes for the simulation of linear seismic waves induced by nonlinear dispersive free-surface water waves | |
| Malkoti et al. | A highly efficient implicit finite difference scheme for acoustic wave propagation |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |