基于機(jī)器學(xué)習(xí)的氣動(dòng)參數(shù)智能修正方法

張家銘 鐘鴻豪 白文艷 孫 友 曹玉騰

北京航天自動(dòng)控制研究所，北京 100854

0 引言

由于飛行器風(fēng)洞試驗(yàn)的局限性以及理論計(jì)算受限于不完善的理論研究，所以應(yīng)用氣動(dòng)參數(shù)辨識技術(shù)從飛行試驗(yàn)數(shù)據(jù)中獲得飛行器模型的各個(gè)氣動(dòng)參數(shù)，已成為飛行器設(shè)計(jì)和研究中不可分割的一部分[1]。

在氣動(dòng)參數(shù)辨識相關(guān)研究方面，國外的起步很早，主要包括美國、俄羅斯(前蘇聯(lián))、歐洲和以色列等傳統(tǒng)航空強(qiáng)國。其中美國從上世紀(jì)中旬開始，研究了從低空亞聲速、高空超聲速到臨近空間高超聲速的、不同構(gòu)型的飛行器空氣動(dòng)力學(xué)現(xiàn)象和氣動(dòng)力參數(shù)，形成了系統(tǒng)全面的研究成果。早在70年代，經(jīng)典的系統(tǒng)辨識方法如極大似然法和最小二乘法剛剛被提出的時(shí)候[2],就在氣動(dòng)數(shù)據(jù)的辨識研究中發(fā)揮了關(guān)鍵作用。例如文獻(xiàn)[3]中，Klein系統(tǒng)研究了最小二乘法在飛行器氣動(dòng)參數(shù)辨識中的應(yīng)用，分析了辨識算法、估計(jì)的準(zhǔn)確性、可辨識性以及模型驗(yàn)證等系統(tǒng)辨識的基本問題，并結(jié)合了多種機(jī)型的例子進(jìn)行了說明。文獻(xiàn)[4-5]介紹了使用極大似然法、最小二乘法以及最小二乘法的變形來估計(jì)氣動(dòng)參數(shù)。在文獻(xiàn)[6]中，德國航空航天中心DLR的Ravindra Jategaonkar博士系統(tǒng)總結(jié)了二十世紀(jì)的飛行器系統(tǒng)辨識的發(fā)展歷程和所應(yīng)用的各種方法，例如方程誤差法、輸出誤差法、濾波器誤差法以及頻域方法、基于濾波的方法、神經(jīng)網(wǎng)絡(luò)方法等等。文獻(xiàn)[7]介紹了飛行器參數(shù)估計(jì)中的輸出誤差方法。文獻(xiàn)[8-9]研究了應(yīng)用在氣動(dòng)力參數(shù)辨識的濾波方法，包括卡爾曼濾波、拓展卡爾曼濾波以及無跡卡爾曼濾波方法。文獻(xiàn)[10]總結(jié)了迭代計(jì)算的幾種系統(tǒng)辨識方法。文獻(xiàn)[11]介紹了一種頻域上的方法。同時(shí)，歐美以及AIAA協(xié)會(huì)也建立了適用于氣動(dòng)參數(shù)辨識的集成工具集和軟件方法庫，并且收集了大量的理論計(jì)算、風(fēng)洞試驗(yàn)和試飛的一手?jǐn)?shù)據(jù)，形成了較為全面的數(shù)據(jù)集。ESTIMA就是這樣一款著名的氣動(dòng)辨識集成軟件工具[12]。在數(shù)據(jù)集方面，文獻(xiàn)[13]總結(jié)了各國無翼的載人航天不同氣動(dòng)外形的返回艙、各種類型的再入彈頭、有翼的航天飛機(jī)、翼身融合體和升力體等不同型號飛行器的氣動(dòng)力數(shù)據(jù)。

