0 引言

我们的地球广泛存在黏弹性和各向异性.地下介质的黏弹性导致高频能量衰减及波形相位变化(李振春和王清振，2007；Qu et al., 2017；李金丽等，2015).因此，过去几十年间，对衰减补偿方法研究的兴趣与日俱增.一种简单直接的方法是使用反Q滤波补偿因介质黏弹性造成的衰减(Bickel and Natarajan, 1985；Gelius，1987；Hargreaves and Calvert, 1991；Varela et al., 1993；裴江云和何樵登，1994；高军等，1996；姚振兴等，2003).然而，由于反Q滤波是基于一维波场反传的，所以无法准确处理实际复杂地质构造(Zhang et al., 2010).为了提高复杂地质构造的成像分辨率，反Q偏移法逐渐得到了发展(崔建军和何继善，2001；杨午阳等，2007；李雪英等，2011).基于射线追踪的反Q偏移法沿射线路径计算衰减效应，在偏移过程中对每个频带进行补偿(Traynin and Reilly, 2008；Xie et al., 2009).基于单程波方程的反Q偏移方法通常在频域中实现，可较好地校正黏弹性的影响(Dai and West, 1994；Mittet and Sollie, 1995；何兵红等，2009).近年来，一些Q补偿的逆时偏移(RTM)方法已被广泛地提出来(Deng and McMechan, 2007；Zhang et al., 2010；Fletcher et al., 2012).

然而，基于传统黏声波方程的逆时偏移，因振幅衰减和相位频散耦合在一起，导致相位频散无法校正(Xue et al., 2018)，因此，许多振幅衰减和相位频散影响解耦的方程被提出，并应用于衰减补偿的逆时偏移中(Zhu et al., 2014；Guo et al., 2016；Qu et al., 2018).对于Q补偿的RTM方法，在衰减补偿过程中，高频成分的指数增长会产生严重的数值噪声，引起不稳定.为了解决这个问题，正则化算子方法(Zhang et al., 2010；Qu et al., 2015a；Wang et al., 2017)、高通滤波方法(Zhu et al., 2014)及黏介质伴随算子方法(李振春等，2014；Guo and McMechan, 2018；李金丽等，2018)可以用以保证衰减补偿过程中的稳定性.

当地下介质存在各向异性时，如果忽略各向异性的影响，会导致成像位置与真实位置存在较大偏差，需要研发能够校正各向异性的成像技术，使得成像位置归位.弹性各向异性方程能够准确模拟地震波在各向异性介质的传播规律(熊煜等，2006；张美根和王妙月，2001；杜启振和秦童，2009；朱峰等，2018)，但是该方法应用于成像时，纵横波无法准确分离且计算量较大，因此，拟声波方程被广泛地应用于各向异性成像中，针对于此，地球物理学者提出了各种拟声波方程.Alkhalifah(1998)将对称轴方向的横波速度赋为零，但是该方法存在不稳定现象，且横波干扰没有压制彻底(Grechka et al., 2004；张岩和吴国忱，2013；李博等，2012).因此多种改进的各向异性拟声波方程被提出来(Duveneck et al., 2008；Yoon et al., 2010；程玖兵等，2013；梁锴等，2009).

基于反演成像理论的最小二乘逆时偏移逐渐得到了国内外地球物理学者的广泛关注(Dai et al., 2012；黄建平等，2014；Gang and Wu, 2015；Zhang et al., 2015；张攀和毛伟建，2018；方修政等，2018).相对于常规偏移算法，最小二乘逆时偏移具有信噪比高、保真度高、分辨率高等优点.对于最小二乘逆时偏移方法早期的研究是基于弹性假设的，但实际上，地下介质不是完全弹性的，存在明显的黏弹性和各向异性.于是，国内外地球物理学家针对复杂介质条件(黏弹性和各向异性等)下的最小二乘逆时偏移进行了研究.为了校正地下介质黏弹性的衰减影响，Dutta and Schuster(2014)及李振春等(2014)通过构建基于广义标准线性固体(GSLS)的不含记忆变量的黏声波动方程，推导并实现了黏声最小二乘逆时偏移方法与技术，在成像过程中校正了介质黏弹性对地震波振幅衰减和相位的影响.杨宇等(2015)发展了基于Kelvin模型的黏声最小二乘逆时偏移方法对品质因子Q进行反演.针对各向异性对地震波运动学和动力学的影响，一些各向异性拟声波最小二乘逆时偏移被发展起来(如李振春等，2017).

