CN103149586A - Tilted layered viscoelasticity dielectric medium wave field forward modelling method - Google Patents

Abstract

The invention discloses a tilted layered viscoelasticity dielectric medium wave field forward modelling method. The method aims at viscoelasticity layered medium, wherein a geologic model of the viscoelasticity layered medium is provided with a tilted interface, a spherical wave stimulated by a point source is decomposed into a plane wave, a character that vector wave equation does not rely on a coordinate system is utilized, under a fixed coordinate system and a moving coordinate system which changes along with the tilted interface, spreading of waves is alternately discussed, firstly, reflection and transmission coefficients of the plane wave of a single tilted interface are deduced through traditional reflection and transmission coefficients, reflection and transmission waves of the interface are obtained, then the reflection and transmission coefficients of the plane wave when passing through multiple tilted interfaces are obtained by recurrence, and finally plane waves after transmitted are combined to obtain a point source excitation wave field. In the process of combination, a rapid oscillatory integral calculation method and an integral path which is suitable for viscoelasticity medium are adopted.

Description

Wave field the Forward Modeling in a kind of inclination stratiform viscoelastic medium

Technical field

The invention belongs to the seismic exploration technique field, relate to a kind of analogy method, especially wave field the Forward Modeling in a kind of inclination stratiform viscoelastic medium, it is the efficient accurate forward simulation technology of seismic wave field in a kind of layered viscoelastic medium for containing tiltedinterface.

Background technology

Wave-field simulation is occupied critical role in oil-gas exploration, although there has been polytype wave-field simulation method, yet for the requirement of the problem that will solve, development wave-field simulation theory and method are very active research fields always.In addition, in fields such as oil-gas explorations, for explanation of seismic figure, or carry out data inversion, usually need only to simulate the ripple of specified type, for example only simulate primary reflection, or analog-converted ripple etc. only.Wave-field simulation method commonly used has polytype, as the method based on ray theory, the numerical solution method, generalized reflection rate method with and promote and obtain anti-/ transmission matrix method etc.These methods respectively have its relative merits.

Ray method can be applicable to more complicated dielectric model, can simulate according in advance given ray code the ripple of interested type, and computing velocity is fast; Weak point is, it is the high-frequency approximation of wave equation, and the seismic response that is used for the Simulation of Complex reservoir has limitation.Adopting the numerical methods of solving wave equation, is the effective way of wave-field simulation in complex dielectrics.Numerical method commonly used has: method of finite difference, and finite element method, pseudo-spectrometry and spectral element method etc., but this class methods simulated earthquake ripple all-wave field, precision is high.But calculated amount is huge, can not look like the ray class methods, selectively only simulates the ripple of the interesting type of our perception.

Layered medium medium wave propagation theory and wave-field simulation technology are subject to people's attention always.The reflectivity method of the people such as Aki and M ü ller development, the all-wave field of dummy level layered viscoelastic medium or the ripple of specified type quickly and accurately.Berkhout has proposed the wave field extrapolation method and can be used for elastic medium wave-field simulation with tiltedinterface.The global generalized of the people such as Chen development is anti-/ transmission matrix method and boundary element method, can simulate with all-wave field in the layered elastic medium of irregular interface.But these methods are not suitable for the layered viscoelastic medium that contains tiltedinterface.

Summary of the invention

For the problems referred to above, the object of the invention is to for the layered viscoelastic medium that contains the inclined wave impedance interface, wave field the Forward Modeling in a kind of inclination stratiform viscoelastic medium is proposed, the method can selectivity be simulated interested various types of ripples, as primary event P ripple, S ripple and transformed wave, and the wave field of simulation is very accurate.

The objective of the invention is to solve by the following technical programs:

Wave field the Forward Modeling in this kind inclination stratiform viscoelastic medium comprises the following steps:

1) spherical wave that will put source forcing is decomposed into plane wave;

2) plane wave is propagated in inclination stratiform viscoelastic medium, is handled as follows when running into the interface: take this interface as

