李文杰

(中國地質科學院 地球物理地球化學勘查研究所， 廊坊 065000)

頻率域航空電磁法在進行礦產勘查時沿測線連續(xù)測量，平均每15 m一個數(shù)據，測線長度十幾千米至上百千米，一個測區(qū)測線總量從幾百測線千米到幾十萬測線千米不等。一般每測線千米數(shù)據中平均會出現(xiàn)2 處～3處明顯的異常，一個測區(qū)則會產生數(shù)千至數(shù)十萬個異常。在如此巨量的航空電磁異常中，排除飛行過載、工業(yè)電、飛行高度變化、地質噪聲等干擾引起的假異常，快速和高效地識別和篩選出有意義的基巖異常，是有效發(fā)揮航空電磁法快速高效勘查優(yōu)勢的關鍵技術。

作者基于中國地質科學院地球物理地球化學勘查研究所(簡稱物化探所，IGGE)研制的國內唯一自主技術固定翼頻率域航空電磁系統(tǒng)—HDY-402三頻航電系統(tǒng)的特有技術，參考國內、外相關技術，借鑒傳統(tǒng)的人工異常識別篩選方法，開發(fā)了計算機軟件異常識別、篩選和解釋算法(Computer-aided Anomaly Detection, CAD)，并編制了應用軟件，應用在內蒙古地區(qū)開展的約42.3萬測線千米航空物探綜合勘查項目中，有效地提高了海量異常識別篩選的工作效率和質量。

1 國內外研究情況

國外在上世紀七十年代初，針對固定翼頻率航電數(shù)據，開發(fā)了基于計算機的交互式異常識別技術[1-2]，提高了異常解釋效率。九十年代末，澳大利亞礦產堪查工業(yè)聯(lián)合研究中心[3](CRC AMET)開發(fā)了EM FLOW數(shù)據處理解釋軟件，可以在時間常數(shù)(tau)域內根據設定的參數(shù)實現(xiàn)航電異常的快速自動識別和反演計算。Valleau[4]系統(tǒng)地介紹了DIGHEM直升機頻率域航空電磁數(shù)據的處理流程，實現(xiàn)了航電異常的計算機提取、篩選及人工輔助判別。Huang[5-6]提出用多個頻率電磁耦合譜的不同特征來識別頻率域航電異常，用于UXO探測。

在國內，八十年代初，周鳳桐[7]總結了固定翼頻率域航空電磁資料的人工處理流程，李文杰[8]研究了固定翼頻率域航空電磁異常軟件識別技術。在航空電磁異常的軟件數(shù)字化自動快速拾取和識別方面，國內未見其他相關報導。

2 方法步驟

這里討論的航電異常數(shù)字化識別篩選，主要是指識別局部異常，篩選出其中的有意義的良導基巖異常。異常的計算機自動識別和交互判別，簡單地說就是將人工識別經驗和技術數(shù)字化為計算機標準識別參數(shù)和程序代碼。航電異常軟件自動識別的過程與人工異常識別大體一致，需要經過如下幾個步驟：

1)異常自動拾?。杭催x出一般意義上的待選異常。

2)異常識別篩選：用約束條件來剔除待選異常中明顯無意義的異常。

3)交互輔助判別：進一步剔除地質噪聲、地形效應異常等難于自動識別的異常。

2.1 異常自動拾取

對經過噪聲濾波、相位改正、零漂改正等處理的三頻航電剖面數(shù)據，取梯度后判別拾取出曲線波峰、波谷，再標記出所有大于三倍噪聲水平的待選異常；然后提取異常信息，計算每個異常的半極值寬度、異常幅值、實虛分量比、視電阻率、視深度、飛行高度、地形及地形匹配系數(shù)、工業(yè)電干擾、過載等相關信息。這些數(shù)字化了的異常信息，是下一步自動識別的基礎信息。

假設航電某頻道曲線為Eij(i為測線號，j為測點號)(圖1)，則它的梯度為

Δe(i,j)=grad(Eij)=E(i,j+1)-E(i,j)

(1)

如果Δe(i,j)=0，則在Eij上出現(xiàn)一個拐點，當在拐點相鄰點的梯度滿足式(2)時

(2)

則拐點對應E(i,j)曲線上的一個谷點(如圖1中的點A、C)。同理，當拐點相鄰點的梯度滿足：

(3)

則拐點對應E(i,j)曲線上的一個峰點(如圖1中的點P)，一個波峰加相鄰兩個波谷就是一個異常。

異常拾取過程中，附帶地可以計算出異常的半極值寬度Pw、異常的峰值高度Pp、實虛分量比，以及峰值點處的視電阻率、視電阻率峰值等信息，以備后續(xù)步驟篩選。

2.2 自動異常識別

根據基巖異常的特征，設定異常峰值、半極值寬度、背景場值、飛行高度、多頻道交集、視電阻率峰值、覆蓋層判別等篩選條件，對待選異常進行自動篩選識別，給出異常分類。

