动态线性系统的卡尔曼滤波

连续线性系统的数学模型

连续线性系统的状态方程和观测方程

为了说明动态系统的数学模型，我们首先举一个例子。卫星在空间飞行，在地球各地的卫星观测站对卫星的瞬时位置——角位置偏差进行观测，可以通过这些观测值来求定和预报卫星飞行轨道的偏差。由于卫星在飞行过程中受大气阻力的影响，要产生阻力加速度，其主要影响是改变卫星沿轨道切线方向的运动，使卫星轨道产生偏差，而不能按标称轨道运行。设在 \(t\) 时刻大气阻力引起的卫星沿轨道切线方向的角位置偏差为 \(X_{1}(t)\)，角速度偏差为 \(X_{2}(t)\)，随机阻力的角加速度为 \(\varOmega(t)\)，则它们的相互关系可以用如下的线性微分方程描述： \[\left.\begin{aligned} \frac{\mathrm{d}X_{1}(t)}{\mathrm{d}t}&=\dot{X}_{1}(t)=X_{2}(t)\\ \frac{\mathrm{d}X_{2}(t)}{\mathrm{d}t}&=\dot{X}_{2}(t)=\varOmega(t) \end{aligned}\right\}\] 令 \[\bm{X}(t)=\begin{bmatrix}X_{1}(t)\\ X_{2}(t)\end{bmatrix},\quad \dot{\bm{X}}(t)=\begin{bmatrix}\dot{X}_{1}(t)\\ \dot{X}_{2}(t)\end{bmatrix}\] 则上面的微分方程可写为 \[\dot{\bm{X}}(t)=\begin{bmatrix}0&1\\ 0&0\end{bmatrix}\bm{X}(t)+\begin{bmatrix}0\\ 1\end{bmatrix}\varOmega(t) \tag{4-1-1}\] 又设角位置偏差的观测值为 \(L(t)\)，相应的观测噪声为 \(\varDelta(t)\)，则观测值 \(L(t)\) 与 \(\bm{X}(t)\) 有关系 \[L(t)=\begin{bmatrix}1&0\end{bmatrix}\bm{X}(t)+\varDelta(t) \tag{4-1-2}\] 因此，求定和预报卫星轨道偏差的问题，也就是由上面的方程 (4-1-1) 和 (4-1-2) 估计随时间不断变化的随机向量 \(\bm{X}(t)\)。

通常称要估计的随时间不断变化的随机向量 \(\bm{X}(t)\) 为动态系统在 \(t\) 时刻的“状态”向量。严格地说，动态系统的“状态”是指全面确定动态系统运动状况的最少的一组参数。一个动态系统的状态是遵循一定的物理规律而变化的，因而不同时刻对状态所作的观测也是互相联系的，对于连续时间系统，在一般情况下，状态随时间 \(t\) 变化的规律，是用一个具有随机初始状态的向量微分方程来描述，称为状态方程或动态方程，如 (4-1-1) 式。

观测向量 \(\bm{L}(t)\) 与状态向量 \(\bm{X}(t)\) 之间一般也存在某种函数关系，仍称为观测方程，也称为输出方程，如 (4-1-2) 式。

动态方程和观测方程一般可表示成 \[\left.\begin{aligned} \dot{\bm{X}}(t)&=f(\bm{X}(t),\bm{U}(t),\bm{\varOmega}(t))\\ \bm{L}(t)&=g(\bm{X}(t),\bm{U}(t),\bm{\Delta}(t)) \end{aligned}\right\}\] 式中 \(f(\cdot)\) 和 \(g(\cdot)\) 分别是已知的 \(m\) 维和 \(n\) 维向量函数，\(\underset{m\times 1}{\bm{X}}(t)\) 是动态系统的状态向量，\(\underset{n\times 1}{\bm{L}}(t)\) 是观测向量或输出向量，\(\underset{r\times 1}{\bm{U}}(t)\) 是 \(r\) 维控制（或输入）向量，\(\underset{p\times 1}{\bm{\varOmega}}(t)\) 是 \(p\) 维动态噪声（或干扰）向量，\(\underset{n\times 1}{\bm{\Delta}}(t)\) 是 \(n\) 维观测噪声（或输出噪声）向量。状态方程和观测方程也可以称为动态系统的函数模型。

通常考虑函数 \(f\)、\(g\) 分别是 \(\bm{X}(t)\)、\(\bm{U}(t)\)、\(\bm{\varOmega}(t)\) 和 \(\bm{X}(t)\)、\(\bm{U}(t)\) 的线性函数，并且观测噪声 \(\bm{\Delta}(t)\) 是一种可加噪声。此时可得到函数模型 \[\left.\begin{aligned} \dot{\bm{X}}(t)&=\bm{A}(t)\bm{X}(t)+\bm{C}(t)\bm{U}(t)+\bm{F}(t)\bm{\varOmega}(t)\\ \bm{L}(t)&=\bm{B}(t)\bm{X}(t)+\bm{G}(t)\bm{U}(t)+\bm{\Delta}(t) \end{aligned}\right\} \tag{4-1-3}\] 式中 \(\underset{m\times m}{\bm{A}}(t)\)、\(\underset{m\times r}{\bm{C}}(t)\)、\(\underset{m\times p}{\bm{F}}(t)\)、\(\underset{n\times m}{\bm{B}}(t)\) 和 \(\underset{n\times r}{\bm{G}}(t)\) 都是随时间连续变化的系数矩阵。一般称由 (4-1-3) 式表示的动态系统为连续线性系统。

如果不考虑系统的确定性输入 \(\bm{U}(t)\)，即当 \(\bm{C}(t)=0\) 和 \(\bm{G}(t)=0\) 时，(4-1-3) 式可简化为 \[\left.\begin{aligned} \dot{\bm{X}}(t)&=\bm{A}(t)\bm{X}(t)+\bm{F}(t)\bm{\varOmega}(t)\\ \bm{L}(t)&=\bm{B}(t)\bm{X}(t)+\bm{\Delta}(t) \end{aligned}\right\} \tag{4-1-4}\] 当 (4-1-3) 式和 (4-1-4) 式中的系数矩阵均与时间 \(t\) 无关时，则变为 \[\left.\begin{aligned} \dot{\bm{X}}(t)&=\bm{A}\bm{X}(t)+\bm{C}\bm{U}(t)+\bm{F}\bm{\varOmega}(t)\\ \bm{L}(t)&=\bm{B}\bm{X}(t)+\bm{G}\bm{U}(t)+\bm{\Delta}(t) \end{aligned}\right\} \tag{4-1-5}\] \[\left.\begin{aligned} \dot{\bm{X}}(t)&=\bm{A}\bm{X}(t)+\bm{F}\bm{\varOmega}(t)\\ \bm{L}(t)&=\bm{B}\bm{X}(t)+\bm{\Delta}(t) \end{aligned}\right\} \tag{4-1-6}\] 称它们为常系数系统的函数模型。

状态方程的解

状态方程 (4-1-3) 是线性微分方程，它的一般解等于它所对应的齐次方程的通解和非齐次方程的特解之和。为此，先求齐次状态方程 \[\dot{\bm{X}}(t)=\bm{A}(t)\bm{X}(t) \tag{4-1-7}\] 的通解。

设 \(\bm{X}(t_{0})\) 为任意 \(m\) 阶向量（\(t\geq t_{0}\)），并设解的形式是 \[\underset{m\times 1}{\bm{X}}(t)=\underset{m\times m}{\overline{\bm{X}}}(t)\,\underset{m\times 1}{\bm{X}}(t_{0}) \tag{4-1-8}\] 其中 \(\overline{\bm{X}}(t)\) 是未知矩阵，则有 \[\dot{\bm{X}}(t)=\dot{\overline{\bm{X}}}(t)\,\bm{X}(t_{0}) \tag{4-1-9}\] 将上两式代入 (4-1-7) 式，得 \[\left[\dot{\overline{\bm{X}}}(t)-\bm{A}(t)\overline{\bm{X}}(t)\right]\bm{X}(t_{0})=0 \tag{4-1-10}\] 因为 \(\bm{X}(t_{0})\) 是任意向量，所以对所有的 \(t\geq t_{0}\)，仅当 \(\overline{\bm{X}}(t)\) 满足微分方程 \[\dot{\overline{\bm{X}}}(t)=\bm{A}(t)\overline{\bm{X}}(t) \tag{4-1-11}\] 时，(4-1-10) 式才能成立，而且当 \(t=t_{0}\) 时，必有 \[\bm{X}(t_{0})=\overline{\bm{X}}(t_{0})\bm{X}(t_{0})\] 或写为 \[\left[\bm{E}-\overline{\bm{X}}(t_{0})\right]\bm{X}(t_{0})=0\] 即有 \[\overline{\bm{X}}(t_{0})=\bm{E} \tag{4-1-12}\] 显然，(4-1-12) 式是 (4-1-11) 的初始条件。因此，齐次方程 (4-1-7) 的通解可写为 \[\bm{X}(t)=\overline{\bm{X}}(t)\bm{X}(t_{0})\] 式中 \(\overline{\bm{X}}(t)\) 应满足 (4-1-10) 和 (4-1-12) 两式。上式也就是 (4-1-8) 式。

下面求 (4-1-3) 式的特解。设求得的特解为 \[\bm{X}(t)=\overline{\bm{X}}(t)\bm{Y}(t) \tag{4-1-13}\] 式中 \(\underset{m\times 1}{\bm{Y}}(t)\) 是未知的 \(m\) 维向量。对上式求导得 \[\dot{\bm{X}}(t)=\dot{\overline{\bm{X}}}(t)\bm{Y}(t)+\overline{\bm{X}}(t)\dot{\bm{Y}}(t) \tag{4-1-14}\] 将上两式代入 (4-1-3) 式，有 \[\dot{\overline{\bm{X}}}(t)\bm{Y}(t)+\overline{\bm{X}}(t)\dot{\bm{Y}}(t) =\bm{A}(t)\overline{\bm{X}}(t)\bm{Y}(t)+\bm{C}(t)\bm{U}(t)+\bm{F}(t)\bm{\varOmega}(t)\] 顾及 (4-1-11) 式，上式可化为 \[\overline{\bm{X}}(t)\dot{\bm{Y}}(t)=\bm{C}(t)\bm{U}(t)+\bm{F}(t)\bm{\varOmega}(t) \tag{4-1-15}\] 可以证明，\(\overline{\bm{X}}(t)\) 是非奇异阵，\(\overline{\bm{X}}^{-1}(t)\) 存在，所以 \[\dot{\bm{Y}}(t)=\overline{\bm{X}}^{-1}(t)\left[\bm{C}(t)\bm{U}(t)+\bm{F}(t)\bm{\varOmega}(t)\right] \tag{4-1-16}\] 对上式积分得 \[\bm{Y}(t)=\int_{t_{0}}^{t}\overline{\bm{X}}^{-1}(\tau)\left[\bm{C}(\tau)\bm{U}(\tau)+\bm{F}(\tau)\bm{\varOmega}(\tau)\right]\mathrm{d}\tau \tag{4-1-17}\] 将 (4-1-17) 式代入 (4-1-13) 式，则得 (4-1-3) 式特解为 \[\bm{X}(t)=\overline{\bm{X}}(t)\int_{t_{0}}^{t}\overline{\bm{X}}^{-1}(\tau)\left[\bm{C}(\tau)\bm{U}(\tau)+\bm{F}(\tau)\bm{\varOmega}(\tau)\right]\mathrm{d}\tau \tag{4-1-18}\]

将 (4-1-8) 与 (4-1-18) 式合并后，可得到一般解为 \[\bm{X}(t)=\overline{\bm{X}}(t)\bm{X}(t_{0}) +\overline{\bm{X}}(t)\int_{t_{0}}^{t}\overline{\bm{X}}^{-1}(\tau)\left[\bm{C}(\tau)\bm{U}(\tau)+\bm{F}(\tau)\bm{\varOmega}(\tau)\right]\mathrm{d}\tau \tag{4-1-19}\] 令 \[\underset{m\times m}{\bm{\varPhi}}(t,\tau)=\overline{\bm{X}}(t)\overline{\bm{X}}^{-1}(\tau) \tag{4-1-20}\] 称 \(\underset{m\times m}{\bm{\varPhi}}(t,\tau)\) 为系统转移矩阵。因 \(\overline{\bm{X}}(t_{0})=\bm{E}\)，所以 \[\bm{\varPhi}(t,t_{0})=\overline{\bm{X}}(t)\overline{\bm{X}}^{-1}(t_{0})=\overline{\bm{X}}(t)\] 故 (4-1-19) 式可写为 \[\bm{X}(t)=\bm{\varPhi}(t,t_{0})\bm{X}(t_{0}) +\int_{t_{0}}^{t}\bm{\varPhi}(t,\tau)\left[\bm{C}(\tau)\bm{U}(\tau)+\bm{F}(\tau)\bm{\varOmega}(\tau)\right]\mathrm{d}\tau \quad(t\geq t_{0}) \tag{4-1-21}\] 上式就是状态方程 (4-1-3) 的解。显然，对于不考虑系统的确定性输入的状态方程 (4-1-4)，其解为

\[\bm{X}(t)=\bm{\varPhi}(t,t_{0})\bm{X}(t_{0}) +\int_{t_{0}}^{t}\bm{\varPhi}(t,\tau)\bm{F}(\tau)\bm{\varOmega}(\tau)\,\mathrm{d}\tau \tag{4-1-22}\]

将 (4-1-20) 式对时间 \(t\) 求导数，可得 \[\dot{\bm{\varPhi}}(t,\tau)=\dot{\overline{\bm{X}}}(t)\overline{\bm{X}}^{-1}(\tau) =\bm{A}(t)\overline{\bm{X}}(t)\overline{\bm{X}}^{-1}(\tau)\] 所以有 \[\dot{\bm{\varPhi}}(t,\tau)=\bm{A}(t)\bm{\varPhi}(t,\tau) \tag{4-1-23}\] 对于一切 \(t\geq t_{0}\)，由 (4-1-20) 式可知 \[\bm{\varPhi}(t,t)=\overline{\bm{X}}(t)\overline{\bm{X}}^{-1}(t)=\bm{E} \tag{4-1-24}\] 因此，以 (4-1-24) 式为初始条件，解微分方程 \[\dot{\bm{\varPhi}}(t,t_{0})=\bm{A}(t)\bm{\varPhi}(t,t_{0}) \tag{4-1-25}\] 再以 \(\tau\) 代换 \(t_{0}\) 就可以求得 \(\bm{\varPhi}(t,\tau)\)。

由状态转移矩阵的定义 (4-1-20) 式，还可得出它的两个重要性质：

(1) 对所有的 \(t_{1}\)，\(t_{2}\)，\(\tau\geq t_{0}\)，有 \[\bm{\varPhi}(t_{2},t_{1})=\bm{\varPhi}(t_{2},\tau)\bm{\varPhi}(\tau,t_{1}) \tag{4-1-26}\]

(2) 对所有的 \(t_{1}\)，\(t_{2}\geq t_{0}\)，有 \[\bm{\varPhi}^{-1}(t_{1},t_{2})=\bm{\varPhi}(t_{2},t_{1}) \tag{4-1-27}\]

例 4-1-1

试求状态方程 (4-1-1) 之解。

为求此状态方程之解，需要利用初始条件 \(\bm{\varPhi}(t,t_{0})=\bm{E}\)，解矩阵微分方程 \[\dot{\bm{\varPhi}}(t,t_{0})=\begin{bmatrix}0&1\\ 0&0\end{bmatrix}\bm{\varPhi}(t,t_{0})\] 也就是解下列微分方程组 \[\left.\begin{aligned} \dot{\varPhi}_{11}(t,t_{0})&=\varPhi_{21}(t,t_{0})\\ \dot{\varPhi}_{12}(t,t_{0})&=\varPhi_{22}(t,t_{0})\\ \dot{\varPhi}_{21}(t,t_{0})&=0\\ \dot{\varPhi}_{22}(t,t_{0})&=0 \end{aligned}\right\}\] 由后面两个方程和初始条件容易得出 \[\left.\begin{aligned} \varPhi_{21}(t,t_{0})&=0\\ \varPhi_{22}(t,t_{0})&=1 \end{aligned}\right\}\] 而前面两个方程可变为 \[\left.\begin{aligned} \dot{\varPhi}_{11}(t,t_{0})&=0\\ \dot{\varPhi}_{12}(t,t_{0})&=1 \end{aligned}\right\}\] 它们的解为 \[\left.\begin{aligned} \varPhi_{11}(t,t_{0})&=1\\ \varPhi_{12}(t,t_{0})&=t-t_{0} \end{aligned}\right\}\] 所以，状态转移矩阵为 \[\bm{\varPhi}(t,\tau)=\begin{bmatrix}1&t-\tau\\ 0&1\end{bmatrix}\] 状态方程之解为 \[\begin{aligned} \bm{X}(t)&=\begin{bmatrix}1&t-t_{0}\\ 0&1\end{bmatrix}\bm{X}(t_{0}) +\int_{t_{0}}^{t}\begin{bmatrix}1&t-\tau\\ 0&1\end{bmatrix}\begin{bmatrix}0\\ 1\end{bmatrix}\varOmega(\tau)\,\mathrm{d}\tau\\ &=\begin{bmatrix}1&t-t_{0}\\ 0&1\end{bmatrix}\bm{X}(t_{0}) +\int_{t_{0}}^{t}\begin{bmatrix}t-\tau\\ 1\end{bmatrix}\varOmega(\tau)\,\mathrm{d}\tau \end{aligned}\]

连续线性系统的随机模型

在实际工作中，常见的连续线性系统是动态噪声和观测噪声均为白噪声，且它们之间完全不相关的情况。对这种白噪声情况下的连续线性系统，通常是假定它有以下随机模型：

(1) 系统的动态噪声 \(\bm{\varOmega}(t)\,(t\geq 0)\) 和观测噪声 \(\bm{\Delta}(t)\,(t\geq t_{0})\) 是零均值白噪声或高斯白噪声随机向量过程，即 \[\left.\begin{aligned} E[\bm{\varOmega}(t)]&=\bm{0}\\ E[\bm{\Delta}(t)]&=\bm{0} \end{aligned}\right\} \tag{4-1-28}\] \[\left.\begin{aligned} \operatorname{cov}[\bm{\varOmega}(t),\bm{\varOmega}(\tau)]&=\bm{D}_{\varOmega}(t)\,\delta(t-\tau)\\ \operatorname{cov}[\bm{\Delta}(t),\bm{\Delta}(\tau)]&=\bm{D}_{\Delta}(t)\,\delta(t-\tau) \end{aligned}\right\} \tag{4-1-29}\] 式中 \(\bm{D}_{\varOmega}(t)\) 是随时间连续变化的对称非负定阵，而 \(\bm{D}_{\Delta}(t)\) 是随时间连续变化的对称正定阵，\(\delta(0)\) 是 Dirac-\(\delta\) 函数。

(2) 动态噪声与观测噪声完全不相关，即 \[\operatorname{cov}[\bm{\varOmega}(t),\bm{\Delta}(\tau)]=\bm{0} \tag{4-1-30}\]

(3) 系统的初始状态 \(\bm{X}(t_{0})\) 是具有某种已知分布或正态分布的随机向量，其均值与方差阵为 \[\left.\begin{aligned} E[\bm{X}(t_{0})]&=\bm{\mu}_{X}(t_{0})\\ \operatorname{var}[\bm{X}(t_{0})]&=\bm{D}_{X}(t_{0}) \end{aligned}\right\} \tag{4-1-31}\]

(4) 动态噪声、观测噪声均与初始状态 \(\bm{X}(t_{0})\) 不相关，即对于一切 \(t\geq t_{0}\) 有 \[\left.\begin{aligned} \operatorname{cov}[\bm{X}(t_{0}),\bm{\varOmega}(t)]&=\bm{0}\\ \operatorname{cov}[\bm{X}(t_{0}),\bm{\Delta}(t)]&=\bm{0} \end{aligned}\right\} \tag{4-1-32}\]

上面是常见的白噪声下的连续线性系统的随机模型，如果将 (4-1-28) 和 (4-1-30) 修改为 \[\left.\begin{aligned} E[\bm{\varOmega}(t)]&=\bm{\mu}_{\varOmega}(t)\\ E[\bm{\Delta}(t)]&=\bm{\mu}_{\Delta}(t) \end{aligned}\right\} \tag{4-1-33}\] \[\operatorname{cov}[\bm{\varOmega}(t),\bm{\Delta}(\tau)]=\bm{D}_{\varOmega\Delta}(t)\,\delta(t-\tau) \tag{4-1-34}\] 即假定动态噪声和观测噪声的均值不为零，且它们在同一时刻是相关的，称为白噪声情况下的一般连续线性系统。

如果动态噪声或观测噪声是有色噪声，即不同时刻的动态噪声或观测噪声是相关的，则称这种系统为有色噪声情况下的连续线性系统。

可以证明，当线性系统 (4-1-3) 具有 (4-1-28)   (4-1-32) 诸式的随机模型时，对于一切 \(\tau>t\geq t_{0}\)，有 \[\left.\begin{aligned} \operatorname{cov}[\bm{X}(t),\bm{\varOmega}(\tau)]&=\bm{0}\\ \operatorname{cov}[\bm{L}(t),\bm{\varOmega}(\tau)]&=\bm{0} \end{aligned}\right\} \tag{4-1-35}\] 当 \(\bm{\Delta}(t)\) 与 \(\bm{\varOmega}(t)\) 完全不相关时，对于一切 \(t,\tau\geq t_{0}\)，有 \[\operatorname{cov}[\bm{X}(t),\bm{\Delta}(\tau)]=\bm{0} \tag{4-1-36}\] 对于一切 \(\tau>t\geq t_{0}\)，有 \[\operatorname{cov}[\bm{L}(t),\bm{\Delta}(\tau)]=\bm{0} \tag{4-1-37}\]

下面讨论状态向量 \(\bm{X}(t)\,(t>t_{0})\) 的均值和自协方差函数。

对于由 (4-1-4) 式表示的线性系统，将状态方程之解 (4-1-22) 的两边取数学期望，并考虑数学期望的运算和积分运算可以交换，则状态向量的均值函数满足以下方程 \[\bm{\mu}_{X}(t)=\bm{\varPhi}(t,t_{0})\bm{\mu}_{X}(t_{0}) +\int_{t_{0}}^{t}\bm{\varPhi}(t,\tau)\bm{F}(\tau)\bm{\mu}_{\varOmega}(\tau)\,\mathrm{d}\tau \tag{4-1-38}\] 上式就是在 (4-1-22) 式中用均值函数代替状态向量的结果。

状态向量的自协方差函数由下式确定 \[\bm{D}_{X}(t_{1},t_{2})=\operatorname{cov}[\bm{X}(t_{1}),\bm{X}(t_{2})] =E\left\{[\bm{X}(t_{1})-\bm{\mu}_{X}(t_{1})]\,[\bm{X}(t_{2})-\bm{\mu}_{X}(t_{2})]^{T}\right\} \tag{4-1-39}\] 因为对于一切 \(\tau\geq t_{0}\)，\(\bm{\varOmega}(\tau)\) 与 \(\bm{X}(t_{0})\) 不相关，且 \[\bm{X}(t)-\bm{\mu}_{X}(t)=\bm{\varPhi}(t,t_{0})\left[\bm{X}(t_{0})-\bm{\mu}_{X}(t_{0})\right] +\int_{t_{0}}^{t}\bm{\varPhi}(t,\tau)\bm{F}(\tau)\left[\bm{\varOmega}(\tau)-\bm{\mu}_{\varOmega}(\tau)\right]\mathrm{d}\tau\] 所以 \[\begin{aligned} \bm{D}_{X}(t_{1},t_{2})={}&\bm{\varPhi}(t_{1},t_{0})\,\bm{D}_{X}(t_{0})\,\bm{\varPhi}^{T}(t_{2},t_{0})+{}\\ &\int_{t_{0}}^{t_{1}}\!\!\int_{t_{0}}^{t_{2}}\bm{\varPhi}(t_{1},\tau_{1})\,\bm{F}(\tau_{1})\,\bm{D}_{\varOmega}(\tau_{1},\tau_{2})\, \bm{F}^{T}(\tau_{2})\,\bm{\varPhi}^{T}(t_{2},\tau_{2})\,\mathrm{d}\tau_{1}\mathrm{d}\tau_{2} \end{aligned} \tag{4-1-40}\] 当 \(\bm{\varOmega}(t)\) 是白噪声过程时，有 \[\begin{aligned} \bm{D}_{X}(t_{1},t_{2})={}&\bm{\varPhi}(t_{1},t_{0})\,\bm{D}_{X}(t_{0})\,\bm{\varPhi}^{T}(t_{2},t_{0})+{}\\ &\int_{t_{0}}^{t^{\prime}}\bm{\varPhi}(t_{1},\tau)\,\bm{F}(\tau)\,\bm{D}_{\varOmega}(\tau)\,\bm{F}^{T}(\tau)\,\bm{\varPhi}^{T}(t_{1},\tau)\,\mathrm{d}\tau \end{aligned} \tag{4-1-41}\] 其中 \(t^{\prime}=\min(t_{1},t_{2})\)。而当 \(t_{1}=t_{2}=t\) 时，有 \[\begin{aligned} \bm{D}_{X}(t)={}&\bm{\varPhi}(t,t_{0})\,\bm{D}_{X}(t_{0})\,\bm{\varPhi}^{T}(t,t_{0})+{}\\ &\int_{t_{0}}^{t}\bm{\varPhi}(t,\tau)\,\bm{F}(\tau)\,\bm{D}_{\varOmega}(\tau)\,\bm{F}^{T}(\tau)\,\bm{\varPhi}^{T}(t,\tau)\,\mathrm{d}\tau \end{aligned} \tag{4-1-42}\] 此时，(4-1-41) 式可改写为 \[\bm{D}_{X}(t_{1},t_{2})= \begin{cases} \bm{\varPhi}(t_{1},t_{2})\,\bm{D}_{X}(t_{2}), & t_{1}>t_{2}\\ \bm{D}_{X}(t_{1})\,\bm{\varPhi}^{T}(t_{2},t_{1}), & t_{1}<t_{2} \end{cases} \tag{4-1-43}\]

离散线性系统的数学模型

离散线性系统的状态方程和观测方程

如果仅在确定的瞬间 \(t_{i}\,(i=0,1,2,\cdots)\) 来研究系统的性能，则把这样的系统叫做离散时间系统。一般来说，它包括两种情况：

(1) 系统本身就是一个离散系统，对于这种系统，只能在离散的瞬间来研究系统的状态，而在这些时刻之间，系统的状态是没有意义的。

(2) 本身是连续系统，但为了研究方便，仅在离散时间内研究其性能。

一个离散线性系统，通常是用一个具有随机初始状态，并带有动态噪声的线性差分方程和离散观测方程来描述。对于上述第二种情况，可以由相应的连续线性系统的状态方程和观测方程离散化得到。下面在离散时刻 \(t_{0},t_{1},\cdots,t_{k},\cdots\) 考虑线性系统 (4-1-3)，即 \[\dot{\bm{X}}(t)=\bm{A}(t)\bm{X}(t)+\bm{C}(t)\bm{U}(t)+\bm{F}(t)\bm{\varOmega}(t) \tag{4-2-1}\] \[\bm{L}(t)=\bm{B}(t)\bm{X}(t)+\bm{G}(t)\bm{U}(t)+\bm{\Delta}(t) \tag{4-2-2}\] 由状态方程之解 (4-1-21) 可得 \[\begin{aligned} \bm{X}(t_{k+1})={}&\bm{\varPhi}(t_{k+1},t_{k})\,\bm{X}(t_{k})+{}\\ &\int_{t_{k}}^{t_{k+1}}\bm{\varPhi}(t_{k+1},\tau)\left[\bm{C}(\tau)\bm{U}(\tau)+\bm{F}(\tau)\bm{\varOmega}(\tau)\right]\mathrm{d}\tau \end{aligned} \tag{4-2-3}\] 用下标 \(k+1\) 和 \(k\) 表示 \(t_{k+1}\) 和 \(t_{k}\)，并设 \[\left.\begin{aligned} \bm{\varPsi}_{k+1,k}&=\int_{t_{k}}^{t_{k+1}}\bm{\varPhi}(t_{k+1},\tau)\,\bm{C}(\tau)\,\mathrm{d}\tau\\ \bm{\varGamma}_{k+1,k}&=\int_{t_{k}}^{t_{k+1}}\bm{\varPhi}(t_{k+1},\tau)\,\bm{F}(\tau)\,\mathrm{d}\tau\\ \bm{\varOmega}_{k}&=\frac{1}{\Delta t}\int_{t_{k}}^{t_{k+1}}\bm{\varOmega}(\tau)\,\mathrm{d}\tau \end{aligned}\right\} \tag{4-2-4}\] 对于观测方程，设 \[\bm{\Delta}_{k+1}=\frac{1}{\Delta t}\int_{t_{k}}^{t_{k+1}}\bm{\Delta}(\tau)\,\mathrm{d}\tau \tag{4-2-5}\] 其中 \(\Delta t=t_{k+1}-t_{k}\)。则 (4-2-3) 式和观测方程 (4-2-2) 可写为 \[\bm{X}_{k+1}=\bm{\varPhi}_{k+1,k}\bm{X}_{k}+\bm{\varPsi}_{k+1,k}\bm{U}_{k}+\bm{\varGamma}_{k+1,k}\bm{\varOmega}_{k} \tag{4-2-6}\] \[\bm{L}_{k+1}=\bm{B}_{k+1}\bm{X}_{k+1}+\bm{G}_{k+1}\bm{U}_{k+1}+\bm{\Delta}_{k+1} \tag{4-2-7}\] 以上两式就是离散线性系统的状态方程和观测方程。

当不考虑系统具有确定性输入时，状态方程和观测方程即为 \[\left.\begin{aligned} \bm{X}_{k+1}&=\bm{\varPhi}_{k+1,k}\bm{X}_{k}+\bm{\varGamma}_{k+1,k}\bm{\varOmega}_{k}\\ \bm{L}_{k+1}&=\bm{B}_{k+1}\bm{X}_{k+1}+\bm{\Delta}_{k+1} \end{aligned}\right\} \tag{4-2-8}\] 对于常系数线性系统，则有 \[\left.\begin{aligned} \bm{X}_{k+1}&=\bm{\varPhi}(\Delta t)\,\bm{X}_{k}+\bm{\varPsi}(\Delta t)\,\bm{U}_{k}+\bm{\varGamma}(\Delta t)\,\bm{\varOmega}_{k}\\ \bm{L}_{k+1}&=\bm{B}\bm{X}_{k+1}+\bm{G}\bm{U}_{k+1}+\bm{\Delta}_{k+1} \end{aligned}\right\} \tag{4-2-9}\] \[\left.\begin{aligned} \bm{X}_{k+1}&=\bm{\varPhi}(\Delta t)\,\bm{X}_{k}+\bm{\varGamma}(\Delta t)\,\bm{\varOmega}_{k}\\ \bm{L}_{k+1}&=\bm{B}\bm{X}_{k+1}+\bm{\Delta}_{k+1} \end{aligned}\right\} \tag{4-2-10}\]

例 4-2-1

写出 [例 4-1-1] 的状态方程和观测方程的离散化形式。

