低频地质雷达屏蔽天线开发及在地下空间工程探测应用

许献磊，张朝阳，彭苏萍，陈令洲 ，张 迪

(1.中国矿业大学(北京) 煤炭资源与安全开采国家重点实验室，北京 100083；2.中国矿业大学(北京) 地球科学与测绘工程学院，北京 100083)

地质雷达(GPR)技术因其探测频率高、精度高，自带发射激励源、现场操作便捷，在地下工程地球物理勘探方法之中具有明显的技术优势[1]，目前已经广泛应用于矿井超前地质构造探测和城市地下空间安全探测中。常规地质雷达探测天线系统受限于复杂地质条件，存在探测深度浅、深部探测精度降低，且不能实现实时处理功能，在实际应用中与矿井智能化开采地质保障及城市地下空间工程探测需求差距较大。

近年来，低频地质雷达屏蔽天线系统的开发一直是行业学术研究的热点和难点之一。在地质雷达天线研发方面，WU Bingheng[2]、郭士礼[3]、吴彦奇[4-5]等在地质雷达的天线参数设计、屏蔽腔内吸波材料填充等方面设计仿真并应用，对于改善探测效果有一定成效，但其大多是针对高频雷达天线。目前，国内外商业化低频地质雷达天线的中心频率集中在100 MHz左右，如美国劳雷GSSI系列，其探测深度一般为5 m。中国矿业大学(北京)研究团队在2015年开发了低频组合雷达天线，中心频率涵盖了12.5、25、50 MHz，均为非屏蔽天线，已经应用在地质勘查和矿井灾害源探测中。低频天线一般尺寸较大，屏蔽设计难度大。此外，一些学者提出加大发射功率来提高探测深度，但这种做法会增强干扰信号[6-7]。当前地质雷达天线技术面临的突出问题有:①探测深度浅；②受城市环境干扰、深部复杂地质岩体结构、电磁波衰减的影响，深部弱反射信号被噪声掩盖，信噪比低，深部地质异常体的探测精度低。

在地质雷达数据处理研究中，目前存在多种方法：基于小波分析法，OUADFEUL S[8]、ZHANG Z[9]等对 GPR 信号进行了降噪处理；基于杂波模型的方法，如 LP法[10]、均值法；基于信号统计特性的算法，如卡尔曼滤波法[11-12]；基于子空间投影的方法如奇异值分解、带通滤波法[13-14]。这些算法在地质雷达上的应用大多基于后处理软件，处理效率较低，不具备实时性，对于大量级的雷达数据尚不具备快速处理的能力，且处理过程严重依赖操作者的经验，智能水平低[15-16]。

随着芯片技术的不断发展，现场可编程门阵列(FPGA)等一些可编程逻辑器件的并行结构特点为高速图像处理提供了新的思路。基于SOPC(System on Programmable Chip，可编程片上系统)技术将并行结构一体化，为地质雷达数据的实时处理提供了可能。图像处理领域的很多先进理论和方法是基于软件方法执行顺序指令，如VC与MATLAB。但随着处理数据量的增大，以及实时性的要求，具有可重构性、并行性、高速性的FPGA发挥出了其独特的优势，可以很好地满足数据处理需求，例如地质雷达信号接收机中量大且复杂的信号处理[17-19]。

笔者及其团队在前期研究的基础上，开发了一种低频地质雷达屏蔽天线系统：通过低频屏蔽及阻抗匹配设计，增强天线的聚焦和辐射特性，进而提高数据的信噪比；同时提出FPGA内嵌NiosⅡ软核处理器的实时图像处理系统设计方案，实现数据的实时处理。研究成果可提高地质雷达的探测深度，为地质结构精细勘查和安全隐患探测提供技术保障[20-22]。

1 超深屏蔽地质雷达天线设计与开发

1.1 地质雷达天线功能

