曲波域統(tǒng)計(jì)量自適應(yīng)閾值探地雷達(dá)數(shù)據(jù)去噪技術(shù)*

李靜和 何展翔 楊俊? 孟淑君 李文杰 廖小倩

1) (桂林理工大學(xué)地球科學(xué)學(xué)院,桂林 541004)

2) (南方科技大學(xué)地球與空間科學(xué)系,深圳 518055)

非線性、非平穩(wěn)探地雷達(dá)數(shù)據(jù)常摻雜各種復(fù)雜噪聲源,其對(duì)精確提取弱反射波信號(hào)、識(shí)別繞射波雙曲線同相軸特征具有嚴(yán)重影響,忽略噪聲影響給探地雷達(dá)探測(cè)數(shù)據(jù)全波形偏移成像及后續(xù)解譯造成較大誤差.采用傳統(tǒng)閾值函數(shù)的曲波變換去噪需要根據(jù)數(shù)據(jù)噪聲水平人為確定合理閾值控制系數(shù).對(duì)此,本文開(kāi)展自適應(yīng)閾值函數(shù)的曲波變換去噪算法研究.引入塊狀復(fù)數(shù)域閾值函數(shù)算法,分析傳統(tǒng)閾值函數(shù)曲波變換去噪的效果隨閾值控制系數(shù)變化的規(guī)律;利用高階統(tǒng)計(jì)量理論,對(duì)曲波變換系數(shù)在尺度、方向上進(jìn)行相關(guān)性疊加,通過(guò)相關(guān)性統(tǒng)計(jì)量自適應(yīng)確定有效信號(hào)在曲波變換系數(shù)分布尺度、旋轉(zhuǎn)方向,由此確定清除噪聲成分閾值范圍,構(gòu)建統(tǒng)計(jì)量自適應(yīng)閾值函數(shù)曲波變換去噪算法.針對(duì)包含隨機(jī)噪聲、相關(guān)噪聲合成探地雷達(dá)數(shù)據(jù)及實(shí)測(cè)探地雷達(dá)數(shù)據(jù),采用傳統(tǒng)閾值函數(shù)曲波變換去噪與本文提出去噪算法處理結(jié)果對(duì)比分析,檢驗(yàn)了本文算法的有效性及可行性.研究成果對(duì)復(fù)雜探地雷達(dá)數(shù)據(jù)精確推斷解譯具有指導(dǎo)意義.

1 前 言

探地雷達(dá)勘探采用天線發(fā)射不同頻率電磁波,利用地下介質(zhì)的電磁性質(zhì)差異,根據(jù)接收回波的幅值、波形等運(yùn)動(dòng)學(xué)及動(dòng)力學(xué)特征推斷目標(biāo)介質(zhì)的空間、物性分布,廣泛應(yīng)用于地質(zhì)災(zāi)害監(jiān)測(cè)、工程與環(huán)境地質(zhì)勘察、水文地質(zhì)勘查及軍事偵查等領(lǐng)域.然而,隨著探地雷達(dá)勘探環(huán)境越來(lái)越復(fù)雜、目標(biāo)探測(cè)要求越來(lái)越精細(xì)[1,2],如何在強(qiáng)能量的各種復(fù)雜干擾噪聲湮沒(méi)環(huán)境中提取微弱的目標(biāo)信號(hào)是一項(xiàng)艱巨的任務(wù)[3,4].

常規(guī)中值濾波法[5]、S變換去噪[6-8]及F-K(頻率-波速)濾波方法[9,10]等探地雷達(dá)勘探數(shù)據(jù)去噪技術(shù)多基于簡(jiǎn)單最優(yōu)化或正交變換算法,構(gòu)建固定濾波窗口,去噪過(guò)程導(dǎo)致時(shí)頻域重疊微弱有效信號(hào)失真,具有特定噪聲處理的局限性.由于各種噪聲源摻雜,探地雷達(dá)數(shù)據(jù)常為非線性、非平穩(wěn)信號(hào)序列[11].基于傅里葉變換理論研發(fā)的信號(hào)處理技術(shù)可有效解決平穩(wěn)信號(hào)分析處理問(wèn)題,但非自適應(yīng)時(shí)頻窗口局限性使其無(wú)法解決探地雷達(dá)復(fù)雜數(shù)據(jù)的去噪問(wèn)題.小波變換用于探地雷達(dá)數(shù)據(jù)去噪的假設(shè)前提是有效信號(hào)與噪聲信號(hào)頻譜分離,通過(guò)采用具有伸縮和平移特性基函數(shù)在特定頻段的稀疏表示,設(shè)置噪聲頻段閾值以達(dá)到保真去噪的目的[12].如利用高頻閾值函數(shù)小波變換去噪理論開(kāi)展探地雷達(dá)弱信號(hào)提取,為推斷異常位置提供依據(jù)[13].改進(jìn)閾值的提升小波變換用于混凝土脫空探地雷達(dá)數(shù)據(jù)去噪,提高了解釋推斷精度[14].可協(xié)調(diào)方向閾值函數(shù)的小波變換用于考古學(xué)和巖土勘察探地雷達(dá)數(shù)據(jù)去噪及雜波壓制,改善了數(shù)據(jù)信噪比等[15,16].因而,高效的閾值函數(shù)可有效提升去噪效果,實(shí)現(xiàn)保真去噪的目的.

