title: 地震反演第三章最速下降法
date: 2022-12-21 14:23:46
description: 最速下降，节选自《Seismic Inversion》
copyright: true
mathjax: true
tags: 最速下降法

categories: Seismic Inversion

第3章最速下降法

本章介绍最速下降优化方法（Steppest descent，SD），其通过沿着最速下降，即负梯度$-\rm \textbf g$方向来迭代地求解优化问题。对于大型优化问题，由于其不需要计算Hessian矩阵的逆，所以其相对于牛顿方法而言计算量更小，但是由于不考虑曲率信息，其Hessian是病态的，所以可能会导致收敛效率较低。然而经验告诉我们，基于预处理和正则化的迭代最速下降法，再加上多尺度反演方法，对于解决大型的地震数据反演问题是非常有效的。

3.1 最速下降法介绍

在多数三维地球物理问题中，求取和存储Hessian矩阵$\rm \textbf H$的逆是不现实的。例如一个共$1000^3=10^9$网格点的速度模型，其Hessian矩阵的逆的求取需要对$10^{18}$个元素进行存储和计算。另一种避免计算大型矩阵的逆，但是可能会导致低收敛效率的方法是最速下降法。回顾公式2.18：

$\rm \textbf x^{(k+1)}=\rm \textbf x^{(k)}-\alpha[\rm \textbf H^{(k)}]^{-1}\rm \textbf g^{(k)} \tag{2.18}$

用单位矩阵代替$\rm \textbf H^{-1}$，即可得到最速下降的公式

$\rm \textbf x^{(k+1)}=\rm \textbf x^{(k)}-\alpha\rm \textbf g^{(k)} \tag{3.1}$

其中，$\alpha$是一个正数，其表示步长。$\rm \textbf x^{(k+1)}$的迭代公式沿着负梯度方向（下坡方向）$-\rm \textbf g^{(k)}$寻找一个更好的解。给定起始点$\rm \textbf x^{(k)}$以及负梯度方向$-\rm \textbf g^{(k)}$，需要进一步计算$\alpha$以使得$f(\rm \textbf x^{(k)}-\alpha\rm \textbf g^{(k)})$最小。与之相等价的是，计算的步长需要使$\rm \textbf x^{(k+1)} - \rm \textbf x^{(k)}$与 $\rm \textbf x^{(k+1)}$的等值线相切或接触（如图3.1.b所示）。当残差矢量（误差函数的均方根）小于某个设定的值后迭代终止。

3.1.1收敛速度

3.2 步长计算

精确线搜索和数值线搜索是两种常用的计算公式3.1中的步长$\alpha$来使以下公式成立的算法：

$\underset{\alpha>0} {\rm min} f(\rm \textbf x^{(k)}-\alpha \Delta \rm \textbf x)\tag{3.7}$

