CN109725346B - Time-space high-order VTI rectangular grid finite difference method and device - Google Patents

Time-space high-order VTI rectangular grid finite difference method and device Download PDF

Info

Publication number
CN109725346B
CN109725346B CN201811542112.4A CN201811542112A CN109725346B CN 109725346 B CN109725346 B CN 109725346B CN 201811542112 A CN201811542112 A CN 201811542112A CN 109725346 B CN109725346 B CN 109725346B
Authority
CN
China
Prior art keywords
order
equation
wave
finite difference
spatial
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.)
Active
Application number
CN201811542112.4A
Other languages
Chinese (zh)
Other versions
CN109725346A (en
Inventor
刘洋
徐世刚
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Petroleum Beijing
China National Petroleum Corp
Original Assignee
China University of Petroleum Beijing
China National Petroleum Corp
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China University of Petroleum Beijing, China National Petroleum Corp filed Critical China University of Petroleum Beijing
Priority to CN201811542112.4A priority Critical patent/CN109725346B/en
Publication of CN109725346A publication Critical patent/CN109725346A/en
Application granted granted Critical
Publication of CN109725346B publication Critical patent/CN109725346B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The embodiment of the specification provides a time-space high-order VTI rectangular grid finite difference method and device. The method comprises the following steps: converting a VTI elastic wave equation into a matrix wave equation; expanding the matrix wave equation to a time fourth-order precision form to obtain a discrete wave equation; calculating a spatial partial derivative based on the discrete wave equation; and solving an optimized finite difference coefficient based on the spatial partial derivative. By the method, the time dispersion precision can be effectively improved, and the wave field extrapolation calculation efficiency is improved while the precision is ensured.

Description

Time-space high-order VTI rectangular grid finite difference method and device
Technical Field
The embodiment of the specification relates to the technical field of wave field extrapolation, in particular to a time-space high-order VTI rectangular grid finite difference method and device.
Background
Compared with an acoustic wave field, an elastic wave field can provide richer and comprehensive underground geological structures, so that the offset imaging and waveform inversion based on an elastic wave equation become leading edges and hot spots of seismic exploration.
The high-precision and high-efficiency wave field extrapolation method is one of the key factors of elastic wave migration imaging and waveform inversion. The finite difference method, especially the staggered grid algorithm, is widely applied to elastic wave number value simulation due to the advantages of simple implementation, small calculation amount, easy parallelism and the like. When the finite difference method is used for carrying out wave field numerical simulation, the main reason influencing the finite difference discrete precision of the elastic wave equation is the numerical dispersion of spatial and time terms.
At present, high-precision finite difference methods for elastic wave equations mainly have four types: the first type is space high-order, time second-order staggered grid discrete, and corresponding differential coefficients are solved through Taylor series expansion, and the method has the advantages that the realization is simple, but the method is easily interfered by space and time frequency dispersion; the second type is that a linear or nonlinear optimization algorithm is adopted to optimize the difference coefficient, so that the spatial precision is obviously improved, but the time dispersion still exists and the stability is reduced; the third type is that a higher-precision difference coefficient is obtained through the frequency dispersion relation of a time-space domain, so that higher time and space simulation precision can be obtained, and more complex and time-consuming nonlinear optimization is realized; the fourth type is that more grid nodes are added near the paraxial of the traditional cross differential template, so that the time high-order precision is realized, the time and space can reach any even-order precision, but the realization is more complex, and the calculation amount is increased due to the introduction of additional grid nodes for calculation.
Because the spatial derivative and the anisotropic parameter are mutually decoupled, the spatial high-order discrete format based on Taylor expansion and optimized difference coefficient can be directly popularized to the numerical simulation of VTI elastic waves. In contrast, the time-space domain method and the time-high order method by adding paraxial grid points are difficult to be extended to the VTI seismic wave equation numerical simulation because the time-domain dispersion is a nonlinear function in which the difference coefficient and the anisotropy parameter are coupled with each other. Therefore, the current wave equation in the VTI medium is still discretized by using a time second order difference format, so that the conventional method suffers from serious time value dispersion and even instability when the time step is large.
Disclosure of Invention
The purpose of the embodiments of the present specification is to provide a finite difference method for a time-space high-order VTI rectangular grid, so as to realize wavefield extrapolation by using the finite difference method more efficiently on the premise of higher precision.
In order to solve the above technical problem, an embodiment of the present application provides a time-space high-order VTI rectangular grid finite difference method and apparatus, which are implemented as follows:
a time-space high-order VTI rectangular grid finite difference method comprises the following steps:
converting a VTI elastic wave equation into a matrix wave equation;
expanding the matrix wave equation to a time fourth-order precision form to obtain a discrete wave equation;
calculating a spatial partial derivative based on the discrete wave equation;
and solving an optimized finite difference coefficient based on the spatial partial derivative.
A time-space high-order VTI rectangular grid finite difference device, comprising:
the matrix wave equation conversion module is used for converting the VTI elastic wave equation into a matrix wave equation;
the discrete wave equation acquisition module is used for expanding the matrix wave equation to a time fourth-order precision form to obtain a discrete wave equation;
the spatial partial derivative calculating module is used for calculating a spatial partial derivative based on the discrete wave equation;
and the finite difference coefficient solving module is used for solving the optimized finite difference coefficient based on the spatial partial derivative.
According to the technical scheme provided by the embodiment of the specification, when the VTI elastic wave equation in the medium is continuously converted to obtain a form beneficial to discrete expansion, the equation is expanded to a discrete wave equation form with time fourth-order precision, then the spatial partial derivative is obtained according to the discrete wave equation, and the finite difference coefficient is further obtained, so that the time-space high-order finite difference in the medium is realized. By the method, wave field extrapolation can be more reliably provided for elastic wave migration imaging and waveform inversion, and geological exploration and other operations can be facilitated.
Drawings
In order to more clearly illustrate the embodiments of the present specification or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly described below, it is obvious that the drawings in the following description are only some embodiments described in the specification, and other drawings can be obtained by those skilled in the art without creative efforts.
FIG. 1 is a flowchart illustrating a method for finite difference of a temporal-spatial high-order VTI rectangular grid according to an embodiment of the present disclosure;
FIG. 2 is a block diagram of a time-space high-order VTI rectangular grid finite difference device according to an embodiment of the present disclosure;
FIG. 3A is a diagram of a snapshot of a wave field at 1.4s in a time step of 2ms in the conventional method in this specification;
FIG. 3B is a schematic diagram of a 2.1s wave field snapshot of the conventional method in 2ms time step in this specification;
FIG. 3C is a diagram of a wave field snapshot at 1.4s in a time step of 1ms in the conventional method in this specification;
FIG. 3D is a schematic diagram of a 2.1s wave field snapshot of the conventional method in a time step of 1 ms;
FIG. 3E is a schematic diagram of a snapshot of the wavefield at 1.4s in one embodiment of the present description at a time step of 2 ms;
FIG. 3F is a schematic diagram of a 2.1s wavefield snapshot in one embodiment of the present description at a time step of 2 ms;
FIG. 4A is a schematic view of a parametric model of compressional velocity in one embodiment of the present disclosure;
FIG. 4B is a graphical illustration of a parametric model of density in one embodiment of the present disclosure;
FIG. 4C is a graphical illustration of a parametric model of the anisotropic parameters in one embodiment of the present disclosure;
FIG. 4D is a graphical illustration of a parametric model of anisotropic parameters in one embodiment of the present disclosure;
FIG. 5A is a schematic diagram of a wave field snapshot at 3.6s under a time step of 1.2ms in the conventional method in this specification;
FIG. 5B is a diagram of a snapshot of a wave field at 3.6s in a time step of 0.6ms in the conventional method in this specification;
FIG. 5C is a schematic diagram of a 3.6s wavefield snapshot at a time step of 1.2ms in one embodiment of the present description;
FIG. 6A is a schematic diagram of an elastic wave forward modeling seismic record under a time step of 1.2ms in the conventional method in this specification;
FIG. 6B is a schematic diagram of an elastic wave forward modeling seismic record under a time step of 0.6ms in the conventional method in this specification;
FIG. 6C is a schematic diagram of an elastic wave forward modeling seismic recording at a time step of 1.2ms in one embodiment of the present description.
Detailed Description
The technical solutions in the embodiments of the present disclosure will be clearly and completely described below with reference to the drawings in the embodiments of the present disclosure, and it is obvious that the described embodiments are only a part of the embodiments of the present disclosure, and not all of the embodiments. All other embodiments obtained by a person of ordinary skill in the art based on the embodiments in the present specification without any creative effort shall fall within the protection scope of the present specification.
An embodiment of a time-space high-order VTI rectangular grid finite difference method according to the present application is described below with reference to fig. 1. The execution subject of the method is the server. The time-space high-order VTI rectangular grid finite difference method comprises the following specific implementation steps:
s100: and converting the VTI elastic wave equation into a matrix wave equation.
With the development of address exploration technology, numerical simulation becomes an important means for researching the kinematic and dynamic characteristics of seismic wave propagation in the earth medium. In many methods of numerical simulation, a wave equation numerical solution is used for discrete gridding of a calculation area, and a wave equation for describing seismic wave propagation is obtained through numerical solution, so that the method is a better implementation method. After the approximate solution with higher precision of the wave equation is obtained by solving, the method can be used for realizing the simulation of the propagation of the seismic wave. The finite difference method is a quick and effective method for solving the wave equation, and realizes the wave field extrapolation by transforming the wave differential equation in the medium into a finite difference form for solving through the discrete gridding of the elastic medium model.
An isotropic medium (isotropic medium) refers to a medium in which the properties of matter (physical, chemical, etc.) are the same in different directions. In isotropic media, the propagation speed of light is independent of the polarization direction of the light, i.e. light has only one refractive index. Air, water, optical glass, ordinary glass, etc. are isotropic media. VTI media, i.e., Transversely isotropic (isotropic) media, is one of the isotropic media that possesses the same physical properties only in the transverse direction. In the VTI medium, the elastic wave can be expressed by an equation including information such as the velocity of the elastic wave, the stress, and the state parameter of the medium. According to the elastic wave equation, through a series of conversion, the time-space high-order finite difference based on the rectangular staggered grid can be realized.
The following is a detailed analysis thereof. For elastic waves in VTI media, there is an equation of the form:
Figure BDA0001908401990000041
where ρ is the density of the medium, (v)x,vz) Is the velocity component (τ)xx,τzz,τxz) In order to be a stress component,
Figure BDA0001908401990000042
and
Figure BDA0001908401990000043
is a VTI elastic modulus parameter, wherein vpAnd vsLongitudinal and transverse wave velocities, respectively, and ε and δ are Thomsen anisotropy parameters. The above equations mainly describe the properties of speed and stress, and specifically map the elastic waves to corresponding state parameters, so that the subsequent solving process is more intuitive and the subsequent process is better realized.
In order to facilitate the later deduction of the velocity component and the stress component, and to facilitate the later solution, a first matrix p is provided for the velocity component (v)x,vz) And stress component (tau)xx,τzz,τxz) Integration is carried out, i.e. let p ═ vx,vz,τxx,τzz,τxz)T. Therefore, the elastic wave equation in the VTI medium can be expressed by using the first matrix, and the subsequent solving process is facilitated.
In the case where the first matrix has been set, the elastic wave equation in the VTI medium can be expressed more conveniently. Based on the set first matrix p ═ (v)x,vz,τxx,τzz,τxz)TEquation of the elastic wave
Figure BDA0001908401990000051
Conversion into a matrix wave equation
Figure BDA0001908401990000052
Wherein, w is set as the second matrix,
Figure BDA0001908401990000053
ρ is the density of the medium, (v)x,vz) Is the velocity component (τ)xx,τzz,τxz) In order to be a stress component,
Figure BDA0001908401990000054
and
Figure BDA0001908401990000055
is the VTI elastic coefficient parameter, vpAnd vsLongitudinal and transverse wave velocities, respectively, and ε and δ are Thomsen anisotropy parameters.
S200: and expanding the matrix wave equation to a time fourth-order precision form to obtain a discrete wave equation.
To achieve a discrete representation of the elastic wave in the medium, the elastic wave equation for the first matrix is expanded into a discrete form according to a rectangular interleaved network, so that extrapolation of the elastic wave field can be achieved later. In a specific example, the matrix wave equation is expanded by adopting Taylor series at a time node t + delta t to reach the time fourth-order discrete precision, so that a discrete wave equation is obtained
Figure BDA0001908401990000056
In the formula (I), the compound is shown in the specification,
Figure BDA0001908401990000061
wherein the content of the first and second substances,
Figure BDA0001908401990000062
Figure BDA0001908401990000063
Figure BDA0001908401990000064
Figure BDA0001908401990000065
Figure BDA0001908401990000066
Figure BDA0001908401990000067
Figure BDA0001908401990000068
Figure BDA0001908401990000069
p is the density of the mixture,
Figure BDA00019084019900000610
Figure BDA00019084019900000611
and
Figure BDA00019084019900000612
as a parameter of the modulus of elasticity, vpIs the velocity of longitudinal wave, vsFor shear wave velocity,. epsilon.and.delta.are Thomsen anisotropy parameters.
Based on the discrete wave equation, the wave field extrapolation can be effectively realized according to the state parameters of the known region, and the data simulation is carried out on the current target measurement region.
S300: and solving a spatial partial derivative based on the discrete wave equation.
To further refine the discrete wave equation, the spatial partial derivatives thereof still need to be calculated.
For the first spatial partial derivative in the second matrix
Figure BDA0001908401990000071
And
Figure BDA0001908401990000072
the discrete format of the first order spatial partial derivatives can be expressed as follows:
Figure BDA0001908401990000073
Figure BDA0001908401990000074
wherein h isxAnd hzRepresenting the spatial separation in the x and z directions, N1To first order spatial operator length, cnIs a first order finite difference coefficient.
For the third order spatial partial derivatives in the third matrix
Figure BDA0001908401990000075
And
Figure BDA0001908401990000076
the discrete format of the third order spatial partial derivatives can be expressed as follows:
Figure BDA0001908401990000077
Figure BDA0001908401990000078
in the formula, hxAnd hzRepresenting the spatial separation in the x and z directions, N2As length of spatial operator, cn' is a third order finite difference coefficient.
For mixed third-order spatial partial derivatives in the third matrix
Figure BDA0001908401990000079
And
Figure BDA00019084019900000710
the discrete format of the mixed third-order spatial partial derivatives can be represented as follows:
Figure BDA00019084019900000711
Figure BDA00019084019900000712
the third-order spatial partial derivative obtained by the calculation
Figure BDA00019084019900000713
And
Figure BDA00019084019900000714
and (4) the calculation result is brought back to the third matrix, and the third matrix can be calculated.
The number of points used for calculating the third-order spatial partial derivative is consistent with the number of points used for calculating the first-order spatial partial derivative. And higher simulation precision can be obtained by adopting more points, namely, the lengths of the space operators are equal. However, similar to calculating the higher order dispersion of the first order spatial partial derivatives, the use of larger spatial operator lengths significantly increases the amount of additional computation when calculating the third order spatial partial derivatives. Therefore, in actual calculation, in order to balance the calculation efficiency and accuracy, the formula N is adopted2＝int(N1And/2) calculating a third-order spatial partial derivative, and improving the calculation efficiency on the premise of ensuring the precision.
S400: and solving an optimized finite difference coefficient based on the spatial partial derivative.
The optimized finite difference coefficient comprises a first order difference coefficient and a third order difference coefficient.
For first order difference coefficient c in first order spatial partial derivativen(n＝1,2,...,N1) Can be obtained by solving a system of linear equations
Figure BDA0001908401990000081
And (4) obtaining. After the first order difference coefficient is obtained through calculation, the first order difference coefficient is brought into a discrete format of a first order spatial partial derivative, and then the first order spatial partial derivative is put into
Figure BDA0001908401990000082
And
Figure BDA0001908401990000083
and substituting the second matrix to obtain the second matrix which can be specifically solved, so as to facilitate the later calculation.
For third-order difference coefficient c 'in third-order spatial partial derivative'n(n＝1,2,...,N2) Can be obtained by solving a system of linear equations
Figure BDA0001908401990000084
And (4) obtaining. Solving the equation, and substituting the obtained third-order difference coefficient into a third-order partial derivative
Figure BDA0001908401990000085
And
Figure BDA0001908401990000086
i.e. for subsequent evaluation of the third matrix.
The difference coefficient obtained by taylor expansion can obtain high precision in the wavelet number range, but the dispersion error gradually increases with the increase of the wave number. In order to further improve the simulation accuracy, the finite difference coefficients can be optimized by using a least square method. And combining a plane wave solution equation and the spatial partial derivative to obtain an intermediate equation, minimizing errors at two ends of the intermediate equation, combining a least square method theory to obtain an optimized equation, and then solving the optimized equation to finally obtain an optimized finite difference coefficient.
A specific embodiment for finding the optimized third-order staggered grid difference coefficient is described below. Solving plane waves
Figure BDA0001908401990000087
Into the equation
Figure BDA0001908401990000088
In the specification, can obtain after simplificationTo
Figure BDA0001908401990000089
Wherein h isxAnd hzRepresenting the spatial separation in the x and z directions, N2For spatial operator length, κn(β)＝2sin[(n-1/2)β]，f(β)＝β3，β＝kxhx，β∈[0,π]，kx、kzIn x and z directions. Minimization equation
Figure BDA0001908401990000091
And establishing an objective function according to errors at two ends as follows:
Figure BDA0001908401990000092
based on the least square theory, obtaining an optimization equation of
Figure BDA0001908401990000093
And solving the optimization equation by using a least square method again to obtain an optimized finite difference coefficient of the third-order spatial partial derivative.
Based on the method described in the above embodiment, the first-order finite difference coefficient can be obtained in the same manner, and details are not repeated again. By utilizing the optimized finite difference coefficient, an experimental result can be better solved, and an elastic wave field simulation result can be more accurately obtained.
The effectiveness and accuracy of the application are proved by applying the method to the specific embodiment and combining the comparison between the effect realized by the traditional method and the effect realized by the traditional method.
In one embodiment, as shown in fig. 3A to 3F, for the homogeneous medium model, the dimension of the model is set to 800 × 800, the grid spacing along the x direction is 15m, and the grid spacing along the z direction is 10 m. The speeds of P wave and S wave are 3100m/S and 1750 m/S respectively, and the medium density is 2300kg/m3The values of the anisotropy parameters were ∈ ═ 0.15 and δ ═ 0.15, respectively. The Rake wavelet with the main frequency of 27Hz is placed in the center of the model to generate vibration. Finite difference method for conventional staggered gridMethod, using spatial operator with length N 18. For the improved time fourth-order staggered grid finite difference method, the length of a space operator is N 18 and N24. Fig. 3A is a schematic view of a wave field snapshot calculated by a conventional method when a time step of 2ms is 1.4s, fig. 3C is a schematic view of a wave field snapshot calculated by a conventional method when a time step of 1ms is 1.4s, and fig. 3E is a schematic view of a wave field snapshot calculated by a method according to the present application when a time step of 2ms is 1.4 s. By comparing the three schematic diagrams, the accuracy of the traditional method is obviously reduced when a larger time step is adopted, and the schematic diagram obtained by the method is more accurate when the same time step is adopted, so that the experimental effect obtained when the traditional method is adopted for a smaller time step is achieved. Under the condition of realizing the same accuracy, the method can adopt larger time step length for calculation, thereby undoubtedly reducing the calculation steps and improving the calculation efficiency. Fig. 3B is a schematic view of a wave field snapshot calculated at 2.1s when a time step length of 2ms is adopted in the conventional method, fig. 3D is a schematic view of a wave field snapshot calculated at 2.1s when a time step length of 1ms is adopted in the conventional method, and fig. 3F is a schematic view of a wave field snapshot at 2.1s when a time step length of 2ms is adopted in the method according to the present application. The accuracy and effectiveness of the application are further proved through the three schematic diagrams.
In order to further verify the effect obtained by the application, a two-dimensional Hess VTI model is adopted to carry out a numerical simulation test. Fig. 4A is a schematic diagram of a parametric model of a longitudinal wave velocity in the Hess VTI model, fig. 4B is a schematic diagram of a parametric model of a density in the Hess VTI model, fig. 4C is a schematic diagram of a parametric model of a Thomsen anisotropic parameter epsilon in the Hess VTI model, and fig. 4D is a schematic diagram of a parametric model of a Thomsen anisotropic parameter delta in the Hess VTI model. The size of the model is 22000m × 9000m, and the model is resampled to 1808 × 750 grid cells using a 12m × 12m grid. The grid step along the x-direction is 10m and the grid step along the z-direction is 15 m. The velocity ratio of longitudinal and transverse waves is 1.7:1, the seismic source is a period of Rake wavelets, the dominant frequency is 12Hz, and the location of the seismic source is (11.28km, 0.12km)The recording time period was 6 s. For the traditional staggered grid finite difference method, the length of a space operator is N 112. For the improved time fourth-order staggered grid finite difference method, the length of a space operator is respectively N 112 and N 26. Fig. 5A is a schematic diagram of an elastic wave forward modeling simulated wave field snapshot calculated at 3.6s when the time step of 1.2ms is adopted in the conventional method, fig. 5B is a schematic diagram of an elastic wave forward modeling simulated wave field snapshot calculated at 3.6s when the time step of 0.6ms is adopted in the conventional method, and fig. 5C is a schematic diagram of an elastic wave forward modeling simulated wave field snapshot calculated at 3.6s when the time step of 1.2ms is adopted in the embodiment of the present application. It can be seen from comparison among the three diagrams that the conventional method can improve the accuracy of the method by shortening the time step, but cannot ensure the accuracy of the wave field calculated by the method when the time step is large, and is easily interfered by time dispersion. However, according to the effect shown by the embodiment obtained by the method of the present application, it can be seen that when the same time step length is adopted, the wave field snapshot obtained by the embodiment of the present application can have the accuracy obtained according to the conventional method with a shorter time, that is, under the condition of the same accuracy limitation, the method adopted by the present application can adopt a smaller time step length, thereby reducing the amount of calculation and improving the calculation efficiency. Further, fig. 6A is a schematic diagram of an elastic wave forward modeling simulated seismic record obtained by the conventional method by using a time step of 1.2ms, fig. 6B is a schematic diagram of an elastic wave forward modeling simulated seismic record obtained by the conventional method by using a time step of 0.6ms, and fig. 6C is a schematic diagram of an elastic wave forward modeling simulated seismic record obtained by the method of the present application by using a time step of 1.2 ms. The accuracy and reliability of the method of the present application are further illustrated by the comparison between the three diagrams.
An embodiment of a time-space high-order VTI rectangular grid finite difference device according to the present application is described below. As shown in fig. 2, the time-space high-order VTI rectangular grid finite difference device includes:
a matrix wave equation conversion module 210, configured to convert the VTI elastic wave equation into a matrix wave equation;
the discrete wave equation obtaining module 220 is configured to expand the matrix wave equation to a time fourth-order precision form to obtain a discrete wave equation;
a spatial partial derivative calculation module 230, configured to calculate a spatial partial derivative based on the discrete wave equation;
and an optimized finite difference coefficient solving module 240, configured to solve the optimized finite difference coefficient based on the spatial partial derivative.
The matrix wave equation conversion module 210 includes:
a parameter matrix setting subunit 211 for setting the parameter matrix p ═ (v) based on the velocity component and the stress component in the elastic wave equationx,vz,τxx,τzz,τxz)TWherein v isx,vzIs the velocity component of the elastic wave, τxx,τzz,τxzIs the stress component of the elastic wave;
a matrix wave equation obtaining subunit 212 for obtaining the elastic wave equation
Figure BDA0001908401990000111
Conversion into a matrix wave equation
Figure BDA0001908401990000112
Wherein w is a second matrix,
Figure BDA0001908401990000113
p is the density of the mixture,
Figure BDA0001908401990000114
and
Figure BDA0001908401990000115
as a parameter of the modulus of elasticity, vpIs the velocity of longitudinal wave, vsFor shear wave velocity,. epsilon.and.delta.are Thomsen anisotropy parameters.
The discrete wave equation obtaining module 220 includes:
discrete wave equationAn expansion subunit 221, configured to expand the matrix wave equation by using Taylor series
Figure BDA0001908401990000116
Discrete wave equation expanded to time fourth order precision
Figure BDA0001908401990000117
In the formula (I), the compound is shown in the specification,
Figure BDA0001908401990000121
wherein the content of the first and second substances,
Figure BDA0001908401990000122
Figure BDA0001908401990000123
Figure BDA0001908401990000124
Figure BDA0001908401990000125
Figure BDA0001908401990000126
Figure BDA0001908401990000127
Figure BDA0001908401990000128
Figure BDA0001908401990000129
p is the density of the mixture,
Figure BDA00019084019900001210
Figure BDA00019084019900001211
and
Figure BDA00019084019900001212
as a parameter of the modulus of elasticity, vpIs the velocity of longitudinal wave, vsFor shear wave velocity,. epsilon.and.delta.are Thomsen anisotropy parameters.
The spatial partial derivative calculating module 230 includes:
a first order spatial partial derivative solving subunit 231 for utilizing the formula
Figure BDA00019084019900001213
Figure BDA00019084019900001214
First order spatial partial derivative is obtained
Figure BDA00019084019900001215
And
Figure BDA00019084019900001216
in the formula, hxAnd hzSpatial separation along the x and z directions, respectively, N1To first order spatial operator length, cnIs a first order finite difference coefficient;
a third order spatial partial derivative solving subunit 232 for utilizing the formula
Figure BDA0001908401990000131
Figure BDA0001908401990000132
Calculating the third-order spatial partial derivative
Figure BDA0001908401990000133
And
Figure BDA0001908401990000134
in the formula, hxAnd hzRepresenting the spatial separation in the x and z directions, N2For three spatial operator lengths, cn' is a third order finite difference coefficient;
hybrid third-order spatial partial derivative solving subunit 233 for utilizing the formula
Figure BDA0001908401990000135
Figure BDA0001908401990000136
Solving for the mixed third-order spatial partial derivative, where hxAnd hzRepresenting the spatial separation in the x and z directions, N2Is the third order spatial operator length.
The optimized finite difference coefficient comprises a first-order finite difference coefficient and a third-order finite difference coefficient; the optimized finite difference coefficient solving module 240 includes:
a first order finite difference coefficient solving subunit 241 for solving a first order equation set by using Taylor expansion based on the first order spatial partial derivative
Figure BDA0001908401990000137
To obtain a first order finite difference coefficient cn(n＝1,2,...,N1) In the formula, N1Is the spatial operator length;
a third-order finite difference coefficient solving subunit 242, configured to solve the third-order equation set by using taylor expansion based on the third-order spatial partial derivative
Figure BDA0001908401990000138
Obtaining a third-order finite difference coefficient c'n(n＝1,2,...,N2) In the formula, N2Is the spatial operator length.
Meanwhile, the optimizing finite difference coefficient solving module 240 may further include:
an intermediate equation solving subunit 243, configured to combine the plane wave solution with the spatial partial derivative to obtain an intermediate equation;
an optimization equation solving subunit 244, configured to minimize errors at two ends of the intermediate equation, and obtain an optimization equation by combining a least square method;
and the finite difference coefficient calculating subunit 245 is configured to solve the optimization equation to obtain the optimized finite difference coefficient.
The finite difference coefficients as described in the finite difference coefficient calculation subunit 245 include at least one of:
a first order finite difference coefficient;
third order finite difference coefficients.
In the 90 s of the 20 th century, improvements in a technology could clearly distinguish between improvements in hardware (e.g., improvements in circuit structures such as diodes, transistors, switches, etc.) and improvements in software (improvements in process flow). However, as technology advances, many of today's process flow improvements have been seen as direct improvements in hardware circuit architecture. Designers almost always obtain the corresponding hardware circuit structure by programming an improved method flow into the hardware circuit. Thus, it cannot be said that an improvement in the process flow cannot be realized by hardware physical modules. For example, a Programmable Logic Device (PLD), such as a Field Programmable Gate Array (FPGA), is an integrated circuit whose Logic functions are determined by programming the Device by a user. A digital system is "integrated" on a PLD by the designer's own programming without requiring the chip manufacturer to design and fabricate a dedicated integrated circuit chip 2. Furthermore, nowadays, instead of manually making an Integrated Circuit chip, such Programming is often implemented by "logic compiler" software, which is similar to a software compiler used in program development and writing, but the original code before compiling is also written by a specific Programming Language, which is called Hardware Description Language (HDL), and HDL is not only one but many, such as abel (advanced Boolean Expression Language), ahdl (alternate Language Description Language), traffic, pl (core unified Programming Language), HDCal, JHDL (Java Hardware Description Language), langue, Lola, HDL, laspam, hardbyscript Description Language (vhr Description Language), and the like, which are currently used by Hardware compiler-software (Hardware Description Language-software). It will also be apparent to those skilled in the art that hardware circuitry that implements the logical method flows can be readily obtained by merely slightly programming the method flows into an integrated circuit using the hardware description languages described above.
The systems, devices, modules or units illustrated in the above embodiments may be implemented by a computer chip or an entity, or by a product with certain functions. One typical implementation device is a computer. In particular, the computer may be, for example, a personal computer, a laptop computer, a cellular telephone, a camera phone, a smartphone, a personal digital assistant, a media player, a navigation device, an email device, a game console, a tablet computer, a wearable device, or a combination of any of these devices.
From the above description of the embodiments, it is clear to those skilled in the art that the present specification can be implemented by software plus a necessary general hardware platform. Based on such understanding, the technical solutions of the present specification may be essentially or partially implemented in the form of software products, which may be stored in a storage medium, such as ROM/RAM, magnetic disk, optical disk, etc., and include instructions for causing a computer device (which may be a personal computer, a server, or a network device, etc.) to execute the methods described in the embodiments or some parts of the embodiments of the present specification.
The embodiments in the present specification are described in a progressive manner, and the same and similar parts among the embodiments are referred to each other, and each embodiment focuses on the differences from the other embodiments. In particular, for the system embodiment, since it is substantially similar to the method embodiment, the description is simple, and for the relevant points, reference may be made to the partial description of the method embodiment.
The description is operational with numerous general purpose or special purpose computing system environments or configurations. For example: personal computers, server computers, hand-held or portable devices, tablet-type devices, multiprocessor systems, microprocessor-based systems, set top boxes, programmable consumer electronics, network PCs, minicomputers, mainframe computers, distributed computing environments that include any of the above systems or devices, and the like.
This description may be described in the general context of computer-executable instructions, such as program modules, being executed by a computer. Generally, program modules include routines, programs, objects, components, data structures, etc. that perform particular tasks or implement particular abstract data types. The specification may also be practiced in distributed computing environments where tasks are performed by remote processing devices that are linked through a communications network. In a distributed computing environment, program modules may be located in both local and remote computer storage media including memory storage devices.
While the specification has been described with examples, those skilled in the art will appreciate that there are numerous variations and permutations of the specification that do not depart from the spirit of the specification, and it is intended that the appended claims include such variations and modifications that do not depart from the spirit of the specification.

