頻域響應辨識優化方法研究

李雅靜, 劉旭華, 賈曉鵬, 焦崗

(1.中國飛行試驗研究院 飛機所, 陜西 西安 710089;2.西安飛豹科技發展有限公司 設計部, 陜西 西安 710089)

0 引言

等效系統在飛機動態特性評價中是很重要的方法,對試驗數據進行時域、頻域轉換是進行低階等效系統擬配的基礎,試飛數據頻域響應分析的質量直接影響擬配效果。

使用快速傅里葉方法進行頻域分析時存在兩方面問題:第一,需要人工選擇窗口寬度[1],且由于單一窗口尺寸不能在整個頻率范圍內產生滿意的相干性,因此單個窗口寬度存在局限性;第二,頻域辨識更適用于長時間數據[2],而試飛中使用脈沖、倍脈沖等短時輸入較多,包含的頻段較窄導致頻域辨識結果較差,降低了采集數據的可用性。

本文針對上述問題提出了頻域響應辨識的優化方法,即組合窗法和頻域辨識修正法。仿真和真實飛行數據的計算結果說明這兩種方法能較好地解決前述問題。

1 組合窗法

1.1 窗口選擇

首先給出兩個公式說明窗口寬度與辨識最小頻率和頻率響應中的隨機誤差的關系。

(1)窗口寬度直接決定傅里葉變換中的最小頻率fmin:

fmin=1/Twin

(1)

式中,Twin為窗口寬度。

(2)頻率響應估計中的隨機誤差計算公式:

(2)

式中,Cε為常數,用于反應窗口重疊效應,增加窗口重疊比能減小隨機誤差,代價是快速增加計算量;γxy為相干函數,飛行試驗應在風和湍流擾動最小時進行,以獲得數據最大信噪比,使相干函數最大化;nd=Trec/Twin為用于求平均的獨立時間歷程數據個數。較短的窗口Twin能增加nd,但會損害低頻范圍的辨識精度。式(1)和式(2)表明,大窗口可以增加低頻的數據信息、加寬動態范圍,而小窗口有助于抑制高頻部分的隨機誤差,所以在不同的頻率點應該選擇不同的窗口。

1.2 組合窗法

選擇單一窗口寬度需要進行折中處理,即為了提高部分頻率點的頻率響應精度,必須在其他頻率點做出犧牲。即使人工找到一個窗口寬度,也僅對所關心的頻率范圍內有限頻段具有最優性。因此需要一種方法能夠去除重復的、人工的窗口尺寸優化工作,在所關心的頻率范圍內生成一個單一的、較準確的頻率響應。為此,本文提出了組合窗法。

根據不同窗口寬度在不同頻段作用,對動作段采用4個常用窗口寬度32 s,16 s,8 s和4 s進行頻域辨識,給出窗口權值曲線如圖1所示。32 s窗口有效頻段為0.1～5.0 rad/s,16 s窗口有效頻段為1～8 rad/s,8 s窗口有效頻段為3～10 rad/s,4 s窗口有效頻段為5～50 rad/s。

圖1 窗口寬度權重Fig.1 Window size weighting coefficient

本文以幅值為例對每個頻率點進行加權,相位、相干函數類似,不再陳述。

(3)

式中,GDBC為綜合后的幅值;GDBi為相應窗口寬度對應的幅值;γi為相應窗口寬度對應的相干系數;Wi為相應窗口寬度對應的權重值。根據上述加權方法可得到最終的估計幅值、相位和相關函數。以波音747模型為例進行仿真,圖2和圖3為模型的輸入、輸出時間歷程;根據數據估計得到不同窗口寬度的頻域響應如圖4所示。

圖2 仿真輸入時間歷程Fig.2 Simulation input time history

圖3 仿真輸出時間歷程Fig.3 Simulation output time history

圖4 頻域響應Fig.4 Frequency response

圖4中r2為相干系數。由仿真結果可以看出,使用小窗口計算得到的頻率響應在低頻段相干性較差,使用大窗口計算得到的結果在高頻段隨機誤差較大,而使用組合窗法計算得到的結果在整個頻段都有較好的相干性且較為平滑。辨識結果與建模傳遞函數頻響比較如圖5所示。可以看出二者十分接近,表明使用組合窗法不需人工選擇窗口寬度,并能得到一個單一的、在整個頻率范圍內較精確的頻域響應。

圖5 真實頻響與辨識結果比較Fig.5 Comparison between true frequency response and identification result

2 頻域辨識修正法

頻域辨識方法適合長時間數據,時域辨識方法適用于短時間數據。試飛中脈沖、倍脈沖等短時間動作較多、數據較短,對時域辨識結果進行頻域分析,既可用于修正頻域辨識結果,也可用于驗證頻域辨識結果,最終能得到較好的辨識結果。

2.1 時域辨識方法

設飛機的高階增穩系統及其等效系統在相同的輸入u(t)作用下的輸出響應分別為yH(t)和yL(t),則兩個系統近似為輸入、輸出等效的條件。這里yL(t)并不是一般意義上的低階等效系統,而是一個若干階的系統。

