基于时域与频域的卡尔曼平滑滤波方法及平滑滤波器技术

本发明专利技术涉及信号平滑滤波技术领域，为基于时域与频域的卡尔曼平滑滤波方法及平滑滤波器，包括获取待滤波信号；采用泰勒级数低阶近似模型，根据信号采样频率和模型维数，确定平滑滤波器的状态方程和量测方程；根据卡尔曼平滑滤波算法，确定平滑滤波器的频率响应函数，进而可得到其幅频特性函数；根据所设计平滑滤波器的幅频特性参数，估算系统噪声和量测噪声的统计特性参数；根据卡尔曼滤波模型及其噪声统计特性参数，采用卡尔曼平滑滤波算法，对所述待滤波信号进行平滑滤波处理。由于采用双向卡尔曼滤波，估值精度优于传统的卡尔曼低通滤波器。在频率域，可同传统的FIR低通滤波器一样，使用方便，且对于低信噪比测量信号，可以获得更高的滤波精度。

【技术实现步骤摘要】
基于时域与频域的卡尔曼平滑滤波方法及平滑滤波器
本专利技术涉及信号平滑滤波
，具体涉及基于时域与频域的卡尔曼平滑滤波方法及平滑滤波器。
技术介绍
卡尔曼滤波技术广泛应用于信号测量、信号传输、智能控制、GPS动态数据处理和惯性导航等领域。由于外界环境、动态变化等因素的影响，经传感器获得的测量信号、以及经数据链路传输的测量信号通常包含大量噪声。噪声通常以高频为主，需要采用平滑滤波技术处理，消除高频噪声干扰，以提取信号中低频中的有效成分信息。特别是，当测量信号的信噪比较低时，采用平滑滤波技术对信号进行滤波处理，十分必要。卡尔曼滤波应用需要建立准确的数学模型，且要求噪声的统计特性是已知的。而实际情况并非如此，模型的偏差及噪声的不确定性,都会对滤波结果和滤波精度产生影响(甚至发散)。
技术实现思路
本专利技术提供了基于时域与频域的卡尔曼平滑滤波方法及平滑滤波器，解决了噪声统计特性参数未知的问题，提出了基于频域的噪声统计特性参数估算方法，使得所设计的卡尔曼平滑滤波器，不仅在时域，也可在频域，对测量信号进行滤波处理。本专利技术提供了基于时域与频域的卡尔曼平滑滤波方法，包括：S1，获取待滤波信号；S2，采用泰勒级数低阶近似模型，根据信号采样频率和模型维数，确定平滑滤波器的状态方程和量测方程；S3，根据卡尔曼平滑滤波算法，确定平滑滤波器的频率响应函数，进而可得到其幅频特性函数；根据所设计平滑滤波器的幅频特性参数，估算系统噪声和量测噪声的统计特性参数；S4，根据卡尔曼滤波模型及其噪声特性参数，采用卡尔曼平滑滤波算法，对所述待滤波信号进行平滑滤波处理。可选的，采用泰勒级数低阶近似模型，滤波器的状态方程为：Xk＝ΦXk-1+ΓWk-1(1)式中，Wk＝[Wk]Xk为tk时刻的状态向量,Φ为tk-1时刻至tk时刻的一步转移矩阵,Γ为系统噪声耦合矩阵,Wk为零均值白噪声或高斯白噪声，其协方差矩阵为Qk；其中，Xk中包含xk和xk的m阶导数(m＝1,2,…,n-1)的状态信息序列，例：表示xk的一阶导数，ts为采样间隔，n为模型维数。可选的，采用泰勒级数低阶近似模型，滤波器的量测方程为Zk＝HXk+Vk(2)式中，Zk＝[x′k]，H＝[100...0]，Vk＝[Vk]，Zk为tk时刻的量测向量，x′k为待滤波信号，H为量测矩阵，Vk为零均值白噪声或高斯白噪声，其协方差矩阵为Rk。可选的，所述卡尔曼平滑滤波算法包括正向滤波算法，正向滤波算法如下：状态一步预测状态估计滤波增益Kk＝Pk/k-1HT(HPk/k-1HT+Rk)-1(5a)或一步预测均方误差Pk/k-1＝ΦPk-1ΦT+ΓQk-1ΓT(6)估计均方误差或Pk＝(I-KkH)Pk/k-1(7b)或式中,I为单位矩阵。给定初始状态向量和初始估计均方误差矩阵P0，按(3)～(7)式，(7)式包括(7a)、(7b)及(7c)，其他的标号也一样，由量测向量{Zk}(k＝1,2,…,M)，递推得到状态向量估计表示滤波器Xk的先验估计，表示滤波器Xk-1的后验估计，Zk为待滤波信号，Kk表示卡尔曼滤波增益，Φ为tk-1时刻至tk时刻的一步转移矩阵。可选的，所述卡尔曼平滑滤波算法还包括反向滤波算法，反向滤波算法如下：平滑方程:平滑均方误差:令在正向滤波基础上，按(8)式和(9)式，由状态向量估计递推得到平滑后的状态向量估计可选的，所述S3具体包括：输入信号为Zk(x′k)，输出信号为即根据卡尔曼平滑滤波器结构和平滑滤波算法，导出的平滑滤波器频率响应函数为：其中，I为单位矩阵，j为虚数单位，ω为圆频率，Φ为tk-1时刻至tk时刻的一步转移矩阵。可选的，在稳态时，根据滤波算法导出P＝(I-PHTR-1H)(ΦPΦT+ΓQΓT)；其中，P表示滤波状态估计误差的协方差矩阵。可选的，按照平滑滤波算法，在稳态条件下，频率响应函数表示为：将包含P矩阵代数方程的解代入频率响应函数H(ejω)，并以此为基础，导出滤波器的幅频特性函数|H(f)|(ω＝2πf),在噪声参数确定的条件下，可求出滤波器对应的截止频率fc。可选的，在量测噪声方差R估算基础上，根据滤波模型、滤波器的频率响应函数H(ejω)和滤波器截止频率fc，估算系统噪声方差Q。本专利技术还提供了用于基于时域与频域的卡尔曼平滑滤波方法的平滑滤波器，包括：待滤波数据获取模块，用于获取待滤波信号及数据采样间隔ts参数；卡尔曼滤波模型和滤波参数确定模块，用于确定系统状态方程、量测方程和滤波器的截止频率；预测方程和更新方程确定模块，用于根据卡尔曼平滑滤波算法,确定出平滑滤波器的预测方程、滤波更新方程和平滑更新方程；噪声统计特性参数确定模块，用于根据卡尔曼平滑滤波算法，确定平滑滤波器的频率响应函数，进而得到其幅频特性函数，并以此为基础，由滤波器截止频率，确定量测噪声序列方差R和系统噪声序列方差Q。有益效果：本专利技术提供了基于时域与频域的卡尔曼平滑滤波方法及平滑滤波器，包括：获取待滤波信号；采用泰勒级数低阶近似模型，根据信号采样频率和模型维数，确定平滑滤波器的状态方程和量测方程；根据卡尔曼平滑滤波算法，确定平滑滤波器的频率响应函数，进而可得到其幅频特性函数；根据所设计平滑滤波器的幅频特性参数，估算系统噪声和量测噪声的统计特性参数；根据滤波模型及其噪声特性参数，采用卡尔曼平滑滤波算法，对所述待滤波信号进行平滑滤波处理。通过采用固定区间平滑滤波算法，利用在整个区间内采样得到的所有测量值，得到每次采样时的最优估计。由于采用双向卡尔曼滤波，估值精度优于传统的卡尔曼低通滤波器。且在频率域，同传统的FIR低通滤波器一样，使用方便，且对于低信噪比测量信号，可以获得更高的滤波精度。上述说明仅是本专利技术技术方案的概述，为了能够更清楚了解本专利技术的技术手段，并可依照说明书的内容予以实施，以下以本专利技术的较佳实施例并配合附图详细说明如后。本专利技术的具体实施方式由以下实施例及其附图详细给出。附图说明此处所说明的附图用来提供对本专利技术的进一步理解，构成本申请的一部分，本专利技术的示意性实施例及其说明用于解释本专利技术，并不构成对本专利技术的不当限定。在附图中：图1为本专利技术基于时域与频域的卡尔曼平滑滤波方法及滤波器的信号处理流程图；图2为本专利技术基于时域与频域的卡尔曼平滑滤波方法及滤波器的卡尔曼平滑滤波器结构示意图；图3为本专利技术基于时域与频域的卡尔曼平滑滤波方法及滤波器的卡尔曼平滑滤波算法(循环计算)示意图；图4为本专利技术基于时域与频域的卡尔曼平滑滤波方法及滤波器的噪声统计特性参数估算框图示意图；图5为本专利技术基于时域与频域的卡尔曼平滑滤波方法及滤波器的卡尔曼平滑滤波器模块结构示意图。具体实施方式以下结合附图对本专利技术的原理和特征本文档来自技高网...

