基于散射矩阵分解的反射系数二阶近似

龚诚诚 吴国忱*② 单俊臻

(①中国石油大学(华东)地球科学与技术学院，山东青岛 266580；②海洋国家实验室海洋矿产资源评价与探测技术功能实验室，山东青岛 266580)

0 引言

AVO技术利用振幅信息研究岩性、检测油气，经过数十年的发展，已经形成了一套完整而严谨的理论体系。Zoeppritz[1]研究了平面波入射在水平界面上的纵横波反射、透射情况，建立了Zoeppritz方程。如何利用地震反射振幅随入射角变化的信息检测油气成为地球物理学家关注的问题，但是由于Zoeppritz方程的结构过于复杂且物理含义并不直观，不易进行数值计算，一直没有得到很好的应用。Muskat等[2]首先提出了AVO概念，即“平面波的反射和透射系数是入射角的函数”；Koefoed[3]将泊松比与反射系数相联系，给出了精确Zoeppritz方程的近似公式，使AVO技术得到快速发展；Bortfeld[4]在界面两侧岩性参数变化小的假设下得到了区分流体、固体的近似公式；Richards等[5]在假设介质参数变化小的情况下根据位移和应力连续条件推导出较直观的反射/透射系数近似表达式；Aki等[6]给出了反射/透射系数精确表达式，并利用射线参数与角度的关系给出了弱参数变化下的反射/透射系数近似式。Ostrander[7]利用AVO技术识别“亮点”型含油气砂岩，使AVO 技术用于实际生产。在前人的基础上，人们从不同角度推导出各类近似表达式，表征纵波反射系数与不同岩石物理参数间的关系。如Shuey[8]给出了具有不同角度项的关于泊松比的纵波反射系数近似式；Smith 等[9]利用Gardner公式，采用加权分析提出了速度随密度变化的经验公式；Hilterman[10]在 Shuey近似的基础上，给出了关于泊松比的另一种近似式；郑晓东[11]、杨绍国等[12]推出了Zoeppritz方程的近似幂级数表达式；Fatti等[13]发布了关于相对波阻抗的纵波反射系数近似式；Goodway等[14]提出了关于拉梅参数的纵波反射系数近似式；Xu等[15]引入剪切模量、体积模量和压缩模量简化了纵波反射系数；Gray等[16]给出了关于剪切模量、体积模量和压缩模量相对变化量的近似表达式；Russell等[17]在多孔流体饱和岩石的前提下提出了关于Russell流体因子的纵波反射系数近似表达式；宗兆云等[18]推出了关于杨氏模量和泊松比纵波反射系数近似方程。

常规AVO反演利用精确Zoeppritz方程的一阶近似式，更加适用于弱介质变化的反射界面、小角度或小炮检距反射问题，其假设条件导致计算误差[19]，不利于准确提取密度参数[20]，不能很好地预测复杂储层或中深部储层，且无法充分利用近临界角数据。因此，针对复杂储层或中深部储层，需要改进反射系数近似方程，使其适应近临界角入射，充分利用大角度地震数据提高反演精度。Xu等[21]基于Aki-Richards近似方法给出了具有归一化速度、密度二阶项的近似公式；Wang[22]对纵波反射系数精确方程的射线参数表达形式进行泰勒展开，得到了关于垂直慢度、射线参数的二阶近似公式；Ramos[23]对纵波反射系数精确方程进行三阶泰勒展开，给出了适用于强反射界面的横波反射系数近似公式；Ursin等[24]利用上、下行波特征向量矩阵及矩阵对称性对反射/透射系数进行泰勒二阶展开，得到二阶近似公式；Charles[25]提出最优Zoeppritz近似，通过研究纵波反射系数精确值与纵横波速度和密度相对变化量各阶泰勒展开的关系得到拟线性近似，提高了近临界角处的反射系数反演精度。

基于反射系数精确表达式利用Aki-Richards近似方法进行高阶近似过于繁琐[22-23]，利用特征向量化简物理意义不够明确[24]。为此，本文直接从P-SV平面波入射/散射矩阵出发，给出了一种求取散射矩阵高阶近似的方法，即利用扰动思想将散射矩阵分解为背景矩阵与一阶、二阶扰动矩阵，求取纵波反射系数背景项与一阶、二阶扰动项，推得纵波反射系数的二阶近似公式。模型对比分析表明，所推公式在中高角度乃至近临界角入射情况下具有较高的精度，对密度参数的敏感性更高，为充分利用大炮检距地震数据准确地反演物性参数提供了基础。

1 基于散射矩阵分解的高阶近似

