數字巖心多分辨率存儲及連通區域分析

彭成

(中國石油化工股份有限公司石油勘探開發研究院， 北京 100083)

隨著計算機技術的發展，運用圖像技術將電鏡或電子計算機斷層掃描(computed tomography，CT)掃描后的二維和三維巖心數據進行處理，進而形成數字巖心[1]。數字巖心可用多種計算機算法進行建模、孔隙度分析、孔喉結構分析等，是目前較為常用的巖心分析手段。王云龍等[2]基于綜合微計算機斷層掃描技術、數字巖心分析技術，實現低滲儲層數字巖心滲吸過程的模擬。宋文輝等[3]對數字巖心進行納米孔隙流動數學模型建立，研究了多尺度孔隙結構特征。

數字巖心的建模方法早期是通過二維圖片，通過一些典型算法建模出三維圖像，例如蒙特卡洛、模擬退火、過程法等。羅家榮[4]采用過程法模擬數字巖心，通過模擬巖石沉積過程建立數字巖心。

隨著計算機處理能力和掃描能力的進步，三維數字巖心圖像現在可以直接通過掃描獲取。在此基礎上，對數字巖心進行相應算法分析，構建孔隙模型、孔喉結構等。王志兵等[5]采用Avizo軟件中先進的數字圖像處理技術提取孔隙網絡模型。柳青兵等[6]采用傳統拼接方法，通過相同孔徑條件下流過孔隙結構的流體體積變化相等的拼接原理確定不同巖性樣品的拼接點位置，進而確定完整孔隙大小。

對于數字巖心存儲，目前多采用單一文件存儲，數據量大，查看局部數據時速度較慢，且沒有區分粗細粒度，應該進行分割存儲并存儲不同粗細級別的數據來方便查看[7-8]。陳國軍等[9]基于分層四叉樹對二維數字巖心進行了多分辨率數字巖心體素模型建立，不過沒有實現三維數字巖心相關的處理。

對于孔隙的連通性，各個連通區域之間的區分程度不夠高，應該將不同連通區域進行分色來幫助顯示[10-11]。對于每個連通區域的面積或體積，范圍，中心位置以及其他一些信息，應該建立連通區域和其屬性的對應關系。

針對上述問題，現通過四叉樹和八叉樹分塊方法對數字巖心進行分布式多級分辨率存取，根據瀏覽范圍選擇粒度層級；通過掃描線、三角剖分、四面體剖分和地圖染色方法識別和劃分連通區域；基于移動立方體算法求取巖心孔隙邊界面。通過上述方法，以期為數字巖心的高效存取、連通區域劃分和性質分析、孔喉結構的建模和形態分析提供支持。

1 數字巖心多分辨率分塊存儲

為提升數字巖心的存取效率，采用四叉樹和八叉樹分別對二維和三維數字巖心進行分塊切割存儲，使得在獲取局部范圍內的數字巖心數據時只需讀取區域對應的分塊即可。將分塊數據逐級向上合并形成多級分辨率，并在顯示較大范圍數字巖心時，根據屏幕像素選擇適合的分辨率級別讀取對應的分塊。

1.1 分塊參數配置

數字巖心分為二維三維兩種類型，格式主要包括RAW和BMP等。對于RAW格式，除巖心數據本身外不帶有其他信息，需要給出巖心的X、Y、Z(長、寬、高)像素數量。對于BMP格式，文件中帶有X、Y(長、寬)兩個方向的像素值，可以直接獲取，Z方向像素數量需要另外給出。接下來指定分塊大小，二維巖心設置分塊的X、Y像素數量，三維巖心設置X、Y、Z像素數量，后續會將數字巖心按照設定的分塊大小切割成塊。

分別選取二維和三維RAW格式文件，其中二維RAW文件的X、Y像素值分別為110 342和104 900，文件大小為10 GB，四叉樹子塊設為長寬1 024，即子塊大小為1 MB；三維RAW文件的X、Y、Z像素值分別為1 200、1 200、1 600，文件大小為2 GB，八叉樹子塊設為長寬高128，即子塊大小為2 MB。子塊的大小由用戶指定，為提升子塊分割效率，依照四叉樹和八叉樹的性質，子塊的長寬高大小選取的依據為：盡量保證分割后各個方向上切割出的數量相近，并接近且小于等于一個2的次冪的個數，這樣可以保證各個方向上切分后的層級相同，且空塊較少。

