杜榮華, 張 翔, 廖文和

(南京理工大學機械工程學院, 江蘇 南京 210094)

0 引 言

僅測角相對導航并不是一個新概念,最早應用于19世紀科學家通過望遠鏡觀測天體并計算其軌道的研究中。當前,僅測角相對導航已經(jīng)在航海[1]、深空探測[2]、衛(wèi)星軌道確定[3]以及編隊飛行[4-5]任務中得到了廣泛應用和研究。

僅測角相對導航的概念相當簡單,即通過單個相機在一段時間內(nèi)測量主動星與空間目標之間的視線角(line of sight,LoS)就可以確定其相對運動狀態(tài)[6]。目前,常用來測量空間物體之間相對運動狀態(tài)的有源主動式傳感器(例如微波雷達、激光雷達等),由于其功耗高、體積大和造價貴,使得其在實際工程應用中受到很大限制,而光學相機等無源被動式傳感器(例如紅外相機、可見光相機等)在這一方面具有很大的優(yōu)勢[7-9],因為可應用于各種衛(wèi)星間距,并且對衛(wèi)星質量、功率的設計影響很小[10],非常適合微小衛(wèi)星平臺的應用。實際上,許多航天器配備了星敏感器,如果方向適當,可用于跟蹤其視野范圍內(nèi)的空間目標,并根據(jù)LoS觀測值執(zhí)行相對導航任務[11]。其中,僅測角相對導航技術代表了多種高級分布式空間系統(tǒng)的明確使能技術,包括自主交會和對接、空間態(tài)勢感知、分布式孔徑科學以及非合作目標的在軌維修等任務[12]。此外,僅測角相對導航技術在天基反導監(jiān)測以及空間碎片探測和清理等領域也有較大的應用前景[13-14]。

在過去十年內(nèi),已經(jīng)有在軌項目利用僅測角相對導航技術進行了非合作目標交會接近的演示驗證。瑞典的原型研究儀器和太空任務技術進步(prototype research instruments and space mission technology advancement,PRISMA)計劃是使用僅測角相對導航技術實現(xiàn)非合作目標遠距離自主交會的典型代表。2011年4月,PRISMA主承包商OHB-Sweden在自主交會(autonomous rendezvous, ARV)實驗中實現(xiàn)了從30 km至50 m的基于LoS的自主會合[15-17]。2011年8月,在PRISMA運行即將結束時,德國宇航中心(Germany Space Operation Center,GSOC)操作追蹤星Mango從60 km抵近目標星TanDEM-X至4 km處,抵近過程仍然使用基于LoS的相對導航濾波器[18-21]。2012年4月,GSOC使用全球定位系統(tǒng)和光學導航進行了高級交會論證(簡稱為ARGON)實驗[22],在ARGON實驗中,地面研發(fā)的專用飛行動力學系統(tǒng)可對星上收集的圖像進行處理,并根據(jù)單目相機測得的LoS矢量估計追蹤星和空間目標之間的相對軌道信息,進而實施一系列的軌道機動,在5天內(nèi)完成了從30 km至3 km的抵近操作。2016年11月,GSOC進行了自主視覺接近導航和目標識別(autonomous vision approach navigation and target identification,AVANTI)實驗[23],與ARGON實驗相比,AVANTI實驗中的圖像處理、相對導航以及軌道機動規(guī)劃均為星載自主實現(xiàn),而ARGON實驗中需要將追蹤星拍攝的圖像數(shù)據(jù)傳回地面,在地面處理完成后再將軌道機動指令發(fā)送到衛(wèi)星上。AVANTI實驗共進行了兩次基于僅測角相對導航的抵近操作,第一次抵近操作中將兩顆衛(wèi)星距離從13 km減少至1 km,然后依靠空氣阻力將衛(wèi)星間距拉大至2.5 km處,在第二抵近操作中將兩顆衛(wèi)星的距離減小至50 m以內(nèi)[23]。

