孟憲海，成文迪，徐 博，楊 欽

（1.北京航空航天大學(xué)計算機學(xué)院，北京 100191；2.中國石油天然氣勘探開發(fā)公司，北京 100034）

1 緒 論

網(wǎng)格生成，即是將復(fù)雜幾何結(jié)構(gòu)離散成一組簡單幾何單元的過程，生成三角形或四面體網(wǎng)格的過程又稱三角化。其中，Delaunay三角形/四面體網(wǎng)格具有很高的質(zhì)量，對復(fù)雜區(qū)域有很好的逼近性，應(yīng)用最為廣泛。

Delaunay三角網(wǎng)格在1934年由俄國數(shù)學(xué)家Delaunay提出[1]，在此類網(wǎng)格中，任何一個三角形/四面體網(wǎng)格單元的外接圓/球中，均不包含網(wǎng)格中的其他頂點。Delaunay三角網(wǎng)格以其優(yōu)良的幾何性質(zhì)，在工程中得到了十分廣泛的應(yīng)用。

點集的Voronoi圖，由每個點的Voronoi單元組成，點a的Voronoi單元，定義為空間中所有到a點的歐氏距離不大于到點集中其他任何點的歐氏距離的所有點的集合。點集的Delaunay三角化結(jié)果，和點集的Voronoi圖互相對偶。即如果點集中任意兩點的Voronoi單元之間存在公共面，則這兩點的連線即為Delaunay三角化中的一條邊。二者是彼此等價的。

傳統(tǒng)的點集Delaunay三角形/四面體網(wǎng)格化算法，可以生成高質(zhì)量的網(wǎng)格，還能處理一定的限定條件[2-3]。但近年來，由于高精度科學(xué)計算，數(shù)值仿真等領(lǐng)域的需求，需要處理更大數(shù)據(jù)量、生成更高精度的網(wǎng)格，還對網(wǎng)格的生成次序、形態(tài)等有特殊的要求，這些因素促使人們不斷尋求新的網(wǎng)格生成思路。

經(jīng)典的Delaunay三角化生成算法已經(jīng)較為成熟。其中以B/W增量算法應(yīng)用最為廣泛[4-5]。在用B/W算法進行點集的Delaunay三角化時，首先用少量點生成初始的粗網(wǎng)格；之后，每插入一個新點時，檢查哪些三角網(wǎng)格單元的外接球包含這個點，把他們刪去，形成一個空洞；最后將新點與組成空洞的每個頂點相連，形成新的網(wǎng)格。不斷重復(fù)這個過程直到所有點都被加入到網(wǎng)格之中，算法終止[6]。

在這種思路下，每插入一個點，都要在現(xiàn)有的整體網(wǎng)格中搜索外接球包含它的網(wǎng)格單元，每一個新網(wǎng)格的生成都與當(dāng)前整體網(wǎng)格相關(guān)，算法具有全局性、連續(xù)性。雖然效率高，但當(dāng)存在某些限定條件、或者無法得到全局信息時，該算法不太適用。但實際上，Delaunay三角網(wǎng)格是具有局部性特點的，本文即是充分利用這個性質(zhì)，提出一種新的Delaunay三角化策略。

2 Delaunay三角網(wǎng)格的局部性原理

2.1 Delaunay三角網(wǎng)格的局部性特點

傳統(tǒng)的B/W增量算法[4-5]是基于整個點集的。但實際上，一個點的三角網(wǎng)格，只與其鄰近有限距離內(nèi)的點相關(guān)。我們可以充分利用這個性質(zhì)，來設(shè)計新的生成策略。

2.2 最小Voronoi鄰近點集

點集的Delaunay三角網(wǎng)格中，與某個點相關(guān)的網(wǎng)格單元，只由鄰近局部范圍內(nèi)的有限點決定。我們把這些點，稱為該生長點的最小Voronoi鄰近點集。

