用于铁路线路基础结构伤损识别的概率密度分布演化方法

徐 磊，翟婉明，陈宪麦，陶伟峰

(1.西南交通大学牵引动力国家重点实验室，四川成都 610031；2.中南大学土木工程学院，湖南长沙410075；3. 同济大学土木工程学院，上海 200092)

我国拥有庞大的铁路线路网，线路长期稳定、安全运营所需的状态监测管理工作重要性日益凸显。国内外一般采用综合检测列车对影响行车的线路异常情况进行实时检测与综合评判。影响车辆-轨道耦合系统(以下简称此系统)的随机因素众多，各因素之间的相互作用机理尚未完全明确。此系统的激扰源具有随机性和演变性(如常见的轨道不平顺)，在理想的平顺状态下能够反映到车辆响应中的线路伤损，极易被线路不平顺产生的响应掩盖。从随机不平顺激扰下的车辆响应中解耦出线路损伤产生的响应亦是困难的，因为除不平顺的随机变化特性外，线路本身的动力参数也是随机的。因此，通过车辆系统的振动响应实时反馈线路基础结构的状态及伤损情况是一个艰巨的研究课题。

现阶段，一般基于计算机或雷达视觉分析轨道系统图形特征，进行病害的诊断[1-2]。本文基于随机分析的一般理论，考虑线路不平顺的空间随机性和长时效演变性，以及轨道结构动力参数的随机性，将车辆-轨道耦合动力学理论[3]与概率密度演化方法[4]结合，提出一种用于铁路线路基础结构伤损识别的概率密度分布演化方法。该方法以车辆-轨道耦合动力模型为核心，在妥善解决用于动力模型激励的轨道不平顺随机性及演变性模拟、线路基础结构动力参数随机模拟与组合等两个问题的基础上，将正常与伤损下的动力响应进行概率密度分布的差分，从而发现线路损伤。

1 轨道不平顺的随机性和演变性模拟

对于某段待监测的铁路线路，采用如下方法模拟此区段轨道不平顺的随机性和长时效演变性。

(1)制定不同线路条件下的累计概率谱。文献[5]指出轨道不平顺谱近似服从自由度为 2 的χ2分布函数。可针对不同的铁路线路区段(路基、桥梁、隧道等)，构建相应的累计概率谱(图1)。

(2)通过长期的线路监测及统计分析，确定不同监测时段内线路区段轨道不平顺谱线的随机波动范围。若以月为监测的时间单位，假定第i月此线路区段轨道随机不平顺谱在累计上、下限概率分别为Ci,u和Ci,d的累计概率谱内波动；同时，用τi→k,u、τi→k,d表示从第i月到第k月的上、下限谱累计概率变化率，一般存在等式

Ck,u|d=Ci,u|d(1±τi→k,u|d)

( 1 )

(3)根据(2)中确定的线路不平顺谱线累计概率范围，采用拉丁超立方抽样(LHS)方法[6]进行此累计概率范围内轨道谱的抽样，对采集的轨道谱进行逆傅里叶变换[3]，生成随机轨道不平顺时域样本。

需要指出的是：

①每条谱样本由服从一定空间频率-谱密度规律的样本点组成。一般而言，频率越低，谱密度值越大，因此采用LHS法对不同频率下谱密度值抽样时，需满足临近低频点的谱密度更大、谱密度值对应的累计概率均在拟定的谱累计概率范围内这两个基本要求。

图1 武广高速铁路普通路基段轨道不平顺的累计概率谱(截止波长：1～120 m)

②计算及实测资料表明[7]，同一路段的轨道不平顺呈现“记忆”重复特性，即其随机波形具有一定的相似性。因此，在用逆傅里叶变换法[3]生成轨道随机不平顺时域样本时，取不同不平顺时域样本对应的独立相位序列ζn相同，即ζn=exp(iΦn)，Φn服从0～2π均匀分布。

2 参数随机的轨道模型构建与计算