1.2 分塊切割

基于四叉樹和八叉樹對數字巖心進行分塊，首先計算樹的層級數，對于二維數字巖心數據，四叉樹的層數需同時滿足關系式如下。

2(n-1)Xblock≥Xfile

(1)

2(n-1)Yblock≥Yfile

(2)

對于三維數字巖心數據，八叉樹的層數還需滿足關系式為

2(n-1)Zblock≥Zfile

(3)

式中：n為層數；Xblock、Yblock、Zblock分別為子塊X、Y、Z像素值；Xfile、Yfile、Zfile分別為文件的X、Y、Z像素值。

滿足這些不等式的最小層數作為樹的層級數，依此計算得到本文選取的二維數字巖心數據四叉樹的層數為8，對應子塊的總個數為16 384個，三維數字巖心數據八叉樹的層數為5，按照八叉樹規則得到對應子塊的總個數為4 096個。基于四叉樹和八叉樹編碼規則對子塊賦予編號[12-13]，如圖1所示，為八叉樹空間結構劃分及線性莫頓編碼示意圖，按照圖中左側所示對三維數字巖心進行切割，每個子塊的編碼如圖1中右側所示。二維數字巖心按照類似方式基于四叉樹進行切割和編碼。

圖1 八叉樹空間結構劃分及線性莫頓編碼Fig.1 Octree spatial structure partition and linear Morton coding

接下來按子塊編號的順序依次讀取對應位置的數字巖心數據，對于靠近邊界區域的子塊，其長寬高可能達不到設定的子塊大小，截取到邊界即可，對于邊界外的子塊，生成空子塊。

數字巖心數據量較大，各分塊相對獨立且處理方式一致，適用于用多線程的方式進行讀取，按照四叉樹和八叉樹的性質，為二維數字巖心設置4的冪級數量的線程，為三維數字巖心設置8的冪級數量的線程，在例子中，設定的二維和三維線程數分別為16和8。為每個線程分配數量相同的子塊進行處理。

1.3 分塊合并

在取子塊數據的過程中，同時進行子塊合并工作。對于二維數字巖心數據，每獲取完4個子塊進行合并，如圖2(a)所示，即在X、Y方向間隔獲取像素進行粗化，紅色格子為所選像素，4個子塊間隔抽取合并出右邊粗化子塊，即在長寬方向上各進行50%的壓縮，合并完的子塊即上一粒度的粗子塊，實際數據效果如圖2(b)所示。合并后將這4個塊所占內存釋放。當上一級子塊數量也達到4個時，將它們進行合并存儲。依此直到最頂層。

圖2 二維數字巖心四叉樹分塊合并示意圖Fig.2 Two dimensional digital core quadtree block merging

對于三維數字巖心數據的合并，與二維類似，每獲取8個子塊進行合并成上一級粗粒度子塊，每級粒度的子塊滿8個時都進行合并，直到所有子塊合并完畢。按照邊讀取邊合并的方式可以大幅度減少內存，正常方式讀取所有子塊需要的內存為數字巖心文件大小，而邊讀取邊合并的方式占用的內存最多為

M=3Mblocklevel

(4)

式(4)中：level為樹層級數；Mblock為單個子塊占用內存大小。

對于多線程的模式，各個線程單獨完成上面的操作，每個線程做完所有的合并操作后，再將各個線程的塊按照同樣方法合并起來，直到最頂層。例子中二維線程數為16，即四叉樹的第二層，再將它們進行兩級合并即可。三維線程數為8，為八叉樹的第1層，再進行一級合并即可，最后將這些線程子塊合并到最上層即樹的第0層。

1.4 讀取粒度與區域內分塊選擇

