榮建鋒，朱力

（1.中廣核鈾業(yè)發(fā)展有限公司，北京 100029；2.核工業(yè)北京地質(zhì)研究院，北京 100029）

白崗巖型鈾礦是納米比亞在世界上獨(dú)有的鈾礦類(lèi)型，納米比亞湖山白崗巖型鈾礦床（Husab Uranium Deposits）是繼羅辛鈾礦床之后在納米比亞發(fā)現(xiàn)的又一個(gè)世界級(jí)特大型鈾礦床，是21 世紀(jì)以來(lái)非洲乃至世界鈾資源勘查領(lǐng)域新發(fā)現(xiàn)的最大鈾礦，其資源量和產(chǎn)量均位居世界第3 位［1-2］，在全球占有重要位置。

將三維建模技術(shù)與地質(zhì)統(tǒng)計(jì)學(xué)方法結(jié)合起來(lái)估算資源儲(chǔ)量是當(dāng)前礦產(chǎn)資源儲(chǔ)量估算的發(fā)展方向［3］，但在國(guó)內(nèi)礦產(chǎn)勘查中應(yīng)用還不普及。由于湖山鈾礦床資源量巨大，礦量集中，礦體主要呈似層狀、脈狀、透鏡狀，礦體數(shù)量多，利用傳統(tǒng)的塊段法估算資源礦體連接、圈定和計(jì)算繁雜，工作量巨大，周期長(zhǎng)易出差錯(cuò)，利用地質(zhì)統(tǒng)計(jì)學(xué)和專(zhuān)業(yè)礦業(yè)軟件估算資源勢(shì)在必行。2016 年以前，湖山鈾礦床的資源量估算和更新均委托國(guó)外公司開(kāi)展。2019 年，中廣核鈾業(yè)公司在完成湖山鈾礦1、2 號(hào)礦帶外圍補(bǔ)充勘探的基礎(chǔ)上，應(yīng)用地質(zhì)統(tǒng)計(jì)學(xué)和專(zhuān)業(yè)礦業(yè)軟件Data Mine，首次自主開(kāi)展資源估算，以驗(yàn)證湖山鈾礦1、2 號(hào)礦體的資源量，并積累經(jīng)驗(yàn)、培養(yǎng)隊(duì)伍。本文以湖山鈾礦北部的1號(hào)礦帶（簡(jiǎn)稱(chēng)Z1）為例探討地質(zhì)統(tǒng)計(jì)法在資源量估算中的應(yīng)用。

1 礦床地質(zhì)特征和地質(zhì)模型

湖山鈾礦床發(fā)現(xiàn)于2006 年，通過(guò)約5×105m的鉆探施工，控制了一條由1～5 號(hào)礦帶組成、南北走向上長(zhǎng)約8 km 的特大型白崗巖型鈾礦帶，其中Z1 礦體規(guī)模最大，資源量占整個(gè)礦床的近一半（圖1）。湖山白崗巖型鈾礦床大地構(gòu)造上位于納米比亞埃龍戈地區(qū)鈾礦省的泛非達(dá)馬拉造山帶的南部中央?yún)^(qū)帶，礦區(qū)出露主要地層自下而上有：前達(dá)馬拉基底的阿巴比斯組雜巖體（Abbabis）、達(dá)馬拉層序的艾杜西斯組（Etusis）、可汗組（Khan）、羅辛組（Rossing）、楚斯組（Chuos）、阿蘭蒂斯組（Arandis）、卡里比布組（Karibib）和卡塞布組（Kuiseb）。含礦主巖為淺灰白色白崗巖，礦化巖體主要侵入羅辛組，極少量侵入可汗組內(nèi)，礦體主要呈似層狀、脈狀、透鏡狀。

圖1 湖山鈾礦地層和構(gòu)造示意圖Fig.1 Schematic map of strata and structure in Husab uranium deposit

Z1 礦體主要侵入于湖山北部穹窿的東翼，走向28°，西翼走向355°，傾角為20°～55°。礦體主要產(chǎn)于羅辛組內(nèi)，礦化白崗巖產(chǎn)狀與地層一致，受地層控制明顯，呈似層狀、脈狀和透鏡狀，總體成群脈狀（圖2）。

