張小艷, 朱圣凱, 楊鑫磊

(西安科技大學計算機科學與技術學院，西安 710000)

傳統(tǒng)的二維圖表在地質構造規(guī)律及拓撲關系表達上受限于其平面化表達方式，存在著空間信息缺失、制圖過程復雜的缺陷，越來越無法滿足人們對地層認識和空間分析的需求。而三維模型由于其直觀的空間關系、可視化的操作分析等優(yōu)點，在表達礦山地質體的地層信息及空間構造特征方面具有巨大優(yōu)勢[1-4]。

近年來隨著計算機技術的發(fā)展，以及數字礦山、智慧礦山等概念的提出，礦山地質體三維建模受到越來越多的關注[5]。中外學者對礦山地質三維建模的進行了大量的研究，吳志春等[6]利用PRB數字填圖技術在GOCAD軟件平臺上構建淺部三維地質模型，所建立的地質界面與實際地質數據較為吻合，但所建三維模型難以在其他軟件平臺上展示，移植性較差；Yucel等[7]借助于Agisoft軟件進行無人機圖像的3D建模，實現(xiàn)了Etili Comakli露天煤礦三維地形模型的構建與可視化，高效地完成了三維地質建模過程中涉及的圖像匹配、地理配準和數字高程建模過程；堅潤堂等[8]基于DIMINE軟件，采用距離冪次反比法對礦體進行品位推估并建立西藏白容—崗講礦區(qū)礦床三維實體模型；魏永勇[9]采用C/C++及IDL混合編程，借助于二次曲面Shepard插值法構建了地層的數字高程模型，實現(xiàn)了三維地層建模，但所建模型中地層交界線缺乏必要的曲線平滑，預測精度需進一步提高；此外楊彪等[10]借助Surpac礦業(yè)軟件結合VRML和HTML語言, 構建礦山地表、地層、斷層等地質模型，實現(xiàn)了礦山在Web環(huán)境下三維可視化與交互。

綜上，已有礦山地質三維建模的研究雖各有特色，但也有其缺陷和不足，借助ARC/INFO、3DMine等C/S結構建模軟件實現(xiàn)的礦體三維模型難以跨平臺使用，無法與礦山企業(yè)現(xiàn)有的Web管理系統(tǒng)相結合；而利用VRML、Flash等Web 3D[11]技術實現(xiàn)的礦體三維模型則需要額外加載瀏覽器插件或者虛擬機才能完成礦山地質體模型的渲染，限制了礦體三維模型在Web端的可視化；此外建模過程中采用單一空間插值預測算法難以處理復雜的地質構造，無法實現(xiàn)切合實際地質情況的三維地質模型?；谝陨戏治觯靡环N可在變異過程中修正變異方向的自適應差分進化算法來優(yōu)化Kriging空間插值算法(updating and modifying adaptive differential evolution Kriging algorithm, UMDE-Kriging)，結合規(guī)則格網法建立精確度更高的數字高程模型(digital elevation model, DEM)，并在此基礎上研究基于Web GL標準、以瀏覽器為載體的三維數字化建模方式，利用Three.js和Catmull-Rom樣條曲線實現(xiàn)無插件、可移植、跨平臺、支持多瀏覽器運行的采煤工作面煤層三維可視化模型，為煤礦采煤工作面的開采掘進提供重要參考。

1 關鍵技術

1.1 Web GL

Web GL(Web graphics library)是一種跨平臺的3D繪圖規(guī)范，通過在OpenGL ES 2.0中添加JavaScript綁定，從而為HTML5中Canvas提供硬件3D加速渲染[12]。與早期的Web 3D實現(xiàn)相比，Web GL能夠在大多數平臺上運行，兼容性好，同時還支持CPU硬件加速，實現(xiàn)三維圖形渲染。

1.2 Three.js介紹

