基于線性協(xié)方差分析的接近天體在軌導(dǎo)航方法研究*

徐 超，黃翔宇，王大軼

(1.北京控制工程研究所，北京 100190；2.空間智能控制技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100190)

對(duì)于接近深空天體的探測(cè)任務(wù)，探測(cè)器的飛行距離遠(yuǎn)、運(yùn)行時(shí)間長(zhǎng)、環(huán)境未知性較大，其導(dǎo)航和控制系統(tǒng)在實(shí)時(shí)性、精度和可靠性方面的要求更高.傳統(tǒng)的基于地面測(cè)控網(wǎng)的導(dǎo)航與控制方法受到精度、傳輸時(shí)間、運(yùn)行成本和可觀測(cè)弧段等方面的限制，往往難以滿足任務(wù)需要[1].

目前采用的接近段導(dǎo)航方式主要有兩類：一類是結(jié)合地面無(wú)線電測(cè)量與光學(xué)成像測(cè)量的地面組合導(dǎo)航方式[2-3]；另一類是單純依靠導(dǎo)航相機(jī)進(jìn)行的光學(xué)自主導(dǎo)航.光學(xué)自主導(dǎo)航完全依賴于探測(cè)器上光學(xué)敏感器提供的觀測(cè)信息，由于受光學(xué)敏感器精度、視場(chǎng)和目標(biāo)天體光照等條件的限制，很難滿足探測(cè)器對(duì)精度和可靠性的要求.融合光學(xué)信息和無(wú)線電信息的地面組合導(dǎo)航方式雖然滿足探測(cè)器對(duì)導(dǎo)航精度和可靠性的要求，但其需要將光學(xué)測(cè)量信息傳回地面站，再將地面站計(jì)算得到的導(dǎo)航信息返回探測(cè)器，這樣無(wú)法滿足實(shí)時(shí)性的要求.因此，為充分發(fā)揮光學(xué)和無(wú)線電導(dǎo)航的優(yōu)勢(shì)，有必要研究在軌融合光學(xué)信息和無(wú)線電信息的自主導(dǎo)航方法，從而提高在軌導(dǎo)航的精度和可靠性[4].

由于光學(xué)成像測(cè)量可以獲得多種觀測(cè)信息[5]，為選擇合適的光學(xué)觀測(cè)信息和無(wú)線電信息在軌融合，本文以火星接近段為背景結(jié)合線性協(xié)方差分析方法和基于B平面參數(shù)的誤差分析方法，分析了不同光學(xué)觀測(cè)信息和無(wú)線電測(cè)量信息在軌融合的自主導(dǎo)航方案.首先，介紹了標(biāo)準(zhǔn)的線性協(xié)方差分析方法；然后給出了接近段基于光學(xué)和無(wú)線電的導(dǎo)航測(cè)量模型，并給出了基于這兩種測(cè)量的自主導(dǎo)航方案；其次，引入了B平面點(diǎn)位誤差橢圓描述導(dǎo)航誤差的置信度；最后，利用基于線性協(xié)方差分析的B平面參數(shù)誤差分析方法對(duì)給出的自主導(dǎo)航方案進(jìn)行了仿真分析.

1 線性協(xié)方差分析方法

線性協(xié)方差分析方法是由Geller和Stastny提出的一種新的導(dǎo)航控制誤差分析工具，該方法只需要運(yùn)行一次即可獲得與蒙特卡羅仿真相似的結(jié)果從而大大節(jié)省仿真時(shí)間[6-7].

利用蒙特卡羅仿真進(jìn)行導(dǎo)航誤差分析是通過(guò)對(duì)N次仿真試驗(yàn)結(jié)果進(jìn)行統(tǒng)計(jì)分析，從而得到導(dǎo)航誤差的統(tǒng)計(jì)信息.真實(shí)狀態(tài)偏差協(xié)方差陣Dtrue、導(dǎo)航狀態(tài)偏差協(xié)方差陣Dnav和導(dǎo)航誤差協(xié)方差Ptrue的計(jì)算公式分別為