本论文首先对最小二乘逆时偏移的基本理论进行描述，然后分别对最小二乘逆时偏移中的黏弹性和各向异性进行校正，最后对黏弹性和各向异性同时校正.在实际资料试算部分，针对渤海湾数据进行测试，验证对最小二乘逆时偏移中的黏弹性和各向异性校正的效果.

1 方法原理 1.1 最小二乘逆时偏移基本理论

最小二乘逆时偏移(Least-Squares Reverse Time Migration, LSRTM)通过最优化算法寻找一个最优的成像结果使得模拟数据与观测数据最佳匹配.定义模拟数据与观测数据残差的2范数为误差函数:

(1)

其中，d_obs表示输入的观测数据，d_cal表示模拟的地震记录，L表示线性正演算子，m表示反射系数模型.在本文中，我们使用预条件的共轭梯度方法(Dai et al., 2012)来求解最小二乘的反演问题.共轭梯度方法使用下述方程更新反射系数模型：

(2)

其中，m^(k)和m^(k+1)表示第k与第k+1次迭代的反射系数模型，α^(k)和z^(k)分别表示共轭梯度法第k次迭代的步长和梯度方向.它们由以下四个方程求得：

(3)

(4)

(5)

(6)

其中，g^(k+1)和β^(k+1)分别为由最速下降法求解获得的梯度方向与步长.

在各向同性声波介质中，地震波动方程为

(7)

其中，v为声波速度，t为时间，x为空间坐标，▽²为拉普拉斯算子，f为震源，u为声波波场值.传统声波最小二乘逆时偏移的反偏移算子为(附录A)：

(8)

(9)

(10)

其中，u₀为背景波场，δu为残差波场，v₀为背景速度，t为时间，x_s为震源坐标，x_r为检波点坐标.根据伴随状态理论(Plessix，2006)，偏移算子可由下式求得：

(11)

其中，u^R为伴随波场，Δd表示数据残差，由下式求得：

(12)

传统最小二乘逆时偏移的梯度公式为

(13)

最小二乘逆时偏移的具体实现流程如下：

(1) 当第一次迭代(k=0)时，d_cal=0；

(2) 采用方程(8)和(11)计算u₀和u^R；

(3) 采用互相关成像条件计算成像剖面：

(14)

(4) 当迭代次数大于1(k>0)时，使用反偏移算子方程(8)—(10)计算合成地震数据d_cal；

(5) 利用方程(12)计算数据残差Δd，判断是否满足收敛条件，如果满足，则进入步骤(9)，否则进入步骤(6)；

(6) 采用方程(8)和(11)计算u₀和u^R；

(7) 利用方程(13)计算梯度，并利用方程(4)—(5)计算共轭梯度方向和步长；

(8) 利用方程(2)更新反射系数，k=k+1，返回步骤(4)；

(9) 输出成像剖面.

1.2 最小二乘逆时偏移中的黏弹性校正

当地下介质存在明显的黏弹性时，需要对黏弹性影响进行校正.基于GSLS模型的一阶波动方程在波场反传时存在一些困难.而且该方程存在记忆变量，因此规则化算子很难应用.为了克服这一缺点，很多不含记忆变量的黏声拟微分波动方程被推导出来，用来实现黏声介质的逆时偏移成像.本论文采用Bai等(2013)提出的黏声拟微分方程：

(15)

其中，▽²为拉普拉斯算子，τ=τ_ε/τ_σ－1由Q值确定(Liu et al., 1976).应力松弛时间τ_σ和应变松弛时间τ_ε可通过τ-Q关系计算得到(Carcione, 2001)：

(16)

(17)

