1、 武汉 大学 硕士学位论文 一 种新的频率选取策略及其在频域波 形反演 中 的应用 研 究 生 姓 名： 学 号： 指导教师姓名、职称： 专 业 名 称 ： 微电子 学与固体电子学 研 究 方 向： 信息 处理与 微 系统 二 XX 年 X 月 摘 要 近 年来，全波形反演已经 成 为估计地下介质参数的有力工具，其具体 实现 方法分为频域和时域两 大 类。频 域 波形反演方法相对于时域反演方法而言，数据的选取更加灵活 ，考 虑到地震数据在频域表示中的冗 余 性， 可 知前 者 具有更高的计算效率。 频域全 波形反演是一种能够获得高分辨率反演结果的反演方法， 但 在模型较大时进行全频 点 反演的
2、计算量 很 大 。 Sirgue 和 Pratt (2004)提出 一种波形反演中的频率点选取的方法，在不影响反演精度的前提下，可以 减 少频域波 形 反演所 需 计算的频率点数 目 ，从而 节省 计算量。 目前 ，大多 数的 频率点选取方法都是基于地下介质速度为均匀分布这一 假设前提的，为了研究在更加接近实际地下介质的背景模型下的情形，本文在线性递增速度模型的基础上，利用 散射 波 场 的 Born 近似 、 WKBJ 近似 以 及地震波的一般传播规律 ， 根据 模型修正量中 竖直 方向上的波数覆盖，提出一种新的频率点选取 策略。 结果表明，在接近实际的地下介质速度分布情况下，可以采用比 S
3、irgue 和 Pratt 方法更少的频点来进行频域反演，而不会影响反演的精度，从而可以进一步提高反演计算的效率。本文针对 一 维速度模型、 Marmousi 模型和 Overthrust 三 个 模型进行反演实验，模拟算例结果验证了该优化的频率点选 取方法的有效性。 关键词 ： 频率点选取 ；波数覆盖 ； 线性 递增模型 ； 频 域 波形反演 ABSTRACT In recent years, waveform inversion has become a powerful tool for estimating the parameters of the underground mediu
4、m, and its implementation is divided into two major categories of frequency domain and time domain method. Compared to time domain inversion method, the frequency domain inversion method has more flexible data selection, considering the redundancy of seismic data in the frequency domain, show that t
5、he frequency domain inversion method has higher computational efficiency. The method of frequency-domain full-waveform inversion is capable of obtaining high resolution results, but it suffered a heavy computation when inverting large-scale models with all frequencies. Sirgue and Pratt (2004) propos
6、ed a discretization strategy that reduces the input frequencies required for inversion and increases efficiency without loss accuracy of the inversion results. At present, the most frequency selection methods are based on the assumption that the velocity of underground medium follows uniform distrib
7、ution, in order to study the situation that the background is more close to the actual underground medium, in this paper we use the velocity of medium linearly increases with depth, apply the Born approximation of scattered field, the WKBJ approximation, and the general law of seismic wave propagati
8、on, according to the coverage of vertical direction wave number in model correction, and proposed a new frequency selection strategy. Experiment results showed that under the more realistic assumption, our method could obtain the same level of accuracy with fewer frequencies than that required by Si
9、rgue and Pratts strategy. In this paper we conducted experiments with 1D model, Marmousi model and Overthrust model, and the results demonstrated the validity of our method. Key words: Frequency selection; Wave-number coverage; Linear increasing velocity; Frequency domain waveform inversion 目 录 摘 要
10、ABSTRACT 1 绪论 1.1 研究 背景和意义 1.2 国 内外发展现状 1.3 主 要研究内容 2 频 域波形反 演 算法简介 2.1 波 形反演的基本步骤 2.2 频 域正演模拟 2.2.1 二 维 频 域波动方程 2.2.2 有 限差分表示 2.2.3 PML 层 边界处理方法 2.3 频 域 波 形反演 3 频 率点的选取策略 3.1 散 射波场的 Born 近似 及 Green 函数求 解 3.2 频 域波形反演中的波数覆盖原理 3.3 线性 递增速度介质中地震波的传播 3.4 频率 点的选取方法 4 实验 验证 4.1 一 维模型合成数据反演 4.2 二 维模型合成数据反演
11、4.2.1 Marmousi 模型 合成数据反演 4.2.2 Overthrust 模型 合成数据反演 5 总结 与展望 参考 文献 1 绪论 随着 人们 对 地下非均匀介质认识的不断深入以及计算机性能的大幅度提升，地 震 波反 演成像技术也相应地得到了非常迅猛地发展。然而 ，在 使用传统频率点选取 策略所得频率点 进 行 频率 从低到高的迭代频域 波形 反演 的 前提下， 在 勘探的目标区域尺度非常大 的情况下 ， 计算 效率就会成 为此 种反演 方法所 面临的挑战。此时 ， 基于 一种新的梯度图像 中 波数覆盖 的 频率点选取策略可能 会 为 地 震波形反演 带 来新的发 展 。 本章 将
12、回顾 相关 的研究背景及现状， 简述 传统的频率点选取策略的原 理 ，阐述它在实际应用中的局限性及其改进的 可 能性，给出本文的总体研究思路和主要研究内容。 1.1 研究 背景和意义 地 震勘探 (seismic exploration)是 地球物理勘探中的一种重要方法，其主要目标 是 使用人工震源激发的 地 震波来定位 地 下介质中的矿藏 (如 石油，天然 气 ，煤 ，矿石，矿物，地热能 源 ，水资源等 )， 获得地下的 岩 石特性 、碳 氢化合物 属性 等地质信息。在 陆 地 地 震勘探中，由地表或近地表的人工震源激发 的地 震波 经由地下介质向地球内部传播，当遇到 地 下介质的 物 性参
13、 数快速 变化 区域 时 地 震波 便 会产生散射波 (包括 反射波，折射波，绕射波等 )， 一部分 返 回到设置在地 表的 高灵敏 度检波器的 散射 波会被此检波器记录下来，所有的这些被不同检波器(hydrophone)记录的 地 震记录构成地震数据集； 同 理，在海洋地震勘探中， 返回到 水面 下方 沿 拖曳 电 缆方向 设置 的 高 灵 敏 度检波器 处 的散射波会被此检波器记录下来， 所 有 这 些不同检波器的散射波记录构成了地震数据集。 对于上述 两种勘探情形下所得的 地 震数据集， 其中 主要记录了地震波传播的旅行时信息 以 及振幅响应信息 。通过 分析地震数据集，人们一方面想建立
14、一种从地震数据集到地下结构概貌的映射 ， 这一过程通常被称作是地震成像或地 震 偏移成像；另一方面想通过 地 震数据集来对目标区域中未知的物理参数 (如 速度，弹性模量，密度，阻抗等 )进行适当估算 ， 这一过程通常被称作是地 震 反演；当然，人们更希望可以通过分析地 震数据 集一次性地获得目标区域的结构概貌和物性参 数 。 由于 地震数据集的复杂性， 要 想完全正确地实现上述的想法依然是目标我们所面临的困难。 地 震数据集中 包括 了所有检波器所能接收的散射波的能量信息， 每 一个能量信息代表 了 从震源到检波器 的地 震 波 所特有的散射 传播 方式，根据散射波在地下介质中的具体传播方式可
15、以将其分为单次反射 波 和多次反射 波 。根据 微 扰理论，可以将地 震 波的单次散射 波表 示为 目标 区域物性参 数扰 动的线性 映射 ，相应地多次反射则是一个非线 性 映射。 考虑 到非线性 反 问题相对于 线性 问题的复杂性，目前大多数 地 震成像和地震反演算法都是基于数 据中只包含地震波单次散射波的假设。 由于 地震数据在频 域 中存在 信息 冗余 (Mulder and Plessix 2004)，因此可以使用所有频率 点 的一个子集来表示 整体 信息 ，考虑 到频 域中 声波方程的有限 差 分的计算量 与 频率点数目成正比 (Plessix et al. 2001)，利用此结论
16、就 可以 将时 间域 中 的波形反演转化到频率域中从而提高 反 演的计 算效率 。 就 从低频率点到高频率点的频率迭代反演来说，其计算用时 线性 正比于 波形反 演 时 的频率点 数 ，因此如何高效 地去除 地 震数据集中的冗余频率点，也即如何 合理 地选择最好 (或最少 )的 频率点组合 ， 就成了 提高 此时计算效率的关键。 现 在常用的方法 有 ， (1)基于均匀介质假设的前提下采用 固定 间隔频率点 ； (2)基于均匀介质假设 的 前提 下 使用 间隔随频率值线性递增的频率点； (3)通过 速度值随深度线性递增模型中地 震 波的最大传播深度 以 及 多 尺度反演 中频率 点 对局部 极
17、点的影响来对方法 (2)中 的线性系数进行简单的修正。 鉴于实际 的地下介质往往是非均匀 的 复杂 情 形， 以 及方法 (3)只 是对方法 (2)的 简单修正，因此我们有必要 研 究一种新的频率点选取策略。本论文 旨在充分借鉴现有的频率点选取策略，发现它们的不 足 之处，提出 一 种 完全基于速度值随深度递增的线性速度模型下 的 新的频率点选取方法 ，为 下一步的研究打下良好的基础。 1.2 国内 外发展现状 人们 对于 地 震数据的处理方法一般分为地 震偏移 成像和地震反演，其中 偏移成像的主要目的是 通过 恢复光滑背景模型上 的 不连续区域来 构建 地下结构概貌，而地震反演的目的则是定量
18、地恢复地下介质的物理参数。 如果 基于最小二乘法的失配函数 (error functional)的成像 方法 所 得的黑塞 (Hessian)矩阵 是一个对角 (或 近对角 )阵 ， 那么此 失配函数的最小值可以通过 对 梯度矩阵前乘黑塞 矩阵 (或伪(pseudo-)黑塞 矩阵 )的逆 矩阵 来 获得 ，此 时的成像方法 可 以用 于 表示 真振幅 偏移成像 (Mulder and Plessix 2004; Plessix and Mulder 2004)，考虑到基于 最小 二 乘法的频域波形反演的情形，我们可以发现地震成像和地震反演有时并没有明确 的 界线 (Tarantola 1986
19、)。 由于 地震数据 集 在频域 表示 时存 在 信息冗余 (Mulder and Plessix 2004; Plessix et al. 2001)， 因此频 域 波形反演中可以使用 全 频率点的一个子集 来提高 反演的效率。 Mulder 和 Plessix (2004)指出在 频率点选 择 时需要考虑两个主要准则： (1)由等价 海森 堡测不准原理 (Heisenbergs principle)所 确定的偏移图像中的 反射 体 分辨 率 ； (2)频域奈 奎斯特抽样定理 (Nyquists theorem)，考虑 到偏移 包括 在给定背景模型上的地 震波 传播 以 及多个震源 检 波器
20、的求和， 因 此尽管有时不符合奈奎斯特抽样定理但依然可以获得一个没有混叠的结果。 在 均匀背景 的 假设下，利用上述的两准则，得到了一个 频率点 间隔为固定值的频率点选取策略。 Sirgue 和 Pratt (2004)根据 频 域最 速 下 降反演方法， 在 均匀介质假设的情 况 下，利用 偏移成像和 频 域波形反演的 等 价性 (Tarantola 1986)， 借助 平面波近似及 波 恩近似 (Born approximation)原理 ， 揭示 了利用多震源检波 距 地震数据进行频域反演时在垂直方向上波数覆盖存在冗余的现象， 据此 提出了一种适用此情形下的频率点间隔随频率线性递增 的
21、选取策略。 同 年， Yokota 和 Matsushima (2004)在 均匀介质假设基础上， 根据观 测数据与正演数据之间的相 位 误差大于 时就 会在结果上出现 跳 周 现象 (Cycle Skipping， 即 一 个异常 高的 渡 越时间记录 )的情形，进而 提出一种相对应的频率点选取策略。 Mulder 和 Plessix (2008)利用速度值 随深度线性 递增模型 中地 震 波传播的最大 深度 以及 多 尺度反演 中频率 点 对局部 极点的 影响 (Bunks et al. 1995)， 对 Sirgue 和 Pratt (2004)所 提出的频率点选取策略 中 的线性系数进
22、行 了 简单修正。 Anagaw 和 Sacchi (2014)以 及 Kim 等 (2011)对 Sirgue 和Pratt (2004)的方法所得的频率点进行简单 地 分组 ，从而 验证不同 的 频率点组合对反演的影响。 目前 来说， 对 于频率点选取策略的研究还 不 是 很 多，前面 所 述的方法 基本 上都是基于均匀背景 (或 只是对均匀背景简单修正 )的 情形 。 随 着 反演理论的发展以及 勘探 的需要，人们自然 要 问是否存在基于相 对于 均匀背景而言更一般的 频率 点选取策略。目前 ，该 课 题 的研究还处于初 级 阶段，对于从事反演研究的学者来说它是一个具有 一 定挑战性的研
23、究方 向 。 1.3 主 要 研究内容 在认真 分析总结 现 有频率点选取策略的基于基础上，本 文结合傅 里 叶 变换 的意义，首先建立 基于 一阶 Born 近似的 散射波场积分表达式 ；然后 在 WKBJ 近似的情形下获得 Green 函数 的解析表达式 ； 在此基础上结合频域波形反演算法通过梯度图像中的波数覆盖情况 推导出 线性递增速度背景模型下 相邻 两个频率点之间的关系 式 ，进而得到新 的 频率点选取策略。 本论文 各章节的内容安排如下： 第 一 章 介绍 了本 论文 的研究目的及意义 ，对相关 研究现状进行 总结 ， 并简要 介绍了本论文的研究内容及安排。 第 二章 回顾 了频
24、域 波形 反 演的基于原理和方法，并对预处理梯度法进行了简要介绍。 第 三 章 针对 线性 递增速度 背景 模型 ，提出 一种基于 散射 波场 Born 近似 及WKBJ 近似 的 新的频率点选取策略，并 在不 同源检 距 组合情形下，推导 了 反 演 中相邻两个频率点之间的关系。 第 四章 对于本 文提出的新的频率点选取策略，在不同的 背景 速度模型下与Sirgue 和 Pratt (2004)提出 的方法进行 模型 测试。详细 对比 两种频率点选取策略下反演 结果 的分辨率及计算用时。 第 五章 总结 本论文的相关研究结果，分析研究中存在的问题，并展望下一步研究的工作重点。 2 频 域波形
25、反演 算法简介 地 震波形反演是一个非线性问题， 为 了降低问题的求解难度，通常会使用基于梯度的线性处理方法来对其进行求 解 。 常用 了 求 解方法有：全牛顿法，高斯牛顿 法，最速下降法， 预 处理梯度法等。 在实际 使用中，由于 黑 塞 矩阵 的尺寸及其较高的计算消耗，一 般 很 难 对 其 进行 直接 计算 (Mulder and Plessix 2004; Pratt et al. 1998; Tarantola 1984)。通常 ， 最速 下降法 具有 较低的收敛速度 (Tarantola 1986)，考虑 到 黑 塞矩阵的 求 解难度， 我 们可以 通过使用 伪黑塞矩阵 预 处理
26、来 加速 梯度 法的收敛 过程 (Mulder and Plessix 2004; Plessix and Mulder 2004)。本章以 频域反演为前提，主要介绍预处理梯度 反 演方法，为后 续 的研究奠定基础。 2.1 波 形反演的基本步骤 目前 地震成像中绝大多数的偏移成像方法都是采用线性化简化， 其思想为 对散射波场进行线性化处理 (也 即 Born 近似 )，然后 将 Born 近似 的积分形式自然地与广义的 Radon 变换 结合 起 来，进而通过逆 广义 Radon 变换来求取 近似地反演解，进而 恢复地 下介质结构。 就 其中的数学部分 而言， Radon 已 于 1917
27、年 首次证明了利用投影函数可以 恢复 出满足条件的唯一的图像函数，因此利用上述的偏移成像方法可以唯一地重建地下介质结构，以达到 地 震勘探的 最 终目的。 而 我们这里要介绍的 波形 反演方法 也遵从 上面所述的理论， 它们都 是 通过 建立 勘探 所得的实际数据与用理论计算得到的模 拟 数据的数据残 差 ，并 将 此残差 与实际地下模型与先验地下模型之间的扰动 之间 形成 映射 ，通过不断修正先验 地 下模型来达到最终获得一个在其上 的模拟数据与勘探数据拟合 度最好 的 地下模型(Anagaw and Sacchi 2014; Gauthier et al. 1986; Mora 1987;
28、 Pratt 1999; Tarantola 1984)。 原理 虽然相同，但是对于具体的实现方法之间则有较 大的区分， 为 了论文的需要这里仅就波形反 演 的方法进行讨论，其实现过程基本可 分 为以下几个部分： (1)在 某一 地 理区域设置勘探的震源及相应的检波器 (也 即接收点 )。 (2)采集 各个震源所对 应 的检波 器 上的波形记录， 形成 实测 (观 测 )地震记录。由 于 实测地震记录是波形反演中的 表 征地下结构的最重要的数据， 故实测 地震记录的 质量显著 地影响了波形反演的效果及可信程度。 (3)根据 实测地震记录中震源及 检波器 的 坐标 位置 以 及 检波 器中的 地
29、 震道数据 来 建立初始 介质 模型， 也 即先验模型。 在 通常的基于梯度的 局部 最优化 反 演方法中，如果先验模型与实际勘探 区域 的 实际 模型 差别过 大，则在反演过程中 陷 入局部极 小 点的可能性大 大 增加 ，在 实际 的 应 用 过 程 中一般可以 使用 走时反演的反演结果作为 初始 模型，而对于基于全局优化的方法则可以减弱对初始模型的 先验限制。 (4)在初始介质 模型上按实际勘探时的震源及检 波 器的设置情 形 进行其上的正演模拟， 并 记录检波器处的波形信息，即预测 (正演 )地 震记录。 (5)把正演 地震记录与观测地震记录 按 每个 检 波器做差得到地震记录残差，
30、由波 恩 近似 可以看出地震记录残差为 先 验模型与目标 模型之间的模型扰动的映射，由 Radon 变换 可知 可用 此 地 震记录残差来 解得 模型的 扰 动也即恢复地下 介质 结构。 将 地震记录残差表示成列矩阵的形式， 其中 的每个元素代表相应的检波器上的 地 震记录的残差， 把 地震残差 矩阵按 一定方法映射到一个标量目标函数 (如 F范数 )。 (6)在 目标函数的形式上， 使用 泛函 分 析的相关理论 获得 目 标 函数极小时所对应 的相 对于初始 介质 模型的模型修正量，用此修正量修正初始模型。把 修正 后的初始 模型当作先验 介质 模型 重 复步骤 (2)-(6)，直到修正 后的初始模型 上 的目 标 函数满足设定的 阈 值时 停止对介质 模型的迭代更新。 概括 起来波形反 演 的基本步骤 可 分为：获取观测地震记录，设置初始 介质 模型，在初始 介质 模型 上 正演获得正演地震记录， 对 正演地震记录与观 测 地 震 记录做差 并 将结果 表示 成 波 场残差矩阵， 然后将 波场残 差 矩阵映射成标 量 目标函数，再通过 泛 函 方法 获得 介质 模型修正量并 按 此对 介质 模型进行更新 ， 具体的流程如图 2.1 所 示。
