孫堯堯，郝立波，趙新運(yùn)，陸繼龍，馬成有，魏俏巧

(吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院，吉林 長(zhǎng)春 130026)

0 引言

在地球化學(xué)找礦工作中，異常的圈定是極為重要的環(huán)節(jié)。目前，在巖性復(fù)雜區(qū)水系沉積物地球化學(xué)異常圈定時(shí)仍存在一些亟待解決的問(wèn)題，其中，巖性背景問(wèn)題最為突出，尤其是在區(qū)域地球化學(xué)勘查中[1-4]。若不消除巖性背景差異的影響而直接采用統(tǒng)一的背景上限進(jìn)行異常圈定，一是可能會(huì)將高背景值區(qū)的一些正常樣品劃分為異常樣品，從而形成虛假地球化學(xué)異常；二是可能會(huì)將低背景值區(qū)的一些異常樣品劃分為正常樣品，從而掩蓋掉部分低弱但重要的地球化學(xué)異常[3-5]。

針對(duì)巖性背景問(wèn)題，學(xué)者們提出了許多處理方法，如多元線性回歸法[6]、滑動(dòng)平均標(biāo)準(zhǔn)化法[1]和一元線性回歸法[3-4]等。然而，這些方法的提出并不是從直接區(qū)分巖性背景的角度出發(fā)的，或多或少存在一些局限性，因此很難從根本上消除巖性背景的影響。Zhao等認(rèn)為巖性背景問(wèn)題的本質(zhì)是多重母體問(wèn)題，因此提出利用Expectation-Maximization(EM)算法來(lái)分解多重母體，進(jìn)而達(dá)到消除巖性背景影響的目的[5]。該方法地質(zhì)意義明確，對(duì)合理圈定巖性復(fù)雜區(qū)水系沉積物地球化學(xué)異常應(yīng)能起到一定的作用。筆者以湖南省某地1∶20萬(wàn)水系沉積物地球化學(xué)數(shù)據(jù)為例，嘗試將EM聚類方法應(yīng)用到地球化學(xué)背景與異常劃分中，旨在為巖性復(fù)雜區(qū)水系沉積物地球化學(xué)異常的圈定提供參考。

1 巖性背景問(wèn)題

研究表明，地球化學(xué)背景與異常劃分的影響因素主要有巖性背景、地球化學(xué)景觀條件以及樣品粒度等[1-2,6]。但在相同或相近的景觀條件下，巖性背景的影響最為突出[3-4]。巖性復(fù)雜區(qū)往往存在多種巖石類型，不同類型巖石的元素豐度通常存在明顯差異。研究表明，巖性單一的地質(zhì)體的元素含量往往服從(近似)正態(tài)分布[7-8]，其數(shù)據(jù)可構(gòu)成一(近似)正態(tài)分布的母體。因此，巖性復(fù)雜區(qū)應(yīng)存在多個(gè)(近似)正態(tài)母體。巖石風(fēng)化形成的水系沉積物對(duì)基巖的化學(xué)成分具有顯著的繼承性[9-10]，因此巖性復(fù)雜區(qū)水系沉積物元素含量數(shù)據(jù)內(nèi)也應(yīng)包含多個(gè)(近似)正態(tài)母體[5]。來(lái)源于不同巖性背景的元素含量數(shù)據(jù)發(fā)生混合，即正態(tài)母體混合或正態(tài)分布疊加，可產(chǎn)生偏態(tài)分布，包括左偏、右偏及多峰等分布型式[7-8]，這是巖性背景對(duì)地球化學(xué)異常圈定的直接影響。由此可見，巖性背景問(wèn)題即多重母體問(wèn)題，而多重母體問(wèn)題主要源于巖性背景差異。

