邢志成，王解先

（同濟(jì)大學(xué) 土木工程學(xué)院 測(cè)量與國(guó)土信息工程系，上海200092）

0 引 言

用戶利用GPS廣播星歷參數(shù)可以計(jì)算出GPS衛(wèi)星的位置，繼而解算用戶的位置；反過來，也可以利用GPS衛(wèi)星的位置，擬合出廣播星歷的參數(shù)［1－2]。GPS廣播星歷就是由地面控制中心在確定衛(wèi)星軌道的基礎(chǔ)上，以衛(wèi)星的位置作為觀測(cè)量擬合出廣播星歷參數(shù)再上傳給衛(wèi)星，然后，衛(wèi)星再轉(zhuǎn)發(fā)這些星歷參數(shù)給用戶［3]。

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

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

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

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

若要計(jì)算時(shí)刻t衛(wèi)星的空間坐標(biāo)，根據(jù)廣播星歷參數(shù)，按如下步驟計(jì)算：

1）求長(zhǎng)半軸A

2）計(jì)算平角速度n0

3）計(jì)算從需要時(shí)刻到參考時(shí)刻的時(shí)間差tk

4）改正平角速度n

5）計(jì)算平近點(diǎn)角Mk

6）按下式迭代計(jì)算偏近點(diǎn)角Ek

7）由下式計(jì)算真近點(diǎn)角

8）計(jì)算緯度參數(shù)φk

9）周期改正項(xiàng)

10）計(jì)算改正后的緯度參數(shù)uk

11）計(jì)算改正后的向徑rk

12）計(jì)算改正后的傾角ik

13）計(jì)算衛(wèi)星在軌道平面內(nèi)的坐標(biāo)（xk，yk）

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

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

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

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

當(dāng)δ足夠小時(shí)，這一導(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個(gè)歷元共11組坐標(biāo)，則可以得到式（4）的參數(shù)的雅可比矩陣B.

1.3 平差模型

對(duì)于給定的參考時(shí)刻t，設(shè)參數(shù)為Xi（i＝1，2，…15），由一個(gè)歷元的衛(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），對(duì)每個(gè)歷元的衛(wèi)星坐標(biāo)可列出3個(gè)誤差方程，若選取11個(gè)歷元，那么一共有33個(gè)方程，對(duì)于15個(gè)參數(shù)，進(jìn)行間接平差迭代求解。迭代的收斂條件為：｜σi＋1－σi｜／σi＜ξ，其中，ξ是預(yù)先給定的閾值，取1.0e－05.是第i次迭代的單位權(quán)方差［2]。

2 算例分析

使用2010年1月1日從0時(shí)45分到3時(shí)15分第02顆GPS衛(wèi)星的共11組IGS精密星歷數(shù)據(jù)（間隔15min）。15參數(shù)中六個(gè)軌道根數(shù)的初值如表2所示，其余的9個(gè)參數(shù)初值設(shè)為0.

表2 六個(gè)軌道根數(shù)初值

為了比較小量δi的取值對(duì)擬合精度和效率的影響，算例根據(jù)1.2節(jié)所述的方法，利用公式Δ＝f（x1，x2，…，xi＋δ…）－f（x1，x2…，xi…）得到的差值Δ所在的量級(jí)不同，采用12種方案確定增量δi的取值，分別為：Δ量級(jí)為千萬(wàn)米級(jí)、百萬(wàn)米級(jí)、十萬(wàn)米級(jí)、萬(wàn)米級(jí)、千米級(jí)；百米級(jí)、十米級(jí)、米級(jí)、分米級(jí)、厘米級(jí)、毫米級(jí)、微米級(jí)。當(dāng)然，增量的取值對(duì)計(jì)算結(jié)果的影響和計(jì)算機(jī)的性能有很大關(guān)系，本算例是在32位操作系統(tǒng)下完成的。根據(jù)這個(gè)方法可以確定增量δi的取值范圍。表3列出了方案8米級(jí)的增量δi的取值。

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

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

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

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

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

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

圖1 方案8米級(jí)小量計(jì)算坐標(biāo)殘差

3 結(jié) 論

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

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

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

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

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

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

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