基于整體最小二乘原理的GPS 軌道坐標(biāo)擬合

謝孟辛，張捍衛(wèi)，李鵬杰，趙東方

（1.河南理工大學(xué) 測(cè)繪與國(guó)土信息工程學(xué)院，河南 焦作 454003；2.河南理工大學(xué) 資源與環(huán)境學(xué)院，河南 焦作 454003）

0 引言

國(guó)際全球衛(wèi)星導(dǎo)航系統(tǒng)（global navigation satellite system,GNSS）服務(wù)組織（International GNSS Service,IGS）提供的全球定位系統(tǒng)（global positioning system,GPS）精密星歷，其采樣間隔為15 min。而實(shí)際觀測(cè)中，GPS 接收機(jī)的采樣率一般低于1 min。因此，用戶須對(duì)GPS 精密星歷進(jìn)行數(shù)值擬合，才能夠獲得任意時(shí)刻衛(wèi)星的軌道坐標(biāo)[1]。

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

針對(duì)目前常用的最小二乘擬合法在高階次精密星歷擬合時(shí)擬合誤差逐漸增大且趨于發(fā)散的問(wèn)題，本文提出一種更為穩(wěn)健的擬合方法即整體最小二乘擬合法，并以正交多項(xiàng)式和常規(guī)多項(xiàng)式為基函數(shù)，分別根據(jù)最小二乘原理和整體最小二乘原理對(duì)精密星歷數(shù)據(jù)進(jìn)行擬合。

1 最小二乘擬合

1.1 切比雪夫多項(xiàng)式擬合

在擬合弧段tΔ 內(nèi)，以n階切比雪夫多項(xiàng)式為基函數(shù)，擬合GPS 衛(wèi)星精密星歷。設(shè)t0為擬合弧段的初始?xì)v元，由切比雪夫多項(xiàng)式的性質(zhì)，使用轉(zhuǎn)換公式[11]為

將變量t歸化到區(qū)間τ∈[? 1,+1 ]上，衛(wèi)星坐標(biāo)可用切比雪夫多項(xiàng)式[13]表示為

式中：n為切比雪夫多項(xiàng)式的階數(shù)；和分別為衛(wèi)星X、Y和Z坐標(biāo)分量的切比雪夫多項(xiàng)式系數(shù)，上標(biāo)T表示切比雪夫多項(xiàng)式擬合系數(shù)。切比雪夫多項(xiàng)式Ti遞推公式[14]為

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

式中：m為擬合過(guò)程所用到的精密星歷歷元個(gè)數(shù)；為使用切比雪夫多項(xiàng)式擬合法計(jì)算的衛(wèi)星X坐標(biāo)與觀測(cè)值xi之差。

將式（4）寫(xiě)成矩陣形式為

根據(jù)最小二乘原理平差準(zhǔn)則VTPV=min，由式（5）可得法方程為

式中P為偽距殘差。

由于各衛(wèi)星軌道觀測(cè)坐標(biāo)（精密星歷）為等精度觀測(cè)值，所以P為單位權(quán)矩陣，因此有

直接求逆矩陣，解得C為

1.2 勒讓德?tīng)柖囗?xiàng)式擬合

與切比雪夫多項(xiàng)式擬合類似，將自變量t歸化處理后，衛(wèi)星軌道坐標(biāo)的勒讓德?tīng)枺↙egendre）多項(xiàng)式擬合函數(shù)為[6]：

式中：n1為勒讓德?tīng)柖囗?xiàng)式的階數(shù)；和分別為X、Y和Z坐標(biāo)分量的勒讓德?tīng)柖囗?xiàng)式系數(shù)，上標(biāo)P表示勒讓德?tīng)柖囗?xiàng)式擬合系數(shù)。勒讓德?tīng)柖囗?xiàng)式Pi遞推公式[16]為

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

1.3 常規(guī)多項(xiàng)式擬合

衛(wèi)星軌道坐標(biāo)擬合時(shí)不再使用正交多項(xiàng)式，而使用一組線性無(wú)關(guān)的常規(guī)多項(xiàng)式{t0,t1,…,tn}作為坐標(biāo)擬合的基函數(shù)。為保證數(shù)值計(jì)算精度，將自變量t正則化處理[1]，即

式中 [t0,t1,…,tm]為擬合過(guò)程用到的精密星歷對(duì)應(yīng)的歷元時(shí)刻；t∈[t0,t1,…,tm]；mean（）與std（）分別表示求均值和標(biāo)準(zhǔn)差。