為突破小波變換采用矩形時(shí)頻窗口去噪的局限性,曲波變換采用小波域伸縮及平移特性,并且引入一個(gè)方向參量從而具有更好的方向識(shí)別能力[17,18].曲波變換自提出便迅速?gòu)V泛應(yīng)用于地球物理數(shù)據(jù)處理領(lǐng)域,特別是在地震數(shù)據(jù)噪聲壓制、多次波分離等領(lǐng)域,研究成果不斷涌現(xiàn)[19-22].其中,Neelamani等[23]將曲波變換應(yīng)用于三維地震數(shù)據(jù),通過(guò)設(shè)定閾值函數(shù)有效地清除隨機(jī)噪聲、相關(guān)噪聲,提高了地震資料的信噪比;張金良等[24]研究快速曲波變換域的SAR (synthetic aperture radar)圖像去噪算法,使用均值閾值濾波器使得圖像的可視化和解譯變得精確;Terrasse等[25]通過(guò)人工挑選探地雷達(dá)數(shù)據(jù)曲波變換域系數(shù),明確需要去除雜亂回波及噪聲源的尺度、方位信息,利用硬閾值函數(shù)實(shí)現(xiàn)去噪處理.可見(jiàn),曲波變換應(yīng)用于探地雷達(dá)數(shù)據(jù)噪聲壓制的前提條件是獲取較為準(zhǔn)確的噪聲閾值參數(shù)及其所屬曲波系數(shù)在尺度、方位上分布范圍.

實(shí)測(cè)復(fù)雜探地雷達(dá)數(shù)據(jù)常包含隨機(jī)噪聲、相關(guān)噪聲及其他未知噪聲類(lèi)型,所需閾值參數(shù)及所屬曲波系數(shù)在尺度、方位上分布范圍亦有所不同.如朱自強(qiáng)等[26]在曲波變換域中將選擇的角度窗函數(shù)方向因子設(shè)置為零,同時(shí)估計(jì)噪聲方差確定閾值函數(shù),實(shí)現(xiàn)去除表層直達(dá)波、噪聲源,在幅值強(qiáng)度較高的直達(dá)波背景下提取了鋼筋層和裂隙水層的弱化反射信號(hào).Bao等[27]認(rèn)為背景噪聲主要能量集中在曲波變換域方向90°區(qū)域附近,而隨機(jī)噪聲會(huì)相對(duì)均勻地分布在整個(gè)曲波域,采用二維濾波器估計(jì)噪聲分布.Tzanis[28]通過(guò)數(shù)值模擬人為地建立了不同裂隙構(gòu)造對(duì)應(yīng)探地雷達(dá)數(shù)據(jù)曲波變換域尺度、方位分布范圍,繼而實(shí)現(xiàn)特定方向發(fā)育裂隙結(jié)構(gòu)探地雷達(dá)反射波信號(hào)提取.因而,如何有效地確定目標(biāo)體引起的探地雷達(dá)數(shù)據(jù)曲波變換域尺度、方位上分布范圍是實(shí)現(xiàn)高效去噪、精確偏移成像的關(guān)鍵[29].當(dāng)前,曲波變換配置研發(fā)了諸多計(jì)算閾值的方法,如分段線性濾波方法、L2標(biāo)準(zhǔn)差法、曲波正變換法、對(duì)角實(shí)數(shù)及復(fù)數(shù)閾值法等,實(shí)測(cè)數(shù)據(jù)去噪需根據(jù)噪聲類(lèi)型選用,總體上無(wú)法滿足自適應(yīng)保真去噪的目的.因而,曲波域自適應(yīng)閾值的構(gòu)建是學(xué)者關(guān)注的研究熱點(diǎn).基于此,結(jié)合曲波變換與高階相關(guān)統(tǒng)計(jì)分析的優(yōu)勢(shì),本文提出在曲波變換多尺度、多方向分析理論基礎(chǔ)上采用高階相關(guān)統(tǒng)計(jì)量構(gòu)建自適應(yīng)閾值算法,探索探地雷達(dá)數(shù)據(jù)曲波域自適應(yīng)保真去噪方法技術(shù).

