波阻抗约束反演方法研究

1、波阻抗约束反演方法研究吕铁良1 刁永波2 王珊珊3（1中国石油大学(华东)，山东东营 257061）（2. 四川石油物探公司，四川成都，611430）（3.胜利石油管理局黄河钻井公司）摘要：本文对波阻抗约束反演中的各种约束反演进行了研究。依据数学算法的不同，反演方法可分为线性反演和非线性反演。依据这个分类标准，介绍了线性约束反演方法中的广义线性反演（GLI）、宽带约束反演（BCI）、地震资料的测井约束多项式拟合外推反演、测井构造约束地震资料目标反演；以及非线性反演方法中的模拟退火法(SA)、遗传算法(GA)和人工神经网络方法(ANN)。关键词：线性反演；非线性反演；约束条件；波阻抗Study

2、of constraint methodology in constrained impedance inversionLü Tieliang1 Diao Yongbo2 Wang Shanshan3(1. China University of Petroleum ( East China ), Dongying, Shandong 257061) ( 2. Sichuan Geophysics Company, Chengdo, Sichuan, 611430) ( 3. Shengli Yellow River Drilling Co. 257046)Abstract: In

3、this paper, several constrained inversion methods in the impedance inversion are studied. Based on the mathematics algorithm, the inversion method can be divided into linear inversion and nonlinear inversion, the linear inversion includes: generalized linear inversion (GLI), broad-band constrained i

4、nversion (BCI), log-constrained polynomial fitting extrapolation inversion of seismic data, target inversion of seismic data under the restrictions of logging and structural data the nonlinear inversion method includes: simulated annealing (SA), genetic algorithm (GA) and artificial neural networks

5、(ANN). Keywords: linear inversion; nonlinear inversion; constraint condition; impedance inversion前 言地球物理反演就是利用观测数据恢复地下地质结构和岩石性质的方法,狭义地说,反演就是从有限频带宽度的地震数据中恢复出宽带波阻抗。在寻找和发现薄的砂泥岩交互储集作者简介：吕铁良，1982年生，男，汉，山东德州人， 2005年在中国石油大学获得学士学位，现在中国石油大学攻读地球探测与信息技术硕士学位。主要从事波阻抗约束反演方法、微地震技术监测在油田开发中的应用等方面的学习与研究。 联系地址：山东省东营市石

6、油大学地球资源与信息学院地球物理系，257061；层过程中,用常规地震记录已显得力不从心。线性反演的理论框架是在19671970年期间，由美国地球物理学家Backus和应用数学家Gilbert奠定的，这一理论简称BG理论，该方法是建立在连续介质这一假设基础上的。由于地震反演问题通常都是强烈非线性的，目标函数的最优解，其可行性完全依赖于初值的选取，当反演参数的初值未知时，就无法判断反演结果的可靠性。因此，对于非线性的地球物理反演问题，常规的算法就受到很大的局限，必须采用基于全局最优的算法才能有效地解决反演问题，提高反演结果的可靠性。然而，在实际的应用过程中，我们在追求非线性全局寻优算法的同时，并

7、不完全放弃对常规线性方法的使用。下面将对线性反演方法和非线性反演方法做一个较为详细的介绍。线性反演方法1广义线性反演方法的基本原理现行的波阻抗约束反演方法,是建立在D. A. Cooke 的广义线性反演( GLI)方法理论基9础之上,以褶积模型的泰勒级数展开为基础的波阻抗反演方法。我们知道在忽略噪声条件下的自激自收反射地震道s(t)可以写为10s(t)=w(t)rj(tTj)j=1n1 (2-1) 式中,w(t)地震子波,(t)为单位脉冲函数,rj为对应Tj时刻的一次波反射振幅。对(21) 式离散取样,令t=it,Tj=jt,t为地震记录的采样间隔。由(21)式不难求出它的偏导数矩阵A 的元素

8、为ai,j=2zj1s(it)w(it)i(j1)=t2zj(zj+zj+1) (2-2) 2zj+1(zj+zj+1)2w(it)(ij)t于是,地震道反演的最小二乘算法可归结为求Tzk=Ak1Ak1+kI()1TAk,2,L (2-3) 1(ssk1) k=1,2,L (2-4) k=1式中,k为迭代次数,zk为第k次迭代波阻抗向量zk的修改向量,sk为由zk得来的合zk=zk+1+µzk 成地震道,s为实际地震道,k为阻尼因子,µ为松驰因子,满足0<k<1, 0<µ10。我们采用多维求根的牛顿法求解方程组(2-3)和(2-4) 。2宽带约束

