高 鑫，連 峰，孔軼男，陳 功，張光華，王文正

(1.西安交通大學(xué)，西安 710049；2.中國(guó)空氣動(dòng)力研究與發(fā)展中心，綿陽(yáng) 621000；3.電子科技大學(xué) 航空航天學(xué)院，成都 611731)

0 引言

飛行器氣動(dòng)力系數(shù)建模是研究飛行器控制系統(tǒng)設(shè)計(jì)、飛行控制、動(dòng)態(tài)分析以及性能評(píng)估的基礎(chǔ)。隨著現(xiàn)代航空航天技術(shù)的不斷發(fā)展，現(xiàn)代戰(zhàn)機(jī)可以在更廣的包絡(luò)線內(nèi)機(jī)動(dòng)。飛行器在飛行過(guò)程中氣動(dòng)環(huán)境復(fù)雜，氣動(dòng)系數(shù)的模型呈現(xiàn)出多變量、高度非線性的特點(diǎn)。

在氣動(dòng)力學(xué)建模領(lǐng)域，國(guó)內(nèi)外已經(jīng)逐漸形成了一套完整的體系，主要分為基于物理機(jī)理建模和基于數(shù)據(jù)建模兩個(gè)方向。基于物理機(jī)理建模的有氣動(dòng)導(dǎo)數(shù)模型、狀態(tài)空間模型、階躍響應(yīng)模型以及非線性微分方程模型[1-8]；這些模型的物理意義明確，適用于飛行器優(yōu)化設(shè)計(jì)與性能評(píng)估?；跀?shù)據(jù)建模主要是運(yùn)用各種智能學(xué)習(xí)算法進(jìn)行建模，例如神經(jīng)網(wǎng)絡(luò)模型、模糊邏輯模型[9-12]等；得益于這些智能算法的強(qiáng)大映射能力，基于數(shù)據(jù)建?？梢苑浅Ｓ行У氐玫椒蔷€性領(lǐng)域內(nèi)的模型。上述氣動(dòng)系數(shù)建模方法通常包含大量重復(fù)的地面測(cè)試試驗(yàn)，如風(fēng)洞試驗(yàn)與計(jì)算機(jī)仿真試驗(yàn)；此外還需進(jìn)行多次精確的飛行試驗(yàn)，利用飛行試驗(yàn)結(jié)果對(duì)地面建模結(jié)果進(jìn)行修正，才能得到精確的全局氣動(dòng)模型。

為了簡(jiǎn)化建模流程，美國(guó)航空航天局NASA 提出“邊學(xué)邊飛”概念（Learn-To-Fly），介紹了一種新的建模方法，以最少的地面測(cè)試、人機(jī)交互和分析時(shí)間，實(shí)現(xiàn)飛行器自主建立實(shí)時(shí)的全局氣動(dòng)模型。相關(guān)研究已經(jīng)表明[13-16]，對(duì)于新型高效機(jī)動(dòng)動(dòng)作[17-18]而言，僅依靠飛行數(shù)據(jù)與先進(jìn)的非線性建模技術(shù)相結(jié)合，即可實(shí)時(shí)獲得目標(biāo)飛行器的全局非線性模型。基于飛行數(shù)據(jù)的實(shí)時(shí)建模是為了在滿足精度的條件下，得到一個(gè)盡可能精簡(jiǎn)的氣動(dòng)力系數(shù)模型，更有利于目標(biāo)飛行器快速分析模型并及時(shí)根據(jù)氣動(dòng)環(huán)境的變化指定控制策略。多元正交函數(shù)建模（multivariate orthogonal function modeling,MOFM）方法是一種基于數(shù)據(jù)的建模方法，其基于實(shí)時(shí)氣動(dòng)數(shù)據(jù)建立氣動(dòng)力和力矩系數(shù)關(guān)于攻角、側(cè)滑角、升降舵偏、方向舵偏以及副翼舵偏的非線性函數(shù)時(shí)，需要確定模型的結(jié)構(gòu)并對(duì)參數(shù)進(jìn)行估計(jì)。一些學(xué)者采用逐步回歸方法確定多項(xiàng)式結(jié)構(gòu)，但是回歸矩陣會(huì)受實(shí)時(shí)數(shù)據(jù)影響，出現(xiàn)病態(tài)矩陣問(wèn)題，無(wú)法準(zhǔn)確確定模型結(jié)構(gòu)。本文采用QR 分解矩陣并逐步迭代實(shí)時(shí)數(shù)據(jù)的方法，能夠有效避免回歸矩陣病態(tài)問(wèn)題；同時(shí)，系數(shù)矩陣維度不隨采樣數(shù)增加而升高，能滿足建模過(guò)程中時(shí)效性的要求。