實(shí)際工程中的飛行試驗(yàn)數(shù)據(jù)往往包含氣動(dòng)參數(shù)的重要信息，但卻未被好好利用。本文將利用氣動(dòng)參數(shù)辨識方法，充分挖掘已有飛行試驗(yàn)數(shù)據(jù)中的氣動(dòng)參數(shù)信息，對氣動(dòng)表進(jìn)行修正。隨著飛行數(shù)據(jù)的增多，準(zhǔn)確性得到不斷提高。然而傳統(tǒng)的氣動(dòng)參數(shù)辨識方法如最小二乘法、極大似然法、拓展卡爾曼濾波法等多基于線性回歸模型，或者對非線性模型進(jìn)行線性化，略去高階非線性項(xiàng)，因此傳統(tǒng)辨識方法多用于線性化較強(qiáng)的模型，一旦模型非線性較強(qiáng)或時(shí)變性較快時(shí)，傳統(tǒng)辨識方法的穩(wěn)定性將難以保證[14-15]。另外，傳統(tǒng)的氣動(dòng)參數(shù)辨識方法直接對氣動(dòng)參數(shù)和飛行狀態(tài)的關(guān)系進(jìn)行建模。對于一些復(fù)雜的飛行器構(gòu)型，由于氣動(dòng)參數(shù)的模型階次較高，其本身繪制出的曲線可能較為復(fù)雜或存在劇烈波動(dòng)及跳變的存在，所以會(huì)嚴(yán)重影響氣動(dòng)參數(shù)的辨識結(jié)果[16]。而BP神經(jīng)網(wǎng)絡(luò)具有逼近任意非線性函數(shù)的功能，其學(xué)習(xí)算法屬于全局逼近算法，具有較強(qiáng)的泛化能力，可以用于氣動(dòng)參數(shù)的辨識。

在德國RLV Phoenix飛行器研究[17]中發(fā)現(xiàn)，對真實(shí)氣動(dòng)數(shù)據(jù)與標(biāo)稱氣動(dòng)數(shù)據(jù)的偏差建模，并將其用于氣動(dòng)參數(shù)的修正，將簡化氣動(dòng)參數(shù)修正過程，避免上述問題的發(fā)生。

本文提出的基于機(jī)器學(xué)習(xí)的氣動(dòng)參數(shù)智能修正方法，利用BP神經(jīng)網(wǎng)絡(luò)對氣動(dòng)數(shù)據(jù)的偏差模型進(jìn)行擬合，根據(jù)擬合后的神經(jīng)網(wǎng)絡(luò)以及標(biāo)稱氣動(dòng)表，即可得到修正后的氣動(dòng)表，并且可以進(jìn)一步得到任意狀態(tài)的氣動(dòng)參數(shù)，實(shí)現(xiàn)離線對網(wǎng)絡(luò)模型進(jìn)行訓(xùn)練，并利用歷史飛行數(shù)據(jù)進(jìn)行模型修正的目的。在本文的第四部分，針對公開的CAV(Common Aero Vehicle)系統(tǒng)的氣動(dòng)參數(shù)進(jìn)行了辨識，驗(yàn)證了方法的有效性。

1 問題描述

使用公開的CAV系統(tǒng)的標(biāo)稱飛行氣動(dòng)參數(shù)表和縱向飛行的彈道[18]作為研究對象。本項(xiàng)目使用的CAV飛行器為類似美國HTV-2高超聲速飛行器的升力體氣動(dòng)外形，外形具體見圖1。

此CAV飛行器只在尾部有兩個(gè)水平可控舵面，依靠同步偏轉(zhuǎn)實(shí)現(xiàn)升降舵偏，利用差動(dòng)偏轉(zhuǎn)進(jìn)行飛行器的滾轉(zhuǎn)控制，并無方向舵。其縱向運(yùn)動(dòng)模型為[19]：

(1)

式中：m為飛行器的質(zhì)量，m為地球引力常量，Jz為飛行器繞彈體系 z 軸的轉(zhuǎn)動(dòng)慣量，V為飛行器的速度，H為飛行高度，α、q、wz分別為攻角、航跡角和俯仰角速率，T為發(fā)動(dòng)機(jī)推力；L、D分別為飛行器所受升力和阻力、推力，Mz為俯仰力矩，且有：

(2)

式中：CL、CD和Cm分別為升力系數(shù)、阻力系數(shù)和俯仰力矩系數(shù)及偏航力矩系數(shù)，ρ為該垂直高度下的大氣密度，S為飛行器的橫截面積，bA飛行器平均氣動(dòng)弦長，S和bA可以認(rèn)為是常數(shù)。

