王 威，徐亞杏，汪振東，張 炯，高一鳴

武漢工程大學(xué)資源與土木工程學(xué)院，湖北武漢 430074

0 引言

三維地質(zhì)建模與可視化技術(shù)［1］是通過(guò)獲取原始地質(zhì)勘探數(shù)據(jù)如地質(zhì)圖、地形線(xiàn)、剖面、鉆孔等，解譯地質(zhì)對(duì)象，利用各種信息形成一個(gè)復(fù)雜三維模型的過(guò)程，從而實(shí)現(xiàn)可視化.我國(guó)對(duì)于研究簡(jiǎn)單地質(zhì)體的建模技術(shù)還是比較成熟的，但對(duì)于復(fù)雜地質(zhì)現(xiàn)象(如斷層、透鏡體、巖石破碎帶、地層倒轉(zhuǎn)等)研究得還不特別深入［2］.

目前，國(guó)內(nèi)使用的斷層建模方式主要有基于地層恢復(fù)的建模技術(shù)(整體法)、基于分區(qū)插值的建模技術(shù)(局部法)、斷層與地層的統(tǒng)一建模技術(shù)三類(lèi).本文基于統(tǒng)一建模技術(shù)方法，分析城市地質(zhì)斷層數(shù)據(jù)的特點(diǎn)，將斷層線(xiàn)等各種邊界線(xiàn)作為約束邊插入，重新進(jìn)行三角剖分，生成Tin(Triangulated Irregular Network)模型.當(dāng)勘察數(shù)據(jù)更新或地質(zhì)認(rèn)識(shí)改變時(shí)，只需重新調(diào)整斷層特征線(xiàn)，便可在調(diào)整部分生成新的斷層模型來(lái)進(jìn)行實(shí)時(shí)更新.通過(guò)實(shí)際應(yīng)用，相較于其它方法，此法更尊重原始地質(zhì)勘察數(shù)據(jù)，特別是在建造大型三維地質(zhì)體模型時(shí)，大大節(jié)省了時(shí)間，同時(shí)保證了模型的準(zhǔn)確性，使模型更加符合真實(shí)地質(zhì)情況.本文最后以武漢市漢陽(yáng)區(qū)為研究對(duì)象建立了三維地質(zhì)模型，分析其中某一小塊斷層模型的生成過(guò)程，驗(yàn)證了此構(gòu)模方法的有效性.

1 含斷層地質(zhì)體的生成原理

含斷層地質(zhì)體的生成主要采用多源數(shù)據(jù)耦合的建模方法(如圖1)，通過(guò)勘探數(shù)據(jù)，如收集的鉆孔資料、地表地質(zhì)圖、基巖地質(zhì)圖、地質(zhì)剖面(利用鉆孔人工生成剖面并通過(guò)實(shí)地考察檢驗(yàn)和地質(zhì)專(zhuān)家的分析來(lái)確定是否合理，如合理，即采用)等信息作為參考，利用帶約束Delaunay三角網(wǎng)建立斷層面、第四系地層面和基巖面等的 TIN模型，再通過(guò)手動(dòng)拼接輪廓線(xiàn)對(duì)某一地層的頂面與底面縫合，形成一個(gè)完整的地質(zhì)體模型.

2 帶約束邊的Delaunay三角剖分法

約束三角剖分就是利用約束條件，影響三角剖分的過(guò)程，或者在三角剖分之后根據(jù)這些約束條件去調(diào)整三角網(wǎng).在地質(zhì)曲面重構(gòu)的過(guò)程中地表模型的斷裂線(xiàn)、地質(zhì)剖面的斷層線(xiàn)和地層邊界、地質(zhì)圖上的地形線(xiàn)等都可作為重要的約束條件.此方法的基本思路是:首先，利用無(wú)約束的點(diǎn)集生成初始Delaunay三角網(wǎng)，再將較復(fù)雜的地質(zhì)中的特征線(xiàn)作為約束線(xiàn)依次插入到初始三角網(wǎng)中來(lái)重新調(diào)整初始三角網(wǎng)，使之符合Delaunay三角網(wǎng)的基本特性，最后就形成了帶約束的Delaunay三角網(wǎng).

圖1 城市三維建模流程圖Fig.1 City 3d modeling flowchart

2.1 無(wú)約束點(diǎn)集生成Delaunay三角網(wǎng)過(guò)程

設(shè)點(diǎn)集S里有n個(gè)離散點(diǎn).