圖2 Z1 礦體7505850 勘探線礦體解譯剖面圖Fig.2 Interpretation of Z1 ore body in section of exploration line 7505850

2 礦床數(shù)學(xué)模型的建立

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

本文以湖山鈾礦床Z1 的鈾礦體為資源量估算對(duì)象，區(qū)域化變量選取主要成礦元素U的品位。該變量來(lái)自于1 223 個(gè)鉆孔，共260 236 個(gè)巖性編錄數(shù)據(jù)，278 623 件樣品分析結(jié)果，在Data Mine 軟件中創(chuàng)建地質(zhì)數(shù)據(jù)庫(kù)，建立如下數(shù)據(jù)表。

孔位表：包括鉆孔的孔口數(shù)據(jù)、鉆探深度及鉆孔類(lèi)型等基本信息；

測(cè)斜表：包括鉆探施工過(guò)程中的測(cè)斜信息；

巖性編錄表：包括巖心/巖屑的巖性記錄、地層代號(hào)等信息；

樣品品位表：包括鉆孔取樣分析結(jié)果及其取樣位置等信息。

2.2 礦體實(shí)體模型

在建立的數(shù)據(jù)庫(kù)基礎(chǔ)上，首先建立地層模型，然后按照工業(yè)指標(biāo)將相鄰工程對(duì)應(yīng)的見(jiàn)礦部位用直線進(jìn)行連接，圈定出剖面礦體解譯線（圖3）。根據(jù)剖面上礦體解譯線，進(jìn)入三維界面將相鄰剖面對(duì)應(yīng)的線框進(jìn)行連接，在三維空間上形成礦體的實(shí)體模型［4］（圖4）。實(shí)體模型直觀地體現(xiàn)了礦體空間展布的幾何形態(tài)，其內(nèi)部包含的原始品位數(shù)據(jù)則是礦塊模型估值的重要基礎(chǔ)。

圖3 數(shù)據(jù)庫(kù)圖形化文件（鉆孔＋地層三維模型）Fig.3 Graphical database file（boreholes and strata model）

圖4 Z1 礦體三維空間模型示意圖Fig.4 Schematic 3D model of Z1 orebody

3 礦床數(shù)學(xué)模型的分析研究

3.1 特高品位識(shí)別和處理方法

特高品位的識(shí)別和處理是資源量估算中必須要考慮的步驟。傳統(tǒng)的識(shí)別方法有3 種：采用品位均值的6～8 倍值作為特異值下限；將累積概率達(dá)到97.5%的樣品品位值作為特異值下限［5］；侯景儒等曾提出根據(jù)樣品品位的變異（化）系數(shù)來(lái)識(shí)別特異值，在不同的變異系數(shù)范圍，選取不同倍數(shù)（2～20 倍）的均值作為特異值下限［6］。國(guó)際咨詢公司地質(zhì)師在進(jìn)行特高樣品的處理時(shí)，通常依據(jù)變化系數(shù)和西舍爾值來(lái)確定累積概率的分位值，一般選擇礦體內(nèi)品位數(shù)據(jù)的97.5%或99%分位數(shù)的數(shù)值作為特高品位的下限值。

經(jīng)統(tǒng)計(jì)，在Z1 礦體中，樣品品位最大值為243 618×10-6，平均品位為458×10-6，標(biāo)準(zhǔn)方差為1 654.3，品位變化系數(shù)為3.612，鈾品位在礦體域中呈偏態(tài)分布，需要進(jìn)行特高品位的處理（圖5）。

圖5 Z1 礦體原始樣品鈾品位統(tǒng)計(jì)直方圖和累積分布曲線Fig.5 Statistical histogram and cumulative distribution curve of uranium grade in Z1 orebody

