彭成

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

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

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

隨著計(jì)算機(jī)處理能力和掃描能力的進(jìn)步，三維數(shù)字巖心圖像現(xiàn)在可以直接通過(guò)掃描獲取。在此基礎(chǔ)上，對(duì)數(shù)字巖心進(jìn)行相應(yīng)算法分析，構(gòu)建孔隙模型、孔喉結(jié)構(gòu)等。王志兵等[5]采用Avizo軟件中先進(jìn)的數(shù)字圖像處理技術(shù)提取孔隙網(wǎng)絡(luò)模型。柳青兵等[6]采用傳統(tǒng)拼接方法，通過(guò)相同孔徑條件下流過(guò)孔隙結(jié)構(gòu)的流體體積變化相等的拼接原理確定不同巖性樣品的拼接點(diǎn)位置，進(jìn)而確定完整孔隙大小。

對(duì)于數(shù)字巖心存儲(chǔ)，目前多采用單一文件存儲(chǔ)，數(shù)據(jù)量大，查看局部數(shù)據(jù)時(shí)速度較慢，且沒(méi)有區(qū)分粗細(xì)粒度，應(yīng)該進(jìn)行分割存儲(chǔ)并存儲(chǔ)不同粗細(xì)級(jí)別的數(shù)據(jù)來(lái)方便查看[7-8]。陳國(guó)軍等[9]基于分層四叉樹(shù)對(duì)二維數(shù)字巖心進(jìn)行了多分辨率數(shù)字巖心體素模型建立，不過(guò)沒(méi)有實(shí)現(xiàn)三維數(shù)字巖心相關(guān)的處理。

對(duì)于孔隙的連通性，各個(gè)連通區(qū)域之間的區(qū)分程度不夠高，應(yīng)該將不同連通區(qū)域進(jìn)行分色來(lái)幫助顯示[10-11]。對(duì)于每個(gè)連通區(qū)域的面積或體積，范圍，中心位置以及其他一些信息，應(yīng)該建立連通區(qū)域和其屬性的對(duì)應(yīng)關(guān)系。

針對(duì)上述問(wèn)題，現(xiàn)通過(guò)四叉樹(shù)和八叉樹(shù)分塊方法對(duì)數(shù)字巖心進(jìn)行分布式多級(jí)分辨率存取，根據(jù)瀏覽范圍選擇粒度層級(jí)；通過(guò)掃描線、三角剖分、四面體剖分和地圖染色方法識(shí)別和劃分連通區(qū)域；基于移動(dòng)立方體算法求取巖心孔隙邊界面。通過(guò)上述方法，以期為數(shù)字巖心的高效存取、連通區(qū)域劃分和性質(zhì)分析、孔喉結(jié)構(gòu)的建模和形態(tài)分析提供支持。

1 數(shù)字巖心多分辨率分塊存儲(chǔ)

為提升數(shù)字巖心的存取效率，采用四叉樹(shù)和八叉樹(shù)分別對(duì)二維和三維數(shù)字巖心進(jìn)行分塊切割存儲(chǔ)，使得在獲取局部范圍內(nèi)的數(shù)字巖心數(shù)據(jù)時(shí)只需讀取區(qū)域?qū)?yīng)的分塊即可。將分塊數(shù)據(jù)逐級(jí)向上合并形成多級(jí)分辨率，并在顯示較大范圍數(shù)字巖心時(shí)，根據(jù)屏幕像素選擇適合的分辨率級(jí)別讀取對(duì)應(yīng)的分塊。

1.1 分塊參數(shù)配置

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

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

1.2 分塊切割

基于四叉樹(shù)和八叉樹(shù)對(duì)數(shù)字巖心進(jìn)行分塊，首先計(jì)算樹(shù)的層級(jí)數(shù)，對(duì)于二維數(shù)字巖心數(shù)據(jù)，四叉樹(shù)的層數(shù)需同時(shí)滿足關(guān)系式如下。

2(n-1)Xblock≥Xfile

(1)

2(n-1)Yblock≥Yfile

(2)

對(duì)于三維數(shù)字巖心數(shù)據(jù)，八叉樹(shù)的層數(shù)還需滿足關(guān)系式為

2(n-1)Zblock≥Zfile

(3)

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

