碳酸鹽巖儲層自適應模型常規測井滲透率預測

徐鵬宇 周懷來 官俊潔 陳 莉 周 捷 宋雨莎

(①成都理工大學地球物理學院，四川成都 610059； ②成都理工大學“油氣藏地質及開發工程”國家重點實驗室，四川成都 610059)

0 引言

準確預測滲透率是碳酸鹽巖儲層常規測井評價的重點內容。碳酸鹽巖儲層通常具有原生孔隙、次生孔隙及各種裂縫共同發育的特點，因受沉積環境、成巖作用及含氣水飽和度等多種因素的影響，非均質性較強，常規地球物理測井解釋存在多解性，需借助成像測井作為參照[1-3]。目前常規測井解釋滲透率主要利用物理模型[4]和擬合模型[5]，常用的物理模型有Kozeny-Carman[6-7]、指數型和冪律型[8]等，常用的擬合模型有多元線性回歸和逐步回歸等模型。在實際應用中，物理模型和擬合模型互為補充，其間的耦合程度還需進一步探索，如物理模型與機器學習模型對滲透率預測精度的影響程度以及離散巖石類型(Discrete Rock Type, DRT)[9]的劃分原則等，都是影響碳酸鹽巖非均質儲層滲透率解釋精度的因素。

為更科學地利用常規測井數據建立儲層參數的擬合模型，人們應用機器學習模型替代儲層參數與常規測井數據間的非線性表達式。目前性能較好的機器學習模型主要有支持向量機(Support Vector Machine, SVM)[10-11]、隨機森林(Random Forest, RF)[12-14]和XGBoost[15-17]等。Zhang等[11]采用6種模型對碳酸鹽巖儲層流動單元指數(Flow Zone Indicator, FZI)進行預測并間接地計算滲透率，認為粒子群優化算法(Particle Swarm Optimization,PSO)-SVM模型性能最佳，測試集決定系數高達0.869； 郭建宏等[13]利用RF模型預測煤層氣含量； 孫予舒等[15]利用XGBoost模型對復雜碳酸鹽巖巖性進行識別。XGBoost[18]在眾多機器學習模型中脫穎而出，得益于其算法和運行方式上的一些改進：在計算過程中加入正則化項，能有效防止過擬合； 在運行方式上開辟多個線程并行處理，能在大規模計算時提高CPU利用率。

面對傳統機器學習模型錯綜復雜的參數調試問題，人們引入遺傳算法(Genetic Algorithm, GA)[19]、粒子群優化算法[20]、模擬退火(Simulated Annealing, SA)算法[21]等替代繁瑣的人工調試手段，同時提高模型的預測精度。谷宇峰等[17]利用PSO算法改進XGBoost并得到了可靠的滲透率解釋模型，但在求解復雜的碳酸鹽巖儲層滲透率非線性問題時效果卻不理想，當滲透率數量級超過10-2～10mD時，孔隙結構更復雜，僅依靠機器學習模型直接預測滲透率難以取得較好的效果。為此，本文從優化模型參數的角度出發，利用改進的粒子群優化算法(IPSO)[22]同時優化物理模型參數及機器學習模型參數，進而提出滲透率預測的自適應模型，并闡述該算法原理及其在非均質碳酸鹽巖儲層滲透率預測中的應用。

1 方法原理

Kozeny[6]和Carman[7]在前人研究的基礎上根據滲流力學原理建立了Kozeny-Carman方程，定量地描述孔隙度與滲透率的復雜關系

(1)

式中：K為滲透率；F為形狀系數，圓柱形管時取2；τ為孔隙介質迂曲度，通常取1.5～3.0；φ是有效總孔隙度；S是單位體積顆粒比表面積，與φ值共同由氣測實驗確定。令

H=Fτ2

(2)

式中：H為卡茲尼常數，在實際應用中通常取值范圍是5～100。

