楊丹，李偉，魏永梁，宋斌

(1.西南交通大學 地球科學與環(huán)境工程學院,四川 成都 611756；2.中國中鐵科學研究院有限公司,四川 成都 610032)

0 引言

為探明掌子面前方地質條件，查清該隧道施工掌子面前方斷層、節(jié)理裂隙分布等工程地質情況，保證施工安全，采用探地雷達進行超前地質預報。但由于掌子面不平整、隧道所在地層地質條件極其復雜，探測條件較差，隧道內施工環(huán)境復雜等原因，導致采集的數(shù)據(jù)中存在強干擾覆蓋有效信號的現(xiàn)象，嚴重影響數(shù)據(jù)采集質量，干擾工程人員的判斷，因此需要進行信號處理技術研究，以改善數(shù)據(jù)質量，增強有效回波信號，消除雜波及多次回波干擾的影響，提高超前地質預報的準確性。

傳統(tǒng)的探地雷達(GPR)去噪方法以傅里葉變換理論為基礎，通過數(shù)學計算實現(xiàn)濾波，能夠在干擾信號和有效信號的頻譜范圍存在明顯分界時有效濾除指定頻率(干擾波頻譜)的信號。但探地雷達信號屬于非平穩(wěn)信號，用傳統(tǒng)的濾波方法去噪的同時也會拓寬波形，損失重要的信號信息，且無法同時體現(xiàn)信號頻率與時域的信息，難以滿足復雜地質條件和復雜工程環(huán)境下的信號處理。近年來，專家、學者將小波變換理論引入探地雷達信號處理中，此舉在繼承和發(fā)展了短時傅里葉變換的局部化思想的同時，又克服了窗口大小不隨頻率變換的缺點，可利用小波變換的調焦功能和頻率域—時間域雙重局部性來壓制噪聲[2-6]。多數(shù)學者選用離散小波變換處理探地雷達信號數(shù)據(jù)[7-10]，但離散小波變換存在平移敏感性、頻率混疊和方向性選擇少的缺點，會影響數(shù)據(jù)處理的準確性，而雙樹復小波變換(dual-tree complex wavelet transform，DTCWT)承襲了傳統(tǒng)離散小波變換優(yōu)越的時頻分析和多分辨率分析能力，同時又克服了傳統(tǒng)離散小波變換的平移敏感性、頻率混疊和缺少方向選擇性等不足[11]。目前，國內外學者利用雙樹復小波對探地雷達信號進行處理的研究還相對較少。Behrooz Oskooi等人在研究利用雙樹復小波對探地雷達信號進行去噪處理時，在閾值的選擇上僅考慮了通用閾值[12]，而通用閾值存在“過扼殺”小波系數(shù)的情況，會丟失有用的高頻信號[13]。

本文以川藏鐵路拉林段某隧道為例，在對探測信號實施雙樹復小波變換的基礎上，采用通用閾值、無偏似然估計閾值、啟發(fā)式閾值、極大極小閾值4種閾值選擇方法，結合軟閾值、硬閾值、半軟閾值、非負消減4種閾值處理方法，對小波高頻系數(shù)進行閾值量化處理后，進行雙樹復小波逆變換重構信號，進而達到對探地雷達信號進行去噪的目的，以提高超前地質預報識別精度。

1 掌子面測線布置

根據(jù)前期地質勘探資料，結合隧道施工需求，在隧道掌子面上布置2條測線(圖1)，采用SIR-20型探地雷達和100 MHz天線,對掌子面前方30 m范圍的斷層、節(jié)理裂隙發(fā)育情況等進行探測。為提高探地雷達探測信號質量，利用雙樹復小波變換進行4層分解，分別采用4種不同的閾值選擇方法并結合4種處理方法對小波系數(shù)進行閾值處理，用雙樹復小波逆變換進行重構，并對重構后的信號進行去噪效果評價。

圖1 探地雷達測線布置示意

2 基于雙樹復小波的信號處理方法

2.1 雙樹復小波變換

探地雷達向被探測的介質發(fā)射寬帶電磁波脈沖信號，其中一部分電磁波在遇到不同介電常數(shù)的介質界面時被反射回來，而這些反射回來的電磁波信號就攜帶了大量的信息，工程人員通過分析這些信息，解讀、判識被探測物體內部信息；但在探測過程中存在大量環(huán)境干擾，這些回波信息中存在各種噪聲信息，在后期數(shù)據(jù)處理時需要壓制和去除干擾信號，解讀有效信號。文中采用雙樹復小波方法進行信號數(shù)據(jù)變換，利用雙樹結構、采用滿足希爾伯特變換的濾波器進行探地雷達隧道超前地質預報信號處理與降噪。

探地雷達接收的雷達信號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個階段：分解、閾值量化和重構，在雙樹復小波變換的基礎上，將含噪信號進行分解，保留最底層的低頻系數(shù)，僅對各個分解尺度下的高頻系數(shù)選擇一個閾值進行量化處理，將處理過的小波系數(shù)通過雙樹復小波逆變換進行重構，進而獲得去除噪聲后的信號(圖4)。這個過程中最關鍵的是如何選擇閾值，并利用該閾值對小波系數(shù)進行量化處理，它關系到信號消噪的質量。

