正交各向异性介质反射系数精确解

张雪莹 孙鹏远 马学军 芦 俊 李梦琦

（①中国地质大学（北京）能源学院，北京100083；②东方地球物理公司物探技术研究中心，河北涿州072751；③中国石化西北油田分公司勘探开发研究院，乌鲁木齐830011；④中国地质大学（北京）地球物理与信息技术学院，北京100083）

0 引言

中国陆相沉积岩大多具有薄互层特征，当受到构造运动的影响时，会发育垂向或近似垂向的高角度裂缝，呈现出典型的正交各向异性特征。作为一种重要的含油气储层，裂缝型储层中发育的裂缝既可以作为油气储集空间，又可以作为渗流通道［1］。研究正交各向异性的AVO 响应对裂缝型储层的精细刻画有重要意义。

裂缝储层的各向异性是指含裂缝地层的弹性性质在不同方向上具有差异，且与诱导裂缝的构造应力有关［2］。国外学者在各向异性领域的研究起步较早，Christoffel［3］提出用弹性张量描述介质的各向异性，因此各向异性波动方程称为Christoffel方程。Crampin［4］首次发现各向异性介质中的地震横波分裂现象，并分析了地球上常见的各向异性介质的种类与特征。在此基础上，Hudson［5-6］推导了由于裂缝存在而产生的一阶和二阶扰动量，给出了含定向裂缝介质的弹性模量的计算方法，但只适用于高频情况。Schoenberg等［7-8］在位移不连续和旋转不变性的假设下，忽略裂缝的形状与结构，提出了线性滑移理论，并对多组裂缝尺度进行了数值模拟。Thomsen［9］提出了弱各向异性理论，给出了表征弱各向异性VTI介质的参数，并据此推导了相速度和群速度的近似表达式。当介质中存在与地层走向垂直的裂隙时，使用正交各向异性模型描述更合理。Tsvankin［10］基于Thomsen的TI介质参数，分别在三个对称面内求解Christoffel方程，给出了正交各向异性介质的参数化方法，并推导了弱正交各向异性介质中qP波相速度近似表达式。基于Tsvankin提出的正交各向异性参数，Bakulin等［11-13］对两种正交各向异性模型进行了研究，分别在VTI背景下嵌入一组垂直裂缝和两组相互垂直的裂缝，运用Schoenberg的柔度建模理论构建刚度系数矩阵。Rokhlin［14］基于Christoffel方程和边界条件提出了一种计算一般各向异性介质反射系数和透射系数的方法；Rüger［15-17］推导了弱各向异性介质中反射系数近似公式，并运用纵波方位AVO 检测裂缝。

国内学者在各向异性领域的研究起步较晚，但是在正交各向异性方面的研究也有一些进展，傅旦丹等［18］研究了伪谱法在正交各向异性介质中的应用，分析了伪谱法的稳定性，并将Backus方法与Hudson理论结合，推出了周期性薄互层（PTL）与广泛扩容各向异性（EDA）组合的等效正交各向异性弹性常数的计算公式。张文生等［19］利用有限差分方法模拟正交各向异性介质中地震波场；魏建新［20］应用物理模拟方法研究了正交各向异性介质的波场特征；李娜［21］基于Tsvankin的各向异性参数研究了正交各向异性介质的相速度和群速度。许茜茹等［22］通过耦合S波射线理论和基于迭代的各向异性相速度与偏振矢量的高阶近似解，得到了适用于正交各向异性介质以qP波入射产生的二阶反射、透射系数的计算公式。张繁昌等［23］针对两组裂缝垂直的正交各向异性介质，基于Schoenberg的柔度建模理论构建刚度矩阵，推导了反射系数公式，并运用Fourier级数预测了裂缝方位。唐杰等［24］以两组垂直裂缝的正交各向异性介质为研究对象，研究了干燥和完全水饱和的地震波响应特征。郭恺等［25］研究了正交各向异性介质的多参数建模方法；单俊臻等［26］推导了一种新的PP 波HTI介质反射系数一阶扰动公式，突破了传统反射系数近似公式仅适用于小炮间距的局限。秦海旭等［27］将两组裂缝垂直斜交的各向异性地层等效为单斜介质，推导了qP波方程并进行了正演模拟。

