联邦卡尔曼滤波

2024-08-24

联邦卡尔曼滤波（精选8篇）

联邦卡尔曼滤波篇1

0 引言

在现代战争实际环境里, 通常测得的目标特征数据是非常有限的, 而目标的方位几乎成了唯一可靠的参数, 因而利用所测得的目标角度信息可以估计目标的运动参数 (位置、速度、加速度等) , 从而对目标实施有效的打击和电子干扰。同时, 多站定位系统通过组合来自多个传感器的数据, 解除了单站观测对传感器运动轨迹的限制, 可以完成单一传感器独自所不能完成的任务, 随着传感器技术和现代网络及通信技术的发展, 站间通信等问题己不再成为制约多站无源定位技术发展的瓶颈, 多站测角传感器被动定位[1,2,3,4]由于它的工程实现简单、测量信息多、作用距离远和定位精度高等特点, 已获得日益广泛的应用, 在侦察、救援、宇航、航空、航海等领域发挥着重要作用。

在现代传感器被动定位研究中, 常常把问题转化为附加噪声的非线性动态系统的滤波问题, 滤波模型可转化为离散时间的一对差分方程[5], 即对目标的状态方程和量测方程。但是, 随着现代传感器被动定位研究的进一步深入, 发现目标跟踪定位的研究经常面临两方面的困难:1) 低运算复杂度下目标跟踪精度提高的问题;2) 未知噪声先验信息下的滤波递推估计问题。

针对第一个问题, 近年来, 学者们提出了多种算法, 如利用一阶Taylor展开进行系统线性化的扩展卡尔曼滤波 (EKF) 算法[6], 基于UT变换的采样 (Sigma) 点的不敏卡尔曼滤波 (UKF) 算法[7]等。其中在运算复杂度相当的前提下, UKF算法的精度要高于EKF算法, 因而成为近年来应用极为广泛的非线性滤波算法。最近, 一种新兴的非线性滤波方法—容积卡尔曼滤波 (CKF) 算法[8]被提出, 它采用一组等权值的容积点集来解决贝叶斯滤波的积分问题, 为非线性估计问题提供了一种新的实现方式。为保证算法的稳定性和正定性, 本文运用CKF的一种鲁棒改进形式—SCKF作为传感器被动定位研究的主要算法。

另外, 在现代战争复杂电磁环境下, 对于传感器噪声的量测往往是非经验的, 如果通过预先假定目标的噪声先验统计信息, 将给实际多站传感器被动定位带来误差。这是因为在先验信息不完全的情况下, 采用与真实值相差较多的噪声统计进行滤波, 将得到不可靠的状态滤波结果。同时, 在实际应用中, 动态系统随时可能受到各种外界干扰, 因此观测噪声的统计特性随时会发生变化, 仅仅依靠有限的先验信息对噪声进行描述是不可靠的, 需要一种能够对噪声进行实时估计的算法。基于这个原因, 本文引入一种经验模态分解 (EMD) 算法[9]对目标噪声的先验统计信息进行实时的估计, 以满足实际传感器被动定位跟踪需要, 同时, 将目标的过程噪声协方差和量测噪声协方差融入到整个循环过程中, 进一步提高了算法精度。

1 贝叶斯滤波过程

对于目标的跟踪其本质是一个贝叶斯滤波循环过程[10], 具体如下图1所示。对于目标的跟踪转化为如何求解目标状态的后验概率密度的问题。为得到时刻k的状态, 用下面的两个基本步骤更新 (k-1) 时刻状态的后验概率密度。

a) 时间更新, 即计算预测概率密度:

这里Dk-1 (28) {u i, zi} (k-1) i (28) 1表示对时刻k-1的输入量测的历史值;p (xkxk-1, uk-1) 表示已知的目标状态转移概率密度;p (xk-1Dk-1) 是时刻k-1的后验概率密度, 由式 (1) 获得。

b) 量测更新, 即计算当前状态的后验概率密度:

使用贝叶斯准则得到:

这里的正态常数ck由下面的式子得到:

建立式 (2) 先验和后验概率的递归关系, 输入需要满足下面的关系:

这里的控制Dk-1要有足够的信息产生输入uk的充分信息。具体的, 输入uk用xˆk|k-1产生。在这种情况下, 可以等价的写出:

因此, 把式 (6) 代入式 (3) 中, 得到:

根据假设

但在实际情况中, 当面对一个非线性滤波问题时, 很难得到贝叶斯后验概率密度的最优解, 因此需要寻求对贝叶斯滤波的次优解。在计算过程中, 对于后验概率密度的次优解可以通过下面两个方法获得:

1) 全局方法, 即利用抽样的方式进行贝叶斯递推。例如, 高斯混合滤波方法, 使用重要性抽样蒙特卡罗运算的粒子滤波[11]等, 都属于这第1) 类的情况, 但是这种全局方法要求的计算量很大。

2) 局部方法。即通过转变后验概率密度形式来达到非线性滤波的目的。例如, 假设它是高斯的非线性滤波方法, 如EKF (扩展卡尔曼滤波) , CDKF (中心差分卡尔曼滤波) [12], UKF (不敏卡尔曼滤波) , 都属于这种类型。对于局部信息的强调使得滤波过程简捷快速。

2 UKF (不敏卡尔曼滤波) 算法

作为局部滤波方法的代表之一, 不敏卡尔曼滤波 (UKF) 是一种建立在高斯域的逼近贝叶斯滤波器, 它的特点是使用一种称为UT变换的方式来进行递推运算, 考虑一个拥有均值为, 协方差为Σ的n维随机变量x, 它的 (2n+1) 个抽样点和权值的集合{i, i}2ni (28) 0, 满足:

在这里, 一个系统分布式点集, 被称为西格码点集, 表述如下:

参数用于尺度化来自先验均值的西格码点, 称为“尺度化参数”。为了紧密获得先验概率密度的峰值, 一般选择 (28) (3-n) 。

随着传感器被动定位方式的广泛应用, 既要实现对目标快速的定位, 又要尽量避免牺牲定位精度, 这就迫切需要挖掘更好的局部滤波方法, 来进行更高精度的传感器被动定位;现今对于UKF技术的研究已经日趋成熟, 人们逐渐把目光转移到追求更高跟踪精度的滤波方法上, 最近新提出的求容积法卡尔曼滤波 (CKF) 作为一种新的滤波方式, 提供了解决问题的新契机。

3 扩展SCKF算法

UKF算法在进行递推运算的过程中可能会面临的一个很大的问题是:如何很好地保持误差协方差以及相关运算的正定性, 同样, 在CKF算法的递推中也会遇到这样的问题, 使用求平方根的形式对CKF算法进行改进, 不但可以进一步提高算法的精度, 而且保证了算法进行过程中的正定性。

EMD (经验模态分解) 算法能够分离信号和噪声, 因此可以用来估计含噪声信号中噪声的标准差。在观测噪声未知的情况下, 可以首先选取一段长度为L的一段信号, 即相当于一个观测区间, 在这个观测区间上进行经验模态分解, 分离出其中的高频噪声, 估算出观测噪声的标准差。

另外, 本文还将状态噪声协方差和量测噪声协方差融入到整个循环过程中, 进一步提高算法准确性。

3.1 时间更新

1) 假设在k-1时刻的后验概率密度函数:p (xk-1Dk-1) (28) N (xk-1k-1, Pk-1k-1) 已知。分解得到:

2) 估计容积点的值 (i (28) , 1, 2, m) :

3) 估计传播的容积点

4) 估计预测状态:

5) 估计预测误差协方差的平方根因子:

其中:Tria () 代表通过QR分解取下三角阵, 这里的SQ, k-1代表Qk-1的一个平方根因子, 这里的Qk-1 (28) SQ, k-1SQT, k-1, 并且权值, 中心 (先验均值被减去) 的矩阵为

3.2 量测更新

1) 估计容积点

2) 估计传播容积点

对量测噪声做EMD估计, 具体步骤如下:

(a) 计算出Zi, k|k-1所有的极值点。

(b) 采用三次样条插值算法[13], 求解所有的极大值点构成的上包络线和所有的极小值点构成的下包络线, 分别记为u1 (t) 和v1 (t) ;为抑制边缘效应, 需使用镜像延拓的方法, 来延拓极值点[14,15]。

