反演原理及公式介绍工科

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解的置信界限是什么？能否用其它方法来评价？第二节 反演的数学基础解超定线性反问题1.

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

4、G-1d(1-2-4)但是由于试验误差，实际数据将不能精确拟合获得，故采用最小平方法求解。解的矩阵表示式为m = GTG广GTd(1-2-5)上式具体计算时可用奇异值分解方法_TG=UA V最后,得m?= (GTG) -1GTd=VA-1UTd(1-2-6)约束线性最小平方反演为了得到最合适的解，通常可在方程 d=Gm中加先验信息，进行约束反演。约束方程为Dm=h (1-2-7)D 一般为只有对角线有值的矩阵，我们希望朝着hj偏置mj使得平最小。=(d-Gm)T(d-Gm) +(Dm-h )T( Dm-h)(1-2-8)如果D是单位矩阵，可以得到约束解n?c= (GTG-+3 2I),(Gd+

5、32h)(1-2-9)式中，3称为Lagrange乘子。三.解非线性反演问题1 .思路在实际工作中许多问题都是非线性的，而非线性问题求解通常比较复杂，这样就产生这 样一个问题，给定一些非线性问题，而它们又不服从简单的线性变换，那么能否用通用的方 法使我们可以用一些线性反演的方法来估算未知模型参数，并最终求得问题的解决呢？答案 是肯定的。2 .初始模型和线性化对于非线性问题d i=fi (m, n2,m>) =fi (m), i=1 , 2,n(1-2-10)设n0为初始模型，则其响应为d0 = f(m°)(1-2-11 )现假定f (限 在m0附近是线性的，从而关于m0的模型响

6、应的微小摄动可以用 Taylor级 数展开为f (m) = fi (m0 、. mi ,m； 、. m2 ,m0 、m3, mP 、- mp)=f/m0)+-f-6 m+-f-5 m +-f-6 m + +-f-6mp+高次项fmifmi2fmh：:mp或简记为p f f( m°)f(m)= f(mo)十fmSmj> + O(|Sm|2)尸 dmljm-,实际情况要考虑噪声d=f(rnj) +e (1-2-12 )mi mj-e = d - f (m): d - f (m0) - ： fi(m ) |j m ； mj令y=d-f(n0), A= %/amj,x = 5m,则有e