氣動(dòng)系數(shù)函數(shù)可表示為：

[CL,CD,Cm]=f(H,Ma,α,δe)

(3)

式中：Ma為馬赫數(shù)，δe為升降舵偏。

針對彈道數(shù)據(jù)，求取氣動(dòng)參數(shù)CL、CD和Cm。

為了求取氣動(dòng)參數(shù)，對真實(shí)氣動(dòng)表和標(biāo)稱氣動(dòng)表間的偏差進(jìn)行建模：

C真實(shí)=C標(biāo)稱+E

(4)

式中：E代表真實(shí)氣動(dòng)表和標(biāo)稱氣動(dòng)表間的偏差。對偏差進(jìn)行建模時(shí)，偏差模型的階數(shù)可由經(jīng)驗(yàn)確定或在相關(guān)文獻(xiàn)中查找。

可以定義修正后的氣動(dòng)數(shù)據(jù)表與標(biāo)稱氣動(dòng)表之間的偏差函數(shù)為：

Δ(H,Ma,α,δe)=

f修正(H,Ma,α,δe)-f標(biāo)稱(H,Ma,α,δe)

(5)

偏差模型擬合修正法就是將擬合氣動(dòng)數(shù)據(jù)表，改為擬合這個(gè)偏差函數(shù)。我們選取偏差模型為：

Δ(H,Ma,α,δe)=a1+a2H+a3Ma+a4α+a5δe

(6)

修正后的氣動(dòng)數(shù)據(jù)表為標(biāo)稱氣動(dòng)表加上使用偏差模型計(jì)算出的偏差，即：

f修正(H,Ma,α,δe)=

(7)

2 基于機(jī)器學(xué)習(xí)的氣動(dòng)參數(shù)智能修正方法

2.1 BP神經(jīng)網(wǎng)絡(luò)

BP神經(jīng)網(wǎng)絡(luò)法是一種被廣泛使用的經(jīng)典非線性系統(tǒng)辨識方法也是最常用最流行最成熟的人工神經(jīng)網(wǎng)絡(luò)之一。在網(wǎng)絡(luò)中，各神經(jīng)元接受上一級輸入，輸出到下一級，在理論上可以逼近任意連續(xù)非線性函數(shù)。BP網(wǎng)絡(luò)是一種2層或2層以上的階層型神經(jīng)網(wǎng)絡(luò)，即輸入層、隱含層(也稱中間層)和輸出層，其原理和過程如圖2所示[19]。

圖2 BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)

BP算法的主要思想是把學(xué)習(xí)過程分為信號的正向傳播與誤差的反向傳播2個(gè)階段。在正向傳播階段，輸入信息從輸入層經(jīng)隱含層傳向輸出層，在輸出端產(chǎn)生輸出信號。在信號的向前傳遞過程中網(wǎng)絡(luò)的權(quán)值固定不變，每一層神經(jīng)元的狀態(tài)只影響下一層神經(jīng)元的狀態(tài)。如果在輸出層不能得到期望的輸出，則轉(zhuǎn)入誤差信號反向傳播。在反向傳播階段，未能滿足精度要求的誤差信號由輸出端開始，以某種方式逐層向前傳播，并將誤差分?jǐn)偨o各層的所有單元，依據(jù)誤差信號動(dòng)態(tài)地調(diào)整各單元層的連接權(quán)重。通過周而復(fù)始的正向傳播與反向調(diào)節(jié)，神經(jīng)元間的權(quán)值得到不斷修正。當(dāng)輸出信號的誤差滿足精度要求時(shí)，停止學(xué)習(xí)。

2.2 基于BP神經(jīng)網(wǎng)絡(luò)的氣動(dòng)參數(shù)修正算法

氣動(dòng)表的非線性函數(shù)關(guān)系可能非常復(fù)雜，使用神經(jīng)網(wǎng)絡(luò)模型對其進(jìn)行擬合可以得到較好的效果。

