臧劍文, 畢曉燁, 金 釗, 劉 浩, 閆 明

(1. 大連理工大學(xué)航空航天學(xué)院, 遼寧 大連 116024; 2. 沈陽(yáng)飛機(jī)設(shè)計(jì)研究所, 遼寧 沈陽(yáng) 110035)

0 引 言

近年來(lái),隨著系統(tǒng)辨識(shí)理論的研究,飛行器氣動(dòng)參數(shù)辨識(shí)在飛行試驗(yàn)中的應(yīng)用得到了飛躍發(fā)展,并在飛行器設(shè)計(jì)中起著越來(lái)越重要的作用,采用參數(shù)辨識(shí)方法來(lái)確定飛機(jī)參數(shù),能夠準(zhǔn)確并迅速地將真實(shí)氣動(dòng)特性從試驗(yàn)結(jié)果中分離出來(lái),對(duì)縮短數(shù)據(jù)處理時(shí)間和減少飛行試驗(yàn)周期等具有重要應(yīng)用價(jià)值[1]。

在大迎角機(jī)動(dòng)時(shí),由于機(jī)身機(jī)翼表面的氣流經(jīng)歷了附著流-旋渦流-渦破裂的歷程,這才使得飛機(jī)氣動(dòng)特性呈現(xiàn)強(qiáng)非線性等特征。傳統(tǒng)的線性辨識(shí)方法對(duì)大迎角機(jī)動(dòng)非線性氣動(dòng)參數(shù)辨識(shí)精度較低,因此線性辨識(shí)模型不適用于大迎角機(jī)動(dòng)情況。

張家銘等[2-3]提出了基于機(jī)器學(xué)習(xí)的氣動(dòng)參數(shù)辨識(shí)方法,但該方法需要足夠多的數(shù)據(jù)樣本,不適用于在線辨識(shí)。李正強(qiáng)等[4]提出基于最小二乘支持向量機(jī)(least squares support vector machines,LS-SVM)飛機(jī)大迎角動(dòng)態(tài)辨識(shí)方法,該方法的弊端在于LS-SVM模型的學(xué)習(xí)時(shí)間較長(zhǎng),支持向量機(jī)的速度較慢,不適用于在線辨識(shí)。喬偉等[5-6]提出基于遞推最小二乘法的飛行器模型參數(shù)在線辨識(shí),但弊端在于無(wú)法直接用于非線性復(fù)雜運(yùn)動(dòng)過(guò)程。蘇振宇[7-8]提出基于試飛數(shù)據(jù)的飛機(jī)大迎角氣動(dòng)力參數(shù)辨識(shí)方法,但弊端在于人為將迎角分割為若干小區(qū)間,不具備自動(dòng)分區(qū)的能力。張婉鑫等[9]提出大迎角非定常氣動(dòng)參數(shù)辨識(shí)方法,該方法在數(shù)據(jù)處理的過(guò)程中要求通過(guò)插值法求解二階導(dǎo)數(shù)、三階導(dǎo)數(shù)等,數(shù)據(jù)存在較大的誤差。

飛機(jī)大幅度機(jī)動(dòng)飛行時(shí),氣動(dòng)特性表現(xiàn)出非線性特征,特別是當(dāng)飛機(jī)具有高迎角姿態(tài)時(shí),數(shù)據(jù)愈加離散,傳統(tǒng)的辨識(shí)方法擬合效果較差[10]。本文基于“局部線性化”的思想,提出大迎角數(shù)據(jù)實(shí)時(shí)自動(dòng)分區(qū)的辨識(shí)方法,將大范圍的迎角數(shù)據(jù)自動(dòng)劃分為幾個(gè)小范圍的子區(qū)域,將一個(gè)復(fù)雜的非線性建模問(wèn)題分解為若干個(gè)局部線性問(wèn)題,進(jìn)而辨識(shí)每個(gè)子區(qū)域內(nèi)的氣動(dòng)參數(shù)。該方法與傳統(tǒng)的參數(shù)辨識(shí)方法相比,辨識(shí)模型結(jié)構(gòu)相對(duì)簡(jiǎn)單,且適用于非線性復(fù)雜運(yùn)動(dòng)等情況[11];與深度學(xué)習(xí)的方法相比,不需大量的數(shù)據(jù)樣本進(jìn)行神經(jīng)網(wǎng)絡(luò)的訓(xùn)練,可應(yīng)用于在線辨識(shí)[12];與插值法求解高階導(dǎo)數(shù)的方法相比,數(shù)據(jù)精度較高,有利于提高后續(xù)的辨識(shí)精度[13]。

