基于Delaunay三角剖分的二維交互建模研究

楊 強(qiáng)

(中國(guó)石化 石油物探技術(shù)研究院，南京 211103)

0 引言

地震勘探中，正演模擬不但在地震數(shù)據(jù)采集中得到應(yīng)用，在地震資料處理和地震資料解釋中也是重要的驗(yàn)證技術(shù)手段，是進(jìn)行地震反演的基礎(chǔ)。而正演模擬的基礎(chǔ)就是需要準(zhǔn)確且合理的地質(zhì)模型，因此建模是地震正演的重要基礎(chǔ)。傳統(tǒng)上地質(zhì)模型都沿用Cenveny提出的模型結(jié)構(gòu)，即所謂層狀結(jié)構(gòu)模型。它要求每一個(gè)分界面都必須從模型體的左邊界貫穿到模型體的右邊界，分界面按順序由上到下依序排列，不得交叉。對(duì)于逆斷層、尖滅、透鏡體等復(fù)雜地層情況，只能人為地簡(jiǎn)化模型，從而滿(mǎn)足層狀模型。

在實(shí)際的地球物理勘探中，層狀結(jié)構(gòu)模型有兩個(gè)嚴(yán)重不足:①在采集、處理、解釋各階段的模型大都有逆斷層、尖滅、透鏡體等復(fù)雜地質(zhì)元素,層狀結(jié)構(gòu)模型無(wú)法準(zhǔn)確描述復(fù)雜的地質(zhì)拓?fù)浣Y(jié)構(gòu)，從而無(wú)法得到網(wǎng)格化模型;②處理中的疊前深度偏移速度建模以及地震解釋構(gòu)造建模等，都要求能夠交互修改模型，交互編輯中的層位修改及移動(dòng)必須遵循層狀規(guī)則，這對(duì)交互式復(fù)雜建模而言很困難，也不便利。蔣先藝[1]提出用點(diǎn)、段、線、面的概念描述二維封閉結(jié)構(gòu)模型，實(shí)現(xiàn)了二維復(fù)雜地質(zhì)結(jié)構(gòu)，其優(yōu)點(diǎn)是地震建模的數(shù)據(jù)層位、斷層、透鏡體都可以用點(diǎn)、線、面集來(lái)描述。這種方法比較復(fù)雜，難點(diǎn)就在于封閉塊體的追蹤。由于建模中的基本元素是點(diǎn)和線段，由這些基本元素需要得到地質(zhì)結(jié)構(gòu)的拓?fù)潢P(guān)系，采用三角剖分是一種理想的處理手段。

1 方法及思路

Delaunay 三角剖分是二維平面內(nèi)的最優(yōu)三角剖分，它在有限元分析、信息可視化、計(jì)算機(jī)圖形學(xué)等應(yīng)用領(lǐng)域有著重要應(yīng)用[2]。Ruppert的二維高質(zhì)量網(wǎng)格生成算法是第一個(gè)理論上保證網(wǎng)格劃分算法在實(shí)踐中真正令人滿(mǎn)意的算法[3]。Refine Delaunay 三角化方法解決了保邊界和內(nèi)嵌邊界的問(wèn)題,該方法往往通過(guò)在保留邊(約束邊)上加入新的節(jié)點(diǎn)以實(shí)現(xiàn)保邊界的目的[4]。筆者通過(guò)對(duì)Delaunay三角剖分的研究，以點(diǎn)和線段為基礎(chǔ)，對(duì)層位、斷層、透鏡體進(jìn)行三角剖分，采用封閉多邊形結(jié)構(gòu)描述復(fù)雜地質(zhì)結(jié)構(gòu)的拓?fù)潢P(guān)系，從而解決正演模擬建模中存在的問(wèn)題。

1.1 方法分析：

經(jīng)典的Delaunay三角剖分算法主要有兩類(lèi)：①Watson算法;②局部變換法。Watson算法又稱(chēng)為Delaunay空洞算法或加點(diǎn)法，從一個(gè)三角形開(kāi)始，每次加一個(gè)點(diǎn)，保證每一步得到的當(dāng)前三角形是局部?jī)?yōu)化的。Delaunay三角剖分加點(diǎn)算法采用點(diǎn)定位方法[5]。點(diǎn)法利用Delaunay空洞性質(zhì)，簡(jiǎn)明地實(shí)現(xiàn)了三角剖分。這種方法的優(yōu)點(diǎn)是在實(shí)現(xiàn)上比局部變換算法相對(duì)容易，而且與空間的維度無(wú)關(guān)。該算法在處理新點(diǎn)加入時(shí)，會(huì)重新計(jì)算三角形單元并判斷其屬性。如果包含新點(diǎn)的三角形單元不再符合Delaunay屬性，則這些三角形單元被刪除，形成Delaunay空洞，然后算法將新點(diǎn)與組成空洞的每一個(gè)頂點(diǎn)相連生成一個(gè)新邊，根據(jù)空球?qū)傩钥梢宰C明這些新邊都是局部Delaunay的，因此新生成的三角網(wǎng)格仍是Delaunay的[6]。