Axle, the intersection point of interface and global coordinate system Z axis are that initial point is set up local coordinate system, Downward perpendicular to the interface, the incident wave in global coordinate system is transformed in local coordinate system; Try to achieve reflection wave and transmitted wave under local coordinate system; Again the reflection wave under local coordinate system and transmitted wave are transformed under global coordinate system, obtain reflection coefficient and transmission coefficient under global coordinate system, and reflection wave and transmitted wave; Choose as required the type of ripple on each interface, the ripple that can choose during the incident of P ripple comprises reflected P P ripple, reflected P S ripple, transmission PP ripple, transmission PS ripple; The ripple that can choose during the incident of S ripple comprises reflection SP ripple, reflection SS ripple, transmission SP ripple, transmission SS ripple; Thereby simulate selected ripple type, the plane wave that obtains simulating;

3) by step 2) wave field that excites according to the following formula resultant source of the plane wave that obtains of simulation:

$P=\frac{1}{2\mathrm{\π}}{\∫}_{\mathrm{\Γ}}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}[R_i\left(p\right)\·\frac{\exp \left[\mathrm{j\ω}(p_{1i}x-q_{1i}z)\right]}{-2\mathrm{jq}}]\mathrm{dp}$

In formula, Γ is path of integration; N is the number of plane wave; R _i(p) be that transmission and the reflection coefficient of plane wave on a plurality of interfaces is comprehensive; J is imaginary unit, and ω is angular frequency; p _1i, q _1iBe respectively i boundary reflection and return the horizontal and vertical slowness of the plane wave of the 1st layer; P, q are respectively the horizontal and vertical slowness of point source place plane wave, and x and z represent position coordinates.

Further, above-mentioned steps 2) in, the primary event P ripple when the simulation focus is the P ripple is established focus and wave detector all at ground floor, and the recursion step that focus place incident wave is successively propagated downwards is:

1. often run into an interface i, namely take this interface as Axle, the intersection point of interface and global coordinate system Z axis are that initial point is set up local coordinate system

Downward perpendicular to the interface; With the incident wave φ in global coordinate system _iBe transformed in local coordinate system and get

2. try to achieve reflection coefficient under local coordinate system And transmission coefficient

And reflection wave

And transmitted wave

3. again the reflection wave under local coordinate system and transmitted wave are transformed under global coordinate system, obtain the reflection coefficient under global coordinate system And transmission coefficient

And reflection wave φ ' _{I, i}With transmitted wave φ _i+1

4. reflection wave φ ' _{I, i}Successively upwards recursion to ground floor, is carried out 1.-3. step, but is only got transmitted wave when running into the interface;

5. transmitted wave φ _i+1Continue downward recursion, get back to step when running into next interface 1..

Further, in above-mentioned steps 2) in, the ripple of simulating other type adopts same procedure and the step of the primary event P ripple when being the P ripple with the simulation focus to carry out.

In above-mentioned steps 3) in, in formula, integration method adopts the GFM integral method that Frazer and Gettrust propose.

The present invention has following some beneficial effect:

(1) the present invention is applicable to the layered viscoelastic medium with tiltedinterface, and the vergence direction at each interface and angle of inclination can be different, only need satisfy the disjoint condition in interface.

(2) the present invention can correctly simulate the decay characteristics of viscoelastic medium medium wave, with respect to method of finite difference, has higher counting yield, is more suitable for layer thicker and violent situation of interlayer parameter variation.

(3) the present invention is very flexible, can simulate respectively various ripples, as primary event P ripple, SV ripple, SH ripple and PS, SP transformed wave, also can simulate our interested multiple reflection, the ripple type that only needs hope is simulated is with level is given with coding form accordingly, and the incident wave incident angle Mintrop wave can occur when surpassing critical angle.

(4) wave field of the present invention's simulation very near analytic solution, can be the numerical method solution contrast reference is provided.So the present invention has significant application value and application prospect in seismic data interpretation and inverting.

Description of drawings

Fig. 1 is the inclination stratiform viscoelastic medium model that the present invention is directed to;

Fig. 2 is reflection and the transmission schematic diagram of ripple when inciding on the interface;

Path of integration schematic diagram when Fig. 3 is midplane ripple resultant source excitation wave field of the present invention.

Embodiment

The geologic model that the present invention is directed to is the Layered Viscoelastic medium that contains tiltedinterface, the spherical wave of point dynamite source is decomposed into plane wave, utilize vector wave equation not rely on this character of coordinate system, under fixed coordinate system and the moving coordinate system with the tiltedinterface variation, wave propagation is discussed alternately.At first derive plane wave reflection and transmission coefficient on single tiltedinterface by traditional reflection and transmission coefficient, obtain reflection and transmitted wave on this interface; Reflection and transmission coefficient when then recursion is tried to achieve plane wave through a plurality of tiltedinterface.Merged by each plane wave after propagating at last and obtain the point dynamite source wave field, adopted the path of integration of oscillatory integral computing method and a suitable viscoelastic medium fast during merging.

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

