1 引言

探地雷达(Ground Penetrating Radar)简称GPR，是一种利用电磁波进行无损探测的浅层地球物理技术，被广泛应用于考古学，冰川学，矿产勘查，以及环境与工程等领域.Norhidayahti(2012)等利用探地雷达在马来西亚的仁岭进行考古研究； Arcone(1991)研究了南极洲冰雪覆盖的3.2 km深度内的地质结构；Moore等(1999)用50 MHz和200 MHz雷达天线对挪威斯瓦尔巴特群岛南部的多热型冰川进行了高精度勘察；Singh(2007)利用探地雷达对大型矿砂矿床资源量进行评估.对于植被严重覆盖的区域，或人类无法到达的危险地区如战场、雷场等，常规的探地雷达则显得无能为力，然而机载探地雷达(Damm et al., 2006；刘四新等，2012)(Airborne GPR)却是一种有效的探测手段.

目前，机载探地雷达的发展和其在工程上的应用仍然受很多因素制约.美军寒地研究及工程实验室(CRREL)通过在直升机上悬吊传统探地雷达系统(Arcone，1991)，在阿拉斯加及南极洲对冰层厚度及冰层中裂隙进行了探测研究.德国地球科学和自然资源研究所(BGR)开发了一套特殊的机载探地雷达系统(Krellmann et al., 2008)，它被安装在AS350直升机上，角反射器天线作为发射和接收天线，操作频率为150 MHz，该雷达系统成功应用于地质调查.美国德克萨斯大学奥斯汀分校地球物理研究所(UTIG)开发了固定翼机载探地雷达系统(Holt et al., 2006)，该平台包括了重力仪、磁力仪、激光高度计、导航系统、数码相机以及冰雷达，该系统为窄带系统，特点是速度快、功率大.

由于机载探地雷达的天线距离地表较高，其电磁波传播要比常规探地雷达复杂.常规探地雷达天线通过与地面耦合向地下发射电磁波，而机载探地雷达天线首先与空气耦合，电磁波能量在空气中由于几何扩散而衰减，到达地表时，一部分能量反射(散射)回空气中，剩下的能量继续向地下传播，遇到地下介质不连续面便产生反射波，该反射波经过地层衰减以及几何扩散，最终到达接收天线.可见，相比于常规探地雷达，机载探地雷达的探测机理变得更为复杂.

影响机载探地雷达探测效果的因素主要有天线极化方向、天线高度以及地表粗糙度等，为了研究这些因素与机载探地雷达探测效果之间的关系，采用三维时间域有限差分(3D-FDTD)方法对沙漠地区地下空洞掩体机载探地雷达探测进行数值模拟，并分析各因素对探测效果的影响.机载探地雷达探测过程中，起伏的地表会导致地表散射杂波的存在，这些散射杂波的能量一般较强，在地下目标体埋深较浅的情况下，目标体信号会被地表杂波掩盖.为了消除地表杂波的影响，张蓓等(2005)提出子空间地表杂波抑制方法，并取得了一定的效果，本文使用逆时偏移成像(Chang and McMechan， 1994；Wu et al., 1996；Sun et al., 2006；龙桂华等，2011；王童奎等，2012；徐兴荣等，2012；孙小东等，2012；王保利等，2012；胡明顺等，2013；张岩和吴国忱，2013)的方法来消除地表杂波的影响.

传统探地雷达探测中，与天线接触的地表一般比较平坦，而机载探地雷达探测过程中遇到的地形情况比较复杂.为了模拟粗糙不平的起伏地表，引入高斯随机粗糙表面(Kobayashi et al., 2002；法文哲和金亚秋，2010)，它的功率谱表示如下：

图 1

Fig. 1

图 1 高斯随机粗糙面 (a)kσ1=0.5;(b)kσ2=1;(c)kσ3=2. 色标为均方根高度，单位m. Fig. 1 Gaussian r and om roughness surface

通过对典型地质模型进行雷达波的正演模拟，可以了解电磁波在介质中的传播过程，有助于对雷达数据的解释.从麦克斯韦微分方程组出发，Yee(1966)提出了时间域有限差分求解法，它将麦克斯韦旋度方程组进行差分离散，在时间上迭代推进求解空间中每一个时刻的电磁场值.

