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

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

(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鐵路應用的重要基礎，本文采用基于方位角的曲率方法對線路的平面線形進行特征識別與分段擬合，通過對青藏線實測線路數據的測試，驗證了算法的可行性及高精度，并且相對于傳統的一系列離散點的地圖表示方法，本算法可以重構出線路的方位角及曲率變化情況，使數字軌道地圖的平面線形更加準確完善，滿足應用需求。