(1)

假設(shè)導(dǎo)航系統(tǒng)真實(shí)狀態(tài)模型為

(2)

其中，w(t)為零均值的高斯白噪聲序列，噪聲方差為E(w(t)wT(t′))=Q(t)δ(t-t′)，Q(t)為噪聲方差陣，δ(t-t′)為Dirac-δ函數(shù).導(dǎo)航敏感器在tk時(shí)刻得到觀測(cè)量為

zk=h(xk,tk)+vk

(3)

圖1 狀態(tài)描述

(4)

(5)

(6)

(7)

(8)

(9)

(10)

由于系統(tǒng)真實(shí)狀態(tài)偏差δx(t)不受測(cè)量更新影響，故有

(11)

(12)

(13)

取擴(kuò)展?fàn)顟B(tài)協(xié)方差為Θ=E(XXT)，則有

(14)

(15)

則式中各協(xié)方差陣可按下式計(jì)算:

Dtrue=E(δx(t)(δx(t))T)

=[In×n0n×n]Θ[In×n0n×n]T

=[0n×nIn×n]Θ[0n×nIn×n]T

=[In×n-In×n]Θ[In×n-In×n]T

綜上，式(8)、(14)、(15)即組成線性協(xié)方差分析方法的協(xié)方差時(shí)間和測(cè)量更新方程.可看出利用線性協(xié)方差分析方法只需運(yùn)行一次即可得到與蒙特卡羅仿真相似的導(dǎo)航性能分析結(jié)果.

對(duì)于導(dǎo)航問(wèn)題，因?yàn)樵诜抡嬷刑綔y(cè)器實(shí)際軌道作為已知量，因此可認(rèn)為實(shí)際軌道和標(biāo)稱軌道相同，則上述求導(dǎo)航誤差協(xié)方差陣的過(guò)程[8]可簡(jiǎn)化為

其中，

(17)

2 接近段導(dǎo)航測(cè)量

接近目標(biāo)天體段，導(dǎo)航系統(tǒng)主要以目標(biāo)天體為觀測(cè)對(duì)象，通過(guò)導(dǎo)航相機(jī)獲取包含背景恒星的目標(biāo)天體光學(xué)圖像，經(jīng)過(guò)圖像處理后得到天體視半徑、天體中心點(diǎn)以及恒星點(diǎn)像素信息.同時(shí)通過(guò)在探測(cè)器上安裝的接收機(jī)接收地面站發(fā)射的無(wú)線電導(dǎo)航信息，在軌獲取探測(cè)器與地面站間的距離和徑向速度等觀測(cè)量.

2.1 基于光學(xué)信息的測(cè)量

2.1.1 目標(biāo)天體中心像素信息測(cè)量

由光學(xué)導(dǎo)航原理可得，目標(biāo)天體中心點(diǎn)像素的觀測(cè)模型為

(18)

其中Cij(i=1,2,3;j=1,2,3)為目標(biāo)天體J2000慣性坐標(biāo)系到相機(jī)坐標(biāo)系的方向余弦矩陣CcI中對(duì)應(yīng)的元素，由姿態(tài)敏感器給出；[x,y,z]T為探測(cè)器在目標(biāo)天體J2000慣性坐標(biāo)系中位置；f為導(dǎo)航相機(jī)焦距；kx、ky為由毫米轉(zhuǎn)換為像素的參數(shù)；p0、l0為導(dǎo)航敏感器中心像元和像線；vp、vl為觀測(cè)噪聲.

2.1.2 星光角距測(cè)量

通過(guò)目標(biāo)天體中心點(diǎn)像素和多個(gè)恒星點(diǎn)像素信息可獲得多個(gè)星光角距信息.星光角距測(cè)量是對(duì)目標(biāo)天體與恒星視線夾角的測(cè)量，其觀測(cè)模型為

ψsk=arccos

(19)

