基于多元正交函數的非線性快速氣動建模

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

(1.西安交通大學，西安 710049；2.中國空氣動力研究與發展中心，綿陽 621000；3.電子科技大學 航空航天學院，成都 611731)

0 引言

飛行器氣動力系數建模是研究飛行器控制系統設計、飛行控制、動態分析以及性能評估的基礎。隨著現代航空航天技術的不斷發展，現代戰機可以在更廣的包絡線內機動。飛行器在飛行過程中氣動環境復雜，氣動系數的模型呈現出多變量、高度非線性的特點。

在氣動力學建模領域，國內外已經逐漸形成了一套完整的體系，主要分為基于物理機理建模和基于數據建模兩個方向。基于物理機理建模的有氣動導數模型、狀態空間模型、階躍響應模型以及非線性微分方程模型[1-8]；這些模型的物理意義明確，適用于飛行器優化設計與性能評估。基于數據建模主要是運用各種智能學習算法進行建模，例如神經網絡模型、模糊邏輯模型[9-12]等；得益于這些智能算法的強大映射能力，基于數據建模可以非常有效地得到非線性領域內的模型。上述氣動系數建模方法通常包含大量重復的地面測試試驗，如風洞試驗與計算機仿真試驗；此外還需進行多次精確的飛行試驗，利用飛行試驗結果對地面建模結果進行修正，才能得到精確的全局氣動模型。

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

本文使用F-16 風洞試驗仿真程序[19]生成建模需要的數據集，通過逐條讀取訓練集中的氣動數據模擬真實飛行過程中的實時采樣數據；采用多元正交函數進行氣動力學系數快速建模，并對比了該方法與徑向基函數（RBF 網絡）、BP 網絡對風洞仿真數據的預測結果。研究結果表明，利用多元正交函數進行氣動力學系數建模能夠快速準確地建立氣動模型，為研究快速建模方法提供了新的思路與方法。

1 多元正交函數建模原理

基于實時飛行數據的多元正交函數非線性氣動學系數建模需要解決兩個問題：1）需要確定模型的結構。這需要從候選函數池的大量非線性項中剔除對模型擬合效果影響較小的項，采用最小預測均方誤差（predicted squared error,PSE）中準則可以滿足上述需求。2）需要根據實時數據辨識出結構參數。傳統的基于最小二乘回歸進行參數估計的方法會由于變量之間的線性相關性導致病態的回歸矩陣，造成參數估計不準。本文將回歸矩陣QR 分解后在R矩陣上迭代數據，解決了上述問題。

1.1 模型初始化

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

對W0進行QR 分解，得到初始化正交函數矩陣Q及系數矩陣R：

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

式中，qj是矩陣Q的第j列，C j是氣動系數C j在正交矩陣Q的正交列向量上的投影向量。這些向量的值表示了Q的正交列也即正交函數與C j的相關度。

計算公式（5）中令代價函數J對 θ的導數為零，且有QTQ=I，求解方程便可得到未知參數 θ的最優估計：

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

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

1.2 參數辨識

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

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

公式（11）是完成一次數據更新的結果，其中帶有上標“′”的符號表示包含新數據的信息，ε表示新數據的殘差，是新數據無法投影到正交基中的殘余信息。進一步迭代時需要將矩陣式的最后一行用新數據代替，重復上述迭代過程，實現對正交函數池的更新。整個迭代過程只需進行簡單的矩陣乘法，不需要重復正交函數的生成過程，極大地提高了迭代效率。

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

引入模型預測方差（prediction squared error,PSE）確定模型的結構：

式中：N為加入模型中的正交函數個數。模型的預測誤差比模型的均方擬合誤差多了一項懲罰項常數是實時數據與已經辨識得到的模型之間的誤差平方的上界估計值：

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

僅依據MSE無法得到模型的最優結構，因為MSE 函數是一個單調遞減函數，影響因子始終為負。因此，僅依靠MSE 一個指標得到的模型會包含所有的正交函數，導致數據過度擬合。而懲罰項始終為正，因此PSE 準則必然有一個全局最小值點。如圖1 是一個測試算例的尋優結果，正交函數數量為10，在n=5處取得PSE 的全局最小值點。測試算例為：

