反演原理及公式介绍

1、第一章 反演理论第一节 基本概念一 反演和正演1 反演反演是一个很广的概念，根据地震波场、地球自由振荡、交变电磁场、重力场以及热学等地球物理观测数据去推测地球内部的结构形态及物质成分，来定量计算各种有关的物理参数，这些都可以归结为反演问题。在地震勘探中，反演的一个重要应用就是由地震记录得到波阻抗。有反演，还有正演。要正确理解反演问题，还要知道正演的概念。2 正演正演和反演相反，它是对一个假设的地质模型，给定某些参数（如速度、层数、厚度）用理论关系式（数学模型）推导出某种可测量的量（如地震波）。在地震勘探中，正演的一个重要应用就是制作合成地震记录。3 例子考虑地球内部的温度分布，假定地球内部的温

2、度随深度线性增加，其关系式可表示成： t(z)=a+bz 正演：给定a和b，求不同深度z的对应温度t(z) 反演：已经在不同点z测得t(z)，求a和b。二反演问题描述和公式表达的几个重要问题1应用哪种参数化方式离散的还是连续的？2地球物理数据的性质是什么？观测中的误差是什么？3问题能不能作为数学问题提出，如果能够，它是不是适定的？4对问题有无物理约束？5能获得什么类型的解，达到什么精度？要求得到近似解、解的范围、还是精确解？6问题是线性的还是非线性的？7问题是欠定的、超定的、还是适定的？8什么是问题的最好解法？9解的置信界限是什么？能否用其它方法来评价？第二节 反演的数学基础一解超定线性反问题

3、 1简单线性回归 可利用最小平方法确定参数a、b使误差的平方和最小。 （1-2-1） 拟合公式为： （1-2-2） 该方法的公式原来只适用于解超定问题，但同样适用于欠定问题，当我们有多个参数时，称为多元回归，在地球物理领域广泛采用这种方法。此过程用矩阵形式表示，则称为广义最小平方法矩阵方演。2非约束最小平方法反演广义矩阵方法由前面讨论可知，参数估计的最小平方方法用矩阵公式表示，所得到的算法等价于一个或多个模型参数的一个或多个数据集反演，步骤为：问题定义矩阵公式最小平方解 线性问题采用广义矩阵形式 d=gm （1-2-3） 对于精确的数据模型，参数m为 m=g-1d （1-2-4） 但是由于试验

4、误差，实际数据将不能精确拟合获得，故采用最小平方法求解。解的矩阵表示式为 （1-2-5）上式具体计算时可用奇异值分解方法 g=uvt 最后，得 =（gtg）-1gtd=v-1utd （1-2-6）二 约束线性最小平方反演为了得到最合适的解，通常可在方程d=gm中加先验信息，进行约束反演。约束方程为dm=h （1-2-7） d一般为只有对角线有值的矩阵，我们希望朝着偏置使得最小。 =（d-gmd-gm）+2（dm-hdm-h） （1-2-8） 如果d是单位矩阵，可以得到约束解 =（gtg+2i）（gtd+2h） （1-2-9） 式中，称为lagrange乘子。三解非线性反演问题1思路在实际工作中

5、许多问题都是非线性的，而非线性问题求解通常比较复杂，这样就产生这样一个问题，给定一些非线性问题，而它们又不服从简单的线性变换，那么能否用通用的方法使我们可以用一些线性反演的方法来估算未知模型参数，并最终求得问题的解决呢？答案是肯定的。2初始模型和线性化对于非线性问题 di=fi（m1，m2，mp）=fi（m）， i=1，2，n （1-2-10）设m0为初始模型，则其响应为 （1-2-11）现假定f（m）在m0附近是线性的，从而关于m0的模型响应的微小摄动可以用taylor级数展开为 或简记为 实际情况要考虑噪声 d=f（m）+e （1-2-12） 令y=d-f（m0），则有 e=d-=y-ax

