基于分層四叉樹的多分辨率數(shù)字巖心的表示與生成*

陳國軍 裴利強 李 勝 姜 朕

（中國石油大學(xué)（華東）計算機(jī)科學(xué)與技術(shù)學(xué)院 青島 266580）

1 引言

數(shù)字巖心技術(shù)是通過特定技術(shù)對巖心進(jìn)行模擬，建立巖石三維數(shù)字化圖像，根據(jù)數(shù)字巖心技術(shù)生成的數(shù)字巖心模型能夠反映真實巖心微觀孔隙結(jié)構(gòu)，模擬巖石物理屬性及實驗，可以精確地分析巖心各項特征和性質(zhì)。相比于傳統(tǒng)的巖石物理實驗，數(shù)字巖心技術(shù)快速簡便，省時省力且結(jié)果精細(xì)等優(yōu)點。

數(shù)字巖心的建模方法［1～5］根據(jù)建模數(shù)據(jù)的來源可以分為兩種：物理重建與數(shù)值重建，物理重建全部使用外部數(shù)據(jù)重建數(shù)字巖心，不通過計算機(jī)模擬巖心的結(jié)構(gòu)元素；數(shù)值重建則部分依靠外部數(shù)據(jù)，然后通過計算機(jī)模擬出數(shù)字巖心。物理重建分為X射線掃描成像法和切片組合法等［6～8］。

X射線CT技術(shù)出現(xiàn)和發(fā)展［9］，促進(jìn)了數(shù)字巖心技術(shù)的發(fā)展，通過CT設(shè)備掃描，能夠獲得分辨率極高的巖心二維圖像，使得建立精細(xì)、準(zhǔn)確的數(shù)字巖心孔隙骨架結(jié)構(gòu)更加簡單便易。

數(shù)字巖心的表現(xiàn)形式有多種，體素模型、面模型、孔隙網(wǎng)絡(luò)模型［10］等，其中體素模型是對真實數(shù)字巖心進(jìn)行數(shù)據(jù)建模的一種直觀表示方式，用一系列規(guī)則尺寸的體元（如立方體和正十二面體等）來模擬巖心空間。建立體素模型數(shù)字巖心的一般步驟：首先將二維巖心切片預(yù)處理，分割巖心骨架和孔隙像素，疊加到三維場景中，根據(jù)每一層切片像素信息生成對應(yīng)的體素。體素模型模擬數(shù)字巖心的特點是簡單、準(zhǔn)確，便于進(jìn)行巖心孔隙骨架空間分布的研究。但體素數(shù)字巖心容易帶來數(shù)據(jù)冗余、占用存儲空間大等問題，難以在普通計算機(jī)上進(jìn)行快速生成和渲染，因此，需要采用和合理的數(shù)據(jù)組織方式。通常，數(shù)字巖心部分結(jié)構(gòu)需要較高分辨率與精度，以獲得更多的細(xì)節(jié)信息，而其余部分則不需高分辨率，難以平衡分辨率問題。張晴等［11］分別對頁巖巖心采用高、低兩個分辨率獲取二維切片圖像，然后分別建立了二者的數(shù)字巖心，在低分辨率下建立的頁巖數(shù)字巖心，為了表征出內(nèi)部的納米級孔隙特征，選取特定的區(qū)域，在高分辨率下建立了相應(yīng)的數(shù)字巖心。崔利凱等［12］利用圖像配準(zhǔn)方法將不同分辨率下的巖心掃描圖像進(jìn)行空間配準(zhǔn)，按照分辨率從高到低的順序進(jìn)行孔隙及骨架分割，構(gòu)建多尺度、多組分的數(shù)字巖心模型。有些研究［13～18］采用八叉樹的方式組織數(shù)據(jù)，但八叉樹仍存在構(gòu)建較為耗時，遍歷復(fù)雜，維護(hù)管理不方便等問題，且在X、Y、Z三個方向上分辨率相同，在表示模型邊界時精確度不高，三個方向上尺寸相差較大的情況下造成了數(shù)據(jù)的冗余。劉亞靜［19］等對八叉樹數(shù)據(jù)結(jié)構(gòu)進(jìn)行擴(kuò)展，將形體的邊界信息單獨存儲，然后對三維礦體進(jìn)行重構(gòu)，提升了模型的邊界精度。

