基于測量船的實時定軌改進方法?

沐俊山??，劉冰，薛國虎，戴正旭

（中國衛(wèi)星海上測控部，江蘇江陰214431）

基于測量船的實時定軌改進方法?

沐俊山??，劉冰，薛國虎，戴正旭

（中國衛(wèi)星海上測控部，江蘇江陰214431）

為了適應多目標任務的測定軌需求，進一步提高測量船定軌的時效性，對時間序列最優(yōu)定軌方法進行了研究，分析了定軌工作中涉及的飛行動力學模型，對序貫定軌的過程進行了貝葉斯統(tǒng)計描述?；诠こ虘玫目尚行裕捎脭U展卡爾曼濾波方法和無跡濾波方法，并對兩者進行了改進，在計算上采用平方根擴展卡爾曼濾波方法和平方根無跡濾波方法，在定軌模式上，采用一種新的混合定軌方法。改進后的方法在保持精度的基礎上提高了濾波器時間更新效率，模擬數(shù)值計算和實戰(zhàn)數(shù)據(jù)驗證的結(jié)果表明，兩種濾波定軌方法均能夠在較短的時間內(nèi)收斂，并達到預期的定軌精度。

測量船；最優(yōu)軌道確定；平方根擴展卡爾曼濾波；平方根無跡濾波

1 引言

航天測量船主要從事對航天器的跟蹤測量、遙測及遙控工作。其中，對目標航天器的彈道狀態(tài)估計和軌道測量與確定是測量船工作的重要內(nèi)容之一。當前，測量船統(tǒng)計定軌方法中，無論是事后處理還是短弧段的準實時處理，均采用基于最小二乘原理的批處理方法。單位矢量法［1］和微分改正方法均屬于批處理方法，此類方法的特點是測量數(shù)據(jù)需要進行離線處理。

隨著衛(wèi)星編隊任務等應用需求，多目標發(fā)射任務越來越多。在這種發(fā)射任務中，火箭不僅攜帶多個有效載荷，而且還包括多個衛(wèi)星之間的承載艙，在各目標分離時，多個目標之間存在一定的碰撞風險，因此需要對各目標的位置、速度狀態(tài)矢量盡可能快速準確地進行估計和預報，并快速地評估和預報碰撞風險。目前的事后定軌方法不僅在時效上滯后，另外無法加入有動力的數(shù)據(jù)，導致有效定軌弧段較短。針對新的任務需求，我們研制并改進了基于時間序列最優(yōu)估計的實時定軌系統(tǒng)。本文就研制過程中涉及的相關技術進行闡述。

2 定軌算法

2.1 動力學模型

飛行器在飛行過程中受到多種力的作用，其中包括地球二體和非球形引力、大氣阻力、太陽光壓、日月和行星的攝動、潮汐等作用。對近地軌道衛(wèi)星所受到的力除了地球引力和大氣阻力外，其他力的作用產(chǎn)生的加速度均低于10－7km2/s［2］。當飛行器在某個時刻的狀態(tài)矢量確定后，其在這個時刻的瞬時密切軌道根數(shù)也就確定了。飛行器狀態(tài)矢量包括其速度r和位置v。該飛行器在其他時刻的狀態(tài)可以使用當前位置、速度基于已知的動力學模型進行外推，目前的軌道計算方法既可以使用分析方法［3］也可以使用數(shù)值方法進行外推，通常這種計算基于某個地心慣性坐標系進行，例如J2000.0。

其中，r和v是飛行器在地心慣性系下的位置、速度矢量；式（2）右邊第一項為地球和飛行器之間理想的二體模型，第二項包括各種攝動作用；anonsphere為地球非球形引力作用，asun/moon為日月攝動，adrag為大氣阻力攝動，asunlight為太陽光壓攝動，atide為固體潮作用；anonsphere是非球形引力勢Unonshpere的梯度，Plm為締合勒讓德多項式形成的系數(shù)，Clm、Slm為引力場系數(shù)。在基于高速率測量數(shù)據(jù)的實時定軌應用中，上述攝動項可以根據(jù)測量精度進行裁剪，軌道外推和誤差狀態(tài)轉(zhuǎn)移矩陣的攝動影響一般分別進行考慮。例如，利用分析模型計算誤差狀態(tài)轉(zhuǎn)移時可以只考慮非球形引力J2項作用［4］，但是在進行軌道外推時，則需要根據(jù)測量數(shù)據(jù)精度來匹配合適的外推模型精度，一般情況下，大氣攝動也考慮在內(nèi)。

