线性高斯反问题的解--长度方法1

最新推荐文章于 2023-09-07 14:33:21 发布

原创最新推荐文章于 2023-09-07 14:33:21 发布 · 578 阅读

0 ·

CC 4.0 BY-SA版权

地球物理专栏收录该内容

12 篇文章

订阅专栏

本文深入探讨了最小二乘法的基本原理，包括长度概念、范数推广及其在解决线性反问题中的应用。通过具体实例，如直线拟合、抛物线拟合和平面拟合，展示了最小二乘法的计算过程。同时，文章分析了最小二乘法在面对过定、恰定和欠定问题时的表现，并提出了解决欠定问题的策略。

3.1 长度

首先，以一个简单的线性回归问题（一条直线拟合数据的问题），来直观的体会“长度”。
a)对数据（z,d）的最小二乘直线拟合；b) $e_i = d_{i}^{obs} - d_{i}^{pre}$

那么，最佳的拟合直线所具有模型参数（截距和斜率）将使总误差E最小，即
$\sum\limits_{i=1}^{n}e_{i}^{2}=\bf {e^Te} = min$
式中，总误差 $E$ 恰恰使向量 $e\bf e$ 的欧几里得长度。

从长度方法的观点来看，最小二乘法是通过寻找预测误差的最小长度所对应的模型参数（截距和斜率）来估计反问题的解。通常来讲在，求解反问题的过程中使用长度方法是最简单的方法，也符合人的直观感受。

3.2 范数==》长度的推广

所谓范数是指对某些长度或大小的度量，其公式如下：
$∥e∥p=(∑i=1nep)1p{\|e\|}_{p} = (\sum\limits_{i=1}^{n}e^p)^{\frac{1}{p}}$
以下三种范数最为常用：
$L_{1}$ norm: $∥e∥1=[∑i∣ei∣1]\quad\|\mathbf{e}\|_{1}=\left[\sum_{i}\left|e_{i}\right|^{1}\right]$
$L_{2}$ norm: $∥e∥2=[∑i∣ei∣2]1/2\quad\|\mathbf{e}\|_{2}=\left[\sum_{i}\left|e_{i}\right|^{2}\right]^{1 / 2}$
$L∞L_{\infty}$ norm: $∥e∥∞=max⁡i∣ei∣\quad\|\mathbf{e}\|_{\infty}=\max _{i}\left|e_{i}\right|$
幂次越高的范数给 $e\mathbf{e}$ 中最大元素更大的权重， $L∞L_{\infty}$ 表明仅对 $e\mathbf{e}$ 的最大元素施加权重
在这里插入图片描述
从上图上看， $∣ e ∣$ 的显著性最明显，而 $e^{10}|$ 最差

那么，最小二乘法究竟该选择哪种长度（即采取哪种范数），这个问题的答案涉及对远离平均趋势的离群数据进行加权的方式,如下图所示：
在这里插入图片描述
图中， $L_1$ 范数给离群点的权重最小。

3.3 线性反问题的最小二乘解

最小二乘能够以一种非常直接的方式推广到一般线性反问题。
$E=eTe=(d−Gm)T(d−Gm)=∑i=1N[di−∑j=1MGijmj][di−∑k=1MGikmk]E=\mathbf{e}^{\mathrm{T}} \mathbf{e}=(\mathbf{d}-\mathbf{G} \mathbf{m})^{\mathrm{T}}(\mathbf{d}-\mathbf{G} \mathbf{m})=\sum_{i=1}^{N}\left[d_{i}-\sum_{j=1}^{M} G_{i j} m_{j}\right]\left[d_{i}-\sum_{k=1}^{M} G_{i k} m_{k}\right]$
为避免混乱，对上式进行重新整理：
$E=∑j=1M∑k=1Mmjmk∑i=1NGijGik−2∑j=1Mmj∑i=1NGijdi+∑i=1NdidiE=\sum_{j=1}^{M} \sum_{k=1}^{M} m_{j} m_{k} \sum_{i=1}^{N} G_{i j} G_{i k}-2 \sum_{j=1}^{M} m_{j} \sum_{i=1}^{N} G_{i j} d_{i}+\sum_{i=1}^{N} d_{i} d_{i}$

令 $∂E∂mp=0\frac{\partial E}{\partial m_{p}}=0$
第一项为：
$∂∂mq[∑j=1M∑k=1Mmjmk∑i=1NGijGik]=∑j=1M∑k=1M[δjqmk+mjδkq]∑i=1NGijGik=2∑k=1Mmk∑i=1NGiqGik\begin{aligned} \frac{\partial}{\partial m_{q}}\left[\sum_{j=1}^{M} \sum_{k=1}^{M} m_{j} m_{k} \sum_{i=1}^{N} G_{i j} G_{i k}\right] &=\sum_{j=1}^{M} \sum_{k=1}^{M}\left[\delta_{j q} m_{k}+m_{j} \delta_{k q}\right] \sum_{i=1}^{N} G_{i j} G_{i k} \\ &=2 \sum_{k=1}^{M} m_{k} \sum_{i=1}^{N} G_{i q} G_{i k} \end{aligned}$

