馬氏距離法在東昆侖東段多元異常圈定中的對比試驗

耿國帥，楊帆

(1.中國地質(zhì)大學(xué)(北京) 地球科學(xué)與資源學(xué)院，北京 100083； 2.中國地質(zhì)調(diào)查局 地球物理調(diào)查中心，河北 廊坊 065000； 3.北京礦產(chǎn)地質(zhì)研究院，北京 100012； 4.中國地質(zhì)調(diào)查局 土地質(zhì)量地球化學(xué)調(diào)查評價研究中心，河北 廊坊 065000)

0 引言

在地球化學(xué)數(shù)據(jù)處理中，許多情況下異常下限并不是某個重要元素的濃度值, 而是樣品中幾種元素的相互關(guān)系[1]，因此，用傳統(tǒng)的統(tǒng)計學(xué)方法確定的單元素異常下限, 可能無法揭示或突出數(shù)據(jù)集中有價值的復(fù)合信息, 還可能誤入歧途。目前，已有多種方法來識別多元素地球化學(xué)異常[2-5]，但應(yīng)用最廣泛的是基于馬氏距離的多元異常識別方法[6-13]。常規(guī)馬氏距離方法是依據(jù)數(shù)據(jù)的算術(shù)平均值和樣本的協(xié)方差矩陣來計算各樣本點離數(shù)據(jù)中心的距離，將計算的馬氏距離與卡方分布的臨界值(通常取卡方分布的97.5%)相比較，大于臨界值的數(shù)據(jù)點為異常點[14]。但算術(shù)平均值和協(xié)方差矩陣是不穩(wěn)健的統(tǒng)計量，部分異常點的存在會造成算術(shù)平均值的偏大或偏小及協(xié)方差矩陣的變化，進而掩蓋一些異常點[15-19]。為了克服不穩(wěn)健的問題，許多科學(xué)家采用了穩(wěn)健的方法來計算馬氏距離[20-21]。應(yīng)用最廣泛的是Rousseeuw提出的基于最小協(xié)方差矩陣行列式(MCD)的快速方法(FMCD)[22]。Majewska認(rèn)為，基于M-估計、S-估計、MM-估計、MVE、MCD 和FMCD的穩(wěn)健馬氏距離在識別異常的效果上已經(jīng)證明比常規(guī)馬氏距離好，在這些穩(wěn)健方法中，F(xiàn)MCD的方法相對較好[23]。Filzmoser等認(rèn)為這種根據(jù)馬氏距離大于某一固定值判定為異常的方法也有很大的缺陷，它沒有考慮樣品的個數(shù)及數(shù)據(jù)的結(jié)構(gòu)，因此提出了基于校正的穩(wěn)健馬氏距離法(adaptive method)[24]。此外Sajesh和Srinivasan提出了用Comedian來代替Covariance的穩(wěn)健馬氏距離法來識別異常(comedian method)[25]。筆者利用常規(guī)馬氏距離方法、FMCD法、Adaptive法和Comedian法等4種方法，分別對東昆侖東段水系沉積物測量地球化學(xué)數(shù)據(jù)進行處理，用于圈定多元素組合異常效果對比，旨在優(yōu)選適宜該區(qū)的方法。

1 研究區(qū)概況

青海省東昆侖東段跨越了古亞洲和特提斯兩大構(gòu)造域，位于南北兩大構(gòu)造域結(jié)合部位，區(qū)內(nèi)可劃分為秦祁昆、特提斯2個一級單元，東昆侖造山帶、巴顏喀拉造山帶2個二級構(gòu)造單元及昆北、昆中、昆南、北巴 4 個三級構(gòu)造單元(圖1)。

區(qū)內(nèi)地層出露齊全，從古元古界到新生界均有出露，其中昆北帶主要出露三疊系鄂拉山組，昆中帶主要出露太古宇古老變質(zhì)巖系，昆南帶出露地層最為復(fù)雜，從古元古界到侏羅系均有出露，北巴帶內(nèi)主要出露三疊系巴顏喀拉群的復(fù)理石沉積。區(qū)內(nèi)巖漿侵入活動強烈，以花崗巖類為主，在4個構(gòu)造單元均有出露，尤以昆中帶最密集，北巴帶最稀少，時代從元古宙、早古生代、晚古生代到中生代均有分布，其中中生代最為發(fā)育。

