基于克里金算法的井田煤層三維建模方法研究

張小艷，王萌娟

(西安科技大學 計算機科學與技術(shù)學院，陜西 西安 710600)

0 引 言

近年來隨著計算機科學技術(shù)與工業(yè)生產(chǎn)傳統(tǒng)技術(shù)的快速融合與發(fā)展，數(shù)字化礦山已經(jīng)成為了礦產(chǎn)行業(yè)的熱點問題[1]。在建立煤層三維模型時，常用到的建模軟件有Surpac Vision[2]、Google Ocean[3]、GOCAD[4]、3DMine礦業(yè)工程軟件[5]、ArcGIS軟件[6]、MICROMINE[7]、CATIA[8]、龍軟[9]等，但采用這些軟件所構(gòu)建的三維模型很難進行跨平臺使用， 同時也無法與煤礦企業(yè)現(xiàn)有的一些Web管理系統(tǒng)進行較好的兼并，不利于企業(yè)對信息的統(tǒng)一管理；同時，建模中常用的方法是單一空間插值算法或曲面擬合算法，不能很好地解決復雜的地質(zhì)構(gòu)造，構(gòu)造的三維模型不太符合實際的煤層分布情況及空間關(guān)系[10]。因此，需要建立一個實現(xiàn)無插件、可移植、跨平臺、準確度高、支持多瀏覽器運行的三維可視化模型。

在實際生產(chǎn)過程中，鉆孔點不是按照規(guī)則化的網(wǎng)格進行獲取的，而是空間中的隨機離散點，鉆孔數(shù)據(jù)沒法做到規(guī)則化網(wǎng)格布點，也無需做如此大規(guī)模的浪費[11]，同時井田范圍大、鉆孔數(shù)據(jù)較少。為了構(gòu)建一個能夠精確地表征煤層分布情況的三維模型，該文首先采用改進的支持向量回歸算法(SVR)優(yōu)化克里金插值算法(Kriging)；然后，建立基于改進SVR的克里金高程估算模型，根據(jù)鉆孔數(shù)據(jù)和煤層高程信息構(gòu)建估算模型，對未知點的高程進行了估算，增加實現(xiàn)三維建模的數(shù)據(jù)點；最后，采用three.js技術(shù)和三角剖分算法[12]實現(xiàn)煤層三維模型，其能夠較為精確地展示出煤層分布范圍以及井田煤層空間關(guān)系。

1 理論基礎

1.1 克里金插值算法

克里金算法是地質(zhì)統(tǒng)計學中的一種空間插值技術(shù)，廣泛用于各類觀測的空間插值，地面風場、降雨、土壤、環(huán)境污染等領(lǐng)域。克里金插值法是根據(jù)采樣點與估算點之間的位置關(guān)系與關(guān)聯(lián)程度，由已知點的數(shù)據(jù)加權(quán)求和來估計未知點的值。它是一種具有無偏、線性特征的最優(yōu)統(tǒng)計方法[13]。設x1,x2,…,xn為區(qū)域內(nèi)的n個采樣點，z(x1),z(x2),…,z(xn)為對應采樣點的觀測值。區(qū)域化變量在x處的估算值z*(x)表示為：

(1)

其中，λi(i=1,2,…,n)為權(quán)重系數(shù)，并且各權(quán)重系數(shù)的取值需滿足兩個條件：無偏性和最優(yōu)性。權(quán)重系數(shù)λi的求解有兩個步驟：首先，求得空間區(qū)域化變量對應的變差函數(shù)；其次，再進行插值[14]。

變差函數(shù)是用來表示區(qū)域化變量的空間結(jié)構(gòu)特征，是克里金插值的基礎。在求解變差函數(shù)模型中要求得三個參數(shù)：變程、塊金值和偏基臺值，這些參數(shù)反映了空間變化特征。在得到實驗變差函數(shù)曲線之后，根據(jù)曲線的趨勢去選擇合適的理論變差函數(shù)模型，常見的三種理論變差函數(shù)模型有：球狀模型、指數(shù)模型、高斯模型。

1.2 支持向量回歸模型