9、反演宽带约束反演是从确定一个初始波阻抗模型开始的(图2l)；模型被参数化为反射系2数和延迟时间，进而形成地震道估计，并与实际地震道比较产生剩余误差道；利用误差道来修正模型参数，直到满意为止。从BCI的实现过程分析，可看出它比合成声波测井等常规方法在理论上有许多优点：第一，常规方法仅仅从地震资料入手，求得的波阻抗只是有限带宽的，缺少真实的低频成分，这不利于进一步求取岩性和孔隙度等参数。然而，BCI求取的是宽频带波阻抗模型(或绝对速度)，可直接用来求岩性等参数。这主要是因为，由测井资料读取的波阻抗曲线是宽频带的，具有真实的低频信息。第二，BCI求取的宽频带波阻抗模型是一个优化的波阻抗模型，它显示出

10、垂向分辨率高于由有限带宽地震数据提供的分辨率。只要给定的初始模型具有较高的垂向分辨率(与采样率也有关)，经过迭代后，就有可能获得一个最终的高分辨率模型。图21 宽带约束反演流程图 Flow chart of the broad width constrained inversion 3地震资料的测井约束多项式拟合外推反演 将具有高垂向分辨率的测井资料作为约束条件，与具有较高横向分辨率的地震资料一起反演，再根据地层变化情况进行横向外推反演计算，这种约束反演方法能大大地减轻反演解的“不适定”程度。本方法采用最小二乘算法，在井旁进行测井资料约束地震反演，在横向上根据多项式拟合方法拟合的相位时间变化的

11、方向，进行横向拟合外推反演。 多项式拟合外推约束反演方法 对于等权约束反演公式和非等权约束反演公式为测井约束单道反演公式，也就是一维约束反演模型。若对多道地震记录进行约束反演，就应在此基础上进行横向外推约束反演。3若地层是水平的，平推出去即可，当地层倾斜或弯曲时，平推会使斜层分段拉平，这是不合理的。俞寿朋先生等（1988）提出的用多项式拟合地震信号相位时间的方法，提供了沿地层产状合理外推的良好方法。该法假设信号在各道上出现时间符合一个待定的多项式，而不采用时间相同或线性变化的假设。对于叠加剖面上在时空域划分出的每一个窗口，用正交多项式来表示窗口的中点时间，所导出的03次多项式p0(x)p3(x

12、)的表达式为p0(x)=1p1(x)=x (2-5a) (2-5b)1p2(x)=x2N(N+1) (2-5c) 31 (2-5d) p3(x)=x3+(3N2+3N1)x 5其中窗口的大小在空间方向为2N+1道，x为相对道序号。在窗口范围内，x=N至x=N。注意到上面导出的p2(x)表达式，与文献6中给出的p2(x)表达式是不同的。当窗口内的道数给定时，正交多项式是固定的。通过对窗口内数据计算归一化多道互相关值，可扫描出窗口内信号的期望相位时间多项式系数，由此确定了信号窗口的形状，即约束条件的外推方向。采用多项式相位时间拟合方法进行外推约束反演的具体做法是从井开始，将由测井资料得到的先验模型

13、参数向量或上道反演出的模型参数向量，沿多项式拟合出的相位变化方向进行外推，参与下一地震道的约束反演，从而实现测井资料约束横向拟合外推反演。2.4测井-构造约束地震资料目标反演对叠后地震剖面进行构造解释后，再从井旁道出发对目标层段进行非等时的外推波阻抗反演，使用相同层位的地震数据能够在相同低频背景的约束条件下进行反演运算。这对于提高解的精度、减少反演的多解性是十分有益的17。反演方法原理假设反射系数序列模型可表示为r(t)=a(i)tit i=1,2,L,L (2-6)i=1L式中：r(t)是有L个反射界面的反射系数序列模型；a(i)是第i个反射界面的反射系数；t是采样间隔。地震道用下式来表示s

14、(t)=r(t)w(t)+n(t) t=1,2,L,M (2-7)式中：s(t)是M个时间取样值的地震道；w(t)是频带有限的地震子波；n(t)是频带有限的附加干扰。利用泰勒级数将式（2-7）在给定的波阻抗初始值附近展开S=S0+G(PP0)+N (2-8)式中：S是实际地震道矢量；S0是由初始波阻抗正演计算得到的合成地震道矢量；P是波阻抗矢量；P0是给定的初始波阻抗矢量；N是附加的干扰矢量；G是雅可比矩阵，其中每个元素为gi,j=2zj1(zj+zj+1)2w(it)i(j1)t2zj+1(zj+zj+1)2w(it)(ij)t式中：Zj为第j层的波阻抗；i=1,2,L,L;j=1,2,L,

