基于地層線的三維地質(zhì)剖切算法

黃 忠， 江巨浪， 張佑生

（1.安慶師范學(xué)院 物理與電氣工程學(xué)院，安徽 安慶 246011；2.合肥工業(yè)大學(xué) 計(jì)算機(jī)與信息學(xué)院，安徽 合肥 230009）

基于地層線的三維地質(zhì)剖切算法

黃 忠1，2， 江巨浪1， 張佑生2

（1.安慶師范學(xué)院 物理與電氣工程學(xué)院，安徽 安慶 246011；2.合肥工業(yè)大學(xué) 計(jì)算機(jī)與信息學(xué)院，安徽 合肥 230009）

文章以地質(zhì)鉆孔為離散點(diǎn)，運(yùn)用空間插值以及由曲面擬合得到的各地質(zhì)結(jié)構(gòu)趨勢面，建立了工程地質(zhì)三棱柱結(jié)構(gòu)模型，在此模型基礎(chǔ)上提出一種任意切面的剖切方法。該方法首先應(yīng)用三維空間求交算法快速求出剖切面與凹凸起伏的地層趨勢面的交點(diǎn)；然后將這些交點(diǎn)按照地層層序連接成三角形并根據(jù)其貫穿的地層屬性進(jìn)行精確著色，從而得到三維的剖切實(shí)體圖、切面圖，為地質(zhì)工作者在3D空間中觀察地質(zhì)結(jié)構(gòu)、分析地質(zhì)特征提供了新的方法和手段，可為地質(zhì)工程項(xiàng)目的設(shè)計(jì)、管理和施工提供科學(xué)的決策依據(jù)。

地質(zhì)結(jié)構(gòu)趨勢面；空間求交；剖面；傾斜剖切；剖切分析

一般地，各種地質(zhì)信息，包括地表地形、地層界面、斷層、地下水位、風(fēng)化層厚度分布以及各種物探化探資料，都可以通過野外測量獲得。這些信息是一些離散的數(shù)據(jù)，地質(zhì)工作者很難直接利用它們分析其在地質(zhì)體中的分布規(guī)律［1］。研究者們面對這些實(shí)測信息，利用趨勢面擬合或空間內(nèi)插的方法獲得地質(zhì)結(jié)構(gòu)趨勢面，再根據(jù)三維地質(zhì)體特征，采用屬性插值、地質(zhì)解譯、評估分析、統(tǒng)計(jì)預(yù)測等方法建立三維地質(zhì)模型，然后運(yùn)用可視化算法，對數(shù)字仿真模型進(jìn)行處理，繪制出具有真實(shí)感的三維地質(zhì)結(jié)構(gòu)圖。由于不同地層面是由規(guī)則網(wǎng)格分解成三角網(wǎng)組成，而這些三角網(wǎng)在水平面上的投影是相同的，這樣使得相鄰層對應(yīng)三角形和豎直方向上互相平行的3條棱邊就構(gòu)成了三棱柱結(jié)構(gòu)［1］，如圖1所示。

建立了工程區(qū)地質(zhì)體三棱柱模型后，將同層的三棱柱著以相同的顏色，不同層的三棱柱著以不同顏色，通過顏色反映地質(zhì)特征的變化，達(dá)到直觀描述地下復(fù)雜的地質(zhì)構(gòu)造情況，形象地表達(dá)地質(zhì)構(gòu)造的形態(tài)特征的目的。但為了展示地質(zhì)體內(nèi)部特征、空間特征并對其地質(zhì)特征進(jìn)行綜合分析，還需要對地質(zhì)體進(jìn)行三維的剖切，以便清晰地展現(xiàn)模型的內(nèi)部細(xì)節(jié)，最大限度地增強(qiáng)地質(zhì)分析的直觀性和準(zhǔn)確性［2］。因此，三維工程地質(zhì)模型剖切的實(shí)現(xiàn)一直是研究的熱點(diǎn)［3］。

圖1 地質(zhì)體三棱柱結(jié)構(gòu)

目前，工程地質(zhì)三維剖切主要有2種途徑：

（1）求構(gòu)成三維地質(zhì)模型的不規(guī)則三角形面與空間剖切平面的交線，這種算法需要搜索所有的面與線，計(jì)算量大、時間長［4］。

（2）采用投影方法，將三維離散鉆孔點(diǎn)集轉(zhuǎn)化為二維平面點(diǎn)集［5］，這種算法違反Delaunay三角剖分“最大最小角”優(yōu)化準(zhǔn)則［6］。

基于上面2種方法的特點(diǎn)，本文提出基于地層線的三維地質(zhì)剖切算法，以實(shí)現(xiàn)任意剖面的三維剖切分析，這需要有效地解決如下2個問題：