本文使用F-16 風(fēng)洞試驗(yàn)仿真程序[19]生成建模需要的數(shù)據(jù)集，通過(guò)逐條讀取訓(xùn)練集中的氣動(dòng)數(shù)據(jù)模擬真實(shí)飛行過(guò)程中的實(shí)時(shí)采樣數(shù)據(jù)；采用多元正交函數(shù)進(jìn)行氣動(dòng)力學(xué)系數(shù)快速建模，并對(duì)比了該方法與徑向基函數(shù)（RBF 網(wǎng)絡(luò)）、BP 網(wǎng)絡(luò)對(duì)風(fēng)洞仿真數(shù)據(jù)的預(yù)測(cè)結(jié)果。研究結(jié)果表明，利用多元正交函數(shù)進(jìn)行氣動(dòng)力學(xué)系數(shù)建模能夠快速準(zhǔn)確地建立氣動(dòng)模型，為研究快速建模方法提供了新的思路與方法。

1 多元正交函數(shù)建模原理

基于實(shí)時(shí)飛行數(shù)據(jù)的多元正交函數(shù)非線性氣動(dòng)學(xué)系數(shù)建模需要解決兩個(gè)問(wèn)題：1）需要確定模型的結(jié)構(gòu)。這需要從候選函數(shù)池的大量非線性項(xiàng)中剔除對(duì)模型擬合效果影響較小的項(xiàng)，采用最小預(yù)測(cè)均方誤差（predicted squared error,PSE）中準(zhǔn)則可以滿足上述需求。2）需要根據(jù)實(shí)時(shí)數(shù)據(jù)辨識(shí)出結(jié)構(gòu)參數(shù)。傳統(tǒng)的基于最小二乘回歸進(jìn)行參數(shù)估計(jì)的方法會(huì)由于變量之間的線性相關(guān)性導(dǎo)致病態(tài)的回歸矩陣，造成參數(shù)估計(jì)不準(zhǔn)。本文將回歸矩陣QR 分解后在R矩陣上迭代數(shù)據(jù)，解決了上述問(wèn)題。

1.1 模型初始化

將辨識(shí)對(duì)象飛行器的空氣動(dòng)力學(xué)系數(shù)C j，如阻力系數(shù)CA、升力系數(shù)CN、側(cè)力系數(shù)CZ、滾轉(zhuǎn)力矩系數(shù)Cmx、俯仰力矩系數(shù)Cmy及偏航力矩系數(shù)Cmz，作為因變量；馬赫數(shù)Ma、攻角α、側(cè)滑角 β 及舵機(jī)偏轉(zhuǎn)角 δi（i=1,2,3,4）等作為自變量。xi（i=1,2,···,m，表示m個(gè)n維自變量），指定一組有序正整數(shù)序列 [k1,k2,···km]，得到自變量的候選單項(xiàng)式集合，也稱回歸因子：

對(duì)W0進(jìn)行QR 分解，得到初始化正交函數(shù)矩陣Q及系數(shù)矩陣R：

式中，θ=[θ1,θ2,···,θn]T，a=Rθ，e為模型誤差。正交建模的最終目的就是得到一組最優(yōu)θ（或a）值，從而使得模型的最小二乘代價(jià)函數(shù)J最小，其中：

式中，qj是矩陣Q的第j列，C j是氣動(dòng)系數(shù)C j在正交矩陣Q的正交列向量上的投影向量。這些向量的值表示了Q的正交列也即正交函數(shù)與C j的相關(guān)度。

