杜菊民，祖輔平，景永波，陳 誠，陳 峰，蔡衛(wèi)東

(1.江蘇省地質(zhì)工程有限公司，江蘇 南京 210018； 2.西安石油大學地球科學與工程學院，陜西 西安 710065； 3.江蘇省地質(zhì)勘查技術(shù)院，江蘇 南京 210049； 4.江蘇省礦產(chǎn)勘查局第一地質(zhì)大隊，江蘇 南京 210018)

隨著地質(zhì)統(tǒng)計學與計算機科學的發(fā)展，三維地質(zhì)建模已在資源量估算、成礦預(yù)測等領(lǐng)域得到廣泛應(yīng)用[1-3]。國內(nèi)受到傳統(tǒng)地質(zhì)勘查規(guī)范的約束[4]，通常采用邊界品位、工業(yè)品位的雙指標來建立礦體三維模型，并進行資源量估算[5-6]。以工業(yè)雙指標體系圈定礦體的方式屬于硬邊界處理，割裂了地質(zhì)數(shù)據(jù)的空間相關(guān)性[7]，不僅影響了三維資源量估算工作的效率及礦床成礦預(yù)測的準確性[8-10]，還阻礙了動態(tài)資源量估算，掩蓋了三維軟件估值的優(yōu)勢[11]。目前國際社會普遍采用低于邊界品位的礦化域指標進行資源量估算，即軟邊界處理方式[7]。實踐證明，采用礦化域指標對斑巖型、熱液蝕變型等礦床進行資源量估算，不僅在邊界品位動態(tài)優(yōu)化、工業(yè)指標論證時，能大幅提高工作效率[9,12]，而且更能客觀反映礦化規(guī)律，有利于進行成礦預(yù)測[8,13]。隨著《資源量估算規(guī)程》(DZ/T 0338—2020)的發(fā)布，采用礦化域模型進行資源量估算必將得到更廣泛的應(yīng)用。

澳大利亞某金礦成礦類型為與侵入巖相關(guān)型，總體品位偏低,礦巖界線呈過渡關(guān)系，因此，可先用相對低的品位圈連出完整性較好的礦化域，再采用礦塊指標體系圈定礦體范圍。本次研究在地質(zhì)統(tǒng)計的基礎(chǔ)上，對礦化域品位邊界值、特異值、理論變異函數(shù)、搜索橢球體等關(guān)鍵參數(shù)進行了優(yōu)化設(shè)置，以普通克里格法進行了資源量估算，并對估值結(jié)果進行全局驗證。本文對相關(guān)資源量估算工作具有一定參考作用。

1 礦床地質(zhì)特征

澳大利亞拉克蘭造山帶(Lachlan Orogen)為典型的大陸邊緣增生造山帶，由西帶、中帶和東帶三個亞帶組成[14]。本文研究的金礦位于拉克蘭造山帶中帶，礦體賦存于晚古生代細粒二長閃長巖株中，該巖株侵入至早古生代強烈褶皺變形的石英片巖。主要蝕變類型包括：硅化、絹云母化、鈉長石化、黃鐵礦化、毒砂化等，圍繞細粒二長閃長巖株蝕變暈分帶明顯，在該巖株外側(cè)還發(fā)育局部角巖化蝕變帶。礦體呈筒狀，直徑約160 m，垂直往下延伸近1 000 m。礦化以金為主，砷含量較高，為低品位大型金礦。

2 數(shù)據(jù)庫建立與樣品基本統(tǒng)計

2.1 鉆孔數(shù)據(jù)庫的建立

鉆孔數(shù)據(jù)庫是進行三維地質(zhì)解譯、品位估值的基礎(chǔ)。本文基于3DMine三維礦業(yè)工程軟件平臺技術(shù)要求，提取鉆孔定位數(shù)據(jù)表、測斜數(shù)據(jù)表、樣品分析數(shù)據(jù)表、地質(zhì)巖性表等，建立鉆孔數(shù)據(jù)庫。經(jīng)數(shù)據(jù)糾錯、完整性和邏輯性檢查后，用于后期樣品統(tǒng)計和地質(zhì)建模工作。鉆孔數(shù)據(jù)根據(jù)需要進行顯示，以方便后續(xù)地質(zhì)解譯(圖1)。礦體由54個鉆孔控制，孔深18.0～1 312.6 m，總進尺12 091 m，采樣6 364件，樣長1～3 m。

圖1 澳大利亞某金礦地表及鉆孔三維效果圖Fig.1 3D renderings of surface and boreholes of a gold mine in Australia

