丁娣，車競，錢煒祺，汪清

1.中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室，綿陽 621000 2. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

飛機飛行過程中的結冰現(xiàn)象是影響飛行安全甚至導致災難性事故的重要原因之一，國內(nèi)外發(fā)生的多起與結冰相關飛行事故表明開展飛機結冰及其防護相關基礎問題研究具有重要而緊迫的意義。飛機結冰探測是結冰防護的前提，從探測方式看主要有直接和間接探測兩種，直接探測以各種傳感器探測為主[1-2]，間接探測主要有以下3種方式[3]：一是基于狀態(tài)觀測的結冰探測技術，通過探測飛機狀態(tài)量的異常變化來對結冰進行診斷；二是新息法(Innovation Approach)探測技術，利用統(tǒng)計學工具分析結冰對新息序列(Innovation Sequence)的統(tǒng)計顯著性進行結冰診斷；三是基于參數(shù)辨識的結冰診斷方法。與直接探測方法相比，間接探測方法一般可以獲得更多的結冰相關信息，二者可以互為補充。NASA結冰研究中心就提出了將傳感器直接探測與參數(shù)辨識技術相結合，建立結冰智能管理系統(tǒng)(Smart Icing System，SIS)[4]的研究思路。

隨著飛機結冰智能管理、容錯飛行和可重構控制等需求的提出，參數(shù)辨識方法在飛機結冰探測領域受到越來越多的關注，通過參數(shù)辨識方法不僅能獲得飛機結冰累積過程中飛機氣動特性參數(shù)的變化，而且能夠獲得參數(shù)化的結冰累積模型。在飛機結冰探測中，要求參數(shù)辨識算法能夠實時快速捕捉由結冰引起的飛機氣動特性參數(shù)的微小變化，且精度高、響應時間短，并具有一定的魯棒性，因此濾波算法是一種較好的選擇。對于飛機結冰參數(shù)辨識問題，需要實時跟蹤飛機結冰氣動參數(shù)，因此屬于非線性濾波，目前應用于飛機結冰探測的參數(shù)辨識方法主要包括擴展卡爾曼濾波[5-6](Extended Kalman Filter，EKF)、H∞濾波和多模型自適應估計[7](Multiple Model Adaptive Estimator，MMAE)等。擴展卡爾曼濾波作為非線性濾波的經(jīng)典方法應用非常廣泛，但由于EKF自身存在的固有缺陷使得人們不斷尋求新的濾波方法[8]，比如專門為魯棒性設計的H∞濾波，相比EKF，H∞算法能夠處理建模誤差和噪聲不確定性[9]。H∞濾波算法首先由Didinsky等[10]提出并證明其漸進收斂性，由于計算效率高，響應快，且具有一定的魯棒性，近年來在飛機結冰參數(shù)在線辨識中應用較多。Melody等[6]利用NASA雙水獺飛機對比了批處理最小二乘算法、擴展卡爾曼濾波和H∞算法在結冰參數(shù)辨識中的性能，仿真發(fā)現(xiàn)在擾動和測量噪聲的影響下，只有H∞算法能夠提供實時準確的辨識結果。在此基礎上，Melody等[11]進一步探討了H∞算法在飛機結冰在線識別中的應用，討論了算法對中度和重度結冰累積過程中飛機氣動參數(shù)連續(xù)變化的跟蹤能力。由于H∞算法的魯棒性較好且在線跟蹤能力較強，Melody認為H∞算法適合應用于結冰管理系統(tǒng)(Ice Management System，IMS)[12]。Schuchard等[13]結合Melody的研究內(nèi)容，在H∞算法參數(shù)辨識的基礎上結合神經(jīng)網(wǎng)絡方法對飛機結冰位置和嚴重性進行分類和識別，Dong和Ai[14-15]后續(xù)也開展了相似的研究。應思斌等[16]也對H∞算法在飛機結冰參數(shù)辨識中的應用進行了探討，分析了將時變和時不變H∞算法靈活應用于飛機結冰探測的策略[17]，并基于H∞參數(shù)辨識和H2控制算法[18]，設計了飛機結冰可重構控制系統(tǒng)。

