基于改進(jìn)模擬退火算法的橫波速度求取

王璞，吳國忱，李偉

（中國石油大學(xué)（華東） 地球科學(xué)與技術(shù)學(xué)院，山東 青島266580）

0 引言

模擬退火算法作為一種隨機性全局最優(yōu)求解方法，具有廣泛的適用性，很多學(xué)者對其進(jìn)行過研究［1］。模擬退火算法最早是由N.Metropolis 等［2］于1953年提出的，常規(guī)的模擬退火算法是指Kirkpatrick 等［3］提出的Metropolis 算法。非線性反演算法往往受到計算效率的制約，為了提高模擬退火的計算效率，Ingber［4］提出的非?？焖倌M退火算法（VFSA）中，擾動模型采用了依賴于溫度的似Cauchy 分布代替了常規(guī)模擬退火算法中的高斯分布；Press 等［5］采用單純形法與模擬退火算法相結(jié)合的算法；Penna［6］由廣義Gibbs 分布給出新的接受概率計算表達(dá)式；姚姚［7］從模擬退火最低溫度的選擇入手，根據(jù)模擬退火算法與統(tǒng)計力學(xué)的Gibbs-Markov 模型之間的關(guān)系導(dǎo)出臨界溫度的近似表達(dá)式；紀(jì)晨等［8］在模擬退火算法中引入均勻設(shè)計方法；文浩［9］采用遺傳模擬退火算法，對阿爾奇公式中的a，m，n 參數(shù)進(jìn)行估計，并通過現(xiàn)場應(yīng)用對測井?dāng)?shù)據(jù)進(jìn)行了模擬計算；閆安菊等［10］借鑒模擬退火的冷卻思想，設(shè)定高斯濾波尺度因子和梯度閥值為時間的函數(shù)，構(gòu)建了各向異性擴散濾波模型。模擬退火算法是建立在隨機搜索方法基礎(chǔ)之上的，因此，擾動模型需要滿足在解分布的最大范圍內(nèi)均勻隨機找值，達(dá)到對目標(biāo)函數(shù)的精度要求。

本文以橫波速度預(yù)測為例進(jìn)行分析，來解決實際工區(qū)中橫波信息缺失的現(xiàn)象。在求取橫波速度時，需要用到儲層孔隙度，然而實際資料中泥巖層段孔隙度是不作解釋的。這就需要基于非線性全局尋優(yōu)算法來重構(gòu)孔隙度，在反演孔隙度時，模擬退火算法是應(yīng)用較多的一種方法。本文基于DEM 模型［11］，在縱波速度約束下，用改進(jìn)的模擬退火算法對某工區(qū)1 口砂泥巖井進(jìn)行應(yīng)用分析，取得了很好的效果。

1 模擬退火算法基本原理

模擬退火算法的思想源于物理中固體物質(zhì)的退火過程與一般組合優(yōu)化問題之間的相似性，它把優(yōu)化問題的可行解看成是材料的各種狀態(tài)，將優(yōu)化目標(biāo)視為材料的能量或者熵。在優(yōu)化過程中，它不僅接受讓目標(biāo)函數(shù)變好的解，對于讓目標(biāo)函數(shù)變差的解也會以某一概率接受，這樣可使算法跳出局部最優(yōu)解，最終獲得全局最優(yōu)解。

模擬退火算法的一般步驟：1） 給定初始溫度T=T0，在解可能取到的范圍內(nèi)隨機產(chǎn)生初始狀態(tài)X=X0，并求取目標(biāo)函數(shù)值F（X0）。2）基于上一個狀態(tài)擾動產(chǎn)生一個新解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），即對解進(jìn)行一定次數(shù)的擾動。5）將溫度T 緩慢退溫。6）重復(fù)步驟2）至5），直至滿足某個收斂準(zhǔn)則為止。

2 模擬退火算法的改進(jìn)

在Ingber［4］給出的非?？焖倌M退火方法中，在模型擾動時，采用依賴于溫度的似Cauchy 分布來產(chǎn)生新的擾動解，在高溫時對應(yīng)大范圍搜索，低溫時對應(yīng)當(dāng)前解附近小范圍搜索，有助于加快收斂速度。擾動模型的表示形式為