其中，Cd為大氣阻尼系數(shù)，A為飛行器垂直與速度矢量方向上的截面面積，ms為飛行器質(zhì)量，ρ是飛行器所在位置的大氣密度，υrel為飛行器相對大氣的飛行速度。

2.2 最優(yōu)實時軌道確定的系統(tǒng)模型

軌道確定中的動力學模型和測量模型均為非線性的模型。其中力學模型為連續(xù)時間模型，可以轉(zhuǎn)化為適合計算機處理的離散模型［5］。而測量模型是自然的離散時間模型。離散的動力學模型和測量模型可以用下式表達：

測量模型和動力學模型分別在上文中描述，其中f為表達動力學模型式（1）和式（2）的過程模型，h為測量模型［6］。xk∈Rn為在離散時刻k的n×1維狀態(tài)矢量，zk∈Rm為該時刻m×1維測量元素矢量。wk∈Rn是n×1維過程噪聲矢量，υ∈Rm為m ×1維測量噪聲矢量。υ和w均假設為期望為0的高斯隨機過程，并且

其中，δ為Kroneck函數(shù)。上述過程和測量模型系統(tǒng)模型可以定義為基于馬爾可夫鏈的貝葉斯統(tǒng)計模型：

在項目研制過程中，我們分別實現(xiàn)了多種非線性濾波方法，其中包括擴展卡爾曼濾波（EKF）［7－8］方法、無跡卡爾曼濾波（UKF）、高斯赫米特濾波（GHF）和粒子濾波（PF）。鑒于工程應用的可行性，我們選擇了計算負荷相對較低的EKF和UKF并對其進行數(shù)值計算的穩(wěn)定性改進。計算字長限制導致計算舍入誤差的積累將導致狀態(tài)協(xié)方差的估計值失去非負定性，負定的狀態(tài)協(xié)方差將使卡爾曼增益的計算失真，從而使從殘差中提取的反饋補償信息失真，最終導致濾波器發(fā)散。解決計算舍入誤差問題的一個主要方式是通過平方根方法，即在濾波的過程中不計算狀態(tài)協(xié)方差，而計算狀態(tài)協(xié)方差的平方根。使用平方根方法不僅能保證狀態(tài)協(xié)方差的非負定性，并且在數(shù)值計算中，使用一半的計算機字長即可達到狀態(tài)協(xié)方差的計算精度。平方根方法是一大類方法，項目主要采用QR和Cholesky兩種矩陣分解技術進行平方根計算，形成平方根擴展卡爾曼濾波（SR-EKF）和平方根無跡濾波方法（SR-UKF）。限于篇幅，這里略去具體的濾波器數(shù)學模型。

3 數(shù)值實例與分析

3.1 數(shù)值計算結(jié)果

在實時定軌系統(tǒng)的實現(xiàn)過程中，還要涉及許多問題，例如：時間系統(tǒng)和坐標系統(tǒng)的實現(xiàn)與轉(zhuǎn)換、測量數(shù)據(jù)在線野值剔除、飛行器外推動力學模型、EKF定軌使用的狀態(tài)轉(zhuǎn)換矩陣的計算方法、濾波器參數(shù)定制等內(nèi)容。例如，在外推動力學模型工作中，我們分別實現(xiàn)了精密的數(shù)值法外推模型和基于平均根數(shù)理論的分析法外推模型等，在數(shù)值方法中，可以根據(jù)測量數(shù)據(jù)精度，人工選擇各種攝動模型及階數(shù)。在狀態(tài)轉(zhuǎn)換矩陣的形成上分別采用了基于受攝力偏微分的分析法、基于單位矢量的分析法和基于有限差分的數(shù)值法。我們分別使用了模擬和實測數(shù)據(jù)對系統(tǒng)進行了測試。

圖1～3為利用一組8min的模擬測量數(shù)據(jù)進行定軌的結(jié)果，定軌采用秒節(jié)點雷達測量數(shù)據(jù)，模擬數(shù)據(jù)中測距隨機誤差標準差設置為100 m，測角隨機誤差標準差設置為0.02°。兩個濾波器的軌道外推模型采用數(shù)值法外推模型，攝動包括2×0階非球型引力和大氣阻力攝動。在定軌模型中，我們采用了一種混合定軌方案：在SR-EKF濾波器中，使用數(shù)值法進行軌道外推，使用分析法并進行有限差分形成誤差狀態(tài)轉(zhuǎn)移矩陣；在SR-UKF濾波器中，在sigma點形成的誤差橢圓中心點利用數(shù)值法外推，而周圍的點使用分析法外推。這種混合定軌方法保持了數(shù)值法定軌的精度，充分利用了分析外推的速度，使濾波器在保證足夠精度的前提下提高時間相應的效率。

