用于航空電磁數據調平處理的導數突變點統計法

陳 雄，曾昭發，孫東明，李 靜

(1.吉林大學，吉林長春 130026;2.吉林省公路勘測設計院，吉林長春 130026)

用于航空電磁數據調平處理的導數突變點統計法

陳 雄1，曾昭發1，孫東明2，李 靜1

(1.吉林大學，吉林長春 130026;2.吉林省公路勘測設計院，吉林長春 130026)

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

航空電磁法;調平;條帶狀誤差;偽切割線法;導數

0 前言

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

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

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

基于Huang［1］提出的在塊狀誤差判定時利用導數判定的原理，作者在本文提出基于導數突變點統計法的調平方法，減小了參數選取的難度，算法簡單，效率高，并在極大程度上保留了真實地質信息。該方法經實際處理應用，取得了良好的效果。

1 導數突變點統計法的原理

導數突變點統計法計算式為:

其中

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

根據式(4)，可判斷出存在零漂線號:

其中 Lj為一行向量，其值為“1”的元素，所對應的航線即為存在零漂異常的航線;LT為閾值，即只有在每條測線上導數突變點的統計個數大于閾值時，才認為這條線存在零漂。

關于這兩個閾值如何選擇，作者將在后面介紹。

在確定了含有零漂的測線之后，我們采用均值濾波對其進行處理:其中 x'i，j為調平處理后數據;xi，j為網格化節點數據;ˉxi，j為誤差線上點一維濾波結果;ˉxi，j－1、ˉxi，j+1分別為與之相鄰測線上的數據點，經過一維濾波后的結果。

此均值濾波器使用方便，實現簡單，運算迅速、穩定。

導數突變點統計法是利用統計學方法，來區分出含噪聲的測線。只要閾值選取適當，就能準確地找出含噪聲的測線。該方法公式簡單，但效果明顯。

2 導數突變點統計法

在沿偽切割線方向，對每一條偽切割線進行求導時。在導數曲線中存在很多突變點，突變點受系統零漂影響，也受大小不一的地質起伏等因素影響。為了更好地突出系統零漂的影響，壓制地質體影響，我們先沿測線方向對數據進行一維濾波(見公式(2))，認為濾波結果為背景值和系統零漂之和，不含地質體影響。然后對一維濾波后的數據沿偽切割線方向進行求導運算，再利用高通濾波公式(1)，濾掉小的由殘余地質起伏影響造成的導數突變點。

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

在整個處理流程中，兩個閾值的選擇是該方法能夠實現的關鍵:

(1)如果高通濾波器閾值過小，會引入很多小的導數突變點，給后面的突變點個數統計造成干擾。

(2)如果閾值過大，將遺漏零漂誤差線，造成的導數突變點，使濾波不徹底。

統計點個數閾值正確選擇的重要意義，與高通濾波器閾值的選擇類似。如果兩個閾值選擇合理，一次濾波便能達到目的。如果發現有含零漂誤差的測線遺漏，可進行多次濾波。

閾值選取原則為:

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

圖1 導數突變點統計法處理流程圖Fig． 1 Flow chart for the derivative statistics methodleveling procedure

(2)統計閾值一般取值以線上點數的1/3到1/2為宜。

第一步的一維濾波也會對調平質量產生一定的影響:

(1)濾波器長度過長，會產生較強的邊界效應，影響數據質量。

(2)濾波器長度過短，不能很好地壓制地質體的影響，對導數突變點的求取構成一定的影響。

在一般情況下，一維濾波器的長度選在測線測點數的1/20～1/10之間為宜。

我們把二維移動平均自動調平方法［9、12］與導數突變點統計法做一下對比。二維移動平均自動調平處理過程是在網格化數據的基礎上，進行二維均值濾波，濾波結果得到區域背景值。原始網格數據與得到的這個背景值之差，便得到包含調平誤差和局部地質異常的中間量。接著再對這個中間量沿測線方向進行一維均值濾波，可以得到調平誤差量。原始網格數據減去調平誤差量，便得到最終調平結果。