Three.js依據Web GL規(guī)范，對底層Web GL代碼進行簡單封裝，可以將二維數據以三維方式呈現(xiàn)出來，是一款開源的Web GL圖形開發(fā)框架。簡單、直觀的封裝Web GL常用的三維對象，通過掩飾細節(jié)來加快開發(fā)效率，減輕開發(fā)者的開發(fā)負擔[13]。由于Three.js采用了JavaScript語言，因此能夠很容易的與其他瀏覽器組件進行交互。

Three.js與其他的Web GL圖形框架一樣，在使用時都需先定義三維世界中所包含的基本元素，主要流程如下：

(1)設置三維場景(scene)：場景是用來存放一切所創(chuàng)建對象的容器，用THREE.Scene來表示。

(2)相機設置(camera)：相機決定了顯示場景中哪一個角度的景色。

(3)設置燈光源(light)：Three.js中提供一系列的光源，包括環(huán)境光、方向光、點光源等。

(4)設置物體(object)：物體就是場景中需要處理的三維對象。Three.js通過加載3D建模工具導出的模型文件或者借助THREE.mesh三維網格動態(tài)創(chuàng)建三維場景中的模型。

(5)添加渲染器(renderer)：渲染器決定了3D空間中的物體對象映射到二維平面上的方式，以及渲染的結果應該繪制到場景中什么元素上。

2 自適應差分進化Kriging算法

2.1 Kriging算法

Kriging插值法又稱空間自協(xié)方差最佳插值法，是依據協(xié)方差函數對隨機過程/隨機場進行空間建模和預測(插值)的回歸算法[14]。根據采樣點與預測點之間的空間位置和相關程度，確定一個影響待插點值的距離范圍。對其中每一個樣品賦予一定的權重，最后進行加權平均來估計待預測點屬性值。

Kriging插值法的一般公式為

(1)

式(1)中：Z*(x)為待求插值點的估計值；Z(xi)為第i個采樣點xi處的數值；λi是第i個采樣點的插值權重系數，并且需滿足無偏性和最優(yōu)性兩個條件。即：

(2)

σ2=E[Z*(x)-Z(x)]2=min

(3)

式中：σ為待求插值點的估計值Z*(x)和真實值Z(x)之間的標準差，在Kriging法中，用來衡量各個樣本點之間空間相關程度的是變差函數。

(4)

式(4)中：h為兩點之間距離；N是由h分開的成對樣本點的數量；Z(xi)是點的屬性值。當滿足本征假設時，可用變差函數來表示Kriging方程組：

(5)

式(5)中：γ(xi，xj)表示兩個采樣點xi和xj的變差函數值；μ為拉格朗日系數；j=1，2,…，n。

2.2 差分進化算法

差分進化算法(DE)思想源于早期提出的遺傳算法(GA)，模擬遺傳學中的雜交、變異、復制來設計遺傳算子，通過優(yōu)化具體參數、值函數解決連續(xù)的非線性優(yōu)化問題[15]。

使用DE算法中的變異、交叉和選擇運算實現(xiàn)種群的進化，其中差分變異公式為

vi(g+1)=xr1(g)+F·[xr2(g)-xr3(g)]，

i=1,2,…,NP, 1≤r1≠r2≠r3≤NP

(6)

式(6)中：vi(g+1)為g+1代變異個體；xr1(g)、xr2(g)、xr3(g)是從t代種群中隨機選取的三個個體；F為縮放比例因子；NP表示種群大小。

由g+1代變異個體vi(g+1)與上一代個體xi(g)交叉得到g+1代子個體ui(g+1)，交叉公式為

(7)

式(7)中，當rand(0,1)函數產生的隨機數小于交叉概率RC時，ui(g+1)中對應參數來自變異個體vi(g+1)，否則該參數來自上一代個體xi(g)。

DE算法采用貪婪算法來選擇進入下一代種群的個體，選擇公式為

(8)

式(8)中，f為待優(yōu)化求解的目標函數。

2.3 UMDE-Kriging算法實現(xiàn)