解在 [例 4-1-1] 中已求得状态转移矩阵为 \[\bm{\varPhi}(t,\tau)=\begin{bmatrix}1&t-\tau\\ 0&1\end{bmatrix}\] 所以 \[\bm{\varPhi}_{k+1,k}=\begin{bmatrix}1&t_{k+1}-t_{k}\\ 0&1\end{bmatrix}\] 而由 (4-2-4) 知 \[\begin{aligned} \bm{\varGamma}_{k+1,k}&=\int_{t_{k}}^{t_{k+1}}\bm{\varPhi}(t_{k+1},\tau)\,\bm{F}(\tau)\,\mathrm{d}\tau =\int_{t_{k}}^{t_{k+1}}\begin{bmatrix}1&t_{k+1}-\tau\\ 0&1\end{bmatrix}\begin{bmatrix}0\\ 1\end{bmatrix}\mathrm{d}\tau\\ &=\int_{t_{k}}^{t_{k+1}}\begin{bmatrix}t_{k+1}-\tau\\ 1\end{bmatrix}\mathrm{d}\tau =\begin{bmatrix}\dfrac{1}{2}(t_{k+1}-t_{k})^{2}\\[8pt] t_{k+1}-t_{k}\end{bmatrix} \end{aligned}\] 所以，离散化的状态方程为 \[\bm{X}_{k+1}=\begin{bmatrix}1&t_{k+1}-t_{k}\\ 0&1\end{bmatrix}\bm{X}_{k} +\begin{bmatrix}\dfrac{1}{2}(t_{k+1}-t_{k})^{2}\\[8pt] t_{k+1}-t_{k}\end{bmatrix}\varOmega_{k}\] 离散化的观测方程为 \[L_{k+1}=\begin{bmatrix}1&0\end{bmatrix}\bm{X}_{k+1}+\varDelta_{k+1}\]

离散线性系统的随机模型

完全不相关的白噪声情况下的离散线性系统具有以下随机模型：

(1) 动态噪声和观测噪声是零均值白噪声或高斯白噪声序列，即有 \[\left.\begin{aligned} E(\bm{\varOmega}_{k})&=\bm{0}\\ E(\bm{\Delta}_{k})&=\bm{0} \end{aligned}\right\} \tag{4-2-11}\] \[\left.\begin{aligned} \operatorname{cov}(\bm{\varOmega}_{k},\bm{\varOmega}_{j})&=\bm{D}_{\varOmega}(k)\,\delta_{kj}\\ \operatorname{cov}(\bm{\Delta}_{k},\bm{\Delta}_{j})&=\bm{D}_{\Delta}(k)\,\delta_{kj} \end{aligned}\right\} \tag{4-2-12}\] 其中 \(\bm{D}_{\varOmega}(k)\) 是 \(\bm{\varOmega}_{k}\) 的非负定方差阵，\(\bm{D}_{\Delta}(k)\) 是 \(\bm{\Delta}_{k}\) 的正定方差阵，\(\delta_{kj}\) 是克罗内克 \(\delta\) 函数。

(2) 动态噪声与观测噪声完全不相关，即 \[\operatorname{cov}(\bm{\varOmega}_{k},\bm{\Delta}_{j})=\bm{0} \tag{4-2-13}\]

(3) 系统的初始状态 \(\bm{X}_{0}\) 是具有正态分布或其他分布的随机向量，其均值和方差阵为 \[\left.\begin{aligned} E(\bm{X}_{0})&=\bm{\mu}_{X}(0)\\ \operatorname{var}(\bm{X}_{0})&=\bm{D}_{X}(0) \end{aligned}\right\} \tag{4-2-14}\] 以上的随机模型也可以说是相应的连续系统随机模型的离散化形式。

对于白噪声情况下一般离散性系统，其随机模型仅将上面的 (4-2-11) 式和 (4-2-13) 式修改为 \[\left.\begin{aligned} E(\bm{\varOmega}_{k})&=\bm{\mu}_{\varOmega}(k)\\ E(\bm{\Delta}_{k})&=\bm{\mu}_{\Delta}(k) \end{aligned}\right\} \tag{4-2-15}\] \[\operatorname{cov}(\bm{\varOmega}_{k},\bm{\Delta}_{j})=\bm{D}_{\varOmega\Delta}(k)\,\delta_{kj} \tag{4-2-16}\] 即假定动态噪声和观测噪声的均值不为零，且它们在同一时刻是相关的。

此外，在有色噪声情况下的离散线性系统中，动态噪声和观测噪声在不同时刻是相关的。

容易证明，对于常见的白噪声情况下的离散线性系统，即具有 (4-2-11)   (4-2-14) 诸式的随机模型的离散线性系统，还具有性质：对于一切 \(j\geq k\geq 0\)，有 \[\operatorname{cov}(\bm{X}_{k},\bm{\varOmega}_{j})=\bm{0} \tag{4-2-17}\] 对于一切 \(j\geq k\geq 0\)，若 \(\operatorname{cov}(\bm{\varOmega}_{k},\bm{\Delta}_{j})=\bm{0}\)，有 \[\operatorname{cov}(\bm{L}_{k},\bm{\varOmega}_{j})=\bm{0} \tag{4-2-18}\] 对于一切 \(j>k\geq 0\)，有 \[\operatorname{cov}(\bm{L}_{k},\bm{\Delta}_{j})=\bm{0} \tag{4-2-19}\]

下面讨论状态向量 \(\bm{X}_{k}\,(k>0)\) 的均值和自协方差函数。由状态方程 (4-2-6) 可写出 \[\left.\begin{aligned} \bm{X}_{k}&=\bm{\varPhi}_{k,k-1}\bm{X}_{k-1}+\bm{\varPsi}_{k,k-1}\bm{U}_{k-1}+\bm{\varGamma}_{k,k-1}\bm{\varOmega}_{k-1}\\ \bm{X}_{k-1}&=\bm{\varPhi}_{k-1,k-2}\bm{X}_{k-2}+\bm{\varPsi}_{k-1,k-2}\bm{U}_{k-2}+\bm{\varGamma}_{k-1,k-2}\bm{\varOmega}_{k-2}\\ &\cdots \end{aligned}\right\}\] 且由状态转移矩阵的性质知 \[\bm{\varPhi}_{k,j}=\bm{\varPhi}_{k,i}\bm{\varPhi}_{i,j}\] 所以 \[\bm{X}_{k}=\bm{\varPhi}_{k,j}\bm{X}_{j} +\sum_{i=j+1}^{k}\bm{\varPhi}_{k,i}\bm{\varPsi}_{i,i-1}\bm{U}_{i-1} +\sum_{i=j+1}^{k}\bm{\varPhi}_{k,i}\bm{\varGamma}_{i,i-1}\bm{\varOmega}_{i-1} \tag{4-2-20}\] 若取 \(j=0\)，则 \[\bm{X}_{k}=\bm{\varPhi}_{k,0}\bm{X}_{0} +\sum_{i=1}^{k}\bm{\varPhi}_{k,i}\bm{\varPsi}_{i,i-1}\bm{U}_{i-1} +\sum_{i=1}^{k}\bm{\varPhi}_{k,i}\bm{\varGamma}_{i,i-1}\bm{\varOmega}_{i-1} \tag{4-2-21}\] 由 (4-2-6) 式和 (4-2-21) 式可得状态向量 \(\bm{X}_{k+1}\) 的均值函数为 \[\bm{\mu}_{X}(k+1)=\bm{\varPhi}_{k+1,k}\bm{\mu}_{X}(k) +\bm{\varPsi}_{k+1,k}\bm{U}_{k}+\bm{\varGamma}_{k+1,k}\bm{\mu}_{\varOmega}(k) \tag{4-2-22}\] \[\bm{\mu}_{X}(k+1)=\bm{\varPhi}_{k+1,0}\bm{\mu}_{X}(0) +\sum_{i=1}^{k+1}\bm{\varPhi}_{k+1,i}\bm{\varPsi}_{i,i-1}\bm{U}_{i-1} +\sum_{i=1}^{k+1}\bm{\varPhi}_{k+1,i}\bm{\varGamma}_{i,i-1}\bm{\mu}_{\varOmega}(i-1) \tag{4-2-23}\] 当 \(\bm{\mu}_{\varOmega}(k)=\bm{0}\) 时，则上式最后一项为零。

对于白噪声情况下的离散线性系统，因为 \(\bm{U}_{k}\,(k\geq 0)\) 是非随机项，故由 (4-2-21) 式按协方差传播律可得状态向量的自协方差函数为 \[\begin{aligned} \bm{D}_{X}(k,j)&=\operatorname{cov}(\bm{X}_{k},\bm{X}_{j}) =\bm{\varPhi}_{k,0}\bm{D}_{X}(0)\,\bm{\varPhi}_{j,0}^{T}+{}\\ &\qquad\sum_{i=1}^{k^{\prime}}\bm{\varPhi}_{k,i}\bm{\varGamma}_{i,i-1}\bm{D}_{\varOmega}(i-1)\,\bm{\varGamma}_{i,i-1}^{T}\bm{\varPhi}_{j,i}^{T} \end{aligned} \tag{4-2-24}\] 其中 \(k^{\prime}=\min(k,j)\)，此时，\(\bm{X}_{k}\) 的方差阵为 \[\bm{D}_{X}(k)=\bm{\varPhi}_{k,0}\bm{D}_{X}(0)\,\bm{\varPhi}_{k,0}^{T} +\sum_{i=1}^{k}\bm{\varPhi}_{k,i}\bm{\varGamma}_{i,i-1}\bm{D}_{\varOmega}(i-1)\,\bm{\varGamma}_{i,i-1}^{T}\bm{\varPhi}_{k,i}^{T} \tag{4-2-25}\] 显然，协方差函数可表示为 \[\bm{D}_{X}(k,j)= \begin{cases} \bm{\varPhi}_{k,j}\bm{D}_{X}(j), & k>j\\ \bm{D}_{X}(k)\,\bm{\varPhi}^{T}(j,k), & k<j \end{cases} \tag{4-2-26}\] 此外，由 (4-2-6) 式按协方差传播律得 \[\bm{D}_{X}(k+1)=\bm{\varPhi}_{k+1,k}\bm{D}_{X}(k)\,\bm{\varPhi}_{k+1,k}^{T} +\bm{\varGamma}_{k+1,k}\bm{D}_{\varOmega}(k)\,\bm{\varGamma}_{k+1,k}^{T} \tag{4-2-27}\]

例 4-2-2

设有离散状态方程 \[X_{k+1}=(k+1)X_{k}+\varOmega_{k}\] 动态噪声是零均值白噪声序列，且与 \(X_{0}\) 独立，而 \(D_{\varOmega}(k)=k^{2}\)，\(\mu_{x}(0)=1\)，\(D_{X}(0)=1\)，试求 \(\mu_{X}(k)\) 和 \(D_{X}(k)\)。

解因 \(\varPhi_{k+1,k}=k+1\)，\(\mu_{\varOmega}(k)=0\)，\(\varGamma_{k+1,k}=1\)，所以 \[\begin{gathered} \varPhi_{k,0}=\varPhi_{k,k-1}\cdots\varPhi_{2,1}\varPhi_{1,0}=k!\\ \varPhi_{k,i}=\varPhi_{k,k-1}\cdots\varPhi_{i+2,i+1}\varPhi_{i+1,i}=\frac{k!}{i!} \end{gathered}\] 故由 (4-2-22)、(4-2-23) 式可得 \[\mu_{X}(k)=\varPhi_{k,0}\mu_{X}(0)=k!\] 而由 (4-2-25) 式可得 \[\begin{aligned} D_{X}(k)&=\varPhi_{k,0}D_{X}(0)\,\varPhi_{k,0}^{T} +\sum_{i=1}^{k}\varPhi_{k,i}\varGamma_{i,i-1}D_{\varOmega}(i-1)\,\varGamma_{i,i-1}^{T}\varPhi_{k,i}^{T}\\ &=(k!)^{2}+\sum_{i=1}^{k}\left(\frac{k!}{i!}\right)^{2}i^{2} \end{aligned}\]

离散线性系统的卡尔曼滤波

本节讨论由 (4-2-6)、(4-2-7) 式所表示的线性系统的状态估计问题，通常称为离散线性系统的卡尔曼（Kalman）滤波。离散线性系统的状态方程和观测方程为 \[\bm{X}_{k+1}=\bm{\varPhi}_{k+1,k}\bm{X}_{k}+\bm{\varPsi}_{k+1,k}\bm{U}_{k}+\bm{\varGamma}_{k+1,k}\bm{\varOmega}_{k} \tag{4-3-1}\] \[\bm{L}_{k+1}=\bm{B}_{k+1}\bm{X}_{k+1}+\bm{Z}_{k+1}+\bm{\Delta}_{k+1} \tag{4-3-2}\] 式中，\(\bm{Z}_{k+1}=\bm{G}_{k+1}\bm{U}_{k+1}\)，其他符号的意义与 4-2 节相同。

所谓离散线性系统的状态估计，就是利用观测向量 \(\bm{L}_{1},\bm{L}_{2},\cdots,\bm{L}_{k}\)，根据其数学模型求定第 \(t_{j}\) 时刻状态向量 \(\bm{X}_{j}\) 的最佳估值。通常把所得到的估计量记为 \(\hat{\bm{X}}(j/k)\)。它可以分为三种情况：

(1) 当 \(j=k\) 时，称 \(\hat{\bm{X}}(k/k)\) 为最佳滤波值，并把 \(\hat{\bm{X}}(k/k)\) 的求定过程称为卡尔曼滤波；

(2) 当 \(j>k\) 时，称 \(\hat{\bm{X}}(j/k)\) 为最佳预测值，而把求定 \(\hat{\bm{X}}(j/k)\) 的过程称为预测或外推；

(3) 当 \(j<k\) 时，称 \(\hat{\bm{X}}(j/k)\) 为最佳平滑值，而把求定 \(\hat{\bm{X}}(j/k)\) 的过程称为平滑或内插。

设系统 (4-3-1)、(4-3-2) 的随机模型就是 4-2 节中给出的最基本的模型，即有 \[\left.\begin{aligned} E(\bm{\varOmega}_{k})&=\bm{0}\\ E(\bm{\Delta}_{k})&=\bm{0} \end{aligned}\right\} \tag{4-3-3}\] \[\left.\begin{aligned} \operatorname{cov}(\bm{\varOmega}_{k},\bm{\varOmega}_{j})&=\bm{D}_{\varOmega}(k)\,\delta_{kj}\\ \operatorname{cov}(\bm{\Delta}_{k},\bm{\Delta}_{j})&=\bm{D}_{\Delta}(k)\,\delta_{kj}\\ \operatorname{cov}(\bm{\varOmega}_{k},\bm{\Delta}_{j})&=\bm{0} \end{aligned}\right\} \tag{4-3-4}\] \[\left.\begin{aligned} E(\bm{X}_{0})&=\bm{\mu}_{X}(0)=\hat{\bm{X}}(0/0)\\ \operatorname{var}(\bm{X}_{0})&=\bm{D}_{X}(0)\\ \operatorname{cov}(\bm{X}_{0},\bm{\varOmega}_{k})&=\bm{0}\\ \operatorname{cov}(\bm{X}_{0},\bm{\Delta}_{k})&=\bm{0} \end{aligned}\right\} \tag{4-3-5}\] 其中动态噪声和观测噪声都是零均值白噪声序列，且它们之间在任何时刻都不相关，因而将对上述模型的卡尔曼滤波，称为完全不相关的白噪声作用下的卡尔曼滤波。

本节仅讨论求 \(\hat{\bm{X}}(k/k)\) 的卡尔曼滤波问题，也就是讨论在已知 \(\bm{X}_{0}\)、\(\bm{\varOmega}_{k}\) 和 \(\bm{\Delta}_{k}\) 的前二阶矩的情况下，对系统的状态向量 \(\bm{X}_{k}\) 进行估计求定 \(\hat{\bm{X}}(k/k)\) 的问题。应该说，这种估计问题可以直接应用第 1 章所介绍的估计方法的基本公式来解决。但是，容易发现，直接应用那些公式在计算上有困难。为此，需要导出滤波的递推公式，应用这些递推公式，既可以减少计算工作量，提高计算速度；又不需要储存大量以前的观测数据，从而满足了实用上的要求。

对于离散线性系统的状态估计问题，可以应用第 1 章所述的最小方差估计或极大验后估计原理来导出求估值 \(\hat{\bm{X}}(k/k)\) 的卡尔曼滤波方程，也可以应用广义最小二乘原理来推导。采用这些方法所得到的结果是相同的。本节阐述应用广义最小二乘原理并采用逐次平差方法导出卡尔曼滤波的过程。

卡尔曼滤波方程

当有观测向量 \(\bm{L}_{1},\bm{L}_{2},\cdots,\bm{L}_{k}\) 时，状态方程和观测方程可写为如下形式： \[\left.\begin{aligned} \bm{X}_{1}&=\bm{\varPhi}_{1,0}\bm{X}_{0}+\bm{\varPsi}_{1,0}\bm{U}_{0}+\bm{\varGamma}_{1,0}\bm{\varOmega}_{0}\\ \bm{L}_{1}&=\bm{B}_{1}\bm{X}_{1}+\bm{Z}_{1}+\bm{\Delta}_{1}\\ \bm{X}_{2}&=\bm{\varPhi}_{2,1}\bm{X}_{1}+\bm{\varPsi}_{2,1}\bm{U}_{1}+\bm{\varGamma}_{2,1}\bm{\varOmega}_{1}\\ \bm{L}_{2}&=\bm{B}_{2}\bm{X}_{2}+\bm{Z}_{2}+\bm{\Delta}_{2}\\ &\cdots\\ \bm{X}_{k}&=\bm{\varPhi}_{k,k-1}\bm{X}_{k-1}+\bm{\varPsi}_{k,k-1}\bm{U}_{k-1}+\bm{\varGamma}_{k,k-1}\bm{Q}_{k-1}\\ \bm{L}_{k}&=\bm{B}_{k}\bm{X}_{k}+\bm{Z}_{k}+\bm{\Delta}_{k} \end{aligned}\right\} \tag{4-3-6}\] 可以看到，随机模型诸式给出的仅仅是初始状态 \(\bm{X}_{0}\)，以及动态噪声 \(\bm{\varOmega}_{k}\) 和观测噪声 \(\bm{\Delta}_{k}\) 的统计性质，而状态 \(\bm{X}_{1},\bm{X}_{2},\cdots\) 的统计性质实际上是未知的。因此，可以说它们都是非随机参数，只有初始状态是随机参数。同时，如果将状态方程 (4-3-1) 中包含动态噪声的项 \(\bm{\varGamma}_{k+1,k}\bm{\varOmega}_{k}\,(k\geq 0)\) 看成是对应于某种观测向量 \(\bm{L}_{s}(k)\) 的观测噪声，则状态方程也可以看成是相应的观测方程，可令 \[\left.\begin{aligned} \bm{\varepsilon}_{k}&=\bm{\varGamma}_{k+1,k}\bm{\varOmega}_{k}\\ \bm{L}_{s}(k)&=-\bm{\varPsi}_{k+1,k}\bm{U}_{k} \end{aligned}\right\} \tag{4-3-7}\] 此时，\(\bm{\varepsilon}_{k}\) 的方差阵为 \[\bm{D}_{\varepsilon}(k)=\bm{\varGamma}_{k+1,k}\bm{D}_{\varOmega}(k)\,\bm{\varGamma}_{k+1,k}^{T} \tag{4-3-8}\] 于是，(4-3-1) 式可写为 \[\bm{L}_{s}(k)=\bm{\varPhi}_{k+1,k}\bm{X}_{k}-\bm{X}_{k+1}+\bm{\varepsilon}_{k} \tag{4-3-9}\] 根据广义最小二乘原理，应将随机参数 \(\bm{X}_{0}\) 的先验期望看成是虚拟观测值。则可以由 (4-3-6) 式和 (4-3-9) 式写出误差方程 \[\left.\begin{aligned} \bm{V}_{X}(0)&=\hat{\bm{X}}_{0}-\hat{\bm{X}}(0/0)\\ \bm{V}_{s}(0)&=\bm{\varPhi}_{1,0}\hat{\bm{X}}_{0}-\hat{\bm{X}}_{1}-\bm{L}_{s}(0)\\ \bm{V}_{1}&=\bm{B}_{1}\hat{\bm{X}}_{1}+\bm{Z}_{1}-\bm{L}_{1}\\ \bm{V}_{s}(1)&=\bm{\varPhi}_{2,1}\hat{\bm{X}}_{1}-\hat{\bm{X}}_{2}-\bm{L}_{s}(1)\\ \bm{V}_{2}&=\bm{B}_{2}\hat{\bm{X}}_{2}+\bm{Z}_{2}-\bm{L}_{2}\\ &\cdots\\ \bm{V}_{s}(k-1)&=\bm{\varPhi}_{k,k-1}\hat{\bm{X}}_{k-1}-\hat{\bm{X}}_{k}-\bm{L}_{s}(k-1)\\ \bm{V}_{k}&=\bm{B}_{k}\hat{\bm{X}}_{k}+\bm{Z}_{k}-\bm{L}_{k} \end{aligned}\right\} \tag{4-3-10}\] 这就是说，根据广义最小二乘原理，可以将上述动态系统的状态估计问题变换为一个最小二乘平差问题。此时 (4-3-10) 式中的 \(\hat{\bm{X}}_{0},\hat{\bm{X}}_{1},\cdots,\hat{\bm{X}}_{k}\) 均看做非随机参数，而 \(\hat{\bm{X}}(0/0)\)，\(\bm{L}_{s}(0)\)，\(\cdots\)，\(\bm{L}_{s}(k-1)\)，\(\bm{L}_{1}\)，\(\cdots\)，\(\bm{L}_{k}\) 是方差阵分别为 \(\bm{D}_{X}(0)\)，\(\bm{D}_{\varepsilon}(0)\)，\(\cdots\)，\(\bm{D}_{\varepsilon}(k-1)\)，\(\bm{D}_{\Delta}(1)\)，\(\cdots\)，\(\bm{D}_{\Delta}(k)\) 的观测值，它们之间的协方差阵均为零。

由于 (4-3-10) 式中所包含的观测值很多，将它们放在一起进行整体平差很不方便。为此，可以采用逐次平差法进行计算，或者是按逐次平差法导出一组递推计算的公式来计算，这组公式称为卡尔曼滤波方程，其推导步骤如下：

(1) 取 (4-3-10) 式中的第一、二式进行平差，即利用 \(\hat{\bm{X}}(0/0)\) 和 \(\bm{L}_{s}(0)\) 进行平差，此时，相应的法方程为 \[\left.\begin{aligned} &\left[\bm{D}_{X}^{-1}(0)+\bm{\varPhi}_{1,0}^{T}\bm{D}_{\varepsilon}^{-1}(0)\bm{\varPhi}_{1,0}\right]\hat{\bm{X}}_{s}(0/0) -\bm{\varPhi}_{1,0}^{T}\bm{D}_{\varepsilon}^{-1}(0)\hat{\bm{X}}(1/0)-{}\\ &\qquad\left[\bm{D}_{X}^{-1}(0)\hat{\bm{X}}(0/0)+\bm{\varPhi}_{1,0}^{T}\bm{D}_{\varepsilon}^{-1}(0)\bm{L}_{s}(0)\right]=0\\ &-\bm{D}_{\varepsilon}^{-1}(0)\bm{\varPhi}_{1,0}\hat{\bm{X}}_{s}(0/0) +\bm{D}_{\varepsilon}^{-1}(0)\hat{\bm{X}}(1/0)+\bm{D}_{\varepsilon}^{-1}(0)\bm{L}_{s}(0)=0 \end{aligned}\right\} \tag{4-3-11}\] 式中 \(\hat{\bm{X}}_{s}(0/0)\) 和 \(\hat{\bm{X}}(1/0)\) 表示由 \(\hat{\bm{X}}(0/0)\) 和 \(\bm{L}_{s}(0)\) 平差求得的 \(\bm{X}_{0}\) 和 \(\bm{X}_{1}\) 的估值。将上式的第二式左乘 \(\bm{\varPhi}_{1,0}^{T}\)，再与第一式相加，即可得 \[\hat{\bm{X}}_{s}(0/0)=\hat{\bm{X}}(0/0)\] 再将它代入 (4-3-11) 的第二式，又可得 \[\hat{\bm{X}}(1/0)=\bm{\varPhi}_{1,0}\hat{\bm{X}}(0/0)-\bm{L}_{s}(0)\] 或写为 \[\hat{\bm{X}}(1/0)=\bm{\varPhi}_{1,0}\hat{\bm{X}}(0/0)+\bm{\varPsi}_{1,0}\bm{U}_{0} \tag{4-3-12}\] 也就是说，利用 \(\hat{\bm{X}}(0/0)\) 和 \(\bm{L}_{s}(0)\) 进行平差时，所得到的平差值 \(\hat{\bm{X}}_{s}(0/0)\) 就等于初始状态 \(\bm{X}_{0}\) 的先验期望 \(\hat{\bm{X}}(0/0)\)；而 \(\hat{\bm{X}}(1/0)\) 等于以 \(\hat{\bm{X}}(0/0)\) 和 \(\bm{\varOmega}_{0}\) 的数学期望 \(E(\bm{\varOmega}_{0})=\bm{0}\) 代入状态方程所得到的结果。实际上，(4-3-10) 式中的第一、二式中包含有 \(2m\) 个参数（\(\bm{X}_{0}\)，\(\bm{X}_{1}\)），而观测值 \(\hat{\bm{X}}(0/0)\) 和 \(\bm{L}_{s}(0)\) 中也只有 \(2m\) 个，就是说，它们没有多余观测，自然会得到这样的结果。在求定 \(\hat{\bm{X}}_{s}(0/0)\) 和 \(\hat{\bm{X}}(1/0)\) 时，不需要组成法方程 (4-3-11)，只需要由 (4-3-12) 式求得 \(\hat{\bm{X}}(1/0)\)，因此，在随机模型中，只要求 \(\bm{D}_{X}(0)\) 和 \(\bm{D}_{\varOmega}(k)\) 是非负定阵。

由 (4-3-12) 式按协方差传播律可得到 \(\hat{\bm{X}}(1/0)\) 的误差方差阵及它与 \(\hat{\bm{X}}(0/0)\) 的协方差阵为 \[\left.\begin{aligned} \bm{D}_{\hat{X}}(1/0)&=\bm{\varPhi}_{1,0}\bm{D}_{X}(0)\bm{\varPhi}_{1,0}^{T} +\bm{\varGamma}_{1,0}\bm{D}_{\varOmega}(0)\bm{\varGamma}_{1,0}^{T}\\ \operatorname{cov}[\hat{\bm{X}}(1/0),\hat{\bm{X}}(0/0)]&=\bm{\varPhi}_{1,0}\bm{D}_{X}(0) \end{aligned}\right\} \tag{4-3-13}\]

(2) 按逐次平差法，以 \(\hat{\bm{X}}(1/0)\) 作为虚拟观测值，与 \(\bm{L}_{1}\) 一起作第二次平差，此时的误差方程为 \[\left.\begin{aligned} \bm{V}_{X}(1/0)&=\hat{\bm{X}}(1/1)-\hat{\bm{X}}(1/0)\\ \overline{\bm{V}}_{1}&=\bm{B}_{1}\hat{\bm{X}}(1/1)+\bm{Z}_{1}-\bm{L}_{1} \end{aligned}\right\} \tag{4-3-14}\] 其法方程为 \[\left[\bm{D}_{X}^{-1}(1/0)+\bm{B}_{1}^{T}\bm{D}_{\Delta}^{-1}(1)\bm{B}_{1}\right]\hat{\bm{X}}(1/1) -\left[\bm{D}_{X}^{-1}(1/0)\hat{\bm{X}}(1/0)+\bm{B}_{1}^{T}\bm{D}_{\Delta}^{-1}(1)\,(\bm{L}_{1}-\bm{Z}_{1})\right]=0 \tag{4-3-15}\] 可解得 \[\hat{\bm{X}}(1/1)=\left[\bm{D}_{X}^{-1}(1/0)+\bm{B}_{1}^{T}\bm{D}_{\Delta}^{-1}(1)\bm{B}_{1}\right]^{-1} \left[\bm{D}_{X}^{-1}(1/0)\hat{\bm{X}}(1/0)+\bm{B}_{1}^{T}\bm{D}_{\Delta}^{-1}(1)\,(\bm{L}_{1}-\bm{Z}_{1})\right] \tag{4-3-16}\] 其方差阵为 \[\bm{D}_{\hat{X}}(1/1)=\left[\bm{D}_{X}^{-1}(1/0)+\bm{B}_{1}^{T}\bm{D}_{\Delta}^{-1}(1)\bm{B}_{1}\right]^{-1} \tag{4-3-17}\] 利用矩阵反演公式可将以上二式变换为 \[\hat{\bm{X}}(1/1)=\hat{\bm{X}}(1/0)+\bm{D}_{\hat{X}}(1/0)\bm{B}_{1}^{T} \left[\bm{B}_{1}\bm{D}_{\hat{X}}(1/0)\bm{B}_{1}^{T}+\bm{D}_{\Delta}(1)\right]^{-1}\cdot \left[\bm{L}_{1}-\bm{Z}_{1}-\bm{B}_{1}\hat{\bm{X}}(1/0)\right] \tag{4-3-18}\] \[\bm{D}_{\hat{X}}(1/1)=\bm{D}_{\hat{X}}(1/0)+\bm{D}_{\hat{X}}(1/0)\bm{B}_{1}^{T} \left[\bm{B}_{1}\bm{D}_{\hat{X}}(1/0)\bm{B}_{1}^{T}+\bm{D}_{\Delta}(1)\right]^{-1}\bm{B}_{1}\bm{D}_{\hat{X}}(1/0) \tag{4-3-19}\]

(3) 再以 \(\hat{\bm{X}}(1/1)\) 作为虚拟观测值，与 \(\bm{L}_{s}(1)\) 一起作第三次平差，此时的误差方程可写为 \[\left.\begin{aligned} \bm{V}_{X}(1/1)&=\hat{\bm{X}}_{s}(1/1)-\hat{\bm{X}}(1/1)\\ \bm{V}_{s}(1)&=\bm{\varPhi}_{2,1}\hat{\bm{X}}_{s}(1/1)-\hat{\bm{X}}(2/1)-\bm{L}_{s}(1) \end{aligned}\right\} \tag{4-3-20}\] 将上式与第一式平差时的误差方程 (4-3-10) 式中的第一、二式比较可知，两式的形式相同，仅相应的序号增加 1，且 (4-3-20) 式中的未知参数 \(\hat{\bm{X}}_{s}(1/1)\) 和 \(\hat{\bm{X}}(2/1)\) 为 \(2m\) 个，观测值 \(\hat{\bm{X}}(1/1)\) 和 \(\bm{L}_{s}(1)\) 也是 \(2m\) 个，因此，可以类似地得到 \[\hat{\bm{X}}_{s}(1/1)=\hat{\bm{X}}(1/1)\] \[\hat{\bm{X}}(2/1)=\bm{\varPhi}_{2,1}\hat{\bm{X}}(1/1)+\bm{\varPsi}_{2,1}\bm{U}_{1} \tag{4-3-21}\] 而 \(\hat{\bm{X}}(2/1)\) 的方差阵及它与 \(\hat{\bm{X}}(1/1)\) 的协方差阵为 \[\left.\begin{aligned} \bm{D}_{\hat{X}}(2/1)&=\bm{\varPhi}_{2,1}\bm{D}_{\hat{X}}(1/1)\bm{\varPhi}_{2,1}^{T} +\bm{\varGamma}_{2,1}\bm{D}_{\varOmega}(1)\bm{\varGamma}_{2,1}^{T}\\ \operatorname{cov}[\hat{\bm{X}}(2/1),\hat{\bm{X}}(1,1)]&=\bm{\varPhi}_{2,1}\bm{D}_{\hat{X}}(1/1) \end{aligned}\right\} \tag{4-3-22}\]

