胡 渤，蒲 軍，茍斐斐

（1.中國石化油田勘探開發(fā)事業(yè)部，北京 100728；2.中國石化石油勘探開發(fā)研究院，北京 101599）

中國致密砂巖油藏資源豐富，初步評估可采儲量達35×108～40×108t［1-6］，約占世界致密油總可采儲量的十分之一，是未來現(xiàn)實的資源接替類型。致密砂巖儲層多發(fā)育微米級孔隙和納米級喉道，孔喉半徑級差大、形態(tài)不規(guī)則、結(jié)構(gòu)復(fù)雜、微觀非均質(zhì)性強，極大地影響了儲層的滲流能力［7-13］，定量描述微觀孔喉結(jié)構(gòu)特征可以為致密砂巖儲層質(zhì)量差異化表征提供依據(jù)，也可以為致密砂巖儲層分類及開發(fā)對策研究奠定基礎(chǔ)，具有重要的意義。但常規(guī)實驗方法（壓汞、核磁、驅(qū)替等）存在實驗周期長、精度低等問題［14-16］，難以準確地確定微觀孔喉結(jié)構(gòu)參數(shù)，制約了致密砂巖孔喉結(jié)構(gòu)特征研究。

近年來，隨著掃描成像設(shè)備與計算機技術(shù)的快速發(fā)展，基于高精度掃描圖像構(gòu)建數(shù)字巖心，建立多孔介質(zhì)模型，開展統(tǒng)計分析，為研究致密砂巖儲層孔喉結(jié)構(gòu)特征提供了新的技術(shù)手段［17-19］。目前業(yè)界內(nèi)主要采用孔隙網(wǎng)絡(luò)模型方法來分割孔喉并計算孔喉特征參數(shù)。該方法是將巖石孔喉空間簡化為球棍模型，用圓球代表孔隙，孔隙之間用細棍代表喉道。常用的孔隙網(wǎng)絡(luò)模型方法主要有最大球法［20-21］和中軸線法［22-26］。最大球法需要計算每個像素的最大球，最終生成的孔隙網(wǎng)絡(luò)模型是對真實孔喉的一種等效假設(shè)，計算量較大，孔喉定量表征精度低。中軸線法受孔喉形態(tài)影響大，致密砂巖孔喉空間復(fù)雜多變，提取的中軸線會產(chǎn)生大量分支，也影響了孔喉分割與表征的精度。

直接在數(shù)字巖心圖像的基礎(chǔ)上開展孔喉結(jié)構(gòu)定量表征，無需對孔喉空間進行簡化處理，可以最大限度地保證孔喉表征的精度與可靠性，成為當(dāng)前研究的熱點與前沿，但目前仍然存在孔喉提取效率低、精細分割難度大等問題。筆者從圖像形態(tài)學(xué)出發(fā)（以SEM-Maps 掃描圖像為例，分辨率為10 nm），開發(fā)了自適應(yīng)孔喉識別與提取方法，研發(fā)了直接基于圖像形態(tài)學(xué)算法的孔喉分割方法及孔喉結(jié)構(gòu)特征參數(shù)計算方法，并以紅河油田長8 儲層典型井致密砂巖巖樣為例，分析了不同儲層的微觀孔喉結(jié)構(gòu)及連通性差異，研究成果可以為致密砂巖儲層質(zhì)量表征提供依據(jù)和借鑒。

1 致密砂巖數(shù)字巖心圖像處理

數(shù)字巖心圖像（包括SEM-Maps 和CT 等）的灰度值為0～255，灰度值較高的像素反映的是相對密度和原子序數(shù)較大的物質(zhì)，如巖石骨架；灰度值較低的像素反映的是相對密度和原子序數(shù)較低的物質(zhì)，如孔喉空間［9-10］。數(shù)字巖心圖像處理的目的是提高圖像清晰度并提取孔喉空間，為下一步孔喉結(jié)構(gòu)定量表征做好準備。

