0 引言

近年来，多波多分量地震勘探不断发展并得到了工业界越来越多的关注.多波的传播属性与地下介质的物性参数有关，可以利用纵波和横波传播的运动学和动力学属性及差异，研究岩石岩性和物性，岩石裂隙发育状况，孔隙度及流体饱和度等(Li，1997；Qian et al., 2007；栗学磊和毛伟建，2016).特别是横波信息，有助于无井区的油气检测，可以降低勘探风险.而多波多分量地震数据的偏移成像在多波多分量勘探中占据着重要地位，是能否有效获取更加准确的地下构造及岩性信息的关键.

传统的多分量数据偏移先将地震记录分解为标量的P波和S波，然后基于声波方程进行标量偏移.由于不同的波型在传播过程中会相互耦合，所以事先分离通常很难将其完全分离开.大多数的波型分离方法，如亥姆霍兹分解和偏振滤波，都是基于波传播和偏振方向等几何特性，没有考虑弹性动力学参数耦合的影响.这样会引起成像结果中不同波型之间的串扰.同样，声波方程也不能准确地描述弹性波的传播特性，进而不能准确实现弹性波的保幅偏移和反演.基于弹性波动方程的多分量偏移始于20世纪80年代，包括弹性Kirchhoff偏移(Kuo and Dai, 1984；Dai and Kuo, 1986；Berkhout and Wapenaar, 1989；Hokstad，2000；Schleicher et al., 2001)和弹性逆时偏移(Chang and McMechan, 1994；Yan and Save, 2008；Du et al., 2012, 2014)等.但它们中的大多数还是标量偏移，对PS转换波的极性翻转问题的解决还不够理想(Li et al., 2018).目前更为先进的技术是不对多波多分量数据进行事先分离，而是在偏移过程中利用矢量波场特征，将其自动分离为P波和S波.这样不仅可以减少波型分离不准确造成的串扰噪声，还能降低处理成本.多波多分量矢量波场成像的优势见栗学磊和毛伟建(2016)、Li等(2018)的论述.

以Kirchhoff偏移为代表的射线类方法具有灵活性和计算效率高的优势，但不能处理复杂地质体，进行陡倾角成像.以逆时偏移为代表的波动方程类方法成像精度高但计算量大.Hill(1990, 2001)提出的高斯束偏移方法兼顾两者优势，可以处理多路径问题，是一种快速且精度较高的成像方法，具有明显的工业应用价值.

另一方面，地下介质的黏弹性会导致地震波发生能量衰减与相位畸变，特别是在类似气云这样的高衰减区域.基于声介质或者弹性介质假设的偏移方法在对这些区域进行成像时，效果会受到影响.而补偿这类衰减效应可以提高地震偏移的成像质量.通常，该效应可以由品质因子Q来定量描述.早期的吸收补偿是直接对地震道做反Q滤波(Wang, 2002)，这种情况无法得到地震波的传播路径，而实际中地震波的能量衰减与其传播路径密切相关(陈雪, 2016).因此，更理想的做法是根据地震波的实际传播路径，在偏移成像的过程中进行补偿.波动方程类偏移方法的衰减补偿通常是通过构造黏弹性波动方程(Deng and McMechan, 2008；Zhu et al., 2014, Zhu and Sun, 2017)或者延拓算子(Mittet，2007)进行的；而对于射线类方法，可用更简便的方法，即通过定义与Q有关的衰减走时影响波形，从而进行相位和振幅补偿(Traynin et al., 2008；Xie et al., 2009, 2010；Cavalca，2013；吴娟等，2015).以上衰减补偿算法的实现，都需要依赖相应的黏弹性吸收衰减模型.黏弹性可以根据介质非弹性吸收本构关系的Boltzmann线性叠加原理，通过不同的物理模式(弹簧和阻尼元件不同方式的组合，如Maxwell体，Kelvin-Voigt体，标准线性体等)或数学模式(如幂律模型，常Q模型等)来表示(张立彬等，2011).理论上，各种模型都可以被用来实现成像过程中的衰减补偿.在地震勘探观测频率范围内，可以认为Q值大小实际不随频率而变化.本文中我们的做法便依赖于这一假设.

