探地雷达数据采集与解释

探地雷达数据采集以及解释
山东大学岩土中心
第1章．探地雷达简介
1．1工作基本原理
探地雷达（Ground Penetrating Radar，简称GPR）是利用频率介于106～109Hz的无线电波来确定地下介质的一种地球物理探测仪器。

随着微电子技术和信号处理技术的不断发展，探地雷达技术被广泛应用于工程地质勘察、建筑结构调查、公路工程质量检测、地下管线探测等众多领域。

探地雷达的基本原理如图1所示。

发射天线将高频短脉冲电磁波定向送入地下，电磁波在传播过程中遇到存在电性差异的地层或目标体就会发生反射和透射，接收天线收到反射波信号并将其数字化，然后由电脑以反射波波形的形式记录下来。

对所采集的数据进行相应的处理后，可根据反射波的旅行时间、幅度和波形，判断地下目标体的空间位置、结构及其分布。

探地雷达是在对反射波形特性分析的基础上来判断地下目标体的，所以其探测效果主要取决于地下目标体与周围介质的电性差异、电磁波的衰减程度、目标体的埋深以及外部干扰的强弱等。

其中，目标体与介质间的电性差异越大，二者的界面就越清晰，表现在雷达剖面图上就是同相轴不连续。

可以说，目标体与周围介质之间的电性差异是探地雷达探测的基本条件。

图1 探地雷达基本原理
1．2电磁波传播特征
探地雷达的电磁脉冲在介质中的传播速度为：
v
=
其中c为电磁波在空气中的传播速度，ε为介质的介电常数，常见介质的介电常数如表1所示。

材质相对介电常数材质相对介电常数
粉质粘土 6 水81
干砂3～5 灰岩4～8
湿砂20～30 花岗岩4～7
金属300 砂岩 6
PVC塑料 3.3 页岩5～15
混凝土 6.4 淤泥5～30
空气 1 海水80
粘土5～40
表1各种常见介质的介电常数
电磁波脉冲在地质界面上的反射系数为：
k=
根据电磁脉冲的传播规律，在地质界面上如果反射系数为负，则相位与发射脉冲相反，若反射系数为正，则相位与反射脉冲一致。

如图2和图3，可以清除看到反射波相位的变化规律。

图2
图3
1．3雷达的分辨率
对于地质雷达的探测方式．它的分辨率也是一个必须了解的内容．地质雷达的分辨率包括垂直分辨率和水平分辨率。

地质雷达的垂直分辨率主要由地质雷达的波长的二分之一决定。

从波的传播规律可知，识别目标体的尺度一般需大于1／2波长，假设垂育最小可分辨的层的厚度为Dm。

则它的计算式为：
其中，c为电磁波在真空中的传播速度。

可见频率越高，介质的介电常数越大，D m越小．即垂直可分辨层的厚度越薄，垂直分辨率越高。

地质雷达的水平分辨率是指地质雷达在水平方向上所能分辨的最小异常体的尺寸的能力。

通常用Fresenel来表示，当反射面的埋深为H，发射，接受天线的间距远小于H时，
R F就是水平分辨率的最小尺度。

从计算公式可以看出．当目标体埋深越大，雷达波频率越低，波长越长，则R F越大，水平分辨率越低．反之，水平分辨率越高。

第2章．雷达数据的采集
雷达数据采集的步骤如下：
1．仪器安装调试；
2．现场地质以及其它情况记录；
3．布置测线；
4．参数设置；
5．参数现场校核，如不合格，重新调试参数；
6．数据采集。

下面将各步骤的注意事项说明。

2．1仪器安装调试
①必须在断电的状态下进行安装和拆卸；
②如果现场空气比较潮湿或者有水，注意防潮防水，以免短路导致损坏电路；
2．2现场地质记录
记录现场的地质条件，以免现场的干扰在雷达记录中造成假象，扰乱最终的解释。

2．3布置测线
如果已知探测对象的大体走向，尽量使测线与其走向正交；
2．4参数设置
参数设置是关系到采集数据质量的关键工作，各项参数如下：
4）：
注：①增益：地质雷达发射的电磁波在介质中传播过程中，在电性（介电常数）分界
面上会发生发射，有一部分电磁波继续向下传播，传播过程中电磁波能量会被介质吸收。

