精密定軌中三類經(jīng)驗(yàn)加速度的有效建模方法

王文彬，高 揚(yáng)

(中國科學(xué)院空間應(yīng)用工程與技術(shù)中心，北京100094)

0 引言

受太陽活動、地磁活動、季風(fēng)等復(fù)雜因素影響以及大氣密度、輻射光壓、重力場等模型誤差影響，很難建立地球低軌衛(wèi)星高精度確定性動力學(xué)模型。目前確定性動力學(xué)模型還無法匹配GPS載波相位的測量精度(毫米級)［1－2］，所以地球低軌衛(wèi)星若只采用確定性動力學(xué)定軌方法較難獲得高精度的定軌結(jié)果。簡化動力學(xué)定軌方法在確定性動力學(xué)模型的基礎(chǔ)上引入經(jīng)驗(yàn)加速度模型吸收未模型化動力學(xué)誤差，使得簡化動力學(xué)模型(包括確定性動力學(xué)模型和經(jīng)驗(yàn)加速度模型)精度能夠匹配GPS載波相位的測量精度，同時(shí)能夠保持動力學(xué)模型對衛(wèi)星軌道的約束作用，保證軌道解的穩(wěn)定性。

目前存在三類典型的經(jīng)驗(yàn)加速度模型:瞬時(shí)速度脈沖(Instantaneous Velocity Pulses，IVP)［2－5］、逐段常量經(jīng)驗(yàn)加速度(Piecewise Empirical Accelerations，PEA)［1，6－8］以及周期性攝動加速度(Once－Per－Revolution Perturbations，OPRP)［9－10］。

1)瞬時(shí)速度脈沖定義在預(yù)先劃分的時(shí)間歷元上，只在有限個(gè)離散的歷元上存在速度的突變，用于改變衛(wèi)星速度狀態(tài)，其他歷元上速度是連續(xù)的。

2)逐段常量經(jīng)驗(yàn)加速度定義在預(yù)先劃分的時(shí)間段上，在任意時(shí)間段內(nèi)經(jīng)驗(yàn)加速度保持常量。方向分別為軌道徑向、切向和法向，切向主要吸收大氣阻力模型誤差，法向吸收未模型化的地球反照光壓誤差、太陽輻射光壓誤差等［1］。

3)周期性攝動加速度模型利用軌道的周期性，假設(shè)未模型化動力學(xué)誤差可用傅立葉級數(shù)的形式描述，每個(gè)軌道周期包括一組傅立葉級數(shù)系數(shù)。

利用瞬時(shí)速度脈沖來對動力學(xué)模型誤差進(jìn)行補(bǔ)償，它的優(yōu)點(diǎn)是模型簡單，無需對含有瞬時(shí)速度脈沖參數(shù)的微分方程進(jìn)行數(shù)值積分，但速度脈沖使衛(wèi)星軌道速度精度下降，且在有限離散歷元上不再連續(xù)。

瞬時(shí)速度脈沖和逐段常量經(jīng)驗(yàn)加速度都稱為偽隨機(jī)參數(shù)模型，模型本身沒有約束性，其值受GPS測量值的影響較大。優(yōu)點(diǎn)是模型可以很好地匹配GPS測量值，當(dāng)測量值的質(zhì)量和精度較高時(shí)定軌精度和穩(wěn)定性較高;缺點(diǎn)是當(dāng)GPS測量值出現(xiàn)間斷或頻繁出現(xiàn)周跳時(shí)，定軌精度則較差(文獻(xiàn)［1］對CHAMP和GRACE星載GPS數(shù)據(jù)存在間斷的時(shí)段定軌表明精度為米級)。

周期性攝動加速度模型采用周期為一個(gè)軌道周期的傅立葉級數(shù)對未模型化動力學(xué)誤差進(jìn)行建模。該模型本身具有約束性，理論上受GPS測量值質(zhì)量和精度的影響相對于偽隨機(jī)參數(shù)模型要小，且待估參數(shù)較少。但有些誤差明顯不具有周期性，無法用該模型來描述，比如姿態(tài)誤差引起的面質(zhì)比變化等［2］。