基于差分進化算法以及Kriging插值算法，結合煤礦采煤工作面實際地質情況以及礦井測量人員經驗公式，設計合適的自適應放縮因子F，通過差分、變異、選擇建立變差函數模型參數最優(yōu)解，進而求解出滿足精度的待預測點高程屬性值。

UMDE-Kriging算法主要步驟如下：

2.3.1 確定變差函數參數

通過式(4)得出實驗變差函數云圖，選定合適的變差函數模型，并確定其參數c0、c、a，其中c0為塊金值，c為偏基臺值，a為變程。

2.3.2 選定合適目標函數

變差函數模型的擬合本質上是表征實驗變差函數與理論變差函數的近似程度。因此，可以將各距離下的實驗變差函數值與理論變差函數值之間的差值平方和作為目標函數。

(9)

式(9)中，γ*(h)為理論變差函數值；γ(h)為實驗變差函數值；m為實驗變差函數值的個數。

2.3.3 初始化種群

在解空間隨機生成M個個體，因為解空間是模型參數c0、c、a的集合域，即每個個體由三維向量組成。

(10)

每個個體中的參數取值如下：

(11)

式(11)中，(C0_max-C0_min)、(Cmax-Cmin)、(Amax-Amin)分別為模型參數c0、c、a的取值區(qū)間，需要人為設定。

2.3.4 確定差分進化算法的有關運行參數

在差分進化算法的變異過程中通過增加自適應縮放因子F，修正變異方向。確定最大進化代數、種群大小和交叉概率。

(12)

式(12)中，F(xiàn)0為設置的縮放因子；fi(pre)為當前個體的上一代個體適應值；fmin(pre)為上一代最優(yōu)個體適應值；fi(cur)為當代個體的適應值。

2.3.5 差分進化運算

在隨機生成的初始種群基礎上，進行變異[式(6)]、交叉[式(7)]和選擇[式(8)]三個過程，多輪進化迭代直至得出最優(yōu)解fmin，求得變差函數模型參數，擬合出相應的模型表達式。

2.3.6 預測值計算

將得到的變差函數模型表達式代入式(5)中求解，得出采樣點權重系數λi，進而通過式(1)得到待預測點屬性值。

3 數據準備

選取某煤礦16402采煤工作面0～700 m深度范圍內的54組巷道鉆孔數據以及切眼地質采樣數據作為樣本數據集合(包括22組巷道鉆孔數據和32組切眼采樣數據)，如表1所示。

表1 采樣點數據

所選工作面區(qū)域橫向寬度為200 m，鉆孔采樣相對煤層底板的高程范圍為0～2 m，樣本點分布如圖1所示。

圖1 采樣點分布

3.1 算法可行性驗證

從樣本數據集合中提取出同屬于分一煤層下表面的控制點數據，選取Kriging插值法中的指數模型作為優(yōu)化對象，利用修正變異方向的自適應差分進化算法(UMDE)進行優(yōu)化。其中，差分進化算法的參數設置如表2所示，模型參數邊界區(qū)間如表3所示。

表2 差分進化算法的參數設置

表3 模型參數的邊界區(qū)間

首先根據54組初始樣本點數據按照采樣編號繪制出高程值分布折線A；之后分別根據UMDE-Kriging算法和普通Kriging算法計算出相應位置的高程理論值，繪制出相應折線B和折線C，如圖2所示。

圖2 采樣點高程分布折線圖

根據54組預測結果和實際高程值進行交叉驗證，分別計算均方根誤差(RMSE)和平均絕對百分誤差(MAPE)。

(13)

(14)

表4 算法預測精度對比

由圖2和表4可知，UMDE-Kriging算法插值的結果要比普通Kriging算法更加接近于與實際采樣點高程值，誤差更小，在高程預測求解過程中表現(xiàn)出更好的有效性和準確性。因此，將UMDE-Kriging插值算法應用于采煤工作面信息預測完全可行。

3.2 數據格網化處理

數字高程模型(DEM)是通過有限的地形高程數據實現(xiàn)地形表面的數字模擬或表面地形的數字表示[16]，其中規(guī)則格網DEM以其數據結構簡單、運算速度快、效率高的優(yōu)點大量應用于實際工程。

