数字轨道地图平面线形特征提取方法研究

陶维杰， 蔡伯根， 王 剑， 刘 江， 上官伟

(1.北京交通大学电子信息工程学院，北京 100044；2.北京市轨道交通电磁兼容与卫星导航工程技术研究中心，北京 100044)

随着卫星导航系统及卫星定位技术的迅速发展，其在水路、航路和道路交通领域的应用也趋于成熟，全球各地对其在铁路交通领域的应用也展开了广泛的研究。相对于传统的轨道电路、应答器的定位方式，卫星定位有利于列车定位由依赖地面设备的被动定位转变为车载设备自主定位，这能极大地减少地面设备的铺设和安装，节省建设和维护成本。此外，列车自主定位有助于推进固定闭塞控制向移动闭塞控制的升级，从而充分增强线路的通过能力，提高列车的运行效率。因此，基于卫星定位的列车定位是下一代列车运行控制系统的重要发展方向，也是列控系统向智能化、自主化发展的必然趋势。列车运行在既定的轨道线路上，线路的地理和几何信息是数字轨道地图的关键组成部分，也是卫星定位应用于列车定位的关键基础，它是将列车的空间位置坐标匹配对应到轨道一维坐标的重要依据。

数字轨道地图数据的主要来源是轨道线路的测量，目前主要采用实时动态差分RTK(Real Time Kinematic)技术对轨道中心线进行测量，获得线路上一系列的离散点。如何利用测量数据高效准确地表示线路特征，并满足存储和地图匹配的应用要求，是数字轨道地图生成的重要问题。文献[1]提出用分段折线近似表示曲线轨道，并将横向误差作为约束条件，采用启发式算法对测量数据进行约简以尽可能简单高效地表示轨道。文献[2]提出数字轨道地图的分层结构，对测量数据进行剔除、约简后选取特征点在低尺度上描述线路，并通过插值在高尺度上细化数字轨道地图。文献[3]根据主曲线、优化理论和二分法等思想提出了三种用数据点和线段表示线路的方法，其中自适应半径算法综合性能较优。上述方法都采用一系列离散的点及其构成的折线描述线路，虽然能在一定程度上描述线路情况，但无法准确地捕捉到线路的几何线形特征，且地图的精度及地图匹配的效率会受制于地图数据点的个数。文献[4]将里程作为参数，采用三次多项式对线路平面的两个方向分别进行动态分段拟合，该算法能保证线路几何的连续与光滑，但其与实际线路的线形不完全一致，且对地图匹配算法的应用具有较高的计算要求。考虑到铁路线路在设计时就采用三种基本线形：直线、圆曲线和缓和曲线，因此可以从这三种线形的基本特点出发，对测量数据进行相应的特征识别及线形分段，进而对不同的线形元素分别进行拟合，尽可能遵循设计的思想完成数字轨道地图的平面几何线形的准确描述。

针对线形元素识别方面，主要有方位角和曲率识别两种方法，文献[5]采用启发式算法对里程-方位角图中的直线、圆曲线和缓和曲线段分别用横直线、二次多项式和斜直线拟合，并在特征点处一阶导数连续作为约束条件，通过不断调整特征点的位置寻找最优解，此方法计算复杂，且需事先选好初始特征点。文献[6]首先用三次样条曲线拟合测点，再对拟合后的曲率进行平滑和截断处理，然后根据曲率识别线形及特征点，该方法计算步骤较多且曲率截断处理需根据实际情况调整参数。文献[7-8]基于方位角计算曲率，根据曲率值分辨直线段和圆曲线段，根据曲率变化率确定缓和曲线段，通过曲率的直方图确定圆曲线的半径，但该方法忽略了线路中存在相同半径圆曲线段的情况，且结果跟直方图选取的区间大小相关。在线形拟合方面，主要有样条曲线法[6，9]、分段多项式曲线法[10]、最小二乘法[7-8]，这些方法都比较成熟且能精确描述线形。

综上所述，传统的用离散点及其构成的折线进行地图描述的方法简单易实现，但无法保证线形的连续与平滑；采用线形描述的方法能准确描述地图，但存在计算步骤过于繁琐、参数调整过多、应用场景受限的问题，因此本文根据数字轨道地图的平面线形描述的需求，提出一种平面线形的特征提取方法。考虑到先拟合再计算测点曲率的步骤复杂，且计算结果易因测点误差造成曲率的突变，而方位角变化相对平缓，算法先通过计算方位角再近似估计曲率，而后通过设定曲率阈值对于线形进行初步识别，根据横向误差的要求对直线段和圆曲线段进行迭代直至特征点确定，最后采用整体最小二乘法分别对直线、圆曲线和缓和曲线拟合得到相应的线形参数，尽可能按照线路的线形特征描述地图。

1 平面线形特征识别方法

