约束于圆弧的路线平面直线段重构算法

彭书航，蒲浩



约束于圆弧的路线平面直线段重构算法

彭书航1，蒲浩2, 3

(1. 西北工业大学 航空学院，陕西 西安 710072；2. 中南大学 土木工程学院，湖南 长沙 410075；3. 高速铁路建造技术国家工程实验室，湖南 长沙 410075)

提出一种利用一组(,)直角坐标点拟合与给定的圆弧相切的直线的算法。该算法遵循总体最小二乘准则来构建约束于圆弧的拟合直线的数学模型，在分析函数单调隔根区间的基础上，利用反拉格朗日插值迭代法计算出函数在定义域内的所有数值根，再经最小二乘准则检验获得数学模型的解，即拟合直线的回归参数。研究结果表明：拟合数学模型简明，基于单调隔根区间的迭代算法稳健，速度快，效果优，可用于重构路线平面的精确几何参数。

路线重构；约束直线拟合；总体最小二乘；隔根区间；反拉格朗日迭代

自动驾驶技术是人工智能研究发展的一个重要方向，其实现机理之一是利用导航定位技术来实时测量车辆的坐标，然后与路线平面的曲线数学方程进行匹配验证，从而决策车辆的运行方向和速度，控制车辆平稳、安全、舒适地运行。根据中华人民共和国国家统计局公布数据，截至2016年，我国道路通车里程达469.63万km，铁路通车里程达14.2万km。其中相当大一部分道路和铁路因建设于非信息化时代，缺乏准确的路线平面解析数据(曲线数学方程)，不能为自动驾驶提供基础的信息服务，因此需要通过现代测量手段来测量路线平面上的一系列坐标点，然后拟合出路线平面的曲线数学方程。Casal等[1−9]对路线中线的优化重构方法进行了广泛研究，它们的共同点是基于所有曲线、直线段均具有同权的拟合结果来构建拟合数学模型。当曲线较长而与其相邻的直线测点少且密时，测量误差对直线拟合效果的影响非常明显，导致其间不能插入合理的缓和曲线，此时宜考虑直线段与曲线段具有不同的拟合权，给予圆弧优先拟合权，直线拟合须与优先拟合的圆弧相切，确保车辆行驶轨迹的连续性。针对这一难点，本文提出一种解算约束于圆弧的直线拟合方法。

1 最小二乘直线拟合

平面直线的最小二乘拟合分为普通最小二乘拟合和总体最小二乘拟合。设直线的回归方程为：

则普通最小二乘拟合的准则是：

总体最小二乘的准则是：

式中：为直线的斜率；为轴的截距；d为测点到拟合直线的垂直距离。和是回归参数。

普通最小二乘计算过程简单，但只考虑了自变量的测量误差，其实质是使因变量坐标差的平方和最小，拟合结果往往不是最优的，直线斜率绝对值越大，拟合结果越差。总体最小二乘是使各测点到拟合直线垂直距离的平方和最小，此准则同时考虑了和变量的测量误差，可以获得直线的最佳拟合效果[10]。表1和表2为某企业铁路专用线上一个直线段RTK测量数据的2种拟合结果(为遵循相关法律，表中坐标数据经过变换处理)。计算结果表明，直线较陡时，总体最小二乘拟合明显优于普通最小二乘拟合。

因此，本文采用总体最小二乘准侧来拟合约束于圆弧的直线。

表1 RTK测量数据

表2 2种准则的直线拟合比较

2 拟合数学模型

车辆行驶中线是由若干直线段和圆弧段(可能有过渡缓和曲线)交互相连构成。当直线段短，测点少且密集，不能通过过渡曲线与圆曲线相连时，宜先拟合出圆弧的参数，再拟合与之相连的直线，使直线与圆弧相切。

刘元朋等[11]给出了引入系数约束条件的圆曲线总体最小二乘拟合方法，解决了其他迭代法运算量大和近似算法拟合效果差的问题。对于约束于圆弧的直线拟合，引入系数约束条件不是最佳方法。Francisco等[14]以航向角为测量对象，拟合结果是航向角最优而非测点到拟合直线的垂直距离平方和最小。

假设圆弧的圆心坐标为(，)，圆的半径为，那么圆的参数方程为：

若拟合直线与圆弧的切点为(00)，圆心到点的方向角为，则

过这一点的切线方程便可以写为：

经过整理，可得：

式(7)为拟合直线的回归方程，方向角为拟合直线回归参数。

于是第个样本点到直线的距离的平方为：

