CN103853930A - Numerical simulation method and device for earthquake vector wave field - Google Patents

Abstract

The invention provides a numerical simulation method and device for an earthquake vector wave field. The numerical simulation method comprises the following steps: performing L-times auto-convolution operation on an original window function, thereby acquiring an auto-convolution window function; performing weighted array operation on the auto-convolution window function and the original window function, thereby acquiring an auto-convolution combined window function; cutting off a space convolution sequence of a pseudo-spectral method by the auto-convolution combined window function, thereby acquiring a finite difference operator; performing numerical simulation on the earthquake vector wave field according to the finite difference operator. According to the numerical simulation method provided by the invention, the technical problems of the prior art that the value dispersion of a numerical simulation result of the finite difference earthquake vector wave field is larger and the numerical simulation precision is low because the cutting of the imminent finite difference operator through the window function cannot consider the spectrum coverage area and precision error stability are solved, and the technical effect of effectively increasing the simulation precision of the numerical simulation for the earthquake vector wave field is achieved.

Description

Earthquake vector wave Numerical Simulation method and apparatus

Technical field

The present invention relates to numerical simulation technology field, particularly a kind of earthquake vector wave Numerical Simulation method and apparatus.

Background technology

Along with deepening continuously that oil gas is surveyed, also more and more urgent for the demand of more high-precision imaging and inverting, and also become particularly important as the seismic wave field numerical simulation technology on imaging and inversion method basis.Seismic wave field numerical simulation is always a popular domain of study of seismology, understand the propagation law of seismic event for people, explain actual seismic data, characterize underground medium domain lithology, and earth resources exploitation etc. all has important theory and realistic meaning.

Conventional method for numerical simulation mainly contains three classes: how much rays methods, integral equation method, Wave Equation Method and method of finite difference, wherein, the analog result of Wave Equation Method is because therefore the kinematics and dynamics feature that has comprised wave field is applied comparatively extensive.

The numerical simulation technology of method of finite difference starts from twentieth century end of the sixties, and after having experienced the development of nearly half a century, the method has been widely applied in earthquake numerical simulation and migration and imaging techniques.And Explicit finite difference method is simple because of its algorithm, it is convenient that programming realizes, and it is lower to assess the cost, and becomes comparatively popular method for numerical simulation.

Up to this point, the inevitable problem that method of finite difference faces is exactly numerical value frequency dispersion, the emphasis of the research of people to finite difference is for many years exactly how to suppress numerical value frequency dispersion, the generation of numerical value frequency dispersion is owing to approaching differentiating operator by difference operator, after blocking, will inevitably cause the generation of error, this is to be determined by the essence of method of finite difference.Generally, in order to reduce numerical value frequency dispersion, can adopt thinner grid or reduce wavelet dominant frequency, but thinner grid means the calculated amount of magnanimity, reduce counting yield, storage and calculated performance are all proposed to challenge, reduced wavelet dominant frequency and can lose radio-frequency component, affected resolution.

The mode of current existing compacting numerical value frequency dispersion mainly contains:

1) utilize flux alignment technique (FCT) to suppress numerical value frequency dispersion, FCT technology has been applied in the 3-component earthquake numerical simulation in anisotropic medium;

2) adopt High-order Difference Methods to suppress numerical value frequency dispersion;

3) utilize optimization method, in the given limits of error, the finite difference coefficient of calculation optimization.For example: conventional optimization method is to utilize least square method, in the given limits of error, obtain larger spectrum coverage, calculation optimization coefficient; Up-to-date is to utilize global optimization method simulation finite difference to optimize operator, makes the optimization operator obtaining have higher spectrum coverage and the less limits of error.

4) in spatial domain, the space convolution sequence of utilizing window function to block pseudo-spectrometry is derived finite difference operator, and the proposition of this mode is mainly that to can be understood as based on finite difference method be the blocking of space convolution of pseudo-spectrometry.Pseudo-spectrometry, because adopted all points, has solved the problem of numerical value frequency dispersion.In other words,, in spatial domain, can adopt different truncated window to go the space convolution sequence of blocking pseudo-spectrometry to derive finite difference operator.Having plenty of at present the Hanning(Hamming that has used broad sense weighting) window blocks the finite difference operator being optimized, to use Gaussian window to block the finite difference operator being optimized in addition, also have plenty of the binomial of use window and unified limit difference method and pseudo-spectrometry, and provided the scheme that a kind of improved binomial window blocks optimization.

In above four kinds of optimization methods, adopt optimization method, in given approximate error limit, the method of trying to achieve optimized coefficients is the most conventional, but the coefficient that this optimization method need to be optimized is more, the selection of the limits of error and objective function is too large on result impact, and Optimizing Flow is not very directly perceived, normally directly calculates the difference coefficient that meets the limits of error.And window function method for cutting is a kind of visual method the most directly perceived, its essence is to adopt the window function space convolution sequence of blocking pseudo-spectrometry, and the precision of approaching as much as possible pseudo-spectrometry can change the parameter of window function at any time, obtains better approximation accuracy.But different window functions blocks different results, because it has different amplitude responses, the main lobe in window function amplitude response and side lobe performance have directly had influence on the precision of approaching, and have also just directly affected the accuracy of numerical simulation.Therefore the window function that, how to confirm is suitable seems and is even more important.

