丁娣，車(chē)競(jìng)，錢(qián)煒祺，汪清

1.中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，綿陽(yáng) 621000 2. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000

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

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

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

1 結(jié)冰研究樣機(jī)

為了滿(mǎn)足大型飛機(jī)的研制需求，分析結(jié)冰對(duì)大型飛機(jī)氣動(dòng)特性的影響，課題組專(zhuān)門(mén)構(gòu)建了與C919、Boeing737和空客A320外形相近的大型飛機(jī)作為該結(jié)冰項(xiàng)目的研究樣機(jī)，該樣機(jī)同時(shí)具備大型飛機(jī)的基本空氣動(dòng)力學(xué)和飛行力學(xué)特征，其外形結(jié)構(gòu)如圖1所示，幾何物理參數(shù)如表1所示。

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

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

表2 結(jié)冰研究樣機(jī)無(wú)冰及結(jié)冰情況下的氣動(dòng)導(dǎo)數(shù)

2 飛機(jī)結(jié)冰氣動(dòng)參數(shù)辨識(shí)方法及算例

首先推導(dǎo)了考慮結(jié)冰累積過(guò)程的飛行動(dòng)力學(xué)模型以及能夠用于時(shí)變參數(shù)辨識(shí)的H∞濾波算法，利用動(dòng)力學(xué)模型和參數(shù)辨識(shí)算法構(gòu)建可用于飛機(jī)結(jié)冰在線(xiàn)探測(cè)的仿真框架，在此基礎(chǔ)上討論了辨識(shí)算法參數(shù)對(duì)辨識(shí)精度的影響，并由仿真算例分析了氣動(dòng)導(dǎo)數(shù)的辨識(shí)誤差及結(jié)冰預(yù)報(bào)延時(shí)。

2.1 飛機(jī)結(jié)冰飛行動(dòng)力學(xué)模型

建立考慮結(jié)冰累積過(guò)程的飛機(jī)縱向運(yùn)動(dòng)方程，在速度坐標(biāo)系下建立質(zhì)心運(yùn)動(dòng)方程，在體坐標(biāo)系下建立轉(zhuǎn)動(dòng)方程，則飛機(jī)的縱向動(dòng)力學(xué)方程可表示為

Vsinαcosβcosθcosφ

(1)

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

(2)

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

(3)

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

(4)

式中:N1和N2為常數(shù)，與不同時(shí)刻的結(jié)冰嚴(yán)重性參數(shù)大小相關(guān)，N2代表式(4)的非線(xiàn)性程度，N2>0為冰增長(zhǎng)過(guò)程，N2<0為冰退化過(guò)程，二者可表示為

(5)

式中:Tcld為結(jié)冰累積時(shí)間。

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

(6)

當(dāng)Tcld、ηice(Tcld/2)和ηice(Tcld)已知時(shí)，由式(3)～式(6)能夠給出飛機(jī)任意氣動(dòng)導(dǎo)數(shù)受結(jié)冰影響隨時(shí)間的變化曲線(xiàn)，其中ηice(Tcld/2)和ηice(Tcld)的取值表征不同結(jié)冰累積時(shí)間對(duì)應(yīng)的結(jié)冰嚴(yán)重性系數(shù)，可以反應(yīng)結(jié)冰累積速度。

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

將式(1)對(duì)應(yīng)的飛行動(dòng)力學(xué)系統(tǒng)寫(xiě)成參數(shù)化形式[10]，即

(7)

(8)

(9)

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

(10)

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

(11)

(12)

(13)

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

(14)

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

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

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

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

表3 測(cè)量噪聲標(biāo)準(zhǔn)差[21]Table 3 Standard deviations of measurement noises[21]

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

表4 雙水獺飛機(jī)結(jié)冰預(yù)報(bào)延時(shí)[11]Table 4 Delays in ice detection of Twin Otter Airplane[11]

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

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

表5 矩陣Q0對(duì)氣動(dòng)導(dǎo)數(shù)辨識(shí)精度的影響

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

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

表6 辨識(shí)結(jié)果誤差及結(jié)冰預(yù)報(bào)延時(shí)Table 6 Identification results errors and ice detection delays

因此，對(duì)于結(jié)冰研究樣機(jī)機(jī)翼結(jié)冰情況，通過(guò)3個(gè)氣動(dòng)導(dǎo)數(shù)的辨識(shí)結(jié)果均可對(duì)結(jié)冰嚴(yán)重性程度進(jìn)行預(yù)報(bào)，但三者的辨識(shí)精度和預(yù)報(bào)延遲差別較大，可以從中選擇預(yù)報(bào)效果最好的參數(shù)辨識(shí)結(jié)果作為結(jié)冰在線(xiàn)探測(cè)的依據(jù)。

3 飛機(jī)機(jī)翼結(jié)冰探測(cè)精度

為了評(píng)估H∞算法在不同工況下的結(jié)冰探測(cè)能力，下面通過(guò)仿真分析不同結(jié)冰累積過(guò)程以及測(cè)量噪聲大小對(duì)辨識(shí)結(jié)果的影響。

3.1 結(jié)冰累積過(guò)程影響分析

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

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

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

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

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

3.2 隨機(jī)噪聲影響下的氣動(dòng)導(dǎo)數(shù)辨識(shí)精度

由于實(shí)際飛行過(guò)程中測(cè)量噪聲的存在，必然會(huì)影響結(jié)冰氣動(dòng)導(dǎo)數(shù)辨識(shí)精度，為了對(duì)實(shí)際飛行過(guò)程中基于H∞的結(jié)冰氣動(dòng)導(dǎo)數(shù)辨識(shí)精度進(jìn)行定量評(píng)估，下面在給定測(cè)量噪聲標(biāo)準(zhǔn)差大小范圍下進(jìn)行多次蒙特卡羅仿真，分析結(jié)冰氣動(dòng)導(dǎo)數(shù)的辨識(shí)精度。

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

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

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

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

表10 不同結(jié)冰水平下氣動(dòng)導(dǎo)數(shù)辨識(shí)延時(shí)統(tǒng)計(jì)Table 10 Delays in aerodynamic derivatives identification in different icing levels

統(tǒng)計(jì)標(biāo)準(zhǔn)差隨結(jié)冰水平增大而增加，最大不超過(guò)2.5 s，Cmα的跟蹤延時(shí)均值最大不超過(guò)9.5 s，統(tǒng)計(jì)標(biāo)準(zhǔn)差最大為2.8 s左右。

4 結(jié) 論

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

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

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

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

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

參 考 文 獻(xiàn)

[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] 占榮輝，張軍. 非線(xiàn)性濾波理論與目標(biāo)跟蹤應(yīng)用[M]. 北京：國(guó)防工業(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)估計(jì):卡爾曼,H∞及非線(xiàn)性濾波[M].張勇剛, 李寧, 奔粵陽(yáng), 譯. 北京：國(guó)防工業(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] 應(yīng)思斌, 葛彤, 艾劍良. 飛機(jī)結(jié)冰時(shí)不變參數(shù)辨識(shí)技術(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] 應(yīng)思斌, 葛彤, 艾劍良. 飛機(jī)結(jié)冰氣動(dòng)參數(shù)綜合檢測(cè)方法研究[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.