顧及GPS軌道凹凸性的標準化方法改進研究*

趙自強 劉志平

(中國礦業大學國土環境與災害監測國家測繪地理信息局重點實驗室，徐州 221116)

1 引言

IGS 提供的SP3 精密星歷產品間隔為15 分鐘，遠低于導航定位應用的采樣率(30 s、15 s 甚至更高)，故需對衛星軌道進行標準化處理［1］。常用的軌道標準化方法有插值算法［2，3］和擬合方法［1，4，5］，在相同階數的情況下，擬合方法較插值算法的標準化結果更為穩定［5］。常規擬合方法主要有Chebyshev 多項式［6］、Legendre 多項式［1］、三角函數［7］等，當采用多項式擬合方法進行軌道標準化時，過高的階數易產生“Runger現象”［8］，對此，許多學者［2，9，10］研究滑動式擬合方法，有效地提高了軌道標準化精度及可靠性。但是，滑動式擬合方法采用逐節點移動處理策略，一般只保留滑動區間中心結果，從而降低了數據處理效率。文獻［4，11］發現，衛星軌道在天球坐標系下具有較好的周期規律，進而提出了顧及衛星軌道周期特征的滑動式擬合方法。但是，該方法需進行天球坐標與地固坐標的相互轉換，軌道標準化精度易受轉換矩陣的影響。此外，最佳擬合階數與軌道弧段特征、節點數及擬合精度緊密相關［12］，但現有研究僅考慮擬合精度并得出擬合階數為9 階［5］、10 ～11 階［1］、11 ～12 階［9］、15階［13］等，缺乏合理統一的定階準則。而且，常規擬合方法中多項式系數的求解多采用等權最小二乘［1，13］，沒有利用IGS 事后精密星歷產品中的坐標中誤差信息。綜上分析，本文提出了顧及地固坐標系下軌道凹凸性特征的快速加權擬合方法，并利用實際數據驗證了該方法的有效性與正確性。

2 常規擬合方法

由于Chebyshev 多項式擬合、Legendre 多項式擬合方法的軌道標準化精度是等效的［1，5］，故以Chebyshev 多項式為例進行討論。

Chebyshev 多項式為最佳一致逼近。要在區間［t0，t0+Δt］進行n 階Chebyshev 多項式擬合，需先對時間變量t∈［t0，t0+Δt］進行變換:

式中，τ∈［－1，1］，t0為初始歷元時間，Δt 為擬合區間長度。

衛星X 坐標分量可用q 階Chebyshev 多項式表示為:

式中，Ci為Chebyshev 多項式系數，Ti(τ)有如下遞推規律:

T0(τ)=1，T1(τ)=τ，…，

Tn(τ)=2τTn－1(τ)－Tn－2(τ)

衛星Y、Z 坐標分量的求解方法與X 坐標分量相同。

3 擬合方法的改進

3.1 凹凸區間

擬合區間的選取與軌道弧段特征緊密相關。鑒此，本文根據衛星軌道運動特征，以衛星坐標二次差代數符號的變化(拐點)為特征點進行區間的劃分。首先，介紹拐點的確定方法:

衛星的X 坐標分量序列為X(i)，(i=1，2，…，n)，每次取出4 個連續歷元的坐標序列X(i－1)，X(i)，X(i+1)，X(i+2)，(i=2，…，n－2)，按式(3)求二次差:

判斷式(3)所得兩個連續二次差結果k(i－1)與k(i)的代數符號，若二者符號相反，則第i 個歷元即為拐點;否則，i=i +1 繼續判斷。衛星Y、Z 坐標分量拐點的求解方法與X 坐標分量相同。圖1 顯示了PRN2 衛星軌道拐點求解的結果。

圖1 PRN2 衛星軌道拐點Fig.1 Inflection points of PRN2 orbit

事實上，上述二次差結果近似表征了坐標分量加速度，表明衛星坐標分量在相鄰兩拐點之內的區間具有同向加速度，而該區間在數學圖形上表現為上凸、上凹兩種特征(圖1)。因此，本文將相鄰兩拐點確定的連續區間稱之為凹凸區間。考慮到“Runger 現象”主要發生在擬合區間兩端［8］，因此在凹凸區間兩端各增加一個節點構成新的區間，本文稱之為擬合區間。多項式擬合在上述整個擬合區間內進行，加密星歷的內插計算則在凹凸區間內進行。

分析圖1 可知，凹凸區間的衛星軌道凹向或凸向同一側，區間內節點具有一致的軌道運動特征(加速度同向)，因此有利于多項式建模描述。此外，凹凸區間的節點數完全由衛星坐標序列拐點確定(不同凹凸區間節點數不盡相同，X、Y 方向節點數為14 ～23，Z 方向為24 ～27)，表明凹凸區間不同于固定節點數的滑動區間，能夠較好地顧及衛星軌道運動特征。需要說明的是，本文所得的凹凸區間一般不能完全覆蓋單天所有歷元節點，在單天開始與結束歷元附近需要搭接相鄰的精密星歷產品。

3.2 加權擬合

