声波最小二乘偏移不同精确伴随算子对的定量关系分析

2022-06-11 01:24王建森任玉晓杨传根许新骥

石油地球物理勘探 2022年3期

王建森任玉晓陈磊严冬杨传根许新骥

(①山东大学齐鲁交通学院，山东济南 250061； ②山东大学岩土与结构工程研究中心，山东济南 250061； ③华能西藏水电安全工程技术研究中心，西藏林芝 860000)

0 引言

在反射波地震勘探中，偏移是利用地震记录对地下结构进行成像的重要手段。近几十年来，地球内部及近地表构造研究的需求推动了偏移技术的迅速发展。然而，当面对地震观测数据不足或不规则、地下结构复杂和波场带宽有限时，传统偏移技术如Kirchhoff偏移和逆时偏移(RTM)很难将数据偏移到准确的位置，难以满足实际地震勘探需求。

随着计算能力的飞速发展，提出了基于最小二乘优化框架的地震偏移方法，在真实速度扰动模型成像方面取得了良好的应用效果。在最小二乘偏移(LSM)中，地震偏移成像结果通常通过最小化数据残差的L2范数获得，该残差表征了正演模拟数据与实际观测数据之间的差异[1]。LSM的思想首先应用于Kirchhoff偏移[2-3]，然后推广到单程波动方程偏移[4-5]。目前，该思想已应用于双程波动方程，形成了最小二乘逆时偏移(LSRTM)[6-11]。与传统的RTM 成像[12-14]相比，LSRTM有助于生成高分辨率、高保真度和较少低频噪声的偏移成像结果[15-17]。

一般来说，常规的偏移方法利用地震正演算子的共轭转置而非真正意义的逆算子成像，导致由成像结果重建的地震数据与原始地震数据之间存在误差。LSM通过基于线性反演迭代的方法逼近逆算子，最终获得同预处理地震数据相匹配的地下构造模型，这种线性反演方法需要通过正演算子(或反偏移算子)和偏移算子建立地下构造模型与地震观测数据之间的线性关系[18-19]，正演—偏移算子对的伴随性是影响迭代收敛和最终偏移成像结果的重要因素。通常，正演算子通过Born近似理论导出[20-21]，偏移算子最常用的是RTM。然而，Xu等[22]证明了使用一阶声波方程时该算子对(Born正演和RTM)经过数值离散计算后，很难满足精确伴随关系[18]。因此，构建满足精确伴随关系的正演—偏移算子对是最小二乘逆时偏移中的一个关键步骤[23-25]，也是提高偏移质量的有效手段。

针对伴随算子对存在的问题，国内外学者主要采用两种途径进行研究：一是基于伴随状态法推导连续伴随波动方程，再谨慎进行离散获得波场延拓伴随算子，保证正演—偏移算子对的精确伴随性[26-27]；二是基于矩阵形式先进行离散后直接推导伴随矩阵的策略。通常采用点积测试[19]作为一种便捷的数值方法检验两个算子之间的伴随性，而点积测试最简单、直接的方法之一是利用算子矩阵的转置。不同于传统的叠后RTM算子，Ji[28]提出了一种适用于时域声波方程有限差分模拟的矩阵计算方法并推导了伴随矩阵； Yao等[29]基于矩阵形式对频率域LSRTM精确伴随算子对进行了研究； Xu等[22]详细描述了基于声波Born近似理论的精确伴随算子对的推导过程，并应用于预条件LSM。许多学者将这种基于Born正演的精确伴随算子对应用于声波[22，30-31]、弹性波[15，32]，甚至各向异性[33-34]、黏滞介质[7，9，35]的LSM。由于RTM是最常用的偏移算子，因此，另一种构造精确伴随算子对的方法是先通过RTM过程的矩阵表示，再将其转置，得到RTM反偏移算子(DeRTM算子)，DeRTM-RTM算子对同样有较好的成像效果[36]。此外，一些学者利用自伴随波动方程将这两种算子对混在一起，并展示了Born-RTM算子对在各种数值算例下的良好应用效果[30，37-39]。然而，这三组伴随算子对之间的区别与联系尚不清楚，需要进一步研究和探索。

为此，本文首先在Ji[28]的研究基础上，将时域二阶声波方程离散成矩阵形式，得到Born正演和RTM过程的矩阵表示。将矩阵转置即产生两组精确伴随算子对，分别表示为Born-AdjBorn算子对和RTM-DeRTM算子对。然后，基于声波方程的自伴随离散化，对这两个算子进行混合匹配，得到了第三组精确伴随算子对，即自伴随Born-RTM算子对。最后进一步推导了这三组精确伴随算子对之间的定量关系，并通过二维数值算例进行了验证。

