高速铁路中，桥梁占线路比重通常较大，京沪高铁桥梁工程占比达86.5%，桥上行车安全成为最受关注的问题之一. 桥上行车安全常采用耦合振动方法进行研究，但受轨道几何不平顺等因素的影响，车-桥耦合系统响应具有明显的随机性. 已有研究表明，同一列车即使以相同速度在同一桥梁上在行驶，车-桥耦合系统的动力响应变异系数也可多达30%^[1]. 因此，研究车-桥耦合系统响应的可靠性具有重要意义.

车-桥耦合振动的求解方法主要有随机振动方法和蒙特卡罗模拟（monte carlo simulation, MCS）法. 随机振动方法是通过求解系统的统计响应以分析系统的可靠性，而MCS法是时域法，通过对大量样本进行统计得到系统响应的统计值.

目前，在车-桥耦合系统中应用的随机振动方法主要有协方差分析法、虚拟激励法及概率密度演化法等. 晋智斌^[2]将轨道不平顺激励转化为白噪声激励，提出了车-桥时变系统竖向随机振动分析的方差递推算法. 张志超等^[3-5]采用虚拟激研究了车-桥耦合系统随机振动响应，得出了有益的结论. 彭勇波等^[6-7]基于概率守恒原理获取系统响应的概率密度信息，从而提出分析车-桥耦合系统随机振动的概率密度演化近似方法. 随机振动方法通常求解效率较高，但通过这种方法得到系统响应统计量，并根据均值和方差进行评价，且多数方法需要对轮轨进行线性化处理.

在实际工程中，常用实际测量或计算机模拟得到桥面或轨道不平顺，用时域法求解非线性车-桥耦合系统的动力响应^[8-11]. 采用直接MCS法进行可靠性评价时需要大量样本，当桥梁为大跨度桁梁桥时，计算耗时会进一步增大. 针对其中关键响应指标，采用动态系统代理模型可显著减少计算时间. 陈皓等^[12]提出了一种基于贝叶斯算法的NARX （nonlinear autoregressive with exogenous input）模型，能较精确地预测列车脱轨系数. 庞学苗^[13]以一种改进的NARX神经网络，提高轮轨相互作用力的预测精度. Han等^[14]用一种NARX-ANN （带外生输入的非线性自回归人工神经网络）代理模型分析了轨道高低不平顺作用下车-桥耦合系统随机振动的竖向响应. Xiang等^[15]用一种子集分裂模拟结合NARX代理模型的方法分析了小概率条件下车桥耦合竖向系统的随机振动. 通过代理模型的方法，可显著减少每个样本的计算时间，提高计算效率，但采用NARX-ANN训练样本时所需时间仍然较长，此外所有的时间滞后项数也较多，单个样本的计算效率有待进一步提高.

ARMAX （autoregression and moving average with exogenous）模型作为多项式线性回归模型^[16]，与传统的ARMA （autoregression and moving average） 相比，既融合了AR （autoregression） 模型和MA （moving average）模型的特点，线性组合现在的干扰值、过去的观测值和干扰值，又引入与输出序列相关的另一序列（即输入序列）作为回归项，提高模型的准确性.

本文采用ARMAX代理模型研究三维列车过桥的行车安全性. 从车辆模型、桥梁模型、轮轨相互作用力及轨道几何不平顺模拟等方面介绍车-桥耦合振动模型. 回顾ARMAX模型的原理，给出基于ARMAX代理模型的车-桥耦合系统随机振动响应分析的算法步骤. 讨论ARMAX代理模型预测列车横向和竖向加速度响应的效率和精度，分析车-桥耦合系统的可靠性.

1.   车-桥耦合振动模型

1.1   车辆模型

列车可视为四轴二系悬挂质量弹簧阻尼车辆^[17]，车体各部件由弹簧元件和阻尼元件连接，包含1个车体、2个转向架和4个轮对共7个刚体部件. 车体和转向架有横移、浮沉、侧滚、点头及摇头5个自由度，轮对有横移和摇头2个自由度，整节车辆共23个自由度. 基于达朗贝尔原理可建立车辆运动方程，见式（1）.

