李 亮，徐娟娟，葉會壽

(1.昆明理工大學(xué) 國土資源工程學(xué)院，云南 昆明 650093)

(2.洛陽欒川鉬業(yè)集團股份有限公司，河南 欒川 471542)

(3.中國地質(zhì)科學(xué)院礦產(chǎn)資源研究所，北京 100073)

0 前言

地質(zhì)統(tǒng)計學(xué)是法國著名統(tǒng)計學(xué)家G.Matheron于20 世紀(jì)60 年代在大量理論研究的基礎(chǔ)上創(chuàng)立并發(fā)展的一門集數(shù)學(xué)與地質(zhì)采礦為一體的邊緣學(xué)科，目前已在國內(nèi)外獲得廣泛應(yīng)用。它是以區(qū)域化變量理論作為理論基礎(chǔ)，以變異函數(shù)為主要工具，采用不同的克立格方法，對在空間上既有隨機性又有結(jié)構(gòu)性的變量(如品位值)進行研究的科學(xué)。經(jīng)過40 多年的發(fā)展，目前已經(jīng)形成了一套較為完整的理論體系［1-2］。本文利用礦業(yè)工程軟件，運用地質(zhì)統(tǒng)計學(xué)原理和方法，對三道莊鉬礦床進行研究。采用最優(yōu)、線性無偏估計研究鉬礦石品位的空間變異性，并計算給定邊界品位下的礦石的儲量，為礦山資源合理開發(fā)利用提供有意義的決策依據(jù)。

1 礦床地質(zhì)概況

三道莊鉬礦床位于河南省欒川縣境內(nèi)，礦床中心地理坐標(biāo)為東經(jīng)111°29'40″、北緯33°55'00″，鉬資源量豐富，為特大型鉬礦床［3］。該礦床位于華北地臺南緣褶皺帶盧氏—欒川多金屬成礦帶中的南泥湖鉬礦田內(nèi)，區(qū)內(nèi)地層、構(gòu)造、變質(zhì)帶及巖漿巖均呈北西西-北西向分布。三道莊鉬礦主要礦體為大厚度的似層狀礦體，呈扇形產(chǎn)于巖體外接觸帶三川組上段的矽卡巖、鈣硅酸角巖及南泥湖中段的黑云母長英角巖中。在主礦體的上、下部位有少數(shù)零星小礦體分布。此外，在主礦體頂、底板有分枝復(fù)合及褶皺構(gòu)造轉(zhuǎn)折處有膨脹狹縮等現(xiàn)象［3-4］。礦體形態(tài)、產(chǎn)狀受巖性、構(gòu)造、巖漿巖控制較為明顯。礦體的分布范圍、產(chǎn)狀與鈣硅酸角巖、矽卡巖及巖漿巖緩傾斜部位的分布范圍、產(chǎn)狀相一致。礦體走向280°～310°，傾向南西，傾角較為平緩。主礦體沿走向在礦區(qū)范圍內(nèi)長度大于1 420 m，沿傾向長度大于1 120 m，厚度一般80～150 m，最大厚度可達364.56 m，厚度變化系數(shù)為44%～60%［3-4］(圖1)。

圖1 三道莊鎢鉬礦區(qū)橫9 線勘探剖面

2 數(shù)據(jù)選取及數(shù)據(jù)庫的建立

根據(jù)研究需要，選取三道莊鉬礦體主要成礦組分鉬為研究對象，共收集了礦體上146 個地表勘探鉆孔共21 107 件樣品，鉬化驗分析數(shù)據(jù)208 97 件。這些樣品大多來自河南省地質(zhì)礦產(chǎn)局第一地質(zhì)調(diào)查大隊在20 世紀(jì)70 年代礦床詳細勘探過程中所采集的樣品，部分為后期礦山基建勘探所采集的樣品。將這些數(shù)據(jù)以特定格式錄入計算機中，利用Surpac專業(yè)軟件將收集到的每一個樣品化驗數(shù)據(jù)連同其對應(yīng)的三維坐標(biāo)及其他屬性(樣號、化驗數(shù)據(jù)、巖性、勘探工程編號、勘探時間等)存入數(shù)據(jù)庫中。

3 數(shù)據(jù)的統(tǒng)計分析

