CN103901472B - Frequency domain forward modeling method and device - Google Patents

Abstract

The invention relates to the technical field of geological exploration, in particular to a frequency domain forward modeling method and device. The method includes the steps that a 17-point-format difference formula is established, and a sparse matrix is constructed; wavelet parameters and a velocity model are read in; single frequency wave fields are acquired through frequency circulation; inverse Fourier transform of all the frequency wave fields is solved; the forward modeling result is obtained. With the method and device, the limitation that dx must be equal to dz in the prior art is overcome through fourth-order high-accuracy forward modeling, the method and device not only can be applied under the condition that dx is equal to dz but also can be applied under the condition that dx is not equal to dz and is more suitable for being applied in actual production, only 2.4 grids are needed for minimum wavelength with the method, and the method is superior to other similar methods.

Description

A kind of frequency field forward modeling method and device

Technical field

The present invention relates to technical field of geological exploration, is a kind of frequency field forward modeling method and device concretely.

Background technology

The Main Means of look for oil and gas has geophysical survey, oil geology method etc., and along with the progress of computer technology, geophysical survey more and more comes into one's own in oil gas is found, and more human and material resources drop into wherein.One of method method of seismic prospecting as geophysical survey has effect of crucial importance in oilfield prospecting developing, is that other geophysical methods are incomparable.

Method of seismic prospecting has a large amount of new technologies to emerge in large numbers every year, and wherein full waveform inversion method is exactly hot issue in recent years.This method can utilize the kinematics of seismic wave field and dynamic information to rebuild subsurface velocities structure, has under disclosing complicated geological background and constructs the potentiality with lithology detailed information.Full waveform inversion can be divided into time domain, frequency field and Laplace domain.Frequency field method enjoys favor with its following advantage: first frequency field method can simulate many big guns data simultaneously, and next is different from time domain method, and it only can carry out inverting to the data of a certain frequency, the 3rd, and this method does not exist cumulative errors.

It is the basis of frequency field full waveform inversion that frequency field is just being drilled, and only carry out frequency field and just drill, Inversion in frequency domain method just can be utilized effectively.In recent years, a lot of scholar is also had to propose variously just to drill difference method.

The key that frequency field is just being drilled balances the calculated amount just drilled precision and just drilling, and higher precision must cause higher calculated amount.Nineteen ninety Pratt and Worthington proposes classical five-point difference scheme in the literature the earliest, but net point required in this difference scheme minimum wavelength be 13 (within relative error 1%, in said minimum wavelength, required net point is all based on this standard below), cause calculated amount excessive, can not well just drill for frequency field.In order to reduce Grid dimension required in minimum wavelength, the people such as Jo publish an article the difference scheme proposed under a kind of 9 rotating coordinate systems newly for 1996, thus effectively reduce Grid dimension required in minimum wavelength, and then decrease calculator memory occupancy volume and Computing amount.Above two kinds of methods are all based on rotating coordinate system, result in them and have a defect, the situation that rate pattern dx (spatial sampling interval of horizontal direction) is not equal to dz (spatial sampling interval of depth direction) can not be applicable to exactly, and in actual exploration production, horizontal sampling interval is often not equal to depth direction sampling interval, and therefore this drawback must overcome.

Summary of the invention

The object of the invention is to solve prior art mainly concentrate on second order frequency territory and just drill, can not meet high precision and just drill demand, part high-order is just being drilled and there is the restriction of applicable elements,

Embodiments provide a kind of frequency field forward modeling method, comprise,

Set up the difference formula of 17 following dot formats,

Wherein

$\begin{array}{c}{\stackrel{\_}{P}}_{m+1,n}={\mathrm{\α}}_2{P}_{m+1,n}+\frac{1-{\mathrm{\α}}_2}{2}(P_{m+1,n+1}+P_{m+1,n+1}),\\ {\stackrel{\_}{P}}_{m-1,n}={\mathrm{\α}}_2{P}_{m-1,n}+\frac{1-{\mathrm{\α}}_1}{2}(P_{m-1,n+1}+P_{m-1,n-1}),\\ {\stackrel{\_}{P}}_{m+2,n}={\mathrm{\α}}_1{P}_{m+2,n}+\frac{1-{\mathrm{\α}}_1}{2}(P_{m+2,n+2}+P_{m+2,n-2}),\\ {\stackrel{\_}{P}}_{m-2,n}={\mathrm{\α}}_1{P}_{m-2,n}+\frac{1-{\mathrm{\α}}_1}{2}(P_{m-2,n+2}+P_{m-2,n-2}),\\ {\stackrel{\_}{P}}_{m,n}={\mathrm{\α}}_3{P}_{m,n}+{\mathrm{\α}}_4(P_{m,n+1}+P_{m,n-1})+\frac{1-2{\mathrm{\α}}_4-{\mathrm{\α}}_3}{2}(P_{m,n+2}+P_{m,n-2}),\end{array}---\left(3\right)$

With

$\begin{array}{c}{\tilde{P}}_{m,n+1}={\mathrm{\β}}_2{P}_{m,n+1}+\frac{1-{\mathrm{\β}}_2}{2}(P_{m+1,n+1}+P_{m-1,n+1}),\\ {\tilde{P}}_{m,n-1}={\mathrm{\β}}_2{P}_{m,n+1}+\frac{1-{\mathrm{\β}}_2}{2}(P_{m+1,n-1}+P_{m-1,n-1}),\\ {\tilde{P}}_{m,n+2}={\mathrm{\β}}_1{P}_{m,n+2}+\frac{1-{\mathrm{\β}}_1}{2}(P_{m+2,n+2}+P_{m-2,n+2}),\\ {\tilde{P}}_{m,n-2}={\mathrm{\β}}_1{P}_{m,n-2}+\frac{1-{\mathrm{\β}}_1}{2}(P_{m+2,n-2}+P_{m-2,n-2}),\\ {\tilde{P}}_{m,n}={\mathrm{\β}}_3{P}_{m,n}+{\mathrm{\β}}_4(P_{m+1,n}+P_{m-1,n})+\frac{1-2{\mathrm{\β}}_4-{\mathrm{\β}}_3}{2}(P_{m+2,n}+P_{m-2,n}),\end{array}---\left(4\right)$

Wherein P _m,n≈ P (m Δ z, n Δ x), Δ x, Δ z represent x direction sampling interval and z direction sampling interval, m and n represents the mesh coordinate in z direction and the mesh coordinate in x direction, α _i, β _i, b, c, d, e, f represent weighting coefficient when differentiate;

Structure sparse matrix;

Read in Wavelet parameter and rate pattern;

Frequency cycle obtains single-frequency wave field;

All frequency wave fields are negated Fourier transform;

Just drilled result.

The embodiment of the present invention additionally provides a kind of frequency field and just drills device, comprises,

Difference block, for setting up the difference formula of 17 following dot formats,

Wherein

$\begin{array}{c}{\stackrel{\_}{P}}_{m+1,n}={\mathrm{\α}}_2{P}_{m+1,n}+\frac{1-{\mathrm{\α}}_2}{2}(P_{m+1,n+1}+P_{m+1,n+1}),\\ {\stackrel{\_}{P}}_{m-1,n}={\mathrm{\α}}_2{P}_{m-1,n}+\frac{1-{\mathrm{\α}}_1}{2}(P_{m-1,n+1}+P_{m-1,n-1}),\\ {\stackrel{\_}{P}}_{m+2,n}={\mathrm{\α}}_1{P}_{m+2,n}+\frac{1-{\mathrm{\α}}_1}{2}(P_{m+2,n+2}+P_{m+2,n-2}),\\ {\stackrel{\_}{P}}_{m-2,n}={\mathrm{\α}}_1{P}_{m-2,n}+\frac{1-{\mathrm{\α}}_1}{2}(P_{m-2,n+2}+P_{m-2,n-2}),\\ {\stackrel{\_}{P}}_{m,n}={\mathrm{\α}}_3{P}_{m,n}+{\mathrm{\α}}_4(P_{m,n+1}+P_{m,n-1})+\frac{1-2{\mathrm{\α}}_4-{\mathrm{\α}}_3}{2}(P_{m,n+2}+P_{m,n-2}),\end{array}---\left(3\right)$

With

