地球物理非线性反演

1、第四章 传统非线性反演方法 （非启发式方法） 地球物理问题大多是非线性的，可以通过参数化和泰勒级数展开法将其转换为线性问题求解。 非线性问题在线性化过程中，由于舍取了高次项，使反演问题的多解性更加严重。 随着计算机发展，非线性反演方法发展很快，除传统的非线性反演方法（非启发式）以外，模拟退火法、遗传算法、神经网络法、小波分析法等启发式反演方法得到了很快发展。第四章 传统非线性反演方法（非启发式方法） 梯度Gradient 法(最速下降the steepest descent法) 共轭梯度conjugate gradient 法 牛顿法Newton 拟牛顿法（变尺度法）一、梯度法(最速下降法)

2、十九世纪中叶由Cauchy提出。 梯度法就是从一个初始模型出发沿负梯度方向收束目标函数(x)极小点的一种最优化方法。 不难理解，沿目标函数负梯度方向收索，只要步长适当，经过反复迭代，最终总可达到目标函数的极小点。一、梯度法(最速下降法)梯度法的关健1)初始模型的选择：合适的初始模型，收敛速度快；不合适则可能导致局部极小。2)步长大小的选择：太小，收敛速度太慢；太大，收敛不稳定。一、梯度法(最速下降法) 负梯度方向： (x)在任一点x的梯度：TNxxx)(,)(,)()()(21xxxxxg 负梯度方向：g(x)g(x)p一、梯度法(最速下降法) 步长大小的选择：设经i次收索后模型位置为x(i)

3、，此时负梯度方向为p(i),求步长 (i)g()g()()()()()()() 1(iiiiiiixxppxx 将目标函数在x(i)，处泰勒展开，并取 x(i)为 (i) p(i)()(2)()()()()()()(! 21)()()(iiTiiTiiiixxxxxxxx一、梯度法(最速下降法) 令(x(i+1)=0，得步长 (i)：TxxNiiNiiiiNiiiiTiiiiixxxxxxxxxxxxxxx)(,)(,)()()(21)()()(2)(1)()()(2)(1)()()()()() 1(这里p 舍去二次项及以上的项，有：)(1)()()()()()()()()()(ikNkkii

4、iTiiipxxxxxp一、梯度法(最速下降法) 按上式计算得到步长 (i)后，可沿最速下降方向计算得到新的模型位置x(i+1)：)g()g()()()()()()() 1(iiiiiiixxppxx)(1)()()()()()()()()()(ikNkkiiiTiiipxxxxxp一、梯度法(最速下降法) 计算(x(i+1),若： 成立，说明校正步长可取，继续进行下一次迭代；否则就减小如减小一半，并重复上述步骤，直至目标函数满足 为止，最终求出(x)极小点。 有时为了简便，保持步长不变。只要 不成立，迭代就结束，认为已经达到最小值。 梯度法的优点是简单，在远离极小点时收敛较快，缺点是在接近极

5、小点时，收敛很慢。 )()()() 1(iixx)()()() 1(iixx)()()() 1(iixx二、共轭梯度法 由初始模型x1沿任一下降方向p1搜索得到x2，根据p1和x2处的梯度方向找到方向p2，使沿p2可搜索得到极小值模型x*。p1p2-g2x2x1X*1、基本思路 为使p2直指极小点， p2必需满足： 满足上式的两个向量，称为共轭向量。利用共轭向量求极小值的办法，称为共轭梯度法。0, 0,0)(! 21)()()(22122112222222*22222222*22222*pHpgpppHgHgHgpTTTTT：xx：，xxxxxxx：xxx：x有并注意两边左乘必需极小泰勒展开设

6、0221pHpT1、基本思路 Hessian海森矩阵：矩阵）（海森）矩阵，（对称为（HessianxxxxxxxxxxxxxxxxxxxxxxxNNNNNNN)()(212122212121112H1、基本思路 共轭向量的求法：11211222112122121222121122121221212122212212112220,ppHppHggppHpgHpgpgppHppHgpHpgHppHpgHppHpHppgpTTTTTTTTTTTT：所以由此解得有两边左乘设1、基本思路 迭代步长求法：2222222222*2222222222222*222222222*222*0)(! 2)()(pp

7、HppgppHppgpHppgpHppgpTTTTTTTTxx：x：dxd：，xx：x：x所以由此可求得必需极小泰勒展开由1、基本思路共轭梯度法步骤：（1）由初始模型x1沿任一下降方向p1（可设为负梯度方向-g1)搜索得到x2:11gp11111111112ppHppgpTTxxx1、基本思路（2）计算p1的共轭向量p2 ，沿p2搜索得到极小值点x*:112122121222ppHpgHpgpgpTT2222222222*ppHppgpTTxxx1、基本思路 对二维二次目标函数，从任意初始点x1出发，沿任一下降方向p1作直线搜索得到x2，再从x2出发沿p1的共轭向量方向p2作直线搜索，所得到的

8、x*必是极小点。 对其它目标函数，由上面2步搜索得到的x*不一定是极小点，还需继续搜索。2、共轭梯度法步骤 共轭梯度法步骤：（1）设k=1，由初始模型x1沿任一下降方向p1（可设为负梯度方向-g1)搜索得到x2:11gp11111111112ppHppgpTTxxx2、共轭梯度法步骤（2）判断(xk+1)？若成立则搜索结束，否则继续下一步（3）计算pk的共轭向量pk+1 ，沿pk+1直线搜索下一个极小值位置xk+2kkkTkkkTkkkkkkppHppHggpgp111111111111111112kkkTkkTkkkkkkxxxppHppgp（4） k=k+1,返回第（2）步2、共轭梯度法步

9、骤 可以证明：对N维二次目标函数 从任意初始点x1出发，经上面N次迭代，最后得到的xN就是目标函数的极小点。所构成的N个搜索方向相互共轭。cxTTxbHxx21)(3、 二次目标函数搜索方向相互共轭 证明：N维二次目标函数N个搜索方向相互共轭031132332233312211222HppHppHppHpppgpHppHpppgpTTTTTT要证明)(1)(21)(121111111122ggHpHpgbpHbHgxbHxxxxcx：TT由3、 二次目标函数搜索方向相互共轭 可证得 因此：p1、 p2、 p3相互共轭，如此下去，可证明N维二次目标函数N个搜索方向相互共轭。)(112311312

10、31313gggHpgHppHpgHppTTTTT0021221313221221112223HppgpgpgpgpHpHpgHpggTTTTT：得以及由122211pgpgp又013HppT可证明：pjTgm=0，jN，则设x1=xk+1,返回（1）。否则计算搜索方向pk+1kkkkkpgggp22111（6）令k=k+1,返回（3）三、牛顿法 共轭梯度法克服了最小二乘法在接近在接近极小点时收敛慢的缺点，缺点是对二次N维目标函数这样简单的问题，也需迭代N次才能达到极小点。那么是否存在一种方法，对二次N维目标函数一次迭代就能达到极小点？牛顿法 将目标函数在x(i)，处泰勒展开，并舍取二次项以上

11、的项)()()()()()()()(2)()()()()()(! 21)()(! 21)()()(iiTiiTiiiiTiiTiiiixxxxxxxxxxxxHg三、牛顿法 对x求导数并令其为零，即可得x)(1)()()()()()()()(0)(iiiiiiiiixxxxxgHHg 新的模型位置x(i+1)：)(1)()()()() 1(iiiiiixxxxgH三、牛顿法 牛顿法的优点是对二次目标函数迭代一次就能达到极小点，这是其它方法不能比拟的；对非二次目标函数，在接近极小点时收敛较快； 缺点是在远离极小点时收敛很慢，并且要计算海森矩阵及逆矩阵，而且其逆往往病态或奇异。 牛顿法往往与梯度法

12、相配合，在远离极小点时用梯度法，接近极小点时用牛顿法。四、变尺度（拟牛顿）法 牛顿法，收敛快，但需计算海森逆矩阵Hk-1 能否用某种容易计算的矩阵Gk近似取代牛顿法中的海森逆矩阵Hk-1？从而摆脱牛顿法计算海森逆矩阵的缺点，又能保持牛顿法收敛快的优点拟牛顿法kkkkkkgHxxxx11kkkkkkgGxxxx11、变尺度法的一般形式要求Gk具有：1)对称正定，以保证搜索方向 为下降方向；2)有简单的迭代形式：3)满足拟牛顿条件。 前式可写成一般形式：kkkkktgGxx1显然Gk=I时，上式变为梯度法。kkkkkkkkkkTkTxxggGEGGgGgpg1111)(01、变尺度法的一般形式 拟

13、牛顿条件kkkkkkkkkk：xxggH)x(xHgg1111111)(得海森逆矩阵满足由因此近似取代海森逆矩阵的矩阵必需满足上述条件kkkkkkkkkkkkkkkkkkk：yGsyEsyGxxsggyxxggG或有拟牛顿条件记111111,)(1、变尺度法的一般形式拟牛顿算法的一般步骤：1)选定初始点x1,计算 1=(x1)和g1= 1,选定近似矩阵G1，如G1=I，置k=12)计算搜索方向pk=-Gkgk3)作直线搜索xk+1= xk+tkpk4)计算 k+1、gk+1、yk和sk5)判别迭代是否结束6)计算Ekyk=sk-Gkyk得Gk+17)计算k= k+1，转（2）2、变尺度法的DF

