基于基本多項式的GPS精密星歷插值方法研究

張 志，廖 瑛，文援蘭，楊雅君，李 志,2

(1. 國防科技大學 航天與材料工程學院，湖南 長沙 410073； 2. 中國空間技術研究院，北京 100094)

一、引 言

在GPS精密定位中，尤其是高精度、動態定位中，經常要用到IGS[1](International GNSS Service)提供的精密星歷，其中三維坐標的最終精度優于2.5 cm，但GPS精密星歷的數據采樣間隔為15 min。而GPS接收機的采樣率一般為30 s或15 s，甚至更密，因此，要想利用某一時刻的衛星位置，就必須對精密星歷進行高精度的插值。在GPS數據后處理中，用戶需要選擇合理的數學模型對GPS精密星歷進行插值或擬合，才能獲得任意時刻的衛星坐標，因而研究如何利用精密星歷快速地插值出觀測瞬間衛星的位置成為實際應用中一項十分必要的工作。

對GPS精密星歷的內插法主要有多項式內插法和三角內插法，采用多項式內插法和三角內插法得到的精度基本相當[2]。對于多項式內插法，要想獲得理想的精度就要合理地選擇多項式的階數和插值的時間間隔。目前，國內用于精密星歷插值較常用的方法有Lagrange插值法、Neville插值法[3-4]及Chebyshev多項式插值法[5-6]等，但對基本多項式插值法的研究很少。Mark Schenewerk[2]和Milan Horemu?[7]均提到了基本多項式插值法用于精密星歷的插值計算，前者只給出了概括性的分析說明，后者雖對多項式插值法在精密星歷的插值過程中出現的問題進行了分析，但相關文獻并未給出具體的插值計算方法。

本文利用Householder變換求解基本多項式擬合系數來避免對病態的法方程系數矩陣直接求逆計算二乘解，并對擬合曲線進行插值計算，為避免高階基本多項式插值過程中出現的龍格現象的影響，給出了多項式擬合區間長度、有效區間長度和基本插值多項式階數的選取準則，最后對插值的結果進行比較分析。

二、插值多項式選取原則

對于插值多項式(插值函數)來說，為了使插值函數的數值計算具有可操作性，通常取插值函數類為有限維子空間。為使插值便于實現，插值函數類子空間的一組基函數應當容易計算。一旦確定了插值函數類子空間的一組基函數，插值問題就可以轉化為求一組組合系數的問題。

由定理說明組合系數的唯一性：設φi,i=1,2,…，n+1是線性無關的多項式，那么任一多項式p可唯一地由φ1t,φ2t,…,φn+1t的線性組合來表示[8]。

設線性無關的一組基函數1,t,t2,…,tn，則插值多項式(插值函數)pt可唯一地由基函數線性組合表示pt=c1tn+c2tn-1+…+cnt+cn+1，此插值多項式稱為基本多項式。同樣，Chebyshev插值多項式和Legendre插值多項式等均可以唯一地由一組線性無關的基函數表示。以下將以基本多項式為例進行分析討論。

三、基本多項式插值方法

給定m個歷元時刻t1,t2,…,tm的觀測點坐標，假設需要計算n階基本多項式系數。為了改進數值計算精度，需要將變量t正則化為

式中，t∈t1,t2,…,tm；mean()、std()分別表示向量的均值和標準差。

則GPS衛星坐標(x,y,z)的基本多項式分別為

(1)

式中，n為基本多項式的階數；cxi、cyi、czi分別為(x,y,z)坐標分量的基本多項式系數，為待求量。

誤差方程的矩陣形式表示為

(2)

為最小。此處給出最小二乘問題的法方程組

(3)

式中，R1是為n+1階上三角矩陣；R中后面m-(n+1)行元素為0。

利用正交矩陣不改變向量長度的性質，又有

這樣用QR分解的最小二乘法求解出基本多項式系數cxi，同理，可求得cyi、czi，利用這些系數，根據式(1)便可計算出t∈t1,tm時間區間內任意時刻的衛星坐標。

四、實例分析

為更好地進行試驗比對，采用Mark Sche-newerk[2]中的數據，包括時間間隔(步長)為15 min的精密星歷ECF_15MI.200，文中稱為試驗數據；作為插值比較的5 min間隔數據文件ECF_5MIN.200，文中稱為真值數據，以上數據均可以從http:∥www.ngs.noaa.gov/gps-toolbox 下載。以試驗數據的第一個歷元點2002年1月1日00∶00∶00.000作為擬合起始歷元，分別采用2 h、3 h、4 h、5 h弧段進行基本多項式擬合，然后對擬合的曲線進行5 min間隔的插值，得到的插值數據與對應歷元真值數據之差稱為插值誤差，以此來說明方法的有效性。插值誤差如圖1—圖4所示。本文給出了PRN01星的比較結果，其他各星結果類似。

