竇立謙，杜苗苗，張秀云，王躍萍

(1. 天津大學(xué)電氣自動(dòng)化與信息工程學(xué)院，天津 300072；2. 飛行控制航空科技重點(diǎn)實(shí)驗(yàn)室(航空工業(yè)自控所)，西安 710065)

可重復(fù)使用運(yùn)載器不僅飛行速度快、高度范圍廣，而且其發(fā)射成本與傳統(tǒng)的天地往返系統(tǒng)相比較低[1-2]，是一種兼具航空器和航天器特點(diǎn)且能在未來(lái)天地往返的重要工具[3].但是，在其從太空返回地球的過(guò)程中，RLV將面臨再入段帶來(lái)的嚴(yán)峻考驗(yàn).在再入返回段中，由于RLV飛行速度快，飛行狀態(tài)變化劇烈，氣動(dòng)參數(shù)表現(xiàn)為不確定和快時(shí)變的特點(diǎn)[4]，導(dǎo)致RLV氣動(dòng)的耦合性和非線性增強(qiáng)，使氣動(dòng)模型難以設(shè)計(jì)和控制，降低了飛行器的穩(wěn)定性[5].因此，根據(jù)RLV再入返回段的量測(cè)數(shù)據(jù)，進(jìn)行氣動(dòng)參數(shù)辨識(shí)研究，對(duì) RLV的氣動(dòng)模型的建立和控制設(shè)計(jì)具有很重要的意義.

在已知模型結(jié)構(gòu)的前提下對(duì)氣動(dòng)參數(shù)進(jìn)行辨識(shí)，是根據(jù)試驗(yàn)數(shù)據(jù)估計(jì)模型中未知參數(shù)的過(guò)程.在實(shí)際飛行中，系統(tǒng)噪聲是往往不可忽略的，這便要求參數(shù)辨識(shí)算法的精確性.用于研究氣動(dòng)參數(shù)辨識(shí)的算法有最小二乘法、極大似然算法、卡爾曼濾波算法、擴(kuò)展卡爾曼濾波算法，其中卡爾曼濾波、擴(kuò)展卡爾曼濾波在參數(shù)辨識(shí)中應(yīng)用較為廣泛，但是，卡爾曼濾波適用于線性系統(tǒng)的狀態(tài)估計(jì)問(wèn)題，而 EKF適用于非線性系統(tǒng)的狀態(tài)及參數(shù)估計(jì)的問(wèn)題，且已廣泛應(yīng)用于飛行器氣動(dòng)參數(shù)辨識(shí)問(wèn)題的研究[6].

眾所周知，用于描述 EKF測(cè)量和過(guò)程噪聲的先驗(yàn)統(tǒng)計(jì)數(shù)據(jù)是需要假設(shè)為已知數(shù)據(jù)[7].這些噪聲特性的設(shè)置影響 EKF算法在參數(shù)估計(jì)時(shí)的成功或失敗.特別是在噪聲統(tǒng)計(jì)特性隨時(shí)間變化的環(huán)境中，EKF算法會(huì)導(dǎo)致較大的估計(jì)誤差甚至?xí)霈F(xiàn)偏差.因此，初始狀態(tài)、測(cè)量噪聲協(xié)方差矩陣和過(guò)程噪聲協(xié)方差矩陣的選擇，在控制 EKF學(xué)習(xí)算法的收斂中是非常重要的.故需要將噪聲的統(tǒng)計(jì)特性也視為未知參數(shù)，這便屬于存在隱含變量的參數(shù)辨識(shí)問(wèn)題.最大期望算法是一種用于含有隱變量參數(shù)模型的迭代的極大似然估計(jì)算法，Bavdekar等[8]曾經(jīng)將EM 算法用于非線性系統(tǒng)的參數(shù)估計(jì)，文獻(xiàn)[9]使用EM 算法對(duì)噪聲協(xié)方差進(jìn)行了精確的估計(jì).因此本文采用 EKF算法對(duì) RLV氣動(dòng)模型的未知?dú)鈩?dòng)參數(shù)進(jìn)行估計(jì)，使用 EM 算法對(duì)噪聲的統(tǒng)計(jì)特性進(jìn)行估計(jì)，并與 Rauch-Tung-Striebel(RTS)平滑器相結(jié)合，完成其最大化的求解過(guò)程，能夠在精確的噪聲特性的基礎(chǔ)上，更好地提高EKF算法參數(shù)估計(jì)的精度.