雖然數(shù)據預處理時進行了噪聲濾波，但系統(tǒng)噪聲在異常篩選階段依然需要被考慮。由于噪聲的存在，異常峰值高度Pp一般會有一個最小限制即Ppmin≥3·vN，vN是系統(tǒng)的噪聲水平，即當異常的峰值高度Pp大于三倍噪聲水平時才被考慮。

圖1 異常的拾取原理

1)半極值寬度。大多數(shù)礦體規(guī)模比較小，尤其硫化礦體，其走向長度大約在200 m～300 m以下的居多。即使是規(guī)模比較大的礦體(如磁鐵礦體或大型硫化礦體)，其走向長度可以大于500 m ～600 m，但往往由于地形或者礦體本身傾伏等原因，并不總是能使礦體整個地位于系統(tǒng)的探測深度范圍之內，因而實際上反映在航電異常曲線的半極值寬度一般都小于500 m。

2)異常背景。有意義的礦體異常一般不會出現(xiàn)在大面積的低阻區(qū)，因此異常波谷的場值可以作為異常背景加以約束，波谷谷點處的場值大于某個限制后，異常可以不被考慮。

3)測量飛行高度。航電接收的二次場隨著測量高度的增加以三次冪或更大的速率衰減。測量高度超高時，航電響應中噪聲、零漂等干擾因素所占比重成倍增加，測得異常的可信度很低。根據實踐經驗，一般可以取70 m～80 m作為異常峰值點測量高度的上限。

4)工業(yè)電干擾。航電接收的二次場響應異常有一部分是由于50 Hz工業(yè)電干擾及飛行過載引起的，這類異?？梢岳幂o助設備提供的輔助參數(shù)來判別和剔除。一般航電系統(tǒng)50 Hz監(jiān)測輸出的信號小于100 mV；當測線上有工業(yè)電高壓線等干擾時，50 Hz監(jiān)測輸出的信號會遠大于100 mV，并隨干擾強度的變化而改變。

5)過載干擾。系統(tǒng)自帶的過載計在飛行相對平穩(wěn)時的輸出通常小于±0.15 g，當飛機遇到氣流或突然的機動動作會使過載計輸出產生正負變化，當過載輸出大于±0.15 g時所引起的異常可以判斷為過載異常，并加以剔除。

6)多頻道信息。HDY-402航電系統(tǒng)共有463 Hz、1 563 Hz、8 333 Hz三個工作頻率，相對而言中頻和低頻的探測深度較大，因此在判斷異常時主要以中低頻為主、高頻為輔。如果中低頻同時具有異常，則是良導基巖異常的可能性較大。結合低、中、高頻的視電阻率、視飛行高度數(shù)據，依據視電阻率與二層大地模型的關系，還可以推斷出二層大地的電阻率高低，對覆蓋層異常的判別幫助很大，也可以作為一類約束條件。

在異常自動識別運算過程中，約束條件如果采用上限和下限雙邊限制，結合特別給定的約束條件組合，就可以對某類特定類型的異常進行識別，而不光是良導基巖異常。

舉例說明，如果將半極值寬度、峰值高度、異常背景、測量高度等約束條件的上限和下限都放寬，即相當于不限制這些條件，而只將輔助信息中的50 Hz工業(yè)干擾下限設得較小，如小于100 mV，上限放寬，則可以得到滿足大于100 mV下限標準的的50 Hz工業(yè)電干擾引起的干擾異常。

2.3 交互輔助識別

數(shù)據解釋人員根據異常信息，參考地質、化探等相關資料，對自動識別拾取的異常逐個進行解釋推斷，剔除地形效應、負地形、覆蓋層等地質噪聲引起的假異常，以及過載、工業(yè)電等引起的干擾異常，對推斷的基巖異常進行標識。同時對未被自動識別拾取的異常進行補查，以免遺漏可能的基巖異常；最后將所有已標識的異常信息進行輸出和匯總。

圖2是兩處在內蒙二連浩特測區(qū)用軟件自動識別出來的航電典型異常。圖2中曲線自上至下依次為異常響應的視電阻率、高頻8 333 Hz(HI/HQ)、中頻1 563 Hz(MI/MQ)、低頻463 Hz(LI/LQ)實虛分量、飛行高度、地形等剖面曲線。圖2(a)為三頻航電在一處已知熱液型銅礦化點上的異常響應曲線，對應有明顯的中低頻視電阻率異常。圖2(b)為一處典型的殘坡積(低洼地)引起的地表覆蓋層異常，對應有高頻視電阻率異常，中低頻視電阻率無明顯異常。

2.4 異常自動識別的軟件流程

1)準備數(shù)據。用于進行異常識別的航電數(shù)據需經過噪聲濾波處理、相位改正及零漂改正。噪聲濾波處理是必須的，其意不言自明。未進行相位改正和零漂改正的航電數(shù)據，在進行異常信息提取時，將得到錯誤的視電阻率、視深度結果，無法保證異常背景、覆蓋層異常判定等約束條件的有效性。