(4) 又以 \(\hat{\bm{X}}(2/1)\) 作为虚拟观测值，与 \(\bm{L}_{2}\) 一起作第四次平差，此时的误差方程可写为 \[\left.\begin{aligned} \bm{V}_{X}(2/1)&=\hat{\bm{X}}(2/2)-\hat{\bm{X}}(2/1)\\ \overline{\bm{V}}_{2}&=\bm{B}_{2}\hat{\bm{X}}(2/2)+\bm{Z}_{2}-\bm{L}_{2} \end{aligned}\right\} \tag{4-3-23}\] 上式与 (4-3-14) 式的形式相同，仅相应的序号增加 1，因此，可类似解得 \[\hat{\bm{X}}(2/2)=\left[\bm{D}_{X}^{-1}(2/1)+\bm{B}_{2}^{T}\bm{D}_{\Delta}^{-1}(2)\bm{B}_{2}\right]^{-1} \left[\bm{D}_{X}^{-1}(2/1)\hat{\bm{X}}(2/1)+\bm{B}_{2}^{T}\bm{D}_{\Delta}^{-1}(2)\,(\bm{L}_{2}-\bm{Z}_{2})\right] \tag{4-3-24}\] 其方差阵为 \[\bm{D}_{\hat{X}}(2/2)=\left[\bm{D}_{X}^{-1}(2/1)+\bm{B}_{2}^{T}\bm{D}_{\Delta}^{-1}(2)\bm{B}_{2}\right]^{-1} \tag{4-3-25}\] 由矩阵反演公式也可得到 \[\resizebox{\textwidth}{!}{$ \left.\begin{aligned} \hat{\bm{X}}(2/2)&=\hat{\bm{X}}(2/1)+\bm{D}_{\hat{X}}(2/1)\bm{B}_{2}^{T} \left[\bm{B}_{2}\bm{D}_{\hat{X}}(2/1)\bm{B}_{2}^{T}+\bm{D}_{\Delta}(2)\right]^{-1} \left[\bm{L}_{2}-\bm{Z}_{2}-\bm{B}_{2}\hat{\bm{X}}(2/1)\right]\\ \bm{D}_{\hat{X}}(2/2)&=\bm{D}_{\hat{X}}(2/1)+\bm{D}_{\hat{X}}(2/1)\bm{B}_{2}^{T} \left[\bm{B}_{2}\bm{D}_{\hat{X}}(2/1)\bm{B}_{2}^{T}+\bm{D}_{\Delta}(2)\right]^{-1}\bm{B}_{2}\bm{D}_{\hat{X}}(2/1) \end{aligned}\right\}$} \tag{4-3-26}\]

(5) 依此类推，就可得到递推计算公式，即卡尔曼滤波方程。由 (4-3-21) 式、(4-3-22) 式和 (4-3-26) 式可得其形式为 \[\left.\begin{aligned} \hat{\bm{X}}(k/k)&=\hat{\bm{X}}(k/k-1) +\bm{J}_{k}\left[\bm{L}_{k}-\bm{Z}_{k}-\bm{B}_{k}\hat{\bm{X}}(k/k-1)\right]\\ \bm{D}_{\hat{X}}(k/k)&=(\bm{E}-\bm{J}_{k}\bm{B}_{k})\,\bm{D}_{\hat{X}}(k/k-1) \end{aligned}\right\} \tag{4-3-27}\] 式中 \[\left.\begin{aligned} \hat{\bm{X}}(k/k-1)&=\bm{\varPhi}_{k,k-1}\hat{\bm{X}}(k-1/k-1)+\bm{\varPsi}_{k,k-1}\bm{U}_{k-1}\\ \bm{D}_{\hat{X}}(k/k-1)&=\bm{\varPhi}_{k,k-1}\bm{D}_{\hat{X}}(k-1/k-1)\bm{\varPhi}_{k,k-1}^{T} +\bm{\varGamma}_{k,k-1}\bm{D}_{\varOmega}(k-1)\bm{\varGamma}_{k,k-1}^{T} \end{aligned}\right\} \tag{4-3-28}\] \[\bm{J}_{k}=\bm{D}_{\hat{X}}(k/k-1)\bm{B}_{k}^{T} \left[\bm{B}_{k}\bm{D}_{\hat{X}}(k/k-1)\bm{B}_{k}^{T}+\bm{D}_{\Delta}(k)\right]^{-1} \tag{4-3-29}\] (4-3-28) 式是利用 \(\hat{\bm{X}}(k-1/k-1)\) 求 \(\bm{X}_{k}\) 的预测值 \(\hat{\bm{X}}(k/k-1)\) 及其方差公式，也就是利用 \(\bm{L}_{1}\)，\(\bm{L}_{2}\)，\(\cdots\)，\(\bm{L}_{k-1}\) 求 \(\bm{X}_{k}\) 的预测值的公式，称此式为一步预测公式；而称 \(\hat{\bm{X}}(k/k-1)\) 为一步预测值。(4-3-27) 的第一式表示滤波值 \(\hat{\bm{X}}(k/k)\) 等于用 \(\bm{J}_{k}\left[\bm{L}_{k}-\bm{Z}_{k}-\bm{B}_{k}\hat{\bm{X}}(k/k-1)\right]\) 修正一步预测值 \(\hat{\bm{X}}(k/k-1)\) 后的结果。因此，通常称 \(\bm{J}_{k}\) 为滤波增益矩阵。

由 (4-3-24)、(4-3-25) 式还可得卡尔曼滤波方程的另一种形式为 \[\left.\begin{aligned} \hat{\bm{X}}(k/k)={}&\left[\bm{D}_{X}^{-1}(k/k-1)+\bm{B}_{k}^{T}\bm{D}_{\Delta}^{-1}(k)\bm{B}_{k}\right]^{-1} \left[\bm{D}_{X}^{-1}(k/k-1)\hat{\bm{X}}(k/k-1)+{}\right.\\ &\left.\bm{B}_{k}^{T}\bm{D}_{\Delta}^{-1}(k)\,(\bm{L}_{k}-\bm{Z}_{k})\right]\\ \bm{D}_{\hat{X}}(k/k)={}&\left[\bm{D}_{X}^{-1}(k/k-1)+\bm{B}_{k}^{T}\bm{D}_{\Delta}^{-1}(k)\bm{B}_{k}\right]^{-1} \end{aligned}\right\} \tag{4-3-30}\] 可以看到，这种滤波递推公式与逐次平差的计算公式的形式相同。

当动态方程和观测方程 (4-3-1)、(4-3-2) 中的非随机控制项为零，即 \(\bm{U}_{k}=\bm{0}\)，\(\bm{Z}_{k}=\bm{0}\) 时，卡尔曼滤波方程为 \[\left.\begin{aligned} \hat{\bm{X}}(k/k)&=\hat{\bm{X}}(k/k-1)+\bm{J}_{k}\left[\bm{L}_{k}-\bm{B}_{k}\hat{\bm{X}}(k/k-1)\right]\\ \bm{D}_{\hat{X}}(k/k)&=(\bm{E}-\bm{J}_{k}\bm{B}_{k})\,\bm{D}_{\hat{X}}(k/k-1) \end{aligned}\right\} \tag{4-3-31}\] 其中 \[\left.\begin{aligned} \hat{\bm{X}}(k/k-1)&=\bm{\varPhi}_{k,k-1}\hat{\bm{X}}(k-1/k-1)\\ \bm{D}_{\hat{X}}(k/k-1)&=\bm{\varPhi}_{k,k-1}\bm{D}_{\hat{X}}(k-1/k-1)\bm{\varPhi}_{k,k-1}^{T} +\bm{\varGamma}_{k,k-1}\bm{D}_{\varOmega}(k-1)\bm{\varGamma}_{k,k-1}^{T}\\ \bm{J}_{k}&=\bm{D}_{\hat{X}}(k/k-1)\bm{B}_{k}^{T} \left[\bm{B}_{k}\bm{D}_{\hat{X}}(k/k-1)\bm{B}_{k}^{T}+\bm{D}_{\Delta}(k)\right]^{-1} \end{aligned}\right\} \tag{4-3-32}\] 以及 \[\left.\begin{aligned} \hat{\bm{X}}(k/k)={}&\left[\bm{D}_{X}^{-1}(k/k-1)+\bm{B}_{k}^{T}\bm{D}_{\Delta}^{-1}(k)\bm{B}_{k}\right]^{-1} \left[\bm{D}_{X}^{-1}(k/k-1)\hat{\bm{X}}(k/k-1)+{}\right.\\ &\left.\bm{B}_{k}^{T}\bm{D}_{\Delta}^{-1}(k)\bm{L}_{k}\right]\\ \bm{D}_{\hat{X}}(k/k)={}&\left[\bm{D}_{X}^{-1}(k/k-1)+\bm{B}_{k}^{T}\bm{D}_{\Delta}^{-1}(k)\bm{B}_{k}\right]^{-1} \end{aligned}\right\} \tag{4-3-33}\]

卡尔曼滤波的特点

卡尔曼滤波方程是一组递推计算公式，其计算过程是一个不断地预测、修正的过程。在求解时不需要储存大量的观测数据，并且当得到新的观测数据时，可随时算得新的滤波值，便于实时处理观测成果。

卡尔曼滤波还具有以下特点：

(1) 由卡尔曼滤波方程 (4-3-27)   (4-3-29) 可以看出，对 \(\bm{D}_{X}(0)\)、\(\bm{D}_{\Delta}(k)\) 和 \(\bm{D}_{\varOmega}(k-1)\) 同乘一个相同的标量时，增益矩阵 \(\bm{J}_{k}\) 之值不变。因此，在实际计算中可以用相应的协因数阵代替它们。当然，在按 (4-3-30) 式进行计算时，也可以用相应的权阵代替式中的方差阵的逆阵。

(2) 容易看到，增益矩阵 \(\bm{J}_{k}\) 与观测值无关，因此，可以预先算出，从而减少实时处理的计算工作量。

(3) 从滤波方程还可以看到，当 \(\bm{D}_{\Delta}(k)\) 增大时，增益矩阵就变小。这在直观上是容易理解的，因为如果观测噪声增大，则滤波的增益就应取小一些，以减弱观测噪声对滤波值的影响。而当 \(\bm{D}_{X}(0)\) 或 \(\bm{D}_{\varOmega}(k-1)\) 变小时，\(\bm{D}_{\hat{X}}(k/k-1)\) 变小，且 \(\bm{D}_{\hat{X}}(k/k)\) 和 \(\bm{J}_{k}\) 也变小，这在直观上也是显然的，因为 \(\bm{D}_{X}(0)\) 变小，表示初始状态较好，\(\bm{D}_{\varOmega}(k-1)\) 变小，表示系统的干扰变小，故增益矩阵也应小一些，以便给予较小的修正。这就是说，增益矩阵 \(\bm{J}_{k}\) 与 \(\bm{D}_{\Delta}(k)\) 成反比，而与 \(\bm{D}_{\varOmega}(k-1)\) 成正比。

(4) 上述卡尔曼滤波方程是根据广义最小二乘原理导出的，在第 1 章已经证明，按广义最小二乘原理得到的结果具有无偏性和方差最小性。因此，卡尔曼滤波值 \(\hat{\bm{X}}(k/k)\) 是由观测值 \(\bm{L}_{1},\bm{L}_{2},\cdots,\bm{L}_{k}\) 得到 \(\bm{X}_{k}\) 的最小方差无偏估计，其误差方差阵 \(\bm{D}_{\hat{X}}(k/k)\) 就是所有估计中的最小方差阵。

例 4-3-1

设状态方程和观测方程为 \[\left.\begin{aligned} X_{k+1}&=0.5X_{k}+\varOmega_{k}\\ L_{k+1}&=X_{k+1}+\varDelta_{k+1} \end{aligned}\right\}\] 式中 \(X\)，\(L\)，\(\varOmega\)，\(\varDelta\) 都是标量，其随机模型是： \[\left.\begin{aligned} &E(\varDelta_{k})=0,\ E(\varOmega_{k})=0,\ E(X_{0})=\hat{X}(0/0)=0\\ &\operatorname{cov}(\varDelta_{k},\varDelta_{j})=0,\ \operatorname{cov}(\varOmega_{k},\varOmega_{j})=0,\ \operatorname{cov}(\varOmega_{k},\varDelta_{j})=0\\ &\operatorname{cov}(X_{0},\varDelta_{k})=0,\ \operatorname{cov}(X_{0},\varOmega_{k})=0\\ &D_{\varDelta}(k)=2,\ D_{\varOmega}(k)=1,\ D_{X}(0)=1 \end{aligned}\right\}\] 又已知两次观测数据 \(L_{1}=4\)，\(L_{2}=2\)，试求 \(\hat{X}(2/2)\) 和 \(D_{\hat{X}}(2/2)\)。

解按卡尔曼滤波方程 (4-3-27)、(4-3-28)、(4-3-29) 式计算：

(1) 计算一步预测值 \[\hat{X}(1/0)=\varPhi_{1,0}\hat{X}(0/0)=0\] \[D_{\hat{X}}(1/0)=\varPhi_{1,0}D_{X}(0)\varPhi_{1,0}^{T}+D_{\varOmega}(0)=0.5^{2}+1=1.25\]

(2) 求增益矩阵 \(J_{1}\) \[\begin{aligned} J_{1}&=D_{\hat{X}}(1/0)B_{1}^{T}\left[B_{1}D_{\hat{X}}(1/0)B_{1}^{T}+D_{\varDelta}(1)\right]^{-1}\\ &=1.25\,(1.25+2)^{-1}=\frac{1.25}{3.25}=0.385 \end{aligned}\]

(3) 计算“新息” \(\bar{L}(1/0)\) \[\bar{L}(1/0)=L_{1}-B_{1}\hat{X}(1/0)=4\]

(4) 计算 \(\hat{X}(1/1)\) 和 \(D_{\hat{X}}(1/1)\) \[\hat{X}(1/1)=\hat{X}(1/0)+J_{1}\bar{L}(1/0)=0+0.385\times 4=1.54\] \[D_{\hat{X}}(1/1)=(1-J_{1}B_{1})\,D_{\hat{X}}(1/0)=(1-0.385)\times 1.25=0.77\]

(5) 依 (1)   (4) 的公式计算 \(\hat{X}(2/1)\)，\(D_{\hat{X}}(2/1)\)，\(J_{2}\)，\(\bar{L}(2/1)\) 和 \(\hat{X}(2/2)\)，\(D_{\hat{X}}(2/2)\)，得 \[\hat{X}(2/1)=\varPhi_{2,1}\hat{X}(1/1)=0.5\times 1.54=0.77\] \[D_{\hat{X}}(2/1)=\varPhi_{2,1}D_{\hat{X}}(1/1)\varPhi_{2,1}^{T}+D_{\varOmega}(1)=0.5^{2}\times 0.77+1=1.19\] \[J_{2}=D_{\hat{X}}(2/1)B_{2}^{T}\left[B_{2}D_{\hat{X}}(2/1)B_{2}^{T}+D_{\varDelta}(2)\right]^{-1} =1.19\times(1.19+2)^{-1}=0.373\] \[\bar{L}(2/1)=L_{2}-B_{2}\hat{X}(2/1)=2-0.77=1.23\] \[\hat{X}(2/2)=\hat{X}(2/1)+J_{2}\bar{L}(2/1)=0.77+0.373\times 1.23=1.23\] \[D_{\hat{X}}(2/2)=(1-J_{2}B_{2})\,D_{\hat{X}}(2/1)=(1-0.373)\times 1.19=0.75\]

例 4-3-2

设有标量系统 \[\begin{gathered} X_{k+1}=X_{k}+\varOmega_{k}\\ L_{k+1}=X_{k+1}+\varDelta_{k+1} \end{gathered}\] \(\varOmega_{k}\) 与 \(\varDelta_{k}\) 均为零均值白噪声序列，且它们完全不相关，并与 \(X_{0}\) 相互独立，已知 \[\begin{gathered} D_{\varOmega}(k)=25,\ D_{\varDelta}(k)=15\\ D_{X}(0)=100,\ E(X_{0})=\hat{X}(0/0)=0 \end{gathered}\] 观测值为 \[L_{1}=10,\ L_{2}=20,\ L_{3}=30,\ L_{4}=40\] 求 \(\hat{X}(k/k)\) 和 \(D_{\hat{X}}(k/k)\)（\(k=1,2,3,4\)）。

解由滤波方程可得 \[\begin{gathered} \hat{X}(k/k)=\hat{X}(k-1/k-1)+J_{k}\left[L_{k}-\hat{X}(k-1/k-1)\right]\\ D_{\hat{X}}(k/k)=(1-J_{k})\,D_{\hat{X}}(k/k-1)\\ J_{k}=D_{\hat{X}}(k/k-1)\left[D_{\hat{X}}(k/k-1)+D_{\varDelta}(k)\right]^{-1}\\ D_{\hat{X}}(k/k-1)=D_{\hat{X}}(k-1/k-1)+D_{\varOmega}(k) \end{gathered}\] 按以上各式计算，其结果列于表 4-1。

动态测量系统的卡尔曼滤波

为说明卡尔曼滤波方法在测量中的应用，本节讨论动态测量系统卡尔曼滤波的基本模型，并举例说明实用模型的建立方法。

动态测量系统的基本模型

动态测量系统一般是由处于运动之中的一组地面点或空间点所构成。这些点的位置可看做时间 \(t\) 的函数。高程监测网是一维网，点的位置是高程；平面监测网是二维网，点的位置是平面坐标；三维地面网点或卫星的位置可以是空间坐标或大地坐标（大地经纬度和大地高），为了估计动态测量系统的运动状态，通常是以点的位置和它们的运动速率为状态向量。设点 \(i\) 在时刻 \(t\) 的位置向量为 \(\bm{\xi}_{i}(t)\)，其瞬时速率为 \(\bm{\lambda}_{i}(t)\)，而将瞬时加速率 \(\bm{\varOmega}_{i}(t)\) 看做是一种随机干扰，即视 \(\bm{\varOmega}_{i}(t)\) 为动态噪声，则有以下微分关系式： \[\left.\begin{aligned} \dot{\bm{\xi}}_{i}(t)&=\bm{\lambda}_{i}(t)\\ \dot{\bm{\lambda}}_{i}(t)&=\bm{\varOmega}_{i}(t) \end{aligned}\right\} \tag{4-4-1}\] 记 \(i\) 点的状态向量为 \(\underset{2d\times 1}{\bm{X}_{i}}(t)\)，全网的状态向量为 \(\bm{X}(t)\)，即令 \[\underset{2d\times 1}{\bm{X}_{i}}(t)=\begin{bmatrix}\bm{\xi}_{i}(t)\\ \bm{\lambda}_{i}(t)\end{bmatrix},\qquad \bm{X}(t)=\begin{bmatrix}\bm{X}_{1}(t)\\ \bm{X}_{2}(t)\\ \vdots\\ \bm{X}_{q}(t)\end{bmatrix}\] 式中 \(d\) 表示网的维数，\(q\) 为全网的待定点数。类似地记 \[\bm{\varOmega}(t)=\begin{bmatrix}\bm{\varOmega}_{1}(t)&\bm{\varOmega}_{2}(t)&\cdots&\bm{\varOmega}_{q}(t)\end{bmatrix}^{T}\] 则 (4-4-1) 式可写为 \[\dot{\bm{X}}_{i}(t)=\begin{bmatrix}\bm{0}&\bm{E}\\ \bm{0}&\bm{0}\end{bmatrix}\bm{X}_{i}(t) +\begin{bmatrix}\bm{0}\\ \bm{E}\end{bmatrix}\bm{\varOmega}_{i}(t) \tag{4-4-2}\] 全网的状态方程为 \[\dot{\bm{X}}(t)=\begin{bmatrix} \bm{0}&\bm{E}&\bm{0}&\bm{0}&\cdots\\ \bm{0}&\bm{0}&\bm{0}&\bm{0}&\cdots\\ \bm{0}&\bm{0}&\bm{0}&\bm{E}&\cdots\\ \bm{0}&\bm{0}&\bm{0}&\bm{0}&\cdots\\ \vdots&\vdots&\vdots&\vdots& \end{bmatrix}\bm{X}(t) +\begin{bmatrix} \bm{0}&\bm{0}&\cdots\\ \bm{E}&\bm{0}&\cdots\\ \bm{0}&\bm{0}&\cdots\\ \bm{0}&\bm{E}&\cdots\\ \vdots&\vdots& \end{bmatrix}\bm{\varOmega}(t) \tag{4-4-3}\] 式中 \(\bm{0}\) 是 \(d\) 阶零矩阵，\(\bm{E}\) 是 \(d\) 阶单位阵。

状态方程 (4-4-2) 是一个微分方程。为了将其化为离散形式，需要求此微分方程的解。可以看到，当 \(d=1\) 时，(4-4-2) 式与 [例 4-1-1] 中的状态方程相同，因此可以按比例的结果直接写出状态方程的解为 \[\bm{X}_{i}(t)=\begin{bmatrix}\bm{E}&(t-t_{0})\bm{E}\\ \bm{0}&\bm{E}\end{bmatrix}\bm{X}_{i}(t_{0}) +\int_{t_{0}}^{t}\begin{bmatrix}\bm{E}&(t-\tau)\bm{E}\\ \bm{0}&\bm{E}\end{bmatrix} \begin{bmatrix}\bm{0}\\ \bm{E}\end{bmatrix}\bm{\varOmega}_{i}(\tau)\,\mathrm{d}\tau\] 即 \[\bm{X}_{i}(t)=\begin{bmatrix}\bm{E}&(t-t_{0})\bm{E}\\ \bm{0}&\bm{E}\end{bmatrix}\bm{X}_{i}(t_{0}) +\int_{t_{0}}^{t}\begin{bmatrix}(t-\tau)\bm{E}\\ \bm{E}\end{bmatrix}\bm{\varOmega}_{i}(\tau)\,\mathrm{d}\tau\] 由 (4-2-4) 式知 \[\begin{aligned} \bm{\varGamma}_{i+1,k}&=\int_{t_{k}}^{t_{k+1}}\bm{\varPhi}_{i}(t_{k+1},\tau)\,\bm{F}(\tau)\,\mathrm{d}\tau\\ &=\int_{t_{k}}^{t_{k+1}}\begin{bmatrix}(t_{k+1}-\tau)\bm{E}\\ \bm{E}\end{bmatrix}\mathrm{d}\tau =\begin{bmatrix}\dfrac{1}{2}(t_{k+1}-t_{k})^{2}\bm{E}\\[8pt] (t_{k+1}-t_{k})\bm{E}\end{bmatrix} \end{aligned}\] 所以，(4-4-2) 式离散化的状态方程为 \[\bm{X}_{i,k+1}=\begin{bmatrix}\bm{E}&\Delta t_{k}\bm{E}\\ \bm{0}&\bm{E}\end{bmatrix}\bm{X}_{i,k} +\begin{bmatrix}\dfrac{1}{2}\Delta t_{k}^{2}\bm{E}\\[6pt] \Delta t_{k}\bm{E}\end{bmatrix}\bm{\varOmega}_{i,k} \tag{4-4-4}\] 式中 \(\Delta t_{k}=t_{k+1}-t_{k}\)。全网的状态方程为 \[\bm{X}_{k+1}=\begin{bmatrix} \bm{E}&\Delta t_{k}\bm{E}&\bm{0}&\bm{0}&\cdots\\ \bm{0}&\bm{E}&\bm{0}&\bm{0}&\cdots\\ \bm{0}&\bm{0}&\bm{E}&\Delta t_{k}\bm{E}&\cdots\\ \bm{0}&\bm{0}&\bm{0}&\bm{E}&\cdots\\ \vdots&\vdots&\vdots&\vdots& \end{bmatrix}\bm{X}_{k} +\begin{bmatrix} \dfrac{1}{2}\Delta t_{k}^{2}\bm{E}&\bm{0}&\cdots\\[6pt] \Delta t_{k}\bm{E}&\bm{0}&\cdots\\[4pt] \bm{0}&\dfrac{1}{2}\Delta t_{k}^{2}\bm{E}&\cdots\\[6pt] \bm{0}&\Delta t_{k}\bm{E}&\cdots\\ \vdots&\vdots& \end{bmatrix}\bm{\varOmega}_{k} \tag{4-4-5}\]

当运动的加速率呈系统性时，还可将加速率 \(\bm{\varOmega}_{i}(t)\) 视作状态，而将加速率的瞬时变化 \(\bm{\zeta}_{i}(t)\) 看做是随机干扰，记 \(i\) 点的状态向量为 \[\underset{1\times 3d}{\bm{X}_{i}^{T}}(t)=\begin{bmatrix}\underset{1\times d}{\bm{\xi}_{i}^{T}}(t)&\underset{1\times d}{\bm{\lambda}_{i}^{T}}(t)&\underset{1\times d}{\bm{\varOmega}_{i}^{T}}(t)\end{bmatrix}\] 则状态方程为 \[\dot{\bm{X}}_{i}(t)=\begin{bmatrix} \bm{0}&\bm{E}&\bm{0}\\ \bm{0}&\bm{0}&\bm{E}\\ \bm{0}&\bm{0}&\bm{0} \end{bmatrix}\bm{X}_{i}(t) +\begin{bmatrix}\bm{0}\\ \bm{0}\\ \bm{E}\end{bmatrix}\bm{\zeta}_{i}(t) \tag{4-4-6}\] 离散化后的状态方程为 \[\bm{X}_{i,k+1}=\begin{bmatrix} \bm{E}&\Delta t_{k}\bm{E}&\dfrac{1}{2}\Delta t_{k}^{2}\bm{E}\\[6pt] \bm{0}&\bm{E}&\Delta t_{k}\bm{E}\\ \bm{0}&\bm{0}&\bm{E} \end{bmatrix}\bm{X}_{i,k} +\begin{bmatrix} \dfrac{1}{6}\Delta t_{k}^{3}\bm{E}\\[6pt] \dfrac{1}{2}\Delta t_{k}^{2}\bm{E}\\[6pt] \Delta t_{k}\bm{E} \end{bmatrix}\bm{\zeta}_{i}(t) \tag{4-4-7}\] 以上讨论的状态方程是在假定地面上各点的动态变化彼此不相关的情况下得到的。可以说，这是一种基本的函数模型。为不同目的建立的各种动态测量系统，其动态过程都受到不同因素的影响，因此，它们的动态变化有着各自不同的特点，各点的动态变化之间可能存在各种不同的规律和关系。为了更有效地描述它们的动态变化，应该在建立状态方程时考虑这些规律和关系。下面以动态水准监测系统的几种情况为例，简述建立动态测量系统状态方程的具体方法。

动态水准监测系统

(1) 以点的高程及其速率为状态方程

设点 \(i\) 的状态向量为 \(\bm{X}_{i}^{T}(t)=\begin{bmatrix}H_{i}(t)&\lambda_{i}(t)\end{bmatrix}^{T}\)，由 (4-4-4) 式可写出离散化后 \(i\) 点在第 \(k+1\) 期的状态方程为 \[\bm{X}_{i,k+1}=\begin{bmatrix}1&\Delta t_{k}\\ 0&1\end{bmatrix}\bm{X}_{i,k} +\begin{bmatrix}\dfrac{1}{2}\Delta t_{k}^{2}\\[6pt] \Delta t_{k}\end{bmatrix}\varOmega_{i,k} \tag{4-4-8}\] 若记全网 \(q\) 个待定点的状态向量和动态噪声向量为 \[\begin{gathered} \underset{1\times 2q}{\bm{X}^{T}}(t)=\begin{bmatrix}H_{1}(t),\cdots,H_{q}(t)&\lambda_{1}(t),&\cdots,&\lambda_{q}(t)\end{bmatrix}\\ \underset{1\times q}{\bm{\varOmega}^{T}}(t)=\begin{bmatrix}\varOmega_{1}(t)&\cdots&\varOmega_{q}(t)\end{bmatrix} \end{gathered}\] 则可由 (4-4-5) 式得到全网在第 \(k+1\) 期的状态方程 \[\underset{2q\times 1}{\bm{X}_{k+1}}=\begin{bmatrix}\bm{E}&\Delta\bm{T}_{k}\\ \bm{0}&\bm{E}\end{bmatrix}\bm{X}_{k} +\begin{bmatrix}\dfrac{1}{2}\Delta\bm{T}_{k}^{2}\\[6pt] \Delta\bm{T}_{k}\end{bmatrix}\bm{\varOmega}_{k} \tag{4-4-9}\] 式中 \(\bm{0}\)，\(\bm{E}\)，\(\Delta\bm{T}_{k}\) 均为 \(q\) 阶方阵，且 \[\Delta\bm{T}_{k}=\Delta t_{k}\cdot\bm{E}\] 如果认为第 \(k+1\) 期的观测高差都是在 \(t_{k+1}\) 时刻观测的，则 \(h_{ij}\) 的观测方程为 \[h_{ij/k+1}=-H_{i,k+1}+H_{j,k+1}+\varDelta_{ij/k+1} \tag{4-4-10}\] 其形式与经典平差相同。如果各个观测高差的观测时刻 \(t_{ij/k+1}\) 不相同，设第 \(k+1\) 期的中心时刻为 \(t_{k+1}\)，则观测方程为 \[h_{ij/k+1}=-H_{i,k+1}+H_{j,k+1}-\Delta t_{ij/k+1}\lambda_{i}+\Delta t_{ij/k+1}\lambda_{j}+\varDelta_{ij/k+1} \tag{4-4-11}\] 式中，\(\Delta t_{ij/k+1}=t_{ij/k+1}-t_{k+1}\)。

(2) 以点的高程和区域速率为状态参数

假定可以将动态水准监测网分为两区（如图 4-1），例如断裂带的两侧，水库坝区的坝上区与坝下区等。此时可以对每个水准点设一个高程参数，每个区仅设一个速率参数。记全网的状态向量和动态噪声向量为 \[\begin{gathered} \bm{X}^{T}(t)=\begin{bmatrix}H_{1}(t)&H_{2}(t)&\cdots&H_{q}(t)&\lambda_{\mathrm{I}}(t)&\lambda_{\mathrm{II}}(t)\end{bmatrix}\\ \bm{\varOmega}^{T}(t)=\begin{bmatrix}\varOmega_{1}(t)&\varOmega_{\mathrm{II}}(t)\end{bmatrix} \end{gathered}\] 则离散化的第 \(k+1\) 期的状态方程为 \[\bm{X}_{k+1}=\begin{bmatrix} \underset{q\times q}{\bm{E}}&\underset{q\times 2}{\Delta\bm{T}_{k}}\\ \underset{2\times q}{\bm{0}}&\underset{2\times 2}{\bm{E}} \end{bmatrix}\bm{X}_{k} +\begin{bmatrix} \underset{q\times 2}{\Delta\bm{T}_{k}^{\prime}}\\ \Delta t\cdot\underset{2\times 2}{\bm{E}} \end{bmatrix}\bm{\varOmega}_{k} \tag{4-4-12}\]

式中 \[\begin{gathered} \Delta\bm{T}_{k}=\Delta t_{k}\begin{bmatrix}\bm{e}_{\mathrm{I}}&\bm{e}_{\mathrm{II}}\end{bmatrix}\\ \Delta\bm{T}_{k}^{\prime}=\frac{1}{2}\Delta t_{k}^{2}\begin{bmatrix}\bm{e}_{\mathrm{I}}&\bm{e}_{\mathrm{II}}\end{bmatrix}\\ \bm{e}_{\mathrm{I}}^{T}=\begin{bmatrix}\underbrace{1\quad\cdots}_{q_{1}\text{ 个}}1&\underbrace{0\quad\cdots}_{q_{2}\text{ 个}}0\end{bmatrix}\\ \bm{e}_{\mathrm{II}}^{T}=\begin{bmatrix}\underbrace{0\quad\cdots}_{q_{1}\text{ 个}}0&\underbrace{1\quad\cdots}_{q_{2}\text{ 个}}1\end{bmatrix} \end{gathered}\] \(q_{1}\) 和 \(q_{2}\) 分别为 \(\mathrm{I}\)、\(\mathrm{II}\) 两个区中的点数，\(q=q_{1}+q_{2}\)。