国内外基于车辆、轨道系统参数确定的动力分析居多[3]，这在线路的初始设计、评估阶段十分有用。随着车-线系统的长期运营，此系统参数和激励源一般会发生随机演变，从而进入随机振动的范畴。本文假定车辆参数确定，主要简述线路参数随机时的模拟与组合方法和线路参数随机时轨道系统的构建及计算方法。

2.1 线路参数随机模拟与组合

以板式轨道为例，阐述动力参数(取轨下胶垫垂横向刚度、阻尼参数，CA砂浆垂横向刚度、阻尼参数)的随机模拟与组合方法。在时间域，假定不同的随机参数服从正态分布[8]且独立；在空间域(即沿线路区段)，同一参数的随机模拟样本间具有相关性，而不同随机参数的模拟序列是独立的。

上文已述及同一结构层的系统部件(如轨下胶垫垂向刚度)参数样本间具有相关性，其相关系数为

( 2 )

通过式( 2 )，可以形成n维相关正态分布随机变量抽样序列的相关系数矩阵

( 3 )

步骤3通过线性变换，令

( 4 )

式中：μ为均值矩阵。

根据Cantor集合映射的思想[12]，依据不同样本的出现概率进行不同随机参数的样本序列组合，实现动力计算样本的降维[13]。

2.2 轨道模型建模与计算

线路参数随机与否对轨道系统的建模与计算影响较大。将2.1节中模拟出的随机参数组合序列均布在轨道模型中，采用有限单元法进行轨道系统建模。

2.2.1 轨道模型

轨道模型[14](图2)采用板式轨道结构，取相邻轨下胶垫间距为一轨段单元，钢轨为连续点支承Euler梁，钢轨垫层为线性刚弹簧及阻尼器单元，轨道板通过CA砂浆与路基连接，CA砂浆模拟为线性均布面弹簧和黏滞阻尼器，每个轨段单元有40个自由度，可表示为

( 5 )

式中：

δ1=[U1LrV1LrW1LrθX1LrθY1LrθZ1LrU1RrV1RrW1RrθX1RrθY1RrθZ1RrV1LsW1LsθX1LsθY1LsθZ1LsW1RsθX1RsθY1Rs]T

( 6 )

δ2=[U2LrV2LrW2LrθX2LrθY2LrθZ2LrU2RrV2RrW2RrθX2RrθY2RrθZ2RrV2LsW2LsθX2LsθY2LsθZ2LsW2RsθX2RsθY2Rs]T

( 7 )

式中：U、V、W分别为沿纵向X、横向Y、垂向Z的位移；下标1、2分别为沿Y轴方向的右侧自由度和左侧自由度；下标L、R分别为沿Z轴方向的左侧自由度和右侧自由度；θ为转角位移；下标 r、s分别为钢轨单元和轨道板单元。

图2 板式轨道模型三维视图

2.2.2 循环计算方法

由于线路参数在时间-空间域具有随机特性，在车辆-轨道耦合动力计算时，定点激励法[3]不能反映线路参数随机时车轨系统振动响应的时变特征，须采用移动车辆激励法考虑车辆与轨道之间的相对运动状态。采用移动车辆激励法时，往往存在计算长度上的限制，特别是采用有限元建模时，其自由度较多，将占用较大的计算机内存并降低计算效率。采用循环计算方法，可有效解决这一类问题。

循环计算方法的实现过程为：

步骤1取M为前期平稳运行长度，N为循环计算长度，L为数值积分时车辆前后截取的轨道模型长度，同时设置全局坐标系xOy和局部坐标系x′O′y′(图3)。在平稳运行阶段，亦可只取长度为M的轨段单元进行数值积分，以减小时域积分算法中等效刚度矩阵的规模，减少求解时间。

图3 循环计算模型