1 可重復(fù)使用運(yùn)載器(動(dòng)力學(xué))模型

考慮類X-37B構(gòu)型的RLV飛行器再入段六自由度方程[10-12]為

式中：x、y、z、v、α、β、φ、θ、ψ、p、q、r分別為 RLV的縱向位移、側(cè)向位移、飛行高度、飛行速度、攻角、側(cè)滑角、滾轉(zhuǎn)角、俯仰角、偏航角、滾轉(zhuǎn)角速度、俯仰角速度、偏航角速度；慣性積分別為沿x、y、z軸的轉(zhuǎn)動(dòng)慣量；L、D、Y分別為升力、阻力和側(cè)力；、、分別為 RLV的滾轉(zhuǎn)力矩、俯仰力矩和偏航力矩；其表達(dá)式為

式中aδ、eδ、δro分別為飛行器的副翼偏轉(zhuǎn)角、升降舵偏轉(zhuǎn)角和方向舵偏轉(zhuǎn)角.

由于飛行器力矩Cl、Cm、Cn的參數(shù)與氣動(dòng)力CL、CD、CY的參數(shù)辨識(shí)類似，因此，本文只對(duì)氣動(dòng)力部分參數(shù)進(jìn)行辨識(shí)，力矩參數(shù)的相關(guān)辨識(shí)過(guò)程不再贅述.

2 算法理論

2.1 非線性狀態(tài)空間模型

考慮非線性系統(tǒng)的狀態(tài)空間模型如下.

狀態(tài)方程和觀測(cè)方程分別為

式中：Xk、Yk分別為系統(tǒng)的狀態(tài)向量和觀測(cè)向量；η為待估計(jì)的未知參數(shù)；u為系統(tǒng)的輸入向量；f(·)為非線性函數(shù)；Γw為σ×p維的系統(tǒng)過(guò)程噪聲輸入矩陣；H為τ×n維的觀測(cè)矩陣；為e維的系統(tǒng)隨機(jī)過(guò)程噪聲序列；為τ維的系統(tǒng)隨機(jī)觀測(cè)噪聲序列.過(guò)程噪聲和觀測(cè)噪聲假設(shè)為滿足高斯零均值分布，且兩者相互獨(dú)立，即

式中：Q為過(guò)程噪聲wk的e× e維對(duì)稱正定方差矩陣；R為觀測(cè)噪聲kχ的τ×τ維對(duì)稱正定方差矩陣.

氣動(dòng)力參數(shù)主要與飛行器模型中式(1)的v、z、α、β相關(guān)，因此辨識(shí)時(shí)只考慮這 4個(gè)狀態(tài)方程即可.假設(shè)系統(tǒng)噪聲只對(duì)狀態(tài)量v、z、α、β產(chǎn)生干擾影響，角度及角速度等均可測(cè).式(1)中關(guān)于v、z、α、β對(duì)狀態(tài)方程變換如下.

狀態(tài)方程為

觀測(cè)方程可表示為

至此，類 X-37B RLV 模型已轉(zhuǎn)換為式(7)和式(8)所示的非線性系統(tǒng)的狀態(tài)空間模型，下文將針對(duì)此狀態(tài)空間模型，對(duì)氣動(dòng)參數(shù)進(jìn)行辨識(shí)研究.

2.2 擴(kuò)展卡爾曼濾波算法

擴(kuò)展卡爾曼濾波(EKF)是一種間接對(duì)參數(shù)進(jìn)行辨識(shí)估計(jì)的算法，其原理是將系統(tǒng)的參數(shù)估計(jì)問(wèn)題轉(zhuǎn)化為對(duì)系統(tǒng)狀態(tài)進(jìn)行估計(jì)，即要對(duì)原始的狀態(tài)向量進(jìn)行擴(kuò)維，得到由未知參數(shù)和原始狀態(tài)向量組成的擴(kuò)維狀態(tài)向量，即，從而得到新的狀態(tài)空間模型為

假設(shè)擴(kuò)維狀態(tài)向量的初值X0服從均值為μ、方差為P0的高斯分布，基于由式(13)和式(14)組成的擴(kuò)維狀態(tài)空間模型，使用擴(kuò)展卡爾曼濾波對(duì)擴(kuò)維狀態(tài)向量進(jìn)行預(yù)測(cè)和更新的過(guò)程如下.

1) EKF預(yù)測(cè)

狀態(tài)預(yù)測(cè)方程為