注意：模型参数是相互独立的变量， $∂mi/∂mj\partial m_{i} / \partial m_{j}$ 形式的导数只有在 $i = j$ 时候等于 $1$ ，而 $\neq j$ 是为0，因此 $∂mi/∂mj\partial m_{i} / \partial m_{j}$ 克罗内克函数（Kronecker delta） $δij\delta_{i j}$ ，包含它的公式可以明显的简化

第二项为：
$\frac{\partial}{\partial m_{q}}\left[\sum_{j=1}^{M} m_{j} \sum_{i=1}^{N} G_{i j} d_{i}\right]=-2 \sum_{j=1}^{M} \delta_{j q} \sum_{i=1}^{N} G_{i j} d_{i}=-2 \sum_{i=1}^{N} G_{i q} d_{i}$

第三项为：
$∂∂mq[∑i=1Ndidi]=0\frac{\partial}{\partial m_{q}}\left[\sum_{i=1}^{N} d_{i} d_{i}\right]=0$

结合三项，可得：
$∂E∂mq=0=2∑k=1Mmk∑i=1NGiqGik−2∑i=1NGiqdi\frac{\partial E}{\partial m_{q}}=0=2 \sum_{k=1}^{M} m_{k} \sum_{i=1}^{N} G_{i q} G_{i k}-2 \sum_{i=1}^{N} G_{i q} d_{i}$

重新写成矩阵形式：
$GTGm−GTd=0\mathbf{G}^{\mathrm{T}} \mathbf{G m}-\mathbf{G}^{\mathrm{T}} \mathbf{d}=0$
假设 $[GTG]−1\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1}$ 存在（有不存在的时候），那么有如下解：
$mest=[GTG]−1GTd\mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1} \mathbf{G}^{\mathrm{T}} \mathbf{d}$

对于维度较大的情况， $GTG\mathbf{G}^{\mathrm{T}} \mathbf{G}$ 的计算成本可能很高，而且 $GTG\mathbf{G}^{\mathrm{T}} \mathbf{G}$ 极少像 $G\mathbf G$ 那样稀疏。在这种情况下，优先考虑迭代的矩阵求解方案，如双共轭梯度算法（biconjugate gradient method）

3.4 一些最小二乘法的例子

1）直线拟合问题：

模型是 $d_i =m_1+m_2z_i$ ,方程 $Gm=d\mathbf{Gm=d}$ 形式如下：
$[1z11z2⋮⋮1zN][m1m2]=[d1d2⋮dN]\left[\begin{array}{cc} 1 & z_{1} \\ 1 & z_{2} \\ \vdots & \vdots \\ 1 & z_{N} \end{array}\right]\left[\begin{array}{c} m_{1} \\ m_{2} \end{array}\right]=\left[\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{N} \end{array}\right]$
==》
$GTG=[11⋯1z1z2⋯zN][1z11z2⋮⋮1zN]=[N∑i=1Nzi∑i=1Nzi∑i=1Nzi2]\mathbf{G}^{\mathrm{T}} \mathbf{G}=\left[\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ z_{1} & z_{2} & \cdots & z_{N} \end{array}\right]\left[\begin{array}{cc} 1 & z_{1} \\ 1 & z_{2} \\ \vdots & \vdots \\ 1 & z_{N} \end{array}\right]=\left[\begin{array}{cc} N & \sum_{i=1}^{N} z_{i} \\ \sum_{i=1}^{N} z_{i} & \sum_{i=1}^{N} z_{i}^{2} \end{array}\right]$

$GTd=[11⋯1z1z2⋯zN][d1d2⋮dN]=[∑i=1Ndi∑i=1Ndizi]\mathbf{G}^{\mathrm{T}} \mathbf{d}=\left[\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ z_{1} & z_{2} & \cdots & z_{N} \end{array}\right]\left[\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{N} \end{array}\right]=\left[\begin{array}{c} \sum_{i=1}^{N} d_{i} \\ \sum_{i=1}^{N} d_{i} z_{i} \end{array}\right]$
$GTd=[11⋯1z1z2⋯zN][d1d2⋮dN]=[∑i=1Ndi∑i=1Ndizi]\mathbf{G}^{\mathrm{T}} \mathbf{d}=\left[\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ z_{1} & z_{2} & \cdots & z_{N} \end{array}\right]\left[\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{N} \end{array}\right]=\left[\begin{array}{c} \sum_{i=1}^{N} d_{i} \\ \sum_{i=1}^{N} d_{i} z_{i} \end{array}\right]$

