第五章 射线路径走时层析

射线路径走时层析

正如第1章所讨论的，地球物理走时层析成像是一种通过反演指定体波的走时来获得地球速度分布的方法。从20世纪70年代初期开始，它已成为地震成像地球内部的主要手段(Aki et al., 1977; Dziewonski and Anderson, 1984; Humphreys et al.,1984; Bishop et al., 1985; Nolet, 1987; Iyer, 1989; Zelt and Smith, 1992; Iyer and Hirahara, 1993)。之前大多数采用射线追踪的高频近似实现走时层析成像，但近来趋势是使用有限频率的层析成像(Harlan, 1990; Luo and Schuster, 1991a; Woodward, 1992; Marquering et al., 1999; Montelli et al., 2004)，其中走时残差沿着波路径（wavepaths）而不是射线路径（raypaths）进行反投影。有限频率波动方程层析成像，也被称为波动方程层析成像(Woodward, 1992)或者波动方程走时层析成像(Luo and Schuster, 1991a; Marquering et al., 1999)，其精度要高于射线路径层析。但是，由于其必须对每个源求波动方程的数值解，所以其成本至少比射线路径层析成像高一个数量级。
走时层析成像还在其他领域中使用，例如医学成像（Duric et al., 2007），非破坏性评估（Leonard et al., 2002）和气象学（Ostashev et al., 2009）。
在本章中，我们介绍射线路径走时层析成像方法。射线理论用于计算源和接收器之间的射线路径，并沿着这些射线路径涂抹加权走时残差以更新模型。层析速度模型是下一次迭代的起始模型，并且由于模型变化会扰动射线，所以射线通常需要重新计算。

5.1 扰动走时积分

3D Eikonal方程（如下所示，推导见公式5.43）：

$$|\nabla t(\mathbf{x})|^2=\left|\frac{\partial t(\mathbf{x})}{\partial x}\right|^2+\left|\frac{\partial t(\mathbf{x})}{\partial y}\right|^2+\left|\frac{\partial t(\mathbf{x})}{\partial z}\right|^2=s(\mathbf{x})^2 \tag{5.1}$$

可以用来推导扰动慢度介质中的走时积分。这里，$t(\mathbf x)$是一个event从源传播到特定检波器的走时，$s(\mathbf x)$是慢度分布。

走时积分方程是走时层析成像中有效更新慢度模型的关键方程。它的推导始于定义背景慢度场$s(\mathbf x)$的慢度扰动$\delta s(\mathbf x)$。相应的扰动走时场定义为$t(\mathbf x) + \delta t(\mathbf x)$，其中$t(\mathbf x)$是背景模型中未扰动的走时场，$\delta t(\mathbf x)$是由慢度扰动引起的走时扰动。扰动走时场满足Eikonal方程

$$\begin{aligned} |\nabla t(\mathbf{x})+\nabla \delta t(\mathbf{x})|^2 & =|\nabla t(\mathbf{x})|^2+2 \nabla \delta t(\mathbf{x}) \cdot \nabla t(\mathbf{x})+|\nabla \delta t(\mathbf{x})|^2, \\ & =[s(\mathbf{x})]^2+2 \delta s(\mathbf{x}) s(\mathbf{x})+[\delta s(\mathbf{x})]^2 . \end{aligned} \tag{5.2}$$

从公式5.2中减去未扰动的Eikonal方程5.1可得

$$2 \nabla t(\mathbf{x}) \cdot \nabla \delta t(\mathbf{x})+|\nabla \delta t(\mathbf{x})|^2=2 \delta s(\mathbf{x}) s(\mathbf{x})+[\delta s(\mathbf{x})]^2 . \tag{5.3}$$

忽略扰动参数的二阶项可得

$$\begin{aligned} \nabla t(\mathbf{x}) \cdot \nabla \delta t(\mathbf{x}) & =|\nabla t(\mathbf{x})| \hat{d \mathbf{l}} \cdot \nabla \delta t(\mathbf{x}) \\ & =\delta s(\mathbf{x}) s(\mathbf{x}) \end{aligned} \tag{5.4}$$

其中，定义$\hat{d \mathbf{l}}$为平行于微扰动射线方向的单位矢量，所以$\nabla t(\mathbf{x})= |\nabla t(\mathbf{x})| \hat{d \mathbf{l}}$。

定义沿着$\hat{d \mathbf{l}}$的方向导数$d/dl= $\hat{d \mathbf{l}} \cdot \nabla$（此方向导数由背景慢度模型决定，而不是扰动模型），对公式5.4除以$|\nabla t(\mathbf x)|=s(\mathbf x)$得：

$$d \delta t(\mathbf{x}) / d l=\delta s(\mathbf{x}) \tag{5.5}$$

重新整理可得走时扰动相对于慢度变化的Frechet导数

$$d \delta t(\mathbf{x}) / \delta s=d l \tag{5.6}$$

Frechet导数是评价数据对于慢度模型中扰动的敏感性的工具。
公式5.5两侧同时乘以$d l$并沿着旧的从源到检波器$\mathbf x$的路径积分可得到扰动走时积分

$$\delta t(\mathbf{x})=\int_{\text {raypath }} \delta s\left(\mathbf{x}^{\prime}\right) d l^{\prime}, \tag{5.7}$$

这个方程在扰动参数上是一阶精确的。被术语“射线路径”表示的积分路径是源位置和接收器位置$\mathbf x$的一个隐函数。
方程5.7表明，由于慢度扰动而引起的走时扰动是通过沿着旧的射线路径进行加权积分来计算的。这是一种节约计算成本的有效方法，因为走时扰动的计算使用旧的射线路径，而不需要通过扰动慢度模型重新跟踪射线。