6、 （1-2-13） e=y-ax这样，非线性问题转化成线性问题，我们可以用线性的方法求出问题的解。四、无约束非线性反演1问题的公式化目标函数： q=ete=(d-f(m)t(d-f(m) （1-2-14） 利用前述结果，上式改写为 q=ete=(y-ax)t(y-ax) （1-2-15）2问题的解法：gauss-newton法对参数摄动的最小平方解 （1-2-16） 将摄动（x=m）应用于起始模型m0，迭代公式如下： (1-2-17) 其中mk为jacobian矩阵a的赋值。3gauss-newton法的局限性当ata病态（本征值很小或近于0）时，计算的解会大到令人难以置信。因此在实践当中，必

7、须对mk做x的微小校正。4最速下降（梯度）法初始模型仅在目标函数q的负梯度方向予以校正，即 （1-2-18） 其中k是合适的常数，进一步推导可得 （1-2-19） 以上方程中以ata-1取代常数因子2k，将变为方程1-2-16所定义的gauss-newton法，k值决定校正步长。但以上方程并不含有任何逆矩阵，因此较gauss-newton法具备更好的起始收敛特征。最速下降法当采用最小平方解法时，其收敛速率将下降，因此不宜在实际反演中应用。5对非稳定性和非收敛性的补救办法当ata是病态时，为防止无界解的增大，levenberg（1944）提出了一种阻尼最小平方的方法，该方法可在taylor近似的

8、逐次应用过程中，阻滞参数摄动的绝对值。levenberg建议应在ata的主对角线上加一个随意选取的正的权因子，并且要显示出当权因子相等时，q2的剩余和的方向导数为最小。这种想法以后为maequardt（1963，1970）用来开发了一种非常有用的非线性算法。该技术称为岭回归（ridge regression）或marquardt-levenberg方法，是地球物理领域最常见的一种反演算法。就其本质来讲，实际上是gauss-newton法和最速下降法之间的内插，一种成功地结合二者有用特性的混合技术。五、约束反演：岭回归或marquardt-levenberg法1目标函数 （1-2-20） 目的：

9、误差和摄动量均取极小。其中摄动量是新增的约束条件，从本质上讲，岭回归法实际上是约束非线性最小平方法。是lagrange乘子，可认为是阻尼因子。如果赋值近于0，则其解近似于gauss-newton解。2问题的求解求解方法与非约束最小平方法相同，最终的解为： （1-2-21） 而后可将解xr用于迭代过程 （1-2-22） 其中a是k+1次迭代对mk求的值 （1-2-23）岭回归法实际上是最速下降法和gauss-newton法二者相结合的混合技术。当初始模型与问题的解相差甚远时，最速下降法起主要作用；而当接近于最终解时，最小平方法起主要作用。六非线性偏置估计对一组既不完整又不准确的数据进行解释时，通

10、常比较明智的做法是寻找一个和先验数据相一致的模型，这些先验数据可以是先前的地球物理研究数据，地质数据、测井数据，这些附加的先验信息可以帮助我们从不准确的实际数据得出的所有的解中求出最可信的一个，附有先验信息的反演问题可在一个统一的偏置估计框架内进行讨论。此方法强调实际过程的简单有效，为清楚起见，在此种方法中将初始模型和先验信息加以区别。1理论基础偏置估计的理论很简单，其基本原理类似于约束线性最小平方反演方法。特别的是除起始（或初始）模型m0外引入了先验信息h。同时，用对角线加权矩阵w=-1i来比例数据方程，使求解过程稳定。2应用先验信息的非线性反演为设有p个参数，h为先验数据，dm=h形式的约

11、束方程可表示为 （1-2-24） 为使相邻物理参数之间的差异降至最小平滑度，需采取twoneytikhonoy平滑度措施。 （1-2-25）我们的目的是要使m偏向于h，不妨将问题简单陈述为：给定一组有限的不准确的观测数据，在所有等效解中求其真解（考虑数据和模型误差）并使之与观测数据相吻合，且满足模型参数的可靠估计。从数学意义来讲，上述问题就等效于对预测误差ete和最终解与特定约束的偏差极小 （1-2-26） 如果f(m)是连续的并且可微，则可用taylor定理将其相对于初始模型m0 展开，从而给出方程(1-2-26)的线性近似（1-2-27） 令b=t，展开上式，并将偏微分置0，最后得偏置解为

12、 (1-2-28) 迭代公式 （1-2-29）如果先验信息有疑义（或不可信），那就需要将约束置为，即h=0，0，0t，而且所有的元素均置为相等的常数（0），这样所有的参数都具有相等的权重。在这种情况下，可以方便地由一单值未定乘子所取代。这样就有参数校正解 （1-2-30） 其迭代公式 （1-2-31） 因为d=1，这里2i用以控制求解的步长，而2mk有助于减小其向零矢量h的位置，我们可以将这种方法称为平滑度约束反演或最小偏置算法。3与标准方法的关系在偏置估计中，如果0，那么上述所有迭代估计公式均会简化为改进了的加权经典最小平方化式 （1-2-32） 偏置估计方法的稳定性和有效性主要取决于和d。

13、方程（1-2-30）与通常的阻尼最小平方或岭回归优化公式 （1-2-33）的不同之处在于一2m0项，唯一目的是要对参数增量的变化范围置一边界。我们可将约束反演问题定义为求反lagrange函数的极值问题： （1-2-34） 要搜寻的是最佳拟合数据的起始模型的有界摄动。在方程（1-2-29）中以量e取代wtw，我们有： （1-2-35）如果b可以统计地解释为先验参数协方差矩阵的逆，则上述方程即等效于jackson和matsuura的bayes估计方法，并类似于tarantola和valette的非线性算法。因此，应用简单代数，我们事实上已经导出一种与基于具有先验数据的概率统计处理的数学上比较严谨

14、的非线性反演法相类似的方法，但是应该注意到，tarantola和valette的里程碑方法中的反演理论和先验信息的使用均与我们的方法不尽相同。我们的主要兴趣在于迫使最终解尽可能与那些先验参数估计相一致，因此方程（1-2-35）右端最后一项不为0，因为在实际情况下，已知的先验参数估计很少。在tarantola和valette算法中，h即是实际起始模型m0。在这种情况下，正如pous等人指出的那样，方程（1-2-35）的最后一项在第一次迭代中应为0。我们将h和m0分开，这点和jackson和matsuura的作法大致相当。但在jackson和matsuura的算法中，方程（1-2-35）括号中的量

15、乘了一个适当选择的因子b(0<b<1)。如此说来，我们的简单方法或许更为通用。第三节 地震勘探中的反演方法一 地震反演的分类地震反演通常分成叠前和叠后反演两大类：叠前反演应用较少，较成熟的是avo反演；叠后反演的大量应用是波阻抗反演，这是当前地震资料处理的重要结果。从反演方法上可将地震反演划分为基于波动理论的波动方程反演和基于robinson褶积模型反演两大类，在实际工作中主要是基于robinson模型的反演。我们通常所说的地震资料波阻抗反演指的是基于robinson褶积模型的叠后地震资料反演，目前常见的有递推反演（recursive inversion），稀疏脉冲反演(spars

16、e-spike inversion)和基于模型的反演(model-based inversion)，后面两种反演方法通常称为宽带反演。递推反演包括道积分、glog、vlog、seislog、块状反演、带限反演和pivt等；稀疏脉冲反演包括多种实现方法，如模方法、最小熵方法、最大似然方法等；基于模型的反演也包括广义线性反演（gli）（cooke，1983）、地震岩性模拟（slim）（gelfand，1984）、鲁棒的速度反演方法（rovim）（fabre，1989）、宽带约束反演（bci）（martinez，1988）、parm和jason等。二. 递推反演方法1 波阻抗递推公式对于两层介质，反

17、射系数为： 分别为上下界面的密度，v1、v2为上下界面的速度。当地下为多层水平介质时，任意第i个界面的反射系数为： （1-3-1） 对应的波阻抗为： （1-3-2） 递推公式： （1-3-3） 如果用经过特殊处理的地震剖面记录道xn代替反射系数rn，则上式可写成 （1-3-4）这就是递推反演的基本公式。2低频补偿地震反射系数剖面的频带是有限的，它缺失的是高频成分和低频成分，对波阻抗反演而言，缺失高频成分只影响分辨率，而缺失低频成分就失去了速度的直流分量及速度斜坡的信息。这种波阻抗剖面通常称为相对波阻抗剖面或剩余波阻抗剖面，剩余波阻抗剖面只能反映波阻抗的相对变化而不能反映波阻抗的真实情况，因此必

18、须在剩余波阻抗剖面上，再加上合理的低频成分，进行低频补偿。（1） 利用声波测井资料补偿低频这是最常用的方法。（2） 利用叠加速度补偿低频 但叠加速度补偿因为实际三维速度场精度有限，会出现低频缺口，造成声测井曲线的值偏高和偏低的振荡。（3） 利用地质模型补偿低频这种方法比较费时。低频缺口在波阻抗反演中是常见的，有时也是比较严重的问题。所幸其横向上速度相对变化通常是正确的，仍然能确定目的层段上有意义的岩性变化。3. 一个简单应用道积分 图1-1 道积分剖面 该方法不做低频补偿，得到的是相对波阻抗。 用连续时间函数表示（1-3-4）式 （1-3-5）如果经反褶积处理后的地震道x（t）的脉冲宽度足够小

19、，认为x（t）与反射函数r（t）成比例 r（t）x（t） 则可近似求出任一时刻t的近似波阻抗 （1-3-6） k为比例系数。 具体实现步骤为： 将地震记录振幅标定到反射系数数量级 计算积分道 将积分结果转换为波阻抗 对转换结果作带通滤波得地层相对波阻抗 图1-1给出了具体的应用例子，处理资料为tjh三维工区一剖面。三 稀疏脉冲反演方法 这种方法假设地下反射系数序列是由一系列大的反射系数叠加在服从高斯分布的小反射系数背景上构成的，主要有：l1范数反褶积、最小熵方法、最大似然方法等。l1范数反褶积最早由barrodale于1973年提出，后经taylor1979年及oldenburg1983年的研

20、究，改进成为一种独特的反褶积方法，它的特点是对子波的各种相位特性都有较好的适应性。常规脉冲反褶积及预测反褶积都要假定子波是最小相位的，并且反射系数是白噪。在这两个条件下，反褶积的求解运算工作只有在最小二乘的意义下（与期输出波形均方根误差最小），才能得到一组toeplitz方程组，才能用莱文森递推法快速求解反褶积因子。误差的最小平方就被称为其范数为2。而l1范数是不做平方的判断，而用误差的绝对值之和作为标准，故称其范数为l1范数。1985年王承曙等又提出lp范数反褶积，即其判断的范数可以不是2，也不是1，而是一个任意的正整数p。由于采用了l1范数，带来的好处是对子波的相位特性放宽了限制，但是在计

21、算中没有脉冲反褶积那样简单了。它一般是从线性规划的理论出发，求解一组超定方程组的最优解。在求解过程中必须反复迭代，或者化为一组非线性方程组，用非线性规划方法迭代求解。最小熵方法由r.a.wiggins首先提出，它以方差模为判断准则，信号的规则性达到最大，熵为最小。wiggins的方差模的定义如下： (1-3-7)式中是地震道数据在某个时窗内的第个数值。当波形很突出时，达到很大振幅值，于是达到极大值，此时认为效果达到最佳。此方法对子波的相位特性不做约束，而且在一定程度上可以把混合相位子波向零相位靠拢。但该方法假设反射系数是稀疏的，只有当有少数大的脉冲存在时效果才很好，所以在具有亮点强波的剖面上，往往得到较好的反演效果。稀疏脉冲反演方法的输出为矩形波阻抗曲线形式，地层边界清晰，对厚层碳酸盐岩地区较为合适。然而其致命的弱点是要求反射系数是稀疏的，而实际上大多数地震道的反射系数是稠密的。四 基于模