非局部均值濾波在重力異常數(shù)據(jù)處理中的應(yīng)用

劉洪波

山東省物化探勘查院，濟(jì)南 250013

0 引言

重力勘探是以地下介質(zhì)密度差異為基礎(chǔ)的一種物探方法[1]，在地質(zhì)構(gòu)造、礦產(chǎn)資源和能源勘探等方面應(yīng)用普遍。重力異常是地下物質(zhì)密度不均勻分布的綜合反映，利用重力異常解決具體地質(zhì)任務(wù)時需對重力異常數(shù)據(jù)進(jìn)行處理，以提取出目標(biāo)地質(zhì)體異常[2]。低通濾波是重力異常數(shù)據(jù)處理的常用方法，能夠消除重力異常中淺層高頻干擾以及非地質(zhì)體引起的高頻干擾，突出異常中的區(qū)域特征。位場邊界特征是位場數(shù)據(jù)解釋非常重要的信息，傳統(tǒng)的低通濾波方法會模糊重力異常的邊界，降低異常分辨率。為了克服低通濾波的這一缺陷，楊高印[3]提出了非線性小子域濾波方法，該方法能以較高的分辨率保留異常的邊界特征，但小子域濾波存在濾波輸出不穩(wěn)定導(dǎo)致的異常曲線扭曲以及附加高頻干擾。馬濤等[4]針對小子域濾波法存在子域劃分重心不穩(wěn)的問題，提出對稱子域劃分方法，使輸出結(jié)果更加穩(wěn)定；張鳳旭等[5]將小子域濾波應(yīng)用于重力異常三方向(X方向、Y方向、XY方向)，提出三方向小子域濾波，提高了小子域濾波法對斷裂平面位置信息的識別；肖鋒等[6]對小子域濾波的小子域剖分方式進(jìn)行了改進(jìn)，提出‘田’型小子域劃分方式，使異常的邊界得到加強(qiáng)；蔣甫玉等[7]在‘田’型子域劃分方式基礎(chǔ)上，加入具有方向信息的四個子域，進(jìn)一步提高小子域濾波的穩(wěn)定性并增強(qiáng)小子域濾波對斷裂水平位置的識別能力；馬國慶等[8，9]提出小子域的‘十’和‘X’型劃分方式，并通過水平導(dǎo)數(shù)改進(jìn)小子域判別準(zhǔn)則，使小子域濾波輸出結(jié)果更穩(wěn)定、合理。

上述對小子域濾波的改進(jìn)主要集中在小子域的劃分方式上，只是針對不同的重力異常數(shù)據(jù)選擇不同的子域劃分方式，能在一定程度上提高小子域濾波的穩(wěn)定性，不可能從根本上改變小子域濾波法導(dǎo)致的重力異常曲線扭曲及高頻干擾引入。通過變換子域劃分方式改進(jìn)小子域濾波的穩(wěn)定性不具有普適性。小子域濾波導(dǎo)致的重力異常曲線扭曲主要表現(xiàn)為高頻成分，可通過低通濾波消除這種高頻干擾，但是傳統(tǒng)低通濾波會模糊異常邊界，這與小子域濾波的初衷相悖。非局部均值濾波[10]是一項新的低通濾波技術(shù)，能充分利用圖像中的冗余信息，在去噪的同時能最大程度地保持圖像的細(xì)節(jié)特征(高頻信息)。因此，筆者提出小子域濾波與非局部均值濾波的聯(lián)合應(yīng)用，既可增強(qiáng)異常邊界又能壓制小子域濾波邊界扭曲效應(yīng)，從而提高小子域濾波的穩(wěn)定性及普適性。

1 基本原理

1.1 小子域濾波

滑動平均可看作是卷積運(yùn)算過程，在濾波過程中卷積因子(卷積核)不變，取滑動窗口內(nèi)不同數(shù)據(jù)的平均值作為窗口中心的輸出結(jié)果，模糊了異常的邊界特征。小子域濾波是對滑動平均法的改進(jìn)，以滑動窗口中心為中心點(diǎn)劃分八個包含中心點(diǎn)的子域(圖1)，并計算每個子域的均方差；以均方差最小的子域平均值作為整個滑動窗口中心的輸出結(jié)果。小子域濾波過程的卷積因子是變化的，是一種非線性濾波。