IGS 事后精密星歷產品不僅給出了軌道坐標，也給出了相應坐標的中誤差［14］。若以加權最小二乘融合上述中誤差信息σi，則所得結果將更加合理。因此，利用IGS 產品中的坐標中誤差構造權陣P，進而獲得加權Chebyshev 多項式擬合系數C 如下:

式中，

其中，σXi表示X 坐標分量中誤差。衛星Y、Z 坐標分量的求解方法與X 坐標分量相同。

3.3 最佳階數

按照上述方法對15 分鐘間隔的衛星軌道數據進行凹凸區間劃分，擬合區間節點數記為m，各區間以8 ～m－1 階進行逐階加權擬合，將擬合中誤差小于0.001 m 時的最小階數作為相應區間的最佳擬合階數。圖2 顯示了32 顆衛星軌道區間節點數與最佳擬合階數的統計關系，表1 為PRN2 衛星一個區間的計算結果。

為保證各區間擬合結果都能滿足上述精度要求，故取圖2各處最大值，得出最佳擬合階數與區間節點總數的經驗關系式(階數曲線如圖2 所示):

圖2 擬合階數與節點數關系Fig.2 Relation between fitting order and node number

表1 PRN2 擬合誤差統計結果Tab.1 Statistical results of PRN2 fitting errors

式(5)中，q 為擬合階數，m 為擬合區間的節點總數，int表示按四舍五入原則取整。

4 應用結果與分析

4.1 外符合精度分析

為驗證本文方法的正確性和可靠性，以文獻［11］提供的2002年1月1日的15 分鐘精密星歷ECF_15MIN.200(其中PRN12、16、19、24、32 無數據)內插為5 分鐘間隔的低頻數據，并與5 分鐘間隔的精密星歷ECF_5MIN.200 比較，以內插誤差和中誤差為指標進行外符合精度分析。其中，PRN2 衛星各分量誤差結果如圖3(a)所示，所有衛星的點位中誤差如圖3(b)所示。

圖3(a)結果顯示，PRN2 號衛星X、Y、Z 方向內插結果的最大誤差均不超過3 mm，且基本服從正態分布，各方向的中誤差分別為0.6、0.5 和0.7 mm。其他衛星的內插誤差分布序列與圖3(a)類似，不再展示。圖3(b)統計了所有衛星內插結果的點位中誤差，其最大值不超過1.2 mm，平均值為0.9 mm。由此，說明本文方法所得衛星位置精度可達到1 mm水平，驗證了外符合精度的可靠性。

圖3 外符合精度統計結果Fig.3 Statistical results of external accuracy

4.2 內符合精度分析

為進一步驗證所提方法的優越性和可靠性，采用不同的Chebyshev 擬合法將ECF_15MIN.200 和ECF_5MIN.200 分別內插為1 s 采樣率的數據，并以兩者對應歷元坐標的最大互差和中誤差為指標進行差異分析。具體設計方案如下:

方案一:采用常規固定長度的Chebyshev 擬合法［1］，每次取3 小時長度的數據段，以10 階內插得到1s 采樣率的數據序列。

方案二:采用常規8 階滑動式Chebyshev 擬合法［2］，其余處理策略同方案一。

方案三:采用本文方法，其余處理策略同方案一。

上述三個方案中，PRN2 衛星坐標分量誤差的統計結果見表2，所有衛星坐標點位誤差的統計結果如圖4 所示。

表2 PRN2 坐標分量誤差統計結果比較Tab.2 Results comparison of PRN2 coordinate component errors

圖4 衛星坐標點位誤差統計結果比較Fig.4 Results comparison of satellite position errors

由表2 可知，以X 坐標分量為例，方案一、二、三的最大互差依次為0.031 2、0.004 7、0.000 7 m，呈現一種明顯的遞減規律，表明所提出方法的最大互差均較其他方法的小;以中誤差為指標，方案三比方案一的精度提高7 倍，比方案二的精度提高3 倍。從而表明，以X 坐標分量進行高頻內插時，本文方法較其他方法更具穩定性。其他坐標分量也呈現類似的統計特性。

圖4(a)中徑向代表最大互差，切向代表衛星的PRN 序列;圖4(b)中徑向代表中誤差，切向代表衛星的PRN 序列。圖4 統計結果顯示，三個方案的最大點位互差及其平均中誤差依次為:方案一分別為5.1 cm 和2.3 mm，方案二分別為8 mm 和1.5 mm，方案三的分別為3 mm 和0.5 mm。對比結果說明，對于不同星歷產品的軌道標準化誤差，固定長度擬合方法存在很大差異，且不同衛星序列間差異顯著，而滑動式擬合算法有所改善，但它們均不及本文所提的方法。由此表明，利用本文方法進行軌道標準化時，所得擬合軌道形態具有較好的一致性，較其他方法具有更高的精度和穩定性。

為說明本文方法標準化效率的優越性，比較分析上述三個方案程序計算所消耗的時間。程序基于MATLAB2009a 平臺，在主頻為2.93 GHz，內存為3.5 G 的計算機中運行，所得CPU 時間如表3 所示。

表3 三個方案程序執行的CPU 時間Tab.3 Execution CPU time of program with three methods

