劉佶林，楊 忠，王孝東，馮光華

(1. 云南華聯(lián)鋅銦股份有限公司，云南 文山663701；2. 昆明理工大學(xué) 國土資源工程學(xué)院，昆明650093)

目前，我國礦山在進(jìn)行資源儲量估算時，普遍采用塊段法和剖面法等傳統(tǒng)儲量估算方法[1]。數(shù)字礦山技術(shù)的運用，首次實現(xiàn)了對開采礦體的規(guī)模、形態(tài)、品位等屬性在空間分布上的數(shù)字化整體描述。近年來，諸如3DMine、DIMine、Surpac等三維礦業(yè)軟件的推廣與應(yīng)用，使得利用地質(zhì)統(tǒng)計學(xué)方法進(jìn)行資源儲量估算已成為國內(nèi)外資源評估的重要手段。

云南都龍銅街曼家寨錫鋅多金屬礦區(qū)(以下簡稱銅曼礦區(qū))位于云南省東南部，礦區(qū)面積約5 km2。礦區(qū)北部和東部以F0斷層為界，寬約1.5 km。礦區(qū)南面延伸至中越國境線，長約8 km，構(gòu)成南北向展布的錫、鋅、銅、銀等多金屬礦帶[2]。礦體走向近南北，向西傾斜，傾角為10°～40°，局部傾角高達(dá)60°。礦區(qū)具有疊瓦狀排列、分支復(fù)合、尖滅再現(xiàn)的特點，隨含礦層同步褶曲。銅曼露天采場內(nèi)地形地質(zhì)情況復(fù)雜，礦巖交錯分布，小礦體數(shù)量較多，礦體邊緣部夾石和夾礦現(xiàn)象明顯，如何充分利用大量的原始勘探數(shù)據(jù)和礦山生產(chǎn)的地質(zhì)編錄數(shù)據(jù)進(jìn)行儲量估算對礦區(qū)的礦產(chǎn)資源評估與礦山開采方案制定具有重要的意義。

1 三維地質(zhì)模型構(gòu)建

1.1 地質(zhì)數(shù)據(jù)庫

地質(zhì)數(shù)據(jù)庫是進(jìn)行地質(zhì)解譯、塊體品位估值、儲量估算等工作的重要基礎(chǔ)[3]，本次研究基于3DMine三維礦業(yè)軟件，充分利用各勘探時期的勘探成果和生產(chǎn)勘探資料，共計738個鉆孔和78條探槽數(shù)據(jù)，建立了三維地質(zhì)數(shù)據(jù)庫，如圖1所示。三維地質(zhì)數(shù)據(jù)庫的建立可將數(shù)字形式的勘探資料用三維圖形形象化、具體化，便于管理和分析利用。在地質(zhì)數(shù)據(jù)庫中可以用三維顯示方式瀏覽所有鉆孔的基本信息，顯示單個或多個工程的地質(zhì)品位、深度、軌跡等數(shù)據(jù)信息，還可根據(jù)需要設(shè)置不同的顯示風(fēng)格來查看鉆孔的空間分布情況。

圖1 三維地質(zhì)數(shù)據(jù)庫Fig.1 3D geologic database

1.2 礦體模型

礦體(實體)模型，是在三維空間內(nèi)由一系列剖面或空間點構(gòu)成的三角網(wǎng)包裹成封閉的實體，最直接的作用就是模擬礦體形態(tài)。本次研究利用3DMine三維礦業(yè)軟件，以銅曼礦區(qū)三維地質(zhì)數(shù)據(jù)庫(如圖1所示)為基礎(chǔ)，按照礦體圈定原則對礦區(qū)進(jìn)行了地質(zhì)解譯，圈連礦體(實體)模型四百余個(如圖2所示)。

圖2 礦體模型Fig.2 The ore body model

1.3 塊體模型

利用上述方法建立的礦體模型，礦體內(nèi)部是空的，沒有任何信息。為了便于后期儲量計算、境界設(shè)計等工作，需要在不規(guī)則的礦體內(nèi)部以及周邊小部分范圍內(nèi)填充規(guī)則的三維等塊狀模型(塊體)，這種塊體集稱為塊體模型。塊體模型是品位估值和儲量估算的基礎(chǔ)，也是大多數(shù)數(shù)學(xué)優(yōu)化方法的基礎(chǔ)。目前的露天開采規(guī)劃優(yōu)化方法、儲量計算、境界設(shè)計、采剝計劃編制等，幾乎都是以塊體模型作為研究手段。塊體模型可根據(jù)生產(chǎn)需要添加諸多屬性，如礦體編號、經(jīng)濟(jì)類型、礦石品位、礦石體重、礦巖類型等，便于在生產(chǎn)實踐中隨著已知信息量的增加或變化不斷更新屬性信息?？紤]到塊體邊界與礦體范圍的吻合度越高，塊體所反映的空間位置與礦體更趨近于一致，因此采用次級模塊的方法對原始塊體進(jìn)行分割，然后利用已知組合樣品點對整個礦體范圍內(nèi)的單元塊的品位進(jìn)行估計，并在此基礎(chǔ)上進(jìn)行儲量估算。