1.1 散射矩阵分解

Aki等[6]研究了P-SV平面波在固体—固体分界面上的反射与透射问题，得到了入射与散射P-SV平面波完整系统(图1)，并基于此系统推导了完整的散射矩阵。

(1)

M=Mb+ΔM1+ΔM2

(2)

N=Nb+ΔN1+ΔN2

(3)

R=Rb+ΔR1+ΔR2

(4)

式中：Mb、Nb、Rb分别为M、N、R的背景矩阵，由各向同性介质参数构成；ΔM1、ΔN1、ΔR1分别为M、N、R的一阶扰动矩阵，由介质参数的一阶变化构成；ΔM2、ΔN2、ΔR2分别为M、N、R的二阶扰动矩阵，由介质参数的二阶变化构成。

图1 入射与散射的P-SV平面波完整系统

1.2 纵波反射系数二阶近似定量表达(附录A)

对式(1)变形得到如下矩阵表达式

MR=N

(5)

首先，仅考虑系数矩阵M、N的背景项Mb、Nb时，可以推得散射矩阵的背景项Rb

(6)

当某种波型入射时，由于背景介质均为各向同性，故背景介质中不存在该波型的反射以及其他波型的反射与透射，在Rb中相同波型的透射系数为1，其他波型的反射、透射系数为0。

假设M、N与R仅由背景矩阵与一阶扰动矩阵构成，则式(5)可表示为

(Mb+ΔM1)(Rb+ΔR1)=(Nb+ΔN1)

(7)

对式(7)进行整理，忽略一阶以上高阶项，得

ΔR1=(Mb)-1(ΔN1-ΔM1Rb)

(8)

假设M、N与R由背景矩阵与一阶扰动矩阵、二阶扰动矩阵共同构成，则式(5)可表示为

(Mb+ΔM1+ΔM2)(Rb+ΔR1+ΔR2)

=Nb+ΔN1+ΔN2

(9)

对式(9)进行整理，忽略二阶以上高阶项，得

ΔR2=(Mb)-1(ΔN2-ΔM1ΔR1-ΔM2Rb)

(10)

整理二阶扰动矩阵及其内部参数，经过计算可以得到ΔN2=ΔM2Rb，则式(10)可以进一步简化为

ΔR2=(Mb)-1(-ΔM1ΔR1)

(11)

(12)

分别计算散射矩阵的背景项、一阶扰动项、二阶扰动项并相加，整理可得入射、散射P-SV平面波完整系统下的纵波反射系数RPP关于纵、横波速度与密度的二阶近似定量表达式

(13)

式中： Δρ、ΔVP、ΔVS分别为界面两侧介质密度、纵波速度、横波速度的差值；ρ、VP、VS分别为界面两侧介质密度、纵波速度、横波速度的平均值；θP、θS分别为界面两侧下、上行纵(θP1、θP2)、横波入射角(θS1、θS2)的平均值；p为射线参数。

1.3 横波因子cosθS

观察式(13)发现，基于散射矩阵得到的纵波反射系数的二阶项受横波入射角、透射角的平均值θS影响，但在实际纵波反演方法中不存在θS。因此，文中将cosθS称为横波因子，并由反射界面的弹性参数近似表示。

由Snell定律可知

(14)

假设界面两侧介质参数变化较小，入射角、反射角与透射角均为实数，则界面两侧介质参数可由p的一阶扰动项及其平均值的形式表示，即

VS1=VS-ΔVS/2

VS2=VS+ΔVS/2

VP1=VP-ΔVP/2

VP2=VP+ΔVP/2

θS1=θS-ΔθS/2

θS2=θS+ΔθS/2

θP1=θP-ΔθP/2

θP2=θP+ΔθP/2

将各参数代入式(14)，得

(15)

进一步化简并忽略高阶扰动项，得

(16)

横波因子可由实测纵、横波速度与入射角定量表示为

(17)

将式(17)代入式(13)，可得到纵波反射系数的另一种表达形式

(18)

式(18)为化简后的关于纵、横波速度以及密度的纵波反射系数二阶近似式。与式(13)相比，用横纵波速度比的形式表示射线参数p与横波因子cosθS，使近似公式的参数形式得到统一，便于选取实际反演参数，更加适用于AVO反演。

2 纵波反射系数二阶近似方程精度分析

2.1 二阶项对精度的影响

若忽略式(18)中的二阶扰动项，便得到关于纵、横波速度以及密度的纵波反射系数一阶近似方程

(19)

(20)