在回歸模型中，給定一組數(shù)據(jù){(x1,y1),(x2,y2),…,(xm,ym)}，其中，xi∈Rd，yi∈R，希望學習到一個模型f(x)=ωTx+b，使得f(x)與y盡可能接近。支持向量機假設能容忍f(x)與y之間最多有ε的偏差，即僅當|f(x)-y|>ε時，才會計算損失[15]。SVR可以將模型轉(zhuǎn)化為：

(2)

SVM算法在處理小樣本、非線性和高維數(shù)等問題時有著許多獨特優(yōu)勢，當模型為非線性回歸時，SVR引入核函數(shù)，構(gòu)造映射函數(shù)φ(x)將樣本從原始空間映射到合適的可線性回歸的特征空間中，從而使其線性回歸。

1.3 可修改變異方向的差分進化算法

差分進化算法首先對種群進行初始化，然后通過對種群中所有個體進行變異、交叉操作，產(chǎn)生中間個體，將中間個體與當前種群中對應個體進行比較，選擇出適應值較好的個體作為下一代的新個體。隨著種群的進化，保留“基因”優(yōu)秀的個體，確保最優(yōu)個體逐步收斂至全局最優(yōu)解[16]。為了提高差分進化算法的搜索效率，選用一種可以在變異過程中修正變異方向的自適應差分進化算法(UMDE)[17]，其中自適應的縮放因子定義為：

(3)

其中，fi(pre)是指當前個體的上一代個體的適應值，fmin(pre)是與上一代最優(yōu)個體的適應值，fi(cur)為當代個體的適應值。如果在第一次變異之后，發(fā)現(xiàn)有個體適應值不如上一代個體適應值，則說明此次變異偏移了最優(yōu)個體適應值的逼近方向，此時就需要適當調(diào)整縮放因子F，使其在下一次進化過程中向全局最優(yōu)個體適應值的方向發(fā)展[17]。

2 UMDE-SVR-Kriging高程估算模型

地質(zhì)統(tǒng)計學中的克里金插值算法已被應用在高程估算中，多位學者也證明了該算法有著較好的估算精度。因此，選擇克里金插值算法對煤層高程進行估算。在克里金插值算法中，變差函數(shù)的求解是克里金插值算法的核心[18]。為了提高克里金算法的估算精度，該文采用UMDE優(yōu)化后的SVR擬合克里金的實驗變差函數(shù)。首先根據(jù)已有的鉆孔坐標信息進行計算，得到實驗變差函數(shù)值，并求得步長進行分組，然后求解SVR最優(yōu)參數(shù)組，進而對變差函數(shù)進行擬合。采用改進的支持向量回歸算法擬合變差函數(shù)能夠克服常規(guī)克里金對變差函數(shù)選擇的依賴性、局限性，以及變差函數(shù)選擇的主觀性和人為性，進而提高克里金插值算法的估算精度。

2.1 基于差分進化算法優(yōu)化的支持向量回歸

在SVR的參數(shù)求解中，運用UMDE算法優(yōu)化支持向量回歸，隨著進化代數(shù)不斷的增加，通過變異、交叉、選擇操作，取得SVR最優(yōu)的參數(shù)對(ε,C,γ)，求解步驟如下：

Step1：確定所求參數(shù)(ε,C,γ)；

Step2：初始化種群，其中差分進化算法的種群數(shù)量為100，縮放因子為0.6，交叉概率為0.65，遺傳代數(shù)為500；

Step3：計算種群中各個個體的適應度值；

Step4：條件判斷，若達到最大進化代數(shù)，則停止，輸出最優(yōu)結(jié)果，否則繼續(xù)執(zhí)行下一步；

Step5：執(zhí)行差分進化操作，進行交叉、變異以及選擇操作，生成新一代；

Step6：進化代數(shù)加一，然后跳轉(zhuǎn)到第4步。

通過以上步驟的運算，便可求解到SVR最優(yōu)參數(shù)對(ε,C,γ)。采用UMDE優(yōu)化SVR，通過種群的不斷迭代，UMDE在交叉、變異、選擇操作后，求解出SVR的最優(yōu)參數(shù)組，提高了SVR的擬合精度，為構(gòu)建基于SVR的克里金高程估算模型奠定了基礎。

2.2 UMDE-SVR-Kriging高程估算模型建立流程

