雙樹復(fù)小波變換在川藏鐵路拉林段某隧道超前地質(zhì)預(yù)報(bào)中的應(yīng)用

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

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

0 引言

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

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

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

1 掌子面測線布置

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

圖1 探地雷達(dá)測線布置示意

2 基于雙樹復(fù)小波的信號處理方法

2.1 雙樹復(fù)小波變換

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

探地雷達(dá)接收的雷達(dá)信號x(t)可以分解為：

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

(1)

雙樹復(fù)小波變換對信號x(t)進(jìn)行兩路離散小波變換，即采用2個(gè)獨(dú)立的小波來計(jì)算信號的實(shí)部與虛部。如圖2所示，樹a為實(shí)部，樹b為虛部，由實(shí)部和虛部2個(gè)樹狀的共軛正交濾波器組構(gòu)成一個(gè)復(fù)小波，h0(n)和h1(n)分別對應(yīng)實(shí)部的低通濾波器和高通濾波器，g0(n)和g1(n)分別對應(yīng)虛部的低通濾波器和高通濾波器，且滿足：

圖2 雙樹復(fù)小波結(jié)構(gòu)示意

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

(2)

式中：G0(ejw)和H0(ejw)分別是g0(n)和h0(n)的傅里葉變換[12]。雙樹復(fù)小波的逆變換如圖3所示。

圖3 雙樹復(fù)小波逆變換

2.2 小波閾值去噪

2.2.1 小波去噪信號處理過程

小波閾值去噪過程分為3個(gè)階段：分解、閾值量化和重構(gòu)，在雙樹復(fù)小波變換的基礎(chǔ)上，將含噪信號進(jìn)行分解，保留最底層的低頻系數(shù)，僅對各個(gè)分解尺度下的高頻系數(shù)選擇一個(gè)閾值進(jìn)行量化處理，將處理過的小波系數(shù)通過雙樹復(fù)小波逆變換進(jìn)行重構(gòu)，進(jìn)而獲得去除噪聲后的信號(圖4)。這個(gè)過程中最關(guān)鍵的是如何選擇閾值，并利用該閾值對小波系數(shù)進(jìn)行量化處理，它關(guān)系到信號消噪的質(zhì)量。

圖4 雙樹復(fù)小波去噪流程示意

2.2.2 閾值選擇及閾值處理

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

1)通用閾值

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

(3)

式中：σ為噪聲的標(biāo)準(zhǔn)差，N為信號長度。

2)無偏似然估計(jì)閾值

無偏似然估計(jì)閾值是一種基于 Stein 的無偏似然估計(jì)原理的自適應(yīng)閾值選擇，是一種軟件閾值估計(jì)器。

3)啟發(fā)式閾值

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

4)極大極小閾值

極值閾值也是一種固定閾值選擇形式，它所產(chǎn)生的是一個(gè)最小均方差的極值，不是無誤差。主要應(yīng)用極大極小原理選擇閾值[3]。

(4)

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

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

1)硬閾值處理

(5)

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

2)軟閾值處理

(6)

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

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

3)半軟閾值處理

(7)

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

4)非負(fù)消減處理

(8)

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

2.3 去噪效果評價(jià)

為評價(jià)雙樹復(fù)小波變換去噪重構(gòu)后的性能，采用以下3個(gè)參數(shù)進(jìn)行評價(jià)。

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

(9)

2)信噪比SNR：

(10)

3)能量比P：

(11)

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

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

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

3.1 1#掌子面應(yīng)用情況

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

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

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

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

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

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

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

3.2 2#掌子面應(yīng)用情況

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

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

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

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

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

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

4 結(jié)論

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