本文借鉴Wang[22]的四种模型，分别利用不同近似表达式给出了不同岩石界面下的纵波反射系数曲线。这四种模型分别表示页岩/砂岩、页岩/石灰岩(或白云岩)、硬石膏/砂岩和硬石膏/石灰岩(或白云岩)界面，模型参数如表1所示。

表1 模型参数

图2为RPP曲线(未考虑横波因子)。由图可见：①当为页岩/砂岩(图2a)、页岩/石灰岩(或白云岩)(图2b)界面时，上层介质的波阻抗小于下层介质(正阻抗界面)，Aki-Richards近似方程仅在小角度入射时得到的RPP具有较高精度，在中高角度处会产生明显误差，而二阶近似方程在中高角度仍能取得很好的近似效果，且本文得到的二阶近似方程(式(18))在近临界角处仍能取得很高的近似效果；②由于页岩/砂岩(图2a)界面两侧密度与纵波速度变化小、横波速度变化大，页岩/石灰岩(或白云岩)(图2b)界面两侧密度与速度变化趋势相反，因此后者的二阶近似方程精度更高，表明二阶RPP反映了横波的影响，即变化小的横波速度提高了二阶RPP的精度；③当为硬石膏/砂岩(图2c)、硬石膏/石灰岩(或白云岩)(图2d)界面时，上层介质波阻抗大于下层介质(负阻抗界面)，在小角度入射时三种近似方程得到的RPP均具很高的精度，但在中高角度处仅有二阶近似方程可以满足精度要求，且本文推导的二阶近似方程(式(18))所产生的误差最小；④由于硬石膏中的密度、纵波速度、横波速度远大于砂岩与石灰岩，介质参数变化大，进行射线参数与横波因子处理时，在大角度情况下会产生较大误差，二阶近似方程的精度随着入射角的增大而逐渐变低。

图2 RPP曲线(未考虑横波因子)

2.2 横波因子对精度的影响

将横波因子表达式(式(17))代入式(13)与式(20)，可得到RPP曲线(考虑横波因子)(图3)。对比未考虑横波因子(图2)与考虑横波因子(图3)的RPP曲线可知，后者的精度较高，确保了近似方程在入射角增大情况下的适用性，如本文提出的二阶近似方程(式(18))曲线与精确方程曲线在近临界角处几乎完全拟合。当界面两侧物性差异较大时(图3c、图3d)，考虑横波因子明显提高了近似方程的精度。

图3 RPP曲线(考虑横波因子)

3 密度敏感性分析

基于AVO反演可以从地震数据中获取相应介质参数，小炮检距(小入射角)地震振幅数据主要与地震波速度相关，密度参数对中高炮检距(中高入射角)地震振幅数据更敏感[26-29]，但是由于纵波反射系数一阶近似方程在中高入射角时无法获得很高的近似效果，因此基于此类方程的AVO反演方法的密度反演精度相对有限[30]。Fatti等[13]提出了相对波阻抗的概念，并在Aki-Richards的反射系数近似式的基础上推导了关于纵波阻抗、横波阻抗、密度的Fatti近似式。由于纵、横波阻抗及纵、横波速度都与密度直接相关，其中密度对阻抗的影响最大[31]。因此为了更清晰地体现本文所提近似方程(式(18))在中高角度地震数据密度反演中的优势，将式(18)

改写为关于纵波阻抗、横波阻抗、密度的近似方程

(21)

图4 Fatti近似公式RPP曲线

4 结束语

本文基于P-SV平面波入射与散射完整系统下的散射矩阵，给出了一种利用矩阵分解推导反射/透射系数高阶近似的方法，将散射矩阵分解为背景矩阵与不同阶数扰动矩阵，得到了基于入射角的纵波反射系数二阶近似公式。与常规AVO分析中的Aki-Richards线性近似[6]、Wang[22]的二阶近似表达式相比，文中所提反射系数二阶近似公式在中高角度入射时与精确方程具有良好的吻合性，且在近临界角入射情况下也具有较高的精度，当界面两侧横波速度与密度变化差异适中时，近似效果可以达到最佳。因此利用本文的近似公式可以充分利用大炮检距地震振幅数据，解决了一阶近似公式未能充分利用临界角数据的问题，提高了中高入射角的反射系数反演精度，并且为近临界角反演提供了一种方法。利用AVO反演可以提取各种弹性、物性参数，其中密度参数对小入射角地震响应的影响弱，在大炮检距地震振幅数据中富含更精细的信息，常规地震方法很难获得稳定的密度参数估计。本文的近似公式增强了对密度的敏感性，提高了富含密度信息的大炮检距地震振幅数据的利用程度，为准确而稳定地密度反演提供了基础。