基于序貫高斯模擬的儲量估算分層更新技術

李 迪，胡乃聯 ，李國清，侯 杰

(北京科技大學金屬礦山高效開采與安全教育部重點實驗室，北京 100083)

礦產資源儲量是礦山的生命線，也是礦業項目評價、礦床開發和礦山生產中最主要的指標之一[1-2]。精確估算礦產儲量對礦山生產管理和規劃具有重要的指導意義，目前礦山通常采用具有無偏估計的克里格法進行儲量估算。

隨著礦山實際生產，將產生海量炮孔巖粉數據，這些高密度數據的加入能夠大幅提高儲量估算精度。為此，迫切需要將儲量估算模型進行更新。更新后的模型用于指導礦山實際作業，如配礦、開采設計、選礦等后續工作，因而需要模型的品位分布有高度準確性。然而，儲量估算是以地質勘探數據為基礎的靜態計算過程，其更新問題極其復雜，原因如下：其一，炮孔巖粉數據與用于初步估值的地質鉆孔數據網度、樣長不同，如何進行有效融合需要系統討論；其二，常規估值方法有諸多弊端，如距離冪次反比法缺乏考慮地質因素，克里格法帶有平滑效應，均不能較好的模擬實測數據的品位分布及其空間自相關性，需尋求新方法；其三，儲量估算流程相對固定和復雜，地質工作者主觀決策因素較多，使得模型便捷更新較為不易。為此，本課題采用常用于儲層建模的條件模擬方法，并以國內某大型露天鉬礦為工程背景嘗試其在礦床儲量模型更新中的應用。

1 常用儲量估算方法及分析

儲量估算方法有很多，大體包括傳統方法、SD法以及地質統計學方法等。近年來，隨著神經網絡、遺傳算法、分形等理論迅速發展，相繼出現結合這些理論的新興估值方法。在眾多估值方法中最常用，并廣泛引入礦業軟件的有距離冪次反比法與克里格法[3]。

距離冪次反比法以其原理簡單、插值過程快速高效的特點，廣泛應用于各領域[4]。國內外諸多軟件引入距離冪次反比法作為內部插值計算方法，如礦業軟件3Dmine、MicroMine、Surpac，遙感軟件ERmapper[5]及統計分析軟件等。距離冪次反比法是一種常用而簡便的空間插值方法，它以插值點與樣本點間的距離為權重進行加權平均。然而該方法單單只考慮了距離與待估值之間的關系，并沒有考慮到空間地質因素所造成的空間相關性，具有一定的弊端。為提高估值精度，人們廣泛應用地質統計學中經典的克里格法進行估值。

克里格法是以礦化空間結構理論為基礎，區域化變量為核心，變異函數為工具的數學地質方法，是地質統計學中的核心方法，在國內外各學科領域有著廣泛的應用。在礦業領域，地質統計學不僅廣泛應用于資源儲量計算，在其他礦業活動，如異常評價[6]、找礦勘探[7]、礦體圈定[8]及采礦設計[9]等方面也具有明顯的優越性。

克里格法遵循無偏估計和最優估計原則，在插值計算中引入區域化變量理論，因而估算精度相對較高，對于整個礦體的儲量估算工作，只要求得到礦床儲量而言，是相對合理的方法。然而，該方法的缺陷是其平滑效應[10]，這種平滑效應使得礦石品位均化，不能完好的再現具體的品位分布及波動，因而對指導礦山具體配礦、開采設計等相關具體作業的儲量更新工作而言，精度欠佳。

2 基于序貫高斯模擬的動態更新技術

條件模擬技術是地質統計學繼克里格族群之后，迅速發展的一個新分支，它的結果是每個待估點的可能取值以及其發生的條件概率，該方法并完全尊重于實測數據。隨著條件模擬各種方法的出現和不斷完善，其技術在許多方面都得到廣泛的應用。條件模擬技術最初由A.Journel教授在他的著作“礦業地質統計學”中提出并詳細地闡述其基本概念和算法，其后在此基礎上很多新算法相繼出現，如LU三角矩陣分解法、序貫高斯模擬、序貫指示模擬及分形模擬等方法。目前，條件模擬主要在石油儲層建模方面的得到廣泛應用，如使用條件模擬方法解決確定沉積相帶的空間分布[11]，以及模擬非均質石油儲層建模等等。另外，條件模擬的發展還體現在了一些非地質領域，如水土科學、森林資源評估等。

