数字轨道地图线路特征分段拟合方法研究

袁 祎，陈光武

（1.兰州交通大学 光电技术与智能控制教育部重点实验室，兰州 730070；2.兰州交通大学 自动化与电气工程学院，兰州 730070；3.甘肃省高原交通信息工程及控制重点实验室，兰州 730070）

由我国北斗卫星导航、欧洲Galileo、美国全球定位系统（GPS，Global Positioning System）、俄罗斯全球导航卫星系统（GLONASS，Global Navigation Satellite System）组成的全球导航卫星系统（GNSS，Global Navigation Satellite System）在铁路交通领域得到了广泛的应用。基于GNSS 的列车定位系统定位方式不同于传统轨道电路及应答器辅助的定位方式，减少了轨旁设备的铺设，降低了建设和运营成本，提高了铁路运输的效率。基于卫星的列车定位方式作为铁路关键安全应用中的一部分，其对可靠性与安全性要求更高[1]。

考虑到铁路轨道可将列车约束至既定线路，许多学者将地图匹配技术作为定位辅助手段，以校正GNSS 定位误差，这种方式很大程度上提高了列车定位的精度，同时也节约了成本。而作为地图匹配方法实现的前提，轨道电子地图对线路地理和几何信息的准确描述直接影响了地图匹配算法的精度，因此，构建高精度的轨道电子地图是非常必要的[2]。在地理信息系统中常通过多条折线逼近的方式来描述曲线段，虽然有效，但会使线路的整体特征无法体现，同时损失了轨迹数据的精度[3]，不适用于铁路列车定位领域。

许多学者就轨道电子地图的构建方法进行了研究，文献[4—5]通过卡尔曼 (Kalman) 滤波算法构建预测方程，对采集到的数据集进行野值数据的剔除，利用角度差法对处理后的轨迹数据进行分段曲线拟合，完成轨道电子地图的构建；文献[6]针对原始GPS 地图测量数据，采用多点加权距离检验及Kalman 估计差值检验实现错误数据剔除，利用三次B 样条曲线反算生成拟合曲线，并按照一定的距离分辨率进行等距离插值，实现轨道电子地图的构建；文献[7]利用启发式算法对轨迹数据进行约简，通过分段折线对轨道线形进行描述；文献[8—9]研究从大量低精度的GPS 轨迹数据融合生成高精度GPS 数据的快速算法，通过对比得出最优的轨迹融合表示结果；文献[10]针对线形的识别进行研究，采用启发式算法对直线、圆曲线和缓和曲线段进行识别拟合，但需要选定初始点，且算法复杂；文献[11]提出了一种低成本的水平对齐拟合方法，通过曲率和方向讨论了点位置误差对道路线形表示噪声的影响，但拟合精度不高；文献[12]通过计算方位角近似估计曲率的方式对线形进行初步识别，通过设定横向误差对直线段和圆曲线段进行迭代，直至特征点确定，并采用最小二乘（LS，Least Squares）方法进行拟合，完成线路描述，但没有考虑到野值对拟合精度的影响。

本文考虑野值数据的影响，采用方位角近似估计曲率的方式对轨迹平面线形进行初步识别分类，设定横向误差为约束条件，对直线及曲线的分界点进行确定，采用引入自适应权值因子的LS 方法进行分段轨迹的拟合，以提高可靠轨迹数据拟合权值，而降低野值数据拟合权值的思路提高轨迹拟合平滑度，完成整体的线路地图构建。

1 铁路平面线形识别

1.1 铁路平面线形介绍

通常，实际铁路线路在设计的过程中受地形、技术等因素的限制，不能以单一长直线进行铁轨的铺设，为提高行车安全及舒适度，在列车转向处往往需要将相邻的两段直线用平面曲线连接起来，以减少对行驶列车的冲击。铁路线路通常划分为直线、圆曲线、缓和曲线3 种类型，其中，缓和曲线分为前缓和曲线与后缓和曲线，主要实现直线与圆曲线之间的过渡，铁路的线形平面示意如图1 所示。

图1 铁路线形平面示意

半径、方位角及曲率作为铁路平面线形的3 个描述指标，其特点如表1 所示，表中，R为圆曲线半径。

表1 铁路线形平面参数特点

1.2 基于方位角曲率特征的平面线形识别