2.2 原始樣品基本統(tǒng)計

提取所有原始樣品點進行基本分析。從原始樣品品位分布直方圖及分布曲線來看，樣品呈正偏態(tài)分布(圖2(a))，經(jīng)對數(shù)處理后，礦化部分樣品服從正態(tài)分布，因此可以采用地質(zhì)統(tǒng)計法進行資源量估算。 原始樣品中大于金檢出限樣品共5 808件，最小值0.01 g/t，最大值59.4 g/t。 其中，平均值(0.806 98 g/t)與西舍爾估值(1.282 06 g/t)偏差過大，偏度(25.267 01)也較大(表1)，說明需進行特異值處理。

圖2 原始樣品品位分布直方圖Fig.2 Histogram of original sample grade distribution

表1 不同類型樣品基本統(tǒng)計信息Table 1 Basic statistical information of different types of samples

2.3 礦化域品位邊界值的確定

原始樣品品位值對數(shù)處理后，品位對數(shù)值-0.6～-0.8之間的樣品數(shù)量分布存在明顯拐點(圖2(b))。對品位對數(shù)值進行分布區(qū)間比例統(tǒng)計，結(jié)果見表2。由表2可知，在對數(shù)值為-0.744 7時，樣品占比為30.0%，對應(yīng)原始樣品品位為0.19 g/t，即大于等于0.19 g/t的樣品占樣品總數(shù)的70.0%，代表了拐點后的樣品礦化較為連續(xù)。因此，本文采用0.19 g/t作為礦化域品位邊界值。

3 礦化域模型的構(gòu)建

3.1 礦化域?qū)嶓w模型

實體模型又稱礦體的表面模型或線框模型，由封閉的三角網(wǎng)組成，可以進行體積估值、三維顯示等[15]。先根據(jù)礦化域品位邊界值圈定大致連續(xù)的礦化范圍，再采用樣條曲線對礦化域邊界進行平滑與漸變處理[11,13]，使實體模型更接近自然實際情況。建立的礦化域模型見圖3(a)。

表2 原始樣品及礦化域樣品品位分布Table 2 Grade distribution of original samples and mineralized domain samples

3.2 礦化域塊體模型

實體模型內(nèi)部為空值，無法進行空間品位、體重、礦巖類型等屬性賦值或估值，需按礦化域走向、傾向等因素，確定基本塊體規(guī)格，創(chuàng)建塊體模型。礦化域?qū)嶓w邊緣設(shè)置次級塊體，以精細模擬礦化域邊界。基本塊體規(guī)格一般設(shè)為工程間距的1/4～1/5，次級礦體規(guī)格一般設(shè)為基本塊體的1/2或1/4?；緣K體太大無法顯示品位變化，塊太小可導致品位估值不可靠。本次工作根據(jù)礦化域模型形態(tài)及鉆孔間距，將礦化域基本塊體規(guī)格設(shè)置為4 m×4 m×10 m，次級塊體設(shè)置為1.0 m×1.0 m×2.5 m。礦化域塊體模型見圖3(b)。

4 特異值處理及合理性檢驗

4.1 特異值的識別與處理

特異值是指特高品位和特低品位，在品位估值中通常指特高品位。由于其超過了期望的品位分布范圍，對品位均值、方差等造成較大影響，會引起品位高估[16]。采用普通克里格法進行估值，進行變異函數(shù)分析及特異值處理時，應(yīng)以礦化域或礦體為單位進行。特異值的識別與處理有多種方法，本文選擇分布衰減或不連續(xù)的點作為特高品位。

利用礦化域?qū)嶓w模型提取礦化域樣品進行基本分析，共提取4 052件樣品，最小值0.19 g/t，最大值59.4 g/t。其中平均值(1.130 05 g/t)與西舍爾估值(1.135 87 g/t)接近，但偏度(27.556 02)較大(表1)，說明存在特高品位。

對礦化域樣品分別繪制原始品位與對數(shù)品位分布直方圖(圖4)。 由圖4(a)可知，在品位大于5 g/t之后樣品數(shù)量大幅減少。 由圖4(b)可知，品位對數(shù)值在0.66附近出現(xiàn)不連續(xù)，此時樣品品位為4.57 g/t。因此，本文將4.57 g/t作為特高品位下限進行替換處理，共替換12件樣品。

圖3 礦化域?qū)嶓w及塊體模型(局部)示意圖Fig.3 Part of 3D mineralized domain entity model and block model of gold deposit

圖4 礦化域樣品品位分布直方圖Fig.4 Histogram of sample grade distribution in mineralized domain