2 方法原理

2.1 曲波變換

擬采用第二代截?cái)嚯x散曲波變換算法[17,18]開(kāi)展去噪算法研究.曲波變換變量包含頻率、尺度及方位(角度),其變換表達(dá)式為

定義頻率窗口U、尺度窗口W和角度窗口V為

式中,j為尺度,m為旋轉(zhuǎn)方向,r和t為空間及時(shí)間域參數(shù).其中,角度窗口V為環(huán)形域,并滿足|x|∈[2j,2j+1]及極坐標(biāo)定義θj,m=2πm·2-［j/2].對(duì)尺度j、旋轉(zhuǎn)角度θj,m及空間位置

對(duì)f(x)∈L2(R2),曲波系數(shù)定義為

曲波域噪聲壓制主要思路為：1)對(duì)探地雷達(dá)數(shù)據(jù)進(jìn)行二維快速傅里葉正變換獲取系數(shù)分布;2)按(2)—(4)式配置頻率、尺度和角度窗口;3)根據(jù)給定閾值函數(shù)截?cái)嗳ピ?或?qū)⒃肼曀鶎傧禂?shù)范圍設(shè)置為零實(shí)現(xiàn)噪聲壓制;4)對(duì)3)部分處理剩余曲波系數(shù)利用二維快速傅里葉反變換算法獲取處理后的數(shù)據(jù)結(jié)果.曲波域噪聲壓制技術(shù)的關(guān)鍵是獲取較為準(zhǔn)確的閾值函數(shù)及其所屬系數(shù)范圍.對(duì)于合成探地雷達(dá)數(shù)據(jù),可根據(jù)加入的噪聲類(lèi)型、強(qiáng)度給出準(zhǔn)確的閾值函數(shù)及確定其所屬系數(shù)范圍;對(duì)于實(shí)測(cè)探地雷達(dá)數(shù)據(jù),常規(guī)操作是根據(jù)實(shí)際情況,人為選取合適的估計(jì)閾值函數(shù)及所屬系數(shù)范圍.

利用第二代曲波變換算法[17,18]開(kāi)展含噪數(shù)據(jù)處理對(duì)比研究,其中曲波變換涉及的閾值窗口在所有尺度設(shè)置為,δ為閾值分布范圍控制系數(shù),可人為設(shè)定;Ec為曲波變換系數(shù)L2范數(shù).δ值越大,過(guò)濾窗口越大,殘余信息越少;反之,δ值越小,噪聲及有效信息量越多.關(guān)于數(shù)據(jù)過(guò)濾判斷范圍的估計(jì),有多種方法供參考,如計(jì)算數(shù)據(jù)L2范數(shù)、塊狀復(fù)數(shù)閾值估計(jì)等[30].對(duì)于復(fù)雜數(shù)據(jù)處理,塊狀復(fù)數(shù)閾值函數(shù)更具有優(yōu)勢(shì),并將其用于曲波變換去噪過(guò)程,如(6)式所示：

式中,Ψ為復(fù)數(shù)域閾值函數(shù),用于估計(jì)判斷濾除范圍,具體推導(dǎo)請(qǐng)參考文獻(xiàn)[30];Γ1為相鄰道曲波系數(shù)權(quán)重值,Γ2為曲波系數(shù)開(kāi)平方歸一化和權(quán)重值.

2.2 曲波域統(tǒng)計(jì)量自適應(yīng)閾值去噪

對(duì)探地雷達(dá)數(shù)據(jù)做二維的曲波正變換獲取不同尺度、不同選擇角度的系數(shù)分布;對(duì)每一道探地雷達(dá)數(shù)據(jù)在不同尺度、不同角度下計(jì)算高階相關(guān)統(tǒng)計(jì)量,根據(jù)統(tǒng)計(jì)量分布明確閾值分布權(quán)重,實(shí)現(xiàn)噪聲濾除;最后采用二維的曲波反變換重建去噪后的探地雷達(dá)數(shù)據(jù).其中,未偏移的三階相關(guān)函數(shù)表示為

歸一化加權(quán)系數(shù)為

式中,i為探地雷達(dá)數(shù)據(jù)接收點(diǎn);t為時(shí)間采樣點(diǎn);q為平移因子,取值為1;為第i道探地雷達(dá)數(shù)據(jù)曲波域第j個(gè)尺度、第θ個(gè)方向正變換系數(shù).