$mest=[GTG]−1GTd=[N∑i=1Nzi∑i=1Nzi∑i=1Nzi2]−1[∑i=1Ndi∑i=1Ndizi]\mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1} \mathbf{G}^{\mathrm{T}} \mathbf{d}=\left[\begin{array}{cc} N & \sum_{i=1}^{N} z_{i} \\ \sum_{i=1}^{N} z_{i} & \sum_{i=1}^{N} z_{i}^{2} \end{array}\right]^{-1}\left[\begin{array}{c} \sum_{i=1}^{N} d_{i} \\ \sum_{i=1}^{N} d_{i} z_{i} \end{array}\right]$

2 抛物线拟合

模型是 $d_{i}=m_{1}+m_{2} z_{i}+m_{3} z_{i}^{2}$ ，方程 $Gm=d\mathbf{Gm=d}$ 形式如下：
$[1z1z121z2z22⋮⋮⋮1zNzN2][m1m2m3]=[d1d2⋮dN]\left[\begin{array}{ccc} 1 & z_{1} & z_{1}^{2} \\ 1 & z_{2} & z_{2}^{2} \\ \vdots & \vdots & \vdots \\ 1 & z_{N} & z_{N}^{2} \end{array}\right]\left[\begin{array}{c} m_{1} \\ m_{2} \\ m_{3} \end{array}\right]=\left[\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{N} \end{array}\right]$
==》
$GTG=[11⋯1z1z2⋯zNz12zN2⋯zN2][1z1z121z2z22⋮⋮⋮1zNzN2]\mathbf{G}^{\mathrm{T}} \mathbf{G}=\left[\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ z_{1} & z_{2} & \cdots & z_{N} \\ z_{1}^{2} & z_{N}^{2} & \cdots & z_{N}^{2} \end{array}\right]\left[\begin{array}{ccc} 1 & z_{1} & z_{1}^{2} \\ 1 & z_{2} & z_{2}^{2} \\ \vdots & \vdots & \vdots \\ 1 & z_{N} & z_{N}^{2} \end{array}\right]$
$GTd=[11⋯1z1z2⋯zNz12zN2⋯zN2][d1d2⋮dN]=[∑i=1Ndi∑i=1Nzidi∑i=1Nzi2di]\mathbf{G}^{\mathrm{T}} \mathbf{d}=\left[\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ z_{1} & z_{2} & \cdots & z_{N} \\ z_{1}^{2} & z_{N}^{2} & \cdots & z_{N}^{2} \end{array}\right]\left[\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{N} \end{array}\right]=\left[\begin{array}{c} \sum_{i=1}^{N} d_{i} \\ \sum_{i=1}^{N} z_{i} d_{i} \\ \sum_{i=1}^{N} z_{i}^{2} d_{i} \end{array}\right]$

$mest=[GTG]−1GTd=[N∑i=1Nzi∑i=1Nzi2∑i=1Nzi∑i=1Nzi2∑i=1Nzi3∑i=1Nzi2∑i=1Nzi3∑i=1Nzi4]−1[∑i=1Ndi∑i=1Nzidi∑i=1Nzi2di]\mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1} \mathbf{G}^{\mathrm{T}} \mathbf{d}=\left[\begin{array}{ccc} N & \sum_{i=1}^{N} z_{i} & \sum_{i=1}^{N} z_{i}^{2} \\ \sum_{i=1}^{N} z_{i} & \sum_{i=1}^{N} z_{i}^{2} & \sum_{i=1}^{N} z_{i}^{3} \\ \sum_{i=1}^{N} z_{i}^{2} & \sum_{i=1}^{N} z_{i}^{3} & \sum_{i=1}^{N} z_{i}^{4} \end{array}\right]^{-1}\left[\begin{array}{c} \sum_{i=1}^{N} d_{i} \\ \sum_{i=1}^{N} z_{i} d_{i} \\ \sum_{i=1}^{N} z_{i}^{2} d_{i} \end{array}\right]$
在这里插入图片描述
上图为开普勒第三定律的验证，它阐述了行星轨道半径的立方正比于轨道周期的平方。a)对于太阳系，数据（红圈）用二次方公式 $d_i=m_1+m_2z_i+m_3z_{i}^{2}$ 实现了最小二乘拟合，其中 $d_i$ 为半径的立方， $z_i$ 为周期。b)拟合误差，将数据和误差独立显示的目的是将误差以恰当的尺度画出。数据源于维基百科（Wikipedia）。

