寇 鵬, 劉永祥, 張 弛, 李瑋杰, 張雙輝, 霍 凱

(1. 國防科技大學電子科學學院, 湖南 長沙 410073; 2. 西安衛(wèi)星測控中心, 陜西 西安 710600)

0 引 言

航天器姿態(tài)是描述其在軌狀態(tài)的主要參數(shù),更是判斷其動作意圖以及故障搶救的重要依據(jù)。航天器姿態(tài)通常定義為軌道坐標系下的偏航角、俯仰角和滾動角。一般而言,大部分正常工作的航天器為三軸或自旋穩(wěn)定姿態(tài),而失效的航天器多為重力梯度穩(wěn)定或自由翻滾等姿態(tài)[1-3]。逆合成孔徑雷達(inverse synthetic aperture radar, ISAR)具有全天時、全天候工作的優(yōu)點,是航天器姿態(tài)估計的重要手段[4-6]。

航天器ISAR像可以看作是其三維結(jié)構(gòu)在成像平面上的投影。ISAR成像不同于光學成像,其成像平面與雷達視線(line of sight, LOS)和目標位置緊密相關(guān)。文獻[7-8]認為成像平面為相干積累時間(coherent processing interval, CPI)內(nèi)LOS掃過的平面,這種假設(shè)忽略了LOS通常為非共面矢量,因而其僅適用于雷達站在航天器軌道面內(nèi)或CPI非常短的特殊情況。通常而言,姿態(tài)穩(wěn)定航天器為了提高方位向的分辨率,往往需要較長CPI,LOS的非共面性一般不能忽略。目前,航天器姿態(tài)估計的方法主要有兩種。第一種是基于三維重建的方法[9-12],這種方法通過散射點相關(guān)性的因子分解方法重構(gòu)目標的幾何結(jié)構(gòu),從而估計目標的姿態(tài)。這種方法本質(zhì)是獲得不同姿態(tài)下散射中心的投影矩陣關(guān)聯(lián),缺點是對視角變化較為敏感,容易受“角閃爍”和噪聲等影響。第二種是典型結(jié)構(gòu)模板匹配法[13-15],這種方法通常需要提取特定的幾何結(jié)構(gòu)(如矩形、線性結(jié)構(gòu)等)進行匹配,在觀測積累數(shù)據(jù)樣本充分的情況下,這種方法具有較高的估計精度[16-17]。但現(xiàn)有方法對幾何結(jié)構(gòu)的提取常需要人工介入[18-19],且對復雜結(jié)構(gòu)航天器姿態(tài)估計的最優(yōu)幾何結(jié)構(gòu)選取問題,目前尚未見到相關(guān)文獻報道。

本文提出了對成像CPI內(nèi)LOS矢量進行奇異值分解(singular value decomposition, SVD)[20],從而得到最優(yōu)ISAR成像平面的方法。在此基礎(chǔ)上,提出了利用Unet3+深度網(wǎng)絡[21]分割出航天器典型部件,進而自動提取典型部件的線性結(jié)構(gòu)的方法。針對實測ISAR像存在散焦、腔體散射和噪聲等造成部件尺寸提取誤差較大的問題,提出了使用不同部件結(jié)構(gòu)建立優(yōu)化函數(shù),通過求解不同結(jié)構(gòu)的優(yōu)化函數(shù)從而得到復雜結(jié)構(gòu)航天器在軌姿態(tài)角的方法。

1 在軌姿態(tài)穩(wěn)定航天器ISAR成像平面分析

1.1 在軌航天器坐標系與姿態(tài)角定義

如圖1所示,建立J2000.0慣性坐標系O-XeYeZe,坐標原點為地球質(zhì)心O,Xe軸在基本平面內(nèi)由地球質(zhì)心指向J2000.0慣性坐標系的平春分點,Ze指向北極方向,Ye軸與Xe軸、Ze軸構(gòu)成右手坐標系。建立軌道坐標系Ot-XoYoZo,坐標原點為航天器質(zhì)心Ot,Xo軸指向航天器運動方向,Yo軸指向軌道面的負法向,Zo軸指向地心。建立本體坐標系Ot-XsYsZs,坐標原點為航天器質(zhì)心Ot,Xs軸指向航天器主體最長軸方向;Ys軸指向航天器帆板方向,Zo軸與Xs軸、Ys軸構(gòu)成右手系。