基于UMDE-SVR的克里金算法高程估算模型建立流程如下：

Step1：準備構(gòu)建模型所需的鉆孔數(shù)據(jù)、煤層高程數(shù)據(jù)、煤層厚度數(shù)據(jù)，并對數(shù)據(jù)進行正態(tài)分布檢驗，若不滿足，對數(shù)據(jù)進行處理，使其正態(tài)分布；

Step2：將觀測點進行兩兩配對，計算每對數(shù)據(jù)點之間的距離以及半方差值。例如有采樣點(x1,x2,x3)，將數(shù)據(jù)點進行兩兩配對之后有數(shù)據(jù)點對(x1,x2)，(x1,x3)，(x2,x3)，計算每對數(shù)據(jù)點之間的距離和半方差值，點對之間距離時選擇歐氏距離，如公式(4)所示，半方差值的公式為公式(5)；

(4)

(5)

Step3：利用采樣點數(shù)據(jù)計算變差函數(shù)值可以得到實驗變差函數(shù)云圖，它以離散且不規(guī)則的形式存在，需要對其進行處理才能使用。因此需要對步驟2中求得的距離進行分組，并計算組內(nèi)變差函數(shù)值的平均值作為該組的值，在分組時選擇平均最近距離作為步長對其進行分組；

Step4：擬合實驗變差函數(shù)，采用2.1小節(jié)中的UMDE-SVR算法擬合變差函數(shù)，求得最優(yōu)的SVR參數(shù)對，最終確定最優(yōu)理論變差函數(shù)；

Step5：根據(jù)求得的變差函數(shù)模型構(gòu)建的高程估算模型，對未知點的高程值進行插值估算。

變差函數(shù)反映的是區(qū)域化變量的空間變化結(jié)構(gòu)，是克里金插值算法的基礎。大多數(shù)變差函數(shù)模型的求解，往往是依據(jù)計算得到的變差函數(shù)值所繪制的實驗變差函數(shù)云圖的圖形曲線，選擇合適的理論變差函數(shù)模型，因此，變差函數(shù)模型的選擇較為局限，影響了模型的精度。采用UMDE-SVR來對克里金插值算法的變差函數(shù)進行擬合，構(gòu)建符合變差函數(shù)的最優(yōu)理論變差函數(shù)模型，避免了選擇常規(guī)的幾種模型的局限性，以及模型選擇的主觀性，從而使得高程估算模型有著較高的估算精度。

3 實驗與結(jié)果分析

3.1 實驗驗證與分析

本次實驗數(shù)據(jù)來源為某煤礦的真實數(shù)據(jù)，包括鉆孔數(shù)據(jù)、煤層數(shù)據(jù)以及各煤層的煤質(zhì)信息，有705個鉆孔數(shù)據(jù)，9層煤層數(shù)據(jù)，其中有8層煤為可采煤層。根據(jù)實驗要求對數(shù)據(jù)進行處理，從鉆孔數(shù)據(jù)表中提取鉆孔編號、鉆孔坐標X和Y；在煤層信息表中，將各個煤層的煤層信息進行整理，其中包括煤層編號、鉆孔編號、煤層底板高程、煤層厚度；將各個煤層的鉆孔編號與鉆孔信息表中的鉆孔編號對應，將煤層編號、鉆孔編號、鉆孔坐標X和Y、煤層底板標高、煤層厚度進行整合，構(gòu)建了各個煤層信息表。本次實驗采用的數(shù)據(jù)主要為鉆孔坐標、煤層底板標高以及煤層厚度。

(1)實驗環(huán)境。

本次實驗的實驗環(huán)境如表1所示。

表1 實驗環(huán)境

(2)檢驗數(shù)據(jù)分布。

根據(jù)正態(tài)QQ圖可以判斷插值對象是否服從正態(tài)分布。在正態(tài)QQ圖中，如果數(shù)據(jù)點圍繞著一條直線上下波動，則服從正態(tài)分布。對于不服從正態(tài)分布的數(shù)據(jù)，可以對其采用對數(shù)變換、冪變換等方法，使其滿足正態(tài)分布。

在ArcGIS軟件中對該5-2煤層的高程數(shù)據(jù)進行空間數(shù)據(jù)探索性分析，得出其分布的正態(tài)QQ圖(見圖1)。

