CN103414450B - 噪声统计特性未知系统的实时多速率h∞融合滤波方法 - Google Patents

Abstract

本发明公开了一种噪声统计特性未知系统的实时多速率H∞融合滤波方法，包括以下步骤：A:在由多个传感器构成的多速率多传感器时变离散系统中，当每个测量值到达时，融合中心首先在计数器Ni记录该滤波周期内已到达融合中心的量测个数；B:利用计数器给出有限域H∞性能指标函数，作为滤波过程的约束指标；C:在上述性能指标的约束下，进行测量值的融合滤波；D:当该滤波周期内所有测量值都到达融合中心时，得到系统待估信号基于全局信息的融合估计结果。本发明通过最优分析的方法分析各传感器以不同速率工作时的系统状态与系统噪声之间的耦合关系，实现多速率多传感器系统的实时融合滤波问题，能够同时适用于时不变系统及时变系统的多速率H∞融合滤波过程。

Description

噪声统计特性未知系统的实时多速率H∞融合滤波方法

技术领域

本发明涉及一种滤波方法，尤其涉及一种针对多个传感器构成的多速率多传感器时变离散系统在噪声统计特性未知情况下的实时多速率H∞融合滤波方法。

背景技术

现代工业工程日趋复杂化，多传感器信息融合技术及在此类系统中的应用，日益受到广泛关注。在现有的多传感器信息融合研究中，研究较多的是同步融合问题,即假设所有传感器具有相同的采样率，且均为同时采样，没有固有延迟与通信延迟，同时把数据送到融合中心。但是，实际系统中所使用的传感器类型异构多样、布设位置差异较大、采样起始时刻不尽相同，因此，亟需解决多传感器异步融合问题。多速率多传感器系统是一类常见的异步融合系统。在多速率多传感器系统中，各传感器的采样频率不尽相同，各传感器的起始采样时刻也可能存在差异。这些特点都给多速率多传感器系统融合滤波过程带来了困难。闫莉萍、胡艳艳等学者提出了基于Kalman滤波器的多速率多传感器系统的融合滤波方法，该方法要求系统噪声必须满足方差均值已知的高斯分布，但是在实际系统中系统噪声不都是服从高斯分布的，且噪声统计特性在短时间内通常都是难以准确获知的；为此，梁彦、陈通文等学者在假设系统噪声在整个时间域内能量有界的前提下，提出了一种基于各传感器采样周期最小公倍数为提升融合周期的多速率多传感器系统H∞融合滤波方法，但是该方法仅适用于线性时不变系统，当系统参数随时间发生变化时，该类方法难以做出有效响应，且需要以提升融合周期为融合处理周期，使得融合滤波过程在实时性方面具有一定不足。此外，该类方法要求系统噪声在[0,+∞)整个时间域内能量有界这一条件意味着在一定时间之后，系统噪声趋于零，这与实际系统略有不符。

发明内容

本发明的目的是基于实际系统中噪声不会呈现为无穷大的这一事实，利用系统噪声在有限时间域内能量有界的性质，提供噪声统计特性未知系统的实时多速率H∞融合滤波方法，能够克服现有技术不足，对时不变系统和时变系统进行多速率H∞融合滤波。

本发明采用下述技术方案：

一种噪声统计特性未知系统的实时多速率H∞融合滤波方法，包括以下步骤：

A:在由多个传感器构成的多速率多传感器时变离散系统中，当每个测量值到达时，融合中心首先在计数器N_i记录该滤波周期内已到达融合中心的量测个数；

B:利用计数器给出有限域H∞性能指标函数，作为滤波过程的约束指标；

C:在上述性能指标的约束下，进行测量值的融合滤波；

D:当该滤波周期内所有测量值都到达融合中心时，即可得到系统待估信号基于全局信息的融合估计结果。

所述的A步骤中，多速率多传感器时变离散系统如下所示

x(k)=F(k,k-1)x(k-1)+w(k,k-1);

y_i(kn_i+j)=H_i(kn_i+j)x(kn_i+j)+v_i(kn_i+j),i=1,2,…,p;

z(km)=L(km)x(km);

其中，x(k)为演变周期为h的离散系统状态；y_i(kn_i+j)表示第i个传感器第k次采样得到的测量值，其采样起始时刻为j时刻，采样间隔为n_i；kn_i+j是k×n_i+j的缩写形式，表示其采样时刻；z(km)表示滤波周期为m×h（缩写为m）的系统待估信号在k×m×h（缩写为km）时刻的表达式；n_i,m均为正整数；F(k,k-1),H_i(kn_i+j),L(km)分别为系统状态从k-1时刻到k时刻的转移矩阵、传感器i第k次采样时的观测矩阵以及系统在第k个滤波周期的待估信号加权矩阵；w(k,k-1)为系统从k-1时刻到k时刻积累的过程噪声，v_i(kn_i)为第i个传感器第k次采样时的测量噪声，两类噪声的统计特性未知；该多速率多传感器时变离散系统含有p个传感器。

所述的C步骤中，进行测量值的融合滤波包括以下步骤

c1:利用测量值建立当前滤波时刻的伪量测矩阵；

c2：估计系统噪声及其与系统状态之间的耦合关系；

c3：对测量值进行滤波实现待估信号估计。

所述的B步骤中，有限域H∞性能指标函数为

$\underset{w,v\&Element;l_2}{\sup }\frac{\underset{\underset{N_i\&NotEqual;0}{i=1}}{\overset{k-1}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N_i}{\mathrm{\&Sigma;}}}{e}_{z,j}^T\left(\mathrm{im}\right)e_{z,j}\left(\mathrm{im}\right)+\underset{j=1}{\overset{l}{\mathrm{\&Sigma;}}}{e}_{z,j}^T\left(\mathrm{km}\right)e_{z,j}\left(\mathrm{km}\right)}{\underset{\underset{N_i\&NotEqual;0}{i=1}}{\overset{k-1}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N_i}{\mathrm{\&Sigma;}}}{v}_{{\mathrm{\&alpha;}}_j^i}^T\left({\mathrm{\&beta;}}_j^i\right)v_{{\mathrm{\&alpha;}}_j^i}\left({\mathrm{\&beta;}}_j^i\right)+\underset{j=1}{\overset{l}{\mathrm{\&Sigma;}}}{v}_{{\mathrm{\&alpha;}}_j^k}^T\left({\mathrm{\&beta;}}_j^k\right)v_{{\mathrm{\&alpha;}}_j^k}\left({\mathrm{\&beta;}}_j^k\right)+\underset{i=1}{\overset{\mathrm{mk}}{\mathrm{\&Sigma;}}}{w}^T(i,i-1)w(i,i-1)+{(x\left(0\right)-{\hat{x}}_0)}^T{P_0}^{-1}(x\left(0\right)-{\hat{x}}_0)}<{\mathrm{\&gamma;}}^2$

其中，e_z,j(ik)表示融合中心在第i个滤波周期利用第j个到达的测量值对待估信号进行滤波时得到的估计误差，γ为给定的有限域H∞性能指标值； 表示在第i个滤波周期第j个到达的测量值对应的测量噪声，表示采集该测量信息的传感器标号，表示该测量信息的采集时刻；e_z,j(im)表示融合中心在第i个滤波周期利用第j个到达的测量值对待估信号z(im)进行滤波时得到的估计误差，；x(0),分别表示初始状态及其估计值，P₀表示偏离x(0)的程度。

所述的c1步骤中，

对应的伪量测噪声为： $v_{y,l}^{*}\left(\mathrm{km}\right)=-H_l^{*}\left(\mathrm{km}\right)w(\mathrm{km},{\mathrm{\&beta;}}_l^k);$

对应的伪量测矩阵为： $H_l^{*}\left(\mathrm{km}\right)=H_{{\mathrm{\&alpha;}}_l^k}\left({\mathrm{\&beta;}}_l^k\right)F({\mathrm{\&beta;}}_l^k,\mathrm{km});$

对应的伪测量信息为： $y_l^{*}\left(\mathrm{km}\right)=y_{{\mathrm{\&alpha;}}_l^k}\left({\mathrm{\&beta;}}_l^k\right)=H_l^{*}\left(\mathrm{km}\right)x\left(\mathrm{km}\right)+v_{y,l}^{*}\left(\mathrm{km}\right);$

