劉剛，邵楠

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

1 引言

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

用于內(nèi)插的方法有很多種，包括:切比雪夫多項式插值、拉格朗日多項式插值、牛頓多項式插值、三角函數(shù)多項式插值等。目前在GPS數(shù)據(jù)處理中，應(yīng)用比較多的是拉格朗日多項式插值和切比雪夫多項式插值。本文首先對這兩種插值方法進行了介紹，然后通過自編程序?qū)崿F(xiàn)了這兩種插值方法并對兩種插值方法的結(jié)果進行比較分析。

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

2.1 切比雪夫多項式插值

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

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

則衛(wèi)星坐標(biāo)X，Y，Z分量可用如下切比雪夫多項式表示:

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

根據(jù)精密星歷，設(shè)衛(wèi)星的Xk坐標(biāo)為觀測值，則誤差方程為:

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

則有向量表達式:

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

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

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

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

mx，my，mz和 mp分別表示衛(wèi)星坐標(biāo) X，Y，Z 方向上的中誤差和點位中誤差。

2.2 拉格朗日多項式插值

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

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

3 計算實例及比較分析

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

不同階次多項式內(nèi)插中心點處精度比較 表1

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

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

圖1 切比雪夫多項式插值結(jié)果

圖2 拉格朗日多項式插值結(jié)果

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

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

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

4 結(jié)語

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

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

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

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

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

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

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