地质雷达天线主要由发射天线和接收天线组成[23]。对地质雷达而言，天线的辐射场通常是向地下耦合传播。发射天线将激励信号转换为一个在已知时空分布的电磁波波场，当电磁波信号在传播过程中遇到波阻抗面(电性不同的介质)时，会在波阻抗界面处发生反射、折射和透射。反射的电磁波信号被同步于发射天线的接收天线接收，可测出矢量电磁场随时间的变化，并将其转化为一个可记录的电信号。接收天线接收到的反射信号强度取决于不同层介质间的电性差异，差异越大，反射信号越强。通过雷达主机精确记录反射回来的电磁波的运动特征，包括接收到波的旅行时间(亦称双程走时)、振幅信息与相位信息等，获得目标体介质的断面扫描图像，通过对扫描图像的进一步处理和解译，可推断电性界面的形态、埋深和构造等特征。

天线辐射的电磁波在有耗介质(如岩石、矿物质等)中的传播因媒介吸收而迅速衰减，限制了地质雷达的探测深度。因此，低频雷达天线设计在满足大带宽、低旁瓣、低互耦电平等要求的同时还要考虑电介质材料与天线之间的激发场对天线辐射的影响[24]。地质雷达系统所采用的超宽带窄脉冲信号要求天线具有良好的辐射特性、宽带特性及优越时域特性。电磁波的穿透能力与地下介质(如岩石、矿物质等)的电学性质、天线的中心频率密切相关：地下介质的主要电学性质为电导率σ和相对介电常数εr，电导率越低，穿透能力越强；开发低压雪崩高频信号采样电路，中心频率越低，穿透能力越强[25]。

1.2 超深屏蔽地质雷达天线系统设计思路

天线探测深度与天线工作频率及地下介质特性有关。在地下介质特性一定的情况下，天线工作频率越小，电磁波在地下介质中的衰减越小，探测深度越深。而天线的工作中心频率与天线的长度相关。一般而言，天线的长度要大于等于所设计工作频率对应电磁波波长的1/4，这给缩小天线尺寸带来了一定的难度。蝶形天线具有体积小、平面结构、宽频带、方向性好等优点，因此采用蝶形样式来设计超深屏蔽地质雷达天线，天线中心频率为50 MHz。

基于时域有限差分法FDTD(Finite-Difference Time-Domain)的三维电磁仿真软件CST进行天线及其屏蔽结构的仿真模拟：

1)通过天线辐射面几何参数的优化和阻抗匹配设计，实现发射天线电磁脉冲低频率、超宽带的性能要求；针对城市环境干扰，设计天线屏蔽装置，进而增强能量辐射的方向性、提高天线辐射效率，屏蔽外界干扰，提高目标信息信噪比。

2)采用宽频带发射机和高精度接收机，开发级联雪崩电路和地纹波保护电路，提高发射脉冲信号的带宽。

3)开发低压雪崩高频信号采样电路，提高接收机灵敏度，同时减少纹波噪声干扰，实现弱信号的识别。

1.3 天线及其屏蔽结构仿真分析

1.3.1 天线辐射面尺寸及无屏蔽模拟

天线辐射面形状如图1(a)所示，天线辐射面长1 700 mm，宽499 mm。CST软件仿真中辐射天线面基板材质均设置为FR-4(loss free)，厚1.8 mm，蝶形辐射片材质设置为PEC，厚0.2 mm。

(a)天线辐射面结构

从无屏蔽的发射天线和接收天线的回波损耗S11图像(图1(b)、(c))中可以看出，在没有任何屏蔽的情况下，接收天线和发射天线在0～100 MHz的带宽中回波损耗严重，均超过-10 dB。

1.3.2 屏蔽设计仿真分析

1)屏蔽装置背腔的高度H仿真及参数优化

将接收天线和发射天线根据地质雷达系统中设计的距离进行组合[26]，在此基础上设计屏蔽装置的尺寸。测试接收天线和发射天线在屏蔽壳不同高度，以及2层屏蔽壳之间电阻加载不同位置的性能。将屏蔽装置背腔的高度H作为变量，考虑到尺寸大小因素，H分别取450、500、600 mm进行仿真分析，结果见图2。

