張小艷，許慧，姜水軍

(1.西安科技大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，陜西 西安 710600; 2.中國(guó)神華神東煤炭集團(tuán)，陜西 榆林 719315)

0 引言

近年來(lái)，隨著數(shù)字礦山概念的提出，三維地質(zhì)建模在行業(yè)領(lǐng)域中備受關(guān)注。因此，對(duì)綜采工作面煤炭賦存形態(tài)空間三維可視化成為煤炭生產(chǎn)管理過(guò)程中的一個(gè)重要環(huán)節(jié)。

目前三維地質(zhì)建模在行業(yè)領(lǐng)域中備受關(guān)注。范文瑤等[1]基于GOCAD軟件對(duì)研究區(qū)域進(jìn)行三維建模與可視化操作；尚福華等[2]為了解決傳統(tǒng)三維地質(zhì)模型中不能有效地分析地層和斷層之間的形態(tài)和結(jié)構(gòu)，構(gòu)建出了一種TIN-Octree混合空間數(shù)據(jù)模型的精確三維地質(zhì)結(jié)構(gòu)方法；宋越等[3]將研究的重點(diǎn)放在了解決邊界不連續(xù)的問(wèn)題層面上，提出了一種基于機(jī)器學(xué)習(xí)的擴(kuò)展紋理樣本庫(kù)的方法。

空間數(shù)據(jù)插值是三維建模過(guò)程中的一個(gè)至關(guān)重要的環(huán)節(jié)。自20世紀(jì)初以來(lái)，就有大量的學(xué)者對(duì)此展開研究并取得了卓越的成效。Kriging算法被廣泛用于各類觀測(cè)的空間插值。本文提出了一種用自適應(yīng)PSO算法優(yōu)化Kriging插值算法(UPSO-Kriging)，將其用于工作面煤層高程值預(yù)測(cè)，為繪制工作面煤層賦存形態(tài)三維圖提供基礎(chǔ)數(shù)據(jù)，并以神東集團(tuán)大柳塔煤礦604工作面的實(shí)際數(shù)據(jù)驗(yàn)證了算法的可行性,進(jìn)而基于規(guī)則網(wǎng)格法構(gòu)建DEM數(shù)字高程模型，以Three.js為平臺(tái)實(shí)現(xiàn)了瀏覽器端綜采工作面三維可視化，為煤炭企業(yè)開采提供科學(xué)依據(jù)。

1 Kriging插值法

Kriging插值法是一種基于協(xié)方差函數(shù)對(duì)空間隨機(jī)過(guò)程建模和預(yù)測(cè)的回歸算法[4]。一般公式可表示為：

(1)

式中:待測(cè)樣本點(diǎn)的值為Z*(x)，已知樣本點(diǎn)的值以及個(gè)數(shù)分別為Z(xi)、n，第i個(gè)已知樣本點(diǎn)相對(duì)于未知樣本點(diǎn)的權(quán)重為λi[5]。

當(dāng)Lagrange乘數(shù)法求取的最小估計(jì)方差滿足本征假設(shè)時(shí)，用變差函數(shù)表示Kriging方程組如式(2)所示：

(2)

式中：μ為L(zhǎng)agrange系數(shù)，γ(xi,xj)為兩個(gè)已知采樣點(diǎn)xi,xj的變差函數(shù)值。

2 UPSO-Kriging插值法

變差函數(shù)模型擬合中的參數(shù)確定是Kriging插值過(guò)程中的關(guān)鍵所在[6-7]。本文針對(duì)傳統(tǒng)Kriging插值法在變差函數(shù)模型擬合中存在的過(guò)擬合或欠擬合而導(dǎo)致算法整體精度下降的問(wèn)題，提出一種自適應(yīng)PSO算法來(lái)優(yōu)化模型中的3個(gè)參數(shù)，提高算法整體預(yù)測(cè)精度。

2.1 自適應(yīng)PSO算法

2.1.1 PSO算法

粒子群優(yōu)化算法(PSO)是目前廣泛運(yùn)用于各個(gè)領(lǐng)域的進(jìn)化算法[8]，具有簡(jiǎn)單、準(zhǔn)確、收斂速度快等特點(diǎn)。標(biāo)準(zhǔn)PSO算法更新公式如下[9]：

c2r2[pg(t)-xi(t)],

(3)

xi(t+1)=xi(t)+vi(t+1),

(4)

2.1.2 UPSO算法

針對(duì)PSO算法存在的搜索時(shí)間過(guò)長(zhǎng)、易陷入局部極值的最優(yōu)化等問(wèn)題[13]，本文采用張小艷等[14]提出的優(yōu)化方法。分別對(duì)算法尋優(yōu)過(guò)程中的粒子位置更新方式、迭代預(yù)選值處理以及基于正態(tài)分布的慣性調(diào)整策略進(jìn)行優(yōu)化。

