碳酸鹽巖儲層自適應(yīng)模型常規(guī)測井滲透率預(yù)測

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

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

0 引言

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

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

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

1 方法原理

Kozeny[6]和Carman[7]在前人研究的基礎(chǔ)上根據(jù)滲流力學(xué)原理建立了Kozeny-Carman方程，定量地描述孔隙度與滲透率的復(fù)雜關(guān)系

(1)

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

H=Fτ2

(2)

式中：H為卡茲尼常數(shù)，在實際應(yīng)用中通常取值范圍是5～100。

為實現(xiàn)對孔隙結(jié)構(gòu)的定量評定，引入流動單元指數(shù)(FZI)，其計算公式為

(3)

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

(4)

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

K=1014.24(FZI*)2φ

(5)

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

DRT=ROUND(2lnFZI+C)

(6)

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

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

設(shè)定IPSO算法中粒子和種群為

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

(7)

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

(8)

(9)

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

圖1 自適應(yīng)模型的計算流程

2 實驗驗證

2.1 數(shù)據(jù)來源及實驗設(shè)計

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

圖2 實驗?zāi)康膶佣蝺宇愋?/p>

通過收集整理現(xiàn)有資料，共計獲取研究區(qū)6口井目的層段425個巖心樣本的孔隙度、滲透率及常規(guī)測井?dāng)?shù)據(jù)。孔隙度和滲透率數(shù)據(jù)均由氣測實驗獲得，其中孔隙度分布范圍為0.30%～12.39%，滲透率分布范圍為0.001～216.236mD。

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

表1 樣本數(shù)據(jù)(部分)

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

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

圖3 實驗?zāi)康膶佣稳⌒姆治隹紫抖扰c滲透率關(guān)系

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

需要注意的是，式(7)中C值取動態(tài)調(diào)試后的最優(yōu)值10.151，該常數(shù)會直接影響DRT的劃分?jǐn)?shù)量及回歸方程的擬合度R2，間接地影響最終計算結(jié)果。

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

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

圖4 實驗?zāi)康膶佣稳⌒臉颖綝RT劃分后不同測井曲線參數(shù)概率密度函數(shù)分布

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

表3 預(yù)測模型主要參數(shù)設(shè)置

2.2 實驗1

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

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

圖5 實驗1靜態(tài)模式下由FZI確定的C值與DRT劃分?jǐn)?shù)量及MSE的關(guān)系

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

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

圖6 實驗1機器學(xué)習(xí)模型與物理模型自適應(yīng)動態(tài)優(yōu)化過程及優(yōu)化結(jié)果

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

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

圖7 實驗1動態(tài)優(yōu)化后機器學(xué)習(xí)模型預(yù)測結(jié)果對比

2.3 實驗2

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

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

圖8 實驗2靜態(tài)模式下由FZI*確定的C值與DRT劃分?jǐn)?shù)量及MSE的關(guān)系

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

圖9 實驗2機器學(xué)習(xí)模型與物理模型動態(tài)自適應(yīng)優(yōu)化過程及優(yōu)化結(jié)果

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

圖10 實驗2動態(tài)優(yōu)化機器學(xué)習(xí)模型預(yù)測結(jié)果對比

2.4 綜合分析及數(shù)據(jù)驗證

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

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

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

為驗證本文自適應(yīng)模型在實際井中的應(yīng)用情況，利用IPSO-XGBoost算法單獨對研究區(qū)M9井目的層段孔隙度進行了連續(xù)預(yù)測，并對M9井目的層段儲層進行了常規(guī)測井綜合評價(圖12)。基于孔隙度及滲透率預(yù)測結(jié)果，根據(jù)有效儲層下限判斷原則在5420～5510m分析得到13個小層。常規(guī)測井解釋結(jié)果與射孔區(qū)段實際產(chǎn)氣量對比分析表明，本文建立的自適應(yīng)模型能夠滿足目的層段儲層評價的要求。

圖12 M9井目的層段5420～5510m測井解釋結(jié)果與產(chǎn)剖面對比

3 結(jié)論與認(rèn)識

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

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

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

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

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