基于二維小波能量閾值的重磁數(shù)據(jù)去噪算法

索 奎,呂曉春,張貴賓,賈正元

(1.華北水利水電大學(xué)地球科學(xué)與工程學(xué)院,河南鄭州 450046;2.中國地質(zhì)大學(xué)(北京)地球物理與信息技術(shù)學(xué)院,北京 100083)

隨著儀器和相關(guān)理論的發(fā)展,實(shí)測重磁數(shù)據(jù)精度越來越高,數(shù)據(jù)量迅速增加。由于地質(zhì)構(gòu)造的復(fù)雜性和儀器的非穩(wěn)定性,數(shù)據(jù)包含了人文噪聲、系統(tǒng)噪聲以及數(shù)據(jù)處理引入的誤差等,這類噪聲具有非平穩(wěn)和非線性特征[1],影響數(shù)據(jù)處理和反演的精度,甚至掩蓋有效信息。線性去噪方法去除非線性噪聲的效果不佳;小波變換能集中有效信號的能量,使噪聲在小波域具有隨機(jī)性,通過對小波系數(shù)施加某種非線性操作實(shí)現(xiàn)去噪,在地球物理數(shù)據(jù)處理領(lǐng)域應(yīng)用廣泛[2-6]。但小波去噪主要面臨小波基函數(shù)選取[7-10],小波信噪分離階數(shù)選取[11-13]和去噪閾值選取[13-15]3個關(guān)鍵問題。選擇分離階數(shù)主要有基于信號特征[16]和基于噪聲特征的方法;閾值選擇方法也多種多樣。雖然各類方法針對不同問題提出了解決方案[17],但處理實(shí)際數(shù)據(jù)時多數(shù)情況仍依賴經(jīng)驗(yàn)選取參數(shù)。由于模型試驗(yàn)參照標(biāo)準(zhǔn)確定,可以采用各種評價標(biāo)準(zhǔn)去評價去噪效果[18],解決上述3個問題;而實(shí)測數(shù)據(jù)沒有參照標(biāo)準(zhǔn),多數(shù)評價標(biāo)準(zhǔn)只能估計或作為參考,導(dǎo)致去噪?yún)?shù)選擇標(biāo)準(zhǔn)不同,結(jié)果因人而異。Lyrio等[19-20]提出了一維小波變換中對3個關(guān)鍵問題的解決方案;索奎[21]將其擴(kuò)展至二維情形,利用二維小波將重磁數(shù)據(jù)分為3個高頻部分分別進(jìn)行處理,從而得到更精確的結(jié)果。筆者重點(diǎn)研究基于能量閾值的二維小波去噪中的3個關(guān)鍵問題,分別針對加性噪聲和乘性噪聲開展模型試驗(yàn)。

1 小波去噪原理與方法

小波變換可將信號分解為不同頻段,再利用其多分辨率的性質(zhì)來研究不同尺度下每個頻段的成分。與傅里葉變換相比,小波變換在時域和頻域都具有良好的局部性質(zhì),能夠聚焦研究對象的細(xì)節(jié)。

1.1 二維小波變換原理

(1)

式中,j為分解尺度。

對于重磁面積數(shù)據(jù)Δg(x,y),也可以通過二維小波將總場分解為局部場(細(xì)節(jié))和區(qū)域場(近似),n階分解的表達(dá)式為

Δg(x,y)=A0f(x,y)=Anf(x,y)+

(2)

二維離散小波變換先利用高通濾波器和低通濾波器沿著水平(行)方向?qū)?shù)據(jù)進(jìn)行一維分解,然后利用同樣的濾波器沿著垂直(列)方向?qū)π蟹纸饨Y(jié)果進(jìn)行分解,二維小波3層分解如圖1所示。經(jīng)過3層分解,得到4種信號:水平、垂直方向均為低頻的信號LL3,水平方向低頻、垂直方向高頻的信號LH1、LH2、LH3,水平方向高頻、垂直方向低頻的信號HL1、HL2、HL3以及水平、垂直方向均為高頻的信號HH1、HH2、HH3。其中LL3稱為逼近部分,其余3種信號稱為細(xì)節(jié)部分。