5.2 射线路径走时成像

一个超定的线性方程组中，方程数量多于未知数的个数，例如一个$4 \times 2$方程组$\mathbf {Ls} = \mathbf {t}$

$$\left[\begin{array}{cc} 0 & 1 \\ 0 & 1 \\ 1 & 1 \\ 0.7 & 1 \end{array}\right] \cdot\left(\begin{array}{l} s_1 \\ s_2 \end{array}\right)=\left(\begin{array}{c} 0.5 \\ 0.4 \\ 1 \\ 1 \end{array}\right) . \tag{5.8}$$

这些方程是不相容（inconsstent）的，因为没有解能够同时满足所有方程。例如，第一个和第二个方程是冲突的。
与方程5.8相关的物理示例是断层图像，成像实验如图5.1所示。每条射线的旅行时由方程5.7所确定，速度模型被离散为N个未知常慢度的单元。走时积分可以近似通过对第$j$个网格中的第$i$条射线路径段$l_{ig}$加权求和得到

$$t_i=\sum_{j=1}^N l_{i j} s_j \tag{5.9}$$

这里我们假设数据集中共拾取了$M$个旅行时。由此形成了$M\times N$的方程组，用矩阵向量可以表示为$\mathbf {Ls} = \mathbf t$，其中$\mathbf t$表示测量得到的大小为$M \times 1$的旅行时矢量，$\mathbf s$表示大小为$N \times 1$的每个网格中的未知慢度矢量。大小为$M \times N$的矩阵$\mathbf L$包含了先的分段长度。
我们的目标是求解方程组5.9并在每个单元格中找到未知的慢度值$s_j$，即如图1.1中的速度断层图所示。几何上，将方程5.8绘制为在图5.2中所示的虚线，没有交叉点意味着这些方程是不相容的。实际上，这种不相容性可能是由于走时拾取错误和/或用于模拟数据的物理机制是不完备的而产生的。解不存在违反了Hadamard（1902）第1章讨论了适定解的准则。

尽管方程5.8没有确切解，但是在线交叉最密集的地方可以找到近似解。这种折中方案即使目标函数$\epsilon$最小化的最小二乘解：

$$\begin{aligned} \epsilon & =1 / 2[\mathbf{L s}-\mathbf{t}]^T[\mathbf{L s}-\mathbf{t}] \\ & =1 / 2 \sum_{i=1}^4 \sum_{j=1}^2\left(l_{i j} s_j-t_i\right)^2 \\ & =1 / 2 \sum_{i=1}^4 r_i^2 \end{aligned} \tag{5.10}$$

其中$r_i=\sum_j l_{i j} s_j-t_i$是第$i$次的残差，即预测走时和观测走时第$i$个分量的差。如果射线弯曲，那么公式5.8中矩阵的元素将取决于未知数$\mathbf s$。在这种情况下，走时方程是非线性的，必须在每次迭代后为新的慢度模型重新计算射线。
根据s1和s2绘制如图5.2损失函数碗（error bowl）。此损失函数碗的底部直接位于最优解$\mathbf s^*$之上，即为最小二乘解。由于这个损失函数碗有底，所以最小二乘解存在。
通过检查找到最优解

图5.2可能适用于方程组只有几个未知数的情况，但当未知数个数变多时是不切实际的。更系统的方法是认识到在误差碗的底部，损失函数对于所有下标$k$均为零的偏导数均为零，即$\partial \epsilon /\partial s_k=0$，其中

$$\begin{aligned} \partial \epsilon / \partial s_k & =\sum_{i=1}^4 \sum_{j=1}^2\left[\left(l_{i j} s_j-t_i\right) l_{i j} \frac{\partial s_j}{\partial s_k}+\left(l_{i j} s_j-t_i\right) s_j \frac{\partial l_{i j}}{\partial s_k}\right] \\ & \approx \sum_{i=1}^4 \sum_{j=1}^2\left(l_{i j} s_j-t_i\right) l_{i j} \delta_{j k}, \end{aligned} \tag{5.11}$$

其中当$j=k$时，Kronecker delta函数$\delta_{jk}=1$，否则其为0。通常要忽略公式5.11中的最右侧项$s_j \partial l_{i j} / \partial s_k\left(\sum_j l_{i j} s_j-t_i\right)$，一是因为其计算成本太高，二是因为当背景慢度模型足够接近实际模型时，偏导数$\partial l_{ij}/\partial s_k$非常小。使用Gauss-Newton（Gill et al., 1987）求解公式5.11所示的超定方程组。如果$\left|s_j \partial l_{i j} / \partial s_k\right|$和$\left|l_{i j} \partial s_j / \partial s_k\right|$具有相同阶数，然后保留第二项即为高斯-牛顿大残差方法（Gill等人，1981）。参见练习13和14以及第26章。
所以Gauss-Newton法的梯度（Gill等人，1981）如下所示。

$$\begin{aligned} \partial \epsilon / \partial s_k & =\sum_{i=1}^4\left(l_{i k} s_k-t_i\right) l_{i k} \\ & =\left[\mathbf{L}^T(\mathbf{L s}-\mathbf{t})\right]_k=0 \quad k=1,2, \end{aligned} \tag{5.12}$$

这里，第$k$个分量的梯度矢量为$\left[\mathbf{L}^T(\mathbf{L s}-\mathbf{t})\right]_k$。

5.2.1 正规方程

5.2.2 病态方程和正则化

5.3 迭代最速下降解