基于Neville算法的GPS精密星歷數(shù)據(jù)加密及精度驗(yàn)證①

曾 超，劉智敏

（山東科技大學(xué)測(cè)繪與工程學(xué)院，山東 青島266590）

0 引 言

精密星歷是按一定的時(shí)間間隔（通常是15 min），由若干衛(wèi)星跟蹤站的觀測(cè)數(shù)據(jù)經(jīng)事后處理得到的衛(wèi)星三維坐標(biāo)、三維運(yùn)動(dòng)速度及衛(wèi)星鐘改正數(shù)等信息。國際地球動(dòng)力學(xué)服務(wù)組織（IGS）組織提供的精密星歷通常需要在測(cè)量約2周后獲得，其提供的數(shù)據(jù)采樣間隔為15min，而實(shí)際觀測(cè)中的采樣率一般為30s、15s甚至更密。為了滿足精密單點(diǎn)定位等一系列高精度衛(wèi)星定位服務(wù)的需求，常采用適當(dāng)?shù)姆椒▽?duì)精密星歷提供的數(shù)據(jù)進(jìn)行加密［1，2]。

精密星歷數(shù)據(jù)的加密主要包括衛(wèi)星在空間的三維坐標(biāo)的加密及鐘差改正數(shù)加密。常用的加密方法有拉格朗日多項(xiàng)式插值（Lagrange Polynomial Interpolation）、切比雪夫多項(xiàng)式擬合（Tshebyshev Polynomial Interpolation）、納威（Neville）算法等［3]。以上三種方法中，拉格朗日多項(xiàng)式插值在展開級(jí)數(shù)較高的情況下容易產(chǎn)生龍格振蕩且耗費(fèi)時(shí)間較多［4]；切比雪夫多項(xiàng)式擬合雖然可以避免擬合區(qū)間兩端產(chǎn)生數(shù)據(jù)振蕩和跳躍，但其數(shù)學(xué)模型相對(duì)復(fù)雜計(jì)算量較大［5]；納維爾算法是一種逐次線性插值方法，數(shù)學(xué)模型簡單計(jì)算量小，內(nèi)插區(qū)間兩端無數(shù)據(jù)跳躍和振蕩現(xiàn)象。根據(jù)相關(guān)成果表明［5]：上述三種方法均能達(dá)到理想精度。

采用納維爾算法對(duì)精密星歷三維坐標(biāo)數(shù)據(jù)和鐘差改正數(shù)據(jù)進(jìn)行加密，將坐標(biāo)數(shù)據(jù)與廣播星歷計(jì)算數(shù)據(jù)進(jìn)行比對(duì)判斷精度，將鐘差加密結(jié)果與IGS提供的30s鐘差信息進(jìn)行比較獲取精度。

1 Neville算法插值原理

Neville算法作是一種逐次線性插值方法。常用的拉格朗日多項(xiàng)式有一個(gè)缺點(diǎn)，即在計(jì)算函數(shù)近似值精度不夠需要增加插值節(jié)點(diǎn)導(dǎo)致原有數(shù)據(jù)不能利用，而且一般的多項(xiàng)式插值不容易估計(jì)誤差，所以，只能靠經(jīng)驗(yàn)估計(jì)確定插值多項(xiàng)式的次數(shù)。Neville算法則能有效地克服上述缺點(diǎn)。

Neville算法進(jìn)行插值計(jì)算的遞推公式為

通過表1可以簡單明了地展現(xiàn)Neville算法的計(jì)算過程［6]。

表1 Neville算法的計(jì)算過程

不難看出，需要增加一個(gè)節(jié)點(diǎn)時(shí)就計(jì)算一行，I0，1，I0，1，2，I0，1，2，3，I0，1，2，3，4依次是1次到4次插值多項(xiàng)式的值。如果插值精度不滿足要求，只需增加一個(gè)節(jié)點(diǎn)即可，前面的計(jì)算結(jié)果完全有效。并可以逐次生成插值結(jié)果，每做一步檢查一次計(jì)算結(jié)果的精度，如不滿足精度要求，則繼續(xù)增加節(jié)點(diǎn)，直到滿足精度要求為止。

2 Neville算法的編程實(shí)現(xiàn)及精度驗(yàn)證

2.1 Neville算法的編程實(shí)現(xiàn)

由于精密星歷的數(shù)據(jù)格式比較固定，包括頭文件和一天之中的96個(gè)歷元數(shù)據(jù)，且每一顆衛(wèi)星的數(shù)據(jù)包含坐標(biāo)信息及鐘差改正信息，容易讀成矩陣數(shù)據(jù)。根據(jù)這一特點(diǎn)可采用Matlab進(jìn)行編程實(shí)現(xiàn)精密星歷的加密。