目前H∞濾波在飛機結冰參數(shù)辨識中的研究大部分都是基于NASA的雙水獺結冰研究樣機，研究中發(fā)現(xiàn)當存在測量噪聲時算法存在跟蹤延時和辨識精度降低的現(xiàn)象，這是由H∞算法本身關于噪聲不確定性假設導致的，針對該問題，本文擬通過調(diào)節(jié)算法中關于噪聲不確定性相關參數(shù)加以解決。在此基礎上，本文針對為滿足國內(nèi)大型飛機結冰防護研究需求而專門設計的大型結冰研究樣機，開展基于H∞算法的飛機機翼結冰氣動參數(shù)辨識方法研究。文中首先介紹了結冰研究樣機無冰、中度和重度結冰下的縱向氣動導數(shù)的CFD擬合結果，建立了考慮結冰累積過程的動力學模型及時變參數(shù)H∞辨識算法，利用動力學模型和H∞算法構建能夠用于結冰在線探測的仿真框架，通過辨識算例分析了H∞算法參數(shù)選擇的合理性以及H∞算法用于結冰研究樣機氣動導數(shù)辨識的有效性，最后分析了H∞算法對81種不同結冰累積過程的辨識能力以及不同測量噪聲對算法辨識精度的影響。

1 結冰研究樣機

為了滿足大型飛機的研制需求，分析結冰對大型飛機氣動特性的影響，課題組專門構建了與C919、Boeing737和空客A320外形相近的大型飛機作為該結冰項目的研究樣機，該樣機同時具備大型飛機的基本空氣動力學和飛行力學特征，其外形結構如圖1所示，幾何物理參數(shù)如表1所示。

結冰研究樣機的氣動特性數(shù)據(jù)由CFD手段得到，基于雷諾平均Navier-Stokes(RANS)方法和多塊結構化網(wǎng)格技術對結冰研究樣機進行基本氣動力計算，得到了干凈外形及機翼中度和重度結冰下的縱向氣動特性，機翼中度和重度結冰對應的冰型如圖2所示。

表1 結冰研究樣機幾何參數(shù)Table 1 Geometric parameters of icing research prototype

表2 結冰研究樣機無冰及結冰情況下的氣動導數(shù)

2 飛機結冰氣動參數(shù)辨識方法及算例

首先推導了考慮結冰累積過程的飛行動力學模型以及能夠用于時變參數(shù)辨識的H∞濾波算法，利用動力學模型和參數(shù)辨識算法構建可用于飛機結冰在線探測的仿真框架，在此基礎上討論了辨識算法參數(shù)對辨識精度的影響，并由仿真算例分析了氣動導數(shù)的辨識誤差及結冰預報延時。

2.1 飛機結冰飛行動力學模型

建立考慮結冰累積過程的飛機縱向運動方程，在速度坐標系下建立質(zhì)心運動方程，在體坐標系下建立轉動方程，則飛機的縱向動力學方程可表示為

Vsinαcosβcosθcosφ

(1)

式中:V為飛機合速度；α為迎角；β為側滑角；p為滾轉角速率；q為俯仰角速率；r為偏航角速率；θ為俯仰角；φ為滾轉角；h為離地高度；q∞為動壓；P=[PxPyPz]為發(fā)動機推力；M=[MxMyMz]為發(fā)動機推力力矩；g=[gxgygz]為重力加速度矢量，下標x、y和z表示相應坐標系下的坐標；CD、CL和Cm分別為阻力、升力和俯仰力矩系數(shù)，其中阻力、升力和俯仰力矩系數(shù)相關氣動導數(shù)取值見表2,表達式為

(2)

飛機結冰會引起阻力、升力和俯仰力矩系數(shù)的變化，描述結冰對飛機氣動特性的影響需要建立氣動導數(shù)在結冰累積過程中的變化模型。伊利諾斯大學的Bragg等[4]提出了一種描述結冰對氣動導數(shù)影響的數(shù)學模型，該模型雖然不能完全精確反應結冰對氣動導數(shù)的影響，但是通過該模型能夠獲得結冰前后氣動導數(shù)的一般變化趨勢，且能夠用于智能結冰系統(tǒng)的開發(fā)，目前已被應用于飛機結冰在線識別技術的討論[7,11,15]。結冰對飛機氣動導數(shù)的影響可表示為

(3)