考慮到透鏡體的網(wǎng)格化問(wèn)題,這里采用了加點(diǎn)法的Constraint Delaunay[7]算法，利用其空?qǐng)A特性和最大化最小角特性，先將模型的原始點(diǎn)集和邊進(jìn)行三角剖分，再結(jié)合文獻(xiàn)[1]提出的二維復(fù)雜地質(zhì)結(jié)構(gòu)，利用原始(非添加)點(diǎn)集和線段集尋找封閉區(qū)域，從而形成系列的多邊形來(lái)描述模型。

1.2 模型的三角剖分

Constraint Delaunay 要求處理的對(duì)象必須是PLSG(Planar straight line graph)，PLSG定義為點(diǎn)和線段組成的集合，要求每個(gè)線段的端點(diǎn)，應(yīng)該在點(diǎn)的序列中，而端點(diǎn)由點(diǎn)的序號(hào)來(lái)標(biāo)示。借鑒文獻(xiàn)[1]的思路和Constraint Delaunay三角剖分，只要在前期處理中把層位、斷層、透鏡體等對(duì)象分解為點(diǎn)和線段集合就可以滿(mǎn)足三角剖分的輸入要求，而對(duì)于每次交互式修改，就相當(dāng)于對(duì)離散點(diǎn)集和線段集合進(jìn)行重構(gòu)，從而克服了層狀結(jié)構(gòu)模型的限制約束。對(duì)模型進(jìn)行Constraint Delaunay 三角剖分的步驟：

1)對(duì)原始輸入數(shù)據(jù)處理，形成滿(mǎn)足PLSG 圖的點(diǎn)和線段的集合，采用加點(diǎn)法進(jìn)行Delaunay 三角剖分。

2)恢復(fù)丟失的線段：用Constriant Delaunay 三角剖分，刪除線段有重疊的三角形，在線段兩側(cè)重新進(jìn)行三角剖分，以此確保在沒(méi)有引進(jìn)任何新點(diǎn)的情況下線段的恢復(fù)。

3)刪除不必要的三角形。主要是刪除“空洞”和“凹”處的三角形。三角剖分的計(jì)算量很大，這種處理有利于提高三角剖分的效率。

Lawson 算法還有一個(gè)步驟是保持三角剖分的Delaunay 性質(zhì)。它通過(guò)插入新點(diǎn)，并對(duì)三角剖分進(jìn)行細(xì)化，直到對(duì)三角形最小角和最大面積的限制條件滿(mǎn)足[8]。保持三角剖分的Delaunay 性質(zhì)本來(lái)是Delaunay三角剖分的一個(gè)核心，但本方法的目的在于形成封閉多邊形結(jié)構(gòu)，對(duì)于多邊形的形狀要求盡可能保持模型本身的形狀，對(duì)于保持三角剖分的Delaunay特性沒(méi)有特殊要求，因而在本方法的剖分中沒(méi)有應(yīng)用。同樣，步驟2)也是為了模型的準(zhǔn)確性以及減少三角化的計(jì)算量，采用不增加新點(diǎn)的方法，從而簡(jiǎn)化后面的封閉拓?fù)錁?gòu)成的工作。

1.3 形成封閉拓?fù)浣Y(jié)構(gòu)

模型三角剖分結(jié)果的結(jié)果是一系列的三角形，對(duì)于建模需要得到封閉體還需要基于三角剖分的結(jié)果(三角形集合)進(jìn)行搜索，形成由原始邊(線段)構(gòu)成的封閉多邊形，其步驟如下：

1)三角剖分后形成若干個(gè)三角形，標(biāo)記不是原拾取點(diǎn)形成的邊以及原拾取點(diǎn)所形成的邊集合。

2)確定一個(gè)起始三角形，預(yù)定義一個(gè)多邊形。

3) 將三角形的原拾取點(diǎn)和邊添加到定義的多邊形，并依據(jù)非原始邊去找新的三角形，依次循環(huán)遞歸查找，直到無(wú)法找到非原始點(diǎn)的邊。

圖1 封閉拓?fù)浣Y(jié)構(gòu)流程Fig.1 Flow chart of closed topology

圖2 Marmousi模型剖分及封閉結(jié)果Fig.2 Results of Marmousi model triangulation

圖3 Marmousi原始速度模型Fig.3 Original Marmousi vp-model

4)利用上步記錄的原始點(diǎn)邊集合添加到多邊形，可組成一個(gè)封閉多邊形(未形成閉合就不紀(jì)錄此多邊形)。

