0 引言

地下介质广泛存在各向异性，在复杂各向异性地区采用传统的各向同性介质成像算法，会引起绕射波收敛不完全，反射波归位不准确等问题，进而会导致成像精度不准确，影响储层预测的可靠性(段新意等，2014；吴国忱，2006；刘强等，2016).近年来，随着地震勘探技术的进步，生产中采集到越来越多的多分量地震资料，同时勘探目标逐渐转向复杂各向异性探区，采用纵波和转换波联合勘探能够降低勘探多解性，提高构造勘探、储层预测和油气藏描述的精度(黄绪德和杨文霞，2008；徐少波等，2014；毕丽飞等，2015；Han et al., 2014；张凯等，2015；刘洋和魏修成，2003).然而，目前对复杂介质转换波成像的研究还比较少，生产中多采用纵波成像方法来处理资料，因而有必要发展一种适应于复杂各向异性介质的转换波成像方法，从而为进一步提高地震资料解释的精度，解决复杂油气藏勘探的问题提供帮助.

20世纪70年代学者们提出了转换波勘探的理论，20世纪80年代到90年代初开始出现转换波资料处理的基本流程(Stewart et al., 1999).在速度横向变化比较大复杂区域，共偏移距域和共炮域共成像点道集中存在一些运动学和动力学假象(Xu et al., 2001; Xu and Parra, 1998).为了消除这些假象，学者开始研究角度域成像方法.目前，角度域成像可以分为波动方程类和射线类两大类.第一类波动方程类角度域成像由De Bruin(1990)、Xie和Wu(2002)、Prucha(2000)、许士勇(2002)等提出，后经过Sava和Fomel(2003)、Fomel(2004)等的发展，逐渐提高了算法的适用性和稳定性，并将其应用到转换波成像中，然而该类算法实现过程中计算比较复杂，计算效率不够理想，生产中具有一定局限性.第二类射线类角度域成像方法由Xu和Parra(1998)、Brandsberg-Dahl(1999)等提出, 学者们多研究角度域Kirchhoff偏移方法，该类方法可直接利用射线追踪中的传播角度信息提取角度域道集，具有较高的计算效率，但Kirchhoff偏移具有一定的局限性，在复杂区块不能得到理想的成像效果.针对波动方程类偏移和普通Kirchhoff类偏移存在的不足，学者们提出了高斯束偏移成像方法，该方法是Kirchhoff射线类方法的改进，兼具Kirchhoff偏移和波动方程类偏移方法的优势，不仅具有较好的成像精度，还具有较高的计算效率(Hill, 1990, 2001; Gray, 2005).李振春等将高斯束偏移逐步推广到角度域，提出了一种高斯束偏移角道集的提取方法(李振春等，2010；岳玉波，2011；秦宁等，2015). 段鹏飞等(2013)在Zhu等(2007)的研究基础上研究了各向异性介质局部角度域高斯束偏移成像方法.随后，Han等(2014)在Alkhalifah(1995)研究的声波各向异性介质追踪算法的基础上，发展了基于弹性参数的各向异性介质转换波高斯束偏移成像方法，然而，该算法实现复杂，存在计算效率不够理想等问题.

本文在Zhu等(2007)的研究基础上，推导了各向异性VTI介质弹性多波追踪方程，发展了转换波射线追踪算法.并将研究的转换波射线追踪应用到高斯束偏移中，实现了各向异性介质角度域转换波偏移成像方法.模型试算的结果表明：同传统的基于弹性参数的各向异性追踪算法和各向同性追踪算法相比，本文研究的方法在保证成像精度的前提下，具有更高的计算效率，本文研究的方法可以对复杂构造有效成像，提取的角道集可以为后期的偏移速度分析提供依据.

1 基本原理 1.1 转换波高斯束偏移成像的基本原理 1.1.1 高斯束波场延拓公式

假设U_m^ν(x, x₀, ω)为x点处接收到的由x₀处ν型波(P波或S波)震源引起位移矢量，利用弹性动力学高斯束，可将其通过一系列x₀点以不同出射角出射的高斯束叠加积分表示(Červený and Pšenčík，1983)：

(1)

式中，p_x(x₀)，p_z(x₀)分别为高斯束在x₀处射线参数的水平和垂直分量；为计算点x处ν型波的高斯束位移矢量；ω为频率；Ψ^ν为权系数，有：