随着深度的增加，电磁波能量减弱，信号幅度相应地减小，不利于信号识别和辨认。

为了能更好地识别信号特征，采用增益（gain）函数来提高信号的幅度，使得信号的细微变化更容易显示和识别。

增益菜单有两个选项：手动/自动Manual/Auto，增益点数Points,GP1、GP2、GP3…GP5。

在采集主机的屏幕的右半部分有一个示波器窗口，除了显示波形之外，还有一条红色的曲线，该曲线就是增益曲线，曲线的转折点点就是增益控制点。

适当的增益函数会提高信号的可视性，但是增益过大会出现削波现象，应该避免削波现象。

②位置确定：雷达剖面的最顶端并不是第一个反射面（比如在地面探测时，剖面的顶端并不是地面），这主要是因为：第一，系统延时，即主机给出发射指令到天线开始发射的延迟时间，第二，直达波，即有发射天线直接到接受天线的电磁波。

为了更加精确的定位，应该去除这两方面的干扰。

但是由于确定零线较为困难。

可以在探测时在地面放置一根电缆与测线正交，天线经过电缆时在剖面上会记录下电缆的位置，通过识别电缆就可以确定零线的位置了。

图4 collect菜单的注解
Playbacck菜单（如图5）
该菜单用于浏览已经采集的文件，在该菜单下有一个SCAN和PROCESS子菜单，通过子菜单可以可以进行一些参数修改和数据处理，以便得到更好的显示效果。

注：这些修改和处理只是用于显示，并不会改变原数据。

图5 Playbacck菜单的注解
Output菜单（见图6）
该菜单的主要功能就是设置数据的显示参数和数据的导出。

图6 Output菜单的注解
System菜单参数设置（见图7）
图7 System菜单的注解
2．5参数现场调试
调试技巧：
①如果出现削波（振幅过大）或者振幅过小的现象，需要调节增益设置；
②如果出现雪花现象（特别是深部），需要增加叠加次数；
③如果存在已知深度的目标体并且该目标体可以在雷达剖面上识别的话，可以通过该目标
体的深度来反算波速，进而求出介电常数。

这样求出来的介电常数比较接近真实值。

如果没有已知深度的目标体，可以打钻确定一个目标体，然后量测其深度，用同样的方法来反算介电常数。

2．6 数据采集
注意事项：
①天线拖动过程中要匀速，并且不要发生跳动现象，跳动现象会造成非常大的干扰。

②打标记要及时，尽量不要重复。

③隧道拱顶或者拱腰部位探测时，可以使用台车（由装载机拖动台车），汽车架子和装载
机架子，如图8。

其中台车最安全，汽车架子最不平稳。

台车（由装载机拖动）汽车架子装载机架子
图8 各种架子
第3章数据分析
数据分析和数据处理的目的就是压制干扰，突出有效信号，提高信噪比，这是进行成果解释的前提，只有进行仔细的处理，才能获得良好的效果。

图9为数据分析和处理的流程图，其中数据分析和数据处理需要进行多次调试对比，才能达到较好的处理效果。

本章将对各个步骤进行详细讲解。

图9 雷达数据分析和处理流程
3．1 绘制测线布置图和尺寸调整
3．1．1 绘制测线布置图
内业工作的第一步就是将现场的测线布置草图绘制成正式文档。

由于现场雷达数据文件是以编号命名的，所以每条测线所对应的文件编号需要记录在文档中，以防混乱。

如下图10，就是雷达探测桥洞的测线布置图，可见在图中标注了测线的方向以及测线所对应的文件编号，应该尽量多的将现场信息反映在布置图上，有助于在解释过程中识别干扰，必要时可以做文字说明。

图10 雷达测线布置图
3．1．2雷达数据剖面的尺寸调整
由于在数据采集过程中存在以下的情况：
①采用时间连续采集模式的时候天线走速不均匀，导致标记（等距离标记）之间的道数不
一样，甚至差别非常大；
②地表起伏比较大，容易在雷达剖上造成假象；
③由于直达波的存在，使得剖面的最顶部并不是地面的反映，使得深度产生误差。

