基于貝葉斯推斷的產量遞減綜合預測新模型

王軍磊 位云生 齊亞東 倪 佳 于 偉 袁 賀 朱漢卿 雷丹鳳

1.中國石油勘探開發研究院 2.中國石油西南油氣田公司頁巖氣研究院 3.德克薩斯大學奧斯汀分校

0 引言

油氣井生產數據分析或產量遞減分析是指油氣藏進入投產階段以后，通過使用合理模型擬合歷史動態數據、預測未來產量，并為早期氣田產能建設和后期方案調整優化提供可靠信息的技術方法。在衰竭式油氣藏開采過程中，生產井呈現緩慢的產量遞減特征，根據建模原理將產量遞減模型分為理論類（如生產數據分析方法，即RTA）[1-3]，經驗類[4-7]（如各類經驗遞減模型）以及數據驅動類（基于機器學習、神經網絡算法等建模）[8-12]3類方法。其中，理論類方法雖然具有明確的物理意義，但難以處理好物理模型的全因素假設和模型模擬效率間的矛盾；經驗類方法在大量生產數據分析基礎上建立經驗式產量遞減模型，但沒有嚴格的滲流理論支撐、適用性差；數據驅動類方法更多的是反映數據到數據的映射關系，預測精度直接取決于訓練數據的質量及算法的適宜性。需要明確的是，以上方法原則上適用于頁巖、致密砂巖等不同類型的非常規油氣藏/井動態數據分析。

盲目使用以上方法應分析實際生產數據會極大增加產量遞減分析預測的不確定性[13-14]，如何合理量化這種不確定性是關鍵，頻率學派和貝葉斯學派是目前2種主要流派。頻率學派依據最大似然估計原理，主要有3種方法：①在產量數據庫基礎上逐井擬合獲得遞減模型參數的概率分布模型，結合隨機模擬獲得氣井產量的概率性預測[15]；②對特定井的生產數據進行有放回地多次抽樣以形成多數據組合，通過擬合獲得產量預測的概率性分布[16]；③假定模型參數概率密度，通過隨機參數抽樣預測產量，根據目標函數設定以篩選概率性產量預測結果[17]。貝葉斯學派則是以最大后驗估計為基礎，克服頻率學派模擬時的隨機性，Gong等[18]首次利用貝葉斯理論量化產量遞減分析的不確定，Fulford等[19]以流態識別為依據，通過使用瞬時雙曲模型和參數分布似然函數改進了Gong提出的方法；Paryani等[20]使用近似的復雜似然函數建立近似貝葉斯概率性方法，大幅度提高了貝葉斯模擬效率，通過多模型協同約束以降低預測不確定性；Holanda等[21]建立具有物理意義的全流態演化模型，結合隨機最大似然原理更新似然函數中的協方差矩陣，有效提高了貝葉斯方法的運行效率。

需要指出的是，以上研究均是針對某種特定的模型，所選模型類型本身也會影響預測結果的不確定性。傳統方法中僅通過對比適用條件和歷史擬合效果而優選“好”模型是不科學的[22-23]，適用性好、擬合效果好的模型并不代表是真正意義上的“好”模型，而是代表該模型是“好”模型的概率（可視為權重）為1，人為區分“好”模型和“差”模型相當于指定了模型概率（即好模型概率為1，差模型概率為0），大大增加產量預測的風險性。本文以6種常見經驗模型為候選模型，使用“貝葉斯概率”量化雙重的不確定性，即單一模型EUR預測的不確定性和多個模型選擇權重的不確定性，充分發揮多種模型間相互兼容、相互制約的技術優勢，有效提高單井EUR預測可信度。該方法不僅適用于經驗式模型，也適用于解析模型和數值模型，具有良好的可擴展性。綜合模型可以有效改進單一模型EUR預測結果的不確定性和風險性，為我國非常規油氣開發提供有益借鑒。

1 產量遞減模型適用性評價

產量遞減模型原理是回歸歷史生產數據。實際應用時，根據實測數據利用各種優化算法，獲得最優擬合效果（即優選最“好”模型），通過反演計算模型參數（向量）以獲得確定性的預測結果。

經驗式產量遞減模型是建立在同一種或幾種流動狀態下動態分析基礎上的，即每種模型對應特定的流動狀態。本文給出了Arps型、冪律指數型（PLE）、擴展指數型（SEPD）、Duong型和邏輯增長型（LGM）等適用于不同流態的6種遞減模型，根據量綱一致性原理將模型轉化為無量綱形式，以便使用典型圖版擬合法進行數據分析。相應模型數學表達式和出處見表1。