1.1 圖像降噪

數(shù)字巖心圖像在數(shù)字化和傳輸過程中常受到成像設(shè)備與外部環(huán)境噪聲的影響，圖像中含有大量噪點，如圖1所示，這些噪點極大地影響了孔喉識別與提取的準確性［11-12］。因此，孔喉空間識別與提取的首要步驟是進行降噪處理。

圖1 數(shù)字巖心圖像噪點Fig. 1 Noises in digital core images

系統(tǒng)分析數(shù)字巖心圖像噪點分布特征（以SEM-Maps 掃描圖像為例），發(fā)現(xiàn)圖像中的噪點主要為高斯噪點（圖2），且不同礦物顆粒和孔喉空間具有不同的噪點標(biāo)準差。圖1 中1—4 位置處的噪點標(biāo)準差分別為7.82，8.48，0 和6.87，選擇合適的降噪算法有效去除噪點是獲取高質(zhì)量數(shù)字巖心圖像的關(guān)鍵。

圖2 SEM-Maps掃描圖像噪點分布（圖1位置1處）Fig. 2 Distribution of noises in SEM-Maps scanning image（location1 on Fig.1）

雙邊降噪算法是在高斯降噪算法中加入圖像灰度權(quán)重項，不但可以有效去除高斯噪點，同時能夠保持孔隙和骨架邊界清晰［13］，成為數(shù)字巖心圖像降噪的首選算法。

設(shè)f(i,j)表示原圖像灰度值，g(i,j)表示降噪后圖像灰度值，則雙邊降噪公式為：

根據(jù)圖像噪點統(tǒng)計結(jié)果，取統(tǒng)計標(biāo)準差的最大值σd=σr=8.48。雙邊降噪前后數(shù)字巖心圖像如圖3所示，礦物顆粒和孔喉空間的噪點大幅度減少，相同灰度閾值下提取的孔隙更加準確（灰度閾值為150，紅色為孔喉）。

圖3 降噪前后數(shù)字巖心圖像孔喉空間提取結(jié)果對比Fig. 3 Comparison of extraction results of pore throat spaces in digital images before and after noise reduction（Threshold is150，and the red part is pore throat）

1.2 孔喉提取

消除噪點后，可以對孔喉空間進行識別和提取。數(shù)字巖心圖像孔喉提取的關(guān)鍵在于確定合適的灰度分割閾值，常規(guī)方法是通過人工觀察確定分割閾值，但準確性較低。筆者提出了自適應(yīng)孔喉識別與提取方法（圖4），通過對降噪后數(shù)字巖心圖像的灰度值進行統(tǒng)計，結(jié)合大津法，自動計算得到合適的灰度分割閾值，避免了人為因素的干擾。

圖4 自適應(yīng)孔喉識別與提取方法流程Fig. 4 Process of adaptive pore throat recognition and extraction

對降噪后數(shù)字巖心圖像中的所有像素點的灰度值進行統(tǒng)計，得到灰度累積分布曲線（圖5）?；叶壤鄯e分布計算方法與面孔率（孔隙度）的計算方法一致，因此，圖5的縱坐標(biāo)可以理解為對應(yīng)灰度值的面孔率（孔隙度）。

圖5 數(shù)字巖心圖像灰度與孔隙度的對應(yīng)關(guān)系Fig. 5 Relationship between gray value and porosity in digital core image

致密砂巖巖心孔隙度通常大于4%，通過數(shù)字巖心圖像觀察，結(jié)合巖心孔隙度分析，確定灰度值小于50 的像素點為孔喉空間（圖5）；當(dāng)灰度值大于165 時，孔隙度急劇增大，這是因為灰度值大于165時，大量巖石礦物顆粒被誤識別為孔喉。因此，合適的孔隙/巖石骨架灰度分割閾值為50～165。

為準確地確定灰度分割閾值，可以對數(shù)字巖心圖像進行灰度拉伸：