為實現對孔隙結構的定量評定，引入流動單元指數(FZI)，其計算公式為

(3)

Amaefule等[23]根據FZI對Kozeny-Carman方程進行變形，得到

(4)

Mirzaei-Paiaman等[24-26]提出一種修正的流動單元指數(FZI*)，并得到

K=1014.24(FZI*)2φ

(5)

Kharrat等[9]提出離散巖石類型(DRT)聚類算法，通過FZI或FZI*值劃分DRT，表達式為

DRT=ROUND(2lnFZI+C)

(6)

式中：C為經驗系數，由FZI或FZI*值的正態分布特性確定，初始值取10.6； ROUND(·)表示取整。

在建立儲層參數模型時需要設定許多經驗參數，如機器學習模型中SVM的懲罰系數、RF的決策樹個數和XGBoost的最大回歸樹深度，物理模型中DRT聚類算法中的經驗系數C等。采用IPSO優化算法能有效地處理多參數極值求取問題，同時解決線性下降(Linear Decreasing Inertia Weight，LDIW)-PSO算法在迭代過程中全局搜索和局部搜索比例不變問題，通過該算法優化可以加強滲透率預測模型中擬合模型和物理模型的耦合程度。

設定IPSO算法中粒子和種群為

X=(x1,x2,···,xN)T∈RN×D

(7)

式中：X為種群，xN為第N個粒子，N為種群中粒子個數；D為優化參數個數，包含機器學習模型中的參數及式(6)中的C。在種群迭代尋優過程中，每個粒子xN通過個體極值和群體極值調整自身的速度和位置

(8)

(9)

綜合前述原理，基于Python語言編譯了孔、滲關系計算方程自動擬合算法。圖1為該自適應模型的計算流程，圖中MSE為均方誤差(Mean Squared Error)，根據測試集所有DRT聚類算法優化預測模型得出的滲透率與實測滲透率計算獲得，計算時對滲透率數據取對數。

圖1 自適應模型的計算流程

2 實驗驗證

2.1 數據來源及實驗設計

研究區目的層為川中北部高石梯—磨溪地區震旦系燈影組燈二段，該沉積期主要發育臺地、臺緣、斜坡及陸棚等沉積相帶，是目前四川盆地發現的最古老的碳酸鹽巖含氣層段[27-28]。根據分析取心及成像測井資料，將儲層分為孔隙型(圖2a)、裂縫—孔隙型(圖2b)、孔洞型(圖2c)及裂縫—孔洞型(圖2d)，裂縫-孔洞型儲層以斜交縫和水平縫為主，發育少量高角度縫，為該區氣藏提供了優質的滲流通道。在同一孔隙結構中，不同產狀裂縫的發育使滲透率呈冪級數增長，因此，得到準確碳酸鹽巖儲層滲透率解釋值的關鍵在于科學合理劃分孔滲等級進而確定儲層的DRT。

圖2 實驗目的層段儲層類型

通過收集整理現有資料，共計獲取研究區6口井目的層段425個巖心樣本的孔隙度、滲透率及常規測井數據?？紫抖群蜐B透率數據均由氣測實驗獲得，其中孔隙度分布范圍為0.30%～12.39%，滲透率分布范圍為0.001～216.236mD。

由式(4)和式(5)分別計算FZI和FZI*。每個樣本包含5種常規測井曲線，分別是聲波時差(AC)、補償中子曲線(CNL)、密度曲線(DEN)、自然伽馬曲線(GR)和原狀地層電阻率曲線(RT)。部分樣本數據見表1。

表1 樣本數據(部分)

對目的層實測孔隙度和滲透率數據進行交會分析，并建立指數模型進行擬合，其擬合度R2僅為0.483(圖3a)，可見單一回歸方程不滿足解釋精度。

根據式(4)和式(7)對復雜的孔滲關系進行劃分得到16種DRT劃分結果，將首尾樣本數量小于4的DRT合并入相鄰集合(圖3b)，對劃分后的樣本建立回歸方程進行擬合，部分計算結果見表2。

