數(shù)字軌道地圖平面線形特征提取方法研究

陶維杰， 蔡伯根， 王 劍， 劉 江， 上官偉

(1.北京交通大學(xué)電子信息工程學(xué)院，北京 100044；2.北京市軌道交通電磁兼容與衛(wèi)星導(dǎo)航工程技術(shù)研究中心，北京 100044)

隨著衛(wèi)星導(dǎo)航系統(tǒng)及衛(wèi)星定位技術(shù)的迅速發(fā)展，其在水路、航路和道路交通領(lǐng)域的應(yīng)用也趨于成熟，全球各地對(duì)其在鐵路交通領(lǐng)域的應(yīng)用也展開(kāi)了廣泛的研究。相對(duì)于傳統(tǒng)的軌道電路、應(yīng)答器的定位方式，衛(wèi)星定位有利于列車定位由依賴地面設(shè)備的被動(dòng)定位轉(zhuǎn)變?yōu)檐囕d設(shè)備自主定位，這能極大地減少地面設(shè)備的鋪設(shè)和安裝，節(jié)省建設(shè)和維護(hù)成本。此外，列車自主定位有助于推進(jìn)固定閉塞控制向移動(dòng)閉塞控制的升級(jí)，從而充分增強(qiáng)線路的通過(guò)能力，提高列車的運(yùn)行效率。因此，基于衛(wèi)星定位的列車定位是下一代列車運(yùn)行控制系統(tǒng)的重要發(fā)展方向，也是列控系統(tǒng)向智能化、自主化發(fā)展的必然趨勢(shì)。列車運(yùn)行在既定的軌道線路上，線路的地理和幾何信息是數(shù)字軌道地圖的關(guān)鍵組成部分，也是衛(wèi)星定位應(yīng)用于列車定位的關(guān)鍵基礎(chǔ)，它是將列車的空間位置坐標(biāo)匹配對(duì)應(yīng)到軌道一維坐標(biāo)的重要依據(jù)。

數(shù)字軌道地圖數(shù)據(jù)的主要來(lái)源是軌道線路的測(cè)量，目前主要采用實(shí)時(shí)動(dòng)態(tài)差分RTK(Real Time Kinematic)技術(shù)對(duì)軌道中心線進(jìn)行測(cè)量，獲得線路上一系列的離散點(diǎn)。如何利用測(cè)量數(shù)據(jù)高效準(zhǔn)確地表示線路特征，并滿足存儲(chǔ)和地圖匹配的應(yīng)用要求，是數(shù)字軌道地圖生成的重要問(wèn)題。文獻(xiàn)[1]提出用分段折線近似表示曲線軌道，并將橫向誤差作為約束條件，采用啟發(fā)式算法對(duì)測(cè)量數(shù)據(jù)進(jìn)行約簡(jiǎn)以盡可能簡(jiǎn)單高效地表示軌道。文獻(xiàn)[2]提出數(shù)字軌道地圖的分層結(jié)構(gòu)，對(duì)測(cè)量數(shù)據(jù)進(jìn)行剔除、約簡(jiǎn)后選取特征點(diǎn)在低尺度上描述線路，并通過(guò)插值在高尺度上細(xì)化數(shù)字軌道地圖。文獻(xiàn)[3]根據(jù)主曲線、優(yōu)化理論和二分法等思想提出了三種用數(shù)據(jù)點(diǎn)和線段表示線路的方法，其中自適應(yīng)半徑算法綜合性能較優(yōu)。上述方法都采用一系列離散的點(diǎn)及其構(gòu)成的折線描述線路，雖然能在一定程度上描述線路情況，但無(wú)法準(zhǔn)確地捕捉到線路的幾何線形特征，且地圖的精度及地圖匹配的效率會(huì)受制于地圖數(shù)據(jù)點(diǎn)的個(gè)數(shù)。文獻(xiàn)[4]將里程作為參數(shù)，采用三次多項(xiàng)式對(duì)線路平面的兩個(gè)方向分別進(jìn)行動(dòng)態(tài)分段擬合，該算法能保證線路幾何的連續(xù)與光滑，但其與實(shí)際線路的線形不完全一致，且對(duì)地圖匹配算法的應(yīng)用具有較高的計(jì)算要求?？紤]到鐵路線路在設(shè)計(jì)時(shí)就采用三種基本線形：直線、圓曲線和緩和曲線，因此可以從這三種線形的基本特點(diǎn)出發(fā)，對(duì)測(cè)量數(shù)據(jù)進(jìn)行相應(yīng)的特征識(shí)別及線形分段，進(jìn)而對(duì)不同的線形元素分別進(jìn)行擬合，盡可能遵循設(shè)計(jì)的思想完成數(shù)字軌道地圖的平面幾何線形的準(zhǔn)確描述。