（1）快速求出剖切面與凹凸起伏的地層面的交點(diǎn)。

（2）根據(jù)不同地層線的交線繪制出反映地層屬性的剖面。

1 地層面與剖切面求交算法

考慮地質(zhì)體的三棱柱結(jié)構(gòu)，對三維地質(zhì)體的剖切變成對三棱柱的剖切，其實(shí)質(zhì)是對剖面與上下三角面片和棱柱側(cè)面求交［7］。但地質(zhì)體中含有大量的三棱柱，每個三棱柱將與剖切面進(jìn)行5次面面求交測試［8］，其計(jì)算量在交互環(huán)境中是無法接受的。因此，當(dāng)求地層面與剖切面的交線時，首先用三角網(wǎng)格逼近地層面，然后用這些近似三角網(wǎng)格與剖面的交線來近似表示地層面和剖切面的交線，如圖2所示［9］。

圖2 地層面與剖切面交線

設(shè)剖切面方程為ax＋by＋cz＋d＝0，某一點(diǎn)M（X，Y，Z）與剖面位置關(guān)系可以通過將其坐標(biāo)值代入剖面方程來判定［10］。設(shè)D＝aX＋bY＋cZ＋d，則有：

（1）若D＞0，M點(diǎn)處在平面正法向一側(cè)。

（2）若D＝0，M點(diǎn)處在平面上。

（3）若D＜0，M點(diǎn)處在平面負(fù)法向一側(cè)。

檢測地層面網(wǎng)格每個三角面片的3個頂點(diǎn)與剖切面的位置關(guān)系，若D值均大于0或均小于0，則三角面片與剖切平面不相交；否則，進(jìn)行求交運(yùn)算，這時分別對該三角面片的邊與剖切平面進(jìn)行求交判定，若一條邊的2個頂點(diǎn)的D值一個大于0而另一個小于0，則由該邊的方程與剖切面方程聯(lián)立求解可得交點(diǎn)P（x，y，z）。

在地層面網(wǎng)格上定義u、v2個方向，并劃分為（u－num＋1）×（v－num＋1）的網(wǎng)格，剖切面與地層面交線在水平面上投影如圖3所示，圖中實(shí)線表示剖切面與地層面網(wǎng)格u方向所有平行線相交，虛線表示剖切面與地層面網(wǎng)格u方向部分平行線相交。當(dāng)u方向平行線與剖切面沒有交點(diǎn)時（如u＝u2等），需要將v方向的平行邊界線與剖切面的交點(diǎn)（圖3中P點(diǎn)）作為u＝u2上交點(diǎn)。具體算法如下。

其次，要注意拓展視聽材料的來源。筆者在一學(xué)期的教學(xué)實(shí)踐中使用的視聽材料大多來自《看視頻學(xué)英語口語-衣食住行一本通》《悅讀聯(lián)播》《英美文化讀本》中的視頻、音頻材料和英語紀(jì)錄片，對英語電影和英語歌曲的使用較少。根據(jù)問卷結(jié)果，學(xué)生希望增加英語電影和英語歌曲的內(nèi)容。因此，在教學(xué)內(nèi)容的選擇上應(yīng)該更多地考慮高中學(xué)生年齡特點(diǎn)、學(xué)習(xí)興趣和認(rèn)知水平。

圖3 剖切面與地層面交線在水平面上投影

1.1 求地層網(wǎng)格與剖面在u方向的交點(diǎn)

對于地層網(wǎng)格u方向上每一條平行線ui（i＝0，1，…，u－num），取v方向上相鄰2點(diǎn)Pui，vj（x1，y1，z1）、Pui，vj＋1（x2，y2，z2），計(jì)算D1、D2及它們的乘積Flag。

如果Flag＞0，表明Pui，vj、Pui，vj＋1位于剖切面同側(cè)，需要繼續(xù)計(jì)算下2個相鄰點(diǎn)。

如果Flag＜0，表明Pui，vj、Pui，vj＋1點(diǎn)位于剖切面的兩側(cè)，可求出線段（Pui，vj，Pui，vj＋1）與剖面的交點(diǎn)P。

如果一條網(wǎng)格線ui上所有點(diǎn)的Flag大于0，表明該網(wǎng)格線與剖切面沒有交點(diǎn)，則計(jì)算v方向上的邊界線與剖切面的交點(diǎn)，即求線段（Pui－1，v－num（x1，y1，z1），Pui，v－num（x2，y2，z2））與剖切面交點(diǎn)，并將此交點(diǎn)作為該網(wǎng)格線的交點(diǎn)Pu。

