基于整體最小二乘原理的GPS 軌道坐標擬合

謝孟辛，張捍衛，李鵬杰，趙東方

（1.河南理工大學 測繪與國土信息工程學院，河南 焦作 454003；2.河南理工大學 資源與環境學院，河南 焦作 454003）

0 引言

國際全球衛星導航系統（global navigation satellite system,GNSS）服務組織（International GNSS Service,IGS）提供的全球定位系統（global positioning system,GPS）精密星歷，其采樣間隔為15 min。而實際觀測中，GPS 接收機的采樣率一般低于1 min。因此，用戶須對GPS 精密星歷進行數值擬合，才能夠獲得任意時刻衛星的軌道坐標[1]。

國內外學者對基于GPS 精密星歷的衛星軌道逼近方法有很多論述，主要分為插值法和擬合法2 個大類[2]。常用的插值方法為拉格朗日（Lagrange）插值，因其數學模型簡單、插值公式易于求得，是最經典最常用的衛星軌道坐標插值方法[3-5]。當用戶需要插值較長弧段時，根據拉格朗日插值定理，須將插值多項式展開到較高的階次；但在此情況下，拉格朗日插值在插值區間端點附近將產生龍格現象，即在插值區間端點附近根據插值函數計算的插值坐標含有較大誤差[6]。為解決這個問題，文獻[7-8]分別使用滑動插值法和切比雪夫（Chebyushev）多項式插值法：前者鑒于長弧段插值時在插值區間端點將產生龍格現象，而使插值點始終位于插值區間的中心；后者依據切比雪夫多項式定理通過重新選取插值節點使插值誤差的最大值最小，借此改善插值區間端點附近衛星軌道坐標計算不準確的問題。目前常用的擬合方法是根據最小二乘原理，以切比雪夫多項式為基函數，對精密星歷進行擬合，此擬合方法不會在擬合區間端點附近出現龍格現象，在擬合弧段長度和擬合精度方面也優于拉格朗日多項式插值[9-12]。因此，現在一般多采用切比雪夫多項式最小二乘擬合法對精密星歷進行數值擬合。文獻[13-14]利用此擬合方法進行具體應用，得到較好的擬合結果。但根據最小二乘原理以切比雪夫多項式為基函數對精密星歷進行數值擬合時，隨著擬合弧段增長和擬合階數增加，法方程將變為病態方程[15]，對精密星歷擬合函數系數的求解產生極大擾動，導致精密星歷擬合結果誤差增大且趨于發散，衛星軌道坐標擬合精度降低。

針對目前常用的最小二乘擬合法在高階次精密星歷擬合時擬合誤差逐漸增大且趨于發散的問題，本文提出一種更為穩健的擬合方法即整體最小二乘擬合法，并以正交多項式和常規多項式為基函數，分別根據最小二乘原理和整體最小二乘原理對精密星歷數據進行擬合。

1 最小二乘擬合

1.1 切比雪夫多項式擬合

在擬合弧段tΔ 內，以n階切比雪夫多項式為基函數，擬合GPS 衛星精密星歷。設t0為擬合弧段的初始歷元，由切比雪夫多項式的性質，使用轉換公式[11]為

將變量t歸化到區間τ∈[? 1,+1 ]上，衛星坐標可用切比雪夫多項式[13]表示為

式中：n為切比雪夫多項式的階數；和分別為衛星X、Y和Z坐標分量的切比雪夫多項式系數，上標T表示切比雪夫多項式擬合系數。切比雪夫多項式Ti遞推公式[14]為

以求解式（2）X坐標分量的切比雪夫多項式系數為例，設xi為觀測值，即IGS 提供的GPS精密星歷坐標數據，目前在求解最優擬合系數時所列誤差方程是在經典的高斯-馬爾可夫（Gause-Markov,G-M）模型下建立的，認為只有觀測值向量含有誤差，誤差方程[13]為

式中：m為擬合過程所用到的精密星歷歷元個數；為使用切比雪夫多項式擬合法計算的衛星X坐標與觀測值xi之差。

將式（4）寫成矩陣形式為

根據最小二乘原理平差準則VTPV=min，由式（5）可得法方程為

式中P為偽距殘差。

由于各衛星軌道觀測坐標（精密星歷）為等精度觀測值，所以P為單位權矩陣，因此有

直接求逆矩陣，解得C為

