工程科学学报.第42卷.第3期：390-398.2020年3月 Chinese Journal of Engineering,Vol.42,No.3:390-398,March 2020 https://doi.org/10.13374/j.issn2095-9389.2019.04.12.001;http://cje.ustb.edu.cn 隧道地质预报探地雷达信号干扰消除方法 刘宗辉2，)，吴一帆2，刘保东，刘毛毛》，蓝日彦)区，孙怀凤) 1)广西新发展交通集团有限公司，南宁5300282)广西大学土木建筑工程学院.南宁5300043)山东大学土建与水利学院，济南250061 4)南宁城建管廊建设投资有限公司，南宁530219 ☒通信作者，E-mail:47573043@qq.com 摘要受探测环境制约，隧道超前地质预报过程中探地雷达反射波往往具有“弱信号，强干扰”的特征，给数据处理和解译 带来极大的困难.将剪切变换(shearlet变换，ST)引入探地雷达信号处理，根据有效信号和干扰信号在剪切域中不同尺度、不 同方向上的能量差异，提出一种基于自适应阀值的随机干扰去除方法，并通过正演模拟数据验证了该方法在随机干扰去除上 的优势：在此基础上针对隧道超前地质预报中常见的能量接近、频率异常干扰信号，以实际数据为例说明小波变换(WT)对 其去除效果：从而进一步提出小波变换与剪切变换联合干扰压制方法，即首先使用小波变换对异常频率干扰进行分离，然后 采用基于自适应阀值的剪切变换对随机干扰进行压制.现场溶洞探测案例应用效果表明.本文所提出的方法能在去除干扰 的同时很好地保留有效信号，根据处理后的波形堆积图可以很好地凸显地质异常区域，从而提高探地雷达资料解译精度 关键词隧道地质预报：探地雷达；剪切变换：小波变换：干扰消除 分类号TU375.4 Research on the interference elimination method of GPR signal for tunnel geological prediction LIU Zong-hui2),WU Yi-fan,LIU Bao-dong,LIU Mao-mao,LAN Ri-yan SUN Huai-feng 1)Guangxi Xinfazhan Communications Group Co Ltd.,Nanning 530028,China 2)College of Civil Engineering and Architecture,Guangxi University,Nanning 530004,China 3)School of Civil Engineering,Shandong University,Jinan 250061,China 4)Nanning Pipe Gallery Construction Investment Co.,Ltd.,Nanning 530219,China Corresponding author,E-mail:47573043@qq.com ABSTRACT Ground-penetrating radar(GPR)has been used in a wide range of shallow detection applications,such as underground geological mapping,highway detection,and hydrogeology survey.In recent years,GPR has been most widely utilized in tunnel geological prediction because it has the advantages of high resolution,intuitionistic results,and fast scanning.In addition,GPR signal is a typical nonstationary and time-varying signal,with its electromagnetic wave exhibiting strong absorption attenuation and dispersion as it propagated in complex surrounding rock.At the same time,the GPR response is often characterized by a weak signal and a strong interference because of numerous system interferences in the tunnel detection environment,which lead to difficulties in data processing and interpretation.Therefore,interference elimination is always a difficult problem when GPR is applied to tunnel geological prediction. In this study,through the introduction of shearlet transform(ST)to GPR signal processing,an adaptive thresholding method is proposed to eliminate random interference on the basis of the energy difference between effective and interference signals in the shearlet domain at different scales and directions.The advantages of this method in random interference removal are verified by forward simulation data. 收稿日期：2019-04-12 基金项目：国家自然科学基金资助项目(51708136)：广西自然科学基金资助项目(2017 GXNSFBA198199):中国博士后基金面上项目 (2019M653310):广西科技基地和人才专项资助项目（桂科AD19245153,桂科AD17129047)
隧道地质预报探地雷达信号干扰消除方法 刘宗辉1,2,3)，吴一帆2)，刘保东4)，刘毛毛2)，蓝日彦1) 苣，孙怀凤3) 1) 广西新发展交通集团有限公司，南宁 530028 2) 广西大学土木建筑工程学院，南宁 530004 3) 山东大学土建与水利学院，济南 250061 4) 南宁城建管廊建设投资有限公司，南宁 530219 苣通信作者，E-mail：47573043@qq.com 摘 要 受探测环境制约，隧道超前地质预报过程中探地雷达反射波往往具有“弱信号，强干扰”的特征，给数据处理和解译 带来极大的困难. 将剪切变换（shearlet 变换，ST）引入探地雷达信号处理，根据有效信号和干扰信号在剪切域中不同尺度、不 同方向上的能量差异，提出一种基于自适应阀值的随机干扰去除方法，并通过正演模拟数据验证了该方法在随机干扰去除上 的优势；在此基础上针对隧道超前地质预报中常见的能量接近、频率异常干扰信号，以实际数据为例说明小波变换（WT）对 其去除效果；从而进一步提出小波变换与剪切变换联合干扰压制方法，即首先使用小波变换对异常频率干扰进行分离，然后 采用基于自适应阀值的剪切变换对随机干扰进行压制. 现场溶洞探测案例应用效果表明，本文所提出的方法能在去除干扰 的同时很好地保留有效信号，根据处理后的波形堆积图可以很好地凸显地质异常区域，从而提高探地雷达资料解译精度. 关键词 隧道地质预报；探地雷达；剪切变换；小波变换；干扰消除 分类号 TU375.4 Research on the interference elimination method of GPR signal for tunnel geological prediction LIU Zong-hui1,2,3) ，WU Yi-fan2) ，LIU Bao-dong4) ，LIU Mao-mao2) ，LAN Ri-yan1) 苣 ，SUN Huai-feng3) 1) Guangxi Xinfazhan Communications Group Co Ltd., Nanning 530028, China 2) College of Civil Engineering and Architecture, Guangxi University, Nanning 530004, China 3) School of Civil Engineering, Shandong University, Jinan 250061, China 4) Nanning Pipe Gallery Construction Investment Co., Ltd., Nanning 530219, China 苣 Corresponding author, E-mail: 47573043@qq.com ABSTRACT Ground-penetrating radar (GPR) has been used in a wide range of shallow detection applications, such as underground geological mapping, highway detection, and hydrogeology survey. In recent years, GPR has been most widely utilized in tunnel geological prediction because it has the advantages of high resolution, intuitionistic results, and fast scanning. In addition, GPR signal is a typical nonstationary and time-varying signal, with its electromagnetic wave exhibiting strong absorption attenuation and dispersion as it propagated in complex surrounding rock. At the same time, the GPR response is often characterized by a weak signal and a strong interference because of numerous system interferences in the tunnel detection environment, which lead to difficulties in data processing and interpretation. Therefore, interference elimination is always a difficult problem when GPR is applied to tunnel geological prediction. In this study, through the introduction of shearlet transform (ST) to GPR signal processing, an adaptive thresholding method is proposed to eliminate random interference on the basis of the energy difference between effective and interference signals in the shearlet domain at different scales and directions. The advantages of this method in random interference removal are verified by forward simulation data. 收稿日期: 2019−04−12 基金项目: 国家自然科学基金资助项目（ 51708136）；广西自然科学基金资助项目（ 2017GXNSFBA198199）；中国博士后基金面上项目 （2019M653310）；广西科技基地和人才专项资助项目（桂科 AD19245153，桂科 AD17129047） 工程科学学报，第 42 卷，第 3 期：390−398，2020 年 3 月 Chinese Journal of Engineering, Vol. 42, No. 3: 390−398, March 2020 https://doi.org/10.13374/j.issn2095-9389.2019.04.12.001; http://cje.ustb.edu.cn
刘宗辉等：隧道地质预报探地雷达信号干扰消除方法 391· On this basis,the interference signal,as well as its energy proximity and frequency anomaly,common in advanced tunnel geological prediction is taken as an example to illustrate the effect of wavelet transform(WT)on its removal.In this manner,WT and ST are combined to suppress interference.First,WT is used to separate abnormal frequency interference.Then,ST based on the adaptive thresholding method is used to suppress random interference.The results of practical engineering cases of karst cave detection in the field show that the method proposed in this study can remove the interference signal,retain the effective signal,and highlight the abnormal geological area on the basis of the processed waveform stacking diagram to improve the interpretation accuracy of GPR data. KEY WORDS tunnel geological prediction;ground-penetrating radar;shearlet transform:wavelet transform:interference suppression 探地雷达(GPR)达近年来已成为遂道超前地 果，找到同时适应尺度和角度的阀值函数是该 质预报中最主要的短距离物探手段，其在不良地 方法成功应用亟需解决的问题.此外，目前探地雷 质探测方面具有分辨率高、结果直观、扫描速度 达信号去噪研究重心多集中在对算法的更新而忽 快等其它物探方法无法比拟的优势-探地雷达 略了对干扰数据本身特征的研究，虽然有学者研 信号是一种典型的非平稳、时变信号)，电磁波在 究总结了隧道超前地质预报中的干扰类型)，但 复杂的隧道围岩中传播时存在强烈的吸收衰减、 没有从信号处理角度给出具有针对性的干扰压制 色散，同时由于隧道探测环境中大量系统干扰， 方法 使得采集到的雷达反射波数据通常具有“弱信号， 剪切变换(shearlet变换，ST)是一种较新的多 强干扰”的特征，给数据的处理和解释带来极大困 尺度多方向时频分析技术，它相比于小波变换 难.因此，干扰消除一直是探地雷达隧道地质预报 (wavelet变换，WT)和其他多尺度几何分析方法有 应用中普遍关注而又没有得到很好解决的难题. 更好的方向敏感性，信号保真度高，该技术在地震 随着复信号分析技术的发展，探地雷达信号 波消噪方面已被证实有较好的应用效果61m.由 去噪已从中值滤波、频域滤波、FK滤波等传统方 于地震波和雷达电磁波的诸多相似性，本文引进 法的基础上发展为小波变换、曲波变换(curverlet 剪切变换，利用数学理论框架将其改进，并将其与 变换)以及经验模态分解(EMD)等方法，如柳钢 小波变换相结合，提出剪切变换与小波变换联合 等、Bao等阿、Gan等m的研究成果.这些方法在 去除干扰方法.实际案例处理效果表明该方法在 一定程度上提高了探地雷达的数据解译精度，但 保证去干扰效果同时能较好地保留有效信号 大都仍存在各自难以克服的缺点.如Ouadfeul等阁、 李才明等四采用的小波分析方法尽管充分利用了 1 基于剪切变换的自适应阀值去噪 小波变换多尺度分析特性，但是该方法的有效性 1.1基本方法原理 依赖于选定的小波基函数和设定的软硬阀值，不 1.1.1小波变换基本原理 能适用于复杂多变的探地雷达实测环境，且常用 小波变换是在傅立叶变换基础上发展起来 的二维小波对二维信号中直线或曲线等边缘特征 的，其最大优势是可以由粗到细逐步观察信号，能 难以精确表达.经验模态分解012)方法尽管充分 很好地表征信号局部特征，对于信号频率具有很 利用了经验模态分解分解的多分辨率和局部时频 高的敏感性 分析特性，但已有的研究成果基本都是针对指定 连续小波变换的表达式为： 探地雷达信号，通过经验人员分析，才能使经验模 态分解分解技术达到探地雷达信号降噪的目的， WTf(a,b)=(f(t),ab(t))= S.)a 同时由于探地雷达信号是超宽带信号，其中目标 (1) 体一次回波信号、多次回波信号、杂波和各种噪 式中：9ab()为小波基函数，通过母函数()平移、 声信号会发生严重的频率混叠，使得经验模态分 伸缩得到，α、b分别为尺度因子和平移因子，可调 解分解不能在探地雷达探测中对接收信号进行有 节时窗和频窗的大小，从而控制时间和频率的分 效增强处理.曲波变换3是在小波变换的基础 辨率，小波变换的时频窗结构能很好地满足超宽 上，增加了一个方位参数，解决了小波变换在处理 频带雷达信号去噪要求 二维信号时的不足，但该方法在干扰信息与有效 连续小波变换逆变换表达式为： 信息方向性一致的情况下，分离效果并不理想.曲 f(t)= WTr(a,b)pab(t)dadb (2) 波变换阀值函数的选取直接关系到曲波降噪效
On this basis, the interference signal, as well as its energy proximity and frequency anomaly, common in advanced tunnel geological prediction is taken as an example to illustrate the effect of wavelet transform (WT) on its removal. In this manner, WT and ST are combined to suppress interference. First, WT is used to separate abnormal frequency interference. Then, ST based on the adaptive thresholding method is used to suppress random interference. The results of practical engineering cases of karst cave detection in the field show that the method proposed in this study can remove the interference signal, retain the effective signal, and highlight the abnormal geological area on the basis of the processed waveform stacking diagram to improve the interpretation accuracy of GPR data. KEY WORDS tunnel geological prediction；ground-penetrating radar；shearlet transform；wavelet transform；interference suppression 探地雷达（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 小波变换基本原理 小波变换是在傅立叶变换基础上发展起来 的，其最大优势是可以由粗到细逐步观察信号，能 很好地表征信号局部特征，对于信号频率具有很 高的敏感性. 连续小波变换的表达式为： WTf(a,b) = ⟨f(t),φa,b(t)⟩ = 1 √ a w R f(t)φ ( t−b a ) dt （1） 式中： φa,b(t) 为小波基函数，通过母函数 φ(t) 平移、 伸缩得到，a、b 分别为尺度因子和平移因子，可调 节时窗和频窗的大小，从而控制时间和频率的分 辨率，小波变换的时频窗结构能很好地满足超宽 频带雷达信号去噪要求. 连续小波变换逆变换表达式为： f(t) = 1 Cφ x R 1 a 2 WTf (a,b)φa,b (t)dadb （2） 刘宗辉等： 隧道地质预报探地雷达信号干扰消除方法 · 391 ·
392 工程科学学报，第42卷，第3期 ufdo<,u)为40的FFT变换 式中：Cp=JRm 较小的阀值以便更好地保护有效信号；当某一分 小波变换处理数据时基函数的选择尤为重 解方向中剪切系数值较小时，则该方向的系数主 要，小波波形越接近待处理的瞬态信号波形，处理 要是噪声系数，应该采取较大的阀值以便能压制 效果越理想.综合时频域的分辨率来看，DB族小 更多的噪声.因此，在经典剪切算法中6，仅设置 波是比较适合分析探地雷达信号的一种小波基 一个随尺度变化的阀值难以满足方向上能量变化 函数啊 的需求，容易产生“过扼杀”、“除不净”的现象 1.1.2剪切变换基本原理 本文在经典尺度阀值函数的基础上，根据有 剪切变换又称剪切波变换，最初是由Guo和 效信号和干扰信号在剪切域中不同尺度、不同方 Easley等l&-l9通过对剪切基函数进行缩放、剪切 向上能量的差异，设置一个随尺度和方向变化的 和平移合成具有膨胀性的仿射系统.当维数为 阀值函数，从而使得在剪切域进行噪声压制时每 2时，该合成膨胀仿射系统定义为： 一个子带都可以根据其能量特征自适应选择最优 阀值，其表达式如下： pAB(）={j(m)=ldetA(BAx-k):leZ,k∈Z2 (3) Thin=oV2logN×2-J-D cos 式中：山∈L2(R2),A、B为2阶的可逆矩阵且 ldet B=l,矩阵A控制尺度参数，可进行伸缩放大 j=1,2,J;1=1,2,.Li 操作：矩阵B控制着方向参数，可进行旋转剪切 (7) 操作 式中：σ为噪声标准差，N为信号长度，ε为任意小 当中AB(W)满足Parseval框架时，称该系统为复 的正数，L为尺度j上方向分解总数，e表示尺度 合小波系统.当A=A,B=B,时，复合小波系统进 j方向1上系数的能量，max(e)表示尺度j上各分 一步简化为剪切波系统.其中 解方向中系数的最大能量 A0=(各9)B=(01) 基于剪切变换的自适应阀值去噪具体过程包 括以下步骤： 对于已知信号f∈L2(R),其连续剪切变换定 (1)信号常规处理，包括：直达波去除、直流去 义为： 漂移、信号增益等； Sril,k=〈f,),(,k)∈S (4) (2)格式转换，将雷达数据格式转换为矩阵形 式中：j>0是尺度参数，s∈R是剪切参数，k∈R2是 式以保证算法的实现： 位置参数 (3)伪极坐标变换，将数据从笛卡尔坐标系 按照以下方式采样，可将剪切变换离散化： 转换到伪极坐标系，并且在伪极坐标系上进行 aj=22j,jeZ:sk=k2,-2j≤1≤2 (5) 多尺度划分，并生成剪切基函数对数据进行窗口 离散剪切系统为： 子带化： {w9(w=232w0(B'Ax-k):j≥0， (4)计算剪切系数，得到系数矩阵之后通过上 述自适应阀值函数对剪切系数矩阵进行干扰压制 -2J≤1≤2J-1,k∈Z2 (6) 处理； 12自适应阀值去噪 (5)根据去噪处理后的剪切系数对探地雷达 含噪雷达信号经过剪切变换后，数据信息被 二维数据进行重构 映射到不同尺度不同方向的剪切系数上，通常有 2 随机干扰消除 效信号会根据自身特点集中在一些特定方向上， 而随机噪声信号不具有方向性，有效信号所映射 2.1正演模拟及干扰设置 的某些特定方向上的剪切系数值往往较大，而随 利用时域有限差分法分别模拟圆形空洞与方 机噪声将在各个方向上分布，其对应的剪切系数 形空洞两种异常体，得到纯净的雷达数据.图1和 值往往较小. 图2分别为异常体几何模型和波形堆积图，从图中 从阀值去噪角度来讲，如果某一分解方向中 可以看出异常体反射界面清晰，同相轴特征明显 剪切系数值较大，则该方向为有效信号系数主要 在正演模拟获得的纯净数据中加入高斯白噪 集中的方向，对于该方向上的系数处理应该采取 声，加噪后图像信噪比为-2.5dB左右，波形堆积
Cφ = r R |Ψ(ω)| 2 |ω| 式中： dω < ∞，Ψ(ω) 为 φ(t) 的 FFT 变换. 小波变换处理数据时基函数的选择尤为重 要，小波波形越接近待处理的瞬态信号波形，处理 效果越理想. 综合时频域的分辨率来看，DB 族小 波是比较适合分析探地雷达信号的一种小波基 函数[15] . 1.1.2 剪切变换基本原理 剪切变换又称剪切波变换，最初是由 Guo 和 Easley 等[18‒ 19] 通过对剪切基函数进行缩放、剪切 和平移合成具有膨胀性的仿射系统. 当维数为 2 时，该合成膨胀仿射系统定义为： ϕAB(ψ) = { ψj,l,k(x) = |detA| j 2 ψ(B lA j x−k) : j,l ∈ Z, k ∈ Z 2 } （3） ψ ∈ L 2 ( R 2 ) |detB| = 1 式 中 ： ， A、 B 为 2 阶 的 可 逆 矩 阵 且 ，矩阵 A 控制尺度参数，可进行伸缩放大 操作；矩阵 B 控制着方向参数，可进行旋转剪切 操作. 当 ϕAB (ψ) 满足 Parseval 框架时，称该系统为复 合小波系统. 当 A=A0，B=B0 时，复合小波系统进 一步简化为剪切波系统. 其中 A0 = ( 4 0 0 2 ) , B0 = ( 1 1 0 1 ) . f ∈ L 2 对于已知信号 (R) ，其连续剪切变换定 义为： S f(j,l, k) = ⟨ f,ψj,l,k ⟩ ,(j,l, k) ∈ S （4） j > 0 s ∈ R k ∈ R 式中： 是尺度参数， 是剪切参数， 2 是 位置参数. 按照以下方式采样，可将剪切变换离散化： aj = 2 −2 j , j ∈ Z;sj,k = k2 −j ,−2 −j ⩽ l ⩽ 2 −j （5） 离散剪切系统为： { ψ (0) j,l,k (x) = 2 3 j/2ψ (0) (B lA j x−k) : j ⩾ 0, −2 j ⩽ l ⩽ 2 j −1, k ∈ Z 2 } （6） 1.2 自适应阀值去噪 含噪雷达信号经过剪切变换后，数据信息被 映射到不同尺度不同方向的剪切系数上，通常有 效信号会根据自身特点集中在一些特定方向上， 而随机噪声信号不具有方向性. 有效信号所映射 的某些特定方向上的剪切系数值往往较大，而随 机噪声将在各个方向上分布，其对应的剪切系数 值往往较小. 从阀值去噪角度来讲，如果某一分解方向中 剪切系数值较大，则该方向为有效信号系数主要 集中的方向，对于该方向上的系数处理应该采取 较小的阀值以便更好地保护有效信号；当某一分 解方向中剪切系数值较小时，则该方向的系数主 要是噪声系数，应该采取较大的阀值以便能压制 更多的噪声. 因此，在经典剪切算法中[16] ，仅设置 一个随尺度变化的阀值难以满足方向上能量变化 的需求，容易产生“过扼杀”、“除不净”的现象. 本文在经典尺度阀值函数的基础上，根据有 效信号和干扰信号在剪切域中不同尺度、不同方 向上能量的差异，设置一个随尺度和方向变化的 阀值函数，从而使得在剪切域进行噪声压制时每 一个子带都可以根据其能量特征自适应选择最优 阀值，其表达式如下： Thj,l =σ √ 2logN ×2 −(J−j) × cos ej,l max 1⩽l⩽Lj ( ej,l ) +ε , j = 1,2,...J;l = 1,2,...Lj （7） σ ε Lj ej,l max(e j,l) 式中： 为噪声标准差，N 为信号长度， 为任意小 的正数， 为尺度 j 上方向分解总数， 表示尺度 j 方向 l 上系数的能量， 表示尺度 j 上各分 解方向中系数的最大能量. 基于剪切变换的自适应阀值去噪具体过程包 括以下步骤： （1）信号常规处理，包括：直达波去除、直流去 漂移、信号增益等； （2）格式转换，将雷达数据格式转换为矩阵形 式以保证算法的实现； （ 3）伪极坐标变换，将数据从笛卡尔坐标系 转换到伪极坐标系，并且在伪极坐标系上进行 多尺度划分，并生成剪切基函数对数据进行窗口 子带化； （4）计算剪切系数，得到系数矩阵之后通过上 述自适应阀值函数对剪切系数矩阵进行干扰压制 处理； （5）根据去噪处理后的剪切系数对探地雷达 二维数据进行重构. 2 随机干扰消除 2.1 正演模拟及干扰设置 利用时域有限差分法分别模拟圆形空洞与方 形空洞两种异常体，得到纯净的雷达数据. 图 1 和 图 2 分别为异常体几何模型和波形堆积图，从图中 可以看出异常体反射界面清晰，同相轴特征明显. 在正演模拟获得的纯净数据中加入高斯白噪 声，加噪后图像信噪比为‒2.5 dB 左右，波形堆积 · 392 · 工程科学学报，第 42 卷，第 3 期
刘宗辉等：隧道地质预报探地雷达信号干扰消除方法 393· (a) (b) 0.2 0.2 0.4 0.6 0.8 0.8 1.0 1.0 0.5 1.0 1.5 2.0 0.51.01.5 2.0 Position/m Position/m 图1正演几何模型.(a)圆形空洞：(b)方形空洞 Fig.I Geometric model of forward simulation:(a)circular hole;(b)square hole 0 0 al (b) 10 15 15 200 400 600 0 100200300400 500 Channel number Channel number 图2正演模拟结果.(a)圆形空洞：(b)方形空洞 Fig.2 Forward simulation results:(a)circular hole;(b)square hole 图如图3所示.从图中可以看出异常体反射界面 下反射面同相轴几乎不可见，方形空洞数据只能 的有效信号已基本被淹没在噪声之中，圆形空洞 观察到水平部分能量特别强的部位. 0 (a) 5 10 10 15 15 0 200 400 600 0 100200300400 500 Channel number Channel number 图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.处理后的波 效的还原出来.从方形空洞处理结果同时可以看 形堆积如图4和图5所示 出剪切变换处理后数据与原数据相比信号强度稍 从图4和图5中可以很直观的看出小波变换 弱，说明存在少量的有效信号被过度剔除
图如图 3 所示. 从图中可以看出异常体反射界面 的有效信号已基本被淹没在噪声之中，圆形空洞 下反射面同相轴几乎不可见，方形空洞数据只能 观察到水平部分能量特别强的部位. 2.2 去噪效果分析 2.2.1 波形堆积图对比 对上述加噪后的数据分别运用小波变换与本 文所提出的基于自适应阀值的剪切变换进行去噪 处理，经过反复对比确定两种方法具体参数为：小 波变换选取 DB4 小波基，分解尺度为 4；剪切变换 采样率设置为 2，分解尺度为 4，方向数设置为可 分解的最多方向，剪切逆变换采用迭代法求逆矩 阵，迭代总数为 10，误差限值为 10−5 . 处理后的波 形堆积如图 4 和图 5 所示. 从图 4 和图 5 中可以很直观的看出小波变换 处理后同相轴信息并没有被完全还原出来，圆形 空洞数据的下反射面同相轴模糊，方形空洞数据 上反射面边缘不清晰，数据整体依然含有较多噪 声，去噪效果并不理想. 而使用本文提出的剪切变 换去噪后的数据图像，目标信号被极大限度的还 原，圆形空洞数据中被噪声掩盖的下反射面同相 轴在经过处理后清晰可见，几乎与原始信号数据 没有差异，方形空洞的起止位置处的信号也被有 效的还原出来. 从方形空洞处理结果同时可以看 出剪切变换处理后数据与原数据相比信号强度稍 弱，说明存在少量的有效信号被过度剔除. Position/m 0.5 1.5 1.0 2.0 0.2 0.4 0.6 0.8 1.0 Depth/m (a) Position/m 0.5 1.5 1.0 2.0 0.2 0.4 0.6 0.8 1.0 Depth/m (b) 图 1 正演几何模型. （a）圆形空洞；（b）方形空洞 Fig.1 Geometric model of forward simulation: (a) circular hole; (b) square hole Channel number 0 600 200 400 0 5 10 15 Time/ns (a) Channel number 0 300 400 100 200 500 0 5 10 15 Time/ns (b) 图 2 正演模拟结果. （a）圆形空洞；（b）方形空洞 Fig.2 Forward simulation results: (a) circular hole; (b) square hole Channel number 0 600 200 400 0 5 10 15 Time/ns (a) Channel number 0 300 400 100 200 500 0 5 10 15 Time/ns (b) 图 3 加噪后数据. （a）圆形空洞；（b）方形空洞 Fig.3 Data with random interference: (a) circular hole; (b) square hole 刘宗辉等： 隧道地质预报探地雷达信号干扰消除方法 · 393 ·
394 工程科学学报，第42卷，第3期 (a) b 10 15 200 400 600 0 100200300400500 Channel number Channel number 图4小波变换处理结果.(a)圆形空洞：(b)方形空洞 Fig.4 Results after Wavelet transform processing:(a)circular hole;(b)square hole 0 0 (a) (b) 5 10 15 15 200 400 600 0. 100200300400500 Channel number Channel number 图5剪切变换处理结果.(a)圆形空洞：(b)方形空洞 Fig.5 Results after shearlet transform processing:(a)circular hole;(b)square hole 利用信噪比(SNR)、峰值信噪比(PSNR)和均 2.2.2单道波对比 方误差(MSE)三个指标对去噪前后波形堆积图数 以圆形空洞的第200道单道波数据为例，对小 据质量进行定量分析.由表1可知，剪切变换处理 波变换与剪切变换两种方法的处理效果做进一步 后数据的信噪比与峰值信噪比均大于小波变换处 对比分析，图6(a)、6(b)和6(c)分别为原始数据 理结果，而均方误差由个位数缩小到01以内，进 与两种方法处理后单道波对比图.从图中可以看 一步说明剪切变换自适应阀值法在处理随机噪声 出加入高强度噪声后原始信号已基本被淹没，在 上的优势，数据噪声残留少 经过DB4小波去噪处理后有效信号波形大致的形 表1小波变换与剪切变换处理前后信噪比、峰值信噪比、均方误差对比表 Table 1 Comparison of SNR,PSNR and MSE before and after wavelet and shearlet transform processing SNR/dB PSNR MSE Model type Noisy data WT ST Noisy data WT ST Noisy data WT ST Circular -2.457 9.559 22.871 14.144 26.170 39.481 3.347 0.2590.041 Square -2.528 9.469 22.987 14.153 26.150 39.667 4.891 0.2490.002 Wavelet transform Shearlet transform 300 150 Clean data (a) [(b) 15 [(e) Clean data 200 100 100 100 50 0 0 100 -50 -50 -200 -500 -500 12 12 Time/ns Time/ns Time/ns 图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
利用信噪比（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 Model type SNR / dB PSNR MSE Noisy data WT ST Noisy data WT ST Noisy data WT ST Circular ‒2.457 9.559 22.871 14.144 26.170 39.481 3.347 0.259 0.041 Square ‒2.528 9.469 22.987 14.153 26.150 39.667 4.891 0.249 0.002 Channel number 0 600 200 400 0 5 10 15 Time/ns (a) Channel number 0 300 400 100 200 500 0 5 10 15 Time/ns (b) 图 4 小波变换处理结果. （a）圆形空洞；（b）方形空洞 Fig.4 Results after Wavelet transform processing: (a) circular hole; (b) square hole Channel number 0 600 200 400 0 5 10 15 Time/ns (a) Channel number 0 300 400 100 200 500 0 5 10 15 Time/ns (b) 图 5 剪切变换处理结果. （a）圆形空洞；（b）方形空洞 Fig.5 Results after shearlet transform processing: (a) circular hole; (b) square hole Time/ns 8 12 300 200 −100 −200 Amplitude/mV 0 4 (a) 100 0 Noisy data Clean data Time/ns 8 12 150 100 −50 −500 Amplitude/mV 0 4 (b) 50 0 Wavelet transform Clean data Time/ns 8 12 150 100 −50 −500 Amplitude/mV 0 4 (c) 50 0 Shearlet transform Clean data 图 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 · 394 · 工程科学学报，第 42 卷，第 3 期