其中，(x,y,z)為探測(cè)器在目標(biāo)天體J2000慣性系下位置，βsk,σsk為識(shí)別的第k顆恒星的赤經(jīng)和赤緯，k=1,2,…,n，n為識(shí)別的恒星數(shù)，本文相關(guān)方案中取n=4.

2.1.3 目標(biāo)天體視半徑測(cè)量

通過(guò)對(duì)導(dǎo)航相機(jī)拍攝的目標(biāo)天體圖像進(jìn)行處理可得到目標(biāo)天體的視半徑信息，其觀測(cè)模型為

(20)

其中，R為目標(biāo)天體參考半徑，r為探測(cè)器相對(duì)目標(biāo)天體中心的距離.

2.2 基于無(wú)線電信息的測(cè)量

無(wú)線電導(dǎo)航信息通常包含多普勒信息和距離測(cè)量信息，多普勒觀測(cè)得到的是探測(cè)器與地面站之間的相對(duì)速度vr，距離測(cè)量得到的是探測(cè)器與地面站之間的距離rr.通過(guò)在探測(cè)器上安裝接收機(jī)，在一定條件下探測(cè)器可接收到地面站發(fā)送的無(wú)線電信息，因此，可在軌獲取探測(cè)器與地面站間的距離和徑向速度信息，其觀測(cè)模型為

(21)

其中

(22)

3 接近段自主導(dǎo)航方案

接近段自主導(dǎo)航系統(tǒng)以探測(cè)器在目標(biāo)天體慣性系中的位置和速度為狀態(tài)量，狀態(tài)方程依據(jù)接近段軌道動(dòng)力學(xué)模型建立，導(dǎo)航濾波算法采用擴(kuò)展卡爾曼濾波.針對(duì)前述幾種觀測(cè)模型，給出下述6種導(dǎo)航方案，其中hi(i=1,2,…，6)為第i種方案對(duì)應(yīng)觀測(cè)量.

3.1 基于無(wú)線電測(cè)量的自主導(dǎo)航方案

探測(cè)器在軌接收地面站發(fā)射的無(wú)線電信息，可獲取的信息主要包括探測(cè)器相對(duì)地面站的徑向速度和距離信息，則基于無(wú)線電的在軌導(dǎo)航方案觀測(cè)方程為

(23)

其中rr、vr如式(21)所示，分別為探測(cè)器相對(duì)地面站的距離和徑向速度.

3.2 基于光學(xué)測(cè)量的自主導(dǎo)航方案

為充分利用光學(xué)導(dǎo)航圖像獲取的信息進(jìn)行高精度導(dǎo)航，光學(xué)導(dǎo)航方案里應(yīng)同時(shí)包含方向和距離信息，因此對(duì)上述3種光學(xué)觀測(cè)量進(jìn)行組合可有以下2種組合導(dǎo)航方式：

(24)

其中，p、l、ρ如式(18)、(20)所示，分別為目標(biāo)天體中心像元、像線及目標(biāo)天體視半徑，ψs=[ψs1ψs2…ψsn]T為觀測(cè)到的n個(gè)星光角距信息.

對(duì)于這兩種基于光學(xué)測(cè)量組合的導(dǎo)航方案，前一種方案由觀測(cè)模型可看出其需要探測(cè)器的姿態(tài)信息和敏感器的安裝信息，因此其導(dǎo)航精度受姿態(tài)誤差和敏感器安裝誤差的影響；而后一種方案雖也是由目標(biāo)天體和恒星點(diǎn)像素信息得到，但在計(jì)算夾角信息時(shí)姿態(tài)誤差和安裝誤差引起的測(cè)量誤差被抵消，因此該方案不受這些因素的影響.

3.3 融合光學(xué)和無(wú)線電信息的導(dǎo)航方案

根據(jù)前述給出的光學(xué)導(dǎo)航方案，給出如下3種融合無(wú)線電和光學(xué)觀測(cè)信息的在軌導(dǎo)航方案：

(25)

4 接近段導(dǎo)航B平面參數(shù)誤差分析