經(jīng)驗(yàn)加速度模型已在簡化動力學(xué)最小二乘批處理中得到了廣泛的應(yīng)用［1，5－7，10－11］。在批處理過程中，偽隨機(jī)參數(shù)模型的參數(shù)達(dá)到400～800個(gè)，敏感矩陣求解規(guī)模和數(shù)據(jù)存儲量較大，同時(shí)由于GPS觀測量多，法方程的構(gòu)建時(shí)間較長。設(shè)計(jì)矩陣中包含大量敏感矩陣，若都通過數(shù)值積分得到，則計(jì)算量非常大。同時(shí)設(shè)計(jì)矩陣大致呈下三角形式，存在大量零元素，傳統(tǒng)建模方法會使法方程構(gòu)建過程增加大量重復(fù)或無效計(jì)算，且占用大量數(shù)據(jù)存儲空間。

本文首先對三類經(jīng)驗(yàn)加速度模型建模。根據(jù)不同模型的特點(diǎn)，分別給出敏感矩陣求解的有效方法，并定性比較該方法和傳統(tǒng)方法的數(shù)值積分規(guī)模以及存儲量。同時(shí)根據(jù)敏感矩陣求解式的特點(diǎn)，將參數(shù)改正方程組表示為分塊子矩陣的形式，從而得到構(gòu)建法方程的有效方法。然后將該方法擴(kuò)展到基于GPS觀測模型的最小二乘批處理定軌過程中。最后應(yīng)用GRACE－A星載GPS觀測數(shù)據(jù)對三種經(jīng)驗(yàn)加速度模型的法方程構(gòu)建效率進(jìn)行測試與分析。

1 軌道動力學(xué)模型

衛(wèi)星狀態(tài)yT(t)=(rT(t)，vT(t))包括位置矢量r(t)和速度矢量v(t)，動力學(xué)方程表示為

式中，a為地球衛(wèi)星加速度矢量和，p為動力學(xué)參數(shù)。

y(t)關(guān)于初始狀態(tài)y0(t)(y0=y(t0))的偏導(dǎo)數(shù)陣記為狀態(tài)轉(zhuǎn)移矩陣Φ(t，t0)，y(t)關(guān)于p的偏導(dǎo)數(shù)陣為敏感矩陣Sp(t)。它們的一階微分方程分別為［12］:

和

式中:

其中，E為單位矩陣。若設(shè)Φ(t，t0)為齊次微分方程組(2)的解，根據(jù)非齊次微分方程組解結(jié)構(gòu)的一般理論［13］，式(3)的解可表示為

簡化動力學(xué)定軌過程中，同時(shí)對動力學(xué)方程(2)和(3)微分方程和進(jìn)行數(shù)值積分，可以得到任意時(shí)刻的衛(wèi)星狀態(tài)、狀態(tài)轉(zhuǎn)移矩陣和敏感矩陣，用于計(jì)算GPS測量值與模型值的殘差、組成最小二乘參數(shù)估計(jì)法方程。

2 經(jīng)驗(yàn)加速度建模

2．1瞬時(shí)速度脈沖

瞬時(shí)速度脈沖通過瞬時(shí)改變衛(wèi)星速度的大小和方向來補(bǔ)償確定性動力學(xué)模型誤差，作用等同于在加速度上施加一個(gè)單位脈沖函數(shù)。設(shè)定子區(qū)間長度為τ，整個(gè)定軌時(shí)間區(qū)間［t0，tn］劃分為nv段，瞬時(shí)速度脈沖vi=(viR，viT，viN)T作用在離散歷元ti上(ti=t0+i·τ，i=1，…，nv－1)，viR viT和viN對應(yīng)的方向分別為軌道徑向eR(t)、切向eT(t)和法向eN(t)。

