多粒度礦體模型不確定性表達與傳遞研究

李 江，徐江嬿

(1.湖北省自然資源廳信息中心,湖北 武漢, 430071；2. 湖北省自然資源廳地質(zhì)勘查基金管理中心,湖北 武漢, 430071)

礦體三維模型是開展地質(zhì)空間分析、礦產(chǎn)資源評價、儲量計算等方面工作的數(shù)據(jù)基礎和對建模區(qū)域地質(zhì)情況的直觀表達[1-2]。受礦體地質(zhì)結構的復雜性、人類對礦體認知的不完備性及礦體采集數(shù)據(jù)的稀疏性等多方面因素影響，不確定性成為礦體三維模型的固有屬性并貫穿于礦體建模過程的始終[3-5]。對礦體建模過程中的誤差來源、表達方式和傳遞機制進行分析，建立模型的不確定性傳遞模型和分析方法，將極大提高其實用價值和應用程度。國內(nèi)外研究人員對礦體三維模型不確定性進行了大量研究。徐華等[6]分析了礦體三維建模各環(huán)節(jié)中的不確定性問題, 建立了基于三維地質(zhì)模擬各環(huán)節(jié)的誤差檢查機制, 并針對不確定性極為突出的斷層模型, 設計了斷層模型質(zhì)量評價方法。Keefer[7]認為模型的不確定性來源主要包括數(shù)據(jù)的準確性和精度、數(shù)據(jù)的數(shù)量和空間分布、地質(zhì)復雜性、地質(zhì)解譯等四個方面。英國地質(zhì)調(diào)查局在DGSM項目中，采用因果關系圖對三維地質(zhì)模型不確定性的來源進行了分析，利用poor、ok、good三個等級相應的分析規(guī)則對模型質(zhì)量不確定性進行了半定量分析[8]。Brdossy等[9]提出了不確定性區(qū)間(uncertainty intervals)、模糊集(fuzzy sets)、概率范圍(probability range)、神經(jīng)網(wǎng)絡(neural networks)、混合算法(hybrid arithmetic)、模糊地質(zhì)統(tǒng)計學(fuzzy geostatistics)等一系列描述地質(zhì)不確定性的數(shù)學方法，以期采用數(shù)學方法解決模型不確定性的評估問題。史文中[10]對三維GIS的位置誤差進行研究，提出了三維地理信息系統(tǒng)中幾何特征(三維點、三維線段和三維線)誤差模型。Wellmann課題組[11-12]依托建模數(shù)據(jù)的概率分布及初始地質(zhì)模型，將地質(zhì)空間劃分成若干 Grid 或 Voxel 單元體，對層狀地質(zhì)體的不確定性及誤差敏感度進行熵值計算，實現(xiàn)了地質(zhì)體三維空間內(nèi)的不確定性分析。盡管概率論、空間統(tǒng)計學、模糊數(shù)學、粗糙集等多種理論和方法被用于礦體三維模型的不確定性傳遞和定量分析，然而不確定性在三維礦體建模整個過程中的傳播和定量描述仍然是尚未解決的難題，尤其是誤差計算方法、修正機制及其實現(xiàn)技術等仍處于探索階段，對礦體三維模型的不確定性表達與傳遞過程的研究仍需進一步提高。為此，本課題組依據(jù)礦山企業(yè)對三維模型在波動性、精確性等方面不同的特征表達需求，提出了以置信度為約束的基于語義的礦體多粒度模型建模方法，并根據(jù)信息熵理論提出了礦體三維模型的不確定性定量分析、模型誤差檢測與修正的技術框架和實現(xiàn)步驟[13]，本文在此研究基礎上，針對礦體多粒度模型的不確定性，應用不確定性推理理論的可信度(Certainty Factor，CF)方法和地質(zhì)統(tǒng)計學中的克里格算法，對礦體模型構建過程中的不確定性分布與傳遞過程進行分析，建立多粒度礦體模型的不確定性表達與傳遞模型框架,以期為礦山企業(yè)利用礦體三維模型實現(xiàn)勘探設計、施工計劃、目標決策等應用方面提供參考。

1 礦體三維建模不確定性來源與傳遞過程

