扩展卡尔曼滤波算法

扩展卡尔曼滤波算法（精选9篇）

扩展卡尔曼滤波算法 第1篇

关键词：公路货运周转量,预测,BP算法,扩展卡尔曼滤波

公路货运受多种因素影响,各因素之间的相互作用很难用精确的数学语言来准确表达,因此可以把公路货运周转量的预测看成是一个复杂的非线性系统问题。神经网络理论中的Kolmogorov连续性定理表明任一连续函数都可以精确地由一个3层神经网络来实现,任何一组时间序列都可以视为一个非线性机制确定的输入输出系统,因此该理论从数学上保证了神经网络用于时间序列预测的可行性。

在货运周转量预测的研究中,文献[1]采用BP算法来进行货运周转量的预测。BP算法采用梯度下降法来确定权值调整量,将确定网络权重和阈值分开计算,使得训练速度加快,但是梯度下降法具有计算负担重的缺点,因而在实际应用中受到一定的限制。针对BP算法的缺点,本文参照参考文献[2,3],根据扩展卡尔曼滤波(extended kalman filter, EKF)算法优化原理,建立了可以应用于EKF的神经网络权值调整模型,并以历年全国公路货运周转量为例,分别采用BP训练算法和EKF训练算法对神经网络进行训练,预测结果对比表明了EKF算法的有效性和准确性。

1 扩展卡尔曼滤波(EKF)

扩展卡尔曼滤波是一种最优估计,它采用状态方程和观测方程组成的非线性随机系统的状态空间模型来描述滤波器,是解决非线性滤波问题常用的和有效的一种方法。它采用状态方程的递推性,按线性无偏最小均方差估计准则,采用递推算法对滤波器的状态变量做最佳估计,从而求得滤掉噪声的有用信号的最佳估计。现考虑以下非线性离散时间状态系统:

$\{\begin{array}{l}\mathit{X}{}_k=f(X{}_{k-1},k-1)+\mathit{W}{}_{k-1}\\ \mathit{{\rm Z}}{}_k=h(X{}_k,k)+\mathit{V}{}_k\end{array}(1)$

式中:Xk和Zk分别为系统在k时刻的状态变量和测量变量,函数f()和h()分别为非线性系统的状态方程和测量方程;Wk和Vk分别为过程噪声和测量噪声,并设有以下性质:

$\{\begin{array}{l}E[\mathit{W}{}_k]=0,E[\mathit{W}{}_k\mathit{W}{}_j^{\text{Τ}}]=\mathit{Q}{}_k\mathit{U}{}_{kj}\\ E[\mathit{V}{}_k]=0,E[\mathit{V}{}_k\mathit{V}{}_j^{\text{Τ}}]=\mathit{R}{}_k\mu {}_{kj}\\ E[\mathit{W}{}_k\mathit{V}{}_j^{\text{Τ}}]\end{array}(2)$

式中:Qk和Rk分别为k时刻系统过程噪声协方差矩阵和观测噪声协方差矩阵。由以上假设可得扩展卡尔曼滤波算法递推过程如下:

1) 时间更新(预测)。向前推算状态变量:

$\widehat{X}{}_{k,k-1}=f(\widehat{X}{}_{k-1},k-1)(3)$

向前推算误差协方差:

${\rm P}{}_{k,k-1}=\phi {}_{k,k-1}{\rm P}{}_{k-1}\mathit{\phi }{}_{k,k-1}^{\text{Τ}}+Q{}_{k-1}(4)$

2) 测量更新(校正)。

① 计算卡尔曼增益:

${\rm K}{}_k={\rm P}{}_{k,k-1}\mathit{{\rm H}}{}_k^{\text{Τ}}(\mathit{{\rm H}}{}_k{\rm P}{}_{k,k-1}\mathit{{\rm H}}{}_k^{\text{Τ}}+R{}_k){}^{-1}(5)$

② 由观测变量Zk更新估计:

$\widehat{X}{}_k=\widehat{X}{}_{k,k-1}+{\rm K}{}_k(\mathit{{\rm Z}}{}_k-h(\widehat{X}{}_{k,k-1},k))(6)$

③ 更新误差协方差:

${\rm P}{}_k=({\rm I}-{\rm K}{}_k\mathit{{\rm H}}{}_k){\rm P}{}_{k,k-1}(7)$

式(3)～(7)中:$\widehat{X}$k,k-1、Pk,k-1分别为在k-1阶段对k阶段状态变量和误差协方差的先验估计,$\widehat{X}$k、Pk分别为第k阶段的状态变量和误差协方差的后验估计,φk.k-1、Hk计算如下:

$\{\begin{array}{l}\mathit{\phi }{}_{k.k-1}=\frac{\partial f(\widehat{X}{}_{k-1},k-1)}{\partial X{}_{k-1}}|{}_{X{}_{k-1}=\widehat{X}{}_{k-1}}\\ \mathit{{\rm H}}{}_k=\frac{\partial h(\widehat{X}{}_{k,k-1},k)}{\partial X{}_k}|{}_{X{}_k=\widehat{X}{}_{k,k-1}}\end{array}(8)$

2 利用EKF算法训练前向神经网络

EKF算法用于神经网络训练的思想是将网络权值作为EKF滤波器的状态变量,通过输入输出数据不断调整状态,直到估计误差满足要求。因此,为利用EKF算法对图1所示的三层前向神经网络进行训练,令神经元之间的连接权值为EKF的状态变量,神经网络的输出为EKF的观测变量。如图1所示网络的滤波器状态变量可表示为:x=[w${}_{11}^{1,2}$w${}_{mn}^{1,2}$w${}_{11}^{2,3}$w${}_{nd}^{2,3}$]T。

根据式 (1),EKF的状态方程和观测方程可以表示为:

$\{\begin{array}{l}x(k+1)=x(k)\\ y{}^d(k)=y(k)+\epsilon (k)=h(x(k),u(k))+\epsilon (k)\end{array}(9)$

式中:y(k)和yd(k)为网络实际输出和期望输出;h()为网络输出、输入和权值间的非线性函数关系;u(k)为网络输入;ε(k)为随机白噪声。将期望输出展开为泰勒级数并略去二次以上项得:

$y{}^d(k)=h(x,(k),u(k))+\mathit{{\rm H}}(k)(x(k)-\widehat{x}(k))(10)$

式中:H(k)为前述h对x的偏导雅可比矩阵,因此可得EKF训练算法对应的计算步骤如下:

1) 随机设置网络初始权值$\widehat{x}$ (1)(滤波器的初始状态变量)以及误差协方差P(1)。

2) 对样本k=1到N,重复如下步骤直至达到设定的精度要求:

根据输入u(k)及权值$\widehat{x}$(k)计算输出h($\widehat{x}$(k),u(k)),并根据期望输出进行状态估计及更新:

