在列车高速运行状态下的车轮冲击及外界复杂环境等因素的影响，铁路轨道几何参数，包括轨距、水平、左/右高低、左/右轨向和三角坑等，会逐渐偏离理想位置，形成轨道不平顺状态。轨道不平顺会直接影响列车安全、平稳运行。铁路部门需要通过维修作业恢复轨道质量，保障行车安全、舒适。随着“状态修” “预防修”的推进，掌握轨道不平顺状态及其未来发展趋势对制订科学、合理的维修作业计划至关重要。

综合检测列车和轨道检测车(以下简称“动轨检车”)是我国铁路工务部门感知轨道状态的重要检测手段之一，基于动轨检车检测数据(以下简称“动轨检数据”)计算轨道质量指数(track quality index，TQI)，实现轨道不平顺状态的量化。在轨道各项几何参数中，左/右高低在维修决策中发挥着重要指导作用，国内外学者采用了不同类型的函数描述轨道高低不平顺劣化模式。文[1-2]基于指数型函数建立了2次机械捣固之间的轨道高低不平顺劣化预测模型；郭然等^[3]提出了2次机械捣固之间25 m单元区段的指数型轨道不平顺发展模型；文[4-6]利用线性函数拟合2次维修之间的轨道高低不平顺劣化过程；Soleimanmeigouni等^[7]建立两阶段分段线性拟合模型来描述轨道几何状态劣化过程；Drozdziel等^[8]基于对数型函数建立了轨道质量状态劣化模型。

探索轨道高低不平顺劣化模式是预测轨道未来发展状态的基础。实际中，受多方因素影响，轨道高低不平顺劣化模式复杂多变。首先，我国地域辽阔，气候复杂，夏季的持续高温、强降雨、东北地区的极寒等极端天气会加剧轨道的劣化；其次，维修作业(如机械捣固、钢轨打磨等)会有效恢复轨道质量状态，同时也会改变轨道刚度和路基支撑刚度，导致轨道高低不平顺劣化率发生变化；再者，动轨检车在检测过程中，不可避免会产生异常检测数据，进而导致计算的TQI异常，干扰对轨道高低不平顺劣化模式的分析。因此，消除异常检测数据的影响，根据维修作业日期划分轨道高低不平顺劣化过程是研究其劣化模式的前提。

为克服异常检测数据的影响，自动感知轨道不平顺状态的变化，识别历史维修作业日期，作者在之前的研究工作中建立了一套动态检测数据驱动的轨道不平顺劣化自适应分段建模方法(minimum-description-length-based rail track deterioration adaptive segmentation framework，MDL-RTDAS)^[9]。在应用中发现其中对于疑似维修作业日期的识别会直接影响维修作业日期的识别准确度和运算效率。为解决这一问题，本文基于Mahalanobis距离提出双层Fast-MCD估计法来改进疑似维修作业日期的识别，同时给出每一层Fast-MCD临界值的确定方法及具体数值。经双层Fast-MCD估计法改进后，维修作业日期的识别准确度和运算效率得到了显著提高。利用改进后的MDL-RTDAS将每个线路基本单元的轨道左/右高低不平顺劣化过程划分为多个维修周期，基于每个维修周期深入探索轨道不平顺劣化模式。

1 改进的轨道不平顺劣化自适应分段建模方法 1.1 MDL-RTDAS

作者基于Rissanen^[10]提出的最优编码理论——最小描述长度(minimum description length，MDL)原则研究建立了MDL-RTDAS，将轨道不平顺劣化过程中的维修作业日期识别问题转换为最优模型选择问题，建立了优化目标方程，并提出相应的求解算法，具体过程简单描述如下。

