黄国娇 孙江兵 白超英 钱 卫

(①河海大学地球科学与工程学院地质科学与工程系，江苏南京 211100； ②江苏省地质资料馆，江苏南京 210018； ③长安大学地质工程与测绘学院地球物理系，陕西西安 710054)

1 引言

地震层析成像技术经过近40年的不断发展和完善，已成为研究地球内部结构以及各向异性的有效技术之一。随着油气等矿产资源的不断开发以及对地壳和上地幔结构的深入研究，对地震成像分辨率的要求也越来越高。

传统的各向异性层析成像是在弱各向异性假设下，主要利用走时信息进行线性反演。如：Hirahara等[1]利用qP波走时反演三维各向异性参数；Jech[2]进行了三维 TI介质走时反演；Chapman等[3]提出了一种适用于各向异性介质中的线性化走时层析成像方法，并应用于井间层析成像。上述研究以弱各向异性假设为前提，即所有正、反演理论均基于各向同性介质背景下的一阶扰动而建立，简化了数学计算的难度，且在反演迭代中不用重复进行射线追踪计算(雅克比矩阵元素固定不变)，从而有利于进行线性反演重建介质属性 。然而，弱各向异性假设忽略了介质横向变化，反演结果可信度不高。

近年来，各向异性层析成像已经发展到利用多种地震信息进行非线性反演。如前人在遗传算法、非线性解析法波动方程反演等方面进行了大量的理论研究[4-6]；Zheng等[7]利用qP波的慢度和质点振动测量值确定弱各向异性参数；阮爱国等[8]推导了弱各向异性介质中的地震波速，同时介绍了反演上地幔各向异性参数的具体算法；郭飚[9]从弱各向异性条件出发构建地震各向异性P波走时层析成像算法；傅旦丹等[10]采用坐标扰动法反演一维正交各向异性介质的模型参数；陈竟一等[11]采用模拟退火方法讨论了利用VSP观测的(水平界面)反射波走时反演椭圆各向异性速度；刘玉柱等[12]实施了VTI介质中的多参数联合走时反演；刘瑞合等[13]讨论了VTI介质中各向异性参数的反演策略；Wang等[14]在介质的对称轴与界面垂直的假设条件下，提出了基于网格射线追踪算法的TTI介质走时层析成像方法；Wang[15]利用修改的牛顿迭代法求解稳定非线性问题，采用最小变化约束策略避免每次迭代时计算雅可比矩阵元素，并在井间观测方式下反演了各向异性介质模型参数；Waheed等[16]提出了针对VTI介质的伴随状态法初至波走时层析成像法，利用梯度法求解非线性问题，避免了直接计算偏导数矩阵。然而，纵观现有的各向异性反演方法，大部分是研究弱各向异性介质或简单的介质模型，如一维正交各向异性介质模型、椭圆各向异性介质及VTI介质。

为实现一般各向异性介质中的非线性反演，Zhou等[17-19]讨论了走时关于各向异性弹性参数的偏导数计算问题，采用初至波(qP、qSV、qSH)走时分别反演了各向异性参数。对于地震层析成像而言分辨率至关重要，而仅用初至波资料在一定程度上限制了成像深度和空间分辨率。因此，如何提高成像分辨率是一个亟待解决的问题。黄光南等[20]利用反射qP波走时实现了二维TI介质层析成像，认为利用多波走时联合有助于改善反演结果。目前，已经实现了二维TI介质多震相旅行时层析成像[21]。

本文研究了三维TI介质中多波走时层析成像技术。将各向同性介质中的分区多步不规则最短路径算法[22，23]推广至三维TI介质中，实现多次透射、反射以及转换波的追踪计算，采用各向同性地震多震相走时层析成像技术[24，25]，结合有关群速度、相速度偏导数计算[17]，通过联合多波走时同时反演5个各向异性弹性参数。数值模拟结果表明，多波走时资料联合反演可有效提高成像的空间分辨率。

2 多波射线追踪计算方法简介

将各向同性介质中分区多步不规则最短路径算法推广至各向异性介质中进行多次波的射线追踪计算，其中在单元节点上采用Zhou等[26]推导的公式计算群速度。在分区多步计算中的分区是指在模型参数化时，可根据反射界面的具体分布情况将模型分成若干个独立的计算区域(相邻的区域由界面连接)。多步计算则是根据所要计算地震震相的种类从震源所在的区域逐区(步)进行计算。具体来说，自震源所在的单元开始进行扫描，按照一定的步长进行波前扩展，若当前区域所有网格单元内节点都计算完毕后，波前就停留在该区的速度离散界面上。若追踪下行透射波(或透射转换波)，则从速度界面中走时最小的点开始(根据惠更斯原理，视为次级震源)，继续新区内的波前扩展和射线追踪；若追踪上行反射波(或反射转换波)，则从速度界面中走时最小的点开始，继续原区内的波前扩展和射线追踪；若存在转换波，则调用相应的速度模型。重复上述步骤则可实现多震相地震射线追踪计算。