式中：C*表示飛機的任意氣動導數(shù)，上標“iced”表示結冰后的氣動導數(shù)，“clean”表示干凈外形的氣動導數(shù)；ηice為結冰嚴重性參數(shù)；KC*為氣動導數(shù)變化斜率，由結冰前后氣動導數(shù)大小決定。ηice=0表示干凈外形的氣動特征，ηice=1表示機翼結冰后飛機的氣動特征，通過構造ηice的變化曲線來近似飛機氣動導數(shù)受結冰影響后的變化過程。Melody等[11]給出了一種反應結冰連續(xù)增長的結冰模型，考慮了結冰過程同時受大氣環(huán)境因素和冰累計量的影響，將ηice的變化記為

(4)

式中:N1和N2為常數(shù)，與不同時刻的結冰嚴重性參數(shù)大小相關，N2代表式(4)的非線性程度，N2>0為冰增長過程，N2<0為冰退化過程，二者可表示為

(5)

式中:Tcld為結冰累積時間。

式(4)為結冰增長模型，dη表征大氣環(huán)境對結冰的影響，這里將其表示為上升余弦函數(shù)，即

(6)

當Tcld、ηice(Tcld/2)和ηice(Tcld)已知時，由式(3)～式(6)能夠給出飛機任意氣動導數(shù)受結冰影響隨時間的變化曲線，其中ηice(Tcld/2)和ηice(Tcld)的取值表征不同結冰累積時間對應的結冰嚴重性系數(shù)，可以反應結冰累積速度。

2.2 H∞時變參數(shù)辨識方法

將式(1)對應的飛行動力學系統(tǒng)寫成參數(shù)化形式[10]，即

(7)

(8)

(9)

將式(4)代入式(9)，得

(10)

在初始時刻ηice=0，那么聯(lián)立式(6)、式(8)、式(10)，得

(11)

(12)

(13)

根據(jù)二次型變換，由式(13)可推導出H∞濾波算法，即

(14)

式中：Σ∈R(n1+n2)×(n1+n2)，n1為系統(tǒng)狀態(tài)量維度，n2為系統(tǒng)未知參數(shù)維度，則Σ1∈n1×n1，當選取P0=I，以及時，對所有Q0>0都有γ*≡1成立。

2.3 算法參數(shù)調(diào)節(jié)及辨識結果分析

利用辨識算例來說明H∞算法在飛機結冰參數(shù)辨識中的有效性，由于目前缺乏飛行試驗數(shù)據(jù)，只能采用仿真數(shù)據(jù)作為系統(tǒng)辨識依據(jù)。

假定飛機在定直平飛情況下機翼出現(xiàn)結冰，結冰累積過程由式(3)描述，其中Tcld=200 s，ηice(Tcld)=1.0，ηice(Tcld/2)=0.7，飛機巡航高度為5 000 m，飛行速度為馬赫數(shù)0.3，升降舵偏輸入信號幅值為3°、周期為40 s的正弦信號?？紤]測量噪聲的影響，假設噪聲為高斯白噪聲，其標準差取值參照A340的傳感器測量精度，如表3[21]所示。在該初始條件下生成考慮測量噪聲的動力學仿真數(shù)據(jù)，利用H∞算法辨識受機翼結冰影響較大的3個與迎角相關的氣動導數(shù)，即CDα、CLα和Cm α，對辨識結果進行分析以評估H∞算法的辨識精度。

表3 測量噪聲標準差[21]Table 3 Standard deviations of measurement noises[21]

由H∞算法推導，算法初始需要確定的參數(shù)包括Q0和γ，且參數(shù)的取值會對辨識結果產(chǎn)生影響。對于擾動衰減水平γ，根據(jù)Didinsky的研究結果，當γ的取值接近且不等于γ*時，算法收斂效果最好[10]；Q0為系統(tǒng)未知參數(shù)初始偏差的權值矩陣，由于這一偏差的存在使得Q0的取值會對辨識結果產(chǎn)生影響。Melody等在對NASA的雙水獺結冰研究樣機進行結冰氣動參數(shù)辨識研究時，將這兩個參數(shù)分別取為γ=1.01，Q0=(1×10-5)I，針對兩種不同結冰累積情況對Cm α變化情況進行跟蹤，發(fā)現(xiàn)辨識結果相對真值存在較大的跟蹤延時，見表4[11]。