因为受地形、地质、技术等因素的限制，铁路线路无法用一条长直线延续始终，其方向需要进行改变。出于行车安全考虑，在转向处需要将相邻的两段直线用曲线连接起来，这种曲线称为平面曲线。平面曲线根据性质分为两种：圆曲线和缓和曲线。圆曲线是具有一定半径的圆弧，其半径大小根据车速的不同有相应的设计规范；为实现曲率半径的逐渐变化，减少对车辆的冲击，直线和圆曲线之间采用缓和曲线平滑过渡，其半径由直线半径(无穷大)逐渐减小为圆曲线半径，设计中常采用回旋曲线。线路的平面线形基本组合为“直线-前缓和曲线-圆曲线-后缓和曲线”。

方位角和曲率是直观描述平面线形的两个量。不同线形元素的方位角存在以下特点：直线段的方位角为一定值，圆曲线的方位角根据转弯方向线性递增或递减，缓和曲线的方位角呈二次抛物线变化。而曲率是切线方位角的一阶导数，因此直线段上曲率恒为零、圆曲线段的曲率是非零常数1/R(半径的倒数)，缓和曲线段上某点的曲率则与其长度成正比，即呈线性变化，见图1。

线路中心线各测点的方位角是从该点正北方向线算起，按照顺时针方向至该点方向线之间的水平夹角，范围为0°～360°。本文采用圆弧切线法计算测点的方位角，首先利用“三点定圆法”计算圆心和半径，然后将中间点在圆弧的切线与正北方向的夹角作为该点的方位角，见图2，公式为

( 1 )

式中：αi为i点的方位角；(xp，yp)为i+1点在切线方向上的投影点坐标；(xi，yi)为i点坐标。

图2 圆弧切线法计算方位角

测点的曲率可以通过方位角的一阶导数近似进行计算。曲率为

( 2 )

式中：ρi为i点的曲率；αi、αi+1分别为i、i+1点的方位角；dis(i，i+1)为i，i+1两点间的距离。

根据直线段曲率恒为零、曲线段曲率非零的特点，选取相应的曲率阈值即可将测点初步分为直线段和曲线段。此外，本文引入横向误差作为约束条件以确保直线段分类的准确性，并进一步确定圆曲线和缓和曲线段的范围。

横向误差是测点到相应拟合线形的正交投影距离，对应数字轨道地图在垂直股道方向上的误差。铁路上平行股道的线路中心线间距为5 m左右，若以中线2.5 m作为列车股道占用识别的分界线，并且作为参考的数字轨道地图精度应高一个数量级，则数字轨道地图本身的横向误差不应超过0.25 m。

算法的流程如下：首先根据轨道线路中心的测点坐标计算方位角，由方位角计算曲率，考虑到测点误差对计算的影响，本文采用滑动平均算法对方位角和曲率进行平滑滤波；根据直线段曲率为零的特点，设定阈值将各测点初步分类为直线段和曲线段；以横向误差为约束条件，采用相应的方法依次对直线段和曲线段进行迭代拟合，确定直线、圆曲线和缓和曲线的起始点、终止点和相关参数，最终得到线路整体的平面线形。算法的实现流程见图3。

图3 平面线形拟合算法流程

2 平面线形拟合方法

确定测点的平面线形之后，需对不同的线形元素进行拟合，求解线形的相关参数。通常根据线形的特征选取相应的数学模型，将测点的坐标作为观测量，线形的参数作为待估计量，一般采用最小二乘法求解。

经典的最小二乘法都以y为因变量，以x为自变量，假设因变量有误差，自变量无误差。这会导致选取的自变量和因变量不同时，得到的拟合结果也不同。在实际线路测量中，由于仪器、模型等存在误差，各测点的横纵坐标都不可避免地存在误差，整体最小二乘法就是在同时考虑x、y误差的情况下求解最优的参数估计，其数学模型称为EIV(Errors in Variables)模型。

EIV的线性数学模型为[11]

Y+eY=(A+eA)·ξ

( 3 )

式中：Y为n维观测向量；eY为Y的随机误差向量；A为n×m维系数矩阵；eA为系数矩阵A的随机误差矩阵；ξ为待求的m×1维参数向量。

当观测向量和系数矩阵为等精度观测时，eY和eA独立且都服从零均值定方差的高斯分布。

将式( 3 )变形为

( 4 )

式中：M=[AY]，e=[eAeY]。

整体最小二乘的约束准则为

min‖eA；eY‖F

( 5 )

式中：‖P‖F为n×m维矩阵P的Frobenius范数，其定义如下

( 6 )

式中：pij为矩阵P第i行第j列的元素；tr(·)为矩阵的迹。

常用的整体最小二乘数值求解方法有奇异值分解法SVD(Singular Value Decomposition)、完全正交法、增广矩阵法、拉格朗日迭代法等[12]。

2.1 直线线形拟合

设直线方程为yi=axi+b，i=1，2，3，…，n，(xi，yi)为n对线路中线的测点坐标，a为斜率，b为截距。考虑到x和y都分别存在误差vx和vy，方程则表示为

yi+vyi=a(x+vxi)+bi=1，2，3，…，n

( 7 )

其误差方程按照EIV模型可以表示为