$\begin{array}{c}{\tilde{P}}_{m,n+1}={\mathrm{\β}}_2{P}_{m,n+1}+\frac{1-{\mathrm{\β}}_2}{2}(P_{m+1,n+1}+P_{m-1,n+1}),\\ {\tilde{P}}_{m,n-1}={\mathrm{\β}}_2{P}_{m,n+1}+\frac{1-{\mathrm{\β}}_2}{2}(P_{m+1,n-1}+P_{m-1,n-1}),\\ {\tilde{P}}_{m,n+2}={\mathrm{\β}}_1{P}_{m,n+2}+\frac{1-{\mathrm{\β}}_1}{2}(P_{m+2,n+2}+P_{m-2,n+2}),\\ {\tilde{P}}_{m,n-2}={\mathrm{\β}}_1{P}_{m,n-2}+\frac{1-{\mathrm{\β}}_1}{2}(P_{m+2,n-2}+P_{m-2,n-2}),\\ {\tilde{P}}_{m,n}={\mathrm{\β}}_3{P}_{m,n}+{\mathrm{\β}}_4(P_{m+1,n}+P_{m-1,n})+\frac{1-2{\mathrm{\β}}_4-{\mathrm{\β}}_3}{2}(P_{m+2,n}+P_{m-2,n}),\end{array}---\left(4\right)$

Wherein P _m,n≈ P (m Δ z, n Δ x), Δ x, Δ z represent x direction sampling interval and z direction sampling interval, m and n represents the mesh coordinate in z direction and the mesh coordinate in x direction, α _i, β _i, b, c, d, e, f represent weighting coefficient when differentiate;

Sparse matrix module, for constructing sparse matrix;

Read in module, for reading in Wavelet parameter and rate pattern;

Single-frequency wave field module, obtains single-frequency wave field for frequency cycle;

Inversefouriertransform module, for Fourier transform of negating to all frequency wave fields;

Result output module, is just drilled result.

By the method and apparatus of above-described embodiment, quadravalence high precision is just drilling the condition restriction overcoming and must equal dz in the past to dx, not only can be applicable to dx=dz, also the situation of dx ≠ dz is gone for, be more suitable for using in actual production, and the grid number required for minimum wavelength only needs 2.4 in the method, is also better than other similar approach.

Accompanying drawing explanation

Read the detailed description to embodiment in conjunction with the following drawings, above-mentioned feature and advantage of the present invention, and extra feature and advantage, will be more readily apparent from.

Figure 1 shows that the process flow diagram of a kind of frequency field forward modeling method of the embodiment of the present invention;

Figure 2 shows that a kind of frequency field of inventive embodiments is just drilling the structural drawing of device;

Figure 3 shows that the particular flow sheet of a kind of frequency field forward modeling method of the embodiment of the present invention;

Fig. 4 a and Fig. 4 b is depicted as embodiment of the present invention coefficient distribution diagram;

Table 1 is when Δ x >=Δ z, for the optimized coefficients table that different Δ x/ Δ z is corresponding;

Table 2 is depicted as when Δ x< Δ z, for the optimized coefficients table that different Δ x/ Δ z is corresponding;

Fig. 5 a is the distribution plan of rate pattern and recording geometry;

Fig. 5 b is the result schematic diagram just drilled;

Fig. 6 a is rate pattern and monitoring system arrangement figure thereof;

Fig. 6 b is single-frequency wave field schematic diagram;

Fig. 6 c be adopt the embodiment of the present invention obtain just drill the schematic diagram of wave field in time domain;

Design sketch shown in Fig. 6 d.

Embodiment

Description below can make any those skilled in the art utilize the present invention.Specific embodiment and the descriptor provided in applying are only example.Various extension and the combination of embodiment as described herein are apparent for those skilled in the art, and when not departing from the spirit and scope of the invention, the rule of the present invention's definition can be applied in other embodiments and application.Therefore, the present invention is not only limited to shown embodiment, and the maximum magnitude consistent with principle shown in this paper and feature is contained in the present invention.

Be illustrated in figure 1 the process flow diagram of a kind of frequency field forward modeling method of the embodiment of the present invention.

Comprise step 101, set up the difference formula of 17 following dot formats,

Wherein

$\begin{array}{c}{\stackrel{\_}{P}}_{m+1,n}={\mathrm{\&alpha;}}_2{P}_{m+1,n}+\frac{1-{\mathrm{\&alpha;}}_2}{2}(P_{m+1,n+1}+P_{m+1,n+1}),\\ {\stackrel{\_}{P}}_{m-1,n}={\mathrm{\&alpha;}}_2{P}_{m-1,n}+\frac{1-{\mathrm{\&alpha;}}_1}{2}(P_{m-1,n+1}+P_{m-1,n-1}),\\ {\stackrel{\_}{P}}_{m+2,n}={\mathrm{\&alpha;}}_1{P}_{m+2,n}+\frac{1-{\mathrm{\&alpha;}}_1}{2}(P_{m+2,n+2}+P_{m+2,n-2}),\\ {\stackrel{\_}{P}}_{m-2,n}={\mathrm{\&alpha;}}_1{P}_{m-2,n}+\frac{1-{\mathrm{\&alpha;}}_1}{2}(P_{m-2,n+2}+P_{m-2,n-2}),\\ {\stackrel{\_}{P}}_{m,n}={\mathrm{\&alpha;}}_3{P}_{m,n}+{\mathrm{\&alpha;}}_4(P_{m,n+1}+P_{m,n-1})+\frac{1-2{\mathrm{\&alpha;}}_4-{\mathrm{\&alpha;}}_3}{2}(P_{m,n+2}+P_{m,n-2}),\end{array}---\left(3\right)$

With

$\begin{array}{c}{\tilde{P}}_{m,n+1}={\mathrm{\&beta;}}_2{P}_{m,n+1}+\frac{1-{\mathrm{\&beta;}}_2}{2}(P_{m+1,n+1}+P_{m-1,n+1}),\\ {\tilde{P}}_{m,n-1}={\mathrm{\&beta;}}_2{P}_{m,n+1}+\frac{1-{\mathrm{\&beta;}}_2}{2}(P_{m+1,n-1}+P_{m-1,n-1}),\\ {\tilde{P}}_{m,n+2}={\mathrm{\&beta;}}_1{P}_{m,n+2}+\frac{1-{\mathrm{\&beta;}}_1}{2}(P_{m+2,n+2}+P_{m-2,n+2}),\\ {\tilde{P}}_{m,n-2}={\mathrm{\&beta;}}_1{P}_{m,n-2}+\frac{1-{\mathrm{\&beta;}}_1}{2}(P_{m+2,n-2}+P_{m-2,n-2}),\\ {\tilde{P}}_{m,n}={\mathrm{\&beta;}}_3{P}_{m,n}+{\mathrm{\&beta;}}_4(P_{m+1,n}+P_{m-1,n})+\frac{1-2{\mathrm{\&beta;}}_4-{\mathrm{\&beta;}}_3}{2}(P_{m+2,n}+P_{m-2,n}),\end{array}---\left(4\right)$

Wherein P _m,n≈ P (m Δ z, n Δ x), Δ x, Δ z represent x direction sampling interval and z direction sampling interval, m and n represents the mesh coordinate in z direction and the mesh coordinate in x direction, α _i, β _i, b, c, d, e, f represent weighting coefficient when differentiate.

Step 102, structure sparse matrix.

Step 103, reads in Wavelet parameter and rate pattern.

Step 104, frequency cycle obtains single-frequency wave field.

Step 105, to negate Fourier transform to all frequency wave fields.

Step 106, is just drilled result.

It is the basis of full waveform inversion that frequency field is just being drilled, it is an important step in full waveform inversion, and it is indispensable, full waveform inversion is the hot issue of current seismic prospecting, it can utilize the kinematics of seismic wave field and dynamic information to rebuild subsurface velocities structure, have under disclosing complicated geological background and construct the potentiality with lithology detailed information, thus more detailed structure and the distribution understanding stratum, be petroleum prospecting predicting hydrocarbon reservoirs.

Above-mentioned steps 102-step 106 is the conventional method that prior art medium frequency territory is just being drilled, and repeats no more in the application's literary composition.

17 dot format difference formulas of above-mentioned steps 101 can for the situation of dx ≠ dz, can also for the situation of dx=dz.

Also comprise after described step 101, ask for the optimal value of described weighting coefficient.