針對(duì)上述問(wèn)題，學(xué)者們開展了大量的研究工作，提出了許多處理方法。研究表明，在使用任何數(shù)據(jù)處理方法劃分地球化學(xué)背景與異常之前，應(yīng)先確定元素含量的概率分布型式，且以(近似)正態(tài)分布為目標(biāo)[11]。為滿足正態(tài)分布的前提，一些數(shù)據(jù)變換方法被應(yīng)用到異常圈定中，如對(duì)數(shù)變換法[12-13]、倒數(shù)變換法[14]和Box-Cox變換法[15]等，其中對(duì)數(shù)變換法最為常用。Stanley研究指出，對(duì)數(shù)變換及類似的數(shù)學(xué)變換方法通常難以消除地球化學(xué)數(shù)據(jù)偏態(tài)分布的特征[16]。在地球化學(xué)數(shù)據(jù)處理工作中，這一現(xiàn)象是普遍存在的，通常是偏態(tài)分布經(jīng)對(duì)數(shù)轉(zhuǎn)換后仍不服從(近似)正態(tài)分布。然而，這些方法目前仍被廣泛使用，可能已經(jīng)造成了嚴(yán)重的錯(cuò)誤。

按地質(zhì)單元界線區(qū)分巖性背景的效果也不是很理想，原因是同一地質(zhì)單元的巖性組成往往并不單一，即同一地質(zhì)單元內(nèi)也可能含有多個(gè)(近似)正態(tài)母體。因此，一些學(xué)者嘗試從回歸分析的角度來(lái)解決巖性背景問(wèn)題。周蒂提出了多元線性回歸分析方法[6]，郝立波等提出了一元線性回歸分析方法[3]。相比于傳統(tǒng)的數(shù)據(jù)變換處理方法，這些方法具有很大的優(yōu)勢(shì)，可以在一定程度上消除巖性背景的影響，較準(zhǔn)確地劃分出成礦元素的地球化學(xué)背景與異常。然而，回歸分析方法也有其局限性，因?yàn)樵擃惙椒▋H適用于那些與巖性存在良好線性關(guān)系的成礦元素，而不適用于那些受巖性背景影響小或不受影響的元素。另外，線性回歸分析要求殘差符合正態(tài)分布，這也是該方法使用受限的原因之一。

李寶強(qiáng)和孫澤坤提出了滑動(dòng)平均標(biāo)準(zhǔn)化方法[1]，其不足之處是窗口大小的確定具有很大的人為性，若窗口過(guò)大，則難以消除元素背景含量的影響，而窗口過(guò)小，又會(huì)縮小各點(diǎn)的元素含量差異，損失異常信息。至于分形/多重分形方法，對(duì)于服從分形分布的元素含量數(shù)據(jù)，該類方法可能是一種很好的選擇[17-18]，然而，該類方法的地質(zhì)和地球化學(xué)意義仍不明確[3-5,19]，另外，該類方法也無(wú)法區(qū)分來(lái)自不同巖性母體的樣品。同樣地，線性判別分析方法[8,20]和概率格紙法[21]等也面臨同樣的問(wèn)題，即無(wú)法判別樣品具體歸屬于哪一類巖性母體[5]。

解決巖性背景問(wèn)題的根本途徑應(yīng)是直接區(qū)分巖性背景或分解多重母體。Zhao等研究認(rèn)為，在巖性復(fù)雜區(qū)可利用能反映巖性變化的主量和(或)親石微量元素，通過(guò)EM聚類算法——一種有效的高斯混合模型聚類方法，對(duì)來(lái)自不同巖性背景的水系沉積物樣品進(jìn)行分類，以達(dá)到區(qū)分巖性背景或分解多重母體的目的[5]?；诖?，筆者以湖南省某地1∶20萬(wàn)水系沉積物地球化學(xué)數(shù)據(jù)為例，嘗試采用EM聚類方法進(jìn)行多重母體分解，然后進(jìn)行地球化學(xué)異常圈定，并將其與傳統(tǒng)的均值+k倍標(biāo)準(zhǔn)偏差方法圈定的異常進(jìn)行比較。

2 方法流程

2.1 樣品分類