針對(duì)線形元素識(shí)別方面，主要有方位角和曲率識(shí)別兩種方法，文獻(xiàn)[5]采用啟發(fā)式算法對(duì)里程-方位角圖中的直線、圓曲線和緩和曲線段分別用橫直線、二次多項(xiàng)式和斜直線擬合，并在特征點(diǎn)處一階導(dǎo)數(shù)連續(xù)作為約束條件，通過(guò)不斷調(diào)整特征點(diǎn)的位置尋找最優(yōu)解，此方法計(jì)算復(fù)雜，且需事先選好初始特征點(diǎn)。文獻(xiàn)[6]首先用三次樣條曲線擬合測(cè)點(diǎn)，再對(duì)擬合后的曲率進(jìn)行平滑和截?cái)嗵幚?，然后根?jù)曲率識(shí)別線形及特征點(diǎn)，該方法計(jì)算步驟較多且曲率截?cái)嗵幚硇韪鶕?jù)實(shí)際情況調(diào)整參數(shù)。文獻(xiàn)[7-8]基于方位角計(jì)算曲率，根據(jù)曲率值分辨直線段和圓曲線段，根據(jù)曲率變化率確定緩和曲線段，通過(guò)曲率的直方圖確定圓曲線的半徑，但該方法忽略了線路中存在相同半徑圓曲線段的情況，且結(jié)果跟直方圖選取的區(qū)間大小相關(guān)。在線形擬合方面，主要有樣條曲線法[6，9]、分段多項(xiàng)式曲線法[10]、最小二乘法[7-8]，這些方法都比較成熟且能精確描述線形。

綜上所述，傳統(tǒng)的用離散點(diǎn)及其構(gòu)成的折線進(jìn)行地圖描述的方法簡(jiǎn)單易實(shí)現(xiàn)，但無(wú)法保證線形的連續(xù)與平滑；采用線形描述的方法能準(zhǔn)確描述地圖，但存在計(jì)算步驟過(guò)于繁瑣、參數(shù)調(diào)整過(guò)多、應(yīng)用場(chǎng)景受限的問(wèn)題，因此本文根據(jù)數(shù)字軌道地圖的平面線形描述的需求，提出一種平面線形的特征提取方法。考慮到先擬合再計(jì)算測(cè)點(diǎn)曲率的步驟復(fù)雜，且計(jì)算結(jié)果易因測(cè)點(diǎn)誤差造成曲率的突變，而方位角變化相對(duì)平緩，算法先通過(guò)計(jì)算方位角再近似估計(jì)曲率，而后通過(guò)設(shè)定曲率閾值對(duì)于線形進(jìn)行初步識(shí)別，根據(jù)橫向誤差的要求對(duì)直線段和圓曲線段進(jìn)行迭代直至特征點(diǎn)確定，最后采用整體最小二乘法分別對(duì)直線、圓曲線和緩和曲線擬合得到相應(yīng)的線形參數(shù)，盡可能按照線路的線形特征描述地圖。

1 平面線形特征識(shí)別方法

因?yàn)槭艿匦巍⒌刭|(zhì)、技術(shù)等因素的限制，鐵路線路無(wú)法用一條長(zhǎng)直線延續(xù)始終，其方向需要進(jìn)行改變。出于行車安全考慮，在轉(zhuǎn)向處需要將相鄰的兩段直線用曲線連接起來(lái)，這種曲線稱為平面曲線。平面曲線根據(jù)性質(zhì)分為兩種：圓曲線和緩和曲線。圓曲線是具有一定半徑的圓弧，其半徑大小根據(jù)車速的不同有相應(yīng)的設(shè)計(jì)規(guī)范；為實(shí)現(xiàn)曲率半徑的逐漸變化，減少對(duì)車輛的沖擊，直線和圓曲線之間采用緩和曲線平滑過(guò)渡，其半徑由直線半徑(無(wú)窮大)逐漸減小為圓曲線半徑，設(shè)計(jì)中常采用回旋曲線。線路的平面線形基本組合為“直線-前緩和曲線-圓曲線-后緩和曲線”。