1 方法原理

1.1 声波波场延拓及其伴随算子矩阵形式

1.1.1 波场延拓算子

一维时空域的二阶声波方程为

(1)

式中：p为压力场；c(z)为声波速度；t和z分别表示时间和深度。根据中心二阶有限差分格式，式(1)的离散形式可以写为

(2)

式中：n=1，2，…，nt，nt为时间方向样点数；i=1，2，…，nz，nz为深度方向样点数；α=Δt2/Δz2，其中Δt和Δz分别为时间间隔和深度网格间隔。将式(2)表示为矩阵形式

(3)

式中：pn=[p1,n,p2,n,…,pnz,n]T为n时刻波场的向量形式；I是单位矩阵；T是一个带状矩阵，具体为

(4)

Bp=s

(5)

(6)

地震记录可以通过观测矩阵D提取各时刻对应的波场。因此，观测数据可以表示为

d=Dp=DB-1s

(7)

二维波场延拓算子与一维类似，详细推导见附录A。

1.1.2 波场延拓伴随算子

如上所述，推导一对伴随算子最直接的方法之一是通过矩阵转置。因此，了解波场延拓算子B-1及其伴随算子B-T对于后续推导正演—偏移算子对有较大帮助。为了简单起见，只给出一维矩阵分析，对于高维可得出类似的结论。

应用线性代数方法，B-1可以显式表示为

B-1=I+A+A2+…+Ant

(8)

定义一系列矩阵算子{Pn}，n∈{1,2,…,nt}，将零时刻波场映射到n时刻波场。考虑式(6)中的定义，令P0=I，可得

(9)

因此，波场延拓p=B-1s可以写成子矩阵求和的形式

(10)

上式可以理解为，在n时刻的波场pn是由此前源项si通过Pn-i映射的所有波场的总和，这符合惠更斯原理和波场叠加原理。其伴随过程p=B-Ts，也可以表示为子阵求和形式

(11)

(12)

通过引入速度加权波场p(z,t)=c2(z)q(z,t)，可以将式(12)转化为常规波动方程

(13)

式(13)表明伴随波场q可以通过在常规波动方程中加载速度加权源项c2s获得。在速度模型不均匀的情况下，忽略速度加权因子将导致伴随波场的计算不准确，但这与二阶声波方程自伴随的性质不一致。为了保持空间离散中的自伴随性，引入另一个中间波场u(z,t)=c-1(z)p(z,t)，满足

(14)

1.2 三组精确伴随算子对

1.2.1 基于Born近似的精确伴随算子对

根据Born近似理论，在小速度扰动假设下，Born正演方程为

(15)

(16)

式中⊙表示Hadamard乘积。因此Born正演获得的观测数据可以表示为

(17)

Born正演的伴随算子是式(17)中正演算子FBorn的转置，即

(18)

1.2.2 基于RTM的精确伴随算子对

经典的RTM过程包括正向波场延拓、逆向波场延拓和成像。利用矩阵形式表示RTM全过程是推导基于RTM精确伴随算子对的关键。正向和逆向波场延拓可表示为

(19)

(20)

(21)

通过零延迟互相关成像条件，偏移过程可表示为

(22)

式中FRTM为逆时偏移算子。相应的反偏移过程可表示为

(23)

1.2.3 基于混合匹配的精确伴随算子对

为了保证混合匹配后算子对的伴随性，简单直接的方法是使用式(14)，而不是传统的声波方程，可以保证FBorn-s=FDeRTM-s和FAdjBorn-s=FRTM-s，即利用自伴随声波方程得到的两组伴随算子对Born-RTM和DeRTM-AdjBorn实际上是相同的。因此，在一般情况下，可以使用自伴随Born-RTM算子对作为LSM中的另一组精确伴随算子对。

综上所述，得到了LSM的三组精确伴随算子对，分别是利用常规离散波动方程的Born-AdjBorn算子对、利用常规离散波动方程的DeRTM-RTM算子对和利用自伴随离散波动方程的Born-RTM算子对。对应的LSM过程分别称为LSBM、LSRTM和自伴随LSBRTM。此外，考虑自伴随波场u、常规离散化波场p及其伴随波场q之间的关系，即u=c-1p=cq，可以得到三组精确伴随算子对成像结果之间的定量关系。