其中，ω为角频率.式(15)中第二项为衰减项.在该黏声拟微分方程中不含记忆变量项.当τ等于零时，黏声拟微分波动方程变为声波方程.使用高阶有限差分法求解方程(15)，其中衰减项在波数域计算，其他项可在波数域或者时空域求解(附录C).

黏声最小二乘逆时偏移的反偏移算子应由式(8)和(9)变换为(附录A)

(18)

(19)

黏声介质残差波场的反向传播方程为

(20)

在波场反传过程中，衰减项的符号应由正号改为负号以补偿衰减.当地震波在衰减介质中进行衰减补偿时，会产生严重的高频分量数值噪音，导致振幅的巨大振荡.在本论文中，我们引入规则化项 (Qu et al., 2015a)来稳定衰减补偿的波形.因此，黏声最小二乘逆时偏移的偏移算子为

(21)

其中，σ规则化参数，本论文中σ=0.02.

1.3 最小二乘逆时偏移中的各向异性校正

对于各向异性介质，Duveneck等(2008)推导的二维常密度VTI拟微分方程被广泛地应用于各向异性介质逆时偏移中，其表达式为

(22)

其中，p和q分别表示波场；v_pz为对称轴方向的相速度；ε和δ为各向异性参数.VTI介质的反偏移方程为(附录B)

(23)

其中，p₀和q₀为背景波场，v_pz⁰为背景相速度，背景波场可由方程(22)计算得到.波场p和q的伴随变量p^R和q^R满足如下的偏移算子方程：

(24)

其中，Δd_p和Δd_q可由下式求得：

(25)

误差泛函f需要同时考虑p和q分量，那么，VTI介质最小二乘逆时偏移的梯度为

(26)

1.4 最小二乘逆时偏移中的衰减和各向异性同时校正

地下介质通常同时存在黏弹性和各向异性，在应用最小二乘逆时偏移时，如果没有对黏弹性和各向异性进行校正，模拟数据与观测数据无法准确地匹配，因此误差函数无法正确收敛，导致成像结果较差.因此，当黏弹性和各向异性影响较强时，需要同时校正地下介质的黏弹性和各向异性.

采用如Suh等(2012)和徐文才等(2016)同样的方法，将方程(22)加入公式(15)中的衰减项，得到黏声VTI介质拟微分方程(Qu et al., 2015b)：

(27)

同样采用混合域差分格式求解方程(27).

黏声VTI介质最小二乘逆时偏移的反偏移算子为(附录B)

(28)

方程(28)的右端项为震源项，其中包括背景波场和反射系数模型.在反偏移方程的右端项，一阶偏导数控制着能量的衰减.当方程(28)变成各向同性波动方程的形式并保留着方程两端的衰减项时，该方程可以简化为黏声各向同性反偏移方程(19)，当去掉方程(28)两端的衰减项时，该方程可以简化为声波各向异性反偏移方程(23)，当移除方程(28)中的衰减项并将其变成各向同性波动方程的形式，该方程可简化为常规的声波各向同性反偏移方程(9).

将数据残差作为方程(27)的震源项，并将衰减项变为负号，同时引入两个规则化项和来稳定衰减补偿的波形，可得到黏声VTI介质最小二乘逆时偏移的偏移方程：

(29)

2 实际资料试算 2.1 偏移参数场和观测系统

对渤海湾采集的地震数据进行偏移成像，测试对最小二乘逆时偏移中黏弹性和各向异性的校正效果.图 1所示的偏移参数场，偏移模型大小为16250 m×6000 m，其中图 1a为偏移速度场，图 1b为品质因子场，通过经验公式由速度场计算得到，图 1c和图 1d为各向异性参数场，偏移速度场和各向异性参数场都是通过多参数层析方法得到.从图 1a所示的偏移速度中可以看出，最小速度为海水的速度，约为1500 m·s^-1，最大速度为4600 m·s^-1.各向异性参数中ε和δ的最大值分别为0.1和0.06，表现出较为明显的各向异性影响.