在地質(zhì)統(tǒng)計學(xué)中，計算實驗半變異函數(shù)之前，首先按照有效數(shù)據(jù)須確定在建立長度的支撐上的要求，即對參與計算的所有樣品數(shù)據(jù)進行重新組合，使同類參數(shù)的地質(zhì)樣品段的長度達到一致［2，5］。因此，必須對樣長不等的鉆孔原始數(shù)據(jù)按一定的長度進行重新組合計算，并對原始數(shù)據(jù)進行等長組合處理。在組合過程中，筆者根據(jù)三道莊礦床鉆孔的特點和礦體研究的技術(shù)要求，采用2 m 的長度進行樣品組合，并對組合后的樣品進行統(tǒng)計分析，鉬品位分布特征結(jié)果參數(shù)見表1。

表1 三道莊鉬礦床組合樣品統(tǒng)計分析結(jié)果

由于有用元素在地質(zhì)環(huán)境中自然富集，造成部分樣品元素品位相對較高，不能真實地反映礦體的實際分布特征，因此需要進行特高品位的處理［5］。根據(jù)礦山勘探報告對特異值的識別及處理原則，筆者選取97.5%的置信區(qū)間處的樣品值0.296%，將其定為特高品位下限，對特高品位進行處理，然后對鉬品位分布進行統(tǒng)計分析。研究發(fā)現(xiàn)，鉬品位分布呈現(xiàn)非正偏斜，表現(xiàn)為右偏，不服從正態(tài)分布(圖1-a)，將組合樣品位進行對數(shù)轉(zhuǎn)換后，采用正態(tài)分布的峰度偏度檢驗法，發(fā)現(xiàn)鉬組合樣品位基本服從三參數(shù)的對數(shù)正態(tài)分布(圖1-b)。因此可以利用轉(zhuǎn)化后的數(shù)據(jù)進行變異函數(shù)計算，并且可以選用克立格方法對礦體進行品位估值和儲量計算。

圖2 鉬品位及鉬品位經(jīng)對數(shù)轉(zhuǎn)換后的分布直方圖

4 變異函數(shù)模型的建立

4.1 實驗變異函數(shù)的計算

在經(jīng)典統(tǒng)計學(xué)中，通常采用均值、方差、標(biāo)準(zhǔn)偏差、變化系數(shù)等一類參數(shù)來概括反映地質(zhì)體的全貌(如各元素總的變化情況)，卻無法反映局部范圍和特定方向上地質(zhì)特征的變化，這是因為經(jīng)典統(tǒng)計學(xué)只考慮數(shù)值大小，不考慮數(shù)據(jù)的空間位置關(guān)系。而在地質(zhì)統(tǒng)計學(xué)中，研究礦床中金屬品位分布特征的一個重要的方面就是掌握品位沿某一個或某幾個方向或距離的變化規(guī)律，通常就采用變異函數(shù)來分析，它兼顧了品位數(shù)值大小和品位變化的空間位置關(guān)系，使得統(tǒng)計的品位分布特征更趨于合理［6-7］。

所謂變異函數(shù)就是假設(shè)在一個空間中的某個域V，V 內(nèi)區(qū)域變化量按一定方向間隔的兩點x，x+h的變化程度，變異函數(shù)反映了區(qū)域變化量的結(jié)構(gòu)特征和空間連續(xù)性。但是在實踐中，往往樣品的數(shù)目總是有限的，把用有限的實際采集的樣品值構(gòu)制成的變異函數(shù)稱為實驗變異函數(shù)，實驗變異函數(shù)就是根據(jù)實際所取數(shù)據(jù)構(gòu)造變異函數(shù)γ(h)的估計值γ*(h)

式中，γ*(h)為實驗變異函數(shù);h 為滯后距，即在一定方向上，距離為| h |的矢量;N(h)表示滯后距為h 時參加實驗變異函數(shù)計算的樣品對數(shù);Z(xi)，Z(xi+h)為區(qū)域化變量在空間點上的品位的測定值［8-9］。

對組合樣品位數(shù)據(jù)進行變異函數(shù)曲線分析，分析其空間變異特征，是研究礦體變化性以及計算礦體儲量必不可少的步驟［6］。由于三道莊鉬品位數(shù)據(jù)服從三參數(shù)的對數(shù)正態(tài)分布，因此在本次計算三道莊鉬品位實驗變異函數(shù)時，首先需將鉬元素組合樣的品位數(shù)據(jù)進行對數(shù)轉(zhuǎn)換，選取經(jīng)過對數(shù)轉(zhuǎn)換處理后的數(shù)據(jù)作為區(qū)域化變量，根據(jù)鉬礦體成礦地質(zhì)特征，選擇分別在礦體走向、傾向、厚度3 個方向，即對應(yīng)為變異函時繪制了實驗變異函數(shù)曲線圖。從圖3、圖4、圖5 中曲線A 可以看出3 個方向的曲線波動性都較小，說明三道莊鉬礦體鉬品位變化較小、分布較均勻。

