基于上下行波分解的分数阶拉普拉斯算子黏滞声波逆时偏移成像方法

王喜鹏，张庆臣，毛伟建*

（1.中国科学院精密测量科学与技术创新研究院计算与勘探地球物理研究中心，湖北武汉 430077；2.大地测量与地球动力学国家重点实验室，湖北武汉 430077； 3.中国科学院大学，北京 100049）

0 引言

近年来，基于双程波的逆时偏移成像技术突破了成像角度的限制，对陡倾构造成像具有显著优势［1-2］。但地震波在传播过程中，由于地下介质通常具有黏滞性，导致其能量衰减及相位畸变，尤其当勘探目标区域存在气云等强衰减构造时，对其做偏移成像会导致分辨率降低和照明不足［3］。为此，学者们进行了大量的研究，模拟地震波在地下介质传播过程中的衰减效应。

起初，对黏性逆时偏移的研究是基于近似常Q模型［4］波动方程展开的，通过单个或多个力学元件的串、并联或叠加（如标准线性体、Maxwell体等）描述地下介质的衰减效应［5-8］，但基于该模型推导的波动方程中包含记忆变量，需较多表征参数，这些参数还需要计算优化选择，额外增加计算量［9］。此外，在这类波动方程中，振幅衰减与相位频散相耦合，在逆时偏移成像过程中不能准确地校正相位，会导致深部目标层成像的位置误差。基于常Q模型［4］推导的波动方程中含有分数阶时间导数，对波动方程计算当前时刻波场值时，需依赖该时刻以前所有时刻的波场，因此其计算成本高、效率低且内存消耗巨大。

Chen 等［10］采用分数阶拉普拉斯算子代替分数阶时间导数，将时间导数转换为空间导数，运用伪谱法在波数域求解黏声波动方程，提高了计算效率。Treeby等［11］基于幂律衰减模型，推导了解耦的分数阶拉普拉斯算子黏滞声波波动方程，实现了振幅衰减与相位频散的分离。刘文卿等［12］基于GPU/CPU 协同系统，通过GPU 实现波场逆时外推，并采用随机速度边界提高算法的并行性，解决了逆时偏移中大规模存储的I/O 问题。Zhu 等［13］基于常Q频散关系，推导了含解耦的分数阶拉普拉斯算子时域黏声波动方程，在波传播过程中具有近似常Q的衰减和频散特征，实现了振幅衰减与相位频散的分离，在逆时偏移过程中对能量补偿、相位校正十分便利，只需在波场外推过程中，反转振幅衰减算子的符号且保持相位频散算子符号不变，就可较准确地校正由地下介质黏滞性引起的振幅衰减与相位频散效应，并基于此方程实现了Q补偿逆时偏移（Q-compensated reverse-time migration，Q-RTM）［14］。

然而，采用Q-RTM 在对地震波振幅进行衰减补偿时，也会增强数据中的高频噪声，在波场外推过程中出现数值不稳定现象。此前，高频噪声通常采用低通Tukey 滤波器压制［14-15］。但该方法是时不变的，不能适应在空间变化的Q值和深度。

Treeby［16］提出一种频域时变Tukey 滤波器以稳定衰减补偿过程，但采用低通滤波法会损失地震数据中部分有用高频信号。Wang［17］通过引入稳定因子，提出一种反Q偏移方法稳定补偿过程，并认为某深度处的地震波吸收衰减效应应是从地表到当前深度吸收衰减的累积； Wang 等［18］在此基础上提出一种新的Q补偿逆时偏移自适应稳定方法，并通过推导分数阶拉普拉斯算子黏滞声波波动方程及其补偿方程和衰减方程的波数域格林函数，从理论上解释了波场外推在高频端不稳定的原因。与传统基于低通滤波的方法相比，该方法具有更强的时间不一致性和Q依赖性，提高了衰减补偿稳定性和抗噪性。陈汉明等［19］应用最小二乘反演理论，推导了分数阶拉普拉斯算子黏滞声波方程对应的Born正演模拟算子和伴随方程，利用反演思路补偿地震波的吸收衰减，解决了补偿不稳定问题。Chen等［20］提出一种隐式稳定策略，将时变滤波器隐式地纳入振幅增强波场外推步骤，并在频谱中加入一个稳定因子，避免了Q-RTM 中任何显式的滤波操作，且该方法不增加额外计算量。

