基于导波多模态融合的无缝钢轨温度应力估计算法

王 嵘，余祖俊，朱力强，许西宁

(1.北京交通大学 机械与电子控制工程学院，北京 100044；2.北京交通大学 载运工具先进制造与测控技术教育部重点实验室，北京 100044)

我国高速铁路规模的不断增大对高速铁路运行安全提出了更高要求。高速铁路通常采用的无缝钢轨线路在长度方向上不能自由伸缩，当钢轨温度发生变化时，其内部将产生较大的温度应力，严重时会导致胀轨、跑道及断轨等事故的发生，直接威胁行车安全。因此，在线监测无缝钢轨内部纵向温度应力并在超限前预警，对保障高速铁路的安全运营具有重要意义。

超声导波对应力敏感，在钢轨中传播时可引起钢轨横截面内全部质点振动，受钢轨表面残余应力的影响较小，已经成为检测钢轨纵向温度应力的一种重要方法。文献[1]通过对杆和平板的研究发现，随着应力增大，导波的相速度增大、群速度减小，通过测量群速度或相速度可对应力进行检测。文献[2]推导一维弹性波导介质在轴向应力作用下的波动方程，并分析应力作用下导波的传播特性。文献[3]获得钢绞线低频率段模态的群速度与拉应力标定曲线，通过测定钢绞线中该模态的群速度就能计算出钢绞线的应力值。文献[4]采用半解析有限元方法分析轴向拉应力作用下钢轨水平弯曲模态、垂直弯曲模态和扭转模态在1 kHz以内相速度的变化趋势。文献[5]选取一种对应力敏感的35 kHz导波模态，获得该模态在不同应力下的群速度标定曲线，通过测定钢轨中该导波的群速度即可计算钢轨温度应力。

以上基于超声导波的应力检测技术研究都假设波导体的弹性模量是常数，实际服役中的无缝钢轨，当环境温度变化时，其弹性模量将发生明显变化。文献[6]推导在温度发生变化时，锁定钢轨弹性模量和温度应力变化对相速度的影响，发现前者的影响比后者大一个数量级。为了准确测量温度应力，需首先获得钢轨的应力-相速度-弹性模量三维标定曲面。这种方法在工程中难以实现，因为获取该曲面的工作量大且较繁琐，在使用时还需要同时测量钢轨弹性模量和相速度。因此，如何消除弹性模量的影响是目前急需解决的技术难题。本文利用钢轨多个导波模态的相速度(或者波数)对应力和弹性模量敏感程度不同的特性，提出一种多模态融合估计算法，克服弹性模量变化对应力检测的影响，在弹性模量未知的条件下估计钢轨温度应力。

1 理论基础

由声弹性原理可知，应力σ与波数ξ以及弹性模量E之间存在函数关系。已知波数ξ和弹性模量E时，利用此函数估计应力σ，记为

σ=f(ξ，E)

( 1 )

此时，应力的误差传递函数为

( 2 )

为了评估此方法的可行性，通过半解析有限元方法对有初始应力的钢轨进行建模，求解函数f及Δσ。

1.1 有初始应力的钢轨半解析有限元建模

对于规则横截面的波导介质，可以建立解析解的频散方程，研究其导波的传播特性。钢轨横截面较复杂，无法建立解析解的频散方程。半解析有限元方法是国际上普遍采用的求解任意横截面波动方程的一种有效方法，其基本原理是在波导介质的横截面上进行有限元离散，沿波导介质传播方向的位移用简谐波的振动方式表示，根据哈密尔顿原理可以推导出导波在其中传播的波动方程[7]。

我国高速铁路广泛使用60 kg/m钢轨，其截面如图1所示，定义其横截面为y-z平面，波传播方向为x方向。

图1 采用三角形单元离散的60 kg/m钢轨y-z横截面

( 3 )

( 4 )

假定导波沿x方向传播时波导介质做简谐振动，其中任意一点的位移可用空间分布函数表示为

( 5 )

式中：ω为导波角频率；t为时间参数。

ε=[εxxεyyεzzεxyεyzεzx]T=

( 6 )

采用半解析有限元法分析钢轨中导波的传播特性，可得到60 kg/m钢轨中超声导波的波动方程[7]为

( 7 )

当波导介质存在较大的初始应力σ(0)时，必须考虑初始应变与位移的二次项，即非线性应变-位移关系。有初始应力时，应变S与位移u的关系为

