ABSTRACT

The size, shape, and ice thickness of thermokarst lakes have obviously changed in the Tibetan Plateau permafrost region due to global climate changing. These changes have strong correlation with ice-permafrost and climate changing. Since thermokarst lakes have widely appeared and covered large area in Tibetan Plateau, it’s difficult to quickly monitor the change of ice thickness using ground-penetrating radar (GPR). Unmanned aerial vehicles (UAV) GPR can adapt to the complex area and realize the quickly and more times monitor, however, the ground surface as a strong reflection interface can generate the multiple waves, which will seriously impact on identification of ice thickness. In this paper, the researchers used the forward and reverse continuation method to eliminate the ground surface’s effect based on the widely appeared the ground surface’s effect for the UAV GPR data. Through comparison with the forward continuation and reverse continuation (FCRC) of zero-offset reverse time migration (RTM) and raw zero-offset RTM profile, the researchers demonstrate the FCRC method can eliminate the ground surface’s effect. From the zero-offset RTM profile of different elevation’s model, the researchers analyze the effectiveness of FCRC method for UAV GPR data with different flight altitude. The results of zero-offset RTM imaging show that the FCRC method can effectively eliminate the surface’s effect, however, with the increase of flight altitude, the resolution of migration imaging will decrease. From the result of forward modeling and zero-offset RTM imaging for models data, the researchers consider that the FCRC method can effectively eliminate multiple waves, better reflect the flat interface information, but cannot image the dip interface. The process of FCRC and flight altitude can decrease the UAV GPR resolution which will go against the high resolution detection; we still need make new progress for the process of UAV GPR data.

Keywords:Unmanned Aerial Vehicles Ground-Penetrating Radar, Multiple Wave, Forward Continuation and Reverse Continuation, Reverse Time Migration

无人机载探地雷达数据消除地面影响 的正反向延拓方法研究

鲁兴林，刘雅琼，宋翱，王玉洁，钱荣毅^*

中国地质大学(北京)，“地下信息探测技术与仪器”教育部重点实验室，北京

收稿日期：2017年1月1日；录用日期：2017年1月14日；发布日期：2017年1月18日

摘 要

随着全球气候的变化，青藏高原冻土区热融湖塘的大小、形状和冰层厚度在区域内存在明显的变化，热融湖塘的变化与区域内富冰冻土层和气候变化息息相关。高原区热融湖塘分布广，覆盖面大，现有的地面探地雷达技术很难快速准确地监测冰层厚度的变化。机载探地雷达能在复杂的地形区域实现快速、重复监测，但由于强反射地表面产生的多次波的干扰，严重影响冰层厚度的识别。本篇文章基于机载探地雷达数据中普遍存在地表面问题，应用正反向延拓方法压制地表面引起的多次波，消除地表面对地下构造识别的影响。通过比较正反向延拓前后的零偏移距逆时偏移剖面，论证了正反向延拓方法在消除地表面影响方面的优势；应用不同飞行高度的模型分析正反向延拓方法对不同飞行高度机载探地雷达数据的效果。零偏移距逆时偏移成像结果显示，正反向延拓方法能有效地去除地表面的影响，但随着飞行高度的增加，偏移成像的分辨率会降低。模型数据的正演模拟和零偏移距逆时偏移成像分析认为，正反向延拓方法能有效地压制多次波干扰，清晰地反映平界面特征，但对倾斜界面成像较差。延拓计算和飞行高度会影响偏移成像的分辨率，不利于高分辨率探测，还需要进一步改善。

关键词 :无人机载探地雷达，多次波，正反向延拓，逆时偏移

Copyright © 2017 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

热融湖塘是青藏高原冻土区典型的地貌特征，其分布范围广，覆盖面积大，并且大小、形状和冰层厚度都有明显的变化。热融湖塘的形成与区域内冻土层和气候变化紧密相关，准确地掌握青藏高原区热融湖塘大小、形状和冰层厚度的变化，对研究高原区冻土层的融化特征和区域内气候变化规律具有重要意义 [1] [2] 。地面探地雷达技术是目前探测冰层厚度最有效的地球物理方法，凭借探测速度快，分辨率高和探测结果直观成像的优势，广泛应用于环境生态保护领域 [3] [4] [5] 。但地面探地雷达探测效率低，覆盖范围小，不适宜在高原区开展大面积探测作业。无人机载探地雷达的问世，可有效地弥补地面探地雷达的缺陷 [6] [7] 。