對灰度拉伸后的圖像（圖6a）采用Ostu 算法［14］，以像素灰度值級差最大為分割依據(jù)，得到孔喉空間與巖石骨架灰度分割閾值界限為141，分割得到的結(jié)果如圖6b 所示，白色為巖石礦物顆粒，黑色為孔喉空間。

圖6 灰度閾值分割后的數(shù)字巖心圖像Fig. 6 Digital core images after grey level threshold segmentation

2 基于數(shù)字巖心圖像的孔喉分割

常規(guī)砂巖孔喉級差小、形態(tài)規(guī)則，采用孔隙網(wǎng)絡(luò)模型分割的孔隙和喉道與實際情況基本相符（圖7）。但致密砂巖巖心孔喉級差大、形態(tài)復(fù)雜多變，采用孔隙網(wǎng)絡(luò)模型易將大孔隙識別為多個小孔隙的組合，且易將喉道中略寬的部分識別成小孔隙，這將極大地影響孔喉定量表征結(jié)果的準確性。

圖7 不同類型砂巖巖心孔隙網(wǎng)絡(luò)模型提取結(jié)果對比Fig. 7 Comparison of extraction results of pore network models in different sandstone cores

從致密砂巖孔喉形態(tài)差異性出發(fā)，基于圖像形態(tài)學(xué)原理，采用腐蝕、膨脹及邏輯算法建立了一套適用于致密砂巖的孔喉分割方法。以致密砂巖數(shù)字巖心典型孔喉圖像（圖8）為例，將孔喉分割的具體流程說明如下。

步驟1：采用腐蝕算子消除喉道。在圖像形態(tài)學(xué)中，設(shè)U為一個平面區(qū)域，設(shè)A1為區(qū)域U內(nèi)部的1個目標(biāo)區(qū)域，S1為指定大小和形狀的結(jié)構(gòu)元素1，定義位于坐標(biāo)(x,y)的S1所表示的區(qū)域為S1(x,y)，則對目標(biāo)區(qū)域A1的腐蝕操作可以表示為：

S1的大小取決于數(shù)字巖心圖像中喉道的大小，通過觀察可以看到：由于喉道可以沿著任意方向延伸，結(jié)構(gòu)元素的形狀通常采用圓盤狀，這樣可以保證所有喉道均可以被腐蝕（圖9）。

圖9 直徑為5個像素的圓盤狀結(jié)構(gòu)元素示意Fig. 9 Disk-shaped structural elements with a diameter of 5 pixels

在圖8 中，A1為原始孔喉，采用S1腐蝕后，不僅所有喉道被腐蝕（圖10a），孔隙周邊像素也會被腐蝕。因此，圖10a中顯示的并非真實孔隙，這里稱之為殘缺孔隙（A2）。

圖8 致密砂巖數(shù)字巖心典型孔喉圖像（A1）Fig. 8 Typical pore throat image of tight sandstone based on digital core technology（A1）

步驟2：采用膨脹算子修復(fù)殘缺孔隙。在圖像形態(tài)學(xué)中，設(shè)U為一個平面區(qū)域，S2為指定大小和形狀的結(jié)構(gòu)元素2，定義位于坐標(biāo)()x,y的S2所表示的區(qū)域為S2()x,y，則對目標(biāo)區(qū)域A2的膨脹運算可以表示為：

對殘缺孔隙（圖10a）進行膨脹運算修復(fù)后，被腐蝕的孔隙邊界像素得到了恢復(fù)（圖10b）；由于膨脹運算并不是腐蝕運算的逆運算，因此這種恢復(fù)并不準確。為了確保膨脹運算后的孔隙能夠恢復(fù)到原始孔隙，膨脹運算的S2的直徑要略大于步驟1 中腐蝕運算的S1，但S2的結(jié)構(gòu)與S1保持一致即可，膨脹運算修復(fù)后的孔隙為A3。

