高采樣率下GPS衛星軌道坐標插值方法比較*

王青平 陳 光 陳超賢 趙文波

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

高采樣率下GPS衛星軌道坐標插值方法比較*

王青平 陳 光 陳超賢 趙文波

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

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

GPS精密星歷;拉格朗日插值;牛頓插值;內維爾插值;多項式系數求解法

1 引言

GPS定位是在GPS衛星的空間坐標已知的情況下，通過測量衛星至天線的距離來計算天線的空間坐標。因此快速而準確地獲取衛星的空間坐標是GPS定位首要解決的關鍵問題。

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

精密星歷文件每15分鐘更新一次衛星的坐標，而GPS接收機的采樣率一般為30秒、1秒甚至更高，需要用插值的方法得到所需歷元的衛星位置。快速而準確地對精密星歷進行插值成為GPS數據處理中的一項重要工作，特別是高頻采樣［4］。

在實際應用中，選用何種插值算法、選用幾次多項式可滿足實際的精度要求，時常成為困擾眾多用戶的問題。為此，本文將給出各插值算法的理論乘/除法運算量，以精密星歷給出的衛星軌道坐標為原始數據，比較三種不同采樣率下各種插值程序算法的實際運行時間，并對其插值結果進行分析。

2 插值算法

2.1 常用的插值算法

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

拉格朗日插值模型簡單，形式對稱，是經典的插值方法，但是當精度不滿足要求而需要增加插值節點時，原來的插值多項式不能使用，需重新構造新的插值多項式;而牛頓插值和內維爾逐次線性插值在增加新的插值節點時，可以重新利用前面的計算結果，常用于給定精度的變階次多項式插值;三角函數插值基于GPS衛星軌道的準周期性，選用三角函數作為基函數，物理涵義清晰，但過程中涉及到大量三角函數運算，運算量較大。

2.2 插值多項式的階次

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

2.3 多項式系數求解法

對于n次多項式:

把t0，t1，…，tn衛星的空間坐標y0，y1，…，yn代入后，可得n+1個關于多項式系數a0，a1，…，an線性方程組，可通過高斯消去法可求解得多項式系數a0，a1，…，an，再利用秦九韶算法［6］可求出插值區間內任意時刻衛星的空間坐標。

多項式系數求解法［6］需解線性方程組，計算量稍大，早期的采樣率較低，故沒有得到廣泛的應用。

3 插值算法的比較

3.1 理論運算量的比較

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

牛頓插值的乘/除法運算量僅為拉格朗日插值的1/4，為內維爾插值的1/3，是一種比較理想的插值方法。但是當GPS接收機采樣率較高時，所需要的插值點m也越多，特別是當m?n時，多項式系數求解法占有絕對的優勢。

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

3.2 理論曲線計算結果的比較

考慮到GPS衛星運行的周期性，選用下列函數作為理論曲線進行比較:

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

1)各算法的絕對誤差量級都是很小的。在精度要求范圍內，各種算法均是可行的;

2)各算法在中央區域的絕對誤差都很小。在實際插值過程中，待插值的時刻最好位于插值節點中央區域附近，使得插值的誤差最小;

3)拉格朗日多項式插值在插值區間兩端附近絕對誤差較大，出現劇烈的龍格震蕩現象。

圖1 插值節點和插值曲線Fig.1 Interpolation nodes and interpolation curve

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

3.3 實際衛星空間坐標計算結果的比較

從 IGS 網站(ftp://garner.ucsd.edu/pub/products/1721/igs17212.sp3.Z)下載 igs17212.sp3 精密星歷文件:歷元始于 2013-01-01T00:00:00，止于2013-01-01T23:45:00。選用前14個歷元序號為1的衛星的橫坐標作為插值節點，對中央區域(1時30分至1時45分，圖3中藍色區域)進行1 Hz的采樣插值。

因無法悉知非插值節點處衛星的真實坐標，只能通過比較四種插值算法間的離散程度。圖4給出四種插值算法的計算結果間的最大誤差分布，很明顯四種插值算法在該區間中差異都很小，并且分布均勻。表明待插值的時刻位于插值節點的中央時，無論使用何種算法，計算結果均可取。

圖3 衛星的精密星歷以及插值節點分布Fig.3 IGS precise ephemeris and interpolation nodes distribution

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

3.4 計算時間的比較

使用各插值算法分別對30 s、1 Hz以及50 Hz采樣率的插值點進行插值，表2給出四種算法的運行時間(代碼通過MinGW GCC編譯，運行環境為I5的CPU、4G內存、Win7操作系統)，不難看出:

1)拉格朗日多項式插值算法的運行效率很低;

2)牛頓插值與內維爾插值的運行時間相當，運行時間約為拉格朗日插值的1/3;

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

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

4 結論

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

2)當待插值的歷元位于插值節點中央時，無論采用哪種算法進行插值，都能取得滿意的結果。當待插值的歷元插值節點區間兩端附近時，牛頓插值和多項式系數求解法最為穩定，而拉格朗日插值極易產生龍格現象;

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

由于早期接收機的采樣率較低，多項式系數求解法并沒有得到廣泛的使用。隨著科技的發展，主流接收機均朝著高采樣率方向發展，相信不久的將來多項式系數求解法將得到廣泛的應用。

1 王曉明，成英燕，劉立.基于階次組合的GPS精密星歷插值研究［J］.大地測量與地球動力學，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軌道實時預報對比研究［J］.大地測量與地球動力學，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精密單點定位快速收斂方法研究［J］.大地測量與地球動力學，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］.大地測量與地球動力學，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 王世儒，等.計算方法［M］.西安:西安電子科技大學出版社，1996.(Wang Shiru，et al.Calculation method［M］.Xi’an:Xidian University Publishing House，1996)

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

8 張守建，等.兩種IGS精密星歷插值方法的比較分析［J］.大地測量與地球動力學，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 邱蕾，廖遠琴，花向紅.基于IGS精密星歷的衛星坐標插值［J］.測繪工程，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精密星歷的衛星位置內插方法比較［J］.城市勘測，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］.測繪工程，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，Fuzhou350003)

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);中國地震局地震科技星火計劃項目(XH13011Y)

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