一種基于鄰域作用力的三維網(wǎng)格平滑方法

王海起, 彭佳琦, 王 藝, 董倩楠, 車 磊, 龔蔚青

(1.中國(guó)石油大學(xué)地球科學(xué)與技術(shù)學(xué)院,山東青島 266580; 2.中國(guó)石化石油勘探開發(fā)研究院,北京 100083)

一種基于鄰域作用力的三維網(wǎng)格平滑方法

王海起1, 彭佳琦1, 王 藝1, 董倩楠1, 車 磊1, 龔蔚青2

(1.中國(guó)石油大學(xué)地球科學(xué)與技術(shù)學(xué)院,山東青島 266580; 2.中國(guó)石化石油勘探開發(fā)研究院,北京 100083)

針對(duì)以不規(guī)則六面體為基元的三維網(wǎng)格模型在對(duì)三維區(qū)域邊界表達(dá)時(shí)存在鋸齒現(xiàn)象,在保持原有網(wǎng)格拓?fù)浣Y(jié)構(gòu)不變的情況下,提出一種基于離散屬性數(shù)據(jù)的鄰域網(wǎng)格相互作用力平滑方法。受磁鐵“同性相斥”物理特性的啟發(fā),將三維網(wǎng)格模型共用頂點(diǎn)的鄰域網(wǎng)格體元素根據(jù)屬性值劃分為不同的體塊,體塊之間具有相斥的作用力,根據(jù)體塊對(duì)頂點(diǎn)的作用合力方向和大小確定網(wǎng)格頂點(diǎn)平滑移動(dòng)的方向和距離。實(shí)際應(yīng)用表明平滑存在著“拐點(diǎn)”現(xiàn)象,平滑系數(shù)取值接近拐點(diǎn)時(shí)平滑效果較好,取值偏大或偏小評(píng)價(jià)指標(biāo)性能下降。

六面體網(wǎng)格; 鋸齒效應(yīng); 三維平滑; 鄰域網(wǎng)格; 作用力

真三維地質(zhì)體建模的模型體元通常采用六面體、四面體、棱柱體等[1],其中,六面體網(wǎng)格具有組織形式規(guī)則、網(wǎng)格單元數(shù)量及重劃分次數(shù)較少、計(jì)算效率高等優(yōu)點(diǎn),在三維地質(zhì)體、油藏體建模中應(yīng)用廣泛,并被大多數(shù)商業(yè)油藏建模和數(shù)值模擬軟件所支持[2-4]。然而,這種源自矩形網(wǎng)格的體模型在對(duì)三維區(qū)域邊界進(jìn)行離散化網(wǎng)格表達(dá)時(shí)存在鋸齒效應(yīng)。高精度的精細(xì)地質(zhì)模型可以降低鋸齒效應(yīng)的影響、提高表征精度,但即使對(duì)于一個(gè)中等規(guī)模的油藏體,精細(xì)的地質(zhì)表征往往具有百萬(wàn)級(jí)甚至千萬(wàn)級(jí)的網(wǎng)格規(guī)模,受限于實(shí)際實(shí)現(xiàn)技術(shù),精細(xì)地質(zhì)模型和油藏?cái)?shù)值模擬器可以支持的網(wǎng)格規(guī)模之間存在著較大差距,對(duì)三維網(wǎng)格進(jìn)行粗化就是將精細(xì)地質(zhì)模型粗化到油藏?cái)?shù)值模擬器能夠接受的網(wǎng)格規(guī)模[5]。然而粗化后的網(wǎng)格六面體體元較大,對(duì)圈閉、油層等區(qū)域邊界的表達(dá)由于鋸齒效應(yīng)的存在不符合實(shí)際的地質(zhì)現(xiàn)象,影響模型的描述精度、可視化效果及后期模擬計(jì)算的準(zhǔn)確性。對(duì)于六面體模型邊界鋸齒現(xiàn)象,可以采用邊界網(wǎng)格重構(gòu)或加密方法解決[6],這種方式可以構(gòu)造出精細(xì)的邊緣形狀,達(dá)到了精細(xì)的平滑效果,但對(duì)于油藏模型組織結(jié)構(gòu)來(lái)說(shuō),這種方法打亂了原有結(jié)構(gòu)的規(guī)則行列網(wǎng)格檢索方式,為查找和數(shù)值模擬計(jì)算增加了難度。另一種解決途徑通過(guò)適當(dāng)移動(dòng)六面體的頂點(diǎn)進(jìn)行平滑,這種方式只能進(jìn)行粗略的平滑,難以刻畫精細(xì)的邊緣形狀,但保留了原有的網(wǎng)格組織形式,網(wǎng)格數(shù)和頂點(diǎn)數(shù)均沒有發(fā)生改變。針對(duì)鋸齒平滑問題,筆者提出一種鄰域網(wǎng)格作用力平滑算法,考慮到實(shí)際應(yīng)用中油藏模型并不需要特別精細(xì)地刻畫邊緣,同時(shí)要保持網(wǎng)格原有的三維拓?fù)浣Y(jié)構(gòu),平滑采用網(wǎng)格頂點(diǎn)位置自適應(yīng)優(yōu)化方式,根據(jù)不同屬性值體塊之間的相互作用力,確定頂點(diǎn)平滑移動(dòng)的方向和距離。

