CN102749096B - 一种对双观测系统量测噪声方差阵的自适应同步估计方法 - Google Patents

Abstract

本发明提出一种对双观测系统量测噪声方差阵的自适应同步估计方法，属于信号处理技术领域，该方法包括：得到两观测系统在各时刻的数据信号；计算两观测系统数据信号Z1(i)的自差分序列以及两个自差分序列的互差分序列；计算两观测系统数据信号在不同窗口长度下自差分序列的方差以及两个观测系统在不同窗口长度下互差分序列的方差，得到两观测系统不同窗口长度下的量测噪声方差阵；使用最优窗口长度计算得当前时刻的量测噪声方差阵。本发明利用双测量系统自差分序列和互差分序列信息，能够同时有效估计出两个测量系统的量测噪声方差阵，可实时调节计算噪声方差阵时的数据窗口长度，提高估计精度。本发明单纯利用测量系统信息，避免了误差耦合，精度较高。

Description

一种对双观测系统量测噪声方差阵的自适应同步估计方法

技术领域

本发明属于信号处理技术领域，具体涉及一种对双观测系统量测噪声方差阵的自适应同步估计方法，可同时有效自适应估计出不同测量系统的量测噪声方差阵。

背景技术

卡尔曼滤波算法是一项广泛应用的信息融合技术，标准卡尔曼滤波滤波算法是建立在系统噪声方差阵已知情况下的最优估计，但在实际情况下，测量系统的测量噪声方差阵未知，如何有效估计出系统测量噪声方差阵对提高滤波精度具有重要意义。

通常，若要获得系统测量噪声方差阵一般通过对测量系统进行长时间观测，根据对系统输出信息的统计量获得量测噪声方差阵，然而这种基于验前信息获得噪声统计特性的方法无法解决时变噪声的问题。在自适应估计系统量测噪声的研究中，方法主要为基于信息序列的自适应估计方法（IAE），主要代表为sage-husa法，它可以自适应估计系统的状态噪声方差阵和量测噪声方差阵，然而此法的信息序列基于状态一步递推值和观测值，当状态估计值计算不准确时，会造成误差耦合，导致量测噪声阵的估计精度下降。本发明利用两种测量参考系统的不同测量特性，自适应调节窗口长度，对系统量测噪声进行有效估计。

发明内容

针对现有技术中存在的问题，本发明为解决对测量系统量测噪声方差阵进行有效估计的问题，提出了一种对双观测系统量测噪声方差阵的自适应同步估计方法。该方法通过计算双观测系统观测值的自差分序列以及互差分序列的方差，得到量测噪声方差阵，并且对同一时刻不同窗口长度下计算得到的噪声方差阵序列进行一次直线拟合，由拟合得到的直线斜率计算出最优窗口长度并得到最终的量测噪声方差阵。

本发明提出一种对双观测系统量测噪声方差阵的自适应同步估计方法，包括以下几个步骤：

步骤一：分别得到观测系统A、观测系统B在各时刻的数据信号Z1(i)、Z2(i)，其中i为观测系统的数据信号测量时刻；

步骤二：分别计算观测系统A数据信号Z1(i)的自差分序列和观测系统B数据信号Z2(i)的自差分序列，并计算两个自差分序列的互差分序列，具体为：

(1)观测系统A数据信号Z1(i)的自差分序列ΔZ1(i)为：

ΔZ1(i)＝Z1(i)-Z1(i-1)

其中i-1和i分别表示观测系统A的数据采集时刻；

(2)观测系统B数据信号Z2(i)的自差分序列ΔZ2(i)为：

ΔZ2(i)＝Z2(i)-Z2(i-1)

其中i-1和i分别表示观测系统B的数据采集时刻；

(3)观测系统A和观测系统B的两个自差分序列的的互差分序列C(i)为：

C(i)＝ΔZ1(i)-ΔZ2(i)

步骤三：分别计算观测系统A数据信号在不同窗口长度下自差分序列的方差，观测系统B数据信号在不同窗口长度下自差分序列的方差，以及两个观测系统在不同窗口长度下互差分序列的方差，并利用方差进行相关计算得到观测系统A、观测系统B不同窗口长度下的量测噪声方差阵R1、R2，具体为：

(1)首先选取最大窗口长度M_Max、最小窗口长度M_Min和窗口长度间隔M_Dis，且窗口长度间隔M_Dis为M_Max与M_Min之差的整数倍，若当前数据信号量测时刻k小于等于M_Max，则直接采用窗口长度累积的方式计算观测系统A自差分序列的方差：

$\left\{\begin{array}{c}{E}_1\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\mathrm{\ΔZ}1\left(i\right)\\ {\mathrm{\σ}}_1\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[\mathrm{\ΔZ}1\left(i\right)-E_1\left(k\right)]}^2\end{array}\right.$

其中E₁(k)表示观测系统A的1到k时刻观测值序列的均值，σ₁(k)表示观测系统A的1到k时刻观测值序列的方差；

若当前数据信号量测时刻k大于M_Max，则计算观测系统A数据信号在不同窗口长度下自差分序列的方差：