動力學(xué)方程需在確定性加速度模型a(t)的基礎(chǔ)上增加瞬時(shí)速度脈沖項(xiàng)［8］:

式中，δ(t)為單位脈沖函數(shù)。vi的敏感矩陣記為Svi(t)，根據(jù)式(6)，其可由下式計(jì)算

當(dāng)t≥ti時(shí)，利用脈沖函數(shù)的性質(zhì)，式(8)右端的積分項(xiàng)為常值矩陣，記為Pi，利用邊值條件得到

此時(shí)Svi(t)表示為Φ(t，t0)和常值矩陣Pi(與i相關(guān))的乘積。最終敏感矩陣的解為

速度脈沖參數(shù)敏感矩陣的微分方程在除有限個(gè)點(diǎn)外是齊次的，敏感矩陣的計(jì)算無需數(shù)值積分。

圖1給出的是與瞬時(shí)速度脈沖參數(shù)相關(guān)的敏感矩陣作用區(qū)域，“灰色”為非零區(qū)域，敏感矩陣可表示為狀態(tài)轉(zhuǎn)移矩陣和常值矩陣的乘積。傳統(tǒng)方法需存儲所有敏感矩陣，存儲規(guī)模為O(nv2)，而對于敏感矩陣的有效計(jì)算方法，只需存儲矩陣Pi(i=1，…，nv－1)即可，存儲規(guī)模為O(nv)，該方法有效降低敏感矩陣的存儲規(guī)模。

圖1 與瞬時(shí)速度脈沖參數(shù)相關(guān)的敏感矩陣作用區(qū)域Fig．1 Activity areas of sensitive matrices related to instantaneous velocity pulses

下面給出瞬時(shí)速度脈沖模型的參數(shù)改正方程和法方程的矩陣分塊建模方法。記zl、hl和εl分別為tl時(shí)刻的測量值、模型值和殘差值，其中tl∈［ti，ti+1)，i=0，1，…，nv－1。假設(shè)動力學(xué)參數(shù)只包括初值狀態(tài)和經(jīng)驗(yàn)加速度參數(shù)(下同)，參數(shù)改正方程為

根據(jù)式(10)，上式進(jìn)一步展開為

其中，bl=zl－h(huán)l，Pm由式計(jì)算。

若考慮［ti，ti+1)內(nèi)所有觀測量hi，參數(shù)改正方程組可表示為子矩陣的形式

其中

若考慮［t0，tn)內(nèi)所有觀測量，參數(shù)改正方程組AvΔxv=bv可表示為分塊子矩陣的形式，設(shè)計(jì)矩陣Av、Δxv和右端項(xiàng)bv分別為

和

最小二乘批處理法方程為

利用矩陣分塊結(jié)構(gòu)，法方程矩陣表示為

和

用矩陣分塊方法構(gòu)建法方程非常有效，且矩陣結(jié)構(gòu)非常清晰:

1)法方程矩陣由Pi(與敏感矩陣相關(guān))以及和(與狀態(tài)轉(zhuǎn)移矩陣相關(guān))組成，i

3)Nmn由左乘，右乘Pn得到;Lm由左乘得到，m，n=0，1，…，n－1。v

2．2逐段常量經(jīng)驗(yàn)加速度

整個(gè)區(qū)間［t0，tn］等間隔劃分為na段，逐段常量經(jīng)驗(yàn)加速度ai=(aiR，aiT，aiN)T在子區(qū)間［t0，tna］上保持常量(ti=t0+iτ，i=0，…，na－1)。ai的方向和瞬時(shí)速度脈沖模型一致。

同理，動力學(xué)方程表示為［8］:

其中

根據(jù)式(6)，敏感矩陣Sai(t)可表示為

當(dāng)t≥ti+1時(shí)，式(22)右端的積分項(xiàng)為常值矩陣，記為Pi+1，其值由式(24)給出。即當(dāng)微分方程為齊次時(shí)，Sai(t)可表示為和常值矩陣Pi+1(與i相關(guān))的乘積。當(dāng)ti≤t＜ti+1時(shí)，則按照微分方程