在上述研究中，基于波动方程的正演模拟方法得到的正交各向异性的地震波场非常复杂，多种模式的波相互混杂，再加上难以完全消除的数值频散与边界效应，无法用于定量分析正交各向异性介质的AVO 响应特征。另外，弱正交各向异性假设下的反射系数近似公式的精度没有得到有效验证，无法适用于裂缝发育程度高或者VTI背景各向异性程度较强的地层。

针对上述问题，本文针对VTI背景介质中发育的一组直立裂缝诱导的正交各向异性，采用Tsvankin提出的各向异性参数构建刚度系数矩阵，再根据Christoffel方程求解出慢度矢量和偏振矢量，并结合弹性波传播理论的边界位移连续条件和应力连续条件，导出了计算精确反射系数和透射系数的方法。使用十二个各向异性模型分析背景VTI介质各向异性强度和裂缝发育程度对正交各向异性介质反射系数的影响。为了验证本文计算结果的正确性，与Rüger的反射系数近似公式的计算结果进行了对比。

1 正交各向异性介质弹性刚度矩阵

本文的正交各向异性介质是指在VTI背景介质中发育有一组垂直裂缝，如图1的下层所示。

图1中x1轴垂直于裂缝面，对称面x2Ox3平行于裂缝面，x1Ox3垂直于裂缝面，x1Ox2平行于水平界面，θ 为入射波传播方向与x3轴的夹角（入射角），φ 为入射平面与x1轴的夹角（方位角）。文中用上标“U、L”分别代表上层介质和下层介质，下标“1、2、3”对应轴x1、x2、x3轴向。相速度传播方向为波前法向，慢度为相速度的倒数。

正交各向异性介质的弹性刚度系数矩阵为

Tsvankin［9］定义的正交各向异性介质中三个不同的正交对称面上的各向异性参数为

式中：VP0为沿着垂直方向（x3方向）传播的P波速度；VS0为沿着垂直方向传播且沿着x2方向偏振的S 波 速 度；ε（1，2）、γ（1，2）、δ（1，2，3）表 示 不 同 平 面 内 的Thomsen型各向异性参数，其中上标“1、2、3”分别代表x2Ox3、x1Ox3和x1Ox2平面［11］。

可根据以上Tsvankin的各向异性参数构建弹性刚度系数矩阵，即

2 正交各向异性介质弹性波反射／透射系数精确解

Christoffel方程是由波动方程导出，用于研究地震波的慢度和偏振矢量等传播特征［28］。基于Christoffel方程及其变形，结合Snell定律求解出反射波和透射波的慢度矢量和偏振矢量，并根据边界的位移连续条件和应力连续条件，求出反射系数和透射系数精确解。

2.1 入射P波的慢度矢量计算

一般均匀完全弹性各向异性介质中地震波时间域Christoffel方程［3］为

即

式中：E 为单位矩阵；ρ为介质的密度；P 为波的偏振矢量。

根据广义虎克定律，可得

式中s1、s2、s3为慢度矢量的三个分量。

当界面上、下介质为正交各向异性或是更高对称性的介质时，假设入射波以图1方式入射，则入射波波前的法向单位矢量为

式（4）两边同乘V2／ρ可得［29］

其展开形式为

此时求解式（8）就等价于相速度V（特征值）和偏振矢量P（特征向量）的求解问题。为满足非零解条件，Christoffel矩阵的行列式值为零，即变成求解V2的一元三次方程

解式（11）可得入射波的相速度VP、VS1、VS2（VP＞VS1＞VS2），则入射P波的慢度为

2.2 反射波和透射波的慢度矢量计算

若用si表示不同波的慢度矢量（i＝0、1、…、6，分别代表入射P 波、反射P、S1、S2波、透射P、S1、S2波），则

根据Snell定律，对于所有的入射、反射和透射波，与界面平行的两个慢度矢量的分量是相等的，即

式（4）的未知解是慢度矢量和偏振矢量，对于反射波和透射波慢度矢量的分量si3，为满足非零解条件，则需求解方程

把上层介质和下层介质的参数以及已知的si1和si2代入式（16），可以解出上、下介质中的六个si3。

2.3 反射波和透射波的偏振矢量计算

