CN103675894A - Method for synthesizing seismic records based on three-dimensional Gaussian beam ray tracing and frequency domain - Google Patents

Abstract

The invention relates to a method for synthesizing seismic records based on three-dimensional Gaussian beam ray tracing and a frequency domain. The method comprises the following steps of: 1) performing kinematics tracking by utilizing a common ray tracing way to obtain a ray path and kinematics characteristics; 2) performing dynamics tracking by utilizing a Gaussian ray beam to calculate dynamics characteristics; 3) synthesizing the seismic records in the frequency domain. As a Futyterman equation based attenuation algorithm, a frequency based green function and a weak anisotropic model of a coupling ray theory can be considered for the synthesis of the seismic records in the frequency domain, the dynamics characteristics of a seismic wave field can be described more delicately in comparison with a conventional method, and the seismic records synthesized by utilizing the ray theory is close to the precision of the wave theory. The method for synthesizing the seismic records based on the three-dimensional Gaussian beam ray tracing and the frequency domain can calculate information such as ray amplitude and phase shift conveniently, can be applied to the fields such as manufacturing of synthesized records, prestack depth migration and waveform tomography inversion.

Description

A kind of based on three-dimensional Gaussian beam ray tracing and frequency field theogram method

Technical field

The present invention relates to geophysical prospecting for oil field, particularly about a kind of based on three-dimensional Gaussian beam ray tracing and frequency field theogram method.

Background technology

Seismic wave field numerical simulation, through development for many years, has obtained a large amount of achievements in research at present, mainly comprises two large class, i.e. Wave Equation Numerical method and ray theory forward modeling methods.Wave Equation Numerical method can disclose ripple and change in amplitude, frequency, the phase place of Propagation, the dynamic characteristic that reflects truly efferent echo, it is the effective means of carrying out wave field research, but it is large to take calculator memory, many while expending machine, sometimes the situation of complex geological structure is also difficult to research.Ray theory forward modeling method is according to model structure and velocity variations, movement locus and the ripple process velocity interface of tracking seismic event in underground medium reflects, transmission changes, and change approximate simulation wave field characteristics by calculating energy loss and the phase place of the ripples such as wave field diffusion, transmission coefficient, reflection coefficient, ray theory forward modeling method is mainly used in studying track, the distribution of reflection spot, the kinematics character of the ripples such as distribution range of wave field that ripple is propagated, and is also the basis of carrying out tomographic inversion.Ray tracing can be divided into kinematics ray tracing and dynamics ray tracing two classes at present.Kinematics ray tracing need to calculate raypath, wavefront surface, information while walking, and dynamics ray tracing also needs to calculate the information such as ray amplitude, phase shift on this basis.Kinematics ray tracing is the basis of dynamics ray tracing.Ray-tracing scheme has extremely important effect in geophysics tomography, and it is one of effective ways of seismic wave propagation problem under research medium arbitrary speed distribution situation.The speed that the path accuracy of ray tracing and tracking are calculated is directly determining quality and the computing velocity of imaging.Therefore, research is ray-tracing scheme fast and accurately, has the practical significance of particular importance for tomography problem.

Numerous solving in raypath and the method while walking of the prior art, traditional method is x ray level algebraic method, the Karal of the method Shi You U.S. and Keller (1959) introduce elastic wave field by electromagnetic achievement in research and derive and take the lead in proposing in the process of ray series solution of the elastokinetics equation of motion, and it is the predecessor of Asymptotic ray theory method.Development in recent years the new methods such as eikonal equation method, wavefront reconstruction method, critical path method (CPM), simulated annealing, genetic algorithm.The concrete grammar of ray tracing of the prior art is a lot, conventionally use trial fire method (shooting method), bending method and wave-front method, wherein, trial fire method utilizes Snell law to determine raypath and acceptance point, constantly adjust the incident angle of ray, make between ray outgoing position and geophone station constantly close, till finally meeting given accuracy requirement, thereby the ray tracing of realizing point-to-point transmission calculates, it is a kind of important method that solves two point tracking problems.Trial fire method need to attempt repeatedly initially firing for adjustment, by continuous iterative modifications ray incident angle, obtains raypath accurately, more consuming time, particularly prominent place under three-dimensional situation.Bending method is fixedly source point and acceptance point, Fermat principle based on minimum ray whilst on tour obtains correct raypath by iteration repeatedly, and to general medium, bending method is to solve Nonlinear Optimization Problem in fact, in computation process, be absorbed in sometimes local convergence, obtain the optimum solution of the non-overall situation.When trial fire method and bending method can solving model multipath be walked, the aspect such as tomography, synthetic record while being generally used for away.Wave-front method, based on Huygens (Huygens-Fresnel) principle, seismic event can form a wave front in communication process, wave front is as new focus divergent-ray, through a plurality of secondary sources, finally propagate into geophone station, first medium to be divided into many net points, require ray must pass through these net points, source point and geophone station are also respectively on net point, ray is from the residing net point of source point, while walking with minimum via each net point, arrive geophone station, these net points of ray process have just formed Shortest ray path, when net point quantity increases, calculated amount is also by proportional increase.

