張 志,廖 瑛,文援蘭,楊雅君,李 志,2
(1. 國(guó)防科技大學(xué) 航天與材料工程學(xué)院,湖南 長(zhǎng)沙 410073; 2. 中國(guó)空間技術(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)星位置,就必須對(duì)精密星歷進(jìn)行高精度的插值。在GPS數(shù)據(jù)后處理中,用戶(hù)需要選擇合理的數(shù)學(xué)模型對(duì)GPS精密星歷進(jìn)行插值或擬合,才能獲得任意時(shí)刻的衛(wèi)星坐標(biāo),因而研究如何利用精密星歷快速地插值出觀測(cè)瞬間衛(wèi)星的位置成為實(shí)際應(yīng)用中一項(xiàng)十分必要的工作。
對(duì)GPS精密星歷的內(nèi)插法主要有多項(xiàng)式內(nèi)插法和三角內(nèi)插法,采用多項(xiàng)式內(nèi)插法和三角內(nèi)插法得到的精度基本相當(dāng)[2]。對(duì)于多項(xiàng)式內(nèi)插法,要想獲得理想的精度就要合理地選擇多項(xiàng)式的階數(shù)和插值的時(shí)間間隔。目前,國(guó)內(nèi)用于精密星歷插值較常用的方法有Lagrange插值法、Neville插值法[3-4]及Chebyshev多項(xiàng)式插值法[5-6]等,但對(duì)基本多項(xiàng)式插值法的研究很少。Mark Schenewerk[2]和Milan Horemu?[7]均提到了基本多項(xiàng)式插值法用于精密星歷的插值計(jì)算,前者只給出了概括性的分析說(shuō)明,后者雖對(duì)多項(xiàng)式插值法在精密星歷的插值過(guò)程中出現(xiàn)的問(wèn)題進(jìn)行了分析,但相關(guān)文獻(xiàn)并未給出具體的插值計(jì)算方法。
本文利用Householder變換求解基本多項(xiàng)式擬合系數(shù)來(lái)避免對(duì)病態(tài)的法方程系數(shù)矩陣直接求逆計(jì)算二乘解,并對(duì)擬合曲線進(jìn)行插值計(jì)算,為避免高階基本多項(xiàng)式插值過(guò)程中出現(xiàn)的龍格現(xiàn)象的影響,給出了多項(xiàng)式擬合區(qū)間長(zhǎng)度、有效區(qū)間長(zhǎng)度和基本插值多項(xiàng)式階數(shù)的選取準(zhǔn)則,最后對(duì)插值的結(jié)果進(jìn)行比較分析。
對(duì)于插值多項(xiàng)式(插值函數(shù))來(lái)說(shuō),為了使插值函數(shù)的數(shù)值計(jì)算具有可操作性,通常取插值函數(shù)類(lèi)為有限維子空間。為使插值便于實(shí)現(xiàn),插值函數(shù)類(lèi)子空間的一組基函數(shù)應(yīng)當(dāng)容易計(jì)算。一旦確定了插值函數(shù)類(lèi)子空間的一組基函數(shù),插值問(wèn)題就可以轉(zhuǎn)化為求一組組合系數(shù)的問(wèn)題。
由定理說(shuō)明組合系數(shù)的唯一性:設(shè)φi,i=1,2,…,n+1是線性無(wú)關(guān)的多項(xiàng)式,那么任一多項(xiàng)式p可唯一地由φ1t,φ2t,…,φn+1t的線性組合來(lái)表示[8]。
設(shè)線性無(wú)關(guān)的一組基函數(shù)1,t,t2,…,tn,則插值多項(xiàng)式(插值函數(shù))pt可唯一地由基函數(shù)線性組合表示pt=c1tn+c2tn-1+…+cnt+cn+1,此插值多項(xiàng)式稱(chēng)為基本多項(xiàng)式。同樣,Chebyshev插值多項(xiàng)式和Legendre插值多項(xiàng)式等均可以唯一地由一組線性無(wú)關(guān)的基函數(shù)表示。以下將以基本多項(xiàng)式為例進(jìn)行分析討論。
給定m個(gè)歷元時(shí)刻t1,t2,…,tm的觀測(cè)點(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)
為最小。此處給出最小二乘問(wèn)題的法方程組
(3)
式中,R1是為n+1階上三角矩陣;R中后面m-(n+1)行元素為0。
利用正交矩陣不改變向量長(zhǎng)度的性質(zhì),又有
這樣用QR分解的最小二乘法求解出基本多項(xiàng)式系數(shù)cxi,同理,可求得cyi、czi,利用這些系數(shù),根據(jù)式(1)便可計(jì)算出t∈t1,tm時(shí)間區(qū)間內(nèi)任意時(shí)刻的衛(wèi)星坐標(biāo)。
為更好地進(jìn)行試驗(yàn)比對(duì),采用Mark Sche-newerk[2]中的數(shù)據(jù),包括時(shí)間間隔(步長(zhǎng))為15 min的精密星歷ECF_15MI.200,文中稱(chēng)為試驗(yàn)數(shù)據(jù);作為插值比較的5 min間隔數(shù)據(jù)文件ECF_5MIN.200,文中稱(chēng)為真值數(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)式擬合,然后對(duì)擬合的曲線進(jìn)行5 min間隔的插值,得到的插值數(shù)據(jù)與對(duì)應(yīng)歷元真值數(shù)據(jù)之差稱(chēng)為插值誤差,以此來(lái)說(shuō)明方法的有效性。插值誤差如圖1—圖4所示。本文給出了PRN01星的比較結(jié)果,其他各星結(jié)果類(lèi)似。
圖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ū)域插值精度為毫米級(jí),因此可以稱(chēng)擬合區(qū)間中間2 h區(qū)域?yàn)橛行^(qū)間,這樣,只要每次插值計(jì)算都進(jìn)入有效區(qū)間,則能有效避免龍格現(xiàn)象的影響,從而使擬合達(dá)到毫米級(jí)的精度。對(duì)2 h的有效區(qū)間,最大插值誤差不超過(guò)5 mm,這對(duì)于厘米級(jí)的精密星歷來(lái)說(shuō),是完全可以接受的。
綜合考慮計(jì)算量和擬合階數(shù),本文選取階數(shù)為16,擬合區(qū)間為4 h的多項(xiàng)式,在2 h有效區(qū)間進(jìn)行插值計(jì)算,如圖5所示,對(duì)試驗(yàn)數(shù)據(jù)進(jìn)行4 h擬合,在一天內(nèi)有效區(qū)間插值數(shù)據(jù)與對(duì)應(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á)到毫米級(jí),誤差小于3 mm。由于精密單點(diǎn)定位對(duì)衛(wèi)星星歷誤差的要求為毫米級(jí),因此這并不影響多項(xiàng)式擬合法插值出的數(shù)據(jù)在實(shí)際中的應(yīng)用。用兩個(gè)統(tǒng)計(jì)量Max和RMS來(lái)反映試驗(yàn)數(shù)據(jù)所有衛(wèi)星基本多項(xiàng)式插值結(jié)果,表1列出了PRN01、PRN15和PRN31號(hào)衛(wèi)星的分析結(jié)果。其中,Max指插值結(jié)果三維位置誤差絕對(duì)值的最大值;RMS指內(nèi)插結(jié)果三維位置誤差的均方根誤差。
由表1可知,采用2 h的有效區(qū)間基本多項(xiàng)式插值方法,試驗(yàn)數(shù)據(jù)3顆衛(wèi)星x、y、z3個(gè)方向的插值誤差最大值均不超過(guò)3 mm,3個(gè)方向的均方根小于0.6 mm。對(duì)于厘米級(jí)的精密星歷來(lái)說(shuō),試驗(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
本文對(duì)插值多項(xiàng)式的選取原則進(jìn)行了探討,根據(jù)選取原則設(shè)計(jì)了一種簡(jiǎn)單的基本插值多項(xiàng)式,然后利用基本多項(xiàng)式擬合得到插值多項(xiàng)式系數(shù),并對(duì)精密星歷進(jìn)行插值測(cè)試。試驗(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)則,對(duì)試驗(yàn)數(shù)據(jù)進(jìn)行一天的插值,給出有效區(qū)間之內(nèi)的插值結(jié)果,插值誤差小于3 mm,能夠滿(mǎn)足用戶(hù)定位精度要求。
參考文獻(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]. 飛行器測(cè)控學(xué)報(bào),2004, 23(4): 43-46.
[4] 王曉明,成英燕,劉立. 基于階次組合的GPS精密星歷插值研究[J]. 大地測(cè)量與地球動(dòng)力學(xué),2011, 31(4): 103-106.
[5] 彭澤泉. GPS精密星歷擬合方法的研究[J]. 測(cè)繪科學(xué),2010, 35(S1): 63-65.
[6] 孔巧麗. 用切貝雪夫多項(xiàng)式擬合GPS衛(wèi)星精密坐標(biāo)[J]. 測(cè)繪通報(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] 金凱德,切尼,王國(guó)榮. 數(shù)值分析[M].3版.北京: 機(jī)械工業(yè)出版社, 2005:221-224.
[11] 吉長(zhǎng)東,徐愛(ài)功,馮磊. GPS精密星歷擬合與插值中龍格現(xiàn)象的處理方法[J]. 測(cè)繪科學(xué),2011, 36(6): 169-171.