該算法的優(yōu)勢(shì)在于不需要估計(jì)曲波變換閾值函數(shù)及其所屬系數(shù)范圍,而是通過(guò)統(tǒng)計(jì)理論建立目標(biāo)體引起的探地雷達(dá)有效信息所屬的尺度及旋轉(zhuǎn)角度范圍.其基本步驟如下.

步驟1,對(duì)原始探地雷達(dá)數(shù)據(jù)做曲波正變換獲取所有尺度j和選擇角度θ系數(shù)分布,并提取大尺度和小尺度的曲波系數(shù)和.

步驟2,采用曲波變換數(shù)值分析確定參與高階相關(guān)統(tǒng)計(jì)量計(jì)算的尺度范圍.

步驟3,Forj= 1,···,Jdo,

Forθ=1,···,?do.

Fori=1,···,M-1do.

end.

i=i+1 ;θ=θ+1;j=j+1 ;

end.

3 模擬試驗(yàn)

3.1 傳統(tǒng)曲波變換去噪

建立理論模型并采用二維有限差分程序計(jì)算合成模擬數(shù)據(jù).矩形目標(biāo)體規(guī)模為2 m × 1 m,X軸分布—1—1 m,Z軸分布為2—3 m;電阻率為10 Ω·m、介電常數(shù)為30 F/m;背景介質(zhì)電阻率為1000 Ω·m、介電常數(shù)為3 F/m;探地雷達(dá)觀測(cè)系統(tǒng)頻率為100 MHz,點(diǎn)距為0.1 m,X軸方向設(shè)置測(cè)點(diǎn)數(shù)為101,時(shí)間域采集點(diǎn)數(shù)為241.文中峰值信噪比(PSNR)定義為,其中MSE為數(shù)據(jù)均方差,表示原始信號(hào)和噪聲(或去噪后)信號(hào)的近似程度,max(s)為原始信號(hào)的峰值.

圖1(a)所示為無(wú)噪聲模擬合成探地雷達(dá)原始時(shí)間域數(shù)據(jù),可見(jiàn)在空間域設(shè)定矩形目標(biāo)模型邊界引起的探地雷達(dá)反射波信號(hào)分布于不同時(shí)間范圍.圖1(b)和圖1(c)分別為加入隨機(jī)噪聲(PSNR= 17 dB)、相關(guān)噪聲后的探地雷達(dá)數(shù)據(jù)(PSNR= 15.51 dB).椒鹽式無(wú)規(guī)律的隨機(jī)噪聲分布于整個(gè)探地雷達(dá)剖面,目標(biāo)弱反射信號(hào)完全被噪聲淹沒(méi);相關(guān)噪聲根據(jù)有效反射信號(hào)范圍比例分布,可見(jiàn)部分假反射信號(hào).采用(6)式閾值函數(shù)曲波變換去噪數(shù)據(jù)的PSNR值隨閾值控制系數(shù)δ的變化如圖1(d)所示.隨著閾值控制系數(shù)δ取值由小變大,去噪數(shù)據(jù)PSNR值趨于極值后逐漸降至原噪聲數(shù)據(jù)PSNR值水平,其中含隨機(jī)噪聲數(shù)據(jù)去噪結(jié)果PSNR值在閾值控制系數(shù)δ= 0.15 (與原始數(shù)據(jù)噪聲水平一致)附近取得PSNR極值(27.33 dB),含相關(guān)噪聲數(shù)據(jù)去噪結(jié)果對(duì)應(yīng)在δ= 0.35附近取得PSNR極值(23.27 dB).上述結(jié)果分析表明,采用閾值函數(shù)曲波變換去噪效果的好壞取決于閾值控制系數(shù)(原始數(shù)據(jù)噪聲水平)取值是否合理.同時(shí),含隨機(jī)噪聲與相關(guān)噪聲數(shù)據(jù)采用閾值函數(shù)曲波變換去噪的合理閾值控制系數(shù)規(guī)律不同,前者閾值控制系數(shù)取噪聲水平值(δ= 0.10—0.20),而后者取略高于閾值控制系數(shù)值(δ= 0.3—0.4).

圖1 含噪聲模擬合成探地雷達(dá)數(shù)據(jù)的曲波變換閾值去噪結(jié)果 (a)原始合成數(shù)據(jù);(b)含隨機(jī)噪聲數(shù)據(jù),PSNR = 17 dB;(c)含相關(guān)噪聲數(shù)據(jù),PSNR = 15.51 dB;(d)去噪數(shù)據(jù)PSNR值隨閾值控制系數(shù)δ的變化Fig.1.Denoised results for the synthetic ground penetrating radar (GPR) data with random and coherent noise using traditional curvelet transform：(a) Original GPR data;(b) data with random noise,PSNR = 17 dB;(c) data with coherent noise,PSNR = 15.51 dB;(d) theδvalue vs.PSNRvalue curves .

