王璞，吳國(guó)忱，李偉

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

0 引言

模擬退火算法作為一種隨機(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)用分析，取得了很好的效果。

1 模擬退火算法基本原理

模擬退火算法的思想源于物理中固體物質(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)則為止。

2 模擬退火算法的改進(jì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ī)賦值。

3 數(shù)值模型分析

為了驗(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ù)變化

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

實(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)后模擬退火算法的有效性。

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

本文對(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.