由于較差的動態(tài)可觀性,AGRON和AVANTI兩個實驗都使用了先驗初始解來初始化僅測角相對導航濾波器,即北美航空航天防御司令部(簡稱為NORAD)提供的兩行根數(shù)(two line elements,TLE)和地面雷達測軌信息來確保僅測角相對導航濾波器的收斂[22-23]。從嚴格意義上講,整個相對導航過程并不是完全自主的。為了提高僅測角相對導航系統(tǒng)的自主性,研究僅測角相對導航初始相對軌道確定(initial relative orbit determination, IROD)方法是實現(xiàn)非合作目標完全自主抵近的關鍵。然而,由于單目相機深度信息的缺失,僅測角相對導航存在星間距離不可觀測的問題[24-25]。目前,已經(jīng)有許多學者研究了提高僅測角相對導航可觀性的方法,例如軌道機動法[26-27]、相機偏置法[28-29]、非線性動力學法[30-31]和多星多傳感器法[32-33],其中相機偏置法只能應用于近距離范圍內(nèi),而軌道機動法會增加衛(wèi)星的燃耗,多星多傳感器法則存在多衛(wèi)星多傳感器配置復雜的缺陷。

針對上述各種方法存在的問題,本文擬在無機動操縱觀測弧內(nèi),引入非線性相對運動模型,該非線性相對運動模型基于相對軌道根數(shù)(relative orbit elements, ROE)建立追蹤星和空間目標之間的相對運動狀態(tài),并將星間距離和相對軌道形狀進行解耦。然后,在線性理論獲得的共線性解附近系統(tǒng)地改變星間距離大小,并執(zhí)行一系列最小二乘擬合,隨后利用二分法或牛頓迭代法快速在全局范圍內(nèi)尋找最小擬合殘差的最優(yōu)解。最后,在4種不同的軌道場景中對提出的方法進行了半物理仿真校驗,仿真結果證明所提出的方法可以較快地獲得高精度的僅測角相對導航IROD的解。

1 問題描述

僅測角相對導航的目的在于推導出空間目標和主動星之間的相對運動狀態(tài),該相對運動狀態(tài)對應于ti時刻獲得的一組LoS觀測值{ui:i=1, 2,…,k},這些觀測值是指向空間目標的單位LoS矢量。一般情況下,僅測角相對導航的問題與系統(tǒng)的可觀性密切相關。Woffinden和Geller[34]已經(jīng)證明,在線性化的測量模型和相對運動模型的假設下,僅測角相對導航系統(tǒng)是不可觀測的。

首先,令x(t)表示空間目標和主動星在時間t處的相對運動狀態(tài),則線性化的相對運動模型的演化方程可表示為

x(t)=Φ(t,t0)x(t0)

(1)

式中,Φ(t,t0)是狀態(tài)轉移矩陣。

假設相對位置矢量r(t)與相對運動狀態(tài)x(t)之間具有線性相關性,其關系可表示為

r(t)=L(t)x(t)

(2)

一組給定的LoS觀測值{ui:i=1, 2,…,k}和相對位置矢量r(t)之間滿足以下關系:

ui×r(ti)=0,i∈[1,k]

(3)

將式(1)和式(2)代入式(3)中,得到

ui×(L(ti)Φ(ti,t0)x(t0))=0,i∈[1,k]

(4)

從式(4)中可以看出,如果x(t0)=x0是式(4)的解,則縮放解φx0也是式(4)的解,從而導致一組給定的LoS觀測值具有無限個解,造成系統(tǒng)不可觀測,這個被稱為Woffinden困境[35]。因此,為了通過僅測角相對導航獲得空間目標和主動星之間的相對運動狀態(tài),必須解決Woffinden困境,提高僅測角相對導航系統(tǒng)的可觀性。

2 算法說明

由第1節(jié)的Woffinden困境可知,采用線性化的相對運動模型和測量模型將導致僅測角相對導航系統(tǒng)不可觀測,因此本文引入了非線性相對運動模型。首先,為捕獲線性模型所忽略的細微差異,重新建立相對運動狀態(tài)更新方程,可記為

x(t)=f(t,x(t0))

(5)

式中,f是非線性函數(shù),對應于包含J2攝動的相對運動方程的數(shù)值積分。

測量模型h可表示為

(6)

式中,r(t,x0)為t時刻空間目標相對于主動星在直角笛卡爾坐標系中的相對位置矢量?？闪顁(t,x0)=g(x(t)),其中g一般為非線性轉換函數(shù)。

