CN104950332A - 一种弹性多层介质中平面波反射系数的计算方法 - Google Patents

Abstract

本发明公开了一种弹性多层介质中平面波反射系数的计算方以下法，步骤如下(1)设平面谐和波从介质n+1向目的层系入射，则在n+1介质内产生反射纵波及横波，在介质1内产生透射纵波及横波，各层介质均可写成标量位与向量位形式；(2)确定位移、应力与位移位存在的关系；(3)把步骤(1)中的标量位和向量位带入到步骤(2)中获得第n+1层内及第1层的位移和应力关系式；(4)根据(3)中获得的第n+1层内及第1层的位移和应力关系式得到包含各层弹性系数的位移、应力传递矩阵，对其进行求解，可获得反射及透射系数。本发明充分考虑了弹性层系多次波和转换波及厚度、频率对反射系数的影响，适宜于薄层AVO分析模拟，计算效率高。

Description

一种弹性多层介质中平面波反射系数的计算方法

技术领域

本发明涉及一种弹性多层介质中平面波反射系数的计算方法，属于薄层AVO模拟技术领域。

背景技术

模拟薄层AVO分为两大类，(1)时间域褶积模型，对某一特定角度情况下，对时间域采样点逐个计算反射系数，利用子波与反射系数进行褶积，形成该角度对应地震道，该方法虽然计算快捷，公式简单，容易实现，但是难以体现薄层效应，模拟精度较低。(2)采用弹性波方程的积分解有限差分法模拟薄层波场，可模拟出薄层效应，但是计算效率较低，受到建模影响较大，工业生产难以实用化。

常规的Zoeppritz方程计算反射系数都是基于单界面的，无厚度变量，虽然也能计算AVO特征，但不能直接分析厚度对各频率分量的影响效果。

发明内容

针对现有技术存在的不足，本发明目的是提供一种弹性多层介质中平面波反射系数的计算方法，充分考虑了弹性层系多次波和转换波及厚度、频率对反射系数的影响，适宜于薄层AVO分析模拟，计算效率高，模拟精度高。

为了实现上述目的，本发明是通过如下的技术方案来实现：

本发明的一种弹性多层介质中平面波反射系数的计算方法，具体包括以下几个步骤：

(1)设平面谐和波P_e从介质n+1向目的层系以入射，则在n+1介质内产生反射纵波和反射横波，其反射系数分别为P_r和S_r，在介质1内产生透射纵波和透射横波，其透射系数分别为P_t和S_t，各层介质均可写成标量位与向量位形式；

(2)在二维情况下，确定位移、应力与位移位之间的关系式；

(3)把步骤(1)中涉及到标量位和向量位带入到步骤(2)中，可以获得第n+1层内及第1层的位移和应力关系式；

(4)根据(3)中获得的第n+1层内及第1层的位移和应力关系式，可以得到包含各层弹性系数的位移、应力传递矩阵，然后对其进行求解，即可获得反射系数、透射系数。

步骤(1)具体包括以下步骤：

设夹层包括n－1个水平层，层序编号为2、3、……n，下部介质为1层；各层中纵波和横波波速分别用带下标的α,β表示，下标为层序号，取x坐标轴与第n层底界面相重合；设一个波自n+1层方向入射到夹层顶面，记纵波入射波为P_e，此时，在n+1层中有一纵波反射波P_r，横波反射波S_r，在1层中有纵波透射波P_t和横波透射波S_t，设纵波和横波到各层入射角

则在n+1层介质中总的位移位为：

${\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)$

其中： $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

其中，ω为频率，φ_e为入射纵波位函数、φ_r反射纵波位函数、为第(n+1)层入射纵波的振幅、e为常数、j为虚数单位、z为方向、为反射纵波的振幅、x为方向、t为时间、ψ_r代表反射横波的位函数、为反射横波的振幅；

C为波沿分界面方向的视速度，根据斯奈尔定律，对分界面上的各个波都是相等的；

在n层介质中纵波和横波位函数表达式分别为：

${\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)$

其中，为入射纵波振幅、为反射纵波振幅、ψ_e为入射横波的位函数、为入射横波振幅、为反射横波振幅、kⁿ第n层纵波的水平波数、Kⁿ第n层横波的水平波数；

则在1层介质中透射波位函数表达式分别为：

${\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)$

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

φ_t为透射纵波位函数、为透射纵波振幅、ψ_t为透射横波位函数、为透射横波振幅、k¹为第一层纵波水平波数、K¹为第一层横波水平波数。

步骤(2)中，在二维情况，位移、应力与位移位的关系式为：