若用Pi表示偏振矢量（i＝0、1、…、6，分别代表入射P波、反射P、S1、S2波、透射P、S1、S2波），即

式（16）已求得慢度矢量，再结合偏振矢量的归一化，对于上层介质有

对于下层介质有

运用Zhou等［30］的方法求解式（18）和式（19），可得到反射波和透射波的偏振矢量。

对于求得的7个偏振矢量，需要通过其与相速度的关系进行偏振矢量符号的校正。根据Chen［31］的校正方法，对于入射、反射和透射的P 波，若偏振矢量P 与慢度矢量s 的点乘大于0，则偏振矢量的符号为正；对于反射和透射的S波，入射面的法向量N 为慢度矢量s 与入射波波前法向矢量n 的叉乘，若偏振矢量P 与入射面法向量N 的点乘为正，则偏振矢量的符号为正。

2.4 反射系数和透射系数的求解

根据弹性波传播理论，当入射P 波入射到界面时，应满足边界的位移连续条件和应力连续条件。根据Rokhlin［14］给出的边界条件可以导出反射系数和透射系数计算公式，即

式中

3 模型分析

为了分析不同类型各向异性参数对反射系数的影响，先将正交各向异性分别退化到用单组裂缝弱度参数描述的HTI介质和用Thomsen参数描述的VTI介质。然后，分两种情况进行数值分析：①固定裂缝弱度参数不变，变化背景介质的VTI各向异性参数；②固定背景VTI各向异性参数不变，变化裂缝弱度参数。

为了验证本文中正交各向异性反射系数计算方法的有效性，将计算结果与Rüger在正交各向异性介质中x1Ox3与x2Ox3面的PP波近似反射系数进行对比；再将退化的HTI和VTI介质的反射系数与Rüger近似公式所得到的反射系数进行对比。

3.1 正交各向异性参数与VTI各向异性参数和裂缝弱度参数的关系

Bakulin等［12］给出了正交各向异性参数与背景介质的VTI各向异性参数以及裂缝弱度参数的关系为式中：下标“b”表示VTI背景介质；g＝V2S0／V2P0；ZN、ZV和ZH是无量纲参数，分别表示裂缝的法向弱度、垂直切向弱度以及水平切向弱度［32］。

3.2 P波从低阻抗介质入射到高阻抗介质

假设模型的上层为低阻抗各向同性介质，下层为高阻抗的正交各向异性介质。根据式（3），刚度系数矩阵可以用P 波垂直速度VP0、S 波垂直速度VS0、密度ρ以及各向异性参数表示。

模型上层介质参数为：VP0＝3150m／s，VS0＝1615m／s，ρ＝2322kg／m3，ε（1）＝ε（2）＝γ（1）＝γ（2）＝δ（1）＝δ（2）＝δ（3）＝0。

模型下层介质参数为：VP0＝3310m／s，VS0＝1697m／s，ρ＝2351kg／m3，其各向异性参数如表1或表2所示。

3.2.1 裂缝弱度参数固定不变，变化背景介质的VTI各向异性参数

首先，固定三个裂缝弱度参数不变，改变背景介质的VTI各向异性参数，分析其对反射系数的影响。表1 参数由式（23）计算得到，裂缝弱度参数ZN、ZV和ZH分别为0.215、0.120、0.090；VTI各向异性参数由弱变强：模型A1的背景为各向同性介质，εb＝0、γb＝0、δb＝0，此时正交各向异性介质退化为HTI介质；模型A2 的背景为弱各向异性VTI介质，εb＝0.1、γb＝0.12、δb＝0.07；模型A3的背景为强各向异性VTI介质，εb＝0.3、γb＝0.25、δb＝0.15。

固定θ＝20°，φ 从0°变化到90°时模型A1、A2、A3的反射系数RPP、RPS1、RPS2曲线如图2所示。可以看出：

（1）PP波反射系数的方位各向异性较弱；从模型A1到模型A3，VTI各向异性的增强导致RPP增大，模型A2比模型A1增大了40%，模型A3比模型A1增大了86%。

（2）反射系数RPS1方位各向异性很强，且随方位角增大而减小；从模型A1到模型A3，VTI各向异性的增强导致RPS1减小，最大时模型A2比模型A1减小了33%，模型A3比模型A1减小了83%；

