地震各向异性研究进展(81届SEG年会)-1st稿件1.doc

地震各向异性研究进展
一、引言
2011年第81届SEG年会“各向异性”专题涉及文章32篇（含张贴报告），讨论的内容主要包括各向异性介质波动方程解耦、地震波在各向异性介质中的传播、各向异性介质偏移和成像方法研究、各向异性介质时差速度反演与各向异性参数建模方法、各向异性弹性参数反演与裂缝预测、近地表各向异性分析等诸多方面。
在本文中，我们对专题论文涉及的主要内容进行了归类，总结了各向异性介质地震波传播、成像与反演方面的研究进展。
二、各向异性介质伪P波正演和逆时偏移
地震波在各向异性介质中传播是一种弹性现象，由于弹性各向异性介质波动方程十分复杂，所以大多数各向异性正演和偏移方法通常使用伪纵波方程。
Chu等（2011）给出了使用伪差分算子的泰勒级数展开近似解耦得到纯P波方程，级数展开产生了一个新的纯纵波各向异性方程如式（1）所示，对于一定范围的波数和各向异性参数都是收敛的。
（1）
使用有限差分法求解导出的方程（1），并用纯P波方程进行数值模拟，其结果与伪纵波方程的结果十分接近，精度较高，并且能保证完全不受横波的影响。图1为BP 倾斜各向异性（TTI）模型的一部分，(a)为速度模型，（b）、（c）、（d）分别为各向异性参数模型，其中各向异性参数变化很大。模拟结果如图2所示，其中（a）为伪声波波动方程（）的模拟结果，（b）为伪声波波动方程（）的模拟结果；（c）为纯声波波动方程模拟结果；（d）为图（c）与图（b）叠加结果。从图2中可以看出当横波速度为零时，用伪声波波动方程产生的结果不稳定，引入非零S波速度后的结果稳定，但是用该方程得到的结果包含了震源产生的S波和P-SV转换波。纯纵波方程能够解决这个问题而不涉及稳定性问题，并且不会产生任何的SV波。

(a) (b)

(c) (d)
图1 BP TTI模型（部分）
(a) (b) (c) (d)

(a) (b)

(c) (d)
图2 用不同的波动方程计算的BP TTI模型1.5秒处的波场快照
(a)伪声波波动方程，其中；(b)伪声波波动方程，其中；
(c)纯声波波动方程；(d)图c(黑色)叠加在图b(红色)之上
Paul J. Fowler（2011）将横向各向同性介质中的伪P波正演和逆时偏移的方法推广到正交各向异性介质中，推导了正交各向异性介质中耦合偏微分方程如式（2）所示，可以准确地模拟正交各向异性介质中伪P波的传播，无需知道准确的横波速度。图3为使用伪P波一阶耦合方程组得到波场快照。其中（a）为横测线方向的波场快照，(b)为纵测线方向的波场快照，（c）为深度切片，模型参数为，，，。从图中可以看出明显的三叉区现象，表明该模型存在强烈的各向异性现象。
（2）

(a) (b)

(c)
图3 使用伪P波一阶耦合方程组得出的波场快照
(a)主测线; (b) 联络测线; (c)水平切片
Zhang（2011）推导了一组声波方程，如式（3）所示，在垂直和倾斜的正交各向异性介质中进行正演和逆时偏移，数值计算表明该方法在更为常见的各向异性介质中可提供稳定且高质量的逆时偏移结果。
（3）
其中是密度函数，是质点位移矢量,是包含三个主应力的矢量，是体力矢量。在方程（7）中，对角矩阵被定义为。是3×3的弹性常数矩阵。，其中
（4）
Zhang计算了TTI介质和倾斜正交各向异性介质的逆时偏移脉冲响应（图4），其中（a）、（d）为纵测线方向的逆时偏移脉冲响应，（b）、（e）为横测线方向的逆时偏移脉冲响应，（c）、（f）为深度为1.5km的波场快照。TTI介质和倾斜正交各向异性介的参数定义为，，，，，倾斜的坐标轴为平面的旋转角度为。无论在平面还是在平面，倾斜正交各向异性介质的脉冲响应与TTI介质的脉冲响应相比都有明显的运动学差别。深度切片对比表明倾斜正交各向异性介质逆时偏移的格林函数更加不对称。

