《1. 引言》

1. 引言

逆时偏移成像（reverse-time migration， RTM）计算成本虽高，但近20年来随着硬件的发展，被广泛应用于近地表地震和探地雷达（ground-penetrating radar， GPR）数据处理中。常用来处理大型地震勘探数据，并成为首选成像方法。利用波的全方位传播特性，避免了其他偏移方法中的近似求解，实现了对陡倾角界面和横向速度变化介质的精准成像。另外，该方法还改善了波的衰减和振幅损失问题[1]。本文主要研究叠后RTM。

近几年来，大量学者将逆时偏移算法引入GPR数据处理中。研究初期，Fisher等[2]利用声波方程和爆炸反射界面模型来进行GPR数据的叠后偏移成像。Sanada和Ashida[3]研究了电导率对基于麦克斯韦方程组的叠后RTM产生的影响。

由于GPR图像的构建严重受制于地形起伏和勘探深度变化，针对基于高程静校正的叠后偏移等常规方法在处理时可能出现较大误差的情况，发展新的处理方法很有必要（图1）[4,5]。Heincke等[6]提出了另一种地形偏移方法，该方法在形态图像处理后，利用衍射相似性来识别复杂地形中的陡倾角反射界面。Zhu和Huang[7]摆脱声波方程的限制，应用经典的Yee氏算法[8]来解决RTM中爆炸反射界面产生的零偏移距问题，并通过沙丘地层模型正演数据检验了算法。研究表明，RTM虽然可以提高沙丘中常见的陡地层界面成像精度，但仍需对野外数据做进一步研究。

《图1》

图1.（ a）沿校正基准面采集到的零偏移距旅行路径剖面图（垂直虚线表示标准高程静校正中采集面到基准面的校正；实线表示校正后的射线传播路径），这种方法导致基于高程静校正基准面偏移结果较差（改编自参考文献[4]）；（b）为沿（a）所示的采集面得到的合成GPR数据（εr是相对介电常数）。

针对高程静校正偏移所存在的缺陷，本文提出了一种基于爆炸反射界面模型的叠后RTM算法，采用直接由地表向下的波场延拓方法，避免了对高程静校正或基准面的计算，从而解决了崎岖地形成像问题。该算法是基于麦克斯韦方程的波动方程解，可用于解释横向复杂速度变化和电导率分布情况，并通过数值模拟和复杂环境下的野外沙丘数据对该算法进行了检验。

《2. RTM原理》

2. RTM原理

RTM是利用负时间步长求解麦克斯韦方程组，并在表面以逆时序列插入雷达图像作为边界条件，通过将记录波场向下延拓来实现的。当时间逆推至零时刻时，即能量全部收敛回模型且向下回归至初始点，满足叠后偏移成像条件。该算法中包含两个熟知的注意事项，首先，必须对速度场进行平滑处理，以防延拓时出现虚假反射，在子波主频下，使用两个波长的三角平滑方法，可以最大限度地降低虚假反射；其次，爆炸反射界面模型只对一次反射波有较好的处理效果，如果数据中存在高幅多次波，则会使成像产生误差。

先求解麦克斯韦方程的解耦和二阶微分形式。在二维（two-dimensional，2D）平面中，电场分量垂直于成像平面[横向电场（transverse electric，TE）模型]，这些方程被简化为阻尼标量波动方程：

式中，E为电场；μ为磁导率；t为时间；x和z为空间变量；ε和σ为介电常数和电导率。本算法中，假设磁导率μ等于自由空间的磁导率μ0，利用损耗角正切σ/(ωε)<<1进行低损耗近似。在这种情况下，速度υ和衰减系数α如下：

在爆炸反射界面模型中，波场近似单向传播，模型速度是真实速度的一半υ′=υ/2，式中，激发参数表示爆炸反射模型中使用的有效值。根据公式（2）求解公式（1）时，需进行变量替换，令ε′=4ε。

