于雪梅，程偉，谷偉巖

(1.北京航空航天大學(xué)航空科學(xué)與工程學(xué)院，北京 100191;2.中航通用飛機(jī)有限責(zé)任公司試飛交付中心，廣東珠海 519015)

引言

飛行器氣動(dòng)參數(shù)辨識(shí)自 Warner和 Norton［1］的早期工作以來，已經(jīng)有九十多年的歷史。以往飛行器氣動(dòng)參數(shù)的確定是通過理論計(jì)算和風(fēng)洞試驗(yàn)進(jìn)行的，而理論計(jì)算有其局限性，風(fēng)洞試驗(yàn)與實(shí)際飛行條件也存在差異，兩者所得到的氣動(dòng)特性均難以準(zhǔn)確反應(yīng)實(shí)際飛行特性。因此，飛行試驗(yàn)是確定飛機(jī)飛行性能最有效的方法，應(yīng)盡可能在飛行試驗(yàn)的基礎(chǔ)上確定飛機(jī)在使用范圍內(nèi)的飛行性能，但這樣不可避免地存在耗資大、周期長的問題。

隨著計(jì)算機(jī)和數(shù)值計(jì)算技術(shù)的發(fā)展，系統(tǒng)辨識(shí)理論在飛行器設(shè)計(jì)和研制過程中的作用越來越大，尤其是經(jīng)過20世紀(jì)60年代后期以來的研究，飛機(jī)氣動(dòng)參數(shù)辨識(shí)已發(fā)展成為一個(gè)較為獨(dú)立的系統(tǒng)辨識(shí)分支［2］。極大似然法是參數(shù)辨識(shí)領(lǐng)域中常用的一種方法，理論上已經(jīng)證明參數(shù)的極大似然估計(jì)是一致漸進(jìn)、有效、無偏的估計(jì)，且具有良好的收斂性［3-4］，因而該方法在飛行器參數(shù)估計(jì)中得到廣泛應(yīng)用。目前，世界發(fā)達(dá)國家對飛機(jī)狀態(tài)估計(jì)和參數(shù)辨識(shí)的研究已發(fā)展到實(shí)用水平［5-6］。我國學(xué)者在飛行器的飛行試驗(yàn)參數(shù)識(shí)別方面也進(jìn)行了許多研究，但總體來說，大部分研究都集中于操縱特性的氣動(dòng)導(dǎo)數(shù)和模態(tài)特性研究上［7-8］，對于飛行性能的氣動(dòng)參數(shù)識(shí)別和總體特性研究并不多見。

本文結(jié)合Y12飛機(jī)雙發(fā)起飛飛行試驗(yàn)實(shí)測數(shù)據(jù)，研究了利用極大似然法進(jìn)行飛機(jī)起飛性能參數(shù)辨識(shí)的問題，建立了Y12飛機(jī)起飛性能的數(shù)學(xué)模型，并從簡化靈敏度導(dǎo)數(shù)計(jì)算角度提出了以飛機(jī)在起飛滑跑過程中速度增量為觀測量的觀測方程［9］，分析給出待辨識(shí)參數(shù)和利用極大似然法進(jìn)行起飛參數(shù)辨識(shí)的方法，采用Matlab軟件進(jìn)行了系統(tǒng)仿真［10］，并針對主要辨識(shí)影響因素變化對結(jié)果的影響進(jìn)行了討論。

極大似然法與以前關(guān)于飛機(jī)起飛性能的數(shù)據(jù)處理方法相比具有較高的處理效率和辨識(shí)精度，因此，將極大似然法應(yīng)用于飛機(jī)起飛性能參數(shù)辨識(shí)具有一定的現(xiàn)實(shí)意義，既可提高工作效率，又可節(jié)約試驗(yàn)成本。

1 飛機(jī)起飛性能建模

起飛性能涉及的內(nèi)容包括:全部發(fā)動(dòng)機(jī)工作正常起飛、一臺(tái)發(fā)動(dòng)機(jī)失效繼續(xù)起飛、一臺(tái)發(fā)動(dòng)機(jī)失效中斷起飛。本文僅研究利用極大似然法對全發(fā)起飛地面滑跑起飛性能的辨識(shí)。

1.1 起飛性能數(shù)學(xué)模型