積分計(jì)算，初值條件Sai(ti)=0。因Sai(t)在整個(gè)區(qū)間上是連續(xù)的，故數(shù)值積分得到的Sai(ti+1)在ti+1時(shí)刻也可表示為Φ(ti+1，t0)與Pi+1的乘積，于是得到常值矩陣

綜上所述，Sai(t)為

逐段常量經(jīng)驗(yàn)加速度的敏感矩陣求解式相比瞬時(shí)速度脈沖式，多了數(shù)值積分項(xiàng)，這是導(dǎo)致敏感矩陣有效區(qū)域，參數(shù)改正方程以及法方程結(jié)構(gòu)與瞬時(shí)速度脈沖參數(shù)模型存在區(qū)別的根本原因。

圖2給出的是與逐段常量經(jīng)驗(yàn)加速度相關(guān)的敏感矩陣有效區(qū)域?！盎疑眳^(qū)域的敏感矩陣可以表示為狀態(tài)轉(zhuǎn)移矩陣和常值矩陣的乘積，無需數(shù)值積分，“×”區(qū)域需進(jìn)行數(shù)值積分，積分和存儲規(guī)模為O(na)。而傳統(tǒng)方法的數(shù)值積分區(qū)域包括所有“×”和灰色區(qū)域，積分和存儲規(guī)模為O(na2)。

圖2 與逐段常量經(jīng)驗(yàn)加速度相關(guān)的敏感矩陣作用區(qū)域Fig．2 Activity areas of sensitive matrices related to piecewise empirical accelerations

逐段常量經(jīng)驗(yàn)加速度的參數(shù)改正方程組也可表示為分塊矩陣的形式

其中，i=0，1，…，na－1，Hi和Pm+1分別見式(14)和(24)。與瞬時(shí)速度脈沖不同，該式多出一項(xiàng):

該項(xiàng)需要數(shù)值積分計(jì)算。

對于所有觀測量，參數(shù)改正方程AaΔxa=ba表示為分塊矩陣的形式，設(shè)計(jì)矩陣Aa和右端項(xiàng)ba分別為

和

法方程為該矩陣表示為兩個(gè)矩陣之和，第一個(gè)矩陣結(jié)構(gòu)與瞬時(shí)速度脈沖模型一致，第二個(gè)矩陣與相關(guān)。設(shè)

其中

子矩陣Nmn和右端項(xiàng)Lm的定義分別與式(18)和(19)相同，m，n=0，1，…，na－1。

矩陣M2a和右端項(xiàng)b2a分別為

和

瞬時(shí)速度脈沖模型的法方程矩陣結(jié)構(gòu)是逐段常量經(jīng)驗(yàn)加速度模型的簡化形式，程序設(shè)計(jì)中可考慮二者間的繼承性。關(guān)于法方程構(gòu)建效率的算例分析見第4節(jié)。

2．3周期性攝動加速度

周期性攝動加速度模型假設(shè)未模型化力具有周期性，可用傅立葉級數(shù)描述，相當(dāng)于對未模型化力做了一種約束，自變量也從時(shí)間(t)變?yōu)樯唤蔷?u(t))。

周期性攝動加速度模型參數(shù)(ciR，siR，ciT，siT，ciN，siN)T在預(yù)先定義的區(qū)間［ti，ti+1)保持常量(ti=t0+iT，i=0，…，ncpr－1)，T為一個(gè)軌道周期，方向與偽隨機(jī)參數(shù)模型一致。動力學(xué)方程為

式中:

周期性攝動加速度模型的敏感矩陣計(jì)算方法、設(shè)計(jì)矩陣和法方程矩陣的分塊結(jié)構(gòu)與逐段常量經(jīng)驗(yàn)加速度完全一致，不再贅述。一般周期性攝動加速度模型參數(shù)個(gè)數(shù)要明顯少于偽隨機(jī)參數(shù)模型，法方程構(gòu)建效率的算例分析見第4節(jié)。

