基于移動區(qū)間的GPS精密星歷內(nèi)插方法①

張晨艷，何秀鳳，常 亮

（河海大學(xué)衛(wèi)星及空間信息應(yīng)用研究所，江蘇 南京210098）

0 引 言

GPS精密星歷由國際衛(wèi)星導(dǎo)航服務(wù)（IGS）、美國國家大地測量局（NGS）等分析中心發(fā)布，給出了衛(wèi)星在國際地球參考架（ITRF）中的三維坐標及衛(wèi)星鐘差，其中最終星歷的精度小于5cm，可以滿足高精度定位的要求［1]。GPS精密星歷的數(shù)據(jù)采樣間隔一般為15min，但在實際應(yīng)用中，GPS接收機的采樣率為30s、15s、5s，甚至更密［2－4]。因此，需要將GPS精密星歷進行內(nèi)插或外推以得到特定歷元的GPS衛(wèi)星的精密軌道。另外，幾乎所有的GPS數(shù)據(jù)處理都需要對GPS精密星歷進行內(nèi)插或外推，無論是GPS相對定位還是精密單點定位，都需要某時刻的衛(wèi)星位置才能順利建立各個歷元的觀測方程［4－5]。

國內(nèi)外比較常用的GPS精密星歷插值方法是Lagrange多項式插值、Neville多項式插值、Chebyshev多項式擬合、Trigonometric多項式插值等方法［5]。由于這些多項式插值方法隨著階數(shù)的增加，出現(xiàn)精度衰減或不穩(wěn)定的問題。對這一問題，本文提出了移動區(qū)間的概念，在精密星歷內(nèi)插的過程中，通過使被插值節(jié)點始終位于移動區(qū)間內(nèi)，提高了上述多項式插值方法的精度，插值精度更加穩(wěn)定。

采用移動區(qū)間的方法比較7家GPS分析中心所提供的精密星歷的質(zhì)量。7家提供GPS精密星歷的分析中心分別是：國際GNSS服務(wù)（IGS）、美國國家大地測量局（NGS）、加拿大自然資源部（NRCan）、歐洲航天局（ESA）、噴氣動力試驗室（JPL）、麻省理工學(xué)院（MIT）、斯科利普斯海洋研究所（SIO），它們提供的數(shù)據(jù)產(chǎn)品的名稱分別為：igs、ngs、emr、esa、jpl、mit、sio.然后采用數(shù)據(jù)質(zhì)量較好的精密星歷，利用基于移動區(qū)間的插值方法，比較隨著階數(shù)的增加、插值方法的變化插值精度的變化情況。

2 幾種GPS精密星歷插值方法

2.1 Lagrange多項式插值

設(shè)有n＋1個節(jié)點歷元t0，t1，t2…，tn，對應(yīng)的精密星歷坐標某項為：x0，x1，x2…，xn，則計算相應(yīng)時刻衛(wèi)星坐標的n階插值多項式為［7－11]：

利用公式（1）分別對衛(wèi)星坐標的三個分量（X，Y，Z）進行插值計算，即可得到觀測時刻衛(wèi)星的坐標。

2.2 Neville插值

Neville插值是一種線性逐次插值，其基本思想是通過低一次多項式的組合來獲得高一次插值多項式。首先求出若干個一次插值多項式的值，用一次多項式兩兩組合得到若干個二次插值多項式的值，利用二次多項式之間兩兩組合得到三次多項式得到值，依次采用迭代方式求出各插值節(jié)點的值。

設(shè)有n＋1個歷元時刻t0，t1，t2…，tn對應(yīng)的精密星歷坐標某項為：x0，x1，x2…，xn，令Ti，0＝xi（i＝0，1，2…n），則有［8]

2.3 三角多項式（Trigonometric）插值

三角多項式插值的函數(shù)為

式中：t為精密星歷的歷元時刻；T（t）為歷元t對應(yīng)的衛(wèi)星的坐標分量（X，Y，Z）；a0、ai、bi（i＝1，2，…，n）為待求系數(shù)；ω為特征因子。根據(jù)衛(wèi)星的運動周期，對于慣性系下的軌道，ω＝2π／0.997 269 566 34.其中，0.997 269 566 34為恒星時與太陽時的比率，而對于地心地固坐標系下的軌道縮小了一倍，即ω＝2π／2×0.997 269 566 34.三角函數(shù)多項式插值中考慮了衛(wèi)星坐標在兩種框架下各自的周期特性［1，2，10]。