①尋找點(diǎn)集周?chē)x質(zhì)心最近的點(diǎn)作為起始點(diǎn)，假設(shè)為點(diǎn)O;②計(jì)算出中心點(diǎn)和點(diǎn)集內(nèi)其余n-1個(gè)點(diǎn)的方向角、距離，升序排列其他所有點(diǎn)的傾斜度、方向角和距離的大小;③將O點(diǎn)和其他n-1個(gè)點(diǎn)連接成扇形依次進(jìn)行徑向掃描;再連接扇形的兩個(gè)角點(diǎn)，形成初始三角網(wǎng)(圖2-A);④從扇形邊的任一點(diǎn)開(kāi)始，按順時(shí)針或逆時(shí)針?lè)较虬紓?cè)連接，如圖2-B所示.假如以a為起點(diǎn)，規(guī)定按沿逆時(shí)針?lè)较蛑鹨徊檎?，a和b的下一個(gè)點(diǎn)分別是b和c，如果c點(diǎn)出現(xiàn)在線(xiàn)段ab方向的左側(cè)，則b點(diǎn)選為當(dāng)前點(diǎn)，從b點(diǎn)繼續(xù)向前查找;如果c點(diǎn)出現(xiàn)在線(xiàn)段ab的右側(cè)，就連接 ac，形成新三角形△abc(如圖2-B)，繼續(xù)向下搜索.同理，d點(diǎn)在ac的右側(cè)，連接ad，又生成新的三角形△acd，繼續(xù)查找d的下一個(gè)點(diǎn)e，發(fā)現(xiàn)點(diǎn)e在ad的左側(cè)，則d仍為當(dāng)前點(diǎn)，繼續(xù)以d點(diǎn)向前查找，直到將凹多邊形全變成凸多邊形為止(如圖2-C).當(dāng)徑向三角形極點(diǎn)被連接并第一次掃描徑向三角形外邊時(shí)，需要用一個(gè)數(shù)組記錄外包邊上的相鄰邊，則下一次掃描時(shí)就可以直接定位到下條邊［3］.⑤以邊界內(nèi)某一三角形為起始三角形，開(kāi)始對(duì)凸四邊形的對(duì)角線(xiàn)依次檢查，這些對(duì)角線(xiàn)是由兩個(gè)相鄰三角構(gòu)成的，最終成功的最優(yōu)三角剖分(如圖2-D).

圖2 徑向掃描算法圖解Fig.2 Diagram of the radial scan algorithm

以上過(guò)程表明:徑向掃描算法的基本思想是先徑向掃描離散點(diǎn)，構(gòu)建初始扇形三角網(wǎng)，然后從扇形的任一點(diǎn)開(kāi)始按一定方向進(jìn)行凹側(cè)連接，直到把外邊界變成凸多邊形為止.用局部?jī)?yōu)化算法調(diào)整到最優(yōu).

2.2 在Delaunay三角網(wǎng)中嵌入約束邊

本文的初始三角網(wǎng)是通過(guò)Delaunay三角剖分算法中的徑向掃描算法獲得的，與帶約束關(guān)系的離散點(diǎn)相關(guān)，下面介紹將帶約束關(guān)系的約束邊嵌入到初始三角網(wǎng)的具體方法.

建立一個(gè)存在的三角形集合T(V;L)，MT?T(V;L)，由MT中三角形的外邊構(gòu)成的多邊形稱(chēng)為影響多邊形 Q={pi，V1，V2，…，VK，pj，…，pi};設(shè)要插入的約束邊為 l，l=pipj，pipj∈V，與約束邊 l相交的三角形所組成的區(qū)域即為約束邊l的影響域;

對(duì)于影響多邊形Q，F(xiàn)loria［4］給出了如下的性質(zhì):

①簡(jiǎn)單多邊形Q被對(duì)角線(xiàn)pipj分成Qu和Qd兩部分，且Qu和Qd也為簡(jiǎn)單多邊形，如圖3所示.

根據(jù)“配第—克拉克”定律，在一個(gè)國(guó)家或地區(qū)的經(jīng)濟(jì)發(fā)展過(guò)程中，符合產(chǎn)業(yè)結(jié)構(gòu)演變的基本規(guī)律是第一產(chǎn)業(yè)在國(guó)民經(jīng)濟(jì)中的比重會(huì)越來(lái)越少，而第二、第三產(chǎn)業(yè)所占比重將不斷加大。農(nóng)業(yè)在國(guó)民經(jīng)濟(jì)中的地位將會(huì)隨著工業(yè)化和現(xiàn)代化進(jìn)程被工業(yè)和服務(wù)業(yè)所逐步取代，這也是一個(gè)國(guó)家和地區(qū)的產(chǎn)業(yè)結(jié)構(gòu)從低級(jí)向高級(jí)優(yōu)化的必然過(guò)程。國(guó)內(nèi)學(xué)者楊琳、李建偉（2002）采用第二、第三產(chǎn)業(yè)增加值之和與當(dāng)年國(guó)民生產(chǎn)總值的比值（ISR）來(lái)衡量一個(gè)國(guó)家或者地區(qū)產(chǎn)業(yè)結(jié)構(gòu)的變化。用公式表示則如下所示：