根据我国《高速铁路线路维修规则》^[11]，TQI为以200 m轨道区段为单元的各项几何参数检测数据的标准差之和。本文以200 m轨道区段作为线路基本单元，每个线路基本单元的左/右高低不平顺分别通过该单元范围内左/右轨高低项检测数据的标准差量化(以下统称为“高低不平顺”，无特别说明不再区分左右)。设在一段时间内动轨检车在某条线路上完成了n次检测，每次检测日期为t₁, t₂, …, t_n；该线路上任意一个线路基本单元在各次检测的高低不平顺值为q₁, q₂, …, q_n，记其集合为q={q₁, q₂, …, q_n}；期间在该单元进行了m次维修作业，根据维修作业日期，将高低不平顺劣化过程分割为m+1个维修周期。设第j个维修周期从第B_j次动态检测开始，该周期内共完成了n_j次动态检测；本文认为各个维修周期内劣化过程相互独立，利用一元线性回归模型分别拟合每个维修周期内的劣化过程，得到q的拟合值$ \hat q$，k_j和b_j分别表示第j个维修周期线性拟合函数的斜率和截距，μ_j和σ_j²分别表示第j个维修周期线性拟合后残差的均值和方差。令L_q表示q的总描述长度，则MDL-RTDAS的优化目标方程为

$ \begin{gathered} \min L_{q}=\ln m+(m+1) \ln n+\frac{n}{2} \ln \sum\limits_{j=1}^{m+1} \sigma_{j}^{2}+ \\ \sum\limits_{j=1}^{m+1} n_{j}+\sum\limits_{j=1}^{m+1} \frac{1}{2 \sigma_{j}^{2}} \sum\limits_{i=B_{j-1}}^{B_{j}-1}\left(q_{i}-\hat{q}_{i}-\mu_{j}\right)^{2}+\frac{n}{2} \ln (2 \pi), \\ 3 \leqslant n_{j}<n, \\ 0 \leqslant m<n, \\ \hat{q}_{i}=k_{j} \times t_{i}+b_{j}, \\ q_{i}-\hat{q}_{i} \sim N\left(\mu_{j}, \sigma_{j}^{2}\right). \end{gathered} $

(1)

考虑到高低不平顺具有数据样本量小且受异常检测干扰较大的特征，MDL-RTDAS提出了一套求解式(1)的方法，共分为2步。

1) 确定疑似维修作业点。

将高低不平顺劣化过程中的异常检测和维修作业后的第1趟检测的检测日期统称为疑似维修作业点。异常检测点和维修作业点的表现特征不同，异常检测点明显偏离于当前劣化趋势，而进行一次维修作业后，高低不平顺状态会被明显改善，不平顺值会下降到标准范围内，当前劣化过程结束，进入到下一个劣化过程。因此，疑似维修作业点所对应的一阶差分值会明显偏低/偏高。基于上述分析，认为当q的一阶差分序列中某点的值大于设定的阈值，认为该点对应的为疑似维修作业点。

2) 确定最优分段拟合模型。

将识别出来的疑似维修作业点进行任意方式的组合，组合方式不同，划分出的维修周期不同，计算得到的L_q也不同，通过动态规划算法求解L_q的最小值，对应的即为最优分段拟合模型。

基于上述分析可知，识别的疑似维修作业点会直接影响算法运行的效率和准确度。通过算法识别的疑似维修作业点应能涵盖所有真实的维修作业点；而识别出来的疑似维修作业点数量过多时，则会导致运行效率降低，搜索最优分段拟合模型所需要的时间过长。在异常检测影响较小、高低不平顺数值变化平稳且波动不大的情况下，通过MDL-RTDAS中的阈值判断方法能够有效识别疑似维修作业点。而实际中，多数线路基本单元的高低劣化过程复杂多变，波动较大，如图 5所示。因此，在节1.2中提出了改进的疑似维修作业点识别方法。

1.2 Fast-MCD

根据上述对疑似维修作业点特征的分析，考虑通过离群点识别算法来改进。在基于距离的离群点识别算法中，Mahalanobis距离具有较强的稳定性和鲁棒性，不受各维数据量纲的影响，排除变量之间相关性的干扰，被广泛应用于实际数据的离群点识别中。

Mahalanobis距离是由印度统计学家Mahalanobis^[12]提出的，反映一个样本与样本总体之间的距离。设在p维空间中的一个样本矢量X=[X₁, X₂, …, X_p]^T，均值向量为μ=[μ₁, μ₂, …, μ_p]^T，其中μ_i为X_i的期望；协方差矩阵S=(s_ij)_p×p，其中s_ij为X_i和X_j的协方差。则某个样本x与样本总体X之间的Mahalanobis距离D_M(x)为