圖3 塊體模型Fig.3 The block model

2 儲量估算

一般來說，在三維礦業(yè)軟件中，利用三維地質(zhì)模型和克里格法進(jìn)行儲量估算的主要步驟[1]如圖4所示。

圖4 儲量估算基本流程Fig.4 The basic flow of reserve estimation

2.1 區(qū)域變化量選擇

應(yīng)用地質(zhì)統(tǒng)計學(xué)進(jìn)行儲量計算時，須根據(jù)礦床的具體情況和特點以及所采用的方法、手段來選取區(qū)域變化量[4]。由于銅曼礦區(qū)礦體形態(tài)復(fù)雜，礦巖交錯現(xiàn)象明顯，礦體形態(tài)、走向、傾向、厚度變化較大，故本次研究先根據(jù)勘探線剖面把礦體的實體模型確定下來，再利用克里格法對所圈定的礦體進(jìn)行品位估值，本次研究主要選擇品位值作為區(qū)域變化量。

2.2 基本統(tǒng)計

在對塊體模型進(jìn)行品位估值之前，要結(jié)合礦床的成因?qū)ΦV區(qū)樣品的品位分布特征進(jìn)行統(tǒng)計分析，為后續(xù)的品位估值提供數(shù)據(jù)基礎(chǔ)，以便根據(jù)礦床自身特點選擇適當(dāng)?shù)墓乐捣桨竅5]。利用3DMine軟件中的地質(zhì)統(tǒng)計學(xué)模塊對銅曼礦區(qū)的樣品點進(jìn)行統(tǒng)計分析，研究樣品點數(shù)據(jù)的分布特征，如數(shù)據(jù)不符合正態(tài)分布則進(jìn)行變換。本次研究主要分析銅曼礦區(qū)的Zn元素，統(tǒng)計結(jié)果如表1和表2所示。

表1 Zn元素樣品點分位數(shù)統(tǒng)計Table 1 Quantile statistics of Zn samples

表2 Zn元素樣品點基本統(tǒng)計Table 2 The basic statistics of Zn samples

統(tǒng)計分析的主要目的是確定礦區(qū)Zn元素樣品點的分布類型，為隨后的變異函數(shù)計算及Zn品位估值提供參考。

從圖5(a)可以看出，Zn元素樣品點不服從正態(tài)分布，故將Zn樣品點進(jìn)行對數(shù)轉(zhuǎn)換，使其服從正態(tài)分布,如圖5(b)所示。由于Zn元素樣品點存在一些特高品位，因此必須對特高品位進(jìn)行處理。本次特高品位處理參照礦區(qū)儲量核實報告，對超過礦床平均品位8倍的特高品位，采用工程平均品位代替。

2.3 樣品組合

地質(zhì)統(tǒng)計學(xué)要求參與估值計算的數(shù)據(jù)的支撐(指樣品的長度或體積)應(yīng)該一致[6]，因此樣品組合就是要將探礦工程中的樣長和品位值量化到離散點上，即每段樣長的中點，只有在工程方向上產(chǎn)生均勻(即等距離)的離散點才能用于資源儲量估算[7]。因此，樣品組合產(chǎn)生的離散點將用于塊體模型估值。

根據(jù)樣品的統(tǒng)計分析，其原始平均樣品長度為1.21，絕大部分樣品的樣長在1 m左右。因此，本次研究采用等距離為1 m的樣品長度進(jìn)行計算分析，最小組合樣長為平均樣長的50%，即0.5 m。

2.4 變異函數(shù)模型

