地基光學圖像中本體-帆板結構航天器姿態估計方法

楊 成, 李 勰, 李志輝, 孫 軍

(1. 北京航天飛行控制中心，北京 100094；2. 航天飛行動力學技術重點實驗室，北京 100094； 3.中國空氣動力研究與發展中心，綿陽 621000)

0 引 言

三維姿態信息是反映空間目標運動狀態的重要參數，對空間態勢感知、目標運動分析、航天器運行控制和故障診斷具有重要意義?？臻g目標的觀測主要采用地基雷達觀測和光學觀測兩種方式[1-3]。雷達觀測向目標發射電磁波并接收返回信號，屬于主動探測，其探測能力與作用距離的四次方成反比，主要作用于近地軌道空間。光學觀測屬于被動觀測，探測能力與距離的平方成反比，具有更好的跟蹤能力。隨著光學成像設備性能的發展，地基光學圖像在空間目標姿態估計的應用越來越多[4-5]。

從圖像信息中提取空間目標姿態主要有兩類方法。一類是模板匹配方法，根據空間目標的幾何外形、表面材質和空間構造等三維模型信息，預先仿真不同姿態下的觀測圖像，建立目標觀測圖像庫，然后將實際觀測圖像進行圖像匹配，根據匹配相似度確定目標的姿態。這類方法需要提前構建匹配圖像庫，在姿態估計過程中進行大量的匹配處理，具有較高的計算復雜度，實際應用中效率較低，且需要定義合適的圖像匹配相似度，對姿態估計的結果有較大影響[6-8]。另一類是直接根據三維空間到二維圖像的投影關系進行解算，一般是通過提取圖像特征，根據計算機視覺成像原理，建立空間三維信息到二維圖像的幾何投影模型，根據特征在二維圖像的像素坐標反算出三維空間坐標，從而確定目標的三維姿態。這類方法避免了建立圖像庫和匹配搜索處理，但空間目標的特征可能被遮擋，無法實現任意姿態的估計[9-11]。

本文提出一種本體-帆板結構航天器的姿態估計方法。針對該類航天器幾何外形結構特點，定義航天器本體主軸和帆板轉動軸兩軸結構作為圖像特征，提出了地基光學圖像邊緣模糊條件下的特征提取方法，根據三維空間到二維圖像投影方程，建立了航天器姿態解算方法，實現了單站序列圖像姿態和角速度估計。

1 地基光學圖像特征提取

1.1 本體-帆板結構航天器外形特點

航天器包括本體、帆板及天線、敏感器等結構和部件。本體一般為方形或者圓柱形，是航天器的主體，同時為其它儀器設備和有效載荷提供安裝框架。帆板具有較大的尺寸，一般為左右兩部分分別連接在本體兩側，為航天器提供能源[12]。天線、敏感器、發動機等載荷安裝在本體表面，其尺寸相對本體和帆板很小，在進行表面力建模和光學圖像分析等處理時可忽略不計。因此，這類航天器在進行幾何外形建模時可簡化為本體和帆板兩部分組成，稱為本體-帆板結構(Box-Wing)航天器。

圖1 航天器本體-帆板結構示意圖Fig.1 Schematic diagram of Box-Wing spacecraft

1.2 圖像特征定義

地基光學雷達對空間航天器進行觀測，根據視覺原理產生圖像，其成像幾何原理和地面景物成像一致。由于航天器距離地面光學設備的成像距離非常大，成像視場角很小，航天器于空間高速運行，因此圖像的信噪比非常低，椒鹽噪聲影響嚴重。同時，成像過程跨越大氣層，導致成像信息的邊緣非常模糊。這些都給圖像特征的準確提取帶來了影響。如圖2所示，在光學圖像中，可根據航天器形狀先驗知識，辨認出帆板、本體結構，而其它載荷部件的細節丟失。

針對航天器本體-帆板結構特點，定義光學圖像中航天器左、右帆板和本體結構的端點作為圖像特征，獲取三點的像素坐標恢復航天器三維姿態。如圖3所示，P(+X)為本體前端面中心點，P(-Y)為左帆板外測邊緣中心點，P(+Y)為右帆板外測邊緣中心點，P(-Y)和P(+Y)連線中點記為P0。

圖3 圖像特征點定義Fig.3 Definition of feature points

1.3 圖像特征提取

利用數字圖像處理方法，從光學圖像提取特征點的像素坐標[13],步驟如下所述。

1)圖像增強。首先對圖像進行中值濾波，消除圖像中椒鹽噪聲影響；然后進行膨脹，恢復圖像中的連通性。

2)邊緣檢測。對圖像進行Harris邊緣檢測，對邊緣檢測結果使用連通性檢測得到最長的邊緣曲線作為航天器的邊緣。待提取的特征點位于檢測出的邊緣曲線上。

圖5 邊緣檢測結果Fig.5 Results of edge detection

3)特征點確定。根據帆板形狀擬合出帆板邊緣位置，獲得帆板的中心點和帆板轉軸方向，確定P(-Y)和P(+Y)兩特征點的位置；然后查找剩下邊緣點的中點，即為P(+X)特征點的位置。

