協(xié)同克立格法在廣西林旺金礦化探數(shù)據(jù)處理中的應(yīng)用

徐銀河，孫春巖，竺瑋煌，趙志濤，楊 奎，劉 鈞

(中國地質(zhì)大學(xué) 工程技術(shù)學(xué)院，北京100083)

協(xié)同克立格法在廣西林旺金礦化探數(shù)據(jù)處理中的應(yīng)用

徐銀河，孫春巖，竺瑋煌，趙志濤，楊 奎，劉 鈞

(中國地質(zhì)大學(xué) 工程技術(shù)學(xué)院，北京100083)

協(xié)同克立格法，同時(shí)兼具化探數(shù)據(jù)多元性及克立格法表征空間屬性的特點(diǎn)?？紤]到地質(zhì)變量兩個(gè)以上的空間屬性，在化探數(shù)據(jù)處理中運(yùn)用協(xié)同克立格法，可以進(jìn)一步提高估計(jì)精度。運(yùn)用協(xié)同克立格法對(duì)廣西林旺礦區(qū)中金元素成礦進(jìn)行預(yù)測(cè)，經(jīng)相關(guān)性分析顯示，在礦區(qū)中元素的基本組合是Au、Ag、As、Hg，協(xié)同克立格插值以Au為主要變量，伴生元素Ag、As、Hg為次要變量。把插值計(jì)算后的結(jié)果與傳統(tǒng)多元統(tǒng)計(jì)方法和普通克立格法的計(jì)算結(jié)果進(jìn)行比較，結(jié)果協(xié)同克立格法得到的估值誤差較小，預(yù)測(cè)精度較高，在成礦元素的預(yù)測(cè)中具有一定程度的優(yōu)越性。

協(xié)同克立格法;相關(guān)分析;變差函數(shù);估計(jì)精度

0 前言

在金屬礦化探數(shù)據(jù)處理中，傳統(tǒng)的多元統(tǒng)計(jì)分析處理實(shí)際問題效果比較好，但其缺點(diǎn)是隨機(jī)理論的一些假設(shè)前提，在地質(zhì)上難以滿足，并且不能分析地質(zhì)變量所特有的空間屬性。在地質(zhì)統(tǒng)計(jì)學(xué)中，最重要的方法，即克立格法是一種線性、無偏、最小估計(jì)方差的估計(jì)方法。它是一種既能保持概率統(tǒng)計(jì)方法的有效性，又更切合地質(zhì)變量特點(diǎn)的方法途徑?？肆⒏穹◤膯我蛔兞砍霭l(fā)的插值計(jì)算，避免了系統(tǒng)誤差，可以做到最優(yōu)的估計(jì)。由于化探數(shù)據(jù)多元性的特點(diǎn)，使得運(yùn)用協(xié)同克立格法可以做到對(duì)多個(gè)相關(guān)變量的克立格估值，提高預(yù)測(cè)精度。作者利用協(xié)同克立格法對(duì)廣西林旺金礦數(shù)據(jù)資料進(jìn)行了分析，與傳統(tǒng)多元統(tǒng)計(jì)及普通克立格法的計(jì)算結(jié)果比較，對(duì)計(jì)算得到的誤差精度進(jìn)行了分析對(duì)比［1、2]。

1 礦區(qū)原生暈地球化學(xué)特征

本區(qū)礦體主要賦存于三疊系百逢組，少量分布于板納組的碎屑巖擠壓、揉皺帶中。礦體多成似層狀，脈狀、透鏡狀分布于揉皺強(qiáng)烈的砂巖、泥巖斷裂、節(jié)理、劈理裂隙帶中，以及褶皺軸部巖層受擠壓的裂隙帶中。

在礦區(qū)內(nèi)，褶皺和斷裂都較發(fā)育，幾乎所有成礦成暈元素都是沿北西向東翼壓扭性斷裂分布。其中，尤以Au、Ag、As、Hg表現(xiàn)最為明顯，說明了北西方向控礦構(gòu)造的存在，Au異常受控礦構(gòu)造的控制。土壤與巖石中的Mo、W、Bi，可能指示礦體線剝蝕程度，而Ni、Co的出現(xiàn)，可能是礦體剝蝕更深的指示。

