基于HHT能量和最大Lyapunov指数的蛇行分类方法

王 敏，宁 静，赵 飞，李艳萍，陈春俊

(西南交通大学 机械工程学院，成都 610031)

0 引言

蛇行运动是车辆动力学系统的核心问题之一，车辆在运行过程中发生收敛较慢的小幅蛇行或者剧烈的大幅蛇行都会严重地影响车辆运行安全，因此对车辆小幅蛇行和大幅蛇行的在线监测都至关重要。学者们通过大量的理论研究总结出轮轨参数和悬挂系统参数对车辆系统蛇行运动的影响规律[1-5]。在这些理论研究的成果上，更多的学者致力于探索如何将这些成果应用到高速列车蛇行监测领域。宋兴武[6]提出横向位移峰值方法(LMP)，通过计算轮对和构架的横向加速度确定轮对横移量的峰值来识别车辆的蛇行运动。朴明伟等[7]发现当车辆系统服从超临界分岔时，随着运行速度的增加，转向架会发生小幅值蛇行，且随着速度的增加，转向架的蛇行幅值变大，最终会出现大幅蛇行失稳。蔡里军[8]总结了大量的实测数据认为当构架横向加速度峰值连续6次达到或超过2 m/s2，车辆发生了蛇行运动。方明宽[9]等考虑到构架横向加速度数据的非平稳性，提出NKJADE方法对多传感器横向加速度数据进行特征融合，从而对车辆的正常运行，小幅蛇行，大幅蛇行状态进行识别。Sun[10]发现当车辆出现蛇行运动时，构架与车体的加速度信号存在相位延迟，通过对车体和构架的横向加速度信号进行互相关分析，将互相关指标作为判别车辆蛇行失稳的阀值。曾元辰、赵飞[11-12]等发现车辆出现蛇行运动时，轮对、构架和车体的横向加速度都具有较强的周期性，提出利用信号周期性特点来识别车辆的蛇行运动。崔万里[13]通过样本熵理论以及等距映射算法(ISOMAP)对高速列车构架横向加速度信号进行特征提取，然后通过LS-SVM进行车辆小幅蛇行异常识别，达到了较好的识别效果。叶运广[14]通过MEEMD进行特征提取与使用LS-SVM对车辆的运行状态进行特征识别亦取得很好的效果。冉伟[15]因而提出了一种基于EEMD-SVD-LTSA的特征提取框架，用于识别高速列车小幅蛇行运动的演变趋势。陈杨[15]通过建立道岔模型，研究高速列车在小幅蛇行的状态下通过道岔的演变规律。宁云志[15]提出基于1D-CNN和CGAN的预测方法对小幅蛇行进行预警，实现了蛇行失稳的提前预警。王晓东[18]利用1D CNN自适应地对构架横向加速度数据进行特征提取，使用LSTM对车辆的正常、小幅蛇行和大幅蛇行状态进行分类，取得了很好的分类效果。但是目前，列车实际运行过程中对蛇行运动的监测还存在以下不足：

问题一：现有的蛇行监测标准将构架横向加速度信号峰值连续6次达到8 m/s2时视为大幅蛇行，但实际过程中即使构架横向加速度没有达到现有的监测标准也会出现大幅的振动。此外，现有标准没有定量反映蛇行失稳对车辆的影响程度，仅用一个固定值来评判是否产生蛇行，具有一定的局限性。

问题二：某些车辆受到部分轨道不平顺周期性的干扰，虽然发生了明显的谐波振动，但只是运行中的极少时刻，且会快速收敛，不会形成稳定的周期蛇行运动，对运行安全性没有较大的影响[19]，将此类状态认为快速蛇行收敛。

本文针对以上问题，提出HHT(Hilbert-Huang transform)能量法从信号频域主频的大小与频谱的集中性以及频率值等方面来判断高速列车是否存在蛇行运动，对两个蛇行运动状态(小幅蛇行和大幅蛇行)的蛇行程度的大小进行定量的评估。利用最大Lyapunov指数进一步表征信号的周期性从而将快速蛇行收敛和蛇行运动(小幅蛇行和大幅蛇行)区分开来。最后通过HHT能量法与最大Lyapunov指数相结合，从高速列车实际运行的需求出发，对车辆系统进行正常运行、小幅蛇行、快速蛇行收敛、大幅蛇行的定性识别和蛇行程度大小的定量分析，以实现对蛇行运动的具体监测。

