0 引言

逆时偏移(Reverse Time Migration，RTM)具有无倾角限制、适应任意复杂速度模型的优点，是目前复杂地下构造成像的重要方法(Yoon et al., 2003; Yoon and Marfurt, 2006; Zhang et al., 2011).相较声波逆时偏移而言，弹性波逆时偏移(Elastic Reverse Time Migration，ERTM)假设地球介质中P波和S波同时传播，更符合实际地质情况，因此利用弹性特征能够获取更准确的成像结果.

与声波逆时偏移(Claerbout, 1971)的实现思路保持一致，弹性波逆时偏移处理多分量地震数据时主要包括两部分：波场构建和成像条件.早期基于笛卡尔分量的弹性波成像算法(Chang and McMechan, 1986, 1994)对波场位移分量进行矢量成像，其存在波模式串扰、物理意义不明确且伪像难以压制的问题(Yan and Sava, 2008).为了压制伪像干扰并明确成像结果物理意义，引入波场分离得到纯波模量以提取PP、PS、SP和SS成像结果(Yan and Sava, 2008).目前针对纯波模量主要存在两种构建思路：标量处理方法和弹性处理方法.其中，标量处理方法是分离得到标量P波和S波波场，并用于标量逆时偏移(Sun and McMechan, 2001; Sun et al., 2006)，该方法可利用声波逆时偏移已有成熟的方法，但忽略了矢量波场固有的弹性特征.弹性处理方法则基于弹性波动方程进行波场延拓，该方法可有效保持波场传播过程中的弹性特征和矢量特征.

其中，波场分离是弹性处理方法的关键环节.基于Helmholtz定理可以分离得到无旋P波波场和无散S波波场(Sun and McMechan, 2001)，从而实现纯波波场成像.但三维转换波成像时，存在标量P波和矢量S波成像困难的问题，因此引入标定方法将矢量S波标量化进行互相关成像(Du et al., 2014)，针对极性反转问题，引入转换波极性校正因子实现转换波极性校正(Balch and Erdemir, 1994; Du et al., 2012)，但是复杂地质情况下转换波成像效果仍不理想.另一方面，基于解耦延拓方程可以实现矢量P波和S波波场分离(马德堂和朱光明，2003；李振春等, 2007).为提高计算效率，引入辅助P波应力和S波应力对解耦延拓的实现策略进行改进(Xiao and Leaney, 2010; Du et al., 2017).与Helmholtz波场分解方法相比，解耦延拓方程计算效率虽有待提高，但是在分离波场相位和振幅保持方面优势明显，且成像结果物理意义更明确.

成像条件则是逆时偏移另一关键环节，可以实现纯波波场中属性提取.最早有激发时间成像条件(Chang and McMechan, 1986, 1994)，即射线追踪确定节点处传播时间，以此提取接收波场值作为成像结果.该成像条件计算和存储量消耗较低，但保幅性较差.以此为基础，提取接收波场和震源波场的比值作为成像结果(Nguyen and McMechan, 2012, 2014)，并成功应用到弹性波逆时偏移中(张智等，2013).该成像条件兼顾了计算效率与保幅性，但是在弹性转换波成像中应用有待商榷.目前常用互相关成像条件(Claerbout, 1971; Loewenthal et al., 1991)，通过接收波场与震源波场互相关作为成像结果，其稳定性高且充分利用所有波场信息，但深部反射能量弱且低波数噪声严重.因此，Kaelin和Guitton(2006)提出归一化互相关成像条件进行照明补偿以增强深部反射层能量.Rocha等(2016)从能量密度出发提出了能量互相关成像条件，有效压制低波数噪声的同时提取出明确表征波场能量的成像结果.但是得到结果中弹性波所有波场能量耦合在一起，难以区分纯波波场信息.

为此，我们基于能量互相关成像条件，引入Helmholtz定理和体应力构建弹性波自解耦能量互相关成像条件：基于Helmholtz原理可以将势能互相关成像条件自解耦为PP波、转换波(C波)及SS波势能成像条件；引入体应力作为P波应力以构建P波波场质点振动速度，从而构建动能解耦互相关成像条件.该自解耦成像条件既可以在振幅相位保持特征基础上获取清晰表征能量的成像结果，同时实现了弹性波波场模式分解，避免了引入额外的分离运算，提高了计算效率.

1 理论基础 1.1 基于弹性波能量密度的成像条件

对于各向同性介质，能量密度(徐芝纶，1990)定义如下：