（3）反射系数RPS2也具有很强的方位各向异性，曲线变化趋势与RPS1沿着方位呈大体镜像关系，但是量值上有差异；从模型A1到模型A3，VTI各向异性的增强导致RPS2减小，最大时模型A2 比模型A1减小了30%，模型A3比模型A1减小了70%。

表1 裂缝弱度参数不变、背景介质的VTI各向异性参数变化时的正交各向异性参数

固定φ＝60°，模型A1、A2、A3的反射系数RPP、RPS1、RPS2随θ的变化曲线如图3所示。由图可知：

（1）小入射角时（θ＜10°），反射系数RPP受背景介质VTI各向异性参数的影响不大。但是，随着θ的增大，模型A1 的RPP随θ 增大而减小，且在θ≈35°时发生极性转变；从模型A1到模型A3，VTI各向异性增强导致大入射角的RPP增大。

（2）随着θ的增大，反射系数RPS1在θ＜30°时都呈增大的趋势，而在30°～45°范围内模型A2、A3开始减小，模型A3在θ≈35°出现了极性转变；从模型A1到模型A3，VTI各向异性参数的增大导致RPS1减小，在θ＝30°时模型A2比模型A1减小了50%，模型A3比模型A1减小了75%。

（3）随着θ的增大，反射系数RPS2与RPS1变化趋势一致；从模型A1到模型A3，VTI各向异性系数的增大导致RPS2减小，在θ＝30°时模型A2 比模型A1减小了37.5%，模型A3比模型A1减小了75%。

φ 在0°～90°内变化、θ在临界角范围以内，P波由各向同性介质入射到正交各向异性模型A1、A2和A3时的反射系数RPP、RPS1和RPS2的等值线图分别如图4、图5、图6所示，可以看出：

（1）随着VTI背景各向异性程度的增强，RPP等值线曲率变化不大，只在反射系数的量值上有整体增大，且靠近临界角时增加幅度最大。

（2）RPS1和RPS2随方位变化的程度远大于RPP，都有随着背景介质VTI各向异性强度增强而减小的趋势，且等值线曲率变化也较大。

由于背景介质中插入的裂缝弱度值不大，裂缝诱导的HTI各向异性属于弱各向异性。针对此类正交各向异性介质的AVO 分析，单纯的PP波对裂缝的发育以及背景介质的VTI各向异性程度的敏感度远小于PS1波和PS2波。

图3 φ＝60°时模型A1、A2和A3的R PP（a）、R PS1（b）、R PS2（c）随θ的变化曲线

图4 入射P波由各向同性到各向异性模型A1（a）、A2（b）、A3（c）反射系数R PP等值线图

图5 入射P波由各向同性到各向异性模型A1（a）、A2（b）、A3（c）反射系数R PS1等值线图

图6 入射P波由各向同性到各向异性模型A1（a）、A2（b）、A3（c）反射系数R PS2等值线图

3.2.2 背景介质VTI各向异性系数固定不变，变化裂缝弱度参数

由式（23）可知，当固定背景VTI介质的三个各向异性参数不变（εb＝0.1，γb＝0.12，δb＝0.07，与表1中模型A2的背景介质VTI各向异性参数取值相同）时，改变三个裂缝弱度值会影响HTI各向异性的强弱，弱度越大则HTI各向异性越强。表2为下层高阻抗介质的背景VTI各向异性参数固定不变、变化裂缝弱度参数计算的正交各向异性参数值。模型B1的裂缝弱度参数ZN＝ZV＝ZH＝0，此时正交各向异性介质退化为VTI各向异性介质；模型B2的裂缝弱度参数ZN＝0.215、ZV＝0.120、ZH＝0.090，呈弱HTI各向异性，等同于表1中的模型A2；模型B3的裂缝弱度参数ZN＝0.215、ZV＝0.120、ZH＝0.090，裂缝发育程度较高，HTI各向异性也较强。

固定θ＝20°，φ 在0°～90°变化时模型B1、B2、B3的反射系数RPP、RPS1、RPS2曲线如图7所示。由图7可知：