表1 常見的6種經驗式產量遞減模型表

以北美地區某口致密氣壓裂井生產數據為例進行分析[26]，相應的地質工程參數為：原始地層壓力為19.07 MPa，地層溫度為30 ℃，儲層有效厚度為4.5 m，孔隙度為8.5%，含氣飽和度為80%，原始地層壓力下氣體黏度（μgi）為0.021 6 mPa·s、偏差因子（Zgi）為0.776 1。建立RTA解析模型擬合生產數據，擬合后的裂縫長度為12 m，地層滲透率為2.67 mD，井控地質儲量（OGIP）為2 049×104m3，設定20年生產周期，EUR預測值為1 185×104m3（該值可視為EUR真實值）。使用經驗式模型進行典型圖版擬合分析（圖1），氣井先后經歷雙線性流+線性流+擬穩態流三個連續生產階段。結果顯示：除Duong模型外，其余模型均可獲得較好的歷史擬合效果，EUR預測結果 為 ：EURArps=1 314×104m3、EURPLE=1 149×104m3、EURSEPD=1 016×104m3、EURLGM=1 387×104m3、EURFDC=1 572×104m3、EURDuong=2 625×104m3。

圖1 產量遞減模型典型圖版擬合效果圖

由EUR預測結果可知，該算例中除Duong模型不適用外，其余5種模型EUR結果介于1 149×104～1 572×104m3之間，包括EUR真實值（EURRTA=1 185×104m3）。其中，FDC模型結果最為樂觀，SEPD模型結果最為保守，PLE模型結果最為接近，Arps和LGM模型結果近似。根據模型適用性條件分析可知：①Arps模型（圖1-a）僅適用某種特定流態，當2＜b＜4時適用雙線性流，當b=2時適用線性流，當b＜1時適用擬穩態流[27]，該井擬穩態生產周期相對整個周期較長，EUR預測值較合理但仍偏高；②PLE模型（圖1-b）通過將遞減指數b修正為關于時間的變量，適用從（雙）線性到擬穩態的整個流態變化過程，該井流態較為清晰，通過擬合該模型獲得最為合理的EUR預測值；③SEPD模型（圖1-c）適用任意非穩態生產數據，但在后期累積產量趨近于界限值，該井較長的擬穩態生產歷史使得模型低估了EUR；④LGM模型（圖1-d）難以擬合整個流態過程，擬合擬穩態過程時需要更新模型參數，適用條件和預測結果與Arps模型類似；⑤Duong模型（圖1-e）適用于（雙）線性流數據，該井后期的擬穩態數據仍解釋為線性流，因而該模型EUR預測值偏大；⑥FDC模型（圖1-f）能較好擬合各種流態，包括過渡流等，該井在擬穩態階段近似指數遞減，而模型衰減速率要低于指數衰減，導致EUR預測結果較高。因此，判斷某種模型是否適用于特定的流動狀態，應以流態識別為依據的同時結合模型適用條件進行具體分析[13]。

總的來看，該口井經歷的流態較多（即生產周期相對較長）、氣井數據質量較高（生產制度穩定、數據噪音小、流態清晰），即使獲得了最優歷史擬合，不同模型EUR預測結果仍有較大不確定性。不確定性主要來自兩個方面：①單個模型本身參數擬合的多解性；②不同模型適用條件的差異性，而且模型參數越多、模型差異越大，預測EUR的不確定性越高。如何用“概率”思維量化雙重不確定性是擴展經驗遞減模型適用范圍的關鍵。

2 貝葉斯推斷原理

貝葉斯原理主要用以描述模型參數、實測數據和模型輸出值之間的關聯，其將概率看成對事件發生的信心，并且保留不確定性。將多維參數向量（θ）視為隨機變量，則存在與θ相關的概率分布函數，給定任何θ取值都能得到相應的概率值。經典貝葉斯推斷定理可表述為根據先驗分布和可能性（似然）分布獲得后驗分布的過程，滿足如下公式[23]：

式中q表示觀測數據，對應歷史生產數據，如產量或壓力等；θ表示模型參數向量；p(θ|q)表示給定實測數據時關于模型參數向量（θ）的后驗分布概率；p(q|θ)表示似然函數；p(θ)表示關于模型參數向量（θ）的先驗分布；p(q)表示邊際似然函數。

對于邊際似然函數，本文分為兩種情況進行計算：