至今，林旺礦區(qū)氧化礦的主要礦段已基本采完，而原生礦尚未完全探明和開采利用，故對(duì)深部原生礦進(jìn)一步的勘查，對(duì)礦區(qū)今后的發(fā)展具有重大的意義。

礦床原生暈最重要的特點(diǎn)之一，是它具有分帶性。造成分帶的原因，主要是成礦溶液在運(yùn)移過程中，其組份濃度的不斷變化，各種元素化學(xué)活動(dòng)性及析出順序的不同，往往導(dǎo)致了同種指示元素的含量及不同指示元素的含量，在空間上有規(guī)律的變化，這種現(xiàn)象被稱為原生暈的分帶性。在以往的工作中，已應(yīng)用原蘇聯(lián)學(xué)者格里戈亮計(jì)算原生暈中元素分帶序列的方法，對(duì)礦山礦段的部份原生暈樣品，計(jì)算了元素的分帶序列，獲得礦山礦段的元素分帶序列為:

(1)前緣(上部)元素組合(As－Sb－Hg－Au)。

(2)近礦(中部)元素組合(Hg－W－Au－Ag－Mo)。

(3)尾部(下部)元素組合(Cu－Co－Bi－Ni－Zn－Pb)。

可以看到，不同指示元素在空間上是按一定的規(guī)律變化分布的［3]。

原生暈中元素的分帶性，指示了元素的空間變化規(guī)律，所以在化探數(shù)據(jù)處理中，需要考慮到元素的空間屬性。而傳統(tǒng)的多元統(tǒng)計(jì)方法基于隨機(jī)理論，不能夠完整地表達(dá)出地質(zhì)變量的空間相關(guān)性。而協(xié)同克立格法，既能反映變量的空間屬性，同時(shí)也考慮了元素多元性的特點(diǎn)，更加貼合地質(zhì)變量的實(shí)際情況。

2 協(xié)同克立格法

克立格法是以變差函數(shù)理論和結(jié)構(gòu)分析為基礎(chǔ)，對(duì)區(qū)域化變量進(jìn)行線性、無偏、最優(yōu)內(nèi)插估計(jì)的方法。它是從研究礦產(chǎn)儲(chǔ)量計(jì)算及其誤差估計(jì)問題而產(chǎn)生和發(fā)展起來的，在化探數(shù)據(jù)處理中，應(yīng)用克立格法考慮了化探異常的空間屬性，可以避免系統(tǒng)誤差，同時(shí)提高估計(jì)精度。

化探數(shù)據(jù)具有多元性的特點(diǎn)，而普通克立格法研究區(qū)域化變量只從單一變量出發(fā)，不能夠很好地滿足化探數(shù)據(jù)處理的需求。用協(xié)同克立格法進(jìn)行分析評(píng)價(jià)，可以很好地兼顧克立格法與數(shù)據(jù)多元性的特點(diǎn)。協(xié)同克立格法是把區(qū)域化變量的最佳估值方法，從單一屬性發(fā)展到二個(gè)以上的協(xié)同區(qū)域化屬性。協(xié)同克立格法在插值時(shí)，用一個(gè)或更多的次要變量，這些次要變量與主要變量有相關(guān)關(guān)系。變量之間的相關(guān)關(guān)系，能用于提高主要變量的預(yù)測(cè)精度［2、4]。

3 研究區(qū)元素統(tǒng)計(jì)特征

通過統(tǒng)計(jì)分析后可以運(yùn)用變差函數(shù)建立模型，進(jìn)而用克立格法估值。在計(jì)算變差函數(shù)之前，必須對(duì)樣品數(shù)據(jù)進(jìn)行分析和處理。樣品特異值(即樣品中的特高含量值)的存在，將嚴(yán)重影響變異函數(shù)的結(jié)構(gòu)性。在對(duì)原始數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析前，將含量異常的特異值剔除掉，在得到的數(shù)據(jù)中對(duì)Au元素含量進(jìn)行對(duì)數(shù)轉(zhuǎn)化后服從對(duì)數(shù)正態(tài)分布，見圖1，所以可以用克立格法進(jìn)行插值估計(jì)。