目前，有关多波多分量矢量波场弹性高斯束偏移方法的研究还不多，至于考虑到黏弹性的方法更是少见的.白敏等(2016)提出了基于吸收衰减补偿的多分量高斯束偏移，但他们首先把波场分解为PP波和PS波，然后分别基于声波方程进行标量方法偏移，并不是真正的基于弹性波方程的矢量波偏移成像.栗学磊和毛伟建(2016)提出了一种在偏移过程中实现分离的机制，将高斯束偏移从声介质推广至弹性介质.为了直接对矢量波进行成像，同时解决极性翻转问题，他们引入了散度/旋度互相关成像条件.但是此方法没有考虑实际地下介质的黏弹性，从而无法补偿能量衰减和校正相位畸变，这使得该方法对一些含高黏弹性地层成像效果不佳.为此，在他们工作的基础上，我们在偏移过程中沿射线路径引入品质因子Q来考虑黏弹性影响，通过复速度对走时和振幅进行修正，将高斯束偏移从弹性介质成像推广至黏弹性衰减补偿成像.

本文首先给出黏弹性高斯束的表达式，然后给出经过衰减补偿修正后的向下延拓和成像公式.同时，给出在矢量波场偏移过程中提取角度域共成像点道集的方法，以便用于成像质量控制，并为后续速度和黏弹性参数反演提供所需的数据.利用2D层状模型和洼陷模型进行方法测试，其成像结果用于验证所提出的黏弹性衰减补偿高斯束偏移方法的可行性和有效性.

1 原理与方法 1.1 黏弹性Q衰减及补偿

根据对应原理(Ben-Menahem and Singh, 2012)，波在黏弹介质中的传播可以认为是复速度的波在弹性介质中的传播.复速度可以定义为(Aki and Richards, 1980)

(1)

式中Q是品质因子，w是频率，ω_r代表参考频率，c₀是在弹性介质中的速度(可以认为是地震波频率为ω_r时对应的速度).从复速度可以推出衰减走时的表达式：

(2)

其中，τ是弹性介质的波走时.τ^*是黏弹性引起的走时变化，也称为全局吸收因子(Červený，2001 )，即

(3)

衰减可以通过公式(2)中依赖频率的复走时进行表达，该式中的后两项分别与品质因子Q引起的振幅衰减和相位畸变有关.

将公式(2)中的后两项取反号即可实现黏弹性吸收衰减补偿，我们将其定义为衰减补偿走时，即

(4)

1.2 黏弹性高斯束

高斯束是波动方程在高频傍轴近似情况下的一种特解. Červený和Pšeník (1983)推导了3D弹性非均匀介质中的高斯束表达式，栗学磊和毛伟建(2016)给出了它的归一化多波表示的统一形式

(5)

其中，ν代表地震波模式(P, S1, S2；在3D各向同性介质中，S波一般情况下为椭圆偏振，为了表示和计算方便，将其分解为S1和S2两种线性偏振横波( Červený, 2001 ))，u_GB^ν是点x处的ν波高斯束位移矢量，r和p^ν分别为高斯束的入射点和入射慢度矢量，q^T=(q₁, q₂)是垂直射线的平面坐标系.为了区分品质因子Q，使用Q_d表示动力学参数，且Q_d是一2×2的矩阵.a^ν为ν波偏振方向矢量，且有关系式

(6)

其中，M_IJ为2×2矩阵M的元素，是波前走时在(q₁, q₂)坐标系的二阶偏导数，满足M=P_dQ_d^-1.e₁，e₂和e₃≡t是3D中心射线坐标基.式(6)包含了与坐标q_J相关的附加项，这是偏振方向矢量的校正项.当x离中心射线很近时，q_J很小，此时校正项可以忽略不计.α、β分别为纵波波速和横波波速.

2D情况下弹性高斯束的表达式退化为