條件模擬在許多領域都取得了一些進展，但目前在礦山方面的應用較少，需要進一步研究。

2.1 序貫高斯模擬的基本原理

條件模擬的方法有很多，其中序貫高斯模擬能正確反映區域化變量的空間波動特點，能再現真實資源特性空間變異曲線的波動，是條件模擬中的首選算法。其特點是隨機抽取樣本模擬后將模擬值最為下次模擬的條件，序貫慮每個網格節點直到所有點模擬結束，下面簡述其基本原理。

正態得分變換:由于序貫高斯模擬嚴格要求數據服從高斯分布，因而需將不滿足條件的數據進行正態得分變換，模擬后再將數據逆變換。

隨機變量Z(u)，給定n個數據條件，則其條件累積分布函數(ccdf )為

Fz(z|(n))=Prob{Z

(1)

序貫高斯模擬嚴格要求數據服從高斯分布，因此FZ(z)需進行正態變換.若其累積分布函數FY(y)為標準高斯累積分布函數G(y)，即有正態變換

(2)

對n個Z(u)采樣的觀測數據以單增排列，則有相應于第k個Z值的累積頻率，Zk的正態轉換就是標準正態分布函數的分位數，即

(3)

模擬原理：設爆區品位數據Z(u)為滿足二階平穩假設的區域化變量，在空間任一點u處的真實值可表示為普通克立格估值Zk×(u)與其偏差之和，即

(4)

兩者具有相同的數學期望、協方差和變差，Zsk×(u)為Zs(u)的克立格估計值.由于Zs(u)與Z(u)同構，Zsk×(u)與Z×(u)同構，則其克立格方程組相同，其權系數解λi也應相同[12]，即有

(5)

則在實測點ui要求模擬值等于觀測值條件下，點u處Z(u)的條件模擬值Zsc(u)可表示為非條件模擬Zs(u)與其在各實測點偏差的克立格估計值之和，即

Zsc(u)=Zs(u)+∑λi[Z(ui)-Zs(ui)]

(6)

其中，Zs(u)和Zs(ui)可通過隨機抽樣方法得到。

2.2 序貫高斯模擬與克里格法的關聯

由上述原理可知，序貫高斯模擬值實際上是克里格插值結果與期望為0，方差等于克里格方差的隨機偏差之和。多次隨機模擬的平均值無限趨近克里格插值結果。然而克里格插值是使用已知點數據估計未知點，而序貫高斯模擬的每一次計算都加將之前模擬的數值作為已知條件加入計算，因而建立了數據的全局相關性，這是條件模擬能再現實測數據分布的根本原因。

克里格法的特點在于無偏性和估值誤差方差最小，但卻導致平滑效應，使得對全局的空間變異程度的估計偏低。因此，本文認為克里格法應主要用于計算平均結果，如初步求解礦體儲量；而序貫高斯模擬保持了全局空間變異性，并在每個待估點提供任意多個可能性，雖然任一次的模擬結果較克里格插值的可信度差，但由于取值的隨機性，多次模擬可提高其估值精度，因而更適合指導具體作業的儲量更新模型，體現數據的整體分布。

3 工程應用

本課題以國內某大型露天鉬礦為工程背景研究序貫高斯模擬在儲量估算模型更新中的應用。該礦區采礦權面積4.6 km2，開采深度由640m至0m標高，為目前國內最大單一鉬礦床。礦區最主要的Ⅰ號礦體工業礦鉬資源礦量79504.75萬t，金屬量717662.90t，平均品位0.090%。針對資源特點，采用露天開采方式，以15m為臺階高度進行分層剝離。

3.1 技術路線圖

工程應用的基本思路如下：首先根據礦山地質鉆孔數據建立礦床資源模型，并采用對數克里格法進行初步儲量估算，隨后加入生產中周期性產生的炮孔數據，將地質鉆孔與生產爆破數據融合，并采用序貫高斯模擬法計算儲量，最后以爆區為單位逐步計算，代替初步儲量估算模型，以達到逐漸更新模型，提高精度的目的。技術路線圖見圖1。

圖1 技術路線圖

3.2 初步儲量估算模型

初步儲量估算以地質勘探鉆孔數據為基礎借助礦業軟件SURPAC進行插值計算。本次儲量估算共收集到173個鉆孔，勘探網度為100m×100m，其中測斜數據736行，巖性數據1596行，品位化驗數據27942行，按照數據格式生成鉆孔數據庫。經統計，鉬品位化驗數據最小值為0，最大值為24790，均值為519.88，標準差為778.684，變異系數為1.498，偏度為6.60，峰度101.127，其頻率基本符合對數正態分布。隨后根據國家相關標準及該礦相關指標圈定礦體實體模型，如圖2所示。初步儲量估算采用無偏估計的克里格法進行計算，首先將實體內鉆孔數據進行組合，組合樣長2m，而后進行變異函數的擬合及對數克里格插值計算最后將對數組合樣進行轉化，估值結果如圖3所示。

圖2 礦體三維實體模型

圖3 估值后礦體品位分布

3.3 基于序貫高斯模擬的分層更新

為提高儲量估算及品位分布精度，必須充分利用礦山化驗數據，如已有地質鉆孔數據及階段性供給的加密炮孔巖粉數據對礦床模型進行更新計算。由于炮孔數據以臺階爆破為單位，樣長均為15m，因此選用15m臺階高度為單位分層對模型進行更新較為合理。這里將原樣長為2m的地質鉆孔以臺階高度標高為基準，15m為樣長進行樣品組合，而后融入炮孔數據，對模型進行更新。由于更新后的模型用于指導礦山后續配礦、境界優化、開采設計等實際工作，要求礦石品位分布的高度準確性，這里采用能夠真實還原全局空間變異性的序貫高斯模擬算法進行礦床三維模型的分層更新。

本研究使用SGeMS軟件中的序貫高斯模擬模塊對生產勘探的炮孔數據進行插值更新。SGeMS(Stanford Geostatistical Modeling Software)[13]是美國斯坦福大學開發的一套地質統計學程序庫，其包含了地質統計學中絕大多數算法，包括多種克里格估值模塊，條件模擬模塊等。

數據以該礦495-510臺階爆破數據為例進行研究，并根據6m×6m爆破網度以及15m臺階高度，建立6m×6m×15m的塊體。經統計，鉬品位化驗數據共227個，均值654.653ppm,方差189755，最大值2960最小值0，中位數580，其空間分布如圖4所示。

圖4 爆區散點分布圖

為得到模擬搜索橢球體參數，需將數據進行變異函數的擬合。搜索展開角選擇22.5°，生成變異函數曲線。由于組合數據為水平層狀數據，因而其短軸沒有數據，次主軸只在0度附近有少量數據。經變異函數擬合后得到塊金值為100958，方位角139.56，傾角0°，變程138.95的橢球體搜索參數,其主軸擬合過程見圖5。序貫高斯模擬的前提是數據服從高斯分布，因此需將數據進行正態得分變換，而后將擬合參數帶入SGeMS運行序貫高斯模擬模塊對數據進行插值計算，模擬啟動源隨機產生為14071788，最后將得到的結果進行正態變換的逆運算。序貫高斯模擬結果如圖6所示。

圖5 變異函數擬合

圖6 序貫高斯模擬結果

最后按照臺階爆破周期將炮孔數據進行序貫高斯模擬，并通過實體約束將模擬結果與初步模型導入Surpac自由變量塊體進行塊體融合與更新。這里得到的更新后模型是由兩種塊度構成的，分別是基于原地質數據的12m×12m×15m的塊體以及基于炮孔數據的6m×6m×15m的塊體，實現了不同塊度、不同數據來源、不同儲量級別的更新塊體，見圖7。隨著臺階爆破的進行，得到逐步精確的儲量值及品位分布。

圖7 更新后的模型

3.4 結果分析

為進一步分析擬合結果，這里針對同一組爆破數據進行克里格估值，并將序貫高斯模擬得到的結果與克里格插值結果進行對比，頻率直方圖如圖8所示。

由圖8可知，序貫高斯模擬與實測數據的頻率直方圖基本一致，體現了序貫高斯模擬與實測數據具有相同的空間自相關結構，更能正確反映品位數據空間波動特點。然而克里格估值使得估值結果過于集中，平滑效應明顯，不能真實再現實測數據的品位分布。從這一角度看，序貫高斯模擬的預測結果可信度更高。

從變異函數角度看(圖9)，序貫高斯模擬的變異函數曲線的變化趨勢更接近實測數據一些，但其變異函數值與實測數據相比有一定差距，即方差較大。相反，克里格估值的變異函數曲線變化趨勢與實測數據有明顯不同，然而其方差卻比較小。該組折線圖反映了序貫高斯模擬方法能更好的顯示品位空間變異性，但與克里格法相比方差較大的特點。

圖8 序貫高斯模擬與克星格法對比直方圖

圖9 序貫高斯模擬與克星格法變異函數對比

4 結論

本文結合礦山實際需求，針對常規儲量估算工作并不能充分利用有效數據的現狀，提出采用序貫高斯模擬算法加入炮孔數據進行礦床模型的分層更新，并以國內某大型露天鉬礦為例進行驗證，得出如下結論。

1) 炮孔巖粉數據較地質鉆孔數據網度大幅加密，將兩種數據融合，充分利用有效數據是提高儲量估算精度的基礎。