针对以上的情况，需要进行尺寸调整，各功能模块如下：
剖
面
尺
寸
调
整
水平尺寸调整
垂直尺寸调整
距离正常化
水平长度的缩放文件拼接
深度正常化
零线调整
距离正常化
距离正常化，该项功能允许你在标记（必须是等距标记）之间建立等长尺寸，也就是要求在标记之间每单位距离上的尺寸等同，或道数一样。

在没有测量轮的情况下以连续模式采集数时，天线移动速度难以保持恒速，这样就需要利用距离正常化的功能。

联合标记或者距离标记（注：关于Mark数据库的性质与使用见软件说明书的的第二章）设定后，该项功能就会通过增或删的方式来修正每个标记之间的道数，当然同时也修正了采集速度（天线运行速度）。

参数输入框
注意：
①在运行该功能之前，必须确保标记信息是正确的（无重复标记，无丢失标记，首尾标记都存在，所有的用户标记都已经转化为距离标记或联合标记）。

②Scans/unit,unit/mark都必须在头文件中设置好，以便运行该功能。

③unit/mark是根据测量时的设定来设置的。

④Scans/unit需要用鼠标来清点每个标记之间的道数（注意：因为在雷达剖面中显示的是已经做了叠加处理的数据，若每个标志之间的道数为3，而叠加次数为16，则实际每个标志之间的道数为二者的乘积48，切记！！）
⑤该功能运行之后，原数据中在第一个标志之前的部分已经被cut off了。

水平尺寸缩放
可以通过水平尺寸缩放中的叠加，去除，添加功能来修改雷达数据。

参数设置对话框
注意每次只能这三项功能中的一个。

叠加（stacking）：使用该功能可以对数据进行简单的滑动平均处理。

该功能就是将所指定的几道平均叠加之后输出一道数据。

去除（skipping）：选定此项，你可以将指定的道去除，比如，你输入参数1，就会每隔一道去除一道数据（因此，数据被压缩到原来的二分之一）
添加（stretching）：选择此项，将拓展水平尺寸。

该功能将计算出每相邻两道的平均值（或指定道数），然后将平均道添加到已有数据中。

拼接文件
有时因为场地条件所限，我们不得不分部采集数据，在后处理中为了将各部分数据连接在一起，就要使用该功能了。

①选择File＞Append File.
②选择所要连接的各个文件。

③点击Done，完成连接，重新命名加以保存。

表面位置调整（非常重要）
由于系统延时和直达波存在，使得整个剖面的最顶部并不是地面的位置。

确定地面的位置，对于精确的深度定位来说非常关键，但是如何确定，目前还没有定论，下面列举三个常用的方法。

①根据RADAN介绍，百分之九十的情况下把直达波的第一个正峰位置作为地面。

②根据华东院资料，将直达波的第二个波瓣作为地面，如下图
③根据经验，可以在探测时在起始部位放置一根电缆，在后处理时在剖面上识别出该电缆，
这样就可以确定地面位置了。

下面以一个探测剖面为例说明一下这三种方法的用法与区别（如图11）：
图11 零线的确定
如图11，其第一个正峰位置是10.34ns，第二个峰值是14ns，而电缆位置是16.34ns，可见三者最大相差6ns,按照介电常数为8计算，深度偏差了30cm，对于超前预报来说该误差可以允许。

而对于衬砌检测来说偏差较大，需要综合三种方法来分析。

表面正常化
在测线布置时会遇到地表起伏较大的情况，这就需要修正地表起伏对数据剖面的影响，进而可以使水平或接近水平状的反射体的反应更接近实际。

通过输入标记的z值就可以实现该功能。

3．2 数据分析
傅立叶谱分析
数据处理是进行数据解释的基础，目前比较常用的处理方法有一维滤波，二维滤波，以及反滤波，这几种方法都是以傅立叶谱分析为基础的，傅立叶谱分析是将雷达数据由时间域转化为频率域，表现的是各种谐波频率的振幅分布，如图12。

