楊開偉 劉涵宇 孫騰達(dá) 孫秀聰 徐明
(1 中國電子科技集團公司第五十四研究所,石家莊 050081)(2 北京航空航天大學(xué)宇航學(xué)院,北京 102206)
近年來,全球衛(wèi)星導(dǎo)航系統(tǒng)(GNSS)迅速發(fā)展,在人類社會生活的方方面面發(fā)揮著不可或缺的作用,涉及到測繪、交通、救援等諸多方面。隨著自動駕駛、物聯(lián)網(wǎng)、5G等新興產(chǎn)業(yè)的發(fā)展,社會生產(chǎn)生活對衛(wèi)星導(dǎo)航服務(wù)提出了更高的要求,對時空信息的需求發(fā)展為精準(zhǔn)、實時、動態(tài)和全球[1-2]。GNSS的不足之處也逐漸顯露出來。GNSS在常規(guī)情況下只能提供10m左右的定位精度,而無法滿足用戶高精度的需求;此外,導(dǎo)航衛(wèi)星信號到達(dá)地面時極其微弱,不但導(dǎo)致其極易受到干擾和欺騙,而且導(dǎo)致其難以穿透物理遮蔽,無法在室內(nèi)、水下等復(fù)雜地形中提供高精度的定位服務(wù)[1,3-4]。
導(dǎo)航增強技術(shù)是解決上述問題、提升導(dǎo)航服務(wù)性能的一種可行方案,能利用低軌導(dǎo)航增強衛(wèi)星產(chǎn)生測距信號并與中高軌GNSS系統(tǒng)協(xié)同提供導(dǎo)航服務(wù)[5]。傳統(tǒng)的導(dǎo)航增強系統(tǒng)采用的是地球靜止軌道(GEO)衛(wèi)星,但GEO衛(wèi)星軌道資源有限,而且用戶接收到的信號弱,容易受到干擾,無法滿足用戶日益增長的高精度定位需求[6]。低軌衛(wèi)星則沒有上述問題,其軌道高度低,信號的自由空間損耗低,落地信號更強,抗干擾性能更好。此外,由于低軌導(dǎo)航增強衛(wèi)星的飛行速度較快,相等時間的位置變化較大,從而使得歷元間量測信息的相關(guān)性下降,能提高系統(tǒng)的可觀測度,加快載波相位模糊度的收斂速度[4]。例如,美國的銥下一代(Iridium NEXT)星座提供的衛(wèi)星授時與定位服務(wù),已經(jīng)能夠脫離GPS單獨提供定位服務(wù),因而低軌導(dǎo)航增強衛(wèi)星可以在GPS受到干擾的情況下提供備份的導(dǎo)航服務(wù)[4]。我國的珞珈一號衛(wèi)星搭載有導(dǎo)航增強載荷,并開展了一系列低軌衛(wèi)星增強相關(guān)試驗[7]。
高精度確定用戶的位置需要導(dǎo)航衛(wèi)星高精度的位置信息,導(dǎo)航衛(wèi)星的位置信息由廣播星歷計算,因此廣播星歷的高精度計算方法得到了廣泛的研究。目前,應(yīng)用廣泛的廣播星歷參數(shù)模型有2種,主要是針對中高軌衛(wèi)星,一種是以GPS為代表的基于軌道根數(shù)的廣播星歷參數(shù)模型;一種是以格洛納斯(GLONASS)為代表的基于軌道狀態(tài)的廣播星歷參數(shù)模型[8]。GPS廣播星歷參數(shù)模型的精度和GLONASS相當(dāng),但前者的預(yù)報能力更強,用戶算法更簡單;后者的優(yōu)勢是其參數(shù)更少,占用導(dǎo)航信號的資源少。衛(wèi)星設(shè)備、通信鏈路和用戶性能等多方面的限制,要求廣播星歷模型具有參數(shù)少、外推能力強、用戶算法簡單等特點,因此GPS廣播星歷參數(shù)模型應(yīng)用更加廣泛[9]。
傳統(tǒng)的GPS系統(tǒng)采用的是包括6個開普勒軌道參數(shù)、3個軌道參數(shù)變化項、6個調(diào)和項系數(shù)和1個參考?xì)v元的16參數(shù)模型[10]。在原來16參數(shù)的基礎(chǔ)上,GPS系統(tǒng)后續(xù)又發(fā)展了18參數(shù)模型。18參數(shù)模型中考慮衛(wèi)星所受攝動力對衛(wèi)星軌道半長軸、平均角速度及升交點赤經(jīng)的影響,相較于16參數(shù)模型更加注重衛(wèi)星軌道的瞬時變化,擬合精度更高[11]。16參數(shù)模型和18參數(shù)模型都是針對中高軌衛(wèi)星設(shè)計的,而低軌衛(wèi)星與中高軌衛(wèi)星所受攝動力有所區(qū)別,如果將這種參數(shù)模型直接應(yīng)用在低軌衛(wèi)星上,將無法滿足高精度的定位需求[12],因此需要根據(jù)低軌衛(wèi)星的特點重新設(shè)計廣播星歷參數(shù)模型。文獻(xiàn)[12]中利用軌道高度約1000km的低軌衛(wèi)星在16參數(shù)模型的基礎(chǔ)上設(shè)計了26參數(shù)的廣播星歷擬合算法,增加了軌道偏心率的變化率、升交角線性變化率和8個調(diào)和改正項振幅共10個參數(shù),其2h軌道弧段的均方差小于10m,局部弧段擬合誤差小于25m,但參數(shù)過多,因此占用的導(dǎo)航信號資源多。文獻(xiàn)[6]中考慮了低軌衛(wèi)星攝動力的短期變化快,在GLONASS廣播星歷參數(shù)模型的基礎(chǔ)上設(shè)計了基于軌道狀態(tài)的22參數(shù)模型。其研究顯示:當(dāng)?shù)蛙壭l(wèi)星軌道高度大于700km時,用戶距離誤差(FURE)精度在5cm以內(nèi);當(dāng)軌道高度降低至500km時,FURE精度優(yōu)于0.1m。雖然相對于26參數(shù)模型來說22參數(shù)模型參數(shù)少,但是用于確定衛(wèi)星位置的用戶算法外推效率較低,因此限制了其被廣泛應(yīng)用。針對26參數(shù)和22參數(shù)占用導(dǎo)航信號資源多、用戶算法復(fù)雜的問題,本文提出了一種兼具擬合精度、計算效率和占用導(dǎo)航信號資源少的模型??紤]到大氣阻力主要改變平面內(nèi)的軌道要素,在18參數(shù)模型的基礎(chǔ)上僅增加偏心率變化率這一參數(shù),設(shè)計了19參數(shù)模型,并針對500~2000km高度的低軌衛(wèi)星進(jìn)行了19參數(shù)模型最小二乘擬合,檢驗了這種星歷模型的有效性,驗證了其擬合精度要優(yōu)于18參數(shù)模型。本文提出的19參數(shù)模型對大傾角(不小于30°)衛(wèi)星的長時弧段(4h)的位置擬合精度與26參數(shù)模型近似,但所需參數(shù)量大大減少,可為低軌導(dǎo)航增強衛(wèi)星的廣播星歷提供一種兼顧精度與計算資源的模型。
18參數(shù)模型采用軌道根數(shù)加攝動參數(shù)的形式表示,包括6個開普勒軌道參數(shù)、5個軌道根數(shù)變化項、6個調(diào)和項系數(shù)和1個參考?xì)v元。19參數(shù)模型在18參數(shù)模型的基礎(chǔ)上添加偏心率變化率這一參數(shù)。表1給出了參考?xì)v元toe時刻19參數(shù)廣播星歷的符號說明。
表1 19參數(shù)廣播星歷參數(shù)定義Table 1 19-parameter broadcast ephemeris parameter definition
用戶算法是指用戶使用廣播星歷計算衛(wèi)星位置的公式。在短時間內(nèi),使用廣播星歷可以較精確地描述衛(wèi)星的位置。在之后的一段時間里,廣播星歷仍能描述衛(wèi)星的位置,只是誤差會隨著外推時間變長而增大。對于一般軌道,19參數(shù)廣播星歷用戶算法如下。
1)計算觀測歷元t到參考?xì)v元t0的時間差tk
tk=t-t0
(1)
式中:k為第k時刻。
2)計算觀測歷元半長軸A
(2)
3)計算衛(wèi)星平均角速度
(3)
式中:μ為地心引力常數(shù);Δn為觀測歷元相比參考?xì)v元的衛(wèi)星平均角速度變化量;n為觀測歷元衛(wèi)星平均角速度。
4)計算觀測瞬間的平近點角Mk
Mk=M0+ntk
(4)
5)計算偏心率e
(5)
6)迭代計算偏近點角Ek
Ek=Mk+esinEk
(6)
7)計算真近點角Vk
(7)
8)計算攝動改正前的緯度輻角Φk
Φk=Vk+ω
(8)
9)計算攝動改正后的緯度輻角uk、衛(wèi)星矢量半徑rk和軌道傾角ik的攝動改正項δuk,δrk,δik
(9)
10)計算經(jīng)過攝動改正的uk,rk,ik
(10)
11)計算衛(wèi)星在軌道平面上的位置(xk,yk)
(11)
12)計算觀測時刻的升交點經(jīng)度(格林威治子午圈到衛(wèi)星軌道升交點)
(12)
13)計算衛(wèi)星在地心固連坐標(biāo)系中的位置(X,Y,Z)
(13)
當(dāng)直接利用常規(guī)的軌道根數(shù)方法對軌道傾角近似為0°的衛(wèi)星進(jìn)行廣播星歷參數(shù)擬合時,會導(dǎo)致法方程系數(shù)矩陣成為病態(tài)矩陣,使方程組的解誤差變大,進(jìn)而造成擬合精度降低。因此,為了解決小軌道傾角造成的衛(wèi)星廣播星歷擬合精度降低甚至擬合失敗的問題,本文采用軌道旋轉(zhuǎn)法[11],先將地心慣性坐標(biāo)系繞X軸(由地心指向春分點)旋轉(zhuǎn)一個角度,在新坐標(biāo)系計算廣播星歷,在用戶預(yù)測衛(wèi)星位置時將慣性坐標(biāo)系旋轉(zhuǎn)相反的角度,得到位置矢量在原慣性坐標(biāo)系下的表達(dá),進(jìn)一步得到地心固連坐標(biāo)系下的位置矢量。對于軌道傾角近似為0°的情況,坐標(biāo)轉(zhuǎn)換的過程見圖1。
19參數(shù)廣播星歷擬合算法流程如圖2所示。首先,輸入被擬合衛(wèi)星一段時間的精密位置和初始時刻的精密速度,設(shè)置擬合時長,并對小傾角情況進(jìn)行坐標(biāo)系旋轉(zhuǎn)。然后,設(shè)計19參數(shù)的迭代初值,除軌道六要素對應(yīng)的參數(shù)外其他參數(shù)置為零。最后,利用最小二乘擬合算法迭代計算19參數(shù),直至迭代改正量足夠小時停止。其中,當(dāng)傾角不大于ε時,認(rèn)為屬于小傾角情況,需要進(jìn)行坐標(biāo)系旋轉(zhuǎn),本文取ε為1°。通過判斷改正量的模是否足夠小(可以取閾值為1×10-2)或已經(jīng)達(dá)到最大迭代次數(shù)(可以取10000次)作為迭代終止條件。下面對擬合算法進(jìn)行詳細(xì)介紹。
圖2 廣播星歷擬合算法流程Fig.2 Broadcast ephemeris fitting algorithm flow
廣播星歷參數(shù)的設(shè)計是實現(xiàn)高精度擬合的前提條件,而擬合算法的設(shè)計則是擬合精度和穩(wěn)定性的重要保證。擬合算法的核心內(nèi)容是參數(shù)估計,經(jīng)典的參數(shù)估計法有最小二乘估計、極大似然估計、極大驗后估計、貝葉斯估計等。其中,經(jīng)典的最小二乘估計法是比較簡單且容易理解的估計方法,在參數(shù)模型合理、觀測值誤差較小的情況下,采用最小二乘擬合算法就可以保證較高的擬合精度。
批處理最小二乘擬合算法的基本原理為是結(jié)合動力學(xué)模型與觀測模型對狀態(tài)量進(jìn)行迭代平差,使得給定時間段內(nèi)測量估計值與實際測量值之差的平方和最小。
對于一個非線性系統(tǒng),有
z=h(x)+ν
(14)
式中:z為低軌衛(wèi)星的精確位置信息;x為19參數(shù)廣播星歷;h(x)為通過用戶算法計算出的低軌衛(wèi)星的位置信息;ν為精確位置和計算位置的誤差。
最小二乘擬合算法是使殘差的平方和最小,即讓評價函數(shù)J最小。
J=(z-h(x))T(z-h(x))
(15)
由于最小二乘擬合算法只能用于線性系統(tǒng),因此需要對系統(tǒng)線性化,在x0處線性化后的殘差可以表示為
Δz-HΔx0
(16)
可以使用線性最小二乘擬合算法求解式(14),得到
Δx0=(HTH)-1(HTΔz)
(17)
進(jìn)而,原問題的最小二乘解為x=x0+Δx0。在線性化過程中舍棄了系統(tǒng)二階及以上的信息,所以不可避免地引入誤差,造成精度損失,即所得到的估計值并不是原問題的最優(yōu)估計值。因此,還需要以所得估計值作為新的初值,重新迭代計算,直到算法收斂,即可得到狀態(tài)量的最優(yōu)估計值。
(18)
式中:rk為第k時刻的衛(wèi)星位置;?rk/?e由文獻(xiàn)[11]給出。
數(shù)值穩(wěn)定修正是為了確保擬合算法能得到正確的結(jié)果,下面對修正方法進(jìn)行具體介紹。
(2)升交點赤經(jīng)、近地點幅角、平近點角。為保證數(shù)值穩(wěn)定性,需要在每次迭代施加改正量之后對升交點赤經(jīng)、近地點幅角、平近點角進(jìn)行歸一化,也就是讓其落在[0,2π],之后再繼續(xù)進(jìn)行下一次迭代。例如,當(dāng)Mk過大時,可能會因為數(shù)值問題無法通過迭代求出Ek。
(3)解線性方程組。最小二乘擬合算法問題需要求解線性方程組(HTH)Δx0=HTΔz,但HTH條件數(shù)很小,甚至接近病態(tài),直接求逆或者直接使用高斯消元法求解誤差較大。因此,本文使用下三角上三角(LU)分解法求解線性方程組。LU分解法可以將1個矩陣分解為1個下三角矩陣和1個上三角矩陣的乘積。如果使用杜利特爾(Doolittle)分解,得到的下三角矩陣為單位下三角矩陣;如果使用克勞特(Crout)分解,得到的上三角矩陣為單位上三角矩陣。首先,對(HTH)Δx0=HTΔz中的HTH進(jìn)行LU分解,即LU=HTH,則(HTH)Δx0=HTΔz則可以寫為LUΔx0=HTΔz,令y=UΔx0,可以先解線性方程組Ly=HTΔz,得到y(tǒng),再解線性方程組y=UΔx0,得到Δx0,從而相比于直接解線性方程組能夠降低系數(shù)矩陣的病態(tài)程度,誤差更小。
本文使用高精度軌道預(yù)報器預(yù)報的衛(wèi)星位置作為精密星歷,高精度軌道預(yù)報器中非球形引力攝動考慮70階70次全球超高階地球重力場模型(EGM2008),大氣阻力模型使用質(zhì)譜儀非相關(guān)散射模型(NRLMSISE-00),同時考慮日月三體引力;參考時刻為2010-09-19T04:00(UTC),參考時刻的偏心率取0,升交點赤經(jīng)、近地點幅角和平近點角均取0°。
下面對不同軌道高度(500km,1000km,1500km,2000km)及不同軌道傾角(0°,30°,60°,90°)的軌道進(jìn)行19參數(shù)廣播星歷最小二乘擬合。設(shè)置擬合初值如下。參考時刻半長軸取為參考時刻衛(wèi)星的半長軸,參考時刻升交點赤經(jīng)變化率取為0,19參數(shù)中的參考時刻、偏心率、升交點赤經(jīng)、近地點幅角和平近點角按照衛(wèi)星參考時刻的實際數(shù)據(jù)設(shè)置,其余參數(shù)設(shè)為0。衛(wèi)星位置采樣時間間隔為1min,即每分鐘獲取1組位置,擬合時長分別考慮20min,30min,60min,2h,3h,4h的情況。
繪制上述仿真場景下的位置均方根誤差(RMSE),計算公式為
(19)
式中:Δx,Δy,Δz為地心固連坐標(biāo)系下的位置誤差。
表2~5中給出了18參數(shù)和19參數(shù)模型擬合的RMSE曲線對各個時刻的均方根,圖3~6中給出了同一軌道傾角、同一擬合時長時不同軌道高度對擬合效果的影響。
圖3 軌道傾角0°、擬合時長60min時不同軌道高度對18參數(shù)和19參數(shù)模型的擬合效果Fig.3 Fitting results of different orbit altitudes on 18-parameter and 19-parameter models when inclination angle is 0° and fitting time is 60 minutes
表2 軌道傾角為0°時不同軌道高度擬合結(jié)果的RMSETable 2 RMSEs of fitting results for different orbit altitudes at an inclination angle of 0° m
從表2和圖3可以看出:對于任何擬合時長,擬合RMSE隨軌道高度的升高而減小,原因是軌道高度越高,同等時間下軌道的非線性越弱,最小二乘法更能發(fā)揮優(yōu)勢;同時,隨著軌道高度的升高及非球形引力攝動和大氣阻力攝動的下降,軌道更接近二體軌道,因此擬合誤差會有所下降。19參數(shù)模型擬合時長為20min時精度達(dá)到厘米級,30min時精度達(dá)到分米級,60min時精度達(dá)到米級。隨著擬合時間的增長,RMSE有所增大,這是因為軌道根數(shù)預(yù)報是基于二體模型的,即使廣播星歷增加了攝動修正項,但其只有6個系數(shù),遠(yuǎn)遠(yuǎn)不足以描述全部的攝動力,因此也只能解決短時間內(nèi)的擬合問題,隨著時間增長,誤差仍然會越來越大。總體上來說,19參數(shù)模型的擬合精度要高于18參數(shù)模型。
從表3和圖4可以看出:軌道傾角為30°時同樣可以得到與軌道傾角為0°時相同的結(jié)論,即對于任何擬合時長,擬合RMSE隨軌道高度升高而減小,整體上19參數(shù)模型的精度要高于18參數(shù)模型。對比軌道傾角為30°和0°的2種場景,可以發(fā)現(xiàn):當(dāng)軌道高度和擬合時長相等時,30°軌道的擬合誤差略大于0°軌道的擬合誤差。
圖4 軌道傾角30°、擬合時長60min時不同軌道高度對18參數(shù)和19參數(shù)模型的擬合效果Fig.4 Fitting results of different orbit altitudes on 18-parameter and 19-parameter models when inclination angle is 30° and fitting time is 60 minutes
表3 軌道傾角為30°時不同軌道高度擬合結(jié)果的RMSETable 3 RMSEs of fitting results for different orbit altitudes at an inclination angle of 30° m
從表4和圖5可以看出:軌道傾角為60°時同樣可以得到與軌道傾角為0°和30°時相同的結(jié)論。軌道傾角60°時的擬合誤差大于軌道傾角0°時的擬合誤差,與30°傾角時近似。
圖5 軌道傾角60°、擬合時長60min時不同軌道高度對18參數(shù)和19參數(shù)模型的擬合效果Fig.5 Fitting results of different orbit altitudes on 18-parameter and 19-parameter models when inclination angle is 60° and fitting time is 60 minutes
表4 軌道傾角為60°時不同軌道高度擬合結(jié)果的RMSETable 4 RMSEs of fitting results for different orbit altitudes at an inclination angle of 60° m
從表5和圖6可以看出:軌道傾角為90°時同樣可以得到與軌道傾角為0°、30°和60°時相同的結(jié)論。
圖6 軌道傾角90°、擬合時長60min時不同軌道高度對18參數(shù)和19參數(shù)模型的擬合效果Fig.6 Fitting results of different orbit altitude on 18-parameter and 19-parameter models when inclination angle is 90° and fitting time is 60 minutes
表5 軌道傾角為90°時不同軌道高度擬合結(jié)果RMSETable 5 RMSEs of fitting results for different orbit altitudes at an inclination angle of 90° m
總的來看,19參數(shù)模型對不同軌道高度(500~1000km)、不同軌道傾角的衛(wèi)星在不同擬合弧段時長下的擬合均具有很好的適用性;同時,19參數(shù)模型的精度要高于18參數(shù)模型。
進(jìn)一步對比19參數(shù)模型和26參數(shù)模型的長弧段(4h)位置擬合精度,結(jié)果如表6所示。從表6可以看出:當(dāng)軌道傾角不小于30°時,在12種不同軌道高度和傾角的情況中,有8種情況是19參數(shù)模型精度更高,另外4種情況是26參數(shù)模型精度更高。這證明了增加軌道傾角變化率這一參數(shù)對長弧段擬合精度的提升是有顯著幫助的,使得19參數(shù)模型對于大傾角衛(wèi)星的長弧段擬合達(dá)到了與26參數(shù)模型近似的精度。同時,26參數(shù)模型在有些情況下擬合精度低于18參數(shù)模型,這是因為26參數(shù)模型是由16參數(shù)模型發(fā)展而來的,其沒有考慮18參數(shù)模型相較16參數(shù)模型增加的軌道高度變化率和角速度變化率。
表6 不同軌道高度、不同傾角的4h擬合結(jié)果RMSETable 6 RMSEs of fitting results for different orbit altitudes and inclination angles when fitting time is 4 hours m
本文在GPS的18參數(shù)模型基礎(chǔ)上增加了偏心率變化率這一參數(shù),仿真結(jié)果證明這一參數(shù)對于提高低軌導(dǎo)航增強衛(wèi)星的位置擬合精度具有重要作用。如何在本文提出的19參數(shù)模型基礎(chǔ)上,在不增加過多參數(shù)前提下進(jìn)一步提高位置擬合精度,值得進(jìn)一步研究。表6為此提供了重要方向,即:對于小傾角(0°)衛(wèi)星的長弧段擬合,26參數(shù)模型的位置擬合精度遠(yuǎn)高于19參數(shù)模型;而對于大傾角(不小于30°)衛(wèi)星,二者精度則近似。因此,導(dǎo)致26參數(shù)模型對小傾角衛(wèi)星的位置擬合精度較高的主要參數(shù)十分值得研究,進(jìn)而可以將這些起主要作用的參數(shù)與19參數(shù)模型結(jié)合,以提高后者對于小傾角衛(wèi)星的位置擬合精度。
本文針對軌道高度為500~2000km的低軌衛(wèi)星,在GPS的18參數(shù)模型基礎(chǔ)上僅增加偏心率變化率這一參數(shù),設(shè)計了19參數(shù)模型,并與18參數(shù)模型的精度進(jìn)行了比較,試驗分析結(jié)論如下。
(1)對于軌道高度為500~2000km的低軌衛(wèi)星,采用19參數(shù)模型擬合精度高。
(2)隨著低軌衛(wèi)星軌道高度的降低,19參數(shù)模型的擬合精度降低,這是因為軌道高度越低,大氣阻力的影響越顯著,導(dǎo)致擬合精度降低。
(3)總體來說,在相同軌道高度、軌道傾角和擬合時長的條件下,19參數(shù)模型的擬合精度要高于GPS的18參數(shù)模型。
此外,本文通過仿真驗證了19參數(shù)模型對于大軌道傾角衛(wèi)星的長弧段位置擬合精度與26參數(shù)模型近似,但所需參數(shù)量大大減少,因此可以提高計算效率和減少存儲資源的需求。因此,本文的研究為低軌導(dǎo)航增強衛(wèi)星的廣播星歷提供了一種兼顧精度高與資源需求少的模型。