第一種連續型變量（即模型參數），對應的分布為概率密度函數，邊際似然函數為積分形式：

第二種離散型變量（即模型類型），對應的分布為概率質量函數，邊際似然函數為求和形式：

假設損失函數滿足標準正態分布N(0,σ2)，似然函數為：

理論上先驗概率分布p(θ)可以任意給定，通常假設隨機變量均勻分布，而最終（穩定）的后驗分布通過式（1）～（5）結合實測數據確定。

2.1 離散型變量隨機模擬

將模型類型視為離散型變量，利用概率關系聯合多個單一模型構建綜合模型。根據貝葉斯推理，式（1）第j個模型mj為“最佳”模型的后驗概率為

式中θj表示第j個模型的參數向量。

離散型參數的邊際似然函數p(q)的積分形式可以寫為離散求和形式。與連續型變量不同，式（6）可以直接計算。

根據先驗分布特征，假設各個模型的先驗概率分布相同，即

則式（6）改寫為：

根據似然函數定義，將式（5）代入式（8），特定模型參數向量條件下的后驗概率為：

式中p(θj,mj|q)表示給定觀測數據q下第j個模型參數θj的后驗分布概率，可以作為第j個模型預測結果的權重；mj表示第j個模型；M表示模型總個數；εij表示第j個模型中第i組數據間的隨機誤差；σj表示第j個模型的數據標準差。

2.2 連續型變量隨機模擬

借助馬爾科夫鏈—蒙特卡洛（MCMC）方法，即從近似分布中獲取一系列的模型參數采樣點，通過校正采樣點獲得后驗分布更好的近似后驗分布。由于隨機變量的后驗分布是未知，需要通過從另一個概率分布中抽樣獲得。MH（Metropolis—Hastings）方法是實現MCMC過程的經典算法，通過建立滿足細致平穩方程的轉移概率矩陣，基于“拒絕采樣”原則沿著馬爾科夫鏈不斷逼近平穩分布，即從任一狀態出發通過不斷進行狀態轉移最終收斂到平穩分布。

MH算法的關鍵是建議分布p(θ|θ*)和接受率分布（α），其中θ表示隨機變量，θ*表示給定的隨機變量。這里設定接受率（α），即接受θ=θ*的概率為α，不接受概率為1-α。設定接受率大于1時為1，規整化處理后滿足如下公式：

通常情況下建議分布滿足對稱假設，即p(θ|θ*)=p(θ*|θ)，式（10）可進一步簡化為 ：

MH算法假設參數間相互獨立導致接受率過低，這里使用AM算法增加變量間的自相關性以提高接受率，基本原理為：每次迭代過程中，在上一步協方差矩陣（Ci）基礎上生成并更新建議分布[27]：

式中i0表示初始周期（更新協方差矩陣的周期）；ε'表示極小值以確保協方差矩陣非奇異，一般取10-12；Id表示d維單位矩陣；sd表示取決于變量維數的換算因子，一般取sd=2.42/d；C0表示初始正定矩陣。

以一維隨機變量模型為例對比MH和AM算法的采樣模擬效果。假設變量目標的概率密度（PD，取值范圍介于0～1）分布滿足正態分布加權平均形式，對應的后驗概率π為：

式中ω1=0.3，ω2=0.7，μ1=0，μ2=10，σ1=2，σ2=2。模擬結果如圖2，可以看出AM算法接受率更高（α=0.354），而且在采樣過程中隨機變量的取值范圍更廣、更均勻，可以獲得更好的目標分布。隨機變量維數越高，AM算法優勢越顯著。

圖2 MH算法與AM算法模擬結果對比圖

3 綜合模型產量遞減不確定性分析

3.1 單一模型分析

量化單一模型擬合結果的不確定性分為3個步驟：①根據最小二乘法原理，即實測數據與模型預測值間的損失函數最小，獲得確定性“最優”模型參數θ0；②求解損失函數最小時對應的近似參數向量（）及對應的協方差矩陣；③以最優擬合參數（θ0）及協方差矩陣為多維高斯概率模型的初始值代入式（12），并根據建議概率分布N(θi,Ci)抽樣參數向量θ*，進而調用MCMC模擬擬合或預測結果的不確定性。

以FDC模型為例說明單一模型預測的不確定性。迭代次數設定為10 000次、每100次更新協方差矩陣（i0=100）。為了提高模擬效率，將產量遞減模型處理為關于最高產量和時間的形式（即三維參數向量降維為二維）：

