航天器在軌運行可視化中姿態數據的處理方法

淡 鵬，王 丹，李 力，李代偉

(1.宇航動力學國家重點實驗室，陜西 西安 710043；2.西安衛星測控中心，陜西 西安 710043)

0 引言

在航天器飛行試驗任務中，地面中心已經越來越多地采用三維可視化技術來表現和展示航天器飛行過程及狀態的變化情況(文獻[1-3]對比進行了研究)。三維可視化技術的使用為工程技術人員提供了更加直觀、更易于理解的飛行過程監視手段，方便了其對任務狀態的有效和正確把握，也為飛控決策的制定提供了有力依據，已經成為一些重大任務不可或缺的軟件配置。

航天器的三維顯示以可視化的形式展示了其位置、姿態、相對關系、部件運動和工作狀況等狀態信息，其中姿態信息的展示是航天器運行狀態表現中非常重要的一部分[4-5]。由于一些任務飛控過程比較復雜，常伴隨著航天器的自主機動控制，以及為了減輕操作人員的負擔，實際任務中常需要使用實時采集的姿態數據來驅動三維顯示，以便更加真實地反映航天器運行狀態的變化。但由于設備異常、采集及傳輸誤碼、數據預處理異常和操作人員過失等原因，導致實時采集的姿態數據常含有異常或跳變數據。在數據處理中，一般將這種不合理或具有較大誤差的數據稱為野值。航天器的在軌飛行三維可視化過程中，野值數據的存在容易使顯示畫面出現跳變，影響實際顯示效果，因而在數據使用時首先需要對野值數據予以剔除。同時，為了應對數據中的丟點和跳變點等問題，還需對數據進行必要的插值和平滑處理，以免影響可視化表達的連續性。

對于姿態數據的野值剔除、平滑及插值等問題，可使用的方法較多，較常見的是采用卡爾曼濾波的方法[6-7]來進行野值剔除和平滑處理。但這類方法實現過程相對比較復雜，并且需要較繁瑣的濾波參數設置[8]。在一些可視化場合，可能不需要進行此類復雜的運算。為此，本文結合工程實際情況，采用一些較易實現的數值計算方法對三維可視化顯示中姿態數據的野值剔除及平滑插值等問題進行了實現。

1 姿態數據及可視化

姿態數據是航天器運行狀態的一項重要參數，描述了航天器本體相對基準坐標系的相對關系。對三軸穩定航天器，其姿態信息通常采用3個歐拉角[9-10](俯仰角、偏航角和滾動角)或四元數[11-12]等形式來表示。

航天器在軌運行過程的歐拉角姿態數據一般建立在質心軌道坐標系或當地東南坐標系下，四元數姿態一般相對于慣性系，這幾類姿態數據的可視化顯示方法可參見文獻[4]，且歐拉角與四元數2類表達形式之間可方便地進行相互轉換[9]。

在三維可視化的實時姿態數據使用過程中，姿態數據由于其不同于星歷、軌道等其他類數據的特殊性(轉序問題、可能不連續和不同表現形式等)，使其預處理較為復雜，下面將針對2類表達形式分別給出解決方法。

2 歐拉角姿態數據的處理

歐拉角是一類常見的姿態表達形式，包含3個角度：滾動角(γ)、偏航角(ψ)和俯仰角(φ)。對其進行可視化數據處理時可對3個角度分量分別進行處理。

2.1 野值剔除

野值數據通常情況下表現為脫離大多數數據點變化趨勢，幅值相差懸殊，點數較少。一般情形下，對歐拉角野值數據的粗處理可采用以下幾種方法。

2.1.1 粗大野值的閾值剔除

給定姿態數據某測量分量ν的閾值范圍[νmin，νmax](例如[-360°，360°]角度范圍)，當測量數據不在該范圍時直接剔除，這樣可將明顯的錯誤數據剔掉，減少后續處理的壓力。

閾值剔除是數據處理的第一步，但在其使用中需合理設置閾值范圍，不能設置過小以免將某些正常數據剔掉，使用時可根據理論飛行狀態及其變化情況進行設置。

2.1.2 差分檢測法

考慮到航天器姿態運動過程中，姿態角速度范圍是有限的，因而可通過前后兩點姿態角變化與姿態角速度的相關性比較來判定野值。設角速度某分量的閾值范圍為[ωmin，ωmax](可為負數)，相應姿態角分量前后兩點測量值為ν1，ν2，兩點間時間間隔為Δt，則應滿足ωmin·Δt≤ν2-ν1≤ωmax·Δt，否則，可認為后點為野值(正常往后遞推剔野時，非起步判斷過程)。