$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)$

其中，u代表位移、σ_zz代表沿z轴的主应变、τ_zx代表剪应力、z为坐标方向、λ和μ拉梅常数。

步骤(3)具体包括以下步骤：

设n层厚度为h，计算n层中的位移分量和应力分量，取其z＝h的值，为第n层顶面上的位移和应力分量值，记为u⁽ⁿ⁾、ω⁽ⁿ⁾、将式(2-1-2)、(2-1-3)、(2-1-4)带入(2-1-6)式，可以得到矩阵：其中p＝dⁿh，Q＝sⁿh，dⁿ为第n层纵波垂直波数、h为厚度、sⁿ为第n层横波垂直波数；

$\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)$

其中

$\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)$

其中，s为第1层横波垂直波数和d为第1层纵波垂直波数、λ_n和μ_n为第n层的拉梅常数 

取位移和应力分量在z＝0处的值，可得n层底面上的位移分量和应力分量，根据分界面位移和应力连续条件，它应与(n－1)层顶面的相应值相等，记为u^(n-1)、ω^(n-1)、 当z＝0时，p＝0，Q＝0，由矩阵(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)层顶面上的位移与应力分量值可表示为：

$\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)$

由上述公式可以建立起(n)层与(n－1)层顶面的位移分量和应力分量之间的关系，为此，求解方程组(2-1-10)可有：

$\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)$

其中(b_ij)表示的逆矩阵：

$\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)$

其中，K为横波水平波数

将式(2-1-11)代入(2-1-7)可得

$\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)$

取记号(a_ij)表示矩阵乘积(B_ij)(b_ij)，同样为4×4阶方阵，各元素为：

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) 

其中，γ＝i_s、θ＝i_d分别表示波在层中的入射角，所以有sinγ＝σ/K，sinθ＝σ/k，α,β为纵横波速度，ρ为介质密度，在式(2-1-12)中建立了(n)层和(n－1)层顶面位移分量和应力分量之间的关系，计算(a_ij)矩阵元素，在公式(2-1-13)将使用(n)层的参 数数值，为此，在式(2-1-12)中系数矩阵用上角码(n)表示，可有：

$\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)$

得到了夹层各分界面上的位移分量和应力分量的递推公式，满足关系式(2-1-14)等价于满足夹层分界面上的位移与应力连续条件。

步骤(4)具体包括以下步骤：

为了确定夹层顶面反射系数和透过夹层在夹层底面以下传播的透射系数，根据公式(2-1-14)建立(n)层顶面和(1)层顶面上的位移分量和应力分量的关系：

$\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)变形得

$\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)$

将式(2-1-1)、式(2-1-5)代入式(2-1-16)得到矩阵方程(2-1-17)，求解后可以获得夹层的反射系数和透射系数；

考虑到A_ij是复数，令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)$

其中，ρ_n+1第n+1层的密度、第n+1层的横波水平波数

其中B称为解矩阵，它的各元素为：

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}$

其中，D_n1为复数；

A_n1的虚部、D_n2为复数A_n2的虚部、d₁为厚度、C_n3为复数A_n3的虚部、ρ₁为第一层的密度、表示第一层横波的水平波数、C_n4为复数A_n4的虚部、s₁第一层横波垂直波数、C_n1为复数A_n1的虚部、C_n2为复数A_n2的虚部解出解矩阵表达式公式(2-2-17)即可得到平面谐和波在弹性层系的纵波反射系数R_pp、横波反射系数V^/、横波透射系数W^/、纵波透射系数W。

本发明基于弹性层系的反射系数谱理论，应用弹性层系的位移和应力递推公式，针对固体多层介质，推导出了纵、横波反射、透射系数公式，该方法可以计算出不同入射角度的平面波反射系数，在石油物探上，可以用弹性层系模型代表地质上的薄(互)层地层结构，因此本发明的方法可以用于计算薄(互)层AVO(Amplitude versus offset)，可以解决多次波、转换波对薄层AVO的影响，同时可以定量模拟不同储层厚度、不同频率变化引起的AVO变化，同时可以用于解释地震道集资料与常规商业化软件基于井zoeppritz方程正演道集不匹配的问题；本发明可以做为目前石油勘探工业商业化软件的一个重要补充。

附图说明

图1为层状介质的反射和透射模型；

图2为单个薄层的反射和透射；

图3为不同储层厚度情况下的储层AVO特征模拟； 

图4为不同入射波频率情况下的储层AVO特征模拟； 

图5为单层介质与传统商业化软件相符。