其中，$-\Delta \rm \textbf x$是预定义的搜索方向，即当前模型的下山方向。精确线搜索方法使用解析公式确定$\alpha$。而数值线搜索方法则在方向$-\Delta \rm \textbf x$上通过$\alpha$的多项式来近似$f(\rm \textbf x)$（Nocedal and Wright，1999。通过对梯度和曲率施加限定条件来避免陷入局部极小。
限制性步长（restricted-step）方法只在一个置信区域（trust region）内搜索最优模型$\rm \textbf x^{(k)}+\alpha \Delta \rm \textbf x$，这个置信区域以$\rm \textbf x$为中心，以$r$为半径，同时置信区域和步长之间的关系满足$||\alpha \Delta \rm \textbf x||<r$。由于根据非线性目标函数估计得到的梯度和曲率只在$\rm \textbf x$附近是合适的，而在较远的变量位置是不合适的，因此这种限制有时是必须的。
与直接使用原目标函数$f(\rm \textbf x)=\rm \textbf x^T\rm \textbf H \rm \textbf x+\rm \textbf g^T\rm \textbf x$相比，也可以通过构造一个阻尼目标函数$\tilde f(\rm \textbf x)=\rm \textbf x^T\rm \textbf H \rm \textbf x+\eta^2\rm \textbf x^T\rm \textbf x+\rm \textbf g^T\rm \textbf x$来使病态Hessian稳定下来。这里,$\eta^2$是一个控制置信区域的正整数。$\eta^2$越大，步长越小，所以目标函数$\Delta \tilde f=\tilde f(\rm \textbf x+\Delta \tilde{\rm \textbf x})-\tilde f(\rm \textbf x)$的减少量将会比$\Delta f= f(\rm \textbf x+\Delta \rm \textbf x)- f(\rm \textbf x)$更小。如果$\Delta \tilde f/ \Delta f \geq 0.5$，阻尼不够大，所以需要增大$\eta^2$直到$\Delta \tilde f/ \Delta f$在0.25到0.5之间，对于大型3D问题，试错法（trail-and-error）的计算效率是非常低的。### 3.2.1 精确线搜索
下面给出精确线搜索的推导。一维线搜索问题定义为，对于给定的$\rm \textbf x$和$\Delta \rm \textbf x$，找到最优的$\alpha$来使$f(\rm \textbf x+\alpha\Delta \rm \textbf x)$最小，这种情况下，公式2.9变为

$f({\rm \textbf x}+\alpha\Delta {\rm \textbf x}) \approx f({\rm \textbf x}) +\alpha {\rm \textbf g}^T \Delta {\rm \textbf x} +\frac{\alpha^2}{2}\Delta {\rm \textbf x}^T \rm \textbf H \Delta \rm \textbf x. \tag{3.8}$

上式对$\alpha$求导，并令导数等于0，可得：

$\frac{\partial f(\rm \textbf x+\alpha \Delta \rm \textbf x)}{\partial \alpha} =\rm \textbf g^T\Delta \rm \textbf x+\alpha \Delta \rm \textbf x^T \rm \textbf H \Delta \rm \textbf x =0\tag{3.9}$

可得到$\alpha$的解析解为：

$\alpha ^\ast = \frac{-\rm \textbf g^T \Delta \rm \textbf x}{\Delta \rm \textbf x^T \rm \textbf H \Delta \rm \textbf x}\tag{3.10}$

由此求得的步长为精确步长。这种线搜索方法为精确线搜索，因为如果$f(\rm \textbf x)$是二次泛函，那么公式3.10中的步长是完全正确的。对于最速下降法，$\Delta \rm \textbf x=-\rm \textbf g$，所以：

$\alpha^{\ast}= \frac{\rm \textbf g^T \rm \textbf g}{ \rm \textbf g^T \rm \textbf H \rm \textbf g} \tag{3.10}$

经验告诉我们，精确线搜索方法应用效果非常不好，而地震反问题中通常使用数值线搜索方法。

3.2.2 非精确牛顿法和非精确线搜索

近似Hessian矩阵逆的方法被称为非精确牛顿法（Gill et al., 1981）。如果有一小部分反演是基于共轭梯度法来迭代近似求解的，那么非精确牛顿法又被称为截断牛顿法（truncated Newton method）。在一定的迭代次数后，共轭梯度法的解将会位于最速下降法和牛顿矢量之间，如图3.2所示。AlTheyab et al., (2013)提出了一种内外双循环的FWI算法，在内部循环中通过线性共轭梯度法来迭代近似Hessian的逆。在章节21中介绍了这种方法，其内部循环迭代得到的速度模型被用作外部循环的背景模型，而外部循环是通过非线性共轭梯度法进行迭代的。
通过设定终止条件来停止线搜索的方法被称为非精确线搜索方法。常用的停止准则有Wolfe准则（Nocedal and Wright, 1999）：