4.2 理論變異函數(shù)曲線擬合

圖3 鉬品位厚度方向?qū)嶒炞儺惡瘮?shù)曲線及其擬合曲線

圖4 鉬品位傾向方向?qū)嶒炞儺惡瘮?shù)曲線及其擬合曲線

圖5 鉬品位走向方向?qū)嶒炞儺惡瘮?shù)曲線及其擬合曲線

從三道莊鉬礦體鉬品位分布得出的實驗變異函數(shù)曲線的特征來看(如圖3、圖4、圖5 中曲線A)，雖然數(shù)據(jù)點的分布規(guī)律性相對較差，但仍可看出變異函數(shù)γ(h)隨距離h 首先增加，然后趨于穩(wěn)定，符合球狀模型的特點，因而采用實驗變異函數(shù)球狀模型用加權(quán)多項式回歸法分別對走向、傾向、厚度3 個方向進行擬合。球狀模型公式為:

式中C0為塊金常數(shù);C0+C 為基臺值;a 為變程;h為滯后距。

從公式(2)中可以看出，γ(h)隨h 的增加而增加，當(dāng)h 達到變程時，γ(h)達到基臺值C0+C;之后γ(h)，便保持常值C0+C。這種特征的物理意義是:當(dāng)樣品之間的距離小于變程時，樣品是相互關(guān)聯(lián)的，關(guān)聯(lián)程度隨間距的增加而減小，或者說，變異程度隨間距的增加而增大;當(dāng)間距達到一定值時，樣品之間的關(guān)聯(lián)性消失，變?yōu)橥耆S機，這時γ(h)即為樣品的方差。因此，變程實際上代表樣品的影響范圍［9-11］。

通過對三道莊鉬礦體鉬品位3 個方向?qū)嶒炞儺惡瘮?shù)的研究，發(fā)現(xiàn)它們都可以用1 個帶有塊金常數(shù)的球狀模型進行擬合，經(jīng)過反復(fù)實驗，確定三道莊鉬礦體鉬品位的理論變異函數(shù)曲線，如圖3、圖4、圖5中B 所示:

從圖3、圖4、圖5 圖中曲線B 可以看出，三道莊鉬礦體鉬品位變異函數(shù)在3 個方向上的基臺值(C0+C)都基本相同，而變程(a)不同，表現(xiàn)出各向異性的特征。同時，鉬的品位在走向方向變程最大，厚度方向最小，說明了鉬礦化的連續(xù)性在走向方向上較其它兩個方向強，而厚度方向的連續(xù)性最弱，這與礦床地質(zhì)成因是相一致的。求得球狀模型的參數(shù)如下:

其在各個方向上的變異函數(shù)分別表示為:

礦體厚度方向:

從以上3 個變異函數(shù)公式表明三道莊鉬礦床鉬品位變異函數(shù)具有方向上的異向性，即在不同方向上表現(xiàn)為基臺值相同，變程不同的幾何異向性，各向異性比為:

在本次研究中，筆者用帶狀各向異性對模型進行套合(帶狀各向異性模型幾乎可用于任意試驗各向異性模型)，變異函數(shù)的套合結(jié)構(gòu)為:

式中sph(250.299)、sph(350.962)和sph(484.328)分別表示變程為250.299 m、350.962 m 和484.328 m 的球狀模型。其相應(yīng)的轉(zhuǎn)換矩陣為:

通過以上對礦體各方向變異函數(shù)的計算及曲線的模擬，可以確定三道莊鉬礦體鉬品位的變異函數(shù)具有幾何異向性，空間各方向的變程(距離)構(gòu)成了變程各向異性。

4.3 變異函數(shù)交叉驗證

為了對結(jié)構(gòu)模型進行檢驗，采用了交叉驗證方法，以判斷變異函數(shù)擬合參數(shù)的選取是否正確。交差驗證是用來檢驗變異函數(shù)曲線擬合的可靠性，并輔助確定變異函數(shù)參數(shù)。其基本方法是:利用所得結(jié)構(gòu)模型和已知樣品去估已知值，然后把這些真值和估計值進行比較，對兩者的殘差(差值)進行統(tǒng)計分析，以判斷變異函數(shù)結(jié)構(gòu)的正確性［8，9，12］。交叉驗證結(jié)果見表2、圖6、圖7。

表2 殘差統(tǒng)計參數(shù)一覽表

圖6 殘差VS 品位統(tǒng)計圖

圖7 鉬品位交叉檢驗殘差分布直方圖