4.2 特異值的空間分布

確定特高品位值后，需觀察特高品位樣品在空間的位置，如空間距離很近且分布有規(guī)律，則需將高品位樣帶單獨圈出進行資源量估算，而不作為特高品位樣品處理[16]。經(jīng)統(tǒng)計與觀察，12件特高品位樣品三維空間分布零散，不具備單獨劃分高品位礦體域的條件。

4.3 特異值處理后的檢查

采用西舍爾估值進行合理性檢驗，處理特異值后樣品的算術(shù)平均值小于且與西舍爾估值接近時，表明特異值處理合理。 對處理完特異值的樣品數(shù)據(jù)進行基本統(tǒng)計，算術(shù)平均值(1.101 92 g/t)小于且與西舍爾估值(1.127 36 g/t)接近，偏度(1.357 51)相對特異值處理前大幅減少，變異系數(shù)(0.611 31)進一步減少(表1)，說明處理后樣品接近正態(tài)分布，礦化相對均勻。因此，以品位4.57 g/t為下限進行特異值處理取值合理。

5 樣品組合與變異函數(shù)分析

5.1 樣品組合

樣品組合是將空間上長度不等的樣品量化到等長的離散點，以保證參數(shù)統(tǒng)計的無偏估計[6,15,17]。 對處理特異值后礦化域樣品進行樣長統(tǒng)計，樣長集中在1 m、2 m和3 m，平均1.7 m(圖5)，所以本次以2 m作為樣品組合長度，最小樣品長度選擇1 m。 樣品等長組合后，共生成3 435件組合樣。與處理特異值后礦化域樣品相比，組合樣品的平均值(1.093 62 g/t)與西舍爾估值(1.118 07 g/t)更為接近(表1)。

圖5 礦化域樣品樣長直方圖Fig.5 Sample length histogram of mineralized domain samples

5.2 變異函數(shù)擬合

變異函數(shù)是地質(zhì)統(tǒng)計學中研究樣品空間分布隨機性和結(jié)構(gòu)性的主要工具。使用塊金值(C0)代表樣品變化的隨機成分，基臺值(C)代表樣品變化的相關(guān)成分，變程(R)代表相關(guān)性存在的最大范圍。樣品的相關(guān)性與礦種、成礦類型、樣品采集等因素有關(guān)。

對組合后的樣品進行變異函數(shù)分析，根據(jù)樣品空間分布位置，設(shè)置16個扇區(qū)，進行主軸變異函數(shù)分析，依次確定主軸、次軸和短軸的實驗變異函數(shù)；采用球狀模型理論變異函數(shù)，通過調(diào)整塊金值、基臺與變程對實驗變異函數(shù)進行擬合。本次工作擬合了三套參數(shù)進行交叉驗證，從中優(yōu)選出最合理的變異函數(shù)模型參數(shù)，各套模型參數(shù)詳見表3。

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

使用構(gòu)建的變異函數(shù)對已知樣品點進行估值，比較估計值和真值的差值，并進行差值統(tǒng)計分析，交叉驗證塊金值、基臺和變程的合理性。真值與估計值的誤差均值趨近于“0”，誤差方差與克里格估計方差比值趨近于“1”，為交叉驗證的理想結(jié)果[15]。

基于以上思路，分別對三套理論變異函數(shù)參數(shù)，設(shè)置主軸搜索半徑范圍在10 m至200 m之間，開展交叉驗證對比， 驗證結(jié)果生成折線圖(圖6)。 隨著主軸搜索半徑的增大，結(jié)果顯示如下幾方面特征：①三套參數(shù)的誤差均值快速增加，并趨于平穩(wěn)；②誤差方差與克里格估計方差比值增加至峰值后，緩慢回落；③參數(shù)一在主軸搜索半徑設(shè)置為50 m時，真值與估計值的誤差均值為0.010 1，誤差方差與克里格估計方差比值為0.927 2；④參數(shù)二在主軸搜索半徑設(shè)置為80 m時，誤差均值為0.010 2，方差比值為0.895 6；⑤參數(shù)三在主軸搜索半徑設(shè)置為70 m時，誤差均值為0.009 7，方差比值為0.999 6。

表3 理論變異函數(shù)參數(shù)組合Table 3 Combinations of theoretical variogram parameters

圖6 不同參數(shù)交叉驗證結(jié)果統(tǒng)計Fig.6 Statistics of cross-validation results of different parameters