$ {D_{\rm{M}}}(x) = \sqrt {{{\left( {x - \mathit{\boldsymbol{\mu }}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{S}}^{ - 1}}\left( {x - \mathit{\boldsymbol{\mu }}} \right)} . $

(2)

当样本的Mahalanobis距离超过某个范围时，认为该样本为离群点。为避免离群点对样本Mahalanobis距离的干扰，Rousseeuw^[13]提出了最小协方差行列式(minimum covariance determinant，MCD)估计方法，为了提高计算效率，Rousseeuw等^[14]提出了改进的Fast-MCD。在Fast-MCD中，构建初始样本集计算初始均值向量和协方差矩阵时，对于p维样本集要随机抽取p+1个样本组成初始样本集。考虑到初始样本的选择会直接影响Fast-MCD的收敛速度和计算效率，因此，此处利用模糊聚类算法代替随机抽取形成初始样本集。改进后Fast-MCD具体流程见图 1。

1.3 基于双层Fast-MCD的改进MDL-RTDAS

根据节1.1中对疑似维修作业点的定义，疑似维修作业点包含了异常检测点和维修作业点。通过分析实际检测数据，不同趟次的动轨检车产生的异常检测点之间不具有相关性，部分异常检测点数值较小，而高低不平顺值本身具有一定的波动性，不是严格随时间单调递增/递减，难以通过设定固定的阈值等方法将异常检测点识别出来；有些异常检测点数值很大(见图 2)，明显偏离其他数值，严重干扰维修作业点的识别，将其称为明显异常点。

基于上述分析，本节建立双层Fast-MCD算法来识别任一线路基本单元高低不平顺劣化过程中的明显异常点和疑似维修作业点。要实现明显异常点和疑似维修作业点的准确定位，需要确定的关键参数是每一层Fast-MCD的临界值D₀。本节利用昌福高速下行方向2013—2020年所有趟次动轨检数据来分别确定每一层Fast-MCD的临界值。将昌福高速下行方向K22+000~K215+000路段中的岔区去掉后，得到887个位于区间正线的线路基本单元。动轨检车的采样间隔为0.025 m，在2013—2020年期间，动轨检车在该路段范围内上采集了超过15万条检测数据。依托实验室研发的“基于动态检测数据的轨道变形分析系统”完成原始动轨检数据的里程偏差修正、TQI计算等数据预处理工作^[15]。

第1层：识别明显异常点。

将每个线路基本单元从2013—2020年的轨道左/右高低不平顺值作为一个样本集，每个样本处于检测时间-轨道高低不平顺值的二维坐标系下，因此通过模糊聚类将样本集分为3类，通过每一类的中心值构建初始协方差矩阵和初始均值向量，根据节1.2中所述方法进行迭代，达到收敛条件后，根据式(2)计算每个样本在当前样本集的稳健Mahalanobis距离。

令D₀₁表示第一层Fast-MCD的临界值，认为样本在当前样本集中的稳健Mahalanobis距离超过D₀₁时为明显异常点。D₀₁的取值要确保识别出每个线路基本单元的明显异常点，同时不能将维修作业点误判为明显异常点。通过人工查看每个线路基本单元左/右高低不平顺的历史变化情况，记录其中远大于/小于其他趟次检测值的检测日期，即为明显异常点。根据所有明显异常点的稳健Mahalanobis距离进行排序，其中最小值为11.053，当D₀₁=11时，未发现将维修作业点误判为明显异常点，因此D₀₁=11。以某一个线路基本单元的右高低不平顺劣化过程为例进行结果展示(如图 2所示)，图中横坐标为各次检测与第一次检测的间隔天数，纵坐标为各次检测对应的右高低不平顺值；以迭代结束后的稳健均值向量为中心点，根据稳健协方差矩阵绘制置信椭圆，根据D₀₁对置信椭圆进行伸缩变换，则置信椭圆范围之外的认为是明显异常点。

第2层：识别疑似维修作业点。

通过第1层将每个线路基本单元高低不平顺劣化过程中的明显异常点去掉后，下一步要解决的问题是识别疑似维修作业点。根据节1.1中提到的，疑似维修作业点对应的一阶差分值会明显偏高/偏低。将每个线路基本单元的左/右高低不平顺值的一阶差分作为一个样本集，根据节1.2中所述方法进行迭代，达到收敛条件后，根据式(2)计算每个样本在当前样本集的稳健Mahalanobis距离。

令D₀₂表示第二层Fast-MCD的临界值，认为一阶差分序列中稳健Mahalanobis距离超过D₀₂的样本为疑似维修作业点。D₀₂的取值要确保将实际维修作业点全部包含进去，同时D₀₂不能过小，过小则会导致判别出的疑似维修作业点过多，求解最优解需要的时间过长。通过人工查看多个线路基本单元左/右高低不平顺的历史变化情况，记录每个单元的维修作业点。根据所有维修作业点的稳健Mahalanobis距离进行排序，最小值约为3.5。因此，设定D₀₂=3.5。识别的疑似维修作业点包含维修作业点和异常检测点，通过求解式(1)判别维修作业点。为进一步缩小疑似维修作业点的范围，根据异常检测点区别于维修作业点的特征增设一条判别规则：若一阶差分样本集中2个相邻样本点的稳健Mahalanobis距离均大于3.5，且对应的一阶差分值的乘积小于0，则认为该点为异常检测点，不属于疑似维修作业点。以某一个线路基本单元的左高低不平顺的一阶差分为例进行结果展示(如图 3所示)，图 3中横坐标为各次检测与第一次检测的间隔天数，纵坐标为高低不平顺一阶差分值；以迭代结束后的稳健均值向量为中心点，根据稳健协方差矩阵绘制置信椭圆，根据D₀₂对置信椭圆进行伸缩变换，找出置信椭圆范围之外的点，根据增设的判别规则筛选出疑似维修作业点。

2 改进前后MDL-RTDAS对比

通过上述双层Fast-MCD算法对MDL-RTDAS中的疑似维修作业点识别方法进行了改进，本节从运算效率、识别精确度和召回率来对比分析改进前后MDL-RTDAS的性能。本节抽取出357个线路基本单元，利用改进前后MDL-RTDAS分别识别每个线路基本单元的左/右高低不平顺劣化过程中的维修作业点，记录运行时间和识别结果。以改进后的MDL-RTDAS在每个样本集上的运行时间为横坐标，改进前的运行时间为纵坐标，对比结果如图 4所示。经统计，改进前平均运行时间为1.27 s，改进后的平均运行时间为0.63 s，运算效率提升了近50%。

将算法识别出的维修作业点与人工分析标定的实际维修作业点进行对比，若两者一致，则认为是被正确识别的维修作业点。对于每个线路基本单元，通过人工分析其高低劣化过程标定出的维修作业点总个数记为N_real，识别出的维修作业点总个数记为N_total，被正确识别的维修作业点个数记为N_tp，将识别精确度定义为100%×N_tp/N_total，召回率为100%×N_tp/N_real。经统计，改进前平均识别精确度为50.89%，平均召回率为65.28%；而改进后平均识别精确度为96.14%，平均召回率为94.71%。在高低不平顺波动幅度较大的情况下，改进前的识别精确度及召回率较低，而改进后的较高。以2个线路基本单元的高低不平顺劣化过程为例对比改进前后识别结果，如图 5和6所示。综上，改进后MDL-RTDAS在运算效率、识别精确度和召回率上都显著优于改进前的。

3 轨道高低不平顺劣化模式

掌握轨道高低不平顺劣化模式对于预测其未来劣化状态十分重要。根据文献调研，国内外学者常用于描述轨道高低不平顺劣化模式的函数有3类：线性、指数和对数函数。为了确定哪类函数更贴近轨道高低不平顺劣化模式，通过改进后的MDL-RTDAS识别出每个线路基本单元高低不平顺劣化过程的维修作业点，将劣化过程划分为若干个维修周期，利用上述3类函数对每个维修周期内的劣化过程进行拟合及预测，通过对比拟合效果和预测能力确定最优拟合函数。

3.1 拟合效果对比

为了便于计算，将每个维修周期的第一趟检测的时间T₁记为1，之后第i次检测对应的时间T_i为该趟检测与第一趟检测的时间差(以d为单位)。用于拟合的线性、对数和指数函数表达式如式(3)—(5)所示。将高低不平顺值进行相应的尺度变换后，利用最小二乘法求得线性函数中的参数k_linear和b_linear、对数函数中的参数k_log和b_log以及指数函数中的参数k_exp和b_exp。

$ {v_{{\rm{linear}}}}\left( T \right) = {k_{{\rm{linear}}}} \times T + {b_{{\rm{linear}}}}, $

(3)

$ {v_{{\rm{log}}}}\left( T \right) = {k_{{\rm{log}}}} \times {\rm{ln}}\;T + {b_{{\rm{log}}}}, $

(4)

$ {v_{{\rm{exp}}}}\left( T \right) = {{\rm{e}}^{{k_{{\rm{exp}}}} \times T + {b_{{\rm{exp}}}}}}. $

(5)

其中：v_linear(T)、v_log(T)、v_exp(T)分别为线性、对数和指数函数得到的高低不平顺拟合值。以每个维修周期为单位，计算利用3类模型拟合后的均方根误差和平均绝对误差。经汇总统计，线性、对数、指数拟合模型的平均绝对误差的均值分别为0.044 mm、0.055 mm、0.045 mm，均方根误差的均值分别为0.056 mm、0.069 mm、0.057 mm。在所有的维修周期中，线性拟合的平均绝对误差小于对数拟合的占比72.968%，小于指数拟合的占比56.611%；线性拟合的均方根误差小于对数拟合的占比74.587%，小于指数拟合的占比69.089%。

平均绝对误差、均方根误差数值越小，认为函数拟合效果越好。通过上述对比，可得对数函数拟合模型的2个指标值均是最大的，线性函数拟合模型的2个指标值略小于指数函数拟合模型的。因此，得到结论：对数函数拟合效果最差，线性函数与指数函数的拟合效果较好，线性函数的拟合效果略优于指数函数。

3.2 预测能力对比

为进一步分析轨道高低不平顺劣化模式更贴近哪类函数，将每个维修周期中的数据分为2部分，一部分数据用于求解模型的参数，另一部分数据用于与模型预测值进行对比，进而比较各类模型的预测能力。本文将每个维修周期中前3/4的轨道高低不平顺值来求解模型参数，利用求解出的模型预测后1/4轨道高低不平顺值，与实际数据进行比较，计算均方根误差和平均绝对误差，对比预测能力。经汇总统计，线性、对数、指数模型预测的平均绝对误差的均值分别为0.048 mm、0.058 mm、0.051 mm，均方根误差的均值分别为0.062 mm、0.073 mm、0.066 mm。在所有的维修周期中，线性模型预测的平均绝对误差小于对数模型预测的占比63.770%，小于指数模型预测的占比62.859%；线性模型预测的均方根误差小于对数模型预测的占比62.316%，小于指数模型预测的占比66.836%。以某个线路基本单元为例，用3类函数分别对每个维修周期内的高低不平顺劣化过程进行拟合和预测，如图 7所示。

预测结果的平均绝对误差、均方根误差越小，认为模型预测能力越强。通过上述对比，可得线性模型的预测值的均方根误差和平均绝对误差均小于其他2类模型。因此，线性模型的预测值与实际值更接近，预测能力更强。

综上，线性函数在拟合效果和预测能力上都优于对数函数和指数函数，因此认为线性函数更贴近轨道高低不平顺劣化模式。同时，《高速铁路线路维修规则》中对允许速度大于250 km/h的线路区段的左/右高低项TQI管理标准值为0.8~0.9 mm^[11]。文中计算出的线性函数拟合平均绝对误差的最大值为0.35 mm，平均值为0.044 mm；预测平均绝对误差的最大值为0.3 mm，平均值为0.048 mm，均远小于左/右高低项TQI管理标准值，认为线性函数拟合及预测的精度满足日常维修工作。

4 结论

由于动态检测数据驱动的轨道不平顺劣化自适应分段建模方法中疑似维修作业点的识别精确度较低，直接影响该方法识别维修作业点的精确度和效率。为解决这一问题，本文基于Mahalanobis距离提出了双层Fast-MCD对其进行改进，通过昌福高速铁路下行方向2013—2020年的动态检测数据确定出每一层Fast-MCD的临界值。经验证，改进后的MDL-RTDAS在运算效率和识别精确度上都得到了显著提升。

通过改进后的MDL-RTDAS识别出昌福高速铁路下行方向上位于区间正线的线路基本单元的维修作业日期，将左/右高低不平顺劣化过程划分为多个维修周期。基于每个维修周期内的劣化过程对比线性函数、指数函数和对数函数的拟合效果和预测能力，认为线性函数更贴近轨道高低不平顺劣化模式。