由于建立在目標(biāo)天體B平面上的參數(shù)與探測(cè)器軌道狀態(tài)參數(shù)之間存在較好的線性關(guān)系，又可用來(lái)描述探測(cè)器目標(biāo)軌道的散布，因此B平面參數(shù)除了可用于深空探測(cè)的軌道設(shè)計(jì)和控制外，還可用于深空接近任務(wù)段導(dǎo)航性能的評(píng)估.

圖2 B平面及其誤差橢圓

如圖2所示，B平面一般定義為過(guò)目標(biāo)天體中心并垂直于進(jìn)入軌道的雙曲線漸近線的平面.B平面坐標(biāo)系定義為，以目標(biāo)天體中心為原點(diǎn)，S軸為進(jìn)入軌道漸近線方向，M為參考方向，此處取為目標(biāo)天體赤道面法線方向，S與M的矢量積為T軸，R軸由S軸和T軸按右手法則確定，B矢量定義為由B平面原點(diǎn)指向雙曲線軌道漸近線與B平面交點(diǎn)的矢量，記為B，其大小為雙曲線短半軸長(zhǎng)度，B平面參數(shù)BT和BR定義為B在T軸和R軸上的分量，即

BT=B·T,

BR=B·R

(26)

B平面靶點(diǎn)的導(dǎo)航不確定性可利用誤差橢圓表達(dá)，有關(guān)B平面參數(shù)和誤差橢圓的詳細(xì)內(nèi)容可參考文獻(xiàn)[9-10].因線性協(xié)方差分析方法可獲得與蒙特卡羅仿真相似的狀態(tài)誤差分布，因此可利用線性協(xié)方差分析方法確定的導(dǎo)航誤差協(xié)方差陣確定B平面誤差橢圓.

利用線性協(xié)方差分析方法可得到探測(cè)器在導(dǎo)航結(jié)束時(shí)刻tf的導(dǎo)航誤差協(xié)方差Pf(x).由于B平面參數(shù)與探測(cè)器軌道狀態(tài)參數(shù)之間有較好的線性關(guān)系，即ΔB=KΔx，其中ΔB=[ΔBTΔBR]，故有

PB=KPf(x)KT

(27)

其中K為B平面參數(shù)敏感矩陣，可通過(guò)數(shù)值方法求解，PB為B平面參數(shù)誤差協(xié)方差陣.令σT、σR分別為BT、BR的誤差標(biāo)準(zhǔn)差，ρRT為相關(guān)系數(shù)，則PB也可寫為

(28)

對(duì)于實(shí)對(duì)稱矩陣PB有

P′=UTPBU

(29)

式中，P′是以PB特征值λ1、λ2為主對(duì)角線元素的對(duì)角陣，U是由對(duì)應(yīng)特征向量的標(biāo)準(zhǔn)正交基構(gòu)成的矩陣.令[x′y′]T=UTΔB則有B平面誤差橢圓方程如下：

(30)

當(dāng)L分別取1、2、3時(shí)，得到橢圓分別為1σ、2σ、3σ誤差橢圓.設(shè)λ1≥λ2，則誤差橢圓長(zhǎng)短半軸及方向角可計(jì)算如下:

(31)

式中Ui1(i=1,2)為U中對(duì)應(yīng)元素.

接近天體點(diǎn)高度標(biāo)準(zhǔn)差可計(jì)算為

(32)

5 仿真分析

以接近火星探測(cè)任務(wù)為例，針對(duì)給出的6種接近天體在軌導(dǎo)航方案和基于線性協(xié)方差分析的B平面參數(shù)誤差分析方法進(jìn)行數(shù)學(xué)仿真.

5.1 仿真條件

1)標(biāo)稱軌道由STK生成，考慮火星80×80階引力場(chǎng)、太陽(yáng)引力、火衛(wèi)星引力和太陽(yáng)光壓攝動(dòng)；仿真采用簡(jiǎn)化動(dòng)力學(xué)模型，只考慮太陽(yáng)引力攝動(dòng)和火星非球形J2項(xiàng)攝動(dòng)，取系統(tǒng)模型誤差方差陣為