本文利用三维时间域有限差分法(刘四新等，2005； 2007；2011；2012；Taflove et al., 2005)对起伏的沙漠地区地下空洞掩体机载探地雷达探测进行了数值模拟.图 2a是建立的三维地质模型，其空间尺度在x方向为116 m，在y方向为60 m，在z方向为30 m.用度量值kσ₂=1高斯随机粗糙表面代 表起伏的沙漠表面，沙漠表面所在的平均高程为15 m，即干沙的厚度约为15 m，其相对节电常数为4，电导率为0.1 mS·m^-1.在干沙地层内部有两个沿y走向空洞掩体，分别用充满空气的两个圆柱体表示.图 2b是y=15 m处的二维切片模型，从图可知，两个 空洞的xz截面中心坐标分别为(48，10.5)和(68，7.5)m，空洞的截面半径均为1.5 m.图中虚线代表天线的飞行轨迹，即合成数据的记录位置，其距离平均地表的平均高度约12 m，水平方向30~90 m，每隔0.5 m记录一道数据.

图 2

Fig. 2

图 2(a)沙漠地区空洞三维模型；(b)相应的二维模型，色标为相对介电常数，无量纲 Fig. 2(a)Three-dimensional model of cavity embedded in the s and ;(b)Corresponding two-dimensional model

整个模型空间离散为0.2 m的立方体单元格，选取中心频率为50 MHz的雷克子波作为激励源，时间采样间隔为0.2 ns，时窗为250 ns.当激励源位于合成记录位置的中点，采用x方向的电场极化，截取不同时刻的波场快照，见图 3，白色星号表示激励源所在位置，白色曲线代表起伏的沙漠地表，两个圆形分别代表地下目标体.观察y=15 m截面不同时刻x方向的电场，图 3a是60 ns时刻的波场快照，可以看到此时波前刚刚传到起伏的沙漠表面，在高程大于平均高程(15 m)的60 m处，已经有部分能量散射到空气中；图 3b是100 ns时刻的波场快照，从图可知，由于起伏地表的影响，一部分能量散射回空气中，并且散射波杂乱交错，剩下的能量穿透沙漠地表继续往前传播，其波前面已经触及到浅部异常体，并且波前面变的不规整；图 3c是140 ns时刻的波场快照，此时波前面已经传到深部异常体，发现两个异常体的反射波能量弱于空气中的散射杂波能量；图 3d是180 ns时刻的波场快照，其波前面可以看做是在起伏地表上的许多二次源的波前面叠加效果，可见起伏地表会造成地表散射杂波，而且其能量往往强于地下目标体反射回波信号.

图 3

Fig. 3

图 3 不同时刻的波场快照 (a)60 ns;(b)100 ns;(c)140 ns;(d)180 ns. Fig. 3 Wave field snapshots at different time

机载雷达探测实际过程中，由于直升机的飞行速度、当地风向及风速的变化等影响，天线的姿态会发生变化，如左右旋转，上下摇摆.根据机载雷达的特殊性，本文研究了四种典型的天线配置方式，如图 4所示，天线的飞行轨迹为沿x方向的带右端箭头的虚线，其中FBy代表天线平行并分别配置于直升机的头部和尾部，天线的极化方向为y方向；FBx代表天线同轴共线配置于直升机的头部和尾部，天线的极化方向为x方向；LRy代表天线共线配置于直升机的左右两翼，天线的极化方向为y方向；LRx代表天线平行配置于直升机的左右两翼，天线的极化方向为x方向.

图 4

Fig. 4

图 4 机载探地雷达四种不同的天线配置方式 Fig. 4 Four different antenna configurations for airborne GPR survey