根據采煤工作面邊界約束，利用規(guī)則格網法生成各地質分層的格網表面，并結合UMDE-Kriging空間插值算法進行任意格網點的屬性值預測和采集，具體步驟如下：

3.2.1 地層劃分

針對16402采煤工作面涉及的煤層地質數據進行綜合分析，將相同的巖性并且垂直方向位置相似的看作同一層，根據地層分界控制點信息將其劃分為分一煤層、夾矸層、分二煤層。

3.2.2 格網劃分

分別計算采煤工作面各地層起始深度和終止深度的采樣點水平坐標(x，y)的最大值和最小值，合理劃分網格單元大小，計算每個格網點的x、y坐標。

3.2.3 格網點插值

利用UMDE-Kriging空間插值算法對各地層格網點的高程、煤質指標進行插值擬合處理，生成具有屬性值的規(guī)則格網模型。

將16402采煤工作面0～700 m范圍區(qū)域劃分為11×11格網，以分一煤下表面為例，結合UMDE-Kriging算法建立規(guī)則格網模型，如圖3所示。

圖3 規(guī)則格網模型

4 煤層三維模型實現(xiàn)

煤層是空間中的一種層狀地質實體，其空間形態(tài)與分布相對于非層狀實體較為簡單，具有較強的全局規(guī)律性[17]，針對采煤工作面的煤層地質情況，可將其分為開采煤分層和夾矸層。

以16402采煤工作面54組樣本數據為基礎，采用JetBrains WebStorm 2018作為開發(fā)工具，利用JavaScript與Three.js進行開發(fā)。借助于規(guī)則格網法和UMDE-Kriging插值法形成規(guī)則網格數據模型，根據點聚成線，線聚成面，面聚成體的原則逐步實現(xiàn)采煤分層模型與夾矸層模型；建立采煤工作面0～700 m范圍內煤層三維體元模型，實現(xiàn)部分煤質信息的動態(tài)展示。其中，所建模型長度為700 m，寬度為200 m，高度為2 m。

4.1 采樣點分布模型

將各煤層規(guī)則網格數據模型中每個數據點視為一個獨立的球體模型，依據其三維坐標屬性值(x，y，z)確定其空間位置，結合THREE.SphereGeometry()構造方法循環(huán)完成各煤分層和夾矸層所有數據點的繪制。其中點半徑選取為30，球體維度方向和經度方向上的分段數均設置為10：var geometry=new THREE.SphereGeometry(30, 10, 10)，采樣點分布模型如圖4(a)所示。

4.2 煤層面平滑曲線模型

Catmull-Rom樣條曲線是一種3次的分段插值樣條，它利用前后的2個節(jié)點直接計算出各點處的切矢量，每相鄰4個節(jié)點的信息就可以構造出一段經過中間兩點的3次多項式函數，并且所構造的每一段樣條能夠保持連續(xù)[18]。Catmull-Rom曲線每個分段擁有四個控制點，假定這四個控制點分別為P0、P1、P2、P3。則Catmull-Rom樣條函數的矩陣表達式為

(15)

式(15)中，浮點坐標t∈[0,1]，其中M為樣條曲線插值矩陣：

(16)

式(16)中，T稱為形狀因子，當樣條用于離散點連續(xù)平滑時，取T=0.5為最佳值。根據式(15)和式(16)可以在P1、P2之間構造出一段平滑的曲線，經過拼接每一段平滑曲線可構成經過所有控制點的完整連續(xù)曲線。

利用所建立的規(guī)則格網模型，分別存儲其中同屬一條線上的控制點信息，借助Catmull-Rom插值算法繪制Catmull-Rom樣條曲線，建立煤層面光滑曲線模型，為煤層面曲面模型繪制打下基礎。如圖4(b)所示，以分一煤層下表面格網模型數據為基礎，分別繪制巷道平行方向和垂直方向Catmull-Rom樣條曲線。