礦體三維模型是利用地質(zhì)剖面、勘探鉆孔、斷層及褶皺構造等多源地質(zhì)數(shù)據(jù)結合插值擬合等算法進行數(shù)學模擬的結果，是對礦區(qū)地質(zhì)和礦產(chǎn)賦存等情況的描述和表示。在建模過程中,不確定性貫穿于地質(zhì)數(shù)據(jù)的獲取、整理、入庫及模型計算的始終，因此礦體三維模型的建立過程也是其不確定性的傳遞過程。根據(jù)建模步驟，礦體三維模型的不確定性來源可歸納為：樣本采集與測量誤差、數(shù)據(jù)不完整性和不一致性、礦體認知不確定性、內(nèi)插或外推隨機不確定性以及建模算法的不確定性[14]。礦體三維建模不確定性的分布如圖1所示。

圖1 礦體三維模型不確定性分布

2 礦體建模不確定性算法實現(xiàn)

2.1 基于可信度的不確定性推理

不確定性的表達和處理是人工智能研究的重要課題之一，屬于非線性復雜問題范疇，主要研究方法包括基于可信度、基于模糊和基于概率三個方面，其中，可信度(CF)方法是在確定性理論的基礎上，結合概率論、非概率和非形式化推理過程提出的一種不確定性推理模型，是在不確定性的知識或現(xiàn)象中引入了可信度概念，使原本模糊的知識或現(xiàn)象變得定量化和明確化。可信度又稱可信度因子或規(guī)則強度，包括證據(jù)可信度和規(guī)則可信度，取值在[-1,1]范圍內(nèi)，其定義一般通過概率給出，假設CF(H,E)表示當證據(jù)E出現(xiàn)時結論H成立的支持程度，CF(H,E)形式化定義為：

(1)

式中，MB(H,E)表示在證據(jù)E下對結論H的信任度的增加量；MD(H,E)表示在證據(jù)E下對結論H不信任度的增加量；P(H)表示結論H發(fā)生的先驗概率；P(H/E)表示當E為真時H發(fā)生的條件概率[15]。基于不確定性推理理論，本文將采用可信度方法中的CF模型對礦體三維模型的不確定性表達及傳遞進行描述，以可信度作為不確定性傳遞過程的統(tǒng)一度量單位。為簡化描述過程，本文將可信度值域[-1,1]歸一化為[0,1]，即可信度取值為0時，結果為“假”，可信度取值為1時，結果為“真”，因此可信度值越大，不確定性越小，結果越接近“真”，反之不確定性越大，結果越接近“假”。不確定性傳遞計算過程中的主要算法包括合成計算和疊加傳遞計算，其中，合成算法是將多個證據(jù)進行組合，證據(jù)間為合取關系或析取關系，即分別以單一證據(jù)的最大值或最小值作為合成結果的可信度；疊加傳遞計算是以證據(jù)可信度為真時推導出的結論可信度為知識可信度，通過與證據(jù)可信度相乘，獲得最終結論的可信度值，疊加傳遞算法描述為：

CF(H)=CF(E)×CF(H，E)

(2)

式中，CF(E )為證據(jù)可信度；CF(H,E)為知識可信度；CF(H)為疊加傳遞計算得到的結論可信度。

根據(jù)圖1中礦體三維模型不確定性分布情況，將建模過程中各種誤差產(chǎn)生的不確定性分為推理型不確定性和計算型不確定性兩種類型。

推理型不確定性是難以通過計算或統(tǒng)計方法確定的誤差，它可由專家根據(jù)地質(zhì)資料或礦體的復雜程度進行解釋和推斷并賦予相應的可信度，如構建礦體表面模型過程中，當?shù)V體形態(tài)簡單、內(nèi)部無夾石分支、礦體輪廓線上工程見礦點擺列規(guī)律且相鄰剖面間見礦點數(shù)目差別不大時，根據(jù)剖面輪廓線建立的礦體表面模型擬合較好，即解譯結果具有較高的可信度，因此該類型礦體圈定可信度可賦予較高值，如0.8～0.9，而對于存在斷層、夾石和分支等情況的復雜礦體，在根據(jù)稀疏分布的采樣數(shù)據(jù)進行剖面輪廓線重構時，該類型礦體表面模型可因人工解譯的經(jīng)驗性和個體性差異而導致模型擬合度差，因此該類型礦體圈定誤差可賦予較低可信度，如0.2～0.3。

