一種邊優先的露天煤礦DEM構建算法

趙景昌,白潤才,劉 威,劉光偉

(1.遼寧工程技術大學礦業學院,遼寧阜新 123000;2.遼寧工程技術大學理學院,遼寧阜新 123000)

?

一種邊優先的露天煤礦DEM構建算法

趙景昌1,白潤才1,劉 威2,劉光偉1

(1.遼寧工程技術大學礦業學院,遼寧阜新 123000;2.遼寧工程技術大學理學院,遼寧阜新 123000)

摘 要:DEM(Digital Elevation Model,數字高程模型)是以高程為地形屬性,用來描述地形表面形態屬性信息的數字表達模型。根據數字化露天煤礦DEM構建需求,提出了一種邊優先的DEM構建算法。該算法基于二維空間直線體素遍歷原理建立邊的格網索引,并采用“最小外接矩形”法搜索DT點。為了確保DEM對地形描述的正確性,采用改進的平坦區域搜索修正算法對相鄰等值線之間的平坦區域進行修正。實驗及應用結果表明,該算法時間效率高,運行穩定,建模精度可靠,能夠滿足露天煤礦DEM構建需要。

關鍵詞:DEM;邊優先;體素遍歷;最小外接矩形;DT點搜索;平坦區域修正

責任編輯:許書閣

趙景昌,白潤才,劉 威,等.一種邊優先的露天煤礦DEM構建算法[J].煤炭學報,2015,40(8):1827-1833.doi:10.13225/ j.cnki.jccs.2014.1026

DEM(Digital Elevation Model,數字高程模型)是以高程為地形屬性,用來描述地形表面形態屬性信息的數字表達模型。DEM最主要的3種表示模型是:規則格網模型、等高線模型和不規則三角網(triangulated irregular network,TIN)模型[1-2]。其中,TIN通過從不規則分布的數據點生成連續三角面來逼近地形表面,與格網模型相比,TIN模型在某一特定分辨率下能用更少的空間和時間更精確地表示更加復雜的表面,因此可以認為TIN為目前進行精細地形表達最好和最常用的方法[3]。對于TIN模型,其基本要求有3點:TIN是惟一的;力求最佳的三角形幾何形狀,即每個三角形盡量接近等邊三角形;保證最鄰近的點構成三角形,即三角形的邊長之和最小。在所有可能的三角網中,Delaunay三角網在地形擬合方面表現最為出色,因此常被用于TIN的生成[4]。

在數字化露天煤礦過程中,基于地質與地形勘測數據、采場與排土場測量驗收數據以及規劃設計數據等構建DEM為礦區時空數據庫建立、生產計劃編制設計、采剝工程量計算、礦區虛擬現實等綜合應用的基礎[5]。構建露天煤礦DEM就是對各種不同來源的建模數據進行三角剖分,在實際建模數據中,除離散點外,經常還存在大量約束邊,如:建模邊界、山脊線、山谷線、斷層線、臺階線坡頂線與坡底線等,對含有約束邊的建模數據進行三角剖分時,應保持其原有的約束關系,否則生成的DEM難以滿足實際應用需要。

對含有約束邊的建模數據進行Delaunay三角剖分即為約束Delaunay三角剖分[6-7]( Constrained Delaunay Triangulation,CDT)。CDT算法主要包括以下5類:分治算法、加密算法、邊優先生長算法、兩步法、掃描線算法[8-12]等。分治算法采用遞歸分塊策略使算法復雜度接近線性,但由于在子塊合并時需要將相鄰子塊兩側凸包邊界采用自下而上的連接算法,涉及大量復雜的判斷,使浮點數誤差錯誤幾率增大,算法穩定性較差;兩步法首先對約束數據集建立非約束Delaunay三角網(初始三角網),然后嵌入約束線段,在約束域較簡單情況下,算法效率很高,但當約束域存在洞、島嶼或約束邊界比較復雜時,該算法在去除多余三角形方面較繁瑣且不夠穩定;加密算法由于加密點的存在破壞了原始數據集,且該算法前期數據處理較煩瑣,使算法整體效率降低;掃描線算法首先將平面多邊形域的所有頂點按掃描方式排列,然后按有序點尋找局部區域的一條邊來構造新的三角形,但在實際應用過程中,為滿足設定的優化標準,該算法需要進行邊翻轉等操作,且不適合在大量離散點存在的情況下應用;邊優先生長算法以非約束Delaunay三角網生長算法為基礎,在建網過程中,對約束邊進行優先擴展,通常結合網格索引搜索擴展邊DT點(能與擴展邊構成Delaunay三角形的“第3點”),以實現快速三角剖分。