本文的修正算法思路為：為根據(jù)式(5)計(jì)算得到偏差函數(shù)，利用神經(jīng)網(wǎng)絡(luò)對偏差函數(shù)進(jìn)行擬合，偏差模型為式(6)的形式，然后根據(jù)擬合得到的偏差模型神經(jīng)網(wǎng)絡(luò)以及標(biāo)稱氣動(dòng)表，即可得到修正后的氣動(dòng)表。

該算法的示意圖如圖3 所示。

圖3 基于BP神經(jīng)網(wǎng)絡(luò)的氣動(dòng)參數(shù)修正算法示意圖

基于BP神經(jīng)網(wǎng)絡(luò)的氣動(dòng)表修正算法的設(shè)計(jì)流程如下：

步驟1：偏差函數(shù)的構(gòu)建。

利用彈道氣動(dòng)數(shù)據(jù)以及標(biāo)稱氣動(dòng)表數(shù)據(jù)，根據(jù)式(5)構(gòu)建偏差函數(shù)。

步驟2：BP神經(jīng)網(wǎng)絡(luò)構(gòu)建。

根據(jù)擬合非線性函數(shù)特點(diǎn)確定BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)。本文研究對象的氣動(dòng)表非線性函數(shù)有4個(gè)輸入?yún)?shù)：高度H，馬赫數(shù)Ma，攻角α，升降舵偏δe；3個(gè)輸出參數(shù)：升力系數(shù)CL，阻力系數(shù)CD，俯仰力矩系數(shù)Cm。

對于神經(jīng)網(wǎng)絡(luò)來說，隱含層的節(jié)點(diǎn)數(shù)過少，網(wǎng)絡(luò)會(huì)欠擬合，即對數(shù)據(jù)的擬合程度不高，網(wǎng)絡(luò)沒有很好地捕捉到數(shù)據(jù)特征。節(jié)點(diǎn)數(shù)過多時(shí)，網(wǎng)絡(luò)會(huì)過擬合。因此，為了確定神經(jīng)網(wǎng)絡(luò)隱含層的層數(shù)和節(jié)點(diǎn)數(shù)，做出如下的規(guī)定：

在單隱含層中，單層測試節(jié)點(diǎn)數(shù)的上限為75，在雙隱含層中，遵循以下原則：

1) 第一隱含層節(jié)點(diǎn)數(shù)不小于第二隱含層節(jié)點(diǎn)數(shù)。

2) 第一隱含層和第二隱含層的節(jié)點(diǎn)數(shù)之和上限為60。

采用等間隔取樣法，從節(jié)點(diǎn)數(shù)為5開始，以5為單位，進(jìn)行單隱含層網(wǎng)絡(luò)和雙隱含層網(wǎng)絡(luò)的測試?；陔[含層節(jié)點(diǎn)信息約束條件獲取不同隱含層層數(shù)以及不同節(jié)點(diǎn)數(shù)下的均方誤差(Mean squared error, MSE)，并將其作為神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)的選取指標(biāo)，選擇網(wǎng)絡(luò)的超參數(shù)。MSE的值越大，表明網(wǎng)絡(luò)的預(yù)測效果越差，因此可以將其當(dāng)作網(wǎng)絡(luò)結(jié)構(gòu)的選取依據(jù)。

之后，根據(jù)不同氣動(dòng)參數(shù)所對應(yīng)的每組結(jié)構(gòu)，選擇出不同氣動(dòng)參數(shù)所對應(yīng)的MSE較小的幾組神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)作為最優(yōu)神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)集合，縮小最優(yōu)神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)的選取范圍。

根據(jù)訓(xùn)練結(jié)果的MSE，對于氣動(dòng)參數(shù)CL，選擇雙層網(wǎng)絡(luò)結(jié)構(gòu)(20,15)作為最優(yōu)網(wǎng)絡(luò)結(jié)構(gòu)；CD和Cm均選擇雙層網(wǎng)絡(luò)結(jié)構(gòu)(30,5)作為最優(yōu)網(wǎng)絡(luò)結(jié)構(gòu)。(注：這些結(jié)構(gòu)僅為取用等間隔采樣法下的最優(yōu)網(wǎng)絡(luò)配置，并非在所有可能的網(wǎng)絡(luò)下的最優(yōu)網(wǎng)絡(luò)配置，縮小間隔可以獲取更準(zhǔn)確的配置信息，但同時(shí)，計(jì)算量也會(huì)增大。)

