0 引言

逆时偏移(Reverse Time Migration，RTM)具有无倾角限制、适应任意复杂模型的优点，是目前复杂地质构造偏移成像的重要方法(Baysal et al., 1983；McMechan, 1983；Whitmore, 1983).相较声波逆时偏移，弹性波逆时偏移(ERTM，Elastic Reverse Time Migration)假设地球介质为弹性介质，允许纵波和横波同时传播，更符合地震波实际传播情况，因此弹性波逆时偏移能够获取更准确的成像结果.

弹性波逆时偏移主要包括两部分：构建震源波场和检波波场以及利用成像条件提取成像信息.早期，弹性波逆时偏移(Chang and McMechan, 1986; Chang, 1987)基于笛卡尔分量构建弹性波位移分量场，以提取到偏移成像结果，该结果中存在多波模式串扰等问题.为避免串扰问题，后续研究多围绕输出物理意义明确的P波和S波成像结果展开.根据是否保持弹性波场的矢量特征，有标量处理方法和矢量处理方法两种处理思路.其中，标量处理方法是基于声波波动方程构建标量波场(Sun and McMechan, 2001; Sun et al., 2006)，虽然充分利用声波逆时偏移已有成熟的处理技术，但是忽略了弹性波波形转换等矢量特征.矢量处理方法则基于弹性波动方程进行波场延拓，在延拓过程中可以有效保持弹性波场的矢量特征，从而构建得到符合地震波实际传播过程的震源波场和检波波场.

构建准确的成像条件则是提取成像结果实现弹性逆时偏移的另一关键环节.最早，激发时间成像条件(Chang and McMechan, 1986, 1994)根据传播时间提取检波波场的波场值作为成像结果，该成像条件计算和存储消耗较低，但振幅保持性较差.为此，提取检波波场和震源波场的比值作为成像结果(Nguyen and McMechan, 2012; 张智等，2013)，该成像条件兼顾了计算效率与振幅保持特性，但是弹性转换波成像的合理性需进一步研究.目前常用的互相关成像条件(Claerbout, 1971; Loewenthal et al., 1991; Du et al., 2017)，通过对相同时刻的源检波场进行互相关运算提取成像结果.该成像条件稳定性高且充分利用全波场信息，但是受双程波动方程的影响，会在波传播路径上产生背向散射.

针对背向散射噪声干扰问题，目前压制技术主要分为两大类(杜启振等，2013)：一是进行波场构建时避免背向散射(Baysal et al., 1983; Loewenthal et al., 1987)；二是选择性成像，利用根据波传播方向所确定的散射角进行小角度成像(Liu et al., 2007)，或者成像后对成像结果进行滤波处理(Sava, 2007).该压制技术虽然可以成功衰减背向散射，但是同样会损害振幅相位等有效信息.然而，目前发展的能量互相关成像条件，基于能量密度将空间域互相关成像条件拓展为时间-空间域成像条件(Rocha et al., 2016a, b, 2017;张晓语等，2018)，提取得到弹性波场反射能量分布信息的同时，可以有效压制背向散射噪声.

但是，目前构建的能量互相关成像条件由位移空间导数及时间导数组成，应用于一阶弹性波方程时需要将振动速度及应力转换为位移变量，难以直接应用.为此，我们针对振动速度及应力，基于满足能量守恒的弹性能量密度表达式，构建速度-应力表征的能量范数.将能量范数拓展为速度-应力表征的能量内积，在震源端与检波端的基矢量正方向保持一致的基础上，构建得到背向散射有效压制的能量互相关成像条件.

1 原理 1.1 一阶弹性波方程

各向同性介质中，一阶弹性波方程为(Levander，1988)：

(1)

式(1)中，υ=(v_x, v_y, v_z)^T，σ=(σ_xx, σ_yy, σ_zz, σ_yz, σ_xz, σ_xy)^T，，.其中，υ是质点振动速度矢量，σ是应力矢量，ρ是密度，t是时间, L是偏导数矩阵，C是刚度矩阵.在各向同性介质中，刚度矩阵C可以由2个独立的弹性模量即拉梅系数λ和μ表述.考虑到弹性波场传播速度有限，认为空间区域Ω足够大从而确保边界∂Ω处满足υ、σ及其导数为0.

一般情况下，二阶位移弹性波动方程假定介质速度及密度的空间导数为零，忽略了介质速度及密度空间变化对弹性波场的影响，尤其是密度变化对弹性波场的影响.相较于二阶弹性波方程，一阶弹性波方程考虑到介质速度及密度空间变化的影响.因此，一阶弹性波方程更适应于非均匀介质情况，尤其是密度非均匀的情况.