當圖像中航天器帆板、主體結構完整清晰時，可用上述圖像處理方法自動獲得特征點的像素坐標；若成像質量差或者處于遮擋嚴重的角度圖像中航天器結構不完整，采用人工標注的方法提取特征點。

2 航天器姿態解算

2.1 投影方程解算

在地基光學設備對空間目標跟蹤成像過程中，按照成像模型，定義光學設備坐標系如圖6所示，坐標原點為光心，Y軸指向地心，Z軸沿主光軸指向空間目標方向，X軸與Y軸、Z軸構成右手坐標系。產生的空間目標光學圖像位于X-Y平面，圖像坐標系像素坐標X軸為水平向右，Y軸為垂直向下。

圖6 光學設備坐標系定義Fig.6 Coordinate of optical telescope

對于本體-帆板結構航天器，根據其結構特點，本體系定義為：本體主軸和帆板轉動軸交點為原點，本體主軸向前進方向為X軸，Y軸為轉動軸，Z軸與X軸、Y軸構成右手坐標系。

記地面光學圖像提取到的特征點像素坐標分別為：

(1)

(2)

(3)

帆板中點P0坐標為：

P0=(P(+Y)+P(-Y))/2

(4)

光學圖像像素分辨率k表示一個像素對應空間中的物理長度，帆板中點至本體前端點在光學圖像的二維投影坐標為P(+X)-P0，則其在光學設備坐標系下X-Y坐標為k(P(+X)-P0)。由于航天器本體坐標系的X軸定義為帆板中點至本體前端點方向，則航天器本體坐標系的X軸在光學設備坐標系下坐標可表示為：

(5)

式中：l1為根據航天器的幾何尺寸中，帆板中點P0至本體前端點P(+X)的長度；由于投影平面為X-Y平面，Z軸方向的坐標不能確定，用未知量z1表示。

同理，航天器本體坐標系的Y軸在光學設備坐標系下坐標可表示為：

(6)

式中：l2為根據航天器的幾何尺寸中，帆板中點P0至右帆板外測邊緣中心點P(+Y)的長度；由于投影平面為X-Y平面，Z軸方向的坐標也不能確定，用未知量z2表示。

本方法沒有在光學圖像中提取航天器本體坐標系Z軸的特征，其在光學設備坐標系下坐標用Zc來表示，則航天器本體坐標系到光學設備坐標系的旋轉矩陣為：

(7)

由Xc,Yc為正交單位向量，可得如下方程：

(8)

(9)

(10)

式中:k,z1,z2為未知量，滿足k>0，|z1|≤1，|z2|≤1。

根據式(8)，z1可變換為k的表達式；根據式(9)，z2也可變換為k的表達式。代入式(10)中，得到關于k的方程，可求得k的可能解。根據約束條件，最終可確定k,z1,z2的兩組解，進而求得Xc和Yc，利用Zc=Yc×Xc，即可得到航天器本體坐標系到光學設備坐標系的旋轉矩陣RI=[Xc,Yc,Zc]3×3。

上述方法通過二維圖像信息，確定了三維空間中航天器對應的兩種可能姿態，這兩種姿態關于相機成像平面成鏡像關系，通過單張圖像不能確定唯一姿態。實際應用中，地面光學設備會間隔一定時間持續對航天器進行觀測，通過序列圖像恢復姿態之間的關聯性，對兩種可能結果進行合理性判斷，確定航天器唯一的姿態結果。

2.2 姿態坐標系轉換

1)光學設備坐標系到測站坐標系的轉換

定義測站坐標系原點為光學設備跟蹤天線的旋轉中心，X-Z平面為經過坐標原點的地球切平面，X軸位于該切平面指向東，Z軸指向北，Y軸與X軸、Z軸形成右手坐標系。在對空間目標進行觀測過程中，需要通過轉動觀測設備的方向來瞄準空間目標，產生方位角A和仰角E，如圖7所示。

圖7 測站坐標系定義Fig.7 Coordinate of observer station

方位角A旋轉矩陣為：

(11)

仰角E旋轉矩陣為：

(12)

圖像坐標系到測站坐標系的旋轉矩陣為：

RA×RE

(13)

2)測站坐標系到地球慣性系的轉換

記測站在地球慣性系的位置坐標為P，可計算出測站坐標系Y軸向量，Y軸與數值單位向量[0，0，1]的叉乘，可得到測站坐標系X軸向量，然后Y軸與X軸叉乘得到Z軸向量。于是得到測站坐標系到地球慣性系的旋轉矩陣RC。

圖8 測站坐標系和地球慣性系示意圖Fig.8 Relationship of observer station and J2000 coordinates

3)地球慣性系到軌道坐標系的轉換

根據航天器在地球慣性系的位置和速度矢量，可算出航天器軌道坐標系與地球慣性系的旋轉矩陣RO。

4)歐拉角計算

可算出航天器本體系在地球慣性系的旋轉矩陣為[14]：

MJ2000=RC×RA×RE×RI

(14)

航天器本體系在軌道系的旋轉矩陣為：