方位角和曲率是直觀描述平面線形的兩個(gè)量。不同線形元素的方位角存在以下特點(diǎn)：直線段的方位角為一定值，圓曲線的方位角根據(jù)轉(zhuǎn)彎方向線性遞增或遞減，緩和曲線的方位角呈二次拋物線變化。而曲率是切線方位角的一階導(dǎo)數(shù)，因此直線段上曲率恒為零、圓曲線段的曲率是非零常數(shù)1/R(半徑的倒數(shù))，緩和曲線段上某點(diǎn)的曲率則與其長(zhǎng)度成正比，即呈線性變化，見(jiàn)圖1。

線路中心線各測(cè)點(diǎn)的方位角是從該點(diǎn)正北方向線算起，按照順時(shí)針?lè)较蛑猎擖c(diǎn)方向線之間的水平夾角，范圍為0°～360°。本文采用圓弧切線法計(jì)算測(cè)點(diǎn)的方位角，首先利用“三點(diǎn)定圓法”計(jì)算圓心和半徑，然后將中間點(diǎn)在圓弧的切線與正北方向的夾角作為該點(diǎn)的方位角，見(jiàn)圖2，公式為

( 1 )

式中：αi為i點(diǎn)的方位角；(xp，yp)為i+1點(diǎn)在切線方向上的投影點(diǎn)坐標(biāo)；(xi，yi)為i點(diǎn)坐標(biāo)。

圖2 圓弧切線法計(jì)算方位角

測(cè)點(diǎn)的曲率可以通過(guò)方位角的一階導(dǎo)數(shù)近似進(jìn)行計(jì)算。曲率為

( 2 )

式中：ρi為i點(diǎn)的曲率；αi、αi+1分別為i、i+1點(diǎn)的方位角；dis(i，i+1)為i，i+1兩點(diǎn)間的距離。

根據(jù)直線段曲率恒為零、曲線段曲率非零的特點(diǎn)，選取相應(yīng)的曲率閾值即可將測(cè)點(diǎn)初步分為直線段和曲線段。此外，本文引入橫向誤差作為約束條件以確保直線段分類的準(zhǔn)確性，并進(jìn)一步確定圓曲線和緩和曲線段的范圍。

橫向誤差是測(cè)點(diǎn)到相應(yīng)擬合線形的正交投影距離，對(duì)應(yīng)數(shù)字軌道地圖在垂直股道方向上的誤差。鐵路上平行股道的線路中心線間距為5 m左右，若以中線2.5 m作為列車股道占用識(shí)別的分界線，并且作為參考的數(shù)字軌道地圖精度應(yīng)高一個(gè)數(shù)量級(jí)，則數(shù)字軌道地圖本身的橫向誤差不應(yīng)超過(guò)0.25 m。

算法的流程如下：首先根據(jù)軌道線路中心的測(cè)點(diǎn)坐標(biāo)計(jì)算方位角，由方位角計(jì)算曲率，考慮到測(cè)點(diǎn)誤差對(duì)計(jì)算的影響，本文采用滑動(dòng)平均算法對(duì)方位角和曲率進(jìn)行平滑濾波；根據(jù)直線段曲率為零的特點(diǎn)，設(shè)定閾值將各測(cè)點(diǎn)初步分類為直線段和曲線段；以橫向誤差為約束條件，采用相應(yīng)的方法依次對(duì)直線段和曲線段進(jìn)行迭代擬合，確定直線、圓曲線和緩和曲線的起始點(diǎn)、終止點(diǎn)和相關(guān)參數(shù)，最終得到線路整體的平面線形。算法的實(shí)現(xiàn)流程見(jiàn)圖3。

圖3 平面線形擬合算法流程

2 平面線形擬合方法

確定測(cè)點(diǎn)的平面線形之后，需對(duì)不同的線形元素進(jìn)行擬合，求解線形的相關(guān)參數(shù)。通常根據(jù)線形的特征選取相應(yīng)的數(shù)學(xué)模型，將測(cè)點(diǎn)的坐標(biāo)作為觀測(cè)量，線形的參數(shù)作為待估計(jì)量，一般采用最小二乘法求解。

經(jīng)典的最小二乘法都以y為因變量，以x為自變量，假設(shè)因變量有誤差，自變量無(wú)誤差。這會(huì)導(dǎo)致選取的自變量和因變量不同時(shí)，得到的擬合結(jié)果也不同。在實(shí)際線路測(cè)量中，由于儀器、模型等存在誤差，各測(cè)點(diǎn)的橫縱坐標(biāo)都不可避免地存在誤差，整體最小二乘法就是在同時(shí)考慮x、y誤差的情況下求解最優(yōu)的參數(shù)估計(jì)，其數(shù)學(xué)模型稱為EIV(Errors in Variables)模型。

