反演原理及公式介绍.docx

第一章反演理论第一节基本概念一. 反演和正演1. 反演反演是一个很广的概念，根据地震波场、地球自由振荡、交变电磁场、重力场以及热学 等地球物理观测数据去推测地球内部的结构形态及物质成分，来定量计算各种有关的物理参 数，这些都可以归结为反演问题。在地震勘探中，反演的一个重要应用就是由地震记录得到 波阻抗。有反演，还有正演。要正确理解反演问题，还要知道正演的概念。2. 正演正演和反演相反，它是对一个假设的地质模型，给定某些参数(如速度、层数、厚度) 用理论关系式(数学模型)推导出某种可测量的量(如地震波】在地震勘探中，正演的一个 重要应用就是制作合成地震记录。3. 例子考虑地球内部的温度分布，假定地球内部的温度随深度线性增加，其关系式可表示成：T(z)=a+bz正演：给定a和b，求不同深度z的对应温度T(z)反演：已经在不同,曲测得T(z)，求a和b。二. 反演问题描述和公式表达的几个重要问题1. 应用哪种参数化方式一一离散的还是连续的？2. 地球物理数据的性质是什么？观测中的误差是什么？3. 问题能不能作为数学问题提出，如果能够，它是不是适定的？4. 对问题有无物理约束？5. 能获得什么类型的解，达到什么精度？要求得到近似解、解的范围、还是精确解？6. 问题是线性的还是非线性的？7. 问题是欠定的、超定的、还是适定的？8. 什么是问题的最好解法？9. 解的置信界限是什么？能否用其它方法来评价？第二节反演的数学基础一. 解超定线性反问题1 .简单线性回归可利用最小平方法确定参麴、b使误差的平方和最小。'_Zy b£x_ L a 一一 y bxn(1-2-1)s ，一 一一、b =(n E-Z x A)nE x2 - (Z x)2拟合公式为：y = a + bx(1-2-2)该方法的公式原来只适用于解超定问题但同样适用于欠定问题，当我们有多个参数时， 称为多元回归，在地球物理领域广泛采用这种方法。此过程用矩阵形式表示，则称为广义最 小平方法矩阵方演。2. 非约束最小平方法反演一一广义矩阵方法由前面讨论可知，参数估计的最小平方方法用矩阵公式表示，所得到的算法等价于一个 或多个模型参数的一个或多个数据集反演，步骤为：问题定义一矩阵公式一最小平方解线性问题采用广义矩阵形式d=Gm(1-2-3)对于精确的数据模型，参数n为m=G-1d(1-2-4)但是由于试验误差，实际数据将不能精确拟合获得，故采用最小平方法求解。解的矩阵 表示式为m = GtG -1 Grd(1-2-5)上式具体计算时可用奇异值分解方法G=UAVT最后，得人，、，、m = (GtG) -1GTd=VA-1UTd(1-2-6)二. 约束线性最小平方反演为了得到最合适的解，通常可在方fS=Gm中加先验信息，进行约束反演。 约束方程为Dm=h(1-2-7)D 一般为只有对角线有值的矩阵，我们希望朝着偏置m使得中最小。 j j中=(d-Gm)( d-Gm) +B2 (Dm-h)r( Dm-h)(1-2-8)如果D是单位矩阵，可以得到约束解m = (GTG+B2I) -1 (GTd+B2h)(1-2-9)c式中，B称为Lagrange乘子。三. 解非线性反演问题1. 思路在实际工作中许多问题都是非线性的，而非线性问题求解通常比较复杂，这样就产生这 样一个问题，给定一些非线性问题，而它们又不服从简单的线性变换，那么能否用通用的方 法使我们可以用一些线性反演的方法来估算未知模型参数，并最终求得问题的解决呢？答案 是肯定的。2. 初始模型和线性化对于非线性问题d=f (m，m，皿)=f (m)，i=1， 2，n(1-2-10)i i 12p i设mo为初始模型，则其响应为d 0 = / (m0)(1-2-11)现假定f (m)在mo附近是线性的，从而关于mo的模型响应的微小摄动可以用Taylor级 数展开为f (m) = f (mo +8 m ,mo +8 m ,mo +8 m ,mo +8 m )=f (mo) + 直8 m + 直8 m +直8 m + f 8 m + 高次项1 dm 问题的公式化目标函数：q=ee=(df(m)r(df(m)(1-2-14)利用前述结果，上式改写为q=Te=(y-Ax)T (y-Ax)(1-2-15) dm 问题的解法：Gauss-Newton法 3m 对参数摄动的最小平方解X =( AtA)-1 Aiy(1-2-16)将摄动(x=5m)应用于起始模型n。，迭代公式如下：3m p123p或简记为'df 帝)、f (m) = f (m ) + 2 f _) I 8 m + 0(118 mll2)o Idmm=mojV.J=1.j实际情况要考虑噪声 d=f(m) +e(1-2-12)'df 帝)、e = d - f (m) = d - f (mo) - 2 (m)I 8 m. >dmm=m0ji j=1.j令 y=d-f(m0)，A = df /dm.,X = 8m，则有e=df (m) =y-Ax(1-2-13)e=y-Ax这样，非线性问题转化成线性问题，我们可以用线性的方法求出问题的解。四、无约束非线性反演mk+i = mk + ( At A)-1 Ary(1-2-17)其中mk为Jacobian矩阵A的赋值。3. Gauss-Newton法的局限性当AtA病态(本征值很小或近于)时，计算的解会大到令人难以置信。因此在实践当 中，必须对mk做x的微小校正。4. 最速下降(梯度)法初始模型仅在目标函麴的负梯度方向予以校正，即工=一 "¥ I(1-2-18)dm J其中k是合适的常数，进一步推导可得尤=-k-2At(d - f (m) = 2kAT(d - f (m) = 2kAry(1-2-19)以上方程中以AtA-1取代常数因子2k，将变为方程1-2-16所定义的Gauss-Newton法, k值决定校正步长。但以上方程并不含有任何逆矩阵因此较Gauss-Newton法具备更好的起 始收敛特征。最速下降法当采用最小平方解法时，其收敛速率将下降，因此不宜在实际反演中应用。5. 对非稳定性和非收敛性的补救办法当AtA是病态时，为防止无界解的增大，Levenberg (1944)提出了一种阻尼最小平方的 方法，该方法可在Taylor近似的逐次应用过程中，阻滞参数摄动的绝对值。Levenberg建议 应在AtA的主对角线上加一个随意选取的正的权因子，并且要显示出当权因子相等时2,的 剩余和的方向导数为最小。这种想法以后Maequardt (1963， 1970)用来开发了一种非常 有用的非线性算法。该技术称为岭回归(RidgeRegression)或Marquardt-Levenberg方法， 是地球物理领域最常见的一种反演算法就其本质来讲，实际上是Gauss-Newton法和最速下 降法之间的内插，一种成功地结合二者有用特性的混合技术。五、约束反演：岭回归或(arquardt-Levenberg法1. 目标函数中=q 邛 q = ere +。(心工一L)(1-2-20)120目的：误差和摄动量均取极小。其中摄动量是新增的约束条件，从本质上讲，岭回归法 实际上是约束非线性最小平方法。B是Lagrange乘子，可认为是阻尼因子。如果B赋值近于 0,则其解近似 于Gauss-Newton解o2. 问题的求解求解方法与非约束最小平方法相同，最终的解为：尤=4A + P T-i A"(1-2-21)r而后可将解x用于迭代过程 rmk+i = mk +AM + p I -i A(1-2-22)其中A是k+1次迭代对mk求的值mk m0 + Xk + Xk-1 + Xk-2 + Xk-3 + Xi ( 1-2-23)岭回归法实际上是最速下降法和Gauss-Newton法二者相结合的混合技术。当初始模型 与问题的解相差甚远时，最速下降法起主要作用；而当接近于最终解时，最小平方法起主要 作用。六. 非线性偏置估计对一组既不完整又不准确的数据进行解释时通常比较明智的做法是寻找一个和先验数 据相一致的模型，这些先验数据可以是先前的地球物理研究数据，地质数据、测井数据，这 些附加的先验信息可以帮助我们从不准确的实际数据得出的所有的解中求出最可信的一个， 附有先验信息的反演问题可在一个统一的偏置估计框架内进行讨论。此方法强调实际过程的 简单有效，为清楚起见，在此种方法中将初始模型和先验信息加以区别。1. 理论基础偏置估计的理论很简单，其基本原理类似于约束线性最小平方反演方法。特别的是除起 始(或初始)模型mo外引入了先验信,tho同时，用对角线加权矩阵W=o-1I来比例数据方程， 使求解过程稳定。2. 应用先验信息的非线性反演为设有p个参数，h为先验数据，Dm=h形式的约束方程可表示为1_1Dm =m 1h 111mh2=2.mlhpp(1-2-24)为使相邻物理参数之间的差异降至最小平滑度，需采WoneyTikhonoy平滑度措施。1 -1Dm =1-11 -1mh11mh.2 _.mhpp(1-2-25)我们的目的是要使皿偏向于h，不妨将问题简单陈述为：给定一组有限的不准确的观测 数据，在所有等效解中求其真解(考虑数据和模型误差)并使之与观测数据相吻合，且满足 模型参数的可靠估计。从数学意义来讲，上述问题就等效于对预测误差e和最终解与特定 约束的偏差极小L =(Wd - Wf (m)T (Wd - Wf (m) + (。Dm - h)r (饥Dm - h)(1-2-26)如果f(m)是连续的并且可微，则可用Taylor定理将其相对于初始模型1。展开，从而给 出方程(1-2-26)的线性近似L = (Wy - WAx)t(Wy - WAx) + D(m0 + x) - hTpp D(mo + x) - h(1-2-27)令B=BtB，展开上式，并将偏微分置)，最后得偏置解为x = (WA)tWA + B-1(WA)TWy + Bh - m。 (1-2-28)迭代公式mk+1 = mk + (WA)tWA + B-1(WA)TWy + Bh - mk (1-2-29)如果先验信息有疑义(或不可信,)那就需要将约束置为，即h=。，。，0t，而且所 有B的元素均置为相等的常数(VBV1),这样所有的参数都具有相等的权重。在这种情 况下，B可以方便地由一单值未定乘子B所取代。这样就有参数校正解X = (WA)t +p 21 -i(WA)TWy-p 2 mo(1-2-3o)其迭代公式mk+1 = mk + (WA)t +P 21 -i(WA)rWy 。2 mk (1-2_31)因为D=1,这里饥I用以控制求解的步长，而命有助于减小其向零矢也的位置，我 们可以将这种方法称为平滑度约束反演或最小偏置算法。3. 与标准方法的关系在偏置估计中，如果BT,那么上述所有迭代估计公式均会简化为改进了的加权经典 最小平方化式mk+1 = mk + (WA)rWA-i(WA)rWy(1-2-32)偏置估计方法的稳定性和有效性主要取决于B和。方程(1-2-30)与通常的阻尼最小平方或岭回归优化公式(1-2-33)x = (WA)tWA + pi -1(WA)rWy的不同之处在于一危m0项，唯一目的是要对参数增量的变化范围置一边界。我们可将约束反演问题定义为求Lagrange函数的极值问题：L = (Wd 一 Wf (m)T(Wd 一 Wf (m) + p(m 一 mo)(m 一 mo)( 1-2-34)要搜寻的是最佳拟合数据的起始模型的有界摄动。在方程(1-2-29)中以量E取代WtW，我们有：mk+1= mk + (AtEA + B-1ArEy + Bh 一mk(1-2-35)如果B可以统计地解释为先验参数协方差矩阵的逆，则上述方程即等效Jackson和 Matsuui的Bayes估计方法，并类似于Tarantola和 Valette的非线性算法。因此，应用简单 代数，我们事实上已经导出一种与基于具有先验数据的概率统计处理的数学上比较严谨的非 线性反演法相类似的方法，但是应该注意到，Tarantola和 Valette的里程碑方法中的反演理论 和先验信息的使用均与我们的方法不尽相同。我们的主要兴趣在于迫使最终解尽可能与那些 先验参数估计相一致，因此方程(1-2-35)右端最后一项不为0,因为在实际情况下，已知的 先验参数估计很少。在Tarantola和 Valette算法中，h即是实际起始模型m。在这种情况下， 正如Pous等人指出的那样，方程1-2-35)的最后一项在第一次迭代中应为。我们将h和 m0分开，这点和Jackson和Matsuura的作法大致相当。但在Jackson和Matsuura的算法中， 方程(1-2-35)括号中的量乘了一个适当选择的因子)(0<b<1>如此说来，我们的简单方法 或许更为通用。第三节地震勘探中的反演方法一. 地震反演的分类地震反演通常分成叠前和叠后反演两大类：叠前反演应用较少，较成熟的莉反演； 叠后反演的大量应用是波阻抗反演，这是当前地震资料处理的重要结果。从反演方法上可将 地震反演划分为基于波动理论的波动方程反演和基Robinson褶积模型反演两大类在实际 工作中主要是基于Robinson模型的反演。我们通常所说的地震资料波阻抗反演指的是基于 Robinson褶积模型的叠后地震资料反演，目前常见的有递推反演(Recursive inversion，稀疏 脉冲反演(Sparse-spike inversioi和I基于模型的反演Model-based inversion,)后面两种反演方法 通常称为宽带反演。递推反演包括道积分GLOG、VLOG、SEISLOG、块状反演、带限反演 和PIVT等；稀疏脉冲反演包括多种实现方法如L模方法、最小熵方法、最大似然方法等； 基于模型的反演也包括广义线性反(GLI) (Cooke，1983)、地震岩性模ft(SLIM) (Gelfand, 1984)、鲁棒的速度反演方法ROVIM) (Fabre，1989)、宽带约束反演(BCI) (Martinez 1988)、PARM和 Jason等。二. 递推反演方法1. 波阻抗递推公式对于两层介质，反射系数为：p V -PVR 2-1-p V +P V2 21 1pP2分别为上下界面的密度，V、v2为上下界面的速度。当地下为多层水平介质时，任意第个界面的反射系数为：(1-3-1)p V _pz -zR a44_+4i_i i+4ii p V +pV z + z i+1 i+1 i ii+1 i对应的波阻抗为:zi+1=z )i 1-Ri(1-3-2)递推公式:zn+1=z n 1+Rn 0 n1 - R n=0n(1-3-3)zn+1如果用经过特殊处理的地震剖面记录道n代替反射系数Rn,则上式可写成(1-3-4)这就是递推反演的基本公式。2. 低频补偿地震反射系数剖面的频带是有限的，它缺失的是高频成分和低频成分，对波阻抗反演而 言，缺失高频成分只影响分辨率，而缺失低频成分就失去了速度的直流分量及速度斜坡的信 息。这种波阻抗剖面通常称为相对波阻抗剖面或剩余波阻抗剖面，剩余波阻抗剖面只能反映 波阻抗的相对变化而不能反映波阻抗的真实情况，因此必须在剩余波阻抗剖面上，再加上合 理的低频成分，进行低频补偿。(1) 利用声波测井资料补偿低频这是最常用的方法。(2) 利用叠加速度补偿低频但叠加速度补偿因为实际三维速度场精度有限会出现低频缺口，造成声测井曲线的值 偏高和偏低的振荡。(3) 利用地质模型补偿低频这种方法比较费时。低频缺口在波阻抗反演中是常见的，有时也是比较严重的问题。所幸其横向上速度相对 变化通常是正确的，仍然能确定目的层段上有意义的岩性变化。3. 一个简单应用道积分该方法不做低频补偿，得到的是相对波阻抗。用连续时间函数表示1-3-4) 式ZQ) = Z(t )exp2tr(t)dt(1-3-5)0 p如果经反褶积处理后的地震遥(t)的脉冲宽度足够小，认为X(t)与反射函数r (t) 成比例r (t)8x (t)则可近似求出任一时刻t的近似波阻抗Z(t) = Z(t )expKtx(t)dt(1-3-6)0们K为比例系数。具体实现步骤为： 将地震记录振幅标定到反射系数数量级 计算积分道A(k) = Y xi i=0 将积分结果转换为波阻抗Z(k) = Z exp A(k) 对转换结果作带通滤波得地层相对波阻抗图1-1给出了具体的应用例子，处理资料为JH三维工区一剖面。三. 稀疏脉冲反演方法这种方法假设地下反射系数序列是由一系列大的反射系数叠加在服从高斯分布的小反 射系数背景上构成的，主要有:L1范数反褶积、最小熵方法、最大似然方法等。L1范数反褶积最早曲arrodale于1973年提出，后经Taylor1979年及Oldenburg198年的 研究，改进成为一种独特的反褶积方法，它的特点是对子波的各种相位特性都有较好的适应 性。常规脉冲反褶积及预测反褶积都要假定子波是最小相位的，并且反射系数是白噪。在这 两个条件下，反褶积的求解运算工作只有在最小二乘的意义下(与期输出波形均方根误差最 小)，才能得到一组Toeplitz方程组，才能用莱文森递推法快速求解反褶积因子误差的最小 平方就被称为其范数为。而L1范数是不做平方的判断，而用误差的绝对值之和作为标准故 称其范数为L1范数。1985年王承曙等又提出范数反褶积，即其判断的范数可以不是，也不p是1，而是一个任意的正整数。由于采用了L1范数，带来的好处是对子波的相位特性放宽了 限制，但是在计算中没有脉冲反褶积那样简单了。它一般是从线性规划的理论出发，求解一 组超定方程组的最优解。在求解过程中必须反复迭代，或者化为一组非线性方程组，用非线 性规划方法迭代求解。最小熵方法由R.A.Wigginst先提出，它以方差模为判断准则，信号的规则性达到最大， 熵为最小。Wiggins的方差模的定义如下：Z X 4V = R -(1-3-7)ari maxZ尤上i式中七是地震道数据在某个时窗内的第个数值。当波形很突出时，X：达到很大振幅 值，于是七:睥、达到极大值，此时认为效果达到最佳。此方法对子波的相位特性不做约束， 而且在一定程度上可以把混合相位子波向零相位靠拢。但该方法假设反射系数是稀疏的，只 有当有少数大的脉冲存在时效果才很好，所以在具有亮点强波的剖面上，往往得到较好的反 演效果。稀疏脉冲反演方法的输出为矩形波阻抗曲线形式，地层边界清晰，对厚层碳酸盐岩地区 较为合适。然而其致命的弱点是要求反射系数是稀疏的，而实际上大多数地震道的反射系数 是稠密的。四. 基于模型的反演1.流程框图模型为基础的方法，或简称模型法，首先构造一个地质模型，并将其与地震资料进行比 较，然后利用比较的结果，迭代地更新模型，直至其与地震资料资料吻合为止。其示意图如 图 1-2。图1-2基于模型的反演示意图这种方法避免了直接对地震资料进行反演，可以有较高的分辨率。然而，随之而来的是 多解性，很可能一个与地震资料吻合得很好的模型却是错的。尽管如此，由于有地震测井和 地质资料的约束，常常可以把多解性降低到最低限度，在储层横向预测和油藏描述中起重要 作用。基于模型的反演包括广义线性反演GLI) (Cooke，1983)、地震岩性模拟(SLIM)(Gelfand, 1984)、鲁棒的速度反演方法RO VIM) (Fabre，1989)、宽带约束反演(BCI) (Martinez 1988)、PARM和Jason等，代表着当前反演的主流和发展趋势。2 .应用实例以Strata软件为例。在软件中有三种反演方法：带限Bnd Limited）反演，Block反演以及稀疏脉冲Sparse Spike）反演。此处应用了带限反演方法，它是一种传统的递推反 演计算。这种方法把地震道看成是经过零相位子波滤波后一系列的反射系数。由于地震数据 中速度的低原成份已滤去，而经过该处理后会恢复，因此要加上一个低步郎艮制，将反射系数 序列转化为声阻抗。此过程是忽略子波影响的，所以声阻抗值很平滑。总之，这种方法是比 较简单的，只需要较短的计算时间，但由于不考虑子波影响，薄层很难分辨。图1-3显示的就是过井g247测线剖面作的地震反演。图1-3 一种基于模型的反演一一带限反演