复杂环境中探地雷达数据的逆时偏移成像

工程4（2018）661研究应用地球物理学-文章复杂环境下从起伏地形到成像探地雷达数据的John H.布拉德福德a，贾纳·普里维特a，大卫·威尔金斯a，理查德·福特ba美国博伊西州立大学地球科学系，ID 83725-1535b美国犹他州奥格登韦伯州立大学地球科学系，邮编84408-2507阿提奇莱因福奥文章历史记录：2017年12月20日收到2018年5月10日修订2018年9月7日接受在线发售2018年保留字：探地雷达逆时偏移沙丘幅度分析A B S T R A C T在探地雷达（GPR）成像中，探测深度通常与地表地形的变化程度相同。在这种情况下，在应用高程静校正后从基准面进行偏移时，偏移会失败。我们介绍了一种逆时偏移（RTM）算法的基础上的二阶解耦形式的麦克斯韦方程组，它只需要计算的电场。波场外推是直接从采集表面计算，而不需要基准面。在一个合成的案例研究中，该算法显着提高了图像精度的处理序列，其中迁移后进行高程静校正。此外，我们在美国犹他州的珊瑚粉沙丘（CPSD）获得了一个现场数据集这些数据是在崎岖的地形上获得的RTM算法显着改善雷达深度图像在这种具有挑战性的环境。©2018 The Bottoms.由爱思唯尔有限公司代表中国工程院和高等教育出版社有限公司出版。这是一个在CC BY-NC-ND许可下的开放获取文章（http://creativecommons.org/licenses/by-nc-nd/4.0/）中找到。1. 介绍逆时偏移（RTM）是计算昂贵的，然而，在过去的二十年中，硬件的进步，使人们有可能经常应用这种方法的许多近地表地震和探地雷达（GPR）的应用。RTM通常应用于大型勘探地震数据集，是首选的成像工具。RTM使用波传播的全部物理过程，因此避免了其他偏移方法的假设。RTM能准确处理反射层倾角和速度非均质性。此外，衰减和传输损耗可以并入RTM框架[1]。在这项研究中，我们专注于RTM在叠后域。以前的一些学者已经将RTM算法应用于GPR数据。在早期的例子中，Fisher et al.[2]利用声波算法和爆炸反射体模型对GPR数据进行叠后偏移Sanada和Ashida[3]根据麦克斯韦方程组开发了一种RTM叠后算法地质雷达图像的构建是具有挑战性的地形变化在同一顺序的调查深度*通讯作者。电子邮件地址：jbradford@mines.edu（J.H.布拉德福德）。[4]的文件。当应用常规处理策略时，可能会导致严重的误差，例如高程静校正之后进行叠后偏移（图1）[4];在这些条件下必须采用替代处理策略[5]。 Heincke等人[6]描述了一种替代的地形偏移方法，该方法利用衍射扫描和随后的形态图像处理来识别复杂地形中的陡倾裂缝。Zhu和Huang[7]通过应用经典Yee算法[8]来解决零偏移距数据的RTM问题（使用通常的爆炸反射体概念），从而超越了声波他们通过将其应用于沙丘地层模型生成的合成数据来测试他们的算法。Zhu和Huang[7]得出结论，RTM可以改善沙丘中常见的陡峭地层的成像，并建议使用野外数据进行进一步研究。在这项研究中，我们提出了一个RTM叠后成像方法在复杂地形的爆炸反射体的概念的基础上。直接从地形表面进行向下波场外推，从而避免了对高程静校正或基准面校正的需要。此外，该算法是基于麦克斯韦方程组的波动方程解，并可以考虑复杂的速度和电导率分布。我们测试的算法与合成的例子，并与现场数据采集的沙丘，包含陡峭的倾斜，复杂的https://doi.org/10.1016/j.eng.2018.09.0042095-8099/©2018 THE COMEORS.由爱思唯尔有限公司代表中国工程院和高等教育出版社有限公司出版。这是一篇基于CC BY-NC-ND许可证的开放获取文章（http://creativecommons.org/licenses/by-nc-nd/4.0/）。可在ScienceDirect上获得目录列表工程杂志主页：www.elsevier.com/locate/eng662J.H. Bradford等人/工程 2018年第4期第661¼¼er4r52e在爆炸反射体模型中，场近似为单向传播，模型速度为真实速度v0v=2，其中，初始参数表示爆炸反射体模型中使用的有效值。从等式（2），很明显，爆炸反射器的解决方案方程。（1）需要变量替换e04e。为了确保初级到达在单向传播期间经历正确量的衰减，等式（3）要求：r0rl0zrrl02z4Fig. 1. (a)沿地形表面采集的零偏移剖面在基准面校正后的行进路径示意图：垂直虚线显示行进路径，这是在标准高程静校正中隐含假设的;实线表示正确的行进路径;假设的路径是其中z是到反射器的深度等式的右侧（4）是双向传播的衰减，而左手边是爆炸反射器模型中单向传播的衰减。在取消公因子后，重新排列，并注意到e0¼4e，我们发现如下：re0在高程静校正后，近似导致从基准面偏移失败（改编自参考文献[4]）。(b)沿（a）中所示的采集表面采集的合成GPR数据。εr是相对介电常数。地层学加上地形变化，与目标深度相同。2. RTM算法RTM在概念上很简单：记录的波场通过求解麦克斯韦方程组而被外推到地下以雷达图作为边界条件现在将初始本构参数e0和r0代入Eq. （1）得到了爆炸反射体模型的修正波动方程。我们使用二阶精度有限差分算法求解波动方程，并使用完全匹配层（PML）吸收边界条件，如公式所示：Zhou等人[9];这导致在边界区域求解耦合各向异性方程组。例如，在z方向上，Zhou等人的“的公式，在边界区域求解的微分方程组是az@D1zrzD1z¼@E在地面上以相反的时间顺序。在堆栈后的情况下，当时钟向后计数到零时，满足成像条件@te0@D2zrz@z@2D1z时间-也就是说，当所有记录的能量都被注入时，的z @te0D2z¼@t@zð6Þ进入模型并向下传播到其原点。@2E@E@2E@2E@D2z该RTM算法中包含两个首先，速度场必须平滑，以防止在外推过程中的人工反射。我们在小波的主频率处使用两个波长上的三角形平滑器;这种平滑程度以最小的精度成本最第二个警告是，爆炸反射体模型仅正确处理了一次反射;如果数据中存在高振幅多次波，则会导致在这里，我们解决了麦克斯韦方程组的解耦，二阶微分形式。在二维（2D）系统中，电场垂直于像平面极化（横向电（TE）模式），这些方程被简化为：阻尼标量波动方程：l0e0@t2l0r0@t¼@z2@x2@t其中，辅助变量D1z和D2z由等式2定义（六）、变量rz和az是PML吸收参数。我们在MATLAB中实现的算法，利用矩阵计算，以获得计算效率，利用多核/处理器系统。波场外推使用Eq.（1），但具有负的时间步长。记录的波场作为边界条件沿记录或地形表面以相反的时间顺序输入。由于该算法是在正方形网格上计算的，因此沿表面的网格不连续性会引入一些误差。由于这些不连续性比波长小得多，因此引入的误差最小（数据中最高频率处的k/10比如说，le@2Elr@E@2E@2E对于100 MHz Ricker小波，1@t2@t¼@x2@z2其中E是电场，l是磁导率，t是时间，x和z是空间变量。变量e和r是介电常数和电导率。在我们的算法中，我们假设磁导率l等于自由空间的磁导率l0。我们利用低损耗近似，其中损耗角正切r=εxε1，其中x是波的角频率在这种情况下，速度v和衰减系数a如下：数据大约为2.5f0或250 MHz。一个典型的马克西-μ m的相对介电常数εr= 9，网格间距为4 cm，与在100 MHz主频处的1 m波长一致。计算网格没有散射是一个很好的指示，表明近似是合理的。3. 合成实施例3.1. 合成实施例1：简单速度结构1l0a¼rrl0ð2Þð3Þ我们首先用简单的分层速度结构研究大地形变化的影响（图1（a））。 该模型具有较大的表面形貌变化相对于调查的深度，这是在GPR调查中经常遇到的情况。gation。 正弦表面形貌具有谷-峰v¼p2e0级2er0¼2J.H. Bradford等人/工程4（2018）661663高度为4 m，不规则反射层深度为1 ~ 5 m。在模拟的探地雷达数据中，存在一个单一的反射界面，它产生了一个复杂的波场（图1）。图1（b））示出了具有多个向斜的反射体的三重特性。我们考虑三种处理顺序：① RTM假设在平面上采集数据，然后进行高程静校正（ 图 2 （ a ）） ; ②高程 静校正， 然后 从基准面进行RTM（图 2（a））。 2（b））;或③ RTM从地形（图。 2（c））。在没有地形校正的情况下迁移数据会导致很大的误差;图像聚焦不良，反射位置不正确（图2（a））。数据似乎过偏移了7至10 m，13至18 m之间的背斜似乎是平坦的。应用标准高程校正，然后进行偏移，可以改善结果，但图像重建仍然包含较大的误差（图2（b））。请注意向斜中心沿剖面距离15 m。它似乎是过度偏移的，这可能导致错误的结论，即速度模型是不正确的（图1）。 2（b））。图像在6到8米之间的区域没有正确聚焦该区域包含地形的最低点和反射器上的陡峭倾斜的短段。聚焦不良导致在该区域中出现两个明显的来自地形的RTM正确地聚焦图像，并将反射边界放置在正确的深度（图1）。2（c））。3.2. 合成实施例2：大横向速度梯度RTM的一个主要优点是它能够处理大的横向速度对比。为了测试我们的算法对横向速度突变和相对于调查深度的严重地形的响应，我们构建了模型，示在图3（a）中。该模型具有与示例1中的模型相同的几何形状，但是我们在模型左侧的第一层中添加了介电常数的突然增加（速度的减小）。所得合成探地雷达数据类似于图1的数据。 1（b），图二、 合成数据来自Fig. 1（b）具有三种 不同的 处理序列。(a)RTM假设在平面上采集数据，然后进行高程静校正;（b）高程静校正，然后根据基准面进行RTM;（c）根据地形进行RTM。（a）和（b）中的图像聚焦不良，而来自形貌的RTM（c）正确地聚焦了图像，并将反射边界放置在正确的深度处，如黑色实线所示。但是在模型的左侧具有一些显著的传播时间延迟，其中速度低于第一模型。本文只考虑两种处理方案：①常规流程高程静校正后再进行RTM，②从地形图进行RTM对于高程静校正，采用较高速度的材料作为置换速度。在高程静校正后，使用正确的速度模型从基准面进行RTM，但图像聚焦不佳，反射波位置不正确（图1）。 3（c））。这张照片与图相似。第2段（b）分段。 来自地形的RTM精确地重建了图像，即使在覆盖层中横向速度对比大的区域也能产生反射边界的清晰图像（图1）。3（d））。值得注意的是，反射器上的陡峭倾斜段（中心在7米的距离处）难以成像。在这种情况下，垂直速度边界折射的大部分能量反射从陡峭的部分远离表面。最终结果是在该位置处的重建图像中的低振幅还值得注意的是，模型左侧的反射系数较小，导致沿界面10 m距离处的振幅较低4. 现场测试：美国犹他州珊瑚粉沙丘珊瑚粉沙丘（CPSD）位于犹他州南部CPSD坐落在纳瓦霍砂岩上，图3.第三章。（ a）第一层10米处横向对比鲜明的合成模型;（b）从（a）生成的合成数据;（c）在高程静校正之后从基准面的RTM在聚焦图像方面做得很差;（d）从地形的RTM即使在具有大的横向速度梯度的情况下也能精确地重建图像，并且避免了基准面校正的需要。664J.H. Bradford等人/工程 2018年第4期第661被Sevier正断层一分为二，该正断层也形成了沿低沙丘区（LDF）东部边界的基岩悬崖。为了检验断层控制的地形沿下伏基岩表面控制沙丘的形成和几何形状的假设，我们进行了地质雷达研究。我们的主要目标是绘制我们收集了超过20公里的配置文件沿25个断面与50和100 MHz的天线。使用连续全球定位系统（GPS）进行高程控制，并在后处理中进行全球定位系统基站位于所有样带的10公里范围内。通过将探地雷达采集时钟与GPS时钟同步，将GPS和探地雷达位置联系起来探地雷达信号的穿透性非常好，我们在某些位置记录了深度超过35米的这为现代沙丘和古沙丘地层提供了出色的图像沿某些断面的露头和/或浅钻孔虽然探地雷达信号质量很好，但地形和地层的复杂性给现场的成像带来了挑战。一些剖面的地表地形变化超过25米，持续梯度大于30°。这个数据集提供了一个很好的现场例子来测试我们的RTM算法。试验了两种处理流程：①高程静校正后的RTM处理流程和②地形RTM处理流程。每个流程的处理步骤见表1。关键的是，我们在一个恒定的介电常数，我们确定通过检查衍射双曲线的崩溃后，在几个不同的许可证值迁移。使用相同的介电常数在整个勘测区域产生了良好的由于唯一的重要复杂因素是地形和地层，因此该数据集是从地形方法评估RTM的理想现场示例4.1. 位置1图图4和图5展示了具有一些特别有趣的特征的典型轮廓。一般来说，我们看到现代沙丘砂位于基岩表面的顶部，基岩由古老的石化沙丘组成。基岩常暴露在沙丘内部。区域1是基岩地层内的陡倾床型带。从基准面RTM，底模差表1高程静校正后RTM与地形RTM对比的处理流程从基准面的RTM（高程静校正后）从地形的几何形状：来自差分校正的GPS的迹线（x，y，z）几何形状：来自差分校正的GPS的迹线（x，y，z）时间零点校正：基于直接无线电波时间零点校正：基于直接无线电波时域带通滤波器：1/2主周期高通，3×主周期时域带通滤波器：1/2主周期高通，3×主周期按t2缩放的增益按t高程静校正：基准面为剖面最高高程，置换速度为0.11m·ns-1从基准面的RTM： er7.4（v= 0.11m·ns-1）的常介电常数偏移（v= 0.11m·ns-1）见图4。 CPSD测量的剖面示例，显示（a）偏移前的预处理数据，以及（b）地形RTM后的数据。J.H. Bradford等人/工程4（2018）661665图五、在区域1中，将（a）来自基准面的RTM与（b）来自地形的RTM进行比较，我们注意到基岩地层中20至30 m深度处距离小于50 m的陡峭反射的不同位置，以及从地形迁移时100至125 m距离之间的改进聚焦在区域2中，比较基准面的RTM（c）和地形的RTM（d），在340和360 m之间有一个向右倾斜的事件，我们将其解释为正断层。从地形上用RTM法可以很好地聚焦断层面，但从基准面偏移时很难解释例如，左边的斜坡似乎穿过了一个更老的、更深的表面，深度为30利用地形的RTM，陡倾地层更好地聚焦，并且看起来定位正确区域2包含一个特征，我们将其解释为基岩内的正断层。从基准面进行RTM，断层面聚焦不好，难以解释地层和断层的相对几何形状利用地形的RTM，断层面被很好地聚焦，我们可以看到复杂的压裂几何形状。请注意，这条断层似乎并没有抵消地表附近的侵蚀接触，因此在现代不太可能活跃4.2. 位置2图6显示了我们在现代沙丘结构中观察到的一些复杂性。剖面图穿过两个7米高的沙丘，方向垂直于沙丘轴。沙丘直接位于基岩表面，沙丘内部区域仅有一层薄薄的沙子覆盖。考虑剖面右侧沙丘的陡峭、在沙丘的陡峭面与近似水平的沙丘内表面相交的地方，基准面的RTM严重错位基岩表面，使基岩看起来就在沙丘底部上方露出;然而，基岩绝对不是见图6。(a)用100兆赫天线从沙丘区中心部分的一条线上采集的预处理数据;（b）从地形上采集的RTM后的数据，其中沙丘的内部地层和基岩界面清晰;（c）从基准面上采集的RTM后显示的右手沙丘;（d）从地形上采集的RTM后的同一沙丘。在从沙丘底部到沙丘表面的过渡处，从基准面偏移时，基岩表面明显错位此外，陡峭倾斜的内部反射器在RTM之后从地形更好地聚焦。666J.H. Bradford等人/工程 2018年第4期第661露头地形的RTM正确地定位了基岩表面，我们发现基岩正好在一层薄薄的沙子覆盖层之下。5. 讨论和结论当对崎岖地形下的复杂地层成像时，记录的波场的运动学显著偏离标准高程静校正中隐含的垂直近地表传播的假设。在这样的校正之后，从基准面进行RTM会显著地使图像失真，从而导致反射的不良聚焦和不正确定位。RTM从地形正确对待波场运动学，并产生一个准确的图像，即使在存在大的横向速度梯度。精确的水平和垂直测量对于生成精确的图像至关重要。爆炸反射器模型，我们在这里实现利用二阶解耦形式的麦克斯韦方程组。由于磁场不是直接计算的，如在Yee算法中，这种方法具有仅计算和存储电场的计算优势。本研究中讨论的成像算法中使用的爆炸反射体模型仅适当地考虑了一次到达;如果数据中包含大的多次波，则会发生显著的成像误差。最后，RTM程序需要平滑速度模型，以避免波场提取期间的人工反射有了这些警告，这里介绍的案例研究表明，当在复杂和恶劣的环境中生成准确的图像至关重要时，叠后RTM是一种有价值的工具。准确的偏移结果依赖于准确的速度模型。当速度结构比CPSD遇到的更复杂时，将需要更严格的速度估计方法对于多偏移距数据，可以利用具有完整双向波动方程的叠前RTM，从而避免爆炸反射体模型中所需的单向传播近似。叠前和叠后RTM GPR成像的其他现场实例将提高我们对优势的理解以及这种计算密集型成像技术的局限性。确认洛桑大学的Herbette基金会为RTM算法的开发提供了支持。遵守道德操守准则John H.布拉德福德、詹娜·普里维特、大卫·威尔金斯和理查德·福特宣布他们没有利益冲突或财务冲突需要披露。引用[1] Deng F，McMechan GA.真振幅叠前深度偏移地球物理2007;72（3）：S155-66.[2] Fisher E，McMechan GA，Annan AP，Cosway SW.单通道探地雷达剖面逆时偏移实例。地球物理1992;57（4）：577-86.[3] Sanada Y，Ashida Y.探地雷达数据的成像算法In：Powers MH，CramerL，BellRS，editors.地球物理学在工程和环境问题中的应用研讨会论文集; 1999年3月14 -18日;美国加利福尼亚州奥克兰。Wheat Ridge：环境与工程地球物理学会，1999年。p.565- 73[4] Lehmann F，Green AG. georadar数据的地形偏移：对采集和处理的影响。地球物理2000;65（3）：836-48.[5] 放大图片作者：J. J.地形起伏地区三维地质雷达测量的采集和处理策略。地球物理2005;70（6）：K53-61.[6] [10]张晓刚，张晓刚，张晓刚.基于半均衡的地形偏移（SBTM）：一种在三维地质雷达数据中识别断裂带的方法。Near Surf Geophys2006;4（2）：79-88.[7] 朱伟，黄强.地质雷达逆时偏移在沙丘内部构造探测中的应用在：探测地球物理学家协会第85届年度国际展览和会议的会议记录;2015年10月18日至23日;新奥尔良，洛杉矶，美国。塔尔萨：勘探地球物理学家协会; 2015年。p. 2269- 74[8] Yee K各向同性介质中麦克斯韦方程初边值问题的数值解。IEEE Trans Propag1966;14（3）：302-7.[9] 周东，黄文萍，徐春林，方东刚，陈波。标量时域有限差分法的完全匹配层边界条件。IEEEPhotonics Technol Lett 2001;13（5）：454-6.