觀測(cè)預(yù)測(cè)方程為

協(xié)方差矩陣預(yù)測(cè)方程為

2) EKF更新

Kalman濾波增益更新公式為

狀態(tài)更新公式為

協(xié)方差更新公式為

由上述的 EKF過(guò)程，可以知道狀態(tài)估計(jì)的精度與擴(kuò)維狀態(tài)初值的分布特點(diǎn)(μ，P0)、噪聲的特性(Q，R)有著密切聯(lián)系，假設(shè)模型的統(tǒng)計(jì)量λ(μ，也是未知的，即在估計(jì)系統(tǒng)參數(shù)η的同時(shí)，需要通過(guò) EM 算法對(duì)λ進(jìn)行估計(jì)[13]，從而得到參數(shù)η更準(zhǔn)確的估計(jì)值.

2.3 EM-EKF算法

完成EM算法對(duì)上述參數(shù)λ的估計(jì)，首先需要獲得完整數(shù)據(jù)的概率函數(shù).假設(shè)系統(tǒng)的初始條件、狀態(tài)向量和觀測(cè)向量的變化均滿足高斯分布.由馬爾可夫鏈可得其表達(dá)式為

則對(duì)數(shù)似然函數(shù)的表達(dá)式為

3) 求期望

EM 是一種不斷迭代的優(yōu)化算法，假設(shè)當(dāng)前的迭代次數(shù)為i，令輔助函數(shù)為

其中，常數(shù)C=[K(σ+τ)+σ]ln(2π).

令

在第2.2節(jié)中所介紹的EKF算法流程為狀態(tài)的前向傳播過(guò)程，為了計(jì)算的值，需要應(yīng)用 Rauch-Tung-Strieber(RTS)平滑器的后向傳播過(guò)程，其原理如下.

式中：Jk為平滑器增益；均值初值為；互協(xié)方差初值為；方差初值為

故

又由于

則

同理可得

故輔助函數(shù)

為了得到參數(shù)λ的最佳估計(jì)值，需要對(duì)輔助函數(shù)極大化，即通過(guò)式(37)對(duì)參數(shù)μ、P0、Q、R求偏導(dǎo)數(shù)，并令其偏導(dǎo)數(shù)等于0，可得表達(dá)式為

下文將在由式(38)～式(41)得到 EKF初始狀態(tài)、測(cè)量和過(guò)程噪聲協(xié)方差矩陣最佳估計(jì)值的基礎(chǔ)上，采用EKF算法對(duì)未知?dú)鈩?dòng)參數(shù)η進(jìn)行精確估計(jì).

3 參數(shù)辨識(shí)結(jié)果與分析

考慮 RLV再入段飛行跨度大，范圍為 20～80km，氣動(dòng)起主要作用的區(qū)域范圍為 30～60km，且飛行器再入階段的攻角變化范圍較大.故以其飛行高度和攻角為基準(zhǔn)，每10km一個(gè)區(qū)間劃分飛行器的飛行階段，分別進(jìn)行氣動(dòng)參數(shù)的在線辨識(shí)，區(qū)間劃分如下.

區(qū)間 1：高度變化范圍為 50～60km，攻角范圍為 40°～45°，速度范圍為 4947～6281m/s.

區(qū)間 2：高度變化范圍為 40～50km，攻角范圍為 35°～40°，速度范圍為 3176～4947m/s.

區(qū)間 3：高度變化范圍為 30～40km，攻角范圍為 30°～35°，速度范圍為 1509～3176m/s.

由于在不同區(qū)間的狀態(tài)量、觀測(cè)量和待辨識(shí)的氣動(dòng)參數(shù)均一樣，只是在不同區(qū)間的初值和參數(shù)值發(fā)生改變，故只需介紹一個(gè)區(qū)間內(nèi)的氣動(dòng)參數(shù)辨識(shí)過(guò)程，下文將以區(qū)間1的氣動(dòng)參數(shù)辨識(shí)過(guò)程為例.

基于方程(11)與(12)可得：RLV 氣動(dòng)參數(shù)辨識(shí)的狀態(tài)向量，觀測(cè)向量Yk=，待辨識(shí)的氣動(dòng)參數(shù)，擴(kuò)維狀態(tài)向量，輸入向量

為了減少初值的設(shè)置對(duì)辨識(shí)結(jié)果的影響，選取待辨識(shí)參數(shù)初始值為 0，擴(kuò)維狀態(tài)向量的協(xié)方差P、過(guò)程噪聲和測(cè)量噪聲的方差矩陣Q、R的初值均為單位陣

