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

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

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

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

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

模擬退火算法能夠較好地實(shí)現(xiàn)全局最優(yōu)化求解?；诜蔷鶆蜃儺愃枷?，文中給出了一個(gè)新的模擬退火算法擾動(dòng)模型，該模型具有一定的靈活性，通過(guò)選取合適的擾動(dòng)模型參數(shù)，可以達(dá)到減少迭代次數(shù)、提高收斂速度的目的。在模擬退火算法退溫前，增加最優(yōu)解附近小范圍內(nèi)尋優(yōu)，以減少對(duì)擾動(dòng)模型的依賴(lài)，通過(guò)改進(jìn)，模擬退火算法得到了優(yōu)化。數(shù)值實(shí)例分析表明，改進(jìn)后的模擬退火算法比常規(guī)模擬退火算法在迭代次數(shù)和擾動(dòng)次數(shù)上有一定的優(yōu)勢(shì)。采用該算法，利用DEM模型進(jìn)行橫波速度求取后，與測(cè)井橫波比較，誤差較小，從而驗(yàn)證了新模擬退火算法的實(shí)用性和有效性。

模擬退火；最優(yōu)解；迭代次數(shù)；橫波

0　引言

模擬退火算法作為一種隨機(jī)性全局最優(yōu)求解方法，具有廣泛的適用性，很多學(xué)者對(duì)其進(jìn)行過(guò)研究［1］。模擬退火算法最早是由N.Metropolis等［2］于1953年提出的，常規(guī)的模擬退火算法是指Kirkpatrick等［3］提出的Metropolis算法。非線性反演算法往往受到計(jì)算效率的制約，為了提高模擬退火的計(jì)算效率，Ingber［4］提出的非常快速模擬退火算法（VFSA）中，擾動(dòng)模型采用了依賴(lài)于溫度的似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ì)，并通過(guò)現(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)模型需要滿(mǎn)足在解分布的最大范圍內(nèi)均勻隨機(jī)找值，達(dá)到對(duì)目標(biāo)函數(shù)的精度要求。

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

1　模擬退火算法基本原理

模擬退火算法的思想源于物理中固體物質(zhì)的退火過(guò)程與一般組合優(yōu)化問(wèn)題之間的相似性，它把優(yōu)化問(wèn)題的可行解看成是材料的各種狀態(tài)，將優(yōu)化目標(biāo)視為材料的能量或者熵。在優(yōu)化過(guò)程中，它不僅接受讓目標(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），直至滿(mǎn)足某個(gè)收斂準(zhǔn)則為止。

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

在Ingber［4］給出的非?？焖倌M退火方法中，在模型擾動(dòng)時(shí)，采用依賴(lài)于溫度的似Cauchy分布來(lái)產(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)解的過(guò)程中，可將非均勻變異思想與溫度結(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)解，該方法具有較好的靈活性。通過(guò)調(diào)節(jié)控制變量b，可以有效提高模擬退火尋優(yōu)速度。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)過(guò)程（見(jiàn)圖1a，1b）為：在某一溫度時(shí)，模擬退火尋優(yōu)過(guò)程如圖1a所示，尋優(yōu)半徑為R1，X*為準(zhǔn)確解位置，通過(guò)尋優(yōu)，最優(yōu)解由X1變?yōu)閄2；尋優(yōu)范圍R1依賴(lài)于溫度T，即溫度T不變時(shí)，尋優(yōu)范圍固定。退溫時(shí)，尋優(yōu)過(guò)程如圖1b所示，搜索范圍由R1變?yōu)镽2，然后循環(huán)執(zhí)行圖1a和圖1b所示流程，直至找到全局最優(yōu)解。其中，由1a到1b為常規(guī)模擬退火流程，由1a到1c，再到1d為改進(jìn)后的流程。

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

和常規(guī)模擬退火算法比較，改進(jìn)后的模擬退火算法擾動(dòng)模型基于非均勻變異思想，大大加快了收斂速度，同時(shí)在退溫前增加最優(yōu)解附近小范圍內(nèi)尋優(yōu)可以減小對(duì)擾動(dòng)模型的依賴(lài)，從而減少了擾動(dòng)模型的隨機(jī)賦值。