圖10 數(shù)字巖心孔喉圖像腐蝕膨脹處理結(jié)果Fig. 10 Processing results of corrosion and expansion in pore throat images

步驟3：采用交集運算提取原始孔隙。采用膨脹運算修復(fù)后的孔隙A3與原始孔喉A1進行交集運算，即可得到準確的孔隙A4（圖11a）。交集運算可以表示為：

步驟4：采用差集運算提取喉道。分割得到準確的孔隙后，從原始孔喉圖像中用差集運算去除孔隙部分，即可得到準確的喉道A5（圖11b）。差集運算可以定義為：

圖11 基于數(shù)字巖心圖像的孔喉分割結(jié)果Fig. 11 Segmentation results of pore throats based on digital core images

3 孔喉定量表征

孔隙和喉道均為幾何圖形，具有相同的幾何特征，例如面積、周長等。此外，孔隙和喉道在形態(tài)上又具有較大的差異，存在特征參數(shù)，例如形狀因子、配位數(shù)等。

3.1 孔喉共有參數(shù)計算

采用GRAY提出的模板匹配方法來計算孔喉的面積和周長［26］。在結(jié)構(gòu)元素模板中，1 表示目標(biāo)像素，0表示空白像素（圖12）。

圖12 結(jié)構(gòu)元素模板Fig. 12 Template of structural element

采用n{ }Q表示孔喉圖像中匹配的模板圖像Q的個數(shù)。孔喉面積和周長的計算公式分別為：

3.2 孔隙特征參數(shù)計算

3.2.1 形狀因子

形狀因子用來描述孔隙接近圓的程度，即圓的形狀因子為1。由于前文已經(jīng)計算得到孔隙的面積和周長，所以這里用面積和周長來描述孔隙的形狀因子。對于周長為Pp的圓，其面積為：

則可以定義孔隙的形狀因子為：

3.2.2 面孔率

面孔率是指數(shù)字巖心孔隙的面積與總面積的比值。假設(shè)數(shù)字巖心的總像素個數(shù)為N，孔隙的像素個數(shù)為M，則面孔率的計算公式為：

3.3 喉道特征參數(shù)計算

3.3.1 長度和平均寬度

已知喉道的面積At和周長Pt，設(shè)喉道長度為L，平均寬度為d，則有：

求解（15）和（16）式，可得喉道的長度和平均寬度分別為：

3.3.2 配位數(shù)

在儲層孔隙結(jié)構(gòu)表征中，配位數(shù)是指每個孔隙連通喉道的數(shù)量，通過對孔喉圖像進行運算，可以得到孔喉連接關(guān)系矩陣cIJ，配位數(shù)計算公式為：

3.4 孔隙連通性計算

巖心中并非所有孔隙都是相互連通的。根據(jù)孔隙的連通關(guān)系，可以將數(shù)字巖心圖像劃分為不同的連通體，而最大連通體的大小直接決定巖心的滲透率，因此可以最大連通體面積占比來表征孔隙連通性：

4 實例應(yīng)用

基于紅河油田長8 儲層物性差異選取3 口典型井，分別采集巖樣進行SEM-Maps掃描，構(gòu)建數(shù)字巖心，定量計算孔喉結(jié)構(gòu)參數(shù)，并對比分析不同儲層微觀孔喉結(jié)構(gòu)的差異性。

4.1 巖樣選取與數(shù)字巖心構(gòu)建

基于沉積與成巖作用研究，結(jié)合巖心物性分析測試結(jié)果，選擇紅河油田長8儲層典型井巖樣，基本信息見表1所示。

表1 紅河油田長8儲層典型井巖樣基本信息Table1 Information of rock samples from Chang8 reservoir in Honghe Oilfield