全部样本点的距离的平方和为：

式(9)中：除了切点方向角，其他所有量都是已知量，即2是的一元单值函数。

设函数()=2，则：

根据整体最小二乘准则，切点方向角应使()取得最小值，即()的一阶导数=0。

设定

于是：

用倍角公式将其展开，可以得到下式：

等号两侧同时平方并整理为sin的表达式：

式(13)换元变形：

令=4(2+2)，=4()，=2+2，24，=22，则：

3 反拉格朗日插值迭代

方程()=0一般是一个4次方程(≠0)，其求根公式极为复杂[12]，有时会因计算机的计算误差而导致求根结果不能满足精度要求，有时会出现虚根。

李惠[13]在给定的区间内采用二分或切线迭代法来计算一元四次方程的近似解。每次人工给定一个区间时，如果区间内有奇数个根，则迭代出其中一个根，当区间内有偶数个根时，则给定区间内被判定为无根(因为区间两端点的函数值乘积大于0)。因此该算法需要不断调整初始区间，效率很低，且不稳健。

可见，迭代法求根的关键是计算回归方程的单调的隔根区间。

求函数零点的迭代方法有很多种，比如牛顿迭代法、二分法等。牛顿法的收敛条件过于苛刻，不具有普遍性[15]；二分法收敛速度较慢，为此本文采用反拉格朗日插值迭代法。

如图1，某函数在隔根区间[0，2]上存在零点，于是选择端点(0，(0))，(2，(2))以及中点(1，(1))为初值。过这3个点可以得到一条开口向下的抛物线，这条抛物线的方程可以由拉格朗日插值公式给出。

图1 反拉格朗日插值方法

拉格朗日插值公式为：

这是一个关于的次多项式。当=2，在求解2()=0时，会出现多根，增加迭代计算难度。为了让问题尽可能的简单化，可以把看成的因变量，构造反拉格朗日插值公式：

这是一条开口朝左或者朝右的抛物线，它有且只有1个根，如图1中的曲线所示。

自此，可在[2(0),1]区间进行下一次迭代：分别以(2(0),(2(0)))和(1，(1))为端点，重复上述操作，直到lim(2(0))→0，或[2(0),1]区间距离小于阀值为止。在达到精度要求的前提下，2(0)便是方程的解。

反拉格朗日插值法并不是一种全局收敛的迭代法，只有在有单根的单调区间上才绝对收敛，证明如下：

使用反拉格朗日插值法解非线性方程的具体步骤如下：

1) 初始隔根区间[0,2]，且(0)*(2)< 0；

以表3和表4中的实测数据来检验(测点数目及坐标同表1)。

表3 10个样本点坐标

表4 拟定圆弧坐标

计算式(15)中的参数：=1.0，=0.010 102 766 20，=−0.000 221 859 0，=−0.000 329 474 4，= −0.000 002 278 0，精度要求为=0.000 000 001。

反拉格朗日插值法迭代8次就获得了区间[−1，0.028 625 458 667 695 315]的近似解，二分法则需要30次迭代才能满足精度要求，可见反拉格朗日迭代法速度快，精度高。

4 回归参数的数值解

()=A+D−≥0，(−)=A+D+2AD≥0

可见，()和(1)都是非负值，在[−1,1]区间内，()有0~4个解，分布在4个单调区间，如图2中的第1个曲线()。

为了求解()的单调区间，可先求出其一阶导数()的所有零点。而()是一个3次函数(式(18))，仍然需要用反拉格朗日插值法来计算数值根。为此可再次求导，计算2阶导函数()的零点。()是一个二次函数(式(19))，它的零点可以通过一元二次方程的求根公式直接计算。

求函数()的隔根区间及数值根的计算过程如下：

1)设()0的2个解分别为21和22，如果21不存在，则令21=−1.0；如果22不存在，则令22=21；

2) 构造()函数的3个单调区间：[−1,21]，[21,22]，[22,1]，如图2中的()曲线。

3) 单调区间内有无根的判别：求区间2个端点的函数值，如果它们的乘积为正，那么在这个区间上没有根；如果乘积为负，那么在这个区间上有1个根。于是逐区间计算()=0的3个根31，32和33：用拉格朗日反插值法计算增区间[−1,21]内的根31(如果无根，则31=−1)，[21,22]内的根32(如果无根，则32=31)，[22,1]内的根33(如果无根，则33=32)。

4) 构造()函数的4个单调区间：[−1,31]，[31,32]，[32,33]，[33,1]。

