基于GNSS定位和衛星影像的鐵路中線數據生成方法研究

馬潤澤 曲以勝 李 克 李立改 何 為

(1.中國科學院微系統與信息技術研究所，上海 201808；2.中國鐵路烏魯木齊局集團有限公司，烏魯木齊 830011；3.中國鐵道科學研究院集團有限公司，北京 100081)

隨著全球導航衛星系統(Global Navigation Satellite System，GNSS)和衛星遙感測量技術的發展，衛星定位和衛星影像技術越來越多應用于工程測繪領域[1-2]。其中，實時動態差分(Real Time Kinematic，RTK)定位精度可達厘米級至毫米級[3]，此外還有操作便捷、測量效率高等優勢，近年來在鐵路工程測繪中得到廣泛應用[4-5]。

鐵路軌道中線數據是鐵路工程測繪中的重要信息，隨著我國北斗衛星導航系統的推廣應用，未來的列車控制系統將采用基于衛星定位的多傳感融合定位技術[6-7]，以減少軌旁設備，如我國青藏鐵路ITCS列控系統采用衛星及軌道電子地圖進行列車定位，此時需要對軌道線路進行測量和電子地圖制作，以實現列車定位的軌道綁定、虛擬閉塞信號觸發等功能[8-9]。

在鐵路未開通前，無法采用車載定位移動獲取數據的方式，多采用人工持衛星定位設備上道采集經緯高數據。梁旺等采用千尋定位GNSS RTK技術測量鐵軌中線，獲得比傳統單基站GNSS RTK更優的定位精度[10]。如果要獲得密度更高的數據點，可以縮短采樣間隔，或者采用數據插值方法填充數據。龍明濤等針對HJ-1衛星影像數據特點，提出一種具有針對性的快速插值算法，使得HJ-1影像的每個像元都具有經緯度信息[11]，該算法是一種平面插值算法，運算時不會使用高度數據。

以下通過離散(間隔幾十米至上百米)的經緯高數據采樣點，采用優化的3D空間曲線擬合插值方法，插值生成所需密度(分辨率)的鐵軌中線數據，再通過衛星影像圖進行校正。最后，通過浩吉鐵路的實際采集數據，驗證所提出方法的有效性。

1 地理數據坐標轉換

衛星定位接收機輸出結果為大地坐標，即經度(λ)、緯度(L)、海拔高度(h)，該數據中，假設地球為旋轉橢球體，地球自轉軸(極軸)與一橢圓短半軸重合，橢圓的橢圓度(扁率)為1/298.26(WGS84坐標系)[12]，橢圓繞其短半軸旋轉構成橢球體的表面，該描述中地球赤道是圓的，旋轉橢球和子午圈橢圓示意見圖1、圖2。

圖1 旋轉橢球基本概念示意

圖2 子午圈橢圓示意

圖中，λ為經度；φ為地心緯度；L為常用的地理緯度(簡稱緯度)；Re和Rp分別為橢圓的長半軸和短半軸。在曲線插值步驟中，需要使用右手直角坐標系，也稱為地心地固坐標系(Earth-Center Earth-Fixed，ECEF)，坐標原點選在地心，oeze為自轉軸且指向北極，oexe軸指向赤道與本初子午線的交點，oeye軸在赤道平面且指向90°經線，ECEF系與地球固連。

通過人工持衛星定位設備上道測量，可以獲得鐵軌道心的順次測量的地理坐標集(λi，Li，hi)i=1，2，…，n，再轉化為地心直角坐標集(xi，yi，zi)i=1，2，…，n；通過曲線插值方法生成更密集的空間直角坐標集(xj，yj，zj)j=1，2，…，N(N?n)，再轉換回地理坐標(λj，Lj，hj)j=1，2，…，N。其坐標相互轉換的關系式為[13]

(1)

式中，e為地球橢圓偏心率；RN為卯酉圈曲率半徑；計算式為

(2)

λ=arctan2(y，x)

(3)

其中，arctan2()為計算給定橫、縱坐標點的反正切函數，取-180°～180°；緯度L通過如下迭代式求解

(4)

迭代初值t0=0，經過5～6次迭代，可得到足夠精度的ti，進而計算緯度L和海拔高度h，有

L=arctan(ti)

(5)

(6)

2 3D曲線插值方法