通過驗證計算，經(jīng)對數(shù)轉(zhuǎn)換后的鉬組合樣品位實際值與估計值之間的誤差均值趨近于零，即誤差的均值趨近于零(圖6);經(jīng)對數(shù)轉(zhuǎn)換后的鉬組合樣品位實際值與估計值之間的單個偏差趨近于零，且趨近于平均克立格方差，即誤差的方差趨近于零且方差比趨近于1;殘差分布屬于正態(tài)分布，且95%置信限位于正負兩倍的克立格方差范圍內(nèi)(圖7)。以上結(jié)果表明模型確定合理，變異函數(shù)參數(shù)對經(jīng)轉(zhuǎn)換后的鉬品位進行估計是無偏的，滿足足區(qū)域化變量內(nèi)蘊假設(shè)，從而驗證了變異函數(shù)模型是正確的，可以用于下一步的礦床鉬品位估值和資源儲量計算。

4.4 克立格法品位估值

完成變異函數(shù)模型結(jié)構(gòu)分析，充分了解了三道莊鉬礦體的空間變異性，下一步工作是設(shè)計克立格估值方案?？肆⒏穹ㄊ歉鶕?jù)代估樣本點(或塊段)有限鄰域內(nèi)若干已測定的樣本點數(shù)據(jù)，考慮了樣本點的形狀、大小和空間相互位置關(guān)系，與待估樣本點的相互空間位置關(guān)系，以及變異函數(shù)提供的結(jié)構(gòu)信息，對待估樣本點值進行的一種線性無偏最優(yōu)估計［14］。

地質(zhì)統(tǒng)計學(xué)主要是在結(jié)構(gòu)分析的基礎(chǔ)上采用各種克立格法來估計和解決實際問題。根據(jù)研究的目的和條件的不同，主要劃分為普通克立格法、泛克立格法、對數(shù)正態(tài)克立格法、指示克立格法等，在不同的條件下應(yīng)選取適當(dāng)?shù)姆椒ㄟM行估值，其中最基本、應(yīng)用最為廣泛的是普通克立格法［10-14］。普通克立格法適用于滿足二階平穩(wěn)(或內(nèi)蘊)假設(shè)的礦床。綜合分析針對三道莊礦床鉬元素品位滿足內(nèi)蘊假設(shè)條件，選擇采用普通克立格法進行礦床品位估值。其克立格主要參數(shù)設(shè)置如下:橢球體定向方位為300°、傾伏角為0°、傾角為18°;各項異性比率為1∶1.38∶1.935;最大距離分別選取110 m，220 m，440 m;最少樣品和最多樣品總數(shù)分別設(shè)置為3 個和16 個。結(jié)合已經(jīng)得到的理論變異函數(shù)參數(shù)值(塊金值、基臺值、變程等)，從而完成對礦床鉬品位的估值。

5 儲量計算

根據(jù)組合樣最低品位值限制，按邊界品位0.03%，礦石體積質(zhì)量取3.2 t/m3進行儲量統(tǒng)計計算。將儲量計算結(jié)果與采用傳統(tǒng)平行斷面法計算結(jié)果進行驗證對比，結(jié)果如圖表3 所示。

從對比結(jié)果可以看出，采用普通克立格法計算礦體資源量與采用傳統(tǒng)平行斷面法計算的礦體資源量基本吻合，其誤差基本在5%左右，說明應(yīng)用普通克立格法對礦體資源量的估算是可靠的。所對應(yīng)的鉬品位—噸位曲線見圖8。

圖8 三道莊鉬礦床鉬品位—噸位曲線

從圖8 可以看出，隨邊界品位的逐漸增高，礦體的平均品位逐漸升高，鉬資源金屬量則不斷減少，因此礦山可根據(jù)當(dāng)前的金屬價格和企業(yè)的生產(chǎn)成本，可以快捷準(zhǔn)確地確定礦體的經(jīng)濟開采品位和噸位，為開采設(shè)計方案的優(yōu)化提供可靠的依據(jù)。

表3 平行斷面法、克立格法儲量計算結(jié)果對比表

隨著地質(zhì)統(tǒng)計學(xué)方法在礦產(chǎn)儲量計算過程中的應(yīng)用與發(fā)展，傳統(tǒng)的儲量計算方法(如地質(zhì)塊段法、平行斷面法等)在計算較復(fù)雜的礦床時，受到了巨大挑戰(zhàn)。通過本文的儲量計算過程可以看出，基于Surpac 軟件的地質(zhì)統(tǒng)計學(xué)方法可以十分快捷地計算礦床儲量，并提交儲量報告。同時，可以根據(jù)不同需求及時繪制出新的指標(biāo)下礦體的形態(tài)，并動態(tài)計算任意邊界的礦產(chǎn)儲量，大大節(jié)約了人力、物力和時間，提高了效率和精度。但是如果按照傳統(tǒng)的方法，在改算礦體的儲量時，必須重新圈定礦體和確定計算邊界，這時工作量是非常大的，效率很低。因此應(yīng)用地質(zhì)統(tǒng)計學(xué)方法進行礦床儲量計算具有十分明顯的優(yōu)勢。

