趙自強(qiáng) 劉志平

(中國礦業(yè)大學(xué)國土環(huán)境與災(zāi)害監(jiān)測國家測繪地理信息局重點實驗室，徐州 221116)

1 引言

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

2 常規(guī)擬合方法

由于Chebyshev 多項式擬合、Legendre 多項式擬合方法的軌道標(biāo)準(zhǔn)化精度是等效的［1，5］，故以Chebyshev 多項式為例進(jìn)行討論。

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

式中，τ∈［－1，1］，t0為初始?xì)v元時間，Δt 為擬合區(qū)間長度。

衛(wèi)星X 坐標(biāo)分量可用q 階Chebyshev 多項式表示為:

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

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

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

衛(wèi)星Y、Z 坐標(biāo)分量的求解方法與X 坐標(biāo)分量相同。

3 擬合方法的改進(jìn)

3.1 凹凸區(qū)間

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

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

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

圖1 PRN2 衛(wèi)星軌道拐點Fig.1 Inflection points of PRN2 orbit

事實上，上述二次差結(jié)果近似表征了坐標(biāo)分量加速度，表明衛(wèi)星坐標(biāo)分量在相鄰兩拐點之內(nèi)的區(qū)間具有同向加速度，而該區(qū)間在數(shù)學(xué)圖形上表現(xiàn)為上凸、上凹兩種特征(圖1)。因此，本文將相鄰兩拐點確定的連續(xù)區(qū)間稱之為凹凸區(qū)間?？紤]到“Runger 現(xiàn)象”主要發(fā)生在擬合區(qū)間兩端［8］，因此在凹凸區(qū)間兩端各增加一個節(jié)點構(gòu)成新的區(qū)間，本文稱之為擬合區(qū)間。多項式擬合在上述整個擬合區(qū)間內(nèi)進(jìn)行，加密星歷的內(nèi)插計算則在凹凸區(qū)間內(nèi)進(jìn)行。

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

3.2 加權(quán)擬合

IGS 事后精密星歷產(chǎn)品不僅給出了軌道坐標(biāo)，也給出了相應(yīng)坐標(biāo)的中誤差［14］。若以加權(quán)最小二乘融合上述中誤差信息σi，則所得結(jié)果將更加合理。因此，利用IGS 產(chǎn)品中的坐標(biāo)中誤差構(gòu)造權(quán)陣P，進(jìn)而獲得加權(quán)Chebyshev 多項式擬合系數(shù)C 如下:

式中，

其中，σXi表示X 坐標(biāo)分量中誤差。衛(wèi)星Y、Z 坐標(biāo)分量的求解方法與X 坐標(biāo)分量相同。

3.3 最佳階數(shù)

按照上述方法對15 分鐘間隔的衛(wèi)星軌道數(shù)據(jù)進(jìn)行凹凸區(qū)間劃分，擬合區(qū)間節(jié)點數(shù)記為m，各區(qū)間以8 ～m－1 階進(jìn)行逐階加權(quán)擬合，將擬合中誤差小于0.001 m 時的最小階數(shù)作為相應(yīng)區(qū)間的最佳擬合階數(shù)。圖2 顯示了32 顆衛(wèi)星軌道區(qū)間節(jié)點數(shù)與最佳擬合階數(shù)的統(tǒng)計關(guān)系，表1 為PRN2 衛(wèi)星一個區(qū)間的計算結(jié)果。

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

圖2 擬合階數(shù)與節(jié)點數(shù)關(guān)系Fig.2 Relation between fitting order and node number

表1 PRN2 擬合誤差統(tǒng)計結(jié)果Tab.1 Statistical results of PRN2 fitting errors

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

4 應(yīng)用結(jié)果與分析

4.1 外符合精度分析

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

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

圖3 外符合精度統(tǒng)計結(jié)果Fig.3 Statistical results of external accuracy

4.2 內(nèi)符合精度分析

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

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

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

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

上述三個方案中，PRN2 衛(wèi)星坐標(biāo)分量誤差的統(tǒng)計結(jié)果見表2，所有衛(wèi)星坐標(biāo)點位誤差的統(tǒng)計結(jié)果如圖4 所示。

表2 PRN2 坐標(biāo)分量誤差統(tǒng)計結(jié)果比較Tab.2 Results comparison of PRN2 coordinate component errors

圖4 衛(wèi)星坐標(biāo)點位誤差統(tǒng)計結(jié)果比較Fig.4 Results comparison of satellite position errors

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

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

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

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

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

5 結(jié)論

顧及地固坐標(biāo)系下GPS 衛(wèi)星軌道運動特征，并利用精密星歷產(chǎn)品的坐標(biāo)中誤差信息，構(gòu)建了一種基于軌道凹凸性的快速加權(quán)標(biāo)準(zhǔn)化方法。所提的加權(quán)準(zhǔn)則與最佳階數(shù)經(jīng)驗公式使得整個凹凸區(qū)間內(nèi)的標(biāo)準(zhǔn)化結(jié)果更準(zhǔn)確可靠，有效地提高了軌道標(biāo)準(zhǔn)化效率。大量計算結(jié)果表明，該方法的內(nèi)外符合精度分別達(dá)到0.5 mm、1 mm，較Chebyshev 擬合方法更加穩(wěn)定可靠。應(yīng)指出的是，本文的凹凸區(qū)間劃分等改進(jìn)研究可直接應(yīng)用于其他擬合類型的軌道標(biāo)準(zhǔn)化方法。

1 馮煒，等.兩種常用GPS 星歷擬合方法的精度分析［J］.大地測量與地球動力學(xué)，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 衛(wèi)星精密星歷和鐘差三種內(nèi)插方法的比較［J］.武漢大學(xué)學(xué)報(信息科學(xué)版)，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］.大地測量與地球動力學(xué)，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］.武漢大學(xué)學(xué)報(信息科學(xué)版)，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 精密星歷內(nèi)插與擬合法精度的比較［J］.大地測量與地球動力學(xué)，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 楊學(xué)鋒，等.利用切比雪夫多項式擬合衛(wèi)星軌道坐標(biāo)的研究［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 精密星歷擬合與插值中龍格現(xiàn)象的處理方法［J］.測繪科學(xué)，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］.全球定位系統(tǒng)，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 常亮，何秀鳳.基于移動區(qū)間的GPS 軌道標(biāo)準(zhǔn)化方法［J］.大地測量與地球動力學(xué)，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 高偉，姜水生.分段曲線擬合與離散度加權(quán)的數(shù)據(jù)誤差處理方法［J］.中國測試技術(shù)，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 衛(wèi)星精密坐標(biāo)［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.