滿足這些不等式的最小層數(shù)作為樹(shù)的層級(jí)數(shù)，依此計(jì)算得到本文選取的二維數(shù)字巖心數(shù)據(jù)四叉樹(shù)的層數(shù)為8，對(duì)應(yīng)子塊的總個(gè)數(shù)為16 384個(gè)，三維數(shù)字巖心數(shù)據(jù)八叉樹(shù)的層數(shù)為5，按照八叉樹(shù)規(guī)則得到對(duì)應(yīng)子塊的總個(gè)數(shù)為4 096個(gè)?；谒牟鏄?shù)和八叉樹(shù)編碼規(guī)則對(duì)子塊賦予編號(hào)[12-13]，如圖1所示，為八叉樹(shù)空間結(jié)構(gòu)劃分及線性莫頓編碼示意圖，按照?qǐng)D中左側(cè)所示對(duì)三維數(shù)字巖心進(jìn)行切割，每個(gè)子塊的編碼如圖1中右側(cè)所示。二維數(shù)字巖心按照類似方式基于四叉樹(shù)進(jìn)行切割和編碼。

圖1 八叉樹(shù)空間結(jié)構(gòu)劃分及線性莫頓編碼Fig.1 Octree spatial structure partition and linear Morton coding

接下來(lái)按子塊編號(hào)的順序依次讀取對(duì)應(yīng)位置的數(shù)字巖心數(shù)據(jù)，對(duì)于靠近邊界區(qū)域的子塊，其長(zhǎng)寬高可能達(dá)不到設(shè)定的子塊大小，截取到邊界即可，對(duì)于邊界外的子塊，生成空子塊。

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

1.3 分塊合并

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

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

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

M=3Mblocklevel

(4)

式(4)中：level為樹(shù)層級(jí)數(shù)；Mblock為單個(gè)子塊占用內(nèi)存大小。

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

1.4 讀取粒度與區(qū)域內(nèi)分塊選擇

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

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

對(duì)于數(shù)字巖心圖像粗細(xì)粒度的選擇，即四叉樹(shù)或八叉樹(shù)級(jí)別的選擇，二維三維數(shù)字巖心所取級(jí)別需分別滿足以下關(guān)系式。

Mapfile/(Pixblock4lev-1)

(5)

Mapfile/(Pixblock8lev-1)

(6)

式中：Mapfile為文件地理坐標(biāo)范圍，式(5)中為面積范圍，式(6)中為體積范圍，在導(dǎo)入時(shí)設(shè)定，默認(rèn)取值為其像素坐標(biāo)范圍；Pixblock在式(5)中為子塊像素面積，在式(6)中為體積；lev為所取級(jí)別；Mapscr和Pixscr分別為屏幕地理坐標(biāo)范圍和像素范圍。式(5)和式(6)分別對(duì)應(yīng)二維和三維數(shù)字巖心。

對(duì)于局部區(qū)域巖心數(shù)據(jù)查看，根據(jù)屏幕所顯示的地理坐標(biāo)范圍，以及數(shù)字巖心的地理坐標(biāo)范圍的位置區(qū)域，得到對(duì)應(yīng)要查看的巖心區(qū)域，讀取區(qū)域內(nèi)子塊數(shù)據(jù)，子塊層級(jí)選擇按前面所述。

1.5 存取效率及粒度選擇實(shí)驗(yàn)對(duì)比

實(shí)驗(yàn)運(yùn)行的環(huán)境為Windows 10，64位操作系統(tǒng)，4核8 GB內(nèi)存，Intel(R) Core(TM) i7-10510U CPU @ 1.80 GHz 2.30 GHz。對(duì)1.1節(jié)中所述的兩個(gè)巖心數(shù)據(jù)進(jìn)行分布式存儲(chǔ)。首先搭建分布式服務(wù)器，如圖4所示，建立GFS Master服務(wù)器負(fù)責(zé)子塊的切分和獲?。唤FS Chunk存儲(chǔ)服務(wù)器來(lái)存放子塊數(shù)據(jù)，最后共建立了3個(gè)存儲(chǔ)服務(wù)器100個(gè)虛擬存儲(chǔ)節(jié)點(diǎn)來(lái)實(shí)現(xiàn)分布式存儲(chǔ)。

圖4 分布式服務(wù)器架構(gòu)Fig.4 Distributed server architecture

按照1.1節(jié)的切塊配置進(jìn)行切塊，分布式存儲(chǔ)到各個(gè)節(jié)點(diǎn)中，各數(shù)據(jù)的切塊結(jié)果如表1所示。

表1 實(shí)驗(yàn)數(shù)據(jù)切塊結(jié)果Table 1 Experimental data segmentation results