At first main symbol definition of the present invention is as shown in table 1:

Main notation convention in table 1 the present invention

Fig. 1 is inclination stratiform viscoelastic medium model schematic diagram, and the density of each layer medium is ρ _i(i=1,2 ..., N-1), quality factor is Q _i, p wave interval velocity is α _i, the S wave velocity is β _i, the inclination angle of each tiltedinterface is θ _i, the Z axis intersection point of each tiltedinterface and vertical direction is z _iPoint source is at the first sublevel, and wave detector can be at an arbitrary position.Take horizontal direction as X-axis, vertical downward direction is Z axis, sets up global coordinate system XOZ.

1) spherical wave that will put source forcing is decomposed into plane wave

According to the Weyl integration, the spherical wave of uniform dielectric Point Source excitation is expressed as the stack of a series of plane waves, and under two-dimensional case, the expression formula of each plane wave is:

φ=exp[jω(px+qz)]

Omitted time factor exp (j ω t) in following formula, j is imaginary unit, and ω is angular frequency, and p is horizontal slowness, and q is vertical slowness, and x and z represent position coordinates.

2) plane wave is propagated in inclination stratiform viscoelastic medium, is handled as follows when running into the interface: take this interface as

Axle, the intersection point of interface and global coordinate system Z axis are that initial point is set up local coordinate system,

Downward perpendicular to the interface, the incident wave in global coordinate system is transformed in local coordinate system; Try to achieve reflection wave and transmitted wave under local coordinate system; Again the reflection wave under local coordinate system and transmitted wave are transformed under global coordinate system, obtain reflection coefficient and transmission coefficient under global coordinate system, and reflection wave and transmitted wave.Can choose as required the type of ripple on each interface, the ripple that can choose during the incident of P ripple comprises reflected P P ripple, reflected P S ripple, transmission PP ripple, transmission PS ripple; The ripple that can choose during the incident of S ripple comprises reflection SP ripple, reflection SS ripple, transmission SP ripple, transmission SS ripple; Thereby only simulate our interested ripple type.For example need simulate the primary event P ripple of focus when being the P ripple, establish focus and wave detector all at ground floor, the recursion step that focus place incident wave is successively propagated downwards is:

1. often run into an interface i, namely take this interface as

Axle, the intersection point of interface and global coordinate system Z axis are that initial point is set up local coordinate system Downward perpendicular to the interface.With the incident wave φ in global coordinate system _iBe transformed in local coordinate system and get

2. try to achieve reflection coefficient under local coordinate system

And transmission coefficient

And reflection wave

And transmitted wave

3. again the reflection wave under local coordinate system and transmitted wave are transformed under global coordinate system, obtain the reflection coefficient under global coordinate system

And transmission coefficient

And reflection wave φ ' _{I, i}With transmitted wave φ _i+1

4. reflection wave φ ' _{I, i}Successively upwards recursion to ground floor, is carried out 1.-3. step, but is only got transmitted wave (upward traveling wave) when running into the interface;

5. transmitted wave φ _i+1Continue downward recursion, get back to step when running into next interface 1..

3) wave field that is excited by the plane wave resultant source;

$P=\frac{1}{2\mathrm{\π}}{\∫}_{\mathrm{\Γ}}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}[R_i\left(p\right)\·\frac{\exp \left[\mathrm{j\ω}(p_{1i}x-q_{1i}z)\right]}{-2\mathrm{jq}}]\mathrm{dp}$

