勒讓德多項(xiàng)式擬合IGS 精密星歷的算法改進(jìn)

申少飛，雷偉偉，李振南

(河南理工大學(xué) 測(cè)繪與國(guó)土信息工程學(xué)院,河南 焦作 454003)

0 引言

在進(jìn)行GPS 精密單點(diǎn)定位(PPP)時(shí)，需要用到衛(wèi)星的精密在軌位置.目前有兩種獲取衛(wèi)星軌道坐標(biāo)的方法，分別是通過(guò)廣播星歷和事后發(fā)布的精密星歷獲得[1].廣播星歷是通過(guò)瞬時(shí)參數(shù)求得衛(wèi)星速度與位置，計(jì)算復(fù)雜且精度較低，精密星歷的衛(wèi)星位置和鐘差精度要高于廣播星歷求得的衛(wèi)星位置和鐘差精度兩個(gè)數(shù)量級(jí)[2].因此，廣播星歷求得的衛(wèi)星軌道坐標(biāo)精度在實(shí)際應(yīng)用中難以滿足高精度用戶的需求.精密星歷由國(guó)際GNSS 服務(wù)(IGS)提供，是由若干衛(wèi)星跟蹤站的觀測(cè)數(shù)據(jù)，經(jīng)事后處理算得的用于衛(wèi)星精密定位等使用的衛(wèi)星軌道信息，其精度可達(dá)5 cm 甚至更高.IGS 發(fā)布的精密星歷采樣數(shù)據(jù)間隔為15 min，而接收機(jī)的采樣率通常為30 s、15 s甚至更短的時(shí)間間隔[3]，因此需采用軌道逼近的算法求取時(shí)間間隔更短的衛(wèi)星軌道坐標(biāo).目前用于軌道逼近的方法主要分為插值法和擬合法[4].

常用的插值方法有拉格朗日插值法、牛頓插值法、Neville 插值法和三角插值[5-7].常用的擬合方法包括切比雪夫多項(xiàng)式擬合和勒讓德多項(xiàng)式擬合[8].目前最常用的插值方法是拉格朗日插值法，但是當(dāng)改變插值點(diǎn)的個(gè)數(shù)時(shí)，需要重新計(jì)算，且計(jì)算量大，為此引入了牛頓插值法和Neville 插值法.文獻(xiàn)[9]說(shuō)明拉格朗日插值法與Neville 插值法實(shí)質(zhì)是一致的，但當(dāng)階數(shù)增加時(shí)，Neville 插值法原有計(jì)算仍然可用，插值公式不需要全部建立.文獻(xiàn)[10]指出當(dāng)取合適的階數(shù)時(shí)，三種插值方法的插值精度都能達(dá)到毫米級(jí).當(dāng)衛(wèi)星插值的軌道弧段較長(zhǎng)，需要將公式展開到較高的階次，此時(shí)，拉格朗日插值容易在區(qū)間兩端出現(xiàn)龍格震蕩現(xiàn)象.為解決這個(gè)問(wèn)題，文獻(xiàn)[11]采用滑動(dòng)式拉格朗日插值方法對(duì)GPS 精密星歷進(jìn)行插值，插值精度較好，即使階數(shù)較高，也不會(huì)出現(xiàn)龍格現(xiàn)象.文獻(xiàn)[12]采用一種新方法，即三角插值，得到的插值精度和滑動(dòng)式拉格朗日插值相當(dāng).文獻(xiàn)[13]采用改進(jìn)的切比雪夫多項(xiàng)式擬合法，在進(jìn)行高階次的軌道擬合時(shí)，也能保持較高的精度和穩(wěn)定性.文獻(xiàn)[14]采用勒讓德多項(xiàng)式擬合法，并與拉格朗日插值法進(jìn)行比較，結(jié)果表明勒讓德多項(xiàng)式擬合能達(dá)到毫米級(jí)精度且不會(huì)出現(xiàn)龍格震蕩之類的現(xiàn)象.文獻(xiàn)[15]指出用擬合法在計(jì)算系數(shù)時(shí)需要對(duì)法方程進(jìn)行求逆，階數(shù)過(guò)高容易使矩陣出現(xiàn)奇異，使法方程變?yōu)椴B(tài)方程，系數(shù)的微小變動(dòng)都會(huì)對(duì)解產(chǎn)生較大的誤差.

