物探數(shù)據(jù)三角網(wǎng)逆生長(zhǎng)網(wǎng)格化方法

李瑞雪，張道軍，黃 龍，席振銖，王 鶴，馮萬(wàn)杰

(1. 中南大學(xué) 有色金屬成礦預(yù)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410083；2. 中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083；2. 長(zhǎng)沙五維地科勘察技術(shù)有限責(zé)任公司，長(zhǎng)沙 410205)

目前，常用于物探數(shù)據(jù)成圖的軟件有 Surfer、Grapher、AutoCAD 和 MapGIS等，其中 Surfer和Grapher屬于科學(xué)繪圖軟件，有良好的成圖功能但圖件的編輯能力相對(duì)較弱；AutoCAD是輔助設(shè)計(jì)軟件，在圖件編輯方面有極強(qiáng)優(yōu)勢(shì)但不具備物探數(shù)據(jù)成圖功能；MapGIS等是地理信息系統(tǒng)軟件，可用于物探數(shù)據(jù)的處理但通常需要與其他軟件結(jié)合使用，這些軟件都不是針對(duì)物探專業(yè)開(kāi)發(fā)。專用于物探數(shù)據(jù)可視化的軟件開(kāi)發(fā)具有重要意義，而數(shù)據(jù)的網(wǎng)格化是其可視化的基礎(chǔ)，所以，很有必要研究適用于物探數(shù)據(jù)的網(wǎng)格化方法。

所謂網(wǎng)格化就是將平面區(qū)域劃分為規(guī)則或不規(guī)則的格網(wǎng)，然后通過(guò)一定的方法得到網(wǎng)格節(jié)點(diǎn)處的物理量值。網(wǎng)格化方法廣泛應(yīng)用于有限元計(jì)算、計(jì)算機(jī)輔助幾何設(shè)計(jì)及可視化分析。早期出現(xiàn)的網(wǎng)格化主要采用矩形網(wǎng)格化方法，但因其網(wǎng)格單元結(jié)構(gòu)簡(jiǎn)單，在對(duì)復(fù)雜地形的擬合方面靈活性較差，逐步發(fā)展到四邊形網(wǎng)格化方法[1]，1970年后開(kāi)始出現(xiàn)了以三角形為網(wǎng)格單元的網(wǎng)格化方法。目前，Delaunay三角網(wǎng)是公認(rèn)的最優(yōu)三角網(wǎng)，已有許多專家學(xué)者對(duì)此進(jìn)行了深入的研究和相應(yīng)的改進(jìn)[2?9]。1975年到 1978年，SHAMOS和 HOEY[2]、LAWSON[3]以及 GREEN 和 SIBSON[4]先后提出了分治算法、逐點(diǎn)插入算法及生長(zhǎng)算法等主流三角網(wǎng)格化算法，1987年CHEW[5]提出了Delaunay細(xì)化算法并成功應(yīng)用于均勻二維網(wǎng)格，改進(jìn)了網(wǎng)格質(zhì)量，但未提及邊界的擬合問(wèn)題，1998年楊欽等[6]解決了平面任意域離散點(diǎn)集的三角網(wǎng)格化問(wèn)題，解決了邊界擬合的問(wèn)題，但計(jì)算過(guò)程復(fù)雜。

在地球物理學(xué)領(lǐng)域，國(guó)內(nèi)外研究學(xué)者越來(lái)越關(guān)注三角網(wǎng)格化方法[10?15]。1986年，WANNAMAKER等[10]在研究有限元二維電磁響應(yīng)時(shí)用三角網(wǎng)格化方法剖分了帶地形的二維模型，1999年，阮百堯等[11]采用三角網(wǎng)格對(duì)激發(fā)極化數(shù)據(jù)進(jìn)行了二維反演，2000年，李桐林等[12]研究了規(guī)則邊界地震模型的三角網(wǎng)格化方法，2006年，湯井田等[13]介紹了任意復(fù)雜地球物理模型的三角剖分?？梢?jiàn)，三角網(wǎng)格化方法在地球物理學(xué)領(lǐng)域的應(yīng)用集中在正反演計(jì)算的模型剖分，而應(yīng)用到物探數(shù)據(jù)的可視化分析則較少。例如，2010年，劉兆平等[16]對(duì)于物探數(shù)據(jù)網(wǎng)格化的研究仍集中在矩形網(wǎng)格化方法上。