定義1 設(shè)S為空間點集，點o∈S，點集M為S的子集。任取p∈S，若p的Voronoi單元和o的Voronoi單元相鄰（Vor(p) ∩Vor(o)≠空集），則p∈M；若p的Voronoi單元和o的Voronoi單元不相鄰（Vor(p) ∩Vor(o)＝空集），則p?M。稱M為生長點o在點集S中的最小Voronoi鄰近點集（Minimum Voronoi Neighbors，簡稱MVN）。

對于點集中的每個點，它作為生長點時的Voronoi單元的形狀只與其最小Voronoi鄰近點集相關(guān)，而且最小Voronoi 鄰近點集中點的個數(shù)必然是有限個，分布在該點的局部有限距離內(nèi)。因此，只要找到其MVN點集，就可以在當(dāng)前生長點的局部構(gòu)建Voronoi單元了。只要構(gòu)造出整個點集的Voronoi圖，就可以通過對偶得到點集的最終Delaunay三角網(wǎng)格。

為了從點集中找到生長點o的MVN，我們給出如下準則：

設(shè)S為空間點集，點o∈S，M為S的一個子集，T為M中的點與o點組成的Delaunay三角形/四面體網(wǎng)格單元集合。如果M滿足：

1）T中的每個三角形/四面體網(wǎng)格單元都以o點為一個頂點。

2）T中的三角形/四面體網(wǎng)格單元組成一個凸包，如果o點不在凸殼上，則T圍繞o點構(gòu)成星形三角形/四面體集合。如果o在凸殼上，則凸殼邊界包含以o點為頂點的網(wǎng)格單元。

3）T中任意網(wǎng)格單元的外接圓/球中均不包含S中的其他點，則M為o點的MVN。

上面給出的三條規(guī)則，可以作為o點的最小Voronoi鄰近點集判定準則。

3 Delaunay局部網(wǎng)格生成算法

3.1 算法流程

3.1.1 最小Voronoi鄰近點集搜索策略

為了尋找點集中某個點的最小Voronoi鄰近點集，需要對點集中的點進行遍歷，并依照空球準則對每個點進行處理。處理過程實際上是一個局部加點過程。

在B/W增量算法中，用于初始化的粗網(wǎng)格是覆蓋全體點集的，因此，每一次被加入的點，一定落在現(xiàn)有的某個三角形/四面體網(wǎng)格單元內(nèi)。但是在本文提出的局部遍歷過程中，初始網(wǎng)格并不是覆蓋全體的，在現(xiàn)有網(wǎng)格單元達到閉合之前，我們需要額外維護一個結(jié)構(gòu)來記錄現(xiàn)有網(wǎng)格的邊界。

二維情況下，需要維護的邊界記錄比較簡單。如圖1(1)所示，在加入初始的兩個點之后，我們得到一個Delaunay三角網(wǎng)格單元，當(dāng)前網(wǎng)格的邊界就是由生長點和兩個加入點組成的兩條直線，平面被分成4個區(qū)域?？梢院唵蔚陌阉鼈冇涗洖椤白簏c”和“右點”。

在局部網(wǎng)格對生長點達到閉合之前，每次加入新點，需要首先判斷它與邊界記錄的關(guān)系，即判斷它落在如圖所示的A、B、C、D四個區(qū)域中的哪一個里：

1）落在A區(qū)域：判斷干涉情況加點，不更新左/右點，如圖1(2)。

2）落在B/C區(qū)域：判斷干涉情況加點，同時更新左/右點，如圖1(3)。

3）落在D區(qū)域：判斷干涉情況加點，達到閉合，如圖1(4)。

達到閉合之后，舍棄左點和右點的記錄。

圖1 二維條件下的局部加點過程示意

三維情況要比二維要復(fù)雜些。首先，三維情況下的邊界記錄是面，空間被劃分成若干區(qū)域。其次，邊界記錄的數(shù)目在加點過程中可能變化。

圖2展示了三維情況下維護邊界記錄的過程。在圖2(1)中，生長點和初始加入的3個點構(gòu)成了第一個Delaunay四面體網(wǎng)格單元。此時，邊界由它的3個包含生長點的面確定。在加入第四個點時，需要首先判斷它落在八個區(qū)域中的哪一個，從而采取不同的策略。

