基于FADS/IMU信息融合的大氣數(shù)據(jù)測量方法研究

肖地波，蔣保睿，劉 鵬

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

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

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

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

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

1 系統(tǒng)設(shè)計

1.1 算法流程概述

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

圖1 算法流程圖

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

1.2 FADS與IMU數(shù)據(jù)預(yù)處理

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

圖2 開孔位置示意圖

FADS測壓孔位置參數(shù)如表1所示。

表 1 FADS測壓孔位置參數(shù)

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

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

(1)

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

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

(2)

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

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

圖3 攻角計算改進流程圖

圖3中TAO24為P2孔與P4孔壓力數(shù)據(jù)的差值。

(3)

(4)

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

1.3 CKF設(shè)計與計算

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

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

(5)

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

(6)

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

(7)

(8)

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

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

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

(9)

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

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

① CKF濾波初始化。

(10)

② 時間更新。

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

(11)

估計容積點為

(12)

式中：

(13)

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

估算狀態(tài)為

(14)

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

(15)

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

③ 量測更新。

將Pk∣k-1分解為

(16)

估計容積點為

(17)

估計量測預(yù)測值為

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 仿真參數(shù)設(shè)置

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

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

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

2.2 仿真結(jié)果分析

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

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

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

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

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

圖6 馬赫數(shù)誤差曲線

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

表2 不同方案估計馬赫數(shù)的誤差結(jié)果對比

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

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

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

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

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

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

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

3 結(jié)束語

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

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