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

权利要求书

1.一种基于α-β-γ滤波和二阶互差分的噪声方差测量方法，包括以下几个步骤：

步骤一：对α-β-γ滤波器进行离线设计，利用包络的方法求得原始测量信号噪声幅度和滤波器输出噪声幅度，通过调整滤波器的参数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)

其中：为反映信号变化的状态量，x₁(t)、x₂(t)、x₃(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)H^T(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)+W

Z(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.对平滑型α-β-γ滤波器输出值进行包络识别，得到滤波器的噪声幅度|V₁|：

|V1|=Σ(upperEnv1k-lowerEnv1k)2m1]]>

其中：m₁为整个原始测量信号的数据数量，upperEnv_1k和lowerEnv_1k为滤波器在k时刻点的上、下包络点；

D.对整个原始测量信号进行包络分析，得到原信号的噪声幅度|V₂|：

|V2|=Σ(upperEnv2k-lowerEnv2k)2m1]]>

其中：upperEnv_2k和lowerEnv_2k为原始测量信号在k时刻点的上、下包络点；

E.调整滤波器h值，重复B与C两步，直到滤波器噪声幅度|V₁|与原始测量信号噪声幅度|V₂|满足：

130≤|V1||V2|≤110]]>

F.取此时的h值为平滑型α-β-γ滤波器的h值h_S，取h_S的10～100倍为跟踪型的α-β-γ滤波器h值h_T；

步骤二：构造出平滑型和跟踪型α-β-γ滤波器后，对传感器采集的实时信号的噪声方差进行测量，传感器实时采集的信号为原始测量信号，使用平滑型α-β-γ滤波器和跟踪型α-β-γ滤波器分别在线对其进行滤波，得到平滑型滤波序列与跟踪型滤波序列以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]]> 2

其中：和分别为k时刻的平滑型滤波器和跟踪型滤波器输出；

步骤三：使用数据选择算法对跟踪型、平滑型α-β-γ滤波器窗口内输出进行筛选，选出缓变信号或者线性变化信号；

数据选择算法筛选出缓变或者线性变化数据的具体过程为：

A.建立线性最小二乘模型：

X^S|kk+m=HlaS(k)bS(k)]]>

X^T|kk+m=HlaT(k)bT(k)]]>

Hl=1121......m1]]>

其中：H_l为最小二乘观测矩阵，a_S(k)、b_S(k)与a_T(k)、b_T(k)为拟合直线x(k+i)＝a·i+b的a、b的估值；

B.使用最小二乘估计a_S(k)、b_S(k)与a_T(k)、b_T(k)：

aS(k)bS(k)=(HlTHl)-1HlTX^S|kk+m]]>

aT(k)bT(k)=(HlTHl)-1HlTX^T|kk+m]]>

C.将窗口内数据按时间分为前后两段，以拟合直线为中心，分别统计每段内的上下平均幅值A₁、A₂、A₃、A₄：

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(a_S(k))-arctan(a_T(k))|

E.判断A₁、A₂、A₃、A₄是否一致和两拟合直线间夹角是否足够小：

|A1A2|,|A1A3|,|A3A4|∈(1-th_amp,1+th_amp)|θ|<th_ang]]>

其中：th_amp和th_ang分别为幅度和角度的阈值；

当A₁、A₂、A₃、A₄绝对值之比在范围(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为遗忘因子，d_i为加权系数，bⁱ⁺¹是遗忘因子b的i+1次幂，Var(X)_i为第i次计算的原信号噪声方差。