Wherein, R _i(p) be that transmission and the reflection coefficient of plane wave on a plurality of interfaces is comprehensive, as when simulating a P-P reflection wave $R_i\left(p\right)=\underset{n=1}{\overset{i-1}{\mathrm{\Π}}}{T}_{\mathrm{pp},n}^u\·R_{\mathrm{pp},i}^d\·\underset{m=1}{\overset{i-1}{\mathrm{\Π}}}{T}_{\mathrm{pp},m}^d;$ P and q are the horizontal and vertical slowness of focus place incident plane wave; p _1iAnd q _1iHorizontal and vertical slowness when arriving wave detector place layer after i boundary reflection for plane wave, they are functions of incident wave horizontal slowness p.Γ is path of integration, and the present invention selects path of integration as shown in Figure 2,1/ α _iWith 1/ β _iBe the slowness of P ripple and S ripple, the real part of the flex point B of path in first quartile should be a bit larger tham max (1/ β _i), to single receiver, θ=tan ^-1(xz _r-z _s|) time Integral Convergence the fastest (x, z _r, z _sRespectively the z axial coordinate in offset distance and receiver and source), a plurality of receivers receive simultaneously in this way, desirable certain compromise value.The inclination angle of line segment OB should make the below of minimum limit in the very close first quartile of OB section.Integration method adopts Frazer and Gettrust to propose GFM (Generalization ofFilon ' s Method) integral method.

Step 2) reflection wave that in, incident wave produces on the interface and the circular of transmitted wave are:

As Fig. 2, being provided with a horizontal slowness is p ₁P ripple φ ₁Propagate downwards among medium I from initial point, amplitude is 1, and this plane wave can be written as:

φ ₁=exp[jω(p ₁x+q ₁z)],

Incident P ripple is transformed into local coordinate system by global coordinate system, obtains local coordinate system Under incident P ripple be:

${\widehat{\mathrm{\φ}}}_1=A_1\exp \left[\mathrm{j\ω}({\hat{p}}_1\hat{x}+{\hat{q}}_1\hat{z})\right],$

Wherein:

A ₁=exp[jω(p ₁x ₀+q ₁z ₀)], $\left[\begin{array}{c}{\hat{p}}_1\\ {\hat{q}}_1\end{array}\right]=C\left[\begin{array}{c}{p}_1\\ \mathrm{<figure-callout\; id="1"\; label="q"\; filenames="HDA00002816807000011.png"\; state="\{\{state\}\}">q</figure-callout>}1\end{array}\right],$ $C=\left[\begin{array}{cc}\cos \mathrm{\&theta;}& \sin \mathrm{\&theta;}\\ -\sin \mathrm{\&theta;}& \cos \mathrm{\&theta;}\end{array}\right].$

At local coordinate system

Under, the P ripple

Incide the interface downwards from medium I, can produce at the interface the up P ripple of reflected back into medium I

With the SV ripple

And be transmitted to downgoing P-wave in medium II

With the SV ripple

The expression formula of each reflection wave and transmitted wave is respectively:

${\widehat{\mathrm{\&phi;}}}_1^{\&prime;}={\hat{R}}_{\mathrm{pp}}^d\&CenterDot;A_1\exp \left[\mathrm{j\&omega;}({\hat{p}}_1\hat{x}-{\hat{q}}_{1p}\hat{z})\right],$

${\widehat{\mathrm{\&psi;}}}_1^{\&prime;}={\hat{R}}_{\mathrm{ps}}^d\&CenterDot;A_1\exp \left[\mathrm{j\&omega;}({\hat{p}}_1\hat{x}-{\hat{q}}_{1s}\hat{z})\right],$

${\widehat{\mathrm{\&phi;}}}_2={\hat{T}}_{\mathrm{pp}}^d\&CenterDot;A_1\exp \left[\mathrm{j\&omega;}({\hat{p}}_1\hat{x}+{\hat{q}}_{2p}\hat{z})\right],$

${\widehat{\mathrm{\&psi;}}}_2={\hat{T}}_{\mathrm{ps}}^d\&CenterDot;A_1\exp \left[\mathrm{j\&omega;}({\hat{p}}_1\hat{x}+{\hat{q}}_{2s}\hat{z})\right],$

Wherein

Reflection and transmission coefficient.Under local coordinate system, each ripple horizontal slowness equates, vertical slowness is provided by following formula:

${\hat{q}}_{\mathrm{ip}}=\sqrt{{\mathrm{\&alpha;}}_i^{-2}-{\hat{p}}_1^2}.$ ${\hat{q}}_{\mathrm{is}}=\sqrt{{\mathrm{\&beta;}}_i^{-2}-{\hat{p}}_1^2}.(i=\mathrm{1,2})$

With respect to local coordinate system, the interface is level, the reflection of this situation and transmission coefficient Computing method can be with reference to pertinent literature (Muller, 1985).All reflection waves and transmitted wave are transformed into global coordinate system by local coordinate system, arrange after conversion the expression formula of various ripples and obtain that under global coordinate system, each reflection wave and transmitted wave are:

${\mathrm{\&phi;}}_1^{\&prime;}=R_{\mathrm{pp}}^d\&CenterDot;\exp \left[\mathrm{j\&omega;}(p_{1p}x-q_{1p}z)\right],$

${\mathrm{\&psi;}}_1^{\&prime;}=R_{\mathrm{ps}}^d\&CenterDot;\exp \left[\mathrm{j\&omega;}(p_{1s}x-q_{1s}z)\right],$

${\mathrm{\&phi;}}_2=T_{\mathrm{pp}}^d\&CenterDot;\exp \left[\mathrm{j\&omega;}(p_{2p}x+q_{2p}z)\right],$

${\mathrm{\&psi;}}_2=T_{\mathrm{ps}}^d\&CenterDot;\exp \left[\mathrm{j\&omega;}(p_{2s}x+q_{2s}z)\right],$

In formula:

$R_{\mathrm{pp}}^d={\hat{R}}_{\mathrm{pp}}^d\&CenterDot;A_1\exp [-\mathrm{j\&omega;}(p_{1p}{x}_0-q_{1p}{z}_0)],$

$R_{\mathrm{ps}}^d={\hat{R}}_{\mathrm{ps}}^d\&CenterDot;A_1\exp [-\mathrm{j\&omega;}(p_{1s}{x}_0-q_{1s}{z}_0)],$

$T_{\mathrm{pp}}^d={\hat{T}}_{\mathrm{pp}}^d\&CenterDot;A_1\exp [-\mathrm{j\&omega;}(p_{2p}{x}_0+q_{2p}{z}_0)],$

$T_{\mathrm{ps}}^d={\hat{T}}_{\mathrm{ps}}^d\&CenterDot;A_1\exp [-\mathrm{j\&omega;}(p_{2s}{x}_0+q_{2s}{z}_0)],$

$\left[\begin{array}{c}{p}_{1p}\\ q_{1p}\end{array}\right]=C\left[\begin{array}{c}{\hat{p}}_1\\ {\hat{q}}_{1p}\end{array}\right].$ $\left[\begin{array}{c}{p}_{1s}\\ q_{1s}\end{array}\right]=C\left[\begin{array}{c}{\hat{p}}_1\\ {\hat{q}}_{1s}\end{array}\right].$ $\left[\begin{array}{c}{p}_{2p}\\ q_{2p}\end{array}\right]=C^{-1}\left[\begin{array}{c}{\hat{p}}_2\\ {\hat{q}}_{2p}\end{array}\right].$ $\left[\begin{array}{c}{p}_{2s}\\ q_{2s}\end{array}\right]=C^{-1}\left[\begin{array}{c}{\hat{p}}_2\\ {\hat{q}}_{2s}\end{array}\right].$

Wherein

Namely be respectively P-P reflection coefficient, P-S reflection coefficient, P-P transmission coefficient and the P-S transmission coefficient of viscoelasticity tiltedinterface under global coordinate system XOZ.

Upwards incided the situation at interface by the below for ripple, establishing has a P ripple φ ' under global coordinate system ₂Upwards propagate in medium II:

φ′ ₂=exp[jω(p ₂x-q ₂z)].

At the downgoing P-wave φ that can produce at the interface reflected back into medium II ₂With SV ripple Ψ ₂, and be transmitted to up P ripple φ ' in medium I ₁With SV ripple Ψ ' ₁Treatment scheme according to the P ripple during from the incident of interface top: first with incident wave φ ' ₂Be transformed under local coordinate system, try to achieve traditional reflective and the transmission coefficient of each reflection wave and transmitted wave in local coordinate system; Again each reflection wave and transmitted wave under local coordinate system are transformed into global coordinate system, obtain that under global coordinate system, the expression formula of reflection wave and transmitted wave is respectively:

${\mathrm{\&phi;}}_2=R_{\mathrm{pp}}^u\&CenterDot;\exp \left[\mathrm{j\&omega;}(p_{2p}x+q_{2p}z)\right],$

${\mathrm{\&psi;}}_2=R_{\mathrm{ps}}^u\&CenterDot;\exp \left[\mathrm{j\&omega;}(p_{2s}x+q_{2s}z)\right],$