对于逆时偏移成像产生的低频噪声，目前主要采用如角度衰减因子成像条件［21-22］、对成像结果进行滤波（最小二乘滤波器、拉普拉斯滤波等）［23-24］、波场方向分解［25-26］等方法压制。杜启振等［27］总结了各类成像方法，认为采用波场分解的成像方法在复杂构造区能更好地消除低频噪声。Liu 等［25-26］提出上下行波场分解的方法，构造新的偏移成像条件，实现不同向的震源波场和检波点波场的互相关成像，从而消除低频噪声。Fei等［28］指出速度模型中存在强速度梯度或反射界面时，震源下行波场和检波点上行波场在逆时偏移成像过程中可能产生假象，成像时需略去这两项。

此外，常规波场分解是在f-k域实现，这需将震源波场和检波点波场进行存储，然后再进行时空域的傅里叶变换，因此会产生巨大的计算量和I/O 量。针对该问题，胡江涛等［29］提出一种基于声波波动方程的解析时间波场外推方法，并将其运用于VSP（Vertical Seismic Profile）数据的逆时偏移成像。该方法在时间外推过程中实现了波传播方向的显式分解，避免了巨大的I/O 量和计算量，实现了低频噪声与有效信号的分离，并能有效地消除图像中由复杂介质引起的假像。Zhou 等［30］采用一步法波场外推实现了Q补偿逆时偏移波场外推，并基于该方法波场为解析信号的性质实现了波场的分离，运用Q补偿逆时偏移和分解的波场成像极大程度提高了偏移成像的质量，但该方法采用的是传统低通滤波器稳定化衰减补偿过程，且一步法外推的计算速度较慢。

目前，Q-RTM 主要采用成像后滤波的方式压制低频噪声，但该方法会改变成像剖面中同相轴振幅大小和波数成分。基于波场分解的成像方法能更有效地压制双程波逆时偏移成像中的低频噪声和速度模型中存在较强速度梯度时产生的假象。

综上所述，本文采用基于常Q模型解耦的分数阶拉普拉斯算子黏声波动方程对震源波场和检波点波场进行外推，避免由于时间历史记忆变量而需存储每个时刻的波场，实现衰减算子和频散算子的分离，可通过反转衰减算子符号和保持频散算子符号不变达到补偿能量和校正频散的目的，并基于该方程实现Q补偿逆时偏移算法。

采用自适应稳定因子方案稳定化衰减补偿过程中数值不稳定问题。对于逆时偏移产生的低频噪声和假象，采用解析时间波场外推方法，在时间外推过程中实现波传播方向的显式分解，选取分解后的震源下行波场和检波点上行波场进行互相关成像，减少由复杂介质及强速度梯度引起的偏移成像中的低频噪声和假象。

1 方法原理

1.1 常Q 模型应力—应变关系

在一维线性非弹性介质中，应力与应变的关系［31］为

式中：ψ1-2γ为弛豫函数；t为时间；t0=1/ω0为参考时间；γ=arctan(1/Q)/π 为无量纲参数，且0＜γ＜1/2，Q值越小介质的衰减效应越强，当Q→∞时，γ→0，方程退化为声波方程；为体积模量，其中ρ0为密度，c0为参考频率ω0处的相速度；*为卷积算子；ε为应变场。

弛豫函数定义［9］为

式中：Γ为欧拉伽马函数；H(t)为Heaviside 阶跃函数。

1.2 解耦的分数阶拉普拉斯算子波动方程

对于均匀声波介质，一阶线性动量守恒方程为

应变速度方程为

两式中：σ为应力场；v为速度。

结合式（3），可得时间分数阶波动方程

式中： ∇2为拉普拉斯算子；。由于时间分数阶导数的存在，计算某时刻波场时需存储该时刻之前每一时刻的波场值，内存需求较大，且数值计算困难。为此，Chen等［10］提出用分数阶拉普拉斯算子(-∇2)α/2代替分数阶时间导数∂α/∂tα，分数阶拉普拉斯算子可在空间域计算，避免较大的内存需求及数值计算难题。

将平面波解exp(iωt-ik·r)代入式（5），可得对应的频散方程

