模拟退火法<sup>[1，]<sup>

模拟退火算法最早在1953年由 Metropolis等人提出。在地球物理中的最早应用是Rothman在1983年利用模拟退火算法处理地震资料的剩余静校正。模拟退火法也是类似于蒙特卡洛法的随机搜索方法。但是在产生模型的过程中引入一些规则，能有效地加快搜索速度，有时又称这类方法为启发式蒙特卡洛法。

模拟退火法概念源于统计物理学，是模拟固体熔化状态逐渐缓慢冷却最终达到能量最小的结晶状态的物理过程。对于一个熔化的金属，当处于某个温度的热平衡状态时，它的每一个分子都有它可能所处的状态，有些分子可能能量高一些，有些分子可能能量低一些，分子处于何种状态的概率由分子所具有的能量决定。设分子所有可能的能级总数为n(微观粒子的能量都是量子化的，不连续的)，则分子处于某种状态的概率满足玻尔兹曼概率分布:

地球物理反演教程

其中:Ei为第i个分子的能量K为玻尔兹曼常数T为绝对温度n为分子所有可能的能级总数，分母称为配分因子pi为第i个分子处于能量Ei的概率。

如果把地球物理反演的模型向量看作分子，把目标函数看作分子的能量，把目标函数的极小值看成分子冷却结晶的最小能量，反演问题(最优化问题)可以模拟式(8.11)金属退火的过程，通过缓慢地减小温度进行反演，使目标函数(能量)逐渐达到极小值，这时所对应的模型(分子状态)就是反演结果。

为了改善于蒙特卡洛法的随机搜索方法，1953年 Metropolis等人在产生模型的过程中引入Metropolis接受准则，模型产生并不是完全随机，而是以前一个模型为基础随机产生。对能量减小的模型完全接受，对能量增加的模型按一定的概率接受，这样能有效地加快搜索速度，同时又有可能跳出局部极小值。具体如下:

设原来模型向量为mi，新的模型为mi+1(在mi基础上随机修改产生)，各自的能量(目标函数)为E(mi)和E(mi+1)。如果E(mi+1)＜E(mi)，则目标函数在减小，新模型可以接受。如果E(mi+1)＞E(mi)，则目标函数在增加，按照一定概率来确定是否接受新的模型。具体规则见式(8.12):

E(mi+1)＜E(mi) 完全接受mi+1为新模型

地球物理反演教程

式(8.12)就是Metropolis接受准则。它使得反演过程可以接受使目标函数增加的模型，因此也就使得模拟退火法有可能跳出局部极小，收敛于全局极小值点。由于玻尔兹曼常数K只是起到尺度因子的作用，在实际计算中K可取为1来简化公式。从式(8.12)可以看出，当温度较低时，pi+1/pi较小，因此接受使能量增加的新模型的可能性较小。而一般温度较低时，目标函数较小，模型比较靠近真实模型，这时基本上只接受使目标函数减小的模型，使模型尽快收敛于极小值点。

在模拟退火反演中，要求温度T随着迭代次数的增加而缓慢降温。常用的温度函数有两种。

(1)指数下降型:

Tk=T0·exp(-ck1/N) (8.13)

式中:k为迭代次数c为衰减因子N为模型参数的个数T0为初始温度。上式也可以改写为

地球物理反演教程

通常选择0.7≤α≤1。在实际应用中可采用0.5或1代替式(8.14)的1/N。图8.4(a)为指数降温曲线。采用参数为:T0=200℃，α=0.99，1/N=0.9。

(2)双曲线下降型:

T=T0αk (8.15)

式中:T0为初始温度k为迭代次数α为衰减因子，通常取0.99。初始温度T0不能取得太高，否则增加计算时间浪费机时T0也不能太低，否则模型选取不能遍及整个模型空间，只是在初始模型附近选取，不能进行全局寻优。所以T0的确定只有通过实验计算得到。图8.4(b)为双曲线降温曲线。采用参数为:T0=200℃，α=0.99。从图8.4可以看出通过对不同温度曲线和相关参数进行选择，可以控制温度下降的方式和速度。

图8.4 模拟退火法降温曲线

模拟退火法主要有三种:

(1)MSA算法(Metropolis Simulated Annealing)

(2)HBSA算法(Heat Bath Simulated Annealing)

(3)VFSA算法(Very Fast Simulated Annealing)。

图8.5 模拟退火MSA算法程序流程图

前面介绍的利用 Metropolis接受准则的算法就是经典的模拟退火法。图8.5为模拟退火 MSA算法的程序流程图。从中可以看出 MSA算法有一套模型修改准则，依次改变模型参数，每次改变都是在原来模型基础上改变一个参数，因此容易保持已有搜索成果，持续不断地向目标函数最小值点接近，因此搜索效率比蒙特卡洛法高。此外，MSA算法允许接受使目标函数增加的模型，这样又易于跳出局部极小，达到全局极小。但 MSA算法在任何温度下和蒙特卡洛法一样都是在模型全空间进行搜索，不能根据当前温度和模型减小搜索空间，此外由于模型的修改全凭运气，所以不可能像前面介绍的最小二乘法那样目标函数基本上持续减小，而是呈不规则振荡在宏观上逐渐减小，因此效率较低。

HBSA算法与 MSA算法的不同之处是在模型的修改上。也是首先随机选择一个初始M维模型向量m0(它具有M个参数)然后限制各个模型参数可能的取值范围，对取值离散化。假设每个模型参数都有N个可能的值，首先固定模型第2个参数m0(2)直到第M个参数m0(M)保持不变，只修改第1个参数m0(1)计算m0(1)的所有取值时的目标函数，然后按式(8.16)计算“概率”，它就是式(8.11)配分因子取1的公式。即

