CN104950332A - Method for calculating plane wave reflection coefficients in elastic multi-layered medium - Google Patents

Abstract

The invention discloses a method for calculating plane wave reflection coefficients in an elastic multi-layered medium. The method comprises the following steps: (1) if plane harmonic waves enter into a target layer series from a medium n+1, generating reflected longitudinal waves and reflected transverse waves in the n+1 medium, generating transmitted longitudinal waves and transmitted transverse waves in a medium 1, and writing each medium layer in a scalar potential form and a vector potential form; (2) determining relation among the displacement, stress and bit shift; (3) loading the scalar potential and vector potential of the step (1) into the step (2) to obtain a displacement and stress relation expression between the (n+1)th layer and the first layer; and (4) obtaining a displacement and stress transfer matrix containing elasticity coefficients of each layer according to the displacement and stress relation expression, which is obtained in the steps (3), between the (n+1)th layer and the first layer, solving, and obtaining reflection and transmission coefficients. According to the method, the influence of multiple waves and transformed waves of the elastic layer series as well as thickness and frequency on the reflection coefficient is considered fully, so that the method is applicable to thin-layer AVO (Amplitude Variation with Offset) analytical simulation and is high in calculation efficiency.

Description

A kind of computing method of Elastic Multi-layer medium midplane wave reflection coefficient

Technical field

The present invention relates to a kind of computing method of Elastic Multi-layer medium midplane wave reflection coefficient, belong to thin layer AVO analogue technique field.

Background technology

Simulation thin layer AVO is divided into two large classes, (1) time domain convolution model, to in a certain special angle situation, to time domain sampled point computational reflect coefficient one by one, wavelet and reflection coefficient is utilized to carry out convolution, form the corresponding seismic trace of this angle, although the method calculates quick, formula is simple, easily realizes, but be difficult to embody thin-bed effect, simulation precision is lower.(2) adopt the Integral Solution method of finite difference simulation thin layer wave field of equations for elastic waves, can simulate thin-bed effect, but counting yield is lower, be subject to modeling impact comparatively large, commercial production is difficult to practical.

Conventional Zoeppritz equation computational reflect coefficient all based on single interface, without gauge variation, although also AVO feature can be calculated, can not Direct Analysis thickness to the impact effect of each frequency component.

Summary of the invention

For the deficiency that prior art exists, the object of the invention is to provide a kind of computing method of Elastic Multi-layer medium midplane wave reflection coefficient, take into full account that elastic layer system multiple reflection and conversion involve thickness, frequency to the impact of reflection coefficient, be suitable for thin layer AVO analysis mode, counting yield is high, and simulation precision is high.

To achieve these goals, the present invention realizes by the following technical solutions:

The computing method of a kind of Elastic Multi-layer medium midplane wave reflection coefficient of the present invention, specifically comprise following step:

(1) plane harmonization P is established _efrom medium n+1 to object series of strata with incidence, then produce reflected P-wave and reflection wave in n+1 medium, and its reflection coefficient is respectively P _rand S _r, in medium 1, produce transmitted P-wave and transmitted shear wave, its transmission coefficient is respectively P _tand S _t, each layer medium all can be write as scalar potential and vector potential form;

(2) under two-dimensional case, the relational expression between displacement, stress and displacement position is determined;

(3) relating to scalar potential in step (1) and vector potential is brought in step (2), can obtain in (n+1)th layer and the displacement of the 1st layer and stress relation formula;

(4) according in (n+1)th layer that obtains in (3) and the displacement of the 1st layer and stress relation formula, the displacement, the Stress transmit matrix that comprise each layer elasticity coefficient can be obtained, then it is solved, reflection coefficient, transmission coefficient can be obtained.

Step (1) specifically comprises the following steps:

If interlayer comprises n-1 flat seam, sequence is numbered 2,3 ... n, bottom medium is 1 layer; In each layer, compressional wave and transverse wave speed use the α of subscripting respectively, and β represents, under be designated as sequence number, get x coordinate axis and n-th layer bottom boundary coincides; If a ripple incides interlayer end face from n+1 layer direction, note compressional wave incident wave is P _e, now, in n+1 layer, have a longitudinal wave reflection ripple P _r, transverse wave reflection ripple S _r, in 1 layer, have compressional wave transmitted wave P _twith shear wave transmitted wave S _tif compressional wave and shear wave are to each layer incident angle

Then total in n+1 layer medium bit shift is:

${\mathrm{\φ}}_{n+1}={\mathrm{\φ}}_e+{\mathrm{\φ}}_r=(A_1^{n+1}{e}^{-{\mathrm{jd}}^{n+1}z}+A_2^{n+1}{e}^{{\mathrm{jd}}^{n+1}z})e^{j(\mathrm{\σ}x-\mathrm{\ω}t)}$

${\mathrm{\ψ}}_{n+1}={\mathrm{\ψ}}_r=B_2^{n+1}{e}^{{\mathrm{js}}^{n+1}z}{e}^{j(\mathrm{\σ}x-\mathrm{\ω}t)}---(2-1-1)$

Wherein: $d^{n+1}=k^{n+1}{\mathrm{cosi}}_d^{n+1}=\sqrt{{\left(k^{n+1}\right)}^2-{\mathrm{\σ}}^2},s^{n+1}=K^{n+1}{\mathrm{cosi}}_s^{n+1}=\sqrt{{\left(K^{n+1}\right)}^2-{\mathrm{\σ}}^2},$ σ＝ω/c， k ⁿ⁺¹＝ω/α _n+1，K ⁿ⁺¹＝ω/β _n+1

Wherein, ω is frequency, φ _efor incident longitudinal wave bit function, φ _rreflected P-wave bit function, be amplitude, the e of (n+1) layer incident longitudinal wave be constant, j is imaginary unit, z is direction, for the amplitude of reflected P-wave, x be direction, t is time, ψ _rrepresent the bit function of reflection wave, for the amplitude of reflection wave;

C is the apparent velocity of ripple along interphase direction, according to Snell's law, is all equal to each ripple on interphase;

In n layer medium, compressional wave and shear wave bit function expression formula are respectively:

${\mathrm{\φ}}_n={\mathrm{\φ}}_e+{\mathrm{\φ}}_r=(A_1^n{e}^{-{\mathrm{jd}}^{n}z}+A_2^n{e}^{{\mathrm{jd}}^{n}z})e^{j(\mathrm{\σ}x-\mathrm{\ω}t)}$

${\mathrm{\ψ}}_n={\mathrm{\ψ}}_e+{\mathrm{\ψ}}_r=B_1^n{e}^{-{\mathrm{js}}^{n}z}{e}^{j(\mathrm{\σ}x-\mathrm{\ω}t)}+B_2^n{e}^{{\mathrm{js}}^{n}z}{e}^{j(\mathrm{\σ}x-\mathrm{\ω}t)}---(2-1-2)$

$d^n=k^n{\mathrm{cosi}}_d^n=\sqrt{{\left(k^n\right)}^2-{\mathrm{\σ}}^2}---(2-1-3)$

$s^n=K^n{\mathrm{cosi}}_s^n=\sqrt{{\left(K^n\right)}^2-{\mathrm{\σ}}^2}---(2-1-4)$

Wherein, for incident longitudinal wave amplitude, for reflected P-wave amplitude, ψ _efor incident shear wave bit function, for incident shear wave amplitude, for reflection wave amplitude, k ⁿthe horizontal wave number of n-th layer compressional wave, K ⁿthe horizontal wave number of n-th layer shear wave;

Then in 1 layer of medium, transmitted wave bit function expression formula is respectively:

${\mathrm{\φ}}_1={\mathrm{\φ}}_t=A_1^1{e}^{-{\mathrm{jd}}^{1}z}{e}^{j(\mathrm{\σ}x-\mathrm{\ω}t)}$

${\mathrm{\ψ}}_1={\mathrm{\ψ}}_t=B_1^1{e}^{-{\mathrm{js}}^{1}z}{e}^{j(\mathrm{\σ}x-\mathrm{\ω}t)}---(2-1-5)$

Wherein, $d^1=\sqrt{{\left(k^1\right)}^2-{\mathrm{\σ}}^2},s^1=\sqrt{{\left(K^1\right)}^2-{\mathrm{\σ}}^2},$

φ _tfor transmitted P-wave bit function, for transmitted P-wave amplitude, ψ _tfor transmitted shear wave bit function, for transmitted shear wave amplitude, k ¹for the horizontal wave number of ground floor compressional wave, K ¹for the horizontal wave number of ground floor shear wave.

In step (2), in two-dimensional case, the relational expression of displacement, stress and displacement position is:

$u=\frac{\∂\mathrm{\φ}}{\∂x}-\frac{\∂\mathrm{\ψ}}{\∂z}$

