CN107942375B - A kind of implicit time-space domain finite difference numerical simulation method of nonlinear optimization based on ACOUSTIC WAVE EQUATION - Google Patents
A kind of implicit time-space domain finite difference numerical simulation method of nonlinear optimization based on ACOUSTIC WAVE EQUATION Download PDFInfo
- Publication number
- CN107942375B CN107942375B CN201711142923.0A CN201711142923A CN107942375B CN 107942375 B CN107942375 B CN 107942375B CN 201711142923 A CN201711142923 A CN 201711142923A CN 107942375 B CN107942375 B CN 107942375B
- Authority
- CN
- China
- Prior art keywords
- time
- difference
- formula
- space
- implicit
- 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
Links
Classifications
-
- 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. analysis, for interpretation, for correction
-
- 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. analysis, for interpretation, for correction
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
-
- 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/51—Migration
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/67—Wave propagation modeling
- G01V2210/673—Finite-element; Finite-difference
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/70—Other details related to processing
- G01V2210/74—Visualisation of seismic data
Abstract
The invention discloses a kind of implicit time-space domain finite difference numerical simulation methods of nonlinear optimization based on ACOUSTIC WAVE EQUATION, comprising the following steps: (1) reads parameter；(2) the high-order discrete differential format of second time derivative is solved, and calculates higher difference coefficient；(3) second order spatial derivative is solved using the implicit discrete scheme in space, solves implicit difference coefficient；(4) the implicit difference recurrence algorithm of time high-order and space is constructed while being had, time-space domain dispersion relation is obtained；(5) time-space domain difference coefficient is solved using nonlinear optimization algorithm；(6) reflection is absorbed using mixed absorbing boundary, carries out recursion according to wave equation, obtain any time wave field and entire earthquake record；(7) wave field snapshot is recorded, earthquake record is exported and is terminated.Method proposed by the present invention can be substantially reduced under the premise of guaranteeing precision calculates the time, efficiently provides precision higher wave field for reverse-time migration and full waveform inversion method.
Description
Technical field
The invention belongs to geophysical exploration technical field of data processing, are related to a kind of based on the high-precision of solution ACOUSTIC WAVE EQUATION
Degree, high efficiency nonlinear optimize implicit time-space domain finite difference numerical simulation method.
Background technique
Currently based on Acoustic Wave-equation reverse-time migration method compared to based on one way wave equation offset imaging method and
Gram western Hough imaging method based on ray theory has shown that big advantage in terms of the complicated structures such as processing salt dome.
Based on the full waveform inversion method of ACOUSTIC WAVE EQUATION also increasingly by the attention of Geophysicist in terms of inverting reservoir parameter.
The main reason is that both methods takes full advantage of the kinematics and dynamic characteristic of seismic wave, consider in imaging and inverting
The various information such as primary wave, multiple wave, back wave and refracted wave, can disclose subsurface structure information and reservoir to the greatest extent
Parameter.Full waveform inversion and reverse-time migration method are directed to the solution of wave equation, and the calculating time is long, and calculation amount is huge, certain
Restrict the development of both technologies in degree.Therefore high-precision, efficient method for numerical simulation just seem and are even more important.Have
Calculus of finite differences is limited compared to other method for numerical simulation such as FInite Element, pseudo- spectrometries, with easy to operate, computational efficiency is high, is easy to simultaneously
Capable and required memory is small to wait many advantages, is widely used in the numerical solution of ACOUSTIC WAVE EQUATION.
Current finite difference calculus major defect is to be easy to frequency dispersion and stability is poor.Frequency dispersion can be divided into according to discrete objects
Time difference frequency dispersion and space parallax frequency dividing dissipate.Time frequency dispersion can be suppressed using lesser time step, but this equally brings
The huge calculating time.And it directlys adopt higher difference time proximity derivative to will lead to communication process unstable.Space frequency dispersion can
To be suppressed by increasing operator length or reducing spatial mesh size.But both strategies, which can all increase, calculates the time.Space frequency dispersion
It can also be suppressed by using implicit difference method.Compared to application of explicit difference method, implicit difference method can use lesser operator
Length come reach with the consistent precision of explicit difference, therefore facilitate reduce calculate the time.However, depositing due to time frequency dispersion
The raising of spatial accuracy can't effectively suppress whole frequency dispersion.
Summary of the invention
The purpose of the present invention is to reduce the problem that traditional finite differential method time frequency dispersion is serious, stability is poor
A kind of implicit time-space domain finite difference numerical simulation method of nonlinear optimization based on ACOUSTIC WAVE EQUATION of offer, it is proposed by the present invention
Method can be substantially reduced under the premise of guaranteeing precision calculates the time, efficiently mentions for reverse-time migration and full waveform inversion method
For the higher wave field of precision.
A kind of implicit time-space domain finite difference numerical simulation method of nonlinear optimization based on ACOUSTIC WAVE EQUATION of the invention, packet
Include following steps:
(1) parameter is read, the parameter includes rate pattern Parameter File, finite difference operator required for forward simulation
Length, nonlinear iteration solve the limits of error, wavelet function and the wavelet dominant frequency, forward modeling needed when difference coefficient and are used
Time and space step-length and earthquake record duration；
(2) the high-order discrete differential format of second time derivative is solved, and calculates higher difference coefficient；
(3) the solution second order spatial derivative implicit difference scheme that Claerbout (1985) propose is used for reference, conventional explicit high
On the basis of order difference, Second-Order Central Difference format is introduced in denominator, construction space implicit difference scheme solves second order spatial and leads
Number, for reducing finite difference operator length；To further decrease frequency dispersion, it is based on space dispersion relation, solves implicit difference system
Number；
(4) it is based on the implicit discrete scheme of the time high-order discrete scheme and space, brings difference scheme into ACOUSTIC WAVE EQUATION,
There is the implicit difference recurrence algorithm of time high-order and space simultaneously；Using plane wave analysis theories, by plane wave equation
Bring dispersion relation into, abbreviation obtains time-space domain dispersion relation；
(5) using the difference coefficient in time and space difference discrete format as variable, by minimization equation dispersion relation,
Establish the non-linear objective function about time and space difference coefficient；The difference system being calculated with step (2) and step (3)
Number is initial value, is iterated solution to time and space difference coefficient using Solvopt nonlinear optimization algorithm；
(6) capable absorption is reflected into boundary using mixed absorbing boundary, is carried out according to the discrete scheme of wave equation
Recursion, obtain any time wave field and entire earthquake record；
(7) wave field snapshot is recorded, earthquake record is exported and is terminated.
In step (2), centered difference is used when solving second time derivative, form is as follows:
Wherein,For pressure field, τ is time step, and t is the time；In order to further increase time discrete
Diamond shape difference operator is introduced into formula (1) by precision, obtains the following discrete formula of time high-order:
Wherein, h is spatial sampling step-length, and v is speed, am,nFor finite difference coefficient；N is that the difference in diamond shape operator is calculated
Sub- length, value is bigger, and time difference approximation quality is higher, but increased calculation amount also increases；N and m is in diamond shape operator
Coordinate position sequence.
Above-mentioned finite difference coefficient am,nBy the way that formula (2) is calculated using Taylor expansion, formula is as follows:
Wherein, j=2,3 ..., N, ξ=1,2 ..., int (j/2), r are the bright number in library, and r=v τ/h, int are rounding operations；
It is discrete to time-derivative progress by using formula (2), and using difference coefficient shown in formula (3), time discrete precision reaches
To high-order, it to be used for pressing time frequency dispersion.
In step (3), approximation is carried out to second order spatial derivative using the implicit discrete scheme of such as down space, form is as follows:
Wherein, a '0、a′mIt is implicit difference coefficient with b, h is spatial mesh size, and x is direction coordinate, and M is difference operator length；
The corresponding space dispersion relation of formula (4) is as follows:
Wherein, β=kxH, kxFor horizontal direction wave number；The corresponding space dispersion relation of minimization formula (4), objective function
For the linear convex function about difference coefficient, optimal difference coefficient is obtained by solving following formula:
Wherein, n=1,2 ..., M+1, d=[d1,d2,…,dM,dM+1]T=[a '1,a′2,…,a′M-1,a′M,b]T, q is to be situated between
Natural number between 0 and 1, controls limit of integration,
In step (4), when solving time-derivative, calculated using formula (2), when solution room derivative, using formula (4)
It calculates, formula (2) and formula (4) is brought into ACOUSTIC WAVE EQUATION and carry out plane wave analysis, process is as follows:
In homogeneous model, ACOUSTIC WAVE EQUATION pressure P meets following plane wave equation,
Wherein, ω is angular frequency,It is imaginary unit, kzFor vertical direction wave number；M, n and l is position and time
Sequence；
It brings formula (8) into formula (2) and carries out abbreviation and obtain following formula:
Equally, space derivation formula is analyzed, obtains following dispersion relation:
Above-mentioned (9), (10) and (11) are brought into ACOUSTIC WAVE EQUATION and simplification and recombination obtain following time-space domain dispersion relation,
Wherein, ω is angular frequency, kxAnd kzRespectively along the wave number of x and z directions.
Step (5) specific steps are as follows:
Firstly, difference coefficient a 'mIt is calculated with b by formula (6), difference coefficient am,0And am,nIt is calculated by formula (3)
It arrives；It is asked using the above-mentioned difference coefficient acquired as initial input using nonlinear optimization algorithm based on time-space domain dispersion relation (12)
Solve time-space domain difference coefficient:
Objective function is defined as follows,
Wherein, θ is propagation angle, and q is the real number between 0 and 1,
Solution in formula (13) is calculated by nonlinear optimization algorithm, and formula is as follows:
(am,0,am,n,a′m,b,r)optimal=Minimize [J2(am,0,am,n,a′m,b,r)] (16)
Above-mentioned nonlinear optimization algorithm is specifically sought using SolvOpt nonlinear optimization algorithm, using following strategy
To reduce operation time required for calculating difference coefficient:
It is assumed that the rate pattern of simulated object (rate pattern is the object of forward simulation, is input in program as input)
Speed value range in parameter is [v1,v2], then when calculating difference coefficient, difference coefficient is calculated at interval of dv, form is such as
Under:
Positioned at [v1,v1+ dv] between the corresponding difference coefficient of speed and v1The difference coefficient at place is consistent.Even now operation one
Determine to have lost computational accuracy in degree, but calculate the time since this method greatly reduces, computational efficiency can be significantly improved.And
The calculating of complex model finds that calculating error caused by this strategy is simultaneously little.
In step (6), entire zoning is divided into interior zone, borderline region and transitional region, specific method
It is as follows:
(6-1) calculates difference coefficient using discrete formula (2) and (4) discrete time and space derivation, using formula (17),
Wave field recursion is carried out, the round trip wave field of any time at all areas is obtained；
(6-2) solves sound wave one way wave equation, calculates the one way wave equation at borderline region and transitional region, calculation formula
And method is as follows:
By taking left margin and the upper left corner as an example, left margin one way wave equation are as follows:
Wherein, x and z represent coordinate；Upper left corner one way wave equation are as follows:
Other boundaries and corner point one-way wave are similar to be obtained, and carries out corresponding solution；Using above-mentioned absorbing boundary condition
And recursion is carried out according to the method described above, obtain any time wave field value；
(6-3) calculates final wave field；
Boundary wave field is one way wave field value, and wave field is two-way wave field value, wave field at middle transition region at interior zone
It is obtained by following linear combination:
P=BiPtwo+(1-Bi)Pone, (20)
Wherein, PtwoFor two-way wave field, PoneFor one way wave field, BiFor weight coefficient, Bi=(i-1)/Nb, i=2 ..., Nb,
NbFor the transitional region absorbing boundary number of plies；Transitional region from inside to outside, BiBy (Nb-1)/NbChange to 1/NbIf Nb=1, then
It is one-way wave absorbing boundary condition that mixed absorbing boundary, which is degenerated, at this time.
Its remarkable advantage is the present invention compared with prior art:
A) time precision is high, and frequency dispersion is small.
B) stability is good, and bigger time step can be used in simulation.
C) method space frequency dispersion is compared traditional Taylor expansion and is preferably suppressed, and can use smaller operator length,
It effectively saves and calculates the time, computational efficiency is high.
Detailed description of the invention
Fig. 1 is the implicit time-space domain finite difference numerical simulation side of a kind of nonlinear optimization based on ACOUSTIC WAVE EQUATION of the invention
Method workflow schematic diagram；
Fig. 2 (a) is conventional implicit difference method dispersion curve step change figure at any time；
Fig. 2 (b) is the method for the present invention dispersion curve step change figure at any time；
Fig. 3 is the method for the present invention stability condition curve and conventional method comparison diagram；
Shock curve is the same as reference at the homogeneous model (1200m, 2250m) that Fig. 4 (a) is calculated for conventional implicit difference method
Solve comparison diagram；
Shock curve is the same as with reference to solution comparison at the homogeneous model (1200m, 2250m) that Fig. 4 (b) is calculated for this method
Figure；
Fig. 4 (c) is that conventional implicit difference method and the method for the present invention simulation shock curve are square with the statistics referred between solution
Root error comparison diagram；
Fig. 5 is SEG/EAGE salt dome rate pattern；
Fig. 6 (a) is SEG/EAGE In A Salt-dome Model 2.4s moment reference wave field snapshot；
Fig. 6 (b) is 2.4s moment wave field snapshot after the conventional implicit difference method of SEG/EAGE In A Salt-dome Model application；
Fig. 6 (c) is 2.4s moment wave field snapshot after SEG/EAGE In A Salt-dome Model application the method for the present invention；
Fig. 7 (a) is SEG/EAGE In A Salt-dome Model using earthquake record at (500m, 4500m) after conventional implicit difference method
Seismogram and reference solution comparison diagram；
Fig. 7 (b) is earthquake record seismogram at (500m, 4500m) after SEG/EAGE In A Salt-dome Model application the method for the present invention
With reference solution comparison diagram.
Specific embodiment
Specific embodiments of the present invention will be described in further detail with reference to the accompanying drawings and examples.
The present invention provides a kind of implicit time-space domain finite difference numerical simulation sides of the nonlinear optimization based on ACOUSTIC WAVE EQUATION
Method, this method on the basis of Conventional Time Second-Order Discrete format, joined diamond shape difference operator first, and having derived has high-order
The time discrete format and difference coefficient of precision.When solution room derivative, the present invention is being protected using the implicit difference method of optimization
Reduce operator length under the premise of card precision.In conjunction with the implicit space difference method of high-order, the present invention has obtained while having had high-order
The recurrence algorithm and its time-space domain dispersion relation of time precision and implicit spatial accuracy.From time-space domain dispersion relation, this hair
It is bright to use nonlinear optimization algorithm, propose the method and strategy for calculating time-space domain difference coefficient.Compared to conventional method, the present invention
Method precision is high, and stability is good.To avoid calculating the difference coefficient at each speed in complex model, the present invention is further mentioned
A kind of selection speed strategy saved and calculate the time is supplied.
As shown in Figure 1, the implicit time-space domain finite difference of a kind of nonlinear optimization based on ACOUSTIC WAVE EQUATION proposed by the present invention
Method for numerical simulation comprises the following specific steps that:
(1) read parameter, including forward simulation required for rate pattern Parameter File, finite difference operator length (when
Between derivative and space derivation), nonlinear iteration solve difference coefficient when need termination condition, wavelet function and wavelet dominant frequency,
The parameters such as time and spatial mesh size used by forward modeling and earthquake record duration；
(2) it is solved for the high-order discrete differential format construction and difference coefficient of time-derivative:
Conventional method uses centered difference when solving 2 rank time-derivative, and form is as follows:
Wherein,For pressure field, τ is time step, and t is the time.In order to further increase time discrete
Diamond shape difference operator is introduced into formula (1) by precision, the present invention, proposes the following discrete formula of time high-order:
Wherein, h is spatial sampling step-length, and v is speed, and N is the difference operator length in diamond shape operator, and value is bigger, when
Between difference approximation precision it is higher, but increased calculation amount also increases.N and m is the coordinate position sequence in diamond shape operator.Finite difference
Divide coefficient am,nIt can be by being obtained to formula (2) using Taylor expansion, formula is as follows:
Wherein, r is the bright number in library, and r=v τ/h, int are rounding operations.By using formula (2) to time-derivative carry out from
It dissipates, and using difference coefficient shown in formula (3), time discrete precision can achieve high-order, effective pressing time frequency dispersion.
(3) implicit difference method solution room derivative is utilized, to further decrease frequency dispersion, space dispersion relation is based on, leads to
Cross Optimization Method implicit difference coefficient:
For second order spatial derivative, to effectively reduce operator length, the present invention uses following implicit difference discrete scheme pair
Space second dervative carries out approximation, and form is as follows:
Wherein, a '0、a′mIt is implicit difference coefficient with b.H is spatial mesh size, and x is direction coordinate, and M is difference operator length.
Conventional method is seeking difference coefficient a '0、a′mWith when b use Taylor expansion, precision is low, is easy to frequency dispersion.Herein using optimization
Method solves, and the corresponding space dispersion relation of formula (4) is as follows:
Wherein, β=kxH, kxFor horizontal direction wave number.The corresponding space dispersion relation of minimization formula (4), objective function
For the linear convex function about difference coefficient, can be obtained by solving following formula:
Wherein, d=[d1,d2,…,dM,dM+1]T=[a '1,a′2,…,a′M-1,a′M,b]T, q is oneself between 0 and 1
So number, controls limit of integration,
(4) the implicit difference recurrence algorithm of time high-order and space is constructed while being had, and obtains time-space domain dispersion relation:
It when solving time-derivative, is calculated using formula (2), when solution room derivative, is calculated using formula (4), it will be public
Formula (2) and formula (4) are brought into ACOUSTIC WAVE EQUATION and carry out plane wave analysis method, available following time-space domain dispersion relation,
Wherein, ω is angular frequency, kxAnd kzRespectively along the wave number of x and z directions.
(5) from time-space domain dispersion relation, the difference coefficient in the above way acquired is initial guess, and use is non-linear
Optimization algorithm is iterated Optimization Solution to time-space domain difference coefficient:
Firstly, difference coefficient a 'mIt is calculated with b by formula (6), difference coefficient am,0And am,nIt is calculated by formula (3)
It arrives.Difference coefficient a is being obtained by the above methodm,0、am,n、a′mAfter b, frequency dispersion error still can pass through nonlinear optimization
Algorithm is further decreased, and objective function is defined as follows,
Wherein, θ is propagation angle, and q is the real number between 0 and 1,
Formula (9) can be calculated by nonlinear optimization algorithm, and formula is as follows:
(am,0,am,n,a′m,b,r)optimal=Minimize [J2(am,0,am,n,a′m,b,r)].(12)
Nonlinear optimization algorithm is specifically sought using SolvOpt nonlinear optimization algorithm, and SolvOpt is a kind of non-thread
Property Local Optimization Algorithm, objective function fast convergence rate facilitate rapid solving nonlinear optimization difference coefficient.
For complex model, above-mentioned operation needs calculate each velocity amplitude, due to being related to iteratively solving, when operation
Between it is very long.Present invention proposition reduces operation time required for calculating difference coefficient using following strategy.
It is assumed that speed value range involved in the rate pattern of simulated object is [v1,v2], then calculating difference coefficient
When, difference coefficient is calculated at interval of dv, form is as follows:
Positioned at [v1,v1+ dv] between the corresponding difference coefficient of speed and v1The difference coefficient at place is consistent, other speed areas
Between it is similar with the interval computation method.By above-mentioned strategy, the time required for calculating difference coefficient is substantially reduced, and is used for simultaneously
The memory of storage difference coefficient also accordingly reduces.Meanwhile complex model shows that the strategy can't significantly reduce difference method
Precision.
(6) absorbing boundary condition appropriate is used, wave field extrapolation is carried out:
After acquiring difference coefficient, recursion is carried out according to wave equation recurrence algorithm, obtains the wave field of any time and entire
Earthquake record.The present invention absorbs reflection using the mixed absorbing boundary that Liu and Sen (2010) is proposed.
(7) wave field snapshot, output earthquake record and terminator are recorded.
The following are one embodiment of the present of invention, illustrate based on a kind of implicit space-time of the nonlinear optimization based on ACOUSTIC WAVE EQUATION
The realization process and final result of domain finite difference numerical simulation method.For the ease of illustrating, hereinafter, conventional method is remembered
It is ISD, method proposed by the present invention is denoted as ITSD-2.
Fig. 2 (a) and Fig. 2 (b) is method proposed by the present invention and the comparative analysis of conventional method frequency dispersion error, ordinate in figure
For numerical simulation speed and true velocity ratio, the precision of method is reflected.When ordinate is 1.0, method is without frequency dispersion at this time.
The position that ordinate deviates 1.0 is remoter, then frequency dispersion is bigger.In figure, speed parameter 1500m/s, spatial sampling step-length is 15m,
Operator length M=6.Dispersion curve when using different time step-length when the direction of propagation is 45 degree is shown in figure.Comparison discovery,
Compared to conventional method, the method proposed in the present invention significantly improves the precision of conventional method, has suppressed numerical solidification.
Fig. 3 is method proposed by the present invention and conventional method Stability Contrast figure, in figure ordinate be stability because
Son, when the bright number in library is greater than the stability factor, difference method is unstable.The value is bigger, and method is more stable.Comparison discovery, passes
System method stability is worst, and stability factor is minimum.Method stability proposed by the present invention is more preferable, and it is bigger to facilitate method use
Time step.
Fig. 4 a and Fig. 4 b are simulation wave field seismogram and reference solution comparison diagram at homogeneous model (1200m, 2250m).Model
Size is 4500m × 4500m, spatial mesh size 15m, time step 3ms, the wavelet used be dominant frequency be 18Hz rake it is sub
Wave, the heart excites in a model.With reference to solution to be calculated by Cagniard-de Hoop (De Hoop, 1960) method.The present invention
Middle method uses M=6 and N=2.Comparison discovery, conventional method frequency dispersion error is obvious, and it is poor to coincide with reference solution.The present invention proposes
Method then significantly improve precision, coincide with reference solution more preferable.Fig. 4 c is the statistics root mean square between each method and reference solution
Error, hence it is evident that method error proposed by the present invention is smaller.
Fig. 5 is commonly used SEG/EAGE In A Salt-dome Model in seismic prospecting numerical simulation.Model size be 5000m ×
15000m, spatial mesh size 25m, the wavelet used are the Ricker wavelet of 15Hz for dominant frequency, are excited at earth's surface center, time step
For 2.5ms.Method uses M=8 and N=2 in the present invention.Fig. 6 (a) to Fig. 6 (c) is each method propagating wavefield at the 2.4s moment
Snapshot is calculated for conventional explicit time-space domain method using the time step and very big operator length of very little with reference to solution.Together
With reference to solution comparison discovery, method obviously has higher precision in the present invention, and white arrow meaning region frequency dispersion compares tradition side
Method is obviously small.Fig. 7 (a) and Fig. 7 (b) is that each method earthquake record seismogram at (500m, 4500m) is compared with reference to solution comparison
Scheme, the error amount shown in figure is the root-mean-square error of statistics, hence it is evident that method root-mean-square error proposed by the present invention is smaller and non-
Linear optimization method is higher compared to linear optimization method precision, and error is smaller.
To reduce frequency dispersion error, traditional implicit difference method use the root-mean-square error that is calculated when 1ms time step for
1.1824e-4, the CPU time of consumption is 403.58s.Method proposed by the present invention when using 2.5ms time step, calculating
Root-mean-square error is 1.0206e-4, the CPU time of consumption is 240.27s.It can be found that method proposed by the present invention is maintaining
It under the premise of precision, significantly reduces and calculates the time, efficiency is greatly enhanced.
This method can significantly improve computational efficiency under the premise of guaranteeing precision, be the huge reverse-time migration of calculation amount
Method and full waveform inversion method provide convenience.The explanation being not directed in a specific embodiment of the invention belongs to known in this field
Technology, can refer to well-known technique and be implemented.
The above specific embodiment and embodiment are to a kind of nonlinear optimization based on ACOUSTIC WAVE EQUATION proposed by the present invention
The specific support of implicit time-space domain finite difference numerical simulation method and technology thought, cannot limit protection model of the invention with this
It encloses, all any equivalent variationss or equivalent according to the technical idea provided by the invention, done on the basis of the technical program
Change still falls within the range of technical solution of the present invention protection.
Claims (7)
1. a kind of implicit time-space domain finite difference numerical simulation method of nonlinear optimization based on ACOUSTIC WAVE EQUATION, which is characterized in that
Including the following steps:
(1) parameter is read, the parameter includes that rate pattern Parameter File, finite difference operator required for forward simulation are long
Degree, nonlinear iteration solve used by the limits of error, wavelet function and wavelet dominant frequency, the forward modeling needed when difference coefficient
Time and space step-length and earthquake record duration；
(2) the high-order discrete differential format of second time derivative is solved, and calculates higher difference coefficient；
(3) it solves second order spatial derivative implicit difference scheme and introduces second-order central in denominator on the basis of explicit higher difference
Difference scheme, construction space implicit difference scheme solves second order spatial derivative, for reducing finite difference operator length；For into one
Step reduces frequency dispersion, is based on space dispersion relation, solves implicit difference coefficient；
(4) it is based on the implicit discrete scheme of the time high-order discrete scheme and space, difference scheme is brought into ACOUSTIC WAVE EQUATION, obtains
There is the implicit difference recurrence algorithm of time high-order and space simultaneously；Using plane wave analysis theories, plane wave equation is brought into
Dispersion relation, abbreviation obtain time-space domain dispersion relation；
(5) it using the difference coefficient in time and space difference discrete format as variable, by minimization equation dispersion relation, establishes
Non-linear objective function about time and space difference coefficient；It is with the difference coefficient that step (2) and step (3) are calculated
Initial value is iterated solution to time and space difference coefficient using Solvopt nonlinear optimization algorithm；
(6) capable absorption is reflected into boundary using mixed absorbing boundary, carries out recursion according to the discrete scheme of wave equation,
Obtain any time wave field and entire earthquake record；
(7) wave field snapshot is recorded, earthquake record is exported and is terminated；
In step (2), centered difference is used when solving second time derivative, form is as follows:
Wherein,For pressure field, t is time step, and t is the time；In order to further increase time discrete precision,
Diamond shape difference operator is introduced into formula (1), obtains the following discrete formula of time high-order:
Wherein, h is spatial sampling step-length, and v is speed, am,nFor finite difference coefficient；N is that the difference operator in diamond shape operator is long
Degree；N and m is the coordinate position sequence in diamond shape operator.
2. the implicit time-space domain finite difference numerical simulation side of the nonlinear optimization according to claim 1 based on ACOUSTIC WAVE EQUATION
Method, which is characterized in that the finite difference coefficient am,nBy the way that formula (2) is calculated using Taylor expansion, formula is such as
Under:
Wherein, j=2,3, L, N, ξ=1,2, L, int (j/2), r are the bright number in library, and r=v τ/h, int are rounding operations；By adopting
Time-derivative is carried out with formula (2) discrete, and using difference coefficient shown in formula (3), time discrete precision reaches high-order,
For pressing time frequency dispersion.
3. the implicit time-space domain finite difference numerical simulation side of the nonlinear optimization according to claim 2 based on ACOUSTIC WAVE EQUATION
Method, which is characterized in that in step (3), approximation is carried out to second order spatial derivative using the implicit discrete scheme of such as down space, form is such as
Under:
Wherein, a '0、a′mIt is implicit difference coefficient with b, h is spatial mesh size, and x is direction coordinate, and M is difference operator length；Formula
(4) corresponding space dispersion relation is as follows:
Wherein, β=kxH, kxFor horizontal direction wave number；The corresponding space dispersion relation of minimization formula (4), objective function are to close
In the linear convex function of difference coefficient, optimal difference coefficient is obtained by solving following formula:
Wherein, n=1,2, L, M+1, d=[d1,d2,L,dM,dM+1]T=[a '1,a′2,L,a′M-1,a′M,b]T, q is between 0 and 1
Between natural number, control limit of integration,
4. the implicit time-space domain finite difference numerical simulation side of the nonlinear optimization according to claim 3 based on ACOUSTIC WAVE EQUATION
Method, which is characterized in that in step (4), when solving time-derivative, calculated using formula (2), when solution room derivative, used
Formula (4) calculates, and formula (2) and formula (4) are brought into ACOUSTIC WAVE EQUATION and carry out plane wave analysis, process is as follows:
In homogeneous model, ACOUSTIC WAVE EQUATION pressure P meets following plane wave equation,
Wherein, ω is angular frequency,It is imaginary unit, kzFor vertical direction wave number；M, n and l is position and time series；
It brings formula (8) into formula (2) and carries out abbreviation and obtain following formula:
Equally, space derivation formula is analyzed, obtains following dispersion relation:
Above-mentioned (9), (10) and (11) are brought into ACOUSTIC WAVE EQUATION and simplification and recombination obtain following time-space domain dispersion relation,
Wherein, ω is angular frequency, kxAnd kzRespectively along the wave number of x and z directions.
5. the implicit time-space domain finite difference numerical simulation side of the nonlinear optimization according to claim 4 based on ACOUSTIC WAVE EQUATION
Method, which is characterized in that step (5) specific steps are as follows:
Firstly, difference coefficient a 'mIt is calculated with b by formula (6), difference coefficient am,0And am,nIt is calculated by formula (3)；Base
In time-space domain dispersion relation (12), using the above-mentioned difference coefficient acquired as initial input,
Time-space domain difference coefficient is solved using nonlinear optimization algorithm:
Objective function is defined as follows,
Wherein, θ is propagation angle, and q is the real number between 0 and 1,
Solution in formula (13) is calculated by nonlinear optimization algorithm, and formula is as follows:
(am,0,am,n,a′m,b,r)optimal=Minimize [J2(am,0,am,n,a′m,b,r)] (16)。
6. the implicit time-space domain finite difference numerical simulation side of the nonlinear optimization according to claim 5 based on ACOUSTIC WAVE EQUATION
Method, which is characterized in that the nonlinear optimization algorithm is specifically sought using SolvOpt nonlinear optimization algorithm, using such as
Lower strategy come reduce calculate difference coefficient required for operation time:
It is assumed that the speed value range in the rate pattern parameter of simulated object is [v1,v2], then when calculating difference coefficient, often
It is spaced dv and calculates difference coefficient, form is as follows:
Positioned at [v1,v1+ dv] between the corresponding difference coefficient of speed and v1The difference coefficient at place is consistent.
7. the implicit time-space domain finite difference numerical simulation side of the nonlinear optimization according to claim 6 based on ACOUSTIC WAVE EQUATION
Method, which is characterized in that in step (6), entire zoning is divided into interior zone, borderline region and transitional region, specifically
Method it is as follows:
(6-1) calculates difference coefficient using discrete formula (2) and (4) discrete time and space derivation, using formula (17), carries out
Wave field recursion obtains the round trip wave field of any time at all areas；
(6-2) solves sound wave one way wave equation, calculates the one way wave equation at borderline region and transitional region, calculation formula and side
Method is as follows:
By taking left margin and the upper left corner as an example, left margin one way wave equation are as follows:
Wherein, x and z represent coordinate；Upper left corner one way wave equation are as follows:
Other boundaries and corner point one-way wave are similar to be obtained, and carries out corresponding solution；Using above-mentioned absorbing boundary condition and press
Recursion is carried out according to the above method, obtains any time wave field value；
(6-3) calculates final wave field；
Boundary wave field is one way wave field value, and wave field is two-way wave field value at interior zone, and wave field passes through at middle transition region
Following linear combination obtains:
P=BiPtwo+(1-Bi)Pone, (20)
Wherein, PtwoFor two-way wave field, PoneFor one way wave field, BiFor weight coefficient, Bi=(i-1)/Nb, i=2, L, Nb, NbFor mistake
Cross the region absorbing boundary number of plies；Transitional region from inside to outside, BiBy (Nb-1)/NbChange to 1/NbIf Nb=1, then it mixes at this time
Closing absorbing boundary condition and degenerating is one-way wave absorbing boundary condition.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711142923.0A CN107942375B (en) | 2017-11-17 | 2017-11-17 | A kind of implicit time-space domain finite difference numerical simulation method of nonlinear optimization based on ACOUSTIC WAVE EQUATION |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711142923.0A CN107942375B (en) | 2017-11-17 | 2017-11-17 | A kind of implicit time-space domain finite difference numerical simulation method of nonlinear optimization based on ACOUSTIC WAVE EQUATION |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107942375A CN107942375A (en) | 2018-04-20 |
| CN107942375B true CN107942375B (en) | 2019-04-30 |
Family
ID=61931635
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201711142923.0A Active CN107942375B (en) | 2017-11-17 | 2017-11-17 | A kind of implicit time-space domain finite difference numerical simulation method of nonlinear optimization based on ACOUSTIC WAVE EQUATION |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107942375B (en) |
Families Citing this family (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107976710B (en) * | 2017-11-17 | 2019-05-28 | 河海大学 | A kind of implicit time-space domain finite difference numerical simulation method of linear optimization based on ACOUSTIC WAVE EQUATION |
| CN110579796A (en) * | 2018-06-08 | 2019-12-17 | 中国科学院地质与地球物理研究所 | wave field simulation method, device and equipment for expanding finite difference stability condition |
| CN110609325B (en) * | 2018-06-15 | 2021-02-09 | 中国石油化工股份有限公司 | Elastic wave field numerical simulation method and system |
| CN110858002B (en) * | 2018-08-23 | 2022-07-15 | 中国科学院地质与地球物理研究所 | Wave field simulation method, device and equipment for expanding explicit differential stability condition |
| CN109725351A (en) * | 2018-12-18 | 2019-05-07 | 中国石油天然气集团有限公司 | A kind of the determination method, apparatus and system of 3D elastic wave mixed absorbing boundary |
| CN110285876A (en) * | 2019-07-01 | 2019-09-27 | 中国人民解放军军事科学院国防科技创新研究院 | A kind of acquisition methods of ocean acoustic field all-wave solution |
| CN112859170B (en) * | 2021-03-24 | 2022-08-02 | 西南石油大学 | High-precision wave field numerical simulation method based on time domain high-order finite difference method |
| CN113987841B (en) * | 2021-12-24 | 2022-04-19 | 苏州浪潮智能科技有限公司 | Method for solving phonon heat transport at interface and storage medium |
| CN114444351B (en) * | 2022-01-11 | 2023-03-31 | 中国空气动力研究与发展中心计算空气动力研究所 | Shock wave noise simulation method based on CCSSR-HW-6-BOO format |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102062875A (en) * | 2010-11-30 | 2011-05-18 | 中国石油集团川庆钻探工程有限公司 | Forward computation method of motion equation of elastic wave on relief surface |
| CN102508296A (en) * | 2011-11-14 | 2012-06-20 | 中国石油天然气股份有限公司 | Method and device for analyzing dispersion and attenuation of unsaturated double-porosity medium earthquake waves |
| CN105911584A (en) * | 2015-09-25 | 2016-08-31 | 中国科学院地质与地球物理研究所 | Implicit staggered-grid finite difference elastic wave numerical simulation method and device |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9864082B2 (en) * | 2008-11-06 | 2018-01-09 | Pgs Geophysical As | Fourier finite-difference migration for three dimensional tilted transverse isotropic media |
| CN102590862B (en) * | 2012-01-19 | 2014-03-19 | 中国科学院地质与地球物理研究所 | Prestack time migration method for compensating absorptive attenuation |
| CN103412328B (en) * | 2013-08-01 | 2016-04-20 | 中国石油天然气集团公司 | Wavenumber domain based on staggering mesh finite-difference algorithm protects amplitude wave field separation method |
-
2017
- 2017-11-17 CN CN201711142923.0A patent/CN107942375B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102062875A (en) * | 2010-11-30 | 2011-05-18 | 中国石油集团川庆钻探工程有限公司 | Forward computation method of motion equation of elastic wave on relief surface |
| CN102508296A (en) * | 2011-11-14 | 2012-06-20 | 中国石油天然气股份有限公司 | Method and device for analyzing dispersion and attenuation of unsaturated double-porosity medium earthquake waves |
| CN105911584A (en) * | 2015-09-25 | 2016-08-31 | 中国科学院地质与地球物理研究所 | Implicit staggered-grid finite difference elastic wave numerical simulation method and device |
Non-Patent Citations (1)
| Title |
|---|
| 声波方程隐式时空域有限差分正演模拟方法;王恩江，等;《2017中国地球科学联合学术年会论文集》;20171031;732-733 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107942375A (en) | 2018-04-20 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107942375B (en) | A kind of implicit time-space domain finite difference numerical simulation method of nonlinear optimization based on ACOUSTIC WAVE EQUATION | |
| CN107976710B (en) | A kind of implicit time-space domain finite difference numerical simulation method of linear optimization based on ACOUSTIC WAVE EQUATION | |
| CN107479092B (en) | A kind of frequency domain high order ACOUSTIC WAVE EQUATION the Forward Modeling based on directional derivative | |
| AU2012220584B2 (en) | Sensitivity kernel-based migration velocity analysis in 3D anisotropic media | |
| CN108535775B (en) | Non-stationary seismic data sound impedance inversion method and device | |
| CN102937720B (en) | Well control improves the method for seismic data resolution | |
| CN107817526B (en) | Prestack seismic gather segmented amplitude energy compensation method and system | |
| CN105425289B (en) | The method and apparatus for determining low frequency wave impedance | |
| CN104268412B (en) | A kind of angle gathers ray tomography migration velocity analysis method and device | |
| CN108693555B (en) | Intelligent time-varying blind deconvolution wideband processing method and processing device | |
| CN108415073B (en) | Angle domain back scattering offset imaging method and device | |
| CN107807390A (en) | The processing method and system of geological data | |
| CN108196304B (en) | A kind of multiple wave drawing method and device | |
| CN109490956A (en) | A kind of Acoustic Wave-equation the Forward Modeling and device based on staggered-mesh | |
| CN102866426B (en) | A kind of method utilizing AVO wide-angle road set analysis rock mass hydrocarbon information | |
| MXPA04009215A (en) | Method for seismic migration using explicit depth extrapolation operators with dynamically variable operator length. | |
| CHANG et al. | 3D ACOUSTIC REVERSE‐TIME MIGRATION 1 | |
| CN108957539B (en) | Ray tracing method and device in chromatography migration velocity analysis | |
| CN106814390A (en) | Staggered-mesh the Forward Modeling based on time-space domain optimization | |
| CN104570090B (en) | The extraction of full waveform inversion noise filter operator and the method filtered using its noise | |
| WO2015155597A2 (en) | Attenuating pseudo s-waves in acoustic anisotropic wave propagation | |
| Shragge et al. | Wave-equation migration from topography | |
| CN107272056A (en) | A kind of method that initial model is set up based on Duo Jing stratum transformation factor | |
| CN110109179A (en) | Bandwidth compensation processing method, device and equipment | |
| CN106547021A (en) | Based on the method and apparatus that individual well convolution algorithm sets up initial model |
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 |