式中qi表示產量最高值（qi=11.74×104m3/d）；ti表示產量最高值對應的時間（ti=1 d）；Eα,1( ) 表示Mittag-Leffler函數；α、γ表示待擬合無量綱參數。

根據Mittag-Leffler函數的性質[25]，可以獲得最優模型參數向量{θ0=[α0,γ0]}對應的雅克比行列式：

式中J表示雅克比行列式；t表示時間序列向量，d；α0、γ0表示最優擬合參數。

獲得的參數向量最優解（α=0.836 1，γ=0.031 0）與手動圖版擬合結果基本一致（圖1-f），二維參數向量的接受率為36.57%，滿足接受率的要求。模型參數向量的后驗概率分布如圖3，參數分布較為集中（α和γ分布區間分別介于0.72～0.95和0.015～0.061），參數間的Pearson's r（皮爾遜相關系數）為-0.962 3，近似反向線性相關。

圖3 遞減模型參數的后驗概率分布及雙參數間關聯度圖

根據模型參數后驗分布，對每種參數組合下的模型進行產量擬合及EUR預測，不確定性模擬結果如圖4所示，其中圖4-a為日產氣量（q），圖4-b為累積產氣量（Gp），圖4-c為EUR概率密度分布（PD）及累積分布（CD）。該井EUR預測區間介于1 462×104～2 133×104m3（最低值—最高值），80%EUR置信域（P10～P90）區間介于1 537×104～1 893×104m3，EURP50=1 629×104m3。圖版最優擬合解（EURFDC=1 572×104m3）位于80%置信域內，近似于P50值，說明確定性圖版擬合法和不確定性法的評價結果具有較高相融性。但FDC模型80%置信域并不包括EUR真實值（1 185×104m3），說明選擇該模型預測EUR的可信度較低（即風險性較高）。

圖4 遞減模型不確定產量擬合及EUR預測結果圖

3.2 綜合模型分析

在單一模型不確定性分析基礎上，根據貝葉斯推斷原理和MCMC采樣，建立基于貝葉斯概率的不確定性產量遞減分析綜合模型，以降低單一模型預測的風險性。建模流程如圖5所示（其中白色框架為MCMC-AM采樣內容，黃色框架為貝葉斯推斷原理內容），每次迭代（i=1～N）時計算每種模型（j=1～M）的后驗（貝葉斯）概率值p(θji,mj|q)記為Pji；以觀測數據（q）為約束，基于模型mj使用參數向量（θji）計算得EUR值記為fji，分布如表2所示。

圖5 綜合模型建模工作流程圖

表2 迭代過程中不同模型的權重值分布表

相應地，第i次迭代時綜合模型的EURi預測值為：

式中Pji表示第i次迭代時第j個模型的后驗概率值；fji表示基于第i次迭代時第j個模型計算的EUR值；M表示模型個數。

對于給定觀測數據（q），相對于其他模型，第j個模型是“優選”模型的相對可信度（即優選模型概率）：

式中mj表示第j個模型；N表示迭代次數；p(q) 表示在迭代過程中獨立樣本，即p(q)=1/N。

以算例數據為例說明工作流程，分兩種情況：

第一種情況：分別對其余4種模型進行不確定性評價（Duong模型除外）。各模型參數向量的后驗概率分布如圖6所示，Arps模型（圖6-a）參數分布范圍最大、接受率最低（30.26%）；而PLE模型（圖6-b）參數分布范圍最小、接受率最高（36.35%），其余兩種模型介于兩者之間。

圖6 不同模型間模型參數向量后驗概率分布圖

圖7-a～b為Arps和PLE兩種模型產量遞減的80%置信域區間（P10～P90），PLE模型具有比Arps模型更小的產量遞減置信域區間，這與模型參數的分布范圍集中程度相關；圖7-c為不同模型預測EUR的累積概率分布，其中Arps模型分布區間最大，PLE模型分布范圍最小，進一步驗證了產量遞減置信域區間范圍。

圖7 Arps模型與PLE模型80%置信域的產量遞減區間圖

第二種情況：評價綜合模型EUR預測結果的不確定性。由于每種模型接受率不同導致對應的預測EUR值個數不同，假設EUR樣本數據獨立同分布，采用自助法（bootstrap）進行有放回抽樣[16]，每種模型獲得相同數量的EUR樣本值（表2）。使用式（17）計算各模型的相對可信度，Arps、PLE、SEPD、LGM和FDC模型依次為13.2%、37.6%、15.9%、21.7%和11.6%。其中，PLE模型相對可信度最高，原因在于該模型能夠很好地擬合連續的雙線性+線性+擬穩態三個流態數據，模型參數約束性較好（即未知參數與流態數量相匹配），預測結果更接近真實值。