計算型不確定性是指可根據(jù)解析法、經(jīng)典誤差理論等不確定性度量模型獲得，并以標準差表征轉(zhuǎn)換而來的可信度，如礦體三維建模過程中的勘探工程、巖性分層中的測量誤差及品位模型插值計算誤差等都可歸入該類型。測量誤差一般采用一階Taylor展開式對其不確定性進行近似描述，由于一階Taylor公式僅適用于線性傳播函數(shù)，對于非線性的礦體三維不確定性明顯存在估算精度不足的問題，如改用Taylor高階展開式又會增加整體計算的復雜性，因此，考慮到測量誤差不確定的計算難度，本文將該類不確定性歸入推理型不確定性范圍，其可信度賦值通過抽樣統(tǒng)計方式獲得。

礦體三維模型的不確定性表達與傳遞是一個復雜的多環(huán)節(jié)關聯(lián)與耦合系統(tǒng)，多個輸入與輸出并存，本文以誤差較為集中的礦體三維屬性模型插值計算為對象，通過對該類計算型不確定性的算法解析，結合礦體三維建模各環(huán)節(jié)不確定性分布狀況，對礦體多粒度模型的不確定性表達與傳遞機理進行探究。

2.2 礦體三維模型不確定性表達算法解析

多粒度礦體三維模型是根據(jù)礦體不同地質(zhì)階段逐步豐富的鉆孔采樣數(shù)據(jù)，以地質(zhì)統(tǒng)計學理論為數(shù)學基礎，分別采用克里格估值和序貫高斯模擬兩種不同插值算法，為滿足礦山企業(yè)不同層次需求所建立的具有多粒度特性的礦體模型，基本原理為：在礦體詳查、勘探等采樣數(shù)據(jù)稀疏且不完整階段，利用序貫高斯模擬插值的整體統(tǒng)計性和空間相關性等特點，采用多次模擬的平均值建立可反映地質(zhì)變量波動性的礦體三維模型；在礦體開拓、回采等鉆孔密集、地質(zhì)數(shù)據(jù)完備階段，根據(jù)克里格估值的無偏估計和方差最優(yōu)計算原則，建立精度較高且可精確估算儲量的礦體三維模型。同時，根據(jù)序貫高斯模擬數(shù)學期望無限趨近于克里格估值的特性，在利用模擬算法建立礦體三維模型過程中，以克里格估算和序貫高斯模擬的平均值分別作為真值和近似值，以兩者置信區(qū)間的置信度α為約束條件，根據(jù)α的不同取值定義不同的語義粒度，采用t檢測方法得到對應的模擬計算次數(shù)，進而得到對應于克里格估值和由多次序貫高斯模擬均值所構建的不同語義粒度礦體三維模型。

(3)

(4)

在依據(jù)樣本方差估計母體方差的過程中，樣本方差即估計方差S2，其定義為[17]：

(5)

(6)

由式(6)可知，利用估計方差對母體方差估值的精確度與樣本容量n存在必然聯(lián)系，設母體方差σ2的誤差區(qū)間為[(1-α)σ2,(1+α)σ2]，P為估計方差落在此誤差區(qū)間的概率，因此χ2(n)分布的概率密度函數(shù)為：

(7)

故a≤χ2≤b的概率為：

P{a≤χ2≤b}=

(8)

式中，Γ(n)為伽馬函數(shù)。

當利用估計方差(S2)對母體誤差(σ2)進行估計時，則有公式[18]：

(1-α)σ2≤S2≤(1+α)σ2

(9)

(1-α)(n-1)≤χ2≤(1+α)(n-1)

(10)

則有：

a=(1-α)(n-1)

b=(1+α)(n-1)

(11)

將式(11)帶入式(8)中，可得：

P{(1-α)(n-1)≤χ2≤(1+α)(n-1)}

(12)

