孫華麗，張政治

（61920部隊(duì)，成都 610505）

1 引言

在全球定位系統(tǒng) （global positioning system，GPS）衛(wèi)星導(dǎo)航定位時(shí)，不管是采用差分模式還是非差精密單點(diǎn)定位，都是利用距離后方交會(huì)測(cè)量原理，即用戶通過(guò)測(cè)量接收機(jī)到衛(wèi)星間的距離來(lái)確定接收機(jī)的位置，于是導(dǎo)航定位的前提是實(shí)時(shí)獲取衛(wèi)星的位置和速度，而衛(wèi)星的位置和速度是由衛(wèi)星的廣播星歷參數(shù)來(lái)計(jì)算得到的［1］。

由接收機(jī)收到的導(dǎo)航電文中開(kāi)普勒軌道參數(shù)和軌道攝動(dòng)參數(shù)等信息按照固的公式計(jì)算得到GPS衛(wèi)星的瞬間坐標(biāo)，即所謂的直接法，文獻(xiàn)［2］詳細(xì)論述了其計(jì)算步驟。用廣播星歷直接計(jì)算GPS坐標(biāo)，計(jì)算量大，過(guò)程復(fù)雜，需要在星歷文件與觀測(cè)文件之間不停地進(jìn)行時(shí)間比對(duì)和衛(wèi)星比對(duì)，效率低。并且由于廣播星歷數(shù)據(jù)更新周期為2h，用相鄰兩組星歷計(jì)算同一時(shí)刻坐標(biāo)時(shí)會(huì)產(chǎn)生不同結(jié)果，這為臨界區(qū)域的坐標(biāo)確定、周跳剔除以及觀測(cè)值殘差分析等帶來(lái)諸多不便。在實(shí)際應(yīng)用中一般采用插值的方法獲取所需時(shí)刻的衛(wèi)星位置坐標(biāo)，因此有必要尋找一種計(jì)算過(guò)程簡(jiǎn)便、占用內(nèi)存小，且滿足精度要求的插值算法將特定時(shí)段衛(wèi)星軌道用某些具有代表性的坐標(biāo)點(diǎn)標(biāo)準(zhǔn)化。

2 三種常用插值算法

目前在利用廣播星歷計(jì)算GPS衛(wèi)星坐標(biāo)這一領(lǐng)域，最常用的插值方法是拉格朗日插值，拉格朗日插值模型比較簡(jiǎn)單，但當(dāng)需要進(jìn)行增減節(jié)點(diǎn)時(shí)，所有的插值基函數(shù)將改變，原有公式必須重新建立，計(jì)算量大。當(dāng)需要提高精度進(jìn)行高階插值時(shí)，可能會(huì)產(chǎn)生龍格現(xiàn)象，并且在插值區(qū)間端點(diǎn)處容易產(chǎn)生震蕩。三種插值算法簡(jiǎn)述為［3］

2.1 牛頓插值

牛頓插值將插值基函數(shù)定義成

其多項(xiàng)式的系數(shù)是各階均差，定義為：

一階均差定義為

二階均差定義為

類(lèi)似地，定義k階均差

根據(jù)均差定義，高階均差是低一階均差的均差。定義牛頓均差插值多項(xiàng)式

在實(shí)際計(jì)算中，當(dāng)增加節(jié)點(diǎn)時(shí)，只需增加一個(gè)被加項(xiàng)，原來(lái)計(jì)算全部有效。這樣就避免了增減節(jié)點(diǎn)引起插值基函數(shù)的改變，導(dǎo)致整個(gè)公式的變化。

2.2 埃爾米特插值

埃爾米特插值也叫指定微商值的插值，是一種光滑的插值方法，不但要求在節(jié)點(diǎn)上函數(shù)值相等，而且為了保證光滑程度還要求對(duì)應(yīng)的導(dǎo)數(shù)值也相等。首先構(gòu)造埃爾米特插值多項(xiàng)式，設(shè)已知節(jié)點(diǎn)

a≤x0＜x1…＜xn≤b上，

yj＝f（xj），mj＝f′（xj），（j＝0，1，…，n），這里給出了2n＋2個(gè)條件，可唯一確定一個(gè)次數(shù)不超過(guò)2n＋1的多項(xiàng)式

其形式為

寫(xiě)成用插值基函數(shù)表示的形式

數(shù)學(xué)上可以證明這種插值函數(shù)具有一致收斂性。