式中：k=kR+ikI为复波数矢量；ω为角频率；r为空间坐标矢量。

由i2γ=cos(πγ)+i sin(πγ)，复波数代入k，经一系列计算及近似，可得新的频散方程

分数阶拉普拉斯算子的傅里叶变换为

式中：k为空间域波数； F表示傅里叶变换。

对式（7）进行时空二维反傅里叶变换，可得

1.3 Q-RTM 原理

逆时偏移有三个主要步骤： ①震源波场外推；②检波点波场外推； ③采用适当的成像条件对外推的震源和检波点波场进行成像，如互相关成像条件［32］

式中：S(t，x)为震源波场；R(t，x)为检波点波场；Tmax为波场外推的最大时间值。

Zhu等［14］基于式（9）实现Q补偿逆时偏移，并给出统一的震源及检波点波场外推的黏声波动方程

式中：p为压力场；c0为声波速度；L和H分别为分数阶拉普拉斯算子； {ηL-∇2}p和Hp分别为频散和振幅衰减算子；L=(-∇2)γ+1；H=(-∇2)γ+1/2。该方程中β1和β2对应项分别为频散算子和振幅损失算子，在正向外推过程中，β2=1 表示振幅衰减，反向外推过程中，时间t被-t代替，衰减项符号反转β2=-1，表示放大振幅。只需保持频散算子符号不变（β1=1），就可校正频散效应。若β1和β2都为0，则式（11）退化为声波方程； 若两者都为1，该方程为黏声波方程。而补偿过程则通过反转衰减算子的符号，使β1=1且β2=-1，式（11）变换为

本文通过式（12）对震源波场和检波点波场进行外推，校正地下介质对地震波造成的能量衰减和频散效应，并采用规则网格有限差分法计算时间导数、伪谱法计算拉普拉斯算子。

1.4 自适应稳定因子

Wang等［18］提出一种新的Q补偿逆时偏移自适应稳定补偿策略，推导过程如下。

Q-RTM 的衰减时间传播算子可表示为

补偿时间传播算子可表示为

根据上两式，振幅衰减算子定义为

振幅补偿算子定义为

式（16）由于指数项的存在，它在时间上是发散的，因此定义如下稳定化的振幅补偿算子

式中σ2为稳定因子。

由于振幅补偿算子已体现在波场外推过程中，则稳定因子可表示为

考虑当前时刻的衰减效应为零时刻到当前时刻衰减效应的累积，则每个时间步长稳定因子系数定义为

第n个时间步稳定因子系数为

1.5 上下行波场分解

由于采用双程波动方程，震源波场和检波点波场中同时包含上行波和下行波，运用成像条件使波路径上不必要的震源波场和检波点波场（震源上行波场与检波点上行波场，震源下行波场与检波点下行波场）互相关会在成像结果上产生低频噪声。此外，Fei等［28］指出当速度模型精度不高，存在较强速度梯度时，在特殊地质构造区，震源上行波场和检波点下行波场成像时可能会产生假象。

为避免上述情况，本文采用的成像条件为

式中： 下标u（upgoing）表示上行波； d（downgoing）表示下行波，如Sd则表示震源下行波场。

在f-k域通过简单的计算很容易分离上行波与下行波分量［33］

式中：S(ω，kz)为震源波场S(t，z)的二维傅里叶变换；kz为波数；Su(ω，kz)和Sd(ω，kz)分别为f-k域分解的上行波和下行波分量。由该式可知，波场方向由频率与波数乘积的正负号决定。

Liu 等［26］通过固定波数符号，将震源和检波点波场分为波数大于零和小于零的四个分量，使波数符号相反的震源波场与检波点波场相乘，再将相乘后的两个分量相加，最后使用互相关成像条件成像（如震源波场波数符号为正，则与之相乘的检波点波场波数符号为负。若频率为正，则分别为震源上行波场及检波点下行波场；若频率为负，则分别为震源下行波场和检波点上行波场），实现传播方向相反的震源和检波点波场的互相关成像。但该方法中上行波场与下行波场互相耦合，不能将波场进行显式地分解。