三维TI介质一般可以用7个物理量m∈{a11,a13,a33,a44,a66,θ0,φ0}(5个弹性模量参数和对称轴角度)来表述。为演示该算法的射线追踪能力，给出了一个两层TI介质中多波地震射线路径追踪结果(图1)。可见，无论是一次反射波(图1a)还是多次反射转换波(图1b)，该算法均能进行追踪计算。

3 反演算法

一般而言，非线性多波(含反射或折射波)走时联合同时反演，可归结为带约束的阻尼最小二乘最优化问题，其反演问题的二范数目标函数为

(1)

式(1)的法方程为[27]

(2)

式(2)可用循环迭代的共轭梯度法进行求解[28]。目前的关键问题是如何求解具有偏导性质的雅克比矩阵A中的元素。

雅克比矩阵元素是走时关于弹性参数的导数，采用Zhou等[19]推导的一阶走时扰动方程求解。该公式不同于传统的基于特征向量的走时扰动方程[29，30]，不需要计算特征向量，能够避免两个准剪切波的奇点问题。

图1 两层各向异性介质中的多波地震射线 路径追踪结果

(a)一次透射、反射波射线路径(黑线为P1P1(上角表示上行波，下角表示下行波)，绿线为P1P2P2P1)；(b)多次透射、反射转换波射线路径(黑线为P1SV1P1，绿线为P1SV1P2SV2P1)

第一层弹性模量参数分别为：a11=9.09(km/s)2；a13=2.98(km/s)2；a33=7.54(km/s)2；a44=2.28(km/s)2；a66=3.84(km/s)2。对称轴角度θ0=0°、φ0=0°。第二层弹性模量参数分别为：a11=20.31(km/s)2；a13=9.58(km/s)2；a33=22.29(km/s)2；a44=8.35(km/s)2；a66=11.36(km/s)2。对称轴角度θ0=45°、φ0=315°

沿着射线路径R的走时一阶扰动方程可以用相速度的导数表示为

(3)

式中：m为模型参数，m∈{a11(x),a13(x),a33(x),a44(x),a66(x),θ0(x),φ0(x)}，x=(x1,x2,x3)为射线路径坐标； ds=|dx|。

(4)

根据式(4)就可以得到雅克比矩阵元素的计算公式

(5)

值得注意的是，无论是射线路径和走时的追踪计算，还是反演方程的建立，或是雅克比矩阵的计算，文中方法都未假设介质为弱各向异性，并且每次迭代都重新追踪射线路径和走时，并计算雅克比矩阵元素。因此，该方法适用于TI介质。

4 数值模拟实例及分析

4.1 断层模型

对断层模型(图3a)模拟井间地震及VSP成像，炮检排列如图2a所示。在反演模拟实验中，采用分区多步最短路径法追踪计算断层模型(图3a)的6类波型(初至波qP、qSV和qSH及模型底界面的一次反射波qPqP、qSVqSV和qSHqSH)的走时作为反演的观测走时。图2为射线路径及走时示意图。由图可见，qP速度最快，qSV和qSH波分离，符合各向异性介质中波的传播特点。将断层上盘介质弹性参数作为初始模型反演断层的弹性模量参数。

图2 射线路径及走时示意图 (a)初至波qP、qSV、qSH射线路径； (b)反射波qPqP、qSVqSV、qSHqSH射线路径； (c)初至波qP、qSV、qSH走时曲线； (d)反射波qPqP、qSVqSV、qSHqSH走时曲线 将121个检波器均匀放置于地表(间距为5m)，40个检波器、40个炮点均置于4口钻井中(模型直立的四条棱，每个钻孔 放置10个检波器，间隔为5m)，且炮点与检波器位置相同，当其中某口井作为激发井时，则另外3口井作为接收井