1.2 確定剖面與地層的剖切情況

計(jì)算出了剖切面跟地層面的交點(diǎn)后，還需要考察它在豎直方向上的相交情況。如果剖切面是豎直的，它與各個地層面交點(diǎn)僅高程值不相同（如圖4），此時只需要將各個層面對應(yīng)交點(diǎn)相連成三角形。

如果剖切面是傾斜的，它與各個地層面的交點(diǎn)在水平面上的投影不相同，不僅要計(jì)算剖切面與地層面的交點(diǎn)，還要計(jì)算剖切面在地層內(nèi)的交點(diǎn)。當(dāng)剖切面與某個地層線沒有交點(diǎn)時，將其與邊界的交點(diǎn)作為地層線交點(diǎn)。如圖5中，剖切面與L3層地層線沒有交點(diǎn)，將邊界線（P1，P2）與地層線的交點(diǎn)P做為L3層地層線交點(diǎn)，而L3以下的地層線不需要再計(jì)算交點(diǎn)。計(jì)算出所有層內(nèi)交點(diǎn)后，其剖面的繪制方式與豎直剖面的連接方式相同。

圖4 豎直剖切（各層交點(diǎn)僅高程值不同）

圖5 傾斜剖切（貫穿部分地層）

2 基于地層線的三角形渲染

2.1 三角形與地層線的位置關(guān)系

當(dāng)確定了地層面與剖面在各個地層中的交點(diǎn)后，繪制這些交點(diǎn)形成的三角形以及準(zhǔn)確對其根據(jù)貫穿的地層的顏色進(jìn)行著色也是進(jìn)行地質(zhì)剖切需要解決的問題［10］。具體分析如下：

（2）由于三角形的頂點(diǎn)可能屬于不同地層，中間可能跨越其他地層，這需要將三角形按照地層線分割成多個部分，使每個部分只屬于1個地層。三角形的3個頂點(diǎn)與地層線的位置關(guān)系可能有如下幾種：①3個頂點(diǎn)都在同一層，這時不需要進(jìn)行三角形分割，如圖6a所示；②3個頂點(diǎn)分處2層，如圖6b所示，這需要將三角形分割成3部分：i層部分、i＋1層至j－1層部分（可能沒有）、j層部分；③3個頂點(diǎn)分處于3層，這時地層線與三角形的位置關(guān)系如圖6c所示，這時需要將三角形分割成5個部分：i層部分、i＋1層至j－1層部分（i、j相鄰時沒有）、j層部分，j＋1層至k－1層部分（j、k相鄰時沒有）、k層部分。

（3）在OpenGL中，三角形內(nèi)部點(diǎn)的顏色值是通過3個頂點(diǎn)的顏色線性插值得到的。但是，對于表示地質(zhì)體剖切面的三角形來說，其內(nèi)部各點(diǎn)都應(yīng)按所屬地層的顏色繪制，而不是通過插值計(jì)算。這要求將三角形按照地層線分割成多個部分，使每個部分屬于同一層。

圖6 3頂點(diǎn)在地層中位置關(guān)系

基于以上考慮，本文提出一種基于地層線的三角形渲染算法。該算法的基本思想是應(yīng)用掃描線法的基本原理，將地層線做為掃描線對三角形進(jìn)行分割，利用掃描線與三角形的交點(diǎn)及交點(diǎn)間的地層屬性來完成掃描線間的顏色填充。

2.2 三角形分割繪制算法

通過上面的分析可知，如何實(shí)現(xiàn)三角形自動分割是需要著力解決的問題。如上所述，三角形最多可能需要分割成5個部分，據(jù)此可設(shè)計(jì)出能夠處理各種情況的三角形分割繪制算法。首先將三角形的3個頂點(diǎn)按照所在地層順序從低到高進(jìn)行排序，如果遇到2個頂點(diǎn)處在同一層情況，則按照高程大小排序。對于當(dāng)前三角形，將排好序后的3個頂點(diǎn)分別記為p1、p2、p3，其對應(yīng)的層號為L1、L2、L3，用flag2和flag3分別表示p2與p3是否已經(jīng)處理。算法描述如下：

（1）繪制三角形（p1，p2，p3）在L1層的部分。判斷（p1，p2）、（p1，p3）是否與L1＋1有交點(diǎn)，若有交點(diǎn)，記為p12、p13，則用L1層地層顏色繪制三角形（p1，p12，p13），如圖7a所示；否則，判斷（p2，p3）與L1＋1是否有交點(diǎn)，如果（p2，p3）與L1＋1沒有交點(diǎn)，即p1、p2、p3在同一層，用L1層地層顏色繪制三角形（p1，p2，p3），轉(zhuǎn)步驟（6），否則，用L1層顏色繪制2個三角形（p1，p2，p23），（p1，p23，p13），flag1＝true（如圖7b所示），轉(zhuǎn)步驟（4）。