3 數據處理實例

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

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

針對這種情況，我們在這里分別利用導數突變點統計法和二維平均自動調平方法，對其進行了處理，以作對照。在處理過程中，每條航線測點數為9967個，因此一維濾波器長度選為599，高通濾波器閾值我們取值為2.5，統計點個數閾值取3300，即測線上數據點個數的三分之一，M/3。二維移動平均處理二維窗口選為7*399，一維濾波窗口選為999。

經處理后的結果平面圖如下頁圖3所示。其中，圖3(a)為二維移動平均自動調平方法處理結果，圖3(a)中還包含部份條帶誤差，圖像質量較差。圖3(b)為導數突變點統計法處理，帶狀誤差已經被濾掉，并且圖像比較光滑，且圖3(b)上地質體處數值比二維移動平均方法大，地質信息損耗更少，效果要比二維移動平均自動調平方法處理得好。

圖4(見下頁)為調平誤差，即原始數據與調平結果之間差值。其中，圖4(a)為二維移動平均自動調平法調平誤差，可以看出，條帶誤差被濾掉，部份區域有部份地質異常被當作誤差濾除了。所以這種方法很難區分與調平誤差相近波長的真實地質信息，造成調平過量。而圖4(b)為導數突變點統計方法的調平誤差，從圖4(b)中可以很清楚地看到，濾波只對有誤差的測線進行，而無誤差的測線則沒有被改動，并且無邊界效應。由于圖4(b)中并沒包含明顯的地質信息，說明濾波較好地保留了原始地質信息。

圖3 航空電磁法頻率870實分量調平結果Fig． 3 Leveling result of 870 Hz in-phase component

圖4 航電頻率870實分量調平誤差Fig． 4 Leveling error of 870 Hz in-phase component

在調平后圖3中，地質體處數值要比二維平均自動調平方法處理后的值要大一些，更接近實際地質體異常信息。

導數突變點統計法運用簡單、高效，但是效果好與壞的關鍵，在于處理過程中閾值大小的選擇。對于處理經驗欠缺的操作者，在處理圖像之前，應盡量選取幾條樣線，大致確定閾值范圍，以減少嘗試周期。

通過利用導數突變點統計法對實際數據進行處理的結果，以及與二維移動平均自動調平法對比后發現，導數突變點統計法參數少，算法簡單，運算快捷。調平只針對存在誤差的測線，誤差更少，更準確，并且不存在邊界效應。

4 結論

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

用本方法與基于窗口濾波的調平方法相比較可知，本方法算法簡單，設置參數少，自動化程度高，并能自動識別含誤差測線，調平結果更加精確。

通過對本方法和二維移動平均自動調平方法的處理結果對比表明，本方法對于處理存在條帶狀誤差的數據具有較好的效果，并能有效地消除條帶誤差的影響，完整保留了微小的局部地質體電磁異常。

［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］ 李文杰.頻率域航空電磁數據處理技術研究［D］.北京:中國地質大學，博士學位論文，2008.

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

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

［9］ 李文杰.用于頻率域航空電磁數據的二維自動調平［J］.成都理工大學學報，2007，34(4):44.

［10］周鳳桐，閻永利，陳本池.航空電磁電阻率填圖方法技術及應用研究［R］.河北廊坊:地質礦產部地球物理地球化學勘查研究所，1991.

［11］孫東明，王衛平，曾昭發，等.用于頻率域航電數據處理的偽切割線自動調平法［J］.物探與化探，2010，34(2):246.

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

P631.3+26

A

1001—1749(2011)05—0522—05

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

陳雄(1986－)，男，碩士，主要研究方向為工程與環境地球物理勘探，主要從事探地雷達正演模擬及電磁法算法研究。