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

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

一、引 言

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

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

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

二、插值多項(xiàng)式選取原則

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

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

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

三、基本多項(xiàng)式插值方法

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

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

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

(1)

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

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

(2)

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

(3)

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

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

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

四、實(shí)例分析

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

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

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

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

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

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

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

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

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

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

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

MaxRMSxyzxyzPRN012.55102.15321.98610.48390.48880.4776PRN151.61571.84242.36110.43260.44210.5290PRN311.73111.73941.86420.41360.44300.4534

五、結(jié)束語

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

參考文獻(xiàn)：

[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]. 飛行器測控學(xué)報(bào)，2004, 23(4): 43-46.

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

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

[6] 孔巧麗. 用切貝雪夫多項(xiàng)式擬合GPS衛(wèi)星精密坐標(biāo)[J]. 測繪通報(bào)，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] 關(guān)治，陸金甫. 數(shù)值方法[M]. 北京: 清華大學(xué)出版社, 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版.北京: 機(jī)械工業(yè)出版社, 2005：221-224.

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