劉宗輝，吳一帆，劉保東，劉毛毛，藍日彥，孫懷鳳

1) 廣西新發(fā)展交通集團有限公司，南寧 530028 2) 廣西大學土木建筑工程學院，南寧 530004 3) 山東大學土建與水利學院，濟南 250061

4) 南寧城建管廊建設投資有限公司，南寧 530219

探地雷達（GPR）達近年來已成為隧道超前地質預報中最主要的短距離物探手段，其在不良地質探測方面具有分辨率高、結果直觀、掃描速度快等其它物探方法無法比擬的優(yōu)勢[1-2]. 探地雷達信號是一種典型的非平穩(wěn)、時變信號[3]，電磁波在復雜的隧道圍巖中傳播時存在強烈的吸收衰減、色散[4]，同時由于隧道探測環(huán)境中大量系統(tǒng)干擾，使得采集到的雷達反射波數據通常具有“弱信號，強干擾”的特征，給數據的處理和解釋帶來極大困難. 因此，干擾消除一直是探地雷達隧道地質預報應用中普遍關注而又沒有得到很好解決的難題.

隨著復信號分析技術的發(fā)展，探地雷達信號去噪已從中值濾波、頻域濾波、F-K濾波等傳統(tǒng)方法的基礎上發(fā)展為小波變換、曲波變換（curverlet變換）以及經驗模態(tài)分解（EMD）等方法，如柳鋼等[5]、Bao 等[6]、Gan 等[7]的研究成果. 這些方法在一定程度上提高了探地雷達的數據解譯精度，但大都仍存在各自難以克服的缺點. 如Ouadfeul等[8]、李才明等[9]采用的小波分析方法盡管充分利用了小波變換多尺度分析特性，但是該方法的有效性依賴于選定的小波基函數和設定的軟硬閥值，不能適用于復雜多變的探地雷達實測環(huán)境，且常用的二維小波對二維信號中直線或曲線等邊緣特征難以精確表達. 經驗模態(tài)分解[10-12]方法盡管充分利用了經驗模態(tài)分解分解的多分辨率和局部時頻分析特性，但已有的研究成果基本都是針對指定探地雷達信號，通過經驗人員分析，才能使經驗模態(tài)分解分解技術達到探地雷達信號降噪的目的，同時由于探地雷達信號是超寬帶信號，其中目標體一次回波信號、多次回波信號、雜波和各種噪聲信號會發(fā)生嚴重的頻率混疊，使得經驗模態(tài)分解分解不能在探地雷達探測中對接收信號進行有效增強處理. 曲波變換[13-14]是在小波變換的基礎上，增加了一個方位參數，解決了小波變換在處理二維信號時的不足，但該方法在干擾信息與有效信息方向性一致的情況下，分離效果并不理想. 曲波變換閥值函數的選取直接關系到曲波降噪效果，找到同時適應尺度和角度的閥值函數是該方法成功應用亟需解決的問題. 此外，目前探地雷達信號去噪研究重心多集中在對算法的更新而忽略了對干擾數據本身特征的研究，雖然有學者研究總結了隧道超前地質預報中的干擾類型[15]，但沒有從信號處理角度給出具有針對性的干擾壓制方法.

剪切變換（shearlet變換，ST）是一種較新的多尺度多方向時頻分析技術，它相比于小波變換（wavelet變換，WT）和其他多尺度幾何分析方法有更好的方向敏感性，信號保真度高，該技術在地震波消噪方面已被證實有較好的應用效果[16-17]. 由于地震波和雷達電磁波的諸多相似性，本文引進剪切變換，利用數學理論框架將其改進，并將其與小波變換相結合，提出剪切變換與小波變換聯(lián)合去除干擾方法. 實際案例處理效果表明該方法在保證去干擾效果同時能較好地保留有效信號.

1 基于剪切變換的自適應閥值去噪

1.1 基本方法原理

1.1.1 小波變換基本原理

小波變換是在傅立葉變換基礎上發(fā)展起來的，其最大優(yōu)勢是可以由粗到細逐步觀察信號，能很好地表征信號局部特征，對于信號頻率具有很高的敏感性.

連續(xù)小波變換的表達式為：

連續(xù)小波變換逆變換表達式為：

小波變換處理數據時基函數的選擇尤為重要，小波波形越接近待處理的瞬態(tài)信號波形，處理效果越理想. 綜合時頻域的分辨率來看，DB族小波是比較適合分析探地雷達信號的一種小波基函數[15].