1 三維平滑算法

對(duì)于網(wǎng)格平滑,經(jīng)典的平滑算法包括Laplacian平滑、Taubin平滑、平均曲率法等。Laplacian平滑的核心是將網(wǎng)格內(nèi)部節(jié)點(diǎn)的位置移動(dòng)到與該節(jié)點(diǎn)共面節(jié)點(diǎn)組成的多面體的體心處,算法簡(jiǎn)單[7]。Taubin算法在Laplacian算法的基礎(chǔ)上引入了濾波器及權(quán)系數(shù),可抑制拉普拉斯算子引起的變形收縮[8]。平均曲率法則遵循曲面曲率變化均勻即為光滑的原則[9]。上述方法各有優(yōu)缺點(diǎn)[9],不少學(xué)者亦在此基礎(chǔ)上進(jìn)行了大量的優(yōu)化改進(jìn),但是對(duì)于三維油藏網(wǎng)格模型,這些平滑方法在應(yīng)用時(shí)會(huì)產(chǎn)生一些問題。

三維油藏網(wǎng)格模型不同于其他建筑或工具的三維模型。圖1所示為某工程零件的三維六面體網(wǎng)格模型,這類模型所有網(wǎng)格屬性相同,只需對(duì)表面因形狀產(chǎn)生的鋸齒進(jìn)行平滑。油藏網(wǎng)格模型表面并不需要平滑,模型內(nèi)部因?yàn)榫W(wǎng)格屬性值不同產(chǎn)生的鋸齒現(xiàn)象需要平滑,若分別將相同屬性值的網(wǎng)格提取出來(lái)對(duì)表面進(jìn)行平滑,平滑后的拓?fù)潢P(guān)系可能被破壞。

綜上,如何能借助六面體網(wǎng)格模型三維點(diǎn)陣的查找優(yōu)勢(shì),避免破壞原有的網(wǎng)格結(jié)構(gòu),又能保證網(wǎng)格之間的三維拓?fù)潢P(guān)系,是三維油藏模型平滑的關(guān)鍵問題。

圖1 六面體網(wǎng)格邊界重構(gòu)加密Fig.1 Boundary reconstruction and densification of hexahedral meshes

2 鄰域網(wǎng)格作用力平滑算法

在三維網(wǎng)格平滑中,對(duì)上述平滑思路進(jìn)行擴(kuò)展,以頂點(diǎn)周圍共用該頂點(diǎn)的8個(gè)鄰域網(wǎng)格體組成平滑窗口,將網(wǎng)格體劃分為不同的體塊,不同體塊對(duì)頂點(diǎn)相斥作用力的大小度量了頂點(diǎn)平滑移動(dòng)的不同權(quán)重,所有體塊作用力的合力決定了頂點(diǎn)最終移動(dòng)的距離和方向。

以圖2為例闡述算法的主要思路,設(shè)頂點(diǎn)v周圍8個(gè)網(wǎng)格體的離散屬性值代表不同的沉積微相類型,如用1、2、3分別代表壩主體、壩緣、前三角洲類型,相同的屬性值指示相同的沉積微相類型,據(jù)此可將8個(gè)鄰域網(wǎng)格體劃分為3個(gè)體塊B1、B2、B3,分別代表上述3種微相類型。為了便于展示,圖2中將頂點(diǎn)v與鄰域網(wǎng)格單元拉開了距離。

圖2 鄰域作用力示意圖Fig.2 Schematic of neighborhood forces

將每個(gè)體塊看作磁鐵的“同極”,各同極對(duì)頂點(diǎn)v具有相斥的作用力,作用力采用基于重力模型改進(jìn)的空間相互作用力SIM模型(spatial interaction model)進(jìn)行度量[10-11]。設(shè)體塊Bk的質(zhì)量為Mk,頂點(diǎn)v質(zhì)量為Mv,體塊Bk對(duì)頂點(diǎn)v的相斥作用力向量Fk大小為

(1)

