带限局部平面波传播算子与偏移方法

2017-11-01刘少勇韩冰凯顾汉明宋桂桥

石油地球物理勘探 2017年5期

关键词：平面波波场射线

刘少勇韩冰凯* 顾汉明宋桂桥

(①中国地质大学地球物理与空间信息学院，湖北武汉 430074; ②地球内部多尺度成像湖北省重点实验室，湖北武汉 430074； ③中国石化油田勘探开发事业部，北京 100728)

带限局部平面波传播算子与偏移方法

刘少勇①②韩冰凯*①②顾汉明①②宋桂桥③

刘少勇，韩冰凯，顾汉明，宋桂桥.带限局部平面波传播算子与偏移方法.石油地球物理勘探，2017,52(5):948-955.

地震信号的带限特点限制了高频近似下的射线类算子在地震波模拟和偏移中的应用。本文基于射线理论构建带限信号局部平面波传播算子，兼顾射线类方法具有时空局部性和实现灵活高效的特点，有效进行带限地震波场的模拟和偏移。即分析传统高频射线在速度界面的传播特征，利用带限信号频带与菲涅耳带的关系，构建带限斯奈尔定律进行中心射线传播；在局部平面波假设下，基于带限中心射线做旁轴近似展开，构建带限局部平面波传播算子；在此基础上发展基于该算子的带限射线束偏移方法。数值实验结果表明，带限射线束较传统基于高频近似的射线束具有更强穿透能力，改善射线阴影区照明，从而提高复杂构造区成像质量。

射线束偏移带限信号传播算子局部平面波

1 引言

基于高频近似的射线类传播算子被广泛应用于描述地震波传播和偏移成像，但实际采集的地震数据是具有一定带宽的数字信号。在高频渐进射线理论下，地震波旅行时通常被表示为慢度沿着射线路径的积分[1]。Woodward[2]提出，在一定频带范围内的地震波旅行时通常是由慢度在炮检点之间射线路径周围一定范围内的空间体积分决定的，该旅行时与地震波频带相关。

Biondi[3]通过对赫姆霍兹方程进行近似，推导出频率相关的程函方程进行旅行时计算，该方法在计算低频部分旅行时时受限于高频端的频带范围。Hogan等[4]讨论了求解频率相关程函方程计算低频波场的条件，即介质的平滑程度。另一类频率相关的旅行时计算方法可通过射线追踪实现，Foreman[5]从赫姆霍兹方程出发推导了其对应的频率相关的射线追踪系统。这种方法相比于传统的射线追踪，计算复杂度显著增加。Protasov等[6]分析该频率相关的射线追踪算法的特点，将其与射线追踪、时间域有限差分正演的波场进行对比分析，指出此方法得到带限射线的射线路径与震源主频有关，能较好地改善射线阴影区问题。Lomax[7]利用波长平滑进行带限波场传播，具体通过对垂直于射线平面的速度加权平均来实现，其加权函数的宽度正比于波长。Lomax等[8]分析了波长平滑算法的原理，指出通过特定的加权函数可控制带限波场的运动学特征，但这种加权平均的方法在速度差异较大的不规则界面处难以准确地刻画射线的传播。Protasov等[9]研究了射线在复杂速度界面传播特征，其核心思想是在复杂界面利用Kirchhoff边界积分对积分区域内射线进行先加权后传播。Protasov等[10]指出利用该方法构造出的带限射线方向对应的是等效波场透过界面的最大能量方向，并提出通过构建界面处菲涅耳带透射波场的最大值问题来求解带限射线路径。Yarman等[11]对Kirchhoff边界积分引入射线级数解，通过积分表达式中相位部分的近似导出积分区间为第一菲涅耳带，并构建该积分区间和波场频带的关系，将此方法称为带限射线追踪。

本文结合带限射线和射线束的特点，从分析传统高频射线与带限射线的传播特征出发，研究适应于复杂介质的带限射线追踪。在此基础上，通过旁轴近似构建带限局部平面波传播算子模拟带限波场，改善传统射线类方法在射线阴影区和中深层的照明及偏移效果。

2 方法原理

2.1 带限局部平面波传播算子

