基于抗差自适应滤波的高速列车融合测速算法

摘要

针对Kalman滤波在高速列车融合测速过程中因观测粗差和动力学模型误差而引起的融合精度下降问题，提出一种基于抗差自适应滤波的高速列车融合测速算法。首先，在Kalman滤波的基础上构建异常检测函数和误差判别统计量，用于检测和区分传感器异常观测导致的观测粗差和动力学模型误差；然后，针对观测粗差和动力学模型误差，分别采用三段式函数和指数函数构造抗差因子和自适应因子，通过2种因子合理调节观测信息和模型信息在状态估计中的权重，从而降低观测粗差和动力学模型误差对融合结果的影响；最后，通过2种运行场景以及算法对比，仿真验证抗差自适应滤波算法性能。仿真结果表明，与基于Kalman滤波的融合测速算法相比，所提出算法无论在观测粗差场景还是在动力学模型误差场景，均具有更高的精度和稳定性。

关键词

铁路运输; 列车测速; 抗差自适应滤波; 高速列车; 信息融合

列车测速系统是列车自动防护（ATP）系统的重要组成部分。车载ATP系统一方面实时监控列车速度以防列车超速运行，另一方面根据列车速度生成速度‒距离曲线来控制列车运行间隔。因此，列车速度是保证ATP系统正常运转的基本参数，其测量精度和可靠性将直接影响列车运行安全和效率。

单一测速技术的精度和可靠性有限，难以满足高速列车运行需求，而多传感器融合测速技术可以结合多个传感器的优点，通过信息冗余和优势互补，为列车控制系统提供更加可靠、精确的列车速度信息^［1］。目前，对于多传感器融合测速，一般采用Kalman滤波实现信息融合^［2］，但Kalman滤波要求精确的动力学模型和观测噪声统计特性^［3］。在列车实际运行过程中，因各种干扰因素的存在（如车轮空转滑行、列车振动、外界天气条件影响）而导致的观测粗差和动力学模型误差问题常常严重影响滤波性能，造成测速精度下降。

为了解决观测粗差和动力学模型误差问题，国内外学者提出了多种改进方法。在观测粗差方面，Muniandi等^［4］通过在滤波估计之前设置验证门来预处理传感器数据，若观测值无法通过验证门则作为异常信息剔除，但其门限值与动力学模型有关，模型异常扰动可能导致正常观测数据被错误剔除。吴昕慧等^［5］采用H_∞滤波处理传感器观测数据，该算法不必对观测噪声进行任何假设，相比于Kalman滤波具有鲁棒性优势，但具有一定的适用范围，只能在约束条件下使用。蔡煊等^［6］在Kalman滤波的基础上引入抗差估计，通过降低传感器异常观测信息的占比有效抑制了观测粗差对融合结果的影响。在动力学模型误差方面，高学泽等^［7］提出了一种自适应交互式多模型算法，通过3种运动模型的切换实现对列车运动状态的跟踪，但是随着模型精细化程度的提高，运算量急剧增加。肖辰彬^［8］使用单一的当前统计模型跟踪列车运动状态，该模型对列车的加速度变化描述较为准确，但需要预先设置机动频率和加速度极值，若参数设置不准确，则模型跟踪精度下降甚至无法收敛。李增科等^［9］对滤波算法进行了改进，将自适应因子引入联邦Kalman滤波器，在保留原算法高计算效率的基础上实现了动力学模型误差的动态调节修正。

以上研究普遍只关注观测粗差或动力学模型误差的单因素影响，然而在实际运行场景中这2种误差往往都存在，因此若能从观测粗差和动力学模型误差两方面进行控制，将有效提升融合测速效果。针对Kalman滤波在高速列车融合测速中存在的问题，首先设计基于轮速传感器、多普勒雷达和加速度计的高速列车融合测速方案；然后提出一种基于抗差自适应滤波的融合测速算法，该算法将抗差滤波与自适应滤波相结合，利用抗差因子和自适应因子分别抑制观测粗差和动力学模型误差，从而降低传感器异常观测对融合结果的影响；最后仿真分析2种场景下的融合测速精度和稳定性，并与其他算法进行比较，验证了抗差自适应滤波算法的有效性。