Early 1990s, a lot of defects of above-mentioned several method also came out gradually, mainly contain along with pre-stack depth migration is realized industrialization gradually: be 1. difficult to treatment media medium velocity and change violent situation; The difficulty when overall situation minimum in while 2. asking for many-valued walk is walked; 3. counting yield is lower; 4. shadow region and complex structural area ray coverage density are low.Therefore, ray-tracing scheme is mainly carried out work around the problems referred to above in recent ten years.The people (1991) such as Vidale (1988,1990) and Podvin are from eikonal equation, by first obtaining distribute the again way of fast descent direction of field computing time of time field, obtain each acceptance point to the raypath of focus.Subsequently, the people such as Qin (1992) have done improvement to the method for Vidale (1988), have proposed wavefront spread method, although the method has reduced the instability of former method, but make the huge increasing of calculated amount.

Summary of the invention

For the problems referred to above, the object of this invention is to provide a kind ofly based on three-dimensional Gaussian beam ray tracing and frequency field theogram method, can obtain meticulous dynamics wave field information.

For achieving the above object, the present invention takes following technical scheme: a kind of based on three-dimensional Gaussian beam ray tracing and frequency field theogram method, it comprises the following steps: the method for 1) utilizing ordinary beam to follow the trail of is carried out kinematics tracking, obtain raypath and kinematics character, i.e. the propagation law of the base area seismic wave raypath that seismic wave is propagated in actual formation definitely; 2) utilize Gaussian ray bundle to carry out dynamics tracking, calculate dynamic characteristic; 3) at frequency domain, carry out the calculating of theogram, detailed process is: 3.1), for the ripple of a certain type, calculate from focus, by ray parameter

the abundant intensive ray system of regulation, and take out with each ray and measure accordingly A (O from the result of initial value or two spots ray tracing and dynamic ray tracing _s), q (O _s), p (O _s), v (O _s) and determine

wherein, A (O _s) be ray O _sthe complex amplitude factor at place, q (O _s) and p (O _s) be that ray is at O _sthe direction factor at place, controls reflection and transmission direction of wave travel, v (O _s) ray is at O _sthe speed at place,

it is the weight function of Gauss Ray; 3.2) same central ray is determined respectively to the ray center coordinate (s, n) at different geophone station S place, determine phase factor T (S, O simultaneously _s), and get its real part and imaginary part; 3.3) given frequency interval Δ f, at a certain frequency range f ₁-f ₂interior calculating corresponding with time Re[T (S, O _s)] discrete spectrum, and be stored in array; 3.4) 3.2) finish after, what obtain is the contribution of a certain Gaussian ray bundle to all relevant geophone stations, form with discrete spectrum is stored in relevant each road respectively, then turn back to step 3.1) middle contribution of calculating another Gaussian ray bundle Dui Ge road, result is stacked up, until same first-harmonic Gaussian ray bundle is finished; Then calculate and carry out the wave field of another kind of first-harmonic again, until all first-harmonic Gaussian ray bundles are all finished, the result of gained Shi Ge road reflection coefficient sequence in frequency domain, is stored in each road with the form of discrete spectrum; 3.5) after all first-harmonic Gaussian ray bundles are finished, input time signal; 3.6) according to the form of input signal, with constant duration Δ t, in scope, it being sampled sometime, obtain discrete-time signal, carry out obtaining signal discrete frequency spectrum after fast Fourier transform; 3.7) by step 3.6) in signal discrete frequency spectrum respectively with step 3.3) after discrete spectrum in Zhong Ge road reflection coefficient sequence multiplies each other, then carry out anti-Fourier transform IFFT, obtain the theogram of forms of time and space.

Described step 2) comprise the following steps: 2.1) by ray tracing, calculate the whilst on tour along ray; 2.2) determine the polarization vector of ray; 2.3) utilize dynamic ray tracing, calculate dynamic characteristic; 2.4) determine the amplitude fading in complex vector direction; 2.5) determine caustic amount.

Described step 3.5) in, the form of input signal adopts a kind of in Gabor signal, Berlage signal, Muler signal and Ricker signal.