为了实现波场的显式分解，需尽量使用一个变量（频率或波数）的符号来确定波场传播的方向。由于在时间域波场外推中，时间步长上波数信息较易获取，若能获取一个其频谱仅为正或负的波场，则波场的方向分解就可通过波数的正、负号确定。为此，胡江涛等［29］提出一种解析时间波场外推方法。解析波场频谱只包含正频率，可仅依靠空间波数的符号实现波场方向的分解。解析波场包括实部和虚部两部分，实部由信号本身构成，虚部由经希尔伯特变换后的原信号构成。

震源的时间解析波场可表示为

基于式（22）和解析波场仅包含正频率的性质，只需根据波数的符号就可达到对波场进行方向分解的目的。震源的上、下行波场表示为

式中： Re表示实部； F-1表示傅里叶反变换；表示解析时间震源波场关于z方向的傅里叶变换。同理，检波点波场的上、下行波表示为

值得注意的是，由于采用解析波场外推分离上下行波场，与传统逆时偏移成像方法相比多了震源波场外推、检波点波场外推及时间切片上的波场傅里叶变换，因此该方法计算时间会高于传统Q补偿逆时偏移成像。实验表明，本文方法的计算时间约为采用拉普拉斯滤波Q-RTM 方法的1.3倍。但与传统的波场分解相比，该方法规避了I/O 量过大的弊端，解决了波场存储问题，且提高了计算效率。

2 数值算例

2.1 两层模型

为了测试本文方法的有效性，首先采用水平层状模型对解析波场分离上下行波的有效性进行测试。图1a 所示速度模型中，第一、第二层层速度分别为2000、3000 m/s，模型网格尺寸为300×300，x和z方向的网格间距为dx=dz=10 m，得到该层状模型对应的Q模型（图1b）。图2a为1.0 s时刻未分离的波场快照，用解析波场分解方法得到该时刻经波场分离后的下行波场（图2b）和上行波场（图2c）。可见本文方法能有效分离上下行波场。

图1 两水平层速度模型（a）及其Q 模型（b）

图2 1.0 s 时刻波场快照

2.2 抽稀的Sigsbee 模型

为验证采用的新成像方法与Q补偿逆时偏移对振幅补偿和相位校正的有效性及对复杂模型的适应性，采用抽稀后的Sigsbee 模型做数值模拟测试。速度模型（图3a）及对应的Q模型（图3b，由经验公式得出）的大小为1067×401，网格间距为dx=dz=10 m，时间采样间隔为1 ms，采样点数为4001，震源采用主频为30 Hz的Ricker子波，总共104 炮，均匀分布在深度50 m 处，炮间距为100 m。

图3 Sigsbee 速度模型（a）及其Q 模型（b）

图4a 为2.0 s 时刻未分离的波场快照，用解析波场分解方法分别得到该时刻经波场分离后的下行波场（图4b）和上行波场（图4c）。

图4 2.0 s 时刻波场快照

图5a 为未经波场分离的传统声波逆时偏移成像结果，分别得到波场分离后的声波逆时偏移成像结果（图5b，参考）、波场分离后的未补偿衰减数据逆时偏移成像结果（图5c）及Q补偿逆时偏移成像结果（图5d）。图5b右侧子图为选取水平距离8000 m 处作为参考的单道信息（黑色实线），图5c、图5d 右侧子图则分别为水平距离8000 m 处（图5c虚线位置）相应的单道信息和参考单道信息（黑色实线）。图6为从上述子图中截取的0.6～1.6 km 层段单道信息。其对应的图5b～图5d的振幅谱如图7a～图7c所示。

图5 偏移成像结果

图6 截取的0.6～1.6 km 层段单道信息

图7 不同逆时偏移的振幅谱

与图5b（声波）相比，图5c（衰减）偏移成像层位及盐丘下边界能量减弱；通过对应的单道信息子图与振幅谱对比，与作为参考的声波数据相比，衰减数据的振幅减小，相位发生畸变，频带变窄。而图5d（Q补偿逆时偏移）中层位及盐丘下边界更清晰； 通过与参考图像（5b）截取的子图及振幅谱对比，可见其损失的振幅与产生的相移都得到一定程度校正，与参考图像吻合度较高，频带得到拓宽。从图5a与图5b 可知，采用传统的逆时偏移成像方法在成像剖面上存在严重的噪声及假象，降低成像精度，而采用上、下行波场分离成像的偏移结果则避免了这种影响，提高了分辨率。