5)查找沒(méi)有用到過(guò)的另一個(gè)三角形，重復(fù)步驟3)、步驟4)，直到找完全部三角形。

圖4 Marmousi模型網(wǎng)格化結(jié)果Fig.4 Results of gridding Marmousi model

圖5 網(wǎng)格模型正演結(jié)果Fig.5 One shot record of gridding Marmousi model

形成封閉拓?fù)浣Y(jié)構(gòu)流程圖(圖1)。

通過(guò)以上流程可以得到一系列的封閉多邊形，在此基礎(chǔ)上根據(jù)網(wǎng)格劃分的步長(zhǎng)，從道方向?qū)W(wǎng)格點(diǎn)的坐標(biāo)進(jìn)行遍歷，從而確定網(wǎng)格點(diǎn)位于屬于哪個(gè)封閉多邊形，這樣就可以進(jìn)行屬性填充，完成模型的網(wǎng)格化從而得到正演所需的網(wǎng)格模型。

2 應(yīng)用實(shí)例

利用本文方法，結(jié)合Qt.4.X(X>2)開(kāi)發(fā)工具進(jìn)行開(kāi)發(fā)，形成了完整的實(shí)用軟件，以下是具體實(shí)例。

2.1 Marmousi模型

Marmousi模型是地球物理經(jīng)典模型，涵蓋了大多數(shù)地質(zhì)模型對(duì)象，圖2給出了三角剖分結(jié)果、封閉結(jié)構(gòu)結(jié)果(塊狀模型)。

圖3是Marmousi模型的原始縱波速度模型，圖4是由剖分結(jié)果按上述網(wǎng)格化方法填充屬性后得到的網(wǎng)格化模型，對(duì)比圖2與圖3反映了剖分的準(zhǔn)確性。通過(guò)圖3和圖4對(duì)比，反映了網(wǎng)格化結(jié)果滿(mǎn)足建模的需求。

在模型中左邊的3個(gè)狹長(zhǎng)的復(fù)雜的透鏡體，三角剖分時(shí)都被狹長(zhǎng)的三角網(wǎng)覆蓋，三角網(wǎng)內(nèi)多邊形保持了Delaunay特性反應(yīng)了剖分算法的正確性，網(wǎng)格化結(jié)果也反映了尋找封閉區(qū)域算法的準(zhǔn)確性。

圖6 實(shí)際模型剖分及封閉結(jié)果Fig.6 Results of real model triangulation

圖7 實(shí)際模型網(wǎng)格化結(jié)果Fig.7 Results of gridding real model

由網(wǎng)格化的數(shù)據(jù)用iSeisWave(自研軟件)軟件進(jìn)行正演模擬，激發(fā)點(diǎn)位于模型2 100 m處，最小炮檢距為1 050 m，道間距為25 m，得到的單炮記錄如圖5所示，可見(jiàn)本方法能滿(mǎn)足地震正演模擬的建模需求。

2.2 惠民凹陷盤(pán)河模型

此模型是勘探生產(chǎn)中的實(shí)際復(fù)雜模型，三角剖分結(jié)果和網(wǎng)格化結(jié)果如圖6、圖7所示。

圖6、圖7兩個(gè)實(shí)例中都存在狹長(zhǎng)三角形，就是因?yàn)槿瞧史种惺÷粤瞬襟E4)-用Lawson 算法保持三角剖分的Delaunay 性質(zhì)。圖7中，中間及右邊的大斷層起點(diǎn)并沒(méi)有與其他層位或斷層相交，三角剖分正常，網(wǎng)格化時(shí)都因非閉合而被忽略，再次驗(yàn)證了形成封閉體算法的正確性。

通過(guò)實(shí)踐，每次交互修改后形成PLSG，可以對(duì)模型進(jìn)行重構(gòu)——直接進(jìn)行三角剖分，三角剖分及形成封閉多邊形在三角網(wǎng)較大數(shù)量級(jí)(10×104級(jí)別)基本可以做到實(shí)時(shí)，對(duì)交互式建模提供了極大的便利(無(wú)需考慮層位和逆斷層的控制問(wèn)題)。同時(shí)在網(wǎng)格化步長(zhǎng)較大(橫向和縱向步長(zhǎng)都為10 m)的情況下，經(jīng)典模型和實(shí)際模型的網(wǎng)格化結(jié)果都達(dá)到了生成要求。

3 結(jié)論

筆者采用的基于三角網(wǎng)的建模方法突破了層狀建模的局限性，同時(shí)比已有的塊狀建模方法更加簡(jiǎn)單，效率高。本方法不僅能建立塊狀模型，還能得到非結(jié)構(gòu)化網(wǎng)格模型，因此有更好的應(yīng)用前景，為地震勘探中復(fù)雜交互式建模提供了有效的方法。