(1)

式中，E(x, t)表示在时刻t点x处的能量密度，ρ为密度，υ_i为振动速度，σ_ij、ε_ij分别为应力、应变，λ、μ为拉梅系数，(∇·U)为波场散度，∇U为波场梯度，上标T表示矩阵的转置，其中为动能密度，为势能密度.在不施加外力情况下，t时刻点x处，即弹性体中能量不发生变化，具有唯一性和确定性(附录A).

对于震源场U(x, t)和检波场V(x, t)，定义能量内积为：

(2)

其中，U和V分别表示震源波场和接收波场，I_E为能量内积(Rocha et al., 2016).在时刻t₁，入射波和反射波同时传播到成像点x₁处，则I_E(t_i, x_i)≠0，反之，I_E(t_i, x_i)=0.能量内积I_E作为成像条件提取反射界面处的机械能，称之为能量互相关成像条件.能量团互相关成像条件由空间导数u_ij(i, j=x, y, z)和时间导数(i=x, y, z)等12个因子分量组成，其中8个空间导数表征波场的势能成分，而4个时间导数表征波场的动能成分.由于12个因子分量皆为延拓计算中的中间变量，因此实现过程中只增加存储量而不增加额外的计算量.Rocha等(2016)依据平面波理论提出，构建接收波场时间延拓的反向方向即可压制背向散射，因此将式(2)改写为背向散射衰减的成像条件：

(3)

基于能量互相关的能量互相关成像条件可以提取出弹性波能量，具有明确的物理意义，而且与基于速度互相关的常规互相关成像条件相比，能量互相关成像条件可有效压制背向散射，避免低波数噪声的干扰.与Rocha等(2016)提出的弹性波场能量成像条件不同，虽然∫_Ωu_{i, j}u_{j, i}-u_{i, i}u_{j, j}=0，但是空间Ω任意点x₁处u_{i, j}u_{j, i}-u_{i, i}u_{j, j}≠0，因此我们认为u_{i, j}u_{j, i}-u_{i, i}u_{j, j}不能省略.

1.2 弹性波自解耦能量团成像条件

实际传播过程中，弹性波是由纯纵波场(PP波)、纯横波场(SS波)和转换波场(C波)组成.其中，转换波场包括常规意义上的PS转换波和SP转换波，由于PS转换波和SP转换波势能难以分离，因此我们将PS和SP波统称为C波.由于基于弹性波能量互相关成像条件(3)得到成像结果将弹性波中所有波场模式的能量耦合在一起，难以区分PP波场、SS波场和C波场.为此，需要针对能量互相关成像条件进行波场模式分解.考虑成像条件中势能由位移场空间导数组成而动能由位移场时间导数组成，我们分别从势能(4)和动能(5)两部分进行解耦.

(4)

1.2.1 弹性波势能解耦

对弹性波势能(4)进行因式分解：

(6)

式中，U和V分别为震源波场和接收波场，(∇·U)、(∇×U)、∇U分别为相应位移场的散度、旋度和梯度，其中上标*表示伴随矩阵.

在各向同性介质中，根据Helmholtz分解定理(孙成禹和李振春，2011)，任何矢量场都可分解为无旋场和无散场，弹性波位移场中可以分解为满足∇×U_P=0的无旋的P波位移场和满足∇·U_S=0的无散的S波位移场.式(6)第一项(λ+2μ)(∇·U)(∇·V)只包含P波分量而不包括S波，即为纵波势能场(PP波)；第二项μ(∇×U)·(∇×V)只包含横波分量而不包括纵波，即为横波势能场(SS波)；根据能量守恒定律，第三项μ∇U:(∇V)^*为P波和S波耦合的势能场，因此我们统称为转换波势能场(C波).

自解耦得到纵波势能、横波势能以及转换波势能互相关成像条件：

(7)

式中，I′_{E, ε, PP}、I′_{E, ε, SS}和I′_{E, ε, C}分别为纵波势能、横波势能以及转换波势能，上标*表示伴随矩阵.自解耦势能互相关成像条件是由位移的空间导数组成，其中纵波势能由位移场的散度组成、横波势能由位移场的旋度组成、转换波势能由位移场的梯度组成.在不需要进行额外波场分离的情况下，可以实现耦合势能的自解耦，得到保幅的纯波势能场，具有明确的物理意义.其中，计算效率方面，空间导数皆为波场延拓的中间变量，不需要增加额外计算量.

1.2.2 弹性波动能解耦

