馬鵬斌 王 丹

（1中國西安衛(wèi)星測控中心 2宇航動力學(xué)國家重點實驗室）

1 引言

我國天宮一號目標(biāo)飛行器和神舟八號飛船將要進(jìn)行空間交會對接。其測量手段有地（海）基測量站USB測量（測距、測速、測角），中繼星4程測距數(shù)據(jù)，目標(biāo)飛行器和神舟八號飛船間的激光測距和雷達(dá)測距等外測手段。另外在飛船上還搭載有高精度的加速度計?？赏ㄟ^這些數(shù)據(jù)融合來快速確定天宮一號目標(biāo)飛行器和神舟八號飛船的軌道。

2 系統(tǒng)方程

基于J2000慣性坐標(biāo)系，對于航天器的運動方程，采用位置矢量r和速度矢量來描述。運動方程為=f（r，，t），可將多個航天器的軌道同時作為待估量，如目標(biāo)飛行器（r1，）和飛船（r2，）兩個航天器可同時求解。另外還可增加中繼星作為待估量。

在軌道確定中需要求解的量往往不限于航天器狀態(tài)r、，如求解大氣阻尼系數(shù)CD，測量系統(tǒng)誤差等。設(shè)其它需要求解的模型參數(shù)為c，因c為常數(shù)，所以有=0。

令

航天器運動狀態(tài)方程可擴(kuò)展成如下形式：

代表了一個由n個非線性一階常微分方程組成的系統(tǒng)[1]。

3 動力學(xué)模型

航天器運動方程基于J2000慣性坐標(biāo)系，位置矢量r和速度矢量描述。運動方程即為初值問題[1]

右函數(shù)f可分為f0和fε兩部分，f0為二體問題下航天器的加速度，fε為其它各種攝動力加速度之和，可表示如下：

針對序貫處理的需求，對運動方程采用了單步法RKF7（8）方法求解。天宮一號和神舟八號可同時積分求解。

一般情況下，各種攝動力的量級如下表所示[2]：

表1 攝動力量級表

兼顧計算效率，只考慮米級以上的攝動影響，包括地球引力場、日月引力、太陽光壓、大氣阻尼攝動。由于神舟八號飛船搭載有高精度的加速度計，可以只計算地球引力場和日月引力，交會對接過程中的軌道機動推力和其他非引力攝動力可直接使用加速度計實測數(shù)據(jù)。

系統(tǒng)動力學(xué)模型的狀態(tài)轉(zhuǎn)移矩陣為：

在計算傳統(tǒng)動力學(xué)模型的狀態(tài)轉(zhuǎn)移矩陣時，考慮地心引力和J2項攝動即可完全滿足精度要求。考慮地心引力和J2項攝動，則有：

4 測量模型

建立量測方程，將系統(tǒng)在某時刻的測量量和系統(tǒng)狀態(tài)量聯(lián)系起來。狀態(tài)量為J2000慣性坐標(biāo)系下的位置矢量、速度矢量 r，。測量數(shù)據(jù)有：陸（海）基測量站USB測量（測距、測速、測角），目標(biāo)飛行器和神舟八號飛船間的激光測距和雷達(dá)測距，中繼星4程距離和等。

4.1 陸（海）基測量站USB測量

陸（海）基測量站USB測量在測站地平坐標(biāo)系中的觀測量為：斜距R、斜距變化率、方位角A、俯仰角E。設(shè)測站地平坐標(biāo)系中的直角坐標(biāo)為（x，y，z，，，），則：

觀測量對狀態(tài)量的偏導(dǎo)數(shù)為：

上述各式中出現(xiàn)的B和L分別為觀測站的大地緯度和經(jīng)度，R測站為J2000慣性坐標(biāo)系下測站地心向徑矢量。

4.2 目標(biāo)飛行器和飛船間距離測量

令航天器 1 的位置矢量為r1=（x1，y1，z1），航天器 2的位置矢量為r2=（x2，y2，z2），則航天器 1 對航天器2的測距為：

ρ對航天器1，2狀態(tài)矢量的偏導(dǎo)分別為：

4.3 中繼星四程距離和測量

已知中繼星精密星歷，可將四程距離和處理為中繼星對所求解航天器的測距數(shù)據(jù)，認(rèn)為中繼星是一個地面站，其測量模型同USB測距模型相同。

5 序貫估計方法

5.1 擴(kuò)展卡爾曼濾波基本原理[3]

設(shè)狀態(tài)方程為（t）=f[x（t），t]+w（t），w為模型噪聲；量測方程為y（t）=h[x（t），t]+R（t），R為觀測噪聲，線性化后可得xk=Φ（tk，tj）xj+wj，yk=Hkxk+Rk。