隨著加點過程的繼續(xù)，邊界記錄中的記錄數(shù)可能發(fā)生變化，考慮圖2(2)，新點D對于當(dāng)前邊界記錄(OAB,OBC,OCA)的位置向量是(-,+,+)，D點的加入會導(dǎo)致OAB所在邊界面被兩個新的邊界面OAD和ODB所取代，邊界記錄的長度變成4。上述過程一直繼續(xù)，直到出現(xiàn)圖2(5)，這時新加入的點P位于所有邊界面的負側(cè)，加入P點之后，當(dāng)前局部Delaunay四面體網(wǎng)格對生長點o達到了閉合，如圖2(6)。此后的加點過程可以用統(tǒng)一的方式來處理。

圖2 三維情況下的局部加點過程示意

在當(dāng)前局部Delaunay四面體網(wǎng)格對生長點o閉合前，維護一個邊界面的列表。在每次加入一個新點時，首先計算它和邊界列表中每個面的位置關(guān)系，得到它的位置向量。對于位置向量中的每個元素：

1）如果為"-"，處理下一個元素。

2）如果為"+"，刪除該元素對應(yīng)的位置向量，并形成兩個新的候選邊界面?？疾煸撛卦谖恢孟蛄恐械那耙粋€元素，如果為"+"，刪除前一個候選面；同時考察該元素在位置向量中的前一個元素，如果為"-"，刪除后一個候選面。

三維情況下判斷干涉加點的過程也與二維情況類似，可以把二維情況下的局部加點算法當(dāng)成三維情況下邊界記錄鏈表長度固定的一個特例。

給定一個點集，我們可以對其中的每個點，利用上述搜索算法，生成其局部Delaunay網(wǎng)格。這個算法稱為Delaunay局部網(wǎng)格生成算法。

3.1.2 引入背景索引網(wǎng)格的輔助

如果處理每個點時，都要遍歷整個點集，效率將會十分低下。但Delaunay網(wǎng)格的局部性特點，保證了每一次遍歷都會在生長點局部的有限范圍內(nèi)中止。當(dāng)兩個點同時插入一個Delaunay三角網(wǎng)格時，如果它們距離較遠，彼此一定互不影響，則稱它們?yōu)榉菦_突點。反之，稱它們?yōu)闆_突點。有如下定理[7-8]：

定理1 如果r是網(wǎng)格中所有單元外接圓半徑的上界，則距離大于4r的兩個點一定是非沖突點。

對于給定的點集，r的距離是一定的。因此，Voronoi鄰近點集一定分布在生長點周圍的有限區(qū)域內(nèi)。但r的具體值很難預(yù)先得到，我們可以采取間接的手段，建立一個二級索引數(shù)據(jù)結(jié)構(gòu)，人為將點集按照位置關(guān)系劃分為許多單元塊，稱為背景索引網(wǎng)格。背景索引網(wǎng)格應(yīng)該盡可能簡單、規(guī)則。這里選擇最簡單的矩形/立方體。

首先將點集中的每個點布置到其所對應(yīng)的三角網(wǎng)格單元中。之后，對當(dāng)前生長點o所在的背景網(wǎng)格單元及其一階鄰接單元進行檢索，而后再逐級擴展至更高階鄰接單元上的網(wǎng)格單元。當(dāng)某個網(wǎng)格單元與當(dāng)前生長點的局部Delaunay三角形/四面體集合中所有三角形/四面體的外接圓/球的交為空，則不再對該網(wǎng)格單元的鄰接單元進行檢索。由于最小Voronoi鄰近點集中的點是有限個，遍歷過程會在有限步驟內(nèi)終止，如圖3所示。

采用鏈表式的方式管理需要遍歷的背景索引網(wǎng)格單元，可以保證遍歷順序是以節(jié)點為中心環(huán)繞進行，逐級向外的，能使得當(dāng)前局部網(wǎng)格盡早達到閉合。