本文針對傳統(tǒng)的數(shù)字巖心建模方法造成的數(shù)據(jù)量過大，耗時過長，難以表示多分辨率等問題，提出了分層四叉樹模型，在此模型基礎(chǔ)上建立多分辨率數(shù)字巖心體素模型，同時結(jié)合MC算法［20］，生成數(shù)字巖心面模型。最后通過實驗給出了分層四叉樹與普通方式以及八叉樹在性能上的對比。

2 CT切片預(yù)處理

通過射線掃描獲得的CT切片為灰度圖，像素灰度值分布為0～255，需要設(shè)置閾值，將孔隙和骨架進(jìn)行分割，以分別獲得骨架和孔隙數(shù)據(jù)。閾值分割方法實際上通過對比閾值將閾值兩邊的像素點灰度值二值化，設(shè)分割閾值為T，像素點（i，j）的灰度值為f（i，j），若f（i，j）>T，則令其灰度值為255，否則為0，計算式如式（1）所示。

閾值分割算法的關(guān)鍵是確定閾值，如果能確定一個合適的閾值就可準(zhǔn)確地將圖像分割開來。閾值分割法主要分為全局和局部兩種，目前應(yīng)用的閉值分割方法都是在此基礎(chǔ)上發(fā)展起來的，比如最小誤差法、最大相關(guān)法、最大嫡法、矩量保持法、最大類間方差法等，而應(yīng)用最廣泛的是最大類間方差法。

本文采用最大類間方差法進(jìn)行計算閾值，再結(jié)合孔隙度，進(jìn)行閾值修正，最后根據(jù)閾值區(qū)分骨架和孔隙像素數(shù)據(jù)。

1）最大類間方差

最大類間方差法OTSU［21～23］是一種自適應(yīng)閾值確定的方法，基于全局的二值化算法，它是根據(jù)圖像的灰度特性，將圖像分為前景和背景兩個部分。當(dāng)取最佳閾值時，兩部分之間的差別應(yīng)該是最大的，在算法中所采用的衡量差別的標(biāo)準(zhǔn)就是較為常見的最大類間方差。

記T為前景與背景的分割閾值，前景點數(shù)占圖像比例為w0，平均灰度為u0；背景點數(shù)占圖像比例為w1，平均灰度為u1，圖像的總平均灰度為u，前景和背景圖像的方差，計算如式（2）所示。

當(dāng)方差g最大時，可以認(rèn)為此時前景和背景差異最大，此時的灰度T是最佳閾值。

2）閾值修正

采用最大類間方差法計算得到的閾值，只是從概率角度計算得到的，用這個閾值進(jìn)行分割像素，導(dǎo)致孔隙像素個數(shù)較多，得到的孔隙度過大，為此需要進(jìn)行修正。

根據(jù)巖心的孔隙度φ和最大類間方差法計算得到的閾值為T0及巖心總的像素個數(shù)C。分別計算T=T0±di時的孔隙像素個數(shù)Ci/C和Ci/C，最終閾值取|φ-Ci/C|最小時的T作為閾值，如式（3）和（4）所示。

3 分層四叉樹模型

四叉樹和八叉樹［24］是分別用來描述二維空間和三維空間的樹狀數(shù)據(jù)結(jié)構(gòu)，八叉樹每個節(jié)點下至多可以有八個節(jié)點，通常把三維空間劃分為八個區(qū)域，區(qū)域的信息存儲在每個節(jié)點中，若細(xì)分區(qū)域中沒有信息或具有一樣的性質(zhì)，則節(jié)點不再細(xì)分。四叉樹則用來表示二維空間。

本文提出分層四叉樹模型，根據(jù)本文實際情況將四叉樹擴(kuò)展到三維空間，采用多層四叉樹的方式組織CT切片數(shù)據(jù)，分層四叉樹結(jié)構(gòu)如圖1所示。

圖1 分層四叉樹示意圖

3.1 分層四叉樹模型定義

每層CT切片生成一棵四叉樹，由一個一維數(shù)組存儲每棵樹的根節(jié)點。四叉樹節(jié)點數(shù)據(jù)結(jié)構(gòu)包含節(jié)點類型、節(jié)點坐標(biāo)位置、節(jié)點的子節(jié)點指針等。