式中：${\ddot {\boldsymbol{y}}_{\text{v}}}$、${\dot {\boldsymbol{y}}_{\text{v}}}$和${\boldsymbol{ y}}_{\mathrm{v}} $分别为车辆的加速度、速度和位移向量；F_bv为列车受到来自桥梁轮轨的相互作用力矩阵；M_v、C_v和K_v分别为车辆质量、阻尼和刚度矩阵.

分析中采用德国ICE3高速列车参数，车速取200 km/h，列车为8车编组. 为讨论方便，8车均采用动车参数进行分析.

1.2   桥梁模型

桥梁采用跨度为32 m的预应力混凝土箱形截面桥^[14]，截面面积为9.031 m²，竖向惯性矩为10.939 m⁴，横向惯性矩为85.920 m⁴.

在商业软件中，采用空间梁单元建立桥梁单元，并通过模态分析得到桥梁结构的自振特性，同时导出桥梁模型的质量矩阵和刚度矩阵，结合瑞利阻尼计算方法，得到桥梁模型的阻尼矩阵. 将桥梁模型的质量、刚度和阻尼矩阵导入MATLAB中可得到桥梁运动方程，见式（2）.

式中：${\ddot {\boldsymbol{y}}_{\text{b}}}$、${\dot {\boldsymbol{y}}_{\text{b}}}$和$ {\boldsymbol{y}}_{\text{b}} $分别为桥梁加速度、速度和位移向量；F_vb为桥梁受到来自列车的轮轨相互作用力矩阵；M_b、C_b和K_b分别为桥梁质量、阻尼和刚度矩阵.

本研究中，为减弱非平稳特性的影响，桥梁设置为50跨，均为简支梁，共长1 600 m，墩高为10 m，墩梁均为混凝土结构，阻尼比为2%.

1.3   轮轨相互作用

轮轨接触力包括法向力和蠕滑力，法向力公式中的部分系数是蠕滑力的函数，蠕滑力公式中的蠕滑系数又是法向力的函数，所以法向力与蠕滑力互为隐函数，分析中轮轨接触法向力采用轮对的竖向和侧滚振动方程进行迭代求解.

轮轨接触法向力为

式中：N_L、N_R分别为左、右轮法向力（L表示左，R表示右，余同）；δ_L、δ_R分别为轮轨左、右接触角；θ_w为侧滚角位移；系数A、B为蠕滑力的函数，见式（4）.

式中：m_w为轮对质量；$ {\ddot y_{\rm{w}}} $为轮对横移加速度；$F_{y{\rm{pw}}}$为一系悬挂系统对轮对的横向力；F_zpw为一系悬挂系统对轮对的垂向力；y_w、ψ_w分别为横移位移和摇头角位移；T_yL、T_yR （T_xL、T_xR）分别为左、右轮的横向（纵向）蠕滑力.

蠕滑力采用Kalker线性理论进行计算，不计自旋蠕滑影响，见式（5）.

式中：T_x、T_y分别为纵向、横向蠕滑力；M_z为旋转蠕滑力矩；f₁₁、f₂₂、f₂₃为蠕滑系数，由法向力、轮轨接触几何参数和Kalker系数等进行计算，可参考文献[18]；ξ_x、ξ_y分别为纵向、横向蠕滑率，见式（6）.

式中：${\lambda }$为踏面斜度；r₀为轮对中心滚动圆半径；y_a为轨道中心线横移位移；θ_c为轨道水平不平顺位移；b为轮对两滚动圆间距之半；$ {\dot \psi _{\rm{w}}} $为轮对摇头角速度；v为轮对中心速度；$ {\dot y_{\rm{w}}} $为轮对横移速度；$ \dot \theta _1^* $为轮对合成转角速度.

轮轨接触力的计算参数会随轨道不平顺变化而变化，具有一定非线性轮轨接触，而代理模型对这类系统参数变化问题适用性相对较差，为减弱系统参数的非线性和时变特性的影响，将式（6）的部分参数进行线性化处理，即：滚动圆半径r₀近似取为0.43；踏面斜度${\lambda } $取为0.04；b取为0.75. 而该简化对列车响应的极值影响相对较小.

由于未考虑轮对的竖向自由度，并假定轮对和轨道为刚体，且与轨道密贴，因此，竖向的轨道不平顺作用会直接传递到转向架，可写成式（7）形式^[19].

