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

2023 年 8 月第 58 卷 第 4 期基于上下行波分解的分数阶拉普拉斯算子黏滞声波逆时偏移成像方法王喜鹏1，2，3，张庆臣1，2，毛伟建*1，2（1.中国科学院精密测量科学与技术创新研究院计算与勘探地球物理研究中心，湖北武汉 430077；2.大地测量与地球动力学国家重点实验室，湖北武汉 430077；3.中国科学院大学，北京 100049）摘要：地下介质通常存在黏滞性，使地震波在传播过程中产生能量衰减及（相位）频散效应，降低了偏移成像照明能量。采用互相关成像条件的双程波动方程逆时偏移技术，当速度模型精度不高或存在较强速度梯度时，可能会产生低频噪声和假象。为此，提出一种基于上下行波分解的解耦分数阶拉普拉斯算子黏滞声波逆时偏移成像（Q 补偿逆时偏移）方法，采用解耦的分数阶拉普拉斯算子黏滞声波方程分别对震源波场和检波点波场进行外推，在此过程中校正地下介质对地震波产生的衰减效应；引入自适应稳定因子应对地震波衰减补偿过程中数值不稳定问题；采用时间解析波场外推方法，在时间切片上分别对震源、检波点波场的传播方向进行分解；最终运用互相关成像条件对下行的震源波场和上行的检波点波场进行成像。数值测试结果表明，所采用的 Q 补偿逆时偏移方法可校正由地下介质黏性特性引起的地震波振幅衰减和（相位）频散效应，该成像方法能显著减少传统逆时偏移成像中的低频噪声与假象。关键词：Q补偿逆时偏移，波场分解，低频噪声，分数阶中图分类号：P631 文献标志码：A doi：10.13810/ki.issn.1000-7210.2023.04.012Viscoacoustic reversetime migration imaging based on fractionalLaplacian with upgoing and downgoing wave decompositionWANG Xipeng1，2，3，ZHANG Qingchen1，2，MAO Weijian1，2（1.Center for Computational and Exploration Geophysics，Innovation Academy for Precision Measurement Science and Technology，Chinese Academy of Sciences，Wuhan，Hubei 430077，China；2.State Key Laboratory of Geodesy and Earths Dynamics，Wuhan，Hubei 430077，China；3.University of Chinese Academy of Sciences，Beijing 100049，China）Abstract：The subsurface media typically exhibit viscosity，which causes energy attenuation and（phase）dispersion effects during seismic wave propagation and in turn reduces the illumination energy of migration imaging.The twoway wave equation reversetime migration with crosscorrelation imaging conditions may produce lowfrequency noise and false images when there is no highly accurate velocity models or strong velocity gradients.To this end，this paper proposes a viscoacoustic reversetime migration method based on decoupled fractional Laplacian with up and downgoing wave decomposition（Q compensated reversetime migration method）.The viscoacoustic wave equation based on decoupled fractional Laplacian is employed to separately extrapolate the source wavefield and receiver wavefield and correct the attenuation effects generated by subsurface media during the process.An adaptive stability factor is introduced to handle numerical instability arising from the attenuation compensation of seismic waves.The timedomain analytical wavefield extrapolation method is adopted to separately decompose the propagation directions of the source and receiver wavefields on each time slice.Finally，the crosscorrelation imaging condition is utilized for imaging of the downgoing source wavefield and upgoing receiver wavefield.Numerical tests show that the Qcompensated reversetime migration method can correct the amplitude attenuation and（phase）dispersion effects caused by the viscosity characteristics of subsurface media.This imaging method significantly reduces the lowfrequency noise and false images in traditional reversetime migration imaging.偏移成像 文章编号：1000-7210（2023）04-0872-11*湖北省武汉市洪山区徐东大街340号中国科学院精密测量科学与技术创新研究院计算与勘探地球物理研究中心，430077。Email： 本文于 2022年 9月 17日收到，最终修改稿于 2023年 5月 15日收到。本项研究受国家自然科学基金重点项目“海底四分量高精度地震偏移成像方法研究”（42130808）资助。第 58 卷 第 4 期王喜鹏，等：基于上下行波分解的分数阶拉普拉斯算子黏滞声波逆时偏移成像方法Keywords：Qcompensated reversetime migration，wavefield decomposition，lowfrequency noise，fractional order王喜鹏，张庆臣，毛伟建.基于上下行波分解的分数阶拉普拉斯算子黏滞声波逆时偏移成像方法 J.石油地球物理勘探，2023，58（4）：872882.WANG Xipeng，ZHANG Qingchen，MAO Weijian.Viscoacoustic reverse time migration imaging based on fractionalLaplacian with upgoing and downgoing wave decomposition J.Oil Geophysical Prospecting，2023，58（4）：872882.0引言近年来，基于双程波的逆时偏移成像技术突破了成像角度的限制，对陡倾构造成像具有显著优势12。但地震波在传播过程中，由于地下介质通常具有黏滞性，导致其能量衰减及相位畸变，尤其当勘探目标区域存在气云等强衰减构造时，对其做偏移成像会导致分辨率降低和照明不足3。为此，学者们进行了大量的研究，模拟地震波在地下介质传播过程中的衰减效应。起初，对黏性逆时偏移的研究是基于近似常 Q 模型4 波动方程展开的，通过单个或多个力学元件的串、并联或叠加（如标准线性体、Maxwell体等）描述地下介质的衰减效应58，但基于该模型推导的波动方程中包含记忆变量，需较多表征参数，这些参数还需要计算优化选择，额外增加计算量9。此外，在这类波动方程中，振幅衰减与相位频散相耦合，在逆时偏移成像过程中不能准确地校正相位，会导致深部目标层成像的位置误差。基于常Q模型4 推导的波动方程中含有分数阶时间导数，对波动方程计算当前时刻波场值时，需依赖该时刻以前所有时刻的波场，因此其计算成本高、效率低且内存消耗巨大。Chen等10 采用分数阶拉普拉斯算子代替分数阶时间导数，将时间导数转换为空间导数，运用伪谱法在波数域求解黏声波动方程，提高了计算效率。Treeby等11 基于幂律衰减模型，推导了解耦的分数阶拉普拉斯算子黏滞声波波动方程，实现了振幅衰减与相位频散的分离。刘文卿等12 基于 GPU/CPU 协同系统，通过GPU实现波场逆时外推，并采用随机速度边界提高算法的并行性，解决了逆时偏移中大规模存储的 I/O问题。Zhu等13 基于常 Q频散关系，推导了含解耦的分数阶拉普拉斯算子时域黏声波动方程，在波传播过程中具有近似常Q的衰减和频散特征，实现了振幅衰减与相位频散的分离，在逆时偏移过程中对能量补偿、相位校正十分便利，只需在波场外推过程中，反转振幅衰减算子的符号且保持相位频散算子符号不变，就可较准确地校正由地下介质黏滞性引起的振幅衰减与相位频散效应，并基于此方程实现了Q补偿逆时偏移（Qcompensated reversetime migration，QRTM）14。然而，采用QRTM在对地震波振幅进行衰减补偿时，也会增强数据中的高频噪声，在波场外推过程中出现数值不稳定现象。此前，高频噪声通常采用低通 Tukey滤波器压制1415。但该方法是时不变的，不能适应在空间变化的Q值和深度。Treeby16 提出一种频域时变 Tukey滤波器以稳定衰减补偿过程，但采用低通滤波法会损失地震数据中部分有用高频信号。Wang17 通过引入稳定因子，提出一种反Q偏移方法稳定补偿过程，并认为某深度处的地震波吸收衰减效应应是从地表到当前深度吸收衰减的累积；Wang等18 在此基础上提出一种新的Q补偿逆时偏移自适应稳定方法，并通过推导分数阶拉普拉斯算子黏滞声波波动方程及其补偿方程和衰减方程的波数域格林函数，从理论上解释了波场外推在高频端不稳定的原因。与传统基于低通滤波的方法相比，该方法具有更强的时间不一致性和 Q 依赖性，提高了衰减补偿稳定性和抗噪性。陈汉明等19 应用最小二乘反演理论，推导了分数阶拉普拉斯算子黏滞声波方程对应的Born正演模拟算子和伴随方程，利用反演思路补偿地震波的吸收衰减，解决了补偿不稳定问题。Chen等20 提出一种隐式稳定策略，将时变滤波器隐式地纳入振幅增强波场外推步骤，并在频谱中加入一个稳定因子，避免了QRTM中任何显式的滤波操作，且该方法不增加额外计算量。对于逆时偏移成像产生的低频噪声，目前主要采用如角度衰减因子成像条件2122、对成像结果进行滤波（最小二乘滤波器、拉普拉斯滤波等）2324、波场方向分解2526 等方法压制。杜启振等27 总结了各类成像方法，认为采用波场分解的成像方法在复杂构造区能更好地消除低频噪声。Liu等2526 提出上下行波场分解的方法，构造新的偏移成像条件，实现不同向的震源波场和检波点波场的互相关成像，从而消除低频噪声。Fei等28 指出速度模型中存在强速度梯度或反射界面873石 油 地 球 物 理 勘 探2023 年时，震源下行波场和检波点上行波场在逆时偏移成像过程中可能产生假象，成像时需略去这两项。此外，常规波场分解是在fk域实现，这需将震源波场和检波点波场进行存储，然后再进行时空域的傅里叶变换，因此会产生巨大的计算量和I/O量。针对该问题，胡江涛等29 提出一种基于声波波动方程的解析时间波场外推方法，并将其运用于 VSP（Vertical Seismic Profile）数据的逆时偏移成像。该方法在时间外推过程中实现了波传播方向的显式分解，避免了巨大的I/O量和计算量，实现了低频噪声与有效信号的分离，并能有效地消除图像中由复杂介质引起的假像。Zhou等30 采用一步法波场外推实现了 Q补偿逆时偏移波场外推，并基于该方法波场为解析信号的性质实现了波场的分离，运用Q补偿逆时偏移和分解的波场成像极大程度提高了偏移成像的质量，但该方法采用的是传统低通滤波器稳定化衰减补偿过程，且一步法外推的计算速度较慢。目前，QRTM主要采用成像后滤波的方式压制低频噪声，但该方法会改变成像剖面中同相轴振幅大小和波数成分。基于波场分解的成像方法能更有效地压制双程波逆时偏移成像中的低频噪声和速度模型中存在较强速度梯度时产生的假象。综上所述，本文采用基于常Q模型解耦的分数阶拉普拉斯算子黏声波动方程对震源波场和检波点波场进行外推，避免由于时间历史记忆变量而需存储每个时刻的波场，实现衰减算子和频散算子的分离，可通过反转衰减算子符号和保持频散算子符号不变达到补偿能量和校正频散的目的，并基于该方程实现Q补偿逆时偏移算法。采用自适应稳定因子方案稳定化衰减补偿过程中数值不稳定问题。对于逆时偏移产生的低频噪声和假象，采用解析时间波场外推方法，在时间外推过程中实现波传播方向的显式分解，选取分解后的震源下行波场和检波点上行波场进行互相关成像，减少由复杂介质及强速度梯度引起的偏移成像中的低频噪声和假象。1方法原理1.1常 Q模型应力应变关系在一维线性非弹性介质中，应力与应变的关系31 为=1-2(t)*ddt=M0t-202t2（1）式中：1-2为弛豫函数；t为时间；t0=1/0为参考时间；=arctan(1/Q)/为无量纲参数，且 0 1/2，Q值越小介质的衰减效应越强，当Q时，0，方程退化为声波方程；M0=0c20cos2(/2)为体积模量，其中0为密度，c0为参考频率0处的相速度；*为卷积算子；为应变场。弛豫函数定义9 为1-2=M0(1-2)(tt0)-2H(t)（2）式中：为欧拉伽马函数；H(t)为 Heaviside 阶跃函数。1.2解耦的分数阶拉普拉斯算子波动方程对于均匀声波介质，一阶线性动量守恒方程为tv=10（3）应变速度方程为t=v（4）两式中：为应力场；v为速度。结合式（3），可得时间分数阶波动方程2-2t2-2=c2-202（5）式中：2为拉普拉斯算子；c2=M00=c20cos2(/2)。由于时间分数阶导数的存在，计算某时刻波场时需存储该时刻之前每一时刻的波场值，内存需求较大，且数值计算困难。为此，Chen等10 提出用分数阶拉普拉斯算子(-2)/2代替分数阶时间导数/t，分数阶拉普拉斯算子可在空间域计算，避免较大的内存需求及数值计算难题。将平面波解exp(it-ikr)代入式（5），可得对应的频散方程2c2=(i)2-202k2（6）式中：k=kR+ikI为复波数矢量；为角频率；r为空间坐标矢量。由i2=cos()+i sin()，复波数代入k，经一系列计算及近似，可得新的频散方程2c2=-202k2cos()+i-202k2sin()（7）874第 58 卷 第 4 期王喜鹏，等：基于上下行波分解的分数阶拉普拉斯算子黏滞声波逆时偏移成像方法分数阶拉普拉斯算子的傅里叶变换为(-2)/2(r，t)Fk(k，)（8）式中：k为空间域波数；F表示傅里叶变换。对式（7）进行时空二维反傅里叶变换，可得1c22t2=(-2)+1+t(-2)+1/2（9）式 中：=-c2-10-20sin()；=-c2-10-20 cos()。该方程为新的含分数阶拉普拉斯算子的黏声波动方程，将时间分数阶导数转换为分数阶拉普拉斯算子，并实现了振幅衰减与相位频散的解耦。1.3QRTM 原理逆时偏移有三个主要步骤：震源波场外推；检波点波场外推；采用适当的成像条件对外推的震源和检波点波场进行成像，如互相关成像条件32 I(x)=0TmaxS(t，x)R(t，x)dt（10）式中：S(t，x)为震源波场；R(t，x)为检波点波场；Tmax为波场外推的最大时间值。Zhu等14 基于式（9）实现Q补偿逆时偏移，并给出统一的震源及检波点波场外推的黏声波动方程 1c202pt2=2p+1L-2p+2tHp（11）式中：p为压力场；c0为声波速度；L和H分别为分数阶拉普拉斯算子；L-2p和tHp分别为频散和振幅衰减算子；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）变换为1c202pt2=2p+L-2p-tHp（12）本文通过式（12）对震源波场和检波点波场进行外推，校正地下介质对地震波造成的能量衰减和频散效应，并采用规则网格有限差分法计算时间导数、伪谱法计算拉普拉斯算子。1.4自适应稳定因子Wang等18 提出一种新的Q补偿逆时偏移自适应稳定补偿策略，推导过程如下。QRTM的衰减时间传播算子可表示为att(k，t)=sin1()k t e-2()k t1()k（13）补偿时间传播算子可表示为comp(k，t)=sin1()k t e2()k t1()k（14）上 两 式 中：1(k)=-2c4|k4+2-4c2|k2+22，2(k)=-c2|k|2+1/2，分别对应解耦的频散关系 式（7）解的实部和虚部。根据上两式，振幅衰减算子定义为(k，t)=e-2(k)t=e12c2-10-20sin()|k2+1t（15）振幅补偿算子定义为(k，t)=-1(k，t)=e2(k)t（16）式（16）由于指数项的存在，它在时间上是发散的，因此定义如下稳定化的振幅补偿算子 (k，t)=(k，t)2(k，t)+2=e2(k)t1+2e22(k)t（17）式中2为稳定因子。由于振幅补偿算子已体现在波场外推过程中，则稳定因子可表示为S(k，t)=11+2e22(k)t（18）考虑当前时刻的衰减效应为零时刻到当前时刻衰减效应的累积，则每个时间步长稳定因子系数定义为n=1ms(k，nt)=S(k，mt)（19）第n个时间步稳定因子系数为s(k，lt)=S(k，nt)Sk，(n-1)t=1+2e22(k)(n-1)t1+2e22(k)nt n=2，3，m1.5上下行波场分解由于采用双程波动方程，震源波场和检波点波场中同时包含上行波和下行波，运用成像条件使波路径上不必要的震源波场和检波点波场（震源上行波场与检波点上行波场，震源下行波场与检波点下行波场）（20）875石 油 地 球 物 理 勘 探2023 年互相关会在成像结果上产生低频噪声。此外，Fei等28 指出当速度模型精度不高，存在较强速度梯度时，在特殊地质构造区，震源上行波场和检波点下行波场成像时可能会产生假象。为避免上述情况，本文采用的成像条件为I(x)=0TmaxSd(t，x)Ru(t，x)dt（21）式中：下标u（upgoing）表示上行波；d（downgoing）表示下行波，如Sd则表示震源下行波场。在fk域通过简单的计算很容易分离上行波与下行波分量33 Sd(，kz)=S(，kz)kz 00 kz 0Su(，kz)=S(，kz)kz 00 kz 0（22）式中：S(，kz)为震源波场S(t，z)的二维傅里叶变换；kz为波数；Su(，kz)和Sd(，kz)分别为 fk域分解的上行波和下行波分量。由该式可知，波场方向由频率与波数乘积的正负号决定。Liu等26 通过固定波数符号，将震源和检波点波场分为波数大于零和小于零的四个分量，使波数符号相反的震源波场与检波点波场相乘，再将相乘后的两个分量相加，最后使用互相关成像条件成像（如震源波场波数符号为正，则与之相乘的检波点波场波数符号为负。若频率为正，则分别为震源上行波场及检波点下行波场；若频率为负，则分别为震源下行波场和检波点上行波场），实现传播方向相反的震源和检波点波场的互相关成像。但该方法中上行波场与下行波场互相耦合，不能将波场进行显式地分解。为了实现波场的显式分解，需尽量使用一个变量（频率或波数）的符号来确定波场传播的方向。由于在时间域波场外推中，时间步长上波数信息较易获取，若能获取一个其频谱仅为正或负的波场，则波场的方向分解就可通过波数的正、负号确定。为此，胡江涛等29 提出一种解析时间波场外推方法。解析波场频谱只包含正频率，可仅依靠空间波数的符号实现波场方向的分解。解析波场包括实部和虚部两部分，实部由信号本身构成，虚部由经希尔伯特变换后的原信号构成。震源的时间解析波场可表示为 -S(t，z)=S(t，z)+iSH(t，z)（23）式中；-S(t，z)为震源时间解析波场；S(t，z)为原始震源波场；SH(t，z)为震源波场经时间维希尔伯特变换的波场。同理，可得检波点的时间解析波场 -R(t，z)=R(t，z)+iRH(t，z)（24）式中：-R(t，z)为检波点时间解析波场；R(t，z)为原始检波点；RH(t，z)为检波点波场经时间维希尔伯特变换的波场。对于时间维希尔伯特变换的震源波场和检波点波场，则是通过分别对子波和地震记录进行时间维希尔伯特变换后由波场外推所获得。基于式（22）和解析波场仅包含正频率的性质，只需根据波数的符号就可达到对波场进行方向分解的目的。震源的上、下行波场表示为 Su(x，t)=ReF-1S+(x，kz)S+(x，kz)=S(x，kz)kz 00 kz 0Sd(x，t)=ReF-1S-(x，kz)S-(x，kz)=S(x，kz)kz 00 kz 0（25）式中：Re表示实部；F-1表示傅里叶反变换；S表示解析时间震源波场-S关于 z方向的傅里叶变换。同理，检波点波场的上、下行波表示为 Ru(x，t)=ReF-1R+(x，kz)R+(x，kz)=R(x，kz)kz 00 kz 0Rd(x，t)=ReF-1R-(x，kz)R-(x，kz)=R(x，kz)kz 00 kz 0（26）式中R表示解析时间检波点波场-R关于z方向的傅里叶变换。将式（26）中的震源下行波场和检波点上行波场代入式（21），即可得到最终的成像结果。值得注意的是，由于采用解析波场外推分离上下行波场，与传统逆时偏移成像方法相比多了震源波场外推、检波点波场外推及时间切片上的波场傅里叶变换，因此该方法计算时间会高于传统Q补偿逆时偏移成像。实验表明，本文方法的计算时间约为采用拉普拉斯滤波QRTM方法的1.3倍。但与传统的波场分解相比，该方法规避了 I/O 量过大的弊端，解决了波场存储问题，且提高了计算效率。876第 58 卷 第 4 期王喜鹏，等：基于上下行波分解的分数阶拉普拉斯算子黏滞声波逆时偏移成像方法2数值算例2.1两层模型为了测试本文方法的有效性，首先采用水平层状模型对解析波场分离上下行波的有效性进行测试。图 1a 所示速度模型中，第一、第二层层速度分别为2000、3000 m/s，模型网格尺寸为 300300，x和 z方向的网格间距为dx=dz=10 m，得到该层状模型对应的Q模型（图1b）。图2a为1.0 s时刻未分离的波场快照，用解析波场分解方法得到该时刻经波场分离后的下行波场（图 2b）和上行波场（图 2c）。可见本文方法能有效分离上下行波场。2.2抽稀的 Sigsbee模型为验证采用的新成像方法与 Q补偿逆时偏移对振幅补偿和相位校正的有效性及对复杂模型的适应性，采用抽稀后的 Sigsbee模型做数值模拟测试。速度模型（图 3a）及对应的 Q模型（图 3b，由经验公式得出）的大小为1067401，网格间距为dx=dz=10 m，时间采样间隔为 1 ms，采样点数为 4001，震源采用主频为30 Hz的Ricker子波，总共104炮，均匀分布在深度50 m处，炮间距为100 m。图 4a为 2.0 s时刻未分离的波场快照，用解析波场分解方法分别得到该时刻经波场分离后的下行波场（图4b）和上行波场（图4c）。图 5a为未经波场分离的传统声波逆时偏移成像结果，分别得到波场分离后的声波逆时偏移成像结果（图 5b，参考）、波场分离后的未补偿衰减数据逆时偏移成像结果（图 5c）及 Q 补偿逆时偏移成像结果 （图5d）。图5b右侧子图为选取水平距离8000 m处作为参考的单道信息（黑色实线），图 5c、图 5d右侧子图则分别为水平距离8000 m处（图5c虚线位置）相应的单道信息和参考单道信息（黑色实线）。图6为从上述子图中截取的 0.61.6 km 层段单道信息。其对应的图5b图5d的振幅谱如图7a图7c所示。与图5b（声波）相比，图5c（衰减）偏移成像层位及盐丘下边界能量减弱；通过对应的单道信息子图与振图 1两水平层速度模型（a）及其 Q模型（b）图 21.0 s时刻波场快照（a）原始波场；（b）下行波场；（c）上行波场图 3Sigsbee速度模型（a）及其 Q模型（b）877石 油 地 球 物 理 勘 探2023 年幅谱对比，与作为参考的声波数据相比，衰减数据的振幅减小，相位发生畸变，频带变窄。而图5d（Q补偿逆时偏移）中层位及盐丘下边界更清晰；通过与参考图像（5b）截取的子图及振幅谱对比，可见其损失的振幅与产生的相移都得到一定程度校正，与参考图像吻合度较高，频带得到拓宽。从图5a与图5b可知，采用传统的逆时偏移成像方法在成像剖面上存在严重的噪声及假象，降低成像精度，而采用上、下行波场分离成像的偏移结果则避免了这种影响，提高了分辨率。2.3抽稀的 Pluto模型为进一步验证 Q补偿逆时偏移对深层能量补偿的效果及与传统的拉普拉斯滤波方法对偏移成像噪声压制的效果，采用抽稀后的 Pluto模型进行数值模拟测试。其速度模型（图 8a）及对应的 Q模型（图 8b，由经验公式得出）的大小为 1160201，网格间距为dx=dz=10 m，时间采样间隔为 0.5 ms，采样点数为4002，震源采用主频为30 Hz的Ricker子波，总共114炮，炮间距为100 m，均匀分布在深度50 m处。图 9a为 1.0 s时刻未分离的波场快照，用上述方法分别得到该时刻经波场分离后的下行波场（图 9b）和上行波场（图9c）波场快照。图10a为未经波场分离的传统声波逆时偏移成像结果，分别得到声波逆时偏移经拉普拉斯滤波后成像结果（图 10c）、Q补偿逆时偏移经拉普拉斯滤波后成像结果（图10e）、波场分离后的声波逆时偏移成像结图 42.0 s时刻波场快照（a）原始波场；（b）下行波场；（c）上行波场图 5偏移成像结果（a）未波场分离的声波逆时偏移；（b）波场分离后的声波逆时偏移成像；（c）波场分离后的衰减数据偏移；（d）波场分离后的 Q补偿黏声逆时偏移878第 58 卷 第 4 期王喜鹏，等：基于上下行波分解的分数阶拉普拉斯算子黏滞声波逆时偏移成像方法图 8Pluto速度模型（a）及其 Q模型（b）果（图10b，参考）、波场分离后的衰减数据逆时偏移成像结果（图10d）及波场分离后的Q补偿逆时偏移成像结果（图10f）。图10b右侧子图为水平距离6000 m处作为参考的单道信息（黑色实线），图 10d、图 10f右侧子图则分别为水平距离 6000 m 处相应的单道信息（红色虚线）和参考单道信息（黑色实线）。图 11a 图11c对应为10b、图10d、图10f 中水平距离58 km、深度1.31.8 km范围的振幅谱。与图10b（参考）相比，图10d（衰减未补偿）因地层的衰减效应，其能量减弱，偏移成像层位不清晰，尤其是盐下边界以下；图 10f 由于采用 Q补偿逆时偏移，其能量得到补偿。对比 10b、图 10d、图 10f 的单道信息与图 11a图 11c 的振幅谱，可见缘于地层的非完全弹性效应，对成像结果造成振幅衰减及相位畸变，而采用Q补偿逆时偏移（图11c）能在一定程度上减缓这种效应的影响，拓宽其振幅谱，补偿衰减的能量，并对畸变的相位进行校正，提高了成像质量。由图10a可知，采用传统的互相关偏移方法，成像结果中浅层的有效信息完全湮没在噪声中，严重干扰成像精度；而采用拉普拉斯滤波后的成像结果（图10c、图10e）中低频噪声受到一定程度压制，但对浅层信息仍有较明显的干扰。综上所述，本文采用的Q补偿逆时偏移方法能有效图 6截取的 0.61.6 km层段单道信息（a）从图 5d截取；（b）从图 5c截取图 7不同逆时偏移的振幅谱（a）声波逆时偏移；（b）衰减数据逆时偏移；（c）Q补偿逆时偏移879石 油 地 球 物 理 勘 探2023 年补偿由地下介质非弹性特性造成的振幅损失和相移，采用该成像方法所得偏移成像结果（图 10d 图 10f）中，有效地压制由传统成像条件所引起的低频噪声及假象。图 91.0s时刻波场快照（a）原始波场；（b）下行波场；（c）上行波场图 10不同的逆时偏移成像结果 （a）未波场分离的声波逆时偏移；（b）波场分离后声波逆时偏移；（c）采用拉普拉斯滤波的声波逆时偏移；（d）波场分离后未补偿黏声逆时偏 移；（e）采用拉普拉斯滤波的 Q补偿逆时偏移；（f）波场分离后 Q补偿黏声逆时偏移880第 58 卷 第 4 期王喜鹏，等：基于上下行波分解的分数阶拉普拉斯算子黏滞声波逆时偏移成像方法图 11图 10d中红框对应位置振幅谱（a）声波逆时偏移；（b）衰减数据逆时偏移；（c）Q补偿逆时偏移值得注意的是，该方法是在进行两次震源波场和检波点波场外推的情况下，还要在时间切片上对波场进行傅里叶变换，所以计算时间稍长于采用拉普拉斯滤波的 QRTM，测试结果（表 1）显示约为其 1.3倍。但该方法不用对波场进行存储，避免了巨大的I/O量。表 1不同方案下的计算时间模型Pluto模型Sigsbee模型方案声波黏声波声波黏声波拉普拉斯滤波上下行分离拉普拉斯滤波上下行分离拉普拉斯滤波上下行分离拉普拉斯滤波上下行分离计算时间s502.05599.70572.41736.26415.11533.43769.371011.453结论本文提出一种上、下行波分解的解耦分数阶拉普拉斯算子黏滞声波逆时偏移成像方法：采用解耦的分数阶拉普拉斯算子黏滞声波方程对震源波场和检波点波场进行外推，能显著减缓由于地下介质的黏滞性所引起的振幅衰减和频散效应，提高了偏移成像的照明能量；采用分数阶拉普拉斯算子代替分数阶时间导数，不需在求解波动方程计算当前时刻波场值时依赖该时刻以前所有时段的波场。采用解析时间波场能显式地分解上下行波场。新的成像条件能避免同向传播的震源波场和检波点波场在应用互相关成像条件时所产生的低频噪声，消除由于偏移速度模型存在较强速度梯度时由震源上行波场和检波点下行波场互相关产生的假象，从而提高偏移成像的分辨率和质量。参 考 文 献1 陈生昌,刘亚楠,李代光.基于波动方程偏移的地震数据地层成像J.石油地球物理勘探,2022,57(5):11051113.CHEN Shengchang,LIU Yanan,LI Daiguang.Stratigraphic imaging with seismic data based on wave equation migrationJ.Oil Geophysical Prospecting,2022,57(5):11051113.2 吴丹,吴海莉,李群,等.应用自适应矩估计的快速最小二乘逆时偏移J.石油地球物理勘探,2022,57(2):386394.WU Dan,WU Haili,LI Qun,et al.Fast leastsquares reverse time migration with adaptive moment estimationJ.Oil Geophysical Prospecting,2022,57(2):386394.3 ZHOU J,BIRDUS S,HUNG B,et al.Compensating attenuation due to shallow gas through Q tomography and QPSDM,a case study in BrazilC.SEG Technical Program Expanded Abstracts,2011,30:33323336.4 KJARTANSSON E.Constant Qwave propagation and attenuationJ.Journal of Geophysical Research:Solid Earth,1979,84(B9):47374748.5 LIU H P,ANDERSON D L,KANAMORI H.Velo city dispersion due to anelasticity;implications for seismology and mantle compositionJ.Geophysical Journal International,1976,47(1):4158.6 EMMERICH H,KORN M.Incorporation of attenuation into timedomain computations of seismic wave fieldsJ.Geophysics,1987,52(9):12521264.7 PAPOULIA K D,PANOSKALTSIS V P,KURUP N V,et al.Rheological representation of fractional order viscoelastic material modelsJ.Rheologica Acta,2010,49(4):381400.8 MOCZO P,KRISTEK J.On the rheological models used for timedomain methods of seismic wave propagationJ.Geophysical Research Letters,2005,32(1):881石 油 地 球 物 理 勘 探2023 年L01306.9 BLANCH J O,ROBERTSSON J O A,SYMES W W.Modeling of a constant Q:Methodology and algorithm for an efficient and optimally inexpensive visco elastic techniqueJ.Geophysics,1995,60(1):176184.10 CHEN W,HOLM S.Fractional Laplacian timespace models for linear and nonlinear lossy media exhibiting arbitrary frequency powerlaw dependencyJ.The Journal of the Acoustical Society of America,2004,115(4):14241430.11 TREEBY B E,COX B T.Modeling power law absorption and dispersion for acoustic propagation using