部分6-全波形反演
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
将其转换到频率域，并对空间求导进行有限差分离散可得到如下的线性方程 组：
SP = F
其中P为离散之后的频域声压波场；S为阻抗矩阵，其为频率、介质参数、 离散网格大小的函数；F为非齐次震源项。
二、频率-空间域全波形反演
若阻抗矩阵的逆可以显式写出，则可得：
P = S -1F
对于多炮情况下，只需求取一次阻抗矩阵的逆便可通过简单的矩阵相乘得
f t
it
F e
it
d
i i t t
F i e
d
F i e
e d f t et
引入复频率后的波形反演即为阻尼波场的波形反演，或称为LaplaceFourier域的波形反演。 物理解释： 1、时域解释：在时间方向进行阻尼，在叠前记录上看，起到了加窗的效果， 减弱了不同时刻同相轴的串扰；在反演结果上看，起到了层剥离的效果，由浅 及深的反演使得反演过程更稳定。 2、频域解释：时间阻尼作为一时变滤波器，改变了原始数据的频率响应，使 得相位谱更趋于线性，振幅谱更平滑，减少了陷入局部极小值的可能性。
这里：
1 1 hmax / z
2
k0
v0
为最大半偏移距； 为目标层深度。
hmax
z
频率选取的准则为低频的最大波数等于 高频的最小波数。
二、频率-空间域全波形反演
Marmousi模型反演结果：
二、频率-空间域全波形反演
5、Laplace-Fourier域全波形反演与近地表速度 Laplace-Fourier域波形反演，对于单一 一个阻尼参数来说，可以看做是一种 层剥离的效果。通过控制不同时间上的波场的权重来控制所能反演的深度， 以此来减少由于不同波形的串扰引起的非线性性。 近地表速度在波形反演中占有重要的作用。表层速度由于影响到所有记录中 的所有波形，固其精确与否直接影响随后波形反演的稳定性。
二、频率-空间域全波形反演
带预处理的最速下降法：
mk 1 mk k (diagHa I)1ςmmO m
其中：
k 为模型参数更新步长；
ς m为平滑算子，由模型 平滑约束推出；
mO m 为梯度向量；
diagHa 为Hessian矩阵的对角线元素。
(diagHa I)1 称为预处理算子，主要校正几何扩散的影响， 进行照明补偿等。
迭代算法选取 计算量过大 牛顿类算法：Gauss-Newton、Full-Newton等， 高效的---计算能力能承受; Hessian矩阵的计算 相位编码、数据压缩等
高精度--- 成像分辨率高
反演流程选取 …… 多尺度反演策略、逐层反演等
一、前 言
全波形反演理论推导的两种方法：
First Optimize Then Discretize (OTD) 涉及到泛函分析等复杂的数学推导； 时域反演理论的推导用的这种方式。 First DiscretizeThen Optimize (DTO) 首先将问题离散化之后然后利用相应的优化理论进行推导，主要 涉及线性代数的相关知识，相对简单易懂； 频率-空间域波形反演的理论推导这样利用这种方式。 对应这两种方法，下面分别讨论频-空域全波形反演与时间域全波形反演。
下面首先讨论阻尼波场的层剥离效果，然后结合阻尼波场减少由于表层速度 不准确引起的反演的不稳定性。
二、频率-空间域全波形反演
二、频率-空间域全波形反演
2、数值实现
有限差分格式：
频 -空域数值模拟利用混合交错网格（交错网格 +旋
转交错网格）。该方法在保证计算精度的情况下使
得阻抗矩阵带宽尽可能小，以提高计算效率。
混合网格差分格式
频散曲线
阻抗矩阵基本形态
二、频率-空间域全波形反演
边界条件： 选取PML或CPML边界条件。在频率域PML边界条件的实现较为简单。 线性方程组求解： 对于 2D情况及小尺度的 3D情况下，选用直接求解方法，如多波前并行 稀 疏 直 接 求 解 法 （ MUltifrontal Massively Parallel Sparse direct Solver）。 对于大尺度 3D情况，一种解决策略为利用迭代解法或者混合求解方法 求解大规模线性方程组；另一种解决策略为通过时域正演记录提取波场 的频域响应。
二、频率-空间域全波形反演
频率-空间域波形反演最终结果
反演速度与真实速度对比
二、频率-空间域全波形反演
反演中间过程中每个频率的最终结果
二、频率-空间域全波形反演
4、Marmousi模型波形反演测试 如右下图所示，修改后的Marmousi模型。主要的修改为对模型进行抽稀且只 截取了一部分，然后在表层加一水层。炮点和检波点位置如图中白线所示。 所用初始模型如左下图所示。
一、前 言
存在的关键问题：
初始模型建立 目标函数选取 主要目标：
旅行时层析、偏移速度分析（成像域）、Laplace域波形反演等 不同测度： L2模、 L1 模、混合模等 可靠的 --有效的避开局部极小值； 不同类型数据残差：相关型、反褶积型、对数型等
稳定的---对数据要求不是很苛刻； 梯度类算法：最速源自文库降、共轭梯度等，预处理方法
diagH a diag J t J * diag H a i
p k k 1 mi
n
pk m i
对比梯度向量的计算可知，即使是 Hessian对角线元素的计算也不可避免需
要显式的计算偏导数波场。对于每个频率每次激发，都需要m+1次的正演过程。
R F S
t
1
R
梯度的物理解释：
波场剩余量的逆时传播与正向波场的零延迟互相关；
隐式求解敏感性矩阵的计算量： 对于一次激发只需两次正演模拟，计算量与模型参数的个数不再是正比关系。
二、频率-空间域全波形反演
Hessian矩阵对角线元素的计算： Hessian矩阵的对角线元素为：
其中n为检波点的个数，m为介质参数的个数。 对频率-空间域离散后的波动方程相对于介质参数进行求导得到：
SP F P S S P mi mi mi mi P S 1f (i ) , mi f (i ) = S P mi
偏导数波场物理解释： 入射波场经过地下散射点散射后接收到的波场； 梯度的物理解释： 偏导数波场与数据剩余量的零延迟互相关； 显式求解敏感性矩阵的计算量： m个模型参数需要m+1次正演；
由于地面地震检波点与炮点位置大量重复， 其个数远小于地下网格点数，利用波场的互易性， 不再计算地下 mi 到地面各个点的波场，而是计 算地面点xj到地下各个点的波场。由于虚震源的 个数明显减少，将明显的提高计算效率。 mi
xj
二、频率-空间域全波形反演
阻尼波场的应用（Laplace-Fourier域全波形反演）： 对于如下形式的傅里叶变换： 引入复频率可得：
二、频率-空间域全波形反演
梯度向量的快速求解 将敏感性矩阵显式写出为：
J [S1f 1 S1f 2
将其代入梯度向量表达式可得：
S1f m ] S1F
t O t 1 m O J R S F R m
F S
t
1 t
15HZ单频快照
二、频率-空间域全波形反演
3、简单模型波形反演测试
速度模型如右图所示；背景 速度场为3500m/s，中间异 常体为3700m/s。炮点与检 波点位置如右图所示。 给定背景 速度 场 ，对 随 机 选取的几个频率进行频率 空间域速 度反 演 ，反 演 的 目标为中 间异 常体 ， 反 演 的频率分别为：4HZ、5HZ、 7HZ 、 10HZ 、 13HZ 、 16HZ、20HZ。
二、频率-空间域全波形反演
反演所用频率分别为： 0.5,1.2,2.5,5.6,12.5 ，这五个频率是利用 Laurent
Sirgue（2004）在1D介质假设下提出的公式计算得到，如右下图所示。 在一维介质假设下，可反演得到的最小和最大空间波数为：
k z min 2k0 k z max 2k0
主要利用一次波
一、前 言
尽管全波形反演在成像精度方面具有不可替代的优势，但
是由于问题的非线性性、数据的不完整性（低频缺失、噪
声）、数学物理模型描述实际问题能力有限性等问题，使 得反演问题复杂性显著增加，较易收敛至局部极小值。 目前，全波形反演方法还存在着像计算量过大、不稳定、 依赖于初始模型等各式各样的问题。
G(x r1 , x sm ) G(x r 2 , x sm ) G(x rn , x sm )
二、频率-空间域全波形反演
基于优化理论的波形反演方法 定义数据剩余量为：
R P d
定义L2模目标函数为：
1 O m = R t R * min 2
其中m为模型参数向量，可以选为声波速度、慢度、压缩模量及其变形等。 模型参数化对随后的反演过程有重要的影响，同时也是反演之前要做的第一 步。 对于局部寻优的梯度类算法，推导的重点是目标函数相对于模型参数的一阶 导数（梯度向量）及二阶导数（Hessian矩阵）。
提 纲
一、前言
二、频率-空间域全波形反演
1、基本原理 2、数值实现 3、模型试算
4、小结
三、时间域全波形反演
四、结论与认识
二、频率-空间域全波形反演
1、基本原理
问题的离散化： 考虑时域声波方程：
2 p x, t 1 p x, t 1 p x, t s x, t 2 2 v x t x x z z 1
到多炮正演记录:
P P s1 P s2 S-1 F s1 F s2
阻抗矩阵逆的每一列为格林函数：
-1 P sn = S F
F sn
G(x r1 , x s1 ) G(x r1 , x s2 ) G(x , x ) G(x , x ) r2 s1 r2 s2 -1 S = G(x rn , x s1 ) G(x rn , x s1 )
二、频率-空间域全波形反演
3、模型试算
1、均匀介质模型频域正演
均匀介质模型：速度为 4000m/s，采样点 为 101×201 ，采样间隔为 40m ，震源点 位于（ x=1000,z=500 ）处，输入子波为 雷克子波。
二、频率-空间域全波形反演
2、Marmousi 模型频域正演试算
为验证正演算法在复杂介质情况下 的稳定性，对 Marmousi 模型进行 测试。 5HZ单频快照
由最速下降法可以看出梯度向量与 Hessian矩阵的计算在反演中的重要作 用，下面重点推导这两者相应的表达式。
二、频率-空间域全波形反演
梯度向量的计算：
m O
O J t R m
其中J为Frechét 求导矩阵（敏感性矩阵），其元素为：
J ij pi , m j i 1, 2,..., n; j 1, 2,..., m
反演方法相比还存在相应的不足。
一、前 言
全波形反演介质参数成像方法
原始波动方程
一次反射波、多次波、 折射波等波形信息
传统介质参数成像方法 低波数
模型的各个波 数成分
缺少中波数信息
高波数
介质参数分解
波动方程线性化
m m0 m
偏移成像 宏观速度模型建立 （旅行时）信息
叠加速度分析； 偏移速度分析； 层析成像；
提
纲
1、叠前时间域速度分析 2、共成像点道集提取 3、叠前深度域偏移速度分析 4、弹性矢量波联合层析速度反演 5、全波形反演
提 纲
一、前言 二、频率-空间域全波形反演 三、时间域全波形反演 四、结论与认识
一、前 言
全波形反演算法最早于上世纪八十年代由 Tarantola 等提出，
当时限于计算能力、野外数据采集方式及质量、数值反演算法 等方面的限制，未能得到大范围的发展。在这个阶段取而代之 是基于 Born 或 Rytov 近似的背景速度场估计（宏观速度模型 建立）及模型扰动量成像（偏移成像）。 虽然传统的成像方法得到了广泛的应用，但在理论上与全波形
SP = F
其中P为离散之后的频域声压波场；S为阻抗矩阵，其为频率、介质参数、 离散网格大小的函数；F为非齐次震源项。
二、频率-空间域全波形反演
若阻抗矩阵的逆可以显式写出，则可得：
P = S -1F
对于多炮情况下，只需求取一次阻抗矩阵的逆便可通过简单的矩阵相乘得
f t
it
F e
it
d
i i t t
F i e
d
F i e
e d f t et
引入复频率后的波形反演即为阻尼波场的波形反演，或称为LaplaceFourier域的波形反演。 物理解释： 1、时域解释：在时间方向进行阻尼，在叠前记录上看，起到了加窗的效果， 减弱了不同时刻同相轴的串扰；在反演结果上看，起到了层剥离的效果，由浅 及深的反演使得反演过程更稳定。 2、频域解释：时间阻尼作为一时变滤波器，改变了原始数据的频率响应，使 得相位谱更趋于线性，振幅谱更平滑，减少了陷入局部极小值的可能性。
这里：
1 1 hmax / z
2
k0
v0
为最大半偏移距； 为目标层深度。
hmax
z
频率选取的准则为低频的最大波数等于 高频的最小波数。
二、频率-空间域全波形反演
Marmousi模型反演结果：
二、频率-空间域全波形反演
5、Laplace-Fourier域全波形反演与近地表速度 Laplace-Fourier域波形反演，对于单一 一个阻尼参数来说，可以看做是一种 层剥离的效果。通过控制不同时间上的波场的权重来控制所能反演的深度， 以此来减少由于不同波形的串扰引起的非线性性。 近地表速度在波形反演中占有重要的作用。表层速度由于影响到所有记录中 的所有波形，固其精确与否直接影响随后波形反演的稳定性。
二、频率-空间域全波形反演
带预处理的最速下降法：
mk 1 mk k (diagHa I)1ςmmO m
其中：
k 为模型参数更新步长；
ς m为平滑算子，由模型 平滑约束推出；
mO m 为梯度向量；
diagHa 为Hessian矩阵的对角线元素。
(diagHa I)1 称为预处理算子，主要校正几何扩散的影响， 进行照明补偿等。
迭代算法选取 计算量过大 牛顿类算法：Gauss-Newton、Full-Newton等， 高效的---计算能力能承受; Hessian矩阵的计算 相位编码、数据压缩等
高精度--- 成像分辨率高
反演流程选取 …… 多尺度反演策略、逐层反演等
一、前 言
全波形反演理论推导的两种方法：
First Optimize Then Discretize (OTD) 涉及到泛函分析等复杂的数学推导； 时域反演理论的推导用的这种方式。 First DiscretizeThen Optimize (DTO) 首先将问题离散化之后然后利用相应的优化理论进行推导，主要 涉及线性代数的相关知识，相对简单易懂； 频率-空间域波形反演的理论推导这样利用这种方式。 对应这两种方法，下面分别讨论频-空域全波形反演与时间域全波形反演。
下面首先讨论阻尼波场的层剥离效果，然后结合阻尼波场减少由于表层速度 不准确引起的反演的不稳定性。
二、频率-空间域全波形反演
二、频率-空间域全波形反演
2、数值实现
有限差分格式：
频 -空域数值模拟利用混合交错网格（交错网格 +旋
转交错网格）。该方法在保证计算精度的情况下使
得阻抗矩阵带宽尽可能小，以提高计算效率。
混合网格差分格式
频散曲线
阻抗矩阵基本形态
二、频率-空间域全波形反演
边界条件： 选取PML或CPML边界条件。在频率域PML边界条件的实现较为简单。 线性方程组求解： 对于 2D情况及小尺度的 3D情况下，选用直接求解方法，如多波前并行 稀 疏 直 接 求 解 法 （ MUltifrontal Massively Parallel Sparse direct Solver）。 对于大尺度 3D情况，一种解决策略为利用迭代解法或者混合求解方法 求解大规模线性方程组；另一种解决策略为通过时域正演记录提取波场 的频域响应。
二、频率-空间域全波形反演
频率-空间域波形反演最终结果
反演速度与真实速度对比
二、频率-空间域全波形反演
反演中间过程中每个频率的最终结果
二、频率-空间域全波形反演
4、Marmousi模型波形反演测试 如右下图所示，修改后的Marmousi模型。主要的修改为对模型进行抽稀且只 截取了一部分，然后在表层加一水层。炮点和检波点位置如图中白线所示。 所用初始模型如左下图所示。
一、前 言
存在的关键问题：
初始模型建立 目标函数选取 主要目标：
旅行时层析、偏移速度分析（成像域）、Laplace域波形反演等 不同测度： L2模、 L1 模、混合模等 可靠的 --有效的避开局部极小值； 不同类型数据残差：相关型、反褶积型、对数型等
稳定的---对数据要求不是很苛刻； 梯度类算法：最速源自文库降、共轭梯度等，预处理方法
diagH a diag J t J * diag H a i
p k k 1 mi
n
pk m i
对比梯度向量的计算可知，即使是 Hessian对角线元素的计算也不可避免需
要显式的计算偏导数波场。对于每个频率每次激发，都需要m+1次的正演过程。
R F S
t
1
R
梯度的物理解释：
波场剩余量的逆时传播与正向波场的零延迟互相关；
隐式求解敏感性矩阵的计算量： 对于一次激发只需两次正演模拟，计算量与模型参数的个数不再是正比关系。
二、频率-空间域全波形反演
Hessian矩阵对角线元素的计算： Hessian矩阵的对角线元素为：
其中n为检波点的个数，m为介质参数的个数。 对频率-空间域离散后的波动方程相对于介质参数进行求导得到：
SP F P S S P mi mi mi mi P S 1f (i ) , mi f (i ) = S P mi
偏导数波场物理解释： 入射波场经过地下散射点散射后接收到的波场； 梯度的物理解释： 偏导数波场与数据剩余量的零延迟互相关； 显式求解敏感性矩阵的计算量： m个模型参数需要m+1次正演；
由于地面地震检波点与炮点位置大量重复， 其个数远小于地下网格点数，利用波场的互易性， 不再计算地下 mi 到地面各个点的波场，而是计 算地面点xj到地下各个点的波场。由于虚震源的 个数明显减少，将明显的提高计算效率。 mi
xj
二、频率-空间域全波形反演
阻尼波场的应用（Laplace-Fourier域全波形反演）： 对于如下形式的傅里叶变换： 引入复频率可得：
二、频率-空间域全波形反演
梯度向量的快速求解 将敏感性矩阵显式写出为：
J [S1f 1 S1f 2
将其代入梯度向量表达式可得：
S1f m ] S1F
t O t 1 m O J R S F R m
F S
t
1 t
15HZ单频快照
二、频率-空间域全波形反演
3、简单模型波形反演测试
速度模型如右图所示；背景 速度场为3500m/s，中间异 常体为3700m/s。炮点与检 波点位置如右图所示。 给定背景 速度 场 ，对 随 机 选取的几个频率进行频率 空间域速 度反 演 ，反 演 的 目标为中 间异 常体 ， 反 演 的频率分别为：4HZ、5HZ、 7HZ 、 10HZ 、 13HZ 、 16HZ、20HZ。
二、频率-空间域全波形反演
反演所用频率分别为： 0.5,1.2,2.5,5.6,12.5 ，这五个频率是利用 Laurent
Sirgue（2004）在1D介质假设下提出的公式计算得到，如右下图所示。 在一维介质假设下，可反演得到的最小和最大空间波数为：
k z min 2k0 k z max 2k0
主要利用一次波
一、前 言
尽管全波形反演在成像精度方面具有不可替代的优势，但
是由于问题的非线性性、数据的不完整性（低频缺失、噪
声）、数学物理模型描述实际问题能力有限性等问题，使 得反演问题复杂性显著增加，较易收敛至局部极小值。 目前，全波形反演方法还存在着像计算量过大、不稳定、 依赖于初始模型等各式各样的问题。
G(x r1 , x sm ) G(x r 2 , x sm ) G(x rn , x sm )
二、频率-空间域全波形反演
基于优化理论的波形反演方法 定义数据剩余量为：
R P d
定义L2模目标函数为：
1 O m = R t R * min 2
其中m为模型参数向量，可以选为声波速度、慢度、压缩模量及其变形等。 模型参数化对随后的反演过程有重要的影响，同时也是反演之前要做的第一 步。 对于局部寻优的梯度类算法，推导的重点是目标函数相对于模型参数的一阶 导数（梯度向量）及二阶导数（Hessian矩阵）。
提 纲
一、前言
二、频率-空间域全波形反演
1、基本原理 2、数值实现 3、模型试算
4、小结
三、时间域全波形反演
四、结论与认识
二、频率-空间域全波形反演
1、基本原理
问题的离散化： 考虑时域声波方程：
2 p x, t 1 p x, t 1 p x, t s x, t 2 2 v x t x x z z 1
到多炮正演记录:
P P s1 P s2 S-1 F s1 F s2
阻抗矩阵逆的每一列为格林函数：
-1 P sn = S F
F sn
G(x r1 , x s1 ) G(x r1 , x s2 ) G(x , x ) G(x , x ) r2 s1 r2 s2 -1 S = G(x rn , x s1 ) G(x rn , x s1 )
二、频率-空间域全波形反演
3、模型试算
1、均匀介质模型频域正演
均匀介质模型：速度为 4000m/s，采样点 为 101×201 ，采样间隔为 40m ，震源点 位于（ x=1000,z=500 ）处，输入子波为 雷克子波。
二、频率-空间域全波形反演
2、Marmousi 模型频域正演试算
为验证正演算法在复杂介质情况下 的稳定性，对 Marmousi 模型进行 测试。 5HZ单频快照
由最速下降法可以看出梯度向量与 Hessian矩阵的计算在反演中的重要作 用，下面重点推导这两者相应的表达式。
二、频率-空间域全波形反演
梯度向量的计算：
m O
O J t R m
其中J为Frechét 求导矩阵（敏感性矩阵），其元素为：
J ij pi , m j i 1, 2,..., n; j 1, 2,..., m
反演方法相比还存在相应的不足。
一、前 言
全波形反演介质参数成像方法
原始波动方程
一次反射波、多次波、 折射波等波形信息
传统介质参数成像方法 低波数
模型的各个波 数成分
缺少中波数信息
高波数
介质参数分解
波动方程线性化
m m0 m
偏移成像 宏观速度模型建立 （旅行时）信息
叠加速度分析； 偏移速度分析； 层析成像；
提
纲
1、叠前时间域速度分析 2、共成像点道集提取 3、叠前深度域偏移速度分析 4、弹性矢量波联合层析速度反演 5、全波形反演
提 纲
一、前言 二、频率-空间域全波形反演 三、时间域全波形反演 四、结论与认识
一、前 言
全波形反演算法最早于上世纪八十年代由 Tarantola 等提出，
当时限于计算能力、野外数据采集方式及质量、数值反演算法 等方面的限制，未能得到大范围的发展。在这个阶段取而代之 是基于 Born 或 Rytov 近似的背景速度场估计（宏观速度模型 建立）及模型扰动量成像（偏移成像）。 虽然传统的成像方法得到了广泛的应用，但在理论上与全波形