The present invention is owing to taking above technical scheme, and it has the following advantages: the method that 1, first the present invention utilizes ordinary beam to follow the trail of is carried out kinematics tracking, obtains raypath and kinematics character; Then utilize Gaussian ray bundle to carry out dynamics tracking, calculate dynamic characteristic; Finally calculating the calculating of carrying out theogram on the basis of dynamic characteristic in frequency field, because frequency field is carried out theogram and can be considered the decay algorithm based on Futterman equation, can consider the Green function based on frequency, also can consider the to be coupled weak anisotropy model of ray theory, so compare the dynamic characteristic of description seismic wave field that can be meticulousr with conventional method, make to utilize the synthetic seismologic record of ray theory close to the precision of wave theory.2, the present invention utilizes classical paraxial ray tracing system or dynamics ray tracing system, can calculate easily the information such as ray amplitude, phase shift, be effectively applied to synthetic record, protect the fields such as width pre-stack depth migration, waveform tomographic inversion, for the situation that has caustic and shade in model, can calculate the amplitude of X-wave place by its improved form-Gaussian beam ray theory.3, current, the numerical solution of fluctuation problem is divided into two large classes, one class is to take elastic oscillation mechanical equation numerical solution as basic direct solving method, as method of finite difference, finite element method etc., another kind of is to take inhomogeneous medium wave equation to get high-frequency approximation as the ray method of basic approximate solution, both compare, the former simulation wave field more complete, precision is high, be suitable for the dielectric model that yardstick is less than advantage wave field, but often need mainframe computer, and expend and often make insufferable computing time of people, and wave field is complicated, be not easy to identification and explain.Although the rays method that wave equation is got high-frequency approximation is accurate not as Solution of Wave Equations, be suitable for the model that yardstick is greater than advantage wavelength, and solution speed is fast, wave field is easy to identification and explains.The present invention is further developing at conventional ray tracing, research shows that conventional paraxial ray tracking is not flawless, when for complex model, still meet difficulty, such as caustic district, Deng Fei regular domain, shadow region, paraxial ray back tracking method will lose efficacy, secondly, paraxial ray is followed the trail of the two spots ray tracing that usually needs to do trouble, and counting yield is low.The present invention carries out on the basis of dynamic ray tracing, when solving dynamic ray tracing system of equations, in how much expansions that solve central ray, can also ask for the whilst on tour on central ray neighbor point, whilst on tour second derivative, how much expansions and ray amplitude, follow the trail of and the position of other ray that definite central ray is contiguous simultaneously, under central ray coordinate system, elastic oscillation mechanics wave equation is got to high-frequency approximation and turn to parabolic equation, derive dynamic ray tracing equation and transmission equation, parabolic equation is got to time harmonic solution and obtain Gaussian beam, high-frequency seism wave field is launched or resolves into the wave field that Gaussian beam system represents.Gaussian beam is similar to the major advantage that holography ray theory compares: the one, can eliminate the singularity in the inefficacy districts such as caustic district, and the 2nd, can avoid two spots ray tracing, greatly improve counting yield.At Gaussian beam ray tracing, calculate on the basis of dynamic characteristic, in frequency field, consider the decay algorithm based on Futterman equation, Green function based on frequency, realize the weak anisotropy model of considering coupling ray theory, so compare the dynamic characteristic of description seismic wave field that can be meticulousr with conventional method, make to utilize the synthetic seismologic record of ray theory close to the precision of wave theory.The present invention can be widely used in synthetic record, protect the fields such as width pre-stack depth migration, waveform tomographic inversion.

Accompanying drawing explanation

Fig. 1 is recording geometry of the present invention and raypath schematic diagram;

Fig. 2 is Z component theogram of the present invention (P ripple);

Fig. 3 is X2 component of the present invention theogram (P ripple);

Fig. 4 is Z component theogram of the present invention (P+SV ripple);

Fig. 5 is X2 component of the present invention theogram (P+SV ripple);

Fig. 6 is Z component theogram of the present invention (many direct waves of P+SV+);

Fig. 7 is X2 component of the present invention theogram (many direct waves of P+SV+);

Fig. 8 is that three-dimensional raypath of the present invention shows schematic diagram;

Fig. 9 is that the raypath of the partially different wave field types of non-zero of the present invention shows schematic diagram.

Embodiment

Below in conjunction with drawings and Examples, the present invention is described in detail.

As shown in Fig. 1～9, of the present invention based on three-dimensional Gaussian beam ray tracing and frequency field theogram method, comprise the following steps:

1, the method for utilizing ordinary beam to follow the trail of is carried out kinematics tracking, obtains raypath and kinematics character, i.e. the propagation law of the base area seismic wave raypath that seismic wave is propagated in actual formation definitely.

In seismology, there are two class seismic ray tracing problem baseds: a class is a point ray tracing, known rays initial point (source point) and initial exit direction, solve earthquake wave trajectory, is again initial value ray tracing problem; Another kind of is two spots ray tracing, i.e. the position of known rays initial point (source point) and another observation point (acceptance point), and the initial exit direction of ray is unknown, solves the raypath between 2.Ray tracing of the present invention had both needed to study initial value ray tracing, also needed to study two spots ray tracing.

By eikonal equation, can derive ray tracing equation:

$\left\{\begin{array}{c}\frac{\mathrm{dx}}{\mathrm{d\τ}}=v^2{p}_x,\frac{\mathrm{dy}}{\mathrm{d\τ}}{v}^2{p}_y,\frac{\mathrm{dz}}{\mathrm{d\τ}}=v^2{p}_z\\ \frac{d{p}_x}{\mathrm{d\τ}}=-\frac{1}{v}\frac{\∂v}{\∂x},\frac{d{p}_y}{\mathrm{d\τ}}=-\frac{1}{v}\frac{\∂v}{\∂y},\frac{d{p}_z}{\mathrm{d\τ}}=-\frac{1}{v}\frac{\∂v}{\∂z}\end{array}\right.---\left(1.1\right)$

In formula, $p_x=\∂\mathrm{\τ}/\∂x=\cos \mathrm{\α}/v,$ $p_y=\∂\mathrm{\τ}/\∂y=\cos \mathrm{\β}/v,$ $p_z=\∂\mathrm{\τ}/\∂z=\cos \mathrm{\γ}/v,$ Wherein, p _x, p _y, p _zbe ray exit direction, wherein, α is that angle, the β of emergent ray and x axle is the angle of emergent ray and y axle, and γ is the angle of emergent ray and z axle, and v represents the speed of medium.

Starting condition can be expressed as follows: as τ=τ ₀time, have:

$\left\{\begin{array}{c}x=x_0,y=y_0,z=z_0\\ p_x=p_{x0},p_y=p_{y0},p_z=p_{z0}\end{array}\right.---\left(1.2\right)$

In formula, starting condition { x ₀, y ₀, z ₀determine the position of focus, { p _x0, p _y0, p _z0the initial exit direction of given ray.

2, utilize Gaussian ray bundle to carry out dynamics tracking, calculate dynamic characteristic, detailed process:

2.1) by ray tracing, calculate the whilst on tour along ray;

Suppose that three-dimensional velocity model is to be defined in a space body M:

$M:x_{\min }^i\≤x^i\≤x_{\max }^i(i=\mathrm{1,2,3})$

In formula,

with

the boundary coordinate that represents three-dimensional model, by certain funtcional relationship, parameter that can given three dimensions medium, for example P-wave And S speed V _p, V _s, density p and decay factor

Ray tracing of the present invention has the following steps: (a) ray tracing and calculating are along the whilst on tour of ray; (b) determine the polarization vector of ray; (c) dynamic ray tracing (for example calculating the propagation model of ray); (d) determine the amplitude fading in complex vector direction.

Ray tracing is comprised of the calculating of ray, for example, along ray node-by-node algorithm coordinate.The parameter of ray is determined by independent variable:

$\mathrm{\σ}={\mathrm{\σ}}_0+{\∫}_{{\mathrm{\τ}}_0}^{\mathrm{\τ}}{v}^{\mathrm{NEXPS}}\mathrm{d\τ}={\mathrm{\σ}}_0+{\∫}_{s_0}^s{v}^{\mathrm{NEXPS}-1}\mathrm{ds}---\left(2.1\right)$

In formula, τ is whilst on tour, and what σ represented is certain ray, σ ₀what represent is the reference position of certain ray.S is the arc length along ray, and v is the wave propagation velocity of respective propagation type, and NEXPS integral indices can be determined in use as required, NEXPS=0, and ± 1, ± 2 .....We think that ray tracing system is the form of six first-order ordinary differential equation systems:

$\left.\begin{array}{c}\frac{d{x}^i}{\mathrm{d\σ}}=v^{2-\mathrm{NEXPS}}{G}^{\mathrm{ij}}{p}_j,\\ \frac{d{p}_i}{\mathrm{d\σ}}=v^{2-\mathrm{NEXPS}}(-v^{-3}\frac{\∂v}{\∂x^i}+{\mathrm{\Γ}}_{\mathrm{ij}}^k{G}^{\mathrm{jl}}{p}_k{p}_l).\end{array}\right\}---\left(2.2\right)$

In formula, x ⁱthe coordinate of putting along on ray, p _ibe slowness vector in the covariant component of an xi, and

metric tensor G wherein ^ijand Cristoffel symbols

be defined as follows: for right hand curve coordinate system (x ¹, x ², x ³), the local characteristics of coordinate system can be by covariant component G _ij=G _ij(x ^k) and contravariant component G ^ij=G ^ij(x ^k) i wherein, j, k=1,2,3, Cristoffel symbols are described

on interface, ray is propagated according to Snell's law.

Equation:

$\mathrm{\τ}={\mathrm{\τ}}_0+{\∫}_{{\mathrm{\σ}}_0}^{\mathrm{\σ}}{v}^{-\mathrm{NEXPS}}\mathrm{d\σ}---\left(2.3\right)$

The real part of whilst on tour (2.3) can be derived and be obtained by formula (2.1), and the imaginary part of whilst on tour is calculated by following formula:

${\mathrm{\τ}}^{\mathrm{IM}}=\frac{1}{2}{t}^{*}={\∫}_{{\mathrm{\τ}}_0}^{\mathrm{\τ}}{\left(2Q\right)}^{-1}{v}^{-\mathrm{NEXPS}}\mathrm{d\σ}---\left(2.4\right)$

In formula, Q is relevant quality factor.

2.2) determine the polarization vector of ray;

At random select a ray to represent with Ω, the complete ray tracing of definite polarization vector along ray Ω under multiple theory all plays very important effect.If determine the direction of the displacement vector of the ripple of propagating along Ω, especially, in S ripple (shear wave) situation, just must know polarization vector.In addition, polarization vector is as the basic point of the ray center coordinate system vector being connected with Ω, this ray center coordinate system under many circumstances of great use, especially when dynamic ray tracing.

Popov and Psencik (1978a, b) introduce the ray center coordinate system of quadrature in seismology.The present invention is by (q ¹, q ², q ³=s) define this ray center coordinate, wherein, s is the arc length along ray Ω, q ^kthe Cartesian coordinates plane of initial point on ray Ω, the upper point of this plane tangent rays Ω s=q ³wavefront surface.By three unit vector e ₁, e ₂, e ₃(polarization vector) forms the base vector of ray center coordinate system, wherein, and e ₃tangent rays Ω, e ₁and e ₂perpendicular ray Ω, v is the speed of seismic wave propagation, it is in order to make ray center coordinate system quadrature that this mode is introduced.

Equation:

$\frac{d{H}_{\mathrm{ij}}}{\mathrm{d\σ}}=v^{1-\mathrm{NEXPS}}[\frac{\∂v}{\∂x^k}{G}^{\mathrm{kl}}{H}_{\mathrm{lj}}{p}_i-p_k{G}^{\mathrm{kl}}{H}_{\mathrm{lj}}\frac{\∂v}{\∂x^i}v{\mathrm{\Γ}}_{\mathrm{ik}}^m{G}^{\mathrm{kl}}{p}_l{H}_{\mathrm{mj}}]---\left(2.5\right)$

In formula, G is that metric tensor, Γ are Cristoffel symbols, H _ijray center base vector e _jcovariant component, (e _j) _i=H _ij.

Equation (2.5) is reduced to

$\frac{d{H}_{\mathrm{iJ}}}{\mathrm{d\σ}}=v^{1-\mathrm{NEXPS}}[\frac{\∂v}{\∂x^k}{G}^{\mathrm{kl}}{H}_{\mathrm{iJ}}{p}_i+v{\mathrm{\Γ}}_{\mathrm{ik}}^m{G}^{\mathrm{kl}}{p}_l{H}_{\mathrm{mJ}}]---\left(2.6\right)$

Vector e _jperpendicular ray Ω.As unit vector e ₃while being tangential to ray Ω (j=3 is in formula (2.5)),

$H_{i3}=v{p}_i=G_{\mathrm{ij}}\frac{d{x}^j}{\mathrm{ds}}---\left(2.7\right)$

Formula (2.5) is just equivalent to the second equation in ray tracing system (2.2) like this $d{p}_i/\mathrm{d\σ}=v^{2-\mathrm{NEXPS}}({-v}^{-3}\mathrm{dv}/x^i+{\mathrm{\Γ}}_{\mathrm{ij}}^k{G}^{\mathrm{jl}}{p}_k{p}_l),$ It is enough used for unique compute vectors e ₁covariant component H _i1(i=1,2,3), because when carrying out ray tracing with formula (2.2), e ₃depend on formula (2.7), and e ₂always perpendicular to e ₁and e ₃:

H _i2＝ε _ijkG ^jmH _m3G ^knH _n1[det(G ^rs)] ^-1/2 （2.8）

In formula, ε ₁₂₃=ε ₂₃₁=ε ₃₁₂=1, ε ₁₃₂=ε ₃₂₁=ε ₂₁₃=-1, ε _ijk=0 is applicable to other all ijk combinations.On interface, ray center coordinate system is around the Vector Rotation perpendicular to plane of incidence, like this after rotation, and according to theory, vector e ₃the ripple of the tangent rays in producing is overlapped.

2.3) utilize dynamic ray tracing, calculate dynamic characteristic;