图4 TTI介质（a，b，c）和倾斜正交各向异性介质（d，e，f）的逆时偏移脉冲响应
(a)(d)主测线; (b)(e) 联络测线; (c)(f) 1.5km处深度切片
图5为倾斜正交各向异性的复杂三维盐丘模型中的逆时偏移脉冲响应。速度模型如图5a、d所示。Thomson参数定义为，定义。图5c为横测线方向的脉冲响应，图5f为纵测线方向的脉冲响应，结果表明方程（3）对于复杂的倾斜各向异性介质成像是稳定的。

图5 在具有倾斜正交各向异性的复杂三维盐丘模型中的逆时偏移脉冲响应
(a)(b)(c) 主测线方向的模型、模型及RTM脉冲响应
(d)(e)(f) 联络测线方向的模型、模型及RTM脉冲响应
三、各向异性NMO时差速度和各向异性参数建模
1、各向异性NMO时差速度
Takanashi等（2011）推导出了在分层弱横向速度变化介质中NMO椭圆的解析表达式，提出了有效估算间隔时差参数的算法，并讨论了引起透镜体畸变的原因。对于三层弱各向异性层，中间层为强横向不均匀介质（如图6），从底部反射的反射波NMO椭圆可以表示为：
（5）
（6）
（7）
（8）
式中为水平坐标向量，是2×2的对称矩阵，代表有效的NMO椭圆，是横向均匀介质中对于CMP道集的NMO椭圆，是从模型底部垂直反射的反射波单程走时，是第层的间隔走时，是间隔NMO椭圆的等效速度（各向同性）。公式是在地层的每一层都有水平对称面（正交或单斜对称）的假设条件下推导出来的。如果速度变化不大，和比较接近，分别为第二层和第三层的厚度。

图6 通过3层水平层的反射波路径
R是未受干扰的自然反射点，Grechka和Tsvankin在1999年提出该模型
Takanashi等（2011）研究了一个拉长的速度透镜体对NMO椭圆的影响，并提出了有效估算间隔时差参数的算法，准确地消除透镜体引起的有效NMO椭圆畸变（图7）。

(a)

(b)
图7 在透镜体中心位置的非双曲时差校正CMP道集
(a)使用透镜体转换走时之前，(b) 使用透镜体转换走时之后
Casasanta等（2011）提出对变换后的三维CMP道集自动的不依赖NMO椭圆的校正方法，通过直接求出数据的局部斜率来获得同相轴几何形态。域处理的优点在于它可以用旅行时剥离方法替换Dix型微分，这会极大地简化层间动校正参数的估计（图8和9）。
Casasanta等得到NMO椭圆剥离方程，它们是Casasanta和Fomel（2010）推导的二维方程拓展到三维（小角度）形式：
（9）
（10）
（11）
图8（a）显示了一个应用时变有效NMO椭圆的逆三维域NMO合成得到的CMP数据。速度矩阵V的平方根分量组成一个正梯度加上一个正弦趋势。，，。图8（b）、（c）是沿和方向得到的横测线和纵测线数据的局部斜率。图9（a）所示是将斜率变换为零斜率时间体，将给定的地震道样点对应到它的零斜率NMO校正时间。因此，无需对每一个地震道进行时移，可以在没有使用NMO椭圆信息的情况下自动校平道集，图9（b）表明使用单一速度的各向同性圆形NMO校正不能校平所有测线方位上的数据，而图9（c）证实了该方法能够成功地解决存在方位各向异性的三维时差问题。

图8 (a)正演合成的得到的域CMP道集
用PWD算法估计得到的横测线(b)和纵测线(c)的局部斜率

