隧道地质预报探地雷达信号干扰消除方法

刘宗辉，吴一帆，刘保东，刘毛毛，蓝日彦，孙怀凤

1) 广西新发展交通集团有限公司，南宁 530028 2) 广西大学土木建筑工程学院，南宁 530004 3) 山东大学土建与水利学院，济南 250061

4) 南宁城建管廊建设投资有限公司，南宁 530219

探地雷达（GPR）达近年来已成为隧道超前地质预报中最主要的短距离物探手段，其在不良地质探测方面具有分辨率高、结果直观、扫描速度快等其它物探方法无法比拟的优势[1-2]. 探地雷达信号是一种典型的非平稳、时变信号[3]，电磁波在复杂的隧道围岩中传播时存在强烈的吸收衰减、色散[4]，同时由于隧道探测环境中大量系统干扰，使得采集到的雷达反射波数据通常具有“弱信号，强干扰”的特征，给数据的处理和解释带来极大困难. 因此，干扰消除一直是探地雷达隧道地质预报应用中普遍关注而又没有得到很好解决的难题.

随着复信号分析技术的发展，探地雷达信号去噪已从中值滤波、频域滤波、F-K滤波等传统方法的基础上发展为小波变换、曲波变换（curverlet变换）以及经验模态分解（EMD）等方法，如柳钢等[5]、Bao 等[6]、Gan 等[7]的研究成果. 这些方法在一定程度上提高了探地雷达的数据解译精度，但大都仍存在各自难以克服的缺点. 如Ouadfeul等[8]、李才明等[9]采用的小波分析方法尽管充分利用了小波变换多尺度分析特性，但是该方法的有效性依赖于选定的小波基函数和设定的软硬阀值，不能适用于复杂多变的探地雷达实测环境，且常用的二维小波对二维信号中直线或曲线等边缘特征难以精确表达. 经验模态分解[10-12]方法尽管充分利用了经验模态分解分解的多分辨率和局部时频分析特性，但已有的研究成果基本都是针对指定探地雷达信号，通过经验人员分析，才能使经验模态分解分解技术达到探地雷达信号降噪的目的，同时由于探地雷达信号是超宽带信号，其中目标体一次回波信号、多次回波信号、杂波和各种噪声信号会发生严重的频率混叠，使得经验模态分解分解不能在探地雷达探测中对接收信号进行有效增强处理. 曲波变换[13-14]是在小波变换的基础上，增加了一个方位参数，解决了小波变换在处理二维信号时的不足，但该方法在干扰信息与有效信息方向性一致的情况下，分离效果并不理想. 曲波变换阀值函数的选取直接关系到曲波降噪效果，找到同时适应尺度和角度的阀值函数是该方法成功应用亟需解决的问题. 此外，目前探地雷达信号去噪研究重心多集中在对算法的更新而忽略了对干扰数据本身特征的研究，虽然有学者研究总结了隧道超前地质预报中的干扰类型[15]，但没有从信号处理角度给出具有针对性的干扰压制方法.

剪切变换（shearlet变换，ST）是一种较新的多尺度多方向时频分析技术，它相比于小波变换（wavelet变换，WT）和其他多尺度几何分析方法有更好的方向敏感性，信号保真度高，该技术在地震波消噪方面已被证实有较好的应用效果[16-17]. 由于地震波和雷达电磁波的诸多相似性，本文引进剪切变换，利用数学理论框架将其改进，并将其与小波变换相结合，提出剪切变换与小波变换联合去除干扰方法. 实际案例处理效果表明该方法在保证去干扰效果同时能较好地保留有效信号.

1 基于剪切变换的自适应阀值去噪

1.1 基本方法原理

1.1.1 小波变换基本原理

小波变换是在傅立叶变换基础上发展起来的，其最大优势是可以由粗到细逐步观察信号，能很好地表征信号局部特征，对于信号频率具有很高的敏感性.

连续小波变换的表达式为：

连续小波变换逆变换表达式为：

小波变换处理数据时基函数的选择尤为重要，小波波形越接近待处理的瞬态信号波形，处理效果越理想. 综合时频域的分辨率来看，DB族小波是比较适合分析探地雷达信号的一种小波基函数[15].

1.1.2 剪切变换基本原理

剪切变换又称剪切波变换，最初是由Guo和Easley等[18-19]通过对剪切基函数进行缩放、剪切和平移合成具有膨胀性的仿射系统. 当维数为2时，该合成膨胀仿射系统定义为：

按照以下方式采样，可将剪切变换离散化：

离散剪切系统为：

1.2 自适应阀值去噪

