CN117150835A - Method for simulating seismic wave value of viscous-acoustic anisotropic medium - Google Patents
Method for simulating seismic wave value of viscous-acoustic anisotropic medium Download PDFInfo
- Publication number
- CN117150835A CN117150835A CN202210566039.4A CN202210566039A CN117150835A CN 117150835 A CN117150835 A CN 117150835A CN 202210566039 A CN202210566039 A CN 202210566039A CN 117150835 A CN117150835 A CN 117150835A
- Authority
- CN
- China
- Prior art keywords
- node
- seismic wave
- viscous
- nodes
- split
- 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
- 238000000034 method Methods 0.000 title claims abstract description 84
- 239000011159 matrix material Substances 0.000 claims abstract description 37
- 238000004088 simulation Methods 0.000 claims abstract description 21
- 238000004364 calculation method Methods 0.000 claims description 13
- 238000000354 decomposition reaction Methods 0.000 claims description 7
- 230000003044 adaptive effect Effects 0.000 claims description 4
- 230000001131 transforming effect Effects 0.000 claims description 4
- 238000004422 calculation algorithm Methods 0.000 claims description 3
- 238000010276 construction Methods 0.000 claims description 2
- 230000008569 process Effects 0.000 abstract description 5
- 238000010586 diagram Methods 0.000 description 6
- 239000006185 dispersion Substances 0.000 description 5
- 238000012545 processing Methods 0.000 description 5
- 238000011161 development Methods 0.000 description 3
- 239000012530 fluid Substances 0.000 description 3
- 239000007787 solid Substances 0.000 description 3
- 238000001514 detection method Methods 0.000 description 2
- 238000006073 displacement reaction Methods 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
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 230000000704 physical effect 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
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/282—Application of seismic models, synthetic seismograms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/58—Media-related
- G01V2210/586—Anisotropic media
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- 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
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Remote Sensing (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention relates to a method for simulating seismic wave values of a viscous-acoustic anisotropic medium, which comprises the following steps: step 1, performing node model subdivision based on a speed model; step 2, extracting parameters and storing related parameters; step 3, constructing a new radial basis function differential matrix; step 4, calculating a differential weighting coefficient; step 5, solving a viscous-acoustic anisotropic equation based on a gridless finite difference method, and calculating a seismic wave field; and 6, outputting a seismic wave field. According to the method, flexible node subdivision is carried out based on a speed model, the influence of deformation parameters is eliminated by improving a differential matrix solving process, and a stable differential coefficient is brought into a viscous-acoustic anisotropic equation to obtain a high-precision seismic wave numerical simulation result.
Description
Technical Field
The invention belongs to the field of seismology, relates to the problem of a seismic wave propagation mechanism, and particularly relates to a numerical simulation method for a viscous-acoustic anisotropic medium seismic wave.
Background
The high-precision seismic wave numerical simulation method is a technology for simulating the propagation rule and waveform characteristics of seismic waves in underground medium by mathematical means of numerical approximation. It is not only an important tool and theoretical basis for seismic migration and inversion, but also occupies a vital place in the development of seismic acquisition, processing and subsequent interpretation work to explore subsurface geologic structures and energy resources.
With the continuous development of exploration development and computer performance in China, the exploration target area is changed from a simple structure to a complex hidden oil and gas reservoir, and the complex surface conditions and the complex underground structure of the exploration target area bring great challenges for seismic exploration. Early seismic wave numerical simulations made a number of approximations and assumptions about subsurface media, and were generally considered isotropic media. Whereas actual real subsurface media are exceptionally complex, most subsurface media are anisotropic, which may be caused by the preferential orientation of anisotropic mineral particles and cracks, or lamellar layer theory of isotropic or anisotropic layers, so wave equations derived based on isotropic media assumptions have significant seismic dynamics and kinematic errors in real media. While media with complex physical properties often coexist with severely undulating terrain configurations, this makes accurate simulation of seismic wave propagation under complex geological conditions a great challenge.
In the field of viscoelastic anisotropic media, the former has made a lot of valuable research work. Thomson (1986) demonstrates that most elastic media are weakly anisotropic, laying a theoretical basis for the approximation of anisotropic wave equations. Tsvankin (1996) describes the P-wave characteristics in a Transverse Isotropic (TI) medium with arbitrary anisotropy strength using Thomsen notation. Tsvenk in (1997) replaces stiffness coefficients with two perpendicular (P and S) velocities and seven anisotropic parameters, simplifying the analytical description of wave propagation in orthogonal media. To further simplify the anisotropic wave equation, alkhalifah (1998) proposes an acoustic approximation of the longitudinal wave, assuming that the velocity of the transverse wave along the axis of symmetry is zero, which greatly simplifies the expression describing the elastic model. With acoustic approximation (Alkhalifah, 1998), alkhalifah (2000) proposes the acoustic wave equation for a perpendicular transverse isotropy (VTI) medium. Likewise, alkhalifah (2003) further proposes an acoustic wave equation media acoustic approximation (Alkhalifah, 1998). Since its discovery, it has found wide application in industry and academia because it accurately describes the kinematics of longitudinal wave propagation in elastic media, significantly improving the efficiency of imaging, inversion and other seismic processing applications (Xu et al 2020). On the basis of the above study, we also note that the conventional finite difference method is more applicable to nodes regularly distributed in square or rectangular lattices, which can accurately simulate conventional orthogonal and symmetric media. However, for the complex geologic model of the anisotropic medium, the nodes which are regularly distributed are not suitable for the problems of strong false scattering, poor stability and the like which are easily caused in the numerical simulation process due to the lack of flexibility of mesh subdivision. In recent years, the gridless finite difference method has been developed greatly. Unlike the finite difference differential coefficients obtained by solving the polynomial constructed differential matrix at nodes of a regular distribution, the grid-free finite difference method differential coefficients are obtained by solving the differential matrix constructed from different radial basis functions at a set of nodes of a disordered distribution. However, the grid-free finite difference method based on radial basis functions is prone to stability problems. The stability of interest refers to the numerical disease problem of the differential matrix caused by different selection of shape parameters when the differential coefficient of the gridless finite difference method is calculated on different subdivision nodes. Therefore, a new seismic wave numerical simulation method with freely discrete nodes, strong stability and high precision and suitable for the viscous-acoustic anisotropic medium must be established.
In application number: in chinese patent application No. cn201310741415.X, a method and apparatus for simulating seismic waves in complex anisotropic media are disclosed, the method comprising: assigning an elastic constant and a density to each grid point of the digitized three-dimensional geological model; acquiring first-order spatial derivatives of velocity components of grid points along three coordinate axes in a three-dimensional geological model; substituting the first-order spatial derivative of the velocity component into a first equation of an elastic wave equation set in an elastic medium, and solving to obtain a stress component; acquiring first-order spatial derivatives of stress components of grid points along three coordinate axes in a three-dimensional geological model; substituting the first-order spatial derivative of the stress component into a second equation in an elastic wave equation set in an elastic medium, and solving to obtain a speed component; obtaining a value of a velocity component at a first moment of the current time according to the velocity component obtained through the second equation, and recording seismic waves at the current moment; and acquiring the record of the seismic wave at the next moment until the wave field value of the seismic wave with the time length t is acquired.
In application number: in the chinese patent application CN202110705392.1, a method, an apparatus and a device for simulating seismic wave values of a viscous anisotropic medium are related, where the method includes: determining a first quality factor matrix of the viscoelastic medium according to the first stiffness matrix of the known elastic medium; simplifying the first rigidity matrix and the first quality factor matrix to obtain a second rigidity matrix and a second quality factor matrix; according to the two stiffness matrixes and the second quality factor matrix, a preset dispersion relation, a motion balance equation and a geometric equation are combined, and a viscous anisotropic medium acoustic wave equation is obtained through calculation; and determining parameters of the underground medium and parameters of the seismic waves, and bringing the parameters of the underground medium and the parameters of the seismic waves into the viscous anisotropic medium acoustic wave equation to calculate and obtain a simulation value of the seismic wave field, so that the change characteristics of the velocity and attenuation of the seismic waves in the anisotropic medium along with the direction can be realized, and the accuracy of the seismic wave numerical simulation is improved.
In application number: in the chinese patent application CN201510473854.6, a method for simulating the three-dimensional TTI biphase medium seismic wave field based on finite difference method is involved, wherein the method comprises obtaining solid and fluid stress tensors and solid and fluid strain tensors, and converting them into constitutive equation; obtaining a geometric equation according to the corresponding relation between stress and displacement; obtaining a motion differential equation according to the structure, the geometric equation and the corresponding relation between the motion of the fluid relative to the solid and the stress and the displacement; taking the divergence of the two sides of the motion differential equation to obtain a first longitudinal wave equation and a second longitudinal wave equation of the seismic wave; the first longitudinal wave equation and the second longitudinal wave equation are subjected to differential dispersion by adopting a 2N-order precision expansion type and a second-order precision center differential format, so that the partial derivative of y is equal to zero, and the first differential equation and the second differential equation are obtained; and performing absorption boundary condition processing on the first and second difference equations to obtain corresponding seismic wave field values. By the method, the real-time propagation simulation of the physical seismic wave field is realized.
The prior art is greatly different from the invention, the technical problem which is needed to be solved by the invention is not solved, and a novel method for simulating the seismic wave value of the viscous-acoustic anisotropic medium is invented.
Disclosure of Invention
Aiming at the defects existing in the prior art, the invention provides the seismic wave numerical simulation method for the viscous-acoustic anisotropic medium, which has the advantages of free node dispersion, strong stability and high precision and is suitable for the viscous-acoustic anisotropic medium.
The aim of the invention can be achieved by the following technical measures: the method for simulating the seismic wave values of the viscous-acoustic anisotropic medium comprises the following steps of:
step 1, performing node model subdivision based on a speed model;
step 2, extracting parameters and storing related parameters;
step 3, constructing a new radial basis function differential matrix;
step 4, calculating a differential weighting coefficient;
step 5, solving a viscous-acoustic anisotropic equation based on a gridless finite difference method, and calculating a seismic wave field;
and 6, outputting a seismic wave field.
The aim of the invention can be achieved by the following technical measures:
in step 1, the seismic velocity of the computational domain is node split using an adaptive node splitting method based on a known seismic velocity model.
In step 1, a layer of nodes at the bottom edge of a calculation domain is set as an initial node, node subdivision is carried out from the bottom edge upwards layer by layer, the space position (x, z) of the node at the moment is obtained, the subdivision radius r of the node and the seismic wave velocity v corresponding to the node are obtained p0 Tilt angle θ and corresponding quality factor Q, anisotropic media parameters epsilon and sigma.
In step 2, the space positions (x, z) of all nodes in the domain are calculated after the subdivision is extracted, the subdivision radius r of the nodes, and the seismic velocity v corresponding to the nodes p0 The dip angle theta and the corresponding quality factor Q, the anisotropic medium parameters epsilon and sigma; and (2) searching N points nearest to each split node in the step (1) in the range of all nodes in the calculation domain by using a K nearest neighbor algorithm, and returning to the corresponding node space position (x) i ,z i ) I=1, 2,..n and distance d from the point i I=1, 2,; finally, the parameters are stored.
In step 2, the stored parameters specifically include the spatial position (x, z) of each split node, the split radius r of the node, and the seismic velocity v corresponding to the node p0 Tilt angle θ and corresponding quality factor Q, anisotropic medium parameters epsilon and sigma, and spatial position (x i ,z i ) I=1, 2,.. i I=1, 2, N, i.e. store information { (x, z), r, v for a current node p0 ,θ,Q,ε,σ,(x i ,z i ),d i }。
In step 3, a uniform seismic wave velocity model is adopted, split nodes and parameters corresponding to the uniform seismic wave velocity model are obtained based on a node splitting method, and a new radial basis function differential matrix is constructed.
In step 3, if a new radial basis function differential matrix of the current node is constructed, first, spatial information (x, z) and (x) in the information are stored for each split node in step 2 i ,z i ) And carrying out uniform scaling, carrying out chebyshev polynomial decomposition on the conventional Gaussian radial basis function, and transforming the conventional Gaussian radial basis function into a polar coordinate system (r, theta).
In step 3, chebyshev polynomial decomposition is performed on a conventional gaussian radial basis function, and the form after the chebyshev polynomial decomposition is transformed into a polar coordinate system (r, θ) is specifically as follows:
wherein d j,m As scale factor, c j,m (x k ) Sum s j,m (x k ) As polynomial coefficients, chebyshev parametersAnd T j-2m (r) is:
the parameter j is more than 0, and the parameter 0 is more than or equal to m and less than or equal to (j-p)/2; when the parameter j is even, the parameter p=0, and when the parameter j is odd, the scale factor d j,m Is that
Wherein, the deformation parameter epsilon is properly valued, and the valued range is [0,1];
parameter c in equation (1) j,m (x k ) Sum parameter s j,m (x k ) In particular to
Wherein parameter b 0 =1, parameter b m =2, m > 0; parameter t 0 ＝1/2，t j ＝1,j＞0；F 2 To go beyond the system of equations, where the parameter alpha j,m = (j-2m+p+1)/2, parameter β j,m ＝[j-2m+1,(j+2m+p+1)/2]The method comprises the steps of carrying out a first treatment on the surface of the Therefore, the construction of a new radial basis function differential matrix by the conventional Gaussian radial basis function is specifically as follows:
wherein C represents a polynomial coefficient C j,m (x k ) Sum s j,m (x k ) The coefficient matrix is composed, D represents the scale factor D j,m The scale matrix is formed by the method,representing the parameter ∈Chebyshev>And T j-2m (r) a chebyshev matrix formed together. In step 4, according to the accuracy requirement, selecting a proper differential node number, based on step 3, utilizing a new radial basis function differential matrix corresponding to the current split node number, and obtaining the split node radius r obtained in step 2 and the nearest N split node space positions (x i ,z i ) Distance d from i I=1, 2, N, calculating accurate differential weighting coefficient omega corresponding to current node 1 i And omega 2 i I=1, 2, N, and further obtaining differential weighting coefficients corresponding to all the subdivision nodes.
In step 4, the differential weighting coefficient is specifically:
in the equation (6), x represents the space coordinates of the seismic wave velocity model subdivision nodes, the subscript "0" represents the 0 th nearest neighbor subdivision node around the current seismic wave velocity model subdivision node differential operation, and the subscript "N" represents the N th nearest neighbor subdivision node around the current seismic wave velocity model subdivision node differential operation; for each split node and N nearest split nodes in the step 2, calculating differential weighting coefficients omega corresponding to the N nearest split nodes around the split node by using an equation (6) 1 i And omega 2 i The method comprises the steps of carrying out a first treatment on the surface of the Obtaining N nearest split node differential weighting coefficients corresponding to all split nodes, L [. Cndot.]Representing the partial derivative operator, L [. Cndot.]Delta and respectivelyWhere delta is the laplace operator of the VTI medium.
In step 5, based on the velocity and the viscous-sound anisotropy parameters { (x, z), r, v stored in step 2 p0 ,ρ,Q,ε,σ,(x i ,z i ),d i Using N nearest split node differential weighting coefficients omega corresponding to all split nodes obtained in the step 4 1 i And omega 2 i I=1, 2, N, forward modeling of the seismic wave field is performed by adopting a viscous-acoustic anisotropic equation, obtaining a numerical simulation qP wave seismic wave field U.
In step 5, forward modeling of the seismic wave field is performed by adopting a viscous-acoustic anisotropy equation, wherein the specific equation of the viscous-acoustic anisotropy equation is as follows:
wherein U in the equation 7 represents the wave field of the qP wave, and t is the propagation time of the seismic wave field; v (V) P0 Is the phase velocity of the qP wave along the direction of the symmetry axis, τ is the relaxation coefficient, X ands is an auxiliary scalar operator; the earthquake wave fields corresponding to all the subdivision nodes of the calculation domain are obtained, and a specific calculation equation is as follows:
the invention relates to a viscous-acoustic anisotropic medium seismic wave numerical simulation method, which has the following advantages: 1) Compared with the conventional seismic wave numerical simulation method, the method has the advantages of free node dispersion, strong stability, high precision and adaptability to the viscous-acoustic anisotropic medium; 2) The invention carries out flexible node subdivision based on a speed model, provides a novel grid-free finite difference method based on a radial basis function, eliminates the influence of deformation parameters by improving a difference matrix solving process, brings a stable difference coefficient into a viscous-acoustic anisotropic equation to obtain a high-precision seismic wave numerical simulation result, and can effectively improve the precision and stability of seismic data processing by using the method, thereby providing important theoretical and technical support for oil-gas resource exploration and deep structure detection of the viscous-acoustic anisotropic medium; 3) The invention can be widely applied to the fields of seismology and exploration seismology, and has important theory and application value for geophysical research involving seismic wave numerical simulation.
Drawings
FIG. 1 is a schematic flow chart of a method for simulating seismic wave values of a viscous-acoustic anisotropic medium;
FIG. 2 is a schematic diagram of an anisotropic media Hess model in accordance with an embodiment of the present invention;
FIG. 3 is a schematic diagram of a seismic wavefield snapshot of the anisotropic media Hess model of FIG. 2, in accordance with an embodiment of the invention;
FIG. 4 is a schematic representation of a simulated seismic recording of the anisotropic media Hess model of FIG. 2 in accordance with an embodiment of the invention;
FIG. 5 is a schematic diagram of a Marmousi-2 model of an anisotropic medium according to an embodiment of the invention;
FIG. 6 is a schematic diagram of a seismic simulation wavefield snapshot of the anisotropic media Marmousi-2 model of FIG. 5, in accordance with an embodiment of the invention;
FIG. 7 is a schematic diagram of a simulated seismic recording of the Marmousi-2 model of the anisotropic medium of FIG. 5 in accordance with an embodiment of the invention.
Detailed Description
The present invention will be described in further detail with reference to the drawings and examples, in order to make the objects, technical solutions and advantages of the present invention more apparent. It should be understood that the specific embodiments described herein are for purposes of illustration only and are not intended to limit the scope of the invention.
The invention relates to a method for simulating seismic wave values of a viscous-acoustic anisotropic medium strain, which comprises the following steps: performing node model subdivision based on the speed model; extracting parameters and storing related parameters; constructing a new radial basis function differential matrix; calculating a differential weighting coefficient; solving a viscous-acoustic anisotropic equation based on a gridless finite difference method, and calculating a seismic wave field; outputting the seismic wave field. The invention provides a novel grid-free finite difference method based on a radial basis function based on flexible node subdivision, which eliminates the influence of deformation parameters by improving a difference matrix solving process, brings stable difference coefficients into a viscous-acoustic anisotropic equation to obtain a high-precision seismic wave numerical simulation result, and can effectively improve the precision and stability of seismic data processing by using the method, thereby providing important theoretical and technical support for oil-gas resource exploration and deep structure detection of the viscous-acoustic anisotropic medium.
The following are several embodiments of the invention
Example 1
In an embodiment 1 to which the present invention is applied, please refer to fig. 1, which is a schematic flow chart of a method for simulating seismic wave values of a viscous-acoustic anisotropic medium according to the present invention. The method for simulating the seismic wave value of the viscous-acoustic anisotropic medium comprises the following steps of:
step S100, performing node model subdivision based on a speed model;
based on known seismic velocity modelsThe seismic velocity of the computational domain is node split using An adaptive node splitting method (Duan et al, 2021-An adaptive node-distribution method for radial-basic-function fine-difference modeling with optimal shape parameter), specifically: setting a layer of nodes at the bottom edge of a calculation domain as initial nodes, carrying out node subdivision from the bottom edge to the top layer by layer to obtain the space position (x, z) of the nodes at the moment, the subdivision radius r of the nodes and the seismic wave velocity v corresponding to the nodes p0 Tilt angle θ and corresponding quality factor Q, anisotropic media parameters epsilon and sigma.
Step S200, extracting parameters and storing related parameters;
after the subdivision is extracted, the spatial positions (x, z) of all nodes of the domain are calculated, the subdivision radius r of the nodes and the seismic wave velocity v corresponding to the nodes are calculated p0 The dip angle theta and the corresponding quality factor Q, the anisotropic medium parameters epsilon and sigma; for each split node described in step S100, using a K-nearest neighbor algorithm, searching N points nearest to the node in all node ranges of the calculation domain, and returning to the corresponding node spatial position (x) i ,z i ) I=1, 2,..n and distance d from the point i I=1, 2,; finally, the parameters are stored, and specifically comprise the spatial positions (x, z) of all the split nodes, the split radius r of the nodes and the seismic wave velocity v corresponding to the nodes p0 Tilt angle θ and corresponding quality factor Q, anisotropic medium parameters epsilon and sigma, and spatial position (x i ,z i ) I=1, 2,.. i I=1, 2, N, i.e. store information { (x, z), r, v for a current node p0 ,θ,Q,ε,σ,(x i ,z i ),d i }。
Step S300, a new radial basis function differential matrix is constructed;
and adopting the uniform seismic wave velocity model, obtaining split nodes and parameters corresponding to the uniform seismic wave velocity model based on the node splitting method, and constructing a new radial basis function differential matrix. If a new radial basis function differential matrix of the current node is constructed, firstly, aiming at the spatial information in the storage information of each split node in the step S200x, z) and (x i ,z i ) And carrying out uniform scaling, carrying out chebyshev polynomial decomposition on the conventional Gaussian radial basis function, and transforming the conventional Gaussian radial basis function into a polar coordinate system (r, theta). The method comprises the steps of carrying out Chebyshev polynomial decomposition on a conventional Gaussian radial basis function, and transforming the conventional Gaussian radial basis function into a polar coordinate system (r, theta):
wherein the Chebyshev parameterAnd T j-2m (r) is:
the parameter j is more than 0, and the parameter 0 is less than or equal to m and less than or equal to (j-p)/2. When the parameter j is even, the parameter p=0, and when the parameter j is odd, the scale factor d j,m Is that
Wherein, the deformation parameter epsilon is properly valued, and the valued range is [0,1].
Parameter c in equation (1) j,m (x k ) Sum parameter s j,m (x k ) In particular to
Wherein parameter b 0 =1, parameter b m =2, m > 0; parameter t 0 ＝1/2，t j ＝1,j＞0；F 2 To go beyond the system of equations, where the parameter alpha j,m = (j-2m+p+1)/2, parameter β j,m ＝[j-2m+1,(j+2m+p+1)/2]The method comprises the steps of carrying out a first treatment on the surface of the Thus, the conventional Gaussian radial basis function structureThe new radial basis function differential matrix is built specifically as follows:
step S400, calculating a differential weighting coefficient;
according to the precision requirement, selecting proper differential node numbers, based on the step S300, utilizing a new radial basis function differential matrix corresponding to the current split node number, and obtaining the split node radius r obtained in the step S200 and the nearest N split node space positions (x i ,z i ) Distance d from i I=1, 2, N, calculating accurate differential weighting coefficient omega corresponding to current node 1 i And omega 2 i I=1, 2, N, and further obtaining differential weighting coefficients corresponding to all the split nodes, wherein the differential weighting coefficients are specifically as follows:
in the equation (6), x represents the space coordinates of the seismic wave velocity model subdivision nodes, the subscript "0" represents the 0 th nearest neighbor subdivision node around the current seismic wave velocity model subdivision node differential operation, and the subscript "N" represents the N th nearest neighbor subdivision node around the current seismic wave velocity model subdivision node differential operation; for each split node and N nearest neighboring split nodes in step S200, calculating differential weighting coefficients omega corresponding to N nearest neighboring split nodes around the split node by using equation (6) 1 i And omega 2 i The method comprises the steps of carrying out a first treatment on the surface of the Obtaining N nearest split node differential weighting coefficients corresponding to all split nodes, L [. Cndot.]Representing the partial derivative operator, L [. Cndot.]Delta and respectivelyWhere delta is the laplace operator of the VTI medium.
Step S500, solving a viscous-acoustic anisotropic equation based on a gridless finite difference method, and calculating a seismic wave field;
based on the velocity and the visco-acoustic anisotropy parameters { (x, z), r, v, which correspond to the respective split nodes stored in step S200 p0 ,ρ,Q,ε,σ,(x i ,z i ),d i Differential weighting coefficient omega of N nearest split nodes corresponding to all split nodes obtained in step S400 1 i And omega 2 i I=1, 2, & gt, N, performing forward modeling of the seismic wave field by adopting a viscous-acoustic anisotropic equation to obtain a numerically-modeled qP wave seismic wave field U;
the seismic wave field forward modeling is performed by adopting a viscous-acoustic anisotropic equation, and the specific equation of the viscous-acoustic anisotropic equation is as follows:
wherein U in the equation 7 represents the wave field of the qP wave, and t is the propagation time of the seismic wave field; v (V) P0 Is the phase velocity of the qP wave along the direction of the symmetry axis, τ is the relaxation coefficient, and X and S are auxiliary scalar operators. The earthquake wave fields corresponding to all the subdivision nodes of the calculation domain are obtained, and a specific calculation equation is as follows:
step S600, outputting a wave field U of the qP wave.
Example 2
In a specific embodiment 2 to which the present invention is applied, fig. 2 is an anisotropic medium Hess model: where fig. 2 (a) is the velocity, fig. 2 (b) is the quality factor, fig. 2 (c) is the anisotropy parameter epsilon, and fig. 2 (d) is the anisotropy parameter sigma. FIG. 3 is a seismic wavefield snapshot section of the anisotropic media Hess model shown in FIG. 2: wherein FIG. 3 (a) is a 1.5s seismic wavefield snapshot taken using the pure qP wave equation and FIG. 3 (b) is a 1.5s seismic wavefield snapshot taken using the method of the present invention; as can be seen from fig. 3 (a) and 3 (b), compared with the seismic record obtained by solving the pure qP wave VTI wave equation, the seismic wave obtained by the method of the present invention propagates stably in the viscous-acoustic anisotropic medium, and the wave front surface has the characteristics of attenuation amplitude and delay phase, which well describes the viscosity and anisotropy of the medium. FIG. 4 is a seismic recording section of the anisotropic media Hess model of FIG. 2: wherein FIG. 4 (a) is a seismic record profile obtained using the pure qP wave equation and FIG. 4 (b) is a seismic record profile obtained using the method of the present invention; the same conclusion can also be drawn from fig. 4 (a) and 4 (b), demonstrating the feasibility and effectiveness of the method of the invention.
Example 3
In a specific embodiment 3 to which the present invention is applied, fig. 5 is a schematic diagram of an anisotropic medium Marmousi-2 model provided by the present invention: wherein fig. 5 (a) is the velocity, fig. 5 (b) is the quality factor, fig. 5 (c) is the tilt angle, fig. 5 (d) is the anisotropy parameter epsilon, and fig. 5 (e) is the anisotropy parameter sigma; FIG. 6 is a seismic wavefield snapshot section of the anisotropic media Marmousi-2 model shown in FIG. 5: wherein FIG. 6 (a) is a 1.0s seismic wavefield snapshot obtained using the method of the present invention, and FIG. 6 (b) is a 1.0s seismic wavefield snapshot obtained using the method of the present invention; from fig. 7, it can be seen that the wavefront surface of the method has obvious energy loss and phase delay at the arrival depth of 1km to 2km, which proves the effectiveness of the method; FIG. 7 is a seismic recording section of the anisotropic media Marmousi-2 model of FIG. 5: wherein FIG. 7 (a) is a seismic record profile obtained using the pure qP wave equation and FIG. 7 (b) is a seismic record profile obtained using the method of the present invention; as can be seen from FIG. 7, the seismic wave field of the method can stably propagate in the viscous-acoustic anisotropic medium, and has the characteristics of free dispersion of nodes, strong stability, high precision and adaptability to the viscous-acoustic anisotropic medium. Further demonstrating the effectiveness of the method of this aspect.
According to the method, flexible node subdivision is carried out based on a speed model, the influence of deformation parameters is eliminated by improving a differential matrix solving process, and a stable differential coefficient is brought into a viscous-acoustic anisotropic equation to obtain a high-precision seismic wave numerical simulation result.
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 (12)
1. The method for simulating the seismic wave value of the viscous-acoustic anisotropic medium is characterized by comprising the following steps of:
step 1, performing node model subdivision based on a speed model;
step 2, extracting parameters and storing related parameters;
step 3, constructing a new radial basis function differential matrix;
step 4, calculating a differential weighting coefficient;
step 5, solving a viscous-acoustic anisotropic equation based on a gridless finite difference method, and calculating a seismic wave field;
and 6, outputting a seismic wave field.
2. A method of modeling seismic wave numerical values of a viscous-acoustic anisotropic medium according to claim 1, wherein in step 1, the seismic wave velocities of the computational domain are node-split using an adaptive node-splitting method based on a known seismic wave velocity model.
3. The method for simulating seismic wave values of a viscous-acoustic anisotropic medium according to claim 2, wherein in step 1, a layer of nodes at the bottom edge of a calculation domain is set as a starting node, and node subdivision is performed layer by layer from the bottom edge to the topObtaining the space position (x, z) of the node at the moment, the subdivision radius r of the node and the seismic wave velocity v corresponding to the node p0 Tilt angle θ and corresponding quality factor Q, anisotropic media parameters epsilon and sigma.
4. The method for simulating seismic wave numerical values of a viscous-acoustic anisotropic medium according to claim 1, wherein in the step 2, the spatial positions (x, z) of all nodes of the calculated domain after the subdivision are extracted, the subdivision radius r of the nodes, and the seismic wave velocity v corresponding to the nodes p0 The dip angle theta and the corresponding quality factor Q, the anisotropic medium parameters epsilon and sigma; and (2) searching N points nearest to each split node in the step (1) in the range of all nodes in the calculation domain by using a K nearest neighbor algorithm, and returning to the corresponding node space position (x) i ,z i ) I=1, 2,..n and distance d from the point i I=1, 2,; finally, the parameters are stored.
5. The method of seismic wave numerical simulation of a viscous-acoustic anisotropic medium according to claim 4, wherein in step 2, the stored parameters include the spatial position (x, z) of each split node, the split radius r of the node, and the seismic wave velocity v corresponding to the node p0 Tilt angle θ and corresponding quality factor Q, anisotropic medium parameters epsilon and sigma, and spatial position (x i ,z i ) I=1, 2,.. i I=1, 2, N, i.e. store information { (x, z), r, v for a current node p0 ,θ,Q,ε,σ,(x i ,z i ),d i }。
6. The method for simulating the seismic wave numerical values of the viscous-acoustic anisotropic medium according to claim 1, wherein in the step 3, a uniform seismic wave velocity model is adopted, split nodes and parameters corresponding to the uniform seismic wave velocity model are obtained based on a node splitting method, and a new radial basis function differential matrix is constructed.
7. According to claim 6A method for simulating seismic wave values of a viscous-acoustic anisotropic medium is characterized in that in step 3, if a new radial basis function differential matrix of a current node is constructed, spatial information (x, z) and (x) in information are stored for each split node in step 2 i ,z i ) And carrying out uniform scaling, carrying out chebyshev polynomial decomposition on the conventional Gaussian radial basis function, and transforming the conventional Gaussian radial basis function into a polar coordinate system (r, theta).
8. The method for simulating the numerical value of a seismic wave of a viscous-acoustic anisotropic medium according to claim 7, wherein in the step 3, a chebyshev polynomial decomposition is performed on a conventional gaussian radial basis function, and the form after the conventional gaussian radial basis function is transformed into a polar coordinate system (r, θ) is specifically as follows:
wherein d j,m As scale factor, c j,m (x k ) Sum s j,m (x k ) As polynomial coefficients, chebyshev parametersAnd T j-2m (r) is:
the parameter j is more than 0, and the parameter 0 is more than or equal to m and less than or equal to (j-p)/2; when the parameter j is even, the parameter p=0, and when the parameter j is odd, the scale factor d j,m Is that
Wherein, the deformation parameter epsilon is properly valued, and the valued range is [0,1];
polynomial coefficient c in equation (1) j,m (x k ) And polynomial coefficient s j,m (x k ) In particular to
Wherein parameter b 0 =1, parameter b m =2, m > 0; parameter t 0 ＝1/2，t j ＝1,j＞0；F 2 To go beyond the system of equations, where the parameter alpha j,m = (j-2m+p+1)/2, parameter β j,m ＝[j-2m+1,(j+2m+p+1)/2]The method comprises the steps of carrying out a first treatment on the surface of the Therefore, the construction of a new radial basis function differential matrix by the conventional Gaussian radial basis function is specifically as follows:
wherein C represents a polynomial coefficient C j,m (x k ) Sum s j,m (x k ) The coefficient matrix is composed, D represents the scale factor D j,m The scale matrix is formed by the method,representing the parameter ∈Chebyshev>And T j-2m (r) a chebyshev matrix formed together.
9. The method according to claim 1, wherein in step 4, an appropriate number of differential nodes is selected according to the accuracy requirement, based on step 3, a new radial basis function differential matrix corresponding to the current number of split nodes is utilized, and the radius r of the split node obtained in step 2 and the nearest N spatial positions (x i ,z i ) Distance d from i I=1, 2, N, calculating accurate differential weighting coefficient omega corresponding to current node 1 i And omega 2 i I=1, 2, N, and further obtaining differential weighting coefficients corresponding to all the subdivision nodes.
10. The method of seismic wave numerical simulation of a viscous-acoustic anisotropic medium according to claim 9, wherein in step 4, the differential weighting coefficients are specifically:
in the equation (6), x represents the space coordinates of the seismic wave velocity model subdivision nodes, the subscript "0" represents the 0 th nearest neighbor subdivision node around the current seismic wave velocity model subdivision node differential operation, and the subscript "N" represents the N th nearest neighbor subdivision node around the current seismic wave velocity model subdivision node differential operation; for each split node and N nearest split nodes in the step 2, calculating differential weighting coefficients omega corresponding to the N nearest split nodes around the split node by using an equation (6) 1 i And omega 2 i The method comprises the steps of carrying out a first treatment on the surface of the Obtaining N nearest split node differential weighting coefficients corresponding to all split nodes, L [. Cndot.]Representing the partial derivative operator, L [. Cndot.]Delta and respectivelyWhere delta is the laplace operator of the VTI medium.
11. The method according to claim 1, wherein in step 5, based on the velocity and the viscosity anisotropy parameters { (x, z), r, v for each split node stored in step 2 p0 ,ρ,Q,ε,σ,(x i ,z i ),d i Using N nearest split node differential weighting coefficients omega corresponding to all split nodes obtained in the step 4 1 i And omega 2 i I=1, 2, N, seismic wave using viscous-acoustic anisotropy equationAnd performing field forward modeling to obtain a qP wave seismic wave field U of numerical modeling.
12. The method for simulating seismic wave numerical values of a viscous-acoustic anisotropic medium according to claim 11, wherein in step 5, a forward modeling of the seismic wave field is performed using a viscous-acoustic anisotropic equation, and the specific equation of the viscous-acoustic anisotropic equation is:
wherein U in the equation 7 represents the wave field of the qP wave, and t is the propagation time of the seismic wave field; v (V) P0 Is the phase velocity of the qP wave along the direction of the symmetry axis, tau is the relaxation coefficient, and X and S are auxiliary scalar operators; the earthquake wave fields corresponding to all the subdivision nodes of the calculation domain are obtained, and a specific calculation equation is as follows:
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210566039.4A CN117150835A (en) | 2022-05-23 | 2022-05-23 | Method for simulating seismic wave value of viscous-acoustic anisotropic medium |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210566039.4A CN117150835A (en) | 2022-05-23 | 2022-05-23 | Method for simulating seismic wave value of viscous-acoustic anisotropic medium |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN117150835A true CN117150835A (en) | 2023-12-01 |
Family
ID=88904745
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210566039.4A Pending CN117150835A (en) | 2022-05-23 | 2022-05-23 | Method for simulating seismic wave value of viscous-acoustic anisotropic medium |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117150835A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN118330732A (en) * | 2024-04-30 | 2024-07-12 | 中国地震局地球物理研究所 | Earthquake positioning method based on three-dimensional TTI medium model |
-
2022
- 2022-05-23 CN CN202210566039.4A patent/CN117150835A/en active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN118330732A (en) * | 2024-04-30 | 2024-07-12 | 中国地震局地球物理研究所 | Earthquake positioning method based on three-dimensional TTI medium model |
| CN118330732B (en) * | 2024-04-30 | 2024-08-27 | 中国地震局地球物理研究所 | Earthquake positioning method based on three-dimensional TTI medium model |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104122585B (en) | Seismic forward simulation method based on elastic wave field resolution of vectors and low-rank decomposition | |
| CN108802813B (en) | A kind of multi-component seismic data offset imaging method and system | |
| Ladopoulos | Elastodynamics for Non-linear Seismic Wave Motion in Real-Time Expert Seismology | |
| Minisini et al. | Local time stepping with the discontinuous Galerkin method for wave propagation in 3D heterogeneous media | |
| Liu et al. | Higher-order triangular spectral element method with optimized cubature points for seismic wavefield modeling | |
| Golubev et al. | Simulation of seismic wave propagation in a multicomponent oil deposit model | |
| CN111948708B (en) | Seismic wave field forward modeling method for dipping in undulating surface of boundary | |
| CN109459787B (en) | coal mine underground structure imaging method and system based on seismic channel wave full-waveform inversion | |
| Muratov et al. | Grid-characteristic method on unstructured tetrahedral meshes | |
| CN114139335A (en) | Viscous sound wave simulation method based on single relaxation time lattice Boltzmann model | |
| CN117150835A (en) | Method for simulating seismic wave value of viscous-acoustic anisotropic medium | |
| Koren et al. | Eigenrays in 3D heterogeneous anisotropic media, Part I: Kinematics | |
| CN109239776B (en) | Seismic wave propagation forward modeling method and device | |
| NO20190489A1 (en) | Seismic modeling | |
| CN111257930B (en) | Visco-elastic anisotropic double-phase medium area variable grid solving operator | |
| US20160034612A1 (en) | Re-ordered Interpolation and Convolution for Faster Staggered-Grid Processing | |
| Favorskaya et al. | Wave processes modelling in geophysics | |
| CN106353801A (en) | Simulation method and device for 3D Laplace domain acoustic wave equation value | |
| Wencai | Reflection Seismology: theory, data processing and interpretation | |
| Lehmann et al. | Seismic hazard analysis with a Factorized Fourier Neural Operator (F-FNO) surrogate model enhanced by transfer learning | |
| CN107807392A (en) | A kind of piecemeal space-time of adaptive anti-frequency dispersion is double to become reverse-time migration method | |
| Jing et al. | An optimized time-space-domain finite difference method with piecewise constant interpolation coefficients for scalar wave propagation | |
| Zakian et al. | Spectral finite element simulation of seismic wave propagation and fault dislocation in elastic media | |
| CN111665550A (en) | Underground medium density information inversion method | |
| Calandra et al. | Recent advances in numerical methods for solving the wave equation in the context of seismic depth imaging |
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 |