張 磊，李昊勃，劉天宇，全 馨，樊冬艷，孫 海，付帥師，楊永飛

基于C#語言的數(shù)字巖心數(shù)值實(shí)驗(yàn)演示系統(tǒng)

張 磊1，李昊勃2，劉天宇2，全 馨2，樊冬艷1，孫 海1，付帥師1，楊永飛1

（1. 中國石油大學(xué)（華東） 石油工程學(xué)院，山東 青島 266580；2. 中國石油大學(xué)（華東） 儲運(yùn)與建筑工程學(xué)院，山東 青島 266580）

基于Visual Studio集成開發(fā)環(huán)境和C#語言，開發(fā)了數(shù)字巖心數(shù)值實(shí)驗(yàn)演示系統(tǒng)。該系統(tǒng)具有數(shù)字巖心掃描圖像的輸入輸出、二值化處理、分形分析、孔徑分布統(tǒng)計(jì)、數(shù)據(jù)保存等功能。系統(tǒng)的軟件設(shè)計(jì)采用大量類實(shí)現(xiàn)，在具體功能實(shí)現(xiàn)上只需調(diào)用相應(yīng)類方法，便于系統(tǒng)功能的增加和擴(kuò)展。數(shù)字巖心圖像的二值化采用Otsu算法以及改進(jìn)的分形算法，可根據(jù)閾值進(jìn)行二值化處理，通過調(diào)整二值化閾值觀察數(shù)字巖心的孔徑分布曲線的變化。進(jìn)行二值化處理和邊緣識別后的圖像與原圖并列顯示以進(jìn)行比較，通過調(diào)整二值化閾值觀察孔徑的變化。

數(shù)字巖心；實(shí)驗(yàn)教學(xué)；演示實(shí)驗(yàn)；C#語言

巖石物性是“油層物理”課程中非常重要的內(nèi)容，巖石的很多物理屬性（如滲透率、孔隙度、毛管力等）均取決于它的微觀結(jié)構(gòu)和流體的性質(zhì)，其中巖石的微觀結(jié)構(gòu)是基礎(chǔ)。研究巖石的微觀結(jié)構(gòu)至關(guān)重要。

數(shù)字巖心[1]是通過CT掃描等技術(shù)得到的油藏儲層巖心樣品的數(shù)字化圖像。在數(shù)字巖心技術(shù)出現(xiàn)以前，巖石的物理屬性參數(shù)通常由實(shí)驗(yàn)室物理實(shí)驗(yàn)得到，巖心內(nèi)部結(jié)構(gòu)是一個(gè)“黑箱”。數(shù)字巖心技術(shù)成功地把“黑箱”變?yōu)椤鞍紫洹?，基于?shù)字巖心的計(jì)算機(jī)模擬技術(shù)從而替代常規(guī)物理實(shí)驗(yàn)，成為研究巖石物性的一種新技術(shù)。筆者基于Visual Studio集成開發(fā)環(huán)境，采用C#語言[2]設(shè)計(jì)開發(fā)了數(shù)字巖心數(shù)值實(shí)驗(yàn)演示系統(tǒng)，并將該系統(tǒng)應(yīng)用于教學(xué)中，可以直觀地展示不同處理算法的效果，方便進(jìn)行對比分析，加深學(xué)生對巖石孔隙結(jié)構(gòu)的認(rèn)識，增加學(xué)習(xí)的趣味性。該系統(tǒng)采用類設(shè)計(jì)，學(xué)生可以在此基礎(chǔ)上進(jìn)行功能擴(kuò)展，進(jìn)行創(chuàng)新項(xiàng)目的探索，提高創(chuàng)新能力。

1 軟件設(shè)計(jì)