根據數字化露天煤礦對DEM構建的具體需求,結合露天煤礦DEM建模數據特點,筆者提出了一種改進的邊優先CDT算法,實驗與實際應用證明,該算法時間效率高,運行穩定,建模精度可靠,能夠較好地滿足數字化露天煤礦DEM構建的需要。

1 算法思想與數據結構

1.1 算法思想

本文算法以露天煤礦DEM建模數據中的約束邊作為優先擴展對象,邊界約束邊向內收縮、邊界內約束邊向外擴展,直至所有邊均飽和(除邊界約束邊外,邊鄰接三角形數2個)為止。上述過程中,擴展邊DT點搜索是決定算法效率的關鍵因素之一,而DT點的搜索效率則取決于DT點搜索范圍以及參與DT點可見性判斷的邊集合規模。本文算法采用分塊技術,并基于二維空間直線體素遍歷原理建立邊的格網索引,使每條邊只關聯于其穿越的格網單元,大大減少了參與DT點可見性判斷計算的邊數,同時,采用“最小外接矩形”法,將每條擴展邊DT點的初始搜索范圍限定于該邊最小外接矩形所覆蓋的格網單元內,從而確保算法具有較高的整體效率。

在基于露天煤礦地質與地形屬性等值線構建DEM時,由大量平三角形(3個頂點高程相等的三角形)構成的平坦區域會使DEM對地形與地質實體的空間形態表達失真,筆者基于DEM拓撲關系,根據平坦區域內相鄰三角形所構成四邊形的凸凹性,分別采用頂點插入與對角邊交換法修正平坦區域,提高了DEM對地形表達的精度。

1.2 數據結構

本文算法主要數據對象包括:頂點(Vertex)、邊(Edge)、Delaunay三角形(DT)、約束Delaunay三角網(CDT)、索引格網(Grid)、索引單元格(GridCell) 等,數據結構如圖1所示。

2 邊優先的露天煤礦DEM構建

2.1 建立點與邊格網索引

建立點與邊的格網索引就是將構建DEM的離散數據域劃分為大小相同的單元格,然后把離散點放到其所在單元格中,把邊放到其所穿越的單元格中,每個單元格分別存儲位于其內部的離散點、穿越或位于該單元格內部的邊,從而建立起單元格、離散點、邊之間的索引關系,目的就是使不規則的DEM離散數據分布“規則化”,提高DEM構建過程中DT點的搜索效率。

索引單元格過大或過小都會增加查詢次數而降低算法效率,筆者經多次實驗確定合理的格網單元邊長cellSize為所有約束邊平均邊長的1.3倍。

建立點的格網索引比較簡單:設某點的平面坐標為(x0,y0),建模數據域最小外接矩形左下角點的坐標為(xmin,ymin),格網單元邊長為cellSize,則該點所在索引單元格的位置(即單元格的行坐標rowId與列坐標colId)可按下式計算:

圖1　本文算法主要數據結構Fig.1 Data structures of the algorithm

建立邊的格網索引時,一種較粗略的方法是用邊最小外接矩形來確定邊的索引單元格,如圖2(a)所示,圖中虛線所示矩形是邊AB的最小外接矩形,則該矩形范圍內的16個單元格即為邊AB的索引單元格,而邊AB實際穿越的單元格僅為圖2(b)中陰影所示7個單元格。

圖2　邊格網索引示意Fig.2 Diagram of edge grid index

可見,用邊的最小外接矩形來確定邊的索引單元格范圍過大,在判斷DT點可見性時,會因為存在大量不必要的相交檢測計算而影響算法效率。為了減少DT點可見性判斷時的計算量,提高算法整體效率,本文算法基于二維空間直線體素遍歷原理[13]建立邊的格網索引,使每條邊僅與其實際穿越的單元格建立索引關系,具體步驟如下:

(1)確定主坐標軸。