（1）裂缝的发育会导致PP 波的反射系数变小，模型B2比模型B1减小了32%，模型B3比模型B1减小了64%；即使在裂缝弱度较强的情况下，PP 波的反射系数的方位各向异性也是比较微弱的。

表2 背景VTI各向异性参数不变、裂缝弱度参数变化时的正交各向异性参数

图7 θ＝20°时模型B1、B2和B3的反射系数R PP（a）、R PS1（b）、R PS2（c）随φ 的变化曲线

（2）反射系数RPS1受φ 的影响较大，在φ＝0°时，即平行裂缝面时，RPS1最大；在φ＝90°时，即垂直裂缝面时，RPS1为0；从模型B1到模型B3，裂缝弱度的增强导致RPS1增大，最大时模型B2比模型B1增加了300%，模型B3 比模型B1 增加了600%。当裂缝不发育时（模型B1），RPS1不存在方位各向异性，等同于PSV 波的反射系数。

（3）反射系数RPS2在φ＝0°时最小，φ＝90°时最大；从模型B1 到模型B3，裂缝弱度的增强导致RPS2增加，最大时模型B3比模型B2增加了100%。

固定φ＝60°时，模型B1、B2、B3 的反射系数RPP、RPS1、RPS2随θ 的 变 化 曲 线 如 图8 所 示，由 图可知：

（1）θ 在0°～45°（临界角以内）范围内，反射系数RPP在θ＜5°时无明显差别，随着θ增大，RPP表现不同，裂缝弱度越大，RPP减小越快；最大时模型B2比模型B1减小了75%，模型B3比模型B1减小了

100%。

（2）随着θ增大，RPS1在θ＜25°时都呈增大的趋势，而在25°～45°范围内模型B1、B2 开始减小，模型B1在θ≈30°出现了极性转变；从模型B1到模型B3，裂缝弱度的增强导致RPS1增大，θ＝25°时模型B2比模型B1增大了100%，模型B3比模型B1增大了300%。

（3）随着θ 增大，当裂缝发育时，RPS2都呈负极性方向增大的趋势，在θ＝0°时最小，在θ＝45°时最大；从模型B1 到模型B3，裂缝弱度的增强导致RPS2增大，最大时模型B3比模型B2增大了100%。

φ 在0°～90°内变化、θ在临界角范围以内，P波由各向同性介质入射到正交各向异性模型B1、B2和B3时的反射系数RPP、RPS1和RPS2的等值线图分别如图9、图10、图11所示，可以看出：随着裂缝弱度的增强，RPP在正极性方向减小，RPS1和RPS2在负极性方向增大。总体上，即使裂缝诱导的HTI各向异性较强，RPP的方位各向异性也比较弱，其对于裂缝发育的敏感度远小于RPS1和RPS2。

图8 φ＝60°时模型B1、B2和B3的反射系数R PP（a）、R PS1（b）、R PS2（c）随θ的变化曲线

图9 入射P波由各向同性到各向异性模型B1（a）、B2（b）和B3（c）的反射系数R PP等值线图

图10 入射P波由各向同性到各向异性模型B1（a）、B2（b）和B3（c）的反射系数R PS1等值线图

图11 入射P波由各向同性到各向异性模型B1（a）、B2（b）和B3（c）的反射系数R PS2等值线图

3.3 P波从高阻抗介质入射到低阻抗介质

为了进一步探讨界面上、下阻抗对反射系数的影响，将A、B系列模型的上、下层介质的纵横波速度和密度交换，下层的各向异性参数保持不变，可得模型C1、C2、C3（表1，裂缝弱度不变，VTI各向异性强度变化）和模型D1、D2、D3（表2，裂缝弱度变化，VTI各向异性参数不变）。

固定θ＝20°，模型C1、C2、C3的反射系数RPP、RPS1、RPS2的随φ 的变化曲线如图12所示；固定φ＝60°，模型C1、C2、C3的反射系数RPP、RPS1、RPS2随θ的变化曲线如图13所示。对比图12、图13与图2、图3可知，与下层介质为高阻抗介质的情况相反，当下层介质为低阻抗介质时，背景介质的VTI各向异性强度越强，则PP 波的反射系数在负极性方向越小，而PS1波与PS2波在正极性方向越大。