已知估值及其協(xié)方差矩陣Pk，對tk+1時刻的一個新觀測數(shù)據(jù),卡爾曼濾波基本方法如下：

計算tk+1時刻的預(yù)測軌道；

計算tk+1時刻的預(yù)測協(xié)方差矩陣

進(jìn)行狀態(tài)更新，得到tk+1時刻的估計值

計算tk+1時刻的誤差協(xié)方差矩陣Pk+1=

對運動方程重新初始化，把運動方程中原參考軌道換成最新得出的tk時刻的最優(yōu)估值。把運動方程積分至tk+1時刻以處理tk+1時刻的觀測數(shù)據(jù)。

在卡爾曼濾波實際應(yīng)用中很可能出現(xiàn)的一個現(xiàn)象是處理了一定量的數(shù)據(jù)后濾波值出現(xiàn)發(fā)散現(xiàn)象。每次在協(xié)方差矩陣上加一個小的噪音項Q（即過程噪聲）可以抑制矩陣變得越來越小[4]。則有：

5.2 野值判斷

野值剔除是序貫定軌中很重要的環(huán)節(jié)，一個大的野值可能導(dǎo)致序貫定軌的發(fā)散。

設(shè)得到t時刻的測量量為O，其測量方差為R。預(yù)報至t時刻的軌道狀態(tài)矢量為X，先驗方差矩陣為P，C為測量量的預(yù)報值。

先驗方差矩陣P的對角線元素可以近似認(rèn)為是軌道精度，則

可近似認(rèn)為：

測距預(yù)報值精度SR為Rp

測角預(yù)報值精度SAE為arctan（Rp/CR）

測速預(yù)報值精度SD為（P11P44+P22P55+P33P66）/CR

O-C為觀測值減計算值。若滿足：

O-C＞n1·R 且 O-C＞n2·S

則測量值為野值。其中，R為測量方差，S為測量值預(yù)報精度。n1，n2為人為設(shè)定門限，若使用5σ剔除，則n1，n2為5。如果某一時刻的測量量只有某一項為野值，則剔除該量，其他量仍可使用。

6 實時軟件流程

實時軟件流程如下：

第一步：計算起步初值：

第二步：得到新的tk+1時刻觀測數(shù)據(jù)后，進(jìn)行軌道外推，計算tk+1時刻的預(yù)測軌道；

第三步：計算tk到tk+1時刻的狀態(tài)轉(zhuǎn)移矩陣Φ（tk+1，tk）；

第四步：計算tk+1時刻的預(yù)測協(xié)方差矩陣=Φ（tk+1，tk）PkΦT（tk+1，tk）；

第五步：計算觀測量計算值；

第六步：計算觀測偏導(dǎo)數(shù)；

第七步：進(jìn)行數(shù)據(jù)質(zhì)量檢驗，判斷O-C是否滿足規(guī)定門限要求，或O-C與預(yù)測O-C相差滿足規(guī)定值，如不滿足則剔除此值，返回第二步；

第九步：進(jìn)行狀態(tài)更新，得到tk+1時刻的估計值

第十步：計算tk+1時刻的誤差協(xié)方差矩陣Pk+1=

第十一步：比較Pk-1、Pk與Pk+1，判斷是否發(fā)散，如連續(xù)發(fā)散則返回第一步；

第十二步：輸出航天器軌道狀態(tài)量，也可預(yù)報下一時刻引導(dǎo)值并輸出；

第十三步：返回第二步，如數(shù)據(jù)結(jié)束則轉(zhuǎn)到下一步；

第十四步：進(jìn)行結(jié)束處理。

7 結(jié)論

本文針對我國空間交會對接任務(wù)的快速軌道計算，融合使用USB測量數(shù)據(jù)、中繼星測量數(shù)據(jù)、目標(biāo)飛行器和飛船間測量數(shù)據(jù)以及船載高精度加速度計數(shù)據(jù)等，利用擴(kuò)展卡爾曼濾波方法來快速確定目標(biāo)飛行器和飛船的軌道?？捎糜谖覈臻g交會對接任務(wù)的實時軌道監(jiān)視。 ◇

［1］李濟(jì)生.人造衛(wèi)星精密軌道確定.北京：解放軍出版社,1995

［2］劉林.人造地球衛(wèi)星軌道力學(xué).北京:高等教育出版社,1992

［3］賈沛漳，朱征桃.最優(yōu)估計理論.科學(xué)出版社.1984.

［4］A.C.Long et al.Goddard Trajectory Determination System.GSFC.1989