圖8對比了不同單一模型和綜合模型EUR不確定性評價結果。在單一模型中，FDC模型80%EUR置信域均高于EUR真實值，SEPD模型80%EUR置信域均低于EUR真實值，其余模型80%EUR置信域包括了EUR真實值。此外，Arps模型EUR預測的不確定性最大（標準差為269.11，對應80%置信域區間最大），PLE模型EUR預測的不確定性最?。藴什顬?47.56，對應80%置信域區間最小），PLE模型EUR預測的可靠性最強，EUR真實值位于50%EUR置信域（P25～P75）、對應的EUR的P50值（即中位線）與真實值最接近，這與圖7-b預測結果以及模型相對可信度（PLE為37.6%）解釋相一致。相對于特定單一模型，綜合模型的80%EUR置信域范圍介于各單一模型置信域范圍之間，抵消了單一模型高估/低估預測的風險，而且EUR真實值位于20%置信域內（P40～P60），相對于PLE模型更為接近。原因為：①PLE模型殘差的方差值最小[據式（5）、（9），PLE 模型對應的p(θj,mj|q)值最大]，導致EURPLE預測值所占權重最大，FDC模型則情況相反；②綜合模型同時考慮了其他模型的概率影響，避免了只選擇特定模型所帶來的“人為”風險。

圖8 不同模型間的EUR預測結果置信域分布范圍圖

按照以上分析方法及流程，選取川南地區25口頁巖氣井生產數據進行分析，生產歷史從14個月到58個月不等，井間經歷的流態差異性較大。井間各個模型成為“好”模型的優選概率模擬結果如圖9，其中為優選模型的SEPD涉及井數為10口、PLE模型為7口、LGM模型為4口、PDC模型為3口、Duong模型為1口、Arps模型為0口，原因在于目前生產周期內氣井均處于非穩態生產階段，Arps模型的適用條件均弱于其他模型。同時可以得到沒有一種模型適用于所有井，當模型間相對可信度值較為接近時，實際應用時需要參考綜合模型EUR預測結果；當某種特定模型相對可信度遠高于其他模型時（如14#、18#、22#井），說明該模型成為優選模型的概率高，在綜合模型EUR預測結果基礎上可重點參考該模型。

圖9 25口頁巖氣井不同模型優選模型概率圖

圖10-a給出了使用單一及綜合模型的25口井EUR預測中值，可以得到井間不同模型EUR預測值差異性較大，主要原因在于井間的流態演化不同、數據質量不同，也證實了選擇單一優選模型時的不確定性。圖10-b給出了25口井的綜合模型不確定性EUR預測結果，所有井的均值都位于50%的EUR置信域內，在集合了各個模型優勢基礎之上有效降低了使用單一優選模型所帶來的風險，對于EUR置信域范圍較小的井（如1#、11#），對應的生產歷史較長、流態特征清晰且個數較多、擬合效果好，各模型EUR置信域近似，這種井EUR預測結果可靠性強；對于EUR置信域范圍較大的井（如3#、6#），其數據質量噪音較大、擬合效果差，各模型EUR置信域差異較大，這種井EUR預測結果可靠性較差。

圖10 25口頁巖氣井單一模型及綜合模型EUR預測結果

由此可見，綜合模型EUR預測值的置信域具有合理分布區間，同時能夠保證EUR真實值位于更小的置信域內（即獲得EUR真實值的概率更高），降低了使用單一優選模型帶來的風險。

4 結論

1）經驗產量遞減模型僅適用特定的某個或多個流動狀態，超出條件即使獲得較好的擬合效果也會導致不可靠的EUR預測結果；受制于模型參數、模型選擇，EUR預測值具有不確定性，動態數據相對生產周期越短、模型參數越多、模型適用性越差，不確定性越高。

2）使用馬爾科夫鏈—蒙特卡洛算法可以獲得各個單一模型參數的后驗概率分布及EUR不確定性預測區間，但模型預測的EUR置信域可能不包括EUR真實值，這取決于是否選定合適的模型。

3）以各單一模型參數的貝葉斯后驗概率值為權重，加權平均獲得的綜合模型不確定性EUR預測置信域分布區間更為合理，提高了獲得EUR真實值的概率，生產數據質量越差、優選模型越困難，使用單一模型預測EUR的不確定性和風險性越高，對應的綜合模型法優勢越明顯。