基于TIN的地質界面三維形態分析方法與應用

毛先成 ，趙瑩 ，唐艷華 ，彭正林 ，陳進 ，鄧浩

(1. 中南大學 有色金屬成礦預測教育部重點實驗室，湖南 長沙，410083；2. 中南大學 地球科學與信息物理學院，湖南 長沙，410083)

界面是指兩相間接觸的交界部分，兩物態間的轉換面[1]。地質體是地質空間中不連續的空間實體，包括沉積成因的巖層、侵入成因的巖漿巖體以及受力變形的構造等[2]。宏觀的地質界面可以理解為地質體的分界面或間斷面，或者稱為面狀地質體[2-3]。在地質學上，界面兩側的巖石具有不同的巖性類型、不同的礦物組成、不同的化學成分、不同的結構構造。地質界面是地球內部能量集中、匯集、傳遞、轉化和釋放的地帶，是地球物質變化、活化、遷移和沉淀富集的部位，也是成礦作用發生的場所。自然界中，地質界面成礦現象普遍存在，礦產資源往往在各種形式的地質界面處超常富集。運用地質界面控礦原理進行找礦預測，關鍵是準確識別并定位地質界面，研究其性質并理解其控礦作用，并融合其他成礦理論進行綜合分析，同時注意地質界面控礦在地質異常解釋中的運用[1]。對地質體數學特征的研究有助于查明地質體成因，追溯地質體形成過程和機理，并有助于發現一些隱蔽規律性，從而更好地解決認識、預測和評價等理論和實際問題[4]。所以，采用空間分析方法對地質體和地質界面三維幾何形態等數學特征進行提取和定量分析，對于發掘有利成礦信息具有重要意義。現有的研究多關注于地質體和地質界面模型的建立和表達[5-6]等，但是，在利用三維空間形態分析提取地質控礦因素方面，見到的研究文獻還比較少。針對控礦地質因素提取的三維形態分析問題，本文作者曾提出采用基于柵格模型的地質體三維形態分析方法進行控礦地質因素的分析和提取[7]，該方法是地質三維建模與數學形態學[8-10]相結合形成的一種基于柵格的三維空間分析方法。基于柵格模型的地質體三維形態分析方法使用體素模型(三維柵格模型)表達地質對象空間數據，具有位置隱含、結構簡單、操作方便的優點，可描述地質現象(地質場)連續變化的特性，特別適應于超覆、彎曲(褶皺)、錯斷等復雜形態的地質體。但是，研究結果分析表明，該方法存在如下缺點：柵格數據模型空間數據量大、空間精度相對較低；在進行邊界表示時不能準確表達地質體的邊界[11]，且生成速度、存儲空間和表達精度難以平衡；數學形態學球形結構元素半徑的選取在量化方面需要改進；不適合對薄片狀地質體或存在狹窄部分的地質體進行趨勢形態分析，當結構元素半徑較大時，會出現斷裂的錯誤。因此，還需提出一種更合適地質界面表達與三維形態分析的方法。基于TIN模型的表達方法可有效地模擬地質界面[12]。不規則三角網(Triangulated Irregular Network, TIN)模型可以根據表面的復雜程度變化三角形面片的大小和數量，消除可視化時的數據冗余，并且具有三維面表示模型的易于采集與構建、表達精度高、便于顯示和更新的優點，適于進行地質體輪廓的幾何建模，建模理論成熟，方法易于實現。基于上述分析，本文作者根據地質界面形態分析的實際需要，將利用TIN模型表達和分析地質界面精度高等方面的優勢，與三維地質建模、三維空間分析等方法結合起來，提出了一種精度更高的地質界面形態分析方法—基于TIN的地質界面三維形態分析方法，可以提取地質界面的幾何形態參數、分析地質界面的多級形態趨勢與起伏等，并建立控礦地質因素場模型，實現控礦地質因素的定量表達與可視化。

1 地質界面三維形態分析方法

1.1 TIN模型與地質界面三維形態分析

