二維物理小波標(biāo)架去除探地雷達(dá)信號(hào)隨機(jī)噪聲

石顯新 楊秋芬 侯彥威

（1.中煤科工集團(tuán)西安研究院，陜西 西安 710077；2.西安文理學(xué)院物理系，陜西 西安 710065）

引 言

探地雷達(dá)（簡稱GPR）是利用高頻無線電波確定介質(zhì)內(nèi)部物質(zhì)分布規(guī)律的一種探測(cè)方法，近年來在工程與環(huán)境等淺層地球物理領(lǐng)域得到廣泛應(yīng)用［1-2］。

在探地雷達(dá)技術(shù)發(fā)展過程中，對(duì)接收信號(hào)中的噪聲處理是一個(gè)重要方面［3］。在探地雷達(dá)接收信號(hào)中存在多種干擾與噪聲，因探測(cè)目標(biāo)反射回波信號(hào)較弱，甚至可能掩沒于干擾與噪聲中，直接影響了目標(biāo)信號(hào)的檢測(cè)與判斷，這給數(shù)據(jù)的處理與解譯造成了很大障礙。為獲得高精度高分辨率的數(shù)據(jù)處理與解譯結(jié)果，需對(duì)采集的數(shù)據(jù)進(jìn)行預(yù)處理，以抑制或去除干擾與噪聲，提高數(shù)據(jù)的質(zhì)量。

Larsoneur和Morlet［4］首先提出用小波變換去噪的思想，從20世紀(jì)90年代開始有大量的將小波變換用于信號(hào)去噪的文章發(fā)表。朱光明等［5］把小波變換用于一維濾波。Chakraborty和David等［6］對(duì)比了短時(shí)Foureir變換、連續(xù)小波變換（CWT）及匹配追蹤算法（MPD）對(duì)地震信號(hào)作時(shí)間-頻率域分析的效果，指出MPD優(yōu)于前兩種，且給出了MPD去噪的方法。Miao和 Moon［7］及羅國安等［8］把小波變換用于消除地震記錄的面波干擾。章珂等［9］把二進(jìn)小波變換用于信號(hào)的分時(shí)分頻去噪處理。

在地震資料數(shù)據(jù)處理中二維物理小波標(biāo)架的成功應(yīng)用，是基于地震炮集數(shù)據(jù)中有效信號(hào)（體波）的空間形態(tài)是雙曲線形狀，二維物理小波函數(shù)可以很好地與之匹配；而主要的噪聲（面波）與體波在空間形態(tài)上差異很大。以二維物理小波為母小波，有效信號(hào)（體波）在小波域比較稀疏，而噪聲（面波）不稀疏，這樣在表示體波的有限區(qū)域重構(gòu)信號(hào)，就可以衰減面波干擾。將二維物理小波標(biāo)架用于去除探地雷達(dá)信號(hào)隨機(jī)噪聲，其原理類似，即探地雷達(dá)點(diǎn)目標(biāo)信號(hào)的二維觀測(cè)數(shù)據(jù)在空間上也呈雙曲線形態(tài)，而隨機(jī)噪聲在空間的分布沒有規(guī)律。若將含有噪聲的探地雷達(dá)信號(hào)變換到二維物理小波域，由于二維物理小波函數(shù)與點(diǎn)目標(biāo)信號(hào)能很好地匹配，小波域表示點(diǎn)目標(biāo)信號(hào)的小波系數(shù)就比較稀疏，而表示隨機(jī)噪聲的小波系數(shù)不稀疏。因此，在表示點(diǎn)目標(biāo)信號(hào)的小波系數(shù)區(qū)域重構(gòu)信號(hào)，就可以衰減隨機(jī)噪聲。

基于以上認(rèn)識(shí)，作者開展了二維物理小波標(biāo)架去除探地雷達(dá)信號(hào)隨機(jī)噪聲的研究，并用二維物理小波標(biāo)架閾值去噪對(duì)兩個(gè)典型信號(hào)進(jìn)行了去噪處理，證明了該方法的有效性。