图9 (a)从和中获取零斜率时间()体，(b)各向同性（）圆形时差动校正的三维CMP道集，(c)利用时间体进行不依赖NMO椭圆的三维域时差动校正后数据
2、各向异性参数建模
各向异性深度偏移成像要求横向各向同性介质（TI）或者倾斜横向各向同性介质（TTI）速度模型，对于TI介质中的P波的传播，需要三个各向异性参数（速度和两个各向异性参数）；对于TTI介质，需要引进与倾角有关的参数。
Li等（2011）为了改进深度定位和构造精度，对墨西哥湾花园海岸东北部建立了各向异性速度模型；通过假设垂直横向各向同性（VTI）并通过比较地震和井数据推出了Thomsen各向异性参数；连井和校验炮数据指示了花园海岸东北部两层的各向异性参数的变化趋势：浅部低值和深部高值。利用这两层模型（图10）和扫描的Thomsen各向异性参数模型，进行了各向异性（VTI）成像。
计算的完整公式为
（12）
其中为浅部低值，深部高值。为了平滑两层之间的过渡，在上作用一个阻尼因子。从水底到拾取的水平层，从0到1呈四次幂增加。在较深层，恒等于1。

图10井-1中三维两层趋势模型的inline剖面
红点是井标注层，蓝点是高值层位的顶部
图11（b）显示过井-1各向异性（VTI）成像剖面，与各向同性成像图11（a）相比，各向异性成像最重要的改善是在提高成像质量的同时增加了深度精度。各向异性（VTI）地震成像与各向同性成像相比，在VTI成像中的地震标注层（破折粉红线）几乎与井标注层（红点）完美连接。

图11(a)井-1各向同性成像剖面。井标注层为红点，破折粉红线显示的是地震中相应的标注层。蓝点是高值层的顶部。红色水平线是水底和质量传输集合体的基线。 (b)是(a)对应剖面的VTI成像剖面(Courtesy of CGGVeritas Sevices Inc., Houston, Texas)
Bachrach等（2011）将岩石物理模型预测与由墨西哥湾18口井经过校验炮数据和实施了地震动校之后共同分析所得到的各向异性参数估计进行了对比，结果表明基于岩石物理的各向异性速度建模是可行的（图12、13、14）。

图12 岩石物理各向异性速度模型建立（VMB）工作流程图. 输入的是速度、孔隙度转换和估计的地震速度. 利用基于压实页岩模型得到的热度信息和各向异性图像孔隙度的压实参数进行内插计算。输出为与页岩压实有关的估计的Thomsen参数. 有关模型参数和方程的完整描述参考Bachrach（2010）

图13 调查区域：顶面-绿峡谷地区70墨西哥湾区块的18口井的覆盖区域
底面：用于输入岩石物理模型的层析成像速度

图14 由岩石物理模型（左）得到的和由内插的校验炮数据（右）得到的Epsilon和delta值
注意在深部区域epsilon值的差异
岩石样品的各向异性参数可以在在实验室中测得。由于在实验室的声波实验是在非各向同性应力场下进行的，对称轴方向可能会改变，这将会违背具有垂直或水平对称轴方向的横向各向同性假设。考虑这些因素，促使人们要考虑对称轴与层理面径向和层理面法向的偏离，这不仅使实验室声波实验更自由，而且可以跟踪测量压力作用下对称坐标可能发生的改变。Dariush Nadri（2011）使用了一种层析型观测系统，用来测量圆柱形岩石样品的超声波射线速度（图15、图16）。利用这些速度，通过快速模拟重退火算法（VFSR）找出真实解，最后估计出Thomsen各向异性参数和对称轴坐标。

图15 在一个圆柱形页岩样品上激发器和接收器的排列中射线传播的原理图
当一个传感器做发射源时，其它的传感器为接收器

图16 射线向量（OA），慢度向量（OB），对称轴向量（OC）的原理图
CB是慢度向量在各向同性平面上的投影，是水平慢度（1997，Tsvankin修改后）
Downton（2011）讨论了不同方位模型的速度估计误差，即多种偏移距-方位角组合关系的误差，包含方位角集和共偏移距向量（COV）集，提出了一种利用全方位角的曲面拟合技术来反演HTI模型的速度，该方法能利用可靠的方位动校正速度来拉平CMP道集。通过对阿尔及利亚WAZ地区的陆上地震资料的实例应用，证明了该技术的优势（图17）。

(a) (b)