圖3 實驗目的層段取心分析孔隙度與滲透率關系

表2 孔隙度滲透率計算公式(部分)

需要注意的是，式(7)中C值取動態調試后的最優值10.151，該常數會直接影響DRT的劃分數量及回歸方程的擬合度R2，間接地影響最終計算結果。

C值的確定與實測樣本FZI值的正態分布特性密切相關，圖4展示了目的層段取心樣本DRT劃分后5種常規測井曲線參數的概率密度函數分布情況。由圖可見，劃分后各DRT類型對不同常規測井曲線參數的敏感程度也近似服從高斯正態分布[9]。

AC曲線主要反映的是原生孔隙，對次生孔隙不敏感，具體表現為少部分DRT劃分類型具有相似的AC概率密度函數(圖4a)，各DRT劃分后AC概率密度函數的均值不同。當孔隙結構中發育低角度裂縫或斜交縫時，AC曲線具有周波跳躍特征，造成其概率密度函數的方差變大。DRT聚類劃分的序列越高表明孔隙結構中裂縫越發育、滲流能力越好。圖4c中劃分后大部分DRT具有相似的DEN概率密度函數，均值和方差變化幅度小，表明不同DRT類型反映的是同一孔隙體積中原始孔隙、次生孔隙和裂縫的不同組合。圖4c和圖4d中DRT(21+20)對應的概率密度函數均明顯異于其他類型，通過完井報告、薄片鑒定等手段核實，為溶蝕孔洞中充填泥質所致。影響CNL和RT曲線的地質因素較多，如地層水礦化度、含氣水飽和度等，DRT劃分結果對CNL和RT概率密度函數不敏感(圖4b、圖4e)。

圖4 實驗目的層段取心樣本DRT劃分后不同測井曲線參數概率密度函數分布

針對Kozeny-Carman物理模型因基礎資料不足而無法應用的問題，本文采用機器學習模型替代并建立常規測井曲線參數與FZI或FZI*值的非線性表達式，通過IPSO算法迭代尋優并自適應地獲取機器學習模型參數和C值的最優解。文中采用SVM、RF和XGBoost等3種機器學習模型對FZI進行了預測，設計了2個實驗驗證該自適應模型的預測能力：實驗1為基于FZI間接預測滲透率； 實驗2為基于FZI*間接預測滲透率。預測模型的主要參數設置見表3。

表3 預測模型主要參數設置

2.2 實驗1

將425個樣本實驗數據按照7∶3的比例隨機劃分為2個集合，得到128個樣本的測試集和297個樣本的訓練集。根據式(4)計算全部樣本的FZI值，進而確定C值，再由式(6)根據不同C值計算DRT劃分數量，在靜態模式下探究C值對DRT劃分數量及滲透率MSE的影響(圖5)。需要說明的是，靜態模式不涉及機器學習模型及常規測井曲線參數。

由圖5a可知，不同C值對應的MSE呈現周期為1的性規律變化，表明C值僅在長度為1的閉區間內對MSE有較大的影響，并在C近似取10.82時獲得MSE極值0.033(圖5b)。圖5c分別展示了不同C值對應的合并前、后DRT劃分數量，合并前最大16類、最小15類，合并后最大為15類、最小為14類，表明不同C值會影響DRT劃分的數量。圖5b和圖5c對比表明，采取合并離散樣本建立物理模型的方式會增加MSE。另外，由于實驗時C值的步長為0.01，圖5b中的極值近似等于該組樣本測試集滲透率計算精度的“物理上限”。

圖5 實驗1靜態模式下由FZI確定的C值與DRT劃分數量及MSE的關系