式中,G、α分別為引力常數(shù)和質(zhì)量指數(shù),這里可取1;Vk為體塊Bk中心與頂點(diǎn)v連線向量,體塊中心可取幾何中心或質(zhì)心,|Vk|為向量模,表示體塊Bk與頂點(diǎn)v的距離;β為距離指數(shù),亦可取1。

Fk方向由Vk確定,圖2中3個(gè)體塊對(duì)頂點(diǎn)v的作用力F1、F2、F3大小分別由式(1)確定,方向分別由向量V1、V2、V3確定,則合力Fsum實(shí)際上是3個(gè)作用力向量之和。

對(duì)于每個(gè)Fk,取G=α=β=1,由運(yùn)動(dòng)定律F=ma可知頂點(diǎn)v的加速度向量ak大小為

(2)

ak方向同樣由Vk確定。借鑒萬(wàn)有引力搜索算法GSA(gravitational search algorithm)的處理方式[12],在不考慮初始速度和時(shí)間的前提下,加速度ak代表了頂點(diǎn)v在作用力Fk作用下的移動(dòng)偏移量Δdk,則在合力Fsum作用下總的移動(dòng)偏移量為各向量Δdk之和。

上述公式亦說(shuō)明Δdk大小與體塊Bk的質(zhì)量Mk成正比,與距離|Vk|成反比。將每個(gè)網(wǎng)格單元的質(zhì)量簡(jiǎn)單地視為1,則Δdk大小實(shí)際上與體塊Bk包含的網(wǎng)格單元個(gè)數(shù)nk成正比,表明當(dāng)體塊包含的網(wǎng)格單元較多時(shí),對(duì)頂點(diǎn)的作用力較大,頂點(diǎn)移動(dòng)距離較長(zhǎng),如圖2中的體塊B3;反之,則對(duì)頂點(diǎn)的作用力較小,頂點(diǎn)移動(dòng)距離較短,如體塊B2。

對(duì)三維區(qū)域邊界的鋸齒進(jìn)行平滑處理,目的是降低邊界表達(dá)的噪聲,使其符合實(shí)際的地質(zhì)現(xiàn)象。對(duì)于有些頂點(diǎn),例如斷層、尖滅等,是不需要進(jìn)行平滑處理的。當(dāng)頂點(diǎn)需進(jìn)行平滑處理時(shí),其鄰域的各體塊總是將頂點(diǎn)外推,目的是使體塊尖銳、鋸齒狀的邊界面和邊界點(diǎn)趨向光滑,這與地層分界面的連續(xù)性是一致的,只不過(guò)質(zhì)量較大的體塊將頂點(diǎn)外推的作用力較大,質(zhì)量較小的體塊將頂點(diǎn)外推的作用力較小,不同方向、大小的作用力形成合力最終將頂點(diǎn)推到相對(duì)平衡的位置,既降低了邊界的鋸齒效應(yīng),又保持了網(wǎng)格的穩(wěn)定性。

綜上,算法的主要思路可概括為:將共用同一頂點(diǎn)的8個(gè)鄰域網(wǎng)格體劃分為不同的體塊,每個(gè)體塊包含具有相同屬性類型的網(wǎng)格。將每個(gè)體塊看作磁鐵的“同極”,同極之間具有相斥的作用力,頂點(diǎn)移動(dòng)的方向?yàn)轶w塊作用力的方向,移動(dòng)的距離以作用力的大小作為權(quán)重,作用力大則移動(dòng)距離大,作用力小則移動(dòng)距離小。所有體塊合力的方向和大小決定了頂點(diǎn)最終移動(dòng)的方向和距離。

通過(guò)該方法以頂點(diǎn)鄰域八網(wǎng)格體為平滑窗口,依次遍歷各頂點(diǎn),若某頂點(diǎn)鄰域八網(wǎng)格體具有相同屬性值,則認(rèn)為該頂點(diǎn)鄰域區(qū)域具有均質(zhì)性,不需要平滑。

下面對(duì)網(wǎng)格體塊分布的幾種具體情況進(jìn)行圖示分析:

(1)頂點(diǎn)鄰域八網(wǎng)格存在兩種不同屬性值,如圖3(a)所示,黃色網(wǎng)格個(gè)數(shù)為1,綠色網(wǎng)格個(gè)數(shù)為7,左圖為原始結(jié)構(gòu),中圖黃綠箭頭為兩種力的作用方向,紅色箭頭為合力的作用方向,右圖為頂點(diǎn)沿合力方向平移后的效果。

(2)頂點(diǎn)鄰域八網(wǎng)格存在4個(gè)黃色網(wǎng)格和4個(gè)綠色網(wǎng)格,如圖3(b)所示,上下兩側(cè)的作用力相互抵消,不需要平滑。