The described optimal value asking for weighting coefficient comprises, further by wave equation (1) set out, plane wave substitute into formula (1), arrangement can obtain dispersion equation, wherein, and P ₀for the initial amplitude of ripple, i are imaginary unit,

Wherein k representative is wave number,

As dx >=dz, definition r=Δ x/ Δ z, therefore above-mentioned formula (4) can be rewritten into following form

When dx >=dz, net point G required in minimum wavelength is defined as G=λ/Δ x, and wherein λ is wavelength,

So far, normalized phase velocity formula is

$\begin{array}{c}\frac{V_{\mathrm{ph}}}{v}=\frac{G}{2\mathrm{\&pi;}}{\left[\frac{E+r^{2}F}{b+2c(A+B)+4\mathrm{dAB}+2e(C+D)+4\mathrm{fCD}}\right]}^{1/2},\\ A=\cos \left(\frac{2\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{G}\right),\\ B=\cos \left(\frac{2\mathrm{\&pi;}\sin \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right),\\ C=\cos \left(\frac{4\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{G}\right),\\ D=\cos \left(\frac{4\mathrm{\&pi;}\sin \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right)\\ E=\frac{1}{6}{\mathrm{\&alpha;}}_1(C-\mathrm{CD})-\frac{8}{3}{\mathrm{\&alpha;}}_2(A-\mathrm{AB})+\frac{5}{2}{\mathrm{\&alpha;}}_3(1-D)+5{\mathrm{\&alpha;}}_4(B-D)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}D,\\ F=\frac{1}{6}{\mathrm{\&beta;}}_1(D-\mathrm{CD})-\frac{8}{3}{\mathrm{\&beta;}}_2(B-\mathrm{AB})+\frac{5}{2}{\mathrm{\&beta;}}_2(1-C)+5{\mathrm{\&beta;}}_4(A-C)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}C,\end{array}----\left(7\right)$

Wherein k _x=kcos (θ), k _z=ksin (θ),

Try to achieve the α of above-mentioned optimum _i, β _i, b, c, d, e, f, make V _phand difference is minimum between v;

As dx<dz, definition r=Δ z/ Δ x, G is defined as G=λ/Δ z;

Above-mentioned formula (4) can be rewritten as

Then normalized phase velocity formula is

$\begin{array}{c}\frac{V_{\mathrm{ph}}}{v}=\frac{G}{2\mathrm{\&pi;}}{\left[\frac{E+r^{2}F}{b+2c(A+B)+4\mathrm{dAB}+2e(C+D)+4\mathrm{fCD}}\right]}^{1/2},\\ A=\cos \left(\frac{2\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right),\\ B=\cos \left(\frac{2\mathrm{\&pi;}\sin \left(\mathrm{\&theta;}\right)}{G}\right),\\ C=\cos \left(\frac{4\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right),\\ D=\cos \left(\frac{4\mathrm{\&pi;}\sin \left(\mathrm{\&theta;}\right)}{G}\right)\\ E=\frac{1}{6}{\mathrm{\&alpha;}}_1(C-\mathrm{CD})-\frac{8}{3}{\mathrm{\&alpha;}}_2(A-\mathrm{AB})+\frac{5}{2}{\mathrm{\&alpha;}}_3(1-D)+5{\mathrm{\&alpha;}}_4(B-D)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}D,\\ F=\frac{1}{6}{\mathrm{\&beta;}}_1(D-\mathrm{CD})-\frac{8}{3}{\mathrm{\&beta;}}_2(B-\mathrm{AB})+\frac{5}{2}{\mathrm{\&beta;}}_2(1-C)+5{\mathrm{\&beta;}}_4(A-C)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}C,\end{array}----\left(7\right)$

Try to achieve the α of above-mentioned optimum _i, β _i, b, c, d, e, f, make V _phand difference is minimum between v.

Optimum α is asked above-mentioned _i, β _i, in the step of b, c, d, e, f, least square method can be adopted or gauss-newton method can also be used to try to achieve optimum α _i, β _i, b, c, d, e, f.

As one embodiment of the present of invention, also comprise after step 103 reads in Wavelet parameter and rate pattern, load perfectly matched layer PML boundary condition on the border of described rate pattern, for the wave field of absorbing boundary reflection, the wave field preventing border from reflecting is folded back to earth's surface.

As one embodiment of the present of invention, described step 102 comprises further, constructs sparse matrix in the following way:

c ₁P _m+1,n+c ₂P _m+1,n+1+c ₂P _m+1,n-1+c ₁P _m-1,n+c ₂P _m-1,n+1+c ₂P _m-1,n-1+c ₃P _m+2,n+c ₄P _m+2,n+2+ (8)c ₄P _m+2,n-2+c ₃P _m-2,n+c ₄P _m-2,n+2+c ₄P _m-2,n-2+c ₅P _m,n+c ₆P _m,n+1+c ₆P _m,n-1+c ₇P _m,n+2+c ₇P _m,n-2＝S,

Coefficient in formula is as follows

Each (m, n) a corresponding equation, total m*n equation, being focus item on the right of equation, be a line number is m*n, and columns is a matrix of 1 row, the left side is a line number m*n, columns is also a matrix of m*n, solves the single-frequency wave field that this system of equations just obtains each position of wave field, i.e. P _m,n.

By above-mentioned method, quadravalence high precision is just drilling the condition restriction overcoming and must equal dz in the past to dx, not only can be applicable to dx=dz, also the situation of dx ≠ dz is gone for, be more suitable for using in actual production, and the grid number required for minimum wavelength only needs 2.4 in the method, is also better than other similar approach.

Be illustrated in figure 2 the structural drawing that a kind of frequency field of inventive embodiments is just drilling device.

Comprise difference block 201, for setting up the difference formula of 17 following dot formats,

Wherein

$\begin{array}{c}{\stackrel{\_}{P}}_{m+1,n}={\mathrm{\&alpha;}}_2{P}_{m+1,n}+\frac{1-{\mathrm{\&alpha;}}_2}{2}(P_{m+1,n+1}+P_{m+1,n+1}),\\ {\stackrel{\_}{P}}_{m-1,n}={\mathrm{\&alpha;}}_2{P}_{m-1,n}+\frac{1-{\mathrm{\&alpha;}}_1}{2}(P_{m-1,n+1}+P_{m-1,n-1}),\\ {\stackrel{\_}{P}}_{m+2,n}={\mathrm{\&alpha;}}_1{P}_{m+2,n}+\frac{1-{\mathrm{\&alpha;}}_1}{2}(P_{m+2,n+2}+P_{m+2,n-2}),\\ {\stackrel{\_}{P}}_{m-2,n}={\mathrm{\&alpha;}}_1{P}_{m-2,n}+\frac{1-{\mathrm{\&alpha;}}_1}{2}(P_{m-2,n+2}+P_{m-2,n-2}),\\ {\stackrel{\_}{P}}_{m,n}={\mathrm{\&alpha;}}_3{P}_{m,n}+{\mathrm{\&alpha;}}_4(P_{m,n+1}+P_{m,n-1})+\frac{1-2{\mathrm{\&alpha;}}_4-{\mathrm{\&alpha;}}_3}{2}(P_{m,n+2}+P_{m,n-2}),\end{array}---\left(3\right)$

With

$\begin{array}{c}{\tilde{P}}_{m,n+1}={\mathrm{\&beta;}}_2{P}_{m,n+1}+\frac{1-{\mathrm{\&beta;}}_2}{2}(P_{m+1,n+1}+P_{m-1,n+1}),\\ {\tilde{P}}_{m,n-1}={\mathrm{\&beta;}}_2{P}_{m,n+1}+\frac{1-{\mathrm{\&beta;}}_2}{2}(P_{m+1,n-1}+P_{m-1,n-1}),\\ {\tilde{P}}_{m,n+2}={\mathrm{\&beta;}}_1{P}_{m,n+2}+\frac{1-{\mathrm{\&beta;}}_1}{2}(P_{m+2,n+2}+P_{m-2,n+2}),\\ {\tilde{P}}_{m,n-2}={\mathrm{\&beta;}}_1{P}_{m,n-2}+\frac{1-{\mathrm{\&beta;}}_1}{2}(P_{m+2,n-2}+P_{m-2,n-2}),\\ {\tilde{P}}_{m,n}={\mathrm{\&beta;}}_3{P}_{m,n}+{\mathrm{\&beta;}}_4(P_{m+1,n}+P_{m-1,n})+\frac{1-2{\mathrm{\&beta;}}_4-{\mathrm{\&beta;}}_3}{2}(P_{m+2,n}+P_{m-2,n}),\end{array}---\left(4\right)$

Wherein P _m,n≈ P (m Δ z, n Δ x), Δ x, Δ z represent x direction sampling interval and z direction sampling interval, m and n represents the mesh coordinate in z direction and the mesh coordinate in x direction, α _i, β _i, b, c, d, e, f represent weighting coefficient when differentiate.

Sparse matrix module 202, for constructing sparse matrix.

Read in module 203, for reading in Wavelet parameter and rate pattern.

Single-frequency wave field module 204, obtains single-frequency wave field for frequency cycle.

Inversefouriertransform module 205, for Fourier transform of negating to all frequency wave fields.

Result output module 206, is just drilled result.

The conventional method that above-mentioned sparse matrix module 202-result output module 206 all can adopt prior art medium frequency territory just drilling, repeats no more in the application's literary composition.

The 17 dot format difference formulas that above-mentioned difference block 201 obtains can for the situation of dx ≠ dz, can also for the situation of dx=dz.

Also comprising weighting coefficient module 207, being connected between described difference block 201 and sparse matrix module 202, for asking for the optimal value of described weighting coefficient.

The optimal value that described weighting coefficient module 207 asks for weighting coefficient comprises, further by wave equation (1) set out, plane wave substitute into formula (1), arrangement can obtain dispersion equation, and wherein, P0 is the initial amplitude of ripple, i is imaginary unit,

Wherein k representative is wave number,

As dx >=dz, definition r=Δ x/ Δ z, therefore above-mentioned formula (4) can be rewritten into following form

When dx >=dz, net point G required in minimum wavelength is defined as G=λ/Δ x, and wherein λ is wavelength,

So far, normalized phase velocity formula is

$\begin{array}{c}\frac{V_{\mathrm{ph}}}{v}=\frac{G}{2\mathrm{\&pi;}}{\left[\frac{E+r^{2}F}{b+2c(A+B)+4\mathrm{dAB}+2e(C+D)+4\mathrm{fCD}}\right]}^{1/2},\\ A=\cos \left(\frac{2\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{G}\right),\\ B=\cos \left(\frac{2\mathrm{\&pi;}\sin \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right),\\ C=\cos \left(\frac{4\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{G}\right),\\ D=\cos \left(\frac{4\mathrm{\&pi;}\sin \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right)\\ E=\frac{1}{6}{\mathrm{\&alpha;}}_1(C-\mathrm{CD})-\frac{8}{3}{\mathrm{\&alpha;}}_2(A-\mathrm{AB})+\frac{5}{2}{\mathrm{\&alpha;}}_3(1-D)+5{\mathrm{\&alpha;}}_4(B-D)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}D,\\ F=\frac{1}{6}{\mathrm{\&beta;}}_1(D-\mathrm{CD})-\frac{8}{3}{\mathrm{\&beta;}}_2(B-\mathrm{AB})+\frac{5}{2}{\mathrm{\&beta;}}_2(1-C)+5{\mathrm{\&beta;}}_4(A-C)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}C,\end{array}----\left(7\right)$

Wherein k _x=kcos (θ), k _z=ksin (θ),

Try to achieve the α of above-mentioned optimum _i, β _i, b, c, d, e, f, make V _phand difference is minimum between v;

As dx<dz, definition r=Δ z/ Δ x, G is defined as G=λ/Δ z;

Above-mentioned formula (4) can be rewritten as

Then normalized phase velocity formula is

$\begin{array}{c}\frac{V_{\mathrm{ph}}}{v}=\frac{G}{2\mathrm{\&pi;}}{\left[\frac{E+r^{2}F}{b+2c(A+B)+4\mathrm{dAB}+2e(C+D)+4\mathrm{fCD}}\right]}^{1/2},\\ A=\cos \left(\frac{2\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right),\\ B=\cos \left(\frac{2\mathrm{\&pi;}\sin \left(\mathrm{\&theta;}\right)}{G}\right),\\ C=\cos \left(\frac{4\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right),\\ D=\cos \left(\frac{4\mathrm{\&pi;}\sin \left(\mathrm{\&theta;}\right)}{G}\right)\\ E=\frac{1}{6}{\mathrm{\&alpha;}}_1(C-\mathrm{CD})-\frac{8}{3}{\mathrm{\&alpha;}}_2(A-\mathrm{AB})+\frac{5}{2}{\mathrm{\&alpha;}}_3(1-D)+5{\mathrm{\&alpha;}}_4(B-D)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}D,\\ F=\frac{1}{6}{\mathrm{\&beta;}}_1(D-\mathrm{CD})-\frac{8}{3}{\mathrm{\&beta;}}_2(B-\mathrm{AB})+\frac{5}{2}{\mathrm{\&beta;}}_2(1-C)+5{\mathrm{\&beta;}}_4(A-C)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}C,\end{array}----\left(7\right)$

Try to achieve the α of above-mentioned optimum _i, β _i, b, c, d, e, f, make V _phand difference is minimum between v.

Described weighting coefficient module 207 asks optimum α above-mentioned _i, β _i, in the process of b, c, d, e, f, least square method can be adopted or gauss-newton method can also be used to try to achieve optimum α _i, β _i, b, c, d, e, f.

As one embodiment of the present of invention, also comprise absorbing boundary module 208, read between module 203 and single-frequency wave field module 204 described in being connected to, for loading perfectly matched layer PML boundary condition on the border of described rate pattern, for the wave field of absorbing boundary reflection, the wave field preventing border from reflecting is folded back to earth's surface.

As one embodiment of the present of invention, described sparse matrix module 202 is also further used for, and constructs sparse matrix in the following way:

c ₁P _m+1,n+c ₂P _m+1,n+1+c ₂P _m+1,n-1+c ₁P _m-1,n+c ₂P _m-1,n+1+c ₂P _m-1,n-1+c ₃P _m+2,n+c ₄P _m+2,n+2+ (8)c ₄P _m+2,n-2+c ₃P _m-2,n+c ₄P _m-2,n+2+c ₄P _m-2,n-2+c ₅P _m,n+c ₆P _m,n+1+c ₆P _m,n-1+c ₇P _m,n+2+c ₇P _m,n-2＝S,

Coefficient in formula is as follows

Each (m, n) a corresponding equation, total m*n equation, being focus item on the right of equation, be a line number is m*n, and columns is a matrix of 1 row, the left side is a line number m*n, columns is also a matrix of m*n, solves the single-frequency wave field that this system of equations just obtains each position of wave field, i.e. P _m,n.

By above-mentioned device, quadravalence high precision is just drilling the condition restriction overcoming and must equal dz in the past to dx, not only can be applicable to dx=dz, also the situation of dx ≠ dz is gone for, be more suitable for using in actual production, and the grid number required for minimum wavelength only needs 2.4 in the method, is also better than other similar approach.

Be illustrated in figure 3 the particular flow sheet of a kind of frequency field forward modeling method of the embodiment of the present invention.

Comprise step 301, set up the difference formula of 17 following dot formats.

The wave equation of frequency field can be expressed as

P represents wave field, represent circular frequency, v is ripple velocity of propagation in media as well.

Utilize the method for finite difference method and mean derivative, the difference formula of 17 dot formats made new advances of deriving, as follows:

Wherein

$\begin{array}{c}{\stackrel{\_}{P}}_{m+1,n}={\mathrm{\&alpha;}}_2{P}_{m+1,n}+\frac{1-{\mathrm{\&alpha;}}_2}{2}(P_{m+1,n+1}+P_{m+1,n+1}),\\ {\stackrel{\_}{P}}_{m-1,n}={\mathrm{\&alpha;}}_2{P}_{m-1,n}+\frac{1-{\mathrm{\&alpha;}}_1}{2}(P_{m-1,n+1}+P_{m-1,n-1}),\\ {\stackrel{\_}{P}}_{m+2,n}={\mathrm{\&alpha;}}_1{P}_{m+2,n}+\frac{1-{\mathrm{\&alpha;}}_1}{2}(P_{m+2,n+2}+P_{m+2,n-2}),\\ {\stackrel{\_}{P}}_{m-2,n}={\mathrm{\&alpha;}}_1{P}_{m-2,n}+\frac{1-{\mathrm{\&alpha;}}_1}{2}(P_{m-2,n+2}+P_{m-2,n-2}),\\ {\stackrel{\_}{P}}_{m,n}={\mathrm{\&alpha;}}_3{P}_{m,n}+{\mathrm{\&alpha;}}_4(P_{m,n+1}+P_{m,n-1})+\frac{1-2{\mathrm{\&alpha;}}_4-{\mathrm{\&alpha;}}_3}{2}(P_{m,n+2}+P_{m,n-2}),\end{array}---\left(3\right)$

With

$\begin{array}{c}{\tilde{P}}_{m,n+1}={\mathrm{\&beta;}}_2{P}_{m,n+1}+\frac{1-{\mathrm{\&beta;}}_2}{2}(P_{m+1,n+1}+P_{m-1,n+1}),\\ {\tilde{P}}_{m,n-1}={\mathrm{\&beta;}}_2{P}_{m,n+1}+\frac{1-{\mathrm{\&beta;}}_2}{2}(P_{m+1,n-1}+P_{m-1,n-1}),\\ {\tilde{P}}_{m,n+2}={\mathrm{\&beta;}}_1{P}_{m,n+2}+\frac{1-{\mathrm{\&beta;}}_1}{2}(P_{m+2,n+2}+P_{m-2,n+2}),\\ {\tilde{P}}_{m,n-2}={\mathrm{\&beta;}}_1{P}_{m,n-2}+\frac{1-{\mathrm{\&beta;}}_1}{2}(P_{m+2,n-2}+P_{m-2,n-2}),\\ {\tilde{P}}_{m,n}={\mathrm{\&beta;}}_3{P}_{m,n}+{\mathrm{\&beta;}}_4(P_{m+1,n}+P_{m-1,n})+\frac{1-2{\mathrm{\&beta;}}_4-{\mathrm{\&beta;}}_3}{2}(P_{m+2,n}+P_{m-2,n}),\end{array}---\left(4\right)$

Wherein P _m,n≈ P (m Δ z, n Δ x), Δ x, Δ z represent x direction sampling interval and z direction sampling interval, the same with dx and the dz implication in literary composition, m and n represents the mesh coordinate in z direction and the mesh coordinate in x direction.α _i, β _i, b, c, d, e, f represent weighting coefficient when differentiate, and concrete coefficient distributes as shown in figures 4 a and 4b.

This new form goes for the situation of dx ≠ dz, adds the applicability of this difference formula.The situation of dx=dz is contained in derivation dx >=dz, and the Optimal Parameters also sought out in table 1 in dx=dz situation, frequency field wave equation difference formula under having had Optimal Parameters just to have dx=dz situation, the frequency field having had this difference formula just can do under dx=dz situation is just drilled.The frequency field in embodiment being below exactly the dx=dz situation made of the present invention is just drilled.

Step 302, asks for the weighting coefficient of optimization.

Based on this difference scheme, carry out frequency field sound wave and just drill, also need the problem solving frequency dispersion, little frequency dispersion can obtain just drills effect preferably.Therefore the derivation of frequency dispersion formula is carried out:

First by wave equation formula (1), plane wave substitute into formula (1), arrangement can obtain dispersion equation, wherein, and P ₀for the initial amplitude of ripple, i are imaginary unit

Wherein k representative is wave number, and other meaning of parameters are pointed out hereinbefore, repeat no more.

In two kinds of situation:

(1) as dx >=dz, definition r=Δ x/ Δ z, therefore above-mentioned formula (4) can be rewritten into following form

When dx >=dz, net point G required in minimum wavelength is defined as G=λ/Δ x (λ is wavelength),

So far, normalized phase velocity formula is

$\begin{array}{c}\frac{V_{\mathrm{ph}}}{v}=\frac{G}{2\mathrm{\&pi;}}{\left[\frac{E+r^{2}F}{b+2c(A+B)+4\mathrm{dAB}+2e(C+D)+4\mathrm{fCD}}\right]}^{1/2},\\ A=\cos \left(\frac{2\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{G}\right),\\ B=\cos \left(\frac{2\mathrm{\&pi;}\sin \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right),\\ C=\cos \left(\frac{4\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{G}\right),\\ D=\cos \left(\frac{4\mathrm{\&pi;}\sin \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right)\\ E=\frac{1}{6}{\mathrm{\&alpha;}}_1(C-\mathrm{CD})-\frac{8}{3}{\mathrm{\&alpha;}}_2(A-\mathrm{AB})+\frac{5}{2}{\mathrm{\&alpha;}}_3(1-D)+5{\mathrm{\&alpha;}}_4(B-D)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}D,\\ F=\frac{1}{6}{\mathrm{\&beta;}}_1(D-\mathrm{CD})-\frac{8}{3}{\mathrm{\&beta;}}_2(B-\mathrm{AB})+\frac{5}{2}{\mathrm{\&beta;}}_2(1-C)+5{\mathrm{\&beta;}}_4(A-C)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}C,\end{array}----\left(7\right)$

Wherein k _x=kcos (θ), k _z=ksin (θ)

Reduce frequency dispersion and find optimal coefficient α exactly _i, β _i, b, c, d, e, f, make V _phand difference is minimum between v, such as least square method can be utilized to try to achieve optimal value, the Gauss-Newton method that also can be proposed in document in 2000 by people such as Min.Thus obtain the α that in formula (2) (3) (4), different dx and dz is corresponding _i, β _i, b, c, d, e, f, specifically see table 1 (seeing below attached list) for when Δ x>=Δ z, for the optimized coefficients table that different Δ x/ Δ z is corresponding.

As dx<dz, r=Δ z/ Δ x, G (net point required in minimum wavelength) is defined as G=λ/Δ z (λ is wavelength), optimized coefficients under different dx and dz ratios when can obtain dx<dz, as table 2 (seeing below attached list) is depicted as when Δ x< Δ z, for the optimized coefficients table that different Δ x/ Δ z is corresponding.

Step 303, structure sparse matrix.

Structure sparse matrix, this step can be explained as follows with formula

c ₁P _m+1,n+c ₂P _m+1,n+1+c ₂P _m+1,n-1+c ₁P _m-1,n+c ₂P _m-1,n+1+c ₂P _m-1,n-1+c ₃P _m+2,n+c ₄P _m+2,n+2+ (8)c ₄P _m+2,n-2+c ₃P _m-2,n+c ₄P _m-2,n+2+c ₄P _m-2,n-2+c ₅P _m,n+c ₆P _m,n+1+c ₆P _m,n-1+c ₇P _m,n+2+c ₇P _m,n-2＝S,

Coefficient in formula is as follows

Above-mentioned formula (8) represents in fact a system of equations, each (m, n) corresponding equation, therefore total m*n equation.When given speed model and recording geometry, above-mentioned equation coefficient is determined, namely that required is P _m,n, being focus item on the right of equation, be a line number is m*n, and columns is a matrix of 1 row, and the left side is a line number m*n, and columns is also a matrix of m*n, solves the single-frequency wave field that this system of equations just obtains each position of wave field, i.e. P _{m, n}, LU decomposition to be used at this and separate this system of equations.

Wherein, it is be the product of a lower triangular matrix and a upper triangular matrix by a matrix decomposition that LU decomposes, and is mainly used in numerical analysis, is used for separating linear equation, the matrix or calculate determinant of negating, concrete formula and program with reference to prior art, can repeat no more in this application.

Step 304, reads in Wavelet parameter and rate pattern.Wavelet parameter is time domain, makes it transform to frequency field, obtain the amplitude that each frequency is corresponding by Fourier's change, the S (matrixes of 1 row) above correspondence in formula (8).

Step 305, loads PML (perfectly matched layer) boundary condition, and namely add the medium that one deck absorbs wave field on border, the wave field preventing border from reflecting is folded back to earth's surface.

Step 306, circulates from low to high to frequency, and the single-frequency wave field obtained that at every turn circulates all is stored in a matrix.

Step 307, the inversefouriertransform carrying out all frequencies to point each in wave field is just converted into the wave field of time domain.

Step 308, obtains forward record by the wave field record extracting earth's surface.

Lineups position in forward record and frequency dispersion situation is obtained in order to what see said method clearly, first use a simple stratified model, stratified model velocity structure is simple, be easy to analyze, the mesh spacing dx=dz=5m of this model, the Grid dimension of horizontal direction and depth direction is 100 (nx=nz=100), for such model, and formula (8) and (9) middle factor alpha _i, β _i, b, c, d, e, f adopt the parameter of r=1 in form 1.Focus is placed on the position of (250m, 20m), and wavelet adopts Ricker wavelet, and dominant frequency is 20Hz, and wave detector is positioned at the position of Z=20m, and time sampling interval is 2ms, and the receiving record time is 1s.Fig. 5 a is the distribution plan of rate pattern and recording geometry, Fig. 5 b is the result schematic diagram just drilled, wherein, Fig. 5 a represents the seismic wave propagation speed of this part region underground structure, and different colours represents different speed, as shown in right block, it is high that color is tending towards red representation speed value, the highest is 2000m/s, and blue representation speed is low, and minimum is 1000m/s; The inverted triangle of white represents the distribution of wave detector, point red in the middle of wave detector represents demolition point position, show that lineups position is correct by Fig. 5 a and Fig. 5 b, and border absorbs better, there is not interference wave, owing to have employed the difference scheme of quadravalence, therefore frequency dispersion is very weak, has suppressed frequency dispersion preferably.

For the validity of further verification algorithm, adopt the higher Marmousi rate pattern of seismic prospecting industry perceive degree just to drill, this model dx=12.5m, dz=4m, therefore dx ≠ dz, we adopt the figure parameters of r=dx/dz=3.125 in form 1.This model we get a wherein part and just drill, grid number is nx=nz=301, and wavelet adopts Ricker wavelet equally, and dominant frequency is 12Hz, is placed on x=500m, the place of z=16m, and wave detector is positioned at the place of z=16m.Fig. 6 a is rate pattern and monitoring system arrangement figure thereof, color in this figure is as described in above-mentioned Fig. 5 a, Fig. 6 b is single-frequency wave field schematic diagram, wherein redder part represents that amplitude is higher, inclined yl moiety represents that amplitude is more weak, redder explanation amplitude is higher, more yellow explanation amplitude is more weak, Fig. 6 c be adopt the embodiment of the present invention obtain just drill the schematic diagram of wave field in time domain, in order to verify its correctness, by the method that time domain is just being drilled, this model is just drilled, obtain the design sketch shown in Fig. 6 d, two figure lineups position consistency, illustrate that this patent method is correct.

By above-mentioned method and device, quadravalence high precision is just drilling the condition restriction overcoming and must equal dz in the past to dx, not only can be applicable to dx=dz, also the situation of dx ≠ dz is gone for, be more suitable for using in actual production, and the grid number required for minimum wavelength only needs 2.4 in the method, is also better than other similar approach.

The present invention can realize in any suitable form, comprises hardware, software, firmware or their combination in any.The present invention can according to circumstances selectively partly realize, and such as software performing is in one or more data processor and digital signal processor.Element and the assembly of each embodiment herein can realize physically, functionally, in logic in any suitable manner.In fact, a function can the part in separate unit, in one group of unit or as other functional units realize.Therefore, this system and method both can realize in separate unit, also can be distributed between different unit and processor physically and functionally.

Technician in the related art will recognize that, embodiments of the invention have many possible amendments and combination, although form is slightly different, still adopts identical fundamental mechanism and method.In order to the object explained, aforementioned description with reference to several specific embodiment.But above-mentioned illustrative discussion is not intended to exhaustive or limits the precise forms of inventing herein.Above, many modifications and variations are possible.Selected and described embodiment, in order to explain principle of the present invention and practical application thereof, in order to the amendment for application-specific, the distortion that enable those skilled in the art utilize the present invention and each embodiment best.

Claims (10)

1. a frequency field forward modeling method, is characterized in that comprising,

Set up the difference formula of 17 following dot formats,

Wherein

${\stackrel{\&OverBar;}{P}}_{m+1,n}={\mathrm{\&alpha;}}_2{P}_{m+1,n}+\frac{1-{\mathrm{\&alpha;}}_2}{2}\left(P_{m+1,n-1}\right),$

${\stackrel{\&OverBar;}{P}}_{m-1,n}={\mathrm{\&alpha;}}_2{P}_{m-1,n}+\frac{1-{\mathrm{\&alpha;}}_2}{2}(P_{m-1,n+1}+P_{m-1,n-1})---\left(3\right)$

${\stackrel{\&OverBar;}{P}}_{m+2,n}={\mathrm{\&alpha;}}_1{P}_{m+2,n}+\frac{1-{\mathrm{\&alpha;}}_1}{2}(P_{m+2,n+2}+P_{m+2,n-2}),$

${\stackrel{\&OverBar;}{P}}_{m-2,n}={\mathrm{\&alpha;}}_1{P}_{m-2,n}+\frac{1-{\mathrm{\&alpha;}}_1}{2}(P_{m-2,n+2}+P_{m-2,n-2}),$

${\stackrel{\&OverBar;}{P}}_{m,n}={\mathrm{\&alpha;}}_3{P}_{m,n}+{\mathrm{\&alpha;}}_4(P_{m,n+1}+P_{m,n-1})+\frac{1-{2\mathrm{\&alpha;}}_4-{\mathrm{\&alpha;}}_3}{2}(P_{m,n+2}+P_{m,n-2}),$

With

${\tilde{P}}_{m,n+1}={\mathrm{\&beta;}}_2{P}_{m,n+1}+\frac{1-{\mathrm{\&beta;}}_2}{2}(P_{m+1,n+1}+P_{m-1,n+1}),$

${\tilde{P}}_{m,n-1}={\mathrm{\&beta;}}_2{P}_{m,n-1}+\frac{1-{\mathrm{\&beta;}}_2}{2}(P_{m+1,n-1}+P_{m-1,n-1}),---\left(4\right)$

${\tilde{P}}_{m,n+2}={\mathrm{\&beta;}}_1{P}_{m,n+2}+\frac{1-{\mathrm{\&beta;}}_1}{2}(P_{m+2,n+2}+P_{m-2,n+2}),$

${\tilde{P}}_{m,n-2}={\mathrm{\&beta;}}_1{P}_{m,n-2}+\frac{1-{\mathrm{\&beta;}}_1}{2}(P_{m+2,n-2}+P_{m-2,n-2}),$

${\tilde{P}}_{m,n}={\mathrm{\&beta;}}_3{P}_{m,n}+{\mathrm{\&beta;}}_4(P_{m+1,n}+P_{m-1,n})+\frac{1-{2\mathrm{\&beta;}}_4-{\mathrm{\&beta;}}_3}{2}(P_{m+2,n}+P_{m-2,n}),$

Wherein P _m,n≈ P (m Δ z, n Δ x), Δ x, Δ z represent x direction sampling interval and z direction sampling interval, m and n represents the mesh coordinate in z direction and the mesh coordinate in x direction, α _i, β _i, b, c, d, e, f represent weighting coefficient when differentiate;

Structure sparse matrix;

Read in Wavelet parameter and rate pattern;

Frequency cycle obtains single-frequency wave field;

All frequency wave fields are negated Fourier transform;

Just drilled result.

2. a kind of frequency field forward modeling method according to claim 1, is characterized in that, above-mentioned set up the difference formula of 17 dot formats after also comprise, ask for the optimal value of described weighting coefficient.

3. a kind of frequency field forward modeling method according to claim 2, is characterized in that, described in ask for weighting coefficient optimal value comprise further, by wave equation (1) set out, plane wave substitute into formula (1), arrangement can obtain dispersion equation, wherein, and P ₀for the initial amplitude of ripple, i are imaginary unit,

A＝cos(k _xΔx),

B＝cos(k _zΔz), （5）

C＝cos(2k _xΔx),

D＝cos(2k _zΔz),

$E=\frac{1}{6}{\mathrm{\&alpha;}}_1(C-\mathrm{CD})-\frac{8}{3}{\mathrm{\&alpha;}}_2(A-\mathrm{AB})+\frac{5}{2}{\mathrm{\&alpha;}}_3(1-D)+{5\mathrm{\&alpha;}}_4(B-D)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}D,$

$F=\frac{1}{6}{\mathrm{\&beta;}}_1(D-\mathrm{CD})-\frac{8}{3}{\mathrm{\&beta;}}_2(B-\mathrm{AB})+\frac{5}{2}{\mathrm{\&beta;}}_3(1-C)+{5\mathrm{\&beta;}}_4(A-C)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}C,$

Wherein k representative is wave number,

As dx >=dz, definition r=Δ x/ Δ z, therefore above-mentioned formula (4) is rewritten into following form

A＝cos(k _xΔx),

B＝cos(k _zΔz),

C＝cos(2k _xΔx), （6）

D＝cos(2k _zΔz),

$E=\frac{1}{6}{\mathrm{\&alpha;}}_1(C-\mathrm{CD})-\frac{8}{3}{\mathrm{\&alpha;}}_2(A-\mathrm{AB})+\frac{5}{2}{\mathrm{\&alpha;}}_3(1-D)+{5\mathrm{\&alpha;}}_4(B-D)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}D,$

$F=\frac{1}{6}{\mathrm{\&beta;}}_1(D-\mathrm{CD})-\frac{8}{3}{\mathrm{\&beta;}}_2(B-\mathrm{AB})+\frac{5}{2}{\mathrm{\&beta;}}_3(1-C)+{5\mathrm{\&beta;}}_4(A-C)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}C,$

r＝Δx/Δz,

When dx >=dz, net point G required in minimum wavelength is defined as G=λ/Δ x, and wherein λ is wavelength,

So far, normalized phase velocity formula is

$\frac{V_{\mathrm{ph}}}{v}=\frac{G}{2\mathrm{\&pi;}}{\left[\frac{E+r^{2}F}{b+2c(A+B)+4\mathrm{dAB}+2e(C+D)+4\mathrm{fCD}}\right]}^{1/2},$

$A=\cos \left(\frac{2\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{G}\right)$

$B=\cos \left(\frac{2\mathrm{\&pi;}\sin \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right),---\left(7\right)$

$C=\cos \left(\frac{4\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{G}\right),$

$D=\cos \left(\frac{4\mathrm{\&pi;}\sin \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right)$

$E=\frac{1}{6}{\mathrm{\&alpha;}}_1(C-\mathrm{CD})-\frac{8}{3}{\mathrm{\&alpha;}}_2(A-\mathrm{AB})+\frac{5}{2}{\mathrm{\&alpha;}}_3(1-D)+{5\mathrm{\&alpha;}}_4(B-D)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}D,$

$F=\frac{1}{6}{\mathrm{\&beta;}}_1(D-\mathrm{CD})-\frac{8}{3}{\mathrm{\&beta;}}_2(B-\mathrm{AB})+\frac{5}{2}{\mathrm{\&beta;}}_3(1-C)+{5\mathrm{\&beta;}}_4(A-C)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}C,$

Wherein

Try to achieve the α of above-mentioned optimum _i, β _i, b, c, d, e, f, make V _phand difference is minimum between v;

As dx<dz, definition r=Δ z/ Δ x, G is defined as G=λ/Δ z;

Above-mentioned formula (4) is rewritten as

A＝cos(k _xΔx),

B＝cos(k _zΔz),

C＝cos(2k _xΔx),

D＝cos(2k _zΔz),

r＝Δz/Δx

$E=\frac{1}{6}{\mathrm{\&alpha;}}_1(C-\mathrm{CD})-\frac{8}{3}{\mathrm{\&alpha;}}_2(A-\mathrm{AB})+\frac{5}{2}{\mathrm{\&alpha;}}_3(1-D)+{5\mathrm{\&alpha;}}_4(B-D)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}D,$

$F=\frac{1}{6}{\mathrm{\&beta;}}_1(D-\mathrm{CD})-\frac{8}{3}{\mathrm{\&beta;}}_2(B-\mathrm{AB})+\frac{5}{2}{\mathrm{\&beta;}}_3(1-C)+{5\mathrm{\&beta;}}_4(A-C)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}C,$

Then normalized phase velocity formula is

$\frac{V_{\mathrm{ph}}}{v}=\frac{G}{2\mathrm{\&pi;}}{\left[\frac{E+r^{2}F}{b+2c(A+B)+4\mathrm{dAB}+2e(C+D)+4\mathrm{fCD}}\right]}^{1/2},$

$A=\cos \left(\frac{2\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right)$

$B=\cos \left(\frac{2\mathrm{\&pi;}\sin \left(\mathrm{\&theta;}\right)}{G}\right),$

$C=\cos \left(\frac{4\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right),$

$D=\cos \left(\frac{4\mathrm{\&pi;}\sin \left(\mathrm{\&theta;}\right)}{G}\right)$

$E=\frac{1}{6}{\mathrm{\&alpha;}}_1(C-\mathrm{CD})-\frac{8}{3}{\mathrm{\&alpha;}}_2(A-\mathrm{AB})+\frac{5}{2}{\mathrm{\&alpha;}}_3(1-D)+{5\mathrm{\&alpha;}}_4(B-D)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}D,$

$F=\frac{1}{6}{\mathrm{\&beta;}}_1(D-\mathrm{CD})-\frac{8}{3}{\mathrm{\&beta;}}_2(B-\mathrm{AB})+\frac{5}{2}{\mathrm{\&beta;}}_3(1-C)+{5\mathrm{\&beta;}}_4(A-C)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}C,$

Try to achieve the α of above-mentioned optimum _i, β _i, b, c, d, e, f, make V _phand difference is minimum between v.

4. a kind of frequency field forward modeling method according to claim 3, is characterized in that, asks optimum α above-mentioned _i, β _i, in the process of b, c, d, e, f, adopt least square method or use gauss-newton method to ask optimum α _i, β _i, b, c, d, e, f.

5. a kind of frequency field forward modeling method according to claim 1, is characterized in that, also comprise after reading in Wavelet parameter and rate pattern, loads perfectly matched layer PML boundary condition, for the wave field of absorbing boundary reflection on the border of described rate pattern.

6. frequency field just drills a device, it is characterized in that comprising,

Difference block, for setting up the difference formula of 17 following dot formats,

Wherein

${\stackrel{\&OverBar;}{P}}_{m+1,n}={\mathrm{\&alpha;}}_2{P}_{m+1,n}+\frac{1-{\mathrm{\&alpha;}}_2}{2}(P_{m+1,n+1}+P_{m+1,n-1}),$

${\stackrel{\&OverBar;}{P}}_{m-1,n}={\mathrm{\&alpha;}}_2{P}_{m-1,n}+\frac{1-{\mathrm{\&alpha;}}_2}{2}(P_{m-1,n+1}+P_{m-1,n-1})---\left(3\right)$

${\stackrel{\&OverBar;}{P}}_{m+2,n}={\mathrm{\&alpha;}}_1{P}_{m+2,n}+\frac{1-{\mathrm{\&alpha;}}_1}{2}(P_{m+2,n+2}+P_{m+2,n-2}),$

${\stackrel{\&OverBar;}{P}}_{m-2,n}={\mathrm{\&alpha;}}_1{P}_{m-2,n}+\frac{1-{\mathrm{\&alpha;}}_1}{2}(P_{m-2,n+2}+P_{m-2,n-2}),$

${\stackrel{\&OverBar;}{P}}_{m,n}={\mathrm{\&alpha;}}_3{P}_{m,n}+{\mathrm{\&alpha;}}_4(P_{m,n+1}+P_{m,n-1})+\frac{1-{2\mathrm{\&alpha;}}_4-{\mathrm{\&alpha;}}_3}{2}(P_{m,n+2}+P_{m,n-2}),$

With

${\tilde{P}}_{m,n+1}={\mathrm{\&beta;}}_2{P}_{m,n+1}+\frac{1-{\mathrm{\&beta;}}_2}{2}(P_{m+1,n+1}+P_{m-1,n+1}),$

${\tilde{P}}_{m,n-1}={\mathrm{\&beta;}}_2{P}_{m,n-1}+\frac{1-{\mathrm{\&beta;}}_2}{2}(P_{m+1,n-1}+P_{m-1,n-1}),---\left(4\right)$

${\tilde{P}}_{m,n+2}={\mathrm{\&beta;}}_1{P}_{m,n+2}+\frac{1-{\mathrm{\&beta;}}_1}{2}(P_{m+2,n+2}+P_{m-2,n+2}),$

${\tilde{P}}_{m,n-2}={\mathrm{\&beta;}}_1{P}_{m,n-2}+\frac{1-{\mathrm{\&beta;}}_1}{2}(P_{m+2,n-2}+P_{m-2,n-2}),$

${\tilde{P}}_{m,n}={\mathrm{\&beta;}}_3{P}_{m,n}+{\mathrm{\&beta;}}_4(P_{m+1,n}+P_{m-1,n})+\frac{1-{2\mathrm{\&beta;}}_4-{\mathrm{\&beta;}}_3}{2}(P_{m+2,n}+P_{m-2,n}),$

Wherein P _m,n≈ P (m Δ z, n Δ x), Δ x, Δ z represent x direction sampling interval and z direction sampling interval, m and n represents the mesh coordinate in z direction and the mesh coordinate in x direction, α _i, β _i, b, c, d, e, f represent weighting coefficient when differentiate;

Sparse matrix module, for constructing sparse matrix;

Read in module, for reading in Wavelet parameter and rate pattern;

Single-frequency wave field module, obtains single-frequency wave field for frequency cycle;

Inversefouriertransform module, for Fourier transform of negating to all frequency wave fields;

Result output module, is just drilled result.

7. a kind of frequency field according to claim 6 just drills device, it is characterized in that, also comprises weighting coefficient module, is connected between described difference block and sparse matrix module, for asking for the optimal value of described weighting coefficient.

8. a kind of frequency field according to claim 7 just drills device, it is characterized in that, the optimal value that described weighting coefficient module asks for weighting coefficient comprises, further by wave equation (1) set out, plane wave substitute into formula (1), arrangement can obtain dispersion equation, wherein, and P ₀for the initial amplitude of ripple, i are imaginary unit,

A＝cos(k _xΔx),

B＝cos(k _zΔz), （5）

C＝cos(2k _xΔx),

D＝cos(2k _zΔz),

$E=\frac{1}{6}{\mathrm{\&alpha;}}_1(C-\mathrm{CD})-\frac{8}{3}{\mathrm{\&alpha;}}_2(A-\mathrm{AB})+\frac{5}{2}{\mathrm{\&alpha;}}_3(1-D)+{5\mathrm{\&alpha;}}_4(B-D)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}D,$

$F=\frac{1}{6}{\mathrm{\&beta;}}_1(D-\mathrm{CD})-\frac{8}{3}{\mathrm{\&beta;}}_2(B-\mathrm{AB})+\frac{5}{2}{\mathrm{\&beta;}}_3(1-C)+{5\mathrm{\&beta;}}_4(A-C)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}C,$