（2）繪制L1＋1到L2之間部分。對于L1＋1到L2之間各層L，求出（p1，p2）與L層、L1＋1層的交點(diǎn)p12、p12′，以及（p1，p3）與L層、L1＋1層的交點(diǎn)p13、p13′，用L層地層顏色繪制2個三角形（p12，p12′，p13）、（p13，p12′，p13′）。

（3）繪制三角形（p1，p2，p3）在L2層的部分。求出（p1，p2）與（p1，p3）與L2層之間的交點(diǎn)p12、p13，判斷（p2，p3）、（p1，p3）是否與L2＋1有交點(diǎn)：若都存在交點(diǎn)p23、p13′，則用L2層顏色繪制3個三角形（p12，p2，p23）、（p12，p23，p13）、（p13，p23，p13′），如圖8a所示；否則，用L2層顏色繪制2個三角形（p12，p2，p3）、（p12，p3，p13），如圖8b所示，轉(zhuǎn)步驟（6）。

圖7 繪制三角形（p1，p2，p3）在L1 層的部分

圖8 繪制三角形（p1，p2，p3）在L2 層的部分

（4）繪制L2＋1到L3之間部分。對于L2＋1到L3之間各層L，求出（p2，p3）與L層、L1＋1層的交點(diǎn)p23、p23′，以及（p1、p3）與L層、L1＋1層的交點(diǎn)p13、p13′，用L層地層顏色繪制三角形（p23，p23′，p13′）、（p23，p13′，p13）。

（5）繪制L3層的點(diǎn)。求出（p1，p3）、（p2，p3）與L3的交點(diǎn)p13、p23，用L3層地層顏色繪制三角形（p23，p3，p13），如圖9所示。

（6）結(jié)束。圖7～圖9中，點(diǎn)虛線表示可能存在的地層，粗虛線表示將多邊形分解成三角形的分界線。

圖9 繪制L3層的點(diǎn)

3 實(shí)驗(yàn)結(jié)果

在PC機(jī)上，采用Visual C＋＋和OpenGL編程，對上述算法進(jìn)行了實(shí)現(xiàn)，并應(yīng)用該算法對某區(qū)工程地質(zhì)模型進(jìn)行了三維剖切。該區(qū)域共收集10個鉆孔，包含8個地層，經(jīng)數(shù)據(jù)處理后10個鉆孔的局部位置坐標(biāo)和每個鉆孔貫穿各個地層面高程數(shù)據(jù)見表1、表2所列。以表中10個鉆孔點(diǎn)為離散點(diǎn)，運(yùn)用空間插值以及由曲面擬合得到的各地質(zhì)結(jié)構(gòu)趨勢面，獲得地質(zhì)體的三棱柱結(jié)構(gòu)模型，然后運(yùn)用可視化算法，對該模型進(jìn)行處理，繪制出具有真實(shí)感的三維地質(zhì)結(jié)構(gòu)圖，如圖10、圖11所示。為了更好地了解地質(zhì)內(nèi)部地質(zhì)結(jié)構(gòu)，對地質(zhì)體結(jié)構(gòu)進(jìn)行三維剖切，圖12、圖13分別為豎直剖切后的剖切面結(jié)構(gòu)和剖面分離地質(zhì)體的結(jié)果圖，圖14所示為貫穿所有地層的傾斜剖切面產(chǎn)生的剖切圖。圖15所示為貫穿部分地層的傾斜剖切面產(chǎn)生的剖切圖。

表1 10個鉆孔的局部位置坐標(biāo)

表2 鉆孔貫穿各個地層面高程數(shù)據(jù)

圖10 地質(zhì)體三棱柱網(wǎng)格

圖13 豎直剖切

圖11 三維地質(zhì)結(jié)構(gòu)圖

圖12 剖面圖

圖14 貫穿所有地層的傾斜剖切

圖15 貫穿部分地層的傾斜剖切

4 結(jié)束語