步驟3：偏差模型神經(jīng)網(wǎng)絡(luò)訓(xùn)練。

利用多組彈道氣動(dòng)數(shù)據(jù)，對所構(gòu)建的BP神經(jīng)網(wǎng)絡(luò)進(jìn)行訓(xùn)練。為了避免訓(xùn)練過度擬合的問題，使用k-fold交叉驗(yàn)證的技術(shù)手段。把訓(xùn)練數(shù)據(jù)隨機(jī)分成一定比例的訓(xùn)練集、驗(yàn)證集和測試集，采用交叉驗(yàn)證的方式訓(xùn)練k次，每次隨機(jī)選數(shù)據(jù)進(jìn)行訓(xùn)練，剩余的驗(yàn)證集作為驗(yàn)證，邊訓(xùn)練，邊驗(yàn)證，從而找出最優(yōu)解，得到偏差模型神經(jīng)網(wǎng)絡(luò)。取k=7。

圖4 k-fold交叉驗(yàn)證示意圖

步驟4：修正后氣動(dòng)表構(gòu)建

根據(jù)標(biāo)稱氣動(dòng)表以及擬合得到的偏差模型神經(jīng)網(wǎng)絡(luò)，構(gòu)建修正后的氣動(dòng)表。

根據(jù)得到的修正后的氣動(dòng)表，可以進(jìn)一步計(jì)算出氣動(dòng)參數(shù)。

3 仿真校驗(yàn)及分析

為了驗(yàn)證所提方法的有效性，針對CAV系統(tǒng)的標(biāo)稱飛行氣動(dòng)參數(shù)表和縱向飛行的彈道數(shù)據(jù)，使用本文的修正方法采用MATLAB R2017b中的神經(jīng)網(wǎng)絡(luò)工具箱Neural Network Toolbox 11.0中的Neural Net Fitting對其氣動(dòng)表進(jìn)行了修正。試驗(yàn)中，選取風(fēng)洞試驗(yàn)得到的標(biāo)稱氣動(dòng)參數(shù)，利用量測彈道數(shù)據(jù)求解得到的氣動(dòng)參數(shù)作為真實(shí)氣動(dòng)參數(shù)，得到氣動(dòng)偏差數(shù)據(jù)，利用插值生成1890組偏差氣動(dòng)數(shù)據(jù)，用于訓(xùn)練修正神經(jīng)網(wǎng)絡(luò)。

本文中BP神經(jīng)網(wǎng)絡(luò)隱含層的非線性環(huán)節(jié)選用“Sigmoid”非線性函數(shù)，輸出層的神經(jīng)元類型為線性輸出神經(jīng)元。本算法采用k-fold交叉驗(yàn)證法，選取70%的訓(xùn)練數(shù)據(jù)作為訓(xùn)練集，驗(yàn)證集和測試集各占15%的數(shù)據(jù)量。在訓(xùn)練前對訓(xùn)練數(shù)據(jù)進(jìn)行歸一化處理。反向傳播的訓(xùn)練算法采用Levenberg-Marquardt反向傳播算法。損失函數(shù)為均方誤差(MSE)準(zhǔn)則。

3.1 氣動(dòng)參數(shù)修正相對誤差率

(8)

基于BP神經(jīng)網(wǎng)絡(luò)的氣動(dòng)表修正方法的氣動(dòng)參數(shù)修正結(jié)果見表1。

表1 算法的修正結(jié)果

其中俯仰力矩相對誤差率曲線如圖5所示，橫坐標(biāo)為俯仰力矩大小。

圖5 俯仰力矩系數(shù)相對誤差率隨俯仰力矩系數(shù)變化