在數字巖心文件導入并切割存儲完成后，可以進行查看，如圖3(a)所示，為本文例子中的二維數字巖心圖像，圖3(b)為三維數字巖心的3個方向切片切出來的立體圖像，可以指定任意方向的切片來查看三維數字巖心。在顯示方面，由于生成了不同粗細粒度的圖像，通過不同的四叉樹和八叉樹級別表示，越上層的圖像越粗。在選擇以何種粗細粒度顯示時，選擇的粒度如果過細，則受屏幕分辨率限制不能完全顯示所有像素，如果選擇的粒度過粗，則展現的像素量太少。選擇粒度的準則是在屏幕分辨率的限制下，盡可能選擇更細的粒度。

圖3 二維三維數字巖心顯示效果Fig.3 2D and 3D digital core display effect

對于數字巖心圖像粗細粒度的選擇，即四叉樹或八叉樹級別的選擇，二維三維數字巖心所取級別需分別滿足以下關系式。

Mapfile/(Pixblock4lev-1)

(5)

Mapfile/(Pixblock8lev-1)

(6)

式中：Mapfile為文件地理坐標范圍，式(5)中為面積范圍，式(6)中為體積范圍，在導入時設定，默認取值為其像素坐標范圍；Pixblock在式(5)中為子塊像素面積，在式(6)中為體積；lev為所取級別；Mapscr和Pixscr分別為屏幕地理坐標范圍和像素范圍。式(5)和式(6)分別對應二維和三維數字巖心。

對于局部區域巖心數據查看，根據屏幕所顯示的地理坐標范圍，以及數字巖心的地理坐標范圍的位置區域，得到對應要查看的巖心區域，讀取區域內子塊數據，子塊層級選擇按前面所述。

1.5 存取效率及粒度選擇實驗對比

實驗運行的環境為Windows 10，64位操作系統，4核8 GB內存，Intel(R) Core(TM) i7-10510U CPU @ 1.80 GHz 2.30 GHz。對1.1節中所述的兩個巖心數據進行分布式存儲。首先搭建分布式服務器，如圖4所示，建立GFS Master服務器負責子塊的切分和獲取；建立GFS Chunk存儲服務器來存放子塊數據，最后共建立了3個存儲服務器100個虛擬存儲節點來實現分布式存儲。

圖4 分布式服務器架構Fig.4 Distributed server architecture

按照1.1節的切塊配置進行切塊，分布式存儲到各個節點中，各數據的切塊結果如表1所示。

表1 實驗數據切塊結果Table 1 Experimental data segmentation results

與本地文件整體存儲相比，分布式分塊存儲一方面減少本地存儲空間開銷，另一方面在訪問局部區域數據時只需要讀取對應子塊，減少文件讀寫開銷。表2對比了整體與分塊存儲在訪問一個子塊區域對應的數據時，文件讀寫的跨度。

由表2可以看出，對于整體存儲模式數據連續存儲，在訪問子塊區域數據時，為讀取子塊區域對應的數據，那么會不斷跨整個寬度和深度數據去訪問下一個位置的數據，造成了文件訪問跨度的增加。而分塊存儲只需要訪問對應子塊文件，不需要額外的跨度。

天津市高等學校師資培訓中心（以下簡稱天津市中心）成立于1990年，受天津市教委和天津師范大學雙重領導。中心下設行政辦公室、崗前培訓辦公室、信息技術部、教師資格認定辦公室，當前業務以崗前培訓、教師資格認定、高校教師網絡培訓為主。

表2 訪問區域數據時文件讀寫跨度對比Table 2 Comparison of file read and write span when accessing area data

在訪問屏幕范圍內數據時，會選擇合適的粗細粒度子塊數據進行讀取，表3顯示了與選擇原始數據讀取相比，選擇適合粒度層級進行讀取在數據量上的對比，選擇的屏幕分辨率為1 024×1 024。

表3 選擇合適讀取粒度與原始數據對比Table 3 Select the appropriate reading granularity and compare it with the original data

可以看出，多級分辨率的巖性數據，根據瀏覽范圍選擇合適的粒度層級展現，依照1.4節中的方法即充分利用屏幕分辨率，又整體上減少了數據讀取量，屏幕顯示范圍越大，減少的數據量越多。

2 數字巖心連通區域識別

數字巖心的孔隙分析、孔喉結構模型計算、孔喉連通區域計算等是數字巖心分析的重要方面。孔隙的識別通過設定的像素閾值進行；對孔喉結構連通性的計算，首先基于掃描線方法，先計算連通線，再將相鄰連通線合并成連通區域，最后將各分塊連通區域合并成整體連通區域。