2.3 抽稀的Pluto 模型

为进一步验证Q补偿逆时偏移对深层能量补偿的效果及与传统的拉普拉斯滤波方法对偏移成像噪声压制的效果，采用抽稀后的Pluto 模型进行数值模拟测试。其速度模型（图8a）及对应的Q模型（图8b，由经验公式得出）的大小为1160×201，网格间距为dx=dz=10 m，时间采样间隔为0.5 ms，采样点数为4002，震源采用主频为30 Hz的Ricker子波，总共114炮，炮间距为100 m，均匀分布在深度50 m 处。

图8 Pluto 速度模型（a）及其Q 模型（b）

图9a 为1.0 s 时刻未分离的波场快照，用上述方法分别得到该时刻经波场分离后的下行波场（图9b）和上行波场（图9c）波场快照。

图9 1.0s 时刻波场快照

图10a为未经波场分离的传统声波逆时偏移成像结果，分别得到声波逆时偏移经拉普拉斯滤波后成像结果（图10c）、Q补偿逆时偏移经拉普拉斯滤波后成像结果（图10e）、波场分离后的声波逆时偏移成像结果（图10b，参考）、波场分离后的衰减数据逆时偏移成像结果（图10d）及波场分离后的Q补偿逆时偏移成像结果（图10f）。图10b右侧子图为水平距离6000 m 处作为参考的单道信息（黑色实线），图10d、图10f 右侧子图则分别为水平距离6000 m 处相应的单道信息（红色虚线）和参考单道信息（黑色实线）。图11a～图11c对应为10b、图10d、图10f 中水平距离5～8 km、深度1.3～1.8 km 范围的振幅谱。

图10 不同的逆时偏移成像结果

图11 图10d 中红框对应位置振幅谱

与图10b（参考）相比，图10d（衰减未补偿）因地层的衰减效应，其能量减弱，偏移成像层位不清晰，尤其是盐下边界以下； 图10f 由于采用Q补偿逆时偏移，其能量得到补偿。对比10b、图10d、图10f 的单道信息与图11a～图11c 的振幅谱，可见缘于地层的非完全弹性效应，对成像结果造成振幅衰减及相位畸变，而采用Q补偿逆时偏移（图11c）能在一定程度上减缓这种效应的影响，拓宽其振幅谱，补偿衰减的能量，并对畸变的相位进行校正，提高了成像质量。由图10a可知，采用传统的互相关偏移方法，成像结果中浅层的有效信息完全湮没在噪声中，严重干扰成像精度；而采用拉普拉斯滤波后的成像结果（图10c、图10e）中低频噪声受到一定程度压制，但对浅层信息仍有较明显的干扰。

综上所述，本文采用的Q补偿逆时偏移方法能有效补偿由地下介质非弹性特性造成的振幅损失和相移，采用该成像方法所得偏移成像结果（图10d～图10f）中，有效地压制由传统成像条件所引起的低频噪声及假象。

值得注意的是，该方法是在进行两次震源波场和检波点波场外推的情况下，还要在时间切片上对波场进行傅里叶变换，所以计算时间稍长于采用拉普拉斯滤波的Q-RTM，测试结果（表1）显示约为其1.3 倍。但该方法不用对波场进行存储，避免了巨大的I/O量。

表1 不同方案下的计算时间

3 结论

本文提出一种上、下行波分解的解耦分数阶拉普拉斯算子黏滞声波逆时偏移成像方法：采用解耦的分数阶拉普拉斯算子黏滞声波方程对震源波场和检波点波场进行外推，能显著减缓由于地下介质的黏滞性所引起的振幅衰减和频散效应，提高了偏移成像的照明能量；采用分数阶拉普拉斯算子代替分数阶时间导数，不需在求解波动方程计算当前时刻波场值时依赖该时刻以前所有时段的波场。采用解析时间波场能显式地分解上下行波场。新的成像条件能避免同向传播的震源波场和检波点波场在应用互相关成像条件时所产生的低频噪声，消除由于偏移速度模型存在较强速度梯度时由震源上行波场和检波点下行波场互相关产生的假象，从而提高偏移成像的分辨率和质量。