(c) (d)
图17 (a)沿层各向异性参数 (b)沿层快速度方位角 (c)沿层快速度 (d)沿层快速度误差
四、各向异性介质偏移成像
Zhou等（2011）利用最新发展的叠前深度偏移技术，包括TTI（倾斜各向同性）模型建立和高保真控制束偏移（HFCBM），对LD区块进行了完整的重新处理。结合浅水去除多次波（SWD）技术，再处理使这个地区的成像质量大幅提高，可以看到在构造图像和信噪比上都得到了较大改善，如下图所示。

图18 Kirchhoff和HFCBM剖面浅部断层区域的深度切片对比

图19 Kirchhoff和HFCBM剖面浅部断层区域的叠加对比
图18为Kirchhoff和HFCBM剖面浅部断层区域的深度切片对比示意图，两种方法输入相同的速度模型，比较左图与右图可以看出对于利用HFCBM求取的结果对于确定浅部油气目标有很大的改善作用。在Kirchhoff剖面中很难解释浅部断层的油气目标，然而在HFCBM结果中，成像结果的质量和信噪比得到了较大改善。图19为Kirchhoff和HFCBM剖面浅部断层区域的叠加对比，比较左图和右图，可以得到同样的结论。
Valenciano等（2011）提出了一种能够对各向异性衰减介质准确成像的单程波波动方程方法。在偏移过程中为了补偿Q的影响，采用扩展的傅里叶有限差分（FFD）算法（Ristow和Ruhl，1994）。该混合域的单程波波动方程算法是灵活且有效的，特别对岩下复杂地质构造的成像是准确的。利用FFD粘声PSDM方法可以很好的控制衰减（图20和图21）。
归一化波函数（二维情况）的频散关系可用下式近似：
（13）
其中是角频率，是垂直波数，和是水平波数，和是参考频率处的品质因子和相速度，，，。可以用来在有限差分中采用奇偶有理函数项来解释TTI各向异性（Valenciano等人，2009），暗示了对于所有Thomsen各向异性参数、、极化方位角，以及的组合，可以计算得到精确频散关系的准确解。有理函数的系数可以利用加权最小均方拟合得到。
对北海3D数据体应用了具有Q补偿的偏移，深度Q模型利用频谱比方法得到。图21显示了是未进行和进行了Q补偿偏移后的结果，可以看到在进行了Q补偿偏移之后，分辨率得到了明显地改善。

图20 北海数据的深度偏移：顶部表示未进行Q补偿，底部表示进行了Q补偿

图21墨西哥湾数据的深度偏移：顶部表示未进行Q补偿，底部表示进行了Q补偿
五、各向异性参数反演和裂缝预测
1、各向异性旅行时参数反演
Alkhalifah（2011）提出了一种水平横向各向同性（HTI）介质的旅行时近似式来对对各向异性参数进行简化，该近似给出了在各向同性背景下的旅行时的精确解析表达式，使用现有的动校正方程，可以将非均匀背景下的椭圆各向异性模型拓展到横向对称轴方向缓慢变化的HTI介质中。Golikov（2011）推导出了TTI介质旅行时参数的准确公式和弱各向异性近似公式。证明弱各向异性近似对于较小的TTI介质是适用的，还实现了弱各向异性近似中旅行时参数的解析反演。
Liu等（2011）提出了对于各向异性介质的三维动态规划方法，与三维各向同性介质一样，每一个当前计算点的新旅行时用4个预先算出的邻点旅行时计算出来（图22~图25）。图22是在中国西南部山区勘探采集到的一个三维真实数据的起伏地表。这个地区的高度差超过了1000m。图23和图24显示了在地表下相同波的传播现象，说明该方法在有地表存在时不会违背地震波传播的规律。图25是当考虑TTI介质倾角时得到的三维旅行时结果。

图22 油田数据的起伏地表

图23 震源位于油田数据地表得到的三维旅行时结果

图24 震源位于平面上得到的三维旅行时结果