圖1為利用SR-EKF定軌結(jié)果與參考軌道在參考軌道RIC（徑向、切向、法向）坐標系下位置分量誤差，圖2為利用SR-EKF定軌結(jié)果與參考軌道在參考軌道RIC坐標系下速度分量誤差，圖3為利用SREKF和SR-UKF進行定軌的結(jié)果與參考軌道在半長軸上的對比。

圖1 SR－EKF定軌結(jié)果RIC位置分量誤差Fig.1 The relative errors in RIC position components under SR－EKF

圖2 SR-EKF定軌結(jié)果RIC速度分量誤差Fig.2 The relative errors in RIC velocity components under SR-EKF

圖3 兩種方法定軌半長軸與理論結(jié)果偏差比較Fig.3 Semi-major axis results generated by the two filters

從各圖可以看出，兩個濾波器從啟動開始計算到60 s，收斂速度較快，后續(xù)收斂速度減慢，其中在RIC坐標系位置分量收斂緩慢，在切向分量上存在一定波動，這是由于飛行器在切向分量上速度最大，估計誤差也最大。在RIC速度分量上，100 s后收斂緩慢而平滑，由于初軌確定工作中最重要的參數(shù)就是半長軸，其直接決定軌道的能量，而速度的準確估計很大程度決定半長軸的精度，與參考軌道相比，半長軸上在60 s后已經(jīng)達到2 km以內(nèi)，在100 s以后，定軌結(jié)果已經(jīng)收斂并保持200m以內(nèi)。

另外，我們還使用歷史任務測量數(shù)據(jù)對兩個濾波方法進行了進一步驗證，并將計算結(jié)果與事后定軌方法結(jié)果進行了比較。通過驗證，兩個改進的實時定軌均能夠處理隨機性不穩(wěn)定的實戰(zhàn)任務測量數(shù)據(jù)而未出現(xiàn)不收斂的情況。另外，實時定軌結(jié)果與單位矢量法的定軌結(jié)果精度上基本相當。

3.2 兩種方法的比較

從兩種方法互相比較的結(jié)果看，對于測量船的測量數(shù)據(jù)精度，兩者定軌結(jié)果差別非常小。由圖3可知，UKF在前100 s比SR-EKF對軌道初值和協(xié)方差的不確定性更為敏感，比SR-EKF偏差較大，在150 s后SR-UKF濾波估計結(jié)果略優(yōu)于EKF，差別在20m以內(nèi)。

我們利用SR-EKF和SR-UKF從功能上實現(xiàn)了對飛行器的短弧段、高密集數(shù)據(jù)的實時軌道確定，時間響應性能和精度均達到了預期的水平。兩種定軌方法存在一定區(qū)別。SR-EKF需要計算誤差狀態(tài)轉(zhuǎn)移矩陣，而SR-UKF無需這一環(huán)節(jié)。完整的誤差狀態(tài)轉(zhuǎn)移矩陣需要分析各攝動力學的雅克比矩陣，對于二階SR-EKF甚至要分析得出海森矩陣，并對其進行積分，不僅分析過程復雜，計算消耗也大。具體應用中，根據(jù)精度需求往往利用數(shù)值法或者分析法進行簡化和裁剪。SR-UKF無需計算誤差狀態(tài)轉(zhuǎn)移矩陣，直接通過sigma點外推計算并加權而得。但是，對于航天器狀態(tài)的六維矢量，需要積分器對每次測量進行13次的軌道和測量外推，對于高精度的數(shù)值法軌道外推，也會造成一定的計算負擔，因此我們在兩種濾波器中均采用了新穎的混合定軌方案，既保證了精度又提高了計算速度。

4 結(jié)束語

為了適應新的任務需求，我們研究并改進了基于測量船的時間序列最優(yōu)定軌方法。模擬和實戰(zhàn)數(shù)據(jù)的測試驗證表明，系統(tǒng)在測量弧段較短的情況下實現(xiàn)了預期的收斂性和精度。后續(xù)需要進一步研究的內(nèi)容包括如何提高原始測量數(shù)據(jù)的精度和進一步融入飛行器有動力段的測量數(shù)據(jù)以增加有效定軌弧段。