$\mathrm{\ω}=\frac{\∂\mathrm{\φ}}{\∂z}+\frac{\∂\mathrm{\ψ}}{\∂x}$

${\mathrm{\σ}}_{zz}=\mathrm{\λ}(\frac{\∂u}{\∂x}+\frac{\∂\mathrm{\ω}}{\∂z})+2\mathrm{\μ}\frac{\∂\mathrm{\ω}}{\∂z}$

${\mathrm{\τ}}_{zx}=\mathrm{\μ}(\frac{\∂u}{\∂z}+\frac{\∂\mathrm{\ω}}{\∂x})---(2-1-6)$

Wherein, u represents displacement, σ _zzrepresent principal strain, the τ along z-axis _zxrepresent shear stress, z is coordinate direction, λ and μ Lame's constant.

Step (3) specifically comprises the following steps:

If n layer thickness is h, calculates the displacement component in n layer and the components of stress, get the value of its z=h, be the displacement on n-th layer end face and components of stress value, be designated as u ⁽ⁿ⁾, ω ⁽ⁿ⁾, bring formula (2-1-2), (2-1-3), (2-1-4) into (2-1-6) formula, can matrix be obtained: wherein p=d ⁿh, Q=s ⁿh, d ⁿfor n-th layer compressional wave vertical wavenumber, h are thickness, s ⁿfor n-th layer shear wave vertical wavenumber;

$\left(\begin{array}{c}{u}^{\left(n\right)}\\ {\mathrm{\ω}}^{\left(n\right)}\\ {\mathrm{\σ}}_{zz}^{\left(n\right)}\\ {\mathrm{\τ}}_{zx}^{\left(n\right)}\end{array}\right)=\left(B_{ij}\right)\left(\begin{array}{c}{A}_2^n+A_1^n\\ A_2^n-A_1^n\\ B_2^n-B_1^n\\ B_2^n+B_1^n\end{array}\right)e^{j(\mathrm{\σ}x-\mathrm{\ω}t)}---(2-1-7)$

Wherein

$\left(B_{ij}\right)=\left(\begin{array}{cccc}j\mathrm{\σ}\cos p& -\mathrm{\σ}\sin p& -js\cos Q& s\sin Q\\ -d\sin p& jd\cos p& -\mathrm{\σ}\sin Q& \mathrm{j\σ}cosQ\\ -({\mathrm{\λ}}_n{k}_0^2+2{\mathrm{\μ}}_n{d}^2)\cos p& -j({\mathrm{\λ}}_n{k}_0^2+2{\mathrm{\μ}}_n{d}^2)\sin p& -2{\mathrm{\μ}}_n\mathrm{\σ}s\cos Q& -2{\mathrm{j\μ}}_n\mathrm{\σ}s\sin Q\\ -jd\mathrm{\σ}\sin p& -d\mathrm{\σ}\cos p& \frac{j}{2}(s^2-{\mathrm{\σ}}^2)\sin Q& \frac{1}{2}(s^2-{\mathrm{\σ}}^2)cosQ\end{array}\right)$

$k_0^2={\mathrm{\σ}}^2+d^2---(2-1-8)$

Wherein, s is the 1st layer of shear wave vertical wavenumber and d is the 1st layer of compressional wave vertical wavenumber, λ _nand μ _nfor the Lame's constant of n-th layer

Fetch bit moves and the value of the components of stress at z=0 place, can obtain the displacement component on n layer bottom surface and the components of stress, and according to interphase displacement and the stress condition of continuity, it should be equal with the analog value of (n-1) layer end face, is designated as u ^(n-1), ω ^(n-1), as z=0, p=0, Q=0, can be obtained by matrix (2-1-8)

${\left(B_{ij}\right)}_{p=0,Q=0}=\left(\begin{array}{cccc}j\mathrm{\σ}& 0& -js& 0\\ 0& jd& 0& j\mathrm{\σ}\\ -({\mathrm{\λ}}_n{k}_0^2+2{\mathrm{\μ}}_n{d}^2)& 0& -2{\mathrm{\μ}}_n\mathrm{\σ}s& 0\\ 0& -d\mathrm{\σ}& 0& \frac{1}{2}(s^2-{\mathrm{\σ}}^2)\end{array}\right)---(2-1-9)$

(n-1) displacement on layer end face and components of stress value can be expressed as:

$\left(\begin{array}{c}{u}^{(n-1)}\\ {\mathrm{\ω}}^{(n-1)}\\ {\mathrm{\σ}}_{zz}^{(n-1)}\\ {\mathrm{\τ}}_{zx}^{(n-1)}\end{array}\right)={\left(B_{ij}\right)}_{Q=0}^{p=0}\left(\begin{array}{c}{A}_2^n+A_1^n\\ A_2^n-A_1^n\\ B_2^n-B_1^n\\ B_2^n+B_1^n\end{array}\right)e^{j(\mathrm{\σ}x-\mathrm{\ω}t)}---(2-1-10)$

Can set up the relation between the displacement component of (n) layer and (n-1) layer end face and the components of stress by above-mentioned formula, for this reason, solving equation group (2-1-10) can have:

$\left(\begin{array}{c}{A}_2^n+A_1^n\\ A_2^n-A_1^n\\ B_2^n-B_1^n\\ B_2^n+B_1^n\end{array}\right)e^{j(\mathrm{\σ}x-\mathrm{\ω}t)}=\left(b_{ij}\right)\left(\begin{array}{c}{u}^{(n-1)}\\ {\mathrm{\ω}}^{(n-1)}\\ {\mathrm{\σ}}_{zz}^{(n-1)}\\ \frac{1}{2{\mathrm{\μ}}_{n-1}}{\mathrm{\τ}}_{zx}^{(n-1)}\end{array}\right)---(2-1-11)$

Wherein (b _ij) represent inverse matrix:

$\left(b_{ij}\right)=\left(\begin{array}{cccc}\frac{-2j\mathrm{\σ}}{K^2}& 0& -\frac{1}{{\mathrm{\μ}}_n{K}^2}& 0\\ 0& \frac{-j(s^2-{\mathrm{\σ}}^2)}{{\mathrm{dK}}^2}& 0& \frac{-2\mathrm{\σ}}{{\mathrm{dK}}^2}\\ \frac{j({\mathrm{\λ}}_n{k}_0^2/{\mathrm{\μ}}_n+2d^2)}{{\mathrm{sK}}^2}& 0& \frac{-\mathrm{\σ}}{{\mathrm{\μ}}_n{K}^{2}s}& 0\\ 0& \frac{-2j\mathrm{\σ}}{K^2}& 0& \frac{2}{K^2}\end{array}\right)$

Wherein, K is the horizontal wave number of shear wave

Formula (2-1-11) is substituted into (2-1-7) can obtain

$\left(\begin{array}{c}{u}^{\left(n\right)}\\ {\mathrm{\ω}}^{\left(n\right)}\\ {\mathrm{\σ}}_{zz}^{\left(n\right)}\\ \frac{1}{2{\mathrm{\μ}}_n}{\mathrm{\τ}}_{zx}^{\left(n\right)}\end{array}\right)=\left(B_{ij}\right)\left(b_{ij}\right)\left(\begin{array}{c}{u}^{(n-1)}\\ {\mathrm{\ω}}^{(n-1)}\\ {\mathrm{\σ}}_{zz}^{(n-1)}\\ \frac{1}{2{\mathrm{\μ}}_{n-1}}{\mathrm{\τ}}_{zx}^{(n-1)}\end{array}\right)---(2-1-12)$

Get mark (a _ij) representing matrix product (B _ij) (b _ij), be similarly 4 × 4 rank square formations, each element is:

a ₁₁＝2sin ²γcos p-cos2γcos Q

a ₁₂＝j(tanθcos2γsin p+sin2γsin Q)

$a_{13}=\frac{jsin\mathrm{\θ}}{\mathrm{\ρ}\mathrm{\α}\mathrm{\ω}}(\cos Q-\cos p)$

$a_{14}=\frac{2\mathrm{\β}}{\mathrm{\ω}}(tan\mathrm{\θ}sin\mathrm{\γ}\sin p-cos\mathrm{\γ}\sin Q)$

$a_{21}=j(\frac{\mathrm{\β}cos\mathrm{\θ}}{\mathrm{\α}cos\mathrm{\γ}}sin2\mathrm{\γ}\sin p-tan\mathrm{\γ}cos2\mathrm{\γ}\sin Q)$

a ₂₂＝cos2γcos p+2sin ²γcos Q

$a_{23}=\frac{1}{\mathrm{\ρ}\mathrm{\α}\mathrm{\ω}}(cos\mathrm{\θ}\sin p+tan\mathrm{\γ}sin\mathrm{\θ}\sin Q)$

$a_{24}=\frac{2j\mathrm{\β}}{\mathrm{\ω}}sin\mathrm{\γ}(\cos Q-\cos p)$