數(shù)字巖心技術(shù)的關(guān)鍵是借助于計(jì)算機(jī)圖像處理技術(shù)[3]，對巖石的掃描圖像進(jìn)行處理，分析其微觀孔隙結(jié)構(gòu)特征，獲取反映孔隙空間分布、連通性等的重要信息。在這一個(gè)過程中，涉及一系列的數(shù)字巖心掃描圖像的處理過程。基于Visual Studio集成開發(fā)環(huán)境，采用C#語言編程，能夠使開發(fā)人員準(zhǔn)確、高效地編寫代碼，并且不會丟失當(dāng)前的文件上下文，可以輕松地放大到詳細(xì)信息，例如調(diào)用結(jié)構(gòu)、相關(guān)函數(shù)、簽入和測試狀態(tài)，還可以方便地進(jìn)行重構(gòu)、識別和修復(fù)代碼。.NET框架為C#提供了一個(gè)強(qiáng)大、易用、邏輯結(jié)構(gòu)一致的設(shè)計(jì)環(huán)境。

數(shù)字巖心數(shù)值實(shí)驗(yàn)演示系統(tǒng)包括數(shù)字巖心圖像的讀入、二值化處理（自動化處理和手動處理、分形自動處理）、數(shù)據(jù)轉(zhuǎn)化、分形參數(shù)分析、孔隙尺寸分析和結(jié)果保存等主要功能，并從數(shù)字巖心角度分為二維、三維數(shù)字巖心的數(shù)值實(shí)驗(yàn)。

數(shù)字巖心數(shù)值實(shí)驗(yàn)演示系統(tǒng)的用戶界面如圖1所示。根據(jù)數(shù)字巖心對象不同，該系統(tǒng)分為不同標(biāo)簽頁，每一個(gè)標(biāo)簽頁的數(shù)據(jù)相互獨(dú)立，每個(gè)標(biāo)簽頁左側(cè)為功能按鈕，右側(cè)上方為原始圖像和處理結(jié)果的對比展示，右側(cè)下方為圖像數(shù)據(jù)信息以及相關(guān)曲線的展示。

圖1 系統(tǒng)操作界面

軟件設(shè)計(jì)采用數(shù)據(jù)類的處理方式，通過按鈕是否可用來實(shí)現(xiàn)數(shù)字巖心圖像處理的流程控制。每個(gè)標(biāo)簽頁中的數(shù)據(jù)雖然相互獨(dú)立，但由于采用了數(shù)據(jù)類的處理，相同的功能不存在代碼的重復(fù)，都是通過實(shí)例化同一個(gè)類，然后調(diào)用其中的成員函數(shù)來實(shí)現(xiàn)。

由于不同的功能之間存在先后關(guān)系，例如只有當(dāng)完成了二值化處理以后才能進(jìn)行圖像的導(dǎo)出或者數(shù)據(jù)的存儲，因此在軟件設(shè)計(jì)中通過控件的Enable屬性來控制每個(gè)功能的可用性，而且在每一個(gè)功能執(zhí)行完成后，重新刷新每一個(gè)按鈕控件的Enable屬性，根據(jù)不同功能之間的先后關(guān)系，確定按鈕是否可用，從而保證數(shù)字巖心數(shù)值實(shí)驗(yàn)演示系統(tǒng)的穩(wěn)健性。

2 實(shí)例分析與演示

根據(jù)功能劃分，數(shù)字巖心數(shù)值實(shí)驗(yàn)演示系統(tǒng)分為輸入模塊、數(shù)值實(shí)驗(yàn)?zāi)K、分析模塊、曲線模塊和輸出功能模塊，各模塊的功能如圖2所示。輸入模塊包括讀原始圖像和優(yōu)化參數(shù)的輸入；數(shù)值實(shí)驗(yàn)?zāi)K包括二值化處理的3種方法和邊緣提??；分析模塊包括分形維數(shù)、孔隙尺寸、迂曲度和孔隙度的分析；曲線模塊包括灰度和孔隙尺寸分布統(tǒng)計(jì)曲線；輸出功能包括數(shù)據(jù)、圖像和曲線的存儲。