設直線起點坐標為(x1,y1),終點坐標為(x2, y2),起點與終點x軸方向的距離dx= |x2-x1|,y軸方向的距離dy= |y2-y1|,若dx>dy,則x軸為主坐標軸,否則,y軸為主坐標軸(dx, dy之間存在另外7種關系,限于篇幅,本文僅以dx>dy且x2> x1,y2>y1為例,其他7種情況原理與此相同)。

(2)用式(1)確定邊起點所在單元格。

(3)計算邊穿越的單元格。

如圖3所示,邊起點V1相對其所在單元格左下角頂點沿x軸方向的距離為xs,沿y軸方向的距離為ys,邊與起點單元格右側邊交點V2相對于左下角頂點的高度hy0為

圖3　計算邊索引單元格示意Fig.3 Diagram of calculating the edge indexed grid cells

比較hy0與起點單元格右上角頂點P0高度ey的關系,若hy0>ey,則邊穿越起點單元格的上鄰單元格與右上鄰單元格,否則,邊穿越右鄰單元格。

在到達終點之前,邊與其所穿越單元格的交點坐標在x軸方向以步長遞增,在y軸方向按下式累加計算:

圖3中,單元格坐標用行與列ID表示,起點單元格坐標r0,c0可用式(1)計算,交點V2高度hy0用式(2)計算,起點單元格(r0,c0)右上角頂點P0高度ey,由于hy0ey,則邊穿越當前單元格(r0,c0+1)的上鄰單元格(r0+1,c0+1)與右上鄰單元格(r0+1,c0+2)。

按上述方法計算直至邊終點所在單元格,即可建立起邊與其實際穿越單元格之間的索引關系。

2.2 DT點快速搜索

建立點與邊格網索引的目的是為了提高DT點搜索效率,此外,DT點搜索范圍也是影響DT點搜索效率的關鍵因素之一。陳學工等在文獻[14]中提出了一種基于“最小搜索圓”的DT點搜索算法,將DT點的搜索范圍控制在以當前擴展邊長度1.5倍為半徑的“最小搜索圓”之內,在構建有10 000個離散點的CDTIN時,算法運行時間為2.15 s。

采用“最小外接矩形”法搜索DT點,實驗證明,該方法的時間效率優于“最小搜索圓”法,具體步驟如下:

(1)從邊集合中取出一條非飽和邊(邊鄰接三角形數小于2個)作為當前擴展邊。

(2)用式(1)分別計算當前擴展邊起點與終點單元格,并據此確定擴展邊“最小外接矩形”單元格范圍。

(3)在步驟(2)確定的單元格范圍內搜索可用DT點,若當前擴展邊為建模邊界邊,則根據頂點排列方向(逆時針或順時針)搜索位于邊左側或右側的頂點作為可用DT點;若當前擴展邊不是邊界邊,則可同時搜索邊左側與右側的可用DT點。

在判斷點的可見性時,可按2.1節所述方法動態建立新邊格網索引,并只對新邊索引單元格所關聯的邊進行相交檢測計算。

(4)計算步驟(3)所有可用DT點中與擴展邊構成DT的頂角,并取頂角最大者與當前擴展邊構成DT。

(5)更新頂點與邊的飽和狀態。點是否飽和可通過計算點角判斷,點角是指與點相連的三角形中,以該點為頂點的內角和,當非邊界點的點角為360°、邊界點角與初始點角相等時(邊界點的初始點角為與該點相連的兩條邊界線段的夾角),則該點飽和;邊的飽和狀態則根據邊鄰接的三角形數確定,當邊鄰接的三角形數為2時,則該邊飽和。飽和點與邊不再作為后續三角剖分中的可用DT點與擴展邊,因此,將其從頂點集合與邊集合中動態刪除,以提高后續三角剖分效率。

(6)重復上述步驟,直至邊集合為空。

2.3 局部LOP優化與平坦區域處理

一般三角剖分算法得到的TIN并不能保證是最優的,通常需要LOP優化,但對于約束Delaunay三角網CDTIN而言,由于約束邊的存在,除了一般LOP優化準則外,還應遵循以下準則[15]:

(1)兩個鄰接三角形若不滿足空圓準則應交換對角線,但如果對應對角線是約束邊可不優化;

(2)若三角形的外接圓中不存在同時滿足從三角形3個頂點都理想可見的頂點,則該三角形仍為DT,不需要優化。