樣品分類主要基于EM算法。該算法由Dempster等提出[22]，是一種重要的統(tǒng)計(jì)學(xué)分析工具，已被廣泛應(yīng)用于含有隱變量的概率模型參數(shù)的極大似然估計(jì)中。假設(shè)在某一復(fù)雜巖性區(qū)采集了n個(gè)水系沉積物樣品，其分析結(jié)果可設(shè)為隨機(jī)變量X，X=(x1,x2,x3,…,xn)，樣品采自不同的巖性背景區(qū)，即X內(nèi)含有多個(gè)正態(tài)母體，但每個(gè)樣品的對(duì)應(yīng)母體或類別是未知的。正態(tài)母體的分布參數(shù)也是未知的，樣品的對(duì)應(yīng)類別可設(shè)為變量Z，Z=(z1,z2,z3,…,zn)。EM算法是先初始化正態(tài)母體的分布參數(shù)θ，然后重復(fù)以下兩個(gè)步驟：①根據(jù)分布參數(shù)θ計(jì)算樣品分類變量Z的期望，Qi(zi)=p(zi|xi;θ)；②將概率模型參數(shù)的似然函數(shù)L(θ)最大化以獲得新的分布參數(shù)θ，即

似然函數(shù)L(θ)收斂時(shí)，迭代停止。最終，可得到樣品的具體分類情況，即在尋找最優(yōu)參數(shù)來(lái)極大化樣本似然函數(shù)的過(guò)程中對(duì)樣本完成最優(yōu)分類[5]。

在實(shí)際問(wèn)題處理中，應(yīng)首先根據(jù)研究區(qū)地質(zhì)概況選取可以反映巖性變化，且元素間相關(guān)性較小的非成礦常量和(或)親石微量元素作為聚類指標(biāo)，然后采用EM算法對(duì)樣品進(jìn)行分類。最優(yōu)分類數(shù)可采用數(shù)學(xué)衡量指標(biāo)和地質(zhì)條件雙重限制進(jìn)行合理判別[5]。數(shù)學(xué)衡量指標(biāo)為赤池信息量準(zhǔn)則AIC(akaike’sinformationcriterion)，一個(gè)衡量模型擬合程度的標(biāo)準(zhǔn)[23]。AIC可以表示為：

AIC=2k-2ln(L)。

式中：k是模型的獨(dú)立參數(shù)個(gè)數(shù)；L是模型的極大似然函數(shù)。按極大似然函數(shù)準(zhǔn)則，L越大，模型的擬合度越高，但實(shí)際上，在L過(guò)大的情況下，可能會(huì)出現(xiàn)過(guò)擬合現(xiàn)象。為了避免過(guò)度擬合，AIC準(zhǔn)則在追求L盡可能大的同時(shí)，還要求模型參數(shù)要盡可能少，即k要盡可能小，換句話說(shuō)，AIC值越小，模型的擬合情況越好。因此，在有限混合正態(tài)分布參數(shù)的極大似然估計(jì)中，當(dāng)AIC值達(dá)到最小時(shí)可獲得最優(yōu)分類數(shù)。采用地質(zhì)條件限制時(shí)，可先根據(jù)分類結(jié)果繪制分類圖，然后將其與地質(zhì)圖進(jìn)行對(duì)比，與地質(zhì)圖吻合度最高的分類圖對(duì)應(yīng)的分類數(shù)可確定為最優(yōu)分類數(shù)。按統(tǒng)計(jì)學(xué)的要求，每個(gè)分類中樣品數(shù)應(yīng)不少于30個(gè)，若不滿足，則需調(diào)整分類，或?qū)悠窋?shù)較少的類合并到與其最近的類中。

2.2 數(shù)據(jù)標(biāo)準(zhǔn)化

樣品分類完成后，對(duì)每個(gè)分類中的樣品成礦元素含量數(shù)據(jù)進(jìn)行3S檢驗(yàn)，將均值±3倍標(biāo)準(zhǔn)偏差范圍外的異常數(shù)據(jù)剔除。然后，利用剩余數(shù)據(jù)計(jì)算每個(gè)分類中成礦元素含量的均值和標(biāo)準(zhǔn)偏差，以此對(duì)成礦元素的含量進(jìn)行Z分?jǐn)?shù)標(biāo)準(zhǔn)化，公式如下：

i=1,2,3,…,n;j=1,2,3,…,m。

2.3 異常圈定

經(jīng)樣品分類和元素含量標(biāo)準(zhǔn)化后，巖性背景差異在一定程度上得到了消除。因此，可在全區(qū)采用統(tǒng)一的背景上限進(jìn)行異常圈定，即直接采用傳統(tǒng)的均值加k倍標(biāo)準(zhǔn)偏差的方法或其他方法進(jìn)行異常圈定，k值的大小可依據(jù)實(shí)際情況進(jìn)行適當(dāng)調(diào)整。