前人研究工作表明，机载探地雷达能有效地探测冰层厚度。美军寒地研究及工程实验室 [8] (CRREL)采用直升机搭载常规探地雷达天线的方法，在阿拉斯加和南极洲附近开展了大量的探测冰层和冰内洞穴的工作，取得显著的成果。美国德克赛斯大学奥斯汀分校的地球物理研究所(UTIG)研制的固定翼飞机搭载的探地雷达，用于研究Vostok湖子冰川的层理，能清晰地反映上千米冰下的反射信息 [9] 。国内十二五期间，将无人机搭载的探地雷达应用于黄河内蒙古包头段的防凌破冰工作，采用机载探地雷达获取河道冰层厚度的空间分布数据 [10] 。

青藏高原区地域广阔，无人机载探地雷达能有效地解决高原区热融湖塘冰层厚度的快速探测问题。由于机载探地雷达天线放置在具有一定高度的飞机上，其数据处理方法与地面探地雷达数据处理方法明显不同。针对机载探地雷达数据的处理问题，很多学者已提出不同的数据处理方法。Catapano等 [11] [12] 参考二维半空间几何关系和散射体模型，提出线性反演散射成像方法处理机载探地雷达数据。傅磊等 [8] 用三维正演模拟方法，模拟不同天线极化方向、天线高度及地表粗糙度时的机载探地雷达剖面，分析各个因素对目标体成像的影响，并提出用逆时偏移方法处理机载探地雷达数据。屈乐乐 [13] 针对机载探地雷达探测过程中电磁波在空气与地表面发生的折射问题，提出一种基于非均匀快速傅里叶变换(NUFFT)的快速BP成像算法，用于补偿地表面的折射损失。机载探地雷达数据采集时，雷达天线辐射出的电磁波经过空气层的扩散和传播，经地表层透射进入地下，遇到地下目标体后产生反射并向上传播，反射波信号被接收天线接收记录。空气与地表面为电磁波的强反射界面，会产生自由表面多次波的干扰，前人的数据处理研究工作都未能较好地处理地表面产生的多次波问题。多次波会在偏移成像过程中产生虚假成像，严重影响对浅地表结构特征的识别。

为提高浅地表层机载探地雷达数据成像分辨率和可靠性，需要消除强反射地表面的对数据处理的负面影响。本文借鉴海洋地震数据处理中消除海底强反射界面的方法 [14] [15] ，采用电磁波TM方程正反向延拓法，经两次延拓压制自由地表面产生的多次波，消除强反射地表面的影响，再用逆时偏移方法处理延拓后的数据 [16] 。通过比较正反向延拓前后的逆时偏移剖面，论证了正反向延拓方法在压制自由地表面多次波方面的优势。应用不同飞行高度模型的正反向延拓前后的逆时偏移剖面，分析正反向延拓方法对不同飞行高度机载探地雷达数据的影响。通过模型数据的正演模拟和偏移成像分析认为，正反向延拓能有效地压制自由表面多次波，但延拓计算和飞行高度会降低数据的分辨率，不利于高分辨率机载探地雷达探测，还需要进一步改善。

2. 正反向延拓原理

正反向延拓法压制多次波，其基本原理是通过两次延拓将电磁波在空气中的传播过程替换成电磁波在地下第一层介质中的传播过程，从而消除强反射地面产生的多次波的干扰。正反向延拓计算过程采用Maxwell方程推导出的TM波方程 [17] ，其方程形式为(式1)：

(1)

式中的表示介质介电系数，单位为法拉/米(F/m)，表示磁导系数，单位为亨利/米(H/m)，表示电导率，单位为西门子/米(S/m)，，分别为磁场的x和y分量，为电场分量。

图1为机载探地雷达电磁波传播示意图，发射天线A点发射电磁波，接收天线B点接收反射点O

图1. 机载探地雷达电磁波传播路径及正反向延拓示意图

点的反射信息。将接收天线接收到的电磁波信号应用TM波方程推导出的逆时有限差分法反向延拓到反射点O，并在地表面上A点的正下方一点C点接收反向延拓的信号，反向延拓时用的是空气层的相对介电常数。然后在地表面C点接收到的反向延拓的信号应用TM波方程推导出的正向有限差分法正向延拓至反射O点，并在接收天线B点接收到正向延拓信号，正向延拓时用的是近地表层的相对介电常数。

3. 逆时偏移方法

