普通克里格法在銅曼礦區(qū)儲(chǔ)量估算中的應(yīng)用

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

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

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

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

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

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

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

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

1.2 礦體模型

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

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

1.3 塊體模型

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

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

2 儲(chǔ)量估算

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

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

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

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

2.2 基本統(tǒng)計(jì)

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

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

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

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

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

2.3 樣品組合

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

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

2.4 變異函數(shù)模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表5 交叉驗(yàn)證表Table 5 Cross validation Table

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

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

2.5 搜索橢球體

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

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

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

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

2.6 克里格法

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

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

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

2.7 估值結(jié)果

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

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

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

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

3 結(jié)論

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