為提高測量效率，經緯高數據的采樣測量間隔可選擇在50 m以上，并在兩個相鄰測量點之間填充數據點，使得這些數據點平滑連續。鐵軌路線可分為直線段和曲線段兩類，據此將這些采樣數據點分為對應的兩類。提出一整套判斷數據點類型和數據點間插值的算法。

2.1 直線點與曲線點判定方法

α1=arctan2(yl-yl-1，xl-xl-1)

(7)

α2=arctan2(yl+1-yl，xl+1-xl)

(8)

(9)

(10)

當滿足如下條件式時

(11)

即可判定Pl屬于直線點，否則判定Pl屬于曲線點。其中，ε1、ε2為判定直線的角度經驗閾值，一般可設置為0.2°；對Pl判定后，再用同樣方法繼續判定Pl+1、Pl+2…等，直至判定所有數據點。

2.2 直線點間插值

對于相鄰兩個直線點間的數據插值，可采用線性插值的方法。假設Pl、Pl+1為相鄰的兩個直線點，坐標分別為(xl，yl，zl)和(xl+1，yl+1，zl+1)，如果要使插入數據點的間距不超過δd，則應先計算Pl、Pl+1的間距d，有

(12)

計算插入數據點的個數n，有

(13)

其中，ceil()為向上取整函數，如ceil(3)=3、ceil(3，1)=4等。進而，得到插入的數據點集(xj，yj，zj)(j=1，2，…，n)，有

(14)

通過該方法，可完成對所有相鄰直線點間的數據擬合插值。

2.3 曲線點間插值

對曲線點間插值，提出一種基于改進3D貝塞爾曲線的插值算法，其示意見圖3。

圖3 曲線點插值示意

假設虛線為鐵軌道心實際連線，圓形點為曲線點，方形點為直線點。由圖3可知，曲線點間插值，就是在P1到PN點間生成所需間距的密集數據點，要求形成的曲線平滑且經過P1到PN各點。設P1到PN各點的坐標為(xi，yi，zi)(i=1，2，…，N)。

首先，對于P1和P2、PN-1和PN間的插值，可采用前述的“直線點間”插值的方法。

P2到PN-1間的數據插值，采用曲線插值方法。貝塞爾曲線插值是應用廣泛的曲線插值方法[14]。給定點T0、T1、…、Tn，其貝塞爾曲線表達式為

(15)

圖4 初始控制點計算示意

(16)

根據上式，計算得到兩交點sj、tj后，以點C0=(sj+tj)/2作為初始控制點，根據C0和點P3到PN-2，得到貝塞爾插值計算的N-4個控制點Ci，即

Ci=C0+(Pi-C0)p，i=3，4，…,N-2

(17)

其中，0≤p≤1，代表控制點Ci距離C0和Pi間的相對遠近程度，是一個用于調節插值曲線形狀的參數；獲得Ci后，將P2、C3、C4、…、CN-2、PN-1(共N-2個點)作為貝塞爾曲線插值的輸入控制點，通過調節參數p，即可獲得所需的曲線插值點。

關于從P2到PN-1間插值點數l的選擇，可先計算P2到PN-1間折線段連線的長度，有

(18)

P2到PN-1間插值曲線的長度較折線段長，取極限長度為2d，假設要求插入數據點的間距不超過δd，則插入點數l為

(19)

其中，ceil()為向上取整函數。如此可保證插入的l點數的插值點，彼此之間的間距不超過δd。

3 衛星影像圖修正

用前述方法獲得的曲線插值點(經度緯度)，可在衛星影像圖上顯示其曲線軌跡，并借助影像圖對插值點做進一步算法參數校正，以獲得更準確的鐵軌中線位置數據。為保障位置修正的準確性，需要衛星影像的分辨率小于插值生成數據的間隔δd，如高分2號衛星影像分辨率達到1 m[16]，可滿足間隔1 m以上的插值數據修正。

3.1 墨卡托投影

目前，主要采用墨卡托投影的方法將經緯度點映射到衛星圖上。又稱為“等角正軸圓柱”投影。其基本原理是假設有一個在赤道與地球相切的圓柱體，先把橢球面映射到圓柱體表面，然后展開圓柱面，即實現了球平轉換。該投影具有等角特性，在保證對象的形狀不會改變的同時，也保證了方向和相互位置的正確性，常應用在航海和航空領域。