逆时偏移由正演、反向延拓和成像三个步骤组成：(a) 正演是波场沿着时间轴的正方向传播到最大时刻并将运算结果保存下来；(b) 反向延拓是从接收波场的最后一个采样点开始，沿着时间轴的负方向延拓至零时刻，并读取正演过程相应时刻的波场值；(c) 应用合适的成像条件对正演波场和反向延拓波场成像，得到真实的地下构造信息。逆时偏移成像过程中，选取零延迟互相关成像条件(式(2)) [18] ，成像后采用拉普拉斯滤波去除近地表处的低频噪音 [19] 。

(2)

其中image (x, z)是点(x, z)的成像结果：S (x, z, t)表示时间域的正演波场；R (x, z, t)表示时间域的反向延拓波场。正演模拟是逆时偏移方法的基础，笛卡尔坐标系TM波方程是由麦克斯韦方程组解耦出的显示的一阶偏导方程 [17] 。

依据探地雷达数据采集方式的不同，逆时偏移方法分为零偏移距和多偏移距两种。考虑到机载探地雷达数据采集时，激发天线和接收天线间距离很小，可以近似为零偏移距采集，所以本文中采用零偏移距逆时偏移方法对模型数据做偏移成像。

4. 模型实验

4.1. 正反向延拓法压制多次波

为论证正反向延拓法在压制多次波方面的优势，设计两层台阶模型(图2)，正演模型分三层(图2(a))，其相对介电常数分别为1.0、6.0、4.0，电导率均为0。图2(b)为正反向延拓后逆时偏移所用的相对介电常数模型，空气层充填为浅地表层第一层的相对介电常数。两个模型(图2(a)，图2(b))的有限差分网格大小为200 × 200，横向和纵向网格间距均为0.02 m，雷克子波的主频为900 MHz，时间步长为0.02 ns，正演数据记录长度为60 ns。

应用有限差分正演模拟方法沿着距地表8 m的飞行高度正演出零偏移距剖面(图3(a))，再将其反向延拓到地表面，得出反向延拓零偏移距剖面(图3(b))，反向延拓过程采用空气层的相对介电常数。最后再将反向延拓剖面正向延拓到实际的飞行高度，得出正向延拓零偏移距剖面(图3(c))，正向延拓过程采用浅地表层的相对介电常数。正演剖面显示(图3(a))，时间为25 ns 和 40 ns为第二层平界面的反射波(图3(a)，红箭头A和B)，C和D为第二层界面产生的自由表面多次波。反向延拓后，第二层平界面的反射波和自由表面多次波的时间都明显减小(图3(b)，红箭头A、B、C、D)。正向延拓剖面能清晰地显示第二层平界面的反射波(图3(c)，红箭头A和B)，自由表面多次波得到有效地压制(图3(c)，C和D)。正反向延拓过程能有效地去除自由表面多次波的影响，但延拓过程会降低数据的分辨率(图4)。

应用有限差分逆时偏移方法计算得出正反向延拓前后的逆时偏移剖面(图5)，延拓前相对介电常数模型用实际正演的相对介电常数模型(图2(a))，延拓后相对介电常数模型用去除地表强反射界面的相对介电常数模型(图2(b))。两个偏移剖面都能清晰地反映平界面的特征，但相比于延拓后的偏移剖面，延拓前的偏移剖面存在很多虚假成像(图5(a)，红箭头A、B、C、D)，而延拓后的偏移剖面能有效地弱化虚假成像(图5(b)，红箭头A、B、C)，能较好地压制自由表面多次波的干扰(图5(b)，红箭头D)，能更清晰地突出平界面的特征。

图2. 相对介电常数台阶模型，(a) 正演模拟及正反向延拓前逆时偏移所用相对介电常数模型；(b) 正反向延拓后逆时偏移所用的相对介电常数模型

图3. (a) 零偏移距探地雷达剖面；(b) 反向延拓探地雷达剖面；(c) 正向延拓探地雷达剖面

4.2. 不同飞行高度逆时偏移成像分析

为分析不同飞行高度时正反向延拓法对逆时偏移成像的影响，设计三层凹陷模型(图6(a))，其相对介电常数

图4. 正反向延拓前后机载探地雷达剖面第150道数据(图3，红线) 频谱分析示意图

图5. (a) 逆时偏移剖面；(b) 正反向延拓的逆时偏移剖面

图6. 相对介电常数凹陷模型，(a) 正演模拟及正反向延拓前逆时偏移所用相对介电常数模型；(b) 正反向延拓后逆时偏移所用的相对介电常数模型