圖3 背景索引網(wǎng)格輔助下的一個點的MVN

3.2 算法描述

Delaunay局部網(wǎng)格生成算法分為兩個步驟，一是背景輔助網(wǎng)格的生成，二是針對每個點的局部網(wǎng)格生成。生成背景輔助網(wǎng)格的過程是串行的，對輸入點的分布進行分析，創(chuàng)建背景索引網(wǎng)格鏈表，把點集中的點布置到對應(yīng)的背景網(wǎng)格單元中，并建立背景網(wǎng)格單元的鄰接關(guān)系。這里不再贅述。

針對每個點的每個局部網(wǎng)格生成，以三維情況為例，描述如下：

算法3.2每個點的最小鄰近點集生成算法

SearchMVNwithBackRect(growPoint,growBox,backGridSet)；

算法輸入：矩形背景網(wǎng)格集合backGridSet，生長點growPoint，生長點所在的背景網(wǎng)格單元growBox；

算法輸出：生長點growPoint的最小鄰近點集；

1）定義backRectList為需要遍歷的背景網(wǎng)格單元集合隊列，并定義整型動態(tài)數(shù)組backRectFlags，存放已經(jīng)加入遍歷隊列的背景網(wǎng)格單元在背景網(wǎng)格集合backGridSet中的下標。首先將growBox加入到backRectList，并把growRect的下標加入backRectFlags；

2）設(shè)i=0，對backRectList中的第i個單元t，對t中的每個點輪流執(zhí)行加點算法addPointToTetraList。通過返回值表明該點是否被加入當(dāng)前局部網(wǎng)格；

3）如果t中任意一點被加入當(dāng)前局部網(wǎng)格，則逐一考察t的鄰接單元。在立方體背景網(wǎng)格的條件下，通過計算下標的方式即可得到t的鄰接單元；

4）查找每個鄰接單元的下標是否存在于backRectFlags中，如果存在，說明該單元已經(jīng)在遍歷隊列中。否則，把該單元加入backRectList，下標加入backRectFlags；

5）如果t中沒有點被加入當(dāng)前局部網(wǎng)格，則判斷當(dāng)前局部網(wǎng)格中所有三角形的外接圓是否與t相交，如果存在任意外接圓與t相交，則再次執(zhí)行步驟3)、4)，逐一考察t的鄰接單元，查找其下標是否存在于backRectFlags中，如果存在，說明該單元已經(jīng)在遍歷隊列中。否則，把該單元加入backRectList，下標加入backRectFlags;

6）i++，如果i小于backRectList的當(dāng)前長度，則處理backRectList中的下一個t。

其中UpdateTriToDTSet為局部Delaunay三角形集合三角形更新算法，UpdatePtToDTSet為局部Delaunay三角形集合點更新算法，均指按照Delaunay準則將三角形newT和點p加入到核心點o的鄰近Delaunay三角形集合中，不再贅述。

最終結(jié)果網(wǎng)格中的每個單元，都在以不同頂點為生長點時重復(fù)生成過。在實際計算中，我們可以給每個點設(shè)置一個“是否進行了局部網(wǎng)格生成”的標志變量。避免重復(fù)加入。

3.3 算法執(zhí)行效率分析

假設(shè)輸入點的個數(shù)為n，背景索引網(wǎng)格的生成過程較為簡單，在對輸入點集進行一次遍歷后即可完成，時間復(fù)雜度為O(n)。

每個點的最小Voronoi鄰近點集生成算法的過程則較為復(fù)雜。在算法中，每個點的最小Voronoi鄰近點集數(shù)目有限，且遍歷過程都是集中在點集附近的，所以對每個點來說，生成其局部網(wǎng)格的時間復(fù)雜度介于O(n)到O(1)之間，我們可以根據(jù)網(wǎng)格中的數(shù)量關(guān)系進行簡單估計。

