第一章 前言

反问题

在野外地震勘探中，我们采集地震数据d的过程实际上是一个正演过程，即：

$${\rm \textbf d}={\rm \textbf {Lm}}$$

其中，$\rm \textbf L$为隐式地依赖模型$\rm \textbf m$的非线性正演算子，其能够通过模型$\rm \textbf m$预测数据$\rm \textbf d$。如果原始的正演算子是非线性的，那么它通常可以通过背景模型来线性化。此时，我们可以将地球物理反问题定义为利用记录的数据矢量$\rm \textbf d$反演模型矢量$\rm \textbf m$。

据Hadamard(1902)，适定数学模型，例如线性化后的${\rm \textbf d}={\rm \textbf {Lm}}$，具有以下特点：

1. 存在解。对于地球物理数据而言，噪声较多引起方程不相同或超定会导致解不存在。一个解决办法是寻找能够使得数据误差和惩罚项的加权和最小的解。
2. 解是唯一的。常规的地球物理勘探中，检波器和震源只布设在有限区域内，这样会导致$\rm \textbf L$存在非空零空间，进而使得$\rm \textbf L^T \rm \textbf T$有零特征值。我们通常使用正则化来缓解这一问题。
3. 解持续地依赖数据，否则解将会不稳定。例如，当数据中噪声水平产生微小变化时，肯恩贵引起模型的不连续变化（当病态矩阵$\rm \textbf L$存在许多接近于零的特征值时，会使得许多模型都能拟合相同的数据）。

我们通常会改写反问题以使我们能够找到一个最优模型，这个最优模型能够使数据误差和惩罚方程的加权和最小。

旅行时成像 Traveltime tomography

旅行时成像是通过拾取的特定同相轴的旅行时来估计地下速度分布。旅行时层析需要六个步骤，这六个步骤对于大多数其他地球物理反演问题也是是适用的。

1. 模型$\rm \textbf m$的离散。地球慢度模型被离散化为一系列未知慢度的网格，每一个网格点上都有一个未知的慢度值。假设网格非常细所以在第$j$个网格中慢度$s_j$为常数。$N$个未知慢度构成了$N \times1$的模型矢量$\rm \textbf m$。
2. 数据$\rm \textbf d$的离散。对于初至旅行时成像，我们首先要拾取每一道的初至并记为第$i_{th}$射线的时间$t_i$。$M$道数据能够形成$M\times1$的数据矢量$\rm \textbf d$。
3. 模型算子$\rm \textbf L$的离散。通过给定一个初始模型和基于正演算子$\rm \textbf L$生成的预测数据来解决反问题。旅行时成像是基于高频近似假设的，因此第$i_{th}$射线的旅行时是所有单元格中旅行时$d_i$的和。 $$\underset{j} \Sigma l_{ij}m_j=d_i \rm \space or \space \rm \textbf {Lm}=\rm \textbf d \tag{1.4}$$ 其中，$l_{ij}$是第$i$条射线在第$j$个单元格中的射线路径段长度，$m_j$是第$j$个单元格中的常慢度$s_j$。这里，$\rm \textbf L$是$M\times N$的矩阵。
4. 线性化。方程$1.4$是非线性的，因为射线路径段长度$l_{ij}$取决于速度模型。这并不奇怪，因为根据$Snell$定律，较大的速度梯度变化将会引起射线路径的较大变化。所以$L$隐式地依赖模型$\rm \textbf m$，所以除非我们的初始慢度模型非常接近实际模型，否则将无法反演出好的结果。解决办法是我们从一个接近于真实模型的模型$m^{(0)}=m_o$出发，对模型与数据之间的关系进行线性化。然后利用数据反演出一个更准确的模型$m^{(1)}$，将其作为新的初始模型，反复迭代直到收敛。
   定义第$j$个模型参数为$m_j$，通过将在$\rm \textbf m$下观测到的数据$d_i(\rm \textbf m)$在初始模型$\rm \textbf m_o$处泰勒一阶展开来将其线性化： $$d_i(\rm \textbf m) \approx d_i (\rm \textbf m_o) + \underset {j} \Sigma [ \frac { \partial d_i (\rm \textbf m)} {\partial m_j}]_{\rm \textbf m_o} \delta m_j, \tag{1.5}$$ 其中，$\rm \textbf m=\rm \textbf m_o+\delta \rm \textbf m$。所以有 $$\delta d_i(\rm \textbf m)= \underset {j} \Sigma [\frac{\partial d_i(\rm \textbf m)} {\partial m_j}] {\rm \textbf m_o} \delta m_j, \tag{1.6}$$ 或者写成矩阵矢量形式 $$\delta \rm \textbf d(\rm \textbf m)=\rm \textbf L(\delta \rm \textbf m),\tag{1.7}$$ 这里，数据残差$\delta d_i=[d_i(\rm \textbf m)-d_i(\rm \textbf m_o)]$为观测数据第$i$个分量$\rm \textbf d_i(\rm \textbf m)$与预测数据第$i$个分量$\rm \textbf d_i(\rm \textbf m_o)$的差，模型扰动$\delta \rm \textbf m=\rm \textbf m-\rm \textbf m_o$为实际模型$\rm \textbf m$与猜测模型$\rm \textbf m_o$的误差。矩阵$\rm \textbf L$为雅可比矩阵，其元素为$l_{ij} = \frac{\partial d_i (\rm \textbf m_o) } {\partial m_j}$，即$Fr\acute echet$导数，表示第$i$个数据相对于模型扰动在第$j$个单元格处的敏感核。
   对于旅行时成像和较小的单元格，$Fr\acute echet$导数可以表示为射线段长度$l_{ij}$在第$j$个单元格第$i$条射线，公式$(1.6)$变为 $$\delta t_i=\Sigma_j l_{ij}\delta s_j \tag {1.8}$$