1 仿真模型建立

1.1 高速列车动力模型

本文使用SIMPACK软件建立国内某型号高速列车的整车动力学模型，该车辆的基本参数如表1所示。该车辆包括1个车体、2个构架、4个轮对、8个轴箱、一系悬挂系统(一系钢簧、一系垂向减振器、转臂式轴箱定位装置)和二系悬挂系统(空气弹簧、二系垂向减振器、二系横向减振器、抗蛇行减振器、抗侧滚扭杆、牵引拉杆、横向止挡)。在建立仿真模型时，忽略了车辆系统各部件的弹性变形，将车体、构架、轮对、轴箱等部件假定为刚体。车体、构架和轮对有6个自由度(纵向、横向、垂向、侧滚、点头、摇头)，轴箱只有1个点头自由度。因此，该车辆动力学模型共有50个自由度。

表1 高速列车动力学参数表

1.2 建立仿真工况类型

关于轮对踏面的选择，踏面类型对高速列车系统极限环分岔影响显著，为了模拟高速列车的不同运行状态，同时使高速列车的运行状态更符合实际，本文分别选用了在一个旋修周期(2.5×105km)内3种不同行驶里程的S1002G实测踏面。

踏面1至踏面3对应行驶的里程依次增大，踏面的磨损程度依次增加。其中，踏面1为列车行驶5×104km的行驶里程左右时，测量轮对外形轮廓获得的实测磨损踏面；踏面2为列车行驶 10×104km 的行驶里程左右时，测量轮对外形轮廓获得的实测磨损踏面；踏面3为列车行驶15×104km 的行驶里程左右时，测量轮对外形轮廓获得的实测磨损踏面。

关于抗蛇行减振器的参数设定，抗蛇行减振器的参数对高速列车系统的临界速度和分岔类型也有很大影响。本文为充分模拟高速列车的不同运行状态，使用3种不同阻尼特性的抗蛇行减振器。抗蛇行减振器的阻尼特性曲线如图1所示，从抗蛇行减振器1到抗蛇行减振器3对应的卸荷力和卸荷速度依次增加。使用不同行驶里程的S1002G实测踏面和不同阻尼特性的抗蛇行减振器模拟高速列车(速度：300～400 km/h)3种工况：

图1 抗蛇行减振器阻尼系数图

1)工况1：使用踏面1和抗蛇行减振器1仿真，模拟正常运行；

2)工况2：使用踏面2和抗蛇行减振器2仿真，模拟小幅蛇行；

3)工况3：使用踏面3和抗蛇行减振器3仿真，模拟大幅蛇行和快速蛇行收敛。

根据以上工况，计算得到各工况下的极限环分岔图，如图2所示。工况1的高速列车动力学模型分岔类型为超临界分岔，在此工况下车辆模型在300～400 km/h区间运行时，车辆处于正常运行状态，如图3所示。工况2的高速列车动力学模型分岔类型同样为超临界分岔，在此工况下车辆模型在300～400 km/h区间运行时，车辆的轮对处于小幅振动状态，此时车辆处于小幅周期蛇行状态，随着速度增加，横向振动幅值也逐渐增大，如图4所示。工况3下高速列车动力学模型分岔类型为亚临界分岔，在此工况下，车辆模型在300～400 km/h区间运行时，车辆系统极限环处于不稳定区域，受轨道不平顺的影响，轨道不平顺激扰大时，车辆系统平衡点收敛为稳定极限环值产生大幅蛇行运动，如图5所示。当轨道不平顺激扰不足以让车辆系统产生稳定的极限环时，车辆系统可能会出现短暂的蛇行运动随后车辆系统平衡点收敛为0，无法形成稳定的周期运动，即快速蛇行收敛，如图6所示。

图2 车辆系统极限环分岔图

图3 正常运行

图4 小幅蛇行