TIN模型將空間無重復點的散亂點按某種規則(如Delaunay規則)進行三角剖分而形成連續卻不重疊的不規則三角面片網，以此來模擬三維物體的表面。它具有固定的結構和簡單的處理過程，除了表達各種元素的空間位置以外，同時還可以較好地表示三角形之間的拓撲關系，能夠保持測量數據的原始性，且在表達地質界面的不連續方面具有優勢，因而目前主要采用TIN模型來模擬地質界面。利用TIN模型來建模時，在模型簡單的區域，三角形較少，而在模型復雜的區域，三角形較多。因此，TIN模型能較好地顧及模型特征，逼真表示模型的高低起伏變化，同時能夠克服模型平坦區域的數據冗余。TIN的基本單元三角形的幾何形狀直接決定著TIN的應用質量。TIN模型由數據結構、三角劃分準則和三角剖分算法等三部分組成。從結構上講，TIN是一種矢量數據結構，對三角形的點、線、面的拓撲描述不同而形成不同的數據結構。近年來，借助于飛速發展的計算機軟硬件技術，在TIN的快速建立、壓縮存儲、數據組織等方面均取得了突破性的進展[13-16]。TIN在地形分析[17]、空間統計分析[18]、幾何分析[19]等空間分析中發揮著重要作用。

地質界面三維形態分析，以控礦因素定量提取為目的，以地質界面三維TIN模型為基礎，利用空間幾何原理計算地質界面的一般幾何形態參數(坡度、夾角等)和距離場，利用空間插值技術和趨勢剩余分析技術，提取出地質界面的局部形態特征和總體趨勢，并以柵格模型來表達控礦地質因素場，實現各種控礦地質因素的定量模擬，為進一步建立礦體的立體定量預測模型提供基礎。本研究提出并實現的地質界面三維形態分析方法包括：(1) 地質界面的幾何形態參數提取；(2) 地質界面的形態趨勢與起伏分析。

1.2 地質界面的幾何形態參數提取

1.2.1 坡度提取

坡度是地質界面的基本屬性。由于坡度的計算都是基于地質界面的TIN 模型，因此對于給定點的坡度計算實質上是求該點所對應的三角面片的坡度。具體過程如下：

TIN模型上任意三角形 ABC可以用如下平面方程表示：

平面上坡度處處相等，可以用下式計算該三角形的坡度：

設三角形ABC的三個頂點的坐標為A(x1，y1，z1)，B(x2，y2，z2)，C(x3，y3，z3)，用向量法可以得出三角形所在平面的方程：

結合式(1)和(3)可以得出：

最后將式(4)和(5)代入式(2)就可以求出α，即對應于每個三角面片的坡度信息。通過距離場算法能獲得地質空間中的立體單元到最近的三角面片的標示ID，該三角面片的坡度信息就是該立體單元所要賦值的坡度。

1.2.2 地質界面間夾角提取

求TIN模型中2個地層界面交線上的某點的夾角就是求經過該點的2個三角面片之間的夾角。但是，在TIN模型中，確定地質界面之間相交的三角面片難度較大，需要對2個界面上的三角面片逐一進行幾何求交計算來獲得實際相交的三角面片，尋找相交的三角面片的計算量與構成2個界面的三角面片的數量的乘積呈正相關，計算量大、時間效率低，所以本文提出了一種簡單高效的方法來求取地質界面之間相交的三角面片，從而計算兩地層界面之間的夾角。該方法首先將地質界面柵格化，獲取三角面片的穿過的塊體單元的網格索引，通過網格索引直接定位有可能實際相交的三角面片，然后采用向量法直接計算這些三角面片之間的夾角，最后通過賦值處理得到地質界面之間的夾角。

(1) 地質界面柵格化。地質界面是以TIN模型表達的，其柵格化的實質就是將地質界面TIN穿過的塊體單元提取出來，形成地質界面的塊體模型。具體步驟如下：

① 將三角面片的 3個頂點向坐標軸 XY平面投影，得到三角形為ABC(如圖1所示)。

② 獲取最小外包矩形。

③ 采用計算機圖形學中的 Liang_Barsky參數化線段裁剪多邊形算法[20]，計算三角形的三條邊穿過的平面網格。

④ 采用掃描線算法計算 XY平面上三角形 ABC內部的網格索引(ix，jy)。

⑤ 重復這些步驟分別求得在XZ平面和YZ平面上相應的內部網格索引，然后將三個投影面上的網格索引的分量求交，最后求得三角面片通過的塊體單元的索引，并且記錄該三角形的ID。

⑥ 通過對所有三角面片循環計算，可以求出地質界面通過的塊體單元，最后生成地質界面塊體模型。

圖1 TIN 模型中某個三角面片在XY平面上投影圖Fig. 1 XY plane projection of a triangular of TIN model

