轨道不平顺是引起列车-轨道或列车-桥梁耦合系统振动的重要激励源之一。在车辆参数确定的情况下,轨道不平顺对列车运行的安全性与平稳性有很大的影响[1]。文献[2]探讨了列车荷载作用下不同轨道不平顺样本对桥梁结构动力响应的影响规律,结果表明,桥梁结构的动力响应随轨道不平顺样本的不同而变化,且有较大的离散性和随机性。目前,常采用单一或少数几个轨道不平顺样本作为激励源,将列车-轨道或列车-桥梁耦合系统动力响应时程的最大值(由单一轨道不平顺样本得到)或最大值的平均值(由多个轨道不平顺样本得到)作为系统安全性和平稳性的评价指标,但轨道不平顺是一个空间随机过程,一个样本仅是统计意义上随机过程的一次实现,因而由一个或少数几个轨道不平顺样本得到的系统动力响应时程可能无法反映响应的概率统计特征。此外,在列车-轨道或列车-桥梁耦合系统随机振动响应的概率统计研究中,响应样本数量的选取是一个非常关键的问题,如若响应样本数量选取过少,则无法反映各类动力响应总体的特征,也就无法建立合理的响应概率统计模型,但如果响应样本数量选取过多,又需要十分庞大的计算量。曾庆元等[3]指出,在没有实测振动响应样本的情况下,如何确定计算样本的数量才能得到更合适的随机响应统计结果是值得探讨的问题。
目前,针对列车-轨道(或桥梁)耦合振动的随机性问题和系统响应的概率统计特征问题,许多学者已开展了一定研究。夏禾等[4]取20个轨道不平顺样本计算了列车通过单跨箱梁时桥梁和列车响应的最大值、平均值、标准差及变异系数,结果表明,20次计算值的离散性很大。余志武等[5]考虑车辆荷载及一系、二系竖向弹簧刚度与阻尼的随机性,对列车-桥梁竖向耦合振动响应开展了研究,计算得到了车桥振动响应的均值、标准差和时变概率密度演化分布。此外,谈遂等[6]基于某重载铁路桥梁的现场试验数据,对不同轴重列车过桥时桥梁动力响应的随机性进行了分析,并结合3倍标准差原理对桥梁动力响应最大值进行了估计。李永乐[7]针对风-列车-桥梁耦合系统的振动特性,指出采用基于单样本分析得到的响应最大值来评价列车和桥梁的振动水平不合适,并提出将桥梁的静风位移、静车位移和考虑峰值因子的动位移均方差叠加作为桥梁位移响应的评价指标,将考虑峰值因子的加速度均方根作为桥梁加速度响应的评价指标。Xin等[8]对列车-桥梁耦合系统的不确定性和敏感度开展了研究,揭示了桥梁参数和轨道不平顺的随机性对列车-桥梁耦合系统响应的变异性有较明显的影响,而不同因素之间的交叉影响作用也不容忽视。朱志辉等[9]基于Monte-Carlo方法,将车速和轨道不平顺作为随机变量,随机选取100个车速和轨道不平顺的组合样本,得到了钢桥最不利疲劳部位的等效应力范围和年疲劳损伤系数的概率分布。
目前对车-桥系统或车-轨系统振动随机性的处理主要有蒙特卡洛(Monte-Carlo)法、虚拟激励法、概率密度演化法、随机模拟法等[10]。其中,蒙特卡洛法具有简单直观、有较大样本容量、适用性较强等优点,且常被用于检验其他分析方法的有效性[11-12]。为可靠地评估列车运行的安全性与平稳性,笔者采用多体动力学软件Simpack建立列车-轨道耦合模型,将通过引入三角级数法模拟得到的轨道不平顺作为随机输入激励,基于Monte-Carlo方法计算得到了列车在行驶过程中加速度振动响应的样本序列,并结合高斯混合(Gaussian Mixture Model, GMM)模型,采用期望最大化(Expectation Maximization, EM)算法对概率模型参数进行最大似然估计,从而对列车加速度振动响应最大值分布的统计规律以及响应样本数量的选取进行研究。
以中国某高速列车动车车辆为原型,运用多体动力学分析软件Simpack建立相应的列车多体动力学精细化模型。模型中忽略车体、转向架、轮对、轴箱转臂等部件的弹性变形将其视为刚体,列车的一、二系弹簧和减振器采用弹簧-阻尼元件进行模拟,并考虑其作用点之间的空间距离以及刚度和阻尼的非线性特性,其中,非线性弹簧、减振器主要有一系垂向减振器、横向减振器、抗蛇行减振器和横向止挡等。模型中采用s1002轮对,钢轨选用T60 kg/m型钢轨。采用准弹性接触模型计算轮轨接触力,轮轨法向力采用Herz理论计算,蠕滑力则采用FASTSIM算法求得。根据相关文献,判断一个复杂的多体动力学模型是否建立正确,主要依据静平衡计算名义力时系统的最大残余加速度是否小于0.01 m/s2。对列车-轨道耦合系统多体动力学进行静平衡名义力计算时,最大残余加速度为8.389 84×10-6 m/s2,小于0.01 m/s2,因此,列车-轨道耦合系统多体动力学模型建立正确。列车的多体动力学模型如图 1所示,列车模型主要参数如表 1所示。
轨道不平顺是列车-轨道耦合系统振动重要的激励源,是引起车辆产生振动的主要原因之一。轨道不平顺主要可划分为水平不平顺、轨距不平顺、垂向不平顺、方向不平顺4类。轨道不平顺具有随机性,可视为由不同幅值、波长和相位角的波叠加而成,通常用功率谱密度函数来表示。在求解动力响应的过程中,需要将轨道不平顺功率谱转换为时域样本或空间样本[13]。
基于中国高速铁路总体技术条件,建议对列车进行平稳性分析时使用德国高速线路轨道谱[14],且高速铁路试验段轨道谱的高低不平顺在30 m波长以上的平顺性基本与德国高速低干扰谱接近[15]。因此,选用德国高速低干扰谱,采用三角级数法对轨道不平顺序列进行模拟。在假设轨道不平顺为平稳遍历的高斯白噪声的前提下,轨道不平顺可看作是不同三角级数的叠加,可通过式(1)得到。
式中:f(x)为模拟得到的轨道不平顺序列;S(ωk)为功率谱密度函数,垂向不平顺单位为m2/rad/m,水平不平顺单位为1/rad/m;ωk为轨道不平顺的空间频率,rad/m;φk为第k个频率的相位,是独立均布于0~2π的随机数。
用三角级数法模拟得到的轨道不平顺功率谱密度(Power spectral density,PSD)与目标谱的吻合情况如图 2所示,由图 2可知,模拟的功率谱与目标谱吻合较好。
以随机轨道不平顺作为输入激励,得到列车以200 km/h速度行驶时的加速度响应时程曲线,如图 3所示。列车加速度测点位置根据《铁道车辆动力学性能评定和实验鉴定规范》(GB 5599—1985)规定设定于转向架中心上方横向1 m的车体地板上。
高斯混合模型广泛应用于统计分析领域[16],其作为一种统计模型,多用于构建概率密度函数,能较好地刻画参数空间中数据的分布及其特征,既具有非参数化方法的灵活性,又保持了参数化方法的精确性。高斯混合模型采用有限个特定概率分布密度函数的加权组合来拟合复杂的概率分布模型,通过选择混合分量的类型和个数,可逼近任何连续的概率分布密度函数。高斯混合模型(GMM)由高斯(正态)分布的加权组合得到,其概率密度函数为
式中:参数μzk、σzk分别为第k个高斯成分的均值和方差;πk是随机变量x取至第k个高斯成分的权重系数,表示观察样本属于第k个高斯分量所代表的聚类的相对率,应满足公式$0 \le {\pi _k} \le 1, \sum\limits_{k = 1}^M {{\pi _k}} = 1$,k=1, 2, 3,…,M[17]。
高斯混合模型可以逼近任何连续的概率分布函数,选择合适的权重系数是得出模型类型数量与模型参数的关键。高斯混合模型是一种“软分类聚类”,是基于假设数据集是由一个潜在的混合概率分布产生的,其中每个高斯分量都表示一个不同的聚类。首先通过估计样本数据集的混合概率密度,然后计算样本源中单个样本之于各个高斯分量的后验概率,最后将单个样本分配到后验概率最大的高斯分量所在的聚类组,从而得到样本数据集所服从的高斯混合分布[18]。
高斯混合模型(GMM)的期望为
式中:πk为随机变量x取至第k个高斯成分的权重系数;N(xi; μk; ∑k)为第k个类别的条件概率密度;μk、∑k分别为该高斯分量的均值和协方差矩阵。
对于高斯混合模型的参数,可用期望最大化(EM)算法进行迭代估计[19]。估计步骤为
1) E步
2) M步
式中:τij(k)为第k个高斯分项生成的期望;πi(k+1)、μi(k+1)、∑i(k+1)为对应的参数估计值。E步计算期望,利用对隐藏变量的现有估计值,计算其最大似然估计值;M步为期望最大化,最大化E步上求得的最大似然值来计算参数的值。
3) 收敛条件
不断迭代E步与M步,至似然函数的变化量小于误差值esp或迭代次数k≤K时,迭代结束,否则返回E步。随着迭代次数的增加,迭代误差越来越小,似然函数不断收敛。可接受的迭代误差esp=2×10-16,最大迭代次数K=500。似然函数为
综上所述,EM是一种迭代算法,也是一种聚类算法,它可以通过迭代求出高斯混合模型的参数,并将样本源中的单个样本通过迭代收敛性进行聚类。高斯混合模型聚类通常采用贝叶斯信息准则(BIC)选择模型,模型的BIC值越大,该模型就越符合实际。
假设列车振动加速度响应的最大值服从高斯混合分布,采用期望最大化算法对该概率模型参数进行最大似然估计,再进行拟合度检验。
为了对比得到的高斯混合模型概率密度函数与由列车加速度响应最大值的样本序列得到的频率直方图的拟合效果,采用拟合优度(Adjusted R2)与均方根误差(Root Mean Squared Error,RMSE)两项指标来对概率密度分布曲线的拟合优劣程度进行评价。拟合优度用于评价概率密度分布曲线与直方图之间的相似程度,该值越接近于1,表示拟合程度越好;均方根误差用于评价概率密度分布曲线与直方图之间的偏离程度,该值越接近于0,表示偏离程度越小,拟合程度越好。
样本量是指从总体中抽取的样本元素的总个数,样本量的大小是选择检验统计量的一个重要要素。由抽样分布理论可知,在大样本条件下,如果总体为正态分布,则样本统计量服从正态分布;如果总体为非正态分布,则样本统计量渐近地服从正态分布[20]。
样本量的计算公式为
式中:n为样本量;α为显著水平;Zα/2为置信区间对应的标准分数;E为估计误差,一般小于0.1;σ为标准差,一般为0.5。
在确定样本量时,取α为0.05,则置信度为95%,经查表,Zα/2为1.96;假定的估计误差为0.05,则最小样本量为n=384。在满足样本最小容量的情况下,增加样本量有助于增加检验统计的精度,提高可靠性。
基于Monte-Carlo方法,分别取400、500、600、700、800、900和1 000个随机生成的轨道不平顺样本,计算得到列车的加速度响应时程样本,并将列车竖向和横向加速度响应的最大值作为随机变量,采用EM算法进行聚类分析,比较聚类为1~4类的BIC值,并选择BIC值最大的一组参数,得出某个确定样本容量下列车加速度响应最大值所服从的高斯混合模型。接着,对比不同样本容量下列车加速度响应最大值的高斯混合模型与其相应频率分布直方图的拟合优度及均方根误差,从中选取合适的高斯混合模型概率密度函数(对应的样本容量记为Nrep)作为代表该车速下列车加速度响应最大值的概率密度模型,此时的样本容量Nrep作为代表该车速下列车加速度响应最大值概率统计特征的最优样本数量。
当列车以200 km/h车速行驶时,以1 000个列车竖向加速度响应样本为例,表 2为列车竖向加速度最大值不同聚类个数的BIC值。由表 2可得,当聚类个数N=2时,BIC值最大。表 3为N=2时模型的参数估计值。将表 3中的参数带入到式(2)中,即可得到列车竖向加速度最大值所服从的高斯混合模型概率密度函数,如式(10)所示。图 4为得到的高斯混合模型概率密度函数与相应样本数量下列车竖向加速度响应最大值的频率分布直方图的对比,从图 4可知,其拟合效果较好。
分别取400、500、600、700、800、900和1 000个随机生成的轨道不平顺样本,计算得到车速为200 km/h时不同样本容量下的列车竖向加速度响应最大值的高斯混合模型概率密度分布与其频率分布直方图的拟合优度及均方根误差,如表 4和图 5所示。由表 4和图 5可知,当样本数量为400~700时,列车竖向加速度响应最大值的高斯混合模型概率密度分布与其频率分布直方图的拟合优度随着样本数量的增加逐渐变大,在样本数量达到700后趋于稳定;而均方根误差随着样本数量的增加逐渐减小,在样本数量达到700后,波动较小,趋于稳定。因而,车速为200 km/h时,代表列车竖向加速度响应最大值概率统计特征的最优样本数量Nrep=700。图 6为不同样本数量下的列车竖向加速度最大值所服从的高斯混合模型的概率密度曲线。
当车速分别为100、150、200 km/h时,不同样本数量下列车竖向加速度响应最大值的高斯混合模型概率密度分布与其频率分布直方图的拟合优度及均方根误差如图 7所示。从图 7中可以看出,在样本数量达到700后,3种不同车速下的列车竖向加速度响应最大值的高斯混合模型概率密度分布与其频率分布直方图的拟合优度和均方根误差均波动较小,趋于稳定。因而当车速分别为100、150、200 km/h时,代表列车竖向加速度响应最大值概率统计特征的最优样本数量Nrep均可取为700。
当列车以200 km/h车速行驶时,仍以1 000个列车横向加速度响应样本为例,表 5为列车横向加速度最大值不同聚类个数的BIC值,由表 5可得,当聚类个数N=3时,BIC值最大。表 6为N=3时模型的参数估计值。将表 6中的参数带入式(2)中即可得到列车横向加速度最大值所服从的高斯混合模型的概率密度函数,如式(11)所示。图 8为得到的高斯混合模型的概率密度函数与相应样本数量下列车横向加速度响应最大值的频率分布直方图的对比,由图 8可知其拟合效果较好。
分别取400、500、600、700、800、900和1 000个随机生成的轨道不平顺样本,计算得到车速为200 km/h时不同样本容量下列车横向加速度响应最大值的高斯混合模型概率密度分布与其频率分布直方图的拟合优度及均方根误差,如表 7和图 9所示。由表 7和图 9可知,列车横向加速度响应最大值的高斯混合模型概率密度分布与其频率分布直方图的拟合优度和均方根误差在样本数量达到800后波动较小,趋于稳定。因而,车速为200 km/h时,代表列车横向加速度响应最大值概率统计特征的最优样本数量Nrep=800。图 10为不同样本数量下列车横向加速度最大值所服从的高斯混合模型的概率密度曲线。
当车速分别为100、150、200 km/h时,不同样本数量下列车横向加速度响应最大值的高斯混合模型概率密度分布与其频率分布直方图的拟合优度及均方根误差如图 11所示。从图 11中可以看出,在样本数量达到800后,3种不同车速下列车横向加速度响应最大值的高斯混合模型概率密度分布与其频率分布直方图的拟合优度和均方根误差均波动较小,趋于稳定。因而当车速分别为100、150、200 km/h时,代表列车横向加速度响应最大值概率统计特征的最优样本数量Nrep均可取为800。
列车分别以车速100、150、200 km/h行驶1 000 m,基于Monte-Carlo方法,取1 000个随机生成的轨道不平顺样本计算得到列车的加速度响应时程样本。通过统计列车的竖向加速度最大值及横向加速度最大值,得到列车不同车速下的加速度最大值箱型图,如图 12所示。由图 12可知,列车的竖向加速度最大值和横向加速度最大值均随着车速的增加逐渐变得离散;列车的竖向加速度最大值和横向加速度最大值的均值随着车速的增加而增加。此外,列车竖向加速度最大值的均值都大于中位数,而列车横向加速度最大值的均值在车速为100、150 km/h时大于中位数,在车速为200 km/h时的均值小于中位数。
通过对列车的竖向加速度最大值样本及横向加速度最大值样本进行分析,进一步得到不同车速下列车加速度最大值的概率密度分布图,如图 13所示。从图 13中可知,列车竖向加速度最大值和横向加速度最大值的概率密度曲线均沿着横坐标向右移动。另外,与车速为100、150 km/h相比,车速为200 km/h时列车竖向加速度与横向加速度最大值的概率密度曲线的分布范围均更大;而车速为100 km/h时列车竖向加速度与横向加速度最大值的概率密度曲线最为高耸。这表明随着车速的增加,列车加速度响应最大值分布的离散性增强。
由于列车振动会使乘车人员产生不适感或疲劳,因而引入平稳性指标来度量乘客乘车时的舒适程度。参考高速铁路客车动力学性能评定资料,中国车体振动加速度的平稳性标准界限值可取为:竖向振动加速度av≤1.3 m/s2,横向振动加速度ah≤1.0 m/s2。
在运行距离为1 000 m的情况下,选择列车运行车速为100、150、200 km/h的工况,得到列车加速度最大值随机变量的累积分布函数(Cumulative distribution function, CDF)曲线,如图 14~图 16所示。从图中可以发现,列车在运行车速为100、150、200 km/h时的竖向加速度与横向加速度均满足平稳性要求。
1) 采用高斯混合模型能够有效地拟合列车加速度响应最大值的分布规律。
2) 当车速分别为100、150、200 km/h时,列车竖向加速度响应的样本数量达到700时可较好地代表列车竖向加速度响应最大值的概率统计特征;而列车横向加速度响应的样本数量达到800时能较好地代表列车横向加速度响应最大值的概率统计特征。
3) 列车的竖向加速度最大值和横向加速度最大值均在车速为100 km/h时分布最为集中;整体来讲,随着车速的增加,列车加速度响应最大值分布的离散性增强。