該區(qū)已發(fā)現(xiàn)金、鐵、銅、鉛、鋅、鎢、錫、鈷等礦床(點)110多處，其中大型礦床4處，都為金礦床；中型礦床7處，為金、鐵、多金屬礦；小型礦床13處，為金、銅、鐵、多金屬礦；礦點46處；礦化點41處。近些年來，該區(qū)一直是找礦工作的熱點地區(qū)之一。

張德全把該區(qū)的礦床類型分為噴氣—沉積礦床組合(VHMS型和SEDEX型)和造山礦床組合(斑巖、矽卡巖、熱液脈型和造山型金礦)共兩個組合、六個類型[26]；丁清峰把該區(qū)的礦床類型分為熱水噴流沉積和與造山作用有關(guān)的斑巖型、矽卡巖型、熱液脈型和疊加改造型[27]；田立明把該區(qū)主要成礦類型歸為海相火山巖型、斑巖型/矽卡巖型和造山型金礦等3種組合[28]。筆者根據(jù)前人的研究成果，把該區(qū)的金屬礦床分為3種成因組合(表1)。

1—主縫合帶；2—次縫合帶；3—新元古代-早古生代結(jié)合帶俯沖方向(一側(cè)有齒者為單向俯沖，兩側(cè)有齒者為雙向俯沖)；4—晚古生代-早中生代縫合帶俯沖方向；5—A型俯沖帶；6—公路；7—研究區(qū)位置；Ⅰ—柴達木地塊；Ⅱ—東昆侖造山帶；Ⅱ1—東昆北早古生代弧后裂陷帶(昆北帶)；Ⅱ2—東昆中巖漿弧帶(昆中帶)；Ⅱ3—東昆南構(gòu)造-混雜巖帶(昆南帶)；Ⅲ—巴顏喀拉造山帶(北巴帶)1—main structure zone; 2—secondary structure zone; 3—Neoproterozoic-early Paleozoic combined belt subduction direction(one-way subduction with teeth on one side and two-way subduction with teeth on both sides); 4—subduction direction of late Paleozoic and early Mesozoic suture belt; 5—A type subduction zones; 6—high way; 7—location of study area; Ⅰ—Qaidam massif; Ⅱ—East Kunlun orogenic belt; Ⅱ1—East Kunbei early paleozoic back-arc rife (Kunbei belt); Ⅱ2—East Kunzhong magmatic arc zone(Kunzhong belt); Ⅱ3—East-Kunnan tectonomagmatic belt(Kunnan belt); Ⅲ—Bayan Kara orogenic belt (Beiba belt)圖1 東昆侖地區(qū)大地構(gòu)造分區(qū)Fig.1 The map of geotectonic units in the study area

表1 研究區(qū)礦床成因類型

2 四種馬氏距離算法和元素組合的選取

2.1 四種馬氏距離算法

馬氏距離是由印度統(tǒng)計學(xué)家P. C.Mahalanobis提出的，表示數(shù)據(jù)的協(xié)方差距離[29]。它是一種廣義距離，因為考慮了各變量間的協(xié)方差，與普通歐氏距離相比,能消除量綱及各變量間相關(guān)性的影響。馬氏距離的計算方法為

式中：Xi為分析數(shù)據(jù)；T為數(shù)據(jù)中心位置估計值；C為數(shù)據(jù)離散的估計值。

1) 在常規(guī)方法中，T為算術(shù)平均值，C為協(xié)方差矩陣。

2) FMCD法，實際上就是挑選一個不含異常點的背景子樣本，使其矩陣行列式的值最小，從而計算出該子樣本的算術(shù)平均值和協(xié)方差矩陣，來代替上式中的T和C。子樣本個數(shù)須適宜，少了代表性差，多了穩(wěn)健性差，樣本個數(shù)一般在(n+p+1)/2～3n/4，其中n為樣本數(shù)，p為元素個數(shù)[22]。

3) Adaptive法，計算馬氏距離的算法同步驟2)，只是判斷異常下限的方法不同，它考慮了樣本容量和元素個數(shù)對異常下限的影響[24]。

4) Comedian法，是利用Comedian來代替Covariance的方法，但直接計算的Comedian是一個非半正定矩陣，因此采用了奇異值分解的方法來重新生成一個半正定的Comedian矩陣和計算出中心位置的估計值[25]。

2.2 元素組合的選擇

馬氏距離圈定多元素異常，通常是針對不同類型礦床的特征元素組合的。本文根據(jù)該區(qū)的礦床(點)個數(shù)、規(guī)模及與元素間的關(guān)系，挑選出下列3種元素組合：

