基于Delaunay三角網(wǎng)礦體和采空區(qū)儲(chǔ)量計(jì)算算法與實(shí)現(xiàn)＊

李 華，寧化展，丁希旺

（濱州市水利勘測(cè)設(shè)計(jì)研究院測(cè)量隊(duì)，山東 濱州256612）

0 引 言

數(shù)字地形模型（DTM）是以離散分布的平面點(diǎn)來(lái)模擬連續(xù)分布的地形，是野外地表勘測(cè)成果的數(shù)字化展現(xiàn)，廣泛的應(yīng)用于地理信息系統(tǒng)各領(lǐng)域中。Delaunay三角網(wǎng)是DTM的主要實(shí)現(xiàn)形式，用一系列互不交叉、重疊的連接在一起的三角形網(wǎng)來(lái)表示地形［1］。

礦產(chǎn)儲(chǔ)量，簡(jiǎn)稱(chēng)儲(chǔ)量，一般是指經(jīng)過(guò)地質(zhì)勘探、基建勘探和生產(chǎn)勘探后，經(jīng)勘查證實(shí)存在礦體，探明其空間分布、產(chǎn)狀、形態(tài)、規(guī)模和質(zhì)量，可以開(kāi)發(fā)利用的原地礦產(chǎn)資源量［2］。根據(jù)各種探礦工程和技術(shù)手段所得到的資料及信息（一般是鉆孔資料），通過(guò)一定的計(jì)算方法計(jì)算礦產(chǎn)的地下埋藏量，這一系列的工作稱(chēng)之為礦產(chǎn)儲(chǔ)量（資源量）計(jì)算。

礦產(chǎn)資源儲(chǔ)量的計(jì)算方法有斷面法、算術(shù)平均法、等值線(xiàn)法等，這些算法都是在礦產(chǎn)資源勘探時(shí)期估算礦產(chǎn)資源儲(chǔ)量的計(jì)算方法，這些算法并不適于計(jì)算礦體采空區(qū)的開(kāi)采量。為此，提出了一種基于礦體區(qū)域的采空區(qū)礦體儲(chǔ)量計(jì)算方法，根據(jù)礦體采空區(qū)范圍及回采率計(jì)算出礦體資源已開(kāi)采量和礦體資源剩余儲(chǔ)量。

1 礦體區(qū)域Delaunay三角網(wǎng)的生成

基本Delaunay三角網(wǎng)生成有以下四個(gè)步驟：凸殼的生成過(guò)程、初始三角網(wǎng)生成過(guò)程、離散點(diǎn)插入初始三角網(wǎng)過(guò)程、局部?jī)?yōu)化過(guò)程。礦產(chǎn)資源的儲(chǔ)量計(jì)算是在已知礦產(chǎn)資源分布邊界，以及區(qū)域鉆孔資料數(shù)據(jù)計(jì)算所得，根據(jù)鉆孔的稀疏程度還要對(duì)離散鉆孔點(diǎn)進(jìn)行插值運(yùn)算。因此論述的Delaunay三角網(wǎng)生成有以下六個(gè)步驟：

1）點(diǎn)差值過(guò)程：由于鉆孔點(diǎn)不多以及區(qū)域邊界太長(zhǎng)，為了使生成的三角網(wǎng)最優(yōu)及礦產(chǎn)量的計(jì)算更準(zhǔn)確需要對(duì)邊界點(diǎn)及鉆孔點(diǎn)進(jìn)行插值運(yùn)算。

2）凸殼生成過(guò)程：凸殼［3］是包含所有離散點(diǎn)的最小多邊形，一般的凸殼構(gòu)造方法只是找出了最少點(diǎn)的多邊形，對(duì)于多點(diǎn)在一條凸殼邊上的情況只是找出了兩個(gè)端點(diǎn)而中間點(diǎn)卻沒(méi)有找出來(lái)，這種情況對(duì)礦體的建模有影響，為此采用一種夾角與距離最?。?］的凸殼構(gòu)造方法。

3）初始三角網(wǎng)生成過(guò)程：利用凸殼點(diǎn)生成初始三角網(wǎng)（如圖1所示）。

4）點(diǎn)插入過(guò)程：將點(diǎn)集中的離散點(diǎn)（不包含凸殼點(diǎn)）依次插入［5］，在初始三角形鏈表中找出三角形外接圓包含插入點(diǎn)的三角形（稱(chēng)為該點(diǎn)的影響三角形），剔除查找出的影響三角形的公共邊，將此插入點(diǎn)同查找到的影響三角形的全部頂點(diǎn)連接起來(lái)，完成一個(gè)點(diǎn)在Delaunay三角形鏈表中的插入，如圖2所示。

