[发明专利]一种基于α‑β‑γ滤波和二阶互差分的噪声方差测量方法

摘要

本发明公开了一种基于α‑β‑γ滤波和二阶互差分的噪声方差测量方法，该方法通过α‑β‑γ滤波构造对测量信号的伪冗余测量，使用数据选择算法选取测量信号中的缓变或者直线变化部分，最终利用二阶互差分的方法计算出测量噪声的方差。本发明能够对任意信号的噪声方差实时估计，在Kalman滤波中，能够有效提高滤波精度。

主权项

一种基于α‑β‑γ滤波和二阶互差分的噪声方差测量方法，包括以下几个步骤：步骤一：对α‑β‑γ滤波器进行离线设计，利用包络的方法求得原始测量信号噪声幅度和滤波器输出噪声幅度，通过调整滤波器的参数h，使平滑型α‑β‑γ滤波器输出的噪声幅度为原信号噪声幅度的1/20～1/10，取跟踪型α‑β‑γ滤波器的h值为平滑型的10～100倍；具体包括：A.建立α‑β‑γ滤波器模型X·(t)=A(t)X(t)+F(t)·W(t)]]>Z(t)＝H(t)X(t)+V(t)其中：为反映信号变化的状态量，x1(t)、x2(t)、x3(t)为状态量中元素；W(t)为系统噪声，Z(t)为观测值，V(t)为观测噪声，为状态矩阵，为噪声输入矩阵，H(t)＝[100]为观测矩阵；在Kalman滤波公式中：X^·(t|t)=A(t)X^(t|t)+K(t)(Z(t)-H(t)X^(t|t))]]>K(t)＝P(t|t)HT(t)R‑1(t)P·(t|t)=A(t)P(t|t)+P(t|t)AT(t)-P(t|t)HT(t)R-1(t)H(t)P(t|t)+F(t)Q(t)F-1(t)]]>其中：为t时刻的状态估计值，为t时刻的状态估计值导数，R(t)为观测噪声协方差阵，Q(t)为系统噪声协方差阵，P(t|t)为滤波噪声方差阵，为滤波噪声方差阵的导数，K(t)为增益矩阵；在滤波达到稳定时，K(t)保持定值K：K=2h32h23h]]>h=qr]]>其中：q为系统噪声方差，r为测量噪声方差，h为系统噪声方差与状态噪声方差之比，是增益矩阵K参数；B.任取h值，使用α‑β‑γ滤波器对信号进行跟踪，构建离散化的α‑β‑γ滤波模型：X(k+1)＝Φ·X(k)+WZ(k+1)＝H·X(k+1)+V其中：X(k+1)、Z(k+1)分别为k+1时刻的状态值和观测值，Φ为状态转移矩阵，W为系统噪声矩阵，H为观测矩阵，V为观测噪声矩阵，值分别为：Φ=1T1T1]]>W=00w]]>H＝[1 0 0]V＝v其中，T为采样间隔，w为噪声方差为q的白噪声，v为噪声方差为r的白噪声；计算α‑β‑γ滤波：X^(k)=Φ(k,k-1)X^(k-1)+K(Z(k)-HΦ(k,k-1)X^(k-1))]]>其中：为k时刻的状态估计值；C.对平滑型α‑β‑γ滤波器输出值进行包络识别，得到滤波器的噪声幅度|V1|：|V1|=Σ(upperEnv1k-lowerEnv1k)2m1]]>其中：m1为整个原始测量信号的数据数量，upperEnv1k和lowerEnv1k为滤波器在k时刻点的上、下包络点；D.对整个原始测量信号进行包络分析，得到原信号的噪声幅度|V2|：|V2|=Σ(upperEnv2k-lowerEnv2k)2m1]]>其中：upperEnv2k和lowerEnv2k为原始测量信号在k时刻点的上、下包络点；E.调整滤波器h值，重复B与C两步，直到滤波器噪声幅度|V1|与原始测量信号噪声幅度|V2|满足：130≤|V1||V2|≤110]]>F.取此时的h值为平滑型α‑β‑γ滤波器的h值hS，取hS的10～100倍为跟踪型的α‑β‑γ滤波器h值hT；步骤二：构造出平滑型和跟踪型α‑β‑γ滤波器后，对传感器采集的实时信号的噪声方差进行测量，传感器实时采集的信号为原始测量信号，使用平滑型α‑β‑γ滤波器和跟踪型α‑β‑γ滤波器分别在线对其进行滤波，得到平滑型滤波序列与跟踪型滤波序列以m为窗口长度，在序列和中分别滑动获得数据段，以k时刻开始的序列段和为：X^S|kk+m=x^S(k)x^S(k+1)...x^S(k+m)T]]>X^T|kk+m=x^T(k)x^T(k+1)...x^T(k+m)T]]>其中：和分别为k时刻的平滑型滤波器和跟踪型滤波器输出；步骤三：使用数据选择算法对跟踪型、平滑型α‑β‑γ滤波器窗口内输出进行筛选，选出缓变信号或者线性变化信号；数据选择算法筛选出缓变或者线性变化数据的具体过程为：A.建立线性最小二乘模型：X^S|kk+m=HlaS(k)bS(k)]]>X^T|kk+m=HlaT(k)bT(k)]]>Hl=1121......m1]]>其中：Hl为最小二乘观测矩阵，aS(k)、bS(k)与aT(k)、bT(k)为拟合直线x(k+i)＝a·i+b的a、b的估值；B.使用最小二乘估计aS(k)、bS(k)与aT(k)、bT(k)：aS(k)bS(k)=(HlTHl)-1HlTX^S|kk+m]]>aT(k)bT(k)=(HlTHl)-1HlTX^T|kk+m]]>C.将窗口内数据按时间分为前后两段，以拟合直线为中心，分别统计每段内的上下平均幅值A1、A2、A3、A4：A1=mean{x^S(k+i)-aS(k)i-bS|x^S(k+i)≥aS(k)i-bS,i∈(0,m/2)}]]>A2=mean{x^S(k+i)-aS(k)i-bS|x^S(k+i)<aS(k)i-bS,i∈(0,m/2)}]]>A3=mean{x^S(k+i)-aS(k)i-bS|x^S(k+i)≥aS(k)i-bS,i∈[m/2,m)}]]>A4=mean{x^S(k+i)-aS(k)i-bS|x^S(k+i)<aS(k)i-bS,i∈[m/2,m)}]]>其中：mean代表求平均值的函数；D.计算两拟合直线之间夹角θ的绝对值θ：|θ|＝|arctan(aS(k))‑arctan(aT(k))|E.判断A1、A2、A3、A4是否一致和两拟合直线间夹角是否足够小：|A1A2|,|A1A3|,|A3A4|∈(1-th_amp,1+th_amp)|θ|<th_ang]]>其中：th_amp和th_ang分别为幅度和角度的阈值；当A1、A2、A3、A4绝对值之比在范围(1‑th_amp,1+th_amp)以内且拟合直线夹角的绝对值θ小于阈值th_ang时，则该数据段是缓变或线性的原始测量信号数据段，否则，数据段为快变过程，舍弃；步骤四：将步骤三得到的缓变或线性的原始测量信号数据段，与相应时刻的平滑型α‑β‑γ滤波器输出进行二阶互差分运算，得到噪声方差；具体的：A.设以时刻k开始，长度为m的缓变或线性原始测量信号数据段以及其平滑型α‑β‑γ滤波器输出将其分别自差分得：Δx(k)＝x(k+1)‑x(k)Δx^S(k)=x^S(k+1)-x^S(k)]]>其中，x(k)为信号k时刻测量值，为平滑型滤波器在k时刻的输出值，Δx(k)和分别为它们的自差分值，并构成自差分序列和B.将得到的自差分序列和进行互差分，求得窗口内原信号数据段的噪声方差Var(X|kk+m)=12(ΔX|kk+m-ΔX^S|kk+m)(ΔX|kk+m-ΔX^S|kk+m)T]]>C.计算原信号的噪声方差，并对计算出的数据进行加权得到：di=(1-b)/(1-bi+1)Var(X)i=(1-di)Var(X|kk+m)+diVar(X)i-1,0<b<1]]>其中：X为原始测量信号序列，b为遗忘因子，di为加权系数，bi+1是遗忘因子b的i+1次幂，Var(X)i为第i次计算的原信号噪声方差。