图12 振幅频谱图
关于不同探测介质的振幅频谱特征，一下有几个结论（摘自杨峰资料）：
（1）水对高频电磁波具有很强的吸收作用，这与水离子导电是密切相关的，离子导电增加了介质的电导，而电磁波传播与电导和频率之间呈指数衰减关系。

（2）花岗岩不但对高频成份具有一定吸收，而且形成的振幅谱比较单一。

（3）在干燥的不均匀介质中，形成的振幅谱不但主频特征不明显，而且在天线的高端会形成一定的杂波信号。

这可能是由于高频电磁波在不均匀介质内形成多次干涉造成的。

干涉现象势必加宽信号的频带特征。

具体的谱图如图13。

空气背景雷达频谱特征水背景雷达频谱特征
花岗岩背景雷达频谱特征干燥碎石背景雷达频谱特征
图13 不同介质的频谱图
希尔伯特变换
由于大地介质的不均匀性，地质雷达发射的高频脉冲电磁波在地下传播过程中将发生强烈的衰减、反射、折射、绕射和散射，这些反射波、折射波、绕射波和散射波相互叠加在一起，为数据处理带来了巨大的困难；同时，为了得到更多的反射波特征，地质雷达通常利用宽频带进行记录，因此不可避免地记录下各种干扰噪声。

如果噪声频率带与反射波频率带重叠或接近，利用傅立叶谱分析技术对这样的信号进行分析，有时难以取得理想的效果，严重
影响了图像解释的可信度和精度，进而影响了地质雷达的探测效果。

而希尔伯特变换可以较好的解决这个问题，希尔伯特变换就是将记录道的信息直接在时间域上转化为瞬时振幅，瞬时相位，瞬时频率的技术。

复信号的瞬时振幅、瞬时相位、瞬时频率这3种瞬时信息，一般是指一个特定的瞬间，而不是一个时间段的平均。

地质雷达信号记录道x(t)的复信号分析与地质雷达信号的傅立叶谱分析分别在时间域和频率域上对地质雷达信号的能量、频率和相位等参数进行分析检测，它们在振幅上无本质差别，而瞬时频率与傅立叶分析的频率不同，前者是分析全部谐波叠加波形的视频率，后者则是分析各谐波频率的振幅分布情况。

两者既有区别，又有一定的内在联系。

复信号分析技术与傅立叶谱分析技术的成果输出不同，它可以将地质雷达记录中的瞬时振幅、瞬时相位和瞬时频率分离出来，得到同一个剖面的3个参数图，因而其解释方法与
质雷达信号总能量的平方根，利用这种特征便于确定特殊岩层的变化。

当地层存在明显介质
波的能量强弱，它的相位都能显示出来，即使是弱振幅有效波在瞬时相位图上也能很好地显示出来。

当电磁波在各向同性均匀介质中传播时，其相位是连续的；当电磁波在有异常存在的介质中传播时，其相位将在异常位置发生显著变化，在剖面图中明显不连续。

因此利用瞬
岩性变化，有助于识别地层，当电磁波通过不同介质界面时，电磁波频率将发生明显变化。

这种变化可以在瞬时频率图像剖面中较为清晰地显示出来，在地下介质发生变化的时候，瞬时频率也会发生显著变化，需要指出的是，在反射层处瞬时频率的大小在数值上与反射波的主频对应的很好，所以可以利用瞬时频率的大小和稳定情况来判断地下介质的稳定性和岩性变化。

对于同一探测对象，3种瞬时信息在同一位置发生明显变化就可能反映探测对象在该处的物性变化。

因为在这3个参数中，瞬时相位谱的分辨率最高，而瞬时频率谱和瞬时振幅谱的变化也较为直观，所以通常根据瞬时频率谱和瞬时振幅谱来确定地下异常或分层的大概位置，然后利用瞬时相位谱精确确定异常位置和分层轮廓线。

有些时候，也可以直接利用瞬时相位谱来确定地下异常的位置。

具体的分析实例见第五章。

第4章数据处理
数据处理是进行数据解释的基础，在RADAN中数据处理的方法非常多，应该在数据分析的基础上决定采取那些处理方法，采取怎样的处理步骤。

