基于圖像處理和Kozeny-Carman方程的砂巖儲層滲透率預測*
2022-06-23 04:48:12潘衛(wèi)國
中國海上油氣 2022年2期
馮 進 石 磊 管 耀 潘衛(wèi)國 楊 清
(中海石油(中國)有限公司深圳分公司 廣東深圳 518054)
滲透率指在一定壓差下，巖石允許流體通過的能力，它是表征巖石本身傳導液體能力的參數(shù)，該參數(shù)對油氣儲層的評價至關重要。目前儲層滲透率評價的常用方法主要有兩大類：第一類是實驗室直接測量[1-3]，即利用滲透率儀測量柱塞巖心的滲透率。該方法的優(yōu)點是準確度高，但需要大量的柱塞巖心進行測量。第二類是測井計算獲取[4-5]，根據(jù)所采用的測井方法不同，可分為常規(guī)測井計算滲透率、核磁共振測井計算滲透率和偶極聲波測井計算滲透率方法。該方法的優(yōu)點是可以連續(xù)計算井筒方向的儲層滲透率，但也需要一定數(shù)量的巖心氣測孔隙度和滲透率數(shù)據(jù)用于方法建模。
當柱塞巖心不易獲取(巖心易碎等原因)，巖心氣測滲透率數(shù)據(jù)缺乏的情況下，儲層的滲透率評價會比較困難。已有研究表明，巖石滲透率的大小與孔隙度和孔隙網(wǎng)絡連通性等微觀特性密切相關。徐艷玲[6]、連會青[7]等對SEM圖像進行處理，獲取了孔隙面積、孔隙周長、孔隙半徑、分形維數(shù)等參數(shù)，并建立滲透率與孔隙面積(或孔隙周長)和分形維數(shù)的擬合關系?？紤]到圖像處理“盒計數(shù)法”獲取分形維數(shù)的局限性，姜濤 等[8]嘗試使用“小島法”對砂巖SEM圖像孔隙的面積和周長進行分析并獲取分形維數(shù)，并開展了滲透率預測。目前，已有的基于地質(zhì)圖像的滲透率預測方法取得了一定的效果，但也存在一些有待改進的地方，例如粘連孔隙準確分割問題、孔隙連通性與滲流通道彎曲度問題、滲透率計算公式的巖石物理意義等。考慮到以上問題，本文擬利用鑄體薄片、掃描電鏡等巖石微觀孔隙圖像(通過小碎塊巖心制備觀察樣品)，獲取巖石微觀參數(shù)，基于Kozeny-Carman方程來計算砂巖滲透率。
1 Kozeny-Carman方程介紹
20世紀20至50年代左右，Kozeny[9]和Carman[10-11]提出了如式(1)所示的半經(jīng)驗、半理論的多孔介質(zhì)滲透率模型：
(1)
式(1)中：k為巖石滲透率，mD；μ為流體黏度，Pa·s；γ為單位轉(zhuǎn)換因子；φ為巖石有效孔隙度，小數(shù)；S為巖石比表面積，1/cm；CK-C為地區(qū)經(jīng)驗系數(shù)。
在式(1)的基礎上，大量學者針對不同特征的多孔介質(zhì)以及假設條件，衍生推導出了一系列聯(lián)系滲透率與孔隙度、顆粒半徑、比表面、彎曲度等參數(shù)的公式，這類結(jié)構(gòu)相似且廣泛運用的公式就是著名的Kozeny-Carman方程。以下分別介紹幾類具有代表意義的Kozeny-Carman方程。
1.1 基于毛管束模型的Kozeny-Carman方程
根據(jù)達西定律，通過巖石的流量Q可表示為式(2)。Kozeny-Carman方程假設巖石可以抽象為一系列相互平行的彎曲毛管束，流體在毛管束中流動(毛管束模型)。假設毛管的截面為圓形，半徑為r，毛管中流體的流動可看作層流，單個毛管的流量Q可表示為式(3)。
(2)
式(2)中：Q為流量，m3/s；A為巖石截面積，m2；L為巖石長度，m；Δp為入口端和出口端的壓力差，Pa。
(3)
式(3)中：r為毛管半徑，m；l為彎曲毛管的實際長度，m。
假設毛管截面形狀為圓形，則巖石孔隙度φ和比表面積S可分別表示為式(4)和式(5)：
(4)
(5)
式(4)、(5)中：τ為彎曲度(此處針對毛管)，τ=l/L，無量綱。
聯(lián)立式(2)、(3)，可得：
(6)
將式(4)、(5)代入式(6)中，可得到基于毛管束模型的Kozeny-Carman方程：
(7)
1.2 基于球狀顆粒模型的Kozeny-Carman方程
公式(7)為根據(jù)毛管束模型推導的Kozeny-Carman方程，但實際巖石的孔隙結(jié)構(gòu)與毛管束模型還是有一定差距的，相對而言，使用球狀顆粒模型來表征巖石孔隙結(jié)構(gòu)會更準確一些。因此假設巖石中球狀顆粒的數(shù)量為n，顆粒直徑為D，則巖石中所有顆粒的體積和表面積分別為nπD3/6和nπD2，顆粒占巖石總體積的百分比為1-φ，巖石總體積可表示為nπD3/[6×(1-φ)]，巖石比表面積可表示為公式(8)：
(8)
式(8)中：D為顆粒直徑，m；n為巖石中顆粒數(shù)量。
將式(8)代入式(7)中，可進一步得到球狀顆粒模型的Kozeny-Carman方程：
(9)
1.3 其他Kozeny-Carman方程
對于黏土，房營光 等[12]采用質(zhì)量比表面的方式來展示Kozeny-Carman方程(式10)：
(10)
式(10)中：Cs為孔隙形狀參數(shù)；ρs為黏土密度，g/m3；Ss為黏土質(zhì)量比表面，m2/g。
Mavko和Nur[10]提出了利用孔隙直徑取代顆粒直徑的Kozeny-Carman方程(式(11))：
(11)
式(11)中：D為顆粒直徑，m。
Mavko等[13]引入孔隙連通性概念，在式(9)的基礎上提出了修改后的Kozeny-Carman方程(式(12))。該公式中，新增加了滲流孔隙度下限φp，即低于該孔隙度的巖石滲透率接近0。
(12)
式(12)中：φp為滲流孔隙度下限值，小數(shù)。
2 基于圖像處理的滲透率模型構(gòu)建
前述幾種Kozeny-Carman方程(式(1)、(7)、(9)～(12))中均使用了孔隙度φ，該參數(shù)主要表征滲流空間占巖石總體積的比例，該參數(shù)獲取相對較容易，并且是滲透率最主要的影響因素之一。另外滲流通道的寬度對巖石滲透率也有非常重要的影響，有些公式采用比表面積S對其進行表征(式(1)、(7)、(10))，有些公式采用顆粒直徑對其進行表征(式(9)、(12))，有些公式則采用孔隙直徑對其進行表征(式(11))，到底哪一參數(shù)最適合？這取決于參數(shù)獲取的難易程度以及Kozeny-Carman方程所應用巖石的特點。一般而言，獲取準確的比表面S相對比較困難，需要開展氣體吸附比表面積檢測(BET)，且該方法主要適用于微孔和介孔發(fā)育的巖石[14-18]，對于宏孔(半徑>0.5 μm)發(fā)育的巖石效果較差；顆粒直徑的獲取方法相對成熟且簡單，可通過篩析法、沉降法、激光法或鑄體薄片觀察法等粒度分析方法獲得[19-20]；孔隙直徑的獲取方法較多，常用的有壓汞法、鑄體薄片觀察法和微米CT掃描法[21-22]。另外，為了表征孔隙網(wǎng)絡的復雜性，大部分公式還使用了孔隙彎曲度τ(式(7)、(9)、(11)、(12))，一部分公式使用了孔隙形狀參數(shù)Cs(式(10))。
本文旨在提出一種僅依賴巖石顯微圖像的滲透率計算模型，因此在綜合考慮各種參數(shù)對滲透率的影響程度和各種參數(shù)獲取難易程度的基礎上，基于公式(9)提出如下形式的Kozeny-Carman方程用于巖石滲透率計算：
(13)
式(13)中：c為綜合了單位轉(zhuǎn)換和地區(qū)經(jīng)驗的系數(shù)；φa為顯微圖像面孔率，小數(shù)。
3 基于圖像處理的巖石微觀參數(shù)確定
在缺乏巖石柱塞樣品的情況下，氣體滲透率測試無法進行，但只需要極少量巖石碎塊或巖屑即可進行巖石微觀孔隙結(jié)構(gòu)觀察并拍攝照片。常用的巖石微觀孔隙結(jié)構(gòu)觀察手段包括鑄體薄片、掃描電鏡、微(納)米CT等。本文以鑄體薄片圖像為例，闡述巖石微觀圖像處理及式(13)中滲透率計算所需參數(shù)的獲取方法。
3.1 面孔率確定
鑄體薄片將液態(tài)染色樹脂在真空加壓條件下注入巖石孔隙中，待染色樹脂固化后磨成薄片，可在透射偏光顯微鏡下清楚的觀察巖石孔隙結(jié)構(gòu)。利用各種圖像處理方法，可以基于鑄體薄片照片獲取微觀尺度下的巖石骨架和孔隙喉道信息。圖1展示了巖石鑄體薄片照片及二值化后孔隙網(wǎng)絡圖像，圖1a中藍色樹脂填充部分代表孔隙和喉道，圖1b為經(jīng)過圖像處理后二值化的孔隙網(wǎng)絡圖像，黑色部分代表孔隙和喉道。圖1b本質(zhì)上為一個二維數(shù)組(x,y)，x和y分別為圖1b橫向和縱向上的像素點個數(shù)，統(tǒng)計圖1b中黑色像素點的個數(shù)n，則圖1b的面孔率可表示為式(14)：
(14)
式(14)中：x為圖像橫向上的像素點個數(shù)；y為圖像縱向上的像素點個數(shù)；n為圖像中孔隙所占的像素點個數(shù)。
本文的圖像處理過程中，會生成兩種圖像，第一種是直接通過閾值分割并二值化的孔隙網(wǎng)絡圖像(圖1b)，第二種是在該圖像基礎上對孔喉進行加強(使用粘連顆粒分割算法)后的孔隙度網(wǎng)絡圖像(圖1c)。對于第一種孔隙網(wǎng)絡圖像，由于鑄體薄片圖像的分辨率有限，直接采用閾值分割，會導致低于圖像分辨率的微孔隙和喉道不計入面孔率中，因此獲得的面孔率偏低。對于第二種孔隙網(wǎng)絡圖像，由于孔喉得到了增強，獲得的面孔率大幅度提高，但也存在一個問題，低孔滲巖石的喉道部分會得到額外增強，使得相對誤差增大。在對比分析這兩種圖像結(jié)果時，發(fā)現(xiàn)如下現(xiàn)象：①第一種面孔率與氣測孔隙度呈較好的線性關系，明顯優(yōu)于第二種面孔率；②第一種面孔率基本反映了對滲透率有主要貢獻的大孔隙；③第一種面孔率所經(jīng)歷的操作較少，確定性高，有利于后期將算法轉(zhuǎn)化為自動化處理軟件。因此綜合考慮以后，論文選擇使用第一種孔隙網(wǎng)絡圖像的面孔率。
圖1 巖石鑄體薄片照片及二值化后孔隙網(wǎng)絡圖像
3.2 顆粒直徑確定
巖石顆粒經(jīng)孔隙和喉道分割以后(圖1c)，可分別對巖石顆粒進行編號(圖2a)，并提取其顆粒直徑。由于真實的巖石顆粒并非圓形，常用的顆粒直徑主要有3種：長軸直徑、短軸直徑和平均直徑。長軸直徑指過顆粒重心(圖2a中紅色數(shù)字編號所在位置即為該顆粒的重心)的最長直徑值，短軸直徑指過顆粒重心的最短直徑值，平均直徑則是指過顆粒重心的多條直徑平均值(圖2a)。將圖像識別得到的顆粒平均直徑按照從大到小排序，選取前10%顆粒，計算其平均直徑的平均值作為巖石顆粒直徑值(圖2b)。
3.3 孔隙彎曲度確定
1937年，Carman[23]最早提出孔隙彎曲度(水力彎曲度)的概念，用來表征多孔介質(zhì)中流動路徑的彎曲特征，此后孔隙彎曲度在多孔介質(zhì)滲流及導電特性等領域被廣泛應用?？紫稄澢鹊墓娇珊唵伪硎緸椋?/p>
τ=l/L
(15)
圖2 巖石顆粒平均直徑提取與代表性顆粒直徑確定示意圖
式(15)中：τ為彎曲度(此處針對孔隙)，無量綱；l為彎曲的孔隙路徑長度，m；L為巖石樣品長度，m。
對于結(jié)構(gòu)規(guī)則的多孔介質(zhì)，孔隙彎曲度可以很容易通過孔隙幾何形態(tài)計算得到。但天然巖石孔隙結(jié)構(gòu)復雜，其孔隙彎曲度的獲取就要困難得多。本文使用的方法是通過孔隙網(wǎng)絡中軸提取并計算最短路徑來獲得孔隙彎曲度。
3.3.1孔隙網(wǎng)絡中軸提取
1967年，Blum[24]提出著名的“火燒草場模型”，形象地闡述了中軸的概念：想象有一塊由二維閉合曲線圍成的草場，在這條二維閉合曲線上(草場邊界)放置火種，同時點燃，火焰以相同的速度向圖形內(nèi)部的各個方向燃燒，火焰相遇時熄滅，火焰熄滅點的集合即為該草場的中軸(也稱為骨架)。
常用的中軸提取方法包括Voronoi方法、演化方法、距離場函數(shù)方法、拓撲細化方法等幾大類。本文采用的是Fouard等[25]在進行血管圖像處理時提出的一種倒角距離變換方法。該方法從圖像邊緣向中心刪除血管的像素點，直到無法繼續(xù)刪除像素點為止，最終形成血管的中軸。巖石孔隙網(wǎng)絡與血管網(wǎng)絡具有一定相似性，因此可以使用該倒角距離變換方法提取巖石孔隙網(wǎng)絡的中軸。
圖3為基于倒角距離變換方法的巖石孔隙中軸提取示意圖。圖3a為倒角距離變換所用的掩模。圖3b為鑄體薄片圖像處理后的孔隙-顆粒二值圖某局部區(qū)域，圖中一個小方塊代表一個像素點，白色小方塊表示巖石骨架，黑色小方塊表示孔隙，Ii和Ij分別表示x方向和y方向的圖像尺寸。通過以下2個步驟的掃描來計算巖石孔隙網(wǎng)絡的距離圖：①從圖像左上角到右下角，即坐標(1,1)到坐標(Ii,Ij)，利用正向掩模進行正向掃描；②從圖像右下角到左上角，即坐標(Ii,Ij)到坐標(1,1)，利用反向掩模進行反向掃描。在正向掩?；蚍聪蜓谀呙柰戤吅?，均需刪除圖中低于某一閾值的像素點，并且重新進行二值化，分別得到正向掩模或反向掩模處理后的圖像。再將這兩張圖像取交集即可得到最終的孔隙網(wǎng)絡中軸(圖3c)。
圖3 基于倒角距離變換的巖石孔隙網(wǎng)絡中軸提取示意圖
3.3.2最短路徑確定
孔隙網(wǎng)絡中軸提取以后，下一步是確定流體流動的最短路徑(最短優(yōu)勢路徑)。最短路徑的定義如下：在網(wǎng)絡圖中任意兩點之間的連接線路存在兩條及以上的情況下，肯定存在一條總權值最小的線路，即這兩點之間的最短路徑。本文求取孔隙網(wǎng)絡最短路徑使用的是Dijkstra算法[26]，具有計算速度快、穩(wěn)定性好、準確率較高等優(yōu)勢，目前廣泛應用于交通、導航、電網(wǎng)、通信等領域。
Dijkstra算法使用類似廣度優(yōu)先搜索的方法解決賦權圖的單源最短路徑問題。對于已提取中軸的巖石鑄體薄片孔隙網(wǎng)絡圖，可以用黑色線條表示孔隙網(wǎng)絡中軸(圖4a)，也可以用彩色線條表示中軸，不同顏色表征孔隙或喉道的半徑相對值(圖4b)。下面在圖4b中選取某一局部區(qū)域(圖4c)演示如何確定最短路徑：將所有線條的交點抽象為節(jié)點集合(圓形，代表孔隙)，兩節(jié)點間使用直線段連接(線條，代表喉道)，該直線段的距離用兩節(jié)點之間孔隙網(wǎng)絡中軸的像素點個數(shù)表示，節(jié)點集合(V)與帶權路徑集合(E)就構(gòu)成一個帶權有向連通圖G=
Dijkstra算法最終可得到孔隙網(wǎng)絡中軸圖起點和終點之間考慮權重(孔隙喉道大小)的最短路徑長度，即流體在孔隙網(wǎng)絡中的最短優(yōu)勢通道的像素點個數(shù)l。對于固定視域的鑄體薄片，其左右邊界距離L已知。利用公式(15)即可求得孔隙彎曲度。
圖4 孔隙網(wǎng)絡中軸的最短路徑確定示意圖
4 應用實例
以珠江口盆地番禺氣田珠江組作為研究對象，其儲層包括巖屑長石砂巖、長石巖屑砂巖和長石石英砂巖等，孔隙度分布在8%～24%，滲透率分布在0.1～5 000 mD。選取珠江口盆地番禺氣田珠江組12個不同深度位置的鑄體薄片樣品。
由于薄片圖像表征了極小尺度的巖石微觀結(jié)構(gòu)，而柱塞氣測滲透率則尺度相對較大，因此需要選擇合適尺度的薄片圖像來進行圖像處理。本文在選擇薄片圖像的視域時，主要考慮了如下2個原則：①以選定視域大小拍攝照片時，不同位置拍攝得到照片的面孔率和顆粒直徑差異較小，即該視域下的照片盡量不受巖石非均質(zhì)性的影響；②以選定視域大小拍攝照片時，要求有足夠的分辨率，孔隙、喉道、顆粒等信息要盡量清晰。由于原則①要求視域盡量大，而原則②要求視域盡量小，兩者相互制約，因此需要盡量平衡原則①和原則②，以獲取最佳效果。
本文對每個樣品均拍攝了多張不同放大倍數(shù)的照片，最終選擇放大倍數(shù)為50倍、視域大小為3.072 mm×3.072 mm的照片(該視域下兼顧了原則①和原則②)，用于圖像處理和巖石微觀參數(shù)的確定。12個砂巖樣品的最短路徑提取結(jié)果如圖5所示。12個砂巖樣品的鑄體薄片照片經(jīng)圖像處理后，提取了面孔率、顆粒直徑和孔隙彎曲度，具體參數(shù)見表1。
本文處理得到的面孔率雖然低于氣測孔隙度，但與氣測孔隙度與較好的線性關系，滿足滲透率計算要求；巖心粒度分析表明樣品最大顆粒直徑主要分布在0.25～0.71 mm，與本文處理得到的前10%最大顆粒直徑平均值(230.22 ～588.50 μm)相吻合。滲透率的計算分兩步：首先，選擇4個砂巖樣品(1-15號樣品，K=207.900 mD；1-13號樣品，K=83.660 mD；3-7號樣品，K=0.994 mD；3-18號樣品，K=0.375 mD)，以氣測滲透率作為標定，利用公式(13)得到這4個樣品的地區(qū)經(jīng)驗系數(shù)c分別為3.53、2.58、4.22和2.57，取4個樣品平均值3.23作為地區(qū)經(jīng)驗系數(shù)c的最終取值(當研究區(qū)暫時無氣測孔隙度數(shù)據(jù)時，一般c可在1～5取值)，因此在珠江口盆地番禺氣田珠江組，基于圖像處理的Kozeny-Carman方程公式(13)可進一步表示為公式(16)；然后，利用剩下8個樣品進行模型的檢驗(表1)。
圖5 珠江口盆地番禺氣田珠江組砂巖最短路徑提取結(jié)果Fig .5 Shortest path extraction results of sandstone of Zhujiang Formation in Panyu gas field, Pearl River Mouth Basin
表1 珠江口盆地番禺氣田珠江組砂巖微觀參數(shù)與滲透率數(shù)據(jù)Table 1 Microscopic parameters and permeability data of sandstone of Zhujiang Formation in Panyu gas field, Pearl River Mouth Basin
(16)
計算結(jié)果表明，利用圖像處理獲取的微觀參數(shù)(面孔率、顆粒直徑、孔隙彎曲度)與根據(jù)式(16)計算得到的珠江口盆地番禺氣田珠江組砂巖滲透率與氣測滲透率吻合性較好(表1、圖6)，表明本文提出的滲透率計算方法效果較好。
圖6 珠江口盆地番禺氣田珠江組砂巖氣測滲透率與 計算滲透率交會圖
Pearl River Mouth Basin
5 結(jié)論
論文提出了一種基于巖石顯微圖像(鑄體薄片照片、掃描電鏡圖片等)獲取巖石微觀結(jié)構(gòu)參數(shù)，并利用Kozeny-Carman方程計算巖石滲透率的方法。該方法適用于標準柱塞巖心缺乏的條件下，可以提供一定準確度的滲透率數(shù)據(jù)，為儲層滲透率的評價提供了新思路和新方法。該方法亦可集成于各類顯微鏡系統(tǒng)配套的圖像處理軟件中，提供多孔介質(zhì)樣品圖像拍攝與滲透率計算的智能一體化解決方案。
