在轨道上运行的列车，主要利用轮轨黏着进行列车牵引或制动，即牵引力或制动力的本质是依靠轮轨黏着力发挥作用，因而轮轨黏着一直是牵引或制动技术领域的研究方向.

经过长期的试验研究分析，人们对轮轨黏着关系已经有了较深的认识，并通过线路试验或者理论公式计算黏着力^[1-3].但这些公式较复杂并且是基于一定的假设而不能兼顾所有的轮轨关系因素^[4]，所以无法直接运用于列车牵引计算、制动计算中.而在制动防滑控制中，尤其是在冰霜天气、低黏着力的条件下，充分利用轮轨黏着力十分重要，如何实时监测黏着力大小便成为防滑控制技术的难点之一.为解决这个问题，顾博川^[5]提出一种基于奇异值分解的改进强跟踪滤波算法用于机车黏着系数的估计；李宁洲等^[6]采用自适应子群协作QPSO(量子粒子群)算法估计黏着力，再利用智能模糊优化控制使机车保持最佳黏着力.但这两种算法都十分复杂，要调节的参数很多，不便应用于工程实践.因此，本文考虑用其他有效且更简单的算法对黏着力进行在线估测.

1 轮对动力学模型

根据黏着力估测的推导以及仿真验证要求，本文先建立轮对动力学模型.为研究分析主要问题和设计思路，将不涉及算法推导的仿真模型作一定的简化，只建立单轮对数学模型.图 1为单轮对受力分析图.

轮对动力学方程为