5）局部?jī)?yōu)化過(guò)程：局部?jī)?yōu)化算法是為了生成Delaunay三角網(wǎng)。算法的基本含義：對(duì)由兩個(gè)公共邊組成的四邊形進(jìn)行判斷，如果其中一個(gè)三角形的外接圓包含第四個(gè)頂點(diǎn)，則這個(gè)四邊形的對(duì)角線(xiàn)互換，如圖3所示。

圖1 初始三角網(wǎng)

圖2 點(diǎn)插入

圖3 局部有化過(guò)程

6）三角形剔除過(guò)程：根據(jù)已知邊界條件，將不符合條件的三角形剔除，如圖4所示。若只是按照三角形三個(gè)點(diǎn)都在邊界上這個(gè)條件來(lái)剔除可能會(huì)剔除一些合法的三角形，因此還要加上內(nèi)切圓圓心是否在邊界內(nèi)這個(gè)條件。算法描述如下：①判斷三角形的三個(gè)頂點(diǎn)是否位于邊界上，②若果不是都在邊界上刪除三角形，如果均位于邊界上，求出內(nèi)切圓的圓心，判斷是否位于邊界內(nèi)，若是三角形保留否則三角形剔除［4］，如圖4所示。

圖4 三角形剔除

2 采空區(qū)域Delaunay三角網(wǎng)的生成

礦山儲(chǔ)量動(dòng)態(tài)管理中，遇到的主要問(wèn)題是兩個(gè)三維實(shí)體：采空區(qū)與礦體［6］。需要管理的內(nèi)容是二者相交部分的礦產(chǎn)量，此礦產(chǎn)量的計(jì)算首先要求采空區(qū)的平面區(qū)域與礦體Delaunay三角網(wǎng)區(qū)域的相交Delaunay三角網(wǎng)區(qū)域，如圖5所示，再求出此Delaunay三角網(wǎng)區(qū)域的礦產(chǎn)量。

圖5 采空區(qū)域

根據(jù)采空區(qū)平面二維邊界求其與礦體表面Delaunay三角網(wǎng)相交部分的Delaunay三角網(wǎng)的算法思想如下：

1）如圖5所示，內(nèi)邊界點(diǎn)也要參與Delaunay三角網(wǎng)的構(gòu)建；

2）以?xún)?nèi)邊界的一條邊開(kāi)始查找此邊影響的三角形（即：與此邊相交的三角形）并依次對(duì)查找的三角形編號(hào)記錄下來(lái)，判斷三角形在邊界外側(cè)的頂點(diǎn)與此三角形跟影響邊的交點(diǎn)組成的是三角形還是四邊形，若是三角形直接連接，若是四邊形根據(jù)Delaunay三角形優(yōu)化原則連接四邊形的一條對(duì)角線(xiàn)構(gòu)成兩個(gè)三角形，如圖6所示。

影響三角形的查找方法如下：以起始邊的一個(gè)端點(diǎn)開(kāi)始，首先查找以此端點(diǎn)為頂點(diǎn)的三角形并判斷那個(gè)三角形與此邊相交求出交點(diǎn)，如此得到第一個(gè)影響三角形；找出與第一個(gè)三角形影響邊共邊的三角形并求出此三角形與起始邊的交點(diǎn)，如此得到第二個(gè)影響三角形，依次循環(huán)直到求出的交點(diǎn)為起始邊的另一端點(diǎn)，這樣就找出了所有影響三角形及它們與影響邊的交點(diǎn)；

圖6 邊界影響三角形

3）依照步驟2）求出其它邊界邊的影響三角形和每個(gè)影響三角形與邊界邊的交點(diǎn)，并按步驟2）的方法生成邊界外側(cè)的三角網(wǎng)；

4）邊界點(diǎn)和邊界上的交點(diǎn)與邊界內(nèi)的點(diǎn)重新生成Delaunay三角網(wǎng)，如圖7所示。

圖7 三角網(wǎng)生成

3 基于Delaunay三角網(wǎng)的儲(chǔ)量計(jì)算方法

根據(jù)Delaunay三角網(wǎng)模型，直接利用Delaunay三角形進(jìn)行體積計(jì)算，體積計(jì)算如下：先求出每個(gè)Delaunay三角形所構(gòu)成立柱體的體積Vi，再把區(qū)域內(nèi)所有立柱體的體積相加，從而得到整個(gè)礦區(qū)內(nèi)礦產(chǎn)體積儲(chǔ)量V，再乘以礦產(chǎn)密度，就是礦體重。因此，地質(zhì)體的體積計(jì)算也就是計(jì)算每個(gè)三棱柱［7］，如圖8所示的體積，然后累積和就是地質(zhì)體的體積。

圖8 礦體線(xiàn)性模型

