基于惯性基准的线路沉降车载检测方法

Research on On-Board Detection Method for Track Settlement Based on Inertial Reference

摘要

为提升线路沉降检测效率，提出基于惯性基准的线路沉降车载检测方法。基于测量原理推算误差来源，采用角速度测量法测量轨道高程，用增量求和法预处理IMU数据，基于列车行驶过程中的动态零速假设建立里程计观测修正模型，基于扩展卡尔曼滤波框架实时估计姿态及陀螺仪零偏，融合惯性推算结果与位移测量值计算轨面高程和高低不平顺，建立局部沉降快速筛查与区段高程测量相结合的检测流程，并在不同线路开展试验验证。结果表明：该方法测量轨道超高重复性偏差的95%分位数最大值为0.9 mm，重复性偏差均值为0.24 mm，对姿态测量准确性的提升较为明显；350 km · h^-1速度等级下测量的1 000 m长波高低不平顺重复性偏差的95%分位数为1.57 mm，均值为0.58 mm；120 km · h^-1速度等级下测量的风险沉降区段的1 000 m高低不平顺最大重复性偏差为0.5 mm，相对高程的最大重复性偏差为5 mm；80 km · h^-1检测速度下测量的400 m内相对高程与绝对测量的最大偏差为1.87 mm，能够满足线路沉降区段的车载快速测量的需求。

Abstract

To improve the efficiency of line settlement detection, an on-board detection method for track settlement based on inertial reference is proposed. The sources of error are deduced based on the measurement principle, and the angular velocity measurement method is proposed to measure. The IMU data are preprocessed by the incremental summation method. The odometer observation correction model is established by utilizing the dynamic zero-speed assumption during train travel, and the attitude and gyroscope zero-bias are estimated in real-time based on the extended Kalman filter. Then, the fusion of inertial results and displacement measurements is used to calculate the track elevation and the vertical irregularity. A detection method combining rapid screening of local settlement and section elevation measurement is built, and tests are conducted on different lines. The results show that the 95% value of the repeatability deviation for measuring track superelevation is 0.9 mm, with a mean value of 0.24 mm, representing a more significant improvement in attitude measurement accuracy. The 95% quartile of the repeatability deviation for 1,000 m long-wave vertical irregularity measured at a speed of 350 km · h^-1 is 1.57 mm, with a mean value of 0.58 mm; the maximum repeatability deviation of the 1,000 m elevation and vertical irregularity of the risky settlement section measured at a speed of 120 km · h^-1 is 0.5 mm, and the maximum repeatability deviation of the relative elevation is 5 mm; the maximum deviation of the relative elevation within 400 m measured at a detection speed of 80 km·h^-1 from absolute measurement is 1.87 mm. These results demonstrate that the proposed method can meet the requirements for on-board rapid measurements in track settlement sections.

Graphical abstract

关键词

线路 / 沉降 / 惯性基准 / 相对高程 / 车载检测 / 扩展卡尔曼滤波

Key words

Railroad line / Settlement / Inertial reference / Relative elevation / On-board detection / Extended Kalman filter

在铁路运营过程中，路基沉降变形不可避免，严重的局部变形会增大轨道不平顺，且导致无砟轨道结构层间离缝、上拱和脱空等病害，进而影响行车安全性和舒适性^［1］。线路沉降监测对保障行车安全、研究沉降发展规律有重要意义^［2-3］。

静态检测是识别与监测线路沉降的主要手段，包括大地水准测量、定点监测等^［4］。其中，大地水准测量应用较多，技术成熟，但连续测量的效率低下。定点监测需在基床表层布设观测桩和沉降板，既干扰施工且易受破坏^［5］。轻型轨检小车发展迅速，且应用广泛，配备使用最高等级全站仪（0.5″级全站仪）的轨检仪，能达到1.5 mm的绝对位置精度^［6］。轨检仪虽可实现轨面高程连续测量，但同样面临测量效率较低、过于依赖CP Ⅲ控制点的问题^［7］。GNSS监测点安装灵活，但垂向测量精度较差，且不适用于山区、隧道等^［8］。D-InSAR技术适用于大面积连续区域，但测量周期长，且对离散点绝对沉降量监测能力不足^［9］。搭载于综合检测列车的轨道几何检测系统可对轨道形变进行周期性检测，并对其几何状态进行评判。若能够提升动态检测高低不平顺的波长范围，则可以实现对线路沉降区段的高效普查，从而作为现有线路沉降识别与监测的重要补充手段。

既有研究能够实现对200 m轨道长波不平顺进行测量^［10］，尚难以满足动辄几百米的沉降区段的监测需求。铁路沿线环境复杂，动态测量缺少外部参考值作修正，无法实现对轨道绝对高程的测量，是限制动态测量结果在现有沉降监测体系中发挥作用的主要因素。近年来，以惯性传感器为核心的测量手段应用至轨道及公路测量^［11］、管线测量^［12］、采煤机位置测量^［13-14］等方面，根据不同应用场景构建观测约束条件，能够提升无外部修正条件下的惯性测量的精度保持能力，得到的相对测量结果同样能满足实际应用需求。

本文提出1种基于惯性基准相对测量的线路沉降车载检测方法。首先，用增量求和法将惯性测量单元（Inertial Measurement Unit，IMU）数据转换至空间序列；其次，基于列车的动态零速假设建立了速度观测修正模型；然后，基于角度积分和位移补偿的方法计算轨面高程；最后，将局部沉降快速筛查与区段高程测量2种方式相结合，实现线路局部沉降的动态检测。

1 检测需求及原理

1.1　线路沉降检测需求