1 高速列车融合测速方案

高速列车普遍采用轮速传感器和多普勒雷达相结合的方式测量列车速度。轮速传感器通过测量车轮转速准确计算列车运行速度，然而轮对一旦发生空转或滑行，测量结果将出现严重偏差^［10］。多普勒雷达不受车轮状态影响，在列车发生空转或滑行时仍能保持测量精度，但多普勒雷达对应用环境和外界天气条件比较敏感^［11］，容易产生异常观测。如果不对传感器数据进行处理而直接用于信息融合，就会严重影响融合精度，因此要在信息融合之前增加滤波环节。此外，考虑到滤波器动力学模型含有加速度信息，故引入加速度计辅助滤波器工作。设计的融合测速方案如图1所示。

图1  融合测速方案示意图

Fig.1  Schematic diagram of fusion speed measurement system scheme

测速传感器和加速度计分别采集列车速度信息和加速度信息，采集的信息经处理后得到当前时刻的列车速度和加速度，并输入对应的局部滤波器。局部滤波器在加速度信息辅助下对测速数据进行滤波处理，得到局部估计速度，将2个滤波器的局部估计速度送入全局信息融合模块进行融合，最后输出列车融合速度。

2 传感器误差特性分析

轮速传感器通过测量车轮转速计算列车运行速度，当列车正常运行时车轮转速与车体走行速度基本相同，此时测速误差主要取决于轮速传感器的自身精度。由传感器自身精度造成的误差被视为正常观测噪声。当车轮发生空转或滑行时，车轮转速与车体实际速度偏差较大，此时轮速传感器测量误差被视为异常观测误差。

多普勒雷达属于非接触测量方式，不受车轮空转或滑行的影响，在良好工作条件下测速误差主要来源于传感器正常观测噪声。多普勒雷达对应用环境比较敏感，列车纵向振动、外界天气条件、轨面状况等因素均会导致异常观测误差。

加速度计测量的是水平运动方向上的加速度，当列车处于坡道时将出现测量误差^［12］。此外，列车受到的冲击和振动可能使加速度测量值发生突变。以上2种误差均被视为异常观测误差。

基于上述3种传感器的误差特性分析，可将传感器误差归为2类：正常观测噪声和异常观测误差。正常观测噪声是传感器自身固有的测量误差，误差小且一直存在，一般可以认为是零均值的高斯白噪声。异常观测误差只有在特定情况下才会出现，误差大且没有准确分布特性。在滤波器工作过程中，异常观测误差不仅严重影响当前时刻滤波结果，还对后续时刻造成影响，所以必须对传感器的异常观测误差加以控制。

3 基于抗差自适应滤波的融合测速 算法

Kalman滤波是建立在状态空间模型和白噪声条件下的最优估计方法，对于传感器正常观测噪声，Kalman滤波具有较好的抑制效果。然而，在列车实际运行过程中，由于传感器异常观测问题的存在，观测粗差和动力学模型误差是不可避免的。因此，对Kalman滤波进行了抗差自适应改进，以达到抵抗传感器异常观测的目的。

3.1　Kalman滤波算法

Kalman滤波的状态空间模型由状态方程和观测方程组成。状态方程主要反映载体的运动规律，也被称为动力学模型。以列车运动时的速度$v$和加速度$a$作为状态变量，构建的状态方程为

$${\left[\begin{array}{l}v\\ a\end{array}\right]}_k=\left[\begin{array}{cc}\mathrm{1}& t\\ \mathrm{0}& \mathrm{1}\end{array}\right]{\left[\begin{array}{l}v\\ a\end{array}\right]}_{k-\mathrm{1}}+\left[\begin{array}{c}\mathrm{0}\\ \mathrm{1}\end{array}\right]\mathrm{\Delta }{a}_k+{\left[\begin{array}{l}{w}_v\\ w_a\end{array}\right]}_k$$