為驗(yàn)證UPSO算法的有效性，本文分別選取4個(gè)函數(shù)作為標(biāo)準(zhǔn)測(cè)試函數(shù)，與PSO算法在python2.7環(huán)境下進(jìn)行仿真測(cè)試。其中包括兩個(gè)單峰函數(shù)：Sphere函數(shù)與Schwefel函數(shù)；兩個(gè)多峰函數(shù)：Ackley函數(shù)與Rastrigrin函數(shù)。單峰函數(shù)用來(lái)測(cè)試算法的收斂速度以檢驗(yàn)其局部搜索能力，多峰函數(shù)用來(lái)檢驗(yàn)算法的全局搜索能力從而衡量算法綜合性能。

圖1 Sphere函數(shù)收斂測(cè)試Fig.1 Sphere function convergence test diagram

圖2 Schwefel函數(shù)收斂測(cè)試Fig.2 Schwefel function convergence test diagram

圖3 Ackley函數(shù)收斂測(cè)試Fig.3 Ackley function convergence test diagram

圖4 Rastrigrin函數(shù)收斂測(cè)試Fig.4 Rastrigrin function convergence test diagram

由上圖可知，不管是在單峰測(cè)試函數(shù)下的局部搜索能力還是多峰函數(shù)下的全局搜索能力，UPSO算法相比PSO算法來(lái)說(shuō),得到最優(yōu)解的收斂速度更快，陷入局部極值的次數(shù)明顯少于原始算法。這均表明UPSO算法整體性能得到了提升。

2.2 UPSO-Kriging插值法

UPSO算法優(yōu)化Kriging插值法，主要通過(guò)UPSO算法迭代尋優(yōu)得到最優(yōu)的模型參數(shù)進(jìn)行變異函數(shù)模型的擬合，進(jìn)而得到預(yù)測(cè)結(jié)果。算法步驟如下:

1)確定待求解問(wèn)題的參數(shù)與目標(biāo)函數(shù)。以變差函數(shù)模型中的3個(gè)模型參數(shù)作為求解參數(shù)，即塊金值c0，偏基臺(tái)值c，變程a，目標(biāo)函數(shù)為不同距離下的實(shí)驗(yàn)與理論變差函數(shù)差值平方和[15],即

(5)

2)初始化粒子群。設(shè)群體規(guī)模為M，將3個(gè)模型參數(shù)作為一個(gè)粒子個(gè)體，即xi=(c0,c,a),i=1,2,…,M。種群中每個(gè)個(gè)體的參數(shù)初始化取值如下：

(6)

其中，Rmax是實(shí)驗(yàn)變差函數(shù)的最大值，Amax為最大滯后距[16]。

3)確定UPSO算法的有關(guān)運(yùn)行參數(shù)。其中包括慣性權(quán)重最大最小取值vmax,vmin;學(xué)習(xí)因子最大最小取值c1s,c1f,c2s,c2f；迭代次數(shù)Iter；速度vmax,vmin；位置xmax,xmin。

4)UPSO優(yōu)化運(yùn)算。優(yōu)化過(guò)程中，通過(guò)目標(biāo)函數(shù)來(lái)評(píng)價(jià)群體中每個(gè)個(gè)體的適合度值是否達(dá)到標(biāo)準(zhǔn)，通過(guò)更新迭代過(guò)程中的速度和位置進(jìn)行迭代尋優(yōu)運(yùn)算。若當(dāng)前迭代中的最優(yōu)解不滿足終止條件則繼續(xù);若滿足，則停止。最終結(jié)果即最優(yōu)變差函數(shù)參數(shù)，由此可擬合出相應(yīng)的變差函數(shù)表達(dá)式。

5)Kriging插值預(yù)測(cè)。將最終得到的變差函數(shù)代入式(2)即可得到權(quán)重系數(shù)，再通過(guò)式(1)計(jì)算出待預(yù)測(cè)點(diǎn)的屬性值。

3 實(shí)驗(yàn)驗(yàn)證

3.1 數(shù)據(jù)準(zhǔn)備

本文運(yùn)用神東集團(tuán)大柳塔煤礦604綜采工作面，寬度300 m，深度3 400 m，煤層厚度在6.24～40.17 m的69組煤層煤樣數(shù)據(jù)作為樣本數(shù)據(jù)集，驗(yàn)證UPSO-Kriging算法應(yīng)用于綜采工作面煤層高程屬性預(yù)測(cè)的可行性。部分樣本數(shù)據(jù)如表1所示。

表1 采樣點(diǎn)數(shù)據(jù)