1) 以基性巖有關(guān)的礦床類型的元素組合：Cu、Co、Cr、Ni、V、Fe；

2) 以酸性巖有關(guān)的礦床類型的元素組合：Cu、Mo、Pb、Zn、Ag、Cd；

3) 以造山型金礦有關(guān)的元素組合：Au、As、Sb。

3 結(jié)果及討論

3.1 數(shù)據(jù)的來源

本文所采用的數(shù)據(jù)來源于青海省地質(zhì)礦產(chǎn)局 1∶50萬水系沉積物地球化學(xué)掃面工作(1個點/4 km2，8～16 km2的組合樣作為一個分析樣)。這項工作為近年來東昆侖地區(qū)的金多金屬資源調(diào)查提供了極為寶貴的基礎(chǔ)資料。張文秦依據(jù)該數(shù)據(jù)，對東昆侖地質(zhì)背景進行探討[30]；安國英對數(shù)據(jù)進行處理，分析了元素的地球化學(xué)特征，總結(jié)了各構(gòu)造單元地球化學(xué)元素的分布特征，并用于金礦靶區(qū)的篩選及評價中[31-32]。筆者旨在利用上述數(shù)據(jù)對比優(yōu)選適宜該區(qū)的多元素組合異常圈定的方法。

3.2 不同方法計算的馬氏距離對比

對研究區(qū)4 001件樣品按照上述3種元素組合計算其馬氏距離，由于Adaptive馬氏距離是利用FMCD的方法計算得出的，因此4種方法僅有3種馬氏距離。把3種元素組合經(jīng)FMCD、基于Comedian計算的馬氏距離分別與常規(guī)方法計算的馬氏距離進行對比(圖2)。從圖2可看出：① 常規(guī)馬氏距離與基于FMCD計算的馬氏距離相關(guān)性較好，如Au、As、Sb中，兩者的相關(guān)性近乎呈一條直線。而基于Comedian計算的馬氏距離與常規(guī)馬氏距離的相關(guān)性較差。② 幾乎所有基于FMCD計算的穩(wěn)健馬氏距離值都大于常規(guī)方法計算的馬氏距離值，而基于Comedian計算的穩(wěn)健馬氏距離值大部分大于常規(guī)方法計算的距離值，少部分小于常規(guī)方法計算的距離值。③ 常規(guī)方法圈定的異常點一定是基于FMCD和Adaptive方法圈定的異常點(通過計算，Adaptive法確定的異常下限都大于FMCD法確定的異常下限，因此FMCD穩(wěn)健距離確定的異常點一定包括Adaptive穩(wěn)健馬氏距離確定的異常點)，但不一定是基于Comedian穩(wěn)健距離圈定的異常點。④基于FMCD計算的馬氏距離最大值(圖2中紅色圓圈所突出的點)與基于Comedian計算的馬氏距離最大值(圖2中藍(lán)色圓圈所突出的點)通常不一致，且相差較大，例如Cd、Pb、Zn、Cu、Mo、Ag元素組合中，用FMCD計算的距離最大值，在Comedian計算的距離中甚至達不到異常點的標(biāo)準(zhǔn)，而用Comedian計算出的最大值點，在FMCD計算的距離中值都較大，在FMCD和Adaptive中都是異常點。

3.3 異常下限及識別的異常點數(shù)對比

對比4種馬氏距離方法確定3種元素組合的異常下限及異常點數(shù)，結(jié)果顯示，常規(guī)方法識別異常能力較弱，識別出的異常點數(shù)較少，而其他3種穩(wěn)健馬氏距離識別異常的能力較強，圈定的異常點數(shù)較多，克服了異常數(shù)據(jù)對不穩(wěn)健馬氏距離的的掩蔽效應(yīng)。但各種穩(wěn)健方法識別出的異常點數(shù)也大不相同，Comedian識別出的Au、As、Sb組合異常點最多，而識別出的Cu、Co、Cr、Ni、V、Fe組合異常的點數(shù)最少；FMCD識別出3種元素組合異常的個數(shù)差別不大，而Adaptive識別出的異常點數(shù)差別介于兩者之間(表2)。為了探究元素組合異常與元素離散度間的關(guān)系，對3種組合中的元素標(biāo)準(zhǔn)化后作箱線圖(圖3)，標(biāo)準(zhǔn)化的方法如下：