铁路线路中心各轨迹点的方位角是指从该观测点正北方向起，绕顺时针方向至该观测点方向线之间所成的水平夹角。参照文献[12]，利用圆弧切线法计算各轨迹点的方位角，如图2 所示，采用“三点定圆法”计算圆心 (x0,y0)和 半径R，将观测点在圆弧的切线与正北方向的夹角作为该点的方位角，方位角表达式为

图2 圆弧切线法确定方位角示意

其中，θi为当前观测点的方位角，且θi2[0◦,360◦]；(xi,yi)为当前点i的坐标；(xp,yp)为后一点i+1在圆弧切线方向的投影坐标。则当前点的近似曲率为

其中，ρi为当前点i的近似曲率；dis(i,i+1)为前后两点坐标的投影距离。

依据平面线形的曲率特点，设定曲率阈值，对直线线形及曲线线形进行划分。将横向误差作为约束条件对圆曲线及缓和曲线进行划分，横向误差通常取0.25 m。

平面线形识别算法流程如图3 所示。（1）通过所测GNSS 轨迹点坐标计算方位角，由方位角计算出该点斜率；（2）出于测量误差会影响计算的考虑，使用滑动平均法对采集轨迹点的方位角和曲率进行平滑滤波；（3）利用直线路段曲率为零的特点，对轨迹进行直线段及曲线路段的初步划分；（4）设定横向测量误差为约束条件，对识别出的路段依据对应模型进行迭代拟合，从而确定出3 种线形的起始点、终止点及相关描述参数，得到轨迹的整体平面线形表示。

图3 平面线形识别算法流程

2 铁路平面线形拟合

2.1 最小二乘平面线形拟合方法

对轨迹点所处平面线形进行确定之后，需根据不同线形的特点选取对应的数学模型进行轨迹点的初步拟合，并求解不同线形模型的相关参数。本文使用常见的曲线拟合法—LS 方法[13]对轨迹点进行拟合。将轨迹点坐标作为观测量，线形模型参数为待估计量，即根据一系列离散点(xi,yi)，求近似曲线函数y=C(x)，将轨迹点到曲线拟合点的残差平方和最小的函数作为LS 拟合曲线函数。LS 形拟合算法的主要步骤如下。

Step1：根据给定的m个观测轨迹点数据(xi,yi)(i=1,2,···,m)，确定多项式曲线为

其中，aj(j=0,1,2,···,n) 为多项式各项的待定系数。

将式（3）用矩阵形式表示为

其中，A=[a0,a1,a2,···,an]为多项式的各项待定系数矩阵。

Step2：记点(xi,yi)到拟合曲线y=C(x)残差为 δi，故曲线拟合残差的平方和为

若要使其残差平方和最小，则曲线拟合的目标函数为

Step3：将曲线拟合问题转化为求解目标函数极小值点，则根据可得

k=0,1,···,n，整理可得

Step4：利用矩阵形式表示为XAT=Y，则可知

根据已知离散轨迹点集(xi,yi)，计算待定参数为

Step5：将式（9）代入式（4），即可获取所拟合曲线的多项式。

2.2 三种平面线形拟合

根据铁路平面线形的特点，可将拟合曲线多项式分为直线线形、圆曲线及缓和曲线3 种，其具体拟合方法如下。

2.2.1 直线线形拟合

设直线参数方程为y=ax+b，根据直线上的m个观测点坐标(xi,yi),i=1,2,···,m，根据LS 原理，可求得其系数为

将式（10）中的a、b代入y=ax+b，即可获得当前目标直线方程。

2.2.2 圆曲线拟合

设圆曲线参数方程为R2=(x−x0)2+(y−y0)2，其中， (x0,y0)为圆心坐标；R为圆曲线半径。根据LS 原理对曲线上的各坐标点 (xi,yi)进行拟合，残差平方和函数Q为

由LS 原理约束条件可知

联立式（11）和式（12），可求解相关参数为

2.2.3 缓和曲线拟合

按2.1 节求解方法，代入4 个轨迹点值即可确定上式a、b、c、d这4 个参数。

2.3 引入自适应权值因子的最小二乘拟合方法

在卫星数据采集过程中，受设备条件及环境因素影响，往往会存在野值数据。在标准的最小二乘拟合方法中，并未消除这种数据的影响，会使拟合线形不够平滑、精度下降。因此，本文引入自适应权值因子来降低野值数据的影响，提高可靠轨迹数据的拟合权值，使拟合线形更符合实际线路情况。

2.3.1 自适应权值因子的确定