Wherein k representative is wave number,

As dx >=dz, definition r=Δ x/ Δ z, therefore above-mentioned formula (4) is rewritten into following form

A＝cos(k _xΔx),

B＝cos(k _zΔz),

C＝cos(2k _xΔx), （6）

D＝cos(2k _zΔz),

$E=\frac{1}{6}{\mathrm{\&alpha;}}_1(C-\mathrm{CD})-\frac{8}{3}{\mathrm{\&alpha;}}_2(A-\mathrm{AB})+\frac{5}{2}{\mathrm{\&alpha;}}_3(1-D)+{5\mathrm{\&alpha;}}_4(B-D)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}D,$

$F=\frac{1}{6}{\mathrm{\&beta;}}_1(D-\mathrm{CD})-\frac{8}{3}{\mathrm{\&beta;}}_2(B-\mathrm{AB})+\frac{5}{2}{\mathrm{\&beta;}}_3(1-C)+{5\mathrm{\&beta;}}_4(A-C)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}C,$

r＝Δx/Δz,

When dx >=dz, net point G required in minimum wavelength is defined as G=λ/Δ x, and wherein λ is wavelength,

So far, normalized phase velocity formula is

$\frac{V_{\mathrm{ph}}}{v}=\frac{G}{2\mathrm{\&pi;}}{\left[\frac{E+r^{2}F}{b+2c(A+B)+4\mathrm{dAB}+2e(C+D)+4\mathrm{fCD}}\right]}^{1/2},$