$\{\begin{array}{l}{\rm K}(k+1)={\rm P}(k)\mathit{{\rm H}}{}^{\text{Τ}}(k+1)\\ [\mathit{{\rm H}}(k+1){\rm P}(k)\mathit{{\rm H}}{}^{\text{Τ}}(k+1)+\\ {\rm P}(k+1)]{}^{-1},\\ \widehat{x}(k+1)=\widehat{x}(k)+{\rm K}(k+1)\\ \left[y{}^d(k)-h(\widehat{x}(k),u(k))\right],\\ {\rm P}(k+1)=\left[{\rm I}-{\rm K}(k+1)\mathit{{\rm H}}{}^{\text{Τ}}(k+1)\right]{\rm P}(k)\end{array}(11)$

3 实例分析

为验证以上算法的有效性,本文以全国历年公路货运周转量为例,考虑影响公路货运周转量的因素有:GDP、公路运输线路长度和民用载货汽车拥有量3个,并以各自的历史数据作为一个3层前向人工神经网络的输入,公路货运周转量作为输出,根据经验取隐层节点数为10,网络结构如图1所示,其输入输出为u1,u2,u3和y,状态变量:

x=[w${}_{11}^{1,2}$w${}_{12}^{1,2}$w${}_{110}^{1,2}$w${}_{31}^{1,2}$w${}_{32}^{1,2}$w${}_{310}^{1,2}$w${}_1^{2,3}$w${}_2^{2,3}$

w${}_{10}^{2,3}$]T ,共有Nk=310+101=40,即40个权值。

取该网络隐层的传递函数为Sigmoid函数,即y=1/(1+exp(-x)),输入层和输出层的传递函数取比例系数为1的线性函数,各层神经元的阈值取0,设uj为输入层第j个神经元的输入,Oi为隐层第i个神经元的输出,对Sigmoid函数求导可得∂y/∂x=y(1-y),所以根据式( 8)可得观测方程中的传递矩阵为:

$\mathit{{\rm H}}=\{\begin{array}{l}{\rm O}{}_i(1-{\rm O}{}_i)w{}_i^{2,3}u{}_j,\text{当}\text{地}\text{时}\text{间}w\text{为}\text{输}\text{入}\text{层}\\ \text{与}\text{隐}\text{层}\text{连}\text{接}\text{权}\text{时}‚\\ i=1,2,\cdots ,10;\\ j=1,2,3\\ {\rm O}{}_i\text{当}w\text{为}\text{隐}\text{层}\text{与}\text{输}\text{出}\text{层}\\ \text{连}\text{接}\text{权}\text{时}‚\\ i=1,2,\cdots ,10\end{array}(12)$

下面便可利用EKF算法对网络进行训练。在滤波器递推过程中,初始状态X(1),P(1)难以准确掌握,但由于EKF在递推过程中不断用新的信息对状态进行修正,所以当滤波时间充分长时,初值X(1)对X(k+1)的影响将衰减至近似于零,P(1)对P(k+1)的影响也将衰减至近似于零,因此,滤波的初始条件可以近似确定。根据以上分析,利用表1所示原始数据,通过Matlab7.0编程,对神经网络分别利用BP算法和EKF算法进行训练,得到图2所示的2种训练算法收敛速度对比图和表2所示的2种训练算法的预测结果。

由计算结果可见,EKF算法是一种有效的算法,它能够在较短的时间内得到较理想的预测结果,提高了算法的数值稳定性,不仅节省了训练时间,也减少了迭代次数、提高了网络输出精度,从而可以应用于公路货运周转量的预测。

4 结束语

本文提出了以EKF算法作为人工神经网络的训练方法并将该算法应用于全国公路货运周转量预测的研究中,该算法预测的结果与BP算法的预测结果对比表明:EKF算法应用于神经网络训练具有较好的运算效果,提高了神经网络的运算速度和精度。另外,文中的算例仅以公路货运周转量预测为例,对于该训练算法在各种交通方式中的其它应用将有待进一步的研究。

参考文献

[1]陈富贵,李远富.基于BP神经网络的铁路货运量预测[J].铁道货运,2005(9):15-17

[2]Yang Huizhong,Li Jiang,Ding Feng.A neuralnetwork learning algorithm of chemical processmodeling based on the extended Kalman filter[J].Neurocomputing,2007,70(4):625-632

[3]缑娜,王睿,付莹.基于自适应卡尔曼滤波器的神经网络算法[J].弹箭与制导学报,2006,26(2):272-274

[4]Dan S.Training radial basis neural networks withthe extended Kalman filter[J].Neurocomputing,2002,48(1):455-475

[5]寇毅.组合预测模型在公路货运周转量预测中的运用[J].中国水运,2007(11):216-217

扩展卡尔曼滤波算法 第2篇

在系统机动性不强的情况下,传统的平台内阻尼算法将系统本身的速度信息通过阻尼网络加到系统中,达到提高姿态角精度的目的.将这种平台内阻尼的思想引入到捷联惯性航姿系统中,在系统加速度较小的情况下,利用加速度计的输出估计系统姿态角,通过卡尔曼滤波的形式补偿系统姿态误差.由于加速度的.大小直接影响滤波器精度,本文设计了模糊自适应卡尔曼滤波算法,根据三轴加速度计的输出调整内阻尼量测误差方差阵,从而避免了滤波器的发散.仿真和实验验证,内阻尼算法可明显抑制舒勒周期振荡和傅科周期振荡,避免了系统姿态漂移,有效提高了捷联惯性航姿系统的精度.

作 者：杜亚玲 刘建业 刘瑞华 祝燕华 DU Ya-ling LIU Jian-ye LIU Rui-hua ZHU Yan-hua 作者单位：杜亚玲,DU Ya-ling(宇航智能控制技术国防科技重点实验室,北京,100854)

刘建业,祝燕华,LIU Jian-ye,ZHU Yan-hua(南京航空航天大学自动化学院导航研究中心,南京,210016)

刘瑞华,LIU Rui-hua(中国民航学院天津市智能信号与图像处理重点实验室,天津,300300)

扩展卡尔曼滤波算法 第3篇

关键词降阶线性卡尔曼；永磁同步电机；无速度传感器

中图分类号TM3文献标识码A文章编号1673-9671-(2011)042-0132-02

卡尔曼滤波算法是由美国学者Rudolf E.Kalman提出的一种最小方差意义上的最优预测估计方法，它提供了直接处理随机噪声干扰的解决方案，将参数误差看作噪声以及把预估计量作为空间状态变量，用递推法将系统及测量随机噪声滤掉，得到准确的空间状态值。

1扩展卡尔曼滤波器

卡尔曼滤波算法是一种线性最优递推滤波算法，能够在系统的状态方程和测量方程中具有噪声时，对系统状态变量的最小方差估计。

在α-β坐标系下PMSM数学模型可以表示为：

（1）

其中：

，（2）

上述数学模型中，系统状态变量x为[iα，iβ，ω，θ]T，输入变量u为[uα，uβ]T，输出变量为[iα，iβ]T。uα，uβ分别为定子α轴、β轴电压，iα，iβ分别为定子α轴、β轴电流，Rs为定子电阻，Ls为同步电感，ψf为转子磁通，ω为转子机械角速度，θ为转子机械位置，w（t）系统过程噪声，

v（t）为系统测量噪声，w（t）与v（t）协方差阵分别为Q、R。

取采样时间为T，将（1）离散化，得到离散化的系统状态方程：

（3）

基于离散化的PMSM状态方程，扩展卡尔曼滤波实现算法如下：

①状态预估：

（4）

（5）

2降阶线性卡尔曼滤波器

EKF估计算法的不足就在于系统的协方差矩阵和增益矩阵必须在每个采样时刻进行更新。如果系统方程能被完全线性化，系统方程的协方差序列在每一个采样时刻是时不变的，那么这种Kalman滤波器称为线性卡尔曼滤波器（LKF）。

如果LKF应用于PMSM无速度传感器控制，首先PSMM电机模型必须线性化，其次系统协方差序列是时不变的。

基于LKF理论的PMSM模型为：

（6）

其中：

，

， （7）

选取系统状态为[θ，ωr，ρ']，θ是转子位置，ωr是转子角速度，ρ是白噪声序列，ρ'是白噪声序列的导数。∑(k)为过程噪声矩阵，Δ(k)为测量噪声矩阵，λ为可调参数，R1为∑(k)的协方差矩阵，R2为Δ(k)协方差矩阵。

如果输出信号y（k）相位相差π/2，则系统的协方差pk会成为一个时不变Riccati差别矩阵，可以通过MATLAB-DLQE计算得到。这使得系统Kalman滤波器增益矩阵可以通过下式计算得到：

（8）

本文选取定子侧电压信号作为输出变量，如式（9）所示，k为信号增益。

（9）

综上所述，应用于PMSM无速度传感器控制的RLKF方程如式（10）所示。

与EKF方程（1）-（5）相比，LKF方程结构大大简化，在运算过程中占用更少的存储空间，更加易于数字化实现。

（10）

3仿真分析

本文的仿真和建模是在Matlab7.0的simulink环境下完成的，图1所示为基于EKF算法的PMSM无速度传感器控制系统的结构图，基于RLKF算法的PMSM无速度传感器控制系统与图1相同，kalman观测器的输入为ualpha、ubeta输出为ωr、θ。

永磁同步电机的运行参数如表1所示。

图2所示为基于EKF（a）和RLKF（b）的电机转子估计转速与实际转速之间的误差。

图2结果显示，两种转速检测方法在低速情况下都存在较大误差，这是由于电机参数变化与低电压信号测量的不确定性引起的。

4结语

本文介绍Kalman滤波算法的基础上，提出了一种新的基于降阶线性卡尔曼滤波算法的永磁同步电机无速度传感器控制方法，并与通过MATLAB仿真，与传统的扩展卡尔曼滤波算法进行了比较，实现结果表明，新算法不仅延续了传统扩展卡尔曼算的优势，而且算法更加简单，减轻了繁重的参数调节任务，易于数字化实现，为永磁同步电机的无速度传感器控制提供了一种新的控制方式。

参考文献

[1]王成元,夏加宽,杨俊友,等.电机现代控制技术[M].北京:机械工业出版社,2006,

[2]江俊,沈艳霞,纪志成.基于EKF的永磁同步电机转子位置和速度估计[J].系统仿真学报,2005,17(7):1704-1707.

作者简介

刘祖全（1983—），男，助理工程师，获硕士学位，2009年毕业于山东大学控制科学与工程学院电力电子与电力传动专业，现从事核电自动化专业设计工作。

扩展卡尔曼滤波算法 第4篇

卡尔曼滤波是20世纪60年代发展起来的一种现代滤波方法, 它的一个重要作用在于系统的状态估计。当噪声是正态分布时, 这种滤波给出了状态的最小方差估计, 当不是正态分布情况时, 这种滤波给出了状态的线性最小方差估计。

1 滤波算法的发展

约200年前, 高斯 (Guass) 提出了参数估计最小二乘法, 它利用一系列观测数据Z (0) , Z (1) , ...Z (k) , 使它与估计值之间的误差平方和最小, 最小二乘法使用简单, 在一些简单的估计问题中, 仍得到广泛应用[2]。1940年左右, 维纳和柯尔莫格洛夫首次用统计方法研究随机系统, 提出了基于最小均方误差准则下的最优线性滤波一维纳滤波。维纳滤波理论奠定了用统计方法研究随机控制问题的基础。1960年左右, 卡尔曼 (Kalman) 和布西 (Bucy) 在维纳滤波的基础上提出了最优线性递推滤波卡尔曼滤波。卡尔曼滤波实际上是一种数据处理的递推算法, 它克服了维纳滤波需要整段数据的缺点, 根据前一时刻的估计和现在的数据可以由递推方程得到新时刻的估计, 使存储量和计算量大大减少。它在线性系统中的应用已经很成熟, 对于非线性系统可采用近似方法求解。由于卡尔曼滤波的上述优点, 它在航天技术、导航、通信和经济管理等方面得到了广泛应用[3]。

2 扩展卡尔曼滤波器

2.1 两相静止坐标系下的感应电机的数学状态空间模型

在两相静止α-β坐标系下, 感应电机的4阶数学模型为:

感应电机的转速机械方程为:

φ (X) 和 (4) 式中的系数定义为式中, K为扭转弹性力矩系数;Tγ为转子时间常数;np为极对数, (2) 中Lm为各相间的互感;Ω为转速;is为定子电流, is=[isα, isβ]T;ψγ为转子磁链, ψγ=[ψαγ, ψβγ]T;σ为漏感系数;Ls为定子自感, (4) 中Lm为定子自感;J为转动惯量;Lγ为转子自感;TL为负载力矩;Rs为定子电阻;Rγ为转子电阻。

2.2 扩展卡尔曼滤波器设计模型描述

模型描述为在两相静止α-β坐标系下, 为了利用卡尔曼滤波算法对电机转子转速进行预测, 需加入转子转速ωγ作为新增的状态变量。在采样周期很短的条件下, 可认为ωγ的变化为零, 即ωγ=0。这样两相静止坐标系下的状态空间模型扩为五阶, 选择定子电流iαs, iβs, 转子磁链ψαγ, ψβγ和转子转速ωγ为状态变量, 感应电机状态空间模型表示[4]:

其中:

式中, Rs为定子电阻;Lm为互感;Ls为定子绕组自感;Lγ为转子绕组自感;ωγ为电机角速度;σ为漏磁系数,

可见α-β两相静止坐标系下, 新的感应电机状态空间模型为一个输入输出由定子电压、电流组成, 系统矩阵由定子和转子参数构成的五阶非线性状态方程。

3 仿真试验

本文用RMxprt软件通过输入相关的电机参数来建立电机结构模型[5], 但是RMxprt软件只能建立单相和三相感应电机的模型, 无法建立多相感应电机的模型。因此, 只能先利用RMxprt建立三相感应电机的模型, 然后再通过改变电机绕组的分布方式, 把三相60°相带感应电机改为六相30°相带的感应电机。图1给出了建立的六相电机的模型图。该电机定子绕组的材料为铜, 转子的导条材料为硅青铜, 定子和转子的冲片的材料为D23, 电机定子的槽数是36, 电机的每极每相槽数为3。图2为在RMxprt软件建立电机模型的基础上, 利用Maxwell软件建立的六相感应电机网格剖分图。

用Maxwell软件进行网格剖分可以大大降低分析人员的工作量, 但是为了进一步提高软件的分析精度, 还需要根据不同的电机部件适当地增加网格剖分密度, 电机模型的网格剖分的各部件的节点数如表1所示。由于电机的气隙、定子齿轭和转子齿轭是电机磁路的组成部分, 因此适当增加节点数目可以得到较高的计算精度。通过求解得到的六相感应电机的计算结果如表2所示。

感应电机A相电流的波形如图3所示, 可见将时间的步长定为0.000 2 s, 计算0.4 s, 仿真得到的电机转矩已经稳定。此时可以得到用扩展卡尔曼滤波器观测到的反电动势波形如图4所示, 由图可以看出电机的反电动势的幅值略低于300 V, 即低于电机额定电压311 V (幅值) 。图5为用扩展卡尔曼滤波器法观测到的六相电机输出转矩特性曲线图, 从图中可以看出在电机进入稳态后, 电机的转矩大小稳定在35 Nm左右。

图6为采用卡尔曼滤波器法观测到的六相电机输出转矩特性曲线图。扩展卡尔曼滤波器对电机低速时的观测效果更加令人满意。

4 结语

本文将卡尔曼滤波用于估计感应电机非线性系统的状态, 设计了扩展卡尔曼滤波器。扩展卡尔曼滤波利用泰勒展开式截断的方式将非线性系统线性化, 线性化以后在形式上同卡尔曼滤波无太多差别。通过对建立的六相感应电机的模型进行电机电磁场的有限元分析, 证明了采用扩展卡尔曼滤波器法能拥有更好的动态跟踪性能。一般来说, 扩展卡尔曼滤波不是最优的, 实际上可以把它作为一种限制复杂性的滤波器。它被限定成具有与线性滤波器类似的结构形式。由于使用了线性逼近, 这种滤波器有可能发散, 在使用过程中尤其要注意这点。

摘要：提出了一种扩展卡尔曼滤波器的设计, 将卡尔曼滤波算法推广至非线性感应电机系统的参数辨识。扩展卡尔曼滤波器采用工程实际中普遍采用的泰勒展开式截断的方法, 对非线性方程进行线性化处理。在仿真中将这种算法应用于六相感应电机模型中, 仿真结果证明扩展卡尔曼滤波器获得了很好的动态跟踪性能以及低速状态下良好的观测效果。

关键词：感应电机,卡尔曼滤波器,扩展卡尔曼滤波器

参考文献

[1]Kailath T.An Innovations Approach to Least-Squares Estimation:Part I.Linear Filtering in Additive White Noise[J].IEEE Trans.Autom.Contr, 1968 (13) :646～655

[2]Kailath T.The Innovations Approach to Detection and Estimaion Theory[J].Proc.IEEE, 1970 (58) :680～695

[3]张贤达.现代信号处理[M].清华大学出版社, 2002

[4]Bose, B.K Power Electronics and AC Drives[M].Prentice-Hall, 1986

扩展卡尔曼滤波算法 第5篇

1 系统描述

定义坐标系如下:

(1) 地理坐标系n系:xn, yn, zn依次指向东北天方向。

(2) 载体坐标系b系:xb, yb, zb依次指向载体右前上方向。

得到姿态矩阵:

分别为运动载体的航向角, 俯仰角和横滚角。其中Cnb与基本旋转顺序有关, 若顺序不同则所得结果不同。

2 欧拉角姿态解算方法

设载体坐标系b系相对地理坐标系n系的角速度为wnb, 再次简记为w, 由wnb在载体坐标系的分量可得:

3 基于扩展卡尔曼滤波器的四元数姿态解算算法

基于卡尔曼滤波器的四元数姿态解算算法是采用扩展卡尔曼方法, 将四元数[q0q1q2q3]T作为状态量, 将加速度计和磁力计的测量作为观测量, 以陀螺测得数据来估计姿态四元数, 以加速度计和磁力计测得数据来修正四元数估计量。其中, 加速度计信息修正侧重水平方向, 磁力计信息修正侧重航向。以下给出该法具体推导。

以四元数与陀螺输出关系建立状态方程, 由四元数与角速率关系式进一步可得到状态方程:

将上式离散化后, 得扩展卡尔曼一步预测方程:

其中t为采样时间。

4 仿真及分析

4.1 仿真条件设置

仿真采用MPU6050惯性测量芯片输出的加速度和角速度做为采样值, 考虑静止状态MPU的角速度和加速度输出。采样间隔为10ms, 仿真时间为1h。得到陀螺输出数据和加速度计输出数据后解算出俯仰角与横滚角。

4.2 仿真结果

图1为解算出的俯仰角与横滚角。分析可知采用基于扩展卡尔曼滤波器的四元数车辆姿态解算算法, 当载体静止或匀速时, 可以将姿态稳定在1°以内, 同时, 由于汽车行驶时, 加速度变化不是非常明显。故该方法适用于对汽车姿态检测。实际系统验证也证明了该法用于汽车姿态检测的可行性。

5 结语

本文采用了基于扩展卡尔曼滤波器的四元数姿态解算算法, 在静止或匀速时能够将解算精度控制在1°以内。可以将该算法应用于行驶车辆获得其姿态信息从而达到对车辆的实时监控。

摘要：测量车辆姿态可以实时跟踪车辆行驶信息, 并通过远程发送车辆姿态信息的方式, 实现对车辆实时的监控。结合车辆姿态测量只考虑俯仰角和横滚角的特殊性, 通过MPU6050芯片采集角速率和加速度数据, 采用基于扩展卡尔曼滤波器的四元数车辆姿态解算算法, 实现对车辆姿态信息的获取。通过采集数据进行MATLAB仿真验证, 结果表明该算法能提高姿态测量精度和系统稳定性。

关键词：扩展卡尔曼滤波器,姿态解算,算法研究

参考文献

[1]张金.基于ARM的车辆姿态测量系统设计[D].北京交通大学, 2008.5.

单雷达航迹滤波与卡尔曼滤波算法 第6篇

单雷达数据处理是空管自动化系统中的重要功能, 它接收原始的雷达数据进行点迹跟踪、航迹处理, 然后把得到的航迹信息进行分发, 供显示终端进行单路雷达航迹显示, 以及多雷达数据处理模块进行航迹融合。正确实现目标/航迹的关联、精确地对航迹进行滤波估计是提高融合效果的关键。本文介绍了单雷达航迹状态滤波原理, 并分析了目前广泛用于航迹滤波的卡尔曼滤波算法。

1航迹状态滤波

在实际的单雷达数据跟踪过程中, 数据相关与状态估计几乎是同时进行的。准确的预测航迹状态是数据相关的前提, 正确的航迹滤波是实现平滑、稳定航迹的关键, 它是单雷达数据处理中的难点。

航迹预测是在本次航迹滤波数值的基础上按照目标运动模型来估计目标未来的状态。滤波 (后验估计值) 用来估计目标当前的运动参数, 如位置、速率和加速度等。它是由本次观测值和预测估值 (先验估计值) 合并处理, 形成的新的航迹状态参数。对那些没有接收到观测的预测航迹来说, 上一次的预测估计就成了滤波估计。要预测航迹未来状态, 就要定义数学模型来描述与预测有关的物理现象, 把下一时刻的目标状态表示为当前时刻目标状态的函数, 当前的目标状态包括目标位置、目标速度、加速度, 其中加速度可以看作是具有随机特性的扰动输入, 目标运动方程如下:

$\mathit{x}{}_{k+1}^{-}=\mathit{A}\mathit{x}{}_k+\mathit{G}\mathit{a}{}_k$

式中:xk是第k次雷达扫描时目标的滤波状态矢量;x-k+1表示根据第k次滤波状态矢量计算出的第k+1次的状态预测;A是转移矩阵;G是噪声增益矩阵;ak是目标的加速度。

一般, 常用的航迹滤波和预测算法有卡尔曼滤波和常系数滤波算法。卡尔曼滤波是根据最小均方误差准则建立起来的估计方法, 它基于包含嗓声的观测值对未知系统状态或参数进行估计, 观测噪声和扰动输入都看作符合高斯过程的白噪声。常用的常系数滤波算法有α-β滤波算法, 它是简化的卡尔曼滤波算法, 其特点为滤波系数是常数而不是动态计算出来的, 因此算法简单运行速度快。卡尔曼滤波有高的准确性, 而常系数滤波算法则有计算上的优势。

这两种算法都可以采用递归调用。过去接收到的数据包含在当前的估计计算中, 因此所有的历史数据都被使用。当有新的观测到来时, 产生“新量”即观测和预测的差, 这种“新量”经过最佳加权去更新目标状态估计, 然后利用动态系统模型和观测模型产生一个最佳的预测, 即在下一次观测到来前作出预测。图1给出了滤波和预测的过程。

图中的增益是用来提供滤波的状态估计和获取下一次扫描时间的预测估计。这些增益可能是固定的, 如常系数滤波算法, 也可能是动态计算出来的值, 如卡尔曼滤波算法。

2卡尔曼滤波算法

卡尔曼滤波算法是根据最小均方误差准则建立起来的估计方法。它用线性递推的方法, 对多个测量数

据和多个信号参数进行处理, 给出无偏差的最小均方误差估计。

卡尔曼滤波器用于估计离散时间过程的状态变量x∈Rn。这个离散时间过程由以下离散随机差分方程描述:

$\mathit{x}{}_k=\mathit{A}\mathit{x}{}_{k-1}+\mathit{B}\mathit{u}{}_{k-1}+\mathit{w}{}_{k-1}(1)$

nn阶增益矩阵A将上一时刻k-1的状态线性映射到当前时刻k的状态。实际中A可能随时间变化, 但在这儿假设为常数。nl阶矩阵B代表可选的控制输入u∈Rl的增益。

定义观测变量z∈Rn, 得到量测方程:

$\mathit{z}{}_k=\mathit{{\rm H}}\mathit{x}{}_k+\mathit{v}{}_k(2)$

式中:nm阶矩阵H表示状态变量xk对测量变量zk的增益。实际中H可能随时间变化, 但在这儿假设为常数。

随机信号wk和vk分别表示过程噪声和观测噪声。假设它们为相互独立, 正态分布的白色噪声, 均值为0, 方差分别为Q和R。

定义先验估计误差的协方差P-k和后验估计误差的协方差Pk:

$\begin{array}{l}\mathit{{\rm P}}{}_k^{-}=E\lfloor (\mathit{x}{}_k-\widehat{\mathit{x}}{}_k^{-})(\mathit{x}{}_k-\widehat{\mathit{x}}{}_k^{-}){}^{{\rm T}}\rfloor (3)\\ \mathit{{\rm P}}{}_k=E\lfloor (\mathit{x}{}_k-\widehat{\mathit{x}}{}_k)(\mathit{x}{}_k-\widehat{\mathit{x}}{}_k){}^{{\rm T}}\rfloor (4)\end{array}$

卡尔曼滤波器中, 先验估计$\widehat{\mathit{x}}{}_k^{-}$和加权的测量变量zk及其预测$\mathit{{\rm H}}\widehat{\mathit{x}}{}_k^{-}$之差的线性组合构成了后验状态估计$\widehat{\mathit{x}}$k:

$\widehat{\mathit{x}}{}_k=\widehat{\mathit{x}}{}_k^{-}+\mathit{{\rm K}}(\mathit{z}{}_k-\mathit{{\rm H}}\widehat{\mathit{x}}{}_k^{-})(5)$

式 (5) 中测量变量及其预测之差$(\mathit{z}{}_k-\mathit{{\rm H}}\widehat{\mathit{x}}{}_k^{-})$称为测量过程的新量或残差, 它的大小反映预测值的准确度, 它的方差反映了这种估计的偏离程度。残差为0表明二者完全吻合。nm阶矩阵K叫做残余的卡尔曼增益或混合因数, 作用是使后验估计误差协方差Pk最小, 即使得后验估计值 (滤波值) $\widehat{\mathit{x}}$k最接近真实值xk。

卡尔曼增益K的计算公式为:

$\begin{array}{l}\mathit{{\rm K}}{}_k=\mathit{{\rm P}}{}_k^{-}\mathit{{\rm H}}{}^{{\rm T}}(\mathit{{\rm H}}\mathit{{\rm P}}{}_k^{-}\mathit{{\rm H}}{}^{{\rm T}}+\mathit{R}){}^{-1}=\frac{\mathit{{\rm P}}{}_k^{-}\mathit{{\rm H}}{}^{{\rm T}}}{\mathit{{\rm H}}\mathit{{\rm P}}{}_k^{-}\mathit{{\rm H}}{}^{{\rm T}}+\mathit{R}}\\ (6)\end{array}$

由式 (6) 可以看出, 观测噪声协方差R越小, 新量的增益K越大;先验估计误差协方差P-k越小, 新量的增益K越小。这可以这样理解:随着测量噪声协方差R趋于零, 测量变量zk的权重越来越大, 而zk的预测$\mathit{{\rm H}}\widehat{\mathit{x}}{}_k^{-}$的权重越来越小;随着先验估计误差协方差P-k趋于零, 测量变量zk的权重越来越小, 而zk的预测$\mathit{{\rm H}}\widehat{\mathit{x}}{}_{\mathit{k}}^{-}$的权重越来越大。

卡尔曼滤波器用反馈控制的方法估计过程状态:滤波器估计过程某一时刻的状态, 然后以含噪声的测量变量的方式获得反馈。因此卡尔曼滤波器可分为两个部分:时间更新方程和测量更新方程。其计算过程的工作原理如图2所示。

时间更新方程负责及时向前推算当前状态变量和误差协方差估计的值, 以便为下一个时间状态构造先验估计。测量更新方程负责反馈, 也就是说, 它将先验估计和新的测量变量结合以构造改进的后验估计。

3结束语

卡尔曼滤波算法既适用于暂态过程也适用于稳态过程, 从跟踪精度上来说它是最好的;从计算时间来说, 它的计算时间虽然较其它滤波方法都长, 但也算是在接受范围之内的一种比较简便的方法。因此在目前空管自动化系统的雷达处理模块中广泛应用, 起到关键的作用。

摘要：单雷达数据处理在民航空管自动化系统中起到基础性的作用, 其航迹质量直接关系着自动化系统的准确性和可靠性。航迹状态滤波是单雷达数据处理中的难点, 是实现平滑、稳定航迹的关键。本文以航迹状态滤波算法为切入点, 介绍了广泛应用于自动化系统中的卡尔曼滤波算法。

关键词：单雷达,航迹状态滤波,卡尔曼滤波算法

参考文献

[1]WELCH G, BISHOP G.An introduction to Kalman filtering[EB/OL].UNC-Chapel Hill, 2006-07-24.http://www.cs.unc.edu/～welch/kalman/.

[2]何友, 王国宏, 等.多传感器信息融合及应用[M].北京:电子工业出版社, 2007:19-25.

基于卡尔曼粒子滤波的目标跟踪算法 第7篇

根据跟踪目标的不同,目标跟踪算法通常可以分为3类[2],即基于偏微分方程的目标跟踪算法、基于Mean Shift的目标跟踪算法和基于滤波理论的目标跟踪方法。文中主要对基于滤波理论的目标跟踪方法进行讨论,这种算法实质是将序列图像的目标跟踪问题建模为目标状态的概率密度函数在时间序列上的传播,通常包括先验密度、观测密度以及后验密度3大要素。

1 基于滤波理论的3种算法描述

1.1 Kalman滤波算法描述

Kalman滤波器[3]是由卡尔曼(R.E.Kalman)于1960年提出的一种递推估计器。Kalman滤波是一种在时域内采用递归滤波对系统状态进行最小均方误差估计的方法,根据系统以前的状态对下一个状态作最优估计,预测时具有无偏、稳定和最优的特点[4]。

Kalman滤波使用要求比较严格,必须满足3个条件[5]:动态系统是线性的;观测是状态的线性组合;噪声是附加的,而且服从正态分布。当这些要求有几个不满足时,Kalman滤波尽管仍然可以使用,但却不能显示其优越性了。

1.2 粒子滤波算法描述

粒子滤波[6]的思想基于蒙特卡洛方法(Monte Carlo Methods),它是利用粒子集来表示概率,也就是通过随机抽取的加权粒子来代替状态的后验概率分布,它可以用于任何形式的状态空间模型上。当随机抽取的粒子数量很大时,结果也就无限接近于状态变量实际的后验密度。

在粒子滤波中,针对许多实际情况,无法直接从后验密度中采样粒子。因此,通常采用与后验密度相近的某一提议密度π(xk|xk-1,y1∶k),从该提议密度中采样新的粒子。在每一离散时刻,新的粒子从提议密度中进行采样,并且计算权重以补偿提议密度和后验密度的差异。将权重归一化之后,后验密度由一组粒子进行近似。

1.3 Kalman粒子滤波算法

由于标准粒子滤波算法选择先验概率密度作为重要密度函数,这种选取方法只有在测量精度要求不高的场合才能取得较好的结果。重要性权重的方差随着时间而随机递增,使得粒子的权重主要集中在少数粒子上,甚至在经过多步递归后,可能只有一个粒子有非零权值,而其他粒子的值很小,几乎可以忽略不计,从而使得大量的计算工作都被浪费在用来更新那些对p(xk|z1∶k)的估计不起任何作用的粒子上,结果粒子集无法表达真实的后验概率分布,这就是粒子滤波算法的退化问题[7]。粒子退化是粒子滤波算法中不可避免的现象,解决方法是选取好的重要性密度函数和粒子重采样。

基于以上原因,文中提出基于卡尔曼粒子滤波器的算法,卡尔曼滤波器选取高斯分布作为提议分布,采用最优提议分布线性化的方法可以得到较优的提议分布。对于每一个粒子,一个单独的卡尔曼滤波器用来产生和传播高斯提议分布。用卡尔曼算法可以得到一个更好的重要密度函数,能够使先验分布朝着高似然度区域移动,在一定程度上避免了粒子退化问题,提高了粒子滤波算法的估计精度。

2 卡尔曼粒子滤波在跟踪中的应用

基于卡尔曼粒子滤波器的跟踪算法流程图如图1所示。

(1)初始化。前文中已经提到,粒子滤波器的基本思想是蒙特罗卡模拟,其中状态的后验密度由一组有权重的粒子$\{(\tilde{x}{}_{k-1}^{(i)},W{}_{k-1}^{(i)}),i=1,\cdots ,{\rm N}\}$来近似。从先验密度p(x${}_0^{(i)}$)中采样,得到初始(令k=1)带权重的粒子$\{(\tilde{x}{}_0^{(i)},{\rm P}{}_0^{(i)}),1/{\rm N},i=1,\cdots ,{\rm N}\}$,其中k表示第k帧,$\tilde{x}{}_{k-1}^{(i)}$为对系统状态x${}_{k-1}^i$的估计,W${}_{k-1}^{(i)}$为相应的权重,p${}_0^{(i)}$是第i个粒子的协方差矩阵。

(2)采样阶段。对于i=1,…,N有

1)卡尔曼预测。

$\stackrel{⌢}{x}{}_{k/k-1}^{(i)}=\mathit{A}\tilde{x}{}_{k-1}^{(i)}+(\mathit{{\rm I}}{}_n-\mathit{A})\overline{x}$ (1)

$\stackrel{⌢}{\mathit{{\rm P}}}{}_k^{(i)}=\mathit{A}\mathit{{\rm P}}{}_{k-1}^{(i)}\mathit{A}{}^{\text{Τ}}+\mathit{Q}$ (2)

其中,A是n×n维状态转移矩阵,由$\tilde{x}{}_{k+1/k}=\mathit{A}\tilde{x}{}_{k/k}$得到;P为协方差矩阵;Q为控制点向量矩阵。

2)基于预测状态x${}_{k/k-1}^{(i)}$进行观测得到观测新息v${}_k^{(i)}$

3)应用新息矩阵H计算卡尔曼更新

