陳 雄，曾昭發(fā)，孫東明，李 靜

(1.吉林大學(xué)，吉林長(zhǎng)春 130026;2.吉林省公路勘測(cè)設(shè)計(jì)院，吉林長(zhǎng)春 130026)

用于航空電磁數(shù)據(jù)調(diào)平處理的導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)法

陳 雄1，曾昭發(fā)1，孫東明2，李 靜1

(1.吉林大學(xué)，吉林長(zhǎng)春 130026;2.吉林省公路勘測(cè)設(shè)計(jì)院，吉林長(zhǎng)春 130026)

航空電磁數(shù)據(jù)中經(jīng)常存在條帶狀誤差，為消除這種誤差，通常采用窗口濾波。但是，窗口參數(shù)及濾波參數(shù)要求極為嚴(yán)格，且不能區(qū)分調(diào)平誤差和其類似大小的地質(zhì)異常，需要操作者多次嘗試以尋求合適的處理參數(shù)。這里提出了導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)法，該方法首先在沿偽切割線方向(垂直飛行方向)，對(duì)網(wǎng)格化的航空電磁數(shù)據(jù)求取偽切割線的導(dǎo)數(shù)，采用高通濾波方法濾掉幅值比較小的導(dǎo)數(shù)突變點(diǎn);然后統(tǒng)計(jì)每條航線上的突變點(diǎn)個(gè)數(shù)，根據(jù)預(yù)設(shè)的閾值，確定出存在條帶誤差的航線;再針對(duì)這條航線上的數(shù)據(jù)，采用均值濾波得到調(diào)平結(jié)果。導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)法能自動(dòng)找出存在誤差的航線，并只針對(duì)誤差線進(jìn)行濾波，能有效消除條帶狀誤差，極大地保留真實(shí)地質(zhì)異常。導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)方法算法簡(jiǎn)單，輸入?yún)?shù)少，數(shù)據(jù)處理速度快，為航空電磁數(shù)據(jù)調(diào)平提供了一種新思路。

航空電磁法;調(diào)平;條帶狀誤差;偽切割線法;導(dǎo)數(shù)

0 前言

航空電磁法(AEM，AirborneElectromagnetic Method)是基于巖石的電性、磁性等物性差異，應(yīng)用電磁感應(yīng)的測(cè)量原理，以飛機(jī)為測(cè)量載體，進(jìn)行探測(cè)的地球物理方法。它是一種快速、高效、經(jīng)濟(jì)、適用性強(qiáng)的勘探方法［1～6，9］。航空電磁法數(shù)據(jù)質(zhì)量和前期處理方法的好壞，對(duì)后期的異常反演、地質(zhì)解釋有著重大的影響［3～5，9］。隨著航空電磁法的發(fā)展，數(shù)據(jù)量越來(lái)越大，儀器也改進(jìn)了很多，信噪比得到較大提升。但這些卻無(wú)法徹底地消除漂移誤差，因此室內(nèi)數(shù)據(jù)調(diào)平，對(duì)提高數(shù)據(jù)質(zhì)量就顯得尤為重要［4、5］。

1991 年，周鳳桐［7～10］提出偽切割線調(diào)平方法，實(shí)現(xiàn)了人機(jī)交互調(diào)平。1994年，朱月娥［8］提出用條件濾波法消除航磁圖像測(cè)線干擾條帶，這種方法主要是結(jié)合中位數(shù)法和區(qū)間鄰點(diǎn)比較法判別異常極值點(diǎn)，能極好地保留有用的地質(zhì)信息。1999年，F(xiàn)UGRO公司的Huang［1］提出適用于DIGHEM吊艙式直升機(jī)頻率域電磁系統(tǒng)的二維自動(dòng)調(diào)平法，適用于零漂較小的直升機(jī)航電視電阻率數(shù)據(jù)。2007年，李文杰［9］提出適用于 HDY402三頻航電系統(tǒng)的二維移動(dòng)平均自動(dòng)調(diào)平方法，這種方法不需要網(wǎng)格化處理，方法簡(jiǎn)單，效率高。2008年 Huang［2］又提出基于線與線之間的航空頻率域測(cè)量數(shù)據(jù)的調(diào)平方法，取得良好的效果。