其中，对第k个滤波周期内利用第N_k=l个到达融合中心的测量值 进行滤波，表示在第i个滤波周期内到达融合中心的第j个量测，第j个量测是由第个传感器在时刻采样得到的；w(km,)表示时刻到km时刻积累的过程噪声；F(km,)表示时刻到km时刻系统状态转移矩阵；F(,km)为其(伪)逆矩阵。

所述的c2步骤中，在估计系统噪声及其与系统状态之间的耦合关系时，首先求出下列辅助参量：

（1）求解状态与测量预测值的互Gramian矩阵

$R_{\mathrm{xyz},l}\left(\mathrm{km}\right)=\left[\begin{array}{cc}{R}_{\mathrm{xy},l}\left(\mathrm{km}\right)& R_{\mathrm{xz},l}\left(\mathrm{km}\right)\end{array}\right]$

$=P_l\left(\mathrm{km}\right)[{\left(H_l^{*}\left(\mathrm{km}\right)\right)}^T,L^T\left(\mathrm{km}\right)]-[(\underset{t={\mathrm{\&beta;}}_l^k}{\overset{\mathrm{km}}{\mathrm{\&Sigma;}}}F(t,t-1)F^T(t,t-1)-P_{\mathrm{xw},l}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)){\left(H_l^{*}\left(\mathrm{km}\right)\right)}^T,0]$

其中，P_xw,l(km,|(k-1)m))为状态与过程噪声的互Gramian矩阵；P_l(km)为在第k个滤波周期利用第l个到达的测量值进行滤波时的Riccati变量，其满足如下递推关系：

$\left\{\begin{array}{c}{P}_{l+1}\left(\mathrm{km}\right)=P_l\left(\mathrm{km}\right)-K_l\left(\mathrm{km}\right)R_{\mathrm{xyz},l}^T\left(\mathrm{km}\right)\\ P_1\left((k+1)m\right)=F((k+1)m,\mathrm{km})(P_{N_k}\left(\mathrm{km}\right)-K_{N_k}\left(k\right)R_{\mathrm{xy},N_k}^T\left(k\right))F^T((k+1)m,\mathrm{km})+Q((k+1)m,\mathrm{km})\end{array}\right.$

其中，Q((k+1)m,km)为km时刻到(k+1)m时刻积累的过程噪声的Gramian矩阵； $Q((k+1)m,\mathrm{km})=\underset{i=\mathrm{km}}{\overset{\mathrm{km}+m-1}{\mathrm{\&Sigma;}}}F(i+1,i)F^T(i+1,i);$ (km)表示在第k个滤波周期利用第N_k个到达的测量值进行滤波时的Riccati变量；(k)为在第k个滤波周期利用第N_k个到达的测量值进行滤波时的滤波增益矩阵；(k)为在第k个滤波周期利用第N_k个到达的测量值进行滤波时的系统状态与测量预测之间的互Gramian矩阵；

（2）求解测量预测值的Gramian矩阵

$R_{\mathrm{eyz},l}\left(\mathrm{km}\right)=\left[\begin{array}{cc}{R}_{\mathrm{ey},l}\left(\mathrm{km}\right)& R_{\mathrm{yz},l}\left(\mathrm{km}\right)\\ R_{\mathrm{zy},l}\left(\mathrm{km}\right)& R_{\mathrm{ez},l}\left(\mathrm{km}\right)\end{array}\right]$

$=\left[\begin{array}{c}{H}_l^{*}\left(\mathrm{km}\right)\\ L\left(\mathrm{km}\right)\end{array}\right]P_l\left(\mathrm{km}\right)[{\left(H_l^{*}\left(\mathrm{km}\right)\right)}^T,L^T\left(\mathrm{km}\right)]$

$+\left[\begin{array}{cc}{H}_l^{*}\left(\mathrm{km}\right)(Q(\mathrm{km},{\mathrm{\&beta;}}_l^k)-R_{\mathrm{ww},l}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)){\left(H_l^{*}\left(\mathrm{km}\right)\right)}^T& 0\\ 0& -{\mathrm{\&gamma;}}^{2}I\end{array}\right]$

$-\left[\begin{array}{cc}{H}_l^{*}\left(\mathrm{km}\right)R_{\mathrm{xw},l}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m){\left(H_l^{*}\left(\mathrm{km}\right)\right)}^T& 0\\ 0& 0\end{array}\right]$

$-{\left[\begin{array}{cc}{H}_l^{*}\left(\mathrm{km}\right)R_{\mathrm{xw},l}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m){\left(H_l^{*}\left(\mathrm{km}\right)\right)}^T& 0\\ 0& 0\end{array}\right]}^T$

其中，R_eyz,l(km)表示在第k个滤波周期利用第l个到达的测量值进行滤波时的测量预测值的Gramian矩阵；R_xw,l(km,|(k-1)m)为在第k个滤波周期利用第l个到达的测量值进行滤波时状态与过程噪声估计值的互Gramian矩阵；R_ww,l(km,|(k-1)m)为在第k个滤波周期利用第l个到达的测量值进行滤波时过程噪声估计误差的Gramian矩阵；Q(km,)为在第k个滤波周期利用第l个到达的测量值进行滤波时对应耦合过程噪声的Gramian矩阵；

$Q(\mathrm{km},{\mathrm{\&beta;}}_l^k)=\underset{i={\mathrm{\&beta;}}_l^k}{\overset{\mathrm{km}-1}{\mathrm{\&Sigma;}}}F(i+1,i)F^T(i+1,i);$

（3）状态与过程噪声估计值的互Gramian矩阵

$\left\{\begin{array}{c}{P}_{\mathrm{xw},l}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)=P_{\mathrm{xw},l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)\\ +R_{\mathrm{xyz},l-1}\left(\mathrm{km}\right)R_{\mathrm{eyz},l-1}^{-1}\left(\mathrm{km}\right)R_{\mathrm{wyz},l-1}^T(\mathrm{km},{\mathrm{\&beta;}}_l^k,\mathrm{km}),j>1\\ P_{\mathrm{xw},1}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)=0\end{array}\right.$ ；

其中，P_xw,l(km,|(k-1)m)在第k个滤波周期利用第l个到达的测量值进行滤波时，状态与过程噪声估计值的互Gramian矩阵；

（4）过程噪声估计值与测量值的互Gramian矩阵

$R_{\mathrm{wyz},l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k,\mathrm{km})=[R_{\mathrm{wy},l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k,\mathrm{km}),R_{\mathrm{wz},l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k,\mathrm{km})]$

$=R_{\mathrm{xw},l-1}^T(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m){\left[\begin{array}{c}{H}_{l-1}^{*}\left(\mathrm{km}\right)\\ L\left(\mathrm{km}\right)\end{array}\right]}^T-[R_{\mathrm{ww},l-1}^T(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m){\left(H_{l-1}^{*}\left(\mathrm{km}\right)\right)}^T,0]$ ；

其中，R_wyz,l-1(km,,km)在第k个滤波周期利用第l个到达的测量值进行滤波时，过程噪声估计值与测量值的互Gramian矩阵；R_ww,l(km,|(k-1)m)为在第k个滤波周期利用第l个到达的测量值进行滤波时过程噪声估计误差的Gramian矩阵；

（5）过程噪声估计误差的Gramian矩阵

$R_{\mathrm{ww},l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k,|(k-1)m)={Q(\mathrm{km},{\mathrm{\&beta;}}_l^k)-P}_{\mathrm{ww},l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k,|(k-1)m)$ ；

其中，P_ww,l(km,|(k-1)m)为在第k个滤波周期利用第l个到达的测量值进行滤波时过程噪声估计值的Gramian矩阵；Q(km,)为在第k个滤波周期利用第l个到达的测量值进行滤波时过程噪声的Gramian矩阵；

$Q(\mathrm{km},{\mathrm{\&beta;}}_l^k)=\underset{i={\mathrm{\&beta;}}_l^k}{\overset{\mathrm{km}-1}{\mathrm{\&Sigma;}}}F(i+1,i)F^T(i+1,i);$