對3 塊巖樣進行表面處理，選擇特征區(qū)域進行SEM-Maps 掃描成像，得到分辨率為10 nm、大小為1 mm×1 mm 的8 位灰度圖（圖13）。其中，46 號巖樣的粒間孔較發(fā)育，且存在大孔隙；而73 號巖樣主要發(fā)育晶間孔，HH1 號巖樣同時發(fā)育粒間孔和晶間孔，但均未見到大孔隙。

圖13 紅河油田長8儲層3塊典型井巖樣SEM-Maps掃描圖像Fig. 13 SEM-Maps scanning images of three typical rock samples from Chang8 reservoir in Honghe Oilfield

按照前文所述方法，對研究區(qū)3 塊典型井巖樣的SEM-Maps掃描圖像進行處理，構(gòu)建數(shù)字巖心，并提取孔隙和喉道。結(jié)果如圖14所示，白色部分為孔喉空間，黑色為礦物顆粒。

圖14 紅河油田長8儲層3塊典型井巖樣數(shù)字巖心孔喉圖像Fig. 14 Pore throat images of three typical rock samples from Chang8 reservoir in Honghe Oilfield

4.2 孔喉結(jié)構(gòu)定量表征

在紅河油田長8儲層典型井巖樣數(shù)字巖心圖像孔喉提取的基礎(chǔ)上，開展定量表征與分析，明確不同油井巖樣的孔隙、喉道發(fā)育特征及孔喉連通性特征，為綜合評價儲層質(zhì)量提供重要依據(jù)。

4.2.1 孔隙發(fā)育特征

基于SEM-Maps 掃描圖像對巖樣進行孔隙分割（圖15），并對孔隙進行定量表征，對比分析3類儲層的平均孔隙半徑、平均形狀因子及面孔率，結(jié)果（圖16）表明，46 號巖樣的平均孔隙半徑最大，為1.92 μm；面孔率最高，為0.068；形狀最規(guī)則，其平均形狀因子為6.56。HH1 號巖樣的平均孔隙半徑最小，為0.91 μm；面孔率最低，為0.038；孔隙不規(guī)則，其平均形狀因子為12.02。73 號巖樣的孔隙定量表征參數(shù)介于46 號和HH1 號之間。因此，46 號巖樣的儲集性能明顯優(yōu)于73號和HH1號巖樣。

圖15 紅河油田長8儲層3塊典型井巖樣孔隙分割圖像Fig. 15 Pore segmentation images of three typical rock samples from Chang8 reservoir in Honghe Oilfield

圖16 紅河油田長8儲層3塊典型井巖樣孔隙特征參數(shù)對比Fig. 16 Comparison of pore characteristic parameters of three typical rock samples from Chang8 reservoir in Honghe Oilfield

4.2.2 喉道發(fā)育特征

基于SEM-Maps 掃描圖像對研究區(qū)巖樣進行喉道分割（圖17），量化對比不同巖樣的平均喉道長度、平均喉道寬度及平均配位數(shù)，結(jié)果（圖18）表明，46 號巖樣的平均喉道長度最短，為14.92 μm；平均喉道寬度最大，為0.35 μm；平均配位數(shù)最大，為3.95。HH1 號巖樣的平均喉道長度最長，為25.05 μm；平均喉道寬度最小，為0.22 μm；平均配位數(shù)最小，為1.72。因此，研究區(qū)46 號巖樣的流體滲流能力遠優(yōu)于73號和HH1號巖樣。

圖18 紅河油田長8儲層3塊典型井巖樣的喉道特征參數(shù)對比Fig. 18 Comparison of throat characteristic parameters of three typical rock samples from Chang8 reservoir in Honghe Oilfield

4.2.3 孔喉連通性特征

基于SEM-Maps 掃描圖像對研究區(qū)3 塊典型井巖樣提取連通體，并用不同顏色標(biāo)記，結(jié)果（圖19）顯示，各連通體之間相互獨立。46 號巖樣最大連通體面積占比為33.04%，基本貫穿整個巖樣，表明流體能夠通過整個樣品；73 號巖樣最大連通體面積占比為27.69%，連通性相對較差；HH1 號巖樣最大連通體面積占比為8.36%，表現(xiàn)為局部連通，整體不連通。因此，46 號巖樣連通性明顯優(yōu)于73 號和HH1號巖樣。