S${}_k^{(i)}$=HP${}_k^{(i)}$HT+σ2Im (3)

$\mathit{{\rm K}}{}_k^{(i)}=\stackrel{⌢}{\mathit{{\rm P}}}{}_k^{(i)}\mathit{{\rm H}}{}^{{\rm T}}(\mathit{S}{}_k^i){}^{-1}$ (4)

$\stackrel{⌢}{x}{}_{k|k}^{(i)}=\stackrel{⌢}{x}{}_{k|k-1}^{(i)}+\mathit{{\rm K}}{}_k^{(i)}v{}_k^{(i)}$ (5)

$\mathit{{\rm P}}{}_k^{(i)}=\stackrel{⌢}{\mathit{{\rm P}}}{}_k^{(i)}-\mathit{{\rm K}}{}_k^{(i)}\mathit{S}{}_k^{(i)}(\mathit{{\rm K}}{}_k^{(i)}){}^{\text{Τ}}$ (6)

其中,S为函数空间,H为新息矩阵,

$\mathit{{\rm H}}=\left[\begin{array}{c}\mathit{B}(\mathit{S}){}^{\text{Τ}}\\ \mathit{B}(\mathit{S}){}^{\text{Τ}}\end{array}\right],\mathit{B}$

为B样条基函数,σ2为方差,K为卡尔曼增益。