a ₃₁＝2jωρβsinγ(cosp-cosQ)cos2γ

$a_{32}=-\mathrm{\ρ}\mathrm{\ω}(\frac{{\mathrm{\αcos}}^{2}2\mathrm{\γ}}{cos\mathrm{\θ}}\sin p+4{\mathrm{\βcos\γsin}}^2\mathrm{\γ}\sin Q)$

a ₃₃＝cos2γcos p+2sin ²γcos Q

a ₃₄＝2jρβ ²(cos2γtanθsin p-sin2γsin Q)

$a_{41}=-\mathrm{\ω}(\frac{2}{\mathrm{\α}}{\sin }^2\mathrm{\γ}cos\mathrm{\θ}sinp+\frac{1}{2\mathrm{\β}}\frac{{\cos }^{2}2\mathrm{\γ}}{\cos \mathrm{\γ}}\sin Q)$

$a_{42}=\frac{j\mathrm{\ω}}{\mathrm{\α}}sin\mathrm{\θ}cos2\mathrm{\γ}(\cos p-\cos Q)$

$a_{43}=\frac{j}{2\mathrm{\ρ}}(\frac{sin2\mathrm{\θ}}{{\mathrm{\α}}^2}\sin p-\frac{cos2\mathrm{\γ}}{{\mathrm{\β}}^2}tan\mathrm{\γ}\sin Q)$

a ₄₄＝2sin ²γcos p+cos2γcos Q (2-1-13)

Wherein, γ=i _s, θ=i _drepresent the incident angle of ripple in layer respectively, so there is sin γ=σ/K, sin θ=σ/k, α, β is p-and s-wave velocity, ρ is Media density, establishes (n) layer and the relation between (n-1) layer end face displacement component and the components of stress in formula (2-1-12), calculates (a _ij) matrix element, will use the parameter values of (n) layer at formula (2-1-13), for this reason, in formula (2-1-12), matrix of coefficients corner brace (n) represents to have:

$\left(\begin{array}{c}{u}^{\left(n\right)}\\ {\mathrm{\ω}}^{\left(n\right)}\\ {\mathrm{\σ}}_{zz}^{\left(n\right)}\\ \frac{1}{2{\mathrm{\μ}}_n}{\mathrm{\τ}}_{zx}^{\left(n\right)}\end{array}\right)=\left(a_{ij}^n\right)\left(\begin{array}{c}{u}^{(n-1)}\\ {\mathrm{\ω}}^{(n-1)}\\ {\mathrm{\σ}}_{zz}^{(n-1)}\\ \frac{1}{2{\mathrm{\μ}}_{n-1}}{\mathrm{\τ}}_{zx}^{(n-1)}\end{array}\right)---(2-1-14)$

Obtain the recursion formula of displacement component on each interphase of interlayer and the components of stress, meet relational expression (2-1-14) and be equivalent to and meet displacement on interlayer interphase and the stress condition of continuity.

Step (4) specifically comprises the following steps:

In order to the transmission coefficient determining interlayer end face reflection coefficient and propagate below interlayer bottom surface through interlayer, set up the relation of displacement component on (n) layer end face and (1) layer end face and the components of stress according to formula (2-1-14):

$\left(\begin{array}{c}{u}^{\left(n\right)}\\ {\mathrm{\ω}}^{\left(n\right)}\\ {\mathrm{\σ}}_{zz}^{\left(n\right)}\\ \frac{1}{2{\mathrm{\μ}}_n}{\mathrm{\τ}}_{zx}^{\left(n\right)}\end{array}\right)=\left(a_{ij}^n\right)\left(a_{ij}^{n-1}\right)\mathrm{...}\left(a_{ij}^2\right)\left(\begin{array}{c}{u}^{\left(1\right)}\\ {\mathrm{\ω}}^{\left(1\right)}\\ {\mathrm{\σ}}_{zz}^{\left(1\right)}\\ \frac{1}{2{\mathrm{\μ}}_1}{\mathrm{\τ}}_{zx}^{\left(1\right)}\end{array}\right)---(2-1-15)$

(2-1-15) is out of shape

$\left(\begin{array}{c}{u}^{\left(n\right)}\\ {\mathrm{\ω}}^{\left(n\right)}\\ {\mathrm{\σ}}_{zz}^{\left(n\right)}\\ {\mathrm{\τ}}_{zx}^{\left(n\right)}\end{array}\right)=\left(\begin{array}{cccc}{A}_{11}& A_{12}& A_{13}& A_{14}\\ A_{21}& A_{22}& A_{23}& A_{24}\\ A_{31}& A_{32}& A_{33}& A_{34}\\ A_{41}& A_{42}& A_{43}& A_{44}\end{array}\right)\left(\begin{array}{c}{u}^{\left(1\right)}\\ {\mathrm{\ω}}^{\left(1\right)}\\ {\mathrm{\σ}}_{zz}^{\left(1\right)}\\ {\mathrm{\τ}}_{zx}^{\left(1\right)}\end{array}\right)---(2-1-16)$

Formula (2-1-1), formula (2-1-5) are substituted into formula (2-1-16) and obtain matrix equation (2-1-17), reflection coefficient and the transmission coefficient of interlayer after solving, can be obtained;

Consider A _ijbe plural number, make A ₁₁=R ₁₁+ iI ₁₁, A ₁₂=R ₁₂+ iI ₁₂, A ₁₃=R ₁₃+ iI ₁₃,

$\left(\begin{array}{c}i\mathrm{\σ}\\ -i{d}^{n+1}\\ -{\mathrm{\ρ}}_{n+1}{\mathrm{\β}}_{n+1}(K_{n+1}^2-2{\mathrm{\σ}}^2)\\ 2{\mathrm{\ρ}}_{n+1}{\mathrm{\β}}_{n+1}^2\mathrm{\σ}{d}^{n+1}\end{array}\right)=\left(\begin{array}{cccc}{B}_{11}& B_{12}& B_{13}& B_{14}\\ B_{21}& B_{22}& B_{23}& B_{24}\\ B_{31}& B_{32}& B_{33}& B_{34}\\ B_{41}& B_{42}& B_{43}& B_{44}\end{array}\right)\left(\begin{array}{c}V\\ V^{/}\\ W^{/}\\ W\end{array}\right)---(2-1-17)$

Wherein, ρ _n+1the density of (n+1)th layer, the horizontal wave number of shear wave of (n+1)th layer

Wherein B is called dematrix, and its each element is:

B ₁₁＝B ₂₂＝-iσ，B ₁₂＝is _n+1

B ₂₁＝id _n+1， $B_{31}=-B_{42}={\mathrm{\ρ}}_{n+1}{\mathrm{\β}}_{n+1}^2(K_{n+1}^2-2{\mathrm{\σ}}^2)$

$B_{32}=2{\mathrm{\ρ}}_{n+1}{\mathrm{\β}}_{n+1}^2{\mathrm{\σs}}_{n+1},B_{41}=2{\mathrm{\ρ}}_{n+1}{\mathrm{\β}}_{n+1}^2{\mathrm{\σd}}_{n+1}$

$B_{n4}=-D_{n1}\mathrm{\σ}+D_{n2}{d}_1-C_{n3}{\mathrm{\ρ}}_1{\mathrm{\β}}_1^2(K_1^2-2{\mathrm{\σ}}^2)+C_{n4}2{\mathrm{\ρ}}_1{\mathrm{\β}}_1^2{\mathrm{\σd}}_1(n=1,2,3,4)$

$\begin{array}{c}{B}_{n3}=\[-D_{n1}{s}_1-D_{n2}\mathrm{\σ}+C_{n3}2{\mathrm{\ρ}}_1{{\mathrm{\β}}_1}^2{\mathrm{\σs}}_1{+C}_{n4}{\mathrm{\ρ}}_1{{\mathrm{\β}}_1}^2(K_1^2-2{\mathrm{\σ}}^2)\]+\\ +i\[C_{n1}{s}_1+C_{n2}\mathrm{\σ}+D_{n3}2{\mathrm{\ρ}}_1{{\mathrm{\β}}_1}^2{\mathrm{\σs}}_1{+D}_{n4}{\mathrm{\ρ}}_1{{\mathrm{\β}}_1}^2(K_1^2-2{\mathrm{\σ}}^2)\](n=1,2,3,4)\end{array}$

Wherein, D _n1for plural number;

A _n1imaginary part, D _n2for plural A _n2imaginary part, d ₁for thickness, C _n3for plural A _n3imaginary part, ρ ₁for ground floor density, represent horizontal wave number, the C of ground floor shear wave _n4for plural A _n4imaginary part, s ₁ground floor shear wave vertical wavenumber, C _n1for plural A _n1imaginary part, C _n2for plural A _n2imaginary part solve dematrix expression formula formula (2-2-17) the longitudinal wave reflection coefficients R of plane harmonization at elastic layer system can be obtained _pp, transverse wave reflection coefficient V ^/, shear wave transmission coefficient W ^/, compressional wave transmission coefficient W.

