基于GNSS定位和衛(wèi)星影像的鐵路中線數(shù)據(jù)生成方法研究

馬潤(rùn)澤 曲以勝 李 克 李立改 何 為

(1.中國(guó)科學(xué)院微系統(tǒng)與信息技術(shù)研究所，上海 201808；2.中國(guó)鐵路烏魯木齊局集團(tuán)有限公司，烏魯木齊 830011；3.中國(guó)鐵道科學(xué)研究院集團(tuán)有限公司，北京 100081)

隨著全球?qū)Ш叫l(wèi)星系統(tǒng)(Global Navigation Satellite System，GNSS)和衛(wèi)星遙感測(cè)量技術(shù)的發(fā)展，衛(wèi)星定位和衛(wèi)星影像技術(shù)越來(lái)越多應(yīng)用于工程測(cè)繪領(lǐng)域[1-2]。其中，實(shí)時(shí)動(dòng)態(tài)差分(Real Time Kinematic，RTK)定位精度可達(dá)厘米級(jí)至毫米級(jí)[3]，此外還有操作便捷、測(cè)量效率高等優(yōu)勢(shì)，近年來(lái)在鐵路工程測(cè)繪中得到廣泛應(yīng)用[4-5]。

鐵路軌道中線數(shù)據(jù)是鐵路工程測(cè)繪中的重要信息，隨著我國(guó)北斗衛(wèi)星導(dǎo)航系統(tǒng)的推廣應(yīng)用，未來(lái)的列車(chē)控制系統(tǒng)將采用基于衛(wèi)星定位的多傳感融合定位技術(shù)[6-7]，以減少軌旁設(shè)備，如我國(guó)青藏鐵路ITCS列控系統(tǒng)采用衛(wèi)星及軌道電子地圖進(jìn)行列車(chē)定位，此時(shí)需要對(duì)軌道線路進(jìn)行測(cè)量和電子地圖制作，以實(shí)現(xiàn)列車(chē)定位的軌道綁定、虛擬閉塞信號(hào)觸發(fā)等功能[8-9]。

在鐵路未開(kāi)通前，無(wú)法采用車(chē)載定位移動(dòng)獲取數(shù)據(jù)的方式，多采用人工持衛(wèi)星定位設(shè)備上道采集經(jīng)緯高數(shù)據(jù)。梁旺等采用千尋定位GNSS RTK技術(shù)測(cè)量鐵軌中線，獲得比傳統(tǒng)單基站GNSS RTK更優(yōu)的定位精度[10]。如果要獲得密度更高的數(shù)據(jù)點(diǎn)，可以縮短采樣間隔，或者采用數(shù)據(jù)插值方法填充數(shù)據(jù)。龍明濤等針對(duì)HJ-1衛(wèi)星影像數(shù)據(jù)特點(diǎn)，提出一種具有針對(duì)性的快速插值算法，使得HJ-1影像的每個(gè)像元都具有經(jīng)緯度信息[11]，該算法是一種平面插值算法，運(yùn)算時(shí)不會(huì)使用高度數(shù)據(jù)。

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

1 地理數(shù)據(jù)坐標(biāo)轉(zhuǎn)換

衛(wèi)星定位接收機(jī)輸出結(jié)果為大地坐標(biāo)，即經(jīng)度(λ)、緯度(L)、海拔高度(h)，該數(shù)據(jù)中，假設(shè)地球?yàn)樾D(zhuǎn)橢球體，地球自轉(zhuǎn)軸(極軸)與一橢圓短半軸重合，橢圓的橢圓度(扁率)為1/298.26(WGS84坐標(biāo)系)[12]，橢圓繞其短半軸旋轉(zhuǎn)構(gòu)成橢球體的表面，該描述中地球赤道是圓的，旋轉(zhuǎn)橢球和子午圈橢圓示意見(jiàn)圖1、圖2。

圖1 旋轉(zhuǎn)橢球基本概念示意

圖2 子午圈橢圓示意

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

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

(1)

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

(2)

λ=arctan2(y，x)

(3)

其中，arctan2()為計(jì)算給定橫、縱坐標(biāo)點(diǎn)的反正切函數(shù)，取-180°～180°；緯度L通過(guò)如下迭代式求解

(4)

迭代初值t0=0，經(jīng)過(guò)5～6次迭代，可得到足夠精度的ti，進(jìn)而計(jì)算緯度L和海拔高度h，有

L=arctan(ti)

(5)

(6)

2 3D曲線插值方法