其中

式中：X（i）為擾動產(chǎn)生的第i 個解；Xmax和Xmin分別為X 可取的最大值和最小值；ε 為［0，1］范圍內(nèi)均勻分布的隨機數(shù)。

在生成新擾動解的過程中，可將非均勻變異思想與溫度結(jié)合創(chuàng)建新的擾動模型。非均勻變異算子要求在算法初始階段能夠均勻地搜索整個解空間，在算法后期則集中于若干小范圍內(nèi)精確搜索。搜索中給出的擾動模型可表示為

其中

式中：T（i）為第i 次退溫的溫度；Tmax為溫度可取最大值；b 是控制搜索范圍在退溫時減小快慢的量。

由式（4）可知，a 既保證了最大范圍內(nèi)尋優(yōu)，又保證了擾動模型在溫度高時對應(yīng)大范圍搜索、 在溫度低時對應(yīng)小范圍搜索的特性。相比似Cauchy 分布產(chǎn)生的擾動解，該方法具有較好的靈活性。通過調(diào)節(jié)控制變量b，可以有效提高模擬退火尋優(yōu)速度。b 值的選擇不能太大也不能太?。篵 值太大，迭代次數(shù)會減少，容易陷入局部最小值；b 值太小，尋優(yōu)時間會增長。

以上2 種產(chǎn)生擾動新解的方法，有個共同的特點，即在高溫時對應(yīng)大范圍搜索，在低溫時對應(yīng)當(dāng)前解附近小范圍搜索。常規(guī)模擬退火尋優(yōu)過程（見圖1a，1b）為：在某一溫度時，模擬退火尋優(yōu)過程如圖1a所示，尋優(yōu)半徑為R1，X*為準(zhǔn)確解位置，通過尋優(yōu)，最優(yōu)解由X1變?yōu)閄2；尋優(yōu)范圍R1依賴于溫度T，即溫度T 不變時，尋優(yōu)范圍固定。退溫時，尋優(yōu)過程如圖1b所示，搜索范圍由R1變?yōu)镽2，然后循環(huán)執(zhí)行圖1a和圖1b所示流程，直至找到全局最優(yōu)解。其中，由1a 到1b 為常規(guī)模擬退火流程，由1a 到1c，再到1d 為改進(jìn)后的流程。

圖1 模擬退火流程

對于圖1a所示溫度不變的過程，經(jīng)過改進(jìn)，分為2 步完成。第1 步先找到最優(yōu)解X2；第2 步，設(shè)定一個較小的搜索范圍R3，即a 取0.01 等較小值（見圖1c）。為了不影響尋優(yōu)速度，在搜索范圍R3內(nèi)只需進(jìn)行較少次數(shù)的搜索。由于該搜索范圍很小，因此，如果R3范圍內(nèi)有更好的解X3，是極易找到的，反之，保持解X2不變。退溫時，為了防止出現(xiàn)局部最優(yōu)解，搜索范圍仍采取依賴于溫度T 的R2，最優(yōu)解X3作為下次循環(huán)的初值（見圖1d）。通過改進(jìn)，擾動產(chǎn)生新解的次數(shù)減少，同時也減少了對初值和擾動模型的依賴。

接受概率采用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)后的模擬退火算法擾動模型基于非均勻變異思想，大大加快了收斂速度，同時在退溫前增加最優(yōu)解附近小范圍內(nèi)尋優(yōu)可以減小對擾動模型的依賴，從而減少了擾動模型的隨機賦值。

3 數(shù)值模型分析

為了驗證該算法改進(jìn)后的尋優(yōu)效果，構(gòu)建了縱波約束下反演孔隙縱橫比的巖石物理模型，采用DEM 模型。基質(zhì)礦物為純砂巖，流體為油水混合，具體參數(shù)如表1所示。

表1 模型參數(shù)