图5 大幅蛇行

图6 快速蛇行收敛

通过建立上述工况，确定了高速列车动力学模型4种运行状态的仿真条件。为构建大量的仿真数据集，将车辆运行速度分为(300 km/h、320 km/h、340 km/h、360 km/h、380 km/h、400 km/h)6个速度级，对上文用到的轨道不平顺分别设置了0.5、1.0、2.0三种幅值比例系数，扩充轨道不平顺数据集。每种运行状态对应6个速度级，3种幅值比例轨道不平顺，共18组仿真条件，4种运行状态共对应72组仿真数据。对每组仿真得到的构架横向加速度信号随机取20个样本，每个样本段取样时长为4 s，采样频率为250 Hz，总计1 440个样本段的仿真数据集，供后续实验使用。

2 理论方法

2.1 HHT能量理论

大量的研究表明蛇行信号和正常信号的主频不同，基于这一特点，通过研究信号所有频率成分在频域段内的变化，可以更加准确地找到正常信号和蛇行信号在频域内的不同特征。再结合信号的时域特征进行分析，可以准确地得知信号的时频信息，从而更加全面地分析信号特征。孙新建[20]通过HHT能量方法对非平稳的爆炸震动信号进行分析，对希尔伯特边际谱进行频域积分计算HHT能量，研究结构体的爆炸震动损伤。HHT能量法结合了信号频域主频的大小与频谱的集中性以及频率值和时域幅值的相关信息，因此本文采用此方法研究蛇行运动，

首先对蛇行信号和正常信号的主频及频域范围进行分析，确定蛇行频域主要范围。然后通过HHT边际谱统计正常信号和蛇行信号在特定频率段范围内的频域信息，最后结合信号的时域幅值和频域信息，通过HHT能量法全面反映蛇行信号和正常信号不同的时频特征，从而对车辆运行状态准确的识别。

HHT具体变换原理如下：对原始信号x(t)进行经验模式分解得到n个IMF(intrinsic mode function)信号分量：

(1)

式中，ci为信号的基本模式分量，rn称为残余函数。对ci进行Hilbert变换得到：

(2)

式中，τ为时间延长间隔。构造解析信号zi(t)：

zi(t)=ci(t)+jGi(t)=ai(t)ejΦi(t)

(3)

式中，幅值函数为ai(t)，相位函数为Φi(t)，ci分量信号的Hilbert谱为Hi(ω，t)：

Hi(ω，t)=RP[ai(t)ejΦi(t)dt]

(4)

式中，RP表示取实部。对Hi(ω，t)进行时域上的积分得到边际谱hi(ω)：

(5)

式中，T为信号时间长度。对hi(ω)(i=1，2，…，n)边际谱分析，对主频在2 Hz以上的边际谱进行叠加求得最终的边际谱h(ω)：

(6)

式中，hk(ω)(k=1，2，...j)为主频在2 Hz以上的边际谱，在频域内对h(ω)进行积分得到最终的HHT能量值。

(7)

2.2 最大李雅普诺夫指数理论

蛇行信号和正常信号除了时频上的差异，其周期性特征也有着明显的区别，所以有必要对信号的周期性特征进行分析。李雅普诺夫方法是研究运动稳定性的重要方法，对于非线性动力学系统，李雅普诺夫指数可以从状态变量的时间序列中计算出来。最大李雅普诺夫指数为相空间相邻轨迹的平均指数发散率的数值特征，常用来判定系统的混沌特性，因此其可以用来量化信号的周期性。最大李雅普诺夫指数计算步骤如下。

对于时间序列{x(t)}构造n维空间Rn，T为延迟时间，T=k△t(k=1，2，…)，△t为时间间隔即：

X(t)=[x(t)，x(t+T)，...，x(t+(n-1)T)]T

(8)

N维相空间中的某一时刻，两条邻近轨迹之间的距离可以分解在n个不同的方向，这个n个不同方向的上的距离增长率是不同的，每一个增长率就是一个李雅普诺夫指数。系统可以写成n个自治一阶微分方程组的形式即：

(9)