數(shù)字巖心掃描圖像一般為灰度圖像，像素點(diǎn)的灰度值范圍為0—255，二值化處理過程就是選擇合適的灰度值（閾值），將圖像像素灰度值根據(jù)大于或者小于該閾值，分成孔隙和骨架兩類。由于巖石的非均質(zhì)性和儀器掃描時(shí)受到光線強(qiáng)弱等因素影響，得到的巖石掃描圖像會存在一些誤差。數(shù)字巖心的二值化處理是數(shù)字巖心模擬中一個(gè)最基礎(chǔ)的過程，因此對掃描圖像進(jìn)行合理、正確的分割至關(guān)重要。

圖2 數(shù)字巖心數(shù)值實(shí)驗(yàn)巖石系統(tǒng)功能框架圖

分割算法有迭代法[4]、簡單統(tǒng)計(jì)法[5]、Otsu方法[6]以及改進(jìn)的Otsu算法[7-8]，本演示系統(tǒng)涉及的算法為Otsu方法和基于Otsu算法改進(jìn)的分形分割算法[8]。其中二值化手動處理方法即給定某一閾值進(jìn)行分割，本文不再贅述。

2.1 Otsu方法原理

Otsu方法是學(xué)者Otsu于1979年提出的最大類間方差法。該方法將像素點(diǎn)按灰度值分割成兩組，一組灰度對應(yīng)孔隙，另一組灰度對應(yīng)骨架，目標(biāo)是使兩組灰度值的類內(nèi)方差最小、類間方差最大；而方差值越大，說明構(gòu)成圖像的兩部分差別越大。

2.2 分形分割算法原理

分形分割算法把Otsu方法和分形理論結(jié)合起來，通過分析巖心孔隙結(jié)構(gòu)的分形特征[9-10]改進(jìn)目標(biāo)函數(shù)，有效提高CT圖像分割的準(zhǔn)確性。假設(shè)以某閾值對數(shù)字巖心掃描圖像進(jìn)行分割，通過統(tǒng)計(jì)孔隙數(shù)量可以得到孔隙度()，通過分析得到最大孔隙尺寸和最小孔隙尺寸分別為max和min，則有孔隙度p和孔隙尺寸分形維數(shù)f以及最大、最小孔隙尺寸的關(guān)系[11]為

式中：e為歐式維數(shù)，對于二維和三維空間分別為e＝2和e＝3；f是分形維數(shù)。對于二維空間，1

結(jié)合Otsu算法的目標(biāo)函數(shù)()，定義新目標(biāo)函數(shù)()：

遍歷所有灰度值，當(dāng)目標(biāo)函數(shù)達(dá)到最大值時(shí)，所對應(yīng)的灰度值即為分割閾值。

2.3 分形維數(shù)計(jì)算原理

數(shù)字巖心數(shù)值實(shí)驗(yàn)演示系統(tǒng)采用計(jì)盒維數(shù)法[12-13]計(jì)算數(shù)字巖心孔隙或者骨架的分形維數(shù)。分形維數(shù)的計(jì)算在二值化處理完成以后的黑白圖像的基礎(chǔ)上進(jìn)行，具體計(jì)算步驟如下：

（1）取邊長為2個(gè)像素點(diǎn)的正方形邊框覆蓋二值圖像，統(tǒng)計(jì)覆蓋到黑色像素點(diǎn)的盒子數(shù)()；

（2）邊長依次選取4、6、8個(gè)像素點(diǎn)，直到邊長等于圖像寬度的一半，統(tǒng)計(jì)對應(yīng)邊長下的盒子數(shù)；

（3）利用最小二乘法線性擬合lg和lg()，擬合所得直線的斜率記為分形維數(shù)。

2.4 結(jié)果演示