為了方便說明，將上述分析的模擬退火算法進(jìn)行了分類表示。方法1 的擾動模型為依賴于溫度的似Cauchy 分布；方法2 的擾動模型基于非均勻變異思想；方法3 為改進(jìn)的方法1，即方法1 中退溫前采取兩步法；方法4 指改進(jìn)的方法2，即方法2 中退溫前采取兩步法。

圖2比較了初始孔隙縱橫比為0.4 時方法1 和方法2 的收斂速度?？梢钥闯?，方法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ù)測精度并沒有受到干擾，能較好地求出孔隙縱橫比。

圖2 孔隙縱橫比隨迭代次數(shù)變化

由圖3可以看出，初始孔隙縱橫比取0.4，b 取5時，4 種方法都能較準(zhǔn)確地找到精確解。在找到全局最優(yōu)解之前，方法3 和方法4 比方法1 和方法2 對解的擾動次數(shù)要少，即溫度不變時，通過增加最優(yōu)解附近小范圍內(nèi)尋優(yōu)，減少了對擾動模型的依賴。但是，由于增加了區(qū)域R3的尋優(yōu)，收斂速度沒有得到改善。

圖3 孔隙縱橫比隨擾動次數(shù)變化

4 井資料實例應(yīng)用

實際工區(qū)儲層解釋時，往往缺失橫波信息，因此采用模擬退火方法2 和方法4，對某地區(qū)1 口井進(jìn)行了橫波求取試算。為驗證方法的有效性，所取井段為有橫波信息的砂泥巖層，采用DEM 模型，用Gassmann 方程進(jìn)行流體替換，模型中模擬退火算法用來重構(gòu)孔隙度，即計算虛擬孔隙度，來解決實際生產(chǎn)中孔隙度解釋不準(zhǔn)確，或者泥巖段沒有孔隙度信息的情況。根據(jù)上述數(shù)值分析，2 種模擬退火方法都能較準(zhǔn)確地找到全局最優(yōu)解。這里運用2 種模擬退火方法，在縱波速度的約束下，由DEM 模型求取橫波速度（見圖4）。

圖4 橫波速度預(yù)測結(jié)果

由圖4a可以看出，通過調(diào)節(jié)孔隙度，計算縱波速度和實測縱波速度，誤差非常小。擾動模型都基于非均勻變異思想，不同的是圖4c考慮了退溫前小范圍內(nèi)尋優(yōu)?？梢钥闯?，2 種方法求取的橫波速度與實際測井橫波速度都吻合較好，誤差較小。圖4b證明了基于非均勻變異思想的模擬退火算法在實際應(yīng)用中的適用性，圖4c證明了同時考慮非均勻變異思想帶來的模型擾動和增加退溫前小范圍尋優(yōu)得到的模擬退火算法也具有很好的適用性，即驗證了改進(jìn)后模擬退火算法的有效性。

5 結(jié)束語

本文對常規(guī)模擬退火算法依賴于溫度的似Cauchy分布擾動模型進(jìn)行了改進(jìn)，改進(jìn)后的擾動模型基于非均勻變異思想，大大提高了收斂速度。另外，模擬退火退溫前，增加最優(yōu)解附近小范圍尋優(yōu)，可以減小對擾動模型的依賴，從而減少了擾動模型的隨機賦值。對改進(jìn)后的模擬退火算法進(jìn)行了數(shù)值實例分析和測井資料實例應(yīng)用，結(jié)果表明，改進(jìn)后模擬退火算法和常規(guī)模擬退火算法相比，在迭代次數(shù)和擾動次數(shù)上具有優(yōu)勢。通過測井實例應(yīng)用，驗證了該方法的實用性和有效性。

［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é)報，1995，38（5）：643-650.

［8］紀(jì)晨，姚振興.用于地震物理反演的均勻設(shè)計優(yōu)化算法［J］.地球物理學(xué)報，1996，39（2）：233-242.

［9］文浩.遺傳模擬退火算法在阿爾奇公式參數(shù)估計中的應(yīng)用［J］.斷塊油氣田，2008，15（1）：105-107.

［10］閆安菊，蔡涵鵬.改進(jìn)的各向異性擴散濾波方法壓制地震數(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.