具体实施方式

为使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施方式，进一步阐述本发明。

常规的Zoeppritz方程计算反射系数都是基于单界面的，无厚度变量，虽然也能计算AVO特征，但不能直接分析厚度对各频率分量的影响效果。弹性多层介质的平面波反射系数充分考虑了弹性层系多次波和转换波及厚度、频率对反射系数的影响，因此更适宜于研究薄层、薄互层的AVO分析及正演模拟。因此本发明详细介绍了弹性多层介质中的平面波反射系数及正演计算方法，并对单层模型进行了推导验证。

弹性多层介质中的平面波反射系数计算方法如下：

在两个半无限弹性介质中间的一个夹层。该夹层由多个水平层次构成。设有一个平面波自上部介质入射夹层顶面，则产生一个夹层反射波在上部介质中传播，同时波亦将透过夹层， 产生一个在下部介质中传播的透射波。计算夹层的反射系数和透射系数。

如图1所示，夹层由n－1个水平层构成，层序编号为2、3、……n；下部介质为1层。各层中纵波和横波波速分别用带下标的α,β表示，下标为层序号。取x坐标轴与第n层底界面相重合。有一个P－SV波***自n+1层方向入射到夹层顶面，记纵波入射波为P_e，此时，在n+1层中有一纵波反射波P_r，横波反射波S_r；在1层中有纵波透射波P_t和横波透射波S_t。设纵波和横波到各层入射角

则在n+1层介质中总的位移位为：

${\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)$

其中： $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为波沿分界面方向的视速度，根据斯奈尔定律，对分界面上的各个波都是相等的。

在n层介质中纵波和横波位函数表达式分别为：

${\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)$

则在1层介质中透射波位函数表达式分别为：

${\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)$

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

在二维情况，位移、应力与位移位的关系式为：

$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)$

设n层厚度为h，计算n层中的位移分量和应力分量，取其z＝h的值，为第n层顶面上的位移和应力分量值，记为u⁽ⁿ⁾、ω⁽ⁿ⁾、将式(2-1-2)、(2-1-3)、(2-1-4)带入(2-1-6)式，可以得到矩阵：其中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)$

其中

$\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)$

令一方面，取位移和应力分量在在z＝0处的值，可得n层底面上的位移分量和应力分量。根据分界面位移和应力连续条件，它应与(n－1)层顶面的相应值相等，记为u^(n-1)、ω^(n-1)、 当z＝0时，p＝0，Q＝0，由矩阵(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)层顶面上的位移与应力分量值可表示为：

$\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)$

由上述公式可以建立起(n)层与(n－1)层顶面的位移分量和应力分量之间的关系。为此，求解方程组(2-1-10)可有：

$\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)$

其中(b_ij)表示的逆矩阵：

$\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)$

将式(2-1-11)代入(2-1-7)可得

$\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)$

取记号(a_ij)表示矩阵乘积(B_ij)(b_ij)，同样为4×4阶方阵，各元素为：

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) 

其中γ＝i_s、θ＝i_d分别表示波在层中的入射角，所以有sinγ＝σ/K，sinθ＝σ/k。其它参数，如α,β为纵横波速度，ρ为介质密度，ω为频率。在式(2-1-12)中建立了(n)层和(n－1)层顶面位移分量和应力分量之间的关系，计算(a_ij)矩阵元素，在公式(2-1-13)将使用(n)层的参数数值。为此，在式(2-1-12)中系数矩阵用上角码(n)表示，可有：

$\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)$

得到了夹层各分界面上的位移分量和应力分量的递推公式。满足关系式(2-1-14)等价于满足夹层分界面上的位移与应力连续条件。

为了确定夹层顶面反射系数和透过夹层在夹层底面以下传播的透射系数，根据公式(2-1-14)建立(n)层顶面和(1)层顶面上的位移分量和应力分量的关系：

$\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)变形得

$\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)$

将式(2-1-1)、式(2-1-5)代入式(2-1-16)得到矩阵方程(2-1-17)，求解后可以获得夹层的反射系数和透射系数。

考虑到A_ij是复数，令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)$

其中B称为解矩阵，它的各元素为：

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}$

解出解矩阵表达式公式(2-2-17)即可得到平面谐和波在弹性层系的纵波反射系数R_pp、横波反射系数V^/、横波透射系数W^/、纵波透射系数W。

