2　系统状态滤波预报方法

2.1　系统滤波递推方法

1.观测滤波方程

根据t=tk时已观测的状态变量Z（tk）=[z（t0），z（t1）…z（tk）]T，按无偏最小方差估计原理，由观测方程得t=tk时的状态变量最优滤波估计值为

和C（tk|tk）的估计误差协方差阵：

其中K（tk）称卡尔曼（Kalman）滤波增益矩阵，依下式计算：

P（tk|tk-1）为tk-1时预报的状态变量的误差协方差阵，计算式为

式中：Q、R分别为系统噪声W（t）和观测噪声M（t）的方差，白噪声序列时，Q、R为常数，通过优选确定；Φ（tk）为状态由C（tk-1）向C（tk）转变的转移矩阵，与A的关系为dΦ/dt=AΦ。

通过对观测值的滤波计算，由式（7）、式（8）求得tk时状态最优滤波估计C（tk|tk）和P（tk|tk），以便进一步由系统方程进行状态预报。

2.状态滤波预报方程

按系统预报误差为无偏估计的要求，由式（5）得状态预报微分方程为：dC^（t）/dt=AC（t）+BU（t），该式为齐次线性微分方程，由积分变换法解得

同时还可求得

当系统参量A、B、G、Q、R在时域（tk-1，tk）间保持不变和观测时距固定的情况下，式（11）、式（12）的形式是方便的，这种情况下，Φ（t）和式（12）右边第二项对所有的观测时域都是相同的，只需计算一次就可以了，如果系统为时变的，那么状态预报方程直接采用微分方程的形式更为有效，即：

实用中，采用吉尔方法求解式（13）、式（14），由tk-1积分到t便得到[tk-1，tk]间的、P中间过程值和P（t|tk-1），其终值即和P（tk|tk-1）。

2.2　状态滤波预报计算过程

按上述原理，滤波递推计算进行水质过程预报的步骤简要归纳如下：

（1）取预报开始时（t=tk-1=0）的滤波估计值和P（tk-1|tk-1）=P（0|0），一般按经验优化选取。

（2）按式（13）预报时段内t时刻和时段未tk时刻的各单元水质浓度、，同时按式（14）计算相应的预报误差协方差阵P（t|tk-1）、P（tk|tk-1）。

（3）按式（9）计算增益矩阵K（tk）。

（4）由式（7）、式（8）计算t=tk时的状态最优滤波估计和相应的协方差阵P（tk|tk）。

（5）将tk赋予tk-1，回到第2步，由和P（tk-1|tk-1）重新开始递推计算及预报。如此不断地进行预报-滤波-再预报-再滤波，从而预报出各单元水质变化过程。

2.3　算法稳定性分析

状态方程和观测方程中各系数矩阵A、B、G、H在一定时域内为常数阵，且取系统噪声和观测噪声为平稳白噪声，Q＞0、R＞0，由柯莱姆矩阵性质可得下面的可控性和可观性条件：

式中：n为状态变量的维数，在这里即整个单元数，对于滇池n=31，上述建立的系统滤波模型能满足式（15）的两个条件，从而保证了该滤波算法的稳定性。