步骤2若已知轨下胶垫间距，则可以计算长度N、L所占的轨段单元数目nN、nL。若第1、4号轮对在全局坐标系中轨段单元编号为n1、n4，车辆在循环计算区段的运行距离为N′，那么前、后轮对所接触的局部轨段单元编号分别为

nf=n1+nL-([N′/N]-1)nN

( 8 )

nb=n1-nL-([N′/N]-1)nN

( 9 )

式中：[·]表示朝0方向取整。

式( 8 )、式( 9 )是实现循环计算的关键，其能保证车辆始终在循环区段N内不断运行。而前、后轮对所接触的全局轨段单元编号分别为

(10)

(11)

根据式( 8 )～式(11)计算出的前、后轮对所在的局部和全局轨段单元，以及L、轨下胶垫间距等已知参数，可以得到每一积分步内整车所占用的轨段单元序列编号。

步骤3前期已经将2.1节中模拟出的随机参数序列按轨段单元序列的形式编号，即全局轨道单元序列与随机参数序列一一对应，因此可以较方便地依据步骤2中计算出的整车所占用的轨段单元序列编号，择取相应的线路随机参数序列参与车轨耦合动力计算。

如图3所示，每循环一次，需重新生成一个长度为N+2L+2(Lc+Lt)(Lc为前后构架中心间距的一半，Lt为轴间距离的一半)的轨道单元矩阵(包括刚度矩阵、阻尼矩阵和质量矩阵)。

3 车辆模型

将车辆视为多刚体系统，考虑车体、转向架及轮对的横移、沉浮、侧滚、点头和摇头运动，即车辆系统共35个自由度。具体建模方法参考文献[3]。

4 概率密度演化方法

文献[4]在结构随机分析方面发展了一类概率密度演化方法，用于解决系统随机激扰-振动响应间的概率密度传递问题。

车辆-轨道耦合系统随机动力方程

(12)

式(12)的解答存在且唯一，并与随机向量有关。因此，可写成如下形式

X(t)=ϑ(Θ,t)

(13)

式中：ϑ为集合符号。

根据概率守恒原理，[X(t)ΘT]T的概密度函数pXΘ(x,θ,t)满足如下广义概率密度演化方程[4]。

(14)

式(14)的初始条件为

pXΘ(x,θ,t)|t=0=δ(x-x0)pΘ(θ)

(15)

联合求解式(14)和式(15)，可得pXΘ(x,θ,t)，随后积分可知X(t)的概率密度函数为

(16)

5 用于线路状态监测的概率密度演化模型

车辆-轨道耦合系统中的轨道不平顺及系统参数均具有时-空随机演变特性(图4)，导致基于车辆系统振动响应的“反馈式”线路状态监测工作困难重重。一般来说，由于线路基础结构伤损过程的渐变性和此系统本质上的随机演变性，系统参数、激励源确定的动力计算理论对于线路随机伤损(尤其是小伤损)的实时监测与预判工作略显乏力。

图4 轨道不平顺和线路基础结构的时-空演变示意

现阶段，国内外大部分线路伤损监测的研究均基于特殊或冲击信号，取一段确定性信息进行特征状态的提取一般可获得满意结果，但往往忽视了线路系统(包括轨道不平顺)在列车循环动荷载及其他因素干扰下的随机演化特性。这些不断发展的车轨激励形态给线路状态监测带来了困难。例如：轨道随机不平顺的存在极有可能掩盖那些较微弱的线路伤损造成的动力响应差异，从而造成识别的困难。

针对线路基础结构的伤损预估和动力性能评判问题，作者认为将随机分析方法引入确定性计算理论[3]，统计此系统在不同的参数或激励源随机域下的动力响应特征，当线路状态恶化时，根据此系统确定随机域下的动力响应特征变化预估线路状态恶化与否不失为一种可行的办法。