7、=d-f (m) =y-Ax (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法对参数摄动的最小平方解x = (AT A)“ATy(1-2-16)将摄动(x=8m)应用于起始模型mi,迭代公式如下：mk 1 = mk (ATA)ATy (1-2-17)其中m为Jacobian矩阵A的赋值。3 . Gauss-Newton

8、法的局限性当ATA病态(本征值很小或近于 0)时，计算的解会大到令人难以置信。因此在实践当 中，必须对m做x的微小校正。4 .最速下降(梯度)法初始模型仅在目标函数q的负梯度方向予以校正，即X = -k3>(1-2-18)m,其中k是合适的常数，进一步推导可得x =一k_2AT(d - f(m) =2kAT(d - f (m) = 2kATy (1-2-19 )以上方程中以ATA-1取代常数因子2k,将变为方程1-2-16所定义的Gauss-Newton法， k值决定校正步长。但以上方程并不含有任何逆矩阵，因此较Gauss-Newton法具备更好的起 始收敛特征。最速下降法当采用最小平方

9、解法时，其收敛速率将下降，因此不宜在实际反演中应用。5 .对非稳定性和非收敛性的补救办法当ATA是病态时，为防止无界解的增大，Levenberg (1944)提出了一种阻尼最小平方的 方法，该方法可在Taylor近似的逐次应用过程中，阻滞参数摄动的绝对值。Levenberg建议 应在ATA的主对角线上加一个随意选取的正的权因子，并且要显示出当权因子相等时，q2的剩余和的方向导数为最小。这种想法以后为Maequardt (1963, 1970)用来开发了一种非常有用的非线性算法。该技术称为岭回归 (Ridge Regression )或Marquardt-Levenberg 方法, 是地球物理领

10、域最常见的一种反演算法。就其本质来讲，实际上是Gauss-Newton法和最速下降法之间的内插，一种成功地结合二者有用特性的混合技术。五、约束反演：岭回归或 Marquardt-Levenberg 法1 .目标函数9=q +P q2 = eTe+ B(xTx l0)(1-2-20)目的：误差和摄动量均取极小。其中摄动量是新增的约束条件，从本质上讲，岭回归法 实际上是约束非线性最小平方法。3是Lagrange乘子，可认为是阻尼因子。如果3赋值近于 0,则其解近似于 Gauss-Newton解。2 .问题的求解求解方法与非约束最小平方法相同，最终的解为：T1 Txr=AA+PI Ay（1-2-21

11、）而后可将解Xr用于迭代过程mk+ = mk + AT A + PI rATy（1-2-22）k 2kA1 _Xr+ Xr+X（1-2-23）Gauss-Newton法二者相结合的混合技术。当初始模型其中A是k+1次迭彳t对力求的值 k0 k kJm =mXrXr岭回归法实际上是最速下降法和 与问题的解相差甚远时，最速下降法起主要作用；而当接近于最终解时，最小平方法起主要 作用。六.非线性偏置估计对一组既不完整又不准确的数据进行解释时，通常比较明智的做法是寻找一个和先验数据相一致的模型，这些先验数据可以是先前的地球物理研究数据，地质数据、测井数据，这 些附加的先验信息可以帮助我们从不准确的实际

12、数据得出的所有的解中求出最可信的一个， 附有先验信息的反演问题可在一个统一的偏置估计框架内进行讨论。此方法强调实际过程的 简单有效，为清楚起见，在此种方法中将初始模型和先验信息加以区别。1 .理论基础偏置估计的理论很简单，其基本原理类似于约束线性最小平方反演方法。特别的是除起 始（或初始）模型m1外引入了先验信息h。同时，用对角线加权矩阵W=r-1i来比例数据方程, 使求解过程稳定。2 .应用先验信息的非线性反演为设有p个参数，h为先验数据，Dm=h式的约束方程可表示为11 Dm =m1几I I lh I m2 |h2991JLmpJ LhpJ(1-2-24 )为使相邻物理参数之间的差异降至最

13、小平滑度，需采取Twone尸Tikhonoy平滑度措施。1 -11-1Dm =m1h1I I L II m2h2.=.(1-2-25)一 . I1-1 mP_p1我们的目的是要使m偏向于h,不妨将问题简单陈述为：给定一组有限的不准确的观测 数据，在所有等效解中求其真解(考虑数据和模型误差)并使之与观测数据相吻合，且满足 模型参数的可靠估计。从数学意义来讲，上述问题就等效于对预测误差eTe和最终解与特定约束的偏差极小L =(Wd -Wf(m)T(Wd -Wf(m) ( Dm - h)T( Dm - h)(1-2-26 )如果f(m)是连续的并且可微，则可用 Taylor定理将其相对于初始模型 n

14、0展开，从而给 出方程(1-2-26)的线性近似L = (Wy -WAx)T(Wy-WAx) D(m0 x) - hT ： T D(m0 x) - h(1-2-27 )令B=3 T3 ,展开上式，并将偏微分置0,最后得偏置解为x =(WA)TWA B(WA)TWy Bh - m° (1-2-28)迭代公式mk+ = mk +(WA)TWA+ B-(WA)TWy+ Bh mk(1-2-29)如果先验信息有疑义(或不可信)，那就需要将约束置为，即 h=0, 0，0T,而且所 有3的元素均置为相等的常数(0< 3 v 1 ),这样所有的参数都具有相等的权重。在这种情 况下，3可以方便

15、地由一单值未定乘子3所取代。这样就有参数校正解xs =(WA)T +P2I(WA)TWy-P2m0(1-2-30)其迭代公式mk+ = mk +(WA)T +P2I(WA)TWy-P2mk(1-2-31)因为D=1,这里3 2I用以控制求解的步长，而3 2精有助于减小其向零矢量h的位置，我 们可以将这种方法称为平滑度约束反演或最小偏置算法。3.与标准方法的关系在偏置估计中，如果3 -0,那么上述所有迭代估计公式均会简化为改进了的加权经典 最小平方化式mk1 = mk(WA)TWA广(WA)TWy(1-2-32)偏置估计方法的稳定性和有效性主要取决于3和D。方程(1-2-30)与通常的阻尼最小平

16、方或岭回归优化公式x = (WA)TWA+ PI (WA)TWy(1-2-33)的不同之处在于一 3 2m0项，唯一目的是要对参数增量的变化范围置一边界。我们可将约束反演问题定义为求反 Lagrange函数的极值问题：L = (Wd -Wf (m)T(Wd -Wf (m) + P(m 一 m0)T(m 一 m0)( 1-2-34)要搜寻的是最佳拟合数据的起始模型的有界摄动。在方程(1-2-29)中以量E取彳tWTW,我们有：k 1kT .-1T.km =m +AEA+B A Ey + Bh m (1-2-35)如果B可以统计地解释为先验参数协方差矩阵的逆，则上述方程即等效于 Jackson和M

17、atsu' ura Bayes估计方法，并类似于 Tarantola和Valette的非线性算法。因此，应用简单 代数，我们事实上已经导出一种与基于具有先验数据的概率统计处理的数学上比较严谨的非 线性反演法相类似的方法， 但是应该注意到，Tarantola和Valette的里程碑方法中的反演理论 和先验信息的使用均与我们的方法不尽相同。我们的主要兴趣在于迫使最终解尽可能与那些先验参数估计相一致，因此方程(1-2-35)右端最后一项不为0,因为在实际情况下，已知的 先验参数估计很少。在 Tarantola和Valette算法中，h即是实际起始模型 m°。在这种情况下， 正如Po

18、us等人指出的那样，方程(1-2-35)的最后一项在第一次迭代中应为 0。我们将h和m0分开，这点和Jackson和Matsu' ura勺作法大致相当。但在Jackson和Matsu' ura勺算法中, 方程(1-2-35)括号中的量乘了一个适当选择的因子b(0<b<1)。如此说来，我们的简单方法或许更为通用。第三节地震勘探中的反演方法1 .地震反演的分类地震反演通常分成叠前和叠后反演两大类：叠前反演应用较少，较成熟的是AVO反演;叠后反演的大量应用是波阻抗反演，这是当前地震资料处理的重要结果。从反演方法上可将 地震反演划分为基于波动理论的波动方程反演和基于Robi

19、nson褶积模型反演两大类，在实际工作中主要是基于 Robinson模型的反演。我们通常所说的地震资料波阻抗反演指的是基于 Robinson褶积模型的叠后地震资料反演，目前常见的有递推反演(Recursive inversion),稀疏脉冲反演(Sparse-spike inversion和基于模型的反演(Model-based inversion),后面两种反演方法 通常称为宽带反演。递推反演包括道积分、GLOG、VLOG、SEISLOG块状反演、带限反演和PIVT等；稀疏脉冲反演包括多种实现方法，如匕模方法、最小嫡方法、最大似然方法等；基于模型的反演也包括广义线性反演 (GLI) (Coo

20、ke, 1983)、地震岩性模拟(SLIM) (Gelfand, 1984)、鲁棒的速度反演方法(ROVIM ) (Fabre, 1989)、宽带约束反演(BCI) (Martinez, 1988)、PARM 和 Jason。2 .递推反演方法1.波阻抗递推公式对于两层介质，反射系数为：R £JVV：2V2 TMP1, P2分别为上下界面的密度，Vi、V2为上下界面的速度。当地下为多层水平介质时，任意第i个界面的反射系数为:9 / 15Riz- Zi 1 iZ- ZJ i 1 J i(1-3-1)对应的波阻抗为:r r / R、Zj = Zi()(1-3-2)1- R递推公式：r r

21、 ' 1 RZn+ = ZoLI (1-3-3)21 - Rn如果用经过特殊处理白地震剖面记录道 Xn代替反射系数Rn,则上式可写成N 1Zn1 = Z0nf1Xn一 Xn(1-3-4 )这就是递推反演的基本公式。2 .低频补偿地震反射系数剖面的频带是有限的，它缺失的是高频成分和低频成分，对波阻抗反演而 言，缺失高频成分只影响分辨率，而缺失低频成分就失去了速度的直流分量及速度斜坡的信 息。这种波阻抗剖面通常称为相对波阻抗剖面或剩余波阻抗剖面，剩余波阻抗剖面只能反映 波阻抗的相对变化而不能反映波阻抗的真实情况，因此必须在剩余波阻抗剖面上，再加上合 理的低频成分，进行低频补偿。(1) 利用

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

23、t) 0cx (t)则可近似求出任一时刻t的近似波阻抗Z(t) =Z(t0)expK：x(t)dt (1-3-6)K为比例系数。具体实现步骤为：将地震记录振幅标定到反射系数数量级计算积分道kA(k)6 i i 0将积分结果转换为波阻抗Z(k) =Z°expA(k)对转换结果作带通滤波得地层相对波阻抗图1-1给出了具体的应用例子，处理资料为 TJH三维工区一剖面。三.稀疏脉冲反演方法这种方法假设地下反射系数序列是由一系列大的反射系数叠加在服从高斯分布的小反 射系数背景上构成的，主要有：L1范数反褶积、最小嫡方法、最大似然方法等。L1范数反褶积最早由Barrodale于197班提出，后经

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

25、于采用了 L1范数，带来的好处是对子波的相位特性放宽了 限制，但是在计算中没有脉冲反褶积那样简单了。它一般是从线性规划的理论出发，求解一 组超定方程组的最优解。在求解过程中必须反复迭代，或者化为一组非线性方程组，用非线 性规划方法迭代求解。最小嫡方法由RAWiggins首先提出，它以方差模为判断准则，信号的规则性达到最大， 嫡为最小。Wiggins的方差模的定义如下：ari maxXi4一卜 X2 2(1-3-7)式中xi是地震道数据在某个时窗内的第i个数值。当波形很突出时，xi4达到很大振幅值，于是V . 达到极大值，此时认为效果达到最佳。此方法对子波的相位特性不做约束， ari max而且

26、在一定程度上可以把混合相位子波向零相位靠拢。但该方法假设反射系数是稀疏的，只 有当有少数大的脉冲存在时效果才很好，所以在具有亮点强波的剖面上，往往得到较好的反 演效果。稀疏脉冲反演方法的输出为矩形波阻抗曲线形式，地层边界清晰，对厚层碳酸盐岩地区 较为合适。然而其致命的弱点是要求反射系数是稀疏的，而实际上大多数地震道的反射系数 是稠密的。13 / 15四.基于模型的反演1 .流程框图模型为基础的方法，或简称模型法，首先构造一个地质模型，并将其与地震资料进行比 较，然后利用比较的结果，迭代地更新模型，直至其与地震资料资料吻合为止。其示意图如 图 1-2。图1-2基于模型的反演示意图这种方法避免了直

27、接对地震资料进行反演，可以有较高的分辨率。然而，随之而来的是 多解性，很可能一个与地震资料吻合得很好的模型却是错的。尽管如此，由于有地震测井和 地质资料的约束，常常可以把多解性降低到最低限度，在储层横向预测和油藏描述中起重要 作用。基于模型的反演包括广义线性反演( GLI) (Cooke, 1983)、地震岩性模拟(SLIM) (Gelfand, 1984)、鲁棒的速度反演方法(ROVIM ) (Fabre, 1989)、宽带约束反演(BCI) (Martinez, 1988)、PARM和Jaso育，代表着当前反演的主流和发展趋势。2 .应用实例以 Strata 软件为例。在软件中有三种反演方法：带限(Band Limited )反演，Block 反演以及稀疏脉冲(Sparse Spike )反演。此处应用了带限反演方法，它是一种传统的递推反演计算。这种方法把地震道看成是经过零相位子波滤波后一系列的反射系数。由于地震数据中速度的低频成份已滤去，而