(7)

式中，(s, n)为二维中心射线坐标系，p_d和q_d是高斯束动力学参数.a^ν为ν波偏振方向矢量，且有关系式

(8)

矢量t和n为2D中心射线坐标基.

前面已经提到，波在黏弹介质中的传播可以认为是复速度的波在弹性介质中的传播.对于复速度的运动学射线追踪，Keers等(2001)证明了在一阶近似下，射线路径与Q无关. Červený(2001) 也指出，可以使用全局吸收因子τ^*(见公式(3))完全控制地震波的衰减.这样黏弹性射线追踪的结果可以通过用复走时取代弹性射线追踪得到的走时而获得.黏弹性动力学射线追踪需要使用复速度，从公式(1)中可看出Q对于复速度的影响在1/Q这个量级.勘探地震成像涉及到的地层Q值一般大于20，从而对速度的影响不足5%，一阶近似下认为其对动力学参数的影响可以忽略.由(6)和(8)式，我们认为非弹性衰减对各种波模式的偏振方向的影响也可以忽略.在此基础上，黏弹性高斯束可以通过只将公式(5)或(7)中的走时项τ替换成公式(2)中的τ_a表示.3D黏弹高斯束的表达式如下：

(9)

其中项代表黏弹性引起的振幅的衰减，项代表黏弹性引起的相位的畸变，二维情况类似.同理，补偿后的高斯束，只需将公式(5)或(7)中的走时项τ替换成公式(4)中的τ_Q.高斯束的走时本身就是复数值，其虚部代表波前面上远离中心射线的能量衰减.引入与Q有关的走时校正项，代表沿射线路径方向的衰减.

图 1，图 2和图 3分别给出了均匀介质中弹性(Q→∞)，Q=100和Q=20三种情况下不同波模式(P, SV和SH)的高斯束.此处SV和SH的划分与上文S1和S2的划分不同，S1和S2的划分与坐标系设置有关，属于人为行为(栗学磊和毛伟建，2016).从图中可以直观看出高斯束在黏弹介质中的衰减情况.对这三种情况，使用了相同的高斯束参数和频率，只对走时进行了替换.对于P波而言，与射线方向相同的高斯束分量称为主分量，而在波前面上的称为次分量.而对于SV波，垂直于射线方向的高斯束分量称为主分量，与射线方向相同的称为次分量.在PS转换波中真正有效的S波是SV波，S1和S2的合并才能完整代表SV波.至于SH波则不存在次分量.主分量与次分量的划分，体现了高斯束的矢量特征.图 1和图 2中高斯束的主分量与次分量在黏弹情况下均发生了衰减，且在Q值越小的情况下，衰减越快.SH波同样满足Q值越小衰减越快的规律.

1.3 含Q的高斯束偏移

栗学磊和毛伟建(2016)给出了波场向下延拓的一般表达式，在考虑衰减补偿后，数学形式如下：

(10)

其中，上标Q代表已经进行过Q补偿，u_n^{ν, Q}为ν波在x处的位移矢量, u_{GB, n}^{ν, Q}是矢量u_GB^{ν, Q}的第n个分量, a为束中心间隔，ω_l和w_l分别为高斯束参考频率和参考有效半宽度，L代表高斯束中心坐标，p为慢度矢量，*代表复共轭, D^ν是多波型分解后的标量局部平面ν波，

(11)

D_n为第n个分量的局部平面波.由于位移分量u_n包含多种波型，因此局部平面波D_n是多种波型的叠加.H_n^ν是一个3×3大小的分解矩阵，与传播和偏振方向以及模型参数有关.通过方程(11)可以自动地把多波型波场分解为P，S1和S2波.公式(10)显示，向下延拓过程中地震波场的振幅和Q有关.

当前，大多数多分量偏移方法均假设模型为自由地表模型.下面直接给出自由地表边界条件下矩阵H_n^ν的表达形式(其具体的推导过程，见栗学磊和毛伟建(2016)一文)：

(12)