表3 結果顯示，三個方案中，方案三計算用時最短，方案一次之，方案二最長。其主要原因是本方法大大延拓了內插窗口長度，減少了多項式擬合系數的求解次數。盡管本文方法需預先判斷拐點，但總體上并未使結果計算時間延長。對比結果表明，本文所提方法較現有標準化方法效率更高。

5 結論

顧及地固坐標系下GPS 衛星軌道運動特征，并利用精密星歷產品的坐標中誤差信息，構建了一種基于軌道凹凸性的快速加權標準化方法。所提的加權準則與最佳階數經驗公式使得整個凹凸區間內的標準化結果更準確可靠，有效地提高了軌道標準化效率。大量計算結果表明，該方法的內外符合精度分別達到0.5 mm、1 mm，較Chebyshev 擬合方法更加穩定可靠。應指出的是，本文的凹凸區間劃分等改進研究可直接應用于其他擬合類型的軌道標準化方法。

1 馮煒，等.兩種常用GPS 星歷擬合方法的精度分析［J］.大地測量與地球動力學，2010，(1):145－149.(Feng Wei，et al.Accuracy assessment of two fitting methods for GPS precise ephemeris［J］.Journal of Geodesy and Geodynamics，2010，30(1):145－149)

2 洪櫻，歐吉坤，彭碧波.GPS 衛星精密星歷和鐘差三種內插方法的比較［J］.武漢大學學報(信息科學版)，2006，31(6):516－518，556.(Hong Ying，Ou Jikun and Peng Bibo.Three interpolation methods for precise ephemeris and clock offset of GPS satellite［J］.Geomatics and Information Science of Wuhan University，2006，31(6):516－518，556)

3 張守建，等.兩種IGS 精密星歷插值方法的比較分析［J］.大地測量與地球動力學，2007，(2):80－83.(Zhang Shoujian，et al.Comparative analysis on two methods for IGS precise ephemeris interpolation［J］.Journal of Geodesy and Geodynamics，2007，(2):80－83)

4 劉偉平，郝金明.一種新的IGS 精密星歷插值算法［J］.武漢大學學報(信息科學版)，2011，(11):1320－1323.(Liu Weiping and Hao Jinming.A new interpolation method for IGS precise ephemeris［J］.Geomatics and Information Science of Wuhan University，2011，(11):1320－1323)

5 李明峰，江國焰，張凱.IGS 精密星歷內插與擬合法精度的比較［J］.大地測量與地球動力學，2008，(2):77－80.(Li Mingfeng，Jiang Guoyan and Zhang Kai.Comparison between the accuracies with interpolating and fitting precise ephemeris［J］.Journal of Geodesy and Geodynamics，2008，(2):77－80)

6 楊學鋒，等.利用切比雪夫多項式擬合衛星軌道坐標的研究［J］.測繪通報，2008，(12):1－3.(Yang Xuefeng，et al.On fitting satellite orbit coordinate by using Chebyshev polynomial［J］.Bulletin of Surveying and Mapping，2008，(12):1－3)

7 Feng Y and Zheng Y.Efficient interpolations to GPS orbits for precise wide area applications［J］.GPS Solutions，2005，9(4):273－282.

8 吉長東，徐愛功，馮磊.GPS 精密星歷擬合與插值中龍格現象的處理方法［J］.測繪科學，2011，36(6):169－171.(Ji Changdong，Xu Aigong and Feng Lei.Treatment of Runge’s phenomenon in fitting ＆ interpolating GPS precise ephemeris［J］.Science of Surveying and Mapping，2011，36(6):169－171)

9 萬亞豪，張書畢，侯東陽.GPS 精密星歷插值法與擬合法的精度分析［J］.全球定位系統，2011，36(2):40－44.(Wan Yahao，Zhang Shubi and Hou Dongyang.Accuracy analysis on the interpolation and fitting of GPS precise ephemeris［J］.GNSS World of China，2011，36(2):40－44)

10 常亮，何秀鳳.基于移動區間的GPS 軌道標準化方法［J］.大地測量與地球動力學，2009，(1):110－113.(Chang Liang and He Xiufeng.Standardization of GPS satellite orbits based on moving intervals［J］.Journal of Geodesy and Geodynamics，2009，(1):110－113)

11 Schenewerk M.A brief review of basic GPS orbit interpolation strategies［J］.GPS Solutions，2003，6(4):265－267.

12 高偉，姜水生.分段曲線擬合與離散度加權的數據誤差處理方法［J］.中國測試技術，2005，31(6):55－56.(Gao Wei and Jiang Shuisheng.Method of deal with data error by subsection curve fitting and discrete degree weight［J］.China Measurement Technology，2005，31(6):55－56)

13 孔巧麗.用切貝雪夫多項式擬合GPS 衛星精密坐標［J］.測繪通報，2006(8):1－3.(Kong Qiaoli.Using Chebyshev polynomial to fit the precise satellite ephemeris［J］.Bulletin of Surveying and Mapping，2006，(8):1－3)

14 The extended standard product 3 orbit format(SP3-c)［EB/OL］.http://www.igscb.jpl.nasa.gov/igscb/data/format/sp3c.txt.