李 娜

(中国石化中原油田分公司物探研究院，河南濮阳 457001)

·偏移成像·

基于Huber范数的多震源最小二乘逆时偏移

李 娜*

(中国石化中原油田分公司物探研究院，河南濮阳 457001)

李娜.基于Huber范数的多震源最小二乘逆时偏移.石油地球物理勘探，2017,52(5):941-947.

将Huber范数推广到最小二乘逆时偏移中，提升了算法的稳健性;同时考虑到最小二乘逆时偏移的计算量较大，将随机最优化相位编码技术引入到基于Huber范数的多震源最小二乘逆时偏移中，可有效提高计算效率,压制串扰噪声;最后通过常规有限差分正演模拟记录试算，验证了算法的有效性和对复杂模型的适用性。

最小二乘逆时偏移 相位编码 Huber范数 多震源 随机最优化

1 引言

近年来，随着油气勘探精度的逐渐提高，地震波成像方法也在不断发展[1]。相比于其他偏移方法，逆时偏移应用双程波波动方程进行波场延拓，避免了对波动方程的近似，无倾角限制，是目前较为精确的常规偏移成像方法。但逆时偏移仍属于常规偏移的范畴，其偏移算子是线性Born正演算子的共轭转置，而不是它的逆[2]。当地震数据采集不足或不规则、地下构造情况复杂以及波场带宽有限时，常规偏移方法只能对地下构造模糊成像，无法满足岩性油气藏勘探开发的需求。

最小二乘偏移将成像看作是最小二乘意义下的反演问题，通过迭代不断拟合观测数据和反偏移数据，可以得到振幅保真性更好、分辨率更高、偏移噪声更少的成像结果。与常规偏移类似，最小二乘偏移主要可以分为三类：射线类最小二乘偏移、单程波类最小二乘偏移和双程波类最小二乘偏移。射线类最小二乘偏移主要指最小二乘Kirchhoff偏移[3-5]和最小二乘高斯束偏移[6]。单程波类最小二乘偏移方法又可具体分为最小二乘分步傅里叶偏移[7,8]、最小二乘傅里叶有限差分偏移[9,10]和最小二乘广义屏偏移[11]等。双程波类最小二乘偏移即最小二乘逆时偏移(Least-Squares Reverse Time Migration,LSRTM)。Dai等[12]首先提出了线性和拟线性的LSRTM，并通过混合震源相位编码技术测试了方法的正确性和有效性；黄建平等[13]将LSRTM应用到近地表保幅成像；郭振波等[14]提出了LSRTM的真振幅成像；李振春等[15]将LSRTM推广到更加符合实际地下情况的黏声介质中。

上述LSRTM的目标泛函都是基于L2范数拟合，考虑到地震数据常含有较多噪声，常用的基于L2范数拟合的LSRTM算法对噪声非常敏感，尤其是当地震数据中含有奇异值时，常规LSRTM结果被噪声严重干扰。L1模对大噪声的容忍性比L2模更好，但由于L1模在0值处不可导，常用Huber模或者混合模替代[16]，目前Huber模已在全波形反演中取得较好效果。为此，本文将Huber模发展到LSRTM中，显著提升了LSRTM算法的稳定性，与常规LSRTM相比，本文算法在地震数据含有脉冲噪声时仍能得到较好的反演结果。

LSRTM的计算量过大，限制了其进一步推广应用。多震源相位编码技术可显著提高LSRTM的计算效率，Romero等[17]首次引入相位编码的概念；Krebs等[18]将相位编码技术应用于全波形反演中；Dai等[12]实现了基于相位编码的多震源LSRTM，有效降低了计算成本; 黄建平等[19]提出了基于平面波静态编码的LSRTM; 李庆洋等[20]提出伪深度域声波数值模拟方法；Moghaddam等[21]首先提出了随机最优化相位编码的全波形反演算法；随后，李庆洋等[22]将随机最优化思想推广到相位编码的LSRTM中，相比传统相位编码，LSRTM方法收敛更快。为此，本文将随机最优化相位编码技术应用于基于Huber泛函的多震源LSRTM中，与常规基于Huber范数的相位编码多震源LSRTM相比收敛更快，在适用于含脉冲噪声数据的同时显著提高了计算效率。