圖2 3種元素組合的常規(guī)(a)和穩(wěn)健(b)馬氏距離對比Fig.2 The comparison of classical and robust mahalanobis distance from three element associations

表2 4種馬氏距離確定的異常下限及異常點數(shù)統(tǒng)計

Z=(X-median(X))/MAD(X),

MAD(X)=1.4826×median|X-median(X)|。

式中：median表示數(shù)據(jù)集的中值，而MAD表示中值的絕對偏差，目的是統(tǒng)一各元素的量綱，使各元素值在統(tǒng)一水平線上，同時使各元素的中值處于箱線圖縱軸的零刻度位置，方便進行比較。

標(biāo)準(zhǔn)化元素箱線圖顯示，該區(qū)Sb、Au、As元素組合中，元素的離散度都較大，高值點較多，形成的異常也應(yīng)該較多；Cd、Cu、Mo、Pb、Zn、Ag元素組合中，Pb、Mo、Cd、Ag離散度也較大，Cu、Zn離散度較??；而Co、Cr、Cu、Ni、V、Fe元素組合中，除Cr、Ni離散度較高外，其他元素離散度都較小。

圖3 元素標(biāo)準(zhǔn)化后的箱線Fig.3 Boxplot for elements of standardized data

從該區(qū)所形成的礦床來看，造山型金礦床無疑是該區(qū)最具找礦潛力的類型，與中酸性巖有關(guān)的矽卡巖型和斑巖型次之，與基性巖有關(guān)的礦床找礦潛力最小。故基于Comedian的方法較其他兩種穩(wěn)健方法更好。

3.4 圈定的異常對比

根據(jù)4種方法確定的異常下限值，圈定了Cu、Co、Cr、Ni、V、Fe，Cd、Cu、Pb、Zn、Mo、Ag和Au、As、Sb組合異常(圖4～圖6)，各異常特征對比如下。

3.4.1 Cu、Co、Cr、Ni、V、Fe組合異常

使用常規(guī)馬氏距離方法圈定的異常面積較小，異常較零星，強度較低，濃集趨勢不明顯，而基于穩(wěn)健馬氏距離圈定的異常，強度較高，濃集趨勢也比較明顯。從3種穩(wěn)健方法所圈定的異?？矗珻omedian圈定的異常主要集中分布在昆南帶內(nèi)，而FMCD和Adaptive所圈定的異常，除了大部分集中在昆南帶內(nèi)外，在昆北帶內(nèi)也有較多的小異常存在，在昆中和北巴帶內(nèi)的異常也較前者多。常規(guī)方法和FMCD計算的馬氏距離最大值點都分布在北巴帶內(nèi)，而Comedian計算的馬氏距離最大值點分布在昆南帶內(nèi)的清水泉鉻鐵礦點附近。相較于FMCD方法，Adaptive方法所圈定的異常面積減少，但濃集趨勢不變，礦床與異常對應(yīng)關(guān)系也沒有改變，表明Adaptive圈定的異常較FMCD合理?？傊珻omedian法更符合地質(zhì)情況，所圈定的異常最真實，其次是Adaptive，常規(guī)方法效果最差。

3.4.2 Cd、Cu、Mo、Pb、Zn、Ag組合異常

使用常規(guī)馬氏距離圈定的異常面積較小，分布較零星，強度較低，濃集趨勢不明顯，而采用穩(wěn)健馬氏距離圈定的異常面積較大，強度較高，濃集趨勢較明顯，與礦床點的對應(yīng)性較好。從3種穩(wěn)健方法所圈定的異?？?，Comedian所圈定的異常強度更高，濃集趨勢更明顯，與礦床點的對應(yīng)性也更好。從馬氏距離最大值所處位置來看，常規(guī)方法和Comedian穩(wěn)健方法計算的最大值點都處在圖幅中間，山根果勒南端，而FMCD所計算的最大值點處于圖幅最東端哈拉深溝的南東。目前在Comedian極值點附近已發(fā)現(xiàn)了沒桑確昂、注斯楞等熱液脈型鉛礦點。相較于FMCD，Adaptive圈定的異常面積減小，但濃集趨勢及與礦床點的關(guān)系沒有改變。

3.4.3 Au、As、Sb組合異常

