劉駿標(biāo),邢會林,姜素華,閆偉超

(1.深海圈層與地球系統(tǒng)教育部前沿科學(xué)中心,海底科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,中國海洋大學(xué) 海洋地球科學(xué)學(xué)院,山東 青島 266100;2.嶗山實(shí)驗(yàn)室,山東 青島 266237;3.中國海洋大學(xué) 海底科學(xué)與工程計(jì)算國際中心,山東 青島 266100;4.嶗山實(shí)驗(yàn)室 海洋礦產(chǎn)資源評價(jià)與探測技術(shù)功能實(shí)驗(yàn)室,山東 青島 266237)

近年來,數(shù)字巖石物理(Digital Rock Physics,DRP)技術(shù)被廣泛應(yīng)用于二維或三維巖石圖像處理中[1-3],以取代復(fù)雜的巖石物理實(shí)驗(yàn)(如滲透率、彈性模量及電導(dǎo)率等巖石物理參數(shù))。得益于CT等成像技術(shù)的發(fā)展,可以得到巖石的不同分辨率的二維或三維圖像(從納米到厘米尺度)[4-5]。裂隙作為巖石非均質(zhì)性的重要組成部分,可在高分辨率CT圖像中被識別出來[6]。若按照目前常用的方法,在進(jìn)行圖像可視化及后續(xù)數(shù)值模擬時(shí),為了精確表示出細(xì)小的裂隙,需要將所有體素點(diǎn)全部轉(zhuǎn)為規(guī)則正方體網(wǎng)格。然而,這種方法不僅將產(chǎn)生一個巨大的網(wǎng)格數(shù)據(jù)集,還將額外增加模擬計(jì)算量。因此,需要提出一種在保留裂隙主要幾何特征的同時(shí)實(shí)現(xiàn)數(shù)據(jù)簡化及生成計(jì)算分析所必需的網(wǎng)格的方法。

基于有限元或有限體積法的數(shù)值模擬分析依賴于可用的合適網(wǎng)格[7]。K.Terada等[8]將所有體素點(diǎn)直接轉(zhuǎn)換為正方體網(wǎng)格,此方法簡單直接,但將產(chǎn)生更大的數(shù)據(jù)集和鋸齒邊界,且鋸齒邊界處的模擬精度難以得到保證并最終影響模擬結(jié)果[9];W.E.Lorensen等[10]、Wu Z.J.等[11]、Zhang Y.J.等[12-13]、D.Boltcheva等[14]分別就移動立方體法及其拓展、基于八叉樹的四面體或六面體網(wǎng)格生成方法、德勞內(nèi)細(xì)化算法進(jìn)行了深入研究,并應(yīng)用于網(wǎng)格生成領(lǐng)域,此類方法通過識別出外輪廓或材料之間的邊界曲面,生成相應(yīng)的體積模型,但是在生成裂隙與圍巖之間的邊界曲面時(shí),由于裂隙厚度較小,需要生成大量小尺寸的網(wǎng)格才能保證材料邊界曲面的質(zhì)量,且其生成的三維網(wǎng)格模型并不適用于仿真分析;A.Paluszny等[15]和Liu Y.等[16-17]通過一系列三維線或面表示裂隙,此類方法能用有限數(shù)量的高質(zhì)量三角形網(wǎng)格來描述裂隙的主要幾何特征,顯著減少了網(wǎng)格數(shù)目,并能以三維裂隙面為約束生成非結(jié)構(gòu)化四面體網(wǎng)格,因此本文在后續(xù)將使用該方法。

本文提出了一種有效的特征提取、簡化及網(wǎng)格生成方法,主要工作為:①利用中心面點(diǎn)提取算法對裂隙的特征進(jìn)行提取,將裂隙數(shù)據(jù)體簡化為一系列面;②利用簡化質(zhì)心Voronoi法來對一系列面生成簡化質(zhì)心Voronoi圖,并基于質(zhì)心點(diǎn)進(jìn)行裂隙面網(wǎng)格生成、面片交叉標(biāo)識、網(wǎng)格修復(fù)及優(yōu)化;③采用對生成的裂隙面網(wǎng)格進(jìn)行圖像數(shù)據(jù)體還原的方式,對所生成的網(wǎng)格的代表性及其精度進(jìn)行體素相似度和形狀相似度評估;④以所生成的裂隙面網(wǎng)格和相應(yīng)的巖石邊界為約束,生成四面體網(wǎng)格模型,并用最小二面角、半徑比、邊長比及體積比等參數(shù)對四面體網(wǎng)格質(zhì)量進(jìn)行評價(jià)。