图3为单个薄层模型所对应的纵波反射系数随角度变化的成果，图中，横坐标为角度，从0度递增到60度，纵坐标为纵波反射系数，薄层厚度5米到30米之间变化，从图中可以看出不同储层厚度对应的薄层AVO形态相似，但四条曲线的截距不同，四条曲线的梯度也不同，这要求我们判断薄层的AVO时要考虑薄层的厚度变化。

图4为单个薄层模型所对应的纵波反射系数随角度变化的成果，图中，横坐标为角度，从0度递增到45度，纵坐标为纵波反射系数，频率由10hz到90hz之间变化，从图中可以看出不同频率对应的薄层AVO形态相似，10条曲线都是单调递增的，但10条曲线的梯度各不相同，这要求我们判断薄层的AVO时要考虑入射波的主频变化，尤其同一套薄储层在埋深变化较大的情况下。

为检验上述理论的正确性，按上面的理论推导垂直入射单个薄层情况下的反射系数。

设两个半无限介质中间有一个厚度为h的固体夹层，计算该夹层的反射系数和透射系数。如图2所示，在1和3两个半无限介质中间夹有一个薄层，其厚度为h。有一纵波自介质3方向垂直入射到薄层顶面，将产生一个在介质3中传播的反射波p_r和在介质1中传播的透 射波p_t，根据图中所示的坐标系，各波表达式为：

${\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)$

对第3介质有位函数：

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

对第1介质有位函数：

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

它们应满足的边界条件，根据(2-1-15)，为

$\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)$

其中考虑了垂直入射的情况下，μ＝0，τ_zx＝0。将式(2-2-2)，式(2-2-3)带入式(2-2-4)式，在z＝0上有：

${\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上可有：

${\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}$

边界条件方程组为：

$\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)$

方程组的解就是所求的反射系数和透射系数，它们是：

$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)$

其中T_PP透射系数取透射波与入射波的位移振幅比。将薄层介质2参数及入射角数据，θ＝0，γ＝0，带入式(2-1-13)，可有：

$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)$

将式(2-2-8)代入式(2-2-6)，式(2-2-7)，可得

$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)$

当h＝0时，可以得到一个分界面时纵波垂直入射时的反射系数和透射系数：

$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)$

其结果很明显是正确的。

图5是利用zoppritz方程计算的单界面情况下的反射系数随角度变化的情况，图中，横坐标为角度，范围从0度到60度，纵坐标为反射系数，整个递减变化。为了验证本发明正确性，我们用本发明的计算方法，把中间薄层厚度取为0米，则本发明的模型退变为单界面模型，利用本发明的程序进行计算反射系数随角度变化情况，经过比较，本发明的方法与zoeppritz方程对于单界面模型两者是相同的，这也验证了本发明的正确性。

以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (5)

1.一种弹性多层介质中平面波反射系数的计算方法，其特征在于，具体包括以下几个步骤：

(1)设平面谐和波P_e从介质n+1向目的层系以入射，则在n+1介质内产生反射纵波和反射横波，其反射系数分别为P_r和S_r，在介质1内产生透射纵波和透射横波，其透射系数分别为P_t和S_t，各层介质均可写成标量位与向量位形式；

(2)在二维情况下，确定位移、应力与位移位之间的关系式；

(3)把步骤(1)中涉及到标量位和向量位带入到步骤(2)中，可以获得第n+1层内及第1层的位移和应力关系式；

(4)根据(3)中获得的第n+1层内及第1层的位移和应力关系式，可以得到包含各层弹性系数的位移、应力传递矩阵，然后对其进行求解，即可获得反射系数、透射系数。

2.根据权利要求1所述的弹性多层介质中平面波反射系数的计算方法，其特征在于，步骤(1)具体包括以下步骤：

设夹层包括n－1个水平层，层序编号为2、3、……n，下部介质为1层；各层中纵波和横波波速分别用带下标的α,β表示，下标为层序号，取x坐标轴与第n层底界面相重合；设一个波自n+1层方向入射到夹层顶面，记纵波入射波为P_e，此时，在n+1层中有一纵波反射波P_r，横波反射波S_r，在1层中有纵波透射波P_t和横波透射波S_t，设纵波和横波到各层入射角

则在n+1层介质中总的位移位为：

其中：

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

其中，ω为频率，φ_e为入射纵波位函数、φ_r反射纵波位函数、为第(n+1)层入射纵波的振幅、e为常数、j为虚数单位、z为方向、为反射纵波的振幅、x为方向、t为时间、ψ_r代表反射横波的位函数、为反射横波的振幅；

C为波沿分界面方向的视速度，根据斯奈尔定律，对分界面上的各个波都是相等的；

在n层介质中纵波和横波位函数表达式分别为：