(A+eA)·Δξ=Y+eY

( 8 )

注意到系数矩阵A中，存在元素固定为1的一列，此问题可以转化成混合最小二乘问题，根据混合最小二乘法求解计算[13]，此处不再展开。

2.2 圆曲线拟合

假设圆的圆心为(xc，yc)，半径为r，测点(xi，yi)到圆弧上的距离di(等于测点与圆心之间的距离减去半径)为

( 9 )

(10)

其误差方程为

(11)

其矩阵表达式为

Aδξ=Y+eY

(12)

式中：δξ为待求参数向量的修正值。且

未知参数向量修正值的解为

δξ=(ATA)-1ATY

(13)

(14)

(15)

2.3 缓和曲线拟合

线路设计中缓和曲线常用回旋线，回旋线的基本公式为rl=B2，B为回旋线的参数。回旋线的终点处，r=R(R为圆曲线半径)，l=Ls(Ls为回旋线总长度)，故RLs=B2。计算时需要建立局部坐标系，将直缓点定为原点，见图4。

图4 回旋曲线

回旋曲线的计算式为

(24)

实际应用时，常用三次抛物线来近似计算

(25)

式中：C=B2=RLs。

局部坐标系的计算结果经过坐标变换后，即可得到整体坐标系下的表达式为

式中：xZH、yZH为直缓点在整体坐标系下的绝对坐标；xlocal、ylocal为缓和曲线上的点在局部坐标系中的坐标；β为两个坐标系之间的旋转角。

在确定完直线段和圆曲线段的范围和相关参数之后，对缓和曲线段建立局部坐标系，按照三次抛物线对其进行拟合，同样采用整体最小二乘法，此处不再赘述。

3 算法测试

3.1 数据描述

本文选取青藏铁路一段线路中心线的GPS测量数据进行算法的测试与验证，数据采用Navcom SF-2050双频差分GPS接收机采集，精度为厘米级。该段线路全长约14.7 km，原始数据经预处理后得到6 100个测点，点间距最大3 m，最小1.5 m，平均2.4 m。

首先将测点经纬度转化为平面坐标，本文利用通用横轴墨卡托UTM(Universal Transverse Mercator)投影，为方便计算将起点转化为坐标原点，该段线路平面线形见图5。

图5 线路平面线形

3.2 线形识别与拟合结果

采用文中的“圆弧切线法”计算各测点的方位角，见图6，由方位角可以初步判断出该段线路包含5个曲线段，6个直线段。局部放大后，可以观察到方位角受测量及计算误差的影响存在噪声，为减少噪声对曲率计算的干扰，采用滑动平均法滤波。

用滑动平均后的方位角进行曲率的计算，同样地，对计算后的曲率进行滤波，见图7。从图7中可以观察到滤波后的曲率噪声明显减少。

(a)整体图

(b)局部放大图图6 测点的方位角变化

图7 测点的曲率变化

图8 直线段和曲线段阈值及初始分段

考虑到线路测量数据的精度达到了厘米级，因此本文设定横向误差的最大值为0.1 m，并以此作为限制条件依次对直线段、曲线段进行直线、圆曲线拟合，对剩余的缓和曲线段进行三次抛物线拟合，最终得到线路的整体线形及相应的参数。

拟合后线路整体的曲率图和平面线形的分布图分别见图9、图10。

3.3 算法的性能指标

横向误差Eacross：所有测点到相应拟合线形的正交投影距离，包括最大值和平均值，其值表征了拟合的平面线形精确度及准确度，越小表明拟合线形在垂直股道方向上偏差越小。

图9 拟合后曲率

图10 拟合后线路的整体平面线形

里程误差Ealong：所有测点的原始累积里程值与拟合后各测点在纵向上的累积里程值的对比，其值表明了拟合线路在纵向上的长度损失。

本次算法的性能指标结果见表1。

表1 算法性能指标结果

由以上结果统计，总结出以下结论：

(1)算法可以有效提取出线路的三种线形并对其进行拟合，通过极少数关键点及相应的参数对整段线路的平面线形进行表达，本文可以用22个线形的分界点及线形的参数表达出6 100个测点的信息，数据的约简率极高。

(2)平面线形的拟合结果表明算法的精度很高，在0.1 m的横向误差约束下，最大的横向误差不超过0.083 m，满足列车定位及股道识别的精度要求；且算法的里程累积误差最大仅为0.015 m，可以忽略不计。

(3)本算法拟合出来的线形不仅可以用来计算线路的坐标，还能准确描述出线路的方位角及曲率变化，丰富了数字轨道地图的内容。

4 结束语

高精度的数字轨道地图是GNSS铁路应用的重要基础，本文采用基于方位角的曲率方法对线路的平面线形进行特征识别与分段拟合，通过对青藏线实测线路数据的测试，验证了算法的可行性及高精度，并且相对于传统的一系列离散点的地图表示方法，本算法可以重构出线路的方位角及曲率变化情况，使数字轨道地图的平面线形更加准确完善，满足应用需求。