高采樣率下GPS衛(wèi)星軌道坐標(biāo)插值方法比較*

王青平 陳 光 陳超賢 趙文波

(福建省地震局，福州 350003)

高采樣率下GPS衛(wèi)星軌道坐標(biāo)插值方法比較*

王青平 陳 光 陳超賢 趙文波

(福建省地震局，福州 350003)

比較拉格朗日、牛頓與內(nèi)維爾3種插值算法的運(yùn)算量、精度和運(yùn)行時(shí)間。結(jié)果表明:在精度要求范圍內(nèi)各算法均是可取的，但拉格朗日插值在插值節(jié)點(diǎn)兩端易產(chǎn)生龍格現(xiàn)象;在50 Hz采樣率插值實(shí)驗(yàn)中，多項(xiàng)式系數(shù)求解法的運(yùn)行時(shí)間僅為拉格朗日插值的1/45，為牛頓和內(nèi)維爾插值的1/15。

GPS精密星歷;拉格朗日插值;牛頓插值;內(nèi)維爾插值;多項(xiàng)式系數(shù)求解法

1 引言

GPS定位是在GPS衛(wèi)星的空間坐標(biāo)已知的情況下，通過(guò)測(cè)量衛(wèi)星至天線的距離來(lái)計(jì)算天線的空間坐標(biāo)。因此快速而準(zhǔn)確地獲取衛(wèi)星的空間坐標(biāo)是GPS定位首要解決的關(guān)鍵問(wèn)題。

獲取GPS衛(wèi)星的軌道坐標(biāo)通常有兩種途徑［1］:一種是根據(jù)廣播星歷計(jì)算出星歷參數(shù)，進(jìn)而求出衛(wèi)星的空間坐標(biāo)，另一種是根據(jù)精密星歷計(jì)算衛(wèi)星軌道坐標(biāo)。廣播星歷每?jī)尚r(shí)更新一次，實(shí)時(shí)發(fā)送給用戶，精度為2 m;精密星歷由IGS等國(guó)際組織發(fā)布，包括最終精密星歷(igs)、快速精密星歷(igr)以及預(yù)報(bào)精密星歷(igu)［2］，目前最終精密星歷有12～18天的延遲，GPS最終精密星歷標(biāo)稱精度為2.5 cm、GLONASS標(biāo)稱精度為5 cm，被廣泛應(yīng)用于事后高精度定位。特別是近些年興起的非差精密單點(diǎn)定位技術(shù)(PPP，Precise Point Positioning)，它利用 IGS提供的精密星歷獲得單臺(tái)測(cè)站厘米級(jí)的定位精度［3］。精密星歷的插值精度以及計(jì)算效率將直接影響到PPP等事后高精度數(shù)據(jù)處理的精度和效率。

精密星歷文件每15分鐘更新一次衛(wèi)星的坐標(biāo)，而GPS接收機(jī)的采樣率一般為30秒、1秒甚至更高，需要用插值的方法得到所需歷元的衛(wèi)星位置?？焖俣鴾?zhǔn)確地對(duì)精密星歷進(jìn)行插值成為GPS數(shù)據(jù)處理中的一項(xiàng)重要工作，特別是高頻采樣［4］。

在實(shí)際應(yīng)用中，選用何種插值算法、選用幾次多項(xiàng)式可滿足實(shí)際的精度要求，時(shí)常成為困擾眾多用戶的問(wèn)題。為此，本文將給出各插值算法的理論乘/除法運(yùn)算量，以精密星歷給出的衛(wèi)星軌道坐標(biāo)為原始數(shù)據(jù)，比較三種不同采樣率下各種插值程序算法的實(shí)際運(yùn)行時(shí)間，并對(duì)其插值結(jié)果進(jìn)行分析。

2 插值算法

2.1 常用的插值算法