根據(jù)伽馬函數(shù)定義，由式(12)可計算出不同樣本容量(n)下估計方差(S2)落在其母體誤差(σ2)指定區(qū)間的概率P，P值即表示估計方差的可靠性。因此，在克里格估值計算過程中，根據(jù)克里格方差和理論變差函數(shù)確定的鄰域范圍內(nèi)的參考點數(shù)目，結合預先劃定的無偏誤差區(qū)間，即可獲得以概率為表達的插值結果精度評定指標，并以此為依據(jù)轉(zhuǎn)化為克里格插值算法可信度。

3 多粒度礦體三維建模的不確定性傳遞過程

礦體三維建模過程中的不確定性傳遞即模型構建中各環(huán)節(jié)不確定性的合并、疊加和傳輸至下一階段的過程，為此，提出基于CF模型描述的多粒度礦體三維建模不確定性傳遞模型推理網(wǎng)絡如圖2所示，利用推理網(wǎng)絡中可信度的傳遞來實現(xiàn)多粒度礦體三維建模的不確定性傳遞的描述。 從圖2中可以看出，推理網(wǎng)絡的初始證據(jù)可信度來源于礦山多源建模數(shù)據(jù)，多源建模數(shù)據(jù)集合中的所有不確定性經(jīng)合成算法生成建模數(shù)據(jù)可信度，建模數(shù)據(jù)可信度與同樣經(jīng)合成算法生成的結構建模過程可信度進行疊加計算，形成礦體三維結構模型可信度。在利用結構模型建立多粒度屬性模型過程中，根據(jù)模型用途分別采用估值建模和模擬建模方式進行插值計算，以此形成估值算法

圖2 多粒度礦體三維模型不確定性傳遞過程

可信度和模擬算法可信度，經(jīng)與結構模型可信度進行疊加計算，最終完成三維屬性模型的可信度表達。由前文可知，模擬建模過程由克里格估值和序貫模擬計算組成，因此其可信度由克里格估值可信度和多粒度約束可信度疊加計算生成，其中，多粒度約束可信度則由預先定義的約束置信度α轉(zhuǎn)換而來。

多粒度礦體三維建模不確定性計算流程如下：

Step1：測量、解譯和整理勘探地質(zhì)剖面、鉆孔樣品等多源礦山建模數(shù)據(jù)，根據(jù)數(shù)據(jù)類型和專家經(jīng)驗，對所有數(shù)據(jù)產(chǎn)生的誤差分別進行不確定性賦值并形成礦山多源數(shù)據(jù)不確定性集合，對集合內(nèi)不確定性進行合成計算，生成建模數(shù)據(jù)可信度(D_CF)。

Step2：對礦體圈定、剖面輪廓線提取和重構、剖面多邊形對應、礦體剖分等構模過程誤差進行不確定性賦值，形成結構建模過程不確定性集合，對集合內(nèi)不確定性進行合成計算，生成建模過程數(shù)據(jù)可信度(SMP_CF)。

Step3：通過對D_CF和SMP_CF的疊加計算，生成結構模型可信度(SM_CF)。

Step4：根據(jù)礦體數(shù)據(jù)基礎和礦山企業(yè)需求的不同，對建立的具有多粒度特征的三維屬性模型按模擬建模和估值建模兩種類型分別計算其過程不確定性，即克里格估值可信度和模擬計算可信度。

(1)克里格估值可信度計算過程為：(a)經(jīng)過實驗變差函數(shù)計算和理論模型函數(shù)擬合，并根據(jù)搜索域確定插值樣本和數(shù)目；(b)對樣本進行克里格插值計算，建立礦體三維屬性模型并得到克里格估計方差；(c)根據(jù)模型精度可接受范圍，確定模型無偏誤差區(qū)間，根據(jù)公式(12)和插值樣本數(shù)計算估計方差可靠性概率P，將P值轉(zhuǎn)換為克里格估值可信度(KC_CF)。

(2)模擬計算可信度的計算過程為：(a)對樣本進行克里格估值，得到模擬計算目標值，同時依據(jù)上述克里格估值可信度計算過程，得到克里格估值可信度(KC_CF)；(b)將預先劃定的置信度α轉(zhuǎn)換為多粒度約束可信度(MGC_CF),同時以模擬計算目標值為真值，以置信度α為約束，利用進行多次序貫高斯模擬計算的平均值構建不同粒度的礦體三維屬性模型；(c) 對KC_CF和MGC_CF進行疊加計算，生成模擬計算可信度(SC_CF)。