(3)頂點(diǎn)鄰域八網(wǎng)格存在3種不同屬性值,作用力效果與平滑效果如圖3(c)所示。

(4)如圖3(d)所示的5×5網(wǎng)格中心存在一凸起網(wǎng)格,左圖紅色箭頭示意各頂點(diǎn)的移動(dòng)方向,右圖為頂點(diǎn)移動(dòng)平滑后的效果。

(5)在剖面情況下,圖4左圖為處理前的網(wǎng)格,紅色箭頭示意相應(yīng)頂點(diǎn)移動(dòng)的方向,右圖為頂點(diǎn)移動(dòng)后的網(wǎng)格,可以看出,該方法同樣也適用于二維網(wǎng)格的平滑。

對(duì)于共用頂點(diǎn)的八網(wǎng)格,若體塊包含的網(wǎng)格數(shù)較少,例如只有1個(gè)網(wǎng)格,如圖3(a)所示的黃色網(wǎng)格,則意味著該網(wǎng)格是孤立的,是一個(gè)凸起,需要平滑凸起邊界;反過(guò)來(lái)看,綠色網(wǎng)格的數(shù)量較多,可以看作成包圍趨勢(shì),則意味著綠色網(wǎng)格存在一個(gè)凹陷,需要平滑凹陷邊界。在本文提出的算法中,借助“同性相斥”的作用力思想,給網(wǎng)格數(shù)量較多的體塊一個(gè)較大的力,給網(wǎng)格數(shù)量較少的體塊一個(gè)較小的力,從而合力的方向和大小決定了頂點(diǎn)移動(dòng)的方向和距離,這種方法的結(jié)果也符合前述的平滑規(guī)律。

對(duì)于具有更多屬性值的網(wǎng)格、或者屬性值交錯(cuò)散亂分布的情況,如圖5所示,每個(gè)網(wǎng)格幾乎都相當(dāng)于一個(gè)異變(凸起或凹陷),這種情況一般不需要進(jìn)行平滑處理,應(yīng)保持其原有的特征。對(duì)應(yīng)到本文的方法,若屬性分布散亂,則會(huì)產(chǎn)生四面八方的力,這些力互相抵消或完全抵消,則頂點(diǎn)僅需做細(xì)微的移動(dòng)、或不移動(dòng),因此在復(fù)雜情況下,本文的方法也符合一般的平滑規(guī)律。實(shí)際上,在油藏模型中,地質(zhì)屬性的分布是有一定規(guī)律的,屬性在一個(gè)頂點(diǎn)周圍散亂分布的情況較為少見。

圖3 不同體塊的作用力平滑F(xiàn)ig.3 Smoothing based on forces of different blocks

圖4 剖面情況下兩種屬性的平滑效果Fig.4 Smoothing effect of two attributes in a profile view

圖5 鄰域網(wǎng)格的其他體塊分布情況Fig.5 Distributions of other blocks for eight-neighbor meshes

3 平滑算法流程

在油藏模型中,有些不需要平滑的特殊點(diǎn),如圖6的示意圖中,模型邊界的角點(diǎn)、邊點(diǎn)不需要進(jìn)行平滑;帶有斷層特征的點(diǎn)不需要進(jìn)行平滑;另外一些具有特殊意義的特征點(diǎn),例如特殊底層的尖滅,應(yīng)在平滑處理前由地質(zhì)工程師進(jìn)行交互操作,圈定或輸入一些特征點(diǎn)或特征區(qū)域,避免平滑時(shí)被破壞。遍歷網(wǎng)格邏輯IJK頂點(diǎn)時(shí),判斷其是否為特征點(diǎn),若不是特征點(diǎn),則進(jìn)行平滑處理。

平滑算法的具體流程如下:

(1)根據(jù)網(wǎng)格IJK集合及頂點(diǎn)索引集合依次遍歷三維模型每個(gè)頂點(diǎn),及共用該頂點(diǎn)的鄰域網(wǎng)格;對(duì)于每個(gè)頂點(diǎn),進(jìn)行下面(2)～(9)步處理。

(2)設(shè)遍歷的當(dāng)前頂點(diǎn)v坐標(biāo)為(xv,yv,zv),則①若頂點(diǎn)為邊界面上的點(diǎn),統(tǒng)計(jì)頂點(diǎn)鄰域網(wǎng)格在邊界面上的屬性值分布情況,若屬性值不相同,在二維面上采用(3)～(9)步進(jìn)行頂點(diǎn)移動(dòng)平滑;②若頂點(diǎn)為三維模型的角點(diǎn)或邊點(diǎn),則不進(jìn)行平滑處理;③若頂點(diǎn)為被標(biāo)注過(guò)的特征點(diǎn),則不進(jìn)行平滑處理;④若頂點(diǎn)在模型內(nèi)部,且鄰域網(wǎng)格屬性值均相同(無(wú)效網(wǎng)格亦被看作一種網(wǎng)格屬性值),則認(rèn)為是均質(zhì)的,不需要處理;⑤若頂點(diǎn)在模型內(nèi)部,且鄰域網(wǎng)格屬性值不相同,則進(jìn)行(3)～(9)步處理。