3 應(yīng)用實(shí)例

3.1 研究區(qū)概況

以湖南省某地1∶20萬(wàn)水系沉積物地球化學(xué)數(shù)據(jù)為例，進(jìn)行了基于多重母體分解的地球化學(xué)異常圈定試驗(yàn)。研究區(qū)面積約為4 300 km2，水系沉積物樣品數(shù)為1 050，采樣密度約為1個(gè)樣品/4 km2。水系沉積物化學(xué)成分分析項(xiàng)目包括Ag、As、Au、B、Ba、Be、Bi、Cd、Co、Cr、Cu、F、Hg、La、Li、Mn、Mo、Nb、Ni、P、Pb、Sb、Sn、Sr、Th、Ti、U、V、W、Y、Zn、Zr、Al2O3、CaO、Fe2O3、K2O、MgO、Na2O及SiO2。研究區(qū)發(fā)育有多個(gè)鎢錫多金屬礦床，成礦地質(zhì)條件優(yōu)越。研究區(qū)地層復(fù)雜，但巖性較為簡(jiǎn)單，按照出露面積從大到小的排列依次為上古生界灰?guī)r、泥灰?guī)r，下古生界淺變質(zhì)砂巖，白堊系砂巖，泥盆系砂巖和震旦系砂質(zhì)板巖。研究區(qū)侵入巖主要為印支期和燕山期二長(zhǎng)花崗巖，分布于研究區(qū)中部偏西處。區(qū)內(nèi)發(fā)育有多條走向近SN和NW的斷裂(圖1)。

圖1 研究區(qū)地質(zhì)簡(jiǎn)圖Fig.1 Schematic geological map of the study area

圖2 水系沉積物樣品Al2O3、K2O、Ti、MgO含量概率分布型式Fig.2 The frequency distributions of Al2O3,K2O, Ti and MgO of stream sediment samples

研究區(qū)巖性組成主要為灰?guī)r、泥灰?guī)r、砂巖和二長(zhǎng)花崗巖等。巖性背景差異導(dǎo)致水系沉積物中元素含量的概率分布型式多呈偏態(tài)(圖2)，如K2O和MgO均呈明顯的右偏分布，Ti呈明顯的左偏分布，且K2O和Ti還具有多峰分布特征；Al2O3雖然呈正態(tài)分布，但可能也是多個(gè)正態(tài)分布疊加的結(jié)果。為了合理地劃分研究區(qū)水系沉積物地球化學(xué)背景與異常，需消除巖性背景差異的影響。

3.2 水系沉積物樣品分類

首先，選擇可以反映巖性變化，且元素間相關(guān)性較小的非成礦的常量和親石微量元素作為分類指標(biāo)?；谏鲜鲈瓌t，本次實(shí)驗(yàn)共挑選出7種元素，即Mn、Ti、V、Zr、CaO、K2O和MgO，對(duì)其他受巖性影響較小的元素和成礦元素W、Sn則不予考慮。利用上述分類指標(biāo)，通過(guò)EM算法對(duì)樣品進(jìn)行分類，分類過(guò)程主要在MATLAB中實(shí)現(xiàn)。

AIC值計(jì)算結(jié)果表明，隨著分類數(shù)的增大，AIC值先是迅速減小，當(dāng)分類數(shù)增大到一定程度后，AIC值的下降幅度變得越來(lái)越小，最后趨于平穩(wěn)，即AIC曲線趨于水平(圖3)。根據(jù)AIC判別原則，AIC值達(dá)到最小時(shí)可獲得最優(yōu)分類數(shù)，即分類數(shù)越多越好，可在最大程度上區(qū)分混合正態(tài)母體。然而，在分類過(guò)程中還需考慮實(shí)際情況，因?yàn)榉诸愡^(guò)細(xì)可能會(huì)使某些分類中的樣品數(shù)過(guò)少而達(dá)不到基本統(tǒng)計(jì)要求。