14、P算法 变尺度法的关健在于校正矩阵Ek的计算，计算校正矩阵有很多方法，其中最有效的方法之一是DFP算法 由Davidon(1950)提出，Fletcher和Powell（1963）改进得到，是无约束非线性反演最有效方法之一。是待定常数和是待定向量和kkkkTkkkTkkkkTkkkTkkkkk，vuvvuuEvvuuGG12、变尺度法的DFP算法 上式也有许多解，DFP取：kkTkkTkkkTkkTkkTkkkkTkkkkkTkkkkkkkTkkkTkkk：，：yGyyvsyyuyyGyvvsyuuyGsyvvuu22)(,)()(得两边同左乘取kkkkk：yGsyE由拟牛顿条件kkTkkkk

15、kkTkkkk，yGyyGvsysu,2、变尺度法的DFP算法 DFP算法的效正矩阵计算公式：kkTkkTkkkkTkkTkkyGyGyyGysssE 可以证明：对于具有对称正定H阵的二次目标函数，由DFP法确定的搜索方向也是相互H共轭的，因此对N维二次目标函数，DFP法经过最多N次迭代可得到极小点。第五章第五章 完全非线性反演方法完全非线性反演方法上章介绍的启发式非线性反演方法的基本思想基本思想：每一次迭代，搜索当前模型下目标函数的下降方向，沿此方向以一定步长前进，求得新模型，再以此新模型为起点进行搜索，不断迭代。如果没有下降方向了，搜索停止得到模型。若初始模型在真实模型附近，搜索能达到最小

16、值（或最大值）对应的真实模型。若离真实模型较远，在某局部极值对应的模型附近，则搜索陷入局部极值，无法再改变。（局部极值对应的不是真实模型）。能否得到真实解，强烈地依赖于初始模型的选择能否得到真实解，强烈地依赖于初始模型的选择第五章第五章 完全非线性反演方法完全非线性反演方法解决办法：解决办法： 1.利用丰富的先验知识和信息，选择较好的初始模型 2.发展不依赖于初始模型的完全非线性反演方法常用的完全非线性反演方法常用的完全非线性反演方法1.蒙特卡洛法2.模拟退火法3.遗传算法4.人工神经网络方法5.多尺度反演 第五章第五章 完全非线性反演方法完全非线性反演方法 最简单、最直接的完全非线性反演方法

17、是彻底搜索法或穷举法。即在一定约束条件下对模型参数的一切可能组合得到的模型进行分析、比较，找到某种可接受标准下的满意解或解集。相当于搜索模型空间中的所有点，即进行模型空间的彻底搜索，因此称为彻底搜索法或穷举法。 优点：只要空间中满足条件的解存在，就必然能被搜索到。 致命弱点：计算上不现实，只能是一种理论上存在的方法。一、蒙特卡洛法（Monte CarloMonte Carlo） 又称又称“尝试和误差尝试和误差”方法，它按一定的先验限制信息，以随机方法，它按一定的先验限制信息，以随机（不系统）的方式搜索模型空间，接受符合先验信息所给出的标准（不系统）的方式搜索模型空间，接受符合先验信息所给出的标

18、准的模型，的模型，“排斥排斥”并并“遗忘遗忘”不符合标准的模型。较为现实，在实不符合标准的模型。较为现实，在实际工作中得到了应用。际工作中得到了应用。优点优点：与穷举法不同，用随机抽样搜索代替了系统搜索，比较现实。 弱点弱点：不能保证搜索的彻底性，影响了它的广泛应用。 解决办法解决办法：改进常规蒙特卡洛方法。 改进的主要思路改进的主要思路：不再进行“盲目”的、完全随机的搜索，而是在一定的先验知识引导下进行随机搜索。这就是所谓的启发式蒙特卡洛反演方法。 应用效果最好的两种启发式蒙特卡洛反演方法应用效果最好的两种启发式蒙特卡洛反演方法： （1）模拟退火法：以统计物理学为基础。 （2）遗传算法：以生

19、物工程为基础。二、模拟退火法模拟退火法模拟退火法（模拟退火法（Simulated AnnealingSimulated Annealing，简称简称SASA）：）： 是Kirpatrick等1983年首先提出的。它模拟物质退火的物理过程：加热，物质温度升高，乃至熔化；冷却，温度降低，最后形成结晶。慢慢冷却，形成稳态晶体，系统处于整体能量极小状系统处于整体能量极小状态态；冷却速度太快，形成非晶体状亚稳态玻璃体，系统系统处于能量局部极小状态处于能量局部极小状态； 因此，Kirpatrick把物理系统的能量看作反演问题的目标函数，把晶体生成看作搜索到目标函数的整体极值，把玻璃体的形成看作错误地搜索到

20、局部极值。系统能量达到最小，则系统处于最优状态。二、模拟退火法模拟退火法 据统计物理学（或称统计热力学），物体中每个分子的状态服从吉布斯概率分布，即为：个分子的概率分布。为第为玻耳兹曼常数；个分子所处的状态为第个分子的能量函数；为第其中irPkirirEkTrEkTrErPiiiMjiii)()()(exp)(exp)(1为方便令k=1 退火过程中，温度降低，分子能量减小，物体从熔态向固态转变。分子取哪种状态，取决于其概率密度（大，则取该状态的可能性大）。若温度下降慢，所有分子同时按同一温度下降，则物体最后变成晶体；若温度下降太快，每个分子的温度不能同时均匀下降，则会出现亚玻璃体。（目标函数陷

21、入局部极小）（1）二、模拟退火法模拟退火法 模拟退火法用法用于反演的基本思路：将待反演的模型参数看作熔化物体的一个分子，将目标函数看作熔化物体的能量函数，通过缓慢地减小模拟温度T（温度参数至关重要），进行迭代反演，使目标函数最终达到极小。 目标函数或方差函数值大时，温度取大值，变化率可以稍大一些；目标函数值小时，温度取小值，变化率小，不至于陷入局部极小。 模拟退火法有Metropolis算法（MSA）和Heat Bath算法（HBSA）两种算法：MSA在全空间自动搜索，模型修改量是随机的；HBSA把模型参数限制在一定范围内，修改量是固定值。二、模拟退火法模拟退火法HBSA算法：1)随机地选择一

