摘 要： 在常規逐點插入算法的基礎上，提出了一種改進的逐點插入構建Delaunay三角網的算法。引入散亂點集有序化、三角形單元分類的方法快速生成Delaunay三角網。
關鍵詞： Delaunay三角網 逐點插入 地震層位模型 OpenGL
在三維地理信息系統中，地震層位模型的可視化研究有著廣泛的應用[1]。如在油藏勘探領域，通過重構地表特征層位、判斷其走向可預測儲層的位置。地震層位模型是構建地表層位表面空間位置與其相關屬性信息的數字化表示，它是對地表層位在地震波采樣數據基礎上的表面重構[2]。
由地震波采樣得到并經過速度場解釋之后的地震數據表現為三維空間的散亂數據點集。針對該類數據的特點，普遍采用基于三角網的建模方法構造層位模型。其中，Delaunay三角網具有良好的形態，在表達地質形態方面表現較為出色。
如何快速、高效地構建Delaunay三角網一直是眾多學者研究和關注的焦點。迄今為止出現了不少成熟的算法，基本算法為逐點插入法、分割-合并法及三角網生長法等[3]。其中三角網生長算法由于算法效率低，目前較少采用；分割-合并算法最為高效，但相對復雜，且深度遞歸，對內存要求高；逐點插入法實現簡單、占用內存小，但其時間復雜度較高，效率低于分割-合并算法。
在油藏勘探過程中，有時須補測部分數據，這會引起地震數據的動態變化。因此需要對已存在的Delaunay三角網進行局部更新，按需插入和刪除某些點。這就要求有一種能夠快速構建三角網的方法。
考慮到工程應用的需要，本文在常規逐點插入算法的基礎上，引入散亂點集有序化、三角形單元分類的方法快速構建Delaunay三角網，有效地解決了地震層位建模的工程問題。實測結果表明，改進后算法的時間復雜度大為降低，從而提高了算法的效率。
1 逐點插入算法流程
假設給定散亂數據點集P={Vi|i=1，……，n}，其中Vi的三維坐標為(xi，yi，zi)，n為點的個數。經過Delaunay三角化之后，輸出Delaunay三角網D={Tj|j=1，……，m}，其中Tj表示第j個三角形，m為三角形的個數。逐點插入法的基本思想為：在某一已存在的三角網中插入一點，遍歷所有三角形，查找外接圓包含該點的三角形集合，然后用LOP(Local Optimization Procedure)法則優化[4]，確保所建三角網符合Delaunay法則。以下用偽碼描述逐點插入算法的流程。
SUBROUTINE：Delaunay Triangulate
輸入：散亂數據點集P={Vi|i=1，……，n}
輸出：Delaunay三角網D={Tj|j=1，……，m}
初始化三角形數組
確定包圍三角形，把包圍三角形的頂點添加到頂點數組的末尾
把包圍三角形添加到三角形數組
FOR i=1 to n DO
對于頂點數組中的頂點Vi，初始化邊緩沖區
FOR j=1 to m DO
對于三角形數組中的三角形Tj，計算其外接圓圓心和半徑
IF Vi位于外接圓內 THEN
把Tj的三條邊添加到邊緩沖區，從三角形數組中刪除Tj
ENDIF
ENDFOR
刪除邊緩沖區中所有重復指定的邊，只保留閉合多邊形的邊把由Vi和閉合多邊形的邊形成的所有新三角形添加到三角形數組
ENDFOR
從三角形數組中刪除所有包含包圍三角形頂點的三角形從三角形數組中刪除包圍三角形
END
2 算法改進
2.1 自定義數據結構
合理的數據結構可以提高算法的執行效率。本文設計了點、邊和三角形三種數據類型來存儲空間數據點集，表達三類元素之間的拓撲關系。
(1)點數據類型用來存儲三角網中所有的數據點坐標。
typedef struct
{
float x，y，z； //空間散亂點坐標
} POINT3D；
(2)邊數據類型存儲組成邊的二個頂點的索引號，記錄三角網中的邊與點之間的拓撲關系。
typedef struct
{
long p1，p2； //兩個頂點的索引號
} EDGE；
(3)三角形數據類型存儲組成三角形的三個頂點的索引號，記錄三角網中三角形與點、邊之間的拓撲關系。
typedef struct
{
long p1，p2，p3； //三個頂點的索引號
} TRIANGLE；
可見，空間數據點集坐標存儲在頂點數組中，邊緩沖區和三角形數組存儲對應點的索引號，這樣既節省存儲空間，又便于數據的直接訪問。而且由于數組元素的隨機訪問速度最快，可提高搜索效率。
2.2 散亂點集有序化
常規逐點插入算法依次從存儲散亂數據點集的數組中提取一點，然后遍歷當前三角形數組中所有的三角形，判斷其外接圓是否包含當前待插入點。這一步是點插入過程中最費時的一步，其計算效率取決于三角網的數據結構和搜索方法。本文在對散亂數據點集進行Delaunay三角化之前，首先對其預處理，使散亂點集有序化。具體處理過程如下：(1)遍歷原始散亂數據點集，計算其X軸坐標范圍(Xmin，Xmax)和Y軸坐標范圍(Ymin，Ymax)，通過比較得出坐標分量變化較大的坐標軸；(2)以此坐標軸的坐標分量為比較依據，按照從小到大的順序對散亂數據點集進行排序，得到有序的空間數據點集；(3)將排序之后的點集投影在XOY平面上，刪除X、Y坐標分量相等的重復點，去除奇異性，保證惟一性。
2.3 三角形單元分類
每次在已生成的三角網中插入一點之后，一些不符合Delaunay準則的三角形被刪除，同時又生成一批新的三角形。常規逐點插入算法對這些三角形不加以區分，每次插入操作都必須遍歷所有三角形，計算其外接圓圓心和半徑，判斷插入點是否位于其外接圓內。而實際上，大部分三角形是不需要進行查找、判斷的，這些冗余操作大大降低了算法執行效率，也是導致常規逐點插入算法時間復雜度高的主要原因。本文針對這一瓶頸，對逐點插入過程中生成的三角形單元分類為已完成三角化和未完成三角化二類并加以標識，從而只需要對那些標識為未完成的三角形進行操作，算法效率大為提高。
詳細的優化步驟：(1)每產生一個新三角形，初始化為未完成三角化并賦以標識；(2)在已形成的三角網中插入一點時，循環判斷當前三角形數組中的所有三角形的分類標識。若標識為已完成則跳過該三角形繼續判斷下一個三角形；否則進一步計算其外接圓圓心與半徑，運用LOP法則進行優化，同時根據三角形單元分類規則更新該三角形的分類標識。
三角形單元分類規則是建立在散亂點集有序化基礎上的。假設散亂點集有序化過程中以X軸的坐標分量為比較依據，則定義分類規則如下：
定義 已知插入點Vi(xi，yi，zi)(i=1，2，3，……，n，n為點的個數)和待分類三角形Tj(p1j，p2j，p3j)(j=1，2，3，……，m，m為三角形的個數)，其中Tj的外接圓圓心為Cj(cxj，cyj，czj)、半徑為rj。若xi-cxj≥rj，則三角形Tj標識為已完成三角化，否則標識為未完成三角化。
分類規則的直觀解釋如圖1所示。若Tj外接圓圓心Cj到經過插入點Vi且垂直于X軸的直線MN的距離大于其半徑rj的長度，則可保證所有X坐標分量不小于插入點Vi的點都位于Tj外接圓外。又因為散亂點集有序化之后按照X坐標分量以升序排列，故對于點Vi之后的插入點Vk(k> i)而言，永遠都不會位于三角形Tj的外接圓內，從而Tj可標識為已完成三角化。
3 算法比較與分析
算法用ANSI C++語言進行了編程實現，微機環境為Pentium(R)4，CPU主頻2.4GHz，內存512MB，實驗數據為程序隨機生成的單精度型點集，實驗結果如表1所示。
由實驗結果可以看出，改進后的算法比原算法在速度上有數倍增長，特別是當點數增加時速度增長的幅度也增大，而占用的內存空間只有少量增長。這充分體現了算法的穩定性。在實際運行時對時間和空間的合理綜合運用提高了算法的效率和性能。
4 工程應用
地震層位模型的可視化已成為油藏勘探領域研究的熱點。將Delaunay三角網建模方法與OpenGL圖形可視化技術相結合，即可快速構建地震層位模型，實現地表特征層表面三維顯示，直觀地表達特征層的空間分布與形態。
OpenGL是當今工業界最先進最開放的圖形API，它提供了強大的繪制二維和三維圖形的能力。如提供了大量的圖形變換函數，可以方便地將三維圖形顯示在屏幕窗口進行放大、縮小、旋轉等操作；OpenGL還提供了線面消隱、著色關照、紋理映射和反走樣等技術的一系列函數，可以方便地對地層表面進行處理，增強圖形的真實感。
本文基于上述方法，以地震波采樣解釋之后的地震數據為數據源，重構了勝利油田某地區的層位模型如圖2、圖3所示。模型顯示結果驗證了本文的算法是行之有效的。
5 結 論
本文在仔細研究常規逐點插入算法的基礎上，對其進行了優化，使得改進之后的算法不但繼承了原算法流程清晰、實現簡單、占用內存較小、可局部動態更新三角網的優點，而且提高了執行效率。該算法運用于實際工程中構建地震層位模型，取得了良好的效果。
當前，包括地震層位建模在內的三維地理信息系統正處于從研究走向實用的過渡階段。在實際應用中，有時還存在著約束三角剖分，也就是部分散亂點之間常常存在著某種約束關系，如地表模型中的斷裂線等。在對此種數據集進行三角網生成時，生成的三角網應保持其原有的約束關系，即約束數據集下的三角剖分。因此，如何在約束條件下快速準確地生成Delaunay三角網，是下一步研究工作的重點。
參考文獻
1 包世泰，夏斌，崔學軍等.地質信息模型研究及其應用.大地構造與成礦學，2004；28(4)
2 胡金星，潘懋，馬照亭等.高效構建Delaunay三角網數字地形模型算法研究.北京大學學報(自學科學版)，2003；39(5)
3 武曉波，王世新，肖春生.Delaunay三角網的生成算法研究. 測繪學報，1999；28(1)
4 Lee D T，Schachter B J.Two Algorithms for Constructing a Delaunay Triangulation.Int J.of Computer and Information Sciences，1980；9(3)