常用的插值算法有拉格朗日多項(xiàng)式插值［1，5－11］、牛頓多項(xiàng)式插值［6，7，10］、內(nèi)維爾逐次線性插值［1，5，9，10］以及三角函數(shù)插值［8］等。

拉格朗日插值模型簡(jiǎn)單，形式對(duì)稱，是經(jīng)典的插值方法，但是當(dāng)精度不滿足要求而需要增加插值節(jié)點(diǎn)時(shí)，原來(lái)的插值多項(xiàng)式不能使用，需重新構(gòu)造新的插值多項(xiàng)式;而牛頓插值和內(nèi)維爾逐次線性插值在增加新的插值節(jié)點(diǎn)時(shí)，可以重新利用前面的計(jì)算結(jié)果，常用于給定精度的變階次多項(xiàng)式插值;三角函數(shù)插值基于GPS衛(wèi)星軌道的準(zhǔn)周期性，選用三角函數(shù)作為基函數(shù)，物理涵義清晰，但過(guò)程中涉及到大量三角函數(shù)運(yùn)算，運(yùn)算量較大。

2.2 插值多項(xiàng)式的階次

選用幾次多項(xiàng)式對(duì)GPS衛(wèi)星的空間坐標(biāo)進(jìn)行插值，使得插值的精度最高，許多學(xué)者對(duì)此做了研究:張守建［8］認(rèn)為10階拉格朗日多項(xiàng)式內(nèi)插精度最高，11階的三角函數(shù)多項(xiàng)式插值內(nèi)插精度最好;邱蕾［9］認(rèn)為當(dāng)插值點(diǎn)位于節(jié)點(diǎn)中央時(shí)，9階以上拉格朗日和內(nèi)維爾插值可獲得與精密星歷相當(dāng)?shù)木?，同時(shí)多項(xiàng)式的次數(shù)最好為奇數(shù)，待插值的時(shí)刻位于節(jié)點(diǎn)的中間，插值的誤差也最小;馬?。?0］認(rèn)為拉格朗日、牛頓以及內(nèi)維爾插值在16階精度最高;何玉晶［11］認(rèn)為9階拉格朗日插值的精度最高;王曉明［1］則提出一種高階插值與低階插值相結(jié)合的插值算法來(lái)提高插值精度。雖說(shuō)眾學(xué)者的看法不盡相同，但均認(rèn)為在進(jìn)行多項(xiàng)式插值時(shí)，適當(dāng)提高多項(xiàng)式的次數(shù)，可以提高插值的精度，但在插值節(jié)點(diǎn)兩端附近高次多項(xiàng)式插值，容易產(chǎn)生龍格(Runge)現(xiàn)象［6，7］，特別是拉格朗日多項(xiàng)式插值［10］，因而會(huì)影響整個(gè)插值的精度，故盲目地提高插值多項(xiàng)式的次數(shù)是不可取的。為簡(jiǎn)單起見(jiàn)，本文只考察13次多項(xiàng)式的各種插值算法。

2.3 多項(xiàng)式系數(shù)求解法

對(duì)于n次多項(xiàng)式:

把t0，t1，…，tn衛(wèi)星的空間坐標(biāo)y0，y1，…，yn代入后，可得n+1個(gè)關(guān)于多項(xiàng)式系數(shù)a0，a1，…，an線性方程組，可通過(guò)高斯消去法可求解得多項(xiàng)式系數(shù)a0，a1，…，an，再利用秦九韶算法［6］可求出插值區(qū)間內(nèi)任意時(shí)刻衛(wèi)星的空間坐標(biāo)。

多項(xiàng)式系數(shù)求解法［6］需解線性方程組，計(jì)算量稍大，早期的采樣率較低，故沒(méi)有得到廣泛的應(yīng)用。

3 插值算法的比較

3.1 理論運(yùn)算量的比較