针对不同的目标，有不同处理方法。

各种处理方法的使用方法和注意事项在下面具体讲述。

4．1 去除水平噪音
所谓的水平干扰信号，就是指水平带状干扰，通常具有低频特征，经常会干扰一些真实
的反映体，如下图：
（a）水平干扰（b）频谱图（可以看出低频干扰较多）
图1 水平带状干扰及其对应的频谱图
图2 水平干扰（可以看出剖面主要被水平信号覆盖）
IIR水平高通滤波
滤波器长度应该先设为数据剖面的最大道数（应该为奇数），如果该数值超过225，则
应该选择225。

这样，在水平方向上长度等于或超过该值的特征将被执行滤波，而长度低于
该值的特征受影响很小。

注意：①因为对于IIR水平高通滤波器，其长度最大为255，所以长度超过该值的特征
都将会被执行滤波，这是不可避免的；
②因为直达波也是水平信号，为了不对直达波产生影响，可以通过设置起始/终止样本
点来圈定滤波区域，避开直达波。

FIR背景去噪
滤波器长度应该先设为数据剖面的最大道数（应该为奇数），如果该数值超过1023，则
应该选择1023。

这样，在水平方向上长度等于或超过该值的特征将被执行滤波，而长度低
于该值的特征受影响很小。

FIR背景去噪滤波器，其长度最大为1023，所以长度超过该值的特征都将会被执行滤波，这是不可避免的；
②因为直达波也是水平信号，为了不对直达波产生影响，可以通过设置起始/终止样本点来圈定滤波区域，避开直达波。

因为水平干扰信号往往具有低频干扰，所以垂直高通滤波器可以进行相应的处理，具体滤波器设计需根据频谱特征来确定。

4．2 去除高频干扰
高频信号经常表现为雪花形状，对数据造成了较大干扰。

可以通过垂直低通滤波，水平低通滤波，滑动平均滤波来进行处理。

垂直低通滤波分为IIR 和FIR形式，可以根据频谱图来确定具体的滤波参数。

其原理是当你输入一个非零值，由该值决定的道数会相加平均并将平均值赋予中间道，依次计算。

所以参数值应该为奇数，一般情况设为5就可以很好的去除高频，平滑数据。

4．3 空间滤波
以上就是F-K滤波的流程图，首先经过快速傅立叶变换将时间域雷达数据转化为二维谱图，经过对二维谱图的分析，选择合适的参数进行快速傅立叶逆变换，得到处理之后的时间域雷达数据剖面。

空间快速傅立叶变换滤波器，是一个二维的频率滤波器，在时空二维域中进行滤波。

经常被称作频率－波数滤波，或f-k滤波（注：k就是波数的意思）。

这种方法可以产生一个二维矩阵，代表了雷达波的相位和振幅。

可以用此滤波器进行二维滤波以削减噪音干扰。

对已经变化的数据矩阵进行傅立叶逆变换，此时的滤波器会滤掉一些噪音。

在技术上，通过逆变换，数据由频率域恢复到时间域。

相对于一维的垂直和水平滤波，F-K滤波的优点有：可以对信号和噪音进行更好的区分。

信号和噪音或许在一维处理中会有所重叠，使得分离它们变得非常困难。

但是二维滤波中的情况好的多。

该功能的对话框如下图：主要显示了二维谱图，滤波器参数设置，以及谱图的显示控制参数。

其中谱图的竖轴代表了信号频率，横轴代表了波数（即每单位长度上波周的数目）。

Recalc按钮会将文件转化为二维谱图。

谱图的显示参数可以控制谱图的显示质量，一旦谱图形成，可以用Gain，Zoom来增强显示效果，可以使用Scans，Samples来选择显示范围（显示范围也可以通过鼠标来控制）
设置滤波器参数可以开始快速傅立叶逆变换。

滤波器参数的意义和选择
·Min Freq最小频率，Max Freq 最高频率，这两项控制着滤波器的竖向分量。

·Alpha 和 Delta Alpha控制着水平分量，可以通过图中的直线在调整。