圖19 紅河油田長8儲層3塊典型井巖樣的孔喉連通性分析Fig. 19 Pore throat connectivity analysis of three typical rock samples from Chang8 reservoir in Honghe Oilfield

5 結(jié)論

基于孔喉形態(tài)學(xué)差異，綜合應(yīng)用圖像腐蝕、膨脹和邏輯運算，提出一種基于圖像形態(tài)學(xué)特征的致密砂巖數(shù)字巖心圖像孔喉分割方法，可以準確地分割孔隙和喉道；建立直接基于數(shù)字巖心圖像的孔隙、喉道以及孔喉連通性定量表征方法，并提出各表征參數(shù)的具體計算方法，實現(xiàn)了致密砂巖微觀孔喉結(jié)構(gòu)定量表征，為分析致密砂巖儲層質(zhì)量提供了依據(jù)。對比分析紅河油田長8 儲層3 口典型井巖樣的孔喉結(jié)構(gòu)特征，結(jié)果表明：物性最好的46 號巖樣的平均孔隙半徑最大、面孔率最大，平均喉道長度最小、配位數(shù)最高，最大連通體面積占比最大；物性最差的HH1 號巖樣的平均孔隙半徑最小、面孔率最小，平均喉道長度最大、配位數(shù)最低，最大連通體面積占比最小。

符號解釋

A——孔喉的面積；

A1——區(qū)域U內(nèi)部的1個目標(biāo)區(qū)域；

A2——腐蝕運算后的殘缺孔隙；

A3——膨脹運算后的殘缺孔隙；

A4——準確的孔隙；

A5——準確的喉道；

Ab——第b個連通體像素個數(shù)，個；

Ap——孔隙的面積；

As——圖像中孔喉的總像素個數(shù)，個；

At——喉道的面積；

b——連通體的個數(shù)，個；

cIJ——孔喉連接關(guān)系矩陣；

CNm——配位數(shù)；

d——喉道的平均寬度；

f（i，j），f（k，l）——原圖像灰度值；

F——孔隙的形狀因子；

g（i，j）——降噪后圖像灰度值；

h（i，j）——灰度拉伸后圖像的灰度值；

i，k——圖像x方向坐標(biāo)；

I——孔喉連接關(guān)系矩陣的行數(shù)；

j，l——圖像y方向坐標(biāo)；

J——孔喉連接關(guān)系矩陣的列數(shù)；

L——喉道的長度；

m——孔隙的個數(shù)，個；

M——孔隙的像素個數(shù)，個；

n{ }Q——孔喉圖像中匹配的模板圖像Q的個數(shù)，個；

N——數(shù)字巖心的總像素個數(shù)，個；

P——孔喉的周長；

Pp——孔隙的周長；

Pt——喉道的周長；

Q1，Q2，Q3，Q4，Q5——序號為1，2，3，4，5的模板圖像；

S1——結(jié)構(gòu)元素1；

S2——結(jié)構(gòu)元素2；

T——與1個喉道連接的孔隙個數(shù)，個；

U——平面區(qū)域；

w——降噪卷積核；

wd——空間域核；

wr——灰度值域核；

w（i，j，k，l）——降噪卷積核w在i，j，k，l處的值；

wd（i，j，k，l）——空間域核wd在i，j，k，l處的值；

wr（i，j，k，l）——灰度值域核wr在i，j，k，l處的值；

x——笛卡爾坐標(biāo)x軸位置；

y——笛卡爾坐標(biāo)y軸位置；

α——面孔率；

β——最大連通體面積占比，%；

σd——圖像噪音空間域標(biāo)準差；

σr——圖像灰度值域標(biāo)準差；

Φ——處理后的剩余元素。