The present invention is based on the reflection coefficient spectral theory of elastic layer system, the displacement of application elastic layer system and stress recursion formula, for solid multi-layer medium, be deduced vertical, transverse wave reflection, transmission coefficient formula, the method can calculate the Reflection Coefficient of Planar Wave of different incidence angles degree, on petroleum exploration, geologic thin (mutually) layer stratal configuration can be represented with elastic layer system model, therefore method of the present invention may be used for calculating thin (mutually) layer AVO (Amplitude versus offset), multiple reflection can be solved, transformed wave is on the impact of thin layer AVO, simultaneously can quantitative simulation different reservoir thickness, different frequency changes the AVO change caused, may be used for explanation of seismic road set information and routine business software based on the unmatched problem of well zoeppritz equation Zheng Yan road collection simultaneously, the present invention can as an important supplement of current oil exploration industry commercial software.

Accompanying drawing explanation

Fig. 1 is reflection and the transmission model of layered medium;

Fig. 2 is reflection and the transmission of single thin layer;

Fig. 3 is the reservoir AVO characteristic simulation under different reservoir depth information;

Fig. 4 is the reservoir AVO characteristic simulation under different incident wave frequency content;

Fig. 5 is that single-layer medium conforms to traditional commerce software.

Embodiment

The technological means realized for making the present invention, creation characteristic, reaching object and effect is easy to understand, below in conjunction with embodiment, setting forth the present invention further.

Conventional Zoeppritz equation computational reflect coefficient all based on single interface, without gauge variation, although also AVO feature can be calculated, can not Direct Analysis thickness to the impact effect of each frequency component.The Reflection Coefficient of Planar Wave of Elastic Multi-layer medium has taken into full account that elastic layer system multiple reflection and conversion involve thickness, frequency to the impact of reflection coefficient, is therefore more suitable for studying thin layer, the AVO of thin interbed analyzes and forward simulation.Therefore the present invention describes Reflection Coefficient of Planar Wave in Elastic Multi-layer medium and forward modelling method in detail, and has carried out deriving checking to single-layer model.

Reflection Coefficient of Planar Wave computing method in Elastic Multi-layer medium are as follows:

An interlayer in the middle of two semiinfinite elastic mediums.This interlayer is made up of multiple horizontal level.Be provided with a plane wave from the incident interlayer end face of upper dielectric, then produce an interlayer reflection wave and propagate in upper dielectric, ripple also through interlayer, will produce a transmitted wave at bottom Propagation simultaneously.Calculate reflection coefficient and the transmission coefficient of interlayer.

As shown in Figure 1, interlayer is made up of n-1 flat seam, sequence is numbered 2,3 ... n; Bottom medium is 1 layer.In each layer, compressional wave and transverse wave speed use the α of subscripting respectively, and β represents, under be designated as sequence number.Get x coordinate axis and n-th layer bottom boundary coincides.Have a P-SV wave system system to incide interlayer end face from n+1 layer direction, note compressional wave incident wave is P _e, now, in n+1 layer, have a longitudinal wave reflection ripple P _r, transverse wave reflection ripple S _r; Compressional wave transmitted wave P is had in 1 layer _twith shear wave transmitted wave S _t.If compressional wave and shear wave are to each layer incident angle

Then total in n+1 layer medium bit shift is:

${\mathrm{\φ}}_{n+1}={\mathrm{\φ}}_e+{\mathrm{\φ}}_r=(A_1^{n+1}{e}^{-{\mathrm{jd}}^{n+1}z}+A_2^{n+1}{e}^{{\mathrm{jd}}^{n+1}z})e^{j(\mathrm{\σ}x-\mathrm{\ω}t)}$

${\mathrm{\ψ}}_{n+1}={\mathrm{\ψ}}_r=B_2^{n+1}{e}^{{\mathrm{js}}^{n+1}z}{e}^{j(\mathrm{\σ}x-\mathrm{\ω}t)}---(2-1-1)$

Wherein: $d^{n+1}=k^{n+1}{\mathrm{cosi}}_d^{n+1}=\sqrt{{\left(k^{n+1}\right)}^2-{\mathrm{\σ}}^2},s^{n+1}=K^{n+1}{\mathrm{cosi}}_s^{n+1}=\sqrt{{\left(K^{n+1}\right)}^2-{\mathrm{\σ}}^2},$ σ＝ω/c， k ⁿ⁺¹＝ω/α _n+1，K ⁿ⁺¹＝ω/β _n+1

C is the apparent velocity of ripple along interphase direction, according to Snell's law, is all equal to each ripple on interphase.

In n layer medium, compressional wave and shear wave bit function expression formula are respectively:

${\mathrm{\φ}}_n={\mathrm{\φ}}_e+{\mathrm{\φ}}_r=(A_1^n{e}^{-{\mathrm{jd}}^{n}z}+A_2^n{e}^{{\mathrm{jd}}^{n}z})e^{j(\mathrm{\σ}x-\mathrm{\ω}t)}$

${\mathrm{\ψ}}_n={\mathrm{\ψ}}_e+{\mathrm{\ψ}}_r=B_1^n{e}^{-{\mathrm{js}}^{n}z}{e}^{j(\mathrm{\σ}x-\mathrm{\ω}t)}+B_2^n{e}^{{\mathrm{js}}^{n}z}{e}^{j(\mathrm{\σ}x-\mathrm{\ω}t)}---(2-1-2)$

$d^n=k^n{\mathrm{cosi}}_d^n=\sqrt{{\left(k^n\right)}^2-{\mathrm{\σ}}^2}---(2-1-3)$

$s^n=K^n{\mathrm{cosi}}_s^n=\sqrt{{\left(K^n\right)}^2-{\mathrm{\σ}}^2}---(2-1-4)$

Then in 1 layer of medium, transmitted wave bit function expression formula is respectively:

${\mathrm{\φ}}_1={\mathrm{\φ}}_t=A_1^1{e}^{-{\mathrm{jd}}^{1}z}{e}^{j(\mathrm{\σ}x-\mathrm{\ω}t)}$

${\mathrm{\ψ}}_1={\mathrm{\ψ}}_t=B_1^1{e}^{-{\mathrm{js}}^{1}z}{e}^{j(\mathrm{\σ}x-\mathrm{\ω}t)}---(2-1-5)$

Wherein $d^1=\sqrt{{\left(k^1\right)}^2-{\mathrm{\σ}}^2},s^1=\sqrt{{\left(K^1\right)}^2-{\mathrm{\σ}}^2}$

In two-dimensional case, the relational expression of displacement, stress and displacement position is:

$u=\frac{\∂\mathrm{\φ}}{\∂x}-\frac{\∂\mathrm{\ψ}}{\∂z}$

$\mathrm{\ω}=\frac{\∂\mathrm{\φ}}{\∂z}+\frac{\∂\mathrm{\ψ}}{\∂x}$

${\mathrm{\σ}}_{zz}=\mathrm{\λ}(\frac{\∂u}{\∂x}+\frac{\∂\mathrm{\ω}}{\∂z})+2\mathrm{\μ}\frac{\∂\mathrm{\ω}}{\∂z}$

${\mathrm{\τ}}_{zx}=\mathrm{\μ}(\frac{\∂u}{\∂z}+\frac{\∂\mathrm{\ω}}{\∂x})---(2-1-6)$

If n layer thickness is h, calculates the displacement component in n layer and the components of stress, get the value of its z=h, be the displacement on n-th layer end face and components of stress value, be designated as u ⁽ⁿ⁾, ω ⁽ⁿ⁾, bring formula (2-1-2), (2-1-3), (2-1-4) into (2-1-6) formula, can matrix be obtained: wherein p=d ⁿh, Q=s ⁿh,

$\left(\begin{array}{c}{u}^{\left(n\right)}\\ {\mathrm{\ω}}^{\left(n\right)}\\ {\mathrm{\σ}}_{zz}^{\left(n\right)}\\ {\mathrm{\τ}}_{zx}^{\left(n\right)}\end{array}\right)=\left(B_{ij}\right)\left(\begin{array}{c}{A}_2^n+A_1^n\\ A_2^n-A_1^n\\ B_2^n-B_1^n\\ B_2^n+B_1^n\end{array}\right)e^{j(\mathrm{\σ}x-\mathrm{\ω}t)}---(2-1-7)$

Wherein