②Qu和Qd能被三角剖分.

③VK∈Q，若 VK為到 pipj的最近點(diǎn)，且(VK≠pi，VK≠pj)，則一定有

piVK∈Q，VKpj∈Q

圖3 pipj作為約束邊的影響域Fig.3 The influence domain of pipj as constraint boundary

在Delaunay三角網(wǎng)內(nèi)插入約束邊pipj的算法如下［5］:

①?gòu)娜蔷W(wǎng)中找到以pi為頂點(diǎn)的三角形t1，且t1與pipj相交.如圖4所示.

圖4 首三角形的確定Fig.4 The determination of the first triangle

②從t1開(kāi)始通過(guò)三角形的拓?fù)湫畔⒌乃阉鳎l(fā)現(xiàn)影響域內(nèi)三角形t1、t2、…，tk，將這些三角形保存在影響域三角形陣列中，同時(shí)生成Qu和Qd影響域的邊界陣列，見(jiàn)圖5.

圖5 pipj邊的影響域的確定Fig.5 The determination of the influence domain of pipj boundary

③為了在影響域內(nèi)生成一個(gè)新的三角形，從一開(kāi)始就pipj使用遞歸算法，同時(shí)將三角形的新一代信息保存在影響域三角形陣列的三角形空間中.

3 漢陽(yáng)區(qū)含斷層地質(zhì)體模型的建立

3.1 工程背景

武漢市漢陽(yáng)圖幅位于經(jīng)度114°～114°15'，緯度30°30'～30°40'.區(qū)內(nèi)(如圖6)第四系分布十分廣泛(圖中地表區(qū)域即為第四系地質(zhì))，占圖幅總面積的90%以上，發(fā)育齊全，成因類(lèi)型多樣，通常按照時(shí)代四分法進(jìn)行劃分，即早更新世、中更新世、晚更新世和全新世.

3.2 漢陽(yáng)區(qū)建模過(guò)程及成果展示

漢陽(yáng)區(qū)含斷層地質(zhì)體三維模型(圖6)主要是由帶約束的Delaunay三角剖分得到的三角面構(gòu)建的.此次模型的構(gòu)建數(shù)據(jù)(如圖1)主要來(lái)自于:①地表地質(zhì)圖;②基巖地質(zhì)圖;③鉆孔;④含斷層線(xiàn)的剖面;⑤地表等高線(xiàn)和地表以下各地層等深線(xiàn)等.

根據(jù)漢陽(yáng)城市地質(zhì)的數(shù)據(jù)特點(diǎn)，此次建模分3個(gè)步驟:a.首先對(duì)建模區(qū)數(shù)據(jù)進(jìn)行預(yù)處理(包括鉆孔數(shù)據(jù)的整理、確定鉆孔標(biāo)準(zhǔn)分層、進(jìn)行剖面交叉處一致性檢查、建模區(qū)域單元格的劃分等)，生成初始數(shù)據(jù)點(diǎn)(如圖7);b.對(duì)數(shù)據(jù)點(diǎn)進(jìn)行Delaunay三角剖分，形成DT剖分三角形(如圖8);c.將地表等高線(xiàn)、地形線(xiàn)、剖面線(xiàn)(含斷層線(xiàn))等線(xiàn)形數(shù)據(jù)插入到第2步生成的三角網(wǎng)中，通過(guò)這些步驟完成對(duì)漢陽(yáng)區(qū)數(shù)據(jù)的約束三角剖分(如圖9)，最終生成漢陽(yáng)區(qū)含斷層地質(zhì)體的三維地質(zhì)模型(如圖6).

圖6 漢陽(yáng)區(qū)三維地質(zhì)體Fig.6 3D geological body in Hanyang District

圖7 地表離散點(diǎn)Fig.7 The surface of the discrete points

圖8 地表無(wú)約束的三角網(wǎng)Fig.8 Unconstrained surface triangulation

圖9 地表嵌入約束邊的三角網(wǎng)Fig.9 Constraint edge embedded in the surface triangulation

圖7～圖9都是以地表嵌入約束邊三角剖分為例，下面以剖面上的約束線(xiàn)為例，詳細(xì)闡述漢陽(yáng)區(qū)地質(zhì)體三維建模過(guò)程:

a.觀察剖面上有各種約束線(xiàn)，通過(guò)剖面可推斷此處有一個(gè)正斷層，如圖10加粗線(xiàn)即為斷層線(xiàn).

b.觀察基巖地質(zhì)圖和剖面圖相交部分(如圖11)，發(fā)現(xiàn)此處單元格①中斷層位于基巖面以下.

圖10 漢陽(yáng)區(qū)部分剖面形成的單元格①Fig.10 Part profile forming cell① in Hanyang

圖11 基巖面與剖面對(duì)照?qǐng)DFig.11 The comparison chart of bedrock surface and profile

將剖面上的約束線(xiàn)插入到初始三角網(wǎng)中，在剖面上人工交互建模，可先構(gòu)造斷層的下盤(pán)(如圖12)，構(gòu)建下盤(pán)時(shí)通過(guò)由外向內(nèi)或由內(nèi)向外層層構(gòu)建，最終整個(gè)斷層體構(gòu)建完畢，(如圖13).通過(guò)以上方法，可以對(duì)每個(gè)地層面逐一構(gòu)建，再通過(guò)手動(dòng)拼接輪廓線(xiàn)對(duì)每一地層的頂面與底面縫合，形成一個(gè)完整的地質(zhì)體模型(如圖14).如果想添加剖面分層數(shù)據(jù)或是更改對(duì)地層的認(rèn)識(shí)，即可將新的特征約束線(xiàn)插入新的初始三角網(wǎng)，重新構(gòu)建此單元格地質(zhì)體.

圖12 斷層下盤(pán)最外圈地質(zhì)體Fig.12 The outer ring geological body in the footwall

圖13 斷層地質(zhì)體三維圖Fig.13 3d geological fault

圖14 單元格①中完整地質(zhì)體Fig.14 Complete geological body in cell①

4 結(jié)語(yǔ)

通過(guò)對(duì)漢陽(yáng)區(qū)城市地質(zhì)數(shù)據(jù)信息的分析，建議選用帶約束的三角剖分先生成TIN(不規(guī)則三角網(wǎng))面，最終生成復(fù)雜三維城市地質(zhì)體的方法.采用這種方法，集合了各種點(diǎn)集和特征線(xiàn)(地表等高線(xiàn)、地表地質(zhì)圖上的地形線(xiàn)、剖面上的地層線(xiàn)和斷層線(xiàn)、基巖地質(zhì)圖上的基巖分界線(xiàn)和斷層線(xiàn)、各地層的等深線(xiàn)等)，并將這些數(shù)據(jù)作為約束點(diǎn)和約束線(xiàn)插入進(jìn)行帶約束的三角剖分.因?yàn)轳詈狭硕嘣磾?shù)據(jù)，所以模型的精度更高.本文算法已在Visual C++6.0中得以編程實(shí)現(xiàn)，并在武漢市漢陽(yáng)區(qū)城市地質(zhì)三維模型建立中得以應(yīng)用.為城市三維地質(zhì)的發(fā)展奠定了基礎(chǔ).

［1］ZHONG D H，LI M C，SONG L G，et al.Enhanced NURBS modeling and visualization for large 3D geoengineering applications:an example from the Jinping first level hydropower engineering project，China ［J］.Computer and Geosciences，2006，32(9):1270-1282.

［2］何趙云，李章林，張夏林.含斷層地質(zhì)體三維建模技術(shù)［J］.軟件導(dǎo)刊，2010，9(8):175-177.HE Zhao-yun，LI Zhang-lin，ZHANG Xia-lin.3d geological faults modeling technology ［J］.Software Guide，2010，9(8):175-177.(in Chinese)

［3］林報(bào)嘉，葉燕林，趙建偉.三角網(wǎng)的生成與簡(jiǎn)化算法［J］.地理學(xué)與國(guó)土研究，2001，17(1):69-72.LIN Bao-jia，YE Yan-lin，ZHAO Jian-wei.The estab-lishment ands simplification algorithm of TIN［J］.Geography and Territorial Research.2001，17(1):69-72.(in Chinese)

［4］Floriai LD.An on-line algorithm for constrained delaunay triangulation［J］.CAGIP:Graphical Models and Image Processing，1992，54(3):290-300.

［5］宋占峰，詹振炎，蒲浩.Delaunay三角網(wǎng)剖分中嵌入約束邊的局部調(diào)整算法［J］.西南交通大學(xué)學(xué)報(bào)，2002，37(4):399-403.SONG Zhan-feng，ZHAN Zhen-yan，PU Hao.A local adjustment algorithm for inserting constrained segments in delaunay triangulation［J］.Journal of Southwest Jiaotong University，2002，37(4):399-403.(in Chinese)