(a)发射天线回波损耗S11

由图2(a)、(b)可以看出，发射天线在屏蔽装置背腔高度H=450、500、600 mm时，在0～200 MHz仿真范围内均满足S11<-10 dB，驻波比VSMR<2,在屏蔽下发射天线性能满足需求。从S11相关图像中也可以看出屏蔽壳高度H的变化对发射天线回波损耗的影响并不明显。接收天线屏蔽装置背腔高度H=450、500、600 mm时，S11差异明显，在(50±25) MHz内S11均在-10 dB左右(见图2(c))，经过测量只有H=600 mm时，接收天线的回波损耗在0～200 MHz仿真范围内均小于-10 dB，驻波比VSMR<2(见图2(d))。优选H=600 mm作为屏蔽壳的高度。

2)阻抗匹配仿真及参数优化

收发天线和屏蔽壳之间加载电阻的位置对各自的S11曲线的影响十分明显，加载电阻的位置分别为：①不加载电阻；②收发天线两侧加载电阻；③收发天线顶部加载电阻。在确定屏蔽壳高度为H=600 mm的情况下，以加载电阻位置作为变量对收发天线进行仿真分析，结果见图3。

(a)发射天线回波损耗S11

由图3(a)、(b)可以看出：在发射天线两侧加载电阻时，在0～200 MHz仿真范围内均满足S11小于-10 dB，驻波比VSMR小于2;发射天线在不加载电阻和顶部加载电阻的情况下回波损耗在(50±25) MHz带宽内不满足小于-10 dB。由图3(c)、(d)可以看出：接收天线在3种位置下，在0～200 MHz仿真范围内均满足S11<-10 dB，驻波比VSMR<2；接收天线不加载电阻和顶部加载电阻的S11曲线几乎重合，展现出的效果几乎相同，相比于两侧加载电阻的情况下，不加载电阻和顶部加载电阻的S11在51 MHz的位置出现谐振，符合预期设计中心频率要求。

通过对收发天线加载不同高度的屏蔽装置，以及加载不同位置电阻进行的仿真测试，分析天线回波损耗S11和驻波比VSMR，从而确定了屏蔽装置的尺寸和电阻加载位置。在最优屏蔽情况下，即屏蔽壳高度H=600 mm、发射天线两侧加载电阻、接收天线不加载电阻时，仿真区间0～200 MHz内收发天线S11均小于-10 dB，发射天线在50 MHz回波损耗低至-24 dB，接收天线在50 MHz回波损耗低至-17 dB。

天线加载屏蔽与无屏蔽仿真对比如图4所示。从远场天线方向图来看，未加载屏蔽的发射天线和接收天线为全向型天线，辐射范围均匀且有前后2个辐射方向；经过屏蔽装置加载后，可以从远场天线方向图清楚地看到接收天线和发射天线均具有了良好的方向性和增益。未加载屏蔽装置的发射天线主瓣幅值为1.77 dBi，加载屏蔽装置的发射天线在垂直极化方向的主瓣宽度为111.9 deg、主瓣幅值为5.29 dBi、副瓣电平为-10.9 dB；未加载屏蔽装置的接收天线主瓣幅值为1.93 dBi，加载屏蔽装置的接收天线在垂直极化方向的主瓣宽度为103.8 deg、主瓣幅值为5.97 dBi、副瓣电平为-10.7 dB。

(a)发射天线水平极化方向

2 基于FPGA雷达信号实时处理技术

2.1 FPGA实时处理架构

基于FPGA数据处理算法实现过程如下：①在Matlab中对数据进行算法验证，生成(.mdl)文件，而后通过Singal Compiler将(.mdl)文件转成相应的硬件描述语言HDL(.vhd)文件和tcl脚本；②通过ModelSim进行RTL级仿真，验证算法在Matlab上与硬件编译结果的相通性；③通过Quartus Ⅱ 进行门级的时序验证，确保算法的可行性。