1.1.2 剪切變換基本原理

剪切變換又稱剪切波變換，最初是由Guo和Easley等[18-19]通過對剪切基函數進行縮放、剪切和平移合成具有膨脹性的仿射系統(tǒng). 當維數為2時，該合成膨脹仿射系統(tǒng)定義為：

按照以下方式采樣，可將剪切變換離散化：

離散剪切系統(tǒng)為：

1.2 自適應閥值去噪

含噪雷達信號經過剪切變換后，數據信息被映射到不同尺度不同方向的剪切系數上，通常有效信號會根據自身特點集中在一些特定方向上，而隨機噪聲信號不具有方向性. 有效信號所映射的某些特定方向上的剪切系數值往往較大，而隨機噪聲將在各個方向上分布，其對應的剪切系數值往往較小.

從閥值去噪角度來講，如果某一分解方向中剪切系數值較大，則該方向為有效信號系數主要集中的方向，對于該方向上的系數處理應該采取較小的閥值以便更好地保護有效信號；當某一分解方向中剪切系數值較小時，則該方向的系數主要是噪聲系數，應該采取較大的閥值以便能壓制更多的噪聲. 因此，在經典剪切算法中[16]，僅設置一個隨尺度變化的閥值難以滿足方向上能量變化的需求，容易產生“過扼殺”、“除不凈”的現象.

本文在經典尺度閥值函數的基礎上，根據有效信號和干擾信號在剪切域中不同尺度、不同方向上能量的差異，設置一個隨尺度和方向變化的閥值函數，從而使得在剪切域進行噪聲壓制時每一個子帶都可以根據其能量特征自適應選擇最優(yōu)閥值，其表達式如下：

基于剪切變換的自適應閥值去噪具體過程包括以下步驟：

（1）信號常規(guī)處理，包括：直達波去除、直流去漂移、信號增益等；

（2）格式轉換，將雷達數據格式轉換為矩陣形式以保證算法的實現；

（3）偽極坐標變換，將數據從笛卡爾坐標系轉換到偽極坐標系，并且在偽極坐標系上進行多尺度劃分，并生成剪切基函數對數據進行窗口子帶化；

（4）計算剪切系數，得到系數矩陣之后通過上述自適應閥值函數對剪切系數矩陣進行干擾壓制處理；

（5）根據去噪處理后的剪切系數對探地雷達二維數據進行重構.

2 隨機干擾消除

2.1 正演模擬及干擾設置

利用時域有限差分法分別模擬圓形空洞與方形空洞兩種異常體，得到純凈的雷達數據. 圖1和圖2分別為異常體幾何模型和波形堆積圖，從圖中可以看出異常體反射界面清晰，同相軸特征明顯.

圖1 正演幾何模型. （a）圓形空洞；（b）方形空洞Fig.1 Geometric model of forward simulation: (a) circular hole; (b) square hole

圖2 正演模擬結果. （a）圓形空洞；（b）方形空洞Fig.2 Forward simulation results: (a) circular hole; (b) square hole

在正演模擬獲得的純凈數據中加入高斯白噪聲，加噪后圖像信噪比為-2.5 dB左右，波形堆積圖如圖3所示. 從圖中可以看出異常體反射界面的有效信號已基本被淹沒在噪聲之中，圓形空洞下反射面同相軸幾乎不可見，方形空洞數據只能觀察到水平部分能量特別強的部位.

圖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?5. 處理后的波形堆積如圖4和圖5所示.

從圖4和圖5中可以很直觀的看出小波變換處理后同相軸信息并沒有被完全還原出來，圓形空洞數據的下反射面同相軸模糊，方形空洞數據上反射面邊緣不清晰，數據整體依然含有較多噪聲，去噪效果并不理想. 而使用本文提出的剪切變換去噪后的數據圖像，目標信號被極大限度的還原，圓形空洞數據中被噪聲掩蓋的下反射面同相軸在經過處理后清晰可見，幾乎與原始信號數據沒有差異，方形空洞的起止位置處的信號也被有效的還原出來. 從方形空洞處理結果同時可以看出剪切變換處理后數據與原數據相比信號強度稍弱，說明存在少量的有效信號被過度剔除.

圖4 小波變換處理結果. （a）圓形空洞；（b）方形空洞Fig.4 Results after Wavelet transform processing: (a) circular hole; (b) square hole