Alpha 代表着滤波器的对称程度，当两条直线关于中间直线对称时，Alpha的值接近于0，当Alpha 很高时，意味着两条射线不对称。

Delta Alpha 代表着射线之间的夹角，与反射体的线性尺寸相关。

滤波器类型的选择：该项决定了使用哪种空间快速傅立叶逆变换，总共有五种：
None: FFT文件不作任何修改被恢复，
High-Cut Horizontal: 仅两射线之间的部分被执行FFT逆变换。

High-Cut Vertical: 仅两射线之外部分被执行FFT逆变换。

High-Cut Vert Symm: 仅两射线之外部分被执行FFT逆变换（对称的）
High-Cut Horz Symm: 仅两射线之间的部分被执行FFT逆变换（对称的）
注：当射线不对称时，使用对称滤波器形式可以起到较好的作用。

使用经验：
①不对称的射线与对称的滤波器类型组合，往往得到比较好的效果。

②射线对称时（即Alpha的值接近于0），水平和竖直特征会被突出。

不对称的射线（即Alpha 的值很高）突出倾斜的特征。

④对于High-Cut Horizontal和High-Cut Horz Symm来说，Delta Alph较小时且对称时，突
出的是水平信号，当Delta较大且对称时，则可以包含各种信号。

如下图。

原始数据Alpha较小Alpha较大
⑤High-Cut Vertical和High-Cut Vert Symm方式，当Delta Alph较小且射线对称时，包含各种信号，当Delta较大且对称时，突出竖直信号。

所以运用这个功能可以去除水平干扰。

原始数据Alpha较小Alpha较大
⑥当射线不对称时，突出的是倾斜信号，如下图，可见High-Cut Horizontal倾斜信号不突出。

High-Cut Vertical模式High-Cut Vert Symm High-Cut Horizontal
4．4 去除多次反射－反卷积
当雷达信号在目标体（如一块金属物或湿粘土层）和天线之间来回反射的时候，往往会出现多次反射界面的现象。

或者在两个反射层面之间发生电磁波的振荡，出现多次反射现象。

这种现象会模糊浅部的（或者深度较小的部分）真实信息。

在实际探测中在扫描地下水层，基岩或空洞的时候就会出现多次重复的情况。

反卷积就是为去除此类噪声干扰而设计的滤波方法，还可以提高垂直分辨率，分解间距较小的层。

RADAN中的反卷积方法叫做预测反卷积，这是一种将尖脉冲反卷积作为一种特例的常见方法。

该方法试图尽力在天线与地面耦合的时候去逼近发射脉冲的形状。

假设一个特定长度的震源子波，也称作滤波器算子长度，当震源子波由数据中被清除时，该滤波器可以预测一定距离之外的数据形状，叫做预测延迟。

这就导致了反射子波被压制。

像天线重复反射等此类预测现象，将被移动到比预测延迟更远的位置，可以有效的消除此类现象。

反卷积参数选择
为了更好的运行反卷积滤波，像滤波算子长度，预测延迟，预白噪声化，增益，起始样本，终止样本等参数应该适当的选择。

滤波算子长度：按照组成一个脉冲的样本点数目，滤波器算子长度设定了滤波器的大小。

·较长的滤波算子长度可以对雷达波进行较好的拟合，并且可以得到较好的结果，但是耗时较长。

·滤波器算子长度应该满足一个完整的雷达子波循环，这样将起到较好的作用。

小于该值的参数会导致不好的结果。

射峰值，将两者的样本点数相减，而算子长度值应该大于或等于该值。

预测延迟：该值将被设定为理想的输出脉冲长度（大约为雷达子波的半个循环）。

小于该值的参数会产生更多噪声。

·利用反卷积来去除重复反射时，延迟值应该等于或小于重复之间的空间。

·参数值为5－1的预测延迟被用作拟合尖脉冲反卷积，但是这将给数据带来更多噪声干扰。

预白噪声化：通过加强白噪声（零延迟）元件，预白噪声化可以调整相应的自相关函数。

从数学意义上来讲，预白噪声化可以是滤波器稳定，并且可以是输出的数据光滑，降低噪声干扰。

0.1～1是普通值，而0.8是一个较好的值。