固定θ＝20°，模型D1、D2、D3的反射系数RPP、RPS1、RPS2随φ 的变化曲线如图14 所示；固定φ＝60°，模型D1、D2、D3的反射系数RPP、RPS1、RPS2随θ的变化曲线如图15所示。对比图14、图15与图7、图8可知，裂缝的发育会导致PP 波的反射系数在负极性方向变大，而PS1波和PS2波的反射系数则在正极性方向变小，与下层介质为高阻抗介质的情况相反。

无论是从低阻抗介质入射，还是高阻抗介质入射，PP波反射系数的方位各向异性都远小于PS1波和PS2波。

3.4 与Rüger近似反射系数公式的对比

为了验证本文反射系数计算结果的正确性，与Rüger近似公式计算的正交各向异性反射系数、HTI介质反射系数和VTI介质反射系数进行对比。

图12 θ＝20°时模型C1、C2和C3的反射系数R PP（a）、R PS1（b）、R PS2（c）随φ 的变化曲线

图13 φ＝60°时模型C1、C2和C3的反射系数R PP（a）、R PS1（b）、R PS2（c）随θ的变化曲线

图14 θ＝20°时模型D1、D2和D3的反射系数R PP（a）、R PS1（b）、R PS2（c）随φ 的变化曲线

图15 φ＝60°时模型D1、D2和D3的反射系数R PP（a）、R PS1（b）、R PS2（c）随θ的变化曲线

3.4.1 正交各向异性介质

对于P波入射，Rüger［17］仅推导了正交各向异性介质RPP在x1Ox3和x2Ox3面的近似反射系数公式。将模型A2的参数代入到其公式计算反射系数曲线，并与本文的反射系数精确解的计算结果进行对比（图16）。在x1Ox3面，两条曲线变化趋势一致，当θ在0°～20°范围内时近乎重合（图16a）；在x2Ox3面，两条曲线变化趋势亦一致，当θ 在0°～25°范围内时近乎重合（图16b）。

3.4.2 退化到HTI介质

依据式（23）可以将正交各向异性介质退化到HTI介质。将模型A1代入Rüger［17］的HTI介质反射系数近似公式计算反射系数曲线，并与本文的反射系数精确解对比（图17）。由图17 可以看出RPP、RPS1和RPS2的Rüger近似解与本文精确解变化趋势一致；RPS2在0°～20°范围内近乎重合。

3.4.3 退化到VTI介质

依据式（23）可以将正交各向异性介质退化到VTI介质。将模型B1 参数代入Rüger的VTI介质反射系数近似公式，并与本文的反射系数精确解对比（图18）。从图18可以看出，RPP和RPS的近似解与精确解变化趋势一致，RPP在0°～25°范围内完全重合，RPS在0°～15°内近乎重合。

图16 正交各向异性介质Rüger近似公式与本文方法计算的反射系数曲线对比

图17 HTI介质Rüger近似公式与本文方法计算的反射系数曲线对比

图18 VTI介质Rüger近似公式与本文方法计算的反射系数曲线对比

由以上分析可知，本文反射系数精确解与Rüger近似公式在中、小角度范围内近乎重合，整体趋势一致。因此，本文反射系数精确解计算方法正确。

4 结论

针对互层背景中发育一组裂缝的正交各向异性介质模型，本文基于Christoffel方程及弹性界面的位移连续条件和应力连续条件，推导了反射系数和透射系数精确求解方法，建立了正交各向异性参数与背景介质VTI各向异性参数系数和裂缝参数之间的关联。通过与Rüger反射系数近似公式计算结果的对比，验证了本文计算结果的正确性。理论模型计算结果表明：

（1）当P 波从低阻抗介质入射至高阻抗介质时，背景介质VTI各向异性参数的变大会导致PP波反射系数增大、PS1 波和PS2 波的反射系数减小；而裂缝弱度的增强会导致PP 波反射系数的减小及PS1波和PS2波的反射系数的增大。

（2）当从高阻抗介质入射至低阻抗介质时，背景介质VTI各向异性系数的变大会导致PP波反射系数减小、PS1波和PS2波的反射系数增大；而裂缝弱度的增强会导致PP 波反射系数的增大、PS1 波和PS2波的反射系数的减小。

（3）PP 波反射系数的方位各向异性要远小于PS1波与PS2波。