實際使用時，角速度閾值取值也應適當放大，以免將合理數據剔除。另外，t應進行適當的限定，即當t>tmaxin(tmaxin為連續判斷的時間間隔最大允許值)時，認為數據中斷，應重新進行野值的起步判斷。

2.1.3 預測值管道剔除法

基于一定時間內的觀測數據，可采用多項式擬合、濾波等方法預測出下一時刻的觀測數據，若實際測量值與預測值偏差較大，則可被視為野值。偏差比較采用閾值比較方式進行，此即預測值管道剔除法。管道法剔除遞推剔除示意圖如圖1所示。

圖1 管道法剔除遞推剔除示意

設由前n個數據預測出的下一測量數據為νF，實測量為ν，管道最大值為Δνmax，則當νF-ν>Δνmax時視為野值。

考慮到可能跨越±180°的情況，上式修改為：Trans(νF-ν)>Δνmax(此處Trans表示規約到±180°范圍內的轉換函數)。

上面幾種方法只是對粗大型野值的處理，實際計算中還需結合角度數據的平滑過程進行進一步處理。

2.2 歐拉角插值及補點處理

對于被當做野值剔除掉的數據點，如不進行數據補齊，就可能造成姿態顯示時畫面跳動，為此必須進行補點處理。同時，在一些平滑性要求更高的場合，可能需要一秒多幀顯示，此時就可能需要將數據插值為每秒多點。

工程上常采用的補點處理方法有2類：一類是進行數據插值；另一類是采用預測值代替(可采用多項式擬合預測、濾波預測等)。

歐拉角的插值主要采用2點線性插值、多點拉格朗日插值和樣條插值等方法。

由于工程計算中歐拉角一般被限定在一定范圍內(如[-180°，180°])變化，這就造成角度數據在形式上是不連續的(盡管實際幾何關系可能連續)，如值在±180°前后的變化過程[10]。因而插值時需要對其做連續性處理，方法為：

方法1：采用加減360°的方法，將數據規整為連續型數列；

方法2：三角函數方法，改用三角函數進行插值，然后反解出相應角度。

2.3 多項式最小二乘平滑處理

多項式擬合最小二乘估計處理是一種兼具預測、平滑及剔野功能的方法。

在歐拉角的多項式最小二乘處理使用前，首先需要采用前面的方法進行粗大型野值處理，以剔除掉明顯超出測量范圍的大野值。然后對數據點使用多項式f(t)進行平滑。可以固定每次平滑的數據點數，以保證算法運算速度不會由于數據的增多而減慢。在一次擬合后，計算出樣本方差σ，即

式中，νF為ν的預測值。

然后，計算出每個觀測數據點yi的偏差Ei=yi-f(ti)，進而對其進行方差檢驗(kσ檢驗，k為倍數，如取3)，剔除掉本次迭代計算的野值。

接下來重新進行擬合迭代計算，直到滿足收斂條件為止，收斂條件取為：

① 迭代次數少于給定值；

② 前后2次計算的擬合值相差小于某給定小量(如取1×10-5)。

需要注意的是，該方法的剔野效果有限，因而主要用來進行擬合預測、數據平滑計算以及被剔除點的補點計算等方面。

2.4 起步處理

在數據中斷一定時間后，或首次接收到數據時，需進行起步處理。起步處理的好壞直接影響了后續的差分剔野、擬合值管道剔野以及平滑處理的效果。故此處可采用多點數據起步。

首先進行粗大野值的閾值剔除，然后計算前后兩點的差分姿態角變化率。當不滿足限制條件ωmin·Δt≤ν2-ν1≤ωmax·Δt時，將ν2點舍棄或將其之前的數據全部舍棄(盡量保留新數據，與前面正常遞推過程差分剔野不同)，并繼續積攢數據，重新進行起步成功判斷。該方法的使用需要給定合理的起步時ωmin，ωmax范圍。

2.5 滑窗計算

在多點數據擬合預測、多項式平滑中，可采用多點數據滑窗計算方法。計算時，需注意數據的中斷情況，當前后兩點時標差超過一定時間(如50 s)時，認為數據中斷，需重新起步處理。

同時，當連續多幀(或一定時間段)數據均被當做野值剔除時，可認為當前剔野參數設置不合理，已經適應不了當前姿態的變化(即管道過小)，需對剔野門限放大(可采用連續多點角變化率平均值的多倍進行自適應計算，并進行適當限幅)。

