肖地波，蔣保睿，劉 鵬

(成都信息工程大學 自動化學院，四川 成都 610225)

對于飛行器而言，傳統(tǒng)的空速管、攻角傳感器等大氣數據測量裝置在高速、高機動性的飛行時會產生較強干擾，且對飛行器的隱身效果有一定的影響。嵌入式大氣數據傳感(Flush Air Data Sensing,FADS)系統(tǒng)依靠機體表面的壓力分布，通過一系列算法間接獲得動靜壓、攻角、側滑角等大氣數據，具有維護成本低、經濟性良好等特點，被廣泛用于航空航天領域[1]。

目前，美國[2]與歐洲[3]對FADS系統(tǒng)開展了大量研究，并已開發(fā)出成熟產品。國內學者王鵬[4]研究了用于尖楔前體飛行器的FADS系統(tǒng)算法，并且對算法進行了驗證和測試。王禹等[5]提出了適合工程應用的飛翼布局飛機的FADS 算法模型。

然而，由于氣壓傳感器的測量噪聲和延遲等問題，導致其單獨使用時對數據的估計精度有限[6]，而慣性測量元件(Inertial Measurement Unit,IMU)的數據的瞬時精度較高、延遲較低，可以用于FADS的輔助計算。以IMU測量的加速度、角速度作為飛行器模型的輸入量，以FADS獲得的角度和速度信息作為觀測量構建卡爾曼濾波是一個常用的方法[7]。因為飛行器飛行時涉及到坐標系轉換和當地聲速變化的計算，所以需要對原有的卡爾曼濾波進行改進，以滿足非線性系統(tǒng)的計算要求。陸辰等[8]提出了一種基于 FADS 的擴展卡爾曼濾波(Extended Kalman Filtering,EKF)實時估計大氣參數的方法，驗證了慣性數據測量大氣數據的可靠性，同時建立了大氣數據傳感信息融合測試平臺，能提高大氣數據系統(tǒng)的精度和穩(wěn)定性。程超[9]采用卡爾曼濾波理論，設計了捷聯(lián)慣性導航系統(tǒng)/大氣數據系統(tǒng)組合導航濾波算法，并對算法的有效性進行了仿真驗證，證明了慣性導航系統(tǒng)和大氣數據系統(tǒng)兩者進行融合是有效的。Nourmohammadi等[10]將容積卡爾曼濾波(Cubature Kalman Filter,CKF)用于慣性導航系統(tǒng)和衛(wèi)星導航系統(tǒng)的信息融合上，獲得了比EKF更高的導航精度。

針對高馬赫數飛行器的大氣數據測量與估計的方法進行研究，提出了基于CKF的FADS/IMU耦合的大氣測量方法，利用大氣與飛行器飛行數據建立濾波算法模型，對飛行器的馬赫數、攻角、側滑角等重要大氣數據進行高精度估計，并完成大氣數據測量仿真分析。

1 系統(tǒng)設計

1.1 算法流程概述

基于CKF的FADS/IMU數據融合大氣數據估計系統(tǒng)包括直接測量大氣數據的FADS解算、IMU數據解算和CKF信息融合估計。算法流程如圖 1所示。

圖1 算法流程圖

CKF是貝葉斯濾波中的一種可廣泛應用于非線性濾波問題的濾波器，相比于UKF算法，CKF避免了因矩陣分解可能會出現的矩陣奇異問題，而且對于高維非線性系統(tǒng),其狀態(tài)估計精度更高。文獻[11]經過分析對比得出：當系統(tǒng)狀態(tài)維數大于3時，CKF算法的估計精度高于UKF算法。

1.2 FADS與IMU數據預處理

研究的FADS算法通過飛行器表面特定區(qū)域的壓力分布反推得到飛行參數，解算順序一般為：動靜壓、攻角側滑角、馬赫數和氣壓高度。構建以5個測壓點為輸入的FADS系統(tǒng)，其中一個點位于機頭中心，其余4個點均勻分布在四周，以球形機頭為例，開孔位置示意圖如圖2所示。

圖2 開孔位置示意圖

FADS測壓孔位置參數如表1所示。

表 1 FADS測壓孔位置參數