圖1 二維小波3層分解示意圖Fig.1 Three-level decomposition of two-dimensional wavelet

1.2 小波去噪基本原理

小波變換的能量壓縮特性能夠使信號變換后的大部分能量集中在少部分小波系數(shù)幅值較高的區(qū)域,即低頻系數(shù),這些區(qū)域一般是圖像邊緣,表示圖像的主要特征。大多數(shù)噪聲具有較高的頻率,經(jīng)過小波變換后的系數(shù)幅值很低。

小波去噪的基本原理是將含噪信號變換后的高頻系數(shù)與某個閾值相比,然后按照特定的閾值函數(shù)進(jìn)行修正。最后將修正后的高頻系數(shù)和低頻系數(shù)進(jìn)行重構(gòu)(逆變換),獲得去噪后的信號?；谀芰块撝档亩S小波去噪方法仍基于上述原理。若考慮野外采集時儀器和環(huán)境影響因素,加性噪聲較為符合重磁數(shù)據(jù)的特點(diǎn);若考慮各項(xiàng)改正引起的噪聲,則乘性噪聲更符合數(shù)據(jù)特點(diǎn)。

本文中假設(shè)噪聲符合高斯模型,以含加性噪聲的理論模型為例進(jìn)行闡述。設(shè)置包含了4個密度異常體的組合模型,參數(shù)見表1。

表1 理論模型參數(shù)Table 1 Synthetic model parameters

地面測網(wǎng)范圍為1 023 m×1 023 m,測點(diǎn)點(diǎn)距為1 m,正演獲得了1 024×1 024個重力數(shù)據(jù),數(shù)據(jù)范圍為(0.004 7～1.191 8)×10-5m/s2,如圖2(a)所示,4個異常體邊界較為清晰,其中2號異常體體積小,異常幅值最高,理論值等值線平滑。添加了均值為0,標(biāo)準(zhǔn)差為0.05×10-5m/s2的高斯噪聲,噪聲范圍為(-0.234 6～0.227 3)×10-5m/s2,添加噪聲后的數(shù)據(jù)如圖2(b)所示,仍可以看出基本異常特征,但2號異常不顯著,且等值線波動較大。

圖2 加性噪聲模型理論重力異常和含噪聲重力異常Fig.2 Synthetic gravity anomalies and gravity anomalies with Gaussian noise

1.3 基于信噪比的基函數(shù)選擇方法

小波變換前應(yīng)選定小波基函數(shù)并確定分離階數(shù)。目前小波基函數(shù)有上百種,合適的小波基函數(shù)是去噪的基礎(chǔ)。本次研究的目標(biāo)是基于能量的重力和磁場的噪聲去除,要求滿足含噪信號能量等于小波變換后系數(shù)能量的特性[20],即

(3)

式中,ω(j,k)為尺度j位置k處的小波系數(shù);g(x)為空間域中位置x處的數(shù)據(jù)。

滿足該特性的小波基函數(shù)必須具有正交性,常見小波基函數(shù)有Haar小波、DB(Daubechies)小波、COIF(Coiflets)小波和SYM(Symlets)小波,本文中選取的小波基函數(shù)是來源于上述幾個小波族的。

多種小波基函數(shù)符合特性,因此需要選定最佳的小波基函數(shù)。選擇最佳小波基函數(shù)基于信噪比,定義一個信噪比(RSN)公式[20]:

(4)

式中,ej為分解階數(shù)j的含噪信號的高頻能量;γ為最佳信噪分離階數(shù);J為最大分解階數(shù)。

信噪比越大,該小波基函數(shù)去噪效果越好,本文中根據(jù)式(4)選擇DB5作為小波基函數(shù),其信噪比為1.595 4,γ為5。

1.4 信噪分離階數(shù)選擇方法

在實(shí)際數(shù)據(jù)計算時,雖然無法得知含噪信號中的噪聲和有效信號的能量,但該含噪信號用同一個小波基函數(shù)進(jìn)行分解時,其能量(式(3))等于噪聲能量和有效信號能量之和[19]。