At present, conventional window function has the Hanning(Hamming of broad sense weighting) window, Gaussian window, improved binomial window etc., but a kind of window function of single selection is due to the character of window function self, in the situation that ensureing main lobe performance, being difficult to ensure the performance of secondary lobe, is not fine thereby make the effect of seismic wave field numerical simulation.For example, for the Hanning window of broad sense weighting, its essence is trigonometric function class window, trigonometric function class window amplitude response has the main lobe of relative narrower, but side lobe attenuation degree is not high enough, use this trigonometric function class window function to block the finite difference operator approaching, it is relatively large that it blocks the trueness error fluctuation approaching, even if there is larger spectrum coverage, but the fluctuation of trueness error is very large on the impact of the effect of approaching, improved binomial window is an adjustable window function, but its main lobe or wide, cause transitional zone long, use this improved binomial window function to block the finite difference operator approaching, its trueness error spectrum coverage is less, especially for low order finite difference operator (4, 8, 12 rank), there is larger fluctuation in its trueness error.Therefore,, when one of How to choose has narrower main lobe, can keep again the window function of the decay of secondary lobe to seem particularly important for the numerical simulation of seismic wave field simultaneously.

Summary of the invention

The embodiment of the present invention provides a kind of earthquake vector wave Numerical Simulation method, and to reach the object of Stability and veracity of effective raising numerical simulation result, the method comprises:

Original window function is done to L auto convolution computing and obtain the window function after auto convolution, wherein, L is positive integer;

Window function after auto convolution and original window function are computed weighted, obtain auto convolution combination window function;

Block to approach by auto convolution combination window function and obtain finite difference operator;

According to described finite difference operator, numerical simulation is carried out in earthquake vector wave field.

In one embodiment, described original window function comprises: Hamming window function, improved binomial window function or Chebyshev's window function.

In one embodiment, described Chebyshev's window function is:

$w_C\left(n\right)=\left\{\begin{array}{c}\frac{1}{N+1}\{\frac{1}{r}+2\underset{i=1}{\overset{N/2}{\mathrm{\Σ}}}{C}_N\left(x_0\cos \left(\frac{\mathrm{i\π}}{N+1}\right)\right)\cos \left(\frac{2\mathrm{ni\π}}{N+1}\right)\},\left|n\right|\≤N/2\\ 0,\left|n\right|>N/2\end{array}\right.$

Wherein, $C_N\left(x\right)=\left\{\begin{array}{c}\cos \left[N{\cos }^{-1}\left(x\right)\right],\left|x\right|\≤1\\ \cosh \left[N{\cosh }^{-1}\left(x\right)\right],\left|x\right|>1\end{array}\right.,x_0=\cosh \left[\frac{1}{N}\cos h^{-1}\left(\frac{1}{r}\right)\right],$ R represents ripple rate, and N+1 represents the length of window, and N is even number.

In one embodiment, said method also comprises:

Adjust the ripple rate of described Chebyshev's window function;

Determine under different ripple rates main lobe performance and the side lobe performance of the amplitude response of Chebyshev's window function;

Chebyshev's window function main lobe performance and side lobe performance being met under the ripple ratio of pre-provisioning request is defined as original window function.

In one embodiment, the window function after auto convolution and original window function are computed weighted, obtain auto convolution combination window function, comprising:

Determining need main lobe performance preferential in the situation that the weighting coefficient of the window function of the weighting coefficient of original window function after higher than auto convolution;

Determining that the weighting coefficient of the function after auto convolution is higher than the weighting coefficient of original window function need side lobe performance preferential in the situation that.

In one embodiment, block to approach obtaining finite difference operator by auto convolution combination window function, comprising:

By the auto convolution combination window function obtaining, pseudo-spectrometry space convolution sequence is carried out to truncation, obtain finite difference operator;

Draw approximate error curve according to the finite difference operator obtaining;

Determine frequency spectrum coverage and the approximation accuracy stability of described approximate error curve;

If frequency spectrum coverage and approximation accuracy stability do not meet predetermined requirement, change the value of L or change the weighting coefficient of ranking operation and redefine auto convolution combination window function, meet pre-provisioning request until definite auto convolution combination window function blocks frequency spectrum coverage and the approximation accuracy stability of approaching the finite difference operator obtaining.

The embodiment of the present invention also provides a kind of earthquake vector wave Numerical Simulation device, and to reach the object of Stability and veracity of effective raising numerical simulation result, this device comprises:

Auto convolution module, obtains the window function after auto convolution for original window function being done to L auto convolution computing, and wherein, L is positive integer;

Weighting block, computes weighted for the window function to after auto convolution and original window function, obtains auto convolution combination window function;

Truncation module, obtains finite difference operator for blocking to approach by auto convolution combination window function;

Numerical simulation module, for carrying out numerical simulation according to described finite difference operator to earthquake vector wave field.

In one embodiment, described original window function comprises: Hamming window function, improved binomial window function or Chebyshev's window function.

In one embodiment, described Chebyshev's window function is:

$w_C\left(n\right)=\left\{\begin{array}{c}\frac{1}{N+1}\{\frac{1}{r}+2\underset{i=1}{\overset{N/2}{\mathrm{\Σ}}}{C}_N\left(x_0\cos \left(\frac{\mathrm{i\π}}{N+1}\right)\right)\cos \left(\frac{2\mathrm{ni\π}}{N+1}\right)\},\left|n\right|\≤N/2\\ 0,\left|n\right|>N/2\end{array}\right.$

Wherein, $C_N\left(x\right)=\left\{\begin{array}{c}\cos \left[N{\cos }^{-1}\left(x\right)\right],\left|x\right|\≤1\\ \cosh \left[N{\cosh }^{-1}\left(x\right)\right],\left|x\right|>1\end{array}\right.,x_0=\cosh \left[\frac{1}{N}\cos h^{-1}\left(\frac{1}{r}\right)\right],$ R represents ripple rate, and N+1 represents the length of window, and N is even number.

In one embodiment, said apparatus also comprises:

Adjusting module, for adjusting the ripple rate of described Chebyshev's window function;

Determination module, for determining under different ripple rates, main lobe performance and the side lobe performance of the amplitude response of Chebyshev's window function;

Select module, be defined as original window function for Chebyshev's window function main lobe performance and side lobe performance being met under the ripple ratio of pre-provisioning request.

In one embodiment, determining need main lobe performance preferential in the situation that the weighting coefficient of the window function of the weighting coefficient of original window function after higher than auto convolution;

Determining that the weighting coefficient of the function after auto convolution is higher than the weighting coefficient of original window function need side lobe performance preferential in the situation that.

In one embodiment, described truncation module comprises:

Block unit, for the auto convolution combination window function by obtaining, pseudo-spectrometry space convolution sequence is carried out to truncation, obtain finite difference operator;

Drawing of Curve unit, for drawing approximate error curve according to the finite difference operator obtaining;

Determining unit, for determining frequency spectrum coverage and the approximation accuracy stability of described approximate error curve;

Change unit, be used for the in the situation that of frequency spectrum coverage and satisfied predetermined requirement of approximation accuracy stability, change the value of L or change the weighting coefficient of ranking operation and redefine auto convolution combination window function, meet pre-provisioning request until definite auto convolution combination window function blocks frequency spectrum coverage and the approximation accuracy stability of approaching the finite difference operator obtaining.

In embodiments of the present invention, selected original window function is carried out to auto convolution computing, auto convolution computing makes the main lobe degradation of original window function, side lobe performance improves, then according to current need to computing weighted to the window function after original window function and auto convolution, thereby obtain a window function meeting the demands, then block to approach by the combination window function finally obtaining and obtain finite difference operator, finally by the finite difference operator obtaining, numerical simulation is carried out in earthquake vector wave field.By the way, solve in prior art, window function blocks the finite difference operator approaching cannot take into account spectrum coverage and trueness error stability, thereby the finite difference earthquake vector wave Numerical Simulation result value frequency dispersion causing is excessive, the technical matters that simulation precision is not high, has reached the technique effect of effective raising earthquake vector wave Numerical Simulation precision.

Brief description of the drawings

Accompanying drawing described herein is used to provide a further understanding of the present invention, forms the application's a part, does not form limitation of the invention.In the accompanying drawings:

Fig. 1 is a kind of method flow diagram of the earthquake vector wave Numerical Simulation method of the embodiment of the present invention;

Fig. 2 is the method flow diagram of definite optimum window function of the embodiment of the present invention;

Fig. 3 is the structured flowchart of the earthquake vector wave Numerical Simulation device of the embodiment of the present invention.

Embodiment

For making the object, technical solutions and advantages of the present invention clearer, below in conjunction with embodiment and accompanying drawing, the present invention is described in further details.At this, exemplary embodiment of the present invention and explanation thereof are used for explaining the present invention, but not as a limitation of the invention.

Inventor considers, every kind of window function all exists some shortcomings of self, for example: trigonometric function class window amplitude response has the main lobe of relative narrower, but the attenuation degree of secondary lobe is not high enough; The attenuation degree of improved binomial window secondary lobe is enough high, but main lobe is relatively wide.But in order to obtain, numerical value frequency dispersion is less, numerical simulation result comparatively accurately, often needs a kind of main lobe narrower, and enough large window functions of the decay of secondary lobe.For above-mentioned this problem, a kind of seismic wave field method for numerical simulation has been proposed in embodiments of the present invention, as shown in Figure 1, comprise the following steps:

Step 101: original window function is done to L auto convolution computing and obtain the window function after auto convolution, wherein, L is positive integer;

Step 102: the window function after auto convolution and original window function are computed weighted, obtain auto convolution combination window function;

Step 103: block to approach obtaining finite difference operator by auto convolution combination window function;

Step 104: numerical simulation is carried out in earthquake vector wave field according to described finite difference operator.

In the above-described embodiments, selected original window function is carried out to auto convolution computing, auto convolution computing makes the main lobe degradation of original window function, side lobe performance improves, then according to current need to computing weighted to the window function after original window function and auto convolution, thereby obtain a window function meeting the demands, then block to approach by the combination window function finally obtaining and obtain finite difference operator, finally by the finite difference operator obtaining, numerical simulation is carried out in earthquake vector wave field.By the way, solve in prior art, window function blocks the finite difference operator approaching cannot take into account the problem of composing coverage and trueness error stability, thereby the finite difference earthquake vector wave Numerical Simulation result value frequency dispersion causing is excessive, the technical matters that simulation precision is not high, has reached the technique effect of effective raising earthquake vector wave Numerical Simulation precision.

In above-mentioned steps 102, original window function in the window function to after auto convolution and selection computes weighted, what obtain auto convolution combination window function and be mainly considering current needs is that main lobe performance is preferential or side lobe performance is preferential, the Properties Control of main lobe the scope of transitional zone, namely frequency spectrum coverage, main lobe is more narrow better comparatively speaking, blocks like this finite difference operator approaching, and low order optimization operator can reach the precision of common high-order operator.The performance of secondary lobe has determined that finite difference operator approaches the extent of deviation of differentiating operator, and side lobe attenuation is larger, and to approach the extent of deviation of differentiating operator less for finite difference operator.Become large although carry out the amplitude response side lobe attenuation of the window function after auto convolution, the width of main lobe has also increased, and this is undesirable.Our object is narrower in order to design in an amplitude response main lobe, the window function that side lobe attenuation is large.Therefore, adopt with auto convolution before window function be not weighted the scheme of combination, under the prerequisite that maintains suitable main lobe width, increase as far as possible side lobe attenuation.Determining need main lobe performance preferential in the situation that, namely need larger spectrum coverage, the weighting coefficient of the window function of the weighting coefficient of original window function after higher than auto convolution; Determining need side lobe performance preferential in the situation that, namely need trueness error stability higher, the weighting coefficient of the function after auto convolution is higher than the weighting coefficient of original window function.

In the concrete process of implementing, the performance of considering the auto convolution combination window function obtaining may be not so good, can before finite difference numerical simulation, the auto convolution combination window function obtaining be tested carrying out, for example: by auto convolution combination window function, pseudo-spectrometry space convolution sequence is carried out to truncation, obtain finite difference operator; Draw approximate error curve according to the finite difference operator obtaining; Determine frequency spectrum coverage and the approximation accuracy stability of described approximate error curve; If frequency spectrum coverage and approximation accuracy stability do not meet predetermined requirement, change the value of L or change the weighting coefficient of ranking operation and redefine auto convolution combination window function; If the stability of frequency spectrum coverage and approximation accuracy meets pre-provisioning request, using this auto convolution combination window function as selected auto convolution combination window function, use this window function to block to approach the finite difference operator obtaining numerical simulation is carried out in earthquake vector wave field.That is to say, first utilize window function to carry out truncation, obtain finite difference operator, then determine frequency spectrum coverage and the approximation accuracy stability of the approximate error curve of this finite difference operator, if can meet the demands, so just calculate according to this auto convolution combination window function, if can not meet the demands, so just change the value of L or change the weights of ranking operation and redefine an auto convolution combination window function, thereby obtain an auto convolution combination window function more meeting the demands.

For the original window function of above-mentioned use can comprise following one of at least: Hamming window function, improved binomial window function, Chebyshev's window function, but it should be noted that, above-mentioned three kinds of window functions are only preferred embodiments of the embodiment of the present invention, for the present invention is described better, can also adopt other similar window function, for example Kaiser window function etc., the present invention is not construed as limiting this.

The expression formula of above-mentioned Chebyshev (Chebyshev) window function is:

$w_C\left(n\right)=\left\{\begin{array}{c}\frac{1}{N+1}\{\frac{1}{r}+2\underset{i=1}{\overset{N/2}{\mathrm{\Σ}}}{C}_N\left(x_0\cos \left(\frac{\mathrm{i\π}}{N+1}\right)\right)\cos \left(\frac{2\mathrm{ni\π}}{N+1}\right)\},\left|n\right|\≤N/2\\ 0,\left|n\right|>N/2\end{array}\right.$

Wherein, $C_N\left(x\right)=\left\{\begin{array}{c}\cos \left[N{\cos }^{-1}\left(x\right)\right],\left|x\right|\≤1\\ \cosh \left[N{\cosh }^{-1}\left(x\right)\right],\left|x\right|>1\end{array}\right.,x_0=\cosh \left[\frac{1}{N}\cos h^{-1}\left(\frac{1}{r}\right)\right],$ R represents ripple rate, and N+1 represents the length of window, is even number, w _c(n) represent truncated window function.

In the parameter that changes Chebyshev (Chebyshev) window function, select in the situation of original window function, can select in the following manner original window function: adjust described Chebyshev's window function ripple rate; Determine under different ripple rates main lobe performance and the side lobe performance of the amplitude response of Chebyshev's window function; Main lobe performance and side lobe performance are met to the Chebyshev's window function under the ripple ratio of pre-provisioning request, be defined as the original window function of selecting.Because affecting Chebyshev window function performance parameter is mainly ripple rate, as long as adjust ripple rate, the performance of determining the window function in each ripple rate situation, the performance that just can determine the window function which ripple ratio is corresponding is best, thereby selects original window function.

Below in conjunction with a specific embodiment, the embodiment of the present invention is specifically described, but it should be noted that this specific embodiment is only in order to describe better the present invention, not form inappropriate limitation of the present invention.

Finite difference operator is derived to more clearly understand finite difference by sinc function interpolation theory below.According to the sampling theory of discrete signal, the continuous signal f (x) of a band limit can be by the signal f with a uniform sampling _nrebuild by sinc function interpolation:

$f\left(x\right)=\underset{n=-\∞}{\overset{\∞}{\mathrm{\Σ}}}{f}_n\frac{\sin \frac{\mathrm{\π}}{\mathrm{\Δx}}(x-\mathrm{n\Δx})}{\frac{\mathrm{\π}}{\mathrm{\Δx}}(x-\mathrm{n\Δx})}$ (formula 1)

Wherein, Δ x represents sampling interval,

represent cut-off wave number.

If the right and left of above-mentioned formula 1 is asked respectively to 1 order derivative and 2 order derivatives, and get the derivative value at x=0 place, can obtain:

$\frac{\∂f}{\∂x}{|}_{x=0}=\frac{1}{\mathrm{\Δx}}\underset{n=-\∞}{\overset{\∞}{\mathrm{\Σ}}}[-\frac{1}{n}\cos \left(\mathrm{n\π}\right)]f\left(n\right)$ (formula 2)

$\frac{{\∂}^{2}f}{\∂x^2}{|}_{x=0}=\frac{1}{\mathrm{\Δ}{x}^2}\underset{n=-\∞}{\overset{\∞}{\mathrm{\Σ}}}[-\frac{1}{n^2}\cos \left(\mathrm{n\π}\right)]f\left(n\right)$ (formula 3)

With a length be the window function that N+1 is ordered, wherein, N is even number, removes to block formula 2 and formula 3, obtains finite difference operator:

$\frac{\∂f}{\∂x}{|}_{x=0}=\frac{1}{\mathrm{\Δx}}\underset{n=-N/2}{\overset{N/2\∞}{\mathrm{\Σ}}}w\left(n\right)[-\frac{1}{n}\cos \left(\mathrm{n\π}\right)]f\left(n\right)$ (formula 4)

$\frac{{\∂}^{2}f}{\∂x^2}{|}_{x=0}=\frac{1}{\mathrm{\Δ}{x}^2}\underset{n=-N/2}{\overset{N/2}{\mathrm{\Σ}}}[-\frac{1}{n^2}\cos \left(\mathrm{n\π}\right)]f\left(n\right)$ (formula 5)

Represent truncated window function with w (n), for conventional finite difference operator:

$w_n^N=\left(\begin{array}{c}N\\ \frac{N}{2}+n\end{array}\right)/\left(\begin{array}{c}N\\ \frac{N}{2}\end{array}\right)$

Because n=0 is a singular point of formula 2 and formula 3, formula 2 and formula 3 can be expressed as:

$\frac{\∂f}{\∂x}{|}_{x=0}=\frac{1}{\mathrm{\Δx}}\underset{n=1}{\overset{\∞}{\mathrm{\Σ}}}{d}_n^1(f_n-f_{-n})$ (formula 6)

$\frac{{\∂}^{2}f}{\∂x^2}{|}_{x=0}=\frac{1}{\mathrm{\Δ}{x}^2}(d_0^2{f}_0+\underset{n=1}{\overset{\∞}{\mathrm{\Σ}}}{d}_n^2(f_n+f_{-n}))$ (formula 7)

Wherein, in formula 6 and formula 7 with

for the f in formula 2 and formula 3 _ncoefficient:

$d_n^1=-\frac{1}{n}\cos \left(\mathrm{n\π}\right),n=\mathrm{1,2,3}...,\∞$

$d_n^2=-\frac{2}{n^2}\cos \left(\mathrm{n\π}\right),n=\mathrm{1,2,3}...,\∞$

$d_0^2=-\underset{n=1}{\overset{\∞}{\mathrm{\Σ}}}(d_n^2+d_{-n}^2)=4*(-\mathrm{\ζ}\left(2\right)/2)=-\frac{{\mathrm{\π}}^2}{3}$

Wherein, ζ represents Riemann zeta function, after windowing w (n) blocks, just has:

$\frac{\∂f}{\∂x}{|}_{x=0}=\frac{1}{\mathrm{\Δx}}\underset{n=1}{\overset{N/2}{\mathrm{\Σ}}}{b}_n(f_n-f_{-n})$ (formula 8)

$\frac{{\∂}^{2}f}{\∂x^2}{|}_{x=0}=\frac{1}{\mathrm{\Δ}{x}^2}(c_0{f}_0+\underset{n=1}{\overset{N/2}{\mathrm{\Σ}}}{c}_n(f_n+f_{-n}))$ (formula 9)

In above-mentioned formula:

$b_n=d_n^{1}w\left(n\right)=-\frac{1}{n}\cos \left(\mathrm{n\π}\right)w\left(n\right),n=\mathrm{1,2}...,N/2$

$c_n=d_n^{2}w\left(n\right)=-\frac{2}{n^2}\cos \left(\mathrm{n\π}\right)w\left(n\right),n=\mathrm{1,2}...,N/2$

$c_0=-2\underset{n=1}{\overset{N/2}{\mathrm{\Σ}}}(c_n+c_{-n}).$

Approach the precision of differentiating operator in order to optimize finite difference operator, can select Chebyshev window function to block, wherein, Chebyshev window function is expressed as follows:

$w_C\left(n\right)=\left\{\begin{array}{c}\frac{1}{N+1}\{\frac{1}{r}+2\underset{i=1}{\overset{N/2}{\mathrm{\Σ}}}{C}_N\left(x_0\cos \left(\frac{\mathrm{i\π}}{N+1}\right)\right)\cos \left(\frac{2\mathrm{ni\π}}{N+1}\right)\},\left|n\right|\≤N/2\\ 0,\left|n\right|>N/2\end{array}\right.$ (formula 10)

In above-mentioned formula 10:

$C_N\left(x\right)=\left\{\begin{array}{c}\cos \left[N{\cos }^{-1}\left(x\right)\right],\left|x\right|\≤1\\ \cosh \left[N{\cosh }^{-1}\left(x\right)\right],\left|x\right|>1\end{array}\right.$

r represents ripple rate, represents the attenuation degree of secondary lobe, C _n(x) represent Chebyshev polynomial expression.

In order to understand better the present invention, main lobe and the impact of side lobe performance on finite difference trueness error of window function amplitude response are first described:

1) main lobe size is relevant with transition band width, and main lobe is narrow, and transition band width is little, and frequency spectrum coverage is large, and main lobe is wide, and transitional zone is roomy, and frequency spectrum coverage is little.Frequency spectrum coverage is less, illustrate for low frequency part, finite difference approximation differentiating operator precision is high, for HFS, approximation accuracy sharply declines, only have frequency spectrum coverage large, the precision that the window function of low order blocks the finite difference operator approaching just can reach the precision of common high-order limited difference operator, and this also just need to have a narrower main lobe.

2) relation of side lobe attenuation and approximation accuracy stability, the decay of secondary lobe has directly had influence on the stability of approximation accuracy, and side lobe attenuation is larger, and stablizing of approximation accuracy is higher, and side lobe attenuation is less, and the stability of approximation accuracy is poorer.The stability of approximation accuracy can be understood as in whole spectral limit, precision curve is around the situation of null value fluctuation, fluctuation Shaoxing opera is strong, stability is poorer, fluctuates less, and stability is high, if approximation accuracy stability is bad, even if spectrum coverage is very large, also can cause the result of numerical simulation undesirable, therefore need to have a secondary lobe that decay is larger.

Therefore, taking above-mentioned Chebyshev's window function as example, can obtain in such a way a main lobe narrower, the window function that side lobe attenuation is larger, as shown in Figure 2, comprises the following steps:

Step 201: the main lobe of more multiple existing window function amplitude response and side lobe performance, select optimum window function.

Step 202: consider that window function auto convolution can cause side lobe attenuation performance boost, also can cause that accordingly main lobe broadens, performance dies down, therefore the optimum window function of selecting in step 201 is to auto convolution L time, produces a new window function, for example, L can select 2, be exactly to carry out auto convolution twice so, the convolution operation of spatial domain is exactly the product calculation of frequency domain, can be after frequency domain multiplies each other again inversefouriertransform obtain the window function after auto convolution.

Step 203: the optimum window function of selecting in the new window function producing in step 202 and step 201 is weighted to combination, obtain auto convolution combination window function, in the process that is weighted combination, the weighting coefficient λ selecting selects to obtain from 0 to 1 interval, can be by selecting λ, determine that main lobe performance is preferential or side lobe performance is preferential,, if wish that main lobe performance is preferential, so just by fixed the coefficient of window function of selecting in step 201 larger, otherwise, just by fixed the coefficient of the auto convolution window function obtaining in step 202 larger.

Step 204: the window function producing in applying step 203 removes to block the finite difference operator that sampled signal is optimized, and introduce approximate error function, calculate and draw approximate error curve, determine the frequency spectrum coverage of approximate error curve and the stability of approximation accuracy, if effect is bad, return to step 202 or 203, re-start circulation, until obtain satisfied result, otherwise execution step 205.

Step 205: obtain preferred window function.

In this example, auto convolution combination window function is ensureing, under the prerequisite of side lobe attenuation, to have narrower main lobe, and such amplitude response characteristic, can ensure that approximation accuracy has larger spectral limit, and trueness error fluctuation is less, and stability is stronger.In the process of implementing, only need to regulate the intrinsic parameter ripple of Chebyshev's window function rate r, auto convolution number of times L and weighted array coefficient lambda, just can obtain good approximate error curve.For example, for Chebyshev's window function, after computing repeatedly and experiment, determine and can select r=60, L=2, λ=0.89, the Chebyshev auto convolution composite window obtaining by said method blocks the finite difference operator optimization method approaching, than conventional finite difference operator, have higher precision, draw the following conclusions through experiment: the precision that the Chebyshev auto convolution composite window on 8 rank blocks the finite difference operator approaching has exceeded the conventional finite difference operator in 12 rank, the precision that the Chebyshev auto convolution composite window on 12 rank blocks the finite difference operator approaching has exceeded the conventional finite difference operator in 24 rank, the trueness error of in addition, blocking by this auto convolution composite window the finite difference operator approaching is within 0.0004.Meanwhile, in this example, by changing 3 parameters, can the more free and precision of blocking the finite difference operator approaching that changes visually directly perceived.

It should be noted that, only the description of carrying out as a specific embodiment using Chebyshev window function in the above-described embodiments, can also adopt other window function to carry out computing, in principle, according to main lobe and sidelobe performance on blocking the impact that approaches finite difference operator trueness error, other higher window functions of approximation accuracy after selecting to block.

Based on same inventive concept, in the embodiment of the present invention, also provide a kind of earthquake vector wave Numerical Simulation device, as described in the following examples.Because the principle that earthquake vector wave Numerical Simulation device is dealt with problems is similar to earthquake vector wave Numerical Simulation method, therefore the enforcement of earthquake vector wave Numerical Simulation device can be referring to the enforcement of earthquake vector wave Numerical Simulation method, repeats part and repeat no more.Following used, term " unit " or " module " can realize the combination of software and/or the hardware of predetermined function.Although the described device of following examples is preferably realized with software, hardware, or the realization of the combination of software and hardware also may and be conceived.Fig. 3 is a kind of structured flowchart of the earthquake vector wave Numerical Simulation device of the embodiment of the present invention, as shown in Figure 3, comprising: auto convolution module 301, weighting block 302 and numerical simulation module 303, describe this structure below.

Auto convolution module 301, obtains the window function after auto convolution for original window function being done to L auto convolution computing, and wherein, L is positive integer;

Weighting block 302, computes weighted for the window function to after auto convolution and original window function, obtains auto convolution combination window function;

Truncation module 303, obtains finite difference operator for blocking to approach by auto convolution combination window function;

Numerical simulation module 304, for carrying out numerical simulation according to described finite difference operator to earthquake vector wave field.

In one embodiment, above-mentioned original window function comprises: Hamming window function, improved binomial window function or Chebyshev's window function.

In one embodiment, above-mentioned Chebyshev's window function is:

$w_C\left(n\right)=\left\{\begin{array}{c}\frac{1}{N+1}\{\frac{1}{r}+2\underset{i=1}{\overset{N/2}{\mathrm{\Σ}}}{C}_N\left(x_0\cos \left(\frac{\mathrm{i\π}}{N+1}\right)\right)\cos \left(\frac{2\mathrm{ni\π}}{N+1}\right)\},\left|n\right|\≤N/2\\ 0,\left|n\right|>N/2\end{array}\right.$

$C_N\left(x\right)=\left\{\begin{array}{c}\cos \left[N{\cos }^{-1}\left(x\right)\right],\left|x\right|\≤1\\ \cosh \left[N{\cosh }^{-1}\left(x\right)\right],\left|x\right|>1\end{array}\right.,x_0=\cosh \left[\frac{1}{N}\cos h^{-1}\left(\frac{1}{r}\right)\right],$ R represents ripple rate, and N+1 represents the length of window, and N is even number.

In one embodiment, above-mentioned seismic wave field numerical simulation device also comprises: adjusting module, for adjusting the ripple rate of described Chebyshev's window function; Determination module, for determining under different ripple rates, main lobe performance and the side lobe performance of the amplitude response of Chebyshev's window function; Select module, be defined as original window function for Chebyshev's window function main lobe performance and side lobe performance being met under the ripple ratio of pre-provisioning request.

In one embodiment, determining need main lobe performance preferential in the situation that the weighting coefficient of the window function of the weighting coefficient of original window function after higher than auto convolution; Determining that the weighting coefficient of the function after auto convolution is higher than the weighting coefficient of original window function need side lobe performance preferential in the situation that.

In one embodiment, truncation module 303 comprises: block unit, for the auto convolution combination window function by obtaining, pseudo-spectrometry space convolution sequence is carried out to truncation, obtain finite difference operator; Drawing of Curve unit, for drawing approximate error curve according to the finite difference operator obtaining; Determining unit, for determining frequency spectrum coverage and the approximation accuracy stability of described approximate error curve; Change unit, be used for the in the situation that of frequency spectrum coverage and satisfied predetermined requirement of approximation accuracy stability, change the value of L or change the weighting coefficient of ranking operation and redefine auto convolution combination window function, meet pre-provisioning request until definite auto convolution combination window function blocks frequency spectrum coverage and the approximation accuracy stability of approaching the finite difference operator obtaining.

In another embodiment, also provide a kind of software, the technical scheme that this software is described for carrying out above-described embodiment and preferred implementation.

In another embodiment, also provide a kind of storage medium, stored above-mentioned software in this storage medium, this storage medium includes but not limited to: CD, floppy disk, hard disk, scratch pad memory etc.

From above description, can find out, the embodiment of the present invention has realized following technique effect: selected original window function is carried out to auto convolution computing, auto convolution computing makes the main lobe degradation of original window function, side lobe performance improves, then according to current need to computing weighted to the window function after original window function and auto convolution, thereby obtain a window function meeting the demands, then block the finite difference operator approaching by the combination window function finally obtaining numerical simulation is carried out in earthquake vector wave field.By the way, solve in prior art, window function blocks the finite difference operator approaching cannot take into account the problem of composing coverage and trueness error stability, thereby the finite difference earthquake vector wave Numerical Simulation result value frequency dispersion causing is excessive, the technical matters that simulation precision is not high, has reached the technique effect of effective raising earthquake vector wave Numerical Simulation precision.

Obviously, those skilled in the art should be understood that, each module of the above-mentioned embodiment of the present invention or each step can realize with general calculation element, they can concentrate on single calculation element, or be distributed on the network that multiple calculation elements form, alternatively, they can be realized with the executable program code of calculation element, thereby, they can be stored in memory storage and be carried out by calculation element, and in some cases, can carry out shown or described step with the order being different from herein, or they are made into respectively to each integrated circuit modules, or the multiple modules in them or step are made into single integrated circuit module to be realized.Like this, the embodiment of the present invention is not restricted to any specific hardware and software combination.

The foregoing is only the preferred embodiments of the present invention, be not limited to the present invention, for a person skilled in the art, the embodiment of the present invention can have various modifications and variations.Within the spirit and principles in the present invention all, any amendment of doing, be equal to replacement, improvement etc., within all should being included in protection scope of the present invention.

Claims (12)

1. an earthquake vector wave Numerical Simulation method, is characterized in that, comprising:

Original window function is done to L auto convolution computing and obtain the window function after auto convolution, wherein, L is positive integer;

Window function after auto convolution and original window function are computed weighted, obtain auto convolution combination window function;

Block to approach by auto convolution combination window function and obtain finite difference operator;

According to described finite difference operator, numerical simulation is carried out in earthquake vector wave field.

2. the method for claim 1, is characterized in that, described original window function comprises: Hamming window function, improved binomial window function or Chebyshev's window function.

3. method as claimed in claim 2, is characterized in that, described Chebyshev's window function is:

$w_C\left(n\right)=\left\{\begin{array}{c}\frac{1}{N+1}\{\frac{1}{r}+2\underset{i=1}{\overset{N/2}{\mathrm{\Σ}}}{C}_N\left(x_0\cos \left(\frac{\mathrm{i\π}}{N+1}\right)\right)\cos \left(\frac{2\mathrm{ni\π}}{N+1}\right)\},\left|n\right|\≤N/2\\ 0,\left|n\right|>N/2\end{array}\right.$

Wherein, $C_N\left(x\right)=\left\{\begin{array}{c}\cos \left[N{\cos }^{-1}\left(x\right)\right],\left|x\right|\≤1\\ \cosh \left[N{\cosh }^{-1}\left(x\right)\right],\left|x\right|>1\end{array}\right.,x_0=\cosh \left[\frac{1}{N}\cos h^{-1}\left(\frac{1}{r}\right)\right],$ R represents ripple rate, and N+1 represents the length of window, and N is even number.

4. method as claimed in claim 3, is characterized in that, also comprises:

Adjust the ripple rate of described Chebyshev's window function;

Determine under different ripple rates main lobe performance and the side lobe performance of the amplitude response of Chebyshev's window function;

Chebyshev's window function main lobe performance and side lobe performance being met under the ripple ratio of pre-provisioning request is defined as original window function.

5. the method for claim 1, is characterized in that, the window function after auto convolution and original window function are computed weighted, and obtains auto convolution combination window function, comprising:

Determining need main lobe performance preferential in the situation that the weighting coefficient of the window function of the weighting coefficient of original window function after higher than auto convolution;

Determining that the weighting coefficient of the function after auto convolution is higher than the weighting coefficient of original window function need side lobe performance preferential in the situation that.

6. the method as described in any one in claim 1 to 5, is characterized in that, blocks to approach obtaining finite difference operator by auto convolution combination window function, comprising:

By the auto convolution combination window function obtaining, pseudo-spectrometry space convolution sequence is carried out to truncation, obtain finite difference operator;

Draw approximate error curve according to the finite difference operator obtaining;

Determine frequency spectrum coverage and the approximation accuracy stability of described approximate error curve;

If frequency spectrum coverage and approximation accuracy stability do not meet predetermined requirement, change the value of L or change the weighting coefficient of ranking operation and redefine auto convolution combination window function, meet pre-provisioning request until definite auto convolution combination window function blocks frequency spectrum coverage and the approximation accuracy stability of approaching the finite difference operator obtaining.

7. an earthquake vector wave Numerical Simulation device, is characterized in that, comprising:

Auto convolution module, obtains the window function after auto convolution for original window function being done to L auto convolution computing, and wherein, L is positive integer;

Weighting block, computes weighted for the window function to after auto convolution and original window function, obtains auto convolution combination window function;

Truncation module, obtains finite difference operator for blocking to approach by auto convolution combination window function;

Numerical simulation module, for carrying out numerical simulation according to described finite difference operator to earthquake vector wave field.

8. device as claimed in claim 7, is characterized in that, described original window function comprises: Hamming window function, improved binomial window function or Chebyshev's window function.

9. device as claimed in claim 8, is characterized in that, described Chebyshev's window function is:

$w_C\left(n\right)=\left\{\begin{array}{c}\frac{1}{N+1}\{\frac{1}{r}+2\underset{i=1}{\overset{N/2}{\mathrm{\Σ}}}{C}_N\left(x_0\cos \left(\frac{\mathrm{i\π}}{N+1}\right)\right)\cos \left(\frac{2\mathrm{ni\π}}{N+1}\right)\},\left|n\right|\≤N/2\\ 0,\left|n\right|>N/2\end{array}\right.$

Wherein, $C_N\left(x\right)=\left\{\begin{array}{c}\cos \left[N{\cos }^{-1}\left(x\right)\right],\left|x\right|\≤1\\ \cosh \left[N{\cosh }^{-1}\left(x\right)\right],\left|x\right|>1\end{array}\right.,x_0=\cosh \left[\frac{1}{N}\cos h^{-1}\left(\frac{1}{r}\right)\right],$ R represents ripple rate, and N+1 represents the length of window, and N is even number.

10. device as claimed in claim 9, is characterized in that, also comprises:

Adjusting module, for adjusting the ripple rate of described Chebyshev's window function;

Determination module, for determining under different ripple rates, main lobe performance and the side lobe performance of the amplitude response of Chebyshev's window function;

Select module, be defined as original window function for Chebyshev's window function main lobe performance and side lobe performance being met under the ripple ratio of pre-provisioning request.

11. devices as claimed in claim 7, is characterized in that, are determining need main lobe performance preferential in the situation that the weighting coefficient of the window function of the weighting coefficient of original window function after higher than auto convolution;

Determining that the weighting coefficient of the function after auto convolution is higher than the weighting coefficient of original window function need side lobe performance preferential in the situation that.

12. devices as described in any one in claim 7 to 11, is characterized in that, described truncation module comprises:

Block unit, for the auto convolution combination window function by obtaining, pseudo-spectrometry space convolution sequence is carried out to truncation, obtain finite difference operator;

Drawing of Curve unit, for drawing approximate error curve according to the finite difference operator obtaining;

Determining unit, for determining frequency spectrum coverage and the approximation accuracy stability of described approximate error curve;

Change unit, be used for the in the situation that of frequency spectrum coverage and satisfied predetermined requirement of approximation accuracy stability, change the value of L or change the weighting coefficient of ranking operation and redefine auto convolution combination window function, meet pre-provisioning request until definite auto convolution combination window function blocks frequency spectrum coverage and the approximation accuracy stability of approaching the finite difference operator obtaining.

Images