$\left(B_{ij}\right)=\left(\begin{array}{cccc}j\mathrm{\σ}\cos p& -\mathrm{\σ}\sin p& -js\cos Q& s\sin Q\\ -d\sin p& jd\cos p& -\mathrm{\σ}\sin Q& \mathrm{j\σ}cosQ\\ -({\mathrm{\λ}}_n{k}_0^2+2{\mathrm{\μ}}_n{d}^2)\cos p& -j({\mathrm{\λ}}_n{k}_0^2+2{\mathrm{\μ}}_n{d}^2)\sin p& -2{\mathrm{\μ}}_n\mathrm{\σ}s\cos Q& -2{\mathrm{j\μ}}_n\mathrm{\σ}s\sin Q\\ -jd\mathrm{\σ}\sin p& -d\mathrm{\σ}\cos p& \frac{j}{2}(s^2-{\mathrm{\σ}}^2)\sin Q& \frac{1}{2}(s^2-{\mathrm{\σ}}^2)cosQ\end{array}\right)$

$k_0^2={\mathrm{\σ}}^2+d^2---(2-1-8)$

On the one hand, fetch bit to move with the components of stress in the value at z=0 place, can obtain the displacement component on n layer bottom surface and the components of stress in order.According to interphase displacement and the stress condition of continuity, it should be equal with the analog value of (n-1) layer end face, is designated as u ^(n-1), ω ^(n-1), as z=0, p=0, Q=0, can be obtained by matrix (2-1-8)

${\left(B_{ij}\right)}_{p=0,Q=0}=\left(\begin{array}{cccc}j\mathrm{\σ}& 0& -js& 0\\ 0& jd& 0& j\mathrm{\σ}\\ -({\mathrm{\λ}}_n{k}_0^2+2{\mathrm{\μ}}_n{d}^2)& 0& -2{\mathrm{\μ}}_n\mathrm{\σ}s& 0\\ 0& -d\mathrm{\σ}& 0& \frac{1}{2}(s^2-{\mathrm{\σ}}^2)\end{array}\right)---(2-1-9)$

(n-1) displacement on layer end face and components of stress value can be expressed as:

$\left(\begin{array}{c}{u}^{(n-1)}\\ {\mathrm{\ω}}^{(n-1)}\\ {\mathrm{\σ}}_{zz}^{(n-1)}\\ {\mathrm{\τ}}_{zx}^{(n-1)}\end{array}\right)={\left(B_{ij}\right)}_{Q=0}^{p=0}\left(\begin{array}{c}{A}_2^n+A_1^n\\ A_2^n-A_1^n\\ B_2^n-B_1^n\\ B_2^n+B_1^n\end{array}\right)e^{j(\mathrm{\σ}x-\mathrm{\ω}t)}---(2-1-10)$

The relation between the displacement component of (n) layer and (n-1) layer end face and the components of stress can be set up by above-mentioned formula.For this reason, solving equation group (2-1-10) can have:

$\left(\begin{array}{c}{A}_2^n+A_1^n\\ A_2^n-A_1^n\\ B_2^n-B_1^n\\ B_2^n+B_1^n\end{array}\right)e^{j(\mathrm{\σ}x-\mathrm{\ω}t)}=\left(b_{ij}\right)\left(\begin{array}{c}{u}^{(n-1)}\\ {\mathrm{\ω}}^{(n-1)}\\ {\mathrm{\σ}}_{zz}^{(n-1)}\\ \frac{1}{2{\mathrm{\μ}}_{n-1}}{\mathrm{\τ}}_{zx}^{(n-1)}\end{array}\right)$

$(2-1-11)$

Wherein (b _ij) represent inverse matrix:

$\left(b_{ij}\right)=\left(\begin{array}{cccc}\frac{-2j\mathrm{\σ}}{K^2}& 0& -\frac{1}{{\mathrm{\μ}}_n{K}^2}& 0\\ 0& \frac{-j(s^2-{\mathrm{\σ}}^2)}{{\mathrm{dK}}^2}& 0& \frac{-2\mathrm{\σ}}{{\mathrm{dK}}^2}\\ \frac{j({\mathrm{\λ}}_n{k}_0^2/{\mathrm{\μ}}_n+2d^2)}{{\mathrm{sK}}^2}& 0& \frac{-\mathrm{\σ}}{{\mathrm{\μ}}_n{K}^{2}s}& 0\\ 0& \frac{-2j\mathrm{\σ}}{K^2}& 0& \frac{2}{K^2}\end{array}\right)$

Formula (2-1-11) is substituted into (2-1-7) can obtain

$\left(\begin{array}{c}{u}^{\left(n\right)}\\ {\mathrm{\ω}}^{\left(n\right)}\\ {\mathrm{\σ}}_{zz}^{\left(n\right)}\\ \frac{1}{2{\mathrm{\μ}}_n}{\mathrm{\τ}}_{zx}^{\left(n\right)}\end{array}\right)=\left(B_{ij}\right)\left(b_{ij}\right)\left(\begin{array}{c}{u}^{(n-1)}\\ {\mathrm{\ω}}^{(n-1)}\\ {\mathrm{\σ}}_{zz}^{(n-1)}\\ \frac{1}{2{\mathrm{\μ}}_{n-1}}{\mathrm{\τ}}_{zx}^{(n-1)}\end{array}\right)---(2-1-12)$

Get mark (a _ij) representing matrix product (B _ij) (b _ij), be similarly 4 × 4 rank square formations, each element is:

a ₁₁＝2sin ²γcos p-cos2γcos Q

a ₁₂＝j(tanθcos2γsin p+sin2γsin Q)

$a_{13}=\frac{jsin\mathrm{\θ}}{\mathrm{\ρ}\mathrm{\α}\mathrm{\ω}}(\cos Q-\cos p)$

$a_{14}=\frac{2\mathrm{\β}}{\mathrm{\ω}}(tan\mathrm{\θ}sin\mathrm{\γ}\sin p-cos\mathrm{\γ}\sin Q)$

$a_{21}=j(\frac{\mathrm{\β}cos\mathrm{\θ}}{\mathrm{\α}cos\mathrm{\γ}}sin2\mathrm{\γ}\sin p-tan\mathrm{\γ}cos2\mathrm{\γ}\sin Q)$

a ₂₂＝cos2γcos p+2sin ²γcos Q

$a_{23}=\frac{1}{\mathrm{\ρ}\mathrm{\α}\mathrm{\ω}}(cos\mathrm{\θ}\sin p+tan\mathrm{\γ}sin\mathrm{\θ}\sin Q)$

$a_{24}=\frac{2j\mathrm{\β}}{\mathrm{\ω}}sin\mathrm{\γ}(\cos Q-\cos p)$

a ₃₁＝2jωρβsinγ(cosp-cosQ)cos2γ

$a_{32}=-\mathrm{\ρ}\mathrm{\ω}(\frac{{\mathrm{\αcos}}^{2}2\mathrm{\γ}}{cos\mathrm{\θ}}\sin p+4{\mathrm{\βcos\γsin}}^2\mathrm{\γ}\sin Q)$

a ₃₃＝cos2γcos p+2sin ²γcos Q

a ₃₄＝2jρβ ²(cos2γtanθsin p-sin2γsin Q)

$a_{41}=-\mathrm{\ω}(\frac{2}{\mathrm{\α}}{\sin }^2\mathrm{\γ}cos\mathrm{\θ}\sin p+\frac{1}{2\mathrm{\β}}\frac{{\cos }^{2}2\mathrm{\γ}}{\cos \mathrm{\γ}}\sin Q)$

$a_{42}=\frac{j\mathrm{\ω}}{\mathrm{\α}}sin\mathrm{\θ}cos2\mathrm{\γ}(\cos p-\cos Q)$

$a_{43}=\frac{j}{2\mathrm{\ρ}}(\frac{sin2\mathrm{\θ}}{{\mathrm{\α}}^2}\sin p-\frac{cos2\mathrm{\γ}}{{\mathrm{\β}}^2}tan\mathrm{\γ}\sin Q)$

a ₄₄＝2sin ²γcos p+cos2γcos Q (2-1-13)

Wherein γ=i _s, θ=i _drepresent the incident angle of ripple in layer respectively, so there is sin γ=σ/K, sin θ=σ/k.Other parameter, if α, β are p-and s-wave velocity, ρ is Media density, and ω is frequency.In formula (2-1-12), establish (n) layer and the relation between (n-1) layer end face displacement component and the components of stress, calculate (a _ij) matrix element, the parameter values of (n) layer will be used at formula (2-1-13).For this reason, in formula (2-1-12), matrix of coefficients corner brace (n) represents to have:

$\left(\begin{array}{c}{u}^{\left(n\right)}\\ {\mathrm{\ω}}^{\left(n\right)}\\ {\mathrm{\σ}}_{zz}^{\left(n\right)}\\ \frac{1}{2{\mathrm{\μ}}_n}{\mathrm{\τ}}_{zx}^{\left(n\right)}\end{array}\right)=\left(a_{ij}^n\right)\left(\begin{array}{c}{u}^{(n-1)}\\ {\mathrm{\ω}}^{(n-1)}\\ {\mathrm{\σ}}_{zz}^{(n-1)}\\ \frac{1}{2{\mathrm{\μ}}_{n-1}}{\mathrm{\τ}}_{zx}^{(n-1)}\end{array}\right)---(2-1-14)$