3 簡化動力學(xué)最小二乘批處理

簡化動力學(xué)最小二乘批處理方法(Batch Leastsquares，Batch LSQ)利用GPS非差消電離組合觀測量定軌，測量模型包含接收機(jī)鐘差以及模糊度參數(shù)，除了初始狀態(tài)和經(jīng)驗(yàn)加速度外，動力學(xué)參數(shù)還包括大氣阻力系數(shù)(CD)和太陽輻射光壓系數(shù)(CR)。

最小二乘批處理的法方程結(jié)構(gòu)可見文獻(xiàn)［7］，其中法方程中與動力學(xué)參數(shù)對應(yīng)的矩陣計(jì)算量最大。對比上節(jié)，雖然動力學(xué)參數(shù)多了CD和CR兩項(xiàng)，但矩陣分塊的建模方法仍是有效的。只需形式上將這兩個(gè)參數(shù)組合到式(13)和(26)的子矩陣Hi中［8］，記作

相應(yīng)地，記

那么式(15)和(28)的分塊矩陣中，Hi和Pi將分別被和所替代，法方程矩陣的結(jié)構(gòu)沒有發(fā)生改變。

4 經(jīng)驗(yàn)加速度有效建模方法的計(jì)算效率分析

利用GRACE－A星載GPS載波相位和偽距觀測數(shù)據(jù)，對應(yīng)用在最小二乘參數(shù)估計(jì)過程中的三類經(jīng)驗(yàn)加速度的建模效率進(jìn)行測試。我們只分析矩陣分塊建模方法與傳統(tǒng)建模方法的計(jì)算效率，關(guān)于定軌精度分析由另文給出。

GRACE－A為低軌近圓軌道衛(wèi)星，平均軌道高度460公里。GPS觀測數(shù)據(jù)來自于美國宇航局JPL實(shí)驗(yàn)室物理海洋學(xué)數(shù)據(jù)分發(fā)存檔中心(PO．DAAC)［14］，選取的日期為2006年1月3日;GPS精密星歷以及鐘差來自于歐洲軌道確定中心(CODE)［15］。解算設(shè)置為:批處理時(shí)長為一天(24小時(shí))，數(shù)據(jù)采樣間隔30秒。設(shè)置典型的區(qū)間長度τ:瞬時(shí)速度脈沖和逐段常量經(jīng)驗(yàn)加速度模型為600秒(模型參數(shù)總數(shù)432個(gè))，周期性攝動加速度模型為一個(gè)軌道周期(5662秒，參數(shù)總數(shù)94個(gè))。選取兩組實(shí)驗(yàn)數(shù)據(jù)，第一組只含有偽距觀測量，第二組包含偽距和載波相位觀測量。因?yàn)閮深愑^測量的累加，使得觀測量增加一倍的同時(shí)還引入了模糊度參數(shù)。

測試的筆記本電腦型號為Lenovo T430S，主頻2．9 GHz，內(nèi)存4Gbytes。

4．1偽距觀測量

前文給出了關(guān)于經(jīng)驗(yàn)加速度敏感矩陣的有效計(jì)算方法，已從數(shù)值積分計(jì)算規(guī)模和敏感矩陣存儲規(guī)模的角度上做了定性分析。下面從法方程的構(gòu)建上分析傳統(tǒng)方法和矩陣分塊方法的計(jì)算效率，不包括法方程逆矩陣求解。傳統(tǒng)計(jì)算方法指不對參數(shù)改正方程(式和式)進(jìn)行分塊，直接利用設(shè)計(jì)矩陣相乘得到法方程矩陣。

2006年1月3日的消電離偽距觀測量經(jīng)預(yù)處理后的數(shù)目為21690個(gè)。因最小二乘批處理過程為數(shù)值迭代過程，每次迭代過程的計(jì)算量幾乎相同，為簡明起見只比較一次迭代的計(jì)算時(shí)間。