1.3 三组精确伴随算子对的定量关系

因此，在大多数采用常规离散波动方程的情况下，当且仅当同一源在均匀速度模型中延拓时，Born-AdjBorn算子对和DeRTM-RTM算子对是相同的。在此条件下，Born近似正演和RTM可以作为一组精确的伴随算子对。

其次，根据附录C中的推导，可以得出以下定量结论。

(1)FBorn=FBorn-s和FAdjBorn=FRTM-s。常规离散波动方程的Born-AdjBorn算子对与自伴随离散波动方程的Born-RTM算子对相同。因此，AdjBorn结果等于自伴随RTM结果，LSBM和自伴随LSBRTM的结果是相同的。

由于式(C-6)中的Hessian矩阵通常是高度不适定的，其逆矩阵很难通过数值计算。因此，在式(C-6)中的最小二乘解通常采用迭代优化算法计算，例如本文中使用的预条件共轭梯度法。然而，由于迭代优化算法难以收敛至极小值，导致数值解可能与最小二乘解并不完全一致，式(C-7)中LSM结果之间的定量关系较难满足。

1.4 基于精确伴随算子对的LSM算法

LSM的目标是速度扰动模型，该模型描述了观测数据所隐含的地下构造信息。因此求解过程通常包括正演算子、偏移算子与最小二乘算法。由于最小二乘算法在大多数LSRTM文献中都有讲述[27，40-41]，因此本文将重点关注精确伴随算子对的正确实现。

通过上述推导，得到了LSBM、LSRTM和自伴随LSBRTM伴随算子对的矩阵形式。然而，它们的矩阵运算是非常昂贵的，矩阵维度为nx×nz×nt。因此，需要将矩阵运算转化可实现的程序语言。表1总结了各算子、矩阵的实际意义。

LSM算法的目标是通过最小化目标函数求解速度扰动模型m

(24)

式中F表示前面推导出的三个正演算子之一。目标函数的第一项表示数据残差，第二项是由阻尼参数λ平衡的模型正则化项，并采用了预条件共轭梯度算法获得LSBM、LSRTM和自伴随LSBRTM成像结果。

表1 伴随算子对中各矩阵或算子的实际意义

2 数值算例

在数值算例中，应用Ricker子波、基于这些算子对的有限差分代码计算合成数据。由于所推导的算子对不考虑任何吸收边界条件，本文通过复制模型的外边界层扩展模型，并在数值边界反射波到达之前停止波场延拓。在预条件共轭梯度算法中，阻尼参数λ设置为0.001。

2.1 层状模型

自伴随LSBRTM、LSBM、LSRTM的一次迭代偏移结果分别如图2a～图2c所示，分别等价于自伴随RTM、传统Born偏移、传统RTM结果。LSM的30次迭代偏移结果如图3所示。

与一次迭代偏移结果相比，多次迭代能够以更高的分辨率、更均衡的振幅和更少的伪影刻画出贴近真实速度扰动模型的成像结果。在波速分布由低速变化至高速的界面位置处，其速度扰动模型由负值变为正值，在一次迭代偏移结果中界面刻画均为正值—负值—正值，在LSM结果中界面均收敛为接近真实速度扰动模型的负值—正值。

图1 层状模型(a)原始速度； (b)平滑速度； (c)速度扰动

图2 层状模型不同方法的一次迭代偏移结果(a)自伴随RTM； (b)AdjBorn； (c)RTM

图3 层状模型四种LSM方法的30次迭代结果(a)自伴随LSBRTM； (b)LSBM；(c)LSRTM； (d)速度加权LSRTM

为了验证三组伴随算子对结果之间的定量关系，计算了不同方法偏移成像结果的差值，如图4所示。图4a、图4b和图4d中极小的振幅证明了传统Born偏移结果(AdjBorn)与自伴随RTM的等价性、速度加权自伴随RTM与传统RTM之间的等价性以及LSBM与自伴随LSBRTM的等价性。对于自伴随LSBRTM与速度加权LSRTM结果的差异，由于存在误差，无法收敛至零，图4e难以证明其对应的定量关系。

此外，图4c中的差值剖面是一个没有明显伪影的层状模型图像。然而，在图2a和图2c中，自伴随RTM和传统RTM结果在地表附近都存在较强的低频伪影，其与速度变化和反射结构无关，而成像结果的差异只与速度变化有关，即目标反射体。因此，这可能是一个有效减少RTM成像结果中低频伪影的策略。