2 方法原理

2.1 线性Born正演

二维常密度声波方程可表示为

(1)

式中:s为慢度场；f为震源项；p是声压场。

(2)

(3)

将式(2)与式(3)相减，并应用Born近似，用背景波场p0替代总波场p0+ps，可得到扰动波场ps的控制方程

(4)

式(4)即为线性化Born正演(反偏移)方程，从中可以看出，扰动波场是背景波场与慢度扰动的相互作用，作为二次震源在背景介质中传播产生的场，与慢度扰动呈线性关系，具有明确的物理含义。即实际编程实现时，线性化正演(反偏移)过程需要求解两次正演方程：通过式(2)得到当前时刻的背景波场，然后再利用式(4)求得当前时刻的扰动波场。

利用背景介质中的格林函数G，可将式(4)中的扰动波场表示为

(5)

式中:Ω为积分域；m=-Δs2， 将其定义为模型参数。

为方便后续的推导，式(5)可写成算子的形式，即

ps=Lm

(6)

2.2 基于Huber模的LSRTM

常规LSRTM基于观测数据与反偏移模拟数据的L2模的拟合，其目标泛函为

(7)

式中：pobs为观测记录；ps=Lm为反偏移模拟记录； ‖·‖2代表向量的L2范数。

上述基于L2范数拟合的LSRTM算法对噪声非常敏感，尤其是当地震数据中含有奇异值时，LSRTM结果被噪声严重干扰。L1范数相比L2范数更加稳健，但当数据误差接近于0时，应用L1范数准则求取的目标函数梯度会不稳定，因而常用Huber模或混合模(hybrid)代替。两种目标泛函的形式如下

(8)

(9)

式中： Δd为模拟数据与观测数据的残差；ε=λ×mean(|pobs|)为控制L1范数与L2范数权重的阈值,λ为0～1之间的数，一般取0.5左右。

采用梯度导引类算法(最速下降或共轭梯度)求解，需要计算目标泛函关于模型参数的梯度，三种不同范数的梯度公式都可统一写成

(10)

rnorm=Wnorm(Lm-pobs)

(11)

其中，不同的范数对应不同的加权因子Wnorm，分别表示为

WL2=1

(12)

(13)

(14)

式中a为稳定因子，为避免分母为0，一般取0.01×mean(|pobs|)左右。

从式(10)可以看出，几种不同范数的梯度公式具有相同的形式，仅残差波场的加权因子不同。梯度求取过程也就是一个逆时偏移(RTM)的过程，与常规RTM不同的是，求梯度时的反传波场为残差记录。考虑到波场存储是RTM的瓶颈，本文计算梯度时采用震源波场重建方法，即每个时刻只在边界存储正传波场信息，反传时将其作为边界条件重构震源波场，大幅度降低了存储量[20]。

通过式(10)计算得到梯度后，模型的更新过程为

mk+1=mk-bkgk

(15)

式中：gk为第k次迭代的梯度；bk为第k次迭代的更新步长，步长可通过线性搜索等方法计算得到。

从梯度公式还可以看出，三种范数的唯一不同之处在于，L2模直接利用波场残差计算梯度，而Huber模和混合模则是利用加权后的残差记录计算梯度，因而三种范数的更新流程相同，常规L2模下的各种加速方法都可直接应用于Huber模和混合模中。

2.3 随机最优化相位编码技术

相位编码技术是目前非常实用的一项提高地震采集和数据处理效率的技术，其基本思想是通过设计编码函数将原本多次采集的炮集数据压缩成一次采集的混叠数据，从而减少数据规模，提高采集和处理效率。Moghaddam等[21]在多震源全波形反演中研究发现，通过加权平均前几次迭代的相位编码梯度可加快收敛速度。本文借鉴李庆洋等[22]提出的优化的多震源LSRTM方法，将随机最优化思想推广到基于Huber模或混合模的多震源LSRTM中。