弹性波位移场满足U=U_P+U_S，弹性波动能(5)可表示为：

(8)

式中，和分别为纵波位移场和横波位移场的时间导数，其中-表示求位移场时间导数的负数，即表示位移场沿时间轴传播的相反方向.其中，式(8)的第一项·(-)为纵波动能场(PP波)，第二项·(-)只包含横波分量，为横波动能场(SS波)，第三项(·(-)+·(-))为转换波动能场(C波).解耦得到动能互相关成像条件：

(9)

式中，I′_{E, k, PP}、I′_{E, k, SS}和I′_{E, k, C}分别表示纵波动能、横波动能以及转换波动能.其中，动能场由矢量位移场的时间导数组成，即速度场.为此，基于解耦延拓方程构建体应力τ^P=(λ+2μ)θ，即可得到纵波位移的时间导数, ，从总速度场中减去纵波速度场就可得到横波位移的时间导数.因此，构建得到的速度场为振幅保持的矢量场，且转换波场中无极性反转问题.在计算效率方面，虽然引入体应力等中间变量来构建纵波和横波的时间导数，但是引入的中间变量无需进入迭代，因此，只额外增加了少量计算量.

1.2.3 弹性波能量成像条件

假设弹性体形变小，无热能等其他能量消耗，弹性波场能量由动能和势能组成，根据能量守恒定律，弹性波自解耦能量互相关成像条件：

(10)

式中，I′_{E, PP}、I′_{E, SS}和I′_{E, C}分别表示PP能量场、SS能量场以及C波能量场，其中弹性波能量即为PP、SS和C场能量之和，即I′_E=I′_{E, PP}+I′_{E, SS}+I′_{E, C}.该自解耦互相关成像条件可以实现弹性波场能量的解耦，得到PP能量场、SS能量场和C波能量场，从振幅特征保持方面，虽然引入Helmholtz原理分解势能波场，但是避免对能量场进行额外改造，最终分解得到的仍是振幅保持的势能场；在计算效率方面，只额外引入少量计算.

2 模型测试

为验证方法的正确性，分别用简单平层模型、二维Marmous2模型对算法进行测试.其中，平层模型测试主要用以验证基于自解耦能量成像条件得到的PP波、C波、SS波成像结果的振幅强弱及背向散射的压制能力，并与能量成像结果进行对比，分析纯波模式成像结果的优劣；通过Marmous2模型主要测试方法对偏移噪声的压制、成像及振幅的恢复能力等，同时对比分析模式分离及常规互相关成像结果的优劣.

2.1 简单模型

设计一个大小为300×300(网格数)、网格大小为10 m×10 m二层简单模型进行测试，图 1为模型示意图.正演模拟时，震源采用主频30 Hz的雷克子波，加载爆炸震源，时间采样间隔为1 ms，接收时间长度为2.5 s.震源位于(0 m，1500 m)处，在地表以10 m间隔均匀分布的检波点进行接收，采用有限差分、时间上二阶精度、空间上十阶精度进行单炮数值模拟.

图 2为单炮能量成像结果，其中图 2a-2c分别为采用自解耦能量互相关成像条件得到的PP波能量成像结果、C波能量成像结果和SS波能量成像结果，图 2d为弹性能量成像结果.从成像结果来看，总能量成像结果为PP波、C波、SS波能量成像之和，解耦之后可以发现PP波能量成像结果同相轴连续且振幅较强，垂直入射时能量最强，但是远偏移处同相轴连续性不如C波成像及SS波成像；相较而言，C波成像在近偏移距处振幅较弱，但远偏移距连续性较好，而且没有极性反转问题；类似于C波成像，SS波成像在近偏移距处振幅较弱，而且高波数成分较强.虽然成像结果中存在串波等噪声的干扰，但是能量弱，对反射轴成像结果影响可忽略不计.

为分析背向散射噪声压制效果，对图 1中模型采用相同震源及观测系统进行激发接收，分别采用自解耦成像条件及常规互相关成像条件进行成像，得到成像结果如图 3所示.其中，(a)和(b)为基于自解耦能量成像条件得到的PP波和SS波成像结果，(c)和(d)为基于自解耦成像条件得到的PP波和SS波成像结果.对比发现，常规互相关成像结果存在严重的背向散射噪声，对分界面处成像振幅造成影响，需要进行滤波；然而，自解耦成像结果不需要引入滤波，就可明显压制背向散射噪声，对分界面处成像振幅影响可忽略不计.因此，自解耦能量成像条件的背向散射压制效果理想.

