全波形反演

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