2.1 連通線掃描

首先進行掃描參數的設定，包括閾值、掃描線最小長度、連通區域面積或體積的下限等。依次讀取最細層級的每個子塊進行掃描，在掃描過程中，二維按照X、Y的順序，三維按照X、Y、Z的順序依次讀取子塊的每個點的像素值，對于滿足閾值的像素點，相鄰的連接起來形成線段，即掃描線，如果掃描線長度小于最小長度則剔除掉。

在例子中，二維數字巖心閾值設定為55，掃描線最小長度為3，連通區域面積下限為22，三維數字巖心閾值設定為55，掃描線最小長度為3，連通區域體積下限為50。即掃描二維數字巖心過程中，像素值小于55的點會記錄下來，連續3個或以上這種點可以組成一條掃描線，若干個相鄰的掃描線組成的連通區域總面積還要大于等于22；掃描三維數字巖心過程中，像素值小于55的點會記錄下來，連續3個或以上這種點可以組成一條掃描線，若干個相鄰的掃描線組成的連通區域總體積還要大于等于50。

2.2 分塊連通區域生成

生成完子塊所有的掃描線后，將掃描線相鄰合并形成連通區域，從而得到這個子塊的若干個連通區域。對于二維數字巖心，如果掃描線的Y方向距離小于等于1且X方向不相離，則相鄰；對于三維數字巖心，如果掃描線的Y、Z方向距離均小于等于1且X方向不相離，則相鄰。依照此判斷規則將所有掃描線組成一個個連通區域，生成一個連通區域的方法為：通過堆棧的方式，每次在未分組的掃描線中選取一條，判斷是否與當前堆棧頂部的掃描線相鄰，如果相鄰則加入連通區域和堆棧。所有掃描線選取完畢后將堆棧彈出一個對象添加到連通區域。接下來進行第二次循環，操作與第一次相同，以此類推，直到堆棧為空，便生成了一組掃描線表示的連通區域。如圖5所示，具有相鄰關系的一組連通線組成一個連通區域，共形成了3個連通區域。

圖5 掃描線合并連通區域示意Fig.5 Schematic diagram of scan line merging and connecting area

2.3 合并相鄰分塊連通區域

將相鄰分塊的連通區域進行合并，對于二維數字巖心，相鄰分塊為上下左右4個，對于三維數字巖心，相鄰分塊為上下左右前后6個。對于連通區域存在相鄰的兩個分塊，將它們的連通區域掃描線進行合并得到新的連通區域。

連通區域是否相鄰的判斷方法為：先判斷兩個連通區域的范圍是否有交集，如果沒有交集，則不相鄰，如果有交集，則進一步判斷：將它們所有的掃描線二維按照Y，三維先按照Z再按照Y進行排序，然后掃描線兩兩比較是否相鄰，如果Y或者Z距離大于1的直接跳過。一種快速的方法是，每次兩邊只取二維Y方向相鄰三行，三維Z方向相鄰三行的數據進行比較。對于三維，進一步對這些數據每次取Y方向相鄰三行進行比較。如果有相鄰的掃描線，則表示兩個連通區域相鄰。所有分塊的相鄰連通區域完成合并之后，最后將合并后的各個分塊的連通區域信息作為整體的連通區域信息。如圖6所示，相鄰兩個分塊的兩個連通區域在邊界處相鄰，子塊合并后這兩個連通區域合并成一個連通區域。

圖6 相鄰分塊連通區域合并示意Fig.6 Merging of adjacent block connected areas

2.4 合并整體連通區域

3 連通區域分色及邊界面計算

為方便查看各個孔喉連通區域形態，需要對其進行染色，基于三角剖分和四面體剖分對二維和三維連通區域分色，并結合地理信息技術將各連通區域作為圖元存儲形成圖層，圖元屬性包括顏色、面積或體積等，方便查詢和管理。基于移動立方體(marching cube，MC)算法[14]對三維數字巖心模型孔隙求取邊界面，以便在三維空間中展示孔喉結構。

3.1 三角和四面體剖分連通區域及分色