2.3 樣條插值

樣條插值能用較低次數(shù)的多項(xiàng)式達(dá)到較高階的光滑度，若S（x）滿足以下條件：

（1）在給定區(qū)間［a，b］上有二階連續(xù)導(dǎo)數(shù)；

（2）給定節(jié)點(diǎn)a＝x0＜x1＜…＜xn＝b，在每個(gè)小區(qū)間 ［xi，xi＋1］上，S（x）是三次多項(xiàng)式。

稱(chēng)滿足以上兩個(gè)條件的函數(shù)S（x）為節(jié)點(diǎn)x0，x1，…，xn的三次樣條函數(shù)。

若在節(jié)點(diǎn)xj上給定函數(shù)值

且S（x）＝y(tǒng)i，i＝ （0，1…，n）成立時(shí)，則稱(chēng)S（x）為節(jié)點(diǎn)x0，x1，…，xn的三次樣條插值函數(shù)，其中x0，x1，…，xn為樣條節(jié)點(diǎn)。三次樣條表達(dá)式為

Mi在力學(xué)上解釋為細(xì)梁在xi截面處的彎矩，稱(chēng)為S（x）的矩。文獻(xiàn) ［4］通過(guò)M，T關(guān)系式法構(gòu)造五次樣條插值函數(shù)。

3 三種算法插值計(jì)算GPS坐標(biāo)

3.1 數(shù)值實(shí)驗(yàn)

利用PRN1衛(wèi)星在2007年5月7日歷元時(shí)刻10時(shí)的廣播星歷參數(shù)文件，直接法可計(jì)算出9時(shí)至10時(shí)30分之間的相關(guān)坐標(biāo)值。本文等距選取插值節(jié)點(diǎn)時(shí)刻為9h：15m：10h30m作為訓(xùn)練集，記作節(jié)點(diǎn)0：15：90，這樣選取的歷元時(shí)刻是整分鐘，可以避免出現(xiàn)小數(shù)位，導(dǎo)致插值時(shí)出現(xiàn)浮點(diǎn)運(yùn)算，降低了計(jì)算復(fù)雜度，這也是采用偶數(shù)階插值的原因。等距選取測(cè)試歷元時(shí)刻為9h：5m：10h30m，記為節(jié)點(diǎn)0：5：90。

由于直接法得到的坐標(biāo)值是離散的，因此在埃爾米特插值和樣條插值中，補(bǔ)充導(dǎo)數(shù)條件不能由直接法給出，現(xiàn)有的研究表明8階拉格朗日插值誤差達(dá)到厘米級(jí)［5-8］，可以滿足GPS普通用戶的應(yīng)用需求，于是這里可以借助8階拉格朗日插值多項(xiàng)式對(duì)插值節(jié)點(diǎn)求各階導(dǎo)數(shù)作為補(bǔ)充條件。在進(jìn)行誤差統(tǒng)計(jì)分析時(shí)，由于正負(fù)誤差效果等同，均先將原始誤差值取絕對(duì)值，并且根據(jù)插值的定義，這里還要剔除誤差曲線插值節(jié)點(diǎn)處的零值。

在牛頓插值中，利用7個(gè)節(jié)點(diǎn)構(gòu)造6階牛頓均差插值。在測(cè)試集點(diǎn)0：5：90上，將6階牛頓均差插值廣播星歷所得坐標(biāo)與直接法計(jì)算出來(lái)的真值做比較，圖1給出了三個(gè)方向上的誤差分布曲線，表1列出了三個(gè)方向上的誤差統(tǒng)計(jì)情況。

圖1 6階牛頓插值誤差曲線

表1 6階牛頓插值誤差統(tǒng)計(jì)/m

同樣的方法，在節(jié)點(diǎn)0：15：90的基礎(chǔ)上，通過(guò)直接計(jì)算追加兩個(gè)節(jié)點(diǎn)，共9個(gè)歷元時(shí)刻作為訓(xùn)練集節(jié)點(diǎn)，構(gòu)造8階牛頓均差插值多項(xiàng)式。在測(cè)試集點(diǎn)0：5：90上，將8階牛頓均差插值廣播星歷所得坐標(biāo)與直接法所得結(jié)果做比較，圖2給出了三個(gè)方向上的誤差分布曲線，表2列出了三個(gè)方向上的誤差統(tǒng)計(jì)情況。

圖2 8階牛頓插值誤差曲線