当第 \(k+1\) 期两区的观测高差可认为是在同一时刻（时间段）\(t_{k+1}\) 观测的，其观测方程仍然是 (4-4-10) 式。当各个观测高差的观测时刻 \(t_{ij/k+1}\) 不相同时，对于同一区内两点间的观测高差的观测方程也与 (4-4-10) 式相同；而当 \(i\) 点在 \(\mathrm{I}\) 区，\(j\) 点在 \(\mathrm{II}\) 区时，其观测方程为 \[h_{ij/i+1}=-H_{i,k+1}+H_{j,k+1}-\Delta t_{ij/k+1}\lambda_{\mathrm{I}}+\Delta t_{ij/k+1}\lambda_{\mathrm{II}}+\varDelta_{ij/k+1} \tag{4-4-13}\]

(3) 以点的高程和呈线性关系的高程速率为状态参数

设某动态水准监测网如图 4-2。假定其重心点的高程速率参数为 \(\lambda_{1}\)，而其它各点的高程变化则与该点到重心点的距离 \(S_{i}\) 呈线性关系。大型矿区和地震危险区就可能是这类情况，此时，对于 \(i\) 点有 \[\dot{H}_{i}(t)=\lambda_{i}(t)=\lambda_{1}(t)+S_{i}\lambda_{2}(t) \tag{4-4-14}\] 若设全网的状态向量和动态噪声向量为 \[\begin{gathered} \bm{X}^{T}(t)=\begin{bmatrix}H_{1}(t)&\cdots&H_{q}(t)&\lambda_{1}(t)&\lambda_{2}(t)\end{bmatrix}\\ \bm{\varOmega}^{T}(t)=\begin{bmatrix}\varOmega_{1}(t)&\varOmega_{2}(t)\end{bmatrix} \end{gathered}\] 则有状态方程 \[\underset{(q+2)\times 1}{\dot{\bm{X}}}(t)=\begin{bmatrix} \underset{q\times q}{\bm{0}}&\underset{q\times 2}{\bm{S}}\\ \underset{2\times q}{\bm{0}}&\underset{2\times 2}{\bm{0}} \end{bmatrix}\bm{X}(t) +\begin{bmatrix} \underset{q\times 2}{\bm{0}}\\ \underset{2\times 2}{\bm{E}} \end{bmatrix}\bm{\varOmega}(t) \tag{4-4-15}\]

式中 \[\bm{S}=\begin{bmatrix} 1&1&\cdots&1\\ S_{1}&S_{2}&\cdots&S_{q} \end{bmatrix}^{T}\] 离散化后的第 \(k+1\) 期的状态方程为 \[\bm{X}_{k+1}=\begin{bmatrix}\bm{E}&\Delta\bm{T}_{k_{1}}\\ \bm{0}&\bm{E}\end{bmatrix}\bm{X}_{k} +\begin{bmatrix}\Delta\bm{T}_{kj}^{\prime}\\ \Delta\bm{T}_{k}\end{bmatrix}\bm{\varOmega}_{k} \tag{4-4-16}\] 式中 \[\begin{gathered} \Delta\bm{T}_{k_{1}}=\Delta t_{k}\bm{S}\\ \Delta\bm{T}_{k}=\Delta t_{k}\cdot\bm{E}\\ \Delta\bm{T}_{k_{1}}^{\prime}=\frac{1}{2}\Delta t_{k}^{2}\bm{S} \end{gathered}\] 当 \(h_{ij}\) 的观测时刻 \(t_{ij/k+1}\) 不等于第 \(k+1\) 期的中心时刻 \(t_{k+1}\) 时，其观测方程为 \[h_{ij/k+1}=-H_{ij/k+1}+H_{j/k+1}+\Delta t_{ij/k+1}(-S_{i}+S_{j})\,\lambda_{2}+\varDelta_{ij/k+1} \tag{4-4-17}\]

离散型卡尔曼滤波的推广

在上面讨论的滤波问题中，假定动态噪声和观测噪声是完全不相关的白噪声序列，在这种情况下导出的卡尔曼滤波方程可以说是卡尔曼滤波的基本方程，本节进一步考虑白噪声作用下的一般线性系统，及有色噪声（或称为相关噪声）的情况，所谓有色噪声，指的是噪声序列中每一时刻的噪声与另一时刻的噪声相关。这里所考虑的是有色噪声序列的一种特殊情况，即只考虑一种由白噪声驱动的线性系统所产生的相关噪声，这对于应用来说已足够了。例如，地面雷达站对运动目标进行跟踪观测时，其观测噪声可以认为是这种有色噪声；惯性导航系统的主要动态噪声——陀螺漂移，也可认为是这种有色噪声。

白噪声作用下一般线性系统的滤波

白噪声作用下一般线性系统的状态方程和观测方程仍然是 (4-3-1)、(4-3-2) 式，即 \[\left.\begin{aligned} \bm{X}_{k+1}&=\bm{\varPhi}_{k+1,k}\bm{X}_{k}+\bm{\varPsi}_{k+1,k}\bm{U}_{k}+\bm{\varGamma}_{k+1,k}\bm{\varOmega}_{k}\\ \bm{L}_{k+1}&=\bm{B}_{k+1}\bm{X}_{k+1}+\bm{Z}_{k+1}+\bm{\Delta}_{k+1} \end{aligned}\right\} \tag{4-5-1}\] 所不同的仅仅是在随机模型中假定 \[\left.\begin{aligned} E(\bm{\varOmega}_{k})&=\bm{\mu}_{\varOmega}(k)\neq\bm{0}\\ E(\bm{\Delta}_{k})&=\bm{\mu}_{\Delta}(k)\neq\bm{0}\\ \operatorname{cov}(\bm{\varOmega}_{k},\bm{\Delta}_{j})&=\bm{D}_{\varOmega\Delta}(k)\,\delta_{kj} \end{aligned}\right\} \tag{4-5-2}\] 即假定动态噪声和观测噪声的均值不为零，且它们在同一时刻相关，其余的性质均与 4-3 节相同。

为了导出白噪声作用下一般线性系统的卡尔曼滤波方程，我们先将 \(\bm{\Delta}_{k}\)，\(\bm{\varOmega}_{k}\) 进行变换，以便得到另一组均值为零且完全不相关的观测噪声和动态噪声序列，假定它们是 \[\left.\begin{aligned} \bm{\Delta}_{k}^{*}&=\bm{\Delta}_{k}-\bm{\mu}_{\Delta}(k)\\ \bm{\varepsilon}_{k}^{*}&=\bm{\varGamma}_{k+1,k}\left[\bm{\varOmega}_{k}-\bm{\mu}_{\varOmega}(k)\right]-\bm{C}_{k}\left[\bm{\Delta}_{k}-\bm{\mu}_{\Delta}(k)\right] \end{aligned}\right\} \tag{4-5-3}\] 式中 \(\bm{C}_{k}\) 为待定系数阵。不难看到，此时有 \[\left.\begin{aligned} E(\bm{\Delta}_{k}^{*})&=\bm{0}\\ E(\bm{\varepsilon}_{k}^{*})&=\bm{0} \end{aligned}\right\} \tag{4-5-4}\] 而欲使 \(\bm{\varepsilon}_{k}^{*}\) 与 \(\bm{\Delta}_{k}^{*}\) 完全不相关，就应使 \[\operatorname{cov}(\bm{\varepsilon}_{k}^{*},\bm{\Delta}_{k}^{*}) =\bm{\varGamma}_{k+1,k}\bm{D}_{\varOmega\Delta}(k)-\bm{C}_{k}\bm{D}_{\Delta}(k)=\bm{0}\] 故得 \[\bm{C}_{k}=\bm{\varGamma}_{k+1,k}\bm{D}_{\varOmega\Delta}(k)\,\bm{D}_{\Delta}^{-1}(k) \tag{4-5-5}\] 此时，由 (4-5-5) 式按协方差传播律还可得 \[\left.\begin{aligned} \bm{D}_{\Delta}^{*}(k)&=\bm{D}_{\Delta}(k)\\ \bm{D}_{\varepsilon}^{*}(k)&=\begin{bmatrix}\bm{\varGamma}_{k+1,k}&-\bm{C}_{k}\end{bmatrix} \begin{bmatrix}\bm{D}_{\varOmega}(k)&\bm{D}_{\varOmega\Delta}(k)\\ \bm{D}_{\Delta\varOmega}(k)&\bm{D}_{\Delta}(k)\end{bmatrix} \begin{bmatrix}\bm{\varGamma}_{k+1,k}^{T}\\ -\bm{C}_{k}^{T}\end{bmatrix}\\ &=\bm{\varGamma}_{k+1,k}\bm{D}_{\varOmega}(k)\,\bm{\varGamma}_{k+1,k}^{T}-\bm{C}_{k}\bm{D}_{\Delta}(k)\,\bm{C}_{k}^{T} \end{aligned}\right\} \tag{4-5-6}\]

顾及观测方程，由 (4-5-3) 式可写出 \[\left.\begin{aligned} \bm{\Delta}_{k}&=\bm{\Delta}_{k}^{*}+\bm{\mu}_{\Delta}(k)\\ \bm{\varGamma}_{k+1,k}\bm{\varOmega}_{k}&=\bm{\varepsilon}_{k}^{*}+\bm{\varGamma}_{k+1,k}\bm{\mu}_{\varOmega}(k) +\bm{C}_{k}\left[\bm{L}_{k}-\bm{B}_{k}\bm{X}_{k}-\bm{Z}_{k}-\bm{\mu}_{\Delta}(k)\right] \end{aligned}\right\} \tag{4-5-7}\] 再将上式代入 (4-5-1) 式，即可得到变换后的状态方程和观测方程 \[\left.\begin{aligned} \bm{X}_{k+1}&=\bm{\varPhi}_{k+1,k}^{*}\bm{X}_{k}+\bm{U}_{k}^{*}+\bm{\varepsilon}_{k}^{*}\\ \bm{L}_{k+1}&=\bm{B}_{k+1}\bm{X}_{k+1}+\bm{Z}_{k+1}^{*}+\bm{\Delta}_{k+1}^{*} \end{aligned}\right\} \tag{4-5-8}\] 式中 \[\left.\begin{aligned} \bm{\varPhi}_{k+1,k}^{*}&=\bm{\varPhi}_{k+1,k}-\bm{C}_{k}\bm{B}_{k}\\ \bm{U}_{k}^{*}&=\bm{\varPsi}_{k+1,k}\bm{U}_{k}+\bm{\varGamma}_{k+1,k}\bm{\mu}_{\varOmega}(k) +\bm{C}_{k}\left[\bm{L}_{k}-\bm{\mu}_{\Delta}(k)-\bm{Z}_{k}\right]\\ \bm{Z}_{k}^{*}&=\bm{Z}_{k}+\bm{\mu}_{\Delta}(k) \end{aligned}\right\} \tag{4-5-9}\] 可以看到，(4-5-8) 式中的动态噪声 \(\bm{\varepsilon}_{k}^{*}\) 与观测噪声 \(\bm{\Delta}_{k}^{*}\,(k\geq 0)\) 是完全不相关的零均值白噪声序列，因此，按 4-3 节的卡尔曼滤波方程，可以得到白噪声作用下一般线性系统的卡尔曼滤波方程为 \[\left.\begin{aligned} \hat{\bm{X}}(k/k)&=\hat{\bm{X}}(k/k-1) +\bm{J}_{k}\left[\bm{L}_{k}-\bm{\mu}_{\Delta}(k)-\bm{Z}_{k}-\bm{B}_{k}\hat{\bm{X}}(k/k-1)\right]\\ \bm{D}_{\hat{X}}(k/k)&=(\bm{E}-\bm{J}_{k}\bm{B}_{k})\,\bm{D}_{\hat{X}}(k/k-1) \end{aligned}\right\} \tag{4-5-10}\] 而一步预测公式和增益矩阵公式为 \[\left.\begin{aligned} \hat{\bm{X}}(k/k-1)&=\bm{\varPhi}_{k,k-1}^{*}\hat{\bm{X}}(k-1/k-1)+\bm{U}_{k-1}^{*}\\ \bm{D}_{\hat{X}}(k/k-1)&=\bm{\varPhi}_{k,k-1}^{*}\bm{D}_{\hat{X}}(k-1/k-1)\bm{\varPhi}_{k,k-1}^{*T}+\bm{D}_{\varepsilon}^{*}(k-1)\\ \bm{J}_{k}&=\bm{D}_{\hat{X}}(k/k-1)\bm{B}_{k}^{T} \left[\bm{B}_{k}\bm{D}_{\hat{X}}(k/k-1)\bm{B}_{k}^{T}+\bm{D}_{\Delta}(k)\right]^{-1} \end{aligned}\right\} \tag{4-5-11}\] 式中 \[\left.\begin{aligned} \bm{\varPhi}_{k,k-1}^{*}&=\bm{\varPhi}_{k,k-1}-\bm{C}_{k-1}\bm{B}_{k-1}\\ \bm{U}_{k-1}^{*}&=\bm{\varPsi}_{k,k-1}\bm{U}_{k-1}+\bm{\varGamma}_{k,k-1}\bm{\mu}_{\varOmega}(k-1) +\bm{C}_{k-1}\left[\bm{L}_{k-1}-\bm{\mu}_{\Delta}(k-1)-\bm{Z}_{k-1}\right]\\ \bm{D}_{\varepsilon}^{*}(k-1)&=\bm{\varGamma}_{k,k-1}\bm{D}_{\varOmega}(k-1)\bm{\varGamma}_{k,k-1}^{T} -\bm{C}_{k-1}\bm{D}_{\Delta}(k-1)\bm{C}_{k-1}^{T}\\ \bm{C}_{k-1}&=\bm{\varGamma}_{k,k-1}\bm{D}_{\varOmega\Delta}(k-1)\,\bm{D}_{\Delta}^{-1}(k-1) \end{aligned}\right\} \tag{4-5-12}\] 也可写出 (4-5-10) 式的另一种形式为 \[\left.\begin{aligned} \hat{\bm{X}}(k/k)={}&\left[\bm{B}_{k}^{T}\bm{D}_{\Delta}^{-1}(k)\bm{B}_{k}+\bm{D}_{X}^{-1}(k/k-1)\right]^{-1} \{\bm{D}_{X}^{-1}(k/k-1)\hat{\bm{X}}(k/k-1)+{}\\ &\bm{B}_{k}^{T}\bm{D}_{\Delta}^{-1}(k)\left[\bm{L}_{k}-\bm{\mu}_{\Delta}(k)-\bm{Z}_{k}\right]\}\\ \bm{D}_{\hat{X}}(k/k)={}&\left[\bm{B}_{k}^{T}\bm{D}_{\Delta}^{-1}(k)\bm{B}_{k}+\bm{D}_{X}^{-1}(k/k-1)\right]^{-1} \end{aligned}\right\} \tag{4-5-13}\] 显然，在上面的滤波方程中，当 \(\bm{\mu}_{\varOmega}(k)\)、\(\bm{\mu}_{\Delta}(k)\) 和 \(\bm{D}_{\varOmega\Delta}(k)\) 为零时，它们就变成了 4-3 节中的卡尔曼滤波方程。

动态噪声是白噪声所驱动的有色噪声的情形

对于动态噪声 \(\bm{\varOmega}_{k}\) 是白噪声所驱动的有色噪声的情形，可以采用扩大状态向量维数的方法来处理。设线性系统的状态方程和观测方程仍是 (4-5-1) 式，但其中的动态噪声 \(\bm{\varOmega}_{k}\) 是由白噪声序列 \(\bm{\eta}_{k}\,(k\geq 0)\) 所驱动的有色噪声，即 \[\bm{\varOmega}_{k}=\bm{F}_{k,k-1}\bm{\varOmega}_{k-1}+\bm{\eta}_{k-1} \tag{4-5-14}\] 且有 \[\left.\begin{aligned} &E(\bm{\eta}_{k})=\bm{0},\ \operatorname{var}(\bm{\eta}_{k})=\bm{D}_{\eta}(k)\\ &E(\bm{\varOmega}_{0})=\bm{0},\ \operatorname{var}(\bm{\varOmega}_{0})=\bm{D}_{\varOmega}(0)\\ &\operatorname{cov}(\bm{\varOmega}_{0},\bm{\eta}_{k})=\bm{0}\\ &\operatorname{cov}(\bm{X}_{0},\bm{\eta}_{k})=\bm{0},\ \operatorname{cov}(\bm{X}_{0},\bm{\varOmega}_{0})=\bm{0} \end{aligned}\right\} \tag{4-5-15}\] 其他性质与 4-3 节相同。

现将 \(\bm{\varOmega}_{k}\) 看成状态向量的一部分，即采用增广状态向量的方法，将 (4-5-1) 式写为 \[\left.\begin{aligned} \begin{bmatrix}\bm{X}_{k+1}\\ \bm{\varOmega}_{k+1}\end{bmatrix} &=\begin{bmatrix}\bm{\varPhi}_{k+1,k}&\bm{\varGamma}_{k+1,k}\\ \bm{0}&\bm{F}_{k+1,k}\end{bmatrix} \begin{bmatrix}\bm{X}_{k}\\ \bm{\varOmega}_{k}\end{bmatrix} +\begin{bmatrix}\bm{\varPsi}_{k+1,k}\\ \bm{0}\end{bmatrix}\bm{U}_{k} +\begin{bmatrix}\bm{0}\\ \bm{E}\end{bmatrix}\bm{\eta}_{k}\\ \bm{L}_{k+1}&=\begin{bmatrix}\bm{B}_{k+1}&\bm{0}\end{bmatrix} \begin{bmatrix}\bm{X}_{k+1}\\ \bm{\varOmega}_{k+1}\end{bmatrix} +\bm{Z}_{k+1}+\bm{\Delta}_{k+1} \end{aligned}\right\} \tag{4-5-16}\] 若记 \[\left.\begin{aligned} \bm{X}_{k}^{*}&=\begin{bmatrix}\bm{X}_{k}\\ \bm{\varOmega}_{k}\end{bmatrix},\quad \bm{\varPhi}_{k+1,k}^{*}=\begin{bmatrix}\bm{\varPhi}_{k+1,k}&\bm{\varGamma}_{k+1,k}\\ \bm{0}&\bm{F}_{k+1,k}\end{bmatrix}\\ \bm{\varPsi}_{k+1,k}^{*}&=\begin{bmatrix}\bm{\varPsi}_{k+1,k}\\ \bm{0}\end{bmatrix},\quad \bm{G}_{k}=\begin{bmatrix}\bm{0}\\ \bm{E}\end{bmatrix},\quad \bm{B}_{k}^{*}=\begin{bmatrix}\bm{B}_{k}&\bm{0}\end{bmatrix} \end{aligned}\right\} \tag{4-5-17}\] 则得到一组新的动态方程和观测方程 \[\left.\begin{aligned} \bm{X}_{k+1}^{*}&=\bm{\varPhi}_{k+1,k}^{*}\bm{X}_{k}^{*}+\bm{\varPsi}_{k+1,k}^{*}\bm{U}_{k}+\bm{G}_{k}\bm{\eta}_{k}\\ \bm{L}_{k+1}&=\bm{B}_{k+1}^{*}\bm{X}_{k+1}^{*}+\bm{Z}_{k+1}+\bm{\Delta}_{k+1} \end{aligned}\right\} \tag{4-5-18}\] 此时，初始状态 \(\bm{X}_{0}^{*}\) 的方差阵即为 \[\operatorname{var}(\bm{X}_{0}^{*})=\begin{bmatrix}\bm{D}_{X}(0)&\bm{0}\\ \bm{0}&\bm{D}_{\varOmega}(0)\end{bmatrix}\] (4-5-18) 式的随机模型的形式与系统 (4-3-1) 是一致的，故可按 4-3 节中的卡尔曼滤波方程进行递推计算。由此可见，用扩大状态向量维数的方法，可以将这种有色噪声化为白噪声的情形。

观测噪声是白噪声所驱动的有色噪声的情形

设系统 (4-5-1) 的动态噪声 \(\bm{\varOmega}_{k}\) 是白噪声序列，但观测噪声 \(\bm{\Delta}_{k}\) 是白噪声 \(\bm{\xi}_{k}\) 所驱动的有色噪声序列，即 \[\bm{\Delta}_{k+1}=\bm{H}_{k+1,k}\bm{\Delta}_{k}+\bm{\xi}_{k} \tag{4-5-19}\] 其中 \(\bm{\xi}_{k}\) 是零均值白噪声序列，且有 \[\left.\begin{aligned} &\operatorname{var}(\bm{\xi}_{k})=\bm{D}_{\xi}(k),\ \operatorname{var}(\bm{\Delta}_{0})=\bm{D}_{\Delta}(0)\\ &\operatorname{cov}(\bm{\xi}_{k},\bm{\Delta}_{0})=\bm{0},\ \operatorname{cov}(\bm{\xi}_{k},\bm{\varOmega}_{j})=\bm{0}\\ &\operatorname{cov}(\bm{X}_{0},\bm{\xi}_{k})=\bm{0},\ \operatorname{cov}(\bm{X}_{0},\bm{\Delta}_{0})=\bm{0} \end{aligned}\right\} \tag{4-5-20}\] 其他性质均与 4-3 节相同。

对于这种情况，当然也可以按对有色动态噪声的处理方法，将观测噪声 \(\bm{\Delta}_{k}\) 当作状态向量的一部分，用扩大状态向量维数的方法处理。但这种方法的缺点是计算工作量增加，且求增益矩阵所要计算的逆阵不一定存在，因此，这种方法不太理想。

下面介绍另一种方法。这种方法的特点是，先将观测向量、观测噪声和观测方程作适当的变换，使原来的有色噪声下的线性系统变为白噪声作用下的一般线性系统，然后再利用滤波方程来导出色观测噪声下的线性系统的滤波方程。

首先令 \[\bm{L}_{k}^{*}=\bm{L}_{k+1}-\bm{H}_{k+1,k}\bm{L}_{k} \tag{4-5-21}\] 将 (4-5-1) 中的观测方程代入上式得 \[\bm{L}_{k}^{*}=\bm{B}_{k+1}\bm{X}_{k+1}+\bm{Z}_{k+1}+\bm{\Delta}_{k+1} -\bm{H}_{k+1,k}(\bm{B}_{k}\bm{X}_{k}+\bm{Z}_{k}+\bm{\Delta}_{k})\] 再将 (4-5-19) 式代入，并顾及 (4-5-1) 式中的状态方程，有 \[\begin{aligned} \bm{L}_{k}^{*}={}&\bm{B}_{k+1}(\bm{\varPhi}_{k+1,k}\bm{X}_{k}+\bm{\varPsi}_{k+1,k}\bm{U}_{k}+\bm{\varGamma}_{k+1,k}\bm{\varOmega}_{k}) +\bm{Z}_{k+1}+{}\\ &\bm{H}_{k+1,k}\bm{\Delta}_{k}+\bm{\xi}_{k} -\bm{H}_{k+1,k}(\bm{B}_{k}\bm{X}_{k}+\bm{Z}_{k})-\bm{H}_{k+1,k}\bm{\Delta}_{k} \end{aligned}\] 再记 \[\left.\begin{aligned} \bm{B}_{k}^{*}&=\bm{B}_{k+1}\bm{\varPhi}_{k+1,k}-\bm{H}_{k+1,k}\bm{B}_{k}\\ \bm{Z}_{k}^{*}&=\bm{B}_{k+1}\bm{\varPsi}_{k+1,k}\bm{U}_{k}+\bm{Z}_{k+1}-\bm{H}_{k+1,k}\bm{Z}_{k}\\ \bm{\Delta}_{k}^{*}&=\bm{B}_{k+1}\bm{\varGamma}_{k+1,k}\bm{\varOmega}_{k}+\bm{\xi}_{k} \end{aligned}\right\} \tag{4-5-22}\] 则得到一个新的线性系统 \[\left.\begin{aligned} \bm{X}_{k+1}&=\bm{\varPhi}_{k+1,k}\bm{X}_{k}+\bm{\varPsi}_{k+1,k}\bm{U}_{k}+\bm{\varGamma}_{k+1,k}\bm{\varOmega}_{k}\\ \bm{L}_{k+1}^{*}&=\bm{B}_{k+1}^{*}\bm{X}_{k+1}+\bm{Z}_{k+1}^{*}+\bm{\Delta}_{k+1}^{*} \end{aligned}\right\} \tag{4-5-23}\] 在这个系统中，由于 \(\bm{\varOmega}_{k}\) 和 \(\bm{\xi}_{k}\) 都是白噪声序列，所以，新的观测噪声 \(\bm{\Delta}_{k}^{*}\) 也是零均值白噪声序列，但此时 \(\bm{\Delta}_{k}^{*}\) 与 \(\bm{\varOmega}_{k}\) 在同一时刻是相关的。由协方差传播律可得 \(\bm{\Delta}_{k}^{*}\) 的方差及 \(\bm{\varOmega}_{k}\) 与 \(\bm{\Delta}_{k}^{*}\) 的协方差为 \[\left.\begin{aligned} \operatorname{var}(\bm{\Delta}_{k}^{*})=\bm{D}_{\Delta}^{*}(k) &=\bm{B}_{k+1}\bm{\varGamma}_{k+1,k}\bm{D}_{\varOmega}(k)\bm{\varGamma}_{k+1,k}^{T}\bm{B}_{k+1}^{T}+\bm{D}_{\xi}(k)\\ \operatorname{cov}(\bm{\varOmega}_{k},\bm{\Delta}_{k}^{*})=\bm{D}_{\varOmega\Delta}^{*}(k) &=(\bm{B}_{k+1}\bm{\varGamma}_{k+1,k}\bm{D}_{\varOmega}(k))^{T} \end{aligned}\right\} \tag{4-5-24}\] 也就是说，系统 (4-5-23) 是一个白噪声作用下一般线性系统。

根据白噪声作用下的一般线性系统的滤波方程，就可以写出观测噪声白噪声驱动的有色噪声情形下的卡尔曼滤波方程为 \[\left.\begin{aligned} \hat{\bm{X}}(k/k)&=\hat{\bm{X}}(k/k-1) +\bm{J}_{k}\left[\bm{L}_{k}^{*}-\bm{Z}_{k}-\bm{B}_{k}^{*}\hat{\bm{X}}(k/k-1)\right]\\ \bm{D}_{\hat{X}}(k/k)&=(\bm{E}-\bm{J}_{k}\bm{B}_{k}^{*})\,\bm{D}_{\hat{X}}(k/k-1) \end{aligned}\right\} \tag{4-5-25}\] 一步预测公式和增益矩阵为 \[\left.\begin{aligned} \hat{\bm{X}}(k/k-1)&=\bm{\varPhi}_{k,k-1}^{*}\hat{\bm{X}}(k-1/k-1)+\bm{U}_{k-1}^{*}\\ \bm{D}_{\hat{X}}(k/k-1)&=\bm{\varPhi}_{k,k-1}^{*}\bm{D}_{\hat{X}}(k-1/k-1)\bm{\varPhi}_{k,k-1}^{*T}+\bm{D}_{\varepsilon}^{*}(k-1)\\ \bm{J}_{k}&=\bm{D}_{\hat{X}}(k/k-1)\bm{B}_{k}^{*T} \left[\bm{B}_{k}^{*}\bm{D}_{\hat{X}}(k/k-1)\bm{B}_{k}^{*T}+\bm{D}_{\Delta}^{*}(k)\right]^{-1} \end{aligned}\right\} \tag{4-5-26}\] 式中 \[\left.\begin{aligned} \bm{\varPhi}_{k,k-1}^{*}&=\bm{\varPhi}_{k,k-1}-\bm{C}_{k-1}\bm{B}_{k-1}^{*}\\ \bm{U}_{k-1}^{*}&=\bm{\varPsi}_{k,k-1}\bm{U}_{k-1}+\bm{C}_{k-1}(\bm{L}_{k-1}^{*}-\bm{Z}_{k-1})\\ \bm{D}_{\varepsilon}^{*}(k-1)&=\bm{\varGamma}_{k,k-1}\bm{D}_{\varOmega}(k-1)\bm{\varGamma}_{k,k-1}^{T} +\bm{C}_{k-1}\bm{D}_{\Delta}^{*}(k-1)\bm{C}_{k-1}^{T}\\ \bm{C}_{k-1}&=\bm{\varGamma}_{k,k-1}\bm{D}_{\varOmega\Delta}^{*}(k-1)\,\bm{D}_{\Delta}^{-1}(k-1) \end{aligned}\right\} \tag{4-5-27}\] 上式中的 \(\bm{B}_{k-1}^{*}\) 按 (4-5-22) 式求定，而 \(\bm{D}_{\Delta}^{*}(k-1)\) 和 \(\bm{D}_{\varOmega\Delta}^{*}(k-1)\) 按 (4-5-24) 式求定。

注意，当 \(k=1\) 时，应有 \[\left.\begin{aligned} \hat{\bm{X}}(1/1)&=\bm{\varPhi}_{1,0}\hat{\bm{X}}_{0}+\bm{J}_{1}(\bm{L}_{1}-\bm{B}_{1}\bm{\varPhi}_{1,0}\hat{\bm{X}}_{0})\\ \bm{D}_{\hat{X}}(1/1)&=(\bm{E}-\bm{J}_{1}\bm{B}_{1})\,\bm{D}_{\hat{X}}(1/0) \end{aligned}\right\} \tag{4-5-28}\] \[\left.\begin{aligned} \bm{D}_{\hat{X}}(1/0)&=\bm{\varPhi}_{1,0}\bm{D}_{X}(0)\bm{\varPhi}_{1,0}^{T} +\bm{\varGamma}_{1,0}\bm{D}_{\varOmega}(0)\bm{\varGamma}_{1,0}^{T}\\ \bm{J}_{1}&=\bm{D}_{\hat{X}}(1/0)\bm{B}_{1}^{T} \left[\bm{B}_{1}\bm{D}_{\hat{X}}(1/0)\bm{B}_{1}^{T}+\bm{D}_{\Delta}^{*}(1)\right]^{-1} \end{aligned}\right\} \tag{4-5-29}\] 若在系统 (4-5-1) 中取 \(k\geq 0\)，且设初始状态 \(\bm{X}_{0}\) 的先验期望和先验方差分别为 \(\bm{\mu}_{0}\) 和 \(\bm{D}_{0}\)，则当 \(k\geq 1\) 时，\(\hat{\bm{X}}(k/k)\) 和 \(\bm{D}_{\hat{X}}(k/k)\) 均按 (4-5-25)   (4-5-27) 式计算，而对于 \(k=0\) 应有 \[\left.\begin{aligned} \hat{\bm{X}}_{0}&=\bm{\mu}_{0}+\bm{D}_{0}\bm{B}_{0}^{T} \left[\bm{B}_{0}\bm{D}_{0}\bm{B}_{0}^{T}+\bm{D}_{\Delta}(0)\right]^{-1}(\bm{L}_{0}-\bm{B}_{0}\bm{\mu}_{0})\\ \bm{D}_{\hat{X}}(0)&=\bm{D}_{0}-\bm{D}_{0}\bm{B}_{0}^{T} \left[\bm{B}_{0}\bm{D}_{0}\bm{B}_{0}^{T}+\bm{D}_{\Delta}(0)\right]^{-1}\bm{B}_{0}\bm{D}_{0}\\ &=\left[\bm{D}_{0}^{-1}+\bm{B}_{0}^{T}\bm{D}_{\Delta}^{-1}(0)\bm{B}_{0}^{T}\right]^{-1} \end{aligned}\right\} \tag{4-5-30}\]

