彭書(shū)航，蒲浩

?

約束于圓弧的路線平面直線段重構(gòu)算法

彭書(shū)航1，蒲浩2, 3

(1. 西北工業(yè)大學(xué) 航空學(xué)院，陜西 西安 710072；2. 中南大學(xué) 土木工程學(xué)院，湖南 長(zhǎng)沙 410075；3. 高速鐵路建造技術(shù)國(guó)家工程實(shí)驗(yàn)室，湖南 長(zhǎng)沙 410075)

提出一種利用一組(,)直角坐標(biāo)點(diǎn)擬合與給定的圓弧相切的直線的算法。該算法遵循總體最小二乘準(zhǔn)則來(lái)構(gòu)建約束于圓弧的擬合直線的數(shù)學(xué)模型，在分析函數(shù)單調(diào)隔根區(qū)間的基礎(chǔ)上，利用反拉格朗日插值迭代法計(jì)算出函數(shù)在定義域內(nèi)的所有數(shù)值根，再經(jīng)最小二乘準(zhǔn)則檢驗(yàn)獲得數(shù)學(xué)模型的解，即擬合直線的回歸參數(shù)。研究結(jié)果表明：擬合數(shù)學(xué)模型簡(jiǎn)明，基于單調(diào)隔根區(qū)間的迭代算法穩(wěn)健，速度快，效果優(yōu)，可用于重構(gòu)路線平面的精確幾何參數(shù)。

路線重構(gòu)；約束直線擬合；總體最小二乘；隔根區(qū)間；反拉格朗日迭代

自動(dòng)駕駛技術(shù)是人工智能研究發(fā)展的一個(gè)重要方向，其實(shí)現(xiàn)機(jī)理之一是利用導(dǎo)航定位技術(shù)來(lái)實(shí)時(shí)測(cè)量車(chē)輛的坐標(biāo)，然后與路線平面的曲線數(shù)學(xué)方程進(jìn)行匹配驗(yàn)證，從而決策車(chē)輛的運(yùn)行方向和速度，控制車(chē)輛平穩(wěn)、安全、舒適地運(yùn)行。根據(jù)中華人民共和國(guó)國(guó)家統(tǒng)計(jì)局公布數(shù)據(jù)，截至2016年，我國(guó)道路通車(chē)?yán)锍踢_(dá)469.63萬(wàn)km，鐵路通車(chē)?yán)锍踢_(dá)14.2萬(wàn)km。其中相當(dāng)大一部分道路和鐵路因建設(shè)于非信息化時(shí)代，缺乏準(zhǔn)確的路線平面解析數(shù)據(jù)(曲線數(shù)學(xué)方程)，不能為自動(dòng)駕駛提供基礎(chǔ)的信息服務(wù)，因此需要通過(guò)現(xiàn)代測(cè)量手段來(lái)測(cè)量路線平面上的一系列坐標(biāo)點(diǎn)，然后擬合出路線平面的曲線數(shù)學(xué)方程。Casal等[1?9]對(duì)路線中線的優(yōu)化重構(gòu)方法進(jìn)行了廣泛研究，它們的共同點(diǎn)是基于所有曲線、直線段均具有同權(quán)的擬合結(jié)果來(lái)構(gòu)建擬合數(shù)學(xué)模型。當(dāng)曲線較長(zhǎng)而與其相鄰的直線測(cè)點(diǎn)少且密時(shí)，測(cè)量誤差對(duì)直線擬合效果的影響非常明顯，導(dǎo)致其間不能插入合理的緩和曲線，此時(shí)宜考慮直線段與曲線段具有不同的擬合權(quán)，給予圓弧優(yōu)先擬合權(quán)，直線擬合須與優(yōu)先擬合的圓弧相切，確保車(chē)輛行駛軌跡的連續(xù)性。針對(duì)這一難點(diǎn)，本文提出一種解算約束于圓弧的直線擬合方法。

1 最小二乘直線擬合

平面直線的最小二乘擬合分為普通最小二乘擬合和總體最小二乘擬合。設(shè)直線的回歸方程為：

則普通最小二乘擬合的準(zhǔn)則是：

總體最小二乘的準(zhǔn)則是：

式中：為直線的斜率；為軸的截距；d為測(cè)點(diǎn)到擬合直線的垂直距離。和是回歸參數(shù)。

