双树复小波变换在川藏铁路拉林段某隧道超前地质预报中的应用

杨丹，李伟，魏永梁，宋斌

(1.西南交通大学 地球科学与环境工程学院,四川 成都 611756；2.中国中铁科学研究院有限公司,四川 成都 610032)

0 引言

为探明掌子面前方地质条件，查清该隧道施工掌子面前方断层、节理裂隙分布等工程地质情况，保证施工安全，采用探地雷达进行超前地质预报。但由于掌子面不平整、隧道所在地层地质条件极其复杂，探测条件较差，隧道内施工环境复杂等原因，导致采集的数据中存在强干扰覆盖有效信号的现象，严重影响数据采集质量，干扰工程人员的判断，因此需要进行信号处理技术研究，以改善数据质量，增强有效回波信号，消除杂波及多次回波干扰的影响，提高超前地质预报的准确性。

传统的探地雷达(GPR)去噪方法以傅里叶变换理论为基础，通过数学计算实现滤波，能够在干扰信号和有效信号的频谱范围存在明显分界时有效滤除指定频率(干扰波频谱)的信号。但探地雷达信号属于非平稳信号，用传统的滤波方法去噪的同时也会拓宽波形，损失重要的信号信息，且无法同时体现信号频率与时域的信息，难以满足复杂地质条件和复杂工程环境下的信号处理。近年来，专家、学者将小波变换理论引入探地雷达信号处理中，此举在继承和发展了短时傅里叶变换的局部化思想的同时，又克服了窗口大小不随频率变换的缺点，可利用小波变换的调焦功能和频率域—时间域双重局部性来压制噪声[2-6]。多数学者选用离散小波变换处理探地雷达信号数据[7-10]，但离散小波变换存在平移敏感性、频率混叠和方向性选择少的缺点，会影响数据处理的准确性，而双树复小波变换(dual-tree complex wavelet transform，DTCWT)承袭了传统离散小波变换优越的时频分析和多分辨率分析能力，同时又克服了传统离散小波变换的平移敏感性、频率混叠和缺少方向选择性等不足[11]。目前，国内外学者利用双树复小波对探地雷达信号进行处理的研究还相对较少。Behrooz Oskooi等人在研究利用双树复小波对探地雷达信号进行去噪处理时，在阈值的选择上仅考虑了通用阈值[12]，而通用阈值存在“过扼杀”小波系数的情况，会丢失有用的高频信号[13]。

本文以川藏铁路拉林段某隧道为例，在对探测信号实施双树复小波变换的基础上，采用通用阈值、无偏似然估计阈值、启发式阈值、极大极小阈值4种阈值选择方法，结合软阈值、硬阈值、半软阈值、非负消减4种阈值处理方法，对小波高频系数进行阈值量化处理后，进行双树复小波逆变换重构信号，进而达到对探地雷达信号进行去噪的目的，以提高超前地质预报识别精度。

1 掌子面测线布置

根据前期地质勘探资料，结合隧道施工需求，在隧道掌子面上布置2条测线(图1)，采用SIR-20型探地雷达和100 MHz天线,对掌子面前方30 m范围的断层、节理裂隙发育情况等进行探测。为提高探地雷达探测信号质量，利用双树复小波变换进行4层分解，分别采用4种不同的阈值选择方法并结合4种处理方法对小波系数进行阈值处理，用双树复小波逆变换进行重构，并对重构后的信号进行去噪效果评价。

图1 探地雷达测线布置示意

2 基于双树复小波的信号处理方法

2.1 双树复小波变换

探地雷达向被探测的介质发射宽带电磁波脉冲信号，其中一部分电磁波在遇到不同介电常数的介质界面时被反射回来，而这些反射回来的电磁波信号就携带了大量的信息，工程人员通过分析这些信息，解读、判识被探测物体内部信息；但在探测过程中存在大量环境干扰，这些回波信息中存在各种噪声信息，在后期数据处理时需要压制和去除干扰信号，解读有效信号。文中采用双树复小波方法进行信号数据变换，利用双树结构、采用满足希尔伯特变换的滤波器进行探地雷达隧道超前地质预报信号处理与降噪。

探地雷达接收的雷达信号x(t)可以分解为：

x(t)=s(t)+n(t)，

(1)

双树复小波变换对信号x(t)进行两路离散小波变换，即采用2个独立的小波来计算信号的实部与虚部。如图2所示，树a为实部，树b为虚部，由实部和虚部2个树状的共轭正交滤波器组构成一个复小波，h0(n)和h1(n)分别对应实部的低通滤波器和高通滤波器，g0(n)和g1(n)分别对应虚部的低通滤波器和高通滤波器，且满足：

图2 双树复小波结构示意

G0(ejw)≈H0(ejw)e-j0.5w，