1 二維物理小波標(biāo)架的數(shù)學(xué)原理

二維小波標(biāo)架原子具有二維小波的形式［10］。這里使用文獻(xiàn)［11］中構(gòu)造的二維物理小波，其數(shù)學(xué)表達(dá)式為

一維墨西哥帽小波表達(dá)式為

在x-t域的雙曲線表達(dá)式為

在x方向加上衰減因子e-ax2后物理小波表達(dá)式變?yōu)?/p>

加入衰減因子的目的在于使這種新構(gòu)造的物理小波滿足小波的允許條件。

式（1）、（4）中的a為正實(shí)數(shù)，反映x方向的衰減；H為雙曲線的頂點(diǎn)，可以根據(jù)要求設(shè)置。

物理小波標(biāo)架的原子具有以下形式為

適當(dāng)?shù)倪x取時(shí)間和偏移兩個(gè)方向的尺度因子，如

把式（6）中的（b）以式（1）代入后結(jié)果為1.證明物理小波適當(dāng)選取尺度位移因子后構(gòu)成的二維小波標(biāo)架，可運(yùn)用于分解重構(gòu)二維信號(hào)。

小波的尺度因子是局部化分析中改變小波頻譜中心和帶寬的參量。同一維小波標(biāo)架對(duì)信號(hào)的分解重構(gòu)一樣，尺度因子選擇應(yīng)使得各尺度小波頻譜對(duì)有效信號(hào)頻譜完全覆蓋且盡量覆蓋均勻，即各尺度小波的頻譜峰值覆蓋區(qū)域保證錯(cuò)開。這相當(dāng)于一濾波器組對(duì)信號(hào)各頻帶分別濾波。小波滿足正交條件時(shí)，各尺度的頻帶完全不重合，重構(gòu)效果好，而本文中介紹的物理小波不滿足正交條件，這就應(yīng)仔細(xì)考察選用的每一組尺度下的小波頻帶覆蓋范圍。在具體實(shí)現(xiàn)過程中，按二進(jìn)情況，在頻域（f-k），t方向尺度st倍增，頻譜在f方向的頻帶范圍壓縮為原來的一半；x方向尺度sx倍增，頻譜在k方向的頻帶范圍壓縮為原來的一半。從用Matlab繪制的各個(gè)尺度的物理小波f-k頻譜圖可以清楚的觀察其二維頻譜形狀、峰值、變化趨勢(shì)，從而選擇合適的尺度范圍，更好地對(duì)信號(hào)進(jìn)行局部化分析去除噪聲，重構(gòu)有效信號(hào)。

當(dāng)尺度增大到一定程度時(shí)，物理小波二維頻譜集中在f-k域中心，范圍很小，具體實(shí)現(xiàn)算法中發(fā)現(xiàn)頻譜峰值重復(fù)覆蓋，低頻信息多次疊加，表現(xiàn)為k低頻分量對(duì)重構(gòu)有效信號(hào)干擾大，在x向雙曲線信號(hào)頂點(diǎn)位置形成一條"直杠"。

通過觀察時(shí)間及空間尺度變化時(shí)二維物理小波在時(shí)間-尺度域及頻率-波數(shù)域的特征，發(fā)現(xiàn)二維物理小波函數(shù)在頻率-波數(shù)域呈扇形分布，這時(shí)二維物理小波具有最高的頻率及波數(shù)分辨率。小波函數(shù)在時(shí)間方向上明顯變寬，也就是說時(shí)間分辨率變低，頻帶變窄。小波函數(shù)在空間方向上明顯變寬，即空間分辨率變低，也就是波數(shù)域展布變窄。

2 二維物理小波標(biāo)架去噪的實(shí)現(xiàn)

2.1 物理小波標(biāo)架去噪的實(shí)現(xiàn)

實(shí)現(xiàn)物理小波標(biāo)架去噪的步驟如下:

1）人工合成點(diǎn)目標(biāo)探地雷達(dá)數(shù)據(jù)。為了逐步驗(yàn)證效果，分別制作了理想的無噪數(shù)據(jù)及含有色噪聲污染的數(shù)據(jù)。2）實(shí)現(xiàn)數(shù)據(jù)的物理小波標(biāo)架分解，得到每一組時(shí)間-偏移尺度下的分解系數(shù)。3）通過實(shí)驗(yàn)選取閾值，實(shí)現(xiàn)小波域閾值去噪。4）從處理后的小波域系數(shù)重構(gòu)有用信號(hào)，完成去噪。

2.2 信號(hào)分解和重構(gòu)公式及實(shí)現(xiàn)方法

在進(jìn)行信號(hào)分解時(shí)，需要固定兩個(gè)尺度因子，計(jì)算如下的二重積分

分別對(duì)信號(hào)f（x，t）和物理小波ψ（x、t）沿兩個(gè)方向采樣，x、t方向采樣分別為M、N點(diǎn)，構(gòu)成M×N階的矩陣，當(dāng)x、t單位都?xì)w一化后，就不必考慮兩個(gè)方向采樣的差異，視為采樣間隔均為1，這樣便可以按數(shù)字信號(hào)來實(shí)現(xiàn)。分解公式變?yōu)?/p>

按照小波變換信息的冗余性及實(shí)際應(yīng)用中的要求，選取尺度因子的一些離散值。按照常用的二進(jìn)小波變換的方法，令

分別固定兩個(gè)尺度因子，小波域系數(shù)不抽取時(shí)，式（8）相當(dāng)于計(jì)算一個(gè)二維相關(guān)函數(shù):

令y′=y（-m，-n），則二維相關(guān)函數(shù)計(jì)算化為二維卷積:

利用二維離散傅里葉變換（DFT）及逆變換（IDFT）的性質(zhì)，可以通過頻域計(jì)算卷積

式（12）中DFT和IDFT可由二維快速傅里葉正、逆變換FFT2和IFFT2實(shí)現(xiàn)。

重構(gòu)信號(hào)時(shí)，要滿足緊標(biāo)架重構(gòu)公式:

記WTsx，st（bx，bt）為一組尺度下的二維小波域系數(shù)矩陣，式（13）離散后變?yōu)?/p>

可分別對(duì)式（14）進(jìn)行二維卷積，重構(gòu)每一組尺度下的二維信號(hào)，再累加，最后就能得到重構(gòu)信號(hào)。

2.3 信號(hào)分解重構(gòu)舉例

圖1（a）為一探地雷達(dá)點(diǎn)目標(biāo)模擬散射信號(hào)，其中空間方向256個(gè)采樣點(diǎn)，采樣間隔1 cm，時(shí)間方向1024個(gè)采樣點(diǎn)，共40ns.可以看出，反射信號(hào)有較明顯的雙曲形態(tài)。每一道信號(hào)的反射子波為Ricker子波，其主頻為800MHz.圖1（b）、（c）是利用式（8）對(duì)圖1（a）作小波分解后由式（14）重構(gòu)的結(jié)果。比較兩圖，可以看出基本上重構(gòu)了信號(hào)。

圖1（e）與圖1（f）是不同距離時(shí)的單道原始信號(hào)與重構(gòu)信號(hào)的對(duì)比，細(xì)實(shí)線代表原始信號(hào)，點(diǎn)劃線代表重構(gòu)信號(hào)。信號(hào)的邊界重構(gòu)不完全，這與一維標(biāo)架類似，因?yàn)槎S物理小波不構(gòu)成緊標(biāo)架，信號(hào)重構(gòu)存在誤差。但在實(shí)際信號(hào)處理中這些誤差是可以接受的。

2.4 小波域閾值去噪的實(shí)現(xiàn)

小波域閾值去噪方法的原理是:分別在各個(gè)尺度下選取去噪的閾值 ，然后對(duì)小波域系數(shù)按以下規(guī)則處理:

物理小波標(biāo)架去噪的關(guān)鍵是當(dāng)數(shù)據(jù)投影到物理小波標(biāo)架上時(shí)，信息是冗余的，雙曲線形狀的有效信號(hào)的二維頻譜與某個(gè)尺度的物理小波標(biāo)架原子二維頻譜相似，故可以映射出絕對(duì)值較大的系數(shù)。因二維頻譜重合區(qū)域相關(guān)度小，使得有色噪聲映射的系數(shù)較小，去噪效果較好。具體實(shí)現(xiàn)中，可以分別從每組尺度下小波域系數(shù)中挑選典型的道，確定該組尺度應(yīng)選取的閾值，從而達(dá)到壓制噪聲的目的。

根據(jù)式（3）、（9）及小波標(biāo)架的算法，結(jié)合小波域閾值去噪的原理，就可以實(shí)現(xiàn)二維物理小波域閾值去噪。

3 模型算例

圖2為一探地雷達(dá)點(diǎn)目標(biāo)模擬散射信號(hào)，其中空間方向256個(gè)采樣點(diǎn)，采樣間隔1cm，時(shí)間方向1024個(gè)采樣點(diǎn)，共40ns.可以看出，反射信號(hào)有十分明顯的雙曲線形態(tài)。每一道信號(hào)的反射子波為Ricker子波，主頻800MHz.需要說明的是，圖2與圖1（a）都是點(diǎn)目標(biāo)散射模擬信號(hào)，但模擬時(shí)參數(shù)選擇不同，圖2的雙曲線形態(tài)更為明顯。

圖2 理想模擬信號(hào)

3.1 信噪比1∶1

在圖2的基礎(chǔ)上加上噪聲，其中信號(hào)與噪聲最大幅值比為1∶1，如圖3（a）所示，信號(hào)受到了較大影響，由于相關(guān)性，仍能看到信號(hào)的基本形態(tài)。為了更清楚地看到噪聲對(duì)信號(hào)的影響程度，從圖3（a）中抽取距離為128cm處的一維信號(hào)與原始無噪信號(hào)重疊顯示，如圖3（b）所示，其中點(diǎn)劃線為含噪聲信號(hào)，細(xì)實(shí)線為原始信號(hào)，可以看出，從一維信號(hào)中無法分辨有效信號(hào)的位置及形態(tài)。另外，因?yàn)樗釉肼暸c有效信號(hào)同頻帶，所以傳統(tǒng)的一維濾波不能去除噪聲。圖3（c）是對(duì)圖3（a）閾值去噪結(jié)果，可以看出:噪聲受到了很好的壓制，同時(shí)信號(hào)也受到了一定的影響，但點(diǎn)目標(biāo)的反射形態(tài)仍然比較清晰。圖3（d）為從圖3（c）中抽取距離為128cm處的一維信號(hào)與原始無噪信號(hào)重疊顯示，其中點(diǎn)劃線為去噪信號(hào)，細(xì)實(shí)線為原始信號(hào)，也可以看出，噪聲受到了較好的壓制。

3.2 信噪比1:3

與圖3（a）相比，圖4（a）中信號(hào)與噪聲最大幅值比為1∶3，信號(hào)已完全淹沒于噪聲中。從圖4（a）中抽取距離為128cm處的一維信號(hào)與原始無噪信號(hào)重疊顯示，如圖4（b）所示，其中點(diǎn)劃線為含噪聲信號(hào)，細(xì)實(shí)線為原始信號(hào)，噪聲明顯大于信號(hào)。圖4（c）是對(duì)圖4（a）閾值去噪結(jié)果，可以看出:噪聲受到了較大的壓制，點(diǎn)目標(biāo)的反射形態(tài)比較清晰。圖4（d）為從圖4（c）中抽取距離為128cm處的一維信號(hào)與原始無噪信號(hào)重疊顯示，其中點(diǎn)劃線為去噪信號(hào)，細(xì)實(shí)線為原始信號(hào)，可以看出，信號(hào)的基本形態(tài)得到了恢復(fù)。與圖4（c）相比，雖然效果要差一些，但這是在信噪比很低，一般認(rèn)為是廢棄資料的情況下，取得這樣的處理結(jié)果，也說明了該方法的有效性。

圖4 信噪比為1∶3時(shí)的去噪效果

