一種頁巖氣儲層有機孔提取及測井計算新模型

歐發(fā)輝，賴富強，夏小雪，李嬌，馬麒，李慧

（重慶科技學院復雜油田勘探開發(fā)重慶市重點實驗室，重慶 401331）

0 引言

頁巖氣資源豐富， 在能源發(fā)展戰(zhàn)略中占有重要地位。頁巖氣通常以游離態(tài)廣泛存在于納米孔隙中，四川盆地川南地區(qū)五峰組—龍馬溪組的海相頁巖儲層發(fā)育大量有機質(zhì)孔和少量無機質(zhì)孔， 是頁巖氣的重要存儲空間［1-2］。

測井作為測量地球物理參數(shù)的手段， 在頁巖儲層基本特征和解釋評價方面發(fā)揮了重要作用， 可對礦物含量、 地化參數(shù)及物性參數(shù)等儲層關鍵參數(shù)進行精細評價［1-3］。 近年來，國內(nèi)外學者針對儲層測井評價模型進行的研究取得一定進展。Luffel 等［4］通過研究巖心總孔隙度與測井參數(shù)關系， 采用概率統(tǒng)計方法建立總孔隙度測井模型，此方法考慮因素單一，結果偏差較大。侯頡等［5］通過建立礦物體積模型及測井響應方程，采用最優(yōu)化算法反演計算出礦物組分含量和總孔隙度，然而此模型中有許多參數(shù)需要確定， 難以準確評估且操作復雜。 Herron 等［6］利用礦物元素質(zhì)量百分比計算礦物骨架參數(shù)值，隨后基于礦物體積模型計算孔隙度，最終建立優(yōu)化的總孔隙度測井模型。 由于建立的測井模型有機質(zhì)因素未考慮在內(nèi)， 導致計算結果精度比較低，普適性不高，不適用于有機質(zhì)孔發(fā)育的地區(qū)。因此，有必要建立一種新的頁巖氣儲層有機質(zhì)孔隙度的評價方法。

圖像處理技術是當前各個行業(yè)的重點研究手段，傳統(tǒng)的閾值分割算法和邊緣檢測算法已經(jīng)在多種行業(yè)中得到應用。 如李應果等［7］通過閾值分割算法有效識別了單板穿孔缺陷現(xiàn)象，齊瑞燕等［8］利用迭代閾值分割算法處理電成像圖像，效果良好。

為更精細地對不同層段有機質(zhì)孔隙度進行計算，本文以川南地區(qū)五峰組—龍馬溪組頁巖（來自取心井X1，X2）為例，對研究區(qū)不同深度的掃描電鏡（SEM）圖像進行有機質(zhì)孔隙提取及面孔率計算， 進而代入面孔率和孔隙度關系式，對有機質(zhì)孔隙度進行計算，最終建立測井計算模型；同時進行核磁實驗，利用實驗結果計算的有機質(zhì)孔隙度與基于SEM 圖像有機質(zhì)孔隙度的計算結果進行對比，驗證該方法的準確性。 結果表明，該項工作有利于頁巖氣儲層有機質(zhì)孔縱橫向展布規(guī)律研究，可完善頁巖氣儲層甜點的精細測井評價。

1 區(qū)域地質(zhì)背景

研究區(qū)主體位于滇黔北坳陷的中西部區(qū)域（川南低陡褶皺帶南緣帶）［6-10］（見圖1）， 該區(qū)五峰組—龍馬溪組地層具有厚度大、總有機碳豐度高的特點，且保存良好，頁巖氣賦存量巨大［9-10］。

圖1 研究區(qū)構造及實驗取心井位置Fig.1 Structure and location of experimental coring wells in the study area

本次研究的32 個樣品均采自研究區(qū)典型取心井（X1，X2），取心井位置見圖1。 為反映儲層物性縱向變化，達到刻度測井解釋目的，基本上等間距取樣，取樣縱深跨度基本為30～35 m，每口井取樣8 個，基本覆蓋了所有層位（五峰組—龍一14層），覆蓋了所有優(yōu)質(zhì)儲層，在底部最優(yōu)質(zhì)儲層段，加密了取樣間隔。