分别为1.0、6.0、4.0，电导率均为0。图6(b)为正反向延拓前后逆时偏移所用的相对介电常数模型，空气层充填为浅地表层第一层的相对介电常数。设计飞行高度为5 m、10 m、15 m和20 m，分别沿着飞行高度正演模拟零偏移距剖面，并分别用正反向延拓前后的数据做逆时偏移成像，分析不同飞行高度对偏移成像的影响。整个计算过程，横向和纵向网格间距均为0.2 m，雷克子波的主频为50 MHz，时间步长为0.2 ns，正演数据记录长度为900 ns。

机载探地雷达飞行高度为5 m时，正反向延拓前后偏移剖面显示(图7(a-1))和(图7(b-1))，正向延拓前偏移

图7. (a-1) 逆时偏移剖面(飞行高度为5 m)；(b-1) 正反向延拓的逆时偏移剖面(飞行高度为5 m)；(a-2) 逆时偏移剖面(飞行高度为10 m)；(b-2) 正反向延拓的逆时偏移剖面(飞行高度为10 m)；(a-3) 逆时偏移剖面(飞行高度为15 m)；(b-3) 正反向延拓的逆时偏移剖面(飞行高度为15 m)；(a-4) 逆时偏移剖面(飞行高度为20 m)；(b-4) 正反向延拓的逆时偏移剖面(飞行高度为20 m)

剖面在平界面底部存在自由表面多次波的成像干扰(图7(a-1)，红箭头A和E)，在凹槽底部存在两个能量较强的交叉反射轴(图7(a-1)，红箭头C)，对凹槽左倾界面的成像为不连续的反射轴(图7(a-1)，红箭头B)，并不能清晰地显示凹槽右倾界面的特征(图7(a-1)，红箭头D)，另外在凹槽下方出现明显的干扰波(图7(a-1)，红圈)。正反向延拓处理后，自由表面多次波的干扰得到较好的压制(图7(b-1)，红箭头A和E)，能较好地去除干扰波的影响，如在凹槽底部(图7(b-1)，红箭头C和D)，凹槽下方(图7(b-1)，红圈)和倾斜界面(图7(b-1)，红箭头B和D)，能更清晰地反映平界面的特征。随着飞行高度的增加，凹槽倾斜界面处的成像效果越来越差(图7(a)，红箭头B和D)，凹槽界面下方的干扰波能量越来越大(图7(a)，红圈)。应用正反向延拓偏移方法能较好地去除自由表面多次波的干扰(图7(b)，红箭头A和E)和凹槽下方的能量较强的干扰波(图7(b)，红圈)，但正反向延拓偏移成像方法并不能对倾斜界面清晰成像(图7(b)，红箭头B和D)，还需要做进一步改进。

5. 结论

本文主要开展机载探地雷达数据处理研究工作，解决机载探地雷达数据中普遍存在的强反射地表面问题，应用正反向延拓方法消除强反射地表面引起的自由表面多次波的干扰，并用零偏移距逆时偏移方法分析正反向延拓的处理效果。根据模型数据的正演模拟和逆时偏移成像分析，总结如下：

(1) 电磁波TM方程正反向延拓方法，能有效消除地表面引起的自由表面多次波的干扰，使偏移成像真实可靠；

(2) 零偏移距逆时偏移方法是地面探地雷达数据处理中精度较高的偏移方法，但对机载探地雷达数据处理时，只能反映平界面特征，并不能对倾斜界面清晰成像；

(3) 正反向延拓计算和机载探地雷达飞行高度的增加，都会降低偏移成像的分辨率，不利于高分辨率探测。

机载探地雷达可在高寒、危险系数高，覆盖范围广的高原冻土区内开展重复作业，具有高效、精准的优势，能有效地弥补地面探地雷达数据量小，采集难度大的劣势。正反向延拓方法能有效地压制地表层引起的自由表面多次波，但在数据分辨率方面仍需要进一步改进。

基金项目

国家重点基金委项目(2145803)。

文章引用

鲁兴林,刘雅琼,宋 翱,王玉洁,钱荣毅. 无人机载探地雷达数据消除地面影响的正反向延拓方法研究
Forward and Reverse Continuation of Unmanned Aerial Vehicles (UAV) Ground-Penetrating Radar (GPR) Data to Eliminate Ground Surface’s Effect[J]. 计算机科学与应用, 2017, 07(01): 27-35. http://dx.doi.org/10.12677/CSA.2017.71004

参考文献 (References)