常規(guī)的勒讓德多項(xiàng)式方法在高階求逆時(shí)會(huì)產(chǎn)生較大的誤差，穩(wěn)定性差.對(duì)此，本文使用改進(jìn)的勒讓德多項(xiàng)式擬合法求解，采用LU 分解(LU Decomposition)算法和奇異值分解(SVD)算法求解軌道擬合，與常規(guī)算法相比具有更高的精度和穩(wěn)定性.

1 改進(jìn)的勒讓德多項(xiàng)式擬合法數(shù)學(xué)模型

1.1 勒讓德多項(xiàng)式擬合法

設(shè)擬合區(qū)間的初始?xì)v元為t0，擬合長(zhǎng)度為 Δt.在時(shí)間段 [t0,t0+Δt] 內(nèi)，以n階勒讓德多項(xiàng)式為基函數(shù)擬合衛(wèi)星軌道坐標(biāo)時(shí)，使用轉(zhuǎn)換公式[14]

將變量t歸化到區(qū)間 τ∈[-1,1]，衛(wèi)星軌道坐標(biāo)的勒讓德多項(xiàng)式擬合函數(shù)為

式中：n為勒讓德多項(xiàng)式的階數(shù)；分別為衛(wèi)星軌道X、Y、Z坐標(biāo)分量的多項(xiàng)式系數(shù).其中，Pi遞推公式為[16-17]

1.2 常規(guī)算法求解多項(xiàng)式系數(shù)

以式(2)中GPS 衛(wèi)星坐標(biāo)X分量中的系數(shù)為例，由式(2)～(3)可求得GPS 衛(wèi)星擬合坐標(biāo).設(shè)IGS提供的GPS 衛(wèi)星X坐標(biāo)分量的精密星歷為Xi，則觀測(cè)值向量的誤差方程為[14]

將式(4)簡(jiǎn)化為矩陣形式

式中：m為計(jì)算過(guò)程中用到的IGS 精密星歷的個(gè)數(shù)；X(τi)為對(duì)應(yīng)歷元時(shí)刻的精密星歷.

根據(jù)最小二乘原理VTPV=min 可得

直接求逆，解得多項(xiàng)式系數(shù)矩陣C為

各歷元坐標(biāo)可視為是等權(quán)觀測(cè)值，所以P為單位矩陣，則求解多項(xiàng)式系數(shù)矩陣C可以寫為

將式(7)中求得的多項(xiàng)式系數(shù)矩陣代入式(2)中，即可求得GPS 衛(wèi)星軌道在X方向坐標(biāo)的擬合多項(xiàng)式.同理可求得GPS 衛(wèi)星在Y方向和Z方向的坐標(biāo)擬合多項(xiàng)式.

1.3 LU 分解法求解多項(xiàng)式系數(shù)

根據(jù)普通最小二乘原理得到法方程式(6)，P為單位矩陣，可化簡(jiǎn)為

設(shè)A=BTB，式(9)可寫為

用高斯消去法對(duì)式(10)求解，可用矩陣表示為[18]

式(11)中，把方陣A分解為一個(gè)下三角矩陣和一個(gè)上三角矩陣，稱為 LU分解.P為初等變換矩陣.

則方程(10)可以改寫為

令b=BTX，則

令UC=Y，則有LY=b.那么由此把原方程的求解變?yōu)榍蠼庀禂?shù)矩陣為三角陣的方程，很容易實(shí)現(xiàn).以上過(guò)程即，首先計(jì)算A的 LU分解，由LY=b可得Y，接著由UC=Y求得多項(xiàng)式系數(shù)向量C.

1.4 SVD 分解法求解多項(xiàng)式系數(shù)

通過(guò)以矩陣SVD 分解為基礎(chǔ)的Moore-Penrose偽逆矩陣法求解GPS 衛(wèi)星軌道坐標(biāo)擬合函數(shù)系數(shù)C，系數(shù)C為

式中，B+為矩陣B的Moore-Penrose 偽逆矩陣，其Moore-Penrose 偽逆矩陣求解方法如下.

1.4.1 對(duì)矩陣B進(jìn)行SVD 分解[18]

式中：

U為(m)×(m)維正交矩陣，V為(n+1)×(n+1)維正交矩陣，ui和vi為各自矩陣內(nèi)部?jī)蓛烧坏膯挝幌蛄?，S為(m)×(n+1)維對(duì)角矩陣，Si為矩陣的奇異值(并且S0≥S1≥···≥Sn≥0)，以上三個(gè)矩陣的維數(shù)是通過(guò)矩陣B的維數(shù)所確定的.

