模拟退火法在地球物理反演中的应用-张明浩

模拟退火法在地球物理反演中的应用简述王占文（中国矿业大学北京，地球科学与测绘工程学院）摘要：模拟退火法源于统计热力物理学，它模拟熔融状态下物体缓慢冷却达到结晶状态的物理过程。模拟退火反演算法的基本思想是：生成一系列参数向量模拟粒子的热运动，通过缓慢地减小一个模拟温度的控制参数，使模拟的系统最终冷却结晶达到系统能量最小值的过程。模拟退火反演算法实质是利用了地球物理反演问题求解过程与熔化固体退火过程的相似性，开辟了地球物理反演的新途径，是非线性反演算法中一种最常用的算法。关键词：模拟退火法；地球物理反演；蒙特卡洛法；一、前言地球物理学中,模拟退火法的使用与高度多解性的资料解释有关。模拟退火不仅能提供很好的拟合模型,而且也能提供有地质意义的额外信息。例如:可以解释资料的那些模型能有什么样的差异？那些可用的模型是什么样子的？它们的相对概率是多少？模拟退火法（SimulatedAnnealingMethod，简称SA），源于统计热力物理学，是模拟物质退火的物理过程，即统计试验的组合优化过程。模拟退火算法的基本思想是:生成一系列参数向量模拟粒子的热运动，通过缓慢地减小一个模拟温度的控制参数，使模拟的热系统最终冷却结晶达到系统能量最小值。模拟退火算法的思想最早由Metropolis等人于1953年提出。Kirkpatrik等人研究证实，冷却后物质变成晶体时的状态，恰好是系统处于整体能量极小的状态，并于1983年成功把它应用于组合优化问题。基于这种类比，Rothman在1985年将退火原理引入到地震资料处理中的剩余静校正问题中，并称之为模拟退火法。紧接着Sen和Stoffa等将其用于一维地震波形反演、一维电测深反演。随后，模拟退火法在大地电磁测深、重磁资料等方法中。模拟退火反演算法与传统的线性反演方法相比该方法具有：不依赖初始模型的选择、能寻找全局最小点而不陷入局部极小、在反演过程中不用计算雅克比偏导数矩阵等优点。所以近年来，模拟退火法在地球物理资料的反演中得到广泛应用与研究，并且取得了一些令人瞩目的研究成果，它已经成为了目前非线性反演的重要方法之一。本文主要讨论了自模拟退火法提出以来的一些研究进展和成果，以及现阶段模拟退火法在地球物理反演中的应用情况，所参考文献有26篇，起止年份自1953年至2012年。二、模拟退火法的原理模拟退火算法(SimulatedAnnealing)源于统计物理学，据统计热力学，物体中的每个分子的状态服从Gibbs分布，即：式中:E(ri)为第i个分子的能量函数；；ri为第i个分子所处的状态；k为波耳兹曼常数，T表示温度；(ri)为第i个分子的概率密度，为了方便起见令k=1。算法源于对实际固体退火过程的模拟，即先将固体加温至充分高，再逐渐冷却。加温时，固体内部粒子变为无序状态，内能增大;而逐渐降温时，粒子趋于有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。因此，算法实际上是将优化问题类比为退火过程中能量的最低状态，也就是温度达到最低点时，概率分布中具有最大概率的状态。模拟退火算法的具体步骤如下:1)给定模型每一个参数变化范围,在这个范围内随机选择一个初始模型m0，并计算相应的目标函数值E(m0)；2)对当前模型进行扰动产生一个新模型m,计算相应的目标函数值E(m),得到E=E(m)-E(m0)；3)若E0,则新模型m被接受;若E0,则新模型m按概率P=exp(-E/T)进行接受,T为温度。当模型被接受时,置m0=m；4)在温度T下,重复一定次数的扰动和接受过程,即重复步骤2)、3)；5)缓慢降低温度T；6)重复步骤2)、5),直至收敛条件满足为止。为了避免最好的解在优化过程中被忽略掉,可以稍做改进,即在整个搜索过程中随时记下最好的值,因为该法的特点决定了最后的优化值并不一定就是最优的值。三、模拟退火法的发展和应用在模拟退火法的应用方面，上世纪90年代出现了一批新技术新应用，但是均处于起步阶段，仍有很多地方并不成熟和完善。和一般地球物理反演问题一样，Rothman首先将静校正问题转换为目标函数的极小化问题，即:minE(x)式中:E是共中心点道集的叠加能谱，x为任意炮点到接收点的静校量所组成的矢量。应该看到，由于炮点和接收点非常多，求解一个全局最优的静校正矢量并不是一件非常容易的事情。Rothman引入了Kirkpatrick等提出的模拟退火法对上述问题进行求解。Sen和Stoffa尝试将模拟退火方法用于一维地震波形反演。一维地震波形反演问题是要求一个最优的模型，使该模型正演的地震波形与实际地震波形的拟合度最大，即求这样一个模型矢量m，使观测波形和理论波形的拟合度E最大。为了加快对模型空间的搜索效率，Sen和Stoffa采用了热浴法(Heat-bath)，即在某个温度上，只变化一个参数而固定其它参数，使该模型参数从最小值依次搜索到最大值，计算不同向量对应的Gibbs条件概率根据该概率值对该模型参数进行寻优。然后变化另一个参数而固定其它参数，使另一个参数也从最小值依次搜索到最大值，根据上式计算Gibbs概率进行寻优，依次类推。整个参数寻优一次后，将温度降温，然后再重复刚才模型寻优的过程。随着温度的降低，反演参数将逐渐向最优值演化。PauloM.Carvalho和ArthurB.Weglein在1994年提出可以使用模拟退火算法为界面多次波衰减提供子波估计，且可提供良好的结果。该方法对于求取资料的比例因子，甚至当它变化很大时特别有效。这里，它将涉及早期在这种情况下试图使用梯度最佳化法的失败。应用与模拟退火法有关的梯度尚需进一步研究，有可能改善最终估计。在叠后层速度计算上，也有学者提出了可以利用模拟退火法，并且在早期的PE机上得到了实现。但是这种方法所得到的反演解的分辨率不高，直接解决地质问题是不现实的。不过模拟退火法得到的反演解是正确的，且比初始地质模型更加真实细致的反应底下地质的平均效应。所以可将模拟退火反演解作为共扼梯度法的初始解，获得了可提供更丰富信息并能解决地质问题的最终解。虽然在地球物理非线性反演中，模拟退火法是一种较为先进的启发式蒙特卡洛法，但是实际地质资料处理中，该方法存在着计算效率不高，甚至有时会失效的弊端。因此，诸多学者提出了对模拟退火法的改进。Rothman，Sen和Stoffa提出的用试验方法估求Tmin。需对系统作若干次迭代运算，然后才能得到合适的Tmin的估计值。姚姚提出由部分局部势能估求Tc*（及Tmin）的方法不仅避免了盲目性，计算效率也大为提高。因为计算局部势能只是系统能量计算中的很少一部分，即使全部局部势能计算出来也只相当于迭代运算中一次系统能量的计算。由于Tc与目标函数超曲面的整体形状密切相关，故降低Tc的最佳方法是利用先验信息对目标函数加以改进。于是提出一种根据模糊数学的思想利用先验信息改进目标函数以降低Tc，提高反演效率和稳定性的模糊化方法。另外有学者讨论了模拟退火算法的基础上，将已知条件转化成为具体的约束，以实现对反演过程的控制。所以，在用模拟退火方法解决全局寻优问题的同时，又能合理利用约束条件提高反演的精度和收敛速度。约束主要体现在两个方面：首先是反演中参数取值范围的确定，可参考测井等资料形成这种约束;其次是利用测井和地震解释资料形成合理的初始模型，这是实现约束的重要内容。为了验证法的效果，进行了理论模型试算，并对试算结果进行了分析。结果表明，该方法用于地震反演是有效的，能适应复杂模型的反演计算。20世纪90年代后期至新世纪，模拟退火法不断得到改进和完善，并取得了一些新进展。郑鸿明、彭立等（2001）提出实际资料处理中并不需要从某一高温态开始进行，因为输入剖面已接近一个局部极小值，也就是剖面上部分共中心点道集已经得到较好的叠加，对应这些共中心点道集所包含的炮点和检波点叠加能量或相关函数曲线，将有一个突出的极大值，这个极大值相对来说是比较可靠的。可以拟定一个标准，若极大值与次大值之比大于一个定值，则认为它就是我们所要获得的延迟所对应的峰值，这样考虑的优点就在于：相应于在迭代时保持这些静校正值为常数，它影响邻近区域的扰动向它所在的极小方向搜寻；对于这些炮点和检波点可采用快速降温方式，使它们的概率分布曲线很快脉冲化，降低迭代次数或增大迭代步长，即减少计算量，节省计算时间。而基于前人所做的工作，现阶段多参数优化的模拟退火法已经应用到波阻抗反演、地层横波各向异性反演中。而在一些复杂地表结构地区，地震资料信噪比低，经常规静校正后仍然存在很大的剩余静校正量，直接影响叠加资料的品质。相对于反射波剩余静校正，应用快速模拟退火静校正后，低信噪比区域地震资料同相轴的连续性和分辨率也有了极大的提高，叠加剖面质量改善明显。四、模拟退火法在地球物理反演中的前景改进的模拟退火算法相对于常规的模拟退火算法,加强了局部搜索功能,更有利于提高算法的收敛速度,数值优化结果证实改进算法的有效性和高效性。在前文中曾提到模拟退火法本质上是一种启发式的蒙特卡洛法，所以反演须在全空间进行搜索，要求得全局极小。理论和实践表明，模拟退火法比蒙特卡洛法计算效率高，但模拟退火反演法仍必须进行大量的正演模拟和反演计算，这就大大地增加了计算时间和本，使它在实际应用中受到了一定的限制，将模拟退火算法与其它线性化反演算法甚至是其它非线性反演算法进行结合，进行混合优化反演，是地球物理反演法值得关注的研究方向。参考文献[1]MetropolisN，RosenbluthA，RosenbluthM，et.al.Equationofstatecalculationsbyfastcomputingmachines[J]，J.Chem.Phys.，1953，21:1087~1092.[2]KirkpatrickS，CelattCD，VecchiMP.Optimizationbysimulatedannealing[J]，Sciences，1983，220:671~680.[3]RothmanDH.Nonlinearinversion，statisticalmechanics，andresidualstaticsestimation[J]，Geophysics，1985，50(12):2784~2796.[4]RothmanDH.Automaticestimationoflargeresidualstaticscorrection[J]，Geophysics，1986，51(2):337~346.[5]SenMK，andStoffaPL.Nonlinearone-dimensionalseismicwaveforminversionusingsimulatedannealing[J]，Geophysics，1991，56(10):1624~1638.[6]JacobM.Pedersen，PeterD.Vestergaard，TheoZimmerman，牛毓荃，模拟退火法地震反演[C]，美国勘探地球物理学家学会第61届年会论文集，1991，sl3.3：501~504[7]SendMK，BhattacharyaBB，StoffaPL.Nonlinearinversionofresistivitysoundingdata[J]，Geophysics，1993，58(4):496~507..[8]PauloM.Carvalho，ArthurB.Weglein.使用模拟退火算法为界面多次波衰减提供子波估计[C]，美国勘探地球物理学家学会第64届年会论文集，1994，278~283[9]陈均，钱泉，唐权钧.用模拟退火算法作叠后层速度计算[J]，石油物探，1994，33（2）：11~20[10]姚姚.地球物理非线性反演模拟退火法的改进[J]，地球物理学报，1995，38（5）：643~650[11]蒋录全.地震数据反演中的Monte-Carlo算法、模拟退火算法（SA）及遗传算法（GA）[J]，西南石油学院学报，1995，17（3）：28~36[12]王山山，李青仁，管叶君.约束模拟退火反演[J]，石油地球物理勘探，1995，30