基于布谷鳥(niǎo)算法的儲(chǔ)層物性參數(shù)同步反演

劉仕友 段治川 周 凡 汪 銳

(①中海石油(中國(guó))有限公司海南分公司，海南?？?570000； ②中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580)

0 引言

儲(chǔ)層物性參數(shù)(如孔隙度、泥質(zhì)含量等)可為儲(chǔ)層評(píng)價(jià)提供依據(jù)。隨著油氣勘探開(kāi)發(fā)目標(biāo)日趨復(fù)雜，對(duì)儲(chǔ)層物性參數(shù)的反演精度的要求越來(lái)越高。

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

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

將智能優(yōu)化算法引入地球物理問(wèn)題的求解時(shí)，部分智能優(yōu)化算法由于較強(qiáng)的局部搜索能力，增強(qiáng)了在某一范圍內(nèi)的搜索精度，但同時(shí)會(huì)導(dǎo)致算法的全局搜索能力不夠，在求解過(guò)程中易陷入局部最優(yōu)問(wèn)題中。

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

1 方法原理

1.1 巖石物理模型構(gòu)建

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

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

(1)

改寫(xiě)為矩陣形式

(2)

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

當(dāng)?shù)卣饛椥詤?shù)與物性參數(shù)的關(guān)系無(wú)法近似為線性時(shí)，可通過(guò)多項(xiàng)式擬合的方法，將其表示為多元高階方程的形式

(3)

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

(4)

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

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

(5)

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

對(duì)式(5)兩邊取對(duì)數(shù)，整理得

lnEI=AlnVP+BlnVS+Clnρ+D

(6)

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

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

(7)

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

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

(8)

引入隨機(jī)誤差項(xiàng)ε，則巖石物理模型可表示為

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

(9)

至此，建立了彈性阻抗與物性參數(shù)之間的函數(shù)關(guān)系，為后續(xù)儲(chǔ)層物性參數(shù)的反演奠定了基礎(chǔ)。

1.2 布谷鳥(niǎo)搜索算法

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

在布谷鳥(niǎo)搜索算法中，設(shè)定了三個(gè)理想狀態(tài)：①每只布谷鳥(niǎo)一次只產(chǎn)一只蛋，并隨即選擇寄生鳥(niǎo)巢； ②在一組鳥(niǎo)巢中，質(zhì)量最好的鳥(niǎo)巢將被保留到下一代； ③鳥(niǎo)巢的數(shù)量是固定的，且宿主發(fā)現(xiàn)外來(lái)鳥(niǎo)蛋的概率即發(fā)現(xiàn)概率Pa∈(0,1)。

(10)

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

(11)

式中：u為隨機(jī)步長(zhǎng)；λ為冪次系數(shù)。在布谷鳥(niǎo)算法實(shí)現(xiàn)過(guò)程中，Levy飛行的隨機(jī)步長(zhǎng)由Mantegna算法計(jì)算得到

(12)

式中：β=3/2；μ和ν為服從正態(tài)分布的隨機(jī)數(shù)，即

(13)

(14)

式中Γ為標(biāo)準(zhǔn)Gamma函數(shù)。

布谷鳥(niǎo)搜索算法通過(guò)上述Levy飛行機(jī)制對(duì)鳥(niǎo)巢的位置更新，將更新后的新解與原來(lái)的解進(jìn)行對(duì)比，選擇較好的結(jié)果并保留； 利用發(fā)現(xiàn)概率Pa對(duì)相對(duì)較差的解利用下式進(jìn)一步更新

(15)

重復(fù)以上的過(guò)程，直到滿足設(shè)定誤差的最優(yōu)解或者達(dá)到限制的最大迭代次數(shù)。圖1為三維空間Levy飛行示意圖[29]，展示了布谷鳥(niǎo)搜索算法在三維空間中隨機(jī)搜索步長(zhǎng)。

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

2 反演方法

2.1 目標(biāo)函數(shù)構(gòu)建

通過(guò)彈性阻抗與儲(chǔ)層物性參數(shù)之間的關(guān)系可構(gòu)建反演目標(biāo)函數(shù)[30]，具體形式為

(16)

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

2.2 基于布谷鳥(niǎo)算法的儲(chǔ)層物性參數(shù)反演

(2)對(duì)比初始鳥(niǎo)巢函數(shù)值，并記錄當(dāng)前最優(yōu)解Gbest、函數(shù)值F(G)。