上述方法，在很大程度上提高了數(shù)據(jù)質(zhì)量。但是這些方法大都采用窗口濾波，而窗口參數(shù)和濾波參數(shù)很難選擇，需要大量嘗試才能取得比較好的效果，處理周期長(zhǎng)，而窗口濾波很難區(qū)分與調(diào)平誤差相近波長(zhǎng)的真實(shí)地質(zhì)信息［9～11］，并且窗口濾波還不可避免地存在邊界效應(yīng)［12］。

基于Huang［1］提出的在塊狀誤差判定時(shí)利用導(dǎo)數(shù)判定的原理，作者在本文提出基于導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)法的調(diào)平方法，減小了參數(shù)選取的難度，算法簡(jiǎn)單，效率高，并在極大程度上保留了真實(shí)地質(zhì)信息。該方法經(jīng)實(shí)際處理應(yīng)用，取得了良好的效果。

1 導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)法的原理

導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)法計(jì)算式為:

其中

其中 i表示網(wǎng)格化后的數(shù)據(jù)矩陣的某一行，其方向與偽切割線方向相同;j為數(shù)據(jù)矩陣的列，即沿測(cè)線方向;xi，j為原始數(shù)據(jù)網(wǎng)格化以后的數(shù)據(jù)矩陣元素;ˉxi，j為數(shù)據(jù)矩陣元素進(jìn)行一維濾波后得到的結(jié)果;N為數(shù)據(jù)矩陣的行數(shù)，即網(wǎng)格化后每條測(cè)線上數(shù)據(jù)點(diǎn)數(shù);M為數(shù)據(jù)矩陣的列數(shù)，即測(cè)線總數(shù);PT為高通濾波閾值，即導(dǎo)數(shù)值大于PT的點(diǎn)被認(rèn)為是導(dǎo)數(shù)突變點(diǎn);Nj為第j條測(cè)線上導(dǎo)數(shù)突變點(diǎn)的統(tǒng)計(jì)值;2l+1為一維濾波器長(zhǎng)度，其值如何選擇將在后面介紹。

根據(jù)式(4)，可判斷出存在零漂線號(hào):

其中 Lj為一行向量，其值為“1”的元素，所對(duì)應(yīng)的航線即為存在零漂異常的航線;LT為閾值，即只有在每條測(cè)線上導(dǎo)數(shù)突變點(diǎn)的統(tǒng)計(jì)個(gè)數(shù)大于閾值時(shí)，才認(rèn)為這條線存在零漂。

關(guān)于這兩個(gè)閾值如何選擇，作者將在后面介紹。

在確定了含有零漂的測(cè)線之后，我們采用均值濾波對(duì)其進(jìn)行處理:其中 x'i，j為調(diào)平處理后數(shù)據(jù);xi，j為網(wǎng)格化節(jié)點(diǎn)數(shù)據(jù);ˉxi，j為誤差線上點(diǎn)一維濾波結(jié)果;ˉxi，j－1、ˉxi，j+1分別為與之相鄰測(cè)線上的數(shù)據(jù)點(diǎn)，經(jīng)過(guò)一維濾波后的結(jié)果。

此均值濾波器使用方便，實(shí)現(xiàn)簡(jiǎn)單，運(yùn)算迅速、穩(wěn)定。

導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)法是利用統(tǒng)計(jì)學(xué)方法，來(lái)區(qū)分出含噪聲的測(cè)線。只要閾值選取適當(dāng)，就能準(zhǔn)確地找出含噪聲的測(cè)線。該方法公式簡(jiǎn)單，但效果明顯。

2 導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)法

在沿偽切割線方向，對(duì)每一條偽切割線進(jìn)行求導(dǎo)時(shí)。在導(dǎo)數(shù)曲線中存在很多突變點(diǎn)，突變點(diǎn)受系統(tǒng)零漂影響，也受大小不一的地質(zhì)起伏等因素影響。為了更好地突出系統(tǒng)零漂的影響，壓制地質(zhì)體影響，我們先沿測(cè)線方向?qū)?shù)據(jù)進(jìn)行一維濾波(見公式(2))，認(rèn)為濾波結(jié)果為背景值和系統(tǒng)零漂之和，不含地質(zhì)體影響。然后對(duì)一維濾波后的數(shù)據(jù)沿偽切割線方向進(jìn)行求導(dǎo)運(yùn)算，再利用高通濾波公式(1)，濾掉小的由殘余地質(zhì)起伏影響造成的導(dǎo)數(shù)突變點(diǎn)。