22、个初始模型m(k)(m1 (k) ，m2 (k) ，mN (k) ，取温度为T（k）2)在保持m2 (k) ，mN (k)不变的情况下，使m1 (k) 从m1 (k) min ， m1 (k) max取出所有可能的取值，按（1）式计算概率密度，并选择概率最大之m1 (k) 为其新值，记为； m1 (k+1)接着又固定m1 (k) ，m3 (k) ，mN (k)不变，按m1 (k)的方法选择m2 (k+1)的新值，依次类推，完成一次迭代得到 m(k+1)(m1 (k+1) ，m2 (k+1) ，mN (k+1) 3)取温度T(k+1) T（k） ，k=k+1返回（2）通常取T（k）=0.99kT

23、0即按指数变化。 三、三、 遗传算法遗传算法 遗传算法（遗传算法（Genetic AlgorithmGenetic Algorithm，简称简称GAGA）：）： John Holland 1975年提出的一种启发式蒙特卡洛反演方法。基于生物系统的自然选择原理和自然遗传机制，模拟自然界中的生命进化过程，解决人工系统中复杂的、特定目标的非线性反演问题。 GA算法的特点： 1.对模型群体进行追踪、搜索，即模型状态通过模型群体传送，具有比模拟退火更大、更复杂的“记忆”，潜力更大。 2.将模型参数进行二进制编码后组成一个 “串”，类似于生物遗传中的关键物质“染色体”，模拟染色体遗传基因的变化来改变模型参

