王文彬，高 揚

(中國科學(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)驗加速度模型吸收未模型化動力學(xué)誤差，使得簡化動力學(xué)模型(包括確定性動力學(xué)模型和經(jīng)驗加速度模型)精度能夠匹配GPS載波相位的測量精度，同時能夠保持動力學(xué)模型對衛(wèi)星軌道的約束作用，保證軌道解的穩(wěn)定性。

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

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

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

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

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

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

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

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

本文首先對三類經(jīng)驗加速度模型建模。根據(jù)不同模型的特點，分別給出敏感矩陣求解的有效方法，并定性比較該方法和傳統(tǒng)方法的數(shù)值積分規(guī)模以及存儲量。同時根據(jù)敏感矩陣求解式的特點，將參數(shù)改正方程組表示為分塊子矩陣的形式，從而得到構(gòu)建法方程的有效方法。然后將該方法擴(kuò)展到基于GPS觀測模型的最小二乘批處理定軌過程中。最后應(yīng)用GRACE－A星載GPS觀測數(shù)據(jù)對三種經(jīng)驗加速度模型的法方程構(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é)定軌過程中，同時對動力學(xué)方程(2)和(3)微分方程和進(jìn)行數(shù)值積分，可以得到任意時刻的衛(wèi)星狀態(tài)、狀態(tài)轉(zhuǎn)移矩陣和敏感矩陣，用于計算GPS測量值與模型值的殘差、組成最小二乘參數(shù)估計法方程。

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

2．1瞬時速度脈沖

瞬時速度脈沖通過瞬時改變衛(wèi)星速度的大小和方向來補(bǔ)償確定性動力學(xué)模型誤差，作用等同于在加速度上施加一個單位脈沖函數(shù)。設(shè)定子區(qū)間長度為τ，整個定軌時間區(qū)間［t0，tn］劃分為nv段，瞬時速度脈沖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ǔ)上增加瞬時速度脈沖項［8］:

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

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

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

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

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

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

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

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

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

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

其中

若考慮［t0，tn)內(nèi)所有觀測量，參數(shù)改正方程組AvΔxv=bv可表示為分塊子矩陣的形式，設(shè)計矩陣Av、Δxv和右端項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)驗加速度

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

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

其中

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

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

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

綜上所述，Sai(t)為

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

圖2給出的是與逐段常量經(jīng)驗加速度相關(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)驗加速度相關(guān)的敏感矩陣作用區(qū)域Fig．2 Activity areas of sensitive matrices related to piecewise empirical accelerations

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

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

該項需要數(shù)值積分計算。

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

和

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

其中

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

矩陣M2a和右端項b2a分別為

和

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

2．3周期性攝動加速度

周期性攝動加速度模型假設(shè)未模型化力具有周期性，可用傅立葉級數(shù)描述，相當(dāng)于對未模型化力做了一種約束，自變量也從時間(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為一個軌道周期，方向與偽隨機(jī)參數(shù)模型一致。動力學(xué)方程為

式中:

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

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

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

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

相應(yīng)地，記

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

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

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

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

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

4．1偽距觀測量

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

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

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

圖3 三類經(jīng)驗加速度模型的法方程構(gòu)建計算時間對比，觀測量為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)星若實現(xiàn)厘米級精密定軌需求，一般采用偽距加載波相位作為觀測量。預(yù)處理后的觀測量數(shù)目為43380個，比只考慮偽距觀測量多了一倍，且待估參數(shù)增加了模糊度參數(shù)(428個)，必然加大了法方程的構(gòu)建時間。

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

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

5 結(jié)論

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

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

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

［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ī)脈沖估計的簡化動力學(xué)衛(wèi)星定軌方法［J］．宇航學(xué)報，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é)報，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é)報，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．