基于上述分析，本文作者擬在三角網(wǎng)生長(zhǎng)法的基礎(chǔ)上，充分考慮物探數(shù)據(jù)既有的分布離散但不趨于平均以及邊界點(diǎn)明確等特點(diǎn)，提出一種三角網(wǎng)逆生長(zhǎng)網(wǎng)格化方法。有別于傳統(tǒng)的三角網(wǎng)生長(zhǎng)法首先生成一個(gè)三角形作為生長(zhǎng)基礎(chǔ)再向外擴(kuò)展構(gòu)建三角網(wǎng)，本算法在構(gòu)網(wǎng)時(shí)首先計(jì)算邊界，再由邊界作為基礎(chǔ)進(jìn)行擴(kuò)展，通過(guò)距離和最小原則保證三角網(wǎng)質(zhì)量，自外向內(nèi)逆向生長(zhǎng)構(gòu)建三角網(wǎng)。

1 三角網(wǎng)逆生長(zhǎng)網(wǎng)格化方法原理

三角網(wǎng)逆生長(zhǎng)網(wǎng)格化方法區(qū)別于普通生長(zhǎng)法的特點(diǎn)是構(gòu)網(wǎng)時(shí)先處理邊界。其實(shí)現(xiàn)步驟是首先根據(jù)物探數(shù)據(jù)的特點(diǎn)提取邊界點(diǎn)，生成邊界邊；然后以這些邊界邊為基礎(chǔ)，逐一選取其最佳第三點(diǎn)，構(gòu)建新邊及三角形；再以新生成的邊為基礎(chǔ)重復(fù)上述過(guò)程，直到所有離散點(diǎn)都成為三角網(wǎng)的節(jié)點(diǎn)，整個(gè)三角網(wǎng)構(gòu)建完畢。

1.1 地球物理數(shù)據(jù)邊界提取

三角網(wǎng)逆生長(zhǎng)網(wǎng)格化方法是以邊界邊為基礎(chǔ)進(jìn)行生長(zhǎng)的，在構(gòu)網(wǎng)時(shí)必須首先提取邊界。

地球物理數(shù)據(jù)的每一個(gè)數(shù)據(jù)點(diǎn)都包括坐標(biāo)信息及屬性信息。物探測(cè)深數(shù)據(jù)在縱向上成列分布，橫向上與地形及地下構(gòu)造相關(guān)，一般淺部數(shù)據(jù)密集，深部相對(duì)稀疏，在整個(gè)空間上分布離散且不均勻；一條測(cè)線上的邊界數(shù)據(jù)點(diǎn)包括最小號(hào)、最大號(hào)測(cè)點(diǎn)的全部數(shù)據(jù)以及其他號(hào)測(cè)點(diǎn)的最淺部和最深部數(shù)據(jù)，這些數(shù)據(jù)點(diǎn)信息明確，易識(shí)別提取。對(duì)于物探平面數(shù)據(jù)，其規(guī)律更加明確，邊界數(shù)據(jù)僅由測(cè)線測(cè)點(diǎn)號(hào)即能識(shí)別。邊界點(diǎn)的提取可以由任意一點(diǎn)開(kāi)始，為了便于生成邊界邊及擴(kuò)展邊界邊，確定第一個(gè)邊界點(diǎn)后，按照順時(shí)針或逆時(shí)針的順序順次提取相鄰的邊界點(diǎn)；得到全部邊界點(diǎn)并依次連接，從而得到整個(gè)離散數(shù)據(jù)集的閉合邊界。