1 飛行器參數(shù)辨識(shí)模型

在參數(shù)辨識(shí)之前需要建立辨識(shí)模型,由于飛行器參數(shù)之間存在強(qiáng)耦合,導(dǎo)致辨識(shí)難度大。本文將飛行器的運(yùn)動(dòng)分解為縱向運(yùn)動(dòng)和橫側(cè)向運(yùn)動(dòng),以縱向運(yùn)動(dòng)模型中俯仰力矩系數(shù)、升阻力系數(shù)為例,建立動(dòng)力學(xué)參數(shù)辨識(shí)模型[14]。

1.1 縱向運(yùn)動(dòng)模型

把飛行器運(yùn)動(dòng)分解為縱向運(yùn)動(dòng)和橫側(cè)向運(yùn)動(dòng),并以兩組相互獨(dú)立的運(yùn)動(dòng)方程組來(lái)描述[15]。通過(guò)對(duì)飛行器縱向受力進(jìn)行分析,得到如下微分方程組[16]:

(1)

式中:m為飛行器質(zhì)量;P為發(fā)動(dòng)機(jī)推力;X為空氣阻力;Y為升力;ωz、Mz為俯仰角速率、俯仰力矩;x、y為x軸、y軸方向上的位移;Jz為轉(zhuǎn)動(dòng)慣量;α、?、θ為迎角、俯仰角和軌跡傾角。

1.2 觀測(cè)模型

在真實(shí)的飛行過(guò)程中,升阻力系數(shù)、俯仰力矩系數(shù)不能通過(guò)傳感器直接測(cè)得,而是根據(jù)其他可觀測(cè)的參數(shù)計(jì)算得到,由此需要建立可觀測(cè)參數(shù)與力矩系數(shù)的觀測(cè)模型,為后續(xù)的辨識(shí)過(guò)程提供待辨識(shí)參數(shù)的觀測(cè)值[17]。

(2)

(3)

空氣動(dòng)力學(xué)的整體模型的理想形式是同時(shí)具有估測(cè)的模型參數(shù)值和相關(guān)不確定性的數(shù)學(xué)模型結(jié)構(gòu),可以將無(wú)量綱的空氣動(dòng)力學(xué)中的力和力矩系數(shù)與可以測(cè)量的飛機(jī)狀態(tài)和控制聯(lián)系起來(lái)。所有的全局辨識(shí)都是基于時(shí)域的方程誤差最小二乘法來(lái)進(jìn)行建模。在此方法中,因變量是無(wú)量綱力或力矩系數(shù)中的某一項(xiàng),它是通過(guò)利用解釋變量(無(wú)量綱的飛機(jī)狀態(tài)量和控制面撓度)計(jì)算出非線性模型項(xiàng)的擴(kuò)展來(lái)進(jìn)行建模的[18]。

因此,可以使用特定模型結(jié)構(gòu)來(lái)表達(dá)俯仰力矩系數(shù)的方程誤差最小二乘問(wèn)題。設(shè)置建模函數(shù)如下:

(4)

式(4)可以寫(xiě)成:

z=Hθ+v

(5)

其中,z代表觀測(cè)模型計(jì)算的觀測(cè)值:

(6)

H代表設(shè)置的建模函數(shù):

(7)

θ代表辨識(shí)模型:

(8)

v代表殘差:

(9)

2 動(dòng)力學(xué)參數(shù)辨識(shí)方法

2.1 基于軸正交決策樹(shù)的迎角自動(dòng)分區(qū)

回歸樹(shù)是一種具有邏輯性的分類方法,可以為局部建模網(wǎng)絡(luò)劃分權(quán)重變量,并分出新的單元結(jié)構(gòu),該過(guò)程將解釋變量歸入不同的區(qū)域,來(lái)表征解釋變量的非線性特征。分區(qū)方法通常采用二進(jìn)制拆分法則,圖1所示二進(jìn)制軸正交決策樹(shù)網(wǎng)絡(luò)結(jié)構(gòu),在具有兩個(gè)權(quán)重變量的情況下,在每個(gè)級(jí)別的單元格中再拆分出一個(gè)新單元格。被劃分的單元稱為母單元,隨后拆分出的單元稱為子單元[19]。

圖1 二進(jìn)制軸正交決策樹(shù)網(wǎng)絡(luò)圖

一個(gè)回歸樹(shù)對(duì)應(yīng)著對(duì)輸入空間的一個(gè)劃分,以及在劃分的單元上的輸出值。假設(shè)把輸入空間劃分為M個(gè)互不相交的單元R1,R2,…,RM,且在每個(gè)單元Rm上都有固定的輸出值Cm(m=1,2,…,M),則回歸樹(shù)模型[20]可表示為