1.2 速度-应力表征的弹性波场能量

基于一阶弹性波方程构建的弹性波振动速度-应力场，可以表征弹性波场能量.在动力学状态下，弹性体处于运动平衡状态，t时刻x处动能密度可表示为

(2)

其中，T为t时刻x处的动能密度，ρ为密度.在静力学状态下，不考虑外力影响，t时刻x处势能密度可表示为

(3)

其中，W、V_e和V_σ为t时刻x处势能密度、应变能密度和应变余能密度，C^-1为柔度矩阵.线性弹性介质中，弹性体势能密度可以由应变能密度或应变余能密度来表征，两者在数值上相等，前者为应变的状态函数，后者为应力的状态函数，因此基于应力采用应变余能密度来表征势能密度(徐芝纶，2006).

那么，动能密度与势能密度之和为弹性波场能量密度(杜世通，1996)，公式为

(4)

区域Ω中弹性波场能量可由能量泛函表征公式为

(5)

式中，E(t)表示t时刻区域Ω中的总弹性波场能量.不考虑外力影响时，区域Ω中弹性波场为保守力场，其总弹性波场能量满足(胡德绥, 1989；Rocha et al., 2016a, b)，即弹性能量在弹性介质区域Ω中传播时满足机械能守恒定律.对于弹性各向同性介质，式(5)可表示为

(6)

为区别能量密度和振动速度，分别采用e和g来表示杨氏模量和泊松比，与拉梅系数λ和μ满足关系和，Θ(x, t)=[σ_x(x, t)+σ_z(x, t)]为体积力.这样，基于质点振动速度和应力可以实现弹性波场能量的表征.其中，第一项为动能密度，由质点振动速度组成，即质点位移的时间导数组成，表示弹性波动力学特征；第二项为势能密度，由应力组成，即质点位移的空间导数组成，表示弹性波静力学特征.

1.3 能量范数及能量互相关成像条件

L2范数可表征弹性波场强度，为此常规互相关成像条件通过求取震源波场与检波波场间内积，以提取弹性波场的反射强度.类比位移表征的能量范数(Tanushev et al., 2009; Rocha et al., 2016a, b, b)，针对基于弹性波场振动速度-应力场U=(υ, σ)，结合能量泛函E(t)，将能量范数定义为

(7)

其中，‖·‖_E²为能量范数，表示弹性波场的能量强度，由弹性波场振动速度及应力表征，i表示虚数.

那么，t₀时刻弹性波场U_S的能量范数，可以提取空间Ω中瞬时弹性波能量分布特征.相应的，弹性波场U_R和U_S间的能量内积定义为

(8)

其中，〈·, ·〉_E表示波场间的能量内积，表征波场间能量强度，由弹性波场的质点振动速度及应力矢量构成.

那么，t₀时刻弹性波场U_R与U_S间能量内积，可以提取空间Ω中弹性波场间瞬时能量分布特征.将能量内积应用于逆时偏移技术中，即对震源波场与检波波场求取能量内积，可以提取t₀时刻空间Ω中反射波场的能量分布特征，沿延拓时间[0, T]进行积分可得到延拓过程中反射能量分布特征，那么能量互相关成像条件可构建为

(9)

其中I_E表示采用能量互相关成像条件得到的成像结果，表示弹性波场的反射能量，υ_{R, i}、υ_{S, i}、σ_{R, j}和σ_{S, j}分别表示震源波场与检波波场的振动速度和应力.震源波场U_S(x, t)和检波波场U_R(x, t)同时传播到反射点x₁处时，〈U_R, U_S〉_E≠0，反之〈U_R, U_S〉_E=0.其中，能量内积非零处为弹性波反射处，非零能量内积即为阻抗分界面x处反射能量.因此，基于振动速度及应力构建的能量成像条件(式9)，可以提取弹性波场在波阻抗分界面处的反射波场能量.

1.4 背向散射压制

基于双程波动方程的逆时偏移采用互相关成像条件得到的成像结果存在背向散射干扰(杜启振等，2013)，然而采用能量互相关成像条件可有效压制背向散射(Rocha et al., 2016a, b)，同样地，采用基于速度-应力表征的能量互相关成像条件可有效压制背向散射.

各向同性弹性介质中，基于一阶弹性波方程构建的振动速度及应力场满足哈密尔顿定理(杜世通，1996)，公式为

(10)