$\overset {actual function} { f(\rm \textbf x^{(k)}) } \leq\overset {linear \space prediction\space of \space function} { f(\rm \textbf x^{(k)} + c_1\alpha\Delta \rm \textbf x^{(k)T} \nabla f(\rm \textbf x^{(k)}) } \tag{3.12}$

其中，$0<c_1<1$，函数在$\rm \textbf x^{(k)} + c_1\alpha \rm \textbf x{(k)}$的预测要大于实际值，否则步长会超过线性近似的有效范围。用户设置一个初始步长值，如果步长满足上述条件，则会计算一个新的梯度并重新搜索；如果步长不满足上述条件，原始步长减半并重新测试是否满足Wolfe准则。
有时较小的步长会很容易满足上述准则，进而导致搜索效率低。为了避免过小的步长，我们也需要通过曲率条件来找到一个更合适的步长：

$\overset {slope \space at \space new \space pt.} { { \Delta \rm \textbf x^{(k)} } \nabla f(\rm \textbf x^{(k)}+\alpha {\rm \textbf x^{(k)})} } \geq \overset {scaled \space slope \space at \space starting \space pt.} { c_2\Delta {\rm \textbf x^{(k)}} \nabla f(\rm \textbf x^{(k)}) } \tag{3.13}$

其中，$0<c_1<c_2<1$。在这一准则中，沿着$\rm \textbf x{(k)}$处的搜索方向的放缩后的斜率要小于或者等于$\rm \textbf x^{(k)}+\alpha\rm \textbf x^{(k)}$处的该值。这样就可以避免步长较小，因为$\rm \textbf x^{(k)}+\alpha\rm \textbf x^{(k)}$处的斜率可能等于起始点$\rm \textbf x{(k)}$处的值。

3.2.3 数值线搜索法

当高度非线性的损失函数会打破局部二次假设，这时确切线搜索方法就不再准确了。取而代之可以在下坡方向的多个点上计算损失函数来估计步长。例如。二次线搜索首先定义搜索方向$\Delta \rm \textbf x$，然后计算多个可能步长，例如$\alpha_1=0.5$和$\alpha_2=1$，然后计算两个目标泛函$f(\rm \textbf x+\alpha_1\Delta \rm \textbf x)$和$f(\rm \textbf x + \alpha_2 \Delta \rm \textbf x)$，看他们是否小于$f(\rm \textbf x)$。然后通过下述公式来对泛函进行插值：

$\begin{align} \notag f(\rm \textbf x+\alpha\Delta \rm \textbf x) & = \frac{(\alpha-\alpha_1) (\alpha-\alpha_2)}{\alpha_1\alpha_2} f(\rm \textbf x) \\ \notag & + \frac{\alpha(\alpha-\alpha_2)}{\alpha_1(\alpha_1-\alpha_2)} f(\rm \textbf x+\Delta \alpha_1 \Delta \rm \textbf x) \\ \notag & + \frac{\alpha(\alpha-\alpha_1)}{\alpha_2(\alpha_2-\alpha_1)} f(\rm \textbf x+\Delta \alpha_2 \Delta \rm \textbf x) 0\leq\alpha \leq \alpha_2. \end{align} \tag{3.14}$

上述函数是一个关于$\alpha$的二次多项式。上述函数最小时$\alpha$的值可以通过求$f(\alpha)$对$\alpha$的导数并令其为零，然后求解$\alpha$。
方程3.14需要求目标函数$f(\rm \textbf x)$在三处的值，所以另一种方法是利用$f(\rm \textbf x), f(\rm \textbf x + \alpha_1 \Delta \rm \textbf x)$以及$f(\rm \textbf x)$在$\rm \textbf x$的方向导数进行二次插值：

$f({\rm \textbf x} + \Delta {\rm \textbf x}) = \frac{\alpha^2}{\alpha_1^2}[f({\rm \textbf x} + \Delta {\rm \textbf x}) - \alpha_1 f ({\rm \textbf x})'-f({\rm \textbf x})] + \alpha f({\rm \textbf x})'+f (\rm \textbf x) \tag{3.15}$

其中，$f(\rm \textbf x)’=\frac{df(\rm \textbf x+\alpha\Delta \rm \textbf x)}{d\alpha}$表示$f(\rm \textbf x)$研制搜索方向的导数。

3.2.4 二维平面最小化

有时我们需要从同一数据中反演两种完全不同的模型参数，这些参数的单位和对数据敏感性可能完全不同。例如在旅行时反演中，可能要同时反演慢度分布和静态延迟（static delays），但静态延迟对模型变化的敏感度要更弱一些。当一个矢量分量是各向同性慢度模型，另一个是各向异性参数时也是一个很好的例子。
为了加快收敛速度，我们可以在一个与静态延迟和慢度模型相关的由矢量$\Delta \rm \textbf x_1$和$\Delta \rm \textbf x_2$构成的二维超平面内搜索。也就是说，搜索方向$\Delta \rm \textbf x$可以分解为：

$\Delta \rm \textbf x=\alpha_1 \Delta \rm \textbf x_1+\alpha_2\Delta \rm \textbf x_2\tag{3.16}$

此时，步长求取问题变为找到值$\alpha_1$和$\alpha_2$使得$f({\rm \textbf x}+\alpha_1 \Delta {\rm \textbf x}_1+\alpha_2 \Delta {\rm \textbf x}_2)$最小。将公式3.16带入2.9可得

$\begin{align} \notag f({\rm \textbf x}+\alpha_1 \Delta {\rm \textbf x}_1+\alpha_2 \Delta {\rm \textbf x}_2 \notag & \approx f({\rm \textbf x})+(\alpha_1 \Delta {\rm \textbf x}_1+\alpha_2 \Delta {\rm \textbf x}_2)^T {\rm \textbf g} \notag & + \frac{1}{2}(\alpha_1 \Delta {\rm \textbf x}_1+\alpha_2 \Delta {\rm \textbf x}_2)^T{\rm \textbf H}(\alpha_1 \Delta {\rm \textbf x}_1+\alpha_2 \Delta {\rm \textbf x}_2). \end{align} \tag{3.17}$

将上式分别对$\alpha_1$和$\alpha_2$求导，并令导数为0，可得

$\begin{bmatrix} {\rm \textbf x}^T_1{\rm \textbf H}\Delta {\rm \textbf x}_1 & {\rm \textbf x}^T_1{\rm \textbf H}\Delta {\rm \textbf x}_2 \\ {\rm \textbf x}^T_2{\rm \textbf H}\Delta {\rm \textbf x}_1 & {\rm \textbf x}^T_2{\rm \textbf H}\Delta {\rm \textbf x}_2 \\ \end{bmatrix} \begin{bmatrix} \alpha_1 \\ \alpha_2 \end{bmatrix} + \begin{bmatrix} \Delta {\rm \textbf x}_1^T {\rm \textbf g}\\ \Delta {\rm \textbf x}_2^T {\rm \textbf g}\\ \end{bmatrix} = \begin{bmatrix} 0\\0 \end{bmatrix} \tag{3.18}$

解以上$2\times 2$矩阵即可求得$\alpha_1$和$\alpha_2$。当不同模型参数的单位不同时，使用二维平面搜索相比一维线搜索能加快收敛效率。这一步骤能够拓展到更多参数、多数据联合的反演中(Kennett et al., 1988; Oldenburg et al., 1993）。
把这一过程看作是预处理最速下降法是很有意思的。这两个步骤（预处理和二维平面搜索）都是想把椭圆目标函数整形为圆形目标函数。如果目标函数对于某个特定分量是拉长的（例如，盆地模型的深部反射系数对地震数据的影响与浅部相比差别不大），所以可以将深部模型设置为$\Delta {\rm \textbf x}_2$，二将浅层模型设置为$\Delta {\rm \textbf x}_1$，这样图3.1中拉长的靶心可以在一步内完成收敛。

3.3 最速下降法和线性方程系统

正定矩阵${\rm \textbf L}$构成的$N \times N$线性方程系统${\rm \textbf {Lx}}={\rm \textbf t}$，其有解且能够最小化二次函数（线性方程组的求解问题可转化为求以下函数最小值问题）

$f({\rm \textbf x})= -{\rm \textbf t}^T{\rm \textbf x}+\frac{1}{2}{\rm \textbf x}^T{\rm \textbf {Lx}} \tag{3.19}$

$f({\rm \textbf x})$的梯度为$({\rm \textbf {Lx}}-{\rm \textbf t})$，$Hessian$为$\rm \textbf L$，所以最速下降公式变为：

$\begin{align} \notag \rm \textbf x^{(k+1)} & =\rm \textbf x^{(k)}-\alpha ({\rm \textbf {Lx}}^{(k)}-{\rm \textbf t}) \\ &= {\rm \textbf x}^{(k)}-\alpha {\rm \textbf r}^{(k)} \end{align} \tag{3.20}$

其中，残差矢量为${\rm \textbf r}^{(k)}= {\rm \textbf {Lx}}^{(k)}-{\rm \textbf t}$。

$M \times N$超定线性方程组${\rm \textbf {Lx}}={\rm \textbf t}(M>N)$也能够利用最速下降法最小化目标函数求解：

$\begin{align} \notag f({\rm \textbf x}) &= \frac{1}{2} ({\rm \textbf {Lx}}-{\rm \textbf t}) ^T ({\rm \textbf {Lx}}-{\rm \textbf t}) \\ &=\frac{1}{2} {\rm \textbf t}^T {\rm \textbf t} -{\rm \textbf t}^T {\rm \textbf {Lx}} + \frac{1}{2}{\rm \textbf x}^T{\rm \textbf L}^T{\rm \textbf {Lx}} \end{align} \tag{3.21}$

对未知数$\rm \textbf x$求导可得$f(\rm \textbf x)$的梯度为

$\nabla f (\rm\textbf x)={\rm\textbf L}^T{\rm\textbf {Lx}}-{\rm\textbf L}^T{\rm\textbf t} \tag{3.22}$

此时最速下降法的解为：

$\begin{align} \notag \rm \textbf x^{(k+1)} & =\rm \textbf x^{(k)}-\alpha {\rm \textbf L}^T ({\rm \textbf {Lx}}^{(k)}-{\rm \textbf t}) \\ &= {\rm \textbf x}^{(k)}-\alpha {\rm \textbf L}^T {\rm \textbf r}^{(k)} \end{align} \tag{3.23}$

其中${\rm \textbf r}^{(k)}= {\rm \textbf {Lx}}^{(k)}-{\rm \textbf t}$为第k次迭代的残差。通过对比3.23和3.20可以发现，与方阵$\rm \textbf L$相比，残差项前多了一个$\rm \textbf L$的转置项。如果$\rm \textbf L$的特征值为$\sigma _i$那么，${\rm \textbf L}^T{\rm \textbf L}$的特征值为$\sigma _i^2$。所以一般方程要比方阵更病态。后面层析和偏移章节中会提到，这个转置矩阵相当于层析中的反传旅行时残差和最小二乘偏移中的偏移地震道。

3.3.1 正则化最速下降

如果$\mathbf H =\mathbf L^T \mathbf L$的条件数非常大，就会出现通过预测多个模型$\mathbf x$都能得到几乎完全相同数据的情况，这样就会违反第1.2章中讨论的Hadamard关于病态问题的描述。为了避免这种情况，我们可以施加约束让解必须很小，即让解在起始位置附近。对于一个线性方程组，这等价于限制性的$(M+N)\times N$方程组：

$\left[\begin{array}{c} \mathbf{L} \\ \eta \mathbf{I} \end{array}\right] \mathbf{x}=\left(\begin{array}{l} \mathbf{t} \\ \mathbf{0} \end{array}\right), \tag{3.24}$

其中，$\mathbf I$是$N\times N$的单位矩阵，$\mathbf 0$是$N \times 1$的零矢量，$\eta$是阻尼参数（正标量），较大的$\eta$使解更加满足此约束，而不是满足$\mathbf {Lx}=\mathbf t$，而较小的$\eta$则使解更满足原方程组。
通过对以上方程组两侧乘以$[\mathbf{L}^T \eta \mathbf I]$可得如下对应方程：

$\left[\mathbf{L}^T \mathbf{L}+\eta^2 \mathbf{I}\right] \mathbf{x}=\mathbf{L}^T \mathbf{t}$

与公式3.22相类似，以上方程的解是可以最小化目标函数$\epsilon$的平稳矢量：

$f(\mathbf{x})=\overbrace{\frac{1}{2}(\mathbf{L} \mathbf{x}-\mathbf{t})^T(\mathbf{L} \mathbf{x}-\mathbf{t})}^{\text {data misfit }}+\overbrace{\frac{\eta^2}{2} \mathbf{x}^T \mathbf{x}}^{\text {penalty term }} \tag{3.26}$

以上公式表明可以将通过向损失函数中添加加权惩罚函数的形式，将正则化引入到一般的非线性优化问题中，正如公式1.9所示的那样。这也是正则化如何实际应用的。随着迭代的继续，当解足够接近$\mathbf x=0$时$\eta$的值逐渐减小到零，这种方法被称为Levenberg-Marquardt，阻尼最小二乘或零阶Tikhonov 法(Aster et al., 2005)，如附录3A所示。附录3B描述了如何选择$\eta$的值。
因此，$f(\mathbf x)$的梯度为

$\nabla f(\mathbf{x})=\mathbf{L}^T(\mathbf{L} \mathbf{x}-\mathbf{t})+\eta^2 \mathbf{x} \tag{3.27}$

正则化最速下降法的迭代公式为：

$\begin{aligned} \mathbf{x}^{(k+1)} & =\mathbf{x}^{(k)}-\alpha\left(\left[\mathbf{L}^T \mathbf{L}+\eta^2 \mathbf{I}\right] \mathbf{x}^{(k)}-\mathbf{L}^T \mathbf{t}\right) \\ & =\mathbf{x}^{(k)}-\alpha\left(\mathbf{L}^T \mathbf{r}^{(k)}+\eta^2 \mathbf{x}^{(k)}\right) \end{aligned}$

在第1章讨论的旅行时层析中，残差矢量$\mathbf{r}^{(k)}=\mathbf{L} \mathbf{x}^{(k)}-\mathbf{t}$是观测旅行时$\mathbf t$和预测旅行时$\mathbf {Lx}^{(k)}$之间的残差。

3.3.3 预处理最速下降

如果$\mathbf P$是Hessian矩阵的逆$\matbf H^{-1}$的低成本近似，那么最速下降法的左-预处理公式为

$\mathbf{x}^{(k+1)}=\mathbf{x}^{(k)}-\alpha \mathbf{P g}^{(k)} \tag{3.29}$

其中$\mathbf {PH}$的条件数应该少于$\mathbf H$。对于地震成像而言，scaled的对角矩阵是较为inexpensive的预处理算子，其中$\mathbf P_{ii}$的元素均为$\mathbf H_{ii}$的倒数。如果$\mathbf H_{ii}$接近0，那么正则化预处理算子为$\mathbf P_{ii}=1/(\mathbf H_{ii}+\eta)$，这里$\eta$为阻尼参数。预处理的作用是重塑椭圆轮廓$f(\mathbf x)$为圆形轮廓，因此我们可以快速迭代到$f(\mathbf x)$的最小值。

3.4 总结

本章给出了求解非线性和线型方程的各种无约束最速下降法的概述。预处理和正则化最速下降法的通用公式如3.29所示。如果$\eta^2=0$并且如果Hessian矩阵的主对角线元素分量都等于值1，其变为标准最速下降公式。如果经过预处理，通常可以实现更快的收敛，尤其是Hessian 矩阵是病态的情况下。病态条件也表明非唯一性问题，可以通过向损失函数中添加正则化项来解决。因此，通过正则化将解引导到模型空间中我们认为更合理的区域，预处理则让我们更快地到达那里。