4)从下面的提议密度中采样得到粒子x${}_k^{(i)}$

$\text{π}(x{}_k^{(i)}|\tilde{x}{}_{k-1}^{(i)},y{}_{1:k})=\mathit{{\rm N}}(\stackrel{⌢}{x}{}_{k|k}^{(i)},{\rm P}{}_k^{(i)})$ (7)

其中,N为提议密度矩阵。

(3)计算新的权重。

(4)输出阶段。

(5)重采样阶段。

(6)k=k+1,若继续跟踪,转到步骤(2);否则,算法结束。

在卡尔曼粒子滤波器中,在采样阶段每个独立的粒子按照高斯过程进化,这个过程同时考虑了状态转移密度和观测密度。高斯提议密度表征了随机粒子的一阶和二阶统计特性,从该提议密度中进行采样得到新粒子与父粒子具有相似的分布,因此继承了父粒子的协方差。每个粒子的状态和协方差可以看作后验密度中的一个高斯分量密度,并按照卡尔曼滤波器进化。

在粒子滤波器中,随着粒子数目的增加,算法的性能逐渐改善,但计算量也变得越来越大。因此,在实际的目标跟踪中,需要对算法的性能和计算效率进行分析,找到临界点。下面就通过临界粒子数目对3种不同的跟踪器进行性能和效率比较。

