于雪梅,程偉,谷偉巖
(1.北京航空航天大學航空科學與工程學院,北京 100191;2.中航通用飛機有限責任公司試飛交付中心,廣東珠海 519015)
飛行器氣動參數(shù)辨識自 Warner和 Norton[1]的早期工作以來,已經(jīng)有九十多年的歷史。以往飛行器氣動參數(shù)的確定是通過理論計算和風洞試驗進行的,而理論計算有其局限性,風洞試驗與實際飛行條件也存在差異,兩者所得到的氣動特性均難以準確反應實際飛行特性。因此,飛行試驗是確定飛機飛行性能最有效的方法,應盡可能在飛行試驗的基礎上確定飛機在使用范圍內(nèi)的飛行性能,但這樣不可避免地存在耗資大、周期長的問題。
隨著計算機和數(shù)值計算技術的發(fā)展,系統(tǒng)辨識理論在飛行器設計和研制過程中的作用越來越大,尤其是經(jīng)過20世紀60年代后期以來的研究,飛機氣動參數(shù)辨識已發(fā)展成為一個較為獨立的系統(tǒng)辨識分支[2]。極大似然法是參數(shù)辨識領域中常用的一種方法,理論上已經(jīng)證明參數(shù)的極大似然估計是一致漸進、有效、無偏的估計,且具有良好的收斂性[3-4],因而該方法在飛行器參數(shù)估計中得到廣泛應用。目前,世界發(fā)達國家對飛機狀態(tài)估計和參數(shù)辨識的研究已發(fā)展到實用水平[5-6]。我國學者在飛行器的飛行試驗參數(shù)識別方面也進行了許多研究,但總體來說,大部分研究都集中于操縱特性的氣動導數(shù)和模態(tài)特性研究上[7-8],對于飛行性能的氣動參數(shù)識別和總體特性研究并不多見。
本文結(jié)合Y12飛機雙發(fā)起飛飛行試驗實測數(shù)據(jù),研究了利用極大似然法進行飛機起飛性能參數(shù)辨識的問題,建立了Y12飛機起飛性能的數(shù)學模型,并從簡化靈敏度導數(shù)計算角度提出了以飛機在起飛滑跑過程中速度增量為觀測量的觀測方程[9],分析給出待辨識參數(shù)和利用極大似然法進行起飛參數(shù)辨識的方法,采用Matlab軟件進行了系統(tǒng)仿真[10],并針對主要辨識影響因素變化對結(jié)果的影響進行了討論。
極大似然法與以前關于飛機起飛性能的數(shù)據(jù)處理方法相比具有較高的處理效率和辨識精度,因此,將極大似然法應用于飛機起飛性能參數(shù)辨識具有一定的現(xiàn)實意義,既可提高工作效率,又可節(jié)約試驗成本。
起飛性能涉及的內(nèi)容包括:全部發(fā)動機工作正常起飛、一臺發(fā)動機失效繼續(xù)起飛、一臺發(fā)動機失效中斷起飛。本文僅研究利用極大似然法對全發(fā)起飛地面滑跑起飛性能的辨識。
從動力學角度分析,在不考慮跑道坡度的情況下,飛機在滑跑過程中承受的力有:飛機重力(G)、發(fā)動機推力(P)、氣動阻力(D)、地面摩擦力(Ff)、地面對機輪的支持力(N)、機翼升力L。飛機起飛滑跑段的運動方程可簡化為[10]:
式中,CL,CD分別為飛機起飛構(gòu)型狀態(tài)下(襟翼在起飛位置、起落架放下等)對應于停機迎角的升力系數(shù)、阻力系數(shù);f為地面滑跑摩擦系數(shù);ρ為大氣密度;SW為機翼面積。一般情況下,發(fā)動機推力P是飛行速度、高度的函數(shù),由發(fā)動機廠家提供。
飛機在滑跑過程中各力隨飛行速度的變化如圖1所示。圖中,Q為總阻力;VLOF為離地真速。在起飛離地點處有:V=VLOF,L=G。
圖1 起飛滑跑過程中各力隨速度的變化曲線
式(1)和式(2)中的速度均為飛行真速(VT)。因為飛機升力、阻力的大小只與真速有關,而起飛滑跑距離的大小卻與地速(VG)有關,因此需考慮風對起飛距離的影響。真速和地速的關系為VG=VT±VW(逆風取“-”,順風取“+”)。起飛滑跑距離可表達為d S=VGd t,積分則可獲得起飛地面滑跑距離為:
將各參數(shù)的表達式帶入式(1),并除以G,有:
式中,CL,LOF為離地瞬間的升力系數(shù)。
由式(3)可見,影響起飛滑跑距離的主要因素有:起飛質(zhì)量、襟翼狀態(tài)、跑道條件(粗糙度、場壓)、外界大氣環(huán)境(場溫、風速風向)等。起飛滑跑距離表達式中有3個未知參數(shù)(CD,CL,f),如取綜合阻力系數(shù) A=(CD-fCL),則待辨參數(shù)變?yōu)?f,A。
根據(jù)極大似然法的辨識原理,在進行參數(shù)辨識前,首先需要選擇觀測量并給出觀測方程[3]。如將式(3)作為觀測方程,將L視為觀測量,則待辨識參數(shù)f和A值位于分母中,計算中發(fā)現(xiàn),此時靈敏度導數(shù)不易計算。因此,分析后選擇起飛滑跑過程的速度增量為觀測量,以簡化觀測方程,便于進行靈敏度導數(shù)計算和參數(shù)辨識。
將起飛滑跑過程分為足夠多的小段,每一小段均可視為勻加速過程,則由式(3)可得速度增量為:
由式(8)可見,代價函數(shù)相當于測量量和觀測量的累計差值,根據(jù)模型要求,其值可以適當變化。一般計算時,當|Jj+1-Jj|/|Jj|<ε(迭代精度)時停止迭代,亦可事先定義代價函數(shù)或迭代次數(shù)進行迭代并檢查收斂情況。
起飛性能飛行試驗測試的目的就是通過準確測量飛機在起飛過程中的滑跑距離、速度、姿態(tài)、發(fā)動機功率等參數(shù)的變化規(guī)律來確定起飛性能。Y12飛機起飛性能的測試參數(shù)包括:飛機起飛質(zhì)量、起飛襟翼偏度、大氣溫度、機場場壓高度、風速、風向、指示空速、水平滑跑距離、發(fā)動機功率(包括螺旋槳轉(zhuǎn)速和扭矩)、起飛松剎車信號、前輪離地信號等。以Y12飛機雙發(fā)起飛性能為例,某一狀態(tài)(起飛質(zhì)量5 300 kg,起飛襟翼0°,大氣溫度 -14℃,機場高度145 m,逆風風速0.5 m/s)的指示空速、高度、發(fā)動機扭矩測試結(jié)果如圖2所示(由于篇幅限制,其它參數(shù)未繪出)。
圖2 飛機起飛測試結(jié)果
以上僅給出一組測試數(shù)據(jù),實際辨識過程中采用取多組測試數(shù)據(jù)進行,將辨識結(jié)果取平均值作為這一狀態(tài)飛行的待辨參數(shù)值。
采用Matlab軟件進行編程,具體辨識步驟如下:
(1)對機載測試數(shù)據(jù)(圖2)進行預處理(主要包括濾波、中值、曲線擬合),并將校正空速轉(zhuǎn)換成真速,選取合適辨識數(shù)據(jù)段,建立測試結(jié)果數(shù)據(jù)文件。
(2)根據(jù)試驗狀態(tài)、飛機參數(shù)和發(fā)動機參數(shù),計算對應試驗狀態(tài)的發(fā)動機推力和大氣密度。
(3)根據(jù)待辨識參數(shù)的理論計算結(jié)果和工作經(jīng)驗,選取待辨識參數(shù)的迭代初值。對于Y12飛機,在起飛地面滑跑過程中,有關參數(shù)的取值范圍為:
地面滑跑摩擦系數(shù):f=0.03~0.04
機翼面積:SW=34.27 m2
阻力系數(shù):CD=0.16~0.20
升力系數(shù):CL=1.2~1.4
計算得綜合阻力系數(shù):A=4.1~5.2
因此,在迭代初始,兩個待辨識參數(shù)f和A的初值應在上述范圍內(nèi)取值。
(4)將飛行速度隨時間的變化轉(zhuǎn)換為速度增量隨時間的變化并作為觀測量,建立數(shù)據(jù)文件。
(5)選擇樣本長度N=100。
(6)合理選擇代價函數(shù)和迭代精度或迭代次數(shù),按極大似然法計算待辨參數(shù)的靈敏度導數(shù)和代價函數(shù),進行迭代計算。最終使辨識參數(shù)經(jīng)過有限次數(shù)的迭代后,收斂于一個合理、可信的結(jié)果。在本文的參數(shù)辨識過程中,經(jīng)過比較,選用了人工控制迭代次數(shù)的方法。
經(jīng)上述步驟,得到辨識結(jié)果如表1所示。
表1 起飛參數(shù)辨識結(jié)果
真值是Y12飛機經(jīng)過多年飛行試驗確定的結(jié)果,比較真實可信。從表中可見:辨識結(jié)果與真值的相對誤差約為1%。將試驗時飛機的狀態(tài)參數(shù)、飛行參數(shù)和辨識值/真值分別代入式(3)中進行計算,可得到起飛滑跑距離的擴展結(jié)果分別為407 m,419 m,從工程應用角度看,該辨識結(jié)果能夠滿足性能試飛精度要求。
辨識參數(shù)迭代結(jié)果如圖3所示。
圖3 辨識參數(shù)迭代結(jié)果曲線(J=0.01)
由圖3可見,經(jīng)過30次的迭代后,辨識參數(shù)已基本趨于穩(wěn)定。迭代過程的變化趨勢說明了迭代的收斂性。由于迭代曲線縱坐標的取值范圍很小,因此,曲線在迭代到30次后雖然看似有一定的波動,但其實際變化范圍已經(jīng)很小,可以認為曲線收斂。
將地面滑跑摩擦系數(shù)f和綜合阻力系數(shù)A的迭代初值分別有意識地放大和縮小,采用同樣的數(shù)據(jù)處理和辨識方法進行辨識。結(jié)果表明,經(jīng)過40次左右的迭代后,迭代參數(shù)的波動范圍已趨正常范圍,因此,本文采用的方法對迭代初值偏差具有良好的適應性。
分別采用J=0.05,J=0.01和J=0.005三種代價函數(shù)進行對比計算,結(jié)果表明:代價函數(shù)越小,迭代次數(shù)越多,但迭代結(jié)果并沒有顯著提高,f基本在0.034 6~0.035 0之間波動,波動幅度為1%,A基本在4.550 5~4.551 3之間波動,波動幅度小于0.1%。因此,選擇合適的代價函數(shù)可以經(jīng)過較少次數(shù)的迭代就可以獲得滿意的結(jié)果。
分別采用50,100和500三種迭代次數(shù)進行計算,結(jié)果表明:迭代次數(shù)增多對迭代結(jié)果沒有顯著影響。因此,在進行參數(shù)辨識時,應選擇合適的迭代次數(shù),可以通過較少次數(shù)的迭代獲得滿意的結(jié)果。
在進行測試數(shù)據(jù)參數(shù)辨識初期,對測試數(shù)據(jù)進行預處理后,參數(shù)辨識結(jié)果很不理想,辨識結(jié)果在一個較大的范圍內(nèi)波動,有時還出現(xiàn)發(fā)散現(xiàn)象。
經(jīng)過對觀測方程式(5)的分析發(fā)現(xiàn),對Y12飛機的狀態(tài)參數(shù)和試驗數(shù)據(jù)而言,方程式(5)中的ΔV和C0兩項為同一量級,C1θ1項比ΔV低一個量級,而C2θ2U項則低至少兩個量級。在開始進行實測數(shù)據(jù)參數(shù)辨識時,觀測量ΔV的有效位數(shù)取值到10-1,這樣,其它參數(shù)的變化就被ΔV的取值誤差吃掉了。因此,觀測量ΔV和C0的精度對辨識結(jié)果具有重要影響。結(jié)果表明,將觀測量ΔV和C0的有效位數(shù)取到10-3,參數(shù)辨識結(jié)果明顯改善。
因此,觀測方程中各個參數(shù)對辨識結(jié)果的影響作用是不一樣的。在參數(shù)辨識前,應初步對觀測方程中的各個參數(shù)進行敏感性分析,以保證測試數(shù)據(jù)的選取有足夠的精度。
對Y12飛機起飛性能測試而言,由于速度傳感器有效工作范圍的限制(小速度指示誤差較大),起飛滑跑段速度的有效范圍基本在50 km/h以上,從該速度到飛機起飛離地,對應的時間歷程約為10 s左右。測試系統(tǒng)的采樣率較高(通常50~100 Hz),雖然理論上選取任意幾秒作為數(shù)據(jù)段測試數(shù)據(jù)點都足以進行參數(shù)辨識,但實際上數(shù)據(jù)段的選取對辨識結(jié)果有一定影響,應當有目的地選取測試數(shù)據(jù)段。
從觀測方程式(5)可見,式中C2θ2U項比觀測量ΔV低一個量級甚至更多,而U為起飛滑跑速度的平方,因此,為使C2θ2U項與觀測量ΔV項具有更接近的量級,試驗數(shù)據(jù)段應盡量在較大速度范圍內(nèi)選取。
極大似然法在理論上希望樣本長度越大越好,然而實際樣本過長不僅會增加計算量,而且數(shù)據(jù)中非理想誤差的增加也會影響參數(shù)辨識結(jié)果。為確定樣本長度對辨識結(jié)果的影響,本文進行了不同樣本長度的對比辨識,在同樣的測試數(shù)據(jù)段內(nèi),分別進行樣本長度N=100和N=500情況的迭代,迭代結(jié)果表明:對Y12飛機起飛性能進行辨識,在同樣迭代次數(shù)的情況下,樣本長度的增加并沒有使迭代收斂增快,在樣本長度N=100的情況下,已具有比較滿意的辨識結(jié)果,當樣本長度N=500時,辨識結(jié)果精度并沒有顯著提高,而迭代計算需要的時間卻顯著增加。
本文研究了利用極大似然法進行飛機起飛性能飛行試驗數(shù)據(jù)參數(shù)辨識的問題,確立了辨識模型和待辨參數(shù),采用Matlab軟件進行計算,結(jié)果表明極大似然法可以快速、準確地辨識出起飛性能所需的參數(shù),并具有較高的精度。從工程實用的角度出發(fā),辨識中所使用的觀測方程應具有簡單的形式,以便于靈敏度導數(shù)的計算。在觀測方程中,各參數(shù)對辨識結(jié)果的影響是不一樣的,為保證辨識結(jié)果的準確性,應進行敏感參數(shù)分析,確保測試數(shù)據(jù)的記錄具有足夠的精度。另外,代價函數(shù)、迭代次數(shù)、樣本長度、數(shù)據(jù)段的選取等對參數(shù)辨識結(jié)果精度都可能有影響,在參數(shù)辨識過程中應注意分析,以盡可能經(jīng)過較少迭代獲得最優(yōu)結(jié)果。與傳統(tǒng)方法相比,極大似然法具有數(shù)據(jù)處理效率高、辨識結(jié)果精度高的特點,在今后的飛機起飛性能數(shù)據(jù)分析中應進一步推廣,并應探索將其應用到其它性能分析中的方法。
[1] Warner E P,Norton F H.Preliminary report on free flight test[R].NASA Report No 70,1919.
[2] 蔡金獅.飛行器系統(tǒng)辨識學[M].北京:國防工業(yè)出版社,2002:287-294.
[3] 葉建華.過程辨識技術[M].上海:上海大學出版社,2007:59-82.
[4] 王志賢.最優(yōu)狀態(tài)估計與系統(tǒng)辨識[M].西安:西北工業(yè)大學出版社,2004:223-243.
[5] Klein V,Mogan DR.Estimation ofbias errors inmeasured airplane responses using maximum likelihood method[R].NASA TM-89059,1987:79-82.
[6] Maine R E.Programmer’smanual for MMLE3,a general FORTRAN program formaximum likelihood parameter estimation[R].NASA TP-1690,1981.
[7] 全昌業(yè).氣動參數(shù)辨識飛機飛行試驗研究回顧[J].飛行試驗,1999,15(2):8-14.
[8] 劉興堂,萬少松.一種估計飛行器運動特性的實用極大似然法[J].飛行力學,1999,17(3):65-70.
[9] 谷偉巖.運十二飛機起飛性能參數(shù)辨識研究[D].北京:北京航空航天大學,2000.
[10] 張靜.Matlab在控制系統(tǒng)中的應用[M].北京:電子工業(yè)出版社,2008:364-83.