圖1 ‘米’型小子域劃分方式Fig.1 Division method of “union jack” small sub-domain

以‘米’型子域劃分為例說明小子域濾波造成異常曲線扭曲的原因。為清晰地說明問題，假設(shè)異常的邊界如圖2a所示，將圖2a離散化，邊界兩側(cè)的重力異常值分別為1，2(圖2b)。圖2b中顯示了根據(jù)小子域濾波原理確定的六個子域，可以看出，小子域在異常邊界處很快由低異常區(qū)轉(zhuǎn)到高異常區(qū)，實現(xiàn)了滑動平均過程對異常邊界的保護(hù)。但由低異常區(qū)到高異常區(qū)的快速轉(zhuǎn)變，導(dǎo)致了濾波過程強(qiáng)烈的非線性。

注：圖中(紅色：1，1，2，2，2，2)表示窗口中心，不同顏色實線三角區(qū)域為確定的不同窗口的小子域.a.異常邊界模型；b.離散化后數(shù)值形式.圖2 小子域濾波邊界扭曲產(chǎn)生原因分析Fig.2 Cause analysis of abnormal boundary distortion in small sub-domain filtering

圖3顯示了六個滑動窗口小子域濾波和滑動平均濾波的結(jié)果，滑動平均在邊界兩側(cè)的濾波輸出近似線性，因此濾波結(jié)果比較穩(wěn)定，而小子域濾波在邊界兩側(cè)表現(xiàn)為強(qiáng)烈的非線性，使異常邊界彎折，非常小的噪聲干擾能引起非常大的濾波誤差，導(dǎo)致濾波結(jié)果不穩(wěn)定，造成異常邊界的扭曲。特別是當(dāng)異常邊界比較復(fù)雜時，非線性放大作用更為明顯，不僅造成異常邊界扭曲，而且會引入附加高頻干擾。

圖3 六個窗口中心不同濾波結(jié)果Fig.3 Different filtering results in six window-center

不同小子域劃分方式及子域確定準(zhǔn)則能一定程度上減弱小子域濾波的這種非線性，但是不能從根本上解決小子域濾波輸出的不穩(wěn)定性。雖然小子域濾波輸出不穩(wěn)定是固有缺陷，但是這種不穩(wěn)定性常表現(xiàn)高頻性質(zhì)。因此，可通過一種既能實現(xiàn)低通濾波又能很好保護(hù)高頻信息的濾波方法消減小子域濾波造成的異常曲線扭曲。非局部均值濾波能很好地實現(xiàn)這一目的。

1.2 非局部均值濾波

非局部均值濾波[11-14]的基本思想是：計算以當(dāng)前像素為中心的像素塊與鄰域內(nèi)其他像素塊之間的高斯加權(quán)歐氏距離；以加權(quán)歐氏距離與濾波系數(shù)為變量構(gòu)建高斯函數(shù)作為圖像塊間的權(quán)值函數(shù)，構(gòu)建此函數(shù)的目的是使與當(dāng)前像素塊相同或者相似的像素塊獲得較大的權(quán)重，不相似的像素塊則得到較小的權(quán)重，通過加權(quán)平均這些權(quán)值得到當(dāng)前像素塊中心的估計像素值。其具體去噪過程為：

假設(shè)含噪圖像v={v(i)|i∈I}，v(i)表示第i個像素的像素值，其估計值vNL(i)可通過圖像內(nèi)所有像素值的加權(quán)平均計算

(1)

式(1)中，w(i,j)為與像素i和像素j相似性相關(guān)的加權(quán)系數(shù)。像素i和像素j的相似性可通過以像素i和j為中心的圖像塊的高斯加權(quán)歐氏距離衡量，即

(2)

式(2)中，Ga表示標(biāo)準(zhǔn)差a的高斯核，Nk表示以像素k為中心的固定大小的圖像塊。加權(quán)系數(shù)w(i,j)可由高斯加權(quán)歐氏距離計算

(3)

式(3)中，h為濾波參數(shù)，Z(i)為歸一化常數(shù),其中：

(4)