根据四种典型的天线配置方式，研究了不同极化方式对机载探地雷达探测地下异常体的影响.选取粗糙度为kσ₂=1的高斯随机表面代表起伏的沙漠表面，固定天线距离沙漠地表高度为12 m，模型参数与图 2所示三维模型保持一致.图 5是四种典型的天线配置情况下获得的机载探地雷达探测地下 空洞掩体的雷达剖面，图中的直达波均已去除.图 5a是天线分别配置于直升机头部和尾部，y方向极化时获得的模拟雷达剖面，从图可知，地表的散射杂波主要集中在110~130 ns，其能量远远强于地下目标体能量，浅部异常体反射信号完全淹没于杂波中，但是来自深部异常体的双曲形态反射信号仍然能被识别；图 5b是天线同轴共线分别位于直升机头部和尾部，x方向极化时获得的模拟雷达剖面，从中可知，地表的散射杂波能量仍然比较强，但是浅部的异常体目标反射信号以及深部的目标反射信号均可以被识别；图 5c是天线分别位于直升机左右两侧，x方向极化时获得的模拟雷达剖面，强地表杂波淹没了浅部目标体信号；图 5d是天线分别位于直升机左右两侧，y方向极化时获得的模拟雷达剖面，图中来自地表的散射杂波能量主要集中在110~130 ns范围内，浅部目标体仍然可以识别.

图 5

Fig. 5

图 5 四种不同天线配置方式下模拟的雷达剖面 (a)FBy;(b)FBx;(c)LRy;(d)LRx. 色标为归一化振幅. Fig. 5 Simulated profile results for four antenna configurations of Airborne GPR

综上可知，当天线的极化方向与异常体的走向平行时，地下浅部目标体的信号被来自地表的强杂波覆盖，地下浅部异常体目标很难被识别，而当天线的极化方向与异常体垂直时，地下异常体目标反射波能量增强，虽然地表杂波能量较强，但是浅部异常体目标仍然可以被识别.因此，在机载探地雷达实际探测过程中，在保证测线与异常体走向垂直的同时，天线飞行过程中应尽量保证天线的极化方向与异常体走向垂直.

机载探地雷达探测过程中，为了保证天线的飞行姿态稳定，天线的飞行轨迹应避开地表植被以及地面建筑物，即天线距离地面高度不能过低.为了研 究天线距离地表高度对探测结果的影响，选取粗糙 度为kσ₂=1的高斯随机粗糙表面代表起伏的沙漠表面，天线位于直升机左右两侧，激励源的极化方向为x方向极化，模型中的其他参数保持不变，改变天线的高度，分别从6 m到8 m、10 m、12 m(如图 6所示)模拟了不同高度情况下沙漠地区地下空洞掩体机载雷达探测雷达剖面，图中直达波均已去除.图 6a是当天线距离地表平均高度为6 m时得到的合成雷达剖面，从图可知，位于浅部的异常体及深部的异常体反射波双曲形态同向轴清晰可见，来自地表的散射杂波能量较强，图 6b是天线距离地表平均高度为8 m时的合成雷达剖面，图 6c是天线距离地表平均高度为10 m时的合成雷达剖面，图 6d是当天线距离地表平均高度为12 m时获得的雷达剖面，从图可知，随着天线距离地表的距离增加，来自地下异常体的反射波的能量越来越弱，其双曲线型同向轴的振幅变小，同向轴越来越模糊，而来自地表的反射杂波的能量仍然很强.对于机载探地雷达探测地下异常体，发射天线发射的电磁波能量在空气中由于几何扩散而衰减，其能量与距离的二次方成反比，随着天线高度增加，到达地表的能量减少，从而降低了信噪比，不利于地下异常体的探测.机载雷达野外实际探测中，为了得到更高分辨率的雷达图像，在保证天线飞行姿态稳定的情况下，天线高度尽可能降低.

图 6

Fig. 6

图 6 不同天线高度下模拟的雷达剖面 (a)6 m;(b)8 m;(c)10 m;(d)12 m.色标为归一化振幅. Fig. 6 Simulated radar profile of different antenna height

传统探地雷达探测过程中，天线与地面耦合，大部分能量直接向地下传播；而机载探地雷达探测过程中，天线位于距离地面较高的空中，发射天线发出的电磁波在空气中由于几何扩撒而损耗能量，当电磁波到达地表时，粗糙不平的地表会导致电磁波的空间各个不同方向的散射，为了研究不同粗糙度表面对探测结果的影响，把天线距离地面的高度固定为12 m，天线位于直升机左右两侧，激励源为x方向极化，其他模型参数保持不变，改变沙漠地表的粗糙度.图 7是沙漠地表不同粗糙度情况下获得的机 载探地雷达探测地下空洞的模拟结果.图中水平方 向是沿x方向的距离，垂直方向是电磁波双程传播时间，图中所有剖面的直达波都已去除.