在不同的分離階數(shù)時,噪聲能量和有效信號能量在高頻總能量的占比發(fā)生變化(圖3),總體趨勢是:分解階數(shù)較低時,高頻能量以噪聲能量為主;隨著分解階數(shù)升高,高頻能量以有效信號能量為主,二者“交叉點(diǎn)”處的階數(shù)為最佳分離階數(shù),此處也為含噪信號的能量最低點(diǎn)[19-21]?；谠撛?進(jìn)行信噪分離可以對含噪信號進(jìn)行高階次分解后,選擇最佳階數(shù)作為分界線,低于該階數(shù)的小波系數(shù)全部舍去,即舍去了以噪聲為主的部分,保留了有效信號為主的部分。

二維小波分解會得到3個高頻信號,其中有兩個信號是經(jīng)過高通和低通濾波的結(jié)果,即LH和HL分量,二者含噪成分相當(dāng);一個信號是經(jīng)過了兩次高通濾波的結(jié)果,即HH分量,噪聲成分更多,信噪分離階數(shù)更高[21]。因此分別根據(jù)3個能量曲線確定的信噪分離階數(shù)會有所不同,但由于對高頻信號進(jìn)行了更加細(xì)致的分解并分別進(jìn)行了處理,因此對信噪分離更加準(zhǔn)確。

圖3為圖2(b)數(shù)據(jù)進(jìn)行二維小波變換后3個高頻分量在1～7階分解情形下的有效信號和理論噪聲能量分布變化趨勢。圖3(a)為含噪信號高頻能量和(LH+HL+HH)隨階數(shù)的變化,圖3(b)～(d)分別是LH、HL和HH分量的能量變化。

以圖3(a)為例,含噪信號能量是有效信號和噪聲之和,在1～4階分解時,其能量逐漸減小,在5～7階分解時又迅速增大;有效信號的能量是逐步增加的,噪聲則相應(yīng)減小。在1～4階時,噪聲能量顯著大于有效信號能量,表明此時含噪信號中以噪聲為主,隨著分解階數(shù)繼續(xù)增加,有效信號能量大于噪聲能量,即有效信號占主要成分。從圖3(b)和(c)可以看出,LH和HL分量特征較為一致,含噪數(shù)據(jù)能量最低的階數(shù)為第4階,HH分量噪聲更多,能量最低階數(shù)為第5階(圖3(d)),此時需要分別對1～4階和1～5階的小波系數(shù)置零,處理5～7階和6～7階小波系數(shù)即可。顯然,分別對3種高頻分量的數(shù)據(jù)進(jìn)行去噪,能夠針對不同的高頻分量分別采用相應(yīng)的閾值,獲得更精確的處理結(jié)果。

圖3 合成模型重力信號和噪聲在小波域不同分解階數(shù)能量分布Fig.3 Gravity synthetic signal and noise of model energy distribution of different scales in wavelet domain

1.5 基于能量的去噪閾值選擇方法

采用確定的小波基函數(shù)和分解階數(shù)對含噪數(shù)據(jù)進(jìn)行分解,再對細(xì)節(jié)部分系數(shù)進(jìn)行閾值去噪。目前有多種閾值選擇方法,如:通用閾值、基于Stein無偏似然估計的SureShrink閾值以及假設(shè)信號概率分布函數(shù)服從Beyes準(zhǔn)則的BayesShrink閾值,除此之外還有多種策略的閾值選取算法,適用于各種情況。

本文中采用基于能量的閾值選取方法。由于信號變換后的大部分能量集中在少部分小波系數(shù)幅值較高、表示圖像主要特征的區(qū)域,理論上只需少部分高幅值系數(shù)即可恢復(fù)大部分有效信號。選定某一個閾值,將幅值系數(shù)從大到小進(jìn)行能量累加,達(dá)到閾值即認(rèn)為恢復(fù)了有效信號。

以圖2(b)所示含噪信號為例,將信號進(jìn)行二維小波變換7階分解后,將小波系數(shù)能量按幅值從大到小累加,能量變化曲線如圖4(a)所示?？梢钥闯?對于HL和LH分量,初期僅需20～40個高幅值系數(shù)的能量累加,總能量迅速上升到較高水平,后面雖然仍有大量系數(shù),但能量很少,表明信號的主要能量都集中在了高幅值系數(shù)中。圖中示例恢復(fù)有效信號分別使用了249和253個HL和LH系數(shù);HH分量使用了102個系數(shù),后續(xù)的能量仍然較大,但多是以噪聲為主。