從動(dòng)力學(xué)角度分析，在不考慮跑道坡度的情況下，飛機(jī)在滑跑過程中承受的力有:飛機(jī)重力(G)、發(fā)動(dòng)機(jī)推力(P)、氣動(dòng)阻力(D)、地面摩擦力(Ff)、地面對機(jī)輪的支持力(N)、機(jī)翼升力L。飛機(jī)起飛滑跑段的運(yùn)動(dòng)方程可簡化為［10］:

式中，CL，CD分別為飛機(jī)起飛構(gòu)型狀態(tài)下(襟翼在起飛位置、起落架放下等)對應(yīng)于停機(jī)迎角的升力系數(shù)、阻力系數(shù);f為地面滑跑摩擦系數(shù);ρ為大氣密度;SW為機(jī)翼面積。一般情況下，發(fā)動(dòng)機(jī)推力P是飛行速度、高度的函數(shù)，由發(fā)動(dòng)機(jī)廠家提供。

飛機(jī)在滑跑過程中各力隨飛行速度的變化如圖1所示。圖中，Q為總阻力;VLOF為離地真速。在起飛離地點(diǎn)處有:V=VLOF，L=G。

圖1 起飛滑跑過程中各力隨速度的變化曲線

式(1)和式(2)中的速度均為飛行真速(VT)。因?yàn)轱w機(jī)升力、阻力的大小只與真速有關(guān)，而起飛滑跑距離的大小卻與地速(VG)有關(guān)，因此需考慮風(fēng)對起飛距離的影響。真速和地速的關(guān)系為VG=VT±VW(逆風(fēng)取“－”，順風(fēng)取“+”)。起飛滑跑距離可表達(dá)為d S=VGd t，積分則可獲得起飛地面滑跑距離為:

將各參數(shù)的表達(dá)式帶入式(1)，并除以G，有:

式中，CL，LOF為離地瞬間的升力系數(shù)。

由式(3)可見，影響起飛滑跑距離的主要因素有:起飛質(zhì)量、襟翼狀態(tài)、跑道條件(粗糙度、場壓)、外界大氣環(huán)境(場溫、風(fēng)速風(fēng)向)等。起飛滑跑距離表達(dá)式中有3個(gè)未知參數(shù)(CD，CL，f)，如取綜合阻力系數(shù) A=(CD－fCL)，則待辨參數(shù)變?yōu)?f，A。

1.2 起飛性能辨識(shí)模型及辨識(shí)方法

根據(jù)極大似然法的辨識(shí)原理，在進(jìn)行參數(shù)辨識(shí)前，首先需要選擇觀測量并給出觀測方程［3］。如將式(3)作為觀測方程，將L視為觀測量，則待辨識(shí)參數(shù)f和A值位于分母中，計(jì)算中發(fā)現(xiàn)，此時(shí)靈敏度導(dǎo)數(shù)不易計(jì)算。因此，分析后選擇起飛滑跑過程的速度增量為觀測量，以簡化觀測方程，便于進(jìn)行靈敏度導(dǎo)數(shù)計(jì)算和參數(shù)辨識(shí)。

將起飛滑跑過程分為足夠多的小段，每一小段均可視為勻加速過程，則由式(3)可得速度增量為:

由式(8)可見，代價(jià)函數(shù)相當(dāng)于測量量和觀測量的累計(jì)差值，根據(jù)模型要求，其值可以適當(dāng)變化。一般計(jì)算時(shí)，當(dāng)|Jj+1－Jj|/|Jj|＜ε(迭代精度)時(shí)停止迭代，亦可事先定義代價(jià)函數(shù)或迭代次數(shù)進(jìn)行迭代并檢查收斂情況。

2 起飛性能飛行試驗(yàn)及數(shù)據(jù)采集

起飛性能飛行試驗(yàn)測試的目的就是通過準(zhǔn)確測量飛機(jī)在起飛過程中的滑跑距離、速度、姿態(tài)、發(fā)動(dòng)機(jī)功率等參數(shù)的變化規(guī)律來確定起飛性能。Y12飛機(jī)起飛性能的測試參數(shù)包括:飛機(jī)起飛質(zhì)量、起飛襟翼偏度、大氣溫度、機(jī)場場壓高度、風(fēng)速、風(fēng)向、指示空速、水平滑跑距離、發(fā)動(dòng)機(jī)功率(包括螺旋槳轉(zhuǎn)速和扭矩)、起飛松剎車信號(hào)、前輪離地信號(hào)等。以Y12飛機(jī)雙發(fā)起飛性能為例，某一狀態(tài)(起飛質(zhì)量5 300 kg，起飛襟翼0°，大氣溫度 －14℃，機(jī)場高度145 m，逆風(fēng)風(fēng)速0.5 m/s)的指示空速、高度、發(fā)動(dòng)機(jī)扭矩測試結(jié)果如圖2所示(由于篇幅限制，其它參數(shù)未繪出)。