当粒子数为500时,卡尔曼滤波和粒子滤波都取得了较好的结果。对于这两种跟踪器,逐渐减少其粒子数,直到得到临界粒子数目。采用一个经典的非线性例子[8]进行仿真,比较结果如表1和图2所示。

状态方程

xk=0.5xk-1+25xk-1/(1+x${}_{k-1}^2$)+8 cos 1.2k-1+rk (8)

观测方程

pk=x${}_{k-1}^2$/20+qk (9)

其中,rk和qk分别代表系统噪声和观测噪声,且都是均值为零的高斯白噪声。

由表1可知,粒子滤波和卡尔曼粒子滤波两种算法在各自的临界数目下跟踪误差大致相同,而卡尔曼粒子滤波临界粒子数远少于其他两种滤波,同时每帧跟踪时间也大幅减少。

由图2(a)可知,此时的卡尔曼滤波算法已经不能得到准确结果,而卡尔曼粒子滤波的结果明显更接近于真实值。图2(b)中卡尔曼粒子滤波的估计值误差远低于其他两种跟踪算法,该算法的执行效率和跟踪误差都优于卡尔曼滤波器和粒子滤波器。图3显示卡尔曼粒子滤波器的估计值都在真实值的95%置信区间内。综上可知,基于卡尔曼粒子滤波器的算法能够满足本系统的实时性和鲁棒性的要求。然而,卡尔曼粒子滤波器在硬件系统中实现起来比其他两种算法更复杂。