并记信号与上、下包络线的均值之差为

(d) 判断h11 () (t) 是否满足IMF (本征模式函数) 的两条性质。若满足, 则h11 () (t) 为IMF。否则, 记h11 () (t) 为Zi, k|k-1, 重复步骤 (a) 至步骤 (d) , 假定k次筛选后得到的h11 () (t) 满足IMF定义, 则Zi, k|k-1的第一阶IMF分量c1 (t) (28) h1 (1) (t) 。

(e) 记Z1 (t) (28) Zi, k|k-1-c 1 (t) 为新的待分析信号。重复步骤 (a) 至步骤 (e) , 以得到第二个IMF, 记为c2 (t) , 余项Z2 (t) (28) Z1 (t) -c2 (t) 。重复上述步骤, 直至得到的余项为一个单调信号, 分解结束。如此, 可得到余项Zm (t) , Zm (t) 即对量测噪声的观测, 从而确定新的Rk。

3) 估计预测量测

4) 估计更新协方差矩阵的平方根形式:

这里的SR, k表示kR的一个平方根因子, 即Rk (28) SR, kSRT, k, 并且权值, 中心矩阵为

5) 估计互协方差矩阵:

这里的权值, 中心矩阵为

6) 估计卡尔曼增益:

7) 估计更新状态

8) 估计对应的误差协方差的平方根因子:

进一步得到Pk|k (28) Sk|kSTk|k。

9) 将过程噪声协方差和量测噪声协方差纳入到整个循环过程中:

令新的协方差矩阵为重新代入整个循环中, 并将式 (18) 、式 (25) 分别更改为

4 仿真实验

4.1 仿真条件

以三雷达传感器的目标定位为例, 如图2所示, 三个传感器之间可以相互通信并随时获得传感器之间的相对位置。

假设的实验坐标系为直角坐标系, 总的实验时间为800 s, 一共做了5次实验, 目标的初始位置设定在 (120, 100) , 状态噪声设定为均值为0, 方差为0.01的高斯白噪声, 三个传感器的初始位置分别为 (-100, 0) , (0, 0) 和 (100, 0) 。量测噪声设定为高斯白噪声, 均值为0, 方差取10-4~10-6的随机值。

假设初始的系统误差方差阵为

对目标和传感器的运动分别做了两次假设:

1) 目标和传感器都做匀速运动, 目标在x, y轴方向上的速度都为0.3 m/s, 三个传感器速度均为3m/s时的纯角度传感器跟踪;

2) 目标和传感器都做匀速运动, 目标在x, y轴方向上的速度都为3 m/s, 三个传感器速度均为30 m/s时的纯角度传感器跟踪。

4.2 仿真结果

这里只显示传感器对目标x方向速度的跟踪结果, 如下图36所示。

其中, 图3、图4代表按照第1) 种假设所做的目标速度跟踪以及均方差比较仿真;图5、图6代表按照第2) 种假设所做的目标速度跟踪以及均方差比较仿真。

由图3、图4可以看出, UKF和扩展SCKF都能对目标x方向的速度进行跟踪, 但扩展SCKF的跟踪精度明显优于UKF算法。

当目标速度和传感器速度都增大时, 目标跟踪的精度将会受到很大的影响, 如图5、图6所示。

由图5、图6可以看出, 在UKF对目标的跟踪出现滤波发散的情况下, 扩展SCKF算法仍可实现对目标的稳定跟踪。

扩展SCKF算法的运算时间主要消耗在对于量测噪声统计信息的EMD估计中, 图7是对两种算法的运算复杂度的对比, 即对目标定位做100次蒙特卡罗仿真的时间比较, 由图可以看出, 扩展SCKF算法所消耗的运算时间与UKF相当。

5 结论

在多传感器被动定位问题中, 对于目标状态的递推跟踪算法的研究一直是人们观注的热点, 其本质是求解贝叶斯过程的后验概率密度问题。本文以三运动测角传感器的目标跟踪为背景, 提出运用一种新的非线性滤波算法CKF的改进形式扩展容积卡尔曼滤波算法对目标的速度进行递推跟踪, 仿真结果表明, 当目标速度与传感器速度选取不同值时, 相对于传统的UKF算法, 在运算复杂度相当的前提下, 扩展容积卡尔曼滤波算法不仅可以对未知量测噪声的目标进行跟踪, 而且显著提高了目标的跟踪精度。

摘要：为提高被动定位技术的精度与环境适应性, 本文提出运用一种新的非线性滤波方法—扩展容积卡尔曼滤波算法进行多角度传感器目标定位;它首先利用EMD (经验模态分解) 算法对目标的量测噪声协方差矩阵进行估计;然后, 将过程噪声协方差和量测噪声协方差融入循环过程;同时, 为保持算法的稳定性和正定性, 利用求平方根的形式对算法改进。通过对扩展容积卡尔曼滤波与UKF (不敏卡尔曼滤波) 算法跟踪目标的结果进行比较, 在运算复杂度与UKF相当的前提下, 扩展容积卡尔曼滤波算法不仅可以对未知量测噪声情况下的目标进行跟踪, 而且显著提高了被动定位的精度。

关键词：测角传感器,容积卡尔曼滤波,经验模态分解,被动定位

联邦卡尔曼滤波篇2

针对采用标准卡尔曼滤波器必须知道系统噪声统计特性的.局限性,研究了一类系统噪声未知情况下的自适应联邦滤波方法,指出了自适应滤波方法应用于联邦结构时应当注意的问题,提出了一种基于信息补偿的自适应联邦滤波算法.SINS/BDS/GPS组合导航系统的仿真结果表明,该方法可以有效抑制系统噪声未知情况下的滤波发散现象,提高了滤波的稳定性和估计性能.

作者：邱恺荣军陈天如吴训忠 Qiu Kai Rong Jun Chen Tianru Wu Xunzhong 作者单位：邱恺,Qiu Kai(北京航空工程技术研究中心,北京,100076;空军工程大学工程学院,西安,710038)

荣军,Rong Jun(空军指挥学院后勤与装备系,北京,100089)

陈天如,吴训忠,Chen Tianru,Wu Xunzhong(空军工程大学工程学院,西安,710038)

卡尔曼滤波算法的硬件实现研究篇3

卡尔曼滤波算法是一种解决离散系统线性滤波问题的递推最优估计算法,它是采用状态空间法在时域内设计的滤波器。20世纪60年代以来,卡尔曼滤波算法得到了深入研究和广泛的工程应用,而今,它已经成功应用在导弹制导和控制、飞机自动驾驶、卫星导航[1,2]与定位等无线电技术的诸多领域,并取得了非凡的成就,为军事、工业以及科学工程的发展做出了重大的贡献。

在实际应用中,需要考察卡尔曼滤波算法的运算精度[3]和运算速度。卡尔曼滤波算法常采用通用处理器通过C语言软件编程实现。由于卡尔曼滤波算法存在大量的计算过程,同时通用的处理器串行地执行软件程序,因此在对实时性要求比较高的场合中该方法难以满足系统对运行速度的要求[4]。针对通用处理器实现卡尔曼滤波算法执行速度较慢的问题,考虑到现场可编程门阵列(Field Programmable Gate Array,FPGA)具有的高度并行性、执行速度快的特点,近些年来对在FPGA中设计卡尔曼滤波算法的硬件实现展开了初步的研究[5,6]。本文借鉴了通用处理器的CPU、RAM以及外设三大部件组成架构,设计了一种专用于处理卡尔曼滤波器的并行处理状态机,提出了一种基于FPGA的卡尔曼滤波算法的硬件实现方法。

1 卡尔曼滤波算法原理

卡尔曼滤波算法以“预测—实测—修正”的过程递推运算,根据系统的测量值来消除随机干扰,得到最优估计值。假设系统的状态方程[7]为:

系统的测量方程为:

式中,X(k)为k时刻的系统状态;U(k)为k时刻对系统的控制量;F和B为系统参数;Z(k)为k时刻的测量值;H为测量系统的参数;W(k)和V(k)分别为过程和测量的高斯白噪声,其协方差分别是Q和R。

利用系统的过程模型,可预测系统的下一状态为:

式中,X(k|k-1)是利用上一状态预测的结果;X(k-1|k-1)是上一状态最优的结果;U(k)为现在状态的控制量。

用P表示协方差,则可预测下一状态的协方差为:

式中,P(k|k-1)为X(k|k-1)对应的协方差;P(k-1|k-1)为X(k-1|k-1)对应的协方差;F'为F的转置矩阵;Q为系统过程的协方差。

根据状态的预测结果,计算得到目前状态的测量值,结合预测值和测量值,可计算得到状态k的最优化估算值X(k|k)表示为:

式中,Kg为卡尔曼增益,其计算方法为:

最后更新k状态下X(k|k)的协方差表示为:

式中,I为单位阵,如此递推下去。

2 基于FPGA的卡尔曼滤波算法设计

FPGA内部主要包括可配置逻辑模块(Configurable Logic Block,CLB)、输入输出模块(Input Output Block,IOB)和内部连线(Interconnect)3部分,而基于FPGA的卡尔曼滤波算法实现主要由这3部分构成,其实现原理如图1所示。

卡尔曼滤波算法数据输入输出由FPGA的IOB模块负责完成,算法的控制、运算和存储通过CLB模块设计实现,指令传达以及数据传输由Interconnect模块完成。

通过在FPGA的CLB模块中设计控制单元、运算单元和RAM存储器,可实现卡尔曼滤波器的递推公式运算过程。卡尔曼滤波算法控制单元,按照时序电路设计,依次完成状态预测、协方差预测、状态更新、卡尔曼滤波增益更新以及协方差更新。运算单元负责完成数值的加减乘除基本运算,其中包括加法器IP核、减法器IP核、乘法器IP核和除法器IP核。存储单元由一块RAM存储器构成,负责存储初始数据、中间数据及结果数据。

在卡尔曼滤波算法运算过程中涉及多个矩阵,其中包括状态更新矩阵、观测矩阵、协方差矩阵、状态向量、状态预测向量、协方差矩阵、协方差预测矩阵和卡尔曼滤波增益矩阵等,它们都需要存在一块RAM存储器中,因此在设计过程中为每一个矩阵分配一定的存储空间及相应的存储地址。

在应用基于FPGA的高维卡尔曼滤波算法的过程中,首先对RAM存储器中的滤波算法相关矩阵进行初始化,然后读取观测量,在控制单元的时序设计下递推完成卡尔曼滤波算法,其状态机执行流程如图2所示。

基于FPGA的卡尔曼滤波算法状态机由6个状态组成,其中S1、S2、S3、S4和S5依次通过控制单元执行时序控制过程完成,每个控制过程所涉及的矩阵运算由S6通过运算单元处理完成,在运算过程所需的计算参数从存储单元读取,并把每一控制过程的运算结果存放在存储单元。

3 实现验证

ISE是使用XILINX的FPGA的必备的设计工具。它可以完成FPGA开发的全部流程,包括设计输入、仿真、综合、布局布线、生成BIT文件、配置以及在线调试等,功能非常强大。本论文的实现验证工作在ISE开发工具中完成[8]。

首先在ISE开发环境中,建立工程,在该工程中,添加卡尔曼滤波算法硬件实现所需的RAM存储器、乘法器IP核、除法器IP核、加法器IP核以及减法器IP核。卡尔曼滤波器实现的FPGA工程如图3所示。

其中,RAM存储器的配置界面如图4所示。

由图4可以看出,RAM存储器采用Block Memory模块,因为在数据运算过程中采用双精度浮点运算,所以定义数据宽度为64位,存储深度依据卡尔曼滤波器维数而定。

运算单元的IP核配置界面如图5所示。

由图5可以看出,运算单元的所有IP核采用双精度浮点Floating-point模块,该模块以占用大量的逻辑资源为代价换取高速的运算速度。

关键的卡尔曼滤波算法采用VHDL语言编程实现[9]。如上描述,卡尔曼滤波算法的运算过程可以理解为以下循环的重复和逻辑控制:①从RAM中读取数据;②把数据送至运算器进行处理;③处理的数据结果写入RAM。如此便可完成卡尔曼滤波器的硬件实现。

卡尔曼滤波算法在通用处理器以及FPGA中实现的运行速度对比结果如表1和图6所示。

表1和图6中,通用处理器中采用C语言设计的卡尔曼滤波算法,工作频率为100 MHz;同时,FPGA中采用VHDL设计的卡尔曼滤波算法,FPGA的工作频率为100 MHz。这保证了通用处理器和FPGA工作在相同的时钟频率下。

由表1中的对比数据可以看出本论文所提出的基于FPGA的卡尔曼滤波算法硬件实现方法切实可行,而且运行速度比在通用处理器中执行速度快。

4 结束语

本文针对卡尔曼滤波算法在通用处理器中串行执行软件程序运行速度低的问题,提出了一种基于FPGA的卡尔曼滤波算法的硬件实现方法,并通过ISE工具进行了开发验证,验证结果随着卡尔曼滤波器维数的增加,基于FPGA的卡尔曼滤波算法将比在通用处理器中的运算速度大大提升。但本方法尚未在工程实践中得到应用,这有待于继续探索。

摘要：卡尔曼滤波算法常采用通用处理器串行执行软件程序的工作方式来实现,但在对实时性要求比较高的场合,这种方式往往不能够满足系统对运行速度的要求。针对此问题,提出了一种基于FPGA的卡尔曼滤波算法的硬件实现方法,该实现方法具有并行运行的特点,大大提高了卡尔曼滤波算法的执行速度,具有良好的运行性能。研究结果表明,该方法切实可行,且运行速度比在通用处理器中执行速度快。

关键词：卡尔曼滤波,FPGA,并行

参考文献

[1]周富相,郑晓晶.基于卡尔曼滤波的RAIM方案研究[J].无线电工程,2011,40(11):36-39.

[2]余小游,高亭亭,孙广富,等.卫星导航弱信号的变维卡尔曼滤波跟踪算法[J].国防科技大学学报,2015,3(3):56-60.

[3]梁民赞,黄子豪,曹占启.曲线拟合与卡尔曼滤波器的滤波精度评估[J].无线电工程,2013,43(3):36-39.

[4]CHUI C K,CHEN G.卡尔曼滤波及其实时应用[M].北京:清华大学出版社,2013:1-3.

[5]冯刚,吕茂庭,覃天.基于MATLAB的卡尔曼滤波仿真研究[J].电子产品可靠性与环境试验,2011,29(6):61-63.

[6]徐经纬,刘士君.交互多模型扩展卡尔曼滤波算法的FPGA实现[J].传感器与微系统,2014,33(1):11-14.

[7]李朝伟,周希元,沈文亮.Kalman滤波算法在被动跟踪中的应用[J].无线电工程,2007,39(5):49-51.

[8]田耘,徐文波.Xilinx FPGA开发实用教程[M].北京:清华大学出版社,2008:14-18.

联邦卡尔曼滤波篇4

针对联邦滤波器子系统传感器故障,提出了一种基于状态递推器的滑动残差检验法,并在故障检测函数的基础上,定义了故障因子函数,用于联邦滤波器的信息分配.根据新的算法,联邦滤波器的信息分配系数可根据故障因子函数进行自动调整,以减少故障通过全局重置对其他无故障子系统的污染程度,从而减少了故障信息的影响程度和持续时间,提高了故障隔离后联邦滤波器的.快速恢复能力.仿真结果表明该算法是有效的.

作者：吴训忠周军邱恺 WU Xun-zhong ZHOU Jun QIU Kai 作者单位：吴训忠,WU Xun-zhong(西北工业大学航天工程学院,西安,710072;空军工程大学工程学院,西安,710038)

周军,ZHOU Jun(西北工业大学航天工程学院,西安,710072)

邱恺,QIU Kai(空军工程大学工程学院,西安,710038)

扩展卡尔曼滤波的影响因素分析篇5

本文介绍了扩展卡尔曼滤波算法,讨论了算法中各因素对于扩展卡尔曼滤波的影响;最后给出其选择原则,旨在全面了解和使用扩展卡尔曼滤波。

1 扩展卡尔曼滤波(EKF)算法