衛(wèi)星軌道坐標(biāo)的常規(guī)多項(xiàng)式擬合函數(shù)為

式中：n2為常規(guī)多項(xiàng)式的階數(shù)；和分別為X、Y和Z坐標(biāo)分量的常規(guī)多項(xiàng)式系數(shù)，上標(biāo)Y表示常規(guī)多項(xiàng)式擬合系數(shù)。式（12）中常規(guī)多項(xiàng)式系數(shù)求解過(guò)程與切比雪夫多項(xiàng)式擬合法和勒讓德?tīng)柖囗?xiàng)式擬合法相同。

2 整體最小二乘擬合

上述依據(jù)最小二乘原理的切比雪夫多項(xiàng)式擬合法、勒讓德?tīng)柖囗?xiàng)式擬合法和常規(guī)多項(xiàng)式擬合法求解的衛(wèi)星軌道坐標(biāo)擬合函數(shù)都是在經(jīng)典的GM 模型下根據(jù)最小二乘原理求解而得；然而在實(shí)際情況中，更多的是觀測(cè)值向量x和系數(shù)矩陣T均存在誤差，G-M 模型是不適用的。由式（8）求得的擬合系數(shù)C并不是最優(yōu)無(wú)偏的，此時(shí)引入變量含誤差模型（errors-in-variables,EIV），此模型將觀測(cè)值向量與系數(shù)矩陣的誤差均引入誤差方程[17]。這里以切比雪夫多項(xiàng)式擬合的誤差方程為例進(jìn)行說(shuō)明，即可得

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

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

在變量含誤差模型（EIV）下，使用普通最小二乘法已無(wú)法解決問(wèn)題，需根據(jù)整體最小二乘原理求解最優(yōu)擬合系數(shù)C。

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

依據(jù)平差準(zhǔn)則式（16），文獻(xiàn)[18]導(dǎo)出的整體最小二乘SVD 算法如下：

1）對(duì)B進(jìn)行SVD 分解，即

其中：U為（m+1）×（m+1）維正交矩陣；V為（n+2）×（n+2）維正交矩陣；ui和vi為各自正交矩陣內(nèi)部?jī)蓛烧坏膯挝幌蛄?；S為（m+1）×（n+2）維對(duì)角矩陣；si為矩陣B的奇異值（且s0≥s1≥ …≥sn+1≥0）。矩陣U、S和V的維數(shù)是通過(guò)矩陣B的維數(shù)確定的。

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

式中vi,n+1表示矩陣V第n+1 列第i個(gè)元素。文獻(xiàn)[18-19]已進(jìn)行嚴(yán)密的數(shù)學(xué)推導(dǎo)，并證明上述奇異值分解算法求解的擬合系數(shù)C在EIV 模型下是最優(yōu)無(wú)偏的。

同理，整體最小二乘法的勒讓德?tīng)柖囗?xiàng)式與常規(guī)多項(xiàng)式擬合也依照上述奇異值分解算法進(jìn)行解算。

通過(guò)第1 節(jié)與第2 節(jié)中敘述的擬合算法均可求得衛(wèi)星軌道坐標(biāo)擬合函數(shù)，即可求出擬合弧段[t0,t0+Δt]內(nèi)任意時(shí)刻的衛(wèi)星坐標(biāo)。

3 實(shí)驗(yàn)與結(jié)果分析

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

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

表2 點(diǎn)位中誤差mp（擬合弧段為6 h） mm

式中：Δx、Δy和Δz為擬合的結(jié)果與對(duì)應(yīng)時(shí)刻IGS提供的精密星歷之間的差值；m為擬合過(guò)程中所用到的歷元數(shù)。在衛(wèi)星軌道坐標(biāo)擬合實(shí)驗(yàn)中，相同弧段同一算法擬合結(jié)果的中誤差mx、m y、m z和mp呈現(xiàn)出相同規(guī)律。由于篇幅限制，這里只列出mx和mp的結(jié)果。