表4 雙水獺飛機結冰預報延時[11]Table 4 Delays in ice detection of Twin Otter Airplane[11]

根據(jù)以上分析，對于本文的問題，可將γ的取值設為與雙水獺飛機算例一致，通過調(diào)節(jié)矩陣Q0對角線上的元素取值來達到提高參數(shù)辨識精度以及減小預報延時的目的。

利用前面推導的飛機動力學模型和參數(shù)辨識算法，構建含動力學仿真和H∞參數(shù)辨識為一體的結冰在線探測框架，將前面所給仿真初始條件和算法參數(shù)代入仿真框架，得到3個受結冰影響較為嚴重的氣動導數(shù)辨識結果，根據(jù)辨識結果與真值之間的平方根誤差(Root Mean Square，RMS)分析矩陣Q0的3個對角線元素{q11,q22,q33}對3個氣動導數(shù)辨識精度的影響，如表5所示。表中對每個矩陣對角線元素選取3個不同點，通過H∞算法對參數(shù)的27種組合分別進行辨識，得到3個氣動導數(shù)相對真值歸一化后的平方根誤差，對每個矩陣元素相同的結果取平均，得到每一矩陣元素取值對應的氣動導數(shù)辨識精度。由表中數(shù)據(jù)可知，q11主要對CDα的辨識精度產(chǎn)生影響，對CL α和Cm α的影響幾乎可以忽略，且q11越小辨識誤差越大；q22主要對CL α和Cm α的辨識精度產(chǎn)生影響，對CD α影響相對較小，q22越小CL α和Cm α的辨識精度越高；q33對Cm α影響相對較大，對CD α影響相對較小，對CL α影響幾乎可以忽略，q33越小Cm α的辨識精度越高。根據(jù)以上分析，這里將矩陣Q0設為Q0=diag(1,10-5,10-5)。

表5 矩陣Q0對氣動導數(shù)辨識精度的影響

將前面確定的H∞算法參數(shù)代入結冰在線探測框架，得到氣動導數(shù)辨識結果，如圖3和圖4所示。圖3中給出了CD α、CL α和Cm α3個氣動導數(shù)在結冰累積影響下的辨識結果與真值之間的對比，圖4中給出了由式(3)將辨識結果表示為結冰嚴重性參數(shù)后與真值的對比，圖中黑色虛線為移動平均濾波器(Moving Average Filter)平滑后的辨識結果。由圖中結果可知，CL α和Cm α的辨識結果較為平滑，CD α在氣動導數(shù)變化階段振蕩相對比較明顯，這說明通過前面調(diào)節(jié)加權矩陣Q0能夠有效減小CL α和Cm α的振蕩，但CD α的振蕩現(xiàn)象難以消除；當結冰累積過程結束后，辨識結果能夠收斂到真值附近，即趨近于重度結冰后的氣動導數(shù)取值；在結冰累積階段，辨識結果存在一定程度上的時間延遲，其中CD α的延遲較為嚴重，Cm α次之，CL α最好。

為了對辨識結果進行評估，定義了平方根誤差、最大絕對誤差和結冰嚴重性參數(shù)分別為0.25、0.50、0.75和0.95對應的時間延遲作為評估的指標，為了保證數(shù)據(jù)的可比性，這里針對歸一化后的結冰嚴重參數(shù)來評估，其中誤差和時間延遲的定義均相對氣動導數(shù)真值給出，時間延遲的計算采用了數(shù)據(jù)平滑后的結果。表6中給出了辨識結果的誤差及結冰嚴重性參數(shù)預報延遲，由表中數(shù)據(jù)可知，無論是辨識誤差還是時間延遲，都是CL α辨識結果最好，Cm α次之，CD α最差。由該算例的辨識結果分析，發(fā)現(xiàn)3個氣動導數(shù)的辨識結果平方根誤差較小，其中辨識效果最差的CD α的平方根誤差也在真值的11%左右，辨識效果最好的CLα僅為真值的1.4%；最大絕對誤差CL α僅為真值的7.8%，Cm α為16.4%，而效果最差的CD α為41%左右，這可能是由辨識結果出現(xiàn)振蕩導致的；如果分別以3個氣動導數(shù)的辨識結果作為結冰探測的依據(jù)，那么CD α最大會有20 s左右的預報延遲，Cm α的預報延遲在10 s以內(nèi)，CL α的預報延遲在3 s以內(nèi)。同時，表6中Cm α的預報延時明顯小于表4中的雙水獺飛機算例，說明通過調(diào)節(jié)算法加權矩陣參數(shù)能夠提高算法對飛機結冰的在線探測能力。