4.3 多煤層面元模型

分層建模的基本思想是按照所劃分的地層順序，根據不同地質層面的控制點依次建立各地層面模型?；谒L制的Catmull-Rom樣條曲線，反算每一條曲線上的點的空間位置信息。依據點在曲線的相對位置，按照弧長比例返回其對應的三維坐標值，并將其作為所屬地質分層曲面控制點；之后借助Three.js自帶的參數化幾何模型THREE.ParametricGeometry依次實現(xiàn)地層面曲面模型的繪制：var geometry=new THREE.ParametricGeometry(function(u,v), slices, stacks)，其中slices用于參數化函數的切片數量u；stacks用于參數化函數的堆棧數量v，循環(huán)調用function()返回所有曲面控制點的三維坐標信息。利用接收到控制點位置信息并繪制各煤層曲面模型，采煤工作面各地質分層(煤層與矸石層)曲面模型如圖4(c)所示。

4.4 體模型繪制

通過獲取采樣點水平方向的空間分布信息和垂直方向的煤層分層信息，在生成多個描述不同煤層地質面模型的基礎上，利用各煤層上下表面邊界線，在這兩條閉合輪廓線之間建立一系列的三角面片，設置當前層的顏色，構建出煤層側面模型；結合各煤層頂板面模型、底板面模型共同組合成封閉立體空間，完成采煤工作面各地質分層模型繪制。如圖4(d)所示，每個地質分層(煤層與矸石層)都是有上下兩個表面模型和若個側面模型“縫合”而成。

圖4 采煤工作面模型

4.5 三維模型交互顯示

4.5.1 人機交互

基于Three.js建立的工作面煤層三維模型與用戶的交互采用THREE.OrbitControls控件實現(xiàn)：controls_camera=new THREE.OrbitControls(camera, renderer.domElement)。

在這種控制模式下，分別利用鼠標左鍵、中鍵、右鍵改變場景中相機的視點，實現(xiàn)三維模型的旋轉、放縮、平移。從不同的角度觀察采煤工作面煤層分布情況，如圖5所示。

圖5 人機交互展示

4.5.2 煤質信息展示

通過給采煤工作面三維模型綁定鼠標響應事件來完成各煤層煤質信息的展示：documenet.addEventListener(“mousemove”,函數,false)。

當鼠標光標移動至三維模型中某個平面時，利用THREE.Raycaster控件完成模型中平面對象的拾取，實現(xiàn)選中煤分層或矸石分層表面采樣點三維坐標和煤質信息展示。如圖6所示，在左側信息顯示區(qū)域展示選中曲面所有控制點的三維坐標以及部分煤質指標。

Td為真密度，Ad為灰分

5 結論

以某礦16402綜采工作面的部分采樣數據為基礎，驗證了UMDE-Kriging空間插值算法可行性，并結合規(guī)則格網法完成虛擬采樣點高程值及煤質指標的插值擬合處理，建立了精確度更高的數字高程模型，并借助Three.js技術和Catmull-Rom樣條曲線構建了能夠反映實際地質特征的16402采煤工作面煤層三維可視化模型。

(1)所建模型克服了基于其他建模技術需要額外加載插件才能完成地質體模型渲染的缺陷，具有可移植、跨平臺、支持多瀏覽器運行的優(yōu)勢。

(2)實現(xiàn)了三維模型的旋轉、平移、放縮和分層顯示，具有良好的人機交互性能。

(3)對工作面的地質空間特性及煤質指標進行了良好的展示，在表面模擬和內部屬性表達兩個方面達到較好的融合。

(4)能夠完整、便捷的嵌入至煤炭企業(yè)Web端信息管理系統(tǒng)中，為礦山開采掘進提供技術支持。

但本模型的構建針對煤層中可能存在的斷層、褶皺等特殊地質區(qū)域還不能很好地呈現(xiàn)，在后續(xù)研究中，將結合煤田地質學以及空間預測算法進行修正，進一步完善。