离散线性系统的预测

如 4-3 节所述，在利用观测向量 \(\bm{L}_{1},\bm{L}_{2},\cdots,\bm{L}_{k}\) 求定在第 \(j\) 时刻的状态向量 \(\bm{X}_{j}\) 的最佳估值时，若 \(j>k\)，则称 \(\hat{\bm{X}}(j/k)\) 为最佳预测值，称求定 \(\hat{\bm{X}}(j/k)\) 和它的误差方差阵的过程为预测。根据 \(k\) 和 \(j\) 的变化情况，预测问题可以分为三种类型：

(1) 固定区间预测：固定 \(k=N\)，而变化 \(j\)，即由固定的观测向量 \(\bm{L}_{1},\bm{L}_{2},\cdots,\bm{L}_{N}\)，求不同时刻状态向量的预测值 \(\hat{\bm{X}}(j/N)\,(j>N)\)；

(2) 固定点预测：固定 \(j=N\)，而变化 \(k\)，即由变化的观测向量 \(\bm{L}_{1},\bm{L}_{2},\cdots,\bm{L}_{k}\)，求某个固定时刻的状态向量的预测值 \(\hat{\bm{X}}(N/k)\,(0<k<j=N)\)；

(3) 固定超前预测：\(j\) 和 \(k\) 都变化，但保持 \(j=k+N\)，即由变化的观测向量 \(\bm{L}_{1},\bm{L}_{2},\cdots,\bm{L}_{k}\) 求预测值 \(\hat{\bm{X}}(k+N/k)\)。

上述预测问题，实际上就是在已知前二阶矩特性的情况下，对一个马尔柯夫或高斯-马尔柯夫序列进行预测的问题。这里仅按广义最小二乘原理来推导 4-3 节的白噪声情形下离散线性系统的预测方程。

按照 4-3 节推导卡尔曼滤波方程的方法，在根据广义最小二乘原理和逐次平差，由观测值 \(\bm{L}_{1},\bm{L}_{2},\cdots,\bm{L}_{k}\) 求得 \(\hat{\bm{X}}(k/k)\) 之后不再增加观测值，还要求定预测值 \(\hat{\bm{X}}(j/k)\,(j>k)\)，也就是由以下误差方程进行最小二乘平差 \[\left.\begin{aligned} \bm{V}_{X}(k/k)&=\hat{\bm{X}}_{\varepsilon}(k/k)-\hat{\bm{X}}(k/k)\\ \bm{V}_{\varepsilon}(k)&=\bm{\varPhi}_{k+1,k}\hat{\bm{X}}_{\varepsilon}(k/k)-\hat{\bm{X}}(k+1/k)-\bm{L}_{\varepsilon}(k)\\ \bm{V}_{\varepsilon}(k+1)&=\bm{\varPhi}_{k+2,k+1}\hat{\bm{X}}(k+1/k)-\hat{\bm{X}}(k+2/k)-\bm{L}_{\varepsilon}(k+1)\\ &\cdots \end{aligned}\right\} \tag{4-6-1}\] 式中 \[\bm{L}_{\varepsilon}(k)=-\bm{\varPsi}_{k+1,k}\bm{U}_{k} \tag{4-6-2}\] 它对应的噪声和方差为 \[\left.\begin{aligned} \bm{\varepsilon}_{k}&=\bm{\varGamma}_{k+1,k}\bm{\varOmega}_{k}\\ \bm{D}_{\varepsilon}(k)&=\bm{\varGamma}_{k+1,k}\bm{D}_{\varOmega}(k)\bm{\varGamma}_{k+1,k}^{T} \end{aligned}\right\} \tag{4-6-3}\] 可以看到，在 (4-6-1) 式中，虚拟观测向量 \(\hat{\bm{X}}(k/k)\)，\(\bm{L}_{\varepsilon}(k)\)，\(\bm{L}_{\varepsilon}(k+1)\)，\(\cdots\) 中观测值的个数，与待定参数向量 \(\hat{\bm{X}}_{\varepsilon}(k/k)\)，\(\hat{\bm{X}}(k+1/k)\)，\(\hat{\bm{X}}(k+2/k)\)，\(\cdots\) 中参数的个数相等，也就是说，它们没有多余观测。因此，此时的改正数向量 \(\bm{V}_{X}(k/k)\)，\(\bm{V}_{\varepsilon}(k)\)，\(\bm{V}_{\varepsilon}(k+1)\)，\(\cdots\) 均为零，\(\hat{\bm{X}}_{\varepsilon}(k/k)=\hat{\bm{X}}(k/k)\)，而预测值为 \[\left.\begin{aligned} \hat{\bm{X}}(k+1/k)&=\bm{\varPhi}_{k+1,k}\hat{\bm{X}}(k/k)+\bm{\varPsi}_{k+1,k}\bm{U}_{k}\\ \hat{\bm{X}}(k+2/k)&=\bm{\varPhi}_{k+2,k+1}\hat{\bm{X}}(k+1/k)+\bm{\varPsi}_{k+2,k+1}\bm{U}_{k+1}\\ &=\bm{\varPhi}_{k+2,k+1}\left[\bm{\varPhi}_{k+1,k}\hat{\bm{X}}(k/k)+\bm{\varPsi}_{k+1,k}\bm{U}_{k}\right] +\bm{\varPsi}_{k+2,k+1}\bm{U}_{k+1}\\ &\cdots \end{aligned}\right\} \tag{4-6-4}\] 考虑 (4-1-26) 和 (4-2-4) 两式即 \[\left.\begin{aligned} \bm{\varPhi}_{j,k}&=\bm{\varPhi}_{j,i}\bm{\varPhi}_{i,k}\\ \bm{\varPsi}_{j,k}&=\int_{t_{k}}^{t_{j}}\bm{\varPhi}(t_{j},\tau)\,\bm{C}(\tau)\,\mathrm{d}\tau \end{aligned}\right\} \tag{4-6-5}\] (4-6-4) 式可变为 \[\left.\begin{aligned} \hat{\bm{X}}(k+1/k)&=\bm{\varPhi}_{k+1,k}\hat{\bm{X}}(k/k)+\bm{\varPsi}_{k+1,k}\bm{U}_{k}\\ \hat{\bm{X}}(k+2/k)&=\bm{\varPhi}_{k+2,k}\hat{\bm{X}}(k/k)+\bm{\varPsi}_{k+2,k}\bm{U}_{k}+\bm{\varPsi}_{k+2,k+1}\bm{U}_{k+1} \end{aligned}\right\}\] 或写为 \[\hat{\bm{X}}(k+i/k)=\bm{\varPhi}_{k+i,k}\hat{\bm{X}}(k/k) +\sum_{r=k+1}^{k+i}\bm{\varPsi}_{k+i,r}\bm{U}_{r}\quad(i=1,2,\cdots) \tag{4-6-6}\] 由 (4-6-1) 式进行平差时，因改正数向量均为零，所以，\(\hat{\bm{X}}(k+1/k)\)，\(\hat{\bm{X}}(k+2/k)\)，\(\cdots\) 的误差方差阵即为 \[\left.\begin{aligned} \bm{D}_{\hat{X}}(k+1/k)&=\bm{\varPhi}_{k+1,k}\bm{D}_{\hat{X}}(k/k)\bm{\varPhi}_{k+1,k}^{T} +\bm{\varGamma}_{k+1,k}\bm{D}_{\varOmega}(k)\bm{\varGamma}_{k+1,k}^{T}\\ \bm{D}_{\hat{X}}(k+2/k)&=\bm{\varPhi}_{k+2,k+1}\bm{D}_{\hat{X}}(k+1/k)\bm{\varPhi}_{k+2,k+1}^{T} +\bm{\varGamma}_{k+2,k+1}\bm{D}_{\varOmega}(k+1)\bm{\varGamma}_{k+2,k+1}^{T}\\ &\cdots \end{aligned}\right\}\] 由此可得 \[\begin{aligned} \bm{D}_{\hat{X}}(k+i/k)={}&\bm{\varPhi}_{k+i,k}\bm{D}_{\hat{X}}(k/k)\bm{\varPhi}_{k+i,k}^{T}+{}\\ &\sum_{r=k+1}^{k+i}\bm{\varPhi}_{k+i,r}\bm{\varGamma}_{r,r-1}\bm{D}_{\varOmega}(r-1)\bm{\varGamma}_{r,r-1}^{T}\bm{\varPhi}_{k+i,r}^{T} \end{aligned} \tag{4-6-7}\] 上面的 (4-6-6) 和 (4-6-7) 两式是预测方程的一般形式。下面根据这种一般形式来导出三种类型的预测方程。

(1) 固定区间预测方程

此时是固定观测向量 \(\bm{L}_{1},\bm{L}_{2},\cdots,\bm{L}_{N}\)，即固定 \(k=N\)，而设 \[j=N+i>N=k\quad(i>0)\] 由一般形式 (4-6-6) 和 (4-6-7) 两式可得 \[\left.\begin{aligned} \hat{\bm{X}}(j/N)&=\bm{\varPhi}_{j,N}\hat{\bm{X}}(N/N)+\sum_{r=N+1}^{j}\bm{\varPsi}_{j,r}\bm{U}_{r}\\ \bm{D}_{\hat{X}}(j/N)&=\bm{\varPhi}_{j,N}\bm{D}_{\hat{X}}(N/N)\bm{\varPhi}_{j,N}^{T} +\sum_{r=N+1}^{j}\bm{\varPhi}_{j,r}\bm{\varGamma}_{r,r-1}\bm{D}_{\varOmega}(r-1)\bm{\varGamma}_{r,r-1}^{T}\bm{\varPhi}_{j,r}^{T} \end{aligned}\right\} \tag{4-6-8}\] 式中，\(\bm{\varPhi}_{j,r}=\bm{\varPhi}(t_{j},t_{r})\)。可以由微分方程 \[\left.\begin{aligned} \dot{\bm{\varPhi}}(t,t_{1})&=\bm{A}(t)\bm{\varPhi}(t,t_{1})\\ \bm{\varPhi}(t_{1},t_{1})&=\bm{E} \end{aligned}\right\} \tag{4-6-9}\] 解得 \(\bm{\varPhi}(t,t_{1})\)，然后使 \(r\cdot\Delta t=t_{1}\)，并将 \(\bm{\varPhi}(t,t_{1})\) 离散化即得到 \(\bm{\varPhi}_{j,r}\)。(4-6-9) 式中的 \(\bm{A}(t)\) 是相应的连续系统的系数阵。

(2) 固定点预测方程

固定点预测是根据观测值 \(\bm{L}_{1},\bm{L}_{2},\cdots,\bm{L}_{k}\) 估计将来的某个固定时刻 \(j=N\) 的状态预测值 \(\hat{\bm{X}}(N/k)\)，此时有 \[j=N=k+i>k\] 由 (4-6-5) 和 (4-6-6) 式可得到 \[\left.\begin{aligned} \hat{\bm{X}}(N/k)&=\bm{\varPhi}_{N,k}\hat{\bm{X}}(k/k)+\sum_{r=k+1}^{N}\bm{\varPsi}_{N,r}\bm{U}_{r}\\ \bm{D}_{\hat{X}}(N/k)&=\bm{\varPhi}_{N,k}\bm{D}_{\hat{X}}(k/k)\bm{\varPhi}_{N,k}^{T} +\sum_{r=k+1}^{N}\bm{\varPhi}_{N,r}\bm{\varGamma}_{r,r-1}\bm{D}_{\varOmega}(r-1)\bm{\varGamma}_{r,r-1}^{T}\bm{\varPhi}_{N,r}^{T} \end{aligned}\right\} \tag{4-6-10}\] 根据转移矩阵的性质，可由微分方程 \[\left.\begin{aligned} \frac{\partial\bm{\varPhi}(t_{2},t)}{\partial t}&=-\bm{\varPhi}(t_{2},t)\bm{A}(t)\\ \bm{\varPhi}(t_{2},t_{2})&=\bm{E} \end{aligned}\right\} \tag{4-6-11}\] 解得 \(\bm{\varPhi}(t_{2},t)\)，然后将 \(\bm{\varPhi}(t_{2},t)\) 离散化得到 \(\bm{\varPhi}_{N,k}\)。

(3) 固定超前预测方程

设固定超前量为 \(N\)，则由 \(\bm{L}_{1},\bm{L}_{2},\cdots,\bm{L}_{k}\) 求得的固定超前预测值为 \[\left.\begin{aligned} \hat{\bm{X}}(k+N/k)&=\bm{\varPhi}_{k+N,k}\hat{\bm{X}}(k/k)+\sum_{r=k+1}^{k+N}\bm{\varPsi}_{k+N,r}\bm{U}_{r}\\ \bm{D}_{\hat{X}}(k+N/k)&=\bm{\varPhi}_{k+N,k}\bm{D}_{\hat{X}}(k/k)\bm{\varPhi}_{k+N,k}^{T}+{}\\ &\quad\sum_{r=k+1}^{k+N}\bm{\varPhi}_{k+N,r}\bm{\varGamma}_{r,r-1}\bm{D}_{\varOmega}(r-1)\bm{\varGamma}_{r,r-1}^{T}\bm{\varPhi}_{k+N,r}^{T} \end{aligned}\right\} \tag{4-6-12}\] 为求得转移矩阵 \(\bm{\varPhi}_{k+N,k}\)，可先求 \(\bm{\varPhi}(t+T,t)\) 的微分方程，将 \(\bm{\varPhi}(t+T,t)\) 对时间 \(t\) 求导，可得到 \[\frac{\partial\bm{\varPhi}(t+T,t)}{\partial t} =\bm{A}(t+T)\bm{\varPhi}(t+T,t)-\bm{\varPhi}(t+T,t)\bm{A}(t) \tag{4-6-13}\] 而边界条件 \(\bm{\varPhi}(t_{0}+T,t_{0})\)，可以从 \(t_{0}\) 到 \(t_{0}+T\) 向前解方程 (4-6-8) 式得到，根据边界条件解微分方程 (4-6-13) 式得到 \(\bm{\varPhi}(t+T,t)\)，然后将它离散化就得到 \(\bm{\varPhi}_{k+N,k}\)。

例 4-6-1

设有状态方程和观测方程 \[\begin{aligned} \bm{X}_{k+1}&=\begin{bmatrix}1&1\\ 0&1\end{bmatrix}\bm{X}_{k}+\bm{\varOmega}_{k}\\ L_{k+1}&=\begin{bmatrix}1&0\end{bmatrix}\bm{X}_{k+1}+\varDelta_{k+1} \end{aligned}\] 其中动态噪声 \(\bm{\varOmega}_{k}\) 和观测噪声 \(\varDelta_{k}\) 都是零均值完全不相关的白噪声序列，它们与 \(\bm{X}_{0}\) 相互独立，且 \[\begin{gathered} \bm{D}_{\varOmega}(k)=\sigma_{\varOmega}^{2}\cdot\bm{E},\quad D_{\varDelta}(k)=\sigma_{\varDelta}^{2}\\ E(\bm{X}_{0})=\bm{\mu}_{0},\quad \operatorname{var}(\bm{X}_{0})=\bm{D}_{X}(0) \end{gathered}\] 试求预测值 \(\hat{\bm{X}}(j/k)\,(j>k)\)。

解由状态方程可知 \[\bm{\varPhi}_{j,k}=\begin{bmatrix}1&j-k\\ 0&1\end{bmatrix}\] 预测值的一般形式为 \[\hat{\bm{X}}(j/k)=\begin{bmatrix}1&j-k\\ 0&1\end{bmatrix}\hat{\bm{X}}(k/k)\] 即有 \[\left.\begin{aligned} \hat{X}_{1}(j/k)&=\hat{X}_{1}(k/k)-(j-k)\hat{X}_{2}(k/k)\\ \hat{X}_{2}(j/k)&=\hat{X}_{2}(k/k) \end{aligned}\right\}\] 固定区间预测为 \[\hat{\bm{X}}(j/N)=\begin{bmatrix}1&j-N\\ 0&1\end{bmatrix}\hat{\bm{X}}(N/N)\] 固定点预测为 \[\hat{\bm{X}}(N/k)=\begin{bmatrix}1&N-k\\ 0&1\end{bmatrix}\hat{\bm{X}}(k/k)\] 固定超前预测为 \[\hat{\bm{X}}(k+N/k)=\begin{bmatrix}1&N\\ 0&1\end{bmatrix}\hat{\bm{X}}(k/k)\] 预测值的误差方差阵可由 (4-6-7) 式得 \[\begin{aligned} \bm{D}_{\hat{X}}(j/k)=&\begin{bmatrix}1&j-k\\ 0&1\end{bmatrix}\bm{D}_{\hat{X}}(k/k) \begin{bmatrix}1&0\\ j-k&1\end{bmatrix}+{}\\ &\sum_{r=k+1}^{j}\begin{bmatrix}1&j-r\\ 0&1\end{bmatrix} \begin{bmatrix}\sigma_{\varOmega}^{2}&0\\ 0&\sigma_{\varOmega}^{2}\end{bmatrix} \begin{bmatrix}1&0\\ j-r&1\end{bmatrix} \end{aligned}\] 式中 \[\bm{D}_{\hat{X}}(k/k)=\begin{bmatrix}D_{11}(k/k)&D_{12}(k/k)\\ D_{21}(k/k)&D_{22}(k/k)\end{bmatrix}\] 展开后可得 \[\begin{aligned} D_{11}(j/k)={}&D_{11}(k/k)+2(j-k)D_{12}(k/k)+(j-k)^{2}D_{22}(k/k)+{}\\ &\left[j-k+\sum_{r=k+1}^{j}(j-r)^{2}\,\right]\sigma_{\varOmega}^{2}\\ D_{12}(j/k)={}&D_{12}(k/k)+(j-k)D_{22}(k/k)+\sum_{r=k+1}^{j}(j-r)\,\sigma_{\varOmega}^{2}\\ D_{22}(j/k)={}&D_{22}(k/k)+(j-k)\,\sigma_{\varOmega}^{2} \end{aligned}\]

离散线性系统的平滑

在利用观测向量 \(\bm{L}_{1},\bm{L}_{2},\cdots,\bm{L}_{k}\) 求定第 \(j\) 时刻的状态向量 \(\bm{X}_{j}\) 的最佳估值时，若 \(j<k\)，则称 \(\hat{\bm{X}}(j/k)\) 为最佳平滑值，称求定 \((j/k)\) 和它的误差方差的过程为平滑。根据 \(k\) 和 \(j\) 的变化情况，也可将平滑问题分为三类：

(1) 固定区间平滑：固定 \(k=N\)，变化 \(j\)，即由固定的观测值 \(\bm{L}_{1},\bm{L}_{2},\cdots,\bm{L}_{N}\) 求不同时刻的状态向量的平滑值 \(\hat{\bm{X}}(j/N)\)，\(j<N=k\)；

(2) 固定点平滑：固定 \(j=N\)，变化 \(k\)，即由变化的观测值求某个固定时刻的状态向量的平滑值 \(\hat{\bm{X}}(N/k)\)，\(j=N<k\)；

(3) 固定滞后平滑：\(j\) 和 \(k\) 都变化，但保持 \(j=k-N\)，即由变化的观测值求平滑值 \(\hat{\bm{X}}(k-N/k)\)。

这些平滑问题也是在已知前二阶矩特性的情况下，对一个马尔柯夫或高斯-马尔柯夫序列进行平滑的问题。这里也按广义最小二乘原理来推导 4-3 节的完全不相关的白噪声情形下离散线性系统的平滑方程。

固定区间平滑

根据广义最小二乘原理和 4-3 节所述可知，由观测值 \(\bm{L}_{1},\bm{L}_{2},\cdots,\bm{L}_{N}\) 虚拟观测值 \(\hat{\bm{X}}(0/0)\) 及动态噪声所对应的观测值 \(\bm{L}_{\varepsilon}(i)=-\bm{\varPsi}_{i+1,i}\bm{U}_{i}\,(i=1,2,\cdots,N)\)，可写出误差方程 (4-3-10)；而按逐次平差原理，在利用 \(\bm{L}_{1},\bm{L}_{2},\cdots,\bm{L}_{j}\,(j<k)\) 平差之后，也可写出误差方程 \[\left.\begin{aligned} \bm{V}_{X}(j/j)&=\hat{\bm{X}}(j/N)-\hat{\bm{X}}(j/j)\\ \bm{V}_{\varepsilon}(j)&=\bm{\varPhi}_{j+1,j}\hat{\bm{X}}(j/N)-\hat{\bm{X}}(j+1,N)-\bm{L}_{\varepsilon}(j)\\ \bm{V}_{j+1}&=\bm{B}_{j+1}\hat{\bm{X}}(j+1/N)+\bm{Z}_{j+1}-\bm{L}_{j+1}\\ \bm{V}_{\varepsilon}(j+1)&=\bm{\varPhi}_{j+2,j+1}\hat{\bm{X}}(j+1/N)-\hat{\bm{X}}(j+2/N)-\bm{L}_{\varepsilon}(j+1)\\ &\cdots\\ \bm{V}_{\varepsilon}(N-1)&=\bm{\varPhi}_{N,N-1}\hat{\bm{X}}(N-1/N)-\hat{\bm{X}}(N/N)-\bm{L}_{\varepsilon}(N-1)\\ \bm{V}_{N}&=\bm{B}_{N}\hat{\bm{X}}(N/N)+\bm{Z}_{N}-\bm{L}_{N} \end{aligned}\right\} \tag{4-7-1}\] 对应的方差阵为 \[\bm{D}_{L}=\begin{bmatrix} \bm{D}_{\hat{X}}(j/j) & \bm{0} & \bm{0} & \cdots & \bm{0} & \bm{0}\\ \bm{0} & \bm{D}_{\varepsilon}(j) & \bm{0} & \cdots & \bm{0} & \bm{0}\\ \bm{0} & \bm{0} & \bm{D}_{\Delta}(j+1) & \cdots & \bm{0} & \bm{0}\\ \vdots & \vdots & \vdots & & \vdots & \vdots\\ \bm{0} & \bm{0} & \bm{0} & \cdots & \bm{D}_{\varepsilon}(N-1) & \bm{0}\\ \bm{0} & \bm{0} & \bm{0} & \cdots & \bm{0} & \bm{D}_{\Delta}(N) \end{bmatrix} \tag{4-7-2}\] 其中 \[\bm{D}_{\varepsilon}(i)=\bm{\varGamma}_{i+1,i}\bm{D}_{\varOmega}(i)\bm{\varGamma}_{i+1,i}^{T} \tag{4-7-3}\] 可得到如下的法方程 \[\resizebox{\textwidth}{!}{$ \begin{bmatrix} \bm{N}_{j} & -\bm{\varPhi}_{j+1,j}^{T}\bm{D}_{\varepsilon}^{-1}(j) & \bm{0} & \cdots & \bm{0}\\ -\bm{D}_{\varepsilon}^{-1}(j)\bm{\varPhi}_{j+1,j} & \bm{N}_{j+1} & -\bm{\varPhi}_{j+2,j+1}^{T}\bm{D}_{\varepsilon}^{-1}(j+1) & \cdots & \bm{0}\\ \bm{0} & -\bm{D}_{\varepsilon}^{-1}(j+1)\bm{\varPhi}_{j+2,j+1} & \bm{N}_{j+2} & \cdots & \bm{0}\\ \vdots & \vdots & \vdots & & \vdots\\ \bm{0} & \bm{0} & \bm{0} & \cdots & \bm{N}_{N} \end{bmatrix} \begin{bmatrix} \hat{\bm{X}}(j/N)\\ \hat{\bm{X}}(j+1/N)\\ \hat{\bm{X}}(j+2/N)\\ \vdots\\ \hat{\bm{X}}(N/N) \end{bmatrix} - \begin{bmatrix} \bm{W}_{j}\\ \bm{W}_{j+1}\\ \bm{W}_{j+2}\\ \vdots\\ \bm{W}_{N} \end{bmatrix} =\bm{0}$} \tag{4-7-4}\] 式中 \[\left.\begin{aligned} \bm{N}_{j}&=\bm{D}_{X}^{-1}(j/j)+\bm{\varPhi}_{j+1,j}^{T}\bm{D}_{\varepsilon}^{-1}(j)\bm{\varPhi}_{j+1,j}\\ \bm{W}_{j}&=\bm{D}_{X}^{-1}(j/j)\hat{\bm{X}}(j/j)+\bm{\varPhi}_{j+1,j}^{T}\bm{D}_{\varepsilon}^{-1}(j)\bm{L}_{\varepsilon}(j)\\ \bm{N}_{N}&=\bm{D}_{\varepsilon}^{-1}(N-1)-\bm{B}_{N}^{T}\bm{D}_{\Delta}^{-1}(N)\bm{B}_{N}\\ \bm{W}_{N}&=\bm{D}_{\varepsilon}^{-1}(N-1)\bm{L}_{\varepsilon}(N-1) +\bm{B}_{N}^{T}\bm{D}_{\Delta}^{-1}(N)\,(\bm{L}_{N}-\bm{Z}_{N}) \end{aligned}\right\} \tag{4-7-5}\] \[\left.\begin{aligned} \bm{N}_{i}&=\bm{D}_{\varepsilon}^{-1}(i-1)+\bm{B}_{i}^{T}\bm{D}_{\Delta}^{-1}(i)\bm{B}_{i} +\bm{\varPhi}_{i+1,i}^{T}\bm{D}_{\varepsilon}^{-1}(i)\bm{\varPhi}_{i+1,i}\\ \bm{W}_{i}&=\bm{D}_{\varepsilon}^{-1}(i-1)\bm{L}_{\varepsilon}(i-1) +\bm{B}_{i}^{T}\bm{D}_{\Delta}(i)\,(\bm{L}_{i}-\bm{Z}_{i}) +\bm{\varPhi}_{i+1,i}^{T}\bm{D}_{\varepsilon}(i)\bm{L}_{\varepsilon}(i) \end{aligned}\right\} \tag{4-7-6}\] 在 (4-7-4) 式中，\(\hat{\bm{X}}(N/N)\) 是状态向量 \(\bm{X}_{N}\) 的滤波值，它可以按前面的滤波方程求得；而 \(\hat{\bm{X}}(j/N)\)，\(\hat{\bm{X}}(j+1/N)\)，\(\cdots\)，\(\hat{\bm{X}}(N-1/N)\) 就是 \(\bm{X}_{j}\)，\(\bm{X}_{j+1}\)，\(\cdots\)，\(\bm{X}_{N-1}\) 的平滑值。此式的系数阵是一个拟三对角阵，因此，由它的第一个方程可解得 \[\hat{\bm{X}}(j/N)=\bm{N}_{j}^{-1}\left[\bm{W}_{j} +\bm{\varPhi}_{j+1,j}^{T}\bm{D}_{\varepsilon}^{-1}(j)\hat{\bm{X}}(j+1/N)\right] \tag{4-7-7}\] 再将 (4-7-5) 中的 \(\bm{N}_{j}\)，\(\bm{W}_{j}\) 代入上式得 \[\begin{aligned} \hat{\bm{X}}(j/N)={}&\left[\bm{D}_{X}^{-1}(j/j)-\bm{\varPhi}_{j+1,j}^{T}\bm{D}_{\varepsilon}^{-1}(j)\bm{\varPhi}_{j+1,j}\right]^{-1}\cdot\\ &\{\bm{D}_{X}^{-1}(j/j)\hat{\bm{X}}(j/j) +\bm{\varPhi}_{j+1,j}^{T}\bm{D}_{\varepsilon}^{-1}(j)\left[\hat{\bm{X}}(j+1/N)-\bm{\varPsi}_{j+1,j}\bm{U}_{j}\right]\} \end{aligned} \tag{4-7-8}\] 应用矩阵反演公式可将上式变换为 \[\hat{\bm{X}}(j/N)=\hat{\bm{X}}(j/j)+\bm{\alpha}(j)\, \{\hat{\bm{X}}(j+1/N)-\bm{\varPhi}_{j+1,j}\hat{\bm{X}}(j/j)-\bm{\varPsi}_{j+1,j}\bm{U}_{j}\} \tag{4-7-9}\] 式中 \[\left.\begin{aligned} \bm{\alpha}(j)&=\bm{D}_{\hat{X}}(j/j)\,\bm{\varPhi}_{j+1,j}^{T}\bm{D}_{\hat{X}}^{-1}(j+1/j)\\ \bm{D}_{\hat{X}}(j+1/j)&=\bm{\varPhi}_{j+1,j}\bm{D}_{\hat{X}}(j/j)\bm{\varPhi}_{j+1,j}^{T} +\bm{\varGamma}_{j+1,j}\bm{D}_{\varOmega}(j)\bm{\varGamma}_{j+1,j}^{T} \end{aligned}\right\} \tag{4-7-10}\]

平滑值 \(\hat{\bm{X}}(j/N)\) 的误差方差阵可由法方程 (4-7-4) 求逆得到。因为有 \[\left.\begin{aligned} \bm{N}_{j}\bm{D}_{\hat{X}}(j/N)-\bm{\varPhi}_{j+1,j}^{T}\bm{D}_{\varepsilon}^{-1}(j)\bm{C}_{j,j+1}&=\bm{E}\\ \bm{N}_{j}\bm{C}_{j,j+1}^{T}-\bm{\varPhi}_{j+1,j}^{T}\bm{D}_{\varepsilon}^{-1}(j)\bm{D}_{\hat{X}}(j+1/N)&=\bm{0} \end{aligned}\right\} \tag{4-7-11}\] 式中 \(\bm{C}_{j,j+1}^{T}\) 表示 \(\hat{\bm{X}}(j/N)\) 和 \(\hat{\bm{X}}(j+1/N)\) 的误差间的协方差阵。所以 \[\bm{D}_{\hat{X}}(j/N)=\bm{N}_{j}^{-1}\left[\bm{E} +\bm{\varPhi}_{j+1,j}^{T}\bm{D}_{\varepsilon}^{-1}(j)\bm{D}_{\hat{X}}(j+1/N)\bm{D}_{\varepsilon}^{-1}(j)\bm{\varPhi}_{j+1,j}\bm{N}_{j}^{-1}\right] \tag{4-7-12}\] 应用矩阵反演公式可将上式变换为 \[\bm{D}_{\hat{X}}(j/N)=\bm{D}_{\hat{X}}(j/j) +\bm{\alpha}(j)\left[\bm{D}_{\hat{X}}(j+1/N)-\bm{D}_{\hat{X}}(j+1/j)\right]\bm{\alpha}^{T}(j) \tag{4-7-13}\] 上面得到的 (4-7-9) 和 (4-7-13) 两式就是固定区间平滑的实用公式。当然，(4-7-7) 和 (4-7-12) 式也可以作为固定区间平滑公式。

固定点平滑