非局部均值濾波充分考慮了像素i和其他所有像素j的關(guān)聯(lián)性，能最大程度地保留圖像的高頻細(xì)節(jié)特征。將非局部均值濾波與小子域濾波相結(jié)合既能改善小子域濾波的邊界扭曲缺陷又能最大程度地保留小子域濾波優(yōu)良的邊界特征，提高小子域濾波的穩(wěn)定性與普適性。

2 重力異常數(shù)據(jù)處理

2.1 模型試算

為了驗證非局部均值濾波與小子域濾波相結(jié)合的方法可行性和有效性，設(shè)計如圖4所示的地質(zhì)模型，包含3個長方地質(zhì)體，其具體參數(shù)如表1所示。

表1 長方體正演參數(shù)

圖4 模擬的地質(zhì)模型Fig.4 Simulated geological model

根據(jù)表1中的各項參數(shù)及圖4所示地質(zhì)模型位置，通過正演模擬獲得的地質(zhì)模型的重力異常如圖5a所示。圖5b為模擬的重力異常加入小功率的高斯白噪聲以仿真實際采集時的高頻干擾，同時檢驗算法穩(wěn)定性[15]。對模擬的重力異常分別采用傳統(tǒng)低通濾波和非局部均值濾波進(jìn)行處理，處理結(jié)果如圖5c，5d所示?？梢钥闯?，傳統(tǒng)低通濾波在壓制高頻干擾時使異常邊界模糊(變寬)，而非局部均值濾波在壓制干擾的同時能很好地保留異常邊界的細(xì)節(jié)信息，不會使異常邊界過渡平滑。

應(yīng)用非局部均值濾波能在壓制小子域濾波產(chǎn)生的高頻干擾的同時保留小子域濾波對異常邊界壓縮的良好特性。圖6分別展示了小子域濾波、小子域濾波與傳統(tǒng)低通濾波、小子域濾波與非局部均值濾波相結(jié)合的濾波效果。圖6a表明小子域濾波具有低通濾波的性能，但小子域濾波在壓制高頻干擾的同時使異常產(chǎn)生扭曲，特別是異常邊界部位，使得重力異常等值線不平滑，產(chǎn)生附加高頻干擾，對重力異常的解釋產(chǎn)生不利影響。圖6b為小子域濾波基礎(chǔ)上經(jīng)參數(shù)最小的傳統(tǒng)低通濾波后的重力異常，傳統(tǒng)低通濾波可以消減小子域濾波產(chǎn)生的高頻干擾，但是小子域濾波良好的異常邊界壓縮特性也被破壞，經(jīng)最小參數(shù)的傳統(tǒng)低通濾波后，幾乎完全模糊了小子域濾波的邊界壓縮特性，加大傳統(tǒng)低通濾波參數(shù)，會對邊界產(chǎn)生更大的模糊作用。因此，傳統(tǒng)低通濾波不適合于消除小子域濾波的固有缺陷。圖6c為小子域濾波與非局部均值濾波相結(jié)合的濾波結(jié)果，非局部均值濾波很好地壓制了小子域濾波產(chǎn)生的高頻干擾，消除了小子域濾波的邊界扭曲缺陷。同時，非局部均值濾波很好地保留了小子域濾波的邊界壓縮特性。由于非局部均值濾波與異常邊界的形態(tài)無關(guān)，小子域濾波與非局部均值濾波相結(jié)合的方法不需要考慮異常邊界的形態(tài)，具有更強(qiáng)的普適性與穩(wěn)定性。

a.地質(zhì)模型產(chǎn)生的重力異常；b.加入小功率高斯白噪聲后的重力異常；c.傳統(tǒng)低通濾波；d.非局部均值濾波.圖5 模擬重力異常及不同濾波方法處理結(jié)果Fig.5 Simulated gravity anomaly and different filtering processing results

a.小子域濾波；b.小子域濾波與傳統(tǒng)低通濾波結(jié)合；c.小子域濾波與非局部均值濾波結(jié)合.圖6 小子域濾波及其聯(lián)合濾波結(jié)果Fig.6 Small sub-domain filtering and joint filtering results

2.2 應(yīng)用實例

通過萊州灣海域?qū)崪y重力異常數(shù)據(jù)處理，分析不同濾波方法的濾波效果，說明小子域濾波與非局部均值濾波相結(jié)合的濾波方法對實際數(shù)據(jù)處理的有效性和優(yōu)越性。

