馬國(guó)慶,孟令順,杜曉娟,張鳳旭
(吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院,吉林長(zhǎng)春 130026)
在重、磁數(shù)據(jù)的處理中,往往在處理前,首先要對(duì)原始數(shù)據(jù)進(jìn)行擴(kuò)邊處理。常用的擴(kuò)邊方法有余弦衰減至零和對(duì)折擴(kuò)邊方法,這二種方法都存在明顯的邊界效應(yīng)[1],而邊界效應(yīng)不利于對(duì)異常進(jìn)行解釋。為了減小邊界效應(yīng),作者在本文提出了三方向擴(kuò)邊方法。
磁測(cè)數(shù)據(jù)受外界因素影響較大,即使一些小的磁性干擾,都會(huì)給實(shí)際數(shù)據(jù)的處理和解釋帶來(lái)困難,所以必須利用相應(yīng)的處理方法進(jìn)行消除。中值濾波具有消除突跳點(diǎn)和保護(hù)邊緣的特征[2、3],作者在本文對(duì)中值濾波方法進(jìn)行了改進(jìn),提出了優(yōu)化中值濾波方法。實(shí)踐證明,利用該方法對(duì)磁測(cè)數(shù)據(jù)進(jìn)行處理,能夠取得較好的效果。
實(shí)際測(cè)量中,某一測(cè)點(diǎn)的數(shù)據(jù)不僅受到目標(biāo)地質(zhì)體的影響,還受到來(lái)自周?chē)欢ǚ秶鷥?nèi)的其它地質(zhì)體的影響,因此該點(diǎn)的數(shù)據(jù)是一個(gè)綜合反映[4]。三方向擴(kuò)邊方法的原理是對(duì)一個(gè)點(diǎn)的三方向同時(shí)進(jìn)行擴(kuò)邊,取擴(kuò)邊值的平均值作為最后的擴(kuò)邊結(jié)果。這樣做可綜合周?chē)刭|(zhì)體的影響,減小邊界效應(yīng)。
進(jìn)行縱向擴(kuò)邊,對(duì)于D點(diǎn)(見(jiàn)下頁(yè)圖1),首先分別計(jì)算出Y2,6→D,Y2,9→D,Y2,12→D三個(gè)方向上相鄰點(diǎn)的異常差平均值;然后根據(jù)邊界值和異常差,分別在三個(gè)方向上計(jì)算出D點(diǎn)的值;最后取三個(gè)方向上D值的平均值,做為D點(diǎn)的擴(kuò)邊值。B、C二點(diǎn)因?yàn)樵跀U(kuò)邊數(shù)據(jù)的邊角處,因此采用特殊的擴(kuò)邊方法。對(duì)于B點(diǎn),由于只在三個(gè)方向當(dāng)中的二個(gè)方向上存在原始數(shù)據(jù),所以只利用這二個(gè)方向的邊界值和差值來(lái)計(jì)算出最后結(jié)果。對(duì)于C點(diǎn),將采用“二長(zhǎng)一短”的方法來(lái)進(jìn)行計(jì)算,因?yàn)樵谌齻€(gè)方向當(dāng)中有二個(gè)方向上的數(shù)據(jù)量充足,而另一個(gè)方向的數(shù)據(jù)量不足,最后的結(jié)果也是三個(gè)方向的結(jié)果的平均。在橫向也采用上述方法進(jìn)行擴(kuò)邊。對(duì)于A點(diǎn),因?yàn)槭墙稽c(diǎn),所以在橫向、縱向擴(kuò)邊計(jì)算完畢以后再對(duì)其進(jìn)行擴(kuò)邊。
圖1 三方向擴(kuò)邊方法Fig.1 Three directions of expansion methods
三方向擴(kuò)邊方法的計(jì)算公式為:
當(dāng)i=0時(shí)為D點(diǎn)的計(jì)算公式:
其余點(diǎn)的計(jì)算公式可以依據(jù)公式(3)來(lái)獲得。
理論模型:一個(gè)平坦的地區(qū)存在二個(gè)位置不同,深度不同的垂直磁化的球體,所產(chǎn)生的疊加磁異常見(jiàn)圖2。
圖2 理論模型計(jì)算結(jié)果Fig.2 The results of theoreticalmodel calculations
采用對(duì)折擴(kuò)邊方法和三方向擴(kuò)邊方法,分別對(duì)理論數(shù)據(jù)進(jìn)行擴(kuò)邊,如圖3所示。
模型試算表明,作者在本文提出的三方向擴(kuò)邊法,基本上不會(huì)引起邊界效應(yīng),并且能使邊界信息得到很好的保留。
圖3 不同擴(kuò)邊方法的結(jié)果對(duì)比Fig.3 The comparison of differentmethods of edge expansion
在實(shí)際測(cè)量中,當(dāng)外界干擾較大,對(duì)異常的處理和解釋產(chǎn)生較大的影響時(shí),必須想辦法將干擾消除掉。
中值濾波是常用的一種非線性平滑濾波。它是一種領(lǐng)域運(yùn)算,類(lèi)似于卷積,但不是加權(quán)求和計(jì)算,而是把領(lǐng)域中的像素按灰度等級(jí)進(jìn)行排序,然后選擇該組的中間值作為輸出像素值。它能減弱或消除傅里葉空間的高頻分量,并且影響低頻分量。因?yàn)楦哳l分量對(duì)應(yīng)圖像中的區(qū)域邊緣的灰度值,具有較大較快變化的部份,所以該濾波可將這些分量濾除,使圖像平滑[5]。
作者在本文以中值濾波為基礎(chǔ),結(jié)合干擾的特征以及異常的屬性,提出了優(yōu)化中值濾波的概念。面積性?xún)?yōu)化中值濾波是以環(huán)帶的形式來(lái)進(jìn)行的(見(jiàn)圖4),環(huán)帶的多少由干擾異常的范圍來(lái)決定。
圖4中的矩形區(qū)域?yàn)楦蓴_區(qū),數(shù)字“1”代表的是第一個(gè)環(huán)帶(包括邊界的數(shù)據(jù))內(nèi)的數(shù)據(jù);數(shù)字“2”環(huán)指的是第二個(gè)環(huán)帶(包括邊界數(shù)據(jù))與第一個(gè)環(huán)帶(不包括邊界數(shù)據(jù))之間的數(shù)據(jù);數(shù)字“3”環(huán)及數(shù)字“4”環(huán)的數(shù)據(jù)獲得方式,與數(shù)字“2”環(huán)是一致的。
中值濾波的點(diǎn)數(shù)不同,處理后的結(jié)果反映的信息也不同。利用這種形式來(lái)進(jìn)行中值濾波的優(yōu)點(diǎn)是:
(1)由于有三個(gè)點(diǎn)以上的異常才稱(chēng)為異常,因此選擇一環(huán)的目的是去除突跳點(diǎn),在無(wú)干擾的區(qū)域不會(huì)使異常幅值變小。
(2)二環(huán)、三環(huán)是一個(gè)過(guò)渡區(qū)域,因?yàn)樵谀承﹨^(qū)域干擾的范圍可能會(huì)小一些,所以加這二個(gè)環(huán)帶的目的,是為了去除范圍小的干擾。
圖4 優(yōu)化中值濾波示意圖Fig.4 Schematic diagram of filter optimization
(3)四環(huán)是主要的濾波區(qū)域,在這個(gè)區(qū)域內(nèi)無(wú)干擾的數(shù)據(jù)多于干擾數(shù)據(jù),對(duì)這個(gè)環(huán)帶內(nèi)的數(shù)據(jù)進(jìn)行中值濾波,能夠較好地去除干擾。
對(duì)于剖面濾波是以不同大小的窗口對(duì)曲線進(jìn)行濾波。
對(duì)數(shù)據(jù)進(jìn)行濾波之后,根據(jù)其異常和干擾的特點(diǎn),加上相應(yīng)的限制條件會(huì)使處理結(jié)果更可靠。
以四川省綿陽(yáng)市射洪縣地區(qū)的一條日變曲線為例,來(lái)驗(yàn)證優(yōu)化中值濾波方法的實(shí)用性。
從圖5中可以看出,在“1”處、“2”處由于受到人為的影響而出現(xiàn)蹦跳。在理論上,一天的日變曲線在沒(méi)有大小磁暴的干擾時(shí),應(yīng)為一條比較平滑的曲線,因此需要將日變干擾造成的蹦跳改正掉。圖6(見(jiàn)下頁(yè))是在不同窗口下,對(duì)實(shí)測(cè)日變曲線進(jìn)行中值濾波后的結(jié)果。
圖5 四川省綿陽(yáng)市射洪縣地區(qū)日變曲線(2009年4月27日)Fig.5 Diurnal variation curve on Shehong county area,Mianyang,Sichuan province
從圖6中可以看出,不同窗口的中值濾波,對(duì)“2”號(hào)干擾都起到了作用。但是在“1”號(hào)干擾處,三個(gè)窗口和五個(gè)窗口的中值濾波效果明顯不如七個(gè)窗口的中值濾波,但是七個(gè)窗口的中值濾波在極值處存在一定壓制。
見(jiàn)圖7(見(jiàn)下頁(yè)),結(jié)合不同窗口對(duì)數(shù)據(jù)的處理效果及干擾為負(fù)值的特點(diǎn),將以上三個(gè)窗口濾波后的數(shù)據(jù)進(jìn)行比較,選取最大值來(lái)作為最終的結(jié)果。
從圖7中可以看出,優(yōu)化中值濾波可以很好地去掉磁測(cè)中的偶然誤差,這種方法可以應(yīng)用到實(shí)測(cè)數(shù)據(jù)處理中。
圖8(見(jiàn)下頁(yè))是內(nèi)蒙古呼圖格地區(qū)的實(shí)測(cè)平面磁法數(shù)據(jù),以此為例來(lái)驗(yàn)證擴(kuò)邊方法和優(yōu)化中值濾波的聯(lián)合應(yīng)用效果。
圖8中左側(cè)的負(fù)異常條帶,是由高壓線引起的。為了對(duì)資料進(jìn)行準(zhǔn)確的地質(zhì)解釋,必須將這種干擾去除掉。利用優(yōu)化中值濾波對(duì)其進(jìn)行處理的步驟[6]如下頁(yè)圖9所示。
由原始圖像可知,干擾異常為負(fù)值且異常寬度為1點(diǎn)~4個(gè)點(diǎn),一環(huán)、二環(huán)、三環(huán)帶的作用主要是去掉突跳點(diǎn)和去除范圍小的干擾,四環(huán)帶對(duì)磁干擾條帶有很好的去除作用,但是對(duì)極值有壓制。由于本地區(qū)進(jìn)行磁法測(cè)量的主要目的是尋找鐵磁性物質(zhì),因此在中值濾波之后,取四個(gè)值的最大值做為最終的計(jì)算結(jié)果(見(jiàn)下頁(yè)圖10)。
圖11(見(jiàn)下頁(yè))是采用與優(yōu)化中值濾波同樣大小窗口的方域中值濾波。
從圖10、圖11二幅圖的對(duì)比中可以看出,利用三方向擴(kuò)邊方法,在邊界區(qū)域未出現(xiàn)邊界效應(yīng),所以這種擴(kuò)邊方法實(shí)用性較好。
從圖10與圖11的對(duì)比中可以看出,優(yōu)化中值濾波的效果比普通中值濾波的效果較好,在將干擾很好的去除的前提下,并沒(méi)有使異常的極值降低很多。從原始圖像中我們可以發(fā)現(xiàn),在負(fù)值區(qū)域內(nèi),坐標(biāo)為(16,42)的位置存在一個(gè)小的正異常,普通的方域中值濾波將這個(gè)異常給模糊掉了,而優(yōu)化中值濾波將其完好的保留。在優(yōu)化中值濾波使異常的處理與解釋更可靠,因此這種方法實(shí)用性較高。
圖6 不同窗口所做的中值濾波結(jié)果對(duì)比Fig.6 Comparative results of median filter in different windows
圖7 經(jīng)過(guò)濾波后的日變曲線與實(shí)測(cè)結(jié)果的對(duì)比Fig.7 Comparison of diurnal variation curve and measured results after filter
圖8 內(nèi)蒙古呼圖格地區(qū)的實(shí)測(cè)磁異常平面圖Fig.8 The measured magnetic map of Hutuge area,Neimenggu
圖9 優(yōu)化中值濾波的計(jì)算步驟Fig.9 Calculation steps of optimization median filter
從實(shí)際資料的處理中可以看出,作者在本文提出的三方向擴(kuò)邊方法不會(huì)引起邊界效應(yīng),并且還可以很好地保留邊界信息。從圖像的對(duì)比中我們會(huì)看出,優(yōu)化中值濾波比普通中值濾波效果好,不僅能將干擾去除,而且還不會(huì)降低幅值,對(duì)信息有很好的保留效果。因此以上方法可以在實(shí)際數(shù)據(jù)的處理中得以廣泛應(yīng)用。
圖10 優(yōu)化中值濾波后的結(jié)果Fig.10 Results of optimization Median filter
圖11 方域中值濾波后的結(jié)果Fig.11 Results of domain median filter
[1] 段本春,徐世浙.磁(重力)異常局部場(chǎng)與區(qū)域場(chǎng)分離處理中的擴(kuò)邊方法研究[J].物探化探計(jì)算技術(shù),1997,19(4):299.
[2] 劉財(cái),王典,劉洋,等.二維多級(jí)中值濾波技術(shù)在隨機(jī)噪聲消除中的應(yīng)用初探[J].石油地球物理勘探,2005,40(2):163.
[3] 邵 敏,何展翔,張立恩.二維中值濾波在電磁測(cè)深曲線去噪中的應(yīng)用[J].物探化探計(jì)算技術(shù),2001,23(3):226.
[4] 劉東甲,程方道.劃分重力區(qū)域場(chǎng)與局部場(chǎng)的多次切割法[J].物探化探計(jì)算技術(shù),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.