邢志成，王解先

（同濟大學 土木工程學院 測量與國土信息工程系，上海200092）

0 引 言

用戶利用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ù)的過程中，偏導數(shù)矩陣的求解有解析法和數(shù)值導數(shù)法兩種方式。前者需要推導位置矢量對待估參數(shù)的偏導數(shù)，形式復雜，若改變參數(shù)形式，偏導數(shù)的形式需要重新推導；與前者相比，數(shù)值導數(shù)計算方法形式簡單，計算方便，便于計算機操作［4]。采用數(shù)值導數(shù)計算方法，利用IGS精密星歷衛(wèi)星坐標數(shù)據(jù)，根據(jù)最小二乘平差原理擬合出GPS廣播星歷的15個參數(shù)；并分析了計算參數(shù)的雅可比矩陣時，增量δ的選取對擬合精度和效率的影響，得出了一些結(jié)論。

1 GPS廣播星歷參數(shù)擬合算法

GPS廣播星歷參數(shù)以開普勒軌道根數(shù)為基礎(chǔ)，加上表示軌道攝動調(diào)和系數(shù)以及隨時間的變化率，一共有15個：6個開普勒軌道根數(shù)ω，Ω0，M）和9個攝動力參數(shù)（Δn，i，Ω，Cuc，Cus，Crc，Crs，Cic，Cis）.15個參數(shù)的具體意義參考文獻［2].

1.1 由廣播星歷計算衛(wèi)星位置

若要計算時刻t衛(wèi)星的空間坐標，根據(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)的坐標（xk，yk）

14）改正升交點的經(jīng)度

15）最后計算衛(wèi)星在地固坐標系中的坐標Rk＝（XkYkZk）

1.2 數(shù)值導數(shù)的概念及雅可比矩陣的確定

若y＝f（x1，x2，…，xi…），y到xi的導數(shù)的定義為［5]

當δ足夠小時，這一導數(shù)可以近似寫成數(shù)值導數(shù)的形式

結(jié)合1.1節(jié)，則一組衛(wèi)星坐標Rk＝（XkYkZk）到Xi（i＝1，2，…15）的雅可比矩陣可表示為

若利用同一顆衛(wèi)星的上下各5個歷元共11組坐標，則可以得到式（4）的參數(shù)的雅可比矩陣B.

1.3 平差模型

對于給定的參考時刻t，設(shè)參數(shù)為Xi（i＝1，2，…15），由一個歷元的衛(wèi)星坐標Rk＝［XkYkZk]T，就可以根據(jù)式（1）建立如下的誤差方程

式（6）中的R0代表R的近似值，X0代表X的近似值，F(xiàn)（X0）＝R0.

則（6）式可以寫為

式（7）中L＝（R－R0）（根據(jù)最小二乘原理，利用間接平差求解，可得

通過式（5），對每個歷元的衛(wèi)星坐標可列出3個誤差方程，若選取11個歷元，那么一共有33個方程，對于15個參數(shù)，進行間接平差迭代求解。迭代的收斂條件為：｜σi＋1－σi｜／σi＜ξ，其中，ξ是預先給定的閾值，取1.0e－05.是第i次迭代的單位權(quán)方差［2]。

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的取值，分別為：Δ量級為千萬米級、百萬米級、十萬米級、萬米級、千米級；百米級、十米級、米級、分米級、厘米級、毫米級、微米級。當然，增量的取值對計算結(jié)果的影響和計算機的性能有很大關(guān)系，本算例是在32位操作系統(tǒng)下完成的。根據(jù)這個方法可以確定增量δi的取值范圍。表3列出了方案8米級的增量δi的取值。

表3 第8種方案“米級”增量δi的值

根據(jù)不同的方案，采用1.3節(jié)所述的算法進行計算。計算迭代的次數(shù)統(tǒng)計如表4所示。

表4 12種方案統(tǒng)計

表5給出了方案8擬合出的廣播星歷15參數(shù)與GPS廣播星歷15參數(shù)的差值。

表5 擬合參數(shù)與廣播星歷參數(shù)之差Δ

圖1給出了利用方案8擬合出來的廣播星歷參數(shù)計算11顆衛(wèi)星的坐標與IGS精密星歷差值。

圖1 方案8米級小量計算坐標殘差

3 結(jié) 論

1）通過表4和圖1可知，在最小二乘原理的基礎(chǔ)上，選取適當?shù)男×?，利用?shù)值導數(shù)計算方法可以有效地擬合出GPS廣播星歷，其擬合結(jié)果GPS位置與IGS精密星歷坐標僅在厘米范圍有差別。

2）當小量對衛(wèi)星位置的影響為萬米到厘米時，可以擬合出廣播星歷15參數(shù)，擬合中迭代計算的次數(shù)較少，擬合出的參數(shù)與GPS廣播星歷給出的參數(shù)值基本一致。當小量超出上述范圍時，擬合精度和效率降低，甚至有可能導致擬合失敗。

3）采用方案4至9對擬合的效率和精度都沒有明顯的差異，說明方案4到方案9都是比較理想的預定小量的方法。

［1]郝金明.利用地面測軌資料擬合GPS廣播星歷［D].鄭州：解放軍測繪學院，1989.

［2]呂志偉，易維勇，曾志林.GPS廣播星歷參數(shù)擬合算法及其分析［J].測繪科學技術(shù)學報，2010，27（2）：83－86.

［3]陳劉成，韓春好，陳金平.廣播星歷參數(shù)擬合算法研究［J].北京：測繪科學，2007，32（3）：12－14.

［4]王解先，徐志京.三種坐標間轉(zhuǎn)換的雅克比矩陣數(shù)值導數(shù)計算方法［J].大地測量與地球動力學，2004，24（4）：19－23.