[0070] 其中：?为最小二乘观测矩阵，a s(k)、bs(k)与aT(k)、bT(k)为拟合直线x(k+i)= a · i+b的a、b的估值。
[0071] B.使用最小二乘估计 as(k)、bs(k)与 aT(k)、bT(k):
[0074] C.将窗口内数据按时间分为前后两段，以拟合直线为中心，分别统计每段内 ?的上下平均幅值A1、A2、A3、A4:
[0079] 其中：mean代表求平均值的函数。
[0080] D.计算两拟合直线之间夹角Θ的绝对值I Θ I :
[0081] Θ I = I arctan (as (k)) -arctan (aT (k))
[0082] E.判断Ap A2、A3、A4是否一致和两拟合直线间夹角是否足够小：
[0084] 其中：th_amp和th_ang分别为幅度和角度的阈值。
[0085] 当Ap A2、A3、A4绝对值之比在范围（l-th_amp, l+th_amp)以内且拟合直线夹角的 绝对值I θ I小于阈值th_ang时，认为'"六4一致且拟合直线间夹角较小，则该数据 段是缓变或线性变化过程的数据段，可以用于计算噪声方差。反之，数据段为快变过程，不 能用于计算噪声方差。
[0086] 步骤四：步骤三筛选出的缓变或者线性变化部分原测量信号，与相应时刻的平滑 型α-β-γ滤波器输出进行二阶互差分运算，得到噪声方差。
[0087] Α.设以时刻k开始，长度为m的缓变或线性原测量信号数据段Χ|Γ以及其平滑 型《 - β - 丫滤波器输出尤 ，将其分别自差分得：
[0090] 其中，x(k)为信号k时刻测量值，?\(/〇为平滑型滤波器在k时刻的输出值， Δ X (k)和⑷分别为它们的自差分值，并构成自差分序列从和Als 。
[0091] B.将得到的自差分序列ΔΧ|Γ和I广进行互差分，求得窗口内原始测量信号 数据段Xpm的噪声方差％:
[0093] C.计算原始测量信号的噪声方差，并对计算出的数据进行加权得到：
[0095] 其中：Χ为原始测量信号序列，b为遗忘因子，Cl1为加权系数，b 1+1是遗忘因子b的 i+Ι次幂，Var (X)1为第i次计算的原始测量信号噪声方差。
[0096] 通过上述方法，构造原始测量信号的虚拟冗余测量序列，选择出其中的缓变或者 线性变化过程，通过消去真实值与相对于真实值的相对变化量，得到反映原始测量信号噪 声的序列，最终计算出原始测量信号噪声的方差。该方差可以用于自适应Kalman滤波算 法，提尚算法精度。
[0097] 实例中：
[0098] 取图2所示信号进行仿真，给信号中加入标准差为5的零均值高斯白噪声作为测 量信号，如图3所示。
[0099] 依据上述步骤，本实例仿真如图4所示，虚线框内数据为选择出的缓变信号。使用 蒙特卡洛法仿真得到的相对误差分布如图5所示。从图5中可以看出，本发明的方法计算 噪声方差误差在20 %以内，90. 2 %的测量精度在10 %以内。可以看出，本发明的方法具有 较高的计算精度。
[0100] 本发明在二阶互差分的基础上，使用α - β -γ滤波器构造冗余测量，通过对数据 和α-β-γ滤波器进行筛选，以达到满足二阶互差分条件的目的，最终使用二阶互差分的 方法计算出测量噪声的方差。
【主权项】
1. 一种基于Ct - β -γ滤波和二阶互差分的噪声方差测量方法，包括以下几个步骤： 步骤一：对α-β-γ滤波器进行离线设计，利用包络的方法求得原始测量信号噪声幅 度和滤波器输出噪声幅度，通过调整滤波器的参数h，使平滑型α-β-γ滤波器输出的噪 声幅度为原信号噪声幅度的1/10~1/20,取跟踪型α-β-γ滤波器的h值为平滑型的 10~100倍； 具体包括： A.建立α - β - γ滤波器模型号变化的状态量，X1U)、X2(t)、x 3(t)为状态量中元 素；W(t)为系统噪声，Z(t)为观测值，V(t)为观测噪声，：为状态矩阵，为噪声输入矩阵，H(t) = [I 0 0]为观测矩阵； 在Kalman滤波公式中：其中：为t时刻的状态估计值，为t时刻的状态估计值导数，R(t)为观测 噪声协方差阵，Q(t)为系统噪声协方差阵，P(t 11)为滤波噪声方差阵，为滤波噪声方 差阵的导数，K (t)为增益矩阵； 在滤波达到稳定时，K (t)保持定值K :其中：q为系统噪声方差，r为测量噪声方差，h为系统噪声方差与状态噪声方差之比， 是增益矩阵K参数； B. 任取h值，使用α-β-γ滤波器对信号进行跟踪，构建离散化的α-β-γ滤波模 型：其中：X(k+l)、Z(k+l)分别为k+Ι时刻的状态值和观测值，Φ为状态转移矩阵，W为系 统噪声矩阵，H为观测矩阵，V为观测噪声矩阵，值分别为：其中，T为采样间隔，w为噪声方差为q的白噪声，V为噪声方差为r的白噪声； 计算α-β-γ滤波：其中：为k时刻的状态估计值； C. 对平滑型α-β-γ滤波器输出值进行包络识别，得到滤波器的噪声幅度IV1I :其中：In1为整个原始测量信号的数据数量，upperEnv ^和IowerEnv lk为滤波器在k时刻 点的上、下包络点； D. 对整个原始测量信号进行包络分析，得到原信号的噪声幅度|V2| :其中：UpperEnva^PI IowerEnva为原始测量信号在k时刻点的上、下包络点； E. 调整滤波器h值，重复B与C两步，直到滤波器噪声幅度IV1I与原始测量信号噪声 幅度IV21满足：F. 取此时的h值为平滑型α - β - γ滤波器的h值hs，取匕的10~100倍为跟踪型的 α - β - γ滤波器h值hT; 步骤二：构造出平滑型和跟踪型α-β-γ滤波器后，对传感器采集的实时信号的噪声 方差进行测量，传感器实时采集的信号为原始测量信号，使用平滑型α-β-γ滤波器和跟 踪型α-β-γ滤波器分别在线对其进行滤波，得到平滑型滤波序列尤与跟踪型滤波序列 以m为窗口长度，在序列中分别滑动获得数据段，以k时刻开始的序列段其中：心斤）和奂分别为k时刻的平滑型滤波器和跟踪型滤波器输出； 步骤三：使用数据选择算法对跟踪型、平滑型α-β-γ滤波器窗口内输出进行筛选， 选出缓变信号或者线性变化信号； 数据选择算法筛选出缓变或者线性变化数据的具体过程为： Α.建立线性最小二乘模型：其中：托为最小二乘观测矩阵，为拟合直线x(k+i) = a *i+b 的a、b的估值； B. 使用最小二乘估计C. 将窗口内数据按时间分为前后两段，以拟合直线为中心，分别统计每段内的 上下平均幅值ApApApA4:其中：mean代表求平均值的函数； D. 计算两拟合直线之间夹角Θ的绝对值I θ I : Θ I = I arctan (as (k)) -arctan (aT (k)) E. 判断A1、A2、A3、A4是否一致和两拟合直线间夹角是否足够小 ：其中：th_amp和th_ang分别为幅度和角度的阈值； 当Ap A2、A3、A4绝对值之比在范围（l-th_amp, l+th_amp)以内且拟合直线夹角的绝对 值I θ I小于阈值th_ang时，则该数据段是缓变或线性的原始测量测量信号数据段，否则， 数据段为快变过程，舍弃； 步骤四：将步骤三得到的缓变或线性的原始测量测量信号数据段，与相应时刻的平滑 型α-β-γ滤波器输出进行二阶互差分运算，得到噪声方差； 具体的： Α.设以时刻k开始，长度为m的缓变或线性原始测量测量信号数据段以及其平滑 型α-β-γ滤波器输出，将其分别自差分得：其中，x(k)为信号k时刻测量值，为平滑型滤波器在k时刻的输出值，AX(k)和分别为它们的自差分值，并构成自差分序列2B. 将得到的自差分序列进行互差分，求得窗口内原信号数据段C. 计算原信号的噪声方差，并对计算出的数据进行加权得到：其中：X为原始测量测量信号序列，b为遗忘因子，Cl1为加权系数，b 1+1是遗忘因子b的 i+Ι次幂，Var (X)1为第i次计算的原信号噪声方差。
【专利摘要】本发明公开了一种基于α-β-γ滤波和二阶互差分的噪声方差测量方法，该方法通过α-β-γ滤波构造对测量信号的伪冗余测量，使用数据选择算法选取测量信号中的缓变或者直线变化部分，最终利用二阶互差分的方法计算出测量噪声的方差。本发明能够对任意信号的噪声方差实时估计，在Kalman滤波中，能够有效提高滤波精度。
【IPC分类】G01D3/032
【公开号】CN105180971
【申请号】CN201510582072
【发明人】张海, 张义昕, 张晓鸥, 郭雷
【申请人】北京航空航天大学
【公开日】2015年12月23日
【申请日】2015年9月14日