，，

(中國飛行試驗研究院 飛機所，西安 710089)

0 引言

電傳操縱系統(tǒng)現(xiàn)在已廣泛應(yīng)用于新型軍民用飛機,極大地改變了飛機動力學(xué)系統(tǒng)的特性。電傳飛機飛行品質(zhì)的評價和驗證面臨著新的技術(shù)挑戰(zhàn)。

飛行品質(zhì)規(guī)范[1]普遍建議使用低階等效系統(tǒng)的方法評定電傳飛機的飛行品質(zhì)。用具有固定形式的低階等效系統(tǒng)代替高階系統(tǒng)，使用優(yōu)化方法給出低階等效系統(tǒng)的參數(shù)，計算飛行品質(zhì)指標，評定飛行品質(zhì)等級。[2-3]采用低階等效系統(tǒng)的方法評定電傳飛機飛行品質(zhì)必須確定電傳飛機的低階等效系統(tǒng)模型，然后尋找合理的尋優(yōu)算法，保證低階等效系統(tǒng)能夠合理有效地反應(yīng)電傳飛機動力學(xué)系統(tǒng)的特性。

對于早期沒有電傳操縱系統(tǒng)的飛機，由于飛機動力學(xué)系統(tǒng)的輸入即為飛行員的操縱，飛機舵面與飛行員操縱基本呈線性關(guān)系，從飛行動力學(xué)方程簡化得到的小擾動方程即可反應(yīng)飛機動力學(xué)系統(tǒng)的模型。小擾動方程為典型的線性系統(tǒng)，模型結(jié)構(gòu)簡單，沒有時間延遲，辨識方法較為成熟，國內(nèi)已經(jīng)傳統(tǒng)機械操縱飛機的飛行品質(zhì)評定中總結(jié)出了有效的辨識算法。但是，電傳操縱系統(tǒng)的使用使得飛行員的操縱輸入與飛機舵面不再是簡單的線性關(guān)系。電傳操縱系統(tǒng)的引入使得飛機動力學(xué)系統(tǒng)不再是簡單的線性系統(tǒng)，反饋控制、濾波器等控制環(huán)節(jié)的引入極大地增加了飛機動力學(xué)系統(tǒng)的階次，時間延遲的引入也增加了系統(tǒng)的非線性。傳統(tǒng)的辨識算法已不能夠滿足電傳飛機飛行品質(zhì)評定的需求，存在時間延遲干擾大、參數(shù)初值敏感性高、辨識結(jié)果可靠性低等問題。

本文在飛機小擾動方程線性系統(tǒng)模型的基礎(chǔ)上，添加時間延遲項，建立以氣動導(dǎo)數(shù)和時間常數(shù)為待辨識參數(shù)的飛機動力學(xué)系統(tǒng)傳遞函數(shù)模型，以某型電傳飛機飛行試驗數(shù)據(jù)為模型輸入和輸出，對高階系統(tǒng)進行頻域分析，綜合運用方程誤差法、輸出誤差法和單純形法等多種頻域?qū)?yōu)算法，辨識飛機的低階等效系統(tǒng)模型參數(shù)，計算電傳飛機飛行品質(zhì)指標，評定飛行品質(zhì)等級，為飛行品質(zhì)評定提供技術(shù)支持。

1 低階等效系統(tǒng)模型

電傳飛機縱向動態(tài)特性包括長周期和短周期模態(tài)，橫航向動態(tài)特性包括荷蘭滾、滾轉(zhuǎn)和螺旋3種模態(tài)。其中縱向長周期模態(tài)和橫航向螺旋模態(tài)均為弱模態(tài)，本文主要研究縱向短周期和橫航向荷蘭滾、滾轉(zhuǎn)等復(fù)雜模態(tài)。