本雷达天线系统中选用的FPGA型号为Xilinx的ZYNQ7000，上位机采集控制系统通过USB端口和“FPGA核心控制系统”进行通信。实时处理过程如下：①“FPGA核心控制系统”接收到“上位机采集控制系统”发送的信息，生成步进延迟和固定延迟等信号，通过轮回机制，控制“雷达天线系统”采集；②“雷达天线系统”中接收机采样保持直接送到“FPGA核心控制系”的采集卡[27-30]；③控制单元向发射子系统和接收子系统启动触发脉冲，回波信号经过接收天线接收，变成和差两路信号，经中频放大，送入信号处理机，完成和差两路中频信号的ADC转换，对2种回波信号进行下变频(DDC)处理操作转换到基带，并将得到的数字化雷达反射波信号数据通过FPGA对基带信号进行脉冲压缩;④通过FPGA图像实时处理模块进行实时处理，后经SDRAM存储器进行缓存；⑤将数据通过VGA控制模块以45FPS(帧/s)的帧率实时呈现处理结果，并经微机系统总线存放到内存中[31-32]。雷达信号处理流程如图5所示。

图5 FPGA雷达信号实时处理流程

2.2 数据实时处理算法

地质雷达数据通过FPGA实时图像处理，可以压制随机、规则信号干扰，突出有用信号，以最大可能分辨率在图像上显示反射波，对信号矩阵进行数学算法运算(数据内插、背景去噪、增益、一维滤波和小波变换)，以提升雷达图像的呈现度。下面以北京海淀区某地铁隧道路段为试验探测地点，探究各数据处理方法对地质雷达数据A-Scan与B-Scan图像的影响。

2.2.1 数据内插

地质雷达数据可通过内插的方法增加采样点数，间接地提高数据带宽，在经过不同的数据处理方法后更好地突出异常体目标信号位置，并提高识别空洞、富水等病害体的位置精度。数据内插采用Resample函数功能，其原理是在原始序列的2个采样点之间等间隔地插入整数倍数值，根据奈圭斯特采样定理可知，原始采样序列x(n)可以无失真地恢复出原始模拟信号。恢复出来的模拟信号为：

(1)

式中：y(t)为恢复出来的模拟信号；T为周期，s；x(n)为原始序列。

采样间隔为T′的信号，即将t=mT′代入式(1)即可。以原始数据、4倍插值、8倍插值及10倍插值处理为例，分析其对探测识别精度的影响。不同倍数插值与识别位置如图6所示。

(a)原始数据

原始数据采样点为512，选取雷达图像方框内顶点及其附近5个点作为层位识别点。通过层位识别可以得出打标点处方框内顶点附近5个点的深度，与实际深度进行对比分析。层位深度计算结果如图7所示。

图7 层位深度计算结果

由图7可以看出，经过插值处理后，图像更加平滑，样点细节变化更加丰富，有利于信号的识别与追踪。其中应用8倍插值数据得出的层位识别结果最接近实际深度，效果最优，平均识别误差为2.91%，相比于应用原始数据得出的与实际深度存在18.73%的识别误差，其在定位精度上有了较大提升。

2.2.2 背景去噪

大窗口内时间样点的Hibert瞬时振幅和：

Asum(ti)=∑k∈NAk(ti)

(2)

第n道各点的权重：

(3)

大窗口内权重平均背景道:

(4)

式中:Asum(ti)为瞬时振幅和；Ak(ti)为大窗口内第k道、第i时间样点的瞬时振幅；k为信道号；wn(ti)为各点权重；An(ti)为Hibert瞬时振幅；xwmean(ti)为平均背景道；N为大窗口总道数；xn(ti)为平均处理的振幅。

背景噪声的去除能够显示出噪声之下的更多有效回波信息并且可以从背景去噪后的B-Scan中明显看出有效抑制了直达波信号，如图8所示。