（6）过程噪声估计值的Gramian矩阵

$\left\{\begin{array}{c}{P}_{\mathrm{ww},l}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m))=P_{\mathrm{ww},l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)\\ +R_{\mathrm{wyz},l-1}^T(\mathrm{km},{\mathrm{\&beta;}}_l^k,\mathrm{km})R_{\mathrm{eyz},l-1}^{-1}\left(\mathrm{km}\right)R_{\mathrm{wyz},l-1}^T(\mathrm{km},{\mathrm{\&beta;}}_l^k,\mathrm{km}),j>1\\ P_{\mathrm{ww},1}(k,{\mathrm{\&beta;}}_l^k|k-1)m)=0\end{array}\right.$ ；

其中，R_ww,l(km,|(k-1)m)为在第k个滤波周期利用第l个到达的测量值进行滤波时过程噪声估计误差的Gramian矩阵；

然后计算出过程噪声估计值为

当l=1时， ${\hat{w}}_1(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)=0;$

当l>1时，

${\hat{w}}_l(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)$

${=\hat{w}}_{l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)+R_{\mathrm{wy},l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k,\mathrm{km}){\left(R_{\mathrm{ey},l-1}\left(\mathrm{km}\right)\right)}^{-1}{e}_{y,l-1}\left(\mathrm{km}\right)$ ；

其中，为在第k个滤波周期利用第l个到达的测量值进行滤波时对过程噪声的估计值；

最终得到伪量测噪声估计值为：

${\hat{v}}_{y,l}^{*}\left(\mathrm{km}\right|(k-1)m)=-H_l^{*}\left(\mathrm{km}\right){\hat{w}}_l(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m);$

当辅助参量R_ey,l(km)满足R_ey,l(km)>0，且

$R_{\mathrm{ez},j}\left(\mathrm{km}\right)-R_{\mathrm{zy},l}\left(\mathrm{km}\right)R_{\mathrm{ey},l}^{-1}\left(\mathrm{km}\right)R_{\mathrm{yz},l}\left(\mathrm{km}\right)<0$ 时，进入步骤c3，否则，需要重新设定性能指标值γ，然后重新执行步骤c2。

所述的c3步骤中，利用伪量测噪声估计值对测量值进行滤波包含以下步骤

（1）计算得出系统状态预测

当l=1时，

${\hat{x}}_1\left(\mathrm{km}\right|(k-1)m)=F(\mathrm{km},(k-1)m){\hat{x}}_{N_{k-1}}\left((k-1)m\right|(k-1)m);$

当l>1时，

${\hat{x}}_l\left(\mathrm{km}\right|(k-1)m)={\hat{x}}_{l-1}\left(\mathrm{km}\right|\mathrm{km});$

（2）计算得出观测预测

${\hat{y}}_l^{*}\left(\mathrm{km}\right|(k-1)m)=H_l^{*}\left(\mathrm{km}\right){\hat{x}}_l\left(\mathrm{km}\right|(k-1)m)+{\hat{v}}_{y,l}^{*}\left(\mathrm{km}\right|(k-1)m);$

（3）计算得出滤波增益

$K_{y,l}\left(\mathrm{km}\right)=R_{\mathrm{xy},l}\left(\mathrm{km}\right)R_{\mathrm{ey},l}^{-1}\left(\mathrm{km}\right);$

（4）计算得出系统状态估计

${\hat{x}}_l\left(\mathrm{km}\right|\mathrm{km})={\hat{x}}_l\left(\mathrm{km}\right|(k-1)m)+K_{y,l}\left(\mathrm{km}\right)(y_l^{*}\left(\mathrm{km}\right)-{\hat{y}}_l^{*}\left(\mathrm{km}\right|(k-1)m))$ ；

（5）计算得出系统待估信号估计值

$\hat{z}\left(\mathrm{km}\right|\mathrm{km})=L\left(\mathrm{km}\right){\hat{x}}_l\left(\mathrm{km}\right|\mathrm{km});$

当下一到达测量值仍为滤波器在第k个滤波周期内采样所得的测量值，即时，N_k=l+1，则返回步骤c1。

所述的D步骤中，当该滤波周期内所有测量值都到达融合中心，可得到系统待估信号基于全局信息的融合估计结果：

$\hat{z}\left(\mathrm{km}\right|\mathrm{km})=L\left(\mathrm{km}\right){\hat{x}}_{N_k}\left(\mathrm{km}\right|\mathrm{km}).$

本发明提供了噪声统计特性未知系统的实时多速率H∞融合滤波方法，通过最优分析的方法来深入分析各传感器以不同速率工作时的系统状态与系统噪声之间的耦合关系，并采用序贯式融合的思想，实现多速率多传感器系统的实时融合滤波问题。不仅约束了不同滤波周期间系统噪声对待估信号估计误差的影响，还在同一滤波周期内，约束了序贯式利用各量测对系统状态及待估信号的估计过程。能够同时适用于时不变系统及时变系统的多速率H∞融合滤波过程。

附图说明

图1为本发明的流程图；

图2为仿真算例场景示意图；

图3为仿真算例中待估信号的估计曲线图；

图4为仿真算例中待估信号的绝对估计误差曲线。

具体实施方式

如图1所示，本发明所述的噪声统计特性未知系统的实时多速率H∞融合滤波方法，包括以下步骤：

A:在由多个传感器构成的多速率多传感器离散系统中，当每个测量值到达时，融合中心首先在计数器N_i记录该滤波周期内已到达融合中心的测量值个数；

多速率多传感器时变离散系统如下所示

x(k)=F(k,k-1)x(k-1)+w(k,k-1);

y_i(kn_i+j)=H_i(kn_i+j)x(kn_i+j)+v_i(kn_i+j),i=1,2,…,p;

z(km)=L(km)x(km);

其中，x(k)为演变周期为h的离散系统状态；y_i(kn_i+j)表示第i个传感器第k次采样得到的测量值，其采样起始时刻为j时刻，采样间隔为n_i；kn_i+j是k×n_i+j的缩写形式，表示其采样时刻；z(km)表示滤波周期为m×h（缩写为m）的系统待估信号在k×m×h（缩写为km）时刻的表达式；n_i,m均为正整数；F(k,k-1),H_i(kn_i+j),L(km)分别为系统状态从k-1时刻到k时刻的转移矩阵、传感器i第k次采样时的观测矩阵以及系统在第k个滤波周期的待估信号加权矩阵；w(k,k-1)为系统从k-1时刻到k时刻积累的过程噪声，v_i(kn_i)为第i个传感器第k次采样时的测量噪声，两类噪声的统计特性未知；该多速率多传感器时变离散系统含有p个传感器。

B:利用计数器给出有限域H∞性能指标函数，作为滤波过程的约束指标；

有限域H∞性能指标函数为

$\underset{w,v\&Element;l_2}{\sup }\frac{\underset{\underset{N_i\&NotEqual;0}{i=1}}{\overset{k-1}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N_i}{\mathrm{\&Sigma;}}}{e}_{z,j}^T\left(\mathrm{im}\right)e_{z,j}\left(\mathrm{im}\right)+\underset{j=1}{\overset{l}{\mathrm{\&Sigma;}}}{e}_{z,j}^T\left(\mathrm{km}\right)e_{z,j}\left(\mathrm{km}\right)}{\underset{\underset{N_i\&NotEqual;0}{i=1}}{\overset{k-1}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N_i}{\mathrm{\&Sigma;}}}{v}_{{\mathrm{\&alpha;}}_j^i}^T\left({\mathrm{\&beta;}}_j^i\right)v_{{\mathrm{\&alpha;}}_j^i}\left({\mathrm{\&beta;}}_j^i\right)+\underset{j=1}{\overset{l}{\mathrm{\&Sigma;}}}{v}_{{\mathrm{\&alpha;}}_j^k}^T\left({\mathrm{\&beta;}}_j^k\right)v_{{\mathrm{\&alpha;}}_j^k}\left({\mathrm{\&beta;}}_j^k\right)+\underset{i=1}{\overset{\mathrm{mk}}{\mathrm{\&Sigma;}}}{w}^T(i,i-1)w(i,i-1)+{(x\left(0\right)-{\hat{x}}_0)}^T{P_0}^{-1}(x\left(0\right)-{\hat{x}}_0)}<{\mathrm{\&gamma;}}^2$