扩展卡尔曼滤波(EKF)是对非线性系统状态方程和量测数据方程进行线性化近似后再采用卡尔曼滤波。通过求过程方程和量测数据方程的偏导数来对量测数据求其线性相近。式(1)表示非线性差分方程和量测数据方程

${\begin{cases} x_{k} = f (x_{k - 1}, k - 1) + w_{k - 1} \\ z_{k} = h (x_{k}, k) + v_{k} \end{cases} (1)$

算法可以分为预测和更新两个阶段:

预测:

${\begin{cases} \hat{x}_{\bar{k}} = f (x_{k - 1}, u_{k - 1}, w_{k - 1}) \\ p_{\bar{k}} = A_{k} Ρ_{k - 1} A_{k}^{Τ} + W_{k} Q_{k - 1} W_{k}^{Τ} \end{cases} (2)$

更新:

${\begin{cases} Κ_{k} = Ρ_{k} - Η_{k}^{Τ} (Η_{k} Ρ_{\bar{k}} + V_{k} R_{k} V_{k}^{Τ})^{- 1} \\ \hat{x}_{k} = \hat{x}_{\bar{k}} - Κ_{k} (z_{k} - h (\hat{x}_{k}, 0)) \\ Ρ_{k} = (Ι - Κ_{k} Η_{k}) Ρ_{\bar{k}} \end{cases} (3)$

式中,Ak,Wk分别是f对x和y求偏导的雅可比矩阵;Hk,Vk分别是对和求偏导的雅可比矩阵;Qk,Rk分别是过程激励噪声和量测数据噪声的协方差矩阵;Kk是卡尔曼滤波增益。扩展卡尔曼滤波算法中Kk表达式中的雅可比矩阵Hk能正确地传递加权量测数据信息中的有用部分。

2 参数分析

由式(3)可以看出,除去一些确定的因素和中间值,影响滤波的主要有5个因素:初始状态 $X_{0} = [\begin{matrix} x_{0} v_{x_{0}} & y_{0} v_{y_{0}} \end{matrix}]^{Τ}$ 、初始状态协方差阵P0、状态转移矩阵A、过程激励噪声协方差阵Q、测量噪声协方差R。下面对影响滤波的5个主要因素进行分析。

首先给出一个用以分析的例子,下面的分析都是建立在此例下:假设一个机动目标的状态变量为 $X_{k} = [\begin{matrix} x_{k} v_{x_{k}} & y_{k} v_{y_{k}} \end{matrix}]^{Τ}$ ,其中第1、3项为k时刻的坐标,第2、4项为k时刻的速度。考虑非线性非高斯系统模型,为简化计算取匀速直线运动模型,即转移矩阵Φ一定;过程噪声和观测噪声都为零均值高斯噪声,物体的状态初始值为 $X_{0} = [\begin{matrix} 100 & - 5 100 & - 5 \end{matrix}]^{Τ}$ ,即位置坐标为(100,100),速度为-5 m/s。物体以5 m/s的速度做了50个周期的匀速运动,观测间隔为1;观测误差期望为0,方差为5的高斯白噪声。

2.1 初始状态X0

实际工程应用中,目前较为精确的GPS定位精度可达15 m以内,而被动定位精度可能更低;由于无法精确得到目标位置的先验信息,而使滤波初值可能有较大的误差。按照定位精度20 m为参考,则在精度20 m范围内,X0相对误差不超过20%;取选取初始状态的相对误差为变量进行滤波,结果如图1所示。

由图1可以看出:初始状态对滤波效果影响较大,为使滤波得到较低的误差,最好能控制初始状态的相对误差在12%以内。误差小的同时还有利于加快滤波算法的稳定。

2.2 初始状态协方差阵P0

$Ρ_{0} = [\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{matrix}] (4)$

协方差是用于衡量两个变量的总体误差;由于定位精度的限制,按照(间隔0.5)选取初始状态协方差,并分析各滤波效果。进行滤波,结果如图2所示。

由图2可以看出:滤波随着协方差的变化很快趋于稳定,所以可以看出初始状态协方差阵对滤波效果影响很小,都能较快收敛,可以任意取一个不为零的矩阵,如取单位对角阵。

2.3 过程激励噪声协方差阵Q

Q=10-3×diag([1]) (5)

过程激励噪声协方差是指用于估计离散时间过程的状态变量。这个离散时间过程中的w被称为过程激励噪声。根据白噪声的假定,Q的非对角线元素均为零;工程上一般通过离线测试来确定对于某过程的最优Q值[3,5]。根据文献可知过程激励噪声协方差一般都较小,所以按照0～1 000(间隔1)选取过程激励噪声协方差阵,并分析各滤波效果,结果如图3所示。

由图3可以看出:滤波结果的误差随着过程激励噪声协方差的增大而增大,且变换较为显著;所以过程激励噪声协方差阵的取值越小越好,可以使用一个非常小但不为零的矩阵,可以取Q=10-3×diag([1])或Q=10-4×diag([1])。

2.4 测量噪声协方差R

测量噪声协方差是滤波过程中的一个重要参数,是滤波中的一个重要输入值。但很多情况下在滤波前未必能够准确获得该值,因此有必要研究R的不同取值对滤波的影响。由于受工程上定位精度的限制,在定位精度范围内,测量噪声协方差一般不会超过30,所以按照0～30(间隔1)选取测量噪声协方差阵,并分析各滤波效果。结果如图4所示。

由图4所示:测量噪声协方差R取值偏小,滤波出现误差突然变大的情况;相反,如果取值偏大那么滤波误差变换较为缓慢。这可能是由于滤波次数的限制导致偏差的R值未能使滤波收敛到有效范围。而随着测量噪声协方差R的变化,滤波误差变化较大;一般工程上是在滤波前先测定噪声协方差,然后再用于后续的滤波。

一般工程上为了简化计算,滤波前都先对测定噪声协方差和过程激励噪声协方差进行估计,然后再用于后续的滤波;而系统的动态噪声和观测噪声为零均值且统计特性已知的白噪声,在复杂的实际工作中未必能够满足,因此可能存在较大误差;为解决上述问题,文献[7]针对传统算法的固定设计在先验信息不充分和动态变化环境中存在的不足,提出了一种基于自适应算法。该算法通过计算各时刻滤波残差的变化,以实时修正过程激励噪声协方差和测量噪声协方差。但这样也会带来算法和时间复杂度的增加,采用何种算法需要按实际情况权衡决定。

3 结束语

根据以上对于滤波参数分析,在有限的滤波次数范围内,可以得出如下结论:

(1)初始状态X0对滤波效果影响较大,为使滤波得到较低的误差,最好能够控制初始状态的误差在12%以内。

(2)滤波随着协方差的变化很快趋于稳定,即初始方差矩阵P0对滤波效果影响很小,都能较快收敛,可以任意取一个不为零的矩阵,如取为单位对角阵。

(3)滤波结果的误差随着过程激励噪声协方差的增大而增大,且变换较为显著,所以过程激励噪声协方差阵的取值越小越好,可以使用一个非常小但不为零的矩阵,可以取Q=10-3×diag([1])或Q=10-4×diag([1])。

(4)测量噪声协方差R的取值变大或变小,其滤波效果不一定好;因此可以在滤波前先测定噪声协方差,用于后续的滤波;带较大偏差的测量噪声协方差将导致算法发散。

摘要：介绍了扩展卡尔曼滤波算法流程,并对影响滤波效果的5个主要参数进行了系统的讨论,最后通过仿真实验方法研究不同的参数变化对于滤波的影响,并给出了参数的选择原则。

关键词：扩展卡尔曼滤波,参数选择,非线性系统状态

参考文献

[1]DAUM F.Nonlinear filters:beyond the Kalman filter[J].IEEE AES Systems Magazine,2005,20(8):57-69.

[2]ZHU Jihua,ZHANG Nanning,YUAN Zejian,et al.A SLAMalgo-rithm based on the central difference Kalman filter[C].Intelligent Vehicles Symposium.IEEE,2009:123-128.

[3]王学斌,徐建宏,张章.卡尔曼滤波器参数分析与应用方法研究[J].计算机应用与软件,2012(6):212-215.

[4]余小琴,沈文苗.扩展卡尔曼滤波算法初值选取方法[J].声学与电子工程,2012(1):12-13.

[5]张谦,景占荣.一种过程噪声自适应调节的卡尔慢滤波算法[J].电子测量技术,2007(5):18-20.

[6]李良群,姬红兵,罗军辉.迭代扩展卡尔曼粒子滤波器[J].西安电子科技大学学报:自然科学版,2007,34(2):233-238.

基于卡尔曼滤波的动态矩阵控制篇6

模型预测控制(Model Predictive Control)是一类基于模型的控制算法总称,它的提出源于工业实践。在各种模型预测控制算法中,动态矩阵控制(DMC)是最为经典的算法,DMC是一种基于阶跃响应系数模型的MPC算法,采用增量式可有效地消除系统的稳态误差,适用于渐进稳定的线性对象。但是DMC在处理不可测干扰时具有明显的局限性,如果系统中存在不可测干扰,闭环系统很难实现无静差控制[1,2],因此,对不可测干扰的抑制,成为DMC取得理想控制效果的关键。

文献[3]运用自校正DMC对气温控制系统仿真,但仅仅考虑了阶跃干扰对系统的影响,并没有讨论干扰不可测的情况;文献[4]运用卡尔曼滤波去除系统中的噪声,对倒立摆系统进行了控制,但只考虑了白噪声对倒立摆的影响;文献[5]应用DMC对气体分离装置进行控制,在干扰可测的情况下,能够实现对可测扰动的理想抑制,但是当扰动不可测时,无法实现气体分离装置的理想控制。

基于此,本文采用状态空间模型对不可测干扰建模,通过卡尔曼滤波估计出该不可测干扰的值,按前馈的方法极大地克服了不可测干扰对DMC控制系统的影响。

1 DMC控制器

DMC的推导方法并不是唯一的,其基本思想是模型预测、在线反馈校正、滚动优化。文献[6]和文献[7]分别给出了不同的DMC推导方法,本文选取文献[7]中的DMC作为电机控制器。

DMC在单位阶跃输入作用下,时不变单入单出(SISO)系统的输出响应为{0,s1,s2,…,sN,SN+1,}。

这里假设系统输出恰好在变化N步后达到稳态,向量s=(s1,s2,…,sN)T称为模型向量,N则称为建模时域。

考虑开环稳定系统,从时刻k开始,在M个控制增量△u(K),△u(k+1)|k),…,△u(k+M-1)|k)的作用下,使被控对象在未来p个时刻的输出预测值为:

为使问题有意义,通常规定M≤P≤N。式中为动态矩阵，是由阶跃响应系数s1组成的P×M阵。为模型输出预测值,Y0(k+1)|k-1)为假设当前和未来时刻控制作用不变时的输出预测值,△U(k)为控制增量向量,分别表示:

由于受模型误差和干扰等影响,系统的输出预测值需在预测模型输出的基础上用实际输出误差进行在线反馈校正,以实现闭环预测。反馈校正后的预测输出为:

假设优化的准则是最小化如下性能指标:

其中e(k+i|k)=ys(k+i)-y(k+i|k)为跟踪误差;为未来输出参考值(设定值);W和R分别为由权系数wi和rj构成的对角阵,称为误差权系数矩阵和控制权系数矩阵。

使J(k)取极小的,可以通过极值必要条件求得(当ATWA+R为可逆矩阵时):

其中。

式(7)给出了△u(k),△u(k+1)|k),…,Du(k+M-1|k)的最优值。但DMC只取其中的即时控制增量△u(k)构成实际控制u(k)=u(k-1)+△u(k)作用于对象,到下一时刻,又求解类似的优化问题,得到△u(k+1)。这就是所谓的“滚动优化”的策略。

2 无刷直流电机模型

无刷直流电机系统的动力学微分方程为:

式中各参数如下[8]:电枢电阻Ra=3.72Ω,电枢电感La=7.83mH,粘滞摩擦系数B=2.59×10-4N m/rad/s,转动惯量,电磁系数Kb=0.31V/ras/d,Va为电枢电压,ia为电枢电流,w为电机角速度,θ为电机转角,Tl为负载转矩,为了方便研究干扰对电机的影响,令Tl=0。选取状态变量X=[i w θ]T,u=Va,输出Y=w由式(8)得到系统的状态空间模型为:

其中:

3 不可测扰动抑制方法设计

3.1 不可测扰动的建模与估计

对线性模型(9)进行离散化,考虑到系统中的不可测干扰d,把d可以作为系统的一个状态变量增广到式(9)中,输入u(k)=Va,状态变量X(k)=[i(k)w(k)θ(k)]T,输出Z(k)=w(k),则带有不可测输出干扰的电机离散模型为:

其中,sp是增广输出干扰状态的个数;Gp决定这些状态对输出的影响。通常取Gp=I,使用公式dk=Z(k)-HX(k)来估计输出扰动。其结果是:对于输出干扰状态,这是一个最小拍观测器;对于系统状态,这是一个开环观测器[1]。取Ts=0.008s,可得:

3.2 基于卡尔曼滤波的状态估计

系统控制方案如图1所示。

根据建立的不可测扰动模型,使用卡尔曼滤波对增广系统状态模型进行估计,考虑到系统中存在的过程噪声和测量噪声,令:

得到带有噪声的电机离散模型为:

其中w(k)为过程噪声,v(k)为测量噪声,它们是均值为零且互不相关的白噪声序列。这样,增广状态的估计可使用如下方法得到:

状态滤波增益k被分解为过程模型状态滤波增益k1和输出扰动滤波增益k2,每个周期的增益由下式得到:

其中,。Q和R分别为过程噪声协方差阵和观测噪声协方差阵。

基于增广的状态空间模型,通过卡尔曼滤波器在每个采样周期估计干扰d的大小之后,充分利用这部分信息进行前馈补偿,降低甚至消除干扰d对系统输出的影响。

4 系统仿真

根据系统控制方案,在MATLAB/simulink中搭建系统仿真模型。仿真过程中,假设过程噪声和测量噪声为具有时变统计特性的平稳高斯白噪声序列,取过程噪声v均值为零,方差Q=0.05,测量噪声均值为零,方差R=2。本文假设电机转速须按照以下规律变化:

DMC控制器参数设置为:建模时域N=10,控制时域M=4,优化时域P=5,误差权阵W=Ip×q,控制权阵R=0.5IM×M,采样周期为0.01秒。当不可测干扰d是幅值为10,频率为1rad/s的方波信号时,仿真结果如图2所示。

图2中的真实值为未加任何噪声下DMC控制电机转速的结果,可见,极好的实现了电机转速的预期目标。当加入噪声,未加Kalman滤波器时,系统受到过程噪声、测量噪声以及方波噪声的影响极大,完全偏离了预期控制目标。由于存在输出方波干扰噪声,加入Kalman滤波器后,虽然能够消除过程噪声、测量噪声的影响,小幅度的消除方波噪声的干扰,但是和控制目标存在着较大的恒值差,可见,由于存在方波干扰,系统无法实现无静差控制。当通过Kalman滤波器估计出方波信号的值,再进行前馈补偿后,很大程度上克服了方波干扰,实现了无静差控制。

估计方波噪声值和真实方波噪声值的比较如图3所示。由图3可以看出,通过Kalman滤波器估计的干扰信号较好的跟踪了实际方波干扰,估计值和实际值存在的误差较小。

5 结论

本文以直流电机的转速为研究对象,以干扰模型为基础,设计了一种基于卡尔曼滤波器的前馈补偿控制方法。该方法有效的克服了不可测干扰对直流电机转速的影响,改善了速度跟踪的精度,实现了无静差控制。仿真结果表明,该方法显著提高了DMC的控制性能。

摘要：动态矩阵控制(DMC)在处理不可测干扰时具有明显的局限性,如果系统中存在不可测干扰,闭环系统就达不到所期望的控制目标。本文将不可测干扰包含在过程模型中,通过卡尔曼滤波估计出不可测干扰,利用前馈方法极大地克服了不可测干扰对闭环系统的影响。将此控制方法应用于直流电机的仿真试验,仿真结果表明,该方法能够有效地抑制不可测干扰,在DMC控制下,电机转速能够快速达到期望值。

关键词：动态矩阵控制,不可测干扰,卡尔曼滤波,直流电机

参考文献

