TIME DOMAIN IDENTIFICATION AND COMPARISON OF VERTICAL WHEEL-RAIL FORCE OF RAIL VEHICLES AND ITS MACHINE LEARNING CORRECTION
-
摘要: 为了尽可能减小轨道车辆垂向轮轨力时域识别中存在的误差, 以时域法为基础, 开展了基于机器学习修正的轨道车辆垂向轮轨力识别研究. 首先建立了车辆动力学仿真模型, 获取了车辆在随机轨道激励下以250 km/h速度行驶时的轴箱加速度响应和垂向轮轨力. 其次, 建立了Green函数法、状态空间法这2种时域法对应的动载荷识别模型, 对状态空间法的初值误差进行了分析, 并引入多项式拟合法修正其趋势项误差, 进而对比分析了2种方法的计算精度和计算效率. 然后, 针对时域法存在的识别误差, 提出采用NARX (nonlinear autoregressive models with exogenous inputs)模型对识别误差进行训练和预测, 用于消减模型中存在的如响应观测不全与观测噪声等因素造成的影响, 进而对时域法识别结果进行修正, 提高识别精度. 最后, 通过一个10自由度轨道车辆垂向动力学模型, 对方法的正确性进行了验证. 研究结果表明: 2种方法对轨道车辆垂向轮轨力均具有较高的识别精度, 对于各轮对的识别精度各有优劣; 在计算效率方面, 状态空间法比Green函数法更优; 经NARX模型修正的2种时域法对轨道车辆垂向轮轨力均具有很好的识别效果, 识别值与正演值的Pearson相关系数大于0.99, 为极强相关. 基于NARX模型的机器学习误差修正方法可有效提高时域识别精度, 可以为后续轨道车辆轮轨力预测提供参考, 具有较强的工程运用价值.Abstract: In order to reduce the identification error of rail vehicle vertical wheel-rail force in time domain, time-domain methods based on machine learning correction are carried out for vertical wheel-rail force identification. Firstly, the vehicle dynamics simulation model is established, and the axle box accelerations and vertical wheel-rail forces are obtained under the speed of 250 km/h with random track irregularity as the excitation. Secondly, the dynamic load identification models corresponding to Green function method and state space method are established. The initial value error of state space method is analyzed, and the polynomial fitting method is introduced to correct the trend term error. The performances of the two methods in calculation accuracy and efficiency are compared. Then, aiming at the identification error existing in the time domain methods, it is proposed to use NARX (nonlinear autoregressive models with exogenous inputs) model to train and predict the identification error, so as to reduce the influence of factors such as incomplete response observation and observation noise in the model, and then correct the identification results of the time-domain methods to improve the identification accuracy. Finally, the correctness of the methods are verified by a 10-degree-of-freedom rail vehicle vertical dynamics model. The results show that the two methods have high identification accuracy for vertical wheel-rail forces of rail vehicles, and for each wheelset, each method has its own advantages and disadvantages; in terms of computational efficiency, state space method is better than Green function method; the two time-domain methods corrected by NARX model are effective in identifying the vertical wheel-rail forces of rail vehicles, and the Pearson correlation coefficients between the simulation and identification values are greater than 0.99, belonging to extremely strong correlation. The machine learning error correction method based on NARX model can effectively improve the time domain identification accuracy, and can provide reference for the subsequent prediction of wheel-rail force of rail vehicles, with strong engineering application value.
-
引言
轨道车辆的轮轨接触力是评判其运行安全性与稳定性的重要依据, 然而轮轨接触力往往难以直接获取. 目前轮轨力的获取主要有直接测量法[1-3]和间接测量法[4-6]. 直接测量法主要是采用测力轮对或测力钢轨, 均需要复杂的标定过程, 且需要后续维护, 较难以推广使用[7]. 间接测量法即使用动载荷识别方法, 利用测量的车辆响应, 识别出精度可以接受的轮轨力值. 在诸多载荷识别方法中, 时域法凭借计算过程直观, 可识别的载荷类型多, 可识别短样本数据(克服了频域法的缺点), 并且还可做到实时识别等优点, 获得了广泛的关注与应用[8-11]. 此外, 机器学习方法由于具有非线性映射能力强等优点, 被广泛应用于复杂机械系统中.
Ronasi 等[12]提出了一种轮轨力识别的优化方法, 将车轮视为二维圆盘, 建立其有限元模型, 通过车轮的径向应变时间历程来识别轮轨力, 并通过加入正则化方法, 提高了模型的抗噪性. Wu等[13]采用Green函数法建立反演模型, 使用TSVD正则化方法结合预测响应误差最小原理, 对轨道车辆轮轨力进行了识别. 朱涛[14]基于动态规划方法, 结合Bellman最优化原理建立了高速列车轮轨力识别模型, 并引入正则化方法消减反问题的不适定性, 获得了较好的识别效果. 郭剑峰[15]提出了一种基于数据建模的轨道车辆轮轨力识别方法, 将经过预处理的实测车辆响应和轮轨力进行特征数据融合, 并作为模型的输入, 采用基于多节点神经网络的方法对轮轨力进行识别, 并与其他类型的神经网络方法进行了对比. Sun等[16]基于数值修正的Newmark-β法, 建立了一个逆向货车动力学模型来识别轮轨接触力, 通过VAMPIRE仿真模型对两个轨道轨顶面和横向不平顺的情况进行了模拟, 然后将模拟的货车部件加速度输入到逆向货车模型中. 对一系悬挂力、二系悬挂力和轮轨接触力的预测与模拟结果进行了比较, 结果趋于一致, 此后, 进一步将逆向货车模型简化为一个3自由度的1/4车垂直模型, 以便于在线应用. 孙善超等[17-18]结合最优控制理论、卡尔曼滤波方法以及轮轨蠕滑理论, 建立了轮轨作用力的全息辨识模型, 将轮轨作用力辨识问题转换为最优控制策略的设计问题, 并结合卡尔曼滤波技术及加速度测试值, 对动载荷预测值进行正向修正, 最终实现了轮轨左、右横向力的识别. 周亚波[19]基于线性卡尔曼滤波算法, 对轨道车辆垂向轮轨力和横向轮轴力进行了辨识, 并验证了在轨道随机不平顺、车轮扁疤及车轮多边形激励下该方法的准确性. Zhu等[20]基于Duhamel积分推导了一种时域动载荷识别方法, 将其应用于车辆通过直线和曲线工况的轮轨力识别, 并将识别结果与动力学仿真结果进行对比, 结果表明该方法具有较高的识别精度. 罗金屯等[21]基于卷积神经网络, 提出了一种数据驱动的轨道车辆轮轨力识别方法, 对轨道车辆直线和曲线工况的轮轨力进行了识别, 并与几种传统的数据驱动方法进行了对比. 曾俊玮等[22]针对胶轮车辆轮胎载荷的识别, 提出了一种卡尔曼滤波器与神经网络相结合的模型, 通过神经网络修正卡尔曼滤波器的识别结果, 提高了轮胎载荷的识别精度.
从已有的研究中可以看出, 以往的研究大多使用单一方法进行轮轨力的识别, 缺乏方法间对比和综合评判; 此外, 运用机器学习修正物理模型, 构建二者的联合模型, 是目前各学科领域研究的新热点. 具体到本文研究领域, 学者们大多单独运用时域法或机器学习方法对轮轨力进行识别, 只有少数研究[22]将二者进行结合. 针对现有研究存在的不足, 本文基于轨道车辆的典型特征, 选取了2种具有代表性的时域动载荷识别方法, 在10自由度轨道车辆垂向动力学模型基础上, 对其垂向轮轨力进行识别; 详细对比2种方法在计算速度、计算精度等方面的优劣; 在此基础上, 针对物理时域方法所存在的识别误差, 引入机器学习算法, 将NARX (nonlinear autoregressive models with exogenous inputs)模型作为误差预测模型, 对时域动载荷识别方法的识别误差进行修正, 进而尝试提高动载荷的识别精度. 本文方法旨在为后续的轨道车辆轮轨力识别, 乃至列车运行安全性评估提供理论支撑.
1. 时域动载荷识别方法
1.1 Green函数法
加速度脉冲响应函数如下式所示[23]
$$ {\boldsymbol{\ddot h}}({{t}}) = \delta ({{t}}){\boldsymbol{\alpha }} - \sum\limits_{n = 1}^\infty {{\omega _n^2}} {{\boldsymbol{A}}_n}\sin ({{\omega}_{{{d,n}}}}{{t + }}{\phi _{{n}}}) $$ (1) 式中, ${{\boldsymbol{A}}_n} = \dfrac{{{{\boldsymbol{\alpha }}_n}}}{{{\omega _{d,n}}}}{{\rm{e}}^{ - {\zeta _n}{\omega _n}t}}$; ${\phi _{{n}}} = {\tan ^{ - 1}}\left( {\dfrac{{2{\zeta _n}\sqrt {1 - {\zeta _n^2}} }}{{1 - 2{\zeta _n^2}}}} \right)$; ${{\boldsymbol{\alpha }}_n} = \dfrac{{{{\boldsymbol{\varphi }}_n}^{}{{\boldsymbol{\varphi }}_n}^{{\rm{T}}}}}{{{m_n}}}$; ${\boldsymbol{\alpha }} = \displaystyle\sum\limits_{n = 1}^\infty {{{\boldsymbol{\alpha }}_n}}$; $\delta ({{t}})$为狄拉克函数; $ {\omega _{d,n}} $为阻尼模态频率; $ {\omega _n} $为模态频率; $ {\zeta _n} $为模态阻尼比; $ {{\boldsymbol{\varphi }}_n} $为模态向量; $ {m_n} $为模态质量.
通过脉冲响应函数公式, 可以得到每一时刻的脉冲响应函数值, 并由此计算得出振动系统的Green函数矩阵. 对于单自由度系统, 其响应正演表达式如下所示
$$ {\boldsymbol{z}} = {\boldsymbol{HF}} $$ (2) 式中, z为系统的响应向量, z = {z1, z2, ···, zm}T; F为系统所受的载荷向量, F = {F1, F2, ···, Fm}T; zi和Fi分别代表第i个时刻的响应值和载荷值; m为采样点的个数; H为系统的Green函数矩阵. 式(2)可写为如下的形式
$$ \left\{ {\begin{array}{*{20}{c}} {{z_1}} \\ {{z_2}} \\ \vdots \\ {{z_m}} \end{array}} \right\} = \Delta t\left[ {\begin{array}{*{20}{c}} {{h_1}}&{}&{}&{} \\ {{h_2}}&{{h_1}}&{}&{} \\ \vdots & \vdots & \ddots &{} \\ {{h_m}}&{{h_{m - 1}}}& \cdots &{{h_1}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{F_1}} \\ {{F_2}} \\ \vdots \\ {{F_m}} \end{array}} \right\} $$ (3) 式中, hi (i = 1, 2, ···, m)为脉冲响应函数值; $ \Delta t $为时间步长.
对于单自由度系统, 其Green函数矩阵的构建方法是, 通过脉冲响应函数公式计算出各个时刻的脉冲响应函数值hi (i = 1, 2, ···, m), 令矩阵主对角线以上的元素的值为0, 主对角线上的元素为$ \Delta t $h1, 从矩阵第1列到第m列的元素, 其值均从主对角线上的$ \Delta t $h1开始, 第1列为$ \Delta t $h1 ~ $ \Delta t $hm, 第2列为$ \Delta t $h1 ~ $ \Delta t $hm-1, 以此类推, 直到第m列, 其非零元素为主对角线元素$ \Delta t $h1.
与单自由度系统不同, 多自由度系统通过脉冲响应函数公式计算出的各个时刻的脉冲响应函数值hi (i = 1, 2, ···, m)均为矩阵, 因此多自由度系统的Green函数矩阵可以看做众多形同单自由度Green函数矩阵的子矩阵的组合, 各子矩阵中的元素分别为hi中对应位置的值与$ \Delta t $之积. 多自由度系统的正演模型如下所示
$$ \left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{z}}_1}} \\ {{{\boldsymbol{z}}_2}} \\ \vdots \\ {{{\boldsymbol{z}}_N}} \end{array}} \right\} = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{H}}_1^1}&{{\boldsymbol{H}}_1^2}& \cdots &{{\boldsymbol{H}}_1^M} \\ {{\boldsymbol{H}}_2^1}&{{\boldsymbol{H}}_2^2}& \cdots &{{\boldsymbol{H}}_2^M} \\ \vdots & \vdots & \vdots & \vdots \\ {{\boldsymbol{H}}_N^1}&{{\boldsymbol{H}}_N^2}& \cdots &{{\boldsymbol{H}}_N^M} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{F}}_1}} \\ {{{\boldsymbol{F}}_2}} \\ \vdots \\ {{{\boldsymbol{F}}_M}} \end{array}} \right\} $$ (4) 式中, zi为第i个位置的响应列阵; Fi为第i个位置的载荷列阵; ${\boldsymbol{H}}_{{j}}^{{i}}$为载荷作用点i和响应测量点j之间的Green函数矩阵; M和N分别为作用在系统上的动载荷个数和响应测量点的个数. 为保证该线性方程组有解, 需要满足$ N \geqslant M $.
在正演模型建立完成后, 通过对Green函数矩阵求逆(或求广义逆)来对动载荷进行识别, 单自由度系统和多自由度系统的动载荷识别公式分别如下所示
$$ {\boldsymbol{F }}={{\boldsymbol{H}}^{ - 1}}{\boldsymbol{z}} $$ (5) $$ \left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{F}}_1}} \\ {{{\boldsymbol{F}}_2}} \\ \vdots \\ {{{\boldsymbol{F}}_M}} \end{array}} \right\} = {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{H}}_1^1}&{{\boldsymbol{H}}_1^2}& \cdots &{{\boldsymbol{H}}_1^M} \\ {{\boldsymbol{H}}_2^1}&{{\boldsymbol{H}}_2^2}& \cdots &{{\boldsymbol{H}}_2^M} \\ \vdots & \vdots & \vdots & \vdots \\ {{\boldsymbol{H}}_N^1}&{{\boldsymbol{H}}_N^2}& \cdots &{{\boldsymbol{H}}_N^M} \end{array}} \right]^{{\boldsymbol{ - }}1}}\left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{z}}_1}} \\ {{{\boldsymbol{z}}_2}} \\ \vdots \\ {{{\boldsymbol{z}}_N}} \end{array}} \right\} $$ (6) 1.2 状态空间法
对于一个n自由度线性定常系统, 其动力学方程如下所示
$$ {\boldsymbol{M\ddot x}} + {\boldsymbol{C\dot x}} + {\boldsymbol{Kx}} = {\boldsymbol{f}} $$ (7) 式中, M, C和K分别为系统质量、阻尼和刚度矩阵; $ {\boldsymbol{\ddot x}} $, $ {\boldsymbol{\dot x}} $和$ {\boldsymbol{x}} $分别为系统加速度、速度和位移向量; $ {\boldsymbol{f}} $为外载荷向量.
将上式转化为状态空间方程的形式, 为
$$ {\boldsymbol{\dot X}} = {\boldsymbol{AX}} + {\boldsymbol{Bf}} $$ (8) 式中, $ {\boldsymbol{X}} $为状态向量, 由位移和速度向量组成; $ {\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} {\boldsymbol{0}}&{\boldsymbol{ I}} \\ { - {{\boldsymbol{M}}^{ - 1}}{\boldsymbol{K}}}&{ - {{\boldsymbol{M}}^{ - 1}}{\boldsymbol{C}}} \end{array}} \right] $; I为单位矩阵; $ {\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} {\boldsymbol{0}} \\ {{{\boldsymbol{M}}^{ - 1}}} \end{array}} \right] $.
由文献[24]可知, 一般的离散状态空间方程只能通过位移和速度对动载荷进行识别, 当已知系统加速度时, 需引入Newmark积分对状态空间方程进行离散, 有如下表达式
$$ {{\boldsymbol{X}}_{{{i}} + 1}} = {\boldsymbol{\varPhi }}{{\boldsymbol{X}}_i} + {\boldsymbol{\varGamma }}{{\boldsymbol{f}}_i} $$ (9) 式中, $ {{\boldsymbol{X}}_i} $为状态变量, 此时的状态变量由位移、速度、加速度向量组成, $ {{\boldsymbol{X}}_i}{{ = }}{\left[ {{{\boldsymbol{x}}_i},{{{\boldsymbol{\dot x}}}_i},{{{\boldsymbol{\ddot x}}}_i}} \right]^{{\rm{T}}} } $.
$$\begin{split} &{\boldsymbol{\varGamma }}{{ = }}{\left[ {\begin{array}{*{20}{c}} {\boldsymbol{I}}&{\boldsymbol{0}}&{ - \beta {\boldsymbol{I}}\Delta {{{t}}^2}} \\ {\boldsymbol{0}}&{\boldsymbol{I}}&{ - \gamma {\boldsymbol{I}}\Delta {{t}}} \\ {\boldsymbol{K}}&{\boldsymbol{C}}&{\boldsymbol{M}} \end{array}} \right]^{ - 1}}\left[ {\begin{array}{*{20}{c}} {\boldsymbol{0}} \\ {\boldsymbol{0}} \\ {\boldsymbol{I}} \end{array}} \right] \\ &{\boldsymbol{\varPhi }}{{ = }}{\left[ {\begin{array}{*{20}{c}} {\boldsymbol{I}}&{\boldsymbol{0}}&{ - \beta {\boldsymbol{I}}\Delta {{{t}}^2}} \\ {\boldsymbol{0}}&{\boldsymbol{I}}&{ - \gamma {\boldsymbol{I}}\Delta {{t}}} \\ {\boldsymbol{K}}&{\boldsymbol{C}}&{\boldsymbol{M}} \end{array}} \right]^{ - 1}}\left[ {\begin{array}{*{20}{c}} {\boldsymbol{I}}&{{\boldsymbol{I}}\Delta {{t}}}&{(0.5 - \beta ){\boldsymbol{I}}\Delta {{{t}}^2}} \\ {\boldsymbol{0}}&{\boldsymbol{I}}&{(1 - \gamma ){\boldsymbol{I}}\Delta {{t}}} \\ {\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{0}} \end{array}} \right] \end{split}$$ 式中, $ \beta $和$ \gamma $为Newmark积分参数, $ \beta $一般在0.25到0.5之间选取, $ \gamma $在0到1之间选取[24].
由式(9)可得每一步的外载荷的求解表达式, 如下所示
$$\begin{split} &\left[ {\begin{array}{*{20}{c}} {\boldsymbol{0}} \\ {\boldsymbol{0}} \\ {{{\boldsymbol{f}}_{i{{ + }}1}}} \end{array}} \right]{{ = }}\left[ {\begin{array}{*{20}{c}} {\boldsymbol{I}}&{\boldsymbol{0}}&{ - \beta {\boldsymbol{I}}\Delta {{{t}}^2}} \\ {\boldsymbol{0}}&{\boldsymbol{I}}&{ - \gamma {\boldsymbol{I}}\Delta {{t}}} \\ {\boldsymbol{K}}&{\boldsymbol{C}}&{\boldsymbol{M}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{x}}_{i + 1}}} \\ {{{{\boldsymbol{\dot x}}}_{i + 1}}} \\ {{{{\boldsymbol{\ddot x}}}_{i + 1}}} \end{array}} \right] -\\ &\qquad \left[ {\begin{array}{*{20}{c}} {\boldsymbol{I}}&{{\boldsymbol{I}}\Delta {{t}}}&{(0.5 - \beta ){\boldsymbol{I}}\Delta {{{t}}^2}} \\ {\boldsymbol{0}}&{\boldsymbol{I}}&{(1 - \gamma ){\boldsymbol{I}}\Delta {{t}}} \\ {\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{0}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{x}}_i}} \\ {{{{\boldsymbol{\dot x}}}_i}} \\ {{{{\boldsymbol{\ddot x}}}_i}} \end{array}} \right]\end{split} $$ (10) 将式(10)展开, 可得到下式
$$ {{\boldsymbol{x}}_{{{i + }}1}}{{ = }}{{\boldsymbol{x}}_i} + \Delta {{t}}{{\boldsymbol{\dot x}}_i} + \Delta {t^2}\left[ {(0.5 - \beta ){{{\boldsymbol{\ddot x}}}_i}{{ + }}\beta {{{\boldsymbol{\ddot x}}}_{i + 1}}} \right] $$ (11) $$ {{\boldsymbol{\dot x}}_{{{i + }}1}}{{ = }}{{\boldsymbol{\dot x}}_i} + \Delta t\left[ {(1 - \gamma ){{{\boldsymbol{\ddot x}}}_i}{{ + }}\gamma {{{\boldsymbol{\ddot x}}}_{i + 1}}} \right] $$ (12) $$ {{\boldsymbol{f}}_{{{i + }}1}}{{ = }}{\boldsymbol{K}}{{\boldsymbol{x}}_{i + 1}} + {\boldsymbol{C}}{{\boldsymbol{\dot x}}_{i + 1}} + {\boldsymbol{M}}{{\boldsymbol{\ddot x}}_{i + 1}} $$ (13) 因此在已知加速度响应的情况下, 通过式(10)求解外载荷在本质上就是通过式(11) ~ 式(13)求解外载荷.
在刚度和阻尼矩阵中, 并不是所有位置的响应都与外载荷的计算相关, 为了提高方法的稳定性和精度, 找到对外载荷${{\boldsymbol{f}}_{{{i + }}1}}$的计算有影响的测点位置(即影响系数法的概念), 把这些测点位置的加速度响应作为已知量, 假设已知i时刻(初始时刻)的位移、速度和加速度, 根据式(11)和式(12)可以得到i + 1时刻的位移和速度, 再根据式(13)即可求得i + 1时刻的动载荷. 对于i时刻位移和速度未知的情况, 由初值引起的误差将在后文中进行讨论. 该方法较适用于多自由度弹簧–阻尼–质量模型, 如轨道车辆动力学模型.
2. 初值误差分析与修正
文献[25]中给出了Wilson-θ反分析法在非0初始条件下由初值引起的误差, 并引入多项式拟合趋势项对位移和速度进行了修正. 由于本文采用的Newmark积分离散的状态空间法同样需要对加速度进行数值积分来间接计算速度和位移, 可能会存在相同的问题, 因此同样需要分析非零初始条件下的初值误差对识别结果的影响.
为了提高方法的稳定性和精度, 同样以对外载荷计算有影响的测点位置的加速度为输入, 若系统为非0初始状态, 由下式所示
$$ {{\boldsymbol{x}}_1}{{ = }}{{\boldsymbol{x}}_0} + \Delta {{t}}{{\boldsymbol{\dot x}}_0} + \Delta {t^2}\left[ {(0.5 - \beta ){{{\boldsymbol{\ddot x}}}_0}{{ + }}\beta {{{\boldsymbol{\ddot x}}}_1}} \right] $$ (14) $$ {{\boldsymbol{\dot x}}_1}{{ = }}{{\boldsymbol{\dot x}}_0} + \Delta t\left[ {(1 - \gamma ){{{\boldsymbol{\ddot x}}}_0}{{ + }}\gamma {{{\boldsymbol{\ddot x}}}_1}} \right] $$ (15) 在1时刻这些测点的位移实际值$ {{\boldsymbol{x}}_1} $可由0时刻的测点的位移$ {{\boldsymbol{x}}_0} $、速度$ {{\boldsymbol{\dot x}}_0} $、加速度$ {{\boldsymbol{\ddot x}}_0} $以及1时刻的加速度$ {{\boldsymbol{\ddot x}}_1} $表示, 1时刻的速度实际值$ {{\boldsymbol{\dot x}}_1} $可由0时刻的速度$ {{\boldsymbol{\dot x}}_0} $、加速度$ {{\boldsymbol{\ddot x}}_0} $以及1时刻的加速度$ {{\boldsymbol{\ddot x}}_1} $表示.
由于$ {{\boldsymbol{x}}_0} $和$ {{\boldsymbol{\dot x}}_0} $未知, 若此时仍将系统视为零初始状态, 则带误差的1时刻的这些测点的位移${{{\boldsymbol{x}'_1}}}$和速度${{{\dot {\boldsymbol{x}}'_1}}}$计算值分别为
$$ {{{{\boldsymbol{x}}'_1}}}{{ = }}\Delta {t^2}\left[ {(0.5 - \beta ){{{\boldsymbol{\ddot x}}}_0}{{ + }}\beta {{{{\ddot {\boldsymbol{x}}_1}}}}} \right] $$ (16) $$ {{{\dot {\boldsymbol{x}}'_1}}}{{ = }}\Delta t\left[ {(1 - \gamma ){{{\boldsymbol{\ddot x}}}_0}{{ + }}\gamma {{{\boldsymbol{\ddot x}}}_1}} \right] $$ (17) 令$ {{\boldsymbol{\delta }}_{11}} $为第1步的位移计算误差, $ {{\boldsymbol{\delta }}_{21}} $为第1步的速度计算误差, 则有
$$ \begin{split} & {{\boldsymbol{\delta }}_{11}}{{ = }}{{\boldsymbol{x}}_1} - {{{{{\boldsymbol{x}}'_1}}}} {{ = }}{{\boldsymbol{x}}_0} + \Delta {{t}}{{{\boldsymbol{\dot x}}}_0} + \Delta {t^2}\left[ {(0.5 - \beta ){{{\boldsymbol{\ddot x}}}_0}{{ + }}\beta {{{\boldsymbol{\ddot x}}}_1}} \right] - \\ &\qquad \Delta {t^2}\left[ {(0.5 - \beta ){{{\boldsymbol{\ddot x}}}_0}{{ + }}\beta {{{\boldsymbol{\ddot x}}}_1}} \right] {{ = }}{{\boldsymbol{x}}_0} + \Delta {{t}}{{{\boldsymbol{\dot x}}}_0} \end{split} $$ (18) $$ \begin{split} &{{\boldsymbol{\delta }}_{21}}{{ = }}{{{\boldsymbol{\dot x}}}_1} - {{{{\dot {\boldsymbol{x}}'_1}}}} {{ = }}{{{\boldsymbol{\dot x}}}_0} + \Delta t\left[ {(1 - \gamma ){{{\boldsymbol{\ddot x}}}_0}{{ + }}\gamma {{{\boldsymbol{\ddot x}}}_1}} \right]- \\ & \qquad \Delta t\left[ {(1 - \gamma ){{{\boldsymbol{\ddot x}}}_0}{{ + }}\gamma {{{\boldsymbol{\ddot x}}}_1}} \right] {{ = }}{{{\boldsymbol{\dot x}}}_0} \end{split}$$ (19) 继续计算第2时刻测点的位移、速度的实际值与计算值
$$ {{\boldsymbol{x}}_2}{{ = }}{{\boldsymbol{x}}_1} + \Delta {{t}}{{\boldsymbol{\dot x}}_1} + \Delta {t^2}\left[ {(0.5 - \beta ){{{\boldsymbol{\ddot x}}}_1}{{ + }}\beta {{{\boldsymbol{\ddot x}}}_2}} \right] $$ (20) $$ {{\boldsymbol{\dot x}}_2}{{ = }}{{\boldsymbol{\dot x}}_1} + \Delta t\left[ {(1 - \gamma ){{{\boldsymbol{\ddot x}}}_1}{{ + }}\gamma {{{\boldsymbol{\ddot x}}}_2}} \right] $$ (21) $$ {{{{\boldsymbol{x}}'_2}}}{{ = }}{{\boldsymbol{x}'_1}} + \Delta {{t}}{{{\dot {\boldsymbol{x}}'_1}}} + \Delta {t^2}\left[ {(0.5 - \beta ){{{\boldsymbol{\ddot x}}}_1}{{ + }}\beta {{{\boldsymbol{\ddot x}}}_2}} \right] $$ (22) $$ {{\boldsymbol{\dot x'}}_2}{{ = }}{{\boldsymbol{\dot x'}}_1} + \Delta t\left[ {(1 - \gamma ){{{\boldsymbol{\ddot x}}}_1}{{ + }}\gamma {{{\boldsymbol{\ddot x}}}_2}} \right] $$ (23) 令$ {{\boldsymbol{\delta }}_{12}} $为第2步的位移计算误差, $ {{\boldsymbol{\delta }}_{22}} $为第2步的速度计算误差, 则有
$$ \begin{split} & {{\boldsymbol{\delta }}_{12}}{{ = }}{{\boldsymbol{x}}_2} - {{{{{\boldsymbol{x}}'_2}}}} {{ = }}{{\boldsymbol{x}}_1} + \Delta {{t}}{{{\boldsymbol{\dot x}}}_1} + \Delta {t^2}\left[ {(0.5 - \beta ){{{\boldsymbol{\ddot x}}}_1}{{ + }}\beta {{{\boldsymbol{\ddot x}}}_2}} \right]- \\ &\qquad {{{{{\boldsymbol{x}}'_1}}}} - \Delta {{t}}{{{{\dot {\boldsymbol{x}}'_1}}}} - \Delta {t^2}\left[ {(0.5 - \beta ){{{\boldsymbol{\ddot x}}}_1}{{ + }}\beta {{{\boldsymbol{\ddot x}}}_2}} \right]= \\ &\qquad {{\boldsymbol{x}}_0} + \Delta {{t}}{{{\boldsymbol{\dot x}}}_0} + \Delta {{t}}{{{\boldsymbol{\dot x}}}_0} {{ = }}{{\boldsymbol{x}}_0} + 2\Delta {{t}}{{{\boldsymbol{\dot x}}}_0} \end{split} $$ (24) $$\begin{split} & {{\boldsymbol{\delta }}_{22}}{{ = }}{{{{\dot {\boldsymbol{x}}_2}}}} - {{{{\dot {\boldsymbol{x}}'_2}}}} {{ = }}{{{{\dot {\boldsymbol{x}}_1}}}} + \Delta t\left[ {(1 - \gamma ){{{\boldsymbol{\ddot x}}}_1}{{ + }}\gamma {{{\boldsymbol{\ddot x}}}_2}} \right] - \\ &\qquad {{{{\dot {\boldsymbol{x}}_1'}}}} - \Delta t\left[ {(1 - \gamma ){{{\boldsymbol{\ddot x}}}_1}{{ + }}\gamma {{{\boldsymbol{\ddot x}}}_2}} \right] {{ = }}{{{\boldsymbol{\dot x}}}_1} - {{{{\dot {\boldsymbol{x}}'_1}}}} {{ = }}{{{\boldsymbol{\dot x}}}_0} \end{split}$$ (25) 至此已经可以发现规律, 即第k步的测点位移计算误差$ {{\boldsymbol{\delta }}_{1{{k}}}} $、速度计算误差$ {{\boldsymbol{\delta }}_{2{{k}}}} $分别为
$$\qquad\qquad\quad {{\boldsymbol{\delta }}_{1{{k}}}}{{ = }}{{\boldsymbol{x}}_0} + {{k}}\Delta {{t}}{{\boldsymbol{\dot x}}_0} $$ (26) $$\qquad\qquad\quad {{\boldsymbol{\delta }}_{2{{k}}}}{{ = }}{{\boldsymbol{\dot x}}_0} $$ (27) 即随着迭代步数的增加, 计算得到的测点位移会出现斜率为$ {{\boldsymbol{\dot x}}_0} $, 常数项为$ {{\boldsymbol{x}}_0} $的一次趋势项, 速度会出现大小为$ {{\boldsymbol{\dot x}}_0} $的偏移项, 而这些测点响应又影响着外载荷的计算, 会给每一步的动载荷识别值引入趋势项误差, 这与文献[25]中推导出的Wilson-θ反分析法的误差具有相同的形式, 因此本文采用的Newmark积分离散的状态空间法可以采用与该文献相同的多项式拟合趋势项的方法修正初值误差, 主要的修正步骤如下[25].
假设积分得到的离散数据点为(ti, vi), i = 0, 1, 2,···, n − 1, 此处vi代指积分所得的含趋势项的速度或位移数据点. 设${{{f}}_m}({{t}}) = \displaystyle\sum\limits_{k = 0}^m {{p_k}{t^k}} \in \varPhi$, $ {p_k} $为多项式系数, $ \varPhi $为所有次数不超过$m\;{{ }}({{m}} \leqslant n - 1)$的多项式构成的函数类, 使得
$$ I = \sum\limits_{{{i}} = 0}^{n - 1} {{{[{{{v}}_{{i}}} - {{{f}}_{{m}}}({t_i})]}^2}} = \sum\limits_{{{i}} = 0}^{n - 1} {{{\left( {\sum\limits_{{{k}} = 0}^m {{{{v}}_{{i}}} - {{{p}}_k}{t_i}^k} } \right)}^2}} = \min $$ (28) 要确定趋势项, 就需要求得一组$ {{{p}}_k} $使得$ I $值达到最小, 由多元函数求极值的必要条件可得
$$ \frac{{\partial I}}{{\partial {{{p}}_j}}} = 2\sum\limits_{i = 0}^{n - 1} {t_i^j\left( {\sum\limits_{k = 0}^m {{v_i} - {p_k}t_i^k} } \right)} = 0 $$ (29) 其中j$ = $0, 1, 2, ···, m, 即
$$ \sum\limits_{k = 0}^m {\left( {\sum\limits_{{{i}} = 0}^{n - 1} {t_i^{j + k}} } \right)} {p_k} = \sum\limits_{i = 0}^{n - 1} {t_i^j{v_i}} $$ (30) 将其表示为矩阵形式, 如下式所示
$$\begin{split} &\left[ {\begin{array}{*{20}{c}} {{n}}&{\displaystyle\sum\limits_{i = 0}^{n - 1} {{t_i}} }& \cdots &{\displaystyle\sum\limits_{i = 0}^{n - 1} {{t_i}^m} } \\ {\displaystyle\sum\limits_{i = 0}^{n - 1} {{t_i}} }&{\displaystyle\sum\limits_{i = 0}^{n - 1} {{t_i}^2} }& \cdots &{\displaystyle\sum\limits_{i = 0}^{n - 1} {{t_i}^{m + 1}} } \\ \vdots & \vdots & \vdots & \vdots \\ {\displaystyle\sum\limits_{i = 0}^{n - 1} {{t_i}^m} }&{\displaystyle\sum\limits_{i = 0}^{n - 1} {{t_i}^{m + 1}} }& \cdots &{\displaystyle\sum\limits_{i = 0}^{n - 1} {{t_i}^{2m}} } \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{p_0}} \\ {{p_1}} \\ \vdots \\ {{p_m}} \end{array}} \right] =\\ &\qquad {\left[ {\sum\limits_{i = 0}^{n - 1} {{v_i},{{ }}\sum\limits_{i = 0}^{n - 1} {{t_{{i}}}{{{v}}_{{i}}},{{ }} \cdots ,{{ }}\sum\limits_{i = 0}^{n - 1} {{t_{{i}}}^{{m}}{{{v}}_{{i}}}} } } } \right]^{{\rm{T}}} } \end{split}$$ (31) 可以证明, 方程组系数矩阵是一个对称正定矩阵, 故存在唯一解. 通过上式求解出pk, 从而可得拟合多项式${{{f}}_m}({{t}}) = \displaystyle\sum\limits_{k = 0}^m {{p_k}{t^k}}$.
除了初值误差引入的趋势项之外, 相关研究[26-27]也表明, 对于由信号本身或其他因素所引起的积分信号的曲线趋势项, 多项式拟合修正的方法同样也是有效的.
3. 基于NARX模型的动载荷识别误差修正
采用基于纯物理模型的时域法对动载荷进行识别时, 难免会因为各种因素的影响(如响应观测不全, 观测噪声等)造成识别误差. 为了减小识别误差, 参考相关文献[22, 28]中提出的物理–机器学习相结合的模型思路, 本节引入NARX模型, 利用其强大的非线性逼近能力, 通过一段已知的输入–输出时间历程对模型进行训练, 对基于物理模型时域法的动载荷识别误差进行预测, 进而尝试对识别误差进行修正, 从而提高动载荷识别精度. 这样既避免了纯机器学习模型缺乏物理解释的问题, 又减小了纯物理模型的识别误差.
NARX模型是线性ARX (autoregressive models with exogenous inputs)模型的扩展, 与静态回归模型相比, 增加了多步时延的环节, 并引入了非线性映射函数, 因此NARX模型具有动态特性, 可以利用输入、输出的历史信息来对具有复杂映射关系的未知时间历程进行预测. 对于一个单输入单输出系统, 其NARX模型可以表示为
$$\begin{split} &{{f}}({{t}}){{ = F}}({{f}}({{t}} - 1),{{f}}({{t}} - 2), \cdots ,{{f}}({{t}} - {{{n}}_a}),\\ &\qquad {{z}}({{t}} - {n_k}), \cdots ,{{z}}({{t}} - {n_k} - {n_b} + 1)) \end{split}$$ (32) 式中, $ {{f}}({{t}}) $为t时刻的系统输出; $ {{z}}({{t}}) $为t时刻的系统输入; $ {{{n}}_a} $为输出的阶数; $ {{{n}}_{{k}}} $为输入相对于输出延迟的阶数; $ {{{n}}_{{b}}} $为输入的阶数; F(·)为非线性映射函数, 常用的有小波网络、S型网络、二叉树回归模型、支持向量机等. 由于小波网络具有结构简单、收敛速度快及计算精度高的优点, 故本文的F(·)选用小波网络[29].
具体到本文而言, 由于要使用轨道车辆轴箱响应来预测垂向轮轨力的识别误差, 因此本文采用的NARX模型的输入为轨道车辆轴箱响应, 输出为垂向轮轨力的识别误差, 定义识别误差为物理模型时域法的识别值减去轮轨力的真实值(可以含随机噪声), 具体的NARX模型可以表示为框图的形式, 如图1所示.
采用NRMSE适应度来评价NARX模型的预测精度, 如下式所示
$$ NRMS E{{ = }}\left(1 - \frac{{\left\| {{\boldsymbol{f}} - {\boldsymbol{\hat f}}} \right\|}}{{\left\| {{\boldsymbol{f}} - {\boldsymbol{\bar f}}} \right\|}}\right) \times 100\% $$ (33) 式中, $ {\boldsymbol{f}} $为目标输出; $ {\boldsymbol{\bar f}} $为目标输出的均值; $ {\boldsymbol{\hat f}} $为模型的预测输出.
在使用训练集的输入输出以及测试集的输入得到NARX模型的识别误差预测值之后, 以测试集的响应输入得到物理模型时域法的识别结果, 再将该识别结果减去识别误差预测值即为修正后的动载荷识别值, 其框图如图2示.
当非线性映射函数确定之后, 输出的阶数$ {{{n}}_a} $、输入相对于输出延迟的阶数$ {{{n}}_{{k}}} $以及输入的阶数$ {{{n}}_{{b}}} $将影响NARX模型的预测效果, 参数的选取将在下一节进行讨论.
4. 轨道车辆垂向轮轨力识别
根据《车辆−轨道耦合动力学》[30], 车辆垂向动力学模型由车体、前后转向架和4个轮对组成. 对于车体以及前后转向架, 考虑浮沉运动和点头运动; 对于4个轮对, 考虑浮沉运动, 总计10个自由度, 如图3所示. 在本研究中, 根据该模型进行动力学建模, 以轴箱加速度响应为输入, 对轮轨垂向力进行识别. 通过建立车辆垂向动力学微分方程, 可以提取出该系统的质量M、刚度和阻尼矩阵K, C, 外载荷向量如下.
$$ {\boldsymbol{M}} = {{\rm{diag}}} (({{{M}}_{{c}}},{{{I}}_{{{cy}}}},{{{M}}_{{t}}},{{{I}}_{{{ty}}}},{{{M}}_{{t}}},{{{I}}_{{{ty}}}},{{{M}}_{{w}}},{{{M}}_{{w}}},{{{M}}_{{w}}},{{{M}}_{{w}}})) $$ (34) $$ {\boldsymbol{K}} = \left[ {\begin{array}{*{20}{c}} {4{K_{sz}}}&0&{ - 2{K_{sz}}}&0&{ - 2{K_{sz}}}&0&0&0&0&0 \\ 0&{4{K_{sz}}{l_c}^2}&{2{K_{sz}}{l_c}}&0&{ - 2{K_{sz}}{l_c}}&0&0&0&0&0 \\ { - 2{K_{sz}}}&{2{K_{sz}}{l_c}}&{4{K_{pz}} + 2{K_{sz}}}&0&0&0&{ - 2{K_{pz}}}&{ - 2{K_{pz}}}&0&0 \\ 0&0&0&{4{K_{pz}}{l_t}^2}&0&0&{2{K_{pz}}{l_t}}&{ - 2{K_{pz}}{l_t}}&0&0 \\ { - 2{K_{sz}}}&{ - 2{K_{sz}}{l_c}}&0&0&{4{K_{pz}} + 2{K_{sz}}}&0&0&0&{ - 2{K_{pz}}}&{ - 2{K_{pz}}} \\ 0&0&0&0&0&{4{K_{pz}}{l_t}^2}&0&0&{2{K_{pz}}{l_t}}&{ - 2{K_{pz}}{l_t}} \\ 0&0&{ - 2{K_{pz}}}&{2{K_{pz}}{l_t}}&0&0&{2{K_{pz}}}&0&0&0 \\ 0&0&{ - 2{K_{pz}}}&{ - 2{K_{pz}}{l_t}}&0&0&0&{2{K_{pz}}}&0&0 \\ 0&0&0&0&{ - 2{K_{pz}}}&{2{K_{pz}}{l_t}}&0&0&{2{K_{pz}}}&0 \\ 0&0&0&0&{ - 2{K_{pz}}}&{ - 2{K_{pz}}{l_t}}&0&0&0&{2{K_{pz}}} \end{array}} \right] $$ (35) $$ {\boldsymbol{C}} = \left[ {\begin{array}{*{20}{c}} {4{C_{sz}}}&0&{ - 2{C_{sz}}}&0&{ - 2{C_{sz}}}&0&0&0&0&0 \\ 0&{4{C_{sz}}{l_c}^2}&{2{C_{sz}}{l_c}}&0&{ - 2{C_{sz}}{l_c}}&0&0&0&0&0 \\ { - 2{C_{sz}}}&{2{C_{sz}}{l_c}}&{4{C_{pz}} + 2{C_{sz}}}&0&0&0&{ - 2{C_{pz}}}&{ - 2{C_{pz}}}&0&0 \\ 0&0&0&{4{C_{pz}}{l_t}^2}&0&0&{2{C_{pz}}{l_t}}&{ - 2{C_{pz}}{l_t}}&0&0 \\ { - 2{C_{sz}}}&{ - 2{C_{sz}}{l_c}}&0&0&{4{C_{pz}} + 2{C_{sz}}}&0&0&0&{ - 2{C_{pz}}}&{ - 2{C_{pz}}} \\ 0&0&0&0&0&{4{C_{pz}}{l_t}^2}&0&0&{2{C_{pz}}{l_t}}&{ - 2{C_{pz}}{l_t}} \\ 0&0&{ - 2{C_{pz}}}&{2{C_{pz}}{l_t}}&0&0&{2{C_{pz}}}&0&0&0 \\ 0&0&{ - 2{C_{pz}}}&{ - 2{C_{pz}}{l_t}}&0&0&0&{2{C_{pz}}}&0&0 \\ 0&0&0&0&{ - 2{C_{pz}}}&{2{C_{pz}}{l_t}}&0&0&{2{C_{pz}}}&0 \\ 0&0&0&0&{ - 2{C_{pz}}}&{ - 2{C_{pz}}{l_t}}&0&0&0&{2{C_{pz}}} \end{array}} \right] $$ (36) $$ {\boldsymbol{F}} = {\left( {0,0,0,0,0,0,{{{F}}_1},{{{F}}_2},{{{F}}_3},{{{F}}_4}} \right)^{\rm{T}}} $$ (37) 式中, Mc为车体质量; Mt为构架质量; Mw为轮对质量; Icy为车体绕Y轴转动惯量; Ity为构架绕Y轴转动惯量; Ksz为转向架一侧二系悬挂垂向刚度(N/m); Kpz为每轴箱一系悬挂垂向刚度(N/m); Csz为转向架一侧二系悬挂垂向阻尼(N·s/m); Cpz为每轴箱一系悬挂垂向阻尼(N·s/m); F为外载荷向量, F1, F2, F3和F4分别为第1轮对、第2轮对、第3轮对和第4轮对轮轨垂向动载荷; lc为车辆定距之半; lt为车辆固定轴距之半.
采用通用动力学仿真软件建立具有同等悬挂参数的轨道车辆多体动力学模型, 采样频率为1000 Hz, 得到仿真模型在一段随机高低不平顺轨道激励下以250 km/h速度行驶时的4个轮对的垂向轮轨动载荷和轴箱垂向加速度响应, 在此以第1轮对为例, 截取非零初始条件下的4000组采样点作为测试集, 另取10000组采样点作为训练集, 如图4所示. 采用Pearson相关系数来评价方法的识别精度, 如下式所示[14]
![]()
图 4 第1轮对垂向轮轨力及加速度时间历程Figure 4. Time history of vertical wheel-rail force and acceleration of the first wheelset$$ p = \frac{{\displaystyle\sum\limits_{i = 1}^n {({{{{\boldsymbol{P}}}}'_i}} - {{\bar {\boldsymbol{P}}}}')({{\boldsymbol{P}}_i} - {{\bar {\boldsymbol{P}}}})}}{{\sqrt {\displaystyle\sum\limits_{i = 1}^n {{{({{\boldsymbol{P}}_i}' - {{\bar {\boldsymbol{P}}}}')}^2}} } \sqrt {\displaystyle\sum\limits_{i = 1}^n {{{({{\boldsymbol{P}}_i} - {{\bar {\boldsymbol{P}}}})}^2}} } }} $$ (38) 式中, n为采样点的数目; ${{\boldsymbol{P}}_i}'$为识别载荷; ${{\bar {\boldsymbol{P}}}}'$为识别载荷的均值; ${{\boldsymbol{P}}_i}$为外载荷; ${{\bar {\boldsymbol{P}}}}$为外载荷的均值.
考虑到以往的相关研究均直接使用仿真正演的加速度响应作为模型输入, 而实际上, 在进行加速度响应的观测时难以避免地会存在观测噪声, 为了模拟实际, 在本文算例的轴箱加速度响应中添加5%的随机噪声. 分别使用Green函数法和状态空间法对轮轨垂向力进行识别, 并对比识别精度与计算速度, 如表1所示, 2种方法对第1轮对的识别结果对比如图5所示.
表 1 2种时域法对垂向轮轨力的识别结果对比Table 1. Comparison of the two time domain methods for identifying vertical wheel-rail forcesMethod Pearson correlation
coefficient ofthe first
wheelsetPearson correlation
coefficient of the second
wheelsetPearson correlation
coefficient of the third
wheelsetPearson correlation
coefficient of the fourth
wheelsetTime consumption/s Green function method 0.948 0.938 0.946 0.932 19.137 state space method 0.958 0.921 0.956 0.912 0.505 ![]()
图 5 2种方法对第1轮对垂向轮轨力的识别结果对比Figure 5. Comparison of the two methods for identifying the vertical wheel-rail force of the first wheelset从表1中可以看出, 在计算时间方面, 状态空间法具有较快的计算速度, 对4 s时间历程的识别计算时间在0.5 s左右, 而Green函数法因为需要建立规模较大的Green函数矩阵并对其进行求逆, 因而计算效率较低, 计算时间达到了19 s左右; 在计算精度方面, 2种方法识别值与仿真值的Pearson相关系数均大于0.91, 对应极强相关性, 说明2种方法均能较好地识别出轮轨垂向力. 但从图5中可以看出, 识别值与仿真值的波形在一些位置的对应较差, 识别效果还有一定的提升空间, 因此引入NARX模型对识别误差进行修正, 并对2种方法的识别误差进行提取. 考虑到篇幅的限制, 在后文中只针对第一轮对进行修正, 其余轮对的修正步骤与之相同, 不再进行展开.
为了使结果更具一般性, 同时考察模型对噪声的抵抗能力, 在训练集和测试集的轴箱加速度与识别误差中也引入5%的随机噪声以模拟实际观测噪声. 对于Green函数法, 由于载荷识别结果仅与加速度有关, 因此使用轴箱加速度作为误差预测模型的输入即可; 而对于状态空间法, 经试算发现直接使用轴箱加速度来预测载荷误差较为困难, 预测效果差, 分析其原因可知, 这种方法的载荷识别结果与加速度、速度和位移均有关, 所用到的速度、位移是由加速度间接计算得到的, 由于轴箱加速度本来就是非平稳信号, 其中还含有噪声, 因此间接计算得到的位移即使在使用更高次的多项式拟合趋势项修正后仍然存在一定的误差, 而经间接计算得到的速度误差较小, 因此位移所含的误差对载荷识别结果的影响相对更大, 故选取计算并修正后得到的轴箱位移作为该方法误差预测模型的输入. 在进行修正之前, 还需要寻找NARX模型的误差预测效果和输入阶数$ {{{n}}_{{b}}} $、输出阶数$ {{{n}}_{{a}}} $以及输入相对于输出延迟的阶数$ {{{n}}_{{k}}} $的关系. 需要进行说明的是, 一般的机器学习研究对于模型超参数的选取常在验证集进行, 而本文研究的重点主要是提供算法思路并验证其可行性, 而不是超参数的选取, 故没有再设置验证集, 而是在测试集进行参数选取. 事实上, 本文测试集起到的是验证集的作用, 后续如需进一步验证模型的泛化能力, 再另外选取新的测试集进行测试即可.
在本例中, 暂不考虑$ {{{n}}_{{k}}} $的影响, 只通过改变输入阶数$ {{{n}}_{{b}}} $、输出阶数$ {{{n}}_{{a}}} $来达到一个较好的误差预测效果. 同样以第一轮对为例, 设置$ {{{n}}_{{a}}} $, $ {{{n}}_{{b}}} $的值以5为步长递增, 为了防止输入输出阶数过大而导致的训练发散, 将其最大值设置为150, 观察测试集的NRMSE适应度随着$ {{{n}}_{{a}}} $和$ {{{n}}_{{b}}} $值增加的变化情况, 如图6所示.
![]()
图 6 测试集预测误差NRMSE适应度随na和nb值的变化关系Figure 6. The relationships between the NRMSE of predicted error and the values of na and nb in the test set从图6中可以看出, 2种方法预测误差的NRMSE适应度随着na和nb值的增大, 整体呈现出增大的趋势, 当na和nb值达到60左右时, NRMSE值增长趋势变缓, 且均大于50%. 根据图中呈现的关系, 选取2种方法的na和nb均为150, 2种方法在测试集的误差原始值(含5%随机噪声)和预测值如图7所示, 对应的NRMSE适应度如表2所示.
![]()
图 7 测试集的误差原始值(含5%随机噪声)和预测值Figure 7. The original error (including 5% random noise) and the predicted error in the test set从图7中可以看出, NARX模型能较好地预测2种方法的误差趋势, 在训练集识别误差带5%随机噪声的情况下, 仍然能较准确地预测出2种方法的识别误差, 具有一定的抗噪能力, NRMSE适应度均大于73%. 将图7中的误差预测值记为err, 用于修正图5中的垂向轮轨力识别值(记为F*), 则有轮轨力最终识别值F为
$$ {\boldsymbol{F}} = {{\boldsymbol{F}}^*} - {\boldsymbol{err}} $$ (39) 经修正后的轮轨力最终识别值F如图8所示.
表 2 预测误差的NRMSE适应度Table 2. NRMSE of the predicted errorMethod NRMSE/% Green function method 78.641 state space method 73.166 从图8中可以看到, 经NARX误差预测模型修正后的2种方法对第一轮对的最终识别结果均能很好地逼近原始载荷, 识别误差得到了很好的消除, 其Pearson相关系数分别为0.998和0.997, 对应极强相关性, 相对于未修正的结果而言, Pearson相关系数分别提升了0.050和0.039. 可见, 本文所采用的时域法结合NARX误差修正的动载荷识别模型能很好地识别轨道车辆垂向轮轨力, 并具有一定的抗噪能力.
5. 结论
本文采用基于NARX模型的机器学习方法对两种常用的动载荷时域识别方法进行修正, 在对两方法计算速度、计算精度对比的基础上, 完成了轨道车辆垂向轮轨力的识别, 主要得到以下结论.
(1) 提出了经典时域法结合机器学习修正的动载荷识别方法, 基于两种方法对轨道车辆垂向轮轨力进行识别, 通过机器学习方法修正后的载荷识别误差明显减小, 显著提高了轮轨力的识别精度.
(2) 基于Green函数法和状态空间法的轨道车辆垂向轮轨力识别值与正演值的Pearson相关系数大于0.91, 但两方法的识别精度均还有提升空间; 在计算效率方面, 状态空间法对4 s样本的识别计算时间在0.5 s左右, 具有较快的计算速度, 而Green函数法由于需要建立较大规模的Green函数矩阵并对其进行求逆, 因而计算效率较低, 计算时间达到了19 s左右.
(3) 针对两种方法的误差产生原因, 引入基于NARX模型的机器学习方法, 通过对两方法识别误差的训练和预测, 修正原始识别结果. 研究结果表明: 采用NARX模型修正时域法识别误差的这一模式具有一定的可行性, 可以较好地消除由响应观测不全、随机噪声等因素引起的轮轨力识别误差. 经修正后, 2种方法对第1轮对的最终识别值与仿真值的Pearson相关系数均大于0.99, 识别精度得到了较大的提高, 说明本文方法为时域法识别精度的进一步提高提供了一种思路, 具有一定的工程应用价值.
在后续的研究中, 还可以考虑对轨道车辆横向轮轨力以及带有部分悬挂参数非线性的轨道车辆模型的轮轨力进行识别, 进一步研究本文方法在更复杂的模型中的载荷识别能力.
-
![]()
图 4 第1轮对垂向轮轨力及加速度时间历程
Figure 4. Time history of vertical wheel-rail force and acceleration of the first wheelset
![]()
图 5 2种方法对第1轮对垂向轮轨力的识别结果对比
Figure 5. Comparison of the two methods for identifying the vertical wheel-rail force of the first wheelset
![]()
图 6 测试集预测误差NRMSE适应度随na和nb值的变化关系
Figure 6. The relationships between the NRMSE of predicted error and the values of na and nb in the test set
![]()
图 7 测试集的误差原始值(含5%随机噪声)和预测值
Figure 7. The original error (including 5% random noise) and the predicted error in the test set
表 1 2种时域法对垂向轮轨力的识别结果对比
Table 1 Comparison of the two time domain methods for identifying vertical wheel-rail forces
Method Pearson correlation
coefficient ofthe first
wheelsetPearson correlation
coefficient of the second
wheelsetPearson correlation
coefficient of the third
wheelsetPearson correlation
coefficient of the fourth
wheelsetTime consumption/s Green function method 0.948 0.938 0.946 0.932 19.137 state space method 0.958 0.921 0.956 0.912 0.505 
下载: 导出CSV
表 2 预测误差的NRMSE适应度
Table 2 NRMSE of the predicted error
Method NRMSE/% Green function method 78.641 state space method 73.166
下载: 导出CSV
-
[1] 丁奥. 地面轮轨力连续测试方法和装置研究. [硕士论文]. 成都: 西南交通大学, 2019 Ding Ao. Research on the continuously testing method and device of wheel-rail force on the ground. [Master Thesis]. Chengdu: Southwest Jiaotong University, 2019 (in Chinese)
[2] Bizic M, Petrovic D. Algorithm for inverse determination of derailment coefficient by using instrumented wheelsets. International Journal of Heavy Vehicle Systems, 2022, 29(5): 503-517 doi: 10.1504/IJHVS.2022.128915
[3] 赖伊雯, 陈建政, 马添翼等. 基于全相位滤波的高频轮轨力测量方法研究. 机械, 2022, 49(12): 53-60 (Lai Yiwen, Chen Jianzheng, Ma Tianyi, et al. Research on high-frequency wheel-rail force measurement method based on all-phase filter. Machinery, 2022, 49(12): 53-60 (in Chinese) doi: 10.3969/j.issn.1006-0316.2022.12.009 Lai Yiwen, Chen Jianzheng, Ma Tianyi, et al. Research on high-frequency wheel-rail force measurement method based on all-phase filter. Machinery, 2022, 49(12): 53-60 (in Chinese) doi: 10.3969/j.issn.1006-0316.2022.12.009
[4] Liu Q, Lei X, Rose JG, et al. Vertical wheel-rail force waveform identification using wavenumber domain method. Mechanical Systems and Signal Processing, 2021, 159(8): 107784
[5] Teng F, Zhu R, Zhou Y, et al. A lightweight model of wheel-rail force inversion for railway vehicles. Concurrency and Computation Practice and Experience, 2023, 35: e6443
[6] 王明猛, 朱涛, 王小瑞等. 一种逆结构滤波法的轨道车辆轮轨力识别. 振动工程学报, 2019, 32(4): 602-608 (Wang Mingmeng, Zhu Tao, Wang Xiaorui, et al. An inverse structural filter method for wheel-rail contact forces identification of railway vehicles. Journal of Vibration Engineering, 2019, 32(4): 602-608 (in Chinese) doi: 10.16385/j.cnki.issn.1004-4523.2019.04.006 Wang Mingmeng, Zhu Tao, Wang Xiaorui, et al. An inverse structural filter method for wheel-rail contact forces identification of railway vehicles. Journal of Vibration Engineering, 2019, 32(4): 602-608 (in Chinese) doi: 10.16385/j.cnki.issn.1004-4523.2019.04.006
[7] Zeng J, Wei L, Wu P. Safety evaluation for railway vehicles using an improved indirect measurement method of wheel-rail forces. Journal of Modern Transportation, 2016, 24(2): 114-123 doi: 10.1007/s40534-016-0107-5
[8] 蔡元奇. 时域内动态载荷识别理论及实施技术研究. [博士论文]. 武汉: 武汉大学, 2004 Cai Yuanqi. Theory and applying technique study on identification method for dynamic loads in time domain. [PhD Thesis]. Wuhan: Wuhan University, 2004 (in Chinese)
[9] Liu J, Sun X, Han X, et al. Dynamic load identification for stochastic structures based on Gegenbauer polynomial approximation and regularization method. Mechanical Systems and Signal Processing, 2015, 56-57: 35-54
[10] Wang L, Huang Y, Xie Y, et al. A new regularization method for dynamic load identification. Science Progress, 2020, 103(3): 0036850420931283
[11] Liu R, Hou Z, Wu P, et al. Dynamic load identification for a power battery pack based on a combined regularization algorithm. Journal of Sound and Vibration, 2022, 529: 116928 doi: 10.1016/j.jsv.2022.116928
[12] Ronasi H, Johansson H, Larsson F. A numerical framework for load identification and regularization with application to rolling disc problem. Computers and Structures, 2011, 89(1-2): 38-47 doi: 10.1016/j.compstruc.2010.07.009
[13] Wu J, Zhu T, Wang Y, et al. TSVD regularization-parameter selection method based on Wilson-θ and its application to vertical wheel-rail force identification of rail vehicles. Shock and Vibration, 2022, 2022: 2598040
[14] 朱涛. 高速列车载荷反演技术及其运用研究. [博士论文]. 成都: 西南交通大学, 2012 Zhu Tao. Load identification technology for high-speed train dynamic force and its application. [PhD Thesis]. Chengdu: Southwest Jiaotong University, 2012 (in Chinese)
[15] 郭剑峰. 基于数据建模的轮轨力载荷辨识理论和应用研究. [博士论文]. 北京: 中国铁道科学研究院, 2015 Guo Jianfeng. Theory and application research on wheel rail force load identification based on data modeling. [PhD Thesis]. Beijing: China Academy of Railway Sciences, 2015 (in Chinese)
[16] Sun Y, Cole C, Spiryagin M. Monitoring vertical WR contact force based on freight wagon inverse modeling. Advances in Mechanical Engineering, 2015, 7(5): 1-11
[17] 孙善超, 刘金朝, 王卫东等. 高速铁路轮轨作用力的全息辨识模型. 工程力学, 2018, 35(11): 190-196 (Sun Shanchao, Liu Jinzhao, Wang Weidong, et al. Holographic identification model of wheel & rail contact force for high-speed railway. Engineering Mechanics, 2018, 35(11): 190-196 (in Chinese) doi: 10.6052/j.issn.1000-4750.2017.08.0603 Sun Shanchao, Liu Jinzhao, Wang Weidong, et al. Holographic identification model of wheel & rail contact force for high-speed railway. Engineering Mechanics, 2018, 35(11): 190-196 (in Chinese) doi: 10.6052/j.issn.1000-4750.2017.08.0603
[18] 孙善超. 轨道−车辆系统轮轨力辨识及应用研究. [博士论文]. 北京: 中国铁道科学研究院, 2016 Sun Shanchao. Theoretical research on wheel/rail contact force identification and its application. [PhD Thesis]. Beijing: China Academy of Railway Sciences, 2016 (in Chinese)
[19] 周亚波. 基于铁道车辆振动加速度的轮轨力反演方法研究. [硕士论文]. 成都: 西南交通大学, 2020 Zhou Yabo. Research on the inverse identification method of wheel-rail force based on railway vehicle’s acceleration. [Master Thesis]. Chengdu: Southwest Jiaotong University, 2020 (in Chinese)
[20] Zhu T, Wang X, Fan Y, et al. A time domain method for wheel-rail force identification of rail vehicles. Vehicle System Dynamics, 2022, 60(3): 790-809 doi: 10.1080/00423114.2020.1838562
[21] 罗金屯, 滕飞, 周亚波等. 数据驱动的高速铁路轮轨作用力反演模型. 南京大学学报(自然科学), 2021, 257(2): 299-308 (Luo Jintun, Teng Fei, Zhou Yabo, et al. A wheel-rail force inversion model for high-speed railway. Journal of Nanjing University (Natural Science), 2021, 257(2): 299-308 (in Chinese) Luo Jintun, Teng Fei, Zhou Yabo, et al. A wheel-rail force inversion model for high-speed railway. Journal of Nanjing University(Natural Science), 2021, 257(2): 299-308 (in Chinese)
[22] 曾俊玮, 季元进, 任利惠等. 卡尔曼滤波器与神经网络串行的轮胎载荷识别模型. 振动与冲击, 2023, 42(11): 262-270, 294 (Zeng Junwei, Ji Yuanjin, Ren Lihui, et al. A serial tire load identification model based on Kalman filter and neural network. Journal of Vibration and Shock, 2023, 42(11): 262-270, 294 (in Chinese) Zeng Junwei, Ji Yuanjin, Ren Lihui, et al. A serial tire load identification model based on Kalman filter and neural network. Journal of Vibration and Shock, 2023, 42(11): 262-270 + 294 (in Chinese)
[23] Iwanaga MK, Brennan MJ, Tang B, et al. Some features of the acceleration impulse response function. Meccanica, 2021, 56: 169-177 doi: 10.1007/s11012-020-01265-4
[24] 朱涛, 肖守讷, 阳光武. 一种新的时域动态载荷识别方法. 西南交通大学学报, 2012, 47(6): 968-973 (Zhu Tao, Xiao Shoune, Yang Guangwu. A new time-domain method for force identification. Journal of Southwest Jiaotong University, 2012, 47(6): 968-973 (in Chinese) Zhu Tao, Xiao Shoune, Yang Guangwu. A new time-domain method for force identification. Journal of Southwest Jiaotong University, 2012, 47(6): 968-973 (in Chinese)
[25] 赵凤遥, 张建成, 葛巍等. 基于多项式拟合初值的动载荷识别修正算法. 振动与冲击, 2021, 40(11): 211-219 (Zhao Fengyao, Zhang Jiancheng, Ge Wei, et al. Dynamic load identification correction algorithm based on polynomial fitting initial value. Journal of Vibration and Shock, 2021, 40(11): 211-219 (in Chinese) Zhao Fengyao, Zhang Jiancheng, Ge Wei, et al. Dynamic load identification correction algorithm based on polynomial fitting initial value. Journal of Vibration and Shock, 2021, 40(11): 211-219 (in Chinese)
[26] 张永强, 宋建江, 屠良尧等. 软件数值积分误差原因分析及改进办法. 机械强度, 2006, 3: 419-423 (Zhang Yongqiang, Song Jianjiang, Tu Liangyao, et al. Error analysis and improvement method when numerical integration with software. Journal of Mechanical Strength, 2006, 3: 419-423 (in Chinese) doi: 10.3321/j.issn:1001-9669.2006.03.023 Zhang Yongqiang, Song Jianjiang, Tu Liangyao, et al. Error analysis and improvement method when numerical integration with software. Journal of Mechanical Strength, 2006(3): 419-423 (in Chinese) doi: 10.3321/j.issn:1001-9669.2006.03.023
[27] 李东文, 熊晓燕, 李博. 振动加速度信号处理探讨. 机电工程技术, 2008, 198(9): 50-52, 126 (Li Dongwen, Xiong Xiaoyan, Li Bo. Research on the processing of vibration acceleration signal. Mechanical &Electrical Engineering Technology, 2008, 198(9): 50-52, 126 (in Chinese) Li Dongwen, Xiong Xiaoyan, Li Bo. Research on the processing of vibration acceleration signal. Mechanical & Electrical Engineering Technology, 2008, 198(9): 50-52 + 126 (in Chinese)
[28] 王琦, 李峰, 汤奕等. 基于物理—数据融合模型的电网暂态频率特征在线预测方法. 电力系统自动化, 2018, 42(19): 1-9 (Wang Qi, Li Feng, Tang Yi, et al. On-line prediction method of transient frequency characteristics for power grid based on physical-statistical model. Automation of Electric Power Systems, 2018, 42(19): 1-9 (in Chinese) Wang Qi, Li Feng, Tang Yi, et al. On-line prediction method of transient frequency characteristics for power grid based on physical-statistical model. Automation of Electric Power Systems, 2018, 42(19): 1-9 (in Chinese)
[29] 杨强. 基于轴箱振动的轨道不平顺估计方法研究. [博士论文]. 成都: 西南交通大学, 2013 Yang Qiang. Research on estimation methods of track irregularities based on axle-box vibration. [PhD Thesis]. Chengdu: Southwest Jiaotong University, 2013 (in Chinese)
[30] 翟婉明. 车辆−轨道耦合动力学, 第4版. 北京: 科学出版社, 2015 Zhai Wanming. Vehicle-track Coupled Dynamics, 4th edition. Beijing: Science Press, 2015 (in Chinese)
-
期刊类型引用(0)
其他类型引用(1)
计量
- 文章访问数: 386
- HTML全文浏览量: 122
- PDF下载量: 64
- 被引次数: 1