其中，e_z,j(im)表示融合中心在第i个滤波周期利用第j个到达的测量值对待估信号z(im)进行滤波时得到的估计误差，γ为给定的有限域H∞性能指标值。 表示在第i个滤波周期第j个到达的测量值对应的测量噪声，表示采集该测量信息的传感器标号，表示该测量信息的采集时刻。N_i表示在第i个滤波周期内有N_i个测量值已经到达融合中心,在这里以N_k=l为例进行滤波过程说明。x(0),分别表示初始状态及其估计值，P₀表示偏离x(0)的程度。

C:在上述性能指标的约束下，按照以下步骤进行测量值的融合滤波

c1:利用测量值 建立当前滤波时刻(km时刻)的伪量测矩阵；

对应的伪量测噪声为： $v_{y,l}^{*}\left(\mathrm{km}\right)=-H_l^{*}\left(\mathrm{km}\right)w(\mathrm{km},{\mathrm{\&beta;}}_l^k);$

对应的伪量测矩阵为： $H_l^{*}\left(\mathrm{km}\right)=H_{{\mathrm{\&alpha;}}_l^k}\left({\mathrm{\&beta;}}_l^k\right)F({\mathrm{\&beta;}}_l^k,\mathrm{km});$

对应的伪测量信息为： $y_l^{*}\left(\mathrm{km}\right)=y_{{\mathrm{\&alpha;}}_l^k}\left({\mathrm{\&beta;}}_l^k\right)=H_l^{*}\left(\mathrm{km}\right)x\left(\mathrm{km}\right)+v_{y,l}^{*}\left(\mathrm{km}\right);$

其中，w(km,)表示时刻到km时刻积累的过程噪声；F(km,)表示时刻到km时刻系统状态转移矩阵，F(,km)为其(伪)逆矩阵。

c2：估计系统噪声及其与系统状态之间的耦合关系；

（1）求解状态与测量预测值的互Gramian矩阵

$R_{\mathrm{xyz},l}\left(\mathrm{km}\right)=\left[\begin{array}{cc}{R}_{\mathrm{xy},l}\left(\mathrm{km}\right)& R_{\mathrm{xz},l}\left(\mathrm{km}\right)\end{array}\right]$

$=P_l\left(\mathrm{km}\right)[{\left(H_l^{*}\left(\mathrm{km}\right)\right)}^T,L^T\left(\mathrm{km}\right)]-[(\underset{t={\mathrm{\&beta;}}_l^k}{\overset{\mathrm{km}}{\mathrm{\&Sigma;}}}F(t,t-1)F^T(t,t-1)-P_{\mathrm{xw},l}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)){\left(H_l^{*}\left(\mathrm{km}\right)\right)}^T,0]$

其中，P_xw,l(km,|(k-1)m))为在第k个滤波周期利用第l个到达的测量值进行滤波时状态与过程噪声估计值的互Gramian矩阵，其计算方法将在（3）状态与过程噪声估计值的互Gramian矩阵中详细介绍；P_l(km)为在第k个滤波周期利用第l个到达的测量值进行滤波时的Riccati变量，其满足如下递推关系

$\left\{\begin{array}{c}{P}_{l+1}\left(\mathrm{km}\right)=P_l\left(\mathrm{km}\right)-K_l\left(\mathrm{km}\right)R_{\mathrm{xyz},l}^T\left(\mathrm{km}\right)\\ P_1\left((k+1)m\right)=F((k+1)m,\mathrm{km})(P_{N_k}\left(\mathrm{km}\right)-K_{N_k}\left(k\right)R_{\mathrm{xy},N_k}^T\left(k\right))F^T((k+1)m,\mathrm{km})+Q((k+1)m,\mathrm{km})\end{array}\right.$

在这里，Q((k+1)m,km)为km时刻到(k+1)m时刻累积的过程噪声的Gramian矩阵； $Q((k+1)m,\mathrm{km})=\underset{i=\mathrm{km}}{\overset{\mathrm{km}+m-1}{\mathrm{\&Sigma;}}}F(i+1,i)F^T(i+1,i);$ (km)表示在第k个滤波周期利用第N_k个到达的测量值进行滤波时的Riccati变量；(k)为在第k个滤波周期利用第N_k个到达的测量值进行滤波时的滤波增益矩阵，其具体求解方法将在下文c3(3)计算得出滤波增益中详细介绍；(k)为在第k个滤波周期利用第N_k个到达的测量值进行滤波时的系统状态与测量预测之间的互Gramian矩阵。

（2）求解测量预测值的Gramian矩阵

$R_{\mathrm{eyz},l}\left(\mathrm{km}\right)=\left[\begin{array}{cc}{R}_{\mathrm{ey},l}\left(\mathrm{km}\right)& R_{\mathrm{yz},l}\left(\mathrm{km}\right)\\ R_{\mathrm{zy},l}\left(\mathrm{km}\right)& R_{\mathrm{ez},l}\left(\mathrm{km}\right)\end{array}\right]$

$=\left[\begin{array}{c}{H}_l^{*}\left(\mathrm{km}\right)\\ L\left(\mathrm{km}\right)\end{array}\right]P_l\left(\mathrm{km}\right)[{\left(H_l^{*}\left(\mathrm{km}\right)\right)}^T,L^T\left(\mathrm{km}\right)]$

$+\left[\begin{array}{cc}{H}_l^{*}\left(\mathrm{km}\right)(Q(\mathrm{km},{\mathrm{\&beta;}}_l^k)-R_{\mathrm{ww},l}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)){\left(H_l^{*}\left(\mathrm{km}\right)\right)}^T& 0\\ 0& -{\mathrm{\&gamma;}}^{2}I\end{array}\right]$

$-\left[\begin{array}{cc}{H}_l^{*}\left(\mathrm{km}\right)R_{\mathrm{xw},l}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m){\left(H_l^{*}\left(\mathrm{km}\right)\right)}^T& 0\\ 0& 0\end{array}\right]$

$-{\left[\begin{array}{cc}{H}_l^{*}\left(\mathrm{km}\right)R_{\mathrm{xw},l}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m){\left(H_l^{*}\left(\mathrm{km}\right)\right)}^T& 0\\ 0& 0\end{array}\right]}^T$ ；

其中，R_eyz,l(km)表示在第k个滤波周期利用第l个到达的测量值进行滤波时的测量预测值的Gramian矩阵；R_xw,l(km,|(k-1)m)为在第k个滤波周期利用第l个到达的测量值进行滤波时状态与过程噪声估计值的互Gramian矩阵，其计算方法将在（3）状态与过程噪声估计值的互Gramian矩阵中详细介绍；R_ww,l(km,|(k-1)m)为在第k个滤波周期利用第l个到达的测量值进行滤波时过程噪声估计误差的Gramian矩阵，其计算方法将在（5）过程噪声估计误差的Gramian矩阵中详细介绍；Q(km,)为在第k个滤波周期利用第l个到达的测量值进行滤波时对应耦合过程噪声的Gramian矩阵；

$Q(\mathrm{km},{\mathrm{\&beta;}}_l^k)=\underset{i={\mathrm{\&beta;}}_l^k}{\overset{\mathrm{km}-1}{\mathrm{\&Sigma;}}}F(i+1,i)F^T(i+1,i).$

（3）状态与过程噪声估计值的互Gramian矩阵

$\left\{\begin{array}{c}{P}_{\mathrm{xw},l}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)=\\ P_{\mathrm{xw},l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)+R_{\mathrm{xyz},l-1}\left(\mathrm{km}\right)R_{\mathrm{eyz},l-1}^{-1}\left(\mathrm{km}\right)R_{\mathrm{wyz},l-1}^T(\mathrm{km},{\mathrm{\&beta;}}_l^k,\mathrm{km}),j>1\\ P_{\mathrm{xw},1}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)=0\end{array}\right.$

；