( a ) ( b )
图25TTI介质（倾角为-45度）中z-x平面的切片
(b) TTI介质（倾角为45度）中z-x平面的切片
2、各向异性方位AVO分析
Yao等（2011）基于Rüger（1998）提出的水平横向各向同性（HTI）介质的概念，结合倾斜界面入射角、倾斜角、方位角、响应水平层的入射角之间的关系，提出了对于两个介质间倾斜反射层的AVO/AVA分析流程，两个介质中的一个是弱HTI介质。
HTI介质中倾斜层对应反射振幅的计算公式（Yao等，2011）为：
（14）
其中是两个新的符号，记作两个地层间倾斜反射层的影响（其中一个地层是HTI地层）
（15）
（16）
Yao等采用中国渤海渤中区域一个气田的一套数据来说明倾斜界面对目的层的影响，该例子的进一步的目的是从原始数据中消除倾角效应，并且恢复与地下岩石物理性质直接相关的真实反射系数。图26所示的是一个探井和一个1225到1230米目的含气区的地震叠加剖面。

图26 井A和井间地震数据
图27给出了最终结果，蓝点线代表由真实数据计算出的反射系数，红点线代表由井数据得到的岩石物理参数所求出的水平反射层反射系数的理论值。蓝点线是考虑倾角影响时，求得的理论反射系数；紫色线是真实的不受倾角影响的反射系数，可见地震反射系数随炮检距的增加而减小。钻前的地质分析表明这一勘探区有成为气藏的潜在构造，有利于油气的聚集。A井数据进一步证明了气的存在，这与消除了地层倾角影响的钻前AVO分析结果一致。

图27 从岩石物理模型和原始测试数据中计算得到的AVO响应
蓝点线代表由真实数据计算出的反射系数
红点线代表由井数据得到的岩石物理参数所求出的水平反射层反射系数的理论值
3、裂缝参数估计与裂缝储层预测
Downton等（2011）利用线性滑移理论（Schoenberg，1980）对傅里叶系数进行解释，基于含有平行对称轴的HTI介质的假设，给出了独立的傅里叶系数分量的物理含义，将多个傅里叶系数结合，不仅可以得到各向异性梯度的无偏估计和对称轴的明确估计，还可以分别得到法向与切向弱度参数（图28）。

图28 傅里叶系数方法与近偏移距Rüger方程计算的各向异性梯度和各向同性平面方位角结果比较
（a）傅里叶系数，（b）加权后的傅里叶系数，（c）使用近偏移距Rüger公式计算的各向异性梯度（d）加权后的傅里叶系数，（e）根据傅里叶系数的相位计算的方位，（f）使用Rüger公式计算的方位
注：数据自加拿大的三维地震数据的一条横测线，傅里叶系数由平均入射角为35度的部分叠加方位剖面计算得到，图（a）相当于35度部分角度叠加剖面，图（b）与图（c）吻合很好，图（d）含有独立于图（b）的信息，但两者仍有许多相似之处，图（e）与图（f）吻合很好
Sun等（2011）介绍了一种使用有限方位分布下的叠前数据进行纵波裂缝预测的算法，并将这个算法应用到塔中45区域，反演得到的结果被证明是可信的，并且与地质统计数据匹配的很好。
传统的基于Ruger纵波反射系数方程（Ruger，1998）的反演方法进一步推导为：
（17）
式中是裂缝密度，表征各向异性程度的是四个未知量。假设共反射点叠加次数是，则式（17）能写成：
（18）
式中是一个的矩阵；是一个包含的列向量；是地震反射振幅（不含噪声）。传统的裂缝反演方法是计算未知量。为了避免噪声干扰，需要使用宽角度数据，存在噪声时会使求解不稳定，Sun等在公式（18）中加入噪声为：
（19）
通过奇异值分解方法，公式中的矩阵A能写成：
（20）
通过奇异值分解法计算广义逆得到：
（21）
通过利用公式（21），有效信号能从含噪信号中分离出来，噪声对反演结果的影响得到限制，使得能够利用有限方位数据进行反演。
图29是塔中45区采集方位分布的概要图，这里东向定义为零度角，方位角按逆时针方向增加。激发源和接受组合分布在有限方位范围内（20度到80度）。根据这一区域偏移距的分布，方位角分为4部分：20~35度，35~50度，50~65度，65~80，计算时用每个角度范围的中值，即27.5，42.5，57.5，72.5度。

图29 塔中45区采集方位分布的概要图。激发源和接收组合分布在有限方位的范围内
将使用有限方位地震数据进行裂缝检测的方法应用于塔中45区奥陶系碳酸盐岩储层。图31中可以看出一些裂缝发育层的特征：方位各向异性在构造的西北侧翼要强，特别是在大的断层之间的区域。这些裂缝不能被传统的三维相干技术识别出来。