(2)

式中：G0(ejw)和H0(ejw)分别是g0(n)和h0(n)的傅里叶变换[12]。双树复小波的逆变换如图3所示。

图3 双树复小波逆变换

2.2 小波阈值去噪

2.2.1 小波去噪信号处理过程

小波阈值去噪过程分为3个阶段：分解、阈值量化和重构，在双树复小波变换的基础上，将含噪信号进行分解，保留最底层的低频系数，仅对各个分解尺度下的高频系数选择一个阈值进行量化处理，将处理过的小波系数通过双树复小波逆变换进行重构，进而获得去除噪声后的信号(图4)。这个过程中最关键的是如何选择阈值，并利用该阈值对小波系数进行量化处理，它关系到信号消噪的质量。

图4 双树复小波去噪流程示意

2.2.2 阈值选择及阈值处理

阈值选择直接关系着去噪结果，过大的阈值导致过多小波系数萎缩，去噪后信号过于平滑；而过小的阈值能最大限度保留信号的形状，却无法抑制噪声。各阈值选择方法各有优缺点，本文采用通用阈值(universal threshold)、无偏似然估计阈值(rigrsure)、启发式阈值(heursure)、极大极小阈值(minimaxi)[3]4种方法进行阈值选择。

1)通用阈值

Donoho在高斯噪声模型下，应用多维独立正态变量决策理论，得到了通用的阈值：

(3)

式中：σ为噪声的标准差，N为信号长度。

2)无偏似然估计阈值

无偏似然估计阈值是一种基于 Stein 的无偏似然估计原理的自适应阈值选择，是一种软件阈值估计器。

3)启发式阈值

启发式阈值是前两种阈值的综合，所选择的是最优预测变量阈值，视信号情况选定通用阈值或无偏似然估计阈值。

4)极大极小阈值

极值阈值也是一种固定阈值选择形式，它所产生的是一个最小均方差的极值，不是无误差。主要应用极大极小原理选择阈值[3]。

(4)

式中:wi为尺度i上各点的小波系数。当噪声为非白噪声时，可以在每个分解尺度上都估计噪声的方差，从而计算不同分解尺度上的阈值。由于探地雷达噪声受环境干扰较大，存在非白噪声干扰，因此对不同分解层次的小波系数采用不同的阈值，以达到优化去噪效果。

在把信号进行分解后，要将得到的高频小波系数w(i)与阈值T比较。为了获得准确的数据，在采用软、硬阈值处理方法的同时，又引入了半软阈值处理(firm shrinkage)[14]、非负消减处理(nonnegative garrote shrinkage)[15]2种改进的阈值处理方法。

1)硬阈值处理

(5)

硬阈值处理将小波系数中绝对值小于阈值T的小波系数置零，因此会产生不连续的情况。

2)软阈值处理

(6)

软阈值处理后的小波系数虽然整体连续性好，但是对所有的小波系数都进行了抑制，造成高频信息损失，因此给重构带了不小的误差[13]。

5) 将 gmchcfg.h 文件放置在C:Tornado221 argethugldrivergraphicsintel 路径下(本文假设. VxWorks 开发环境安装在 C:Tornado221 路径下);

3)半软阈值处理

(7)

Gao和Bruce于1997年提出半软阈值处理[14]，其中T2=2T。

4)非负消减处理

(8)

Breiman于1995年提出了非负消减处理[15]，这种方法在软、硬阈值处理之间做了折中。

2.3 去噪效果评价

为评价双树复小波变换去噪重构后的性能，采用以下3个参数进行评价。

1)归一化均方根误差NRMSE：

(9)

2)信噪比SNR：

(10)

3)能量比P：

(11)

式中：x(i)是探地雷达采集到的原始信号，x′(i)为双树复小波去噪重构后的信号，ux为原始信号x(i)的平均值。归一化均方根误差描述了去噪信号x′(i)与原始信号x(i)之间的相似程度[5]，信噪比指的是有效波与干扰波相对强弱的比较，能量比为降噪后的信号在原信号所占的能量比例；归一化均方根误差越小，信噪比越大，能量比越大，说明信号失真越小，信号去噪效果越好[4]。

本文对探地雷达数据进行双树复小波变换去噪，采用4种阈值选择和4种阈值处理方法两两组合，通过对去噪效果的评价，并结合去噪前后波形图及堆积图的对比，找到最佳的阈值选择和阈值处理方法的组合。

3 实测数据处理效果

3.1 1#掌子面应用情况

在隧道1#掌子面上水平布置2条测线，其中测线2数据共计415道，单道数据采样点数1 024，记录长度600 ns。对测线2数据进行双树复小波变换去噪并重构，以第250道数据为例进行评价分析。