能量閾值并非越高越好,這是因?yàn)楹胄盘柕哪芰恐邪擞行盘柡驮肼暤哪芰?因此若選定過高的閾值既會增加計算量,也會使得恢復(fù)的信號中包含過多的噪聲[21]。

圖4 能量閾值與均方誤差關(guān)系Fig.4 Relationship between energy threshold and mean square error

以理論信號為標(biāo)準(zhǔn),保持近似小波系數(shù)不變,用不同數(shù)量的細(xì)節(jié)小波系數(shù)恢復(fù)的信號與理論信號進(jìn)行比對,重構(gòu)信號與有效信號之間的均方誤差如圖4(b)所示:重構(gòu)均方誤差開始為5.085 6×10-5m/s2,隨著重構(gòu)小波系數(shù)的增加,細(xì)節(jié)信息不斷增加,誤差急劇減小到0.008 8×10-5m/s2,表明此時恢復(fù)的以有效信號為主,隨著系數(shù)個數(shù)的繼續(xù)增加,誤差開始逐漸上升,直至2.588 3×10-5m/s2,表明此時恢復(fù)的信號中包含了噪聲,因此與有效信號相比誤差反而增大了。

從圖4(b)可以看出理論最佳閾值應(yīng)該在誤差最小處(左側(cè)虛線),但處理實(shí)際數(shù)據(jù)時由于不知道有效信號的準(zhǔn)確值,無法精確確定該位置,根據(jù)式(4)確定的最佳信噪比中的能量作為閾值來確定參與重構(gòu)的細(xì)節(jié)系數(shù)個數(shù),其位置如圖4(b)中的右側(cè)虛線,此時誤差并非最小,但與最佳值很接近。

基于能量閾值的二維小波去噪方法主要步驟可以總結(jié)為:①依據(jù)公式選擇信噪比最大的小波基函數(shù);②確定最佳分離階數(shù),計算不同階數(shù)的能量,能量最低點(diǎn)即為最佳分離階數(shù);③確定能量閾值,最佳分離階數(shù)及以上的能量和即為該閾值。利用本文方法對數(shù)據(jù)進(jìn)行去噪,既可以得到有效信號利于進(jìn)一步數(shù)據(jù)處理,也可以更準(zhǔn)確估計噪聲水平,為提高三維反演結(jié)果的精度提供依據(jù)[22]。

2 模型試驗(yàn)

2.1 加性噪聲模型試驗(yàn)

為了進(jìn)一步對比,本文中還對添加了標(biāo)準(zhǔn)差為0.10×10-5和0.2×10-5m/s2高斯噪聲的信號進(jìn)行去噪,以對比不同噪聲水平的去噪結(jié)果。

針對含噪信號,本文中選擇巴特沃斯低通濾波器、常規(guī)二維小波和本文方法進(jìn)行去噪,去噪結(jié)果如圖5所示,均方誤差見表2。為了對比方便,圖5插值方法及參數(shù)相同,繪圖參數(shù)設(shè)置相同,繪圖設(shè)置保持一致。

圖5(a)～(c)為巴特沃斯低通濾波器濾波結(jié)果,均方誤差分別為0.027 5×10-5、0.813 0×10-5和0.992 5×10-5m/s2,呈顯著上升趨勢。從圖5可以看出,圖5(a)具有較好的去噪效果,均方誤差與小波去噪相差不大,但在1號異常體處,原本近似方形的異常等值線在去噪后變成了圓形,圖5(b)和(c)也有相同情況,且均方誤差與小波去噪結(jié)果已經(jīng)不在同一個數(shù)量級,這表明了巴特沃斯低通濾波器局限之處。

圖5(d)～(f)是二維小波去噪結(jié)果。本文中選用DB4小波,全局硬閾值去噪,此時的均方誤差最低,分別為0.014 1×10-5、0.034 2×10-5和0.098 9×10-5m/s2,均顯著低于巴特沃斯濾波結(jié)果。

