陳國(guó)軍 裴利強(qiáng) 李 勝 姜 朕

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

1 引言

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

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

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

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

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

2 CT切片預(yù)處理

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

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

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

1）最大類(lèi)間方差

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

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

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

2）閾值修正

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

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

3 分層四叉樹(shù)模型

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

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

圖1 分層四叉樹(shù)示意圖

3.1 分層四叉樹(shù)模型定義

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

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

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

圖2 單層CT切片四叉樹(shù)示意圖

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

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

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

4.1 多分辨率擴(kuò)展

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

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

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

圖3 模板割分示意圖

4.2 體素?cái)?shù)字巖心

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

等值面的梯度分量在沿面的切線(xiàn)方向?yàn)榱?，所以，該點(diǎn)的梯度矢量的方向就代表了該點(diǎn)的法向方向［15］。為了方便后續(xù)等值面法線(xiàn)的計(jì)算，體元各頂點(diǎn)法向量通過(guò)差分公式計(jì)算：

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

4.3 面數(shù)字巖心

4.3.1 基于分層四叉樹(shù)的MC算法

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

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

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

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

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

4.3.2 截面頂點(diǎn)坐標(biāo)及法線(xiàn)計(jì)算

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

5 實(shí)驗(yàn)結(jié)果分析

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

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

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

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

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

圖8 體素?cái)?shù)字巖心總生成時(shí)間對(duì)比

表1 不同方法性能對(duì)比

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

6 結(jié)語(yǔ)

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