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

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

一、引 言

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

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

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

二、插值多項式選取原則

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

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

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

三、基本多項式插值方法

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

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

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

(1)

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

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

(2)

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

(3)

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

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

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

四、實例分析

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

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

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

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

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

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

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

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

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

圖5 擬合區(qū)間和有效區(qū)間選取準則

圖6 一天有效區(qū)間插值結果與真值數(shù)據(jù)差值

MaxRMSxyzxyzPRN012.55102.15321.98610.48390.48880.4776PRN151.61571.84242.36110.43260.44210.5290PRN311.73111.73941.86420.41360.44300.4534

五、結束語

本文對插值多項式的選取原則進行了探討，根據(jù)選取原則設計了一種簡單的基本插值多項式，然后利用基本多項式擬合得到插值多項式系數(shù)，并對精密星歷進行插值測試。試驗結果表明，在插值邊緣部分出現(xiàn)明顯的龍格現(xiàn)象， 用基本多項式插值精密星歷時，4 h擬合區(qū)間、2 h有效區(qū)間、多項式擬合階數(shù)選取16能夠有效消除龍格現(xiàn)象的影響。根據(jù)以上選取準則，對試驗數(shù)據(jù)進行一天的插值，給出有效區(qū)間之內的插值結果，插值誤差小于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] 劉迎春，林寶軍，張曉坤. 一種衛(wèi)星精密星歷的插值方法[J]. 飛行器測控學報，2004, 23(4): 43-46.

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

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

[6] 孔巧麗. 用切貝雪夫多項式擬合GPS衛(wèi)星精密坐標[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] 關治，陸金甫. 數(shù)值方法[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] 金凱德，切尼，王國榮. 數(shù)值分析[M].3版.北京: 機械工業(yè)出版社, 2005：221-224.

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