[1]Kenneth R.Muske,Thomas A.Badgwell.Disturbancemodeling for offset-free linear model predictive control[J].Journal of Process Control,2002,12(5):617-632.

[2]邹涛,丁宝苍,张端.模型预测控制工程应用导论[M].北京:化学工业出版社,2010:140-154.

[3]柳其亮,张福斌.自校正DMC-PID过热汽温控制系统仿真[J].计算机仿真,2008,25(12):254-257.

[4]祁虔,李祖枢,谭智.基于卡尔曼滤波的倒立摆控制系统噪声抑制[J].控制与决策,2010,8(25):1144-1148.

[5]牛建璋.基于Matlab/Simulink的多变量气体分离装置多变量预测控制[J].化工自动化及仪表,2010,37(5):29-34.

[6]席裕庚.预测控制[M].北京:国防工业出版社,1993:10-18.

[7]丁宝苍.预测控制的理论与方法[M].北京:机械工业出版社,2008:42-51.

[8]R.Krishnan.Electric Motor Drives:Modeling,Analysis,andControl[M].Englewood Cliffs,NJ:Prentice-Hall,2001.

联邦卡尔曼滤波篇7

在现实的目标运动分析中, 目标运动方程的建立一般都不是线性的, 特别是针对红外寻的制导导弹, 导引头仅能够测量角度或角速度信息, 无法测量出弹目相对距离、相对速度和目标加速度等先进制导律所需要的制导参数问题。尽管有些可以近似看成线性系统, 但是大多数的系统不仅不能用线性微分方程描述, 而且其非线性因素还不能忽略。此外有时为了更加精确地得到滤波结果, 也必须应用反映实际系统的非线性模型, 因此对于非线性系统的滤波是一个必须要解决的问题。

由于一般的非线性系统在理论上难以找到严格的递推滤波公式, 因此目前只能采用近似方法研究, 而线性化是用近似方法来研究非线性滤波问题的重要途径之一, 这就是扩展卡尔曼滤波算法 (EKF) 。但是EKF线性化过程中忽略了二阶以上的分量, 因此在滤波精度上存在着较大误差。为了弥补EKF方法的不足, 人们希望通过对非线性函数的概率密度分布近似, 来代替对非线性函数的近似, 这样就可以利用采样逼近的方法来解决非线性问题。1997年, Juliear S.J.和Uhlman J.K.提出了一种新的非线性滤波方法——Unscented卡尔曼滤波 (UKF) 。UKF不需要对非线性系统的状态方程和观测方程进行线性化, 而是利用Unscented变换 (UT) 方法来近似非线性函数的概率密度分布, 因此UKF方法在计算精度上要高于EKF方法, 并且不需要计算状态转移矩阵的雅可比矩阵, 这使得其应用范围更加广泛。

1 UT

UT是UKF的基础, UT的思想是用固定数量的参数去近似一个高斯分布, 这比近似任意的非线性函数或变换更容易。

UT的具体公式如下:

式中:λ = α 2 ( n + k ) n是一个比例因子。 α 决定x ˉ 周围 χ 点的分布状态, 调节 α 以使高阶项的影响达到最小, 通常选择0 ≤ α ≤1 。对高斯分布的情况, 当状态变量为单变量时, 选择k = 2 ;当状态变量为多变量时, 选择) i是矩阵的第i列 (当P =ATA 时, 取的第i行;当P = AAT的第i列) 。 χ 点的选取就是选取尽可能代表X分布点, 而这些点的分布程度取决于的大小。

y的均值和方差可以通过下述公式获得:

其中Wim和Wip分别对应着计算y的均值和方差的加权系数, 要求

2 UKF滤波算法

设非线性系统的状态方程和观测方程为:

式中:x (k) 为n维状态向量;z (k) 为m维的观测向量; w (k) 为系统噪声, v (k) 为观测噪声, 假设它们是均值为零的高斯白噪声, 且互不相关;φ (⋅) 为n维向量方程, 是x (k) 、w (k) 和k的非线性函数;h (⋅) 为m维向量方程, 是x (k + 1) 、v (k + 1) 和k + 1的非线性函数。具体算法如下:

(1) 设置初值

(2) 当k > 1, 计算2n + 1个 χ点:

(3) 时间更新

(4) 测量更新

扩展卡尔曼滤波是通过对非线性方程进行线性化变化得到线性部分, 经过泰勒展开式可以得出这种方法的精度为一阶水平, 而UKF算法则可以使均值精确到非线性部分泰勒展开式的三阶水平, 方差精确到二阶水平。

3 UKF对目标位置预测的应用

3.1模型

对目标飞机的运动进行建模时, 可以将目标飞机看成一个质点, 由于假设目标作直线运动的论文比较多, 在此本文主要做目标机动的建模仿真。飞机最常见的机动动作是盘旋, 而飞机匀速圆周运动的定常盘旋最具有代表性, 故以此来建模仿真。其速度与过载的关系如下:

式中:VT表示目标飞机速度;θ 表示倾角;ψ 表示偏角; XT, YT, ZT为目标飞机在地面坐标系中的位置。

飞机在遭遇导弹攻击以后, 一般会进行机动飞行来逃避导弹攻击, 下面建立飞机等过载机动模型:

式中:nx表示飞机切向过载;ny, nz表示法向过载在y轴和z轴上的分量。

3.2仿真及结果分析

3.2.1不同过载下目标机动

设目标飞机以0.6 Ma≈200 m/s的速度做nf= 2到nf= 9的盘旋机动, 跟据公式 (12) 计算出匀速状态下飞机各个过载的盘旋半径和周期, 观测噪声的方差矩阵R = [0.01], 滤波误差的协方差初始矩阵为P (0) = [2 500, 0, 0;0, 2 500, 0;0, 0, 0.000 1], 离散时间T = 0.1 s, 分别进行仿真, 仿真结果如图1~ 图7所示。在UKF中 α = 0.5, β = 2。

从图1, 图2中可以看到使用EKF和UKF都能对目标的位置进行预测, 但是当误差稳定时, UKF的误差要小于EKF的误差, 同时从图2中还能看到UKF的收敛速度要高于EKF。

图3~图6分别是飞机在过载为4和过载为7的盘旋条件下进行的弹目仿真。从图上可以看到随着飞机的过载加大, 飞机的机动性增强, 机动半径减小, 机动时间缩短, 跟踪算法由于目标的机动性增大而相应误差增大, 跟踪时间变长, 但是UKF算法在飞机做大机动的前提下仍优于EKF算法。飞机在不同过载下分别被跟踪的数据见表1。

从表1中可以看到, 随着过载的增大, EKF和UKF的跟踪效果都在降低, 主要原因在于目标机动的过载变大, 时间变短, 机动的半径变得较小, 跟踪误差增大。但是两种算法相比较, UKF仍优于EKF算法。

3.2.2测量干扰较大的目标机动

从上述的实验仿真确定了在相同的条件下UKF算法的精度和收敛速度要高于EKF算法。但是在现实中, 导弹的探测信息中存在多种噪声, 如视线角速度的量测误差, 导引头量测误差, 失调角零位的测量误差等等, 这些因素都会引起信号量测值的波动。因此算法的抗干扰能力对于导弹的跟踪效果显得尤为重要。假设目标以0.6 Ma的速度匀速盘旋, 观测噪声方差阵R =[0.01], 陀螺回转中心与位标器质心不重合引起的漂移服从0.01*N (0, 1) , 滤波误差的协方差初始矩阵P (0) = [2 500, 0, 0;0, 2 500, 0;0, 0, 0.000 1], 离散时间T = 0.1 s, 仿真如图7所示。

从图7中可以看到UKF算法的抗干扰能力更强, 更适合在干扰较强的环境中运用。

4结语

本文介绍了UT和UKF的概念, 并将其应用到飞行目标的跟踪, 在观测值为角度的情况下, 对目标的状态进行了较好的估计, 同时在比较UKF和EKF算法中, 反映出UKF能对所有高斯输入向量的非线性函数进行近似, 均值精确到三阶, 方差精确到二阶, 并且不需要计算雅可比矩阵来对非线性函数作近似变换, 能处理非可导的非线性函数, 计算量与EKF相当。理论分析和仿真结果均表明, UKF算法对于目标方位跟踪领域较其他以往的滤波算法更加稳定, 精度更高。