(2)

式中，i为虚数单位；ω_r为高斯束的参考频率；ω₀为高斯束的初始宽度；v_ν(x₀)为x₀处ν型波的相速度；ρ(x₀)为x₀处介质密度.

本文研究的各向异性介质中，将震源位移波场通过高斯束位移矢量表示，可得到正向延拓后的波场公式：

(3)

式中，v_p(x_s)为震源x_s处qp波相速度；为震源x_s处qp波的高斯束位移矢量.

根据Hill等(1990, 2001)的研究，高斯束的初始波前为平面，据此可将地震记录分解为不同波型、不同初始方向的局部平面波，利用相应的高斯束位移矢量波场进行波场延拓，经过相移校正和高斯窗分解(岳玉波，2011；黄建平等，2015)，令ρ(x_r)≈ρ(x_L)，得到x_L处出射的qsv波的反向延拓的位移公式：

(4)

式中，ΔL为其水平间隔；N_L为高斯窗的数目；x_r为接收点位置；x_L为束中心位置；p_x^qsv(x_L)，p_z^qsv(x_L)分别为x_L处qsv波在水平和垂直方向上的慢度矢量的分量；D^qsv(x_L, p_x^qsv, ω)为以x_L为中心的高斯窗内地震波场的倾斜叠加：

(5)

1.1.2 各向异性介质qpqsv波成像公式

根据Claerbout(1971)成像条件和岳玉波(2011)给出的角度域成像公式，结合式(4)所示的反向延拓的波场公式，可以得到qpqsv波互相关成像公式：

(6)

式中，I_qpqsv为qpqsv波单炮成像值；*表示复值共轭；v_p(x_s)表示x_s处qp波的相速度.

由于转换波偏振方向不同，成像时会出现极性反转现象.本文根据qp波入射角度的正负，引入符号函数对其进行校正.图 1中给出了各向异性介质中qpqsv转换波的传播过程，入射qp波具有不同符号的入射角，导致R₁点产生的qpqsv转换波位移分量为正水平方向，而R₂点产生的qpqsv转换波位移分量为负水平方向，因而使得R₁，R₂点记录到的x分量地震记录极性相反，因此，可以根据qp波入射角α的正负直接对成像结果进行校正(岳玉波，2011)，得到最终的qpqsv转换波成像公式为

(7)

1.2 各向异性介质qpqsv转换波追踪算法

实现各向异性qpqsv转换波高斯束偏移的关键在于运动学和动力学射线追踪.本文在Zhu等(2007)研究的各向异性介质qp波射线追踪算法的基础上，通过推导各向异性介质转换波运动学和动力学射线追踪方程，发展了各向异性介质qpqsv转换波射线追踪算法.各向异性介质qpqsv转换波与qpqp波追踪算法的区别在于，对qpqsv转换波进行射线追踪时，震源处需使用qp波射线追踪，而检波点处需使用qsv波射线追踪.

1.2.1 运动学射线追踪

Červený(1972)推导了一种基于弹性参数的各向异性介质运动学射线追踪方程，然而其计算过程中需要求解Christoffel方程特征值问题，计算过程复杂.为此Zhu等(2007)推导了一种基于相速度的运动学射线追踪方程.本文将Zhu等(2007)研究的各向异性介质qpqp波运动学射线追踪算法推广到各向异性介质qpqsv波，得到如下公式：

(8)

(9)

式中，V_Pi, V_Si分别为qp波和qsv波群速度(i=x, z)，Tsvankin(2001)的研究中给出了计算群速度的近似表达式；v_P和v_S分别为qp波和qsv波的相速度，各向异性VTI介质中，有：

(10)

(11)

式中，v_P0和v_S0分别为P波和S波垂直速度；ε和δ为Thomosen(1986)表示的各向异性参数；σ= 为相速度角.

运动学射线追踪得到qp波和qsv波中心射线的走时和路径.

1.2.2 动力学射线追踪

Hanyga(1988)推导了一种基于弹性参数的各向异性介质动力学射线追踪方程，其求解过程更为复杂.在Zhu等(2007)的研究基础上，可得到射线中心坐标系(y_M, y_N, τ)下各向异性动力学射线追踪方程：

(12)

(13)