3 移動區(qū)間的精密星歷插值方法

采用上述的多項式插值方法進行精密星歷的內(nèi)插時，隨著階數(shù)的增加，插值會發(fā)生衰減或不穩(wěn)定現(xiàn)象。例如，Lagrange多項式插值隨著階數(shù)的增加會出現(xiàn)龍格（Runge）振蕩現(xiàn)象，插值結(jié)果在端點處誤差比較大，基于此提出了基于移動區(qū)間的GPS精密星歷插值方法。

移動區(qū)間的定義：以一個指定大小的區(qū)間以固定的步長在二維空間中移動，移動時，即可采用某種多項式插值計算每個區(qū)間內(nèi)的待插值節(jié)點對應(yīng)的坐標值，此區(qū)間即稱為移動區(qū)間。在GPS精密星歷內(nèi)插的過程中，當(dāng)所有插值節(jié)點按遞增或遞減排列時，移動區(qū)間的選取為位于插值節(jié)點序列中間的兩個插值節(jié)點構(gòu)成的區(qū)間，也就是說插值節(jié)點始終位于插值范圍的中心［4]。例如我們想要得到6階多項式，每次插值需要7個節(jié)點完成插值，而被插值點始終是位于第三個到第四個點之間。第一個點到第七個點就是第一個插值區(qū)間，第三個節(jié)點到第四個節(jié)點稱為移動區(qū)間，僅用來插值位于移動區(qū)間之間的被插值節(jié)點。插值區(qū)間和移動區(qū)間在該次插值過程中區(qū)間的大小始終保持不變，每次將窗口向后移動一個歷元的距離，僅用來插值位于移動區(qū)間兩個歷元之間的被插值點［6]。

基于移動區(qū)間的精密星歷內(nèi)插的實現(xiàn)步驟如下所示：

1）移動區(qū)間的起始節(jié)點t＿start和終止節(jié)點t＿end的確定

①Trigonometric插值起始節(jié)點和終止節(jié)點的

確定：m階插值需要2m＋2個插值節(jié)點t1，t2，…，t2m＋1，t2m＋2，所對應(yīng)的坐標值為 X1，X2，…，X2m＋1，X2m＋2，設(shè)｛Xi｝單調(diào)遞增（i＝1，…，2m＋1，2m＋2）。起始節(jié)點和終止節(jié)點分別為：［t＿start，t＿end]＝［tm＋1，tm＋2]，［t1，t2m＋2]為插值區(qū)間。

②Lagrange插值和Neville插值起始節(jié)點和終止節(jié)點的確定：m階插值需要m＋1個插值節(jié)點t1，t2，…，tm＋1，所 對 應(yīng) 的 坐 標 值 為 X1，X2，…，Xm＋1，設(shè) ｛Xi｝單調(diào)遞增（i＝1，…，m，m＋1）.

當(dāng)m＝2n＋1時，［t＿start，t＿end]＝［tn＋1，tn＋2]，n＝1，2，…，插值區(qū)間為［t1，tm].

當(dāng)m＝2n時，［t＿start，t＿end]＝［tn，tn＋1]，n＝1，2，…，插值區(qū)間為［t1，tm].

2）選取插值區(qū)間節(jié)點［t1，t2m＋2]或［t1，tm]處的坐標值為［X1，X2m＋2]或［X1，Xm].

3）采用移動區(qū)間的方法進行精密星歷的插值計算。

① Trigonometric插值：插值區(qū)間為［t1，t2m＋2]，用 Trigonometric插值求出［tm＋1.tm＋2]內(nèi)各節(jié)點所對應(yīng)的坐標值。

②Lagrange插值、Neville插值：當(dāng)m＝2n＋1時，插值區(qū)間為［t1，tm]，可用 Lagrange插值或Neville插值求出［tn＋1，tn＋2]內(nèi)各節(jié)點所對應(yīng)的坐標值；當(dāng) m＝2n 時，插值區(qū)間［t1，tm]，可 用Lagrange插值或 Neville插值求出［tn，tn＋1]內(nèi)各節(jié)點所對應(yīng)的坐標值。