受三角形幾何特性與約束邊數據特征限制,在基于地質屬性等值線或地形等高線構建DEM時,會存在由大量平三角形(三角形的3個頂點屬性值相等)構成的平坦區域,如圖4所示。平坦區域的存在使DEM對露天煤礦地形表達失真,因此,為了獲得正確的DEM,除了進行局部LOP優化外,還應對平坦區域進行修正。

圖4　平坦區域示意Fig.4 Diagram of flat region

對于露天煤礦DEM而言,在一條閉合約束線(地質屬性等值線、地形等高線、臺階坡頂或坡底線)內出現平三角形的主要原因是在閉合約束線內缺少離散的特征點,如果以地形等高線或地質屬性等值線為約束構建DEM,可以通過手動增加特征點來解決閉合等高線內平三角形問題;若以臺階坡頂或坡底線為約束構建露天煤礦DEM,在閉合臺階線內出現平三角形表示該閉合臺階線內高程無明顯變化,可認為是對露天煤礦地形的正確描述。因此,本文只考慮對兩條約束線之間的平坦區域進行修正。

平坦區域修正算法主要分為3類:數據概化修正法、特征修正法和平坦區域搜索修正法[16-18]。數據概化修正法通過減少約束線上的采樣點數量并增加采樣點之間的距離,從而最大限度地減少平三角形的出現,此方法損失了采樣數據信息量,同時也不能完全消除平三角形;特征修正法主要是通過插入特征點、特征線消除平三角形,然而特征線的提取是一個比較復雜的問題,而且此類算法需要插入大量特征點,導致三角網重構工作量較大,修正效率較低;平坦區域搜索修正法[18]首先搜索約束線間存在的平坦區域,然后根據相鄰三角形所構成四邊形的凸凹性,采用插入點或交換邊實現平坦區域的處理,此類算法在修正平坦區域時,在遇到內部平三角形(3條邊都不是約束線的平三角形)后會出現路徑二義性問題,對算法修正速度有一定影響。本文對平坦區域搜索修正法進行了改進,通過引入緩存棧,較好地解決了搜索路徑二義性而導致的算法效率問題。

如圖5(a)中陰影所示平坦區域(其中,tv為與平坦區域鄰接的內部非平三角形;t1~t7為平三角形;t5為內部平三角形),應用本文改進的平坦區域搜索修正法進行修正的具體步驟如下:

圖5　平坦區域修正示意Fig.5 Diagram of flat region correcting

(1)內部非平三角形tv與平三角形t1構成凹四邊形,因此,在tv與t1公共邊上插入中點(該點高程用距離冪次反比法計算),將凹四邊形分裂為4個三角形,修正后DEM如圖5(b)所示;

(2)分裂后的4個三角形中,與t2鄰接的三角形為與殘留平坦區域鄰接的內部非平三角形tv,tv與t2構成凸四邊形,因此,交換t1與tv公共邊,修正后DEM如圖5(c)所示;

(3)圖5(d),(e)修正方法同步驟(2);

(4)圖5(e)中與tv鄰接的三角形t5為內部平三角形,交換tv與t5公共邊后得到兩個分別與平三角形t6,t7鄰接的內部非平三角形tv與tv′(存在二義性搜索路徑),此時建立一個內部非平三角形緩存棧, 將tv與tv′壓入緩存棧;

(5)依次彈出緩存棧棧頂三角形tv′與tv,交換tv′與t7,tv與t6公共邊,結果如圖5(g),(h)所示,此時平坦區域所有平三角形已全部修正,算法結束。

3 算法實驗與應用實例

3.1 算法實驗

本文算法已基于Visual Studio.NET平臺用C#語言實現,表1所列為本文算法與文獻[3,14]算法在相同實驗環境(Windows 7中文旗艦版SP1操作系統,Intel(R) Core(TM) i7-2720QM CPU 2.20 GHz, 4 GB DDR3內存)下,用5組數據構建DEM的運行時間(其中第4、第5組數據是基于地形等高線的DEM數據,無離散的特征高程點)。

表1　算法運行時間Table 1 Running time of the algorithm