采集地震数据的观测系统为：采用海洋拖揽进行采集，即在海水面处激发，在海水面处接收.与海底电缆方式不同的是，该采集方式无法采集到横波信息.采集到的炮集总数为399炮，每炮共1300道.时间采样率为2 ms，时间采样点数为3072，最大偏移距为3564 m，最大深度为10000 m，相邻检波器之间的距离为12.5 m，相邻炮点距离为37 m.采用左侧激发方式，记录得到的观测数据如图 2所示.在最小二乘逆时偏移的模拟部分，采用纵横向网格间距为12.5 m×16 m，网格点数为1301×401.时间采样间隔为1.0 s，采样点数为4144，采用CFS-PML(Kuzuoglu and Mittra, 1996)边界条件压制人工边界反射，差分精度为10阶，使用25 Hz主频的雷克子波进行模拟.在2.2节中，将对该渤海湾地震数据进行成像.

2.2 最小二乘逆时偏移结果

采用同时校正衰减和各向异性的最小二乘逆时偏移得到的成像剖面如图 3所示，其中图 3a为第一次迭代的成像结果(相当于逆时偏移的成像结果)，图 3b为第30次迭代的成像结果，这里的结果均未进行低频噪音的滤波(如Laplacian滤波等).从两图对比可以看出，采用最小二乘逆时偏移迭代30次之后的成像结果(图 3a)与逆时偏移的成像结果(图 3b)相比，浅层的低频噪音得到了明显地压制，如图虚线椭圆所示.中深层能量得到了补偿，如图实心箭头所示，潜山顶部构造更加连续且清晰，如图虚线箭头所示，潜山内部构造得到了凸显，如图空心箭头所示，除此之外，在中浅层区域的断层构造也得到了很好地刻画.在本例中，潜山内幕构造为需要精细成像的地质目标体.众所周知，潜山内幕构造成像困难，如图 3a中空心箭头所示，潜山内幕构造几乎没有成像能量，而经过本论文最小二乘逆时偏移30次迭代过程中对衰减和各向异性的校正，潜山内幕构造得到了准确地成像.也就是说，通过最小二乘逆时偏移的迭代，得到了更高质量的成像结果.图 4为误差泛函的收敛曲线(即数据误差收敛曲线).从图中可以看出，随着迭代次数的增加，误差泛函首先急剧下降，随后下降速度逐渐变缓，最终稳定在0.2以下.图 5为合成地震记录与采集的地震数据，其中图 5a为采集的地震数据，图 5b为初始模型正演得到的炮记录，图 5c为第一次迭代反偏移炮记录，图 5d为第30次迭代反偏移炮记录.从图中可以看出，初始模型正演得到的炮记录只有几个弱的反射波形，第一次迭代反偏移炮记录出现了较多反射同相轴，走势与采集的地震数据基本吻合，但浅层只包含少量的波形信息(实心箭头所示)，而且深层的波形能量较弱(空心箭头所示)，中部区域的反射轴分辨率较低(虚心椭圆所示).从第30次迭代的反偏移炮记录中可以看出，浅层区域出现了更多的反射同相轴，与输入的地震记录更加匹配，而且中部区域的反射轴分辨率得到了一定程度的改善，深部反射能量大幅度提升，概括来讲，30次迭代后的反偏移炮记录更加匹配输入的地震记录.以上结果证明了同时校正衰减和各向异性的最小二乘逆时偏移的正确性，该方法在应用于实际数据时也有非常好的效果.

接下来，我们对采集到的地震数据分别采用常规最小二乘逆时偏移、只校正黏弹性的最小二乘逆时偏移、只校正各向异性的最小二乘逆时偏移进行试算，以验证最小二乘逆时偏移在实际生产应用中对黏弹性和各向异性校正的优越性和必要性.这里，偏移使用的观测系统与前文相同，当采用常规最小二乘逆时偏移时，只输入偏移速度场(图 1a)，当采用只校正各向异性的最小二乘逆时偏移时，输入的是偏移速度场(图 1a)和品质因子场(图 1b)，当采用只校正黏弹性的最小二乘逆时偏移时，输入的是偏移速度场(图 1a)和各向异性场(图 1c和图 1d).