图4 层状模型不同方法偏移结果的差值剖面(a)传统Born偏移与自伴随RTM； (b)速度加权自伴随RTM与传统RTM； (c)传统RTM与自伴随RTM；(d)自伴随LSBRTM与LSBM； (e)速度加权LSRTM与自伴随LSBRTM

在均匀偏移速度模型(2500m/s)下测试了LSBM与LSRTM，偏移成像结果如图5所示。由于存在速度误差，第二个界面位置存在偏差，但两个偏移成像结果之间的差值为零，验证了本文的推论，即当偏移速度模型是均匀的时候，LSBM和LSRTM会呈现相同的成像结果。

在速度存在误差、数据含有噪声、数据有缺失的情况下，对三组精确伴随算子对进行敏感性分析。

应用平滑半径为100个网格距的偏移速度模型(图6a)，LSBM、LSRTM和自伴随LSBRTM成像结果如图7所示，可见，在速度误差的影响下，LSBM、LSRTM和自伴随LSBRTM均可对地下层状构造较准确成像。

在原始地震数据中添加白噪声，数据信噪比为-20dB(图6b)。三组精确伴随算子对的LSM成像结果如图8所示，受噪声的影响，均产生了较多的伪影。

图5 均匀偏移速度模型下两种方法的LSM结果(a)LSBM； (b)LSRTM

图6 低精度的偏移速度速度模型(a)、震源位于1000m的含噪声数据(b)及随机缺失87.5%的道集数据(c)

每炮地震数据随机缺失175道，即缺失87.5%(图6c)，三组精确伴随算子对的LSM成像结果如图9所示。数据缺失对LSBM、LSRTM和自伴随LSBRTM偏移结果影响一致，均会在浅层产生些许低频噪声。

在速度存在误差、数据含有噪声、数据有缺失的情况下，LSBM与自伴随LSBRTM之间差值极小，在10-4数量级； LSBM与LSRTM、自伴随LSBRTM与LSRTM之间的差值较大，与成像结果的幅值在同一数量级。因此，速度误差、数据噪声以及数据缺失对三组精确伴随算子对之间的定量关系基本没有影响。

2.2 Marmousi模型

应用更复杂的Marmousi模型(图10)测试三种的精确伴随算子对。模型网格数为1000×300，网格间距为10m。50个震源均匀分布在模型表面，震源主频为10Hz。检波器以10m间距均匀分布在模型表面。一次迭代偏移结果如图11所示，经过10次迭代，最终的LSM结果如图12所示。

图7 低精度速度模型的三种方法的LSM结果(a)LSBM； (b)LSRTM； (c)自伴随LSBRTM

图8 数据含有噪声时三种方法的LSM结果(a)LSBM； (b)LSRTM； (c)自伴随LSBRTM

图9 数据有缺失时三种方法的LSM结果(a)LSBM； (b)LSRTM； (c)自伴随LSBRTM

与层状模型结果类似，图11和图12中所有偏移结果符合预期，LSM结果展现了更少的低频噪声和更清晰的界面信息。图12中红色矩形和箭头突出了成像结果之间的差异，通过与速度扰动模型进行对比可以看出，LSBM与自伴随LSBRTM更关注浅层的精细成像； LSRTM相较于LSBM和自伴随LSBRTM对深层(2km以下)成像效果更好，红

图10 Marmousi模型(a)真实速度； (b)平滑速度； (c)速度扰动

图11 Marmousi模型不同方法的一次迭代偏移结果(a)自伴随RTM； (b)AdjBorn； (c)RTM

色矩形中的地质构造更加清晰；并且速度加权项LSRTM对深层的照明更强，红色箭头处的地层信息与速度扰动模型基本一致。说明在Marmousi模型算例中，显式加入速度加权项后会增强深层的成像效果。

在定量关系上，图13中极小的振幅分别验证了一次迭代偏移结果和LSM结果相应的定量关系，与层状模型算例类似，由于存在迭代误差，难以证明速度加权LSRTM与LSBM(或自伴随LSBRTM)之间的定量关系。图14的数据收敛曲线表明三组伴随算子对的收敛速度快且相似，经过10次迭代后，三种LSM方法的归一化数据残差都小于0.001。

图12 Marmousi模型四种LSM方法的10次迭代结果(a)自伴随LSBRTM； (b)LSBM； (c)LSRTM； (d)速度加权LSRTM

图13 Marmousi模型不同方法偏移结果的差值剖面(a)传统Born偏移和自伴随RTM； (b)速度加权自伴随RTM和传统RTM； (c)自伴随LSBRTM与LSBM