1.2 勒讓德爾多項式擬合

與切比雪夫多項式擬合類似，將自變量t歸化處理后，衛星軌道坐標的勒讓德爾（Legendre）多項式擬合函數為[6]：

式中：n1為勒讓德爾多項式的階數；和分別為X、Y和Z坐標分量的勒讓德爾多項式系數，上標P表示勒讓德爾多項式擬合系數。勒讓德爾多項式Pi遞推公式[16]為

式（9）中多項式系數同樣是在G-M 模型下根據最小二乘原理求解，求解算法與切比雪夫多項式擬合法相同。

1.3 常規多項式擬合

衛星軌道坐標擬合時不再使用正交多項式，而使用一組線性無關的常規多項式{t0,t1,…,tn}作為坐標擬合的基函數。為保證數值計算精度，將自變量t正則化處理[1]，即

式中 [t0,t1,…,tm]為擬合過程用到的精密星歷對應的歷元時刻；t∈[t0,t1,…,tm]；mean（）與std（）分別表示求均值和標準差。

衛星軌道坐標的常規多項式擬合函數為

式中：n2為常規多項式的階數；和分別為X、Y和Z坐標分量的常規多項式系數，上標Y表示常規多項式擬合系數。式（12）中常規多項式系數求解過程與切比雪夫多項式擬合法和勒讓德爾多項式擬合法相同。

2 整體最小二乘擬合

上述依據最小二乘原理的切比雪夫多項式擬合法、勒讓德爾多項式擬合法和常規多項式擬合法求解的衛星軌道坐標擬合函數都是在經典的GM 模型下根據最小二乘原理求解而得；然而在實際情況中，更多的是觀測值向量x和系數矩陣T均存在誤差，G-M 模型是不適用的。由式（8）求得的擬合系數C并不是最優無偏的，此時引入變量含誤差模型（errors-in-variables,EIV），此模型將觀測值向量與系數矩陣的誤差均引入誤差方程[17]。這里以切比雪夫多項式擬合的誤差方程為例進行說明，即可得

式中E為系數矩陣T的誤差矩陣，維數與T相同。將式（13）化為為

式中I為單位向量。設B=(T x)，D=(EV)，Z=(C?I)T，誤差方程式（13）可簡寫為

在變量含誤差模型（EIV）下，使用普通最小二乘法已無法解決問題，需根據整體最小二乘原理求解最優擬合系數C。

本文使用奇異值分解（singular value decomposition,SVD）算法求解整體最小二乘問題，SVD 算法的平差準則為

依據平差準則式（16），文獻[18]導出的整體最小二乘SVD 算法如下：

1）對B進行SVD 分解，即

其中：U為（m+1）×（m+1）維正交矩陣；V為（n+2）×（n+2）維正交矩陣；ui和vi為各自正交矩陣內部兩兩正交的單位向量；S為（m+1）×（n+2）維對角矩陣；si為矩陣B的奇異值（且s0≥s1≥ …≥sn+1≥0）。矩陣U、S和V的維數是通過矩陣B的維數確定的。

2）擬合系數C的整體最小二乘解（s n>sn+1）為

式中vi,n+1表示矩陣V第n+1 列第i個元素。文獻[18-19]已進行嚴密的數學推導，并證明上述奇異值分解算法求解的擬合系數C在EIV 模型下是最優無偏的。

同理，整體最小二乘法的勒讓德爾多項式與常規多項式擬合也依照上述奇異值分解算法進行解算。

通過第1 節與第2 節中敘述的擬合算法均可求得衛星軌道坐標擬合函數，即可求出擬合弧段[t0,t0+Δt]內任意時刻的衛星坐標。

3 實驗與結果分析

為分析上述各種擬合方法的衛星軌道坐標擬合結果，本文使用IGS 提供的2019-10-05的精密星歷數據進行實驗分析。實驗中GPS 衛星軌道擬合弧段為6 和12 h，采樣間隔為15 min，擬合弧段初始時刻為2019-10-05 T 00:00:00。將以上擬合算法擬合結果的中誤差列于表1、表2 和圖1～圖3中，其中mx、my、mz和mp分別為GPS 衛星X、Y和Z坐標分量擬合中誤差和點位中誤差，單位為mm，其計算公式為

表1 X 方向擬合中誤差m x（擬合弧段為6 h） mm

