基于迭代再加权最小二乘的地震资料稀疏反演方法—硕士学位论文.doc

中文摘要 iii 中文摘要 地震数据反演的目的在于通过求解反射系数序列来反映地下介 质的分布规律，这就需从地震数据中可靠地消除地震子波的影响。 然而受到有限带宽的影响，往往缺失反射系数序列中的低频和高频 信息；同时由于受到不可避免的噪声的影响，反演得到的结果分辨 率通常比较低。 本论文的目的在于使用迭代再加权的寻优方法由地震数据得到 高分辨率的反射系数序列；该方法通过再加权的方法将稀疏优化问 题转化为一系列权不断变化的线性方程组，在保证结果精度的前提 下增加了求解的速度。对寻优过程中的线性方程组求解采用针对性 的预条件化，并在权重的计算中引入伸缩变换，优化了整个寻优过 程的稳定性和求解速度。经过大量的数值试验，得出了合理的参数 选择范围，通过与真值的比较以及实际资料处理验证了本方法的实 用性。 关键词：波阻抗反演，迭代再加权，gauss-siedel 迭代，共轭梯度 法，预条件化，伸缩变换 abstract iv abstract the basic goal impedance inversion is to recover the layered structure of the surface as a function of depth from observed seismic data, which usually expressed as series of reflection coefficients. such a goal could be achieved when the effects of seismic wavelet is eliminated in an effective way. however, it is hard to get high resolution reflection coefficients series because of the bandlimited effect of seismic data with the resulting loss of low and high frequency information about subsurface structure. the unavoidable noise reduced the resolution,too. the purpose of this dissertation is to invert a high resolution reflection coefficients series from seismic data using iterative reweighted technique. this technique transform the sparse optimization into a series of wighted linear function, speeds up the computing while keep signal- to-noise rate high. in order to stabilize and accelerate the solving,pre-conditioning according to the linear function is used and shrink operator is introduced in computing the weighing. through lots of numerical experimentations, the range of parameters is settled, and the practicability is conformed together with the real data processing. 中文摘要 iii keywords: impedance inversion, iterative reweighted, gauss-siedel iteration,conjugate gradient method, pre-conditioning, shrink 目 录 v 目 录 第一章 引 言 1 1.1 选题依据和研究意义 .1 1.2 波阻抗反演研究现状 .1 1.3 稀疏反演方法的研究情况 .3 1.4 论文结构与主要内容 .4 1.5 论文的主要创新点 .4 第二章 波阻抗与稀疏反演概述 .5 2.1 基本概念 .5 2.1.1 常规波阻抗反演技术的基本假设前提 .5 2.1.2 波阻抗反演方法的基本原理 .5 2.1.3 波阻抗反演的分类 10 2.2 稀疏反演方法 11 2.2.1 匹配追踪 11 2.2.2 稀疏反演和压缩感知 12 第三章 迭代再加权最小二乘地震资料稀疏反演方法 17 3.1 迭代再加权方法原理 17 3.2 反射系数稀疏反演的迭代再加权方法 20 3.3 解线性方程组的迭代方法 21 3.3.1 gauss-seidel 迭代方法 .21 3.3.2 共轭梯度法 22 3.4 收缩变换 26 3.5 预条件化 28 3.5 流程图 30 基于迭代再加权最小二乘的地震资料稀疏反演方法 vi 第四章 数值模拟和实际资料 31 4.1 不同噪声水平下的原始地震数据 31 4.2 不同主频子波反演结果 32 4.3 噪音水平不同时反演结果 34 4.4 权重因子，噪声水平不同时反演结果 37 4.5 范数逼近程度不同时反演结果 39 4.6 实际资料处理 40 第五章 结论与认 识 41 参考文献 .43 作者简介 .47 攻读硕士学位期间已发表及待发表的论文 .48 攻读硕士学位期间参加的科研项目 .48 第一章 引 言 1 第一章 引 言 1.1 选题依据和研究意义 油气资源关系民生和国家安全，与国民经济的大部分产业有直 接或间接的联系，对经济发展有着十分轻重的影响，也是今后相当 长时间内难以替代的主要能源。所以加强国内油气的勘探能力具有 十分重要的意义。 地震、测井、钻井等方法是石油工作者了解地下地质构造、地 层、岩性、物性、含油气性等重要的信息来源。虽然地震方法获得 资料的分辨率远远不及测井、钻井，但是随着地震勘探技术进一步 地发展，地震资料的信噪比、分辨率、成像准确性都获得了很大的 提升。同时地震资料包含大量地下地质信息，且覆盖面积广，具有 三维特性，所以地震勘探愈来愈受到重视。 地震勘探的主要原理是从地表某一点位激发地震子波，子波在 向下传播的过程中遇到波阻抗界面后会发生反射，传回一部分能量， 被传感器接收。随着这一激发、传播、反射、接收的过程不断进行， 传感器记录下一系列按时间排列的地震波，即地震记录。地震勘探 的基础是地下不同地层存在波阻抗差异，因为地震波只有传播到有 波阻抗差异的地层分界面时才会发生反射。 地震资料反演其基本目的是充分利用测井、钻井、地质资料提 供丰富的构造、层位、岩性等信息，由常规的地震剖面推导出地下 地层的波阻抗、密度、速度等信息，为勘探开发提供重要的依据。 而石油勘探开发的形式在不断地发展，计算机技术也不断地提高， 使得地震资料反演技术成为了当前热门的研究课题。其中波阻抗反 基于迭代再加权最小二乘的地震资料稀疏反演方法 2 演（acoustic impedance，ai）是反演技术中很重要的一个部分， 它利用地震资料反演反射系数、进而得到地层波阻抗信息的技术。 同时，地震资料反演技术本身具有的不适定性，应用条件的限 制以及储层的复杂性，常规叠后地震资料品质不高且分辨率低等因 素，一定程度上影响了该技术更广泛的应用。 1.2 波阻抗反演研究现状 地震反演的基本目的是利用地震波在地下介质中的传播规律， 通过数据采集、处理与解释等流程，推测地下岩层结构和物性参数 的空间分布，为勘探开发提供重要的依据 13。波阻抗反演是地震反 演的一个重要组成部分，由于波阻抗信息是联系地质和地球物理的 一座桥梁，叠后计算数据量相对较少，在实际生产中应用方便而且 效果明显。因此，波阻抗反演在地震反演中具有特殊的地位。侠义 的地震反演概念指的就是波阻抗反演 4。 波阻抗反演是利用地震资料反演地层波阻抗（或速度）的地震 特殊处理解释技术，它产生于 20 世纪 70 年代，80 年代开始得到蓬 勃发展。1983 年，cooke 介绍了地震资料广义线性反演方法，从而 揭开了波阻抗反演技术的新篇章 5。周竹生等人在 90 年代初期提出 了综合利用地质、地震和测井资料进行约束反演，克服单一的线性 反演方法的缺陷 6。90 年代中期，李宏兵在国内提出了一条走递推 反演与宽带约束反演结合的道路，此方法为我们解决了从单道出发 的反演方法不能在跟不上消除噪音的困惑 7。在此基础之上，有人 进行了无井多道反演和有井多道反演的研究 8,9，使波阻抗反演方法 更加完善。90 年代至今，围绕一维波阻抗反演的各类算法及应用成 果层出不穷，随着研究的开展，在 1997 年左右开始出现了一些反思 第一章 引 言 3 的文章，指出了波阻抗反演中存在的一些缺陷，并提出了一些解决 方案 1012。bp amoco 公司的 connolly 在 1999 年正式发表了弹性波 阻抗反演方法的论文，该方法最早是在 19931994 年由 bp amoco 公 司发展起来，随后在 2000 年的 seg 年会上同时出现了四篇论文对 ei(elastic impedance)进行了研究，arco 公司介绍了他们申请专 利的弹性波阻抗反演方法，认为在求取的反射系数的稳定性方面要 好于 connolly 方法，而且计算的 ei 和 ai 数值在一个尺度下，同时 bp amoco 公司在会上提出了扩充弹性波阻抗方法，可以用于流体和 岩性的预测。此外，paradigm 公司在其商业软件 vanguard 的最新 版块中也有有关弹性波阻抗反演功能的出现。jason 公司则推出了 rock trace 弹性反演模块，以纵波波阻抗和横波波阻抗的概念来区 别其以往软件中的波阻抗概念。这些进展说明弹性波阻抗已经成为 波阻抗反演进一步发展的方向之一，地震反演的发展正走向 ai 和 ei 结合、ai 和 avo 结合的道路 13。这个发展过程大致上经历了从 直接反演到模型反演、从叠后反演到叠前反演、从线性反演到非线 性反演发展的过程 16,33。近十几年来，非线性反演方法得到了有效 利用。实践表明，非线性反演比线性反演更接近真实。非线性优化 方法有基于导数的最陡下降法、牛顿法、共轭梯度法等；基于非导 数的神经网络方法、遗传算法（ga） 、模拟退火方法（sa） 、随机搜 索等 17。为了提高反演的精度和收敛速度，针对非线性反演方法， 地球物理工作者将各种最优化方法引入到了波阻抗反演中，各种各 样的混合优化算法也陆续出现，力求在解决计算精度的同时提高算 法收敛速度，减少反演多解性 18。尽管已经作了大量的研究工作， 波阻抗反演方法仍然存在许多的不足，如算法复杂、实现困难、抗 干扰能力差(从而导致分辨率降低)、收敛速度慢等，因此仍需要进 基于迭代再加权最小二乘的地震资料稀疏反演方法 4 一步研究和改进这些优化反演算法。在理论指导下，国内外相继研 制出一些商业化程度高的波阻抗反演方法及软件，如 stratr、glog、jason 等 19。目前国内使用的波阻抗反演软件基本 上都是引进软件，国产软件较少。96 年以前的软件一般都是单井控 制下的地震资料反演，只能解决构造简单没有断层或岩性、岩相横 向变化小的问题，在多井条件下和地质条件复杂的地区难以推广使 用。国外从 80 年代末提出了利用声波测井资料作为约束条件，对过 井地震剖面进行正反演联合迭代求取地下波阻抗的方法，将反演方 法推向非线性，这种新方法利用了测井资料的高频信息，因而表现 了强劲的发展势头，目前成为各国外石油公司重点发展的对象，如 美国 hgs 公司宽带约束反演 bci 技术，法国 cgg 公司波阻抗反演模 拟 rovim 技术，俄罗斯的 parm,加拿大的 strata，丹麦的 isis，以 及荷兰推出的 jason 技术。今后地震反演的技术，是在现在的基础 上，充分利用测井、地震、地质的信息，继续探索先进的算法，使 储层横向预测更符合实际。这些方法在实际中取得了一定的应用效 果，但是还不能够完全满足实际需要。不同反演方法结果差别较大， 不同人使用同一种方法结果也有差别，往往导致错误的地质解释。 究其原因，除了地震资料品质、子波提取不准确、垂直入射假设与 实际有误差等因素外，还有另外两个方面的主要影响因素：其一， 由于地震数据具有带限性质，直接反演只能获得波阻抗模型中的中 频成分，缺少应有的低频和高频成分。其二，hardmard 意义下，波 阻抗反演是一个不适定问题，其反演结果会出现不稳定性和多解性。 目前，解决带限问题最有效的手段是对解加以约束进而进行反演 （如测井数据的约束） 。而解决不适定问题主要使用正则化方法并辅 之以适当的最优化技巧 15。通常利用测井资料在地质知识约束的插 第一章 引 言 5 值方法，来弥补地震资料中仍缺乏的低频和高频信息。 1.3 稀疏反演方法的研究情况 线性问题的稀疏解在很多领域都有研究。线性问题可以分为正 定问题，欠定问题和超定问题。正定问题的解存在并且唯一。欠定 问题存在无穷多解，因此一般研究其最简单的解，即稀疏解。超定 问题不存在严格的解，一般求出满足最小二乘意义的近似解，当然 也可先验约束研究其稀疏近似解。目前的稀疏算法主要针对欠定的 线性问题求解。 在实际应用中求稀疏解的最常用方法是匹配追踪类方法（davis et al，1997；donoho et al，2006；tropp ，2004） ，主要有匹配追 踪，正交匹配追踪和其他改进算法。这类方法的执行简单，但是计 算量巨大导致计算速度慢，尤其在解不是严格稀疏和数据存在误差 的情况下该问题尤其突出。 稀疏解问题可以描述为一个目标函数为 l1范数最小化的约束优 化问题，被称为基追踪问题（chen et al，1998） 。这样就可以用常 见的凸优化算法求解。基追踪问题可以变为线性规划和二次规划问 题，然后采用标准的内点法求解（chen et al，1998；kim et al，2007；wang et al，2009b） 。 1.4 论文结构与主要内容 本文针对基于迭代加权的反演问题，系统的介绍了迭代加权的 波阻抗反演方法，具体内容有： （1）对波阻抗反演方法涉及到的基本概念、原理进行阐述，并 介绍了几种常用的反演方法； （2）从介绍波阻抗反演的反褶积分析入手，分析并导出了基于 迭代加权的波阻抗反演的基本问题，并介绍了迭代方法以及迭代的 基于迭代再加权最小二乘的地震资料稀疏反演方法 6 效果。 1.5 论文的主要创新点 本文对迭代波阻抗反演的反褶积问题进行了进一步的分析，从 而导出了基于迭代再加权的波阻抗反演方法，并且应用了几种迭代 方法来解方程组，通过理论模拟和实际资料处理证明了新方法在收 敛速度上的优越性。 (1) 对波阻抗反演方法进行了系统分类，对稀疏反演方法以及 与之相关的匹配追踪和压缩反演进行了论述； (2) 从反褶积问题入手，建立了一种基于迭代加权的波阻抗反 演方法，分析了加权因子的选取，并研究了 gauss-seidel 迭代法、 共轭梯度法和预条件加速技巧； (3) 对以上方法，用不同噪音水平的地震数据进行了反演计算， 计算结果验证了算法的有效性。 第二章 波阻抗与稀疏反演概述 7 第二章 波阻抗与稀疏反演概述 2.1 基本概念 2.1.1 常规波阻抗反演技术的基本假设前提 目前常用的波阻抗反演软件所采用的方法基本上都是在褶积模 型的基础上建立的，因此，要求资料都要满足褶积模型的假设前提， 可以概括为以下四个方面。 （1）地震模型：假设地层是水平层状介质、地震波为平面波法向入 射，其地震剖面为正入射剖面，并假设地震道为地震子波与地层反 射系数的褶积。 （2）反射系数序列：普通递推反演中，假设反射系数为完全随机的 序列；在约束稀疏脉冲反演中，假设反射系数为由一系列大的反射 系数叠加在高斯分布的小反射系数的背景上构成。 （3）噪音分量：通常假设波阻抗反演输入的地震数据的振幅信息反 映了地下波阻抗变化的情况，地震剖面没有多次波和绕射波等噪音 分量。因此在资料处理时，可以考虑的处理流程是反褶积、噪音剔 正，尤其是多次波，处理的最终目标是得到高分辨率的保幅波场。 （4）地震子波：假设反射系数剖面中的每一道都可以看作地下反射 系数与一个零相位子波的褶积。实践情况下地震子波往往是混合相 位或最小相位的，因此需要对地震剖面进行相位校正处理。 2.1.2 波阻抗反演方法的基本原理 波阻抗波在某种介质中传播的速度与介质密度的乘积，在本 论文中其符号和定义为 ， z 表示波阻抗， 表示密度， v 表示v 基于迭代再加权最小二乘的地震资料稀疏反演方法 8 波速。根据弹性波理论，当一个纵波入射到波阻抗不同的两种介质 的分界面时，一般要产生反射纵波、反射横波、透过纵波、透过横 波四种新的波动。一般主要研究纵波反射波；在水平层状介质、地 震波为平面波法向入射的假设下，问题会得到简化，只考虑产生的 反射纵波和透射纵波、反射横波和透射横波的情况。 地震子波震源附近地震波以冲击波的形式传播，当传播到一 定距离时，波形逐渐稳定，此时的地震波被称为地震子波。进行反 演时，求取地震子波的方法 一般有统计法和参数法。统计法是从实际地震记录中统计出地震子 波的振幅谱和相位谱，从而生成地震子波；参数法是给出某种地震 子波类型的参数而生成地震子波（如 ricker）子波。ricker 子波是 最常用的零相位子波，在时间域表示为 （2-1-1） 2)(21()tvmettf 其中， 是主频。零相位子波分辨力较高，最有利于解释。mv 由定义可见，只有层状介质的波阻抗存在差异时反射系数才不 为 0。当子波遇到波阻抗界面后，会发生反射，从而被接收到。由 于是按照时间序列排列的，以下分别记 为接收到的地震)(,)(trwtd 数据、子波、反射系数按照时间的序列，并已经过离散化，则按照 定义，可以得到： )1().)1(2)(1)( ltrwtrtrwtb l 表示子波的长度。数学上定义此运算为褶积；下记“*”为褶积运 第二章 波阻抗与稀疏反演概述 9 算符，则上式可以表示为 （2-1-2）rwb* 因此，地震记录实际上是由反射系数与子波的褶积得到的。如 果将子波作为输入信号，地下介质作为一个线性时不变的“系统” （或者滤波器） ，接收信号为信号通过系统后的输出，则该褶积可以 看做是一个系统的响应过程（陆基孟，1993） 。该响应不但包括一次 反射还包括多次放射。假定震源波形在向地下传播的过程中不变， 因此褶积模型不包括固有的吸收。地震记录可以看作是由震源子波 与地下反射率函数、多次反射、仪器等诸多因素的褶积的过程。虽 然可以给合成地震记录加上随机噪音，但是这里的用来建立反褶积 滤波器的褶积模型暂时不包括随机噪音。 图 2-1 模拟的反射系数、子波和合成的地震记录 如果能求得一个算子 是子波 的逆，即： 。则)(tb)(tw)(1twty 基于迭代再加权最小二乘的地震资料稀疏反演方法 10 它与地震记录 褶积时，可由（2.1）式得)(td )()()(1 trtrtwby （2-1-3） 即地震记录 和子波的逆进行褶积可以得到反射系数，这个过程)(t 叫做反褶积。 由以上过程可以看到，反射系数序列中含有波阻抗随时间变化的 信息，这就提供了速度和密度随时间变化的信息，因此可得到地层、 岩性及构造在地下空间分布的信息，在有利条件下还可得到岩石孔 隙率、渗透率、孔隙流体性质（油、气、水）乃至地层压力的信息。 波阻抗反演需要求出反射系数，因而其反演就是已知地震记录 ，)(tb 由提取的子波 反褶积得到反射系数序列 。)(tw)(tr 在使用炸药作为震源时，震源实际是一个 脉冲，由于地层的 滤波作用，震源信号在地下传播过程中高频被衰减，能量被吸收， 变成一个子波，所以在地震勘探中观测的信号其峰值频率往往小于 100hz，这就导致了记录信号的分辨率降低，记录到的反射是减弱了 强度的宽脉冲，并且通常随着反射时间的增加信号越宽和越弱。子 波在地层中传播，携带着反射系数序列这种有用的地质信息返回地 面，只有通过反褶积消除子波才能恢复反射系数序列；同时反褶积 可以压缩地震子波，压制交混回响和短周期多次波，从而提高时间 分辨率。反褶积可应用于叠前资料，也可应用于叠后资料。理想的 反褶积应该压缩子波并消除多次波，在地震道内只留下地层反射系 数。子波压缩可以通过将反滤波器作为反褶积算子来实现，反滤波 器可以将地震子波转变成尖脉冲，每个脉冲的强弱与界面的反射系 数的大小成正比，脉冲的极性反映界面反射系数的符号。当应用于 地震合成记录时，反滤波输出应为地层脉冲响应。 第二章 波阻抗与稀疏反演概述 11 反射系数根据弹性波理论，当一个纵波入射到波阻抗不同的 两种介质分界面时，一般要产生反射纵波、反射横波、透过纵波、 透过横波四种新的波动。在实际生产工作中，主要研究纵波反射波。 界面倾角不大时，可近似认为是垂直入射。在垂直入射的情况下， 将会使问题简化，纵波入射时我们将只考虑产生的反射纵波和透射 纵波、反射横波和透射横波的情况。这时界面的反射系数定义为 （2-1-4）1212vzarf 其中 af， ar分别表示反射振幅和入射振幅。而 ； ； 分21,z21,21,v 别对应上层和下层的波阻抗、密度、速度。当平面波非垂直入射时， 反射系数将随着入射角的变化而变化，这种变化与界面两边介质的 波阻抗和各种弹性参数有关，特别是与介质的泊松比有很密切的关 系（zoppritz 方程，入射角等于零是它的特例） 。 褶积模型根据介质模型假设不同，波阻抗反演分为基于波动 方程的反演和基于褶积模型的反演两大类。前者由于算法结构复杂、 计算量大、抗干扰差，很难获得一个稳定解，在实际中未得到广泛 应用。目前主要使用基于褶积模型的直接反演方法，它算法简单， 对地震噪音敏感性小，但一般情况下并不一定能得到一个稳定的解。 （1）褶积：在时域中为了便于求得线性时不变系统的零状态响应， 如果能够将任意激励信号分解为单元信号，并使每一单元信号在系 统中产生的零状态响应易于求得，那么根据系统的线性时不变特性， 就可以利用叠加原理很容易求得原信号在系统中产生的零状态响应。 系统对单位冲击函数 的零状态响应为 ，对于线性时不变系统，)(t)(th 当激励是强度为 的延时冲击函数 时，其零状态响应为k0k ，任意信号 可表示为)(0tkh)(te 基于迭代再加权最小二乘的地震资料稀疏反演方法 12 （2-1-5）dtetktketet )()()(lim)(0 作用下系统响应 为)(te)(tr (2-dthetkthetrkt )()(li)(0 1-6) 其中 为积分变量， 为参变量，上式褶积积分又可记作：t （2-1-7）)(*)(her （2）模型产生：水平层状介质模型是共中心点（cmp）叠加方法的 基础，在 cmp 叠加中，对分享同一炮检中点的一些地震道进行正常 时差校正求和，产生一个逼近 1d 层状介质垂直入射平面波响应的求 和道。本着这种数据，形成了一种地震道模型。假定从波阻抗不连 续界面反射（或透射）的子波与入射子波具有同样的波形，一个地 震道顺序记录这样一系列连续子波，也就记录了波阻抗不连续面。 地震道可以简单地看成是每一个独立的反射界面叠加的结果，这里 每个反射界面都认为是一个滤波器。地震记录可认为是震源子波与 介质反射系数（即地下垂直入射反射系数序列）的褶积。 1）时间域表达： 连续情形： (2-)()()()()( tedtrwtetrtwtb 1-8) 其中 为地震子波， 为地震反射序列， 为地震记录， 为随机噪r e 音。 矩阵-向量情形： ewrb （2-1-9） 第二章 波阻抗与稀疏反演概述 13 其中 为地震反射系数向量， 为噪音列向量，nrr,10e 为地震波记录向量， 为子波阵，子波长度为 ,),(bb wl llllnl wwww 21212121)1( 2）频率域表达： 连续情形： )()(erwb （2-1-10） 各表示傅立叶变化下的： -地震记录， -子波， -反射系数， -br 噪音。 2.1.3 波阻抗反演的分类 波阻抗信息是联系地质和地球物理的重要媒介，其叠后计算数 据量相对较少，在实际生产中应用方便、效果明显，因此通常的地 震反演概念指的就是波阻抗反演。将地震剖面转换成波阻抗剖面， 不仅便于解释人员将地震资料与测井资料连接对比，而且能有效地 对储层物性参数的变化进行研究，从而得到物性参数在空间上的分 布规律，对油气的勘探和开发有很重要的指导意义。对于波阻抗反 演的重要性，一个很经典的论述是著名的地球物理学家李庆忠院士 所说的“交到地质人员手中的地震资料应是作了反演的波阻抗剖面” 、 “波阻抗反演是高分辨率地震资料处理的最终表达形式” 25,26。而 基于迭代再加权最小二乘的地震资料稀疏反演方法 14 常常被人们引用的美国阿莫克石油公司在墨西哥湾海上寻找非背斜 油藏不断取得成功的例子 29也印证了这一说法。 波阻抗反演是目前应用最广泛的储层预测和油藏描述手段之一， 其发展日新月异，计算方法多种多样，这些方法由于其基于的数学 物理模型不同，各有优缺点，利用的基础资料和能够达到的反演效 果是不同的。波阻抗反演方法有多种 （seislog,velog,delog,pivt,bci,rovim,slim,parm 和 strata） 14, 概括起来有两大类：基于反射系数逆公式的直接反演和基于正演模 型的迭代反演。迭代反演又可以分为无井条件下的广义性反演 （gli）和有井条件下的宽带约束反演（bci） 33。还有一些属于非 线性反演范畴的方法，如混沌理论、神经网络等。 （1）依据算法的实现方式分为：基于反褶积的反演，包括道积 分、递归、广义线性、稀疏脉冲反演等；基于波动方程的反演，包 括 born 反散射反演等；基于随机过程的反演，包括随机反演、随机 模拟、模拟退火反演等；基于特征分析的反演，包括特征反演、神 经网络反演等；基于动力学特征的反演，包括混沌反演等 （2）按测井和地震的相对作用：无井反演，包括道积分等；地 震为主的测井约束类反演，包括递归反演、广义线性、宽带约束反 演、稀疏脉冲、模拟退火反演等；测井为主的统计类反演，包括随 机模拟、随机反演等；井-震约束联合反演，包括特征反演、神经网 络反演等。这种分类方案体现了测井资料对反演约束作用的强弱， 便于研究人员根据测井和地震基础资料评价各种反演方法的反演效 果。 第二章 波阻抗与稀疏反演概述 15 （3）依数学算法不同可分为：线性反演，包括道积分、递归反 演、广义线性、宽带约束反演、稀疏脉冲等；非线性反演，包括随 机模拟、随机反演等。 2.2 稀疏反演方法 常规反褶积得到的地震资料分辨率低，由于地震子波仅有有限 宽度的频谱，而反射系数序列具有无限宽度的频谱，使得利用地震 数据求取反射系数序列具有高度的不适定性，特别是具有无穷多个 解能够拟合地震数据。此外地震数据往往含有随机噪声，会使得地 震反演具有不稳定性。为了提高地震数据的分辨率，更客观地反映 地下介质的分布规律，需从地震数据中可靠地消除地震子波的影响， 这就要求对地震反演问题进行正则化处理或加入先验知识加以约束。 考虑到地球介质沉积过程的成层特征，可合理假设反射系数序列具 有一定的稀疏性，而稀疏性正是各沉积地层厚度不同特征的数学表 征，即利用等时间间隔表示反射系数序列时，反射系数序列中某些 时刻的反射系数肯定为零。利用反射系数序列稀疏性，实现地震数 据的反演，可在地震数据含有噪声条件下，有效地消除地震子波的 影响，提高地震数据分辨能力，提高对地下介质真实分布情况的刻 画能力。可以说，地震数据的稀疏反演构成当前提高地震分辨能力 的主要途径。 2.2.1 匹配追踪 常规反褶积得到的地震资料分辨率低，子波及反射系数准确度 不高。因此为了提高分辨率，要求反射系数能够稀疏化。但同时， 稀疏化后的反射系数再经过褶积过程会与原始数据产生偏差，导致 基于迭代再加权最小二乘的地震资料稀疏反演方法 16 信噪比降低。因此，需要寻求一种方法，在提高分辨率的同时信噪 比不会显著降低，使分离的子波及反射系数更接近真实情况。 对于稀疏解的求取，一类常用的方法是匹配追踪算法 （matching pursuit, mp） （davis et al，1997；mallat et al，1993） 。对于常见的信号处理过程，通常可以表达为 bar 的形式，其中是 a 一个 变换矩阵， a 的每个列向量称为一个原nm 子。 b 是时空域中的信号， r 是在变换域中的系数，是求解的对象。 通常 a 是欠定的，所以存在无穷多满足上式。匹配追踪法作为一种 贪婪算法，能够迭代的求出最稀疏的。因为 b 表示为 a 的列向量的 线性组合，所以最稀疏的表示方式是唯一的。令 是 a 的ni,.21 , 列向量，将上式展开得 。nxxab21 匹配追踪法的每次迭代包括两步： 1. 寻找当前剩余的相关最大的列，然后求出当前列的系数。 2. 升级剩余，即由当前的剩余，减去选出的列，得到新的剩余。 首先令初始剩余 ，令初始 r 为 0 向量。在第 k 步迭代，选br0 出与当前剩余 相关绝对值最大的列，即满足：1k iknik ara,maxrg11 之后求出第 k 步的剩余 ，求出新原子的系数kkkr,11 ，然后继续下一步的迭代。由此可见，匹配追踪法kar,1 算法简单，易行。缺点是循环量大，即使某个原子被选出来，但是 其系数没有经过优化，所以，在以后的迭代中还会选出这个原子。 这样就使得迭代步数增加。 第二章 波阻抗与稀疏反演概述 17 匹配追踪类方法还包括正交匹配追踪方法及其变形（davis et al，1997；donoho et al，2006；tropp，2004） 。正交匹配追踪法 比匹配追踪法多一个正交化的过程，一旦选出上述过程中要求的原 子，就执行正交化过程求出最优系数，以后的迭代就不会再选这个 原子，因此比匹配追踪法效率要高。 匹配追踪方法保证每次迭代使得解更加逼近真实解。如果解非 常稀疏，则匹配追踪类算法的效率很高。因此对于解严格稀疏或者 维数很小的问题匹配追踪类方法效率较高。在存在适当噪音的情况 下匹配追踪类方法也仍然适用（donoho et al，2006） 。 2.2.2 稀疏反演和压缩感知 求解线性方程组的稀疏解最直观的数学表示为如下优化问题 ， bar （2-2-1） 其中 代表非零元素的个数。问题（2-10）是一个组合优化问题，0 不存在多项式解法，需要计算所有的解才能找到最稀疏的解，但是 这样做是不现实的，在实现上意义不大（candes et al，2005） 。 一般采用以下问题求稀疏解（chen et al，1998）： ， 21.minbarts （2-2-2） 其中 代表绝对值（l 1范数） ，该问题便是基追踪问题（basis 1 pursuit，bp） ，是问题（2-10）的松弛，并且在一定的情况下，问 题（2-2-2）的解就是问题（2-2-1）的解。这样就把一个组合优化 问题变为一个连续的凸优化问题求解。基于 l1范数的算法是稀疏优 基于迭代再加权最小二乘的地震资料稀疏反演方法 18 化方法的主流，也正是本文求解的问题。 问题（2-2-2）的更一般的情形是假设采样的数据存在噪音 ，n 即 ，nard 因此相应的问题变为 ， 21. mibrts （2-2-3） 其中 是非负的实参数，用来控制噪音 的大小。问题（2-2-3）被 n 称为基追踪去噪问题（basis pursuit denoise, bpdn） 。由于问题 （2-2-2）可以变为线性规划，问题（2-2-3）可以变二阶锥优化问 题。因此问题（2-2-1）和（2-2-3）都可采用内点类算法求解 （chen et al，1998） 。 问题（2-2-3）可以变为无约束优化问题 21mindar （2-2-4） 问题（2-2-4）是一个惩罚的最小二乘问题，其目标函数是拉格朗日 函数，其中 称为拉格朗日乘子。根据对偶理论，当选取参数 和 合适时，问题（2-2-3）和（2-2-4）是等价的。问题（2-2-4）还 存在另一个关系紧密的问题 1 2 .minrtsba （2-2-5） 问题（2-2-5）在统计学中被称为 lasso 问题（least absolute shrinkage and selection operator） 。在统计理论中，最初是由文 章（tibshirani，1996）提出该问题并提出了一种组合二次规划迭 第二章 波阻抗与稀疏反演概述 19 代算法。很明显，如果参数 ， 和 选的合适，问题（2-2-3） ， （2-2-4） ， （2-2-5）是等价的，因此可以通过这三个问题中的任意 一个求出稀疏解。问题（2-2-5）可以变为二次规划问题，然后采用 内点法或投影法求解（chen et al，1998；wang et al，2009） 。 由于问题（2-2-2）目标函数不是光滑的，且只有噪声水平已知 时才有吸引力，因此其应用受到限制。问题（2-2-5）因为很少情况 下给出约束的上界，因此应用也受到了限制，更多的算法是将其等 价变形为问题（2-2-2） ，采用拉格朗日法或者内点法求解。 由于问题（2-2-2）是无约束问题，因此很多算法对其求解。这 些方法包括内点法（chen et al，1998；nocedal et al，1999） ， 坐标下降法（fu，1998；shevade et al，2003） ，投影梯度类方法 （schmidt et al, 2007） ，积极集法（perkins et al，2003） ，序 列二次规划（sqp） （lee et al，2006） ，迭代软阈值法，迭代重赋 权法（gersztenkorn et al，1986）等等。 文（kim et el，2007）中的内点法将问题（2-2-2）改为二次 规划并且用预条件共轭梯度法求解。文（wang et al，2005；wang et al，2009b）利用内点法来解遥感中的稀疏优化问题。 以上问题都可以归结为 模型：qpl 0,minprbarq （2-2-6） 可见问题（2-2-2）是问题（2-2-6）的特例（wang et al，2009a） 。 稀疏优化算法种类繁多，但是主要是基于 l1范数的凸优化方法。 沿用上文记号，褶积的矩阵向量乘积形式为： 基于迭代再加权最小二乘的地震资料稀疏反演方法 20 bnar （2-2-7） 其中 是子波延拓成的矩阵， 为反射系数， 为噪音， 为地震记ar d 录。假设噪音是白噪声，则该方程组的求解可以变为最优化问题 2minbar （2-2-8） 假设地震子波 已知或反射系数序列 可以准确估计，通过褶积wr 模型，这样也可直接估计出反射系数序列 ，使得地震记录与褶积结 果达到最大程度的拟合。这里假设噪声近似为零，在时域中反褶积 的目的是使残差最小。但是由于要求反射系数是稀疏的，因此需要 在模型中增加稀疏约束或采用其它策略求出稀疏解。 若反射系数是稀疏的，通过引进 l0范数可以将解方程问题变为 一个组合优化问题： 0minr .tsba （2-2-9） 其中 用来衡量噪音水平的大小。 这个组合优化问题在一定条件下可以变为一个等价的 l1范数优 化问题求解 1 minr .tsba （2-2-10） 但是由于目标函数不是连续可微的，因此需要变形后采用内点法， 投影梯度法，同伦法等方法求解（chen et al，1998；donoho et al，2006；figueiredo et al，2007） 。求解线性问题的稀疏解的方 第二章 波阻抗与稀疏反演概述 21 法有很多种，在前面已经做了简单的介绍，在本文子波已知的前提 下，这些方法可以都可用来解稀疏约束的反褶积问题。在地震勘探 中，稀疏反褶积有以下方法：基于 l1范数约束的方法，基于 cauchy 范数的方法，基于 huber 范数的方法和基于 l2范数的加权迭代最小 二乘法。 针对问题（2-2-10） ，可以采用一个 l1范数的连续可微的近似函 数来作为目标函数。其中最常用的有 cauchy 函数和 huber 函数。基 于 cauchy 范数的模型为 niicauchyrfbar12min （2-2-11） 其中 为正则参数， 越大，解越稀疏；刘喜武等（ 2006）针对问 题（2-2-11）基于预条件共轭梯法来求反射系数序列，提高了运算 速度，降低了内存使用量。此外基于模型： 同样可22mindrbar 以通过可以采用基于共轭梯度法的迭代重赋权最小二乘法得到稀疏 解（sacchi，1995） 。 本文采用的方法，是将 (2-2-12)1minr.ts2bar 等价的转化为 (2-2-13)21 ir 反射系数的稀疏化遇到的最直接的一个问题就是当原始数据不 满足 nyquist 采样定理时是否能够得到完整的反射系数。受到实际 条件的约束，采样时间和采样个数经常受到限制，采集的数据量往 往不满足采样定理。因此需要打破采样定理的限制。 如果采样过程是线性的，则采样过程可以写为 基于迭代再加权最小二乘的地震资料稀疏反演方法 22 bar (2-2-14) 其中是 r 原始数据， b 为采样的数据， a 是采样矩阵。地震数据 采集中，当采样个数满足采样定理时， ，则 。若采样的数ibr 据是不完整的，则该问题成为一个欠定问题，因此是病态的反问题。 而压缩感知理论可以由采样不足的数据重建出原始数据。该理 论基于两个假设条件： 第一个条件是原始数据在某个变换下是稀疏的。即假设存在变 换 使得 ，其中 s 的非零元素的个数很少。于是可以先求出稀sr 疏的 s 再恢复原始信号。即我们先解以下问题 (2-2-15)bsb 其中 ，然后求出 。abr 第二个条件是必须满足约束等距性（restricted isometry property） （candes et al，2005） 。如果满足 b 约束等距性，且采 样个数满足一定条件时，则稀疏解可以求出。如果以上两个条件满 足，我们就可以由上式求出稀疏解 s（candes，2006） 。 压缩感知理论给出了采样不足时稀疏反褶积的理论保证。 总之，对于压缩感知与稀疏反演而言，对于含噪声的数据处理 时，均可将待优化的目标函数定义为下式： 0,min21qpxbaxqp (2-2-16) 对于 p、 q 取值不同的上述优化问题，人们已开展大量的研究工 作，特别是压缩感知理论与技术的出现，直接推动着计算科学、信 息技术、通讯技术和编码与解码技术等多个领域的发展，可以说压 第二章 波阻抗与稀疏反演概述 23 缩感知已成为当前最热点的研究领域。就比较成熟的优化算法而言， 主要存在两种优化技术或策略，即将上述优化问题转化为等价的有 约束优化和无约束优化问题，有约束优化的主要算法是原始-对偶内 点法，无约束优化的主要算法是迭代再加权方法。这两种方法各有 不同优势与不足，但无论在研究领域和实际应用领域中，始终处于 相互竞争与相互促进的过程，本文重点对迭代再加权方法在地震数 据反演开展研究。 第三章 迭代再加权最小二乘地震资料稀疏反演方法 25 第 3 章 迭代再加权最小二乘地震资料稀疏反演方法 迭代再加权方法是解无约束优化问题的重要方法，其核心思想 是通过计算一系列权函数不断更新的加权最小二乘问题得到目标函 数的最优解。虽然在迭代过程中要求解一系列的线性方程组，但相 比内点法等其他方法其计算量要小很多。因此迭代再加权方法的主 要优点在于计算效率高，占用内存小，编程实现较为容易。 3.1 迭代再加权方法原理 一般的讲，对于稀疏反和压缩感知演而言，目标函数可以定义 为下式： (3-1-1)21 minbaxxqq 常见的情况是 ：20,1 (3-1-2)2 inbxxq ；考虑在该 范数下的最优化问题： njjx122 l (3-1-3)22 minbaxxq 其中 是一个足够大的参数，用于平衡地震数据所含噪声与0 稀疏性，对于该问题进行下述正则化处理： (3-1-4)nqrxbax, in2,2 其中 ； 是逐渐趋于 0 的，以逼近 。给定 njjqx12/, )( 0qx 后，记 是上述问题的一个特定解，对于该无约束问题，它满, 基于迭代再加权最小二乘的地震资料稀疏反演方法 26 足梯度方程： 0)(2)(2 ,1,2/1, baxxxqx qtnjqjjq 这是保证最小化的一个必要条件。由于其不是线性的，因此需使用 迭代的方法对上述方程进行线性化处理，方能求解上述方程。采用 再加权技术实现线性化，即设起始点为 ，通过以下迭代过程求解)1(x ：)1(kx (3-1-0)(2)(2 11)(2/1)1( baxxxqktnjkjkjk 5） 或者通过 （3-ba xnjxqdiagt kjqkjt )1(2/1)( ,.,21 1-6） 来求解。在迭代过程中，若记.3,21k 2/1)(1qkjkjxw 则 为 的权重，则（3-1-5）和（3-1-6）两式可以改写为：iw)1(kx 0)(22 11)(1)( baxxktnjkjj （3-1-7） 第三章 迭代再加权最小二乘地震资料稀疏反演方法 27 ba xnjwdiagt kjkj )1(1,.,21 （3-1-8） 可以看到 是随着迭代的进行而不断变化的，同时，为了满足对1kjw 的逼近