转向架浮沉振动荷载：

转向架侧滚振动荷载：

转向架摇头振动荷载：

式中：m = ( j−1) × 2 + l，控制轮对求和，当j = 1时表示前转向架，i = 1，2，当j = 2时表示后转向架，i = 3，4；z_wi和${\dot {\textit{z}}_{{\rm{w}}i}}$为轮对中心竖向移位和速度；${\theta _{{\rm{w}}i}}$和${\dot \theta _{{\rm{w}}i}}$分别为轮对中心侧滚移位和速度；K_1y、K_1z、K_1θ分别为轮对横向、竖向、侧滚刚度；C_1y、C_1z、C_1θ分别为轮对轴箱悬挂横向、竖向、侧滚阻尼；l₁为轮对轴距之半；h₄为转向架构架中心到车轴中心线的高度.

在得到轮对的法向力和蠕滑力后，根据轮轨接触几何关系，可得到车辆和桥梁的相互作用力.

1.4   轨道几何不平顺的模拟

采用AR法模拟轨道几何不平顺，该方法是一种线性自回归模型，也是全极点模型^[20]. AR模型作为一种线性滤波器，能近似将白噪声信号转化为具有指定谱特征的离散随机过程. AR模型可用差分方程表示，见式（8）.

式中：W(t)为自回归序列；p为模型阶数，取值为4；A(k)为自回归系数；k为1 ~ p中的任意一阶模型阶数；t为时间；Δt为时间步长；U(t)为均值取0、方差取σ²的白噪声序列.

基于维纳-辛钦定理可得到轨道不平顺功率密度谱的自相关函数，见式（9）.

式中：S_w( f )为轨道不平顺功率密度谱函数；f为频率；R_w( dΔt )为S_w( f )的自相关函数；d = 0，1，$\cdots，p $.

再根据尤拉-沃克方程建立自相关函数与AR模型自回归系数A(k)的关系式，见式（10）

联立式（9）和式（10a）求得自回归系数A(k)，并代入式（10b）中求得R₀；根据式（11）求得白噪声序列U(t).

式中：n(t)为满足标准正态分布的白噪声.

生成轨道几何不平顺时，取空间步长为0.2 m，S_w( f )采用德国低干扰谱.

针对式（6）中被简化的参数，采用与文献[21-22]相同的车辆参数，对比在正弦波轨道不平顺^[21]和采用德国低干扰谱作为随机轨道不平顺输入时（按式（8）进行模拟），简化前、后车辆的加速度响应如表1所示. 由于文献[21-22]中车辆模型与本文模型有一定差异（文献中车辆为35个自由度，并考虑轨道的影响），但考虑在激励相同时，计算结果基本吻合（见表1）；虽然在正弦激励下，简化后的车体横向加速度明显变大，但考虑在随机激励作用时，模型参数简化对车体加速度响应影响有限，后续研究均采用随机激励.

表  1  加速度响应结果对比

Table  1.  Comparison of acceleration responses m/s²

项目		正弦激励			随机激励
		竖向	横向		竖向	横向
文献[21]		0.775	0.490
文献[22]		0.735	0.480
本文结果	简化前	0.818[23]	0.535[23]		0.369	0.271
	简化后	0.836	0.687		0.369	0.281

下载: 导出CSV 

| 显示表格

2.   动态系统的ARMAX代理模型

代理模型类似于一个黑箱，在不了解复杂系统内部详细结构的情况下，通过输入和输出数据来训练代理模型，再次输入数据便可预测动态响应. 该模型是针对某一个指标，通过在时间轴上迭代来计算响应，计算效率会明显提高.

2.1   ARMAX模型

ARMAX模型是带外源输入的自回归滑动平均模型，也是一种线性输入-输出多项式模型，可当做一个黑箱模型. 单输入和单输出的ARMAX 模型结构^[24]如图1所示. 图中：p(t−n_k)为输入信号，n_k为输入的延迟阶数；v(t)为输出信号；e(t)为均值为0，方差不为0的白噪声信号；Z(q)、H(q)与G(q)分别为模型输出、输入和噪声项，见式（12）.

图  1  ARMAX模型结构

Figure  1.  ARMAX model structure

