邢志成,王解先
(同濟(jì)大學(xué) 土木工程學(xué)院 測量與國土信息工程系,上海200092)
用戶利用GPS廣播星歷參數(shù)可以計算出GPS衛(wèi)星的位置,繼而解算用戶的位置;反過來,也可以利用GPS衛(wèi)星的位置,擬合出廣播星歷的參數(shù)[1-2]。GPS廣播星歷就是由地面控制中心在確定衛(wèi)星軌道的基礎(chǔ)上,以衛(wèi)星的位置作為觀測量擬合出廣播星歷參數(shù)再上傳給衛(wèi)星,然后,衛(wèi)星再轉(zhuǎn)發(fā)這些星歷參數(shù)給用戶[3]。
利用最小二乘原理擬合廣播星歷參數(shù)的過程中,偏導(dǎo)數(shù)矩陣的求解有解析法和數(shù)值導(dǎo)數(shù)法兩種方式。前者需要推導(dǎo)位置矢量對待估參數(shù)的偏導(dǎo)數(shù),形式復(fù)雜,若改變參數(shù)形式,偏導(dǎo)數(shù)的形式需要重新推導(dǎo);與前者相比,數(shù)值導(dǎo)數(shù)計算方法形式簡單,計算方便,便于計算機(jī)操作[4]。采用數(shù)值導(dǎo)數(shù)計算方法,利用IGS精密星歷衛(wèi)星坐標(biāo)數(shù)據(jù),根據(jù)最小二乘平差原理擬合出GPS廣播星歷的15個參數(shù);并分析了計算參數(shù)的雅可比矩陣時,增量δ的選取對擬合精度和效率的影響,得出了一些結(jié)論。
GPS廣播星歷參數(shù)以開普勒軌道根數(shù)為基礎(chǔ),加上表示軌道攝動調(diào)和系數(shù)以及隨時間的變化率,一共有15個:6個開普勒軌道根數(shù)ω,Ω0,M)和9個攝動力參數(shù)(Δn,i,Ω,Cuc,Cus,Crc,Crs,Cic,Cis).15個參數(shù)的具體意義參考文獻(xiàn)[2].
若要計算時刻t衛(wèi)星的空間坐標(biāo),根據(jù)廣播星歷參數(shù),按如下步驟計算:
1)求長半軸A
2)計算平角速度n0
3)計算從需要時刻到參考時刻的時間差tk
4)改正平角速度n
5)計算平近點角Mk
6)按下式迭代計算偏近點角Ek
7)由下式計算真近點角
8)計算緯度參數(shù)φk
9)周期改正項
10)計算改正后的緯度參數(shù)uk
11)計算改正后的向徑rk
12)計算改正后的傾角ik
13)計算衛(wèi)星在軌道平面內(nèi)的坐標(biāo)(xk,yk)
14)改正升交點的經(jīng)度
15)最后計算衛(wèi)星在地固坐標(biāo)系中的坐標(biāo)Rk=(XkYkZk)
若y=f(x1,x2,…,xi…),y到xi的導(dǎo)數(shù)的定義為[5]
當(dāng)δ足夠小時,這一導(dǎo)數(shù)可以近似寫成數(shù)值導(dǎo)數(shù)的形式
結(jié)合1.1節(jié),則一組衛(wèi)星坐標(biāo)Rk=(XkYkZk)到Xi(i=1,2,…15)的雅可比矩陣可表示為
若利用同一顆衛(wèi)星的上下各5個歷元共11組坐標(biāo),則可以得到式(4)的參數(shù)的雅可比矩陣B.
對于給定的參考時刻t,設(shè)參數(shù)為Xi(i=1,2,…15),由一個歷元的衛(wèi)星坐標(biāo)Rk=[XkYkZk]T,就可以根據(jù)式(1)建立如下的誤差方程
式(6)中的R0代表R的近似值,X0代表X的近似值,F(xiàn)(X0)=R0.
則(6)式可以寫為
式(7)中L=(R-R0)(根據(jù)最小二乘原理,利用間接平差求解,可得
通過式(5),對每個歷元的衛(wèi)星坐標(biāo)可列出3個誤差方程,若選取11個歷元,那么一共有33個方程,對于15個參數(shù),進(jìn)行間接平差迭代求解。迭代的收斂條件為:|σi+1-σi|/σi<ξ,其中,ξ是預(yù)先給定的閾值,取1.0e-05.是第i次迭代的單位權(quán)方差[2]。
使用2010年1月1日從0時45分到3時15分第02顆GPS衛(wèi)星的共11組IGS精密星歷數(shù)據(jù)(間隔15min)。15參數(shù)中六個軌道根數(shù)的初值如表2所示,其余的9個參數(shù)初值設(shè)為0.
表2 六個軌道根數(shù)初值
為了比較小量δi的取值對擬合精度和效率的影響,算例根據(jù)1.2節(jié)所述的方法,利用公式Δ=f(x1,x2,…,xi+δ…)-f(x1,x2…,xi…)得到的差值Δ所在的量級不同,采用12種方案確定增量δi的取值,分別為:Δ量級為千萬米級、百萬米級、十萬米級、萬米級、千米級;百米級、十米級、米級、分米級、厘米級、毫米級、微米級。當(dāng)然,增量的取值對計算結(jié)果的影響和計算機(jī)的性能有很大關(guān)系,本算例是在32位操作系統(tǒng)下完成的。根據(jù)這個方法可以確定增量δi的取值范圍。表3列出了方案8米級的增量δi的取值。
表3 第8種方案“米級”增量δi的值
根據(jù)不同的方案,采用1.3節(jié)所述的算法進(jìn)行計算。計算迭代的次數(shù)統(tǒng)計如表4所示。
表4 12種方案統(tǒng)計
表5給出了方案8擬合出的廣播星歷15參數(shù)與GPS廣播星歷15參數(shù)的差值。
表5 擬合參數(shù)與廣播星歷參數(shù)之差Δ
圖1給出了利用方案8擬合出來的廣播星歷參數(shù)計算11顆衛(wèi)星的坐標(biāo)與IGS精密星歷差值。
圖1 方案8米級小量計算坐標(biāo)殘差
1)通過表4和圖1可知,在最小二乘原理的基礎(chǔ)上,選取適當(dāng)?shù)男×?,利用?shù)值導(dǎo)數(shù)計算方法可以有效地擬合出GPS廣播星歷,其擬合結(jié)果GPS位置與IGS精密星歷坐標(biāo)僅在厘米范圍有差別。
2)當(dāng)小量對衛(wèi)星位置的影響為萬米到厘米時,可以擬合出廣播星歷15參數(shù),擬合中迭代計算的次數(shù)較少,擬合出的參數(shù)與GPS廣播星歷給出的參數(shù)值基本一致。當(dāng)小量超出上述范圍時,擬合精度和效率降低,甚至有可能導(dǎo)致擬合失敗。
3)采用方案4至9對擬合的效率和精度都沒有明顯的差異,說明方案4到方案9都是比較理想的預(yù)定小量的方法。
[1]郝金明.利用地面測軌資料擬合GPS廣播星歷[D].鄭州:解放軍測繪學(xué)院,1989.
[2]呂志偉,易維勇,曾志林.GPS廣播星歷參數(shù)擬合算法及其分析[J].測繪科學(xué)技術(shù)學(xué)報,2010,27(2):83-86.
[3]陳劉成,韓春好,陳金平.廣播星歷參數(shù)擬合算法研究[J].北京:測繪科學(xué),2007,32(3):12-14.
[4]王解先,徐志京.三種坐標(biāo)間轉(zhuǎn)換的雅克比矩陣數(shù)值導(dǎo)數(shù)計算方法[J].大地測量與地球動力學(xué),2004,24(4):19-23.