結(jié)合累積概率曲線來(lái)看，鈾品位大致在99%的分位數(shù)附近品位分布明顯少，在6 000×10-6處對(duì)數(shù)累積頻率曲線有明顯的轉(zhuǎn)折（圖6）。按6 000×10-6下限值處理樣品后，品位均值變化不大，品位變化系數(shù)從3.61 降至1.86，西舍爾值與平均值比較接近，標(biāo)準(zhǔn)差、偏度、峰度值也有大幅改善，故將6 000×10-6定為Z1 特異值的下限（表1）。

表1 Z1 礦體特異值處理前后礦化域內(nèi)樣品統(tǒng)計(jì)數(shù)據(jù)對(duì)比Table 1 Statistic of the samples and grade before and after outliers treatment in Z1 orebody

圖6 U3O8品位累積頻率曲線及截止品位Fig.6 The cumulative frequency curve of U3O8 and the top-cut grade

3.2 樣品組合與統(tǒng)計(jì)分析

樣品組合是將不等長(zhǎng)的樣品長(zhǎng)度和品位統(tǒng)一到等長(zhǎng)的樣品上。根據(jù)地質(zhì)統(tǒng)計(jì)學(xué)原理，為確保得到各參數(shù)的最優(yōu)無(wú)偏估計(jì)量，所有的樣品數(shù)據(jù)應(yīng)落在相同的承載上，即同一類(lèi)參數(shù)的地質(zhì)樣品段的長(zhǎng)度應(yīng)該一致［7-9］。因此，在估算之前，必須按照一定長(zhǎng)度對(duì)長(zhǎng)度不等的樣品數(shù)據(jù)進(jìn)行組合，并統(tǒng)計(jì)分析其分布特征。樣品組合有多種方式，如按鉆孔長(zhǎng)度組合、按樣品長(zhǎng)度組合、按臺(tái)階高度組合、混合組合等［10］。經(jīng)統(tǒng)計(jì)，Z1 礦體內(nèi)含有樣品72 053 件，平均樣長(zhǎng)0.82 m，1 m 樣長(zhǎng)占比65%，處理特異值后按組合樣品長(zhǎng)度分別設(shè)為1 m、3 m 進(jìn)行統(tǒng)計(jì)和比較。兩者的統(tǒng)計(jì)學(xué)分布特征，兩種組合方案的樣品特征基本不變，但是3 m 組合方案的變化系數(shù)、偏度和峰度值均有所降低，說(shuō)明樣品數(shù)據(jù)更符合地質(zhì)統(tǒng)計(jì)學(xué)對(duì)參與估值數(shù)據(jù)的要求；且樣品總數(shù)大幅減少有助于節(jié)省軟件處理數(shù)據(jù)花費(fèi)的時(shí)間。綜上，判定3 m 樣品組合是較優(yōu)方案（表2）。

表2 Z1 礦體特異值處理后樣品組合統(tǒng)計(jì)數(shù)據(jù)對(duì)比Table 2 Statistic of the composite grade after outliers treatment in Z1 orebody

對(duì)3 m 樣品組合后的數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析，獲取U3O8樣品的品位分布直方圖呈正偏分布，對(duì)品位數(shù)據(jù)取對(duì)數(shù)后，樣品基本呈正態(tài)分布，無(wú)明顯的多峰分布特征，可以采用普通克里格進(jìn)行品位估值（圖7、8）。

圖7 Z1 礦體組合樣U3O8品位分布直方圖Fig.7 Histogram of U3O8 grades in a composite sample of Z1 orebody

3.3 變異函數(shù)分析

在進(jìn)行品位變異函數(shù)分析時(shí)，根據(jù)礦體產(chǎn)出特征，分別沿礦體走向（主軸）、傾向（次軸）、厚度（最小軸）3 個(gè)方向進(jìn)行變異函數(shù)分析，從而確定各方向?qū)嶒?yàn)變異函數(shù)計(jì)算時(shí)角度及距離參數(shù)（表3）［11］。

表3 實(shí)驗(yàn)變異函數(shù)參數(shù)表Table 3 Parameters of experimental variation functions

3.4 結(jié)構(gòu)分析