圖1 在軌航天器坐標系Fig.1 In orbit spacecraft coordinate system

如圖2所示,航天器在軌姿態(tài)定義為本體坐標系與軌道坐標系的夾角,即偏航角α、俯仰角β和滾動角γ。姿態(tài)角可以理解為軌道坐標系按Z-Y-X順序旋轉(zhuǎn)α、β和γ到本體坐標系。一般而言,姿態(tài)穩(wěn)定航天器的α、β和γ保持不變。

圖2 在軌航天器姿態(tài)角定義Fig.2 Definition of attitude angle of spacecraft in orbit

1.2 姿態(tài)穩(wěn)定航天器ISAR等效成像平面

如圖3所示,當雷達站不在航天器軌道面內(nèi)時,在一個CPI內(nèi),各時刻雷達LOS不在同一平面內(nèi),通常不能簡單地認為雷達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)時刻,雷達在J2000.0慣性坐標系下位置矢量為ORti,航天器在J2000.0慣性坐標系下位置矢量為OTti,則雷達LOS矢量RTti為

RTti=ORti-OTti

(1)

定義CPI內(nèi)LOS矢量矩陣LOST為

LOST=[RTt0,RTt1,…,RTtN][RTt0,RTt1,…,RTtN]T

(2)

利用SVD分解將LOST進行分解:

LOST=UΔVT

(3)

式中:Δ∈R3×3為對角矩陣,其對角元素為LOST的奇異值;U∈R3×3、V∈R3×3均為列正交矩陣。記Δ最小奇異值對應的行號為minl,U的第minl行向量為umin l,則等效成像平面的距離向矢量aR可以表示為

(4)

等效成像平面的距離向矢量aC可表示為

aC=aR×umin l

(5)

等效成像平面的法向量aL可以表示為

aL=aC×aR

(6)

2 航天器部件分割與線性結(jié)構(gòu)提取

2.1 基于Unet3+神經(jīng)網(wǎng)絡的航天器部件分割

Unet3+中的每一個解碼器都結(jié)合了全部編碼器的特征,這些不同尺度的特征能夠獲取細粒度的細節(jié)和粗粒度的語義,能有效彌補Unet[22]和Unet++[23-24]不能精確分割圖像中部件位置和邊界的缺點,因此本文采用Unet3+網(wǎng)絡進行部件分割,網(wǎng)絡結(jié)構(gòu)如圖4所示[21]。網(wǎng)絡損失函數(shù)采用Focal LOSS[25]與IOU(intersection over union)LOSS[26]的加權(quán)求和,其中Focal LOSS主要是為了解決目標檢測中正負樣本比例嚴重失衡的問題。IOU損失表示預測框A和真實框B之間交并比的差值,反映預測框的檢測效果。

圖4 Unet3+網(wǎng)絡模型Fig.4 Unet3+network model

2.2 部件線性結(jié)構(gòu)提取

3 在軌航天器ISAR成像投影關(guān)系與姿態(tài)求解

(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)化方法可以最終實現(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ù)。事實上,圖5為德國TIRA雷達公布的實測航天器ISAR像[29]。由于存在散焦、腔體散射和噪聲等因素的影響,航天器ISAR像邊界輪廓往往確定困難,進而提取的部件尺寸LRj和LCj往往誤差較大,而部件在ISAR圖像中的矢量方向φj不受圖像邊界輪廓影響,因此實測序列ISAR圖像可采用部件線性結(jié)構(gòu)矢量或僅使用矢量方向來估計姿態(tài)。

圖5 德國TIRA雷達實測航天器ISAR像Fig.5 ISAR images of spacecraft measured by Germany TIRA radar

4 本文算法流程

