摘   要： 在常規(guī)逐點(diǎn)插入算法的基礎(chǔ)上，提出了一種改進(jìn)的逐點(diǎn)插入構(gòu)建Delaunay三角網(wǎng)的算法。引入散亂點(diǎn)集有序化、三角形單元分類的方法快速生成Delaunay三角網(wǎng)。
關(guān)鍵詞： Delaunay三角網(wǎng)  逐點(diǎn)插入  地震層位模型  OpenGL

　　在三維地理信息系統(tǒng)中，地震層位模型的可視化研究有著廣泛的應(yīng)用[1]。如在油藏勘探領(lǐng)域，通過重構(gòu)地表特征層位、判斷其走向可預(yù)測儲層的位置。地震層位模型是構(gòu)建地表層位表面空間位置與其相關(guān)屬性信息的數(shù)字化表示，它是對地表層位在地震波采樣數(shù)據(jù)基礎(chǔ)上的表面重構(gòu)[2]。
　　由地震波采樣得到并經(jīng)過速度場解釋之后的地震數(shù)據(jù)表現(xiàn)為三維空間的散亂數(shù)據(jù)點(diǎn)集。針對該類數(shù)據(jù)的特點(diǎn)，普遍采用基于三角網(wǎng)的建模方法構(gòu)造層位模型。其中，Delaunay三角網(wǎng)具有良好的形態(tài)，在表達(dá)地質(zhì)形態(tài)方面表現(xiàn)較為出色。
　　如何快速、高效地構(gòu)建Delaunay三角網(wǎng)一直是眾多學(xué)者研究和關(guān)注的焦點(diǎn)。迄今為止出現(xiàn)了不少成熟的算法，基本算法為逐點(diǎn)插入法、分割-合并法及三角網(wǎng)生長法等[3]。其中三角網(wǎng)生長算法由于算法效率低，目前較少采用；分割-合并算法最為高效，但相對復(fù)雜，且深度遞歸，對內(nèi)存要求高；逐點(diǎn)插入法實(shí)現(xiàn)簡單、占用內(nèi)存小，但其時(shí)間復(fù)雜度較高，效率低于分割-合并算法。
　　在油藏勘探過程中，有時(shí)須補(bǔ)測部分?jǐn)?shù)據(jù)，這會引起地震數(shù)據(jù)的動態(tài)變化。因此需要對已存在的Delaunay三角網(wǎng)進(jìn)行局部更新，按需插入和刪除某些點(diǎn)。這就要求有一種能夠快速構(gòu)建三角網(wǎng)的方法。
考慮到工程應(yīng)用的需要，本文在常規(guī)逐點(diǎn)插入算法的基礎(chǔ)上，引入散亂點(diǎn)集有序化、三角形單元分類的方法快速構(gòu)建Delaunay三角網(wǎng)，有效地解決了地震層位建模的工程問題。實(shí)測結(jié)果表明，改進(jìn)后算法的時(shí)間復(fù)雜度大為降低，從而提高了算法的效率。
1  逐點(diǎn)插入算法流程
　　假設(shè)給定散亂數(shù)據(jù)點(diǎn)集P={Vi|i=1，……，n}，其中Vi的三維坐標(biāo)為(xi，yi，zi)，n為點(diǎn)的個(gè)數(shù)。經(jīng)過Delaunay三角化之后，輸出Delaunay三角網(wǎng)D={Tj|j=1，……，m}，其中Tj表示第j個(gè)三角形，m為三角形的個(gè)數(shù)。逐點(diǎn)插入法的基本思想為：在某一已存在的三角網(wǎng)中插入一點(diǎn)，遍歷所有三角形，查找外接圓包含該點(diǎn)的三角形集合，然后用LOP(Local Optimization Procedure)法則優(yōu)化[4]，確保所建三角網(wǎng)符合Delaunay法則。以下用偽碼描述逐點(diǎn)插入算法的流程。
　　SUBROUTINE：Delaunay Triangulate
　　輸入：散亂數(shù)據(jù)點(diǎn)集P={Vi|i=1，……，n}
　　輸出：Delaunay三角網(wǎng)D={Tj|j=1，……，m}
　　初始化三角形數(shù)組
　　確定包圍三角形，把包圍三角形的頂點(diǎn)添加到頂點(diǎn)數(shù)組的末尾
　　把包圍三角形添加到三角形數(shù)組
　　FOR i=1 to n DO
      　　對于頂點(diǎn)數(shù)組中的頂點(diǎn)Vi，初始化邊緩沖區(qū)
      　　FOR j=1 to m DO
            　　對于三角形數(shù)組中的三角形Tj，計(jì)算其外接圓圓心和半徑
