轨道车辆垂向轮轨力时域识别对比及其机器学习修正1)

朱 涛 吴佳欣 王小瑞 肖守讷 阳光武 杨 冰

(西南交通大学轨道交通运载系统全国重点实验室,成都 610031)

引言

轨道车辆的轮轨接触力是评判其运行安全性与稳定性的重要依据,然而轮轨接触力往往难以直接获取.目前轮轨力的获取主要有直接测量法[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]

通过脉冲响应函数公式,可以得到每一时刻的脉冲响应函数值,并由此计算得出振动系统的Green函数矩阵.对于单自由度系统,其响应正演表达式如下所示

式中,z为系统的响应向量,z={z1,z2,···,zm}T;F为系统所受的载荷向量,F={F1,F2,···,Fm}T;zi和Fi分别代表第i个时刻的响应值和载荷值;m为采样点的个数;H为系统的Green 函数矩阵.式(2)可写为如下的形式

式中,hi(i=1,2,···,m)为脉冲响应函数值;Δt为时间步长.

对于单自由度系统,其Green 函数矩阵的构建方法是,通过脉冲响应函数公式计算出各个时刻的脉冲响应函数值hi(i=1,2,···,m),令矩阵主对角线以上的元素的值为0,主对角线上的元素为 Δth1,从矩阵第1 列到第m列的元素,其值均从主对角线上的 Δth1开始,第1 列为 Δth1～Δthm,第2 列为 Δth1～Δthm-1,以此类推,直到第m列,其非零元素为主对角线元素 Δth1.

与单自由度系统不同,多自由度系统通过脉冲响应函数公式计算出的各个时刻的脉冲响应函数值hi(i=1,2,···,m)均为矩阵,因此多自由度系统的Green 函数矩阵可以看做众多形同单自由度Green 函数矩阵的子矩阵的组合,各子矩阵中的元素分别为hi中对应位置的值与 Δt之积.多自由度系统的正演模型如下所示

式中,zi为第i个位置的响应列阵;Fi为第i个位置的载荷列阵;为载荷作用点i和响应测量点j之间的Green 函数矩阵;M和N分别为作用在系统上的动载荷个数和响应测量点的个数.为保证该线性方程组有解,需要满足N≥M.

在正演模型建立完成后,通过对Green 函数矩阵求逆(或求广义逆)来对动载荷进行识别,单自由度系统和多自由度系统的动载荷识别公式分别如下所示

1.2 状态空间法

对于一个n自由度线性定常系统,其动力学方程如下所示

式中,M,C和K分别为系统质量、阻尼和刚度矩阵;,和x分别为系统加速度、速度和位移向量;f为外载荷向量.

将上式转化为状态空间方程的形式,为

由文献[24]可知,一般的离散状态空间方程只能通过位移和速度对动载荷进行识别,当已知系统加速度时,需引入Newmark 积分对状态空间方程进行离散,有如下表达式

式中,Xi为状态变量,此时的状态变量由位移、速度、加速度向量组成,Xi=[xi,,]T.

式中,β 和 γ 为Newmark 积分参数,β 一般在0.25到0.5 之间选取,γ 在0 到1 之间选取[24].

由式(9)可得每一步的外载荷的求解表达式,如下所示

将式(10)展开,可得到下式

因此在已知加速度响应的情况下,通过式(10)求解外载荷在本质上就是通过式(11)～式(13)求解外载荷.

在刚度和阻尼矩阵中,并不是所有位置的响应都与外载荷的计算相关,为了提高方法的稳定性和精度,找到对外载荷fi+1的计算有影响的测点位置(即影响系数法的概念),把这些测点位置的加速度响应作为已知量,假设已知i时刻(初始时刻) 的位移、速度和加速度,根据式(11)和式(12)可以得到i+1 时刻的位移和速度,再根据式(13)即可求得i+1时刻的动载荷.对于i时刻位移和速度未知的情况,由初值引起的误差将在后文中进行讨论.该方法较适用于多自由度弹簧–阻尼–质量模型,如轨道车辆动力学模型.