Claims (8)

1. A time-space high-order VTI rectangular grid finite difference method is characterized by comprising the following steps:
converting a VTI elastic wave equation into a matrix wave equation;
expanding the matrix wave equation to a time fourth-order precision rectangular staggered grid differential form to obtain a discrete wave equation;
calculating a spatial partial derivative based on the discrete wave equation; the spatial partial derivatives include the following three types: a first order spatial partial derivative, a third order spatial partial derivative, and a mixed third order spatial partial derivative; wherein, include: using formulas
Figure FDA0002967599420000011
Figure FDA0002967599420000012
First order spatial partial derivative is obtained
Figure FDA0002967599420000013
And
Figure FDA0002967599420000014
wherein p (i, j) is a plane wave component in a discrete wave equation, i, j are grid component indexes of the plane wave component along different spatial directions, hxAnd hzSpatial separation along the x and z directions, respectively, N1To first order spatial operator length, cnIs a first order finite difference coefficient; or, using a formula
Figure FDA0002967599420000015
Figure FDA0002967599420000016
Calculating the third-order spatial partial derivative
Figure FDA0002967599420000017
And
Figure FDA0002967599420000018
in the formula, hxAnd hzRepresenting the spatial separation in the x and z directions, N2For three spatial operator lengths, cn' is a third order finite difference coefficient; using formulas
Figure FDA0002967599420000019
Figure FDA00029675994200000110
Solving for the mixed third-order spatial partial derivative, where hxAnd hzRepresenting the spatial separation in the x and z directions, N2Is the length of the three-order spatial operator;
based on the spatial partial derivative, calculating an optimized finite difference coefficient; the optimized finite difference coefficient comprises a first-order finite difference coefficient and a third-order finite difference coefficient; wherein, include: solving a first order equation set by using Taylor expansion based on a first order spatial partial derivative
Figure FDA0002967599420000021
To obtain a first order finite difference coefficient cnWherein c isnN in (1, 2., N)1In the formula, N1Is the first order spatial operator length; or, based on the third-order spatial partial derivative, the third-order equation set is solved by using Taylor expansion
Figure FDA0002967599420000022
Obtaining the third order finite difference coefficient cn', wherein, cnN of' 1,22In the formula, N2Is the third order spatial operator length.
2. The method of claim 1, wherein converting the VTI elastic wave equation to a matrix wave equation comprises:
setting a parameter matrix p ═ (v) based on the velocity component and the stress component in the VTI elastic wave equationx,vz,τxx,τzz,τxz)TWherein v isx,vzIs the velocity component of the elastic wave, τxx,τzz,τxzIs the stress component of the elastic wave;
equation of elastic wave
Figure FDA0002967599420000023
Conversion into a matrix wave equation
Figure FDA0002967599420000024
Wherein w is a second matrix,
Figure FDA0002967599420000025
p is the density of the mixture,
Figure FDA0002967599420000026
Figure FDA0002967599420000027
and
Figure FDA0002967599420000028
as a parameter of the modulus of elasticity, vpIs the velocity of longitudinal wave, vsFor shear wave velocity,. epsilon.and.delta.are Thomsen anisotropy parameters.
3. The method of claim 1, wherein the expanding the matrix wave equation into a time fourth order precision rectangular staggered grid differential form yields a discrete wave equation, comprising:
using Taylor series expansion to transform the matrix wave equation
Figure FDA0002967599420000031
Discrete wave equation expanded to time fourth order precision
Figure FDA0002967599420000032
In the formula (I), the compound is shown in the specification,
Figure FDA0002967599420000033
wherein the content of the first and second substances,
Figure FDA0002967599420000034
Figure FDA0002967599420000035
Figure FDA0002967599420000036
Figure FDA0002967599420000037
Figure FDA0002967599420000038
Figure FDA0002967599420000039
Figure FDA00029675994200000310
Figure FDA00029675994200000311
Figure FDA00029675994200000312
p is the density of the mixture,
Figure FDA00029675994200000313
Figure FDA0002967599420000041
and
Figure FDA0002967599420000042
as a parameter of the modulus of elasticity, vpIs the velocity of longitudinal wave, vsFor shear wave velocity,. epsilon.and. delta.are Thomsen anisotropy parameters, t is a time node, and. delta.t is a preset time interval.
4. The method of claim 1, wherein the deriving optimized finite difference coefficients based on the spatial partial derivatives comprises:
combining the plane wave solution with the spatial partial derivative to obtain an intermediate equation;
minimizing errors at two ends of the intermediate equation, and combining a least square method to obtain an optimized equation;
and solving the optimization equation to obtain an optimized finite difference coefficient.
5. A time-space high-order VTI rectangular grid finite difference device, comprising;
the matrix wave equation conversion module is used for converting the VTI elastic wave equation into a matrix wave equation;
the discrete wave equation acquisition module is used for expanding the matrix wave equation to a time fourth-order precision rectangular staggered grid differential form to obtain a discrete wave equation;
the spatial partial derivative calculating module is used for calculating a spatial partial derivative based on the discrete wave equation; the spatial partial derivatives include the following three types: a first order spatial partial derivative, a third order spatial partial derivative, and a mixed third order spatial partial derivative; wherein, include: using formulas
Figure FDA0002967599420000043
Figure FDA0002967599420000044
First order spatial partial derivative is obtained
Figure FDA0002967599420000045
And
Figure FDA0002967599420000046
wherein p (i, j) is a plane wave component in a discrete wave equation, i, j are grid component indexes of the plane wave component along different spatial directions, hxAnd hzSpatial separation along the x and z directions, respectively, N1To first order spatial operator length, cnIs a first order finite difference coefficient; or, using a formula
Figure FDA0002967599420000047
Figure FDA0002967599420000048
Calculating the third-order spatial partial derivative
Figure FDA0002967599420000049
And
Figure FDA00029675994200000410
in the formula, hxAnd hzRepresenting the spatial separation in the x and z directions, N2For three spatial operator lengths, cn' is a third order finite difference coefficient; using formulas
Figure FDA00029675994200000411
Figure FDA0002967599420000051
Solving for the mixed third-order spatial partial derivative, where hxAnd hzRepresenting the spatial separation in the x and z directions, N2Is the length of the three-order spatial operator;
the optimized finite difference coefficient calculating module is used for calculating an optimized finite difference coefficient based on the spatial partial derivative; the optimized finite difference coefficient comprises a first-order finite difference coefficient and a third-order finite difference coefficient; wherein, include: solving a first order equation set by using Taylor expansion based on a first order spatial partial derivative
Figure FDA0002967599420000052
To obtain a first order finite difference coefficient cnWherein c isnN in (1, 2., N)1In the formula, N1Is the first order spatial operator length; or, based on the third-order spatial partial derivative, the third-order equation set is solved by using Taylor expansion
Figure FDA0002967599420000053
Obtaining the third order finite difference coefficient cn', wherein, cnN of' 1,22In the formula, N2Is the third order spatial operator length.
6. The apparatus of claim 5, wherein the matrix wave equation transformation module comprises:
a parameter matrix setting subunit for setting a parameter matrix p ═ (v) based on the velocity component and the stress component in the elastic wave equationx,vz,τxx,τzz,τxz)TWherein v isx,vzIs the velocity component of the elastic wave, τxx,τzz,τxzIs the stress component of the elastic wave;
a matrix wave equation obtaining subunit for obtaining an elastic wave equation
Figure FDA0002967599420000054
Conversion into a matrix wave equation
Figure FDA0002967599420000061
Wherein w is a second matrix,
Figure FDA0002967599420000062
p is the density of the mixture,
Figure FDA0002967599420000063
and
Figure FDA0002967599420000064
as a parameter of the modulus of elasticity, vpIs the velocity of longitudinal wave, vsFor shear wave velocity,. epsilon.and.delta.are Thomsen anisotropy parameters.
7. The apparatus of claim 5, wherein the discrete wave equation acquisition module comprises:
a discrete wave equation expansion subunit for expanding the matrix wave equation by Taylor series
Figure FDA0002967599420000065
Discrete wave equation expanded to time fourth order precision
Figure FDA0002967599420000066
In the formula (I), the compound is shown in the specification,
Figure FDA0002967599420000067
wherein the content of the first and second substances,
Figure FDA0002967599420000068
Figure FDA0002967599420000069
Figure FDA00029675994200000610
Figure FDA00029675994200000611
Figure FDA00029675994200000612
Figure FDA0002967599420000071
Figure FDA0002967599420000072
Figure FDA0002967599420000073
Figure FDA0002967599420000074
p is the density of the mixture,
Figure FDA0002967599420000075
Figure FDA0002967599420000076
and
Figure FDA0002967599420000077
as a parameter of the modulus of elasticity, vpIs the velocity of longitudinal wave, vsFor shear wave velocity,. epsilon.and. delta.are Thomsen anisotropy parameters, t is a time node, and. delta.t is a preset time interval.
8. The apparatus of claim 5, wherein the optimized finite difference coefficient solving module comprises:
the intermediate equation solving subunit is used for combining the plane wave solution and the spatial partial derivative to obtain an intermediate equation;
the optimization equation solving subunit is used for minimizing errors at two ends of the intermediate equation and combining a least square method to obtain an optimization equation;
and the optimization finite difference coefficient calculation subunit is used for solving the optimization equation to obtain the optimization finite difference coefficient.
CN201811542112.4A 2018-12-17 2018-12-17 Time-space high-order VTI rectangular grid finite difference method and device Active CN109725346B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811542112.4A CN109725346B (en) 2018-12-17 2018-12-17 Time-space high-order VTI rectangular grid finite difference method and device

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811542112.4A CN109725346B (en) 2018-12-17 2018-12-17 Time-space high-order VTI rectangular grid finite difference method and device

