馮剛 楊東春 顏根廷

1.上海宇航系統(tǒng)工程研究所,上海,201109

跟瞄設(shè)備是空間編隊(duì)飛行中普遍使用的相對測量敏感器.目前,在相對導(dǎo)航濾波器設(shè)計(jì)中,普遍采用卡爾曼濾波以及擴(kuò)展卡爾曼濾波作為濾波估計(jì)算法.當(dāng)系統(tǒng)為線性并且噪聲統(tǒng)計(jì)特性滿足高斯分布時(shí),卡爾曼濾波是最小方差估計(jì)器.當(dāng)噪聲不滿足高斯分布,特別是噪聲中含有粗差時(shí),卡爾曼濾波的估計(jì)精度會大幅下降.由Huber提出的以其名字命名的濾波算法Huber濾波器,從魯棒統(tǒng)計(jì)的角度出發(fā),通過結(jié)合L1和L2范數(shù)估計(jì)器的特點(diǎn),其估計(jì)性能在測量誤差含有粗差的情況下具有較好的魯棒性.本文結(jié)合UKF(Unscented Kalman Filter)和Huber濾波算法作為編隊(duì)飛行中非線性狀態(tài)估計(jì)問題的相對導(dǎo)航濾波器.

空間攻防與維護(hù)過程中,追蹤航天器可能距離目標(biāo)航天器較近的情況下進(jìn)行編隊(duì)飛行,為了提高軌道控制的安全性,本文對追蹤航天器的軌控加速度偏差進(jìn)行了建模,采用Twisting算法設(shè)計(jì)了相對軌道運(yùn)動二階滑模變結(jié)構(gòu)控制器,該控制器對推力故障引起的控制加速度偏差具有一定容忍度,且算法比較簡單,具有工程參考價(jià)值.

1 相對軌道動力學(xué)描述

為解決空間操控時(shí)相對路徑規(guī)劃控制的需要,在動力學(xué)模型建立中,引入“編隊(duì)向量”參數(shù),如圖1所示.

圖1 編隊(duì)飛行相對幾何關(guān)系

圖1中是編隊(duì)向量,它描述了在目標(biāo)LVLH坐標(biāo)系(當(dāng)?shù)卮怪彼阶鴺?biāo)系)下的編隊(duì)圖案的變化,將“編隊(duì)向量”設(shè)計(jì)為時(shí)間變量的函數(shù),可形成滿足任務(wù)需求的標(biāo)稱運(yùn)動軌跡,與自主閉環(huán)控制方法相結(jié)合,可實(shí)現(xiàn)對目標(biāo)航天器的逼近、伴飛、繞飛等不同的機(jī)動任務(wù);T、S分別表示目標(biāo)航天器和在追蹤航天器在J2000坐標(biāo)系下的位置矢量.

在目標(biāo)航天器LVLH坐標(biāo)系內(nèi),相對軌道動力學(xué)方程可寫為:

式(1)～式(3)可以寫成如下形式

2 跟瞄相對導(dǎo)航濾波器設(shè)計(jì)

假設(shè)跟瞄測量坐標(biāo)系與追蹤航天器的體軸系重合,測量輸出為:ρ?目標(biāo)航天器相對于追蹤航天器的視線距離;˙ρ?視線距離變化率;α?視線在跟瞄測量坐標(biāo)系中的高低角,定義為視線與其在測量坐標(biāo)系xy平面的投影間的夾角,抬頭為正;˙α?視線高低角變化率;β?視線在跟瞄測量坐標(biāo)系中的方位角,定義為視線在測量坐標(biāo)系xy平面的投影與x軸的夾角,偏向+y軸方向?yàn)檎?˙β?視線方位角變化率.

令目標(biāo)航天器在追蹤航天器體軸系中的相對位置坐標(biāo)為(xb,yb,zb),其時(shí)間變化率為(,,).考慮到跟瞄測量坐標(biāo)系與追蹤航天器的體軸系重合,相對位置和速度在追蹤航天器體軸系下的分量與跟瞄測量值之間的換算關(guān)系為:

圖2 追蹤航天器星載跟瞄設(shè)備測量坐標(biāo)系示意圖

令向量Xb=表示目標(biāo)航天器在追蹤航天器體軸系下的位置和速度,令向量

相對軌道動力狀態(tài)方程模型中定義的狀態(tài)量X與向量Xb之間有如下關(guān)系成立:

考慮到實(shí)際情況,跟瞄測量的˙α、˙β精度不高,因此將跟瞄測量信號中的ρ、α、β、˙ρ作為觀測量,得到觀測方程

式中wρ、wα、wβ、w˙ρ為測量噪聲,綜合上式可以得到觀測方程