文献[25] 中给出了Wilson-θ反分析法在非0 初始条件下由初值引起的误差,并引入多项式拟合趋势项对位移和速度进行了修正.由于本文采用的Newmark 积分离散的状态空间法同样需要对加速度进行数值积分来间接计算速度和位移,可能会存在相同的问题,因此同样需要分析非零初始条件下的初值误差对识别结果的影响.

为了提高方法的稳定性和精度,同样以对外载荷计算有影响的测点位置的加速度为输入,若系统为非0 初始状态,由下式所示

在1 时刻这些测点的位移实际值x1可由0 时刻的测点的位移x0、速度、加速度以及1 时刻的加速度表示,1 时刻的速度实际值可由0 时刻的速度、加速度以及1 时刻的加速度表示.

由于x0和未知,若此时仍将系统视为零初始状态,则带误差的1 时刻的这些测点的位移和速度计算值分别为

令 δ11为第1 步的位移计算误差,δ21为第1 步的速度计算误差,则有

2 初值误差分析与修正

继续计算第2 时刻测点的位移、速度的实际值与计算值

令 δ12为第2 步的位移计算误差,δ22为第2 步的速度计算误差,则有

至此已经可以发现规律,即第k步的测点位移计算误差 δ1k、速度计算误差 δ2k分别为

要确定趋势项,就需要求得一组pk使得I值达到最小,由多元函数求极值的必要条件可得

其中j=0,1,2,···,m,即

将其表示为矩阵形式,如下式所示

可以证明,方程组系数矩阵是一个对称正定矩阵,故存在唯一解.通过上式求解出pk,从而可得拟合多项式.

除了初值误差引入的趋势项之外,相关研究[26-27]也表明,对于由信号本身或其他因素所引起的积分信号的曲线趋势项,多项式拟合修正的方法同样也是有效的.

3 基于NARX 模型的动载荷识别误差修正

采用基于纯物理模型的时域法对动载荷进行识别时,难免会因为各种因素的影响(如响应观测不全,观测噪声等)造成识别误差.为了减小识别误差,参考相关文献[22,28]中提出的物理–机器学习相结合的模型思路,本节引入NARX 模型,利用其强大的非线性逼近能力,通过一段已知的输入–输出时间历程对模型进行训练,对基于物理模型时域法的动载荷识别误差进行预测,进而尝试对识别误差进行修正,从而提高动载荷识别精度.这样既避免了纯机器学习模型缺乏物理解释的问题,又减小了纯物理模型的识别误差.

NARX 模型是线性ARX (autoregressive models with exogenous inputs)模型的扩展,与静态回归模型相比,增加了多步时延的环节,并引入了非线性映射函数,因此NARX 模型具有动态特性,可以利用输入、输出的历史信息来对具有复杂映射关系的未知时间历程进行预测.对于一个单输入单输出系统,其NARX 模型可以表示为

式中,f(t) 为t时刻的系统输出;z(t) 为t时刻的系统输入;na为输出的阶数;nk为输入相对于输出延迟的阶数;nb为输入的阶数;F(·)为非线性映射函数,常用的有小波网络、S 型网络、二叉树回归模型、支持向量机等.由于小波网络具有结构简单、收敛速度快及计算精度高的优点,故本文的F(·)选用小波网络[29].

具体到本文而言,由于要使用轨道车辆轴箱响应来预测垂向轮轨力的识别误差,因此本文采用的NARX 模型的输入为轨道车辆轴箱响应,输出为垂向轮轨力的识别误差,定义识别误差为物理模型时域法的识别值减去轮轨力的真实值(可以含随机噪声),具体的NARX 模型可以表示为框图的形式,如图1 所示.

图1 NARX 模型示意图Fig.1 Schematic diagram of NARX model

采用NRMSE适应度来评价NARX 模型的预测精度,如下式所示

