基于迎角分區的全局飛行器氣動參數辨識方法

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

(1. 大連理工大學航空航天學院, 遼寧 大連 116024; 2. 沈陽飛機設計研究所, 遼寧 沈陽 110035)

0 引 言

近年來,隨著系統辨識理論的研究,飛行器氣動參數辨識在飛行試驗中的應用得到了飛躍發展,并在飛行器設計中起著越來越重要的作用,采用參數辨識方法來確定飛機參數,能夠準確并迅速地將真實氣動特性從試驗結果中分離出來,對縮短數據處理時間和減少飛行試驗周期等具有重要應用價值[1]。

在大迎角機動時,由于機身機翼表面的氣流經歷了附著流-旋渦流-渦破裂的歷程,這才使得飛機氣動特性呈現強非線性等特征。傳統的線性辨識方法對大迎角機動非線性氣動參數辨識精度較低,因此線性辨識模型不適用于大迎角機動情況。

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

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

1 飛行器參數辨識模型

在參數辨識之前需要建立辨識模型,由于飛行器參數之間存在強耦合,導致辨識難度大。本文將飛行器的運動分解為縱向運動和橫側向運動,以縱向運動模型中俯仰力矩系數、升阻力系數為例,建立動力學參數辨識模型[14]。

1.1 縱向運動模型

把飛行器運動分解為縱向運動和橫側向運動,并以兩組相互獨立的運動方程組來描述[15]。通過對飛行器縱向受力進行分析,得到如下微分方程組[16]:

(1)

式中:m為飛行器質量;P為發動機推力;X為空氣阻力;Y為升力;ωz、Mz為俯仰角速率、俯仰力矩;x、y為x軸、y軸方向上的位移;Jz為轉動慣量;α、?、θ為迎角、俯仰角和軌跡傾角。

1.2 觀測模型

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

(2)

(3)

空氣動力學的整體模型的理想形式是同時具有估測的模型參數值和相關不確定性的數學模型結構,可以將無量綱的空氣動力學中的力和力矩系數與可以測量的飛機狀態和控制聯系起來。所有的全局辨識都是基于時域的方程誤差最小二乘法來進行建模。在此方法中,因變量是無量綱力或力矩系數中的某一項,它是通過利用解釋變量(無量綱的飛機狀態量和控制面撓度)計算出非線性模型項的擴展來進行建模的[18]。

因此,可以使用特定模型結構來表達俯仰力矩系數的方程誤差最小二乘問題。設置建模函數如下:

(4)

式(4)可以寫成:

z=Hθ+v

(5)

其中,z代表觀測模型計算的觀測值:

(6)

H代表設置的建模函數:

(7)

θ代表辨識模型:

(8)

v代表殘差:

(9)

2 動力學參數辨識方法

2.1 基于軸正交決策樹的迎角自動分區

回歸樹是一種具有邏輯性的分類方法,可以為局部建模網絡劃分權重變量,并分出新的單元結構,該過程將解釋變量歸入不同的區域,來表征解釋變量的非線性特征。分區方法通常采用二進制拆分法則,圖1所示二進制軸正交決策樹網絡結構,在具有兩個權重變量的情況下,在每個級別的單元格中再拆分出一個新單元格。被劃分的單元稱為母單元,隨后拆分出的單元稱為子單元[19]。

圖1 二進制軸正交決策樹網絡圖

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

(10)

可迭代的二進制軸正交回歸樹算法已成功用于各種非線性靜態和動態建模應用程序中。

假設:

(11)

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

(12)

2.2 飛行器在線辨識激勵信號設計

實時氣動建模時需要迅速地收集到所有剛體自由度的數據信息,而在一系列飛行條件下,將自動正交優化的多正弦擾動施加在控制舵面上可以很有效地收集相關氣動數據[23]。使用擾動輸入來激勵飛行過程,該輸入通過將設計好的激勵信號與控制指令中的執行器命令相加,以此作為飛行器的舵面指令。

設定應用在第j個控制面的擾動輸入為uj,其為具有單個相移的正弦波諧波φk的總和[24]。

(13)

式中:M是可用的諧波相關頻率的總數;T是激勵的時間長度;Ak是正弦波分量的振幅;t是時間分量;nj個輸入中的每一個都是從M個諧波正弦波中選擇具有頻率ωk=2πk/T的分量,k=1,2,…,M。ωm=2πM/T代表激勵輸入的頻帶上限。區間[ω1-ωm]代表了指定期望飛機動力學的頻率范圍。

為了實現均勻的功率分配,Ak被設置為

(14)

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

2.3 時域遞推最小二乘辨識方法

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

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

(15)

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

(16)

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

(17)

z=Hθ4×N

(18)

由式(18)可知:

(19)

2.4 基于加權函數的數據平滑方法

當每個局部模型辨識氣動參數以后,需要對全局范圍內氣動參數進行數據平滑整合,該部分基于加權函數對局部模型的臨界數據進行數據平滑[27]。具體步驟如下:

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

(20)

Ck,j是第j個權重變量在該單元格內的中心點的高斯函數,計算公式如下[28]:

(21)

式中:μ代表在該單元內的該分區變量的均值;σ2代表方差。值得注意的是,由于權重變量是隨著數據點的不同而不斷變化的,而在每段區域中心處求該區域中心點的高斯中心值是固定的,因此根據上述公式得當某點離某一區域的高斯中心越遠時,該區域在這點處有效性越小[29]。所以,比如對第一個區域內的某一數據點來說,主要是第一個區域在該點有效性比較大,離其他區域越遠,則其在該點有效性越小。

(2) 根據上述公式得到各個區域的有效性后,通過式(22)來計算第k個單元的權重[30]。

(22)

(3) 利用計算出的各個單元格的權重,進行全局平滑建模。對每個點的整體建模,采用各個單元內區域建模的加權加和得到[31]:

(23)

首先,平滑的依據有兩個:① 在斷點處的值相等;② 在斷點處導數相等。

(24)

則在臨界點的平滑計算值為

(25)

證畢

以上平滑的兩點依據均滿足,因此可以認為該曲線已經平滑。通過上述過程,可以將各個分區的建模結果進行平滑加權加和,從而得到一個平滑后的有效的整體模型。

3 仿真實例

為驗證本文方法的有效性,針對飛行高度7 500 m、Ma0.4的某飛行器飛行數據,建立縱向運動模型、辨識觀測模型以及遞推最小二乘辨識模型,總仿真時間為ts=25 s,仿真步長為tk=0.01 s。飛行器的相關參數如表1所示。飛行器的飛行狀態數據如圖2所示。

表1 飛行器相關參數

圖2 大迎角機動仿真數據

3.1 仿真計算結果及分析

圖3展示仿真流程示意圖,仿真流程主要包括3個方面:基于激勵信號的數據收集,建立觀測模型、空氣動力學模型、最小二乘辨識模型,基于加權函數的數據平滑處理。

圖3 仿真流程示意圖

本文將激勵信號設計為具有特定諧波頻率,優化相移和指定功率分布之后的正弦波之和[32],如圖4所示。

圖4 正弦激勵信號

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

(26)

選擇特定的分量頻率以覆蓋飛行器姿態運動包含的頻帶[33]。另外,將激勵信號設計為多元輸入,保證其在時域和頻域中相互正交,以使其在短時間內所有軸上都具有較高的數據信息含量[34]。

在充分獲取數據的基礎上,建立俯仰力矩系數、升力系數的觀測模型如式(27)。取迎角變化范圍為5°,基于回歸樹迎角自動分區結果如圖5所示。

圖5 回歸樹迎角自動分區結果

(27)

根據每個分區的辨識情況,總結了計算俯仰力矩系數、升力系數各項參數的辨識結果。由圖6和圖7可知,氣動導數在5 s左右已經收斂,且收斂速度較快。氣動導數的辨識精度分別為2.1%、2.9%、1.6%,詳見表2。

表2 俯仰力矩系數辨識精度

圖6 俯仰靜穩定系數和升降舵效率辨識結果

由圖8和圖9可見,升力系數的氣動導數辨識收斂較快,辨識精度分別為0.6%、0.5%,詳見表3。

表3 升力系數辨識精度

圖8 迎角系數和升降舵系數辨識結果

圖9 升力系數辨識結果

為充分證明迎角分區時域辨識精度較傳統全局辨識精度高,下文闡述迎角分區與傳統全局的辨識結果對比,如圖10所示。

圖10 俯仰靜穩定系數和升降舵效率辨識結果

由表4可知,利用迎角分區辨識的方法,辨識俯仰靜穩定系數、升降舵操縱效率以及俯仰阻尼力矩系數,其相對誤差均在3%以內。與傳統全局辨識結果相比,待辨識參數的平均相對誤差明顯降低,精度更高。

表4 不同辨識方法誤差對比

3.2 不同迎角區間長度的離線辨識時間對比

為研究分區數量與辨識時間的關系,本文分別進行5°迎角區間、3°迎角區間、2°迎角區間以及傳統全局不分區仿真實驗,記錄辨識時間如表5所示。

表5 不同迎角區間長度的離線辨識時間

根據表5能得到如下結論:迎角區間越小,分區的數量越多,離線辨識的時間越短。由此有以下建議:若迎角數據充分且對精度要求不高,可以適量減小迎角區間長度,增大分區的數量,由此可減小程序運行時間,提高計算效率。

4 結 論

本文基于遞推最小二乘法,圍繞局部線化代替整體非線性的思想,采用按回歸樹迎角分區的方法,對飛行器動力學模型中未知參數的時域辨識方法進行了研究。

首先,根據飛行器的縱向動力學方程,建立了適用于仿真驗證的觀測模型。之后,根據時域局部平滑法得到觀測模型所需要的觀測數據。在此基礎上,按照迎角分區方法,將所得數據劃分多個區域,區域間是相互獨立的,且能夠連續更新。并在每個不同的區域內采用遞推最小二乘法,分別進行時域辨識。最后,將每個區域的辨識結果利用加權函數整合,重構俯仰力矩系數。

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