4 結(jié) 論

綜上所述，利用二維物理小波函數(shù)模擬與探地雷達(dá)點(diǎn)目標(biāo)散射信號(hào)匹配的小波函數(shù)，對(duì)含噪聲的探地雷達(dá)模擬信號(hào)進(jìn)行小波域閾值去噪，在信噪比較低的情況下仍能取得較好的效果，證明了該方法的有效性。

［1］陳文超，師振盛，汪文秉，等.小波變換在去除探地雷達(dá)信號(hào)直達(dá)波的應(yīng)用［J］.電波科學(xué)學(xué)報(bào)，2000，15（3）:352-357.CHEN Wenchao，SHI Zhensheng，WANG Wenbing，et al.Suppressing the direct wave noise in GPR data using the 2-D continue directional wavelets［J］.Chinese Joural of Radio Science，2000，15（3）:352-357.（in Chinese）

［2］郭 晨，劉 策，張安學(xué).探地雷達(dá)超寬帶背腔蝶形天線設(shè)計(jì)與實(shí)現(xiàn)［J］.電波科學(xué)學(xué)報(bào)，2010，25（2）:221-226.GUO Chen，LIU Ce，ZHANG Anxue.Design and implement of an UWB bowtie antenna with backcovity for ground penetrating radar［J］.Chinese Journal of Radio Science，2010，25（2）:221-226.（in Chinese）

［3］周輝林，田 茂，熊俊志.探地雷達(dá)回波信號(hào)的特征提取和分類［J］.電波科學(xué)學(xué)報(bào)，2006，21（4）:586-589.ZHOU Huilin，TIAN Mao，XIONG Junzhi.Feature extraction and classification of ground penetrating radar echo signals［J］.Chinese Journal of Radio Science，2006，21（4）:586-589.（in Chinese）

［4］LARSONNEUR J L，MORLET J.Wavelets and seismic interpretation［C］／／ Wavelets，Time-Frequency Methods and Phase Space，1st Int Wavelets Conf Marseille，December 14-18，1989:126-141.

［5］朱光明，高靜懷，王玉貴.小波變換及其在一維濾波中的應(yīng)用［J］.石油物探，1993，32（1）:1-10.ZHU Guangming，GAO Jinghuai，WANG Yugui.Wavelet transform and its application to 1-D filtering［J］.Geophysical Prospecting for Petrole，1993，32（1）:1-10.（in Chinese）

［6］CHAKRABORTY A.Application of wavelet transform to seismic data［C］／／64th SEG Meeting American.Los Angeles，October 23-28，1994:725-728.

［7］MIAO X G.Application of the wavelet transform in seismic data processing［C］／／62th SEG meeting A-merican，1994:1461-1464.

［8］羅國安，杜世通.小波變換及信號(hào)重建在壓制面波中的應(yīng)用［J］.石油地球物理勘探，1996，31（3）:338-349.LUO Guoan，DU Shitong.Application of wavelet transform and signal reconstruction in surface wave elimination［J］.Oil Geophysical Prospecting，1996，31（3）:338-349.（in Chinese）

［9］章 珂，劉貴忠.二進(jìn)小波變換的地震信號(hào)分時(shí)分頻去噪處理［J］.地球物理學(xué)報(bào)，1996，39（2）:265-270.ZHANG Ke，LIU Guizhong，ZOU Dawen，et al.Seismic data time-frequency domain denosing based on dyadic wavelet transform［J］.Chinese Journal of Geophysics，1996，39（2）:265-270.（in Chinese）

［10］ZHANG Rongfeng，ULTYCH T J.Physical wavelet frame denoising［J］.Goephysics，2003，68（1）:225-231.

［11］楊秋芬，石顯新.探地雷達(dá)信號(hào)直達(dá)波的去除［J］.煤炭學(xué)報(bào)，2008，33（7）:770-774.YANG Qiufen，SHI Xianxin.Removing for the direct wave noise in ground penetrating radar［J］.Journal of China Coal Society，2008，33（7）:770-774.（in Chinese）