3 实验结果分析

实验素材选取于一段事先拍摄好的视频。由图4的跟踪结果可以看出,检测运动目标的标识符为白色“十”字符,而该 “十”字符中心点基本能够与目标的中心点重合,达到了预期的目标。由此可知,本文采用的基于卡尔曼粒子滤波的跟踪算法能较好地实现运动目标的跟踪,有效达到预期目的,有良好的实际应用价值和推广价值。

4 结束语

主要研究了基于滤波理论的目标跟踪方法,首先简单介绍了目前常用的基本算法,并对基于滤波理论的3种算法进行了比较分析,最后在Matlab平台上进行了仿真实验,实验结果表明Kalman 粒子滤波算法能够很好地对目标进行跟踪,并在实时性、鲁棒性等方面较为突出,比其他算法更有效。

参考文献

[1]刘卫光,李广鑫.一种通用的视频目标跟踪系统设计[J].计算机技术与发展,2009,19(10):110-113.

[2]NIXON M S,HARRIS C J.Statistical gait description viatemporal moments[C].Proceedings of Proc 4th IEEE South-west Symposium on Image Analysis and Interpretation,2000:291-295.

[3]王利循.无迹Kalman滤波在IMU和GPS组合导航系统中的应用研究[D].南昌:南昌大学,2010.