由表1可以看到，對于升力系數(shù)來說，該算法的最大誤差率和平均誤差率均小于1%；對于阻力系數(shù)來說，該算法的最大誤差率和平均誤差率均小于0.2%，取得了優(yōu)良的修正效果。而對俯仰力矩系數(shù)來說，其最大誤差率和平均誤差率稍大，分別達(dá)到了2.955%和2.333%。根據(jù)圖5分析可知，俯仰力矩系數(shù)本身的數(shù)量級較小，個(gè)別數(shù)據(jù)近似為0，所以誤差率較大，但平均誤差率均在3%以內(nèi)，仍取得了良好的修正效果。

3.2 修正前后飛行結(jié)果對比

基于BP神經(jīng)網(wǎng)絡(luò)的氣動(dòng)表修正方法可以得到修正后的氣動(dòng)表，利用修正后的氣動(dòng)表，可以進(jìn)一步得到修正后的彈道、俯仰角、俯仰角速度以及攻角曲線，將其與修正前以及真實(shí)的曲線進(jìn)行對比，繪制如圖6～9所示。

圖6 彈道曲線對比圖

圖7 俯仰角曲線對比圖

圖8 俯仰角速度曲線對比圖

圖9 攻角隨時(shí)間變化曲線對比圖

根據(jù)圖6可知，修正后的彈道飛行距離誤差由16%降到小于3%；同樣的，根據(jù)圖7～9，可以看到，經(jīng)過本論文的方法修正，修正后曲線與真實(shí)曲線的誤差較小，證明了本論文方法的有效性。

3.3 最小二乘法與神經(jīng)網(wǎng)絡(luò)結(jié)果對比

將傳統(tǒng)方法最小二乘方法與人工神經(jīng)網(wǎng)絡(luò)方法修正結(jié)果比較，以氣動(dòng)參數(shù)CL為例，假設(shè)CL=CL0+CLαα+CLδδ+CLωzωz，待辨識參數(shù)為θ={CLα,CLδ,CLωz,CL0}，設(shè)：

(9)

用最小二乘辨識算法θ=(φTφ)-1φTCL得出：θ={0.116891747195653,-0.179589910444440,

-0.196493513367346,0.287336006156336}.代入待修正的氣動(dòng)參數(shù)，得到其修正平均相對誤差為3.475%，遠(yuǎn)大于人工神經(jīng)網(wǎng)絡(luò)修正平均相對誤差0.712%。

3.4 數(shù)據(jù)量對氣動(dòng)參數(shù)修正誤差的影響分析

彈道數(shù)據(jù)對氣動(dòng)參數(shù)修正相對誤差的影響如圖10所示。

圖10 彈道數(shù)量(數(shù)據(jù)量)對修正CL平均相對誤差的影響

從圖中可以看出，以氣動(dòng)參數(shù)CL為例，隨著彈道數(shù)量的增多，修正誤差越小。即數(shù)據(jù)量的增大，修正CL的平均相對誤差是逐漸減少的。說明隨著飛行數(shù)據(jù)的增多，氣動(dòng)數(shù)據(jù)修正得更精確，往后的發(fā)射就更精準(zhǔn)。

4 結(jié)論

提出了一種基于機(jī)器學(xué)習(xí)的氣動(dòng)參數(shù)智能修正方法，利用彈道數(shù)據(jù)和標(biāo)稱數(shù)據(jù)，擬合出較為準(zhǔn)確的偏差模型網(wǎng)絡(luò)，結(jié)合標(biāo)稱氣動(dòng)表得到修正后的氣動(dòng)參數(shù)。另外，結(jié)合公開的CAV縱向氣動(dòng)數(shù)據(jù)，計(jì)算了基于機(jī)器學(xué)習(xí)的氣動(dòng)參數(shù)智能修正方法的修正結(jié)果，驗(yàn)證了方法的有效性。該方法充分挖掘已有飛行試驗(yàn)數(shù)據(jù)中的氣動(dòng)參數(shù)信息，對氣動(dòng)參數(shù)表進(jìn)行修正，得到較為準(zhǔn)確的氣動(dòng)模型，并隨著彈道數(shù)據(jù)量增多，氣動(dòng)修正精度越高，發(fā)射一發(fā)比一發(fā)更準(zhǔn)確，在工程上具有實(shí)用價(jià)值。