基于布谷鳥算法的儲層物性參數同步反演

劉仕友 段治川 周 凡 汪 銳

(①中海石油(中國)有限公司海南分公司，海南海口 570000； ②中國石油大學(華東)地球科學與技術學院，山東青島 266580)

0 引言

儲層物性參數(如孔隙度、泥質含量等)可為儲層評價提供依據。隨著油氣勘探開發目標日趨復雜，對儲層物性參數的反演精度的要求越來越高。

目前，常規的儲層物性參數預測大多是基于貝葉斯理論框架，利用物性參數的先驗分布信息，構建儲層彈性參數與物性參數之間的似然函數，尋找最大后驗概率解并將其作為反演結果。Mukerji等[1-2]在貝葉斯理論框架下，將統計巖石物理模型與疊前地震反演相結合，實現了巖性識別及流體預測。Eidsvik等[3]將巖石物理模型作為先驗信息，將地震數據與儲層參數之間的關系作為似然函數，實現了儲層及流體預測。Bachrach[4]基于Gassmann方程建立巖石物理模型，實現了孔隙度與含水飽和度的聯合反演。鄧繼新等[5]考慮孔隙度、飽和度與儲層巖石彈性特征的相關性，在巖石物理模型基礎上，結合貝葉斯統計反演，實現了孔隙度與含氣飽和度的聯合反演。Grana等[6]將統計巖石物理模型與貝葉斯反演相結合，在不同條件下對儲層參數進行概率估計，提高了反演精度。印興耀等[7]基于統計巖石物理模型、蒙特卡洛模擬及期望最大化算法，提出了基于彈性阻抗反演的儲層物性參數預測方法。Gui等[8]針對不同彈性參數之間精度的差異，構建了包含權重系數的儲層物性參數與彈性參數之間的加權關系，在貝葉斯框架下，結合馬爾科夫鏈蒙特卡洛(MCMC)隨機模擬技術，提高了儲層物性參數的反演效率。李志勇等[9]在貝葉斯反演框架下，對目標函數的梯度用泰勒公式近似展開，并將儲層彈性參數對物性參數的梯度用差分形式表示，通過共軛梯度算法實現了儲層物性參數與彈性參數的同步反演。楊培杰[10]通過Simon模型建立砂泥巖物性參數與彈性參數之間的定量關系，并基于貝葉斯理論構建多信息聯合約束的反演目標函數，聯合蒙特卡洛模擬與遺傳算法求解，實現了孔隙度與含水飽和度的同步反演。凌東明等[11]將巖相約束下線性化近似的Xu-White模型引入貝葉斯框架中，結合最小二乘優化算法實現了儲層物性參數反演。張佳佳等[12]對復雜巖石物理模型進行泰勒展開，得到線性化巖石物理模型，利用阻尼最小二乘算法求解精確解析解，預測了儲層物性參數。李坤等[13]提出差分進化MCMC隨機模型的相約束疊前地震概率化反演方法，實現了碎屑巖儲層物性參數同步預測。時磊等[14]針對致密砂巖儲層建模及參數反演難度大的問題，利用巖石物理敏感性分析作為指導，通過核貝葉斯判別法反演了致密砂巖儲層孔隙度、滲透率等參數。張健等[15]基于貝葉斯理論得到儲層彈性參數、物性參數及巖相后驗概率分布的解析表達式，由此提出一種構造約束聯合概率反演方法，引入構造信息和井信息，提高了反演結果的橫向連續性及分辨率。但地球物理反問題的非線性特征易導致常規算法在求解過程中陷入局部極值。因此，如何避免這種情況，是提高預測精度的關鍵點之一。