針對(duì)含噪模擬數(shù)據(jù)閾值函數(shù)曲波變換去噪問(wèn)題,預(yù)先獲取的準(zhǔn)確閾值控制系數(shù)范圍可有效用于含噪數(shù)據(jù)處理.但實(shí)測(cè)數(shù)據(jù)常常被不同類(lèi)型、不同強(qiáng)度的噪聲源污染,合理的閾值控制系數(shù)范圍難以確定.為了有效使用閾值函數(shù)曲波變換去噪算法,需要通過(guò)特定方式人為估計(jì)閾值控制系數(shù)范圍,如L2標(biāo)準(zhǔn)方差算法.圖2所示為含噪聲數(shù)據(jù)(圖1(b)和圖1(c))采用L2標(biāo)準(zhǔn)方差估計(jì)方法確定閾值控制系數(shù)為0.08時(shí)的曲波變換去噪結(jié)果.相比原始含噪數(shù)據(jù)(圖1(b)和圖1(c)),閾值函數(shù)曲波變換去噪清除了部分噪聲成分信號(hào),但矩形目標(biāo)模型引起的有效反射波信號(hào)僅依稀可辯,無(wú)法有效用于全波形反演成像計(jì)算及后續(xù)推斷解譯過(guò)程(見(jiàn)圖2).

3.2 統(tǒng)計(jì)量自適應(yīng)閾值曲波變換去噪

為解決傳統(tǒng)曲波變換去噪需要估計(jì)合理閾值函數(shù)范圍的問(wèn)題,采用統(tǒng)計(jì)量自適應(yīng)閾值曲波變換去噪方法對(duì)圖1所示含隨機(jī)噪聲、相關(guān)噪聲數(shù)據(jù)進(jìn)行處理.去噪結(jié)果如圖3所示,基于統(tǒng)計(jì)量自適應(yīng)閾值曲波變換去噪方法不僅可有效衰減隨機(jī)噪聲(PSNR值提高8.3 dB)、相關(guān)噪聲(PSNR值提高6.41 dB),同時(shí)較好地還原剖面中有效反射波信號(hào)時(shí)空域分布特征.在噪聲處理過(guò)程中,本文提出的去噪算法無(wú)需估計(jì)閾值函數(shù)范圍,而是采用(8)式計(jì)算統(tǒng)計(jì)量自適應(yīng)閾值權(quán)重,達(dá)到了保真去噪的目的.其中,隨機(jī)噪聲數(shù)據(jù)去噪效果優(yōu)于相關(guān)噪聲數(shù)據(jù)去噪效果,究其緣由,低信噪比相關(guān)噪聲數(shù)據(jù)(PSNR= 15.51 dB)導(dǎo)致噪聲信號(hào)相關(guān)性統(tǒng)計(jì)量計(jì)算誤差,殘留部分相關(guān)性較強(qiáng)的噪聲信號(hào).總體上,統(tǒng)計(jì)量自適應(yīng)閾值曲波變換去噪結(jié)果可有效地用于探地雷達(dá)數(shù)據(jù)全波形偏移成像處理及后續(xù)解釋推斷過(guò)程.

圖2 δ = 0.08時(shí)含噪聲數(shù)據(jù)(圖1(b)和圖1(c))曲波變換去噪結(jié)果 (a)隨機(jī)噪聲數(shù)據(jù)去噪結(jié)果,PSNR = 20.23 dB,提高3.23 dB;(b)相關(guān)噪聲數(shù)據(jù)去噪結(jié)果,PSNR = 16.75 dB,提高1.24 dBFig.2.Denoised results for the synthetic GPR data (Fig.1.(b) and Fig.1.(c)) withδ = 0.08 using traditional curvelet transform：(a) Result for random noise,PSNR = 20.23 dB;(b) result for coherent noise,PSNR = 16.75 dB.

圖3 含噪聲數(shù)據(jù)(圖1(b)和圖1(c))的統(tǒng)計(jì)量自適應(yīng)閾值曲波變換去噪結(jié)果 (a)隨機(jī)噪聲數(shù)據(jù)去噪結(jié)果,PSNR = 25.3 dB,提高8.3 dB;(b)相關(guān)噪聲數(shù)據(jù)去噪結(jié)果,PSNR = 21.92 dB,提高6.41 dBFig.3.Denoised results for the synthetic GPR data (Fig.1.(b) and Fig.1.(c)) using curvelet transform with statistical self-adaption：(a) Result for random noise,PSNR = 25.3 dB;(b) result for coherent noise,PSNR = 21.92 dB.