假設ui=[ux(ti),uy(ti),uz(ti)]是由相機采集的圖像獲得的空間目標在ti時刻的LoS單位矢量,該LoS單位矢量也可用一組角度測量值(即方位角α和俯仰角ε)進行表示,其中i=1, 2,…,k,表示第i個測量時刻。理想情況下,矢量ui與h(ti,x0)在方向上是平行的,即

ui×h(ti,x0)=0,i∈[1,k]

(7)

由于ui中含有LoS矢量測量誤差,h(ti,x0)含有建模誤差,兩者實際有偏差,導致ui×h≠0。因此,IROD問題可以描述為最優(yōu)化問題:

(8)

式中,J為損失函數(shù),表示前k個時刻LoS矢量誤差的平方和。

由于僅測角導致的弱可觀性,如何設計優(yōu)化算法保證式(8)收斂到全局最小值是找到精確IROD解的關鍵。在線性化的相對運動模型中,空間目標相對于主動星之間的相對運動狀態(tài)可以表示為

x(t)=Φ(t,t0)x0

(9)

假設相對位置矢量r(t,x0)和x(t)之間具有線性相關性,可記為

r(t,x0)=C(t)x(t)

(10)

式中,C(t)為r(t)與x(t)之間的轉換矩陣。

對于k個LoS矢量測量值ui,由式(7)可構造以下線性方程組:

ui×(C(ti)Φ(ti,t0)x0)=0,i∈[1,k]

(11)

(13)

比較式(8)與式(13)的尋優(yōu)過程,可見尋優(yōu)搜索空間從6個維度的x0減少到只有一個維度的φ。圖1給出了該IROD研究方法的流程圖。需要指出的是,縮放解φx0中φ的上限和下限取決于主動星和空間目標之間的距離范圍。對于本研究的應用場景,主動星和空間目標的間隔大約介于1～100 km之間,通過在1～100 km范圍內(nèi)改變比例因子φ,就可以找到使擬合殘差均方根σ(φ)最小的比例因子φ,即為僅測角相對導航IROD的全局最優(yōu)解。

圖1 僅測角相對導航IROD方法流程框圖

證明σ(φ)為凸函數(shù)

實際上,可以通過使用一組ROE或一組笛卡爾坐標系的坐標分量來參數(shù)化相對運動狀態(tài)x。在這種情況下,相對運動模型f描述了相對運動狀態(tài)x的時間演變,而函數(shù)g將相對運動狀態(tài)x映射為直角笛卡爾坐標系中的坐標分量。相對運動狀態(tài)x的建?？刹捎们€型Hill-Clohessy-Wiltshire(HCW)模型[36]和Gim-Alfriend模型[37]。后者基于ROE,但提供了一個線性映射L,以將ROE轉換為曲線笛卡爾坐標系表示的相對運動狀態(tài)(即x(t)=L(t)O(t)ω(t,t0)O(t)L(t)x0=Φ((t,t0)x0),其中ω為基于平均ROE的狀態(tài)轉移矩陣,O為平均軌道根數(shù)和密切軌道根數(shù)之間的轉換矩陣),可以根據(jù)Brouwer[38]和Lyddane[39]中提出的理論進行平均軌道根數(shù)至密切軌道根數(shù)的一階非線性映射。

根據(jù)式(3),可以描述一般非線性的測量方程為

ui×r(ti)=ui×(g(ti,f(ti,x0)))=0,i∈[1,k]

(14)

進一步簡化為

ui×(g(ti,f(ti,x0)))=0,i∈[1,k]

(15)

假設曲線型相對運動狀態(tài)在軌道坐標系中進行表示,并且分量按徑向-切向-法向序列排序。對于正在考慮的問題(即相對運動主要為沿軌道方向的星間距離),直線相對位置r可以由曲線相對位置rcur得到

(16)

式中,R是圓形軌道的半徑;rcur,2表示相對位置沿軌道方向的分量,且R≥rcur,2,所以式(16)可以簡化為

(17)

