磁法數據處理中的擴邊和優化中值濾波方法的研究

馬國慶,孟令順,杜曉娟,張鳳旭

(吉林大學地球探測科學與技術學院,吉林長春 130026)

0 前言

在重、磁數據的處理中,往往在處理前,首先要對原始數據進行擴邊處理。常用的擴邊方法有余弦衰減至零和對折擴邊方法,這二種方法都存在明顯的邊界效應[1],而邊界效應不利于對異常進行解釋。為了減小邊界效應,作者在本文提出了三方向擴邊方法。

磁測數據受外界因素影響較大,即使一些小的磁性干擾,都會給實際數據的處理和解釋帶來困難,所以必須利用相應的處理方法進行消除。中值濾波具有消除突跳點和保護邊緣的特征[2、3],作者在本文對中值濾波方法進行了改進,提出了優化中值濾波方法。實踐證明,利用該方法對磁測數據進行處理,能夠取得較好的效果。

1 三方向擴邊方法

1.1 三方向擴邊方法的原理

實際測量中,某一測點的數據不僅受到目標地質體的影響,還受到來自周圍一定范圍內的其它地質體的影響,因此該點的數據是一個綜合反映[4]。三方向擴邊方法的原理是對一個點的三方向同時進行擴邊,取擴邊值的平均值作為最后的擴邊結果。這樣做可綜合周圍地質體的影響,減小邊界效應。

1.2 三方向擴邊方法計算步驟

進行縱向擴邊,對于D點(見下頁圖1),首先分別計算出Y2,6→D,Y2,9→D,Y2,12→D三個方向上相鄰點的異常差平均值;然后根據邊界值和異常差,分別在三個方向上計算出D點的值;最后取三個方向上D值的平均值,做為D點的擴邊值。B、C二點因為在擴邊數據的邊角處,因此采用特殊的擴邊方法。對于B點,由于只在三個方向當中的二個方向上存在原始數據,所以只利用這二個方向的邊界值和差值來計算出最后結果。對于C點,將采用“二長一短”的方法來進行計算,因為在三個方向當中有二個方向上的數據量充足,而另一個方向的數據量不足,最后的結果也是三個方向的結果的平均。在橫向也采用上述方法進行擴邊。對于A點,因為是交點,所以在橫向、縱向擴邊計算完畢以后再對其進行擴邊。

圖1 三方向擴邊方法Fig.1 Three directions of expansion methods

三方向擴邊方法的計算公式為:

當i=0時為D點的計算公式:

其余點的計算公式可以依據公式(3)來獲得。

1.3 模型數據驗算

理論模型:一個平坦的地區存在二個位置不同,深度不同的垂直磁化的球體,所產生的疊加磁異常見圖2。

圖2 理論模型計算結果Fig.2 The results of theoreticalmodel calculations

采用對折擴邊方法和三方向擴邊方法,分別對理論數據進行擴邊,如圖3所示。

模型試算表明,作者在本文提出的三方向擴邊法,基本上不會引起邊界效應,并且能使邊界信息得到很好的保留。

圖3 不同擴邊方法的結果對比Fig.3 The comparison of differentmethods of edge expansion

2 優化中值濾波在磁法數據處理中的應用

在實際測量中,當外界干擾較大,對異常的處理和解釋產生較大的影響時,必須想辦法將干擾消除掉。

中值濾波是常用的一種非線性平滑濾波。它是一種領域運算,類似于卷積,但不是加權求和計算,而是把領域中的像素按灰度等級進行排序,然后選擇該組的中間值作為輸出像素值。它能減弱或消除傅里葉空間的高頻分量,并且影響低頻分量。因為高頻分量對應圖像中的區域邊緣的灰度值,具有較大較快變化的部份,所以該濾波可將這些分量濾除,使圖像平滑[5]。

作者在本文以中值濾波為基礎,結合干擾的特征以及異常的屬性,提出了優化中值濾波的概念。面積性優化中值濾波是以環帶的形式來進行的(見圖4),環帶的多少由干擾異常的范圍來決定。

圖4中的矩形區域為干擾區,數字“1”代表的是第一個環帶(包括邊界的數據)內的數據;數字“2”環指的是第二個環帶(包括邊界數據)與第一個環帶(不包括邊界數據)之間的數據;數字“3”環及數字“4”環的數據獲得方式,與數字“2”環是一致的。

中值濾波的點數不同,處理后的結果反映的信息也不同。利用這種形式來進行中值濾波的優點是:

(1)由于有三個點以上的異常才稱為異常,因此選擇一環的目的是去除突跳點,在無干擾的區域不會使異常幅值變小。

(2)二環、三環是一個過渡區域,因為在某些區域干擾的范圍可能會小一些,所以加這二個環帶的目的,是為了去除范圍小的干擾。

圖4 優化中值濾波示意圖Fig.4 Schematic diagram of filter optimization

(3)四環是主要的濾波區域,在這個區域內無干擾的數據多于干擾數據,對這個環帶內的數據進行中值濾波,能夠較好地去除干擾。