［1］陳務深，陸本魁，馬靜遠.單位矢量法的數(shù)學模型（MMUVM）極其簡化形式（PUVM1）的迭代法的收斂性［J］.天文學報，2007，48（3）：343－354. CHENWu-sheng，LU Ben-kui，MA Jing-yuan.The convergence prerformance ofMMUVM and PUVM1［J］.Astronomy Journal，2007，48（3）：343－354.（in Chinese）

［2］Vallado D A.Fundamentals ofastrodynamics and applications［M］.3rd ed.New York：Springer，2007.

［3］劉林.航天器軌道理論［M］.北京：國防工業(yè)出版社，2000. LIU Lin.The Orbit Theory of Spacescraft［M］.Beijing：National Defense Industry Press，2000.（in Chinese）

［4］馬鵬斌，王丹.基于EKF和多信息源融合的空間交會對接過程實時定軌方法［J］.載人航天，2001（3）：27－30. MA Peng-bing，WANG Dan.Real-time Orbit Determination for Space Rendezvous and Docking Based on EKF and Multisource Information Fusion［M］.Manned Spaceflight，2001（3）：27－30.（in Chinese）

［5］Jazwinski A H.Stochastic processes and filtering theory［M］. New York：Academic Press，1970.

［6］李紅艷，沐俊山，傅敏輝，等.星箭分離前后的聯(lián)合統(tǒng)計定軌方法［J］.電訊技術，2011，51（12）：48－52. LIHong-yan，MU Jun-shan，F(xiàn)UMin-hui，etal.A New United Statistical Orbit Determination Method for TT＆C Ships［J］. Telecommunication Engineering，2011，51（12）：48－52.（in Chinese）

［7］Montenbruck O，Eberhard G.Satellite orbits：models，methods，and applications［M］.Berlin：Springer，2000.

［8］Tapley BD.Statisticalorbitdetermination［M］.London：Elsevier Academic Press，2004.

沐俊山（1973—），男，安徽人，碩士，高級工程師，主要研究方向為航天測控和軌道力學；

MU Jun-shan was born in Anhui Province，in 1973.He is now a senior engineerwith theM.S.degree.His research interests include maritime space TT＆C technology and satellite orbitmechanics.

Email：mujunshan＠gmail.com

劉冰（1968—），男，江蘇人，碩士，高級工程師，主要研究方向為航天測控總體；

LIU Bing was born in Jiangsu Province，in 1968.He is now a senior engineer with the M.S.degree.His research concerns space TT＆C system.

薛國虎（1974—），男，江蘇人，碩士，高級工程師，主要研究方向為航天測控總體；

XUEGuo-huwas born in Jiangsu Province，in 1974.He isnow a senior engineer with the M.S.degree.His research concerns maritime space TT＆C system.

戴正旭（1983—），男，江蘇泰興人，2006年于南京航天航空大學獲學士學位，現(xiàn)為助理工程師，主要從事彈道數(shù)據(jù)處理、精度分析等方面的研究工作。

DAI Zheng-xu was born in Taixing，Jiangsu Province，in 1983.He received the B.S.degree from Nanjing University of Aeronautics and Astronautics in 2006.He is now an assistant engineer.His research concerns trajectory data processing and precision analysis.

Improvement of Optimal Sequential Orbit Determ ination M ethod for TT＆C Ship

MU Jun-shan，LIU Bing，XUEGuo-hu，DAIZheng-xu
（China Satellite Maritime Tracking and Controlling Department，Jiangyin 214431，China）

To satisfy the requirementofmulti-target launchmission and improve the time performance oforbitdetermination of spacecraft for TT＆C ship，new sequential optimal estimation methods are proposed and studied. The dynamic processmodel in spacecraft orbit determination is analyzed，and a Bayesian description for the Markov process is described.For the engineering application of theoreticalmodels，improved methods SR-EKF and SR-ULF based on two typical nonlinear filters are proposed.Furthermore，a new hybrid orbit determination model is adopted which utilizes both high precision numerical propagator and high speed analytical propagator. The simulation and historicalmeasurement data verification shows that both filters have achieved the expected accuracy，speed and convergence.

TT＆C ship；real-time orbit determination；square root extended Kalman filter；square root unscented Kalman filter

V556

A

1001－893X（2013）02－0177－05

10.3969/j.issn.1001－893x.2013.02.013

2012－08－01；

2012－10－17 Received date：2012－08－01；Revised date：2012－10－17

國家高技術研究發(fā)展計劃（863計劃）項目

Foundation Item：The National High－tech R＆D Program（863 Program）of China

??通訊作者：mujunshan＠gmail.com Corresponding author：mujunshan＠gmail.com