3.2 基于UPSO-Kriging的高程值預(yù)測(cè)

基于樣本數(shù)據(jù)集合中同一煤層下的數(shù)據(jù)，分別優(yōu)化UPSO-Kriging插值中的球面、指數(shù)以及高斯模型，并與PSO-Kriging以及Kriging的插值預(yù)測(cè)效果進(jìn)行比較，其中PSO算法相關(guān)參數(shù)設(shè)置如表2所示。變差函數(shù)模型參數(shù)的邊界區(qū)間如表3所示。

表2 PSO算法的參數(shù)設(shè)置

表3 變差參數(shù)的邊界設(shè)置

將實(shí)際采樣點(diǎn)高程與預(yù)測(cè)點(diǎn)之間的均方根誤差(RMSE)和平均絕對(duì)百分誤差(MAE)作為算法評(píng)價(jià)指標(biāo)進(jìn)行分析，其中MAE和RMSE公式分別為：

(7)

(8)

由表4可知，在UPSO-Kriging插值方法下將指數(shù)模型作為變差函數(shù)理論模型時(shí)，效果最佳。結(jié)合69組初始樣本點(diǎn)數(shù)據(jù)繪制高程值分布折線A；然后分別根據(jù)UPSO-Kriging算法、PSO-Kriging算法計(jì)算出待預(yù)測(cè)點(diǎn)的高程值，繪制出相應(yīng)高程預(yù)測(cè)折線B及折線C，最終實(shí)驗(yàn)對(duì)比圖如圖5所示。

由表4和圖5可以看出， UPSO-Kriging算法在各種評(píng)價(jià)指標(biāo)下插值的結(jié)果都優(yōu)于其他算法。因此，基于UPSO-Kriging算法的工作面煤層高程值預(yù)測(cè)完全可行。

表4 估算精度比較

圖5 采樣點(diǎn)高程分布對(duì)比實(shí)驗(yàn)Fig.5 Broken line diagram of elevation distribution of sampling points

4 工作面煤層三維模型

4.1 基于規(guī)則網(wǎng)格的DEM模型建立

數(shù)字高程模型是利用有限的地形數(shù)據(jù)對(duì)地形進(jìn)行數(shù)值模擬，通常分為規(guī)則網(wǎng)格DEM模型、TIN模型、等高線模型等[17]。本文結(jié)合煤層實(shí)際地形特點(diǎn)，基于matlab仿真實(shí)現(xiàn)的單層DEM高程模型以及多層DEM高程模型分別如圖6,7所示。

圖6 單層DEM模型Fig.6 Single layer DEM model

4.2 基于Three.js的三維模型建立

Three.js 是一種基于javascript編寫的第三方WEBGL庫(kù)[18]。由于其速度快、易用性強(qiáng)、易于交互等特點(diǎn)，近年來(lái)被廣泛應(yīng)用于三維建?？梢暬I(lǐng)域。本文運(yùn)用Three.js平臺(tái)實(shí)現(xiàn)瀏覽器端三維模型建立，并實(shí)現(xiàn)人機(jī)交互操作。

人機(jī)交互功能主要采用Three.js 中的THREE. OrbitControls控件實(shí)現(xiàn)，在這種實(shí)現(xiàn)模式下，用戶可以通過(guò)鼠標(biāo)對(duì)界面圖形進(jìn)行放大、縮小、旋轉(zhuǎn)、平移等相關(guān)操作，從不同的角度觀察工作面煤層分布情況。如圖8所示，其中黃色為矸石，深色為煤層。

圖7 多層DEM模型Fig.7 Multi-layer DEM model

圖8 人機(jī)交互Fig.8 Human-computer interaction

5 結(jié)論

論文以神東集團(tuán)大柳塔煤礦604工作面的實(shí)際數(shù)據(jù)，運(yùn)用本文提出的UPSO-Kriging插值法對(duì)工作面分層結(jié)構(gòu)的高程進(jìn)行預(yù)測(cè)，進(jìn)而基于規(guī)則網(wǎng)格法構(gòu)建DEM數(shù)字高程模型，以Three.js為平臺(tái)實(shí)現(xiàn)了瀏覽器端綜采工作面三維可視化，并與傳統(tǒng)的Kriging算法和PSO-Kriging算法進(jìn)行對(duì)比，結(jié)果表明，利用UPSO-Kriging方法求得的空間插值數(shù)據(jù)更能反映實(shí)際的地質(zhì)特性，能夠比較準(zhǔn)確地反映工作面地質(zhì)結(jié)構(gòu)。從實(shí)用意義上說(shuō),工作面煤層賦存狀態(tài)三維可視化有著十分重要的意義，可以為煤炭企業(yè)的透明開采、智能開采、分質(zhì)開采提供依據(jù)。