3 四元數姿態數據的處理

歐拉角姿態存在旋轉不唯一、有轉序設定問題、旋轉不均勻以及運動方程存在奇點[10，13]等缺點。因而在角度變化較大的姿態計算中，常采用另一類表達形式，即四元數數據。

3.1 剔野處理

四元數為一個四維向量，其姿態旋轉軸由四元數矢部決定，旋轉角度由標部決定[14]。

設航天器在慣性系下的單位姿態四元數qm的4個分量為qm0，qm1，qm2，qm3，其中qm0為四元數的標部；qm1，qm2，qm3為矢部。則野值剔除可通過以下步驟進行：

① 計算四元數的模

若模不為1，則判定為異常值。用公式表示為：

qm-1>η(η為比較閾值，如取0.01)。

② 計算四元數矢部的模，若為0(可取為小于1×10-6)，則說明旋轉軸數據異常，將該數據剔除。

③ 計算前一個接收到的正常姿態數據qt-1到當前四元數qt的增量四元數(或稱作機動四元數)Δq：

計算增量四元數對應的轉角Δα：

Δα=2arccos(Δq[0])，

式中，Δq[0]表示取Δq的標部。

將此轉角與角速度閾值進行比較，設時間增量為Δt，則比較式為ωmin·Δt≤Δα≤ωmax·Δt，若不在此范圍，說明連續2個姿態數據的變化超出正常值，可判定當前四元數為異常值。

④ 序列數據的角速度擬合判斷

對一時間連續未發生中斷的數據序列，由上面計算的每一點與前一點的轉角計算其角速度，可得到一角速度序列，可將最后一點(最新點)計算的角速度與序列多項式擬合值進行閾值比較，判斷是否要將最新點剔除，剔除后用插值結果作為輸出(但原始測量值仍參與下一點數據剃野判斷，以減少機動時的誤判)。

3.2 插值計算

四元數姿態插值的一種簡易算法是先將其轉換為歐拉角姿態再進行插值。但這種算法最大缺陷在于若轉序選擇不好，可能導致歐拉角的奇異問題的發生[9，15]，從而放棄了四元數的優越性，且該方法會導致計算量增大，因而不推薦使用。

四元數插值[16]最簡單的方法是進行線性插值，對于2個四元數q1，q2，其線性插值公式為：

qt=q1+λ·(q2-q1)，0≤λ≤1，

式中，λ為q1變換到q2間的插值參數。插值計算中需注意的是：q1，q2應為單位四元數，且插值后的值也應轉換為單位四元數，將公式修改為：

該算法計算過程簡單，但缺點在于忽略了旋轉的球面特性，其插值過程不是勻速進行的，有加速和減速過程。為此，本文主要采用球面線性插值方法。

Ken Shoemaker[17]采用單位球面上大段圓弧類比歐式空間中直線的插值方法，提出了2個姿態間的球面線性插值(Spherical Linear Interpolation，Slerp插值)，其最初是為了模擬3D旋轉在四元數插值(Quaternion Interpolation)基礎上改進的方法，當給定末端且在0～1之間插值時，Slerp指的是沿著單位半徑大圓的常速運動。

給定t時刻左右兩端的四元數q1，q2(均為單位四元數，q1分量為q10，q11，q12，q13；q2分量為q20，q21，q22，q23)，對t時刻姿態(設為qλ)進行Slerp插值，基本計算原理如圖2所示。

圖2 Slerp插值原理示意

由圖2中幾何關系可推導出插值公式為：

式中，λ為插值參數，且λ∈[0，1]；θ為2個四元數之間的夾角：

θ= arccos(q1·q2)=

arccos(q10q20+q11q21+q12q22+q13q23)。

4 實驗結果

對某航天器姿態機動過程中某段時間內的實測姿態數據進行剔野[18-19]和平滑[20]處理，處理前后的偏航角曲線如圖3所示。

圖3 某段實測數據處理前后對比

從圖3中可見，大多數野值點均被有效剔除，且處理后的曲線較平滑。

原型軟件(采用VC與OpenGL編寫)對某深空探測器四元數姿態可視化效果如圖4所示。

圖4 原型軟件對探測器飛行過程姿態可視化效果

5 結束語

對航天器可視化中姿態數據的處理方法進行了闡述，實驗結果表明該方法是可行的。應該看到，野值剔除是一個比較困難的命題，各種處理方法均有一定的適用范圍，使用中需根據航天器運動情況適時調整參數才能得到更好的結果。下一步將重點在姿態數據的濾波處理方法上進行研究。