劉洪波

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

0 引言

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

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

1 基本原理

1.1 小子域?yàn)V波

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

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

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

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

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

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

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

1.2 非局部均值濾波

非局部均值濾波[11-14]的基本思想是：計(jì)算以當(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)前像素塊中心的估計(jì)像素值。其具體去噪過程為：

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

(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)歐氏距離計(jì)算

(3)

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

(4)

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

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

2.1 模型試算

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

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

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

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

應(yīng)用非局部均值濾波能在壓制小子域?yàn)V波產(chǎn)生的高頻干擾的同時(shí)保留小子域?yàn)V波對(duì)異常邊界壓縮的良好特性。圖6分別展示了小子域?yàn)V波、小子域?yàn)V波與傳統(tǒng)低通濾波、小子域?yàn)V波與非局部均值濾波相結(jié)合的濾波效果。圖6a表明小子域?yàn)V波具有低通濾波的性能，但小子域?yàn)V波在壓制高頻干擾的同時(shí)使異常產(chǎn)生扭曲，特別是異常邊界部位，使得重力異常等值線不平滑，產(chǎn)生附加高頻干擾，對(duì)重力異常的解釋產(chǎn)生不利影響。圖6b為小子域?yàn)V波基礎(chǔ)上經(jīng)參數(shù)最小的傳統(tǒng)低通濾波后的重力異常，傳統(tǒng)低通濾波可以消減小子域?yàn)V波產(chǎn)生的高頻干擾，但是小子域?yàn)V波良好的異常邊界壓縮特性也被破壞，經(jīng)最小參數(shù)的傳統(tǒng)低通濾波后，幾乎完全模糊了小子域?yàn)V波的邊界壓縮特性，加大傳統(tǒng)低通濾波參數(shù)，會(huì)對(duì)邊界產(chǎn)生更大的模糊作用。因此，傳統(tǒng)低通濾波不適合于消除小子域?yàn)V波的固有缺陷。圖6c為小子域?yàn)V波與非局部均值濾波相結(jié)合的濾波結(jié)果，非局部均值濾波很好地壓制了小子域?yàn)V波產(chǎn)生的高頻干擾，消除了小子域?yàn)V波的邊界扭曲缺陷。同時(shí)，非局部均值濾波很好地保留了小子域?yàn)V波的邊界壓縮特性。由于非局部均值濾波與異常邊界的形態(tài)無關(guān)，小子域?yàn)V波與非局部均值濾波相結(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.小子域?yàn)V波；b.小子域?yàn)V波與傳統(tǒng)低通濾波結(jié)合；c.小子域?yàn)V波與非局部均值濾波結(jié)合.圖6 小子域?yàn)V波及其聯(lián)合濾波結(jié)果Fig.6 Small sub-domain filtering and joint filtering results

2.2 應(yīng)用實(shí)例

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

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

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

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

a.小子域?yàn)V波與傳統(tǒng)低通濾波聯(lián)合；b.小子域?yàn)V波與非局部均值濾波聯(lián)合.圖8 小子域?yàn)V波與傳統(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ù)時(shí)，對(duì)小子域的附加高頻干擾異常壓制比較明顯，但對(duì)小子域?yàn)V波引起異常邊界扭曲改善不明顯(圖9)。當(dāng)采用較強(qiáng)濾波能力的非局部均值濾波參數(shù)時(shí)，小子域?yàn)V波與非局部均值濾波聯(lián)合濾波對(duì)異常邊界產(chǎn)生模糊作用。為了能夠兼顧高頻異常干擾與邊界扭曲改善，即使選取較強(qiáng)的濾波參數(shù)，仍然能很好地保留小子域?yàn)V波的邊界增強(qiáng)作用。

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

3 結(jié)論

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

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

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