圖6 三維模型表面頂點(diǎn)不同位置分布示意圖Fig.6 Distribution schematic of different positions of surface vertexes of 3D model

(3)計(jì)算當(dāng)前頂點(diǎn)各鄰域網(wǎng)格質(zhì)心坐標(biāo)。將每個(gè)六面體鄰域網(wǎng)格單元視為均質(zhì),設(shè)網(wǎng)格單元各頂點(diǎn)的坐標(biāo)為(xi,yi,zi),其中i=1,2,…,8,i為頂點(diǎn)序號(hào),網(wǎng)格單元質(zhì)心坐標(biāo)為(xc,yc,zc),質(zhì)心計(jì)算公式為

(3)

(4)計(jì)算體塊質(zhì)心坐標(biāo)。將屬性值相同(將無(wú)效網(wǎng)格看作同一類屬性值)的鄰域網(wǎng)格看作一個(gè)體塊,則頂點(diǎn)(xv,yv,zv)的8個(gè)鄰域網(wǎng)格單元可分為多個(gè)體塊{Bk}。

(4)

(5)計(jì)算體塊Bk對(duì)當(dāng)前頂點(diǎn)作用力的方向向量。利用體塊質(zhì)心與當(dāng)前頂點(diǎn)(xv,yv,zv)的連線確定方向向量Vk,計(jì)算公式為

Vk=(Xk-xv,Yk-yv,Zk-zv).

(5)

(6)計(jì)算體塊Bk對(duì)當(dāng)前頂點(diǎn)作用力的移動(dòng)距離。將體塊Bk包含的每個(gè)網(wǎng)格單元質(zhì)量視為1,體塊質(zhì)量即為體塊包含的網(wǎng)格單元個(gè)數(shù)nk。體塊質(zhì)量大,則對(duì)頂點(diǎn)的作用力大,頂點(diǎn)移動(dòng)距離長(zhǎng);反之,則作用力小,頂點(diǎn)移動(dòng)距離短。

對(duì)于六面體模型,任一頂點(diǎn)(xv,yv,zv)的鄰域網(wǎng)格單元個(gè)數(shù)最多為8個(gè),當(dāng)所有鄰域網(wǎng)格均屬于同一個(gè)體塊時(shí),則意味著鄰域網(wǎng)格屬性值均相同,鄰域網(wǎng)格單元是均質(zhì)的,頂點(diǎn)不需要進(jìn)行平滑處理,即前述第(2)步中④ 說(shuō)明的情況,因此體塊Bk的最大質(zhì)量為7,則nk/7可以作為相對(duì)于最大質(zhì)量作用的最長(zhǎng)移動(dòng)距離,體塊質(zhì)量為nk時(shí)的移動(dòng)距離權(quán)重。設(shè)六面體網(wǎng)格對(duì)角線長(zhǎng)度Ldiag為頂點(diǎn)移動(dòng)距離的上限,則頂點(diǎn)(xv,yv,zv)在體塊Bk作用下的移動(dòng)距離Δdk(Δdk為移動(dòng)偏移向量Δdk的數(shù)值)可表達(dá)為

(6)

式中,γ為平滑系數(shù),用于對(duì)最大移動(dòng)距離進(jìn)行約束。

(7)根據(jù)體塊Bk的方向向量Vk、移動(dòng)距離Δdk,則在體塊Bk作用下頂點(diǎn)(xv,yv,zv)的移動(dòng)距離向量為

(7)

其中

(8)所有體塊{Bk}合力作用下的移動(dòng)分量為

(8)

(9)

算法流程如圖7所示。

圖7 平滑算法流程Fig.7 Flow chart of smoothing algorithm

4 算法分析及應(yīng)用

4.1 平滑系數(shù)γ分析

當(dāng)鄰域網(wǎng)格劃分為體塊后,公式(8)中每個(gè)體塊的質(zhì)心位置(Xk,Yk,Zk)、方向向量Vk、包含的網(wǎng)格個(gè)數(shù)nk及網(wǎng)格對(duì)角線長(zhǎng)度Ldiag都是確定的,因此平滑系數(shù)γ的選擇對(duì)平滑效果至關(guān)重要。