下载: 全尺寸图片 幻灯片

式中：n_z为自回归项阶数；n_h为回归项阶数；n_g为噪声项阶数；[z₁ z₂$\cdots\; {{\textit{z}}_{{n_{\text{z}}}}}$]、[h₁ h₂$\cdots\; {h_{{n_{\text{h}}}}}$]和[g₁ g₂$\cdots\;{g_{{n_{\text{g}}}}}$]为模型参数；q为后向移位算子，可由式（13）求得.

式中：n取值为n_z、n_h + 1、n_g.

因此，ARMAX模型可表达为

式中：Z(q)v(t)为自回归项，即AR部分；G(q)e(t)为滑动平均项，即MA部分；H(q)p(t−n_k)为外源输入（也可称为回归项），即X部分.

联立式（12） ~ （14），即可得到ARMAX模型的差分方程为

分析中，延迟阶数通过试算法得到. 在试算过程中，可观察系数A、B、C随延迟系数的变化情况，若延迟系数非常小，则可删除相应的项.

2.2   分析过程

对车-桥耦合系统模型而言，列车所受荷载除轨道几何不平顺外，还有桥梁位移. 但在分析中，桥梁位移与轨道不平顺有一定关系，考虑计算的效率，针对轨道几何不平顺为0时列车过桥时各个轮对处的桥梁响应为输入，与轨道几何不平顺叠加，形成等效的轨道不平顺. 这个过程相当于仅进行了一次全过程迭代，但桥梁的刚度通常较大，一次迭代的桥梁响应已与实际曲线较为吻合. 具体分析时，采用以下步骤：

步骤1　轨道几何不平顺为0时，提取列车过桥时每个轮对对应处的振动位移和加速度.

步骤2　任意生成一条随机的轨道不平顺，进行车-桥耦合振动计算，得到车辆的响应，将轨道几何不平顺与步骤1中的桥梁响应叠加，得到输入和输出数据.

步骤3　试算代理模型延迟阶数，预测响应，基于均方根值和最大绝对误差，判断预测准确性，直到确定代理模型参数为止.

步骤4　随机生成轨道几何不平顺，与步骤1中的桥梁响应叠加，预测响应，直到所需数量. 寻找响应绝对最大值所对应的样本，以相应的轨道不平顺为输入，进行车-桥耦合振动计算，验证代理模型预测的准确性.

步骤5　若代理模型预测准确性差，则以最大值所对应的输入和输出重新训练代理模型，重新进行步骤4的计算，并将两次预测结果进行平均值处理.

步骤6　统计所有预测样本的信息，形成超越概率曲线，评价行车安全.

不同的车辆响应输入有所差异. 对于竖向响应，主要由竖向的不平顺决定，如车体竖向加速度，主要和式（7）有关. 因此，采用代理模型计算竖向加速度时，在生成轨道几何不平顺后，按式（7a）计算竖向的荷载，得到相应的输入用于训练代理模型和预测响应. 对于横向加速度，主要受横向和竖向轨道几何不平顺的影响，但轮轨横向力（见式（3））不能直接计算得到，因此可直接将4个轮对处的横向轨道不平顺和式（7b）作为输入.

对于步骤5的处理方式，文献[14]直接采用最大值样本的处理方式，但本文中的ARMAX代理默认的基函数为线性函数，采用最大值时仍有一定误差，通过对比发现取平均值的精度会更好. 当然，若任意产生的样本即为所模拟样本集中的最大值，则需要再生成一个响应最大值较小的组样本重新训练代理模型.

分析中，取一条样本中绝对最大值进行统计分析，假如有1000个样本，最大值样本是绝对最大值所对应的时程曲线.

3.   结果分析

为讨论ARMAX代理模型的准确性和计算效率，针对车-桥耦合振动模型，采用直接蒙特卡罗法计算了1 000个样本，每个样本的计算时间约5 min. 由计算结果的统计特征观察到，在车辆参数相同，8辆车编组时，各车辆的车体竖向加速度a_v和横向加速度a_h绝对最大值的均值、1%和99%分位处的响应差异均较小（见图2），因此，后续主要讨论第1辆车的响应.

图  2  8辆车的车体加速度响应