(10)

可迭代的二進(jìn)制軸正交回歸樹(shù)算法已成功用于各種非線性靜態(tài)和動(dòng)態(tài)建模應(yīng)用程序中。

假設(shè):

(11)

對(duì)于任意的t,總要消除y(t)對(duì)其中一個(gè)xi(t)的依賴性,比如xi(t)取xp(t)。通過(guò)分區(qū),可以在{x1(t)x2(t) …xn(t)}的真子集上重新定義y(t)[21]為

(12)

2.2 飛行器在線辨識(shí)激勵(lì)信號(hào)設(shè)計(jì)

實(shí)時(shí)氣動(dòng)建模時(shí)需要迅速地收集到所有剛體自由度的數(shù)據(jù)信息,而在一系列飛行條件下,將自動(dòng)正交優(yōu)化的多正弦擾動(dòng)施加在控制舵面上可以很有效地收集相關(guān)氣動(dòng)數(shù)據(jù)[23]。使用擾動(dòng)輸入來(lái)激勵(lì)飛行過(guò)程,該輸入通過(guò)將設(shè)計(jì)好的激勵(lì)信號(hào)與控制指令中的執(zhí)行器命令相加,以此作為飛行器的舵面指令。

設(shè)定應(yīng)用在第j個(gè)控制面的擾動(dòng)輸入為uj,其為具有單個(gè)相移的正弦波諧波φk的總和[24]。

(13)

式中:M是可用的諧波相關(guān)頻率的總數(shù);T是激勵(lì)的時(shí)間長(zhǎng)度;Ak是正弦波分量的振幅;t是時(shí)間分量;nj個(gè)輸入中的每一個(gè)都是從M個(gè)諧波正弦波中選擇具有頻率ωk=2πk/T的分量,k=1,2,…,M。ωm=2πM/T代表激勵(lì)輸入的頻帶上限。區(qū)間[ω1-ωm]代表了指定期望飛機(jī)動(dòng)力學(xué)的頻率范圍。

為了實(shí)現(xiàn)均勻的功率分配,Ak被設(shè)置為

(14)

式中:n是包含在等式(13)中的正弦分量的數(shù)量;A是輸入的uj幅度。

2.3 時(shí)域遞推最小二乘辨識(shí)方法

使用遞推最小二乘法進(jìn)行參數(shù)辨識(shí)時(shí),首先要利用已知觀測(cè)量和輸出量計(jì)算遞推初值。在已知k時(shí)刻之前所有的觀測(cè)和輸出時(shí),記B(k)=HT(k)H(k),則前k時(shí)刻所計(jì)算的遞推初值為

θ(k)=B(-1)(k)HT(k)z(k)

(15)

在實(shí)際使用的過(guò)程中,隨著信息的增加,信息矩陣的正定性不斷減小,對(duì)新的信息改進(jìn)作用逐步等于零,這一現(xiàn)象稱為數(shù)據(jù)飽和現(xiàn)象[25]。為解決這一問(wèn)題,提出了引入遺忘因子的遞推算法。遞推公式[26]如下:

(16)

經(jīng)式(16)遞推計(jì)算得到θ4×N矩陣:

(17)

z=Hθ4×N

(18)

由式(18)可知:

(19)

2.4 基于加權(quán)函數(shù)的數(shù)據(jù)平滑方法

當(dāng)每個(gè)局部模型辨識(shí)氣動(dòng)參數(shù)以后,需要對(duì)全局范圍內(nèi)氣動(dòng)參數(shù)進(jìn)行數(shù)據(jù)平滑整合,該部分基于加權(quán)函數(shù)對(duì)局部模型的臨界數(shù)據(jù)進(jìn)行數(shù)據(jù)平滑[27]。具體步驟如下:

(1) 對(duì)第k個(gè)區(qū)域進(jìn)行有效性分析:根據(jù)之前確定的權(quán)重變量,φ=[φ1,φ2,…,φnφ]利用下式計(jì)算該區(qū)域的有效性:

(20)

Ck,j是第j個(gè)權(quán)重變量在該單元格內(nèi)的中心點(diǎn)的高斯函數(shù),計(jì)算公式如下[28]:

(21)