1.2 點(diǎn)線位置關(guān)系判斷

對(duì)三角網(wǎng)中的任意一邊進(jìn)行擴(kuò)展時(shí)，并不是離散點(diǎn)集內(nèi)的任意一點(diǎn)都有成為該邊最佳第三點(diǎn)的可能，最佳的第三點(diǎn)只存在于一定的范圍內(nèi)。如圖1所示，邊P1P2已存在于三角形P1P2P3中，因此，有可能成為該邊的最佳第三點(diǎn)的點(diǎn)只可能存在于 P1P2所在直線的右側(cè)。在擴(kuò)展過(guò)程中，對(duì)任何一條邊的擴(kuò)展都需要先通過(guò)判斷點(diǎn)線位置關(guān)系求取符合這一范圍要求的點(diǎn)。目前常用的點(diǎn)線位置關(guān)系的判斷通常利用直線判別正負(fù)區(qū)的原理[12]，即首先建立直線判別式，然后將待擴(kuò)展邊的兩端點(diǎn)代入方程判別式計(jì)算系數(shù)，最后將第三點(diǎn)代入直線方程判別式，根據(jù)所得值的符號(hào)判斷該點(diǎn)位于待擴(kuò)展邊的哪一側(cè)。該判別方法計(jì)算過(guò)程比較復(fù)雜，而利用向量方法則可以大大簡(jiǎn)化這一過(guò)程。

圖1 向量叉積判斷點(diǎn)線位置關(guān)系Fig. 1 Judgment of position relationship of point and edge using vector cross product

對(duì)于某一內(nèi)部邊，需要判斷一點(diǎn)是否位于該邊所在三角形的對(duì)頂點(diǎn)的異側(cè)，對(duì)于邊界邊，需要判斷一點(diǎn)是否位于該邊的內(nèi)側(cè)。利用向量叉積，可以簡(jiǎn)單地實(shí)現(xiàn)二者的判斷。如圖1所示，判斷任意點(diǎn)P與點(diǎn)P3位于邊 P1P2的同側(cè)或異側(cè)，只需將邊 P1P2、P1P3和P1P分別看作向量，首先根據(jù)式(1)和式(2)計(jì)算向量P1P2與P1P3以及向量P1P2與P1P的叉積，根據(jù)向量叉積的右手法則，可知向量P1P2與P1P3的叉積為負(fù)，向量P1P2與P1P的叉積為正；然后將兩叉積的結(jié)果相乘，所得數(shù)值如果小于0，則點(diǎn)P位于點(diǎn)P3的異側(cè)，大于0則點(diǎn)P位于點(diǎn)P3的同側(cè)，若等于0，則點(diǎn)P位于邊P1P2所在的直線上。對(duì)于邊界邊的擴(kuò)展，同樣可以利用該方法判斷點(diǎn)線位置，如果構(gòu)造邊界時(shí)是順時(shí)針?lè)较蜻B接邊界，則取點(diǎn)時(shí)只取位于邊界右側(cè)的點(diǎn)。在圖1中假設(shè)邊P1P2是離散數(shù)據(jù)集的一段邊界，取其左側(cè)的點(diǎn)，假設(shè)為P3，根據(jù)式(1)計(jì)算向量P1P2與P1P3的叉積，為負(fù)值就可以滿足要求。

1.3 最佳三角形的構(gòu)建

盡管對(duì)不均勻的離散數(shù)據(jù)進(jìn)行網(wǎng)格化只能生成不規(guī)則的三角網(wǎng)，但是在實(shí)際應(yīng)用中仍希望網(wǎng)格單元盡量趨近于飽滿，即每個(gè)三角形盡量接近正三角形。網(wǎng)格質(zhì)量的好壞將直接影響后續(xù)計(jì)算的精度。三角形外接圓半徑與內(nèi)切圓半徑的比值稱為三角形的橫縱比(Aspect radio)。用三角形的橫縱比來(lái)衡量三角形的質(zhì)量，具有大的橫縱比的三角形至少有一個(gè)小角且三個(gè)頂點(diǎn)接近共線。如圖2所示，有兩種質(zhì)量差的三角形單元，沒(méi)有短邊的Blade和有一個(gè)短邊的Dagger。