其中W_n^ν是一个3×3的权重矩阵

(13)

在二维情况下，该矩阵退化为

(14)

式中，λ和μ为拉梅系数.p₁，p₂和p₃分别为z，x和y方向上的慢度分量.

黏弹性波场与弹性波场一样，为矢量波场.常规的互相关成像条件不适用于矢量波成像，并且PS转换波存在极性翻转问题.为了解决矢量波成像问题，本文利用了栗学磊和毛伟建(2016)提出的散度/旋度互相关成像条件.成像表达式如下：

(15a)

(15b)

式中，算子▽_x下标表示作用于点x，u_r^Q(x, t)和u_s^Q(x, t)分别代表检波点和炮点经过衰减补偿后的地震波场，函数Sgn(x)为一单位矢量函数，并且有表达式

(16)

其中，p_s和p_r分别为炮点射线慢度矢量和检波点射线慢度矢量，Sgn(x)为同时与正向延拓和反向延拓方向垂直的单位矢量，它有具体的物理意义，代表转换波的一种旋转方向.当Sgn(x)和波场的旋度同向时提取出正标量振幅，当Sgn(x)和波场的旋度反向时提取出负标量振幅, 由此，极性翻转在矢量波场中通过这一符号得以解决(栗学磊和毛伟建，2016).

利用公式(10)和(15)，可以实现黏弹性衰减补偿偏移成像.公式(10)到(16)显示将偏移推广至黏弹性衰减补偿后，整体的架构不发生变化，只是需要重新构造高斯束表达式.实际偏移过程中，黏弹性衰减补偿偏移与弹性偏移最大的区别在于需要对复走时进行与Q有关的校正.对于射线追踪中的每一个步长，认为该段Q值是一定的，对该段走时信息进行Q补偿校正，最终累加获得经过Q校正的总走时.由于高斯束本身就具有复走时，因此引入与衰减有关的走时对于程序实现而言并没有带来太多额外的计算量.

1.4 角度参数计算

在已知射线束走时场的情况下，可以通过构造走时场梯度和地下反射张角的关系提取地下角度域道集.在三维情况下，走时梯度方向和炮检点出发的局部平面波矢量的关系可以表示为

(17)

(18)

其中，p_s和p_r分别为炮点射线慢度矢量和检波点射线慢度矢量，t_s和t_r分别为炮点和检波点端的单向走时.

角度域的成像是通过控制地下射线的角度选定参与成像的数据采样来实现的.本文所定义的角度域主要指偏移张角，对于PP波是反射角的两倍，对于PS波是入射角和反射角之和.在某个成像点上，结合炮点射线慢度矢量p_s及检波点射线的慢度矢量p_r即可求得角度信息.偏移张角可以通过如下运算得到(段鹏飞等，2013)

(19)

在求得偏移张角后，便可将成像值归到角度域共成像点道集对应的位置.

2 数值试验

为了验证本文所提出的方法的正确性，我们设计了两个模型来进行试验.模型记录由Jan Thorbecke提供的二维有限差分波场正演开源程序计算而得到.

2.1 简单水平层状模型

模型参数如图 4所示.反射界面在深度z=1000 m处，炮点位于(x=2000 m，z=0 m)处，检波点位于x=0~4000 m，检波点间隔为10 m.震源采用主频为22.5 Hz的雷克子波.地震记录长度为2.5 s.Q_P和Q_S不随空间变化且均为常数20.图 5展示的是弹性和黏弹性单炮的正演结果.对比(a)和(c)可以看出，虽然PS波在x分量上，PP波在z分量非常明显，但它们在另一分量上也是清晰可见的.左边的(a)、(c)与右边的(b)、(d)对比，可以看出右侧的黏弹性记录反射波的同相轴能量较左侧弹性的弱，表明介质的黏滞性对地震波的吸收和衰减作用明显.图 6和图 8分别是PP波和PS波在三种情况下(参考，补偿前和补偿后)偏移之后的结果.图 6与图 8中的(b)相较于(a)，界面处的能量明显弱很多，且界面位置也有所移动，而(c)中的成像结果，与(a)中的非常接近，证明补偿效果良好.为了更好地展示，分别从图 6和图 8中的(a)、(b)、(c)中各抽取一道(横向距离2000 m)波形进行对比(见图 7和图 9).可以看出，黏弹记录偏移后得到的波形的相位比弹性的要提前，能量也要弱很多.经过Q补偿之后，相位畸变得到了很好的校正，同时能量也得到了恢复.值得注意的是：图 6中的划弧是由于测线两端数据的突然截断引起的边界效应，图 8中炮点正下方界面能量接近于零是由于垂直入射时没有转换波生成.