本文用概率密度分布表达这种动力响应特征，动力响应的概率密度分布基于概率密度演化方法将众多具有概率性质的确定性计算结果融合，是在确定的随机域激扰下此系统响应的全面反映。线路状态恶化时，动力响应的概率密度分布会发生相应的变化，通过此变化即可识别线路伤损，这就是用于铁路线路状态监测的概率密度演化方法，其基本模型或分析流程如图5所示。

图5 用于铁路线路状态监测的概率密度分布演化模型

6 算例分析

本文主要采用动力仿真手段，验证所提方法的有效性。行车速度为350 km/h，轨道不平顺截止波长范围1～120 m。

将线路区间的监测时间长度取为1年，将“月”作为监测时间单位，将有12次线路状态的监测评判。基本计算参数如下：

(1)将线路区段随机不平顺激扰谱的初始累计概率范围取为0.5～0.6，终止累计概率范围取为0.1～0.9，假定此线路区间不同时段间的激扰谱累计概率服从均匀分布。

(2)由于本文的研究重点为线路基础部件的伤损识别，故取较小的线路系统参数随机域，即：各参数正态分布中的标准差较小，同一基础结构层的线路参数相关系数区间为0.7～0.9。

(3)动力计算里程为1 km。本文分别在线路里程K0+400、K0+425、K0+615三处的两侧承台设置了轨下胶垫横向动力系数(即横向刚度、阻尼系数)的伤损；在里程K0+137、K0+330、K0+520三处的两侧承台设置了轨下胶垫垂向动力系数(即垂向刚度、阻尼系数)的伤损。基础部件伤损表现为动力参数值的减小，本文假定损伤部件在不同监测时段间的动力参数值λi可以通过参数μi∈[0,1](i=1,2,…,12)进行调节，若初始参数值为λ0，则λi=λ0(1-μi)。

限于篇幅，图6、图7仅给出i=2，8(μi=0.09，0.73)时通过车辆系统轮对横向加速度进行轨下胶垫伤损状态监测时的概率密度分布差分图，主要步骤参考图5。

图6 μi=0.09时轨下胶垫伤损识别结果

图7 μi=0.73时轨下胶垫伤损识别结果

从图6、图7可知，通过概率密度分布差分图能较好地检测出轨下胶垫在不同比率下的伤损情况，且定位准确；伤损程度不同，对应的概率密度分布差分值亦不同，一般伤损越严重，差分值越大。如图6所示，即使是轨下胶垫的小伤损状态(动力系数仅减小了0.09)，仍引起轮轨横向加速度概率密度分布较大的变化，继而可以通过此变化评判线路基础结构的伤损状态。

此外，图8给出了不同阶段(1～12月)的概率密度差分上、下极限值。

图8 不同阶段的概率密度差分极限值

由图8可知，随着轨下胶垫的动力参数耗损及随机不平顺的演化，伤损位置处动力响应的概率密度差分值范围逐渐变大，基本与伤损程度成正比，这为确定不同损伤状态下的概率密度差分值标准奠定了基础。

7 结论

激励源和系统参数确定时，利用车辆-轨道耦合动力学方法能较好计算出线路基础部件伤损时的系统响应形态。由于车辆-轨道系统的伤损位置、时间具有未知性和渐变性，为了防患于未然，需对线路基础部件的伤损状态进行跟踪测算与评估。

本文根据随机分析的一般理论，将此系统的不平顺激励源和线路参数均做随机处理，考虑此系统的时-空随机演变特性，提出了从车辆系统动力响应的概率密度分布差分结果中发现、定位线路伤损程度及位置的基本设想及方法。计算结果表明：

(1)采用本文给出的概率密度演化差分方法，能够准确定位轨下胶垫的垂横向损伤，表明了方法的可行性；同时，概率密度差分值可用于表征线路基础结构伤损程度。

(2)以往的研究工作多取确定性的实测或仿真信号进行伤损识别，忽视了线路在循环随机荷载等因素作用下的时变本质。本文重点考虑轨道不平顺和线路参数的随机演化特性，提出更通用的一般性方法，可供借鉴。