Obtain the recursion formula of displacement component on each interphase of interlayer and the components of stress.Meet relational expression (2-1-14) to be equivalent to and to meet displacement on interlayer interphase and the stress condition of continuity.

In order to the transmission coefficient determining interlayer end face reflection coefficient and propagate below interlayer bottom surface through interlayer, set up the relation of displacement component on (n) layer end face and (1) layer end face and the components of stress according to formula (2-1-14):

$\left(\begin{array}{c}{u}^{\left(n\right)}\\ {\mathrm{\ω}}^{\left(n\right)}\\ {\mathrm{\σ}}_{zz}^{\left(n\right)}\\ \frac{1}{2{\mathrm{\μ}}_n}{\mathrm{\τ}}_{zx}^{\left(n\right)}\end{array}\right)=\left(a_{ij}^n\right)\left(a_{ij}^{n-1}\right)\mathrm{...}\left(a_{ij}^2\right)\left(\begin{array}{c}{u}^{\left(1\right)}\\ {\mathrm{\ω}}^{\left(1\right)}\\ {\mathrm{\σ}}_{zz}^{\left(1\right)}\\ \frac{1}{2{\mathrm{\μ}}_1}{\mathrm{\τ}}_{zx}^{\left(1\right)}\end{array}\right)---(2-1-15)$

(2-1-15) is out of shape

$\left(\begin{array}{c}{u}^{\left(n\right)}\\ {\mathrm{\ω}}^{\left(n\right)}\\ {\mathrm{\σ}}_{zz}^{\left(n\right)}\\ {\mathrm{\τ}}_{zx}^{\left(n\right)}\end{array}\right)=\left(\begin{array}{cccc}{A}_{11}& A_{12}& A_{13}& A_{14}\\ A_{21}& A_{22}& A_{23}& A_{24}\\ A_{31}& A_{32}& A_{33}& A_{34}\\ A_{41}& A_{42}& A_{43}& A_{44}\end{array}\right)\left(\begin{array}{c}{u}^{\left(1\right)}\\ {\mathrm{\ω}}^{\left(1\right)}\\ {\mathrm{\σ}}_{zz}^{\left(1\right)}\\ {\mathrm{\τ}}_{zx}^{\left(1\right)}\end{array}\right)---(2-1-16)$

Formula (2-1-1), formula (2-1-5) are substituted into formula (2-1-16) and obtain matrix equation (2-1-17), reflection coefficient and the transmission coefficient of interlayer after solving, can be obtained.

Consider A _ijbe plural number, make A ₁₁=R ₁₁+ iI ₁₁, A ₁₂=R ₁₂+ iI ₁₂, A ₁₃=R ₁₃+ iI ₁₃,

$\left(\begin{array}{c}i\mathrm{\σ}\\ -i{d}^{n+1}\\ -{\mathrm{\ρ}}_{n+1}{\mathrm{\β}}_{n+1}(K_{n+1}^2-2{\mathrm{\σ}}^2)\\ 2{\mathrm{\ρ}}_{n+1}{\mathrm{\β}}_{n+1}^2\mathrm{\σ}{d}^{n+1}\end{array}\right)=\left(\begin{array}{cccc}{B}_{11}& B_{12}& B_{13}& B_{14}\\ B_{21}& B_{22}& B_{23}& B_{24}\\ B_{31}& B_{32}& B_{33}& B_{34}\\ B_{41}& B_{42}& B_{43}& B_{44}\end{array}\right)\left(\begin{array}{c}V\\ V^{/}\\ W^{/}\\ W\end{array}\right)---(2-1-17)$

Wherein B is called dematrix, and its each element is:

B ₁₁＝B ₂₂＝-iσ，B ₁₂＝is _n+1

B ₂₁＝id _n+1， $B_{31}=-B_{42}={\mathrm{\ρ}}_{n+1}{\mathrm{\β}}_{n+1}^2(K_{n+1}^2-2{\mathrm{\σ}}^2)$

$B_{32}=2{\mathrm{\ρ}}_{n+1}{\mathrm{\β}}_{n+1}^2{\mathrm{\σs}}_{n+1},B_{41}=2{\mathrm{\ρ}}_{n+1}{\mathrm{\β}}_{n+1}^2{\mathrm{\σd}}_{n+1}$

$B_{n4}=-D_{n1}\mathrm{\σ}+D_{n2}{d}_1-C_{n3}{\mathrm{\ρ}}_1{\mathrm{\β}}_1^2(K_1^2-2{\mathrm{\σ}}^2)+C_{n4}2{\mathrm{\ρ}}_1{\mathrm{\β}}_1^2{\mathrm{\σd}}_1(n=1,2,3,4)$

$\begin{array}{c}{B}_{n3}=\[-D_{n1}{s}_1-D_{n2}\mathrm{\σ}+C_{n3}2{\mathrm{\ρ}}_1{\mathrm{\β}}_1^2{\mathrm{\σs}}_1+C_{n4}{\mathrm{\ρ}}_1{\mathrm{\β}}_1^2(K_1^2-2{\mathrm{\σ}}^2)\]+\\ +i\[-C_{n1}{s}_1-C_{n2}\mathrm{\σ}+D_{n3}2{\mathrm{\ρ}}_1{\mathrm{\β}}_1^2{\mathrm{\σs}}_1+D_{n4}{\mathrm{\ρ}}_1{\mathrm{\β}}_1^2(K_1^2-2{\mathrm{\σ}}^2)\](n=1,2,3,4)\end{array}$

Solve dematrix expression formula formula (2-2-17) and the longitudinal wave reflection coefficients R of plane harmonization at elastic layer system can be obtained _pp, transverse wave reflection coefficient V ^/, shear wave transmission coefficient W ^/, compressional wave transmission coefficient W.

The achievement that the longitudinal wave reflection coefficient of Fig. 3 corresponding to single thin film model changes with angle, in figure, horizontal ordinate is angle, be incremented to 60 degree from 0 degree, ordinate is longitudinal wave reflection coefficient, changes between thickness of thin layer 5 meters to 30 meters, as can be seen from the figure the thin layer AVO plesiomorphism that different reservoir thickness is corresponding, but four the intercept of curve is different, and the gradient of four curves is also different, and this requires that we will consider the variation in thickness of thin layer when judging the AVO of thin layer.

The achievement that the longitudinal wave reflection coefficient of Fig. 4 corresponding to single thin film model changes with angle, in figure, horizontal ordinate is angle, 45 degree are incremented to from 0 degree, ordinate is longitudinal wave reflection coefficient, frequency changes by between 10hz to 90hz, as can be seen from the figure the thin layer AVO plesiomorphism that different frequency is corresponding, article 10, curve is all monotonically increasing, but 10 the gradient of curve is different, this requires that we will consider the dominant frequency change of incident wave when judging the AVO of thin layer, and especially same set of thin reservoir is when buried depth changes greatly.

For checking the correctness of above-mentioned theory, by the reflection coefficient in the single thin layer situation of theory deduction vertical incidence above.

If have a thickness to be the solid interlayer of h in the middle of two half infinite mediums, calculate reflection coefficient and the transmission coefficient of this interlayer.As shown in Figure 2, in the middle of 1 and 3 two half infinite medium, accompany a thin layer, its thickness is h.A compressional wave is had to impinge perpendicularly on thin layer end face from medium 3 direction, by the reflection wave p that generation one is propagated in medium 3 _rwith the transmitted wave p propagated in medium 1 _t, according to the coordinate system shown in figure, each ripple expression formula is:

${\mathrm{\φ}}_e={\mathrm{Ae}}^{-j\mathrm{\ω}t}{e}^{-j\mathrm{\ω}(z-h)/{\mathrm{\α}}_3}$

${\mathrm{\φ}}_r={\mathrm{Be}}^{-j\mathrm{\ω}t}{e}^{j\mathrm{\ω}(z-h)/{\mathrm{\α}}_3}$

${\mathrm{\φ}}_t={\mathrm{Ce}}^{-j\mathrm{\ω}t}{e}^{-j\mathrm{\ω}z/{\mathrm{\α}}_1}---(2-2-1)$

Bit function is had to the 3rd medium:

φ ₃＝φ _e+φ _r(2-2-2)

Bit function is had to the 1st medium:

φ ₁＝φ _t(2-2-3)

The boundary condition that they should meet, according to (2-1-15), for