與本地文件整體存儲(chǔ)相比，分布式分塊存儲(chǔ)一方面減少本地存儲(chǔ)空間開(kāi)銷(xiāo)，另一方面在訪問(wèn)局部區(qū)域數(shù)據(jù)時(shí)只需要讀取對(duì)應(yīng)子塊，減少文件讀寫(xiě)開(kāi)銷(xiāo)。表2對(duì)比了整體與分塊存儲(chǔ)在訪問(wèn)一個(gè)子塊區(qū)域?qū)?yīng)的數(shù)據(jù)時(shí)，文件讀寫(xiě)的跨度。

由表2可以看出，對(duì)于整體存儲(chǔ)模式數(shù)據(jù)連續(xù)存儲(chǔ)，在訪問(wèn)子塊區(qū)域數(shù)據(jù)時(shí)，為讀取子塊區(qū)域?qū)?yīng)的數(shù)據(jù)，那么會(huì)不斷跨整個(gè)寬度和深度數(shù)據(jù)去訪問(wèn)下一個(gè)位置的數(shù)據(jù)，造成了文件訪問(wèn)跨度的增加。而分塊存儲(chǔ)只需要訪問(wèn)對(duì)應(yīng)子塊文件，不需要額外的跨度。

天津市高等學(xué)校師資培訓(xùn)中心（以下簡(jiǎn)稱天津市中心）成立于1990年，受天津市教委和天津師范大學(xué)雙重領(lǐng)導(dǎo)。中心下設(shè)行政辦公室、崗前培訓(xùn)辦公室、信息技術(shù)部、教師資格認(rèn)定辦公室，當(dāng)前業(yè)務(wù)以崗前培訓(xùn)、教師資格認(rèn)定、高校教師網(wǎng)絡(luò)培訓(xùn)為主。

表2 訪問(wèn)區(qū)域數(shù)據(jù)時(shí)文件讀寫(xiě)跨度對(duì)比Table 2 Comparison of file read and write span when accessing area data

在訪問(wèn)屏幕范圍內(nèi)數(shù)據(jù)時(shí)，會(huì)選擇合適的粗細(xì)粒度子塊數(shù)據(jù)進(jìn)行讀取，表3顯示了與選擇原始數(shù)據(jù)讀取相比，選擇適合粒度層級(jí)進(jìn)行讀取在數(shù)據(jù)量上的對(duì)比，選擇的屏幕分辨率為1 024×1 024。

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

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

2 數(shù)字巖心連通區(qū)域識(shí)別

數(shù)字巖心的孔隙分析、孔喉結(jié)構(gòu)模型計(jì)算、孔喉連通區(qū)域計(jì)算等是數(shù)字巖心分析的重要方面?？紫兜淖R(shí)別通過(guò)設(shè)定的像素閾值進(jìn)行；對(duì)孔喉結(jié)構(gòu)連通性的計(jì)算，首先基于掃描線方法，先計(jì)算連通線，再將相鄰連通線合并成連通區(qū)域，最后將各分塊連通區(qū)域合并成整體連通區(qū)域。

2.1 連通線掃描

首先進(jìn)行掃描參數(shù)的設(shè)定，包括閾值、掃描線最小長(zhǎng)度、連通區(qū)域面積或體積的下限等。依次讀取最細(xì)層級(jí)的每個(gè)子塊進(jìn)行掃描，在掃描過(guò)程中，二維按照X、Y的順序，三維按照X、Y、Z的順序依次讀取子塊的每個(gè)點(diǎn)的像素值，對(duì)于滿足閾值的像素點(diǎn)，相鄰的連接起來(lái)形成線段，即掃描線，如果掃描線長(zhǎng)度小于最小長(zhǎng)度則剔除掉。

在例子中，二維數(shù)字巖心閾值設(shè)定為55，掃描線最小長(zhǎng)度為3，連通區(qū)域面積下限為22，三維數(shù)字巖心閾值設(shè)定為55，掃描線最小長(zhǎng)度為3，連通區(qū)域體積下限為50。即掃描二維數(shù)字巖心過(guò)程中，像素值小于55的點(diǎn)會(huì)記錄下來(lái)，連續(xù)3個(gè)或以上這種點(diǎn)可以組成一條掃描線，若干個(gè)相鄰的掃描線組成的連通區(qū)域總面積還要大于等于22；掃描三維數(shù)字巖心過(guò)程中，像素值小于55的點(diǎn)會(huì)記錄下來(lái)，連續(xù)3個(gè)或以上這種點(diǎn)可以組成一條掃描線，若干個(gè)相鄰的掃描線組成的連通區(qū)域總體積還要大于等于50。