圖1 2 h擬合區間，8階基本多項式插值誤差

圖2 3 h擬合區間，12階基本多項式插值誤差

圖3 4 h擬合區間，16階基本多項式插值誤差

圖4 5 h擬合區間，20階基本多項式插值誤差

由圖1—圖4可以看出，高階基本多項式插值結果并不穩定，插值結果存在明顯的龍格現象[11]，即插值的兩邊緣部分的精度較低，且存在一定的震蕩性。但圖2—圖4清晰地表明，基本多項式擬合區間的中間2 h區域插值精度為毫米級，因此可以稱擬合區間中間2 h區域為有效區間，這樣，只要每次插值計算都進入有效區間，則能有效避免龍格現象的影響，從而使擬合達到毫米級的精度。對2 h的有效區間，最大插值誤差不超過5 mm，這對于厘米級的精密星歷來說，是完全可以接受的。

綜合考慮計算量和擬合階數，本文選取階數為16，擬合區間為4 h的多項式，在2 h有效區間進行插值計算，如圖5所示，對試驗數據進行4 h擬合，在一天內有效區間插值數據與對應歷元真值數據的差值如圖6所示。

由圖6可以看出，試驗數據在時間段00∶00∶00.000—23∶45∶00.000內的插值計算結果表明，有效區間插值數據只能得到1～21 h的數據結果，若希望得到試驗數據中24 h的有效插值結果，則需要提供前一天和后一天的精密星歷數據。試驗結果顯示，插值數據精度可達到毫米級，誤差小于3 mm。由于精密單點定位對衛星星歷誤差的要求為毫米級，因此這并不影響多項式擬合法插值出的數據在實際中的應用。用兩個統計量Max和RMS來反映試驗數據所有衛星基本多項式插值結果，表1列出了PRN01、PRN15和PRN31號衛星的分析結果。其中，Max指插值結果三維位置誤差絕對值的最大值；RMS指內插結果三維位置誤差的均方根誤差。

由表1可知，采用2 h的有效區間基本多項式插值方法，試驗數據3顆衛星x、y、z3個方向的插值誤差最大值均不超過3 mm，3個方向的均方根小于0.6 mm。對于厘米級的精密星歷來說，試驗結果是比較理想的。

圖5 擬合區間和有效區間選取準則

圖6 一天有效區間插值結果與真值數據差值

MaxRMSxyzxyzPRN012.55102.15321.98610.48390.48880.4776PRN151.61571.84242.36110.43260.44210.5290PRN311.73111.73941.86420.41360.44300.4534

五、結束語

本文對插值多項式的選取原則進行了探討，根據選取原則設計了一種簡單的基本插值多項式，然后利用基本多項式擬合得到插值多項式系數，并對精密星歷進行插值測試。試驗結果表明，在插值邊緣部分出現明顯的龍格現象， 用基本多項式插值精密星歷時，4 h擬合區間、2 h有效區間、多項式擬合階數選取16能夠有效消除龍格現象的影響。根據以上選取準則，對試驗數據進行一天的插值，給出有效區間之內的插值結果，插值誤差小于3 mm，能夠滿足用戶定位精度要求。

參考文獻：

[1] GERHARD B, MOORE A W, MUELLER I I. The International Global Navigation Satellite Systems Service (IGS):Development and Achievements[J]. Journal of Geodesy，2009, 83(3-4): 297-307.

[2] SCHENEWERK M. A Brief Review of Basic GPS Orbit Interpolation Strategies[J]. GPS Solutions，2003, 6(4): 265-267.

[3] 劉迎春，林寶軍，張曉坤. 一種衛星精密星歷的插值方法[J]. 飛行器測控學報，2004, 23(4): 43-46.

[4] 王曉明，成英燕，劉立. 基于階次組合的GPS精密星歷插值研究[J]. 大地測量與地球動力學，2011, 31(4): 103-106.

[5] 彭澤泉. GPS精密星歷擬合方法的研究[J]. 測繪科學，2010, 35(S1): 63-65.

[6] 孔巧麗. 用切貝雪夫多項式擬合GPS衛星精密坐標[J]. 測繪通報，2006(8): 1-3.

[7] HOREMU? M,ANDERSSON J V. Polynomial Interpolation of GPS Satellite Coordinates[J]. GPS Solutions，2006, 10(1): 67-72.

[8] 關治，陸金甫. 數值方法[M]. 北京: 清華大學出版社, 2006： 192-193.

[9] CRASSIDIS J L,JUNKINS J L. Optimal Estimation of Dynamic Systems[M].Boca Raton,Florida: CRC Press, 2004：7-13.

[10] 金凱德，切尼，王國榮. 數值分析[M].3版.北京: 機械工業出版社, 2005：221-224.

[11] 吉長東，徐愛功，馮磊. GPS精密星歷擬合與插值中龍格現象的處理方法[J]. 測繪科學，2011, 36(6): 169-171.