TI介质的5个弹性参数被分成两部分，一是与qP、qSV波相关的4个弹性参数(a11、a13、a33和a44)，二是与qSH波相关的2个弹性参数(a44和a66)。当采用某一种波反演时，只能得到与之对应的弹性参数。图3为初至qP波反演弹性参数的剖面结果。从整体上分析，y=0剖面反演结果(图3b)明显好于y=15m剖面(图3c)，这主要是因为y=0剖面上既有震源、又有接收点，射线平面即为y=0剖面，该剖面上弹性参数的反演退化为二维反演，射线覆盖相对较好。其中a44的反演效果最好，而其他三个弹性参数(a11、a13、a33)虽然恢复了断层形态和位置，但是未完全反演出断层下盘介质的弹性参数，说明qP波对a44最敏感。

当联合初至波qP、qSV和qSH走时反演弹性参数时，能够反演5个弹性参数(图4)。对比初至qP波反演弹性参数的剖面结果(图3)与初至波qP、qSV和qSH联合反演弹性参数的剖面结果(图4)发现：后者的y=0剖面的5个弹性参数反演结果(图4b)几乎与真实模型(图4a)完全一致，无论断层的起伏、形态和位置，还是断层上、下盘介质的弹性参数值均得到了很好的重建；后者的y=15m剖面的5个弹性参数的反演效果不是很理想(图4c)，但较前者的y=15m剖面(图3c)有所改善，依稀可见断层的形态和位置。

图5为初至波qP、qSV、qSH和反射波qPqP、qSVqSV、qSHqSH联合反演弹性参数的剖面结果。由图可见，当加入反射波后，成像分辨率显著提高，特别是分辨率较差的剖面(图4c)的断层形态和位置也得到很好的重建(图5c)。这主要是因为反射波qPqP、qSVqSV、qSHqSH的射线路径(图2b)与初至波qP、qSV和qSH的射线路径(图2a)存在较大区别，当加入反射波之后，射线对模型的覆盖率及角度覆盖明显增加，进而成像分辨率显著提高。因此，联合多波提高走时成像分辨率，实际上是通过增加有效的射线覆盖和角度覆盖实现的。

4.2 高、低速异常体模型

图6为真实各向异性介质模型，各向异性弹性模量参数及相应的Thomsen参数如表1所示。模拟井间地震及VSP成像时的观测系统同断层模型(图2a)。实验中采用6类波型的走时信息，即初至波qP、qSV和qSH与模型底界面的一次反射波qPqP、 qSVqSV和qSHqSH。将背景介质弹性参数作为初始模型反演高、低速异常体的弹性模量参数。

图3 初至波qP反演弹性参数的剖面结果 (a)真实剖面；(b)y=0剖面；(c)y=15m剖面

断层模型尺度为50m×50m×50m，断层上盘介质的弹性模量参数为a11=9.09(km/s)2、a13=2.981(km/s)2、a33=7.535(km/s)2、a44=2.274(km/s)2、a66=3.843(km/s)2,对称轴角度θ0=45°,φ0=0°。断层下盘介质的弹性模量参数为a11=13.84(km/s)2、a13=4.245(km/s)2、a33=11.34(km/s)2、a44=3.345(km/s)2、a66=5.051(km/s)2，对称轴角度为θ0=45°、φ0=0°

表1 高、低速异常体弹性参数及相应的Thomsen参数表

图7为初至波qP弹性参数反演结果。由图可见，反演的弹性参数能大致反映高、低速异常体的位置和形态，且高速异常体对应的弹性参数的反演效果明显好于低速异常体，这主要是因为高速异常体对地震射线具有汇集作用； 在反演的4个弹性参数中，a44的反演分辨率最高，在10(图7a)、15m切片(图7b)和15m剖面(图7d)上的高速异常体的形态和幅度恢复得较好，这说明a44对qP波群(相)速度起决定作用； 随着深度的增加，成像分辨率越来越低，如20m切片(图7c)只能反演异常位置，不能较好地恢复异常形状和幅度，不同位置的剖面图(图7d、图7e)也证实了这种现象，即无论是高速异常体还是低速异常体，均未较好地重建下部的异常轮廓。

图4 初至波qP、qSV和qSH联合反演弹性参数的剖面结果 (a)真实剖面； (b)y=0剖面； (c)y=15m剖面

为了有效提高走时成像分辨率，进行了初至波qP、qSV和qSH联合反演成像(图8)及初至波qP、qSV和qSH与反射波qPqP、qSVqSV、qSH qSH联合反演成像(图9)。分析初至波qP反演弹性参数的结果(图7)和初至波qP、qSV和qSH联合反演弹性参数的结果(图8)发现，后者不仅能反演5个弹性参数(a11、a13、a33、a44和a66)，而且反演分辨率明显提高，很好地重建了高、低速异常的形态、幅度及深度(只是20m深度切片的结果差些)，特别是也能大致看出前者中分辨率较差的a33(图7)的异常形态和位置。