把狀態(tài)轉移矩陣拆分成兩塊,得到Φ=[Φ1-3,1-6,Φ4-6,1-6]T。讓Φ6,1-6定義矩陣Φ的第6行,且r=g(x)=g(Φx0)=Φ1-3,1-6x0。因此,式(15)可以采取以下形式進行表示:

(18)

因此,對式(15)的一般非線性方程,該問題已簡化為二次函數(shù)的最小化問題,為簡單起見,引入與衛(wèi)星軌道曲率模型相對應的二次部分:

(19)

ui×(Φ1-3,1-6x0)=0,i∈[1,k]

(20)

并擴展成

(21)

式中,Φ1-3,1-6=[Φ1-3,1-5,Φ1-3,6]已拆分為多個塊。

累加k次測量后,獲得線性系統(tǒng)(由表示線性問題的下標L標識):

(22)

其五維解可由最小二乘法給出:

(23)

0,i∈[1,k]

(24)

使用式(21)中采用的相同方法將問題限制在五維空間中:

(25)

(26)

在累加k個測量值之后,再次獲得另一個線性系統(tǒng):

(27)

Γ(φ)T(I-P(φ))Γ(φ)

(28)

為簡化起見,引入投影矩陣P(φ)=A(φ)(A(φ)T·A(φ))-1A(φ)T,得出任意固定為φ時損失函數(shù)J的最小二乘意義的解析公式,但其凸性仍有待證明?？梢酝ㄟ^考慮以下因素來簡化式(25)以便近似得

(29)

(30)

進而矢量?？梢员磉_成以下形式:

Γ= sign(φ)(Γ1+φΓ2)

(31)

(32)

(33)

3 數(shù)值分析

3.1 相對運動狀態(tài)參數(shù)化

引入一組無量綱的ROE來進行相對運動狀態(tài)的參數(shù)化[40]:

δα=[δa,δex,δey,δix,δiy,δλ]T

(34)

式中,δa是無量綱的相對半長軸;δλ代表相對平均經(jīng)度;δe=[δex,δey]T和δi=[δix,δiy]T分別為相對偏心矢量和相對傾斜矢量。

與笛卡爾坐標系表示的相對運動狀態(tài)參數(shù)相比,通過ROE進行參數(shù)化具有以下優(yōu)點:可以快速了解相對運動的幾何形狀,如圖2所示,描述了局部笛卡爾坐標系(徑向、切向和法向或RTN)中表示的相對運動狀態(tài),其單位矢量定義如下:eR和eN分別與主動星的絕對位置和軌道角動量方向對齊,eT與eR和eN滿足右手定則。

此外,從圖2可以看出[35],相對運動的大小可以通過ROE的尺寸(即軌道半長軸a)來描述[41]。

圖2 笛卡爾坐標系和ROE表示的相對運動

這組ROE特別適用于僅測角相對導航問題,因為弱可觀測的星間距離幾乎與aδλ分量是相等的,從而將相對運動的形狀和星間距離從幾何上進行了解耦。

3.2 相對運動模型

需要建立主動星和空間目標相對運動的精確模型,以捕獲系統(tǒng)微弱的非線性效應,提高僅測角相對導航系統(tǒng)的可觀性。本文使用第3.1節(jié)中介紹的ROE(即x=δα)對相對運動模型進行參數(shù)化。在這種情況下,g不再是線性的。為了獲得笛卡爾坐標系中的相對位置,首先需要將平均ROE轉換為密切ROE,從而恢復J2攝動引起的長期和短期效應,然后將密切相對軌道根數(shù)映射到笛卡爾坐標系表示的相對位置,并且允許在相對運動模型中引入其他經(jīng)驗參數(shù),例如反映差動阻力擾動效應的參數(shù)[42-43]。

在數(shù)學上,該模型通過狀態(tài)轉移矩陣Φ將時間t處的相對運動狀態(tài)δα(t)與時間t0處的相對運動狀態(tài)δα(t0)相關聯(lián),即

δα(t)=Φ(t,t0)δα(t0)

(35)

(36)

(37)

(38)

3.3 IROD的快速求解方法

由于使用ROE進行相對運動狀態(tài)的參數(shù)化,因此矩陣C表示的是ROE和直角笛卡爾坐標系表示的相對位置之間的轉換關系[35]:

C(t)=

(39)

式中,u代表主動星的緯度幅角;i代表主動星的軌道傾角。

在最小二乘擬合殘差意義上對應的特定值S=aδλ=φ的線性解可推導如下:

(40)

Γ=-SA6

(41)

(42)

(43)

式中,κi代表κ的第i個分量。

3.3.1 二分法

3.3.2 牛頓迭代法

牛頓迭代法的思想是將非線性函數(shù)(原方程)線性化(切線方程),以線性方程的解逐步逼近非線性方程的解,其原理圖如圖3所示。

圖3 牛頓迭代法原理圖

其原理公式為

(44)

牛頓迭代法的具體步驟如下:

4 仿真驗證

4.1 地面半物理仿真平臺搭建

在完成快速僅測角相對導航IROD算法的設計后,搭建驗證該算法的地面半物理仿真平臺,其結構框圖如圖4所示。

圖4 地面半物理仿真平臺結構框圖

該地面半物理仿真平臺主要包括以下部分。

(1)軌道/姿態(tài)發(fā)生器

該模塊內(nèi)部裝有Matlab/Simulink軟件和STK軟件,主要用于生成主動星和空間目標的軌道/姿態(tài)數(shù)據(jù)。

(2)星光星圖模擬計算機

該計算機作為星光星圖模擬器,內(nèi)部裝有導航星庫、星圖生成軟件等,主要用來模擬恒星和空間目標在星圖中的分布情況及星光等級。

(3)星圖識別和目標檢測計算機

該計算機內(nèi)部裝有星圖識別軟件、星體質心定位軟件和空間目標自主檢測軟件等,主要功能是在原始圖像中找到空間目標,并輸出空間目標相對于攝像機的LoS測量信息。

(4)僅測角相對導航計算機

該計算機的主要功能是完成僅測角相對導航算法的實時運算,并輸出空間目標相對于主動星的相對運動狀態(tài)信息。

(5)實時監(jiān)視顯示計算機

該計算機主要用于監(jiān)視和顯示軌道/姿態(tài)發(fā)生器模擬的主動星和空間目標的軌道/姿態(tài)數(shù)據(jù)以及僅測角相對導航計算機的輸出結果等。

搭建的實際地面半物理仿真平臺如圖5所示。在原有設備的基礎上,配備了三軸精密轉臺和三軸直線導軌等機械微調裝置,可調整星圖顯示裝置和星圖采集模塊同軸。

圖5 地面半物理仿真平臺

4.2 實驗測試與結果分析

在搭建的地面半物理仿真平臺上對該方法進行校驗,主要仿真參數(shù)如表1所示。

表1 主要仿真參數(shù)

該仿真校驗在4種不同的軌道場景中進行,包括遠距離保持點相對運動軌跡(場景1)、自然漂移相對運動軌跡(場景2)、遠近距切換保持點相對運動軌跡(場景3)和主從隨動相對運動軌跡(場景4)。這4個軌道場景剛好包含了主動星從遠距離交會至抵近空間目標的各個過程,其相對運動軌跡如圖6所示。出于模擬的需要,主動星和空間目標的參考運動使用20×20階的重力場模型的數(shù)值積分傳播得到,其中包括第三體引力、太陽輻射壓力和大氣阻力擾動等,因此可以得到LoS測量的參考值以及主動星和空間目標之間無機動的觀測弧,最終根據(jù)主動星和空間目標模擬的狀態(tài)創(chuàng)建一組觀測值。

圖6 4種不同軌道場景的相對運動軌跡圖

此外,為了評估快速僅測角相對導航IROD算法的運行時間,僅測角相對導航計算機選用的是常規(guī)配置的Intel Core i5系列計算機。測試結果如表2～表3和圖7～圖14所示,為該算法在4種不同軌道場景中的性能。其中,表2為不同方法單次運算時耗;表3為采樣周期為30 s,采樣次數(shù)為5次時,不同方法的百分比誤差;圖7～圖14為不同采樣周期和觀測次數(shù)情況下,各方法在不同軌道場景中的IROD性能。

圖7 不同觀測次數(shù)下的IROD方法的性能(場景1,T=10 s)