$\left(\begin{array}{c}{\mathrm{\ω}}^{\left(2\right)}\\ {\mathrm{\σ}}_{zz}^{\left(2\right)}\end{array}\right)=\left(\begin{array}{cc}{a}_{22}& a_{23}\\ a_{32}& a_{33}\end{array}\right)\left(\begin{array}{c}{\mathrm{\ω}}^{\left(1\right)}\\ {\mathrm{\σ}}_{zz}^{\left(1\right)}\end{array}\right)---(2-2-4)$

When it considers vertical incidence, μ=0, τ _zx=0.By formula (2-2-2), formula (2-2-3) brings formula (2-2-4) formula into, and z=0 has:

${\mathrm{\ω}}^{\left(1\right)}=-j\frac{\mathrm{\ω}}{{\mathrm{\α}}_1}{\mathrm{ce}}^{-j\mathrm{\ω}t}$

${\mathrm{\σ}}_{zz}^{\left(1\right)}=-{\mathrm{\ρ}}_1{\mathrm{\ω}}^2{\mathrm{ce}}^{-j\mathrm{\ω}t}$

Z=h can have:

${\mathrm{\ω}}^{\left(3\right)}=-j\frac{\mathrm{\ω}}{{\mathrm{\α}}_3}(A-B)e^{-j\mathrm{\ω}t}$

${\mathrm{\σ}}_{zz}^{\left(3\right)}=-{\mathrm{\ρ}}_3{\mathrm{\ω}}^2(A+B)e^{-j\mathrm{\ω}t}$

Boundary condition equation group is:

$\left(\begin{array}{cc}\frac{j}{{\mathrm{\α}}_3}& a_{22}\frac{j}{{\mathrm{\α}}_1}+a_{23}{\mathrm{\ρ}}_1\mathrm{\ω}\\ -{\mathrm{\ρ}}_3\mathrm{\ω}& a_{32}\frac{j}{{\mathrm{\α}}_1}+a_{33}{\mathrm{\ρ}}_1\mathrm{\ω}\end{array}\right)\left(\begin{array}{c}\frac{B}{A}\\ \frac{C}{A}\end{array}\right)=\left(\begin{array}{c}\frac{j}{{\mathrm{\α}}_3}\\ {\mathrm{\ρ}}_3\mathrm{\ω}\end{array}\right)---(2-2-5)$

Solution of equations is exactly required reflection coefficient and transmission coefficient, and they are:

$R_{PP}=\frac{B}{A}=\frac{-a_{32}+{\mathrm{j\ω\ρ}}_1{\mathrm{\α}}_1{a}_{33}-{\mathrm{j\ω\ρ}}_3{\mathrm{\α}}_3{a}_{22}-{\mathrm{\ρ}}_3{\mathrm{\α}}_3{\mathrm{\ρ}}_1{\mathrm{\α}}_1{\mathrm{\ω}}^2{a}_{23}}{-a_{32}+{\mathrm{j\ω\ρ}}_1{\mathrm{\α}}_1{a}_{33}+{\mathrm{j\ω\ρ}}_3{\mathrm{\α}}_3{a}_{22}+{\mathrm{\ρ}}_3{\mathrm{\α}}_3{\mathrm{\ρ}}_1{\mathrm{\α}}_1{\mathrm{\ω}}^2{a}_{23}}---(2-2-6)$

$T_{PP}=\frac{{\mathrm{\α}}_3}{{\mathrm{\α}}_1}\frac{C}{A}=\frac{2{\mathrm{j\ω\ρ}}_3{\mathrm{\α}}_3}{-a_{32}+{\mathrm{j\ω\ρ}}_1{\mathrm{\α}}_1{a}_{33}+{\mathrm{j\ω\ρ}}_3{\mathrm{\α}}_3{a}_{22}+{\mathrm{\ρ}}_3{\mathrm{\α}}_3{\mathrm{\ρ}}_1{\mathrm{\α}}_1{\mathrm{\ω}}^2{a}_{23}}---(2-2-7)$

Wherein T _pPtransmission coefficient gets the displacement amplitude ratio of transmitted wave and incident wave.By thin-layered medium 2 parameter and incident angle data, θ=0, γ=0, bring formula (2-1-13) into, can have:

$a_{22}=cos\left(\frac{\mathrm{\ω}}{{\mathrm{\α}}_2}h\right)$

$a_{23}=\frac{1}{{\mathrm{\ρ}}_2{\mathrm{\α}}_2}\frac{1}{\mathrm{\ω}}sin\left(\frac{\mathrm{\ω}}{{\mathrm{\α}}_2}h\right)$

$a_{32}=-{\mathrm{\ρ}}_2{\mathrm{\α}}_2\mathrm{\ω}sin\left(\frac{\mathrm{\ω}}{{\mathrm{\α}}_2}h\right)$

$a_{33}=cos\left(\frac{\mathrm{\ω}}{{\mathrm{\α}}_2}h\right)---(2-2-8)$

Formula (2-2-8) is substituted into formula (2-2-6), and formula (2-2-7), can obtain

$R_{PP}=\frac{{\mathrm{\ρ}}_2{\mathrm{\α}}_2\sin \frac{\mathrm{\ω}}{{\mathrm{\α}}_2}h+{\mathrm{j\ρ}}_1{\mathrm{\α}}_1\cos \frac{\mathrm{\ω}}{{\mathrm{\α}}_2}h-{\mathrm{j\ρ}}_3{\mathrm{\α}}_3\cos \frac{\mathrm{\ω}}{{\mathrm{\α}}_2}h-\frac{{\mathrm{\ρ}}_3{\mathrm{\α}}_3{\mathrm{\ρ}}_1{\mathrm{\α}}_1}{{\mathrm{\ρ}}_2{\mathrm{\α}}_2}\sin \frac{\mathrm{\ω}}{{\mathrm{\α}}_2}h}{{\mathrm{\ρ}}_2{\mathrm{\α}}_2\sin \frac{\mathrm{\ω}}{{\mathrm{\α}}_2}h+{\mathrm{j\ρ}}_1{\mathrm{\α}}_1\cos \frac{\mathrm{\ω}}{{\mathrm{\α}}_2}h+{\mathrm{j\ρ}}_3{\mathrm{\α}}_3\cos \frac{\mathrm{\ω}}{{\mathrm{\α}}_2}h+\frac{{\mathrm{\ρ}}_3{\mathrm{\α}}_3{\mathrm{\ρ}}_1{\mathrm{\α}}_1}{{\mathrm{\ρ}}_2{\mathrm{\α}}_2}\sin \frac{\mathrm{\ω}}{{\mathrm{\α}}_2}h}---(2-2-9)$

$T_{PP}=\frac{2{\mathrm{j\ρ}}_3{\mathrm{\α}}_3}{{\mathrm{\ρ}}_2{\mathrm{\α}}_{2}sin\frac{\mathrm{\ω}}{{\mathrm{\α}}_2}h+{\mathrm{j\ρ}}_1{\mathrm{\α}}_{1}cos\frac{\mathrm{\ω}}{{\mathrm{\α}}_2}h+{\mathrm{j\ρ}}_3{\mathrm{\α}}_{3}cos\frac{\mathrm{\ω}}{{\mathrm{\α}}_2}h+\frac{{\mathrm{\ρ}}_3{\mathrm{\α}}_3{\mathrm{\ρ}}_1{\mathrm{\α}}_1}{{\mathrm{\ρ}}_2{\mathrm{\α}}_2}sin\frac{\mathrm{\ω}}{{\mathrm{\α}}_2}h}---(2-2-10)$

As h=0, reflection coefficient when can obtain an interphase during compressional wave vertical incidence and transmission coefficient:

$R_{PP}=\frac{{\mathrm{\ρ}}_1{\mathrm{\α}}_1-{\mathrm{\ρ}}_3{\mathrm{\α}}_3}{{\mathrm{\ρ}}_1{\mathrm{\α}}_1+{\mathrm{\ρ}}_3{\mathrm{\α}}_3}$

$T_{PP}=\frac{2{\mathrm{\ρ}}_3{\mathrm{\α}}_3}{{\mathrm{\ρ}}_1{\mathrm{\α}}_1+{\mathrm{\ρ}}_3{\mathrm{\α}}_3}---(2-2-11)$

Its result is clearly correct.

Fig. 5 is the situation that the reflection coefficient under the single interface conditions utilizing zoppritz equation to calculate changes with angle, and in figure, horizontal ordinate is angle, and scope is from 0 degree to 60 degree, and ordinate is reflection coefficient, whole change of successively decreasing.In order to verify correctness of the present invention, we use computing method of the present invention, thin intermediate thickness is taken as 0 meter, then model regression of the present invention is single INTERFACE MODEL, program of the present invention is utilized to carry out computational reflect coefficient with angle situation of change, through comparing, method of the present invention and zoeppritz equation are identical for single INTERFACE MODEL, and this also demonstrates correctness of the present invention.