表6 辨識結果誤差及結冰預報延時Table 6 Identification results errors and ice detection delays

因此，對于結冰研究樣機機翼結冰情況，通過3個氣動導數(shù)的辨識結果均可對結冰嚴重性程度進行預報，但三者的辨識精度和預報延遲差別較大，可以從中選擇預報效果最好的參數(shù)辨識結果作為結冰在線探測的依據(jù)。

3 飛機機翼結冰探測精度

為了評估H∞算法在不同工況下的結冰探測能力，下面通過仿真分析不同結冰累積過程以及測量噪聲大小對辨識結果的影響。

3.1 結冰累積過程影響分析

考慮結冰累積時間Tcld和ηice(Tcld/2)取值的變化，可以得到不同的結冰累積過程，利用H∞算法對不同結冰累積過程進行辨識，通過辨識結果誤差可評估算法對結冰累積過程的辨識能力。這里將結冰累積時間Tcld的取值區(qū)間定為[50,450]s，ηice(Tcld/2)的取值區(qū)間定為[0.1,0.9]，不考慮測量噪聲的影響，分析辨識結果與真值之間的誤差以及不同結冰水平對應的時間延遲，如圖5～圖7和表7所示。

圖5和圖6給出了不同結冰累積過程下的平方根誤差和最大絕對誤差結果，由于平方根誤差在結冰累積時間較長時會顯著增大，為了便于分析，因此將結果分成了兩部分，見圖5(a)和圖5(b)。圖5和圖6中分別給出了81種不同結冰累積過程對應的CD α、CL α和Cm α辨識結果與真值之間的誤差，為了便于對比3個氣動導數(shù)辨識誤差,分別用不同顏色的標記區(qū)分。對于平方根誤差，當Tcld≥350 s且ηice(Tcld/2)≤0.7時，3個氣動導數(shù)均顯著增大，最大可達真值的50%左右，其余結冰累積情況下平方根誤差均較小，最大約為真值的10%左右，這說明當結冰時間過長且累積速度緩慢時，比較不容易得到準確的辨識結果；對于這81種不同的結冰情況，其中CL α的平方根誤差最小的情況占95%，CD α最小的情況占5%，而CD α在60%的情況下誤差最大，Cm α在40%的情況下誤差最大；且當ηice(Tcld/2)取值在0.4～0.7之間時，CD α辨識誤差相對較大，其余取值情況Cm α辨識誤差相對較大，說明當結冰累積速度較慢或較快時Cm α辨識結果更不準確。對于最大絕對誤差，其結果與平方根誤差的分析結果近似，同樣對于結冰時間較長且結冰速度較慢的情況,辨識最大誤差絕對值會較大；而約有89%的情況下CL α的辨識最大絕對誤差最小，6%的情況下CDα最小，5%的情況下Cm α最?。幌鄳腃D α有62%的情況辨識最大絕對誤差最大，Cm α有38%的情況辨識最大絕對誤差最大。

對辨識結果關于不同結冰水平的時間延遲進行分析，發(fā)現(xiàn)當結冰累積時間較長時存在算法無法跟蹤得到相應結冰水平的情況。圖7中給出了不同結冰累積情況下的3個氣動導數(shù)辨識結果分別在0.25、0.50、0.75和0.95這4個結冰水平上的跟蹤失敗情況，由圖可知跟蹤失敗的情況主要集中在結冰累積時間300～450 s之間。若算法無法跟蹤結冰水平為0.25的情況，那么意味著大于0.25的情況也無法跟蹤，因此圖中4個不同結冰水平的跟蹤失敗標記呈從高水平向低水平覆蓋的狀態(tài)，且逐漸向(450, 0.1)邊界點收縮，其中CDα、CLα和Cmα在不同結冰水平下的預報延時和跟蹤失敗次數(shù)統(tǒng)計見表7，由表中數(shù)據(jù)可知CLα的跟蹤失敗次數(shù)最少，延時均值最小，說明CLα跟蹤更為準確，表中延時標準差在結冰水平為0.25和0.95時較大，這可能是由于辨識結果振蕩造成延時散布較大。