若观测轨迹点在t时刻以前满足测试数据的噪声分布，而在t时刻出现较大偏差，则可将t时刻的观测数据视作是叠加在常规观测值上的一个野值，此时观测值为

其中，ω为与野值相关的偏差分量，yt为平稳状态下的理论观测值。此时残差估计值为

铁路平面线形拟合是已知整体轨迹数据的后处理过程，与轨迹观测点的时序无关，故定义当前测点与拟合曲线的相似耦合度因子 τ(t)为

其中，ρ(t)为两轨迹序列残差的相关系数为包含当前数据点的轨迹序列不含当前数据点的相邻轨迹序列；L为轨迹序列长度，依据采样周期进行相应调节。τ(t)反映了相邻序列间的相似程度，相关系数越高，则该观测点序列平滑，相关系数低，即判断当前点为野值数据。

对残差进行归一化为

其中，jδjmax为残差序列最大值的绝对值，jδjmin为残差序列最小值的绝对值。故权值因子 λt与残差δt的关系为

当观测点(xt,yt)与对应拟合曲线l的残差δt越大时，观测点与拟合轨迹相似耦合程度越低，权值因子 λt随之减小，降低观测点的可信度；反之，当残差 δt较小时，λt增大，提高观测点的可信度，因此，可通过λt对野值进行自适应修正。

2.3.2 改进最小二乘线形拟合算法流程

自适应权值因子改进后的最小二乘目标拟合函数为

自适应权值因子改进的最小二乘算法流程如下。

（1）按照标准最小二乘算法步骤拟合初始曲线；

（2）对观测轨迹点Ot的拟合残差序列δt(t=1,2,···,m)进行计算；

（3）根据式（17）、（18）求解相似耦合度因子 τ(t)及归一化系数η，并由式（19）计算得到自适应权值因子 λt；

（4）根据自适应权值因子 λt改进拟合参数，根据2.1 节拟合方程求解，获取轨迹分段线形。

3 试验验证及分析

3.1 试验数据准备

为验证本文拟合方法的有效性，采用PwrPak7-E1 组合导航系统，测得兰州交通大学至北环路的车辆GNSS 轨迹数据，并进行验证，线路全长约5.7 km，共计800 个测点，对数据进行经纬度坐标的UTM（Universal Transverse Mercator）平面投影转换，转换后的轨迹平面线形如图4 所示。

图4 车辆原始轨迹平面线形

3.2 线形识别及拟合结果

参照文献[12]，采用滑动平均滤波法对其轨迹各测点的方位角进行处理，如图5 所示，已知直线路段曲率恒为零，设定度/m）为直线及曲线的划分阈值，将直线路段属性设定为1，曲线路段属性设为0。可对车辆轨迹的分段进行初步划分，该车辆轨迹可分为17 段，包含9 个曲线段，8 个直线段，如图6 所示。根据2.2 节的拟合公式，分别对划分好的轨迹段进行拟合，可得到相应的线形及参数。拟合后的轨迹曲率特征如图7 所示，拟合后的平面线形的分布如图8 所示。

图5 车辆轨迹点曲率变化示意

图6 直线、曲线初始分段属性

图7 轨迹拟合后的曲率示意

图8 车辆轨迹线形分段拟合示意

拟合误差表示测量数据在拟合轨迹和真实轨迹中投影点的残差，其值代表了平面线形的拟合精度。本文将拟合误差绝对值的均值（AVG）、最大值（MAX）、方差（VAR）、均方误差（MSE）作为衡量指标，与传统最小二乘拟合方法进行对比，拟合误差相关参数如表2 所示。

表2 算法性能指标对比

由表2 可知，引入自适应权值因子改进的LS 拟合方法，在拟合精度上更具优势，可有效抑制野值数据对拟合平滑度的影响，使轨迹平面线形与真实线路更为接近。算法可对线路的直线、缓和曲线及圆曲线3 种平面线形进行有效识别并拟合，通过线形分界点及对应线形拟合参数对轨迹数据进行约简表示的方法适用于数字轨道地图的构建过程。

4 结束语

本文基于铁路线路平面线形的特点，采用方位角近似估计曲率的方法对轨迹平面线形进行初步识别分类；引入自适应权值因子改进的LS 方法进行分段轨迹的曲线拟合；并通过实测车辆轨迹数据进行验证。本文方法能够有效对3 种铁路平面线形进行识别，改进的LS 方法可降低野值数据拟合权值，从而提高轨迹拟合平滑度，比传统LS 方法拟合精度更高。