普通最小二乘計(jì)算過(guò)程簡(jiǎn)單，但只考慮了自變量的測(cè)量誤差，其實(shí)質(zhì)是使因變量坐標(biāo)差的平方和最小，擬合結(jié)果往往不是最優(yōu)的，直線斜率絕對(duì)值越大，擬合結(jié)果越差?？傮w最小二乘是使各測(cè)點(diǎn)到擬合直線垂直距離的平方和最小，此準(zhǔn)則同時(shí)考慮了和變量的測(cè)量誤差，可以獲得直線的最佳擬合效果[10]。表1和表2為某企業(yè)鐵路專(zhuān)用線上一個(gè)直線段RTK測(cè)量數(shù)據(jù)的2種擬合結(jié)果(為遵循相關(guān)法律，表中坐標(biāo)數(shù)據(jù)經(jīng)過(guò)變換處理)。計(jì)算結(jié)果表明，直線較陡時(shí)，總體最小二乘擬合明顯優(yōu)于普通最小二乘擬合。

因此，本文采用總體最小二乘準(zhǔn)側(cè)來(lái)擬合約束于圓弧的直線。

表1 RTK測(cè)量數(shù)據(jù)

表2 2種準(zhǔn)則的直線擬合比較

2 擬合數(shù)學(xué)模型

車(chē)輛行駛中線是由若干直線段和圓弧段(可能有過(guò)渡緩和曲線)交互相連構(gòu)成。當(dāng)直線段短，測(cè)點(diǎn)少且密集，不能通過(guò)過(guò)渡曲線與圓曲線相連時(shí)，宜先擬合出圓弧的參數(shù)，再擬合與之相連的直線，使直線與圓弧相切。

劉元朋等[11]給出了引入系數(shù)約束條件的圓曲線總體最小二乘擬合方法，解決了其他迭代法運(yùn)算量大和近似算法擬合效果差的問(wèn)題。對(duì)于約束于圓弧的直線擬合，引入系數(shù)約束條件不是最佳方法。Francisco等[14]以航向角為測(cè)量對(duì)象，擬合結(jié)果是航向角最優(yōu)而非測(cè)點(diǎn)到擬合直線的垂直距離平方和最小。

假設(shè)圓弧的圓心坐標(biāo)為(，)，圓的半徑為，那么圓的參數(shù)方程為：

若擬合直線與圓弧的切點(diǎn)為(00)，圓心到點(diǎn)的方向角為，則

過(guò)這一點(diǎn)的切線方程便可以寫(xiě)為：

經(jīng)過(guò)整理，可得：

式(7)為擬合直線的回歸方程，方向角為擬合直線回歸參數(shù)。

于是第個(gè)樣本點(diǎn)到直線的距離的平方為：

全部樣本點(diǎn)的距離的平方和為：

式(9)中：除了切點(diǎn)方向角，其他所有量都是已知量，即2是的一元單值函數(shù)。

設(shè)函數(shù)()=2，則：

根據(jù)整體最小二乘準(zhǔn)則，切點(diǎn)方向角應(yīng)使()取得最小值，即()的一階導(dǎo)數(shù)=0。

設(shè)定

于是：

用倍角公式將其展開(kāi)，可以得到下式：

等號(hào)兩側(cè)同時(shí)平方并整理為sin的表達(dá)式：

式(13)換元變形：

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

3 反拉格朗日插值迭代

方程()=0一般是一個(gè)4次方程(≠0)，其求根公式極為復(fù)雜[12]，有時(shí)會(huì)因計(jì)算機(jī)的計(jì)算誤差而導(dǎo)致求根結(jié)果不能滿足精度要求，有時(shí)會(huì)出現(xiàn)虛根。

李惠[13]在給定的區(qū)間內(nèi)采用二分或切線迭代法來(lái)計(jì)算一元四次方程的近似解。每次人工給定一個(gè)區(qū)間時(shí)，如果區(qū)間內(nèi)有奇數(shù)個(gè)根，則迭代出其中一個(gè)根，當(dāng)區(qū)間內(nèi)有偶數(shù)個(gè)根時(shí)，則給定區(qū)間內(nèi)被判定為無(wú)根(因?yàn)閰^(qū)間兩端點(diǎn)的函數(shù)值乘積大于0)。因此該算法需要不斷調(diào)整初始區(qū)間，效率很低，且不穩(wěn)健。

可見(jiàn)，迭代法求根的關(guān)鍵是計(jì)算回歸方程的單調(diào)的隔根區(qū)間。

求函數(shù)零點(diǎn)的迭代方法有很多種，比如牛頓迭代法、二分法等。牛頓法的收斂條件過(guò)于苛刻，不具有普遍性[15]；二分法收斂速度較慢，為此本文采用反拉格朗日插值迭代法。

如圖1，某函數(shù)在隔根區(qū)間[0，2]上存在零點(diǎn)，于是選擇端點(diǎn)(0，(0))，(2，(2))以及中點(diǎn)(1，(1))為初值。過(guò)這3個(gè)點(diǎn)可以得到一條開(kāi)口向下的拋物線，這條拋物線的方程可以由拉格朗日插值公式給出。