图30 展示振幅各向异性的方位角道集

图31 塔中45区反演得到的裂缝密度和方位的结合（图中白线是叠加的解释得到的断层）
六、近地表各向异性的影响
Samik Sil（2011）研究了两个合成数据的例子来分析近地表各向异性对深部各向异性目标的影响，结果表明由于近地表各向异性的存在，目的层反射波的旅行时数据（或VVAZ信号）畸变很严重（图32），但是对目的层反射波的AVAZ信号的影响是可以忽略的。
为了描述近地表各向异性介质对深部各向异性油藏的影响，使用了根据表1和表2中数据建立的六层模型。在模型1中，层3和层5是各向异性介质（HTI介质）。层3表示近地表各向异性介质。层5是一个各向异性油藏。与另一个模型（模型2）对比，模型2的上覆层（层3）是各向同性的。
表1 模型1参数
宽度(km)
Vp(km/s)
Vs(km/s)
(gm/cc)
E

0.5
3.0
1.5
2.45
0
0
0
N/A
0.5
3.5
1.7
2.45
0
0
0
N/A
0.6
4.0
2.0
2.45
0.05
0.05
0.08
0
0.5
4.3
2.15
2.5
0
0
0
N/A
0.6
4.7
2.2
2.55
0.2
0.1
0.15
90

5.2
3.0
2.6
0
0
0
N/A
表2 模型2参数
宽度(km)
Vp(km/s)
Vs(km/s)
(gm/cc)
E

0.5
3.0
1.5
2.45
0
0
0
N/A
0.5
3.5
1.7
2.45
0
0
0
N/A
0.6
4.0
2.0
2.45
0
0
0
N/A
0.5
4.3
2.15
2.5
0
0
0
N/A
0.6
4.7
2.2
2.55
0.2
0.1
0.15
90

5.2
3.0
2.6
0
0
0
N/A

图32（a）是由模型1得到的合成方位角道集，图32（b）是由模型2得到的相似的合成数据。通过这些图可以清楚地发现与5个界面相对应的5个同相轴。以下是对两个模型每一个同相轴的分析。模型1和2对应的第一个同相轴由第一层和第二层各向同性介质间的界面产生，这些同相轴很平，在这个同相轴的旅行时和振幅上没有观测到任何典型的各向异性余弦形变化；图32中的第二个同相轴是模型1和2中第二层和第三层介质间界面产生的，由于模型1中的第三层是各向异性介质，对照图32（a）和（b）中的第二个同相轴，可以清楚地看到由方位各向异性产生的典型振幅余弦形变化。图32（a）中第三个同相轴是模型1中各向异性-各向同性界面产生的，图32（b）中第三个同相轴是模型2中各向同性-各向同性界面产生的，可以看到地震波穿过各向异性介质后，模型1中第三条同相轴的旅行时和振幅都出现了典型的方位变化，而由于模型2中的界面是各向同性到各向同性界面，故合成地震剖面上没有这些变化；同相轴4在模型1和2中都是来自一个各向同性到各向异性的界面，对于模型1可以看到一个与同相轴3相似的旅行时正弦形变化同相，注意到尽管界面上层介质是各向同性的，观测到的旅行时变化使各向同性地层显出各向异性的特征，这个明显的各向异性现象是由于纵波穿过一个各向异性地层到达检波器产生的。正如期望的那样，由于没有上覆各向异性地层，在模型2中没有观测到任何明显的各向异性现象。

图32 炮检距范围为2.41到2.44km对应的方位角道集图左排使用上覆各向异性层的模型1。图右排是上覆各向同性层的模型2。五个同相轴由表1和表2描述的六层模型产生
结束语
在地震资料的处理和解释过程中，忽略各向异性将会导致异常大的偏差，考虑地下介质的各向异性则可以显著地改善地震资料处理的质量，因此，地震各向异性的研究已经成为地球物理工作者必须掌握的基础技术手段。其中的关键是如何正确估计各向异性参数，对成像参数与反演参数的深入研究是各向异性研究的重要课题，该项研究的深入开展将会进一步推动地震各向异性处理和解释技术的进步。
22