根據(jù)以上分析，CL α的辨識精度相對其余兩個氣動導數(shù)較高，下面根據(jù)CL α的平方根誤差和最大絕對誤差大小，分析81種結冰累積過程的相對辨識精度，如圖8所示。圖中以不同顏色對81種結冰累積過程的辨識精度進行區(qū)分，偏藍為辨識精度較高，偏紅為辨識精度較差。由圖可知，結冰累積時間對辨識精度影響較明顯，累積時間較長(>300 s)和較短(50 s)辨識精度稍差，累積時間適中(100～300 s)辨識精度較高，對于ηice(Tcld/2)，累積時間較短和較長對辨識精度有一定影響，當累積時間較長時(>300 s)，結冰累積速度越快辨識精度越高，當累積時間較短時(50 s)，ηice(Tcld/2)取值適中(約為0.3),辨識精度較高。這一結果與H∞算法對飛機結冰氣動參數(shù)變化的捕捉能力相關，參數(shù)變化率太大或太小均會影響算法辨識精度，只有參數(shù)變化率適中才能得到精度較高的辨識結果。

表7 不同結冰水平下預報延時及失敗次數(shù)統(tǒng)計Table 7 Delay and failure statistics in different icing levels

3.2 隨機噪聲影響下的氣動導數(shù)辨識精度

由于實際飛行過程中測量噪聲的存在，必然會影響結冰氣動導數(shù)辨識精度，為了對實際飛行過程中基于H∞的結冰氣動導數(shù)辨識精度進行定量評估，下面在給定測量噪聲標準差大小范圍下進行多次蒙特卡羅仿真，分析結冰氣動導數(shù)的辨識精度。

為了盡量覆蓋實際可能出現(xiàn)的噪聲情況，以表3中A340的傳感器精度為基礎，將測量噪聲標準差變化范圍取為從0到表中對應的10倍大小，在該范圍內(nèi)隨機均勻生成大量測量噪聲標準差作為辨識精度分析的輸入條件，根據(jù)每一標準差生成高斯隨機噪聲代入仿真辨識框架，對每一飛機結冰氣動導數(shù)辨識結果精度進行統(tǒng)計，得到給定噪聲標準差變化范圍下氣動導數(shù)辨識精度。蒙特卡羅仿真次數(shù)的選擇需要分析氣動導數(shù)辨識結果誤差的收斂性，表8中給出了CDα、CLα和Cmα不同蒙特卡羅仿真次數(shù)對應平方根誤差的統(tǒng)計均值，由表中數(shù)據(jù)可知對于4 000、5 000和6 000次仿真，3個氣動導數(shù)的平方根誤差的統(tǒng)計均值幾乎沒有變化，已接近于收斂，因此根據(jù) 6 000次蒙特卡羅仿真結果分析3個氣動導數(shù)的H∞算法辨識精度具有一定的意義。

表8 不同蒙特卡羅仿真次數(shù)下的氣動導數(shù)平方根誤差統(tǒng)計均值