表1為不同幾種插值算法對(duì)m個(gè)插值點(diǎn)進(jìn)行插值的理論乘/除法運(yùn)算量［5－7］。拉格朗日插值、牛頓插值以及內(nèi)維爾插值的時(shí)間復(fù)雜度均為O(mn2)。多項(xiàng)式系數(shù)求解法先用O(n3)的運(yùn)算量計(jì)算出多項(xiàng)式的系數(shù)，后續(xù)使用秦九韶算法對(duì)m個(gè)插值點(diǎn)求值僅需O(mn)的運(yùn)算量，總運(yùn)算量為O(n3+mn)。

牛頓插值的乘/除法運(yùn)算量?jī)H為拉格朗日插值的1/4，為內(nèi)維爾插值的1/3，是一種比較理想的插值方法。但是當(dāng)GPS接收機(jī)采樣率較高時(shí)，所需要的插值點(diǎn)m也越多，特別是當(dāng)m?n時(shí)，多項(xiàng)式系數(shù)求解法占有絕對(duì)的優(yōu)勢(shì)。

表1 各插值方法理論乘/除法運(yùn)算量的比較Tab.1 Comparison of computation of different interpolation methods

3.2 理論曲線計(jì)算結(jié)果的比較

考慮到GPS衛(wèi)星運(yùn)行的周期性，選用下列函數(shù)作為理論曲線進(jìn)行比較:

式中b0，b1，…，bp為系數(shù)，ω 為圓頻率。張守建［8］認(rèn)為11階的三角函數(shù)多項(xiàng)式與衛(wèi)星空間坐標(biāo)擬合較好。隨機(jī)產(chǎn)生一組b0，b1，…，b11進(jìn)而產(chǎn)生14個(gè)插值節(jié)點(diǎn)(圖1中的紅色標(biāo)記)，進(jìn)行各種多項(xiàng)式插值(n=13)，并將插值結(jié)果與式(2)的理論值進(jìn)行比較，誤差分布如圖2所示。不難看出:

1)各算法的絕對(duì)誤差量級(jí)都是很小的。在精度要求范圍內(nèi)，各種算法均是可行的;

2)各算法在中央?yún)^(qū)域的絕對(duì)誤差都很小。在實(shí)際插值過(guò)程中，待插值的時(shí)刻最好位于插值節(jié)點(diǎn)中央?yún)^(qū)域附近，使得插值的誤差最小;

3)拉格朗日多項(xiàng)式插值在插值區(qū)間兩端附近絕對(duì)誤差較大，出現(xiàn)劇烈的龍格震蕩現(xiàn)象。

圖1 插值節(jié)點(diǎn)和插值曲線Fig.1 Interpolation nodes and interpolation curve

圖2 不同插值算法與理論值的絕對(duì)誤差Fig.2 Absolute error of different interpolation algorithms with the theoretical value

3.3 實(shí)際衛(wèi)星空間坐標(biāo)計(jì)算結(jié)果的比較

從 IGS 網(wǎng)站(ftp://garner.ucsd.edu/pub/products/1721/igs17212.sp3.Z)下載 igs17212.sp3 精密星歷文件:歷元始于 2013-01-01T00:00:00，止于2013-01-01T23:45:00。選用前14個(gè)歷元序號(hào)為1的衛(wèi)星的橫坐標(biāo)作為插值節(jié)點(diǎn)，對(duì)中央?yún)^(qū)域(1時(shí)30分至1時(shí)45分，圖3中藍(lán)色區(qū)域)進(jìn)行1 Hz的采樣插值。

因無(wú)法悉知非插值節(jié)點(diǎn)處衛(wèi)星的真實(shí)坐標(biāo)，只能通過(guò)比較四種插值算法間的離散程度。圖4給出四種插值算法的計(jì)算結(jié)果間的最大誤差分布，很明顯四種插值算法在該區(qū)間中差異都很小，并且分布均勻。表明待插值的時(shí)刻位于插值節(jié)點(diǎn)的中央時(shí)，無(wú)論使用何種算法，計(jì)算結(jié)果均可取。

圖3 衛(wèi)星的精密星歷以及插值節(jié)點(diǎn)分布Fig.3 IGS precise ephemeris and interpolation nodes distribution

圖4 各算法間的最大偏差分布Fig.4 Maximum error distribution among the interpolation algorithms

3.4 計(jì)算時(shí)間的比較

使用各插值算法分別對(duì)30 s、1 Hz以及50 Hz采樣率的插值點(diǎn)進(jìn)行插值，表2給出四種算法的運(yùn)行時(shí)間(代碼通過(guò)MinGW GCC編譯，運(yùn)行環(huán)境為I5的CPU、4G內(nèi)存、Win7操作系統(tǒng))，不難看出:

1)拉格朗日多項(xiàng)式插值算法的運(yùn)行效率很低;