計(jì)算公式（5）中令代價(jià)函數(shù)J對(duì) θ的導(dǎo)數(shù)為零，且有QTQ=I，求解方程便可得到未知參數(shù) θ的最優(yōu)估計(jì)：

公式（7）的矩陣形式為：

式中，rij（i=1,2,···,n；j=1,2,···,n）表示矩陣R的各個(gè)元素值。

1.2 參數(shù)辨識(shí)

氣動(dòng)系數(shù)快速建模的第二個(gè)部分是根據(jù)實(shí)時(shí)數(shù)據(jù)迭代更新正交函數(shù)池，同時(shí)更新的還有公式（8）右側(cè)的影響因子。在獲取新的數(shù)據(jù) [ξ1,ξ2,···,ξn]后，將新數(shù)據(jù)添加到公式（8）的最后一行，可得：

其中 ζ為新的因變量數(shù)據(jù)。為了保證滿足后續(xù)QR 分解條件，須將參數(shù)矩陣R′進(jìn)行變換，使其最后一行全為零。這里使用Givens 變換方法對(duì)R′進(jìn)行上三角化，定義Givens 矩陣如下：

公式（11）是完成一次數(shù)據(jù)更新的結(jié)果，其中帶有上標(biāo)“′”的符號(hào)表示包含新數(shù)據(jù)的信息，ε表示新數(shù)據(jù)的殘差，是新數(shù)據(jù)無(wú)法投影到正交基中的殘余信息。進(jìn)一步迭代時(shí)需要將矩陣式的最后一行用新數(shù)據(jù)代替，重復(fù)上述迭代過(guò)程，實(shí)現(xiàn)對(duì)正交函數(shù)池的更新。整個(gè)迭代過(guò)程只需進(jìn)行簡(jiǎn)單的矩陣乘法，不需要重復(fù)正交函數(shù)的生成過(guò)程，極大地提高了迭代效率。

模型均方擬合誤差（mean squared error,MSE）寫(xiě)為：

引入模型預(yù)測(cè)方差（prediction squared error,PSE）確定模型的結(jié)構(gòu)：

式中：N為加入模型中的正交函數(shù)個(gè)數(shù)。模型的預(yù)測(cè)誤差比模型的均方擬合誤差多了一項(xiàng)懲罰項(xiàng)常數(shù)是實(shí)時(shí)數(shù)據(jù)與已經(jīng)辨識(shí)得到的模型之間的誤差平方的上界估計(jì)值：

采用截止頻率為2 Hz 的二階巴特沃斯高通濾波器處理因變量數(shù)據(jù)，得到剔除噪聲的因變量數(shù)據(jù)vi；依據(jù)公式（14）得到實(shí)時(shí)的噪聲方差估計(jì)值以及的保守估計(jì)值。使用的保守估計(jì)值保證PSE 指標(biāo)比新數(shù)據(jù)的實(shí)際平方誤差更高。因此PSE 準(zhǔn)則保守估計(jì)了預(yù)測(cè)情況下的平方誤差。

僅依據(jù)MSE無(wú)法得到模型的最優(yōu)結(jié)構(gòu)，因?yàn)镸SE 函數(shù)是一個(gè)單調(diào)遞減函數(shù)，影響因子始終為負(fù)。因此，僅依靠MSE 一個(gè)指標(biāo)得到的模型會(huì)包含所有的正交函數(shù)，導(dǎo)致數(shù)據(jù)過(guò)度擬合。而懲罰項(xiàng)始終為正，因此PSE 準(zhǔn)則必然有一個(gè)全局最小值點(diǎn)。如圖1 是一個(gè)測(cè)試算例的尋優(yōu)結(jié)果，正交函數(shù)數(shù)量為10，在n=5處取得PSE 的全局最小值點(diǎn)。測(cè)試算例為：

圖1 測(cè)試算例尋優(yōu)結(jié)果Fig.1 Optimization results of a test example