编码函数多种多样，考虑到简便易实现的优点，本文选用震源极性编码，编码函数为

(16)

式中：Nnit为总迭代次数；s为当前炮号；k为当前迭代次数；γ为随机的+1或-1，在每次迭代中随机生成。

随机最优化思想通过加权平均前期得到的迭代梯度减小了梯度的随机波动，本文同样采用指数衰减加权平均，优化后的梯度表达式为

(17)

式中：j为加权的前期迭代次数，综合考虑效果和效率，一般取为10；c为衰减因子，取为0.5。相比于常规LSRTM，基于相位编码的LSRTM只需要修改线性正演算子和观测数据即可，具体的计算流程与常规LSRTM相同，在此不再赘述。

3 模型试算

3.1 简单层状模型

首先以一个简单的层状模型为例验证本文算法的正确性及有效性。模型速度场如图1左所示，模型中既有大尺度的背斜、向斜构造，也有90°高陡构造，非常适合检验成像方法的优劣，对应的扰动模型如图1右所示。网格数为401×271，纵、横向网格间距均为10m。合成观测数据所用计算参数为：在地表以50m间隔均匀分布81炮，每炮都是401个检波点全接收，震源为主频30Hz的雷克子波，时间采样间隔为1ms，最大记录时间为2.5s。

图1 速度模型(左)及反射率模型(右)

LSRTM的正算子为线性Born正演算子，而一般生成炮记录的正算子为常规有限差分正演算子，因而基于波形匹配策略的LSRTM需要首先测试两种算子是否匹配。图2左为常规有限差分正演模拟第41炮单炮记录，图2右为不含噪声的线性Born正演模拟单炮记录。可以看出，线性Born正演模拟记录与常规有限差分正演模拟记录相比，除了缺少明显的直达波外，基本相同，振幅相位都基本一致。为了更加清晰地对比二者的相似度，从图2左和图2右中分别取出5道单道波形进行对比，如图3所示。可以看出本文线性Born正演模拟算子与常规有限差分正演模拟算子匹配得较好，为其线性近似。

图2 有限差分单炮记录(左)及线性Born正演单炮记录(右)

将81炮利用震源极性编码方式组合成一个超道集，使计算量相当于单炮情形，缓解了计算需求。图4为图2左数据不同方法第50次迭代成像结果对比。可以看出，图4b的结果被噪声严重污染，而图4c结果消除了脉冲噪声的影响。对比图4c和图4d可以看出，本文方法可以有效改善成像效果，不仅能更好地压制串扰噪声，且收敛速度更快，成像结果分辨率更高。

图5左为四种算法归一化的数据残差曲线。可以看出，在含脉冲噪声的情况下常规L2模LSRTM不收敛，而基于Huber泛函的LSRTM仍然能较好地收敛。随机最优化编码方案相比常规相位编码的收敛速度更快，更加接近不含噪声情况下的收敛速度。由图5右常规相位编码和随机最优化编码策略的模型残差曲线同样可以看出随机最优化策略的优势。

需要注意的是，图4d虽然没有脉冲噪声的干扰，但仍然含有较多由相位编码引入的高频串扰噪声。为了压制串扰噪声，提高超道集的数目，将81炮利用震源极性编码方式组合成9个超道集，得到的第50次迭代成像结果如图6所示。对比图4d和图6可以看出，超道集数目的增加可以明显压制高频串扰噪声，得到同相轴更加清晰的成像结果。由于炮检距较大，成像结果中仍然含有一些低频噪声，这是由于线性Born正演模拟与常规有限差分模拟结果大炮检距不如小炮检距吻合得好，因而不能完全压制低频噪声，可以借助高通滤波、拉普拉斯滤波等进一步消除。

图3 图2抽取5道单道波形对比

图4 四种方法图2左数据切除直达波后第50次迭代成像结果对比

图5 数据残差收敛曲线(左)及模型残差收敛曲线(右)

图6 采用9个超道集的含噪声随机最优化Huber范数LSRTM结果

3.2 Marmousi模型