公式（1）时，需进行变量替换，令ε′=4ε。为了确保初至波在单向传播中正常衰减，公式（3）需满足：

式中，z为反射深度。公式（4）右侧是双程传播时的衰减，左侧是爆炸反射界面模型单程传播时间的衰减。左右两边消除公因子、移项并将ε′=4ε代入整理后，结果如下：

将本构参数εꞌ和σꞌ代入公式（1），得到修正后的爆炸反射界面模型波动方程。

利用二阶精确有限差分算法求解波动方程，并使用Zhou等[9]提出的完全匹配层（perfectly matched layer，PML）作为吸收边界条件，来计算边界区域中的耦合各向异性方程组。例如，z方向上，用Zhou等的公式，在边界区域求出的微分方程组为

式中，辅助变量D1z和D2z由公式（6）定义。变量σz和az是PML吸收参数。

本文通过MATLAB编程实现该算法，并利用矩阵计算和多核/多处理器系统来提高计算效率。该方法沿着记录面或地表，逆时输入记录波场作为边界条件，使用负时间步长对公式（1）进行波场延拓。由于算法使用正方形网格进行剖分，因此沿着起伏地表的网格不连续性会引起少量误差。但这些不连续单元尺寸远小于波长，所以引入的误差较小（数据中最高频率处产生的误差小于λ/10）。例如，对于频率为100MHz的Ricker子波，数据中包含的最高频率约为2.5f0或250MHz。与典型的最大相对介电常数εr=9、网格间距4cm相比，主频为100MHz时波长为1m。计算网格中的散射缺失是判断近似合理很好的指标。

《3. 数值模拟》

3. 数值模拟

《3.1.正演模型1——简单速度模型》

3.1.正演模型1——简单速度模型

首先，使用简单的层状速度模型来研究地形起伏较大时产生的影响[图1（a）]。相对研究深度而言，该模型横向地形变化较大，这是GPR调查中经常遇到的情况。模型地表为正弦形态，从底部到顶部的高度为4m，不规则反射层深度变化范围为1~5m。在模拟的GPR数据中，存在单个反射界面，产生了复杂的波场[图1（b）]，形成了由多个向斜叠加而成的三重反射特征。使用三种常规处理方法：①基于近似基准面的逆时偏移高程静校正[图2（a）]；②基于高程静校正的基准面逆时偏移[图2（b）]；③地形逆时偏移[图2（c）]。图2（a）中，在没有地形校正的情况下进行偏移，数据会产生较大的误差，图像聚焦不良，反射不能正确归位。在7~10m，数据出现了过偏移现象；而13~18m，背斜几乎是平坦的。图2（b）进行标准高程静校正后，偏移结果得到了改善，但图像重建仍然包含大量误差。沿剖面看，向斜集中在15m处，且出现了过偏移现象，这可能会得出速度模型错误的结论。在6~8m该区地形最低点和陡峭反射短面上，由于图像不能正确聚焦，从而产生了两个明显的界面。相较其他方法而言，图2（c）中的地形逆时偏移法能很好地聚焦图像，并将反射界面归位到正确的深度位置上。

《图2》

图2. 图1（b）中正演数据的三种不同的处理方法。（a）基于近似基准面的逆时偏移高程静校正；（b）基于高程静校正的基准面逆时偏移；（c）地形逆时偏移[（a）和（b）中的图像聚焦不良，而经过地形逆时偏移处理后的（c）图，能正确地聚焦图像并将反射界面归位到正确的深度处，如（c）中黑线所示]。

《3.2.正演模型2——横向速度突变模型》

3.2.正演模型2——横向速度突变模型

逆时偏移的主要优势之一是能够处理横向速度变化较大的资料。为了检验该算法对横向速度突变且地形起伏较大模型的处理效果，本文建立了如图3所示的模型。该模型的结构与模型1相同，但第一层从右往左介电常数陡增（速度降低）。正演结果与图1（b）相似，但模型左边的旅行时存在明显延迟，说明该处的速度比模型1要低。