圖1 反拉格朗日插值方法

拉格朗日插值公式為：

這是一個(gè)關(guān)于的次多項(xiàng)式。當(dāng)=2，在求解2()=0時(shí)，會(huì)出現(xiàn)多根，增加迭代計(jì)算難度。為了讓問(wèn)題盡可能的簡(jiǎn)單化，可以把看成的因變量，構(gòu)造反拉格朗日插值公式：

這是一條開(kāi)口朝左或者朝右的拋物線，它有且只有1個(gè)根，如圖1中的曲線所示。

自此，可在[2(0),1]區(qū)間進(jìn)行下一次迭代：分別以(2(0),(2(0)))和(1，(1))為端點(diǎn)，重復(fù)上述操作，直到lim(2(0))→0，或[2(0),1]區(qū)間距離小于閥值為止。在達(dá)到精度要求的前提下，2(0)便是方程的解。

反拉格朗日插值法并不是一種全局收斂的迭代法，只有在有單根的單調(diào)區(qū)間上才絕對(duì)收斂，證明如下：

使用反拉格朗日插值法解非線性方程的具體步驟如下：

1) 初始隔根區(qū)間[0,2]，且(0)*(2)< 0；

以表3和表4中的實(shí)測(cè)數(shù)據(jù)來(lái)檢驗(yàn)(測(cè)點(diǎn)數(shù)目及坐標(biāo)同表1)。

表3 10個(gè)樣本點(diǎn)坐標(biāo)

表4 擬定圓弧坐標(biāo)

計(jì)算式(15)中的參數(shù)：=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次就獲得了區(qū)間[?1，0.028 625 458 667 695 315]的近似解，二分法則需要30次迭代才能滿足精度要求，可見(jiàn)反拉格朗日迭代法速度快，精度高。

4 回歸參數(shù)的數(shù)值解

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

可見(jiàn)，()和(1)都是非負(fù)值，在[?1,1]區(qū)間內(nèi)，()有0~4個(gè)解，分布在4個(gè)單調(diào)區(qū)間，如圖2中的第1個(gè)曲線()。

為了求解()的單調(diào)區(qū)間，可先求出其一階導(dǎo)數(shù)()的所有零點(diǎn)。而()是一個(gè)3次函數(shù)(式(18))，仍然需要用反拉格朗日插值法來(lái)計(jì)算數(shù)值根。為此可再次求導(dǎo)，計(jì)算2階導(dǎo)函數(shù)()的零點(diǎn)。()是一個(gè)二次函數(shù)(式(19))，它的零點(diǎn)可以通過(guò)一元二次方程的求根公式直接計(jì)算。

求函數(shù)()的隔根區(qū)間及數(shù)值根的計(jì)算過(guò)程如下：

1)設(shè)()0的2個(gè)解分別為21和22，如果21不存在，則令21=?1.0；如果22不存在，則令22=21；

2) 構(gòu)造()函數(shù)的3個(gè)單調(diào)區(qū)間：[?1,21]，[21,22]，[22,1]，如圖2中的()曲線。

3) 單調(diào)區(qū)間內(nèi)有無(wú)根的判別：求區(qū)間2個(gè)端點(diǎn)的函數(shù)值，如果它們的乘積為正，那么在這個(gè)區(qū)間上沒(méi)有根；如果乘積為負(fù)，那么在這個(gè)區(qū)間上有1個(gè)根。于是逐區(qū)間計(jì)算()=0的3個(gè)根31，32和33：用拉格朗日反插值法計(jì)算增區(qū)間[?1,21]內(nèi)的根31(如果無(wú)根，則31=?1)，[21,22]內(nèi)的根32(如果無(wú)根，則32=31)，[22,1]內(nèi)的根33(如果無(wú)根，則33=32)。

4) 構(gòu)造()函數(shù)的4個(gè)單調(diào)區(qū)間：[?1,31]，[31,32]，[32,33]，[33,1]。

5) 計(jì)算各單調(diào)區(qū)間的根，如果有根，則加入()=0的解集1，否則忽略該單調(diào)區(qū)間。

根據(jù)式(14)，還需要對(duì)()=0所有的解求反正弦值。設(shè)sol是從1到集合S的一個(gè)映射，并且1與S的關(guān)系為：

反正弦函數(shù)在[0, 2π]上是多值函數(shù)，一個(gè)值對(duì)應(yīng)2個(gè)值。在從式(12)到式(13)的平方變形過(guò)程中，會(huì)產(chǎn)生增根，因此S中的元素并不一定都是方程(11)的解。為了排除增根，還應(yīng)該把S所有的元素帶入方程(11)檢驗(yàn)，當(dāng)方程的值在誤差限范圍內(nèi)時(shí)，可以被認(rèn)為是方程的一個(gè)解。依次檢驗(yàn)S的每一個(gè)元素，就可以得到方程(11)的解集0。