觀察發(fā)現(xiàn)，分類數(shù)8前后的AIC值變化幅度較大，因此可以考慮8及以上的分類數(shù)；在分類過(guò)程中，分類數(shù)達(dá)到8時(shí)部分類別中的樣品數(shù)過(guò)少，甚至只有個(gè)位數(shù)樣品；另外，研究區(qū)巖性組成相對(duì)簡(jiǎn)單，因此分類數(shù)不宜過(guò)大?；谏鲜隼碛?，在分類數(shù)8的基礎(chǔ)上將一些樣品數(shù)過(guò)少的類合并到了附近的類中，最終將樣品分為5類(圖4)。由于EM算法在計(jì)算過(guò)程中容易出現(xiàn)局部最優(yōu)解的情況，因此在每次分類過(guò)程中進(jìn)行了多次重復(fù)計(jì)算。將分類結(jié)果成圖，比較分類圖與地質(zhì)圖(圖1)的吻合程度，發(fā)現(xiàn)分類結(jié)果較為合理。整體上，分類界線與研究區(qū)主要地質(zhì)單元界線吻合較好，但也存在同一地質(zhì)單元內(nèi)出現(xiàn)多種分類的現(xiàn)象，如圖幅中部的上古生界灰?guī)r、泥灰?guī)r分布區(qū)出現(xiàn)了第1類和第3類樣品，這是同一地質(zhì)單元內(nèi)具有不同巖性單元的真實(shí)反映。另外，還存在同一分類區(qū)域?qū)?yīng)著不同地質(zhì)單元的現(xiàn)象，如圖幅中第4類樣品分布區(qū)對(duì)應(yīng)著白堊系砂巖、泥盆系砂巖和下古生界淺變質(zhì)砂巖?？偟膩?lái)說(shuō)，將水系沉積物樣品分為5類較符合實(shí)際情況，其中，第1類樣品對(duì)應(yīng)灰?guī)r，第2類樣品對(duì)應(yīng)燕山期二長(zhǎng)花崗巖，第3類樣品對(duì)應(yīng)泥灰?guī)r，第4類樣品對(duì)應(yīng)砂巖，第5類樣品對(duì)應(yīng)印支期二長(zhǎng)花崗巖。

圖3 分類數(shù)與AIC值關(guān)系曲線Fig.3 Relationship between number of clusters and AIC values

分類后絕大多數(shù)母體中的樣品元素含量都服從(近似)正態(tài)分布(圖5)，說(shuō)明分解多重母體取得了較好的效果，在確定地球化學(xué)背景與異常的過(guò)程中可直接采用傳統(tǒng)的均值加k倍標(biāo)準(zhǔn)偏差的方法。值得注意的是，某些區(qū)域內(nèi)可能存在水系沉積物混合物，是由來(lái)自不同巖性背景區(qū)的水系沉積物混合形成的。研究表明，多重母體分解后，這些水系沉積物混合物中的多數(shù)元素含量仍不服從(近似)正態(tài)分布[5]。針對(duì)這種情況，目前還沒(méi)有較好的處理方法。在這種情況下可作取對(duì)數(shù)、平方根、倒數(shù)等變換或直接當(dāng)作正態(tài)分布處理，這是一個(gè)策略選擇問(wèn)題。

圖4 分類結(jié)果Fig.4 Optimal classification map

圖5 分類后Al2O3、K2O、Ti、MgO含量概率分布型式Fig.5 The frequency distributions of Al2O3,K2O, Ti and MgO of separated samples

3.3 異常圈定及效果評(píng)價(jià)

為了便于比較，筆者采用EM方法和傳統(tǒng)的均值加k倍標(biāo)準(zhǔn)偏差的方法分別圈定了W和Sn的地球化學(xué)異常(圖6)。由于W和Sn的概率分布型式為非正態(tài)，因此采用傳統(tǒng)方法處理之前先對(duì)原始數(shù)據(jù)進(jìn)行了對(duì)數(shù)變換，然后采用3S檢驗(yàn)法對(duì)所有樣品數(shù)據(jù)進(jìn)行檢驗(yàn)，剔除異常樣品后再計(jì)算均值和標(biāo)準(zhǔn)偏差，以此進(jìn)行Z分?jǐn)?shù)標(biāo)準(zhǔn)化。數(shù)據(jù)標(biāo)準(zhǔn)化后，將背景上限確定為均值加2倍標(biāo)準(zhǔn)偏差，然后進(jìn)行異常圈定。EM方法計(jì)算步驟為：對(duì)每個(gè)子類的樣品數(shù)據(jù)先采用3S檢驗(yàn)方法剔除異常樣品，計(jì)算均值和標(biāo)準(zhǔn)偏差，然后對(duì)各子類中的數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化，以消除巖性背景差異的影響，合并標(biāo)準(zhǔn)化之后的各個(gè)子類數(shù)據(jù)按統(tǒng)一背景上限進(jìn)行異常圈定。