五峰組—龍馬溪組， 巖性整體由硅質(zhì)頁巖逐漸過渡為灰質(zhì)頁巖到灰?guī)r（見圖2），表明五峰組—龍馬溪組沉積期，水體逐漸變淺，鹽度逐漸變高，氧化程度逐漸增強，總有機碳質(zhì)量分數(shù)逐漸減?。?1-17］。 研究區(qū)五峰組—龍馬溪組一段局部地區(qū)發(fā)育深水砂泥質(zhì)陸棚濁積巖砂體。龍馬溪組二段沉積期，水體隨深度減少而變咸且富含氧氣，水動力條件也隨之增強，總體上不利于有機質(zhì)的形成與保存。 滇黔北地區(qū)中西部以淺水泥質(zhì)陸棚沉積為主，五峰組—龍馬溪組，巖性顏色逐漸變淺，灰質(zhì)含量增加， 巖性整體上由硅質(zhì)頁巖逐漸過渡為灰質(zhì)頁巖到灰?guī)r。龍一段主要發(fā)育深水泥質(zhì)陸棚、灰泥質(zhì)陸棚及深水硅泥質(zhì)陸棚微相。 龍馬溪組二段灰質(zhì)含量增加，筆石化石含量相對減少，黃鐵礦含量減少，主要發(fā)育為淺水泥質(zhì)陸棚及淺水灰質(zhì)陸棚微相。因此，微相控制了頁巖的脆性礦物含量， 由龍一11亞段—龍一13亞段硅質(zhì)頁巖向上逐漸過渡為泥質(zhì)頁巖、灰質(zhì)頁巖。優(yōu)質(zhì)頁巖氣儲層主要發(fā)育在五峰組—龍一11亞段， 該層段水平層理發(fā)育［11］。

2 有機質(zhì)面孔率計算及核磁實驗標定

有機質(zhì)孔隙是四川盆地川南地區(qū)五峰組—龍馬溪組海相頁巖儲層非常重要的儲集空間［12］，所以對該地區(qū)的有機質(zhì)孔隙進行識別尤為重要。 在SEM 圖像中，有機質(zhì)孔與無機質(zhì)孔在亮度上有差異， 導致其成像后的閾值也呈現(xiàn)出明顯差異?；谶@種情況，本研究以川南地區(qū)X1，X2 井不同井段的掃描電鏡圖像閾值差異為出發(fā)點，利用Matlab 軟件，構建改進Canny-Otsu 算法，對有機質(zhì)孔進行提取并計算出其面孔率，同時搭建 核磁實驗，計算有機質(zhì)孔隙占比。

2.1 有機質(zhì)孔隙的識別及相關參數(shù)計算

2.1.1 閾值分割法原理

閾值分割主要采用Otsu 閾值分割法［13-14］，基本思路為：根據(jù)巖石掃描電鏡的實驗結果，利用Matlab 分析軟件編制相應的命令流， 通過計算得到掃描電鏡圖像不同區(qū)域的類間方差σ2。當圖像分割閾值在灰度范圍內(nèi)順序取值時，計算SEM 圖像的σ2；當計算的σ2值最大時，此時的分割閾值即為最佳分割閾值。 σ2計算方法見式（1）。

式中：ω0，ω1分別為不同灰度值的概率之和；uT為掃描電鏡圖像灰度均值；u0，u1分別為掃描電鏡圖像不同區(qū)域灰度均值。

利用Otsu 閾值分割算法對有機質(zhì)孔隙進行識別及參數(shù)表征的步驟為：1）對掃描電鏡圖像進行獲?。ㄒ妶D3a）；2）對獲取的掃描電鏡圖像，利用Otsu 閾值分割法對有機質(zhì)孔與有機質(zhì)進行閾值提取， 并利用灰度值繪制灰度直方圖［15-16］（見圖3b）；3）利用掃描電鏡圖像，通過特定的閾值，對有機質(zhì)進行提取（見圖3c）；4）利用掃描電鏡圖像， 通過特定的閾值對有機質(zhì)孔隙進行提?。ㄒ妶D3d）；5）通過有機質(zhì)孔隙的像素點和有機質(zhì)的像素點數(shù)量的比值進行面孔率求取。