图 7

Fig. 7

图 7 机载探地雷达于沙漠地区探测地下空洞的模拟结果 (a)表面水平情况；(b)表面粗糙度为kσ1=0.5；(c)表面粗糙度为kσ2=1；(d)表面粗糙度为kσ3=2.色标为归一化振幅. Fig. 7 Radar profile of airborne GPR detecting cavities at desert area (a)With smooth surface;(b)With kσ1=0.5 roughness surface;(c)With kσ2=1 roughness surface;(d)With kσ3=2 roughness surface.

图 7a是沙漠地表水平情况下得到的雷达剖面，从图可知，水平同向轴是来自沙漠地表的反射波，而两条双曲线分别代表来自地下空洞的反射波，两条双曲线的顶点位置反映了异常体所在的水平位置，从剖面可知，第一个异常体水平位置为48 m，第二个异常体的水平位置为68 m，这与模型是一致的.图 7b是沙漠表面粗糙度为kσ₁=0.5时的雷达剖面，相比于图 7a，可以发现，虽然来自沙漠表面以及地下空洞的反射波均能清晰地识别，但是雷达剖面在110~130 ns之间出现了很多地表散射杂波，这些杂波的存在势必给剖面的解释带来困惑，另外来自地下空洞的反射波的振幅稍有变小.图 7c是沙漠表面粗糙度为kσ₂=1时的雷达剖面，仍然可以清晰地识别来自沙漠表面的反射波，然而由于沙漠表面粗糙度的增加，导致电磁波在地面的散射更严重，穿透到地下的电磁波能量更少，所以来自地下空洞的反射同向轴变得模糊，位于110~130 ns的散射波能量变的更强.图 7d是沙漠表面粗糙度为kσ₃=2的雷达剖面，来自地表的散射杂波能量进一步增强，而来自地下异常体的反射波的双曲线形态已经很难被识别.从该剖面中，我们很难定位地下异常体的位置.

由此可见，粗糙的地表会导致电磁波的散射，地表散射而损耗的电磁波能量越多，散射杂波的能量越强，同时，穿透到地下的电磁波能量相应的就越少，来自地下目标体的反射波能量就越弱.反映在雷达剖面中图像分辨率降低，不利于机载探地雷达探测的开展.

雷达电磁波与地震波在运动学上存在相似性，反射地震上的偏移技术可以借鉴到探地雷达数据处理中.偏移的目的是为了使雷达剖面中倾斜界面共深度映像聚焦，使绕射波归位.为了消除机载探地 雷达探测中由于起伏地形产生的地表散射杂波影响，借鉴反射地震中的逆时偏移成像(RTM)技术.

地震上逆时偏移可分为叠前逆时偏移和叠后逆时偏移.叠前偏移是把共炮点道集记录或共偏移距道集记录中的反射波归位到产生它们的反射界面上并使绕射波收敛到产生它的绕射点上，最后得到能够反映界面反射系数特点的并正确归位了的地震波形剖面.叠后偏移是在水平叠加剖面的基础上进行的，针对水平叠加剖面上存在的倾斜反射层不能正确地归位和绕射波不能完全收敛的问题，采用了爆炸反射面的概念来实现倾斜反射层的正确归位和绕射波的完全收敛.相对其他方法而言，逆时偏移用全程波方程对波场延拓，避免对波动方程的近似，因此没有倾角限制，原理上可以利用转换波、棱镜波、或者多次反射波成像，并获得更精确的振幅等动力学信息，实现保幅成像，还可以更好地对复杂速度成像进行更细化更精确的估计，成像方法不受介质速度变化的影响，能够对复杂区域进行较准确的成像.

