寇 鵬, 劉永祥, 張 弛, 李瑋杰, 張雙輝, 霍 凱
(1. 國防科技大學(xué)電子科學(xué)學(xué)院, 湖南 長沙 410073; 2. 西安衛(wèi)星測控中心, 陜西 西安 710600)
航天器姿態(tài)是描述其在軌狀態(tài)的主要參數(shù),更是判斷其動作意圖以及故障搶救的重要依據(jù)。航天器姿態(tài)通常定義為軌道坐標(biāo)系下的偏航角、俯仰角和滾動角。一般而言,大部分正常工作的航天器為三軸或自旋穩(wěn)定姿態(tài),而失效的航天器多為重力梯度穩(wěn)定或自由翻滾等姿態(tài)[1-3]。逆合成孔徑雷達(dá)(inverse synthetic aperture radar, ISAR)具有全天時、全天候工作的優(yōu)點(diǎn),是航天器姿態(tài)估計的重要手段[4-6]。
航天器ISAR像可以看作是其三維結(jié)構(gòu)在成像平面上的投影。ISAR成像不同于光學(xué)成像,其成像平面與雷達(dá)視線(line of sight, LOS)和目標(biāo)位置緊密相關(guān)。文獻(xiàn)[7-8]認(rèn)為成像平面為相干積累時間(coherent processing interval, CPI)內(nèi)LOS掃過的平面,這種假設(shè)忽略了LOS通常為非共面矢量,因而其僅適用于雷達(dá)站在航天器軌道面內(nèi)或CPI非常短的特殊情況。通常而言,姿態(tài)穩(wěn)定航天器為了提高方位向的分辨率,往往需要較長CPI,LOS的非共面性一般不能忽略。目前,航天器姿態(tài)估計的方法主要有兩種。第一種是基于三維重建的方法[9-12],這種方法通過散射點(diǎn)相關(guān)性的因子分解方法重構(gòu)目標(biāo)的幾何結(jié)構(gòu),從而估計目標(biāo)的姿態(tài)。這種方法本質(zhì)是獲得不同姿態(tài)下散射中心的投影矩陣關(guān)聯(lián),缺點(diǎn)是對視角變化較為敏感,容易受“角閃爍”和噪聲等影響。第二種是典型結(jié)構(gòu)模板匹配法[13-15],這種方法通常需要提取特定的幾何結(jié)構(gòu)(如矩形、線性結(jié)構(gòu)等)進(jìn)行匹配,在觀測積累數(shù)據(jù)樣本充分的情況下,這種方法具有較高的估計精度[16-17]。但現(xiàn)有方法對幾何結(jié)構(gòu)的提取常需要人工介入[18-19],且對復(fù)雜結(jié)構(gòu)航天器姿態(tài)估計的最優(yōu)幾何結(jié)構(gòu)選取問題,目前尚未見到相關(guān)文獻(xiàn)報道。
本文提出了對成像CPI內(nèi)LOS矢量進(jìn)行奇異值分解(singular value decomposition, SVD)[20],從而得到最優(yōu)ISAR成像平面的方法。在此基礎(chǔ)上,提出了利用Unet3+深度網(wǎng)絡(luò)[21]分割出航天器典型部件,進(jìn)而自動提取典型部件的線性結(jié)構(gòu)的方法。針對實(shí)測ISAR像存在散焦、腔體散射和噪聲等造成部件尺寸提取誤差較大的問題,提出了使用不同部件結(jié)構(gòu)建立優(yōu)化函數(shù),通過求解不同結(jié)構(gòu)的優(yōu)化函數(shù)從而得到復(fù)雜結(jié)構(gòu)航天器在軌姿態(tài)角的方法。
如圖1所示,建立J2000.0慣性坐標(biāo)系O-XeYeZe,坐標(biāo)原點(diǎn)為地球質(zhì)心O,Xe軸在基本平面內(nèi)由地球質(zhì)心指向J2000.0慣性坐標(biāo)系的平春分點(diǎn),Ze指向北極方向,Ye軸與Xe軸、Ze軸構(gòu)成右手坐標(biāo)系。建立軌道坐標(biāo)系Ot-XoYoZo,坐標(biāo)原點(diǎn)為航天器質(zhì)心Ot,Xo軸指向航天器運(yùn)動方向,Yo軸指向軌道面的負(fù)法向,Zo軸指向地心。建立本體坐標(biāo)系Ot-XsYsZs,坐標(biāo)原點(diǎn)為航天器質(zhì)心Ot,Xs軸指向航天器主體最長軸方向;Ys軸指向航天器帆板方向,Zo軸與Xs軸、Ys軸構(gòu)成右手系。
圖1 在軌航天器坐標(biāo)系Fig.1 In orbit spacecraft coordinate system
如圖2所示,航天器在軌姿態(tài)定義為本體坐標(biāo)系與軌道坐標(biāo)系的夾角,即偏航角α、俯仰角β和滾動角γ。姿態(tài)角可以理解為軌道坐標(biāo)系按Z-Y-X順序旋轉(zhuǎn)α、β和γ到本體坐標(biāo)系。一般而言,姿態(tài)穩(wěn)定航天器的α、β和γ保持不變。
圖2 在軌航天器姿態(tài)角定義Fig.2 Definition of attitude angle of spacecraft in orbit
如圖3所示,當(dāng)雷達(dá)站不在航天器軌道面內(nèi)時,在一個CPI內(nèi),各時刻雷達(dá)LOS不在同一平面內(nèi),通常不能簡單地認(rèn)為雷達(dá)LOS掃過的即為成像平面,因此這里提出通過SVD分解尋找最優(yōu)等效成像平面的方法。
圖3 姿態(tài)穩(wěn)定航天器ISAR等效成像平面Fig.3 ISAR equivalent imaging plane of attitude stabilized spacecraft
若ti(i=0,1,…,N)時刻,雷達(dá)在J2000.0慣性坐標(biāo)系下位置矢量為ORti,航天器在J2000.0慣性坐標(biāo)系下位置矢量為OTti,則雷達(dá)LOS矢量RTti為
RTti=ORti-OTti
(1)
定義CPI內(nèi)LOS矢量矩陣LOST為
LOST=[RTt0,RTt1,…,RTtN][RTt0,RTt1,…,RTtN]T
(2)
利用SVD分解將LOST進(jìn)行分解:
LOST=UΔVT
(3)
式中:Δ∈R3×3為對角矩陣,其對角元素為LOST的奇異值;U∈R3×3、V∈R3×3均為列正交矩陣。記Δ最小奇異值對應(yīng)的行號為minl,U的第minl行向量為umin l,則等效成像平面的距離向矢量aR可以表示為
(4)
等效成像平面的距離向矢量aC可表示為
aC=aR×umin l
(5)
等效成像平面的法向量aL可以表示為
aL=aC×aR
(6)
Unet3+中的每一個解碼器都結(jié)合了全部編碼器的特征,這些不同尺度的特征能夠獲取細(xì)粒度的細(xì)節(jié)和粗粒度的語義,能有效彌補(bǔ)Unet[22]和Unet++[23-24]不能精確分割圖像中部件位置和邊界的缺點(diǎn),因此本文采用Unet3+網(wǎng)絡(luò)進(jìn)行部件分割,網(wǎng)絡(luò)結(jié)構(gòu)如圖4所示[21]。網(wǎng)絡(luò)損失函數(shù)采用Focal LOSS[25]與IOU(intersection over union)LOSS[26]的加權(quán)求和,其中Focal LOSS主要是為了解決目標(biāo)檢測中正負(fù)樣本比例嚴(yán)重失衡的問題。IOU損失表示預(yù)測框A和真實(shí)框B之間交并比的差值,反映預(yù)測框的檢測效果。
圖4 Unet3+網(wǎng)絡(luò)模型Fig.4 Unet3+network model
(7)
(8)
若航天器等效線性結(jié)構(gòu)矢量為
M=Pi-Pj(i≠k)
則M的ISAR像TM為
(9)
可見TM中不包含圖像中心[u0,v0]T,即消除了序列ISAR圖像中心偏移對姿態(tài)估計的偏差。利用三維姿態(tài)在成像平面的投影誤差建立優(yōu)化函數(shù),通過經(jīng)典粒子群優(yōu)化(particle swarm optimization, PSO)算法[28]等優(yōu)化方法可以最終實(shí)現(xiàn)在軌姿態(tài)的優(yōu)化求解估計。若航天器第j個部件的線性結(jié)構(gòu)矢量為Μj=[m1j,m2j,m3j],Mj在tj時刻的成像平面投影可以寫為
(10)
若總觀測時長tN時間內(nèi)共得到N幅ISAR像,姿態(tài)估計問題可以轉(zhuǎn)化為求解代價函數(shù)f(α,β,γ)的最小值:
(11)
利用PSO即可求解該優(yōu)化函數(shù)。事實(shí)上,圖5為德國TIRA雷達(dá)公布的實(shí)測航天器ISAR像[29]。由于存在散焦、腔體散射和噪聲等因素的影響,航天器ISAR像邊界輪廓往往確定困難,進(jìn)而提取的部件尺寸LRj和LCj往往誤差較大,而部件在ISAR圖像中的矢量方向φj不受圖像邊界輪廓影響,因此實(shí)測序列ISAR圖像可采用部件線性結(jié)構(gòu)矢量或僅使用矢量方向來估計姿態(tài)。
圖5 德國TIRA雷達(dá)實(shí)測航天器ISAR像Fig.5 ISAR images of spacecraft measured by Germany TIRA radar
本文算法流程如圖6所示。首先根據(jù)雷達(dá)回波數(shù)據(jù),采用傳統(tǒng)RD算法進(jìn)行ISAR序列成像;然后對得到的ISAR圖像序列進(jìn)行部件分割和線性結(jié)構(gòu)提取,獲取其在RD成像平面內(nèi)的長度、方向等觀測信息;接著通過航天器軌道及雷達(dá)站位置信息計算各幀圖像對應(yīng)的成像平面,利用成像投影幾何關(guān)系進(jìn)行圖像線性結(jié)構(gòu)的提取匹配估計姿態(tài)角。
圖6 本文算法流程圖Fig.6 Flowchart of the proposed algorithm
在仿真實(shí)驗中,設(shè)置雷達(dá)站址為北緯107.616°,東經(jīng)33.126°,高程為10 m。雷達(dá)中心頻率為10 GHz,帶寬為1 GHz,脈沖重復(fù)頻率為100 Hz,脈寬為5 μs。第5.1節(jié)驗證了本文所提等效成像平面的準(zhǔn)確性;第5.2節(jié)分析了不同軌道航天器ISAR成像CPI時間內(nèi)LOS非共面性;第5.3節(jié)分析了不同線性結(jié)構(gòu)對航天器在軌姿態(tài)估計誤差的影響。
圖7為利用42個點(diǎn)模型模擬包含圓柱主體和太陽能帆板的航天器。航天器軌道高度為375 km,傾角為43.012°,繞地周期為92.193 min,仿真時長為500 s。首先,根據(jù)雷達(dá)站對點(diǎn)模型回波進(jìn)行仿真,利用距離多普勒(range Doppler, RD)算法得到序列ISAR圖像,典型結(jié)果如圖8(a)和圖8(d)所示。圖8(a)中CPI內(nèi)LOS在軌道坐標(biāo)系內(nèi)的變化曲線如圖9所示,可見一般情況下LOS為非共面矢量。圖8(b)和圖8(e)使用Otsu算法[30-31]提取序列ISAR圖像的散射點(diǎn)位置。根據(jù)推導(dǎo)的等效成像平面在J2000.0慣性坐標(biāo)系中的位置,將目標(biāo)的三維坐標(biāo)投影到成像平面上,得到序列“投影圖像”如圖8(c)和圖8(f)所示。最后,逐點(diǎn)比較序列ISAR圖像散射點(diǎn)和三維點(diǎn)模型在不同成像平面“投影圖像”的位置誤差,驗證本文所推導(dǎo)的等效成像平面的準(zhǔn)確性。
圖7 點(diǎn)模型三維分布Fig.7 Three-dimensional distribution of point model
圖8 點(diǎn)模型的ISAR像、散射中心及投影像Fig.8 ISAR image, scattering center and projection image of point model
如圖8(c)、圖8(f)所示,以CPI起止幀LOS構(gòu)成的成像平面與圖8(b)、圖8(e)所示散射中心逐點(diǎn)匹配后,求得平均歐式距離的位置誤差分別為26.5 cm和10.4 cm,而本文提出的等效成像平面平均歐式距離的位置誤差分別為7.7 cm和3.6 cm,說明本文提出的等效成像平面更加準(zhǔn)確、合理。
實(shí)驗首先分析了在相同CPI時間內(nèi),LOS非共面性;然后采用點(diǎn)模型驗證ISAR等效成像平面的準(zhǔn)確性。選取近圓形軌道高度分別為300 km、500 km和700 km,軌道傾角為55°、98°和120°的共9顆航天器進(jìn)行分析,利用STK(sa-tellite tool kit)仿真軟件及SGP4(simplified general perturbation 4)模型計算航天器位置,利用雷達(dá)站的經(jīng)緯度和高度計算LOS序列。成像CPI定義為成像累積時間取方位像分辨率與距離像分辨率相等的時間。在第1.2節(jié)圖3中,以本文SVD分解法構(gòu)成成像平面A,以RTt0和RTtN構(gòu)成平面B,表1給出了246圈次的平面A與平面B法向量夾角的平均值。
表1 平面A與平面B法向量夾角的平均值
從北京時間2021年5月2日4時0分0秒至2021年5月5日4時0分0秒,雷達(dá)對9顆航天器共可見246圈。圖10給出了不同軌道參數(shù)的平面A與平面B法向量的最大夾角,可見圈次雷達(dá)最大俯仰角變化情況,可見CPI時間內(nèi)LOS是非共面角度,隨軌道高度降低而增加,且與雷達(dá)最高俯仰角基本是二次拋物線關(guān)系。
圖10 不同軌道參數(shù)平面A與平面B法向量的最大夾角Fig.10 Maximum angle between normal vectors of plane A and B with different orbital parameters
如圖11所示,截至2021年6月19日,中國空間站在軌模型為“王”字形結(jié)構(gòu),模型可以簡化成主體、帆板1、帆板2和帆板3四個線性結(jié)構(gòu)。實(shí)驗首先使用3dmax軟件建立中國空間站三維模型,根據(jù)中國空間站TLE(two-line element)根數(shù)計算航天器在J2000.0慣性坐標(biāo)系下的位置和速度,然后用物理光學(xué)方法計算模型雷達(dá)散射截面(radar cross section, RCS),進(jìn)而仿真雷達(dá)回波數(shù)據(jù)。根據(jù)雷達(dá)回波,利用RD算法成像后得到序列ISAR像,典型ISAR像結(jié)果如圖12(a)、圖12(g)和圖12(m)所示。實(shí)驗仿真了5個圈次共163幅ISAR像,選取4個圈次共117幅ISAR像作為訓(xùn)練集,選取一個圈次共46幅ISAR像作為測試集。
(1) 部件分割與線性結(jié)構(gòu)提取實(shí)驗
部件分割典型結(jié)果如圖12(b)~圖12(e)、圖12(h)~圖12(k)和圖12(n)~圖12(q)所示,同時利用Radon變換檢測的直線與部件分割結(jié)果輪廓最遠(yuǎn)交點(diǎn)可以計算部件的ISAR圖像尺寸,根據(jù)徑向尺寸與橫向尺寸比值可以計算得到ISAR圖像的部件方向,計算結(jié)果如表2所示。整個測試集部件的徑向尺寸、橫向尺寸和線性結(jié)構(gòu)方向誤差如表3所示,可見部件線性結(jié)構(gòu)的尺寸誤差相對較大,而線性結(jié)構(gòu)方向尺寸誤差相對較小。這其中的原因主要是ISAR圖像邊緣輪廓較難確定,且散焦、高旁瓣和斑點(diǎn)噪聲等因素造成尺寸估計誤差較大,而多數(shù)散射點(diǎn)相對位置關(guān)系分布則較為穩(wěn)定,所以部件線性結(jié)構(gòu)的方向誤差相對較小。
表2 中國空間站典型ISAR像部件尺寸與方向
表3 測試圈部件線性結(jié)構(gòu)矢量平均誤差
(2) 不同線性結(jié)構(gòu)姿態(tài)估計誤差實(shí)驗
實(shí)驗對測試圈次46幅ISAR像,分別選取“天和核心艙”的“十”字型、“天和天舟”組合體的“士”字型和“天和天舟神舟十二”組合體的“王”字型3種線性結(jié)構(gòu)組合,不同線性結(jié)構(gòu)矢量對姿態(tài)估計的精度如圖13所示。
圖13 3種線性結(jié)構(gòu)矢量姿態(tài)估計誤差曲線Fig.13 Error curves of three kinds of linear structure vector attitude estimation
定義平均姿態(tài)角誤差φ=(α+β+γ)/3,當(dāng)ISAR像數(shù)量為5幅、25幅和45幅時,φ的取值見表4所示。當(dāng)ISAR像的數(shù)量相同時,“王”字形的φ均值最小,“十”字形的φ均值最大,說明在估計復(fù)雜結(jié)構(gòu)航天器姿態(tài)角時,應(yīng)該充分選取盡量多的線性結(jié)構(gòu)以減小估計誤差。同時,隨著使用ISAR像數(shù)量增多,φ減小,當(dāng)ISAR像為45幅時,3種線性結(jié)構(gòu)矢量的平均姿態(tài)角誤差φ基本相同,均在1°以內(nèi)。
表4 3種線性結(jié)構(gòu)矢量的平均姿態(tài)角誤差φ
實(shí)測ISAR像部件線性結(jié)構(gòu)矢量尺寸提取往往存在較大誤差,但部件線性結(jié)構(gòu)矢量方向通常誤差較小,因此僅使用3種線性結(jié)構(gòu)的方向來估計姿態(tài)。當(dāng)ISAR像數(shù)量為5幅、15幅和45幅時,僅使用3種線性結(jié)構(gòu)方向的平均姿態(tài)角誤差φ如表5所示。如圖14所示,當(dāng)ISAR像的數(shù)量多于15幅時,3種結(jié)構(gòu)的姿態(tài)估計誤差幾乎相等。同時,當(dāng)ISAR像數(shù)量增多至45幅時,平均姿態(tài)角φ誤差同樣可以小于2°。實(shí)驗說明在ISAR像數(shù)量較多時,僅使用線性結(jié)構(gòu)矢量的方向估計姿態(tài)仍然具有較高的精度。
表5 3種線性結(jié)構(gòu)方向的平均姿態(tài)角誤差φ
圖14 3種線性結(jié)構(gòu)方向姿態(tài)估計誤差曲線Fig.14 Error curves of orientation and attitude estimation of three linear structures
本文從ISAR成像原理出發(fā),分析航天器軌道參數(shù)與ISAR成像平面之間的聯(lián)系,通過雷達(dá)LOS矢量SVD分解得到姿態(tài)穩(wěn)定航天器的等效ISAR成像平面。然后,通過深度神經(jīng)網(wǎng)絡(luò)及圖像處理方法自動獲取復(fù)雜結(jié)構(gòu)航天器線性部件線性結(jié)構(gòu)矢量。最后,分析了不同線性結(jié)構(gòu)對姿態(tài)估計誤差的影響。實(shí)驗表明,等效線性結(jié)構(gòu)越多,姿態(tài)估計誤差越小;使用的ISAR圖像數(shù)量越多,姿態(tài)估計誤差越小。同時,僅采用部件線性結(jié)構(gòu)的方向可有效彌補(bǔ)實(shí)測ISAR像中部件尺寸估計誤差較大的缺點(diǎn),進(jìn)而提升了方法的魯棒性。