圖1 測試算例尋優結果Fig.1 Optimization results of a test example

利用正交函數進行氣動力學系數快速建模，PSE準則能夠評估各個正交建模函數對模型的貢獻度，選取PSE 準則的全局最優點處的模型結構既能保證模型的精確度，又能獲得最簡單的模型結構。采用PSE 準則篩選之外，本文還引入了決定系數R2對qj進一步篩選。進行R2篩選是為了排除部分對模型的MSE 和PSE 影響均很小的正交函數。由公式（12）和公式（13）可知，當正交函數的時，就會被選入模型中，無法排除正交函數對均方擬合誤差貢獻很小的函數。其中：

式中，Q′是從候選函數池中篩選出的m（m≤n）個正交函數，實際辨識過程中m?n，這也正是本文研究方法的優勢所在。將公式（11）代入公式（16），可以得到多項式模型：

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

2 數據采集與處理

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

圖2 CxT 線性插值結果Fig.2 Linear interpolation of CxT

文獻[21]中的仿真程序使用線性插值，插值數據誤差較大，使用MOFM 得到的模型預測結果誤差很大（如圖3），決定系數R2=0.841 2。

圖3 CxT 線性插值模型預測結果Fig.3 Prediction of linear interpolated data

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

圖4 CxT 樣條函數插值與建模結果Fig.4 Spline interpolation and prediction results of CxT

模型初始化前，線性相關性較強的數據會導致初始化矩陣為非滿秩矩陣，無法進行QR 分解。在飛行器配平狀態下控制面沒有輸入擾動時便開始采樣，或采樣頻率高、飛機狀態短時間內不會突變等，均會導致采集到的初始數據具有極高的線性相關性。因此，選取前nc條數據初始化前需要進行線性相關性分析，剔除線性相關性較強的數據。本文使用皮爾森相關系數（R＜0.99）判斷數據之間的相關性，保證初始化矩陣為滿秩矩陣。

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

3 建模結果與分析

對目標飛行器進行氣動力學系數建模，需要分別建立6 個無量綱氣動力和力矩系數的模型。在低馬赫數（Ma＜0.5）飛行條件下，氣動力系數有攻角α、側滑角 β、升降舵機偏轉角 δe、副翼舵機偏轉角 δa、方向舵機偏轉角 δr、滾轉角速度p、俯仰角速度q和偏航角速度r的非線性函數[22]。目標飛行器的部分參數以及飛行條件見表1。

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

進行目標飛行器實時氣動力系數建模時，為了更加有效地激發目標飛行器的運動模態[23]，需要使其在指定飛行高度與飛行速度下配平后保持平飛，并對3 個舵面同時輸入正交優化復合正弦激勵[17,23]（圖5 為平飛配平狀態下加入的正交激勵信號），然后依據動態響應結果進行建模。

圖5 正交優化復合正弦激勵Fig.5 Orthogonal optimized multi-sine perturbation

將采集到的數據集經隨機采樣分為訓練集與驗證集；將訓練集中的數據逐條輸入算法中，進而驗證本文正交函數建模方法的實時性。通過正交函數建模得到的6 個無量綱氣動力系數的多項式模型見表2，建模結果與仿真數據對比結果如圖6。

表2 氣動力模型結構Table 2 Aerodynamic models

圖6 氣動力系數建模結果Fig.6 Aerodynamic coefficients predictions

實驗結果表明，使用正交多項式建立的氣動力系數模型有很好的預測結果。氣動力系數模型誤差的主要原因是氣動數據集考慮了前緣襟翼舵機偏轉角的影響，而在MOFM 方法的建模變量中沒有考慮其對橫向氣動力系數的影響。此外，由于仿真程序中風洞數據是在自變量空間內均勻選點試驗獲得，然后進行樣條函數插值構建數據集，所以真實試驗數據與數據集有誤差。另外，在縱向力系數建模過程中沒有考慮橫向力的影響，也是模型誤差的部分原因。

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

表3 建模結果Table 3 Results of modeling

4 結論

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