圖3 有機質(zhì)孔隙提取及參數(shù)獲取Fig.3 Organic pore extraction and parameter acquisition

統(tǒng)計樣品的圖像數(shù)據(jù)， 可以進一步獲得樣品所處頁巖段的微觀孔隙結構特征， 但在有機質(zhì)孔隙提取和面孔率計算的過程中，會有誤差產(chǎn)生。產(chǎn)生誤差的原因主要有： 掃描電鏡圖像中有機質(zhì)孔隙與有機質(zhì)、 無機質(zhì)的閾值區(qū)分不明顯， 使得閾值分割無法達到理想效果； 在有機質(zhì)孔的提取中， 圖像中無機質(zhì)孔與有機質(zhì)孔的閾值過于近似， 使用閾值分割算法無法剔除裂縫和無機質(zhì)孔的影響， 對有機質(zhì)面積和有機質(zhì)孔面積計算造成誤差； 在計算面孔率時， 由于有機質(zhì)孔邊緣存在高亮區(qū)域，導致有機質(zhì)提取時存在空白區(qū)域，進而計算有機質(zhì)的面積偏小。

當前軟件還不能解決以上問題， 為了盡可能精確表征孔隙微細結構，需要將灰度變換增強［17］、邊緣檢測算法、閾值分割算法進行結合改進，才能使算法更適用于有機質(zhì)孔面孔率的計算，有效解決存在的問題。

2.1.2 Canny 邊緣檢測算法

Canny 邊緣檢測算法是當前較為出色的多邊緣檢測算法，在多個領域已經(jīng)得到廣泛應用，基本思想是通過圖像的幅度值與方向的極大值提取圖片的邊緣信息。 主要步驟為：1）將圖像灰度化后進行高斯濾波，以此達到去除高頻噪聲的目的；2）利用算子計算圖像梯度和強度，表征邊緣強度與方向；3）根據(jù)梯度方向?qū)D像梯度強度進行非極大值抑制，得到單像素邊緣點；4）設置雙閾值保留強邊界和潛在邊界，并連接邊緣。

引入邊緣檢測算法的目的是為了有效避免無機質(zhì)中的孔隙對于有機質(zhì)孔隙提取工作的影響。 但是傳統(tǒng)的Canny 算法對SEM 圖像進行處理會存在有機質(zhì)提取不準確、邊緣不閉合的問題，需要根據(jù)實驗需求進行拓展與改進。

2.1.3 有機質(zhì)孔隙識別及參數(shù)表征

因SEM 掃描電鏡圖像中地質(zhì)信息過多，存在孔隙邊緣高亮問題和無機質(zhì)孔隙的影響， 現(xiàn)對傳統(tǒng)Canny算法進行相應的改進，并結合Otsu 閾值分割算法和灰度變換增強，形成一種改進Canny-Otsu 算法。 利用改進Canny-Otsu 算法，對SEM 掃描電鏡圖像（X2 井龍一12亞段樣品（深度為2 248.88 m）有機質(zhì)孔掃描電鏡圖）進行處理（見圖4），具體步驟為7 步。

1）利用Otsu 閾值分割算法對SEM 圖像進行處理，繪制灰度直方圖， 獲取圖像有機質(zhì)與有機質(zhì)孔隙的閾值。 對原圖觀察分析看出，掃描電鏡圖像有機質(zhì)、有機質(zhì)孔隙區(qū)分不明顯。對直方圖進行分析得出，閾值基本為70～150，閾值集中，難以區(qū)分有機質(zhì)孔隙與有機質(zhì)。