對比三套理論變異函數(shù)交叉驗證結(jié)果(圖6)，采用參數(shù)三時，具有最趨近于“0”的誤差均值，以及最趨近于“1”的誤差方差與克里格估計方差比值。因此，參數(shù)三是最為合理的理論變異函數(shù)參數(shù)，參數(shù)三的主軸、次軸及短軸方向變異函數(shù)見圖7。

圖7 參數(shù)三主軸、次軸及短軸方向變異函數(shù)示意圖Fig.7 Schematic diagram of the variation function from major axis,minor axis and minor axis of the parameter three

6 資源量估值、驗證及分類

6.1 普通克里格估值

根據(jù)礦體特征、工程分布及變異函數(shù)(參數(shù)三)，確定本次普通克里格估值參數(shù)(表4)，對塊體模型的Au品位屬性進行估值。主軸搜索半徑參考勘探工程間距及變程，設(shè)置為25 m及其偶數(shù)倍。最少樣品數(shù)、最多樣品數(shù)及每孔最多選擇樣品數(shù)的設(shè)置可以減弱樣品數(shù)據(jù)的叢聚效應(yīng)，提高估值精度。

表4 普通克里格法估值參數(shù)Table 4 Estimation parameters of ordinary Kringing method

同時，將到樣品點最近距離、到樣品點平均距離、樣品點數(shù)目、克里格方差、克里格效率、負權(quán)重數(shù)目、拉格郞日因子等信息寫到塊體中，以便對估值后的塊體進行統(tǒng)計分類。

6.2 估值結(jié)果驗證

采用全局驗證的方法對估值結(jié)果進行檢驗。在平面圖上按20 m的間距布設(shè)9條南北向勘探線，共劃分9個區(qū)域(圖8)。分別計算每一個區(qū)域內(nèi)參與估值的組合樣品平均值與塊體模型估值的平均值，結(jié)果繪制成兩條曲線圖，其下方為對應(yīng)的樣品數(shù)(圖9)。由圖9可知，塊體估值結(jié)果與實際樣品點平均品位較為接近，估值曲線擬合了實際品位變化，整體而言，兩條曲線吻合度較高，說明本次估值較為合理可靠。在5號勘探線、7號勘探線、8號勘探線，估值品位比樣點平均品位低0.1 g/t，可以發(fā)現(xiàn)這三個區(qū)域樣品點數(shù)相對較少，但同時估值的塊體較多，估值距離增大，是造成估值結(jié)果偏差增加的原因。

6.3 資源量分類

地質(zhì)統(tǒng)計法估算的資源量，可根據(jù)塊體屬性中的最近距離、平均距離、工程數(shù)、樣品數(shù)等記錄內(nèi)容，或是搜索半徑等進行資源量分類[18]。本文依據(jù)塊體估值獲取得到樣品點的平均距離，按0～50 m及50 m以上，分別對相應(yīng)塊體屬性設(shè)置為控制和推斷類別。

礦化域模型不僅估算了一般工業(yè)品位以上的資源量，還估算了礦化域品位邊界值至一般工業(yè)品位區(qū)間的礦石量、平均品位、金屬量等信息，并可按品位變化進行統(tǒng)計(圖10)，因此，在工業(yè)指標論證[9]、礦化特征分析[19]、成礦預(yù)測[8]等方面具有優(yōu)勢。

圖8 全局驗證示意圖Fig.8 Schematic diagram of global verification

圖9 全局驗證曲線圖Fig.9 Global verification diagram

圖10 礦化域品位-資源量統(tǒng)計圖Fig.10 Combination diagram of mineralized domain grade and resource changes

7 結(jié) 論

本文基于三維地質(zhì)軟件平臺及礦化域估算理論，通過對澳大利亞某大型金礦樣品的統(tǒng)計與分析，確立了礦化域品位邊界值，建立了三維礦化域?qū)嶓w與塊體模型；進一步對礦化域內(nèi)樣品進行地質(zhì)統(tǒng)計分析，進行特異值處理及變異函數(shù)優(yōu)化，最終以普通克里格法進行資源量估算。通過上述研究，主要得出以下結(jié)論。

1) 礦化域模型進行資源量估算，不僅估算了邊界品位以上礦體，對低于邊界品位的礦化部分也進行了估值，在品位指標論證時更為便捷高效。

2) 礦化域模型可以最大限度保留礦化的連續(xù)性，反映礦化的自然規(guī)律，展現(xiàn)礦化體的真實形態(tài)，有利于后期成礦預(yù)測。

3) 基于地質(zhì)統(tǒng)計理論，對礦化域品位邊界值、特異值、變異函數(shù)等參數(shù)進行優(yōu)化，是保證估算質(zhì)量的重要步驟。