下面以Marmousi模型为例验证算法对复杂模型的适用性。Marmousi模型速度场如图7所示。模型参数如下：横向网格间距为12.5m，网格点数为481，纵向网格间距为8m，网格点数为375，速度范围为1.5～5.5km/s。采用时间二阶空间八阶有限差分正演模拟得到观测记录，观测系统模拟标准数据的观测系统，共有120炮，每一炮有101个接收道，炮间距为50m，道间距为25m，采用单边放炮，最小炮检距为50m，最大炮检距为2550m，震源为主频是20Hz的雷克子波，时间采样间隔为0.6ms，最大记录时间为3s。

图7 Marmousi速度模型

图8为四种方法逆时偏移成像结果对比。可以看出，不含噪声数据的常规RTM成像剖面中各个同相轴基本清晰可见，但振幅能量不均衡、分辨率较低。对比图8a和图8b可以看出，LSRTM的成像结果有明显的改善，不仅浅层的断层刻画更加清晰，而且深部的背斜也得到较好的凸显，整体能量更加均衡、分辨率明显提高(图8b)。图8c被噪声严重污染，而Huber泛函相比L2模具有更高的噪声容忍度，因而图8d结果显著消除了脉冲噪声的影响，验证了本文算法对复杂模型的适用性。

图8 四种方法逆时偏移结果对比

4 结论

(1)Huber模或者混合模相比L2范数更稳健，能较好地消除地震数据中的脉冲噪声，且计算流程与常规LSRTM相同，简便易实现。

(2)相位编码技术可显著提升LSRTM的计算效率，但收敛较慢。随机最优化的动态相位编码技术相比传统相位编码算法收敛更快，且可有效压制串扰噪声，进一步节省计算成本、提高计算效率，降低LSRTM的计算成本(同常规逆时偏移的计算成本相同)。

[1] 李振春.地震偏移成像技术研究现状与发展趋势.石油地球物理勘探,2014,49(1):1-21. Li Zhenchun.Research status and development trends for seismic migration technology.OGP,2014,49(1):1-21.

[2] Claerbout J F.Earth Soundings Analysis: Processing versus Inversion.Blackwell Scientific Publications,1992.

[3] Nemeth T,Wu C J,Schuster G T.Least-squares migration of incomplete reflection data.Geophysics,1999,64(1),208-221.

[4] 黄建平，李振春，孔雪等.碳酸盐岩裂缝型储层最小二乘偏移成像方法研究.地球物理学报,2013,56(5):1716-1725. Huang Jianping,Li Zhenchun,Kong Xue et al.A study of least-squares migration imaging method for fractured-type carbonate reservoir.Chinese Journal of Geo-physics,2013,56(5):1716-1725.

[5] 刘玉金,李振春,吴丹等.局部倾角约束最小二乘偏移方法研究.地球物理学报,2013,56(3):1003-1011. Liu Yujin,Li Zhenchun,Wu Dan et al.The research on local slope constrained least-squares migration.Chinese Journal of Geophysics,2013,56(3):1003-1011.

[6] Hu H,Liu Y,Zheng Y et al.Least-squares Gaussian beam migration.Geophysics,2016,81(3):S87-S100.

[7] Kuehl H,Sacchi M D.Least-squares wave-equationmigration for AVP/AVA inversion.Geophysics,2003,68(1):262-273.

[8] 黄建平,孙陨松,李振春等.一种基于分频编码的最小二乘裂步偏移方法.石油地球物理勘探,2014,49(4):702-707. Huang Jianping,Sun Yunsong,Li Zhenchun et al.Least-squares split-step migration based on frequencydivision encoding.OGP,2014,49(4):702-707.

[9] 黄建平，薛志广，步长城等.基于裂步DSR的最小二乘偏移方法.吉林大学学报(地球科学版),2014,44(1):369-374. Huang Jianping,Xue Zhiguang,Bu Changcheng et al.The study of least-squares migration method based on split-step DSR.Journal of Jilin University(Earth Science Edition),2014,44(1):369-374.

