馬俊，唐詩華，黃鷹，王文杰

(桂林理工大學(xué)土木與建筑工程學(xué)院，廣西桂林 541004)

基于IGS精密星歷的衛(wèi)星位置內(nèi)插方法比較

馬俊?，唐詩華，黃鷹，王文杰

(桂林理工大學(xué)土木與建筑工程學(xué)院，廣西桂林 541004)

在高精度的GPS測量數(shù)據(jù)處理過程中，獲取高精度衛(wèi)星軌道是重要的環(huán)節(jié)，內(nèi)插是獲取任意歷元的精密軌道信息的重要手段，本文利用IGS提供的精密星歷用MATLAB作為工具編程計算，用實(shí)例分析和比較了拉格朗日插值法、牛頓插值法、內(nèi)維爾逐次插值法以及最小二乘擬合在求衛(wèi)星坐標(biāo)當(dāng)中的應(yīng)用效果。前三種方法在16階時內(nèi)插精度最高，內(nèi)維爾插值和牛頓插值有很多相似的特性，方便靈活，龍格現(xiàn)象沒有拉格朗日插值那么劇烈，將會得到廣泛應(yīng)用。

GPS精密星歷；拉格朗日插值；牛頓插值；內(nèi)維爾插值；最小二乘擬合

1 引 言

在GPS定位中，將GPS衛(wèi)星的位置作為已知值，通過計算衛(wèi)星至接收機(jī)的距離來計算測站點(diǎn)的坐標(biāo)。因此，衛(wèi)星的軌道信息是定位的基本要素，而衛(wèi)星的軌道信息是通過星歷給出的，星歷誤差將成為一種起算數(shù)據(jù)誤差。這種誤差直接影響著單點(diǎn)定位的精度，對絕對定位也有一定的影響。為了降低星歷誤差對單點(diǎn)定位和相對定位的影響，提高目標(biāo)在地心坐標(biāo)系的絕對定位精度，IGS(International Geodetic Service)免費(fèi)給出了若干精度等級和不同時間延遲的精密星歷。精密星歷是按一定時間間隔(通常為15 min)給出相應(yīng)歷元時刻GPS衛(wèi)星在地心坐標(biāo)系的三維坐標(biāo)、三維速度及衛(wèi)星鐘改正數(shù)等信息。

GPS精密數(shù)據(jù)處理需要得到間隔1 s的衛(wèi)星坐標(biāo)。GPS接收機(jī)的采樣率一般為30 s或者15 s甚至更密，因此，要想得到某瞬時的衛(wèi)星坐標(biāo)就必須對精密星歷和進(jìn)行高精度、快速的內(nèi)插或擬合。其中較常用的插值方法是拉格朗日多項(xiàng)式插值法，和牛頓多項(xiàng)式插值法。近來也有一些學(xué)者和論述中提到內(nèi)維爾逐次線性插值和切比雪夫多項(xiàng)式擬合法。本文著重分析了拉格朗日多項(xiàng)式插值、牛頓多項(xiàng)式插值、內(nèi)維爾逐次插值以及最小二乘擬合算法，將其在精密星歷衛(wèi)星位置內(nèi)插應(yīng)用方面做了詳細(xì)的比較，得出了一些有益的結(jié)論。

2 插值數(shù)學(xué)模型

2.1 拉格朗日多項(xiàng)式插值

一般情況，對于給定的n＋1個插值節(jié)點(diǎn)，x0

于是所求n次拉格朗日插值多項(xiàng)式為:

拉格朗日多項(xiàng)式函數(shù)模型簡單，插值效率高，收斂速度快，是經(jīng)典的插值方法，但是要增加或者刪除一個插值節(jié)點(diǎn)需要重新構(gòu)造多項(xiàng)式。

2.2 牛頓多項(xiàng)式插值

為f(x)在k＋1點(diǎn)x0，x1，x2，…xk的k階差商。

得到多項(xiàng)式牛頓插值多項(xiàng)式:

實(shí)際計算中我們可以把計算過程中的各階差商保存為一個差商表，如表1所示。這個差商表是個下三角陣，當(dāng)增加節(jié)點(diǎn)時只需要在差商表最下邊增加一行就行了，原先的計算全部有效，差商表計算完成后，只需將差商表中對角線上的元素帶入牛頓多項(xiàng)式，就可計算出插值點(diǎn)結(jié)果。

牛頓多項(xiàng)式插值差商表 表1