15、M。将式（2-42）简化为S=GP+N式中：S=SS0,P=PP0。 (2-9)利用一套有约束的线性反演方法，把给定的初始波阻抗值与地震资料结合起来，反演采用随机算法。用此算法使初始值P0产生摄动，以得出一个经修改的估算值P。P的估算值由下式给出T1TP=P0+GG+CCGS NP1 (2-10)式中：CN是噪声的协方差矩阵；CP是P的协方差矩阵。在测井资料较多的情况下，为提高反演解的精度和分辨率，对P与观测数据间的最大似然解估计做一些修改，得出下列关系式1TP=P0+GG+CNC(GS+CNCPm0) T11P (2-11)其中m0为先验信息。非线性反演方法近年来，非线性理论在自然科学的各个

16、领域成为研究前沿，受到人们的普遍关注。大多数自然现象都是非线性的，地球物理反问题也是这样。 实践表明，非线性反演要比线性反演更接近实际。非线性优化方法有基于导数的最陡下降法、牛顿法、共轭梯度法等,基于非导数的非线性优化方法包括遗传算法(GA)、模拟退火法(SA)、人工神经网络方法(ANN)、随机搜索、下山单纯形搜索等等。1模拟退火算法模拟退火法(Simulated Annealing，简称SA)是Kirkpatrick等1983年首先提出的，自问世以来很快受到广大地球物理学家的密切注意。近年来，它在可靠性和有效性等方面都得到了很大的发展，成为一种十分受欢迎的非线性多参数联合反演方法（1）、模拟

17、退火算法的物理描述模拟退火法是一种启发式蒙特卡洛反演方法。它模拟退火的物理过程：物质先被熔化，11。然后逐渐冷却。在冷却过程中，有可能产生非晶体状的亚稳态玻璃体，也有可能产生稳态的晶体。晶体相应于该物理系统能量最小的基本状态；玻璃体相应于其能量达到次极小的亚稳态。把物理系统的能量模拟成反演问题的目标函数；把晶体的生成模拟成搜索到目标函数的整体极值；把玻璃体的形成模拟成错误地搜索到局部极值，就形成能有效地求解非线性反演问题、得到相应于整体极值的某种意义下的“最佳”解的模拟退火法。（2）、退火过程的统计描述统计物理从物质是由大量微观粒子组成这一事实出发，认为物质的宏观性质是大量微观粒子热运动的平均

18、结果，宏观量是微观量的统计平均。因此，它研究的是由大量微观粒子组成的宏观系统的统计特性。统计物理学的基本研究成果是得到一个处于平衡状态下的宏观系统的统计分布，其中最重要的一种分布是吉布斯(Gibbs)分布。它认为系统处于某一种状态x下的概率由下式确定 P(x)=E(x)1exp (3-1) ZKBT式中E(x)表示系统处于状态x下的能量，KB为玻尔兹曼(Boltzman)常数，T为绝对温度，Z为正规化常数E(x)Z=exp (3-2) KTxB对于平衡的系统，吉布斯分布函数描述了系统状态的期望扰动。这种扰动既可能增加系统能量，也可能减少能量。向增加能量方向扰动的可能性大还是向减少能量方向扰动的

19、可能性大由吉布斯分布决定。应当注意的是，绝对温度T这一参数对吉布斯分布的影响很大，因而对期望扰动的影响很大。当系统温度T很高时，使系统能量增加的扰动与使系统能量减少的扰动都有差不多的可能性。但是，当系统温度变小时，吉布斯分布逐渐给低能量状态以较大的概率。在极限情况T0时，吉布斯分布只允许向能量减小的方向扰动，系统进入基态。基态相当于最规则的晶体状态。但是，为了达到基态，系统必须慢慢冷却。因为若冷却太快则可能形成非晶体状的亚稳态玻璃体。将晶体形成模拟成搜索到整体极值，玻璃体形成模拟成搜索到局部极值，就可以利用吉布斯分布指导随机搜索。搜索时缓慢降温，保持平衡，就可以求出整体极值解。2遗传算法遗传算