采用三棱柱建模方法構(gòu)造了真實(shí)感的三維地層模型，并應(yīng)用三維空間求交算法和基于地層線的三角形繪制算法，可由任意切面對三維地質(zhì)模型進(jìn)行剖切，得到三維的剖切實(shí)體圖、切面圖，使不熟悉地質(zhì)結(jié)構(gòu)和構(gòu)造復(fù)雜性的人對地質(zhì)空間關(guān)系有一個直觀的認(rèn)識，并提高對難以想象的復(fù)雜地質(zhì)條件的理解和判別，對實(shí)際的工程地質(zhì)分析和解釋具有指導(dǎo)意義，可以很好地滿足層狀地質(zhì)結(jié)構(gòu)實(shí)際應(yīng)用的需要。該算法對地層相互穿插、出現(xiàn)斷層等復(fù)雜地質(zhì)結(jié)構(gòu)的三維剖切效果較差，且未考慮復(fù)雜地質(zhì)構(gòu)造（如侵入體、復(fù)雜褶皺等）的影響與作用，因此仍需隨著實(shí)際應(yīng)用的深入而作進(jìn)一步研究。

［1］劉 喜.三維地質(zhì)建模關(guān)鍵算法研究［D］.合肥：中國科學(xué)技術(shù)大學(xué)，2011.

［2］朱良峰，潘 信，吳信才.地質(zhì)斷層三維可視化模型的構(gòu)建方法與實(shí)現(xiàn)技術(shù)［J］.軟件學(xué)報，2008，19（8）：2004－2017.

［3］李 楠，曾俊鋼，肖克炎，等.一種線框矢量模型表示的地質(zhì)體折面剖切方法［J］.測繪科學(xué)；2011，36（6）：28－31.

［4］梁師俊.基于虛擬鉆孔的工程地質(zhì)三維剖切的實(shí)現(xiàn)［J］.地理與地理信息科學(xué)，2010，26（1）：11－14.

［5］郭加樹，張 勇，畢世普，等.一種三維地層模型可視化及剖切分析方法［J］.海洋地質(zhì)動態(tài)，2010，26（11）：51－54.

［6］鄒 偉.三維地質(zhì)勘查軟件核心技術(shù)研究與實(shí)現(xiàn)［D］.北京：中國地質(zhì)科學(xué)院，2010.

［7］De Ricazo F，Oreste P，Pelizza S.Subsurface geologicalgeotecbnical modelling to sustain underground civil planning［J］.Engineering Geology，2008，96（1）：187－204.

［8］Smirnoff A，Boisvert E，Paradis S J.Support vector machine for 3Dmodelling from sparse geological information of various origins［J］.Computers ＆ Geosciences，2008，34（2）：127－143.

［9］黃地龍，鄧 飛.復(fù)雜地層結(jié)構(gòu)模型三維重構(gòu)與可視化方法研究［J］.成都理工大學(xué)學(xué)報，2008，35（5）：553－558.

［10］陳國良，曹衛(wèi)群，黃心淵.一種由等高線模型生成規(guī)則格網(wǎng)模型 的 算 法 ［J］.中 國 圖 象 圖 形 學(xué) 報，2007，12（6）：1110－1113.

3D sectioning algorithm of geology based on stratum line

HUANG Zhong1，2， JIANG Ju－lang1， ZHANG You－sheng2

（1.School of Physics and Electronic Engineering，Anqing Normal University，Anqing 246011，China；2.School of Computer and Information，Hefei University of Technology，Hefei 230009，China）

Three－prism structure model of engineering geology is built by taking geological drillings as discrete points and applying each trend surface of geological structures，which are calculated by the spatial interpolation and surface fitting algorithm.On this basis，a sectioning method of arbitrary cutting is proposed.First the 3Dspace intersection algorithm is used to quickly calculate the intersections of cutting plane and rugged trend surface，and then in terms of stratigraphic sequence these intersections create triangles，which are colored accurately according to cross－cutting strata properties.As a result，the entity diagram and section diagram of 3Dcutting are obtained，which can give new means for geologists when they observe geological structures and analyze geological features in 3Dspace，thus providing a scientific basis for the decision－making of the design，management and construction of geological projects.

trend surface of geological structure；space intersection；section plane；inclined cutting；sectioning analysis

TP391.41

A

1003－5060（2013）02－0202－06

10.3969／j.issn.1003－5060.2013.02.016

2012－09－21；

2012－10－31

國家自然科學(xué)基金資助項(xiàng)目（60575023）；安徽省自然科學(xué)基金資助項(xiàng)目（090412065）和安慶師范學(xué)院青年科研基金資助項(xiàng)目（KJ201007）

黃 忠（1981－），男，安徽安慶人，安慶師范學(xué)院講師，合肥工業(yè)大學(xué)博士生；

江巨浪（1966－），男，安徽安慶人，博士，安慶師范學(xué)院教授；

張佑生（1941－），男，湖南瀏陽人，合肥工業(yè)大學(xué)教授，博士生導(dǎo)師.

（責(zé)任編輯 馬國鋒）