Publications (2)

Publication Number Publication Date
CN109725346A CN109725346A (en) 2019-05-07
CN109725346B true CN109725346B (en) 2021-06-18

Family

ID=66296222

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811542112.4A Active CN109725346B (en) 2018-12-17 2018-12-17 Time-space high-order VTI rectangular grid finite difference method and device

Country Status (1)

Country Link
CN (1) CN109725346B (en)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112526605B (en) * 2020-12-24 2022-09-02 广州海洋地质调查局 Method for simulating and exploring natural gas hydrate by adopting seismic numerical value
CN113536638B (en) * 2021-07-22 2023-09-22 北京大学 High-precision seismic wave field simulation method based on intermittent finite element and staggered grid

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103091711A (en) * 2013-01-24 2013-05-08 中国石油天然气集团公司 Method and device for full-wave-shape inversion
CN103630933A (en) * 2013-12-09 2014-03-12 中国石油天然气集团公司 Nonlinear optimization based time-space domain staggered grid finite difference method and device

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103091711A (en) * 2013-01-24 2013-05-08 中国石油天然气集团公司 Method and device for full-wave-shape inversion
CN103630933A (en) * 2013-12-09 2014-03-12 中国石油天然气集团公司 Nonlinear optimization based time-space domain staggered grid finite difference method and device