3、平面拟合

模型： $d_{i}=m_{1}+m_{2} x_{i}+m_{3} y_{i}$ ，方程 $Gm=d\mathbf{Gm=d}$ 形式如下：
$[1x1y11x2y2⋮⋮⋮1xNyN][m1m2m3]=[d1d2⋮dN]\left[\begin{array}{ccc} 1 & x_{1} & y_{1} \\ 1 & x_{2} & y_{2} \\ \vdots & \vdots & \vdots \\ 1 & x_{N} & y_{N} \end{array}\right]\left[\begin{array}{c} m_{1} \\ m_{2} \\ m_{3} \end{array}\right]=\left[\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{N} \end{array}\right]$
==》
$GTG=[11⋯1x1x2⋯xNy1y2⋯yN][1x1y11x2y2⋮⋮⋮1xNyN]\mathbf{G}^{\mathrm{T}} \mathbf{G}=\left[\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ x_{1} & x_{2} & \cdots & x_{N} \\ y_{1} & y_{2} & \cdots & y_{N} \end{array}\right]\left[\begin{array}{ccc} 1 & x_{1} & y_{1} \\ 1 & x_{2} & y_{2} \\ \vdots & \vdots & \vdots \\ 1 & x_{N} & y_{N} \end{array}\right]$
$=[N∑i=1Nxi∑i=1Nyi∑i=1Nxi∑i=1Nxi2∑i=1Nxiyi∑i=1Nyi∑i=1Nxiyi∑i=1Nyi2]=\left[\begin{array}{ccc} N & \sum_{i=1}^{N} x_{i} & \sum_{i=1}^{N} y_{i} \\ \sum_{i=1}^{N} x_{i} & \sum_{i=1}^{N} x_{i}^{2} & \sum_{i=1}^{N} x_{i} y_{i} \\ \sum_{i=1}^{N} y_{i} & \sum_{i=1}^{N} x_{i} y_{i} & \sum_{i=1}^{N} y_{i}^{2} \end{array}\right]$

$GTd=[11⋯1x1x2⋯xNy1y2⋯yN][d1d2⋮dN]=[∑i=1Ndi∑i=1Nxidi∑i=1Nyidi]\mathbf{G}^{\mathrm{T}} \mathbf{d}=\left[\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ x_{1} & x_{2} & \cdots & x_{N} \\ y_{1} & y_{2} & \cdots & y_{N} \end{array}\right]\left[\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{N} \end{array}\right]=\left[\begin{array}{c} \sum_{i=1}^{N} d_{i} \\ \sum_{i=1}^{N} x_{i} d_{i} \\ \sum_{i=1}^{N} y_{i} d_{i} \end{array}\right]$

$mest=[GTG]−1GTd=[N∑i=1Nxi∑i=1Nyi∑i=1Nxi∑i=1Nxi2∑i=1Nxiyi∑i=1Nyi∑i=1Nxiyi∑i=1Nxiyi2]−1[∑i=1Ndi∑i=1Nxidi∑i=1Nyidi]\mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1} \mathbf{G}^{\mathrm{T}} \mathbf{d}=\left[\begin{array}{ccc} N & \sum_{i=1}^{N} x_{i} & \sum_{i=1}^{N} y_{i} \\ \sum_{i=1}^{N} x_{i} & \sum_{i=1}^{N} x_{i}^{2} & \sum_{i=1}^{N} x_{i} y_{i} \\ \sum_{i=1}^{N} y_{i} & \sum_{i=1}^{N} x_{i} y_{i} & \sum_{i=1}^{N} x_{i} y_{i}^{2} \end{array}\right]^{-1}\left[\begin{array}{c} \sum_{i=1}^{N} d_{i} \\ \sum_{i=1}^{N} x_{i} d_{i} \\ \sum_{i=1}^{N} y_{i} d_{i} \end{array}\right]$

在这里插入图片描述
上图是使用平面拟合来检测地质断层上发生地震的例子，（圆圈）西北太平洋千岛群岛俯冲带发生的地震。 $x$ 坐标轴指向北， $y$ 坐标轴指向东。地震散布在一个由最小二乘确定的倾斜平面附近。数据来自美国地质调查局（United States Geological Survey, USGS）。

最小二乘法的是为了求解没有精确解的反问题。该方法是为了寻找最佳近似解，从数学公式来看： $E=eTe\mathbf {E = e^Te}$ ，寻找的是 $∂E∂m=0\mathbf{\frac{\partial E}{\partial m}=0}$ 的解（ $m^{est}$ ）即最小化 $L_2$ 误差长度（此时，或许 $e≠0e\neq0$ ）,而不是在求 $E=0\mathbf{E = 0}$ 的解（即误差 $e = 0$ ）,所以说最小二乘法的解是近似解（最佳的），而不是精确解