表2 點位中誤差mp（擬合弧段為6 h） mm

式中：Δx、Δy和Δz為擬合的結果與對應時刻IGS提供的精密星歷之間的差值；m為擬合過程中所用到的歷元數。在衛星軌道坐標擬合實驗中，相同弧段同一算法擬合結果的中誤差mx、m y、m z和mp呈現出相同規律。由于篇幅限制，這里只列出mx和mp的結果。

從表中可看出：擬合弧段為6 h，根據最小二乘原理以切比雪夫多項式和勒讓德爾多項式（正交多項式）為基函數擬合精密星歷時，在擬合階數小于等于20的區間內，X坐標分量擬合中誤差mx與點位中誤差mp隨擬合階數的增加而減小，擬合誤差在20 階時達到最小；切比雪夫多項式擬合的X坐標分量擬合的中誤差為0.090 61 mm，點位中誤差為0.219 35 mm，勒讓德爾多項式擬合的X坐標分量擬合的中誤差為0.082 72 mm，點位中誤差為0.215 81 mm。在擬合階數大于20 至最大擬合階數范圍，擬合中誤差隨著擬合階數的增加而增大，軌道坐標擬合精度降低。擬合弧段為6 h，根據最小二乘原理以常規多項式為基函數擬合精密星歷時，在擬合階數小于等于12的區間內，X坐標分量擬合中誤差mx與點位中誤差mp隨擬合階數的增加而減小，軌道坐標擬合精度在擬合階數為12 階時達到最高；X坐標分量擬合中誤差為20.285 64 mm，點位中誤差為23.664 46 mm。當擬合階數大于12 階時，隨著擬合階數的增加，軌道擬合中誤差急劇增大，趨于發散。

衛星軌道擬合弧段為6 h，根據最小二乘原理使用切比雪夫多項式、勒讓德爾多項式（正交多項式）和常規多項式進行軌道坐標擬合時，3 種擬合方法擬合結果的X方向的中誤差mx和點位中誤差mp的數值變化規律相似：在一定階數范圍內，隨著擬合階數的增加，擬合結果的中誤差減小，軌道擬合精度升高；擬合階數超過此區間，隨著擬合階數增加，擬合結果的中誤差增大且趨于發散，軌道擬合精度降低。不同之處在于：相比于常規多項式擬合，正交多項式擬合在更高階次才出現誤差發散現象，軌道擬合結果較穩定，且使用正交多項式進行軌道坐標擬合的精度要優于常規多項式擬合。

擬合弧段為6 h，根據整體最小二乘原理以切比雪夫多項式、勒讓德爾多項式和常規多項式為基函數擬合精密星歷，在所能達到的最高擬合階數范圍內，擬合精度隨著擬合階數的增加而升高，在高階次擬合時能夠保證擬合結果的穩定性，不會出現最小二乘擬合在高階次擬合時軌道擬合誤差增大并趨于發散的現象。并且根據整體最小二乘原理得到的擬合結果精度更高，特別是以常規多項式為基函數根據整體最小二乘原理進行軌道坐標擬合時，擬合結果精度得到了極大的改善。

如圖1 與圖2 所示，擬合弧段為12 h，根據最小二乘原理和整體最小二乘原理，以切比雪夫多項式和勒讓德爾多項式（正交多項式）為基函數進行衛星軌道坐標擬合，其擬合結果與6 h 擬合弧段擬合結果的精度變化規律類似：擬合弧段為12 h，根據最小二乘原理進行軌道擬合時，2 種正交多項式擬合結果在一定階數范圍內，隨著擬合階數的增加，擬合結果的中誤差減小，軌道擬合精度升高；擬合階數超過此區間，隨著擬合階數的增加，擬合結果的中誤差增大且趨于發散，擬合精度降低；根據整體最小二乘原理進行軌道坐標擬合，2 種正交多項式擬合在所能達到的最高擬合階數范圍內，擬合精度隨著擬合階數的增加而升高，在高階次擬合時能夠保證擬合結果的穩定性。

圖1 以切比雪夫多項式為基函數的擬合點位中誤差mp變化曲線（擬合弧段12 h）

圖2 以勒讓德爾多項式為基函數的擬合點位中誤差mp變化曲線（擬合弧段12 h）