对上述模型使用两种方法进行处理：①传统的基于高程静校正的逆时偏移；②地形逆时偏移。图3（c）中，在高程静校正时高速层速度会替换原始地层速度，所以通过基准面逆时偏移后会得到一个修正的速度模型，导致图像无法聚焦，反射不能正确归位，与图2（b）类似。而对这种表层横向速度突变模型，地形逆时偏移能精准成像，并得到清晰的反射边界[图3（d）]。值得注意的是，较陡的反射层（主要在7m附近）很难成像，因为绝大部分能量不会返回到地表，致使偏移出来的地层振幅较弱。另外，由于模型左边的反射系数较小，使得10m之前的反射界面振幅都很弱。

《图3》

图3.（ a）第一层10 m处速度突变模型（；b）由（a）模型得到的正演图形；（c）基于高程静校正的基准面逆时偏移图像；（d）地形RTM图。

《4. 现场试验——美国犹他州珊瑚粉沙丘》

4. 现场试验——美国犹他州珊瑚粉沙丘

珊瑚粉沙丘（Coral Pink Sand Dunes, CPSD）位于犹他州南部，是大盆地与科罗拉多高原之间的过渡地带中最大的沙丘之一。沙丘位于纳瓦霍基岩之上，并被塞维尔正断层分割成两块，因此沿沙丘下部（lower dune field， LDF）的东部边界形成了基岩悬崖。为了检验“沿下伏基岩表面受断层控制的地形决定了沙丘的构造和形态”这一假设，我们开展了一项基于GPR的研究，主要目的是绘制沙丘和基岩面交界处的地形图并对基岩内部的构造特征进行成像。

现场测试时，使用50MHz和100MHz的天线对25个样区内超过20km的剖面进行数据采集，通过连续全球定位系统（global positioning system， GPS）进行高程控制，并在之后的处理中进行差分校正。GPS基站距离所有样区均不超过10km，且通过将GPR采集时钟与GPS时钟同步来链接GPS和GPR位置。地质雷达信号的穿透力很强，在有些地区可以记录到深度超过35m的反射，因此能够对地表现代沙丘和地下古代沙丘地层清晰成像。样区内的一些露头和浅层钻孔提供了沙丘和基岩接触的真实情况。虽然GPR信号质量非常好，但样区地形和地层复杂性限制了成像精度。在样区某些剖面上，地表起伏达到了25m，坡度超过了30°，这样的地形结构能很好地检验逆时偏移算法。

对采集数据使用两种方法进行处理：①基于高程静校正的逆时偏移；②地形逆时偏移。两种方法的处理步骤如表1所示。值得注意的是，通过对几种不同介电常数模型偏移后的衍射双曲线衰减情况进行检查后，最终确定了一个恒定的介电常数，并应用在整个区域得到了较好的偏移效果。由于地形和地层是非常重要的影响因素，可利用该区域理想的数据从地形角度来评估逆时偏移效果。

《表1》

表1 高程静校正逆时偏移和地形逆时偏移处理流程对比

《4.1.场地1》

4.1.场地1

图4和图5为一组典型特征剖面，一般而言，我们看到的现代沙丘都是覆盖在古代沙丘岩化形成的基岩面上，而基岩往往暴露在沙丘内部区域。区域1是基岩地层中的一块陡峭河床，在基准面逆时偏移后，河床反射界面没有聚焦并被错误归位，如左边的斜坡沉积在30~35m深度穿过了一个更古老更深的界面。而在地形逆时偏移后，这块陡峭的地层被较好地聚焦和归位。区域2的基岩中含有一个正断层。基准面逆时偏移没有很好地聚焦断层面，很难解释地层和断层的相对几何关系。在地形逆时偏移中，断层面聚焦很好并且能够看到复杂的破碎结构。值得注意的是，这条断层没有切断近地表的侵蚀接触面，因此在近代是没有发育的断层。