从表1可以看出: 去噪后的归一化均方差NRMSE值都比较小，在0.1左右，但采用无偏似然估计或启发式阈值选择方法，结合非负消减阈值处理方法时，得到的归一化均方差值最小，仅为0.04;去噪后的SNR信号信噪比较高，均在75～86.7 dB之间，但采用无偏似然估计或启发式阈值选择方法，结合非负消减阈值处理方法时，得到的信噪比值最大；去噪后的信号P值信号能量比相差较大，采用无偏似然估计阈值选择方法，结合非负消减阈值处理方法时，得到的信号能量比最大，达到99.06%。图5为第250道采用无偏似然估计阈值选择方法和非负消减阈值处理方法重构的信号波形，图中显示多次波及大部分噪声被明显去除。

图5 第250道单道原始信号与双树复小波变换去噪信号对比

表1 第250道数据去噪效果评价

通过归一化均方根误差、信噪比、能量比的结果分析，采用无偏似然估计阈值选择方法和非负消减阈值处理方法对1#掌子面信号进行双树复小波变换后的去噪效果最佳。

从整个掌子面探测剖面信号的处理情况看(图6)，采用无偏似然估计阈值选择方法，结合非负消减阈值处理方法进行双树复小波变换去噪处理之后，噪声干扰被有效消除，消减了多次反射波，异常体的边界及位置清晰可见，比传统滤波方法去噪效果更好。图6中可以清晰看到250～350道之间的强反射区波形存在明显畸变，同相轴连续性较差，结合地质勘察资料与掌子面情况，可以推断该段围岩较破碎，节理裂隙较发育，稳定性较差，施工时应加强支护。

图6 1#掌子面传统滤波与双树复小波去噪结果对比

3.2 2#掌子面应用情况

在隧道2#掌子面上水平布置2条测线，其中测线1数据共计800道，单道数据采样点数1 024，记录长度600 ns。对测线1数据进行双树复小波变换去噪并重构，以第700道数据重构信号为例进行去噪效果评价。

从表2可以看出:采用无偏似然估计或启发式阈值选择方法，结合非负消减阈值处理方法得到的归一化均方差NRMSE值最小，还不到0.03，达到了较好的去噪效果,而采用其他方法时，归一化均方差NRMSE值相对较大，去噪效果不理想；采用无偏似然估计阈值或启发式阈值选择方法，结合非负消减阈值处理方法时，去噪后的信噪比较大，可达到60 dB以上，信号能量远大于噪声能量，去噪效果较好，而其他方法去噪后的信噪比只有27～30 dB，去噪效果一般；采用无偏似然估计阈值选择方法，结合非负消减阈值处理方法时，能量比值P最大，达到98.88%，去噪效果较好，而采用其他方法得到的能量比值P相对较小，去噪效果不理想。采用无偏似然估计阈值选择方法和非负消减阈值处理方法重构信号的波形图如图7所示，图中显示多次波及大部分噪声被明显去除。

图7 第700道原始信号与双树复小波变换去噪信号对比

表2 第700道数据去噪效果评价

通过归一化均方根误差、信噪比、能量比的结果分析，采用无偏似然估计阈值选择方法和非负消减阈值处理方法对2#掌子面信号进行双树复小波变换后的去噪效果最佳。从整个掌子面探测剖面信号的处理情况看(图8)，采用无偏似然估计阈值选择方法并结合非负消减阈值处理方法进行双树复小波变换去噪处理之后，噪声干扰被有效消除，消减了多次反射波，异常体的边界及位置清晰可见，比传统滤波方法去噪效果更好。从图8中可以清晰看到600～750道之间的强反射区波形存在明显畸变，同相轴连续性较差；结合地质勘察资料与掌子面情况，可以推断该段围岩较破碎，节理裂隙较发育，稳定性较差，施工时应加强支护。

图8 2#掌子面传统滤波与双树复小波去噪结果对比

4 结论

受川藏铁路极其复杂的地质条件和施工环境的影响，探地雷达地质超前预报技术采集的信号噪声干扰较大，给信号处理与数据分析带来较大困难。在对探测信号实施双树复小波变换的基础上，采用通用阈值、无偏似然估计阈值、启发式阈值、极大极小阈值4种阈值选择方法，结合软阈值、硬阈值、半软阈值、非负消减4种阈值处理方法，对小波高频系数进行阈值量化处理后进行双树复小波逆变换重构信号，达到了较好的去噪效果。对比分析归一化均方差、信噪比、能量比3个评价因子后发现，采用无偏似然估计阈值选择方法，结合非负消减阈值处理方法，可以获得最佳的去噪效果，信号得到很大的改善，噪声残留较少，多次反射波被减弱，提高了超前地质预报中异常体的识别精度，能够为川藏铁路隧道施工提供更好的安全保障。