取两条邻近轨迹L1和L2，起始点分别为x0和y0两起始点的距离为d0=y0-x0，经过△t时间后分别运动到x1和y1，此时距离为d1=y1-x1，如此循环下去经过m△t后得到m个di(i=1，2，…，m)，di=yi-xi，最大李雅普诺夫指数λ1即：

(10)

3 基于仿真数据的蛇行状态研究

3.1 研究步骤

由于实测数据较少，本文利用仿真数据进行分析，再通过实测数据进行验证，通过HHT能量法和最大Lyapunov指数法对蛇行状态进行分类且定量反映蛇行程度的大小，其具体操作如下：

1)仿真得到不同运行状态下的构架横向加速度信号数据，采用0.5～10 Hz对构架横向加速度进行滤波，提取时间长度为4 s的信号进行分析。

2)通过对大量的仿真数据进行EMD(empirical mode decomposition)分析，得出蛇行信号主频和能量频带基本在2 Hz以上。

3)对主频在2 Hz以上IMF(intrinsic mode function)分量信号的Hilbert边际谱进行叠加求和。

4)在频域内对最终得到的Hilbert边际谱进行积分得到最终的HHT能量值。

5)利用SVM(support vector machines)分类得到HHT能量阈值，根据HHT能量值将正常运行、小幅蛇行、大幅蛇行区分开，且通过HHT能量值的大小反映蛇行程度的大小。

6)利用最大Lyapunov指数对高速列车不同运行状态信号的周期性进行分析。

7)利用SVM分类得到最大李雅普诺夫指数阈值，根据最大Lyapunov指数值将快速蛇行收敛与蛇行运动(小幅蛇行和大幅蛇行)区分开来。

3.2 HHT能量分析

对图3、4、5、6的构架横向加速度信号进行示例分析。如表2所示，正常运行信号的主频和能量频带分布比较均匀，主频2 Hz以上IMF信号分量能量占比62%，2 Hz以下IMF信号分量能量值占比38%。如表3和表4所示，蛇行信号(小幅蛇行、大幅蛇行)主频高的IMF分量能量占比高，蛇行频率主频和能量频带基本在2 Hz以上，信号分量能量值均占比达到90%以上。如表5所示，快速蛇行收敛信号因为包含蛇行特征，IMF分量主频在2 Hz以上的占比也高达80%。对主频在2 Hz以上IMF分量信号边际谱进行叠加，得到最终的Hilbert边际谱，最后在频域内对Hilbert边际谱进行积分得到最终的HHT能量值。

表2 正常运行(2 Hz以上HHT能量所占百分比为62.1%)

表3 小幅蛇行(2 Hz以上HHT能量所占百分比为95.87%)

表4 大幅蛇行(2 Hz以上HHT能量所占百分比为99.5%)

表5 快速蛇行收敛(2 Hz以上HHT能量所占百分比为82.9%)

对仿真数据集(1 440个样本)的正常运行、小幅蛇行和大幅蛇行3类数据的HHT能量值进行分析，利用SVM方法计算得出其分类阈值，分别为0.82和2.34，如图7所示。当信号HHT能量值小于0.82时，将其视为正常运行，当信号HHT能量值大于0.82小于2.34时视为小幅蛇行，当信号HHT能量值大于2.34时视为大幅蛇行。HHT能量值的大小反映了蛇行程度的大小，即HHT能量值越大蛇行程度越剧烈，HHT能量值越小蛇行程度越小。如图8所示，在此3类数据上加入快速蛇行收敛数据，发现快速蛇行收敛HHT能量值和小幅蛇行、大幅蛇行HHT能量值重叠，表明HHT能量值指标无法区分这3种运行状态，因此需要补充新的指标将其区分。

图7 三类仿真数据HHT能量统计图

图8 四类仿真数据HHT能量统计图

大幅蛇行和小幅蛇行的最大Lyapunov指数值信号周期性较强，正常运行和快速蛇行收敛的信号周期性较弱。因此可通过最大Lyapunov指数值将快速蛇行收敛和正常运行归为周期性较弱的一类，将小幅蛇行和大幅蛇行归为周期性较强的一类。