式中, Q和P为动力学射线追踪参数，相关系数由下式给出：

(14)

(15)

式中, y_M和y_N为射线中心坐标系下的坐标；q_M=∂τ/∂y_M；V_PN和V_SN分别为射线中心坐标系中qp波和qsv波群速度矢量V的分量.

动力学射线追踪得到复值的动力学射线参数，从而计算出射线的振幅和波前.

1.3 二维转换波角度域道集提取

根据岳玉波(2011)的研究可知，Hill(2001)所选择的高斯束的初始值，可保证高斯束的走时和振幅在一个波场范围内缓慢变化，因而可以直接利用高斯束的实值走时信息来计算传播角度.

对于任意射线Ω，建立二维射线中心坐标系(如图 2)，坐标系的基矢量分别为同射线Ω相切的单位切向量t和同Ω垂直并指向Ω同一侧的单位法向量n₀，R点处高斯束的实值走时满足如下公式：

(16)

式中，T(R)、T(Q)分别为R和Q点处高斯束的走时；M(Q)为Q点处走时的二阶偏导数.

对上式求关于x方向的偏导数，可得

(17)

式中，n为R在Q点为原心的射线中心坐标系中法向量方向的坐标，其满足下式：

(18)

式中，t_x，t_z为Q点射线单位切向量t沿直角坐标轴x，z方向的分量；x_Q，z_Q为Q点的坐标.将式(18)代入式(17)，且，可得

(19)

同理可知：

(20)

最终可得到网格点的传播角度(与正z轴方向的夹角)为

(21)

由式(21)可以计算出震源以及接收点高斯束的传播角度，进而求得偏移张角，这样便可以将成像值叠加到ADCIGs对应的位置.

利用上述原理直接提取qpqsv转换波ADCIGs，其关键是计算地下成像点处反射角.

如图 3，假设在地下成像点处，qp波和qsv波高斯束同正z轴方向的夹角分别为β和α，n为界面法线，qp波和qsv波反射角分别为θ_P，θ_S.则由图中的几何关系可知

(22)

根据斯奈尔定律可得：

(23)

令φ=β-α，将式(22)代入式(23)，化简可得qsv波的反射角：

(24)

在求得qsv波的反射角后，便可以进行转换波ADCIGs的提取.

2 模型试算

为了说明本文研究方法的正确性和有效性，本文通过各向异性VTI介质断块模型和各向异性VTI介质复杂构造模型进行测试.

2.1 各向异性VTI介质断块模型

为了验证本文方法的正确性，本文设计了一个各向异性VTI介质断块模型进行测试.图 4为各向异性VTI介质断块模型速度场及参数场，其大小为18010 m×3000 m.本文采用各向异性介质高斯束正演模拟方法进行正演，震源子波采用雷克子波，主频25 Hz，总炮数为281，炮间隔为50，检波器位于地表，每炮201道接收，道间距为20.总记录时间为3 s，采样间隔大小为1 ms.图 5为各向异性高斯束正演模拟得到的地震记录，图 5a为x分量，图 5b为z分量.本文用各向同性介质转换波高斯束偏移方法、基于弹性参数的各向异性转换波偏移方法(Han et al., 2014)与本文研究方法进行了对比测试，试算结果如图 6所示.图 6a为运用各向同性介质算法得到的ps转换波偏移成像的结果，图 6b为运用基于弹性参数的方法得到的qpqsv转换波偏移成像的结果，图 6c为运用本文研究的方法得到的qpqsv转换波偏移成像的结果.

由图 6可知，与各向同性介质转换波高斯束偏移方法相比，运用基于弹性参数的各向异性VTI介质转换波高斯束偏移方法和本文研究的方法得到的成像结果，信噪比高，同相轴能量聚焦性好，构造位置准确.而运用各向同性方法得到的成像结果，由于不考虑各向异性因素的影响，导致同相轴能量聚焦性不好(箭头所指部分)，同相轴两侧出现上翘的假象(方框所示区域).各向异性VTI介质断块模型的试算结果说明了本文研究方法的正确性和有效性.

为了更好地说明本文的追踪算法在计算效率方面的优势，对各向同性高斯束偏移方法、基于弹性参数的各向异性高斯束偏移方法以及本文方法的计算效率进行测试分析.表 1中给出了各向同性方法、基于弹性参数的偏移方法以及本文研究方法之间的计算效率对比.由此可见，本文研究的方法在保证成像精度的前提下，在计算效率方面具有一定优势.