2）對SEM 圖像的閾值進行灰度變換增強處理。已經(jīng)灰度化的SEM 圖像按一定變換關系逐點改變源圖像中每一個像素灰度值，使圖像中的有機質(zhì)、無機質(zhì)和孔隙的閾值特征更加明顯， 這有利于后期的有機質(zhì)和有機質(zhì)孔隙的提取?；叶茸儞Q增強處理得到的圖像（見圖4c）與原圖（見圖4a）對比發(fā)現(xiàn)，有機質(zhì)與有機質(zhì)孔隙變得更為突出。 由增強處理的灰度直方圖（見圖4d）可以看出，增強處理后的閾值更加趨近于兩端，處理后增強的灰度直方圖顯示圖像量化恰當， 使有機質(zhì)與非有機質(zhì)的閾值差異更大，在后續(xù)工作中，根據(jù)所求閾值范圍自適應提取有機質(zhì)面積更加精準。

3）采用自適應中值濾波代替高斯濾波平滑圖像，去除噪點。中值濾波是一種非線性平滑技術，采用濾波窗口對待處理的圖像矩陣進行處理， 將窗口中的圖像像素點進行排序，并采用中值替代窗口中的像素值，以達到平滑圖像的效果［18］。 中值濾波對比高斯濾波可以在兼顧保留像素信息的基礎上不對原圖像產(chǎn)生較大影響，但是中值濾波受濾波窗口大小的影響較大，易造成圖像不連續(xù)。 自適應中值濾波是在中值濾波的基礎上引入了自適應濾波算法的思想， 可以根據(jù)不同情況設好不同條件，動態(tài)改變中值濾波器的窗口尺寸。

4）對圖像進行非極大值抑制。

5）利用Otsu 算法確定雙閾值，進行閾值篩選。 利用第1）步提取到的有機質(zhì)與有機質(zhì)孔隙的閾值范圍［0，L］確定有機質(zhì)與無機質(zhì)閾值分割值t：有機質(zhì)與有機質(zhì)孔隙的閾值分割值為t0，［t+1，L］為無機質(zhì)閾值范圍；［0，t］為有機質(zhì)閾值范圍［19］。 令有機質(zhì)像素點的全幅比值為a0，平均灰度級為u0；無機質(zhì)像素點比值為a1，平均灰度級為u1，圖像總平均灰度為u，類間方差為σ。

根據(jù)Otsu 算法，當σ 取最大值時，此時的灰度級為最佳閾值，即為高閾值Th。 低閾值T1為

6）邊緣膨脹使有機質(zhì)的邊緣連接成封閉區(qū)域。SEM圖像中少數(shù)孔隙的閾值范圍接近于有機質(zhì)閾值范圍，影響有機質(zhì)提取（見圖5）。 所以，需要進行邊緣膨脹操作，經(jīng)過膨脹后的有機質(zhì)邊緣會形成閉合區(qū)域［20］。

圖5 無機質(zhì)孔放大圖像Fig.5 Enlarged image of inorganic pores

7）連通域面積過濾。 經(jīng)過膨脹后的有機質(zhì)邊緣會形成閉合，將閉合區(qū)域進行面積篩選，保留大面積，去除小面積， 可以有效規(guī)避無機質(zhì)孔隙與無機質(zhì)閉合區(qū)域的影響。 改進Canny-Otsu 算法的流程見圖6。

圖6 改進Canny-Otsu 算法計算有機質(zhì)提取流程Fig.6 Flow chart of organic matter extraction calculated by improved Canny-Otsu algorithm

利用改進Canny-Otsu 算法對實例圖進行處理，可以獲得有機質(zhì)的區(qū)域， 通過計算區(qū)域面積可以獲得有機質(zhì)的面積， 然后通過對有機質(zhì)區(qū)域進行二值化處理提取出有機質(zhì)孔隙并計算出有機質(zhì)孔隙面積。 傳統(tǒng)Otsu 閾值分割算法處理后結果與改進Canny-Otsu 算法的處理后結果對比見圖7、圖8。

圖7 有機質(zhì)區(qū)域提取區(qū)域Fig.7 Extraction results of organic matter regions

圖8 有機質(zhì)孔區(qū)域提取區(qū)域Fig.8 Extraction results of organic pore regions