從表中可看出：擬合弧段為6 h，根據(jù)最小二乘原理以切比雪夫多項(xiàng)式和勒讓德?tīng)柖囗?xiàng)式（正交多項(xiàng)式）為基函數(shù)擬合精密星歷時(shí)，在擬合階數(shù)小于等于20的區(qū)間內(nèi)，X坐標(biāo)分量擬合中誤差mx與點(diǎn)位中誤差mp隨擬合階數(shù)的增加而減小，擬合誤差在20 階時(shí)達(dá)到最?。磺斜妊┓蚨囗?xiàng)式擬合的X坐標(biāo)分量擬合的中誤差為0.090 61 mm，點(diǎn)位中誤差為0.219 35 mm，勒讓德?tīng)柖囗?xiàng)式擬合的X坐標(biāo)分量擬合的中誤差為0.082 72 mm，點(diǎn)位中誤差為0.215 81 mm。在擬合階數(shù)大于20 至最大擬合階數(shù)范圍，擬合中誤差隨著擬合階數(shù)的增加而增大，軌道坐標(biāo)擬合精度降低。擬合弧段為6 h，根據(jù)最小二乘原理以常規(guī)多項(xiàng)式為基函數(shù)擬合精密星歷時(shí)，在擬合階數(shù)小于等于12的區(qū)間內(nèi)，X坐標(biāo)分量擬合中誤差mx與點(diǎn)位中誤差mp隨擬合階數(shù)的增加而減小，軌道坐標(biāo)擬合精度在擬合階數(shù)為12 階時(shí)達(dá)到最高；X坐標(biāo)分量擬合中誤差為20.285 64 mm，點(diǎn)位中誤差為23.664 46 mm。當(dāng)擬合階數(shù)大于12 階時(shí)，隨著擬合階數(shù)的增加，軌道擬合中誤差急劇增大，趨于發(fā)散。

衛(wèi)星軌道擬合弧段為6 h，根據(jù)最小二乘原理使用切比雪夫多項(xiàng)式、勒讓德?tīng)柖囗?xiàng)式（正交多項(xiàng)式）和常規(guī)多項(xiàng)式進(jìn)行軌道坐標(biāo)擬合時(shí)，3 種擬合方法擬合結(jié)果的X方向的中誤差mx和點(diǎn)位中誤差mp的數(shù)值變化規(guī)律相似：在一定階數(shù)范圍內(nèi)，隨著擬合階數(shù)的增加，擬合結(jié)果的中誤差減小，軌道擬合精度升高；擬合階數(shù)超過(guò)此區(qū)間，隨著擬合階數(shù)增加，擬合結(jié)果的中誤差增大且趨于發(fā)散，軌道擬合精度降低。不同之處在于：相比于常規(guī)多項(xiàng)式擬合，正交多項(xiàng)式擬合在更高階次才出現(xiàn)誤差發(fā)散現(xiàn)象，軌道擬合結(jié)果較穩(wěn)定，且使用正交多項(xiàng)式進(jìn)行軌道坐標(biāo)擬合的精度要優(yōu)于常規(guī)多項(xiàng)式擬合。

擬合弧段為6 h，根據(jù)整體最小二乘原理以切比雪夫多項(xiàng)式、勒讓德?tīng)柖囗?xiàng)式和常規(guī)多項(xiàng)式為基函數(shù)擬合精密星歷，在所能達(dá)到的最高擬合階數(shù)范圍內(nèi)，擬合精度隨著擬合階數(shù)的增加而升高，在高階次擬合時(shí)能夠保證擬合結(jié)果的穩(wěn)定性，不會(huì)出現(xiàn)最小二乘擬合在高階次擬合時(shí)軌道擬合誤差增大并趨于發(fā)散的現(xiàn)象。并且根據(jù)整體最小二乘原理得到的擬合結(jié)果精度更高，特別是以常規(guī)多項(xiàng)式為基函數(shù)根據(jù)整體最小二乘原理進(jìn)行軌道坐標(biāo)擬合時(shí)，擬合結(jié)果精度得到了極大的改善。

如圖1 與圖2 所示，擬合弧段為12 h，根據(jù)最小二乘原理和整體最小二乘原理，以切比雪夫多項(xiàng)式和勒讓德?tīng)柖囗?xiàng)式（正交多項(xiàng)式）為基函數(shù)進(jìn)行衛(wèi)星軌道坐標(biāo)擬合，其擬合結(jié)果與6 h 擬合弧段擬合結(jié)果的精度變化規(guī)律類似：擬合弧段為12 h，根據(jù)最小二乘原理進(jìn)行軌道擬合時(shí)，2 種正交多項(xiàng)式擬合結(jié)果在一定階數(shù)范圍內(nèi)，隨著擬合階數(shù)的增加，擬合結(jié)果的中誤差減小，軌道擬合精度升高；擬合階數(shù)超過(guò)此區(qū)間，隨著擬合階數(shù)的增加，擬合結(jié)果的中誤差增大且趨于發(fā)散，擬合精度降低；根據(jù)整體最小二乘原理進(jìn)行軌道坐標(biāo)擬合，2 種正交多項(xiàng)式擬合在所能達(dá)到的最高擬合階數(shù)范圍內(nèi)，擬合精度隨著擬合階數(shù)的增加而升高，在高階次擬合時(shí)能夠保證擬合結(jié)果的穩(wěn)定性。