對(duì)于三維地質(zhì)體平滑,當(dāng)遇到“凸起”或“凹陷”時(shí),需要對(duì)其平滑,而不是將其完全“削掉”或“填平”,故對(duì)頂點(diǎn)移動(dòng)時(shí),應(yīng)小于網(wǎng)格對(duì)角線長(zhǎng)度的一半。直觀上理解,如圖8所示的二維規(guī)則網(wǎng)格,當(dāng)移動(dòng)距離約為對(duì)角線長(zhǎng)度的1/4時(shí),恰好將斜邊平滑近似為一條直線,較為符合常理。

圖8 兩種屬性交界的斜邊平滑F(xiàn)ig.8 Hypotenuse smoothing of borderlines of two attributes

為了確定最佳平滑系數(shù),對(duì)γ不同取值的平滑結(jié)果進(jìn)行比較,評(píng)價(jià)指標(biāo)為平滑頂點(diǎn)與體塊分界面上所有相鄰頂點(diǎn)連線形成的夾角余弦值的平均值。

圖9為二維視角下兩個(gè)鄰域體塊分界線上的平滑頂點(diǎn)與處于分界線上的相鄰頂點(diǎn)形成的向量夾角及相應(yīng)的余弦值。

圖9 不同夾角對(duì)應(yīng)的余弦值Fig.9 Different angle and its cosine

從圖中可以看出,兩向量夾角越接近于180°,即夾角連線越接近于直線,余弦值越接近于-1,平滑效果越好。

三維視角下,對(duì)于分界面上的平滑頂點(diǎn),需要分別計(jì)算與XOY、XOZ、YOZ三個(gè)剖面上處于分界線上的相鄰頂點(diǎn)連線形成的向量夾角余弦值,取其平均值作為評(píng)價(jià)指標(biāo)評(píng)價(jià)整體的平滑效果,其值越接近于-1,平滑效果越好。

采用的測(cè)試網(wǎng)格為二維規(guī)則格網(wǎng),網(wǎng)格劃分為63×63,網(wǎng)格大小為1 m×1 m×1 m,平滑系數(shù)γ不同取值的平滑效果及評(píng)價(jià)指標(biāo)變化情況分別如圖10、圖11所示。

圖10 不同平滑系數(shù)γ的平滑效果Fig.10 Smoothing effect of different smoothing factor γ

分析平滑效果及變化曲線圖,平滑系數(shù)γ為0表示沒有進(jìn)行平滑,如圖10(a)所示;隨著γ值逐漸增大,平滑發(fā)生作用,評(píng)價(jià)指標(biāo)也趨向-1;當(dāng)平滑系數(shù)γ增大到某個(gè)值(本算例為0.275)時(shí),曲線到達(dá)拐點(diǎn)位置,評(píng)價(jià)指標(biāo)達(dá)到最為接近-1的值(本算例為-0.3013),此時(shí)平滑效果最好,如圖10(b)所示;之后,隨著γ取值繼續(xù)增大,平滑性能下降,從圖10(c)可以看出,邊界出現(xiàn)過(guò)度平滑現(xiàn)象。

對(duì)規(guī)則六面體格網(wǎng),采用其他算例進(jìn)行相同的計(jì)算,最佳平滑系數(shù)(曲線拐點(diǎn)位置)變化不大,γ基本在0.26～0.29變化。對(duì)比實(shí)際圖形可發(fā)現(xiàn),當(dāng)區(qū)域邊界出現(xiàn)較多的連續(xù)臺(tái)階狀鋸齒現(xiàn)象時(shí),如圖12(a)所示,拐點(diǎn)位置更接近于γ=0.26;當(dāng)區(qū)域邊界出現(xiàn)較多的直角轉(zhuǎn)彎現(xiàn)象時(shí),如圖12(b)所示,拐點(diǎn)位置更接近于γ=0.29。實(shí)際應(yīng)用中,對(duì)于規(guī)則六面體格網(wǎng),可針對(duì)實(shí)際模型進(jìn)行計(jì)算以確定最佳的平滑系數(shù)γ值,也可簡(jiǎn)單地取區(qū)間[0.26, 0.29]中值0.275為γ近似值。

對(duì)于不規(guī)則六面體網(wǎng)格,由于實(shí)際地層分界面和地質(zhì)模型的多樣性,上述的區(qū)間并不適用,需針對(duì)具體情況進(jìn)行計(jì)算以獲取最優(yōu)的γ值。但拐點(diǎn)現(xiàn)象是存在的,因?yàn)樵u(píng)價(jià)指標(biāo)為平滑頂點(diǎn)與相鄰頂點(diǎn)連線形成的夾角余弦值或余弦平均值,夾角總是在[0, 360°]之間變化,其余弦曲線總是具有拐點(diǎn)位置(圖9)。