【技术保护点】
1.基于时域与频域的卡尔曼平滑滤波方法，其特征在于，包括：/nS1，获取待滤波信号；/nS2，采用泰勒级数低阶近似模型，根据信号采样频率和模型维数，确定平滑滤波器的状态方程和量测方程；/nS3，根据卡尔曼平滑滤波算法，确定平滑滤波器的频率响应函数，进而可得到其幅频特性函数；根据所设计平滑滤波器的幅频特性参数，估算系统噪声和量测噪声的统计特性参数；/nS4，根据卡尔曼滤波模型及其噪声统计特性参数，采用卡尔曼平滑滤波算法，对所述待滤波信号进行平滑滤波处理。/n

【技术特征摘要】
1.基于时域与频域的卡尔曼平滑滤波方法，其特征在于，包括：
S1，获取待滤波信号；
S2，采用泰勒级数低阶近似模型，根据信号采样频率和模型维数，确定平滑滤波器的状态方程和量测方程；
S3，根据卡尔曼平滑滤波算法，确定平滑滤波器的频率响应函数，进而可得到其幅频特性函数；根据所设计平滑滤波器的幅频特性参数，估算系统噪声和量测噪声的统计特性参数；
S4，根据卡尔曼滤波模型及其噪声统计特性参数，采用卡尔曼平滑滤波算法，对所述待滤波信号进行平滑滤波处理。


