基于圖像的三維裂隙巖體特征提取、簡化及網格生成方法

劉駿標,邢會林,姜素華,閆偉超

(1.深海圈層與地球系統教育部前沿科學中心,海底科學與探測技術教育部重點實驗室,中國海洋大學 海洋地球科學學院,山東 青島 266100;2.嶗山實驗室,山東 青島 266237;3.中國海洋大學 海底科學與工程計算國際中心,山東 青島 266100;4.嶗山實驗室 海洋礦產資源評價與探測技術功能實驗室,山東 青島 266237)

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

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

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

1 裂隙特征提取簡化網格生成及評價方法

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

1.1 裂隙特征提取、簡化及其網格生成

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

(1)

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

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

(2)

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

(3)

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

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

(4)

式中:Cvol為裂隙的總體素點數目,t為裂隙的平均厚度。r為自定義的Voronoi半徑,決定著Voronoi單元的傳播距離。因為Cvol和t是已知的,所以Voronoi半徑r為簡化質心Voronoi圖法的唯一變量。

利用簡化質心Voronoi圖法生成n個質心點后,本文將基于這些質心點生成三角形網格,并進行面片交叉標識、網格修復及優化[16],用上述方法對79×79×1的裂隙中心面點數據處理后結果如圖2所示。

圖2 三角形網格生成流程Fig.2 Generation process of triangular mesh

1.2 相似度計算

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

a.對所有裂隙點進行判斷:當一個點的上下左右前后有一個體素點為基質時,這個點為外邊界點;

b.算出每個三角形網格的重心與所有外圍點的最小距離r;

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

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

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

b.根據每個三角形網格的厚度將每個三角形推成三棱柱;

c.標記此三棱柱內部的所有體素點,這些被標記的體素點就是簡化后網格所能還原出來的體素點。

數據還原過程如圖3所示。

圖3 數據還原過程Fig.3 Data restoration process

數據還原后,將進行以下2種相似度判斷[17]:

a.體素相似度:

(5)

b.形狀相似度:

(6)

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

2 實際應用

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

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

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

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

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

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

表1 網格數據Table 1 Mesh data

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

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

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

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

表2 不同采樣率時的圖像數據Table 2 Image data at different sampling rates

圖7 共有9 137節點和46 486個四面體網格的四面體網格模型Fig.7 A tetrahedral mesh model with a total of 9 137 nodes and 46 486 tetrahedral meshes

表3 網格質量Table 3 Mesh quality

3 結 論

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

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

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

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