智能算法可避免陷入局部最優解，提高求解結果的準確性。近年來，部分學者將智能優化算法引入地球物理問題求解中。田景文等[16]針對模擬退火算法和BP算法的不足，將模擬退火算法與Powell算法有機組合并引入BP算法中，提高了算法的收斂速度和精度，實現了薄互層物性參數預測。Fattahi等[17]將智能優化方法，如粒子群算法優化的支持向量機回歸(SVR-PSO)、減法聚類及自適應神經模糊推理系統(ANFIS-SCM)等，與傳統貝葉斯方法預測的孔隙度、含水飽和度進行對比，驗證了智能優化方法的可行性與準確性。張冰等[18]基于巖石物理模型的反演流程，引入模擬退火優化粒子群算法，解決了多參數同時反演問題。邢凱[19]聯合變權系數混合高斯模型與布谷鳥算法，解決了在沒有考慮地質結構及初始巖相比例時得到的結果與實際地層分布不符合的問題，針對傳統隨機反演方法的不確定性，引入信息熵模型，并與布谷鳥混合高斯MCMC反演方法耦合，增強了反演結果的穩定性，提高了反演的精度。王鵬飛等[20]針對傳統全局尋優算法收斂速度慢、易陷入局部極值的缺點，將布谷鳥算法與單純形法、粒子群算法相結合，提高了大地電磁數據反演的精度。王峣鈞等[21]將布谷鳥算法與MCMC方法相融合，避免了在高維數據多參數同步反演時傳統MCMC方法因抽樣隨機性而陷入局部最優的問題，實現了不同巖相條件下儲層參數的精細描述。

將智能優化算法引入地球物理問題的求解時，部分智能優化算法由于較強的局部搜索能力，增強了在某一范圍內的搜索精度，但同時會導致算法的全局搜索能力不夠，在求解過程中易陷入局部最優問題中。

本文在前人研究的基礎上，以彈性阻抗與儲層物性參數關系為基礎，建立儲層物性參數反演所需的目標函數，針對地球物理反問題非線性特點以及部分智能優化算法全局搜索能力較弱的問題，引入一種新型元啟發式算法——布谷鳥算法，對反演目標函數進行求解，從而預測儲層物性參數。相比于其他智能優化算法，布谷鳥算法作為一種新的全局優化算法，其Levy飛行機制使算法有相對較高的概率生成大步長，增強了在求解過程中的全局搜索能力，較大概率避免了算法陷入局部極值，適用于解決地球物理反演非線性問題，同時也解決了部分優化算法全局搜索能力較弱的問題。通過模型試算及實際資料測試，驗證了本文方法能夠有效實現儲層物性參數預測。

1 方法原理

1.1 巖石物理模型構建

巖石物理方法是聯系儲層彈性參數與物性參數的橋梁，對地震正演和反演結果的定量解釋有著十分重要的意義[22-23]。經典巖石物理模型(如Xu-White模型等)因其自身假設條件(如巖石礦物組分、孔隙類型不夠豐富，巖石及組成成分均為各向同性、線性和彈性介質等)，無法兼具不同地區參數性質的多樣性，因而導致經典巖石物理模型在實際資料應用中受限。

在實際資料處理過程中，地震彈性參數與物性參數之間的關系可用曲線擬合方式刻畫。在誤差允許的條件下，可通過線性近似將二者的關系表示為

(1)

改寫為矩陣形式

(2)

式中：VP為縱波速度；VS為橫波速度；ρ為密度；φ為孔隙度；Vsh為泥質含量；Sw為含水飽和度； 系數ai、bi、ci(i=1,2,3)和常數項di(i=1,2,3)可由測井資料多元回歸得到。

當地震彈性參數與物性參數的關系無法近似為線性時，可通過多項式擬合的方法，將其表示為多元高階方程的形式

(3)

式中系數ai,j、bi,j、ci,j(i=1,2,…,n；j=1,2,3，…,n)和常數項dj(j=1,2,3)皆通過最小二乘擬合得出。通過式(1)或式(3)可建立彈性參數與物性參數之間的聯系，結合阻抗方程可進一步建立這種聯系。Connolly[24]對彈性阻抗方程進行了推導，具體為

(4)

式中：θ為入射角度；K=(VS/VP)2，一般令其等于0.25。

Whitcombe[25]在式(4)基礎上，引入歸一化常數，方程變換為

(5)

式中VP0、VS0、ρ0分別為測井資料中縱波速度、橫波速度、密度的均值。

對式(5)兩邊取對數，整理得

lnEI=AlnVP+BlnVS+Clnρ+D

(6)

其中A=1+tan2θ,B=-8Ksin2θ,C=1-4Ksin2θ，D=-tan2θ·lnVP0+8Ksin2θ·lnVS0+4Ksin2θ·lnρ0。