因建模過程中測試集樣本的FZI值由機器學習模型與常規測井曲線建模預測獲得，故預測值和實測值必然存在偏差。樣本測試集滲透率MSE受多個參數的影響，需采用IPSO算法動態確定這些參數。因周期外的C值只影響DRT序列值的變化而不影響其他參數，為方便數據展示，C參數值取10～11(表3)，并將測試集預測滲透率與實測滲透率計算得到的MSE設定為適應度值并迭代尋優。圖6為3種機器學習模型與物理模型的自適應動態優化過程及最終優化結果。

相對于傳統的PSO算法，IPSO算法在迭代次數和尋優精度方面明顯處于優勢，3種自適應模型在迭代次數20以內的優化結果均超越PSO算法，甚至還有更大的尋優潛力(圖6a)。在擬合精度方面，XGBoost模型性能最佳，MSE僅為0.292，表明XGBoost模型耦合能力最強，能較敏感地響應物理模型。圖6b顯示3種自適應模型在迭代100次后C值存在多解性，表明不同機器學習模型與物理模型的耦合程度不同，僅依靠單一經驗參數C較難推廣。

圖6 實驗1機器學習模型與物理模型自適應動態優化過程及優化結果

需要說明的是，IPSO算法在迭代次數超過100后仍有上升空間，但已經能有效完成滲透率預測，從側面說明了IPSO迭代尋優算法的優越性。圖7為3種機器學習模型與物理模型最終自適應建模結果，不同機器學習模型的預測性能主要體現在其對10-3～1mD數量級范圍內滲透率樣本的擬合度效果。

IPSO-SVM的自適應建模結果(圖7a)顯示，該模型對10-3～1mD數量級范圍的滲透率擬合效果較差，這是由于物理模型在該數量級區間分布最緊密，且為碳酸鹽巖原生孔隙和溶蝕孔隙共同發育區域，非均質性極強，表明SVM模型與其耦合能力不足。IPSO-RF的自適應建模結果(圖7b)顯示，該模型提升了對滲透率10-3～1mD數量級的擬合度，MSE由0.401降至0.338。圖7c中建模結果顯示在該范圍內樣本點的擬合度最高，MSE進一步降至0.292。

圖7 實驗1動態優化后機器學習模型預測結果對比

2.3 實驗2

同樣，將425個樣本實驗數據劃分為2個集合，樣本構成要素同實驗1。由式(5)計算全部樣本的FZI*值，進而確定C值，并根據式(6)基于不同C值計算DRT類型，圖8為靜態模式下C值對DRT劃分數量及滲透率MSE的影響。

由圖8a可知，不同C值對應的MSE仍呈周期性規律變化，周期為1，并在C約為16.46時獲得MSE極值0.015(圖8b)。圖8c展示了不同C值計算的原始DRT劃分數量及合并后DRT數量，合并后最大為12類、最小為11類。對比實驗1靜態調試結果(圖5)，實驗2的DRT劃分數量和MSE均有所下降，且相對于MSE，DRT劃分數量下降幅度更大。在保證精度的情況下減少了計算工作量，表明FZI*計算方式更簡潔，在實際應用中更容易被推廣[24-26]。實驗2中最小MSE由原始DRT劃分數量13類合并至12類計算得到，而實驗1中最小MSE是根據原始DRT劃分數量15類計算得到，表明合并離散樣本的方式在DRT劃分數量較多的情況下會增加MSE，但在DRT數量較少的情況下相反，說明DRT劃分數量與近似“物理上限”沒有直接的聯系。為方便數據展示，C值取值范圍設定為16～17。圖9為3種機器學習模型與物理模型的自適應過程及最終優化結果。

圖8 實驗2靜態模式下由FZI*確定的C值與DRT劃分數量及MSE的關系

對比圖6實驗1的優化過程及優化結果，圖9a中3種自適應模型的MSE均有所下降，單從MSE曲線看，基于FZI*的物理模型與機器學習模型耦合程度更高。圖9b顯示3種自適應模型在迭代100次后C值仍存在多解性，再一次驗證了采用IPSO算法優化C值的必要性。