2.3 內(nèi)維爾插值

內(nèi)維爾插值是一種線性逐次插值，主要是通過低一次多項(xiàng)式的組合來獲得高一次插值多項(xiàng)式。具體做法是，首先求出若干個一次插值多項(xiàng)式的值，利用它們之間的兩兩組合得到若干個二次插值多項(xiàng)式，進(jìn)一步利用這些二次多項(xiàng)式間的兩兩組合來得到若干個三次插值多項(xiàng)式，依此類推。給定n＋1個已知點(diǎn)(xi，f(xi))令0階內(nèi)插結(jié)果為，Ti，0＝f(xi)，下一階插值結(jié)果是上一階線性內(nèi)插得到的，第j階插值的公式為:

可列出各階內(nèi)維爾的插值表，如表2所示。它是一個下三角矩陣。事實(shí)上后一項(xiàng)就是前一項(xiàng)和左上角一項(xiàng)的線性內(nèi)差，內(nèi)插表中的任何一個元素都是一次內(nèi)插的結(jié)果。內(nèi)維爾插值同樣解決了增加節(jié)點(diǎn)的問題，其數(shù)學(xué)模型簡單，當(dāng)?shù)趇階插值結(jié)果與第i－1階插值結(jié)果小于限差，終止插值。

內(nèi)維爾插值流程表 表2

2.4 最小二乘擬合

對于給定的數(shù)據(jù)(xi，yi)(i＝0，1，2，…m)，要求在給定的函數(shù):

使s?(x)滿足:

取得極小值。

解方程組，可得a0，a1，a2…an，帶入上式中，即可得所求函數(shù)。

3 實(shí)例分析和比較

3.1 精度分析

在IGS官方網(wǎng)站上下載SP3精密星歷文件，本文選取起始?xì)v元在2011年2月23日0時0分0秒，終止歷元時刻在2011年2月23日23時45分0秒的SP3精密星歷文件作為實(shí)驗(yàn)數(shù)據(jù)。本文選用衛(wèi)星序號為14的歷元值作為實(shí)驗(yàn)數(shù)據(jù)，取擬合的時間間隔為30 min，采用對稱內(nèi)插的方法，計算時取30 min間隔的歷元作為內(nèi)插點(diǎn)計算各個內(nèi)插時間段中間時刻的衛(wèi)星位置，以便將內(nèi)插得到的衛(wèi)星位置同精密星歷給出的衛(wèi)星位置進(jìn)行比較。在Matlab 7.1下進(jìn)行實(shí)驗(yàn)，得到14號衛(wèi)星6階～25階的內(nèi)插結(jié)果，IGS向全球GPS用戶提供的GPS精密星歷，其誤差一般小于5 cm，所以這幾種方法內(nèi)插的精度至少要小于5 cm。

表3給出了拉格朗日插值法，牛頓多項(xiàng)式插值法，內(nèi)維爾插值法的6階～25階的標(biāo)準(zhǔn)差，表4給出最小二乘擬合6階～50階的標(biāo)準(zhǔn)差。為了能更好地比較插值效果，筆者根據(jù)表3做出了圖1～圖4，圖1～圖3是前三種方法不同階次插值的標(biāo)準(zhǔn)差曲線圖，圖4是這前三種方法在x方向不同階次插值效果比較圖，以下將根據(jù)這些數(shù)據(jù)圖做詳細(xì)分析。

三種方法不同階次插值標(biāo)準(zhǔn)差 表3

最小二乘擬合各階次的標(biāo)準(zhǔn)差 表4

圖1 拉格朗日不同階次插值標(biāo)準(zhǔn)差

圖2 牛頓插值不同階次插值標(biāo)

圖4 三種方法不同階次標(biāo)準(zhǔn)差

對于拉格朗日插值，從圖1可以看出，當(dāng)階數(shù)過低或過高時，插值精度都較差，當(dāng)階數(shù)取12階～21階之間時，x、y、z方向的標(biāo)準(zhǔn)差均小于50 mm，完全可以滿足各種應(yīng)用。其中插值階數(shù)取14階～18階之間時，三個方向的標(biāo)準(zhǔn)差均小于10 mm，并且三個方向其中取16階時最小，由表1可知三方向的標(biāo)準(zhǔn)差分別是1.17 mm，0.78 mm及0.37 mm。

