基于扩展卡尔曼滤波的轨道垂向不平顺估计

?

基于扩展卡尔曼滤波的轨道垂向不平顺估计

王　贵1，邢宗义1，王晓浩2，陈岳剑2

(1.南京理工大学自动化学院，南京　210094;2.南京理工大学机械工程学院， 南京　210094)

摘　要：轨道不平顺是影响列车平稳性和舒适度的关键因素，因此及时掌握线路状态对保证列车的运行安全具有重要意义。针对采用单个惯性量较难达到对不同波段不平顺的检测，通过观测多个惯性量，运用扩展卡尔曼滤波解决非线性离散系统的最优估计原理，根据车辆轨道耦合状态空间方程计算递推雅克比矩阵，并结合线性观测方程得到最优状态估计，实现轨道不平顺估计。在Matlab平台下，进行了实测轨道不平顺激励作用下的仿真，将仿真得到的观测值采用本文提出的方法进行轨道垂向不平顺估计，结果表明该算法具有很好的精确性。

关键词：轨道不平顺；扩展卡尔曼滤波；最优估计

轨道不平顺是引起车辆与轨道产生振动的主要激励源之一，是轨道质量状态的集中体现。轨道不平顺不仅会增加车辆与轨道间的相互作用，缩短车辆与轨道的使用寿命，影响乘坐舒适度，当形变累积到一定程度时还会使车辆产生倾斜和侧滚运动，严重威胁列车的安全运行，因此对运行线路的轨道不平顺状态进行实时检测和估计具有重要意义[1-2]。

目前已见众多轨道不平顺在线检测和估计的研究。朱文发等[3]提出了基于捷联惯性系统实现轨道长波不平顺的检测方法，并采用小波去噪和积分滤波器将加速度信号处理得到轨道长波不平顺位移。J Real等[4]在轴箱上安装垂向加速度计，获得轴箱振动加速度信号，并对信号进行二次积分、高通滤波、相位补偿和振动模型逆输入的处理，实现对轨道垂向不平顺的检测。Lee J S等[5]利用轴箱加速度通过卡尔曼滤波和波长带通滤波来估计轨道不平顺。

轨道垂向不平顺包括波长较广范围内的不平顺，常用检测方法中，轴箱振动较难同时检测局部不平顺和长波不平顺；构架振动由于一系簧滤波特性和几何滤波而丢失了较多波长成分；构架倾角对短波不平顺存在空间滤波特性；因此仅以轴箱垂向振动、构架垂向振动或构架点头角速度中单个惯性量较难达到对不同波段不平顺的检测[6-7]。本文中综合考虑振动加速度、构架角速度等多个惯性量，构建以多个惯性响应为输入，以轨道垂向不平顺为输出的数学模型，以达到精确的轨道垂向不平顺估计。

惯性响应与轨道垂向不平顺的关系是复杂多输入非线性动态[8]。常用的非线性滤波方法有扩展卡尔曼滤波(EKF)[9]、无迹卡尔曼滤波(UKF)[10]、粒子滤波(PF)[11]等。扩展卡尔曼滤波具有实时性强、精确度高等优点[12]，对非线性系统滤波能得到近似的可接受的解，因此采用扩展卡尔曼滤波算法实现多惯性监测量与轨道不平顺之间的滤波。

1　扩展卡尔曼滤波算法

对于如下离散非线性系统模型

(1)

(2)

(3)

式(1)为系统的状态方程，式(2)为系统的输出观测方程。f()、h()为非线性传递函数，wk和vk代表过程模型激励噪声和观测噪声，两个噪声对应的协方差分别为Qk，Rk。xk为系统的状态变量，yk为系统的观测变量，uk为系统的驱动函数。

将状态方程在和wk-1=0点处、观测方程在和vk=0点处进行泰勒级数展开，忽略高阶项

(4)

(5)

对偏微分雅克比矩阵Fk-1、Lk-1、Hk、Mk进行计算，状态估计和估计误差协方差的时间更新如下

(6)

(7)

状态估计的测量更新和估计误差协方差更新如下

(8)

(9)

(10)

时间更新方程(6)、(7)将状态估计和协方差估计从k-1时刻向前推算到k时刻。测量更新方程(8)，首先计算卡尔曼增益Kk，其次获得测量输出yk，然后按式(9)产生状态的后验估计，最后按式(10)估计状态的后验协方差。整个过程在计算完时间更新方程和测量更新方程后重复迭代。

