兩種IGS精密星歷多項式插值方法的比較分析

劉剛，邵楠

(1.重慶市勘測院，重慶 400020; 2.沈陽市勘察測繪研究院，遼寧 沈陽 110004)

1 引言

利用IGS提供的精密星歷和鐘差參數，可以實現高精度的定位。IGS發布的SP3格式的精密星歷通常給出的是15 min或者5 min等間隔的衛星位置和鐘差信息，而在實際應用中可能需要任意歷元時刻的位置和鐘差信息，這就需要利用精密星歷提供的等間隔時刻的信息進行內插方法得到特定的歷元時刻的衛星位置和鐘差，然后再進行后續的計算。因此，高精度、快速實現對精密星歷和鐘差的內插或擬合就成為了GPS數據處理中的一項基礎而重要的內容。

用于內插的方法有很多種，包括:切比雪夫多項式插值、拉格朗日多項式插值、牛頓多項式插值、三角函數多項式插值等。目前在GPS數據處理中，應用比較多的是拉格朗日多項式插值和切比雪夫多項式插值。本文首先對這兩種插值方法進行了介紹，然后通過自編程序實現了這兩種插值方法并對兩種插值方法的結果進行比較分析。

2 切比雪夫多項式插值與拉格朗日多項式插值

2.1 切比雪夫多項式插值

切比雪夫多項式擬合就是根據給定的數據擬合出一個函數，使其在給定點的函數值與給定值之間的方差和最小，且該函數是以切比雪夫多項式為基函數構成的函數［1，2］。

假定欲將在時間間隔［t0，t0+△t］的衛星星歷用n階切比雪夫多項式進行逼近(或標準化)，其中t0和△t分別是開始歷元和擬合時間區間的長度。為了用切比雪夫多項式對衛星軌道標準化，首先時間t∈［t0，t0+△t］變換成 τ∈［－1，1］，變換公式為［3］:

則衛星坐標X，Y，Z分量可用如下切比雪夫多項式表示:

在式(2)中，n為切比雪夫多項式的系數，Cxi，Cyi，Czi分別為X坐標分量、Y坐標分量、Z坐標分量切比雪夫多項式的系數。切比雪夫多項式Ti用以下遞推公式確定:

根據精密星歷，設衛星的Xk坐標為觀測值，則誤差方程為:

誤差方程的矩陣展開式為:

則有向量表達式:

應用最小二乘法平差中的間接平差，可以得出法方程為:

式(6)中，N=BTPB，fex=BTPfx，此處，權矩陣 P 為單位陣，因此，N=BTB，fex=BTfx，則方程(5)的解為:

此時，X各分量便為切比雪夫多項式擬合系數Cxi。同理，對Y、Z分量完成上述類似的計算便可得到關于某顆GPS衛星在［t0，t0+△t］區間內各坐標分量的切比雪夫多項式擬合系數 Cxi，Cyi，Czi(i=0，1，2，…，n)。通過這三組系數，便可利用(1)～(3)式分別計算得到［t0，t0+△t］時間區間內任意時刻的衛星坐標。

為了檢查擬合的質量，還需要對擬合結果進行質量評定。一般來說是利用中誤差來評定的，評定公式如下:

mx，my，mz和 mp分別表示衛星坐標 X，Y，Z 方向上的中誤差和點位中誤差。

2.2 拉格朗日多項式插值

設 y=f(x)是區間［a，b］上的一個實函數，xi(i=0，1，…，n)是［a，b］上 n+1 個互異實數，且 y=f(x)在 xi處的值為yi=f(xi)，則區間［a，b］上任意一點x的n階拉格朗日插值多項式的代數表達式為［4，5］

點xi(i=0，1，…，n)稱為插值節點，包含插值節點的區間［a，b］稱為插值區間。

3 計算實例及比較分析

為了比較兩種方法的插值精度，作者根據上述插值原理利用VC編制了程序來實現兩種插值方法。文中選取了IGS發布的2009年3月29日的精密星歷數據作為算例來研究不同階數多項式下兩種插值方法的精度。由于精密星歷給出的衛星位置的時間間隔為15 min，因此本文計算中以30 min為時間間隔選取內插點內插計算得到各個內插時間段中間時刻的衛星位置，以便將兩種方法內插得到的結果與精密星歷給出的衛星位置進行比較，獲得兩種插值方法的精度。從上述的兩種插值原理中可以看出，當采用n階多項式插值時，至少需要n+1個節點的坐標數據，為了檢驗坐標內插的質量，選取n+1個節點中間已知坐標的n個點作為插值點，計算出這n個插值點的坐標，與已知值進行比較，檢核其精度。作者利用自編程序用兩種插值方法分別計算了32顆衛星在不同階次下各插值點處的平均點位誤差，表1給出了不同階次多項式下兩種插值方法在內插時段中心點出的平均點位誤差。