Dynamic characteristic in the present invention can refer to the propagation model of ray, and dynamic ray tracing is comprised of four solutions along the linear first-order differential equation group of ray Ω, and component is 4 * 1 matrix W,

$\frac{\mathrm{dW}}{\mathrm{d\σ}}=v^{2-\mathrm{NEXPS}}\mathrm{SW}---\left(2.9\right)$

Wherein

$S=\left[\begin{array}{cc}0& I\\ {-v}^{-2}V& 0\end{array}\right]$

0 and I be 0 and 11 * 1 matrix, V is that 2 * 2 components are:

$V_{\mathrm{IJ}}=(\frac{{\∂}^{2}v}{\∂x^k\∂x^l}-{\mathrm{\Γ}}_{\mathrm{kl}}^m\frac{\∂v}{\∂x_m})G^{\mathrm{kr}}{G}^{\mathrm{ls}}{H}_{\mathrm{rI}}{H}_{\mathrm{sJ}}---\left(2.10\right)$

Matrix.The first two component of W is the ray center coordinate (q of paraxial ray ¹, q ²) (single order Taylor expansion), what two components of other of W were corresponding is the ray center component of the slow vector of paraxial ray.

Dynamic ray tracing formula (2.9) has four linear independent solutions, adopts Π (σ, σ ₀) as the initial point σ=σ at starting condition ray ₀, Π (σ, σ ₀basis 4 * 4 matrixes of the Line independent solution under)=I.Here I is 4 * 4 unit matrix, basis matrix Π (σ, σ ₀) be also ray propagates model, or the propagation model of dynamic ray tracing system.Along directions of rays it meet below relation:

${\mathrm{\Π}}^T\mathrm{\Σ\Π}=\mathrm{\Σ},\mathrm{with\Σ}=\left[\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ -1& 0& 0& 0\\ 0& -1& 0& 0\end{array}\right]---\left(2.11\right)$

This is also the coupling of Π.

Any solution W (σ) of dynamic ray tracing system equations (2.9) can represent by following form:

W(σ)＝Π(σ,σ ₀)W(σ ₀) (2.12)

The solution of ray tracing system equations (2.9) has very important application in seismology.The ray of assumed calculation meets two parameter γ of ray system ¹, γ ².Like this, dynamic ray tracing system can be used for determining 2 * 2 matrix Q along directions of rays ^rand P ^r, matrix unit is:

$\left.\begin{array}{c}{Q}_{\mathrm{IJ}}^R={\left[\frac{\&PartialD;q^I}{\&PartialD;{\mathrm{\&gamma;}}^J}\right]}_{{\mathrm{<figure-callout\; id="1"\; label="q"\; filenames="HDA0000444945350000011.png"\; state="\{\{state\}\}">q</figure-callout>}}^1={\mathrm{<figure-callout\; id="2"\; label="q"\; filenames="HDA0000444945350000023.png"\; state="\{\{state\}\}">q</figure-callout>}}^2=0},\\ P_{\mathrm{IJ}}^R={\left[\frac{{\&PartialD;}^2\mathrm{\&tau;}}{\&PartialD;q^I\&PartialD;{\mathrm{\&gamma;}}^J}\right]}_{{\mathrm{<figure-callout\; id="1"\; label="q"\; filenames="HDA0000444945350000011.png"\; state="\{\{state\}\}">q</figure-callout>}}^1={\mathrm{<figure-callout\; id="2"\; label="q"\; filenames="HDA0000444945350000023.png"\; state="\{\{state\}\}">q</figure-callout>}}^2=0}={\left[\frac{\&PartialD;p_I^{\left(q\right)}}{\&PartialD;{\mathrm{\&gamma;}}^J}\right]}_{{\mathrm{<figure-callout\; id="1"\; label="q"\; filenames="HDA0000444945350000011.png"\; state="\{\{state\}\}">q</figure-callout>}}^1={\mathrm{<figure-callout\; id="2"\; label="q"\; filenames="HDA0000444945350000023.png"\; state="\{\{state\}\}">q</figure-callout>}}^2=0}\end{array}\right\}---\left(1.13\right)$