2)設定異常篩選條件。設定異常峰值、半極值寬度、背景、測量高度、多頻道交集、視電阻率峰值、覆蓋層判別等篩選條件。上述異常篩選的約束條件因為都用上、下限數(shù)值賦值，因此可以交集使用，或通過將某一項或幾項約束條件的上、下限數(shù)值設定為完全滿足而不參與異常篩選的條件約束。

3)異常拾取。計算航電各頻道實、虛分量曲線的梯度，根據正、負峰組合判斷拾取波峰、波谷，一個波峰及相鄰兩個波谷構成一個被選異常。

4)異常信息提取。計算被選異常的峰值、半極值寬度、視電阻率、視電阻率背景、測量高度、多頻道交集、視電阻率峰值、視飛行高度等信息，并根據半極值寬度等條件識別雙峰異常。

圖2 已知礦化點異常和殘坡積引起的覆蓋層異常綜合剖面圖Fig．2 Composite profiles of an ore anomaly and a residual sediment anomaly (a)礦化點異常；(b)覆蓋層異常

5)應用約束條件篩選異常。應用異常篩選條件，篩選出符合各約束條件交集的異常。如上所述，約束條件可根據任務需要進行設定，因而可以區(qū)分出50Hz工業(yè)電干擾異常、覆蓋層異常等特定類型異常。

6)輸出異常篩選結果。將異常結果輸出為文本文件，其中包含異常的各種信息，及位置坐標、異常種類、評級等，供方法解釋人員進一步篩選。

7)輔助信息解釋。根據地形剖面曲線等輔助條件對異常篩選結果進行進一步分類，剔除負地形異常及地形效應異常。

將異常識別算法用軟件代碼實現(xiàn)，然后再依據以上流程就可以編織出航電異常識別程序(圖3)。

圖3 異常自動識別流程Fig．3 Flow chart of the digital anomaly detection and classification

3 應用實例

圖4為內蒙古某測區(qū)航電異常自動拾取的結果。圖中彩色底圖為視電阻率平面圖(暖色代表視低電阻率，冷色代表高視電阻率)，紫色三角點為自動拾取的異常位置，黑色箭頭所指藍色圓點為已知礦或礦化點?？梢钥吹皆趫D中間最上部①、④號已知礦化點與自動異常拾取的位置是完全對應的。①號異常為已知的熱液型銅礦化點，自動識別在8 040線上有明顯航電響應異常與之對應(圖5)。該異常控制長度2 km，半極值寬度297 m，飛行高度48 m，中頻和低頻實虛分量比分別為1.88和1.76，異常清晰。④號已知熱液型銅礦化點在8 170、8 180線間穿過，在8 180線礦化點附近自動識別找到了航電異常，但異常峰值較小。②號已知熱液型銅礦化點在6 870線上方通過，但測線上未見明顯異常，在相鄰的6 850線、6 860線航電異常明顯。③號已知熱液型錳礦化點位于8 240線、8 250線間，為熱液型銅礦化點，在測線上未發(fā)現(xiàn)明顯異常，但在相鄰的8 270、8 280、8 290線上自動識別出了異常條帶。

圖4 內蒙某測區(qū)航電異常自動識別篩選結果Fig．4 Digital detection and classification results of a survey area in Inner Mongolia

圖5 內蒙某測區(qū)8040線上的已知礦上航電異常響應

4 結論

頻率域航空電磁異常的數(shù)字化自動識別篩選技術補充完善了國產HDY-402三頻航空電磁系統(tǒng)數(shù)據處理技術，開發(fā)的異常識別程序實用、高效，可根據解釋人員定義的不同約束條件檢出不同類型的航電異常，提高了數(shù)據處理的效率和質量，填補了國內該領域的技術空白。但是受到數(shù)據質量不佳及數(shù)據信息量有限等條件的限制，航電異常的自動識別篩選結果中包含有干擾異常，僅能實現(xiàn)異常的初步篩選，還不能完全替代解釋人員的人工判別。

參考文獻：

[1] DAVID FOUNTAIN, RICHARD SMITH. 55 Years of AEM; 50 Years of KEGS! Fugro Web Site, 2003, www.fugro.com.au.

[2] PALACKY G J, WEST G F. Computer Processing of Airborne EM Data[J]. Geophysical Prospecting, 1974, 22:490-509.

[3] EM FLOW. User Guide and Reference Manua[M]. Encom Technology Pty Ltd.

[4] NICHOLAS VALLEAU. HEM Data Processing-A Practical Overview[M]. Proc. 14th ASEG Conf and Exhibition, 2000, Perth, Australia.

[5] HAOPING HUANG, WON I J. Identification of Mineral Deposits using airborne electromagnetic spectra[M]. SEG Annual Meeting, 2002.

[6] HOAPING HUANG. Automated anomaly picking from broadband electromagnetic data in an unexploded ordance (UXO) survey [J] .Geophysics, 2003, 68(6), p1870-1876.

[7] 周鳳桐. 補償式航空電磁法[J]. 物化探研究報道，1982(11):1-78.

[8] 李文杰. 頻率域航空電磁數(shù)據處理技術研究 [D] .北京：中國地質大學, 2008.