EIV的線性數(shù)學(xué)模型為[11]

Y+eY=(A+eA)·ξ

( 3 )

式中：Y為n維觀測(cè)向量；eY為Y的隨機(jī)誤差向量；A為n×m維系數(shù)矩陣；eA為系數(shù)矩陣A的隨機(jī)誤差矩陣；ξ為待求的m×1維參數(shù)向量。

當(dāng)觀測(cè)向量和系數(shù)矩陣為等精度觀測(cè)時(shí)，eY和eA獨(dú)立且都服從零均值定方差的高斯分布。

將式( 3 )變形為

( 4 )

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

整體最小二乘的約束準(zhǔn)則為

min‖eA；eY‖F(xiàn)

( 5 )

式中：‖P‖F(xiàn)為n×m維矩陣P的Frobenius范數(shù)，其定義如下

( 6 )

式中：pij為矩陣P第i行第j列的元素；tr(·)為矩陣的跡。

常用的整體最小二乘數(shù)值求解方法有奇異值分解法SVD(Singular Value Decomposition)、完全正交法、增廣矩陣法、拉格朗日迭代法等[12]。

2.1 直線線形擬合

設(shè)直線方程為yi=axi+b，i=1，2，3，…，n，(xi，yi)為n對(duì)線路中線的測(cè)點(diǎn)坐標(biāo)，a為斜率，b為截距。考慮到x和y都分別存在誤差vx和vy，方程則表示為

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

( 7 )

其誤差方程按照EIV模型可以表示為

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

( 8 )

注意到系數(shù)矩陣A中，存在元素固定為1的一列，此問(wèn)題可以轉(zhuǎn)化成混合最小二乘問(wèn)題，根據(jù)混合最小二乘法求解計(jì)算[13]，此處不再展開(kāi)。

2.2 圓曲線擬合

假設(shè)圓的圓心為(xc，yc)，半徑為r，測(cè)點(diǎn)(xi，yi)到圓弧上的距離di(等于測(cè)點(diǎn)與圓心之間的距離減去半徑)為

( 9 )

(10)

其誤差方程為

(11)

其矩陣表達(dá)式為

Aδξ=Y+eY

(12)

式中：δξ為待求參數(shù)向量的修正值。且

未知參數(shù)向量修正值的解為

δξ=(ATA)-1ATY

(13)

(14)

(15)

2.3 緩和曲線擬合

線路設(shè)計(jì)中緩和曲線常用回旋線，回旋線的基本公式為rl=B2，B為回旋線的參數(shù)?；匦€的終點(diǎn)處，r=R(R為圓曲線半徑)，l=Ls(Ls為回旋線總長(zhǎng)度)，故RLs=B2。計(jì)算時(shí)需要建立局部坐標(biāo)系，將直緩點(diǎn)定為原點(diǎn)，見(jiàn)圖4。

圖4 回旋曲線

回旋曲線的計(jì)算式為

(24)

實(shí)際應(yīng)用時(shí)，常用三次拋物線來(lái)近似計(jì)算

(25)

式中：C=B2=RLs。

局部坐標(biāo)系的計(jì)算結(jié)果經(jīng)過(guò)坐標(biāo)變換后，即可得到整體坐標(biāo)系下的表達(dá)式為

式中：xZH、yZH為直緩點(diǎn)在整體坐標(biāo)系下的絕對(duì)坐標(biāo)；xlocal、ylocal為緩和曲線上的點(diǎn)在局部坐標(biāo)系中的坐標(biāo)；β為兩個(gè)坐標(biāo)系之間的旋轉(zhuǎn)角。

在確定完直線段和圓曲線段的范圍和相關(guān)參數(shù)之后，對(duì)緩和曲線段建立局部坐標(biāo)系，按照三次拋物線對(duì)其進(jìn)行擬合，同樣采用整體最小二乘法，此處不再贅述。

3 算法測(cè)試

3.1 數(shù)據(jù)描述

本文選取青藏鐵路一段線路中心線的GPS測(cè)量數(shù)據(jù)進(jìn)行算法的測(cè)試與驗(yàn)證，數(shù)據(jù)采用Navcom SF-2050雙頻差分GPS接收機(jī)采集，精度為厘米級(jí)。該段線路全長(zhǎng)約14.7 km，原始數(shù)據(jù)經(jīng)預(yù)處理后得到6 100個(gè)測(cè)點(diǎn)，點(diǎn)間距最大3 m，最小1.5 m，平均2.4 m。