1.4.2B的Moore-Penrose 偽逆矩陣

式(17)中的S+為矩陣S的偽逆矩陣，是矩陣S的非零元素取倒數(shù)之后再轉(zhuǎn)置得到的，所以S+為(n+1)×(m)的矩陣，其形式為

至此根據(jù)式(15)和(17)即可求得GPS 衛(wèi)星軌道坐標(biāo)擬合函數(shù)系數(shù)C.將求得的系數(shù)C代入式(2)中，可求得任意時(shí)刻的衛(wèi)星坐標(biāo).

2 勒讓德多項(xiàng)式擬合改進(jìn)算法分析

本文從IGS 官網(wǎng)下載2019 年9 月29 日的衛(wèi)星精密星歷作為計(jì)算數(shù)據(jù)，歷元時(shí)刻為00:00:00—23:45:00.采樣間隔為15 min，以擬合時(shí)段3 h 為例，采樣點(diǎn)所對(duì)應(yīng)時(shí)刻為0 min、15 min、30 min、45 min、···、180 min.為了驗(yàn)證擬合精度，選取采樣點(diǎn)相同時(shí)刻的擬合點(diǎn)，運(yùn)用多項(xiàng)式求出軌道坐標(biāo)擬合結(jié)果，并與擬合點(diǎn)處歷元所對(duì)應(yīng)的衛(wèi)星坐標(biāo)進(jìn)行作差比較，求出誤差絕對(duì)值的最大值和誤差中誤差，并求出每種擬合方法的點(diǎn)位中誤差.

選取3 h、6 h、12 h 三個(gè)擬合時(shí)段進(jìn)行擬合分析，采樣初始時(shí)刻為00:00:00，對(duì)3 種算法的擬合結(jié)果進(jìn)行對(duì)比分析，運(yùn)算環(huán)境在MATLAB 9.8 下進(jìn)行.IGS 精密星歷提供的坐標(biāo)誤差一般小于5 cm，因此3 種算法的擬合精度至少要小于5 cm.表1、表2 中差值是指擬合結(jié)果與對(duì)應(yīng)歷元衛(wèi)星坐標(biāo)X方向差的絕對(duì)值的最大值；均方根(RMS)是衛(wèi)星坐標(biāo)X分量擬合結(jié)果的中誤差；點(diǎn)位中誤差是衛(wèi)星坐標(biāo)X、Y、Z方向中誤差和的平方根.

2.1 擬合時(shí)段為3 h 的勒讓德多項(xiàng)式擬合精度分析

表1 給出勒讓德多項(xiàng)式擬合時(shí)段為3 h 的X方向擬合誤差絕對(duì)值的最大值和中誤差，可以看出3 種擬合方法無(wú)論是誤差絕對(duì)值的最大值還是中誤差結(jié)果一致，都隨著階數(shù)的增大而減小.在第8 階時(shí)，3 種擬合方法誤差絕對(duì)值的最大值均為69.01 mm，中誤差均為44.52 mm.在第9 階時(shí)，誤差絕對(duì)值的最大值和中誤差精度都迅速提高，達(dá)到了毫米級(jí).隨著階數(shù)的增大，在11 階時(shí)誤差絕對(duì)值的最大值和中誤差達(dá)到最小，分別為0.08 mm 和0.04 mm.擬合時(shí)段為3 h的Y方向勒讓德多項(xiàng)式擬合誤差絕對(duì)值的最大值與中誤差和Z方向勒讓德多項(xiàng)式擬合誤差絕對(duì)值的最大值與中誤差，與X方向的規(guī)律大致一樣.圖1 給出了3 種擬合方法的點(diǎn)位中誤差變化曲線，和表1 一樣，3 種擬合方法的精度變化完全一樣，都隨著擬合階數(shù)的增大而不斷提高.這說(shuō)明在擬合時(shí)段較小，階數(shù)較低時(shí)，3 種擬合方法的精度完全一樣.

表1 X 方向勒讓德多項(xiàng)式擬合(擬合時(shí)段3 h) mm

圖1 點(diǎn)位中誤差變化曲線(擬合時(shí)段3 h)

2.2 擬合時(shí)段為6 h 的勒讓德多項(xiàng)式擬合精度分析