6 結(jié)論

(1)筆者綜合運用地質(zhì)統(tǒng)計學(xué)理論，建立了三道莊鉬礦床數(shù)學(xué)模型，通過對鉬品位變異函數(shù)分析，發(fā)現(xiàn)厚度、傾向、走向3 個方向的曲線波動性均較小，且均有較大的變程，說明三道莊鉬礦體鉬品位空間分布具有較好的連續(xù)性，礦體內(nèi)部品位分布較均勻。

(2)通過對比分析普通克立格法和傳統(tǒng)平行斷面法計算礦體的資源量，結(jié)果表明:所建立的三道莊鉬礦床數(shù)學(xué)模型可靠，所選用的區(qū)域化變量和模型估值方法合理，所采用的變異函數(shù)模型及其參數(shù)計算正確，可用于輔助礦山進行資源評估、采礦設(shè)計以及計劃編制等工作。

(3)地質(zhì)統(tǒng)計學(xué)方法與傳統(tǒng)的儲量計算方法相比具有明顯的優(yōu)點。它能最大限度地利用現(xiàn)有的勘探工程所提供的資料，使計算的礦石品位和礦石儲量更為精準(zhǔn)和快捷。

［1］侯景儒.實用地質(zhì)統(tǒng)計學(xué)［M］.北京:地質(zhì)出版社，1998:2-6.

［2］侯景儒，黃競先.地質(zhì)統(tǒng)計學(xué)及其在礦產(chǎn)儲量計算中的應(yīng)用［M］.北京:地質(zhì)出版社，1982:31-46.

［3］甕紀(jì)昌，張云政，黃超勇，等.欒川三道莊特大型鉬鎢礦床地質(zhì)特征及礦床成因［J］.地質(zhì)與勘探，2010，46(1):41-48.

［4］汪江河.河南三道莊鎢鉬礦區(qū)深部下層礦的找礦意義［J］.地質(zhì)調(diào)查與研究，2009，32(4):292-296.

［5］鄭文寶，黎楓佶，唐菊興，等.基于Micromine 軟件下地質(zhì)統(tǒng)計學(xué)在甲瑪矽卡巖型銅多金屬礦儲量計算中的應(yīng)用［J］.地質(zhì)與勘探，2011，47(4):726-736.

［6］孫洪泉.地質(zhì)統(tǒng)計學(xué)及其應(yīng)用［M］.北京:中國礦業(yè)大學(xué)出版社，1990:59.

［7］Surpac software 國際軟件公司.Surpac Vision 地質(zhì)統(tǒng)計學(xué)培訓(xùn)資料［R］.北京:2012.

［8］羅周全，劉曉明，吳亞斌，等.地質(zhì)統(tǒng)計學(xué)在多金屬礦床儲量計算中的應(yīng)用研究［J］.地質(zhì)與勘探，2007，43(3):83-87.

［9］余海軍，李文昌，尹光候，等.普朗銅礦床銅品位分布地質(zhì)統(tǒng)計學(xué)研究［J］.地質(zhì)與勘探，2009，45(4):437-443.

［10］陳愛兵，秦德先，張學(xué)書，等.基于MICROMINE 礦床三維立體模型的應(yīng)用［J］.地質(zhì)與勘探，2004，40(5):77-80.

［11］賈明濤，潘長良，王李管.克服地質(zhì)統(tǒng)計學(xué)礦床建模中主觀因素影響技術(shù)研究［J］.地質(zhì)與勘探，2003，39(4):73-77.

［12］余海軍，李文昌，曾普勝，等.地質(zhì)統(tǒng)計學(xué)在羊拉銅礦儲量計算中的應(yīng)用［J］.地球?qū)W報，2009，30(5):684-690.

［13］徐 靜，胡乃聯(lián).Surpac 軟件在某金銅礦山的應(yīng)用研究［J］.黃金科學(xué)技術(shù)，2007，15(1):54-58.

［14］馮超東，楊 鵬，胡乃聯(lián).克立格法在SURPAC 軟件中的實現(xiàn)及應(yīng)用［J］.金屬礦山，2007，370(4):55-58.