2.2 洼陷模型

在本例中，我们验证衰减补偿偏移算法在多炮叠加情况下的有效性，同时进行角度域的成像.速度和Q值模型如图 10所示，高衰减区域主要在楔形内.介质密度ρ基于Gardner关系式ρ=0.31·V_P^0.25确定.合成的记录共有201炮，炮间距为20 m.每一炮有401个检波点接收信号，其中检波点间隔为10 m.震源采用主频为10 Hz雷克子波.地震记录长度为6.1 s，时间采样间隔为0.001 s.图 11和图 13分别是PP波和PS波在三种情况下(参考，补偿前和补偿后)偏移之后的结果.图 11b楔形左下侧处相较于a，能量很弱，界面的延续性较差，显示出地震波在经过较长距离的高衰减区域后会影响以下地层的成像效果，进一步说明了进行黏弹性衰减补偿的必要性.图 11和图 13中洼陷的底部及最下面弧状基底，均能看出(b)中该处同相轴的能量中比(a)中的弱，且成像位置不准确，而(c)则基本恢复到与(a)中的一致.为了更好地展示，分别从图 11和图 13中的(a)、(b)、(c)中各抽取一道(横向距离2000 m)波形进行对比，得到图 12和图 14中的结果，同样可以看出黏弹记录偏移后得到的界面(特别是深度为1800 m，2800 m和3300 m处)明显发生了错位，能量也比弹性参考下的弱很多，但经过Q补偿之后，相位畸变得到了很好的校正，而能量也得到了恢复.

此外，我们还计算了经衰减补偿后的地震波在某一空间位置处(x=2000 m处)的角度域共成像点道集(ADCIGs)，如图 15和图 16所示.每一道包含了多个偏移界面在某一角度时的成像点结果.可以看到道集中的同相轴均为平直的，这也证明了该角度域成像算法的准确性.图 15中角度域成像结果的第一层在800 m处，与实际模型一致，经计算此处的临界角大致在60°，因此超过120°的部分没有意义.第一层往下，因受观测系统限制，照明角度越来越小，同相轴大体随深度增加而变短.图 15b中第三个界面处的同相轴的能量明显比(a)中的要弱很多，这是因为黏弹情况下的地震波刚经过一段高衰减区.图 16是PS波的角度域共成像点道集.横波的能量比纵波弱很多.PS波的角道集在角度很小时几乎没有能量，说明垂直入射时没有转换波生成.从图 15和图 16中同样可以看出，能量的衰减和反射界面位置的错位，在经过黏弹性补偿之后得到很好地恢复.

3 结论

在本文中，我们将黏弹理论与高斯束偏移成像方法相结合，在3D高斯束矢量波成像的框架下，通过引入与品质因子Q有关的走时校正，实现了黏弹性衰减补偿偏移.相较于波型分离之后再分别对PP波和PS波进行标量偏移的处理技术，我们的方法可以做到在处理过程中考虑波传播的矢量特征，因此对于转换波有更好的成像效果.本文利用二维层状模型和洼陷模型进行了方法测试.PP波和PS波的成像结果表明我们所提出的方法能有效补偿多波多分量地震数据的衰减效应，改进成像效果.角度域的成像，其结果同样可以验证衰减补偿方法的有效性.三维情况可以在与二维类似的偏移框架下实现，我们将在下一步工作中详细讨论.