為提高測(cè)量效率，經(jīng)緯高數(shù)據(jù)的采樣測(cè)量間隔可選擇在50 m以上，并在兩個(gè)相鄰測(cè)量點(diǎn)之間填充數(shù)據(jù)點(diǎn)，使得這些數(shù)據(jù)點(diǎn)平滑連續(xù)。鐵軌路線可分為直線段和曲線段兩類(lèi)，據(jù)此將這些采樣數(shù)據(jù)點(diǎn)分為對(duì)應(yīng)的兩類(lèi)。提出一整套判斷數(shù)據(jù)點(diǎn)類(lèi)型和數(shù)據(jù)點(diǎn)間插值的算法。

2.1 直線點(diǎn)與曲線點(diǎn)判定方法

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

(7)

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

(8)

(9)

(10)

當(dāng)滿(mǎn)足如下條件式時(shí)

(11)

即可判定Pl屬于直線點(diǎn)，否則判定Pl屬于曲線點(diǎn)。其中，ε1、ε2為判定直線的角度經(jīng)驗(yàn)閾值，一般可設(shè)置為0.2°；對(duì)Pl判定后，再用同樣方法繼續(xù)判定Pl+1、Pl+2…等，直至判定所有數(shù)據(jù)點(diǎn)。

2.2 直線點(diǎn)間插值

對(duì)于相鄰兩個(gè)直線點(diǎn)間的數(shù)據(jù)插值，可采用線性插值的方法。假設(shè)Pl、Pl+1為相鄰的兩個(gè)直線點(diǎn)，坐標(biāo)分別為(xl，yl，zl)和(xl+1，yl+1，zl+1)，如果要使插入數(shù)據(jù)點(diǎn)的間距不超過(guò)δd，則應(yīng)先計(jì)算Pl、Pl+1的間距d，有

(12)

計(jì)算插入數(shù)據(jù)點(diǎn)的個(gè)數(shù)n，有

(13)

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

(14)

通過(guò)該方法，可完成對(duì)所有相鄰直線點(diǎn)間的數(shù)據(jù)擬合插值。

2.3 曲線點(diǎn)間插值

對(duì)曲線點(diǎn)間插值，提出一種基于改進(jìn)3D貝塞爾曲線的插值算法，其示意見(jiàn)圖3。

圖3 曲線點(diǎn)插值示意

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

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

P2到PN-1間的數(shù)據(jù)插值，采用曲線插值方法。貝塞爾曲線插值是應(yīng)用廣泛的曲線插值方法[14]。給定點(diǎn)T0、T1、…、Tn，其貝塞爾曲線表達(dá)式為

(15)

圖4 初始控制點(diǎn)計(jì)算示意

(16)

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

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

(17)

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

關(guān)于從P2到PN-1間插值點(diǎn)數(shù)l的選擇，可先計(jì)算P2到PN-1間折線段連線的長(zhǎng)度，有

(18)

P2到PN-1間插值曲線的長(zhǎng)度較折線段長(zhǎng)，取極限長(zhǎng)度為2d，假設(shè)要求插入數(shù)據(jù)點(diǎn)的間距不超過(guò)δd，則插入點(diǎn)數(shù)l為

(19)

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

3 衛(wèi)星影像圖修正

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

3.1 墨卡托投影

目前，主要采用墨卡托投影的方法將經(jīng)緯度點(diǎn)映射到衛(wèi)星圖上。又稱(chēng)為“等角正軸圓柱”投影。其基本原理是假設(shè)有一個(gè)在赤道與地球相切的圓柱體，先把橢球面映射到圓柱體表面，然后展開(kāi)圓柱面，即實(shí)現(xiàn)了球平轉(zhuǎn)換。該投影具有等角特性，在保證對(duì)象的形狀不會(huì)改變的同時(shí)，也保證了方向和相互位置的正確性，常應(yīng)用在航海和航空領(lǐng)域。

墨卡托投影對(duì)經(jīng)緯度的轉(zhuǎn)化方法[17]是，地球表面上某點(diǎn)A(φ，λ)經(jīng)過(guò)墨卡托投影得到新坐標(biāo)點(diǎn)A′(x，y)。其中，φ0為標(biāo)準(zhǔn)緯度；λ0為標(biāo)準(zhǔn)經(jīng)度；e為第一偏心率；e′為第二偏心率；a為長(zhǎng)半軸；b為短半軸，投影公式為

(20)

(21)

3.2 插值曲線參數(shù)優(yōu)化

浩吉鐵路北起內(nèi)蒙，南達(dá)江西，全長(zhǎng)1 814 km，是我國(guó)重要的貨運(yùn)鐵路。為推動(dòng)北斗導(dǎo)航衛(wèi)星技術(shù)應(yīng)用于鐵路CTC系統(tǒng)，國(guó)鐵集團(tuán)組織多家單位在試驗(yàn)區(qū)段驗(yàn)證衛(wèi)星定位技術(shù)，其中，對(duì)鐵道中線的GNSS衛(wèi)星定位測(cè)量場(chǎng)景見(jiàn)圖5。采用在浩吉鐵路坪田到瀏陽(yáng)東區(qū)段實(shí)際采集的高精度衛(wèi)星差分定位數(shù)據(jù)，測(cè)試所提出曲線插值算法的實(shí)際效果。