More than show and describe ultimate principle of the present invention and principal character and advantage of the present invention.The technician of the industry should understand; the present invention is not restricted to the described embodiments; what describe in above-described embodiment and instructions just illustrates principle of the present invention; without departing from the spirit and scope of the present invention; the present invention also has various changes and modifications, and these changes and improvements all fall in the claimed scope of the invention.Application claims protection domain is defined by appending claims and equivalent thereof.

Claims (5)

1. computing method for Elastic Multi-layer medium midplane wave reflection coefficient, is characterized in that, specifically comprise following step:

(1) plane harmonization P is established _efrom medium n+1 to object series of strata with incidence, then produce reflected P-wave and reflection wave in n+1 medium, and its reflection coefficient is respectively P _rand S _r, in medium 1, produce transmitted P-wave and transmitted shear wave, its transmission coefficient is respectively P _tand S _t, each layer medium all can be write as scalar potential and vector potential form;

(2) under two-dimensional case, the relational expression between displacement, stress and displacement position is determined;

(3) relating to scalar potential in step (1) and vector potential is brought in step (2), can obtain in (n+1)th layer and the displacement of the 1st layer and stress relation formula;

(4) according in (n+1)th layer that obtains in (3) and the displacement of the 1st layer and stress relation formula, the displacement, the Stress transmit matrix that comprise each layer elasticity coefficient can be obtained, then it is solved, reflection coefficient, transmission coefficient can be obtained.

2. the computing method of Elastic Multi-layer medium midplane wave reflection coefficient according to claim 1, it is characterized in that, step (1) specifically comprises the following steps:

If interlayer comprises n-1 flat seam, sequence is numbered 2,3 ... n, bottom medium is 1 layer; In each layer, compressional wave and transverse wave speed use the α of subscripting respectively, and β represents, under be designated as sequence number, get x coordinate axis and n-th layer bottom boundary coincides; If a ripple incides interlayer end face from n+1 layer direction, note compressional wave incident wave is P _e, now, in n+1 layer, have a longitudinal wave reflection ripple P _r, transverse wave reflection ripple S _r, in 1 layer, have compressional wave transmitted wave P _twith shear wave transmitted wave S _tif compressional wave and shear wave are to each layer incident angle

Then total in n+1 layer medium bit shift is:

Wherein:

σ＝ω/c， k ⁿ⁺¹＝ω/α _n+1，K ⁿ⁺¹＝ω/β _n+1

Wherein, ω is frequency, φ _efor incident longitudinal wave bit function, φ _rreflected P-wave bit function, be amplitude, the e of (n+1) layer incident longitudinal wave be constant, j is imaginary unit, z is direction, for the amplitude of reflected P-wave, x be direction, t is time, ψ _rrepresent the bit function of reflection wave, for the amplitude of reflection wave;

C is the apparent velocity of ripple along interphase direction, according to Snell's law, is all equal to each ripple on interphase;

In n layer medium, compressional wave and shear wave bit function expression formula are respectively:

wherein, for incident longitudinal wave amplitude, for reflected P-wave amplitude, ψ _efor incident shear wave bit function, for incident shear wave amplitude, for reflection wave amplitude, k ⁿthe horizontal wave number of n-th layer compressional wave, K ⁿthe horizontal wave number of n-th layer shear wave;

Then in 1 layer of medium, transmitted wave bit function expression formula is respectively:

Wherein,

φ _tfor transmitted P-wave bit function, for transmitted P-wave amplitude, ψ _tfor transmitted shear wave bit function, for transmitted shear wave amplitude, k ¹for the horizontal wave number of ground floor compressional wave, K ¹for the horizontal wave number of ground floor shear wave.

3. the computing method of Elastic Multi-layer medium midplane wave reflection coefficient according to claim 2, it is characterized in that, in step (2), in two-dimensional case, the relational expression of displacement, stress and displacement position is:

Wherein, u represents displacement, σ _zzrepresent principal strain, the τ along z-axis _zxrepresent shear stress, z is coordinate direction, λ and μ Lame's constant.

4. the computing method of Elastic Multi-layer medium midplane wave reflection coefficient according to claim 3, it is characterized in that, step (3) specifically comprises the following steps:

If n layer thickness is h, calculates the displacement component in n layer and the components of stress, get the value of its z=h, be the displacement on n-th layer end face and components of stress value, be designated as u ⁽ⁿ⁾, ω ⁽ⁿ⁾, bring formula (2-1-2), (2-1-3), (2-1-4) into (2-1-6) formula, can matrix be obtained: wherein p=d ⁿh, Q=s ⁿh, d ⁿfor n-th layer compressional wave vertical wavenumber, h are thickness, s ⁿfor n-th layer shear wave vertical wavenumber;

Wherein

Wherein, s is the 1st layer of shear wave vertical wavenumber, and d is the 1st layer of compressional wave vertical wavenumber, λ _nand μ _nfor the Lame's constant of n-th layer;

Fetch bit moves and the value of the components of stress at z=0 place, can obtain the displacement component on n layer bottom surface and the components of stress, and according to interphase displacement and the stress condition of continuity, it should be equal with the analog value of (n-1) layer end face, is designated as u ^(n-1), ω ^(n-1), as z=0, p=0, Q=0, can be obtained by matrix (2-1-8)

(n-1) displacement on layer end face and components of stress value can be expressed as:

Can set up the relation between the displacement component of (n) layer and (n-1) layer end face and the components of stress by above-mentioned formula, for this reason, solving equation group (2-1-10) can have:

Wherein (b _ij) represent inverse matrix:

Wherein, K is the horizontal wave number of shear wave

Formula (2-1-11) is substituted into (2-1-7) can obtain

Get mark (a _ij) representing matrix product (B _ij) (b _ij), be similarly 4 × 4 rank square formations, each element is:

a ₁₁＝2sin ²γcosp-cos2γcosQ

a ₁₂＝j(tanθcos2γsinp+sin2γsinQ)

a ₂₂＝cos2γcosp+2sin ²γcosQ

a ₃₁＝2jωρβsinγ(cosp-cosQ)cos2γ

a ₃₃＝cos2γcosp+2sin ²γcosQ

a ₃₄＝2jρβ ²(cos2γtanθsinp-sin2γsinQ)

a ₄₄＝2sin ²γcosp+cos2γcosQ (2-1-13)

Wherein, γ=i _s, θ=i _drepresent the incident angle of ripple in layer respectively, so there is sin γ=σ/K, sin θ=σ/k, α, β is p-and s-wave velocity, ρ is Media density, establishes (n) layer and the relation between (n-1) layer end face displacement component and the components of stress in formula (2-1-12), calculates (a _ij) matrix element, will use the parameter values of (n) layer at formula (2-1-13), for this reason, in formula (2-1-12), matrix of coefficients corner brace (n) represents to have:

Obtain the recursion formula of displacement component on each interphase of interlayer and the components of stress, meet relational expression (2-1-14) and be equivalent to and meet displacement on interlayer interphase and the stress condition of continuity.

5. the computing method of Elastic Multi-layer medium midplane wave reflection coefficient according to claim 4, it is characterized in that, step (4) specifically comprises the following steps:

In order to the transmission coefficient determining interlayer end face reflection coefficient and propagate below interlayer bottom surface through interlayer, set up the relation of displacement component on (n) layer end face and (1) layer end face and the components of stress according to formula (2-1-14):

(2-1-15) is out of shape

Formula (2-1-1), formula (2-1-5) are substituted into formula (2-1-16) and obtain matrix equation (2-1-17), reflection coefficient and the transmission coefficient of interlayer after solving, can be obtained;

Consider A _ijbe plural number, make A ₁₁=R ₁₁+ iI ₁₁, A ₁₂=R ₁₂+ iI ₁₂, A ₁₃=R ₁₃+ iI ₁₃,

Wherein, ρ _n+1the density of (n+1)th layer, the horizontal wave number of shear wave of (n+1)th layer

Wherein B is called dematrix, and its each element is:

B ₁₁＝B ₂₂＝-iσ，B ₁₂＝is _n+1

B ₂₁＝id _n+1，

Wherein, D _n1for plural number;

A _n1imaginary part, D _n2for plural A _n2imaginary part, d ₁for thickness, C _n3for plural A _n3imaginary part, ρ ₁for ground floor density, represent horizontal wave number, the C of ground floor shear wave _n4for plural A _n4imaginary part, s ₁ground floor shear wave vertical wavenumber, C _n1for plural A _n1imaginary part, C _n2for plural A _n2imaginary part

Solve dematrix expression formula formula (2-2-17) and the longitudinal wave reflection coefficients R of plane harmonization at elastic layer system can be obtained _pp, transverse wave reflection coefficient V ^/, shear wave transmission coefficient W ^/, compressional wave transmission coefficient W.