圖3給出的是三種經(jīng)驗(yàn)加速度模型的法方程構(gòu)建計(jì)算時(shí)間對比，縱軸為CPU計(jì)算時(shí)間。同屬于偽隨機(jī)參數(shù)模型的瞬時(shí)速度脈沖和逐段常量經(jīng)驗(yàn)加速度模型的計(jì)算時(shí)間比較接近，且明顯高于周期性攝動加速度。偽隨機(jī)參數(shù)模型的矩陣分塊計(jì)算方法將計(jì)算效率提高了1倍，而周期性攝動加速度模型幾乎沒有提高。事實(shí)上，偽隨機(jī)參數(shù)模型數(shù)目是周期性攝動加速度的4．6倍，待估參數(shù)越多，利用法方程子矩陣的遞歸關(guān)系減少重復(fù)計(jì)算的效果越明顯，而周期性攝動加速度的參數(shù)少，計(jì)算效率的提升并不明顯。

圖3 三類經(jīng)驗(yàn)加速度模型的法方程構(gòu)建計(jì)算時(shí)間對比，觀測量為GPS偽距消電離組合Fig．3 The computing time comparisons of setting up normal equations(without inversion)for three empirical acceleration models using pseudorange observations

4．2偽距加載波相位觀測量

低軌衛(wèi)星若實(shí)現(xiàn)厘米級精密定軌需求，一般采用偽距加載波相位作為觀測量。預(yù)處理后的觀測量數(shù)目為43380個(gè)，比只考慮偽距觀測量多了一倍，且待估參數(shù)增加了模糊度參數(shù)(428個(gè))，必然加大了法方程的構(gòu)建時(shí)間。

圖4為利用偽距加載波相位消電離組合觀測量的法方程構(gòu)建計(jì)算時(shí)間對比。相比于只利用偽距觀測量(圖3)，無論傳統(tǒng)方法還是矩陣分塊方法，計(jì)算時(shí)間增加約2～3倍。對于偽隨機(jī)脈沖模型，矩陣分塊建模方法將效率提高3．5倍，同樣周期性攝動加速度建模的效率沒有明顯提高。

圖4 三類經(jīng)驗(yàn)加速度模型的法方程構(gòu)建計(jì)算時(shí)間對比，觀測量為偽距加載波相位消電離組合Fig．4 The processing time comparison of setting up normal equation(without inversion)for three empirical acceleration models

5 結(jié)論

本文介紹了地球低軌衛(wèi)星精密定軌的三類常用經(jīng)驗(yàn)加速度模型的有效建模方法，包括瞬時(shí)速度脈沖、逐段常量經(jīng)驗(yàn)加速度以及周期性攝動加速度。因經(jīng)驗(yàn)加速度模型參數(shù)一般較多，傳統(tǒng)建模方法的敏感矩陣數(shù)值積分規(guī)模較大，且構(gòu)建法方程的效率較低。

當(dāng)經(jīng)驗(yàn)加速度敏感矩陣的微分方程是齊次方程時(shí)(除有限歷元或有限區(qū)間外)，敏感矩陣可表示為狀態(tài)轉(zhuǎn)移矩陣和某一常量矩陣的乘積(式(10)和式(25))。瞬時(shí)速度脈沖參數(shù)敏感矩陣無需通過數(shù)值積分計(jì)算，另外兩類模型的敏感矩陣計(jì)算規(guī)模是O(n)(n為參數(shù)個(gè)數(shù))，而傳統(tǒng)方法的計(jì)算規(guī)模是O(n2)。在敏感矩陣的有效求解方法基礎(chǔ)上，用分塊子矩陣構(gòu)建法方程效率會有明顯提高。