MOrbit=RO×RC×RA×RE×RI

(15)

采取“偏航-滾動-俯仰”的順序來進行姿態旋轉，旋轉矩陣M3-1-2和偏航(ψ)-滾動(φ)-俯仰(θ)歐拉角之間的轉換關系為：

(16)

由式(16)，根據航天器旋轉矩陣可反推出航天器偏航-滾動-俯仰角度。

3 試 驗

3.1 實驗設置

利用我國地面光學設備獲取的空間在軌航天器圖像對本文方法進行驗證。航天器為艙體-帆板結構，如圖1所示，艙體為圓柱狀，前端面為錐形，帆板連接于艙體后半部分。本體坐標系X軸沿艙體指向前錐面，Y軸指向右帆板，Z軸與X軸、Y軸形成右手坐標系。航天器帆板展開后長度約30 m，軌道高度約400 km，處于慢速自旋狀態。

地面觀測站位于我國西南地區，如圖9所示，觀測弧段內，航天器沿圖中軌跡左上點往右下方向飛行，光學設備間隔5 s獲得航天器光學圖像。本次觀測過程中，成像時間在凌晨5點左右，成像距離約600 km，成像方位角約-10°～70°，俯仰角約20°，共獲得14張圖像，用于本文姿態估計方法實驗。

圖9 航天器經過觀測站上空的路線Fig.9 The routine of spacecraft over ground station

3.2 實驗步驟

對每張圖像處理過程如下所述。

1)使用第1.3節的方法或者人工方式提取圖像中航天器三個特征點的像素坐標，代入第2.1節方程(8)～(10)中進行投影方程解算，獲得航天器本體坐標系到光學設備坐標系的旋轉矩陣RI。

2)根據圖像對應的觀測設備方位角A和仰角E，計算光學設備坐標系到測站坐標系的旋轉矩陣RA×RE。

3)根據航天器軌道根數和觀測站的經緯度，獲得該幅圖像觀測時刻航天器在地球慣性系的位置和速度，以及觀測站在地球慣性系的位置。按照第2.2節中的方法，得到航天器在地球慣性系和軌道坐標系的旋轉矩陣，并計算相應的姿態歐拉角。

3.3 姿態估計結果

為便于分析，以“偏航-滾動-俯仰”順序的歐拉角來表示姿態。每幅圖像都得到兩個可能的解，對應兩種不同的航天器姿態。對14幅序列圖像進行處理，得到了兩組可能姿態。

第一組姿態結果如表1所示，單位為度(°)。

表1 第一組可能的姿態結果Table 1 Estimation of the first probable attitude

圖10和圖11為第一組姿態角度和擬合情況，實線為根據圖像估計的角度數據，虛線為根據序列圖像估計結果進行一次擬合的直線。

圖11 航天器軌道系姿態結果(第一組)Fig.11 Estimation of attitude in orbit coordination (The first group)

第二組姿態結果如表2所示。

表2 第二組可能的姿態結果Table 2 Estimation of the second probable attitude

圖12和圖13為第二組姿態角度和擬合情況，實線為根據圖像估計的角度數據，虛線為根據序列圖像估計結果進行一次擬合的直線。

圖12 航天器地球慣性系姿態結果(第二組)Fig.12 Estimation of attitude in J2000 coordination (The second group)

圖13 航天器軌道系姿態結果(第二組)Fig.13 Estimation of attitude in orbit coordination (The second group)

本次實驗14幅圖像數據間隔5 s，航天器在此期間處于無控飛行狀態，其姿態變化為勻速旋轉。通過對姿態角度進行線性擬合，得出偏航、滾動、俯仰角度的旋轉角速度。利用擬合的一次曲線作為理論姿態角度，可得出本文方法估計結果的均方誤差如下：

(17)

式(17)中θ′i為本方法估計的姿態角度，θ′i為一次曲線擬合姿態角度，對應均方差反應了本文方法估計姿態角度的一致性。兩組可能姿態的角速度和結果方差如表3所示。

從圖10～13和表3可以看出，第二組姿態結果均方差更小，對航天器慢速旋轉狀態擬合較好，確定第二組姿態為估計結果。圖14為第二組結果對應的航天器姿態情況，在序列圖像成像時刻航天器的位置，繪制航天器本體系坐標軸，可直觀地看出航天器姿態變化情況。

表3 兩組姿態估計結果Table 3 Two probable results of estimated attitude

圖14 估計出的航天器姿態Fig.14 Visualization of the estimated attitude

4 結 論

本文提出了一種從光學圖像估計航天器姿態的方法。針對艙體-帆板結構航天器的幾何外形特點，定義了三個易于提取的特征點，并設計了基于圖像邊緣檢測的自動特征提取算法；建立了特征點二三維投影方程，實現了對三維姿態的求解。本方法無需建立大量的模型匹配庫，避免了復雜的匹配搜索過程，可直接根據特征點獲得姿態結果。對在軌航天器進行了光學圖像姿態估計實驗，姿態角進行線性擬合后的平均估計誤差在3°以內，驗證了本方法的有效性。