申少飛，雷偉偉，李振南

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

0 引言

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

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

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

1 改進的勒讓德多項式擬合法數(shù)學(xué)模型

1.1 勒讓德多項式擬合法

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

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

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

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

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

將式(4)簡化為矩陣形式

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

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

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

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

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

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

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

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

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

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

則方程(10)可以改寫為

令b=BTX，則

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

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

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

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

1.4.1 對矩陣B進行SVD 分解[18]

式中：

U為(m)×(m)維正交矩陣，V為(n+1)×(n+1)維正交矩陣，ui和vi為各自矩陣內(nèi)部兩兩正交的單位向量，S為(m)×(n+1)維對角矩陣，Si為矩陣的奇異值(并且S0≥S1≥···≥Sn≥0)，以上三個矩陣的維數(shù)是通過矩陣B的維數(shù)所確定的.

1.4.2B的Moore-Penrose 偽逆矩陣

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

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

2 勒讓德多項式擬合改進算法分析

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

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

2.1 擬合時段為3 h 的勒讓德多項式擬合精度分析

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

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

圖1 點位中誤差變化曲線(擬合時段3 h)

2.2 擬合時段為6 h 的勒讓德多項式擬合精度分析

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

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

圖2 點位中誤差變化曲線(擬合時段6 h)

2.3 擬合時段為12 h 的勒讓德多項式擬合精度分析

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

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

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

圖5 點位中誤差變化曲線(擬合時段12 h)

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

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

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

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

3 結(jié)束語

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