$\left\{\begin{array}{c}{E}_1(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}\mathrm{\ΔZ}1\left(i\right)\\ {\mathrm{\σ}}_1(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}{[\mathrm{\ΔZ}1\left(i\right)-E_1(M,k)]}^2\end{array}\right.M=M_{\mathrm{Min}},M_{\mathrm{Min}}+M_{\mathrm{Dis}}...M_{\mathrm{Max}}-M_{\mathrm{Dis}},M_{\mathrm{Max}}$

其中，k为观测系统A当前数据信号量测时刻，i为数据信号量测时刻，ΔZ1(i)为观测系统A自差分序列在i时刻的值，E₁(M,k)为k时刻窗口长度为M时观测系统A自差分序列的均值，σ₁(M,k)为k时刻窗口长度为M时自差分序列的方差，M为窗口长度序列；

(2)若当前数据量测时刻k小于等于M_Max，则直接采用窗口长度累积的方式计算观测系统B自差分序列的方差：

$\left\{\begin{array}{c}{E}_2\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\mathrm{\ΔZ}2\left(i\right)\\ {\mathrm{\σ}}_2\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[\mathrm{\ΔZ}2\left(i\right)-E_2\left(k\right)]}^2\end{array}\right.$

其中E₂(k)表示观测系统B的1到k时刻观测值序列的均值，σ₂(k)表示观测系统B的1到k时刻观测值序列的均值；

若当前数据量测时刻k大于M_Max，则计算观测系统B数据信号在不同窗口长度下的自差分序列的方差：

$\left\{\begin{array}{c}{E}_2(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}\mathrm{\ΔZ}2\left(i\right)\\ {\mathrm{\σ}}_2(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}{[\mathrm{\ΔZ}2\left(i\right)-E_2(M,k)]}^2\end{array}\right.M=M_{\mathrm{Min}},M_{\mathrm{Min}}+M_{\mathrm{Dis}}...M_{\mathrm{Max}}-M_{\mathrm{Dis}},M_{\mathrm{Max}}$

其中，k为观测系统B当前数据信号量测时刻，i为数据信号量测时刻，ΔZ2(i)为观测系统B自差分序列在i时刻的值，E₂(M,k)为k时刻窗口长度为M时观测系统B自差分序列的均值，σ₂(M,k)为k时刻窗口长度为M时的自差分序列方差，M为窗口长度序列；

(3)若当前数据量测时刻k小于等于M_Max，计算两个观测系统互差分序列的方差σ_C(k)：

$\left\{\begin{array}{c}C\left(i\right)=\mathrm{\ΔZ}1\left(i\right)-\mathrm{\ΔZ}2\left(i\right)\\ E_C\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}C\left(i\right)\\ {\mathrm{\σ}}_C\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[C\left(i\right)-E_C\left(k\right)]}^2\end{array}\right.$

其中，C(i)表示观测系统A和观测系统B自差分序列在i时刻的互差分序列；E_C(k)表示k时刻互差分序列的均值；

若当前系统数据量测时刻k大于M_Max，则按下式计算两个量测系统互差分序列的方差：

$\left\{\begin{array}{c}C\left(i\right)=\mathrm{\ΔZ}1\left(i\right)-\mathrm{\ΔZ}2\left(i\right)\\ E_C(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}C\left(i\right)\\ {\mathrm{\σ}}_C(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}{[C\left(i\right)-E_C(M,k)]}^2\end{array}\right.M=M_{\mathrm{Min}},M_{\mathrm{Min}}+M_{\mathrm{Dis}}...M_{\mathrm{Max}}-M_{\mathrm{Dis}},M_{\mathrm{Max}}$

其中，k为测量系统当前数据信号量测时刻，i为数据信号量测时刻，C(i)表示观测系统A和观测系统B自差分序列在i时刻的互差分序列，E_C(M,k)为k时刻窗口长度为M时互差分序列的均值，σ_C(M,k)为k时刻窗口长度为M时互差分序列的方差，M为窗口长度序列；

（4）利用上述（1）、（2）、（3）步骤中得到的各个方差计算观测系统A和观测系统B的量测噪声方差阵；

若当前数据信号量测时刻k小于等于M_Max，量测噪声方差阵为：

$\left\{\begin{array}{c}{R}_1\left(k\right)=\frac{{\mathrm{\σ}}_c\left(k\right)+({\mathrm{\σ}}_1\left(k\right)-{\mathrm{\σ}}_2\left(k\right))}{4}\\ R_2\left(k\right)=\frac{{\mathrm{\σ}}_c\left(k\right)-({\mathrm{\σ}}_1\left(k\right)-{\mathrm{\σ}}_2\left(k\right))}{4}\end{array}\right.$

R₁(k)表示数据量测时刻k小于M_Max时观测系统A的量测噪声方差阵；R₂(k)表示数据量测时刻k小于M_Max时观测系统B的量测噪声方差阵；

若当前数据信号量测时刻k大于M_Max，量测噪声方差阵为：

$\left\{\begin{array}{c}{R}_1(M,k)=\frac{{\mathrm{\σ}}_c(M,k)+({\mathrm{\σ}}_1(M,k)-{\mathrm{\σ}}_2(M,k))}{4}\\ R_2(M,k)=\frac{{\mathrm{\σ}}_c(M,k)-({\mathrm{\σ}}_1(M,k)-{\mathrm{\σ}}_2(M,k))}{4}\end{array}\right.M=M_{\mathrm{Min}},M_{\mathrm{Min}}+M_{\mathrm{Dis}}...M_{\mathrm{Max}}-M_{\mathrm{Dis}},M_{\mathrm{Max}}$

其中，k为系统当前数据量测时刻，σ₁(k)、σ₂(k)分别为观测系统A、观测系统B自差分序列在k时刻的方差，σ_C(k)为互差分序列在k时刻的方差，M为窗口长度，R₁(M,k)表示数据量测时刻k大于M_Max时观测系统A的量测噪声方差阵；R₂(M,k)表示数据量测时刻k大于M_Max时观测系统B的量测噪声方差阵；

步骤四：根据两个观测系统在不同窗口长度下计算得到的量测噪声方差阵信息，得到当前时刻最优窗口长度，使用最优窗口长度计算得当前时刻的量测噪声方差阵，具体为：

（1）当k小于等于M_Max时，直接按照步骤三（4）中k小于等于M_Max情况下的公式计算量测噪声方差阵：

$\left\{\begin{array}{c}{R}_1\left(k\right)=\frac{{\mathrm{\σ}}_c\left(k\right)+({\mathrm{\σ}}_1\left(k\right)-{\mathrm{\σ}}_2\left(k\right))}{4}\\ R_2\left(k\right)=\frac{{\mathrm{\σ}}_c\left(k\right)-({\mathrm{\σ}}_1\left(k\right)-{\mathrm{\σ}}_2\left(k\right))}{4}\end{array}\right.$

（2）当k大于M_Max时，将k时刻不同窗口下计算得到的噪声方差阵序列使用一次线形函数进行拟合，并记录其斜率；

$\left\{\begin{array}{c}R(M_{\mathrm{Min}},k)=a*(\frac{M_{\mathrm{Min}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1)+b\\ R(M_{\mathrm{Min}}+M_{\mathrm{Dis}},k)=a*(\frac{M_{\mathrm{Min}}+M_{\mathrm{Dis}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1)+b\\ \begin{array}{cc}\begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\end{array}\\ R(M_{\mathrm{Max}}-M_{\mathrm{Dis}},k)=a*(\frac{M_{\mathrm{Max}}-M_{\mathrm{Dis}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1)+b\\ R(M_{\mathrm{Max}},k)=a*(\frac{M_{\mathrm{Max}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1)+b\end{array}\right.$

$\left[\begin{array}{cc}\frac{M_{\mathrm{Min}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1& 1\\ \frac{M_{\mathrm{Min}}+M_{\mathrm{Dis}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1& 1\\ .& .\\ .& .\\ .& .\\ \frac{M_{\mathrm{Max}}-M_{\mathrm{Dis}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1& 1\\ \frac{M_{\mathrm{Max}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1& 1\end{array}\right]{[a,b]}^T=\left[\begin{array}{c}R(M_{\mathrm{Min}},k)\\ R(M_{\mathrm{Min}}+M_{\mathrm{Dis}},k)\\ .\\ .\\ .\\ R(M_{\mathrm{Max}}-M_{\mathrm{Dis}},k)\\ R(M_{\mathrm{Max}},k)\end{array}\right]$

其中，k为当前时刻，a为拟合直线的斜率，b为拟合直线的截距：

则拟合直线的斜率和截距为：

其中R(M_Min，k)、R(M_Min+M_Dis,k)、R(M_Max-M_Dis,k)、R(M_Max,k)分别表示窗口长度为M_Min，(M_Min+M_Dis),(M_Max-M_Dis)，M_Max时的量测噪声方差阵；

（3）将拟合直线计算的斜率带入下式方程得到当前时刻k应使用的数据窗口长度M'，若计算得到的窗口长度M'小于M_Min，则取M'为M_Min，且M_Min取值为100：

其中：fix(·)表示取整，a_max为所拟合的直线达到的最大斜率；

（4）利用最优窗口内的序列信息计算得到最终噪声方差阵的估计值：

$\left\{\begin{array}{c}{R}_{1f}(M,k)=\frac{{\mathrm{\σ}}_c(M^{\′},k)+({\mathrm{\σ}}_1(M^{\′},k)-{\mathrm{\σ}}_2(M^{\′},k))}{4}\\ R_{2f}(M,k)=\frac{{\mathrm{\σ}}_c(M^{\′},k)-({\mathrm{\σ}}_1(M^{\′},k)-{\mathrm{\σ}}_2(M^{\′},k))}{4}\end{array}\right.$

R_1f(M,k)、R_2f(M,k)分别表示观测系统A、观测系统B在k时刻估计得到的量测噪声方差阵。

通过本发明提出的自适应同步估计方法，通过计算两个测量系统自差分序列及互差分序列的方差，计算不同窗口长度下的量测噪声方差阵，得到最优窗口长度，从而计算出最终的噪声方差阵。

本发明的优点在于：

（1）本发明提出了一种对双观测系统量测噪声方差阵的自适应同步估计方法，利用双测量系统自差分序列和互差分序列信息，能够同时有效估计出两个测量系统的量测噪声方差阵。

（2）本发明提出了一种对双观测系统量测噪声方差阵的自适应同步估计方法，可实时调节计算噪声方差阵时的数据窗口长度，提高估计精度。

（3）本发明提出了一种对双观测系统量测噪声方差阵的自适应同步估计方法，单纯利用测量系统信息，避免了误差耦合，精度较高。

附图说明

图1是本发明一种对双观测系统量测噪声方差阵的自适应同步估计方法的流程图；

图2是在测量噪声方差不变情况下，对两个量测系统测量噪声方差阵的估计情况；

图3a在测量噪声法阵为时变情况下，使用固定窗口对量测系统测量噪声方差阵的估计情况；

图3b在测量噪声法阵为时变情况下，使用自适应调节窗口对量测系统测量噪声方差阵的估计情况；

具体实施方式

下面将结合附图对本发明作进一步的详细说明。

本发明提出一种对双观测系统量测噪声方差阵的自适应同步估计方法，如图1所示，包括以下几个步骤：

步骤一：分别得到观测系统A、观测系统B在各时刻的数据信号Z1(i)、Z2(i)，其中i为观测系统的数据信号测量时刻。

步骤二：分别计算观测系统A数据信号Z1(i)的自差分序列和观测系统B数据信号Z2(i)的自差分序列，并计算两个自差分序列的互差分序列，具体包括记下几个步骤：

(1)观测系统A数据信号Z1(i)的自差分序列ΔZ1(i)为：

ΔZ1(i)＝Z1(i)-Z1(i-1)

其中i-1和i分别表示观测系统A的数据采集时刻。

(2)观测系统B数据信号Z2(i)的自差分序列ΔZ2(i)为：

ΔZ2(i)＝Z2(i)-Z2(i-1)

其中i-1和i分别表示观测系统B的数据采集时刻。

(3)观测系统A和观测系统B的两个自差分序列的的互差分序列C(i)为：

C(i)＝ΔZ1(i)-ΔZ2(i)

步骤三：分别计算观测系统A数据信号在不同窗口长度下自差分序列的方差，观测系统B数据信号在不同窗口长度下自差分序列的方差，以及两个观测系统在不同窗口长度下互差分序列的方差，并利用方差进行相关计算得到观测系统A、观测系统B不同窗口长度下的量测噪声方差阵R1、R2，具体包括以下几个步骤：

(1)首先选取一个最大窗口长度M_Max，最小窗口长度M_Min，窗口长度间隔M_Dis，且窗口长度间隔为M_Max与M_Min之差的整数倍，若当前数据信号量测时刻k小于等于M_Max，则直接采用窗口长度累积的方式计算观测系统A自差分序列的方差：

$\left\{\begin{array}{c}{E}_1\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\mathrm{\ΔZ}1\left(i\right)\\ {\mathrm{\σ}}_1\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[\mathrm{\ΔZ}1\left(i\right)-E_1\left(k\right)]}^2\end{array}\right.$

其中E₁(k)表示观测系统A的1到k时刻观测值序列的均值，σ₁(k)表示观测系统A的1到k时刻观测值序列的方差。

若当前数据信号量测时刻k大于M_Max，则计算观测系统A数据信号在不同窗口长度下自差分序列的方差：

$\left\{\begin{array}{c}{E}_1(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}\mathrm{\ΔZ}1\left(i\right)\\ {\mathrm{\σ}}_1(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}{[\mathrm{\ΔZ}1\left(i\right)-E_1(M,k)]}^2\end{array}\right.M=M_{\mathrm{Min}},M_{\mathrm{Min}}+M_{\mathrm{Dis}}...M_{\mathrm{Max}}-M_{\mathrm{Dis}},M_{\mathrm{Max}}$

其中，k为观测系统A当前数据信号量测时刻，i为数据信号量测时刻，ΔZ1(i)为观测系统A自差分序列在i时刻的值，E₁(M,k)为k时刻窗口长度为M时观测系统A自差分序列的均值，σ₁(M,k)为k时刻窗口长度为M时自差分序列的方差，M为窗口长度序列。

(2)若当前数据量测时刻k小于等于M_Max，则直接采用窗口长度累积的方式计算观测系统B自差分序列的方差：

$\left\{\begin{array}{c}{E}_2\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\mathrm{\ΔZ}2\left(i\right)\\ {\mathrm{\σ}}_2\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[\mathrm{\ΔZ}2\left(i\right)-E_2\left(k\right)]}^2\end{array}\right.$

其中E₂(k)表示观测系统B的1到k时刻观测值序列的均值，σ₂(k)表示观测系统B的1到k时刻观测值序列的均值。

若当前数据量测时刻k大于M_Max，则计算观测系统B数据信号在不同窗口长度下的自差分序列的方差：

$\left\{\begin{array}{c}{E}_2(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}\mathrm{\ΔZ}2\left(i\right)\\ {\mathrm{\σ}}_2(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}{[\mathrm{\ΔZ}2\left(i\right)-E_2(M,k)]}^2\end{array}\right.M=M_{\mathrm{Min}},M_{\mathrm{Min}}+M_{\mathrm{Dis}}...M_{\mathrm{Max}}-M_{\mathrm{Dis}},M_{\mathrm{Max}}$

其中，k为观测系统B当前数据信号量测时刻，i为数据信号量测时刻，ΔZ2(i)为观测系统B自差分序列在i时刻的值，E₂(M,k)为k时刻窗口长度为M时观测系统B自差分序列的均值，σ₂(M,k)为k时刻窗口长度为M时的自差分序列方差，M为窗口长度序列。

(3)若当前数据量测时刻k小于等于M_Max，计算两个观测系统互差分序列的方差σ_C(k)：

$\left\{\begin{array}{c}C\left(i\right)=\mathrm{\ΔZ}1\left(i\right)-\mathrm{\ΔZ}2\left(i\right)\\ E_C\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}C\left(i\right)\\ {\mathrm{\σ}}_C\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[C\left(i\right)-E_C\left(k\right)]}^2\end{array}\right.$

其中，C(i)表示观测系统A和观测系统B自差分序列的互差分序列；E_C(k)表示k时刻互差分序列的均值。

若当前系统数据量测时刻k大于M_Max，则按下式计算两个量测系统互差分序列的方差：

$\left\{\begin{array}{c}C\left(i\right)=\mathrm{\ΔZ}1\left(i\right)-\mathrm{\ΔZ}2\left(i\right)\\ E_C(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}C\left(i\right)\\ {\mathrm{\σ}}_C(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}{[C\left(i\right)-E_C(M,k)]}^2\end{array}\right.M=M_{\mathrm{Min}},M_{\mathrm{Min}}+M_{\mathrm{Dis}}...M_{\mathrm{Max}}-M_{\mathrm{Dis}},M_{\mathrm{Max}}$

其中，k为测量系统当前数据信号量测时刻，i为数据信号量测时刻，C(i)为互差分序列在i时刻的值，E_C(M,k)为k时刻窗口长度为M时互差分序列的均值，σ_C(M,k)为k时刻窗口长度为M时互差分序列的方差，M为窗口长度序列。

（4）利用上述（1）、（2）、（3）步骤中得到的各个方差计算观测系统A和观测系统B的量测噪声方差阵。

若当前数据信号量测时刻k小于等于M_Max，量测噪声方差阵为：

$\left\{\begin{array}{c}{R}_1\left(k\right)=\frac{{\mathrm{\σ}}_c\left(k\right)+({\mathrm{\σ}}_1\left(k\right)-{\mathrm{\σ}}_2\left(k\right))}{4}\\ R_2\left(k\right)=\frac{{\mathrm{\σ}}_c\left(k\right)-({\mathrm{\σ}}_1\left(k\right)-{\mathrm{\σ}}_2\left(k\right))}{4}\end{array}\right.$

R1(k)表示数据量测时刻k小于M_Max时观测系统A的量测噪声方差阵；R₂(k)表示数据量测时刻k小于M_Max时观测系统B的量测噪声方差阵。

若当前数据信号量测时刻k大于M_Max，量测噪声方差阵为：

$\left\{\begin{array}{c}{R}_1(M,k)=\frac{{\mathrm{\σ}}_c(M,k)+({\mathrm{\σ}}_1(M,k)-{\mathrm{\σ}}_2(M,k))}{4}\\ R_2(M,k)=\frac{{\mathrm{\σ}}_c(M,k)-({\mathrm{\σ}}_1(M,k)-{\mathrm{\σ}}_2(M,k))}{4}\end{array}\right.M=M_{\mathrm{Min}},M_{\mathrm{Min}}+M_{\mathrm{Dis}}...M_{\mathrm{Max}}-M_{\mathrm{Dis}},M_{\mathrm{Max}}$

其中，k为系统当前数据量测时刻，σ₁(k)、σ₂(k)分别为观测系统A、观测系统B自差分序列在k时刻的方差，σ_C(k)为互差分序列在k时刻的方差，M为窗口长度。R₁(M,k)表示数据量测时刻k大于M_Max时观测系统A的量测噪声方差阵；R₂(M,k)表示数据量测时刻k大于M_Max时观测系统B的量测噪声方差阵。

步骤四：根据两个观测系统在不同窗口长度下计算得到的量测噪声方差阵信息，得到当前时刻最优窗口长度，使用最优窗口长度计算得当前时刻的量测噪声方差阵，具体包括以下几个步骤：

（1）当k小于等于M_Max时，直接按照步骤三（4）中k小于等于M_Max情况下的公式计算量测噪声方差阵。直接使用下式计算量测噪声方差阵。

$\left\{\begin{array}{c}{R}_1\left(k\right)=\frac{{\mathrm{\σ}}_c\left(k\right)+({\mathrm{\σ}}_1\left(k\right)-{\mathrm{\σ}}_2\left(k\right))}{4}\\ R_2\left(k\right)=\frac{{\mathrm{\σ}}_c\left(k\right)-({\mathrm{\σ}}_1\left(k\right)-{\mathrm{\σ}}_2\left(k\right))}{4}\end{array}\right.$

（2）当k大于M_Max时，将k时刻不同窗口下计算得到的噪声方差阵序列使用一次线形函数进行拟合，并记录其斜率。

$\left\{\begin{array}{c}R(M_{\mathrm{Min}},k)=a*(\frac{M_{\mathrm{Min}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1)+b\\ R(M_{\mathrm{Min}}+M_{\mathrm{Dis}},k)=a*(\frac{M_{\mathrm{Min}}+M_{\mathrm{Dis}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1)+b\\ \begin{array}{cc}\begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\end{array}\\ R(M_{\mathrm{Max}}-M_{\mathrm{Dis}},k)=a*(\frac{M_{\mathrm{Max}}-M_{\mathrm{Dis}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1)+b\\ R(M_{\mathrm{Max}},k)=a*(\frac{M_{\mathrm{Max}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1)+b\end{array}\right.$

$\left[\begin{array}{cc}\frac{M_{\mathrm{Min}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1& 1\\ \frac{M_{\mathrm{Min}}+M_{\mathrm{Dis}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1& 1\\ .& .\\ .& .\\ .& .\\ \frac{M_{\mathrm{Max}}-M_{\mathrm{Dis}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1& 1\\ \frac{M_{\mathrm{Max}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1& 1\end{array}\right]{[a,b]}^T=\left[\begin{array}{c}R(M_{\mathrm{Min}},k)\\ R(M_{\mathrm{Min}}+M_{\mathrm{Dis}},k)\\ .\\ .\\ .\\ R(M_{\mathrm{Max}}-M_{\mathrm{Dis}},k)\\ R(M_{\mathrm{Max}},k)\end{array}\right]$

其中，k为当前时刻，a为拟合直线的斜率，b为拟合直线的截距，记：

则拟合直线的斜率和截距即为：

其中R(M_Min，k)、R(M_Min+M_Dis,k)、R(M_Max-M_Dis,k)、R(M_Max,k)分别表示窗口长度为M_Min，(M_Min+M_Dis),(M_Max-M_Dis)，M_Max时的量测噪声方差阵。

（3）将拟合直线计算的斜率带入下式方程得到当前时刻k应使用的数据窗口长度M',并且如果计算得到的窗口长度M'小于M_Min（M_Min取值为100），则取M'为M_Min。

式中：fix(·)表示取整，a_max为所拟合的直线能达到的最大斜率，它反应了量测噪声方差阵的抖动程度，a_max一般取3～7。

（4）利用最优窗口内的序列信息计算得到最终噪声方差阵的估计值。

$\left\{\begin{array}{c}{R}_{1f}(M,k)=\frac{{\mathrm{\σ}}_c(M^{\′},k)+({\mathrm{\σ}}_1(M^{\′},k)-{\mathrm{\σ}}_2(M^{\′},k))}{4}\\ R_{2f}(M,k)=\frac{{\mathrm{\σ}}_c(M^{\′},k)-({\mathrm{\σ}}_1(M^{\′},k)-{\mathrm{\σ}}_2(M^{\′},k))}{4}\end{array}\right.$

R_1f(M,k)、R_2f(M,k)分别表示观测系统A、观测系统B在k时刻估计得到的量测噪声方差阵。通过上述方法，基于不同测量参考系统的自差分序列及互差分序列信息，自适应调节窗口长度，可有效完成对测量噪声方差阵的估计值。

实施例：

应用本发明提出的一种对双观测系统量测噪声方差阵的自适应同步估计方法，分别对观测系统A、观测系统B的测量值加方差值不变白噪声、方差值时变白噪声，使用本发明方法对噪声值进行估计，验证其有效性。

如图2所示，对加了固定方差值白噪声的观测系统A、观测系统B，利用本发明对其量测噪声方差阵进行估计，其中观测系统A所加白噪声方差值为10，观测系统B所加白噪声方差值为5，结果显示，估计值可以有效逼近真值。如图3a和图3b所示，对加了时变方差白噪声的观测系统A和加了固定值白噪声的观测系统B，图3a为使用固定窗口两个系统量测噪声值的估计结果，图3b为使用本发明提出方法的估计结果，其中观测系统A所加白噪声方差为10，在3000s至4000s时间内由10变为50,4000s至5500s时间内为50，5500s至6500s时间内由50变为10,6500s至9000s为10，观测系统B所加白噪声方差为20，本实施例的估计中，使用的最大窗口长度M_Max为1000，最小窗口长度M_Min为100，窗口长度间隔M_Dis为50，根据结果可知，本发明方法的估计结果误差小，抖动小，精度高，能够更加有效逼近真值。

Claims (2)

1.一种对双观测系统量测噪声方差阵的自适应同步估计方法，其特征在于：包括以下几个步骤：

步骤一：分别得到观测系统A、观测系统B在各时刻的数据信号Z1(i)、Z2(i)，其中i为观测系统的数据信号测量时刻；

步骤二：分别计算观测系统A数据信号Z1(i)的自差分序列和观测系统B数据信号Z2(i)的自差分序列，并计算两个自差分序列的互差分序列，具体为：

(1)观测系统A数据信号Z1(i)的自差分序列ΔZ1(i)为：

ΔZ1(i)＝Z1(i)-Z1(i-1)

其中i-1和i分别表示观测系统A的数据采集时刻；

(2)观测系统B数据信号Z2(i)的自差分序列ΔZ2(i)为：

ΔZ2(i)＝Z2(i)-Z2(i-1)

其中i-1和i分别表示观测系统B的数据采集时刻；

(3)观测系统A和观测系统B的两个自差分序列的互差分序列C(i)为：

C(i)＝ΔZ1(i)-ΔZ2(i)

步骤三：分别计算观测系统A数据信号在不同窗口长度下自差分序列的方差，观测系统B数据信号在不同窗口长度下自差分序列的方差，以及两个观测系统在不同窗口长度下互差分序列的方差，并利用方差进行相关计算得到观测系统A、观测系统B不同窗口长度下的量测噪声方差阵R1、R2，具体为：

(1)首先选取最大窗口长度M_Max、最小窗口长度M_Min和窗口长度间隔M_Dis，且M_Max与M_Min之差为窗口长度间隔M_Dis的整数倍，若当前数据信号量测时刻k小于等于M_Max，则直接采用窗口长度累积的方式计算观测系统A自差分序列的方差：

$\left\{\begin{array}{c}{E}_1\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\mathrm{\ΔZ}1\left(i\right)\\ {\mathrm{\σ}}_1\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[\mathrm{\ΔZ}1\left(i\right)-E_1\left(k\right)]}^2\end{array}\right.$

其中E₁(k)表示观测系统A的1到k时刻观测值序列的均值，σ₁(k)表示观测系统A的1到k时刻观测值序列的方差；

若当前数据信号量测时刻k大于M_Max，则计算观测系统A数据信号在不同窗口长度下自差分序列的方差：

$\left\{\begin{array}{c}{E}_1(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}\mathrm{\ΔZ}1\left(i\right)\\ {\mathrm{\σ}}_1(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}{[\mathrm{\ΔZ}1\left(i\right)-E_1(M,k)]}^2\end{array}\right.M=M_{\mathrm{Min}},M_{\mathrm{Min}}+M_{\mathrm{Dis}}...M_{\mathrm{Max}}-M_{\mathrm{Dis}},M_{\mathrm{Max}}$

其中，k为观测系统A当前数据信号量测时刻，i为数据信号量测时刻，ΔZ1(i)为观测系统A自差分序列在i时刻的值，E₁(M,k)为k时刻窗口长度为M时观测系统A自差分序列的均值，σ₁(M,k)为k时刻窗口长度为M时自差分序列的方差，M为窗口长度序列；

(2)若当前数据量测时刻k小于等于M_Max，则直接采用窗口长度累积的方式计算观测系统B自差分序列的方差：

$\left\{\begin{array}{c}{E}_2\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\mathrm{\ΔZ}2\left(i\right)\\ {\mathrm{\σ}}_2\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[\mathrm{\ΔZ}2\left(i\right)-E_2\left(k\right)]}^2\end{array}\right.$

其中E₂(k)表示观测系统B的1到k时刻观测值序列的均值，σ₂(k)表示观测系统B的1到k时刻观测值序列的方差；

若当前数据量测时刻k大于M_Max，则计算观测系统B数据信号在不同窗口长度下的自差分序列的方差：

$\left\{\begin{array}{c}{E}_2(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}\mathrm{\ΔZ}2\left(i\right)\\ {\mathrm{\σ}}_2(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}{[\mathrm{\ΔZ}2\left(i\right)-E_2(M,k)]}^2\end{array}\right.M=M_{\mathrm{Min}},M_{\mathrm{Min}}+M_{\mathrm{Dis}}...M_{\mathrm{Max}}-M_{\mathrm{Dis}},M_{\mathrm{Max}}$

其中，k为观测系统B当前数据信号量测时刻，i为数据信号量测时刻，ΔZ2(i)为观测系统B自差分序列在i时刻的值，E₂(M,k)为k时刻窗口长度为M时观测系统B自差分序列的均值，σ₂(M,k)为k时刻窗口长度为M时的自差分序列方差，M为窗口长度序列；

(3)若当前数据量测时刻k小于等于M_Max，计算两个观测系统互差分序列的方差σ_C(k)：

$\left\{\begin{array}{c}C\left(i\right)=\mathrm{\ΔZ}1\left(i\right)-\mathrm{\ΔZ}2\left(i\right)\\ E_C\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}C\left(i\right)\\ {\mathrm{\σ}}_C\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[C\left(i\right)-E_C\left(k\right)]}^2\end{array}\right.$

其中，C(i)表示观测系统A和观测系统B自差分序列在i时刻的互差分序列；E_C(k)表示k时刻互差分序列的均值；

若当前系统数据量测时刻k大于M_Max，则按下式计算两个量测系统互差分序列的方差：

$\left\{\begin{array}{c}C\left(i\right)=\mathrm{\ΔZ}1\left(i\right)-\mathrm{\ΔZ}2\left(i\right)\\ E_C(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}C\left(i\right)\\ {\mathrm{\σ}}_C(M,k)=\underset{i=k-M}{\overset{k}{\mathrm{\Σ}}}{[C\left(i\right)-E_C(M,k)]}^2\end{array}\right.M=M_{\mathrm{Min}},M_{\mathrm{Min}}+M_{\mathrm{Dis}}...M_{\mathrm{Max}}-M_{\mathrm{Dis}},M_{\mathrm{Max}}$

其中，k为测量系统当前数据信号量测时刻，i为数据信号量测时刻，C(i)表示观测系统A和观测系统B自差分序列在i时刻的互差分序列，E_C(M,k)为k时刻窗口长度为M时互差分序列的均值，σ_C(M,k)为k时刻窗口长度为M时互差分序列的方差，M为窗口长度序列；

(4)利用上述(1)、(2)、(3)步骤中得到的各个方差计算观测系统A和观测系统B的量测噪声方差阵；

若当前数据信号量测时刻k小于等于M_Max，量测噪声方差阵为：

$\left\{\begin{array}{c}{R}_1\left(k\right)=\frac{{\mathrm{\σ}}_c\left(k\right)+({\mathrm{\σ}}_1\left(k\right)-{\mathrm{\σ}}_2\left(k\right))}{4}\\ R_2\left(k\right)=\frac{{\mathrm{\σ}}_c\left(k\right)-({\mathrm{\σ}}_1\left(k\right)-{\mathrm{\σ}}_2\left(k\right))}{4}\end{array}\right.$

R₁(k)表示数据量测时刻k小于等于M_Max时观测系统A的量测噪声方差阵；R₂(k)表示数据量测时刻k小于等于M_Max时观测系统B的量测噪声方差阵；

若当前数据信号量测时刻k大于M_Max，量测噪声方差阵为：

$\left\{\begin{array}{c}{R}_1(M,k)=\frac{{\mathrm{\σ}}_c(M,k)+({\mathrm{\σ}}_1(M,k)-{\mathrm{\σ}}_2(M,k))}{4}\\ R_2(M,k)=\frac{{\mathrm{\σ}}_c(M,k)-({\mathrm{\σ}}_1(M,k)-{\mathrm{\σ}}_2(M,k))}{4}\end{array}\right.M=M_{\mathrm{Min}},M_{\mathrm{Min}}+M_{\mathrm{Dis}}...M_{\mathrm{Max}}-M_{\mathrm{Dis}},M_{\mathrm{Max}}$

其中，k为系统当前数据量测时刻，σ₁(k)、σ₂(k)分别为观测系统A、观测系统B自差分序列在k时刻的方差，σ_C(k)为互差分序列在k时刻的方差，M为窗口长度，R₁(M,k)表示数据量测时刻k大于M_Max时观测系统A的量测噪声方差阵；R₂(M,k)表示数据量测时刻k大于M_Max时观测系统B的量测噪声方差阵；

步骤四：根据两个观测系统在不同窗口长度下计算得到的量测噪声方差阵信息，得到当前时刻最优窗口长度，使用最优窗口长度计算得当前时刻的量测噪声方差阵，具体为：

(1)当k小于等于M_Max时，直接按照步骤三(4)中k小于等于M_Max情况下的公式计算量测噪声方差阵：

$\left\{\begin{array}{c}{R}_1\left(k\right)=\frac{{\mathrm{\σ}}_c\left(k\right)+({\mathrm{\σ}}_1\left(k\right)-{\mathrm{\σ}}_2\left(k\right))}{4}\\ R_2\left(k\right)=\frac{{\mathrm{\σ}}_c\left(k\right)-({\mathrm{\σ}}_1\left(k\right)-{\mathrm{\σ}}_2\left(k\right))}{4}\end{array}\right.$

(2)当k大于M_Max时，将k时刻不同窗口下计算得到的噪声方差阵序列使用一次线形函数进行拟合，并记录其斜率；

$\left\{\begin{array}{c}R(M_{\mathrm{Min}},k)=a*(\frac{M_{\mathrm{Min}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1)+b\\ R(M_{\mathrm{Min}}+M_{\mathrm{Dis}},k)=a*(\frac{M_{\mathrm{Min}}+M_{\mathrm{Dis}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1)+b\\ \begin{array}{cc}.& .\end{array}\\ \begin{array}{cc}.& .\end{array}\\ \begin{array}{cc}.& .\end{array}\\ R(M_{\mathrm{Max}}-M_{\mathrm{Dis}},k)=a*(\frac{M_{\mathrm{Max}}-M_{\mathrm{Dis}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1)+b\\ R(M_{\mathrm{Max}},k)=a*(\frac{M_{\mathrm{Max}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1)+b\end{array}\right.$

$\left[\begin{array}{cc}\frac{M_{\mathrm{Min}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1& 1\\ \frac{M_{\mathrm{Min}}+M_{\mathrm{Dis}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1& 1\\ .& .\\ .& .\\ .& .\\ \frac{M_{\mathrm{Max}}-M_{\mathrm{Dis}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1& 1\\ \frac{M_{\mathrm{Max}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1& 1\end{array}\right]{[a,b]}^T=\left[\begin{array}{c}R(M_{\mathrm{Min}},k)\\ R(M_{\mathrm{Min}}+M_{\mathrm{Dis}},k)\\ .\\ .\\ .\\ R(M_{\mathrm{Max}}-M_{\mathrm{Dis}},k)\\ R(M_{\mathrm{Max}},k)\end{array}\right]$

其中，k为当前时刻，a为拟合直线的斜率，b为拟合直线的截距：

$K=\left[\begin{array}{cc}\frac{M_{\mathrm{Min}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1& 1\\ \frac{M_{\mathrm{Min}}+M_{\mathrm{Dis}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1& 1\\ .& .\\ .& .\\ .& .\\ \frac{M_{\mathrm{Max}}-M_{\mathrm{Dis}}-M_{\mathrm{Min}}}{M_{\mathrm{Dis}}}+1& 1\\ \frac{M_{\mathrm{Max}}-M_{\min }}{M_{\mathrm{Dis}}}+1& 1\end{array}\right]$

则拟合直线的斜率和截距为：

其中R(M_Min,k)、R(M_Min+M_Dis,k)、R(M_Max-M_Dis,k)、R(M_Max,k)分别表示窗口长度为M_Min,(M_Min+M_Dis),(M_Max-M_Dis),M_Max时的量测噪声方差阵；

(3)将拟合直线计算的斜率带入下式方程得到当前时刻k应使用的数据窗口长度M'，若计算得到的窗口长度M'小于M_Min，则取M'为M_Min，且M_Min取值为100：

其中：fix(·)表示取整，a_max为所拟合的直线达到的最大斜率；

(4)利用最优窗口内的序列信息计算得到最终噪声方差阵的估计值：

$\left\{\begin{array}{c}{R}_{1f}(M,k)=\frac{{\mathrm{\σ}}_c(M^{\′},k)+({\mathrm{\σ}}_1(M^{\′},k)-{\mathrm{\σ}}_2(M^{\′},k))}{4}\\ R_{2f}(M,k)=\frac{{\mathrm{\σ}}_c(M^{\′},k)-({\mathrm{\σ}}_1(M^{\′},k)-{\mathrm{\σ}}_2(M^{\′},k))}{4}\end{array}\right.$

R_1f(M,k)、R_2f(M,k)分别表示观测系统A、观测系统B在k时刻估计得到的量测噪声方差阵。

2.根据权利要求1所述的一种对双观测系统量测噪声方差阵的自适应同步估计方法，其特征在于：所述的最大斜率a_max取值为3～7。