2.1 递推线性最小方差估计框架

考虑如下非线性系统

式中，f（·，·）∈Rn为已知的状态函数，x（k）∈Rn为k时刻系统状态，h（·，·）∈Rm为已知的传感器观测函数，z（k）∈Rm为传感器观测数据，为系统噪声，为传感器观测噪声。假设w（k）和v（k）是零均值、方差阵分别为Q_w和R且相互独立噪声，即

式中，E为均值号，T为转置号，δ_tk=0（t≠k），δ（·）是狄拉克函数（Dirac Delta function）。

问题是根据已知观测数据Z_0～k={z（0）～z（k）}，求解状态x（k）的估计。

2.1.1 射影定理

定义2.1[14]基于m×1维随机变量z∈Rm的对n×1维随机变量x∈Rn的线性估计记为

若估计极小化性能指标为

则称为随机变量x的线性最小方差估计，式中E为均值号，T为转置号。

定理2.1[14] 基于z∈Rm对随机变量x∈Rn的线性最小方差估计公式为

其中假设Ex、Ez、P_xz、P_zz均存在。

证明：将式（2-4）代入式（2-5）有

令有

所以有

将式（2-9）代入式（2-7）并定义

可有关系

令，应用矩阵迹求导公式[152]，并整理有

证毕。

推论2.1[14]

证明：由式（2-6）有

证毕。

推论2.2[14]

证明：将式（2-6）代入式（2-16）左边，有

证毕。

推论2.3[14] 与z不相关。

证明：

证毕。

定义2.2[14] 与z不相关称为与z正交（垂直），记为，并称为x在z上的射影，记为。

定义2.3[14] 由随机变量z∈Rm张成的线性流形（线性空间）定义为如下形式的随机变量z∈Rn的集合

推论2.4[14]

证明：

证毕。

定义2.4[14] 设随机变量x∈Rn，随机变量z（1），…，z（k）∈Rm，引入合成随机变量ϖ为

由z（1），…，z（k）∈Rm张成的线性流形L（z（1），…，z（k））定义为

引入分块矩阵

则有

定义2.5[14] 基于随机变量z（1），…，z（k）∈Rm对随机变量x∈Rn的线性最小方差估计定义为

也称为x在线性流形或者L（z（1），…，z（k））上的射影。

推论2.5[14] 设x∈Rn为零均值随机变量，z（1），…，z（k）∈Rm为零均值、互不相关（正交）的随机变量，则有

证明：

推论2.6[14]设随机变量x∈Rp，y∈Rq，随机变量（Ax+By）∈Rn，A∈Rn×p，B∈Rn×q，随机变量z∈Rm，则有

推论2.7[14] 设随机变量x∈Rn，随机变量z∈Rm，则有关系

其中x的分量形式为

证明：

即得到式（2-28）。证毕。

2.1.2 新息序列

定义2.6[14] 设z（1），…，z（k），…∈Rm是存在2阶矩的随机序列，它的新息序列定义为

其中z（k）的一步最优预报估值为

因而新息序列定义为

其中规定，这样可以保证E{ε（1）}=0。

定理2.2[14]新息序列ε（k）是零均值白噪声。

证明：由新息序列定义式（2-33），有

由推论2.1，可得

设i≠j，可以设i＞j，又由于ε（i）⊥L（z（1），…，z（i-1）），且有L（z（1），…，z（j））⊂L（z（1），…，z（i-1）），因此ε（i）⊥L（z（1），…，z（j））。

又因为ε（j）=z（j）-zˆ（j|j-1）∈L（z（1），…，z（j）），因而ε（i）⊥ε（j），即

故ε（i）是白噪声。证毕。

定理2.3[14]新息序列ε（k）与原序列z（k）含有相同的统计信息，即（z（1），…，z（k））与（ε（1），…，ε（k））张成相同的线性流形，即

证明：由式（2-6）和式（2-32），每个ε（k）是z（1），…，z（k）的线性组合，这里引出

从而有

下面用数学归纳法证明z（k）∈L（ε（1），…，ε（k））。

由式（2-32），有

故有

从而有式（2-37）成立。证毕。

推论2.8[14] 设随机变量x∈Rn，则有

定理2.4[14]（递推射影公式）设随机变量x∈Rn，随机序列z（1），…，z（k），…∈Rm，且它们存在2阶矩，则有递推射影公式

证明：引入合成向量

有E{ε（i）}=0。

由式（2-42）和式（2-6），有

证毕。

2.1.3 递推线性最小方差滤波框架