2　基于EKF的轨道垂向不平顺估计

通过车辆轨道耦合模型的动力学方程可以改写得到系统的状态方程，结合传感器的观测方程，经过EKF滤波迭代后得到最优状态，然后根据最优状态逆向计算轨道不平顺。

2.1　车辆轨道耦合系统状态方程

车辆轨道耦合模型的动力学方程可描述为

(11)

式中，M、C、K为车辆的质量、阻尼和弹簧矩阵与轨道的质量、阻尼和弹簧矩阵结合；Q为车辆与轨道动力学系统激励的组合；状态向量x为车辆与轨道动力学系统状态向量的组合[8]

(12)

式中，Zc、Zt1、Zt2、Zw1、Zw2、Zw3、Zw4、Zsj、Zbi分别为车体、前构架、后构架、第1～4轮对、第i轨枕、道床的浮沉运动；βc、βt1、βt2分别为车体、前构架、后构架的点头运动；NM为钢轨模态阶数；N为轨枕支点总数；qk(t)为钢轨正则振型坐标。

对于车辆轨道垂向耦合模型，在仿真计算轨道长度l取70 m的情况下，N=100，根据文献[8]，NM合理取值为90，整个系统的自由度为300，即x为300×1的系统状态向量，M，C，K为300×300的常量矩阵。

将式(11)改写成式(13)

(13)

式中，；；u(t)为轨道不平顺；t为时间；w(t)为模型噪声，A为600×600的常系数矩阵。

2.2　传感器观测方程

系统测量模型采用在车体、转向架构架上安装振动加速度传感器，获得车体、构架的垂向振动加速度t；在车体、构架上安装陀螺仪，获得车体点头角速度c，构架点头角速度t2，观测方程如式(14)所示。

(14)

式中，v(t)为测量噪声；H为测量矩阵；X(t)为t时刻状态向量及其导数的组合。

注意到式(14)包含额外的重力加速度，这部分可视为测量噪声并非为零均值的，本文采用对检测到的加速度需先减去g的处理方式。

2.3　滤波迭代方程配置

对系统状态方程式(13)求取雅克比矩阵Fk-1

为k-1时刻的状态后验估计fk-1为k-1时刻的X(t)和u(t)的函数关系式，忽略下标 ，定义

(16)

式(15)中A为已知的常数矩阵，推导可得Au

(17)

其中

(18)

状态估计和估计误差协方差的时间更新如下

(19)

(20)

式(19)中，Pk为状态协方差；为k时刻先验协方差估计；Fk-1为k-1时刻雅克比矩阵；为k-1时刻后验协方差估计；Qk-1为k-1时刻轮轨力。

式(20)状态估计时间更新采用经典的4阶龙格-库塔积分实现[13]。

式(20)中，为状态的先验估计；uk-1为k-1时刻的轮轨接触力为k-1时刻的X(t)和u(t)的函数关系式。

对观测方程(14)求取雅克比矩阵

式中，Hk和Mk为hk关于x和v的偏微分；为k时刻状态的先验估计。

状态估计的量测更新和估计误差协方差的更新如下

(21)

(22)

(23)

式中，Kk为卡尔曼增益；Rk为量测噪声的协方差；为k时刻状态的后验估计；为k时刻状态的先验估计；为k时刻后验状态协方差估计。

2.4　由最优状态计算轨道不平顺

经过EKF滤波迭代后得到的各时刻最优状态，利用状态估计的时间更新方程的式(20)，求逆后即可得到轨道不平顺激励k，即

(24)

式中，uk-1为k-1时刻的轮轨接触力；k-1为k-1时刻的状态估计；k为k时刻的状态估计。

改写状态转移方程式(13)如下

(25)

式中，τ为常数矩阵。由于轨道不平顺在不同车轮的激励输入仅是时间延迟，因此只需得到首个车轮Zw1受到的激励即可。提取第307个行元素，建立等式如下

(26)

最终，轨道不平顺的反推方程如下

(27)

3　仿真验证

基于EKF算法的轨道不平顺估计研究框架汇总见图1，通过车辆轨道耦合动力学模型仿真得到真实的状态响应，采集可观测的状态如车体加速度、车体角速度、构架加速度和构架角速度等，利用扩展卡尔曼滤波，结合观测方程、状态方程和由估计状态逆向计算轨道不平顺实现轨道不平顺的最优估计。

