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

工程科学学报.第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期 ufdo0是尺度参数，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ω 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 期

刘宗辉等：隧道地质预报探地雷达信号干扰消除方法 395· 状已经开始显现，但与原始信号曲线对比，依然存 带通滤波，去除频率异常部分的信号小波系数：最 在很多杂波的干扰，曲线的平滑度不够，有效信号 后再进行剩余小波系数分量重构，从而达到频率 没有被突出.而经剪切变换处理后，单道波波形图 异常信号消除 更为平滑，仅在部分位置存在有少量的杂波，有效 以岑溪大隧道左洞掌子面DK7+511处探地雷 信号被保留并且被凸显出来，噪声被有效的压制. 达实测数据来说明上述两种方法对能量接近的不 综上所述，本文提出基于自适应阀值的剪切 同频率成分信号混叠干扰的去除效果.现场探测 变换方法可以很好地适应探地雷达的数据结构， 采用意大利IDS公司K2雷达，天线频率为100MHz, 对随机噪声有很好的压制能力，可以提高含噪数 时窗600ns,采样点数512. 据信噪比 图7(a)为常规方法处理后的探地雷达图像， 从图中可以看到明显贯穿整个剖面的强反射波 3频率异常干扰消除 组，对探测范围内地质异常解释存在较大干扰 探地雷达脉冲信号是一种宽带电磁波，具有 图8(a)为该测线8.3m处单道波的广义S变换后 非平稳性、非线性衰减等特点.在隧道掌子面探 的时频分布20，从图中也可以看出沿时间深度存 测时，由于外界信号干扰、地下不同介质的吸收、 在明显低频成分.结合现场探测环境以及实际开 反射等因素，导致仪器所接收到的反射回波往往 挖情况可判定此处低频成分为干扰信号.图7(b) 是多种频率成分的信号叠加.小波变换具有良好 和7(c)分别为剪切变换和小波变换处理后二维雷 的时频局部化性质，其可变的时窗结构对一维信 达剖面，处理过程中具体参数与模型试验一致，对 号具有较强的频率分辨率.本文采用的小波去除 比两幅图可以看出剪切变换能将杂波信号很好去 频率异常信号的基本原理是：首先对信号进行若 除，但对低频干扰信号处理效果并不明显；而小波 干尺度层小波分解，获得各尺度层不同频率的小 变换处理后波长异常区域得到明显改善，但数据 波系数；再结合时频分析结果，取合适的阀值进行 浅部仍存在较多的杂波信号 0 a (b) c) 100 100 200 200 200 300 300 400 400 400 500 500 500 8 4 8 12 0 8 12 Position/m Position/m Position/m 图7频率异常信号干扰消除前后灰度图.()常规方法：(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.3m处单道 4小波变换与剪切变换联合法去干扰 波的时频分布.对比两幅图可以更直观看出小波 变换能很好地将低频成分分离，对频率异常信号 根据干扰信号特点，可以将隧道内常见的干 去除优势明显 扰分为随机干扰和频率异常干扰两种类型，其中

状已经开始显现，但与原始信号曲线对比，依然存 在很多杂波的干扰，曲线的平滑度不够，有效信号 没有被突出. 而经剪切变换处理后，单道波波形图 更为平滑，仅在部分位置存在有少量的杂波，有效 信号被保留并且被凸显出来，噪声被有效的压制. 综上所述，本文提出基于自适应阀值的剪切 变换方法可以很好地适应探地雷达的数据结构， 对随机噪声有很好的压制能力，可以提高含噪数 据信噪比. 3    频率异常干扰消除 探地雷达脉冲信号是一种宽带电磁波，具有 非平稳性、非线性衰减等特点. 在隧道掌子面探 测时，由于外界信号干扰、地下不同介质的吸收、 反射等因素，导致仪器所接收到的反射回波往往 是多种频率成分的信号叠加. 小波变换具有良好 的时频局部化性质，其可变的时窗结构对一维信 号具有较强的频率分辨率. 本文采用的小波去除 频率异常信号的基本原理是：首先对信号进行若 干尺度层小波分解，获得各尺度层不同频率的小 波系数；再结合时频分析结果，取合适的阀值进行 带通滤波，去除频率异常部分的信号小波系数；最 后再进行剩余小波系数分量重构，从而达到频率 异常信号消除. 以岑溪大隧道左洞掌子面 DK7+511 处探地雷 达实测数据来说明上述两种方法对能量接近的不 同频率成分信号混叠干扰的去除效果. 现场探测 采用意大利 IDS 公司 K2 雷达，天线频率为 100 MHz， 时窗 600 ns，采样点数 512. 图 7（a）为常规方法处理后的探地雷达图像， 从图中可以看到明显贯穿整个剖面的强反射波 组，对探测范围内地质异常解释存在较大干扰. 图 8（a）为该测线 8.3 m 处单道波的广义 S 变换后 的时频分布[20]，从图中也可以看出沿时间深度存 在明显低频成分. 结合现场探测环境以及实际开 挖情况可判定此处低频成分为干扰信号. 图 7（b） 和 7（c）分别为剪切变换和小波变换处理后二维雷 达剖面，处理过程中具体参数与模型试验一致，对 比两幅图可以看出剪切变换能将杂波信号很好去 除，但对低频干扰信号处理效果并不明显；而小波 变换处理后波长异常区域得到明显改善，但数据 浅部仍存在较多的杂波信号. 图 8（b）和 8（c）分别为对应测线 8.3 m 处单道 波的时频分布. 对比两幅图可以更直观看出小波 变换能很好地将低频成分分离，对频率异常信号 去除优势明显. 4    小波变换与剪切变换联合法去干扰 根据干扰信号特点，可以将隧道内常见的干 扰分为随机干扰和频率异常干扰两种类型，其中 Time/ns 0 100 200 300 400 500 Position/m 0 4 8 12 (a) Time/ns 0 100 200 300 400 500 Position/m 0 4 8 12 (b) Time/ns 0 100 200 300 400 500 Position/m 0 4 8 12 (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 刘宗辉等： 隧道地质预报探地雷达信号干扰消除方法 · 395 ·

396 工程科学学报，第42卷，第3期 200 (a) 200 (b) 200 (c) ZHW/Xouanbaid 50 150 100 100 50 0 0 0 200 400 600 200 400 600 200 400 600 Time/ns Time/ns Time/ns 图8颜率异常干扰消除前后单道波广义S变换时频分布.()常规方法：(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 随机干扰具有频率随机分布、波形杂乱的特点，而 5 实际工程案例 频率异常干扰往往具有某些特定频率特征、在波 现场试验场地位于广西壮族自治区融水县至 形堆积图上往往具有某些特定规律.通过对正演 河池市在建高速公路罗城段，探测对象为路基边 模拟和现场实际数据的处理可以发现，剪切变换 与小波变换对干扰的压制都有着各自的优势：剪 缘三处形态不规则裸露的溶洞，其中1#为软塑黏 土夹碎石充填型溶洞、2#为干土充填型溶洞、3#为 切变换对信号能量敏感，对于随机噪声、机械噪声 以及电磁信号等频率随机、能量异常的干扰压制 空腔型溶洞.路基为微风化灰岩、岩质坚硬、结构 效果较好，而小波变换对频率异常、能量相近的干 面较发育，周边未见地表水.现场情况、测线与溶 扰压制效果较好 洞平面示意图分别如图10和图11所示，试验坑 结合两种方法各自的优势，进一步提出小波 长11m,深1.5m,实际测线长度10m. 变换与剪切变换相结合去除干扰方法，即先用小 3# 波变换对异常频率信号进行分离，再使用剪切变 换对随机干扰进行压制.具体实现流程图如图9 所示 Data processed by conventional methods Comprehensive analysis based on on-site detection environment,waveform image 图10现场情况 characteristics,and spectral characteristics Fig.10 Field conditions Preliminary estimation of interference types Slope Cavity3 Filled with dry soil Interference with No Filled with abnormal frequency soft plastic clay 12.5m 1# 7.5m Yes m Direction Processing with wavelet transform 10m Testpit Rock subgrade Interference with No 图11测线与溶洞平面示意图 random noise Fig.11 Layout diagram of karst caves and survey line Yes Shearlet transform based 采用意大利IDS公司K2雷达探测，天线频率 onadaptive threshold value 为100MHz,时窗400ns,采样点数1024.通过人为 Displays the outcome of image processing 设置机械电磁噪声和金属干扰来模拟隧道干扰环境 图9联合法去除干扰流程图 图12(a)为常规方法处理后的波形堆积图，从 Fig.9 Flow chart of the combined methods for interference removal 图中可以明显看出240s以下深度中存在强烈同

随机干扰具有频率随机分布、波形杂乱的特点，而 频率异常干扰往往具有某些特定频率特征、在波 形堆积图上往往具有某些特定规律. 通过对正演 模拟和现场实际数据的处理可以发现，剪切变换 与小波变换对干扰的压制都有着各自的优势：剪 切变换对信号能量敏感，对于随机噪声、机械噪声 以及电磁信号等频率随机、能量异常的干扰压制 效果较好，而小波变换对频率异常、能量相近的干 扰压制效果较好. 结合两种方法各自的优势，进一步提出小波 变换与剪切变换相结合去除干扰方法，即先用小 波变换对异常频率信号进行分离，再使用剪切变 换对随机干扰进行压制. 具体实现流程图如图 9 所示. 5    实际工程案例 现场试验场地位于广西壮族自治区融水县至 河池市在建高速公路罗城段，探测对象为路基边 缘三处形态不规则裸露的溶洞，其中 1#为软塑黏 土夹碎石充填型溶洞、2#为干土充填型溶洞、3#为 空腔型溶洞. 路基为微风化灰岩、岩质坚硬、结构 面较发育，周边未见地表水. 现场情况、测线与溶 洞平面示意图分别如图 10 和图 11 所示，试验坑 长 11 m，深 1.5 m，实际测线长度 10 m. 采用意大利 IDS 公司 K2 雷达探测，天线频率 为 100 MHz，时窗 400 ns，采样点数 1024. 通过人为 设置机械电磁噪声和金属干扰来模拟隧道干扰环境. 图 12（a）为常规方法处理后的波形堆积图，从 图中可以明显看出 240 ns 以下深度中存在强烈同 Time/ns 400 600 200 50 0 Frequency/MHz 0 200 (a) 150 100 Time/ns 400 600 200 50 0 Frequency/MHz 0 200 (b) 150 100 Time/ns 400 600 200 50 0 Frequency/MHz 0 200 (c) 150 100 图 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 Data processed by conventional methods Comprehensive analysis based on on-site detection environment, waveform image characteristics, and spectral characteristics Preliminary estimation of interference types Interference with abnormal frequency Processing with wavelet transform Yes No Interference with random noise Yes Shearlet transform based onadaptive threshold value No Displays the outcome of image processing 图 9    联合法去除干扰流程图 Fig.9    Flow chart of the combined methods for interference removal 1# 2# 3# 图 10    现场情况 Fig.10    Field conditions Slope Filled with dry soil Direction 1# 2# Cavity 3# 7.5 m 12.5 m 2.5 m 10 m 11 m Testpit Rock subgrade Filled with soft plastic clay 图 11    测线与溶洞平面示意图 Fig.11    Layout diagram of karst caves and survey line · 396 · 工程科学学报，第 42 卷，第 3 期

刘宗辉等：隧道地质预报探地雷达信号干扰消除方法 397· 相轴异常干扰，3#溶洞反射波信号被淹没于干扰 分.图13(a)为测线3m处单道波时频分布，从图 信号中，根据信号特征可判定该区域同时含有随 中可清晰看出干扰信号频率成分.因此，为凸显异 机干扰和频率异常成分干扰.此外，从图中可以看 常体空间位置以及后续进一步开展属性分析，有 出240ns以上的浅部数据也存在波形杂乱无章的 必要采用本文所提出的联合去噪方法提取异常体 随机噪声干扰，左侧1#和中间2#溶洞位置难以区 反射信号 0 0 (a (b) (c) 50 100 100 100 150 150 200 200 200 250 250 250 300 300 300 350 350 350 400 4 10 400 4 400 0 2 6 8 10 4 6 10 Position/m Position/m Position/m 图12去噪效果对比.(a)常规方法：(b)小波变换：(c)联合算法 Fig.12 Comparison of different denoizing methods:(a)conventional method;(b)wavelet transform;(c)joint algorithm 300(a) 300(b) 300(c) ZHW/ 200 200 ZHW/A 200 100 0 50 0 0 100200 300 400 100200300 400 100200300400 Iime/ns me/ns Time/ns 图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)分别为使用小波变换以及联 分布图中同样可进一步清晰看出干扰信号频率成 合法处理后的波形堆积图，对处理后的雷达数据 分去除效果.通过该案例可以进一步说明小波变 选取测线3m处单道波进行时频分析，时频分布 换与剪切变换联合法去干扰的有效性与必要性 见图13(b)和13(c).从图12(b)中可以明显看出 6结语 低频成分的强反射同相轴已得到较好地去除，但 随机干扰并未得到有效去除.而图12(c)处理结果 (1)利用剪切变换将探地雷达数据转换到剪 显示整个数据的图像质量有了很大的改善，同相 切域，可以更加细致的对信号进行多尺度多方向 轴变得清晰连续，深层以及浅层的随机噪声和低 划分，通常有效信号会根据自身特点集中在一些 频干扰信号都能很好的去除.根据处理后的波形 特定的方向上，在此基础上提出的基于自适应阀 堆积图数据可以很好地将三处地质异常进行圈 值去噪方法可对随机干扰有很好去除效果 定，并能与实际情况对应.图13(b)和13(c)时频 (2)小波变换具有良好的时频局部化性质，其

相轴异常干扰，3#溶洞反射波信号被淹没于干扰 信号中，根据信号特征可判定该区域同时含有随 机干扰和频率异常成分干扰. 此外，从图中可以看 出 240 ns 以上的浅部数据也存在波形杂乱无章的 随机噪声干扰，左侧 1#和中间 2#溶洞位置难以区 分. 图 13（a）为测线 3 m 处单道波时频分布，从图 中可清晰看出干扰信号频率成分. 因此，为凸显异 常体空间位置以及后续进一步开展属性分析，有 必要采用本文所提出的联合去噪方法提取异常体 反射信号. 图 12（b）和 12（c）分别为使用小波变换以及联 合法处理后的波形堆积图，对处理后的雷达数据 选取测线 3 m 处单道波进行时频分析，时频分布 见图 13（b）和 13（c）. 从图 12（b）中可以明显看出 低频成分的强反射同相轴已得到较好地去除，但 随机干扰并未得到有效去除. 而图 12（c）处理结果 显示整个数据的图像质量有了很大的改善，同相 轴变得清晰连续，深层以及浅层的随机噪声和低 频干扰信号都能很好的去除. 根据处理后的波形 堆积图数据可以很好地将三处地质异常进行圈 定，并能与实际情况对应. 图 13（b）和 13（c）时频 分布图中同样可进一步清晰看出干扰信号频率成 分去除效果. 通过该案例可以进一步说明小波变 换与剪切变换联合法去干扰的有效性与必要性. 6    结语 （1）利用剪切变换将探地雷达数据转换到剪 切域，可以更加细致的对信号进行多尺度多方向 划分，通常有效信号会根据自身特点集中在一些 特定的方向上，在此基础上提出的基于自适应阀 值去噪方法可对随机干扰有很好去除效果. （2）小波变换具有良好的时频局部化性质，其 Time/ns 0 100 50 150 200 250 300 400 350 Position/m 0 2 4 6 8 10 (a) Time/ns 0 100 50 150 200 250 300 400 350 Position/m 0 2 4 6 8 10 (b) Time/ns 0 100 50 150 200 250 300 400 350 Position/m 0 2 4 6 8 10 (c) 1# 2# 3# 图 12    去噪效果对比. （a）常规方法；（b）小波变换；（c）联合算法 Fig.12    Comparison of different denoizing methods: (a) conventional method; (b) wavelet transform; (c) joint algorithm Time/ns 200 300 400 100 0 Frequency/MHz 0 100 (a) 300 200 Time/ns 200 300 400 50 0 Frequency/MHz 0 100 (b) 300 200 Time/ns 200 300 400 50 0 Frequency/MHz 0 100 (c) 300 200 图 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 刘宗辉等： 隧道地质预报探地雷达信号干扰消除方法 · 397 ·

·398 工程科学学报，第42卷，第3期 可变的时频窗结构对于一维信号有较强的频率分 [8]Ouadfeul S A,Aliouane L.Multiscale analysis of noisy 3D GPR 辨率，通过多尺度分解后可将频率异常信号分离， data using the directional continuous wavelet transform /2012 从而保留有效信号.但小波变换对于随机分布在 14th International Conference on Ground Penetrating Radar (GPR).Shanghai,2012:257 整个频率域且能量异常的信号去除效果不如剪切 [9]Li C M,Wang L S,Xu M J,et al.Objects recognition of ground 变换 penetrating radar in karst regions using wavelet energy spectrum (3)隧道空间环境复杂，使用探地雷达进行隧 analysis.Chin Geophys,2006,49(5):1499 道超前地质预报时往往会遇到各种成分的干扰混 (李才明，王良书，徐鸣洁，等.基于小波能谱分析的岩溶区探地 叠，采用单一干扰去除方法难以获得满意效果.提 雷达目标识别.地球物理学报，2006,49(5)：1499) 出的剪切变换与小波变换联合方法可以同时对能 [10]Addison A D,Battista B M,Knapp C C.Improved 量异常的随机干扰以及频率异常的干扰信号进行 hydrogeophysical parameter estimation from empirical mode 压制，并且可以保证处理后的数据有着较高信噪比. decomposition processed ground penetrating radar data./Emiron (4)根据具体的干扰类型以及干扰数据特征 Eng Geophys,2009,14(4:171 [11]Chen C S,Jeng Y.Nonlinear data processing method for the signal 选择合适的阀值系数，直接关系到千扰去除和有 enhancement of GPR data.JApp/Geoplys,2011,75(1):113 效信号保留效果，自适应阀值函数仍有许多不确 [12]Zhang Z Y,Zhang X D,Yu HY,et al.Noise suppression based on 定性.实际工程应用中需要不断总结隧道中常见 a fast discrete curvelet transform.Geophys Eng,2010,7(1):105 各类干扰的剪切系数能量特征以及小波系数频率 [13]Zhu Z Q,Zhu H,Lu G Y,et al.Processing of GPR data in tunnel 特征，并形成干扰信号样本库 fissure water based on Curvelet transform.Comput Tech Geophys Geochem Explor,2014,36(5):571 参考文献 (朱自强，朱贺，鲁光银，等.基于Curvelet?变换的隧道裂隙水 [1] Xue Y G,Li S C,Su M X,et al.Study of geological prediction GPR数据处理研究.物探化探计算技术，2014,36(5)：571) implementation method in tunnel construction.Rock Soil Mech [14]Zhou L,Li S C,Xu Z H,et al.Interpretation and treatment of 2011,32(8):2416 interfering factors in advance geological prediction by ground (薛翊国，李术才，苏茂鑫，等.隧道施工期超前地质预报实施方 penetrating radar of tunnel construction.Tunnel Constr,2016. 法研究.岩土力学，2011,32(8)：2416) 36(12):1517 [2]Gao Y T,Xu J,Wu S C,et al.An intelligent identification method (周轮，李术才，许振浩，等.隧道施工期超前预报地质雷达异常 to detect tunnel defects based on the multidimensional analysis of 千扰识别及处理.隧道建设，2016,36(12)：1517) GPR reflections.Chin J Eng,2018,40(3):293 [15]Liu C M.Wang DL,Wang T,et al.Random seismic noise (高永涛，徐俊，吴顺川，等.基于GR反射波信号多维分析的隧 attention based on the Shearlet transform.Acta Petrol Sin,2014 道病害智能辨识.工程科学学报，2018,40(3)：293) 35(4):692 [3]Jol H M.Ground Penetrating Radar:Theory and Applications (刘成明，王德利，王通，等.基于Shearlet?变换的地震随机噪声压 Elsevier,2009 制.石油学报，2014,35(4)：692) [4]Yang G,Liu D W.Frequency spectrum characteristic of quartz [16]Liu C M,Wang D L,Hu B,et al.Seismic date interpolation based sandstone rock mass with GPR.ChinJ Eng,2015,37(11):1397 on sparse in Shearlet domain.J Jilin Univ Earth Sci Ed,2016. (杨光，刘敦文.石英砂岩体的地质雷达波频谱特征.工程科学 46(6):1855 学报，2015,37(11)：1397) (刘成明，王德利，胡斌，等.Shearlet域稀疏约束地震数据重建 [5]Liu G,Li S C,Xue Y G,et al.GPR signal processing approach 吉林大学学报：地球科学版，2016,46(6)：1855) under low signal to noise ratio based on wavelet transforms and its [17]Guo K,Labate D.The construction of smooth Parseval frames of application.Geotech Imvest Surv,2009,37(9):85 shearlets.Math Modell Nat Phenom,2013,8(1):82 (柳刚，李术才，薛翊国，等.基于小波变换的雷达低信噪比信号 [18]Easley G,Labate D,Lim W Q.Sparse directional image 处理技术及应用研究.工程勘察，2009,37(9)：85) representations using the discrete shearlet transform.App/Comput [6]Bao QZ,Li Q C,Chen W C.GPR data noise attenuation on the Harmon Anal,2008,25(1):25 curvelet transform.Appl Geoplys,2014,11(3):301 [19]Zhang X W,Gao Y Z,Fang G Y,et al.Application of generalized [7]Gan L,Zhou L,You X G,et al.The instantaneous frequency S transform with low-pass filtering to layer recognition of Ground extraction of GPR B-scan data based on HHT method /2012 Penetrating Radar.Chin J Geoplrys,2013,56(1):309 International Conference on Machine Learning and Cybernetics (张先武，高云泽，方广有.带有低通滤波的广义$变换在探地雷 (CMLC).Xi'an,2012:982 达层位识别中的应用.地球物理学报，2013,56(1)：309)

可变的时频窗结构对于一维信号有较强的频率分 辨率，通过多尺度分解后可将频率异常信号分离， 从而保留有效信号. 但小波变换对于随机分布在 整个频率域且能量异常的信号去除效果不如剪切 变换. （3）隧道空间环境复杂，使用探地雷达进行隧 道超前地质预报时往往会遇到各种成分的干扰混 叠，采用单一干扰去除方法难以获得满意效果. 提 出的剪切变换与小波变换联合方法可以同时对能 量异常的随机干扰以及频率异常的干扰信号进行 压制，并且可以保证处理后的数据有着较高信噪比. （4）根据具体的干扰类型以及干扰数据特征 选择合适的阀值系数，直接关系到干扰去除和有 效信号保留效果，自适应阀值函数仍有许多不确 定性. 实际工程应用中需要不断总结隧道中常见 各类干扰的剪切系数能量特征以及小波系数频率 特征，并形成干扰信号样本库. 参    考    文    献 Xue  Y  G,  Li  S  C,  Su  M  X,  et  al.  Study  of  geological  prediction implementation  method  in  tunnel  construction. Rock Soil Mech, 2011, 32（8）: 2416 （薛翊国, 李术才, 苏茂鑫, 等. 隧道施工期超前地质预报实施方 法研究. 岩土力学, 2011, 32（8）：2416） [1] Gao Y T, Xu J, Wu S C, et al. An intelligent identification method to detect tunnel defects based on the multidimensional analysis of GPR reflections. Chin J Eng, 2018, 40（3）: 293 （高永涛, 徐俊, 吴顺川, 等. 基于GPR反射波信号多维分析的隧 道病害智能辨识. 工程科学学报, 2018, 40（3）：293） [2] Jol  H  M. Ground Penetrating Radar: Theory and Applications. Elsevier, 2009 [3] Yang  G,  Liu  D  W.  Frequency  spectrum  characteristic  of  quartz sandstone rock mass with GPR. Chin J Eng, 2015, 37（11）: 1397 （杨光, 刘敦文. 石英砂岩体的地质雷达波频谱特征. 工程科学 学报, 2015, 37（11）：1397） [4] Liu  G,  Li  S  C,  Xue  Y  G,  et  al.  GPR  signal  processing  approach under low signal to noise ratio based on wavelet transforms and its application. Geotech Invest Surv, 2009, 37（9）: 85 （柳刚, 李术才, 薛翊国, 等. 基于小波变换的雷达低信噪比信号 处理技术及应用研究. 工程勘察, 2009, 37（9）：85） [5] Bao Q Z, Li Q C, Chen W C. GPR data noise attenuation on the curvelet transform. Appl Geophys, 2014, 11（3）: 301 [6] Gan  L,  Zhou  L,  You  X  G,  et  al.  The  instantaneous  frequency extraction  of  GPR  B-scan  data  based  on  HHT  method  //  2012 International Conference on Machine Learning and Cybernetics (ICMLC). Xi’an, 2012: 982 [7] Ouadfeul S A, Aliouane L. Multiscale analysis of noisy 3D GPR data  using  the  directional  continuous  wavelet  transform  //  2012 14th International Conference on Ground Penetrating Radar (GPR). Shanghai, 2012: 257 [8] Li C M, Wang L S, Xu M J, et al. Objects recognition of ground penetrating  radar  in  karst  regions  using  wavelet  energy  spectrum analysis. Chin J Geophys, 2006, 49（5）: 1499 （李才明, 王良书, 徐鸣洁, 等. 基于小波能谱分析的岩溶区探地 雷达目标识别. 地球物理学报, 2006, 49（5）：1499） [9] Addison  A  D,  Battista  B  M,  Knapp  C  C.  Improved hydrogeophysical  parameter  estimation  from  empirical  mode decomposition processed ground penetrating radar data. J Environ Eng Geophys, 2009, 14（4）: 171 [10] Chen C S, Jeng Y. Nonlinear data processing method for the signal enhancement of GPR data. J Appl Geophys, 2011, 75（1）: 113 [11] Zhang Z Y, Zhang X D, Yu H Y, et al. Noise suppression based on a fast discrete curvelet transform. J Geophys Eng, 2010, 7（1）: 105 [12] Zhu Z Q, Zhu H, Lu G Y, et al. Processing of GPR data in tunnel fissure water based on Curvelet transform. Comput Tech Geophys Geochem Explor, 2014, 36（5）: 571 （朱自强, 朱贺, 鲁光银, 等. 基于Curvelet变换的隧道裂隙水 GPR数据处理研究. 物探化探计算技术, 2014, 36（5）：571） [13] Zhou  L,  Li  S  C,  Xu  Z  H,  et  al.  Interpretation  and  treatment  of interfering  factors  in  advance  geological  prediction  by  ground penetrating  radar  of  tunnel  construction. Tunnel Constr,  2016, 36（12）: 1517 （周轮, 李术才, 许振浩, 等. 隧道施工期超前预报地质雷达异常 干扰识别及处理. 隧道建设, 2016, 36（12）：1517） [14] Liu  C  M,  Wang  D  L,  Wang  T,  et  al.  Random  seismic  noise attention  based  on  the  Shearlet  transform. Acta Petrol Sin,  2014, 35（4）: 692 （刘成明, 王德利, 王通, 等. 基于Shearlet变换的地震随机噪声压 制. 石油学报, 2014, 35（4）：692） [15] Liu C M, Wang D L, Hu B, et al. Seismic date interpolation based on  sparse  in  Shearlet  domain. J Jilin Univ Earth Sci Ed,  2016, 46（6）: 1855 （刘成明, 王德利, 胡斌, 等. Shearlet域稀疏约束地震数据重建. 吉林大学学报: 地球科学版, 2016, 46（6）：1855） [16] Guo K, Labate D. The construction of smooth Parseval frames of shearlets. Math Modell Nat Phenom, 2013, 8（1）: 82 [17] Easley  G,  Labate  D,  Lim  W  Q.  Sparse  directional  image representations using the discrete shearlet transform. Appl Comput Harmon Anal, 2008, 25（1）: 25 [18] Zhang X W, Gao Y Z, Fang G Y, et al. Application of generalized S transform with low-pass filtering to layer recognition of Ground Penetrating Radar. Chin J Geophys, 2013, 56（1）: 309 （张先武, 高云泽, 方广有. 带有低通滤波的广义S变换在探地雷 达层位识别中的应用. 地球物理学报, 2013, 56（1）：309） [19] · 398 · 工程科学学报，第 42 卷，第 3 期