在圖1中，可以看到數(shù)據(jù)點近似分布為一條直線，該煤層的高程值十分接近標準正態(tài)分布，滿足插值對象為正態(tài)分布的前提條件。

圖1 5-2煤層高程值正態(tài)QQ圖

(3)實驗過程。

第一步：計算實驗變差函數(shù)值，繪制變差函數(shù)云圖。

當對5-2煤的鉆孔數(shù)據(jù)(172組)進行兩兩配對計算其距離以及半方差值后，會得到如圖2所示的實驗變差函數(shù)云圖，x軸表示點對之間的距離，也稱滯后距，y軸表示其相應的實驗變差函數(shù)值，即計算得到的半方差值。

圖2 5-2煤實驗變差函數(shù)云圖

第二步：計算步長并進行分組。

在圖2中，可以看到圖中的數(shù)據(jù)點分布散亂，很難對其進行求解，因此在實際的應用中，需要對其按照步長進行分組，求得組內(nèi)平均值并繪制變差函數(shù)云圖，然后再進行求解。本次實驗中采用平均最鄰近距離作為分組步長，根據(jù)上一步計算得到的點與點之間的距離，計算得到平均最近鄰距離值為640.720 37，在具體實驗中對其進行取整，即取分組步長為640 m，對圖2中的點對進行分組，數(shù)據(jù)點可分為27組，得到分組后的實驗變差函數(shù)云圖，如圖3所示。

圖3 5-2煤分組后實驗變差函數(shù)云圖

第三步：利用UMDE優(yōu)化支持向量回歸，得到支持向量參數(shù)。

根據(jù)分組后得到的實驗變差函數(shù)值、距離值，采用自適應差分進化算法優(yōu)化SVR，得到SVR最優(yōu)參數(shù)組，其中差分進化算法的種群數(shù)量為100，縮放因子為0.6，交叉概率為0.65，遺傳代數(shù)為500。

第四步：采用UMDE-SVR擬合變差函數(shù)并估算高程值。

根據(jù)第三步得到SVR最優(yōu)參數(shù)組(ε,C,γ)為(1.81e-12，100.0，0.62)，擬合變差函數(shù)，構(gòu)建高程估算模型。設搜索點數(shù)為30，對樣本點高程進行插值估算。

為了合理評價UMDE-SVR-Kriging模型的估算效果，選用5-2、1-2煤層數(shù)據(jù)，與UMDE-Kriging的指數(shù)、高斯、球狀模型，基于改進粒子群的SVR(Ipso-SVR)克里金法，以及基于留一法的SVR(Loo-SVR)的克里金法，在同一數(shù)據(jù)集上進行實驗，通過均方誤差(MAE)與均方根誤差(RMSE)兩個指標來評價模型的估算效果，證明UMDE-SVR-Kriging模型的高估算精度，結(jié)果如表2和表3所示。

表2 5-2煤不同高程估算模型實驗結(jié)果比較

表3 1-2煤不同高程估算模型實驗結(jié)果比較

根據(jù)表2和表3中的實驗結(jié)果，可以明顯看到，相較于UMDE-Kriging算法的指數(shù)模型、高斯模型、球狀模型，Ipso-SVR-Kriging以及Loo-SVR-Kriging，UMDE-SVR-Kriging的MAE和RMSE值較小，說明該模型估算的高程值相較于其他模型更為精準，證明了采用基于差分進化算法優(yōu)化的SVR的克里金高程估算模型在煤層高程估算方面有著較高的估算精度，以及模型在井田煤層高程估算的可行性。

3.2 三維可視化

由于現(xiàn)有井田鉆孔數(shù)據(jù)較少，畫圖無法真實反映煤層的分布情況，因此該文采用改進的SVR的克里金高程估算模型對未知點的高程值進行估算，增加三維建模數(shù)據(jù)，從而使得繪制出的煤層三維圖更符合現(xiàn)實情況，并在3.1小節(jié)中通過實驗證明了UMDE-SVR-Kriging高程估算模型的較高估算精度以及可行性。