铁路线路沉降包含路基工后沉降、路基与桥涵等结构过渡区域差异沉降、无砟轨道及桥梁等工程结构变形导致的沉降等多种不同情况，最终均会不同程度地反映为轨道几何形态的变化^［15］。其中，路基工后的不均匀沉降主要包括施工期残余沉降和运营期列车动载引起的动变形等，所引起轨道不平顺的波长覆盖几十米至几公里。路桥过渡段沉降主要是因为路基和桥梁的刚度差异明显，在过渡点附近容易产生变形差而导致轨面弯折，长度一般为几十米。桥梁结构形变与其跨度相关，除少数大跨度桥梁外，大多数桥梁的跨度在120 m内。由于大多数线路沉降变形在1 km范围内，因此在养护维修环节，研究重点应放在1 km范围内沉降区段累计变形的监测与评估上。

1.2　动态检测原理及误差分析

轨道动态检测系统搭载于动态检测列车平台，可实时测量轨道几何参数，若能够改善低信噪比下的测量精度，提升轨道高低不平顺的测量范围，可实现线路沉降的快速筛查。

图1为轨道检测系统结构示意图。可见，视觉测距单元与惯性测量单元均安装在检测梁体上，检测梁悬挂于构架端部，编码器安装于轮对一侧。定义IMU本体坐标系o_bx_by_bz_b，其三轴正向分别为前、左、上，原点为IMU中心；定义导航坐标系o_nx_ny_nz_n，其三轴正向分别为东、北、上，原点为起始测量位置；定义载体坐标系o_wx_wy_wz_w，其x轴沿轨道纵向方向，y轴沿轨平面向左，z轴垂直于轨平面向上，原点与本体坐标系原点重合。

传统惯性基准法的高低不平顺计算式可简化为^［16］

\widehat{Z}=\iint (a_{\mathrm{m}\mathrm{b}z}^{}-\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{c}\mathrm{o}\mathrm{s}\psi G_{\mathrm{0}})\mathrm{d}t\mathrm{d}t+d_{\mathrm{V}}（1）

式中：\widehat{Z}为高低不平顺；a_{\mathrm{m}\mathrm{b}z}^{}为z轴加速度计的测量值；\theta和\psi分别为测量平台侧滚角和俯仰角；G_{\mathrm{0}}为重力常量；t为时间；d_{\mathrm{V}}为垂向位移计的测量值。

由式（1）可知垂向加速度计是核心传感器之一，该传感器的测量误差可表示为

a_{\mathrm{m}\mathrm{b}z}^{}=a_{\mathrm{b}z}^{}+{\kappa }_z+W_{\mathrm{a}z}（2）

式中：a_{\mathrm{b}z}^{}为本体坐标系下垂向加速度；{\kappa }_z为z轴加速度计的零偏；W_{\mathrm{a},z}为z轴加速度计的测量噪声。

加速度计零偏所导致的高程测量误差e可以表示为

e=\frac{\mathrm{1}}{\mathrm{2}}{\kappa }_z{t}^{\mathrm{2}}（3）

若行车速度为360 km · h^-1，基于不同零偏加速度计测量的高程积分漂移随里程变化曲线如图2所示。由图2可知：在加速度计零偏的作用下，高程积分漂移迅速增长至厘米级，难以满足应用需求。

文献［16］对比了加速度测量和角速度测量2种方法，在现有传感器选型下，角速度测量法在测量长波不平顺时具有测量信噪比更高的优势。据此，文献［10］中通过观测轨平面角位移变化实现了200 m高低不平顺的测量。在用该方法计算高程时，与加速度计相似，陀螺仪零偏经积分所产生的误差构成了主要的测量误差。因此，准确估计传感器零偏，提升平台的姿态测量精度，是实现高程精准测量的关键。

2 线路沉降检测方法

2.1　数据预处理

为保证测量结果能够按照等距离输出，需将固定频率输出的IMU数据转换至空间相关性序列，本文采用增量求和法进行处理。

系统以每行驶固定距离\mathrm{\Delta }x（本文\mathrm{\Delta }x取0.25 m）进行1次采样和计算，提取2个采样时间间隔内的IMU数据序列，以计算单个采样周期内的角增量和速度增量。以加速度计为例，采样的时序关系如图3所示。图中：k代表第k次采样，k=2，…，N₀，其中N₀为总采样次数；k-1次至k次采样时刻内共接收了N组加速度数据（ a₁， a₂，……， a_N ），{\mathit{a}}_{\mathrm{0}}为{\mathit{a}}_{\mathrm{1}}的前1组数据，{\mathit{a}}_{N+\mathrm{1}}为{\mathit{a}}_N的后1组数据，1组数据包含3个轴向加速度；\mathrm{\delta }{t}_{\mathrm{1}}为第k-1次采样时刻与前1组加速度接收时刻的时间差；\mathrm{\delta }{t}_{\mathrm{2}}为k次采样时刻与前1组加速度接收时刻的时间差。

则第k个采样时刻的角增量向量\mathrm{\Delta }{\mathit{\theta }}_k和速度增量向量\mathrm{\Delta }{\mathit{v}}_k可以表示为

\mathrm{\Delta }{\mathit{\theta }}_k={\displaystyle \sum _{i=\mathrm{1}}^N}\left({\mathit{\omega }}_i\mathrm{\Delta }t\right)+\mathrm{\delta }{t}_{\mathrm{2}}{\mathit{\omega }}_{N+\mathrm{1}}+\left(\mathrm{\Delta }t-\mathrm{\delta }{t}_{\mathrm{1}}\right){\mathit{\omega }}_{\mathrm{0}}（4）

\mathrm{\Delta }{\mathit{v}}_k={\displaystyle \sum _{i=\mathrm{1}}^N}\left({\mathit{a}}_i\mathrm{\Delta }t\right)+\mathrm{\delta }{t}_{\mathrm{2}}{\mathit{a}}_{N+\mathrm{1}}+\left(\mathrm{\Delta }t-\mathrm{\delta }{t}_{\mathrm{1}}\right){\mathit{a}}_{\mathrm{0}}（5）

式中：{\mathit{\omega }}_i为三轴陀螺仪输出的角速度向量；{\mathit{a}}_i为加速度计输出的加速度向量；角标为IMU数据编号；\mathrm{\Delta }t为IMU数据接收的固定时间间隔。