利用正交函數(shù)進(jìn)行氣動(dòng)力學(xué)系數(shù)快速建模，PSE準(zhǔn)則能夠評(píng)估各個(gè)正交建模函數(shù)對(duì)模型的貢獻(xiàn)度，選取PSE 準(zhǔn)則的全局最優(yōu)點(diǎn)處的模型結(jié)構(gòu)既能保證模型的精確度，又能獲得最簡(jiǎn)單的模型結(jié)構(gòu)。采用PSE 準(zhǔn)則篩選之外，本文還引入了決定系數(shù)R2對(duì)qj進(jìn)一步篩選。進(jìn)行R2篩選是為了排除部分對(duì)模型的MSE 和PSE 影響均很小的正交函數(shù)。由公式（12）和公式（13）可知，當(dāng)正交函數(shù)的時(shí)，就會(huì)被選入模型中，無(wú)法排除正交函數(shù)對(duì)均方擬合誤差貢獻(xiàn)很小的函數(shù)。其中：

式中，Q′是從候選函數(shù)池中篩選出的m（m≤n）個(gè)正交函數(shù)，實(shí)際辨識(shí)過(guò)程中m?n，這也正是本文研究方法的優(yōu)勢(shì)所在。將公式（11）代入公式（16），可以得到多項(xiàng)式模型：

本文研究的多元正交化函數(shù)建模方法中，根據(jù)候選函數(shù)池中各函數(shù)對(duì)模型均方擬合誤差的貢獻(xiàn)度和預(yù)測(cè)誤差，確定模型結(jié)構(gòu)。本方法是在時(shí)域進(jìn)行建模，如果由于儀器遙測(cè)故障或計(jì)算延遲而丟失數(shù)據(jù)點(diǎn)，建模算法并不受影響（丟失數(shù)據(jù)點(diǎn)包含較多模型信息的情況除外）。因此，本算法對(duì)采樣數(shù)據(jù)的要求較低，對(duì)實(shí)時(shí)數(shù)據(jù)流中的間隙（時(shí)間）不敏感；但是數(shù)據(jù)大量丟包會(huì)對(duì)建模產(chǎn)生嚴(yán)重不利影響，這也是時(shí)域建模方法存在的普遍問(wèn)題。

2 數(shù)據(jù)采集與處理

本文實(shí)驗(yàn)數(shù)據(jù)是由F-16 風(fēng)洞試驗(yàn)仿真程序生成的。根據(jù)文獻(xiàn)[20]的介紹，仿真程序是對(duì)風(fēng)洞數(shù)據(jù)進(jìn)行線性插值從而得到氣動(dòng)力數(shù)據(jù)集，其風(fēng)洞數(shù)據(jù)從文獻(xiàn)[21]中整理得到（例如，圖2 是對(duì)軸向力系數(shù)風(fēng)洞數(shù)據(jù)插值的結(jié)果）；CxT是總的氣動(dòng)力系數(shù)，由Cx和Cxq組成。文獻(xiàn)[20] 中涉及的氣動(dòng)數(shù)據(jù)分別保存在Cx和Cxq兩張數(shù)據(jù)表中，故這里對(duì)兩組數(shù)據(jù)分別插值。

圖2 CxT 線性插值結(jié)果Fig.2 Linear interpolation of CxT

文獻(xiàn)[21]中的仿真程序使用線性插值，插值數(shù)據(jù)誤差較大，使用MOFM 得到的模型預(yù)測(cè)結(jié)果誤差很大（如圖3），決定系數(shù)R2=0.841 2。

圖3 CxT 線性插值模型預(yù)測(cè)結(jié)果Fig.3 Prediction of linear interpolated data

本文將仿真程序中的插值方法換為三次樣條函數(shù)插值（如圖4（a）），并使用MOFM 進(jìn)行建模，建模結(jié)果如圖4（b）。模型預(yù)測(cè)結(jié)果明顯提高，R2=0.999 4。同理，使用三次樣條函數(shù)插值方法插值得到了完整的氣動(dòng)力數(shù)據(jù)集。