Non-Patent Citations (9)

* Cited by examiner, † Cited by third party
Title
2D各向同性和各向异性介质中高阶交错网格有限差分波场模拟;王馨;《中国优秀硕士学位论文全文数据库 基础科学辑》;20150215(第02期);第5-16页 *
3D acoustic wave modeling with a time-space-domain temporal high-order finite difference scheme;Shigang Xu et al.;《Journal of Geophysics and Engineering》;20180618;第1963-1976页 *
An improved rotated staggered-grid finite-difference method with fourth-order temporal accuracy for elastic-wave modeling in anisotropic media;Kai Gao et al.;《Journal of Computational Physics》;20170831;第363-368页 *
Pseudoacoustic tilted transversely isotropic modeling with optimal k-space operator-based implicit finite-difference schemes;Shigang Xu et al.;《GEOPHYSICS》;20180630;第83卷(第3期);第T139-T157页 *
一阶声波方程时间四阶精度差分格式的伪谱法求解;唐怀谷等;《石油地球物理勘探》;20170228;第52卷(第1期);第72-73页 *
一阶弹性波方程交错网格高阶差分解法;董良国等;《地球物理学报》;20000531;第43卷(第3期);第412-415页 *
交错网格任意阶导数有限差分格式及差分系数推导;杨庆节等;《吉林大学学报（地球科学版）》;20140131;第44卷(第1期);第376-382页 *
基于 k 空间算子补偿的 VTI 介质弹性波交错网格有限差分正演模拟;徐世刚等;《中国地球科学联合学术年会 2018》;20181031;第907-908页 *
横向各向同性介质有限差分法波场模拟方法研究;李宾;《中国优秀硕士学位论文全文数据库 基础科学辑》;20100315(第03期);第10-22、30-31页 *