已知三棱柱上表面三角形三個(gè)點(diǎn)A，B，C和下表面三角形三個(gè)點(diǎn)D，E，F(xiàn)的（X，Y，Z）坐標(biāo)，求三棱柱的體積算法思想如下：

1）找出上表面三角形點(diǎn)Y值最小的點(diǎn)A和下表面三角形點(diǎn)Y值最大的點(diǎn)E，并以這兩點(diǎn)為基準(zhǔn)做橫切面把三棱柱分成一個(gè)直三棱柱和兩個(gè)四棱錐體。

2）根據(jù)A，B，C三點(diǎn)的坐標(biāo)值可求出BB1、CC1、C1B1、AC1、AB1的長(zhǎng)度，進(jìn)而求出四邊形CC1B1B的面積S.在三角形AB1C1中根據(jù)求出的三邊長(zhǎng)可求出A到邊B1C1的高h(yuǎn)，即四棱錐的高。上表面四棱錐的體積v＝S＊h／3.

3）A，E兩點(diǎn)的高程差即是直三棱柱的H，根據(jù)第二步求出的三角形AB1C1三個(gè)邊長(zhǎng)可求出直三棱柱的橫截面面積S，則直三棱柱的體積V＝S＊H.

4）按照步驟②的方法求出下表面四棱錐的體積，三個(gè)體積累計(jì)求和即為此三棱柱的體積，多個(gè)三棱柱體積的累計(jì)即得到地質(zhì)體的體積。

圖9 礦體面狀模型

圖9、10示出了以滕州曹莊煤礦三層煤數(shù)據(jù)及采空區(qū)數(shù)據(jù)生成的煤層模型和采空區(qū)模型。選中煤層或任一采空區(qū)都可計(jì)算出其體積、表面積及投影面積。體積及面積信息在窗體狀態(tài)欄左下角顯示，如圖11所示。

圖10所示淺色區(qū)域?yàn)檫x中的礦體采空區(qū)區(qū)域，其體積為2 894 621.0m3與ARCGIS計(jì)算的體積2 888 418.6m3相似率達(dá)0.997 8，實(shí)踐證明此算法是正確有效的。

圖10 礦體面狀模型

圖11 狀態(tài)欄

4 結(jié) 論

根據(jù)鉆孔的坐標(biāo)及屬性數(shù)據(jù)庫(kù)，利用java3D可視化技術(shù)，通過(guò)對(duì)Delaunay三角形生成過(guò)程中點(diǎn)集進(jìn)行封閉性判斷，生成任意地區(qū)Delaunay三角網(wǎng)，結(jié)合鉆孔屬性數(shù)據(jù)庫(kù)，對(duì)任意區(qū)域礦產(chǎn)儲(chǔ)量進(jìn)行計(jì)算，為計(jì)算任意區(qū)域礦產(chǎn)儲(chǔ)量與可視化提供了一種新思路，該方法不僅克服算術(shù)平均法誤差大及其它方法難以自動(dòng)實(shí)現(xiàn)的缺點(diǎn)，而且為礦山管理提供可視化查詢(xún)與分析功能。

［1］ 袁 翰，李偉波，陳婷婷.對(duì)構(gòu)建Delaunay三角網(wǎng)中凸殼算法的研究與改進(jìn)［J］.計(jì)算機(jī)工程，2007，33（7）：70－72.

［2］ 蔣 銳，宋煥斌，朱杰勇.基于三維柵格數(shù)據(jù)的露天礦山儲(chǔ)量動(dòng)態(tài)監(jiān)測(cè)研究與應(yīng)用［J］.礦產(chǎn)與地質(zhì)，2009，23（5）：469－472.

［3］ 劉永和，王潤(rùn)懷，齊永安.一種非凸包邊界約束不規(guī)則三角網(wǎng)生成算法［J］.測(cè)繪科學(xué)，2008，33（3）：79－81.

［4］ 寧化展，田茂義.帶有邊界條件的Delaunay三角網(wǎng)生成算法的研究與實(shí)現(xiàn)［J］.全球定位系統(tǒng)，2010，35（4）：39－42.

［5］ 徐道柱，劉海硯.大量約束邊條件下Delaunay三角網(wǎng)的快速生成［J］.測(cè)繪工程，2007，16（3）：6－10.

［6］ 張尹軍.礦產(chǎn)資源儲(chǔ)量動(dòng)態(tài)監(jiān)測(cè)的問(wèn)題及建議［J］.商品與質(zhì)量，2010，（7）：74－75.

［7］ 張渭軍，孔金玲，翁曉鵬.基于鉆孔數(shù)據(jù)Delaunay三角網(wǎng)的礦產(chǎn)儲(chǔ)量計(jì)算［J］.測(cè)繪科學(xué)，2010，35（4）：57－58.