圖4 CxT 樣條函數(shù)插值與建模結(jié)果Fig.4 Spline interpolation and prediction results of CxT

模型初始化前，線性相關(guān)性較強(qiáng)的數(shù)據(jù)會(huì)導(dǎo)致初始化矩陣為非滿秩矩陣，無(wú)法進(jìn)行QR 分解。在飛行器配平狀態(tài)下控制面沒(méi)有輸入擾動(dòng)時(shí)便開(kāi)始采樣，或采樣頻率高、飛機(jī)狀態(tài)短時(shí)間內(nèi)不會(huì)突變等，均會(huì)導(dǎo)致采集到的初始數(shù)據(jù)具有極高的線性相關(guān)性。因此，選取前nc條數(shù)據(jù)初始化前需要進(jìn)行線性相關(guān)性分析，剔除線性相關(guān)性較強(qiáng)的數(shù)據(jù)。本文使用皮爾森相關(guān)系數(shù)（R＜0.99）判斷數(shù)據(jù)之間的相關(guān)性，保證初始化矩陣為滿秩矩陣。

此外，需要對(duì)采集的數(shù)據(jù)進(jìn)行歸一化處理。氣動(dòng)數(shù)據(jù)的特殊性決定了不同變量擁有不同的量綱與量綱單位，在進(jìn)行氣動(dòng)數(shù)據(jù)分析建模時(shí)，為達(dá)到快速收斂目的，需要進(jìn)行歸一化處理，消除不同指標(biāo)變量之間的量綱影響。本文采用式（20）給出的最大最小標(biāo)準(zhǔn)歸一化方法對(duì)數(shù)據(jù)進(jìn)行歸一化處理：

3 建模結(jié)果與分析

對(duì)目標(biāo)飛行器進(jìn)行氣動(dòng)力學(xué)系數(shù)建模，需要分別建立6 個(gè)無(wú)量綱氣動(dòng)力和力矩系數(shù)的模型。在低馬赫數(shù)（Ma＜0.5）飛行條件下，氣動(dòng)力系數(shù)有攻角α、側(cè)滑角 β、升降舵機(jī)偏轉(zhuǎn)角 δe、副翼舵機(jī)偏轉(zhuǎn)角 δa、方向舵機(jī)偏轉(zhuǎn)角 δr、滾轉(zhuǎn)角速度p、俯仰角速度q和偏航角速度r的非線性函數(shù)[22]。目標(biāo)飛行器的部分參數(shù)以及飛行條件見(jiàn)表1。

表1 F-16 部分參數(shù)及飛行條件Table 1 Parameters of F-16 and flight conditions

進(jìn)行目標(biāo)飛行器實(shí)時(shí)氣動(dòng)力系數(shù)建模時(shí)，為了更加有效地激發(fā)目標(biāo)飛行器的運(yùn)動(dòng)模態(tài)[23]，需要使其在指定飛行高度與飛行速度下配平后保持平飛，并對(duì)3 個(gè)舵面同時(shí)輸入正交優(yōu)化復(fù)合正弦激勵(lì)[17,23]（圖5 為平飛配平狀態(tài)下加入的正交激勵(lì)信號(hào)），然后依據(jù)動(dòng)態(tài)響應(yīng)結(jié)果進(jìn)行建模。

圖5 正交優(yōu)化復(fù)合正弦激勵(lì)Fig.5 Orthogonal optimized multi-sine perturbation

將采集到的數(shù)據(jù)集經(jīng)隨機(jī)采樣分為訓(xùn)練集與驗(yàn)證集；將訓(xùn)練集中的數(shù)據(jù)逐條輸入算法中，進(jìn)而驗(yàn)證本文正交函數(shù)建模方法的實(shí)時(shí)性。通過(guò)正交函數(shù)建模得到的6 個(gè)無(wú)量綱氣動(dòng)力系數(shù)的多項(xiàng)式模型見(jiàn)表2，建模結(jié)果與仿真數(shù)據(jù)對(duì)比結(jié)果如圖6。

表2 氣動(dòng)力模型結(jié)構(gòu)Table 2 Aerodynamic models