3　數(shù)值模型分析

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

表1　模型參數(shù)

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

圖2比較了初始孔隙縱橫比為0.4時(shí)方法1和方法2的收斂速度?？梢钥闯?，方法2通過(guò)調(diào)節(jié)參數(shù)b，可以提高收斂速度。根據(jù)式（4），隨著b值的增大，在退溫過(guò)程中，尋優(yōu)范圍減小的速度加快，理論上減小了大范圍搜索全局最優(yōu)解的次數(shù)，因此，參數(shù)b不宜取值太大，以免出現(xiàn)局部最優(yōu)解。通過(guò)數(shù)值分析，可以看出b值取6左右，既可以有效地減少迭代次數(shù)，又可以避免迭代次數(shù)過(guò)少而產(chǎn)生局部最小值。精確b值的確定，還有待進(jìn)一步研究。由圖2還可以看出，方法2雖然減少了迭代次數(shù)，但是其預(yù)測(cè)精度并沒(méi)有受到干擾，能較好地求出孔隙縱橫比。

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

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

圖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)行流體替換，模型中模擬退火算法用來(lái)重構(gòu)孔隙度，即計(jì)算虛擬孔隙度，來(lái)解決實(shí)際生產(chǎn)中孔隙度解釋不準(zhǔn)確，或者泥巖段沒(méi)有孔隙度信息的情況。根據(jù)上述數(shù)值分析，2種模擬退火方法都能較準(zhǔn)確地找到全局最優(yōu)解。這里運(yùn)用2種模擬退火方法，在縱波速度的約束下，由DEM模型求取橫波速度（見(jiàn)圖4）。

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

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

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

本文對(duì)常規(guī)模擬退火算法依賴(lài)于溫度的似Cauchy分布擾動(dòng)模型進(jìn)行了改進(jìn)，改進(jìn)后的擾動(dòng)模型基于非均勻變異思想，大大提高了收斂速度。另外，模擬退火退溫前，增加最優(yōu)解附近小范圍尋優(yōu)，可以減小對(duì)擾動(dòng)模型的依賴(lài)，從而減少了擾動(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ì)。通過(guò)測(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.

（編輯石愛(ài)萍）

Shear-wave velocity prediction based on revised and simulated annealing algorithm

Wang Pu，Wu Guochen，Li Wei
（College of Geosciences and Technology，China University of Petroleum，Qingdao 266580，China）

Simulated annealing algorithm is a good way to get the optimal solution.Based on the idea of non-uniform mutation，a new simulated annealing perturbation model is given in this paper.The perturbation model has a certain flexibility in its application.By choosing an appropriate parameter，the iterative number is reduced and the convergence speed is improved.In addition，before the temperature is decreased in the simulated annealing algorithm，the dependence on perturbation model can be reduced by adding a small search scope near the optimal solution.The results show that the revised method is better than traditional method from iterative and perturbing number.A rock physics model is constructed with DEM model to estimate shear-wave velocity under constraint of P-wave velocity.Compared with logging data，the computed shear-wave velocities have high accuracy.

simulated annealing；optimal solution；iterative number；shear-wave

國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃（973計(jì)劃）項(xiàng)目“多元多尺度地球物理信息映射關(guān)系及反演”（2013CB228604）；國(guó)家科技重大專(zhuān)項(xiàng)專(zhuān)題“復(fù)雜儲(chǔ)層及流體識(shí)別方法與預(yù)測(cè)技術(shù)研究”（2011ZX05009-003-004）

TE319

A

10.6056/dkyqt201503013

2014-11-13；改回日期：2015-03-06。

王璞，男，1991年生，在讀碩士研究生，研究方向?yàn)榈卣饚r石物理。E-mail：826378939@qq.com。

引用格式：王璞，吳國(guó)忱，李偉.基于改進(jìn)模擬退火算法的橫波速度求?。跩］.斷塊油氣田，2015，22（3）：330-333. Wang Pu，Wu Guochen，Li Wei.Shear-wave velocity prediction based on revised and simulated annealing algorithm［J］.Fault-Block Oil&Gas Field，2015，22（3）：330-333.