仿真退火反演

5.3.1 仿真退火与系统能量最小状态

退火是众所周知的一种金属加工工艺，指金属从高温逐渐冷却的一种热处理方式。高温时，分子相对运动的热活动性大，在冷却时热活动性逐渐消失，最后形成原子的有序排列，处于热动力系统的最小能量状态。快速的突然冷却常不能达到最小能量状态，只有在缓慢冷却的条件下，在原子失去热活动性的同时有足够的时间让它们重新分布，才能达到最小能量状态。

退火这一物理过程的计算机仿真就成为一种求泛函极小的最优化算法。与退火过程类似，这种算法要求温度缓慢下降，但这与降低成本相矛盾。

仿真退火算法将非线性优化问题与统计力学中的热平衡问题进行类比，找出了优化问题求解与固体退火过程之间的相似性，从而另辟优化问题数值求解的新途径。

对于非线性优化问题，设非负目标函数为φ（x），要求

φ（v*）=minφ（vi），vi∈V （5.3.1）

其中vi为某一个可能的解，所有vi的集合记为V=｛vi｝，i=1，2，…，N；v*为V中使φ取极小的一个解估计。在与退火过程类比时，φ视为系统的内能，要使它由高到低，最后取极小而终止。vi为解空间的一个点，对应为固体内的一种粒子状态，它的扰动范围随内能的减小而逐渐减小。温度T用控制参数T′=kT类比，其中k为波尔兹曼系数。仿真退火算法在求解优化问题时先取很大的T′，T′可称为仿真温度，然后逐渐减小T′，计算时不同的迭代步对应着时间与温度T′的变化。当系统温度高内能大时粒子的随机扰动大，相应于系统状态的vi扰动的自由度也大。随迭代进行温度降低，当前解vi不断更新而形成序列，构成了对解空间的搜索。

按照统计力学中关于退火过程的理论，仿真退火算法遵循的基本规律是：某种状态vi在当前出现的概率服从波尔兹曼分布

地球物理数据处理教程

由（5.3.2）式可见，当仿真温度T′=kT→0和目标函数φ取极小时（vi=v*），系统处于内能最低的基态，这时概率f*取极大而接近于1，说明系统最稳定，熵最小。因此，仿真温度降至某一阀值以下，而且系统的熵达到最小可作为停止迭代的准则。

仿真退火算法的关键在于：在给定仿真温度T′i和当前状态vi之下，如何确定下一个状态vi+1，即如何构造迭代修改序列以便使空间中的搜索合理地进行，而最后找到使φ取全局极小的v*。仿真退火算法的创始人N.Metropolis于1953年提出了一种抽样算法。他说，假定当前的状态为v，利用随机扰动产生一种新扰动v′，其对应系统内能的变化为Δφ=φ（v′）-φ（v）。根据（5.3.2）式，这种新状态可以被接受的概率为

p=exp（-Δφ/T′） （5.3.3）

如果φ（v′）比φ（v）减小或相同，即Δφ≤0，则有概率p≥1，这种新状态毫无疑问要被接受。如果Δφ＞0，由（5.3.3）式计算p＜1，但也不一定为零。这时新状态虽然使目标函数增大，但为保证解空间的搜索是全局性的，不至于落入局部极小而无机会逃出，给予p＞0的机会是正确的，它给予从局部极小逃出的机会。根据（5.3.3）式建立了Metropolis内循环抽样算法，将在下面详述。

根据仿真退火原理，可以把算法的实现步骤大致归结如下。

第一步，选取i=0时的初始模型及仿真温度。随机选取初始地球模型v0，并依照经验选取初始温度T′0。也可通过实验选取T′0，例如，随机抽样一组模型 ｛vi｝，计算对应的目标函数 ｛φj｝，取 ｛φj｝ 的方差或最大Δφj的若干倍为T′0。

第二步，在给定vi和仿真温度T′i之下用Metropolis抽样对解空间进行搜索，具体由以下内循环迭代过程构成：

（1）令k=0时当前解估计v（0）=vi，参数T=T′i。

（2）按某一半随机准则在当前解v（k）附近产生一个邻近子集Vn（v（k））?V，它不能为空集，与v（k）的距离（或半径）与参数T成比例。在子集中随机地取出一个新的扰动模型v′作为下一步迭代的候选解，计算

Δφ=φ（v′）-φ（v（k）） （5.3.4）

（3）如Δφ≤0，则接收v′为下步迭代的当前解，即令v（k+1）=v′。如若Δφ＞0，则按（5.3.3）式计算概率p′，并计算一个在（0，1）之间的随机数p*。如果p′＞p*，则接受v′为下一步迭代的当前解，否则令v（k+1）=v（k），即保留原来的当前解估计。

（4）令k=k+1，如果k+1不超过某个规定的次数，返回（5.3.2）继续抽样搜索，否则用当前解估计v（k+1）跳出内循环，转入下述主程序的第三步。