表9和表10中給出了CDα、CLα和Cmα3個氣動導數(shù)辨識結果的平方根誤差、最大絕對誤差和4種不同結冰水平預報延時按正態(tài)分布擬合得到的統(tǒng)計均值μ、標準差σ以及二者的95%置信區(qū)間范圍μci和σci。由表中數(shù)據(jù)可知，對于平方根誤差和最大絕對誤差，從統(tǒng)計均值和標準差可知，都有CDα誤差最大，Cmα次之，CLα誤差最小，其中辨識效果最好的CLα平方根誤差均值為1.8%左右，最大絕對誤差均值為8%左右，Cmα的平方根誤差均值為4%左右，最大絕對誤差均值為19%，而CDα的平方根誤差均值為52%，最大絕對誤差均值高達250%，這說明在測量噪聲影響下，CDα的辨識結果可信性較低，而CLα和Cmα能夠得到較為準確的辨識結果。由表10中數(shù)據(jù)可知，對于不同結冰水平下的預報延時統(tǒng)計結果，由于CDα辨識精度較差，6 000次仿真中出現(xiàn)了6次利用CDα辨識結果無法跟蹤結冰系數(shù)為0.95的情況，且統(tǒng)計均值出現(xiàn)了負值，說明CDα跟蹤結果在結冰累積后期出現(xiàn)了較大的振蕩，根據(jù)辨識結果跟蹤該結冰系數(shù)出現(xiàn)了早于真值的情況；對4種不同結冰水平，3個氣動導數(shù)的預報延時統(tǒng)計均值除CLα明顯較小之外，CDα和Cmα相當，但統(tǒng)計標準差CDα明顯大于其余二者，說明由于CDα的辨識精度較低導致結冰預報時間延時散布較大；根據(jù)統(tǒng)計結果，4種不同結冰水平下CLα的預報延時均值最大為3 s左右，

表9 氣動導數(shù)辨識結果平方根誤差和最大絕對誤差統(tǒng)計Table 9 Statistics of RMS and maximum absolute errors of identification results of aerodynamic derivatives

表10 不同結冰水平下氣動導數(shù)辨識延時統(tǒng)計Table 10 Delays in aerodynamic derivatives identification in different icing levels

統(tǒng)計標準差隨結冰水平增大而增加，最大不超過2.5 s，Cmα的跟蹤延時均值最大不超過9.5 s，統(tǒng)計標準差最大為2.8 s左右。

4 結 論

1) 調(diào)節(jié)H∞算法中的加權矩陣會對辨識精度產(chǎn)生影響，通過選擇合適的算法參數(shù)使得H∞算法能夠對結冰研究樣機氣動導數(shù)結冰累積過程進行有效辨識，其中辨識效果最差的CDα的歸一化平方根誤差為真值的11%左右，辨識效果最好的CLα僅為真值的1.4%。

2) 當參數(shù)變化率太快或太慢時H∞算法辨識精度會降低，當參數(shù)變化率適中，即結冰累積時間適中(100～300 s)時，算法辨識精度較高，這與H∞算法對時變參數(shù)變化的捕捉能力相關。

3) 在測量噪聲影響下，CDα的辨識結果可信度較低，而CLα和Cmα能夠得到較為準確的辨識結果，在文中給定的測量噪聲標準差變化范圍內(nèi)，CLα歸一化后的平方根誤差均值為1.8%左右，Cmα的平方根誤差均值為4%左右，4種不同結冰水平下CLα的預報延時均值最大為3 s左右，Cmα的預報延時均值最大不超過9.5 s。

4) 由大量辨識結果分析可知，對于受機翼結冰影響較為嚴重的3個氣動導數(shù)，在結冰累積時間適中的情況下，都有CLα和Cmα辨識結果較好，CDα較差，這與飛機本身的氣動特性相關，由于CDα相對CLα和Cmα為小量，因此H∞時變參數(shù)模型建模偏差和測量噪聲等擾動導致CDα辨識結果容易出現(xiàn)振蕩，辨識精度降低，在實際結冰探測中，可以選擇法向氣動導數(shù)和俯仰力矩導數(shù)作為辨識對象。

本文的研究工作可以作為大型飛機結冰防護系統(tǒng)研究的基礎，為開展虛擬飛行結冰探測提供算法支持，得到的算法辨識精度定量分析結論可為后續(xù)工作的開展提供重要的評估依據(jù)。下一步還可以結合可重構控制開展基于結冰研究樣機的結冰防護相關研究。

參 考 文 獻

[1] GORAJ Z.An overview of the deicing and anti-icing technologies with prospects for the future[C]∥24th International Congress of Aeronautical Sciences, 2004.

[2] JARVINEN P. Aircraft ice detection method: AIAA-2007-696[R].Reston, VA: AIAA, 2007.

[3] CALISKAN F, HAJIYEV C. A review of inflight detection and identification of aircraft icing and reconfigurable control[J]. Progress in Aerospace Sciences, 2013, 60: 12-34.

[4] BRAGG M B, BASAR T, PERKINS W R, et al. Smart icing systems for aircraft icing safety: AIAA-2002-0813[R]. Reston, VA: AIAA, 2002.