圖1 林旺金礦Au含量對(duì)數(shù)統(tǒng)計(jì)直方圖Fig．1 Logarithmic histogram of gold element

首先，對(duì)林旺金礦主要的十三種元素進(jìn)行統(tǒng)計(jì)分析，其元素含量的相關(guān)矩陣見下頁表1。通過計(jì)算相關(guān)系數(shù)，可以得出:

(1)主成礦元素Au與Ag、As、W、Hg元素在含量變化上呈正相關(guān)特性。

(2)Au與Sb、Mo元素在含量變化上呈相對(duì)微弱的正相關(guān)特性。

(3)Au與Bi、Co、Cu、Ni、Pb、Zn元素呈負(fù)相關(guān)特性。

Au元素的異常組合，可以通過正交旋轉(zhuǎn)主因子解進(jìn)一步確定，見下頁表2。從表2中可以看出，Au元素與Ag、As和Hg元素關(guān)系密切［5]。

4 結(jié)果分析

4．1 建立和計(jì)算變差函數(shù)模型

協(xié)同克立格法是研究一個(gè)變量取樣值時(shí)，借助其它與之相關(guān)的變量樣品信息，來提高對(duì)該變量的估計(jì)精度。協(xié)同克立格法利用了更多的樣品信息，比較普通克立格法而言，計(jì)算量更大一些。通過統(tǒng)計(jì)分析得出，Au元素與Ag、As和Hg元素關(guān)系密切，進(jìn)而運(yùn)用協(xié)同克立格法對(duì)其估值。

統(tǒng)計(jì)分析的結(jié)果表明，Au元素的含量服從對(duì)數(shù)正態(tài)分布，所以在計(jì)算變差函數(shù)之前，要對(duì)Au元素的含量數(shù)據(jù)進(jìn)行對(duì)數(shù)轉(zhuǎn)換。在統(tǒng)計(jì)分析的基礎(chǔ)上，首先要建立變差函數(shù)模型，變差函數(shù)的重要性在于獲得了表征區(qū)域化變量空間特性的三個(gè)參數(shù)，即變程、塊金效應(yīng)、基臺(tái)值，這三個(gè)參數(shù)可用于求取克立格計(jì)算時(shí)所需的權(quán)系數(shù)，所以變差函數(shù)能較好地反映成礦異常的空間結(jié)構(gòu)性和隨機(jī)性。

表1 林旺礦區(qū)十三種元素含量的相關(guān)矩陣Tab．1 Correlation matrix of 13 elements

表2 正交旋轉(zhuǎn)主因子矩陣Tab．2 Rotated componentmatrix

協(xié)同克立格變差函數(shù)如下:

協(xié)同克立格法的任務(wù)，是應(yīng)用估計(jì)領(lǐng)域內(nèi)的有效數(shù)據(jù)，來估計(jì)中心點(diǎn)在x0的待估域V(x0)上，變量k0的平均值Zvk0的估計(jì)量

在插值中，要考慮到領(lǐng)域樣品的搜索形式及搜索半徑的選擇，通常的搜索方式有橢球體和球體二種。由于克立格法的適用條件需滿足平穩(wěn)性假設(shè)與內(nèi)蘊(yùn)假設(shè)條件，所以選擇用于區(qū)域化變量在空間平穩(wěn)分布的球狀模型。這就要求在應(yīng)用時(shí)，對(duì)所研究的礦床有一定了解，取樣數(shù)據(jù)點(diǎn)不能太少［6]。