狀態(tài)向量的初始值選取如下：

則擴(kuò)維狀態(tài)向量初始值為

系統(tǒng)常值參數(shù)如表1所示.

表1 系統(tǒng)常值參數(shù)Tab.1 System constant parameters

分別采用 EM-EKF算法(迭代次數(shù)為 2000)與極大似然(ML)算法對(duì)區(qū)間 1中的氣動(dòng)參數(shù)進(jìn)行辨識(shí)，其結(jié)果如表2所示.

表2 區(qū)間1參數(shù)辨識(shí)結(jié)果Tab.2 Parameter identification results of interval 1

同理可得，其他兩個(gè)區(qū)間的氣動(dòng)參數(shù)辨識(shí)值如表3和表4所示.

表3 區(qū)間2參數(shù)辨識(shí)結(jié)果Tab.3 Parameter identification result of interval 2

以機(jī)理模型的原始數(shù)據(jù)、EM-EKF(濾除噪聲)和極大似然(未濾除噪聲)算法(ML)的氣動(dòng)參數(shù)辨識(shí)值為基礎(chǔ)，基于區(qū)間 1分別分析攻角、側(cè)滑角、副翼偏轉(zhuǎn)角、升降舵偏轉(zhuǎn)角和方向舵偏轉(zhuǎn)角對(duì)升力、阻力、側(cè)力、俯仰力矩、滾轉(zhuǎn)力矩和偏航力矩系數(shù)的影響，具體對(duì)比效果如圖1～圖6所示.

表4 區(qū)間3參數(shù)辨識(shí)結(jié)果Tab.4 Parameter identification result of interval 3

從圖1～圖6的曲線可以看出，雖然三者之間的誤差不是很大，但是以 EM-EKF算法的氣動(dòng)參數(shù)辨識(shí)值為基礎(chǔ)的系數(shù)曲線能夠較好地與原始數(shù)據(jù)的曲線相擬合，即在系統(tǒng)噪聲存在的情況下，采用 EMEKF算法對(duì)氣動(dòng)參數(shù)具有較好的辨識(shí)效果.

圖2 阻力系數(shù)值的比較Fig.2 Comparison of the drag coefficient values

圖3 側(cè)力系數(shù)值的比較Fig.3 Comparison of the lateral force coefficient values

圖4 偏航力矩系數(shù)值的比較Fig.4 Comparison of the yawing-moment coefficient values

圖5 滾轉(zhuǎn)力矩系數(shù)值的比較Fig.5 Comparison of the rolling-moment coefficient values

圖6 俯仰力矩系數(shù)值的比較Fig.6 Comparison of the pitching moment coefficient values

4 結(jié) 語(yǔ)

本文針對(duì)RLV飛行器再入段氣動(dòng)特性的強(qiáng)耦合性和非線性等問(wèn)題，在系統(tǒng)噪聲存在的前提下，提出了一種基于 EM-EKF算法的氣動(dòng)參數(shù)辨識(shí)方法.首先，以飛行高度和速度為基準(zhǔn)，將RLV再入段分為3個(gè)飛行區(qū)間；然后，分別對(duì) 3個(gè)區(qū)間的氣動(dòng)參數(shù)進(jìn)行參數(shù)辨識(shí)；在辨識(shí)的過(guò)程中，考慮系統(tǒng)噪聲的存在，采用EKF算法對(duì)狀態(tài)和氣動(dòng)參數(shù)進(jìn)行估計(jì)以濾除噪聲，之后為了減少測(cè)量和過(guò)程噪聲的先驗(yàn)統(tǒng)計(jì)量的設(shè)置對(duì)EKF辨識(shí)結(jié)果帶來(lái)的影響，在EKF辨識(shí)算法的基礎(chǔ)上，結(jié)合 EM 算法更精確的估計(jì) EKF先驗(yàn)統(tǒng)計(jì)量，能夠更好地提高 EKF算法對(duì)氣動(dòng)參數(shù)估計(jì)的精度；最后，通過(guò)將提出的EM-EKF結(jié)合的參數(shù)辨識(shí)算法與極大似然方法的辨識(shí)值對(duì)各種氣動(dòng)系數(shù)影響進(jìn)行仿真對(duì)比，驗(yàn)證了辨識(shí)方法的準(zhǔn)確性.