24、数。三、三、 遗传算法遗传算法遗传算法的基本过程：遗传算法的基本过程： 1.模型参数进行二进制编码（组成染色体），并随机选择一组初始模型群体，； 2.通过 “选择”、“交换” 和 “变异”得到新的模型群体3.重复第二步直至模型群体变得“一致”为止。所谓群体“一致”，即群体目标函数（或后验概率）的方差或标准差很小，或者群体目标函数的均值接近于最大值。三、三、 遗传算法遗传算法 例：设m是正整数，且0m127，求方程(m)=m2的极大值。 1、模型参数的编码： 选单位为1，0127用二进制编码需7位。这样模型参数化结果有128个可选模型：1111111000001000000010000000ma

25、x1261272312min1mmmmmmmmm也可选别的单位，如=2，则可选模型减少到64个。三、三、 遗传算法遗传算法 2.2.初始模型群体产生：初始模型群体产生：随机产生初始模型群体，原则是个体在模型空间中分布均匀，最好是模型空间各代表区域均有成员。初始模型群体大，有利于搜索，但太大会增加计算量。 本例选初始模型数Q=4，即有2对父母：imi二进制(m)P (mi)11500011112250.012270100011049000.26531011100101102010.533456011100031360.170AVG4616三、三、 遗传算法遗传算法 3.3.选择：选择：选择模型配