為了同時反演孔隙度、泥質含量、含水飽和度三個參數，需要引入三個角度的彈性阻抗[26]。將三個角度θ1、θ2、θ3所對應的彈性阻抗分別記為EI1、EI2、EI3，則式(6)可以表示為

(7)

將式(1)或式(3)代入式(7)，可得表征彈性阻抗與儲層物性參數之間關系的函數，即

[EI1,EI2,EI3]=fRPM(φ,Vsh,Sw)

(8)

引入隨機誤差項ε，則巖石物理模型可表示為

[EI1,EI2,EI3]=fRPM(φ,Vsh,Sw)+ε

(9)

至此，建立了彈性阻抗與物性參數之間的函數關系，為后續儲層物性參數的反演奠定了基礎。

1.2 布谷鳥搜索算法

布谷鳥搜索算法(Cuckoo Search Algorithm)是由Yang等[27-28]受布谷鳥孵育寄生現象啟發而提出的一種新型元啟發式算法。相對于其他啟發式算法，該算法具有參數少、全局搜索能力強的優點。布谷鳥搜索算法的主要思路是布谷鳥的巢寄生行為(育雛行為)與鳥類的Levy飛行機制兩種策略的結合。根據巢寄生行為，如果宿主發現布谷鳥的卵，則放棄該巢，并利用Levy飛行機制尋找新的位置。鳥類的Levy飛行行為作為一種典型的隨機游走機制，具有高度隨機性，由較長時間的短步長與較短時間的長步長組成。隨機游走步長滿足重尾分布，偶爾的長步長能擴大搜索范圍，增加種群多樣性，有效避免算法陷入局部極值。

在布谷鳥搜索算法中，設定了三個理想狀態：①每只布谷鳥一次只產一只蛋，并隨即選擇寄生鳥巢； ②在一組鳥巢中，質量最好的鳥巢將被保留到下一代； ③鳥巢的數量是固定的，且宿主發現外來鳥蛋的概率即發現概率Pa∈(0,1)。

(10)

Levy(β)～u=t-λ1<λ≤3

(11)

式中：u為隨機步長；λ為冪次系數。在布谷鳥算法實現過程中，Levy飛行的隨機步長由Mantegna算法計算得到

(12)

式中：β=3/2；μ和ν為服從正態分布的隨機數，即

(13)

(14)

式中Γ為標準Gamma函數。

布谷鳥搜索算法通過上述Levy飛行機制對鳥巢的位置更新，將更新后的新解與原來的解進行對比，選擇較好的結果并保留； 利用發現概率Pa對相對較差的解利用下式進一步更新

(15)

重復以上的過程，直到滿足設定誤差的最優解或者達到限制的最大迭代次數。圖1為三維空間Levy飛行示意圖[29]，展示了布谷鳥搜索算法在三維空間中隨機搜索步長。

圖1 三維空間Levy飛行示意圖P1、P2、P3分別表示三維坐標系的坐標軸

2 反演方法

2.1 目標函數構建

通過彈性阻抗與儲層物性參數之間的關系可構建反演目標函數[30]，具體形式為

(16)

式中：參數向量m=(φ,Vsh,Sw)是待反演的儲層物性參數，分別為孔隙度、泥質含量、含水飽和度； EIreal,i(m)為利用實際數據計算得到的彈性阻抗； EIcal,i(m)為利用式(9)計算得到的彈性阻抗；i=1,2,3表示三個不同角度。

2.2 基于布谷鳥算法的儲層物性參數反演

(2)對比初始鳥巢函數值，并記錄當前最優解Gbest、函數值F(G)。

圖2 基于布谷鳥算法的儲層物性參數反演流程圖

(4)產生一組在(0，1)均勻分布的隨機數r，并與發現概率Pa比較。若r>Pa，則利用式(15)對該位置進一步更新，反之不變。計算更新后鳥巢位置的函數值，與之前對應的鳥巢位置進行對比、選擇，進一步更新最優解Gbest、函數值F(G)。