圖8 采樣周期下的IROD方法的性能(場景1, 觀測次數(shù)為5)

圖9 不同觀測次數(shù)下的IROD方法的性能(場景2, T=10 s)

圖10 不同采樣周期下的IROD方法的性能(場景2,觀測次數(shù)為5)

圖11 不同觀測次數(shù)下的IROD方法的性能(場景3, T=10 s)

圖12 不同采樣周期下的IROD方法的性能(場景3,觀測次數(shù)為5)

圖13 不同觀測次數(shù)下的IROD方法的性能(場景4, T=10 s)

圖14 不同采樣周期下的IROD方法的性能(場景4,觀測次數(shù)為5)

表3 IROD方法的百分比誤差

從表2中可以看出,在相同軌道場景中,批量搜索法單次運算時耗大于牛頓迭代法和二分法,并且牛頓迭代法單次運算時耗最低。在不同軌道場景中,場景4的單次運算時耗明顯高于其他3種軌道場景,其次是場景2,場景1和場景3兩種軌道場景的單次運算時間相差不大,這是由于場景1和場景3的軌道場景相似,只是距離不同導致的。

表2 單次運算時耗

為了評估算法性能,定義均方根誤差、標準差和百分比誤差分別為

(45)

(46)

(47)

本文將相對平均經(jīng)度aδλ的均方根誤差Mx、標準差σx和百分比誤差Px作為評價算法性能的指標,因為僅測角相對導航的不可觀測性主要集中在該參數(shù)上[45]。

分析表3、圖7、圖9、圖11和圖13,在采樣周期相同但觀測次數(shù)不同的情況下,在相同軌道場景中,牛頓迭代法的IROD性能最好,其相對定軌誤差可達所測距離的6%左右,其次是二分法和牛頓迭代法。此外,不同的觀測次數(shù)也會影響算法的IROD性能,但并不是次數(shù)越多性能越好,在觀測次數(shù)大于一定數(shù)量后,算法性能趨于穩(wěn)定。

分析表3、圖8、圖10、圖12和圖14,在觀測次數(shù)相同但采樣周期不同的情況下,采樣周期會影響算法的IROD性能,其中采樣周期越大,算法的IROD性能越好。這是由于大的采樣周期可以增加系統(tǒng)的非線性和擾動效應,并抵消相機噪聲產(chǎn)生的影響,從而增加了僅測角相對導航的可觀性。比較不同軌道場景中算法的IROD性能,在場景4的軌道場景中,各方法的IROD性能明顯低于其他3種軌道場景,其相對定軌誤差都在所測距離的30%以上。這是由于場景4的系統(tǒng)可觀性較差導致的,這說明主動星和空間目標之間的相對運動軌跡也會影響所提算法的IROD性能。

5 結 論

本文針對僅測角相對導航IROD問題,提出了一種快速僅測角相對導航IROD方法,并通過地面半物理仿真平臺測試了該方法在4種不同軌道場景中的性能,得到以下結論。

(1)采用基于ROE建立的非線性相對動力學模型,通過最小二乘擬合殘差的方法可以獲得較高精度的僅測角相對導航IROD的解,且性能受采樣周期和觀測次數(shù)的影響。增大采樣周期可以提高所提方法的IROD的性能,這是由于大的采樣周期可以增加系統(tǒng)的非線性和擾動效應,并抵消相機噪聲產(chǎn)生的影響,從而提高系統(tǒng)的可觀性。增加觀測次數(shù)也可以提高所提方法的IROD的性能,但并不是次數(shù)越多性能越好,在觀測次數(shù)大于一定數(shù)量后,算法性能趨于穩(wěn)定。因此,在實際工程應用中,應該合理選擇采樣周期和觀測次數(shù)。

(2)提出的二分法和牛頓迭代法可以快速地獲得僅測角相對導航IROD的解。但是,所提出的方法在主從隨動軌道場景中的IROD性能明顯低于其他3種軌道場景,這是由于主從隨動軌道場景中的系統(tǒng)可觀性較差所導致的,所以未來將繼續(xù)開展僅測角相對導航可觀性最優(yōu)的軌跡設計的研究。