圖2 兩種畸形三角形單元Fig. 2 Two kinds of deformed triangle element

在進(jìn)行網(wǎng)格化時(shí)，利用距離和最小準(zhǔn)則能保證離散點(diǎn)集中最鄰近的三點(diǎn)構(gòu)成三角形，避免生成畸形三角形單元及三角形的交叉。第三點(diǎn)的確定，是通過(guò)使用三角形邊角關(guān)系的余弦定理來(lái)進(jìn)行的：三角形中任意兩邊的平方和減去此兩邊與其夾角的余弦值等于第三邊的平方，即：由式(3)可得

對(duì)于固定的 c，C 角越大，cosC 值越小，?2ab·(1+cosC)越大，a+b也就越小。如圖1所示，如果某點(diǎn)P的角C在每個(gè)滿足條件的點(diǎn)與待擴(kuò)展邊所構(gòu)成的三角形中是最大的，那么點(diǎn)P距離待擴(kuò)展邊兩端點(diǎn)P1、P2的距離之和最小，該三點(diǎn)是離散點(diǎn)集中最鄰近的三點(diǎn)，選該點(diǎn)P為第三點(diǎn)滿足最小距離和準(zhǔn)則，所生成的三角形符合最佳三角形條件。

2 三角網(wǎng)逆生長(zhǎng)網(wǎng)格化方法實(shí)現(xiàn)

本算法的實(shí)現(xiàn)采用C Sharp語(yǔ)言，算法流程圖如圖3所示，生長(zhǎng)示意圖如圖4所示。

圖3 程序流程圖Fig. 3 Program flow chart

具體實(shí)現(xiàn)步驟：

Step 1: 如圖 4(a)所示，對(duì)于給定的點(diǎn)集P= ｛ p｝n，首先按次序提取其邊界點(diǎn)并加入邊界點(diǎn)

圖 4 三角網(wǎng)逆生長(zhǎng)示意圖: (a) 離散點(diǎn)集; (b) 邊界提取;(c) 邊界邊生長(zhǎng); (d) 內(nèi)部邊生長(zhǎng)Fig. 4 Triangular mesh inverse growth sketch: (a) Discrete points set; (b) Boundary extraction; (c) Boundary growth;(d) Inner boundary growth

Step 2: 如圖4(c)所示，對(duì)邊界執(zhí)行生長(zhǎng)算法，判斷邊界邊的三角形數(shù)量，若等于0，對(duì)該邊進(jìn)行擴(kuò)展，查找其最佳第三點(diǎn)，構(gòu)造新邊，判斷該邊是否存在，若不存在，加入內(nèi)部邊集 E = ｛ ei｝k。2 i=0

Step 3: 如圖4(d)所示，對(duì)內(nèi)部邊執(zhí)行生長(zhǎng)算法，判斷內(nèi)部邊的三角形數(shù)量，若等于 1，對(duì)該邊進(jìn)行擴(kuò)展，查找其最佳第三點(diǎn)，構(gòu)造新邊，判斷該邊是否存在，若不存在，加入內(nèi)部邊集 E = ｛ e2i ｝ik=0。

3 實(shí)例

應(yīng)用某測(cè)區(qū)的EH4電磁測(cè)深數(shù)據(jù)，分別用傳統(tǒng)三角網(wǎng)生長(zhǎng)法及三角網(wǎng)逆生長(zhǎng)法進(jìn)行處理，得到如圖 5和6所示實(shí)驗(yàn)結(jié)果。由圖5可見(jiàn)，圖5(a)和圖5(b)所示分別為應(yīng)用傳統(tǒng)三角網(wǎng)生長(zhǎng)法和三角網(wǎng)逆生長(zhǎng)法得到的不規(guī)則三角網(wǎng)，圖 5(a)中虛線部分為未被正確擬合的邊界。由圖6可見(jiàn)，圖6(a)和圖6(b)所示分別為應(yīng)用兩個(gè)不規(guī)則三角網(wǎng)進(jìn)行可視化分析得到的等值線圖，圖 6(a)中，①處的等值線在地表處并不存在，是虛假等值線，②處的等值線超出數(shù)據(jù)區(qū)域，不合理。