（1）

即：

$${\mathit{X}}_k=\mathit{\Phi }{\mathit{X}}_{k-\mathrm{1}}+\mathit{B}\mathrm{\Delta }{a}_k+{\mathit{W}}_k$$

（2）

式中：$k$为当前时刻；$t$为$(k-\mathrm{1})$至$k$时刻的时间步长；${\mathit{X}}_k$为状态向量；$\mathit{\Phi }$为状态转移矩阵；$\mathit{B}$为输入控制矩阵；$\mathrm{\Delta }{a}_k$为加速度计$k$与$(k-\mathrm{1})$时刻的测量值之差，通过引入加速度信息来描述$(k-\mathrm{1})$至$k$时刻的外力作用；${\mathit{W}}_k={\left[w_v{w}_a\right]}^{\mathrm{T}}$为系统噪声向量，其中$w_v$、$w_a$分别为影响速度和加速度的系统噪声，假设${\mathit{W}}_k$是均值为零、方差为${\mathit{Q}}_k$的白噪声序列。

以测速传感器的观测值为滤波对象，测速传感器的观测方程表示为

$${\mathit{Z}}_k=\mathit{H}{\mathit{X}}_k+{\mathit{V}}_k=\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\end{array}\right]{\left[\begin{array}{c}v\\ a\end{array}\right]}_k+{\mathit{V}}_k$$

（3）

式中：${\mathit{Z}}_k$为测速传感器$k$时刻观测值；$\mathit{H}$为观测矩阵；${\mathit{V}}_k$为测速传感器观测噪声，假设其是均值为零、方差为${\mathit{R}}_k$的白噪声序列。多普勒雷达的观测方程与式（3）相似，只是噪声方差取值不同。

在建立式（2）的状态方程和式（3）的观测方程之后，即可对测量数据进行滤波处理。假设$k$时刻状态最优估计值为${\widehat{\mathit{X}}}_k$，估计误差协方差矩阵为${\mathit{P}}_k$，则Kalman滤波工作过程如下^［13］：

$${\widehat{\mathit{X}}}_{k,k-\mathrm{1}}=\mathit{\Phi }{\widehat{\mathit{X}}}_{k-\mathrm{1}}+\mathit{B}\mathrm{\Delta }{a}_k$$

（4）

$${\mathit{P}}_{k,k-\mathrm{1}}=\mathit{\Phi }{\mathit{P}}_{k-\mathrm{1}}{\mathit{\Phi }}^{\mathrm{T}}+{\mathit{Q}}_k$$

（5）

$${\mathit{K}}_k={\mathit{P}}_{k,k-\mathrm{1}}{\mathit{H}}^{\mathrm{T}}(\mathit{H}{\mathit{P}}_{k,k-\mathrm{1}}{\mathit{H}}^{\mathrm{T}}+{\mathit{R}}_k{)}^{-\mathrm{1}}$$

（6）

$${\widehat{\mathit{X}}}_k={\widehat{\mathit{X}}}_{k,k-\mathrm{1}}+{\mathit{K}}_k({\mathit{Z}}_k-\mathit{H}{\widehat{\mathit{X}}}_{k,k-\mathrm{1}})$$

（7）

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

（8）

式中：${\widehat{\mathit{X}}}_{k,k-\mathrm{1}}$为状态预测值；${\mathit{P}}_{k,k-\mathrm{1}}$为预测误差协方差矩阵；${\mathit{K}}_k$为滤波增益矩阵；$\mathit{I}$为单位矩阵。由式（4）~（8）可以看出，Kalman滤波工作过程就是根据前一时刻的估计值和当前时刻的观测值不断递推处理的过程。

3.2　抗差自适应滤波

在列车运行过程中，测速传感器的异常观测导致观测信息${\mathit{Z}}_k$出现粗差现象，加速度计的异常观测导致模型预测信息${\widehat{\mathit{X}}}_{k,k-\mathrm{1}}$出现偏差。由式（7）可知，以上2种误差将对状态估计${\widehat{\mathit{X}}}_k$产生严重干扰。针对上述问题，将抗差因子和自适应因子引入Kalman滤波，增强算法抵抗观测粗差和动力学模型误差的能力。