決定函數(shù) Neville（x1，y1，x）算法精度的是它的階數(shù)，在本程序中就是選取的節(jié)點(diǎn)時(shí)間段數(shù)x1。通過實(shí)驗(yàn)數(shù)據(jù)反復(fù)計(jì)算可知，當(dāng)階數(shù)為11時(shí)，Neville算法計(jì)算的結(jié)果精度達(dá)到最佳；y1為與x1相對(duì)應(yīng)的節(jié)點(diǎn)數(shù)據(jù)，可以是精密星歷提供的衛(wèi)星坐標(biāo)數(shù)據(jù)或鐘差改正數(shù)據(jù)；x為需要加密的時(shí)間段。在本程序中，為了匹配精密星歷數(shù)據(jù)15min采樣率的特點(diǎn)，將加密的區(qū)間固定為15min，這樣也能最大化利用精密星歷提供的數(shù)據(jù)。

2.2 衛(wèi)星坐標(biāo)加密精度驗(yàn)證

根據(jù)已有的文獻(xiàn)資料對(duì)于精密星歷衛(wèi)星坐標(biāo)加密的精度驗(yàn)證來看，基本上是采用幾種加密方法加密出與精密星歷提供的15min采樣率相一致的數(shù)據(jù)進(jìn)行比對(duì)［7]，得出相關(guān)結(jié)論。需要指出的是，如果和精密星歷提供的數(shù)據(jù)進(jìn)行精度對(duì)比，則需要剔除提供對(duì)比數(shù)據(jù)的該歷元，計(jì)算時(shí)取30min間隔的星歷作為內(nèi)插點(diǎn)內(nèi)插計(jì)算各個(gè)內(nèi)插時(shí)間段中間時(shí)刻的衛(wèi)星位置［8]，這樣就無法充分利用精密星歷的數(shù)據(jù)。

對(duì)此，通過對(duì)廣播星歷計(jì)算衛(wèi)星坐標(biāo)和精密星歷提供的整點(diǎn)坐標(biāo)之間的關(guān)系研究，得出了一種驗(yàn)證Neville算法加密結(jié)果精度的方法。

對(duì)于如何通過廣播星歷計(jì)算瞬時(shí)衛(wèi)星坐標(biāo)，可以查閱文獻(xiàn)［1，9]，同樣采用 Matlab進(jìn)行編程實(shí)現(xiàn)。廣播星歷采用IGS網(wǎng)站提供的2010年6月16日的文件brdc1670.10n，精密星歷則選用同天的jpl15883.sp3文件。時(shí)間區(qū)間選取2010年6月16日8時(shí)至9時(shí)，表2即為選取5顆衛(wèi)星分別在每隔15min的節(jié)點(diǎn)上Δs的計(jì)算結(jié)果，其中Δx＝｜x廣－x精｜，Δy＝｜y廣－y精｜，Δz＝｜z廣－z精｜，Δs＝單位均為m.

由于Neville插值屬于線性插值，再由表2數(shù)據(jù)可以推測(cè)，Δs值與時(shí)間序列存在某種線性關(guān)系。利用Matlab的最小二乘法［10]對(duì)此相關(guān)數(shù)據(jù)進(jìn)行處理得到如圖1的結(jié)果。

表2 PRN2、7、9、15、22號(hào)衛(wèi)星的相關(guān)計(jì)算結(jié)果

根據(jù)表1提供的數(shù)據(jù)得出，當(dāng)采用4階多項(xiàng)式擬合時(shí)，殘差最小。利用廣播星歷計(jì)算的衛(wèi)星坐標(biāo)對(duì)精密星歷加密數(shù)據(jù)進(jìn)行精度驗(yàn)證的思想：利用已知5個(gè)節(jié)點(diǎn)數(shù)據(jù)進(jìn)行最小二乘擬合，選取最優(yōu)擬合多項(xiàng)式，再利用該多項(xiàng)式求出對(duì)應(yīng)加密間隔的值λ，并與Δs值進(jìn)行作差處理（需剔除五個(gè)節(jié)點(diǎn)數(shù)據(jù)）ε＝｜λ－Δs｜，通過分析差值數(shù)據(jù)ε的中誤差σ間接驗(yàn)證利用Neville方法加密坐標(biāo)數(shù)據(jù)的精度。