2)牛頓插值與內(nèi)維爾插值的運(yùn)行時(shí)間相當(dāng)，運(yùn)行時(shí)間約為拉格朗日插值的1/3;

表2 各插值方法的運(yùn)行時(shí)間(ms)Tab.2 Running time of different interpolation algorithms(unit:ms)

3)多項(xiàng)式系數(shù)求解法的效率相當(dāng)高，與其他三種算法相比均有一個(gè)量級(jí)的提升，且采樣率越高性能提升得明顯。特別是50 Hz采樣率運(yùn)行時(shí)間僅為拉格朗日插值的1/45，為牛頓插值以及內(nèi)維爾插值的1/15;

4 結(jié)論

1)拉格朗日插值、牛頓插值和內(nèi)維爾插值算法的時(shí)間復(fù)雜度均為O(mn2)，而多項(xiàng)式系數(shù)求解法的時(shí)間復(fù)雜度為O(n3+mn)。

2)當(dāng)待插值的歷元位于插值節(jié)點(diǎn)中央時(shí)，無(wú)論采用哪種算法進(jìn)行插值，都能取得滿意的結(jié)果。當(dāng)待插值的歷元插值節(jié)點(diǎn)區(qū)間兩端附近時(shí)，牛頓插值和多項(xiàng)式系數(shù)求解法最為穩(wěn)定，而拉格朗日插值極易產(chǎn)生龍格現(xiàn)象;

3)多項(xiàng)式系數(shù)求解法的運(yùn)行效率最高。50 Hz采樣率下，其效率約為牛頓插值和內(nèi)維爾插值的15倍;約為拉格朗日插值的45倍。

由于早期接收機(jī)的采樣率較低，多項(xiàng)式系數(shù)求解法并沒(méi)有得到廣泛的使用。隨著科技的發(fā)展，主流接收機(jī)均朝著高采樣率方向發(fā)展，相信不久的將來(lái)多項(xiàng)式系數(shù)求解法將得到廣泛的應(yīng)用。

1 王曉明，成英燕，劉立.基于階次組合的GPS精密星歷插值研究［J］.大地測(cè)量與地球動(dòng)力學(xué)，2011，(4):103－106.(Wang Xiaoming，Cheng Yingyan and Liu Li.Research on GPS precise ephemeris interpolation based on high/low order integrated［J］.Journal of Geodesy and Geodynamics，2011，(4):103 －106)

2 劉偉平，等.GPS軌道實(shí)時(shí)預(yù)報(bào)對(duì)比研究［J］.大地測(cè)量與地球動(dòng)力學(xué)，2012，(1):94 －96.(Liu Weiping，et al.Comparison and study on real-time prediction of GPS orbit［J］.Journal of Geodesy and Geodynamics，2012，(1):94 －96)

3 郝明，王慶良，崔篤信.GPS精密單點(diǎn)定位快速收斂方法研究［J］.大地測(cè)量與地球動(dòng)力學(xué)，2009，(2):88－91.(Hao Ming，Wang Qingliang and Cui Duxin.Study on fast convergence method in precise point positioning［J］.Journal of Geodesy and Geodynamics，2009，(2):88 －91)