以人造巖心為例進(jìn)行展示。圖3為巖心經(jīng)過CT掃描后得到的灰度圖像。從圖像上可以看到巖心中存在不同尺寸的顆粒以及膠結(jié)物；圖4為數(shù)字巖心掃描圖像灰度分布圖，可以看出圖像具有明顯的雙峰分布，說明圖像中孔隙（黑色）和骨架（白色）的區(qū)分還是比較明顯的。

采用Otsu算法和改進(jìn)的分形算法分別進(jìn)行二值化處理，處理結(jié)果如圖5所示。可以看到：改進(jìn)的分形算法處理結(jié)果更接近真實(shí)情況，能夠準(zhǔn)確地分辨出數(shù)字巖心圖像中的固體顆粒。在此基礎(chǔ)上，可以對巖心進(jìn)行孔隙尺寸分析。圖6為巖石的孔隙尺寸分布曲線圖，可以看出，巖石孔隙主要寬度為5個(gè)像素。

圖3 數(shù)字巖心原始掃描圖像

圖4 原始圖像灰度分布圖

圖5 二值化處理結(jié)果

圖6 孔隙分布曲線結(jié)果

通過參數(shù)分析結(jié)果，可以得到處理過程中的所有參數(shù)設(shè)置，以及通過數(shù)字巖心分析得到的巖石物性參數(shù)，從而進(jìn)一步了解巖石的物理性質(zhì)。

3 結(jié)語

基于Visual Studio集成開發(fā)環(huán)境開發(fā)的數(shù)字巖心數(shù)值實(shí)驗(yàn)演示系統(tǒng)，實(shí)現(xiàn)了輸入功能、數(shù)值實(shí)驗(yàn)、分析功能和輸出功能等，具有良好的用戶交互界面，方便對數(shù)字巖心的掃描圖像進(jìn)行二值化處理、分形分析等一系列操作，最終得到巖心的孔隙度、孔徑分布等物理屬性，對處理結(jié)果圖像和相關(guān)曲線進(jìn)行了實(shí)時(shí)顯示。該系統(tǒng)在油層物理巖石物理屬性的教學(xué)中達(dá)到了以下教學(xué)目的：

（1）掌握數(shù)字巖心掃描圖像二值化處理的基本原理和方法，原始圖像和結(jié)果圖像的實(shí)時(shí)對比演示展現(xiàn)了不同算法的實(shí)際效果，同時(shí)給出孔隙度等參數(shù)，使學(xué)生對算法和巖石的基本性質(zhì)加深理解；

（2）通過對處理結(jié)果圖像的分析，給出孔隙半徑分布曲線，有助于加深對巖石孔隙結(jié)構(gòu)的認(rèn)識

（3）通過對二值化圖像的分形分析，使學(xué)生對巖石孔隙結(jié)構(gòu)的分形特征有更深入了解。

將該演示系統(tǒng)應(yīng)用于油層物理的課程教學(xué)和實(shí)驗(yàn)教學(xué)中，對學(xué)生深入理解基礎(chǔ)理論知識和掌握相應(yīng)的實(shí)驗(yàn)方法有積極的作用，可以激發(fā)學(xué)生的科研興趣和創(chuàng)新意識，提高學(xué)生的創(chuàng)新能力。目前，在此教學(xué)演示系統(tǒng)的基礎(chǔ)上，已培養(yǎng)我校校級大學(xué)生創(chuàng)新項(xiàng)目兩項(xiàng)，取得了較好的實(shí)驗(yàn)教學(xué)效果。

[1] 姚軍，趙秀才.數(shù)字巖心及孔隙級滲流模擬理論[M].北京：石油工業(yè)出版社，2010.

[2] 匡松，張淮鑫. C#開發(fā)寶典[M].北京：中國鐵道出版社，2010.

[3] 岡薩雷斯.數(shù)字圖像處理[M].北京：電子工業(yè)出版社，2011: 62–190.