按照數(shù)據(jù)分析流程，需對(duì)區(qū)域化變量進(jìn)行實(shí)驗(yàn)變異函數(shù)計(jì)算，在此基礎(chǔ)上擬合出理論變異函數(shù)。借助結(jié)構(gòu)分析可以更全面的掌握區(qū)域化變量的結(jié)構(gòu)性和變異性，從而保證實(shí)驗(yàn)變異函數(shù)計(jì)算的精度。實(shí)驗(yàn)變異函數(shù)計(jì)算完成后，需選擇合適的數(shù)學(xué)模型表達(dá)變異函數(shù)的曲線特征。用一個(gè)球狀模型來(lái)表示各個(gè)方向的變異性，從而實(shí)現(xiàn)理論變異函數(shù)的擬合，這是業(yè)內(nèi)廣泛使用的方法［12］。

圖8 Z1 礦體組合樣U3O8品位對(duì)數(shù)值分布Fig.8 Logarithmic distribution of U3O8 grades in a composite sample of Z1 orebody

球狀變異函數(shù)模型標(biāo)準(zhǔn)形式為：

式中：γ（h）—變異函數(shù)；h—滯后距離；C0—塊金常數(shù)；C—躍遷常數(shù)；C0+C—基臺(tái)值；a—變程［13］。

根據(jù)礦體空間賦存特征，結(jié)合成礦規(guī)律研究成果，使用Data Mine 的分析模擬功能調(diào)算各種參數(shù)，最終得到礦體走向、傾向、厚度方向變異函數(shù)擬合曲線（圖9）。所擬合出的球狀模型理論變異函數(shù)最佳結(jié)構(gòu)參數(shù)見(jiàn)表4。

表4 Z1 礦體最佳理論變異函數(shù)結(jié)構(gòu)參數(shù)表Table 4 Main parameters of the best variogram of Z1 orebody

圖9 Z1 礦體三維方向變異函數(shù)曲線Fig.9 The variogram of grade in 3D of Z1 orebody

3.5 理論變異函數(shù)交叉驗(yàn)證

在進(jìn)行理論變異函數(shù)的擬合后，用交叉驗(yàn)證的方法對(duì)模型設(shè)置的參數(shù)進(jìn)行檢驗(yàn)。交叉驗(yàn)證的統(tǒng)計(jì)結(jié)果顯示，均值差較小接近0，且方差/標(biāo)準(zhǔn)差趨近于1，良好的交叉驗(yàn)證結(jié)果證明變異函數(shù)結(jié)構(gòu)模型可靠［14］，可作為普通克里格估值的參數(shù)（表5）。因此，前面分析得到的理論變異函數(shù)參數(shù)是可靠的，可用于該礦床U3O8品位估算。

表5 各項(xiàng)異性變異函數(shù)參數(shù)的交叉驗(yàn)證結(jié)果Table 5 Results of cross validation of each hetero variogram parameters

4 資源量估算

4.1 礦塊模型的建立

在確定單元塊尺寸時(shí)主要考慮采礦方法、最小開(kāi)采單元大小、礦體形態(tài)、勘探線網(wǎng)度、最小可采厚度、夾石剔除厚度以及運(yùn)算效率等因素［15］。該礦床采用露天開(kāi)采方式，最小采礦臺(tái)階高7.5 m，設(shè)定單元塊尺寸為12.5 m×12.5 m×7.5 m，在礦體邊緣部次級(jí)塊設(shè)為6.25 m×6.25 m×2.5 m。

4.2 礦塊模型估值

礦塊模型估值就是求出塊模型中每個(gè)塊中心位置的品位值，這就需要在具有相關(guān)性的空間范圍內(nèi)搜索已知點(diǎn)，并通過(guò)給予一定的權(quán)系數(shù)實(shí)現(xiàn)“最佳、線性、無(wú)偏”估計(jì)［16］。本次工作采用普通克里格法分三次進(jìn)行估值，完成估值后的塊體模型可以在三維空間詳細(xì)列示不同品位級(jí)別的礦石的分布情況（表6，圖10）。

圖10 Z1 礦體品位估值塊體模型圖Fig.10 Block model diagram of Z1 orebody

表6 Z1 礦體橢球體參數(shù)表Table 6 Parameters of search volume of Z1 orebody