先考慮二維情況下，根據(jù)歐拉定理，對于三角網(wǎng)格中的數(shù)量關(guān)系，可以得到下述結(jié)論[7]：

定理2 平面三角網(wǎng)格中，三角網(wǎng)格單元數(shù)T、邊數(shù)E和點數(shù)n是同一數(shù)量級，即O(T)=O(E)=O(n)。

對于一個點來說，以它為生長點的三角網(wǎng)格單元數(shù)是O(1)級的，根據(jù)經(jīng)驗，在點集分布較為均勻時，所需遍歷到的點的數(shù)目一般是最小Voronoi鄰近點集大小的5倍左右，由于背景索引網(wǎng)格的存在，這個基數(shù)不會太大。我們可以把二維Delaunay三角網(wǎng)格局部生成算法的時間復(fù)雜度歸納為O(kn)，其中k是個較大的常數(shù)。

在三維情況下，每個點最小Voronoi鄰近點集的基數(shù)較大。從歐拉定理可知[7]：

定理3 三維四面體網(wǎng)格的邊數(shù)E、面數(shù)F和四面體數(shù)T的下限是O(n)，上限是O(n2)。

在三維情況下，對于一個點來說，以它為生長點的四面體網(wǎng)格單元數(shù)是下限是O(1)，上限是O(n)。整個算法的下限是O(n)，上限是O(n2)。

不過在大多數(shù)實際問題中，三維四面體網(wǎng)格的數(shù)目并不會達到上限，該算法在一般情況下仍然具有可行性和實用價值。

4 算法的應(yīng)用

4.1 地質(zhì)斷裂層面上的PEBI網(wǎng)格生成

在油藏數(shù)值模擬計算中，常需要精確描述地質(zhì)層面的網(wǎng)格，PEBI網(wǎng)格是較好的選擇[9]。PEBI網(wǎng)格在定義上等價于Voronoi圖，生成PEBI網(wǎng)格的問題也就是確定生長點集并生成Voronoi圖的問題。

在實際問題中，所研究的地質(zhì)層面經(jīng)常具有復(fù)雜的地質(zhì)斷裂構(gòu)造，傳統(tǒng)的二維限定PEBI網(wǎng)格生成算法，是從整體生成點集的Delaunay圖，然后通過對偶的方式獲得Voronoi圖。這種方法對限定條件的處理復(fù)雜、繁瑣，通用性差。特別是當(dāng)其應(yīng)用于逆斷層附近時，由于其該局部在垂直方向上的重疊，從而會造成上下盤層面上的點之間相互干擾，生成錯誤的Voronoi單元。

圖4 地質(zhì)層面模型上的正斷層與逆斷層

在這種情況下，本文提出的基于索引的局部Delaunay三角化算法可以派上用場。利用搜索最小Voronoi鄰近點集的方法在生長點局部構(gòu)造Voronoi單元可以有效的避免全局搜索所帶來的逆斷層上下盤點之間相互干擾的問題，只需合理地設(shè)置邊界附近的生長點，對斷裂邊界的恢復(fù)問題也能夠得到很好的解決。圖5和圖6顯示了利用Delaunay對包含逆斷層的地質(zhì)層面進行PEBI網(wǎng)格生成后的結(jié)果。

圖5 利用本文算法生成的PEBI網(wǎng)格

圖6 圖5中逆斷層附近的PEBI網(wǎng)格

4.2 Delaunay三角化算法的并行化處理

子步驟之間的獨立性，是該算法的另一個優(yōu)勢，使得它十分適合進行共享內(nèi)存、無需通信的并行化處理。適合在同一臺機器上進行多核并行，也可以在支持分布式共享內(nèi)存的機群上運行。

表1顯示了三維條件下算法在不同線程數(shù)下并行執(zhí)行所耗的時間，表格中還包含了本算法與傳統(tǒng)的B/W增量算法的執(zhí)行效率的比較。從結(jié)果可以看出，隨著點集的規(guī)模逐漸增大，和并行度的增加，并行算法的優(yōu)勢逐漸顯露，其時間復(fù)雜度優(yōu)于傳統(tǒng)的未經(jīng)優(yōu)化的B/W增量算法。