含噪雷达信号经过剪切变换后，数据信息被映射到不同尺度不同方向的剪切系数上，通常有效信号会根据自身特点集中在一些特定方向上，而随机噪声信号不具有方向性. 有效信号所映射的某些特定方向上的剪切系数值往往较大，而随机噪声将在各个方向上分布，其对应的剪切系数值往往较小.

从阀值去噪角度来讲，如果某一分解方向中剪切系数值较大，则该方向为有效信号系数主要集中的方向，对于该方向上的系数处理应该采取较小的阀值以便更好地保护有效信号；当某一分解方向中剪切系数值较小时，则该方向的系数主要是噪声系数，应该采取较大的阀值以便能压制更多的噪声. 因此，在经典剪切算法中[16]，仅设置一个随尺度变化的阀值难以满足方向上能量变化的需求，容易产生“过扼杀”、“除不净”的现象.

本文在经典尺度阀值函数的基础上，根据有效信号和干扰信号在剪切域中不同尺度、不同方向上能量的差异，设置一个随尺度和方向变化的阀值函数，从而使得在剪切域进行噪声压制时每一个子带都可以根据其能量特征自适应选择最优阀值，其表达式如下：

基于剪切变换的自适应阀值去噪具体过程包括以下步骤：

（1）信号常规处理，包括：直达波去除、直流去漂移、信号增益等；

（2）格式转换，将雷达数据格式转换为矩阵形式以保证算法的实现；

（3）伪极坐标变换，将数据从笛卡尔坐标系转换到伪极坐标系，并且在伪极坐标系上进行多尺度划分，并生成剪切基函数对数据进行窗口子带化；

（4）计算剪切系数，得到系数矩阵之后通过上述自适应阀值函数对剪切系数矩阵进行干扰压制处理；

（5）根据去噪处理后的剪切系数对探地雷达二维数据进行重构.

2 随机干扰消除

2.1 正演模拟及干扰设置

利用时域有限差分法分别模拟圆形空洞与方形空洞两种异常体，得到纯净的雷达数据. 图1和图2分别为异常体几何模型和波形堆积图，从图中可以看出异常体反射界面清晰，同相轴特征明显.

图1 正演几何模型. （a）圆形空洞；（b）方形空洞Fig.1 Geometric model of forward simulation: (a) circular hole; (b) square hole

图2 正演模拟结果. （a）圆形空洞；（b）方形空洞Fig.2 Forward simulation results: (a) circular hole; (b) square hole

在正演模拟获得的纯净数据中加入高斯白噪声，加噪后图像信噪比为-2.5 dB左右，波形堆积图如图3所示. 从图中可以看出异常体反射界面的有效信号已基本被淹没在噪声之中，圆形空洞下反射面同相轴几乎不可见，方形空洞数据只能观察到水平部分能量特别强的部位.

图3 加噪后数据. （a）圆形空洞；（b）方形空洞Fig.3 Data with random interference: (a) circular hole; (b) square hole

2.2 去噪效果分析

2.2.1 波形堆积图对比

对上述加噪后的数据分别运用小波变换与本文所提出的基于自适应阀值的剪切变换进行去噪处理，经过反复对比确定两种方法具体参数为：小波变换选取DB4小波基，分解尺度为4；剪切变换采样率设置为2，分解尺度为4，方向数设置为可分解的最多方向，剪切逆变换采用迭代法求逆矩阵，迭代总数为10，误差限值为10−5. 处理后的波形堆积如图4和图5所示.

从图4和图5中可以很直观的看出小波变换处理后同相轴信息并没有被完全还原出来，圆形空洞数据的下反射面同相轴模糊，方形空洞数据上反射面边缘不清晰，数据整体依然含有较多噪声，去噪效果并不理想. 而使用本文提出的剪切变换去噪后的数据图像，目标信号被极大限度的还原，圆形空洞数据中被噪声掩盖的下反射面同相轴在经过处理后清晰可见，几乎与原始信号数据没有差异，方形空洞的起止位置处的信号也被有效的还原出来. 从方形空洞处理结果同时可以看出剪切变换处理后数据与原数据相比信号强度稍弱，说明存在少量的有效信号被过度剔除.

图4 小波变换处理结果. （a）圆形空洞；（b）方形空洞Fig.4 Results after Wavelet transform processing: (a) circular hole; (b) square hole

图5 剪切变换处理结果. （a）圆形空洞；（b）方形空洞Fig.5 Results after shearlet transform processing: (a) circular hole; (b) square hole