圖6 氣動(dòng)力系數(shù)建模結(jié)果Fig.6 Aerodynamic coefficients predictions

實(shí)驗(yàn)結(jié)果表明，使用正交多項(xiàng)式建立的氣動(dòng)力系數(shù)模型有很好的預(yù)測(cè)結(jié)果。氣動(dòng)力系數(shù)模型誤差的主要原因是氣動(dòng)數(shù)據(jù)集考慮了前緣襟翼舵機(jī)偏轉(zhuǎn)角的影響，而在MOFM 方法的建模變量中沒(méi)有考慮其對(duì)橫向氣動(dòng)力系數(shù)的影響。此外，由于仿真程序中風(fēng)洞數(shù)據(jù)是在自變量空間內(nèi)均勻選點(diǎn)試驗(yàn)獲得，然后進(jìn)行樣條函數(shù)插值構(gòu)建數(shù)據(jù)集，所以真實(shí)試驗(yàn)數(shù)據(jù)與數(shù)據(jù)集有誤差。另外，在縱向力系數(shù)建模過(guò)程中沒(méi)有考慮橫向力的影響，也是模型誤差的部分原因。

建模的初始化過(guò)程、單條數(shù)據(jù)迭代過(guò)程以及模型參數(shù)辨識(shí)過(guò)程，耗時(shí)均在毫秒級(jí)（如表3），因此飛行試驗(yàn)中的采樣頻率可以達(dá)到500 Hz，模型辨識(shí)頻率也能達(dá)到200 Hz；而目標(biāo)飛行器的狀態(tài)在毫秒級(jí)內(nèi)不會(huì)有突變，完全滿足實(shí)時(shí)建模的要求。將基于多元正交函數(shù)建立的多項(xiàng)式模型與徑向基網(wǎng)絡(luò)（RBF 網(wǎng)絡(luò)）和具有兩層中間層的BP 網(wǎng)絡(luò)進(jìn)行比較可以發(fā)現(xiàn)：前者在預(yù)測(cè)精度上有明顯提升，6 個(gè)氣動(dòng)系數(shù)的預(yù)測(cè)誤差均小于RBF 網(wǎng)絡(luò)與BP 網(wǎng)絡(luò)預(yù)測(cè)值的誤差。

表3 建模結(jié)果Table 3 Results of modeling

4 結(jié)論

本文主要研究了利用多元正交函數(shù)對(duì)非線性氣動(dòng)力參數(shù)建模的方法，利用目標(biāo)飛行器在不同配平狀態(tài)下加入正交優(yōu)化復(fù)合正弦激勵(lì)得到的氣動(dòng)數(shù)據(jù)集驗(yàn)證了模型的有效性。本文的建模方法使用QR 分解矩陣后，利用Givens 矩陣進(jìn)行數(shù)據(jù)迭代，極大地減少了迭代過(guò)程消耗的時(shí)間，滿足實(shí)時(shí)建模的要求。而且由于R矩陣為上三角矩陣，在R矩陣上逐條迭代數(shù)據(jù)可以有效避免因數(shù)據(jù)相關(guān)性較強(qiáng)而導(dǎo)致的迭代過(guò)程中出現(xiàn)的病態(tài)系數(shù)矩陣。實(shí)驗(yàn)結(jié)果表明MOFM方法對(duì)高度非線性的氣動(dòng)力系數(shù)有較好的建模能力，說(shuō)明了該方法是處理多變量、非線性氣動(dòng)力問(wèn)題的一種真實(shí)可行的建模方法。針對(duì)一些較為復(fù)雜的氣動(dòng)模型，如多變量高階的多項(xiàng)式模型；采用本文方法建立氣動(dòng)力模型時(shí)需要的氣動(dòng)數(shù)據(jù)較多，這會(huì)對(duì)本文方法產(chǎn)生一定限制；后續(xù)研究中可與因子分析類方法結(jié)合，對(duì)變量與單項(xiàng)進(jìn)行篩選，可以減少正交函數(shù)池的規(guī)模，再進(jìn)行MOFM 建模。