根據(jù)小擾動理論，飛機的縱向短周期模態(tài)可以用式(1)表示[6][12],其中，ɑ為迎角，q為俯仰角速率，Nz為法向過載，δe為飛機的縱向桿位移。由于迎角的測量誤差較大，在飛機的縱向短周期模型中，將迎角作為系統(tǒng)的狀態(tài)量，將俯仰角速率和法向過載作為系統(tǒng)的輸出量，縱向桿位移為系統(tǒng)的輸入量。將式(1)所示的狀態(tài)空間模型進行拉普拉斯變換，可以獲得飛機縱向短周期模態(tài)的傳遞函數(shù)模型，以[a0，a1，b0，b1，c0，c1]代替復(fù)雜的拉普拉斯變量系數(shù)，增加時間延遲項，可以得到如式(2)所示的縱向短周期模態(tài)傳遞函數(shù)模型。

(1)

(2)

同理，飛機的橫航向小擾動方程可以用式(3)表示[2-3]，其中β為側(cè)滑角，p為滾轉(zhuǎn)角速率，r為偏航角速率，Ny為側(cè)向過載，δa和δr分別為橫向桿位移和腳蹬位移。對式(3)進行相似的處理，以滾轉(zhuǎn)角速率和偏航角速率為系統(tǒng)輸出，可以得到如式(4)和(5)的荷蘭滾、滾轉(zhuǎn)模態(tài)的傳遞函數(shù)模型。

(3)

(4)

(5)

基于以上傳遞函數(shù)模型計算其特征根，從而得到飛機的縱向短周期無阻尼自然頻率、阻尼比，荷蘭滾無阻尼自然頻率、阻尼比，滾轉(zhuǎn)模態(tài)時間常數(shù)等飛行品質(zhì)指標。用特征根計算飛行品質(zhì)指標的計算方法詳見參考文獻6，本文不再贅述。

2 頻域相干性分析

頻域參數(shù)辨識方法需要對系統(tǒng)的輸入輸出進行頻譜分析，而相關(guān)函數(shù)則是一種重要的檢查輸入輸出因果關(guān)系的方法。相干性分析可以為頻域參數(shù)辨識方法確定辨識的頻率范圍。對于已知的系統(tǒng)輸入和系統(tǒng)輸出，根據(jù)快速傅立葉變換(FFT)方法，可以得到輸入輸出的頻譜[4]：

X(f)=FFT(x),Y(f)=FFT(y)

(6)

自相關(guān)函數(shù)定義如下：

(7)

其中：T表示記錄的時間段長度，互相關(guān)函數(shù)定義如下：

(8)

根據(jù)以上公式，可以得到各個頻率點的相干系數(shù)定義如下：

(9)

相干系數(shù)分析可以確定辨識的頻譜信息范圍，以下列輸入輸出信息為例，如圖1～2所示。

圖1 輸入輸出時間歷程

圖2 相關(guān)系數(shù)

根據(jù)圖1和圖2所示，在2～15 rad/s的頻率范圍內(nèi)，輸入輸出的相關(guān)度較高，而飛行品質(zhì)關(guān)心的頻率范圍為0.1～10 rad/s，因此確定頻域參數(shù)辨識的頻率范圍為2～10 rad/s。

3 高精度有限傅立葉變換

高精度有限傅立葉變換[5]可以對有限長度的時域信號在指定的頻率范圍內(nèi)進行分解[2]。結(jié)合本文第二節(jié)的相關(guān)性分析結(jié)果，可以得到指定頻率范圍內(nèi)的精確的頻譜信息。

假設(shè)所選擇的頻率范圍為f∈[f0,f1)，并且取定M個離散頻率點

fk=f0+kΔf,k=0,1,2,…,M-1

(10)

其中：f=(f1-f0)/M。因此離散傅立葉變換可以寫為如下形式：

k=0,1,2,…,M-1

(11)

定義

Φ0=2πf0Δt,ΔΦ=2πΔfΔt

A=ejΦ0,Z=ejΔΦ

代入上式可以得到：

(12)

隨著k的增加，AZk表示在Z平面上單位矢量畫出的輪廓線。Ф0表示單位矢量的初始位置的角度和位置，和初始頻率f0有關(guān)；△ Ф表示沿著單位圓隨著頻率的增加而導(dǎo)致的角度和位置增量，和△f與k有關(guān)。對以上形式使用參考文獻[2]中的修正方法，即可得到準確的高精度有限傅立葉變換。