由于變異函數(shù)計算直接影響到變異函數(shù)的擬合及克里格法估值的精度[8-9]，因此變異函數(shù)是克里格法儲量估算的重要組成部分。地質(zhì)統(tǒng)計學(xué)中擬合各向異性的基本思路是求三個相互垂直方向(主軸、次軸、短軸)的變異函數(shù)，這三個方向上的變程的比值就是各向異性中軸的比例。通常，對于大多數(shù)金屬礦床，可以根據(jù)礦體走向、傾向、厚度進(jìn)行變異函數(shù)的分析。本次研究在3DMine軟件中對Zn元素樣品進(jìn)行走向(主軸)、傾向(次軸)、厚度(短軸)3個方向的實驗變異函數(shù)計算。結(jié)合銅曼礦區(qū)勘探工程間距為80 m×80 m的實際情況，在計算實驗變異函數(shù)時的基本滯后距離取勘探工程的1/2(即40 m)，滯后距誤差限為勘探工程間距的1/4(即20 m)，變異函數(shù)計算方向如表3所示。

圖5 Zn元素原始樣品分布直方圖Fig.5 The original sample distribution histogram of Zn samples

圖6 銅曼礦區(qū)全部樣品點樣長統(tǒng)計Fig.6 The point sample length statistics total samples of Tongman mine表3 Zn元素品位變異函數(shù)計算方向Table 3 The calculation direction of variation function for Zn samples

由于Zn元素樣品點在各個方向上的影響半徑各不相同，因此需要找到每個方向上影響距離的比率，即各向異性。在3DMine軟件地質(zhì)統(tǒng)計模塊中，雙擊主軸函數(shù)圖，選擇模型中的“指數(shù)模型”進(jìn)行擬合。調(diào)整指數(shù)模型曲線上的紅點，使指數(shù)模型的曲線與變異曲線盡可能形態(tài)一致，調(diào)整完主軸后再依次調(diào)整次軸和短軸。

調(diào)整曲線的過程中，右側(cè)變異函數(shù)參數(shù)及各向異性參數(shù)都在發(fā)生變化。根據(jù)經(jīng)驗判斷一個合理的變異函數(shù)的基本原則是，隨著距離增大，伽瑪值(Gamma)不斷上升，變程也越大。在所有扇區(qū)中選擇一個最符合正態(tài)分布的方向設(shè)為“主變異函數(shù)方向”，即搜索橢球體的主軸，如圖7所示。

圖7 Zn品位主軸(走向方向)變異函數(shù)曲線Fig.7 The variation function curve of main axial for Zn grade

主軸確定后，在垂直于主軸的方向上將生成一個平面，該平面分為8個扇區(qū)，在這8個扇區(qū)中，又以相同的方式，找到一個最符合正態(tài)分布的方向并確定為次軸變異函數(shù)方向，如圖 8所示。

圖8 Zn品位次軸(傾向方向)變異函數(shù)曲線Fig.8 The variation function curve of secondary axial for Zn grade

當(dāng)確定主軸和次軸方向后，短軸方向?qū)⒆詣哟_定，如圖9所示。

圖9 Zn品位短軸(厚度方向)變異函數(shù)曲線Fig.9 The variation function curve of short axial for Zn grade

擬合完畢后，將此工程保存起來，計算獲得的參數(shù)將用于為Zn品位模型估值。擬合結(jié)果如表4所示。

表4 Zn元素理論變異函數(shù)擬合參數(shù)Table 4 The fitting parameters of variation function for Zn samples

變異函數(shù)的擬合參數(shù)主要有塊金值、基臺值和變程。其中，塊金值由于變量空間分布的不均勻性和測量誤差的存在，在最小采樣尺度下變量的變異性，可反映出區(qū)域變化量隨機(jī)性的大小?；_值表示變量在空間中的總變異性，即h大于變程時變差函數(shù)的返回值。變程是指區(qū)域變化量在空間上具有相關(guān)性的范圍。在變程范圍之內(nèi)，數(shù)據(jù)具有相關(guān)性，但在變程之外，數(shù)據(jù)之間的相關(guān)性減弱直至消失，用超出變程之外的數(shù)據(jù)對未知點進(jìn)行估值等同于數(shù)學(xué)平均。

理論變異函數(shù)參數(shù)將用于估算塊體模型中的Zn品位，這對Zn品位估值的準(zhǔn)確性有很大的影響。因此，在估值前應(yīng)當(dāng)對變異函數(shù)的參數(shù)進(jìn)行交叉驗證，即對應(yīng)用這些參數(shù)進(jìn)行品位估值時的可靠性進(jìn)行初步判斷[10]。理論變異函數(shù)參數(shù)的可靠性通?；谝韵聝蓚€方面的交叉驗證結(jié)果來判斷[11]：

1)交叉驗證的原始均值和估計均值趨于相等，交叉驗證的誤差均值應(yīng)趨近于0。

2)誤差方差和誤差平均值趨于相等且盡可能小。

表5 交叉驗證表Table 5 Cross validation Table