1 裂隙特征提取簡化網(wǎng)格生成及評價(jià)方法

裂隙特征提取是指從巖石數(shù)字圖像中提取裂隙的形態(tài)、大小、分布等幾何特征的過程,旨在對巖石的微觀結(jié)構(gòu)進(jìn)行精確分析和建模。本文首先會對原始數(shù)字巖石圖像數(shù)據(jù)進(jìn)行二值化、拓展[17]、膨脹與腐蝕等預(yù)處理操作。其中,二值化操作后數(shù)據(jù)將被分為基質(zhì)和裂隙,拓展操作則能將裂隙面上的小孔洞補(bǔ)全,而先腐蝕后膨脹操作可以消除3D圖像上的細(xì)小噪聲并平滑裂隙邊界。在預(yù)處理之后,本文將使用以下方法對三維巖石圖像裂隙進(jìn)行特征提取、簡化、網(wǎng)格生成、還原及評估,然后以生成的面為約束生成四面體網(wǎng)格。

1.1 裂隙特征提取、簡化及其網(wǎng)格生成

為將裂隙提取為一系列三維面從而得到裂隙的幾何特征,本文采用中心面點(diǎn)提取算法提取出裂隙的中心面點(diǎn)[16]。某個裂隙體素點(diǎn)v的26鄰域如圖1所示,分為8個區(qū)域(每個小立方體為一獨(dú)立區(qū)域),如果這8個區(qū)域都符合公式(1),則認(rèn)為這個裂隙體素點(diǎn)v為中心面點(diǎn)。

(1)

式中:N(v)為裂隙體素點(diǎn)v的26鄰域,Ni(v)為N(v)第i個分區(qū),|·|代表每個分區(qū)中裂隙體素點(diǎn)的數(shù)目。例如,對一個81×81×3的裂隙數(shù)據(jù)進(jìn)行中心面點(diǎn)提取算法后,得到79×79×1的裂隙中心面點(diǎn)數(shù)據(jù)。

Vi={x∈Ω|dis(x,zi)

(2)

式中:dis(·)為RN中的歐幾里得范數(shù),生成器zi也可被稱為Voronoi單元的質(zhì)心:

(3)

式中:ρ(x)是Vi的密度函數(shù)。

簡化質(zhì)心Voronoi圖法對Voronoi圖法做了2點(diǎn)限制:①假設(shè)各組裂隙被提取為一系列的平面,厚度變?yōu)?,裂隙由一組具有相同體積和密度的體素表示,因此公式(3)中的ρ(x)=1;②由于裂隙被提取為一系列的平面,在構(gòu)建Voronoi圖時(shí),采用廣度優(yōu)先搜索方法(Breadth First Search Method)傳播Voronoi單元[19]。

(4)

式中:Cvol為裂隙的總體素點(diǎn)數(shù)目,t為裂隙的平均厚度。r為自定義的Voronoi半徑,決定著Voronoi單元的傳播距離。因?yàn)镃vol和t是已知的,所以Voronoi半徑r為簡化質(zhì)心Voronoi圖法的唯一變量。

利用簡化質(zhì)心Voronoi圖法生成n個質(zhì)心點(diǎn)后,本文將基于這些質(zhì)心點(diǎn)生成三角形網(wǎng)格,并進(jìn)行面片交叉標(biāo)識、網(wǎng)格修復(fù)及優(yōu)化[16],用上述方法對79×79×1的裂隙中心面點(diǎn)數(shù)據(jù)處理后結(jié)果如圖2所示。

圖2 三角形網(wǎng)格生成流程Fig.2 Generation process of triangular mesh

1.2 相似度計(jì)算

本文采用對生成的網(wǎng)格進(jìn)行圖像數(shù)據(jù)體還原的方式,對所生成的網(wǎng)格的代表性及其精度進(jìn)行相似度評估。在進(jìn)行圖像數(shù)據(jù)體還原前,需要在生成的裂隙網(wǎng)格上添加厚度屬性,厚度的計(jì)算方法如下:

a.對所有裂隙點(diǎn)進(jìn)行判斷:當(dāng)一個點(diǎn)的上下左右前后有一個體素點(diǎn)為基質(zhì)時(shí),這個點(diǎn)為外邊界點(diǎn);

b.算出每個三角形網(wǎng)格的重心與所有外圍點(diǎn)的最小距離r;

c.將最小距離r乘以2作為厚度。

在賦予三角形網(wǎng)格厚度屬性后,本文會把帶厚度屬性的三角形網(wǎng)格還原回三維圖像數(shù)據(jù),數(shù)據(jù)還原方法如下:

a.當(dāng)三角形網(wǎng)格不平行于XY、XZ、YZ平面時(shí),利用三角形網(wǎng)格的重心與所有外圍點(diǎn)的最小距離求厚度值方法一般會導(dǎo)致厚度值偏小,所以本文對所有三角形網(wǎng)格都加上一個額外厚度Δd。在本文中,額外厚度Δd設(shè)置為1;

b.根據(jù)每個三角形網(wǎng)格的厚度將每個三角形推成三棱柱;

c.標(biāo)記此三棱柱內(nèi)部的所有體素點(diǎn),這些被標(biāo)記的體素點(diǎn)就是簡化后網(wǎng)格所能還原出來的體素點(diǎn)。

數(shù)據(jù)還原過程如圖3所示。

圖3 數(shù)據(jù)還原過程Fig.3 Data restoration process

數(shù)據(jù)還原后,將進(jìn)行以下2種相似度判斷[17]:

a.體素相似度:

(5)

b.形狀相似度:

(6)

式中:Cmesh代表的是還原后裂隙點(diǎn)的數(shù)目,Cimage代表的是原始數(shù)據(jù)裂隙點(diǎn)的數(shù)目。Cmesh,in代表的是還原后的裂隙與輸入圖像裂隙相同點(diǎn)的數(shù)目,Cmesh,out表示還原后裂隙與輸入圖像裂隙中不同點(diǎn)的數(shù)目,即Cmesh,out=Cmesh-Cmesh,in,Svoxel的取值范圍為[0,1],Sshape的取值范圍為[-∞,1]。

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

為證明本文方法能有效地應(yīng)用于具有復(fù)雜裂隙的三維巖石圖像,本文選取了一個128×128×60的二值化巖石裂隙圖像數(shù)據(jù),圖像數(shù)據(jù)如圖4所示。

圖4 數(shù)字巖石圖像與裂隙Fig.4 Digital rock image and fractures

從圖4-A中可看出此圖像數(shù)據(jù)已被分為基質(zhì)(藍(lán)色)和裂隙部分(紅色),圖4-B為該數(shù)據(jù)的裂隙部分。用本文方法對數(shù)據(jù)的裂隙部分進(jìn)行特征提取、簡化及網(wǎng)格生成,以不同Voronoi半徑為參數(shù)生成的網(wǎng)格如圖5所示。

從圖5可看出,隨著Voronoi半徑的增大,裂隙面網(wǎng)格越來越大,圖像也越加簡化并丟失更多信息,尤其是寬度小或長度小的裂隙更容易被簡化掉。Voronoi半徑與相似度的關(guān)系如圖6所示,網(wǎng)格數(shù)據(jù)如表1所示。

從圖6可以看出,隨著Voronoi半徑的增大,體素相似度與形狀相似度都降低,且體素相似度與形狀相似度的差基本不變。原始數(shù)據(jù)中裂隙點(diǎn)數(shù)目為110 467,占總體素點(diǎn)數(shù)目11.24%,表1中的簡化后網(wǎng)格數(shù)目表明本文方法能用少量的網(wǎng)格代表大量的裂隙點(diǎn)。此外,當(dāng)形狀相似度高于60%時(shí),本文方法生成的網(wǎng)格就能較好地代表原始裂隙。例如,當(dāng)Voronoi半徑選為3時(shí),形狀相似度為64.61%。從圖4-B與圖5-A中可以看出,生成的網(wǎng)格基本保留了裂隙的主要幾何特征。

圖5 基于不同Voronoi半徑生成的裂隙面網(wǎng)格Fig.5 Generated fracture surface meshes based on different Voronoi radius

表1 網(wǎng)格數(shù)據(jù)Table 1 Mesh data