由表2 可知，在擬合時(shí)段為6 h 的勒讓德多項(xiàng)式擬合中，采用常規(guī)算法的誤差絕對(duì)值的最大值與中誤差隨著階數(shù)的不斷升高呈現(xiàn)先減小后增大的趨勢(shì).在10～21 階時(shí)，誤差絕對(duì)值的最大值與中誤差不斷減小，在第10 階時(shí)分別為402.82 mm 和205.97 mm，誤差精度較低.隨著階數(shù)的升高，在第21 階時(shí)精度達(dá)到最高，分別為0.27 mm 和0.15 mm，精度達(dá)到亞毫米級(jí).接著隨著階數(shù)的升高精度不斷降低，在23 階時(shí)，常規(guī)算法的誤差絕對(duì)值的最大值和中誤差分別為319.34 mm 和102.68 mm，精度達(dá)到了分米級(jí)，達(dá)不到衛(wèi)星軌道坐標(biāo)的精度要求.LU 分解法和SVD 分解法的誤差絕對(duì)值的最大值與中誤差則隨著階數(shù)的升高不斷減小，且精度幾乎一致，在第10 階時(shí)，分別為402.82 mm 和205.97 mm.隨著階數(shù)的不斷增大，兩種算法的精度不斷提高，在第23 階時(shí)精度達(dá)到了最高，LU 分解法誤差絕對(duì)值的最大值與中誤差分別為0.21 mm 和0.08 mm，SVD 分解法誤差絕對(duì)值的最大值與中誤差分別為0.20 mm 和0.08 mm，精度均較高.可以看出在擬合階數(shù)較高時(shí)，LU 分解法與SVD分解法明顯比常規(guī)算法精度高.擬合時(shí)段為6 h 的Y方向勒讓德多項(xiàng)式擬合誤差絕對(duì)值的最大值與中誤差和Z方向勒讓德多項(xiàng)式擬合誤差絕對(duì)值的最大值與中誤差，與X方向的規(guī)律大致一樣.圖2 給出了3 種擬合方法的點(diǎn)位中誤差變化曲線，其結(jié)果和表2中3 種方法的規(guī)律相似，也進(jìn)一步驗(yàn)證了LU 分解法與SVD 分解法在高階時(shí)比常規(guī)算法更有優(yōu)勢(shì).

表2 X 方向勒讓德多項(xiàng)式擬合(擬合時(shí)段6 h) mm

圖2 點(diǎn)位中誤差變化曲線(擬合時(shí)段6 h)

2.3 擬合時(shí)段為12 h 的勒讓德多項(xiàng)式擬合精度分析

由圖3～4 可知，當(dāng)擬合時(shí)段為12 h，擬合階數(shù)較低時(shí)，3 種擬合方法在X方向誤差絕對(duì)值最大值和中誤差變化曲線一致，隨著階數(shù)的增大而不斷減小，趨近于0.從34 階開始，常規(guī)算法誤差絕對(duì)值的最大值和中誤差則隨著擬合階數(shù)的增大精度不斷降低，LU分解法和SVD 分解法的變化曲線則繼續(xù)趨近于0.從42 階開始，LU 分解法的精度開始降低，曲線趨于發(fā)散.而SVD 分解法不會(huì)出現(xiàn)類似的情況，隨著階數(shù)的增大，誤差絕對(duì)值的最大值和中誤差不斷減小，擬合精度越來(lái)越高，在擬合的最大階數(shù)處，精度最高，擬合曲線趨于收斂.圖5 是擬合時(shí)段為12 h 的點(diǎn)位中誤差，擬合曲線變化趨勢(shì)和圖3～4 大致一樣.從3 張圖中可以看出，擬合精度最高最穩(wěn)定的是SVD 分解法.

圖3 X 方向誤差最大值變化曲線(擬合時(shí)段12 h)

圖4 X 方向中誤差變化曲線(擬合時(shí)段12 h)

圖5 點(diǎn)位中誤差變化曲線(擬合時(shí)段12 h)

2.4 從條件數(shù)角度對(duì)結(jié)果進(jìn)行分析