[5] WENZ A, JOHANSEN T A. Icing detection for small fixed wing UAVs using inflight aerodynamic coefficient estimation[C]∥IEEE Conference on Control Applications. Piscataway, NJ: IEEE Press, 2016.

[6] MELODY J W, BASAR T, PERKINS W R, et al. Parameter identification for inflight detection and characterization of aircraft icing[J]. Control Engineering Practice, 2000, 8(9): 985-1001.

[7] CRISTOFARO A, JOHANSEN T A, AGUIAR A P. Icing detection and identification for unmanned aerial vehicles: Multiple model adaptive estimation[C]∥European Control Conference, 2015.

[8] 占榮輝，張軍. 非線性濾波理論與目標跟蹤應用[M]. 北京：國防工業(yè)出版社，2013: 8-11.

ZHAN R H, ZHANG J. Nonlinear filtering theory with target tracking application[M]. Beijing: National Defense Industry Press, 2013: 8-11(in Chinese).

[9] SIMON D. 最優(yōu)狀態(tài)估計:卡爾曼,H∞及非線性濾波[M].張勇剛, 李寧, 奔粵陽, 譯. 北京：國防工業(yè)出版社，2015: 255-269.

SIMON D. Optimal state estimation: Kalman,H∞and nonlinear approaches[M]. ZHANG Y G, LI N, BEN Y Y, translated. Beijing: National Defense Industry Press, 2015: 255-269(in Chinese).

[10] DIDINSKY G, PAN Z, BASAR T. Parameter identification for uncertain plants usingH∞methods[J]. Automatica, 1995, 31(9): 1227-1250.

[11] MELODY J W, HILLBRAND T, BASAR T, et al.H∞parameter identification for inflight detection of aircraft icing: The time-varying case[J]. Control Engineering Practice, 2001, 9(12): 1327-1335.

[12] MELODY J W. Inflight characterization of aircraft icing[D]. Illinois Urbana: University of Illinois at Urbana-Champaign Graduate College, 2004: 1-6.

[13] SCHUCHARD E A, MELODY J W, BASAR T, et al. Detection and classification of aircraft icing using Neural Networks: AIAA-2000-0361[R]. Reston, VA: AIAA, 2000.

[14] DONG Y Q, AI J L. Research on inflight parameter identification and icing location detection of the aircraft[J]. Aerospace Science and Technology, 2013, 29(1): 305-312.

[15] DONG Y Q, AI J L. Inflight parameter identification and icing location detection of the aircraft: The time-varying case[J/OL]. Journal of Control Science and Engineering, 2014: 1-11. [2014-07-10].http://dx.doi.org/10.1155/2014/396532.

[16] 應思斌, 葛彤, 艾劍良. 飛機結冰時不變參數(shù)辨識技術[J]. 指揮控制與仿真, 2012, 34(4): 55-60.

YING S B, GE T, AI J L. Time invariant parameter identification of inflight aircraft icing[J]. Command Control & Simulation, 2012, 34(4): 55-60(in Chinese).

[17] 應思斌, 葛彤, 艾劍良. 飛機結冰氣動參數(shù)綜合檢測方法研究[J]. 指揮控制與仿真, 2012, 34(5): 128-133.

YING S B, GE T, AI J L. Research on comprehensive parameter identification of inflight aircraft icing[J]. Command Control & Simulation, 2012, 34(5): 128-133(in Chinese).

[18] YING S B, GE T, AI J L.H∞parameter identification and H2feedback control synthesizing for inflight aircraft icing[J]. Journal of Shanghai Jiaotong University, 2013, 18(3): 317-325.

[19] RATVASKY T P, RANAUDO R J. Icing effects on aircraft stability and control determined from flight data: AIAA-1993-0398[R]. Reston, VA: AIAA, 1993.

[20] BRAGGM B, HUTCHISON T, MERRET J, et al. Effect of ice accretion on aircraft flight dynamics: AIAA-2000-0360[R]. Reston, VA: AIAA, 2000.

[21] AYKAN R, HAJIYEV C, CALISKAN F. Aircraft icing detection, identification and reconfigurable control based on Kalman filtering and neural networks: AIAA-2005-6220[R]. Reston, VA: AIAA, 2005.