首先需要判断观测信息或模型信息是否出现异常。由于预测残差能够较好地反映观测异常和模型异常情况，因此基于预测残差进行判断。预测残差及其协方差矩阵分别为：

$${\mathit{\epsilon }}_k={\mathit{Z}}_k-\mathit{H}{\widehat{\mathit{X}}}_{k,k-\mathrm{1}}$$

（9）

$$\mathrm{E}\{{\mathit{\epsilon }}_k{\mathit{\epsilon }}_k^{\mathrm{T}}\}=\mathit{H}{\mathit{P}}_{k,k-\mathrm{1}}{\mathit{H}}^{\mathrm{T}}+{\mathit{R}}_k$$

（10）

在传感器均无异常观测的情况下，预测残差服从均值为零的高斯分布，即：

$${\mathit{\epsilon }}_k~N(\mathrm{0},\mathit{H}{\mathit{P}}_{k,k-\mathrm{1}}{\mathit{H}}^{\mathrm{T}}+{\mathit{R}}_k)$$

（11）

基于预测残差定义异常检测函数为

$${\lambda }_k={\mathit{\epsilon }}_k^{\mathrm{T}}(\mathit{H}{\mathit{P}}_{k,k-\mathrm{1}}{\mathit{H}}^{\mathrm{T}}+{\mathit{R}}_k{)}^{-\mathrm{1}}{\mathit{\epsilon }}_k$$

（12）

根据式（12），若${\mathit{\epsilon }}_k$服从高斯分布，则${\lambda }_k$服从自由度为$m$的${\chi }^{\mathrm{2}}$分布，$m$为观测维度。当观测信息或模型信息异常时，异常检测函数不再服从上述分布，因此可以通过${\chi }^{\mathrm{2}}$分布的上分位点性质进行检验。例如，在$m=\mathrm{1}$且误检率为10^-3时，对应检测门限$T_{\mathrm{D}}$=10.83，若${\lambda }_k>T_{\mathrm{D}}$，则判定有异常，否则认为无异常。然而，${\lambda }_k$无法具体区分是观测信息异常还是模型信息异常，因此对两者分别构造统计量：