若取固定的 \(j=N\)，则由 (4-7-9) 式可得到 \[\begin{aligned} \hat{\bm{X}}(N/k)&=\hat{\bm{X}}(N/N)+\bm{\alpha}(N)\left[\hat{\bm{X}}(N+1/k)-\bm{\varPhi}_{N+1,N}\hat{\bm{X}}(N/N)-\bm{\varPsi}_{N+1,N}\bm{U}_{N}\right]\\ \hat{\bm{X}}(N/k-1)&=\hat{\bm{X}}(N/N)+\bm{\alpha}(N)\left[\hat{\bm{X}}(N+1/k-1)-\bm{\varPhi}_{N+1,N}\hat{\bm{X}}(N/N)-\bm{\varPsi}_{N+1,N}\bm{U}_{N}\right] \end{aligned}\] 两式相减即得 \[\hat{\bm{X}}(N/k)-\hat{\bm{X}}(N/k-1) =\bm{\alpha}(N)\left[\hat{\bm{X}}(N+1/k)-\hat{\bm{X}}(N+1/k-1)\right] \tag{4-7-14}\] 同理，对于 \(i=N+1,\cdots,k-1\)，可得 \[\hat{\bm{X}}(i/k)-\hat{\bm{X}}(i/k-1) =\bm{\alpha}(i)\left[\hat{\bm{X}}(i+1/k)-\hat{\bm{X}}(i+1/k-1)\right]\] 则 (4-7-14) 式可写为 \[\hat{\bm{X}}(N/k)-\hat{\bm{X}}(N/k-1) =\prod_{i=N}^{k-1}\bm{\alpha}(i)\left[\hat{\bm{X}}(k/k)-\hat{\bm{X}}(k/k-1)\right] \tag{4-7-15}\] 于是，可以得到固定点平滑公式 \[\hat{\bm{X}}(N/k)=\hat{\bm{X}}(N/k-1) +\bm{\beta}(N,k)\left[\hat{\bm{X}}(k/k)-\hat{\bm{X}}(k/k-1)\right] \tag{4-7-16}\] 式中 \[\left.\begin{aligned} \bm{\beta}(N,k)&=\prod_{i=N}^{k-1}\bm{\alpha}(i)=\bm{\beta}(N,k-1)\bm{\alpha}(k-1)\\ \bm{\beta}(N,N)&=\bm{E} \end{aligned}\right\} \tag{4-7-17}\] 因为由滤波方程知 \[\hat{\bm{X}}(k/k)=\hat{\bm{X}}(k/k-1) +\bm{J}_{k}\left[\bm{L}_{k}-\bm{Z}_{k}-\bm{B}_{k}\hat{\bm{X}}(k/k-1)\right]\] 故可将 (4-7-16) 式写为 \[\hat{\bm{X}}(N/k)=\hat{\bm{X}}(N/k-1) +\bm{\beta}(N,k)\bm{J}_{k}\left[\bm{L}_{k}-\bm{Z}_{k}-\bm{B}_{k}\hat{\bm{X}}(k/k-1)\right] \tag{4-7-18}\] 根据 (4-7-13) 式可写出 \(\hat{\bm{X}}(N/k)\) 的误差方差阵为 \[\bm{D}_{\hat{X}}(N/k)=\bm{D}_{\hat{X}}(N/N) +\bm{\alpha}(N)\left[\bm{D}_{\hat{X}}(N+1/k)-\bm{D}_{\hat{X}}(N+1/N)\right]\bm{\alpha}^{T}(N)\] 而 \(\hat{\bm{X}}(N/k-1)\) 的误差方差阵为 \[\bm{D}_{\hat{X}}(N/k-1)=\bm{D}_{\hat{X}}(N/N) +\bm{\alpha}(N)\left[\bm{D}_{\hat{X}}(N+1/k-1)-\bm{D}_{\hat{X}}(N+1/N)\right]\bm{\alpha}^{T}(N)\] 两式相减得 \[\bm{D}_{\hat{X}}(N/k)-\bm{D}_{\hat{X}}(N/k-1) =\bm{\alpha}(N)\left[\bm{D}_{\hat{X}}(N+1/k)-\bm{D}_{\hat{X}}(N+1/k-1)\right]\bm{\alpha}^{T}(N)\] 同理，对于 \(i=N+1,\cdots,k-1\) 也可得到类似的结果，因此有 \[\bm{D}_{\hat{X}}(N/k)=\bm{D}_{\hat{X}}(N/k-1) +\bm{\beta}(N,k)\left[\bm{D}_{\hat{X}}(k/k)-\bm{D}_{\hat{X}}(k/k-1)\right]\bm{\beta}^{T}(N,k) \tag{4-7-19}\] 又因为 \[\bm{D}_{\hat{X}}(k/k)=(\bm{E}-\bm{J}_{k}\bm{B}_{k})\,\bm{D}_{\hat{X}}(k/k-1)\] 所以 \[\bm{D}_{\hat{X}}(N/k)=\bm{D}_{\hat{X}}(N/k-1) -\bm{\beta}(N,k)\bm{J}_{k}\bm{B}_{k}\bm{D}_{\hat{X}}(k/k-1)\bm{\beta}^{T}(N,k) \tag{4-7-20}\] (4-7-16) 和 (4-7-19) 各式就是固定点平滑公式。

还可以将 (4-7-18) 和 (4-7-20) 两式改写为 \[\hat{\bm{X}}(N/k)=\bm{X}(N/k-1) +\bar{\bm{\beta}}(N,k)\bm{B}_{k}^{T}\bm{D}_{\Delta}^{-1}(k) \left[\bm{L}_{k}-\bm{Z}_{k}-\bm{B}_{k}\hat{\bm{X}}(k/k-1)\right] \tag{4-7-21}\] \[\begin{aligned} \bm{D}_{\hat{X}}(N/k)={}&\bm{D}_{\hat{X}}(N/k-1) -\bar{\bm{\beta}}(N,k)\left[\bm{B}_{k}^{T}\bm{D}_{\Delta}^{-1}(k)\bm{B}_{k}\bm{D}_{\hat{X}}(k/k-1)\bm{B}_{k}^{T}\bm{D}_{\Delta}(k)\bm{B}_{k}+{}\right.\\ &\left.\bm{B}_{k}^{T}\bm{D}_{\Delta}^{-1}(k)\bm{B}_{k}\right]\bar{\bm{\beta}}(N,k) \end{aligned} \tag{4-7-22}\] 式中 \[\bar{\bm{\beta}}(N,k)=\bar{\bm{\beta}}(N,k-1)\bm{\varPhi}^{T}(k,k-1) \left[\bm{E}-\bm{B}_{k}^{T}\bm{D}_{\Delta}^{-1}(k)\bm{B}_{k}\bm{D}_{\hat{X}}(k/k)\right] \tag{4-7-23}\] 固定点平滑的计算在工程实践中有着广泛的应用。例如，求卫星在 \(\lambda\) 轨点时刻 \(N\) 的状态向量 \(\bm{X}(N)\) 的估值，就是固定点平滑问题。

固定滞后平滑

为了由变化的观测值 \(\bm{L}_{1},\bm{L}_{2},\cdots,\bm{L}_{N+k}\) 求固定滞后平滑值 \(\hat{\bm{X}}(k/N+k)\)，可在固定点平滑公式 (4-7-18) 中以 \(k\) 代替 \(N\)，以 \(N+k\) 代替 \(k\)，得 \[\begin{aligned} \hat{\bm{X}}(k/N+k)={}&\hat{\bm{X}}(k/N+k-1) +\bm{\beta}(k,N+k)\bm{J}_{N+k}\left[\bm{L}_{N+k}-\bm{Z}_{N+k}-{}\right.\\ &\left.\bm{B}_{N+k}\hat{\bm{X}}(N+k/N+k-1)\right] \end{aligned} \tag{4-7-24}\] 式中 \[\left.\begin{aligned} \bm{\beta}(k,N+k)&=\prod_{i=k}^{N+k-1}\bm{\alpha}(i)\\ &=\bm{\alpha}^{-1}(k-1)\bm{\beta}(k-1,N+k-1)\bm{\alpha}(N+k-1)\\ \bm{\beta}(0,N)&=\prod_{i=0}^{N-1}\bm{\alpha}(i) \end{aligned}\right\} \tag{4-7-25}\] 而由固定区间平滑公式 (4-7-9) 知，当 \(\bm{\alpha}(j)\) 为非奇异阵时，有 \[\hat{\bm{X}}(j+1/N)=\bm{\varPhi}_{j+1,j}\hat{\bm{X}}(j/j)+\bm{\varPsi}_{j+1,j}\bm{U}_{j} +\bm{\alpha}^{-1}(j)\left[\hat{\bm{X}}(j/N)-\hat{\bm{X}}(j/j)\right]\] 以 \(k\) 和 \(N+k-1\) 分别代替上式中的 \(j+1\) 和 \(N\)，得 \[\begin{aligned} \hat{\bm{X}}(k/N+k-1)={}&\bm{\varPhi}_{k,k-1}\hat{\bm{X}}(k-1/k-1)+\bm{\varPsi}_{k,k-1}\bm{U}_{k-1}+{}\\ &\bm{\alpha}^{-1}(k-1)\left[\hat{\bm{X}}(k-1/N+k-1)-\hat{\bm{X}}(k-1/k-1)\right] \end{aligned} \tag{4-7-26}\] 或写为 \[\hat{\bm{X}}(k/N+k-1)=\bm{\varPhi}_{k,k-1}\hat{\bm{X}}(k-1/N+k-1)+\bm{\varPsi}_{k,k-1}\bm{U}_{k-1}+{}\] \[\bm{\gamma}(k-1)\left[\hat{\bm{X}}(k-1/N+k-1)-\hat{\bm{X}}(k-1/k-1)\right] \tag{4-7-27}\] 式中 \[\bm{\gamma}(k-1)=\bm{\alpha}^{-1}(k-1)-\bm{\varPhi}_{k,k-1} =\left[\bm{D}_{\hat{X}}(k-1/k-1)\bm{\varPhi}_{k,k-1}^{T}\bm{D}_{\hat{X}}^{-1}(k/k-1)\right]^{-1}-\bm{\varPhi}_{k,k-1} \tag{4-7-28}\] 顾及 \[\bm{D}_{\hat{X}}(k/k-1)=\bm{\varPhi}_{k,k-1}\bm{D}_{\hat{X}}(k-1/k-1)\bm{\varPhi}_{k,k-1}^{T} +\bm{\varGamma}_{k,k-1}\bm{D}_{\varOmega}(k-1)\bm{\varGamma}_{k,k-1}^{T}\] 可知 \[\bm{\gamma}(k-1)=\bm{\varGamma}_{k,k-1}\bm{D}_{\varOmega}(k-1)\bm{\varGamma}_{k,k-1}^{T} \bm{\varPhi}_{k,k-1}^{-T}\bm{D}_{\hat{X}}^{-1}(k-1/k-1) \tag{4-7-29}\] 将 (4-7-26) 式代入 (4-7-24) 式中，即得 \[\begin{aligned} \hat{\bm{X}}(k/N+k)={}&\bm{\varPhi}_{k,k-1}\hat{\bm{X}}(k-1/N+k-1) +\bm{\beta}(k,N+k)\bm{J}_{N+k}\left[\bm{L}_{N+k}-\bm{Z}_{N+k}-{}\right.\\ &\left.\bm{B}_{N+k}\hat{\bm{X}}(N+k/N+k-1)\right] +\bm{\gamma}(k-1)\left[\hat{\bm{X}}(k-1/N+k-1)-{}\right.\\ &\left.\hat{\bm{X}}(k-1/k-1)\right]+\bm{\varPsi}_{k,k-1}\bm{U}_{k-1} \end{aligned} \tag{4-7-30}\] 在求得 \(\hat{\bm{X}}(0/N)\) 之后，就可按上式依次计算固定滞后平滑值 \(\hat{\bm{X}}(1/N+1)\)，\(\hat{\bm{X}}(2/N+2)\)，\(\cdots\)

固定滞后估值 \(\hat{\bm{X}}(k/N+k)\) 的误差方差阵的递推公式也不难导出。由固定区间平滑公式 (4-7-13) 可知 \[\begin{aligned} \bm{D}_{\hat{X}}(k/N+k-1)={}&\bm{D}_{\hat{X}}(k/k-1) +\bm{\alpha}^{-1}(k-1)\left[\bm{D}_{\hat{X}}(k-1/N+k-1)-{}\right.\\ &\left.\bm{D}_{\hat{X}}(k-1/k-1)\right]\bm{\alpha}^{-T}(k-1) \end{aligned} \tag{4-7-31}\] 同时，由固定点平滑公式 (4-7-20) 得 \[\begin{aligned} \bm{D}_{\hat{X}}(k/N+k)={}&\bm{D}_{\hat{X}}(k/N+k-1) -\bm{\beta}(k,N+k)\bm{J}_{N+k}\bm{B}_{N+k}\bm{D}_{\hat{X}}(N+k/N+k-1)\cdot\\ &\bm{\beta}^{T}(k,N+k) \end{aligned} \tag{4-7-32}\] 将 (4-7-31) 式代入 (4-7-32) 式，有 \[\resizebox{\textwidth}{!}{$ \begin{aligned} \bm{D}_{\hat{X}}(k/N+k)={}&\bm{D}_{\hat{X}}(k/k-1) -\bm{\beta}(k,N+k)\bm{J}_{N+k}\bm{B}_{N+k}\bm{D}_{\hat{X}}(N+k/N+k-1)\bm{\beta}^{T}(k,N+k)-{}\\ &\bm{\alpha}^{-1}(k-1)\left[\bm{D}_{\hat{X}}(k-1/k-1)-\bm{D}_{\hat{X}}(k-1/N+k-1)\right]\bm{\alpha}^{-T}(k-1) \end{aligned}$} \tag{4-7-33}\] 在求得 \(\bm{D}_{\hat{X}}(0/N)\) 之后，可由此式计算。

例 4-7-1

题同 [例 4-3-2]，试求：

(1) 在 \(N=4\) 时最佳区间平滑值；

(2) 初始点的平滑值；

(3) 固定滞后 \(N=1\) 情况下的平滑值。

解(1) 因系统的状态方程和观测方程为 \[\left.\begin{aligned} X_{k+1}&=X_{k}+\varOmega_{k}\\ L_{k+1}&=X_{k+1}+\varDelta_{k+1} \end{aligned}\right\}\] 由固定区间平滑公式 (4-7-9) 和 (4-7-13) 可得 \[\hat{X}(j/4)=\hat{X}(j/j)+\alpha(j)\left[\hat{X}(j+1/4)-\hat{X}(j/j)\right]\] \[\alpha(j)=\frac{D_{\hat{X}}(j/j)}{D_{\hat{X}}(j/j)+D_{\varOmega}(j)}=\frac{D(j/j)}{D(j/j)+25}\] \[D_{\hat{X}}(j/4)=D_{\hat{X}}(j/j)+\alpha^{2}(j)\left[D_{\hat{X}}(j+1/4)-D_{\hat{X}}(j+1/j)\right]\] 按以上各式计算所得结果列于表 4-2。表 4-2 中的 \(D_{\hat{X}}(j/j-1)\)、\(J_{j}\)、\(\hat{X}(j/j)\)、\(D_{\hat{X}}(j/j)\) 取自 [例 4-3-2] 的表 4-1。

(2) 由固定点平滑公式 (4-7-21) 式可得 \[\bar{\beta}(0,k)=\bar{\beta}(0,k-1)\left[1-\frac{1}{15}D_{\hat{X}}(k/k)\right]\] \[\hat{X}(0/k)=\hat{X}(0/k-1)+\bar{\beta}(0,k)\cdot\frac{1}{15}\left[L_{k}-\hat{X}(k/k-1)\right]\] \[D_{\hat{X}}(0/k)=D_{\hat{X}}(0/k-1)-\frac{1}{15}\bar{\beta}^{2}(0,k)\left[\frac{1}{15}D_{\hat{X}}(k/k-1)+1\right]\] 按上式计算得初始点的平滑值列于表 4-3。

(3) 由固定滞后平滑公式 (4-7-30) 和 (4-7-32) 式，取 \(N=1\)，有 \[\gamma(k-1)=\frac{25}{D_{\hat{X}}(k-1/k-1)}\] \[\beta(k,1+k)=\alpha(k)\] \[\beta(0,1)=\alpha(0)=0.80\] \[\begin{aligned} \hat{X}(k/k+1)={}&\hat{X}(k-1/k)+\beta(k,1+k)J_{k+1}\left[L_{k+1}-\hat{X}(k+1/k)\right]+{}\\ &\gamma(k-1)\left[\hat{X}(k-1/k)-\hat{X}(k-1/k-1)\right]\\ D_{\hat{X}}(k/k+1)={}&D_{\hat{X}}(k/k-1)-\alpha^{2}(k)J_{k+1}D_{\hat{X}}(k+1/k)-{}\\ &\frac{1}{\alpha^{2}(k-1)}\left[D_{\hat{X}}(k-1/k-1)-D_{\hat{X}}(k-1/k)\right]\\ D_{\hat{X}}(0/1)={}&28.60 \end{aligned}\] 由以上各式计算所得的结果列于表 4-4。

线性确定系统的能观性和能控性

在研究动态线性系统的最佳估计问题时，人们会提出以下问题：

(1) 在一般适当的时间内，由观测值（输出值）\(L(t)\) 能否求得系统的状态 \(\bm{X}(t)\)？

(2) 如果已知系统的状态 \(\bm{X}(t)\)，加入适当的控制项（输入项），能否使系统转移到预期的状态？

为了解决这两个问题，就要引用能观性和能控性的概念。能观性是描述系统的状态与观测量的关系的概念，而能控性则是描述系统的控制项对系统运动控制能力的概念。本节讨论它们的明确定义和判断系统是否具有能观性和能控性的具体方法。

讨论线性确定系统的能观性和能控性，是假定系统在理想状态下，或者说，假定动态噪声和观测噪声为零。

离散系统的能观性

考虑以下没有动态噪声和观测噪声的理想状态下的离散线性系统 \[\left.\begin{aligned} \bm{X}_{k+1}&=\bm{\varPhi}_{k+1,k}\bm{X}_{k}+\bm{\varPsi}_{k+1,k}\bm{U}_{k}\\ \bm{L}_{k+1}&=\bm{B}_{k+1}\bm{X}_{k+1} \end{aligned}\right\} \tag{4-8-1}\] 这一系统具有能观性的定义是：如果存在某个正整数 \(N\)，可以由 \(N\) 次观测 \(\bm{L}_{k}\)，\(\bm{L}_{k+1}\)，\(\cdots\)，\(\bm{L}_{k+N-1}\) 及已知的控制项系列 \(\bm{U}_{k}\)，\(\bm{U}_{k+1}\)，\(\cdots\)，\(\bm{U}_{k+N-1}\) 来唯一确定 \(k\) 时刻的状态 \(\bm{X}_{k}\)，则称系统在 \(k\) 时刻是能观的，如果系统在任何时刻都是能观的，则称系统是一个完全能观的系统。

显然，如果观测向量 \(\underset{n\times 1}{\bm{L}_{k}}\) 与状态向量 \(\underset{m\times 1}{\bm{X}_{k}}\) 的维数相同（即 \(n=m\)），则当 \(\bm{B}_{k}\) 是非奇异阵时，通过同一时刻的一组观测就唯一确定了当时的状态，即 \[\bm{X}_{k}=\bm{B}_{k}^{-1}\bm{L}_{k}\] 在一定情况下，系统 (4-8-1) 式在 \(k\) 时刻具有能观性的充要条件是：存在正整数 \(N\)，使 \[\operatorname{rk}(\bm{W}_{G})=m \tag{4-8-2}\] 其中 \(\bm{W}_{G}\) 是一个 \(N_{n}\times m\) 矩阵，即 \[\bm{W}_{G}=\begin{pmatrix} \bm{B}_{k}\\ \bm{B}_{k+1}\bm{\varPhi}_{k+1,k}\\ \vdots\\ \bm{B}_{k+N-1}\bm{\varPhi}_{k+N-1,k+N} \end{pmatrix} \tag{4-8-3}\] 以上充要条件也等价于： \[\bm{W}_{G}^{T}\bm{W}_{G}=\sum_{j=0}^{N-1}\bm{\varPhi}_{k+j,k}^{T}\bm{B}_{k+j}^{T}\bm{B}_{k+j}\bm{\varPhi}_{k+j,k}>0 \tag{4-8-4}\] 即 \(\bm{W}_{G}^{T}\bm{W}_{G}\) 是对称正定阵。下面来证明之。

因为与 \(N\) 次观测相应的观测方程是 \[\left.\begin{aligned} \bm{L}_{k}&=\bm{B}_{k}\bm{X}_{k}\\ \bm{L}_{k+1}&=\bm{B}_{k+1}\bm{X}_{k+1} \end{aligned}\right\} \tag{4-8-5}\] 由 \(\bm{X}_{k}\) 转移 \(\bm{X}_{i}\,(i=k,k+1,\cdots,k+N-1)\) 的表达式为 \[\bm{X}_{i}=\bm{\varPhi}_{i,k}\bm{X}_{k}+\sum_{j=k+1}^{i}\bm{\varPhi}_{i,j}\bm{\varPsi}_{j,j-1}\bm{U}_{j-1} \tag{4-8-6}\] 将它代入 (4-8-5) 式，得 \[\left.\begin{aligned} \bm{L}_{k}&=\bm{B}_{k}\bm{\varPhi}_{k,k}\bm{X}_{k}\\ \bm{L}_{k+1}&=\bm{B}_{k+1}\bm{\varPhi}_{k+1,k}\bm{X}_{k}+\bm{B}_{k+1}\bm{\varPhi}_{k+1,k+1}\bm{\varPsi}_{k+1,k}\bm{U}_{k}\\ &\cdots\\ \bm{L}_{k+N-1}&=\bm{B}_{k+N-1}\bm{\varPhi}_{k+N-1,k}\bm{X}_{k} +\sum_{j=k+1}^{k+N-1}\bm{B}_{k+N-1}\bm{\varPhi}_{k+N-1,j}\bm{\varPsi}_{j,j-1}\bm{U}_{j-1} \end{aligned}\right\} \tag{4-8-7}\] 式中 \(\bm{\varPhi}_{k,k}=\bm{\varPhi}_{k+1,k+1}=\bm{E}\)，若记 \[\bm{Z}_{i}=\bm{L}_{i}-\sum_{j=i+1}^{i}\bm{B}_{i}\bm{\varPhi}_{i,j}\bm{\varPsi}_{j,j-1}\bm{U}_{j-1} \tag{4-8-8}\] 则 (4-8-7) 式可写为 \[\begin{pmatrix} \bm{B}_{k}\bm{\varPhi}_{k,k}\\ \bm{B}_{k+1}\bm{\varPhi}_{k+1,k}\\ \vdots\\ \bm{B}_{k+N-1}\bm{\varPhi}_{k+N-1,k} \end{pmatrix}\bm{X}_{k} =\begin{pmatrix} \bm{Z}_{k}\\ \bm{Z}_{k+1}\\ \vdots\\ \bm{Z}_{k+N-1} \end{pmatrix} \tag{4-8-9}\] 这是一个线性方程组。因为该方程组中的 \(\bm{X}_{k}\) 有唯一解的充要条件就是 \(\operatorname{rk}(\bm{W}_{G})=m\)，因此，它也就是系统 (4-8-1) 具有能观性的充要条件。

对于定常线性系统 \[\left.\begin{aligned} \bm{X}_{k+1}&=\bm{\varPhi}\bm{X}_{k}+\bm{\varPsi}\bm{U}_{k}\\ \bm{L}_{k+1}&=\bm{B}\bm{X}_{k+1} \end{aligned}\right\} \tag{4-8-10}\] 它具有能观性的充要条件是 \[\operatorname{rk}(\bm{W}_{G})=\operatorname{rk}\begin{pmatrix} \bm{B}\\ \bm{B}\bm{\varPhi}\\ \vdots\\ \bm{B}\bm{\varPhi}^{m-1} \end{pmatrix}=m \tag{4-8-11}\]

例 4-8-1

已知定常线性系统的动态方程和观测方程（理想状态）为 \[\left.\begin{aligned} \bm{X}_{k+1}&=\begin{bmatrix} -a_{3} & -a_{2} & a_{1}\\ a_{4} & 0 & 0\\ 0 & a_{5} & a_{6} \end{bmatrix}\bm{X}_{k}+\begin{bmatrix}1\\ 0\\ 0\end{bmatrix}\bm{U}_{k}\\ \bm{L}_{k+1}&=\begin{bmatrix}0&a_{7}&0\end{bmatrix}\bm{X}_{k+1}\quad(a_{7}\neq 0) \end{aligned}\right\}\] 试判别其能观性。

解因为 \[\bm{W}_{G}=\begin{bmatrix}\bm{B}\\ \bm{B}\bm{\varPhi}\\ \bm{B}\bm{\varPhi}^{2}\end{bmatrix} =\begin{bmatrix} 0 & a_{7} & 0\\ a_{4}a_{7} & 0 & 0\\ -a_{3}a_{4}a_{7} & -a_{2}a_{4}a_{7} & a_{1}a_{4}a_{7} \end{bmatrix}\] 所以 \[\operatorname{rk}(\bm{W}_{G})=\operatorname{rk}\begin{bmatrix} 0 & a_{7} & 0\\ a_{4}a_{7} & 0 & 0\\ 0 & 0 & a_{1}a_{4}a_{7} \end{bmatrix} =\begin{cases} 1, & \text{当}\ a_{1}\neq 0,\ a_{4}=0,\\ 2, & \text{当}\ a_{1}=0,\ a_{4}\neq 0,\\ 3, & \text{当}\ a_{1}\neq 0,\ a_{4}\neq 0. \end{cases}\] 可见，只有当 \(a_{1}\neq 0\) 且 \(a_{4}\neq 0\) 时，该系统才具有能观性。

离散系统的能控性

系统 (4-8-1) 具有能控性的定义是：如果对 \(k\) 时刻的给定状态 \(\bm{X}_{k}=\bm{X}^{*}\)，存在某个正整数 \(N\)，只要选取适当的控制序列 \(\bm{U}_{k},\bm{U}_{k+1},\cdots,\bm{U}_{k+N-1}\)，就可以使 \(\bm{X}_{k+N}\) 达到零状态，则称状态 \(\bm{X}^{*}\) 在 \(k\) 时刻是能控的；如果任意状态 \(\bm{X}_{k}\) 在 \(k\) 时刻都能控，则称系统在 \(k\) 时刻是完全能控的；如果系统在时间定义域上的每一时刻都完全能控，则称该系统是一个完全能控的系统。

系统 (4-8-1) 在 \(k\) 时刻完全能控的充要条件是，存在正整数 \(N\)，使 \[\operatorname{rk}(\bm{W}_{K})=m \tag{4-8-12}\] 式中 \[\bm{W}_{K}=\begin{bmatrix}\bm{\varPsi}_{k+N,k+N-1}&\vdots&\bm{\varPhi}_{k+N,k+N-1}\bm{\varPsi}_{k+N-1,k+N-2}&\vdots&\cdots&\vdots&\bm{\varPhi}_{k+N,k+1}\bm{\varPsi}_{k+1,k}\end{bmatrix} \tag{4-8-13}\] 现证明如下。

因为由 (4-8-1) 式可将 \(\bm{X}_{k+N}\) 表示为 \[\bm{X}_{k+N}=\bm{\varPhi}_{k+N,k}\bm{X}_{k}+\sum_{j=k+1}^{k+N}\bm{\varPhi}_{k+N,j}\bm{\varPsi}_{j,j-1}\bm{U}_{j-1} \tag{4-8-14}\] 根据能控性的定义，系统 (4-8-1) 的完全能控性，相当于对于任何 \(\bm{X}_{k}\)，存在 \(N\)，可以找到 \(\bm{U}_{k}\)，\(\bm{U}_{k+1}\)，\(\cdots\)，\(\bm{U}_{k+N-1}\) 使 \(\bm{X}_{k+N}=0\)。即有 \[\sum_{j=k+1}^{k+N}\bm{\varPhi}_{k+N,j}\bm{\varPsi}_{j,j-1}\bm{U}_{j-1}+\bm{\varPhi}_{k+N}\bm{X}_{k}=0 \tag{4-8-15}\] 显然，对于给定的 \(\bm{X}_{k}\)，序列 \(\bm{U}_{k},\bm{U}_{k+1},\cdots,\bm{U}_{k+N-1}\) 就是 (4-8-15) 这个线性方程组的解。将 (4-8-15) 式写为 \[\resizebox{\textwidth}{!}{$ \begin{bmatrix}\bm{\varPsi}_{k+N,k+N-1}&\vdots&\bm{\varPhi}_{k+N,k+N-1}\bm{\varPsi}_{k+N-1,k+N-2}&\vdots&\cdots&\vdots&\bm{\varPhi}_{k+N,k+1}\bm{\varPsi}_{k+1,k}\end{bmatrix} \begin{bmatrix}\bm{U}_{k+N-1}\\ \bm{U}_{k+N-2}\\ \vdots\\ \bm{U}_{k}\end{bmatrix} +\bm{\varPhi}_{k+N,k}\bm{X}_{k}=0$} \tag{4-8-16}\] 其系数阵就是 \(\bm{W}_{K}\)。因为线性方程组 (4-8-15) 中的 \(\bm{U}_{k},\bm{U}_{k+1},\cdots,\bm{U}_{k+N-1}\) 有唯一解的充要条件就是 \(\operatorname{rk}(\bm{W}_{K})=m\)，且当 \(N\) 充分大时，总可使 \(\bm{W}_{K}\) 的列数等于它的行数 \(m\)。这就证明了 \(\operatorname{rk}(\bm{W}_{K})=m\) 是系统 (4-8-1) 具有能控性的充要条件。

又因为 \(\operatorname{rk}(\bm{W}_{K})=m\) 等价于 \(\bm{W}_{K}\bm{W}_{K}^{T}\) 正定，即等价于 \[\bm{W}_{K}\bm{W}_{K}^{T}=\sum_{j=1}^{N}\bm{\varPhi}_{k+N,k+j}\bm{\varPsi}_{k+j,k+j-1}\bm{\varPsi}_{k+j,k+j-1}^{T}\bm{\varPhi}_{k+N,k+j}^{T}>0 \tag{4-8-17}\] 故有时也以 (4-8-17) 式作为系统 (4-8-1) 具有完全能控性的充要条件。

对于定常系统 (4-8-10)，它具有完全能控性的充要条件是 \[\operatorname{rk}(\bm{W}_{K})=\operatorname{rk}\begin{bmatrix}\bm{\varPsi}&\bm{\varPhi}\bm{\varPsi}&\bm{\varPhi}^{2}\bm{\varPsi}&\cdots&\bm{\varPhi}^{N-1}\bm{\varPsi}\end{bmatrix}=m \tag{4-8-18}\]

例 4-8-2

试判别 [例 4-8-1] 的线性系统的完全能控性。

解因为 \[\bm{W}_{K}=\begin{bmatrix}\bm{\varPsi}&\bm{\varPhi}\bm{\varPsi}&\bm{\varPhi}^{2}\bm{\varPsi}\end{bmatrix} =\begin{bmatrix} 1 & -a_{3} & a_{3}^{2}-a_{2}a_{4}\\ 0 & a_{4} & -a_{3}a_{4}\\ 0 & 0 & a_{4}a_{5} \end{bmatrix}\] \[\operatorname{rk}(\bm{W}_{K})=\operatorname{rk}\begin{bmatrix} 1 & 0 & 0\\ 0 & a_{4} & 0\\ 0 & 0 & a_{4}a_{5} \end{bmatrix} =\begin{cases} 1, & \text{当}\ a_{4}=0,\\ 2, & \text{当}\ a_{5}=0,\ a_{4}\neq 0,\\ 3, & \text{当}\ a_{4}\neq 0,\ a_{5}\neq 0. \end{cases}\] 所以，只有当 \(a_{4}\neq 0\) 且 \(a_{5}\neq 0\) 时，该系统才是完全能控的。结合例 4-8-1 的结论可知，只有当 \(a_{1}\neq 0\)，\(a_{4}\neq 0\) 且 \(a_{5}\neq 0\) 时，该系统才是完全能观和完全能控的系统。