表2 8階牛頓插值誤差統(tǒng)計(jì)/m

在埃爾米特插值計(jì)算中，本文考慮插值精度的要求和便于比較，筆者在訓(xùn)練集節(jié)點(diǎn)上分別構(gòu)造了三點(diǎn)5次和四點(diǎn)4次埃爾米特插值多項(xiàng)式，并應(yīng)用于衛(wèi)星軌道插值。在測(cè)試集點(diǎn)0：5：90上分別將分段5次和7次埃爾米特插值所得結(jié)果與直接法計(jì)算出來(lái)的真值做比較。圖3給出了5次埃爾米特插值在三個(gè)方向上的誤差分布曲線，表3列出了5次埃爾米特插值在三個(gè)方向上的誤差統(tǒng)計(jì)情況，圖4和表4給出了7次的情形。

圖3 5次埃爾米特插值誤差曲線

表3 5次埃爾米特插值誤差統(tǒng)計(jì)/m

圖4 7次埃爾米特插值誤差曲線

表4 7次埃爾米特插值誤差統(tǒng)計(jì)/m

在樣條插值中，筆者嘗試了3次樣條插值和5次樣條插值，在測(cè)試集點(diǎn)0：5：90上，將3次樣條和5次樣條插值結(jié)果與直接法計(jì)算出的真值做比較。圖5給出了3次樣條插值在3個(gè)方向上的誤差分布曲線，表5列出了3次樣條插值3個(gè)方向上的誤差統(tǒng)計(jì)情況，圖6和表6給出了5次的情形。

圖5 3次樣條插值誤差曲線

表5 3次樣條插值誤差統(tǒng)計(jì)/m

圖6 五次樣條插值誤差曲線

表6 五次樣條插值誤差統(tǒng)計(jì)/m

為了考察三種插值方法插值結(jié)果與直接法所得真值的離散程度，綜合考慮階數(shù)和精度的要求，下面選取6階牛頓插值、5次埃爾米特插值、五次樣條插值來(lái)計(jì)算在節(jié)點(diǎn)處的均方根 （root mean square，RMS），其公式是

在同一坐標(biāo)系中畫(huà)出三種方法的RMS分布圖如圖7所示。

圖7 三種插值方法的RMS分布圖

3.2 精度說(shuō)明

由圖1和表1可知，6階牛頓插值在Z方向上精度最高，達(dá)到厘米級(jí)，Y方向次之，為分米級(jí)，X方向僅為米級(jí)。由圖2和表2可知，8階牛頓插值在Y、Z方向上都達(dá)到了毫米級(jí)，X方向達(dá)到了厘米級(jí)。8階有小幅震蕩，6階的震蕩現(xiàn)象比較嚴(yán)重。三個(gè)方向上在插值區(qū)間中間部分精度都比端點(diǎn)出高一個(gè)量級(jí)，與插值的計(jì)算特點(diǎn)相符，下面埃爾米特插值和樣條插值的結(jié)果也有類(lèi)型的現(xiàn)象。

由圖3和表3可知，五次埃爾米特插值在X、Y、Z方向精度分別為分米、厘米、毫米。由圖4和表4可知，7次埃爾米特插值誤差在X方向?yàn)槔迕准?jí)，Y、Z方向達(dá)到毫米級(jí)。三個(gè)方向在插值區(qū)間左端點(diǎn)處仍有較明顯的震蕩現(xiàn)象。

由圖5和表5可知，3次樣條插值的計(jì)算精度在十米級(jí)，不能滿足用戶的精度需求，但是沒(méi)有震蕩現(xiàn)象發(fā)生。由圖6和表6可知，5次樣條插值精度在X、Y方向?yàn)榉置准?jí)，在Z方向?yàn)槔迕准?jí)，在插值區(qū)間端點(diǎn)處震蕩現(xiàn)象較不明顯。

由圖7三種插值方法的RMS分布圖知，區(qū)間內(nèi)精度較高，區(qū)間端點(diǎn)處會(huì)產(chǎn)生震蕩。牛頓插值誤差斜率變化比較大，震蕩現(xiàn)象最為嚴(yán)重，而隨著插值階數(shù)的提高，震蕩也會(huì)更嚴(yán)重。埃爾米特插值比較好的處理了區(qū)間端點(diǎn)處的震蕩問(wèn)題。

3.3 結(jié)果分析