Step5：將Step3得到的結構模型可信度(SM_CF)分別與模擬建模和估值建模過程產(chǎn)生的克里格估值可信度(KC_CF)及多粒度約束可信度(MGC_CF)進行疊加計算，得到多粒度礦體三維屬性模型可信度(AM_CF)。

4 應用實例與分析

在江西省城門山銅礦地質(zhì)勘查數(shù)據(jù)的基礎上，以Cu品位模型建立過程中的克里格插值計算為例，對誤差較為集中的屬性建模中克里格估值不確定性計算的有效性進行驗證。

城門山銅礦為Cu、S、Mo共生的綜合性礦床，伴生Au、Ag等稀有元素，主要礦體產(chǎn)狀包括似層狀、豆莢狀、透鏡狀、帶狀及席狀等類型。礦區(qū)勘查數(shù)據(jù)來源于礦山地質(zhì)過程中鉆孔資料，鉆孔樣品數(shù)共18 421件，結合樣品中Cu、S、Mo等元素化學成分分析結果，完成20條勘探線和159個鉆探工程數(shù)據(jù)的編錄工作。以iExploration-EM(固體礦產(chǎn)資源量估算與礦體三維建模系統(tǒng))為數(shù)據(jù)處理和建模平臺，經(jīng)單工程礦體圈定、剖面礦體圈定、輪廓線拼接等人工解譯過程，形成實驗礦區(qū)Cu礦體結構模型如圖3所示。從圖3中可看出，礦區(qū)共包含三個礦體，分別為礦體Ⅰ、礦體Ⅱ、礦體Ⅲ。經(jīng)鉆孔數(shù)據(jù)的特高品位值處理、鉆孔樣品組合、實驗變差函數(shù)計算、理論變差函數(shù)結構分析與套合等步驟，建立Cu屬性模型,在建模過程中，選取礦化較集中的12至14號勘探線區(qū)間大小為80 m×120 m×290 m的區(qū)域作為驗證區(qū)塊，該區(qū)塊包含于礦體Ⅰ內(nèi)，在驗證區(qū)塊內(nèi)隨機選取20個采樣點作為驗證塊體，采樣點Cu品位平均值為0.853，方差σ2為0.962。 將全樣本的誤差區(qū)間分別設定為[0.9σ2,1.1σ2]和[0.8σ2,1.2σ2]兩個區(qū)域， 在驗證區(qū)塊內(nèi)，首先利用不同數(shù)目的樣本對選定的20個采樣點經(jīng)30次克里格法插值

圖3 城門山礦區(qū)Cu礦體結構模型

計算分別得到其估值方差及其落在預定誤差區(qū)間的概率，再利用公式(12)計算出相對應樣本數(shù)下的計算方差落在預定誤差區(qū)間的概率P，比較兩種方法的計算結果如表1所示。從表1中可以看出，在相同的樣本數(shù)下，克里格估值方差在預定區(qū)間的概率與利用公式(12)計算得到的計算方差落在指定區(qū)間的概率非常接近，表明兩種方法得到的結果偏差較小，因此，在對克里格估值可信度進行確定的過程中，可以通過其插值樣本的數(shù)目，由式(12)計算獲得在指定誤差區(qū)間內(nèi)的以概率為表達的插值結果精度評定指標。

表1 克里格估值方差與數(shù)值計算方差落在指定誤差區(qū)間概率比較(單位：%)Table 1 Comparison of the falling probability between Kriging interpolation variance and numerical calculated variance within the specified error interval

5 結語

(1) 對多粒度礦體三維模型的不確定性采用克里格插值計算方法，并對克里格方差進行可信度轉(zhuǎn)換,提出了基于樣本數(shù)的樣本方差對母體方差的精度評估方法，并以此得到克里格估值計算過程中基于概率表達的精度評估指標，實現(xiàn)多粒度礦體三維模型的不確定性的定量評價。

(2)通過對試驗區(qū)礦體樣本方差與母體方差精確度的比較以及對樣本克里格估值方差在指定誤差區(qū)間內(nèi)的概率進行統(tǒng)計分析，驗證了多粒度礦體三維模型的不確定性的定量評價方法是有效的。