由于導(dǎo)數(shù)值較大的突變點(diǎn)，是由沿測(cè)線方向條帶誤差引起的，所以在高通濾波之后的導(dǎo)數(shù)突變點(diǎn)數(shù)據(jù)上，我們沿測(cè)線方向統(tǒng)計(jì)每一條測(cè)線的導(dǎo)數(shù)突變點(diǎn)數(shù)量。只要我們選取一個(gè)合理的統(tǒng)計(jì)點(diǎn)數(shù)量閾值，便能得到存在零漂誤差的測(cè)線(見公式(4))。對(duì)每一條誤差線進(jìn)行均值濾波(見公式(5))，最終得到調(diào)平結(jié)果，具體流程圖見下頁(yè)圖1。

在整個(gè)處理流程中，兩個(gè)閾值的選擇是該方法能夠?qū)崿F(xiàn)的關(guān)鍵:

(1)如果高通濾波器閾值過(guò)小，會(huì)引入很多小的導(dǎo)數(shù)突變點(diǎn)，給后面的突變點(diǎn)個(gè)數(shù)統(tǒng)計(jì)造成干擾。

(2)如果閾值過(guò)大，將遺漏零漂誤差線，造成的導(dǎo)數(shù)突變點(diǎn)，使濾波不徹底。

統(tǒng)計(jì)點(diǎn)個(gè)數(shù)閾值正確選擇的重要意義，與高通濾波器閾值的選擇類似。如果兩個(gè)閾值選擇合理，一次濾波便能達(dá)到目的。如果發(fā)現(xiàn)有含零漂誤差的測(cè)線遺漏，可進(jìn)行多次濾波。

閾值選取原則為:

(1)先將原始數(shù)據(jù)成圖，仔細(xì)觀察誤差條帶與相鄰條帶之間的數(shù)據(jù)差異，以此為高通閾值，或進(jìn)行微調(diào)。

圖1 導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)法處理流程圖Fig． 1 Flow chart for the derivative statistics methodleveling procedure

(2)統(tǒng)計(jì)閾值一般取值以線上點(diǎn)數(shù)的1/3到1/2為宜。

第一步的一維濾波也會(huì)對(duì)調(diào)平質(zhì)量產(chǎn)生一定的影響:

(1)濾波器長(zhǎng)度過(guò)長(zhǎng)，會(huì)產(chǎn)生較強(qiáng)的邊界效應(yīng)，影響數(shù)據(jù)質(zhì)量。

(2)濾波器長(zhǎng)度過(guò)短，不能很好地壓制地質(zhì)體的影響，對(duì)導(dǎo)數(shù)突變點(diǎn)的求取構(gòu)成一定的影響。

在一般情況下，一維濾波器的長(zhǎng)度選在測(cè)線測(cè)點(diǎn)數(shù)的1/20～1/10之間為宜。

我們把二維移動(dòng)平均自動(dòng)調(diào)平方法［9、12］與導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)法做一下對(duì)比。二維移動(dòng)平均自動(dòng)調(diào)平處理過(guò)程是在網(wǎng)格化數(shù)據(jù)的基礎(chǔ)上，進(jìn)行二維均值濾波，濾波結(jié)果得到區(qū)域背景值。原始網(wǎng)格數(shù)據(jù)與得到的這個(gè)背景值之差，便得到包含調(diào)平誤差和局部地質(zhì)異常的中間量。接著再對(duì)這個(gè)中間量沿測(cè)線方向進(jìn)行一維均值濾波，可以得到調(diào)平誤差量。原始網(wǎng)格數(shù)據(jù)減去調(diào)平誤差量，便得到最終調(diào)平結(jié)果。