图4. CPSD示例成像剖面。（a）偏移前的预处理数据；（b）地形逆时偏移数据。

《图5》

图5.（ a）区域1的基准面逆时偏移图像；（b）区域1的地形逆时偏移图像[通过比较（a）和（b），发现在基岩横向50 m内深度为20 ～ 30 m的一组陡峭地层，以及横向100 ～ 125 m的地层，在进行地形逆时偏移后同相轴聚焦情况得到改善]；（c）区域2的基准面逆时偏移图像；（d）区域2的地形逆时偏移图像[通过比较（c）和（d），发现在横向上340~360 m有一条陡峭同相轴，解释为正断层。在地形逆时偏移图像中，这条断层面聚焦得很好，但在基准面逆时偏移中很难被识别]。

《4.2.场地2》

4.2.场地2

图6为现代沙丘的复杂结构剖面图。剖面垂直于沙丘轴线，并穿过了两座7m高的沙丘。这两座沙丘位于基岩表面，且在沙丘表面只覆盖了一层薄沙，整体向左倾斜，右侧面为下风向。在基准面逆时偏移后，陡峭沙丘面和近似水平沙丘底界面相交处，基岩面错位严重，令其看起来像是位于沙丘底部的一个露头，这在地质上是不合理的。地形逆时偏移能将其正确地归位，最后发现基岩在一层薄沙之下。

《图6》

图6.（ a）沙丘中央使用100 MHz天线在一条测线上采集的预处理数据；（b）地形逆时偏移后沙丘内部地层和基岩表面成像图；（c）右边沙丘在基准面逆时偏移后的图像；（d）右边沙丘在地形逆时偏移后的图像（在基准面逆时偏移后，基岩表面位置在沙丘斜面和沙丘内部表面相交处错位严重，而地形逆时偏移能对内部陡峭的反射界面更好地聚焦）。

《5. 讨论和结论》

5. 讨论和结论

在起伏地表下复杂地层成像时，地震波在实际记录中的传播方式与在高程静校正中假设的近地表垂直传播方式明显不同。因此，图像在进行高程静校正逆时偏移后会严重扭曲，导致反射面同相轴不聚焦或者被偏移归位于错误的位置。而地形逆时偏移校正了地震波的传播路径，使得在速度横向突变时也能得到准确的图像。精确的水平和垂直地层位置对成像至关重要。

本文在建立爆炸反射界面模型时用到了二阶解耦形式的麦克斯韦方程，这种方法不用直接计算磁场，而是类似于Yee氏算法，只计算和保存电场值。研究中介绍的基于爆炸反射界面模型偏移算法，仅适用于一次波，如果数据中存在高幅多次波，则会产生严重的成像误差。最后，在逆时偏移时需要对速度模型进行平滑，以免波场延拓产生虚假反射。结合上述正演模拟和实例分析可以看出，在对崎岖地形和复杂地层精准成像时，叠后RTM是一种非常有效的处理方法。另外，精确的偏移结果需要精确的速度模型，当遇到速度结构比CPSD更复杂时，就需要更严格的速度估算方法。此时可利用多通道、多偏移距的采集方法，通过全波场双程波动方程实现叠前RTM，避免爆炸反射界面模型中的近似单向传播对成像效果的影响。而在其他测线上采集的雷达资料的叠前、叠后逆时偏移图像，有助于加深对这种密集型计算方法的优势和不足的理解。

《Acknowledgements 》

Acknowledgements

The Herbette Foundation at the University of Lausanne provided support for the development of the RTM algorithm. Compliance with ethics guidelines.

《Compliance with ethics guidelines 》

Compliance with ethics guidelines

John H. Bradford, Janna Privette, David Wilkins, and Richard Ford declare that they have no conflict of interest or financial conflicts to disclose.