Q=diag{0,0,0,10-5,10-5,10-5}

2)探測(cè)器在火星J2000慣性系下初始位置、速度：

X0=[227326.99km;-160020.36km;

126502.87km;-2.58386631km/s;

1.8939265km/s;-1.5008243km/s]

3)導(dǎo)航初始誤差：位置各方向?yàn)?00km，速度各方向10m/s.

4)觀測(cè)誤差：無(wú)線電距離測(cè)量誤差100m(1σ)；無(wú)線電徑向速度測(cè)量誤差0.1m/s(1σ)；光學(xué)敏感器焦距677mm，光學(xué)觀測(cè)誤差0.1像素(1σ)，光學(xué)角度測(cè)量誤差0.01°(1σ).

5)姿態(tài)誤差各軸方向5″(1σ).

6)濾波方法EKF，濾波周期30s，導(dǎo)航仿真時(shí)間20h.

5.2 仿真結(jié)果及分析

圖3、4分別給出了方案4由線性協(xié)方差分析方法和蒙特卡羅仿真得到的導(dǎo)航誤差分布和B平面參數(shù)誤差分布結(jié)果.圖3中虛線為由線性協(xié)方差分析方法得到的狀態(tài)誤差標(biāo)準(zhǔn)差包絡(luò)線，實(shí)線為200次蒙特卡羅仿真結(jié)果，圖4中3個(gè)橢圓是由線性協(xié)方差分析方法得到的對(duì)應(yīng)1σ、2σ、3σ的B平面誤差橢圓，小圓圈代表1000次蒙特卡羅仿真得到的B平面打靶結(jié)果，由仿真結(jié)果可看出線性協(xié)方差分析方法得到結(jié)果與蒙特卡羅仿真結(jié)果基本一致.表1給出了無(wú)線電導(dǎo)航、光學(xué)導(dǎo)航及融合無(wú)線電和光學(xué)信息的導(dǎo)航方案B平面參數(shù)精度對(duì)比.圖5給出了不同導(dǎo)航方案的B平面誤差橢圓，圖6比較了不同姿態(tài)誤差下方案4的B平面誤差橢圓與方案5的B平面誤差橢圓.

表1 不同導(dǎo)航方案對(duì)應(yīng)的B平面參數(shù)誤差

由仿真結(jié)果可看出，光學(xué)導(dǎo)航精度遠(yuǎn)優(yōu)于無(wú)線電導(dǎo)航，融合光學(xué)信息和無(wú)線電信息的導(dǎo)航方案精度高于光學(xué)或無(wú)線電信息單獨(dú)進(jìn)行導(dǎo)航的精度；當(dāng)姿態(tài)誤差較小時(shí)，基于火心像素和火星視半徑信息的光學(xué)導(dǎo)航方式優(yōu)于基于星光角距和火星視半徑信息的光學(xué)導(dǎo)航，相應(yīng)的融合無(wú)線電信息后的導(dǎo)航精度前者也優(yōu)于后者，但隨著姿態(tài)誤差增大，融合無(wú)線電、星光角距和火星視半徑信息的導(dǎo)航方式精度逐漸高于融合無(wú)線電、火心像素和火星視半徑信息的導(dǎo)航方式精度，見(jiàn)圖6，這是由于火心像素觀測(cè)信息對(duì)姿態(tài)誤差極其敏感，而基于星光角距和火星視半徑信息的光學(xué)導(dǎo)航不受姿態(tài)誤差的影響，與前面分析一致；由圖5的結(jié)果可看出方案5比方案6的精度稍高但并不明顯，這主要是因?yàn)榛鹦且暟霃叫畔⒃诰嗷鹦禽^遠(yuǎn)時(shí)對(duì)探測(cè)器狀態(tài)變化不敏感所致.

圖3 線性協(xié)方差與蒙特卡羅仿真結(jié)果比較