5. 正则化解。记录数据含有的噪声会导致一组不一致的(inconsistent)（当我们使用的模拟算子无法模拟数据所需要的所有的物理信息时，也会出现不一致问题）超定方程。此外，由于许多模型都可以满足相同的数据，所以解可能是不稳定的。
   为了一定程度上解决这种问题，我们通常寻找能够使得目标函数$\epsilon$最小的模型，目标函数由一个惩罚项和数据残差的$p$范数的$p$次方构成： $$\epsilon=\frac{1}{p} ||\rm \textbf L\delta \rm \textbf m-\delta \rm \textbf d||^p_p+ \eta^2g(\rm \textbf m), \tag {1.9}$$ 其中$\eta^2$为一个非常小的正数，$g(\rm \textbf m)$是惩罚因子项，当估计模型接近实际模型的先验估计时其会变得很小。$p$是正整数，残差矢量$\rm \textbf L \delta \rm \textbf m-\delta \rm \textbf d$是观测数据和预测数据之差。当$p=2$时，残差矢量的平方长度即为误差函数。
   有时惩罚项会被表示为$g(\rm \textbf m)=\frac{1}{p} ||m-m^{‘}||^p_p$（模型残差），其中$m^{‘}$为先验模型。通常设置$p=2$，即平方数据残差和约束模型惩罚之和$||\rm \textbf L \delta \rm \textbf m -\delta \rm \textbf d ||^2+\eta ^2||\delta \rm \textbf m||^2$作为目标函数。这种情况下，正则化阻尼最小二乘解为: $$\delta \rm \textbf m=[\rm \textbf L^{\rm \textbf T}\rm \textbf T + \eta^2\rm \textbf I]^{-1}\rm \textbf L ^{\rm \textbf T}\delta \rm \textbf d, \tag {1.10}$$ 其中$I$为单位矩阵，当元素均为复数时转置算子表示复共轭。对于复数，我们使用伴随符号$\dagger$代替转置符号。
6. 迭代正则化解。数据与模型之间的关系通常是非线性的，所以公式$1.10$的解常由迭代的形式来得到: $$\rm \textbf m^{(k+1)}=\rm \textbf m^{(k)}-\alpha [\rm \textbf L^{\rm \textbf T}\rm \textbf T + \eta^2\rm \textbf I]^{-1}\rm \textbf L ^{\rm \textbf T}\delta \rm \textbf d^{(k)},\tag {1.11}$$ 其中，$\alpha$为步长，$\rm \textbf m^{(k)}$（$\delta \rm \textbf d^{(k)}$）表示第$k$次迭代的模型（数据残差）。与矩阵$\rm \textbf L$相关的射线是由第$k$次的速度模型计算得到的，不断迭代得到后续的速度模型知道数据残差到我们能够接受的水平。矩阵$[\rm \textbf L^T\rm \textbf L]$的求取、存储和反演成本巨大，因此我们通常使用其主对角元素$[\rm \textbf L^T\rm \textbf L]_{ij} ≈[\rm \textbf L^T\rm \textbf L]_{ii}\delta _{ij}$来获得预处理后的最速下降梯度解。 $$\delta m_i^{(k+1)}=\delta m_i^{(k)}-\frac{\alpha [\rm \textbf L ^{\rm \textbf T}\delta \rm \textbf d^{(k)}]} {[\rm \textbf L^T \rm \textbf L]_{ii}+\eta^2},\tag {1.12}$$ 非线性地震反演的主要问题是目标函数常常包含许多局部小值；所以，迭代解通常会陷入局部极小并且永远无法达到全局最优，或者无法反演到真实模型。为了避免这一问题，我们通常在前面的迭代中只反演重要的特征（骨骼化数据，Luo and Schuster，1991a；1991b），即反演长波长信息，后续的反演中慢慢添加高阶信息。
   常见的骨骼化的方法包括multi-scale反演，early-arrivals反演，near-offset反演等。在后续的反演中逐渐添加远偏移距、长记录、以及更复杂的物理信息。