圖6 Voronoi半徑與相似度的關(guān)系Fig.6 The relationship between Voronoi radius and similarity

考慮到原始數(shù)據(jù)是三維圖像數(shù)據(jù),對體素點(diǎn)進(jìn)行均勻采樣可以降低圖像分辨率、簡化數(shù)據(jù)和減少數(shù)據(jù)集的規(guī)模,但這種方法存在相似度急劇下降和體素相似度與形狀相似度之差增大的問題,如表2所示。表1中顯示,采用本文方法計(jì)算的體素相似度與形狀相似度之差最大僅為2.05%。與之相比,表2中均勻采樣法的體素相似度與形狀相似度之差最小就高達(dá)9.8%。因此,相比于均勻采樣,本文方法能保持體素相似度與形狀相似度的差基本不變,說明本文所提簡化方案可以較好地代表原始裂隙。

本文以裂隙面網(wǎng)格(圖5-A)和相應(yīng)的巖石邊界為約束,生成四面體網(wǎng)格模型(圖7)。從圖7中黑框區(qū)域可看出,四面體網(wǎng)格模型是基于裂隙面網(wǎng)格的節(jié)點(diǎn)約束而生成的。

通過采用最小二面角、半徑比、邊長比和體積比等4種網(wǎng)格質(zhì)量度量參數(shù)[20],可以評價(jià)圖7中四面體網(wǎng)格模型的網(wǎng)格質(zhì)量,計(jì)算結(jié)果如表3所示。從表3中的質(zhì)量數(shù)據(jù)范圍及其平均數(shù)據(jù)可以看出,生成的四面體網(wǎng)格平均質(zhì)量接近于理想的正四面體網(wǎng)格,網(wǎng)格質(zhì)量較好。

表2 不同采樣率時(shí)的圖像數(shù)據(jù)Table 2 Image data at different sampling rates

圖7 共有9 137節(jié)點(diǎn)和46 486個四面體網(wǎng)格的四面體網(wǎng)格模型Fig.7 A tetrahedral mesh model with a total of 9 137 nodes and 46 486 tetrahedral meshes

表3 網(wǎng)格質(zhì)量Table 3 Mesh quality

3 結(jié) 論

本文以一個128×128×60的二值化巖石裂隙圖像數(shù)據(jù)為例子,將裂隙提取為一系列平面并進(jìn)行網(wǎng)格生成與相似度評估,然后以所生成的裂隙面網(wǎng)格與相應(yīng)的巖石邊界為約束,利用自行研制的四面體網(wǎng)格生成程序生成了質(zhì)量較好的四面體網(wǎng)格模型。得到如下主要結(jié)論:

a.提出了一種既能保留裂隙主要特征又能顯著降低數(shù)據(jù)量的網(wǎng)格生成方法,通過中心面點(diǎn)提取算法、網(wǎng)格生成和厚度計(jì)算,達(dá)到用有限數(shù)量的高質(zhì)量三角形網(wǎng)格來描述裂隙的主要幾何特征的目的,并探討了不同Voronoi半徑參數(shù)與方法效果的關(guān)系。

b.在示例中,形狀相似度高于60%的裂隙面網(wǎng)格就能較好地代表原始裂隙。當(dāng)Voronoi半徑為3時(shí),基于本文方法生成的裂隙面網(wǎng)格能較好地代表原始裂隙,示例的體素相似度與形狀相似度為66.14%和64.61%。相比于原始數(shù)據(jù),簡化后的裂隙面網(wǎng)格數(shù)與初始裂隙體素點(diǎn)數(shù)之比為約為1：39,生成的四面體網(wǎng)格數(shù)與圖像總體素點(diǎn)數(shù)之比約為1：21;相比于所有體素點(diǎn)轉(zhuǎn)換為正方體網(wǎng)格,生成的四面體網(wǎng)格與正方體網(wǎng)格占用存儲空間之比約為1：59。

c.隨著Voronoi半徑的增大,雖然數(shù)據(jù)量能被進(jìn)一步降低,但一些寬度或長度較短的裂隙將更容易被簡化掉,從而造成形狀相似度進(jìn)一步下降。當(dāng)Voronoi半徑大于等于5后,示例的形狀相似度已經(jīng)小于60%,基于本文方法生成的裂隙面網(wǎng)格僅能代表原始裂隙中較為明顯的裂隙。