式中,f为目标输出;为目标输出的均值;为模型的预测输出.

在使用训练集的输入输出以及测试集的输入得到NARX 模型的识别误差预测值之后,以测试集的响应输入得到物理模型时域法的识别结果,再将该识别结果减去识别误差预测值即为修正后的动载荷识别值,其框图如图2 示.

当非线性映射函数确定之后,输出的阶数na、输入相对于输出延迟的阶数nk以及输入的阶数nb将影响NARX 模型的预测效果,参数的选取将在下一节进行讨论.

4 轨道车辆垂向轮轨力识别

根据《车辆−轨道耦合动力学》[30],车辆垂向动力学模型由车体、前后转向架和4 个轮对组成.对于车体以及前后转向架,考虑浮沉运动和点头运动;对于4 个轮对,考虑浮沉运动,总计10 个自由度,如图3 所示.在本研究中,根据该模型进行动力学建模,以轴箱加速度响应为输入,对轮轨垂向力进行识别.通过建立车辆垂向动力学微分方程,可以提取出该系统的质量M、刚度和阻尼矩阵K,C,外载荷向量如下.

图3 轨道车辆垂向动力学模型Fig.3 Vertical dynamic model of rail vehicles

式中,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 组采样点作为测试集,另取10 000 组采样点作为训练集,如图4 所示.采用Pearson 相关系数来评价方法的识别精度,如下式所示[14]

图4 第1 轮对垂向轮轨力及加速度时间历程Fig.4 Time history of vertical wheel-rail force and acceleration of the first wheelset

式中,n为采样点的数目;为识别载荷;为识别载荷的均值;Pi为外载荷;为外载荷的均值.

考虑到以往的相关研究均直接使用仿真正演的加速度响应作为模型输入,而实际上,在进行加速度响应的观测时难以避免地会存在观测噪声,为了模拟实际,在本文算例的轴箱加速度响应中添加5%的随机噪声.分别使用Green 函数法和状态空间法对轮轨垂向力进行识别,并对比识别精度与计算速度,如表1 所示,2 种方法对第1 轮对的识别结果对比如图5 所示.

表1 2 种时域法对垂向轮轨力的识别结果对比Table 1 Comparison of the two time domain methods for identifying vertical wheel-rail forces

图5 2 种方法对第1 轮对垂向轮轨力的识别结果对比Fig.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 模型的误差预测效果和输入阶数nb、输出阶数na以及输入相对于输出延迟的阶数nk的关系.需要进行说明的是,一般的机器学习研究对于模型超参数的选取常在验证集进行,而本文研究的重点主要是提供算法思路并验证其可行性,而不是超参数的选取,故没有再设置验证集,而是在测试集进行参数选取.事实上,本文测试集起到的是验证集的作用,后续如需进一步验证模型的泛化能力,再另外选取新的测试集进行测试即可.

在本例中,暂不考虑nk的影响,只通过改变输入阶数nb、输出阶数na来达到一个较好的误差预测效果.同样以第一轮对为例,设置na,nb的值以5 为步长递增,为了防止输入输出阶数过大而导致的训练发散,将其最大值设置为150,观察测试集的NRMSE适应度随着na和nb值增加的变化情况,如图6 所示.

图6 测试集预测误差NRMSE 适应度随na 和nb 值的变化关系Fig.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 所示.

表2 预测误差的NRMSE 适应度Table 2 NRMSE of the predicted error

图7 测试集的误差原始值(含5%随机噪声)和预测值Fig.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为

经修正后的轮轨力最终识别值F如图8 所示.

图8 最终识别值对比Fig.8 Comparison of the final identification values

从图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,识别精度得到了较大的提高,说明本文方法为时域法识别精度的进一步提高提供了一种思路,具有一定的工程应用价值.

在后续的研究中,还可以考虑对轨道车辆横向轮轨力以及带有部分悬挂参数非线性的轨道车辆模型的轮轨力进行识别,进一步研究本文方法在更复杂的模型中的载荷识别能力.