26、对（亲本）。选择时，可取适值大的模型（即目标函数好或概率大）作为父本，淘汰适值小的模型。 本例取繁殖概率为： 其中AVG为目标函数的平均，目标函数小于它的模型不传代，大于它的可以传代。 对第一代，为了防止基因丢失，可先不舍去概率小的模型，让它与概率大的模型配对。如：本例中1与3配对，2与4配对。AVGkkQmmp)()(三、三、 遗传算法遗传算法4.4.交换：交换：将配对的两亲本模型的部分染色体相互交换，形成二个新的模型（子代）。两个染色体遗传基因的交换过程，是遗传算法的“繁殖”过程，是亲本的重组。其中交换点可随机选择。imi交换点子代mi(m)P (mi)115000-11110000101

27、5250.00092101110-01011101111111123210.463701-0001101111000120144000.53784560-11100000001106360.0013AVG6696三、三、 遗传算法遗传算法 5.5.更新更新：目标函数小于平均的模型不传代，大于平均的传代，重新选择配对：i配对变异交换点子代mi111011111101111110111-11101110110211110001111000111100-01111001121311011111101111110-1111110101070411110001111010111-101011111111

28、27AVG 6.6.变异：变异：将模型参数二进制编码的某一位由1变为0或由0变为1，使染色体基因发生变异。三、三、 遗传算法遗传算法 7.7.收敛：收敛：模型群体经多次选择、交换和变异后，大小不变，平均目标函数值逐渐变大，最后聚集在模型空间中一个小范围内。则找到了整体极大值对应的解。四、人工神经网络 人工神经网络（人工神经网络（Artificial Neural Network，简称简称ANN）：）：是一个由大量简单的处理单元组成的高度复杂的大规模非线性自适应系统。模拟人脑处理信息的模拟人脑处理信息的功能，功能，是一种非线性处理系统，十分适于进行非线性反演。 四、人工神经网络1. ANN历史历

29、史2.2. 人工神经元的基本构成人工神经元的基本构成3.3. 激活函数激活函数(Activation Function)4. 常用的人工神经网络常用的人工神经网络Hopfield神经网络 误差回传神经网络（Back-Propagation Neural Network B-P NN）1、 ANN历史n萌芽期（20世纪40年代） 人工神经网络的研究最早可以追溯到人类开始研究自己的智能的时期，到1949年止。 1943年，心理学家McCulloch和数学家Pitts建立起了著名的阈值加权和模型，简称为M-P模型。发表于数学生物物理学会刊Bulletin of Methematical Biophy

30、sics 1949年，心理学家D. O.Hebb提出神经元之间突触联系是可变的假说Hebb学习律。1、 ANN历史 第一高潮期（19501968） 人们乐观地认为几乎已经找到了智能的关键。许多部门都开始大批地投入此项研究，希望尽快占领制高点。 以Marvin Minsky，Frank Rosenblatt，Bernard Widrow等为代表人物，代表作是单级感知器（Perceptron）。n可用电子线路模拟。1、 ANN历史 反思期（19691982） M. L.Minsky和S.Papert，Perceptron，MIT Press，1969年 “异” “或”运算不可表示 二十世纪70年代

31、和80年代早期的研究结果 认识规律：认识实践再认识1、 ANN历史 第二高潮期（19831990）n1）1982年，J.Hopfield提出循环网络用Lyapunov函数作为网络性能判定的能量函数，建立ANN稳定性的判别依据 阐明了ANN与动力学的关系用非线性动力学的方法来研究ANN的特性 n指出信息被存放在网络中神经元的联接上1、 ANN历史2）1984年，J.Hopfield设计研制了后来被人们称为Hopfield网的电路。较好地解决了著名的TSP问题，找到了最佳解的近似解，引起了较大的轰动。 3）1985年，Hinton、Sejnowsky、Rumelhart等人所在的并行分布处理（PD

32、P）小组的研究者在Hopfield网络中引入了随机机制，提出所谓的Boltzmann机。 4）1986年，并行分布处理小组的Rumelhart等研究者重新独立地提出多层网络的学习算法多层网络的学习算法BP算法，较好地解决了多层网络的学习问题。（Paker1982和Werbos1974年） 国内首届神经网络大会1990年12月在北京举行1、 ANN历史 再认识与应用研究期（1991）n问题： n1）应用面还不够宽 n2）结果不够精确 n3）存在可信度的问题 n研究： n1）开发现有模型的应用，并在应用中根据实际运行情况对模型、算法加以改造，以提高网络的训练速度和运行的准确度。 n2）希望在理论上

33、寻找新的突破，建立新的专用/通用模型和算法。 n3）进一步对生物神经系统进行研究，不断地丰富对人脑的认识2、人工神经元的基本构成 生物神经网六个基本特征1）神经元及其联接； 2）神经元之间的联接强度决定信号传递的强度 3）神经元之间的联接强度是可以随训练改变 4）信号可以是起刺激作用的，也可以是起抑制作用的； 5）一个神经元接受的信号的累积效果决定该神经元的状态； 6) 每个神经元可以有一个“阈值”。2、人工神经元的基本构成 人工神经元是构成神经网络的最基本单元，应该具有生物神经元的六个基本特性。输入：X=（x1，x2，xn） 联接权：W=（w1，w2，wn）T 网络输入：net=xiwi；

34、向量形式：net=XW3、 “阈值”激活函数(Activation Function) 神经元的输出由激励函数（也称活化函数）决定：out=f(net)。 具有“阈值”性质，当net大于“阈值”时被激励；小于“阈值”时被抑制3、 激活函数(Activation Function) 常用的激活函数a、线性函数（Liner Function）b、非线性斜面函数(Ramp Function)c、阈值函数（Threshold Function）阶跃函数d、压缩函数（Squashing Function）f（net）=a+b/(1+exp(-d*net) 4、BP神经网络 误差回传神经网络（Back-Propagation Neural Network B-P NN）：一种无反馈前向网络，由输入层、输出层和隐蔽层组成。每一层中的每一个神经元均接收