$A=\cos \left(\frac{2\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{G}\right)$

$B=\cos \left(\frac{2\mathrm{\&pi;}\sin \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right),---\left(7\right)$

$C=\cos \left(\frac{4\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{G}\right),$

$D=\cos \left(\frac{4\mathrm{\&pi;}\sin \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right)$

$E=\frac{1}{6}{\mathrm{\&alpha;}}_1(C-\mathrm{CD})-\frac{8}{3}{\mathrm{\&alpha;}}_2(A-\mathrm{AB})+\frac{5}{2}{\mathrm{\&alpha;}}_3(1-D)+{5\mathrm{\&alpha;}}_4(B-D)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}D,$

$F=\frac{1}{6}{\mathrm{\&beta;}}_1(D-\mathrm{CD})-\frac{8}{3}{\mathrm{\&beta;}}_2(B-\mathrm{AB})+\frac{5}{2}{\mathrm{\&beta;}}_3(1-C)+{5\mathrm{\&beta;}}_4(A-C)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}C,$

Wherein

Try to achieve the α of above-mentioned optimum _i, β _i, b, c, d, e, f, make V _phand difference is minimum between v;

As dx<dz, definition r=Δ z/ Δ x, G is defined as G=λ/Δ z;

Above-mentioned formula (4) is rewritten as

A＝cos(k _xΔx),

B＝cos(k _zΔz),

C＝cos(2k _xΔx),

D＝cos(2k _zΔz),

r＝Δz/Δx

$E=\frac{1}{6}{\mathrm{\&alpha;}}_1(C-\mathrm{CD})-\frac{8}{3}{\mathrm{\&alpha;}}_2(A-\mathrm{AB})+\frac{5}{2}{\mathrm{\&alpha;}}_3(1-D)+{5\mathrm{\&alpha;}}_4(B-D)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}D,$

$F=\frac{1}{6}{\mathrm{\&beta;}}_1(D-\mathrm{CD})-\frac{8}{3}{\mathrm{\&beta;}}_2(B-\mathrm{AB})+\frac{5}{2}{\mathrm{\&beta;}}_3(1-C)+{5\mathrm{\&beta;}}_4(A-C)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}C,$