　　            IF Vi位于外接圓內(nèi) THEN
        　　　      把Tj的三條邊添加到邊緩沖區(qū)，從三角形數(shù)組中刪除Tj
　　            ENDIF
      　　ENDFOR
　　      刪除邊緩沖區(qū)中所有重復(fù)指定的邊，只保留閉合多邊形的邊把由Vi和閉合多邊形的邊形成的所有新三角形添加到三角形數(shù)組
　　ENDFOR
　　從三角形數(shù)組中刪除所有包含包圍三角形頂點(diǎn)的三角形從三角形數(shù)組中刪除包圍三角形
　　END
2  算法改進(jìn)
2.1 自定義數(shù)據(jù)結(jié)構(gòu)
　　合理的數(shù)據(jù)結(jié)構(gòu)可以提高算法的執(zhí)行效率。本文設(shè)計(jì)了點(diǎn)、邊和三角形三種數(shù)據(jù)類型來存儲空間數(shù)據(jù)點(diǎn)集，表達(dá)三類元素之間的拓?fù)潢P(guān)系。
　　(1)點(diǎn)數(shù)據(jù)類型用來存儲三角網(wǎng)中所有的數(shù)據(jù)點(diǎn)坐標(biāo)。
　　typedef struct
　　　{
        　　float x，y，z；        //空間散亂點(diǎn)坐標(biāo)
　　} POINT3D；
　　(2)邊數(shù)據(jù)類型存儲組成邊的二個(gè)頂點(diǎn)的索引號，記錄三角網(wǎng)中的邊與點(diǎn)之間的拓?fù)潢P(guān)系。
　　typedef struct
　　{
        　　long p1，p2；        //兩個(gè)頂點(diǎn)的索引號
　　} EDGE；
　　(3)三角形數(shù)據(jù)類型存儲組成三角形的三個(gè)頂點(diǎn)的索引號，記錄三角網(wǎng)中三角形與點(diǎn)、邊之間的拓?fù)潢P(guān)系。
　　typedef struct
　　{
　　        long p1，p2，p3；    //三個(gè)頂點(diǎn)的索引號
　　} TRIANGLE；
　　可見，空間數(shù)據(jù)點(diǎn)集坐標(biāo)存儲在頂點(diǎn)數(shù)組中，邊緩沖區(qū)和三角形數(shù)組存儲對應(yīng)點(diǎn)的索引號，這樣既節(jié)省存儲空間，又便于數(shù)據(jù)的直接訪問。而且由于數(shù)組元素的隨機(jī)訪問速度最快，可提高搜索效率。
2.2 散亂點(diǎn)集有序化
　　常規(guī)逐點(diǎn)插入算法依次從存儲散亂數(shù)據(jù)點(diǎn)集的數(shù)組中提取一點(diǎn)，然后遍歷當(dāng)前三角形數(shù)組中所有的三角形，判斷其外接圓是否包含當(dāng)前待插入點(diǎn)。這一步是點(diǎn)插入過程中最費(fèi)時(shí)的一步，其計(jì)算效率取決于三角網(wǎng)的數(shù)據(jù)結(jié)構(gòu)和搜索方法。本文在對散亂數(shù)據(jù)點(diǎn)集進(jìn)行Delaunay三角化之前，首先對其預(yù)處理，使散亂點(diǎn)集有序化。具體處理過程如下：(1)遍歷原始散亂數(shù)據(jù)點(diǎn)集，計(jì)算其X軸坐標(biāo)范圍(X_min，X_max)和Y軸坐標(biāo)范圍(Y_min，Y_max)，通過比較得出坐標(biāo)分量變化較大的坐標(biāo)軸；(2)以此坐標(biāo)軸的坐標(biāo)分量為比較依據(jù)，按照從小到大的順序?qū)ιy數(shù)據(jù)點(diǎn)集進(jìn)行排序，得到有序的空間數(shù)據(jù)點(diǎn)集；(3)將排序之后的點(diǎn)集投影在XOY平面上，刪除X、Y坐標(biāo)分量相等的重復(fù)點(diǎn)，去除奇異性，保證惟一性。
2.3 三角形單元分類
　　每次在已生成的三角網(wǎng)中插入一點(diǎn)之后，一些不符合Delaunay準(zhǔn)則的三角形被刪除，同時(shí)又生成一批新的三角形。常規(guī)逐點(diǎn)插入算法對這些三角形不加以區(qū)分，每次插入操作都必須遍歷所有三角形，計(jì)算其外接圓圓心和半徑，判斷插入點(diǎn)是否位于其外接圓內(nèi)。而實(shí)際上，大部分三角形是不需要進(jìn)行查找、判斷的，這些冗余操作大大降低了算法執(zhí)行效率，也是導(dǎo)致常規(guī)逐點(diǎn)插入算法時(shí)間復(fù)雜度高的主要原因。本文針對這一瓶頸，對逐點(diǎn)插入過程中生成的三角形單元分類為已完成三角化和未完成三角化二類并加以標(biāo)識，從而只需要對那些標(biāo)識為未完成的三角形進(jìn)行操作，算法效率大為提高。
　　詳細(xì)的優(yōu)化步驟：(1)每產(chǎn)生一個(gè)新三角形，初始化為未完成三角化并賦以標(biāo)識；(2)在已形成的三角網(wǎng)中插入一點(diǎn)時(shí)，循環(huán)判斷當(dāng)前三角形數(shù)組中的所有三角形的分類標(biāo)識。若標(biāo)識為已完成則跳過該三角形繼續(xù)判斷下一個(gè)三角形；否則進(jìn)一步計(jì)算其外接圓圓心與半徑，運(yùn)用LOP法則進(jìn)行優(yōu)化，同時(shí)根據(jù)三角形單元分類規(guī)則更新該三角形的分類標(biāo)識。
　　三角形單元分類規(guī)則是建立在散亂點(diǎn)集有序化基礎(chǔ)上的。假設(shè)散亂點(diǎn)集有序化過程中以X軸的坐標(biāo)分量為比較依據(jù)，則定義分類規(guī)則如下：
　　定義 已知插入點(diǎn)V_i(x_i，y_i，z_i)(i=1，2，3，……，n，n為點(diǎn)的個(gè)數(shù))和待分類三角形T_j(p_1j，p_2j，p_3j)(j=1，2，3，……，m，m為三角形的個(gè)數(shù))，其中T_j的外接圓圓心為C_j(c_xj，c_yj，c_zj)、半徑為r_j。若xi-c_xj≥r_j，則三角形T_j標(biāo)識為已完成三角化，否則標(biāo)識為未完成三角化。
　　分類規(guī)則的直觀解釋如圖1所示。若T_j外接圓圓心C_j到經(jīng)過插入點(diǎn)V_i且垂直于X軸的直線MN的距離大于其半徑r_j的長度，則可保證所有X坐標(biāo)分量不小于插入點(diǎn)Vi的點(diǎn)都位于Tj外接圓外。又因?yàn)樯y點(diǎn)集有序化之后按照X坐標(biāo)分量以升序排列，故對于點(diǎn)Vi之后的插入點(diǎn)V_k(k> i)而言，永遠(yuǎn)都不會位于三角形Tj的外接圓內(nèi)，從而Tj可標(biāo)識為已完成三角化。