利用GRACE－A星載GPS觀測數(shù)據(jù)測試，在典型的精密定軌設(shè)置下，定量分析矩陣分塊建模方法的法方程構(gòu)建效率。觀測量為偽距時(shí)，偽隨機(jī)參數(shù)模型的效率提高1倍;觀測量為偽距加載波時(shí)，效率提高3．5倍。觀測量越多，待估參數(shù)越多，矩陣分塊的建模方法越有效。周期性攝動加速度模型的構(gòu)建效率沒有明顯提高，其參數(shù)數(shù)目明顯小于偽隨機(jī)參數(shù)模型。此外，逐段常量經(jīng)驗(yàn)加速度和周期性攝動加速度模型的法方程是在瞬時(shí)速度脈沖模型的基礎(chǔ)上增加了一個(gè)方程，建模方法具有繼承性。

［1］Kroes R．Preciserelative positioning of formation flying spacecraft using GPS［D］．Delft:Delft University of Technology，2006．

［2］Beutler G，J?ggi A，Hugentobler U，et al．Efficient satellite orbit modelling using pseudo－stochastic parameters［J］．J．Geod．，2006，80(7):353－372．

［3］Beutler G，Brockmann E，Gurtner W，et al．Extended orbit modeling techniques at the CODE processing center of the international GPS service for geodynamics(IGS):theory and initial results［J］．Manuscr Geod，1994，19:367－386．

［4］J?ggi A，Hugentobler U，Beutler G．Pseudo－stochastic orbitmodelingtechniques for low earth orbiters［J］．J．Geod．，2006，80(1):47–60．

［5］ 趙春梅，程鵬飛，益鵬舉．基于偽隨機(jī)脈沖估計(jì)的簡化動力學(xué)衛(wèi)星定軌方法［J］．宇航學(xué)報(bào)，2011，32(4):762－766．［Zhao Chun－mei，Cheng Peng－fei，Yi Peng－ju．Reduced－dynamics satellite orbit determination based on pseudo－stochastic pulse estimation［J］．Journal of Astronautics，2011，32(4):762－766．］

［6］Visser P，Van den Ijssel J．Aiming at a 1cm orbit for low Earth orbiters:reduced－dynamic and kinematic precise orbit determination［J］．Space Science Reviews，2003，108:27－36．

［7］Montenbruck O，Van Helleputte T，Kroes R，et al．Reduced dynamic orbit determination using GPS code and carrier measurements［J］．Aerosp．Sci．Tech．，2005，9(3):261－271．

［8］Swatschina P．Dynamic and reduced－dynamic precise orbit determination of satellites in low earth orbits［D］．Vienna:Vienna University of Technology，2009．

［9］Colombo O L．The dynamics of global positioning orbits and the determination of precise ephemerides［J］．J．Geophys．Res:Solid Earth(1978－2012)．1989，94(B7):9167－9182．

［10］ 朱駿，王家松，陳建榮，等．HY－2衛(wèi)星DORIS厘米級精密定軌［J］．宇航學(xué)報(bào)，2013，34(2):163－169．［Zhu Jun，Wang Jia－song，Chen Jian－rong，et al．Centimeter precise orbit determination for HY－2 via DORIS［J］．Journal of Astronautics，2013，34(2):163－169．］

［11］ 王文彬，劉榮芳．基于雙頻GPS觀測的簡化動力學(xué)最小二乘批處理精密定軌［J］．空間科學(xué)學(xué)報(bào)，2014，34(4):460－467．［Wang Wen－bin，Liu Rong－fang．Precise orbit determination based on reduced dynamic batch LSQ estimation method using dual－frequency GPS observations［J］．Chin．J．Space Sci．，2014，34(4):460－467．］

［12］Montenbruck O，Gill E．Satellite orbits:models，methods and applications［M］．Heidelberg:Springer Verlag，2000．

［13］ 王高雄，周之銘，等．常微分方程(第二版)［M］．北京:高等教育出版社，1982:185－196．

［14］NASA Jet Propulsion Laboratory PO．DAAC．GRACE Level 1B products［OL］．［2013－11－21］．ftp://podaac．jpl．nasa．gov/allData/grace/．

［15］Center for Orbit Determination in Europe(CODE)．High rate clocks and orbits［OL］．［2014－01－08］．ftp://ftp．unibe．ch/aiub/CODE．