文章目录
笔者按
{~~~~~~~}
正所谓“温故而知新”。最近,我又开启对“集合xxx滤波”的回顾和学习。
{~~~~~~~}
其实关于集合数据同化我是有这样一个心路历程:
{~~~~~~~}
刚开始接触集合同化,便将一切“集合xxx滤波”归为“EnKF”,觉得理所当然也无可厚非;
{~~~~~~~}
后来随着对“EAKF”的研究,知晓了EnKF和EAKF的区别在于“观测是否扰动”以及“随机性和确定性”,但我始终声称“EAKF是一种确定性的EnKF”并乐此不疲;
{~~~~~~~}
直到最近被审稿人指出“EAKF是一种确定性的EnKF”的提法令他困惑,我才终于意识到是时候、有必要好好掰扯掰扯这些“集合xxx滤波”了。
{~~~~~~~}
于是在导师的帮助下,参考Tippett等(2003),在这篇分享中对卡尔曼滤波(Kalman filter)、均方根滤波(square root filter,SRF)、集合卡尔曼滤波(ensemble Kalman filter,EnKF)、集合均方根滤波(ensemble square root filter,EnSRF)、集合变换卡尔曼滤波(ensemble transform Kalman filter,ETKF)和集合调整卡尔曼滤波(ensemble adjustment Kalman filter,EAKF)的公式推导以及它们的关系作以整理。
{~~~~~~~}
这次的分享可能对初学者不太友好,要求读者对于数据同化、集合滤波、线性代数等数理知识有较好的基础。其实每过一段时间温习一遍类似的数学公式推导(如数理方程、数值模式、同化算法等等)对于沉淀知识净化心灵很有效果(哈哈)。后面有时间我会更新一些入门级的分享。这也是我第一次发布学习笔记,有不恰当之处我会及时改正,请读者多多包涵。欢迎感兴趣的同行一起交流!
1 前言
{~~~~~~~} 集合数据同化方法利用状态估计方法以及预测和分析误差协方差矩阵的低秩形式来同化观测。这类方法的一个关键要素是采用适当的统计手段将预测集合转换为分析集合。这种转换可以通过将观测视为随机变量来“随机性”实现(Houtekamer和Mitchell,1998;Burgers等,1998;Evensen,2003),这类方法完善了Evensen(1994)最初提出的集合卡尔曼滤波(EnKF);也可以要求更新后的分析扰动满足卡尔曼滤波分析误差协方差的方程来“确定性”实现。“随机性”更新分析集合的EnKF方法和后续发展的“确定性”更新分析集合的方法都可被视作卡尔曼均方根滤波(Kalman SRF)的范畴(Bierman,1977;Maybeck,1982;Heemink等,2001;Tippett等,2003)。第2部分依次介绍了卡尔曼滤波、SRF框架下的“随机性”EnKF,并提出“确定性”Kalman SRF的概念。第3部分介绍了四种不同的“确定性”Kalman SRF(即EnSRF)变体。第4部分是简单的总结。
2 卡尔曼均方根滤波(Kalman SRF)
2.1 卡尔曼滤波误差协方差的演化
{~~~~~~~}
为方便后续的表达,首先给出一些常用符号的含义。
{~~~~~~~}
n
n
n,状态变量维度;
{~~~~~~~}
p
p
p,观测个数;
{~~~~~~~}
P
k
f
{\bf{P}}_k^f
Pkf,预报误差协方差矩阵,维度为
n
×
n
n \times n
n×n;
{~~~~~~~}
P
k
a
{\bf{P}}_k^a
Pka,分析误差协方差矩阵,维度为
n
×
n
n \times n
n×n;
{~~~~~~~}
M
k
{\bf{M}}_k^{}
Mk,切线性动力模型;
{~~~~~~~}
Q
k
{{\bf{Q}}_k}
Qk,模型误差协方差矩阵,维度为
n
×
n
n \times n
n×n;
{~~~~~~~}
H
k
{\bf{H}}_k^{}
Hk,观测投影算子,维度为
p
×
n
p \times n
p×n;
{~~~~~~~}
R
k
{\bf{R}}_k^{}
Rk,观测误差协方差矩阵,维度为
p
×
p
p \times p
p×p;
{~~~~~~~}
K
k
{\bf{K}}_k^{}
Kk,卡尔曼增益矩阵,维度为
n
×
p
n \times p
n×p;
{~~~~~~~}
卡尔曼滤波的核心是预报误差协方差和分析误差协方差的演化,它们的演化方程以及卡尔曼增益矩阵的定义如下:
P
k
f
=
M
k
P
k
−
1
a
M
k
T
+
Q
k
(
1
)
{\bf{P}}_k^f = {\bf{M}}_k^{}{\bf{P}}_{k - 1}^a{\bf{M}}_k^{\rm{T}} + {{\bf{Q}}_k}~~~~~~~~(1)
Pkf=MkPk−1aMkT+Qk (1)
P
k
a
=
(
I
−
K
k
H
k
)
P
k
f
(
2
)
{\bf{P}}_k^a = \left( {{\bf{I}} - {\bf{K}}_k^{}{\bf{H}}_k^{}} \right){\bf{P}}_k^f~~~~~~~~(2)
Pka=(I−KkHk)Pkf (2)
K
k
≡
P
k
f
H
k
T
(
H
k
P
k
f
H
k
T
+
R
k
)
−
1
(
3
)
{\bf{K}}_k^{} \equiv {\bf{P}}_k^f{\bf{H}}_k^{\rm{T}}{\left( {{\bf{H}}_k^{}{\bf{P}}_k^f{\bf{H}}_k^{\rm{T}} + {\bf{R}}_k^{}} \right)^{ - 1}}~~~~~~~~(3)
Kk≡PkfHkT(HkPkfHkT+Rk)−1 (3)其中,
k
k
k表示当前分析步。可以看出,在以模型误差作为强迫的动力作用下,分析误差的传播由式(1)给出;式(2)表示在最优的数据同化方案中同化观测得到的分析误差协方差总是小于预报误差协方差。
2.2 矩阵均方根及其与扰动矩阵的关系
{~~~~~~~} 为了在卡尔曼SRF的框架下分析EnKF和EnSRF,有必要清楚矩阵均方根是什么及其与扰动矩阵(扰动向量集合)的关系。对于一个维度为 n × n n \times n n×n的协方差矩阵 (对称正定矩阵),如果它的秩为 m m m,那么总有一个维度为 n × m n \times m n×m的矩阵 P = Z Z T {\bf{P}} = {\bf{Z}}{{\bf{Z}}^{\rm{T}}} P=ZZT这里允许 m ≪ n m \ll n m≪n。不难知道, P k f {\bf{P}}_k^f Pkf和 P k a {\bf{P}}_k^a Pka都是对称正定矩阵,因此可以分解为 P k f = Z k f Z k f T {\bf{P}}_k^f = {\bf{Z}}_k^f{\bf{Z}}_k^{f{\rm{T}}} Pkf=ZkfZkfT P k a = Z k a Z k a T {\bf{P}}_k^a = {\bf{Z}}_k^a{\bf{Z}}_k^{a{\rm{T}}} Pka=ZkaZkaT其中 Z k f {\bf{Z}}_k^f Zkf和 Z k a {\bf{Z}}_k^a Zka分别是 P k f {\bf{P}}_k^f Pkf和 P k a {\bf{P}}_k^a Pka的矩阵均方根。根据集合卡尔曼滤波(EnKF)的定义, P k a {\bf{P}}_k^a Pka可以利用维度为 n × m n \times m n×m的扰动矩阵(扰动向量集合) S k a {\bf{S}}_k^a Ska计算,即 P k a = ( m − 1 ) − 1 S k a S k a T {\bf{P}}_k^a = {\left( {m - 1} \right)^{ - 1}}{\bf{S}}_k^a{\bf{S}}_k^{a{\rm{T}}} Pka=(m−1)−1SkaSkaT这里 m m m表示集合数。于是,容易得出 Z k a {\bf{Z}}_k^a Zka和 S k a {\bf{S}}_k^a Ska的关系是 Z k a = ( m − 1 ) − 1 2 S k a {\bf{Z}}_k^a = {\left( {m - 1} \right)^{ - \frac{1}{2}}}{\bf{S}}_k^a Zka=(m−1)−21Ska二者仅相差一个系数。因此,个人理解,无妨认为 Z k a {\bf{Z}}_k^a Zka是广义的分析向量的集合;同理, Z k f {\bf{Z}}_k^f Zkf是广义的预报集合。
2.3 集合卡尔曼滤波(EnKF)
{~~~~~~~}
下面利用
Z
k
f
{\bf{Z}}_k^f
Zkf和
Z
k
a
{\bf{Z}}_k^a
Zka重新计算式(1)和式(2)。如果忽略
Q
k
{{\bf{Q}}_k}
Qk,则有
Z
k
f
Z
k
f
T
=
M
k
Z
k
−
1
a
Z
k
−
1
a
T
M
k
T
=
M
k
Z
k
−
1
a
(
M
k
Z
k
−
1
a
)
T
{\bf{Z}}_k^f{\bf{Z}}_k^{f{\rm{T}}} = {\bf{M}}_k^{}{\bf{Z}}_{k - 1}^a{\bf{Z}}_{k - 1}^{a{\rm{T}}}{\bf{M}}_k^{\rm{T}} = {\bf{M}}_k^{}{\bf{Z}}_{k - 1}^a{\left( {{\bf{M}}_k^{}{\bf{Z}}_{k - 1}^a} \right)^{\rm{T}}}
ZkfZkfT=MkZk−1aZk−1aTMkT=MkZk−1a(MkZk−1a)T于是容易得到预报集合的演化如下
Z
k
f
=
M
k
Z
k
−
1
a
(
4
)
{\bf{Z}}_k^f = {\bf{M}}_k^{}{\bf{Z}}_{k - 1}^a~~~~~~~~(4)
Zkf=MkZk−1a (4)从集合同化的角度看,式(4)意味着
Z
k
f
{\bf{Z}}_k^f
Zkf的每个列向量是
M
k
{\bf{M}}_k^{}
Mk作用到
Z
k
−
1
a
{\bf{Z}}_{k-1}^a
Zk−1a的每个列向量的结果,因此无妨认为
Z
k
f
{\bf{Z}}_k^f
Zkf和
Z
k
a
{\bf{Z}}_k^a
Zka的不同列就代表着不同的集合成员。事实上,模式误差的影响是不可以忽略的。如果有维度为
n
×
q
n \times q
n×q的
Z
k
d
{\bf{Z}}_k^d
Zkd满足
Q
k
=
Z
k
d
Z
k
d
T
{{\bf{Q}}_k} = {\bf{Z}}_k^d{\bf{Z}}_k^{d{\rm{T}}}
Qk=ZkdZkdT那么
P
k
f
{\bf{P}}_k^f
Pkf的一个均方根矩阵可以表示为
Z
k
f
=
M
k
Z
k
−
1
a
Z
k
d
{\bf{Z}}_k^f = {\bf{M}}_k^{}{\bf{Z}}_{k - 1}^a{\bf{Z}}_k^d
Zkf=MkZk−1aZkd其维度是
n
×
m
n \times m
n×m,这里
m
=
m
e
+
q
m = {m_e} + q
m=me+q,
m
e
{m_e}
me表示动力演化的预报扰动的数目。可以通过计算集合的奇异值分解和丢弃方差较小的部分来限制集合的大小(Heemink等,2001)。在分析阶段可以使用具有演化分析误差和模型误差的较大的集合,在动态预报阶段可以使用较小的集合。
{~~~~~~~}
下面重点讨论分析集合。由Evensen(1994)最早提出的集合卡尔曼滤波(EnKF)规定对于每个集合成员都利用相同的同化算法做更新,即
Z
k
a
=
(
I
−
K
k
H
k
)
Z
k
f
(
5
)
{\bf{Z}}_k^a = \left( {{\bf{I}} - {\bf{K}}_k^{}{\bf{H}}_k^{}} \right){\bf{Z}}_k^f~~~~~~~~(5)
Zka=(I−KkHk)Zkf (5)于是,最初的EnKF的分析误差协方差矩阵可表示为
P
k
a
=
Z
k
a
Z
k
a
T
=
(
I
−
K
k
H
k
)
Z
k
f
Z
k
f
T
(
I
−
K
k
H
k
)
T
=
(
I
−
K
k
H
k
)
P
k
f
(
I
−
K
k
H
k
)
T
(
6
)
\begin{aligned} {\bf{P}}_k^a &= {\bf{Z}}_k^a{\bf{Z}}_k^{a{\rm{T}}}\\ &= \left( {{\bf{I}} - {\bf{K}}_k^{}{\bf{H}}_k^{}} \right){\bf{Z}}_k^f{\bf{Z}}_k^{f{\rm{T}}}{\left( {{\bf{I}} - {\bf{K}}_k^{}{\bf{H}}_k^{}} \right)^{\rm{T}}}\\ &= \left( {{\bf{I}} - {\bf{K}}_k^{}{\bf{H}}_k^{}} \right){\bf{P}}_k^f{\left( {{\bf{I}} - {\bf{K}}_k^{}{\bf{H}}_k^{}} \right)^{\rm{T}}}\end{aligned}~~~~~~~~(6)
Pka=ZkaZkaT=(I−KkHk)ZkfZkfT(I−KkHk)T=(I−KkHk)Pkf(I−KkHk)T (6)很明显,式(6)中的
P
k
a
{\bf{P}}_k^a
Pka相较于(2)中的
P
k
a
{\bf{P}}_k^a
Pka被“低估”了
I
−
K
k
H
k
{\bf{I}} - {\bf{K}}_k^{}{\bf{H}}_k^{}
I−KkHk倍。Houtekamer和Mitchell(1998)和Burgers等(1998)不约而同地利用“随机性”的方法对最初的EnKF(Evensen,1994)作了修订,其关键是引入了一组在统计上反应观测误差的“观测集合”。它们的方法相当于将式(5)中的
Z
k
a
{\bf{Z}}_k^a
Zka修订成
Z
k
a
=
(
I
−
K
k
H
k
)
Z
k
f
+
K
k
W
k
(
7
)
{\bf{Z}}_k^a = \left( {{\bf{I}} - {\bf{K}}_k^{}{\bf{H}}_k^{}} \right){\bf{Z}}_k^f + {\bf{K}}_k^{}{\bf{W}}_k^{}~~~~~~~~(7)
Zka=(I−KkHk)Zkf+KkWk (7)其中,
W
k
{\bf{W}}_k^{}
Wk是维度为
p
×
m
p \times m
p×m的矩阵,它由
m
m
m个长度为
p
p
p、方差为
R
k
m
{{{\bf{R}}_k^{}} \over m}
mRk的随机数列向量组成。利用(7)中的
Z
k
a
{\bf{Z}}_k^a
Zka计算分析误差协方差矩阵有
P
k
a
=
Z
k
a
Z
k
a
T
=
(
I
−
K
k
H
k
)
P
k
f
(
I
−
K
k
H
k
)
T
+
K
k
R
k
K
k
T
(
8
)
\begin{aligned} {\bf{P}}_k^a &= {\bf{Z}}_k^a{\bf{Z}}_k^{a{\rm{T}}}\\ &= \left( {{\bf{I}} - {\bf{K}}_k^{}{\bf{H}}_k^{}} \right){\bf{P}}_k^f{\left( {{\bf{I}} - {\bf{K}}_k^{}{\bf{H}}_k^{}} \right)^{\rm{T}}} + {\bf{K}}_k^{}{\bf{R}}_k^{}{\bf{K}}_k^{\rm{T}}\end{aligned}~~~~~~~~(8)
Pka=ZkaZkaT=(I−KkHk)Pkf(I−KkHk)T+KkRkKkT (8)若式(8)中的
P
k
a
{\bf{P}}_k^a
Pka的期望等于式(2)中的
P
k
a
{\bf{P}}_k^a
Pka则可以弥补式(7)中
P
k
a
{\bf{P}}_k^a
Pka的低估。但是,这种扰动观测的方法引入了一个额外的抽样误差来源,这会降低分析误差协方差的准确性,增加低估分析误差协方差的可能性(Mitchell和Houtekamer,2002)。
2.4 “确定性”Kalman SRF
{~~~~~~~}
卡尔曼SRF提供了一个“确定性”更新分析误差协方差矩阵的思路,即将式(2)改写成以下形式
P
k
a
=
Z
k
a
Z
k
a
T
=
[
I
−
P
k
f
H
k
T
(
H
k
P
k
f
H
k
T
+
R
k
)
−
1
H
k
]
P
k
f
=
[
I
−
Z
k
f
Z
k
f
T
H
k
T
(
H
k
Z
k
f
Z
k
f
T
H
k
T
+
R
k
)
−
1
H
k
]
Z
k
f
Z
k
f
T
=
Z
k
f
Z
k
f
T
−
Z
k
f
Z
k
f
T
H
k
T
(
H
k
Z
k
f
Z
k
f
T
H
k
T
+
R
k
)
−
1
H
k
Z
k
f
Z
k
f
T
=
Z
k
f
[
I
−
Z
k
f
T
H
k
T
(
H
k
Z
k
f
Z
k
f
T
H
k
T
+
R
k
)
−
1
H
k
Z
k
f
]
Z
k
f
T
=
Z
k
f
{
I
−
(
H
k
Z
k
f
)
T
[
H
k
Z
k
f
(
H
k
Z
k
f
)
T
+
R
k
]
−
1
H
k
Z
k
f
}
Z
k
f
T
=
Z
k
f
(
I
−
V
k
D
k
−
1
V
k
T
)
Z
k
f
T
(
9
)
\begin{aligned} {\bf{P}}_k^a &= {\bf{Z}}_k^a{\bf{Z}}_k^{a{\rm{T}}}\\ &= \left[ {{\bf{I}} - {\bf{P}}_k^f{\bf{H}}_k^{\rm{T}}{{\left( {{\bf{H}}_k^{}{\bf{P}}_k^f{\bf{H}}_k^{\rm{T}} + {\bf{R}}_k^{}} \right)}^{ - 1}}{\bf{H}}_k^{}} \right]{\bf{P}}_k^f\\ &= \left[ {{\bf{I}} - {\bf{Z}}_k^f{\bf{Z}}_k^{f{\rm{T}}}{\bf{H}}_k^{\rm{T}}{{\left( {{\bf{H}}_k^{}{\bf{Z}}_k^f{\bf{Z}}_k^{f{\rm{T}}}{\bf{H}}_k^{\rm{T}} + {\bf{R}}_k^{}} \right)}^{ - 1}}{\bf{H}}_k^{}} \right]{\bf{Z}}_k^f{\bf{Z}}_k^{f{\rm{T}}}\\ &= {\bf{Z}}_k^f{\bf{Z}}_k^{f{\rm{T}}} - {\bf{Z}}_k^f{\bf{Z}}_k^{f{\rm{T}}}{\bf{H}}_k^{\rm{T}}{\left( {{\bf{H}}_k^{}{\bf{Z}}_k^f{\bf{Z}}_k^{f{\rm{T}}}{\bf{H}}_k^{\rm{T}} + {\bf{R}}_k^{}} \right)^{ - 1}}{\bf{H}}_k^{}{\bf{Z}}_k^f{\bf{Z}}_k^{f{\rm{T}}}\\ &= {\bf{Z}}_k^f\left[ {{\bf{I}} - {\bf{Z}}_k^{f{\rm{T}}}{\bf{H}}_k^{\rm{T}}{{\left( {{\bf{H}}_k^{}{\bf{Z}}_k^f{\bf{Z}}_k^{f{\rm{T}}}{\bf{H}}_k^{\rm{T}} + {\bf{R}}_k^{}} \right)}^{ - 1}}{\bf{H}}_k^{}{\bf{Z}}_k^f} \right]{\bf{Z}}_k^{f{\rm{T}}}\\ &= {\bf{Z}}_k^f\left\{ {{\bf{I}} - {{\left( {{\bf{H}}_k^{}{\bf{Z}}_k^f} \right)}^{\rm{T}}}{{\left[ {{\bf{H}}_k^{}{\bf{Z}}_k^f{{\left( {{\bf{H}}_k^{}{\bf{Z}}_k^f} \right)}^{\rm{T}}} + {\bf{R}}_k^{}} \right]}^{ - 1}}{\bf{H}}_k^{}{\bf{Z}}_k^f} \right\}{\bf{Z}}_k^{f{\rm{T}}}\\ &= {\bf{Z}}_k^f\left( {{\bf{I}} - {\bf{V}}_k^{}{\bf{D}}_k^{ - 1}{\bf{V}}_k^{\rm{T}}} \right){\bf{Z}}_k^{f{\rm{T}}}\end{aligned}~~~~~~~~(9)
Pka=ZkaZkaT=[I−PkfHkT(HkPkfHkT+Rk)−1Hk]Pkf=[I−ZkfZkfTHkT(HkZkfZkfTHkT+Rk)−1Hk]ZkfZkfT=ZkfZkfT−ZkfZkfTHkT(HkZkfZkfTHkT+Rk)−1HkZkfZkfT=Zkf[I−ZkfTHkT(HkZkfZkfTHkT+Rk)−1HkZkf]ZkfT=Zkf{I−(HkZkf)T[HkZkf(HkZkf)T+Rk]−1HkZkf}ZkfT=Zkf(I−VkDk−1VkT)ZkfT (9)其中,定义维度为
m
×
p
m \times p
m×p的矩阵
V
k
≡
(
H
k
Z
k
f
)
T
{\bf{V}}_k^{} \equiv {\left( {{\bf{H}}_k^{}{\bf{Z}}_k^f} \right)^{\rm{T}}}
Vk≡(HkZkf)T和维度为
p
×
p
p \times p
p×p的矩阵
D
k
≡
V
k
T
V
k
+
R
k
{\bf{D}}_k^{} \equiv {\bf{V}}_k^{\rm{T}}{\bf{V}}_k^{} + {\bf{R}}_k^{}
Dk≡VkTVk+Rk。于是,分析集合(分析误差协方差矩阵均方根)可以这样计算
Z
k
a
=
Z
k
f
X
k
U
k
(
10
)
{\bf{Z}}_k^a = {\bf{Z}}_k^f{\bf{X}}_k^{}{\bf{U}}_k^{}~~~~~~~~(10)
Zka=ZkfXkUk (10)其中,
X
k
X
k
T
=
I
−
V
k
D
k
−
1
V
k
T
{\bf{X}}_k^{}{\bf{X}}_k^{\rm{T}} = {\bf{I}} - {\bf{V}}_k^{}{\bf{D}}_k^{ - 1}{\bf{V}}_k^{\rm{T}}
XkXkT=I−VkDk−1VkT,
U
k
{\bf{U}}_k^{}
Uk为任意维度为
m
×
m
m \times m
m×m的正交矩阵。从式(10)可以看出,
Z
k
a
{\bf{Z}}_k^a
Zka是
Z
k
f
{\bf{Z}}_k^f
Zkf的一个线性组合,其求解需要计算
D
k
{\bf{D}}_k^{}
Dk的逆以及
I
−
V
k
D
k
−
1
V
k
T
{\bf{I}} - {\bf{V}}_k^{}{\bf{D}}_k^{ - 1}{\bf{V}}_k^{\rm{T}}
I−VkDk−1VkT的均方根矩阵
X
k
{\bf{X}}_k^{}
Xk,这也是集合均方根滤波(EnSRF)类同化方法需要解决的问题。
3 集合均方根滤波(EnSRF)
3.1 直接(direct)EnSRF
{~~~~~~~} 一种直接求解 X k {\bf{X}}_k^{} Xk的方法需要借助维度为 p × m p \times m p×m的矩阵 Y k {\bf{Y}}_k^{} Yk,其满足 ( H k P k f H k T + R k ) Y k = H k Z k f ( 11 ) \left( {{\bf{H}}_k^{}{\bf{P}}_k^f{\bf{H}}_k^{\rm{T}} + {\bf{R}}_k^{}} \right){\bf{Y}}_k^{} = {\bf{H}}_k^{}{\bf{Z}}_k^f~~~~~~~~(11) (HkPkfHkT+Rk)Yk=HkZkf (11)当 R k {\bf{R}}_k^{} Rk的逆可求时,利用Sherman-Morrison-Woodbury方法(Golub和Van Loan,1996)可以求得 Y k = ( H k P k f H k T + R k ) − 1 H k Z k f = { R k − 1 − R k − 1 H k Z k f [ I + ( H k Z k f ) T R k − 1 ( H k Z k f ) ] − 1 ( H k Z k f ) T R k − 1 } H k Z k f ( 12 ) \begin{aligned}{\bf{Y}}_k^{} &= {\left( {{\bf{H}}_k^{}{\bf{P}}_k^f{\bf{H}}_k^{\rm{T}} + {\bf{R}}_k^{}} \right)^{ - 1}}{\bf{H}}_k^{}{\bf{Z}}_k^f\\ &= \left\{ {{\bf{R}}_k^{ - 1} - {\bf{R}}_k^{ - 1}{\bf{H}}_k^{}{\bf{Z}}_k^f{{\left[ {{\bf{I}} + {{\left( {{\bf{H}}_k^{}{\bf{Z}}_k^f} \right)}^{\rm{T}}}{\bf{R}}_k^{ - 1}\left( {{\bf{H}}_k^{}{\bf{Z}}_k^f} \right)} \right]}^{ - 1}}{{\left( {{\bf{H}}_k^{}{\bf{Z}}_k^f} \right)}^{\rm{T}}}{\bf{R}}_k^{ - 1}} \right\}{\bf{H}}_k^{}{\bf{Z}}_k^f\end{aligned}~~~~~~~~(12) Yk=(HkPkfHkT+Rk)−1HkZkf={Rk−1−Rk−1HkZkf[I+(HkZkf)TRk−1(HkZkf)]−1(HkZkf)TRk−1}HkZkf (12)于是,可以求得 I − V k D k − 1 V k T = I − ( H k Z k f ) T Y k {\bf{I}} - {\bf{V}}_k^{}{\bf{D}}_k^{ - 1}{\bf{V}}_k^{\rm{T}} = {\bf{I}} - {\left( {{\bf{H}}_k^{}{\bf{Z}}_k^f} \right)^{\rm{T}}}{\bf{Y}}_k^{} I−VkDk−1VkT=I−(HkZkf)TYk从而进一步求得 X k {\bf{X}}_k^{} Xk和 Z k a {\bf{Z}}_k^a Zka。
3.2 顺序(serial)EnSRF
{~~~~~~~}
当观测误差是不相关的,观测可以被逐个地或顺序地同化(Houtekamer和Mitchell,2001)。对于单个观测,
p
=
1
p = 1
p=1,
V
k
{\bf{V}}_k^{}
Vk是一个列向量,
D
k
{\bf{D}}_k^{}
Dk是一个标量,那么以下变换成立
I
−
V
k
D
k
−
1
V
k
T
=
I
−
D
k
−
1
V
k
V
k
T
=
(
I
−
β
k
V
k
V
k
T
)
(
I
−
β
k
V
k
V
k
T
)
T
(
13
)
{\bf{I}} - {\bf{V}}_k^{}{\bf{D}}_k^{ - 1}{\bf{V}}_k^{\rm{T}} = {\bf{I}} - {\bf{D}}_k^{ - 1}{\bf{V}}_k^{}{\bf{V}}_k^{\rm{T}} = \left( {{\bf{I}} - \beta _k^{}{\bf{V}}_k^{}{\bf{V}}_k^{\rm{T}}} \right){\left( {{\bf{I}} - \beta _k^{}{\bf{V}}_k^{}{\bf{V}}_k^{\rm{T}}} \right)^{\rm{T}}}~~~~~~~~(13)
I−VkDk−1VkT=I−Dk−1VkVkT=(I−βkVkVkT)(I−βkVkVkT)T (13)
其中
β
k
=
[
D
k
±
(
R
k
D
k
)
1
/
2
]
−
1
\beta _k^{} = {\left[ {{\bf{D}}_k^{} \pm {{\left( {{\bf{R}}_k^{}{\bf{D}}_k^{}} \right)}^{1/2}}} \right]^{ - 1}}
βk=[Dk±(RkDk)1/2]−1每当一个观测被同化后的分析集合就变成
Z
k
a
=
Z
k
f
(
I
−
β
k
V
k
V
k
T
)
(
14
)
{\bf{Z}}_k^a = {\bf{Z}}_k^f\left( {{\bf{I}} - \beta _k^{}{\bf{V}}_k^{}{\bf{V}}_k^{\rm{T}}} \right)~~~~~~~~(14)
Zka=Zkf(I−βkVkVkT) (14)
{~~~~~~~}
习惯上,把Whitaker和Mitchell(2002)提出的方法称为EnSRF(但EnSRF实际上代表一类同化方法),其做法是定义一个新的卡尔曼增益
K
~
k
{\bf{\tilde K}}_k^{}
K~k,满足
(
I
−
K
~
k
H
k
)
P
k
f
(
I
−
K
~
k
H
k
)
T
=
P
k
a
=
(
I
−
K
k
H
k
)
P
k
f
(
15
)
\left( {{\bf{I}} - {\bf{\tilde K}}_k^{}{\bf{H}}_k^{}} \right){\bf{P}}_k^f{\left( {{\bf{I}} - {\bf{\tilde K}}_k^{}{\bf{H}}_k^{}} \right)^{\rm{T}}} = {\bf{P}}_k^a = \left( {{\bf{I}} - {\bf{K}}_k^{}{\bf{H}}_k^{}} \right){\bf{P}}_k^f~~~~~~~~(15)
(I−K~kHk)Pkf(I−K~kHk)T=Pka=(I−KkHk)Pkf (15)对于单个观测,解得
K
~
k
=
[
1
+
R
k
D
k
]
−
1
K
k
=
β
k
Z
k
f
V
k
{\bf{\tilde K}}_k^{} = {\left[ {1 + \sqrt {\frac{{{\bf{R}}_k^{}}}{{{\bf{D}}_k^{}}}} } \right]^{ - 1}}{\bf{K}}_k^{} = \beta _k^{}{\bf{Z}}_k^f{\bf{V}}_k^{}
K~k=[1+DkRk]−1Kk=βkZkfVk,于是
Z
k
a
=
(
I
−
K
~
k
H
k
)
Z
k
f
=
(
I
−
β
k
Z
k
f
V
k
H
k
)
Z
k
f
=
Z
k
f
(
I
−
β
k
V
k
V
k
T
)
(
16
)
\begin{aligned}{\bf{Z}}_k^a &= \left( {{\bf{I}} - {\bf{\tilde K}}_k^{}{\bf{H}}_k^{}} \right){\bf{Z}}_k^f\\ &= \left( {{\bf{I}} - \beta _k^{}{\bf{Z}}_k^f{\bf{V}}_k^{}{\bf{H}}_k^{}} \right){\bf{Z}}_k^f\\ &= {\bf{Z}}_k^f\left( {{\bf{I}} - \beta _k^{}{\bf{V}}_k^{}{\bf{V}}_k^{\rm{T}}} \right)\end{aligned}~~~~~~~~(16)
Zka=(I−K~kHk)Zkf=(I−βkZkfVkHk)Zkf=Zkf(I−βkVkVkT) (16)
3.3 集合变换卡尔曼滤波(ETKF)
{~~~~~~~} 如果 R k − 1 {\bf{R}}_k^{ - 1} Rk−1可求,则利用Sherman-Morrison-Woodbury方法可以对 I − V k D k − 1 V k T {\bf{I}} - {\bf{V}}_k^{}{\bf{D}}_k^{ - 1}{\bf{V}}_k^{\rm{T}} I−VkDk−1VkT变换形式为 ( I + Z k f T H k R k − 1 H k Z k f ) − 1 {\left( {{\bf{I}} + {\bf{Z}}{{_k^f}^{\rm{T}}}{\bf{H}}_k^{}{\bf{R}}_k^{ - 1}{\bf{H}}_k^{}{\bf{Z}}_k^f} \right)^{ - 1}} (I+ZkfTHkRk−1HkZkf)−1,因此可以避免计算 D k − 1 {\bf{D}}_k^{ - 1} Dk−1。Bishop等(2001)提出的集合变换卡尔曼滤波(ETKF)利用了这个变换,并给出分析集合的公式为 Z k a = Z k f C k ( Γ k + I ) − 1 2 ( 17 ) {\bf{Z}}_k^a = {\bf{Z}}_k^f{\bf{C}}_k^{}{\left( {{\bf{\Gamma }}_k^{} + {\bf{I}}} \right)^{ - \frac{1}{2}}}~~~~~~~~(17) Zka=ZkfCk(Γk+I)−21 (17)其中, C k Γ k C k T {\bf{C}}_k^{}{\bf{\Gamma }}_k^{}{\bf{C}}_k^{\rm{T}} CkΓkCkT是 Z k f T H k R k − 1 H k Z k f {\bf{Z}}_k^{f{\rm{T}}}{\bf{H}}_k^{}{\bf{R}}_k^{ - 1}{\bf{H}}_k^{}{\bf{Z}}_k^f ZkfTHkRk−1HkZkf的特征值分解,那么 C k ( Γ k + I ) − 1 C k T {\bf{C}}_k^{}{\left( {{\bf{\Gamma }}_k^{} + {\bf{I}}} \right)^{ - 1}}{\bf{C}}_k^{\rm{T}} Ck(Γk+I)−1CkT就是 I − V k D k − 1 V k T {\bf{I}} - {\bf{V}}_k^{}{\bf{D}}_k^{ - 1}{\bf{V}}_k^{\rm{T}} I−VkDk−1VkT的特征值分解,所以 C k ( Γ k + I ) − 1 2 {\bf{C}}_k^{}{\left( {{\bf{\Gamma }}_k^{} + {\bf{I}}} \right)^{ - \frac{1}{2}}} Ck(Γk+I)−21实际上就是 I − V k D k − 1 V k T {\bf{I}} - {\bf{V}}_k^{}{\bf{D}}_k^{ - 1}{\bf{V}}_k^{\rm{T}} I−VkDk−1VkT的一个均方根矩阵,也就是 X k {\bf{X}}_k^{} Xk。
3.4 集合调整卡尔曼滤波(EAKF)
{~~~~~~~} Anderson(2001)给出了另一种分析更新的公式 Z k a = A k Z k f ( 18 ) {\bf{Z}}_k^a = {\bf{A}}_k^{}{\bf{Z}}_k^f~~~~~~~~(18) Zka=AkZkf (18)其中, A k {\bf{A}}_k^{} Ak被称为调整矩阵,其定义为 A k ≡ F k G k C ~ k ( Γ ~ k + I ) − 1 2 G k − 1 F k T {\bf{A}}_k^{} \equiv {\bf{F}}_k^{}{\bf{G}}_k^{}{\bf{\tilde C}}_k^{}{\left( {{\bf{\tilde \Gamma }}_k^{} + {\bf{I}}} \right)^{ - \frac{1}{2}}}{\bf{G}}_k^{ - 1}{\bf{F}}_k^{\rm{T}} Ak≡FkGkC~k(Γ~k+I)−21Gk−1FkT这里 F k G k 2 F k T {\bf{F}}_k^{}{\bf{G}}_k^2{\bf{F}}_k^{\rm{T}} FkGk2FkT是 P k f {\bf{P}}_k^f Pkf的特征值分解,正交矩阵 C ~ k {\bf{\tilde C}}_k^{} C~k应满足 C ~ k T G k F k T H k T R k − 1 H k F k G k C ~ k = Γ ~ k {\bf{\tilde C}}_k^{\rm{T}}{\bf{G}}_k^{}{\bf{F}}_k^{\rm{T}}{\bf{H}}_k^{\rm{T}}{\bf{R}}_k^{ - 1}{\bf{H}}_k^{}{\bf{F}}_k^{}{\bf{G}}_k^{}{\bf{\tilde C}}_k^{} = {\bf{\tilde \Gamma }}_k^{} C~kTGkFkTHkTRk−1HkFkGkC~k=Γ~k其中 Γ ~ k {\bf{\tilde \Gamma }}_k^{} Γ~k是对角矩阵。如果令 Γ ~ k = Γ k {\bf{\tilde \Gamma }}_k^{} = {\bf{\Gamma }}_k^{} Γ~k=Γk C ~ k = G k − 1 F k T Z k f C k {\bf{\tilde C}}_k^{} = {\bf{G}}_k^{ - 1}{\bf{F}}_k^{\rm{T}}{\bf{Z}}_k^f{\bf{C}}_k^{} C~k=Gk−1FkTZkfCk那么调整矩阵可以表示为 A k = Z k f C k ( Γ k + I ) − 1 2 G k − 1 F k T {\bf{A}}_k^{} = {\bf{Z}}_k^f{\bf{C}}_k^{}{\left( {{\bf{\Gamma }}_k^{} + {\bf{I}}} \right)^{ - \frac{1}{2}}}{\bf{G}}_k^{ - 1}{\bf{F}}_k^{\rm{T}} Ak=ZkfCk(Γk+I)−21Gk−1FkT于是,式(18)就变成 Z k a = Z k f C k ( Γ k + I ) − 1 2 G k − 1 F k T Z k f ( 19 ) {\bf{Z}}_k^a = {\bf{Z}}_k^f{\bf{C}}_k^{}{\left( {{\bf{\Gamma }}_k^{} + {\bf{I}}} \right)^{ - \frac{1}{2}}}{\bf{G}}_k^{ - 1}{\bf{F}}_k^{\rm{T}}{\bf{Z}}_k^f~~~~~~~~(19) Zka=ZkfCk(Γk+I)−21Gk−1FkTZkf (19)比较式(17)和式(19)会发现,EAKF的分析集合相当于在ETKF的分析集合上右乘 G k − 1 F k T Z k f {\bf{G}}_k^{ - 1}{\bf{F}}_k^{\rm{T}}{\bf{Z}}_k^f Gk−1FkTZkf。事实上,这个 G k − 1 F k T Z k f {\bf{G}}_k^{ - 1}{\bf{F}}_k^{\rm{T}}{\bf{Z}}_k^f Gk−1FkTZkf是正交矩阵,也是 Z k f {\bf{Z}}_k^f Zkf的右奇异向量。因此, C k ( Γ k + I ) − 1 2 G k − 1 F k T Z k f {\bf{C}}_k^{}{\left( {{\bf{\Gamma }}_k^{} + {\bf{I}}} \right)^{ - \frac{1}{2}}}{\bf{G}}_k^{ - 1}{\bf{F}}_k^{\rm{T}}{\bf{Z}}_k^f Ck(Γk+I)−21Gk−1FkTZkf同样是 I − V k D k − 1 V k T {\bf{I}} - {\bf{V}}_k^{}{\bf{D}}_k^{ - 1}{\bf{V}}_k^{\rm{T}} I−VkDk−1VkT的一个均方根矩阵。
4 总结
{~~~~~~~} 集合预测/同化方法使用预测和分析误差协方差矩阵的低秩集合表示,这些集合是误差协方差矩阵的矩阵平方根,因此集合数据同化方法可以统一地看作是平方根滤波(SRF)。同化观测后,可以“随机性”地构建分析集合(如EnKF)或“确定性”构建分析集合(如EnSRF)。二者最大的区别是前者需要构造观测集合,而后者消除了生成观测集合带来的抽样误差的影响,而且EnSRF在某些应用实例中可以得到比EnKF更精确的状态估计(Whitaker和Hamill,2002;Anderson,2001)。估计理论中SRF的不唯一性已被广泛用于设计具有理想计算和数值性质的滤波方法。一个悬而未决的问题是,是否存在一种特定的SRF,其得到的集合特性能够优于其他方法,或许比较的关键就是计算成本。例如,可以选择一种保留预测集合的高阶非高斯统计量的分析更新方案。这个问题只能通过在现实环境中对不同方法的详细比较来回答。在现实环境中,同化系统的其他细节,如系统误差建模或数据质量控制,可能被证明同样重要。