2.根据权利要求1所述的基于时域与频域的卡尔曼平滑滤波方法，其特征在于，采用泰勒级数低阶近似模型，滤波器的状态方程为：
Xk＝ΦXk-1+ΓWk-1(1)
式中，Wk＝[Wk]
Xk为tk时刻的状态向量,Φ为tk-1时刻至tk时刻的一步转移矩阵,Γ为系统噪声耦合矩阵,Wk为零均值白噪声或高斯白噪声，其协方差矩阵为Qk；
其中，Xk中包含xk和xk的m阶导数(m＝1,2,…,n-1)的状态信息序列，表示xk的一阶导数，ts为采样间隔，n为模型维数。


3.根据权利要求1所述的基于时域与频域的卡尔曼平滑滤波方法，其特征在于，采用泰勒级数低阶近似模型，滤波器的量测方程为
Zk＝HXk+Vk(2)
式中，Zk＝[x′k]，H＝[100...0]，Vk＝[Vk]，Zk为tk时刻的量测向量，x′k为待滤波信号，H为量测矩阵，Vk为零均值白噪声或高斯白噪声，其协方差矩阵为Rk。


4.根据权利要求1所述的基于时域与频域的卡尔曼平滑滤波方法，其特征在于，所述卡尔曼平滑滤波算法包括正向滤波算法，正向滤波算法如下：
状态一步预测
状态估计
滤波增益Kk＝Pk/k-1HT(HPk/k-1HT+Rk)-1(5a)
或
一步预测均方误差Pk/k-1＝ΦPk-1ΦT+ΓQk-1ΓT(6)
估计均方误差
或Pk＝(I-KkH)Pk/k-1(7b)
或
式中,I为单位矩阵，给定初始状态向量和初始估计均方误差矩阵P0，按(3)～(7)式，(7)式包括(7a)、(7b)及(7c)，其他的标号也一样，由量测向量{Zk}(k＝1,2,…,M)，递推得到状态向量估计表示滤波器Xk的先验估计，表示滤波器Xk-1的后验估计，Zk为待滤波信号，Kk表示卡尔曼滤波增...

【专利技术属性】
技术研发人员：王静波，熊盛青，罗锋，郭志宏，周锡华，王冠鑫，
申请(专利权)人：北方工业大学，中国自然资源航空物探遥感中心，
类型：发明
国别省市：北京;11