表1中a為圓周角；b為圓錐角。

由測壓孔的位置參數和每個孔的壓力值，可列出方程求解動靜壓：

(1)

式中：pi為第i個測壓孔在當前時刻所測壓力；σi為第i個孔對應的氣流方向；孔的數量為5，即n=5；f為一個單調增函數，其自變量為動靜壓之比，因變量為修正系數，可由實驗或牛頓理論確定，用于修正高馬赫數下測壓孔數據的偏差。利用5個方程組成的方程組迭代即可求解qc與P∞兩個未知數。

為了減少計算量，可使用矩陣偽逆的方法將式(1)化簡為

(2)

化簡后的方程組在迭代求解qc與P∞時計算更為簡單。

基于均勻分布的測壓孔，可以使用“三點法”求解攻角和側滑角數據，但由于三點法涉及反三角函數的計算，在大攻角(α>45°)與小攻角(α<45°)的計算方法不同，原有經典方法[12]僅說明了不同攻角范圍時的結算步驟，但未給出相應的判斷方法，可以使用如圖3所示的攻角計算改進流程進行判斷與校正。

圖3 攻角計算改進流程圖

圖3中TAO24為P2孔與P4孔壓力數據的差值。

(3)

(4)

式中:l11為L第1行第1列的數;l12為L第1行第2列的數;以此類推。

1.3 CKF設計與計算

本文提出一種融合IMU數據與FADS數據的大氣數據參數估計算法。算法以飛行器飛行速度、加速度、姿態(tài)角、角速度為狀態(tài)量建立系統(tǒng)狀態(tài)方程；以FADS所測的空速為狀態(tài)量，建立系統(tǒng)量測方程；考慮到狀態(tài)方程和量測方程均有較強的非線性特性,采用CKF對馬赫數、攻角、側滑角進行估計，結合測量的動、靜壓獲取較為全面的大氣參數。

選取飛行器的三軸速度分量作為狀態(tài)量X，即：

(5)

大氣的基本風場可以分為平均風、大氣紊流、風切變和突風4種形式，其中前兩者最為普遍。平均風是風速的基準值，表現為無風或風速、風向不變，此時地速與空速的加速度一致，大氣紊流的時間短、速度及其改變量小，也可以認為地速與真空速的加速度基本一致[13]，即：

(6)

經整理可得狀態(tài)方程：

(7)

(8)

式中：L為式(3)中的旋轉矩陣。

由于地速與真空速的加速度一致，Z可直接由IMU的數據輸出，即IMU可提供所需的量測信息。其關系如下：

Z=[X(1),X(2),X(3)]T+v

(9)

式中:X(1)、X(2)、X(3)分別為狀態(tài)量X的第1～第3項；v為量測噪聲向量。因為IMU的數據經慣導系統(tǒng)修正后通常都可以提供較高的角速度，所以此處旋轉矩陣直接使用狀態(tài)量中的(θ,φ,ψ)。

CKF 濾波采用三階容積法則，用數值積分來近似高斯加權積分，利用一組等權值容積點加權求和來代替加權高斯問題，尤其在高維非線性系統(tǒng)中，可以獲得較高的估計精度[14]。針對式(5)～式(9)的非線性系統(tǒng)，CKF濾波算法具體流程如下[15]。

① CKF濾波初始化。

(10)

② 時間更新。

對于k-1時刻的狀態(tài)濾波誤差陣，將其因式分解為

(11)

估計容積點為

(12)

式中：

(13)

m為系統(tǒng)向量維數的2倍，這里即為18。

估算狀態(tài)為

(14)

估計誤差協(xié)方差預測值為

(15)

式中：Qk-1為第k-1時刻的系統(tǒng)過程噪聲的協(xié)方差矩陣。

③ 量測更新。

將Pk∣k-1分解為

(16)

估計容積點為

(17)

估計量測預測值為

Zi,k∣k-1=h(Xi,k∣k-1)

(18)

式中：h(*)為觀測方程，用于取狀態(tài)量X的第1～第3維。

估計新息協(xié)方差矩陣為

(19)

式中：Rk為第k時刻的測量噪聲的協(xié)方差矩陣。

估計互協(xié)方差矩陣與增益為