其中，P_xw,l(km,|(k-1)m)在第k个滤波周期利用第l个到达的测量值进行滤波时，状态与过程噪声估计值的互Gramian矩阵。

（4）过程噪声估计值与测量值的互Gramian矩阵

$R_{\mathrm{wyz},l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k,\mathrm{km})=[R_{\mathrm{wy},l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k,\mathrm{km}),R_{\mathrm{wz},l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k,\mathrm{km})]$

$=R_{\mathrm{xw},l-1}^T(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m){\left[\begin{array}{c}{H}_{l-1}^{*}\left(\mathrm{km}\right)\\ L\left(\mathrm{km}\right)\end{array}\right]}^T-[R_{\mathrm{ww},l-1}^T(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m){\left(H_{l-1}^{*}\left(\mathrm{km}\right)\right)}^T,0]$ ；

其中，R_wyz,l-1(km,,km)在第k个滤波周期利用第l个到达的测量值进行滤波时，过程噪声估计值与测量值的互Gramian矩阵。R_ww,l(km,|(k-1)m)为在第k个滤波周期利用第l个到达的测量值进行滤波时过程噪声估计误差的Gramian矩阵，其计算方法将在（5）过程噪声估计误差的Gramian矩阵中详细介绍。

（5）过程噪声估计误差的Gramian矩阵

$R_{\mathrm{ww},l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k,|(k-1)m)={Q(\mathrm{km},{\mathrm{\&beta;}}_l^k)-P}_{\mathrm{ww},l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k,|(k-1)m)$ ；

其中，P_ww,l(km,|(k-1)m)为在第k个滤波周期利用第l个到达的测量值进行滤波时过程噪声估计值的Gramian矩阵，其计算方法将在（6）过程噪声估计值的Gramian矩阵中详细介绍；Q(km,)为在第k个滤波周期利用第l个到达的测量值进行滤波时过程噪声的Gramian矩阵；

$Q(\mathrm{km},{\mathrm{\&beta;}}_l^k)=\underset{i={\mathrm{\&beta;}}_l^k}{\overset{\mathrm{km}-1}{\mathrm{\&Sigma;}}}F(i+1,i)F^T(i+1,i).$

（6）过程噪声估计值的Gramian矩阵

$\left\{\begin{array}{c}{P}_{\mathrm{ww},l}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m))=P_{\mathrm{ww},l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)\\ +R_{\mathrm{wyz},l-1}^T(\mathrm{km},{\mathrm{\&beta;}}_l^k,\mathrm{km})R_{\mathrm{eyz},l-1}^{-1}\left(\mathrm{km}\right)R_{\mathrm{wyz},l-1}^T(\mathrm{km},{\mathrm{\&beta;}}_l^k,\mathrm{km}),j>1\\ P_{\mathrm{ww},1}(k,{\mathrm{\&beta;}}_l^k|k-1)=0\end{array}\right.$

;

其中，R_ww,l(km,|(k-1)m)为在第k个滤波周期利用第l个到达的测量值进行滤波时过程噪声估计误差的Gramian矩阵；

然后计算出过程噪声估计值为

当l=1时， ${\hat{w}}_1(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)=0$

当l>1时，

${\hat{w}}_l(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)$

${=\hat{w}}_{l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)+R_{\mathrm{wy},l-1}(\mathrm{km},{\mathrm{\&beta;}}_l^k,\mathrm{km}){\left(R_{\mathrm{ey},l-1}\left(\mathrm{km}\right)\right)}^{-1}{e}_{y,l-1}\left(\mathrm{km}\right)$

其中，为在第k个滤波周期利用第l个到达的测量值进行滤波时对过程噪声的估计值；

经过计算，可得到最终伪量测噪声估计值为：

${\hat{v}}_{y,l}^{*}\left(\mathrm{km}\right|(k-1)m)=-H_l^{*}\left(\mathrm{km}\right){\hat{w}}_l(\mathrm{km},{\mathrm{\&beta;}}_l^k|(k-1)m)$

当辅助参量R_ey,l(km)满足R_ey,l(km)>0，且

$R_{\mathrm{ez},j}\left(\mathrm{km}\right)-R_{\mathrm{zy},l}\left(\mathrm{km}\right)R_{\mathrm{ey},l}^{-1}\left(\mathrm{km}\right)R_{\mathrm{yz},l}\left(\mathrm{km}\right)<0$ 时，进入步骤c3，否则，需要重新设定性能指标值γ，然后重新执行步骤c2。

c3：按照以下步骤对测量值进行滤波，实现待估信号估计：

（1）计算得出系统状态预测

当l=1时，

${\hat{x}}_1\left(\mathrm{km}\right|(k-1)m)=F(\mathrm{km},(k-1)m){\hat{x}}_{N_{k-1}}\left((k-1)m\right|(k-1)m);$

当l>1时，

${\hat{x}}_l\left(\mathrm{km}\right|(k-1)m)={\hat{x}}_{l-1}\left(\mathrm{km}\right|\mathrm{km});$

（2）计算得出观测预测

${\hat{y}}_l^{*}\left(\mathrm{km}\right|(k-1)m)=H_l^{*}\left(\mathrm{km}\right){\hat{x}}_l\left(\mathrm{km}\right|(k-1)m)+{\hat{v}}_{y,l}^{*}\left(\mathrm{km}\right|(k-1)m);$

（3）计算得出滤波增益

$K_{y,l}\left(\mathrm{km}\right)=R_{\mathrm{xy},l}\left(\mathrm{km}\right)R_{\mathrm{ey},l}^{-1}\left(\mathrm{km}\right);$

（4）计算得出系统状态估计

${\hat{x}}_l\left(\mathrm{km}\right|\mathrm{km})={\hat{x}}_l\left(\mathrm{km}\right|(k-1)m)+K_{y,l}\left(\mathrm{km}\right)(y_l^{*}\left(\mathrm{km}\right)-{\hat{y}}_l^{*}\left(\mathrm{km}\right|(k-1)m))$ ；

（5）计算得出系统待估信号估计值

$\hat{z}\left(\mathrm{km}\right|\mathrm{km})=L\left(\mathrm{km}\right){\hat{x}}_l\left(\mathrm{km}\right|\mathrm{km});$

当下一到达测量值仍为滤波器在第k个滤波周期内采样所得的测量值，即时，N_k=l+1，则返回步骤c1；

D:当下一到达测量值为滤波器在第k+1个滤波周期内采样所得的测量值，即第k个滤波周期内所有测量值已全部到达融合中心，待估信号估计值在第k个滤波周期基于全局信息的融合估计结果，即 $\hat{z}\left(\mathrm{km}\right|\mathrm{km})=L\left(\mathrm{km}\right){\hat{x}}_{N_k}\left(\mathrm{km}\right|\mathrm{km}).$

下面以由2个传感器构成的多速率多传感器时变离散系统为例，具体阐述本发明所述噪声统计特性未知系统的实时多速率H∞融合滤波方法。

由2个传感器构成的多速率多传感器时变离散系统如下所示