Figure  2.  Acceleration response of 8 car bodies of the train

下载: 全尺寸图片 幻灯片

图2的加速度响应值较小，这在类似分析算例中也有相似的结果^[25-26]. 这是因为分析中计算车速远小于其最高运营速度，加之桥梁刚度较大，因此车辆响应较小. 由于轨道不平顺对桥梁响应影响较小，因此后续不再讨论桥梁响应.

针对第一辆车的车体竖向加速度进行代理模型训练，建立代理模型过程中，通过试算得n_z = 4，n_h = 5，n_g = 0，n_k = 0，此时训练的代理模型准确性高. 相比NARX等其他代理模型，此时ARMAX代理模型的延迟阶数明显较少，因此，计算效率会明显提高，结果显示每个样本的计算时间不超过0.2 s. 竖向系统不需进行2.2节中的步骤5便可得到预测结果. 图3给出了第1辆车的车体竖向加速度响应的理论值和代理模型的预测值对比. 图中：横轴为目标值，由式（1）和式（2）进行直接MCS得到；纵轴为预测值，由代理模型预测得到. 由图3可见：ARMAX模型对车-桥耦合竖向系统车辆加速响应的预测具有良好精度和效率.

图  3  竖向加速度的目标值和预测值

Figure  3.  Target value vs predicted value ofvertical acceleration

下载: 全尺寸图片 幻灯片

列车横向系统为非线性系统，而ARMAX代理模型是线性模型，在模拟横向系统时会有一定误差. 图4给出了横向加速度响应的最大绝对值. 建立代理模型过程中，通过试算得到n_a = 50，n_b = 50，n_c = 10，n_k = 0. 由图4可知：初始样本训练得到的代理模型预测值总体较目标值偏大，而最大值样本所得到的代理模型预测值较总体目标值偏小，取两者平均值后准确性明显提高.

图  4  横向加速度的目标值和预测值

Figure  4.  Target value vs predicted value oflateral acceleration

下载: 全尺寸图片 幻灯片

为进一步讨论代理模型预测车体横向加速度响应的准确性，图5（a）给出了车体横向加速度目标时程和预测时程的对比曲线，图5（b）给出了目标值和预测值的最大误差时程曲线. 由图5（b）可见：最大误差的幅值小于0.03 m/s²，所以利用ARMAX代理模型预测非线性的横向系统响应时仍有一定准确性. 同样，车体横向加速度每条样本的计算时间小于0.2 s，代理模型的精度和效率是可接受的.

图  5  车体横向加速度时程曲线

Figure  5.  Lateral acceleration time-history curves ofcar body

下载: 全尺寸图片 幻灯片

铁路线路有固定的养护周期，每天运行列车的对数是确定值，因此，评价行车安全的概率水平也应确定，这表明从概率角度考虑行车安全可行. 为进一步研究行车安全性，图6（a）、（b）给出车体竖向和横向加速度的超越概率曲线，当超越概率为10⁻³时，列车加速度响应仍较小，这在文献[18]中也有类似的结果，同时加速度响应较小也与计算速度较低有一定关系，这表明行车安全具有较高的安全储备. 此外，在代理模型基础上，按MCS统计超越概率时，在保证变差系数为0.1的条件下，模拟的次数为超越概率倒数100倍，由图6可见：在超越概率为10⁻²时，预测值与目标值已较为吻合. 这表明基于ARMAX代理可近似车桥系统响应.

图  6  车体加速度超越概率

Figure  6.  Exceedance probability of car body acceleration

下载: 全尺寸图片 幻灯片

4.   结　论

本文采用MATLAB编写车-桥耦合系统的分析程序，提出一种ARMAX代理模型近似得到车-桥耦合系统中的车辆竖向和横向响应，并与MCS计算结果进行对比，探究车-桥耦合系统可靠性分析的效率和精度，得出如下结论：

1） ARMAX代理模型预测车辆竖向和横向加速度响应的时间短，对车-桥耦合系统进行可靠性分析的效率较高；

2） 就车-桥耦合系统车辆加速度响应而言，竖向响应的预测结果与MCS结果吻合较好，横向响应的预测精度虽然略低于竖向，但平均处理后的超越概率曲线与直接MCS得到的结果吻合较好.