3  算法比較與分析
　　算法用ANSI C++語言進(jìn)行了編程實(shí)現(xiàn)，微機(jī)環(huán)境為Pentium(R)4，CPU主頻2.4GHz，內(nèi)存512MB，實(shí)驗(yàn)數(shù)據(jù)為程序隨機(jī)生成的單精度型點(diǎn)集，實(shí)驗(yàn)結(jié)果如表1所示。

　　由實(shí)驗(yàn)結(jié)果可以看出，改進(jìn)后的算法比原算法在速度上有數(shù)倍增長，特別是當(dāng)點(diǎn)數(shù)增加時(shí)速度增長的幅度也增大，而占用的內(nèi)存空間只有少量增長。這充分體現(xiàn)了算法的穩(wěn)定性。在實(shí)際運(yùn)行時(shí)對時(shí)間和空間的合理綜合運(yùn)用提高了算法的效率和性能。
4  工程應(yīng)用
　　地震層位模型的可視化已成為油藏勘探領(lǐng)域研究的熱點(diǎn)。將Delaunay三角網(wǎng)建模方法與OpenGL圖形可視化技術(shù)相結(jié)合，即可快速構(gòu)建地震層位模型，實(shí)現(xiàn)地表特征層表面三維顯示，直觀地表達(dá)特征層的空間分布與形態(tài)。
　　OpenGL是當(dāng)今工業(yè)界最先進(jìn)最開放的圖形API，它提供了強(qiáng)大的繪制二維和三維圖形的能力。如提供了大量的圖形變換函數(shù)，可以方便地將三維圖形顯示在屏幕窗口進(jìn)行放大、縮小、旋轉(zhuǎn)等操作；OpenGL還提供了線面消隱、著色關(guān)照、紋理映射和反走樣等技術(shù)的一系列函數(shù)，可以方便地對地層表面進(jìn)行處理，增強(qiáng)圖形的真實(shí)感。
　　本文基于上述方法，以地震波采樣解釋之后的地震數(shù)據(jù)為數(shù)據(jù)源，重構(gòu)了勝利油田某地區(qū)的層位模型如圖2、圖3所示。模型顯示結(jié)果驗(yàn)證了本文的算法是行之有效的。