图8 背景去噪处理后的雷达图像

2.2.3 一维滤波

一维滤波可以压制有效信号频带以外的干扰信号，以达到突出有效信号、提高剖面信噪比的目的。如果单位脉冲响应是一个有限长序列，这种系统称为“有限长单位脉冲响应系统”，即FIR滤波器。FIR数字滤波主要根据有效信号频带设置滤波器，常见滤波方式有低通、高通、带通和陷通。根据不同滤波参数可以得到相应的滤波因子和频响的振幅谱。FIR滤波器的设计可以采用窗口法、频率采样法，进行最优化设计。FIR数字滤波器系统函数为：

使用AMOS进行结构方程建模及分析，以找出船员疲劳与疲劳影响因子之间的联系，船员的负荷矩阵见表5，基本研究方法见图2。

(5)

将式(5)分解为二阶实系数因子形式，得到FIR数字滤波器的级联型结构，即：

(6)

经过有限冲击响应系统H(z)后，离散信号x(n)即完成数字滤波处理，变为y(n):

(7)

式中：x(n-k)为离散信号；y(n)为数字滤波处理后信号；h(k)为有限脉冲响应序列；β0、β1k、β2k均为权重参数。

一维FIR滤波处理后的雷达图像如图9所示，可以看出，有效信号频带以外的干扰信号经FFT振幅谱分析后得到去除。

图9 一维FIR滤波处理后的雷达图像

2.2.4 增益

增益是地质雷达图像显示的关键数学运算，常见增益分为时变增益和自动增益。

为了弥补地下反射信号振幅随时间及深度的衰减，使深部的信号经放大后显示清楚，设计了按反射信号在地下介质中传播过程的衰减规律进行补偿，即时变增益。典型的时变增益有扩散和指数补偿SEC。时变增益通用表达式如下：

A′x,t=Ax,t(αt+eβt)

(8)

式中：α为线性增益；β为指数增益。

自动增益即AGC增益。当信号较强时增益自动降低，而当信号较弱时增益自动提高，从而保证强弱信号的均匀性，便于有效波的追踪。自动增益依靠雷达记录乘以随时间变化的因子来实现，算法如下：

x′(t)=P(t)x(t)

(9)

式中:x(t)为自动增益前的雷达记录；P(t)为自动增益权函数；x′(t)为自动增益后雷达记录。

经过增益处理后的雷达图像如图10所示，可以看出，空洞位置与周围振幅由图谱显示明显不同，雷达信号经过有效放大后，图像更易于识别。

图10 增益处理后的雷达图像

2.2.5 小波变换

如果ψ(t)L2(R)，傅里叶变换满足允许条件则称ψ(t)为一个基本小波，将母函数ψ(t)经过伸缩和平移后，就可以得到小波序列：

(10)

任意的函数f(t)∈L2(R)的连续小波变换定义为：

(11)

(12)

式(12)即为小波变换的逆变换[33]。小波变换可以更好地去除杂波，压制噪声，使图像更加精细化，如图11所示。由图11可以看出，经过小波变换后去除了杂波，雷达图像的分辨率得到进一步提高，有用信号进一步得到展现，更易分辨出空洞与地铁隧道的位置。

图11 小波变换处理后的雷达图像

2.3 实时处理效果评价

为了验证基于FPGA数据实时处理与算法的有效性，选取信息熵和平均梯度2个指标对雷达数据FPGA实时处理结果与常规处理结果进行对比评价。

1)信息熵：利用统计学中信息熵的概念可以评价图像信息量丰富度。信息熵越大，所含信息越完备。通过比较实时处理与常规处理的信息熵，评价基于FPGA实时处理的效果。其计算公式如下：

(13)

式中:L表示图像的灰度级数；pi表示频率为i的灰度值。

2)平均梯度：平均梯度是图像边界或影线两侧附近灰度值的变化率，可以作为雷达图像分辨率的一种评判指标，平均梯度越大，图像整体分辨率越高。其计算公式如下：