拉格朗日多項(xiàng)式在增加新的節(jié)點(diǎn)時(shí)，原有基函數(shù)需要重新計(jì)算，而牛頓法通過(guò)構(gòu)造差商，在加入新節(jié)點(diǎn)后，只需計(jì)算一次更高階的差商，原有差商結(jié)果仍可繼續(xù)使用，解決了動(dòng)態(tài)增加節(jié)點(diǎn)的問(wèn)題，使迭代方式具有承襲性，程序?qū)崿F(xiàn)靈活，但其精度不高。

用分段低次埃爾米特插值可以做到降低插值階數(shù)情況下提高插值精度，這樣可以降低計(jì)算復(fù)雜度，降低了插值多項(xiàng)式階數(shù)還可以進(jìn)一步避免了龍格震蕩現(xiàn)象的發(fā)生。增加1個(gè)節(jié)點(diǎn)時(shí)，埃爾米特插值可以提高2個(gè)階數(shù)，并且由于采用分段低次處理插值問(wèn)題，可以較好的避免拉格朗日多項(xiàng)式在增加新的節(jié)點(diǎn)時(shí)，原有基函數(shù)需要重新計(jì)算的問(wèn)題，只需局部改變多項(xiàng)式基函數(shù)，模型變化小，易于程序移植。

5次樣條插值方法模型簡(jiǎn)單，精度高，可克服分段插值函數(shù)整體光滑性較差的缺點(diǎn)，同時(shí)計(jì)算復(fù)雜度比牛頓插值和埃爾米特插值降低了一個(gè)階［4］。5次樣條插值方法采用分段處理插值問(wèn)題，各段都是相對(duì)獨(dú)立的，并且不涉及插值基函數(shù)，插值函數(shù)表達(dá)式形式固定，只需解方程求出區(qū)間內(nèi)樣條節(jié)點(diǎn)的二階和四階導(dǎo)數(shù)值，較好的解決了動(dòng)態(tài)增減節(jié)點(diǎn)的問(wèn)題，程序?qū)崿F(xiàn)靈活。

4 結(jié)論

三種插值方法各有特點(diǎn)，牛頓插值承襲性較強(qiáng)，但是精度沒(méi)有得到明顯提高，適應(yīng)于衛(wèi)星有限、點(diǎn)位快速計(jì)算且無(wú)需高精度；埃爾米特插值精度高，在節(jié)點(diǎn)處具有一階光滑度、整體逼近效果好并且一致收斂，適用于較多衛(wèi)星位置計(jì)算或局部軌道計(jì)算；樣條插值在穩(wěn)定性、震蕩性方面有其獨(dú)特的優(yōu)勢(shì)，適合于軌道弧段邊界點(diǎn)的位置計(jì)算。在利用GPS廣播星歷對(duì)衛(wèi)星軌道進(jìn)行插值計(jì)算中，上述三種方法在一定程度上可以替代傳統(tǒng)的拉格朗日插值方法，在一定范圍內(nèi)值得推廣。

［1］魏子卿，葛茂榮.GPS相對(duì)定位的數(shù)學(xué)模型［M］.北京：測(cè)繪出版社，1998.

［2］HOFNM ANN－WELLENHOF B，LICHTENEGGER H，COLLINS J.GPS Theory and Practice.Wien：Springer－Verlag，1997.

［3］李慶揚(yáng)，王能超，易大義.數(shù)值分析［M］.4版.北京：清華大學(xué)出版社，2005.

［4］唐建國(guó)，蔣九林.五次樣條函數(shù)的構(gòu)造［J］.鄭州大學(xué)學(xué)報(bào)：理學(xué)版，2005，37（3）：15－18.

［5］江國(guó)焰，李明峰，朱振宇，等.GPS衛(wèi)星廣播星歷的Lagrange等距插值算法［J］.南京工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，2008，30（1）：34－38.

［6］萬(wàn)亞豪，張書(shū)畢，侯東陽(yáng).基于GPS廣播星歷的衛(wèi)星位置擬合精度分析［J］.測(cè)繪工程，2011，20（3）：36－40.

［7］陳劉成，韓春好，陳金平.廣播星歷參數(shù)擬合算法研究［J］.測(cè)繪科學(xué)，2007，32（3）：12－15.

［8］萬(wàn)亞豪，張書(shū)畢，侯東陽(yáng).基于GPS廣播星歷的衛(wèi)星位置擬合精度分析［J］.測(cè)繪工程，2011，20（3）：36－40.