${\mathrm{\&phi;}}_1^{\&prime;}=T_{\mathrm{pp}}^u\&CenterDot;\exp \left[\mathrm{j\&omega;}(p_{1p}x-q_{1p}z)\right],$

${\mathrm{\&psi;}}_1^{\&prime;}=T_{\mathrm{ps}}^u\&CenterDot;\exp \left[\mathrm{j\&omega;}(p_{1s}x-q_{1s}z)\right],$

Wherein:

$R_{\mathrm{pp}}^u={\hat{R}}_{\mathrm{pp}}^u\&CenterDot;A^{\&prime;}\exp [-\mathrm{j\&omega;}(p_{2p}{x}_0+q_{2p}{z}_0)],$

$R_{\mathrm{ps}}^u={\hat{R}}_{\mathrm{ps}}^u\&CenterDot;A^{\&prime;}\exp [-\mathrm{j\&omega;}(p_{2s}{x}_0+q_{2s}{z}_0)],$

$T_{\mathrm{pp}}^u={\hat{T}}_{\mathrm{pp}}^u\&CenterDot;A^{\&prime;}\exp [-\mathrm{j\&omega;}(p_{1p}{x}_0-q_{1p}{z}_0)],$

$T_{\mathrm{ps}}^u={\hat{T}}_{\mathrm{ps}}^u\&CenterDot;A^{\&prime;}\exp [-\mathrm{j\&omega;}(p_{1s}{x}_0-q_{1s}{z}_0)],$

A′=exp[jω(p ₂x ₀-q ₂z ₀)]

$\left[\begin{array}{c}{p}_{2p}\\ q_{2p}\end{array}\right]=C^{-1}\left[\begin{array}{c}{\hat{p}}_2\\ {\hat{q}}_{2p}\end{array}\right],$ $\left[\begin{array}{c}{p}_{2s}\\ q_{2s}\end{array}\right]=C^{-1}\left[\begin{array}{c}{\hat{p}}_2\\ {\hat{q}}_{2s}\end{array}\right],$ $\left[\begin{array}{c}{p}_{1p}\\ q_{1p}\end{array}\right]=C\left[\begin{array}{c}{\hat{p}}_1\\ {\hat{q}}_{1p}\end{array}\right],$ $\left[\begin{array}{c}{p}_{1s}\\ q_{1s}\end{array}\right]=C\left[\begin{array}{c}{\hat{p}}_1\\ {\hat{q}}_{1s}\end{array}\right],$

${\hat{q}}_{\mathrm{ip}}=\sqrt{{\mathrm{\&alpha;}}_i^{-2}-{\hat{p}}_1^2},$ ${\hat{q}}_{\mathrm{is}}=\sqrt{{\mathrm{\&beta;}}_i^{-2}-{\hat{p}}_1^2},$ (i=1，2) ${\hat{p}}_1={\hat{p}}_2,$ $\left[\begin{array}{c}{\hat{p}}_2\\ {\hat{q}}_2\end{array}\right]=C^{-1}\left[\begin{array}{c}{p}_2\\ {\mathrm{<figure-callout\; id="2"\; label="q"\; filenames="HDA00002816807000011.png"\; state="\{\{state\}\}">q</figure-callout>}}_2\end{array}\right].$

Wherein

Namely be respectively global coordinate system have a down dip P-P reflection coefficient, P-S reflection coefficient, P-P transmission coefficient and the P-S transmission coefficient at interface.

In like manner, for the situation of SV ripple from tiltedinterface top (or below) incident, it is identical when finding the solution with the incident of P ripple.

Above-mentioned steps 3) in, the GFM integration method is:

In order to utilize the GFM integral method, first integration type in step 3) is rewritten into following form (only provide here the summation number in one):

P=∫ _Γf(p)exp[sg(p)]dp

Wherein:

$f\left(p\right)=\frac{{-R}_i\left(p\right)}{2\mathrm{\&pi;}\&CenterDot;2\mathrm{jq}},$ s=jω·max(x,z), $g\left(p\right)=\frac{p_{1i}x-q_{1i}z}{\max (x,z)}.$

The GFM trapezoidal rule is:

$\underset{a}{\overset{b}{\&Integral;}}f\left(p\right)\exp \left[\mathrm{sg}\left(p\right)\right]\mathrm{dp}=\left\{\begin{array}{cc}\frac{\mathrm{\&delta;p}}{\mathrm{s\&delta;g}}[\mathrm{\&delta;}\left({\mathrm{fe}}^{\mathrm{sg}}\right)-\frac{\mathrm{\&delta;}\left(f\right)\mathrm{\&delta;}\left(e^{\mathrm{sg}}\right)}{\mathrm{s\&delta;}\left(g\right)}]& \mathrm{\&delta;g}\&NotEqual;0\\ \frac{\mathrm{\&delta;p}}{2}[f\left(a\right)e^{\mathrm{sg}\left(a\right)}+f\left(b\right)e^{\mathrm{sg}\left(b\right)}]& \mathrm{\&delta;g}=0\end{array}\right..$

Here b-a equals integration step, to arbitrary function h, and δ h=h (b)-h (a).

Claims (4)

1. wave field the Forward Modeling in an inclination stratiform viscoelastic medium, is characterized in that, comprises the following steps:

1) spherical wave that will put source forcing is decomposed into plane wave;

2) plane wave is propagated in inclination stratiform viscoelastic medium, is handled as follows when running into the interface: take this interface as

Axle, the intersection point of interface and global coordinate system Z axis are that initial point is set up local coordinate system,

Downward perpendicular to the interface, the incident wave in global coordinate system is transformed in local coordinate system; Try to achieve reflection wave and transmitted wave under local coordinate system; Again the reflection wave under local coordinate system and transmitted wave are transformed under global coordinate system, obtain reflection coefficient and transmission coefficient under global coordinate system, and reflection wave and transmitted wave; Choose as required the type of ripple on each interface, the ripple that can choose during the incident of P ripple comprises reflected P P ripple, reflected P S ripple, transmission PP ripple, transmission PS ripple; The ripple that can choose during the incident of S ripple comprises reflection SP ripple, reflection SS ripple, transmission SP ripple, transmission SS ripple; Thereby simulate selected ripple type, the plane wave that obtains simulating;

3) by step 2) wave field that excites according to the following formula resultant source of the plane wave that obtains of simulation:

$P=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{\mathrm{\&Gamma;}}\underset{i=1}{\overset{N-1}{\mathrm{\&Sigma;}}}[R_i\left(p\right)\&CenterDot;\frac{\exp \left[\mathrm{j\&omega;}(p_{1i}x-q_{1i}z)\right]}{-2\mathrm{jq}}]\mathrm{dp}$

In formula, Γ is path of integration; N is the number of plane wave; R _i(p) be that transmission and the reflection coefficient of plane wave on a plurality of interfaces is comprehensive; J is imaginary unit, and ω is angular frequency; p _1i, q _1iBe respectively i boundary reflection and return the horizontal and vertical slowness of the plane wave of the 1st layer; P, q are respectively the horizontal and vertical slowness of point source place plane wave, and x and z represent position coordinates.

2. wave field the Forward Modeling in inclination stratiform viscoelastic medium according to claim 1, it is characterized in that step 2) in, the primary event P ripple when the simulation focus is the P ripple, if focus and wave detector be all at ground floor, the recursion step that focus place incident wave is successively propagated downwards is:

1. often run into an interface i, namely take this interface as

Axle, the intersection point of interface and global coordinate system Z axis are that initial point is set up local coordinate system Downward perpendicular to the interface; With the incident wave φ in global coordinate system _iBe transformed in local coordinate system and get

2. try to achieve reflection coefficient under local coordinate system

And transmission coefficient

And reflection wave

And transmitted wave

3. again the reflection wave under local coordinate system and transmitted wave are transformed under global coordinate system, obtain the reflection coefficient under global coordinate system

And transmission coefficient And reflection wave φ ' _{I, i}With transmitted wave φ _i+1

4. reflection wave φ ' _{I, i}Successively upwards recursion to ground floor, is carried out 1.-3. step, but is only got transmitted wave when running into the interface;

5. transmitted wave φ _i+1Continue downward recursion, get back to step when running into next interface 1..

3. wave field the Forward Modeling in inclination stratiform viscoelastic medium according to claim 1, is characterized in that step 2) in, the ripple of simulating other type adopts same procedure and the step of the primary event P ripple when being the P ripple with the simulation focus to carry out.

4. wave field the Forward Modeling in inclination stratiform viscoelastic medium according to claim 1, is characterized in that, in step 3), in formula, integration method adopts the GFM integral method that Frazer and Gettrust propose.

Images