對于牛頓插值法，從圖2可以看出，和拉格朗日插值法一樣，當(dāng)階數(shù)過低或過高時，插值精度都較差，當(dāng)階數(shù)取12階～21階之間時，x、y、z方向的標(biāo)準(zhǔn)差均小于50 mm，完全可以滿足各種應(yīng)用。其中插值階數(shù)取14階～19階之間時，三個方向的標(biāo)準(zhǔn)差均小于10 mm，取16階時均達(dá)到最小值，由表1可知三方向的標(biāo)準(zhǔn)差分別是1.88 mm，0.78 mm及0.37 mm，幾乎和拉格朗日插值在16階時的精度一樣，從圖2還可以看出從14階～17階之間，三個方向的曲線幾乎重合。

對于內(nèi)維爾逐次插值，從圖3可以看出，從10階～19階x、y、z三個方向的標(biāo)準(zhǔn)差曲線圖與牛頓插值的標(biāo)準(zhǔn)差曲線圖走勢相同，當(dāng)插值取12階～21階時，內(nèi)插精度完全滿足各種應(yīng)用需要。由表4可以看出從10階～19階，牛頓插值與內(nèi)維爾插值在三個方向的標(biāo)準(zhǔn)差相同，所以維爾插值在16階時三個方向標(biāo)準(zhǔn)差均達(dá)到最小值，分別是1.88 mm，0.78 mm及0.37 mm。

由圖1至圖3中還可以看出從10階～20階之間，z方向插值的精度要比x、y方向的精度高。

對于最小二乘擬合，由表4可以看出，從16階～41階最小標(biāo)準(zhǔn)差都小于100 m，但是精度最小也只能達(dá)到分米級，小于16階或者大于41階的，精度都超過100 m，6階和50階都超過萬米，因此最小二乘擬合法不能作為衛(wèi)星位置內(nèi)插的選擇，因此沒有再為其作圖與其他方法進(jìn)行詳細(xì)比較，下文就不再將其詳述

圖4給出了拉格朗日插值，牛頓插值及內(nèi)維爾插值在x方向各階次的標(biāo)準(zhǔn)差曲線，從圖4中我們可以看出從9階～20階，牛頓插值和內(nèi)維爾插值重合，從12階～21階三種方法標(biāo)準(zhǔn)差都小于50 mm，其中16階，17階時三條曲線重合，標(biāo)準(zhǔn)差同時達(dá)到最小值，但是從13階～18階拉格朗日插值標(biāo)準(zhǔn)差波動較大，19階～21階之間，拉格朗日插值法與其他兩種方法相比，標(biāo)準(zhǔn)差較小。

3.2 龍格現(xiàn)象

為了詳細(xì)研究比較龍格現(xiàn)象對拉格朗日插值、牛頓插值和內(nèi)維爾插值的影響，筆者分別作出其在標(biāo)準(zhǔn)差較好的插值階數(shù)下部分歷元的誤差曲線圖，如圖5是拉格朗日插值法部分歷元的誤差曲線圖，圖6，圖7分別是牛頓插值法、內(nèi)維爾插值法的部分歷元的插值曲線圖。

圖5 拉格朗日插值法16階～18階誤差曲線

圖6 牛頓插值法15階～17階誤差曲線

圖7 內(nèi)維爾值法14階～17階誤差曲線

從圖5～圖7中可以看出龍格現(xiàn)象對三種插值方法的影響主要集中在前兩個內(nèi)插節(jié)點(diǎn)和后兩個內(nèi)插節(jié)點(diǎn)，三種方法在16階，17階時插值精度最好，其他階次在前兩個和后兩個插值節(jié)點(diǎn)震蕩都較明顯。對于拉格朗日插值法16階精度最好的內(nèi)插節(jié)點(diǎn)是從第9個～35個歷元，精度都在10 mm以內(nèi)；牛頓插值法16階精度最好的內(nèi)插節(jié)點(diǎn)也是從第9個～35個歷元，精度同拉格朗日；對于內(nèi)維爾插值法，精度小于10 mm的內(nèi)插節(jié)點(diǎn)是從第7個歷元～第39個歷元?？偟膩碚f，內(nèi)維爾逐次插值法在16階次精度小于10 mm的內(nèi)插節(jié)點(diǎn)最多，兩端龍格震蕩最為平緩。

4 結(jié) 論