第三步，按一定准则降温。此时令 i=i+1，Ti+1=λiTi，λi＜1。例如，λi可在（0.2，0.99）之间变化，但随i增大而增大。

第四步，判断是否停止外循环迭代，如果不是，则返回第二步，以当前vi+1=v（k+1）及Ti+1重复Metropolis抽样搜索。外循环停止计算的准则是T′小于某一规定阀值和系统的熵达到极小，或者由（5.3.3）式计算当前解的概率f→1等。

仿真退火算法的收敛性可以用马尔可夫链来研究，因为非平稳可数马尔可夫链可恰当地描述解空间的搜索序列。总的来说，在仿真退火算法中，如果仿真温度T′的变化等参数满足一定的条件收敛是有保证的，但收敛速度较慢。如果T′0足够高而且下降充分慢，以及对每一次下降，Metropolis抽样都很稳定，则能达到φ（vi）的全局极小值。这些条件正好要求加大计算成本，这是仿真退火固有的缺点。

5.3.2 仿真退火地震反演

进行地震反演时，目标泛函可选类似式（5.3.1）的任何形式。例如，对不考虑正则化的情况，以拟解选

地球物理数据处理教程

为计算第k次迭代时的目标泛函值，其中r为地震道反射系数序列，s为地震道数据，sk为第k次迭代时用 k-1次迭代结果计算合成的地震道。

初始温度T0可选为10000～20000。每一次迭代，温度下降为Tk+1=0.92Tk，这与混沌反演λk下降一致，几何级数递减。每一次迭代随机试验拟解的修改最大次数为3000。设当前第j个模型参数为rj，随机扰动的公式为

r′j=rj+Rand（I）·βTk （5.3.6）

式中Rand（I）为一个［-0.5，0.5 ］之间取值的随机数，β为与当前迭代温度有关的权。温度Tk越低时，权β应越小才能提高成功率。

按式（5.3.6）修改模型参数后作地震合成并代入（5.3.5）式求φk（r′），如果φk减小则修改成功，否则按Metropolis算法确定是否用r′代替r。接下去再作扰动试验。

对每个温度级（即每个k），规定最大扰动成功次数（如500），在此级成功修改多次后即降低温度进入（k+1）次迭代。

迭代停止的准则是：

（1）在k个温度级修改r3000次，成功率为零；

（2）温度已降到室温；

（3）|φk+1-φk|＜ε，ε为规定的小数。

按以上参数准则，利用仿真退火进行了与混沌反演相似的地震道反演，无误差结果示意如图5.7，子波为11点零相位Ricker子波，频率为30 Hz。

由图可见，仿真退火反演的分辨率是不错的，如中段低阻抗薄层有清晰显示；但是，方差较大，即抖动得厉害。

从目标泛函的下降曲线看，它呈现阶梯形状，与混沌反演的阶段性类似。在迭代230次以后φk已接近全局最小，最终在k=265时达到全局极小“山谷”，φmin=4.92×10-5。

以上是用双精度计算的反演结果。当用单精度计算时，迭代次数和φ曲线下降台阶减少，迭代次数在238次停止，φmin=4.88×10-4。

图5.7 仿真退火地震反演试验（精确数据）

（A）地层波阻抗模型，虚线为初始模型；（B）反演最终取得的波阻抗序列；（C）目标泛函随迭代的下降曲线

如果在数据中出现微小高斯噪音，则迭代226次，φmin=1.59×10-3，没有明显的全局极小“山谷”出现，此时计算时仍比混沌反演多一倍。

总之，仿真退火反演是一种非常稳定的非线性反演方法，但与混沌反演相比，优点并不突出。

二项式定理公式tk+1=Cnkan-kbk。

二项展开式的特点

1、项数展开式有共n+1项；系数：都是组合数，依次为Cn°，Cn，Cn2，Cn3等，指数的特点：a的指数由n一0（降幂）；b的指数由0一n（升幂）；a和b的指数和为n；利用二项式定理和展开式的通项公式可以求某些特殊项，如含某个幂的项、常数项、有理项、最火项等问题。

2、二项展开式中的各项的“二项式象数”与“条数”的区别，这是两个不同的概念，“二项式象数”仅指Cn0、Cn、.Cn这些组合数而言，不包括字母a、b所表示式子中的条数。通项Ckan-kbk是展开式中的第k+1项，而不是第k项。要灵活性、正确的应用二项展开式的通项公式。

通过探索二项式定理，感受由特殊到一般地认识事物的规律；在探究过程中，培养观察分析和综合、判断的能力。激发发现规律的积极性，鼓励勇于探索的精神。学生能够借助问题的引导，猜想发现、归纳并证明二项式定理，准确复述二项式定理的定义，并利用二项式定理准确展开式子。