3 數(shù)據(jù)處理實(shí)例

圖2為原始航空電磁法頻率870Hz實(shí)分量異常圖。從圖2中可以看出，存在很明顯的條帶狀誤差，這些誤差的存在，掩蓋了很多真實(shí)的地質(zhì)信息，如果不加以處理，將會(huì)給后面的處理解釋帶來(lái)很大的困難和誤差。

圖2 頻率為870Hz實(shí)分量航空電磁測(cè)量原始異常圖Fig.2 Raw anomaly map of 870 Hz in-phase component

針對(duì)這種情況，我們?cè)谶@里分別利用導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)法和二維平均自動(dòng)調(diào)平方法，對(duì)其進(jìn)行了處理，以作對(duì)照。在處理過(guò)程中，每條航線測(cè)點(diǎn)數(shù)為9967個(gè)，因此一維濾波器長(zhǎng)度選為599，高通濾波器閾值我們?nèi)≈禐?.5，統(tǒng)計(jì)點(diǎn)個(gè)數(shù)閾值取3300，即測(cè)線上數(shù)據(jù)點(diǎn)個(gè)數(shù)的三分之一，M/3。二維移動(dòng)平均處理二維窗口選為7*399，一維濾波窗口選為999。

經(jīng)處理后的結(jié)果平面圖如下頁(yè)圖3所示。其中，圖3(a)為二維移動(dòng)平均自動(dòng)調(diào)平方法處理結(jié)果，圖3(a)中還包含部份條帶誤差，圖像質(zhì)量較差。圖3(b)為導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)法處理，帶狀誤差已經(jīng)被濾掉，并且圖像比較光滑，且圖3(b)上地質(zhì)體處數(shù)值比二維移動(dòng)平均方法大，地質(zhì)信息損耗更少，效果要比二維移動(dòng)平均自動(dòng)調(diào)平方法處理得好。

圖4(見下頁(yè))為調(diào)平誤差，即原始數(shù)據(jù)與調(diào)平結(jié)果之間差值。其中，圖4(a)為二維移動(dòng)平均自動(dòng)調(diào)平法調(diào)平誤差，可以看出，條帶誤差被濾掉，部份區(qū)域有部份地質(zhì)異常被當(dāng)作誤差濾除了。所以這種方法很難區(qū)分與調(diào)平誤差相近波長(zhǎng)的真實(shí)地質(zhì)信息，造成調(diào)平過(guò)量。而圖4(b)為導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)方法的調(diào)平誤差，從圖4(b)中可以很清楚地看到，濾波只對(duì)有誤差的測(cè)線進(jìn)行，而無(wú)誤差的測(cè)線則沒有被改動(dòng)，并且無(wú)邊界效應(yīng)。由于圖4(b)中并沒包含明顯的地質(zhì)信息，說(shuō)明濾波較好地保留了原始地質(zhì)信息。

圖3 航空電磁法頻率870實(shí)分量調(diào)平結(jié)果Fig． 3 Leveling result of 870 Hz in-phase component

圖4 航電頻率870實(shí)分量調(diào)平誤差Fig． 4 Leveling error of 870 Hz in-phase component

在調(diào)平后圖3中，地質(zhì)體處數(shù)值要比二維平均自動(dòng)調(diào)平方法處理后的值要大一些，更接近實(shí)際地質(zhì)體異常信息。

導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)法運(yùn)用簡(jiǎn)單、高效，但是效果好與壞的關(guān)鍵，在于處理過(guò)程中閾值大小的選擇。對(duì)于處理經(jīng)驗(yàn)欠缺的操作者，在處理圖像之前，應(yīng)盡量選取幾條樣線，大致確定閾值范圍，以減少嘗試周期。

通過(guò)利用導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)法對(duì)實(shí)際數(shù)據(jù)進(jìn)行處理的結(jié)果，以及與二維移動(dòng)平均自動(dòng)調(diào)平法對(duì)比后發(fā)現(xiàn)，導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)法參數(shù)少，算法簡(jiǎn)單，運(yùn)算快捷。調(diào)平只針對(duì)存在誤差的測(cè)線，誤差更少，更準(zhǔn)確，并且不存在邊界效應(yīng)。