In formula, q ⁱray parameter γ ¹, γ ²the ray center coordinate of definite adaxial ray;

the ray center component of corresponding paraxial slow vector,

matrix Q ^ralso be called geometrical attenuation model, | detQ ^r| be geometrical attenuation.It represents expansion and the contraction of ray channel.Geometrical attenuation depends on the parameter of ray.

Introduce ray coordinates (γ ¹, γ ², γ ³), γ ¹and γ ²ray parameter, γ ³it is the parameter (as arc length s) along ray.Matrix Q like this ^ralso can be understood as from ray coordinates (γ ¹, γ ²) to the ray center coordinate (q along ray Ω ¹, q ²) transition matrix.Similar, matrix P ^rcan be regarded as from ray coordinates (γ ¹, γ ²) to phase space coordinate

transition matrix, this matrix X of 4 * 2 ^rbe decided to be:

$X^R=\left[\begin{array}{c}{Q}^R\\ P^R\end{array}\right]---\left(2.14\right)$

Meet along ray Ω direction

$\frac{d{X}^R}{\mathrm{d\&sigma;}}=v^{2-\mathrm{NEXPS}}{\mathrm{SX}}^R---\left(2.15\right)$

Dynamic ray tracing system.The solution of formula (2.15) can be expressed as following form:

X ^R(σ)＝Π(σ,σ ₀)X ^R(σ ₀) (2.16)

If Q ^rand P ^rbe known, also can determine whilst on tour field M ^r(σ) 2 * 2 matrixes of second-order partial differential coefficient, matrix element is:

$M_{\mathrm{IJ}}^R={\left[\frac{{\&PartialD;}^2\mathrm{\&tau;}}{\&PartialD;q^I\&PartialD;q^J}\right]}_{{\mathrm{<figure-callout\; id="1"\; label="q"\; filenames="HDA0000444945350000011.png"\; state="\{\{state\}\}">q</figure-callout>}}^1={\mathrm{<figure-callout\; id="2"\; label="q"\; filenames="HDA0000444945350000023.png"\; state="\{\{state\}\}">q</figure-callout>}}^2=0}---\left(2.17\right)$

Q ^r, P ^rand M ^rrelation be

M ^R＝P ^R(Q ^R) ^-1 (2.18)

The conversion solution of the dynamic ray tracing system at nearly interface is determined by the Snell's law of adaxial ray.

2.4) determine the amplitude fading in complex vector direction;

The amplitude of vector decay:

$A_i=A_i^R{\left(\mathrm{\&rho;v}\right|\det Q^R\left|\right)}^{1/2}---\left(2.19\right)$

It is the ray amplitude representing at ray center coordinate system wave impedance (density p and speed ν's is long-pending) and geometrical attenuation are multiplied by the complex vector skew of (comprising the phase deviation causing due to frequency dispersion) | detQ ^r| the square root of product.Suppose the amplitude A of vector decay _istarting point at ray is unit one.If this starting point is positioned at Free Surface or any internal interface, inverting coefficient is included in A _ithe inside.

Definition 3 * 3 matrix A _ijby corresponding to S ripple at starting point e ₁the amplitude A of the vector decay of the polarization in direction _i1, corresponding to S ripple at starting point e ₂the amplitude A of the vector decay of the polarization in direction _i2with the vector decay amplitude A corresponding to starting point P ripple _i3form.The component A of vector decay amplitude _ijat ray center coordinate system, represent.Ray unit in single thing piece inside is constant, except frequency dispersion, because at this moment according to their will be multiplied by-i or-1 of the type of frequency dispersion.On interface, the amplitude of vector decay can convert by reducing reflection/transmission (R/T) coefficient.

$R^E=R{\left[\frac{\widetilde{\mathrm{\&rho;}}\tilde{v}\left|\cos \widetilde{\mathrm{\&alpha;}}\right|}{\mathrm{\&rho;v}\left|\cos \mathrm{\&alpha;}\right|}\right]}^{1/2}---\left(2.20\right)$

In formula, R is plane wave reflection/transmission coefficient.These amounts with conjugate of symbol are corresponding to relevant reflection/transmission ripple, and those do not have the corresponding incident wave of amount of conjugate of symbol; In two kinds of situations of incidence point, α and reflection/transmission ripple incident angle separately, wherein,

what represent is the impedance of medium.

2.5) determine caustic amount;