在分析變差函數(shù)特點(diǎn)的基礎(chǔ)上，確定用球狀模型進(jìn)行理論變差函數(shù)擬合。變差函數(shù)曲線反映了一個(gè)采樣點(diǎn)與其相鄰采樣點(diǎn)的空間關(guān)系，對(duì)異常采樣點(diǎn)具有很好的探測(cè)作用。在圖2(見下頁)中，點(diǎn)圖為實(shí)驗(yàn)變差函數(shù)，通過計(jì)算樣品值的變差函數(shù)，得到了理論變差函數(shù)。在圖2中的曲線，即為理論變差函數(shù)曲線。從圖2中可以看出，實(shí)驗(yàn)變差函數(shù)與理論變差函數(shù)曲線擬合情況良好。

Au元素樣品點(diǎn)分布不具有趨勢(shì)，擬合的變差函數(shù)不具有幾何各向異性，即各個(gè)方向的變差函數(shù)均具有相同的基臺(tái)值和變程。當(dāng)采樣點(diǎn)間的距離增大時(shí)，變異函數(shù)從初始的塊金值達(dá)到一個(gè)相對(duì)穩(wěn)定的常數(shù)時(shí)，該常數(shù)值稱為基臺(tái)值。當(dāng)變差函數(shù)值超過基臺(tái)值時(shí)，即函數(shù)值不隨采樣點(diǎn)間隔距離而改變時(shí)，空間相關(guān)性不存在。當(dāng)變差函數(shù)的取值由初始的塊金值達(dá)到基臺(tái)值時(shí)，采樣點(diǎn)的間隔距離被稱為變程。變程表示在某種觀測(cè)尺度下，空間相關(guān)性的作用范圍，其大小受觀測(cè)尺度的限定。塊金值就是在二采樣點(diǎn)非常接近時(shí)，它們的變差函數(shù)值不為0，即存在塊金值。球狀模型的變程是0．233 833，偏基臺(tái)值是0．518 81，塊金值是1.161 9，滯后距是0．040 936，滯后數(shù)是12。

4．2 變差函數(shù)的交叉驗(yàn)證

用克立格法估值時(shí)，需要將原始數(shù)據(jù)分成二部份，即計(jì)算數(shù)據(jù)和檢驗(yàn)數(shù)據(jù)。其原理是用選定的理論變差函數(shù)的參數(shù)，對(duì)已知數(shù)據(jù)利用其周圍的數(shù)據(jù)進(jìn)行克立格估值，將所得的這套估計(jì)值與原始值進(jìn)行比較，并對(duì)其結(jié)果進(jìn)行統(tǒng)計(jì)分析。良好的檢驗(yàn)結(jié)果將是二套數(shù)據(jù)的均值之差很小，方差也很小，以及有良好的正相關(guān)性，從而判斷參數(shù)的正確性。

驗(yàn)證之后的數(shù)據(jù)見表3和圖3。

從表3可以看出，標(biāo)準(zhǔn)平均值非常接近于0，均方根預(yù)測(cè)誤差也很小，平均標(biāo)準(zhǔn)誤差與均方根預(yù)測(cè)誤差非常接近，均符合最優(yōu)變差函數(shù)模型的要求?？梢缘贸鼋Y(jié)論:計(jì)算變差函數(shù)模型是最優(yōu)的;對(duì)樣品含量的估計(jì)是無偏的;滿足區(qū)域化變量的內(nèi)蘊(yùn)假設(shè);驗(yàn)證該變差函數(shù)的模擬是正確的。

圖2 實(shí)驗(yàn)變差函數(shù)與理論變差函數(shù)擬合圖Fig．2 Experimental variogram and theoretic variogram

表3 交叉驗(yàn)證模型參數(shù)Tab．3 Parameters of cross validation

圖3 交叉驗(yàn)證的預(yù)測(cè)圖Fig．3 Cross validation prediction