圖11平滑系數(shù)和評(píng)價(jià)指標(biāo)的關(guān)系曲線Fig.11 Curve between smoothing factor and evaluation index

圖12 兩種鋸齒現(xiàn)象Fig.12 Two kinds of jaggies

4.2 算法應(yīng)用

采用C++及OSG(OpenSceneGraph)三維引擎,以自主研發(fā)的Simtools油藏模擬分析軟件為平臺(tái),對(duì)基于鄰域網(wǎng)格作用力的三維網(wǎng)格平滑方法進(jìn)行了實(shí)現(xiàn),并加載了實(shí)際的三角洲沉積相數(shù)據(jù),加載數(shù)據(jù)的網(wǎng)格規(guī)模為97×105×8,網(wǎng)格平均大小為200 m×200 m×30 m。圖13～15分別為平滑前后俯視圖效果對(duì)比、俯視圖局部細(xì)節(jié)對(duì)比及剖面效果對(duì)比,可以看出,通過(guò)對(duì)沉積相邊界的平滑處理,模型的可視化效果有了較大的改善。

圖13 平滑前后俯視圖效果對(duì)比Fig.13 Effect contrast of top view before and after smoothing

圖14 平滑前后局部細(xì)節(jié)對(duì)比Fig.14 Local details contrast of top view before and after smoothing

圖15 平滑前后三維模型部分剖面對(duì)比Fig.15 Part profile contrast of 3D model before and after smoothing

5 結(jié)束語(yǔ)

基于鄰域網(wǎng)格作用力的三維網(wǎng)格體平滑方法，充分利用六面體網(wǎng)格數(shù)據(jù)組織形式,以IJK對(duì)頂點(diǎn)及其鄰域網(wǎng)格進(jìn)行遍歷,簡(jiǎn)單、快速。算法對(duì)需平滑頂點(diǎn)的鄰域網(wǎng)格根據(jù)屬性值劃分為不同的體塊,根據(jù)體塊對(duì)頂點(diǎn)的作用合力方向確定頂點(diǎn)平滑移動(dòng)的方向,通過(guò)移動(dòng)方向的鄰域網(wǎng)格對(duì)角線長(zhǎng)度限制頂點(diǎn)移動(dòng)的最大距離,使頂點(diǎn)平滑距離不超出鄰域網(wǎng)格的范圍,從而在平滑過(guò)程中完全保持了網(wǎng)格原有的拓?fù)潢P(guān)系。本方法基于“互斥作用力”的平滑思路不僅可以用于六面體網(wǎng)格的二、三維平滑,也可以推廣到其他類型基元構(gòu)成的三維模型,如四面體網(wǎng)格、三棱柱網(wǎng)格等。三維地質(zhì)體、油藏體模型使用的角點(diǎn)網(wǎng)格需對(duì)尖滅、缺失、斷層等地質(zhì)現(xiàn)象建模,導(dǎo)致網(wǎng)格模型存在著無(wú)效網(wǎng)格或鄰域網(wǎng)格不共點(diǎn)等復(fù)雜情況。下一步工作將具體分析特殊情況下鄰域網(wǎng)格的查找方式,對(duì)三維網(wǎng)格體平滑方法進(jìn)行完善和優(yōu)化。

[1] 武強(qiáng),徐華.三維地質(zhì)建模與可視化方法研究[J].中國(guó)科學(xué)(D輯):地球科學(xué),2004,34(1):54-60.

WU Qiang, XU Hua. On three-dimensional geological modeling and visualization[J]. Science in China(Ser D): Earth Sciences, 2004,47(8):739-748.

[2] 趙永軍,舒曉,胡勇,等.一種復(fù)雜曲流帶儲(chǔ)層三維構(gòu)型建模新方法[J].中國(guó)石油大學(xué)學(xué)報(bào)(自然科學(xué)版),2015,39(1):1-7.

ZHAO Yongjun, SHU Xiao, HU Yong, et al. A new 3D architecture modeling method of complex meander belt reservoir[J]. Journal of China University of Petroleum (Edition of Natural Science), 2015,39(1):1-7.

[3] 林承焰,董春梅,任麗華,等.油藏描述技術(shù)發(fā)展及啟示[J].中國(guó)石油大學(xué)學(xué)報(bào)(自然科學(xué)版),2013,37(5):22-27.

LIN Chengyan, DONG Chunmei, REN Lihua, et al. Development of reservoir characterization and its stimulation[J]. Journal of China University of Petroleum (Edition of Natural Science), 2013,37(5):22-27.

[4] 毛小平,張志庭,錢真.用角點(diǎn)網(wǎng)格模型表達(dá)地質(zhì)模型的剖析及在油氣成藏過(guò)程模擬中的應(yīng)用[J].地質(zhì)學(xué)刊,2012,36(3):265-273.