图 6所示为采用四种最小二乘逆时偏移第一次迭代的成像结果，其中图 6a为普通声波LSRTM的成像结果，图 6b为只校正衰减性的LSRTM的成像结果，图 6c为只校正各向异性的LSRTM的成像结果，图 6d同时校正衰减性和各向异性的LSRTM的成像结果.相比于同时校正衰减性和各向异性的LSRTM的成像结果(图 6d)，只校正黏弹性的LSRTM的成像结果(图 6b)在中浅层区域(虚线椭圆区域)同相轴发生错乱，连续性变差，成像位置没有准确归位，主要断层没有成像出来；只校正各向异性的LSRTM的成像结果(图 6c)中深层能量较弱(实心箭头所示)，而且深部区域的能量几乎没有(如图空心箭头所示)，因此，衰减的能量没有得到补偿.普通声波LSRTM的成像结果中的黏弹性和各向异性能量都没有得到校正，因此同相轴出现明显错乱，而且深部能量较弱.而同时校正衰减性和各向异性的LSRTM的成像质量相比于图 6a、6b和6c得到了明显的改善.同时校正衰减性和各向异性的LSRTM成像对于目的层(潜山顶部和潜山内幕)构造得到了更加精准地成像.

图 7为四种最小二乘逆时偏移第30次迭代的成像结果，其中图 7a为普通声波LSRTM的成像结果，图 7b为只校正衰减性的LSRTM的成像结果，图 7c为只校正各向异性的LSRTM的成像结果，图 7d同时校正衰减性和各向异性的LSRTM的成像结果.普通声波LSRTM和只校正黏弹性的LSRTM的成像结果(图 7a和7b)，因为未校正各向异性的影响，除了中浅层区域同相轴错乱以外，随着迭代次数的增加，同相轴位置错乱的情况没有得到改善，且分辨率也没有得到提升，反而分辨率变差了(由图 7a和7b中的实线箭头所示).普通声波LSRTM和只校正各向异性的LSRTM的成像结果(图 7a和7c)，因为没有补偿衰减的影响，随着迭代次数的增加，潜山顶部和潜山内幕构造的成像能量没有得到明显地改善，同相轴依旧较弱，而且不连续(如图 7a和7c中空心箭头所示).图 8为四种最小二乘逆时偏移的误差泛函收敛曲线，其中红色、绿色、蓝色和粉色曲线分别表示传统声波LSRTM、只校正各向异性的LSRTM、只校正黏弹性影响的LSRTM和同时校正黏弹性和各向异性的LSRTM的误差泛函收敛曲线.从图中可以看出，采用本论文同时校正黏弹性和各向异性的LSRTM具有更快的收敛速度，而且稳定的误差泛函最小值也更低.传统声波LSRTM的收敛速度最慢，而且归一化误差泛函最小值在0.3左右，为四种方法的最大值.其他只单独校正一种性质的LSRTM收敛速度介于两者之间，但是该例子不具有一般性，因为不同探区各向异性和黏弹性的影响不同.

综上所述：同时校正黏弹性和各向异性的LSRTM在低频噪音压制、震源效应压制、深部能量改善、分辨率提升、目的层刻画等方面都明显优于其他三种LSRTM的成像结果.

3 讨论与结论

地下介质存在黏弹性和各向异性，而非完全弹性.最小二乘逆时偏移受黏弹性和各向异性的影响比其他成像方法(如逆时偏移)更严重，因为最小二乘逆时偏移本身就是通过不断迭代使得观测数据与模拟数据的误差不断降低的过程，如果不校正黏弹性和各向异性的影响，在匹配观测数据和模拟数据时，无法准确匹配，会导致成像误差的逐渐累加.