On the ray calculating, the decay of compute vectors amplitude needs the information of diacaustic point, and this diacaustic point does not have the character of single ray, but has the information of ray field.It depends on the mutual alignment of being calculated ray and polarization ray.If calculated ray at starting point σ=σ ₀2 * 2 model M of second order local derviation of whilst on tour field ^iNIT=M ^r(σ ₀) be known, their mutual alignment can be determined so.This model M ^iNITcan be according to matrix Q ^iNIT=Q ^r(σ ₀) and P ^iNIT=P ^r(σ ₀) be defined as M ^iNIT=P ^iNIT(Q ^iNIT) ^-1.

3, at frequency domain, carry out theogram;

Because theogram is carried out after ray tracing and dynamic ray tracing complete, thereby arbitrfary point O on needed ray in synthetic _sthe amount at place is as τ (O _s), M (O _s) and A (O _s) etc. to every central ray, be all known, by these known amounts, utilize following formula:

Get final product theogram.In formula

by known

determine, N is the ray number that participates in stack, and weight function is:

And

$\mathrm{Re}{\{-i[M\left(O_s\right)-M^R\left(O_s\right)\left]\right\}}^{\frac{1}{2}}>0---\left(3.3\right)$

Q ^r(O _s), M ^r(O _s) be q (O _s) and M (O _s) real part, can be determined by dynamic ray tracing.What in the result due to ray tracing, provide is the component of displacement vector, in 3-D situation, and U (S, ω)={ U _x, U _y, U _z, so synthesizing must be to U _x(S, ω), U _y(S, ω), U _z(S, ω) carries out respectively, so just can obtain the three-component seismologic record of x, y, z, and concrete process is:

3.1) for the ripple of a certain type, calculate from focus, by ray parameter

the abundant intensive ray system of regulation, and take out with each ray and measure accordingly A (O from the result of initial value (or 2 points) ray tracing and dynamic ray tracing _s), q (O _s), p (O _s), v (O _s) and determine

a (O wherein _s) be ray O _sthe complex amplitude factor at place, q (O _s) and p (O _s) be that ray is at O _sthe direction factor at place, controls reflection and transmission direction of wave travel, v (O _s) ray is at O _sthe speed at place, τ (O _s) be that ray is at O _sthe whilst on tour at place,

it is the weight function of Gauss Ray;

3.2) same central ray is determined respectively to the ray center coordinate (s, n) at different geophone station S place, determine phase factor T (S, O simultaneously _s) and get its real part and imaginary part, real part is the time of arrival of ripple, imaginary part is perpendicular to the decay factor of geophone station S place amplitude in central ray direction.

Above-mentioned centre coordinate (s, n) is must know while determining the amplitude of additional components:

n(S,O _s)＝[x(S)-x(O _s)]t _z(O _s)-[z(s)-z(O _s)]t _x(O _s) (3.4)

s(S,O _s)＝[x(S)-x(O _s)]t _x(O _s)-[z(s)-z(O _s)]t _z(O _s) (3.5)

In formula, t _x(O _s), t _z(O _s) be at an O _sx component and the z component of the unit vector t that place and ray are tangent.

Phase factor:

$T(S,O_s)=\mathrm{\&tau;}\left(O_s\right)+v^{-1}\left(O_s\right)I^T\left(O_s\right)X(S,O_s)+\frac{1}{2}{X}^T(S,O_s)W\left(O_s\right)X(S,O_s)---\left(3.6\right)$

Wherein: $X(S,O_s)=\left[\begin{array}{c}x\left(S\right)-x\left(O_s\right)\\ z\left(S\right)-z\left(O_s\right)\end{array}\right]---\left(3.7\right)$

$I\left(O_s\right)=\left[\begin{array}{c}{t}_x\left(O_s\right)\\ t_z\left(O_s\right)\end{array}\right]---\left(3.8\right)$

W(O _s)＝Ω ^T(O _s)A(O _s)Ω(O _s) （3.9）

$\mathrm{\&Omega;}\left(O_s\right)=\left[\begin{array}{cc}{t}_z\left(O_s\right)& -t_x\left(O_s\right)\\ t_x\left(O_s\right)& t_z\left(O_s\right)\end{array}\right]---\left(3.10\right)$

$A\left(O_s\right)=v^{-2}\left(O_s\right)\left[\begin{array}{cc}{v}^2\left(O_s\right)M\left(O_s\right)& -v_n\left(O_s\right)\\ -v_n\left(O_s\right)& -v_s\left(O_s\right)\end{array}\right]---\left(3.11\right)$

$\left\{\begin{array}{c}{v}_n\left(O_s\right)=v_x\left(O_s\right)t_z\left(O_s\right)-v_z\left(O_s\right)t_x\left(O_s\right)\\ v_s\left(O_s\right)=v_x\left(O_s\right)t_x\left(O_s\right)+v_z\left(O_s\right)t_z\left(O_s\right)\end{array}\right.---\left(3.12\right)$

$M\left(O_s\right)=\frac{P\left(O_s\right)}{q\left(O_s\right)}---\left(3.13\right)$

3.3) given frequency interval Δ f, at a certain frequency range f ₁-f ₂interior calculating corresponding with time Re[T (S, O _s)] discrete spectrum, and be stored in array.

3.4) 3.2) finish after, what obtain is the contribution of a certain Gaussian ray bundle to all relevant geophone stations, form with discrete spectrum is stored in relevant each road respectively, then turn back to step 3.1) middle contribution of calculating another Gaussian ray bundle Dui Ge road, result is stacked up, until same first-harmonic Gaussian ray bundle is finished, then calculate again the wave field that carries out another kind of first-harmonic, until all first-harmonic Gaussian ray bundles are all finished, the result of gained Shi Ge road reflection coefficient sequence in frequency domain, is stored in each road with the form of discrete spectrum.