對于剖面濾波是以不同大小的窗口對曲線進行濾波。

對數據進行濾波之后,根據其異常和干擾的特點,加上相應的限制條件會使處理結果更可靠。

3 實際數據驗算

以四川省綿陽市射洪縣地區的一條日變曲線為例,來驗證優化中值濾波方法的實用性。

從圖5中可以看出,在“1”處、“2”處由于受到人為的影響而出現蹦跳。在理論上,一天的日變曲線在沒有大小磁暴的干擾時,應為一條比較平滑的曲線,因此需要將日變干擾造成的蹦跳改正掉。圖6(見下頁)是在不同窗口下,對實測日變曲線進行中值濾波后的結果。

圖5 四川省綿陽市射洪縣地區日變曲線(2009年4月27日)Fig.5 Diurnal variation curve on Shehong county area,Mianyang,Sichuan province

從圖6中可以看出,不同窗口的中值濾波,對“2”號干擾都起到了作用。但是在“1”號干擾處,三個窗口和五個窗口的中值濾波效果明顯不如七個窗口的中值濾波,但是七個窗口的中值濾波在極值處存在一定壓制。

見圖7(見下頁),結合不同窗口對數據的處理效果及干擾為負值的特點,將以上三個窗口濾波后的數據進行比較,選取最大值來作為最終的結果。

從圖7中可以看出,優化中值濾波可以很好地去掉磁測中的偶然誤差,這種方法可以應用到實測數據處理中。

圖8(見下頁)是內蒙古呼圖格地區的實測平面磁法數據,以此為例來驗證擴邊方法和優化中值濾波的聯合應用效果。

圖8中左側的負異常條帶,是由高壓線引起的。為了對資料進行準確的地質解釋,必須將這種干擾去除掉。利用優化中值濾波對其進行處理的步驟[6]如下頁圖9所示。

由原始圖像可知,干擾異常為負值且異常寬度為1點～4個點,一環、二環、三環帶的作用主要是去掉突跳點和去除范圍小的干擾,四環帶對磁干擾條帶有很好的去除作用,但是對極值有壓制。由于本地區進行磁法測量的主要目的是尋找鐵磁性物質,因此在中值濾波之后,取四個值的最大值做為最終的計算結果(見下頁圖10)。

圖11(見下頁)是采用與優化中值濾波同樣大小窗口的方域中值濾波。

從圖10、圖11二幅圖的對比中可以看出,利用三方向擴邊方法,在邊界區域未出現邊界效應,所以這種擴邊方法實用性較好。

從圖10與圖11的對比中可以看出,優化中值濾波的效果比普通中值濾波的效果較好,在將干擾很好的去除的前提下,并沒有使異常的極值降低很多。從原始圖像中我們可以發現,在負值區域內,坐標為(16,42)的位置存在一個小的正異常,普通的方域中值濾波將這個異常給模糊掉了,而優化中值濾波將其完好的保留。在優化中值濾波使異常的處理與解釋更可靠,因此這種方法實用性較高。

圖6 不同窗口所做的中值濾波結果對比Fig.6 Comparative results of median filter in different windows

圖7 經過濾波后的日變曲線與實測結果的對比Fig.7 Comparison of diurnal variation curve and measured results after filter

圖8 內蒙古呼圖格地區的實測磁異常平面圖Fig.8 The measured magnetic map of Hutuge area,Neimenggu

圖9 優化中值濾波的計算步驟Fig.9 Calculation steps of optimization median filter

4 結論

從實際資料的處理中可以看出,作者在本文提出的三方向擴邊方法不會引起邊界效應,并且還可以很好地保留邊界信息。從圖像的對比中我們會看出,優化中值濾波比普通中值濾波效果好,不僅能將干擾去除,而且還不會降低幅值,對信息有很好的保留效果。因此以上方法可以在實際數據的處理中得以廣泛應用。

圖10 優化中值濾波后的結果Fig.10 Results of optimization Median filter

圖11 方域中值濾波后的結果Fig.11 Results of domain median filter

[1] 段本春,徐世浙.磁(重力)異常局部場與區域場分離處理中的擴邊方法研究[J].物探化探計算技術,1997,19(4):299.

[2] 劉財,王典,劉洋,等.二維多級中值濾波技術在隨機噪聲消除中的應用初探[J].石油地球物理勘探,2005,40(2):163.

[3] 邵 敏,何展翔,張立恩.二維中值濾波在電磁測深曲線去噪中的應用[J].物探化探計算技術,2001,23(3):226.

[4] 劉東甲,程方道.劃分重力區域場與局部場的多次切割法[J].物探化探計算技術,1997,19(1):32.

[5] SPECTOR A.Comments and errata to classic papers[J].Geophysics,1985,50(11):2278.

[6] ARCE G R,FOSTER R E.Detail-preserving rankedor-der based filter for image processing[J].IEEE Transon ASSP,1987,37(1):83.