(2) 夾角的計算。首先將兩地層界面的塊體模型進行求交，其實質就是計算塊體模型中重合的塊體單元的索引值，同時記錄穿過重合塊體單元的三角面片的標識 ID值，然后計算穿過同一塊體單元的三角面片之間的夾角。其中，需要對一種情況進行特殊處理：當相交處的塊體單元為三角面片的頂點時，地質界面的TIN模型就會有多個三角面片同時穿過相應的塊體單元，這種情況下，對2個地質界面TIN模型中穿過相同塊體單元的三角面片兩兩求夾角，最后將這些夾角的平均值作為最后結果。

2個三角形的夾角計算經常采用向量法，首先規定三角形的法向量向上，然后求2個三角形的法向量，法向量的夾角即為2個平面的夾角。

設 2個三角形的法向量分為 V1(x1，y1，z1)和V2(x2，y2，z2)，兩法向量的夾角為α，則：

(3) 夾角的賦值。通過上面步驟，獲得了兩地層界面相交處的塊體集合，并且該集合中的塊體包含了兩者的夾角信息。因此，對于整個地質空間中的立體單元就通過求取與上述塊體集合中的塊體的最近距離來獲取相應的夾角信息。求得的塊體集合中最短距離所對應的塊體的夾角信息就是該立體單元所需要賦值的夾角。

1.3 地質界面的形態趨勢與起伏分析

在礦床地質空間中，相關地質體的主體形態將對周圍成礦起到一定因素的影響。因此，對地質體形態的分析將是控礦因素定量提取的重要步驟。

地質體形態可以通過地質界面的波狀起伏來描述。借鑒傳統的曲面的趨勢形態分析方法，結合地質界面的實際情況，本文采用地質界面的原始TIN模型，利用距離平方反比法對地質界面進行趨勢形態分析。地質界面TIN模型的趨勢形態分析實際上就是利用估測點周圍一定距離之內的樣品點對估測點的z重新賦值，然后利用計算得到的新的z生成地質界面的趨勢面。具體過程如下。

(1) 確定搜索距離 d，只有與估測點距離小于 d的樣品點才能參與計算，其中距離是指兩點在 XY平面上的相對距離，而不是兩點的絕對距離。

(2) 利用距離平方反比法重新計算樣品點的所有z，公式如下：

(3) 利用計算得到的新的 z，建立新的地質界面的TIN模型，形成地質界面的趨勢面。

地質體起伏形態可以通過地質界面的波狀起伏來描述。這種波狀起伏可以通過地質界面相對于其趨勢平面的起伏幅度來表達。這種起伏幅度又叫波幅，反映地質界面的彎曲程度和在空間中的波狀起伏幅度變化。

這種波狀起伏的建模可以通過前文所述的趨勢形態分析來實現，即將地質界面標高z的觀測值的變化分解為 2個部分：(1) 觀測值總體變化重建的一個曲面，即趨勢面T(x，y)；(2) 觀測值的局部變化，稱謂為信號，即剩余趨勢面 R(x，y)；則剩余趨勢面R(x，y)：

對于地質曲面上的一個點 z(x0，y0)，可通過R(x，y)絕對值來描述該點所處的局部的形態起伏程度，通過R(x，y)的符號來描述該點所處的局部狀態，有：

對于R(x，y)所表示的局部起伏，可視為不同級別起伏的疊加，有：

其中：R1(x，y)，R2(x，y)，…，Rn(x，y)表示n級局部起伏變化，各級局部變化的篩分可用趨勢剩余分析方法。

利用上述方法可實現形態起伏的分級提取。因為插值搜索半徑d決定可濾除的波形的大小，所以通過改變d可提取出不同級別的起伏。

地質界面的剩余趨勢面通過波幅很好地反映了地質界面在空間中的起伏變化程度，為了定量描述地質界面上任何一點的波幅，需要計算地質界面與趨勢面上具有相同x和y的點的z坐標的差，實質上是通過x，y計算地質界面和趨勢面的z。計算過程如下：

地質界面TIN模型中任意三角形ABC(A(x1，y1，z1)，B(x2，y2，z2)，C(x3，y3，z3))，三角形 ABC 所在平面中的任意點M(x，y，z)，平面方程同式(3)。

可得：

其中：a=(x2-x1)(y3-y1)-(y2-y1)(x3-x1)，b=(z2-z1)×(y3-y1)(x-x3)+(y-y1)(x2-x1)(z3-z1)-(x-x1)(y2-y1)(z3-z1)-(y-y1)(z2-z1)(x3-x1)。