圖2 飛機(jī)起飛測試結(jié)果

以上僅給出一組測試數(shù)據(jù)，實(shí)際辨識(shí)過程中采用取多組測試數(shù)據(jù)進(jìn)行，將辨識(shí)結(jié)果取平均值作為這一狀態(tài)飛行的待辨參數(shù)值。

3 極大似然法參數(shù)辨識(shí)

采用Matlab軟件進(jìn)行編程，具體辨識(shí)步驟如下:

(1)對機(jī)載測試數(shù)據(jù)(圖2)進(jìn)行預(yù)處理(主要包括濾波、中值、曲線擬合)，并將校正空速轉(zhuǎn)換成真速，選取合適辨識(shí)數(shù)據(jù)段，建立測試結(jié)果數(shù)據(jù)文件。

(2)根據(jù)試驗(yàn)狀態(tài)、飛機(jī)參數(shù)和發(fā)動(dòng)機(jī)參數(shù)，計(jì)算對應(yīng)試驗(yàn)狀態(tài)的發(fā)動(dòng)機(jī)推力和大氣密度。

(3)根據(jù)待辨識(shí)參數(shù)的理論計(jì)算結(jié)果和工作經(jīng)驗(yàn)，選取待辨識(shí)參數(shù)的迭代初值。對于Y12飛機(jī)，在起飛地面滑跑過程中，有關(guān)參數(shù)的取值范圍為:

地面滑跑摩擦系數(shù):f=0.03～0.04

機(jī)翼面積:SW=34.27 m2

阻力系數(shù):CD=0.16～0.20

升力系數(shù):CL=1.2～1.4

計(jì)算得綜合阻力系數(shù):A=4.1～5.2

因此，在迭代初始，兩個(gè)待辨識(shí)參數(shù)f和A的初值應(yīng)在上述范圍內(nèi)取值。

(4)將飛行速度隨時(shí)間的變化轉(zhuǎn)換為速度增量隨時(shí)間的變化并作為觀測量，建立數(shù)據(jù)文件。

(5)選擇樣本長度N=100。

(6)合理選擇代價(jià)函數(shù)和迭代精度或迭代次數(shù)，按極大似然法計(jì)算待辨參數(shù)的靈敏度導(dǎo)數(shù)和代價(jià)函數(shù)，進(jìn)行迭代計(jì)算。最終使辨識(shí)參數(shù)經(jīng)過有限次數(shù)的迭代后，收斂于一個(gè)合理、可信的結(jié)果。在本文的參數(shù)辨識(shí)過程中，經(jīng)過比較，選用了人工控制迭代次數(shù)的方法。

經(jīng)上述步驟，得到辨識(shí)結(jié)果如表1所示。

表1 起飛參數(shù)辨識(shí)結(jié)果

真值是Y12飛機(jī)經(jīng)過多年飛行試驗(yàn)確定的結(jié)果，比較真實(shí)可信。從表中可見:辨識(shí)結(jié)果與真值的相對誤差約為1%。將試驗(yàn)時(shí)飛機(jī)的狀態(tài)參數(shù)、飛行參數(shù)和辨識(shí)值/真值分別代入式(3)中進(jìn)行計(jì)算，可得到起飛滑跑距離的擴(kuò)展結(jié)果分別為407 m，419 m，從工程應(yīng)用角度看，該辨識(shí)結(jié)果能夠滿足性能試飛精度要求。

辨識(shí)參數(shù)迭代結(jié)果如圖3所示。

圖3 辨識(shí)參數(shù)迭代結(jié)果曲線(J=0.01)

由圖3可見，經(jīng)過30次的迭代后，辨識(shí)參數(shù)已基本趨于穩(wěn)定。迭代過程的變化趨勢說明了迭代的收斂性。由于迭代曲線縱坐標(biāo)的取值范圍很小，因此，曲線在迭代到30次后雖然看似有一定的波動(dòng)，但其實(shí)際變化范圍已經(jīng)很小，可以認(rèn)為曲線收斂。