参考文献

[1]JULIER S J, UHLMANN J K.A new extension of the Kalman filter to nonlinear systems[C]//Proceedings of the 11th International Symposium on Aerospace/Defence Sensing, Simulation and Controls.Orlando, Florida:IEEE, 1997:1-12.

[2]ZHOU Hong-ren, KUMAN K S P.A current statistical model and adaptive algorithm for estimating maneuvering targets AIAA[J].Journal of Guidance, Control and Dynamics, 1984, 7 (5) :596-602.

[3]高磊, 严卫生, 任章, 等.改进型自适应推广卡尔曼滤波器的应用[J].火力与指挥控制, 2000, 25 (3) :13-16.

[4]刘铮.UKF算法及改进算法研究[D].长沙:中南大学, 2009.

[5]张红梅, 邓正隆.UKF方法在陆地车辆组合导航中的应用[J].中国惯性技术学报, 2004, 12 (4) :20-23.

[6]管旭军, 芮国胜.基于UKF的单站无源定位算法[J].电光与控制, 2004, 11 (1) :34-36.

[7]吴玲, 卢发兴, 刘忠.UKF算法及其在目标被动跟踪中的应用[J].系统工程与电子技术, 2005, 27 (1) :49-51.

[8]秦永元, 张洪钺, 汪叔华.卡尔曼滤波与组合导航原理[M].西安:西北工业大学出版社, 1998.

联邦卡尔曼滤波篇8

近年来,组合导航技术发展很快,由于组合导航系统充分有效地利用了各子导航系统的导航信息,在降低成本的同时还提高了系统的精度和可靠性,因而在许多领域都得到了广泛的应用。随着我国北斗二号卫星导航系统开始逐步提供服务,国内采用BD-2/SINS组合导航的车载系统的研究开始受到关注[1,2]。但当BD-2/SINS车载组合导航系统工作于恶劣条件时,随时都有可能出现故障或者工作不正常的状况,因此,故障诊断与容错技术是车载BD-2/SINS组合导航系统设计和实现的关键技术。研究能够实时确定BD-2/SINS车载组合导航各子系统信息输出合理性的方法以及保证系统信息输出能够达到可靠性要求的系统结构,已经成为车载BD-2/SINS组合导航系统研究的重要方向。

1 车载BD-2/SINS组合导航的故障诊断算法

车载BD-2/SINS组合导航由SINS和BD-2系统构成。捷联惯性制导系统(SINS)具有完全自主性和全天候工作能力,可连续提供实时的全参数(位置、速度、姿态)导航信息,具有抗干扰能力强以及短时导航精度高等优点,通常作为主导航系统。北斗二号卫星导航系统(BD-2)具有更高的导航定位精度且不随时间发散,但由于系统工作需要依赖外部信息,当卫星信号由于遮挡或者受到干扰时,会出现丧失导航定位能力或者定位误差急剧增大的情况,因此将BD-2作为辅助导航系统。BD-2/SINS组合导航系统的多源冗余信息为提高系统的可靠性与容错性创造了条件,通过合理设计系统结构、算法和组合模式,并采用故障诊断技术与系统重构方法,能够大大提高系统的无缝导航能力。对于组合导航的系统级故障检测,目前应用较多的是基于系统模型的统计检验方法,主要思想就是查验系统输出信号是否在统计意义上服从高斯分布并且具有假设的均值和协方差。按照构成检验统计量方法的不同,可分为两种方法:残差χ2检验和状态χ2检验[3]。

1.1 残差χ2诊断原理

残差χ2检验法通过对卡尔曼滤波的残差所具有的统计特性进行假设检验。

定义第k时刻系统状态的递推值 $\hat{X} (k | k - 1)$

$\hat{X} (k | k - 1) = Φ (k, k - 1) \hat{X} (k - 1) (1)$

k时刻系统量测的预测 $\hat{Ζ} (k | k - 1)$

$\hat{Ζ} (k | k - 1) = Η (k) \hat{X} (k | k - 1) (2)$

则卡尔曼滤波器的残差d(k)可表示成

$d (k) = Ζ (k) - \hat{Ζ} (k | k - 1) (3)$

其协方差矩阵为

S(k)=H(k)P(k|k-1)HT(k)+R(k) (4)

定义距离函数D(k)

D(k)=dT(k)S-1(k)d(k) (5)

根据卡尔曼滤波器残差序列的统计特性,D(k)应服从自由度为m(观测向量z的维数)的χ2分布。若故障发生,残差d(k)将不再是零均值白噪声过程,则D(k)将变大。因此可采用这一思路实现对故障的诊断。取D(k)大于某一门限M的概率为α,即

P{D(k)>M}=α (6)

式(6)中α为允许的虚警概率。

此时检测过程为