2.2 Marmous2模型

对Marmious2原始模型(图 4)进行正演模拟获得地震记录，对平滑后的Marmous2模型(图 5)作为背景速度场进行逆时偏移.Marmious2原始模型中既有平层发育又有陡倾角断层发育，平滑后平层及断层处反射界面不再明显.实际模拟时采用1 ms时间采样间隔，进行时间长度为6 s的采样，中心放炮双边接收，炮间距为20 m，共156炮.

图 6-8分别为采用自解耦能量互相关成像条件得到的PP波偏移结果、SS波偏移结果及C波偏移结果.从成像结果可以看出，在复杂模型中可以得到纯波波场模式的能量成像结果，而且低波数噪声较为有效地得到压制.其中，PP波成像结果中地质构造能够清晰成像，包括中浅部的大倾角构造、断裂及深部的角度不整合地层、突变体均得到了比较好的成像.与PP波成像结果相比，SS波成像在不整合地层和突变体处能量更突出，但也存在突变体处噪声扰动及深层成像相对较弱的问题，这与横波传播速度慢接收时间有限有关，可以通过延长接收时间及衰减补偿等方法对深部能量进行补偿.相对比而言，转换波成像中背景噪声较大，这与模型复杂绕射散射现象突出等问题有关.

为分析基于自解耦能量成像条件的成像效果，对图 4中模型采用相同震源及观测系统进行激发接收，采用常规互相关成像条件进行成像，得到滤波前PP波成像结果如图 9所示.发现常规成像结果中，浅层背向散射噪声干扰严重，导致深层反射信息被掩盖.因此，引入Laplace滤波来压制浅层背向散射噪声，得到滤波后PP波成像结果如图 10所示.与自解耦成像条件得到的PP波成像结果(图 6)对比发现，滤波后常规互相关成像结果可压制浅层背向散射噪声.而自解耦成像条件得到的成像结果不需要额外引入滤波，即可压制背向散射噪声.因此，自解耦能量成像条件得到的成像结果在低波数噪声压制方面效果理想.

3 结论

基于能量密度构建的能量互相关成像条件可以得到能量成像结果，但是其将所有弹性波能量耦合在一起，难以区分纯波波场能量.为此，针对能量成像条件难以分离纯波能量的问题，本文引入Helmholtz原理和体应力构建出自解耦能量成像条件.该成像条件既可以得到纯波模式能量成像结果，而且避免波场分离过程对波场振幅和相位特征进行改造，从而避免了极性反转问题及额外计算量问题.数值实验验证该方法可得到低波数噪声有效压制的PP波、SS波和转换波成像结果，并刻画出地质构造中纯波波场的能量分配.以此为基础，可将自解耦能量互相关成像条件应用于偏移成像等方面，以便更好进行地震数据处理.

附录A 能量密度函数

假定弹性体在受力作用的过程中始终保持平衡，即弹性体动能恒定，而且弹性体非机械能恒定，则外力势能(也就是外力所做的功)完全转变为弹性体的形变势能(称之为应变能)存储于弹性体内部.根据能量守恒定律，形变势能与弹性体受力次序无关，完全取决于应力及形变的最终大小.因此，设弹性体受6个应力分量σ_x, σ_y, σ_z, τ_yz, τ_xz, τ_xy和6个形变分量ε_x, ε_y, ε_z, ε_yz, ε_xz, ε_xy，得到应变能密度定义如下(徐芝纶，1990)：

(A1)

弹性体的动能密度表示为：

(A2)

在有限空间域Ω的能量分布可以用能量密度函数表示：

(A3)

则有限空间中某时刻的总能量即：

(A4)

对于各向同性介质c_ijkl=λδ_ijδ_kl+μδ_ikδ_jl+μδ_ilδ_jk而言，利用σ_ij=c_ijklε_kl，则应力势能密度可以表示为：

(A5)

基于本构方程，则基于位移的应力势能密度可以表示为：

(A6)

对有限空间中总能量求时间偏导：

(A7)

根据链式法则，等式右边后三项分别表示为：

(A8)

(A9)

(A10)

将式(A8)(A9)(A10)代入式(A7)有：

(A11)

根据散度定理，边界处衰减为0，则式(A7)可以变换为：

(A12)

当没有外力作用，即f=0时，，根据能量守恒定律，在没有外力作用下，弹性体中总能量不发生变化，即具有唯一性和稳定性.