3.5 应用最小二乘的一些缺陷（存在性、唯一性问题）

公式 $mest=[GTG]−1GTd\begin{array}{c} \mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1} \mathbf{G}^{\mathrm{T}} \mathbf{d} \end{array}$ 中，隐含地假设仅存在唯一的最佳近似解。如果反问题存在多个解，并且这些解具有相同的最小预测误差，那么最小二乘法将会失效。

反例

在这里插入图片描述
如上图，考虑仅具有一个数据点的直线拟合问题，显然该问题的解是非唯一的，从图上来看，可以画出可以穿过该点的无穷条直线，且各直线都具有零预测误差（e=0）。
用公式 $mest=[GTG]−1GTd\mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1} \mathbf{G}^{\mathrm{T}} \mathbf{d}$ 去求解的话，则有
$\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1}=\left[\begin{array}{cc} N & \sum_{i=1}^{N} z_{i} \\ \sum_{i=1}^{N} z_{i} & \sum_{i=1}^{N} z_{i}^{2} \end{array}\right]^{-1} \rightarrow\left[\begin{array}{cc} 1 & z_{1} \\ z_{1} & z_{1}^{2} \end{array}\right]^{-1}$
由于 $A−1=A∗∣A∣A^{-1}=\frac{A^{*}}{|A|}$
然而
$\mathbf{|G^TG|} = \frac{1}{z_1^2-z_1^2}=\infty$
显然， $[GTG]−1\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1}$ 是奇异的，也就意味着最小二乘法的公式失效了！！！！！

从反例来看，显然最小二乘法不能解决所有反问题，从这个角度上看，等式 $G m = d$ 没有提供了足够信息来唯一的确定模型参数 $m\mathbf m$ ，基于这一点（也可以说，基于最小二乘法的解是否唯一），可以对反问题进行分类。

1）恰定问题

$Gm=d\mathbf{Gm=d}$ 恰好存在足够的信息来唯一的确定模型参数，且具有零预测误差（e=0）。

2 ）超定问题

当 $Gm=d\mathbf{Gm=d}$ 拥有过多的信息来确定一个精确解时，则为超定问题（overdetermined）。这时，最小二乘法可以找到一个最佳近似解。
超定问题拥有的数据通常多于未知量（ $N > M$ ）。但并非是绝对的，有些问题即使当数据少于未知量（ $N < M$ ）,在一定程度上依然是超定的。

3）欠定问题

当 $Gm=d\mathbf{Gm=d}$ 不能提供足够的信息来唯一地确定模型参数是，则为欠定问题（underdetermined）。就如上如的反例，欠定发生在几个解通常具有零预测误差的情况。这时，最小二乘法是失效的。欠定问题通常未知量多于数据（ $M > N$ ）。但也存在一些问题，当 $M < N$ 时，它们在一定程度上也是欠定的，**在很多情况下，数据唯一的确定了某些模型参数，但不能确定其他模型参数。**如下图所示：
在这里插入图片描述
上图为一个声学实验，由于对第二块砖没有测定慢度，这个模型参数完全不受数据的约束。而第一块砖因多次测量慢度，是超定的。
$h\left[\begin{array}{cc} 1 & 0 \\ 1 & 0 \\ \vdots & \vdots \\ 1 & 0 \end{array}\right]\left[\begin{array}{c} s_{1} \\ s_{2} \end{array}\right]=\left[\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{N} \end{array}\right]$
式中： $s_i$ 为第 $i$ 块砖的慢度； $h$ 为砖的宽度； $d_n$ 为第n次测量的旅行时。在用最小二乘法求解时， $[GTG]−1\mathbf{{[G^TG]}^{-1}}$ 是奇异的，此时即使 $M < N$ ,依然是欠定问题。这时常见的仅有部分模型参数是欠定的情况，在实际的情况中，还有很多欠定的形式更加不易察觉。
另外，并非所有的欠定问题都具有零预测误差。
为此，对欠定问题进一步分类：

将具有非零预测误差的欠定问题称为混定问题（mixed-determined problems）;
将具有零预测误差的欠定问题称为纯欠定问题（purely undetermined problems）;

从以上分类来看，单纯的最小二乘法无法应对欠定问题，需要其他的求解方法来应对欠定问题

3.6 纯欠定问题（purely undetermined problems）