4 實驗結(jié)果與分析

以30min間隔的衛(wèi)星位置插值15min間隔的衛(wèi)星位置。采用2011年1月29日PRN2的15 min間隔的精密星歷，精密星歷數(shù)據(jù)是由以下分析中心提供的：IGS、NGS、ESA、NRCan、ESA、JPL、MIT、SIO。利用相關(guān)程序提取各分析中心30min間隔的衛(wèi)星位置作為已知的插值節(jié)點，將各分析中心的衛(wèi)星位置求平均值作為15min間隔的衛(wèi)星位置的真值，采用基于移動區(qū)間的插值方法進行插值計算，比較出各分析中心的數(shù)據(jù)質(zhì)量，采用數(shù)據(jù)質(zhì)量比較好的分析中心提供的數(shù)據(jù)進行插值計算，分析隨著插值階數(shù)的增加，插值方法的變化，插值精度的變化情況。

4.1 不同數(shù)據(jù)產(chǎn)品的數(shù)據(jù)質(zhì)量比較

分別將7家分析中心2011年1月29日PRN2的30min間隔的數(shù)據(jù)采用基于移動區(qū)間的Lagrange多項式插值、Neville多項式插值、Trigonometric多項式插值方法進行插值計算，插值到11階。將插值結(jié)果與各分析中心的平均值即真值之間進行求差，得到的結(jié)果如下圖所示。

圖1、圖2和圖3中各分析中心的插值結(jié)果與真值的坐標差是基于移動區(qū)間的Lagrange多項式、Neville多項式、Trigonometric多項式在11階時與真值的坐標差。在圖1、圖2、圖3中畫兩條額外的深虛線，分析各圖可知：由IGS提供的精密星歷數(shù)據(jù)經(jīng)過內(nèi)插之后與真值之差全部在±2cm之內(nèi)，而其他分析中心的數(shù)據(jù)產(chǎn)品的插值結(jié)果與真值之差有部分數(shù)據(jù)落在虛線外。因此，IGS分析中心提供的精密星歷質(zhì)量比較好，在將來的計算中，需要特定時間間隔的衛(wèi)星位置時，首選IGS提供的精密星歷數(shù)據(jù)進行插值計算。

圖1 不同數(shù)據(jù)產(chǎn)品基于移動區(qū)間的Lagrange插值結(jié)果與真值之差的比較圖

4.2 各插值方法精度比較

經(jīng)過上面的實驗分析，IGS分析中心提供的精密星歷數(shù)據(jù)質(zhì)量比較好。將IGS提供的2011年1月29日PRN2的30min間隔的精密星歷分別采用基于移動區(qū)間的Lagrange插值、Neville插值、Trigonometric插值內(nèi)插到15min，并與IGS提供的15min間隔的數(shù)據(jù)（稱為真值）進行比較。計算結(jié)果如表1所示。

進行分析可知：

1）Lagrange插值方法和Neville插值方法的插值精度相同。

2）采用基于移動區(qū)間的插值方法收斂比較迅速，并且隨著階數(shù)的增加，插值精度比較穩(wěn)定，在特定的階數(shù)下Neville多項式插值和Trigonometric多項式插值方法均能達到亞厘米級的精度。

3）基于移動區(qū)間的Neville多項式插值在11階時可以滿足精度要求，并且在11階以后收斂，可以達到亞厘米級的精度，而Trigonometric收斂速度更快，在4階時就能滿足精度要求，并且在4階以后收斂，同樣，也能達到亞厘米級的精度要求且精度保持穩(wěn)定。

4）基于移動區(qū)間的Neville插值方法的收斂速度較Trigonometric插值方法慢，而Trigonometric插值方法理論知識比Neville插值方法復(fù)雜，但是Trigonometric插值方法考慮了衛(wèi)星坐標分量的周期性，本身精度就很高，加上移動區(qū)間的方法，使其精度更加穩(wěn)定，因此，Trigonometric多項式插值方法是實現(xiàn)衛(wèi)星位置內(nèi)插計算的一種行之有效的方法。