結(jié)果表明，兩種方法圈定的地球化學(xué)異常存在明顯差異，消除巖性背景后圈定的異常無(wú)論從強(qiáng)度上還是形態(tài)上均發(fā)生了顯著變化?？偟膩?lái)說(shuō)，分類后圈定的W和Sn的異常與已知礦床(點(diǎn))的位置對(duì)應(yīng)關(guān)系更好，而常規(guī)方法圈定的異常有偏移現(xiàn)象。另外，本文方法圈定的異常強(qiáng)度得到了明顯增強(qiáng)，不僅發(fā)現(xiàn)了一些傳統(tǒng)方法沒(méi)有發(fā)現(xiàn)的低弱異常，如B1、B2、B3、B4和B5，還消除了巖性背景高值區(qū)形成的虛假異常，如二長(zhǎng)花崗巖分布區(qū)的A1和A2。這些虛假異常是由巖性背景造成的。統(tǒng)計(jì)發(fā)現(xiàn)，不同巖性分布區(qū)的元素背景值差異較大(表1)，砂巖分布區(qū)的元素背景值最低，其次是灰?guī)r和泥灰?guī)r分布區(qū)，而燕山期二長(zhǎng)花崗巖分布區(qū)的元素背景值最高，印支期二長(zhǎng)花崗巖次之，若全區(qū)統(tǒng)一進(jìn)行異常圈定，顯然會(huì)壓低砂巖分布區(qū)的異常，放大二長(zhǎng)花崗巖分布區(qū)的異常。

受巖性背景影響較大的W和Sn，兩種方法圈定的異常差異較大。由于二長(zhǎng)花崗巖分布區(qū)的元素背景值明顯高于其他巖性分布區(qū)，因此傳統(tǒng)方法圈定的W、Sn異常多位于二長(zhǎng)花崗巖分布區(qū)，而在灰?guī)r、泥灰?guī)r及砂巖分布區(qū)幾乎沒(méi)有顯示。相比較而言，本文方法消除了巖性背景差異的影響，圈定出的W、Sn異常主要分布于二長(zhǎng)花崗巖巖體的周邊，且在砂巖、灰?guī)r及泥灰?guī)r分布區(qū)也有顯示。另外，本文方法圈定的異常規(guī)模也較大，強(qiáng)度也較高，與已知礦床(點(diǎn))的對(duì)應(yīng)關(guān)系更好，說(shuō)明利用該方法圈定的異常更為合理。

圖6 傳統(tǒng)方法(左側(cè))與EM方法(右側(cè))圈定的W、Sn地球化學(xué)異常Fig.6 Geochemical anomalies of W and Sn delineated by traditional method(left) and the EM method(right)

表1 分類樣品中W、Sn背景值

4 結(jié)論

1) 采用EM分類方法能夠快速、有效地區(qū)分來(lái)自不同巖性背景的水系沉積物樣品，可在一定程度上消除巖性背景的影響。該方法以水系沉積物對(duì)基巖化學(xué)成分顯著的繼承性為基礎(chǔ)，具有明確的地質(zhì)意義。

2) 消除巖性背景影響后，不僅可以圈定出低弱地球化學(xué)異常，還可以削弱虛假地球化學(xué)異常。

3) 分類完成后，對(duì)各類樣品元素含量進(jìn)行正態(tài)分布檢驗(yàn)，若服從(近似)正態(tài)分布，則可直接采用傳統(tǒng)的均值加k倍標(biāo)準(zhǔn)偏差的方法或其他方法圈定地球化學(xué)異常，若不服從(近似)正態(tài)分布，則可作取對(duì)數(shù)、平方根、倒數(shù)等變換或直接當(dāng)作正態(tài)分布處理，這是一個(gè)策略選擇問(wèn)題。