地壳形变反演的方法比较张永奇长安大学地质工程和测绘学院，陕西(710054)摘要：本文主要介绍了几种在地壳形变反演中的方法，有传统的蒙特卡洛法，模拟退火法， 还有比较少用的共轭梯度法介绍了三种方法的基本原理以及计算步骤，比较了他们的优缺 点，在应用时的注意事项关键字：蒙特卡洛法； 模拟退火法； 共轭梯度法引言从概率的观点来看，非线性反演方法可以分成即统计方法和确定性方法两类第一类又 称为启发式反演，如蒙特卡洛法、模拟退火法、遗传算法等就属于这一类确定性方法为第 二类，又称为非启发式反演，主要是通过迭代过程来完成非线性反演的确定性非线性反演方法采用的反演策略是非线性问题的线性化，主要利用目标函数的梯 度信息，通过反复迭代，寻找反演的最优解梯度类方法包括最速下降法、牛顿法、共轭梯 度法(ConjugateGradient, CG)、变尺度法等，也包括纯粹利用梯度信息的最小二乘、广义 逆等变种的反演方法梯度类方法的优点是收敛速度快，缺点是容易陷入局部极值，其解严 重依赖于初始猜测虽然统计类反演的最大优点是不完全依赖于初始猜测，理论上在反演过程中不会陷入局 部极值但是，这类方法的缺点也是十分明显的，就是计算工作量巨大，效率很低。
地球物 理反演中涉及的模型参数成百上千，就目前的计算条件来说，统计类方法仍然不能满足大规 模地球物理反演的要求梯度类方法计算效率高，特别是共轭梯度法经过不断的改进，并且 和统计类反演方法结合形成了统计加迭代的组合反演方法，消除了依赖于初始猜测的缺点， 成了一种广受欢迎的反演方案本文将介绍最常用的梯度类方法-共轭梯度法1 蒙特卡洛反演算法的基本原理1.1 什么是蒙特卡洛法为纪念著名的赌城一Mont eCarl人们将反演过程中任何一个阶段，用随机(或伪随 机)生器产生模型、以实现模型全空问搜索的方法统称为蒙特卡洛反演法对蒙特卡洛法之所以有兴趣，是因为它可以解决其它反演方法难以解决的许多问题诸 如，多参数反演、多个极值的高次非线性反演等假如，我们已知待求模型的参数的上下界 限，ma < mo < ma (a g IM) (1.1)inf sup式中，ma代表第a个模型参数的下界，而ma代表第a个模型参数的上界inf sup这里有两种方法对模型空间进行搜索，一种是彻底地搜索法，把模型空间允许的范围都搜索到，看那一个模型，或那一组模型的计算值(d(m))和观测数据d(m)拟合最好这种方法也叫穷举法；另一种搜索法是在模型空问允许的范围内随机地搜索，对每一个随机产生 的模型计算其理论值并把它与观测值进行比较，看其是否可以接受，这就是传统的蒙特卡洛 法。
显然，最简单最直接的完全非线性反演法是彻底搜索法，或穷举法穷举法是在一定条件下， 对模型参数进行一切可能的组合，得到大量不同的模型，并依次对这些模型进行正演计算， 得到相应的理论数据，并对这些理论数据与实际观测数据进行比较，如符合预先设定的条件， 则模型被接受，否则进行下一模型计算如此反复，直到所有模型都被验证，找到符合要求 的模型为止这种方法的优点是，只要在模型空间存在着满足条件的解，必然会被搜索出来 其致命问题是计算时间很长，长到无现实意义的地步因此穷举法只有理论上的意义，算法 上要实现，实际上是不可能的传统的蒙特卡洛法是穷举法的改进它不是对模型空间进行彻底的搜索，而是随机搜索 随机搜索能达到彻底的搜索的效果吗?回答，可以根据概率论中的大数定律，只要随机抽 样数n足够大，就可以用随机抽样的搜索结果来代替穷举法的搜索结果 蒙特卡洛法可分为传统的蒙特卡洛法和现代的蒙特卡洛法传统的蒙特卡洛法又称为‘尝试 法'（trial and error），这种方法在计算中按一定的先验信息，随机地产生大量可选择的 模型，并对这些模型进行计算，将其结果与实际观测的结果进行比较并根据预先给定的先 验信息来确定该模型是接受还是放弃。
如符合给定的先验信息，则模型被接受，否则就放弃， 并重新进行下一模型计算如此反复，直到所有模型都被验证，找到符合要求的空间模型为 止显然，传统的蒙特卡洛法不同于穷举法，虽然它只在模型空间随机地进行搜索、随机地 选择模型，但实践证明，在众多的情况下它仍能达到穷举法进行全局寻优的目的至于现代的蒙特卡洛法，如模拟退火法，遗传算法，以及我们将要讲述的其它类似的其 它方法，它们和传统的蒙特卡洛法不同，不是随机地选择模型，而是在一定的原则下，有指 导的选择模型，因此我们称它为启发式蒙特卡洛法1.2 传统的蒙特卡洛反演法的主要计算步骤1） 输人观测数据和必要的判别指数；2） 在一定的原则下，有指导的随机的生成模型；3） 对初始模型进行正演计算，并检验它是否合符要求保留合格解，撇弃不合格解；4） 如是不合格解，则返回2），如合格，则终止程序1.3 随机生成模型参数的方法简述在蒙特卡洛反演法中，随机数是在计算机中用某些计算公式产生的它占用的内存少， 速度快，便于计算但是这并不满足真正随机数的要求，严格讲，这不是随机数，而是一种 伪随机数在实际应用中只要选择得好，这样的伪随机数是可以满足要求的经过验证，这 些方法可以得到近似于［0,1］区间上均匀分布的随机序列。
常用的伪随机序列发生器有Von Neumann建议的'平方取中法'，和平方取中法的修正 ‘乘积取中法'，以及应用相当广泛的‘乘同余法'等在大多数计算机中，都有产生均匀 分布的伪随机数的内部函数或库文件可以调用，使用极其方便，这里不再赘述2模拟退火反演算法的基本原理2.1 物理退火过程金属加热到高温熔融状态成为液体，然后退火冷却成为固体的过程可以用热力物理学的微观分子运动论来进行描述在高温状态下，液体物质的大量分子或原子彼此之间进行着相 对自由移动如果该液体慢慢冷却，分子或原子的可动性就会消失，大量分子或原子常常能 够自行排列成行，形成一个纯净的晶体对于这个系统来说，晶体状态是能量最低状态；如 果液体被迅速冷却，就会达到一种具有较高能量的多晶体状态或非结晶状态(玻璃体)对于 这个系统来说，多结晶或非结晶状态不是能量最低状态2.2 Metropolis 准则固体在恒温下达到的热平衡态的过程可以用Monte Carlo算法进行模拟Monte Carlo 模拟方法比较简单，必须大量随机采样才能得到比较精确的结果，因而计算量很大为此， Metropolis等于1953年提出了效率更高的采样法，这种采样方法以概率来接受新的状态， 即新状态的接受与否要依据粒子处于该状态的几率来判断，这就是著名的Metropolis接受 准则，具体如下：设粒子的在温度T时的初始(当前)状态为i，该状态系统的能量为E。
然后假设粒子产i生一个随机扰动得到一个新的状态j，新状态系统的能量为E如果E 5，则接受新状态j，否则，舍弃新状态 j 。
这就是著名的 Metropolis 接受准则，即p (E T E )=ijE - Eexp(卄)bfE < Eijif Ei < E jij(2.4)2.3 模拟退火反演算法的基本思想模拟退火反演算法的基本思想，就是将待反演的模型的每个参数看作是熔化物体的每一 个分子，将目标函数看作是熔化物体的能量函数，通过缓慢减小一个模拟温度的控制参数来 进行迭代反演，使目标函数最终达到全局极值点不失一般性，先讨论L 2范数情况下非线性的地球物理反演问题非线性的地球物理反演问题的目标函数为①(m)二[dobs — F(m)]H + [dobs — F(m)] (2.5)式中：dobs表示观测数据向量；m为地球物理的模型参数向量；上标H表示共轭转置, 当观测数据为实数时，共轭转置退化为常规转置；F(m)表示正演问题，它满足F (m)二 d obs将非线性地球物理反演问题的每一个模型参数向量m等效为物体的某种状态r，将地 ii球物理反演问题的目标函数①(m)等效为物体的能量函数E，引入一个随迭代次数变化而i变化的控制参数T模拟物体的温度，就可以得到非线性地球物理反演的Metropolis接受准 则：P(m. T m」= 1exp[①纬-(m) j ]I bif ①(m ) < ①(m )j . (2.6)if ①(m ) >O(m ) ()ji式中： k 为玻尔兹曼常数，在地球物理反演过程中可设为1。
b至此，可以完全利用热退火的步骤进行模拟退火反演，使模型参数向量逐渐向最优模型 参数向量演化，直到地球物理反演问题的目标函数达到最小值模拟退火反演算法实质是利 用了地球物理反演问题求解过程与熔化固体退火过程的相似性，开辟了地球物理反演的新途 径2.4 温度的选取温度T的选取在模拟退火反演中是最重要的问题从式(2.6)可以看出，温度T实质上 是一种非线性加权当T值取得很大时，它对于能量差异小的模型并不能突出出来，对于能 量差异大的模型也不能完全压制因此，当温度取值很大时，无论是能量差异小还是差异大的模型选取的概率接近相等，这样就允许反演过程中，在比较大的模型空间中随机选取模型 当T值逐渐变小时，能量差异大的模型被选取的几率减小，反演只能在一个较小的模型空间 内随机选取，并随着温度的继续变小，最终得到一个反演问题的最优解模拟退火反演中 要求温度T随着迭代次数的增加而缓慢降温常用的温度退火机制有如下2种：1) 指数下降型T 二 T exp(—ck 1/n )k0式中：k为迭代次数；c为衰减因子；N为状态的个数；T为初始温度02) 双曲线下降型T 二 T (0.99)kk0式中： k 为迭代次数； T 为初始温度。
0初始温度不能取得太高，否则增加计算时间浪费机时； T 也不能太低，否则模型选取0不能遍及整个模型空间，只是在初始模型附近选取，不能进行全局寻优所以的确定只有通 过实验计算得到3共轭梯度反演算法的基本原理对于一个非线性反演问题，要寻找目标函数的极值以获取对应的极小值解作为反演的结 果通常情况下，目标函数都可以近似表示成一个多维的多项高次代数式，要求取这样目标 函数的极小值对应的解，可以从任意的初始猜测开始，沿目标函数的负梯度方向，即最速下 降方向，向前搜索即可快速到达目标函数的极小点处沿目标函数的负梯度方向搜索的方法， 称为最速下降法(Steepest Descent)理论上说，最速下降法沿着负梯度方向搜索应能够很 快找到极小解，但是在实际。