本文算法流程如圖6所示。首先根據(jù)雷達回波數(shù)據(jù),采用傳統(tǒng)RD算法進行ISAR序列成像;然后對得到的ISAR圖像序列進行部件分割和線性結(jié)構(gòu)提取,獲取其在RD成像平面內(nèi)的長度、方向等觀測信息;接著通過航天器軌道及雷達站位置信息計算各幀圖像對應的成像平面,利用成像投影幾何關(guān)系進行圖像線性結(jié)構(gòu)的提取匹配估計姿態(tài)角。

圖6 本文算法流程圖Fig.6 Flowchart of the proposed algorithm

5 仿真實驗

在仿真實驗中,設(shè)置雷達站址為北緯107.616°,東經(jīng)33.126°,高程為10 m。雷達中心頻率為10 GHz,帶寬為1 GHz,脈沖重復頻率為100 Hz,脈寬為5 μs。第5.1節(jié)驗證了本文所提等效成像平面的準確性;第5.2節(jié)分析了不同軌道航天器ISAR成像CPI時間內(nèi)LOS非共面性;第5.3節(jié)分析了不同線性結(jié)構(gòu)對航天器在軌姿態(tài)估計誤差的影響。

5.1 航天器ISAR等效成像平面準確性驗證

圖7為利用42個點模型模擬包含圓柱主體和太陽能帆板的航天器。航天器軌道高度為375 km,傾角為43.012°,繞地周期為92.193 min,仿真時長為500 s。首先,根據(jù)雷達站對點模型回波進行仿真,利用距離多普勒(range Doppler, RD)算法得到序列ISAR圖像,典型結(jié)果如圖8(a)和圖8(d)所示。圖8(a)中CPI內(nèi)LOS在軌道坐標系內(nèi)的變化曲線如圖9所示,可見一般情況下LOS為非共面矢量。圖8(b)和圖8(e)使用Otsu算法[30-31]提取序列ISAR圖像的散射點位置。根據(jù)推導的等效成像平面在J2000.0慣性坐標系中的位置,將目標的三維坐標投影到成像平面上,得到序列“投影圖像”如圖8(c)和圖8(f)所示。最后,逐點比較序列ISAR圖像散射點和三維點模型在不同成像平面“投影圖像”的位置誤差,驗證本文所推導的等效成像平面的準確性。

圖7 點模型三維分布Fig.7 Three-dimensional distribution of point model

圖8 點模型的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)所示散射中心逐點匹配后,求得平均歐式距離的位置誤差分別為26.5 cm和10.4 cm,而本文提出的等效成像平面平均歐式距離的位置誤差分別為7.7 cm和3.6 cm,說明本文提出的等效成像平面更加準確、合理。

5.2 ISAR成像CPI時間內(nèi)LOS非共面性分析