高频射线在速度界面的传播可由斯奈尔定律控制，带限射线则对应于带限波场传播的最大能量方向[9]。带限斯奈尔定律可控制带限射线的传播，但没有传统斯奈尔定律的简洁数学表达式。带限射线通过界面时的影响范围不再是一个点，而是由菲涅耳带控制的一个范围，在数学上可表达为求解透射波场能量最大的优化问题对应的信任区间。在速度界面处，中心射线透过界面上一点x0的透射波场可表示为

(1)

式中：x为地下成像点；A(x)为与格林函数有关的振幅项；T(x)为透射系数；积分区间FZx0为中心射线路径坐标x0附近的菲涅耳带的范围，由频率f和与旅行时有关的相位项Δτ(x,x0)决定。如图1所示，中心射线穿过速度界面时，在其附近的菲涅耳带FZx0内构建局部平面Γx0。局部平面内，以高频射线描述等效局部平面波的传播，入射射线参数pin(x)和出射射线参数pout(x)遵循运动学射线方程组

(2)

式中：p表示射线的慢度矢量；v表示射线路径上的速度；s和τ分别表示弧长和射线的旅行时；x∈Γx0。由高频射线的出射射线参数pout(x)加权得到等效局部平面波的出射方向pe(x0)，对应透射波场u(x0)的最大能量方向，即带限中心射线在速度界面处以带限斯奈尔定律传播

|u(x0)|

(3)

基于旁轴近似理论，由带限中心射线可构建带限局部平面波传播算子。局部平面波宽度范围的旅行时由中心射线旁轴近似展开[24]

图1 带限射线在速度界面的传播示意图

(4)

式中： Δs为x在中心射线上投影点与x0间距离；r为x到中心射线的距离(图2a)；mr表征波前曲率

(5)

(6)

式中：v0表示中心射线初射点的速度； Δθ0表示中心射线初射角度间隔(图2b)。

图2 笛卡尔坐标系下中心射线旅行时旁轴近似展开(a)与射线束宽度(b)

2.2 带限射线束偏移

射线束偏移区别于Kirchhoff积分偏移单道输入的特点，其对应的输入数据也是由地表地震记录分解的局部平面波。传统τ-p变换可将地表地震记录分解成不同方向的局部平面波。这种线性Randon变换方法由于其基函数不完备，往往存在能量泄露而降低局部平面波的分辨率。平面波分解通常可统一在反演理论框架下进行，通过构建反问题求解以获得更高的分辨率。本文通过在反问题中引入L0范数的约束，在压缩感知框架下求解局部平面波问题。基于分解后的局部平面波进行带限射线束的成像。

在频率域中，构建平面波分解的正问题描述为：局部平面波源经过传播矩阵得到地表局部射线束范围内地震数据

d(x,ω)=A(x,pe,ω)D(pe,ω)

(7)

式中：d(x,ω)表示地表局部地震数据；D(pe,ω)表示待求解的局部平面波；A(x,pe,ω)表示局部平面波源和地表地震数据之间的投影矩阵。通过求解对应反问题获得射线束控制范围内的局部平面波源，即给定相移矩阵和空间局部信号求解局部平面波。该反问题的实现通常借助于求解如下最小二乘泛函实现

J=‖A(x,pe,ω)D(pe,ω)-d(x,ω)‖2

(8)

矩阵A(x,pe,ω)通常为欠定矩阵，造成上述反演问题不稳定，可通过加入模型约束求最小二乘解，但是其结果通常受到能量泄露的影响。叠前地震数据可以通过在变换域中实现数据的稀疏表达，求解欠定问题比较有效的方法是假设模型具有一定的稀疏性[28]，对应在局部平面波分解中，该稀疏性的物理解释为：在局部范围内仅存在少数几个平面波源。这种假设在局部空间窗数据中通常成立，可借助于对模型的L1和L0范数约束来实现[29]。更进一步，为了得到更稀疏的局部平面波解，可将模型估计改为L0范数约束

min ‖D(pe,ω)‖0使得

‖A(x,pe,ω)D(pe,ω)-d(x,ω)‖2<ε

(9)

该模型的L0范数约束代表解的个数极小[28]，在射线束合成中代表局部平面波的个数较少。模型L0范数约束是模型稀疏约束比较极端的情况，通常该泛函借助一些贪婪算法实现，比如匹配追踪等方法。

基于以上局部平面波数据，利用带限局部平面波传播算子进行炮检端的波传播加上成像条件就可以进行成像处理。射线束成像条件类似于传统Kirchhoff积分偏移的等时成像条件，只是在成像过程前已经对数据进行了方向筛选。传统的Kirchhoff积分偏移成像公式可写成[30]