圖4 雙樹復小波去噪流程示意

2.2.2 閾值選擇及閾值處理

閾值選擇直接關系著去噪結果，過大的閾值導致過多小波系數(shù)萎縮，去噪后信號過于平滑；而過小的閾值能最大限度保留信號的形狀，卻無法抑制噪聲。各閾值選擇方法各有優(yōu)缺點，本文采用通用閾值(universal threshold)、無偏似然估計閾值(rigrsure)、啟發(fā)式閾值(heursure)、極大極小閾值(minimaxi)[3]4種方法進行閾值選擇。

1)通用閾值

Donoho在高斯噪聲模型下，應用多維獨立正態(tài)變量決策理論，得到了通用的閾值：

(3)

式中：σ為噪聲的標準差，N為信號長度。

2)無偏似然估計閾值

無偏似然估計閾值是一種基于 Stein 的無偏似然估計原理的自適應閾值選擇，是一種軟件閾值估計器。

3)啟發(fā)式閾值

啟發(fā)式閾值是前兩種閾值的綜合，所選擇的是最優(yōu)預測變量閾值，視信號情況選定通用閾值或無偏似然估計閾值。

4)極大極小閾值

極值閾值也是一種固定閾值選擇形式，它所產生的是一個最小均方差的極值，不是無誤差。主要應用極大極小原理選擇閾值[3]。

(4)

式中:wi為尺度i上各點的小波系數(shù)。當噪聲為非白噪聲時，可以在每個分解尺度上都估計噪聲的方差，從而計算不同分解尺度上的閾值。由于探地雷達噪聲受環(huán)境干擾較大，存在非白噪聲干擾，因此對不同分解層次的小波系數(shù)采用不同的閾值，以達到優(yōu)化去噪效果。

在把信號進行分解后，要將得到的高頻小波系數(shù)w(i)與閾值T比較。為了獲得準確的數(shù)據(jù)，在采用軟、硬閾值處理方法的同時，又引入了半軟閾值處理(firm shrinkage)[14]、非負消減處理(nonnegative garrote shrinkage)[15]2種改進的閾值處理方法。

1)硬閾值處理

(5)

硬閾值處理將小波系數(shù)中絕對值小于閾值T的小波系數(shù)置零，因此會產生不連續(xù)的情況。

2)軟閾值處理

(6)

軟閾值處理后的小波系數(shù)雖然整體連續(xù)性好，但是對所有的小波系數(shù)都進行了抑制，造成高頻信息損失，因此給重構帶了不小的誤差[13]。

5) 將 gmchcfg.h 文件放置在C:Tornado221 argethugldrivergraphicsintel 路徑下(本文假設. VxWorks 開發(fā)環(huán)境安裝在 C:Tornado221 路徑下);

3)半軟閾值處理

(7)

Gao和Bruce于1997年提出半軟閾值處理[14]，其中T2=2T。

4)非負消減處理

(8)

Breiman于1995年提出了非負消減處理[15]，這種方法在軟、硬閾值處理之間做了折中。

2.3 去噪效果評價

為評價雙樹復小波變換去噪重構后的性能，采用以下3個參數(shù)進行評價。

1)歸一化均方根誤差NRMSE：

(9)

2)信噪比SNR：

(10)

3)能量比P：

(11)

式中：x(i)是探地雷達采集到的原始信號，x′(i)為雙樹復小波去噪重構后的信號，ux為原始信號x(i)的平均值。歸一化均方根誤差描述了去噪信號x′(i)與原始信號x(i)之間的相似程度[5]，信噪比指的是有效波與干擾波相對強弱的比較，能量比為降噪后的信號在原信號所占的能量比例；歸一化均方根誤差越小，信噪比越大，能量比越大，說明信號失真越小，信號去噪效果越好[4]。

本文對探地雷達數(shù)據(jù)進行雙樹復小波變換去噪，采用4種閾值選擇和4種閾值處理方法兩兩組合，通過對去噪效果的評價，并結合去噪前后波形圖及堆積圖的對比，找到最佳的閾值選擇和閾值處理方法的組合。

3 實測數(shù)據(jù)處理效果

3.1 1#掌子面應用情況

在隧道1#掌子面上水平布置2條測線，其中測線2數(shù)據(jù)共計415道，單道數(shù)據(jù)采樣點數(shù)1 024，記錄長度600 ns。對測線2數(shù)據(jù)進行雙樹復小波變換去噪并重構，以第250道數(shù)據(jù)為例進行評價分析。

從表1可以看出: 去噪后的歸一化均方差NRMSE值都比較小，在0.1左右，但采用無偏似然估計或啟發(fā)式閾值選擇方法，結合非負消減閾值處理方法時，得到的歸一化均方差值最小，僅為0.04;去噪后的SNR信號信噪比較高，均在75～86.7 dB之間，但采用無偏似然估計或啟發(fā)式閾值選擇方法，結合非負消減閾值處理方法時，得到的信噪比值最大；去噪后的信號P值信號能量比相差較大，采用無偏似然估計閾值選擇方法，結合非負消減閾值處理方法時，得到的信號能量比最大，達到99.06%。圖5為第250道采用無偏似然估計閾值選擇方法和非負消減閾值處理方法重構的信號波形，圖中顯示多次波及大部分噪聲被明顯去除。