(1)最小二乘擬合精度最好的時候只能達(dá)到分米級，因此不能作為GPS衛(wèi)星位置內(nèi)插的方法。而將拉格朗日多項(xiàng)式插值、牛頓多項(xiàng)式插值和內(nèi)維爾逐次線性插值用于GPS精密星歷插值處理，無論使用哪一種方法，插值階數(shù)都不應(yīng)少于12階，最高不能超過21階，更高階數(shù)對精度提高的貢獻(xiàn)非常有限，甚至?xí)a(chǎn)生更大的差值，并且插值過程中盡量將插值點(diǎn)位置選在節(jié)點(diǎn)的中間，并盡可能選取偶數(shù)個節(jié)點(diǎn)。在進(jìn)行16階，17階的插值時，三種插值方法精度最高，三個方向的標(biāo)準(zhǔn)差最大不超過3 mm。隨著節(jié)點(diǎn)數(shù)的增加，插值精度反而下降。

(2)拉格朗日插值模型簡單、結(jié)構(gòu)緊湊，利用基函數(shù)就很容易得到插值函數(shù)，是經(jīng)典的插值方法。但是當(dāng)插值節(jié)點(diǎn)增加、減少或其位置變化時，整個插值多項(xiàng)式的結(jié)構(gòu)都會改變，需要重新計算增加了算法的復(fù)雜度，且當(dāng)增大插值階數(shù)時容易出現(xiàn)龍格現(xiàn)象。

(3)內(nèi)維爾插值和牛頓插值具有十分相似的特性，在14階～19階，三個方向的精度均達(dá)到毫米級。模型簡單、計算量小，在內(nèi)插節(jié)點(diǎn)兩端龍格震蕩現(xiàn)象比拉格朗日插值平緩；程序方便靈活，當(dāng)加入新的節(jié)點(diǎn)時，可保留原有計算結(jié)果而不降低程序的效率，對于短時間內(nèi)插衛(wèi)星位置無疑是一個好的選擇，是一種較實(shí)用的內(nèi)插方法。

[1] 陳兆林，朱彥武，張書畢.一種GPS精密星歷的多項(xiàng)式插值方法[J].工程勘察，2007(4)

[2] 卓寧.GPS數(shù)據(jù)處理中衛(wèi)星精密星歷插值方法分析[J].戰(zhàn)術(shù)導(dǎo)彈技術(shù)，2010(4)

[3] 陳遠(yuǎn)鴻.基于精密星歷的計算衛(wèi)星位置的方法[J].測繪信息與工程，2008(2)

[4] 李明峰，江國焰，張凱.IGS精密星歷內(nèi)插與擬合法精度的比較[J].大地測量與地球動力學(xué)，2008，28(2)

[5] 李勇軍，丁士軍.不同星歷誤差對靜態(tài)單點(diǎn)定位精度的影響與分析[J].測繪工程，2010，6(19)

[6] 王磊，羅年學(xué).GPS精密星歷的軌道內(nèi)插方法比較[J].地理空間信息，2009，7(5)

[7] 馬昌鳳，林偉川.現(xiàn)代數(shù)值計算方法(MATLAB版)[M].北京:科學(xué)出版社，2008

[8] 張德豐.MATLAB數(shù)值分析與應(yīng)用[M].北京:國防工業(yè)出版社，2007

[9] Jan Kouba，Pierre Héroux.Precise Point Positioning Using IGS Orbit and Clock Products[J].GPS Solutions，2001，5(2)

IGS Precise Ephemeris Based on the Comparison of the Ways for the Satellite Position Interpolation

Ma Jun，Tang Shihua，Huang Ying，Wang Wenjie
(Guilin University of Technology，Construction Engineering College，Guilin 541004，China)

In the high－precision GPS measurement data processing，to obtain high－precision satellite orbit is an important link.Interpolation is an important means for obtaining any precise orbit information in any epoch.This article provided by IGS precise ephemeris using MATLAB as a tool for programming comparatives and analyses the ways which contain Lagrange interpolation，Newton interpolation，Neville interpolation and least square fitting with the example that these ways is used to obtain the position of satellite.The first three in the 16－order time interpolation have the highest accuracy，Neville interpolation and Newton Interpolation has many similar features，it is convenient and flexible.Its runge phenomenon is less severe than Lagrange interpolation.So it will be widely used.

GPS Precise Ephemeris；Lagrange Interpolation；Newton interpolation；Neville interpolation；leastsquares approximation

2011—04—01

馬俊(1985—)，男，碩士研究生，主要研究方向:GPS單點(diǎn)定位技術(shù)。

1672－8262(2011)05－89－05

P228

A