4 結(jié)論

針對(duì)航空電磁測(cè)量中存在的條帶狀測(cè)量誤差，作者在本文提出了導(dǎo)數(shù)突變點(diǎn)統(tǒng)計(jì)方法。該方法能根據(jù)數(shù)據(jù)的特點(diǎn)，進(jìn)行誤差消除處理。在方法的實(shí)施中，通過(guò)沿偽切割線方向(即垂直于飛行方向)求取沿線數(shù)據(jù)的導(dǎo)數(shù)，并統(tǒng)計(jì)導(dǎo)數(shù)變化較大點(diǎn)的數(shù)量，以判斷是否包含測(cè)量調(diào)平誤差。在獲得含噪聲的測(cè)線后，采用濾波方法對(duì)其進(jìn)行調(diào)平處理。

用本方法與基于窗口濾波的調(diào)平方法相比較可知，本方法算法簡(jiǎn)單，設(shè)置參數(shù)少，自動(dòng)化程度高，并能自動(dòng)識(shí)別含誤差測(cè)線，調(diào)平結(jié)果更加精確。

通過(guò)對(duì)本方法和二維移動(dòng)平均自動(dòng)調(diào)平方法的處理結(jié)果對(duì)比表明，本方法對(duì)于處理存在條帶狀誤差的數(shù)據(jù)具有較好的效果，并能有效地消除條帶誤差的影響，完整保留了微小的局部地質(zhì)體電磁異常。

［1］ HAOPING HUANG，DOUGLAS C． FRASER． Airborneresistivity data leveling［J］． GEOPHYSICS，1999，64( 2) : 378．

［2］ HAOPING HUANG． Airborne geophysical data levelingbased on line to line correlations ［J］． GEOPHYSICS，2008， 73( 3) : 83．

［3］ EIRIK MAURING， OLA KIHLE． Leveling aerogeophysicaldata using a moving differential median filter［J］． GEOPHYSICS，Jan － Feb， 2006， 71( 1) : L5．

［4］ MAJID BEIKI， MEHRDAD BASRANI， LAUST B．Pedersen． Leveling HEM and aeromagnetic data usingdifferential polynomial fitting［J］． GEOPHYSICS，Jan－ Feb，2010， 75( 1) : L13．

［5］ J． BRADLEY NELSON． Leveling total-field aerornagneticdata with measured horizontal gradients［J］． GEOPHYSICS，1994，59( 8) : 1166

［6］ 李文杰.頻率域航空電磁數(shù)據(jù)處理技術(shù)研究［D］.北京:中國(guó)地質(zhì)大學(xué)，博士學(xué)位論文，2008.

［7］ 周鳳桐，陳本池，閻永利.航空電磁法數(shù)據(jù)處理與圖示技術(shù)［J］.物探與化探，1997(21):348.

［8］ 朱月娥.用條件濾波法消除航磁圖像測(cè)線干擾條帶［J］.物探化探計(jì)算技術(shù)，1994，16(1):44.

［9］ 李文杰.用于頻率域航空電磁數(shù)據(jù)的二維自動(dòng)調(diào)平［J］.成都理工大學(xué)學(xué)報(bào)，2007，34(4):44.

［10］周鳳桐，閻永利，陳本池.航空電磁電阻率填圖方法技術(shù)及應(yīng)用研究［R］.河北廊坊:地質(zhì)礦產(chǎn)部地球物理地球化學(xué)勘查研究所，1991.

［11］孫東明，王衛(wèi)平，曾昭發(fā)，等.用于頻率域航電數(shù)據(jù)處理的偽切割線自動(dòng)調(diào)平法［J］.物探與化探，2010，34(2):246.

［12］孫東明，頻率域航空電磁法數(shù)據(jù)處理方法研究［D］.吉林大學(xué)，碩士論文，2010.

P631.3+26

A

1001—1749(2011)05—0522—05

2011－01－16 改回日期:2011－06－07

陳雄(1986－)，男，碩士，主要研究方向?yàn)楣こ膛c環(huán)境地球物理勘探，主要從事探地雷達(dá)正演模擬及電磁法算法研究。