CN104878668A - 基于扩展卡尔曼滤波的轨道垂向不平顺估计方法及系统 - Google Patents

Abstract

本发明公开了一种基于扩展卡尔曼滤波的轨道垂向不平顺估计方法及系统。该系统在车体和构架上均安装振动加速度传感器和陀螺仪，将振动加速度传感器、陀螺仪连同车速传感器输出端均接入中央处理单元，以车速信号实现车体垂向振动、车体点头角速度、构架垂向振动和构架点头角速度的等空间采样。根据车辆轨道耦合模型的动力学方程构建车辆轨道耦合系统状态方程；根据所有传感器采集到的信号，构建传感器观测方程；配置滤波迭代方程，得到状态估计和估计误差协方差的时间更新方程和量测更新方程；根据扩展卡尔曼滤波得到的各时刻最优状态，利用状态估计的时间更新方程求逆后最终得到轨道不平顺估计。本发明具有成本低、工程实施性好、操作简单和在线实时检测的优点。

Description

基于扩展卡尔曼滤波的轨道垂向不平顺估计方法及系统

技术领域

本发明涉及轨道垂向不平顺在线监测的技术领域，特别是一种基于扩展卡尔曼滤波的轨道垂向不平顺估计方法及系统。

背景技术

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

朱文发公开了一种(朱文发，柴晓冬，郑树彬等基于捷联惯性系统的轨道长波不平顺检测。城市轨道交通研究，2012)基于捷联惯性系统实现轨道长波不平顺的检测方法，该方法使用采用小波去噪和积分滤波器将加速度信号处理得到轨道长波不平顺位移。但是在车体安装惯性系统，经过了一系簧和二系簧的传递，很难具体的评估轨道状态。LeeJ S公开了一种(Lee J S,Choi S,Kim S S,et al.Estimation of Rail Irregularity by Axle-BoxAccelerometers on a High-Speed Train[M]//Noise and Vibration Mitigation for RailTransportation Systems.Springer Japan,2012)利用轴箱加速度通过卡尔曼滤波和波长带通滤波来估计轨道不平顺的方法。J Real公开了一种(J Real,P Salvador,L Montalbán,MBueno,Determination of rail vertical profile through inertial methods.Proceedings of theInstitution of Mechanical Engineers,Part F:Journal of Rail and Rapid Transit.June 2010.)基于轴箱加速度检测轨道垂向不平顺的方法，该方法通过在轴箱上安装垂向加速度计，获得轴箱振动加速度信号，并对信号进行二次积分、高通滤波、相位补偿和振动模型逆输入的处理，实现对轨道垂向不平顺的检测。因为轨道不平顺包括波长较广范围内的不平顺，常见轨道不平顺检测方法中，轴箱振动较难同时检测局部不平顺和长波不平顺。牛小骥公开了一种基于INS/GNSS的轨道不平顺检测系统及方法，以INS/GNSS组合测量系统作为核心测量设备，固定安装于移动支架上并随移动支架一起在轨道上移动，测量获得轨道的三维位置坐标序列和姿态角序列。根据轨道的三维位置坐标、姿态角及轨距测量值，并结合轨道设计曲线参数，采用基于坐标法的轨道不平顺评估方法评估轨道不平顺。但该方法在系统测量单元成本高且复杂，在滤波算法上使用卡尔曼滤波并不能有效滤除非线性系统的噪声。

发明内容

本发明的目的在于提供一种成本低、工程实施性好的基于扩展卡尔曼滤波的轨道垂向不平顺估计方法及系统，基于车体、转向架构架上的振动加速度和车体、转向架的点头角速度对轨道垂向不平顺进行估计。

实现本发明目的的技术解决方案是：一种基于扩展卡尔曼滤波的轨道垂向不平顺估计方法，包括以下步骤：

步骤1，在车体、转向架构架上安装振动加速度传感器，获得车体、构架的垂向振动加速度，在车体、构架上安装陀螺仪，获得车体点头角速度及构架点头角速度，以车速信号实现车体垂向振动、车体点头角速度、构架垂向振动和构架点头角速度的等空间采样；

步骤2，根据车辆轨道耦合模型的动力学方程构建车辆轨道耦合系统状态方程；

步骤3，根据步骤1中所述振动加速度传感器和陀螺仪采集到的信号，构建传感器观测方程；

步骤4，配置滤波迭代方程，得到状态估计的时间更新方程和量测更新方程、估计误差协方差的时间更新方程和量测更新方程；

步骤5，根据步骤4滤波迭代方程进行滤波迭代后，得到各时刻最优状态，利用状态估计的时间更新方程求逆后得到轨道不平顺估计。

一种基于扩展卡尔曼滤波的轨道垂向不平顺估计系统，包括陀螺仪、振动加速度传感器和中央处理单元，其中陀螺仪安装于列车车体和走行部一系弹簧上方的构架上，振动加速度传感器安装于车体和走行部二系弹簧下方的构架上，所述振动加速度传感器和陀螺仪的输出端均接入中央处理单元，中央处理设置于车体内。

本发明与现有技术相比，其显著优点是：(1)成本低，节约了轨检车的运行、维护和调度成本；(2)工程实施性好，陀螺仪、振动和速度传感器能承受现场的恶劣环境，且安装方便；(3)在线实时监测，可及时发现轨道垂向不平顺状态变化，从而提供及时的维护预警，为铁路维护进行更有效的规划。

附图说明

图1为本发明基于扩展卡尔曼滤波的轨道垂向不平顺估计方法的流程图。

图2为车辆-轨道垂向耦合动力学模型及其参数的示意图。

图3为实测轨道不平顺空间波形图和功率谱图，其中(a)为实测轨道不平顺空间波形图，(b)为实测轨道不平顺功率谱与德国高干扰谱的比较图。

图4为观测量空间域波形图，其中(a)为构架角速度波形图，(b)为构架振动加速度波形图，(c)为车体角速度波形图，(d)为车体振动加速度波形图。

图5为最优状态量与真实值比较空间域波形和功率谱图，(a)为车轮加速度的波形图，(b)为车轮加速度的功率谱图，(c)为构架加速度的波形图，(d)为构架加速度的功率谱图，(e)为车体加速度的波形图，(f)为车体加速度的功率谱图，(g)为车轮速度的波形图，(h)为车轮速度的功率谱图，(i)为构架角速度的波形图，(j)构架角速度的功率谱图，(k)为车体角速度的波形图，(l)为车体角速度的功率谱图。

图6为滤波后轨道不平顺与实际激励值比较空间域波形和功率谱图，其中(a)为波形图，(b)为功率谱图。

具体实施方式

下面结合附图及具体实施例对本发明作进一步详细说明。

结合图1，本发明垂向轨道长波不平顺估计方法，包含以下步骤：

步骤1，在车体、转向架构架上安装振动加速度传感器，获得车体、构架的垂向振动加速度，在车体、构架上安装陀螺仪，获得车体点头角速度及构架点头角速度，以车速信号实现车体垂向振动、车体点头角速度、构架垂向振动和构架点头角速度的等空间采样，具体为：

按照轨道短波不平顺0.5～5m的波长范围和列车1～360km/h的典型车速，加速度计的响应特征频率f_a需包含0～200Hz；同样，按照轨道长波不平顺5～200m的波长范围，陀螺仪的响应特征频率f_θ需包含0～20Hz；由采样定律知，为复现0.5～200m波长的轨道不平顺，令空间间隔为Δs，则Δs<0.25m，f_a>400Hz，f_θ>40Hz，实际应用过程中Δs<0.1m，f_a>1000Hz，f_θ>100Hz；若在时域获得列车实时运行过程中轴箱振动加速度信号a(t)、构架点头陀螺仪输出ω(t)和列车速度v(t)，则需通过a(t)、ω(t)分别除以v(t)得到空间域序列a(s)、ω(s)，从而将a(t)、ω(t)对应到轨道的各个位置；由于目前列车轴端加装的测速传感器为旋转脉冲速度传感器，该传感器基于霍尔原理，采集车轴端固定个数的金属孔产生的脉冲，将采集到的脉冲对时间差分得到速度。因此利用脉冲信号的跳变沿触发陀螺仪信号的采样，即可达到等空间采样的目的，获得以空间间隔为Δs的空间域序列a(s)和ω(s)。

步骤2，根据车辆轨道耦合模型的动力学方程构建车辆轨道耦合系统状态方程，具体为：

车辆轨道耦合模型的动力学方程为：

$M_z\stackrel{\&CenterDot;\&CenterDot;}{x}+C\stackrel{\&CenterDot;}{x}+Kx=Q$

式中，M_z为车辆的质量矩阵与轨道的质量矩阵结合、C为车辆的阻尼矩阵与轨道的阻尼矩阵结合、K为车辆的弹簧矩阵与轨道的弹簧矩阵结合；状态向量x为车辆与轨道动力学系统状态向量的组合：

x＝[Z_c β_c Z_t1 β_t1 Z_t2 β_t2 Z_w1 Z_w2 Z_w3 Z_w4 q_k(t) Z_sj Z_bi]^T(k＝1～NM,j＝1～N)式中，Z_c、Z_t1、Z_t2、Z_w1、Z_w2、Z_w3、Z_w4、Z_sj、Z_bi依次分别为车体、前构架、后构架、第1～4轮对、第i轨枕、道床的浮沉运动，β_c、β_t1、β_t2依次分别为车体、前构架、后构架的点头运动，NM为钢轨模态阶数，N为轨枕支点总数；q_k(t)表示钢轨正则振型坐标_，激励矩阵Q为车辆与轨道动力学系统激励的组合，由下式表示：

$\begin{array}{c}Q={\left[\begin{array}{ccccccccccccc}{M}_{c}g& 0& M_{t}g& 0& M_{t}g& 0& M_{w}g-2p_1\left(t\right)& M_{w}g-2p_2\left(t\right)& M_{w}g-2p_3\left(t\right)& M_{w}g-2p_4\left(t\right)& \underset{j=1}{\overset{4}{\&Sigma;}}{p}_j\left(t\right)Y_k\left(x_{wj}\right)& 0_{Nx1}& 0_{N\&times;1}\end{array}\right]}^T\\ \left(k=1~NM\right)\end{array}---\left(1\right)$

式中，M_c、M_t、M_w分别为车体、构架、轮对的质量，g为重力加速度，Y_k(x_wj)中Y_k为正则振型函数，x_wj为各车轮运动坐标，p₁(t)～p₄(t)为4个车轮的轮轨接触力；由上式可知Q为状态向量x和轨道不平顺导致的轮轨接触力的非线性函数，即Q＝F(x,u(t))，u(t)表示轨道不平顺。

将车辆轨道耦合模型的动力学方程式改写得：

$\left[\begin{array}{c}\stackrel{\&CenterDot;}{x}\\ \stackrel{\&CenterDot;\&CenterDot;}{x}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -{\&lsqb;M\&rsqb;}^{-1}\&lsqb;K\&rsqb;& -{\&lsqb;M\&rsqb;}^{-1}\&lsqb;C\&rsqb;\end{array}\right]\left[\begin{array}{c}x\\ \stackrel{\&CenterDot;}{x}\end{array}\right]+\left\{\begin{array}{c}0\\ {\&lsqb;M\&rsqb;}^{-1}\&lsqb;Q\&rsqb;\end{array}\right\}$

令：

$X=\left[\begin{array}{c}x\\ \stackrel{\&CenterDot;}{x}\end{array}\right],A=\left[\begin{array}{cc}0& 1\\ -{\&lsqb;M\&rsqb;}^{-1}\&lsqb;K\&rsqb;& -{\&lsqb;M\&rsqb;}^{-1}\&lsqb;C\&rsqb;\end{array}\right],f(x\left(t\right),u\left(t\right))=\left\{\begin{array}{c}0\\ {\&lsqb;M\&rsqb;}^{-1}\&lsqb;Q\&rsqb;\end{array}\right\}$

得车辆轨道耦合系统状态方程为：

$\stackrel{\&CenterDot;}{X}\left(t\right)=AX\left(t\right)+f(X\left(t\right),u\left(t\right))+w\left(t\right)---\left(2\right)$

式中u(t)是轨道不平顺，t表示时间，w(t)为模型噪声，A为常系数矩阵。

步骤3，根据步骤1中所述振动加速度传感器和陀螺仪采集到的信号，构建传感器观测方程，具体为

根据安装在车体上的振动加速度传感器获取车体的垂向振动加速度转向架构架上的振动加速度传感器获取构架的前构架垂向振动加速度和后构架垂向振动加速度车体上的陀螺仪获取的车体点头角速度转向架上的陀螺仪获取的前构架点头角速度和后构架点头角速度构建传感器观测方程：

$y\left(t\right)=\left[\begin{array}{c}{\stackrel{\&CenterDot;\&CenterDot;}{Z}}_c\\ {\stackrel{\&CenterDot;}{\mathrm{\&beta;}}}_c\\ {\stackrel{\&CenterDot;\&CenterDot;}{Z}}_{t1}\\ {\stackrel{\&CenterDot;}{\mathrm{\&beta;}}}_{t1}\\ {\stackrel{\&CenterDot;\&CenterDot;}{Z}}_{t2}\\ {\stackrel{\&CenterDot;}{\mathrm{\&beta;}}}_{t2}\end{array}\right]=\&lsqb;H\&rsqb;X\left(t\right)+\left[\begin{array}{c}g\\ 0\\ g\\ 0\\ g\\ 0\end{array}\right]+v\left(t\right)---\left(3\right)$

式(3)中，包含有额外的重力加速度g，视为测量噪声并非为零均值的，对检测到的加速度需先减去g，v(t)为测量噪声，H为量测矩阵，X(t)为t时刻状态向量及其导数的组合。

步骤4，配置滤波迭代方程，得到状态估计的时间更新方程和量测更新方程、估计误差协方差的时间更新方程和量测更新方程，具体为：

对车辆轨道耦合系统状态方程式(2)求取雅克比矩阵F_k-1：

$F_{k-1}=\frac{\&part;f_{k-1}}{\&part;X}{|}_{{\hat{x}}_{k-1}^{+}}=\frac{\&part;}{\&part;X}\&lsqb;AX\left(t\right)+f\left(X\left(t\right),u\left(t\right)\right)\&rsqb;{|}_{{\hat{x}}_{k-1}^{+}}=A+\frac{\&part;}{\&part;X}f\left(X\left(t\right),u\left(t\right)\right){|}_{{\hat{x}}_{k-1}^{+}}$

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

$A_u=\frac{\&part;}{\&part;X}f(X\left(t\right),u\left(t\right))$

按照车辆轨道耦合系统状态方程式(2)求取的对于状态的雅克比矩阵F_k-1整理为：

$F_{k-1}=(A+A_u){|}_{{\hat{x}}_{k-1}^{+}}$

状态方程式(2)求取的对于系统噪声w_k的雅克比矩阵L_k-1为单位阵：

$L_{k-1}=\frac{\&part;f_{k-1}}{\&part;w}{|}_{{\hat{x}}_{k-1}^{+}}=I$

估计误差协方差的时间更新方程如式(4)，状态估计的时间更新方程如式(5)：

$P_k^{-}=F_{k-1}{P}_{k-1}^{+}{F}_{k-1}^T+Q_{k-1}---\left(4\right)$

${\hat{x}}_k^{-}=f_{k-1}({\hat{x}}_{k-1}^{+},u_{k-1},0)---\left(5\right)$

式(4)中，P_k表示状态协方差，为k时刻先验协方差估计，F_k-1为k-1时刻的雅克比矩阵，为k-1时刻后验协方差估计，Q_k-1为k-1时刻的轮轨力，式(5)所示状态估计的时间更新采用4阶龙格-库塔积分实现，式(5)中，为状态的先验估计，u_k-1为k-1时刻的轮轨接触力，f_k-1为k-1时刻X(t)和u(t)的函数关系式；

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

$H_k=\frac{\&part;h_k}{\&part;x}{|}_{{\hat{x}}_{\stackrel{\&OverBar;}{k}}}=H,M_k=\frac{\&part;h_k}{\&part;v}{|}_{{\hat{x}}_{\stackrel{\&OverBar;}{k}}}=I$

其中，H_k、M_k分别是h_k关于x、v的偏微分，为k时刻状态的先验估计。

状态估计的量测更新如式(7)，估计误差协方差的更新如式(8)：

$K_k=P_k^{-}{H}_k^T{(H_k{P}_k^{-}{H}_k^T-R_k)}^{-1}---\left(6\right)$

${\hat{x}}_k^{+}={\hat{x}}_k^{-}+K_k(y_k-H_k{\hat{x}}_k^{-})---\left(7\right)$

$P_k^{+}=(I-K_k{H}_k)P_k^{-}---\left(8\right)$

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

步骤5，根据步骤4滤波迭代方程进行滤波迭代后，得到各时刻最优状态，利用状态估计的时间更新方程求逆后得到轨道不平顺估计。

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

$u_{k-1}=f_{k-1}^{-1}({\hat{x}}_{k-1},{\hat{x}}_k,0)$

其中，u_k-1为k-1时刻的轮轨接触力，是k-1时刻的状态估计，是k时刻的状态估计。

改写状态转移方程式如下：

$f(X_{k-1},u_{k-1})=\frac{X_k-X_{k-1}}{T}-{\mathrm{AX}}_{k-1}=\mathrm{\&tau;}$

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

$g-\frac{2}{M_w}{p}_1\left(t\right)={\mathrm{\&tau;}}_{\left(307\right)}$

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

$\begin{array}{c}{z}_{k-1}^{01}=-{\left(M_w{G}^{1.5}\frac{g+{\mathrm{\&tau;}}_{\left(307\right)}}{2}\right)}^{\frac{2}{3}}+Z_{w1}-Z_{r1}\\ =-{\left(M_w{G}^{1.5}\frac{g+{\mathrm{\&tau;}}_{\left(307\right)}}{2}\right)}^{\frac{2}{3}}+X_{k-1}^{\left(7\right)}+\underset{i=1}{\overset{NM}{\&Sigma;}}{Y}_i\left(x_{w1}\right)X_{k-1}^{(i+10)}\end{array}$

下面结合具体实施例对本发明作进一步详细说明。

本发明基于扩展卡尔曼滤波的轨道垂向不平顺估计系统，包括陀螺仪、振动加速度传感器和中央处理单元，其中陀螺仪安装于列车车体和走行部一系弹簧上方的构架上，振动加速度传感器安装于车体和走行部二系弹簧下方的构架上，所述振动加速度传感器和陀螺仪的输出端均接入中央处理单元，中央处理设置于车体内。所述陀螺仪的敏感轴与列车的轮对车轴平行。所述陀螺仪采用美国KVH公司的DSP-3000，中央处理单元采用北京阿尔泰科技的EPC-9351。

实施例1

通过建立车辆-轨道垂向耦合动力学模型，采集可观测的状态如车体加速度、车体角速度、构架加速度和构架角速度等，利用扩展卡尔曼滤波，结合观测方程、状态方程和由估计状态逆向计算轨道不平顺实现轨道不平顺的最优估计。

结合图2，建立翟婉明提出的(翟婉明著.车辆-轨道耦合动力学(第三版).科学出版社.2007)车辆-轨道垂向耦合动力学模型，包括车辆子模型、轨道子模型、轮轨作用力、Hertz非线性接触理论和新型预测-校正显式积分方法。车辆动力学子模型包括车体、前后构架的沉浮、俯仰运动，四个轮对的沉浮运动；轨道子模型采用连续支撑的欧拉梁力学模型，扣件、枕木和道床支撑视为连续支撑的弹簧阻尼模型；用Hertz非线性接触理论可以计算轮轨作用力；模型求解采用翟婉明提出的新型预测-校正显式积分方法，将所有的振动方程对时间步长进行积分，求取模型的动力学响应。对于车辆轨道垂向耦合模型，在仿真计算轨道长度l取70m的情况下，N＝100，根据上述文献，NM合理的取值为0.9N，因此本文NM＝90。表1为我国高速线路HST60基本动力参数，作为模型固定参数。

表1车辆-轨道垂向动力学模型参数

结合图3，该段实测轨道不平顺数据是以0.5m为间隔的离散采样序列，长度2km。图3绘制了该段实测数据的空间波形图和功率谱图。图3(a)为实测轨道不平顺空间波形图，可见该段数据包含的不平顺幅值小于5mm。图3(b)绘制了德国高干扰轨道谱线，可见该段数据包含高于德国高干谱的6-300m波长不平顺。由于采用的实测数据以0.5m为间隔，因此对数据进行三次样条插值处理，得到空间步长0.001m的轨道不平顺数据作为模型输入。仿真过程中，设定列车运行速度为20m/s，迭代空间步长0.001m，迭代时间间隔Δt＝0.05ms。

结合图4，利用MATLAB强大的矩阵运算能力，使得动力学仿真模型的编程实现相对于其他实现平台更有效而简洁。模型求解编程实现经过模型参数矩阵配置、初始化、轮轨力计算、钢轨振动力计算和迭代循环等过程。模型输出的车体振动加速度、车体点头角速度、构架振动加速度和构架点头角速度如图4所示，其中图4(a)为构架角速度波形图，图4(b)为构架振动加速度波形图，图4(c)为车体角速度波形图，图4(d)为车体振动加速度波形图，构架角速度的幅值响应范围为±5°/s，构架振动加速度幅值响应范围为±2m/s²，车体角速度的幅值响应范围为±0.2°/s，车体振动加速度幅值响应范围为±0.5m/s²。目前光纤陀螺仪分辨率<10°/h，可达±5°/s/h，完全能满足角速度检测要求，该方案具有工程可实施性。为模拟实际振动加速度计和陀螺仪观测过程中的观测噪声，加入方差为仿真响应最大幅值20％的高斯白噪声，其中加速度的单位为m/s²，角速度的单位为°/s。由于无法直接观测到过程信号x_k，过程激励噪声协方差的Q_k值通常较难确定，Q_k的意义是给过程信号“注入”不确定度，从而在获得可接受精度的前提下建立一个简单的状态模型。可见Q_k的选取是一个模型逼近的过程，车辆轨道耦合的动力学模型往往也存在一定的不确定度，实际现场往往需要通过实验的方式确定Q_k，由于是理论仿真，Q_k通过实验最终确定为0.05I_600×600。

结合图5，根据步骤2～3建立的车辆轨道耦合系统状态方程和传感器观测方程，以及步骤4配置的滤波迭代式，对仿真得到的惯性量进行最优滤波估计。首先将车辆轨道耦合模型的动力学方程改写为车辆轨道耦合系统状态方程，其次结合传感器观测方程，再后进行滤波迭代。滤波器迭代过程中，设定状态估计的积分步长为0.1ms，滤波器迭代步长1ms，空间步长0.1m。滤波后车辆系统状态量响应与动力学仿真输出的真实状态对比见图5，图5(a)为车轮加速度的波形图，图5(b)为车轮加速度的功率谱图，图5(c)为构架加速度的波形图，图5(d)为构架加速度的功率谱图，图5(e)为车体加速度的波形图，图5(f)为车体加速度的功率谱图，图5(g)为车轮速度的波形图，图5(h)为车轮速度的功率谱图，图5(i)为构架角速度的波形图，图5(j)构架角速度的功率谱图，图5(k)为车体角速度的波形图，图5(l)为车体角速度的功率谱图。可见有效的抑制了观测噪声带来的影响，得到了状态量可接受精度的最优估计。

对步骤四处理后的最优状态数据进行步骤五处理，得到最终的轨道不平顺估计。在步骤5中经过滤波迭代后得到的各时刻最优状态利用状态估计的时间更新方程式(5)求逆后即可得到轨道不平顺激励即

$u_{k-1}=f_{k-1}^{-1}({\hat{x}}_{k-1},{\hat{x}}_k,0)$

改写状态转移方程式如下：

$f(X_{k-1},u_{k-1})=\frac{X_k-X_{k-1}}{T}-{\mathrm{AX}}_{k-1}=\mathrm{\&tau;}$

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

$g-\frac{2}{M_w}{p}_1\left(t\right)={\mathrm{\&tau;}}_{\left(307\right)}$

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

$\begin{array}{c}{z}_{k-1}^{01}=-{\left(M_w{G}^{1.5}\frac{g+{\mathrm{\&tau;}}_{\left(307\right)}}{2}\right)}^{\frac{2}{3}}+Z_{w1}-Z_{r1}\\ =-{\left(M_w{G}^{1.5}\frac{g+{\mathrm{\&tau;}}_{\left(307\right)}}{2}\right)}^{\frac{2}{3}}+X_{k-1}^{\left(7\right)}+\underset{i=1}{\overset{NM}{\&Sigma;}}{Y}_i\left(x_{w1}\right)X_{k-1}^{(i+10)}\end{array}$

得到的轨道不平顺与真实轨道不平顺激励对比见图6，其中图6(a)为波形图，图6(b)为功率谱图。可见最优状态计算得到的轨道不平顺在各个波段均逼近真实值，尤其在2m-300m波长范围内功率谱值与真实不平顺表现了很强的一致性。

综上所述，本发明基于车体振动、构架振动、构架角速度和车速信号的检测，利用扩展卡尔曼滤波实现轨道垂向不平顺的最优估计，具有成本低、工程实施性好和在线实时监测的优点。

Claims (8)

1.一种基于扩展卡尔曼滤波的轨道垂向不平顺估计方法，其特征在于，包括以下步骤：

步骤1，在车体、转向架构架上安装振动加速度传感器，获得车体、构架的垂向振动加速度，在车体、构架上安装陀螺仪，获得车体点头角速度及构架点头角速度，以车速信号实现车体垂向振动、车体点头角速度、构架垂向振动和构架点头角速度的等空间采样；

步骤2，根据车辆轨道耦合模型的动力学方程构建车辆轨道耦合系统状态方程；

步骤3，根据步骤1中所述振动加速度传感器和陀螺仪采集到的信号，构建传感器观测方程；

步骤4，配置滤波迭代方程，得到状态估计的时间更新方程和量测更新方程、估计误差协方差的时间更新方程和量测更新方程；

步骤5，根据步骤4滤波迭代方程进行滤波迭代后，得到各时刻最优状态，利用状态估计的时间更新方程求逆后得到轨道不平顺估计。

2.根据权利要求1所述的轨道垂向不平顺估计方法，其特征在于，步骤2所述根据车辆轨道耦合模型的动力学方程构建车辆轨道耦合系统状态方程，具体为：

车辆轨道耦合模型的动力学方程为：

$M_z\stackrel{\&CenterDot;\&CenterDot;}{x}+C\stackrel{\&CenterDot;}{x}+Kx=Q$

式中，M_z为车辆的质量矩阵与轨道的质量矩阵结合、C为车辆的阻尼矩阵与轨道的阻尼矩阵结合、K为车辆的弹簧矩阵与轨道的弹簧矩阵结合；状态向量x为车辆与轨道动力学系统状态向量的组合：

x＝[Z_c β_c Z_t1 β_t1 Z_t2 β_t2 Z_w1 Z_w2 Z_w3 Z_w4 q_k(t) Z_sj Z_bi]^T(k＝1～NM,j＝1～N)

式中，Z_c、Z_t1、Z_t2、Z_w1、Z_w2、Z_w3、Z_w4、Z_sj、Z_bi依次分别为车体、前构架、后构架、第1～4轮对、第i轨枕、道床的浮沉运动，β_c、β_t1、β_t2依次分别为车体、前构架、后构架的点头运动，NM为钢轨模态阶数，N为轨枕支点总数；q_k(t)表示钢轨正则振型坐标，激励矩阵Q为车辆与轨道动力学系统激励的组合，由下式表示：

$\begin{array}{c}Q={\left[\begin{array}{ccccccccccccc}{M}_{c}g& 0& M_{t}g& 0& M_{t}g& 0& M_{w}g-2p_1\left(t\right)& M_{w}g-2p_2\left(t\right)& M_{w}g-2p_3\left(t\right)& M_{w}g-2p_4\left(t\right)& \underset{j=1}{\overset{4}{\&Sigma;}}{p}_j\left(t\right)Y_k\left(x_{wj}\right)& 0_{N\&times;1}& 0_{N\&times;1}\end{array}\right]}^T\\ (k=1~NM)\end{array}---\left(1\right)$

式中，M_c、M_t、M_w分别为车体、构架、轮对的质量，g为重力加速度，Y_k(x_wj)中Y_k为正则振型函数，x_wj为各车轮运动坐标，p₁(t)～p₄(t)为4个车轮的轮轨接触力；由式(1)知Q为状态向量x和轨道不平顺导致的轮轨接触力的非线性函数，即Q＝F(x,u(t))，u(t)表示轨道不平顺；

将车辆轨道耦合模型的动力学方程式改写得：

$\left[\begin{array}{c}\stackrel{\&CenterDot;}{x}\\ \stackrel{\&CenterDot;\&CenterDot;}{x}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -{\&lsqb;M\&rsqb;}^{-1}\&lsqb;K\&rsqb;& -{\&lsqb;M\&rsqb;}^{-1}\&lsqb;C\&rsqb;\end{array}\right]\left[\begin{array}{c}x\\ \stackrel{\&CenterDot;}{x}\end{array}\right]+\left\{\begin{array}{c}0\\ {\left[M\right]}^{-1}\&lsqb;Q\&rsqb;\end{array}\right\}$

令：

$X=\left[\begin{array}{c}x\\ \stackrel{\&CenterDot;}{x}\end{array}\right],A=\left[\begin{array}{cc}0& 1\\ -{\&lsqb;M\&rsqb;}^{-1}\&lsqb;K\&rsqb;& -{\&lsqb;M\&rsqb;}^{-1}\&lsqb;C\&rsqb;\end{array}\right],f\left(x\right(t),u\left(t\right))=\left\{\begin{array}{c}0\\ {\left[M\right]}^{-1}\&lsqb;Q\&rsqb;\end{array}\right\}$

得车辆轨道耦合系统状态方程为：

$\stackrel{\&CenterDot;}{X}\left(t\right)=AX\left(t\right)+f\left(X\right(t),u\left(t\right))+w\left(t\right)---\left(2\right)$

式中，u(t)是轨道不平顺，t表示时间，w(t)为模型噪声，A为常系数矩阵。

3.根据权利要求1所述的轨道垂向不平顺估计方法，其特征在于，步骤3所述构建传感器观测方程，具体为：

根据安装在车体上的振动加速度传感器获取车体的垂向振动加速度转向架构架上的振动加速度传感器获取构架的前构架垂向振动加速度和后构架垂向振动加速度车体上的陀螺仪获取的车体点头角速度转向架上的陀螺仪获取的前构架点头角速度和后构架点头角速度构建传感器观测方程：

$y\left(t\right)=\left[\begin{array}{c}{\stackrel{\&CenterDot;\&CenterDot;}{Z}}_c\\ {\stackrel{\&CenterDot;}{\mathrm{\&beta;}}}_c\\ {\stackrel{\&CenterDot;\&CenterDot;}{Z}}_{t1}\\ {\stackrel{\&CenterDot;}{\mathrm{\&beta;}}}_{t1}\\ {\stackrel{\&CenterDot;\&CenterDot;}{Z}}_{t2}\\ {\stackrel{\&CenterDot;}{\mathrm{\&beta;}}}_{t2}\end{array}\right]=\&lsqb;H\&rsqb;X\left(t\right)+\left[\begin{array}{c}g\\ 0\\ g\\ 0\\ g\\ 0\end{array}\right]+v\left(t\right)---\left(3\right)$

式(3)中，包含有额外的重力加速度g，视为测量噪声并非为零均值的，对检测到的加速度需先减去g，v(t)为测量噪声，H为量测矩阵，X(t)为t时刻状态向量及其导数的组合。

4.根据权利要求1所述的轨道垂向不平顺估计方法，其特征在于，步骤4所述配置滤波迭代方程，得到状态估计的时间更新方程和量测更新方程、估计误差协方差的时间更新方程和量测更新方程，具体为：

对车辆轨道耦合系统状态方程式(2)求取雅克比矩阵F_k-1：

$F_{k-1}=\frac{\&part;f_{k-1}}{\&part;X}{|}_{{\hat{x}}_{k-1}^{+}}=\frac{\&part;}{\&part;X}\&lsqb;AX\left(t\right)+f\left(X\right(t),u\left(t\right))\&rsqb;{|}_{{\hat{x}}_{k-1}^{+}}=A+\frac{\&part;}{\&part;X}f\left(X\right(t),u\left(t\right)){|}_{{\hat{x}}_{k-1}^{+}}$

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

$A_u=\frac{\&part;}{\&part;X}f\left(X\right(t),u\left(t\right))$

按照车辆轨道耦合系统状态方程式(2)求取的对于状态的雅克比矩阵F_k-1整理为：

$F_{k-1}=(A+A_u){|}_{{\hat{x}}_{k-1}^{+}}$

状态方程式(2)求取的对于系统噪声w_k的雅克比矩阵L_k-1为单位阵：

$L_{k-1}=\frac{\&part;f_{k-1}}{\&part;w}{|}_{{\hat{x}}_{k-1}^{+}}=I$

估计误差协方差的时间更新方程如式(4)，状态估计的时间更新方程如式(5)：

$P_k^{-}=F_{k-1}{P}_{k-1}^{+}{F}_{k-1}^T+Q_{k-1}---\left(4\right)$

${\hat{x}}_k^{-}=f_{k-1}({\hat{x}}_{k-1}^{+},u_{k-1},0)---\left(5\right)$

式(4)中，P_k表示状态协方差，为k时刻先验协方差估计，F_k-1为k-1时刻的雅克比矩阵，为k-1时刻后验协方差估计，Q_k-1为k-1时刻的轮轨力，式(5)所示状态估计的时间更新采用4阶龙格-库塔积分实现，式(5)中，为状态的先验估计，u_k-1为k-1时刻的轮轨接触力，f_k-1为k-1时刻X(t)和u(t)的函数关系式；

对观测方程式(3)求取雅克比矩阵：

$H_k=\frac{\&part;h_k}{\&part;x}{|}_{{\hat{x}}_k^{-}}=H,M_k=\frac{\&part;h_k}{\&part;v}{|}_{{\hat{x}}_k^{-}}=I$

其中，H_k、M_k分别是h_k关于x、v的偏微分，为k时刻状态的先验估计。

状态估计的量测更新如式(7)，估计误差协方差的更新如式(8)：

$K_k=P_k^{-}-H_k^T{(H_k{P}_k^{-}{H}_k^T-R_k)}^{-1}---\left(6\right)$

${\hat{x}}_k^{+}={\hat{x}}_k^{-}+K_k(y_k-H_k{\hat{x}}_k^{-})---\left(7\right)$

$P_k^{+}=(I-K_k{H}_k)P_k^{-}---\left(8\right)$

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

5.根据权利要求1所述的轨道垂向不平顺估计方法，其特征在于，步骤5所述根据步骤4滤波迭代方程进行滤波迭代后，得到各时刻最优状态，利用状态估计的时间更新方程求逆后得到轨道不平顺估计，其中：

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

$u_{k-1}=f_{k-1}^{-1}({\hat{x}}_{k-1},{\hat{x}}_k,0)$

其中，u_k-1为k-1时刻的轮轨接触力，是k-1时刻的状态估计，是k时刻的状态估计。

6.一种基于扩展卡尔曼滤波的轨道垂向不平顺估计系统，其特征在于，包括陀螺仪、振动加速度传感器和中央处理单元，其中陀螺仪安装于列车车体和走行部一系弹簧上方的构架上，振动加速度传感器安装于车体和走行部二系弹簧下方的构架上，所述振动加速度传感器和陀螺仪的输出端均接入中央处理单元，中央处理设置于车体内。

7.根据权利要求6所述的基于扩展卡尔曼滤波的轨道垂向不平顺估计系统，其特征在于，所述陀螺仪的敏感轴与列车的轮对车轴平行。

8.根据权利要求6所述的基于扩展卡尔曼滤波的轨道垂向不平顺估计系统，其特征在于，所述陀螺仪采用美国KVH公司的DSP-3000，中央处理单元采用北京阿尔泰科技的EPC-9351。