《地壳形变反演的方法比较.docx》由会员分享，可在线阅读，更多相关《地壳形变反演的方法比较.docx（14页珍藏版）》请在三一办公上搜索。
1、地壳形变反演的方法比较弓长永奇长安大学地质工程和测绘学院，陕西(710054)摘要：本文主要介绍了几种在地壳形变反演中的方法，有传统的蒙特卡洛法，模拟退火法， 还有比较少用的共轭梯度法。介绍了三种方法的基本原理以及计算步骤，比较了他们的优缺 点，在应用时的注意事项。关键字：蒙特卡洛法；模拟退火法；共轭梯度法引言从概率的观点来看，非线性反演方法可以分成即统计方法和确定性方法两类。第一类又 称为启发式反演，如蒙特卡洛法、模拟退火法、遗传算法等就属于这一类。确定性方法为第 二类，又称为非启发式反演，主要是通过迭代过程来完成非线性反演的。确定性非线性反演方法采用的反演策略是非线性问题的线性化，主要利用
2、目标函数的梯 度信息，通过反复迭代，寻找反演的最优解。梯度类方法包括最速下降法、牛顿法、共轭梯 度法(Conjugate Gradient，CG)、变尺度法等，也包括纯粹利用梯度信息的最小二乘、广义 逆等变种的反演方法。梯度类方法的优点是收敛速度快，缺点是容易陷入局部极值，其解严 重依赖于初始猜测。虽然统计类反演的最大优点是不完全依赖于初始猜测，理论上在反演过程中不会陷入局 部极值。但是，这类方法的缺点也是十分明显的，就是计算工作量巨大，效率很低。地球物 理反演中涉及的模型参数成百上千，就目前的计算条件来说，统计类方法仍然不能满足大规 模地球物理反演的要求。梯度类方法计算效率高，特别是共轭梯度
3、法经过不断的改进，并且 和统计类反演方法结合形成了统计加迭代的组合反演方法，消除了依赖于初始猜测的缺点， 成了一种广受欢迎的反演方案。本文将介绍最常用的梯度类方法共轭梯度法。1蒙特卡洛反演算法的基本原理1.1什么是蒙特卡洛法为纪念著名的赌城一MonteCarlo,人们将反演过程中任何一个阶段，用随机(或伪随 机)生器产生模型、以实现模型全空问搜索的方法统称为蒙特卡洛反演法。对蒙特卡洛法之所以有兴趣，是因为它可以解决其它反演方法难以解决的许多问题。诸 如，多参数反演、多个极值的高次非线性反演等。假如，我们已知待求模型的参数的上下界 限，ma mo ma(a g IM)(1.1)式中，m f代表第
4、a个模型参数的下界，而msu代表第a个模型参数的上界。这里有两种方法对模型空间进行搜索，一种是彻底地搜索法，把模型空间允许的范围都 搜索到，看那一个模型，或那一组模型的计算值(d(m)和观测数据d(m)拟合最好。这种 方法也叫穷举法；另一种搜索法是在模型空问允许的范围内随机地搜索，对每一个随机产生 的模型计算其理论值并把它与观测值进行比较，看其是否可以接受，这就是传统的蒙特卡洛 法。显然，最简单最直接的完全非线性反演法是彻底搜索法，或穷举法。穷举法是在一定条件下， 对模型参数进行一切可能的组合，得到大量不同的模型，并依次对这些模型进行正演计算， 得到相应的理论数据，并对这些理论数据与实际观测数
5、据进行比较，如符合预先设定的条件， 则模型被接受，否则进行下一模型计算。如此反复，直到所有模型都被验证，找到符合要求 的模型为止。这种方法的优点是，只要在模型空间存在着满足条件的解，必然会被搜索出来。 其致命问题是计算时间很长，长到无现实意义的地步。因此穷举法只有理论上的意义，算法 上要实现，实际上是不可能的。传统的蒙特卡洛法是穷举法的改进。它不是对模型空间进行彻底的搜索，而是随机搜索。 随机搜索能达到彻底的搜索的效果吗?回答，可以。根据概率论中的大数定律，只要随机抽 样数n足够大，就可以用随机抽样的搜索结果来代替穷举法的搜索结果。蒙特卡洛法可分为传统的蒙特卡洛法和现代的蒙特卡洛法。传统的蒙特
6、卡洛法又称为尝试 法 （trial and error），这种方法在计算中按一定的先验信息，随机地产生大量可选择的 模型，并对这些模型进行计算，将其结果与实际观测的结果进行比较。并根据预先给定的先 验信息来确定该模型是接受还是放弃。如符合给定的先验信息，则模型被接受，否则就放弃， 并重新进行下一模型计算。如此反复，直到所有模型都被验证，找到符合要求的空间模型为 止。显然，传统的蒙特卡洛法不同于穷举法，虽然它只在模型空间随机地进行搜索、随机地 选择模型，但实践证明，在众多的情况下它仍能达到穷举法进行全局寻优的目的。至于现代的蒙特卡洛法，如模拟退火法，遗传算法，以及我们将要讲述的其它类似的其 它方
7、法，它们和传统的蒙特卡洛法不同，不是随机地选择模型，而是在一定的原则下，有指 导的选择模型，因此我们称它为启发式蒙特卡洛法。1.2传统的蒙特卡洛反演法的主要计算步骤1）输人观测数据和必要的判别指数；2）在一定的原则下，有指导的随机的生成模型；3）对初始模型进行正演计算，并检验它是否合符要求。保留合格解，撇弃不合格解；4）如是不合格解，则返回2），如合格，则终止程序。1.3随机生成模型参数的方法简述在蒙特卡洛反演法中，随机数是在计算机中用某些计算公式产生的。它占用的内存少， 速度快，便于计算。但是这并不满足真正随机数的要求，严格讲，这不是随机数，而是一种 伪随机数。在实际应用中只要选择得好，这样
8、的伪随机数是可以满足要求的。经过验证，这些方法可以得到近似于0,1区间上均匀分布的随机序列。常用的伪随机序列发生器有Von Neumann建议的平方取中法，和平方取中法的修正 乘积取中法，以及应用相当广泛的乘同余法等。在大多数计算机中，都有产生均匀 分布的伪随机数的内部函数或库文件可以调用，使用极其方便，这里不再赘述。2模拟退火反演算法的基本原理2.1物理退火过程金属加热到高温熔融状态成为液体，然后退火冷却成为固体的过程可以用热力物理学的微观分子运动论来进行描述。在高温状态下，液体物质的大量分子或原子彼此之间进行着相 对自由移动。如果该液体慢慢冷却，分子或原子的可动性就会消失，大量分子或原子常
9、常能 够自行排列成行，形成一个纯净的晶体。对于这个系统来说，晶体状态是能量最低状态；如 果液体被迅速冷却，就会达到一种具有较高能量的多晶体状态或非结晶状态(玻璃体)。对于 这个系统来说，多结晶或非结晶状态不是能量最低状态。2.2 Metropolis 准则固体在恒温下达到的热平衡态的过程可以用Monte Carlo算法进行模拟。Monte Carlo 模拟方法比较简单，必须大量随机采样才能得到比较精确的结果，因而计算量很大。为此， Metropolis等于1953年提出了效率更高的采样法，这种采样方法以概率来接受新的状态， 即新状态的接受与否要依据粒子处于该状态的几率来判断，这就是著名的Met
10、ropolis接受 准则，具体如下：设粒子的在温度T时的初始(当前)状态为，该状态系统的能量为E。然后假设粒子产 生一个随机扰动得到一个新的状态J，新状态系统的能量为E j。如果Ej 5，则接受新状态j，否则，舍弃新状态j。这就是著名的Metropolis接受准则，即exp(ifE. EifEi E：(2.4)2.3模拟退火反演算法的基本思想模拟退火反演算法的基本思想，就是将待反演的模型的每个参数看作是熔化物体的每一 个分子，将目标函数看作是熔化物体的能量函数，通过缓慢减小一个模拟温度的控制参数来 进行迭代反演，使目标函数最终达到全局极值点。不失一般性，先讨论L 2范数情况下非线性的地球物理反
11、演问题。非线性的地球物理反 演问题的目标函数为中(m) = dobs - F(m)H + dobs - F(m)(2.5)式中：dobs表示观测数据向量；m为地球物理的模型参数向量；上标H表示共轭转置，当观测数据为实数时，共轭转置退化为常规转置；F(m)表示正演问题，它满足F (m) = dobs将非线性地球物理反演问题的每一个模型参数向量mi等效为物体的某种状态r，将地 球物理反演问题的目标函数中(m)等效为物体的能量函数E,，引入一个随迭代次数变化而 变化的控制参数T模拟物体的温度，就可以得到非线性地球物理反演的Metropolis接受准则：1P叫-mj)= jexpkbT中(m) (m)
12、 ij (2.6)if (m ) (m )在地球物理反演过程中可设为1。式中：kb为玻尔兹曼常数至此，可以完全利用热退火的步骤进行模拟退火反演，使模型参数向量逐渐向最优模型 参数向量演化，直到地球物理反演问题的目标函数达到最小值。模拟退火反演算法实质是利 用了地球物理反演问题求解过程与熔化固体退火过程的相似性，开辟了地球物理反演的新途 径。2.4温度的选取温度T的选取在模拟退火反演中是最重要的问题。从式(2.6)可以看出，温度T实质上 是一种非线性加权。当T值取得很大时，它对于能量差异小的模型并不能突出出来，对于能 量差异大的模型也不能完全压制。因此，当温度取值很大时，无论是能量差异小还是差异
13、大 的模型选取的概率接近相等，这样就允许反演过程中，在比较大的模型空间中随机选取模型。 当T值逐渐变小时，能量差异大的模型被选取的几率减小，反演只能在一个较小的模型空间 内随机选取，并随着温度的继续变小，最终得到一个反演问题的最优解。模拟退火反演中， 要求温度T随着迭代次数的增加而缓慢降温。常用的温度退火机制有如下2种：1) 指数下降型T = T exp(-ck 1/ n )式中：k为迭代次数；c为衰减因子；N为状态的个数；孔为初始温度。2) 双曲线下降型T = T (0.99)k式中：k为迭代次数；t为初始温度。初始温度不能取得太高，否则增加计算时间浪费机时；t也不能太低，否则模型选取 不能
14、遍及整个模型空间，只是在初始模型附近选取，不能进行全局寻优。所以的确定只有通 过实验计算得到。3共轭梯度反演算法的基本原理对于一个非线性反演问题，要寻找目标函数的极值以获取对应的极小值解作为反演的结 果。通常情况下，目标函数都可以近似表示成一个多维的多项高次代数式，要求取这样目标 函数的极小值对应的解，可以从任意的初始猜测开始，沿目标函数的负梯度方向，即最速下 降方向，向前搜索即可快速到达目标函数的极小点处。沿目标函数的负梯度方向搜索的方法， 称为最速下降法(Steepest Descent)。理论上说，最速下降法沿着负梯度方向搜索应能够很 快找到极小解，但是在实际的计算中发现，当迭代接近于极
15、值点时，会增加成百上千次的搜 索也难于到达极值点，这主要缘于两个问题：在极值点附近，由于梯度值逐渐减小，迭代 步长越来越短，将迅速增加迭代次数；在极值点附近，由于梯度值逐渐减小，很小的数值 误差，都有可能改变指向极值点的迭代方向，从而增加迭代次数。这两个因素使得最速下降 法变得不具有实用性。有没有一种方案使得在有限迭代次数内，搜索到目标函数的极值点呢?回答是肯定的。 当目标函数是n维二次多项式时，如果迭代的搜索方向沿共轭方向搜索时，能够在n次迭代 过程中，到达极值点。下面我们来介绍共轭梯度法的原理及实现。3.1下降法的搜索方向设有如下的反演目标函数/、1，中(x) = xtHx - bTx +
16、 c(3 1)2其中x表示模型参数，x和b均是n维向量，c是一个常数，h是一个n x n的正定对称矩阵。我们要求取的解就是目标函数对应极值点的解x *。式(3.1)是一个标准的二次型，它对应的极值解就是(3.2)中，(x) = V(x) = Hx - b = 0即Hx = bx* = H -ib我们假设在点x0处开始沿负梯度方向r0 = r(x) = -p(x)搜索，到达点x 1，即x =x +ar为了使搜索能够快速到达极值点，我们这样选取a使甲(x1)=甲(x0 + ar0)达到最小， 即取?甲(x ) = 0,或者：d 甲(x ) = d 甲(x +ar ) *(x) r = 0daida
17、idao o i o上式说明，如果a的取值，使得中(气)与ro正交，即前一次搜索的方向必与下一次的 搜索方向正交，这样使得最速下降的搜索路径成空间锯齿形，且下降方向是一组交替平行的 梯度方向(图3.1)。这启示我们，虽然最速下降的搜索过程有成百上千次，但实际上只需要 一组n个(维)彼此正交的梯度方向就可以搜索到极值解。然而，在设计实际算法时，从一组 n个(维)彼此正交的梯度方向，直接计算出每一方向上的搜索步长是相当困难的。但是，如 果把正交搜索方向改为共轭方向，算法就很容易实现。3.2共轭方向设P1和P2为二个n维向量，如果存在一个n x n的正定对称矩阵H，使得：p：Hp2 = 0(3.3)
18、就称P1和P2是关于矩阵H共轭的。如果矩阵H是一个单位阵，即H = I，则pTHp = pTlp = pT p = 0(3.4)121212这就是我们熟知的正交向量的定义。所以，我们可以把共轭看成是正交的推广，或者把 正交看成是共轭的特殊情况。我们把彼此互相共轭的p1和p2所代表的n维空间的两个方向 称为共轭方向。图3.1二维空间情形下的最速下降法搜索路径3.3共轭向量的构造沿共轭方向的搜索，关键是要构造出n个彼此共轭的方向向量。为了说明，共轭梯度法 的原理，我们先考虑形式如式(3.1)表示的目标函数。设有一组n维彼此关于nx n的正定对 称矩阵H共轭的向量p,P1,LL,P1，能够使我们分别
19、沿着这n个共轭向量所指的方向各搜 索一次，就可以达到极值点x *。设有一组线性无关的向量g0, g1,LL, gn1，可以通过对它线性组合构造出一组n个彼此 共轭的向量。不妨设p o = g ，为此构造与p 0共轭的向量P1:P1 = g +a p(3.5)其中a1为常数。由于p与p1关于H共轭的向量，故：pTHp = pTH(g +Op ) = pTHg +a prHp = 0(3.6)00000 00所以，prHg_ PH00(3.7)因而P = g - ；H 1 P11 pT Hp 0(3.8)由于g 0=p0，上式也可以写成:p =g - g*1 pi 1g THg 0 0在构造了共轭
20、向量p0 与 p 1的基础上，设我们已经求出了前k个彼此共轭的向量P , P,LL,P- 1，构造一个向量pk与p,%LL,P- 1都H共轭，并先设:k-1p = g +乙 a pi=0(3.9)欲使 p 与 P , P, LL, P 都 H 共轭，便有：pTHp = 0, (i = 0,1, LL, k -1)k 0 1n-1i k将式(9 )代人上式，得-piHgk , (i = 0,1, LL, k -1) p T Hppk=gk-歹1 ElppT Hp ii=1 i i（3.10）（3.11）依此类推，就可以求得所有一组彼此关于丑共轭的所有向量% PLL, Pn-1。3.4共轭梯度反演
21、方法在上一节中，我们看到构造一组彼此共轭的向量并不难，只需要找到一组线性无关的向量g ,g ,LL,g 即可。在3.1节中，我们知道，最速下降法在搜索的过程中，可以比较容 0 1n-1易地计算出一组n个互相正交的梯度方向向量，我们不妨取g ,g ,LL,g 作为这一组n个 0 1n-1互相正交的梯度方向向量，来计算出共轭向量。为此，我们构造沿共轭方向搜索的反演算法 如下：(3.12)（%）=所（%）= （IF，* LL, /）12n甚至，我们可以从上式构造一组正交向量，不失一般性，取:g o=(罗,0, LL,0)力平g1= (0,产,LL,0)2g = (0,0, LL,-)n-1oxn很显
22、然，上述的g0,g1,LL,gn-1向量，是一组线性无关的正交向量组，符合构造共轭向 量组的条件。首先，令0= g 02) 置循环变量以= 1,2, LL, n -1)，按照公式(3.9)计算出共轭向量P 1到P- 1，3) 从初始搜索点x0开始，沿共轭向量P0,P1,LL,P所指向的方向依次搜索。每步搜索 的计算公式是：x = x + P p ,(k = 1,2, LL, n -1)(3.13)kk-1k-1 k-1每次迭代的修改量为Ax = x -x =P p 。kk-1k-1 k-1最终，经过n步搜索，到达极值点x* = x = x +罗1 P p(3.14)i=0为了求取P k，(k
23、= 0,1,LL,n -1)，我们根据求极值的定义，在每次沿p搜索时，有000p*(气)=函、+ Ppk-1)*(xk)Tpk-1 = 0根据式(3.2)有：中(x )t p = (Hx b)t p = 0(H (x + P p ) b) t p = 0k-1k-1 k-1k-1从而P = (% f) % =叮 1 Pk 1(3.15)k(pT Hp ) pT Hpk 1 k 1k 1 k 1如果把式(3.14)，代人(3.2)可以证明x*就是目标函数(x)的极小值，同时也是方程Hx = b的解。在n维目标函数为标准二次型的情况下，共轭梯度法经过n次迭代，能获取 精确解的性质被称为二次截止性。
24、3.4非线性共轭梯度反演当目标函数是高于二次的连续函数(即目标函数的梯度存在)时，其对应的解方程是非线 性方程，非线性问题的目标函数可能存在局部极值，并且破坏了二次截止性，共轭梯度法需 要在两个方面加以改进后，仍然可以用于实际的反演计算，但共轭梯度法不能确保收敛到全 局极值。1) 首先是共轭梯度法不能在维空间内依靠n步搜索到达极值点，需要重启共轭梯度 法，继续迭代，以完成搜索极值点的工作。2) 在目标函数复杂，用式(3.10)计算a时，由于需要局部线性化，需计算海赛(Hessian) 矩阵H，且计算工作量比较大，矩阵H也有可能是病态的。一些研究者给出了不同的a计 算方案，但几乎都抛弃了矩阵H的
25、计算，其中以Fletcher和Reeves(1964)的方案最为常 用，具体形式如下：a -gLi gk-1(3.16) 式中gk-1和g分别为第k -1和看k次搜索时计算出来的目标函数的梯度。以Fletcher和Reeves计算方案为基础的共轭梯度反演算法如下：1)给定n维目标函数中3)，初值X = %，最大迭代次数max，误差8 V L计算 ；计算 8 newi V i82)右 max，且执行newg0=-VX0)，置i = ig = g，p = g82g僧执行(3)，否则退出，并输出计算结果Xg = (V甲(x) t p3)计算pT V 2 甲(X)pg =-Vp (x), 6 dd =
26、84).8,8= gTg, ad = new, p = g +adp, k = k +1new new8dd如果是k = n或者gTp - ，置p = g, k = 5) = +1，返回第(2)步。4存在的问题及进一步的展望本文只是简单的介绍了三种在地壳形变中的反演方法的原理，以及在进行反演时的主要 步骤。从而对三种反演方法有了一些初步的认识，但是本文没有结合地壳形变监测进行具体 的算例分析，在今后的工作中要将这三种方法应用在实际中，从而在理性和感性上重新认识 它们。这些方法现在在地球物理反演中都有广泛的应用，在测绘领域尤其是地壳形变反演中也 有重要的作用，所以掌握这些方法在今后的学习和工作中会事倍功半。参考文献1 王家映.地球物理反演理论M.北京：高教出版社，2002.2 姚姚.蒙特卡洛非线性反演方法及其应用M.北京：冶金工业出版社，1997.3 王家映.地球物理反演问题概述3.工程地球物理学报，2007, 4(1)： 13.4 杨文采.线性地球物理反演方法回顾与展望M.地球物理学进展，2002, 17(2)： 255 261.5 王家映.蒙特卡洛法3.工程地球物理学报，2007，4(2)： 8185.6 朱培民，王家映.随机共轭梯度反演法J .石油地球物理勘探，2000，35(2)： 208 213.