圖9 實驗2機器學習模型與物理模型動態自適應優化過程及優化結果

圖10為本實驗3種機器學習模型與物理模型自適應優化后的最終建模預測結果。相比實驗1動態建模結果(圖7)，實驗2動態建模結果擬合效果均有提升，模型中的樣本都更靠近45°線。仔細對比可見，IPSO-SVM模型中雖然大于0.1mD的樣本都基本靠近45°線，但對10-3～10-1mD范圍內的滲透率預測能力仍然不足，表明該模型在與物理模型的耦合方面難以推廣應用(圖10a)。對于IPSO-RF，有少部分10-3～10-1mD范圍的滲透率預測值與實測值不相符，樣本點散落在45°線外側，表明該模型能有效地適應不同的物理模型，且預測能力較穩定(圖10b)。對于IPSO-XGBoost，該范圍內的滲透率預測值與實測值均為同一數量級，樣本點更靠近45°線，表明該模型不僅穩定且預測能力最強(圖10c)。

圖10 實驗2動態優化機器學習模型預測結果對比

2.4 綜合分析及數據驗證

綜合分析實驗1、實驗2靜態調試結果，表明實際取心過程受多種非地質因素的影響，FZI及FZI*值難以呈現理論上的正態分布特征，故單依靠固定的經驗參數C值難以準確計算滲透率。綜合分析實驗1、實驗2動態調試結果，表明IPSO算法可結合樣本自身特性動態優化各種參數，自適應地耦合物理模型和擬合模型，最大程度地逼近“物理上限”。實驗1、實驗2綜合結果分析表明，物理模型決定預測精度上限，擬合模型決定預測精度下限，3種機器學習模型中，XGBoost模型預測性能最佳。

將建模以外的P1井79個取心樣本作為驗證集，驗證本文自適應模型的外推能力。圖11展示了IPSO-XGBoost算法基于不同預測方式對驗證集的計算結果，同樣以MSE評價預測精度。由圖可知，基于FZI和FZI*計算方式得到的結果均優于直接預測方式，經測算MSE由1.286分別降至0.386和0.296。從預測效果看，直接計算方式得到的結果數值跳躍誤差更大(黑色虛線框內)，表明自適應模型中DRT聚類算法有效地卡控了滲透率預測的數量級?；?種不同流動單元指數計算方法得到的結果均能為后續準確識別有效儲層提供合理支撐。

圖11 P1井目的層段取心樣本IPSD-XGBoost模型滲透率計算結果

為驗證本文自適應模型在實際井中的應用情況，利用IPSO-XGBoost算法單獨對研究區M9井目的層段孔隙度進行了連續預測，并對M9井目的層段儲層進行了常規測井綜合評價(圖12)?；诳紫抖燃皾B透率預測結果，根據有效儲層下限判斷原則在5420～5510m分析得到13個小層。常規測井解釋結果與射孔區段實際產氣量對比分析表明，本文建立的自適應模型能夠滿足目的層段儲層評價的要求。

圖12 M9井目的層段5420～5510m測井解釋結果與產剖面對比

3 結論與認識

本文從優化擬合模型、物理模型參數出發，利用IPSO算法同時優化物理模型參數及機器學習模型參數，得到碳酸鹽巖儲層滲透率預測的自適應模型，并基于流動單元指數的2種不同計算方式分別進行了實驗，得到以下認識：

(1)針對建模時物理模型基礎資料不充足問題，采用擬合模型作為補充，兩者耦合提高了模型的預測能力及適用性；

(2)DRT聚類算法中經驗系數C值在應用過程中存在多解性，依靠單一經驗值適用性差；

(3)IPSO算法在應對多參數求極值問題方面明顯優于傳統PSO算法，迭代次數更少，尋優精度更高，具有很強的推廣性；

(4)物理模型決定滲透率預測精度上限，擬合模型決定滲透率預測精度下限。