圖5 三角網(wǎng)格化結(jié)果對(duì)比: (a) 生長(zhǎng)法三角網(wǎng)格化結(jié)果; (b)逆生長(zhǎng)法三角網(wǎng)格化結(jié)果Fig. 5 Comparison of triangulation results: (a) Triangulation result of growth algorithm; (b) Triangulation result of inverse growth algorithm

圖6 可視化結(jié)果對(duì)比: (a) 生長(zhǎng)法可視化結(jié)果; (b) 逆生長(zhǎng)法可視化結(jié)果Fig. 6 Comparison of visualization results: (a) Visualization result of growth algorithm; (b) Visualization result of inverse growth algorithm

通過(guò)對(duì)比分析可知，三角網(wǎng)逆生長(zhǎng)網(wǎng)格化方法得到的不規(guī)則三角網(wǎng)邊界擬合更準(zhǔn)確，效果更好，應(yīng)用于可視化分析所繪制的等值線均位于數(shù)據(jù)范圍以內(nèi)，在地表以上沒(méi)有形成虛假的等值線，深部也未存在不合理的等值線。

4 結(jié)論

1) 三角網(wǎng)逆生長(zhǎng)網(wǎng)格化方法思路簡(jiǎn)單，算法穩(wěn)定，生成的三角網(wǎng)符合Delaunay準(zhǔn)則，質(zhì)量最優(yōu)且唯一，能簡(jiǎn)單而精確地?cái)M合物探數(shù)據(jù)的復(fù)雜邊界。

2) 相比于常用的矩形網(wǎng)格化方法，三角網(wǎng)逆生長(zhǎng)網(wǎng)格化方法得到的網(wǎng)格節(jié)點(diǎn)值全部是原始數(shù)據(jù)，避免了大量的數(shù)值計(jì)算以及由此引起的誤差；相比于傳統(tǒng)的三角網(wǎng)生長(zhǎng)算法，恰當(dāng)?shù)貞?yīng)用了物探數(shù)據(jù)的特點(diǎn)進(jìn)行邊界處理，避免了復(fù)雜的邊界計(jì)算，應(yīng)用于物探數(shù)據(jù)處理及資料解釋符合實(shí)際情況，以此為基礎(chǔ)的可視化計(jì)算得到的等值線圖不產(chǎn)生非法等值線。

[1] 李海生. Delaunay三角剖分理論及可視化應(yīng)用研究[M].哈爾濱: 哈爾濱工業(yè)大學(xué)出版社, 2010: 4?7.

LI Hai-sheng. Study of Delaunay triangulation theory and visualization applications [M]. Harbin: Harbin Institute of Technology Press, 2010: 4?7.

[2] SHAMOS M I, HOEY D. Closest-point problems[C]//IEEE Computer Society. Proceedings of the 16th Annual Symposium on the Foundations of Computer Science. New York: ACM Press,1975: 151?162.

[3] LAWSON C L. Software for C1 surface interpolation [C]// RICE J R. Mathematical software Ⅲ. New York: Academic Press,1977: 161?194.

[4] GREEN P J, SIBSON R. Computing dirichlet tessellations in the plane [J]. The Computer Journal, 1978, 21(2): 168?173.

[5] CHEW L P. Constrained Delaunay triangulations [C]//SIGGRAPH. Proceedings of the Third Annual Symposium on Computational Geometry. New York: ACM Press, 1987: 215?222.