實驗首先分析了在相同CPI時間內(nèi),LOS非共面性;然后采用點模型驗證ISAR等效成像平面的準確性。選取近圓形軌道高度分別為300 km、500 km和700 km,軌道傾角為55°、98°和120°的共9顆航天器進行分析,利用STK(sa-tellite tool kit)仿真軟件及SGP4(simplified general perturbation 4)模型計算航天器位置,利用雷達站的經(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秒,雷達對9顆航天器共可見246圈。圖10給出了不同軌道參數(shù)的平面A與平面B法向量的最大夾角,可見圈次雷達最大俯仰角變化情況,可見CPI時間內(nèi)LOS是非共面角度,隨軌道高度降低而增加,且與雷達最高俯仰角基本是二次拋物線關(guān)系。

圖10 不同軌道參數(shù)平面A與平面B法向量的最大夾角Fig.10 Maximum angle between normal vectors of plane A and B with different orbital parameters

5.3 線性結(jié)構(gòu)選取對“中國空間站”姿態(tài)估計精度影響分析

如圖11所示,截至2021年6月19日,中國空間站在軌模型為“王”字形結(jié)構(gòu),模型可以簡化成主體、帆板1、帆板2和帆板3四個線性結(jié)構(gòu)。實驗首先使用3dmax軟件建立中國空間站三維模型,根據(jù)中國空間站TLE(two-line element)根數(shù)計算航天器在J2000.0慣性坐標系下的位置和速度,然后用物理光學方法計算模型雷達散射截面(radar cross section, RCS),進而仿真雷達回波數(shù)據(jù)。根據(jù)雷達回波,利用RD算法成像后得到序列ISAR像,典型ISAR像結(jié)果如圖12(a)、圖12(g)和圖12(m)所示。實驗仿真了5個圈次共163幅ISAR像,選取4個圈次共117幅ISAR像作為訓練集,選取一個圈次共46幅ISAR像作為測試集。

(1) 部件分割與線性結(jié)構(gòu)提取實驗

部件分割典型結(jié)果如圖12(b)～圖12(e)、圖12(h)～圖12(k)和圖12(n)～圖12(q)所示,同時利用Radon變換檢測的直線與部件分割結(jié)果輪廓最遠交點可以計算部件的ISAR圖像尺寸,根據(jù)徑向尺寸與橫向尺寸比值可以計算得到ISAR圖像的部件方向,計算結(jié)果如表2所示。整個測試集部件的徑向尺寸、橫向尺寸和線性結(jié)構(gòu)方向誤差如表3所示,可見部件線性結(jié)構(gòu)的尺寸誤差相對較大,而線性結(jié)構(gòu)方向尺寸誤差相對較小。這其中的原因主要是ISAR圖像邊緣輪廓較難確定,且散焦、高旁瓣和斑點噪聲等因素造成尺寸估計誤差較大,而多數(shù)散射點相對位置關(guān)系分布則較為穩(wěn)定,所以部件線性結(jié)構(gòu)的方向誤差相對較小。

表2 中國空間站典型ISAR像部件尺寸與方向

表3 測試圈部件線性結(jié)構(gòu)矢量平均誤差

(2) 不同線性結(jié)構(gòu)姿態(tài)估計誤差實驗

實驗對測試圈次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,當ISAR像數(shù)量為5幅、25幅和45幅時,φ的取值見表4所示。當ISAR像的數(shù)量相同時,“王”字形的φ均值最小,“十”字形的φ均值最大,說明在估計復雜結(jié)構(gòu)航天器姿態(tài)角時,應該充分選取盡量多的線性結(jié)構(gòu)以減小估計誤差。同時,隨著使用ISAR像數(shù)量增多,φ減小,當ISAR像為45幅時,3種線性結(jié)構(gòu)矢量的平均姿態(tài)角誤差φ基本相同,均在1°以內(nèi)。

表4 3種線性結(jié)構(gòu)矢量的平均姿態(tài)角誤差φ

實測ISAR像部件線性結(jié)構(gòu)矢量尺寸提取往往存在較大誤差,但部件線性結(jié)構(gòu)矢量方向通常誤差較小,因此僅使用3種線性結(jié)構(gòu)的方向來估計姿態(tài)。當ISAR像數(shù)量為5幅、15幅和45幅時,僅使用3種線性結(jié)構(gòu)方向的平均姿態(tài)角誤差φ如表5所示。如圖14所示,當ISAR像的數(shù)量多于15幅時,3種結(jié)構(gòu)的姿態(tài)估計誤差幾乎相等。同時,當ISAR像數(shù)量增多至45幅時,平均姿態(tài)角φ誤差同樣可以小于2°。實驗說明在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

6 結(jié) 論

本文從ISAR成像原理出發(fā),分析航天器軌道參數(shù)與ISAR成像平面之間的聯(lián)系,通過雷達LOS矢量SVD分解得到姿態(tài)穩(wěn)定航天器的等效ISAR成像平面。然后,通過深度神經(jīng)網(wǎng)絡及圖像處理方法自動獲取復雜結(jié)構(gòu)航天器線性部件線性結(jié)構(gòu)矢量。最后,分析了不同線性結(jié)構(gòu)對姿態(tài)估計誤差的影響。實驗表明,等效線性結(jié)構(gòu)越多,姿態(tài)估計誤差越小;使用的ISAR圖像數(shù)量越多,姿態(tài)估計誤差越小。同時,僅采用部件線性結(jié)構(gòu)的方向可有效彌補實測ISAR像中部件尺寸估計誤差較大的缺點,進而提升了方法的魯棒性。