eiωτdωdΩ

(10)

式中：Ω表示包含炮点xs和中心检波点xr0的成像孔径；ps和pr分别表示炮点、检波点射线参数；W为加权函数；D表示频率域局部平面波；τ为旅行时。在声学介质中，旅行时表达是常数，D在频率域的积分可以通过与δ函数卷积得到时间域表达的局部平面波d。由此，通过引入激励时间成像条件去除式(10)的频率域积分，得到时空域成像公式

pr(xr0),τ]dΩ|τ=τs+τr

(11)

3 带限射线束偏移的实现

地震数据偏移的实现通常有两个核心步骤：炮检波场的传播和成像条件的时间。射线束偏移一般采用局部平面波传播算子，本节从方法的具体实现出发，对带限局部平面波传播算子的构建和带限射线束的偏移流程进行说明。

如图3所示，时空域共炮集记录的带限射线束偏移流程中带限局部平面波传播算子的构建分为三个步骤： ①构建局部平面波，当中心射线经过模型空间的非均匀区域时，在射线束宽度范围内的局部平面内以一定密度高频射线近似局部平面波的传播； ②以带限射线追踪计算射线路径和旅行时，在局部平面内加权计算等效射线参数作为中心射线的参数，以带限斯奈尔定律描述中心射线的传播； ③带限中心射线旅行时和振幅信息由旁轴近似展开分配到射线束宽度范围内的成像网格点上。

图3 带限射线束偏移实现流程图

结合当前计算集群多节点、多核心共享存储的特点，为了提高带限射线束偏移的效率，其实现过程需遵循如下原则： ①炮点端各方向格林函数依次计算并保存到内存中，避免后续炮检点射线束相交成像过程中对炮点正问题的重复计算； ②对炮循环利用消息传递接口(MPI)进行并行，内部空间方向成像利用OpenMP线程并行，利用计算机单节点多核的特点节省内存； ③一次输入单炮数据，成像结果记录于当前计算节点的本地盘中，最后通过并行规约(Reduce)统一抽取成像道集，提高节点本地存储的I/O利用率。

4 数值算例

首先对水平海底模型进行带限射线传播的测试，其中海底以下为高速(4000m/s)介质，用水平海底模拟速度跃变界面。在该模型坐标(1200m，400m)处设置震源，从震源出发进行带限射线束传播(射线出射角度为25°)，并以有限差分正演模拟的波场快照作为参考(主频为20Hz，时刻为1.0s)。如图4所示，传统的高频射线无法穿过速度跃变海底，带限射线却能穿过速度跃变的海底，带限射线宽度由绿色实线表示，带限射线的1.0s旅行时对应坐标位置和局部波前分别由红色小圆圈和红色实线表示。带限射线穿透能力较高频射线强，且随着带限射线频率的降低，其穿透能力增强。对比波场快照，带限射线旅行时与有限差分计算的波前能准确地吻合。

图4 水平海底模型高频射线与带限射线路径比较(射线初射角度为25°)

(a)高频射线；(b)、(c)、(d)依次代表高、中、低频带限射线。红色实线表示中心射线，绿色实线标注射线束宽度，红色小圆圈及其对应红色实线标注旅行时1．0s的局部波前位置，对应20Hz有限差分正演模拟1．0s波场快照

进一步，为验证带限射线束进行波传播的有效性和准确性，同时验证带限射线束成像的效果，设计一个崎岖海底模型进行测试。在图5a中，以正弦函数型界面模拟速度跃变的崎岖海底，对于高频射线，崎岖海底下方是典型的射线阴影区。如图5b和图5c所示，在模型坐标(4000m，0)处设置震源，高频射线和带限射线的射线路径叠加在有限差分正演模拟波场快照(主频为20Hz，时刻为1.2s)的背景上。对比崎岖海底模型中的射线分布可以看出，带限射线较传统高频射线具有更强的穿透能力，改善了崎岖海底下方的照明，带限射线也能更好地与波场快照吻合。图5d和图5e分别为传统射线束和带限射线束的成像结果，相比于基于高频射线的射线束偏移结果，带限射线束偏移结果提高了目标层的连续性，其中红色箭头标注的水平目标层的成像能量强且均匀。在偏移效率方面，带限射线束和传统射线束由于在相同的实现框架下，因此理论上有类似的计算效率。其执行效率的差别源于正算子中计算格林函数效率不同，单炮测试中，带限射线束传播算子相比传统射线束传播算子增加计算量达373.7%。但由于该模块并不是偏移过程中的热点，因此其导致的偏移效率计算量只增加了不到18%，详细数据如表1所示。