現(xiàn)對(duì)PRN2衛(wèi)星在8時(shí)0分0秒至9時(shí)0分0秒這一區(qū)間加密的衛(wèi)星坐標(biāo)數(shù)據(jù)按照上述方法進(jìn)行精度驗(yàn)證，結(jié)果如圖2所示。經(jīng)過計(jì)算各衛(wèi)星的差值數(shù)據(jù)ε的中誤差σ值參見表3.

圖1 五顆衛(wèi)星Δs值與時(shí)間擬合曲線及殘差

表3 各衛(wèi)星的差值數(shù)據(jù)ε的中誤差σ

由表3數(shù)據(jù)顯示最大中誤差為6.07mm，最小為1.71mm，通過多天多顆衛(wèi)星的大量同類實(shí)驗(yàn)數(shù)據(jù)計(jì)算表明，利用Neville算法進(jìn)行30s坐標(biāo)內(nèi)插的數(shù)據(jù)精度符合精密單點(diǎn)定位計(jì)算要求。

2.3 衛(wèi)星鐘差改正數(shù)精度驗(yàn)證

將一小時(shí)區(qū)間內(nèi)加密所得到的鐘差改正數(shù)據(jù)與IGS提供的30s鐘差改正數(shù)文件的數(shù)據(jù)進(jìn)行對(duì)比，結(jié)果如表4所示。

表4 加密鐘差與IGS提供鐘差的差值比較

表4中的數(shù)據(jù)顯示，加密鐘差數(shù)據(jù)與IGS提供的30s鐘差數(shù)據(jù)的最大差值未超過1ns，最大精度可達(dá)10－3納秒級(jí)。同樣對(duì)多天多歷元衛(wèi)星的鐘差加密數(shù)據(jù)進(jìn)行相同處理后發(fā)現(xiàn)：Neville算法加密30s鐘差的精度可以滿足精密單點(diǎn)定位的需求。

4 結(jié) 論

利用Neville算法編程實(shí)現(xiàn)精密星歷的衛(wèi)星坐標(biāo)和鐘差加密計(jì)算，根據(jù)Neville算法的線性插值特性，對(duì)已有加密坐標(biāo)精度的驗(yàn)證方法進(jìn)行改進(jìn)，結(jié)合廣播星歷計(jì)算瞬時(shí)衛(wèi)星坐標(biāo)，通過最小二乘擬合實(shí)現(xiàn)加密坐標(biāo)精度的驗(yàn)證，得到毫米級(jí)的精度驗(yàn)證結(jié)果。不足之處在于此種驗(yàn)證方法只適用于Neville算法，對(duì)于其他加密方法如拉格朗日內(nèi)插、切比雪夫多項(xiàng)式擬合等非線性方法并不適用。對(duì)于鐘差改正數(shù)的加密驗(yàn)證結(jié)果得到亞納秒級(jí)的精度，均符合精密單點(diǎn)定位計(jì)算的要求。

［1]李征航，黃勁松.GPS測(cè)量與數(shù)據(jù)處理［M].武漢：武漢大學(xué)出版社，2005.

［2]BOCK H.Efficient methods for determining precise orbits of low earth orbiters using the global positioning system ［J].Geodatisch－geophysikalische Arbeiten in der Schweiz，2003，65（18）：52－55.

［3]魏二虎，柴 華.GPS精密星歷插值方法的比較研究［J].全球定位系統(tǒng)，2006，31（5）：13－15.

［4]洪 櫻，歐吉坤，彭碧波，等.GPS衛(wèi)星精密星歷和鐘差三種內(nèi)插方法的比較［J].武漢大學(xué)學(xué)報(bào)，2006，6（31）：516－518.

［5]邱 蕾，廖遠(yuǎn)琴.基于IGS精密星歷的衛(wèi)星坐標(biāo)插值［J].測(cè)繪工程，2008，4（17）：15－18.

［6]李慶揚(yáng)，王能超，易大義，等.數(shù)值分析［M].武漢：華中科技大學(xué)出版社，2006.

［7]楊善婷，陳藝軍，鄭家柱，等.GPS衛(wèi)星星歷插值擬合方法研究［J].森林工程，2010，1（26）：80－83.

［8]任迎華.IGS精密星歷內(nèi)插GPS衛(wèi)星位置研究［J].礦山測(cè)量，2010（3）：54－56.

［9]丁文武，劉萬科，李征航.基于廣播星歷的GPS衛(wèi)星坐標(biāo)的多項(xiàng)式擬合算法［J].海洋測(cè)繪，2007，6（27）：4－8.

［10]呂喜明，李明遠(yuǎn).最小二乘曲線擬合的 MATLAB實(shí)現(xiàn)［J].內(nèi)蒙古民族大學(xué)學(xué)報(bào)，2009，3（24）：125－127.