[4]程娟.复杂背景下运动目标识别算法研究[D].武汉:武汉理工大学,2008.

[5]ZHOU Feng,MENG Xiuyun.Adaptive kalman filter of transferalignment with un-modeled wing flexure of aircraft[J].Jour-nal of Beijing Institute of Technology,2008,17(4):434-438.

[6]朱娟.蒙特卡洛滤波算法在目标跟踪中的应用[D].长春:中国科学院研究生院,2010.

[7]LI P H,ZHANG T W,ARTHUR E C P.Visual contour track-ing based on particle filters[J].Image and Vision Compu-ting,2003(21):111-123.

基于卡尔曼滤波的相对定位算法研究 第8篇

随着移动通信系统[1]的不断发展,用户对移动通信平台的功能需求也逐渐提高。利用现有移动通信平台的测距能力,可以进行进一步功能升级,使移动平台与固定基站之间除了实现稳定、可靠的通信功能之外,还可实现它们之间的相对定位功能。

为了实现移动平台与固定基站之间的相对定位,本文结合卡尔曼滤波技术在定位方面的优点,研究了基于卡尔曼滤波的相对定位算法,给出了算法设计的详细过程,并通过构建仿真场景,对算法进行了仿真验证。

二、相对定位原理

相对定位需要先建立相对坐标系。首先指定一个固定基站作为中心站。相对坐标系的原点由它指定。然后,它估计出其在相对坐标系中的坐标,同时将包括相对坐标、地理位置等数据的位置消息广播出去。已入网的其他站点接收到该位置消息后,可以估计自身在相对坐标系中的坐标。但此坐标位置还需要不断的校正,以提高其精度,校正的过程如下:

进行相对定位的移动平台会周期的接收到固定基站的位置消息。移动平台选择其中几个基站,利用接收到的位置消息中基站的位置、自身的位置估计数据,计算出移动平台与基站之间的距离估计值。而移动平台在接收位置消息时,能得到准确的测距值(TOA),从而计算得到与基站之间的观测距离。将距离估计值与观测距离进行比较,可以计算出定位误差,从而改善位置的估计值。

三、相对定位算法设计

相对定位算法采用以测距值(TOA)作为观测值的扩展卡尔曼滤波算法[2]。该算法包括2个基本的方程:a)状态方程b):观测方程:

其中,Xk是状态向量,Φk,k-1是系统的状态转移矩阵,ΓK,K-1是系统噪声输入矩阵,WK是系统噪声序列,Zk是观测噪声序列,Hk是观测矩阵,Vk是观测噪声序列。

观测量是通过TOA计算得到的伪距R0即移动平台与基站的观测距离。观测模型可以表示为:

其中,c为光速,TOA为相对于移动平台本地时钟而言,观测到的到达时间,Rc为根据对自身位置的最佳估值,以及基站发送的位置消息中包含的基站位置信息计算得到的距离,bt为基站的时钟偏差,b为移动平台自身的时钟偏差,N为所有的观测噪声和,xt、yt、zt为基站的位置坐标。

R0是关于状态向量的非线形方程,用R0对状态向量取偏导即可得到线形化观测矩阵H。

除了上述2个基本方程外,扩展卡尔曼滤波算法还包括以下变量:

a)状态向量的预测值b)状态向量的估计值

基于卡尔曼滤波的相对定位算法[3]的基本流程为:a)为状态向量、误差协方差矩阵指定初值;b)预测下一时刻的状态向量、误差协方差矩阵;c)下一时刻到来时,用收到的位置消息中基站的位置、预测的移动平台位置计算预测的Rc,并求得观测矩阵;d)用基站的位置误差方差扩充预测的误差协方差矩阵,并结合观测矩阵计算滤波器增益矩阵;e)计算R0的观测值和预测值:f)计算状态向量和误差协方差矩阵的估计值。

上述b)至f)即为卡尔曼滤波器的一次状态更新过程,不断重复该过程即可使各种状态得到周期性的更新。

四、算法仿真

4.1仿真场景

仿真场景包括宽范围和窄范围场景。宽范围场景中的各个源相对距离较大,位置关系较好。窄范围场景中的各个源相对距离较小,位置关系较差。下面分别进行描述。

宽范围场景的条件如下:场景中包括4个固定基站G1、G2、G3、G4,以及4个移动平台A21、A23、A24、A32,4个移动平台速度恒定,A21的运动轨迹为匀速圆周运动,速度为300ft/s。A32的运动轨迹为匀速直线运动,速度为455ft/s。每个平台周期发送位置消息,每个平台的初始时钟相位偏差值平均分布在0~1.8ms之间,初始时钟频率偏差值服从高斯分布,标准差为10-8。窄范围场景跟宽范围场景相比,除了成员位置分布和运动轨迹不同,其他条件一样。它们的平台位置分布如图1所示。

4.2仿真结果

根据上述宽范围仿真场景,得到的仿真结果如下:

A32的仿真曲线图如2~3所示。在图中,RMSx曲线图的纵轴是滤波器每次更新后得到的x坐标均方根误差值,横轴是更新次数。RMXy曲线图的纵轴是滤波器每次更新后得到的y坐标均方根误差值,横轴是更新次数。

A21的仿真曲线图如图4~5所示:

根据上述窄范围仿真场景,得到的仿真结果如下:

A32仿真曲线图如图6~7所示:

A21仿真曲线图如图8~9所示:

五、结束语

本文针对移动通信平台与固定基站之间的相对定位功能需求,对基于卡尔曼滤波的相对定位算法进行了较深入的研究,详细介绍了算法的具体设计过程,并通过特定的应用场景,对算法进行了仿真分析,验证了算法的有效性,结果表明平台间位置关系对定位精度有较大的影响,位置关系较好时得到的定位精度比位置关系较差时的精度要高,从而为实际项目提供了一定的理论参考。

摘要：根据移动通信网络对平台之间相对定位的应用需求,研究了基于卡尔曼滤波的相对定位算法,给出了具体的算法设计过程,并对算法进行了仿真分析,验证了算法的有效性,仿真结果为实际工程提供了一定的理论参考。

关键词：相对定位,卡尔曼滤波

参考文献

[1]韩斌杰.《GSM原理及其网络优化》.机械工业出版社,2003.

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

扩展卡尔曼滤波算法 第9篇

Paliwal提出的时域卡尔曼(Kalman)滤波语音增强[1]方法通过状态空间模型描述含噪语音的进化过程,对语音信号预测与估计获得MMSE估计值,非常适合语音信号的非平稳特性,不存在音乐噪声,自然度高。近年来,出现了多种基于Kalman滤波的语音增强方法。时域Kalman滤波需要的线性预测阶数多,算法复杂度高。文献[2,3]提出子带Kalman滤波,降低了一定的阶数,计算量仍然比较大。Zavarehei等人尝试用AR模型描述短时傅立叶变换(DFT)域的频谱进化过程,并提出在短时DFT域对含噪语音进行Kalman滤波[3]。该方法需要的阶数较少,但要同时对DFT的实部和虚部系数运用Kalman滤波。文献[4]先用Ephraim的MMSE方法估计语音信号的对数DFT幅度谱,再串接Kalman滤波进行,显然该方法复杂度高。