代價函數[3]為:

(4)

yH(t)可以通過直接測量得到;yL(t)在第k個采樣時刻的數值應滿足如下差分方程:

A(q-1)yL(k)=B(q-1)u(k)+e(k)

(5)

其中:

A(q-1)=1+a1q-1+…+anq-n

(6)

B(q-1)=b1q-1+…+bnq-n+1

(7)

式中,e(k)為方程殘差,由剩余高階特性、隨機風干擾及測量噪聲等引起。

MATLAB自帶的系統辨識工具箱[4]使得辨識方法十分簡便,上述時域辨識方法只是其中的一種ARX,還包括PEM,BJ,OE,ARMAX等。時域法辨識采用多個模型結構和不同的辨識算法。使用不同辨識函數進行辨識后,比較辨識系統輸出與真實系統輸出的相似度并選出最優辨識結果,即選擇相似度最大的一組進行頻域響應計算,得到時域辨識系統的頻響。

2.2 時域辨識對頻域辨識修正法

使用時域方法得到的頻響修正直接頻域辨識結果。修正方法為:分別對數據進行時域辨識和頻域辨識,時域辨識得到狀態空間方程或傳遞函數,再對該傳遞函數進行頻域分析,得到系統頻域響應;頻域辨識直接使用快速傅里葉變換進行頻域分析;將兩個結果進行比較并加權得到最終的頻響。

首先使用桿輸入激勵由系統辨識工具箱得到的時域模型,若輸出與飛機真實輸出相似度大于85%,則說明該模型可以代替飛機真實模型;若相似度小于85%則不修正。修正量計算式為:

(8)

(9)

Wf(i)=1-Wt(i)

(10)

修正公式為:

G(i)=Wt(i)Gt(i)+Wf(i)Gf(i)

(11)

P(i)=Wt(i)Pt(i)+Wf(i)Pf(i)

(12)

式中,r2(i)為頻域辨識的相干系數;fitmax為時域辨識的相似度;Wt(i),Wf(i)分別為時域辨識頻響與頻域辨識頻響所占比例;Gt(i),Gf(i)分別為各頻率點時域辨識與頻域辨識幅值響應;Pt(i),Pf(i)分別為各頻率點時域辨識與頻域辨識相位響應。由式(9)可見，加權系數由頻域辨識的相干系數和時域辨識的相似度組成,哪個大,則響應辨識方法的結果所占權重大。

一般情況下,使用相干函數的大小評價頻域計算結果,圖6為某機俯仰角速度頻域辨識結果。直接辨識結果的低頻段相干函數小,擬配時一般適當縮小擬配頻率范圍,縮小到0.7～10.0 rad/s[5]。圖7所示的擬配結果并不好,因為該機的頻率范圍較低,從0.7 rad/s開始擬配頻率范圍較窄,導致擬配結果差。

圖6 俯仰角速度頻域辨識結果Fig.6 Frequency domain identification results of rate of pitch

圖7 時域修正前俯仰角速度擬配結果Fig.7 Matching results of rate of pitch before time domain

在低頻段辨識結果差、縮小擬配范圍也無效的情況下,使用時域結果或綜合結果可以得到較好的擬配結果,如圖8所示。大量試飛數據計算結果表明,時域辨識對頻域辨識的修正方法大大提高了試飛數據的利用率。此外時域結果可以對頻域結果進行驗證,如圖9所示。在計算過程中,如果時域估計模型的頻響與直接計算得到的頻響較接近,一般情況下可以得到很好的擬配結果。

圖8 時域修正后俯仰角速度擬配結果Fig.8 Matching results of rate of pitch after time domain

圖9 俯仰角速度頻響Fig.9 Frequency response identification result of rate of pitch

3 結束語

本文從實際出發,提出了組合窗法、時域辨識對頻域辨識修正的方法,很好地解決了頻域響應辨識過程中存在的問題,通過實例驗證了方法的有效性。將本文方法應用于實際飛機動態特性指標計算,得到了很好的效果,使辨識結果更準確、擬配結果更合理,而且大大提高了試飛數據的利用率。

參考文獻：

[1] Tichler M B,Remple R K.Aircraft and rotorcraft system identification[M].USA:American Institute of Aeronautics and Astronautics,2006:145-167.

[2] Klyde D H,Bachelder E N,Thompson P M,et al.Flying qualities parameter identification using short duration flight test inputs[R].AIAA-2007-6385,2007.

[3] 馬維金,王俊元,李鳳蘭,等.基于ARX模型的控制系統辨識及穩定性分析[J].中北大學學報,2010,31(1):9-13.

[4] 徐昕,李濤,伯曉晨.MATLAB工具箱應用指南:控制工程篇[M].北京:電子工業出版社,2000:28-41.

[5] 包立平,李春錦.評定增穩飛機飛行品質的縱向頻域等效系統法[J].北京航空學院學報,1984,(4):31-43.