20、法，由于它不需要计算目标函数的导数，是一种全局优化算法，其非线性优化的实现机制完全是随机的(统计的)，不存在任何确定性的规律来指导参数空间的搜索。在遗传算法中，大多数模型及复杂的模型用简单的二进制字符串来表示，这些字符串的变化方式用生物系统从遗传到繁殖的演化来类比。在找到拟合更好的模型时，代表模型的字符串中某些优生的字节被保留，以使参数空间的搜索效率大大提高。因此，遗传算法是在地6球物理反演中应用得较多的一种。1、遗传算法概述与模拟退火法相似，遗传算法(Genetic A1gorithm，简称GA)也是一种启发式蒙特卡洛反演方法，即有指导地而不是盲目地随机搜索的方法。它可以解决复杂的大尺度、多

21、变量和非线性反演问题。与模拟退火法模拟物理系统的结晶过程不同，遗传算法基于生物系统的自然选择原理和自然遗传机制，模拟自然界中的生命进化过程，在人工系统中解决复杂的、特定目标的非线性反演问题。遗传算法用经二进制编码后的模型参数进行工作。模型参数经二进制编码后组成一个“串”，类似于生物遗传中的关键物质“染色体”。遗传算法模拟生物遗传中染色体遗传基因的变化来改变模型参数。图3-3 遗传算法基本流程图2、遗传算法的实现方法（1）编码和初始群体形成首先要确定参数的变化范围及所要求的分辨率。通常遗传算法对模型参数的二进制编码进行工作，所以遗传算法的首要步骤是对模型参数进行二进制编码。按照遗传学的术语，这个

22、二进制的“染色体”中的每一个二进制位称为一个“基因”，只能取0或1两个值。需要注意的是，参数编码并不意味着从十进制换算成二进制这种简单的变换。同一个参数7可以编成不同的码，取决于参数的取值范围(最大值max，最小值min)和要求的精度。如果该染色体有n位的话，则min为该染色体n位全部为0所表示的参数值，而max=min+2(n1)，其余值可类推得到。（2）选择和再生选择是产生新的模型群体的过程中的第一步。它从群体中挑选模型配成对(亲本)以进行交换。选择的基本思想为群体中的每个成员都有合理的繁殖机会，但较优秀的成员(即适值较大的成员)应有更多的机会。因此，选择是建立在群体中各模型适值大小的基础

23、上进行的，是按适值算出的选择概率进行随机采样得到的。适值的大小与模型对应的目标函数(或后验概率)密切相关。若求极小，则目标函数(或后验概率)越小适值越大；反之，若求极大，则目标函数(或后验概率)越大适值越小。按适值计算选择概率可以有多种方法(只要能保证适值越大，概率越大即可)。一种最常用也最简单的选择概率计算公式为PS(xi)=f(xi)式中f(xi)为模型xi的适值。（3）基因交换一旦亲本选择出来之后，就可以进入交换这一步。交换是遗传算法的“繁殖”过程，是遗传算法的内在力。交换为亲本模型的重组，即将两个亲本模型拷贝的片段剪接在一起构成后代子本模型。显然，这种交换完全模拟遗传过程中两个染色体遗

24、传基因的交换过程。最基本也最简单的交换方式为一点交换。在染色体内部随机地选择一个交换点，将一个亲本染色体在此点前的第一段与另一个亲本染色体在此点后的第二段结合在一起构成子本后代染色体，从而得到两个子本后代。例如，对于某一成员A1和A2： f(x) (3-3) iiA1 = 0 1| 1 0 1A2 = 1 1| 0 0 0随机选择的交换位置如竖线所示，则基因交换如图3-4所示： 基因交换所产生的下一代两个新成员为A1和A2A1= 0 1 0 0 0 8= 1 1 1 0 1 A2因为交换的实质是在模型空间中进行大范围的搜索，搜索的空间区域很可能与原先的采样区域相距较远，故这种搜索属于非邻近区域搜索过程。它可以产生一个十分有效的模型空间普查。线性化或拟线性反演方法、模拟退火法都属于邻近区域搜索过程，搜索的强度和内在的潜能远不如遗传算法。（4）基因变异变异是对偶然的(按较低的变异概率随机选择的)后代中的一个或多个随机选择的基因作随机摄动。变异在遗传过程中是十分重要的。因为若不存在变异，则子本模型不可能获得群体中不存在的染色体基因，因而也就不可能出现强有力的进化，出现超过前代的变化。最简单的变异方法就是将模型参数二进制编码的某一位由1变为0或由0变为1。发生变异的情况应当少(即变异概率Pm应当低)，但必须有。生物学中变异是保证物种不会退化的重要手段。在非线性反演中，变异是使得模