[6] 楊 欽, 徐永安, 陳其明, 譚建榮. 任意平面域離散點(diǎn)集的三角化方法[J]. 軟件學(xué)報(bào), 1998, 9(4): 241?245.

YANG Qin, XU Yong-an, CHEN Qi-ming, TAN Jian-rong.Triangulation algorithm of scattered data on arbitrary planar domain [J]. Journal of Software, 1998, 9(4): 241?245.

[7] 劉少華, 程朋根, 趙寶貴. 約束數(shù)據(jù)域的 Delaunay三角剖分算法研究及應(yīng)用[J]. 計(jì)算機(jī)應(yīng)用研究, 2006(3): 26?28.

LIU Shao-hua, CHENG Peng-gen, ZHAO Bao-gui. A study on algorithm of Delaunay triangulation for the constrained data set and application [J]. Application Research of Computers, 2006(3):26?28.

[8] 高曉楓. 2D-Delaunay三角網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)與遍歷[J].天津理工大學(xué)學(xué)報(bào), 2006, 22(2): 66?69.

GAO Xiao-feng. Data structure and traverse of 2D-Delaunay triangulation [J]. Journal of Tianjin University of Technology,2006, 22(2): 66?69.

[9] 周佳文, 薛之昕, 萬(wàn) 施. 三角剖分綜述[J]. 計(jì)算機(jī)與現(xiàn)代化,2010(7): 75?79.

ZHOU Jia-wen, XUE Zhi-xin, WAN Shi. Survey of triangulation methods [J]. Computer and Modernization, 2010(7): 75?79.

[10] WANNAMAKER P E, SYODT J A, RIJO L. Two dimensional topographic responses in magnetollurics modeled using finite elements [J]. Geophysics, 1986, 51(11): 2131?2144.

[11] 阮百堯, 村上裕, 徐世浙. 激發(fā)極化數(shù)據(jù)的最小二乘二維反演方法[J]. 地球科學(xué): 中國(guó)地質(zhì)大學(xué)學(xué)報(bào), 1999, 24(6):619?624.

RUAN Bai-yao, YUTAKA Murakami, XU Shi-zhe. Least square 2-D inversion method for induced polarization data [J]. Earth Science: Journal of China University of Geosciences, 1999,24(6): 619?624.

[12] 李桐林, 齊守奎, 趙海珍. 有限元金屬礦地震波場(chǎng)模擬中三角網(wǎng)的自動(dòng)剖分與模型編輯[J]. 鈾礦地質(zhì), 2000, 16(1): 42?45.

LI Tong-lin, QI Shou-kui, ZHAO Hai-zhen. Automatic area dividing and editing techniques for finite element ore deposit seismic modeling [J]. Uranium Geology, 2000, 16(1): 42?45.

[13] 湯井田, 任政勇, 化希瑞. 任意地球物理模型的三角形和四面體有限單元剖分[J]. 地球物理學(xué)進(jìn)展, 2006, 21(4):1272?1280.TANG Jing-tian, REN Zheng-yong, HUA Xi-rui. Triangle and tetrahedral finite element meshing fromarbitrary geophysical model data [J]. Progress in Geophysics, 2006, 21(4): 1272?1280.

[14] KERRY K, CHESTER W. Adaptive finite-element modeling using unstructured grids: The 2D magnetotelluric example [J].Geophysics, 2006, 71(6): 291?299.

[15] TANG Jing-tian, WANG Fei-yan, REN Zheng-yong. 2.5-D DC resistivity by adaptive finite-element method with unstructured triangulation [J].Chinese Journal of Geophysics, 2010, 53(3):708?716.

[16] 劉兆平, 楊 進(jìn), 武 煒. 地球物理數(shù)據(jù)網(wǎng)格化方法的選取[J]. 物探與化探, 2010, 34(1): 93?97.

LIU Zhao-ping, YANG Jin, WU Wei. The Choice of gridding methods for geophysical data [J]. Geophysical and Geochemical Exploration, 2010, 34(1): 93?97.