表1 三維條件下算法并行執(zhí)行時間

4.3 無線傳感器網(wǎng)絡(luò)中的定位問題

無線傳感器網(wǎng)絡(luò)（wireless sensor networks,WAN）是一個由大量簡單、廉價的傳感器集成無線通信接口所組成的網(wǎng)絡(luò)，廣泛應(yīng)用于環(huán)境監(jiān)測、災(zāi)難救助等領(lǐng)域。在無線傳感器網(wǎng)絡(luò)的應(yīng)用中，傳感節(jié)點自我定位是一個基礎(chǔ)問題[10]。

目前的定位算法通常是依靠傳感器收集通信或者跳數(shù)信息，來計算節(jié)點位置的。但是，在一個大型的無線傳感器網(wǎng)絡(luò)中，收集全局的信息是非常困難的，所以采用分布式的定位算法，而且消息的復(fù)雜度不能過高。上文提出的基于最小Voronoi鄰近點集的Voronoi圖生成方法，用于節(jié)點定位是很合適的。在定位某個節(jié)點S時，通過廣播消息收集周圍節(jié)點的拓撲信息，生成自己的局部Voronoi鄰近點集，每個節(jié)點不需要知道全局的節(jié)點分布信息也可以正確定位。

5 結(jié)論與展望

本文以Delaunay三角網(wǎng)格的局部性原理為理論基礎(chǔ)，研究并實現(xiàn)一種基于Voronoi最小鄰近點集的Delaunay局部網(wǎng)格生成算法，它是基于節(jié)點的，子過程具有局部性和獨立性，因此可以用于解決一些傳統(tǒng)算法不適用的問題，包括地質(zhì)斷裂層面上的PEBI網(wǎng)格生成，Delaunay網(wǎng)格并行化生成，無線傳感器網(wǎng)絡(luò)定位等。

實現(xiàn)點集的網(wǎng)格生成對于整個研究過程還只是一個開端，后續(xù)工作主要在以下兩個方面進行。一是改善算法本身，使其能夠應(yīng)用于更多的實際問題。二是進一步研究網(wǎng)格細化、質(zhì)量改善、簡單限定條件的恢復(fù)、復(fù)雜限定條件恢復(fù)等問題。

[1]Delaunay B N.Sur la sphere vide.Bulletin [M].Acade′mie des Sciences URSS.1934: 793-800.

[2]Hazlewood C.Approximating constrained tetrahedrizations [J].Computer Aided Geometric Design, 1993, (1): 67-87.

[3]Weatherill N P, Hassan O.Efficient three dimensional Delannay triangulation with automatic point creation and imposed boundary constraints [J].International Journal for Numerical Methods in Engineering, 1994,37(2): 2005-2039.

[4]Bowyer A.Computing Dirichlet tessellations [J].The Computer Journal, 1981, 24(2): 162-166.

[5]Watson D F.Computing the Delaunay tessellation with applications to Voronoi Polytopes [J].The Computer Journal, 1981, 24(2): 167-172.

[6]Thompson J F.Handbook of grid generation [M].CRC Press LLC.1999: 644-661.

[7]閔衛(wèi)東, 唐 澤.圣三角網(wǎng)格中的數(shù)量關(guān)系[J].計算機輔助設(shè)計與圖形學(xué)學(xué)報, 1996, 8(2): 81-86.

[8]Toomre A.On the distribution of matter within highly flattened galaxies [J].The Astrophysical Journal, 1963,138: 385-392.

[9]孫 邁.面向油藏數(shù)值模擬的PEBI網(wǎng)格生成技術(shù)研究[D].北京: 北京航空航天大學(xué), 2007.

[10]王繼春, 劉黃生.基于Voronoi圖的無需測距的無線傳感器網(wǎng)絡(luò)節(jié)點定位算法[J].計算機研究與發(fā)展,2008, 45(1): 119-125.