沐俊山??，劉冰，薛國(guó)虎，戴正旭

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

基于測(cè)量船的實(shí)時(shí)定軌改進(jìn)方法?

沐俊山??，劉冰，薛國(guó)虎，戴正旭

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

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

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

1 引言

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

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

2 定軌算法

2.1 動(dòng)力學(xué)模型

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

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

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

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

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

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

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

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

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

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

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

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

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

圖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 兩種方法定軌半長(zhǎng)軸與理論結(jié)果偏差比較Fig.3 Semi-major axis results generated by the two filters

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

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

3.2 兩種方法的比較

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

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

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

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

［1］陳務(wù)深，陸本魁，馬靜遠(yuǎn).單位矢量法的數(shù)學(xué)模型（MMUVM）極其簡(jiǎn)化形式（PUVM1）的迭代法的收斂性［J］.天文學(xué)報(bào)，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］.北京：國(guó)防工業(yè)出版社，2000. LIU Lin.The Orbit Theory of Spacescraft［M］.Beijing：National Defense Industry Press，2000.（in Chinese）

［4］馬鵬斌，王丹.基于EKF和多信息源融合的空間交會(huì)對(duì)接過(guò)程實(shí)時(shí)定軌方法［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ì)定軌方法［J］.電訊技術(shù)，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—），男，安徽人，碩士，高級(jí)工程師，主要研究方向?yàn)楹教鞙y(cè)控和軌道力學(xué)；

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—），男，江蘇人，碩士，高級(jí)工程師，主要研究方向?yàn)楹教鞙y(cè)控總體；

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.

薛國(guó)虎（1974—），男，江蘇人，碩士，高級(jí)工程師，主要研究方向?yàn)楹教鞙y(cè)控總體；

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年于南京航天航空大學(xué)獲學(xué)士學(xué)位，現(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

國(guó)家高技術(shù)研究發(fā)展計(jì)劃（863計(jì)劃）項(xiàng)目

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

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