5) 计算各单调区间的根，如果有根，则加入()=0的解集1，否则忽略该单调区间。

根据式(14)，还需要对()=0所有的解求反正弦值。设sol是从1到集合S的一个映射，并且1与S的关系为：

反正弦函数在[0, 2π]上是多值函数，一个值对应2个值。在从式(12)到式(13)的平方变形过程中，会产生增根，因此S中的元素并不一定都是方程(11)的解。为了排除增根，还应该把S所有的元素带入方程(11)检验，当方程的值在误差限范围内时，可以被认为是方程的一个解。依次检验S的每一个元素，就可以得到方程(11)的解集0。

图2 函数的隔根区间与导数0点的依赖关系

Fig. 2 Mapping between the function’s interval containing a null point and the null point of the derivative

0中的元素可以使′()=0，但并不一定能使()取到最小值。因此，还要根据式(9)求出0中每一个元素的函数值()，并比较各函数值的大小。能使得()取最小值的解才是本文问题的最优解，即与给定圆弧相切的拟合直线的回归参数，拟合直线的斜率：

5 实例应用

为验证本文的算法，用Microsoft Visual C++ 2008编程进行测试。这个测试中输入数据为10个RTK测点，以及1个已经拟合的圆弧(见第3节 表2)。

运行中间成果数据即式(11)的4个系数：=−159 329.407 064 198，=−1 572 003.362 484 304，=159 400.781 899 597，=3 155 740.414 248 4(式(15)的5个系数值见第3节的检验)。

有效根即回归参数=6.276 137 776 1，拟合直线斜率=141.891 314 641 1，拟合残差的平方和为0.008 698 153 2，稍大于总体最小二乘拟合残差的平方和，这是因为切圆使拟合直线发生了偏转。

6 结论

1) 理论分析及实际验证表明，将测点少、没有缓和曲线的短小直线拟合到与圆弧相切，会使路线平面的重构结果更接近实际工况。

2) 普通最小二乘的拟合效果与测点的平面分布相关，当拟合直线的斜率绝对值越大时，拟合残差越大。总体最小二乘拟合与测点的平面分布无关，无论直线如何倾斜，拟合残差都相同，因此总体最小二乘的拟合效果最佳。

3) 遵循总体最小二乘准则来建立约束于圆弧的直线拟合数学模型，经换元变化，得到一个关于(sin)的一元四次方程，方程系数均是已知量的二次函数，因此计算过程简单，易于编程实现。

4) 根据函数一阶导数的0点是函数极值点原理，可以准确构建函数的隔根区间，从而降低迭代计算工作量，算法简明实用。

[1] Casal G, Santamarina D, Vázquez-Méndez M E. Optimization of horizontal alignment geometry in road design and reconstruction[J]. Transportation Research Part C Emerging Technologies, 2017, 74(1): 261−274

[2] 刘威, 胡光常, 麻丁一, 等. 铁路既有线平面自动重构方法研究[J]. 高速铁路技术, 2015, 6(6): 79−84. LIU Wei, HU Guangchang, MA Dingyi, et al. Research on existing railway plane line auto-reconstruction[J]. High Speed Railway Technology, 2015, 6(6): 79−84.

[3] 郭良浩, 刘成龙, 宋韬, 等. 铁路既有线平面和竖面线形精确分段方法研究[J]. 铁道工程学报, 2014, 31(7): 48−52. GUO Lianghao, LIU Chenglong, SONG Tao, et al. Research on the new method for accurate linear segmentation of plane and vertical curve type in existing railway[J]. Journal of Railway Engineering Society, 2014, 31(7): 48−52.

[4] 张航, 黄云, 龚良甫. 基于三次样条函数拟合公路平面线形方法研究[J]. 武汉理工大学学报(交通科学与工程版), 2007, 31(5): 925−927. ZHANG Hang, HUANG Yun, GONG Liangpu. Study on fitting highway alignment based on the cubic spline function[J]. Journal of Wuhan University of Technology (Transportation Science & Engineering), 2007, 31(5): 925−927.

[5] 刘苏, 王文强, 查旭东, 等. 基于法线偏差的旧路平面线形拟合精度评估方法[J]. 中国公路学报, 2007, 20(5): 36−40. LIU Su, WANG Wenqiang, ZHA Xudong, et al. Estimate method for accuracy of fitting plane linear in old highway based on normal error[J]. China Journal of Highway and Transport, 2007, 20(5): 36−40.