但是，當三角形ABC垂直于XY平面時，相同的x和y會有無數個z對應，對于此種情況，將z確定為最小值與最大值和的一半，即：

2 地質界面三維形態分析應用實例

實例：福建省尤溪縣丁家山鉛鋅礦床的不整合面。

丁家山鉛鋅礦床是屬于梅仙鉛鋅礦田的一個大型礦床。其編號為Ⅰ，Ⅱ和Ⅲ的三層礦體分布于龍北溪組上段(Pt2-3l3)綠片巖地層，產狀平緩，礦體產狀與含礦地層產狀協調一致，礦床空間分布特征與該含礦巖性段總體延伸特征保持一致。在侏羅系長林組(J3c)地層中產出鉛鋅礦體(Ⅳ號礦體)，呈長條狀或透鏡狀，礦體傾向與不整合面產狀基本一致[21]。

龍北溪組上段綠片巖系為主要賦礦層位，燕山期石英斑巖、花崗斑巖脈在形成過程中礦區地層和含礦巖系是主要的控礦地質因素。某些鉛鋅礦體賦存于巖層構造裂隙中，礦體受斷裂構造控制，產狀與斷裂構造產狀一致，呈脈狀，受裂隙構造控制的礦體的形成和分布可能與后期的熱液疊加改造有關[21]。

不整合面在礦區內分布廣泛，是控制礦體分布的主要構造。產于龍北溪組上段(Pt2-3l3)綠片巖地層中的礦體同時還受不整合構造控制。當不整合面較陡，且與下伏綠片巖地層層面呈較大角度斜交時，則在不整合面兩側一定范圍內均有利于成礦。一方面，在不整合面構造兩側，礦體的產狀受不整合面產狀的影響，且礦體的厚度和品位有增大變富的趨勢；另一方面，在不整合面構造中，還產出有充填型的脈狀礦體。

從定性的角度，可以發現丁家山鉛鋅礦床不整合面的空間形態與成礦存在一定聯系，然而對這種關系的認識是有限的，而且存在主觀性，于是采用三維形態分析來對其進行定量化研究成為必然。又由于丁家山鉛鋅礦床不整合面屬于正常產狀，產狀平緩，因此采用本方法進行形態分析。首先，采用Datamine軟件建立地質界面三維TIN數據模型，并且建立礦化空間的柵格模型，然后利用研究開發的三維形態提取算法和程序，對不整合面地層界面進行三維形態分析，如對不整合面的形態趨勢提取以及形態起伏分級提取,不整合面坡度提取，不整合面的夾角提取，不整合面距離場因素指標提取等，最后利用礦化空間柵格模型來表達各種控礦地質因素場。

圖2所示為不整合面形態起伏分級提取。不整合面趨勢-起伏因素分析主要是揭示不整合面起伏對不整合面周圍地質空間的控礦作用影響。經過反復試驗，選擇100 m和200 m作為一級濾波和二級濾波的插值搜索范圍半徑，對不整合面原始TIN模型進行一級形態濾波和二級形態濾波，對應的得到不整合面的一級趨勢與一級起伏(圖 2(a))、二級趨勢與二級起伏(圖2(b))。

waU和 wbU分別代表不整合面一級起伏和二級起伏的起伏程度，用單元到不整合面的最近距離處的不整合面趨勢部分的歐式距離 Einner或 Eouter來度量。顯然，若單元到不整合面的最近距離處屬于不整合面相應級別起伏的外凸部分，則waU和wbU大于0；若單元到不整合面的最近距離處屬于不整合面相應級別起伏的內凹部分，則waU和wbU小于0。計算結果如圖3所示。

陡傾不整合構造是層控型鉛鋅礦床后期改造最有利的構造，陡傾不整合構造旁側礦體均強烈顯示變厚變富的趨勢，因此不整合面坡度因素gU的提取具有重要意義。計算結果的柵格模型表達如圖4所示。

不整合面與下伏綠片巖地層層面呈較大角度斜交的這種對成礦有利的情況可以利用不整合面夾角因素aU_S來表達。不整合面與多個地層界面有相交，定義單元的不整合面夾角指標值取值為單元到不整合面的最近距離處對應的最近的夾角。計算結果如圖5所示。陡傾不整合構造旁側礦體均強烈顯示變厚變富的趨勢，并且產于火山巖蓋層及其他層位的“充填型”礦體，受裂隙構造控制，這種構造因素可以用不整合面的距離場模型描述。從丁家山鉛鋅礦床不整合面三維TIN模型獲取不整合面的線框模型，以單元到不整合面的最小距離作為距離場值來表達，即可得到不整合面距離場因素指標dU，如圖6所示。為區分不整合面上下單元的距離場值的差異，將位于不整合面分布范圍之上的單元場值置為正(正距離場)，位于不整合面分布范圍之下的單元場值為負(負距離場)。