对仿真数据集(1 440个样本)的正常运行、小幅蛇行、大幅蛇行、快速蛇行收敛数据的最大Lyapunov指数进行分析，如图9所示。利用SVM方法计算得出其分类阈值10.26，李雅普诺夫指数值低于10.26时，将其认为属于正常运行和快速蛇行收敛一类，最大Lyapunov指数值大于10.26时将其认为小幅蛇行和大幅蛇行一类。

图9 仿真数据最大Lyapunov指数统计图

通过最大Lyapunov指数法和HHT能量法结合，将正常运行、小幅蛇行、大幅蛇行、快速蛇行收敛区分开来。当信号HHT能量值小于0.82时，将其视为正常运行，当信号HHT能量值大于0.82小于2.34且最大Lyapunov指数值小于10.26时视为小幅蛇行，当信号HHT能量值大于2.34且最大Lyapunov指数值小于10.26时视为大幅蛇行，当信号HHT能量值大于0.82且最大Lyapunov指数值大于10.26时视为快速蛇行收敛。此方法既考虑了蛇行的周期性特点，也考虑了蛇行信号的频域特性包括频率主频的大小、频谱的集中性以及频率值等特点，从而可以定性识别车辆系统不同的运行状态，且能够定量表示蛇行信号的能量值从而反映蛇行程度的大小。

4 实测数据验证

为验证本文所提方法的正确性，采用实测数据对其进行验证。图10为某高速列车速度在300～400 km/h区间的部分实测构架横向加速度数据，时长为1 220 s，单个样本数据长度为4 s，样本总数为305。

图10 实测构架横向加速度信号

用HHT-最大Lyapunov指数方法对实测数据进行分析，图11为实测数据的HHT能量值，图12为实测数据的最大Lyapunov指数图。图13～15为局部信号监测放大图。

图11 实测数据HHT能量图

图12 实测数据Lyapunov指数图

图13 小幅-大幅监测图

从图13的识别结果中可以观察到，704～708 s这一段信号，传统的方法会将其识别为正常运行，但这段信号出现了明显的周期性特征，本文所提方法将其识别为小幅蛇行状态，这段小幅蛇行识别可以起到预警作用。708～712 s这一段信号，其加速度峰值未达到目前所使用的蛇行报警标准(连续6个周期加速度峰值达到8 m/s2)，传统的方法会将其识别为正常运行，但通过观察可以看出，该信号峰值(6～8 m/s2之间)接近目前所使用的蛇行报警标准，且出现了明显的周期性，本文所提方法将其识别为大幅蛇行状态，故本文认为其识别结果合理。

从图14的识别结果可以看到712～716 s这一段信号，传统的方法无法识别其小幅蛇行状态，但其具有明显的周期性特征，716～720 s这一段信号本方法识别为快速蛇行收敛，从图中可以看出其周期性特征逐渐消失，证明了本方法的合理性。

图14 小幅-快速蛇行收敛监测图

从图15的识别结果看到，1 084～1 092 s这一段信号从正常运行变为小幅蛇行的状态。其信号周期性特征逐渐变强，幅值不断上升，车辆系统出现小幅蛇行。传统的方法会将这一段信号整体识别为正常运行，本方法可以准确地监测这一变化过程。

图15 正常-小幅监测图

本文所提方法准确地识别了车辆蛇行运动的不同状态，能够及时采取措施预防列车发生严重的蛇行运动，保证列车运行安全。综上，本文利用实际数据验证了方法具有可行性，能够在车辆运行过程中实现更为准确的在线监测。

5 结束语

现有的高速列车蛇行失稳标准仅针对大幅蛇行且单一的从信号时域的角度分析。本文将HHT能量值和最大Lyapunov指数值两个指标进行了结合，综合考虑了信号的时域、频域和周期性特点，以此作为区分高速列车不同蛇行状态的参考指标。本方法重点并非确定不同蛇行运行状态的准确阈值，由于仿真数据的局限性，本文通过仿真数据所确立的阈值标准存在一定的误差。阈值的准确性问题可通过台架试验或线路实测等方式对真实数据进行补充，再结合迁移学习等机器学习算法，提高实际监测中阈值的准确性。