4 吳繼忠，吳文壇.利用高頻GPS進(jìn)行地震動(dòng)態(tài)變形分析及地震定位［J］.大地測(cè)量與地球動(dòng)力學(xué)，2012，(2):20－23.(Wu Jizhong and Wu Wentan.Kinematic analysis of coseismic deformation and seismic location using high-rate GPS［J］.Journal of Geodesy and Geodynamics，2012，(2):20－23)

5 William H Press，et al.Numerical recipes:The art of scientific computing(3rd edition)［M］:New York:Cambridge University Press.2007.

6 王世儒，等.計(jì)算方法［M］.西安:西安電子科技大學(xué)出版社，1996.(Wang Shiru，et al.Calculation method［M］.Xi’an:Xidian University Publishing House，1996)

7 徐萃薇，孫繩武.計(jì)算方法引論［M］.北京:高等教育出版社，2002.(Xu Cuiwei and Sun Shengwu.Introduction to computational methods［M］.Beijing:Higher Education Press，2002)

8 張守建，等.兩種IGS精密星歷插值方法的比較分析［J］.大地測(cè)量與地球動(dòng)力學(xué)，2007，(2):80－83.(Zhang Shoujian，et al.Comparative analysis on two methods for igs precise ephemeris interpolation［J］.Journal of Geodesy and Geodynamics，2007，(2):80 －83)

9 邱蕾，廖遠(yuǎn)琴，花向紅.基于IGS精密星歷的衛(wèi)星坐標(biāo)插值［J］.測(cè)繪工程，2008，17(4):15 － 18.(Qiu Lei，Liao Yuanqin and Hua Xianghong.Interpolation methods for precision ephemeris based on IGS［J］.Engineering of Surveying and Mapping，2008，17(4):15 －18)

10 馬俊，等.基于IGS精密星歷的衛(wèi)星位置內(nèi)插方法比較［J］.城市勘測(cè)，2011，5:89 － 93.(Ma Jun，et al.IGS precise ephemeris based on the comparison of the ways for the satellite position interpolation［J］.Urban Geotechnical Investigation ＆ Surveying，2011，5:89 －93)

11 何玉晶，楊力.基于拉格朗日插值方法的GPS IGS精密星歷插值分析［J］.測(cè)繪工程，2011，20(5):60－ 62.(He Yujing and Yang Li.Analysis of interpolation results on GPS IGS precise ephemeris based on Lagrange interpolation［J］.Engineering of Surveying and Mapping，2011，20(5):60－62)

COMPARISON ON INTERPOLATION METHOD FOR HIGH SAMPLING RATE OF GPS SATELLITE ORBIT

Wang Qingping，Chen Guang，Chen Chaoxian and Zhao Wenbo
(Earthquake Administration of Fujian Province，F(xiàn)uzhou350003)

The computation and the accuracy and running time of Lagrange interpolation，Newton interpolation，Neville successive linear interpolation are compared.The result indicates that each algorithm is feasible within the required accuracy range，but it is easy to produce Runge phenomenon in Lagrange interpolation.Coefficients of the interpolating polynomial algorithm is about 45 times faster than Lagrange interpolation and 15 times faster than Newton and Neville interpolations at 50 Hz sampling rate interpolation experiments.

GPS precise ephemeris;Lagrange polynomial interpolation;Newton polynomial interpolation;Neville successive linear interpolation;polynomial coefficients solving method

P207

A

1671-5942(2013)05-0049-04

2013-02-09

福建省地震局青年科技基金(Y201201、Y201003);中國(guó)地震局地震科技星火計(jì)劃項(xiàng)目(XH13011Y)

王青平，男，1984年生，博士，工程師，主要從事GPS形變監(jiān)測(cè)研究.E－mail:wacs5@126.com