3.5) after all first-harmonic Gaussian ray bundles are finished, input time signal, wherein, except can people for providing discrete signal form, can time signal form for you to choose have following several:

1. Gabor signal

In formula, f _mfor dominant frequency f _m, γ, t ₀,

it is all input parameter.

2. Berlage signal

f(t)＝exp{-β(t-t ₀)}sin[2πf _m(t-t ₀)](t-t ₀) ^α

In formula, f _mfor dominant frequency f _m, β, α, t ₀it is all input parameter.

3. Muler signal

$\left\{\begin{array}{c}f\left(t\right)=\sin \left(\frac{\mathrm{n\&pi;t}}{t_p}\right)-\left(\frac{n}{n+2}\right)\sin \left[(n+2)\mathrm{\&pi;t}\right]t_p,(0\&le;t\&le;t_p)\\ f\left(t\right)=0,(t<0,t>t_p)\end{array}\right.$

In formula, t _pwith n be input parameter.

4. Ricker signal

f(t)＝{1-2[β(t-t ₀)] ²}exp{-[β(t-t ₀)] ²}

In formula, β, t ₀it is input parameter.

3.6) according to the form of input signal, with the constant duration Δ t(sampling interval of seismic trace namely) in scope, it being sampled sometime, obtain discrete-time signal, carry out obtaining signal discrete frequency spectrum after fast Fourier transform.

3.7) by step 3.6) in signal discrete frequency spectrum respectively with step 3.3) after discrete spectrum in Zhong Ge road reflection coefficient sequence multiplies each other, then carry out anti-Fourier transform IFFT, obtain the theogram of forms of time and space.

The various embodiments described above are only for illustrating the present invention, and wherein each step of implementation method etc. all can change to some extent, and every equivalents of carrying out on the basis of technical solution of the present invention and improvement, all should not get rid of outside protection scope of the present invention.

Claims (3)

1. based on three-dimensional Gaussian beam ray tracing and a frequency field theogram method, it comprises the following steps:

1) method of utilizing ordinary beam to follow the trail of is carried out kinematics tracking, obtains raypath and kinematics character, i.e. the propagation law of the base area seismic wave raypath that seismic wave is propagated in actual formation definitely;

2) utilize Gaussian ray bundle to carry out dynamics tracking, calculate dynamic characteristic;

3) at frequency domain, carry out the calculating of theogram, detailed process is:

3.1) for the ripple of a certain type, calculate from focus, by ray parameter

the abundant intensive ray system of regulation, and take out with each ray and measure accordingly A (O from the result of initial value or two spots ray tracing and dynamic ray tracing _s), q (O _s), p (O _s), v (O _s) and determine

wherein, A (O _s) be ray O _sthe complex amplitude factor at place, q (O _s) and p (O _s) be that ray is at O _sthe direction factor at place, controls reflection and transmission direction of wave travel, v (O _s) ray is at O _sthe speed at place,

it is the weight function of Gauss Ray;

3.2) same central ray is determined respectively to the ray center coordinate (s, n) at different geophone station S place, determine phase factor T (S, O simultaneously _s), and get its real part and imaginary part;

3.3) given frequency interval Δ f, at a certain frequency range f ₁-f ₂interior calculating corresponding with time Re[T (S, O _s)] discrete spectrum, and be stored in array;

3.4) 3.2) finish after, what obtain is the contribution of a certain Gaussian ray bundle to all relevant geophone stations, form with discrete spectrum is stored in relevant each road respectively, then turn back to step 3.1) middle contribution of calculating another Gaussian ray bundle Dui Ge road, result is stacked up, until same first-harmonic Gaussian ray bundle is finished; Then calculate and carry out the wave field of another kind of first-harmonic again, until all first-harmonic Gaussian ray bundles are all finished, the result of gained Shi Ge road reflection coefficient sequence in frequency domain, is stored in each road with the form of discrete spectrum;

3.5) after all first-harmonic Gaussian ray bundles are finished, input time signal;

3.6) according to the form of input signal, with constant duration Δ t, in scope, it being sampled sometime, obtain discrete-time signal, carry out obtaining signal discrete frequency spectrum after fast Fourier transform;

3.7) by step 3.6) in signal discrete frequency spectrum respectively with step 3.3) after discrete spectrum in Zhong Ge road reflection coefficient sequence multiplies each other, then carry out anti-Fourier transform IFFT, obtain the theogram of forms of time and space.

2. as claimed in claim 1 a kind of based on three-dimensional Gaussian beam ray tracing and frequency field theogram method, it is characterized in that: described step 2) comprise the following steps:

2.1) by ray tracing, calculate the whilst on tour along ray;

2.2) determine the polarization vector of ray;

2.3) utilize dynamic ray tracing, calculate dynamic characteristic;

2.4) determine the amplitude fading in complex vector direction;

2.5) determine caustic amount.

3. as claimed in claim 1 or 2 a kind of based on three-dimensional Gaussian beam ray tracing and frequency field theogram method, it is characterized in that: described step 3.5), the form of input signal adopts a kind of in Gabor signal, Berlage signal, Muler signal and Ricker signal.

Images