图1　基于EKF算法的轨道不平顺估计研究框架

图2　高铁实测轨道不平顺

通过一段实测轨道高低不平顺数据对车辆-轨道垂直耦合动力学模型进行数值求解。该段数据来自武广高速铁路DK1 534～DK1 536的运营线路，利用GJ-6型轨检车基于激光检测法检测输出。该段数据是以0.1 m为间隔的离散采样序列，长度2 km，如图2所示。动力学模型仿真可直接获取车体振动加速度、车体点头角速度、构架振动加速度和构架点头角速度，仿真响应如图3所示。

图4　最优状态量与仿真值比较空间域波形和功率谱(一)

图3　观测量空间域波形

滤波器初始化过程中，定义=0。由于初始协方差P0的选择并不关键，通常P0≠0的任何数值均可，因此本文P0=I。此外，注意到式(20)的状态估计时间更新包含轨道不平顺激励uk-1的输入，uk-1可设为0。在实际监测过程中，uk-1可以是上一次的轨道不平顺值，在设备刚投入使用阶段，无法避免要面临uk-1未知的情况，因此本文设定式(20)的状态估计时间更新中uk-1=0，这会给式(20)的状态估计带来一定的不确定度。

滤波器迭代过程中，设定式(20)的状态估计的积分步长为0.1 ms，滤波器迭代步长1 ms，空间步长0.1 m。滤波后车辆系统状态量响应与动力学仿真输出的真实状态对比见图4和图5。图4为车轮加速度，构架加速度和车体加速度最优状态量与仿真值比较空间域波形和功率谱(一)。图5为车轮速度，构架角速度，车体角速度最优状态量与仿真值比较空间域波形和功率谱(二)。

图5　最优状态量与仿真值比较空间域波形和功率谱(二)

为了量化衡量图4和图5中滤波后车辆系统状态量响应与动力学仿真输出的真实状态逼近程度，计算两者之间的相关度系数R，标准差SD和归一化均方误差NMSE。仿真过程之前的状态量均添加了协方差为Rk的观测噪声，因此为更直观地展示拓展卡尔曼滤波的精度，对真实状态值添加噪声后计算相关度系数R，标准差SD和归一化均方误差NMSE。表1给出了滤波前后，车轮加速度等6个状态变量的R，SD和NMSE。对比表1中各个状态变量在滤波前后R，SD和NMSE，可发现滤波后相关系数R明显增大，均高于0.98，表现了估计值和仿真值之间很高的相关程度；标准差SD和归一化均方误差NMSE均有减小，且SD多数小于0.1 m/s2，表明估计值自身的变化程度很小，NMSE均低于-10，表明估计值与仿真值的误差很小。综合滤波后各个状态变量的R，SD和NMSE，各个状态量都表现出了很高的逼近程度，有效地抑制了观测噪声带来的影响，得到了状态量可接受精度的最优估计。

表1　最优状态量与真实值逼近程度量化指标

状态变量滤波前滤波后RSDNMSERSDNMSE车轮加速度0.6474.048m/s21.440.9970.339m/s2-19.89构架加速度0.6620.474m/s21.080.9830.118m/s2-15.33车体加速度0.7860.099m/s2-2.030.9980.069m/s2-20.93车轮速度0.5870.041m/s2.830.9970.0031m/s-19.42构架角速度0.7460.562°/s-0.980.9980.056°/s-10.44车体角速度0.8160.110°/s-3.020.9990.083°/s-21.28

按式(24)由EKF得到的最优状态，计算轨道不平顺。得到的轨道不平顺与真实轨道不平顺激励对比见图6，可见最优状态计算得到的轨道不平顺在各个波段均逼近真实值，尤其在2～300 m波长范围内功率谱值与真实不平顺表现了很强的一致性。滤波后相关系数R为0.984 3，标准差SD为0.171 mm，归一化均方误差NMSE为-27.2，表明运营列车在广泛的车速区间内均能实现垂向不平顺的预测，验证了所提扩展卡尔曼滤波方法的有效性。

图6　滤波后轨道不平顺与实际激励值比较空间域波形和功率谱

4　结论