连续线性系统的能观性和能控性

设理想状态下的连续线性系统为 \[\left.\begin{aligned} \dot{\bm{X}}(t)&=\bm{A}(t)\bm{X}(t)+\bm{C}(t)\bm{U}(t)\\ \bm{L}(t)&=\bm{B}(t)\bm{X}(t) \end{aligned}\right\} \tag{4-8-19}\] 这里给出它具有能观性和能控性的定义，并且不加证明地给出这种系统具有能观性和能控性的充要条件。

系统 (4-8-19) 具有能观性的定义是，在有限区间 \([t_{1},t_{2}]\) 内，如果通过观测值 \(\bm{L}(t)\)，能够唯一地确定系统在 \(t_{1}\) 时刻的状态 \(\bm{X}(t_{1})\)，就称该系统在 \(t_{1}\) 时刻是完全能观的。如果在任何时刻 \(t_{1}\) 都是完全能观的，则称该系统是完全能观的。

系统 (4-8-19) 具有完全能观性的充要条件是：对于时刻 \(t_{2}>t_{1}\)，当对称矩阵 \[\bm{W}_{G}(t_{1},t_{2})=\int_{t_{1}}^{t_{2}}\bm{\varPhi}^{T}(t,t_{1})\bm{B}^{T}(t)\bm{B}(t)\,(\bm{\varPhi}t,t_{1})\,\mathrm{d}t \tag{4-8-20}\] 是满秩阵，即当 \[\operatorname{rk}\left[\bm{W}_{G}(t_{1},t_{2})\right]=m \tag{4-8-21}\] 时，系统 (4-8-19) 是完全能观的。

对于定常连续线性系统 \[\left.\begin{aligned} \dot{\bm{X}}(t)&=\bm{A}\bm{X}(t)+\bm{C}\bm{U}(t)\\ \bm{L}(t)&=\bm{B}\bm{X}(t) \end{aligned}\right\} \tag{4-8-22}\] 它具有完全能观性的充要条件是 \[\operatorname{rk}(\bm{W}_{G})=\operatorname{rk}\begin{pmatrix} \bm{B}\\ \bm{B}\bm{A}\\ \vdots\\ \bm{B}\bm{A}^{N-1} \end{pmatrix}=m \tag{4-8-23}\]

系统 (4-8-19) 具有能控性的定义是：对时刻 \(t_{1}\) 的任意状态 \(\bm{X}(t_{1})\)，如果存在控制项 \(\bm{U}(t)\)，在有限的区间 \([t_{1},t_{2}]\) 内，使得 \(\bm{X}(t_{2})=0\)，则称该系统在 \(t_{1}\) 时刻是能控的。如果系统在任何时刻都是能控的，则称系统是完全能控的。

系统 (4-8-19) 完全能控的充要条件是：对于时刻 \(t_{2}>t_{1}\)，当对称阵 \[\bm{W}_{K}(t_{1},t_{2})=\int_{t_{1}}^{t_{2}}\bm{\varPhi}(t_{1},t)\,\bm{C}(t)\,\bm{C}^{T}(t)\,\bm{\varPhi}^{T}(t_{1},t)\,\mathrm{d}t \tag{4-8-24}\] 是满秩阵，即当 \[\operatorname{rk}\left[\bm{W}_{K}(t_{1},t_{2})\right]=m \tag{4-8-25}\] 时，系统 (4-8-19) 是完全能控的。

对于定常连续系统 (4-8-22)，它具有完全能控性的充要条件是 \[\operatorname{rk}(\bm{W}_{K})=\operatorname{rk}\begin{bmatrix}\bm{C}&\bm{A}\bm{C}&\bm{A}^{2}\bm{C}&\cdots&\bm{A}^{N-1}\bm{C}\end{bmatrix}=m \tag{4-8-26}\]

卡尔曼滤波的稳定性

在 4-3 节中推导卡尔曼滤波方程 (4-3-27) 时，曾假定初始状态 \(\bm{X}(0)\) 的均值 \(\bm{\mu}_{X}(0)\) 和方差阵 \(\bm{D}_{X}(0)\) 都是已知的，但在实际问题中，\(\bm{\mu}_{X}(0)\) 和 \(\bm{D}_{X}(0)\) 一般是不准确的，甚至是完全不了解的。那么，所取的这种初始值的不准确性会对滤波的精度带来什么影响呢？本节将作出结论性的说明，如果对系统加上一定的要求，可以使得无论怎样选取初始值，只要递推时间充分长，就能保证其滤波值与最优滤波值任意接近，即在时间充分长之后，初始值的影响可以忽略。这就是滤波的稳定性问题。为叙述方便，我们仅介绍离散型卡尔曼滤波的稳定性。

滤波稳定性的概念

由 4-3 节导出的离散性卡尔曼滤波方程可得 \[\hat{\bm{X}}(k/k)=\left[\bm{E}-\bm{J}_{k}\bm{B}_{k}\right]\bm{\varPhi}_{k,k-1}\hat{\bm{X}}(k-1/k-1) +\bm{J}_{k}(\bm{L}_{k}-\bm{Z}_{k})+\bm{\varPsi}_{k,k-1}\bm{U}_{k-1} \tag{4-9-1}\] 若记 \[\bm{\varPhi}_{KC}(i,j)=(\bm{E}-\bm{J}_{i}\bm{B}_{i})\,\bm{\phi}_{i,j} \tag{4-9-2}\] \[\bm{Y}_{k}=\bm{J}_{k}(\bm{L}_{k}-\bm{Z}_{k})+\bm{\varPsi}_{k,k-1}\bm{U}_{k-1} \tag{4-9-3}\] 则方程 (4-9-1) 可改写为 \[\hat{\bm{X}}(k/k)=\bm{\varPhi}_{KC}(k,k-1)\hat{\bm{X}}(k-1/k-1)+\bm{Y}_{k} \tag{4-9-4}\] 不难看到，卡尔曼滤波方程可以看成为一个线性系统，其中 \(\bm{Y}_{k}\) 就是它的输入项。

一般称 \(\bm{\varPhi}_{KC}(i,j)\) 为系统 (4-9-4) 的状态转移阵，而 \(\bm{\varPhi}_{kC}(k,k-1)\) 是它的一步转移阵。对于 \(\bm{\varPhi}_{KC}(i,j)\) 有 \[\bm{\varPhi}_{KC}(i,j)= \begin{cases} \bm{E}, & i=j\\ \bm{\varPhi}_{KC}(i,i-1)\,\bm{\varPhi}_{KC}(i-1,i-2)\cdots\bm{\varPhi}_{KC}(j+1,j), & i>j\\ \bm{\varPhi}_{KC}^{-1}(j,i), & i<j \end{cases}\] \(\bm{\varPhi}_{KC}(i,j)\) 的可逆性由 (4-9-2) 式可以验证，因为该式中的 \(\bm{\varPhi}_{k,k-1}\) 是可逆阵，且 \((\bm{E}-\bm{J}_{k}\bm{B}_{k})\) 的逆阵也存在，即有 \[(\bm{E}-\bm{J}_{k}\bm{B}_{k})^{-1} =\bm{E}+\bm{D}_{\hat{X}}(k/k-1)\bm{B}_{k}^{T}\bm{D}_{\Delta}^{-1}(k)\bm{B}_{k}\] 因此，讨论滤波方程稳定性问题就归结为讨论方程 (4-9-4) 的稳定性问题。

设初始状态有两个不同的初始值 \(\hat{\bm{X}}_{1}(0)\) 和 \(\hat{\bm{X}}_{2}(0)\)，相应的方差阵为 \(\bm{D}_{\hat{X}_{1}}(0)\) 和 \(\bm{D}_{\hat{X}_{2}}(0)\)，所得到的滤波值分别为 \(\hat{\bm{X}}_{1}(k)\) 和 \(\hat{\bm{X}}_{2}(k)\)，则有以下定义：

(1) 对于任意给定的 \(\varepsilon>0\)，若存在相应的 \(\delta=\delta(\varepsilon)>0\)，当 \(\|\hat{\bm{X}}_{1}(0)-\hat{\bm{X}}_{2}(0)\|<\delta\) 时，使得对一切 \(k=1,2,\cdots\) 均能保证 \[\|\hat{\bm{X}}_{1}(k)-\hat{\bm{X}}_{2}(k)\|<\varepsilon \tag{4-9-5}\] 则称系统 (4-9-4) 是稳定的。其中 \(\|\cdot\|\) 表示向量的范数，如 \(\|\bm{a}\|=\sqrt{\bm{a}^{T}\bm{a}}\)。

(2) 如果系统 (4-9-4) 是稳定的，而且对于任何初始值 \(\hat{\bm{X}}_{1}(0)\) 和 \(\hat{\bm{X}}_{2}(0)\) 均有 \[\lim_{k\rightarrow\infty}\|\hat{\bm{X}}_{1}(k)-\hat{\bm{X}}_{2}(k)\|=0, \tag{4-9-6}\] 则称系统 (4-9-4) 是渐近稳定的。

(3) 如果系统 (4-9-4) 是渐近稳定的，且当输入项 \(\bm{Y}_{k}\) 为有界量时，\(\hat{\bm{X}}(k)\) 必定也是有界量，则称此系统是一致渐近稳定的。

对于稳定的系统来说，只要初始值足够精确，就能保证以后的滤波值达到任意指定的精度。在实际应用中仅仅只有稳定性是不够的，渐近稳定的系统则具有更优的性质，即不论初始值有多大的偏差，只要滤波时间充分长，就能自行消除初始值偏差的影响，保证滤波值的准确性。而且有一致渐近稳定性的系统，在实际应用中才是有效的。

上述稳定性的定义是很直观的，它有以下充分必要条件。

(1) 系统是稳定的充要条件是：存在常数 \(c>0\)，使得所有的 \(k>0\)，均有 \[\|\bm{\varPhi}_{KC}(k,0)\|\leq c \tag{4-9-7}\]

(2) 系统是渐近稳定的充要条件是：当 \(k\rightarrow\infty\) 时，有 \[\lim_{k\rightarrow\infty}\|\bm{\varPhi}_{KC}(k,0)\|=0 \tag{4-9-8}\]

(3) 系统是一致渐近稳定的充要条件是：存在常数 \(c_{1}>0\)，\(c_{2}>0\)，使得对所有的 \(i\geq j>0\)，有 \[\|\bm{\phi}_{KC}(i,j)\|\leq c_{2}e^{-c_{1}(i-j)} \tag{4-9-9}\] 上面的 \(\|\cdot\|\) 表示矩阵的范数，对于矩阵 \(\underset{m\times n}{\bm{A}}\) 有 \[\|\bm{A}\|=\sqrt{\sum_{i=1}^{m}\sum_{j=1}^{n}|a_{ij}|^{2}}\] 由以上充要条件可见，系统 (4-9-4) 的稳定性取决于 \(\bm{\phi}_{KC}(k,k-1)\) 的性质。为了找出什么样的状态方程和观测方程所对应的滤波系统是稳定的，下面讨论随机系统的能观性和能控性。

随机系统的能观性和能控性

在 4-8 节中讨论了确定系统的能观性和能控性的概念，这里仍讨论 4-3 节的随机系统，即 \[\left.\begin{aligned} \bm{X}_{k+1}&=\bm{\varPhi}_{k+1,k}\bm{X}_{k}+\bm{\varGamma}_{k+1,k}\bm{\varOmega}_{k}\\ \bm{L}_{k+1}&=\bm{B}_{k+1}\bm{X}_{k+1}+\bm{\Delta}_{k+1} \end{aligned}\right\} \tag{4-9-10}\] 式中 \(\bm{\varOmega}_{k}\) 和 \(\bm{\Delta}_{k+1}\) 均为完全不相关的零均值白噪声序列，具有 4-3 节中所给出的随机模型。对于这种随机系统，其能观性的定义是：

对于离散随机系统 (4-9-10)，如果对时刻 \(k\)，存在正整数 \(N\)，使得 \[\bm{C}_{G}(k-N+1,k)=\sum_{i=k-N+1}^{k}\bm{\varPhi}_{i,k}^{T}\bm{B}_{i}^{T}\bm{D}_{\Delta}^{-1}(i)\bm{B}_{i}\bm{\varPhi}_{i,k}>0 \tag{4-9-11}\] 则称该系统是完全能观的。而对于该系统，如果存在正整数 \(N\) 和与 \(k\) 无关的正整数 \(\alpha\)，\(\beta\)，使得对所有的 \(k\geq N\)，有 \[0<\alpha\cdot\bm{E}\leq\bm{C}_{G}(k-N+1,k)<\beta\cdot\bm{E} \tag{4-9-12}\] 则称该系统是一致完全能观的。

离散随机系统 (4-9-10) 的能控性定义是：对于离散随机系统 (4-9-10)，如果对时刻 \(k\)，存在正整数 \(N\)，使得 \[\bm{C}_{K}(k-N+1,k)=\sum_{i=k-N+1}^{k}\bm{\varPhi}_{k,i}\bm{\varGamma}_{i,i-1}\bm{D}_{\varOmega}(i-1)\bm{\varGamma}_{i,i-1}^{T}\bm{\varPhi}_{k,i}^{T}>0 \tag{4-9-13}\] 则称该系统是完全能控的。而对于该系统，如果存在正整数 \(N\) 以及与 \(k\) 无关的正数 \(\alpha\)、\(\beta\)，使得所有的 \(k\geq N\)，有 \[0<\alpha\cdot\bm{E}\leq\bm{C}_{K}(k-N+1,k)\leq\beta\cdot\bm{E} \tag{4-9-14}\] 则称该系统是一致能控的。

能观性的概念与状态方程、观测方程有关，从形式上看，随机系统能观性的定义式 (4-9-14) 与确定性系统能观性充要条件 (4-8-4) 相比较，(4-9-14) 式中仅多了 \(\bm{D}_{\Delta}^{-1}(i)\)。如果 \(\bm{C}_{G}(k-N+1,k)>0\)，再令 \[\bm{X}^{*}(k)=\bm{C}_{G}^{-1}(k-N+1,k)\sum_{i=k-N+1}^{k}\bm{\varPhi}_{i,k}^{T}\bm{B}_{i}^{T}\bm{D}_{\Delta}^{-1}(i)\bm{L}_{i}\] 它是 \(N\) 组观测 \(\bm{L}_{k-N+1},\cdots,\bm{L}_{k-1},\bm{L}_{k}\) 的线性函数，其数学期望为 \[\begin{aligned} E[\bm{X}^{*}(k)]&=\bm{C}_{G}^{-1}(k-N+1,k)\sum_{i=k-N+1}^{k}\bm{\varPhi}_{i,k}^{T}\bm{B}_{i}^{T}\bm{D}_{\Delta}^{-1}(i)E(\bm{L}_{i})\\ &=\bm{C}_{G}^{-1}(k-N+1,k)\sum_{i=k-N+1}^{k}\bm{\varPhi}_{i,k}^{T}\bm{B}_{i}^{T}\bm{D}_{\Delta}^{-1}(i)E(\bm{B}_{i},\bm{X}_{i}) \end{aligned}\] 因为 \[\bm{X}_{i}=\bm{\varPhi}_{i,k}\bm{X}_{k}-\sum_{j=i+1}^{k}\bm{\varPhi}_{i,j}\bm{\varGamma}_{j,j-1}\bm{\varOmega}_{j-1}\] 所以 \[E[\bm{X}^{*}(k)]=\bm{C}_{G}^{-1}(k-N+1,k)\sum_{i=k-N+1}^{k}\bm{\varPhi}_{i,k}^{T}\bm{B}_{i}^{T}\bm{D}_{\Delta}^{-1}(i)\bm{B}_{i}\bm{\varPhi}_{i,k}E(\bm{X}_{k})=E(\bm{X}_{k})\] 可见，\(\bm{X}^{*}(k)\) 是状态 \(\bm{X}_{k}\) 基于观测值 \(\bm{L}_{k-N+1},\cdots,\bm{L}_{k-1},\bm{L}_{k}\) 的线性无偏估计。反过来可以验证，如果存在基于观测值 \(\bm{L}_{k-N+1},\cdots,\bm{L}_{k-1},\bm{L}_{k}\)，而不依赖 \(\bm{X}_{k}\) 的验前知识的 \(\bm{X}_{k}\) 的线性无偏估计，则 \(\bm{C}_{k}(k-N+1,k)>0\)。也就是说，离散随机系统的完全能观性说明，基于 \(\bm{L}_{k-N+1},\cdots,\bm{L}_{k-1},\bm{L}_{k}\) 而不依赖于其验前知识的 \(\bm{X}_{k}\) 的线性无偏估计是存在的。

将系统的状态方程改写为 \[\bm{X}_{k}=\bm{\varPhi}_{k,k-N}\bm{X}_{k-N}+\sum_{i=k-N+1}^{k}\bm{\varPhi}_{k,i}\bm{\varGamma}_{i,i-1}\bm{\varOmega}_{i-1}\] 则有 \[\begin{aligned} &E\left[(\bm{X}_{k}-\bm{\varPhi}_{k,k-N}\bm{X}_{k-N})\,(\bm{X}_{k}-\bm{\varPhi}_{k,k-N}\bm{X}_{k-N})^{T}\right]\\ ={}&E\left[\left(\sum_{i=k-N+1}^{k}\bm{\varPhi}_{k,i}\bm{\varGamma}_{i,i-1}\bm{\varOmega}_{i-1}\right) \left(\sum_{i=k-N+1}^{k}\bm{\varPhi}_{k,i}\bm{\varGamma}_{i,i-1}\bm{\varOmega}_{i-1}\right)^{T}\right]\\ ={}&\sum_{i=k-N+1}^{k}\bm{\varPhi}_{k,i}\bm{\varGamma}_{i,i-1}\bm{D}_{\varOmega}(i-1)\bm{\varGamma}_{i,i-1}^{T}\bm{\varPhi}_{k,i}^{T} =\bm{C}_{k}(k-N+1,k) \end{aligned}\] 如果所考虑的动态噪声序列是正态分布的，则 \(\bm{X}_{k}-\bm{\varPhi}_{k,k-N}\bm{X}_{k-N}\) 服从 \(m\) 维正态分布。当系统是完全能控时，其概率密度处处大于零，就说明从任意给定的状态 \(\bm{X}_{k-N}\) 出发，能够以正概率达到 \(k\) 时刻的任何状态 \(\bm{X}_{k}\)。当然，这也是确定性系统能控性概念的推广。

滤波误差方差阵的上、下界

对于具有一致完全能观性和一致完全能控性的线性系统 (4-9-10)，它的滤波误差方差阵 \(\bm{D}_{\hat{X}}(k/k)\) 对于 \(k\) 有一致的有限上界和正定下界。也就是说，可以证明滤波误差方差阵具有以下性质：

(1) 设系统 (4-9-10) 为一致完全能观和一致完全能控，即 (4-9-12) 和 (4-9-14) 两式成立，则对于所有的 \(k\geq N\)，有 \[\bm{D}_{\hat{X}}(k/k)\leq\frac{1+(m\beta)^{2}}{\alpha}\cdot\bm{E} \tag{4-9-15}\]

(2) 在以上条件下，如果初始方差阵 \(\bm{D}_{\hat{X}}(0)\) 是正定的，即 \(\bm{D}_{\hat{X}}(0)>0\)，则对所有的 \(k\geq N\)，有 \[\bm{D}_{\hat{X}}(k/k)\geq\frac{\alpha}{1+(m\beta)^{2}}\cdot\bm{E} \tag{4-9-16}\]

(3) 设系统 (4-9-10) 是一致完全能控的，即式 (4-9-14) 成立，则对于所有的 \(k\geq N\)，恒有 \[\bm{D}_{\hat{X}}(k/k)>0 \tag{4-9-17}\] 性质 (3) 中去掉了性质 (2) 中的 \(\bm{D}_{\hat{X}}(0)>0\) 的限制。

将以上性质归结起来，可以得出结论：对于一致完全能观和一致完全能控的线性系统，它的滤波误差方差阵 \(\bm{D}_{\hat{X}}(k/k)\)，当 \(k\geq N\) 时，恒满足不等式 \[\frac{\alpha}{1+(m\beta)^{2}}\cdot\bm{E}\leq\bm{D}_{\hat{X}}(k/k)<\frac{1+(m\beta)^{2}}{\alpha}\cdot\bm{E} \tag{4-9-18}\] 其中 \(N\) 为正整数，\(\alpha>0\)，\(\beta>0\)，由 (4-9-12) 和 (4-9-14) 两式给出。

滤波误差方差阵的上、下界的性质是这类线性系统的基本特性，也是证明滤波稳定性要用到的重要性质。

稳定性定理

这里只结论性地给出两个滤波稳定性定理。

定理一：设系统 (4-9-10) 是一致完全能观和一致完全能控的，则它的最优线性滤波系统 (4-9-4) 是一致渐近稳定的，即存在常数 \(c_{1}>0\) 和 \(c_{2}>0\)，使得对一切 \(i\geq j\geq 0\)，有 \[\|\bm{\varPhi}_{KC}(i,j)\|\leq c_{2}e^{-c_{1}(i-j)}\] 这个定理表明，对于一致完全能观和一致完全能控的线性系统，当时间充分长以后，它的最优线性滤波值将渐近地不依赖于滤波初值的选取，而且有界的观测输入将得到有界的滤波输出。

定理二：在定理一的条件下，如果 \(\bm{D}_{\hat{X}}^{(1)}(0)\) 和 \(\bm{D}_{\hat{X}}^{(2)}(0)\) 是两个不同的初始误差方差阵，则对于 \(k\) 时刻相应的滤波误差方差阵 \(\bm{D}_{\hat{X}}^{(1)}(k/k)\) 和 \(\bm{D}_{\hat{X}}^{(2)}(k/k)\)，存在常数 \(c_{3}>0\) 和 \(c_{4}>0\)，使得对一切 \(k\geq j>0\)，有 \[\|\bm{D}_{\hat{X}}^{(2)}(k/k)-\bm{D}_{\hat{X}}^{(1)}(k/k)\| \leq c_{4}e^{-c_{3}(k-j)}\,\|\bm{D}_{\hat{X}}^{(2)}(j/j)-\bm{D}_{\hat{X}}^{(1)}(j/j)\|\] 上式表明，对于一致完全能观和一致完全能控的线性系统，当 \(k\rightarrow\infty\) 时，有 \[\|\bm{D}_{\hat{X}}^{(2)}(k/k)-\bm{D}_{\hat{X}}^{(1)}(k/k)\| \leq c_{4}e^{-c_{3}k}\,\|\bm{D}_{\hat{X}}^{(2)}(0)-\bm{D}_{\hat{X}}^{(1)}(0)\|\rightarrow 0\] 即当时间充分长以后，它的滤波误差方差阵也将渐近地不依赖于初始方差阵的选取。

定常系统的情形

如果线性系统 (4-9-10) 是定常系统，即设系统为 \[\left.\begin{aligned} \bm{X}_{k+1}&=\bm{\varPhi}\bm{X}_{k}+\bm{\varGamma}\bm{\varOmega}_{k}\\ \bm{L}_{k+1}&=\bm{B}\bm{X}_{k+1}+\bm{\Delta}_{k} \end{aligned}\right\} \tag{4-9-19}\] 且设 \[\begin{gathered} \bm{D}_{\varOmega}(k)=\bm{D}_{\varOmega}>0\\ \bm{D}_{\Delta}(k)=\bm{D}_{\Delta}>0 \end{gathered}\] 容易看出，此时对于 \(k\geq N\)，有 \[\begin{aligned} \bm{C}_{G}(k-N+1,k)&=\sum_{j=k-N+1}^{k}(\bm{\varPhi}^{j-k})^{T}\bm{B}^{T}\bm{D}_{\Delta}^{-1}\bm{B}\bm{\varPhi}^{j-k}\\ &=(\bm{\varPhi}^{-N+1})^{T}\sum_{i=0}^{N-1}(\bm{\varPhi}^{i})^{T}\bm{B}^{T}\bm{D}_{\Delta}^{-1}\bm{B}\bm{\varPhi}^{i}\,(\bm{\varPhi}^{-N+1}) \end{aligned}\] \[\begin{aligned} \bm{C}_{K}(k-N+1,k)&=\sum_{i=k-N+1}^{k}(\bm{\varPhi}^{k-i})\,\bm{\varGamma}\bm{D}_{\varOmega}\bm{\varGamma}^{T}(\bm{\varPhi}^{k-i})^{T}\\ &=\sum_{i=0}^{N-1}(\bm{\varPhi}^{i})\,\bm{\varGamma}\bm{D}_{\varOmega}\bm{\varGamma}^{T}(\bm{\varPhi}^{i})^{T} \end{aligned} \tag{4-9-20}\] 可以证明，定常系统 (4-9-19) 的一致完全能观和一致完全能控的充分必要条件分别是 \[\sum_{i=0}^{m-1}(\bm{\varPhi}^{i})^{T}\bm{B}^{T}\bm{B}(\bm{\varPhi}^{i})>0 \tag{4-9-21}\] \[\sum_{i=0}^{m-1}(\bm{\varPhi}^{i})\,\bm{\varGamma}\bm{\varGamma}^{T}(\bm{\varPhi}^{i})^{T}>0 \tag{4-9-22}\] 式中 \(m\) 是状态的维数。

由此可知，对定常系统而言，一致完全能观性和一致完全能控性，分别与完全能观性、完全能控性是一样的。

定常系统的渐近性质还有以下定理：

定理三：对于完全能观和完全能控的定常线性系统 (4-9-19)，存在唯一的正定阵 \(\bm{D}\) 和定常矩阵 \(\bm{J}\)，使得从任意的初始方差阵 \(\bm{D}_{\hat{X}}(0)\) 出发，当 \(k\rightarrow\infty\) 时，恒有 \[\left.\begin{aligned} \lim_{k\rightarrow\infty}\bm{D}_{\hat{X}}(k/k)&=\bm{D}\\ \lim_{k\rightarrow\infty}\bm{J}_{k}&=\bm{J} \end{aligned}\right\} \tag{4-9-23}\]

由这个定理可知，对于完全能观和完全能控的定常线性系统，不论如何选取滤波初值，当时间充分长以后，它的滤波误差方差将趋于一个唯一确定的正定矩阵，从而它的增益阵也将趋于一个唯一确定的增益矩阵。这称为滤波达到稳态。定理还说明，一个渐近稳定的滤波系统，其估计误差方差既不会无限增大，也不会无限制地精确到零。所以，我们对于滤波的精度不可要求过高，也不必担心失之过大。

模型误差分析

在讨论卡尔曼滤波的性质时，一般都认为系统的数学模型是完全正确的。但实际工作中将具体问题模型化时，往往由于对问题的本质认识不足，或由于为了简化计算而包含某些近似，使得所采用的模型与实际情况存在差别，这种差别就称为模型误差。显然，模型误差对任何具体问题都是存在的，因此，分析模型误差对滤波所产生的影响是非常必要的。

设有实际的离散系统为 \[\left.\begin{aligned} \bm{X}_{k+1}&=\bm{\varPhi}_{k+1,k}\bm{X}_{k}+\bm{U}_{k}+\bm{\varGamma}_{k+1,k}\bm{\varOmega}_{k}\\ \bm{L}_{k+1}&=\bm{B}_{k+1}\bm{X}_{k+1}+\bm{\Delta}_{k+1} \end{aligned}\right\} \tag{4-10-1}\] 其中 \[\left.\begin{aligned} &E(\bm{\varOmega}_{k})=\bm{0},\quad E(\bm{\Delta}_{k})=\bm{0},\quad E(\bm{X}_{0})=\bm{\mu}_{X}(0),\\ &\operatorname{cov}(\bm{\varOmega}_{k},\bm{\varOmega}_{j})=\bm{D}_{\varOmega}(k)\delta_{kj},\quad \operatorname{cov}(\bm{\Delta}_{k},\bm{\Delta}_{j})=\bm{D}_{\Delta}(k)\delta_{kj}\\ &\operatorname{var}(\bm{X}_{0})=\bm{D}_{X}(0),\quad \operatorname{cov}(\bm{\varOmega}_{k},\bm{\Delta}_{j})=\bm{0},\\ &\operatorname{cov}(\bm{X}_{0},\bm{\varOmega}_{k})=\bm{0},\quad \operatorname{cov}(\bm{X}_{0},\bm{\Delta}_{k})=\bm{0} \end{aligned}\right\} \tag{4-10-2}\] 但在滤波时，将它模型化为 \[\left.\begin{aligned} \overline{\bm{X}}_{k+1}&=\overline{\bm{\varPhi}}_{k+1,k}\overline{\bm{X}}_{k}+\overline{\bm{U}}_{k}+\overline{\bm{\varGamma}}_{k+1,k}\overline{\bm{\varOmega}}_{k}\\ \bm{L}_{k+1}&=\overline{\bm{B}}_{k+1}\overline{\bm{X}}_{k+1}+\overline{\bm{\Delta}}_{k+1} \end{aligned}\right\} \tag{4-10-3}\] 而其中 \[\left.\begin{aligned} &E(\overline{\bm{\varOmega}}_{k})=\bm{0},\quad E(\overline{\bm{\Delta}}_{k})=\bm{0},\ E(\overline{\bm{X}}_{0})=\overline{\bm{\mu}_{X}}(0)\\ &\operatorname{cov}(\overline{\bm{\varOmega}}_{k},\overline{\bm{\varOmega}}_{j})=\overline{\bm{D}}_{\varOmega}(k)\delta_{kj},\quad \operatorname{cov}(\overline{\bm{\Delta}}_{k},\overline{\bm{\Delta}}_{j})=\overline{\bm{D}}_{\Delta}(k)\delta_{kj}\\ &\operatorname{var}(\overline{\bm{X}}_{0})=\overline{\bm{D}}_{X}(0),\quad \operatorname{cov}(\overline{\bm{\varOmega}}_{k},\overline{\bm{\Delta}}_{j})=\bm{0}\\ &\operatorname{cov}(\overline{\bm{X}}_{0},\overline{\bm{\Delta}}_{k})=\bm{0},\quad \operatorname{cov}(\overline{\bm{X}}_{0},\overline{\bm{\varOmega}}_{k})=\bm{0} \end{aligned}\right\} \tag{4-10-4}\] 当采用模型 (4-10-3) 进行卡尔曼滤波时，其相应的滤波方程是 \[\hat{\overline{\bm{X}}}(k/k)=\hat{\overline{\bm{X}}}(k/k-1) +\overline{\bm{J}}_{k}\left[\bm{L}_{k}-\overline{\bm{B}}_{k}\hat{\overline{\bm{X}}}(k/k-1)\right] \tag{4-10-5}\] \[\overline{\bm{D}}_{\hat{X}}(k/k)=(\bm{E}-\overline{\bm{J}}_{k}\overline{\bm{B}}_{k})\,\overline{\bm{D}}_{\hat{X}}(k/k-1) \tag{4-10-6}\] \[\left.\begin{aligned} \hat{\overline{\bm{X}}}(k/k-1)&=\overline{\bm{\varPhi}}_{k,k-1}\hat{\overline{\bm{X}}}(k-1/k-1)+\overline{\bm{U}}_{k-1}\\ \overline{\bm{D}}_{\hat{X}}(k/k-1)&=\overline{\bm{\varPhi}}_{k,k-1}\overline{\bm{D}}_{\hat{X}}(k-1/k-1)\overline{\bm{\varPhi}}_{k,k-1}^{T} +\overline{\bm{\varGamma}}_{k,k-1}\overline{\bm{D}}_{\varOmega}(k-1)\overline{\bm{\varGamma}}_{k,k-1}^{T}\\ \overline{\bm{J}}_{k}&=\overline{\bm{D}}_{\hat{X}}(k/k-1)\overline{\bm{B}}_{k}^{T} \left[\overline{\bm{B}}_{k}\bm{D}_{\hat{X}}(k/k-1)\overline{\bm{B}}_{k}^{T}+\overline{\bm{D}}_{\Delta}(k)\right]^{-1} \end{aligned}\right\} \tag{4-10-7}\] 或 \[\hat{\overline{\bm{X}}}(k/k)=\overline{\bm{D}}_{\hat{X}}(k/k) \left[\overline{\bm{D}}_{\hat{X}}^{-1}(k/k-1)\hat{\overline{\bm{X}}}(k/k-1) +\overline{\bm{B}}_{k}^{T}\overline{\bm{D}}_{\Delta}^{-1}(k)\bm{L}_{k}\right] \tag{4-10-8}\] \[\overline{\bm{D}}_{\hat{X}}(k/k)=\left[\overline{\bm{D}}_{\hat{X}}^{-1}(k/k-1) +\overline{\bm{B}}_{k}^{T}\overline{\bm{D}}_{\Delta}^{-1}(k)\overline{\bm{B}}_{k}\right]^{-1} \tag{4-10-9}\] 不难理解，按 (4-10-5) 或 (4-10-8) 式求得的 \(\hat{\overline{\bm{X}}}(k/k)\) 不是真实的最优滤波值，而由 (4-10-6) 或 (4-10-9) 式求得的 \(\overline{\bm{D}}_{\hat{X}}(k/k)\) 即不是 \(\bm{X}_{k}\) 的真实滤波值的误差方差阵，也不是 \(\hat{\overline{\bm{X}}}(k/k)\) 的实际均方误差阵。当然，可以导出计算 \(\hat{\overline{\bm{X}}}(k/k)\) 的实际均方误差阵的递推公式，但这些公式实际上是不能用于具体计算的，因为问题的真正模型并不知道。