图14 Mairmousi模型三种LSM方法的数据残差对比(a)第15号炮； (b)第35号炮； (c)所有50炮

3 讨论

一般认为Born正演和RTM分别是LSM的正演和偏移算子。然而，在数值离散化后，它们可能无法准确地通过点积测试。当采用非精确伴随算子对进行LSM时，会导致收敛性差或累积误差较大的偏移成像结果。因此，在声波方程矩阵表达式和矩阵运算的基础上，本文从理论上推导了基于Born近似理论和RTM过程的三组精确伴随算子对，并找出了它们之间的等价条件和定量关系。同时，开发了一套无矩阵代码，通过数值算例验证该推论。

正如在精确伴随算子对推导中所证明的，LSBM和LSRTM成像结果之间的主要区别源于相应算子对中的波场延拓矩阵B-1和B-T，即式(1)及其伴随波动方程(式(14))。后者可以进一步扩展为

(25)

式(25)意味着算子B-T实际上是传统的声波方程B-1加上与速度空间扰动有关的两项。如果偏移速度模型为均匀速度模型，即∂c2/∂z=0、∂2c2/∂z2=0，则LSBM和LSRTM会得到相同的结果。因此，可以从不同的角度对Born-AdjBorn和DeRTM-RTM伴随算子对进行对比，并得到相同的结论。并且，在二阶声波方程中(以一维情况为例)，若将速度项移至方程左侧，即变为慢度项，即

(26)

将该声波方程进行离散，并加载震源项s，其结果与加载速度加权源项c2s的式(13)相似。这意味着，在LSM过程中，将慢度扰动模型作为收敛目标即可利用该声波方程的自伴随性质进行波场延拓模拟。因此，在最小二乘框架下的LSM方法中，波动方程的选择及模型参数化对于迭代收敛性能和最终的偏移成像结果影响较大。

此外，不论对于连续波动方程还是波动方程的离散形式，由于Born-AdjBorn和DeRTM-RTM算子对都是与一次反射波对应的线性化算子，实际观测数据中可能包含更复杂的波场分量，如多次波或绕射波等。因此，通过正演模拟(Born正演和DeRTM过程)很难得到与实际观测数据准确匹配的合成数据，所以不应追求获得非常小的数据残差，特别是对于复杂的Marmousi模型。这使得验证LSRTM与LSBM(或自伴随LSBRTM)结果之间的定量关系(式(C-7))变得更加困难。

在实际地震勘探过程中，地震观测数据存在较严重的噪声，且实际子波与波场延拓模拟子波之间存在差异，正演地震数据与实际观测数据本身具有一定误差。若采用非精确伴随算子对进行最小二乘偏移，会引入新的误差，且在迭代收敛过程中容易累计误差，难以保证迭代收敛，导致实际资料偏移结果存在伪影、假象。因此，在偏移速度模型较为准确的前提下，采用精确伴随算子对进行最小二乘偏移成像可有效提升实际地震资料的成像质量。

4 结论

最小二乘偏移需要正演和偏移两个过程实现成像结果的迭代更新，并且正演和偏移过程需满足精确伴随关系。基于Born近似理论和RTM过程，本文首先在常规离散化声波方程的基础上推导出了两对理论伴随算子。采用自伴随方法进行波动方程离散时，这两组算子对会产生另一组精确伴随算子对。这三组算子对均可通过无矩阵算法编程实现，并能准确通过点积测试。此外，还推导了这三组伴随算子对之间的定量关系，并应用二维模型进行了验证。定量关系为：采用同一源在均匀速度模型中进行波场延拓时，Born-AdjBorn算子对和DeRTM-RTM算子对是相同的；速度加权的LSRTM结果等于LSBM和自伴随LSBRTM的结果。

附录A 二维声波波场延拓算子

对于时空域二维声波方程

(A-1)

假设Δz=Δx，其有限差分格式为

pi+1,j,n+pi,j-1,n+pi,j+1,n)-pi,j,n-1

(A-2)

将波场表示为矩阵形式

pn=[p1,1,n,p2,1,n,…,pnz,1,n,p1,2,n,p2,2,n,…,

pnz,2,n,…,p1,nx,n,p2,nx,n,…,pnz,nx,n]T

(A-3)

式中nx、nz分别为x和z方向的样点数。在二维情况下，矩阵T变成了大小为(nx×nz)×(nx×nz)的三对角块矩阵

(A-4)