Also Published As

Publication number Publication date
CN109725346A (en) 2019-05-07

Similar Documents

Publication Publication Date Title
Anderson et al. Time-reversal checkpointing methods for RTM and FWI
Min et al. Improved frequency-domain elastic wave modeling using weighted-averaging difference operators
Komatitsch et al. Spectral-element simulations of global seismic wave propagation—II. Three-dimensional models, oceans, rotation and self-gravitation
Casadei et al. A mortar spectral/finite element method for complex 2D and 3D elastodynamic problems
Alkhalifah Scanning anisotropy parameters in complex media
Liu et al. A mixed-grid finite element method with PML absorbing boundary conditions for seismic wave modelling
Nissen‐Meyer et al. A two‐dimensional spectral‐element method for computing spherical‐earth seismograms–I. Moment‐tensor source
Askan et al. Full waveform inversion for seismic velocity and anelastic losses in heterogeneous structures
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
CN109725346B (en) Time-space high-order VTI rectangular grid finite difference method and device
Fang et al. Elastic full-waveform inversion based on GPU accelerated temporal fourth-order finite-difference approximation
Long et al. A temporal fourth-order scheme for the first-order acoustic wave equations
Le Bouteiller et al. A discontinuous Galerkin fast-sweeping eikonal solver for fast and accurate traveltime computation in 3D tilted anisotropic media
Ma et al. Calculating Ray Paths for First‐Arrival Travel Times Using a Topography‐Dependent Eikonal Equation Solver
Zou et al. 3D elastic waveform modeling with an optimized equivalent staggered-grid finite-difference method
Di Bartolo et al. High-order finite-difference approximations to solve pseudoacoustic equations in 3D VTI media
Huang et al. A novel hybrid method based on discontinuous Galerkin method and staggered‐grid method for scalar wavefield modelling with rough topography
Abedi et al. A new acoustic assumption for orthorhombic media
Ma et al. A stable auxiliary differential equation perfectly matched layer condition combined with low-dispersive symplectic methods for solving second-order elastic wave equations
Xu et al. Time-space-domain temporal high-order staggered-grid finite-difference schemes by combining orthogonality and pyramid stencils for 3D elastic-wave propagation
Adcroft Numerical algorithms for use in a dynamical model of the ocean
He et al. A numerical dispersion-dissipation analysis of discontinuous Galerkin methods based on quadrilateral and triangular elements
Shukla et al. Modeling the wave propagation in viscoacoustic media: An efficient spectral approach in time and space domain
Xu et al. Pseudoacoustic tilted transversely isotropic modeling with optimal k-space operator-based implicit finite-difference schemes

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
GR01 Patent grant
GR01 Patent grant