擬合弧段為12 h，根據最小二乘原理以常規多項式為基函數進行軌道坐標擬合，所能達到的最高擬合精度只有米級，如擬合階數為15 階時，擬合結果的點位中誤差mp最小，為2 548.578 48 mm；擬合階數大于15，擬合結果中誤差隨著擬合階數的增加而急劇增大且趨于發散。即使根據整體最小二乘原理，其在高階次軌道坐標擬合時也無法保證擬合結果的穩定性，軌道坐標擬合中誤差在高階次擬合時增大且趨于發散，擬合結果精度降低，如圖3 所示。

圖3 以常規多項式為基函數的擬合點位中誤差mp變化曲線（擬合弧段12 h）

綜上所述，相比于最小二乘擬合，整體最小二乘擬合的軌道擬合精度更高，在高階次軌道擬合時擬合結果更穩定，其原因在于：在軌道擬合過程中，求解擬合函數系數C時會遇到病態方程求解問題，特別是最小二乘擬合；隨著擬合弧段增長，擬合階數增加，法方程式（7）TTTC=TTx系數矩陣TTT（PTP,YTY）的條件數（通常以條件數衡量矩陣的病態性）會急劇增加。如表3 與表4 所示。

表中：cond（TTT）、cond（PTP）和cond（YTY）分別表示切比雪夫多項式擬合、勒讓德爾多項式擬合和常規多項式擬合的法方程系數矩陣的條件數；cond（T）、cond（P）和cond（Y）分別表示切比雪夫多項式擬合、勒讓德爾多項式擬合和常規多項式擬合的誤差方程系數矩陣的條件數。法方程式（7）將變為嚴重的病態方程，求解病態法方程而得到的擬合函數系數C含有極大的誤差。針對解決病態方程求解不準確的問題，整體最小二乘擬合有2 個方面的優勢：①整體最小二乘擬合無須通過法方程而解得擬合系數C，其是對矩陣T的增廣矩陣(T x)進行SVD 分解，由分解而得右奇異矩陣V的元素求解擬合系數C，在此可只考慮矩陣T的病態性，由表3 與表4 可以看出，矩陣T比矩陣TTT的條件數數量級小得多，即使在高階次擬合時方程仍為良態方程；②整體最小二乘擬合是在EIV 模型下求解擬合函數系數C，計算時將觀測值向量誤差與系數矩陣誤差均納入誤差方程，依據平差準則式（16），使誤差方程中系數矩陣誤差與觀測值向量誤差最小，由此可使病態誤差方程中系數矩陣誤差與觀測值向量誤差對擬合系數C的求解影響降至最小，而最小二乘擬合只考慮觀測值含有誤差。

表3 法方程與誤差方程系數矩陣條件數（擬合弧段為6 h）

表4 法方程與誤差方程系數矩陣條件數（擬合弧段為12 h）

但整體最小二乘法也無法解決條件數極大的病態方程求解的問題。當軌道擬合弧段為12 h 時，對于常規多項式的整體最小二乘擬合，其在高階次軌道坐標擬合時方程系數矩陣Y的條件數數量級達到1016以上，如表4 所示，此時根據整體最小二乘原理求得的擬合系數C含有較大的誤差，衛星軌道坐標擬合結果的誤差在高階次增大且趨于發散，擬合精度降低。正交多項式（切比雪夫多項式、勒讓德爾多項式）無論在最小二乘原理還是在整體最小二乘原理下進行軌道坐標擬合，其系數矩陣條件數的數量級明顯小于常規多項式系數矩陣條件數的數量級，因此其軌道擬合精度與高階次擬合時擬合結果的穩定性要優于常規多項式，這正是目前多采用切比雪夫等正交多項式進行軌道坐標擬合的原因。

4 結束語

本文分別基于最小二乘原理和整體最小二乘原理，以正交多項式和常規多項式為基函數，對6 和12 h 弧段的GPS 精密星歷進行擬合，并對擬合結果精度進行分析。結果表明：

1）整體最小二乘擬合解決了最小二乘擬合在高階次精密星歷擬合時擬合誤差逐漸增大且趨于發散的問題，在擬合結果的穩定性和擬合精度方面均優于最小二乘擬合；

2）在擬合基函數選擇方面，鑒于正交多項式（切比雪夫多項式、勒讓德爾多項式）在擬合過程的矩陣計算中具有更小的條件數，以正交多項式為基函數進行擬合，擬合結果的穩定性和擬合精度更高。