最后，选择更一般的Salt Bag模型[31]进行盐下成像的测试，图6a所示模型包含不规则盐丘体，盐丘上界面为一崎岖界面。如图6b和图6c所示，在模型坐标(5000m，0)处设置震源，高频射线和带限射线的射线路径叠加在有限差分正演模拟波场快照(主频为20Hz，时刻为1.2s)的背景上。对比高频射线与带限射线的射线路径可知，带限射线能更好地穿透盐丘，对盐下区域有更好的照明，这是利用改善射线类方法对盐下成像的基础。图6d和图6e分别为对应的传统射线束和带限射线束的偏移成像结果，对比红色矩形框内的成像区域，带限射线束偏移对盐丘下界面和盐丘下方的目标层位成像效果优于传统高频射线束偏移。

图5 崎岖海底模型高频射线与带限射线传播对比以及对应射线束成像结果

(a)崎岖海底模型，海底下方为高速(4500m/s)盐丘体； (b)高频射线路径叠加于20Hz有限差分正演模拟1.2s波场快照，红色实线表示射线路径，绿色实心圆点射线路径对应1.2s旅行时的波前，图c同； (c)带限射线路径叠加于20Hz有限差分正演模拟1.2s波场快照； (d)基于高频射线的射线束成像剖面； (e)带限射线束成像剖面

图6 Salt Bag模型高频射线与带限射线传播对比以及对应射线束成像结果

(a)Salt Bag速度模型； (b)高频射线路径叠加于20Hz有限差分正演模拟1.2s波场快照，红色实线表示射线路径，绿色实心圆点射线路径对应1.2s旅行时的波前，图c同； (c)带限射线路径叠加于20Hz有限差分正演模拟1.2s波场快照； (d)基于高频射线的射线束成像剖面； (e)带限射线束偏移成像剖面

表1 崎岖海底模型传统射线束与带限射线束偏移效率对比

5 结论与讨论

本文在射线理论框架下，基于带限中心射线，利用旁轴近似构建带限局部平面波传播算子，进行带限地震波场的传播和偏移。该算子在理论上保留了传统射线束传播算子灵活高效和具有局部时空性的优点，同时兼顾波动类算子能够刻画带限信号的特点。在实际应用中带限局部平面波传播能提高射线束的穿透能力，改善传统高频射线理论对应的射线阴影区等问题，有利于改善复杂构造的照明。基于该带限传播算子的成像方法提高了崎岖海底和盐丘等复杂构造的偏移成像质量。数值实验结果也表明本文提出的带限局部平面波传播算子能够准确有效地描述带限波场的传播，改善射线阴影区和盐下构造等复杂探区的成像质量。本文研究为基于带限射线束的层析和反演等后续研究提供了基础。

[1] Babich V M and Buldyrev V S.Asymptotic Methods in Problems of Diffraction of Short Waves.Nauka，Moscow,1972.

[2] Woodward M J.Wave-equation tomography.Geo-physics，1992，57(1)：15-26.

[3] Biondi B.Solving the frequency dependent eikonal equation.SEG Technical Program Expanded Abstracts，1992，11：1315-1319.

[4] Hogan C M，Margrave G F.Ray-tracing and eikonal solutions for low-frequency wavefields.CREWES Research Report，2007.

[5] Foreman T.A Frequency Dependent Ray Theory[D].The University of Texas at Austin，1987.

[6] Protasov M，Gadylshin K.Exact frequency dependent rays on the basis of Helmholtz solver.SEG Technical Program Expanded Abstracts，2015，34：3739-3743.

[7] Lomax A.The wavelength-smoothing method for approximating broad-band wave propagation through complicated velocity structures.Geophysical Journal International，1994，117(2)：313-334.

[8] Lomax A，Snieder R.Estimation of finite-frequency waveforms through wavelength-dependent averaging of velocity.Geophysical Journal International，1996，126(2)：369-381.