分析初至波qP、qSV、qSH和反射波qPqP、qSVqSV、qSH qSH联合反演弹性参数的结果(图9)可知，反演的5个弹性参数均很好地重建了高、低速异常体，所恢复的高、低速异常体的形状、位置、幅度接近于真实模型，异常边界也清晰可见。因此，联合多波走时资料可以有效提高成像分辨率。

模型尺度为50m×50m×50m，模型中包含两个速度异常体(红色为高速异常体，蓝色为低速异常体)，大小均为10m×10m×10m，深度范围为10～20m，背景和异常体均为TI介质，对称轴角度都为(45°,0°)

图7 初至波qP弹性参数反演结果 (a)z=10m切片； (b)z=15m切片； (c)z=20m切片； (d)y=15m剖面； (e)y=35m剖面

图8 初至波qP、qSV和qSH弹性参数联合反演结果 (a)z=10m切片； (b)z=15m切片； (c)z=20m切片； (d)y=15m剖面； (e)y=35m剖面

图9 初至波qP、qSV、qSH与反射波qPqP、qSVqSV、qSHqSH弹性参数联合反演结果 (a)z=10m切片； (b)z=15m切片； (c)z=20m切片； (d)y=15m剖面； (e)y=35m剖面

5 结束语

地球内部地震各向异性介质广泛存在，前人在研究地球内部地震各向异性时都假设介质为弱各向异性，且多为VTI或HTI介质，从而简化了算法，并根据扰动理论更新模型参数，即不需要每次迭代都进行射线追踪及雅克比矩阵元素的计算。但是这种处理忽略了介质的横向不均匀性，结果往往很难与实际情况吻合。为了进行TI介质的多波走时层析成像，本文利用一阶走时扰动方程计算雅克比矩阵元素，并结合分区多步不规则最短路径射线追踪算法和共轭梯度法求解带约束的阻尼最小二乘问题，形成了一套可用于TI介质中多波走时联合反演各向异性弹性参数的技术。数值模拟实验证明该方法不仅能有效地反演各向异性弹性参数，而且还能提高成像分辨率，是一种实用的反演算法。

[1] Hirahara K and Yuzo I.Traveltime inversion for three-dimensional P-wave velocity anisotropy.Journal of Physics of the Earth,1984,32(3):197-218.

[2] Jech J.Three-dimensional inversion problem for inhomogeneous transversely isotropic media.Studia Geophysica et Geodetica,1988,32(2):136-143.

[3] Chapman C H and Pratt G R.Traveltime tomography in anisotropic media—Ⅰ.Theory.Geophysical Journal International,1992,109(1):1-19.

[4] 周辉,何樵登.非线性各向异性波形反演方法.石油地球物理勘探,1995,30(6):725-735. Zhou Hui and He Qiaodeng.Nonlinear waveform inversion in anisotropic medium.OGP,1995,30(6):725-735.

[5] 张文生,何樵登.利用VSP和跨孔走时反演横向各向同性介质中的弹性常数.石油地球物理勘探,1997,32(3):345-356. Zhang Wensheng,He Qiaodeng.Deriving the elastic constants of transversal isotropic medium from VSP and crosshole travel times.OGP,1997,32(3):345-356.

[6] 张秉铭,张中杰.一种新的地层弹性参数直接反演方法.地震学报,2000,22(6):654-660. Zhang Bingming,Zhang Zhongjie.A new way of directly elastic parameters inversion.Acta Seismologica Sinica,2000,22(6):654-660.

[7] Zheng X Y and Psenick I.Local determination of weak anisotropy parameters from qP-wave slowness and particle motion measurements.Pure and Applied Geophysics,2002,159(7-8):1881-1905.

[8] 阮爱国,王椿镛.上地幔各向异性的反演方法.西北地震学报,2002,24(2):104-112. Ruan Aiguo,Wang Chunyong.The inversion methods of the upper mantle anisotropy.Northwestern Seismological Journal,2002,24(2):104-112.

[9] 郭飚.非均匀各向异性介质的地震P波走时层析成像研究[学位论文].北京:中国地震局地震研究所,2010.

[10] 傅旦丹,何樵登,刘一峰等.一种新的全局最优化反演计算方法.石油地球物理勘探,2000,35(4):536-542. Fu Dandan,He Qiaodeng,Liu Yifeng et al.A novel global-optimized inversion method.OGP,2000,35(4):536-542.