首先將測(cè)點(diǎn)經(jīng)緯度轉(zhuǎn)化為平面坐標(biāo)，本文利用通用橫軸墨卡托UTM(Universal Transverse Mercator)投影，為方便計(jì)算將起點(diǎn)轉(zhuǎn)化為坐標(biāo)原點(diǎn)，該段線路平面線形見(jiàn)圖5。

圖5 線路平面線形

3.2 線形識(shí)別與擬合結(jié)果

采用文中的“圓弧切線法”計(jì)算各測(cè)點(diǎn)的方位角，見(jiàn)圖6，由方位角可以初步判斷出該段線路包含5個(gè)曲線段，6個(gè)直線段。局部放大后，可以觀察到方位角受測(cè)量及計(jì)算誤差的影響存在噪聲，為減少噪聲對(duì)曲率計(jì)算的干擾，采用滑動(dòng)平均法濾波。

用滑動(dòng)平均后的方位角進(jìn)行曲率的計(jì)算，同樣地，對(duì)計(jì)算后的曲率進(jìn)行濾波，見(jiàn)圖7。從圖7中可以觀察到濾波后的曲率噪聲明顯減少。

(a)整體圖

(b)局部放大圖圖6 測(cè)點(diǎn)的方位角變化

圖7 測(cè)點(diǎn)的曲率變化

圖8 直線段和曲線段閾值及初始分段

考慮到線路測(cè)量數(shù)據(jù)的精度達(dá)到了厘米級(jí)，因此本文設(shè)定橫向誤差的最大值為0.1 m，并以此作為限制條件依次對(duì)直線段、曲線段進(jìn)行直線、圓曲線擬合，對(duì)剩余的緩和曲線段進(jìn)行三次拋物線擬合，最終得到線路的整體線形及相應(yīng)的參數(shù)。

擬合后線路整體的曲率圖和平面線形的分布圖分別見(jiàn)圖9、圖10。

3.3 算法的性能指標(biāo)

橫向誤差Eacross：所有測(cè)點(diǎn)到相應(yīng)擬合線形的正交投影距離，包括最大值和平均值，其值表征了擬合的平面線形精確度及準(zhǔn)確度，越小表明擬合線形在垂直股道方向上偏差越小。

圖9 擬合后曲率

圖10 擬合后線路的整體平面線形

里程誤差Ealong：所有測(cè)點(diǎn)的原始累積里程值與擬合后各測(cè)點(diǎn)在縱向上的累積里程值的對(duì)比，其值表明了擬合線路在縱向上的長(zhǎng)度損失。

本次算法的性能指標(biāo)結(jié)果見(jiàn)表1。

表1 算法性能指標(biāo)結(jié)果

由以上結(jié)果統(tǒng)計(jì)，總結(jié)出以下結(jié)論：

(1)算法可以有效提取出線路的三種線形并對(duì)其進(jìn)行擬合，通過(guò)極少數(shù)關(guān)鍵點(diǎn)及相應(yīng)的參數(shù)對(duì)整段線路的平面線形進(jìn)行表達(dá)，本文可以用22個(gè)線形的分界點(diǎn)及線形的參數(shù)表達(dá)出6 100個(gè)測(cè)點(diǎn)的信息，數(shù)據(jù)的約簡(jiǎn)率極高。

(2)平面線形的擬合結(jié)果表明算法的精度很高，在0.1 m的橫向誤差約束下，最大的橫向誤差不超過(guò)0.083 m，滿足列車定位及股道識(shí)別的精度要求；且算法的里程累積誤差最大僅為0.015 m，可以忽略不計(jì)。

(3)本算法擬合出來(lái)的線形不僅可以用來(lái)計(jì)算線路的坐標(biāo)，還能準(zhǔn)確描述出線路的方位角及曲率變化，豐富了數(shù)字軌道地圖的內(nèi)容。

4 結(jié)束語(yǔ)

高精度的數(shù)字軌道地圖是GNSS鐵路應(yīng)用的重要基礎(chǔ)，本文采用基于方位角的曲率方法對(duì)線路的平面線形進(jìn)行特征識(shí)別與分段擬合，通過(guò)對(duì)青藏線實(shí)測(cè)線路數(shù)據(jù)的測(cè)試，驗(yàn)證了算法的可行性及高精度，并且相對(duì)于傳統(tǒng)的一系列離散點(diǎn)的地圖表示方法，本算法可以重構(gòu)出線路的方位角及曲率變化情況，使數(shù)字軌道地圖的平面線形更加準(zhǔn)確完善，滿足應(yīng)用需求。