對各個連通區域的中點進行三角/四面體剖分連接起來。二維采用三角剖分算法，三維采用四面體剖分算法[15]。如圖7所示為例子中二維數字巖心連通區域進行三角剖分后的結果。

圖7 二維數字巖心三角剖分結果Fig.7 Two dimensional digital core triangulation results

利用三角/四面體剖分的結果，可以得到各個連通區域中心點之間是否有邊相連，以及每個點所連的邊的數量，如果兩個點有邊相連，則表示這兩個連通區域距離較近，需要染不同顏色。具體每個連通區域染哪種色的計算方法如圖8所示。

圖8 連通區域染色流程Fig.8 Dyeing process of connected areas

按照地圖染色算法，統計各個點連接的邊數量，從大到小排序，從邊數量最多的點開始染色，每染色一個點，將其四周有邊相連的點的邊計數減去1，然后將剩余點按照所連接的邊的數量從大到小排序，繼續從最多邊數量點染色，每次染色要與周圍有邊相連的點的顏色不同。然后依此類推循環下去，直到所有點都被染色。

接下來對于二維數字巖心，將染色結果用地理信息圖層表示。每個連通區域作為圖層中一個圖元對象，并賦予所染顏色。圖元對應屬性存放這個連通區域的中點坐標、面積、范圍、掃描線集合等。在例子中，二維數字巖心生成的地理信息圖層如圖9所示。

圖9 二維數字巖心連通區域識別Fig.9 Recognition of connected area of two-dimensional digital core

對于三維數字巖心，構建三維模型并將各個連通區域作為三維圖元對象，存放的屬性信息除掃描線集合、中點、范圍外，還存儲了連通區域的像素體積和地理上的體積。在例子中，三維數字巖心的連通區域染色效果如圖10所示，展示了3個方向上的切片。

圖10 三維數字巖心連通區域識別Fig.10 Recognition of connected area of 3D digital core

3.2 連通區域邊界面計算

三維數字巖心的連通區域表示孔洞，將孔洞剔除掉，剩余的部分代表喉道，構建喉道三維模型的流程如圖11所示。

圖11 構建喉道三維模型流程Fig.11 Process of constructing laryngeal three-dimensional model

對數字巖心的三維空間范圍中每個像素點，如果其屬于某個連通區域內，則賦值為0，如果不屬于連通區域內，則賦值為1，然后用MC算法，對每個立方體進行處理，得到各個立方體的01分界面，然后合并成整體的01分界面即孔洞分界面[11]。具體的方法為，按Z方向每次取相鄰的兩層切片，切片中夾了若干個立方體，數量為

n=(Xpix-1)(Ypix-1)

(7)

式(7)中：Xpix為長像素值；Ypix為寬像素值。

對于立方體的8個頂點，01分界面共有256種情況，可以分為15類，如圖12所示，列出了不同情況的等值面構建方法。將每個立方體形成的等值面連接起來，即是孔喉的分界面。連接時，與連通區域分組類似，將具有相連關系的分界面劃分為一組，屬于同一個連通孔喉，互不相連的分界面屬于不同孔喉。

圖12 移動立方體等勢面構建情況組合Fig.12 Combination of construction of equipotential surface of moving cube

例子中將三維數字巖心的等值面生成后的孔喉結構效果如圖13所示，可以用兩種顯示方式，圖13(a)～圖13(d)是顯示等值面，圖13(e)～圖13(h)是只顯示等值面的邊即鏤空模式，從左到右分別是從小到大設定不同閾值得到的結果，閾值分別為50、70、90、110。閾值設定較小，得到的連通區域小，剔除后剩余的部分較多；閾值設定大，得到的連通區域大，剔除后剩余的部分較少。

圖13 三維數字巖心連通區域邊界面Fig.13 Edge interface of 3D digital core connected area

4 結論

(1)分別以四叉樹和八叉樹模式分塊存儲二維和三維數字巖心數據，生成多級分辨率的巖性數據，根據瀏覽范圍選擇合適的粒度層級展現，整體上減少了數據讀取量。

(2)通過掃描線合并構建連通區域，基于移動立方體方法獲取孔隙分界面，得到了孔喉結構的形態及連通性質。通過三角剖分和四面體剖分對連通區域分色，方便區分顯示各連通區域。