$\left\{\begin{array}{c}x\left(k\right)=F(k,k-1)x(k-1)+w(k,k-1)\\ y_1(2k+1)=H_1(2k+1)x(2k+1)+v_1(2k+1)\\ y_2(3k+2)=H_2(3k+2)x(3k+2)+v_2(3k+2)\\ z\left(3k\right)=L\left(3k\right)x\left(3k\right)\end{array}\right.;$

其中，时变离散系统的演变周期为h，y₁(2k+1)为第一个传感器采集的观测值，其采样周期为2h，该传感器在第1时刻就开始采样；y₂(3k+2)是由第二个传感器采集的测量值，其采样周期为3h，从系统开始后的第2个时刻才开始采样。系统的滤波更新周期为3h，这样就更成了一个典型的简化多速率多传感器融合滤波系统，且系统干扰w(k,k-1),v₁(2k+1),v₂(3k+2)的统计特性未知，但能量在有限域内有界。

在此多速率多传感器时变离散系统中，当每个测量值到达时，融合中心首先在计数器N_i记录该滤波周期内已到达融合中心的测量值个数；

在该仿真算例中，第一个传感器的采样周期为2h，第二个传感器的采样周期为3h，融合中心的滤波更新周期为3h，那么奇数滤波周期内测量信息到达融合中心的情况相同，偶数滤波周期内测量信息到达融合中心的情况也相同，但相邻滤波周期内测量信息到达融合中心的情况各不相同。各传感器的采样时刻与融合中心的滤波更新时刻如图2所示。

在取系统初始值为x₀=[500,1],P₀=[0.1,0;0,0.1]，系统参数为H₁=[5,1],H₂(k)=[1,10],L(k)=[1,0]，F(k,k-1)=[0.95,1;0,0.98]时，首先模拟各传感器的测量输出。然后利用上述滤波方法按照递推循环的方式依次处理融合中心在各时刻获取的测量值，仿真系统待估信号前100时刻的估计曲线与绝对估计误差曲线如图3和图4所示。

下面以一次仿真实现的前两个时刻来说明本发明公布的滤波过程。由于滤波更新周期为3，时刻1与时刻2都属于第1个滤波周期，即k=1

在时刻1，到达融合中心的是第一个传感器得到的测量值，在仿真算例中通过y₁(1)=H₁(1)x(1)+v₁(1)=H₁(1)F(1，0)x₀=[5,1][0.95,1;0,0.98][500,1]+v₁(1)得到，v₁(1)由命令randn(1)生成随机数得到，本实施例中，第一个传感器在时刻1所得测量值的一个随机实现为2386，利用该测量值建立时刻3的伪量测矩阵可通过下式得到 $H_1^{*}\left(3\right)=H_1\left(1\right)F\left(\mathrm{1,3}\right)=\left[\begin{array}{cc}5.5402& -10.0922\end{array}\right],$ 伪量测噪声估计值为： ${\hat{v}}_{y,1}^{*}\left(3\right|0)=-H_1^{*}\left(3\right){\hat{w}}_1\left(\mathrm{3,1}\right|0)=-H_1^{*}\left(3\right)\&times;0=0,$ 利用其对测量值进行滤波，过程如下：

（1）系统状态预测

${\hat{x}}_1\left(3\right|0)=F\left(\mathrm{3,0}\right){\hat{x}}_0=F\left(\mathrm{3,2}\right)F\left(\mathrm{2,1}\right)F\left(\mathrm{1,0}\right){\hat{x}}_0=[431.4814;0.9412]$ ；

（2）观测预测 ${\hat{y}}_1^{*}\left(3\right|0)=H_1^{*}\left(3\right){\hat{x}}_1\left(3\right|0)+{\hat{v}}_{y,1}^{*}\left(3\right|0)=2381$

（3）状态与测量值互Gramian矩阵

$R_{\mathrm{xyz},1}\left(3\right)=\left[\begin{array}{cc}{R}_{\mathrm{xy},1}\left(3\right)& R_{\mathrm{xz},1}\left(3\right)\end{array}\right]=\left[\begin{array}{cc}14.7105& 8.2960\\ -12.8321& 3.0965\end{array}\right]$

（4）测量值的Gramian矩阵

$R_{\mathrm{eyz},1}\left(3\right)=\left[\begin{array}{cc}{R}_{\mathrm{ey},1}\left(3\right)& R_{\mathrm{yz},1}\left(3\right)\\ R_{\mathrm{zy},1}\left(3\right)& R_{\mathrm{ez},1}\left(3\right)\end{array}\right]=\left[\begin{array}{cc}391.1725& 14.7105\\ 14.7105& -0.7040\end{array}\right]$

（5）滤波增益 $K_{y,1}\left(3\right)=R_{\mathrm{xy},1}\left(3\right)R_{\mathrm{ey},1}^{-1}\left(3\right)=[0.0376;-0.0328]$

（6）状态估计：

${\hat{x}}_1\left(3\right|3)={\hat{x}}_1\left(3\right|0)+K_{y,1}\left(3\right)(y_1^{*}\left(3\right)-{\hat{y}}_1^{*}\left(3\right|0))=[431.6560;0.7889]$

（7）待估信号估计值： ${\hat{z}}_1\left(3\right|3)=L\left(3\right){\hat{x}}_1\left(3\right|3)=431.6560$

（8）Riccati变量： $P_1\left(3\right)=\left[\begin{array}{cc}8.2960& 3.0965\\ 3.0965& 2.9714\end{array}\right]$

在时刻2，到达融合中心的是第二个传感器得到的测量值，在仿真算例中通过y₂(2)=H₂(2)x(2)+v₂(2)=H₂(2)F(2,1)x(1)+v₂(2)=[1,10][0.95,1;0,0.98]x(1)+v₂(2)实现，v₂(2)由命令randn(1)生成随机数得到，本实施例中，第二个传感器在时刻2所得测量值的一个随机实现为486.1855，利用该测量值建立时刻3的伪测量矩阵可以通过下式得到 $H_2^{*}\left(3\right)=H_2\left(2\right)F\left(\mathrm{2,3}\right)=\left[\begin{array}{cc}1.0526& 9.13\end{array}\right],$ 伪量测噪声估计值为： ${\hat{v}}_{y,2}^{*}\left(3\right|0)=-H_2^{*}\left(3\right){\hat{w}}_2\left(\mathrm{3,1}\right|0)=0,$ 利用其对测量值进行滤波，过程如下：

（1）系统状态预测 ${\hat{x}}_2\left(3\right|0)={\hat{x}}_1\left(3\right|3)=[431.6560;0.7889];$

（2）观测预测 ${\hat{y}}_2^{*}\left(3\right|3)=H_2^{*}\left(3\right){\hat{x}}_2\left(3\right|0)+{\hat{v}}_{y,2}^{*}\left(3\right|0)=462.7840;$

（3）状态与测量值的互Gramian矩阵

$R_{\mathrm{xyz},2}\left(3\right)=\left[\begin{array}{cc}{R}_{\mathrm{xy},2}\left(3\right)& R_{\mathrm{xz},2}\left(3\right)\end{array}\right]=\left[\begin{array}{cc}293.6780& 54.4287\\ 134.8309& 25.6219\end{array}\right];$

（4）测量值的Gramian矩阵

$R_{\mathrm{eyz},2}\left(3\right)=\left[\begin{array}{cc}{R}_{\mathrm{ey},2}\left(3\right)& R_{\mathrm{yz},2}\left(3\right)\\ R_{\mathrm{zy},2}\left(3\right)& R_{\mathrm{ez},2}\left(3\right)\end{array}\right]=\left[\begin{array}{cc}1701.3& 292.3\\ 292.3& 46.4\end{array}\right];$

（5）滤波增益 $K_{y,2}\left(3\right)=R_{\mathrm{xy},2}\left(3\right)R_{\mathrm{ey},2}^{-1}\left(3\right)=[0.1726;0.0793];$

（6）状态估计：

${\hat{x}}_2\left(3\right|3)={\hat{x}}_2\left(3\right|0)+K_{y,2}\left(3\right)(y_2^{*}\left(3\right)-{\hat{y}}_2^{*}\left(3\right|0))=[435.6956;2.6435];$

（7）待估信号估计值： ${\hat{z}}_2\left(3\right|3)=L\left(3\right){\hat{x}}_2\left(3\right|3)=435.6956;$

（8）Riccati变量： $P_2\left(3\right)=\left[\begin{array}{cc}8.9141& 4.9319\\ 4.9319& 3.6522\end{array}\right].$

Claims (4)

1.一种噪声统计特性未知系统的实时多速率H∞融合滤波方法，其特征在于，包括以下步骤：

A:在由多个传感器构成的多速率多传感器时变离散系统中，当每个测量值到达时，融合中心首先在计数器N_i记录该滤波周期内已到达融合中心的量测个数；

所述的A步骤中，多速率多传感器时变离散系统如下所示：

x(k)＝F(k,k-1)x(k-1)+w(k,k-1)；

y_i(kn_i+j)＝H_i(kn_i+j)x(kn_i+j)+v_i(kn_i+j),i＝1,2,…,p；

z(km)＝L(km)x(km)；

其中，x(k)为演变周期为h的离散系统状态；y_i(kn_i+j)表示第i个传感器第k次采样得到的测量值，其采样起始时刻为j时刻，采样间隔为n_i；kn_i+j是k×n_i+j的缩写形式，表示其采样时刻；z(km)表示滤波周期为m×h(缩写为m)的系统待估信号在k×m×h(缩写为km)时刻的表达式；n_i,m均为正整数；F(k,k-1),H_i(kn_i+j),L(km)分别为系统状态从k-1时刻到k时刻的转移矩阵、传感器i第k次采样时的观测矩阵以及系统在第k个滤波周期的待估信号加权矩阵；w(k,k-1)为系统从k-1时刻到k时刻积累的过程噪声，v_i(kn_i)为第i个传感器第k次采样时的测量噪声，两类噪声的统计特性未知；该多速率多传感器时变离散系统含有p个传感器；

B:利用计数器给出有限域H∞性能指标函数，作为滤波过程的约束指标；

所述的B步骤中，有限域H∞性能指标函数为

$\underset{w,v\&Element;l_2}{sup}\frac{\underset{N_i\&NotEqual;0}{\overset{k-1}{\underset{i=1}{\&Sigma;}}}\underset{j=1}{\overset{N_i}{\&Sigma;}}{e}_{z,j}^T\left(im\right)e_{z,j}\left(im\right)+\underset{j=1}{\overset{l}{\&Sigma;}}{e}_{z,j}^T\left(km\right)e_{z,j}\left(km\right)}{\underset{N_i\&NotEqual;0}{\overset{k-1}{\underset{i=1}{\&Sigma;}}}\underset{j=1}{\overset{N_i}{\&Sigma;}}{v}_{{\mathrm{\&alpha;}}_j^i}^T\left({\mathrm{\&beta;}}_j^i\right)v_{{\mathrm{\&alpha;}}_j^i}\left({\mathrm{\&beta;}}_j^i\right)+\underset{j=1}{\overset{l}{\&Sigma;}}{v}_{{\mathrm{\&alpha;}}_j^k}^T\left({\mathrm{\&beta;}}_j^k\right)v_{{\mathrm{\&alpha;}}_j^k}\left({\mathrm{\&beta;}}_j^k\right)+\underset{i=1}{\overset{mk}{\&Sigma;}}{w}^T(i,i-1)w(i,i-1)+{\left(x\right(0)-{\hat{x}}_0)}^T{P}_0^{-1}\left(x\right(0)-{\hat{x}}_0)}<{\mathrm{\&gamma;}}^2;$

其中，e_z,j(ik)表示融合中心在第i个滤波周期利用第j个到达的测量值对待估信号进行滤波时得到的估计误差，γ为给定的有限域H∞性能指标值；表示在第i个滤波周期第j个到达的测量值对应的测量噪声，表示采集该测量信息的传感器标号，表示该测量信息的采集时刻；e_z,j(im)表示融合中心在第i个滤波周期利用第j个到达的测量值对待估信号z(im)进行滤波时得到的估计误差，分别表示初始状态及其估计值，P₀表示偏离x(0)的程度；

C:在上述性能指标的约束下，进行测量值的融合滤波；

所述的C步骤中，进行测量值的融合滤波包括以下步骤：

c1:利用测量值建立当前滤波时刻的伪量测矩阵；

所述的c1步骤中，

对应的伪量测噪声为： $v_{y,l}^{*}\left(km\right)=-H_l^{*}\left(km\right)w(km,{\mathrm{\&beta;}}_l^k);$

对应的伪量测矩阵为： $H_l^{*}\left(km\right)=H_{{\mathrm{\&alpha;}}_l^k}\left({\mathrm{\&beta;}}_l^k\right)F({\mathrm{\&beta;}}_l^k,km);$

对应的伪测量信息为： $y_l^{*}\left(km\right)=y_{{\mathrm{\&alpha;}}_l^k}\left({\mathrm{\&beta;}}_l^k\right)=H_l^{*}\left(km\right)x\left(km\right)+v_{y,l}^{*}\left(km\right);$

其中，对第k个滤波周期内利用第N_k＝l个到达融合中心的测量值进行滤波，表示在第i个滤波周期内到达融合中心的第j个量测，第j个量测是由第个传感器在时刻采样得到的；表示时刻到km时刻积累的过程噪声；表示时刻到km时刻系统状态转移矩阵；为其(伪)逆矩阵；

c2：估计系统噪声及其与系统状态之间的耦合关系；

c3：对测量值进行滤波实现待估信号估计；

D:当该滤波周期内所有测量值都到达融合中心时，即可得到系统待估信号基于全局信息的融合估计结果。

2.根据权利要求1所述的噪声统计特性未知系统的实时多速率H∞融合滤波方法，其特征在于：所述的c2步骤中，在估计系统噪声及其与系统状态之间的耦合关系时，首先求出下列辅助参量：

(1)求解状态与测量预测值的互Gramian矩阵

$\begin{array}{c}{R}_{xyz,l}\left(km\right)=\left[\begin{array}{cc}{R}_{xy,l}\left(km\right)& R_{xz,l}\left(km\right)\end{array}\right]\\ =P_l\left(km\right)\&lsqb;{\left(H_l^{*}\right(km\left)\right)}^T,L^T\left(km\right)\&rsqb;-\&lsqb;\left(\underset{t={\mathrm{\&beta;}}_l^k}{\overset{km}{\&Sigma;}}F\right(t,t-1\left)F^T\right(t,t-1)-P_{xw,l}(km,{\mathrm{\&beta;}}_l^k|\left(k-1\right)m\left)\right){\left(H_l^{*}\right(km\left)\right)}^T,0\&rsqb;\end{array};$

其中，为状态与过程噪声的互Gramian矩阵；P_l(km)为在第k个滤波周期利用第l个到达的测量值进行滤波时的Riccati变量，其满足如下递推关系：

$\left\{\begin{array}{c}{P}_{l+1}\left(km\right)=P_l\left(km\right)-K_l\left(km\right)R_{xyz,l}^T\left(km\right)\\ P_1\left(\right(k+1\left)m\right)=F\left(\right(k+1)m,km)\left(P_{N_k}\right(km)-K_{N_k}(k\left)R_{xy,N_k}^T\right(k\left)\right)F^T\left(\right(k+1)m,km)+Q\left(\right(k+1)m,km)\end{array}\right.;$

其中，Q((k+1)m,km)为km时刻到(k+1)m时刻积累的过程噪声的Gramian矩阵； $Q\left(\right(k+1)m,km)=\underset{i=km}{\overset{km+m-1}{\&Sigma;}}F(i+1,i)F^T(i+1,i);$ 表示在第k个滤波周期利用第N_k个到达的测量值进行滤波时的Riccati变量；为在第k个滤波周期利用第N_k个到达的测量值进行滤波时的滤波增益矩阵；为在第k个滤波周期利用第N_k个到达的测量值进行滤波时的系统状态与测量预测之间的互Gramian矩阵；

(2)求解测量预测值的Gramian矩阵

$\begin{array}{c}{R}_{eyz,l}\left(km\right)=\left[\begin{array}{cc}{R}_{ey,l}\left(km\right)& R_{yz,l}\left(km\right)\\ R_{zy,l}\left(km\right)& R_{ez,l}\left(km\right)\end{array}\right]\\ =\left[\begin{array}{c}{H}_l^{*}\left(km\right)\\ L\left(km\right)\end{array}\right]P_l\left(km\right)\&lsqb;{\left(H_l^{*}\right(km\left)\right)}^T,L^T\left(km\right)\&rsqb;\\ +\left[\begin{array}{cc}{H}_l^{*}\left(km\right)\left(Q\right(km,{\mathrm{\&beta;}}_l^k)-R_{ww,l}(km,{\mathrm{\&beta;}}_l^k|\left(k-1\right)m\left)\right){\left(H_l^{*}\right(km\left)\right)}^T& 0\\ 0& -{\mathrm{\&gamma;}}^{2}I\end{array}\right]\\ -\left[\begin{array}{cc}{H}_l^{*}\left(km\right)R_{xw,l}\left(km,{\mathrm{\&beta;}}_l^k|\left(k-1\right)m\right){\left(H_l^{*}\right(km\left)\right)}^T& 0\\ 0& 0\end{array}\right]\\ -{\left[\begin{array}{cc}{H}_l^{*}\left(km\right)R_{xw,l}\left(km,{\mathrm{\&beta;}}_l^k|\left(k-1\right)m\right){\left(H_l^{*}\right(km\left)\right)}^T& 0\\ 0& 0\end{array}\right]}^T\end{array};$

其中，R_eyz,l(km)表示在第k个滤波周期利用第l个到达的测量值进行滤波时的测量预测值的Gramian矩阵；为在第k个滤波周期利用第l个到达的测量值进行滤波时状态与过程噪声估计值的互Gramian矩阵；为在第k个滤波周期利用第l个到达的测量值进行滤波时过程噪声估计误差的Gramian矩阵；为在第k个滤波周期利用第l个到达的测量值进行滤波时对应耦合过程噪声的Gramian矩阵；

$Q(km,{\mathrm{\&beta;}}_l^k)=\underset{i={\mathrm{\&beta;}}_l^k}{\overset{km-1}{\&Sigma;}}F(i+1,i)F^T(i+1,i);$

(3)状态与过程噪声估计值的互Gramian矩阵

$\{\begin{array}{c}{P}_{xw,l}(km,{\mathrm{\&beta;}}_l^k|\left(k-1\right)m)=P_{xw,l-1}(km,{\mathrm{\&beta;}}_l^k|\left(k-1\right)m)\\ +R_{xyz,l-1}\left(km\right)R_{eyz,l-1}^{-1}\left(km\right)R_{wyz,l-1}^T(km,{\mathrm{\&beta;}}_l^k,km)\\ P_{xw,1}(km,{\mathrm{\&beta;}}_l^k|\left(k-1\right)m)=0\end{array},j>1;$

其中，在第k个滤波周期利用第l个到达的测量值进行滤波时，状态与过程噪声估计值的互Gramian矩阵；

(4)过程噪声估计值与测量值的互Gramian矩阵

$\begin{array}{c}{R}_{wyz,l-1}(km,{\mathrm{\&beta;}}_l^k,km)\&lsqb;R_{wy,l-1}(km,{\mathrm{\&beta;}}_l^k,km),R_{wz,l-1}(km,{\mathrm{\&beta;}}_l^k,km)\&rsqb;\\ =R_{xw,l-1}^T(km,{\mathrm{\&beta;}}_l^k|\left(k-1\right)m){\left[\begin{array}{c}{H}_{l-1}^{*}\left(km\right)\\ L\left(km\right)\end{array}\right]}^T-\&lsqb;R_{ww,l-1}^T(km,{\mathrm{\&beta;}}_l^k|\left(k-1\right)m){\left(H_{l-1}^{*}\right(km\left)\right)}^T,0\&rsqb;\end{array};$

其中，在第k个滤波周期利用第l个到达的测量值进行滤波时，过程噪声估计值与测量值的互Gramian矩阵；为在第k个滤波周期利用第l个到达的测量值进行滤波时过程噪声估计误差的Gramian矩阵；

(5)过程噪声估计误差的Gramian矩阵

$R_{ww,l-1}(km,{\mathrm{\&beta;}}_l^k|\left(k-1\right)m)=Q(km,{\mathrm{\&beta;}}_l^k)-P_{ww,l-1}(km,{\mathrm{\&beta;}}_l^k|\left(k-1\right)m);$

其中，为在第k个滤波周期利用第l个到达的测量值进行滤波时过程噪声估计值的Gramian矩阵；为在第k个滤波周期利用第l个到达的测量值进行滤波时过程噪声的Gramian矩阵；

$Q(km,{\mathrm{\&beta;}}_l^k)=\underset{i={\mathrm{\&beta;}}_l^k}{\overset{km-1}{\&Sigma;}}F(i+1,i)F^T(i+1,i);$

(6)过程噪声估计值的Gramian矩阵

$\left\{\begin{array}{c}{P}_{ww,l}(km,{\mathrm{\&beta;}}_l^k|\left(k-1\right)m))=P_{ww,l-1}(km,{\mathrm{\&beta;}}_l^k|\left(k-1\right)m)\\ +R_{wyz,l-1}^T(km,{\mathrm{\&beta;}}_l^k,km)R_{eyz,l-1}^{-1}\left(km\right)R_{wyz,l-1}^T(km,{\mathrm{\&beta;}}_l^k,km)\\ P_{ww,l}(k,{\mathrm{\&beta;}}_l^k|k-1)=0\end{array}\right.,j>1;$

其中，为在第k个滤波周期利用第l个到达的测量值进行滤波时过程噪声估计误差的Gramian矩阵；

然后计算出过程噪声估计值为

当l＝1时， ${\hat{w}}_1(km,{\mathrm{\&beta;}}_l^k|\left(k-1\right)m)=0;$

当l＞1时，

$\begin{array}{c}{\hat{w}}_l(km,{\mathrm{\&beta;}}_l^k|\left(k-1\right)m)\\ ={\hat{w}}_{l-1}(km,{\mathrm{\&beta;}}_l^k|\left(k-1\right)m)+R_{wy,l-1}(km,{\mathrm{\&beta;}}_l^k,m){\left(R_{ey,l-1}\right(km\left)\right)}^{-1}{e}_{y,l-1}\left(km\right)\end{array};$

其中，为在第k个滤波周期利用第l个到达的测量值进行滤波时对过程噪声的估计值；

最终得到伪量测噪声估计值为：

${\hat{v}}_{y,l}^{*}\left(km\right|\left(k-1\right)m)=-H_l^{*}\left(km\right){\hat{w}}_l(km,{\mathrm{\&beta;}}_l^k|\left(k-1\right)m);$

当辅助参量R_ey,l(km)满足R_ey,l(km)＞0，且

$R_{ez,j}\left(km\right)-R_{zy,l}\left(km\right)R_{ey,l}^{-1}\left(km\right)R_{yz,l}\left(km\right)<0$ 时，进入步骤c3，否则，需要重新设定性能指标值γ，然后重新执行步骤c2。

3.根据权利要求2所述的噪声统计特性未知系统的实时多速率H∞融合滤波方法，其特征在于：所述的c3步骤中，利用伪量测噪声估计值对测量值进行滤波包含以下步骤：

(1)计算得出系统状态预测

当l＝1时，

${\hat{x}}_1\left(km\right|\left(k-1\right)m)=F(km,(k-1\left)m\right){\hat{x}}_{N_{k-1}}\left(\right(k-1\left)m\right|\left(k-1\right)m);$

当l＞1时，

${\hat{x}}_l\left(km\right|\left(k-1\right)m)={\hat{x}}_{l-1}\left(km\right|km);$

(2)计算得出观测预测

${\hat{y}}_l^{*}\left(km\right|\left(k-1\right)m)=H_l^{*}\left(km\right){\hat{x}}_l\left(km\right|\left(k-1\right)m)+{\hat{v}}_{y,l}^{*}\left(km\right|\left(k-1\right)m);$

(3)计算得出滤波增益

$K_{y,l}\left(km\right)=R_{xy,l}\left(km\right)R_{ey,l}^{-1}\left(km\right);$

(4)计算得出系统状态估计

${\hat{x}}_l\left(km\right|km)={\hat{x}}_l\left(km\right|\left(k-1\right)m)+K_{y,l}\left(km\right)\left(y_l^{*}\right(km)-{\hat{y}}_l^{*}(km|\left(k-1\right)m\left)\right);$

(5)计算得出系统待估信号估计值

$\hat{z}\left(km\right|km)=L\left(km\right){\hat{x}}_l\left(km\right|km);$

当下一到达测量值仍为滤波器在第k个滤波周期内采样所得的测量值，即时，N_k＝l+1，则返回步骤c1。

4.根据权利要求3所述的噪声统计特性未知系统的实时多速率H∞融合滤波方法，其特征在于：所述的D步骤中，当该滤波周期内所有测量值都到达融合中心，可得到系统待估信号基于全局信息的融合估计结果： $\hat{z}\left(km\right|km)=L\left(km\right){\hat{x}}_{N_k}\left(km\right|km).$