頁巖儲層微觀結構非均質性圖像法定量表征

荊 玥，劉登科

（西安交通大學，陜西西安 710000）

非均質性是油氣儲層內在固有屬性，頁巖儲層的微觀非均質性影響了儲層的各項物性、孔隙分布以及礦物分布等。我國陸相頁巖微觀結構差異較大，有機質與無機礦物的分布呈紋層、似層狀以及分散等分布結構，具有極強的非均質性[1]。

頁巖儲層的微觀非均質性表征依賴測量方法，國內外研究者利用鑄體薄片[2]、掃描電鏡[3]、微納米CT[4]等圖像技術定性描述不同類別的儲層微觀結構，或者利用低溫氮氣吸附[5]、二氧化碳吸附[6]、高壓壓汞[7]等流體注入實驗定量計算分析儲層的礦物、孔隙分布、含油氣性等差異。但是，基于圖像法定量表征頁巖儲層微觀結構非均質性的研究較少。本文基于圖像處理礦物識別技術定量研究頁巖儲層微觀結構的非均質性，以期豐富儲層微觀結構非均質性研究方法。

本文首先研究頁巖儲層微觀結構對巖石薄片圖像特征參數的影響，其次進行了特征參數的篩選與降維，最后定量分析頁巖樣品的微觀結構類型。

1 圖像處理及參數提取方法

在頁巖儲層的研究中，鑄體薄片鑒定、場發射環境掃描電鏡、微納米CT 掃描等圖像技術能直觀可視研究儲層性質。本文選取鄂爾多斯盆地三疊系的延長組長7 頁巖氣儲層的薄片圖像數據進行分析，16 組薄片圖像按照礦物微觀結構被人工劃分為分散狀、層狀以及塊狀分布等三種類型（圖1）。本節首先介紹薄片圖像的礦物識別分割處理方法，之后介紹礦物結構非均質性及紋理特征參數的提取方法。

1.1 圖像分割處理

頁巖薄片的彩色圖像中每個像素由0～255 范圍的RGB 三個分量組成，對這三個分量使用加權平均法計算可得到像素灰度值，所得灰度圖像見圖1d～f。通常灰度圖像選取合適的閾值可分割出感興趣的研究區域進行下一步分析，但是頁巖薄片圖像的消光特征使得不同角度下礦物的灰度差異過大，在研究礦物微觀結構時灰度閾值分割并不準確。

礦物識別分割可以提高圖像分析的精度，本文使用開源工具包ilastik1.3.3 對研究區的16 個薄片圖像進行礦物識別，ilastik 工具包使用隨機森林算法基于頁巖圖像中不同礦物的顏色、強度及紋理等特征在不同的尺度上計算每個像素屬于幾個類之一的概率。使用工具包的像素分類工作流可以獲得較理想的頁巖微觀結構圖像（圖1j～i），同時，之后進行特征參數提取工作時還可以將分割圖像進行二值化處理，突出微觀礦物在有機質中的分布特征。

圖1 頁巖薄片圖像

1.2 圖像特征參數的提取

礦物含量的非均質性特征采用代表視域面積與變異系數兩種方式描述。

1.2.1 代表視域面積 代表視域面積是一種基于圖像數據評價儲層非均質性的方法[8]，統計從小到大不同尺度下的視域中評價要素所占的比重，當比重穩定時的視域面積就是這種因素的代表視域，代表視域越大，表示該因素非均質性越強。

本文將礦物分布二值化圖從9 種路徑依次擴大視域，統計相對應路徑和邊長視域的礦物占視域面積的百分含量，當礦物含量穩定時提取參數。

1.2.2 變異系數 變異系數是數據標準差與平均數的比，可以比較數據離散程度大小，也是評價儲層非均質性的方法[3]，可以通過水平和垂直兩個方向上的變異系數分別描述兩個方向上的非均質性，通過式（1）計算水平方向變異系數，式（2）計算垂直方向變異系數：

將礦物在有機質中的分布結構看作圖像紋理，通過灰度共生矩陣（Gray Level Co-occurrence Matrix）計算得到薄片圖像的紋理特征參數。對于灰度級為L 的圖像，通過式（3）可生成L×L 的灰度共生矩陣：

式中：P（i，j，d，θ）-坐標為（x，y）與（x+a，y+b）的像素點，在θ 方向d 距離上出現灰度級組合為（i，j）的頻率；θ-灰度共生矩陣的生成方向，一般選取0°、45°、90°、135°進行計算。本節首先將薄片圖像分為八個灰度級從灰度共生矩陣提取出對比度、逆差矩、相關性及熵四個紋理特征。

1.2.2.1 對比度 對比度描述像素點與其領域像素的亮度對比情況，它反映了圖像中同尺度變化的強烈程度，在頁巖薄片圖像中，黃鐵礦等圖像亮度較高的礦物與圖象亮度低的有機質之間的遞變越多，灰度共生矩陣的對比度越大。

1.2.2.2 相關性 相關性描述了灰度共生矩陣元素在行或列方向上的相似程度，它反映了水平或垂直方向上圖像紋理的情況，即如果在某方向上礦物和有機物的紋理越豐富，則此方向的相關度值越大。