[11] 陈竟一,张中杰.利用有偏VSP反射波旅行时反演椭圆各向异性速度—模拟退火方法.中国地球物理学会第二十届年会论文集,2004,501.

[12] 刘玉柱,王光银,董良国等.VTI介质多参数联合走时层析成像方法.地球物理学报,2014,57(10):3402-3410. Liu Yuzhu,Wang Guangyin,Dong Liangguo et al.Joint inversion of VTI parameters using nonlinear traveltime tomography.Chinese Journal of Geophy-sics,2014,57(10):3402-3410.

[13] 刘瑞合,赵金玉,印兴耀等.VTI介质各向异性参数层析反演策略与应用.石油地球物理勘探,2017,52(3):484-490. Liu Ruihe,Zhao Jinyu,Yin Xingyao et al.Strategy of anisotropic parameter tomography inversion in VTI medium.OGP,2017,52(3):484-490.

[14] Wang X X,Tsvankin I.Ray-based gridded tomography for tilted transversely isotropic media.Geophysics,2013,78(1):C11-C23.

[15] Wang Y H.Seismic ray tracing in anisotropic media:A modified Newton algorithm for solving highly nonlinear systems.Geophysics,2014,79(1):T1-T7.

[16] Waheed B U,Flagg G and Yarman C E.First-arrival traveltime tomography for anisotropic media using the adjoint-state method.Geophysics,2016,81(4):R147-R155.

[17] Zhou B and Greenhalgh S A.Analytic expressions for the velocity sensitivity to the elastic moduli for the most general anisotropic media.Geophysical Prospecting,2005,53(4):619-641.

[18] Zhou B and Greenhalgh S A.Nonlinear traveltime inversion scheme for crosshole seismic tomography in titled transversly isotropic media.Geophysics,2008,73(4):D17-D33.

[19] Zhou B and Greenhalgh S A.Nonlinear traveltime inversion for 3D seismic tomography in strongly anisotropic media.Geophysical Journal International,2008,172(1):383-394.

[20] 黄光南,Zhou Bing,邓居智等.各向异性TI介质qP反射波走时层析成像.地球物理学报,2015,58(6):2035-2045. Huang Guangnan,Zhou Bing,Deng Juzhi et al.Tra-veltime tomography of qP reflection waves in anisotropic TI media.Chinese Journal of Geophysics,2015,58(6):2035-2045.

[21] 黄国娇,白超英,钱卫.一般TI介质中多震相旅行时层析成像——以井间成像为例.石油地球物理勘探,2016,51(1):115-126. Huang Guojiao,Bai Chaoying,Qian Wei.Multi-phase traveltime tomography in general anisotropic TI media:an example of crosshole tomography.OGP,2016,51(1)：115-126.

[22] Bai C Y,Tang X P and Zhao R.2D/3D multiple transmitted,converted and reflected arrivals in complex layered media with the modified shortest path method.Geophysical Journal International,2009,179(1):201-214.

[23] Bai C Y,Huang G J and Zhao R.2D/3D irregular shortest-path raytracing for multiple arrivals and its applications.Geophysical Journal International,2010,183(3):1596-1612.

[24] 黄国娇,白超英.二维复杂层状介质中地震多波走时联合反演成像.地球物理学报,2010,53(12):2972-2981. Huang Guojiao,Bai Chaoying.Simultaneous inversion with multiple traveltimes within 2-D complex layered media.Chinese Journal of Geophysics,2010,53(12):2972-2981.

[25] 白超英,黄国娇,李忠生.三维复杂层状介质中多震相走时联合反演成像.地球物理学报,2011,54(1):182-192. Bai Chaoying,Huang Guojiao,Li Zhongsheng.Simultaneous inversion combining multiple-phase travel times within 3D complex layered media.Chinese Journal of Geophysics,2011,54(1):182-192.

[26] Zhou B and Greenhalgh S A.Raypath and traveltime computations for 2D transversely isotropic media with dipping symmetry axes.Exploration Geophysics,2006,37(2):150-159.

[27] Menke W.Geophysical Data Analysis:Discrete Inverse Theory.Academic Press Inc,San Diego,California,1984,40-126.

[28] Zhou B,Greenhalgh S A,Sinadinovski C.Iterative algorithm for the damped minimum norm,least-squares and constrained problem in seismic tomography.Exploration Geophysics,1992,23(3):497-505.

[29] Crampin S.Effective anisotropic elastic constants for wave-propagation through cracked solids.Geophysical Journal of the Royal Astronomical Society,1984,76(1):135-145.