2.1.2节，在最小方差意义下，递推射影定理被给出。本节我们将给出一种具体的滤波估计框架。

定理2.5[116] 对系统式（2-1）和式（2-2），局部滤波器有如下递推滤波框架

其中滤波增益为

滤波误差方差阵为

其中

预报误差方差阵P（k+1|k）为

证明：根据最小方差估计理论，一步预测是状态的条件数学期望，即有

可以得到式（2-48）。

即有

然后可以得到式（2-49）。

由预报误差协方差阵P_xz（k+1|k）定义有

因为假设v（k）是具有零均值且独立的Gauss噪声，所以得到式（2-50）。

由观测误差协方差矩阵P_z（k+1|k）定义有

类似于P_xz（k+1|k），式（2-56）可以写为式（2-51）。

由预报误差方差阵P（k+1|k）定义有

可得式（2-52）。

将式（2-45）代入滤波误差协方差矩阵定义式，整理得

基于最小方差估计准则，，可以得到式（2-46）。证毕。

2.1.4 Kalman滤波器

滤波是去除噪声还原真实数据的一种数据处理技术。Kalman滤波在观测方差已知的情况下能够从一系列存在观测噪声的数据中，估计动态系统的状态。由于它便于计算机编程实现，并能够对现场采集的数据进行实时更新和处理，因此Kalman滤波是目前应用最为广泛的滤波算法，在通信、导航、制导与控制等领域得到了较好的应用。

考虑如下多传感器定常线性随机系统[14]

其中x（k）∈Rn为状态，z（k）∈Rm为第j个传感器的观测，为观测白噪声，w（k）∈Rr为输入白噪声，Φ、Γ、Η为已知的适当维常阵。

假设1w（k）∈Rr为相互独立的，方差阵各为Q_w和R的互不相关的白噪声，且噪声均值和方差统计为

假设2x（0）不相关于w（k）和v（k）

Kalman滤波问题是：基于观测Z_0～k={z（0）～z（k）}，求解状态x（j）的线性最小方差估计，它极小化性能指标为

对于j=k，j＞k，j＜k，分别称为Kalman滤波器、预报器和平滑器。下面应用射影定理推导Kalman滤波器。

定理2.6[14] 系统式（2-59）和式（2-60）在假设1和假设2下，经典Kalman滤波方程组如下：

证明：由递推射影公式（2-43）有

令

则有式（2-65）成立。称K（k+1）为Kalman滤波增益。对式（2-59）两边取射影有

由式（2-59）迭代有

将式（2-75）代入式（2-60）有

引出如下关系

由假设1、假设2和式（2-77）有

应用式（2-6）和E{w（k）}=q可得

于是得到式（2-67）成立。

对式（2-60）取射影有

由假设1、假设2和式（2-77）有

应用式（2-6）和E{v（k）}=r可得

于是有

将式（2-83）代入式（2-33），得到新息表达式（2-66）成立。

记滤波和预报估值误差及方差阵为

则由式（2-60）和式（2-84）有新息表达式

且由式（2-59）和式（2-67）有

由式（2-65）有

将式（2-88）代入式（2-90）有

其中I_n为n×n单位阵。因为

故有

这里引出

于是由式（2-89）得到式（2-69）成立。

又因为

故有

这引出

于是由式（2-88）有

且由式（2-91）可得

由式（2-88）有

由射影正交性有

且存在关系，于是由式（2-100）有

将式（2-102）和式（2-98）代入式（2-73），则式（2-68）成立。

将式（2-68）代入式（2-99）并化简整理得式（2-70）成立。证毕。

Kalman滤波递推算法框图如图2-1所示。

2.1.5 ARMA新息模型

由式（2-59）和式（2-60）有

其中I_n为n×n单位阵，q-1为单位滞后算子，q-1x（k）=x（k-1）。引入左素分解

其中多项式矩阵A（q-1）和B（q-1）有形式

将式（2-104）代入式（2-103）引出自回归滑动平均（Autoregressive Moving Aerage，ARMA）新息模型

其中D（q-1）是稳定的，新息ε（k）∈Rm是零均值、方差阵为Q_ε的白噪声，D（q-1）和Q_ε可用G-W（Gevers-Wouters）算法[14]求得。

2.1.6 基于ARMA新息模型的稳态Kalman滤波器

定理2.7[14] 系统式（2-59）和式（2-60）在假设1和假设2下，基于现代时间序列的稳态Kalman滤波算法如下：

其中M_i可递推计算为

其中规定M₀ =I _m，M_t=0（t＜0），D_t=0（t＞n_d）。

证明：见文献[14]。