四叉樹任一節(jié)點的子節(jié)點為四個或零個，若節(jié)點范圍內(nèi)存在像素值不同的像素點，即不滿足式（15），則節(jié)點細(xì)分為四個子節(jié)點，子節(jié)點每個維度上 的 邊 長 分 別 為L1=（xmax-xmin）/2、L2=（ymaxymin）/2，若范圍內(nèi)像素值一致，即滿足式（5），則節(jié)點為葉節(jié)點，節(jié)點中存儲范圍坐標(biāo)、子節(jié)點指針及節(jié)點像素屬性等信息，對于四叉樹，所有葉節(jié)點滿足節(jié)點內(nèi)像素屬性一致。

單層CT切片生成的四叉樹如圖2所示。

圖2 單層CT切片四叉樹示意圖

3.2 分層四叉樹模型擴(kuò)展

分層四叉樹由于結(jié)構(gòu)特點，其管理較為快速簡便。遍歷：通過根節(jié)點數(shù)組直接獲得根節(jié)點指針，然后通過節(jié)點存儲的子節(jié)點指針進(jìn)行遍歷，直至滿足條件，對于樹T，其深度depth（T）=max｛log2rmax，log2cmax｝；拆分：將待拆分節(jié)點根據(jù)節(jié)點范圍坐標(biāo)一分為四，得到新的四個葉節(jié)點，對新的葉節(jié)點進(jìn)行賦值；合并：由于四叉樹節(jié)點子節(jié)點數(shù)要么為0要么為4，因此合并即為4個子節(jié)點合并，直接將4個子節(jié)點刪除，保留其父節(jié)點即可。

4 數(shù)字巖心的多分辨率模型生成

4.1 多分辨率擴(kuò)展

首先由低分辨率CT切片數(shù)據(jù)建立四叉樹模型。

通過高分辨率模板對低分辨率CT切片進(jìn)行修改時，低分辨率CT切片孔隙、骨架像素點可能會發(fā)生翻轉(zhuǎn)變化，因此需要對四叉樹進(jìn)行相應(yīng)的修改擴(kuò)展。設(shè)擴(kuò)展切片號為n，模板所在行列值范圍為rmin～rmax，cmin～cmax。

根據(jù)切片編號i從根節(jié)點數(shù)組中取出對應(yīng)四叉樹根節(jié)點RootNode=Slice（n）；從根節(jié)點通過子節(jié)點指針進(jìn)行遍歷，比較子節(jié)點Si坐標(biāo)范圍與模板所在行列范圍，若滿足rmin>xmin，rmaxymin，cmax

圖3 模板割分示意圖

4.2 體素數(shù)字巖心

依據(jù)根節(jié)點數(shù)組，遍歷每一層切片（即每一棵四叉樹），根據(jù)節(jié)點屬性由式（6）判斷節(jié)點類型，取出葉節(jié)點中存儲的坐標(biāo)信息，由xmin、xmax、ymin、ymax、z、z+1構(gòu)建立方體體元，同時根據(jù)節(jié)點屬性分別賦予骨架體元和孔隙體元不同渲染顏色。

等值面的梯度分量在沿面的切線方向為零，所以，該點的梯度矢量的方向就代表了該點的法向方向［15］。為了方便后續(xù)等值面法線的計算，體元各頂點法向量通過差分公式計算：

式中：Δx、Δy、Δz分別代表體元三個方向上的間距，此處Δx=Δy=Δz=1，gx、gy、gz分別代表坐標(biāo)點（i，j，k）三個方向上的梯度值，法向量為N。

4.3 面數(shù)字巖心

4.3.1 基于分層四叉樹的MC算法

Marching Cubes算法是一種經(jīng)典的等值面提取算法，它的原理是在由一系列二維圖像組成的三維圖像庫中選取上下連續(xù)兩層相鄰8個像素點作為定點組成的立方體體元中，根據(jù)8個頂點的像素值生成三角形截面，即等值面。對體元中8個角點的像素值進(jìn)行判斷，如果頂點值為0，則說明該頂點為孔隙點，將該點標(biāo)記為0，如果頂點值為255，則說明該點為骨架點，標(biāo)記為1，遍歷體元8個頂點，便得到孔隙、骨架點的分布情況，如果體元一條邊的兩個頂點屬性不一樣，那么截面一定經(jīng)過此邊，由此規(guī)律，可以找到體元中與截面三角形相交的邊，進(jìn)而分析出體元中截面的分布情況。體元頂點標(biāo)記只有0和1兩種情況，所以總共存在28=256種不同的體元。由旋轉(zhuǎn)和對稱變換可知256種體元可以總結(jié)為15種。

