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

李 江，徐江嬿

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

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

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

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

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

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

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

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

(1)

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

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

(2)

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

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

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

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

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

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

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

(3)

(4)

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

(5)

(6)

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

(7)

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

P{a≤χ2≤b}=

(8)

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

當(dāng)利用估計(jì)方差(S2)對(duì)母體誤差(σ2)進(jìn)行估計(jì)時(shí)，則有公式[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)可計(jì)算出不同樣本容量(n)下估計(jì)方差(S2)落在其母體誤差(σ2)指定區(qū)間的概率P，P值即表示估計(jì)方差的可靠性。因此，在克里格估值計(jì)算過程中，根據(jù)克里格方差和理論變差函數(shù)確定的鄰域范圍內(nèi)的參考點(diǎn)數(shù)目，結(jié)合預(yù)先劃定的無偏誤差區(qū)間，即可獲得以概率為表達(dá)的插值結(jié)果精度評(píng)定指標(biāo)，并以此為依據(jù)轉(zhuǎn)化為克里格插值算法可信度。

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

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

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

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

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

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

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

Step3：通過對(duì)D_CF和SMP_CF的疊加計(jì)算，生成結(jié)構(gòu)模型可信度(SM_CF)。

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

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

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

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

4 應(yīng)用實(shí)例與分析

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

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

圖3 城門山礦區(qū)Cu礦體結(jié)構(gòu)模型

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

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

5 結(jié)語

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

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