因此當(dāng)探測(cè)器姿態(tài)誤差較小時(shí)，應(yīng)以目標(biāo)天體中心像素為主要的光學(xué)觀測(cè)信息，當(dāng)姿態(tài)誤差較大時(shí)，則應(yīng)以星光角距為主要的光學(xué)觀測(cè)信息；當(dāng)距離目標(biāo)天體較近時(shí)，可增加目標(biāo)天體視半徑為主要光學(xué)觀測(cè)信息.

圖4 方案4的B平面誤差橢圓和蒙特卡羅打靶結(jié)果

圖5 不同導(dǎo)航方案的B平面誤差橢圓

圖6 在不同姿態(tài)誤差下方案4與方案5 B平面誤差橢圓比較

6 結(jié) 論

本文結(jié)合線性協(xié)方差分析方法和基于B平面參數(shù)的誤差分析方法，分析了在軌融合無(wú)線電測(cè)量信息和光學(xué)觀測(cè)信息的自主導(dǎo)航方法.由線性協(xié)方差分析得到的B平面參數(shù)誤差表明，在軌融合無(wú)線電測(cè)量信息和光學(xué)觀測(cè)信息的導(dǎo)航方法可有效提高預(yù)測(cè)的B平面參數(shù)和目標(biāo)天體近心距精度；在考慮姿態(tài)誤差影響時(shí)，選擇合適的光學(xué)觀測(cè)信息和無(wú)線電測(cè)量信息進(jìn)行融合可提高系統(tǒng)導(dǎo)航精度和可靠性.

參 考 文 獻(xiàn)

[1]王大軼,黃翔宇.深空探測(cè)自主導(dǎo)航與控制技術(shù)綜述[J].空間控制技術(shù)與應(yīng)用,2009,35(3):6-12

Wang D Y, Huang X Y. Survey of autonomous navigation and control for deep space exploration[J]. Aerospace Control and Application, 2009, 35(3):6-12

[2]Bhaskaran S, Riedel J E, Synnott S P, et al. The deep space 1 autonomous navigation system: a post-flight analysis[C]. AIAA/AAS Astrodynamics Specialist Conference, Denver, CO, USA, Aug 2000

[3]Ryne M, Jefferson D, Craig D, et al. Ground-based orbit determination for deep impact[C]. The 16thAnnual AAS/AIAA Space Flight Mechanics Meeting, San Diego, CA, 2006

[4]黃翔宇,王大軼.基于光學(xué)和無(wú)線電信息的接近天體在軌導(dǎo)航方法研究[C].全國(guó)第十五屆空間及運(yùn)動(dòng)體控制技術(shù)學(xué)術(shù)會(huì)議,長(zhǎng)春,2012

Huang X Y, Wang D Y. Study on autonomous navigation based on optical and radio measurements for the approach phase of interplanetary mission[C]. The 15thAcademic Conference on Space and Movement Body Control Technology, Changchun, 2012

[5]黃翔宇.深空自主導(dǎo)航系統(tǒng)的可觀性分析[J].宇航學(xué)報(bào),2006,27(3):332-337

Huang X Y. Observability analysis of deep space autonomous navigation system[J]. Journal of Astronautics, 2006, 27(3): 332-337

[6]Geller D K, Rose M B, Woffinden D C. Event triggers in linear covariance analysis with applications to orbital rendezvous[J]. Journal of Guidance Control and Dynamics, 2009, 32(1): 102-111

[7]Geller D K. Linear covariance techniques for orbital rendezvous analysis and autonomous onboard mission planning[J]. Journal of Guidance Control and Dynamics, 2006, 29(6): 1404-1414

[8]Stastny N B, Geller D K. Autonomous optical navigation at Jupiter: a linear covariance analysis[J]. Journal of Spacecraft and Rockets, 2008, 45(2): 290-298

[9]Kizner W. A method of describing miss distances for lunar and interplanetary trajectories [R]. NASA CR-63483, 1959

[10]Byron D T, Bob E S, George H B. Statistical orbit determination[M]. USA: Elsevier Academic Press, 2004, 251-261