肖地波，蔣保睿，劉 鵬

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

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

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

然而，由于氣壓傳感器的測(cè)量噪聲和延遲等問(wèn)題，導(dǎo)致其單獨(dú)使用時(shí)對(duì)數(shù)據(jù)的估計(jì)精度有限[6]，而慣性測(cè)量元件(Inertial Measurement Unit,IMU)的數(shù)據(jù)的瞬時(shí)精度較高、延遲較低，可以用于FADS的輔助計(jì)算。以IMU測(cè)量的加速度、角速度作為飛行器模型的輸入量，以FADS獲得的角度和速度信息作為觀測(cè)量構(gòu)建卡爾曼濾波是一個(gè)常用的方法[7]。因?yàn)轱w行器飛行時(shí)涉及到坐標(biāo)系轉(zhuǎn)換和當(dāng)?shù)芈曀僮兓挠?jì)算，所以需要對(duì)原有的卡爾曼濾波進(jìn)行改進(jìn)，以滿足非線性系統(tǒng)的計(jì)算要求。陸辰等[8]提出了一種基于 FADS 的擴(kuò)展卡爾曼濾波(Extended Kalman Filtering,EKF)實(shí)時(shí)估計(jì)大氣參數(shù)的方法，驗(yàn)證了慣性數(shù)據(jù)測(cè)量大氣數(shù)據(jù)的可靠性，同時(shí)建立了大氣數(shù)據(jù)傳感信息融合測(cè)試平臺(tái)，能提高大氣數(shù)據(jù)系統(tǒng)的精度和穩(wěn)定性。程超[9]采用卡爾曼濾波理論，設(shè)計(jì)了捷聯(lián)慣性導(dǎo)航系統(tǒng)/大氣數(shù)據(jù)系統(tǒng)組合導(dǎo)航濾波算法，并對(duì)算法的有效性進(jìn)行了仿真驗(yà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)航精度。

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

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

1.1 算法流程概述

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

圖1 算法流程圖

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

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

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

圖2 開(kāi)孔位置示意圖

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

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

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

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

(1)

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

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

(2)

化簡(jiǎn)后的方程組在迭代求解qc與P∞時(shí)計(jì)算更為簡(jiǎn)單。

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

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

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

(3)

(4)

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

1.3 CKF設(shè)計(jì)與計(jì)算

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

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

(5)

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

(6)

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

(7)

(8)

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

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

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

(9)

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

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

① CKF濾波初始化。

(10)

② 時(shí)間更新。

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

(11)

估計(jì)容積點(diǎn)為

(12)

式中：

(13)

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

估算狀態(tài)為

(14)

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

(15)

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

③ 量測(cè)更新。

將Pk∣k-1分解為

(16)

估計(jì)容積點(diǎn)為

(17)

估計(jì)量測(cè)預(yù)測(cè)值為

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

(18)

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

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

(19)

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

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

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

(20)

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

(21)

④ 濾波值輸出。

(22)

2 仿真實(shí)驗(yàn)

2.1 仿真參數(shù)設(shè)置

使用MATLAB對(duì)前述CKF算法進(jìn)行仿真實(shí)驗(yàn)，將濾波前測(cè)量結(jié)果的誤差、EKF算法[16]誤差和本文提出的CKF濾波的誤差進(jìn)行對(duì)比實(shí)驗(yà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°，初始時(shí)刻飛行器坐標(biāo)系與地面坐標(biāo)系朝向相同，飛行器坐標(biāo)系0點(diǎn)位于地面坐標(biāo)系的(0,0,-1000 m)處。飛行總時(shí)長(zhǎng)1000 s，采樣周期T=1 s。

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

2.2 仿真結(jié)果分析

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

圖4 不同濾波方法估計(jì)速度對(duì)比

圖5 不同濾波方法估計(jì)速度誤差對(duì)比

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

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

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

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

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

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

圖7 不同方案估計(jì)攻角的值與誤差

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

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

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

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

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

3 結(jié)束語(yǔ)

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

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