王璞,吳國(guó)忱,李偉
(中國(guó)石油大學(xué)(華東) 地球科學(xué)與技術(shù)學(xué)院,山東 青島266580)
模擬退火算法作為一種隨機(jī)性全局最優(yōu)求解方法,具有廣泛的適用性,很多學(xué)者對(duì)其進(jìn)行過研究[1]。模擬退火算法最早是由N.Metropolis 等[2]于1953年提出的,常規(guī)的模擬退火算法是指Kirkpatrick 等[3]提出的Metropolis 算法。非線性反演算法往往受到計(jì)算效率的制約,為了提高模擬退火的計(jì)算效率,Ingber[4]提出的非常快速模擬退火算法(VFSA)中,擾動(dòng)模型采用了依賴于溫度的似Cauchy 分布代替了常規(guī)模擬退火算法中的高斯分布;Press 等[5]采用單純形法與模擬退火算法相結(jié)合的算法;Penna[6]由廣義Gibbs 分布給出新的接受概率計(jì)算表達(dá)式;姚姚[7]從模擬退火最低溫度的選擇入手,根據(jù)模擬退火算法與統(tǒng)計(jì)力學(xué)的Gibbs-Markov 模型之間的關(guān)系導(dǎo)出臨界溫度的近似表達(dá)式;紀(jì)晨等[8]在模擬退火算法中引入均勻設(shè)計(jì)方法;文浩[9]采用遺傳模擬退火算法,對(duì)阿爾奇公式中的a,m,n 參數(shù)進(jìn)行估計(jì),并通過現(xiàn)場(chǎng)應(yīng)用對(duì)測(cè)井?dāng)?shù)據(jù)進(jìn)行了模擬計(jì)算;閆安菊等[10]借鑒模擬退火的冷卻思想,設(shè)定高斯濾波尺度因子和梯度閥值為時(shí)間的函數(shù),構(gòu)建了各向異性擴(kuò)散濾波模型。模擬退火算法是建立在隨機(jī)搜索方法基礎(chǔ)之上的,因此,擾動(dòng)模型需要滿足在解分布的最大范圍內(nèi)均勻隨機(jī)找值,達(dá)到對(duì)目標(biāo)函數(shù)的精度要求。
本文以橫波速度預(yù)測(cè)為例進(jìn)行分析,來解決實(shí)際工區(qū)中橫波信息缺失的現(xiàn)象。在求取橫波速度時(shí),需要用到儲(chǔ)層孔隙度,然而實(shí)際資料中泥巖層段孔隙度是不作解釋的。這就需要基于非線性全局尋優(yōu)算法來重構(gòu)孔隙度,在反演孔隙度時(shí),模擬退火算法是應(yīng)用較多的一種方法。本文基于DEM 模型[11],在縱波速度約束下,用改進(jìn)的模擬退火算法對(duì)某工區(qū)1 口砂泥巖井進(jìn)行應(yīng)用分析,取得了很好的效果。
模擬退火算法的思想源于物理中固體物質(zhì)的退火過程與一般組合優(yōu)化問題之間的相似性,它把優(yōu)化問題的可行解看成是材料的各種狀態(tài),將優(yōu)化目標(biāo)視為材料的能量或者熵。在優(yōu)化過程中,它不僅接受讓目標(biāo)函數(shù)變好的解,對(duì)于讓目標(biāo)函數(shù)變差的解也會(huì)以某一概率接受,這樣可使算法跳出局部最優(yōu)解,最終獲得全局最優(yōu)解。
模擬退火算法的一般步驟:1) 給定初始溫度T=T0,在解可能取到的范圍內(nèi)隨機(jī)產(chǎn)生初始狀態(tài)X=X0,并求取目標(biāo)函數(shù)值F(X0)。2)基于上一個(gè)狀態(tài)擾動(dòng)產(chǎn)生一個(gè)新解X1,求取目標(biāo)函數(shù)F(X1),進(jìn)而得到ΔF=F(X0)-F(X1)。3)若ΔF<0,則接受新解X1;若ΔF>0,則新解以概率P=exp(-ΔF/T)進(jìn)行接受。4)在同一溫度T0下,重復(fù)步驟2)和3),即對(duì)解進(jìn)行一定次數(shù)的擾動(dòng)。5)將溫度T 緩慢退溫。6)重復(fù)步驟2)至5),直至滿足某個(gè)收斂準(zhǔn)則為止。
在Ingber[4]給出的非??焖倌M退火方法中,在模型擾動(dòng)時(shí),采用依賴于溫度的似Cauchy 分布來產(chǎn)生新的擾動(dòng)解,在高溫時(shí)對(duì)應(yīng)大范圍搜索,低溫時(shí)對(duì)應(yīng)當(dāng)前解附近小范圍搜索,有助于加快收斂速度。擾動(dòng)模型的表示形式為
其中
式中:X(i)為擾動(dòng)產(chǎn)生的第i 個(gè)解;Xmax和Xmin分別為X 可取的最大值和最小值;ε 為[0,1]范圍內(nèi)均勻分布的隨機(jī)數(shù)。
在生成新擾動(dòng)解的過程中,可將非均勻變異思想與溫度結(jié)合創(chuàng)建新的擾動(dòng)模型。非均勻變異算子要求在算法初始階段能夠均勻地搜索整個(gè)解空間,在算法后期則集中于若干小范圍內(nèi)精確搜索。搜索中給出的擾動(dòng)模型可表示為
其中
式中:T(i)為第i 次退溫的溫度;Tmax為溫度可取最大值;b 是控制搜索范圍在退溫時(shí)減小快慢的量。
由式(4)可知,a 既保證了最大范圍內(nèi)尋優(yōu),又保證了擾動(dòng)模型在溫度高時(shí)對(duì)應(yīng)大范圍搜索、 在溫度低時(shí)對(duì)應(yīng)小范圍搜索的特性。相比似Cauchy 分布產(chǎn)生的擾動(dòng)解,該方法具有較好的靈活性。通過調(diào)節(jié)控制變量b,可以有效提高模擬退火尋優(yōu)速度。b 值的選擇不能太大也不能太小:b 值太大,迭代次數(shù)會(huì)減少,容易陷入局部最小值;b 值太小,尋優(yōu)時(shí)間會(huì)增長(zhǎng)。
以上2 種產(chǎn)生擾動(dòng)新解的方法,有個(gè)共同的特點(diǎn),即在高溫時(shí)對(duì)應(yīng)大范圍搜索,在低溫時(shí)對(duì)應(yīng)當(dāng)前解附近小范圍搜索。常規(guī)模擬退火尋優(yōu)過程(見圖1a,1b)為:在某一溫度時(shí),模擬退火尋優(yōu)過程如圖1a所示,尋優(yōu)半徑為R1,X*為準(zhǔn)確解位置,通過尋優(yōu),最優(yōu)解由X1變?yōu)閄2;尋優(yōu)范圍R1依賴于溫度T,即溫度T 不變時(shí),尋優(yōu)范圍固定。退溫時(shí),尋優(yōu)過程如圖1b所示,搜索范圍由R1變?yōu)镽2,然后循環(huán)執(zhí)行圖1a和圖1b所示流程,直至找到全局最優(yōu)解。其中,由1a 到1b 為常規(guī)模擬退火流程,由1a 到1c,再到1d 為改進(jìn)后的流程。
圖1 模擬退火流程
對(duì)于圖1a所示溫度不變的過程,經(jīng)過改進(jìn),分為2 步完成。第1 步先找到最優(yōu)解X2;第2 步,設(shè)定一個(gè)較小的搜索范圍R3,即a 取0.01 等較小值(見圖1c)。為了不影響尋優(yōu)速度,在搜索范圍R3內(nèi)只需進(jìn)行較少次數(shù)的搜索。由于該搜索范圍很小,因此,如果R3范圍內(nèi)有更好的解X3,是極易找到的,反之,保持解X2不變。退溫時(shí),為了防止出現(xiàn)局部最優(yōu)解,搜索范圍仍采取依賴于溫度T 的R2,最優(yōu)解X3作為下次循環(huán)的初值(見圖1d)。通過改進(jìn),擾動(dòng)產(chǎn)生新解的次數(shù)減少,同時(shí)也減少了對(duì)初值和擾動(dòng)模型的依賴。
接受概率采用P=exp(-ΔF/T)[1],降溫方式采用T(n)=T0αnm。其中:n 為迭代次數(shù),T0為初始值,α 取值為[0,1],m 取0.5 或1[1]。為了減少循環(huán)次數(shù),將每次循環(huán)尋優(yōu)得到的最優(yōu)解保存,并將它們進(jìn)行比較;如果一定次數(shù)內(nèi)沒有得到改善,循環(huán)終止。
和常規(guī)模擬退火算法比較,改進(jìn)后的模擬退火算法擾動(dòng)模型基于非均勻變異思想,大大加快了收斂速度,同時(shí)在退溫前增加最優(yōu)解附近小范圍內(nèi)尋優(yōu)可以減小對(duì)擾動(dòng)模型的依賴,從而減少了擾動(dòng)模型的隨機(jī)賦值。
為了驗(yàn)證該算法改進(jìn)后的尋優(yōu)效果,構(gòu)建了縱波約束下反演孔隙縱橫比的巖石物理模型,采用DEM 模型?;|(zhì)礦物為純砂巖,流體為油水混合,具體參數(shù)如表1所示。
表1 模型參數(shù)
為了方便說明,將上述分析的模擬退火算法進(jìn)行了分類表示。方法1 的擾動(dòng)模型為依賴于溫度的似Cauchy 分布;方法2 的擾動(dòng)模型基于非均勻變異思想;方法3 為改進(jìn)的方法1,即方法1 中退溫前采取兩步法;方法4 指改進(jìn)的方法2,即方法2 中退溫前采取兩步法。
圖2比較了初始孔隙縱橫比為0.4 時(shí)方法1 和方法2 的收斂速度??梢钥闯觯椒? 通過調(diào)節(jié)參數(shù)b,可以提高收斂速度。根據(jù)式(4),隨著b 值的增大,在退溫過程中,尋優(yōu)范圍減小的速度加快,理論上減小了大范圍搜索全局最優(yōu)解的次數(shù),因此,參數(shù)b 不宜取值太大,以免出現(xiàn)局部最優(yōu)解。通過數(shù)值分析,可以看出b值取6 左右,既可以有效地減少迭代次數(shù),又可以避免迭代次數(shù)過少而產(chǎn)生局部最小值。精確b 值的確定,還有待進(jìn)一步研究。由圖2還可以看出,方法2 雖然減少了迭代次數(shù),但是其預(yù)測(cè)精度并沒有受到干擾,能較好地求出孔隙縱橫比。
圖2 孔隙縱橫比隨迭代次數(shù)變化
由圖3可以看出,初始孔隙縱橫比取0.4,b 取5時(shí),4 種方法都能較準(zhǔn)確地找到精確解。在找到全局最優(yōu)解之前,方法3 和方法4 比方法1 和方法2 對(duì)解的擾動(dòng)次數(shù)要少,即溫度不變時(shí),通過增加最優(yōu)解附近小范圍內(nèi)尋優(yōu),減少了對(duì)擾動(dòng)模型的依賴。但是,由于增加了區(qū)域R3的尋優(yōu),收斂速度沒有得到改善。
圖3 孔隙縱橫比隨擾動(dòng)次數(shù)變化
實(shí)際工區(qū)儲(chǔ)層解釋時(shí),往往缺失橫波信息,因此采用模擬退火方法2 和方法4,對(duì)某地區(qū)1 口井進(jìn)行了橫波求取試算。為驗(yàn)證方法的有效性,所取井段為有橫波信息的砂泥巖層,采用DEM 模型,用Gassmann 方程進(jìn)行流體替換,模型中模擬退火算法用來重構(gòu)孔隙度,即計(jì)算虛擬孔隙度,來解決實(shí)際生產(chǎn)中孔隙度解釋不準(zhǔn)確,或者泥巖段沒有孔隙度信息的情況。根據(jù)上述數(shù)值分析,2 種模擬退火方法都能較準(zhǔn)確地找到全局最優(yōu)解。這里運(yùn)用2 種模擬退火方法,在縱波速度的約束下,由DEM 模型求取橫波速度(見圖4)。
圖4 橫波速度預(yù)測(cè)結(jié)果
由圖4a可以看出,通過調(diào)節(jié)孔隙度,計(jì)算縱波速度和實(shí)測(cè)縱波速度,誤差非常小。擾動(dòng)模型都基于非均勻變異思想,不同的是圖4c考慮了退溫前小范圍內(nèi)尋優(yōu)??梢钥闯?,2 種方法求取的橫波速度與實(shí)際測(cè)井橫波速度都吻合較好,誤差較小。圖4b證明了基于非均勻變異思想的模擬退火算法在實(shí)際應(yīng)用中的適用性,圖4c證明了同時(shí)考慮非均勻變異思想帶來的模型擾動(dòng)和增加退溫前小范圍尋優(yōu)得到的模擬退火算法也具有很好的適用性,即驗(yàn)證了改進(jìn)后模擬退火算法的有效性。
本文對(duì)常規(guī)模擬退火算法依賴于溫度的似Cauchy分布擾動(dòng)模型進(jìn)行了改進(jìn),改進(jìn)后的擾動(dòng)模型基于非均勻變異思想,大大提高了收斂速度。另外,模擬退火退溫前,增加最優(yōu)解附近小范圍尋優(yōu),可以減小對(duì)擾動(dòng)模型的依賴,從而減少了擾動(dòng)模型的隨機(jī)賦值。對(duì)改進(jìn)后的模擬退火算法進(jìn)行了數(shù)值實(shí)例分析和測(cè)井資料實(shí)例應(yīng)用,結(jié)果表明,改進(jìn)后模擬退火算法和常規(guī)模擬退火算法相比,在迭代次數(shù)和擾動(dòng)次數(shù)上具有優(yōu)勢(shì)。通過測(cè)井實(shí)例應(yīng)用,驗(yàn)證了該方法的實(shí)用性和有效性。
[1]張霖斌,姚振興.快速模擬退火算法及應(yīng)用[J].石油地球物理勘探,1997,32(5):654-660.
[2]Metropolis N,Rosenbluth A W,Rosenbluth M N,et al.Equation of state calculations by fast computing machines[J].The journal of chemical physics,1953,21(6):1087-1092.
[3]Kirkpatrick S.Optimization by simulated annealing: Quantitative studies[J].Journal of statistical physics,1984,34(5/6):975-986.
[4]Ingber L.Very fast simulated re-annealing[J].Mathematical and computer modelling,1989,12(8):967-973.
[5]Press W H,Teukolsky S A.Simulated annealing optimization over continuous spaces[J].Computers in Physics,1991,5(4):426-429.
[6]Penna T J P.Traveling salesman problem and Tsallis statistics[J].Physical Review E,1995,51(1):R1-R3.
[7]姚姚.地球物理非線性反演模擬退火法的改進(jìn)[J].地球物理學(xué)報(bào),1995,38(5):643-650.
[8]紀(jì)晨,姚振興.用于地震物理反演的均勻設(shè)計(jì)優(yōu)化算法[J].地球物理學(xué)報(bào),1996,39(2):233-242.
[9]文浩.遺傳模擬退火算法在阿爾奇公式參數(shù)估計(jì)中的應(yīng)用[J].斷塊油氣田,2008,15(1):105-107.
[10]閆安菊,蔡涵鵬.改進(jìn)的各向異性擴(kuò)散濾波方法壓制地震數(shù)據(jù)噪聲[J].斷塊油氣田,2012,19(5):592-595.
[11]Berryman J G.Single-scattering approximations for coefficients in Biot′s equations of poroelasticity [J].The Journal of the Acoustical Society of America,1992,91(2):551-571.