利用信噪比（SNR）、峰值信噪比（PSNR）和均方误差（MSE）三个指标对去噪前后波形堆积图数据质量进行定量分析. 由表1可知，剪切变换处理后数据的信噪比与峰值信噪比均大于小波变换处理结果，而均方误差由个位数缩小到0.1以内，进一步说明剪切变换自适应阀值法在处理随机噪声上的优势，数据噪声残留少.

2.2.2 单道波对比

以圆形空洞的第200道单道波数据为例，对小波变换与剪切变换两种方法的处理效果做进一步对比分析，图 6（a）、6（b）和 6（c）分别为原始数据与两种方法处理后单道波对比图. 从图中可以看出加入高强度噪声后原始信号已基本被淹没，在经过DB4小波去噪处理后有效信号波形大致的形状已经开始显现，但与原始信号曲线对比，依然存在很多杂波的干扰，曲线的平滑度不够，有效信号没有被突出. 而经剪切变换处理后，单道波波形图更为平滑，仅在部分位置存在有少量的杂波，有效信号被保留并且被凸显出来，噪声被有效的压制.

表1 小波变换与剪切变换处理前后信噪比、峰值信噪比、均方误差对比表Table 1 Comparison of SNR, PSNR and MSE before and after wavelet and shearlet transform processing

图6 第 200 道单道波数据去噪前后对比. （a）原始数据；（b）小波变换；（c）剪切变换Fig.6 Comparison before and after denoizing of the 200th A-scan: (a) raw data; (b) wavelet transform; (c) shearlet transform

综上所述，本文提出基于自适应阀值的剪切变换方法可以很好地适应探地雷达的数据结构，对随机噪声有很好的压制能力，可以提高含噪数据信噪比.

3 频率异常干扰消除

探地雷达脉冲信号是一种宽带电磁波，具有非平稳性、非线性衰减等特点. 在隧道掌子面探测时，由于外界信号干扰、地下不同介质的吸收、反射等因素，导致仪器所接收到的反射回波往往是多种频率成分的信号叠加. 小波变换具有良好的时频局部化性质，其可变的时窗结构对一维信号具有较强的频率分辨率. 本文采用的小波去除频率异常信号的基本原理是：首先对信号进行若干尺度层小波分解，获得各尺度层不同频率的小波系数；再结合时频分析结果，取合适的阀值进行带通滤波，去除频率异常部分的信号小波系数；最后再进行剩余小波系数分量重构，从而达到频率异常信号消除.

以岑溪大隧道左洞掌子面DK7+511处探地雷达实测数据来说明上述两种方法对能量接近的不同频率成分信号混叠干扰的去除效果. 现场探测采用意大利IDS公司K2雷达，天线频率为100 MHz，时窗600 ns，采样点数512.

图7（a）为常规方法处理后的探地雷达图像，从图中可以看到明显贯穿整个剖面的强反射波组，对探测范围内地质异常解释存在较大干扰.图8（a）为该测线8.3 m处单道波的广义S变换后的时频分布[20]，从图中也可以看出沿时间深度存在明显低频成分. 结合现场探测环境以及实际开挖情况可判定此处低频成分为干扰信号. 图7（b）和7（c）分别为剪切变换和小波变换处理后二维雷达剖面，处理过程中具体参数与模型试验一致，对比两幅图可以看出剪切变换能将杂波信号很好去除，但对低频干扰信号处理效果并不明显；而小波变换处理后波长异常区域得到明显改善，但数据浅部仍存在较多的杂波信号.

图7 频率异常信号干扰消除前后灰度图. （a）常规方法；（b）剪切变换；（c）小波变换Fig.7 Image of before and after elimination of abnormal frequency signal interference: (a) conventional method; (b) shearlet transform; (c) wavelet transform

图8（b）和 8（c）分别为对应测线 8.3 m 处单道波的时频分布. 对比两幅图可以更直观看出小波变换能很好地将低频成分分离，对频率异常信号去除优势明显.

4 小波变换与剪切变换联合法去干扰

图8 频率异常干扰消除前后单道波广义S变换时频分布. （a）常规方法；（b）剪切变换；（c）小波变换Fig.8 Generalized S transform (GST) spectrogram of GPR A-scan before and after eliminating the interference: (a) conventional method; (b) shearlet;(c) wavelet

根据干扰信号特点，可以将隧道内常见的干扰分为随机干扰和频率异常干扰两种类型，其中随机干扰具有频率随机分布、波形杂乱的特点，而频率异常干扰往往具有某些特定频率特征、在波形堆积图上往往具有某些特定规律. 通过对正演模拟和现场实际数据的处理可以发现，剪切变换与小波变换对干扰的压制都有着各自的优势：剪切变换对信号能量敏感，对于随机噪声、机械噪声以及电磁信号等频率随机、能量异常的干扰压制效果较好，而小波变换对频率异常、能量相近的干扰压制效果较好.