由圖7a，7b 看出，有機質(zhì)中去除了無機孔的影響，并且有機質(zhì)孔邊緣高亮區(qū)域的影響也得到了去除。 由圖8a，8b 可以清晰看出，圖7 中的無機質(zhì)孔已經(jīng)去除。

2.1.4 表征結果

在X1 井和X2 井的電鏡圖像中，選取代表性的8張（具有統(tǒng)一比例尺）圖像進行有機質(zhì)和有機質(zhì)孔提取，計算有機質(zhì)孔面積與有機質(zhì)面積的比值，獲取面孔率。 最終計算結果見表1。

表1 X1 井和X2 井面孔率計算Table 1 Surface porosity calculation of Well X1 and X2

2.2 基于核磁實驗有機質(zhì)孔和無機質(zhì)孔比例計算

飽和水、 飽和油核磁對比實驗的一個目標是判斷有機質(zhì)孔、無機質(zhì)孔的比例［21］。 如圖9 所示，藍色數(shù)據(jù)點為壓力在20.68 MPa、 溫度在90 ℃狀態(tài)下飽和鹽水的核磁T2譜減去來樣狀態(tài)的核磁T2譜， 目的是消除黏土束縛水的影響［20-28］。 紅色數(shù)據(jù)點為壓力在20.68 MPa、溫度在90 ℃狀態(tài)下飽和油的核磁T2譜減去來樣狀態(tài)的核磁T2譜。

對比2 套T2譜的包絡面積可以分別得到有機質(zhì)孔（親油孔）和無機質(zhì)孔（親水孔）的孔隙度。X1 井2 個樣品有機質(zhì)孔與無機質(zhì)孔的比例分別為56.7∶43.3 和54.2∶45.8。 X2 井2 個樣品有機質(zhì)孔與無機質(zhì)孔的比例分別為54.6∶45.4 和51.9∶48.1。

綜合分析得出，研究區(qū)有機質(zhì)孔、無機質(zhì)孔的比例約為55∶45。

3 測井計算模型構建

利用研究區(qū)礦物元素類型及質(zhì)量分數(shù)構建的儲層礦物組分體積模型見圖10。 利用最優(yōu)化反演計算儲層的有機質(zhì)骨架干酪根、 無機質(zhì)骨架礦物質(zhì)量分數(shù)及孔隙度，為構建有機質(zhì)孔孔隙度計算模型做準備。

圖10 頁巖氣儲層礦物組分體積模型Fig.10 Volume model of mineral components in shale gas reservoirs

3.1 頁巖儲層礦物組分的確定

X1 井和X2 井均進行了QEMSCAN 礦物成分分析，與LithoScanner 巖性掃描測井解釋結果（巖心刻度后）對比結果見圖11。 由圖11 可以看出，整體上巖性掃描測井的礦物質(zhì)量分數(shù)與巖心測試結果一致性較強，為計算模型的建立提供了基礎。

圖11 X2 井QEMSCAN 礦物成分結果與標定后巖性掃描測井處理結果對比Fig.11 Comparison of QEMSCAN mineral composition results with lithological scanning logging results after calibration of Well X2

3.2 總有機質(zhì)骨架干酪根體積分數(shù)計算方法

1）總有機碳質(zhì)量分數(shù)計算方法。 將密度測井響應特征納入總有機碳質(zhì)量分數(shù)計算影響因素，可表示為

式中：TOC為總有機碳質(zhì)量分數(shù)，g/g；R為深側向電阻率，Ω·m；駐t為聲波時差，μs/m；ρ 為測井密度，g/cm3；A，B，C為無量綱擬合系數(shù)。

2）總有機質(zhì)骨架干酪根體積計算模型［22］。 依據(jù)密度測井，假設ρ 為主要測量無機礦物骨架密度值，TOC也可為

式中：ρtoc為有機碳密度，g/cm3；Vtoc為有機質(zhì)骨架干酪根體積分數(shù)。

根據(jù)式（5），可得：

井眼環(huán)境及含氣性的一些影響導致測井值偏低，或黃鐵礦的影響導致測井值偏高，因此，需要對其進行適當?shù)男Ｕ?，則校正后干酪根體積分數(shù)為［23］