$$\left\{\begin{array}{l}{\mathit{\epsilon }}_{\mathrm{1}}={\mathit{Z}}_k-\mathit{H}{\widehat{\mathit{X}}}_{k-\mathrm{1}}\\ {\mathit{\epsilon }}_{\mathrm{2}}=\mathit{H}({\widehat{\mathit{X}}}_{k,k-\mathrm{1}}-{\widehat{\mathit{X}}}_{k-\mathrm{1}})\end{array}\right.$$

（13）

$$\left\{\begin{array}{l}{t}_{\mathrm{1}}={\mathit{\epsilon }}_{\mathrm{1}}^{\mathrm{T}}(\mathit{H}{\mathit{P}}_{k-\mathrm{1}}{\mathit{H}}^{\mathrm{T}}+{\mathit{R}}_k{)}^{-\mathrm{1}}{\mathit{\epsilon }}_{\mathrm{1}}\\ t_{\mathrm{2}}={\mathit{\epsilon }}_{\mathrm{2}}^{\mathrm{T}}(\mathit{H}({\mathit{P}}_{k-\mathrm{1}}+{\mathit{P}}_{k,k-\mathrm{1}}){\mathit{H}}^{\mathrm{T}}{)}^{-\mathrm{1}}{\mathit{\epsilon }}_{\mathrm{2}}\end{array}\right.$$

（14）

式中：$t_{\mathrm{1}}$为观测判别统计量；$t_{\mathrm{2}}$为模型判别统计量。当判定系统出现异常后，分别计算$t_{\mathrm{1}}$和$t_{\mathrm{2}}$，如果是观测粗差引发的异常，就会出现$t_{\mathrm{1}}\ge t_{\mathrm{2}}$，此时应当进行抗差滤波；如果是动力学模型误差引发的异常，就会出现$t_{\mathrm{1}}<t_{\mathrm{2}}$，此时应当进行自适应滤波。

抗差滤波通过抗差因子放大观测噪声方差以调节增益矩阵${\mathit{K}}_k$的值，从而降低观测信息对状态估计的权重，即将式（6）改为

$${\mathit{K}}_k={\mathit{P}}_{k,k-\mathrm{1}}{\mathit{H}}^{\mathrm{T}}{\left(\mathit{H}{\mathit{P}}_{k,k-\mathrm{1}}{\mathit{H}}^{\mathrm{T}}+\frac{\mathrm{1}}{r_k}{\mathit{R}}_k\right)}^{-\mathrm{1}}$$

（15）

式中，$r_k$为抗差因子。抗差因子构造方法主要有Huber法、丹麦法和三段式函数法^［14-15］。选用三段式函数法构造抗差因子，其函数表达式为：

$$r_k=\left\{\begin{array}{l}\mathrm{1},\left|{\overline{\mathit{\epsilon }}}_k\right|\le k_{\mathrm{0}}\\ \frac{k_{\mathrm{0}}}{\left|{\overline{\mathit{\epsilon }}}_k\right|}{\left(\frac{k_{\mathrm{1}}-\left|{\overline{\mathit{\epsilon }}}_k\right|}{k_{\mathrm{1}}-k_{\mathrm{0}}}\right)}^{\mathrm{2}},k_{\mathrm{0}}<\left|{\overline{\mathit{\epsilon }}}_k\right|\le k_{\mathrm{1}}\\ {\mathrm{10}}^{-\mathrm{6}},\left|{\overline{\mathit{\epsilon }}}_k\right|>k_{\mathrm{1}}\end{array}\right.$$

（16）

式中：${\overline{\mathit{\epsilon }}}_k$为标准化预测残差；$k_{\mathrm{0}}$、$k_{\mathrm{1}}$为常数，通常$k_{\mathrm{0}}$取1.0~2.5，$k_{\mathrm{1}}$取3.5~8.0。当观测信息出现粗差时，$r_k$取值减小，${\mathit{K}}_k$也减小，使得观测信息在状态估计中的占比降低，而模型信息占比增大；当观测粗差比较大时，$r_k$取值几乎为零，此时观测信息不参与状态估计，实现了粗差淘汰。

自适应滤波通过自适应因子放大预测误差协方差矩阵来降低动力学模型误差对状态估计的影响。自适应滤波的增益矩阵计算方法为

$${\mathit{K}}_k=\frac{\mathrm{1}}{{\alpha }_k}{\mathit{P}}_{k,k-\mathrm{1}}{\mathit{H}}^{\mathrm{T}}{\left(\frac{\mathrm{1}}{{\alpha }_k}\mathit{H}{\mathit{P}}_{k,k-\mathrm{1}}{\mathit{H}}^{\mathrm{T}}+{\mathit{R}}_k\right)}^{-\mathrm{1}}$$

（17）

式中，${\alpha }_k$为自适应因子。基于指数函数构造自适应因子的两段式函数模型^［16-17］：

$${\alpha }_k=\left\{\begin{array}{l}\mathrm{1},\mathrm{\Delta }{\mathit{\epsilon }}_k\le c\\ {\mathrm{e}}^{-(\mathrm{\Delta }{\mathit{\epsilon }}_k{-c)}^{\mathrm{2}}},\mathrm{\Delta }{\mathit{\epsilon }}_k>c\end{array}\right.$$

（18）

$$\mathrm{\Delta }{\mathit{\epsilon }}_k={\left(\frac{{\mathit{\epsilon }}_k^{\mathrm{T}}{\mathit{\epsilon }}_k}{\mathrm{t}\mathrm{r}(\mathit{H}{\mathit{P}}_{k,k-\mathrm{1}}{\mathit{H}}^{\mathrm{T}}+{\mathit{R}}_k)}\right)}^{\mathrm{1}/\mathrm{2}}$$

（19）

式中：$c$为常数，取值为1.0~3.0；$\mathrm{t}\mathrm{r}$为矩阵的迹。当动力学模型出现误差时，${\alpha }_k$取值减小，${\mathit{K}}_k$增大，使得模型预测信息在状态估计中的占比降低。

综上所述，本文抗差自适应滤波算法相比于传统Kalman滤波算法，不仅能控制观测信息的异常误差影响，还对动力学模型误差也有一定的抑制能力，其算法流程如图2所示。

图2  抗差自适应滤波算法流程

Fig.2  Flow chart of robust adaptive filter algorithm

滤波开始时，首先获取当前时刻的速度和加速度测量信息，然后根据加速度测量值计算模型预测信息${\widehat{\mathit{X}}}_{k,k-\mathrm{1}}$，并与速度观测信息${\mathit{Z}}_k$作差得到预测残差${\mathit{\epsilon }}_k$。根据${\mathit{\epsilon }}_k$计算异常检测函数${\lambda }_k$，若${\lambda }_k\le T_{\mathrm{D}}$，则认为没有异常，此时选择Kalman滤波；若${\lambda }_k>T_{\mathrm{D}}$，则判定有异常，然后计算观测判别统计量$t_{\mathrm{1}}$和模型判别统计量$t_{\mathrm{2}}$，比较$t_{\mathrm{1}}$和$t_{\mathrm{2}}$的值，由此选择抗差滤波或自适应滤波。最终得到状态估计${\widehat{\mathit{X}}}_k$和误差协方差矩阵${\mathit{P}}_k$。

3.3　全局信息融合方法

将2个滤波器的局部估计信息输入全局信息融合模块进行融合，得到全局估计信息。2个滤波器的估计信息分别表示为${\widehat{\mathit{X}}}_{\mathrm{1},k}$和${\mathit{P}}_{\mathrm{1},k}$、${\widehat{\mathit{X}}}_{\mathrm{2},k}$和${\mathit{P}}_{\mathrm{2},k}$。全局估计${\widehat{\mathit{X}}}_{\mathrm{g}}$为局部估计的线性组合，即：

$${\widehat{\mathit{X}}}_{\mathrm{g}}={\mathit{W}}_{\mathrm{1}}{\widehat{\mathit{X}}}_{\mathrm{1},k}+{\mathit{W}}_{\mathrm{2}}{\widehat{\mathit{X}}}_{\mathrm{2},k}$$

（20）

式中，${\mathit{W}}_{\mathrm{1}}$和${\mathit{W}}_{\mathrm{2}}$为权重矩阵。若全局估计${\widehat{\mathit{X}}}_{\mathrm{g}}$是最优无偏的，则${\mathit{W}}_{\mathrm{1}}$和${\mathit{W}}_{\mathrm{2}}$满足

$$\left\{\begin{array}{c}{\mathit{W}}_{\mathrm{1}}=({\mathit{P}}_{\mathrm{1},k}^{-\mathrm{1}}+{\mathit{P}}_{\mathrm{2},k}^{-\mathrm{1}}{)}^{-\mathrm{1}}{\mathit{P}}_{\mathrm{1},k}^{-\mathrm{1}}\\ {\mathit{W}}_{\mathrm{2}}=({\mathit{P}}_{\mathrm{1},k}^{-\mathrm{1}}+{\mathit{P}}_{\mathrm{2},k}^{-\mathrm{1}}{)}^{-\mathrm{1}}{\mathit{P}}_{\mathrm{2},k}^{-\mathrm{1}}\end{array}\right.$$

（21）

将式（21）代入式（20），可得

$${\widehat{\mathit{X}}}_{\mathrm{g}}=({\mathit{P}}_{\mathrm{1},k}^{-\mathrm{1}}+{\mathit{P}}_{\mathrm{2},k}^{-\mathrm{1}}{)}^{-\mathrm{1}}({\mathit{P}}_{\mathrm{1},k}^{-\mathrm{1}}{\widehat{\mathit{X}}}_{\mathrm{1},k}+{\mathit{P}}_{\mathrm{2},k}^{-\mathrm{1}}{\widehat{\mathit{X}}}_{\mathrm{2},k})$$

（22）

全局信息融合模块按式（22）对滤波器的输出信息进行处理后，可以得到列车状态的全局最优估计值，其速度分量即为列车融合速度。

4 仿真验证与分析

针对本文提出的基于抗差自适应滤波的列车融合测速算法，通过2种仿真场景以及与其他算法的对比，验证其抗差能力和自适应能力。

4.1　仿真参数设置

以CRH380A型动车组为仿真车型，车重取定员时的429 t，回转质量系数取0.08，在不考虑附加阻力的情况下，通过列车牵引计算生成列车运行数据。仿真产生的列车速度如图3所示，列车运行时间共计800 s，运行过程包括加速、匀速、惰行和减速等阶段。传感器参数设置为：轮速传感器输出频率20 Hz，速度测量误差0.2 m·s^-1；多普勒雷达输出频率2 Hz，速度测量误差0.2 m·s^-1；加速度计输出频率100 Hz，加速度测量误差0.03 m·s^-2。信息融合周期设为0.5 s。

图3  CRH380A型动车组速度仿真曲线

Fig.3  Speed simulation curve of CRH380A train

仿真设置2种运行场景：场景1，在正常观测噪声的基础上，对轮速传感器和多普勒雷达加入异常观测数据（模拟观测粗差）；场景2，对加速度计加入异常观测数据（模拟动力学模型误差）。在3.3节全局信息融合方法的基础上，在局部滤波器中分别采用Kalman滤波（KF）算法、文献［18］的抗差无迹Kalman滤波（抗差UKF）算法以及本文抗差自适应滤波算法进行融合测速仿真对比。其中，抗差因子常数$k_{\mathrm{0}}$、$k_{\mathrm{1}}$设为2.5、4.0，自适应因子常数c设为2.0。

4.2　仿真验证

4.2.1　抗差性能验证

在场景1中，对本文算法的抗差性能进行验证。根据工作特点对测速传感器加入异常观测数据，在列车加速/减速阶段对轮速传感器加入空转/滑行误差，每次空转/滑行持续时间15 s，对多普勒雷达每隔100 s加入突变测量误差，用于模拟列车振动和运行环境影响。图4为轮速传感器和多普勒雷达相对于列车速度仿真数据的测量误差。

图4  测速传感器速度测量误差

Fig.4  Measurement errors of speed sensors

图5给出了局部滤波过程中抗差因子的变化。可以看出，2个局部滤波器的抗差因子都能准确识别异常观测数据，并且大幅降低自身数值，使得异常观测数据无法参与滤波计算，这是滤波算法能够抵御观测粗差的主要原因。

图5  滤波过程中抗差因子变化

Fig.5  Change of robust factor during filtering

以仿真速度为参考值，场景1融合速度误差如图6所示。可以看出，由于不具备抗差能力，因此KF算法在空转和滑行期间出现了大幅误差；抗差UKF算法能够抑制大部分的观测粗差，但在部分时刻仍然出现了较大误差，原因是其采用了基于验后残差的抗差因子，计算存在滞后性，影响了抗差能力的发挥；本文算法使用的是预测残差，不存在滞后问题，因此在整个融合过程中都能稳定工作。

图6  场景1融合速度误差

Fig.6  Fusion speed errors of scene 1

为了进一步比较算法性能，采用平均绝对误差（e_MAE）衡量精度，采用均方根误差（e_RMSE）衡量稳定性。表1为3种算法在场景1下的平均绝对误差和均方根误差统计结果。在精度方面，本文算法的测速精度为0.038 4 m·s^-1，相比于KF和抗差UKF算法分别提升了52%和10%。在稳定性方面，本文算法的均方根误差最小，说明其稳定性最强，相比于KF和抗差UKF算法分别提升了72%和13%。因此，本文算法相比于其他2种算法，在观测粗差环境下的精度和稳定性均有明显提升，其中在稳定性方面的提升更为明显。

表1  3种算法在场景1下的平均绝对误差和均方根误差

Tab.1  e_MAE and e_RMSE of three algorithms in scene 1

算法	eMAE/（m·s-1）	eRMSE/（m·s-1）
KF算法	0.079 2	0.171 2
抗差UKF算法	0.042 6	0.055 8
本文算法	0.038 4	0.048 5

4.2.2　自适应性能验证

在场景2中，验证本文算法对动力学模型误差的自适应能力。在部分时刻对2个加速度计分别加入突变测量误差，用于模拟动力学模型误差。图7为2个加速度计相对于列车加速度仿真数据的测量误差。

图7  加速度计加速度测量误差

Fig.7  Measurement errors of accelerometers

图8给出了局部滤波过程中自适应因子变化。可以看出，2个局部滤波器的自适应因子能够准确识别模型异常信息并且大幅降低模型信息在状态估计中的权重。

图 8  滤波过程中自适应因子变化

Fig.8  Change of adaptive factor during filtering

图9为不同算法场景2融合速度误差对比。可以看出，由于不具备自适应能力，KF算法和抗差UKF算法均出现了较大误差，而本文算法则有效控制了动力学模型误差对融合结果的影响。此外，抗差UKF算法在部分时刻的误差大于KF算法，原因是动力学模型误差导致残差变大，抗差因子进一步增大模型信息权重，使得滤波结果更加不可靠。

图 9  场景2融合速度误差

Fig.9  Fusion speed errors of scene 2

表2为3种算法在场景2下的平均绝对误差和均方根误差统计结果。在精度方面，本文算法的测速精度为0.039 3 m·s^-1，相比于KF和抗差UKF算法分别提升了18%和25%。在稳定性方面，本文算法的均方根误差最小，稳定性最强，相比于KF和抗差UKF算法分别提升了30%和43%。因此，本文算法相比于其他2种算法在模型异常扰动环境下的精度和稳定性均有明显提升，在稳定性方面的提升更为明显。

表2  3种算法在场景2下的平均绝对误差和均方根误差

Tab.2  e_MAE and e_RMSE of three algorithms in scene 2

算法	eMAE/（m·s-1）	eRMSE/（m·s-1）
KF算法	0.048 1	0.084 2
抗差UKF算法	0.052 6	0.104 3
本文算法	0.039 3	0.059 1

综合2种场景下的仿真结果可以得出：在有观测粗差或动力学模型误差影响的情况下，本文算法的抗差因子和自适应因子均能准确识别异常信息并降低异常信息对融合结果的干扰，最终的融合精度分别达到了0.038 4 m·s^-1和0.039 3 m·s^-1 ，说明本文算法具有较强的抗差能力和自适应能力；与KF算法相比，本文算法在2种场景下的融合精度分别提升了52%和18%，在抗差性能上的改进更为明显；与抗差UKF算法相比，本文算法在2种场景下的融合精度分别提升了10%和25%，在自适应性能上的改进更为明显。

5 结语

在高速列车融合测速过程中，由于传感器异常观测问题的存在，采用传统Kalman滤波往往容错性能不佳，导致融合精度下降。针对Kalman滤波存在的缺陷，提出一种基于抗差自适应滤波的列车融合测速算法。该算法将抗差滤波与自适应滤波相结合，利用抗差因子和自适应因子分别解决由传感器异常观测导致的观测粗差和动力学模型误差问题。仿真结果表明，该方法能有效改善观测粗差和动力学模型误差对融合结果的影响，提高了列车融合测速的精度和稳定性。

作者贡献声明

王小敏：总体构想，算法设计，论文撰写。

贾钰林：算法仿真，论文撰写。

张亚东：算法设计，论文修改。

魏维伟：测速方案修改。

何 静：实验分析。