圖2 基于布谷鳥(niǎo)算法的儲(chǔ)層物性參數(shù)反演流程圖

(4)產(chǎn)生一組在(0，1)均勻分布的隨機(jī)數(shù)r，并與發(fā)現(xiàn)概率Pa比較。若r>Pa，則利用式(15)對(duì)該位置進(jìn)一步更新，反之不變。計(jì)算更新后鳥(niǎo)巢位置的函數(shù)值，與之前對(duì)應(yīng)的鳥(niǎo)巢位置進(jìn)行對(duì)比、選擇，進(jìn)一步更新最優(yōu)解Gbest、函數(shù)值F(G)。

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

3 模型試算

為檢驗(yàn)本文基于布谷鳥(niǎo)搜索算法的儲(chǔ)層物性參數(shù)預(yù)測(cè)方法的效果，選用一維模型數(shù)據(jù)進(jìn)行試算。

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

利用本文全局尋優(yōu)反演方法對(duì)目標(biāo)函數(shù)進(jìn)行求解，反演目標(biāo)層段儲(chǔ)層物性參數(shù)。圖5為利用布谷鳥(niǎo)算法反演得到的儲(chǔ)層物性參數(shù)與原始曲線對(duì)比。由圖可知，利用本文方法反演得到的物性參數(shù)與原始數(shù)據(jù)整體趨勢(shì)基本吻合，僅在1760～1765m處結(jié)果變化幅度較大，但反演參數(shù)誤差在合理區(qū)間內(nèi)，驗(yàn)證了本文方法的可行性。

圖3 一維模型曲線

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

圖5 不同方法物性參數(shù)反演結(jié)果對(duì)比

4 應(yīng)用實(shí)例

選取D工區(qū)A井實(shí)際資料對(duì)本文方法進(jìn)行測(cè)試。首先進(jìn)行一維實(shí)際數(shù)據(jù)(圖6)測(cè)試。在反演儲(chǔ)層物性參數(shù)之前，需要分析實(shí)際測(cè)井?dāng)?shù)據(jù)，構(gòu)建統(tǒng)計(jì)巖石物理模型。通過(guò)統(tǒng)計(jì)巖石物理模型由孔隙度、泥質(zhì)含量、含水飽和度計(jì)算彈性阻抗曲線(10°、20°和30°)，并與實(shí)際彈性阻抗結(jié)果進(jìn)行對(duì)比(圖7)。由圖可見(jiàn)，實(shí)際曲線與計(jì)算結(jié)果整體變化趨勢(shì)相同。

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

圖6 A井實(shí)際測(cè)井曲線

圖7 A井統(tǒng)計(jì)巖石物理模型與阻抗方程計(jì)算結(jié)果對(duì)比

圖8 A井基于布谷鳥(niǎo)算法的儲(chǔ)層物性參數(shù)反演結(jié)果

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

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

5 結(jié)束語(yǔ)

針對(duì)當(dāng)前常規(guī)物性參數(shù)反演容易陷入局部極值的問(wèn)題，本文引入元啟發(fā)式算法——布谷鳥(niǎo)算法，提出了一種基于布谷鳥(niǎo)算法的儲(chǔ)層物性參數(shù)反演方法。采用布谷鳥(niǎo)算法對(duì)目標(biāo)函數(shù)求解，反演儲(chǔ)層物性參數(shù)，兼具大量局部搜索及必要的全局搜索。對(duì)比利用布谷鳥(niǎo)算法預(yù)測(cè)得到的儲(chǔ)層物性參數(shù)反演結(jié)果與實(shí)際測(cè)井物性參數(shù)可知，預(yù)測(cè)結(jié)果與實(shí)際測(cè)井?dāng)?shù)據(jù)具有較高的吻合度，誤差較小，證明了本文算法的可行性和有效性。

Levy分布的重尾分布特征保證了布谷鳥(niǎo)算法的全局尋優(yōu)能力，但也同樣降低了該算法的計(jì)算效率?？紤]全局尋優(yōu)算法的計(jì)算效率下降問(wèn)題，將來(lái)可進(jìn)一步采取布谷鳥(niǎo)算法與常規(guī)全局尋優(yōu)算法的結(jié)合策略，在求解精度允許范圍內(nèi)提高算法的計(jì)算效率，進(jìn)一步提高本文方法的實(shí)用性。