式中：Vkero為校正后干酪根體積分數(shù)；λ 為校正因子。

式（7）中，對于Vtoc偏低情形時，λ 一般為1.0～2.0；對于Vtoc偏高情況下，λ 為0.5～1.0。 根據(jù)前人的研究成果，ρtoc一般為1.3 g/cm3，范圍為［1.2，1.5］。

3.3 有機孔孔隙度計算方法

本文構建的模型參數(shù)， 主要包括泥質(zhì)體積分數(shù)（Vsh）、石英長石砂巖體積分數(shù)（VQFM）、碳酸鹽巖體積分數(shù)（Vcar）、黃鐵礦體積分數(shù)（Vpyr）、校正后干酪根體積分數(shù)（Vkero），還有孔隙度（?）。 通過最優(yōu)化反演計算出這些礦物體積分數(shù)之后，則總孔隙度?t為

基于本次研究所涉及的基于閾值分割算法所求得的有機質(zhì)孔面孔率和利用核磁實驗所求得的有機質(zhì)孔隙度占比，可分別得到不同的計算模型（見式（9）、式（10）、式（11））。

式中：β 為圖像處理獲得的不同小層的有機質(zhì)孔面孔率均值；ε 為核磁實驗確定的有機質(zhì)孔占比；?organic為有機質(zhì)孔孔隙度；?inorganic為無機質(zhì)孔孔隙度。

3.4 有機質(zhì)孔隙度計算模型實例與分析

利用圖像分割算法和基于最優(yōu)化反演與核磁實驗結合的有機質(zhì)孔孔隙度計算式（9）、式（10），分別得到了X1 井（見圖12）、X2 井（見圖13）的計算有機質(zhì)孔孔隙度和核磁實驗有機質(zhì)孔孔隙度。

圖12 X1 井巖心實驗測得的有機質(zhì)孔孔隙度與測井計算模型結果對比Fig.12 Comparison of organic porosity measured by core experiment with results of logging calculation model of Well X1

圖13 X2 井巖心實驗測得的有機質(zhì)孔孔隙度與測井計算模型結果對比Fig.13 Comparison of organic porosity measured by core experiment with results of logging calculation model of Well X2

2 條數(shù)據(jù)曲線對比可以看出，數(shù)據(jù)曲線基本吻合，只有局部有所偏差。 造成此現(xiàn)象的原因可能是改進Canny-Otsu 算法所處理的圖像樣本不足。 計算的實驗有機質(zhì)孔隙度和核磁實驗分析的有機質(zhì)孔隙度的絕對誤差均在10%以內(nèi)，X1 井樣品的平均絕對誤差為6.30%，X2 井樣品的平均絕對誤差為5.02%。

4 結論

1）通過對掃描電鏡圖像進行圖像處理與分析，構建出適用于SEM 圖像的改進Canny-Otsu 算法。 分別對X1 和X2 井的樣品圖像進行圖像處理，計算出面孔率為11.68%～21.72%，平均面孔率為15.55%。 飽和鹽水、 飽和油的核磁對比實驗得出，X1 井和X2 井樣品的有機質(zhì)孔、無機質(zhì)孔比例約55∶45。

2）依據(jù)最優(yōu)化反演計算得到了儲層的有機質(zhì)骨架干酪根及總孔隙度， 結合圖像分割算法得到的研究區(qū)面孔率與核磁實驗得到的有機質(zhì)孔和無機質(zhì)孔比例，構建了基于研究區(qū)面孔率的有機質(zhì)孔孔隙度計算模型和基于有機質(zhì)孔占比的有機質(zhì)孔孔隙度計算模型。

3）基于圖像分割算法的有機質(zhì)孔孔隙度計算方法和基于核磁實驗的有機質(zhì)孔孔隙度計算方法所得到的有機質(zhì)孔孔隙度絕對誤差均在10%以內(nèi)。 這表明，本文的研究成果為頁巖氣儲層有機質(zhì)孔隙度測井計算提供了一種可借鑒的方法和思路， 為海相頁巖氣儲層甜點評價和預測奠定了基礎。