将角增量向量和速度增量向量分别除以2次采样的时间间隔\mathrm{\Delta }{T}_k，可得到平均加速度向量和平均角速度向量。由于列车速度不断变化，该计算过程可认为是IMU数据经过了1个可变窗长的矩形窗滤波器，窗长近似为\mathrm{\Delta }{T}_k/\mathrm{\Delta }t。参考文献［17］中的信号处理过程，画出不同速度下滤波器的幅频特性曲线，如图4所示。

由图4可知：滤波过程空间域截止波长约为1 m；滤波器的零点位于0.25 m，与采样间距相等；滤波过程的幅频响应并不随速度变化而变化，保证了不同速度下检测结果的一致性。

经过上述预处理之后，IMU数据得以转换至空间相关性序列。

2.2　惯性基准测量方法

为准确估计传感器的零偏，提升姿态测量精度，本文利用列车行驶过程中动态零速假设，推算系统观测与修正方程，然后基于扩展卡尔曼滤波算法框架，实现里程计与IMU的组合姿态测量。

2.2.1　动态零速修正测量模型

扩展卡尔曼滤波算法主要包括3个步骤：状态更新、观测修正和误差反馈。状态更新通过状态方程进行推导；观测修正根据观测方程计算误差状态的估计值；误差反馈则将误差状态回馈到惯性递推的状态值中。

1）状态更新

状态更新过程需定义状态量，并得到状态量的更新方程。定义系统待估计的状态量为姿态角、陀螺仪零偏和速度。参考文献［11-12］，状态方程可表示为

{\dot{\mathit{C}}}_{\mathrm{n}\mathrm{b}}={\mathit{C}}_{\mathrm{n}\mathrm{b}}{\left({\mathit{\omega }}_{\mathrm{m}}-{\mathit{\epsilon }}_{\mathrm{b}}\right)}^{\wedge }（6）

{\dot{\mathit{\epsilon }}}_{\mathrm{b}}=\mathrm{0}（7）

{\dot{\mathit{v}}}_{\mathrm{n}}={\mathit{C}}_{\mathrm{n}\mathrm{b}}{\mathit{a}}_{\mathrm{m}}-{\mathit{G}}_{\mathrm{n}}（8）

其中，

{\mathit{G}}_{\mathrm{n}}={\left[\mathrm{0}\mathrm{0}{G}_{\mathrm{0}}\right]}^{\mathrm{T}}

式中：{\mathit{C}}_{\mathrm{n}\mathrm{b}}为姿态矩阵；{\mathit{\omega }}_{\mathrm{m}}为三轴陀螺仪输出的角速度向量；{\mathit{\epsilon }}_{\mathrm{b}}为陀螺仪的零偏向量；{\left(·\right)}^{\wedge }为反对称矩阵；{\mathit{v}}_{\mathrm{n}}为导航坐标系下的速度向量；{\mathit{a}}_{\mathrm{m}}为三轴加速度计输出的加速度向量；{\mathit{G}}_{\mathrm{n}}为导航坐标系下的重力向量。

在惯性测量中，通过泰勒级数展开将状态方程进行局部线性化，用误差状态值进行迭代计算，可减小非线性系统离散化过程所产生的误差，有助于提升状态估计的精度。设定误差状态向量\mathrm{\delta }\mathit{\xi }表示为

\mathrm{\delta }\mathit{\xi }=(\mathrm{\delta }{\mathit{\phi }}^{\mathrm{T}}\mathrm{\delta }{\mathit{\epsilon }}_b^{\mathrm{T}}\mathrm{\delta }{\mathit{v}}_n^{\mathrm{T}}{)}^{\mathrm{T}}（9）

式中：\mathrm{\delta }\mathit{\phi }为状态积分递推过程产生的失准角向量；\mathrm{\delta }{\mathit{\epsilon }}_{\mathrm{b}}为零偏误差向量；\mathrm{\delta }{\mathit{v}}_{\mathrm{n}}为速度误差向量。

误差状态与真实状态的关系式为

{\mathit{C}}_{\mathrm{n}\mathrm{b},\mathrm{t}}^{}=\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{\delta }\mathit{\phi }){\mathit{C}}_{\mathrm{n}\mathrm{b}}（10）

{\mathit{\epsilon }}_{\mathrm{b},\mathrm{t}}={\mathit{\epsilon }}_{\mathrm{b}}+\mathrm{\delta }{\mathit{\epsilon }}_{\mathrm{b}}（11）

{\mathit{v}}_{\mathrm{n},\mathrm{t}}^{}={\mathit{v}}_{\mathrm{n}}-\mathrm{\delta }{\mathit{v}}_{\mathrm{n}}（12）

式中：{\mathit{C}}_{\mathrm{n}\mathrm{b},\mathrm{t}}^{}为真实姿态矩阵；{\mathit{\epsilon }}_{\mathrm{b},\mathrm{t}}为真实零偏向量；{\mathit{v}}_{\mathrm{n},\mathrm{t}}^{}为真实速度向量。

将式（10）—式（12）带入到式（6）—式（8），推导误差状态的微分方程，离散化后的迭代更新方程可表示为

\mathrm{\delta }{\mathit{\xi }}_k^{}={\mathit{F}}_{k-\mathrm{1}}\mathrm{\delta }{\mathit{\xi }}_{k-\mathrm{1}}+{\mathit{V}}_{k-\mathrm{1}}（13）

其中，