5  結(jié)  論
　　本文在仔細(xì)研究常規(guī)逐點(diǎn)插入算法的基礎(chǔ)上，對其進(jìn)行了優(yōu)化，使得改進(jìn)之后的算法不但繼承了原算法流程清晰、實(shí)現(xiàn)簡單、占用內(nèi)存較小、可局部動態(tài)更新三角網(wǎng)的優(yōu)點(diǎn)，而且提高了執(zhí)行效率。該算法運(yùn)用于實(shí)際工程中構(gòu)建地震層位模型，取得了良好的效果。
當(dāng)前，包括地震層位建模在內(nèi)的三維地理信息系統(tǒng)正處于從研究走向?qū)嵱玫倪^渡階段。在實(shí)際應(yīng)用中，有時(shí)還存在著約束三角剖分，也就是部分散亂點(diǎn)之間常常存在著某種約束關(guān)系，如地表模型中的斷裂線等。在對此種數(shù)據(jù)集進(jìn)行三角網(wǎng)生成時(shí)，生成的三角網(wǎng)應(yīng)保持其原有的約束關(guān)系，即約束數(shù)據(jù)集下的三角剖分。因此，如何在約束條件下快速準(zhǔn)確地生成Delaunay三角網(wǎng)，是下一步研究工作的重點(diǎn)。
參考文獻(xiàn)
1   包世泰，夏斌，崔學(xué)軍等.地質(zhì)信息模型研究及其應(yīng)用.大地構(gòu)造與成礦學(xué)，2004；28(4)
2   胡金星，潘懋，馬照亭等.高效構(gòu)建Delaunay三角網(wǎng)數(shù)字地形模型算法研究.北京大學(xué)學(xué)報(bào)(自學(xué)科學(xué)版)，2003；39(5)
3   武曉波，王世新，肖春生.Delaunay三角網(wǎng)的生成算法研究. 測繪學(xué)報(bào)，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)