$ \left\{ \begin{array}{l} \dot v = - \frac{{{F_{{\rm{bw}}}}}}{{{m_{\rm{w}}}}}\\ \dot \omega = \frac{{\left( {{F_{{\rm{bw}}}} - {K_{{\varphi _{\rm{K}}}}}} \right){R_{\rm{w}}}}}{{{I_{\rm{w}}}}} \end{array} \right. $

(1)

式中：v为轮对纵向速度；ω为轮对的旋转角速度；K为电制动和空气制动作用在轮对的等效夹紧力; φ_K为闸瓦与车轮间的摩擦因数；F_bw为轮轨黏着力.模型考虑用轴重代替轮对质量，即可忽略架构对轮对的作用力.

车轮在钢轨上做滚动运动时，轮轨接触面上会产生微小滑动，不是纯滚动，是一种“滚中带滑”的状态，称之为蠕滑.也就是说在制动过程中，机车的速度v要大于车轮的线速度ωR_w.一般用滑移率s描述滑动的大小，如式(2)所示.

$ s = \frac{{v - \omega {R_{\rm{w}}}}}{v} $

(2)

轮轨黏着力的计算可按照Polach黏着力计算公式计算^[7]，如式(3)所示.

$ \left\{ \begin{array}{l} F = \frac{{2Qf}}{{\rm{ \mathsf{ π} }}}\left( {\frac{{{k_A}{\xi _1}}}{{1 + {{\left( {{k_A}{\xi _1}} \right)}^2}}} + \arctan \left( {{k_{\rm{s}}}{\xi _1}} \right)} \right)\\ {\xi _1} = \frac{2}{3}\frac{{C{\rm{ \mathsf{ π} }}{a^2}b}}{{Qf}}s\\ C = \frac{3}{8}\frac{G}{a}{C_{jj}} \end{array} \right. $

(3)

式中：Q为正压力；C_jj为Kalker系数；a、b为椭圆接触斑半轴长；G为刚性模量；f为轮轨间摩擦因数；k_A、k_s为调节参数.

2 在线估测计算方法

根据轮对动力学方程，可以直接测量或者间接计算的信号有等效夹紧力K和车轮角速度ω.本节利用K和ω设计多种算法在线估测轮轨黏着力F_bw.

2.1 微分相减法

观察轮对旋转动力学方程，如果能得到轮对角速度的微分信号$ \mathit{\dot{\omega }}$，那么黏着力可直接根据该方程代数运算得到.

对于一个有噪声的信号，直接微分容易扩大噪声.这里采用离散最速反馈系统^[8]跟踪角速度的微分信号$ \mathit{\dot{\omega }}$.

$ \left\{ \begin{array}{l} {F_{\rm{h}}} = {f_{\rm{h}}}\left( {\left( {\hat \omega \left( k \right) - \omega \left( k \right),\hat {\dot \omega} \left( k \right),r,h} \right)} \right.\\ \hat \omega \left( {k + 1} \right) = \hat \omega \left( k \right) + h\hat {\dot \omega} \left( k \right)\\ \hat {\dot \omega} \left( {k + 1} \right) = \hat {\dot \omega} \left( k \right) + h{F_{\rm{h}}} \end{array} \right. $

(4)

式中：f_h为fhan函数，在参考文献[8]中被称为最速综合函数；h为计算周期；r为可调节参数; 速度信号ω(k)为输入，输出为$ \mathit{\hat{\omega }}\left( \mathit{k} \right)、\mathit{\hat{\dot{\omega }}}\left( \mathit{k} \right)$，其中$ \mathit{\hat{\omega }}\left( \mathit{k} \right)$以r的加速度跟踪ω(k), $ \mathit{\hat{\dot{\omega }}}\left( \mathit{k} \right)$近似为ω(k)的微分.

得到角速度微分的估计值后，轮轨黏着力可由下式估测：

$ {{\hat F}_{{\rm{bw}}}} = K{\varphi _{\rm{K}}} + \frac{{{I_{\rm{w}}}\hat {\dot \omega} }}{{{R_{\rm{w}}}}} $

(5)

这种方法关键在于如何减低噪声对微分的影响，除了上述的离散最速反馈系统，还有很多滤波方法可以使用，算法自由，实现简单；不过通常实时性较差，一般应用在没有实时性要求的场合.

2.2 卡尔曼滤波法^[9]

根据轮对旋转动力学方程，其离散状态方程和观测方程分别为

$ \left\{ \begin{array}{l} \omega \left( {k + 1} \right) = \omega \left( k \right) - \frac{{h{R_{\rm{w}}}}}{{{I_{\rm{w}}}}}\left( {{\varphi _{\rm{K}}}K\left( k \right) - {F_{{\rm{bw}}}}\left( k \right)} \right)\\ y\left( {k + 1} \right) = \omega \left( k \right) + \chi \left( k \right) \end{array} \right. $

(6)

式中：y(k+1)为k时刻车轮角速度的测量值；χ(k)为观测噪声，均值为0，协方差为R.

黏着力估测值的迭代过程如下：

$ \begin{array}{l} \hat \omega \left( {k + 1\left| k \right.} \right) = \hat \omega \left( {k\left| k \right.} \right) - \\ \;\;\;\;\frac{{h{R_{\rm{w}}}}}{{{I_{\rm{w}}}}}\left( {{\varphi _{\rm{K}}}K\left( k \right) - {{\hat F}_{{\rm{bw}}}}\left( k \right)} \right) \end{array} $

$ {K_1}\left( {k + 1} \right) = P\left( {k\left| k \right.} \right){\left( {P\left( {k\left| k \right.} \right) + R} \right)^{ - 1}} $

$ \begin{array}{l} \hat \omega \left( {k + 1\left| {k + 1} \right.} \right) = \hat \omega \left( {k + 1\left| k \right.} \right) + \\ \;\;\;\;{K_1}\left( {k + 1} \right)\left( {y\left( {k + 1} \right) - \hat \omega \left( {k + 1\left| k \right.} \right)} \right) \end{array} $

$ P\left( {k + 1\left| {k + 1} \right.} \right) = \left( {I - {K_1}\left( {k + 1} \right)} \right)P\left( {k\left| k \right.} \right) $

$ {d_k} = \left( {1 - b} \right){\left( {1 - {b^{k + 1}}} \right)^{ - 1}} $

$ \begin{array}{l} {\varepsilon _1}\left( {k + 1} \right) = \hat \omega \left( {k + 1\left| {k + 1} \right.} \right) - \\ \;\;\;\;\;\hat \omega \left( {k\left| k \right.} \right) + \frac{{h{\varphi _{\rm{K}}}{R_{\rm{w}}}}}{{{I_{\rm{w}}}}}K\left( k \right) \end{array} $

$ {{\hat F}_{{\rm{bw}}}}\left( {k + 1} \right) = \left( {1 - {d_k}} \right){F_{{\rm{bw}}}}\left( k \right) + \frac{{{d_k}{I_{\rm{w}}}}}{{h{R_{\rm{w}}}}}{\varepsilon _1}\left( {k + 1} \right) $

(7)

式中：$ {\mathit{\hat{\omega }}}$(k+1|k)为角速度的预测值(先验值)；K₁(k+1)为卡尔曼增益；$ {\mathit{\hat{\omega }}}$(k+1|k+1)为角速度的估计值，初始值$ {\mathit{\hat{\omega }}}$(0|0)取角速度测量值y(0)效果最好；P(k+1|k+1)为估计值与真实值之间的协方差，初始值P(0|0)取10；b为遗忘因子；d_k为黏着力迭代权值；ε₁(k+1)为过程变量；$ {{{\mathit{\hat{F}}}}_{\rm{bw}}}$(k+1)为黏着力的估计值，初始值$ {{{\mathit{\hat{F}}}}_{\rm{bw}}}$(0)可取0.

卡尔曼滤波每一步迭代都极大利用上一步估计值和当前测量值，迭代效率极高；但该方法理论性强，且更强调对信号的滤波作用.

2.3 扩张状态观测器(ESO)法

将轮轨黏着力视为一个独立的状态变量，则轮对旋转运动状态方程可写作

$ \left\{ \begin{array}{l} \dot \omega = \frac{{{R_{\rm{w}}}}}{{{I_{\rm{w}}}}}{F_{{\rm{bw}}}} - \frac{{{\varphi _{\rm{K}}}{R_{\rm{w}}}}}{{{I_{\rm{w}}}}}K\\ {{\dot F}_{{\rm{bw}}}} = \zeta \end{array} \right. $

(8)

式中：ζ为黏着力的变化率.对式(8)描述的系统建立状态观测器为

$ \left\{ \begin{array}{l} \varepsilon = \hat \omega - \omega \\ \hat {\dot \omega} = \frac{{{R_{\rm{w}}}}}{{{I_{\rm{w}}}}}{{\hat F}_{{\rm{bw}}}} - \frac{{{\varphi _{\rm{K}}}{R_{\rm{w}}}}}{{{I_{\rm{w}}}}}K - {\beta _1}f\left( {\varepsilon ,{\alpha _1},{h_2}} \right)\\ {{\hat {\dot F}}_{{\rm{bw}}}} = - {\beta _2}f\left( {\varepsilon ,{\alpha _2},{h_2}} \right) \end{array} \right. $

(9)

其中，f(ε, α, h₂)一般取f_al函数，表达式为

$ {f_{{\rm{al}}}}\left( {x,\alpha ,\delta } \right) = \left\{ \begin{array}{l} \frac{x}{{{\delta ^{\left( {1 - \alpha } \right)}}}},\left| x \right| \le \delta \\ {\mathop{\rm sgn}} \left( x \right){\left| x \right|^\alpha },\left| x \right| > \delta \end{array} \right. $

(10)

式中:α和δ为函数的可调节参数，见参考文献[8].

当f(ε, α, h₂)取f_al函数时，上述观测器叫做非线性扩张状态观测器(NESO).当f(ε, α, h₂)不取f_al函数，而直接用ε代替时，观测器又叫做线性扩张状态观测器(LESO).

扩张状态观测器可以一次性分离出同一状态的不同变量——已知量和未知量(F_bw)，在引入特殊非线性函数提高滤波效果的同时，也保证了算法收敛时间.

2.4 负梯度迭代法^[10]

将轮对旋转运动微分方程进行拉式变换，再对两边同时进行一阶滤波处理，如下：

$ \frac{s}{{s + \lambda }}\omega \left( s \right) - \frac{1}{{s + \lambda }}\frac{{{R_{\rm{w}}}}}{{{I_{\rm{w}}}}}{F_{{\rm{bw}}}} = - \frac{1}{{s + \lambda }}\frac{{{\varphi _{\rm{K}}}{R_{\rm{w}}}}}{{{I_{\rm{w}}}}}K\left( s \right) $

(11)

式中：λ为一阶惯性环节的转折频率.再利用卷积定理，化简得

$ \begin{array}{l} \omega = \frac{{1 - {{\rm{e}}^{ - \lambda t}}}}{\lambda }\frac{{{R_{\rm{w}}}}}{{{I_{\rm{w}}}}}{F_{{\rm{bw}}}} - \frac{{{\varphi _{\rm{K}}}{R_{\rm{w}}}}}{{{I_{\rm{w}}}}}\int_0^t {{{\rm{e}}^{ - \lambda \left( {t - r} \right)}}K{\rm{d}}r} + \\ \;\;\;\;\;\lambda \int_0^t {{{\rm{e}}^{ - \lambda \left( {t - r} \right)}}\omega {\rm{d}}r} + {{\rm{e}}^{ - \lambda t}}\omega \left( 0 \right) \end{array} $

(12)

轮对角速度的估计误差ε为

$ \varepsilon = \hat \omega - \omega = \frac{{1 - {{\rm{e}}^{ - \lambda t}}}}{\lambda }\frac{{{R_{\rm{w}}}}}{{{I_{\rm{w}}}}}\left( {{{\hat F}_{{\rm{bw}}}} - {F_{{\rm{bw}}}}} \right) $

(13)

为使ε逼近0，F_bw的迭代方向应该是ε^Tε下降最快的方向(负梯度方向)，即

$ {{\hat {\dot \omega}}_{{\rm{bw}}}} = - p\frac{{\partial \left[ {{\varepsilon ^{\rm{T}}}\varepsilon } \right]}}{{\partial {{\hat F}_{{\rm{bw}}}}}} = - p\frac{{1 - {{\rm{e}}^{ - \lambda t}}}}{\lambda }\frac{{{R_{\rm{w}}}}}{{{I_{\rm{w}}}}}\varepsilon $

(14)

式中：p为负梯度迭代的可调增益.该式为黏着力估计值的微分方程，迭代时利用梯形公式转换成差分方程即可，初始值$ {{{\mathit{\hat{F}}}}_{\rm{bw}}}$(0)取0.

该方法利用拉式变换和卷积定理，将一阶滤波器融入到算法里，同时利用取负梯度方向的微分形式，使得该方法是几种方法中唯一一种收敛方式近似二阶系统的方法.不过该算法理论性强，推导复杂.

3 仿真对比分析 3.1 仿真参数设置

某型列车车辆参数如下：轴重m_w=4 797 kg；轮对转动惯量I_w=145 kg·m²；车轮半径R_w=0.42 m；闸瓦摩擦因数φ_K=0.3.本文借用Simulink仿真平台，根据上述参数搭建单轮对模型并编写算法，将轮对初始速度设为20 m·s^-1.

为更真实地模拟实际环境，测量的信号都加入噪声功率为10^-5的白噪声.另外某型列车采用MVB(multifunction vehicle bus)网络和CAN(controller area network)总线传输数据，最长传输时间为32 ms，本文用将e^-τs模拟测量信号的传输延迟，传输时间τ都设为0.05 s.

各算法的参数经过调试后，取值如表 1所示.

3.2 级位不变工况下的仿真

某型列车制动减速度为1.1 m·s^-2时，轮对等效夹紧力K为21 kN.如果轨面条件和制动级位都不变，车轮也没有滑行时，根据Polach黏着力公式，黏着力建立后将保持不变.此时各算法对黏着力的估测效果如图 2所示.

由图 2可知，算法对黏着力的估测都有一定的响应时间(第一次达到实际黏着力的时间)，同时在噪声的影响下，算法即便稳定收敛后也有一定的波动.现在以该工况下5种算法的响应时间和稳定后的最大误差来对比分析算法的估计效果，结果如表 2所示.

综合来看，NESO法的估计算法效果最佳.

3.3 黏着力变化时的仿真

实际运用由于轨面条件、车速等发生变化，黏着力也会变化，为反映算法能跟随黏着力的变化，现在模拟一种极端恶劣的情况——黏着力正(余)弦变化(频率为0.2 Hz).图 3为NESO算法对黏着力的变化跟踪效果.

从图 3可知，黏着力估测值相对实际值虽然有一定的幅值减弱和相位滞后，但整体还是可以在这种极端的情况下对黏着力变化进行估算的.

4 实例分析

本文以广佛地铁二期线路为例，用试验数据进行分析验证.广佛地铁二期线上车辆制动系统是EP2002架控系统，试验是在纯空气制动时进行的，测量的数据为制动缸压力和轮速，测量方式及测量仪器如图 4所示.

已知制动缸压力p_bc，可以根据以下公式计算得到轮对等效夹紧力：

$ K = n\eta \left( {\frac{{{\rm{ \mathsf{ π} }}{d^2}}}{4}{p_{{\rm{bc}}}} - {F_{{\rm{spring}}}}} \right) $

(15)

式中：n为制动盘数量，取2；η为夹钳效率，取0.97；d为制动缸活塞直径，取0.129 5 m；F_spring为制动缸复原弹簧恢复力，取1 450 N.

实际的轮轨黏着力利用公式直接计算得到.需要说明的是，此方法需要对信号进行高强度的滤波处理，实时性差，不适合在线估计，也不适合用其去判断其他估计算法的响应快慢.这里只用其去验证NESO估计算法对黏着力变化估测的准确性.

图 5是在不滑行工况下测得的制动缸压力和列车轮速变化曲线.

图 5中，0~2 s列车处于非制动状态，因不知列车牵引信息，这一段不做分析.2 s之后，制动施加，制动缸压力缓慢上升，列车开始减速.制动过程的黏着力变化估测如图 6所示.

图 6显示本文的NESO算法对黏着力的估测值十分接近黏着力非实时计算值，验证了NESO算法对黏着力变化估测的准确性.

5 结语

本文基于轮对动力学模型设计了微分相减、卡尔曼滤波、NESO、LESO和负梯度迭代等算法，对黏着力进行在线估测.仿真结果显示了算法的可行性，并表明NESO算法对黏着力变化的估测效果最好，并且通过实测数据验证了其准确性.

由于黏着力不易测量，本文设计的算法有利于列车牵引或制动时的黏着力控制.更重要的是，对于列车制动防滑控制，不仅需要防止车轮抱死，也希望充分利用黏着以减小制动距离，因此，实时估测黏着力将对轮轨作用关系、牵引和制动有很大帮助，所以本文的算法未来可用于列车防滑控制.