4 辨識(shí)結(jié)果影響因素分析

4.1 迭代初值影響

將地面滑跑摩擦系數(shù)f和綜合阻力系數(shù)A的迭代初值分別有意識(shí)地放大和縮小，采用同樣的數(shù)據(jù)處理和辨識(shí)方法進(jìn)行辨識(shí)。結(jié)果表明，經(jīng)過40次左右的迭代后，迭代參數(shù)的波動(dòng)范圍已趨正常范圍，因此，本文采用的方法對迭代初值偏差具有良好的適應(yīng)性。

4.2 代價(jià)函數(shù)影響

分別采用J=0.05，J=0.01和J=0.005三種代價(jià)函數(shù)進(jìn)行對比計(jì)算，結(jié)果表明:代價(jià)函數(shù)越小，迭代次數(shù)越多，但迭代結(jié)果并沒有顯著提高，f基本在0.034 6～0.035 0之間波動(dòng)，波動(dòng)幅度為1%，A基本在4.550 5～4.551 3之間波動(dòng)，波動(dòng)幅度小于0.1%。因此，選擇合適的代價(jià)函數(shù)可以經(jīng)過較少次數(shù)的迭代就可以獲得滿意的結(jié)果。

4.3 迭代次數(shù)影響

分別采用50，100和500三種迭代次數(shù)進(jìn)行計(jì)算，結(jié)果表明:迭代次數(shù)增多對迭代結(jié)果沒有顯著影響。因此，在進(jìn)行參數(shù)辨識(shí)時(shí)，應(yīng)選擇合適的迭代次數(shù)，可以通過較少次數(shù)的迭代獲得滿意的結(jié)果。

4.4 敏感參數(shù)影響

在進(jìn)行測試數(shù)據(jù)參數(shù)辨識(shí)初期，對測試數(shù)據(jù)進(jìn)行預(yù)處理后，參數(shù)辨識(shí)結(jié)果很不理想，辨識(shí)結(jié)果在一個(gè)較大的范圍內(nèi)波動(dòng)，有時(shí)還出現(xiàn)發(fā)散現(xiàn)象。

經(jīng)過對觀測方程式(5)的分析發(fā)現(xiàn)，對Y12飛機(jī)的狀態(tài)參數(shù)和試驗(yàn)數(shù)據(jù)而言，方程式(5)中的ΔV和C0兩項(xiàng)為同一量級(jí)，C1θ1項(xiàng)比ΔV低一個(gè)量級(jí)，而C2θ2U項(xiàng)則低至少兩個(gè)量級(jí)。在開始進(jìn)行實(shí)測數(shù)據(jù)參數(shù)辨識(shí)時(shí)，觀測量ΔV的有效位數(shù)取值到10－1，這樣，其它參數(shù)的變化就被ΔV的取值誤差吃掉了。因此，觀測量ΔV和C0的精度對辨識(shí)結(jié)果具有重要影響。結(jié)果表明，將觀測量ΔV和C0的有效位數(shù)取到10－3，參數(shù)辨識(shí)結(jié)果明顯改善。

因此，觀測方程中各個(gè)參數(shù)對辨識(shí)結(jié)果的影響作用是不一樣的。在參數(shù)辨識(shí)前，應(yīng)初步對觀測方程中的各個(gè)參數(shù)進(jìn)行敏感性分析，以保證測試數(shù)據(jù)的選取有足夠的精度。

4.5 試驗(yàn)數(shù)據(jù)段的選取

對Y12飛機(jī)起飛性能測試而言，由于速度傳感器有效工作范圍的限制(小速度指示誤差較大)，起飛滑跑段速度的有效范圍基本在50 km/h以上，從該速度到飛機(jī)起飛離地，對應(yīng)的時(shí)間歷程約為10 s左右。測試系統(tǒng)的采樣率較高(通常50～100 Hz)，雖然理論上選取任意幾秒作為數(shù)據(jù)段測試數(shù)據(jù)點(diǎn)都足以進(jìn)行參數(shù)辨識(shí)，但實(shí)際上數(shù)據(jù)段的選取對辨識(shí)結(jié)果有一定影響，應(yīng)當(dāng)有目的地選取測試數(shù)據(jù)段。