本文中根據(jù)式(4)選擇DB5作為小波基函數(shù),最佳分離階數(shù)為5、5和6,去噪結(jié)果如圖5(g)～(i)?？梢钥闯龆鄶?shù)等值線曲線平滑,個別等值線有一些波動,在異常體1的附近波動最大,這是因?yàn)榇颂幍奶荻戎底畲?等值線最密集,總體看對異常識別沒有影響。查看計算得到的噪聲,符合高斯白噪聲的特征,標(biāo)準(zhǔn)差為0.05×10-5m/s2,與理論值一致,這表明噪聲計算較為準(zhǔn)確。從表2可以看到3種情形本文方法均方誤差均為最小,且隨著噪聲水平的增加,增長最慢,表明了方法的穩(wěn)定性。

該試驗(yàn)驗(yàn)證了二維小波閾值去噪算法具有更小的均方誤差,去除不同水平的加性噪聲具有穩(wěn)定性和準(zhǔn)確性。

圖5 模型試驗(yàn)加性噪聲去噪結(jié)果Fig.5 Denoising result of synthetic data with additive noise

表2 不同方法模型試驗(yàn)去噪結(jié)果均方誤差對比Table 2 Comparison of mean square error of denoising results of different methods

2.2 乘性噪聲模型試驗(yàn)

乘性噪聲模型與加性噪聲模型一致,添加了3%的乘性高斯隨機(jī)噪聲,噪聲范圍(-0.135～0.135)×10-5m/s2,根據(jù)式(4)選擇DB4作為小波基函數(shù),對含噪數(shù)據(jù)進(jìn)行三分量閾值去噪。最佳分離階數(shù)為4、4和5,去噪后結(jié)果如圖6。HL分量有效系數(shù)227個,LH分量有效系數(shù)263個,HH分量有效系數(shù)128個,噪聲均方誤差為0.007 3×10-5m/s2。

從圖6(a)中可以看出,異常幅值越高,噪聲幅值越高,符合乘性噪聲的特征;去噪后如圖6(b)所示,總體上較好地恢復(fù)了原始數(shù)據(jù),大部分等值線曲線平滑,等值線0.05有一些波動。圖6(c)和(d)分別為理論噪聲和本文計算的噪聲,異常體1處的噪聲差別較大,此外在異常體2～4的邊緣區(qū)域噪聲差異也相對明顯,主要集中在梯度值較大的區(qū)域,其他區(qū)域則較為“干凈”,沒有因?yàn)槿ピ胨惴◣砣藶樵肼暋?/p>

圖6 模型試驗(yàn)乘性噪聲去噪結(jié)果對比Fig.6 Denoising result of synthetic data with multiplicative noise

3 某區(qū)域?qū)崪y磁異常去噪

本文中選擇的實(shí)測磁異常數(shù)據(jù)位于包含多個隆起、凹陷的構(gòu)造單元匯聚區(qū),面積約1 100 km2,磁異常變化范圍大,有異常變化劇烈區(qū)和相對平靜區(qū),部分區(qū)域有人類活動引起的人文噪聲。

磁異常數(shù)據(jù)分辨率為250 m×500 m,磁異常值范圍為-654.77～1 477.21 nT,原始數(shù)據(jù)如圖7(a)所示,可以看出,圖右側(cè)有本區(qū)域面積最大、幅值最高的磁異常區(qū),中部和左側(cè)有部分小范圍磁異常,磁噪聲使部分區(qū)域形成了小范圍假異常區(qū),等值線較為復(fù)雜。本文中分別采用巴特沃斯低通濾波、二維小波去噪、本文去噪方法3種方法分別進(jìn)行去噪處理。3種方法去噪結(jié)果如圖7(b)～(d)所示,為方便對比,采用相同的繪圖參數(shù)。

