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

摘要

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

主权项

一种对双观测系统量测噪声方差阵的自适应同步估计方法，其特征在于：包括以下几个步骤：步骤一：分别得到观测系统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}{\Sigma }}\mathrm{\Delta Z}1\left(i\right)\\ {\sigma }_1\left(k\right)=\underset{i=1}{\overset{k}{\Sigma }}{[\mathrm{\Delta 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}{\Sigma }}\mathrm{\Delta Z}1\left(i\right)\\ {\sigma }_1(M,k)=\underset{i=k-M}{\overset{k}{\Sigma }}{[\mathrm{\Delta 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}{\Sigma }}\mathrm{\Delta Z}2\left(i\right)\\ {\sigma }_2\left(k\right)=\underset{i=1}{\overset{k}{\Sigma }}{[\mathrm{\Delta 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}{\Sigma }}\mathrm{\Delta Z}2\left(i\right)\\ {\sigma }_2(M,k)=\underset{i=k-M}{\overset{k}{\Sigma }}{[\mathrm{\Delta 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{\Delta Z}1\left(i\right)-\mathrm{\Delta Z}2\left(i\right)\\ E_C\left(k\right)=\underset{i=1}{\overset{k}{\Sigma }}C\left(i\right)\\ {\sigma }_C\left(k\right)=\underset{i=1}{\overset{k}{\Sigma }}{[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{\Delta Z}1\left(i\right)-\mathrm{\Delta Z}2\left(i\right)\\ E_C(M,k)=\underset{i=k-M}{\overset{k}{\Sigma }}C\left(i\right)\\ {\sigma }_C(M,k)=\underset{i=k-M}{\overset{k}{\Sigma }}{[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{{\sigma }_c\left(k\right)+({\sigma }_1\left(k\right)-{\sigma }_2\left(k\right))}{4}\\ R_2\left(k\right)=\frac{{\sigma }_c\left(k\right)-({\sigma }_1\left(k\right)-{\sigma }_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{{\sigma }_c(M,k)+({\sigma }_1(M,k)-{\sigma }_2(M,k))}{4}\\ R_2(M,k)=\frac{{\sigma }_c(M,k)-({\sigma }_1(M,k)-{\sigma }_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{{\sigma }_c\left(k\right)+({\sigma }_1\left(k\right)-{\sigma }_2\left(k\right))}{4}\\ R_2\left(k\right)=\frac{{\sigma }_c\left(k\right)-({\sigma }_1\left(k\right)-{\sigma }_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{{\sigma }_c(M^{\prime },k)+({\sigma }_1(M^{\prime },k)-{\sigma }_2(M^{\prime },k))}{4}\\ R_{2f}(M,k)=\frac{{\sigma }_c(M^{\prime },k)-({\sigma }_1(M^{\prime },k)-{\sigma }_2(M^{\prime },k))}{4}\end{array}\right.$R_1f(M,k)、R_2f(M,k)分别表示观测系统A、观测系统B在k时刻估计得到的量测噪声方差阵。