由分層四叉樹模型構(gòu)建體元：依據(jù)CT切片行列值坐標(biāo)i∈（rmin，rmax）、j∈（cmin，cmax），依次構(gòu)建體元的八個頂點，對于體元頂點p（i，j，k），對分層四叉樹k進(jìn)行遍歷，得到葉節(jié)點N，若滿足xmin

使用高分辨率模板對分層四叉樹模型進(jìn)行擴(kuò)展后，造成局部節(jié)點細(xì)分，分辨率變高，體元的構(gòu)建亦須相應(yīng)的改變。如圖4所示，節(jié)點N為擴(kuò)展節(jié)點，將其進(jìn)行標(biāo)記，視為高分辨率節(jié)點，構(gòu)建體元時須提高分辨率，然后遍歷其鄰接節(jié)點N1、N2、N3、N4、N5、N6，判斷節(jié)點N與鄰接節(jié)點交界處像素值是否一致，若與鄰接節(jié)點Ni像素不一致，則節(jié)點N與Ni間有等值面產(chǎn)生，構(gòu)建體元時須細(xì)分鄰接節(jié)點節(jié)點Ni，提高分辨率，將Ni進(jìn)行標(biāo)記，同時對Ni與其鄰接節(jié)點執(zhí)行上述操作；若與鄰接節(jié)點Ni像素一致，則停止標(biāo)記。至此，N節(jié)點及其周圍被標(biāo)記節(jié)點即為構(gòu)建體元時須細(xì)分節(jié)點。

圖4 局部擴(kuò)展示意圖

標(biāo)記完成后，進(jìn)行體元構(gòu)建時，若p（i，j，k）位于被標(biāo)記節(jié)點起始位置，即高分辨率節(jié)點，則此節(jié)點在構(gòu)建體元時須細(xì)分，此時以模板分辨率為體元邊長進(jìn)行體元頂點的搜索構(gòu)建體元，查找等值面，至節(jié)點結(jié)束位置。

4.3.2 截面頂點坐標(biāo)及法線計算

本文體元頂點只有孔隙和骨架兩種屬性，為了計算方便，截面頂點直接取所在邊中點，同時法線即所在邊頂點法線均值，如式（11）所示，其中N（i1，j1，k1），N（i2，j2，k2）分別為截面頂點所在邊兩端頂點法向量。

5 實驗結(jié)果分析

為了檢驗算法效果及時間、存儲上的性能，本文的實驗在Windows10平臺進(jìn)行，機(jī)器配置為E5-1620CPU，16G內(nèi)存，NVIDIA Quadro K2000圖形顯卡，開發(fā)工具有VS2015，OpenCV開發(fā)庫，OpenGL圖形庫，使用的數(shù)據(jù)為吉林某油藏儲層巖心CT掃描切片，圖像分辨率為995×968，數(shù)量為60張。采用本文方生成的數(shù)字巖心效果如圖5、圖6所示。

圖5 數(shù)字巖心模型

圖6 局部多分辨率擴(kuò)展前后對比圖

除了采用上述提出算法進(jìn)行實例建模，實驗還采用普通方法以及傳統(tǒng)的八叉樹對相同數(shù)據(jù)進(jìn)行建模實驗，與本文方法性能對比如表1、圖7、圖8所示。

圖7 不同模型內(nèi)存占用對比

圖8 體素數(shù)字巖心總生成時間對比

表1 不同方法性能對比

結(jié)果表明，與普通方式相比，采用八叉樹以及分層四叉樹，生成多尺度模型，能夠使體素量減少95%以上。與八叉樹相比，本文的分層四叉樹模型存儲空間占用減少30%左右，同時縮短了數(shù)字巖心模型重建耗時。

6 結(jié)語

實驗證明，本文采用的方法建立的多分辨率數(shù)字巖心模型能夠精確地模擬顯示真實巖心的空間特征，同時優(yōu)化了傳統(tǒng)建模方法存在的數(shù)據(jù)冗余，存儲空間消耗大，建模耗時等缺點。多尺度數(shù)字巖心模型通過分層四叉樹組織，利于后續(xù)模型的維護(hù)及擴(kuò)展等操作。