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

分类号分类号 密级密级 UDC 编号编号 中国科学院研究生院中国科学院研究生院 硕士学位论文硕士学位论文 基于迭代再加权最小二乘的地震资料稀疏反演方法基于迭代再加权最小二乘的地震资料稀疏反演方法 申请学位级别申请学位级别 硕士硕士 学科专业名称学科专业名称 固体地球物理学固体地球物理学 Seismic Data Sparse Inversion Method Based on Iterative Reweighted Technique By Yuan Yi ming A Dissertation Submitted to Graduate University of Chinese Academy of Sciences In partial fulfillment of the requirement For the degree of Master of Science Lanzhou Center for Oil and Gas Resources Institute of Geology and Geophysics Chinese Academy of Science May 26th 2012 中国科学院地质与地球物理研究所中国科学院地质与地球物理研究所 兰州油气资源研究中心兰州油气资源研究中心 研究生学位论文独创性声明研究生学位论文独创性声明 本人声明所呈交的学位论文是本人在 导师指导下进 行的研究工作及取得的研究成果 除了文中特别加以标注和致谢的 内容外 论文中不包含其他人已经发表或撰写过的研究成果 也不 包含为获得中国科学院地质与地球物理研究所兰州油气资源中心或中国科学院地质与地球物理研究所兰州油气资源中心或 其他教育机构其他教育机构的学位或证书而使用过的材料 与我一同工作的同志 对本研究所做的任何贡献均已在论文中做了明确的说明并表示谢意 论文与资料若有不实之处 本人承担一切相关责任 学位论文作者签名 签字日期 年 月 日 学位论文版权使用授权书学位论文版权使用授权书 本学位论文作者完全了解中国科学院地质与地球物理研究所兰中国科学院地质与地球物理研究所兰 州油气资源研究中心州油气资源研究中心有关知识产权的规定 研究所有权保留并向国 家有关部门或机构送交本论文及论文的复印件和电子文档 允许论 文被查阅和借阅 中国科学院地质与地球物理研究所兰州油气资源中国科学院地质与地球物理研究所兰州油气资源 研究中心研究中心可以将学位论文的全部或部分内容编入有关数据库进行检 索和传播 可以采用影印 缩印或扫描等复制手段保存 汇编学位 论文 保密的学位论文在解密后适用本授权书 保密论文注释 本学位论文属于保密范畴 保密期限 年 解 密后适用本授权书 学位论文作者签名 导师签名 签字日期 年 月 日 签字日期 年 月 日 致 谢 衷心感谢我的导师 研究员 自从 2009 年入学以来 导师为我的培养倾 注了大量心血 从课程学习 论文的选题 研究思路到论文的撰写 倾注了老 师的心血 恩师的深入浅出的讲解 将复杂深奥的理论变得津津有味 细致入 微的关怀和帮助 使我顺利完成学业 恩师对物理 数学和地球物理等多门学 科都富有深厚的造诣而且具有广博的学识 其学术思想活跃且有卓越的前瞻洞 察力及敏锐的思维力 恩师高尚的品格 渊博的知识和创新的精神 严谨的治 学态度和废寝忘食的工作精神等等都是我的榜样 对我以后的科研 学习工作 将产生深刻的影响和启迪 谨向尊敬杨老师表示深深的敬意和由衷的感谢 并 祝身体健康 衷心感谢 等老师对我在学业上的许多指导和帮助 老师们渊博的知识 严谨的治学态度 以诚待人的高尚人品 无不给我留下深 刻的印象 感谢中科院地质与地球物理研究所兰州油气研究中心的 老师在三年来给 予我的关怀和帮助 老师乐于任事 关心学生 是一位受学生尊敬的良师益友 感谢研究生部的 老师 老师 老师 老师 老师 感谢他 们对我学习的关心和帮助 感谢中科院地质与地球物理研究所兰州油气研究中心的学长和同学们 正 是有了你们的诸多帮助 才让学习和科研的气氛变得更加地和谐和快乐 谨向 对我给予诸多帮助的油气中心的兄弟姐妹们致以诚挚的谢意 子曰 三人行必 有吾师也 从他们身上我也学到了很多东西 并且一起度过了愉快的研究生时 光 最后向所有给予我关心和帮助的朋友 以及对我研究工作给予直接或间接 帮助的论文作者 专家和学者致以诚挚的谢意 在以后的岁月里 无论身在何处 我都将保持严谨的学风 诚实做人 勤 奋工作 乐于助人 用所学知识造福社会 让世界更变得更美好一些 中文摘要 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 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 第四章第四章 数值模拟和实际资料数值模拟和实际资料 31 4 1 不同噪声水平下的原始地震数据 31 4 2 不同主频子波反演结果 32 基于迭代再加权最小二乘的地震资料稀疏反演方法 VI 4 3 噪音水平不同时反演结果 34 4 4 权重因子 噪声水平不同时反演结果 37 4 5 范数逼近程度不同时反演结果 39 4 6 实际资料处理 40 第五章第五章 结论与认结论与认识识 41 参考文献参考文献 43 作者简介作者简介 47 攻读硕士学位期间已发表及待发表的论文攻读硕士学位期间已发表及待发表的论文 48 攻读硕士学位期间参加的科研项目攻读硕士学位期间参加的科研项目 48 第一章 引 言 1 第一章 引 言 1 1 选题依据和研究意义 油气资源关系民生和国家安全 与国民经济的大部分产业有直接或间接的 联系 对经济发展有着十分轻重的影响 也是今后相当长时间内难以替代的主 要能源 所以加强国内油气的勘探能力具有十分重要的意义 地震 测井 钻井等方法是石油工作者了解地下地质构造 地层 岩性 物性 含油气性等重要的信息来源 虽然地震方法获得资料的分辨率远远不及 测井 钻井 但是随着地震勘探技术进一步地发展 地震资料的信噪比 分辨 率 成像准确性都获得了很大的提升 同时地震资料包含大量地下地质信息 且覆盖面积广 具有三维特性 所以地震勘探愈来愈受到重视 地震勘探的主要原理是从地表某一点位激发地震子波 子波在向下传播的 过程中遇到波阻抗界面后会发生反射 传回一部分能量 被传感器接收 随着 这一激发 传播 反射 接收的过程不断进行 传感器记录下一系列按时间排 列的地震波 即地震记录 地震勘探的基础是地下不同地层存在波阻抗差异 因为地震波只有传播到有波阻抗差异的地层分界面时才会发生反射 地震资料反演其基本目的是充分利用测井 钻井 地质资料提供丰富的构 造 层位 岩性等信息 由常规的地震剖面推导出地下地层的波阻抗 密度 速度等信息 为勘探开发提供重要的依据 而石油勘探开发的形式在不断地发 展 计算机技术也不断地提高 使得地震资料反演技术成为了当前热门的研究 课题 其中波阻抗反演 Acoustic Impedance AI 是反演技术中很重要的一个 部分 它利用地震资料反演反射系数 进而得到地层波阻抗信息的技术 同时 地震资料反演技术本身具有的不适定性 应用条件的限制以及储层 的复杂性 常规叠后地震资料品质不高且分辨率低等因素 一定程度上影响了 该技术更广泛的应用 1 2 波阻抗反演研究现状 地震反演的基本目的是利用地震波在地下介质中的传播规律 通过数据采 集 处理与解释等流程 推测地下岩层结构和物性参数的空间分布 为勘探开 基于迭代再加权最小二乘的地震资料稀疏反演方法 2 发提供重要的依据 1 3 波阻抗反演是地震反演的一个重要组成部分 由于波阻 抗信息是联系地质和地球物理的一座桥梁 叠后计算数据量相对较少 在实际 生产中应用方便而且效果明显 因此 波阻抗反演在地震反演中具有特殊的地 位 侠义的地震反演概念指的就是波阻抗反演 4 波阻抗反演是利用地震资料反演地层波阻抗 或速度 的地震特殊处理解 释技术 它产生于 20 世纪 70 年代 80 年代开始得到蓬勃发展 1983 年 Cooke 介绍了地震资料广义线性反演方法 从而揭开了波阻抗反演技术的新篇 章 5 周竹生等人在 90 年代初期提出了综合利用地质 地震和测井资料进行约 束反演 克服单一的线性反演方法的缺陷 6 90 年代中期 李宏兵在国内提出 了一条走递推反演与宽带约束反演结合的道路 此方法为我们解决了从单道出 发的反演方法不能在跟不上消除噪音的困惑 7 在此基础之上 有人进行了无 井多道反演和有井多道反演的研究 8 9 使波阻抗反演方法更加完善 90 年代至 今 围绕一维波阻抗反演的各类算法及应用成果层出不穷 随着研究的开展 在 1997 年左右开始出现了一些反思的文章 指出了波阻抗反演中存在的一些缺 陷 并提出了一些解决方案 10 12 BP Amoco 公司的 Connolly 在 1999 年正式发 表了弹性波阻抗反演方法的论文 该方法最早是在 1993 1994 年由 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 为了提高反演的精度 第一章 引 言 3 和收敛速度 针对非线性反演方法 地球物理工作者将各种最优化方法引入到 了波阻抗反演中 各种各样的混合优化算法也陆续出现 力求在解决计算精度 的同时提高算法收敛速度 减少反演多解性 18 尽管已经作了大量的研究工作 波阻抗反演方法仍然存在许多的不足 如算法复杂 实现困难 抗干扰能力差 从而导致分辨率降低 收敛速度慢等 因此仍需要进一步研究和改进这些优 化反演算法 在理论指导下 国内外相继研制出一些商业化程度高的波阻抗反 演方法及软件 如 STRATR GLOG JASON 等 19 目前国内使用的波阻抗反 演软件基本上都是引进软件 国产软件较少 96 年以前的软件一般都是单井控 制下的地震资料反演 只能解决构造简单没有断层或岩性 岩相横向变化小的 问题 在多井条件下和地质条件复杂的地区难以推广使用 国外从 80 年代末提 出了利用声波测井资料作为约束条件 对过井地震剖面进行正反演联合迭代求 取地下波阻抗的方法 将反演方法推向非线性 这种新方法利用了测井资料的 高频信息 因而表现了强劲的发展势头 目前成为各国外石油公司重点发展的 对象 如美国 HGS 公司宽带约束反演 BCI 技术 法国 CGG 公司波阻抗反演模 拟 ROVIM 技术 俄罗斯的 PARM 加拿大的 STRATA 丹麦的 ISIS 以及荷兰 推出的 JASON 技术 今后地震反演的技术 是在现在的基础上 充分利用测井 地震 地质的信息 继续探索先进的算法 使储层横向预测更符合实际 这些 方法在实际中取得了一定的应用效果 但是还不能够完全满足实际需要 不同 反演方法结果差别较大 不同人使用同一种方法结果也有差别 往往导致错误 的地质解释 究其原因 除了地震资料品质 子波提取不准确 垂直入射假设 与实际有误差等因素外 还有另外两个方面的主要影响因素 其一 由于地震 数据具有带限性质 直接反演只能获得波阻抗模型中的中频成分 缺少应有的 低频和高频成分 其二 Hardmard 意义下 波阻抗反演是一个不适定问题 其 反演结果会出现不稳定性和多解性 目前 解决带限问题最有效的手段是对解 加以约束进而进行反演 如测井数据的约束 而解决不适定问题主要使用正则 化方法并辅之以适当的最优化技巧 15 通常利用测井资料在地质知识约束的插 值方法 来弥补地震资料中仍缺乏的低频和高频信息 1 3 稀疏反演方法的研究情况 线性问题的稀疏解在很多领域都有研究 线性问题可以分为正定问题 欠 基于迭代再加权最小二乘的地震资料稀疏反演方法 4 定问题和超定问题 正定问题的解存在并且唯一 欠定问题存在无穷多解 因 此一般研究其最简单的解 即稀疏解 超定问题不存在严格的解 一般求出满 足最小二乘意义的近似解 当然也可先验约束研究其稀疏近似解 目前的稀疏 算法主要针对欠定的线性问题求解 在实际应用中求稀疏解的最常用方法是匹配追踪类方法 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 从介绍波阻抗反演的反褶积分析入手 分析并导出了基于迭代加权的 波阻抗反演的基本问题 并介绍了迭代方法以及迭代的效果 1 5 论文的主要创新点 本文对迭代波阻抗反演的反褶积问题进行了进一步的分析 从而导出了基 于迭代再加权的波阻抗反演方法 并且应用了几种迭代方法来解方程组 通过 理论模拟和实际资料处理证明了新方法在收敛速度上的优越性 1 对波阻抗反演方法进行了系统分类 对稀疏反演方法以及与之相关的 匹配追踪和压缩反演进行了论述 2 从反褶积问题入手 建立了一种基于迭代加权的波阻抗反演方法 分 析了加权因子的选取 并研究了 Gauss Seidel 迭代法 共轭梯度法和预条件加 速技巧 3 对以上方法 用不同噪音水平的地震数据进行了反演计算 计算结果 验证了算法的有效性 第二章 波阻抗与稀疏反演概述 5 第二章 波阻抗与稀疏反演概述 2 1 基本概念 2 1 1 常规波阻抗反演技术的基本假设前提 目前常用的波阻抗反演软件所采用的方法基本上都是在褶积模型的基础上 建立的 因此 要求资料都要满足褶积模型的假设前提 可以概括为以下四个 方面 1 地震模型 假设地层是水平层状介质 地震波为平面波法向入射 其地震 剖面为正入射剖面 并假设地震道为地震子波与地层反射系数的褶积 2 反射系数序列 普通递推反演中 假设反射系数为完全随机的序列 在约 束稀疏脉冲反演中 假设反射系数为由一系列大的反射系数叠加在高斯分布的 小反射系数的背景上构成 3 噪音分量 通常假设波阻抗反演输入的地震数据的振幅信息反映了地下波 阻抗变化的情况 地震剖面没有多次波和绕射波等噪音分量 因此在资料处理 时 可以考虑的处理流程是反褶积 噪音剔正 尤其是多次波 处理的最终目 标是得到高分辨率的保幅波场 4 地震子波 假设反射系数剖面中的每一道都可以看作地下反射系数与一个 零相位子波的褶积 实践情况下地震子波往往是混合相位或最小相位的 因此 需要对地震剖面进行相位校正处理 2 1 2 波阻抗反演方法的基本原理 波阻抗 波在某种介质中传播的速度与介质密度的乘积 在本论文中其符 号和定义为 Z 表示波阻抗 表示密度 V 表示波速 根据弹性波理VZ 论 当一个纵波入射到波阻抗不同的两种介质的分界面时 一般要产生反射纵 波 反射横波 透过纵波 透过横波四种新的波动 一般主要研究纵波反射波 在水平层状介质 地震波为平面波法向入射的假设下 问题会得到简化 只考 虑产生的反射纵波和透射纵波 反射横波和透射横波的情况 地震子波 震源附近地震波以冲击波的形式传播 当传播到一定距离时 基于迭代再加权最小二乘的地震资料稀疏反演方法 6 波形逐渐稳定 此时的地震波被称为地震子波 进行反演时 求取地震子波的 方法 一般有统计法和参数法 统计法是从实际地震记录中统计出地震子波的振幅谱 和相位谱 从而生成地震子波 参数法是给出某种地震子波类型的参数而生成 地震子波 如 Ricker 子波 Ricker 子波是最常用的零相位子波 在时间域表 示为 2 1 1 2 222 21 tv M M etvtf 其中 是主频 零相位子波分辨力较高 最有利于解释 M v 由定义可见 只有层状介质的波阻抗存在差异时反射系数才不为 0 当子 波遇到波阻抗界面后 会发生反射 从而被接收到 由于是按照时间序列排列 的 以下分别记为接收到的地震数据 子波 反射系数按照时间 trtwtD 的序列 并已经过离散化 则按照定义 可以得到 1 1 2 1 LtrLwtrwtrwtb L 表示子波的长度 数学上定义此运算为褶积 下记 为褶积运算符 则上 式可以表示为 2 1 2 rwb 因此 地震记录实际上是由反射系数与子波的褶积得到的 如果将子波作 为输入信号 地下介质作为一个线性时不变的 系统 或者滤波器 接收信 号为信号通过系统后的输出 则该褶积可以看做是一个系统的响应过程 陆基 孟 1993 该响应不但包括一次反射还包括多次放射 假定震源波形在向地下 传播的过程中不变 因此褶积模型不包括固有的吸收 地震记录可以看作是由 震源子波与地下反射率函数 多次反射 仪器等诸多因素的褶积的过程 虽然 可以给合成地震记录加上随机噪音 但是这里的用来建立反褶积滤波器的褶积 模型暂时不包括随机噪音 第二章 波阻抗与稀疏反演概述 7 图 2 1 模拟的反射系数 子波和合成的地震记录 如果能求得一个算子是子波的逆 即 则它与地震 tb tw 1 twty 记录褶积时 可由 2 1 式得 tD 2 1 1 trtrttrtwtwtbty 3 即地震记录和子波的逆进行褶积可以得到反射系数 这个过程叫做反褶积 tD 由以上过程可以看到 反射系数序列中含有波阻抗随时间变化的信息 这 就提供了速度和密度随时间变化的信息 因此可得到地层 岩性及构造在地下 空间分布的信息 在有利条件下还可得到岩石孔隙率 渗透率 孔隙流体性质 油 气 水 乃至地层压力的信息 波阻抗反演需要求出反射系数 因而其 反演就是已知地震记录 由提取的子波反褶积得到反射系数序列 tb tw tr 在使用炸药作为震源时 震源实际是一个 脉冲 由于地层的滤波作用 震源信号在地下传播过程中高频被衰减 能量被吸收 变成一个子波 所以在 地震勘探中观测的信号其峰值频率往往小于 100Hz 这就导致了记录信号的分 辨率降低 记录到的反射是减弱了强度的宽脉冲 并且通常随着反射时间的增 加信号越宽和越弱 子波在地层中传播 携带着反射系数序列这种有用的地质 基于迭代再加权最小二乘的地震资料稀疏反演方法 8 信息返回地面 只有通过反褶积消除子波才能恢复反射系数序列 同时反褶积 可以压缩地震子波 压制交混回响和短周期多次波 从而提高时间分辨率 反 褶积可应用于叠前资料 也可应用于叠后资料 理想的反褶积应该压缩子波并 消除多次波 在地震道内只留下地层反射系数 子波压缩可以通过将反滤波器 作为反褶积算子来实现 反滤波器可以将地震子波转变成尖脉冲 每个脉冲的 强弱与界面的反射系数的大小成正比 脉冲的极性反映界面反射系数的符号 当应用于地震合成记录时 反滤波输出应为地层脉冲响应 反射系数 根据弹性波理论 当一个纵波入射到波阻抗不同的两种介质分 界面时 一般要产生反射纵波 反射横波 透过纵波 透过横波四种新的波动 在实际生产工作中 主要研究纵波反射波 界面倾角不大时 可近似认为是垂 直入射 在垂直入射的情况下 将会使问题简化 纵波入射时我们将只考虑产 生的反射纵波和透射纵波 反射横波和透射横波的情况 这时界面的反射系数 定义为 2 1 4 1122 1122 12 12 VV VV ZZ ZZ A A R R F 其中 AF AR分别表示反射振幅和入射振幅 而 分别对应 21 Z Z 21 21 V V 上层和下层的波阻抗 密度 速度 当平面波非垂直入射时 反射系数将随着 入射角的变化而变化 这种变化与界面两边介质的波阻抗和各种弹性参数有关 特别是与介质的泊松比有很密切的关系 Zoppritz 方程 入射角等于零是它的 特例 褶积模型 根据介质模型假设不同 波阻抗反演分为基于波动方程的反演 和基于褶积模型的反演两大类 前者由于算法结构复杂 计算量大 抗干扰差 很难获得一个稳定解 在实际中未得到广泛应用 目前主要使用基于褶积模型 的直接反演方法 它算法简单 对地震噪音敏感性小 但一般情况下并不一定 能得到一个稳定的解 1 褶积 在时域中为了便于求得线性时不变系统的零状态响应 如果能够将 任意激励信号分解为单元信号 并使每一单元信号在系统中产生的零状态响应 易于求得 那么根据系统的线性时不变特性 就可以利用叠加原理很容易求得 原信号在系统中产生的零状态响应 系统对单位冲击函数的零状态响应为 t 对于线性时不变系统 当激励是强度为的延时冲击函数时 thk 0 ttk 第二章 波阻抗与稀疏反演概述 9 其零状态响应为 任意信号可表示为 0 ttkh te 2 1 5 dtetktttkete k t lim 0 作用下系统响应为 te tr 2 1 6 dthetktthtketr k t lim 0 其中为积分变量 为参变量 上式褶积积分又可记作 t 2 1 7 thtetr 2 模型产生 水平层状介质模型是共中心点 CMP 叠加方法的基础 在 CMP 叠加中 对分享同一炮检中点的一些地震道进行正常时差校正求和 产生 一个逼近 1D 层状介质垂直入射平面波响应的求和道 本着这种数据 形成了一 种地震道模型 假定从波阻抗不连续界面反射 或透射 的子波与入射子波具 有同样的波形 一个地震道顺序记录这样一系列连续子波 也就记录了波阻抗 不连续面 地震道可以简单地看成是每一个独立的反射界面叠加的结果 这里 每个反射界面都认为是一个滤波器 地震记录可认为是震源子波与介质反射系 数 即地下垂直入射反射系数序列 的褶积 1 时间域表达 连续情形 2 1 8 tedtrwtetrtwtb 其中为地震子波 为地震反射序列 为地震记录 为随机噪音 wrde 矩阵 向量情形 2 1 eWrb 9 其中为地震反射系数向量 为噪音列向量 n rrrr 10 e 为地震波记录向量 为子波阵 子波长度为 10n bbbb WL 基于迭代再加权最小二乘的地震资料稀疏反演方法 10 L L L L L L L nLn w w ww www www www www ww ww w W 2 12 12 12 12 12 12 1 1 2 频率域表达 连续情形 2 1 erWb 10 各表示傅立叶变化下的 地震记录 子波 反射系数 噪音 bWre 2 1 3 波阻抗反演的分类 波阻抗信息是联系地质和地球物理的重要媒介 其叠后计算数据量相对较 少 在实际生产中应用方便 效果明显 因此通常的地震反演概念指的就是波 阻抗反演 将地震剖面转换成波阻抗剖面 不仅便于解释人员将地震资料与测 井资料连接对比 而且能有效地对储层物性参数的变化进行研究 从而得到物 性参数在空间上的分布规律 对油气的勘探和开发有很重要的指导意义 对于 波阻抗反演的重要性 一个很经典的论述是著名的地球物理学家李庆忠院士所 说的 交到地质人员手中的地震资料应是作了反演的波阻抗剖面 波阻抗反 演是高分辨率地震资料处理的最终表达形式 25 26 而常常被人们引用的美国 阿莫克石油公司在墨西哥湾海上寻找非背斜油藏不断取得成功的例子 29 也印证 了这一说法 波阻抗反演是目前应用最广泛的储层预测和油藏描述手段之一 其发展日 新月异 计算方法多种多样 这些方法由于其基于的数学物理模型不同 各有 优缺点 利用的基础资料和能够达到的反演效果是不同的 波阻抗反演方法有 多种 Seislog Velog Delog PIVT BCI ROVIM SLIM PARM 和 Strata 14 概括 起来有两大类 基于反射系数逆公式的直接反演和基于正演模型的迭代反演 迭代反演又可以分为无井条件下的广义性反演 GLI 和有井条件下的宽带约束 第二章 波阻抗与稀疏反演概述 11 反演 BCI 33 还有一些属于非线性反演范畴的方法 如混沌理论 神经网 络等 1 依据算法的实现方式分为 基于反褶积的反演 包括道积分 递归 广义线性 稀疏脉冲反演等 基于波动方程的反演 包括 Born 反散射反演等 基于随机过程的反演 包括随机反演 随机模拟 模拟退火反演等 基于特征 分析的反演 包括特征反演 神经网络反演等 基于动力学特征的反演 包括 混沌反演等 2 按测井和地震的相对作用 无井反演 包括道积分等 地震为主的测 井约束类反演 包括递归反演 广义线性 宽带约束反演 稀疏脉冲 模拟退 火反演等 测井为主的统计类反演 包括随机模拟 随机反演等 井 震约束联 合反演 包括特征反演 神经网络反演等 这种分类方案体现了测井资料对反 演约束作用的强弱 便于研究人员根据测井和地震基础资料评价各种反演方法 的反演效果 3 依数学算法不同可分为 线性反演 包括道积分 递归反演 广义线 性 宽带约束反演 稀疏脉冲等 非线性反演 包括随机模拟 随机反演等 2 2 稀疏反演方法 常规反褶积得到的地震资料分辨率低 由于地震子波仅有有限宽度的频谱 而反射系数序列具有无限宽度的频谱 使得利用地震数据求取反射系数序列具 有高度的不适定性 特别是具有无穷多个解能够拟合地震数据 此外地震数据 往往含有随机噪声 会使得地震反演具有不稳定性 为了提高地震数据的分辨 率 更客观地反映地下介质的分布规律 需从地震数据中可靠地消除地震子波 的影响 这就要求对地震反演问题进行正则化处理或加入先验知识加以约束 考虑到地球介质沉积过程的成层特征 可合理假设反射系数序列具有一定的稀 疏性 而稀疏性正是各沉积地层厚度不同特征的数学表征 即利用等时间间隔 表示反射系数序列时 反射系数序列中某些时刻的反射系数肯定为零 利用反 射系数序列稀疏性 实现地震数据的反演 可在地震数据含有噪声条件下 有 效地消除地震子波的影响 提高地震数据分辨能力 提高对地下介质真实分布 情况的刻画能力 可以说 地震数据的稀疏反演构成当前提高地震分辨能力的 主要途径 基于迭代再加权最小二乘的地震资料稀疏反演方法 12 2 2 1 匹配追踪 常规反褶积得到的地震资料分辨率低 子波及反射系数准确度不高 因此 为了提高分辨率 要求反射系数能够稀疏化 但同时 稀疏化后的反射系数再 经过褶积过程会与原始数据产生偏差 导致信噪比降低 因此 需要寻求一种 方法 在提高分辨率的同时信噪比不会显著降低 使分离的子波及反射系数更 接近真实情况 对于稀疏解的求取 一类常用的方法是匹配追踪算法 Matching pursuit MP Davis et al 1997 Mallat et al 1993 对于常见的信号处理过程 通常 可以表达为 bAr 的形式 其中是 A 一个变换矩阵 A 的每个列向量称为一个原子 b 是时nm 空域中的信号 r 是在变换域中的系数 是求解的对象 通常 A 是欠定的 所 以存在无穷多满足上式 匹配追踪法作为一种贪婪算法 能够迭代的求出最稀 疏的 因为 b 表示为 A 的列向量的线性组合 所以最稀疏的表示方式是唯一的 令是 A的列向量 将上式展开得 niAi 2 1 nnx AxAxAb 2211 匹配追踪法的每次迭代包括两步 1 寻找当前剩余的相关最大的列 然后求出当前列的系数 2 升级剩余 即由当前的剩余 减去选出的列 得到新的剩余 首先令初始剩余 令初始 r 为 0 向量 在第 k 步迭代 选出与当前bR 0 剩余相关绝对值最大的列 即满足 1 k R iknik ARA maxarg 11 之后求出第 k 步的剩余 求出新原子的系数 kkkkk AARRR 11 然后继续下一步的迭代 由此可见 匹配追踪法算法简单 kkkk Arrr 1 易行 缺点是循环量大 即使某个原子被选出来 但是其系数没有经过优化 所以 在以后的迭代中还会选出这个原子 这样就使得迭代步数增加 匹配追踪类方法还包括正交匹配追踪方法及其变形 Davis et 第二章 波阻抗与稀疏反演概述 13 al 1997 Donoho et al 2006 Tropp 2004 正交匹配追踪法比匹配追踪法 多一个正交化的过程 一旦选出上述过程中要求的原子 就执行正交化过程求 出最优系数 以后的迭代就不会再选这个原子 因此比匹配追踪法效率要高 匹配追踪方法保证每次迭代使得解更加逼近真实解 如果解非常稀疏 则 匹配追踪类算法的效率很高 因此对于解严格稀疏或者维数很小的问题匹配追 踪类方法效率较高 在存在适当噪音的情况下匹配追踪类方法也仍然适用 Donoho et al 2006 2 2 2 稀疏反演和压缩感知 求解线性方程组的稀疏解最直观的数学表示为如下优化问题 2 2 bAr 1 其中代表非零元素的个数 问题 2 10 是一个组合优化问题 不存在多项 0 式解法 需要计算所有的解才能找到最稀疏的解 但是这样做是不现实的 在 实现上意义不大 Candes et al 2005 一般采用以下问题求稀疏解 Chen et al 1998 2 2 2 1 min bArts r 2 其中代表绝对值 l1范数 该问题便是基追踪问题 Basis pursuit BP 1 是问题 2 10 的松弛 并且在一定的情况下 问题 2 2 2 的解就是问题 2 2 1 的解 这样就把一个组合优化问题变为一个连续的凸优化问题求解 基于 l1范数的算法是稀疏优化方法的主流 也正是本文求解的问题 问题 2 2 2 的更一般的情形是假设采样的数据存在噪音 即 n nArD 因此相应的问题变为 2 2 2 1 min bArts r 3 基于迭代再加权最小二乘的地震资料稀疏反演方法 14 其中是非负的实参数 用来控制噪音的大小 问题 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 可以变为无约束优化问题 2 2 2 21 minDArr 4 问题 2 2 4 是一个惩罚的最小二乘问题 其目标函数是拉格朗日函数 其中 称为拉格朗日乘子 根据对偶理论 当选取参数和合适时 问题 2 2 3 和 2 2 4 是等价的 问题 2 2 4 还存在另一个关系紧密的问题 2 2 1 2 2 min rts bAr 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 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 改为二次规划并且用预 第二章 波阻抗与稀疏反演概述 15 条件共轭梯度法求解 文 Wang et al 2005 Wang et al 2009b 利用内点法 来解遥感中的稀疏优化问题 以上问题都可以归结为模型 qp ll 2 2 0 min qprbAr q q p p 6 可见问题 2 2 2 是问题 2 2 6 的特例 Wang et al 2009a 稀疏优化算法种类繁多 但是主要是基于 l1范数的凸优化方法 沿用上文记号 褶积的矩阵向量乘积形式为 2 2 bnAr 7 其中是子波延拓成的矩阵 为反射系数 为噪音 为地震记录 假设ArnD 噪音是白噪声 则该方程组的求解可以变为最优化问题 2 2 2 minbAr 8 假设地震子波已知或反射系数序列 可以准确估计 通过褶积模型 这wr 样也可直接估计出反射系数序列 使得地震记录与褶积结果达到最大程度的r 拟合 这里假设噪声近似为零 在时域中反褶积的目的是使残差最小 但是由 于要求反射系数是稀疏的 因此需要在模型中增加稀疏约束或采用其它策略求 出稀疏解 若反射系数是稀疏的 通过引进 l0范数可以将解方程问题变为一个组合优 化问题 0 min r 2 2 ts bAr 9 其中用来衡量噪音水平的大小 这个组合优化问题在一定条件下可以变为一个等价的 l1范数优化问题求解 1 minr 2 2 ts bAr 基于迭代再加权最小二乘的地震资料稀疏反演方法 16 10 但是由于目标函数不是连续可微的 因此需要变形后采用内点法 投影梯度法 同伦法等方法求解 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 n i i Cauchy r fbAr 1 2 min 11 其中为正则参数 越大 解越稀疏 刘喜武等 2006 针对问题 2 2 11 基于预条件共轭梯法来求反射系数序列 提高了运算速度 降低了内存使 用量 此外基于模型 同样可以通过可以采用基于共轭梯 22 mindrbAr 度法的迭代重赋权最小二乘法得到稀疏解 Sacchi 1995 本文采用的方法 是将 2 2 12 1 min r ts 2 bAr 等价的转化为 2 2 13 2 21 minbArr 反射系数的稀疏化遇到的最直接的一个问题就是当原始数据不满足 Nyquist 采样定理时是否能够得到完整的反射系数 受到实际条件的约束 采样时间和 采样个数经常受到限制 采集的数据量往往不满足采样定理 因此需要打破采 样定理的限制 如果采样过程是线性的 则采样过程可以写为 2 2 14 bAr 其中是 r 原始数据 b 为采样的数据 A 是采样矩阵 地震数据采集中 当 第二章 波阻抗与稀疏反演概述 17 采样个数满足采样定理时 则 若采样的数据是不完整的 则该IA br 问题成为一个欠定问题 因此是病态的反问题 而压缩感知理论可以由采样不足的数据重建出原始数据 该理论基于两个 假设条件 第一个条件是原始数据在某个变换下是稀疏的 即假设存在变换使得 其中 s 的非零元素的个数很少 于是可以先求出稀疏的 s 再恢复原始信sr 号 即我们先解以下问题 2 2 15 Bsb 其中 然后求出 AB sr 第二个条件是必须满足约束等距性 restricted isometry property Candes et al 2005 如果满足 B 约束等距性 且采样个数满足一定条件时 则稀疏解可 以求出 如果以上两个条件满足 我们就可以由上式求出稀疏解 s Candes 2006 压缩感知理论给出了采样不足时稀疏反褶积的理论保证 总之 对于压缩感知与稀疏反演而言 对于含噪声的数据处理时 均可将 待优化的目标函数定义为下式 2 2 16 0 min 2 2 2 1 1 1 qpxxbAx q q q q p p 对于 p q 取值不同的上述优化问题 人们已开展大量的研究工作 特别是 压缩感知理论与技术的出现 直接推动着计算科学 信息技术 通讯技术和编 码与解码技术等多个领域的发展 可以说压缩感知已成为当前最热点的研究领 域 就比较成熟的优化算法而言 主要存在两种优化技术或策略 即将上述优 化问题转化为等价的有约束优化和无约束优化问题 有约束优化的主要算法是 原始 对偶内点法 无约束优化的主要算法是迭代再加权方法 这两种方法各有 不同优势与不足 但无论在研究领域和实际应用领域中 始终处于相互竞争与 相互促进的过程