在能量不泄露的保守力场中，动能与势能之差的时间积分为0.求取震源波场和检波波场能量内积时，为保证震源端与检波端的基矢量正方向保持一致，检波波场进行逆时延拓时，时间传播方向应与正向延拓时的时间传播方向的正方向保持一致.由于逆向延拓的时间传播方向与正向延拓时间的传播方向相反，所以用-υ_{S, i}表示逆时延拓的质点振动速度.那么，震源波场U_R=(υ_R, σ_R)和检波波场U_S'=(-υ_S, σ_S)的能量内积有：

(11)

对于传播路径相同及极化方向一致的低波数噪声，σ_{R, j}=σ_{S, j}, υ_{R, i}=υ_{S, i}，利用式(10)有〈U_R, U_S'〉_E=0，即能量内积为零.反射点处，震源波场与检波波场的传播路径及极化方向不再保持一致，则〈U_R, U_S'〉_E≠0.那么，将能量互相关成像条件(式9)改进为

(12)

其中I_E'为采用能量互相关成像条件得到的成像结果，该能量互相关成像条件(式12)可得到表征反射能量的成像结果，同时对于满足传播路径相同及极化方向一致的背向散射即σ_{R, j}=σ_{S, j}, υ_{R, i}=υ_{S, i}，其能量内积为零，那么可有效压制传播路径及极化方向相同的背向散射.

2 数值模拟

用简单平层模型和二维Marmousi2模型，对成像条件的稳定性、背向散射压制问题、成像及振幅恢复能力等进行测试.

2.1 简单模型

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

首先分析弹性波场能量传播特点，提取600 ms时正向延拓的能量快照(图 2a)和速度快照(图 2b)，以及600 ms时反向延拓的能量快照(图 2c)和速度快照(图 2d).对比发现，弹性波场能量传播伴随弹性波场传播，其传播速度保持一致，但是相位特征不同，其中能量快照相位极性始终保持一致，但位移快照中相位极性有正有负.然后分析简单模型的单炮偏移成像(图 3)，发现基于能量互相关成像条件得到的偏移结果中，同相轴连续性好，虽然存在串扰噪声，但是噪声能量弱，对反射轴成像结果影响可忽略不计.

不同于能量互相关成像条件，常规互相关成像条件对弹性波速度场进行互相关从而提取成像结果.基于相同震源及观测系统进行激发接收，采用常规互相关成像条件对简单模型进行偏移成像，得到成像结果如图 4所示.其中，(a)为滤波前单炮成像结果, (b)为进行Laplace滤波后的成像结果.对比发现，互相关成像结果中存在背向散射噪声，Laplace滤波后可有效压制背向散射噪声，但是反射同相轴振幅难以保持，而且串扰噪声干扰相对突出.然而，能量互相关成像结果中背向散射噪声压制效果明显，不需要滤波，而且串扰噪声相对较弱.提取x=1500 m处单道能量成像结果及滤波后互相关成像结果(图 5)，发现远偏移距处成像能力相当，但是近偏移距处能量成像结果中振幅强度强于后者.综合来说，在成像能力及背向散射压制方面，能量互相关成像条件优势明显.

2.2 Marmousi2模型

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

图 8为基于能量互相关成像条件得到的弹性波能量成像结果，图 9为基于常规互相关成像条件得到的弹性波成像结果，图 10为Laplace滤波后互相关成像结果.在能量互相关成像结果(图 8)中，中浅部的大倾角构造及断裂构造清晰；气藏区成像结果虽然存在些许绕射能量，干扰了气藏下储层成像结果，但仍可以准确追踪到其上下储层界面的构造特征；深部的角度不整合地层及高陡地质体也实现了比较好的构造成像.在背向散射压制方面，对比发现互相关成像结果(图 9)中背向散射突出，浅层强反射界面处背向散射尤为严重，会影响到深层成像效果，滤波后的互相关成像结果(图 10)中背向散射噪声可有效压制.相对于互相关成像结果，能量互相关成像结果(图 8)不需要滤波，就可有效地压制背向散射噪声的干扰，而且串扰噪声的干扰较弱、同相轴连续性强.因此，能量成像条件的成像能力及背向散射衰减能力优势明显.

3 结论

针对能量互相关成像条件难以直接应用于一阶弹性波方程的问题，本文基于一阶弹性波方程及能量密度，构建以振动速度和应力为参数的能量范数，实现弹性波场能量的表征.将能量范数拓展为震源波场与检波波场的能量内积，以提取弹性波场反射能量.结合逆时偏移中震源端与检波端的基矢量方向，构建得到可有效压制背向散射的弹性波能量成像条件.数值模拟结果表明：该成像条件可以得到物理意义明确、背向散射压制、振幅有效保持的能量成像结果.