因為煤層有著一定的厚度，煤層上下有著不同的變化規(guī)律，現(xiàn)有數(shù)據(jù)為鉆孔坐標、煤層底板高度及煤層厚度，根據(jù)底板高度和煤層厚度，計算得到某個鉆孔坐標煤層上的高度，之后分別對煤層上下的高程分別建立估算模型，實現(xiàn)對未知點的煤層上下高度的估算。

采用2.2小節(jié)中介紹的UMDE-SVR-Kriging模型分別建立煤層上下的高程值的估算模型，實現(xiàn)對未知點的高程值的估算。在數(shù)據(jù)點選擇中，以鉆孔坐標X、Y的最大值和最小值為最大范圍，隨機選取一定個數(shù)的未知點坐標，利用構(gòu)建的煤層上下高程值的估算模型對同一坐標點煤層上下高程值進行插值估算，求得未知點煤層上下的高程值。根據(jù)估算得到某點煤層上下的高程值可以得到該點煤層厚度，已知每層煤層厚度范圍，通過對估算得到的煤層厚度進行篩選，得到有效的估算點進行繪圖。

為了建立無插件、可移植、跨平臺、準確度高、支持多瀏覽器運行的煤層三維可視化系統(tǒng)，以Visual Studio為開發(fā)工具，采用vue框架等構(gòu)建了井田煤層三維可視化系統(tǒng)，three.js三維引擎和三角剖分算法，以及原始數(shù)據(jù)如表4所示，UMDE-SVR-Kriging高程估算模型估算得到的5-2煤的高程數(shù)據(jù)如表5所示，實現(xiàn)web端的5-2煤的三維模型，如圖4所示。

表4 5-2煤原始部分數(shù)據(jù)

表5 5-2煤估算得到的部分數(shù)據(jù)

圖4 5-2煤層圖

根據(jù)圖4，可以看出5-2煤層的大致范圍和分布情況，同時可對圖形進行旋轉(zhuǎn)，查看不同方向的煤層情況。

已知某井田內(nèi)可采煤層共有8層，為1-2上、1-2、2-2、3-1、4-2、4-3、4-4、5-2煤層。分別根據(jù)各個煤層的鉆孔坐標及高程數(shù)據(jù)，基于改進的SVR的克里金高程估算模型建立不同的估算模型并進行估算。根據(jù)各煤層估算得到的有效數(shù)據(jù)以及原始數(shù)據(jù)，采用Three.js三維引擎和三角剖分算法實現(xiàn)該井田8層煤的三維可視化模型，如圖5所示，從上往下分別是1-2上、1-2、2-2、3-1、4-2、4-3、4-4、5-2煤層。

圖5 某井田煤層圖

由圖5可以看出該井田的8層煤的分布和范圍情況以及空間特性，并且在該模型中，可以隨機選擇煤層號，展示該煤層的具體情況。

4 結(jié)束語

由于在實際生產(chǎn)過程中，井田原始鉆孔數(shù)據(jù)稀疏，由此數(shù)據(jù)繪制的煤層三維圖無法真實地展示出煤層分布情況。該文構(gòu)建了基于克里金插值算法的高程估算模型，對未知點高程值進行估算，增加繪圖數(shù)據(jù)點。為了提高克里金的估算精度，采用改進的支持向量回歸對克里金進行了優(yōu)化，構(gòu)建了基于UMDE-SVR優(yōu)化的克里金高程估算模型，采用某煤礦實際數(shù)據(jù)，將該模型與差分進化算法優(yōu)化的克里金的指數(shù)、球狀、高斯模型，基于改進粒子群的SVR克里金法，以及基于留一法的SVR的克里金法的估算精度進行了對比，分別對某礦5-2、1-2煤層的高程值進行估算。實驗結(jié)果表明，UMDE-SVR-Kriging模型的MAE與RMSE值都較優(yōu)，證明了該模型有著較高的估算精度，可為井田煤層三維模型的實現(xiàn)提供技術(shù)支持。通過UMDE-SVR-Kriging模型構(gòu)建的某井田煤層的三維模型，能夠更為真實地展現(xiàn)出該井田煤層的分布情況，并且可以與煤礦企業(yè)現(xiàn)有的Web管理系統(tǒng)進行良好的契合，對煤礦企業(yè)的生產(chǎn)起到了很好的指導作用。