式中:μ代表在該單元內(nèi)的該分區(qū)變量的均值;σ2代表方差。值得注意的是,由于權(quán)重變量是隨著數(shù)據(jù)點(diǎn)的不同而不斷變化的,而在每段區(qū)域中心處求該區(qū)域中心點(diǎn)的高斯中心值是固定的,因此根據(jù)上述公式得當(dāng)某點(diǎn)離某一區(qū)域的高斯中心越遠(yuǎn)時(shí),該區(qū)域在這點(diǎn)處有效性越小[29]。所以,比如對(duì)第一個(gè)區(qū)域內(nèi)的某一數(shù)據(jù)點(diǎn)來(lái)說(shuō),主要是第一個(gè)區(qū)域在該點(diǎn)有效性比較大,離其他區(qū)域越遠(yuǎn),則其在該點(diǎn)有效性越小。

(2) 根據(jù)上述公式得到各個(gè)區(qū)域的有效性后,通過(guò)式(22)來(lái)計(jì)算第k個(gè)單元的權(quán)重[30]。

(22)

(3) 利用計(jì)算出的各個(gè)單元格的權(quán)重,進(jìn)行全局平滑建模。對(duì)每個(gè)點(diǎn)的整體建模,采用各個(gè)單元內(nèi)區(qū)域建模的加權(quán)加和得到[31]:

(23)

首先,平滑的依據(jù)有兩個(gè):① 在斷點(diǎn)處的值相等;② 在斷點(diǎn)處導(dǎo)數(shù)相等。

(24)

則在臨界點(diǎn)的平滑計(jì)算值為

(25)

證畢

以上平滑的兩點(diǎn)依據(jù)均滿足,因此可以認(rèn)為該曲線已經(jīng)平滑。通過(guò)上述過(guò)程,可以將各個(gè)分區(qū)的建模結(jié)果進(jìn)行平滑加權(quán)加和,從而得到一個(gè)平滑后的有效的整體模型。

3 仿真實(shí)例

為驗(yàn)證本文方法的有效性,針對(duì)飛行高度7 500 m、Ma0.4的某飛行器飛行數(shù)據(jù),建立縱向運(yùn)動(dòng)模型、辨識(shí)觀測(cè)模型以及遞推最小二乘辨識(shí)模型,總仿真時(shí)間為ts=25 s,仿真步長(zhǎng)為tk=0.01 s。飛行器的相關(guān)參數(shù)如表1所示。飛行器的飛行狀態(tài)數(shù)據(jù)如圖2所示。

表1 飛行器相關(guān)參數(shù)

圖2 大迎角機(jī)動(dòng)仿真數(shù)據(jù)

3.1 仿真計(jì)算結(jié)果及分析

圖3展示仿真流程示意圖,仿真流程主要包括3個(gè)方面:基于激勵(lì)信號(hào)的數(shù)據(jù)收集,建立觀測(cè)模型、空氣動(dòng)力學(xué)模型、最小二乘辨識(shí)模型,基于加權(quán)函數(shù)的數(shù)據(jù)平滑處理。

圖3 仿真流程示意圖

本文將激勵(lì)信號(hào)設(shè)計(jì)為具有特定諧波頻率,優(yōu)化相移和指定功率分布之后的正弦波之和[32],如圖4所示。

圖4 正弦激勵(lì)信號(hào)

公式表達(dá)如式(26),其中T=5 s。

(26)

選擇特定的分量頻率以覆蓋飛行器姿態(tài)運(yùn)動(dòng)包含的頻帶[33]。另外,將激勵(lì)信號(hào)設(shè)計(jì)為多元輸入,保證其在時(shí)域和頻域中相互正交,以使其在短時(shí)間內(nèi)所有軸上都具有較高的數(shù)據(jù)信息含量[34]。

在充分獲取數(shù)據(jù)的基礎(chǔ)上,建立俯仰力矩系數(shù)、升力系數(shù)的觀測(cè)模型如式(27)。取迎角變化范圍為5°,基于回歸樹(shù)迎角自動(dòng)分區(qū)結(jié)果如圖5所示。

圖5 回歸樹(shù)迎角自動(dòng)分區(qū)結(jié)果

(27)

根據(jù)每個(gè)分區(qū)的辨識(shí)情況,總結(jié)了計(jì)算俯仰力矩系數(shù)、升力系數(shù)各項(xiàng)參數(shù)的辨識(shí)結(jié)果。由圖6和圖7可知,氣動(dòng)導(dǎo)數(shù)在5 s左右已經(jīng)收斂,且收斂速度較快。氣動(dòng)導(dǎo)數(shù)的辨識(shí)精度分別為2.1%、2.9%、1.6%,詳見(jiàn)表2。

表2 俯仰力矩系數(shù)辨識(shí)精度

圖6 俯仰靜穩(wěn)定系數(shù)和升降舵效率辨識(shí)結(jié)果