纯欠定问题， $Gm=d\mathbf{Gm=d}$ 方程的数据的数量少于未知模型参数的数量，有多个解使预测误差为0( $E=0\mathbf{E}=0$ )。对于这样的问题，我们必须使用某些手段，从无穷多个解中挑选出一个合适的解（ $m^{est}$ ）。为此，我们必须添加一些没有包含在 $Gm=d\mathbf{Gm=d}$ 中的额外信息，即先验信息。

3.6.1 先验信息举例

先验信息可以采取任何形式，但在任何形式下，它仅描述了期望解有什么样的特征，它并非基于实际数据(即与实际数据无关)
例子1：
在拟合通过一个数据点的直线例子中，加入"期望直线穿过远点"的先验信息。此时，这个先验信息提供了足够的信息来找到唯一的解（因为，两点确定一条直线）。
例子2：
期望模型参数拥有给定的符号（如均为正/负）或位于一个给定的范围。
如假设模型参数代表了地球内部不同点的密度。即使没有测量，人们也可以肯定密度是正值。而且，由于地球内部岩石的密度必然处于一个密度范围之内，如 $\sim 10000 kg/m^3$ 。使用这样的先验信息可以极大的缩小解的范围，甚至使之唯一。

对于欠定问题，不得不将先验信息添加到反演问题，从而才能寻找到一个合适的解。但这种方法也存在着一些不确定性。
这些信息来自何处？它们的确定程度有多大？
这些问题没有固定答案，一些情况下，人们也许能找到合适的先验信息，有些时候则不能。如果人们想要论证估计值的唯一性，那么先验信息的有效性至关重要。

3.6.2 最小长度解

最为普遍的先验信息就是，希望反问题的解是简单的，这里所谓的“简单”是通过解的长度在某种度量下来量化的。其中采用欧几里得空间长度最为常见，即解的长度可定义为
$L=mTm=∑mi2\mathbf{L}=\mathbf{m^Tm}=\sum{m_{i}^2}$
我们所期望的解，其解的长度要最小。
虽然，这个度量也许不是对简单程度的实际度量，但它有时是有用的。在一些特定的问题， $L\mathbf{L}$ 可以是一个物理系统的能量（离散化后则代表功率–王家映）
使用先验信息来约束纯欠定反问题可表述如下：
寻找 $mest\mathbf{m}^{est}$ ,以在 $e=Gm−d=0\mathbf{e=Gm-d=0}$ 的约束下，使 $L=mTm\mathbf{L=m^Tm}$ 最小化。
该问题可通过拉格朗日乘数法来求解：
$\Phi(\mathbf{m})=L+\sum_{i=1}^{N} \lambda_{i} e_{i}=\sum_{i=1}^{M} m_{i}^{2}+\sum_{i=1}^{N} \lambda_{i}\left[d_{i}-\sum_{j=1}^{M} G_{i j} m_{j}\right]$
式中， $λi\lambda_i$ 为拉格朗日乘数，求导可得
$\frac{\partial \Phi}{\partial m_{q}}=\sum_{i=1}^{M} 2 \frac{\partial m_{i}}{\partial m_{q}} m_{i}-\sum_{i=1}^{N} \lambda_{i} \sum_{j=1}^{M} G_{i j} \frac{\partial m_{j}}{\partial m_{q}}=2 m_{q}-\sum_{i=1}^{N} \lambda_{i} G_{i q}$
令 $∂Φ∂mq=0\frac{\partial \Phi}{\partial m_{q}}=0$ ，并将其重新写成矩阵形式，得到方程
$2m=GTλ2\mathbf{m}=\mathbf{G^T\lambda}$
将其带入 $d=Gm\mathbf{d}=\mathbf{Gm}$ 中，得到
$d=Gm=G[GTλ/2]\mathbf{d}=\mathbf{Gm}=\mathbf{G[G^T\lambda/2]}$
式中 $GGT\mathbf{GG^T}$ 是一个 $\times N$ 的方阵，若其存在逆矩阵，那么有 $λ=2[GGT]−1d\mathbf{\lambda}=2\mathbf{[GG^T]^{-1}d}$
将其回代到 $2m=GTλ2\mathbf{m}=\mathbf{G^T\lambda}$ 中，可得到
$mest=GT[GGT]−1d\mathbf{m}^{\mathrm{est}}=\mathbf{G}^{\mathrm{T}}\left[\mathbf{G} \mathbf{G}^{\mathrm{T}}\right]^{-1} \mathbf{d}$

注意，上式的一个条件是 $Gm=d\mathbf{Gm=d}$ 是纯欠定的。

3.7 混定问题（mixed-determined problems）

