普通克里格法在銅曼礦區儲量估算中的應用

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

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

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

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

1 三維地質模型構建

1.1 地質數據庫

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

圖1 三維地質數據庫Fig.1 3D geologic database

1.2 礦體模型

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

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

1.3 塊體模型

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

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

2 儲量估算

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

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

2.1 區域變化量選擇

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

2.2 基本統計

在對塊體模型進行品位估值之前，要結合礦床的成因對礦區樣品的品位分布特征進行統計分析，為后續的品位估值提供數據基礎，以便根據礦床自身特點選擇適當的估值方案[5]。利用3DMine軟件中的地質統計學模塊對銅曼礦區的樣品點進行統計分析，研究樣品點數據的分布特征，如數據不符合正態分布則進行變換。本次研究主要分析銅曼礦區的Zn元素，統計結果如表1和表2所示。

表1 Zn元素樣品點分位數統計Table 1 Quantile statistics of Zn samples

表2 Zn元素樣品點基本統計Table 2 The basic statistics of Zn samples

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

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

2.3 樣品組合

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

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

2.4 變異函數模型

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

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

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

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

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

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

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

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

當確定主軸和次軸方向后，短軸方向將自動確定，如圖9所示。

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

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

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

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

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

1)交叉驗證的原始均值和估計均值趨于相等，交叉驗證的誤差均值應趨近于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元素品位模型結構模型較好，變異函數擬合科學合理。從圖 10可看出，Zn元素的大部分樣品點的散點品位1%～8%，且分布在殘差0線(虛線)附近，說明殘差值較小，與2.2節中Zn元素的樣品點分位數區間統計結果相吻合，故可以獲得較好的估值效果。

2.5 搜索橢球體

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

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

結合銅曼礦區礦體形態特點、變異函數模擬結果，搜索橢球體參數設置如表6所示。

表6 搜索橢球體參數Table 6 Ellipsoid’s parameters

2.6 克里格法

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

克里格法估值參數如表7所示。

表7 克里格法估值參數Table 7 Estimation parameters of Kriging

2.7 估值結果

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

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

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

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

3 結論

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