Then normalized phase velocity formula is

$\frac{V_{\mathrm{ph}}}{v}=\frac{G}{2\mathrm{\&pi;}}{\left[\frac{E+r^{2}F}{b+2c(A+B)+4\mathrm{dAB}+2e(C+D)+4\mathrm{fCD}}\right]}^{1/2},$

$A=\cos \left(\frac{2\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right),$

$C=\cos \left(\frac{4\mathrm{\&pi;}\cos \left(\mathrm{\&theta;}\right)}{\mathrm{rG}}\right),$

$D=\cos \left(\frac{4\mathrm{\&pi;}\sin \left(\mathrm{\&theta;}\right)}{G}\right)$

$E=\frac{1}{6}{\mathrm{\&alpha;}}_1(C-\mathrm{CD})-\frac{8}{3}{\mathrm{\&alpha;}}_2(A-\mathrm{AB})+\frac{5}{2}{\mathrm{\&alpha;}}_3(1-D)+{5\mathrm{\&alpha;}}_4(B-D)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}D,$

$F=\frac{1}{6}{\mathrm{\&beta;}}_1(D-\mathrm{CD})-\frac{8}{3}{\mathrm{\&beta;}}_2(B-\mathrm{AB})+\frac{5}{2}{\mathrm{\&beta;}}_3(1-C)+{5\mathrm{\&beta;}}_4(A-C)-\frac{8}{3}\mathrm{AB}+\frac{1}{6}\mathrm{CD}+\frac{5}{2}C,$

Try to achieve the α of above-mentioned optimum _i, β _i, b, c, d, e, f, make V _phand difference is minimum between v.

9. a kind of frequency field according to claim 8 just drills device, it is characterized in that, described weighting coefficient module asks optimum α above-mentioned _i, β _i, in the process of b, c, d, e, f, adopt least square method or use gauss-newton method to ask optimum α _i, β _i, b, c, d, e, f.

10. a kind of frequency field according to claim 1 just drills device, it is characterized in that, also comprise absorbing boundary module, read between module and single-frequency wave field module described in being connected to, for loading perfectly matched layer PML boundary condition on the border of described rate pattern, for the wave field of absorbing boundary reflection.