進(jìn)一步,設(shè)計(jì)雙矩形目標(biāo)體規(guī)模為1 m × 1 m,水平位置分別為—2至—3 m及2至3 m,埋深范圍分別為2至3 m及3至4 m;其他物性參數(shù)及探地雷達(dá)觀測(cè)系統(tǒng)與前述單個(gè)矩形目標(biāo)模型一致.數(shù)值模擬的雙矩形目標(biāo)體探地雷達(dá)時(shí)間剖面如圖4(a)所示,加入隨機(jī)噪聲及相關(guān)噪聲數(shù)據(jù)如圖4(b)和圖4(c)所示.含隨機(jī)噪聲數(shù)據(jù)(PSNR= 16.04 dB)及相關(guān)噪聲數(shù)據(jù)(PSNR= 15.7 dB)完全淹沒(méi)雙矩形目標(biāo)體引起的有效反射波信號(hào),無(wú)法用于雙矩形目標(biāo)體時(shí)空域分布推斷解釋.

圖4(d)和圖4(e)為采用本文提出的去噪算法處理結(jié)果.由去噪結(jié)果可見(jiàn),含隨機(jī)噪聲時(shí)間剖面內(nèi)雙矩形目標(biāo)體引起復(fù)雜有效反射波信號(hào)被完全恢復(fù),PSNR值較原始含噪數(shù)據(jù)提高7.93 dB;含相關(guān)噪聲數(shù)據(jù)處理結(jié)果基本重建目標(biāo)體引起有效反射波信號(hào)時(shí)空分布,但仍殘余部分相關(guān)噪聲成分,PSNR值較原始含噪數(shù)據(jù)提高6.65 dB.抽取圖4所示去噪結(jié)果第50道數(shù)據(jù),繪制原始無(wú)噪聲數(shù)據(jù)、含噪數(shù)據(jù)及去噪數(shù)據(jù)信號(hào)局部細(xì)節(jié)對(duì)比曲線,如圖5所示.可見(jiàn),本文提出的去噪算法有效清除了全局椒鹽式隨機(jī)噪聲分布(圖5(a))、可有效分辨相關(guān)噪聲環(huán)境下微弱有效反射波信號(hào)(圖5(b)),去噪數(shù)據(jù)曲線與原始無(wú)噪聲數(shù)據(jù)曲線吻合較好,驗(yàn)證了本文提出的去噪算法有效性及可行性.

圖4 雙矩形目標(biāo)模型含噪聲數(shù)據(jù)統(tǒng)計(jì)量自適應(yīng)閾值曲波變換去噪結(jié)果 (a)原始合成數(shù)據(jù);(b)含隨機(jī)噪聲數(shù)據(jù),PSNR = 16.04 dB;(c)含相關(guān)噪聲數(shù)據(jù),PSNR = 15.7 dB;(d)隨機(jī)噪聲數(shù)據(jù)去噪結(jié)果,PSNR = 23.97 dB;(e)相關(guān)噪聲數(shù)據(jù)去噪結(jié)果,PSNR = 21.05 dBFig.4.Denoised results for the synthetic GPR data of complex model with random and coherent noise using curvelet transform with statistical self-adaption：(a) Original GPR data;(b) data with random noise,PSNR = 16.04 dB;(c) data with coherent noise,PSNR = 15.7 dB;(d) result for random noise,PSNR = 23.97 dB;(e) result for coherent noise,PSNR = 21.05 dB.

數(shù)值模擬實(shí)驗(yàn)結(jié)果表明,傳統(tǒng)曲波變換去噪效果依賴于閾值控制系數(shù)選取范圍是否符合含噪數(shù)據(jù)噪聲水平范圍,而不同信噪比含噪數(shù)據(jù)噪聲水平各異,因此采用傳統(tǒng)閾值函數(shù)的曲波去噪算法抗噪性能與估計(jì)閾值控制系數(shù)范圍準(zhǔn)確度具有直接關(guān)系.對(duì)于含有弱信號(hào)成分探地雷達(dá)數(shù)據(jù)的去噪處理,普遍應(yīng)用于傳統(tǒng)曲波變換去噪算法理論與實(shí)踐研究的L2標(biāo)準(zhǔn)差估計(jì)的噪聲水平導(dǎo)致偽異常信號(hào)問(wèn)題.本文提出的統(tǒng)計(jì)量自適應(yīng)閾值曲波變換去噪算法可有效確定含噪數(shù)據(jù)噪聲水平范圍,在單一以及多個(gè)異常信號(hào)疊加的探地雷達(dá)數(shù)據(jù)處理方面顯示了較好的應(yīng)用效果.需要指出的是,本文考慮數(shù)據(jù)噪聲水平為15%,遠(yuǎn)遠(yuǎn)高于地球物理勘探數(shù)據(jù)采集質(zhì)量控制誤差水平5%,數(shù)值算例結(jié)果可見(jiàn),含隨機(jī)噪聲、相關(guān)噪聲數(shù)據(jù)處理結(jié)果可有效地用于后續(xù)成像及解譯;而高于15%含噪數(shù)據(jù)通常視為無(wú)效數(shù)據(jù),其相應(yīng)去噪處理分析研究不在本文討論范圍.