由圖8和圖9可見(jiàn),升力系數(shù)的氣動(dòng)導(dǎo)數(shù)辨識(shí)收斂較快,辨識(shí)精度分別為0.6%、0.5%,詳見(jiàn)表3。

表3 升力系數(shù)辨識(shí)精度

圖8 迎角系數(shù)和升降舵系數(shù)辨識(shí)結(jié)果

圖9 升力系數(shù)辨識(shí)結(jié)果

為充分證明迎角分區(qū)時(shí)域辨識(shí)精度較傳統(tǒng)全局辨識(shí)精度高,下文闡述迎角分區(qū)與傳統(tǒng)全局的辨識(shí)結(jié)果對(duì)比,如圖10所示。

圖10 俯仰靜穩(wěn)定系數(shù)和升降舵效率辨識(shí)結(jié)果

由表4可知,利用迎角分區(qū)辨識(shí)的方法,辨識(shí)俯仰靜穩(wěn)定系數(shù)、升降舵操縱效率以及俯仰阻尼力矩系數(shù),其相對(duì)誤差均在3%以內(nèi)。與傳統(tǒng)全局辨識(shí)結(jié)果相比,待辨識(shí)參數(shù)的平均相對(duì)誤差明顯降低,精度更高。

表4 不同辨識(shí)方法誤差對(duì)比

3.2 不同迎角區(qū)間長(zhǎng)度的離線辨識(shí)時(shí)間對(duì)比

為研究分區(qū)數(shù)量與辨識(shí)時(shí)間的關(guān)系,本文分別進(jìn)行5°迎角區(qū)間、3°迎角區(qū)間、2°迎角區(qū)間以及傳統(tǒng)全局不分區(qū)仿真實(shí)驗(yàn),記錄辨識(shí)時(shí)間如表5所示。

表5 不同迎角區(qū)間長(zhǎng)度的離線辨識(shí)時(shí)間

根據(jù)表5能得到如下結(jié)論:迎角區(qū)間越小,分區(qū)的數(shù)量越多,離線辨識(shí)的時(shí)間越短。由此有以下建議:若迎角數(shù)據(jù)充分且對(duì)精度要求不高,可以適量減小迎角區(qū)間長(zhǎng)度,增大分區(qū)的數(shù)量,由此可減小程序運(yùn)行時(shí)間,提高計(jì)算效率。

4 結(jié) 論

本文基于遞推最小二乘法,圍繞局部線化代替整體非線性的思想,采用按回歸樹(shù)迎角分區(qū)的方法,對(duì)飛行器動(dòng)力學(xué)模型中未知參數(shù)的時(shí)域辨識(shí)方法進(jìn)行了研究。

首先,根據(jù)飛行器的縱向動(dòng)力學(xué)方程,建立了適用于仿真驗(yàn)證的觀測(cè)模型。之后,根據(jù)時(shí)域局部平滑法得到觀測(cè)模型所需要的觀測(cè)數(shù)據(jù)。在此基礎(chǔ)上,按照迎角分區(qū)方法,將所得數(shù)據(jù)劃分多個(gè)區(qū)域,區(qū)域間是相互獨(dú)立的,且能夠連續(xù)更新。并在每個(gè)不同的區(qū)域內(nèi)采用遞推最小二乘法,分別進(jìn)行時(shí)域辨識(shí)。最后,將每個(gè)區(qū)域的辨識(shí)結(jié)果利用加權(quán)函數(shù)整合,重構(gòu)俯仰力矩系數(shù)。

仿真過(guò)程和結(jié)果表明,本文基于遞推最小二乘法的時(shí)域分區(qū)辨識(shí)方法,可以較為準(zhǔn)確地給出局部氣動(dòng)參數(shù)的估計(jì)值,具有辨識(shí)精度高、計(jì)算效率高的特點(diǎn),且所有的參數(shù)辨識(shí)誤差小于3%,有利于了解全局氣動(dòng)模型中的局部氣動(dòng)特性。另外,通過(guò)比較不同迎角區(qū)間長(zhǎng)度,得到“可適量減小迎角區(qū)間長(zhǎng)度,增大分區(qū)數(shù)量,提高計(jì)算效率”的結(jié)論。然而,如果各區(qū)間迎角變化范圍取值較小,容易導(dǎo)致各區(qū)間仿真數(shù)據(jù)缺少足夠用于辨識(shí)的有效信息。進(jìn)而遞推初值計(jì)算誤差大,影響辨識(shí)精度,具有一定的局限性。因此,需要提供大量且充分激勵(lì)的仿真數(shù)據(jù)。