[10] 杨其强，张叔伦.最小二乘傅里叶有限差分法偏移.地球物理学进展,2008,23(2):433-437. Yang Qiqiang,Zhang Shulun.Least-squares Fourier finite-difference migration.Progress in Geophysics,2008,23(2):433-437.

[11] 周华敏，陈生昌，任浩然等.基于照明补偿的单程波最小二乘偏移.地球物理学报，2014,57(8)：2644-2655. Zhou Huamin,Chen Shengchang,Ren Haoran et al.One-way wave equation least-squares migration based on illumination compensation.Chinese Journal of Geophysics,2014,57(8):2644-2655.

[12] Dai W,Fowler P,Schuster G T.Multi-source least-squares reverse time migration.Geophysical Prospecting,2012,60(4):681-695.

[13] 黄建平，曹晓莉，李振春等.最小二乘逆时偏移在近地表高精度成像中的应用.石油地球物理勘探,2014,49(1):107-112. Huang Jianping,Cao Xiaoli,Li Zhenchun et al.Least square reverse time migration in high resolution imaging of near surface.OGP,2014,49(1):107-112.

[14] 郭振波，李振春.最小平方逆时偏移真振幅成像.石油地球物理勘探,2014,49(1):113-120. Guo Zhenbo,Li Zhenchun.True-amplitude imaging based on least-squares reverse time migration.OGP,2014,49(1):113-120.

[15] 李振春，郭振波，田坤.黏声介质最小平方逆时偏移.地球物理学报，2014，57(1):214-228. Li Zhenchun,Guo Zhenbo,Tian Kun.Least-squares reverse time migration in visco-acoustic medium.Chinese Journal of Geophysics,2014，57(1):214-228.

[16] Brossier R,Operto S and Virieux J.Which data resi-dual norm for robust elastic frequency-domain full waveform inversion.Geophysics,2010,75(3):R37-R46.

[17] Romero L A,Ghiglia D C,Ober C C et al.Phase encoding of shot records in prestack migration.Geophy-sics,2000,65(2):426-436.

[18] Krebs J R,Anderson J E,Hinkley.Fast full-wavefield seismic inversion using encoded sources.Geophysics,2009,74(6):WCC177-WCC188.

[19] 黄建平,李闯,李庆洋等.一种基于平面波静态编码的最小二乘逆时偏移方法.地球物理学报,2015,58(6):2046-2056. Huang Jianping,Li Chuang,Li Qingyang et al.Least-square reverse time migration with static plane-wave encoding.Chinese Journal of Geophysics,2015,58(6):2046-2056.

[20] 李庆洋,黄建平,李振春等.伪深度域声波数值模拟方法及应用.石油地球物理勘探,2015,50(2):283-289. Li Qingyang,Huang Jianping,Li Zhenchun et al.Acoustic wave numerical simulation in pseudo-depth domain.OGP,2015,50(2):283-289.

[21] Moghaddam P P,Keers H,Herrmann F J et al.A new optimization approach for source-encoding full-waveform inversion.Geophysics,2013,78(3):R125-R132.

[22] 李庆洋,黄建平,李振春等.优化的多震源最小二乘逆时偏移.石油地球物理勘探,2016,51(2):334-341. Li Qingyang,Huang Jianping,Li Zhenchun et al.Optimized multi-source least-squares reverse time migration.OGP,2016,51(2):334-341.

(本文编辑：金文昱)

李娜 在站博士后，1985年生；2007年毕业于山东理工大学数学统计学专业，获理学学士学位；2010年毕业于苏州大学概率论与数理统计专业，获理学硕士学位；2014年获中国石油大学(华东)地质资源与地质工程专业工学博士学位；现在中原油田物探研究院从事地震波传播正演模拟、偏移方法以及裂缝储层预测和描述方法研究。

1000-7210(2017)05-0941-07

P631

A

10.13810/j.cnki.issn.1000-7210.2017.05.006

*河南省濮阳市华龙区庆山街3号中原油田物探研究院，457001。Email：lina19202@163.com

本文于2016年8月29日收到，最终修改稿于2017年6月23日收到。

本项研究受国家油气重大专项项目(2016ZX05006-004)资助。