圖5 剪切變換處理結果. （a）圓形空洞；（b）方形空洞Fig.5 Results after shearlet transform processing: (a) circular hole; (b) square hole

利用信噪比（SNR）、峰值信噪比（PSNR）和均方誤差（MSE）三個指標對去噪前后波形堆積圖數據質量進行定量分析. 由表1可知，剪切變換處理后數據的信噪比與峰值信噪比均大于小波變換處理結果，而均方誤差由個位數縮小到0.1以內，進一步說明剪切變換自適應閥值法在處理隨機噪聲上的優(yōu)勢，數據噪聲殘留少.

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

圖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

綜上所述，本文提出基于自適應閥值的剪切變換方法可以很好地適應探地雷達的數據結構，對隨機噪聲有很好的壓制能力，可以提高含噪數據信噪比.

3 頻率異常干擾消除

探地雷達脈沖信號是一種寬帶電磁波，具有非平穩(wěn)性、非線性衰減等特點. 在隧道掌子面探測時，由于外界信號干擾、地下不同介質的吸收、反射等因素，導致儀器所接收到的反射回波往往是多種頻率成分的信號疊加. 小波變換具有良好的時頻局部化性質，其可變的時窗結構對一維信號具有較強的頻率分辨率. 本文采用的小波去除頻率異常信號的基本原理是：首先對信號進行若干尺度層小波分解，獲得各尺度層不同頻率的小波系數；再結合時頻分析結果，取合適的閥值進行帶通濾波，去除頻率異常部分的信號小波系數；最后再進行剩余小波系數分量重構，從而達到頻率異常信號消除.

以岑溪大隧道左洞掌子面DK7+511處探地雷達實測數據來說明上述兩種方法對能量接近的不同頻率成分信號混疊干擾的去除效果. 現場探測采用意大利IDS公司K2雷達，天線頻率為100 MHz，時窗600 ns，采樣點數512.

圖7（a）為常規(guī)方法處理后的探地雷達圖像，從圖中可以看到明顯貫穿整個剖面的強反射波組，對探測范圍內地質異常解釋存在較大干擾.圖8（a）為該測線8.3 m處單道波的廣義S變換后的時頻分布[20]，從圖中也可以看出沿時間深度存在明顯低頻成分. 結合現場探測環(huán)境以及實際開挖情況可判定此處低頻成分為干擾信號. 圖7（b）和7（c）分別為剪切變換和小波變換處理后二維雷達剖面，處理過程中具體參數與模型試驗一致，對比兩幅圖可以看出剪切變換能將雜波信號很好去除，但對低頻干擾信號處理效果并不明顯；而小波變換處理后波長異常區(qū)域得到明顯改善，但數據淺部仍存在較多的雜波信號.

圖7 頻率異常信號干擾消除前后灰度圖. （a）常規(guī)方法；（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.3 m 處單道波的時頻分布. 對比兩幅圖可以更直觀看出小波變換能很好地將低頻成分分離，對頻率異常信號去除優(yōu)勢明顯.

4 小波變換與剪切變換聯(lián)合法去干擾

圖8 頻率異常干擾消除前后單道波廣義S變換時頻分布. （a）常規(guī)方法；（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

根據干擾信號特點，可以將隧道內常見的干擾分為隨機干擾和頻率異常干擾兩種類型，其中隨機干擾具有頻率隨機分布、波形雜亂的特點，而頻率異常干擾往往具有某些特定頻率特征、在波形堆積圖上往往具有某些特定規(guī)律. 通過對正演模擬和現場實際數據的處理可以發(fā)現，剪切變換與小波變換對干擾的壓制都有著各自的優(yōu)勢：剪切變換對信號能量敏感，對于隨機噪聲、機械噪聲以及電磁信號等頻率隨機、能量異常的干擾壓制效果較好，而小波變換對頻率異常、能量相近的干擾壓制效果較好.

結合兩種方法各自的優(yōu)勢，進一步提出小波變換與剪切變換相結合去除干擾方法，即先用小波變換對異常頻率信號進行分離，再使用剪切變換對隨機干擾進行壓制. 具體實現流程圖如圖9所示.

圖9 聯(lián)合法去除干擾流程圖Fig.9 Flow chart of the combined methods for interference removal

5 實際工程案例