圖5 第250道單道原始信號與雙樹復小波變換去噪信號對比

表1 第250道數(shù)據(jù)去噪效果評價

通過歸一化均方根誤差、信噪比、能量比的結果分析，采用無偏似然估計閾值選擇方法和非負消減閾值處理方法對1#掌子面信號進行雙樹復小波變換后的去噪效果最佳。

從整個掌子面探測剖面信號的處理情況看(圖6)，采用無偏似然估計閾值選擇方法，結合非負消減閾值處理方法進行雙樹復小波變換去噪處理之后，噪聲干擾被有效消除，消減了多次反射波，異常體的邊界及位置清晰可見，比傳統(tǒng)濾波方法去噪效果更好。圖6中可以清晰看到250～350道之間的強反射區(qū)波形存在明顯畸變，同相軸連續(xù)性較差，結合地質勘察資料與掌子面情況，可以推斷該段圍巖較破碎，節(jié)理裂隙較發(fā)育，穩(wěn)定性較差，施工時應加強支護。

圖6 1#掌子面?zhèn)鹘y(tǒng)濾波與雙樹復小波去噪結果對比

3.2 2#掌子面應用情況

在隧道2#掌子面上水平布置2條測線，其中測線1數(shù)據(jù)共計800道，單道數(shù)據(jù)采樣點數(shù)1 024，記錄長度600 ns。對測線1數(shù)據(jù)進行雙樹復小波變換去噪并重構，以第700道數(shù)據(jù)重構信號為例進行去噪效果評價。

從表2可以看出:采用無偏似然估計或啟發(fā)式閾值選擇方法，結合非負消減閾值處理方法得到的歸一化均方差NRMSE值最小，還不到0.03，達到了較好的去噪效果,而采用其他方法時，歸一化均方差NRMSE值相對較大，去噪效果不理想；采用無偏似然估計閾值或啟發(fā)式閾值選擇方法，結合非負消減閾值處理方法時，去噪后的信噪比較大，可達到60 dB以上，信號能量遠大于噪聲能量，去噪效果較好，而其他方法去噪后的信噪比只有27～30 dB，去噪效果一般；采用無偏似然估計閾值選擇方法，結合非負消減閾值處理方法時，能量比值P最大，達到98.88%，去噪效果較好，而采用其他方法得到的能量比值P相對較小，去噪效果不理想。采用無偏似然估計閾值選擇方法和非負消減閾值處理方法重構信號的波形圖如圖7所示，圖中顯示多次波及大部分噪聲被明顯去除。

圖7 第700道原始信號與雙樹復小波變換去噪信號對比

表2 第700道數(shù)據(jù)去噪效果評價

通過歸一化均方根誤差、信噪比、能量比的結果分析，采用無偏似然估計閾值選擇方法和非負消減閾值處理方法對2#掌子面信號進行雙樹復小波變換后的去噪效果最佳。從整個掌子面探測剖面信號的處理情況看(圖8)，采用無偏似然估計閾值選擇方法并結合非負消減閾值處理方法進行雙樹復小波變換去噪處理之后，噪聲干擾被有效消除，消減了多次反射波，異常體的邊界及位置清晰可見，比傳統(tǒng)濾波方法去噪效果更好。從圖8中可以清晰看到600～750道之間的強反射區(qū)波形存在明顯畸變，同相軸連續(xù)性較差；結合地質勘察資料與掌子面情況，可以推斷該段圍巖較破碎，節(jié)理裂隙較發(fā)育，穩(wěn)定性較差，施工時應加強支護。

圖8 2#掌子面?zhèn)鹘y(tǒng)濾波與雙樹復小波去噪結果對比

4 結論

受川藏鐵路極其復雜的地質條件和施工環(huán)境的影響，探地雷達地質超前預報技術采集的信號噪聲干擾較大，給信號處理與數(shù)據(jù)分析帶來較大困難。在對探測信號實施雙樹復小波變換的基礎上，采用通用閾值、無偏似然估計閾值、啟發(fā)式閾值、極大極小閾值4種閾值選擇方法，結合軟閾值、硬閾值、半軟閾值、非負消減4種閾值處理方法，對小波高頻系數(shù)進行閾值量化處理后進行雙樹復小波逆變換重構信號，達到了較好的去噪效果。對比分析歸一化均方差、信噪比、能量比3個評價因子后發(fā)現(xiàn)，采用無偏似然估計閾值選擇方法，結合非負消減閾值處理方法，可以獲得最佳的去噪效果，信號得到很大的改善，噪聲殘留較少，多次反射波被減弱，提高了超前地質預報中異常體的識別精度，能夠為川藏鐵路隧道施工提供更好的安全保障。