從觀測方程式(5)可見，式中C2θ2U項(xiàng)比觀測量ΔV低一個(gè)量級(jí)甚至更多，而U為起飛滑跑速度的平方，因此，為使C2θ2U項(xiàng)與觀測量ΔV項(xiàng)具有更接近的量級(jí)，試驗(yàn)數(shù)據(jù)段應(yīng)盡量在較大速度范圍內(nèi)選取。

4.6 樣本長度影響分析

極大似然法在理論上希望樣本長度越大越好，然而實(shí)際樣本過長不僅會(huì)增加計(jì)算量，而且數(shù)據(jù)中非理想誤差的增加也會(huì)影響參數(shù)辨識(shí)結(jié)果。為確定樣本長度對辨識(shí)結(jié)果的影響，本文進(jìn)行了不同樣本長度的對比辨識(shí)，在同樣的測試數(shù)據(jù)段內(nèi)，分別進(jìn)行樣本長度N=100和N=500情況的迭代，迭代結(jié)果表明:對Y12飛機(jī)起飛性能進(jìn)行辨識(shí)，在同樣迭代次數(shù)的情況下，樣本長度的增加并沒有使迭代收斂增快，在樣本長度N=100的情況下，已具有比較滿意的辨識(shí)結(jié)果，當(dāng)樣本長度N=500時(shí)，辨識(shí)結(jié)果精度并沒有顯著提高，而迭代計(jì)算需要的時(shí)間卻顯著增加。

5 結(jié)束語

本文研究了利用極大似然法進(jìn)行飛機(jī)起飛性能飛行試驗(yàn)數(shù)據(jù)參數(shù)辨識(shí)的問題，確立了辨識(shí)模型和待辨參數(shù)，采用Matlab軟件進(jìn)行計(jì)算，結(jié)果表明極大似然法可以快速、準(zhǔn)確地辨識(shí)出起飛性能所需的參數(shù)，并具有較高的精度。從工程實(shí)用的角度出發(fā)，辨識(shí)中所使用的觀測方程應(yīng)具有簡單的形式，以便于靈敏度導(dǎo)數(shù)的計(jì)算。在觀測方程中，各參數(shù)對辨識(shí)結(jié)果的影響是不一樣的，為保證辨識(shí)結(jié)果的準(zhǔn)確性，應(yīng)進(jìn)行敏感參數(shù)分析，確保測試數(shù)據(jù)的記錄具有足夠的精度。另外，代價(jià)函數(shù)、迭代次數(shù)、樣本長度、數(shù)據(jù)段的選取等對參數(shù)辨識(shí)結(jié)果精度都可能有影響，在參數(shù)辨識(shí)過程中應(yīng)注意分析，以盡可能經(jīng)過較少迭代獲得最優(yōu)結(jié)果。與傳統(tǒng)方法相比，極大似然法具有數(shù)據(jù)處理效率高、辨識(shí)結(jié)果精度高的特點(diǎn)，在今后的飛機(jī)起飛性能數(shù)據(jù)分析中應(yīng)進(jìn)一步推廣，并應(yīng)探索將其應(yīng)用到其它性能分析中的方法。

［1］ Warner E P，Norton F H.Preliminary report on free flight test［R］.NASA Report No 70，1919.

［2］ 蔡金獅.飛行器系統(tǒng)辨識(shí)學(xué)［M］.北京:國防工業(yè)出版社，2002:287-294.

［3］ 葉建華.過程辨識(shí)技術(shù)［M］.上海:上海大學(xué)出版社，2007:59-82.

［4］ 王志賢.最優(yōu)狀態(tài)估計(jì)與系統(tǒng)辨識(shí)［M］.西安:西北工業(yè)大學(xué)出版社，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è).氣動(dòng)參數(shù)辨識(shí)飛機(jī)飛行試驗(yàn)研究回顧［J］.飛行試驗(yàn)，1999，15(2):8-14.

［8］ 劉興堂，萬少松.一種估計(jì)飛行器運(yùn)動(dòng)特性的實(shí)用極大似然法［J］.飛行力學(xué)，1999，17(3):65-70.

［9］ 谷偉巖.運(yùn)十二飛機(jī)起飛性能參數(shù)辨識(shí)研究［D］.北京:北京航空航天大學(xué)，2000.

［10］ 張靜.Matlab在控制系統(tǒng)中的應(yīng)用［M］.北京:電子工業(yè)出版社，2008:364-83.