实际中的大多数反问题，既不是完全超定的，也不是完全欠定的。如下例所示。
在这里插入图片描述
在X射线层析成像中，模型离散为多个小方格，有些方格有多条射线穿过（上图A），显然就这一个网格来说X射线不透明度（模型参数），是超定的。也有一些网格，完全没有X射线穿过（上图B），这些方格是完全欠定的。另外，还存在一些方格，因为射线从两个（或多个）方格中穿过（上图C），所以这些方格不能单独求取，这些方格也是欠定的，因为仅能确定它们的平均不透明度。

方法1

如果可能（理想情况下），我们希望将未知的模型参数分为两组：超定的和欠定的。这样将 $Gm=d==>G′m′=d′\mathbf{Gm=d} ==>\mathbf{G^\prime m^\prime=d^\prime}$ ,其中 $m′\mathbf{m^\prime}$ 分为超定的部分 $m′o\mathbf{m^{\prime o}}$ 和欠定部分 $m′u\mathbf{m^{\prime u}}$ :
$\left[\begin{array}{cc} \mathbf{G}^{\prime 0} & 0 \\ 0 & \mathbf{G}^{\prime \mathrm{u}} \end{array}\right]\left[\begin{array}{c} \mathbf{m}^{\prime 0} \\ \mathbf{m}^{\prime \mathrm{u}} \end{array}\right]=\left[\begin{array}{l} \mathbf{d}^{\prime 0} \\ \mathbf{d}^{\prime \mathrm{u}} \end{array}\right]$
这样上半部分求最小二乘解，下半部分求最小长度解。这样我们将找到一个解，而且给反问题添加了尽量少的先验信息。这个分割过程可以通过数据核的奇异值分解来实现（矩阵 $G$ 的奇异值分解就是把 $m′o\mathbf{m^{\prime o}}$ 和 $m′u\mathbf{m^{\prime u}}$ 分开的一种途径）

方法2

在不分割 $m\mathbf{m}$ 的情况下，阻尼最小二乘法是最常用的。这时我们不仅要求预测误差（E）最小化，还要求解的长度（L）最小化，所以目标函数为二者的某种组合：
$\Phi(\mathbf{m})=E+\varepsilon^{2} L=\mathbf{e}^{\mathrm{T}} \mathbf{e}+\varepsilon^{2} \mathbf{m}^{\mathrm{T}} \mathbf{m}$
式中， $ε2\varepsilon^2$ 决定了预测误差和解的长度的相对重要性。若 $ε\varepsilon$ 足够大，则极小化的过程主要针对上式中第二项进行，这样可以明显的削弱解的欠定部分。但第一项权重相对减少，解估计的预测误差加大（不是最小的E）。结果就是，获得的解既不是最小二乘解，也不是最小长度解。若 $ε\varepsilon$ 设定为0，则最小化过程仅针对上式中的第一项，没有加入先验信息来挑选欠定部分的解。
与最小二乘的推导过程相似，可以给出阻尼最小二乘解的形式：
$\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}+\varepsilon^{2} \mathbf{I}\right] \mathbf{m}^{\mathrm{est}}=\mathbf{G}^{\mathrm{T}} \mathbf{d} \quad \text { or } \quad \mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}+\varepsilon^{2} \mathbf{I}\right]^{-1} \mathbf{G}^{\mathrm{T}} \mathbf{d}$
阻尼最小二乘法可以挑选出 $ε\varepsilon$ 的一些折中值，以期在近似最小化解中欠定部分长度的基础上，近似最小化预测误差。然而，并没有直接的方法来求取折中值，它必须通过试错法来确定（实验/经验）。另外，值得注意的是，误差中不仅包含预测误差，还包含了拟合先验信息误差。

有时候，用 $L=mTm\mathbf{L=m^Tm}$ 来衡量解的简单程度并不合理，例如，假设求解海洋中海水密度扰动这一反问题。我们不希望找到一个最接近0的解，而是想要找到一个接近于某一特定数值（如海水的平均密度）情况下的解。那么，对 $L\mathbf{L}$ 自然可以推广为：
$L=(m−<m>)T(m−<m>)\mathbf{L}=\mathbf{(m-<m>)}^{T}\mathbf{(m-<m>)}$
式中， $< m >$ 为模型参数的先验值（在例子中为已知海水的平均值）。

3.8 利用先验信息加权

