《基于迭代再加权最小二乘的地震资料稀疏反演方法—硕士学位论文.doc》由会员分享，可在线阅读，更多相关《基于迭代再加权最小二乘的地震资料稀疏反演方法—硕士学位论文.doc（70页珍藏版）》请在三一办公上搜索。
1、中文摘要地震数据反演的目的在于通过求解反射系数序列来反映地下介质的分布规律，这就需从地震数据中可靠地消除地震子波的影响。然而受到有限带宽的影响，往往缺失反射系数序列中的低频和高频信息；同时由于受到不可避免的噪声的影响，反演得到的结果分辨率通常比较低。本论文的目的在于使用迭代再加权的寻优方法由地震数据得到高分辨率的反射系数序列；该方法通过再加权的方法将稀疏优化问题转化为一系列权不断变化的线性方程组，在保证结果精度的前提下增加了求解的速度。对寻优过程中的线性方程组求解采用针对性的预条件化，并在权重的计算中引入伸缩变换，优化了整个寻优过程的稳定性和求解速度。经过大量的数值试验，得出了合理的参数选择范
2、围，通过与真值的比较以及实际资料处理验证了本方法的实用性。关键词：波阻抗反演，迭代再加权，Gauss-Siedel迭代，共轭梯度法，预条件化，伸缩变换ABSTRACTThe 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
3、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
4、 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 th
5、e 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, an
6、d the practicability is conformed together with the real data processing.KEYWORDS:impedance inversion, iterative reweighted, Gauss-Siedel iteration,conjugate gradient method, pre-conditioning, shrink目 录第一章 引 言11.1选题依据和研究意义11.2 波阻抗反演研究现状11.3 稀疏反演方法的研究情况31.4 论文结构与主要内容41.5 论文的主要创新点4第二章 波阻抗与稀疏反演概述52.1 基
7、本概念52.1.1 常规波阻抗反演技术的基本假设前提52.1.2 波阻抗反演方法的基本原理52.1.3波阻抗反演的分类102.2 稀疏反演方法112.2.1 匹配追踪112.2.2 稀疏反演和压缩感知12第三章 迭代再加权最小二乘地震资料稀疏反演方法173.1 迭代再加权方法原理173.2 反射系数稀疏反演的迭代再加权方法203.3 解线性方程组的迭代方法213.3.1 Gauss-Seidel迭代方法213.3.2 共轭梯度法223.4 收缩变换263.5 预条件化283.5 流程图30第四章 数值模拟和实际资料314.1 不同噪声水平下的原始地震数据314.2 不同主频子波反演结果324.
8、3 噪音水平不同时反演结果344.4 权重因子，噪声水平不同时反演结果374.5 范数逼近程度不同时反演结果394.6 实际资料处理40第五章 结论与认识41参考文献43作者简介47攻读硕士学位期间已发表及待发表的论文48攻读硕士学位期间参加的科研项目48第一章 引 言1.1选题依据和研究意义油气资源关系民生和国家安全，与国民经济的大部分产业有直接或间接的联系，对经济发展有着十分轻重的影响，也是今后相当长时间内难以替代的主要能源。所以加强国内油气的勘探能力具有十分重要的意义。地震、测井、钻井等方法是石油工作者了解地下地质构造、地层、岩性、物性、含油气性等重要的信息来源。虽然地震方法获得资料的分
9、辨率远远不及测井、钻井，但是随着地震勘探技术进一步地发展，地震资料的信噪比、分辨率、成像准确性都获得了很大的提升。同时地震资料包含大量地下地质信息，且覆盖面积广，具有三维特性，所以地震勘探愈来愈受到重视。地震勘探的主要原理是从地表某一点位激发地震子波，子波在向下传播的过程中遇到波阻抗界面后会发生反射，传回一部分能量，被传感器接收。随着这一激发、传播、反射、接收的过程不断进行，传感器记录下一系列按时间排列的地震波，即地震记录。地震勘探的基础是地下不同地层存在波阻抗差异，因为地震波只有传播到有波阻抗差异的地层分界面时才会发生反射。地震资料反演其基本目的是充分利用测井、钻井、地质资料提供丰富的构造、
10、层位、岩性等信息，由常规的地震剖面推导出地下地层的波阻抗、密度、速度等信息，为勘探开发提供重要的依据。而石油勘探开发的形式在不断地发展，计算机技术也不断地提高，使得地震资料反演技术成为了当前热门的研究课题。其中波阻抗反演（Acoustic Impedance，AI）是反演技术中很重要的一个部分，它利用地震资料反演反射系数、进而得到地层波阻抗信息的技术。同时，地震资料反演技术本身具有的不适定性，应用条件的限制以及储层的复杂性，常规叠后地震资料品质不高且分辨率低等因素，一定程度上影响了该技术更广泛的应用。1.2 波阻抗反演研究现状地震反演的基本目的是利用地震波在地下介质中的传播规律，通过数据采集、
11、处理与解释等流程，推测地下岩层结构和物性参数的空间分布，为勘探开发提供重要的依据13。波阻抗反演是地震反演的一个重要组成部分，由于波阻抗信息是联系地质和地球物理的一座桥梁，叠后计算数据量相对较少，在实际生产中应用方便而且效果明显。因此，波阻抗反演在地震反演中具有特殊的地位。侠义的地震反演概念指的就是波阻抗反演4。波阻抗反演是利用地震资料反演地层波阻抗（或速度）的地震特殊处理解释技术，它产生于20世纪70年代，80年代开始得到蓬勃发展。1983年，Cooke介绍了地震资料广义线性反演方法，从而揭开了波阻抗反演技术的新篇章5。周竹生等人在90年代初期提出了综合利用地质、地震和测井资料进行约束反演，
12、克服单一的线性反演方法的缺陷6。90年代中期，李宏兵在国内提出了一条走递推反演与宽带约束反演结合的道路，此方法为我们解决了从单道出发的反演方法不能在跟不上消除噪音的困惑7。在此基础之上，有人进行了无井多道反演和有井多道反演的研究8,9，使波阻抗反演方法更加完善。90年代至今，围绕一维波阻抗反演的各类算法及应用成果层出不穷，随着研究的开展，在1997年左右开始出现了一些反思的文章，指出了波阻抗反演中存在的一些缺陷，并提出了一些解决方案1012。BP Amoco公司的Connolly在1999年正式发表了弹性波阻抗反演方法的论文，该方法最早是在19931994年由BP Amoco公司发展起来，随后
13、在2000年的SEG年会上同时出现了四篇论文对EI(elastic impedance)进行了研究，ARCO公司介绍了他们申请专利的弹性波阻抗反演方法，认为在求取的反射系数的稳定性方面要好于Connolly方法，而且计算的EI和AI数值在一个尺度下，同时BP Amoco公司在会上提出了扩充弹性波阻抗方法，可以用于流体和岩性的预测。此外，Paradigm公司在其商业软件Vanguard的最新版块中也有有关弹性波阻抗反演功能的出现。Jason公司则推出了Rock Trace弹性反演模块，以纵波波阻抗和横波波阻抗的概念来区别其以往软件中的波阻抗概念。这些进展说明弹性波阻抗已经成为波阻抗反演进一步发展
14、的方向之一，地震反演的发展正走向AI和EI结合、AI和AVO结合的道路13。这个发展过程大致上经历了从直接反演到模型反演、从叠后反演到叠前反演、从线性反演到非线性反演发展的过程16,33。近十几年来，非线性反演方法得到了有效利用。实践表明，非线性反演比线性反演更接近真实。非线性优化方法有基于导数的最陡下降法、牛顿法、共轭梯度法等；基于非导数的神经网络方法、遗传算法（GA）、模拟退火方法（SA）、随机搜索等17。为了提高反演的精度和收敛速度，针对非线性反演方法，地球物理工作者将各种最优化方法引入到了波阻抗反演中，各种各样的混合优化算法也陆续出现，力求在解决计算精度的同时提高算法收敛速度，减少反演
15、多解性18。尽管已经作了大量的研究工作，波阻抗反演方法仍然存在许多的不足，如算法复杂、实现困难、抗干扰能力差(从而导致分辨率降低)、收敛速度慢等，因此仍需要进一步研究和改进这些优化反演算法。在理论指导下，国内外相继研制出一些商业化程度高的波阻抗反演方法及软件，如STRATR、GLOG、JASON等19。目前国内使用的波阻抗反演软件基本上都是引进软件，国产软件较少。96年以前的软件一般都是单井控制下的地震资料反演，只能解决构造简单没有断层或岩性、岩相横向变化小的问题，在多井条件下和地质条件复杂的地区难以推广使用。国外从80年代末提出了利用声波测井资料作为约束条件，对过井地震剖面进行正反演联合迭代
16、求取地下波阻抗的方法，将反演方法推向非线性，这种新方法利用了测井资料的高频信息，因而表现了强劲的发展势头，目前成为各国外石油公司重点发展的对象，如美国HGS公司宽带约束反演BCI技术，法国CGG公司波阻抗反演模拟ROVIM技术，俄罗斯的PARM,加拿大的STRATA，丹麦的ISIS，以及荷兰推出的JASON技术。今后地震反演的技术，是在现在的基础上，充分利用测井、地震、地质的信息，继续探索先进的算法，使储层横向预测更符合实际。这些方法在实际中取得了一定的应用效果，但是还不能够完全满足实际需要。不同反演方法结果差别较大，不同人使用同一种方法结果也有差别，往往导致错误的地质解释。究其原因，除了地震
17、资料品质、子波提取不准确、垂直入射假设与实际有误差等因素外，还有另外两个方面的主要影响因素：其一，由于地震数据具有带限性质，直接反演只能获得波阻抗模型中的中频成分，缺少应有的低频和高频成分。其二，Hardmard意义下，波阻抗反演是一个不适定问题，其反演结果会出现不稳定性和多解性。目前，解决带限问题最有效的手段是对解加以约束进而进行反演（如测井数据的约束）。而解决不适定问题主要使用正则化方法并辅之以适当的最优化技巧15。通常利用测井资料在地质知识约束的插值方法，来弥补地震资料中仍缺乏的低频和高频信息。1.3 稀疏反演方法的研究情况线性问题的稀疏解在很多领域都有研究。线性问题可以分为正定问题，欠
18、定问题和超定问题。正定问题的解存在并且唯一。欠定问题存在无穷多解，因此一般研究其最简单的解，即稀疏解。超定问题不存在严格的解，一般求出满足最小二乘意义的近似解，当然也可先验约束研究其稀疏近似解。目前的稀疏算法主要针对欠定的线性问题求解。在实际应用中求稀疏解的最常用方法是匹配追踪类方法（Davis et al，1997；Donoho et al，2006；Tropp，2004），主要有匹配追踪，正交匹配追踪和其他改进算法。这类方法的执行简单，但是计算量巨大导致计算速度慢，尤其在解不是严格稀疏和数据存在误差的情况下该问题尤其突出。稀疏解问题可以描述为一个目标函数为l1范数最小化的约束优化问题，被称
19、为基追踪问题（Chen et al，1998）。这样就可以用常见的凸优化算法求解。基追踪问题可以变为线性规划和二次规划问题，然后采用标准的内点法求解（Chen et al，1998；Kim et al，2007；Wang et al，2009b）。1.4 论文结构与主要内容本文针对基于迭代加权的反演问题，系统的介绍了迭代加权的波阻抗反演方法，具体内容有：（1） 对波阻抗反演方法涉及到的基本概念、原理进行阐述，并介绍了几种常用的反演方法；（2） 从介绍波阻抗反演的反褶积分析入手，分析并导出了基于迭代加权的波阻抗反演的基本问题，并介绍了迭代方法以及迭代的效果。1.5 论文的主要创新点本文对迭代波阻
20、抗反演的反褶积问题进行了进一步的分析，从而导出了基于迭代再加权的波阻抗反演方法，并且应用了几种迭代方法来解方程组，通过理论模拟和实际资料处理证明了新方法在收敛速度上的优越性。(1) 对波阻抗反演方法进行了系统分类，对稀疏反演方法以及与之相关的匹配追踪和压缩反演进行了论述；(2) 从反褶积问题入手，建立了一种基于迭代加权的波阻抗反演方法，分析了加权因子的选取，并研究了Gauss-Seidel迭代法、共轭梯度法和预条件加速技巧；(3) 对以上方法，用不同噪音水平的地震数据进行了反演计算，计算结果验证了算法的有效性。第二章 波阻抗与稀疏反演概述2.1 基本概念2.1.1 常规波阻抗反演技术的基本假设
21、前提目前常用的波阻抗反演软件所采用的方法基本上都是在褶积模型的基础上建立的，因此，要求资料都要满足褶积模型的假设前提，可以概括为以下四个方面。（1） 地震模型：假设地层是水平层状介质、地震波为平面波法向入射，其地震剖面为正入射剖面，并假设地震道为地震子波与地层反射系数的褶积。（2） 反射系数序列：普通递推反演中，假设反射系数为完全随机的序列；在约束稀疏脉冲反演中，假设反射系数为由一系列大的反射系数叠加在高斯分布的小反射系数的背景上构成。（3） 噪音分量：通常假设波阻抗反演输入的地震数据的振幅信息反映了地下波阻抗变化的情况，地震剖面没有多次波和绕射波等噪音分量。因此在资料处理时，可以考虑的处理流
22、程是反褶积、噪音剔正，尤其是多次波，处理的最终目标是得到高分辨率的保幅波场。（4） 地震子波：假设反射系数剖面中的每一道都可以看作地下反射系数与一个零相位子波的褶积。实践情况下地震子波往往是混合相位或最小相位的，因此需要对地震剖面进行相位校正处理。2.1.2 波阻抗反演方法的基本原理波阻抗波在某种介质中传播的速度与介质密度的乘积，在本论文中其符号和定义为，Z表示波阻抗，表示密度，V表示波速。根据弹性波理论，当一个纵波入射到波阻抗不同的两种介质的分界面时，一般要产生反射纵波、反射横波、透过纵波、透过横波四种新的波动。一般主要研究纵波反射波；在水平层状介质、地震波为平面波法向入射的假设下，问题会得
23、到简化，只考虑产生的反射纵波和透射纵波、反射横波和透射横波的情况。地震子波震源附近地震波以冲击波的形式传播，当传播到一定距离时，波形逐渐稳定，此时的地震波被称为地震子波。进行反演时，求取地震子波的方法一般有统计法和参数法。统计法是从实际地震记录中统计出地震子波的振幅谱和相位谱，从而生成地震子波；参数法是给出某种地震子波类型的参数而生成地震子波（如Ricker）子波。Ricker子波是最常用的零相位子波，在时间域表示为 （2-1-1） 其中，是主频。零相位子波分辨力较高，最有利于解释。 由定义可见，只有层状介质的波阻抗存在差异时反射系数才不为0。当子波遇到波阻抗界面后，会发生反射，从而被接收到。
24、由于是按照时间序列排列的，以下分别记为接收到的地震数据、子波、反射系数按照时间的序列，并已经过离散化，则按照定义，可以得到：L表示子波的长度。数学上定义此运算为褶积；下记“*”为褶积运算符，则上式可以表示为 （2-1-2） 因此，地震记录实际上是由反射系数与子波的褶积得到的。如果将子波作为输入信号，地下介质作为一个线性时不变的“系统”（或者滤波器），接收信号为信号通过系统后的输出，则该褶积可以看做是一个系统的响应过程（陆基孟，1993）。该响应不但包括一次反射还包括多次放射。假定震源波形在向地下传播的过程中不变，因此褶积模型不包括固有的吸收。地震记录可以看作是由震源子波与地下反射率函数、多次反
25、射、仪器等诸多因素的褶积的过程。虽然可以给合成地震记录加上随机噪音，但是这里的用来建立反褶积滤波器的褶积模型暂时不包括随机噪音。图2-1 模拟的反射系数、子波和合成的地震记录 如果能求得一个算子是子波的逆，即：。则它与地震记录褶积时，可由（2.1）式得 （2-1-3）即地震记录和子波的逆进行褶积可以得到反射系数，这个过程叫做反褶积。由以上过程可以看到，反射系数序列中含有波阻抗随时间变化的信息，这就提供了速度和密度随时间变化的信息，因此可得到地层、岩性及构造在地下空间分布的信息，在有利条件下还可得到岩石孔隙率、渗透率、孔隙流体性质（油、气、水）乃至地层压力的信息。波阻抗反演需要求出反射系数，因而
26、其反演就是已知地震记录，由提取的子波反褶积得到反射系数序列。在使用炸药作为震源时，震源实际是一个脉冲，由于地层的滤波作用，震源信号在地下传播过程中高频被衰减，能量被吸收，变成一个子波，所以在地震勘探中观测的信号其峰值频率往往小于100Hz，这就导致了记录信号的分辨率降低，记录到的反射是减弱了强度的宽脉冲，并且通常随着反射时间的增加信号越宽和越弱。子波在地层中传播，携带着反射系数序列这种有用的地质信息返回地面，只有通过反褶积消除子波才能恢复反射系数序列；同时反褶积可以压缩地震子波，压制交混回响和短周期多次波，从而提高时间分辨率。反褶积可应用于叠前资料，也可应用于叠后资料。理想的反褶积应该压缩子波
27、并消除多次波，在地震道内只留下地层反射系数。子波压缩可以通过将反滤波器作为反褶积算子来实现，反滤波器可以将地震子波转变成尖脉冲，每个脉冲的强弱与界面的反射系数的大小成正比，脉冲的极性反映界面反射系数的符号。当应用于地震合成记录时，反滤波输出应为地层脉冲响应。 反射系数根据弹性波理论，当一个纵波入射到波阻抗不同的两种介质分界面时，一般要产生反射纵波、反射横波、透过纵波、透过横波四种新的波动。在实际生产工作中，主要研究纵波反射波。界面倾角不大时，可近似认为是垂直入射。在垂直入射的情况下，将会使问题简化，纵波入射时我们将只考虑产生的反射纵波和透射纵波、反射横波和透射横波的情况。这时界面的反射系数定义
28、为 （2-1-4）其中AF，AR分别表示反射振幅和入射振幅。而；分别对应上层和下层的波阻抗、密度、速度。当平面波非垂直入射时，反射系数将随着入射角的变化而变化，这种变化与界面两边介质的波阻抗和各种弹性参数有关，特别是与介质的泊松比有很密切的关系（Zoppritz方程，入射角等于零是它的特例）。褶积模型根据介质模型假设不同，波阻抗反演分为基于波动方程的反演和基于褶积模型的反演两大类。前者由于算法结构复杂、计算量大、抗干扰差，很难获得一个稳定解，在实际中未得到广泛应用。目前主要使用基于褶积模型的直接反演方法，它算法简单，对地震噪音敏感性小，但一般情况下并不一定能得到一个稳定的解。（1） 褶积：在时
29、域中为了便于求得线性时不变系统的零状态响应，如果能够将任意激励信号分解为单元信号，并使每一单元信号在系统中产生的零状态响应易于求得，那么根据系统的线性时不变特性，就可以利用叠加原理很容易求得原信号在系统中产生的零状态响应。系统对单位冲击函数的零状态响应为，对于线性时不变系统，当激励是强度为的延时冲击函数时，其零状态响应为，任意信号可表示为 （2-1-5）作用下系统响应为 (2-1-6)其中为积分变量，为参变量，上式褶积积分又可记作： （2-1-7）（2）模型产生：水平层状介质模型是共中心点（CMP）叠加方法的基础，在CMP叠加中，对分享同一炮检中点的一些地震道进行正常时差校正求和，产生一个逼近
30、1D层状介质垂直入射平面波响应的求和道。本着这种数据，形成了一种地震道模型。假定从波阻抗不连续界面反射（或透射）的子波与入射子波具有同样的波形，一个地震道顺序记录这样一系列连续子波，也就记录了波阻抗不连续面。地震道可以简单地看成是每一个独立的反射界面叠加的结果，这里每个反射界面都认为是一个滤波器。地震记录可认为是震源子波与介质反射系数（即地下垂直入射反射系数序列）的褶积。1）时间域表达：连续情形： (2-1-8)其中为地震子波，为地震反射序列，为地震记录，为随机噪音。 矩阵-向量情形： （2-1-9）其中为地震反射系数向量，为噪音列向量，为地震波记录向量，为子波阵，子波长度为, 2） 频率域表
31、达： 连续情形： （2-1-10）各表示傅立叶变化下的：-地震记录，-子波，-反射系数，-噪音。2.1.3波阻抗反演的分类波阻抗信息是联系地质和地球物理的重要媒介，其叠后计算数据量相对较少，在实际生产中应用方便、效果明显，因此通常的地震反演概念指的就是波阻抗反演。将地震剖面转换成波阻抗剖面，不仅便于解释人员将地震资料与测井资料连接对比，而且能有效地对储层物性参数的变化进行研究，从而得到物性参数在空间上的分布规律，对油气的勘探和开发有很重要的指导意义。对于波阻抗反演的重要性，一个很经典的论述是著名的地球物理学家李庆忠院士所说的“交到地质人员手中的地震资料应是作了反演的波阻抗剖面”、“波阻抗反演是
32、高分辨率地震资料处理的最终表达形式”25,26。而常常被人们引用的美国阿莫克石油公司在墨西哥湾海上寻找非背斜油藏不断取得成功的例子29也印证了这一说法。波阻抗反演是目前应用最广泛的储层预测和油藏描述手段之一，其发展日新月异，计算方法多种多样，这些方法由于其基于的数学物理模型不同，各有优缺点，利用的基础资料和能够达到的反演效果是不同的。波阻抗反演方法有多种（Seislog,Velog,Delog,PIVT,BCI,ROVIM,SLIM,PARM和Strata）14,概括起来有两大类：基于反射系数逆公式的直接反演和基于正演模型的迭代反演。迭代反演又可以分为无井条件下的广义性反演（GLI）和有井条件
33、下的宽带约束反演（BCI）33。还有一些属于非线性反演范畴的方法，如混沌理论、神经网络等。（1） 依据算法的实现方式分为：基于反褶积的反演，包括道积分、递归、广义线性、稀疏脉冲反演等；基于波动方程的反演，包括Born反散射反演等；基于随机过程的反演，包括随机反演、随机模拟、模拟退火反演等；基于特征分析的反演，包括特征反演、神经网络反演等；基于动力学特征的反演，包括混沌反演等（2） 按测井和地震的相对作用：无井反演，包括道积分等；地震为主的测井约束类反演，包括递归反演、广义线性、宽带约束反演、稀疏脉冲、模拟退火反演等；测井为主的统计类反演，包括随机模拟、随机反演等；井-震约束联合反演，包括特征反
34、演、神经网络反演等。这种分类方案体现了测井资料对反演约束作用的强弱，便于研究人员根据测井和地震基础资料评价各种反演方法的反演效果。（3） 依数学算法不同可分为：线性反演，包括道积分、递归反演、广义线性、宽带约束反演、稀疏脉冲等；非线性反演，包括随机模拟、随机反演等。2.2 稀疏反演方法 常规反褶积得到的地震资料分辨率低，由于地震子波仅有有限宽度的频谱，而反射系数序列具有无限宽度的频谱，使得利用地震数据求取反射系数序列具有高度的不适定性，特别是具有无穷多个解能够拟合地震数据。此外地震数据往往含有随机噪声，会使得地震反演具有不稳定性。为了提高地震数据的分辨率，更客观地反映地下介质的分布规律，需从地
35、震数据中可靠地消除地震子波的影响，这就要求对地震反演问题进行正则化处理或加入先验知识加以约束。考虑到地球介质沉积过程的成层特征，可合理假设反射系数序列具有一定的稀疏性，而稀疏性正是各沉积地层厚度不同特征的数学表征，即利用等时间间隔表示反射系数序列时，反射系数序列中某些时刻的反射系数肯定为零。利用反射系数序列稀疏性，实现地震数据的反演，可在地震数据含有噪声条件下，有效地消除地震子波的影响，提高地震数据分辨能力，提高对地下介质真实分布情况的刻画能力。可以说，地震数据的稀疏反演构成当前提高地震分辨能力的主要途径。2.2.1 匹配追踪常规反褶积得到的地震资料分辨率低，子波及反射系数准确度不高。因此为了
36、提高分辨率，要求反射系数能够稀疏化。但同时，稀疏化后的反射系数再经过褶积过程会与原始数据产生偏差，导致信噪比降低。因此，需要寻求一种方法，在提高分辨率的同时信噪比不会显著降低，使分离的子波及反射系数更接近真实情况。对于稀疏解的求取，一类常用的方法是匹配追踪算法（Matching pursuit, MP）（Davis et al，1997；Mallat et al，1993）。对于常见的信号处理过程，通常可以表达为的形式，其中是A一个变换矩阵，A的每个列向量称为一个原子。b是时空域中的信号，r是在变换域中的系数，是求解的对象。通常A是欠定的，所以存在无穷多满足上式。匹配追踪法作为一种贪婪算法，能
37、够迭代的求出最稀疏的。因为b表示为A的列向量的线性组合，所以最稀疏的表示方式是唯一的。令是A的列向量，将上式展开得。匹配追踪法的每次迭代包括两步：1. 寻找当前剩余的相关最大的列，然后求出当前列的系数。2. 升级剩余，即由当前的剩余，减去选出的列，得到新的剩余。首先令初始剩余，令初始r为0向量。在第k步迭代，选出与当前剩余相关绝对值最大的列，即满足： 之后求出第k步的剩余，求出新原子的系数，然后继续下一步的迭代。由此可见，匹配追踪法算法简单，易行。缺点是循环量大，即使某个原子被选出来，但是其系数没有经过优化，所以，在以后的迭代中还会选出这个原子。这样就使得迭代步数增加。匹配追踪类方法还包括正交
38、匹配追踪方法及其变形（Davis et al，1997；Donoho et al，2006；Tropp，2004）。正交匹配追踪法比匹配追踪法多一个正交化的过程，一旦选出上述过程中要求的原子，就执行正交化过程求出最优系数，以后的迭代就不会再选这个原子，因此比匹配追踪法效率要高。匹配追踪方法保证每次迭代使得解更加逼近真实解。如果解非常稀疏，则匹配追踪类算法的效率很高。因此对于解严格稀疏或者维数很小的问题匹配追踪类方法效率较高。在存在适当噪音的情况下匹配追踪类方法也仍然适用（Donoho et al，2006）。2.2.2 稀疏反演和压缩感知求解线性方程组的稀疏解最直观的数学表示为如下优化问题 ，
39、 （2-2-1）其中代表非零元素的个数。问题（2-10）是一个组合优化问题，不存在多项式解法，需要计算所有的解才能找到最稀疏的解，但是这样做是不现实的，在实现上意义不大（Candes et al，2005）。 一般采用以下问题求稀疏解（Chen et al，1998）： ， （2-2-2）其中代表绝对值（l1范数），该问题便是基追踪问题（Basis pursuit，BP），是问题（2-10）的松弛，并且在一定的情况下，问题（2-2-2）的解就是问题（2-2-1）的解。这样就把一个组合优化问题变为一个连续的凸优化问题求解。基于l1范数的算法是稀疏优化方法的主流，也正是本文求解的问题。 问题（2-
40、2-2）的更一般的情形是假设采样的数据存在噪音，即 ，因此相应的问题变为 ， （2-2-3）其中是非负的实参数，用来控制噪音的大小。问题（2-2-3）被称为基追踪去噪问题（Basis Pursuit Denoise, BPDN）。由于问题（2-2-2）可以变为线性规划，问题（2-2-3）可以变二阶锥优化问题。因此问题（2-2-1）和（2-2-3）都可采用内点类算法求解（Chen et al，1998）。 问题（2-2-3）可以变为无约束优化问题 （2-2-4）问题（2-2-4）是一个惩罚的最小二乘问题，其目标函数是拉格朗日函数，其中称为拉格朗日乘子。根据对偶理论，当选取参数和合适时，问题（2-
41、2-3）和（2-2-4）是等价的。问题（2-2-4）还存在另一个关系紧密的问题 （2-2-5）问题（2-2-5）在统计学中被称为LASSO问题（Least Absolute Shrinkage and Selection Operator）。在统计理论中，最初是由文章（Tibshirani，1996）提出该问题并提出了一种组合二次规划迭代算法。很明显，如果参数，和选的合适，问题（2-2-3），（2-2-4），（2-2-5）是等价的，因此可以通过这三个问题中的任意一个求出稀疏解。问题（2-2-5）可以变为二次规划问题，然后采用内点法或投影法求解（Chen et al，1998；Wang et a
42、l，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），序列二次规划（SQ
43、P）（Lee et al，2006），迭代软阈值法，迭代重赋权法（Gersztenkorn et al，1986）等等。文（Kim et el，2007）中的内点法将问题（2-2-2）改为二次规划并且用预条件共轭梯度法求解。文（Wang et al，2005；Wang et al，2009b）利用内点法来解遥感中的稀疏优化问题。以上问题都可以归结为模型： （2-2-6）可见问题（2-2-2）是问题（2-2-6）的特例（Wang et al，2009a）。 稀疏优化算法种类繁多，但是主要是基于l1范数的凸优化方法。沿用上文记号，褶积的矩阵向量乘积形式为： （2-2-7）其中是子波延拓成的矩阵，为
44、反射系数，为噪音，为地震记录。假设噪音是白噪声，则该方程组的求解可以变为最优化问题 （2-2-8）假设地震子波已知或反射系数序列可以准确估计，通过褶积模型，这样也可直接估计出反射系数序列，使得地震记录与褶积结果达到最大程度的拟合。这里假设噪声近似为零，在时域中反褶积的目的是使残差最小。但是由于要求反射系数是稀疏的，因此需要在模型中增加稀疏约束或采用其它策略求出稀疏解。 若反射系数是稀疏的，通过引进l0范数可以将解方程问题变为一个组合优化问题： （2-2-9）其中用来衡量噪音水平的大小。 这个组合优化问题在一定条件下可以变为一个等价的l1范数优化问题求解 （2-2-10）但是由于目标函数不是连续
45、可微的，因此需要变形后采用内点法，投影梯度法，同伦法等方法求解（Chen et al，1998；Donoho et al，2006；Figueiredo et al，2007）。求解线性问题的稀疏解的方法有很多种，在前面已经做了简单的介绍，在本文子波已知的前提下，这些方法可以都可用来解稀疏约束的反褶积问题。在地震勘探中，稀疏反褶积有以下方法：基于l1范数约束的方法，基于Cauchy范数的方法，基于Huber范数的方法和基于l2范数的加权迭代最小二乘法。针对问题（2-2-10），可以采用一个l1范数的连续可微的近似函数来作为目标函数。其中最常用的有Cauchy函数和Huber函数。基于Cauchy范数的模型为 （2-2-11）其中为正则参数，越大，解越稀疏；刘喜武等（2006）针对问题（2-2-11）基于预条件共轭梯法来求反射系数序列，提高了运算速度，降低了内存使用量。此外基于模型：同样可以通过可以采用基于共轭梯度法的迭代重赋权最小二乘法得到稀疏解（Sacchi，1995）。本文采用的方法，是将 (2-2-12)等价的转化为 (2