[4] PEREZ A, GONZALEZ R C. An iterative thresholding algori-thm for image segmentation[J]. IEEE transactions on pattern analysis and machine intelligence, 1987, 9(6): 742–751.

[5] KITTLER J, ILLING W J, FOGLEIN J. Threshold selection ba-s---ed on a simple image static[J]. Computer Vision, Graphics and Image Procession, 1985, 30: 125–147.

[6] OTSU N. A threshold selection method from gray-level histog-ram[J]. IEEE Trans, 1979, 9(1): 62–66.

[7] 趙秀才，姚軍，房克榮.合理分割巖心微觀結(jié)構(gòu)圖像的新方法[J].中國石油大學(xué)學(xué)報(bào)（自然科學(xué)版），2009, 33(1): 64–67, 72.

[8] 張磊，姚軍，王鋒，等.考慮分形特征的巖心掃描圖像合理分割方法[J].科學(xué)技術(shù)與工程，2015, 15(24): 57–60.

[9] 趙明，郁伯銘.數(shù)字巖心孔隙結(jié)構(gòu)的分形表征及滲透率預(yù)測[J].重慶大學(xué)學(xué)報(bào)，2011, 34(4): 88–94.

[10] 郁伯銘.多孔介質(zhì)輸運(yùn)性質(zhì)的分形分析研究進(jìn)展[J].力學(xué)進(jìn)展，2003, 33(3): 333–346.

[11] YU Boming, LI Jianhua. Some Fractal Characters of Porous Media[J]. Fractals, 2001, 9(3): 365–372.

[12] 張懷亮，孫栓輝，劉舜堯.數(shù)字分形圖像分形維數(shù)計(jì)算[J].計(jì)算機(jī)仿真，2010, 27(10): 215–217, 254.

[13] 彭瑞東，楊彥從，鞠楊，等.基于灰度CT圖像的巖石孔隙分形維數(shù)計(jì)算[J].科學(xué)通報(bào)，2011, 56(26): 2256–2266.

Demonstration system of digital core numerical experiment based on C# language

ZHANG Lei1, LI Haobo2, LIU Tianyu2, QUAN Xin2, FAN Dongyan1, SUN Hai1, FU Shuaishi1, YANG Yongfei1

(1. School of Petroleum Engineering, China University of Petroleum, Qingdao 266580, China; 2. College of Storage, Transportation and Architectural Engineering, China University of Petroleum, Qingdao 266580, China)

Based on the Visual Studio integrated development environment and C# language, a demonstration system of the digital core numerical experiment is developed. This system realizes the functions of the input and output, binary processing, fractal analysis, pore size distribution statistics, data preservation of digital core scanning image, etc. This system uses a large number of class realizations in software design and only calls the corresponding class methods in the realization of the specific function, which is convenient for the increase and expansion of system functions. The Otsu algorithm and improved fractal algorithm are used for binarization of digital core image. According to the threshold value, binarization can be carried out, and the change of pore size distribution curve of digital core can be observed by adjusting the threshold value. The image after binarization and edge recognition is displayed side by side with the original image for comparison.

digital core; experimental teaching; demonstration experiment; C# language

TE319

A

1002-4956(2019)10-0127-04

10.16791/j.cnki.sjg.2019.10.030

2019-03-18

國家自然科學(xué)基金面上項(xiàng)目（51674280）；中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)資金項(xiàng)目（18CX02031A）；中國石油大學(xué)（華東）教學(xué)改革項(xiàng)目（QN201706，JY-B201812）；山東省本科教改項(xiàng)目重點(diǎn)項(xiàng)目（Z2018B118），山東省高等學(xué)校教學(xué)改革立項(xiàng)項(xiàng)目（B2016M008）

張磊（1984—），男，山東濱州，博士，講師，主要研究方向?yàn)閿?shù)字巖心技術(shù)及微觀流動模擬。E-mail: zhlei84@163.com