墨卡托投影對經緯度的轉化方法[17]是，地球表面上某點A(φ，λ)經過墨卡托投影得到新坐標點A′(x，y)。其中，φ0為標準緯度；λ0為標準經度；e為第一偏心率；e′為第二偏心率；a為長半軸；b為短半軸，投影公式為

(20)

(21)

3.2 插值曲線參數優化

浩吉鐵路北起內蒙，南達江西，全長1 814 km，是我國重要的貨運鐵路。為推動北斗導航衛星技術應用于鐵路CTC系統，國鐵集團組織多家單位在試驗區段驗證衛星定位技術，其中，對鐵道中線的GNSS衛星定位測量場景見圖5。采用在浩吉鐵路坪田到瀏陽東區段實際采集的高精度衛星差分定位數據，測試所提出曲線插值算法的實際效果。

圖5 鐵軌衛星定位測點

使用諾瓦泰衛導板卡(OEM718D)采集經緯度高度差分定位數據(精度達到1 cm)，以這些經緯高數據作為輸入，執行前述3D曲線插值算法，得到足夠密集的鐵道中線數據插值點集，將這些插值點按墨卡托投影計算式投影到谷歌衛星影像上，與圖中的鐵道衛星影像進行對比，并調節曲線參數p，以得到與衛星影像最接近的插值曲線。

其中，在一段弧形鐵路上采集了長約1 km的數據。基于這些經緯高測點數據，應用前述的3D坐標轉換算法、曲線插值算法等，得到更密集的數據點。插值計算結果見圖6。

圖6 數據插值結果

由圖6可知，19個采樣測點平均間隔約50 m，數據插值點設置為間距不超過5 m，p取0.9。此時，插值點形成的擬合曲線平滑穿過了各間隔采樣點。

然后，選取不同的p值分別進行插值計算，并將這些采樣點和數據插值點投影到衛星影像上，各插值結果見圖7。

圖7 不同參數值下的插值結果

圖7中，紅色箭頭代表原始的經緯高測點數據(間隔50 m)，綠色箭頭代表按前述算法得到的數據插值結果。由圖7可知，當選擇不同的參數p時，將插值結果與衛星影像進行對比，可發現p=0.9時兩者最為接近，即可認為0.9是最優參數值。再從定量角度計算不同p值下的插值數據誤差情況，誤差結果見表1。

表1 不同p值下的插值數據誤差

由表1可知，p值為0.9時，最小誤差達到0.01 m，平均誤差0.56 m，誤差中位值0.24 m，優于其他p值。不難看出，通過衛星影像觀察到的最優p值，通過定量誤差分析后仍為是最優。此誤差計算涉及大量的空間三維點距離計算，對于小批量數據點計算尚可接受，對大批量數據點的比較計算，則是沉重的負擔，而通過衛星影像的值觀比較修正方式，可以避免這種情況的發生。

此外，如果經緯高衛星測點采樣間距拉大(如間隔100,150 m等)，也可以用同樣方法進行數據插值，插值結果見圖8、圖9。

圖8 間隔100 m(紅色)與間隔50 m(綠色)插值結果對比

圖9 間隔150 m(紅色)與間隔50 m(綠色)插值結果對比

由圖8、圖9可知，對原始采樣數據點，間隔50，100，150 m的插值結果基本吻合，再進行定量誤差計算，結果見表2。

表2 不同衛星定位采樣間隔下的插值數據誤差 m

由表2可知，在不同的采樣間隔條件下，插值結果誤差相差不大，這說明拉大衛星采樣測點的間距可以節省測繪時間、提高測繪效率。

4 結論

(1)與傳統平面經緯度插值方法相比，基于衛星定位和影像校正的鐵路地理數據生成方法的主要特點在于將經緯高數據轉化到直角三維坐標系再做插值，從而避免在經緯度二維平面做插值時不能采用高度信息等問題。

(2)采用該方法，改進了貝塞爾曲線插值方法，采用在三維空間中計算多個控制點的方式，從而可以使插值結果曲線平滑的通過所有經緯度測點。

(3)采用該方法，對間隔更大的原始數據測點也可以獲得較好的數據插值結果，在采樣間隔100 m情況下達到0.56 m平均插值誤差，在間隔150 m情況下達到0.83 m平均誤差，有利于減少衛星經緯高測點的工作量，提高測繪工作效率。