1.2.2.3 逆差矩 紋理逆差矩描述了圖像紋理的局部變化，它反映了圖像紋理的同質性，在圖像的不同區域中的礦物分層結構越像，紋理逆差矩越大。

1.2.2.4 熵 熵度量圖像所具有的紋理信息的隨機性，它表示了圖像中紋理的復雜程度，即圖像中礦物分布越分散越復雜，灰度共生矩陣的熵越大。

此外，對于二值圖像，灰度共生矩陣有特殊的三個參數可以描述礦物的分層結構。對于二值圖像，灰度只有兩級，所以每個方向上的共生矩陣均為一個2×2 的矩陣。

若將礦物賦值為1，則主對角上的元素分別表示連續的礦物和有機質的個數，副對角元素分別表示礦物與有機質之間的遞變。在同尺度下，礦物結構的周長是所有方向上的副對角線元素之和La，該值越大頁理發育越復雜。0°和90°方向下的副對角元素之比表征礦物結構的長徑比Lc，該值越小礦物結構越扁平。同一方向的副對角元素之差Ld 表征該方向下礦物和有機質的遞變程度，該值越大，遞變程度越強，且符號代表了遞變方向。

2 結果與討論

2.1 礦物含量的非均質性特征參數的分析

代表視域面積與變異系數兩種非均質性特征參數定量礦物含量非均質性特征，結果見圖2。由圖2 可以看出，代表視域面積參數可以較好的區分層狀分布和塊狀分布，分散狀分布的參數值變化范圍介于層狀與塊狀參數之間，可見，代表視域面積可部分定量礦物含量非均質性。同時，水平變異系數和垂直變異系數聯合應用也可以很好描述塊狀分布結構在兩個方向上的礦物含量非均質性，但是層狀分布與分散狀分布參數有重合的特征。所以本文選擇代表視域面積作為礦物含量非均質性特征參數。

圖2 三類分布礦物含量非均質性參數

2.2 紋理特征參數的篩選

上節從圖像中提取出7 個可以定量描述礦物分層結構的特征參數，雖然增加特征參數的個數有利于更加全面的描述頁巖微觀結構的改變帶來的圖像視覺形態變化，但是，參數維度的增加帶來更大的計算量且可能過分關注局部細節，從而造成更大的統計誤差，不利于后續的聚類分析，因此需要對上述特征參數進行篩選與降維。

通過對所提取的定量描述礦物分層結構的特征參數之間的相關關系進行研究，得到圖3 所示的相關性矩陣圖，篩選去除強正相關的參數有利于后續分析，將相關性系數超過0.8 的參數視為強正相關，對比度、逆差矩、周長均與其他參數強相關，因此剔除這三個參數。

圖3 圖像紋理特征相關性矩陣圖

2.3 特征參數的降維

通過篩選去除極相似的特征參數，得到一個五維向量可以用來表征礦物的分層結構，但是，五維向量在計算中仍然復雜，需要進一步降維。為此，引入主成分分析方法對多維特征向量進行簡化。主成分分析方法作為一種經典的降維方法，能夠將原來的多個指標重新組合成維度更低的一組互相無關的指標，來衡量圖像的特征。

對特征參數之間的協方差矩陣進行特征值分解，可以得到5 個特征值以及對應的特征向量。特征值在特征值總和中所占的比例就是該主成分的貢獻率。按照從大到小的順序排列，具體的貢獻率見圖4。前兩個主成分即可體現出原圖像特征向量超過98%的有效信息。采用主成分分析方法將五維向量簡化映射到二維向量，降低了多維特征向量的復雜度，可以更加直觀地描述礦物分層結構。兩個主成分分量通過式（8）計算，形成了一個描述礦物分層結構特征的二維向量。

圖4 各主成分貢獻率

通過計算礦物分層結構特征向量之間的歐式距離，使用最短距離法進行hierarchy 層次聚類，首先將數據中的每個樣本看作初始聚類簇，算法運行中找出距離最近的兩個聚類簇進行合并，該過程不斷重復直到預設的聚類簇個數。通過聚類分析得到三種類型的頁巖微觀礦物結構。二維向量具有較好的聚類效果，并且聚類結果可以很好的與儲層微觀結構的人工分類相對應。

3 結論

本文對鄂爾多斯長7 組頁巖氣儲層的微觀礦物結構非均質性進行了研究分析，通過提取巖樣薄片圖像的形狀及紋理特征，經過參數降維、定量聚類研究了三種類型的頁巖微觀礦物結構，具體結論如下：

（1）對于薄片鑒定圖像，鄂爾多斯盆地長7 儲層頁巖的微觀礦物結構有分散狀、層狀以及塊狀分布三種類型，每種類型在視覺上有較大差異。

（2）對于圖像法，代表視域面積作為礦物含量非均質性特征參數比水平和垂直兩方向上的變異系數更直接的表征了礦物含量的非均質性。

（3）通過相關系數篩選和主成分降維可以得到量化頁巖儲層微觀結構非均質性的二維向量，使用該向量定量聚類的結果可以有效的與薄片鑒定結果對應。