F_{k-\mathrm{1}}=\left(\begin{array}{ccc}{\mathit{I}}_{\mathrm{3}}& -{\mathit{C}}_{\mathrm{n}\mathrm{b},k-\mathrm{1}}^{}\mathrm{\Delta }{T}_k& {\mathit{O}}_{\mathrm{3}}\\ {\mathit{O}}_{\mathrm{3}}& {\mathit{I}}_{\mathrm{3}}& {\mathit{O}}_{\mathrm{3}}\\ ({\mathit{C}}_{\mathrm{n}\mathrm{b},k-\mathrm{1}}^{}{\mathit{a}}_{\mathrm{m},k}{)}^{\wedge }\mathrm{\Delta }{T}_k& {\mathit{O}}_{\mathrm{3}}& {\mathit{I}}_{\mathrm{3}}\end{array}\right)

V_{k-\mathrm{1}}=\left(\begin{array}{c}{\mathit{\tau }}_{\mathrm{1},k-\mathrm{1}}\\ \begin{array}{l}{\mathit{\tau }}_{\mathrm{2},k-\mathrm{1}}\\ {\mathit{\tau }}_{\mathrm{3},k-\mathrm{1}}\end{array}\end{array}\right)

式中：\mathrm{\delta }{\mathit{\xi }}_k和\mathrm{\delta }{\mathit{\xi }}_{k-\mathrm{1}}分别为第k次和第k-1次采样时刻的误差状态向量；{\mathit{F}}_{k-\mathrm{1}}为第k-1次采样时刻的状态转移矩阵；{\mathit{V}}_{k-\mathrm{1}}为第k-1次采样时刻的噪声向量；{\mathit{O}}_{\mathrm{3}}为3×3的零矩阵； I₃为3×3的单位矩阵；{\mathit{\tau }}_{\mathrm{1},k-\mathrm{1}}为第k-1次采样时刻的角度预测噪声向量；{\mathit{\tau }}_{\mathrm{2},k-\mathrm{1}}为第k-1次采样时刻的零偏预测噪声向量；{\mathit{\tau }}_{\mathrm{3},k-\mathrm{1}}为第k-1次采样时刻的速度预测噪声向量。若无特殊说明，角标k和k-1均表示第k次和第k-1次采样时刻的物理量。

当k=1时，误差状态量初始化为0，初始姿态矩阵由初始对准过程得到。参考文献［14］，扩展卡尔曼滤波算法一般用预测协方差矩阵表示噪声带来的影响，根据式（13）推算第k-1次采样时刻的预测协方差矩阵{\mathit{Q}}_{k-\mathrm{1}}为

{\mathit{Q}}_{k-\mathrm{1}}=\left(\begin{array}{ccc}{\sigma }_{\mathrm{N},\mathrm{g}}^{\mathrm{2}}\mathrm{\Delta }{T}_k{\mathit{I}}_{\mathrm{3}}& & \\ & {\sigma }_{\mathrm{K},\mathrm{g}}^{\mathrm{2}}\mathrm{\Delta }{T}_k{\mathit{I}}_{\mathrm{3}}& \\ & & {\sigma }_{\mathrm{N},\mathrm{a}}^{\mathrm{2}}\mathrm{\Delta }{T}_k{\mathit{I}}_{\mathrm{3}}\end{array}\right)（14）

式中：{\sigma }_{\mathrm{N},\mathrm{g}}^{}为陀螺仪的角度随机游走参数；{\sigma }_{\mathrm{K},\mathrm{g}}^{}为陀螺仪的角速度随机游走参数；{\sigma }_{\mathrm{N},\mathrm{a}}^{}为加速度计的速度随机游走参数。

通过使用Allan方差对IMU进行标定，可以确定随机噪声参数的数量级，并根据实际效果对参数进行优化调整。

2）观测修正

观测修正过程需得到观测量对应的观测方程，并参与扩展卡尔曼滤波计算。列车在轨道的约束下向前行驶时，忽略轮轨间的蠕滑效应，载体前向速度与里程计测速近似一致。列车在行进过程中航向角速率较小，测量平台与转向架固定安装，若忽略转向架的蛇形游走、横向摆动和沉浮运动的影响，那么横向和垂向的速度近似为零。利用该近似可获得载体坐标系下的三向速度，称之为动态零速假设。

根据动态零速假设，里程计测速的观测方程表示为

{\mathit{v}}_{\mathrm{w},k}\approx {\mathit{C}}_{\mathrm{w}\mathrm{b}}^{}{\mathit{C}}_{\mathrm{n}\mathrm{b},k}^{\mathrm{T}}{\mathit{v}}_{\mathrm{n},k}+{\mathit{W}}_{\mathrm{w}\mathrm{v},k}（15）

其中，

{\mathit{v}}_{\mathrm{w},k}=(v_{\mathrm{f},k}{\mathrm{0}\mathrm{0})}^{\mathrm{T}}

式中：{\mathit{v}}_{\mathrm{w},k}为里程计测量的载体坐标系下的速度向量；v_{\mathrm{f},k}为速度测量值；{\mathit{W}}_{\mathrm{w}\mathrm{v},k}为速度噪声向量；{\mathit{C}}_{\mathrm{w}\mathrm{b},k}为本体坐标系到载体坐标系下的旋转矩阵，由检测梁与轮对间的安装角所决定，需要准确标定。

对观测方程式（15）作误差分析，推算第k次采样时刻的残差，即

{\tilde{\mathit{z}}}_k=-{\mathit{C}}_{\mathrm{w}\mathrm{b}}{\mathit{C}}_{\mathrm{n}\mathrm{b},k-\mathrm{1}}^{\mathrm{T}}({\mathit{v}}_{\mathrm{n},k-\mathrm{1}}{)}^{\wedge }\mathrm{\delta }{\mathit{\theta }}_{k-\mathrm{1}}+{\mathit{C}}_{\mathrm{w}\mathrm{b}}{\mathit{C}}_{\mathrm{n}\mathrm{b},k-\mathrm{1}}^{\mathrm{T}}\mathrm{\delta }{\mathit{v}}_{\mathrm{n},k-\mathrm{1}}+{\mathit{W}}_{\mathrm{w}v}（16）

式中：{\tilde{\mathit{z}}}_k为观测残差向量。

