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

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

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

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

1 礦床地質(zhì)特征

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 礦化域塊體模型

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

4 特異值處理及合理性檢驗(yàn)

4.1 特異值的識(shí)別與處理

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

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

對(duì)礦化域樣品分別繪制原始品位與對(duì)數(shù)品位分布直方圖(圖4)。 由圖4(a)可知，在品位大于5 g/t之后樣品數(shù)量大幅減少。 由圖4(b)可知，品位對(duì)數(shù)值在0.66附近出現(xiàn)不連續(xù)，此時(shí)樣品品位為4.57 g/t。因此，本文將4.57 g/t作為特高品位下限進(jìn)行替換處理，共替換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ī)律，則需將高品位樣帶單獨(dú)圈出進(jìn)行資源量估算，而不作為特高品位樣品處理[16]。經(jīng)統(tǒng)計(jì)與觀察，12件特高品位樣品三維空間分布零散，不具備單獨(dú)劃分高品位礦體域的條件。

4.3 特異值處理后的檢查

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

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

5.1 樣品組合

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

圖5 礦化域樣品樣長(zhǎng)直方圖Fig.5 Sample length histogram of mineralized domain samples

5.2 變異函數(shù)擬合

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

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

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

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

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

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

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

對(duì)比三套理論變異函數(shù)交叉驗(yàn)證結(jié)果(圖6)，采用參數(shù)三時(shí)，具有最趨近于“0”的誤差均值，以及最趨近于“1”的誤差方差與克里格估計(jì)方差比值。因此，參數(shù)三是最為合理的理論變異函數(shù)參數(shù)，參數(shù)三的主軸、次軸及短軸方向變異函數(shù)見(jiàn)圖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 資源量估值、驗(yàn)證及分類

6.1 普通克里格估值

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

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

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

6.2 估值結(jié)果驗(yàn)證

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

6.3 資源量分類

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

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

圖8 全局驗(yàn)證示意圖Fig.8 Schematic diagram of global verification

圖9 全局驗(yàn)證曲線圖Fig.9 Global verification diagram

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

7 結(jié) 論

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

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

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

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