这里考虑两种较简单的情况。

\(\bm{D}_{X}(0)\)、\(\bm{D}_{\varOmega}(k)\) 和 \(\bm{D}_{\Delta}(k)\) 有误差的情况

当初始状态的方差阵 \(\bm{D}_{X}(0)\)、动态噪声和观测噪声的方差阵 \(\bm{D}_{\varOmega}(k)\) \(\bm{D}_{\Delta}(k)\) 有误差时，(4-10-5)   (4-10-7) 式可写为 \[\hat{\overline{\bm{X}}}(k/k)=\hat{\overline{\bm{X}}}(k/k-1) +\overline{\bm{J}}_{k}\left[\bm{L}_{k}-\bm{B}_{k}\hat{\overline{\bm{X}}}(k/k-1)\right] \tag{4-10-10}\] \[\overline{\bm{D}}_{\hat{X}}(k/k)=(\bm{E}-\overline{\bm{J}}_{k}\bm{B}_{k})\,\overline{\bm{D}}_{\hat{X}}(k/k-1) \tag{4-10-11}\] \[\left.\begin{aligned} \hat{\overline{\bm{X}}}(k/k-1)&=\bm{\varPhi}_{k,k-1}\hat{\overline{\bm{X}}}(k-1/k-1)+\bm{U}_{k}\\ \overline{\bm{D}}_{\hat{X}}(k/k-1)&=\bm{\varPhi}_{k,k-1}\overline{\bm{D}}_{\hat{X}}(k-1/k-1) +\bm{\varGamma}_{k,k-1}\overline{\bm{D}}_{\varOmega}(k-1)\bm{\varGamma}_{k,k-1}^{T}\\ \overline{\bm{J}}_{k}&=\overline{\bm{D}}_{\hat{X}}(k/k-1)\bm{B}_{k}^{T} \left[\bm{B}_{k}\overline{\bm{D}}_{\hat{X}}(k/k-1)\bm{B}_{k}^{T}+\overline{\bm{D}}_{\Delta}(k)\right]^{-1} \end{aligned}\right\} \tag{4-10-12}\] 而 \(\hat{\overline{\bm{X}}}(k/k)\) 和 \(\hat{\overline{\bm{X}}}(k/k-1)\) 的实际均方误差阵应为 \[\left.\begin{aligned} \hat{\overline{\bm{D}}}_{\hat{X}}(k/k)&=(\bm{E}-\overline{\bm{J}}_{k}\bm{B}_{k})\,\hat{\overline{\bm{D}}}_{\hat{X}}(k/k-1)\\ \hat{\overline{\bm{D}}}_{\hat{X}}(k/k-1)&=\bm{\varPhi}_{k,k-1}\hat{\overline{\bm{D}}}_{\hat{X}}(k-1/k-1)\bm{\varPhi}_{k,k-1}^{T} +\bm{\varGamma}_{k,k-1}\bm{D}_{\varOmega}(k-1)\bm{\varGamma}_{k,k-1}^{T} \end{aligned}\right\} \tag{4-10-13}\] 现令 \[\left.\begin{aligned} \Delta\bm{D}_{\hat{X}}(k/k)&=\hat{\overline{\bm{D}}}_{\hat{X}}(k/k)-\overline{\bm{D}}_{\hat{X}}(k/k)\\ \Delta\bm{D}_{\hat{X}}(k/k-1)&=\hat{\overline{\bm{D}}}_{\hat{X}}(k/k-1)-\overline{\bm{D}}_{\hat{X}}(k/k-1) \end{aligned}\right\} \tag{4-10-14}\] 则由 (4-10-10)   (4-10-13) 式可得 \[\Delta\bm{D}_{\hat{X}}(k/k)=(\bm{E}-\overline{\bm{J}}_{k}\bm{B}_{k})\,\Delta\bm{D}_{\hat{X}}(k/k-1)\,(\bm{E}-\overline{\bm{J}}_{k}\bm{B}_{k})^{T} \tag{4-10-15}\] \[\begin{aligned} \Delta\bm{D}_{\hat{X}}(k/k-1)={}&\bm{\varPhi}_{k,k-1}\Delta\bm{D}_{\hat{X}}(k-1/k-1)\bm{\varPhi}_{k,k-1}^{T}+{}\\ &\bm{\varGamma}_{k,k-1}\left[\bm{D}_{\varOmega}(k-1)-\overline{\bm{D}}_{\varOmega}(k-1)\right]\bm{\varGamma}_{k,k-1}^{T} \end{aligned} \tag{4-10-16}\]

由上两式可见，如果在选取 \(\overline{\bm{D}}_{\varOmega}(k)\) 和 \(\overline{\bm{D}}_{\Delta}(k)\) 时，采用保险的办法，使得 \[\begin{gathered} \overline{\bm{D}}_{\varOmega}(k)\geq\bm{D}_{\varOmega}(k)\\ \overline{\bm{D}}_{\Delta}(k)\geq\bm{D}_{\Delta}(k) \end{gathered}\] 则对于一切 \(k>0\)，有 \(\overline{\bm{D}}_{\hat{X}}(k/k)\geq\hat{\overline{\bm{D}}}_{\hat{X}}(k/k)\)。

可以看到，对于某些实际问题，虽然不能确定真实的 \(\bm{D}_{X}(0)\)、\(\bm{D}_{\varOmega}(k)\) 和 \(\bm{D}_{\Delta}(k)\)，但在滤波时，可以采用它们的保险值。这样，滤波值的真实均方误差阵必然小于计算得到的误差方差阵。

此外，若还设系统 (4-10-3) 是一致完全可观和一致完全可控的，则由 4-9 节可知，存在常数 \(r>0\)，使得对所有的 \(k\)，有 \[\overline{\bm{D}}_{\hat{X}}(k/k)\leq r\cdot\bm{E}\] 所以，对于我们所讨论的这种情况有 \[\hat{\overline{\bm{D}}}_{\hat{X}}(k/k)\leq\overline{\bm{D}}_{\hat{X}}(k/k)\leq r\cdot\bm{E} \tag{4-10-17}\] 即 \(\hat{\overline{\bm{D}}}_{\hat{X}}(k/k)\) 有一致的上界。

\(\bm{D}_{X}(0)\)、\(\bm{D}_{\varOmega}(k)\)、\(\bm{D}_{\Delta}(k)\) 和 \(\bm{U}_{k}\) 存在误差的情况

非随机项 \(\overline{\bm{U}}_{k}\) 的不准确，一般是由于忽略模型的非线性等原因所产生的。此时，由 (4-10-5) 式可得 \(\hat{\overline{\bm{X}}}(k/k)\) 和 \(\hat{\overline{\bm{X}}}(k/k-1)\) 的估计误差为 \[\widetilde{\overline{\bm{X}}}(k/k)=\bm{X}_{k}-\hat{\overline{\bm{X}}}(k/k) =(\bm{E}-\overline{\bm{J}}_{k}\bm{B}_{k}^{T})\,\widetilde{\overline{\bm{X}}}(k/k-1)-\overline{\bm{J}}_{k}\Delta k \tag{4-10-18}\] \[\begin{aligned} \widetilde{\overline{\bm{X}}}(k/k-1)&=\bm{X}_{k}-\hat{\overline{\bm{X}}}(k/k-1) =\bm{\varPhi}_{k,k-1}\widetilde{\overline{\bm{X}}}(k-1/k-1)+{}\\ &\quad\Delta\bm{U}_{k-1}+\bm{\varGamma}_{k,k-1}\bm{\varOmega}_{k-1} \end{aligned} \tag{4-10-19}\] 式中 \(\Delta\bm{U}_{k-1}=\bm{U}_{k-1}-\overline{\bm{U}}_{k-1}\)。相应的均方误差阵 (4-10-6) 为 \[\hat{\overline{\bm{D}}}_{\hat{X}}(k/k)=(\bm{E}-\overline{\bm{J}}_{k}\bm{B}_{k})\,\hat{\overline{\bm{D}}}_{\hat{X}}(k/k-1)\,(\bm{E}-\overline{\bm{J}}_{k}\bm{B}_{k})^{T} +\overline{\bm{J}}_{k}\bm{D}_{\Delta}(k)\,\overline{\bm{J}}_{k}^{T} \tag{4-10-20}\] \[\hat{\overline{\bm{D}}}_{\hat{X}}(k/k-1)=\bm{\varPhi}_{k,k-1}\hat{\overline{\bm{D}}}_{\hat{X}}(k-1/k-1)\bm{\varPhi}_{k,k-1}^{T} +\bm{\varGamma}_{k,k-1}\bm{D}_{\varOmega}(k-1)\bm{\varGamma}_{k,k-1}^{T}+\Delta\bm{U}_{k-1}\Delta\bm{U}_{k-1}^{T} \tag{4-10-21}\] 假定系统 (4-10-3) 是一致完全能观和一致完全能控的，且设 \[\bm{F}_{k}=(\bm{E}-\overline{\bm{J}}_{k}\bm{B}_{k}) \left[\bm{\varGamma}_{k,k-1}\bm{D}_{\varOmega}(k-1)\bm{\varGamma}_{k,k-1}^{T}+\Delta\bm{U}_{k-1}\Delta\bm{U}_{k-1}^{T}\right] +\overline{\bm{J}}_{k}\bm{D}_{\Delta}(k)\,\overline{\bm{J}_{k}^{T}}\] 有一致的上界，则由 4-9 节的稳定性定理已知，对于 \[\overline{\bm{\varPhi}}_{KC}(k,k-1)=(\bm{E}-\overline{\bm{J}}_{k}\bm{B}_{k})\,\bm{\phi}_{k,k-1}\] 存在常数 \(c_{1}>0\) 和 \(_{2}>0\)，使得一切 \(i\geq j\geq 0\)，有 \[\|\overline{\bm{\varPhi}}_{KC}(i,j)\|\leq c_{2}e^{-c_{1}(i-j)} \tag{4-10-22}\] 将 (4-10-20) 式写为 \[\begin{aligned} \hat{\overline{\bm{D}}}_{\hat{X}}(k/k)&=\overline{\bm{\varPhi}}_{KC}(k,k-1)\,\hat{\overline{\bm{D}}}_{\hat{X}}(k-1/k-1)\,\overline{\bm{\varPhi}}_{KC}^{T}(k,k-1)+\bm{F}_{k}\\ &=\overline{\bm{\varPhi}}_{KC}(k,0)\,\bm{D}_{\hat{X}}(0)\,\overline{\bm{\varPhi}}_{KC}(k,0) +\sum_{i=1}^{k}\overline{\bm{\varPhi}}(k,i)\bm{F}_{i} \end{aligned} \tag{4-10-23}\] 设存在常数 \(c>0\)，对一切 \(k\)，\(\|\bm{F}_{k}\|\leq c\)，则由上式可得 \[\begin{aligned} \hat{\overline{\bm{D}}}_{\hat{X}}(k/k)&\leq\|\overline{\bm{\varPhi}}_{KC}(k,0)\|^{2}\,\|\bm{D}_{\hat{X}}(0)\| +\sum_{i=1}^{k}\|\bm{\varPhi}_{KC}(k,i)\|^{2}\,\|\bm{F}_{i}\|\\ &\leq\|\bm{D}_{\hat{X}}(0)\|\,c_{2}^{2}e^{2c_{1}k}+cc_{2}\sum_{i=1}^{k}e^{-2c_{1}(k-i)}\\ &\leq\|\bm{D}_{\hat{X}}(0)\|\,c_{2}^{2}+cc_{2}^{2}/(1-e^{-2c_{1}}) \end{aligned}\] 于是可得结论：对于 \(\bm{D}_{\hat{X}}(0)\)、\(\bm{D}_{\varOmega}(k)\)、\(\bm{D}_{\Delta}(k)\) 和 \(\bm{U}_{k}\) 有误差的情况，若系统 (4-10-3) 是一致完全能观和一致完全能控的，且 \(\bm{F}_{k}\) 有一致的上界，则 \(\hat{\overline{\bm{D}}}_{\hat{X}}(k/k)\) 也有一致的上界。

还可以证明，对于所考虑的情况，如果系统 (4-10-3) 是一致完全能观和一致完全能控的，\(\|\Delta\bm{U}_{k}\|\) 有一致的上界，且存在常数 \(c_{4}>0\) 和 \(c_{5}>0\)，使得对一切 \(k\)，有 \[\bm{D}_{\varOmega}(k)<c_{4}\overline{\bm{D}}_{\varOmega}(k),\qquad \bm{D}_{\Delta}(k)<c_{5}\overline{\bm{D}}_{\Delta}(k)\] 则滤波的真正均方误差阵 \(\hat{\overline{\bm{D}}}_{\hat{X}}(k/k)\) 必有一致的上界。

由此可见，在滤波时，为了防止由于模型不准确而使 \(\hat{\overline{\bm{D}}}_{\hat{X}}(k/k)\) 无上界，即防止出现发散现象，可以根据这个结论来适当选择 \(\overline{\bm{D}}_{\varOmega}(k)\) 和 \(\overline{\bm{D}}_{\Delta}(k)\)，并使模型化的系统具有一致完全能观性和一致完全能控性，这是一个很基本的手段。

滤波的发散现象和克服发散的方法

从理论上讲，随着观测数据的增加，滤波应该得到更精确的状态估计。但是，在实际应用中有时会发现，滤波得到的状态估计与实际状态之间的误差，远远超出了按公式计算的范围；也就是说，按公式计算的误差方差阵可以逐渐地趋于零，而实际的估计误差却可以趋于无穷大。这种现象称为滤波的发散现象。显然，当滤波发散时，滤波就失去了它的最佳作用，因此，在实际应用中必须克服这种现象。

引起滤波发散的主要原因是：

(1) 由于对具体问题的物理过程了解得不准确，从而使物理系统的数学模型与实际过程不能较好的符合；或者虽然对物理过程了解得比较好，但其数学模型极为复杂（例如是非线性模型或状态维数过高），在简化数学模型时（如线性化或降维处理）时带来模型误差。

(2) 对于动态噪声或观测噪声的统计性质缺乏了解，因而取得不合适。

(3) 由于递推计算是在有穷字长的数字计算机上实现的，每步都有舍入误差，可能使所计算的误差方差阵逐渐失去正定性或对称性，造成计算值与理论值之间的偏差越来越大。

现用一个例子来说明发散现象。

在利用卡尔曼滤波方法估计气球以恒定速度 \(a\) 上升的高度，此时的实际状态方程和观测方程应为 \[\left.\begin{aligned} X_{k+1}&=X_{k}+a\\ L_{k+1}&=X_{k+1}+\varDelta_{k+1} \end{aligned}\right\} \tag{4-11-1}\] 式中 \(X_{k}\) 表示气球的高度。现在忽略了气球上升的速度，而采用了气球为静止时的状态方程，将状态方程和观测方程取为 \[\left.\begin{aligned} \overline{X}_{k+1}&=\overline{X}_{k}\\ L_{k+1}&=\overline{X}_{k+1}+\overline{\varDelta}_{k+1} \end{aligned}\right\} \tag{4-11-2}\] 且假设了 \(\overline{\varDelta}_{k}\) 为零均值白噪声序列，即设 \[E(\overline{\varDelta}_{k})=0,\quad \operatorname{cov}(\overline{\varDelta}_{k},\overline{\varDelta}_{j})=\overline{\sigma_{\varDelta}^{2}}\delta_{kj}\] 并设 \[\begin{gathered} E(\overline{X}_{0})=\mu_{X}(0),\quad \operatorname{var}(X_{0})=\overline{D}_{X}(0)\\ \operatorname{cov}(\overline{X}_{0},\overline{\varDelta}_{k})=0,\quad \overline{D}_{\varOmega}(k)=0 \end{gathered}\] 根据滤波方程可以解得 \[\left.\begin{aligned} \hat{\overline{X}}(k/k)&=\frac{\overline{\sigma_{\varDelta}^{2}}}{\sigma_{\varDelta}^{2}+k\overline{D}_{X}(0)}\mu_{X}(0) +\frac{\overline{D}_{X}(0)}{\sigma_{\varDelta}^{2}+k\overline{D}_{X}(0)}\sum_{i=1}^{k}L_{i}\\ \overline{D}_{\hat{X}}(k/k)&=\frac{\overline{D}_{X}(0)\,\overline{\sigma_{\varDelta}^{2}}}{\sigma_{\varDelta}^{2}+k\overline{D}_{X}(0)} =\frac{\overline{\sigma_{\varDelta}^{2}}}{\dfrac{\overline{\sigma_{\varDelta}^{2}}}{\overline{D}_{X}(0)}+k}\\ \overline{J}_{k}&=\frac{\overline{D}_{X}(0)}{\overline{\sigma_{\varDelta}^{2}}+K\overline{D}_{X}(0)} =\frac{\overline{\sigma_{\varDelta}^{2}}}{\dfrac{\overline{\sigma_{\varDelta}^{2}}}{\overline{D}_{X}(0)}+k} \end{aligned}\right\} \tag{4-11-3}\] 因为由 (4-11-1) 式可知 \[L_{k}=X_{K-1}+a+\varDelta_{k}=X_{0}+ka+\varDelta_{k}\] 所以有 \[\hat{\overline{X}}(k/k)=\frac{\overline{\sigma_{\varDelta}^{2}}}{\sigma_{\varDelta}^{2}+k\overline{D}_{X}(0)}\mu_{X}(0) +\frac{\overline{D}_{X}(0)}{\sigma_{\varDelta}^{2}+k\overline{D}_{X}(0)} \left[kX_{0}+\frac{1}{2}k(k+1)\alpha+\sum_{i=1}^{k}\varDelta_{i}\right] \tag{4-11-4}\] 又由于气球的真实高度是 \[X_{k}=X_{k-1}+\alpha=X_{0}+k\alpha\] 故 \(\hat{\overline{X}}(k/k)\) 的实际估计误差为 \[\begin{aligned} \widetilde{\overline{X}}(k/k)&=X_{k}-\hat{\overline{X}}(k/k)\\ &=\frac{1}{\dfrac{\overline{\sigma_{\varDelta}^{2}}}{\overline{D}_{X}(0)}+k} \left\{\frac{\overline{\sigma_{\varDelta}^{2}}}{\overline{D}_{X}(0)}\left[X_{0}-\mu_{X}(0)+k\alpha\right] +\frac{1}{2}\alpha k(k-1)-\sum_{i=1}^{k}\varDelta_{i}\right\} \end{aligned} \tag{4-11-5}\] 若设 \(\overline{D}_{X}(0)=\infty\)，即设初始状态的先验统计性质未知，则实际的估计误差为 \[\widetilde{\overline{X}}(k/k)=\frac{1}{2}\alpha(k-1)-\frac{1}{k}\sum_{i=1}^{k}\varDelta_{i} \tag{4-11-6}\]

而它的均方误差为 \[\hat{\overline{D}}_{\hat{X}}(k/k)=E\left[\widetilde{\overline{X}}(k/k)\,\widetilde{\overline{X}}^{T}(k/k)\right] =\frac{\alpha^{2}}{4}(k-1)^{2}+\frac{1}{k}\sigma_{\varDelta}^{2} \tag{4-11-7}\] 式中 \(\sigma_{\varDelta}^{2}\) 是实际观测误差 \(\varDelta_{k}\) 的方差。

由 (4-11-3) 和 (4-11-7) 式可以看到，当 \(\overline{D}_{X}(0)=\infty\) 时，按滤波公式算得的 \(\overline{D}_{\hat{X}}(k/k)=\dfrac{\overline{\sigma_{\varDelta}^{2}}}{k}\) 随着 \(k\) 的增加而减小；但实际的均方误差 \(\overline{D}_{\hat{X}}(k/k)\) 却随着 \(k\) 的增大而增大，直至无穷，即滤波是发散的。而且，即使没有完全忽略 \(\alpha\)，而只是取得不精确，也会出现发散现象。

下面根据引起滤波发散的原因来讨论克服发散现象的几种方法。

(1) 限定下界法

由上面的例子可知，引起滤波发散的一个直接原因，是 \(\overline{D}_{\hat{X}}(k/k)\) 随着 \(k\) 的增加而很快的减小，并趋于零，且 (4-11-3) 式中的增益阵 \(\overline{J}_{k}\) 也随着 \(k\) 的增加而迅速变小，使得新得到的观测数据的作用越来越小，从而产生数据饱和现象，使滤波发散。

根据引起发散的这个原因，可得到一种克服这类发散现象的方法，就是强制地使增益阵 \(\overline{J}_{k}\) 在过了一定时间之后不再下降，称为限定下界法。

例如，在上面的例子中，取一个适当的正数 \(M\)，使得当 \(k>M\) 时，新的增益阵 \(\overline{J_{k}^{*}}\) 不再下降，即令 \[\overline{J_{k}^{*}}= \begin{cases} 1\left/\left[\dfrac{\overline{\sigma_{\varDelta}^{2}}}{\overline{D}_{X}(0)}+k\right]\right., & k\leq M\\[14pt] \dfrac{1}{M}, & k>M \end{cases} \tag{4-11-8}\] 若 \(\overline{D}_{X}(0)=\infty\)，则令 \[\overline{J_{k}^{*}}= \begin{cases} \dfrac{1}{k}, & k\leq M\\[10pt] \dfrac{1}{M}, & k>M \end{cases} \tag{4-11-9}\] 此时，可得新的滤波值为 \[\hat{\overline{X}}^{*}(k/k)=(E-\overline{J_{k}^{*}})\,\hat{\overline{X}}^{*}(k-1/k-1)+\overline{J_{k}^{*}}L_{k} \tag{4-11-10}\] 如果初始状态 \(\hat{X}(0)=0\)，则得 \[\hat{\overline{X}}^{*}(k/k)=\sum_{i=1}^{k-1}(E-\overline{J_{k}^{*}})\,(E-\overline{J_{k-1}^{*}})\cdots(E-\overline{J_{i+1}^{*}})\,\overline{J_{i}^{*}}L_{i} +\overline{J_{k}^{*}}L_{k} \tag{4-11-11}\] 可以导出它的均方误差为 \[\hat{\overline{D}}_{\hat{X}}^{*}(k/k)=(M-1)^{2}\left[1-\frac{1}{2}\left(1-\frac{1}{M}\right)^{k-M}\right]^{2}\alpha^{2} +\frac{\sigma_{\varDelta}^{2}}{2M-1}\left[1+\left(1-\frac{1}{M}\right)^{2(k-M)+1}\right] \tag{4-11-12}\] 当 \(k\rightarrow\infty\) 时，有 \[\lim_{k\rightarrow\infty}\hat{\overline{D}}_{\hat{X}}^{*}(k/k)=(M-1)^{2}\alpha^{2}+\frac{\sigma_{\varDelta}^{2}}{2M-1} \tag{4-11-13}\] 这就克服了发散现象。显然，这种滤波方法并不是最佳的。为了起到真正的滤波作用，必须使滤波值的均方误差小于观测噪声的方差，即应选取 \(M>0\)，使得它满足 \[(M-1)^{2}\alpha^{2}+\frac{\sigma_{\varDelta}^{2}}{2M-1}<\sigma_{\varDelta}^{2}\] 因此，\(M\) 只能在下列范围内选取 \[1<M<1+\frac{\sigma_{\varDelta}}{\alpha} \tag{4-11-14}\] 可见，只有当被忽略的 \(\alpha\) 显著地小于 \(\sigma_{\varDelta}\) 时，这个方法才有实用意义。

(2) 扩充状态法（状态增广法）

将状态方程中被忽略的未知部分（如以上例子中的 \(\alpha\)）看做是一种状态，将它扩充到原有的状态向量中，如在例中取 \[\alpha_{k+1}=\alpha_{k}+\xi_{k} \tag{4-11-15}\] 在式中设 \(\xi_{k}\) 为零均值白噪声序列，即设 \[\left.\begin{aligned} E(\xi_{k})&=0\\ \operatorname{cov}(\xi_{k},\xi_{j})&=\sigma_{\xi}^{2}\delta_{kj} \end{aligned}\right\} \tag{4-11-16}\] 再将原来的状态表示为 \(X_{k}^{\prime}\)，令扩充后的状态向量为 \(X_{k}=\begin{bmatrix}X_{k}^{\prime}&\alpha_{k}\end{bmatrix}^{T}\)，则扩充后的状态方程和观测方程为 \[\left.\begin{aligned} X_{k+1}^{*}&=\begin{bmatrix}1&1\\ 0&1\end{bmatrix}X_{k}+\begin{bmatrix}0\\ 1\end{bmatrix}\xi_{k}\\ L_{k+1}&=\begin{bmatrix}1&0\end{bmatrix}X_{k+1}+\Delta k \end{aligned}\right\} \tag{4-11-17}\] 这样，问题就转化为如何适当地选择 \(\sigma_{\xi}\)，才使滤波不发散。因为此时有（\(m=2\)）： \[\begin{gathered} \sum_{i=0}^{m-1}(\bm{\varPhi}^{2})\,\bm{\varGamma}\bm{\varGamma}^{T}(\bm{\varPhi}^{i})^{T} =\begin{bmatrix}1&1\\ 1&2\end{bmatrix}>0\\ \sum_{i=0}^{m-1}(\bm{\varPhi}^{2})^{T}\bm{B}^{T}\bm{B}(\bm{\varPhi}^{i}) =\begin{bmatrix}2&1\\ 1&1\end{bmatrix}>0 \end{gathered}\] 所以，以上系统是一致完全能控和一致完全能观的，故其滤波的误差方差阵有一致的上界，即不会发散。

(3) 渐消记忆滤波

为了消除模型误差所引起的滤波发散现象，还可以设法加大新观测数据的作用，而相对地减小过老数据对滤波的影响。为了达到这个目的，可以采用将以前的 \(D_{X}(0)\)、\(D_{\varOmega}(k)\)、\(D_{\Delta}(k)\) 适当地逐渐加大的办法。

考虑线性系统模型 \[\left.\begin{aligned} X_{k+1}&=\bm{\varPhi}_{k+1,k}X_{k}+\bm{\varGamma}_{k+1,k}\varOmega_{k}\\ L_{k+1}&=B_{k+1}X_{k+1}+\varDelta_{k+1} \end{aligned}\right\} \tag{4-11-18}\] 其随机模型也与 4-3 节相同。为了加大新观测数据的作用，我们可以将初始状态的方差阵 \(D_{X}(0)\) 换为 \(\exp\left(\sum\limits_{i=1}^{N-1}c_{i}\right)D_{X}(0)\)，将动态噪声的方差阵 \(D_{\varOmega}(0)\)、\(D_{\varOmega}(1)\)，\(\cdots\)，\(D_{\varOmega}(N-1)\) 换为 \(\exp\left(\sum\limits_{i=1}^{N-1}c_{i}\right)\cdot D_{\varOmega}(0)\)，\(\exp\left(\sum\limits_{i=2}^{N-1}c_{i}\right)\cdot D_{\varOmega}(1)\)，\(\cdots\)，\(e^{c_{N-1}}D_{\varOmega}(N-2)\)，\(D_{\varOmega}(N-1)\)，而将观测噪声的方差阵 \(D_{\Delta}(1)\)，\(D_{\Delta}(2)\)，\(\cdots\)，\(D_{\Delta}(N)\) 换为 \(\exp\left(\sum\limits_{i=1}^{N-1}c_{i}\right)\cdot D_{\Delta}(1)\)，\(\exp\left(\sum\limits_{i=2}^{N-1}c_{i}\right)D_{\Delta}(2)\)，\(\cdots\)，\(e^{c_{N-1}}D_{\Delta}(N-1)\)，\(D_{\Delta}(N)\)。这种代换也相当于给出一个新的模型 \[\left.\begin{aligned} X_{k+1}^{*}&=\bm{\varPhi}_{k+1,k}X_{k}^{*}+\bm{\varGamma}_{k+1,k}\varOmega_{k}^{*}\\ L_{k+1}&=B_{k+1}X_{k+1}^{*}+\varDelta_{k+1}^{*} \end{aligned}\right\} \tag{4-11-19}\] 不难证明，模型 (4-11-19) 在 \(N\) 时刻的滤波值就是模型 (4-11-18) 在 \(N\) 时刻的渐消记忆滤波。可以看到，对于 \(k\leq N\) 的每一步计算，由于一步预测值 \(\hat{X}^{*}(k/k-1)\) 的误差方差阵 \(D_{X}^{*}(k/k-1)\) 表示式的第一项乘上了一个大于 1 的因子 \(e^{c_{N-1}}\)，因此 \(\hat{X}^{*}(k/k)\) 的误差方差阵 \(D_{X}^{*}(k/k)\) 也增大了，从而使每步的增益阵 \(J_{k}^{*}\) 也增大了，这就加大了每次新观测数据 \(L_{k}\) 的作用。至于常数 \(c_{i}\) 到底如何选择，则要根据具体情况决定。

(4) 限定记忆滤波

为了消除“太老”数据对现时状态估计的影响，在求 \(X_{k}\) 的滤波值 \(\hat{X}(k/k)\) 时，只利用 \(k\) 时刻的前 \(N\) 个时刻的观测数据 \(L_{k-N+1}\)，\(L_{k-N+2}\)，\(\cdots\)，\(L_{k}\)，而把其余的数据全部去掉，这里的 \(N\) 是根据具体对象选定的记忆长度，通常称这种方法为限定记忆滤波。

(5) 自适应滤波

所谓自适应滤波，就是在利用观测数据进行滤波时，不断地对未知的或不确切知道的模型参数以及噪声的统计性质进行估计和修正，以减少模型误差。因此，它也是克服滤波发散的途径之一。自适应滤波是一个需要专门论述的课题，这里就不讨论了。