楊 強(qiáng)

(中國石化 石油物探技術(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)模型。它要求每一個分界面都必須從模型體的左邊界貫穿到模型體的右邊界，分界面按順序由上到下依序排列，不得交叉。對于逆斷層、尖滅、透鏡體等復(fù)雜地層情況，只能人為地簡化模型，從而滿足層狀模型。

在實(shí)際的地球物理勘探中，層狀結(jié)構(gòu)模型有兩個嚴(yán)重不足:①在采集、處理、解釋各階段的模型大都有逆斷層、尖滅、透鏡體等復(fù)雜地質(zhì)元素,層狀結(jié)構(gòu)模型無法準(zhǔn)確描述復(fù)雜的地質(zhì)拓?fù)浣Y(jié)構(gòu)，從而無法得到網(wǎng)格化模型;②處理中的疊前深度偏移速度建模以及地震解釋構(gòu)造建模等，都要求能夠交互修改模型，交互編輯中的層位修改及移動必須遵循層狀規(guī)則，這對交互式復(fù)雜建模而言很困難，也不便利。蔣先藝[1]提出用點(diǎn)、段、線、面的概念描述二維封閉結(jié)構(gòu)模型，實(shí)現(xiàn)了二維復(fù)雜地質(zhì)結(jié)構(gòu)，其優(yōu)點(diǎn)是地震建模的數(shù)據(jù)層位、斷層、透鏡體都可以用點(diǎn)、線、面集來描述。這種方法比較復(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)格生成算法是第一個理論上保證網(wǎng)格劃分算法在實(shí)踐中真正令人滿意的算法[3]。Refine Delaunay 三角化方法解決了保邊界和內(nèi)嵌邊界的問題,該方法往往通過在保留邊(約束邊)上加入新的節(jié)點(diǎn)以實(shí)現(xiàn)保邊界的目的[4]。筆者通過對Delaunay三角剖分的研究，以點(diǎn)和線段為基礎(chǔ)，對層位、斷層、透鏡體進(jìn)行三角剖分，采用封閉多邊形結(jié)構(gòu)描述復(fù)雜地質(zhì)結(jié)構(gòu)的拓?fù)潢P(guān)系，從而解決正演模擬建模中存在的問題。

1.1 方法分析：

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

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

1.2 模型的三角剖分

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

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

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

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

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

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

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

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

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

3) 將三角形的原拾取點(diǎn)和邊添加到定義的多邊形，并依據(jù)非原始邊去找新的三角形，依次循環(huán)遞歸查找，直到無法找到非原始點(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)邊集合添加到多邊形，可組成一個封閉多邊形(未形成閉合就不紀(jì)錄此多邊形)。

5)查找沒有用到過的另一個三角形，重復(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)。

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

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

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

2.1 Marmousi模型

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

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

在模型中左邊的3個狹長的復(fù)雜的透鏡體，三角剖分時都被狹長的三角網(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所示，可見本方法能滿足地震正演模擬的建模需求。

2.2 惠民凹陷盤河模型

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

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

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

3 結(jié)論

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