0 一组测量的均值
假设,在一个实验中对一模型参数
m
1
m_1
m1在有噪声的情况下测量了N次,每次都采集到一个独立数据
d
i
d_i
di。那么,如果随机变量是高斯型的,那么它们的联合概率密度函数能够以方差
σ
2
\sigma^2
σ2和均值
μ
\mu
μ的形式来描述,即
p
(
d
)
=
σ
−
N
(
2
π
)
−
N
/
2
exp
[
−
1
2
σ
−
2
∑
i
=
1
N
[
d
i
−
μ
]
2
]
p(\mathbf{d})=\sigma^{-N}(2 \pi)^{-N / 2} \exp \left[-\frac{1}{2} \sigma^{-2} \sum_{i=1}^{N}\left[d_{i}-{\mu}\right]^{2}\right]
p(d)=σ−N(2π)−N/2exp[−21σ−2i=1∑N[di−μ]2]
数据
d
o
b
s
d^{obs}
dobs可以用N维空间中的一个点来表示,这个空间的坐标轴维
d
1
,
d
2
,
.
.
.
,
d
N
d_1,d_2,...,d_N
d1,d2,...,dN(如下图)。
同样地,数据的概率密度函数也可以通过在这一空间中表示出来(如下图)。
由于假设所有数据具有相同的均值和方差,概率密度函数是球对称的,且直线
d
1
=
d
2
=
.
.
.
=
d
N
d_1=d_2=...=d_N
d1=d2=...=dN穿过该“球形云”。
那么,假设我们先给出一组猜测的均值和方差,从而固定了概率密度函数的中心和半径。之后,我们利用上述的公式可以计算数据的概率密度函数值 p ( d o b s ) p(d^{obs}) p(dobs)。如果猜测的均值和方差的值接近正确值,那么 p ( d o b s ) p(d^{obs}) p(dobs)应该是一个相对较大的值。
我们可以想象,沿着直线滑动上图中的概率云,并调节它的直径,直至它在数据点
d
o
b
s
d^{obs}
dobs处的概率被最大化。这一步骤形象的描述了最大似然法,它估计了分布的未知参数。该方法断言,参数的最佳值最大化了观测数据在实际中能够被观测到的概率。换句话说,概率密度函数在点
d
o
b
s
d^{obs}
dobs的值要尽可能的大。可通过对
p
(
d
o
b
s
)
p(d^{obs})
p(dobs)求导得到,即
∂
p
∂
μ
=
0
\frac{\partial{p}}{\partial{\mu}}=0
∂μ∂p=0
∂
p
∂
σ
=
0
\frac{\partial{p}}{\partial{\sigma}}=0
∂σ∂p=0
为了方便计算,一般使用
log
(
p
)
\log(p)
log(p)来计算,因为它是
p
p
p的单调函数,最大化
log
p
(
d
o
b
s
)
\log{p(d^{obs})}
logp(dobs)与最大化
p
(
d
o
b
s
)
p(d^{obs})
p(dobs)结果是相同的。因此,我们计算似然函数
L
=
log
p
(
d
o
b
s
)
L=\log{p(d^{obs})}
L=logp(dobs)的导数(如下图)。
忽略系数
(
2
π
)
−
N
/
2
(2\pi)^{-N/2}
(2π)−N/2,我们有
L
=
log
(
p
(
d
obs
)
)
=
−
N
log
(
σ
)
−
1
2
σ
−
2
∑
i
=
1
N
(
d
i
obs
−
m
1
)
2
∂
L
∂
m
1
=
0
=
1
2
σ
−
2
2
m
1
∑
i
=
1
N
(
d
i
obs
−
m
1
)
∂
L
∂
σ
=
0
=
−
N
σ
+
σ
−
3
∑
i
=
1
N
(
d
i
obs
−
m
1
)
2
\begin{aligned} L &=\log \left(p\left(\mathbf{d}^{\text {obs }}\right)\right)=-N \log (\sigma)-\frac{1}{2} \sigma^{-2} \sum_{i=1}^{N}\left(d_{i}^{\text {obs }}-m_{1}\right)^{2} \\ \frac{\partial L}{\partial m_{1}} &=0=\frac{1}{2} \sigma^{-2} 2 m_{1} \sum_{i=1}^{N}\left(d_{i}^{\text {obs }}-m_{1}\right) \\ \frac{\partial L}{\partial \sigma} &=0=-\frac{N}{\sigma}+\sigma^{-3} \sum_{i=1}^{N}\left(d_{i}^{\text {obs }}-m_{1}\right)^{2} \end{aligned}
L∂m1∂L∂σ∂L=log(p(dobs ))=−Nlog(σ)−21σ−2i=1∑N(diobs −m1)2=0=21σ−22m1i=1∑N(diobs −m1)=0=−σN+σ−3i=1∑N(diobs −m1)2
求解方程,可获得均值和方差的估计值(所要的解)
μ
e
s
t
=
1
N
∑
i
=
1
N
d
i
o
b
s
{\mu}^{\mathrm{est}}=\frac{1}{N} \sum_{i=1}^{N} d_{i}^{\mathrm{obs}}
μest=N1i=1∑Ndiobs
σ
e
s
t
=
[
1
N
∑
i
=
1
N
(
d
i
o
b
s
−
μ
e
s
t
)
2
]
1
/
2
\sigma^{\mathrm{est}}=\left[\frac{1}{N} \sum_{i=1}^{N}\left(d_{i}^{\mathrm{obs}}-{\mu}^{\mathrm{est}}\right)^{2}\right]^{1 / 2}
σest=[N1i=1∑N(diobs−μest)2]1/2
上述两式中
μ
\mu
μ的估计就是采样均值的一般公式,对
σ
\sigma
σ的估计是均方根误差。我们注意到,这些估计是数据拥有高斯分布这个假设的直接产物。
1 应用于反问题中的最大似然性
1.1 最简单情况
在反问题
G
m
=
d
\mathbf{Gm=d}
Gm=d,其中
d
\mathbf{d}
d为一系列的观测数据,那么按照概率来看,
p
(
d
)
p(\mathbf{d})
p(d)是一个多变量高斯概率密度函数
P
(
d
i
)
=
1
2
π
σ
exp
[
−
(
d
i
−
d
)
2
2
σ
2
]
P\left({d_i}\right)=\frac{1}{\sqrt{2 \pi} \sigma} \exp \left[-\frac{\left( {d_i} - d\right)^{2}}{2 \sigma^{2}}\right]
P(di)=2πσ1exp[−2σ2(di−d)2]
其中
d
i
d_i
di为测量数据,也是概率密度的平均值,
σ
\sigma
σ为相应的标准差。
由于n个互不相关的数据的联合概率分布等于n个独立的高斯分布概率密度的乘积,则数据向量
d
\mathbf{d}
d的联合概率分布为:
P
(
d
)
=
∏
i
=
1
n
(
2
π
)
−
1
/
2
σ
−
1
exp
[
−
(
d
i
−
d
)
2
2
σ
2
]
=
(
2
π
)
−
n
/
2
σ
−
n
exp
[
−
1
2
σ
2
∑
i
=
1
n
(
d
i
−
d
)
2
]
\begin{aligned} P(\mathbf{d}) &=\prod_{i=1}^{n}(2 \pi)^{-1 / 2} \sigma^{-1} \exp \left[-\frac{(d_i -{d})^{2}}{2 \sigma^{2}}\right] \\ &=(2 \pi)^{-n / 2} \sigma^{-n} \exp \left[- \frac{1}{2 \sigma^{2}} \sum_{i=1}^{n}\left(d_{i}-{d}\right)^{2}\right] \end{aligned}
P(d)=i=1∏n(2π)−1/2σ−1exp[−2σ2(di−d)2]=(2π)−n/2σ−nexp[−2σ21i=1∑n(di−d)2]
可简化为:
p
(
d
)
∝
exp
[
−
1
2
(
d
−
G
m
)
T
(
d
−
G
m
)
]
p(\mathbf{d}) \propto \exp \left[-\frac{1}{2}(\mathbf{d}-\mathbf{G m})^{\mathrm{T}}(\mathbf{d}-\mathbf{G m})\right]
p(d)∝exp[−21(d−Gm)T(d−Gm)]
若数据之间具有相关性,则
p
(
d
)
∝
exp
[
−
1
2
(
d
−
G
m
)
T
[
cov
d
]
−
1
(
d
−
G
m
)
]
p(\mathbf{d}) \propto \exp \left[-\frac{1}{2}(\mathbf{d}-\mathbf{G m})^{\mathrm{T}}[\operatorname{cov} \mathbf{d}]^{-1}(\mathbf{d}-\mathbf{G m})\right]
p(d)∝exp[−21(d−Gm)T[covd]−1(d−Gm)]
应用于地球物理的反问题中,数据的多次观测可以得到数据的平均值(
d
o
b
s
d^{obs}
dobs)和方差(
σ
\sigma
σ),而模型参数(
m
m
m)未知。模型参数的最佳值有最大化
p
(
d
o
b
s
)
p(\mathbf{d^{obs}})
p(dobs)决定,相对于最小化
log
p
(
d
o
b
s
)
\log{p(\mathbf{d^{obs}})}
logp(dobs),可简化为下式的最小化问题:
(
d
obs
−
G
m
)
T
[
cov
d
]
−
1
(
d
obs
−
G
m
)
\left(\mathbf{d}^{\text {obs }}-\mathbf{G m}\right)^{\mathrm{T}}[\operatorname{cov} \mathbf{d}]^{-1}\left(\mathbf{d}^{\text {obs }}-\mathbf{G m}\right)
(dobs −Gm)T[covd]−1(dobs −Gm)
从上式可以看出模型参数的最大似然估计相当于加权最小二乘解,其中权矩阵式数据协方差矩阵的逆。
其中权矩阵可记为:
W
e
=
[
cov
d
]
−
1
\mathbf{W}_{e}=[\operatorname{cov} \mathbf{d}]^{-1}
We=[covd]−1
如果,数据是不相关的,且数据集中各数据恰好具用相同的方差(赶巧了!!),那么
[
c
o
v
d
]
=
σ
d
2
I
[cov\mathbf{d}]={\sigma_d}^2\mathbf{I}
[covd]=σd2I,那么最大似然解对应简单的最小二乘解。
如果,数据不相关,但具有不同的方差,那么可以利用方差归一化预测误差(见附录2):
E
=
∑
i
=
1
N
σ
d
i
−
2
e
i
2
E=\sum_{i=1}^{N} \sigma_{d i}^{-2} e_{i}^{2}
E=i=1∑Nσdi−2ei2
其最大似然解对应加权最小二乘解。
1.2 先验分布
当线性问题欠定时,最小二乘解不存在。以概率论的观点看待这个问题就是,这时数据的概率密度函数
p
(
d
o
b
s
)
p(\mathbf{d^{obs}})
p(dobs)并不拥有关于模型参数变化的明确极大值。即便是最好的情况,它也仅拥有一个概率最大的脊状结构(如下图所示)。
上图A,有明显的尖峰,可以通过最小二乘法求解。图B对应欠定问题,仅存在概率最大的脊状结构。
为了求解欠定问题,我们必须增加先验信息以使分布拥有明确的尖峰。为了这一目标,首先将关于模型参数的先验信息写为概率密度函数 p A ( m ) p_A(\mathbf{m}) pA(m),其中下标 A A A代表先验的。那么这个概率密度函数的均值使模型向量的期望值。
模型参数的先验分布可以取各种形式。例如,如果我们希望模型参数接近 < m > <\mathbf{m}> <m>,可以使用具有均值 < m > <\mathbf{m}> <m>和方差的高斯分布,如下图。
上图用概率密度函数
p
(
m
1
,
m
2
)
p(m_1,m_2)
p(m1,m2)代表模型参数
m
1
m_1
m1和
m
2
m_2
m2的先验信息。概率密度函数的宽度反映了模型参数的确定度:图A较为确定;图B相对不确定。
如果一个模型参数的先验信息比其他参数值更加确定,那么我们可以对不同模型参数使用不同的方差,如下图
上图中可以看到,模型参数在
<
m
>
<\mathbf{m}>
<m>附近,对
m
1
m_1
m1的不确定度要小于
m
2
m_2
m2的不确定度。
一般高斯情况是:
p
A
(
m
)
∝
exp
[
−
1
2
(
m
−
⟨
m
⟩
)
T
[
cov
m
]
A
−
1
(
m
−
⟨
m
⟩
)
]
p_{\mathrm{A}}(\mathbf{m}) \propto \exp \left[-\frac{1}{2}(\mathbf{m}-\langle\mathbf{m}\rangle)^{\mathrm{T}}[\operatorname{cov} \mathbf{m}]_{\mathrm{A}}^{-1}(\mathbf{m}-\langle\mathbf{m}\rangle)\right]
pA(m)∝exp[−21(m−⟨m⟩)T[covm]A−1(m−⟨m⟩)]
等式约束可以用一个含有脊状结构的分布来实现,如下图所示
图A表示
m
1
m_1
m1和
m
2
m_2
m2的值未知,但认为
m
1
m_1
m1和
m
2
m_2
m2是相关的;图B使用特定方差的高斯概率密度函数对图A的近似。可以看到,这个分布是非高斯型的,然而如果模型参数的期望区间较小,这个分布可以用一个具有非零方差的高斯分布来近似。
对于不等式约束,也可以用一个先验分布来代表,但它本质上是非高斯型的
上图中,模型参数
m
1
m_1
m1和
m
2
m_2
m2的值未知,但认为有
m
1
≤
m
2
m_1{\leq}m_2
m1≤m2这样的关系。
同理,也可以使用一个先验概率密度函数
p
A
(
d
)
p_A(\mathbf{d})
pA(d)来描述关于测量的信息。它恰当地描述了观测,因此它的均值为
d
o
b
s
\mathbf{d^{obs}}
dobs,且它的方差是数据的先验协方差
[
c
o
v
d
]
[cov\mathbf{d}]
[covd]:
p
A
(
d
)
∝
exp
[
−
1
2
(
d
−
d
o
b
s
)
T
[
cov
d
]
−
1
(
d
−
d
o
b
s
)
]
p_{\mathrm{A}}(\mathbf{d}) \propto \exp \left[-\frac{1}{2}\left(\mathbf{d}-\mathbf{d}^{\mathrm{obs}}\right)^{\mathrm{T}}[\operatorname{cov} \mathbf{d}]^{-1}\left(\mathbf{d}-\mathbf{d}^{\mathrm{obs}}\right)\right]
pA(d)∝exp[−21(d−dobs)T[covd]−1(d−dobs)]
p
A
(
d
)
p_{\mathrm{A}}(\mathbf{d})
pA(d)和
p
A
(
m
)
p_{\mathrm{A}}(\mathbf{m})
pA(m),分别描述了数据
d
\mathbf{d}
d和模型参数
m
\mathbf{m}
m的信息
在求解反问题之前,我们可以通过以下两个步骤对关于反问题的了解情况进行总结:首先对数据定义一个先验概率密度函数
p
A
(
d
)
p_{\mathrm{A}}(\mathbf{d})
pA(d),然后将其与模型的先验密度函数
p
A
(
m
)
p_{\mathrm{A}}(\mathbf{m})
pA(m)。先验数据概率密度函数总结了观测,其均值为
d
o
b
s
\mathbf{d^{obs}}
dobs,且它的方差见附录2。由于先验模型概率密度函数完全独立于数据的实际值,我们可以将两者通过简单地相乘来构建其联合概率密度函数:
p
A
(
m
,
d
)
=
p
A
(
m
)
p
A
(
d
)
p_{\mathrm{A}}(\mathbf{m}, \mathbf{d})=p_{\mathrm{A}}(\mathbf{m}) p_{\mathrm{A}}(\mathbf{d})
pA(m,d)=pA(m)pA(d)
这个概率密度函数可以被形象地描述为以观测数据和先验模型参数为中心的概率云(cloud of probability),它的宽度反映了这些量的确定程度,如下图所示,其分布函数的尖峰位于均值
m
a
p
m^{ap}
map和
d
o
b
s
d^{obs}
dobs处。
2 最大似然性解
2.1 准确理论
假设模型与数据的关系式为
g
(
m
)
=
d
\mathbf{g(m)=d}
g(m)=d(线性/非线性方程),该方程在模型参数和数据所张成的空间中呈一个曲面,而方程的解就位于这个曲面上。那么,最大似然问题就转化为在曲面
d
=
g
(
m
)
\mathbf{d=g(m)}
d=g(m)上,寻找联合概率分布
p
A
(
m
,
d
)
p_A(\mathbf{m,d})
pA(m,d)(或
log
p
A
(
m
,
d
)
\log{p_A(\mathbf{m,d})}
logpA(m,d))的最大值,可抽象描述为:
maximize
log
[
p
(
m
,
d
)
]
\log [p(\mathbf{m}, \mathbf{d})]
log[p(m,d)]
with the constrant
g
(
m
)
−
d
=
0
\mathbf{g}(\mathbf{m})-\mathbf{d}=\mathbf{0}
g(m)−d=0
如下图所示,图A中模型参数
m
m
m和数据
d
d
d的先验联合概率密度
p
A
(
m
,
d
)
p_A(\mathbf{m,d})
pA(m,d)表达的思想是模型参数接近它的先验值
m
a
p
\mathbf{m^{ap}}
map,且数据接近它的观测值
d
o
b
s
\mathbf{d^{obs}}
dobs(白色圆圈)。在这里,认为模型参数和数据是通过准确的理论公式
d
=
g
(
m
)
\mathbf{d=g(m)}
d=g(m)相关联的(白色曲线)。估计模型参数
m
e
s
t
\mathbf{m^{est}}
mest和预测数据
d
p
r
e
\mathbf{d^{pre}}
dpre落在曲线的最大概率点处(黑点)。图B展示了沿曲线计算的概率密度
p
p
p。
注意,上图中模型参数的先验概率密度函数比观测数据的概率密度函数更加确定(
σ
m
<
σ
d
\sigma_m<\sigma_d
σm<σd),那么模型参数的估计(最大似然点)趋近于先验模型参数,如下图解
m
e
s
t
\mathbf{m^{est}}
mest更加接近于
m
a
p
\mathbf{m^{ap}}
map,但可能会远离
d
o
b
s
d^{obs}
dobs。
反之,若
σ
d
<
σ
m
\sigma_d<\sigma_m
σd<σm,模型参数的估计主要体现了数据所包含的信息,如下图,解接近于
d
o
b
s
d^{obs}
dobs,但有可能会远离
m
a
p
m^{ap}
map
准确理论的最大似然法与长度方法与有很好的对应关系,关于
m
\mathbf{m}
m最小化
Φ
(
m
)
=
L
(
m
)
+
E
(
m
)
\Phi(\mathbf{m})=L(\mathbf{m})+E(\mathbf{m})
Φ(m)=L(m)+E(m),其中:
L
(
m
)
=
(
m
−
⟨
m
⟩
)
T
[
cov
m
]
A
−
1
(
m
−
⟨
m
⟩
)
E
(
m
)
=
(
G
m
−
d
o
b
s
)
T
[
cov
d
]
−
1
(
G
m
−
d
o
b
s
)
\begin{array}{l} L(\mathbf{m})=(\mathbf{m}-\langle\mathbf{m}\rangle)^{\mathrm{T}}[\operatorname{cov} \mathbf{m}]_{\mathrm{A}}^{-1}(\mathbf{m}-\langle\mathbf{m}\rangle) \\ E(\mathbf{m})=\left(\mathbf{G m}-\mathbf{d}^{\mathrm{obs}}\right)^{\mathrm{T}}[\operatorname{cov} \mathbf{d}]^{-1}\left(\mathbf{G m}-\mathbf{d}^{\mathrm{obs}}\right) \end{array}
L(m)=(m−⟨m⟩)T[covm]A−1(m−⟨m⟩)E(m)=(Gm−dobs)T[covd]−1(Gm−dobs)
对比长度方法的加权阻尼最小二乘问题,它具有
ε
2
W
m
=
[
cov
m
]
A
−
1
and
W
e
=
[
cov
d
]
−
1
\varepsilon^{2} \mathbf{W}_{m}=[\operatorname{cov} \mathbf{m}]_{\mathrm{A}}^{-1} \quad \text { and } \quad \mathbf{W}_{e}=[\operatorname{cov} \mathbf{d}]^{-1}
ε2Wm=[covm]A−1 and We=[covd]−1
因此,它的解是方程组
F
m
=
f
\mathbf{Fm=f}
Fm=f的最小二乘解,即
F
T
F
m
=
F
T
f
\mathbf{F^TFm=F^Tf}
FTFm=FTf,其中:
F
=
[
[
cov
d
]
−
1
/
2
G
[
cov
m
]
A
−
1
/
2
I
]
and
f
=
[
[
cov
d
]
−
1
/
2
d
o
b
s
[
cov
m
]
A
−
1
/
2
⟨
m
⟩
]
\mathbf{F}=\left[\begin{array}{c} {[\operatorname{cov} \mathbf{d}]^{-1 / 2} \mathbf{G}} \\ {[\operatorname{cov} \mathbf{m}]_{\mathrm{A}}^{-1 / 2} \mathbf{I}} \end{array}\right] \quad \text { and } \quad \mathbf{f}=\left[\begin{array}{c} {[\operatorname{cov} \mathbf{d}]^{-1 / 2} \mathbf{d}^{\mathrm{obs}}} \\ {[\operatorname{cov} \mathbf{m}]_{\mathrm{A}}^{-1 / 2}\langle\mathbf{m}\rangle} \end{array}\right]
F=[[covd]−1/2G[covm]A−1/2I] and f=[[covd]−1/2dobs[covm]A−1/2⟨m⟩]
矩阵
[
cov
d
]
−
1
/
2
{[\operatorname{cov} \mathbf{d}]^{-1 / 2}}
[covd]−1/2和
[
cov
m
]
−
1
/
2
{[\operatorname{cov} \mathbf{m}]^{-1 / 2} }
[covm]−1/2可以看作是
d
o
b
s
\mathbf{d^{obs}}
dobs和
<
m
>
<\mathbf{m}>
<m>的确定度,它们的数值这其不确定较小时较大。因此方程
F
m
=
f
\mathbf{Fm=f}
Fm=f的上半部分是由数据确定度加权的数据方程
G
m
=
d
e
s
t
\mathbf{Gm=d^{est}}
Gm=dest,下部分是由模型参数确定度加权的先验方程
m
=
<
m
>
\mathbf{m=<m>}
m=<m>。加权最小二乘问题的矩阵
W
m
W_m
Wm和
W
e
W_e
We拥有明确的概率论含义。
2.2 不准确理论
在许多实际问题中,存在与理论有关的误差。有些理论的假设是不显示的,或者现存的理论计算方法仅是复杂的准确理论的近似形式。这这种情况种,方程
g
(
m
)
=
d
\mathbf{g(m)=d}
g(m)=d不能再用一个简单曲面来代表。它已经变得"模糊不清(fuzzy)",因为存在与之相关的误差。除了曲面外,可以想象一个以
g
(
m
)
=
d
\mathbf{g(m)=d}
g(m)=d为中心的概率密度函数
p
g
(
m
,
d
)
p_g(\mathbf{m,d})
pg(m,d),其宽度正比于理论的不确定度。这时,我们不应该在曲面上寻找
p
A
(
m
,
d
)
p_A(\mathbf{m,d})
pA(m,d)的最大似然点,而是将
p
A
(
m
,
d
)
p_A(\mathbf{m,d})
pA(m,d)和
p
g
(
m
,
d
)
p_g(\mathbf{m,d})
pg(m,d)合并为一个单独的分布,并在全部体积内寻找最大似然点。下图中图A代表
p
A
(
m
,
d
)
p_A(\mathbf{m,d})
pA(m,d)的分布。图B表示不准去理论由以准确理论(白线)为中心的条件概率密度
p
g
(
m
,
d
)
p_g(\mathbf{m,d})
pg(m,d)。图C为先验信息和不准确理论的合并。
其合并的形式如下:
p
T
(
m
,
d
)
=
p
A
(
m
,
d
)
p
g
(
m
,
d
)
p_{\mathrm{T}}(\mathbf{m}, \mathbf{d})=p_{\mathrm{A}}(\mathbf{m}, \mathbf{d}) p_{\mathrm{g}}(\mathbf{m}, \mathbf{d})
pT(m,d)=pA(m,d)pg(m,d)
3 以相对熵为指导原则求解反问题
3.1 相对熵、信息增益
从信息论的观点来看,模型参数
m
\mathbf{m}
m和数据
d
\mathbf{d}
d的概率密度函数包含了两者的信息,可以通过信息增益来量化这些信息量,定义为:
S
[
p
A
(
m
)
]
=
∫
p
A
(
m
)
log
[
p
A
(
m
)
p
N
(
m
)
]
d
M
m
S
[
p
A
(
d
)
]
=
∫
p
A
(
d
)
log
[
p
A
(
d
)
p
N
(
d
)
]
d
N
d
\begin{array}{l} S\left[p_{\mathrm{A}}(\mathbf{m})\right]=\int p_{\mathrm{A}}(\mathbf{m}) \log \left[\frac{p_{\mathrm{A}}(\mathbf{m})}{p_{\mathrm{N}}(\mathbf{m})}\right] \mathrm{d}^{M} \mathbf{m} \\ S\left[p_{\mathrm{A}}(\mathbf{d})\right]=\int p_{\mathrm{A}}(\mathbf{d}) \log \left[\frac{p_{\mathrm{A}}(\mathbf{d})}{p_{\mathrm{N}}(\mathbf{d})}\right] \mathrm{d}^{N} \mathbf{d} \end{array}
S[pA(m)]=∫pA(m)log[pN(m)pA(m)]dMmS[pA(d)]=∫pA(d)log[pN(d)pA(d)]dNd
式中,零概率密度函数
p
N
(
m
)
p_N(\mathbf{m})
pN(m)和
p
N
(
d
)
p_N(\mathbf{d})
pN(d)是模型参数和数据完全未知状态。
信息增益可以用作求解反问题的指导原则。其思想是寻找某个概率密度函数
p
T
(
m
)
p_{T}(\mathbf{m})
pT(m), 使
p
T
(
m
)
p_{T}(\mathbf{m})
pT(m)相对于先验概率密度函数
p
A
(
m
)
p_{A}(\mathbf{m})
pA(m)的信息增益
S
\mathbf{S}
S最小化。在
S
\mathbf{S}
S最小化过程中需要添加约束,否则解将变成
p
T
(
m
)
=
p
A
(
m
)
p_{T}(\mathbf{m}) =p_{A}(\mathbf{m})
pT(m)=pA(m)。其中一个必不可少的约束是
p
T
(
m
)
p_{T}(\mathbf{m})
pT(m)积分的面积为1。其他约束的选择依赖于特定问题的具体情况。
3.2 例子-求解欠定问题
minimize
:
S
=
∫
p
T
(
m
)
log
(
p
T
(
m
)
p
A
(
m
)
)
d
M
m
\operatorname{minimize}: S=\int p_{\mathrm{T}}(\mathbf{m}) \log \left(\frac{p_{\mathrm{T}}(\mathbf{m})}{p_{\mathrm{A}}(\mathbf{m})}\right) \mathrm{d}^{M} m
minimize:S=∫pT(m)log(pA(m)pT(m))dMm
使用约束:
∫
p
T
(
m
)
d
M
m
=
1
and
∫
p
T
(
m
)
(
d
−
G
m
)
d
M
m
=
0
\int p_{\mathrm{T}}(\mathbf{m}) \mathrm{d}^{M} m=1 \quad \text { and } \quad \int p_{\mathrm{T}}(\mathbf{m})(\mathbf{d}-\mathbf{G m}) \mathrm{d}^{M} m=0
∫pT(m)dMm=1 and ∫pT(m)(d−Gm)dMm=0
式中,最终分布
p
T
(
m
)
p_{\mathrm{T}}(\mathbf{m})
pT(m)是未知的,且先验分布
p
A
(
m
)
p_{\mathrm{A}}(\mathbf{m})
pA(m)是已经指定好的(已知的)。注意,第二个约束表明,误差
e
=
d
−
G
m
\mathbf{e=d-Gm}
e=d−Gm的平均(期望)值是0.
利用拉格朗日乘数法,可得:
Φ
(
m
)
=
p
T
log
(
p
T
)
−
p
T
log
(
p
A
)
+
λ
0
p
T
+
λ
T
(
d
−
G
m
)
p
T
\Phi(\mathbf{m})=p_{\mathrm{T}} \log \left(p_{\mathrm{T}}\right)-p_{\mathrm{T}} \log \left(p_{\mathrm{A}}\right)+\lambda_{0} p_{\mathrm{T}}+\boldsymbol{\lambda}^{\mathrm{T}}(\mathbf{d}-\mathbf{G m}) p_{\mathrm{T}}
Φ(m)=pTlog(pT)−pTlog(pA)+λ0pT+λT(d−Gm)pT
∂
Φ
∂
p
T
=
0
=
log
(
p
T
)
+
1
−
log
(
p
A
)
+
λ
0
+
λ
T
(
d
−
G
m
)
\frac{\partial \Phi}{\partial p_{\mathrm{T}}}=0=\log \left(p_{\mathrm{T}}\right)+1-\log \left(p_{\mathrm{A}}\right)+\lambda_{0}+\lambda^{\mathrm{T}}(\mathbf{d}-\mathrm{Gm})
∂pT∂Φ=0=log(pT)+1−log(pA)+λ0+λT(d−Gm)
或
p
T
(
m
)
=
p
A
(
m
)
exp
{
−
(
1
+
λ
0
)
−
λ
T
(
d
−
G
m
)
}
p_{\mathrm{T}}(\mathbf{m})=p_{\mathrm{A}}(\mathbf{m}) \exp \left\{-\left(1+\lambda_{0}\right)-\mathbf{\lambda}^{\mathrm{T}}(\mathbf{d}-\mathbf{G m})\right\}
pT(m)=pA(m)exp{−(1+λ0)−λT(d−Gm)}
又因为
p
A
(
m
)
∝
exp
{
−
1
2
(
m
−
⟨
m
⟩
)
T
[
cov
m
]
A
−
1
(
m
−
⟨
m
⟩
)
}
p_{\mathrm{A}}(\mathbf{m}) \propto \exp \left\{-\frac{1}{2}(\mathbf{m}-\langle\mathbf{m}\rangle)^{\mathrm{T}}[\operatorname{cov} \mathbf{m}]_{A}^{-1}(\mathbf{m}-\langle\mathbf{m}\rangle)\right\}
pA(m)∝exp{−21(m−⟨m⟩)T[covm]A−1(m−⟨m⟩)}
那么,
p
T
(
m
)
∝
exp
{
−
A
(
m
)
}
p_{\mathrm{T}}(\mathbf{m}) \propto \exp \{-A(\mathbf{m})\}
pT(m)∝exp{−A(m)}
其中
A
(
m
)
=
1
2
(
m
−
⟨
m
⟩
)
T
[
cov
m
]
A
−
1
(
m
−
⟨
m
⟩
)
−
(
1
+
λ
0
)
−
λ
T
(
d
−
G
m
)
A(\mathbf{m})=\frac{1}{2}(\mathbf{m}-\langle\mathbf{m}\rangle)^{\mathrm{T}}[\operatorname{cov} \mathbf{m}]_{\mathrm{A}}^{-1}(\mathbf{m}-\langle\mathbf{m}\rangle)-\left(1+\lambda_{0}\right)-\mathbf{\lambda}^{\mathrm{T}}(\mathbf{d}-\mathbf{G m})
A(m)=21(m−⟨m⟩)T[covm]A−1(m−⟨m⟩)−(1+λ0)−λT(d−Gm)
要寻找的模型参数最佳估计
m
e
s
t
\mathbf{m^{est}}
mest是这个分布的均值,也是这个分布的最大似然点。这个点出现在
A
(
m
)
\mathbf{A(m)}
A(m)取极小值的位置:
∂
A
∂
m
q
=
0
or
0
=
[
cov
m
]
A
−
1
(
m
e
s
t
−
⟨
m
⟩
)
+
G
T
λ
\frac{\partial A}{\partial m_{q}}=0 \quad \text { or } \quad 0=[\operatorname{cov} \mathbf{m}]_{\mathrm{A}}^{-1}\left(\mathbf{m}^{\mathrm{est}}-\langle\mathbf{m}\rangle\right)+\mathbf{G}^{\mathrm{T}} \mathbf{\lambda}
∂mq∂A=0 or 0=[covm]A−1(mest−⟨m⟩)+GTλ
左乘
G
[
c
o
v
m
]
A
\mathbf{G}[cov\mathbf{m}]_A
G[covm]A,回代如
d
=
G
m
\mathbf{d=Gm}
d=Gm,可得
λ
=
{
G
[
cov
m
]
A
G
T
}
−
1
{
d
−
G
⟨
m
⟩
}
\boldsymbol{\lambda}=\left\{\mathbf{G}[\operatorname{cov} \mathbf{m}]_{\mathrm{A}} \mathbf{G}^{\mathrm{T}}\right\}^{-1}\{\mathbf{d}-\mathbf{G}\langle\mathbf{m}\rangle\}
λ={G[covm]AGT}−1{d−G⟨m⟩}
最终得到解:
m
e
s
t
−
⟨
m
⟩
=
[
cov
m
]
A
G
T
{
G
[
cov
m
]
A
G
T
}
−
1
{
d
−
G
⟨
m
⟩
}
\mathbf{m}^{\mathrm{est}}-\langle\mathbf{m}\rangle=[\operatorname{cov} \mathbf{m}]_{\mathrm{A}} \mathbf{G}^{\mathrm{T}}\left\{\mathbf{G}[\operatorname{cov} \mathbf{m}]_{\mathrm{A}} \mathbf{G}^{\mathrm{T}}\right\}^{-1}\{\mathbf{d}-\mathbf{G}\langle\mathbf{m}\rangle\}
mest−⟨m⟩=[covm]AGT{G[covm]AGT}−1{d−G⟨m⟩}
4 三种观点的等效
对于同一线性反问题的一般解,三种观点是等效的
观点1:解是通过最小化
L
2
L_2
L2预测误差和
L
2
L_2
L2解的简单程度的线性组合来获得的,即最小化:
e
T
W
e
e
+
ε
2
[
m
−
⟨
m
⟩
]
T
W
m
[
m
−
⟨
m
⟩
]
\mathbf{e}^{\mathrm{T}} \mathrm{W}_{e} \mathbf{e}+\varepsilon^{2}[\mathbf{m}-\langle\mathbf{m}\rangle]^{\mathrm{T}} \mathrm{W}_{m}[\mathbf{m}-\langle\mathbf{m}\rangle]
eTWee+ε2[m−⟨m⟩]TWm[m−⟨m⟩]
观点2:解是通过最小化模型分辨率,数据分辨率,和模型方差大小三项的线性组合获得的,最小化:
α
1
spread
(
R
)
+
α
2
spread
(
N
)
+
α
3
size
(
[
cov
u
m
]
)
\alpha_{1} \operatorname{spread}(\mathbf{R})+\alpha_{2} \operatorname{spread}(\mathbf{N})+\alpha_{3} \operatorname{size}\left(\left[\operatorname{cov}_{u} \mathbf{m}\right]\right)
α1spread(R)+α2spread(N)+α3size([covum])
观点3:解是通过最大化数据,先验模型参数和理论的联合高斯分布的似然性来获得的,最大化:
L
=
log
p
T
(
m
,
d
)
L=\log p_{\mathrm{T}}(\mathbf{m}, \mathbf{d})
L=logpT(m,d)
附录
1 方差与协方差
把一个有误差的数据
d
i
d_i
di当作随机变量处理时,它的方差定义为:
Var
(
d
i
)
=
lim
L
→
∞
1
L
∑
i
=
1
L
(
d
i
l
−
d
i
^
)
2
\operatorname{Var}\left({d_i}\right)=\lim _{L \rightarrow \infty} \frac{1}{L} \sum_{i=1}^{L}\left(d_{il}-\hat{d_{i}}\right)^{2}
Var(di)=L→∞limL1i=1∑L(dil−di^)2
式中,
L
L
L为测量次数,
d
i
^
\hat{d_{i}}
di^为随机变量
d
i
d_i
di的数学期望
E
(
d
i
)
E(d_i)
E(di),也就是精确的数据,而
d
i
l
d_{il}
dil为第
l
l
l次的测量值。两个带误差的数据
d
i
d_i
di和
d
j
d_j
dj在统计上相关时它们的互协方差定义为以下乘积的数学期望
Cov
(
d
i
,
d
j
)
=
E
{
(
d
i
−
d
i
^
)
(
d
j
−
d
j
^
)
}
=
lim
L
→
∞
1
L
∑
i
=
1
L
(
d
i
l
−
d
i
^
)
(
b
j
l
−
d
j
^
)
\begin{aligned} \operatorname{Cov}\left(d_i, {d}_{j}\right) &=E\left\{\left({d}_{i}-\hat{d_{i}}\right)\left({d}_{j}-\hat{d_{j}}\right)\right\} \\ &=\lim _{L \rightarrow \infty} \frac{1}{L} \sum_{i=1}^{L}\left({d}_{i l}-\hat{d_{i}}\right)\left(b_{jl}-\hat{d_{j}}\right) \end{aligned}
Cov(di,dj)=E{(di−di^)(dj−dj^)}=L→∞limL1i=1∑L(dil−di^)(bjl−dj^)
对于整个数据集(向量
d
\mathbf{d}
d),其协方差矩阵为
C
d
=
[
Var
(
d
1
)
Cov
(
d
1
,
d
2
)
⋯
Cov
(
d
1
,
d
n
)
Cov
(
d
2
,
d
1
)
Var
(
d
2
)
⋯
Cov
(
d
2
,
d
^
n
)
⋮
⋮
⋮
Cov
(
d
n
,
d
1
)
Cov
(
d
n
,
d
2
)
⋯
Var
(
d
n
)
]
\boldsymbol{C}_{d}=\left[\begin{array}{cccc} \operatorname{Var}\left(d_1\right) & \operatorname{Cov}\left({d}_{1}, {d}_{2}\right) & \cdots & \operatorname{Cov}\left({d}_{1}, {d}_{n}\right) \\ \operatorname{Cov}\left({d}_{2}, {d}_{1}\right) & \operatorname{Var}\left({d}_{2}\right) & \cdots & \operatorname{Cov}\left({d}_{2}, \hat{d}_{n}\right) \\ \vdots & \vdots & & \vdots \\ \operatorname{Cov}\left({d}_{n}, {d}_{1}\right) & \operatorname{Cov}\left({d}_{n}, {d}_{2}\right) & \cdots & \operatorname{Var}\left({d}_{n}\right) \end{array}\right]
Cd=⎣⎢⎢⎢⎢⎡Var(d1)Cov(d2,d1)⋮Cov(dn,d1)Cov(d1,d2)Var(d2)⋮Cov(dn,d2)⋯⋯⋯Cov(d1,dn)Cov(d2,d^n)⋮Var(dn)⎦⎥⎥⎥⎥⎤
值得注意的是,在进行实际地球物理数据的反演时,要求出上式的数据协方差,这要求多次重复测量,增加了数据采集成本。在数据统计上不相关而且是由多台不同仪器观测取得时,可以利用不同仪器在同一基点处重复观测的数据确定每台仪器的方差,并把此方差当作某台仪器在某天观测的数据方差,用来对各仪器的数据做归一化处理。
2 归一化,为什么?
参数化后的地球物理模型可以表述为下面的形式:
m
T
=
(
ρ
1
,
⋯
,
ρ
4
,
k
1
,
⋯
⋅
k
4
,
μ
1
,
⋯
,
μ
4
,
r
1
,
r
2
,
r
3
)
m^{T}=\left(\rho_{1}, \cdots, \rho_{4}, k_{1}, \cdots \cdot k_{4}, \mu_{1}, \cdots, \mu_{4}, r_{1}, r_{2}, r_{3}\right)
mT=(ρ1,⋯,ρ4,k1,⋯⋅k4,μ1,⋯,μ4,r1,r2,r3)
其中多个物理量具有不同量纲。由于量纲的原因,它们的数值可以相差几个数量级,此时jacobian矩阵的列向量的模也差别极大,以致性状很坏。因此,在进行数值反演计算之前,应该把所有模型参数 转化为无量纲的量,即归一化。对于数据也有类似的情况,就数据而言,不同的数据具有不同的方差,因此可以用方差或其他统计参数是数据和模型归一化。
以数据为例,假设各个数据在统计上是独立的,标准差为
σ
i
,
i
=
1
,
2
,
⋯
,
N
\sigma_i,i=1,2,\cdots,N
σi,i=1,2,⋯,N(
N
N
N为数据的个数),则
d
i
/
σ
i
d_i/\sigma_i
di/σi是一个无量纲的量。即对角阵
E
=
[
1
/
σ
1
⋱
⋱
⋱
1
/
σ
n
]
\boldsymbol{E}=\left[\begin{array}{cccc} 1 / \sigma_{1} & & & \\ & \ddots & & \\ & & \ddots & \\ & & & \ddots \\ & & & & 1 / \sigma_{n} \end{array}\right]
E=⎣⎢⎢⎢⎢⎡1/σ1⋱⋱⋱1/σn⎦⎥⎥⎥⎥⎤
那么归一化的数据为
d
′
=
E
d
=
E
G
m
\mathbf{d^{\prime}}=\mathbf{Ed}=\mathbf{EGm}
d′=Ed=EGm
或写为:
G
′
m
=
d
′
,
G
′
=
E
G
,
d
′
=
E
d
\boldsymbol{G}^{\prime} \boldsymbol{m}=\boldsymbol{d}^{\prime}, \boldsymbol{G}^{\prime}=\boldsymbol{E} \boldsymbol{G}, \boldsymbol{d}^{\prime}=\boldsymbol{E} \boldsymbol{d}
G′m=d′,G′=EG,d′=Ed
若用最小二乘法求解方程组,则此时需要把预测误差
r
r
r与对角阵
E
\mathbf{E}
E乘积的
L
2
L_2
L2范数最小化:
∥
E
r
∥
2
=
r
T
E
T
E
r
=
r
T
W
e
−
1
r
\| \boldsymbol{E} r \|_{2}=\boldsymbol{r}^{\boldsymbol{T}} \boldsymbol{E}^{\boldsymbol{T}} \boldsymbol{E} \boldsymbol{r}=\boldsymbol{r}^{\boldsymbol{T}} \boldsymbol{W}_{e}^{-1} \boldsymbol{r}
∥Er∥2=rTETEr=rTWe−1r