圖5 鐵軌衛(wèi)星定位測(cè)點(diǎn)

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

其中，在一段弧形鐵路上采集了長(zhǎng)約1 km的數(shù)據(jù)。基于這些經(jīng)緯高測(cè)點(diǎn)數(shù)據(jù)，應(yīng)用前述的3D坐標(biāo)轉(zhuǎn)換算法、曲線插值算法等，得到更密集的數(shù)據(jù)點(diǎn)。插值計(jì)算結(jié)果見(jiàn)圖6。

圖6 數(shù)據(jù)插值結(jié)果

由圖6可知，19個(gè)采樣測(cè)點(diǎn)平均間隔約50 m，數(shù)據(jù)插值點(diǎn)設(shè)置為間距不超過(guò)5 m，p取0.9。此時(shí)，插值點(diǎn)形成的擬合曲線平滑穿過(guò)了各間隔采樣點(diǎn)。

然后，選取不同的p值分別進(jìn)行插值計(jì)算，并將這些采樣點(diǎn)和數(shù)據(jù)插值點(diǎn)投影到衛(wèi)星影像上，各插值結(jié)果見(jiàn)圖7。

圖7 不同參數(shù)值下的插值結(jié)果

圖7中，紅色箭頭代表原始的經(jīng)緯高測(cè)點(diǎn)數(shù)據(jù)(間隔50 m)，綠色箭頭代表按前述算法得到的數(shù)據(jù)插值結(jié)果。由圖7可知，當(dāng)選擇不同的參數(shù)p時(shí)，將插值結(jié)果與衛(wèi)星影像進(jìn)行對(duì)比，可發(fā)現(xiàn)p=0.9時(shí)兩者最為接近，即可認(rèn)為0.9是最優(yōu)參數(shù)值。再?gòu)亩拷嵌扔?jì)算不同p值下的插值數(shù)據(jù)誤差情況，誤差結(jié)果見(jiàn)表1。

表1 不同p值下的插值數(shù)據(jù)誤差

由表1可知，p值為0.9時(shí)，最小誤差達(dá)到0.01 m，平均誤差0.56 m，誤差中位值0.24 m，優(yōu)于其他p值。不難看出，通過(guò)衛(wèi)星影像觀察到的最優(yōu)p值，通過(guò)定量誤差分析后仍為是最優(yōu)。此誤差計(jì)算涉及大量的空間三維點(diǎn)距離計(jì)算，對(duì)于小批量數(shù)據(jù)點(diǎn)計(jì)算尚可接受，對(duì)大批量數(shù)據(jù)點(diǎn)的比較計(jì)算，則是沉重的負(fù)擔(dān)，而通過(guò)衛(wèi)星影像的值觀比較修正方式，可以避免這種情況的發(fā)生。

此外，如果經(jīng)緯高衛(wèi)星測(cè)點(diǎn)采樣間距拉大(如間隔100,150 m等)，也可以用同樣方法進(jìn)行數(shù)據(jù)插值，插值結(jié)果見(jiàn)圖8、圖9。

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

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

由圖8、圖9可知，對(duì)原始采樣數(shù)據(jù)點(diǎn)，間隔50，100，150 m的插值結(jié)果基本吻合，再進(jìn)行定量誤差計(jì)算，結(jié)果見(jiàn)表2。

表2 不同衛(wèi)星定位采樣間隔下的插值數(shù)據(jù)誤差 m

由表2可知，在不同的采樣間隔條件下，插值結(jié)果誤差相差不大，這說(shuō)明拉大衛(wèi)星采樣測(cè)點(diǎn)的間距可以節(jié)省測(cè)繪時(shí)間、提高測(cè)繪效率。

4 結(jié)論

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

(2)采用該方法，改進(jìn)了貝塞爾曲線插值方法，采用在三維空間中計(jì)算多個(gè)控制點(diǎn)的方式，從而可以使插值結(jié)果曲線平滑的通過(guò)所有經(jīng)緯度測(cè)點(diǎn)。

(3)采用該方法，對(duì)間隔更大的原始數(shù)據(jù)測(cè)點(diǎn)也可以獲得較好的數(shù)據(jù)插值結(jié)果，在采樣間隔100 m情況下達(dá)到0.56 m平均插值誤差，在間隔150 m情況下達(dá)到0.83 m平均誤差，有利于減少衛(wèi)星經(jīng)緯高測(cè)點(diǎn)的工作量，提高測(cè)繪工作效率。