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

李 迪，胡乃聯(lián) ，李國清，侯 杰

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

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

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

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

儲量估算方法有很多，大體包括傳統(tǒng)方法、SD法以及地質(zhì)統(tǒng)計學方法等。近年來，隨著神經(jīng)網(wǎng)絡(luò)、遺傳算法、分形等理論迅速發(fā)展，相繼出現(xiàn)結(jié)合這些理論的新興估值方法。在眾多估值方法中最常用，并廣泛引入礦業(yè)軟件的有距離冪次反比法與克里格法[3]。

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

克里格法是以礦化空間結(jié)構(gòu)理論為基礎(chǔ)，區(qū)域化變量為核心，變異函數(shù)為工具的數(shù)學地質(zhì)方法，是地質(zhì)統(tǒng)計學中的核心方法，在國內(nèi)外各學科領(lǐng)域有著廣泛的應(yīng)用。在礦業(yè)領(lǐng)域，地質(zhì)統(tǒng)計學不僅廣泛應(yīng)用于資源儲量計算，在其他礦業(yè)活動，如異常評價[6]、找礦勘探[7]、礦體圈定[8]及采礦設(shè)計[9]等方面也具有明顯的優(yōu)越性。

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

2 基于序貫高斯模擬的動態(tài)更新技術(shù)

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

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

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

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

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

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

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

(1)

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

(2)

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

(3)

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

(4)

兩者具有相同的數(shù)學期望、協(xié)方差和變差，Zsk×(u)為Zs(u)的克立格估計值.由于Zs(u)與Z(u)同構(gòu)，Zsk×(u)與Z×(u)同構(gòu)，則其克立格方程組相同，其權(quán)系數(shù)解λi也應(yīng)相同[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 序貫高斯模擬與克里格法的關(guān)聯(lián)

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

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

3 工程應(yīng)用

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

3.1 技術(shù)路線圖

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

圖1 技術(shù)路線圖

3.2 初步儲量估算模型

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

圖2 礦體三維實體模型

圖3 估值后礦體品位分布

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

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

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

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

圖4 爆區(qū)散點分布圖

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

圖5 變異函數(shù)擬合

圖6 序貫高斯模擬結(jié)果

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

圖7 更新后的模型

3.4 結(jié)果分析

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

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

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

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

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

4 結(jié)論

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

1) 炮孔巖粉數(shù)據(jù)較地質(zhì)鉆孔數(shù)據(jù)網(wǎng)度大幅加密，將兩種數(shù)據(jù)融合，充分利用有效數(shù)據(jù)是提高儲量估算精度的基礎(chǔ)。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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