4 實(shí)測(cè)算例

圖6(a)所示為某地探地雷達(dá)實(shí)測(cè)某測(cè)線原始時(shí)間剖面圖,由于受場(chǎng)地環(huán)境影響,椒鹽式隨機(jī)噪聲遍布整個(gè)剖面數(shù)據(jù);部分采集道數(shù)據(jù)受局部不均勻體影響,出現(xiàn)相鄰觀測(cè)道數(shù)據(jù)幅值畸變.原始數(shù)據(jù)剖面淺部(100—250時(shí)間采樣點(diǎn))介質(zhì)的反射波能量非常微弱,依稀可見(jiàn)幾處繞射波同相軸,但受噪聲淹沒(méi),雙曲線特征不明顯;在250—300時(shí)間采樣范圍內(nèi)出現(xiàn)較強(qiáng)幅值的似平行反射波組能量,同相軸斷斷續(xù)續(xù)并不連續(xù);強(qiáng)幅值反射波組下部(300—400時(shí)間采樣點(diǎn))出現(xiàn)較弱的零散反射波組,難以識(shí)別強(qiáng)幅值反射波組下邊界.

利用上述傳統(tǒng)曲波變換去噪處理方法對(duì)圖6(a)進(jìn)行數(shù)據(jù)處理,通過(guò)L2標(biāo)準(zhǔn)方差估計(jì)閾值范圍,去噪獲得的探地雷達(dá)時(shí)間剖面如圖6(b)所示.由去噪結(jié)果可見(jiàn),部分地下反射波組能量信號(hào)得以加強(qiáng),清除了部分隨機(jī)噪聲信號(hào).但原始數(shù)據(jù)內(nèi)有效繞射波雙曲線特征、錯(cuò)動(dòng)反射波組依然無(wú)法有效識(shí)別,去噪數(shù)據(jù)信噪比提高并不顯著.究其原因,實(shí)測(cè)數(shù)據(jù)中隨機(jī)噪聲、相關(guān)噪聲及其他類(lèi)型噪聲源摻雜,同時(shí),原始數(shù)據(jù)包含幅值強(qiáng)度各異的反射波組,L2標(biāo)準(zhǔn)方差估計(jì)閾值范圍無(wú)法有效涵蓋上述噪聲特點(diǎn)的數(shù)據(jù)去噪范圍.人工試錯(cuò)法確定最佳閾值范圍或能有效清除噪聲,但最佳去噪效果評(píng)價(jià)缺乏科學(xué)依據(jù).因此,基于傳統(tǒng)閾值函數(shù)曲波變換去噪數(shù)據(jù)處理后的探地雷達(dá)時(shí)間剖面進(jìn)行分析解釋,無(wú)法保證資料解譯的準(zhǔn)確度.

采用本文提出的統(tǒng)計(jì)量自適應(yīng)閾值曲波變換去噪算法對(duì)圖6(a)原始數(shù)據(jù)進(jìn)行處理,去噪結(jié)果如圖6(c)所示.由于隨機(jī)噪聲源在曲波變換系數(shù)尺度、方向上并不具備統(tǒng)計(jì)相關(guān)特性;相關(guān)噪聲源變換系數(shù)在尺度、方向上雖具備一定相關(guān)性,但相對(duì)有效反射波組信號(hào)相關(guān)性仍然較小,因而統(tǒng)計(jì)量自適應(yīng)閾值可有效分辨反射波組能量,去除噪聲信號(hào),達(dá)到保真去噪的目的.由圖6(c)可見(jiàn),地下反射波能量得到顯著增強(qiáng),特別是淺部多個(gè)繞射波雙曲線同相軸特征的信噪比得到有效提高;中深度多層介質(zhì)的反射組同相軸連續(xù)、分界面明顯;且分布于強(qiáng)幅值反射波組下部的弱幅值反射波組清晰可辯.經(jīng)過(guò)統(tǒng)計(jì)量自適應(yīng)閾值曲波變換數(shù)據(jù)處理后的探地雷達(dá)時(shí)間剖面進(jìn)行分析解釋,有助于資料解譯的準(zhǔn)確度.

