何 畏,吳文鸝
(中國(guó)地質(zhì)科學(xué)院 地球物理地球化學(xué)勘查研究所,河北廊坊 065000)
目前,地球物理數(shù)據(jù)解釋以二維剖面為主[1、2],而實(shí)際地質(zhì)體在空間賦存為三維結(jié)構(gòu)。為了深入認(rèn)識(shí)地質(zhì)體的構(gòu)造特征,需要將二維剖面信息轉(zhuǎn)化為三維信息來(lái)完成地質(zhì)體的推斷解釋,即利用連接截面輪廓線模擬實(shí)際地質(zhì)體的幾何形狀[9、10],從而提高數(shù)據(jù)處理、分析,以及解釋工作的效率與準(zhǔn)確性。
自 1993年 Simon WHoulding[5]提出三維地質(zhì)建模這一概念以來(lái),其實(shí)現(xiàn)方式已呈現(xiàn)多樣化,如基于 Delaunay四面體模型、剖面模型、三角形多面體模型等建模方式[6~8],輪廓線三維地質(zhì)建模方式也應(yīng)用到了地質(zhì)領(lǐng)域。
對(duì)于眾多輪廓線建模的研究成果,由 H.Fuchs提出的“最小表面積”和 E.Kepple主張的“最大體積”連接算法[3、4]最為有名。針對(duì)以往重構(gòu)表面算法的復(fù)雜、計(jì)算量大等不利因素,作者在文中提出了最優(yōu)路徑 (optimalpath,OP)算法,從而實(shí)現(xiàn)了三角形多面體截面輪廓線三維地質(zhì)體建模。
輪廓線建模的基本思想是,用一系列彼此不相交,互相不重疊的三角面片,將相鄰截面的輪廓線連接起來(lái)的方法。H.Fuchs[3]指出,只有滿足下列二個(gè)條件的三角面片集合,連接才是合理的。
(1)由于每一個(gè)輪廓線的線段必須并且只能在一個(gè)基本三角面片中出現(xiàn)。因此,二輪廓線連接的三維表面模型包括 m+n個(gè)基本三角形 (m、n分別為二輪廓線上節(jié)點(diǎn)數(shù))。
(2)如果一個(gè)跨距在某一基本三角面中為左跨距,則該跨距是而且僅是另一基本三角面片的右跨距,如下頁(yè)圖 1(a)所示。
見(jiàn)下頁(yè),圖 1(a)是位于不同平面的 P、Q二輪廓線頂點(diǎn)按一定規(guī)則連接形成的三維形體效果圖,圖 1(b)為輪廓線上點(diǎn)列連接關(guān)系有向圖。假設(shè)圖 1(a)上、下二輪廓線點(diǎn)列分別為 P0,P1,…,Pm-1及 Q0,Q1,…,Qn-1,則可以用一個(gè) m行、n列的有向圖來(lái)表示點(diǎn)列之間的連接關(guān)系。如有向圖G(V,E)(如圖 1(b)所示 ),其中 Vij表示點(diǎn) Pi與點(diǎn)Qj之間的跨距,Eij表示點(diǎn) Pi與點(diǎn) Qj之間對(duì)應(yīng)的邊。
圖1 輪廓線連接和矩陣圖Fig.1 Map of contour data connection and matrix
從圖 1(b)中可以看出,一種輪廓線建模的三角形連接方式對(duì)應(yīng)于從點(diǎn) V00到點(diǎn) Vm-1,n-1的一條路徑,因此可以滿足上述二個(gè)條件的連接表面總數(shù)等于由點(diǎn) V00到 Vm-1,n-1時(shí)可能存在的路徑總數(shù),記為A[m,n]。
由式 (1)知,代表路徑總數(shù)的 A[m,n]隨著m、n的增加而快速增加[11],然而在這些路徑中,并不是所有路徑對(duì)應(yīng)的連接方式都能滿足連接要求。為了從圖 1(b)中找到合適的路徑來(lái)滿足連接要求,H.Fuchs采用最小表面積,E.Kepple主張最大體積法來(lái)搜索全局最優(yōu)路徑,并取得了較為理想的連接效果。
從理論上講,在圖 1(b)中采用窮舉搜索法,一定可以找出一條可以接受且滿足輪廓線連接要求的全局最優(yōu)路徑,但當(dāng)二輪廓線節(jié)點(diǎn)較多或多輪廓線連接時(shí),搜索時(shí)間必然大幅增加,這勢(shì)必降低工作效率。
作者在本文提出的 OP算法,具有局部計(jì)算和決策啟發(fā)式特點(diǎn),能解決上述出現(xiàn)的不足。該算法不要求實(shí)現(xiàn)全局最優(yōu),而是基于局部計(jì)算來(lái)決定當(dāng)前的選擇。OP算法的“最優(yōu)”主要體現(xiàn)以下二個(gè)方面:①二輪廓線的初始連接節(jié)點(diǎn)是所有跨距最短的兩節(jié)點(diǎn);②在輪廓線節(jié)點(diǎn)矩陣圖中,由“最短對(duì)角線法則”來(lái)決定連接路徑的下一節(jié)點(diǎn)。具體算法描述如下:
(1)賦予圖 1(b)中節(jié)點(diǎn) Nij的初值為 Vij,Vij表示點(diǎn) Pi與點(diǎn) Qj之間的跨距 (i=0,1,…,m-1;j=0,1,…,n-1),求出 V中的最小值 Vopti_i,opti_j,并將Vopti_i,opti_j作為新的V′00值,在P、Q輪廓線上按逆時(shí)針對(duì)節(jié)點(diǎn)排序,形成新的有向路徑圖 G′(由于 G和G′形式一樣,故將圖 1(b)作為新的有向路徑圖G′)和節(jié)點(diǎn)跨距 V′ij。該步保證了在 P、Q輪廓線上第一條連接線直線距離最短,實(shí)現(xiàn)局部最優(yōu)。
(2)對(duì) G′中的每個(gè)節(jié)點(diǎn)賦一方向?qū)傩?對(duì)第 m-1行節(jié)點(diǎn)賦“向右”屬性,對(duì)第 n-1列節(jié)點(diǎn)賦“向下 ”屬性 ,對(duì) V′m-1,n-1節(jié)點(diǎn)賦“終止 ”屬性 ,其余節(jié)點(diǎn)賦“向右”和“向下”屬性。
(3)根據(jù) G′中節(jié)點(diǎn) N′ij的方向?qū)傩?獲取V ′i+1,j和 V ′i,j+1,并比較二者大小 ,將較小值對(duì)應(yīng)的節(jié)點(diǎn)作為路徑的下一節(jié)點(diǎn) (如下頁(yè)圖 2(a)所示),將 P、Q相應(yīng)節(jié)點(diǎn)編號(hào)存入容器 V_P、V_Q中,以備后續(xù)三角形連接使用 (如下頁(yè)圖 2(b)所示)。該步遵循了“最短對(duì)角線”局部最優(yōu)原則。
(4)判斷執(zhí)行步驟 (3)后的新路徑節(jié)點(diǎn)是否為終止節(jié)點(diǎn)V′m-1,n-1,如果為“真”則轉(zhuǎn)至步驟(5);如果為“假”則轉(zhuǎn)向步驟 (3),根據(jù)節(jié)點(diǎn)的方向?qū)傩?判斷其“向右”或“向下”節(jié)點(diǎn)的 V′值。重復(fù)步驟(3)~步驟(4),直至路徑終止節(jié)點(diǎn)V′m-1,n-1。
(5)將在上述過(guò)程中裝入容器 V_P或 V_Q中的節(jié)點(diǎn)編號(hào),與下一位置中的節(jié)點(diǎn)編號(hào)做“相等”判斷,并根據(jù)判斷結(jié)果,選擇 PPQ或 PQQ形式形成輪廓線之間的連接三角面片,算法結(jié)束。
此外,根據(jù) V_P和 V_Q容器形成的三角形個(gè)數(shù)比二輪廓線側(cè)面所需連接三角面片個(gè)數(shù)少二個(gè)的特點(diǎn),需要在步驟 (5)的基礎(chǔ)上追加二個(gè)三角形。
最優(yōu)路徑 (OP)算法的特點(diǎn),是基于連接的二條輪廓線大小、形狀和中心距相差不遠(yuǎn)時(shí),會(huì)取得滿意的連接效果。但當(dāng)二輪廓線不滿足其特點(diǎn)時(shí),可能會(huì)造成連接失敗,如圖 3所示出現(xiàn)連接失敗形成二棱錐體的情況。對(duì)于這種情況,作者給出如下二種解決方案。
圖2 最短對(duì)角線連接及節(jié)點(diǎn)編號(hào)存儲(chǔ)圖Fig.2 Map of the shortest diagonal connection and storageof node number
圖3 最優(yōu)算法連接失敗示意圖Fig.3 Sketch map of the failure of op timal algorithmconnection
(1)在連接失敗的二截面之間,增加一個(gè)或多個(gè)截面,使截面輪廓線中心差減小,從而滿足 OP算法的特點(diǎn)?;蛘呃萌藱C(jī)交互方式移動(dòng)二輪廓線上的節(jié)點(diǎn)及邊,使二輪廓線所形成區(qū)域的形狀、面積大小盡可能的接近,使中心差減小,最終取得較好的建模效果。
(2)在構(gòu)造三角面片之前,將該 P、Q輪廓線變換至以同一原點(diǎn)為中心的單位正方形之內(nèi),從而保證在單位正方形內(nèi)二輪廓線大小和形狀相近,并使對(duì)中情況較好[12]。由于將二輪廓線變換至單位正方形內(nèi),OP算法調(diào)用以及將各輪廓線變換到原來(lái)位置的反變換所需工作量大,所以本方案有時(shí)間復(fù)雜度較大的特點(diǎn)。
本文的算法在 C++編程環(huán)境下,借助于 Open-GL圖形庫(kù),編寫了基于截面輪廓線人機(jī)交互的建模程序。
圖4 截面輪廓線三維建模示意圖Fig.4 Sketch map of 3D modeling of cross - sectioncontour data
圖 4為六個(gè)截面 (在從左到右方向上,第二個(gè)截面只剩下透鏡體的截面,其余地質(zhì)屬性體截面已刪去)的部份輪廓線調(diào)用最優(yōu)路徑 (OP)算法形成的三維地質(zhì)體模型。圖 4(a)是六個(gè)截面形成的3D效果圖,每個(gè)截面由三種不同灰度的區(qū)域構(gòu)成,其中某區(qū)域含有一個(gè)透鏡體截面 (在三維空間中,Z方向?yàn)榇怪狈较?X方向?yàn)闁|西方向,Y軸所在平面為水平面)。在圖 4(b)中,某透鏡體在第一截面和第六截面 (在從左到右方向上)上處于尖滅狀態(tài),“透鏡體輪廓線連接圖”是對(duì)透鏡體截面輪廓線調(diào)用OP算法的三維建模圖,取得了較為理想的三角形連接效果,形象地模擬了透鏡體在巖層中的空間分布。
三維截面輪廓線人機(jī)交互建模是 3D建模的發(fā)展方向,雖然研究難度較大,但它可以在人與數(shù)據(jù)、人與圖形之間進(jìn)行“對(duì)話”,特別是與地球物理正演、反演技術(shù)結(jié)合,可提高地質(zhì)數(shù)據(jù)的處理、解釋效率,具有很好的應(yīng)用前景。
[1] 丁艷紅,夏東嶺,李志勇,等.對(duì)二維地震剖面的解釋與成圖方法的認(rèn)識(shí) [J].石油物探,2002,41(2):211.
[2] 張亞敏,張書法.二維地震疊加偏移剖面“二步法偏移”地質(zhì)解釋方法[J].西安石油大學(xué)學(xué)報(bào) (自然科學(xué)版),2009,24(1):34.
[3] FUCHS H, KEDEM ZM, USELTON S P. Op timal Surface Reconstruction from Planar Contours[ J ]. Communicationof the ACM, 1977, 20 (10) : 693.
[4] KEPPEL E. App roximating Comp lex Surface by Triangulation of Contour L ines [ J ]. IBM Journal Res Develop1975, 1: 2.
[5] HOULD ING SW. 3D GeosciencesModeling - Computer Techniques for Geological Characterization [M ]. BerlinHeidelberg: Sp ringer– Verlag, 1994.
[6]V ICTOR JD TSA I. Delaunay Triangulations in TIN Creation:An overview and a linear - time algorithm [ J ]. International Journal of Geographical Information Systems,1993, 7 (6) : 501
[7] 屈紅剛,潘懋,王勇,等.基于含拓?fù)淦拭娴娜S地質(zhì)建模[J].北京大學(xué)學(xué)報(bào) (自然科學(xué)版),2006,42(6):717.
[8] 吳文鸝,管志寧,高艷芳,等.重磁異常數(shù)據(jù)三維人機(jī)聯(lián)作模擬[J].物探化探計(jì)算技術(shù),2005,27(3):227.
[9] 李光亮,肖海紅,徐遵義,等.平行輪廓線構(gòu)建復(fù)雜斷層地質(zhì)模型研究[J].煤田地質(zhì)與勘探,2007,35(2):22.
[10]陳華根,吳健生,王家林,等.剖面疊置思想在地球物理反演建模中的應(yīng)用 [J].同濟(jì)大學(xué)學(xué)報(bào),2003,31(11):1309.
[11]唐澤圣.三維數(shù)據(jù)場(chǎng)可視化 [M].北京:清華大學(xué)出版社,2000.
[12]CHR ISTIANSEN H N, SEDERBERG TW. Conversion of Comp lex Contour L ine Definitions into Polygonal ElementMosaics [ J ]. Computer Graphics, 1978, 12 ( 2 ) :187.