圖7顯示了實測重力異常及3種低通濾波結(jié)果。小子域濾波增強(qiáng)了異常邊界形態(tài)，使異常邊界更加清晰，但其使異常邊界發(fā)生扭曲，引入嚴(yán)重的其他高頻干擾，降低了重力異常數(shù)據(jù)的質(zhì)量(圖7b)；圖7c為滑動窗5×5傳統(tǒng)滑動平均濾波，可以看出傳統(tǒng)低通濾波能很好地濾除高頻干擾，對異常平滑作用明顯，但邊界形態(tài)模糊嚴(yán)重(圖7c①、②、③、④)，與原始異常數(shù)據(jù)相比，異常邊界等值線變得稀疏，異常分辨率降低；圖7d為非局部均值濾波結(jié)果，非局部均值濾波能實現(xiàn)傳統(tǒng)低通濾波效果，同時能夠保護(hù)異常邊界形態(tài)(圖7d①、②、③、④)，與原始異常數(shù)據(jù)相比，幾乎不改變異常邊界形態(tài)，能最大程度保護(hù)異常邊界。

a.原始數(shù)據(jù)；b.小子域濾波；c.傳統(tǒng)低通濾波；d.非局部均值濾波.圖7 原始重力異常及不同低通濾波結(jié)果Fig.7 Original gravity anomaly and different low-pass filtering results

通過圖7分析可知，小子域濾波能夠增強(qiáng)異常邊界使異常保持很高的分辨率，但小子域濾波會產(chǎn)生嚴(yán)重的異常失真及邊界扭曲，而非局部均值濾波能在實現(xiàn)低通濾波的同時很好地保護(hù)異常邊界信息。圖8a、b分別展示了小子域濾波與傳統(tǒng)低通濾波和非局部均值濾波相結(jié)合的濾波效果，其中傳統(tǒng)低通濾波的參數(shù)為3×3滑動窗口。圖8a中可以看出，即使傳統(tǒng)低通濾波采用最小窗口，其對異常邊界也有很大的模糊作用，幾乎沒有保留小子域濾波的良好邊界保護(hù)特性，但其對邊界的保護(hù)比直接進(jìn)行傳統(tǒng)濾波要好；圖8b可看出，非局部均值濾波與小子域濾波相結(jié)合的濾波方法，不但能消除小子域濾波的異常失真，而且能最大程度保留小子域濾波對異常邊界增強(qiáng)特性，使得異常的邊界相對于直接非局部均值濾波更加清晰。

a.小子域濾波與傳統(tǒng)低通濾波聯(lián)合；b.小子域濾波與非局部均值濾波聯(lián)合.圖8 小子域濾波與傳統(tǒng)低通濾波、非局部均值濾波聯(lián)合濾波Fig.8 Combined filtering with small domain filtering and traditional low-pass filtering and non-local mean filtering

當(dāng)采用較弱非局部濾波參數(shù)時，對小子域的附加高頻干擾異常壓制比較明顯，但對小子域濾波引起異常邊界扭曲改善不明顯(圖9)。當(dāng)采用較強(qiáng)濾波能力的非局部均值濾波參數(shù)時，小子域濾波與非局部均值濾波聯(lián)合濾波對異常邊界產(chǎn)生模糊作用。為了能夠兼顧高頻異常干擾與邊界扭曲改善，即使選取較強(qiáng)的濾波參數(shù)，仍然能很好地保留小子域濾波的邊界增強(qiáng)作用。

圖9 小子域濾波與弱參數(shù)非局部均值濾波聯(lián)合Fig.9 Combination of small domain filtering and weak parameter nonlocal mean filtering

3 結(jié)論

(1)小子域濾波產(chǎn)生的異常邊界扭曲難以通過子域劃分方式的改進(jìn)加以克服，可通過非局部均值濾波加以改善。

(2)小子域濾波與非局部均值濾波聯(lián)合，既增強(qiáng)了異常邊界又壓制了小子域濾波邊界扭曲效應(yīng)，從而提高了小子域濾波的穩(wěn)定性。

(3)小子域濾波與非局部均值濾波聯(lián)合濾波不需要考慮異常局部形態(tài)與子域形態(tài)的關(guān)系，因而具有更強(qiáng)的普適性，濾波輸出更合理，應(yīng)用效果更好。