式（16）对\mathrm{\delta }{\mathit{\xi }}_{k-\mathrm{1}}求偏导，得到的雅克比矩阵 H_k-1为

{\mathit{H}}_{k-\mathrm{1}}=(-{\mathit{C}}_{\mathrm{w}\mathrm{b}}\begin{array}{ccc}{\mathit{C}}_{\mathrm{n}\mathrm{b},k-\mathrm{1}}^{\mathrm{T}}({\mathit{v}}_{\mathrm{n}k-\mathrm{1}}{)}^{\wedge }& {\mathit{O}}_{\mathrm{3}}& {\mathit{C}}_{\mathrm{w}\mathrm{b}}{\mathit{C}}_{\mathrm{n}\mathrm{b},k-\mathrm{1}}^{\mathrm{T}}\end{array})（17）

在得到观测方程的雅可比矩阵与观测残差计算式之后，即可代入至扩展卡尔曼滤波算法框架求解\mathrm{\delta }{\mathit{\xi }}_k的估计值，完成观测修正过程。

3）误差反馈

根据式（10）—式（12），将估计的误差状态反馈至递推状态，即可实现1次当前姿态的准确估计。

2.2.2　观测性分析

观测方程对不同状态量的观测性并不一致，对系统展开观测性分析，有利于评估计算过程中引入的误差项^［18］。将本体坐标系下的速度转换至导航坐标系，可表示为

{\mathit{v}}_{\mathrm{n}}={\mathit{C}}_{\mathrm{n}\mathrm{b}}^{}{\mathit{v}}_{\mathrm{b}}（18）

对其求导可表示为

{\dot{\mathit{v}}}_{\mathrm{n}}={\mathit{C}}_{\mathrm{n}\mathrm{b}}^{}{\left({\mathit{\omega }}_{\mathrm{b},\mathrm{n}\mathrm{b}}^{}\right)}^{\wedge }{\mathit{v}}_{\mathrm{b}}+{\mathit{C}}_{\mathrm{n}\mathrm{b}}^{}{\dot{\mathit{v}}}_{\mathrm{b}}（19）

式中：{\mathit{\omega }}_{\mathrm{b},\mathrm{n}\mathrm{b}}^{}为本体坐标系相对于导航坐标系的角速度。

若忽略安装角偏差，根据动态零速假设有{\mathit{v}}_{\mathrm{b}}\approx {\mathit{v}}_{\mathrm{w}}，将导航坐标系的加速度转换至本体坐标系可表示为

{\mathit{a}}_{\mathrm{b}}={\mathit{C}}_{\mathrm{n}\mathrm{b}}^{\mathrm{T}}{\dot{\mathit{v}}}_{\mathrm{b}}（20）

则有如下推导

\begin{array}{l}{\mathit{a}}_{\mathrm{b}}={\left({\mathit{\omega }}_{\mathrm{n},\mathrm{n}\mathrm{b}}^{}\right)}^{\wedge }{\mathit{v}}_{\mathrm{b}}+{\dot{\mathit{v}}}_{\mathrm{b}}\\ ={\left({\begin{array}{ccc}{\dot{v}}_{\mathrm{f}}& {\omega }_z{v}_{\mathrm{f}}& -{\omega }_{y}v\end{array}}_{\mathrm{f}}\right)}^{\mathrm{T}}\end{array}（21）

式中：v_{\mathrm{f}}为里程计测量的速度；{\omega }_y为{\mathit{\omega }}_{\mathrm{b},\mathrm{n}\mathrm{b}}^{}的y轴分量；{\omega }_z为{\mathit{\omega }}_{\mathrm{b},\mathrm{n}\mathrm{b}}^{}的z轴分量。

加速度计的测量值可分解为

{\mathit{a}}_{\mathrm{m}}^{}={\mathit{a}}_{\mathrm{b}}+{\mathit{C}}_{\mathrm{n}\mathrm{b}}^{\mathrm{T}}{\mathit{G}}_{\mathrm{n}}+\mathit{\kappa }（22）

式中：\mathit{\kappa }为加速度计的零偏向量。

将式（21）代入到式（22）并展开可得

\left\{\begin{array}{l}-a_{\mathrm{m}\mathrm{b}x}^{}+{\dot{v}}_{\mathrm{f}}-G_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}\psi +{\kappa }_x=\mathrm{0}\\ -a_{\mathrm{m}\mathrm{b}y}^{}+v_{\mathrm{f}}{\omega }_z+G_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}\psi \mathrm{s}\mathrm{i}\mathrm{n}\theta +{\kappa }_y=\mathrm{0}\\ -a_{\mathrm{m}\mathrm{b}z}^{}-v_{\mathrm{f}}{\omega }_y+G_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{c}\mathrm{o}\mathrm{s}\psi +{\kappa }_z=\mathrm{0}\end{array}\right.（23）

式中：a_{\mathrm{m}\mathrm{b}x}^{}为x轴加速度计的测量值；a_{\mathrm{m}\mathrm{b}x}^{}为y轴加速度计的测量值；{\kappa }_x，{\kappa }_y和{\kappa }_z分别为x轴、y轴和z轴加速度计的零偏；\psi为俯仰角；\theta为侧滚角。

由式（23）可以看出动态零速修正模型对各状态量的观测效果。由于地球重力常量已知，该方法可对水平姿态角直接观测，而对航向角不可观测。在滤波估计的状态量反馈更新过程中，为避免对航向角作错误修正，设定航向角误差不反馈。

2.2.3　轨道高程与长波高低不平顺测量

轨面高程可通过惯性测量参数与视觉测距单元测量值结合计算而得到。可通过引入约束点高程并进行平差处理提升高程的测量精度，降低惯性参数积分所产生的漂移。

基于航迹递推原理，推算IMU的位移迭代更新计算式为

{\mathit{p}}_k={\mathit{p}}_{k-\mathrm{1}}+{\mathit{C}}_{\mathrm{n}\mathrm{b}}{\left(\begin{array}{ccc}\mathrm{\Delta }x& \mathrm{0}& \mathrm{0}\end{array}\right)}^{\mathrm{T}}（24）