同時(shí)，利用普通克立格法對(duì)礦區(qū)化探數(shù)據(jù)進(jìn)行分析評(píng)價(jià)，從表4可以看到，用協(xié)同克立格法對(duì)原始數(shù)據(jù)進(jìn)行估值時(shí)所得到的預(yù)測(cè)誤差，要比普通克立格法得到的預(yù)測(cè)誤差減少11.11%，標(biāo)準(zhǔn)誤差減少23.38%。由此可見，協(xié)同克立格法與普通克立格法相比，在一定程度上，協(xié)同克立格法是較為優(yōu)越的，這對(duì)于提高地質(zhì)統(tǒng)計(jì)學(xué)的預(yù)測(cè)精度，起到了一定的作用。但該方法唯一缺點(diǎn)就是計(jì)算量比較大，不過可以通過先進(jìn)的軟件來解決。

通過對(duì)變差函數(shù)模型的計(jì)算以及得到的精度，可以對(duì)林旺礦區(qū)的Au元素進(jìn)行估計(jì)，得到的預(yù)測(cè)圖見圖4［4、7]。

圖4 林旺礦區(qū)Au元素異常預(yù)測(cè)圖(1∶10 000，ArcGIS)Fig．4 Gold element abnormity predication

表4 Au元素協(xié)同克立格法與普通克立格法交叉驗(yàn)證結(jié)果比較Tab．4 Cokriging and normal kriging of gold element in cross validation

5 討論

(1)協(xié)同克立格法結(jié)合了化探數(shù)據(jù)處理的多元性，以及克立格法反映空間屬性的特點(diǎn)，與傳統(tǒng)多元統(tǒng)計(jì)分析相比，較能夠反映出地質(zhì)上區(qū)域化變量的空間特征。

(2)在化探數(shù)據(jù)處理中，運(yùn)用協(xié)同克立格法估值，在一定程度上提高了預(yù)測(cè)的精度，得到了較為精確的估計(jì)精度。與普通克立格法相比較，協(xié)同克立格滿足了多變量估值的要求，更加符合化探數(shù)據(jù)統(tǒng)計(jì)的要求，缺點(diǎn)就是計(jì)算量比較大，需要運(yùn)用先進(jìn)的軟件來處理。

(3)運(yùn)用協(xié)同克立格法研究區(qū)域化變量時(shí)，要滿足平穩(wěn)假設(shè)和內(nèi)蘊(yùn)假設(shè)條件，所以需要對(duì)所研究礦床有一定的了解，取樣數(shù)據(jù)點(diǎn)不能太少。

［1]余先川，王世稱，王桂安．穩(wěn)健協(xié)同克里格因子分析及其在化探中的應(yīng)用［J]．地球科學(xué)—中國地質(zhì)大學(xué)學(xué)報(bào)，1998，23(2):1．

［2]王仁鐸，胡光道．線性地質(zhì)統(tǒng)計(jì)學(xué)［M]．武漢:武漢地質(zhì)學(xué)院教材科，1984．

［3]吳錫生．廣西樂業(yè)縣林旺金礦勘查地球化學(xué)深部找礦工作報(bào)告［D]．北京，2009．

［4]伍偉，秦德先，胡志軍，等．地質(zhì)統(tǒng)計(jì)學(xué)在曼家寨錫鋅多金屬礦床中的應(yīng)用［J]．有色金屬，2008，60(2):2．

［5]葉水盛，劉光盛，馬生忠．克立格方法在區(qū)域化探數(shù)據(jù)處理中的應(yīng)用［J]．長春科技大學(xué)學(xué)報(bào)，2000，30(2):1．

［6]楊建宇，秦德先，康澤寧等．地質(zhì)統(tǒng)計(jì)學(xué)在北衙金礦的應(yīng)用［J]．礦業(yè)研究與發(fā)展，2006，26(2):3．

［7]湯國安，楊昕．ArcGIS地理信息系統(tǒng)空間分析實(shí)驗(yàn)教程［M]．北京:科學(xué)出版社，2006．

P 632

A

1001—1749(2011)02—0179—05

廣西國土資源廳普查項(xiàng)目(4500000330364)

2010－09－03改回日期:2010－12－08

徐銀河(1985－)，女，碩士，主要從事基于GIS的數(shù)學(xué)地質(zhì)方法在礦產(chǎn)資源評(píng)價(jià)中的應(yīng)用。