(14)

式中:M、N分别表示图像的横向和纵向的行、列数；F(i,j)为图像的第i行、第j列的灰度值。

不同处理方式的雷达图像对比如图12所示，信息熵和平均梯度计算结果见图13。

(a)常规处理结果

(a)信息熵

通过以上公式关于信息熵与平均梯度的评价指标计算，结合图13可知，相比于普通数据处理，经过FPGA实时处理得到的雷达图像5条测线数据信息熵平均增加了13.72%，说明通过FPGA实时处理的各种算法，图像信息丰富度提高，使雷达图像解释更具合理性、准确度更高，并且平均梯度增幅为12.44%，说明雷达图像的分辨率得到提高，图像更具可读性，在提高数据处理效率的同时提高了质量。

3 工程探测试验

3.1 试验概况

为验证该天线系统的探测能力及实时处理的有效性，应用设计的50 MHz屏蔽天线系统对北京某地铁施工透水位置附近的地面道路进行探测，如图14所示。

(a)现场数据采集

本次探测试验共计4条测线，采样点数为1 024，叠加次数为10次，时窗根据不同探测位置分别设置。通过FPGA控制采集雷达数据图像并且进行实时处理。

3.2 结果分析

测线1探测结果如图15所示，从中可以发现，测线1含有2处富含水区域(方框区域)，富水区Ⅰ和富水区Ⅱ深度分别为17.2～23.0 m和16.8～22.5 m。根据现场打标位置，富水区Ⅰ的位置与地下突水的位置重合，可以判定该突水是由富水区Ⅰ引起；富水区Ⅱ为预测富水区，在后续的作业中得到验证。本次探测为突水事故救援提供了有效数据支持，对地铁隧道施工安全提供技术保障。

图15 测线1数据处理结果

本次工程探测的地质雷达系统中所用信号处理板卡的FPGA芯片采用200 MHz的系统时钟，通过各个处理流程的时钟延迟可以得出FPGA实时处理模块各步骤的处理时间。测算各条测线数据处理的平均值，其中背景去噪处理平均需要4 213个时钟延迟，一维滤波处理平均需要6 857个时钟延迟，增益处理平均需要5 236个时钟延迟，小波变换平均需要7 236个时钟延迟。因本地质雷达系统脉冲间隔为432 μs，在此期间可以完成对数据的全部处理步骤，满足实时性要求，实现了数据的快速处理。

4 结论与展望

为了提高传统低频地质雷达探测深度及精度，并能够实现雷达数据实时处理，开发了一种低频地质雷达屏蔽天线系统。首先，通过低频屏蔽和阻抗匹配设计，优化了天线辐射面结构和屏蔽结构参数，增强了天线的聚焦和辐射特性，可有效屏蔽环境干扰信号，提高数据的信噪比；其次，提出了基于FPGA的数据实时处理框架，应用可编程片上系统开发实现了数据内插、去噪、增益等处理算法的集成，实现了数据的实时处理，相比于常规处理方法，雷达图像信息熵平均增加了13.72%，平均梯度平均增幅为12.44%。相比于中国矿业大学(北京)研究团队研发的低频组合天线，本次进一步优化了天线的屏蔽效果，使屏蔽复杂信号干扰能力更强，探测精度进一步提升，并且增加了数据处理的实时性。低频地质雷达天线系统的验证试验结果表明，该低频地质雷达天线系统(频率为50 MHz)在城市环境下其探测深度可达30 m，基于FPGA的处理模块实现了数据的实时处理。

试验研究结果为特殊地质结构的勘探信息准确性、城市地下空间的安全排查和隐患探测提供了新的理论和技术支撑。下一步将对FPGA实时处理的算法进行优化，同时开展更多的现场探测试验，测试天线系统的稳定性和可靠性，并在此过程中对天线硬件系统进行优化。将该研究成果应用于地下矿山开采前的勘探，可进一步对待采煤矿地质信息进行精细化研究。