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

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

1 引言

數(shù)字巖心技術(shù)是通過特定技術(shù)對巖心進(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é)果精細(xì)等優(yōu)點(diǎn)。

數(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í)數(shù)字巖心進(jìn)行數(shù)據(jù)建模的一種直觀表示方式，用一系列規(guī)則尺寸的體元（如立方體和正十二面體等）來模擬巖心空間。建立體素模型數(shù)字巖心的一般步驟：首先將二維巖心切片預(yù)處理，分割巖心骨架和孔隙像素，疊加到三維場景中，根據(jù)每一層切片像素信息生成對應(yīng)的體素。體素模型模擬數(shù)字巖心的特點(diǎn)是簡單、準(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ò)展，將形體的邊界信息單獨(dú)存儲，然后對三維礦體進(jìn)行重構(gòu)，提升了模型的邊界精度。

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

2 CT切片預(yù)處理

通過射線掃描獲得的CT切片為灰度圖，像素灰度值分布為0～255，需要設(shè)置閾值，將孔隙和骨架進(jìn)行分割，以分別獲得骨架和孔隙數(shù)據(jù)。閾值分割方法實(shí)際上通過對比閾值將閾值兩邊的像素點(diǎn)灰度值二值化，設(shè)分割閾值為T，像素點(diǎn)（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為前景與背景的分割閾值，前景點(diǎn)數(shù)占圖像比例為w0，平均灰度為u0；背景點(diǎn)數(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é)點(diǎn)下至多可以有八個節(jié)點(diǎn)，通常把三維空間劃分為八個區(qū)域，區(qū)域的信息存儲在每個節(jié)點(diǎn)中，若細(xì)分區(qū)域中沒有信息或具有一樣的性質(zhì)，則節(jié)點(diǎn)不再細(xì)分。四叉樹則用來表示二維空間。

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

圖1 分層四叉樹示意圖

3.1 分層四叉樹模型定義

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

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

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

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

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

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

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

4.1 多分辨率擴(kuò)展

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

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

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

圖3 模板割分示意圖

4.2 體素數(shù)字巖心

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

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

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

4.3 面數(shù)字巖心

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

Marching Cubes算法是一種經(jīng)典的等值面提取算法，它的原理是在由一系列二維圖像組成的三維圖像庫中選取上下連續(xù)兩層相鄰8個像素點(diǎn)作為定點(diǎn)組成的立方體體元中，根據(jù)8個頂點(diǎn)的像素值生成三角形截面，即等值面。對體元中8個角點(diǎn)的像素值進(jìn)行判斷，如果頂點(diǎn)值為0，則說明該頂點(diǎn)為孔隙點(diǎn)，將該點(diǎn)標(biāo)記為0，如果頂點(diǎn)值為255，則說明該點(diǎn)為骨架點(diǎn)，標(biāo)記為1，遍歷體元8個頂點(diǎn)，便得到孔隙、骨架點(diǎn)的分布情況，如果體元一條邊的兩個頂點(diǎn)屬性不一樣，那么截面一定經(jīng)過此邊，由此規(guī)律，可以找到體元中與截面三角形相交的邊，進(jìn)而分析出體元中截面的分布情況。體元頂點(diǎ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)建體元的八個頂點(diǎn)，對于體元頂點(diǎn)p（i，j，k），對分層四叉樹k進(jìn)行遍歷，得到葉節(jié)點(diǎn)N，若滿足xmin

使用高分辨率模板對分層四叉樹模型進(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)建體元時須提高分辨率，然后遍歷其鄰接節(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)建體元時須細(xì)分鄰接節(jié)點(diǎn)節(jié)點(diǎn)Ni，提高分辨率，將Ni進(jìn)行標(biāo)記，同時對Ni與其鄰接節(jié)點(diǎn)執(zhí)行上述操作；若與鄰接節(jié)點(diǎn)Ni像素一致，則停止標(biāo)記。至此，N節(jié)點(diǎn)及其周圍被標(biāo)記節(jié)點(diǎn)即為構(gòu)建體元時須細(xì)分節(jié)點(diǎn)。

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

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

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

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

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

為了檢驗(yàn)算法效果及時間、存儲上的性能，本文的實(shí)驗(yàn)在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)行實(shí)例建模，實(shí)驗(yàn)還采用普通方法以及傳統(tǒng)的八叉樹對相同數(shù)據(jù)進(jìn)行建模實(shí)驗(yàn)，與本文方法性能對比如表1、圖7、圖8所示。

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

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

表1 不同方法性能對比

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

6 結(jié)語

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