用来衡量解的简单程度并不只有长度这一度量，平直度（flat）和光滑度（smooth）也是经常用到的两种方式。
一个空间连续函数的平直度可以通过一阶导数（即陡度（stepness）（平直度的反义词））的范数来量化。对于离散化的模型参数，可以用模型空间（物理）上相邻的模型参数的差值作为一阶导数的近似：
那么向量 $m\mathbf{m}$ 的陡度 $I\mathbf{I}$ 是
$\mathbf{I}=\frac{1}{\Delta x}\left[\begin{array}{cccc} -1 & 1 & & \\ & -1 & 1 & \\ & & \ddots & \ddots \\ & & & -1 \end{array}\right]\left[\begin{array}{c} m_{1} \\ m_{2} \\ \vdots \\ m_{M} \end{array}\right]=\mathbf{D}{\mathbf{m}}$
解的光滑度（smoothness）可以通过使用二阶导数量化粗糙度（roughness）(光滑度的反义词)来实现。
于是，解的总陡度/总粗糙度是如下长度：
$L=\mathbf{I}^{\mathrm{T}} \mathbf{I}=[\mathbf{D m}]^{\mathrm{T}}[\mathbf{D m}]=\mathbf{m}^{\mathrm{T}} \mathbf{D}^{\mathrm{T}} \mathbf{D} \mathbf{m}=\mathbf{m}^{\mathrm{T}} \mathbf{W}_{\mathrm{m}} \mathbf{m}$
矩阵 $Wm=DTD\mathbf{W_m=D^TD}$ 可以当作参与计算模型参数向量 $m\mathbf{m}$ 长度的权重因子。
因此，解的简单程度的度量 $L$ 可以推广为：
$L=[\mathbf{m}-\langle\mathbf{m}\rangle]^{\mathrm{T}} \mathbf{W}_{\mathrm{m}}[\mathbf{m}-\langle\mathbf{m}\rangle]$

另外，预测误差的加权度量也是有用的。有一些观测经常比另一些做的更准确，在这种情况下，我们希望较精确观测的预测误差 $e_i$ 在总误差 $E$ 的量化过程中比不精确观测具有更大的权重。为了实现这样的权重，我们定义一个广义的预测误差：
$E=\mathbf{e}^{\mathrm{T}} \mathbf{W}_{\mathrm{e}} \mathbf{e}$
式中，矩阵 $We\mathbf{W_e}$ 定义了各单独测量误差对总预测误差的相对贡献。通常，该矩阵为一对角阵。
例如，共有5个测量数据（N=5）且第三次观测精度是其他几次观测精度的2倍，那么可以使用：
$\mathbf{W}_{\mathrm{e}}=\left[\begin{array}{cccc} 1 & & & \\ & 1 & & \\ & & 2 & \\ & & 1 & \\ & & & 1 \end{array}\right]$

这样，可以通过加权的方式，进一步改进反问题的解：

1）加权最小二乘==》完全超定问题
通过最小化广义观测误差 $E=eTWee\mathbf{E = e^TW_ee}$ 来求解：
$\mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{W}_{\mathrm{e}} \mathbf{G}\right]^{-1} \mathbf{G}^{\mathrm{T}} \mathbf{W}_{\mathrm{e}} \mathbf{d}$
2）加权最小长度==》纯欠定问题
通过最小化广义长度 $L=[m−⟨m⟩]TWm[m−⟨m⟩]TL=[\mathbf{m}-\langle\mathbf{m}\rangle]^{\mathrm{T}} \mathbf{W}_{\mathrm{m}}[\mathbf{m}-\langle\mathbf{m}\rangle]^{\mathrm{T}}$ 来求解：
$\mathbf{m}^{\mathrm{est}}=\langle\mathbf{m}\rangle+\mathbf{W}_{\mathrm{m}}^{-1} \mathbf{G}^{\mathrm{T}}\left[\mathbf{G} \mathbf{W}_{\mathrm{m}}^{-1} \mathbf{G}^{\mathrm{T}}\right]^{-1}[\mathbf{d}-\mathbf{G}\langle\mathbf{m}\rangle]$
3）加权阻尼最小二乘法==》混定问题（稍微欠定）
通过最小化预测误差和解的长度的组合 $Φ(m)=E+ε2L\Phi(\mathbf{m})=E+\varepsilon^{2} L$ 来求解：
$\left[\mathbf{G}^{\mathrm{T}} \mathbf{W}_{\mathrm{e}} \mathbf{G}+\varepsilon^{2} \mathbf{W}_{\mathrm{m}}\right] \mathbf{m}^{\mathrm{est}}=\mathbf{G}^{\mathrm{T}} \mathbf{W}_{\mathrm{e}} \mathbf{d}+\varepsilon^{2} \mathbf{W}_{\mathrm{m}}\langle\mathbf{m}\rangle$
或
$\mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{W}_{\mathrm{e}} \mathbf{G}+\varepsilon^{2} \mathbf{W}_{\mathrm{m}}\right]^{-1}\left[\mathbf{G}^{\mathrm{T}} \mathbf{W}_{\mathrm{e}} \mathbf{d}+\varepsilon^{2} \mathbf{W}_{\mathrm{m}}\langle\mathbf{m}\rangle\right]$