圖7(b)為采用巴特沃斯低通濾波器進(jìn)行低通濾波結(jié)果,參數(shù)為8階,截止波長1 000 m?？梢钥闯?線性濾波器的處理后等值線變得圓滑,高頻噪聲得到了較好的去除,數(shù)據(jù)范圍-496.72～1 475.98 nT,表明了對顯著的高值有較好的保幅效果,但對異常范圍小的負(fù)值則有較為明顯的“損失”,在箭頭指示區(qū)域?qū)Σ糠之惓５膱A滑有一些過度,這與圖5(a)～(c)中對1號異常圖產(chǎn)生的圓滑效果類似。去噪結(jié)果對于定性認(rèn)識夠用,但隨著三維反演算法的發(fā)展,對原始數(shù)據(jù)處理的要求是在去噪的同時盡可能保留“原始信息”,顯然巴特沃斯低通濾波器在濾除噪聲的同時損失了過多的細(xì)節(jié)信息。

圖7 某地實(shí)測磁異常去噪結(jié)果對比Fig.7 Comparison of denoising results of measured magnetic anomalies in a certain area

圖7(c)為采用全局硬閾值二維小波去噪結(jié)果。與圖7(b)相比,小波濾波結(jié)果更接近原始圖像,對數(shù)據(jù)的圓滑效果不顯著,數(shù)據(jù)范圍-584～1 439 nT,除了對負(fù)值的“保護(hù)”有了顯著提升外,對顯著高值也有了較好的保幅效果。當(dāng)然,對正負(fù)異常幅值的“保護(hù)”并不能直接說明去噪效果的優(yōu)劣,但從反演計算對數(shù)據(jù)的要求角度看,圖7(c)顯然保留了更多的細(xì)節(jié)信息。同時存在一些問題,在部分區(qū)域特別是箭頭指向區(qū)域,與原始數(shù)據(jù)的差異較為明顯,對個別小范圍異常出現(xiàn)了濾除,可能是閾值選取的原因。

圖7(d)為使用本文方法去噪結(jié)果。使用DB3小波進(jìn)行處理,進(jìn)行二維小波變換前x、y方向數(shù)據(jù)需延拓到128個。本文中采用點(diǎn)對稱邊界延拓方法,邊界處理的目的是保留小波的重組性,盡量減少邊塊效應(yīng)。前人研究結(jié)果顯示,相較于廣泛應(yīng)用的對稱延拓、周期延拓,點(diǎn)對稱延拓精度高于對稱延拓[23]。3個高頻分量最佳分離階數(shù)為2、2和3階,數(shù)據(jù)范圍-617～1 462 nT。與圖7(c)較為一致,從數(shù)據(jù)范圍看對正負(fù)異常幅值“保護(hù)”更優(yōu),較為顯著的是圖中3個箭頭指向的區(qū)域,在去除高頻噪聲的同時與原始數(shù)據(jù)保持較好的一致性,表明了本文提出的方法對于二維重磁數(shù)據(jù)具有更好的適用性。

從實(shí)際區(qū)域磁數(shù)據(jù)去噪和模型對比試驗(yàn)結(jié)果來看:對于大范圍顯著異常,多數(shù)算法的濾波結(jié)果差異不大,能夠在去除噪聲的同時保持原始異常形態(tài);但對于弱異常的處理則有較大差異,線性濾波通常會出現(xiàn)過度圓滑的情形;二維小波更能適應(yīng)實(shí)際數(shù)據(jù)非線性非平穩(wěn)的特點(diǎn),則能夠在保持異常幅值的同時去除大部分噪聲,這對于反演計算有較為重要的意義。

4 結(jié) 論

(1)二維小波能量閾值去噪算法去除噪聲的相對誤差小于常規(guī)二維小波和線性濾波器,且能夠保持異常形態(tài)不失真,有利于弱異常信號的識別和提取。

(2)含噪信號經(jīng)二維小波分解后,在3個高頻成分中噪聲隨分解階數(shù)的增加能量降低,有效信號則相反。

(3)二維小波能量閾值去噪算法中,小波基函數(shù)和去噪閾值的確定依據(jù)信噪比公式,最佳信噪分離階數(shù)為信號的能量最低階數(shù),3個關(guān)鍵參數(shù)選取有客觀依據(jù),避免了依據(jù)經(jīng)驗(yàn)選擇參數(shù)。