對表1為了進一步比較三種方法的插值精度，將多項式插值到12階的結(jié)果與真值之間求差，由于Lagrange與Neville插值精度相同，本文只比較Neville與Trigonometric的插值精度。結(jié)果如圖4所示。

表1 基于移動區(qū)間的Lagrange插值、Neville插值和Trigonometric插值結(jié)果

圖4 兩種插值方法結(jié)果比較

圖4為Neville多項式、Trigonometric多項式插值到12階時與真值之差。分析圖4可知：

1）對于Neville插值方法X和Y方向的誤差在±1cm，Z方向的誤差在±5cm之內(nèi)，對于Trigonometric插值方法，X和Y方向的誤差在±5mm，Z方向的誤差在±2mm之內(nèi)。因此，基于移動區(qū)間的Trigonometric多項式插值精度比基于移動區(qū)間的Neville多項式插值精度高。

2）圖中畫出的是兩種插值方法插值到12階時的情形，結(jié)合表1可以得出，當(dāng)插值階數(shù)上升到一定的程度，兩種插值方法的精度相差不大。但就收斂速度而言，基于移動區(qū)間的Trigonometric多項式插值方法的收斂速度比其他兩種插值方法都要快，并且插值精度整體上優(yōu)于Neville多項式插值和Lagrange多項式插值的插值精度。

5 結(jié) 論

1）在精密星歷內(nèi)插的過程中，單獨使用Lagrange、Neville、Trigonometric等插值方法，隨著階數(shù)的增加，插值精度會發(fā)生衰減或不穩(wěn)定的現(xiàn)象，本文提出的基于移動區(qū)間的插值方法有效的解決了這一問題，插值精度較高且插值結(jié)果穩(wěn)定。

2）利用提出的基于移動區(qū)間的插值方法對7家分析中心的數(shù)據(jù)進行插值計算，IGS提供的精密星歷數(shù)據(jù)質(zhì)量比較好。

3）用IGS提供的30min間隔的數(shù)據(jù)進行插值15min間隔的數(shù)據(jù)采用基于移動區(qū)間的12階的Lagrange插值、12階Neville插值或基于移動區(qū)間的6階Trigonometric插值，這些方法都可以達到亞厘米級的精度，滿足精密定位的要求。在今后的應(yīng)用中，15min插值5min、5min插值30s等均可采用本文提出的移動區(qū)間方法進行精密星歷的插值計算。

［1]FENG Yan－ming，ZHENG Yi.Efficient interpolations to GPS orbits for precise wide area applications［J].GPS Solutions，2005（9）：273－282.

［2]張守建，李建成，邢樂林，等.兩種IGS精密星歷插值方法的比較分析［J].大地測量與地球動力學(xué)，2007，27（2）：80－83.

［3]彭澤泉.GPS精密星歷擬合方法的研究［J].測繪科學(xué)，2010，35（supp1）：63－65.

［4]常 亮，何秀鳳.基于移動區(qū)間的GPS軌道標準化方法［J].大地測量與地球動力學(xué)，2009，29（1）：110－113.

［5]SCHENEWERK M.A brief review of basic GPS orbit interpolation strategies［J].GPS Solutions，2003（6）：265－267.

［6]洪 櫻，歐吉坤，彭碧波.GPS衛(wèi)星精密星歷和鐘差三種內(nèi)插方法的比較［J].武漢大學(xué)學(xué)報·信息科學(xué)版，2006，31（6）：516－518.

［7]李慶楊，王能超，易大義.數(shù)值分析［M].北京：清華大學(xué)出版社，2001.

［8]鄧建中，劉之行.計算方法［M].西安：西安交通大學(xué)出版社，2001.

［9]邱 蕾，廖遠琴，花向紅.基于IGS精密星歷的衛(wèi)星坐標插值［J].測繪工程，2008，17（4）：15－18.

［10]唐華林，鄧翠美，匡翠林.一種基于三角函數(shù)多項式的GPS軌道標準化方法［J].測繪與空間地理信息，2007，30（4）：85－88.

［11]HU Hui，F(xiàn)ANG Lian.Interpolation and fitting algorithms for GPS satellite orbit［C]∥Proceeding of SPIE，2009.