4.3 可靠性驗(yàn)證

為驗(yàn)證估值品位與原始數(shù)據(jù)的一致性，開(kāi)展全局檢驗(yàn)統(tǒng)計(jì)分析兩組數(shù)據(jù)在三維空間分布特征是必要的。Swath 曲線就是用于刻畫(huà)在特定方向上樣品與估值模型品位的吻合情況。本 文從X 軸、Y 軸 和Z 軸3 個(gè)方 向，分別 對(duì) 比 樣品品位曲線和估值模型的曲線的品位變化趨勢(shì)（圖11）。整體而言，兩組品位曲線較為一致，只是在礦體兩端尖滅部位有較明顯的偏差，這是由于樣品數(shù)量較少，以及克里格估值方法在估值過(guò)程中產(chǎn)生的平滑效應(yīng)引起的。

圖11 Z1 礦體全局檢驗(yàn)Swath 曲線對(duì)比圖Fig.11 Swath plot of global inspection of Z1 orebody

同時(shí)，將礦體內(nèi)等長(zhǎng)組合樣品的品位均值與估算品位均值進(jìn)行比較，評(píng)估平均品位的變化，驗(yàn)證估值過(guò)程是否存在明顯的偏差，從而實(shí)現(xiàn)對(duì)估算結(jié)果的可靠性驗(yàn)證。統(tǒng)計(jì)結(jié)果表明，塊體模型中估值元素經(jīng)過(guò)插值，品位產(chǎn)生平滑效應(yīng)，平均品位有所降低（表7）。礦體等長(zhǎng)樣品位均值與塊體品位插值均值比為92%，說(shuō)明數(shù)據(jù)插值擬合度較高，沒(méi)有明顯偏差。

表7 礦體等長(zhǎng)樣品位統(tǒng)計(jì)與礦塊體模型品位統(tǒng)計(jì)對(duì)比表Table 7 Statistics of samples and interpolation grade

為保證估值質(zhì)量，本次工作還采用距離冪次方反比法對(duì)Z1 資源模型的工業(yè)礦部分進(jìn)行資源量估算結(jié)果驗(yàn)證。經(jīng)對(duì)比，兩種估算方法的礦石量相差1.5%，平均品位相差1.5%，金屬量相差0.2%。兩種估算方法的礦石量、金屬量及品位相差均在5%以內(nèi)，說(shuō)明本次研究采用Data Mine 軟件普通克里格法估算資源量是適當(dāng)?shù)?，估算結(jié)果是可靠的。

與2012 年資源量估算結(jié)果進(jìn)行對(duì)比，Z1 礦體的控制級(jí)以上資源量減少5%，推斷級(jí)資源量也有一定程度的減少，本次工作達(dá)到了資源量核實(shí)的目的。

5 結(jié)論

本次研究綜合運(yùn)用地質(zhì)統(tǒng)計(jì)學(xué)理論，在Data Mine 礦業(yè)軟件平臺(tái)上建立了納米比亞湖山鈾礦床的數(shù)學(xué)模型，在研究和處理特高品位數(shù)據(jù)的基礎(chǔ)上，使用球狀模型計(jì)算和擬合了不同方向?qū)嶒?yàn)變異函數(shù)并進(jìn)行了結(jié)構(gòu)分析，成功構(gòu)建出符合礦床地質(zhì)特征的結(jié)構(gòu)模型，交叉驗(yàn)證各結(jié)構(gòu)參數(shù)合理可靠；最后采用普通克立格法對(duì)鈾礦石品位及礦石量進(jìn)行了估算?？煽啃则?yàn)證及對(duì)比結(jié)果表明，本次研究所選用的區(qū)域化變量和模型估值方法合理，符合礦床空間展布特征和分布規(guī)律。地質(zhì)統(tǒng)計(jì)學(xué)理論和Data Mine 礦業(yè)軟件在資源量估算中快速、動(dòng)態(tài)和直觀的特性得到充分體現(xiàn)，其成果可用于礦山資源評(píng)估、露天采礦設(shè)計(jì)以及生產(chǎn)計(jì)劃編制等工作。