MAO Xiaoping, ZHANG Zhiting, QIAN Zhen. Analysis of geological model expressed by corner point grid model and application in process of simulation of hydrocarbon accumulation process[J]. Journal of Geology, 2012,36(3):265-273.

[5] 劉慧卿,馬代鑫,張松亭.精細(xì)油藏描述數(shù)據(jù)的粗化方法研究[J].石油大學(xué)學(xué)報(bào)(自然科學(xué)版),2003,27(6):36-38.

LIU Huiqing, MA Daixin, ZHANG Songting. Up-scaling technique for handling data of fine reservoir description[J]. Journal of the University of Petroleum, China (Edition of Natural Science),2003,27(6):36-38.

[6] KWAK D Y, IM Y T. Remeshing for metal forming simulations-Part II: three-dimensional hexahedral mesh generation[J]. International Journal for Numerical Methods in Engineering, 2002,53(11):2501-2528.

[7] HENGL T. A practical guide to geostatistical mapping[M].Amsterdam: University of Amsterdam, 2009:117-128.

[8] 李光明,田捷,何暉光,等.基于距離均衡化的網(wǎng)格平滑算法[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),2002,14(9):820-823.

LI Guangming, TIAN Jie, HE Huiguang, et al. A mesh smoothing algorithm based on distance equalization[J]. Journal of Computer-Aided Design & Computer Graphics, 2002,14(9):820-823.

[9] 佟一飛.基于頂點(diǎn)擴(kuò)散的三維網(wǎng)格平滑技術(shù)研究[D].大連:遼寧師范大學(xué),2008.

TONG Yifei. Research on the smoothing techniques for 3-D mesh based on vertex diffusion[D]. Dalian: Liaoning Normal University, 2008.

[10] KIM M, LEE J. A data transformation method for visualizing the statistical information based on the grid[J]. Journal of Korea Spatial Information Society, 2015,23(5):31-40.

[11] KORDI M, FOTHERINGHAM A S. Spatially weighted interaction models (SWIM)[J]. Annals of the American Association of Geographers, 2016,106(5):990-1012.

[12] RASHEDI E, NEZAMABADI-POUR H, SARYAZDI S. GSA: a gravitational search algorithm[J]. Information Sciences, 2009,179(13):2232-2248.

A3Dmeshessmoothingmethodbasedonneighborhoodforces

WANG Haiqi1, PENG Jiaqi1, WANG Yi1, DONG Qiannan1, CHE Lei1, GONG Weiqing2

(1.SchoolofGeosciencesinChinaUniversityofPetroleum,Qingdao266580,China;2.SINOPECResearchInstituteofPetroleumExploration&Development,Beijing100083,China)

Hexahedral meshes are widely used to build three-dimensional (3D) geological and reservoir models. However, this voxel model, derived from rectangular grid, has jagged effects when modeling the boundary of a 3D region, which is not consistent with actual geological phenomena. To address the jagged-edge phenomena, under the condition of keeping original topology structure invariable, a smoothing algorithm is proposed based on the forces of neighborhood meshes. The algorithm is inspired by "poles repel" derived from physical properties of a magnet, and the neighboring meshes of a smoothed vertex are divided into different blocks according to their attributes. The blocks have repulsive force. The direction and magnitude of resultant force determine the vertex direction and distance of smoothing movement. The practical application shows that there is an inflection point for smoothing performance and the smoothing performs better when the smoothing coefficient closes to the inflection point; or conversely, the performance is decreasing when smoothing coefficient is away from the inflection point.

hexahedral mesh; jagged effect; 3D smoothing; neighborhood mesh; force

2017-02-02

國(guó)家自然科學(xué)基金項(xiàng)目(41471322); 山東省自然科學(xué)基金項(xiàng)目(ZR2012DM010)

王海起(1972-), 男, 副教授, 博士, 研究方向?yàn)榈乩硇畔⑾到y(tǒng)與數(shù)字油田。E-mail: wanghaiqi@upc.edu.cn。

1673-5005(2017)06-0080-08

10.3969/j.issn.1673-5005.2017.06.009

TP 391.9

A

王海起,彭佳琦,王藝,等.一種基于鄰域作用力的三維網(wǎng)格平滑方法[J].中國(guó)石油大學(xué)學(xué)報(bào)(自然科學(xué)版),2017,41(6):80-87.

WANG Haiqi, PENG Jiaqi, WANG Yi, et al. A 3D meshes smoothing method based on neighborhood forces[J].Journal of China University of Petroleum(Edition of Natural Science), 2017,41(6):80-87.

(編輯 修榮榮)