${\begin{cases} D (k) > Μ ‚ 系统故障发生 \\ D (k) \leq Μ ‚ 系统 ⌶ 作正常 \end{cases} (7)$

1.2 状态χ2诊断原理

研究表明,残差χ2检验法对于硬故障检测效果比较好,但对于变化比较缓慢的软故障则不容易检测出来,对这类故障经常采用状态χ2检验法进行检测。状态χ2检验法利用由量测值Z(k)经卡尔曼滤波得到的 $\hat{X}$ 1(k)和由先验信息递推计算而得的 $\hat{X}$ 2(k)这两个状态进行估计。其中 $\hat{X}$ 1(k)和测量信息有关,会受到故障影响, $\hat{X}$ 2(k)由先验信息递推计算而得到,与测量信息无关,不受故障影响,利用二者之间的差异即可对故障进行检测和隔离。定义估计误差

${\begin{cases} e_{1} (k) = \hat{X}_{1} (k) - X (k) \\ e_{2} (k) = \hat{X}_{2} (k) - X (k) \end{cases} (8)$

定义

$β (k) = e_{1} (k) - e_{2} (k) = \hat{X}_{1} (k) - \hat{X}_{2} (k) (9)$

当系统发生故障时,由于估计 $\hat{X}$ 2(k)与测量值Z(k)无关,故仍是无偏估计;而估计 $\hat{X}$ 1(k)因受故障影响变成了有偏估计。则

$E β (k) = E \hat{X}_{1} (k) - E \hat{X}_{2} (k) \neq 0 (10)$

对β(k)作二元假设:

H0:无故障Eβ(k)=0;Eβ(k)βT(k)=T(k)

H1:有故障 Eβ(k)=μ;E{[β(k)-μ][βT(k)-μ]}=T(k)

由于β(k)是高斯随机向量,故有以下条件概率密度函数:

$\begin{array}{l} p_{r} (β / Η_{0}) = \frac{1}{\sqrt{2 π} | Τ (Κ) |^{\frac{1}{2}}} \times \\ \exp [- \frac{1}{2} β^{Τ} (k) Τ^{- 1} (k) β (k)] (11) \end{array}$

$\begin{array}{l} p_{r} (β / Η_{1}) = \frac{1}{\sqrt{2 π} | Τ (Κ) |^{\frac{1}{2}}} \times \\ \exp {- \frac{1}{2} [β^{Τ} (k) - μ^{Τ}] Τ^{- 1} (k) [β (k) - μ]} (12) \end{array}$

pr(β/H1)和pr(β/H0)的对数似然比为

$\begin{array}{l} Λ (k) = \ln \frac{p_{r} (β / Η_{1})}{p_{r} (β / Η_{0})} = \frac{1}{2} {β^{Τ} (k) Τ^{- 1} (k) β (k) - \\ [β^{Τ} (k) - μ^{Τ}] Τ^{- 1} (k) [β (k) - μ]} (13) \end{array}$

式(13)中的μ是未知的,可用其极大似然估计 $\hat{μ}$ 代替。求 $\hat{μ}$ 使Λ(k)达到极大,得

$\hat{μ} (k) = β (k) (14)$

将其代入式(13),得到故障检测函数

λ(k)=βT(k)T-1(k)β(k) (15)

由于β(k)是高斯随机向量,故λ(k)服从自由度为n的χ2分布。故障判断准则为:

${\begin{cases} λ (k) > Τ_{D} ‚ 有故障 \\ λ (k) \leq Τ_{D} ‚ 无故障 \end{cases} (16)$

式(16)中TD是预先设置的门限,它决定了故障检测的性能。

1.3 基于卡尔曼滤波器组的故障检测方法

状态χ2检验的基本原理是要求系统具备一个高精度的参考系统,从而通过检验卡尔曼滤波器状态估计与参考系统的一致性来检测故障,因此故障检测的性能很大程度上取决于参考系统的精度。但用状态递推器作为参考系统进行χ2检验时,由于状态递推器并不进行量测更新,误差将使状态递推值越来越偏离真实值,导致系统处于无故障状态时的β(k)也会越来越大,使故障检测的灵敏度降低。现采用一种基于辅助卡尔曼滤波器的方法,为状态递推器引入量测量,提高状态递推值的精度,基本结构如图1所示。

图中,主滤波器M使用系统的全部量测值作为输入,当系统无故障时,得到系统状态的最优估计 $\hat{X}$ M。辅助滤波器A1和A2分别使用系统量测的一个子集作为量测量,其状态估计为 $\hat{X}$ A1和 $\hat{X}$ A2,用来作为故障检测的参考。不管系统发生什么故障,A1和A2其中之一总会受到影响,通过检验主滤波器M和辅助滤波器A1和A2状态估计之间的一致性便可以进行故障检测。由于辅助滤波器引入了量测更新,使参考系统的精度得到保证,因而不会出现状态递推器误差逐渐变大的情况。

2 BD-2/SINS组合导航联邦滤波器方案及仿真分析

2.1 系统联邦滤波器的结构

联邦滤波器是一种具有两级结构的分散化滤波方法,它能够通过各滤波器之间的信息分配而得到全局最优估计,具有结构灵活、计算简单及较好的容错性等特点[4]。基于故障诊断和容错为基本需求的BD-2/SINS组合导航联邦滤波器的基本结构可设计如图2所示。

图2中,公共参考系统为SINS系统,其输出一方面直接给主滤波器,另一方面它给BD-2局部滤波器做为其测量值。BD-2子系统输出给局部滤波器,其量测量为SINS和BD-2位置信息输出之差,其局部估计X1及协方差阵P1送入主滤波器和主滤波器的估计和协方差阵一起进行融合得到全局最优估计 $\hat{X}$ g和协方差阵Pg。而后Pg被放大为β $_{1}^{- 1}$ Pg(β1≤1),反馈到局部滤波器和主滤波器中,并以 $\hat{X}$ g重置局部滤波器和主滤波器的估计值。当系统运行时,首先在局部滤波器中采用残差χ2检验法对测量故障进行检测,然后再进行状态χ2检验,如果残差χ2检验能够检测出故障,则不再进行状态χ2检验。在进行状态χ2检验时,设计一个参考滤波器,它实质上是以SINS先验模型为状态方程的状态递推器。在滤波和故障检测的过程中,参考滤波器只根据SINS的初始状态和初始方差进行时间更新,不进行测量更新。这样的联邦滤波结构能够很容易地实现子系统的出现故障诊断和容错,当发现BD-2子系统出现故障时,主滤波器立即拒绝该子滤波器的输出,即可实现故障的隔离,系统仅采用SINS进行导航定位,达到无缝导航的目的。

2.2 系统联邦滤波器故障诊断和容错的仿真研究

根据前面分析,进行BD-2/SINS组合导航系统及其故障检测算法的仿真研究,其中仿真工具选用Matlab Simulink。

建立各导航子系统的误差仿真模型。把导航信息误差作为系统状态向量,各导航子系统输出的导航信息之差作为观测量。其中,SINS的误差传播方程采用文献[5]所述。将BD-2各种误差因素造成的总的定位误差用一阶马尔可夫过程等效,表示为:

${\begin{cases} \dot{ϕ}_{g} = - \frac{1}{Τ_{g}} ϕ_{g} + w_{g} \\ \dot{λ}_{g} = - \frac{1}{Τ_{g}} λ_{g} + w_{g} \end{cases} (17)$

式(17)中,Tg为相关时间,wg为激励马尔可夫过程的白噪声。

设计卡尔曼滤波模块和故障检测模块,仿真模型如图3所示。

图3中,SINS、BD-2分别为捷联惯导和BD-2接收机仿真模型,输出相应的导航误差信息,车辆行驶的具体参数由m文件载入。其中二个卡尔曼滤波器所使用的状态方程相同,但量测量不同。Kalman filter M使用全部位置和速度量测,而Kalman filter A1仅使用位置作为量测。Fault detection模块对二个卡尔曼滤波器所计算的状态估计和协方差阵进行一致性检验,并输出相应的故障检测结果。

2.3 仿真结果及分析

设定无故障、硬故障(突变)和软故障(慢变)三种不同情况,分别采用残差χ2检验法和状态χ2检验法对三种不同情况的观测量(BD-2的位置信息)进行故障检测。

图4(a)和图4(b)给出了观测量无故障时残差χ2检验法和状态χ2检验法的故障检测函数曲线。图中,无故障时两种检测方法的故障检测函数值一直远远小于1,观测量无故障,可直接作为卡尔曼滤波器的观测量使用。

在600 s处设置硬件故障,具体表现为BD-2观测值突然增大。图5(a)和图5(b)给出了采用残差χ2检验法和状态χ2检验法的故障检测函数曲线。

可以看出,当硬故障出现时,采用残差χ2检验法进行故障检测,函数值λ迅速增大且远远大于1。而采用状态χ2检验法时,函数值λ经过一段时间延迟后其值才变大并超过1,这说明状态χ2检验法对于硬件故障的检测中具有一定的延迟。

在600 s到1 275 s时设置软故障,设定故障幅值变化为130(t-600)m。分别采用残差χ2检验法和状态χ2检验法进行故障诊断,得到的函数曲线如图图6(a)和图6(b)所示。

图中可以看出,当软故障出现时,残差χ2检验法的故障检测函数值有所增加,但一直都保持小于1的状态,可见,残差χ2检验法对软故障的诊断不灵敏。而采用状态χ2检验法进行故障诊断,虽然在故障初始阶段检测函数值小于1,无法诊断故障的出现,但经过一段时间延迟后,故障诊断函数值就明显增大且超过1,进而实现了对故障的诊断。

3 结论

对于采用联邦滤波的车载BD-2/SINS组合导航的故障诊断与容错,残差χ2检验法对故障的响应告警延迟量小,灵敏度明显强于状态χ2检验法,但当系统出现软故障时存在漏检的风险,而状态χ2检验法能有效地检测出系统出现的硬故障和软故障,但其告警延迟量大,实际应用中,将两者加以结合,即可较好的实现对车载BD-2/SINS组合导航的故障诊断和容错处理。

摘要：论述了车载BD-2/SINS组合导航的系统级故障诊断方法,包括残差χ2检验方法和状态χ2检验方法。针对状态χ2检验原理的不足,提出一种基于卡尔曼滤波器组的改进的状态χ2检验方法。文中构建了基于联邦滤波的车载BD-2/SINS组合导航的故障诊断系统结构,将残差χ2诊断和改进的状态χ2诊断方法加以整合,能够实现对BD-2/SINS组合导航硬故障和软故障的诊断,并针对无故障、硬故障和软故障三种不同情况分别进行了仿真研究和分析。

关键词：BD-2/SINS组合导航,故障诊断,卡尔曼滤波,联邦滤波,容错

参考文献

[1]高社生,胡攀,王海维.SINS/SAR/北斗组合导航系统.中国惯性技术学报,2008;(2):178—182

[2]戴邵武,马长里,代海霞.北斗双星/SINS组合导航中的捷联惯导算法研究.计算机与数字工程,2010;38(2):1—3,10

[3]韩瑞,秦红磊,丛丽,等.JIDS/SINS/GPS组合导航系统两级故障检测结构设计.航空电子技术,2008;39(3):43—49

[4] Lee Tae-Gyoo.Centralized Kalman filter with adaptive measurementfusion:its application to a GPS/SDINS integration system with an ad-ditional sensor.International Journal of Control,Automation,andSystems,2003;(4):December,444—452

【联邦卡尔曼滤波】推荐阅读：

卡尔曼滤波法07-29

卡尔曼粒子滤波09-13

自适应卡尔曼滤波06-16

联邦主义06-08

联邦设计08-02