地球物理反演教程

选择“概率”最大的为模型第1个参数的修改值。照此依次对所有模型参数进行修改完成依次迭代计算。在每次迭代计算中保持温度不变。随着迭代次数增加，温度降低，最终达到稳定状态，获得最小能量解。这种方法的计算由于要计算某个参数的所有可能值，所以计算量也是很大的。

1989年Ingber提出了VFSA算法，由于速度较快，最为常用。它使得模拟退火法从理论走向了实际应用。VFSA算法在流程上与传统的模拟退火法相同，但是在模型修改、接受概率以及降温曲线上有所改进。

(1)模型修改:常规模拟退火法采用高斯随机分布修改模型，在任何温度下都是在模型全空间进行搜索。而Ingber提出采用依赖于温度的似cauchy分布产生新的模型。即

地球物理反演教程

yi=Tsgn(u-0.5)[(1+1/T|2u-1|-1](8.18)

其中:mi为当前模型第i个参数，m'i为修改后的模型参数u为[0，1]的随机数[Ai，Bi]为mi和m'i的取值范围sgn( )为符号函数。

采用以上方式能在高温下进行大范围的搜索，低温时在当前模型附近搜索，而且由于似cauchy分布具有平坦的“尾巴”，使其易于迅速跳出局部极值。这一改进大大加快了模拟退火法的收敛速度。

(2)接收概率:当E(mi+1)＞E(mi)时，VFSA算法采用如下概率接受公式:

地球物理反演教程

上式当h→1时变为式(8.12)。h通过实验获得。

(3)降温曲线(退火计划):Ingber在1989年采用式(8.13)得出指数降温曲线。从图8.4可知，温度下降较快。

总之，VFSA算法在模型修改、接受概率以及降温曲线上的改进使得模拟退火算法收敛速度大大加快。后人在此基础上还有很多的改进，读者可以参考相关文献。

模拟退火法的优点:由于不需要计算偏导数矩阵，不需要解线性方程组(当然正演计算的除外)，结构简单，易于编程此外，由于它搜索范围大，能接受较差模型，因此易于达到全局极小。缺点:随机搜索，计算量巨大，往往要计算成百上千次正演，这与前面的最小二乘法十几次的正演计算相比反演时间太长，因此一般应用在一维反演之中，在二维、三维等高维反演中应用较少。

蒙特卡洛： 大量随机抽样下的比对，最后结果就是在当前抽样数量下筛选出的一定是最想要的那个结果。举例：假如篮子里有1000个苹果(你定的测试集)，让你 每次闭着眼睛找一个最大的，可以不限制挑选次数；于是，你可以闭着眼随机拿了一个，然后再随机拿一个与第一个比，留下大的；再随机拿一个，与前次留下的比较，又 可以留下大的；循环往复这样：拿的次数越多，挑出最大苹果的可能性也就越大!但除非你把1000个苹果都挑一遍，否则你无法肯定最终挑出来的就是最大的一个。如果有 10000个苹果的话，继续如此说不定就能找到更大的！

模拟退火 ：“渐渐”清楚自己的目标是什么！并不断朝“越发”明确的目标迈进，“越来越”不被诱惑干扰。举例：为了找出地球上最高的山，一只兔子在开始并没有 合适的策略，它随机地跳了很长时间！在这期间，它可能走向高处，也可能踏入平地或沟壑。但是，随着时间的流逝，它“渐渐清醒”! 并“直直地”朝着最高的方向跳去， 最后就到达了珠穆朗玛峰。

粒子群 ：信息的社会共享，以一个团队的形式来搜索！团队里成员信息共享，共同进步；避免一个人工作时出现目光短浅，没有全局意识。举例：就像下围棋，只 专注于一个角落的战斗不一定能获取最终的胜利，只有放眼全局，把所有己方的棋子都盘活，相互间彼此帮助，才能获得最后胜利。

蚁群 ：和粒子群算法有些相似，都是靠团队的力量共同去找目标！蚁群算法中特殊的是它的"信息素"挥发! 这个效果是其他算法中没有的！

以上所有的最优化算法都很难做到极高的精度，这是必然的： 一是 因为全局搜索已经耗费了大量的时间和资源，再过分强调精度有些不经济； 二是 因为全局搜索得到的最值可以理解为一精确最值的一个准确范围！即进入这个范围再进行精确的搜索一定可以找到精确最值；但是，全局最优的核心是随机/概率，当进入一个准确范围时，这个范围肯定是很小的，如果之后精确搜索还用全局搜索的概率参数(此时来说波动范围太大了)，很可能又会跳出这个好不容易找到的精确区域！

因此： 全局最优算法与局部最优算法是要相结合的 ！全局最优算法负责划定最值所在的一个精确的、较小的范围内，即告诉局部最优算法在这个范围内继续找一定可以找到精确解；局部最优算法按照较小的步长、较高的精度继续搜索精确最值。

常用全局最优算法：蒙特卡洛(MC)、模拟退火(SA)、粒子群(PSO)、蚁群(AG)；

常用局部最优算法：梯度下降法、牛顿法、阻尼牛顿法、共轭梯度法；

推荐搭配1：蒙特卡洛

推荐搭配2：粒子群 + 梯度下降

推荐搭配3：蚁群 + 梯度下降 + 重检机制

以上提到算法的 “程序 + 详细使用说明” 参考以下地址：

优化算法

---