圖10 Zn元素樣品點估值的殘差圖Fig.10 The residual chart of Zn samples’ estimation

從表5可以看出，原始均值0.501 0與估計均值0.501 9之差趨于0，誤差平方均值0.180 7和誤差方差0.180 7相等，說明所選的Zn元素品位模型結(jié)構(gòu)模型較好，變異函數(shù)擬合科學(xué)合理。從圖 10可看出，Zn元素的大部分樣品點的散點品位1%～8%，且分布在殘差0線(虛線)附近，說明殘差值較小，與2.2節(jié)中Zn元素的樣品點分位數(shù)區(qū)間統(tǒng)計結(jié)果相吻合，故可以獲得較好的估值效果。

2.5 搜索橢球體

搜索橢球體是一個代表各向異性的體，它與待估塊體的中心重合。一般情況下，橢球體的長軸方向與礦化域走向一致，次軸和短軸則分別與礦化域的傾向方向和厚度方向一致[1]。搜索橢球體的大小由各個軸向的直徑?jīng)Q定，而各個軸向的直徑由兩個因素決定，即與待估塊相鄰的最近工程之間的距離和礦化域各個方向的變異性。

圖11 橢球體示意圖Fig.11 Schematic diagram of ellipsoid

結(jié)合銅曼礦區(qū)礦體形態(tài)特點、變異函數(shù)模擬結(jié)果，搜索橢球體參數(shù)設(shè)置如表6所示。

表6 搜索橢球體參數(shù)Table 6 Ellipsoid’s parameters

2.6 克里格法

1951年，Krige提出了一種為樣品點賦值，讓塊體品位成為樣品分析結(jié)果的線性組合，使樣品分析結(jié)果與權(quán)重漸近的方法，1963年Matheron將這種漸近的估值方法概括命名為“克里格法”[12]。一般來說，克里格法是一種尋求最優(yōu)、線性、無偏內(nèi)插值估計的方法。它是在考慮了樣品信息的大小、形狀與待估值點之間的空間分布位置等特征以及區(qū)域化變量的空間結(jié)構(gòu)信息的前提下，給每個樣品值分別賦予一定的權(quán)重系數(shù)后，用加權(quán)平均法評估待估值點的方法[13]?？死锔穹ㄊ且粋€獲得未知變量的估計方差最小化(最佳線性無偏估計)的隨機(jī)過程?？死锔穹ü乐抵校胀死锔穹ㄊ琴Y源儲量估算中較為常用的方法[1]。普通克里格法的估算公式如下：

克里格法估值參數(shù)如表7所示。

表7 克里格法估值參數(shù)Table 7 Estimation parameters of Kriging

2.7 估值結(jié)果

應(yīng)用普通克里格估值法對銅曼礦區(qū)進(jìn)行儲量估算，截至2016年末，銅曼礦區(qū)內(nèi)工業(yè)礦石共計7 979.50萬t，平均Zn品位4.19%。為驗證估值結(jié)果的可靠性，將普通克里格法的儲量估算結(jié)果與傳統(tǒng)儲量估算結(jié)果進(jìn)行對比，對比情況如下。

表8 普通克里格法與傳統(tǒng)儲量計算法估值結(jié)果對比Table 8 Estimation results comparison between ordinary Kriging and traditional method

圖12 Zn元素品位—礦石量曲線圖Fig.12 Grade-reserves curve diagram of Zn samples

由表8可看出，普通克里格法的儲量估算結(jié)果與傳統(tǒng)儲量估算結(jié)果相對誤差在合理范圍內(nèi)，無顯著性差異。從圖12可看出，普通克里格法的估值結(jié)果中Zn元素的品位—礦石量曲線圖與圖5(a)中Zn元素樣品點分布規(guī)律圖基本吻合，說明克里格法估值結(jié)果可靠。

3 結(jié)論

本文基于3DMine三維礦業(yè)軟件，對銅曼礦區(qū)建立了三維地質(zhì)模型，闡述了普通克里格法儲量估算應(yīng)用過程中樣品點統(tǒng)計、參數(shù)的選取、變異函數(shù)擬合、橢球體參數(shù)的確定過程。利用普通克里格法對銅曼礦區(qū)礦體儲量進(jìn)行了估算，估算結(jié)果表明基于3DMine軟件的普通克里格法儲量估算結(jié)果相對準(zhǔn)確，可以作為礦山資源儲量管理和開發(fā)利用的依據(jù)。此外，普通克里格法估算方法能充分利用樣品信息，相比傳統(tǒng)的手工計算方式，可提高工程技術(shù)人員工作效率，值得推廣應(yīng)用。