( 8 )

单位体积的应变能φ可表示为

( 9 )

(10)

式中：Sxx为x方向的非线性应变。式(10)中省略部分为高阶无穷小量。

(11)

(12)

初始应力作用下波导介质总体积应变能为

(13)

因此，含有轴向初始应力的钢轨波动方程[9]为

(14)

1.2 波数与弹性模量测量精度对应力测量精度的影响

将式(14)移项，并令

G(ξ，σ，E)=[ξ2·(K2+K(0))+ξ·K1+

(15)

式中：σ为纵向温度应力；G(ξ，σ，E)为关于ξ、σ、E的函数。

根据隐函数求导公式，求解出应力对波数的偏导数为

(16)

根据隐函数求导公式，应力对弹性模量的偏导数为

(17)

σ=α1f1(ξ1，E1)+α2f2(ξ2，E2)+…+

(18)

多模态融合估计应力误差传递函数为

(19)

仿真试验数据表明，波数测量误差导致的应力误差为0.243 5 MPa。若使应力估计误差在2 MPa以内，根据式(19)可知，弹性模量的测量误差应在0.011 13 GPa以内，相对误差为0.005 3%。这个量级的测量精度在对弹性模量的实际测量中是难以达到的。因此，采用传统的补偿弹性模量方式不可取，本文利用不同模态对弹性模量的敏感度不同这一特点，抵消弹性模量对应力估计结果的影响。

2 基于多模态融合的应力估计方法

通过半解析有限元程序计算特定频率下各个模态对温度应力和弹性模量的敏感性，定量分析两者对相速度的影响。利用应力和弹性模量对相速度敏感程度不同这一特性，基于应力估计值最小原则选出特定的模态对应力进行估计。最后，通过数据验证方法的正确性。

2.1 多模态对应力与弹性模量敏感度分析

在特定频率下，各个模态对温度应力和弹性模量敏感程度的判断有两种方法。一是通过导波模态的振型判断。现有研究结果表明，水平弯曲模态和扭转弯曲模态对应力敏感，伸展模态对弹性模量敏感。因此只要利用有限元方法获得导波模态的振型，根据主要振动形态即可大致估计其对应力与弹性模量的敏感程度。二是通过对比应力和弹性模量变化前后的频散曲线判断。例如，对比应力变化前后，特定频率下各个模态相速度的变化率，前后状态变化率最大的是对应力最敏感的模态；模态对弹性模量变化的敏感程度同理。第二种方法可以定量分析各模态对温度应力和弹性模量的敏感性，因此本文采用第二种方法论述。

对于60 kg/m钢轨，温度变化1 ℃大约等效于弹性模量变化0.05 GPa，锁定状态下温度应力变化2.5 MPa[9]。假设室温为20 ℃时，弹性模量为210 GPa，应力为0 MPa。若温度的变化范围为-20～60 ℃，则应力变化范围为-100～100 MPa，弹性模量变化范围为208～212 GPa。以35 kHz为中心频率作为激励信号，计算出3种状态下各模态的相速度(当频率一定时，相速度与波数成反比，求解相速度等价于求解波数)：第一种状态为弹性模量210 GPa、应力0 MPa；第二种状态为弹性模量210 GPa、应力100 MPa；第三种状态为弹性模量212 GPa、应力0 MPa。3种状态下不同模态相速度对比见表1。从表1可以看出，1号模态对温度应力最敏感，20号模态最不敏感；19号和20号模态对弹性模量最敏感，1号模态最不敏感。说明当弹性模量一定时，随着相速度值增大，模态对温度应力敏感度减小；当温度应力一定时，随着相速度值增大，模态对弹性模态敏感度增大。从总体上看，不同模态对温度应力和弹性模量的敏感程度不同。第一种与第三种状态相速度变化率明显大于第一种与第二种状态相速度变化率，即在温度变化工况下弹性模量比温度应力对波数的影响更大。

表1 35 kHz导波在不同弹性模量状态和不同温度应力作用下在钢轨中传播的相速度对比

2.2 多模态估计方法

为了消除弹性模量测量对应力测量的影响，提高应力估计精度，采用应力估计最小误差的方法选取多模态组合对应力进行估计。假设波数ξ与温度应力σ和弹性模量E之间存在线性函数关系g，波数ξ表达式为

ξ=g(σ，E)=aσ+bE+η

(20)

式中：a、b、η为系数，其表达式分别为

(21)

(22)

当温度应力σ和弹性模量E一定时，相应的参数a、b、ξ均可通过式(21)、式(22)计算出来，则参数η为

η=ξ-aσ-bE

(23)

当选取了k个模态，式(20)可改写为

ξi=gi(σ，E)=aiσ+biE+ηi

i=1，2，…，k

(24)

将式(24)改写成矩阵形式为

(25)

式中：ξ=[ξ1ξ2…ξk]T，A=[a1a2…ak]T，B=[b1b2…bk]T，η=[η1η2…ηk]T。

根据最小二乘法，温度应力σ以及弹性模量E的估计量表达式为

(26)

表2 不同弹性模量状态下不同模态组合对零应力的估计值

3 仿真试验验证

基于以上理论，实际测量温度应力分为两步。第一步，获得波数-应力标定函数。利用半解析有限元法建立60 kg/m钢轨模型。在钢轨模型的特定节点施加激励信号，如正信号或汉宁窗调制的正弦信号，通过解析求解得到远处某一点的振动响应。通过改变模型中温度应力和弹性模量参数获取不同温度应力与不同温度状态下导波振动响应，采用二次加权算法计算各个模态波数。通过BP神经网络建立波数与应力的标定函数。限定温度参考的范围可以提高应力估计精度，因此引入轨温作为特征参数。BP神经网络的输入信号是各模态的波数和轨温，输出信号是温度应力。第二步，在测量应力时，采集阵列数据和轨温，计算各模态波数，代入标定函数，估计温度应力。为了验证本文理论的可行性，采用仿真试验进行验证。

3.1 钢轨建模

导波在波导体中的传播特性一般用频散曲线c-f描述，其中f为导波信号的频率，c为导波相速度Cp。对于钢轨这类具有复杂不规则截面的波导体，频散曲线无法通过解析式描述，只能借助仿真方法(如半解析有限元法[7-8，10])获得数值解。

图2为采用半解析有限元法获得的60 kg/m钢轨在零应力状态下的相速度频散曲线，其中，钢轨弹性模量为E=210 GPa，泊松比为ν=0.3，密度ρ=7 800 kg/m3。从图2可以看出，在低频段，钢轨中主要存在4个导波模态；随着频率的增加，能够在钢轨中传播的导波模态也越来越多。图2表明，钢轨中的导波具有频散特性，即对于同一导波模态，频率变化时，导波的速度也发生变化，这将造成波形在传播过程中发生畸变。

图2 零应力状态下60 kg/m钢轨相速度频散曲线

当钢轨中存在拉应力或压应力时，导波的传播速度将发生改变。例如，钢轨中拉应力为100 MPa时，在如图3所示的低频段范围内(0～1 kHz)，各模态的相速度发生了微小变化，且各模态的变化程度(对温度应力的敏感度)不同。相速度值增加时，其对温度应力的敏感性逐渐减弱。

图3 弹性模量为210 GPa时不同应力状态下的相速度频散曲线

当钢轨弹性模量变化时，导波的传播速度将发生改变。例如，在零应力状态下，钢轨的弹性模量由210 GPa变为212 GPa时，在如图4所示的低频段范围内(0～1 kHz)，各模态的相速度发生了微小变化，且各模态的变化程度(对弹性模量的敏感度)不相同。相速度增加时，其对弹性模量的敏感性逐渐增强。

图4 零应力时不同弹性模量状态下的相速度频散曲线

3.2 模态激励与波数测量方法

在35 kHz频率下，首先挑选对温度应力敏感的模态和对弹性模量敏感的模态，同时应便于激励与接收、信号模态相对单一。通过模态振型分析，确定每一个模态的最佳激励与接收位置、方式，获取不同温度应力与不同弹性模量状态下的导波信号。最后，通过二次加权算法计算各个模态波数。在35 kHz频率下选择了模态1、模态3、模态7、模态20。设置11种温度应力状态，即在-100～100 MPa范围内，每间隔20 MPa设置一个状态点；设置5种弹性模量，在206～214 GPa范围内，每间隔2 GPa设置一个状态点。

模态1为水平弯曲模态，波长约为0.056 12 m，温度应力变化时在20个模态中群速度与相速度变化率最大，对温度应力最敏感。通过振型分析，选择在轨底4号节点施加z方向的激励，4号节点x方向设置信号接收点。激励信号设置成以35 kHz为中心频率、汉宁窗调制的5个周期正弦波。从距离激励位置2 m处开始，每间隔0.02 m设置一个数据采集点，共有128个同步采集点。每个采集点的采样频率为1.4 MHz，采样点数为4 096。2、3、4、5 m处接收到的位移信号时域波形如图5所示。从图5可以看出，随着传播距离的增加，信号的波包峰值降低，波形宽度变大，波形在传播过程中发生了畸变。

(a)距激励位置2 m处

(b)距激励位置3 m处

(c)距激励位置4 m处

(d)距激励位置5 m处图5 4号节点不同位置位移信号时域波形

模态3为扭转弯曲模态，波长约为0.064 81 m，温度应力变化时在20个模态中群速度与相速度变化率排第三位。通过振型分析，选择在轨腰36号节点施加y方向的激励，36号节点x方向设置信号接收点。模态7的振型比较复杂，其波长约为0.077 15 m，相对于其他模态更容易被激励出来。通过振型分析，选择在轨顶41号节点施加x方向的激励，41号节点z方向设置信号接收点。模态20为伸展模态，波长约为0.223 2 m，弹性模量变化时在20个模态中群速度与相速度变化率最大，对弹性模量最敏感。通过振型分析，选择在轨顶64号节点施加x方向的激励，64号节点x方向设置信号接收点。

分别对模态1、模态3、模态7、模态20的128点数据做二维傅里叶变换，在弹性模量210 GPa、应力0 MPa状态下处理结果如图6所示，通过频率和波数计算各模态的相速度约分别为1 960.19、2 252.96、2 731.06、7 510.40 m/s，与理论计算的模态速度相符。

(a)模态1 (b)模态3

(c)模态7 (d)模态20图6 零应力下35 kHz导波不同模态的二维傅里叶变换结果

二维傅里叶变换后的波数结果存在量化误差，大致估算的相速度值可用于模态辨识，应力检测则需要高精度的波数。因此，通过频率和二次加权修正波数计算该模态的波数[11]。根据类似于三角形重心坐标的加权公式，频率-波数加权点位置如图7所示。选取距离激励频率f最近的两个傅里叶变换量化频率f1和f2，选取波数ξ附近3个波数ξ1、ξ2、ξ3，同时可以得到6个点s11、s12、s13、s21、s22、s23对应的等高线幅值l11、l12、l13、l21、l22、l23。

图7 频率-波数加权点位置示意

使用等高线幅值归一化对单个频率下波数进行加权修正得到

(27)

(28)

(29)

3.3 多模态融合应力估计结果

模型标定数据选取5种弹性模量和11种温度应力状态下的波数。假设温度变化1 ℃约等效于弹性模量变化0.05 GPa，室温为20 ℃时，弹性模量为210 GPa。钢轨温度变化范围为-20～60 ℃，为了避免估计应力时产生边缘效应，将温度的变化范围扩大至-60～100 ℃，其等效的弹性模量变化范围为206～214 GPa，每间隔2 GPa设置一个状态点。温度应力变化范围是-100～100 MPa，每间隔20 MPa设置一个状态点。除去弹性模量212 GPa状态下，温度应力值为-20、-40、-60、-80、-100 MPa的数据点(这些数据用于测试)，共有50种状态数据。由于测试样本不足，采用添加白噪声形式获得更多样本，防止数据过度拟合。各模态时域信号引入5%的高斯白噪声，弹性模量引入1%的高斯白噪声，每一个状态波数生成50个样本，一共生成2 500组特征信号。

由于试验条件的限制，目前尚不具备进行现场实物验证的条件(这需要一段锁定状态下应力已知的钢轨以及相应的相控阵采集处理系统)。因此，本文通过仿真获取测量数据。测量数据的产生基于半解析有限元法建立的60 kg/m钢轨模型。为了模拟现场实物试验场景，本文使仿真测量的钢轨模型与训练标定函数时的钢轨模型存在差异：设置参数时，弹性模量引入1%的高斯白噪声；同时，传感器采集的阵列信号存在随机测量噪声，时域信号引入5%的高斯白噪声。这些差异和噪声已达到实际应用可能存在的上限。

模型测量数据选取2种弹性模量和5种温度应力状态下的波数：在弹性模量210 GPa状态下，选取温度应力值10～90 MPa，每间隔20 MPa设置一个状态点；在弹性模量212 GPa状态下，选取温度应力值-100～-20 MPa，每间隔20 MPa设置一个状态点，共计10种状态数据。每一个状态生成50个样本，共计500个测量样本。模型预测应力误差柱状分布如图8所示，总体估计应力误差均值为0.210 1 MPa，标准差为0.726 4 MPa。图9为10种状态的应力估计误差，由于每个状态下应力估计标准差均小于1 MPa，图形的纵坐标刻度为20 MPa，因此对应力估计的误差线做放大20倍处理。图9中最小的应力估计误差是弹性模量210 GPa、拉应力10 MPa处，标准差为0.481 MPa；最大的应力估计误差是弹性模量212 GPa、拉应力-40 MPa处，标准差为0.995 MPa。

图8 BP神经网络预测应力误差柱状分布

图9 模型测量数据估计应力及标准差

4 结束语

无缝钢轨纵向温度应力的检测与监测一直是一个技术难题，尽管近年来提出的超声导波无损检测原理得到了广泛应用，但钢轨弹性模量的不确定性制约着超声导波检测应力的测量精度。钢轨弹性模量的不确定性一方面来源于温度变化，另一方面来自不同批次钢轨的差异性或者材料性能随时间的衰变特性。本文基于钢轨温度应力检测的误差传递模型，分析波数、弹性模量与应力的影响关系，发现不同模态对应力和弹性模量的敏感度存在差异，提出利用多模态导波融合的弹性模型补偿算法，实现弹性模量未知条件下钢轨温度应力的精确估计。仿真试验结果表明，本文提出的多模态融合算法能够克服弹性模量变化对温度应力检测的影响，为进一步研究无缝线路温度应力的测量方法提供依据。

参考文献：

[1]CHEN F,WILCOX P D.The Effect of Load on Guided Wave Propagation[J].Ultrasonics,2007,47(1/4):111-122.

[2]LOVEDAY P W.Semi-analytical Finite Element Analysis of Elastic Waveguides Subjected to Axial Loads[J].Ultrasonics,2009,49(3):298-300.

[3]刘增华,刘溯,吴斌,等.预应力钢绞线中超声导波声弹性效应的试验研究[J].机械工程学报,2010,46(2):22-27.

LIU Zenghua,LIU Su,WU Bin,et al.Experimental Research on Acoustoelastic Effect of Ultrasonic Guided Waves in Prestressing Steel Strand[J].Journal of Mechanical Engineering,2010,46(2):22-27.

[4]BARTOLI I,COCCIA S,PHILLIPS R,et al.Stress Dependence of Ultrasonic Guided Waves in Rails[C]//Proceedings of SPIE:The International Society for Optical Engineering,2010,7650:765021.

[5]许西宁,叶阳升,江成,等.钢轨应力检测中超声导波模态选取方法研究[J].仪器仪表学报,2014,35(11):2473-2483.

XU Xining,YE Yangsheng,JIANG Cheng,et al.Research on Method for Mode Selection of Guided Ultrasonic Waves in Stress Measurement of Rails[J].Chinese Journal of Scientific Instrument,2014,35(11):2473-2483.

[6]LOVEDAY P W,WILCOX P D.Guided Wave Propagation as a Measure of Axial Loads in Rails[C]//Proceedings of SPIE:The International Society for Optical Engineering,2010,7650:765023.

[7]HAYASHI T,SONG W J,ROSE J L.Guided Wave Dispersion Curves for a Bar with an Arbitrary Cross-section,a Rod and Rail Example[J].Ultrasonics,2003,41(3):175-183.

[8]BARTOLI I,MARZANI A,LANZA D S F,et al.Modeling Wave Propagation in Damped Waveguides of Arbitrary Cross-section[J].Journal of Sound and Vibration,2006,295(3/5):685-707.

[9]卢耀荣.无缝线路研究与应用[M].2版.北京:中国铁道出版社,2010:37-48.

[10]许西宁.基于超声导波的无缝线路钢轨应力在线监测技术应用基础研究[D].北京:北京交通大学,2014.

[11]王嵘, 余祖俊,朱力强,等.基于导波速度的无缝钢轨应力检测方法[J].中国铁道科学,2018,39(2):18-27.

WANG Rong,YU Zujun,ZHU Liqiang,et al.Detection Method for Stress in Continuously Welded Rail Based on Guided Wave Velocity[J].China Railway Science,2018,39(2):18-27.