(5)判斷搜索結果是否滿足本文算法的預設條件(Max_iter、Tol)，若沒有達到預設條件，則返回步驟(3)； 若滿足預設條件，停止迭代，輸出最優解Gbest。

3 模型試算

為檢驗本文基于布谷鳥搜索算法的儲層物性參數預測方法的效果，選用一維模型數據進行試算。

模型數據曲線如圖3所示。圖4是利用孔隙度、泥質含量、含水飽和度計算得到的三個角度(10°、20°和30°)的彈性阻抗與利用彈性參數表征的彈性阻抗方程計算得到的三個角度的真實彈性阻抗的對比。由圖可以看出，非線性條件下建立的統計巖石物理模型(式(9))計算得到的彈性阻抗曲線與利用彈性阻抗公式計算(式(5))得到的結果吻合程度較高，說明本文所建立的統計巖石物理模型是準確的、合理的。

利用本文全局尋優反演方法對目標函數進行求解，反演目標層段儲層物性參數。圖5為利用布谷鳥算法反演得到的儲層物性參數與原始曲線對比。由圖可知，利用本文方法反演得到的物性參數與原始數據整體趨勢基本吻合，僅在1760～1765m處結果變化幅度較大，但反演參數誤差在合理區間內，驗證了本文方法的可行性。

圖3 一維模型曲線

圖4 不同角度、不同模型真實(黑線)與計算(紅虛線)彈性阻抗對比不同角度線性關系下得到的確定性巖石物理模型：(a)10°； (b)20°； (c) 30°。不同角度線性關系下得到的統計性巖石物理模型：(d)10°； (e)20°； (f)30°。不同角度非線性關系下得到的統計性巖石物理模型：(g)10°； (h)20°； (i) 30°

圖5 不同方法物性參數反演結果對比

4 應用實例

選取D工區A井實際資料對本文方法進行測試。首先進行一維實際數據(圖6)測試。在反演儲層物性參數之前，需要分析實際測井數據，構建統計巖石物理模型。通過統計巖石物理模型由孔隙度、泥質含量、含水飽和度計算彈性阻抗曲線(10°、20°和30°)，并與實際彈性阻抗結果進行對比(圖7)。由圖可見，實際曲線與計算結果整體變化趨勢相同。

利用彈性阻抗數據反演儲層物性參數(孔隙度、泥質含量、含水飽和度)，結果如圖8所示。由圖可見，基于布谷鳥算法反演得到的儲層物性參數與測井曲線基本一致，整體變化趨勢相同。其中，孔隙度相關系數為75.76%，泥質含量相關系數為80.31%，含水飽和度相關系數為84.55%。

圖6 A井實際測井曲線

圖7 A井統計巖石物理模型與阻抗方程計算結果對比

圖8 A井基于布谷鳥算法的儲層物性參數反演結果

然后將本文方法拓展應用于二維測線，對圖9a進行儲層物性參數反演，結果如圖9b～圖9d所示。由圖可見，物性參數(孔隙度、泥質含量、含水飽和度)在橫向上具有較好的連續性，大套層組間的界面較為清晰； 并且在A、B兩井處，與測井巖性解釋結果及物性參數曲線吻合較好(黑色圓圈處)。

圖9 利用本文方法的不同物性參數反演結果(a)原始數據(以中角度彈性阻抗為例)； (b)孔隙度； (c)泥質含量； (d)含水飽和度

5 結束語

針對當前常規物性參數反演容易陷入局部極值的問題，本文引入元啟發式算法——布谷鳥算法，提出了一種基于布谷鳥算法的儲層物性參數反演方法。采用布谷鳥算法對目標函數求解，反演儲層物性參數，兼具大量局部搜索及必要的全局搜索。對比利用布谷鳥算法預測得到的儲層物性參數反演結果與實際測井物性參數可知，預測結果與實際測井數據具有較高的吻合度，誤差較小，證明了本文算法的可行性和有效性。

Levy分布的重尾分布特征保證了布谷鳥算法的全局尋優能力，但也同樣降低了該算法的計算效率。考慮全局尋優算法的計算效率下降問題，將來可進一步采取布谷鳥算法與常規全局尋優算法的結合策略，在求解精度允許范圍內提高算法的計算效率，進一步提高本文方法的實用性。