使用常規(guī)方法圈定的Au、As、Sb組合異常，較常規(guī)方法所圈定的前兩種元素組合異常，分布更零星，面積更小，強度更低，也說明了常規(guī)方法不穩(wěn)健所帶來的掩蔽效果對Au、As、Sb組合更嚴(yán)重。而從3種穩(wěn)健方法所圈定的異常來看，用Comedian所圈定的異常強度更高，濃集趨勢更明顯，且從馬氏距離最大值所處的位置來看，Comedian所計算的馬氏距離最大值處于五龍溝礦床，而其他兩種方法計算的馬氏距離都處于研究區(qū)的最南端。

3.5 Adaptive和Comedian識別出的異常點分布

對Adaptive和Comedian所識別出的異常點分別進行統(tǒng)計(見表3)，從所識別的異常點數(shù)看，兩種方法識別的Co、Cr、Cu、Ni、V、Fe組合異常點數(shù)區(qū)別較小，共同識別的異常點數(shù)占總異常點數(shù)近50%，僅被一種方法所識別的異常點數(shù)各占總異常點數(shù)的近25%。而在識別其他兩種元素組合的異常點數(shù)時差別較大，在Cd、Cu、Mo、Pb、Zn、Ag元素組合中，兩者共同識別出的異常點數(shù)僅占1/3左右，而Comedian所識別的異常點數(shù)遠(yuǎn)高于Adaptive所識別的異常點數(shù)；在Au、As、Sb元素組合中，兩者共同識別出的異常點數(shù)大于總異常點數(shù)的40%，Comedian所識別的異常點數(shù)也遠(yuǎn)遠(yuǎn)大于Adaptive所識別的異常點數(shù)。

異常識別方法的優(yōu)劣，取決于異常點個數(shù)及其與地質(zhì)及礦床點分布的吻合程度。因此，筆者選擇兩種方法所識別異常點數(shù)區(qū)別較小的Co、Cr、Cu、Ni、V、Fe組合，比較兩種方法所識別的異常點與礦床點的吻合程度(見圖7)。

從該區(qū)所收集到的25處該類型的礦床點與異常點的分布情況看，除了1處礦床點周圍無異常點外，其他24處礦床點都處在異常點上或異常點附近，且大都處于兩者共同識別出的異常點上，但從礦床點與兩者單獨識別出的異常點的關(guān)系看，9處礦床點處于Comdeian識別出的異常點附近，而只有3處礦床點處于Adaptive識別出的異常點附近，故Comedian識別出的異常點更合理。

表3 Adaptive和Comedian異常點統(tǒng)計

圖5 4種馬氏距離圈定的Cd、Cu、Mo、Pb、Zn、Ag組合異常Fig.5 Multivariate outlier delineation for Cd，Cu，Mo,Pb,Zn and Ag from four mahalanobis diatance methods in study area

圖6 4種馬氏距離圈定的Au,As,Sb組合異常Fig.6 Multivariate outlier delineation for Au,As,Sb from mahalanobis distance methods in study area

圖7 Adaptive和Comedian識別的Co、Cr、Cu、Ni、V、Fe組合異常點分布Fig.7 Outliers distribution of Co,Cr,Cu,Ni,V,Fe association using Adaptive and Comedian methods

4 結(jié)論

常規(guī)馬氏距離法、FMCD法、Adaptive法、Comedian法在東昆侖東段異常識別和圈定中的實際效果顯示：

1) 常規(guī)馬氏距離容易受到特異值的影響，而造成異常的掩蔽效應(yīng)，識別出的異常點數(shù)較少，圈定的異常面積小，分布零星，強度低。

2) 3種穩(wěn)健方法由于消除了特異值的掩蔽效應(yīng)，識別出的異常點數(shù)多，圈定的異常面積大，強度高，濃集趨勢明顯。但由于FMCD沒有考慮數(shù)據(jù)結(jié)構(gòu)，圈定的異常只與變量個數(shù)有關(guān)，因此與實際情況不太相符；而Adaptive雖然考慮了數(shù)據(jù)結(jié)構(gòu)，但馬氏距離的計算是基于FMCD的方法，Maronna等指出，F(xiàn)MCD估計值的偏差會隨著變量數(shù)的增加而增大[32-33]； Comedian是基于中值和中值的絕對偏差，因此具有高穩(wěn)健性和高截斷點，但由于Comedian是非半正定矩陣，因此不能直接計算，而要采用奇異值分解的方法來重構(gòu)Comdian矩陣。

3) 綜上所述，4種方法在東昆侖東段多元素異常圈定效果依次為Comedian方法、Adaptive法、FMCD法、常規(guī)方法。