式中：{\mathit{p}}_k为第k次采样时刻IMU所在的位置向量。

由式（24）提取IMU的高程{\tilde{z}}_{\mathrm{I},k}，可表示为

{\tilde{z}}_{\mathrm{I},k}={\displaystyle \sum _{j=\mathrm{1}}^k}(\mathrm{s}\mathrm{i}\mathrm{n}{\psi }_j\mathrm{\Delta }x)（25）

式（25）表明IMU的高程可近似为俯仰角随距离的积分。结合视觉传感器测量的位移信息，并根据传感器间空间关系推算左、右轨高程表示为

\left\{\begin{array}{l}{\tilde{z}}_{\mathrm{l},k}={\tilde{z}}_{\mathrm{I},k}+\mathrm{\Delta }z+d_{\mathrm{V},\mathrm{l},k}\\ {\tilde{z}}_{\mathrm{r},k}={\tilde{z}}_{\mathrm{I},k}-\mathrm{\Delta }z+d_{\mathrm{V},\mathrm{r},k}\end{array}\right.（26）

式中：{\tilde{z}}_{\mathrm{l},k}和{\tilde{z}}_{\mathrm{r},k}分别为左右轨高程，角标l和r分别表示左轨和右轨；\mathrm{\Delta }z为侧滚方向的补偿修正项。

在缺少绝对基准修正的条件下，系统测量的高程漂移误差不可避免地随时间增长。由式（23）和（25）可知，z轴加速度计的零偏、安装角偏差均会导致俯仰角测量常值偏差，进而带来高程测量的线性增加的误差。通过线性平差处理可消除测量误差的1阶项，进一步提升动态测量高程的精度。相对高程与测量高程的关系示意如图5所示，图中2个约束点为A和B，经平差处理后的高程为约束点间的相对高程。由图5可知：平差过程不改变测量结果的趋势，相对高程同样能够反应区段内的线路沉降变化。

在惯性基准测量中，通常用高通滤波器对高程进行滤波处理，得到高低不平顺^［17］。基于长波高低不平顺在理论上可对线路沉降区段进行大范围筛查。但轨道几何参数中的高低不平顺包含了轨道不平顺和线路设计值，例如某高速线路的200 m长波高低不平顺如图6所示，图中K74，K76和K78处均为竖曲线区段。由图6可知：波形中包含了竖曲线所带来的结构性不平顺，与沉降区段的波形结构相近，显然会干扰线路沉降区段的识别。

因此利用台账信息减去轨道不平顺中的设计线形，是线路沉降识别的必要步骤。由于综合检测列车的检测里程与现场里程存在一定偏差，可利用线路坡度台账与实测平台俯仰角间的相关性配准里程，减去设计值之后再根据式（25）和（26）计算高程与长波高低不平顺。为尽可能覆盖不同范围的沉降区段，设定长波滤波器的截止波长为1 000 m。

2.3　线路沉降检测流程

线路沉降检测流程可分为以下3个步骤。

步骤1：根据线路实际情况设定沉降管理阈值，利用1 000 m长波高低不平顺筛查线路中可能的沉降区段。

步骤2：对于某一沉降区段，选定能够包含沉降区段的约束点，利用轨道的相对高程对沉降区间的最大沉降值进行监测。

步骤3：利用静态测量手段对该区段复核测量，并制定维修计划。

3 实测数据分析

为验证本文提出的线路沉降识别与监测方法的可行性，某型号综合检测列车分别于2022年11月和2022年3月在环形铁道试验线和郑州试验线开展中低速和高速试验。试验中使用的惯性测量单元（IMU）包括光纤陀螺仪和石英挠性加速度计，陀螺仪零偏稳定性指标为0.05° · h^-1，加速度计零偏稳定性指标为50 μg。

3.1　惯性参数准确性

姿态测量是轨面高程测量的核心环节，因此首先对惯性解算参数进行分析。通常用轨道超高的重复性来评估姿态解算精度。在环形铁道试验线开展的试验中，检测列车以不同的速度进行了6次重复测量（test1—test6），轨道超高重复测量结果如图7所示。由图7可知：多次测量的轨道超高表现出较好的重复性。

从6组数据任选2组，分别采用传统方法和本文方法计算轨道超高重复性偏差，重复性指标见表1。表中：对2次测量重复性偏差（绝对值）做从小到大的排序，取95%处的数值为95%分位数；重复性偏差的平均值简称为均值。由表1可知：基于本文方法测量轨道超高重复性偏差的95%分位数最大值为0.9 mm，基于传统方法测量的95%分位数最大值为1.82 mm，说明了本文方法的测量重复性更优；本文方法和传统方法的重复性偏差均值分别为0.24和0.38 mm，重复性相比于传统方法提升了37%左右，说明了该方法对姿态测量准确性的提升较为明显。

3.2　线路沉降检测准确性

结合试验数据对本文提出的线路沉降检测方法进行重复性与准确性分析。试验高铁线路总长度约为140 km，速度等级为350 km · h^-1。综合检测列车在该线路进行了多次往返测量工作，对该组试验数据开展分析。以6 km为窗长，统计了K20—K152内1 000 m长波高低不平顺的重复性偏差95%分位数，如图8所示。由图8可知：1 000 m高低不平顺的重复性偏差95%分位数在3.5 mm以内，满足线路沉降检测的精度要求。

测量平台的俯仰角可以反映线路的坡度^［19］，利用线路纵断面台账可生成轨道高低不平顺和坡度的设计值。先用坡度信息修正测量数据的里程，再减去高低不平顺的设计值，得到1 000 m截止波长高低不平顺偏差。减去设计值后，2次检测（test7和test8）的1 000 m高低不平顺和偏差如图9所示，统计重复性指标见表2。由图9和表2可知：减去设计值后，1 000 m高低不平顺偏差的95%分位数为1.57 mm，均值为0.58 mm，最大值为4.4 mm，重复性指标相较于图8更优，说明该处理方式有助于提升检测结果的重复性；减去设计值后，高低不平顺偏差更能直观反映局部沉降区段，基本不受到图6中所展示的线路竖曲线的干扰。

综合检测列车位于试验线路开展了多次重复检测试验，线路速度等级为120 km · h^-1，检测速度覆盖60~110 km · h^-1。设定该速度等级的线路以—120 mm作为筛查沉降区段的阈值，并对风险沉降区段进行重复测量，结果如图10所示。图10（a）为全线的1 000 m高低不平顺，图10（b）为沉降区段的1 000 m高低不平顺，图10（c）为沉降区段500 m范围内的相对高程。

由图10可知：风险沉降区段的里程在K2—K3范围内；设定两侧约束点的距离为500 m，测量K2+100—K2+600范围内的相对高程，峰值处表现出较好的重复性。

统计图10中沉降区段多次检测的最大沉降，见表3。由表3可知：K2+253附近4次检测1 000 m高低不平顺的最大重复性偏差为0.5 mm，相对高程的最大重复性偏差为5 mm；当行车速度为65 km · h^-1时，系统测量的最大沉降与更高速下的偏差较大，说明系统在低速下的测量准确性有所下降，但仍在可接受范围内。以上数据说明了本文的方法可对线路局部沉降进行识别和监测，且能够适用于不同速度等级的应用场景。

试验线某区段布设了轨道控制网（CPⅢ），使用安伯格惯导型轨检仪配合轨道控制网可测量轨道的高程信息，仪器的标称精度为1 mm。下面以绝对测量数据作为参考值，对动态测量的高程进行准确性评估。列车的检测速度为80 km · h^-1，试验区段的长度为400 m，设计坡度为0.5‰，以起始和终止点作为约束点，计算系统测量和绝对测量的相对高程，2者偏差如图11所示。图11（a）为左轨相对高程，图11（b）为系统测量和绝对测量的偏差。计算相对高程偏差的各项指标见表4。

由图11和表4可知：在100~300 m区段存在1处较为明显的沉降区段，最大沉降值近35 mm；系统测量的相对高程与绝对测量的相对高程趋势接近，短波不平顺一致性较好；在80 km · h^-1的行车速度下，相对高程偏差的均值为0.54 mm，95%分位数为1.39 mm，最大值为1.87 mm，能够满足一般沉降变形的精度要求。根据惯性基准测量原理，行车速度越快，单位距离内的累积误差越小，系统能够实现准确测量的范围越大，因此系统在中低速下的精度保持能力更具代表性。

综上，本文所提出的方法测量的相对高程能够反映轨道在竖向的真实形变，能够满足一定距离范围内的线路沉降变形区段的车载快速测量的需求，具有广泛的应用前景。

4 结论

（1）基于本文方法测量轨道超高重复性偏差的95%分位数的最大值为0.9 mm，均值为0.24 mm，重复性偏差均值相比于传统方法提升约37%。本文所提出的方法对姿态测量准确性的提升较为明显。

（2）在350 km · h^-1速度等级线路，1 000 m高低不平顺偏差的95%分位数为1.57 mm，均值为0.58 mm。高低不平顺偏差能直观反映局部沉降区段，基本不受到线路竖曲线的干扰。

（3）在120 km · h^-1速度等级线路，风险沉降区段的1 000 m高低不平顺最大重复性偏差为0.5 mm，相对高程的最大重复性偏差为5 mm。本文的方法可对线路局部沉降进行识别和测量，且能够适用于不同速度等级的应用场景。

（4）在80 km · h^-1的检测速度下，系统测量的400 m范围内相对高程与绝对测量结果变化趋势接近，相对高程偏差的95%分位数为1.39 mm，最大偏差为1.87 mm。相对高程能够反映轨道在竖向的真实形变，能够满足线路沉降区段的车载快速测量的需求。

参考文献

[1] 易强，王继军，刘伟斌，等.高速铁路无砟轨道结构基础变形分级评价［J］.中国铁道科学，2022，43（6）：20-28.

[2] YI Qiang， WANG Jijun， LIU Weibin， et al. Classification Evaluation for Foundation Deformation of Ballastless Track Structure in High-Speed Railway ［J］. China Railway Science， 2022， 43 （6）： 20-28. in Chinese

[3] 罗林，张格明，吴旺青，等.轮轨系统轨道平顺状态的控制［M］北京：中国铁道出版社，2006.

[4] LUO Lin， ZHANG Geming， WU Wangqing， et al. Control of Track Smooth State of Wheel-Rail System ［M］. Beijing： China Railway Press， 2006. in Chinese

[5] 赵国堂.高速铁路无砟轨道—路基变形计算模型的研究［J］.中国铁道科学，2016，37（4）：1-8.

[6] ZHAO Guotang. Calculation Model of Ballastless Track and Subgrade Deformation on High Speed Railway ［J］. China Railway Science， 2016， 37 （4）： 1-8. in Chinese

[7] 王平，肖杰灵.高速铁路轨道平顺性检测关键理论与技术［M］.北京：上海科学技术出版社，2019.

[8] WANG Ping， XIAO Jielin. Key Theories and Technologies of Track Irregularities Detection for High-Speed Railways ［M］. Beijing： China Railway Press， 2006. in Chinese

[9] 张玉芝，杜彦良，孙宝臣，等.基于液力测量的高速铁路无砟轨道路基沉降变形监测方法［J］.北京交通大学学报，2013，37（1）：80-84.

[10] ZHANG Yuzhi， DU Yanliang， SUN Baochen， et al. High-Speed Railway Ballastless Track Roadbed Settlement Monitoring Method Based on Fluid Pressure Differential ［J］. Journal of Beijing Jiaotong University， 2013， 37 （1）： 80-84. in Chinese

[11] 杨友涛，孔延花，孔书祥.高速铁路轨检仪绝对测量方法应用研究［J］.铁道建筑，2010，12（3）：97-99.

[12] YANG Youtao， KONG Yanhua， KONG Shuxiang. Study on Application of Absolute Measuring Method with Rail Detector to High-Speed Railway Track Inspection Instrument ［J］. Railway Engineering， 2010， 12 （3）： 97-99. in Chinese

[13] 魏晖，朱洪涛，万坚.绝对测量在无砟轨道的轨向控制中的精度分析［J］.铁道工程学报，2012，29（5）：1-5.

[14] WEI Hui， ZHU Hongtao， WAN Jian. Analysis of Accuracy of Absolute Survey in Control of Horizontal Regularity of Ballastless Track ［J］. Journal of Railway Engineering Society， 2012， 29 （5）： 1-5. in Chinese

[15] CHEN Q J， NIU X J， ZHANG Q， et al. Railway Track Irregularity Measuring by GNSS/INS Integration ［J］. Navigation， 2015， 62 （1）： 83-93.

[16] 何秀凤，高壮，肖儒雅，等.InSAR与北斗/GNSS综合方法监测地表形变研究现状与展望［J］.测绘学报，2022，51（7）：1338-1355.

[17] HE Xiufeng， GAO Zhuang， XIAO Ruya， et al. Application and Prospect of the Integration of InSAR and BDS/GNSS for Land Surface Deformation Monitoring ［J］. Acta Geodaetica et Cartographica Sinica， 2022， 51 （7）： 1338-1355. in Chinese

[18] 陈仕明，魏世斌，李颖，等.基于互补滤波的轨道不平顺动态测量方法［J］.中国铁道科学，2022，43（1）：52-62.

[19] CHEN Shiming， WEI Shibin， LI Ying， et al. Dynamic Measurement Method of Track Irregularity Based on Complementary Filter ［J］. China Railway Science， 2022， 43 （1）： 52-62. in Chinese

[20] 陈起金.A-INS组合导航的铁路轨道几何状态精密测量技术研究［D］.武汉：武汉大学，2016.

[21] CHEN Qijin. Research on Precise Measurement Technology of Railway Track Geometric State Based on A-INS Integrated Navigation ［D］. Wuhan： Wuhan University， 2016. in Chinese

[22] 牛小骥，旷俭，陈起金.采用MEMS惯导的小口径管道内检测定位方案可行性研究［J］.传感技术学报，2016，29（1）：40-44.

[23] NIU Xiaoji， KUANG Jian， CHEN Qijin. Study on the Posibility of the PIG Positioning Using MEMS-Based IMU ［J］. Chinese Journal of Sensors and Actuators， 2016， 29 （1）： 40-44. in Chinese

[24] 江一夫，李四海，谢波，等.采煤机三维轨迹的惯性基厘米级相对测量方法［J］.中国惯性技术学报，2021，29（2）：178-183.

[25] JIANG Yifu， LI Sihai， XIE Bo， et al. Inertial-Based Centimeter-Order Relative Measurement Method of Shearer 3-Dimensional Path ［J］. Journal of Chinese Inertial Technology， 2021， 29 （2）： 178-183. in Chinese

[26] 黎奇，白征东，陈波波，等.GNSS/INS多传感器组合高速铁路轨道测量系统［J］.测绘学报，2020，49（5）：569-579.

[27] LI Qi， BAI Zhengdong， CHEN Bobo， et al. A Novel Track Measurement System Based on GNSS/INS and Multisensor for High-Speed Railway ［J］. Acta Geodaetica et Cartographica Sinica， 2020， 49 （5）： 569-579. in Chinese

[28] 翟婉明，赵春发，夏禾，等.高速铁路基础结构动态性能演变及服役安全的基础科学问题［J］.中国科学 · 技术科学，2014，44（7）：645-660.

[29] ZHAI Wanming， ZHAO Chunfa， XIA He， et al. Basic Scientific Issues on Dynamic Performance Evolution of the High-Speed Railway Infrastructure and Its Service Safety ［J］. Scientia Sinica： Technologica， 2014， 44： 645-660. in Chinese

[30] 陈仕明，侯智雄，王昊，等.轨道几何动态惯性测量误差评估方法［J］.铁道建筑，2022，62（5）：36-41.

[31] CHEN Shiming， HOU Zhixiong， WANG Hao， et al. Evaluation Method of Dynamic Inertial Measurement Error of Track Geometry ［J］. Railway Engineering， 2022， 62 （5）： 36-41. in Chinese

[32] 陈仕明.基于视觉惯性融合的轨道长波不平顺动态检测系统研究及应用［D］.北京：中国铁道科学研究院，2022.

[33] CHEN Shiming. Research and Application of Dynamic Inspection System for Track Longwave Irregularities Based on Vision/Inertial Fusion ［D］. Beijing： China Academy of Railway Sciences， 2022. in Chinese

[34] DISSANAYAKE G， SUKKARIEH S， NEBOT E， et al. The Aiding of a Low-Cost Strapdown Inertial Measurement Unit Using Vehicle Model Constraints for Land Vehicle Applications ［J］. IEEE Transactions on Robotics and Automation， 2001， 17 （5）： 731-747.

[35] 王琰，陈仕明，魏世斌，等.高速铁路线路轨面坡度及竖曲线半径动态测量方法及应用研究［J］.铁道学报，2024，46（3）：176-183.

[36] WANG Yan， CHEN Shiming， WEI Shibin， et al. Study on Dynamic Measurement Method and Application of Rail Surface Slope and Vertical Curve Radius for High Speed Railways ［J］. Journal of the China Railway Society， 2024， 46 （3）： 176-183. in Chinese

基金资助

国家自然科学基金资助项目(52278465)

中国国家铁路集团有限公司科技研究开发计划课题(K2022T002)

中国铁道科学研究院集团有限公司院基金课题(2021YJ218)