2) 加入巖粉數據的模型分層更新結果用于指導礦山的配礦等實際生產，因而采用更能反映數據全局變異性及波動性的序貫高斯模擬算法更為恰當。

3) 序貫高斯模擬相比于克里格算法更能體現實測數據的空間自相關性與變異性，即全局數據分布的準確性，但其局部方差略高于克里格算法，多次模擬可提高精度。

[1] 談樹成,虎雄崗,金艷珠,等.基于地統計學與GIS的固體礦產資源儲量估算方法研究[J].中國礦業，2013,3(22):101-105.

[2] 陳國旭,田宜平,劉剛等.資源儲量估算、圖表編制一體化與可視化系統研究[J].金屬礦山，2009,394(4):102-108.

[3] 張伊麗.基于CUDA的泛克里格算法的研究與設計[D].北京：中國地質大學(北京)，2012.

[4] 牛文杰.基于變異函數的距離加權反比法[J].東華大學學報：自然科學版，2011,37(3):362-367.

[5] 丁建華,肖克炎,婁德波,等.大比例尺三維礦產預測[J].地質與勘探，2009,45 (6):729-734.

[6] 韓天成,張振飛,呂新彪,等.東天山中段區域化探異常評價方法研究[J].地質與勘探，2011,47(5):885-893.

[7] 陳道貴,李建鋒,郭秀峰.基于克里格方差的勘探/取樣網度優化[J].現代礦業，2012,27(7):43-46.

[8] 尤超,張唯,劉修國.指示克里格法在礦體圈定中的應用研究[J].信陽師范學院學報：自然科學版，2012,25(1):130-132.

[9] 王斌,劉保順,王濤.基于Surpac的礦山三維可視化地質模型的研究與應用[J].中國礦業，2011,20(2):106-109.

[10] 李曉暉,袁峰,賈蔡,等.基于地統計學插值方法的局部奇異性指數計算比較研究[J].地理科學，2012,32(2):136-141.

[11] 于正軍,董冬冬,宋維琪,等.相帶控制下協克里金方法孔隙度預測[J].地球物理學進展，2012,27(4):1581-1587.

[12] 張銳,王瑞和,王建冬.應用序貫高斯模擬法預測儲集層潛在傷害[J].大慶石油學院學報，2009,33(2):58-62.

[13] 瞿明凱.幾種地統計學方法在縣域土壤空間信息處理上的應用與研究[D].武漢：華中農業大學，2012.