圖1 以切比雪夫多項(xiàng)式為基函數(shù)的擬合點(diǎn)位中誤差mp變化曲線（擬合弧段12 h）

圖2 以勒讓德?tīng)柖囗?xiàng)式為基函數(shù)的擬合點(diǎn)位中誤差mp變化曲線（擬合弧段12 h）

擬合弧段為12 h，根據(jù)最小二乘原理以常規(guī)多項(xiàng)式為基函數(shù)進(jìn)行軌道坐標(biāo)擬合，所能達(dá)到的最高擬合精度只有米級(jí)，如擬合階數(shù)為15 階時(shí)，擬合結(jié)果的點(diǎn)位中誤差mp最小，為2 548.578 48 mm；擬合階數(shù)大于15，擬合結(jié)果中誤差隨著擬合階數(shù)的增加而急劇增大且趨于發(fā)散。即使根據(jù)整體最小二乘原理，其在高階次軌道坐標(biāo)擬合時(shí)也無(wú)法保證擬合結(jié)果的穩(wěn)定性，軌道坐標(biāo)擬合中誤差在高階次擬合時(shí)增大且趨于發(fā)散，擬合結(jié)果精度降低，如圖3 所示。

圖3 以常規(guī)多項(xiàng)式為基函數(shù)的擬合點(diǎn)位中誤差mp變化曲線（擬合弧段12 h）

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

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

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

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

但整體最小二乘法也無(wú)法解決條件數(shù)極大的病態(tài)方程求解的問(wèn)題。當(dāng)軌道擬合弧段為12 h 時(shí)，對(duì)于常規(guī)多項(xiàng)式的整體最小二乘擬合，其在高階次軌道坐標(biāo)擬合時(shí)方程系數(shù)矩陣Y的條件數(shù)數(shù)量級(jí)達(dá)到1016以上，如表4 所示，此時(shí)根據(jù)整體最小二乘原理求得的擬合系數(shù)C含有較大的誤差，衛(wèi)星軌道坐標(biāo)擬合結(jié)果的誤差在高階次增大且趨于發(fā)散，擬合精度降低。正交多項(xiàng)式（切比雪夫多項(xiàng)式、勒讓德?tīng)柖囗?xiàng)式）無(wú)論在最小二乘原理還是在整體最小二乘原理下進(jìn)行軌道坐標(biāo)擬合，其系數(shù)矩陣條件數(shù)的數(shù)量級(jí)明顯小于常規(guī)多項(xiàng)式系數(shù)矩陣條件數(shù)的數(shù)量級(jí)，因此其軌道擬合精度與高階次擬合時(shí)擬合結(jié)果的穩(wěn)定性要優(yōu)于常規(guī)多項(xiàng)式，這正是目前多采用切比雪夫等正交多項(xiàng)式進(jìn)行軌道坐標(biāo)擬合的原因。

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

本文分別基于最小二乘原理和整體最小二乘原理，以正交多項(xiàng)式和常規(guī)多項(xiàng)式為基函數(shù)，對(duì)6 和12 h 弧段的GPS 精密星歷進(jìn)行擬合，并對(duì)擬合結(jié)果精度進(jìn)行分析。結(jié)果表明：

1）整體最小二乘擬合解決了最小二乘擬合在高階次精密星歷擬合時(shí)擬合誤差逐漸增大且趨于發(fā)散的問(wèn)題，在擬合結(jié)果的穩(wěn)定性和擬合精度方面均優(yōu)于最小二乘擬合；

2）在擬合基函數(shù)選擇方面，鑒于正交多項(xiàng)式（切比雪夫多項(xiàng)式、勒讓德?tīng)柖囗?xiàng)式）在擬合過(guò)程的矩陣計(jì)算中具有更小的條件數(shù)，以正交多項(xiàng)式為基函數(shù)進(jìn)行擬合，擬合結(jié)果的穩(wěn)定性和擬合精度更高。