在擬合時(shí)段為3 h、6 h 和12 h 時(shí)，3 種擬合方法的擬合精度并不完全一樣.在擬合時(shí)段為3 h 時(shí)，3 種擬合方法的精度完全一致.在擬合時(shí)段為6 h 時(shí)，常規(guī)算法的誤差在階數(shù)較高時(shí)開始增大，而另外兩種擬合方法則能保持較高的精度.在擬合階數(shù)較高，即擬合時(shí)段為12 h 時(shí)，LU 分解法的精度在高階時(shí)也開始降低，SVD 分解法無(wú)論低階還是高階都能保持較高的精度.這是因?yàn)樵跀M合時(shí)段較短，即3 h 時(shí)，法方程系數(shù)矩陣BTB維數(shù)較低，不易成為奇異矩陣，因此用常規(guī)算法對(duì)BTB求逆時(shí)，精度較高，與LU 分解法和SVD 分解法的精度相當(dāng).當(dāng)擬合時(shí)段為6 h 時(shí)，法方程系數(shù)矩陣BTB維數(shù)超過(guò)了20，此時(shí)矩陣很容易成為奇異矩陣，因此再用常規(guī)算法進(jìn)行多項(xiàng)式系數(shù)求解時(shí)，很容易出現(xiàn)較大的誤差.采用LU 分解法和SVD分解法不需要對(duì)矩陣求逆，而是分別對(duì)BTB和B矩陣進(jìn)行分解，因此兩者精度都遠(yuǎn)遠(yuǎn)高于用常規(guī)算法求得的系數(shù)C的精度.

在擬合時(shí)段為12 h 時(shí)，此時(shí)擬合階數(shù)較高，矩陣B和BTB的維數(shù)變得很大，因此在擬合過(guò)程中除了要解決滿秩矩陣，還要解決病態(tài)矩陣的問(wèn)題.由于在計(jì)算病態(tài)性矩陣方程的解時(shí)誤差幾乎是不可避免的，所以在可能的時(shí)候識(shí)別并避免病態(tài)性的矩陣是重要的，一般以條件數(shù)來(lái)衡量矩陣的病態(tài)性.由表3 可知，隨著階數(shù)的升高，矩陣B和BTB的條件數(shù)越來(lái)越大.在第40 階時(shí)，BTB的條件數(shù)已經(jīng)達(dá)到了1013，按照雙精度，我們會(huì)得到16-13=3 位正確數(shù)字的解C[18].LU 分解法是通過(guò)對(duì)矩陣BTB進(jìn)行LU 分解來(lái)求得多項(xiàng)式系數(shù)C，因此擬合階數(shù)過(guò)高，精度就會(huì)出現(xiàn)一定程度的下降.而SVD 分解法實(shí)則是對(duì)矩陣B進(jìn)行SVD 分解來(lái)求得系數(shù)C，所以即使在第47 階時(shí)，條件數(shù)也只是達(dá)到了1011，還能保持較高的精度.因此用SVD 分解法無(wú)論是低階還是高階都能保持較高的精度.

表3 矩陣條件數(shù)(擬合時(shí)段12 h)

3 結(jié)束語(yǔ)

計(jì)算表明，3 種勒讓德多項(xiàng)式擬合衛(wèi)星軌道坐標(biāo)的求解方法在不同擬合時(shí)段和擬合階數(shù)精度各不相同.在擬合時(shí)段為3 h，即擬合階數(shù)較低時(shí)，3 種擬合方法精度相當(dāng).在擬合時(shí)段為6 h 時(shí)，由于擬合階數(shù)較高，常規(guī)算法已不能滿足衛(wèi)星軌道坐標(biāo)的精度，改進(jìn)的LU 分解法和SVD 分解法對(duì)奇異矩陣都有較好的解決方法，因此精度相當(dāng).在擬合時(shí)段為12 h 時(shí)，擬合階數(shù)進(jìn)一步提高，此時(shí)在擬合過(guò)程中不僅存在奇異矩陣，還存在病態(tài)矩陣的問(wèn)題.SVD 分解法對(duì)兩種問(wèn)題都能較好的解決，所以在高階擬合時(shí)，無(wú)論是在精度還是穩(wěn)定性方面都要優(yōu)于LU 分解法和常規(guī)算法，具有較好的優(yōu)勢(shì).數(shù)字?jǐn)M合精度越高，表明擬合求得的多項(xiàng)式越逼近衛(wèi)星軌道坐標(biāo)，則對(duì)原星歷的精度影響越小.因此以后在采用勒讓德多項(xiàng)式對(duì)GPS 衛(wèi)星軌道坐標(biāo)進(jìn)行擬合時(shí)，可以優(yōu)先采用SVD 分解法進(jìn)行軌道坐標(biāo)擬合.