4 參數(shù)辨識算法

本文使用頻域方法計算低階等效系統(tǒng)的參數(shù)值。根據(jù)試驗數(shù)據(jù)的高精度有限傅立葉變換得到的頻譜信息，利用頻域方程誤差法和單純形法可以得到合理的待辨識參數(shù)的近似值，進一步使用輸出誤差法進行迭代，可以得到最終的參數(shù)值。

根據(jù)時域極大似然法原理[6]，對于飛行器動力學(xué)系統(tǒng)而言，極大似然函數(shù)形式：

(13)

其中:R為加權(quán)矩陣，無特殊要求時，無差別的定義R為單位陣。根據(jù)傅立葉變換的帕斯瓦爾定理，以縱向短周期模型為例，可以定義頻域極大似然函數(shù)如下：

(14)

其中:w為加權(quán)系數(shù)，用以平衡似然函數(shù)中俯仰角速率與法向過載之間的權(quán)重。一般情況下，我們定義加權(quán)系數(shù)如式(14)所示。加權(quán)系數(shù)可以根據(jù)擬配計算的實際情況進行調(diào)整以獲取更好的擬配效果。

(15)

4.1 方程誤差法

頻域方程誤差法計算方法簡單，可以為輸出誤差法提供迭代待初始值。以縱向單擬配為例，將式(2)縱向短周期傳遞函數(shù)模型中的拉普拉斯變量s替換為jω，改寫為：

(16)

(17)

(18)

代價函數(shù)為:

(19)

?θ=(XXT)-1XTY

(20)

方程誤差法需要給出時間延遲的初始值。實際計算中，可以取輸入數(shù)據(jù)和輸出數(shù)據(jù)第一個峰值之間的時間差作為時間延遲的初值。方程誤差法可以通過簡單的一步迭代獲取待辨識參數(shù)的初步估計值，提高后續(xù)迭代計算的效率。

4.2 單純形法

單純形法在多參數(shù)優(yōu)化函數(shù)最值問題中有廣泛的應(yīng)用。單純形法的計算原理在參考文獻[4]中有非常詳細的介紹。單純形法的優(yōu)點是對于初值的選取不敏感，且迭代收斂特性好；缺點是迭代收斂速度比較慢，計算量大，對于未知參數(shù)較多的情況，單純形法往往不能得到合理的解。綜合使用單純形法和方程誤差法，可以為輸出誤差法提供非常接近最優(yōu)解的初值，為保證輸出誤差法的迭代收斂提供很好的前提。

4.3 輸出誤差法

頻域輸出誤差法是利用牛頓-拉夫遜法[7]對極大似然函數(shù)優(yōu)化的算法。極大似然函數(shù)一般是待辨識參數(shù)的函數(shù)，假設(shè)其形式為J(θ)，對J(θ)進行泰勒展開如下：

(21)

θk+1=θk+Δθ

(22)

(23)

(24)

其中:θ為待辨識模型參數(shù)向量，k表示當(dāng)前的迭代步數(shù)，M和N分別是在θk的基礎(chǔ)上計算得到的系數(shù)矩陣。式(24)即為牛頓-拉夫遜法給出的待辨識參數(shù)迭代算法。其中M矩陣為極大似然函數(shù)J(θ)對待辨識參數(shù)向量的二階導(dǎo)數(shù)矩陣，其形式復(fù)雜，數(shù)值計算耗時，不適用于數(shù)值迭代。比較實用的迭代算法是一種改進的算法，這種算法將M矩陣做了近似處理，以縱向短周期模型為例，其最大似然函數(shù)為：

(25)

(26)

(27)

(28)