结合两种方法各自的优势，进一步提出小波变换与剪切变换相结合去除干扰方法，即先用小波变换对异常频率信号进行分离，再使用剪切变换对随机干扰进行压制. 具体实现流程图如图9所示.

图9 联合法去除干扰流程图Fig.9 Flow chart of the combined methods for interference removal

5 实际工程案例

现场试验场地位于广西壮族自治区融水县至河池市在建高速公路罗城段，探测对象为路基边缘三处形态不规则裸露的溶洞，其中1#为软塑黏土夹碎石充填型溶洞、2#为干土充填型溶洞、3#为空腔型溶洞. 路基为微风化灰岩、岩质坚硬、结构面较发育，周边未见地表水. 现场情况、测线与溶洞平面示意图分别如图10和图11所示，试验坑长11 m，深1.5 m，实际测线长度10 m.

图10 现场情况Fig.10 Field conditions

图11 测线与溶洞平面示意图Fig.11 Layout diagram of karst caves and survey line

采用意大利IDS公司K2雷达探测，天线频率为100 MHz，时窗400 ns，采样点数1024. 通过人为设置机械电磁噪声和金属干扰来模拟隧道干扰环境.

图12（a）为常规方法处理后的波形堆积图，从图中可以明显看出240 ns以下深度中存在强烈同相轴异常干扰，3#溶洞反射波信号被淹没于干扰信号中，根据信号特征可判定该区域同时含有随机干扰和频率异常成分干扰. 此外，从图中可以看出240 ns以上的浅部数据也存在波形杂乱无章的随机噪声干扰，左侧1#和中间2#溶洞位置难以区分. 图13（a）为测线3 m处单道波时频分布，从图中可清晰看出干扰信号频率成分. 因此，为凸显异常体空间位置以及后续进一步开展属性分析，有必要采用本文所提出的联合去噪方法提取异常体反射信号.

图12 去噪效果对比. （a）常规方法；（b）小波变换；（c）联合算法Fig.12 Comparison of different denoizing methods: (a) conventional method; (b) wavelet transform; (c) joint algorithm

图13 广义 S 变换时频分布结果对比. （a）常规方法；（b）小波变换；（c）联合算法Fig.13 Comparison of GST results obtained using different denoizing methods: (a) conventional method; (b) wavelet transform; (c) joint algorithm

图12（b）和 12（c）分别为使用小波变换以及联合法处理后的波形堆积图，对处理后的雷达数据选取测线3 m处单道波进行时频分析，时频分布见图 13（b）和 13（c）. 从图 12（b）中可以明显看出低频成分的强反射同相轴已得到较好地去除，但随机干扰并未得到有效去除. 而图12（c）处理结果显示整个数据的图像质量有了很大的改善，同相轴变得清晰连续，深层以及浅层的随机噪声和低频干扰信号都能很好的去除. 根据处理后的波形堆积图数据可以很好地将三处地质异常进行圈定，并能与实际情况对应. 图 13（b）和 13（c）时频分布图中同样可进一步清晰看出干扰信号频率成分去除效果. 通过该案例可以进一步说明小波变换与剪切变换联合法去干扰的有效性与必要性.

6 结语

（1）利用剪切变换将探地雷达数据转换到剪切域，可以更加细致的对信号进行多尺度多方向划分，通常有效信号会根据自身特点集中在一些特定的方向上，在此基础上提出的基于自适应阀值去噪方法可对随机干扰有很好去除效果.

（2）小波变换具有良好的时频局部化性质，其可变的时频窗结构对于一维信号有较强的频率分辨率，通过多尺度分解后可将频率异常信号分离，从而保留有效信号. 但小波变换对于随机分布在整个频率域且能量异常的信号去除效果不如剪切变换.

（3）隧道空间环境复杂，使用探地雷达进行隧道超前地质预报时往往会遇到各种成分的干扰混叠，采用单一干扰去除方法难以获得满意效果. 提出的剪切变换与小波变换联合方法可以同时对能量异常的随机干扰以及频率异常的干扰信号进行压制，并且可以保证处理后的数据有着较高信噪比.

（4）根据具体的干扰类型以及干扰数据特征选择合适的阀值系数，直接关系到干扰去除和有效信号保留效果，自适应阀值函数仍有许多不确定性. 实际工程应用中需要不断总结隧道中常见各类干扰的剪切系数能量特征以及小波系数频率特征，并形成干扰信号样本库.