肖地波，蔣保睿，劉 鵬

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

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

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

然而，由于氣壓傳感器的測量噪聲和延遲等問題，導(dǎo)致其單獨使用時對數(shù)據(jù)的估計精度有限[6]，而慣性測量元件(Inertial Measurement Unit,IMU)的數(shù)據(jù)的瞬時精度較高、延遲較低，可以用于FADS的輔助計算。以IMU測量的加速度、角速度作為飛行器模型的輸入量，以FADS獲得的角度和速度信息作為觀測量構(gòu)建卡爾曼濾波是一個常用的方法[7]。因為飛行器飛行時涉及到坐標(biāo)系轉(zhuǎn)換和當(dāng)?shù)芈曀僮兓挠嬎?，所以需要對原有的卡爾曼濾波進(jìn)行改進(jìn)，以滿足非線性系統(tǒng)的計算要求。陸辰等[8]提出了一種基于 FADS 的擴(kuò)展卡爾曼濾波(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)航濾波算法，并對算法的有效性進(jìn)行了仿真驗證，證明了慣性導(dǎo)航系統(tǒng)和大氣數(shù)據(jù)系統(tǒng)兩者進(jìn)行融合是有效的。Nourmohammadi等[10]將容積卡爾曼濾波(Cubature Kalman Filter,CKF)用于慣性導(dǎo)航系統(tǒng)和衛(wèi)星導(dǎo)航系統(tǒng)的信息融合上，獲得了比EKF更高的導(dǎo)航精度。

針對高馬赫數(shù)飛行器的大氣數(shù)據(jù)測量與估計的方法進(jìn)行研究，提出了基于CKF的FADS/IMU耦合的大氣測量方法，利用大氣與飛行器飛行數(shù)據(jù)建立濾波算法模型，對飛行器的馬赫數(shù)、攻角、側(cè)滑角等重要大氣數(shù)據(jù)進(jìn)行高精度估計，并完成大氣數(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)估計精度更高。文獻(xiàn)[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)，其中一個點位于機(jī)頭中心，其余4個點均勻分布在四周，以球形機(jī)頭為例，開孔位置示意圖如圖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所示的攻角計算改進(jìn)流程進(jìn)行判斷與校正。

圖3 攻角計算改進(jìn)流程圖

圖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è)滑角進(jìn)行估計，結(jié)合測量的動、靜壓獲取較為全面的大氣參數(shù)。

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

(5)

大氣的基本風(fēng)場可以分為平均風(fēng)、大氣紊流、風(fēng)切變和突風(fēng)4種形式，其中前兩者最為普遍。平均風(fēng)是風(fēng)速的基準(zhǔn)值，表現(xiàn)為無風(fēng)或風(fēng)速、風(fēng)向不變，此時地速與空速的加速度一致，大氣紊流的時間短、速度及其改變量小，也可以認(rèn)為地速與真空速的加速度基本一致[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算法進(jìn)行仿真實驗，將濾波前測量結(jié)果的誤差、EKF算法[16]誤差和本文提出的CKF濾波的誤差進(jìn)行對比實驗，參數(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°，初始時刻飛行器坐標(biāo)系與地面坐標(biāo)系朝向相同，飛行器坐標(biāo)系0點位于地面坐標(biāo)系的(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，標(biāo)準(zhǔn)差約為20 Pa。

2.2 仿真結(jié)果分析

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

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

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

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

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

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

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

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

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

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

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

從仿真結(jié)果可以看出，當(dāng)攻角范圍在±50°，側(cè)滑角范圍在±20°的條件下，CKF濾波算法對飛行器的攻角和側(cè)滑角的估計效果更好，尤其在機(jī)動狀態(tài)時具有更強的穩(wěn)定性。在此條件下進(jì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è)滑角測量不準(zhǔn)確的問題，設(shè)計高精度和高可靠性為特色的濾波方法，并采用飛行曲線數(shù)據(jù)對其進(jìn)行測試。結(jié)果表明該算法能快速準(zhǔn)確地計算出各大氣參數(shù)：攻角誤差小于±0.1°，側(cè)滑角誤差小于±0.1°，馬赫數(shù)誤差小于±0.01,滿足飛行器大氣數(shù)據(jù)估計的基本精度要求。

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