由表1中數據可知,在不同特征、不同規模的數據條件下,本文算法時間效率均高于文獻[3,14]算法:其中文獻[3]算法與本文算法DT點搜索范圍是一致的,但在建立邊格網索引時,文獻[3]為根據邊最小外接矩形來確定邊的索引單元格范圍,在后續建網過程中判斷DT點相對于擴展邊的可見性時,由于存在大量不必要的相交檢測計算而影響了算法整體效率;文獻[14]算法基于“最小搜索圓”確定DT點搜索范圍,較本文“最小外接矩形”法范圍大,DT點搜索時參與計算的點與邊數較多、計算量較大而導致算法效率較低。

3.2 應用實例

本文算法已應用在筆者開發的露天煤礦剝采計劃CAD軟件系統中,圖6為應用本文算法生成的地形DEM二維視圖,圖7為應用本文算法生成的露天煤礦現勢DEM與采場計劃DEM三維渲染圖。

圖6　地形DEM二維視圖Fig.6 2D view of terrain DEM

圖7　露天煤礦現勢DEM與計劃DEM三維渲染圖Fig.7 3D rendering of open-pit coal mine’s present situation and plan DEM

實際應用表明,本文算法運行穩定,時間效率高,建模精度可靠,可廣泛應用于各種數據條件下露天煤礦DEM構建。

4 結 語

DEM構建算法是數字化露天煤礦重要的基礎算法之一。本文算法基于二維直線體素遍歷原理建立邊格網索引,并采用“最小外接矩形”法進行DT點快速搜索,此外,針對基于等值線構建DEM時存在大量平三角形而影響模型精度問題,提出了一種改進的等值線間平坦區域修正算法。實驗及應用結果表明,本文算法時間效率高,運行穩定,建模精度可靠,能夠滿足露天煤礦DEM構建需要。

隨著多核與網絡時代的到來,基于多核的多線程并發編程以及分布式并行編程將是未來算法性能提升的主要途徑,筆者將繼續深入研究多核條件下的DEM并行構建算法,以滿足海量數據露天煤礦DEM快速構建的需要。

參考文獻:

[1]鄔 倫,劉 瑜,張 晶,等.地理信息系統:原理、方法和應用[M].北京:科學出版社,2001:195-200.

Wu Lun,Liu Yu,Zhang Jing,et al.Geographic information system: Theory,method and application[M].Beijing:Science Press,2001: 195-200.

[2]胡金星,潘 懋,馬照亭,等.高效構建Delaunay三角網數字地形模型算法研究[J].北京大學學報(自然科學版),2003, 39(5):736-741.

Hu Jinxing,Pan Mao,Ma Zhaoting,et al.Study on faster algorithm for constructing Delaunay triangulations DTM[J].Acta Scientiarum Naturalium Universitatis Pekinensis,2003,39(5):736-741.

[3]湯國安,劉學軍,閭國年.數字高程模型及地學分析的原理與方法[M].北京:科學出版社,2006:107-110.

Tang Guo’an,Liu Xuejun,Lü Guonian.Principles and methods of digital elevation model and analysis[M].Beijing:Science Press, 2006:107-110.

[4]李志林,朱 慶.數字高程模型[M].武漢:武漢測繪科技大學出版社,2000:34-35.

Li Zhilin,Zhu Qing.Digital elevation model[M].Wuhan:Wuhan Technical University of Surveying and Mapping Press,2000:34-35.

[5]隋 心,徐愛功,宋偉東.露天礦精細DEM模型的建立及更新方法[J].測繪科學,2013,38(3):148-150.

Sui Xin,Xu Aigong,Song Weidong.Establishment and update of fine DEM for strip mine[J].Science of Surveying and Mapping,2013, 38(3):148-150.

[6]劉學軍,龔健雅.約束數據域的Delaunay三角剖分與修改算法[J].測繪學報,2001,30(1):83-88.

Liu Xuejun,Gong Jianya.Delaunay triangulation of constrained data set[J].Acta Geodaetica et Cartographica Sinica,2001,30(1):83-88.

[7]Floriani L D.An online algorithm for constrained Delaunay triangulation[J].Graphical Models and Image Processing,1992,54 (3): 290-300.

[8]Chew L P.Constrained Delaunay triangulations[J].Algorithmica, 1989,4(1):97-108.

[9]Boissinnat J D.Shape reconstruction from planar cross sections[J].Computer Vision,Graphics,and Image Processing,1988,44(1):1-29.