現場試驗場地位于廣西壯族自治區(qū)融水縣至河池市在建高速公路羅城段，探測對象為路基邊緣三處形態(tài)不規(guī)則裸露的溶洞，其中1#為軟塑黏土夾碎石充填型溶洞、2#為干土充填型溶洞、3#為空腔型溶洞. 路基為微風化灰?guī)r、巖質堅硬、結構面較發(fā)育，周邊未見地表水. 現場情況、測線與溶洞平面示意圖分別如圖10和圖11所示，試驗坑長11 m，深1.5 m，實際測線長度10 m.

圖10 現場情況Fig.10 Field conditions

圖11 測線與溶洞平面示意圖Fig.11 Layout diagram of karst caves and survey line

采用意大利IDS公司K2雷達探測，天線頻率為100 MHz，時窗400 ns，采樣點數1024. 通過人為設置機械電磁噪聲和金屬干擾來模擬隧道干擾環(huán)境.

圖12（a）為常規(guī)方法處理后的波形堆積圖，從圖中可以明顯看出240 ns以下深度中存在強烈同相軸異常干擾，3#溶洞反射波信號被淹沒于干擾信號中，根據信號特征可判定該區(qū)域同時含有隨機干擾和頻率異常成分干擾. 此外，從圖中可以看出240 ns以上的淺部數據也存在波形雜亂無章的隨機噪聲干擾，左側1#和中間2#溶洞位置難以區(qū)分. 圖13（a）為測線3 m處單道波時頻分布，從圖中可清晰看出干擾信號頻率成分. 因此，為凸顯異常體空間位置以及后續(xù)進一步開展屬性分析，有必要采用本文所提出的聯(lián)合去噪方法提取異常體反射信號.

圖12 去噪效果對比. （a）常規(guī)方法；（b）小波變換；（c）聯(lián)合算法Fig.12 Comparison of different denoizing methods: (a) conventional method; (b) wavelet transform; (c) joint algorithm

圖13 廣義 S 變換時頻分布結果對比. （a）常規(guī)方法；（b）小波變換；（c）聯(lián)合算法Fig.13 Comparison of GST results obtained using different denoizing methods: (a) conventional method; (b) wavelet transform; (c) joint algorithm

圖12（b）和 12（c）分別為使用小波變換以及聯(lián)合法處理后的波形堆積圖，對處理后的雷達數據選取測線3 m處單道波進行時頻分析，時頻分布見圖 13（b）和 13（c）. 從圖 12（b）中可以明顯看出低頻成分的強反射同相軸已得到較好地去除，但隨機干擾并未得到有效去除. 而圖12（c）處理結果顯示整個數據的圖像質量有了很大的改善，同相軸變得清晰連續(xù)，深層以及淺層的隨機噪聲和低頻干擾信號都能很好的去除. 根據處理后的波形堆積圖數據可以很好地將三處地質異常進行圈定，并能與實際情況對應. 圖 13（b）和 13（c）時頻分布圖中同樣可進一步清晰看出干擾信號頻率成分去除效果. 通過該案例可以進一步說明小波變換與剪切變換聯(lián)合法去干擾的有效性與必要性.

6 結語

（1）利用剪切變換將探地雷達數據轉換到剪切域，可以更加細致的對信號進行多尺度多方向劃分，通常有效信號會根據自身特點集中在一些特定的方向上，在此基礎上提出的基于自適應閥值去噪方法可對隨機干擾有很好去除效果.

（2）小波變換具有良好的時頻局部化性質，其可變的時頻窗結構對于一維信號有較強的頻率分辨率，通過多尺度分解后可將頻率異常信號分離，從而保留有效信號. 但小波變換對于隨機分布在整個頻率域且能量異常的信號去除效果不如剪切變換.

（3）隧道空間環(huán)境復雜，使用探地雷達進行隧道超前地質預報時往往會遇到各種成分的干擾混疊，采用單一干擾去除方法難以獲得滿意效果. 提出的剪切變換與小波變換聯(lián)合方法可以同時對能量異常的隨機干擾以及頻率異常的干擾信號進行壓制，并且可以保證處理后的數據有著較高信噪比.

（4）根據具體的干擾類型以及干擾數據特征選擇合適的閥值系數，直接關系到干擾去除和有效信號保留效果，自適應閥值函數仍有許多不確定性. 實際工程應用中需要不斷總結隧道中常見各類干擾的剪切系數能量特征以及小波系數頻率特征，并形成干擾信號樣本庫.