不同階次多項式內插中心點處精度比較 表1

從表1中可以看出，隨著插值多項式階次的提高，切比雪夫多項式插值的精度先是逐漸提高，當階次達到20次以后，插值精度迅速降低，這稱為高次插值的“龍格”現象。而采用拉格朗日多項式插值，隨著多項式階次的提高，插值精度先是提高，然后趨于穩定。

按照切比雪夫多項式插值時間歸化的原理，兩種多項式插值結果反映為時間與精度的關系如圖1、圖2所示，式中T為所取插值點位于整個插值區段的位置，取中心節點處為0，整個插值區段為(－1，1)。

圖1 切比雪夫多項式插值結果

圖2 拉格朗日多項式插值結果

圖1和圖2分別表示了用不同階次切比雪夫多項式和拉格朗日多項式插值在整個插值區間的精度。從圖中可以看出，兩種插值方法均表現出兩頭大中間小的誤差分布情況，特別是在10階、20階和22階時，兩端處的誤差跳躍較大。從圖中也可以看出，兩種插值方式從10階上升到20階時，走勢基本一致，當多項式階次上升到22階時，兩種方法表現出不同的走勢:切比雪夫多項式插值精度明顯變差，隨時間變化的趨勢也發生了改變，出現了不規則的波動，這是因為階數太高，使數據噪聲納入了擬合模型，影響了模型的精度［6］。而拉格朗日插值法則保持了之前的走勢，離插值中心點較近的地方仍舊保持了較高的精度，在兩端處精度迅速下降。為了更好地分析兩種插值在高階時的差異，抽取圖形局部放大，如圖3所示，可看出在20階時切比雪夫多項式插值精度顯著降低，波動較大，差于拉格朗日插值法。

圖3 20階多項式插值精度對比

由此可以看出，當多項式階數上升到一定程度后，利用切比雪夫多項式插值的結果不再可靠，而拉格朗日插值法仍舊保持高精度的插值。實驗結果顯示，利用兩種方法進行插值時，采用一定階次的多項式進行插值，在插值區間范圍內，80%左右的插值區間均能取得很好的插值效果。根據這一特性，采用切比雪夫多項式進行插值時，可適當提高階次以提高插值區間覆蓋的范圍，在有效范圍內，只需一次計算多項式系數進行存儲即可直接使用，可以免去多項式系數的重復計算。對于一天的數據，需要銜接前后一天的數據，分為幾個時段，分別存儲幾組系數，在應用中根據所需時間直接調用系數進行計算即可。而對于拉格朗日插值法，在各個不同的時間點上進行插值時需要單獨進行運算。

4 結語

本文通過程序實現了利用切比雪夫和拉格朗日多項式對IGS精密星歷進行插值的功能，分析了不同階次多項式對于插值精度的影響，并對兩種方法進行了比較。實驗數據表明，按照 30 min間隔選取插值節點，在多項式階次達到10階以上時，兩種方法均能取得毫米級的插值精度，在一定階次范圍內，兩種插值方法插值精度相當。當多項式階次上升到20階以上時，切比雪夫插值法精度迅速下降，而拉格朗日插值法仍維持較高的精度，說明高次插值時，切比雪夫多項式插值方法“龍格”現象更為突出。根據實驗結果，兩種插值方法的運算都較簡單，精度也較高。在選取插值方法以及多項式階次時，應根據實際需要擬合時間長度來選擇，對于切比雪夫多項式插值法可適當提高多項式階次以實現較長時段的高精度擬合，而對于拉格朗日插值法，在達到要求的插值精度后，再提高多項式的階次已無益于插值精度的提高，因此無需采用過高的階次以提高運算效率。

［1］邱蕾，廖遠琴，華向紅.基于IGS精密星歷的衛星坐標插值［J］.測繪工程，2008，17(4):15～18

［2］劉剛.基于精密星歷的切比雪夫多項式衛星軌道坐標擬合研究［J］.城市勘測，2010(1):53～55

［3］陳正陽，易重海.用切比雪夫多項式進行GPS衛星軌道標準化［J］.礦山測量，2002(2):5～7

［4］曹德欣，曹瓔珞.計算方法(第二版)［M］.徐州:中國礦業大學出版社，2001

［5］朱方生，李大美，李素貞.計算方法［M］.武漢:武漢大學出版社，2003

［6］楊學峰，程鵬飛，方愛平等.利用切比雪夫多項式擬合衛星軌道坐標的研究［J］.測繪通報，2008(12):1～3