這種方法稱為改進的牛頓-拉夫遜方法。在數(shù)值迭代過程中，需要不斷的對M矩陣求逆，盡管M為對稱矩陣，但是由于數(shù)值誤差的累積作用，迭代過程中，會出現(xiàn)M矩陣的各個特征根量級差別較大，產(chǎn)生由于數(shù)值誤差而引起的奇異現(xiàn)象，導(dǎo)致迭代發(fā)散。因此在程序中還應(yīng)當(dāng)采取迭代控制來保證M矩陣求逆順利進行。采用Levenberg-Marquardt方法可以解決這種數(shù)值發(fā)散問題。該方法通過增大M矩陣的主對角元素值的大小來改善矩陣的奇異性，使矩陣求逆更加容易。具體方法是：

M-1=(M0+kA)-1

(28)

其中:M0為原信息矩陣，k為非負標量，A為一個正定矩陣，通常取為單位矩陣。k的取法有很多種。在A取為單位矩陣的情況下，建議取為M0矩陣所有元素的平均值的絕對值，即：

(29)

采用這種方法改進后，輸出誤差法在保證迭代效率的前提下，解決了迭代過程容易發(fā)散的問題，取得了較好的效果。

5 算法應(yīng)用

頻域輸出誤差法具有良好的迭代性能，但對初值的敏感性較高，使用方程誤差法和單純形法得到的初值進行迭代，輸出誤差法可計算得到精確的最優(yōu)參數(shù)解。計算流程如圖3所示。

圖3 計算流程圖

以某型電傳飛機氣壓高度11 km 、馬赫數(shù)0.9的試飛數(shù)據(jù)為例進行辨識。該型飛機的典型試飛數(shù)據(jù)時間歷程如圖4所示。按式(2)的縱向短周期模型，使用方程誤差法、單純形法和輸出誤差法進行參數(shù)辨識計算，獲取[a0，a1，b0，b1，c0，c1，τ1，τ2]等參數(shù)值，再進一步求取式(2)中傳遞函數(shù)的特征根，即可獲得該型飛機在氣壓高度11 km、馬赫數(shù)0.9的縱向短周期無阻尼自頻率ωsp、阻尼比ζsp和時間延遲τe。擬配得到的結(jié)果如圖5所示，飛機的縱向短周期響應(yīng)和低階等效系統(tǒng)的響應(yīng)重合度較高。飛行品質(zhì)指標計算結(jié)果如表1所示，飛行品質(zhì)指標與設(shè)計值相符，辨識結(jié)果合理有效。

以某型電傳飛機氣壓高度5 km 、馬赫數(shù)0.80的飛行狀態(tài)下的試飛數(shù)據(jù)為例，進行橫航向低階等效系統(tǒng)參數(shù)辨識計算，辨識結(jié)果如圖7所示。飛機的橫航向響應(yīng)和低階等效系統(tǒng)的響應(yīng)重合度較高。根據(jù)辨識結(jié)果計算飛機的橫航向荷蘭滾模態(tài)無阻尼自然頻率ωdr、阻尼比ζdr、滾轉(zhuǎn)模態(tài)時間常數(shù)τR等，計算結(jié)果如表2所示，飛行品質(zhì)指標與設(shè)計值相符，辨識結(jié)果合理有效。

圖4 某型飛機縱向短周期模態(tài)時間歷程

圖5 某型飛機縱向短周期模態(tài)擬合結(jié)果

圖6 某型飛機橫航向時間歷程

圖7 某型飛機橫航向辨識結(jié)果

飛行品質(zhì)指標設(shè)計值參數(shù)辨識值ωsp(rad/s)2.0～3.02.75ζsp0.75～0.850.80τe(ms)100～150118

表2 橫航向飛行品質(zhì)辨識結(jié)果

6 結(jié)論

本文提出了實用的低階等效系統(tǒng)模型，采用頻域相關(guān)性分析方法確定辨識頻率范圍，采用高精度有限傅立葉變換方法對試驗數(shù)據(jù)進行頻譜分析，綜合運用頻域方程誤差法、輸出誤差法和單純形法，對某型電傳飛機的低階等效系統(tǒng)模型進行參數(shù)辨識，準確地計算了該型飛機的飛行品質(zhì)指標。本文采用的方法為電傳飛機飛行品質(zhì)的評定提供了有效的技術(shù)支持[8-11]。