本文提出利用离散余弦变换(discrete cosine transformation,DCT)去除帧内语音相关性,然后在每个语音通道分别对DCT系数进行Kalman滤波获得语音。由于DCT系数为不相关的实数,省去了文献[5]在虚部也要Kalman滤波的运算,同时用于描述状态空间方程所需要的阶数少的特点得以保留,适合并行处理。符合统一计算设备架构(compute unified device architecture,CUDA)计算的图形处理单元(graphics processing unit,GPU)具有上百个并行计算单元,大大缩短运算时间。为了提高该方法的工程实用性,本文实现了基于CUDA平台的并行设计。实验结果表明:DCT-Kalman语音增强方法可以有效去除噪声还原语音信号,基于CUDA平台的实现可以有效地缩减DCT-Kalman语音增强程序的运算时间。

1 DCT-Kalman滤波语音增强

对时域被加性噪声污染的带噪语音信号进行分帧后,对应的带噪语音信号的DCT系数、纯净语音信号DCT系数与观察噪声的DCT系数分别设为Y(j,k)、X(j,k)与V(j,k),其中j为系数索引,k为时间帧索引。(为了方便表示,省去j。)

由于DCT具有正交分解特性,每帧语音的DCT系数是相互独立的,而信号帧之间仍然存在较强的相关性,可用AR模型描述语音DCT系数的进化过程

其中,ω(K)是某一频率通道,第k帧的均值为0,方差为的过程噪声;a(k)=[a1(k),…,ap(k)]是DCT系数的AR系数向量,由YuleWalker方程可得

其中,Rx(k)是X(k)的自相关矩阵,r(k)是由X(k)组成的向量。

Kalman滤波方法利用状态方程和测量方程描述随机系统的时间演进,使得该滤波方法能够适用于非平稳条件,结合语音AR模型,用状态空间描述含噪语音、测量噪声与纯净语音之间的关系如下

其中,F为AR系数构成的状态转移矩阵;

在同一频率通道中,S(k)是K时刻的语音DCT系数向量,表示为S(k)=[X(k-p+1),…,X(k)]T;H=GT=[0,…,1]1×p;Y(k)是含噪语音DCT系数。

用于语音增强的DCT-Kalman滤波的状态估计与更新过程可以用式(5)、(6)、(7)、(8)、(9)表示

其中,S(k|K-1)与S(k|k)分别表示纯净语音的预测与估计;分别表示过程噪声方差和环境噪声方差;K(k)表示增益;P(k|k-1)表示预测误差协方差矩阵;P(k|k)表示估计误差协方差矩阵。

最后,对估计出来的语音估计S(k)做N点的反离散余弦变换再去窗函数,去重叠得到增强语音x(n)。

2 DCT-Kalman滤波的CUDA并行分析与实现

CUDA模型如图1所示,由多个SP组成一个多处理器(MultiProcessor),同一个多处理器里的所有SP在每一步都执行相同的运算操作。多处理器上设置有数千个读写速度较快的寄存器,还有容量小,读写速度比寄存器慢一些的共享内存。GPU架构不像CPU架构那样具有高速缓冲器,所以CUDA程序必须让线程以有规律和线性的方式访问显存,否则会导致读写显存的速度很慢。此外还有一块容量小、只读的常量内存,用于存放程序中不变的常量,可以加快读取速度。

CUDA程序需要主机函数和核函数。主机函数由CPU执行,负责初始化,传输数据到GPU端,启动核函数,最后传输数据回CPU端。核函数则以线程块的形式在GPU的多处理器上执行,即每个SP上执行一个线程。CUDA并行处理调用了CUDA库的DCT核函数实现转换,对转换后的每个频率通道的语音信号,调用Kalman核函数对其处理。N个通道信号需要调用N个Kalman核函数处理,每个Kalman核函数1个线程,共N个线程执行,分别占用一个多处理器里面的N个SP。这样就不需要外循环来控制频率通道,实现了并行处理。由于式(5)、(6)、(7)、(8)、(9)描述的Kalman滤波在每个频率通道内随时间进行迭代,其中的若干个向量与矩阵变量需要不断更新,为了提高这类变量的读写速度,使用共享内存存放。此外,为了获得内存的合并访问,必须对每个频率通道的相同变量都开辟对应的存储空间。如:p=2的2×1的S(k|k)存放在2行N列的共享内存数组里,其中列对应频率通道索引:

2×2的矩阵F存放在2行(p*N+1)列的共享内存数组里,存储示意(数字的含义同上):

其中每连续p列依次对应一个频率通道F,最后1列不存储任何数据,仅用来错开数据存放的bank,避免一个warp(32个线程)访问共享内存时产生bank conflict。变量中,向量都是类似S(k|k)一样存储,矩阵都是类似行F一样存储。迭代过程中不变的常量,如式(6)的G矩阵、式(7)的H矩阵则存在常量内存上。过程噪声方差的值则存放在寄存器上。最后Kalman核函数的伪代码如下所示:

3 实验与结果分析

实现Kalman语音增强需同时获得测量噪声的方差,语音的AR系数向量a(k)与过程噪声方差。实际应用时间必须结合语音活动检测、噪声估计方法和谱减法以获得以上参数的估计才能实现去噪。本文主要研究DCT-Kalman与时域Kalman的去噪效果与时间复杂度,因此没有研究以上3个参数在含噪语音中的获得方法。在实验中,直接从纯净语音中获得AR系数与过程噪声方差;直接用已知噪声求出测量噪声的方差作为实验所用。

实验用Intel I3-530CPU和NVIDIA GTX650TI BOOST显卡硬件平台,在Matlab 2013环境下分别运行基于CPU的时域Kalman滤波算法、与DCT域Kalman滤波算法和基于GPU的DCT域Kalman滤波算法,比较它们的去噪性能与运行时间。实验中的语音与带噪语音均为16 kHz采样,选取汉明窗为窗函数,帧长N取128,帧移为50%,即64。时域Kalman滤波程序的语音AR模型阶数p取13,并在一帧信号进行估计;DCT-Kalman滤波的语音AR模型阶数P取2。并在3帧信号中估计。噪声估计的平滑系数设为0.95。

由图2、图3与图4最下面和中间的圈可以比较出,DCT域理想Kalman滤波后的语谱图纹路比时域理想Kalman滤波后更为清晰,由最上面的圈可以看出,DCT域比时域的能量稍微低一些,恢复得更接近纯净语音。

从表1分段信噪比来说,同一语音下,DCT域理想Kalman滤波的分段信噪比比时域理想Kalman滤波要高,且所需的阶数仅为2。从表2来看,用CUDA来做DCT-Kalman的速度比CPU快200倍。

4 结语

CUDA平台适用于并行计算。利用CUDA架构设计Kalman语音增强方法,运算速度提高200倍,充分发挥了GPU的并行处理能力,节省运算时间。依靠GPU的并行计算,语音增强这一类计算复杂度高的科学计算问题将得到具有实际工程意义的解决。

摘要：提出一种卡尔曼滤波语音增强算法的统一计算设备架构并行实现方案。该方案通过离散余弦变换把含噪语音分解为不相关的DCT系数,使原来的时域串行处理转化为统一计算设备架构并行处理,只需较少的线性预测阶数,节省运算时间。实验结果表明:与时域卡尔曼滤波比较,该方法有更高的输出分段信噪比,其统一计算设备架构加速方案可有效缩短语音增强运算的时间。

关键词：语音增强,卡尔曼,离散余弦变换,统一计算设备架构

参考文献

[1]Paliwal K,Basu A.A speech enhancement method based on Kalman filtering[C].Proc.IEEE ICASSP,1987:177-180.

[2]Wu Wen-Rong,Chen Po-Cheng.Subband Kalman filtering for speech enhancement[J].IEEE Transactions on Analog and Digital Signal Processing,1998,45(8):1072-1083.

[3]章旭景,李辉,陆伟.基于子带卡尔曼滤波的语音增强方法[J].信号处理,2009,25(9):1474-1478.

[4]齐贺,张雪英,武奕峰.频域内的卡尔曼滤波语音增强算法[J].电声技术,2011,35(10):55-58.