[10]Piegl L A.Algorithm and data structure for triangulation multiply connected polygonal domains[J].Computer and Graphics,1993,17(5):563-574.

[11]Sloan S W.A fast algorithm for generating constrained Delaunay triangulations[J].Computers & Structures,1993,47(3):441-450.

[12]Domiter V,?alik B.Sweep-line algorithm for constrained Delaunay triangulation[J].International Journal of Geographical Information Science,2008,22(4):449-462.

[13]劉勇奎,云 健,王曉強,等.沿三維直線的非單位體素遍歷的多步整數算法[J].計算機輔助設計與圖形學學報,2006, 18(6):812-818.

Liu Yongkui,Yun Jian,Wang Xiaoqiang,et al.A multi-step integer algorithm for non-unit voxel traversing along a 3D line[J].Journal of Computer-Aided Design & Computer Graphics,2006,18(6): 812-818.

[14]陳學工,馬金金,黃 偉,等.一種基于最小搜索圓平面多邊形域約束Delaunay三角剖分算法[J].小型微型計算機系統, 2011,32(2):374-377.

Chen Xuegong,Ma Jinjin,Huang Wei,et al.Constrained delauny triangulation algorithm for planar polygonal domains based on minimum search circle[J].Journal of Chinese Computer Systems,2011,32(2):374-377.

[15]周曉云,劉慎權.實現約束Delaunay三角剖分的健壯算法[J].計算機學報,1996,19(8):615-624.

Zhou Xiaoyun,Liu Shenquan.A robust algorithm for constrained Delaunay triangulation[J].Chinese Journal of Computers,1996, 19(8):615-624.

[16]Mark Ware J.A procedure for automatically correcting invalid flat triangles occurring in triangulated contour data[J].Computer & Geosciences,1998,24(2):141-150.

[17]陳學工,黃晶晶.基于等高線建立的TIN中平坦區域的修正算法[J].計算機應用,2007,27(7):1644-1646.

Chen Xuegong,Huang Jingjing.Corrective algorithm of flat areas occurring in TIN constructed from contours[J].Computer Applications,2007,27(7):1644-1646.

[18]解愉嘉,劉學軍,胡加佩.無平三角形處理的等高線數據三角化方法[J].南京師大學報(自然科學版),2012,35(4):106-111.

Xie Yujia,Liu Xuejun,Hu Jiapei.Triangulating the contour data without flat triangle treatment[J].Journal of Nanjing Normal University (Natural Science Edition),2012,35(4):106-111.

Zhao Jingchang,Bai Runcai,Liu Wei,et al.An edge-prior DEM construction algorithm for open-pit coal mine[J].Journal of China Coal Society,2015,40(8):1827-1833.doi:10.13225/ j.cnki.jccs.2014.1026

An edge-prior DEM construction algorithm for open-pit coal mine

ZHAO Jing-chang1,BAI Run-cai1,LIU Wei2,LIU Guang-wei1

(1.School of Mining,Liaoning Technical University,Fuxin 123000,China;2.College of Science,Liaoning Technical University,Fuxin 123000,China)

Abstract:DEM(Digital Elevation Model) is a digital representation model of the terrain surface morphology property information using the elevation as the terrain property.According to the demand for DEM construction in digitalizing open-pit coal mine,an edge-prior DEM construction algorithm was put forward in this paper.The edge grid index was established based on two dimensional linear voxel traversal principle,and the“minimal circumscribed rectangle”method was adopted to search the DT point.In addition,in order to ensure the true terrain description of DEM,the flat regions between adjacent contour lines were corrected with the improved flat region searching and correcting algorithm.Test and application results show that the algorithm is of high efficiency,stable,reliable and accurate,and can meet the requirements of open-pit coal mine DEM construction.

Key words:DEM;edge-prior;voxel traversal;minimal circumscribed rectangle;DT point searching;flat region correcting

作者簡介:趙景昌(1974—),男,內蒙古寧城人,講師,博士研究生。E-mail:lntuzjc@126.com

基金項目:國家自然科學基金資助項目(51304104);遼寧省教育廳科學研究一般資助項目(L2011051)

收稿日期:2014-08-06

中圖分類號:TD176

文獻標志碼:A

文章編號:0253-9993(2015)08-1827-07