[6] LI Wei, PU Hao, Paul Schonfeld, et al. Mountain railway alignment optimization with bidirectional distance transform and genetic algorithm[J]. Computer-Aided Civil and Infrastructure Engineering, 2017, 32(8): 691− 709

[7] LI Wei, PU Hao, Paul Schonfeld, et al. Methodology for optimizing constrained 3-dimensional railway alignments in mountainous terrain[J]. Transportation Research Part C, 2016, 7(68): 549−565

[8] 李伟. 铁路数字选线关键技术研究与应用[D]. 长沙: 中南大学, 2014. LI Wei. Research and application for key technologies of digitalizing railway alignment design[D]. Changsha: Central South University, 2014.

[9] Robert T, Gunner L, Sarbaz O, et al. Using naturalistic field operational test data to identify horizontal curves[J]. Journal of Transportation Engineering, 2012, 138(9): 1151−1160.

[10] 丁克良, 沈云中, 欧吉坤. 整体最小二乘法直线拟合[J]. 辽宁工程技术大学学报(自然科学版), 2010, 29(1): 44−47. DING Keliang, SHEN Yunzhong, OU Jikun. Methods of line-fitting based on total least-square[J]. Journal of Liaoning Technical University (Natural Science), 2010, 29(1): 44−47.

[11] 刘元朋, 张定华, 桂元坤, 等. 用带约束的最小二乘拟合平面圆曲线[J]. 计算机辅助设计与图形学学报, 2001, 10(16): 1382−1385. LIU Yuanpeng, ZHANG Dinghua, GUI Yuankun, et al. Fitting planar circles with constrained least squares[J]. Journal of Computer-Aided Design & Computer Graphics, 2001, 10(16): 1382−1385.

[12] 吴佐慧, 廖军, 徐行忠, 等. 一元三次、四次方程根的行列式解法[J]. 湖北大学学报(自然科学版), 2014, 4(36): 381−388.WU Zuohui, LIAO Jun, XU Xingzhong, et al. A solution of the root of the cubic and quartic equationusing determinant[J]. Journal of Hubei University (Natural Science), 2014, 4(36): 381−388.

[13] 李惠. 用C程序求一元四次方程的近似解[J]. 中国科教创新导刊, 2012, 9(25): 98−100. LI Hui. Find the approximate solution of a quartic equation with one unknown in C program[J]. China Education Innovation Herald, 2012, 9(25): 98−100.

[14] Francisco J C, Ana M P, José M C, et al. Use of heading direction for recreating the horizontal alignment of an existing road[J]. Computer-Aided Civil and Infrastructure Engineering, 2015, 30(4): 282−299.

[15] 喻文健. 数值分析与算法[M]. 北京: 清华大学出版社, 2012: 41−45. YU Wenjian. Numerical analysis and algorithm[M]. Beijing: Tsinghua University Press, 2012: 41−45.

A recreating algorithm of straight segment of the horizontal alignment constrained by circular arc

PENG Shuhang1, PU Hao2, 3

(1. School of Areonautics, Northwestern Polytechnical University School of Areonautics, Xi’an 710072, China; 2. School of Civil Engineering, Central South University, Changsha 410075, China; 3. National Engineering Laboratory for High Speed Railway Construction, Changsha 410075, China)

On the basis of a set of (,) points, this paper proposes a method for fitting a regression line that is tangent to the given arc. The algorithm, based on Holistic Least Square Method, firstly established mathematical model fitting a line constrained by circular arc. On the basis of analysis of the intervals on which a function’s null point exists, the Inverse Lagrange’s Interpolation Method was used to calculate all of the function’s numerical solutions on the domain of definition. Additionally, the least squares criterion decided which of the numerical solutions can satisfy the mathematical model which had been established. Numerical solutions satisfying the mathematical model were exactly regression parameters of the regression line. Theoretical derivation and practical application indicate that this is a fast and stable algorithm with accurate results. This feature allows it to reconstruct precise geometric parameters of the Horizontal Alignment.

recreating the horizontal alignment; constrained line fitting; holistic least square method; intervals containing a null point; the inverse Lagrange’s interpolation method

10.19713/j.cnki.43−1423/u.2019.04.011

U412.3

A

1672 − 7029(2019)04 − 0915 − 07

2018−05−05

国家自然科学基金资助项目(51608543)；国家重点研发计划项目(2017YFB1201102)

蒲浩(1973−)，男，四川南充人，教授，博士，从事铁路线站数字化设计理论与方法研究；E−mail：haopu@csu.edu.cn

(编辑 涂鹏)