圖2 不整合面形態起伏分級提取Fig. 2 Hierarchical extraction of undulate shapes of unconformity surface

圖3 不整合面形態起伏因素(waU和wbU)分布柵格模型(-150～-200 m標高范圍)Fig. 3 Raster model of shape factors of unconformity surface (waU and wbU)

圖4 不整合面坡度因素(gU)分布柵格模型(-150～-200 m標高范圍)Fig. 4 Raster model of slope factors of unconformity surface (gU)

圖5 地層界面與不整合面的夾角因素(aU_S)分布柵格模型(-150～-200 m標高范圍)Fig. 5 Raster model of angle factors of stratum interface and unconformity surface (aU_S)

圖6 不整合面距離因素(dU)分布柵格模型(-150～-200 m標高范圍)Fig. 6 Raster model of distance factors of unconformity surface (dU)

3 結論

(1) 結合三維地質建模、空間插值、地質場和空間分析等理論，提出了基于TIN的地質界面三維形態分析方法，通過該方法，可以快速有效的進行地質界面坡度和夾角等幾何形態參數提取、地質界面距離場分析、地質界面形態趨勢與起伏分析等，克服了柵格模型和數學形態學用于地質界面形態分析的缺點，是一種精度更高、用途更廣泛的新方法。

(2) 以福建省尤溪縣丁家山鉛鋅礦床為例，利用提出的方法對地層界面和不整合面等地質界面進行形態分析，提取了與鉛鋅成礦的各種地質形態因素，可以合理地定量化描述地質現象，進而在成礦預測中發揮作用，如陡傾不整合構造是成礦有利因素，而通過形態分析提取的不整合面坡度因素gU，可以便于定位到坡度陡傾的位置，或者便于建立坡度陡傾與成礦的非線性函數關系。

(3) 礦產資源的三維立體定量預測及三維評價系統是危機礦山可接替資源評價與找礦創新體系建立過程中亟待探索和實踐的新理論、新方法。采用地質界面三維形態分析方法，以地質界面三維TIN模型的建立和可視化為基礎，提取出地質界面空間形態指標，從而得到控礦地質因素場模型，進而用于控礦地質因素的定量分析，在隱伏礦體三維可視化立體定量預測中具有重要意義。

[1] 張善明, 呂新彪, 鄧國祥, 等. 地質界面控礦原理及其運用要點[J]. 地質科技情報, 2009, 28(6): 51-56.ZHANG Shanming, Lü Xinbiao, DENG Guoxiang, et al.Ore-controlling mechanism of geological interface and the key of application[J]. Geological Science and Technology Information, 2009, 28(6): 51-56.

[2] 朱志澄, 宋鴻林. 構造地質學[M]. 武漢: 中國地質大學出版社, 1990: 5-14.ZHU Zhicheng, SONG Honglin. Structural geology[M]. Wuhan:China University of Geosciences Press, 1990: 5-14.

[3] 鄧浩. 面向隱伏礦體預測的三維地質建模與空間分析若干技術研究[D]. 長沙: 中南大學地球科學與信息物理學院, 2008:1-17.DENG Hao. Study on 3D geological modeling and spatial analysis for prediction of buried ore bodies[D]. Changsha:Central South University. School of Geosciences and Info-Physics, 2008: 1-17.

[4] 趙鵬大. 試論地質體的數學特征[J]. 地球科學, 1982(1):145-155.ZHAO Pengda. On the mathematical characteristics of geological bodies[J]. Earth Science-Journal of China University of Geoscience, 1982(1): 145-155.

[5] 周訪濱, 劉學軍. 基于柵格DEM自動劃分微觀地貌形態的研究[J]. 武漢理工大學學報: 信息與管理工程版, 2008, 30(2):172-175.ZHOU Fangbin, LIU Xuejun. Research on the automated classification of micro landform based on grid DEM[J]. Journal of Wuhan University of Technology: Information &Management Engineering, 2008, 30(2): 172-175.