圖2 函數(shù)的隔根區(qū)間與導(dǎo)數(shù)0點(diǎn)的依賴(lài)關(guān)系

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

0中的元素可以使′()=0，但并不一定能使()取到最小值。因此，還要根據(jù)式(9)求出0中每一個(gè)元素的函數(shù)值()，并比較各函數(shù)值的大小。能使得()取最小值的解才是本文問(wèn)題的最優(yōu)解，即與給定圓弧相切的擬合直線的回歸參數(shù)，擬合直線的斜率：

5 實(shí)例應(yīng)用

為驗(yàn)證本文的算法，用Microsoft Visual C++ 2008編程進(jìn)行測(cè)試。這個(gè)測(cè)試中輸入數(shù)據(jù)為10個(gè)RTK測(cè)點(diǎn)，以及1個(gè)已經(jīng)擬合的圓弧(見(jiàn)第3節(jié) 表2)。

運(yùn)行中間成果數(shù)據(jù)即式(11)的4個(gè)系數(shù)：=?159 329.407 064 198，=?1 572 003.362 484 304，=159 400.781 899 597，=3 155 740.414 248 4(式(15)的5個(gè)系數(shù)值見(jiàn)第3節(jié)的檢驗(yàn))。

有效根即回歸參數(shù)=6.276 137 776 1，擬合直線斜率=141.891 314 641 1，擬合殘差的平方和為0.008 698 153 2，稍大于總體最小二乘擬合殘差的平方和，這是因?yàn)榍袌A使擬合直線發(fā)生了偏轉(zhuǎn)。

6 結(jié)論

1) 理論分析及實(shí)際驗(yàn)證表明，將測(cè)點(diǎn)少、沒(méi)有緩和曲線的短小直線擬合到與圓弧相切，會(huì)使路線平面的重構(gòu)結(jié)果更接近實(shí)際工況。

2) 普通最小二乘的擬合效果與測(cè)點(diǎn)的平面分布相關(guān)，當(dāng)擬合直線的斜率絕對(duì)值越大時(shí)，擬合殘差越大。總體最小二乘擬合與測(cè)點(diǎn)的平面分布無(wú)關(guān)，無(wú)論直線如何傾斜，擬合殘差都相同，因此總體最小二乘的擬合效果最佳。

3) 遵循總體最小二乘準(zhǔn)則來(lái)建立約束于圓弧的直線擬合數(shù)學(xué)模型，經(jīng)換元變化，得到一個(gè)關(guān)于(sin)的一元四次方程，方程系數(shù)均是已知量的二次函數(shù)，因此計(jì)算過(guò)程簡(jiǎn)單，易于編程實(shí)現(xiàn)。

4) 根據(jù)函數(shù)一階導(dǎo)數(shù)的0點(diǎn)是函數(shù)極值點(diǎn)原理，可以準(zhǔn)確構(gòu)建函數(shù)的隔根區(qū)間，從而降低迭代計(jì)算工作量，算法簡(jiǎn)明實(shí)用。

[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] 劉威, 胡光常, 麻丁一, 等. 鐵路既有線平面自動(dòng)重構(gòu)方法研究[J]. 高速鐵路技術(shù), 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]. 鐵道工程學(xué)報(bào), 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] 張航, 黃云, 龔良甫. 基于三次樣條函數(shù)擬合公路平面線形方法研究[J]. 武漢理工大學(xué)學(xué)報(bào)(交通科學(xué)與工程版), 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] 劉蘇, 王文強(qiáng), 查旭東, 等. 基于法線偏差的舊路平面線形擬合精度評(píng)估方法[J]. 中國(guó)公路學(xué)報(bào), 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] 李偉. 鐵路數(shù)字選線關(guān)鍵技術(shù)研究與應(yīng)用[D]. 長(zhǎng)沙: 中南大學(xué), 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]. 遼寧工程技術(shù)大學(xué)學(xué)報(bào)(自然科學(xué)版), 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]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 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]. 湖北大學(xué)學(xué)報(bào)(自然科學(xué)版), 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]. 中國(guó)科教創(chuàng)新導(dǎo)刊, 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] 喻文健. 數(shù)值分析與算法[M]. 北京: 清華大學(xué)出版社, 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

國(guó)家自然科學(xué)基金資助項(xiàng)目(51608543)；國(guó)家重點(diǎn)研發(fā)計(jì)劃項(xiàng)目(2017YFB1201102)

蒲浩(1973?)，男，四川南充人，教授，博士，從事鐵路線站數(shù)字化設(shè)計(jì)理論與方法研究；E?mail：haopu@csu.edu.cn

(編輯 涂鵬)