[Xi,k∣k-1-Zi,k∣k-1]T

(20)

估計誤差協(xié)方差為

(21)

④ 濾波值輸出。

(22)

2 仿真實驗

2.1 仿真參數設置

使用MATLAB對前述CKF算法進行仿真實驗，將濾波前測量結果的誤差、EKF算法[16]誤差和本文提出的CKF濾波的誤差進行對比實驗，參數設置如下。

飛行高度起始值為1 km，0～200 s(上升段)內上升至巡航高度20 km并保持500 s，最后700～1000 s(下降段)下降至1 km；馬赫數上升段從0.1提升至7，巡航段保持不變，下降段由7降為0.1；初始航向角為0°，初始時刻飛行器坐標系與地面坐標系朝向相同，飛行器坐標系0點位于地面坐標系的(0,0,-1000 m)處。飛行總時長1000 s，采樣周期T=1 s。

為驗證論文算法，設置仿真條件中IMU和FADS的噪聲：加速度計0.06 m/s2(3σ)；陀螺儀0.03°/s2(3σ)；各孔壓力數值由馬赫數、攻角、側滑角等參數的實際值求解并添加噪聲而得；噪聲由加性噪聲和乘性噪聲組成，服從一階高斯-馬爾科夫過程，相關時間系數0.5，標準差約為20 Pa。

2.2 仿真結果分析

對于上述仿真模型，分別采用前述CKF和文獻所述的EKF進行參數估計，圖 4和圖 5為各濾波方法下濾波前后飛行器三軸速度分量和誤差分量的對比。

圖4 不同濾波方法估計速度對比

圖5 不同濾波方法估計速度誤差對比

兩種濾波方法對速度均有很好的估計能力。但隨著時間增加，EKF在線性化過程中泰勒展開導致非線性部分數據丟失，出現輸出結果抖動較大的情況，在100 s時即出現較強的振蕩。

兩種濾波算法估計的馬赫數誤差曲線如圖 6所示。在整個仿真過程中，CKF的效果都更好。EKF 由于采用了一階近似的泰勒逼近方法，只能對非線性的系統(tǒng)做到粗略的近似，損失部分精度，且此現象在飛行器高機動性飛行時尤其顯著，所以在圖6中100 s附近EKF估計的馬赫數產生了較大誤差。

圖6 馬赫數誤差曲線

在馬赫數為0.1～7的變化范圍下，進行100次的仿真實驗，統(tǒng)計各方法估計馬赫數誤差的最大誤差、平均誤差和誤差標準差，結果對比如表 2所示。

表2 不同方案估計馬赫數的誤差結果對比

圖7和圖8分別為不同濾波算法對攻角與側滑角的估計效果。

圖7 不同方案估計攻角的值與誤差

圖8 不同方案估計側滑角的值與誤差

從仿真結果可以看出，當攻角范圍在±50°，側滑角范圍在±20°的條件下，CKF濾波算法對飛行器的攻角和側滑角的估計效果更好，尤其在機動狀態(tài)時具有更強的穩(wěn)定性。在此條件下進行100次的仿真實驗，每次實驗持續(xù)1000 s，采樣周期1 s，表 3和表 4為所有采樣點得到的誤差的統(tǒng)計結果。

表3 EKF估計α和β的誤差結果對比(100次)

表4 CKF估計α和β的誤差結果對比(100次)

從表4中可以看出，利用CKF濾波算法估計的攻角最大誤差和側滑角最大誤差，比EKF所得最大誤差分別減少了39%和47%。

3 結束語

以容積卡爾曼濾波的FADS和IMU為研究對象，針對現有飛行器面臨的馬赫數、攻角、側滑角測量不準確的問題，設計高精度和高可靠性為特色的濾波方法，并采用飛行曲線數據對其進行測試。結果表明該算法能快速準確地計算出各大氣參數：攻角誤差小于±0.1°，側滑角誤差小于±0.1°，馬赫數誤差小于±0.01,滿足飛行器大氣數據估計的基本精度要求。

但是濾波算法在使用時忽略了風場變化等外界擾動帶來的不確定性誤差，后續(xù)需要進一步研究此類偏差存在時的估計方法以提高算法的普遍適用性。