[6] 姜棟, 趙文吉, 朱紅春, 等. DEM地形信息提取對比研究: 以坡度為例[J]. 測繪科學, 2008, 33(5): 177-179.JIANG Dong, ZHAO Wenji, ZHU Hongchun, et al. Comparison of landform information extracted from DEMs: A case study of slope[J]. Science of Surveying and Mapping, 2008, 33(5):177-179.

[7] 毛先成, 唐艷華, 鄧浩. 地質體的三維形態分析方法與應用[J]. 中南大學學報: 自然科學版, 2012, 43(2): 588-595.MAO Xiancheng, TANG Yanhua, DENG Hao.Three-dimensional morphological analysis method for geologic bodies and its application[J]. Journal of Central South University:Science and Technology, 2012, 43(2): 588-595.

[8] Serra J. Introduction to mathematical morphology[J]. Computer Vision, Graphics, and Image Processing, 1986, 35: 283-305.

[9] Serra J. Image analysis and mathematical morphology[M]. San Diego: Academic, 1982: 15-58.

[10] Jahjah M, Ulivieri C, Santilli G. Automatic archaeological feature extraction from satellite VHR images by mathematical morphology functions[C]//International Astronautical Federation-59th International Astronautical Congress 2008. Paris:IAF, 2008: 2751-2755.

[11] 吳立新, 史文中, Christopher G. 3DGIS與3DGMS中的空間構模技術[J]. 地理與地理信息科學, 2003, 19(1): 5-11.WU Lixin, SHI Wenzhong, Christopher G. Spatial modeling technologies for 3D GIS and 3D GMS[J]. Geography and Geo-Information Science, 2003, 19(1): 5-11.

[12] 彭儀普. 地形三維可視化及其實時繪制技術研究[D]. 成都:西南交通大學土木工程學院, 2002: 1-19.PENG Yipu. Studying on three-dimension visualization and real-time rendering of terrain[D]. Chengdu: Southwest Jiaotong University. School of Civil Engineering, 2002: 1-19.

[13] WU Hangbin, LIU Chun. Compression of discrete data based on the triangulate irregular network[C]//Proceedings of the SPIE-The International Society for Optical Engineering.Bellingham: SPIE, 2006: 6420.

[14] LU Linlin, GUO Huadong. Visualization of a digital elevation model[J]. Data Science Journal, 2007, 6: 481-484.

[15] Kothuri R, Ravada S, Fisher E. Triangulated irregular network:US, 12069089[P]. 2010-08-10.

[16] ZHU Qing, ZHANG Yeting, LI Fengchun. Three-dimensional TIN algorithm for digital terrain modeling[J]. Geo-spatial Information Science, 2008, 6(11): 79-85.

[17] YU Ming, CHEN Xiaoyu, AI Tinghua. Application research of terrain based on DEM and data mining[C]//CHEN Jingming, PU Yingxia. Geoinformatics 2007: Geosputial information Science.Bellingham, Wash: SPIE, 2007: 67531w.

[18] 鐘登華, 郭享, 李明超, 等. 基于三維地質模型的地下洞室參數化設計與方案優選[J]. 天津大學學報, 2007, 40(5): 519-524.ZHONG Denghua, GUO Xiang, LI Mingchao, et al. Parametric design and schemes optimization for underground structure based on 3D geological model[J]. Journal of Tianjin University,2007, 40(5): 519-524.

[19] 許麗敏, 薛安. 基于Delaunay三角網與Voronoi圖聯合提取等高線骨架的地形重建算法研究[J]. 北京大學學報自然科學版,2009, 45(4): 647-652.XU Limin, XUE An. Terrain reconstruction from contours by skeleton extraction using Delaunay triangulation and Voronoi diagram[J]. Acta Scientiarum Naturalium Universitatis Pekinensis, 2009, 45(4): 647-652.

[20] 孫家廣.計算機圖形學[M].北京:清華大學出版社,1998:166-215.SUN Jiaguang. Computer graphics[M]. Beijing: Tsinghua University Press, 1998: 166-215.

[21] 吳建設, 黃仁生. 福建省尤溪峰巖鉛鋅銀資源潛力及找礦方向探討[J]. 中國地質, 2001, 28(12): 13-18.WU Jianshe, HUANG Rensheng. Potentials of lead-zinc and silver resources and their prospecting direction at Fengyan,Youxi, Fujian[J]. Chinese Geology, 2001, 28(12): 13-18.