本文借鉴反射地震叠前偏移技术，采用和勘探地震中类似的共炮集观测方式合成雷达数据.为了合成共炮集数据，采用1个发射天线及8个接收天线，接收天线间距为1 m，发射天线与第一个接收天线的偏移距为3 m，发射天线每次移动距离为1 m，这样可以获得一个共炮集数据集.合成数据时，天线距离地面的高度固定为12 m，激励源为x方向极化，其他模型参数保持不变，改变沙漠地表的粗糙度，获得不同粗糙度情况下共炮集雷达数据集.

采用二维RTM算法，对共炮集雷达数据集进行逆时偏移成像，同时为了对比说明RTM的成像效果，本文给出了基尔霍夫偏移成像(Duquet et al., 2000；Eaton et al., 2002；Zhu and Lines， 1998)结果，如图 8所示.图中左边一栏四个图是来自基尔霍夫偏移成像的结果，右边一栏是来自二维逆时偏移成像的结果.图 8a是水平地表情况下基尔霍夫偏移成像结果，其地下目标体的反射波聚焦较好，聚焦位置基本反映了地下目标体的空间位置，图 8e是地表情况下逆时偏移成像结果，从中可知，地下目标体的空间位置、几何形态及尺寸与实际模型基本一致；图 8b是在kσ₁=0.5的粗糙表面情况下获得的基尔霍夫偏移成像结果，图中散射杂波仍然杂乱无章，散射杂波并没有得到聚焦，杂波的存在仍然会影响地下异常体的判断，图 8f是对应的逆时偏移结果，来自地表的散射杂波归位，地下目标体清晰可见，地下目标体的空间位置、几何形态、尺寸以及起伏地表与实际模型基本一致；同样的，图 8g和图 8f中的偏移成像结果明显好于图 8c和图 8d的.

图 8

Fig. 8

图 8 不同地表粗糙度情况下雷达剖面偏移成像结果 (a)地表水平时基尔霍夫偏移结果；(b)kσ1=0.5时基尔霍夫偏移结果；(c)kσ2=1时基尔霍夫偏移结果；(d)kσ3=2时基尔霍夫偏移结果； (e)地表水平时逆时偏移结果；(f)kσ1=0.5时逆时偏移结果；(g)kσ2=1时逆时偏移结果；(h)kσ3=2时逆时偏移结果. Fig. 8 GPR migration results for different roughness surface (a)Kirchhoff migration result at smooth surface；(b)Kirchhoff migration result at kσ1=0.5；(c)Kirchhoff migration result at kσ2=1；(d)Kirchhoff migration result at kσ3=2；(e)Reverse time migration result at smooth surface；(f)Reverse time migration result at kσ1=0.5；(g)Reverse time migration result at kσ2=1；(h)Reverse time migration result at kσ3=2.

基尔霍夫偏移技术只能让地下反射目标体信号聚焦，聚焦位置大致反映目标体空间位置，无法聚焦的散射杂波仍然会对数据解释造成困扰，随着地表粗糙度增加，地下目标体聚焦能量减弱，而逆时偏移 成像技术不但可以使地下目标体反射波聚焦，刻画 地下目标体空间位置、几何形态及尺寸大小，而且可以让地表散射杂波归位，同时起伏的地形亦得到刻画.

(1)数值模拟试验表明天线的极化方向与异常体走向垂直时更利于探测地下目标体，因此机载雷达野外探测时，为了得到更高分辨率的雷达图像，测线与异常体走向垂直的同时，尽可能地保证天线的极化方向与异常体走向垂直.

(2)随着天线高度的增加，从发射天线发射出的电磁波能量由于几何扩散而衰减更多，不利于探测，实际野外探测中，保证天线姿态稳定的情况下，飞行高度尽可能地降低.

(3)地表粗糙度的增加会导致地表散射杂波增强，造成传播到地下的电磁波能量减少，削弱地下目标体信号，雷达剖面的分辨率降低，不利于地下异常体目标识别.

(4)成功借鉴反射地震上逆时偏移成像技术，使地下目标体反射波聚焦，聚焦位置很好地刻画了异常体的空间位置、几何形态及尺寸大小，同时，可以使地表散射杂波归位，起伏的地形在偏移后的剖面中得到刻画；随着地表粗糙度的增加，基尔霍夫偏移成像结果更差，而逆时偏移成像结果仍然能够很好地刻画起伏地形及异常体.