本论文以渤海地区数据为例，校正最小二乘逆时偏移过程中的黏弹性和各向异性，得到了较好的成像结果.本论文是针对二维理论和数据进行测试的，方法可以很容易地发展到三维.本论文方法除了具备传统最小二乘逆时偏移压制低频噪音、提高信噪比、增加保幅性、压制采集脚印、提高振幅均衡性和补偿深部能量等优点，还具备以下几点优势：(1)相比于未补偿黏弹性的最小二乘逆时偏移，深部区域衰减的能量得到了补偿；(2)相比于未校正各向异性的最小二乘逆时偏移，成像位置得到了准确地归位，随着迭代次数的增加，分辨率逐渐得到改善；(3)相比与其他三种传统方法，收敛速度更快.

但本论文的方法也存在一些需要解决的问题：(1)本论文使用的品质因子Q是通过经验公式求得的，采用更准确的Q能够得到更准确的成像结果；(2)该方法没有考虑密度的影响；(3)假设各向异性参数和品质因子不存在扰动，还不是真正的黏声各向异性LSRTM；(4)最小二乘逆时偏移自身问题需要进一步解决，如对偏移速度的依赖性，计算量大的问题，抗噪性差等；(5)本论文对各向异性的补偿是针对VTI介质的，地球介质还表现为方位各向异性、极化各向异性及正交各向异性等特征，采用本论文方法无法对这些性质的影响进行校正.

附录A 声波和黏声最小二乘逆时偏移的反偏移公式

黏声介质拟微分方程为

(A1)

在声波介质中，τ=0，方程(A1)退化为声波方程，因此声波最小二乘的反偏移算子也可通过黏声拟微分方程推导出来，这里不做单独推导.

基于伯恩近似理论，参数扰动δv=v－v₀能够引起如下所示的波场扰动：

(A2)

式中，u₀表示背景波场，v₀表示背景速度，δu表示扰动波场，δv表示扰动速度.背景波场u₀可由下式确定：

(A3)

采用泰勒展开近似得到

(A4)

其中，O(δv)表示v的高阶项.u(x, t)=u₀(x, t)+δu(x, t)代入方程(A1)中，然后减去方程(A3)可得

(A5)

我们定义反射系数模型为

(A6)

将方程(A6)代入方程(A5)并忽略高阶项可得黏声最小二乘逆时偏移反偏移方程：

(A7)

在声波介质中，τ=0，方程(A7)简化为声波最小二乘逆时偏移反偏移方程：

(A8)

附录B 各向异性和黏声各向异性最小二乘逆时偏移的反偏移公式

黏声VTI介质拟微分方程为

(B1)

根据伴随状态理论，相速度v_pz包括两部分：背景速度v_pz0和速度扰动δv_pz.那么，真实的速度模型可由下式描述：

(B2)

为了简化，本论文假设参数ε, δ和τ是准确的，且不存在扰动，即dε=0，dδ=0，dτ=0.那么，波场p和q可以被分成两部分：背景波场(p₀，q₀)和扰动波动(dp，dq).背景波场和扰动波场的关系可以由下式表示：

(B4)

背景波场可由下式计算得到

(B5)

从方程(B1)减去方程(B5)，并使用(A4)所示的泰勒展开近似，可得

(B6)

反射系数模型定义为

(B7)

将反射系数(B7)代入方程(B6)并忽略高阶项可得黏声VTI介质最小二乘逆时偏移的反偏移算子：

(B8)

声波VTI介质最小二乘逆时偏移的反偏移算子，只需要将(B8)中的τ变为0即可：

(B9)

附录C 混合域差分格式

我们使用时间上二阶，空间上2M阶的混合域差分格式求解方程(15)和(27).以方程(15)为例，对于变量u，时间二阶偏导数为

(C1)

衰减项的空间偏导数在波数域中求解：

(C2)

其中，k_x和k_z为分别为水平分量和垂直分量的波数.F为傅里叶变换，F^－1为傅里叶反变换，Δt为时间采样间隔.其他项在时间域求解：

(C3)

(C4)

其中Δx和Δz分别为水平方向和垂直方向的空间采样间隔；c为差分系数.对于变量p，上标和下标分别为时间和空间坐标.则方程(15)的递推公式为

(C5)

致谢  感谢中国石油大学(华东)地震波传播与成像实验室的技术支持.感谢中海油天津分公司的数据支持.