基于 GIS的高質量約束Delaunay三角網格剖分

趙 曉 東,晏 小 寶,沈 永 明,王 亮

(1.大連大學院士創業園中日地層環境科學研究中心,遼寧大連 116622;2.大連理工大學海岸和近海工程國家重點實驗室,遼寧大連 116023）

基于 GIS的高質量約束Delaunay三角網格剖分

趙 曉 東1,晏 小 寶1,沈 永 明2,王 亮2

(1.大連大學院士創業園中日地層環境科學研究中心,遼寧大連 116622;2.大連理工大學海岸和近海工程國家重點實驗室,遼寧大連 116023）

在分析現有非結構化網格剖分算法的基礎上,提出了一種 GIS支持下的改進分治算法實現約束Delaunay三角網格剖分。該方法利用了 GIS的空間拓撲關系對算法輸入數據進行預處理,基于三角形的統一數據結構實現了網格細化,對輸出剖分網格進行準確的拓撲和約束條件的檢查,并基于推進陣面算法思想,結合空間鄰近拓撲關系實現了三角剖分節點和網格的重新編號,方便了實際問題中開邊界條件的賦值,提高了計算效率。實例應用表明,該方法大大簡化了數值模型非結構化網格剖分的前處理過程,集成了幾種綜合算法的優點,在保證原分治算法時間復雜度的基礎上,提高了約束條件下Delaunay三角網格生成的質量。

網格剖分;GIS;約束Delaunay三角剖分

0 引言

在應用有限元、有限差分和有限體積法對力學問題進行數值計算的前處理中,網格自動剖分占有重要的位置。目前,實現網格剖分的算法很多,而且不斷有新的算法推出,其相關研究領域不僅針對數值模擬計算,對 GIS數據表達、地學分析、計算機視覺、表面目標重構等眾多領域也是一項重要的應用技術[1-3]。網格剖分可分為結構化網格(Structured Grid）和非結構化網格(U nstructured Grid）兩種。非結構化網格易于控制網格單元大小、形狀及網格點位置,可以實現合理分布網格的密度,提高計算精度,因此具有比結構化網格更大的靈活性和對復雜邊界更強的適應性。生成二維非結構化網格的常用方法[4]有四叉樹法(Quadtree）、Delaunay三角剖分法(Delaunay Triangulation,DT）和推進陣面法(Advancing Front M ethod,A FM）以及幾種方法的綜合和改進。由于Delaunay三角網是Vo ronoi圖的直線對偶圖,具有空外接圓和最大最小角特性,還可以盡可能地避免病態三角形的出現,實現約束條件下的三角剖分,因此,在數值計算非結構化網格剖分和GIS三角網(TIN）的數據格式表達中廣泛使用。

目前常見的構建Delaunay三角網的算法有:分治算法(Divide-and-conquer）、逐點插入算法(Incremental Insertion）、生長算法和掃描線算法(Sweepline）[3,5-7],其中以分治算法效率最高,時間復雜度為O(N logN）[5]。由于實際的不同需要,除幾何約束條件外,對三角剖分的角度、大小和節點物理量的控制也都有相應約束,目前的剖分算法尚不能滿足以上全部約束條件進行剖分。本文基于改進的分治算法,提供了梯度變化和針對實際問題的單元和節點重編號,并利用 GIS強大的空間數據處理和分析功能對算法做數據前后處理,保證了高質量約束Delaunay三角網格剖分的生成。

1 約束Delaunay三角化方法

分治算法[5]的基本思路是把輸入點集分割為數個較小的點集,在各個子點集內生成小三角網,再逐級合并的一個遞歸過程。子點集最終只有兩點或者三點,形成一條邊、共線邊或者三角形。算法必須將點集連接為凸區域(凸殼）,對凹區域和多連通區域會產生域外三角形,無法保證約束邊的存在。約束Delaunay三角化的分治算法的關鍵過程及實現方法如下。

1.1 點集劃分

原分治算法只有Y軸方向的分割,這里改進為交替分割[8],即點域X軸方向的長度大于Y軸方向的長度,則以X軸方向對半分割點集,否則以Y軸方向對半分割點集。該過程為遞歸調用的分割過程。

1.2 約束特征線重構

在Delaunay三角剖分上,采用逐次加入約束特征線的方法進行約束重構,可以證明該重構方法是收斂的[4],具體步驟如下:1）對于任一條約束直線,查找與之相交的三角形,即該約束線穿越的所有三角形組成的影響域;2）刪除該影響域內的所有三角形,形成一個多邊形空腔,稱為恢復域;3）以該約束直線為界將恢復域一分為二,形成兩個多邊形子域;4）對這兩個多邊形子域進行Delaunay三角剖分,得到包含約束直線的影響域內的三角剖分;5）重復上述步驟,直至所有的約束線都被加入。

1.3 基于三角形的統一數據結構

Delaunay三角剖分涉及三角形和線段兩類幾何數據,如果將線段表示為退化的三角形勢必增加算法特殊情況判斷的冗余。如圖1a所示,利用“影子”層表示完整的三角形結構,凸殼上的每一線段有一個頂點在無窮遠處(算法中采用空指針）的“影子”三角形,特征線段有兩個“影子”三角形。圖1b中下端虛線三角形為子網合并中新建的“影子”三角形,陰影的三角形采用局部LOP優化算法[7],保證三角形的Delaunay特性。

圖1 分治算法的三角形數據結構和子網合并[9]Fig.1 Data structure of triangle and halfway through a merge step of the dividing and conquering algorithm

1.4 網格細化修正

剖分細化是在網格中插入頂點細化網格,使得網格最大和最小角滿足約束。本文采用Shew chuk的細化算法[9,10],它是對 Ruppert和 Chew提出算法[11,12]的綜合改進算法,頂點插入使用LOP算法[7],保證Delaunay特性,遵循兩個規則:

規則1:一條線段的直徑圓是唯一且最小的包含該線段的圓。如圖2所示,如果有一點位于直徑園內,則稱該線段被侵入。任何被侵入線段在其中點位置插入頂點進行分割。

圖2 線段被遞歸細分直至不被侵入為止[9]Fig.2 Segmentsare split recursively until no segmentsare encroached

規則2:一個三角形的外接圓是唯一經過三角形3個頂點的圓。如果三角形的一個角度過小或面積過大的情況下滿足了約束條件,則稱該三角形為壞三角形。如圖3所示,以外心頂點(外接圓的中心）作為插入點進行分割。

圖3 每個壞三角形線在外心頂點插入分割[9]Fig.3 Each bad triangle is split by inserting a vertex at its circumcenter

剖分細化具體算法如下:1）取約束線段集合中的每一條約束線段,檢查它是否被其他點侵入,是則轉第2步;否則轉第3步。2）按規則1將約束線段分割,將新的約束線段加入約束線段集合,轉第1步。3）取剖分三角形集合的三角形進行檢查,全部合格,轉第5步;否則,對于不滿足約束條件的三角形按規則2插入點。4）如果插入點不侵入約束線段,將該點插入網格中,轉第3步;否則,轉第2步。5）算法收斂停止。

2 GIS支持下約束Delaunay三角網格剖分

高質量Delaunay三角網格意味著必須滿足一定的約束條件,而約束條件是根據工程需要或者數值模擬的條件給定的,如流體動力學問題一般在邊界、結構物和擴散源附近采取網格加密方法,這些都要對輸入點集進行預處理;輸出的剖分三角形內角不能太大也不能太小,面積梯度變化適中。因此,要實現本文提出的非結構網格生成算法,需要提取計算區域邊界點的坐標,包括外邊界和內邊界上各節點光滑、間距的處理、內部約束特征線的設置等。如圖4所示,這些預處理都可以在 GIS的支持下完成,同時可以實現空間數據的管理,保證數據的拓撲準確性和剖分三角網格的正確性。

2.1 邊界點的輸入

剖分算法輸入的為點集的坐標值,非實際的拓撲數據,而且邊界點組織直接影響到剖分輸出的結果,需要在輸入前進行必要的預處理。如圖4所示,按照拓撲關系,將輸入數據劃分為Polygon、Polyline和Point 3種矢量數據,其涵蓋了流體動力學問題中的海岸邊界、島嶼、開邊界和水深等高線等所有輸入數據。其中,海岸邊界采用Douglas-Peucker算法做簡化處理后,根據實際計算需要,對邊界和島嶼邊界線上節點進行光滑處理。

圖4 GIS下Delaunay三角化網格剖分流程Fig.4 Flow diagram illustrating Delaunary triangulation mesh generation in GIS

2.2 Delaunay三角網格剖分

采用本文改進的分治算法實現約束Delaunay三角網格剖分,并將算法集成在A rcGIS平臺中。如圖5所示,直接輸入經過預處理的邊界控制線Polyline、邊界控制區域 Polygon和控制點 Point,可直接輸出Delaunay三角網格剖分網格及節點。

2.3 三角剖分網格拓撲檢查及編輯

如果初始輸入的約束之間沒有形成尖銳角度,算法能夠生成最小角度30°的Delaunay剖分網格,并具有較好的梯度。但在初始約束不當和特殊約束條件下,算法還不能保證滿足所有約束條件,因此對算法生成的剖分三角形的拓撲檢查是必要的。如圖6所示,對局部沒有滿足約束條件的三角網格,可通過GIS下的拓撲編輯修正節點以滿足其約束條件。

2.4 剖分三角節點及網格的重編號

圖5 GIS下的Delaunay三角網格剖分Fig.5 Delaunay triangulation mesh generation in GIS

圖6 三角剖分網格拓撲檢查Fig.6 Topological check on triangulation meshes

算法生成的Delaunay三角網格是沒有規律的,給實際數值模擬帶來不便,如水動數值模型的邊界條件是從開邊界開始設定的,因此需要對剖分網格節點及網格重新編號。這樣不僅可以和數值模型開邊界的編號完全耦合,還使得相鄰網格之間和相鄰節點間的編號跨度減小,從而減少計算中調用各網格信息的時間花費,提高計算效率。重編號算法采用AFM法的思想,以開邊界為起始邊界,按照空間鄰接關系向前推進,直到遍歷所有剖分節點,具體步驟如下[13]:1）如圖7所示,給定起始邊界線(這里為開邊界線）、三角剖分 Polygon和剖分節點 Point數據,設定排序方法;2）根據起始邊界線空間檢索與之相交的節點Point,檢索結果按照設定方法排序;3）根據檢索到的節點Point空間檢索與之相交的三角剖分Polygon,以其中心點為基點對未排序三角剖分Polygon按照設定方法排序;4）根據檢索到的三角剖分Polygon空間檢索與之相交剖分節點Point,判斷是否有未排序節點,如果沒有,算法結束,否則,按照設定方法排序后返回第3步,直到算法收斂。

圖7 三角剖分網格重編號Fig.7 Re-numbering of triangulation meshes

3 算例及應用

采用本文算法分別使用隨機離散點集和實例計算海域進行了測試和應用,測試硬件為 Intel Co re Duo CPU 1.8 GHz,內存2 GB,操作系統 W indow s XP SP3,算法時間不包括輸入/輸出時間。

表1為幾種主要算法在隨機離散點集下的測試結果,從時間復雜度看,改進的分治算法基本與原分治算法保持一致,在使用了“影子”三角形的統一數據結構下,時間上明顯優于其他幾種傳統算法。本文提出的算法在保持原分治算法時間復雜度的前提下,滿足了約束條件以適應實際問題計算的要求,其結果令人滿意。表2為日本博多灣(圖8）和中國渤海區域的海洋數值模擬中的應用算例,Delaunay三角剖分能很好地滿足特征約束,計算速度快,能滿足水動數值模擬非結構網格的計算要求。

表1 隨機離散點集算例Table 1 Examples using random points s

表2 應用算例計算輸入參數及輸出結果Table 2 Input parametersand output results of examples in case study

圖8 日本博多灣Delaunay三角剖分Fig.8 Delaunay triangulation of Hakata Bay,Japan

4 結語

本文在 GIS平臺下基于改進的分治算法實現了高質量約束Delaunay三角網格剖分,改進的分治算法保持了原算法的時間復雜度,并利用 GIS空間數據處理功能完成了數據輸入和輸出剖分網格的準確拓撲檢查,提高了非結構化網格的質量和準確度;同時針對實際問題的需要,基于推進陣面算法的思想,結合空間鄰近拓撲關系實現了三角剖分節點和網格的重新編號算法。算例測試和實例應用表明,本文提出的方法集成了幾種綜合算法的優點,算法收斂速度快,網格生成質量高,滿足實際問題的約束條件,應用效果較好。

[1] 芮一康,王結臣.Delaunay三角形構網的分治掃描線算法[J].測繪學報,2007,36(3）:358-362.

[2] 胡金星,馬照亭,吳煥萍,等.基于網格劃分的海量數據Delaunay三角剖分[J].測繪學報,2004,33(2）:163-167.

[3] 武曉波,王世新,肖春生.Delaunay三角網的生成算法研究[J].測繪學報,1999,28(1）:28-35.

[4] 王盛璽,宋松和,鄒正平.基于約束Delaunay三角化的二維非結構網格生成方法[J].計算物理,2009,26(3）:335-348.

[5] LEED T,SCHACHTER B J.Two algorithms for constructing a Delaunay triangulation[J].Computer and Information Sciences,1980,9(3）:219-242.

[6] FORTUNE S.A sweepline algorithm for VoronoiDiagrams[J].Algorithmica,1987,2(2）:153-174.

[7] LAWSON C L.Software for C1 surface interpolation[A].RICE J R.Mathematical Software III[C].New York:Academic Press,1977.161-194.

[8] DW YER R A.A Faster divide-and-conquer algorithm for constructing Delaunay triangulations[J].Algorithmica,1987,2(2）:137-151.

[9] SHEWCHU K J R.Delaunay refinement algo rithm s fo r triangular mesh generation[J].Computational Geometry,2002,22(1-3）:21-74.

[10] SHEWCHUK JR.Triangle:Engineering a 2D qualitymesh generator and Delaunay triangulator[A].L IN M C,MANOCHA D.Applied Computational Geometry:Towards Geometric Engineering[C].First ACM Workshop on Applied Computational Geometry,Lecture Notes in Computer Science,Vol.1148.Berlin:Springer-Verlag,1996.203-222.

[11] RUPPERT J.A Delaunay refinement algorithm for quality 2-dimensionalmesh generation[J].A lgorithms,1995,18(3）:548-585.

[12] CHEW L P.Guaranteed-quality mesh generation for curved surfaces[A].Proceedings of the Ninth Annual ACM Symposium on Computational Geometry[C].San Diego,CA:ACM Press,1993.274-280.

[13] 趙曉東,王海龍,沈永明.基于 GIS的二維非結構化剖分網格優化[J].地理與地理信息科學,2010,26(2）:46-48.

High-Quality Constrained Delaunay Triangulation Based on GIS

ZHAO Xiao-dong1,YAN Xiao-bao1,SHEN Yong-ming2,WANG Liang2
(1.China-JapanResearchCenterforGeo-environmentalScience,PioneerParkofAcademician,DalianUniversity,Dalian 116622;2.StateKeyLaboratoryofCoastalandOffshore Engineering,DalianUniversityofTechnology,Dalian116023,China）

The unstructured mesh generation isoneof the key technical issues in many fields such asmechanical computation and numerical simulation.Based on the analysis of existing unstructured mesh generation algo rithms,the imp roved divide-and-conquer algo rithm suppo rted by GIS isp roposed to dealwith constrained Delaunay triangulation.Thismethod makes useof the GIS spatial topological relations to handle the p re-p rocessing of input algo rithm data,imp lements Delaunay refinement using a triangle-based data structure,and checks the output meshes with accurate topology and constraints.Based on the idea of advancing f ront algo rithm and spatial topological relations between the mesh nodes and triangulation,the re-numbering algo rithm is also p roposed to facilitate the assignment of theopen boundary conditionsaswell as imp roving the computational efficiency.The app lication to the p ractical simulation show s that the p roposed method can simp lify numerical model p re-p rocessing of unstructured mesh generation,p reserve the advantages of several algorithm s w ith the original running time,and imp rove the quality of constrained Delaunay triangulations.

mesh generation;GIS;constrained Delaunay triangulation

P208

A

1672-0504(2010）05-0024-05

2010-04- 21;

2010-06-10

國家自然科學基金重點項目(50839001）;國家自然科學基金項目(50874021、50779006）;遼寧省高等學?？蒲许椖坑媱?L20100321）

趙曉東(1969-）,男,博士,教授,從事 GIS在采礦、巖土工程應用和水動模型耦合的研究。E-mail:xdong.zhao@gmail.com