反演原理及公式介绍工科 第一章反演理论 第一节基本概念 一.反演和正演 1.反演 反演是一个很广的概念,根据地震波场、地球自由振荡、交变电磁场、重力场以及热学等地球物理观测数据去推测地球内部的结构形态及物质成分,来定量计算各种有关的物理参数,这些都可以归结为反演问题在地震勘探中,反演的一个重要应用就是由地震记录得到波阻抗 有反演,还有正演要正确理解反演问题,还要知道正演的概念 2.正演 正演和反演相反,它是对一个假设的地质模型,给定某些参数(如速度、层数、厚度)用理论关系式(数学模型)推导出某种可测量的量(如地震波)在地震勘探中,正演的一个重要应用就是制作合成地震记录 3.例子 考虑地球内部的温度分布,假定地球内部的温度随深度线性增加,其关系式可表示成:T(z)=a+bz 正演:给定a和b,求不同深度z的对应温度T(z) 反演:已经在不同点z测得T(z),求a和b 二.反演问题描述和公式表达的几个重要问题 1.应用哪种参数化方式——离散的还是连续的? 2.地球物理数据的性质是什么?观测中的误差是什么? 3.问题能不能作为数学问题提出,如果能够,它是不是适定的? 4.对问题有无物理约束? 5.能获得什么类型的解,达到什么精度?要求得到近似解、解的范围、还是精确解? 6.问题是线性的还是非线性的? 7.问题是欠定的、超定的、还是适定的? 8.什么是问题的最好解法? 9.解的置信界限是什么?能否用其它方法来评价? 第二节反演的数学基础 一.解超定线性反问题 1.简单线性回归 可利用最小平方法确定参数a 、b 使误差的平方和最小。
??? ? ???∑-∑∑∑-∑=-=∑∑-=2 2)()(x x n y x xy n b x b y n x b y a (1-2-1) 拟合公式为: bx a y +=? (1-2-2) 该方法的公式原来只适用于解超定问题,但同样适用于欠定问题,当我们有多个参数时,称为多元回归,在地球物理领域广泛采用这种方法此过程用矩阵形式表示,则称为广义最小平方法矩阵方演 2.非约束最小平方法反演——广义矩阵方法 由前面讨论可知,参数估计的最小平方方法用矩阵公式表示,所得到的算法等价于一个或多个模型参数的一个或多个数据集反演,步骤为: 问题定义→矩阵公式→最小平方解 线性问题采用广义矩阵形式 d=Gm (1-2-3) 对于精确的数据模型,参数m 为 m=G -1d (1-2-4) 但是由于试验误差,实际数据将不能精确拟合获得,故采用最小平方法求解解的矩阵表示式为 d G G G m T T 1][?-= (1-2-5) 上式具体计算时可用奇异值分解方法 G=U ∧V T 最后,得 m ?=(G T G )-1G T d=V ∧-1U T d (1-2-6) 二. 约束线性最小平方反演 为了得到最合适的解,通常可在方程d=Gm 中加先验信息,进行约束反演。
约束方程为 Dm=h (1-2-7) D 一般为只有对角线有值的矩阵,我们希望朝着j h 偏置j m 使得?最小 ?=(d-Gm ()T d-Gm )+β2(Dm-h ()T Dm-h ) (1-2-8) 如果D 是单位矩阵,可以得到约束解 c m ?=(G T G+β2I )1 -(G T d+β2h ) (1-2-9) 式中,β称为Lagrange 乘子 三.解非线性反演问题 1.思路 在实际工作中许多问题都是非线性的,而非线性问题求解通常比较复杂,这样就产生这样一个问题,给定一些非线性问题,而它们又不服从简单的线性变换,那么能否用通用的方法使我们可以用一些线性反演的方法来估算未知模型参数,并最终求得问题的解决呢?答案是肯定的 2.初始模型和线性化 对于非线性问题 d i =f i (m 1,m 2,…m p )=f i (m ), i=1,2,…n (1-2-10) 设m 0 为初始模型,则其响应为 )(00 m f d = (1-2-11) 现假定f (m )在m 0 附近是线性的,从而关于m 0 的模型响应的微小摄动可以用Taylor 级数展开为 高次项+??+ +??+??+??+=++++=p p i i i i i p p i m m f m m f m m f m m f m f m m m m m m m m f m f δδδδδδδδ 3322110 030 3202101)() ,,,()( 或简记为 )||(|||)()()(21000 m O m m m f m f m f p j j m m j i δ+? ?? ???∑δ??+=== 实际情况要考虑噪声 d=f (m )+e (1-2-12) ? ?????∑δ??--=-===p j j m m j i m m m f m f d m f d e 10 0.|)()()(0 令y=d-f (m 0),m x m f A j ij δ=??=,/,则有 e=d-)(m f =y-Ax (1-2-13) e=y-Ax 这样,非线性问题转化成线性问题,我们可以用线性的方法求出问题的解。
四、无约束非线性反演 1.问题的公式化 目标函数: q=e T e=(d-f(m))T (d-f(m)) (1-2-14) 利用前述结果,上式改写为 q=e T e=(y-Ax)T (y-Ax) (1-2-15) 2.问题的解法:Gauss-Newton 法 对参数摄动的最小平方解 y A A A x T T 1)(-= (1-2-16) 将摄动(x=δm )应用于起始模型m 0 ,迭代公式如下: y A A A m m T T k k 11 )(-++= (1-2-17) 其中m k 为Jacobian 矩阵A 的赋值 3.Gauss-Newton 法的局限性 当A T A 病态(本征值很小或近于0)时,计算的解会大到令人难以置信因此在实践当中,必须对m k 做x 的微小校正 4.最速下降(梯度)法 初始模型仅在目标函数q 的负梯度方向予以校正,即 ? ?? ?????-=m q k x (1-2-18) 其中k 是合适的常数,进一步推导可得 y A k m f d kA m f d A k x T T T ]2[))((2))}((2{=-=---= (1-2-19) 以上方程中以[A T A]-1 取代常数因子2k ,将变为方程1-2-16所定义的Gauss-Newton 法,k 值决定校正步长。
但以上方程并不含有任何逆矩阵,因此较Gauss-Newton 法具备更好的起始收敛特征 最速下降法当采用最小平方解法时,其收敛速率将下降,因此不宜在实际反演中应用 5.对非稳定性和非收敛性的补救办法 当A T A 是病态时,为防止无界解的增大,Levenberg (1944)提出了一种阻尼最小平方的方法,该方法可在Taylor 近似的逐次应用过程中,阻滞参数摄动的绝对值Levenberg 建议应在A T A 的主对角线上加一个随意选取的正的权因子,并且要显示出当权因子相等时,q 2 的剩余和的方向导数为最小这种想法以后为Maequardt (1963,1970)用来开发了一种非常有用的非线性算法该技术称为岭回归(Ridge Regression )或Marquardt-Levenberg 方法,是地球物理领域最常见的一种反演算法就其本质来讲,实际上是Gauss-Newton 法和最速下降法之间的内插,一种成功地结合二者有用特性的混合技术 五、约束反演:岭回归或Marquardt-Levenberg 法 1.目标函数 ) (2021 L x x e e q q T T -β+=β+=? (1-2-20) 目的:误差和摄动量均取极小。
其中摄动量是新增的约束条件,从本质上讲,岭回归法实际上是约束非线性最小平方法β是Lagrange 乘子,可认为是阻尼因子如果β赋值近于0,则其解近似于Gauss-Newton 解 2.问题的求解 求解方法与非约束最小平方法相同,最终的解为: y A I A A x T T r 1][-β+= (1-2-21) 而后可将解x r 用于迭代过程 y A I A A m m T T k k 11 ][-+β++= (1-2-22) 其中A 是k+1次迭代对m k 求的值 ][1 3210r k r k r k r k r k x x x x x m m ++++++=--- (1-2-23) 岭回归法实际上是最速下降法和Gauss-Newton 法二者相结合的混合技术当初始模型与问题的解相差甚远时,最速下降法起主要作用;而当接近于最终解时,最小平方法起主要作用 六.非线性偏置估计 对一组既不完整又不准确的数据进行解释时,通常比较明智的做法是寻找一个和先验数据相一致的模型,这些先验数据可以是先前的地球物理研究数据,地质数据、测井数据,这些附加的先验信息可以帮助我们从不准确的实际数据得出的所有的解中求出最可信的一个,附有先验信息的反演问题可在一个统一的偏置估计框架内进行讨论。
此方法强调实际过程的简单有效,为清楚起见,在此种方法中将初始模型和先验信息加以区别 1.理论基础 偏置估计的理论很简单,其基本原理类似于约束线性最小平方反演方法特别的是除起始(或初始)模型m 0 外引入了先验信息h 同时,用对角线加权矩阵W=σ-1 I 来比例数据方程,使求解过程稳定 2.应用先验信息的非线性反演 为设有p 个参数,h 为先验数据,Dm=h 形式的约束方程可表示为 ? ??? ? ???????=????????????????????? ???=p p h h h m m m Dm 2121111 (1-2-24) 为使相邻物理参数之间的差异降至最小。