抽取圖6所示實(shí)測(cè)探地雷達(dá)數(shù)據(jù)第200道數(shù)據(jù),繪制傳統(tǒng)曲波變換去噪結(jié)果與本文提出的去噪結(jié)果局部細(xì)節(jié)對(duì)比曲線,如圖7所示.由曲線對(duì)比圖可見(jiàn),在0—100時(shí)間采樣范圍內(nèi)地面回波、弱幅值繞射波(100—150時(shí)間采樣)、強(qiáng)幅值繞射波(150—230時(shí)間采樣)及強(qiáng)幅值反射波組(230—350時(shí)間采樣)受到隨機(jī)噪聲及相關(guān)噪聲污染,在第350時(shí)間采樣范圍以外的弱反射波信號(hào)受到較強(qiáng)幅值的上述噪聲源淹沒(méi).相比傳統(tǒng)曲波變換去噪結(jié)果,本文提出的去噪算法可有效恢復(fù)弱幅值繞射波、強(qiáng)幅值繞射波及強(qiáng)幅值反射波組的同相軸特征,同時(shí)可清晰分辨弱幅值反射波信號(hào),驗(yàn)證了本文提出的去噪算法在實(shí)測(cè)探地雷達(dá)時(shí)間剖面數(shù)據(jù)處理應(yīng)用中的可行性及有效性.

5 結(jié) 論

1)針對(duì)傳統(tǒng)曲波變換探地雷達(dá)數(shù)據(jù)去噪需估計(jì)閾值函數(shù)問(wèn)題,對(duì)探地雷達(dá)數(shù)據(jù)做二維曲波正變換獲取不同尺度、不同選擇角度下的系數(shù)分布;對(duì)每一道探地雷達(dá)數(shù)據(jù)在不同尺度、不同角度下計(jì)算高階相關(guān)統(tǒng)計(jì)量,根據(jù)統(tǒng)計(jì)量分布明確閾值分布權(quán)重,開(kāi)展了統(tǒng)計(jì)量自適應(yīng)函數(shù)算法研究,實(shí)現(xiàn)無(wú)需估計(jì)曲波變換閾值函數(shù)及其所屬系數(shù)范圍,而是通過(guò)統(tǒng)計(jì)理論明確目標(biāo)體探地雷達(dá)有效信息所屬的尺度及旋轉(zhuǎn)角度范圍,從而構(gòu)建了統(tǒng)計(jì)量自適應(yīng)閾值函數(shù)曲波變換去噪算法.

2)基于塊狀復(fù)數(shù)域閾值函數(shù)理論,設(shè)計(jì)簡(jiǎn)單矩形及埋深不同雙矩形模型,并合成隨機(jī)噪聲、相關(guān)噪聲合成探地雷達(dá)數(shù)據(jù).由傳統(tǒng)閾值函數(shù)曲波變換去噪分析表明：塊狀復(fù)數(shù)域閾值范圍的計(jì)算需給定閾值控制系數(shù),去噪效果的優(yōu)劣取決于與含噪數(shù)據(jù)是否匹配的閾值控制系數(shù);合成含噪數(shù)據(jù)可通過(guò)給定的估計(jì)閾值函數(shù)實(shí)現(xiàn)高效去噪,但實(shí)測(cè)含噪數(shù)據(jù)無(wú)法準(zhǔn)確估計(jì)閾值函數(shù),因而難以實(shí)現(xiàn)高效保真的去噪目的.

3)對(duì)復(fù)雜噪聲源環(huán)境下采集的探地雷達(dá)實(shí)測(cè)數(shù)據(jù)進(jìn)行去噪分析處理,提取了微弱幅值繞射波雙曲線同相軸特征,恢復(fù)了中深部不同幅值似平行非連續(xù)性反射波組及弱幅值錯(cuò)動(dòng)反射波組,獲得了相應(yīng)去噪結(jié)果.將相應(yīng)結(jié)果與傳統(tǒng)閾值函數(shù)曲波變換去噪結(jié)果對(duì)比,驗(yàn)證了統(tǒng)計(jì)量自適應(yīng)閾值曲波變換去噪算法在復(fù)雜噪聲背景下探地雷達(dá)弱信號(hào)提取方面的獨(dú)特優(yōu)勢(shì),有助于對(duì)探地雷達(dá)探測(cè)資料進(jìn)行可靠、準(zhǔn)確的解譯.