2三角剖分算法介紹

易讀文庫易讀文庫易讀文庫2三角剖分算法介紹2.1Delaunay三角剖分理論基礎(chǔ)2.1.1VORONOI圖1908年，俄國數(shù)學(xué)家M.G.Voronoi提出了Voronoi圖[1]，它是計算幾何學(xué)上非常重要的工具之一，被廣泛應(yīng)用于各個領(lǐng)域中。Voronoi圖是由平面區(qū)域中連接兩鄰點的線段的中垂線所形成的區(qū)域。它又叫Dirichlet圖或Thiessen多邊形[2]。Voronoi圖是自然界中的宏觀物體和微觀物體以空間距離相互作用所形成的一種網(wǎng)格結(jié)構(gòu)，應(yīng)用的范圍相當(dāng)廣泛，特別是在計算幾何學(xué)領(lǐng)域的應(yīng)用[3]。Voronoi圖是一種關(guān)于平面或空間區(qū)域劃分的基礎(chǔ)數(shù)據(jù)結(jié)構(gòu)。100多年來，它被廣泛應(yīng)用在與幾何信息相關(guān)的各個領(lǐng)域。隨著計算機科學(xué)技術(shù)的不斷普及和發(fā)展，Voronoi圖的應(yīng)用領(lǐng)域也在不斷地擴大。對于Voronoi圖的應(yīng)用，以90年代以來的應(yīng)用更為突出。 Voronoi圖的本質(zhì)屬性是由空間實體幾何唯一確定的，不是通過其他方法強加上去的。從Voronoi圖的數(shù)字幾何角度來看，它是針對平面n個離散點而言的，把平面分為若干區(qū)域，每一個區(qū)域包含一個點，該點所在的區(qū)域就是到該點距離最近的點的集合。 設(shè)p1，p2是平面上兩點，L是線段p1p2的中垂線，L將平面分成兩部分LL和LR，位于LL內(nèi)的點pl具有特性：d(pl,p1)<d(pl,p2)，其中d(pl,pi)表示pl與pi之間的歐幾里得距離，i=1，2。這意味著，位于LL內(nèi)的點比平面上其他點更接近點p1，換句話說，LL內(nèi)的點是比平面上其他店更接近于p1的點的軌跡，記為V(p1)，如圖2.1所示。如果用H(p1,p2)表示半平面LL，而LR=H(p2,p1)，則有V(p1)=H(p1,p2)，V(p2)=H(p2,p1)。 給定平面上n個點的點集S，S={p1,p2,…,pn}。定義： V(pi)=即V(pi)表示比其他點更接近pi的點的軌跡是n-1個半平面的交，它是一個不多于n-1條邊的凸多邊形區(qū)域，稱為關(guān)于pi的Voronoi多邊形或關(guān)于pi的Voronoi域。對于S中的每個點都可以作一個Voronoi多邊形，這樣的n個Voronoi多邊形組成的圖稱為Voronoi圖。Voronoi多邊形的每條邊是S中某兩點的連線的垂直平分線，所有這樣的兩點連線構(gòu)成一個圖，稱為Voronoi圖的直線對偶圖。如圖2.2所示，虛線表示voronoi圖，實線表示其對偶圖。圖2.1平面中兩離散點圖2.2Voronoi圖 Voronoi圖有以下一些性質(zhì)：性質(zhì)1.n個點的點集S的Voronoi圖至多有2n-5個頂點和3n-6條邊。性質(zhì)2.每個Voronoi點點好是三條Voronoi邊的交點。性質(zhì)3.Voronoi圖的對偶圖實際上時點集的一種那個三角剖分，該三角剖分就是Delaunay三角剖分[4]。 性質(zhì)3實際上告訴我們可以用構(gòu)造Voronoi圖的方法來求平面點集的三角剖分，但實際上這種方法較少使用，因為算法的效率不好。2.1.2Delaunay三角剖分的基本概念1、基本概念（1）域分割 給定平面上的n個不相重的散亂數(shù)據(jù)點，對每個散亂數(shù)據(jù)點構(gòu)造一個域，使該域內(nèi)的任一點離此散亂點比離其他散亂點更近，這種域分割就是上節(jié)介紹的Voronoi圖，也稱Dirichlet域分割，由上節(jié)可知，域邊界其實就是連接兩個相鄰散亂點的直線的垂直平分線。 （2）Delaunay三角化 對平面上的散亂數(shù)據(jù)點進(jìn)行域分割后，將具有公共域邊界的散亂點對相連形成的三角剖分稱為Delaunay三角剖分。 （3）優(yōu)化 對三角網(wǎng)格進(jìn)行優(yōu)化，就是要使三角網(wǎng)格整體上盡量均勻，避免出現(xiàn)狹長三角形，也就是獲得Delaunay三角化。2、三角剖分優(yōu)化準(zhǔn)則 在三角剖分過程中，人們往往先用一種簡單的方法構(gòu)造散亂點的初始三角剖分，然后對其進(jìn)行優(yōu)化以獲得Delaunay剖分。優(yōu)化的方法取決于所采用的優(yōu)化準(zhǔn)則，平面三角剖分最常用的優(yōu)化準(zhǔn)則有Thiessen區(qū)域準(zhǔn)則、最小內(nèi)角最大準(zhǔn)則、圓準(zhǔn)則。Sibson證明了這三個準(zhǔn)則的等價性，并指出符合這三個準(zhǔn)則的三角剖分只有一個，即Delaunay三角剖分。最小內(nèi)角最大準(zhǔn)則對一個嚴(yán)格凸的四邊形進(jìn)行三角化時，有兩種選擇，最小內(nèi)角最大準(zhǔn)則就是要保證對角線兩側(cè)兩個三角形中的最小內(nèi)角為最大，如圖2.3所示。圖2.3對角線優(yōu)化圓準(zhǔn)則嚴(yán)格凸四邊形中的三個頂點確定一個圓，如果第四個頂點落在圓內(nèi)，則將第四個頂?shù)着c其相對的頂點相連，否則將另外兩個頂點相連，這個準(zhǔn)則稱為圓準(zhǔn)則。也就是說，符合圓準(zhǔn)則的三角剖分中，任一三角形的外接圓內(nèi)不應(yīng)該包含其他點，如圖2.4所示。圖2.4空外接圓準(zhǔn)則局部優(yōu)化局部優(yōu)化是指對任意一個凸四邊形的對角線，依據(jù)某種優(yōu)化準(zhǔn)則做交還測試后所得到的三角剖分。全局優(yōu)化當(dāng)三角形剖分T中每一條內(nèi)邊上的兩個三角形所形成的凸四邊形都滿足局部優(yōu)化標(biāo)準(zhǔn)時，稱該三角剖分T滿足全局優(yōu)化。退化前面已經(jīng)講過Delaunay三角剖分滿足圓準(zhǔn)則，即任一三角形的外接圓內(nèi)不能包含其他的點，如果規(guī)定三角剖分中任一三角形的外接圓內(nèi)和外接圓上都不能有其余頂點，則稱為標(biāo)準(zhǔn)的Delaunay三角剖分，在實際應(yīng)用中，往往存在四點共圓的情況。我們把四點或四點以上的散亂數(shù)據(jù)點共圓的情況稱為退化，相應(yīng)的三角剖分稱為準(zhǔn)Delaunay三角剖分。3、二維任意區(qū)域內(nèi)點集的Delaunay三角劃分概念及基本定理 將所有共邊的Voronoi多邊形的中心連接起來所形成的三角網(wǎng)稱為Delaunay三角剖分[5]。 在區(qū)域D內(nèi)的點集V的三角劃分（記為T(D,V)）是指將V中的點用互不相交的直線段連接起來，使得域D內(nèi)的每一個區(qū)域都是三角形，同時這些連線和三角形都在區(qū)域D內(nèi)，包含在區(qū)域D內(nèi)的邊稱為內(nèi)邊。 二維任意區(qū)域點集的Delaunay三角劃分[5]（簡記為DATA）是指所有的內(nèi)邊都滿足局部優(yōu)化的T(D,V)。4、基本定理 下面的定理是閔衛(wèi)東在文獻(xiàn)[6]中提出來的，其證明參見文獻(xiàn)[6]。定理2.1：設(shè)B，P點在AC邊的同側(cè)，如果角APC小于角ABC，則三角形ABC的外接圓內(nèi)不包含點P，參見圖2.5。圖2.5三角形外接圓圖2.6局部優(yōu)化邊AB定理2.2：設(shè)AB邊為三角形ABC和三角形ABE所共有，如果四邊形ACBE的內(nèi)角CBE≥180°或CAE≧180°，則AB邊是局部優(yōu)化的，見圖2.6。定理2.3：設(shè)三角形ABC是DATA（D，V）的一個三角形，P是V中的一點且不是ABC的三個頂點，如果線段PC在域D內(nèi)且它與AB邊有交點，則ABC的外接圓內(nèi)部不包含點P。定理2.4：（DATA的Circle準(zhǔn)則）三角剖分T（D，V）是DATA（D，V）的充要條件是，對于T（D，V）的任何一個三角形ABC，V中的任何點P，如果它與A，B，C相連的線段PA，PB，PC都在域D內(nèi)，則P不包含在三角形ABC的外接圓內(nèi)部。定理2.5：任何一個T（D，V）都可以通過有限次的局部優(yōu)化操作將其轉(zhuǎn)化為DATA（D，V）。定理2.6：DATA（D，V）在所有的T（D，V）中平均形態(tài)比最大。 以上幾個定理實際上作為Delaunay三角劃分呢算法的理論依據(jù)在許多以Delaunay原理為基礎(chǔ)的網(wǎng)格剖分算法中經(jīng)常用到。5、三角剖分時需要滿足的幾個約束條件 在進(jìn)行網(wǎng)格剖分過程中，生成三角形時，必須滿足一下三個約束條件：三角形相互之間是不相交的，即兩個三角形除端點外不應(yīng)該有別的交點，此條件稱為三角形相交約束。三角形相互之間是互不包含的，即任意一個三角形不能完全包含其他的三角形，此條件稱為三角形包含約束。三角形全部落在區(qū)域之內(nèi)，此條件稱為三角形合法約束。2.2Delaunay三角剖分算法介紹Delaunay三角剖分算法是最常用的一種平面點集的剖分算法，實際上，Delaunay三角剖分算法是基于Delaunay原理的一類三角剖分算法的總稱。總體來說，Delaunay三角剖分算法可分為以下幾種。2.2.1翻邊算法1977年，Lawson基于最小內(nèi)角最大優(yōu)化準(zhǔn)則提出了一種二維點集Delaunay三角剖分翻邊算法[7]。該算法的思想是：給定一個二維離散點集，首先對其進(jìn)行一次初始三角剖分，然后根據(jù)最小內(nèi)角最大化準(zhǔn)則判斷初始三角剖分中形成凸四邊形的共邊三角形是否滿足優(yōu)化準(zhǔn)則；如果不滿足，就交換四邊形的兩條對角線。一次循環(huán)直到所偶的三角形都滿足最小角最大優(yōu)化準(zhǔn)則為止。通過計算可以得到，翻邊算法的時間復(fù)雜度為O(n2)。根據(jù)前面的理論分析知道，運用翻邊算法撲粉出來的所有三角形邊都是Delaunay三角剖分的邊，那么最終得到的三角網(wǎng)格一定是Delaunay三角網(wǎng)格。 目前，二維的翻邊算法發(fā)展比較完善，利用翻邊法對二維離散點集進(jìn)行Delaunay三角剖分原理比較簡單、得到的網(wǎng)格質(zhì)量也較好。1989年，BarryJoe基于空外接球準(zhǔn)則提出了一種三維Delaunay三角剖分的翻面法[8]。翻面法的原理跟翻邊法的原理是一樣的，只不過在三維空間利用翻面法會比較復(fù)雜，有待進(jìn)一步研究。2.2.2逐點插入算法隨著Delaunay三角剖分的不斷發(fā)展，大量剖分算法不斷被提出來。1981年，一種逐點插入算法被Bowyer和Watson提出來[9][10]。該算法的基本思想是：首先構(gòu)造點集的一個初始三角剖分，然后逐一地在當(dāng)前三角剖分中插入一個新點；在插入過程中，要根據(jù)Delaunay三角剖分空外接圓住著呢進(jìn)行三角網(wǎng)格優(yōu)化，直到整個點集為空集為止。 該算法首先構(gòu)造一個極大三角形，使得要剖分的點集中所有點都落在該三角形里面；其次，逐一地插入各個新點，每插入一個點，必須搜索位于外接圓內(nèi)部的點的三角形；然后，將這些三角形從三角形隊列中刪除，可以形成多邊形空腔；最后，將插入點和多邊形空腔連接起來，構(gòu)成若干個以插入點為共同頂點的新三角形。 定位新插入點的位置是逐點插入算法的一個關(guān)鍵步驟，直接關(guān)系到剖分網(wǎng)格的好壞和剖分速度的快慢。由于包含插入點的三角形的外接圓位于插入點的附近區(qū)域，只要能夠準(zhǔn)確定位新插入點在當(dāng)前三角剖分中的位置，便可快速檢索到外接圓包圍它的三角形?；诙ㄎ徊迦朦c位置的重要性，現(xiàn)在大量改進(jìn)的逐點插入算法主要都體現(xiàn)在定位新插入點的位置上。 由于逐點插入算法比較簡單、易懂，可以將其推廣到三維或者更高維空間的三角剖分。唯一不足的就是該算法時間復(fù)雜度較高，很多專家學(xué)者在研究算法的時候也會將該因素重點考慮進(jìn)去，從而設(shè)計出簡單、高效的三角剖分算法。2.2.3分割合并算法Hoey和Shamos提出了一種用于生成Voronoi圖的分割合并算法，由于Voronoi圖與Delaunay圖互為對偶圖，進(jìn)而可以間接得到Delaunay三角網(wǎng)格[11]。后來Dwyer、Katajaine、L.J.Guibas等人提出了用于直接生成Delaunay三角網(wǎng)格的分割合并算法。 分割合并算法主要根據(jù)數(shù)學(xué)遞歸思想，遞歸將點集分割成規(guī)模相當(dāng)?shù)膬刹糠肿蛹?，然后對分割出來的點集進(jìn)行Delaunay三角剖分，再遞歸地將兩個相鄰Delaunay三角剖分進(jìn)行合并組合，從而生成了整個點集的Delaunay三角網(wǎng)格。改進(jìn)后的分割合并算法時間復(fù)雜度為O(nlogn)，目前最好的分割合并算法是Dewall算法。該算法不僅可以應(yīng)用于二維空間，還可對三維或更高維空間進(jìn)行三角剖分，雖然分割合并算法實現(xiàn)的時間效率很高，但它的編程卻非常復(fù)雜。 以上三種算法都是比較典型的Delaunay三角剖分算法，應(yīng)用范圍比較廣泛。除了這幾種經(jīng)典算法外，其實還有很多種其他的Delaunay三角剖分算法。在選擇三角剖分算法的時候，必須從難易程度、效率高低和貼近實際與否等方面選擇，其中逐點插入算法是目前最為簡單、最為流行的一種Delaunay三角剖分算法，該算法思想簡單、易懂，可推廣到維數(shù)更高空間的三角剖分。2.2.4Bowyer-Watson算法Bowyer-Watson算法是一種先形成粗原始網(wǎng)格，再逐步插入新點進(jìn)行細(xì)化的算法，是一種逐點插入算法。假定已有初始剖分為T，向已有網(wǎng)格內(nèi)插入新點的Bowyer-Watson算法如下： Step1：插入一個新點P到現(xiàn)有的Delaunay三角剖分中。 Step2：尋找并刪除所有外接圓包含P點的三角元素，形成一個Delaunay空腔。 Step3：連接P點和Delaunay腔壁上所有各點，形成新的Delaunay三角剖分。 應(yīng)用Bowyer-Watson算法進(jìn)行網(wǎng)格剖分的整個過程簡述如下：Step1：家里一個覆蓋整個計算區(qū)域的粗原始網(wǎng)格。Step2：采用Bowyer-Watson算法將邊界點逐點插入原始網(wǎng)格。邊界點由人工事先給定，并假設(shè)邊界點的分布式合理的。Step3：刪除計算區(qū)域以外的三角形，并確保邊界面的正確三角剖分（進(jìn)行拓?fù)湎嗳菪蕴幚恚?。Step4：用Bowyer-Watson算法逐步向計算區(qū)域內(nèi)插入新點（對于給定點集，按一定順序插入即可，若內(nèi)點是自動生成的，則需要按一定的策略生成內(nèi)點，直到所有內(nèi)點都被插入或者網(wǎng)格達(dá)到一定的撲粉要求為止。Step5：對所生成的網(wǎng)格進(jìn)行拓?fù)湎嗳菪詸z查、網(wǎng)格光順等工作。3礦井涌水量地質(zhì)模型的實現(xiàn)3.1初始粗網(wǎng)格的建立礦井涌水量地質(zhì)模型的建立，首先是要將煤層頂?shù)装宓那骐x散成三角形網(wǎng)格。本文中采用逐點插入法作為主要的算法。 逐點插入法的第一步是要建立一個包含所有計算域頂點的初始粗網(wǎng)格。常用的建立初始網(wǎng)格的方法是包容盒法，一般可建立兩種形式的包容盒，一種是建立矩形的包容盒[12]，另一種是三角形包容盒[13]。無論是給定點集或者是給定區(qū)域的三角剖分均可采用該方法。如果只給出邊界區(qū)域，則內(nèi)點需要自動產(chǎn)生。本文所建立的三角剖分是給定區(qū)域邊界，自動生成內(nèi)點。這種方法在使用過程中需要對生成的三角形進(jìn)行邊界拓?fù)湎嗳菪詸z測。因為拓?fù)湎嗳菪詸z測的步驟，本文給出一種無需建立三角形粗網(wǎng)格或矩形粗網(wǎng)格的方法，直接將給定區(qū)域的多邊形劃分成三角形初始粗網(wǎng)格。 該方法首先將給定邊界按順時針或者逆時針的方向?qū)⑦吔珥旤c進(jìn)行排列，然后計算每一個頂點的角度大小，作為判斷頂點角度的指標(biāo)。隨后每次選擇最小的頂點，將其與相鄰的頂點連成一個三角形，重新計算該頂點兩側(cè)的頂點角度大小，然后再選擇最小角進(jìn)行處理。 由于多邊形并非是完全的凸多邊形，邊界上很可能會出現(xiàn)凹頂點。在這種情況下，需要判斷頂點的凹凸性。 本文使用矢量面積法[14]來判斷多邊形頂點的凹凸性。如果多邊形邊界以逆時針方向排列，那么為了判斷某一頂點的凹凸性，可以用與該頂點相鄰的兩個頂點，共三個頂點來判斷該頂點的凹凸性。計算這三個頂點形成的三角形的矢量面積，如果結(jié)果為正，那么該頂點為凸頂點，如果為負(fù)，那么該頂點為凹頂點。計算公式可簡單描述如下： 設(shè)三個頂點為：A、B、C，那么： 判斷完頂點的凹凸性之后，需要計算頂點的角度（弧度）?？梢允褂孟噜徣齻€頂點形成的兩個向量來計算頂點的角度。、的數(shù)量積可以計算出兩向量的夾角，再結(jié)合頂點凹凸性的判斷，加上或者減去180°便是該頂點的角度。Delaunay剖分的初始化步驟可用以下算法描述：算法3.1Step1：判斷頂點凹凸性，計算每一個頂點的角度。Step2：從頂點數(shù)組中找到角度最小的頂點，連接該頂點兩側(cè)的頂點，形成一個三角形。并計入三角形數(shù)組。Step3：重新計算該頂點兩側(cè)新頂點的角度，并從頂點數(shù)組中刪去該頂點。Step4：從新的頂點數(shù)組中找到新的最小角度頂點，重復(fù)step2至step3。直到頂點數(shù)組中的頂點個數(shù)等于3個。Step5：連接最后三個頂點形成一個三角形，計入三角形數(shù)組。按照以上算法，最后會形成一個初始的粗三角形網(wǎng)格，或者稱為多邊形的一個粗三角劃分。如圖3.1，3.2。圖3.1任意多邊形區(qū)域圖3.2初始網(wǎng)格劃分3.2簡單多邊形的Delaunay剖分在建立了多邊形的初始粗網(wǎng)格以后，便可以開始向計算區(qū)域內(nèi)插入點進(jìn)行Delaunay三角剖分。本文采用自動布點技術(shù)向計算區(qū)域內(nèi)插入點。采用的方法是：選擇面積最大三角形的最長邊的中點最為插入點，建立三角網(wǎng)格。這樣需要設(shè)置一個面積的控制量，以便控制剖分的程度。即三角形的個數(shù)或者三角網(wǎng)格的疏密程度。這是一個可以改變的量，可以方便地控制三角網(wǎng)格的大小。當(dāng)插入一個點以后，需要搜索整個三角形數(shù)組，進(jìn)行空外接圓檢測，如果新插入點落在該三角形的外接圓內(nèi)（包含邊界），那么該三角形被添加到需要優(yōu)化的三角形數(shù)組中，所有三角形檢測完后，符合要求的三角形會形成一個多邊形空腔，刪除這些三角形的公共邊，連接插入點與多邊形空腔的各個頂點，形成新的多個三角形，加入三角形數(shù)組，刪除被合并的三角形。這樣一直到三角形面積符合控制要求結(jié)束。便完成了簡單多邊形的Delaunay剖分。算法可描述如下：算法3.2 Step1：設(shè)定三角形面積控制單元，計算三角形數(shù)組中所有三角形的面積大小。 Step2：選擇面積最大的三角形的最長邊作為新的插入點，并將該點加到頂點數(shù)組的后面。 Step3：搜索整個三角形數(shù)組，找出所有插入點落在三角形外接圓中的三角形，計入優(yōu)化三角形數(shù)組，并將三角形的邊計入邊數(shù)組，并且不能有重復(fù)邊，若有重復(fù)邊，便要從數(shù)組中刪除該邊。Step4：連接插入點與邊數(shù)組中每條邊的兩個頂點，形成新的三角形，加入三角形數(shù)組。形成新的三角形數(shù)組。 Step5：刪除優(yōu)化三角形數(shù)組中的三角形，更新整個三角形數(shù)組，完成一次點的插入。 Step6：重復(fù)step2、step3、step4，直到三角形的面積符合控制要求，結(jié)束點的插入，這樣就完成了簡單多邊形的Delaunay剖分。剖分實例如圖3.3。圖3.3Delaunay剖分在剖分的過程中需要注意幾個問題。如果插入點落在多邊形邊界上，便會形成同一直線上的兩條不同的邊被同時計入邊數(shù)組，當(dāng)連接三角形的時候，便會出現(xiàn)三個頂點落在一條直線上的三角形，所以需要對這種三角形進(jìn)行特別的處理。不然會出現(xiàn)錯誤。處理方法很簡單，在連接三角形的時候，需要計算插入點與兩個頂點連接形成的直線的斜率，比較這兩條直線的斜率，如果斜率相等，那么表示這三個點處在同一條直線上，這樣便不能連接三角形。算法表示如下：算法3.3 Step1：計算邊數(shù)組中將要連接的那條邊的兩個頂點與插入點形成的直線的斜率。 Step2：判斷兩條直線的斜率。 Step3：如果斜率相等，那么那條邊的連接性置為0，如果斜率不等那么直線的連接性置為1。Step4：連接邊的兩個頂點與插入點，形成一個三角形，計入三角形數(shù)組。在計算斜率的時候，由于數(shù)據(jù)類型的限制，三角形的斜率不一定嚴(yán)格相等，所以當(dāng)三點一線的時候斜率會有微小的誤差，這一誤差需要用一常量來作為判斷斜率是否相等的限制，可以將這一誤差設(shè)定為1°，用來幫助判斷斜率。3.3用VB實現(xiàn)的Delaunay剖分算法 本文使用VisualBasic6.0作為開發(fā)環(huán)境來實現(xiàn)Delaunay剖分算法。分為四個模塊：公共數(shù)據(jù)模塊初始化模塊剖分模塊數(shù)據(jù)輸出模塊3.3.1公共數(shù)據(jù)模塊 公共數(shù)據(jù)模塊定義有程序中使用的數(shù)據(jù)類型，包括：頂點類，存儲有頂點的x、y、z坐標(biāo)數(shù)據(jù)，以及頂點是否是邊界格點的判斷變量。三角形類，存儲有每個三角形的三個頂點編號，編號既是三角形頂點在頂點數(shù)組中的序號。臨時頂點類，該數(shù)據(jù)類型主要用于初始化時存儲臨時頂點數(shù)據(jù)，除了包括x、y、z坐標(biāo)外，還有弧度值與該頂點在原頂點數(shù)組中的序號。邊類，該數(shù)據(jù)類主要用于剖分時存儲被優(yōu)化三角形的邊信息，數(shù)組中的邊將成為新生成的三角形的邊。模塊中還定義有全局頂點數(shù)組和全局三角形數(shù)組，全局頂點個數(shù)和全局三角形個數(shù)。全局頂點數(shù)組是一個只增不減的數(shù)組，頂點數(shù)會隨著剖分的進(jìn)行而不斷增加，且數(shù)據(jù)不會改變。全局三角形數(shù)組是一個動態(tài)數(shù)組，數(shù)組元素的個數(shù)和元素的取值都會隨著程序的運行而改變。初始化時已事先為兩個數(shù)組定義了元素個數(shù)。可隨剖分的細(xì)分程度修改。3.3.2初始化模塊 初始化模塊用于對多邊形區(qū)域進(jìn)行初始的網(wǎng)格化。形成初始的三角形存入三角形數(shù)組。該模塊中先定義臨時頂點數(shù)組，該臨時頂點數(shù)組繼承全局頂點數(shù)組的全部數(shù)據(jù)。并存儲每一個頂點在全局頂點數(shù)組中的序號。 隨后計算每一個頂點的角度，存入臨時頂點數(shù)組，使用算法3.1判斷頂點的凹凸性并計算頂點的角度。經(jīng)測試矢量面積法能判斷頂點的凹凸性，由此可以計算每一個頂點的角度大小，在程序中，角度值轉(zhuǎn)換為弧度值計算。 該模塊中包括判斷頂點凹凸性的函數(shù)，初始化網(wǎng)格的過程和繪制網(wǎng)格的過程。這三個函數(shù)或過程在初始化過程中被調(diào)用。 初始化網(wǎng)格的過程中，先找到頂點角度最小的角，將該頂點兩側(cè)的頂點連接形成一個三角形，插入全局三角形數(shù)組。并在臨時頂點數(shù)組中刪除該頂點。重新計算該頂點兩側(cè)的頂點角度大小，并更新臨時頂點數(shù)組。在計算的過程中，需要考慮序號為1、2、n-1、n這四個頂點。因為數(shù)組數(shù)據(jù)類型的限制，當(dāng)最小角頂點序號是這四個值時，頂點兩側(cè)的頂點序號并不連續(xù)，所以需要分開計算。用判斷語句分別判斷這四種情況，分別計算頂點兩側(cè)的頂點角度。其余情況可用循環(huán)語句計算。執(zhí)行完初始化網(wǎng)格過程，全局三角形數(shù)組包含最初的三角形網(wǎng)格，全局頂點數(shù)組不變。調(diào)用繪制網(wǎng)格過程繪制三角形。3.3.3剖分模塊該模塊用于將多邊形剖分成需要的Delaunay三角形網(wǎng)格。模塊首先建立邊數(shù)組，用于存儲需要連接的邊信息。首先計算現(xiàn)有三角形的面積，找到面積最大的三角形。然后取該三角形最長邊的中點作為新的插入點，將插入點計入全局頂點數(shù)組。搜索整個全局三角形數(shù)組，找到插入點落在其外接圓中的三角形，并記錄該三角形在全局三角形數(shù)組中的序號，將其三條邊送入邊數(shù)組等待優(yōu)化。在將邊送入邊數(shù)組的過程中，需要判斷邊是否已經(jīng)存在邊數(shù)組中，如果已經(jīng)存在則要在邊數(shù)組中刪除該邊，如果不存在則存入邊數(shù)組。搜索完全部三角形后，形成了一個記錄需要優(yōu)化三角形的序號數(shù)組和邊數(shù)組，邊數(shù)組中的邊形成一個連續(xù)的多邊形空腔，插入點可能位于空腔內(nèi)部，或者空腔的某一條邊上。如果位于邊上，則需要檢測插入點落在哪條邊上，而這條邊則不能作為新三角形的邊。使用插入點與邊兩個頂點的連線的斜率來判斷插入點是否落在邊上，這樣便可除去三點一線的三角形。按順序連接插入點與邊數(shù)組中的每條邊，形成新的數(shù)個三角形，將這些三角形送入全局三角形數(shù)組。如果符合要求的三角形只有一個，則要專門計算插入三角形的個數(shù)，形成兩個新的三角形，刪除一個舊的三角形。隨后檢測需要刪除的三角形，從數(shù)組尾端開始向前逐個刪除需要刪除的三角形。這樣便完成了一個點在多邊形區(qū)域中的插入，更新了全局頂點數(shù)組和全局三角形數(shù)組。重新計算數(shù)組中每個三角形的面積，找到面積最大的三角形，如果面積值大于某一設(shè)定值，則繼續(xù)向其中插入點，直到滿足要求為止。這樣便完成了多邊形區(qū)域的Delaunay剖分。3.3.4數(shù)據(jù)輸出模塊 該模塊主要用于頂點數(shù)據(jù)和三角形數(shù)據(jù)的輸出，以便在隨后的計算中使用。將需要的數(shù)據(jù)輸出到一個文本文件中，留待后用。程序界面如圖3.4所示：圖3.4程序運行界面3.4對剖分算法的評價該算法可以將簡單的多邊形剖分成Delaunay三角網(wǎng)格。并且算法省去了建立初始大三角形或者矩形的步驟，直接將多邊形劃分成初始粗網(wǎng)格。然后進(jìn)行剖分。這樣就省去了進(jìn)行邊界拓?fù)湎嗳菪詸z查的步驟。和其他算法相比，該算法在時間復(fù)雜度上并無優(yōu)勢，剖分本質(zhì)上是逐點插入法，時間復(fù)雜度較高。 該算法的一個缺點是，當(dāng)三角形長邊位于邊界上時，不能保證該三角形是銳角三角形，也是就該三角形有可能是鈍角三角形。目前尚無有效的辦法解決此問題，這也是選擇布點策略時產(chǎn)生的結(jié)果。3.5網(wǎng)格節(jié)點編號優(yōu)化3.5.1節(jié)點編號優(yōu)化算法有限差分分析時經(jīng)常會碰到具有稀疏矩陣剛度的方程式。減少稀疏矩陣的帶寬可以減少有限差分法的計算時間，和減少占用計算機的存儲空間。在有限差分法的計算過程中，可通過節(jié)點重新編號來減少帶寬。對于節(jié)點重編號的研究較多，最典型的有：Cuthill和Mckee[15]基于圖論的節(jié)點分層編號方法，可在一定程度上減少剛度矩陣的帶寬，這個方法后來被稱為CM方法。Gibbs，Poole和Stockmeyer[16]提出了更加有效的基于圖論的減少稀疏矩陣帶寬和外形的方法，此方法后來被稱為GPS方法。Akhras和Dhatt[17]通過計算節(jié)點商方法來縮小帶寬，該方法編程方便，此方法后來被稱為AD方法。歐陽興[18]提出了前沿法與矩形法。 CM方法和GPS方法，采用了復(fù)雜的圖論方法[19]，沒有估算算法本身的復(fù)雜度。實現(xiàn)比較困難，如果考慮節(jié)點編號算法本身時間消耗大，就會降低節(jié)點編號算法的效率。AD法雖然實現(xiàn)比較方便，但是在網(wǎng)格重劃和網(wǎng)格加密時節(jié)點編號優(yōu)化速度慢。前沿法和矩形法對于板料單元和平面二維單元特別適用。 黃志超，包忠詡，周天瑞[20]等從分析最佳最簡單的矩形出發(fā)，總結(jié)出一種新的節(jié)點重編號方法，是一種改進(jìn)的AD法。 如圖3.5是最典型的最佳二維網(wǎng)格節(jié)點編號，由此來分析其特性，總結(jié)出一般規(guī)律。圖3.5二維網(wǎng)格15節(jié)點優(yōu)化 該15節(jié)點序號是最佳的編號法，它的最優(yōu)同單元節(jié)點編號的最大差值為4。由該圖的編號法可總結(jié)出以下規(guī)則：節(jié)點商是從小到大升序排列的；最小最大節(jié)點編號和也是從小到大升序排列的；最大節(jié)點基本上是按升序排列的，不過有編號重復(fù)的情況。 由此可總結(jié)出一種新的節(jié)點重編號方法。算法3.4：Step1： 計算出每個節(jié)點的相關(guān)單元節(jié)點編號總和；Step2： 計算出該節(jié)點所相關(guān)單元的最大節(jié)點編號與最小節(jié)點編號；Step3： 計算出最小節(jié)點編號與最大節(jié)點編號之和；Step4：計算出該節(jié)點的節(jié)點商，將該節(jié)點的相關(guān)單元節(jié)點編號總和除以該節(jié)點的相關(guān)單元個數(shù)所得的值為節(jié)點商；Step5： 把節(jié)點商作為該節(jié)點重編號的依據(jù)，節(jié)點商小的編號也小，節(jié)點商大的編號也大，如果節(jié)點商相等則把最小最大節(jié)點編號和作為重編號的依據(jù)；step6： 重復(fù)（1）到（5）直到帶寬不能再減小為止。經(jīng)計算該算法能減少50%左右的帶寬。圖3.6為一495個節(jié)點的網(wǎng)格形成的稀疏矩陣每一行的帶寬分布圖。圖3.6節(jié)點帶寬分布圖3.5.2節(jié)點編號優(yōu)化算法的實現(xiàn) 本文使用VB語言實現(xiàn)節(jié)點編號算法。輸入數(shù)據(jù)為剖分完成后的節(jié)點數(shù)組和三角形數(shù)組。該算法對應(yīng)程序的計算水頭模塊。 在公共數(shù)據(jù)模塊中定義一新的數(shù)據(jù)類OptVertex。該數(shù)據(jù)類包含節(jié)點在頂點數(shù)組中的序號number，節(jié)點商NodeTrade和最大最小節(jié)點編號之和maxmin三個數(shù)據(jù)。在計算水頭模塊定義一該類型的數(shù)組作為節(jié)點優(yōu)化數(shù)組，命名為OptVer。數(shù)組長度和全局頂點數(shù)組相同。 算法首先將全局頂點數(shù)組的序號存入節(jié)點優(yōu)化數(shù)組，然后開始循環(huán)。計算節(jié)點優(yōu)化數(shù)組的每一個節(jié)點，搜索全局三角形數(shù)組如果三角形中包含該節(jié)點的編號，便將該節(jié)點相關(guān)單元數(shù)加1，并根據(jù)三角形的三個頂點在全局頂點數(shù)組中的序號從節(jié)點優(yōu)化數(shù)組中找到這三個頂點，并將其在節(jié)點優(yōu)化數(shù)組中對應(yīng)的序號累加存入一變量中。同時要記錄相關(guān)單元中節(jié)點序號的最大值與最小值。三角形數(shù)組搜索完成后，便得到了該節(jié)點的相關(guān)單元數(shù)、相關(guān)單元節(jié)點編號總和、最大節(jié)點編號和最小節(jié)點編號。接著計算該節(jié)點的節(jié)點商和最大最小節(jié)點編號之和，分別存入該節(jié)點的對應(yīng)變量中。這樣依次計算所有節(jié)點優(yōu)化數(shù)組中所有節(jié)點的相關(guān)數(shù)據(jù)存入數(shù)組中。 對計算完成的節(jié)點進(jìn)行排序，使用選擇排序法，將第一個位置的節(jié)點和后面所有的節(jié)點做比較，節(jié)點商大的序號也大，節(jié)點商小的序號也小。如果節(jié)點商相等，比較節(jié)點的最大最小節(jié)點編號之和，小的序號也小，大的序號也大。然后比較第二個位置的節(jié)點和后面的所有節(jié)點，依次排序下去。這樣就完成了一次節(jié)點編號的優(yōu)化。 接著計算每一個節(jié)點的帶寬，初始帶寬為節(jié)點的個數(shù)。然后計算每一個節(jié)點在稀疏矩陣中的帶寬，計算方法為，從三角形數(shù)組中找到包含該節(jié)點的三角形，計算每一個節(jié)點和該節(jié)點之差的絕對值，也就是在稀疏矩陣的一行中包含數(shù)據(jù)的元素到該行對角線的距離的絕對值。計算出最大的帶寬，然后和初始帶寬比較，如果比初始帶寬小，說明稀疏矩陣帶寬減小了，應(yīng)當(dāng)進(jìn)行下一次節(jié)點編號優(yōu)化。將當(dāng)前帶寬作為下一次比較的帶寬，繼續(xù)循環(huán)重復(fù)以上步驟對節(jié)點編號進(jìn)行優(yōu)化排序。直到當(dāng)前帶寬不再改變?yōu)橹埂?這樣就完成了節(jié)點編號的優(yōu)化。3.6三角形網(wǎng)格有限差分法3.6.1三角形均衡網(wǎng)格差分方程的建立三角形網(wǎng)格差分方程是計算網(wǎng)格中格點水頭的重要步驟。也是有限差分方法的核心。差分方程形成的系數(shù)矩陣是主要的數(shù)據(jù)結(jié)構(gòu)。 考慮以格點i為中心的多邊形網(wǎng)格Di作為均衡區(qū)（如圖3.7），下面可根據(jù)達(dá)西定律和水均衡原理建立格點i的差分方程[21]。圖3.7多邊形網(wǎng)格均衡區(qū) 首先求出單位時間內(nèi)通過的周邊流入內(nèi)的水量。 由圖3.7可見，我們可以先求出通過三角形ijk內(nèi)的邊和流入內(nèi)的水量。由達(dá)西定律有：（3.1）式中：、為流段ij和ik的平均導(dǎo)水系數(shù)；、、、為線段長度；為從第e號三角形通過兩線段和流入內(nèi)的水量。根據(jù)（式3.1）的計算方法，對格點i周圍的所有三角形作類似計算，并求和得：（3.2） 另外，由于本文所研究的計算區(qū)域為定水頭邊界，承壓水，不考慮垂向補給量。 所以在單位時間內(nèi)，流入均衡區(qū)內(nèi)的總水量為：Q側(cè)。由水均衡原理得均衡方程：（i=1,2,…N） （3.3）式中：為的面積，為給水度。 這就是多邊形網(wǎng)格的水均衡方程。 由以上公式可以看出，在建立方程時，必須對格點i周圍的每個三角形計算、值，而在通常情況下，只知道格點的坐標(biāo)，因此，為了方便，下面將、以及均用格點坐標(biāo)標(biāo)出。（3.4）（3.5）并且記（3.6）（3.7）其中，=（3.8） 接下來計算多邊形的面積，可以將的計算分解成若干小三角形，分別計算每個小三角形的面積。計算三角形和三角形的面積。（3.9）同理可得（3.10） 類似此法，在格點i周圍的各個三角形中重復(fù)上述工作，即可求得的值。 由上述推導(dǎo)過程可見，要建立格點i的差分方程，只須對格點i周圍的三角形逐一做類似的計算即可。3.6.2三角形網(wǎng)格差分方程的解法 由（式3.2）和（式3.3）可知，在一個差分方程中，將涉及到格點i及其周圍幾個格點的水頭值。在本文的算例中，由于一類邊界格點是定水頭邊界，所以水頭值是已知的。其余內(nèi)格點是未知的，所以必須聯(lián)立方程組求解。 按前述方法，對每個內(nèi)格點建立差分方程，若內(nèi)格點有N1個，則共有N1個方程，而這些方程包含的未知數(shù)也正好是N1個，所以方程數(shù)和未知水頭數(shù)一樣多。此外，如果把格點i的差分方程中的項的系數(shù)排列在系數(shù)矩陣的主對角線上，其系數(shù)矩陣為對角占優(yōu)，故方程組存在唯一解，可用迭代法求解。 SOR（SuccessiveOver-Relaxation）法，即超松弛迭代法[22]，是目前解大型線型方程組的一種最常用的迭代加速方法[23]。它是在高斯-塞德爾迭代法[24]的基礎(chǔ)上進(jìn)行加速的，將前一步的結(jié)果與高斯-賽德爾迭代法的迭代值適當(dāng)?shù)募訖?quán)平均，期望獲得更好的近似值。其迭代公式[25][26]如下： i=1,2,……,n;k=0,1,2,… (3.11) SOR法中ω[27]的取值對迭代公式的收斂速度影響很大[28]，它的好壞直接影響到加速的快慢。為了保證迭代過程的收斂，必須要求0<ω<2，超松弛迭代法取1<ω<2[29]。算法3.4Step1： 將內(nèi)格點原始編號存入一個新的數(shù)組，獲得系數(shù)矩陣大小。Step2： 搜索三角形數(shù)組，找到包含當(dāng)前內(nèi)格點的三角形，計算當(dāng)前內(nèi)格點的多邊形區(qū)域面積，周圍格點的系數(shù)，存入數(shù)組中，直到所有內(nèi)格點都處理完。Step3： 將數(shù)組中的系數(shù)數(shù)據(jù)按照內(nèi)格點順序存入系數(shù)矩陣中，如果格點為邊界格點，將其計算結(jié)果放入右端已知項矩陣中。Step4： 使用SOR法解線性方程組，得到水頭值，存入數(shù)組中。3.6.3計算水頭模塊 三角形網(wǎng)格差分方程的解法在程序的“計算水頭”模塊中實現(xiàn)。該模塊前半部分為網(wǎng)格節(jié)點編號的優(yōu)化。后半部分為三角形網(wǎng)格有限差分法的計算。該部分主要的程序代碼為系數(shù)矩陣的形成。由于研究區(qū)域為定水頭邊界，水頭值已知，且不隨時間的變化而變化。所以這部分格點應(yīng)作為已知量，放到線性方程組右端已知項中。因此，用于計算的未知水頭值個數(shù)即為內(nèi)格點的個數(shù)。在已經(jīng)優(yōu)化的節(jié)點編號中需要除去邊界格點。僅保留內(nèi)格點。重新計算過的節(jié)點數(shù)組大小將只是內(nèi)格點的個數(shù)。以此數(shù)組的大小來定義系數(shù)矩陣，右端已知項矩陣，初始水頭值數(shù)組的大小。計算中包含的重要數(shù)據(jù)結(jié)構(gòu)包括：Matrix（points，points），即系數(shù)矩陣，BMatrix（points），即右端已知項矩陣，x（points），即初始水頭值向量，其中points即為內(nèi)格點的個數(shù)。以下描述程序的具體實現(xiàn)過程。首先形成新的節(jié)點數(shù)組，該數(shù)組不包含邊界格點，僅包含內(nèi)格點。這樣，也即確定了系數(shù)矩陣的大小，重新定義系數(shù)矩陣的大小，方程右端已知項的大小。搜索全局三角形數(shù)組，找到包含當(dāng)前節(jié)點的三角形，逐個計算當(dāng)前節(jié)點周圍節(jié)點的系數(shù)，存入當(dāng)前節(jié)點的系數(shù)數(shù)組中，計算當(dāng)前節(jié)點的多邊形區(qū)域面積，保存在變量中，待三角形全部搜索完畢，即獲得當(dāng)前節(jié)點多邊形區(qū)域的面積和周圍各個節(jié)點的系數(shù)。將數(shù)組中的各個節(jié)點的系數(shù)按照節(jié)點編號順序存入系數(shù)矩陣對應(yīng)位置，如果該節(jié)點為邊界節(jié)點，便獲得其水頭值，將水頭值與系數(shù)相乘加到右端已知項中。隨后將已知項存入已知項矩陣對應(yīng)位置。循環(huán)計算內(nèi)格點的系數(shù)和右端已知項，存入系數(shù)矩陣和右端已知項矩陣對應(yīng)位置，這樣便完成了系數(shù)矩陣和右端已知項矩陣的形成。取得各個格點的初始水頭值，存入初始向量數(shù)組中，將其和系數(shù)矩陣以及右端已知項矩陣作為SOR法的輸入數(shù)據(jù)。調(diào)用SOR法解方程組，迭代完成后即得到該時段末各個內(nèi)格點的水頭值。將其存入相應(yīng)數(shù)組中。以上，便是按照水均衡法原理[30]設(shè)計的礦井涌水量程序，程序獲得水頭值以后便可以交給相關(guān)部門預(yù)測涌水量，以及進(jìn)行礦井水害防治等工作。6 10-108工作面涌水量計算本文以霍州礦區(qū)白龍煤礦480水平，10-108工作面作為研究區(qū)域。將區(qū)域做了簡化，具體研究區(qū)域見圖6.1。并且假定工作面已經(jīng)形成，然后開始涌水，工作面不隨時間變化。圖6.1研究區(qū)域邊界及工作面位置外部多邊形區(qū)域為本文研究區(qū)域的邊界，邊界為定水頭邊界，即水頭不隨時間變化，為一定值550。含水層為K2灰?guī)r，厚度為10m，給水度1/10000，滲透系數(shù)K為2m/天。內(nèi)部矩形區(qū)域為含水層下部巷道的位置，含水層中的水即會從矩形區(qū)域向下涌入巷道，該區(qū)域內(nèi)部也為一定水頭邊界，為一定值480。滲透系數(shù)K為5m/天，下滲路徑長度為5m，矩形區(qū)域面積為100144m2。涌水量計算這也便是本文所要研究的主要內(nèi)容。首先對研究區(qū)域進(jìn)行初始化及剖分。結(jié)果如圖6.2、6.3所示。圖6.2研究區(qū)初始劃分圖6.3研究區(qū)Delaunay剖分 隨后建立差分方程組計算各個內(nèi)部格點的水頭值，然后應(yīng)用達(dá)西定律得到涌水量。達(dá)西定律如下： （6.1）式中：Q為涌水量，K為滲透系數(shù)，F(xiàn)為過水?dāng)嗝?，t為時間，h為總水頭損失，L為滲流路徑長度。 計算中得到的一些數(shù)據(jù)：多邊形區(qū)域邊界點個數(shù)：12工作面面積：100144m2剖分節(jié)點總個數(shù)：661三角形單元個數(shù)：1236內(nèi)部節(jié)點個數(shù)：649工作面內(nèi)節(jié)點個數(shù)：85 應(yīng)用程序得到的涌水量曲線如圖6.4，x軸為時間