其中，为入射纵波振幅、为反射纵波振幅、ψ_e为入射横波的位函数、为入射横波振幅、为反射横波振幅、kⁿ第n层纵波的水平波数、Kⁿ第n层横波的水平波数；

则在1层介质中透射波位函数表达式分别为：

其中，

φ_t为透射纵波位函数、为透射纵波振幅、ψ_t为透射横波位函数、为透射横波振幅、k¹为第一层纵波水平波数、K¹为第一层横波水平波数。

3.根据权利要求2所述的弹性多层介质中平面波反射系数的计算方法，其特征在于，步骤(2)中，在二维情况，位移、应力与位移位的关系式为：

其中，u代表位移、σ_zz代表沿z轴的主应变、τ_zx代表剪应力、z为坐标方向、λ和μ拉梅常数。

4.根据权利要求3所述的弹性多层介质中平面波反射系数的计算方法，其特征在于，步骤(3)具体包括以下步骤：

设n层厚度为h，计算n层中的位移分量和应力分量，取其z＝h的值，为第n层顶面上的位移和应力分量值，记为u⁽ⁿ⁾、ω⁽ⁿ⁾、将式(2-1-2)、(2-1-3)、(2-1-4) 带入(2-1-6)式，可以得到矩阵：其中p＝dⁿh，Q＝sⁿh，dⁿ为第n层纵波垂直波数、h为厚度、sⁿ为第n层横波垂直波数；

其中

其中，s为第1层横波垂直波数，d为第1层纵波垂直波数，λ_n和μ_n为第n层的拉梅常数；

取位移和应力分量在z＝0处的值，可得n层底面上的位移分量和应力分量，根据分界面位移和应力连续条件，它应与(n－1)层顶面的相应值相等，记为u^(n-1)、ω^(n-1)、 当z＝0时，p＝0，Q＝0，由矩阵(2-1-8)可以得到

(n－1)层顶面上的位移与应力分量值可表示为：

由上述公式可以建立起(n)层与(n－1)层顶面的位移分量和应力分量之间的关系，为此，求解方程组(2-1-10)可有：

其中(b_ij)表示的逆矩阵：

其中，K为横波水平波数

将式(2-1-11)代入(2-1-7)可得

取记号(a_ij)表示矩阵乘积(B_ij)(b_ij)，同样为4×4阶方阵，各元素为：

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) 

其中，γ＝i_s、θ＝i_d分别表示波在层中的入射角，所以有sinγ＝σ/K，sinθ＝σ/k，α,β为纵横波速度，ρ为介质密度，在式(2-1-12)中建立了(n)层和(n－1)层顶面位移分量和应力分量之间的关系，计算(a_ij)矩阵元素，在公式(2-1-13)将使用(n)层的参数数值，为此，在式(2-1-12)中系数矩阵用上角码(n)表示，可有：

得到了夹层各分界面上的位移分量和应力分量的递推公式，满足关系式(2-1-14)等价于满足夹层分界面上的位移与应力连续条件。

5.根据权利要求4所述的弹性多层介质中平面波反射系数的计算方法，其特征在于，步骤(4)具体包括以下步骤：

为了确定夹层顶面反射系数和透过夹层在夹层底面以下传播的透射系数，根据公式(2-1-14)建立(n)层顶面和(1)层顶面上的位移分量和应力分量的关系：

对(2-1-15)变形得

将式(2-1-1)、式(2-1-5)代入式(2-1-16)得到矩阵方程(2-1-17)，求解后可以获得夹层的反射系数和透射系数；

考虑到A_ij是复数，令A₁₁＝R₁₁+iI₁₁，A₁₂＝R₁₂+iI₁₂，A₁₃＝R₁₃+iI₁₃，…

其中，ρ_n+1第n+1层的密度、第n+1层的横波水平波数

其中B称为解矩阵，它的各元素为：

B₁₁＝B₂₂＝-iσ，B₁₂＝is_n+1

B₂₁＝id_n+1，

其中，D_n1为复数；

A_n1的虚部、D_n2为复数A_n2的虚部、d₁为厚度、C_n3为复数A_n3的虚部、ρ₁为第一层的密度、表示第一层横波的水平波数、C_n4为复数A_n4的虚部、s₁第一层横波垂直波数、C_n1为复数A_n1的虚部、C_n2为复数A_n2的虚部

解出解矩阵表达式公式(2-2-17)即可得到平面谐和波在弹性层系的纵波反射系数R_pp、横波反射系数V^/、横波透射系数W^/、纵波透射系数W。