李文杰
(中國地質(zhì)科學(xué)院 地球物理地球化學(xué)勘查研究所, 廊坊 065000)
頻率域航空電磁法在進(jìn)行礦產(chǎn)勘查時(shí)沿測(cè)線連續(xù)測(cè)量,平均每15 m一個(gè)數(shù)據(jù),測(cè)線長度十幾千米至上百千米,一個(gè)測(cè)區(qū)測(cè)線總量從幾百測(cè)線千米到幾十萬測(cè)線千米不等。一般每測(cè)線千米數(shù)據(jù)中平均會(huì)出現(xiàn)2 處~3處明顯的異常,一個(gè)測(cè)區(qū)則會(huì)產(chǎn)生數(shù)千至數(shù)十萬個(gè)異常。在如此巨量的航空電磁異常中,排除飛行過載、工業(yè)電、飛行高度變化、地質(zhì)噪聲等干擾引起的假異常,快速和高效地識(shí)別和篩選出有意義的基巖異常,是有效發(fā)揮航空電磁法快速高效勘查優(yōu)勢(shì)的關(guān)鍵技術(shù)。
作者基于中國地質(zhì)科學(xué)院地球物理地球化學(xué)勘查研究所(簡稱物化探所,IGGE)研制的國內(nèi)唯一自主技術(shù)固定翼頻率域航空電磁系統(tǒng)—HDY-402三頻航電系統(tǒng)的特有技術(shù),參考國內(nèi)、外相關(guān)技術(shù),借鑒傳統(tǒng)的人工異常識(shí)別篩選方法,開發(fā)了計(jì)算機(jī)軟件異常識(shí)別、篩選和解釋算法(Computer-aided Anomaly Detection, CAD),并編制了應(yīng)用軟件,應(yīng)用在內(nèi)蒙古地區(qū)開展的約42.3萬測(cè)線千米航空物探綜合勘查項(xiàng)目中,有效地提高了海量異常識(shí)別篩選的工作效率和質(zhì)量。
國外在上世紀(jì)七十年代初,針對(duì)固定翼頻率航電數(shù)據(jù),開發(fā)了基于計(jì)算機(jī)的交互式異常識(shí)別技術(shù)[1-2],提高了異常解釋效率。九十年代末,澳大利亞礦產(chǎn)堪查工業(yè)聯(lián)合研究中心[3](CRC AMET)開發(fā)了EM FLOW數(shù)據(jù)處理解釋軟件,可以在時(shí)間常數(shù)(tau)域內(nèi)根據(jù)設(shè)定的參數(shù)實(shí)現(xiàn)航電異常的快速自動(dòng)識(shí)別和反演計(jì)算。Valleau[4]系統(tǒng)地介紹了DIGHEM直升機(jī)頻率域航空電磁數(shù)據(jù)的處理流程,實(shí)現(xiàn)了航電異常的計(jì)算機(jī)提取、篩選及人工輔助判別。Huang[5-6]提出用多個(gè)頻率電磁耦合譜的不同特征來識(shí)別頻率域航電異常,用于UXO探測(cè)。
在國內(nèi),八十年代初,周鳳桐[7]總結(jié)了固定翼頻率域航空電磁資料的人工處理流程,李文杰[8]研究了固定翼頻率域航空電磁異常軟件識(shí)別技術(shù)。在航空電磁異常的軟件數(shù)字化自動(dòng)快速拾取和識(shí)別方面,國內(nèi)未見其他相關(guān)報(bào)導(dǎo)。
這里討論的航電異常數(shù)字化識(shí)別篩選,主要是指識(shí)別局部異常,篩選出其中的有意義的良導(dǎo)基巖異常。異常的計(jì)算機(jī)自動(dòng)識(shí)別和交互判別,簡單地說就是將人工識(shí)別經(jīng)驗(yàn)和技術(shù)數(shù)字化為計(jì)算機(jī)標(biāo)準(zhǔn)識(shí)別參數(shù)和程序代碼。航電異常軟件自動(dòng)識(shí)別的過程與人工異常識(shí)別大體一致,需要經(jīng)過如下幾個(gè)步驟:
1)異常自動(dòng)拾取:即選出一般意義上的待選異常。
2)異常識(shí)別篩選:用約束條件來剔除待選異常中明顯無意義的異常。
3)交互輔助判別:進(jìn)一步剔除地質(zhì)噪聲、地形效應(yīng)異常等難于自動(dòng)識(shí)別的異常。
對(duì)經(jīng)過噪聲濾波、相位改正、零漂改正等處理的三頻航電剖面數(shù)據(jù),取梯度后判別拾取出曲線波峰、波谷,再標(biāo)記出所有大于三倍噪聲水平的待選異常;然后提取異常信息,計(jì)算每個(gè)異常的半極值寬度、異常幅值、實(shí)虛分量比、視電阻率、視深度、飛行高度、地形及地形匹配系數(shù)、工業(yè)電干擾、過載等相關(guān)信息。這些數(shù)字化了的異常信息,是下一步自動(dòng)識(shí)別的基礎(chǔ)信息。
假設(shè)航電某頻道曲線為Eij(i為測(cè)線號(hào),j為測(cè)點(diǎn)號(hào))(圖1),則它的梯度為
Δe(i,j)=grad(Eij)=E(i,j+1)-E(i,j)
(1)
如果Δe(i,j)=0,則在Eij上出現(xiàn)一個(gè)拐點(diǎn),當(dāng)在拐點(diǎn)相鄰點(diǎn)的梯度滿足式(2)時(shí)
(2)
則拐點(diǎn)對(duì)應(yīng)E(i,j)曲線上的一個(gè)谷點(diǎn)(如圖1中的點(diǎn)A、C)。同理,當(dāng)拐點(diǎn)相鄰點(diǎn)的梯度滿足:
(3)
則拐點(diǎn)對(duì)應(yīng)E(i,j)曲線上的一個(gè)峰點(diǎn)(如圖1中的點(diǎn)P),一個(gè)波峰加相鄰兩個(gè)波谷就是一個(gè)異常。
異常拾取過程中,附帶地可以計(jì)算出異常的半極值寬度Pw、異常的峰值高度Pp、實(shí)虛分量比,以及峰值點(diǎn)處的視電阻率、視電阻率峰值等信息,以備后續(xù)步驟篩選。
根據(jù)基巖異常的特征,設(shè)定異常峰值、半極值寬度、背景場(chǎng)值、飛行高度、多頻道交集、視電阻率峰值、覆蓋層判別等篩選條件,對(duì)待選異常進(jìn)行自動(dòng)篩選識(shí)別,給出異常分類。
雖然數(shù)據(jù)預(yù)處理時(shí)進(jìn)行了噪聲濾波,但系統(tǒng)噪聲在異常篩選階段依然需要被考慮。由于噪聲的存在,異常峰值高度Pp一般會(huì)有一個(gè)最小限制即Ppmin≥3·vN,vN是系統(tǒng)的噪聲水平,即當(dāng)異常的峰值高度Pp大于三倍噪聲水平時(shí)才被考慮。
圖1 異常的拾取原理
1)半極值寬度。大多數(shù)礦體規(guī)模比較小,尤其硫化礦體,其走向長度大約在200 m~300 m以下的居多。即使是規(guī)模比較大的礦體(如磁鐵礦體或大型硫化礦體),其走向長度可以大于500 m ~600 m,但往往由于地形或者礦體本身傾伏等原因,并不總是能使礦體整個(gè)地位于系統(tǒng)的探測(cè)深度范圍之內(nèi),因而實(shí)際上反映在航電異常曲線的半極值寬度一般都小于500 m。
2)異常背景。有意義的礦體異常一般不會(huì)出現(xiàn)在大面積的低阻區(qū),因此異常波谷的場(chǎng)值可以作為異常背景加以約束,波谷谷點(diǎn)處的場(chǎng)值大于某個(gè)限制后,異??梢圆槐豢紤]。
3)測(cè)量飛行高度。航電接收的二次場(chǎng)隨著測(cè)量高度的增加以三次冪或更大的速率衰減。測(cè)量高度超高時(shí),航電響應(yīng)中噪聲、零漂等干擾因素所占比重成倍增加,測(cè)得異常的可信度很低。根據(jù)實(shí)踐經(jīng)驗(yàn),一般可以取70 m~80 m作為異常峰值點(diǎn)測(cè)量高度的上限。
4)工業(yè)電干擾。航電接收的二次場(chǎng)響應(yīng)異常有一部分是由于50 Hz工業(yè)電干擾及飛行過載引起的,這類異常可以利用輔助設(shè)備提供的輔助參數(shù)來判別和剔除。一般航電系統(tǒng)50 Hz監(jiān)測(cè)輸出的信號(hào)小于100 mV;當(dāng)測(cè)線上有工業(yè)電高壓線等干擾時(shí),50 Hz監(jiān)測(cè)輸出的信號(hào)會(huì)遠(yuǎn)大于100 mV,并隨干擾強(qiáng)度的變化而改變。
5)過載干擾。系統(tǒng)自帶的過載計(jì)在飛行相對(duì)平穩(wěn)時(shí)的輸出通常小于±0.15 g,當(dāng)飛機(jī)遇到氣流或突然的機(jī)動(dòng)動(dòng)作會(huì)使過載計(jì)輸出產(chǎn)生正負(fù)變化,當(dāng)過載輸出大于±0.15 g時(shí)所引起的異??梢耘袛酁檫^載異常,并加以剔除。
6)多頻道信息。HDY-402航電系統(tǒng)共有463 Hz、1 563 Hz、8 333 Hz三個(gè)工作頻率,相對(duì)而言中頻和低頻的探測(cè)深度較大,因此在判斷異常時(shí)主要以中低頻為主、高頻為輔。如果中低頻同時(shí)具有異常,則是良導(dǎo)基巖異常的可能性較大。結(jié)合低、中、高頻的視電阻率、視飛行高度數(shù)據(jù),依據(jù)視電阻率與二層大地模型的關(guān)系,還可以推斷出二層大地的電阻率高低,對(duì)覆蓋層異常的判別幫助很大,也可以作為一類約束條件。
在異常自動(dòng)識(shí)別運(yùn)算過程中,約束條件如果采用上限和下限雙邊限制,結(jié)合特別給定的約束條件組合,就可以對(duì)某類特定類型的異常進(jìn)行識(shí)別,而不光是良導(dǎo)基巖異常。
舉例說明,如果將半極值寬度、峰值高度、異常背景、測(cè)量高度等約束條件的上限和下限都放寬,即相當(dāng)于不限制這些條件,而只將輔助信息中的50 Hz工業(yè)干擾下限設(shè)得較小,如小于100 mV,上限放寬,則可以得到滿足大于100 mV下限標(biāo)準(zhǔn)的的50 Hz工業(yè)電干擾引起的干擾異常。
數(shù)據(jù)解釋人員根據(jù)異常信息,參考地質(zhì)、化探等相關(guān)資料,對(duì)自動(dòng)識(shí)別拾取的異常逐個(gè)進(jìn)行解釋推斷,剔除地形效應(yīng)、負(fù)地形、覆蓋層等地質(zhì)噪聲引起的假異常,以及過載、工業(yè)電等引起的干擾異常,對(duì)推斷的基巖異常進(jìn)行標(biāo)識(shí)。同時(shí)對(duì)未被自動(dòng)識(shí)別拾取的異常進(jìn)行補(bǔ)查,以免遺漏可能的基巖異常;最后將所有已標(biāo)識(shí)的異常信息進(jìn)行輸出和匯總。
圖2是兩處在內(nèi)蒙二連浩特測(cè)區(qū)用軟件自動(dòng)識(shí)別出來的航電典型異常。圖2中曲線自上至下依次為異常響應(yīng)的視電阻率、高頻8 333 Hz(HI/HQ)、中頻1 563 Hz(MI/MQ)、低頻463 Hz(LI/LQ)實(shí)虛分量、飛行高度、地形等剖面曲線。圖2(a)為三頻航電在一處已知熱液型銅礦化點(diǎn)上的異常響應(yīng)曲線,對(duì)應(yīng)有明顯的中低頻視電阻率異常。圖2(b)為一處典型的殘坡積(低洼地)引起的地表覆蓋層異常,對(duì)應(yīng)有高頻視電阻率異常,中低頻視電阻率無明顯異常。
1)準(zhǔn)備數(shù)據(jù)。用于進(jìn)行異常識(shí)別的航電數(shù)據(jù)需經(jīng)過噪聲濾波處理、相位改正及零漂改正。噪聲濾波處理是必須的,其意不言自明。未進(jìn)行相位改正和零漂改正的航電數(shù)據(jù),在進(jìn)行異常信息提取時(shí),將得到錯(cuò)誤的視電阻率、視深度結(jié)果,無法保證異常背景、覆蓋層異常判定等約束條件的有效性。
2)設(shè)定異常篩選條件。設(shè)定異常峰值、半極值寬度、背景、測(cè)量高度、多頻道交集、視電阻率峰值、覆蓋層判別等篩選條件。上述異常篩選的約束條件因?yàn)槎加蒙?、下限?shù)值賦值,因此可以交集使用,或通過將某一項(xiàng)或幾項(xiàng)約束條件的上、下限數(shù)值設(shè)定為完全滿足而不參與異常篩選的條件約束。
3)異常拾取。計(jì)算航電各頻道實(shí)、虛分量曲線的梯度,根據(jù)正、負(fù)峰組合判斷拾取波峰、波谷,一個(gè)波峰及相鄰兩個(gè)波谷構(gòu)成一個(gè)被選異常。
4)異常信息提取。計(jì)算被選異常的峰值、半極值寬度、視電阻率、視電阻率背景、測(cè)量高度、多頻道交集、視電阻率峰值、視飛行高度等信息,并根據(jù)半極值寬度等條件識(shí)別雙峰異常。
圖2 已知礦化點(diǎn)異常和殘坡積引起的覆蓋層異常綜合剖面圖Fig.2 Composite profiles of an ore anomaly and a residual sediment anomaly (a)礦化點(diǎn)異常;(b)覆蓋層異常
5)應(yīng)用約束條件篩選異常。應(yīng)用異常篩選條件,篩選出符合各約束條件交集的異常。如上所述,約束條件可根據(jù)任務(wù)需要進(jìn)行設(shè)定,因而可以區(qū)分出50Hz工業(yè)電干擾異常、覆蓋層異常等特定類型異常。
6)輸出異常篩選結(jié)果。將異常結(jié)果輸出為文本文件,其中包含異常的各種信息,及位置坐標(biāo)、異常種類、評(píng)級(jí)等,供方法解釋人員進(jìn)一步篩選。
7)輔助信息解釋。根據(jù)地形剖面曲線等輔助條件對(duì)異常篩選結(jié)果進(jìn)行進(jìn)一步分類,剔除負(fù)地形異常及地形效應(yīng)異常。
將異常識(shí)別算法用軟件代碼實(shí)現(xiàn),然后再依據(jù)以上流程就可以編織出航電異常識(shí)別程序(圖3)。
圖3 異常自動(dòng)識(shí)別流程Fig.3 Flow chart of the digital anomaly detection and classification
圖4為內(nèi)蒙古某測(cè)區(qū)航電異常自動(dòng)拾取的結(jié)果。圖中彩色底圖為視電阻率平面圖(暖色代表視低電阻率,冷色代表高視電阻率),紫色三角點(diǎn)為自動(dòng)拾取的異常位置,黑色箭頭所指藍(lán)色圓點(diǎn)為已知礦或礦化點(diǎn)??梢钥吹皆趫D中間最上部①、④號(hào)已知礦化點(diǎn)與自動(dòng)異常拾取的位置是完全對(duì)應(yīng)的。①號(hào)異常為已知的熱液型銅礦化點(diǎn),自動(dòng)識(shí)別在8 040線上有明顯航電響應(yīng)異常與之對(duì)應(yīng)(圖5)。該異常控制長度2 km,半極值寬度297 m,飛行高度48 m,中頻和低頻實(shí)虛分量比分別為1.88和1.76,異常清晰。④號(hào)已知熱液型銅礦化點(diǎn)在8 170、8 180線間穿過,在8 180線礦化點(diǎn)附近自動(dòng)識(shí)別找到了航電異常,但異常峰值較小。②號(hào)已知熱液型銅礦化點(diǎn)在6 870線上方通過,但測(cè)線上未見明顯異常,在相鄰的6 850線、6 860線航電異常明顯。③號(hào)已知熱液型錳礦化點(diǎn)位于8 240線、8 250線間,為熱液型銅礦化點(diǎn),在測(cè)線上未發(fā)現(xiàn)明顯異常,但在相鄰的8 270、8 280、8 290線上自動(dòng)識(shí)別出了異常條帶。
圖4 內(nèi)蒙某測(cè)區(qū)航電異常自動(dòng)識(shí)別篩選結(jié)果Fig.4 Digital detection and classification results of a survey area in Inner Mongolia
圖5 內(nèi)蒙某測(cè)區(qū)8040線上的已知礦上航電異常響應(yīng)
頻率域航空電磁異常的數(shù)字化自動(dòng)識(shí)別篩選技術(shù)補(bǔ)充完善了國產(chǎn)HDY-402三頻航空電磁系統(tǒng)數(shù)據(jù)處理技術(shù),開發(fā)的異常識(shí)別程序?qū)嵱?、高效,可根?jù)解釋人員定義的不同約束條件檢出不同類型的航電異常,提高了數(shù)據(jù)處理的效率和質(zhì)量,填補(bǔ)了國內(nèi)該領(lǐng)域的技術(shù)空白。但是受到數(shù)據(jù)質(zhì)量不佳及數(shù)據(jù)信息量有限等條件的限制,航電異常的自動(dòng)識(shí)別篩選結(jié)果中包含有干擾異常,僅能實(shí)現(xiàn)異常的初步篩選,還不能完全替代解釋人員的人工判別。
參考文獻(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] 周鳳桐. 補(bǔ)償式航空電磁法[J]. 物化探研究報(bào)道,1982(11):1-78.
[8] 李文杰. 頻率域航空電磁數(shù)據(jù)處理技術(shù)研究 [D] .北京:中國地質(zhì)大學(xué), 2008.