本文提出了一种基于扩展卡尔曼滤波的轨道垂向不平顺估计方法。由于采用单个惯性量较难达到对不同波段不平顺的检测，本文通过观测多个惯性量，并结合车辆轨道耦合状态空间方程，基于扩展卡尔曼滤波实现了轨道垂向不平顺的估计。仿真结果表明本文所提方法的精确性和鲁棒性。

参考文献：

[1]　高春雷，王发灯.利用激光准直技术检测线路的长波不平顺[J].铁道建筑，2009(1):81-85.

[2]　张念，朱焱，王小敏，等.铁路轨道不平顺预测模型研究与应用[J].铁道标准设计，2015，59(9)：56-63.

[3]　朱文发，柴晓冬，郑树彬，等.基于捷联惯性系统的轨道长波不平顺检测[J].城市轨道交通研究，2012(11)：87-90.

[4]　J Real， P Salvador， L Montalbn， M Bueno.Determination of rail vertical profile through inertial methods[J]. Proceedings of the Institution of Mechanical Engineers， Part F: Journal of Rail and Rapid Transit， 2011，225(1):14-23.

[5]　Lee J S， Choi S， Kim S S， et al. Estimation of Rail Irregularity by Axle-Box Accelerometers on a High-Speed Train[M]//Noise and Vibration Mitigation for Rail Transportation Systems. Springer Japan， 2012:571-578.

[6]　Jun Seok Lee. Sunghoon Choi，Sang-Soo Kim，Choonsoo Park，Young Guk Kim. A Mixed Filtering Approach for Track Condition monitoring using acceleration on the axle box and bogie[J]. IEEE transactions on instrumentation and measurement， 2012，61(3):749-758.

[7]　Real JL， Montalban L， Real T， et al. Development of a system to obtain vertical track geometry measuring axle-box accelerations from in-service trains[J]. Journal of Vibroengineering， 2012，14(2):813-826.

[8]　翟婉明.车辆-轨道耦合动力学[M].北京：科学出版社，2007.

[9]　雷雨龙，付尧，刘科，等.基于扩展卡尔曼滤波的车辆质量与道路坡度估计[J].农业机械学报，2014，44(11):9-13.

[10]解少博，林程.基于无迹卡尔曼滤波的车辆状态与参数估计[J].农业机械学报，2012，42(12):6-12.

[11]陈龙.基于粒子滤波的城轨车辆状态监测方法研究[D].上海：上海工程技术大学，2013.

[12]Simon D. Optimal state estimation: Kalman， H infinity， and nonlinear approaches[M]. John Wiley & Sons，2006.

[13]成棣，刘金朝，王成国.基于Kalman滤波器的故障诊断方法及其在铁道车辆中的应用[J].铁道机车车辆，2007(1):9-13.

收稿日期：2015-09-14； 修回日期：2015-11-26

基金项目：中央高校基本科研业务费专项资金项目(30920130132002)

作者简介：王　贵(1990—)，男，硕士研究生，主要从事轨道交通控制与安全研究，E-mail:[e-mail]。通信作者：邢宗义(1974—)，男，副教授，主要从事交通运输信息工程与安全保障研究，E-mail:[e-mail]。

文章编号：1004-2954(2016)07-0014-05

中图分类号：U213.2

文献标识码：A　　

DOI:10.13238/j.issn.1004-2954.2016.07.004

Prediction of Vertical Track Irregularities Based on Extended Kalman Filter

WANG Gui1， XING Zong-yi1， WANG Xiao-hao2， CHEN Yue-jian2

(1.School of Automation， Nanjing University of Science and Technology， Nanjing 210094， China;2.School of Mechanical Engineering，Nanjing University of Science and Technology， Nanjing 210094， China)

Abstract： Track irregularities are main factors affecting stability and comfort of trains and it is important to understand the status of line to ensure safe operation of trains. Due the difficulties to detect different irregularities of different bands with a single inertial value， the method based on the observation of multi inertia values to predict vertical track irregularity is proposed. The optimal estimation principle of extended Kalman filter for nonlinear discrete systems is used to estimate track irregularities according to Recursive Jacobi matrix， the vehicle track coupling status-space equation and the optimal state estimation obtained with linear measurement equation. On Matlab platform， the simulation of actual track irregularities measured under stimulation is conducted. The simulation results show that this algorithm is accurate to estimate track irregularities．

Key words： Track irregularities; Extended Kalman Filter; Optimal estimation