從相對軌道動力學(xué)模型及跟瞄測量方程可見,當(dāng)目標(biāo)軌道為橢圓軌道時(shí),相對軌道動力學(xué)模型和跟瞄測量方程是時(shí)變、強(qiáng)非線性的.Unscented濾波(UKF)是用一系列確定的采樣點(diǎn)來逼近狀態(tài)的后驗(yàn)概率密度,在高斯白噪聲環(huán)境下,對非線性系統(tǒng)有良好的跟蹤性能.假設(shè)UKF經(jīng)過時(shí)間更新后,已經(jīng)得到交互協(xié)方差PXZ和狀態(tài)預(yù)測協(xié)方差PK+1|K,定義K+1時(shí)刻狀態(tài)真值和預(yù)測值的關(guān)系如下

預(yù)測誤差δXK+1的方差為PK+1|K.對觀測方程進(jìn)行線性化,其斜率矩陣為

量測方程可近似為

構(gòu)造線性回歸模型

定義下述各量

ξ的協(xié)方差陣為單位矩陣,并有

定義Huber方法的代價(jià)函數(shù)

其中γ為可調(diào)參數(shù),一般取1.345.

Huber濾波器取決于下式求解:

其中

具體由下式給出

容易得到

采用迭代方法求解,具體求解過程為

其中迭代初值由下式確定:

迭代初值也可以采用UKF測量更新后的估計(jì)值.一般只需要迭代一步就可以達(dá)到滿意的收斂精度.迭代初值結(jié)束后求得估計(jì)的方差為

其中Ψ取迭代最后一步的值.

3 繞飛標(biāo)稱軌跡設(shè)計(jì)

以編隊(duì)飛行時(shí)的主動繞飛為例,繞飛的標(biāo)稱參考軌跡用“編隊(duì)向量”表示.如圖3所示,假設(shè)以目標(biāo)航天器為圓心的一段弧線AD,從LVLH坐標(biāo)系y軸負(fù)方向開始度量,弧線起始點(diǎn)的對應(yīng)的角度為Θ0,弧線終端點(diǎn)的對應(yīng)的角度為Θf,繞飛總時(shí)間為T.

圖3 LVLH坐標(biāo)系下繞飛弧線軌跡

繞飛過程中,把運(yùn)動軌跡分為3段:加速段AB、勻速段BC和減速段CD.

圖4 繞飛各段路徑示意圖

加速段對應(yīng)的起始角θ(t01),對應(yīng)的終止角θ(tf1),起始時(shí)刻t01=0,加速段運(yùn)行的結(jié)束時(shí)刻為tf1=0.2T.加速段的圓心角邊值條件如下:

勻速段對應(yīng)的起始角θ(t02),對應(yīng)的終止角θ(tf1),起始時(shí)刻t02=0.2T,結(jié)束時(shí)刻為tf2=0.8T,勻速段運(yùn)行的角速率為ω=(θf2?θ02)/tf2.勻速段的圓心角邊值條件如下:

減速段對應(yīng)的起始角θ(t03),對應(yīng)的終止角θ(tf3),減速段起始時(shí)刻t03=0.8T,結(jié)束時(shí)刻為tf3=T.減速段的圓心角邊值條件如下:

選擇加速段和減速段的θ(t)的形式為如下多項(xiàng)式形式:

其中

可以驗(yàn)證,上述多項(xiàng)式形式滿足2點(diǎn)邊值條件.

選擇勻速段的θ(t)的形式為如下線性形式:

求出各段的圓心角變化函數(shù)后,目標(biāo)LVLH坐標(biāo)系下的編隊(duì)向量設(shè)計(jì)為

式中,ρ為繞飛半徑,θ為表征繞飛位置的圓心角度參數(shù).

上述分段邊值函數(shù)設(shè)計(jì),保證了繞飛時(shí)從起始點(diǎn)逐漸加速,建立起繞飛速度,在到達(dá)終端點(diǎn)時(shí)候又逐漸減速為0,這樣可以讓繞飛各階段時(shí)的軌跡光滑銜接.

4 相對軌道二階Twisting變結(jié)構(gòu)控制器設(shè)計(jì)

前面基于編隊(duì)向量設(shè)計(jì)了繞飛時(shí)的標(biāo)稱運(yùn)動軌跡,控制器的任務(wù)就是消除實(shí)際軌跡與標(biāo)稱軌跡之間的偏差.考慮編隊(duì)飛行動力學(xué)系統(tǒng)式(1)～式(3),假設(shè)追蹤航天器的推力加速度偏差時(shí),相對軌道動力表示為:

其中,

動力學(xué)方程增加的γ(t?T)項(xiàng)描述了推力故障情況下引起的加速度偏差和擾動加速度是有界的.

其中,εi、ηi為某個已知的正常數(shù).

階躍函數(shù)γ(t?T)描述了推力加速度偏差發(fā)生的時(shí)間,γ(t?T)的表達(dá)式如下:

令?是追蹤航天器對目標(biāo)進(jìn)行編隊(duì)飛行時(shí),相對編隊(duì)點(diǎn)的位置偏差矢量,令=?表示跟蹤誤差.設(shè)計(jì)切換函數(shù)為

滑動模態(tài)的形式為

其中,s:R+×R3→R3.

假設(shè)目標(biāo)無機(jī)動,忽略攝動加速度,不考慮推力故障時(shí),計(jì)算s的二階導(dǎo)數(shù)得到

其中,

把式(62)代入到式(60)中,得到

滿足如下條件

式(64)和式(65)表明

基于Twisting算法,設(shè)計(jì)編隊(duì)飛行的二階滑模變結(jié)構(gòu)控制為(當(dāng)推力器無故障情況)為:

其中,滿足收斂的條件為

證明過程如下:

根據(jù)式(68),有

結(jié)合式(67),可知

式(71)對應(yīng)的等效微分方程為

圖5 s?相平面收斂軌跡示意圖

如圖5,假設(shè)滑模態(tài)變量在t=0時(shí)刻的初始值為si=si0>0、=0.在相平面上,軌跡會進(jìn)入<0的下半平面,由微分方程式(72),相軌跡會進(jìn)入<0的下半平面的軌跡為其中,si1為相軌跡接下來與軸=0相交的坐標(biāo).顯然在ki2>ηi條件下,有

同樣的,相軌跡進(jìn)入si<0的左半平面或si>0的右半平面時(shí),有不等式成立

綜合式(38)、式 (39),表明在有限時(shí)間內(nèi)s→0→0.

進(jìn)行積分,可以得到估算收斂時(shí)間的計(jì)算式為顯然,在收斂條件ki1?ki2>ηi成立情況下,式(77)是有意義的.

當(dāng)考慮追蹤航天器的推力故障時(shí),與前面推導(dǎo)過程類似,可得到

按照同樣的過程,設(shè)計(jì)對推力故障具有容忍性的控制器為

其中,滿足收斂的條件為

為了削弱滑模變結(jié)構(gòu)控制器的抖振影響,引入飽和函數(shù)來代替式(79)中的符號函數(shù),其中飽和函數(shù)的定義為:

圖6 目標(biāo)LVLH坐標(biāo)系下相對標(biāo)稱軌跡的偏差

圖7 跟瞄相對導(dǎo)航誤差

圖8 追蹤飛行器體軸系下控制加速度

其中,?為正的常數(shù),稱為邊界層的厚度.

圖9 軌道面內(nèi)的繞飛軌跡

5 數(shù)值仿真驗(yàn)證

5.1 仿真輸入

假設(shè)目標(biāo)飛行器的初始軌道根數(shù)如表1所示.

表1 目標(biāo)飛行器的初始軌道根數(shù)

假設(shè)追蹤航天器對目標(biāo)進(jìn)行近距離軌道面內(nèi)快速繞飛,繞飛半徑r=200m,繞飛周期T=2000s;跟瞄測量噪聲均方差為σρ=2m、σα=0.06?、σβ=0.06?、σ˙ρ=0.1m/s;設(shè)計(jì)Twisting滑?？刂破鞯膮?shù)為ki1=2.0×10?3,ki2=1.0×10?3,ki3=4.0×10?3,ki4=2.0×10?3;飽和函數(shù)的邊界層的厚度設(shè)置為:位置φr=0.5,速度φ˙r=0.05.假設(shè)由于推力器故障,推力加速度下降了40%.為了對控制效果進(jìn)行對比,本文還用PD控制器進(jìn)行了仿真.

圖10 目標(biāo)LVLH坐標(biāo)系下相對標(biāo)稱軌跡的偏差

圖11 跟瞄相對導(dǎo)航誤差

圖12 追蹤飛行器體軸系下控制加速度

5.2 PD控制器結(jié)果

PD控制器位置控制誤差12m,繞飛共需速度增量?V=9.023m/s.

5.3 Twisting變結(jié)構(gòu)控制器結(jié)果

Twisting控制器位置控制誤差2m,繞飛共需速度增量?V=17.65m/s.

圖13 軌道面內(nèi)的繞飛軌跡

6 結(jié)論

本文以近程編隊(duì)飛行時(shí)的繞飛任務(wù)為例,進(jìn)行了詳細(xì)的跟瞄相對導(dǎo)航和軌跡控制方案設(shè)計(jì),并通過數(shù)學(xué)仿真進(jìn)行了驗(yàn)證.仿真表明:本文設(shè)計(jì)的相對導(dǎo)航和軌跡控制方案滿足繞飛任務(wù)的功能和性能要求,在控制加速度大小下降40%的情況下,PD控制的位置誤差約12m,而Twisting滑模變結(jié)構(gòu)控制器的位置誤差約2m,且Twisting滑模變結(jié)構(gòu)控制器需要的速度增量比PD控制器增加不多,但從安全性上更好,這對近距離編隊(duì)飛行尤為重要.