2.2 各向异性VTI介质复杂构造模型

为了验证本文方法对复杂各向异性介质的适应性，本文设计了一个各向异性VTI介质复杂构造模型进行测试.图 7为各向异性VTI介质复杂构造模型速度场及参数场，模型大小为12010 m×3500 m.该复杂构造模型中存在各向异性介质平层、逆冲断层、洼陷等构造.采用各向异性高斯束正演模拟方法进行正演，震源子波采用雷克子波，主频为25 Hz，总炮数为161，炮间隔为50，检波器位于地表，每炮201道接收，道间距为20.总记录时间为3.5 s，采样间隔大小为1 ms.图 8为地震正演模拟得到的地震记录，图 8a为x分量，图 8b为z分量.本文运用各向同性介质转换波高斯束偏移成像方法、基于弹性参数的各向异性介质转换波高斯束偏移成像方法和本文中的方法分别进行了模型测试，结果如图 9所示，图 9a为各向同性介质转换波高斯束偏移得到的ps转换波成像结果，图 9b为基于弹性参数的高斯束偏移得到的qpqsv转换波成像结果，图 9c为本文研究的方法得到的qpqsv转换波成像结果.

由图 9可知，运用各向同性转换波高斯束偏移方法处理各向异性数据，导致同相轴位置不准确(图中箭头所指部分)，并且逆冲断层以及深部洼陷同相轴能量不聚焦，整体成像质量不够理想.而运用基于弹性参数的各向异性转换波成像方法(Han et al., 2014)和本文研究方法得到的结果，同相轴归位准确，绕射波得到很好的收敛，整体成像质量明显提高.同时，对比图 9b，9c可以发现，二者在成像位置上都比较准确，但受动力学射线追踪的影响，采用基于弹性参数的各向异性转换波偏移方法(Han et al., 2014)得到的结果在部分区域能量聚焦性不理想(图中方框所指区域).

为了进一步验证本文研究方法的正确性和有效性，我们在成像过程中根据反射张角提取了角度范围为0°~15°的角道集(图 10所示).图 10给出了CDP为601时qpqsv转换波偏移角道集结果，图 10a为运用各向同性高斯束偏移方法得到的角道集，图 10b为运用基于弹性参数的各向异性高斯束偏移方法得到的角道集，图 10c为运用本文方法得到的角道集.对比可以发现，运用各向同性方法得到的角道集同相轴没有拉平，而运用基于弹性参数的各向异性方法和本文研究方法得到的角道集同相轴较为平直.角道集的结果说明了本文研究方法的正确性和有效性.

综上所述，本文中研究的各向异性介质角度域转换波高斯束偏移成像方法，更能够对复杂各向异性介质地质构造有效成像，模型试算的结果说明了本文研究方法的正确性和有效性.

3 结论

本文通过修改各向异性介质运动学和动力学射线追踪方程，发展了各向异性介质转换波射线追踪算法，并将其应用到转换波高斯束偏移中，实现了一种各向异性VTI介质角度域转换波高斯束偏移成像方法.通过引入符号函数，较好的解决了转换波成像中的极性校正问题.该方法可以利用横波和纵波波场信息，解决复杂各向异性地下构造的成像问题，从而为提高储层预测的精度，解决复杂油气藏勘探的问题提供帮助.另外，本文中的角度域成像方法提取的角道集可以用于后期的偏移速度分析.

各向异性介质断块模型测试结果表明，与各向同性介质转换波高斯束偏移方法相比，本文中的方法可以对各向异性介质准确成像，与基于弹性参数的各向异性介质转换波高斯束偏移方法相比，本文研究的方法在计算效率方面具有一定优势.二维各向异性介质复杂构造模型试算结果表明，本文研究的方法可以对复杂各向异性地质构造有效成像，进一步说明了本文方法对复杂地质构造的适应性和有效性.下一步拟将本文研究的方法推广到广义TI介质，提高算法的适应性, 并进一步增强方法对更为复杂的地质构造的适应性.

致谢

感谢中国石油大学(华东)地震传播与成像(SWPI)课题组的支持与帮助！感谢审稿专家对本文给出的指导和建议！