[9] Protasov M I，Yarman C E，Nichols D et al.Frequency-dependent ray-tracing through rugose interfaces.SEG Technical Program Expanded Abstracts，2011，30：2992-2996.

[10] Protasov M I，Osypov K S.Frequency dependent ray tracing for irregular boundaries.Seismic Technology,2014,11(3):1-11.

[11] Yarman C E，Cheng X，Osypov K et al.Band-limited ray tracing.Geophysical Prospecting，2013，61(6)：1194-1205.

[13] Hill N R.Gaussian beam migration.Geophysics，1990，55(11)：1416-1428.

[14] Hill N R.Prestack Gaussian-beam depth migration.Geophysics，2001，66(4)：1240-1250.

[15] Hale D.Migration by the Kirchhoff，slant stack，and Gaussian beam methods.CWP Annual Project Review Meeting，Colorado，1992.

[16] Gray S H.Gaussian beam migration of common-shot records.Geophysics，2005，70(1)：953-959.

[17] 邓飞，刘超颖，赵波等.高斯射线束正演与偏移.石油地球物理勘探，2009，44(3)：265-269. Deng Fei，Liu Chaoying，Zhao Bo et al.Gaussian beams forward simulation and migration.OGP，2009，44(3)：265-269.

[18] 李振春，岳玉波，郭朝斌等.高斯波束共角度保幅深度偏移.石油地球物理勘探，2010，43(3)：360-365. Li Zhenchun，Yue Yubo，Guo Chaobin et al.Gaussian beam common angle preserved-amplitude migration.OGP，2010，45(3)：360-365.

[19] Zhu T，Gray S H，Wang D.Prestack Gaussian-beam depth migration in anisotropic media.Geophysics，2007，72(3)：S133-S138.

[20] 韩建光，王赟，张晓波等.VTI介质高斯束叠前深度偏移.石油地球物理勘探，2015，50(2)：267-273. Han Jianguang，Wang Yun，Zhang Xiaobo et al.Gaussian beam prestack depth migration in VTI media.OGP，2015，50(2)：267-273.

[21] 刘强，张敏，李振春等.各向异性介质共炮域高斯束偏移.石油地球物理勘探，2016，51(5)：930-937. Liu Qiang，Zhang Min，Li Zhenchun et al.Common-shot domain Gaussian beam migration in anisotropic media.OGP，51(5)：930-937.

[22] 代福材，黄建平，李振春等.角度域黏声介质高斯束叠前深度偏移方法.石油地球物理勘探，2017，52(2)：283-293. Dai Fucai，Huang Jianping，Li Zhenchun et al.Angle domain prestack Gaussian beam migration for visco-acoustic media.OGP，2017，52(2)：283-293.

[23] Hu C S，Stoffa P L.Slowness-driven Gaussian-beam prestack depth migration for low-fold seismic data.Geophysics，2009，74(6)：WCA35-WCA45.

[24] Liu J，Palacharla G.Multiarrival Kirchhoff beam migration.Geophysics，2011，76(5)：WB109-WB118.

[25] Bube K P，Washbourne J K.Wave tracing：Ray tracing for the propagation of band-limited signals：Part 1 —Theory.Geophysics，2008，73(5)：VE377-VE384.

[26] Washbourne J K，Bube K P，Carillo P et al.Wave tra-cing：Ray tracing for the propagation of band-limited signals：Part 2 — Applications.Geophysics，2008，73(5)：VE385-VE393.

[27] Gray S H，Xie Y，Notfors C et al.Taking apart beam migration.The Leading Edge，2009，28(9)：1098-1108.

[28] Elad M.Sparse and Redundant Representations：from Theory to Applications in Signal and Image Processing.Springer Science & Business Media，2010.

[29] Liu S Y，Gu H and Feng B.Characteristic wave imaging in TTI media via compressed sensing.SEG Technical Program Expanded Abstracts，2016，35：4424-4429.

[30] Schneider W A.Integral formulation for migration in two and three dimensions.Geophysics，1978，43(1)：49-76.

[31] Popov M M，Semtchenok N M，Popov P M et al.Depth migration by the Gaussian beam summation method.Geophysics，2010，75(2)：S81-S93.

(本文编辑：朱汉东)

1000-7210(2017)05-0948-08

P631

10.13810/j.cnki.issn.1000-7210.2017.05.007