張 志,廖 瑛,文援蘭,楊雅君,李 志,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.