2.2 分塊連通區(qū)域生成

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

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

2.3 合并相鄰分塊連通區(qū)域

將相鄰分塊的連通區(qū)域進(jìn)行合并，對(duì)于二維數(shù)字巖心，相鄰分塊為上下左右4個(gè)，對(duì)于三維數(shù)字巖心，相鄰分塊為上下左右前后6個(gè)。對(duì)于連通區(qū)域存在相鄰的兩個(gè)分塊，將它們的連通區(qū)域掃描線進(jìn)行合并得到新的連通區(qū)域。

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

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

2.4 合并整體連通區(qū)域

3 連通區(qū)域分色及邊界面計(jì)算

為方便查看各個(gè)孔喉連通區(qū)域形態(tài)，需要對(duì)其進(jìn)行染色，基于三角剖分和四面體剖分對(duì)二維和三維連通區(qū)域分色，并結(jié)合地理信息技術(shù)將各連通區(qū)域作為圖元存儲(chǔ)形成圖層，圖元屬性包括顏色、面積或體積等，方便查詢和管理?；谝苿?dòng)立方體(marching cube，MC)算法[14]對(duì)三維數(shù)字巖心模型孔隙求取邊界面，以便在三維空間中展示孔喉結(jié)構(gòu)。

3.1 三角和四面體剖分連通區(qū)域及分色

對(duì)各個(gè)連通區(qū)域的中點(diǎn)進(jìn)行三角/四面體剖分連接起來(lái)。二維采用三角剖分算法，三維采用四面體剖分算法[15]。如圖7所示為例子中二維數(shù)字巖心連通區(qū)域進(jìn)行三角剖分后的結(jié)果。

圖7 二維數(shù)字巖心三角剖分結(jié)果Fig.7 Two dimensional digital core triangulation results

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

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

按照地圖染色算法，統(tǒng)計(jì)各個(gè)點(diǎn)連接的邊數(shù)量，從大到小排序，從邊數(shù)量最多的點(diǎn)開(kāi)始染色，每染色一個(gè)點(diǎn)，將其四周有邊相連的點(diǎn)的邊計(jì)數(shù)減去1，然后將剩余點(diǎn)按照所連接的邊的數(shù)量從大到小排序，繼續(xù)從最多邊數(shù)量點(diǎn)染色，每次染色要與周?chē)羞呄噙B的點(diǎn)的顏色不同。然后依此類推循環(huán)下去，直到所有點(diǎn)都被染色。

接下來(lái)對(duì)于二維數(shù)字巖心，將染色結(jié)果用地理信息圖層表示。每個(gè)連通區(qū)域作為圖層中一個(gè)圖元對(duì)象，并賦予所染顏色。圖元對(duì)應(yīng)屬性存放這個(gè)連通區(qū)域的中點(diǎn)坐標(biāo)、面積、范圍、掃描線集合等。在例子中，二維數(shù)字巖心生成的地理信息圖層如圖9所示。

圖9 二維數(shù)字巖心連通區(qū)域識(shí)別Fig.9 Recognition of connected area of two-dimensional digital core

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

圖10 三維數(shù)字巖心連通區(qū)域識(shí)別Fig.10 Recognition of connected area of 3D digital core

3.2 連通區(qū)域邊界面計(jì)算

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

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

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

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

(7)

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

對(duì)于立方體的8個(gè)頂點(diǎn)，01分界面共有256種情況，可以分為15類，如圖12所示，列出了不同情況的等值面構(gòu)建方法。將每個(gè)立方體形成的等值面連接起來(lái)，即是孔喉的分界面。連接時(shí)，與連通區(qū)域分組類似，將具有相連關(guān)系的分界面劃分為一組，屬于同一個(gè)連通孔喉，互不相連的分界面屬于不同孔喉。

圖12 移動(dòng)立方體等勢(shì)面構(gòu)建情況組合Fig.12 Combination of construction of equipotential surface of moving cube

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

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

4 結(jié)論

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

(2)通過(guò)掃描線合并構(gòu)建連通區(qū)域，基于移動(dòng)立方體方法獲取孔隙分界面，得到了孔喉結(jié)構(gòu)的形態(tài)及連通性質(zhì)。通過(guò)三角剖分和四面體剖分對(duì)連通區(qū)域分色，方便區(qū)分顯示各連通區(qū)域。