王保印， 張曙光,*， 賈曉鵬

1.北京航空航天大學 交通科學與工程學院， 北京 100083 2.中國飛行試驗研究院 飛機所， 西安 710089

一種基于短時機動試飛數(shù)據(jù)的動力學辨識優(yōu)化方法

王保印1， 張曙光1,*， 賈曉鵬2

1.北京航空航天大學 交通科學與工程學院， 北京 100083 2.中國飛行試驗研究院 飛機所， 西安 710089

飛行動力學辨識算法的一個關(guān)鍵問題是，如何通過簡單的機動獲取所關(guān)心頻率范圍的響應特性。短時倍脈沖是一種易于實施的激勵信號，兼顧試飛安全性與經(jīng)濟性，但與頻域辨識法通常使用的掃頻輸入激勵相比，短時機動頻譜范圍窄、信噪比低，一般難以得到準確的辨識結(jié)果。對如何基于短時機動飛行試驗數(shù)據(jù)，提高辨識結(jié)果準確性的問題進行了研究。首先分析了經(jīng)典Welch譜估計進行時域-頻域轉(zhuǎn)換過程中，影響非參數(shù)模型辨識精度的主要因素，提出了削減窗函數(shù)邊緣縮減效應的數(shù)據(jù)預處理方法，并結(jié)合多窗口綜合技術(shù)，提高頻域特性辨識結(jié)果的精度。在參數(shù)化模型辨識過程中，針對有限頻譜范圍，提出了利用相干函數(shù)和功率譜密度加權(quán)綜合，確定等效擬配的頻率范圍和頻率節(jié)點的自適應方法，使得低階等效擬配與輸入激勵信號高度相關(guān)，提高參數(shù)化模型辨識的精度、一致性和適應性。通過不同類型飛機的大量短時機動和少量掃頻飛行試驗數(shù)據(jù)模型辨識的工程應用示例，驗證了動力學辨識優(yōu)化方法算法穩(wěn)定、結(jié)果準確，可滿足飛行品質(zhì)模態(tài)特性評價等應用需求。

辨識； 飛行試驗； 短時機動； 有限頻譜； Welch譜估計； 自適應； 低階等效

對飛行器特定輸入下的飛行試驗數(shù)據(jù)進行飛行動力學系統(tǒng)辨識，獲得應用于風洞試驗數(shù)據(jù)校準、飛行控制系統(tǒng)設(shè)計和飛行品質(zhì)評價等多個方面的支持數(shù)據(jù)，是飛行器研制工程中的一項重要工作。隨著現(xiàn)代控制理論的發(fā)展和電傳操縱系統(tǒng)的廣泛應用，使得高增益飛機飛行動力學系統(tǒng)越來越復雜，也帶來了動力學特性辨識的新問題。

時域辨識法是最早應用于飛行動力學模型系統(tǒng)辨識的方法，其主要優(yōu)點是可以直接利用時域的飛行試驗數(shù)據(jù)進行辨識，不存在中間過程誤差。20世紀60年代，時域的輸出誤差法、最大似然法、卡爾曼濾波法在動力學系統(tǒng)辨識，特別是氣動參數(shù)辨識中得到了廣泛的應用，在早期的系統(tǒng)辨識研究中，一直占據(jù)著主導地位[1-2]。飛行器頻域辨識法，其最早可追溯至1945年B-25J飛機傳遞函數(shù)辨識，隨著快速傅里葉變換(FFT)算法的發(fā)展和計算能力的提升，從20世紀70年代，頻域辨識法也有了很大的發(fā)展[3-6]。2006年，Klein和Morelli[4]系統(tǒng)地總結(jié)了頻域辨識法中的極大似然法、輸出誤差法和方程誤差法等。頻域辨識法的主要優(yōu)點是將線性微分方程或積分方程轉(zhuǎn)化為代數(shù)方程而簡化計算，而且可將被測數(shù)據(jù)選擇在感興趣的頻率范圍，具有較強的抗噪性能。雖然頻域辨識法主要應用于線性系統(tǒng)辨識，但因其獨特的優(yōu)勢，目前，依然是廣泛應用的一種飛行器系統(tǒng)辨識方法。比如，美國系統(tǒng)技術(shù)公司(STI)和NASA歷經(jīng)較長周期開發(fā)了基于FFT技術(shù)的多種飛行動力學系統(tǒng)辨識軟件[7-9]，特別是CIFER軟件，采用Chirp-z變換和組合窗處理技術(shù)，提升頻域辨識結(jié)果的精度，融入了多種數(shù)據(jù)處理方法，實現(xiàn)了多輸入多輸出的頻域辨識，得到了廣泛的工程應用[5]。

目前，由于大部分飛機依然是操縱姿態(tài)獲得軌跡改變的“經(jīng)典”飛行方式，因而其動力學描述仍可以用經(jīng)典的模型結(jié)構(gòu)，亦即仍然用“等效模型”的概念進行頻域擬配獲取系統(tǒng)模型，并進行飛行品質(zhì)評價，也就是所謂的“三步法”頻域辨識，包括數(shù)據(jù)預處理、時域-頻域轉(zhuǎn)換和頻域等效擬配等步驟[10]。其中，從辨識算法的角度，時域-頻域轉(zhuǎn)換是其核心，常轉(zhuǎn)變?yōu)楣β首V密度的計算。經(jīng)典譜估計以傅里葉變換為基礎(chǔ)，具有計算結(jié)果穩(wěn)定、效率高的優(yōu)點，其主要缺陷是方差性能差，易受頻率分辨率和窗口效應的影響。針對經(jīng)典譜估計存在的問題，現(xiàn)代譜估計以模型為基礎(chǔ)，利用采樣數(shù)據(jù)建立模型，對數(shù)據(jù)進行外推，努力改善譜估計的分辨率，對短時數(shù)據(jù)具有更強的適應性，但也存在頻率分辨率與信號的信噪比相關(guān)，譜估計質(zhì)量受模型階次影響和譜線分裂現(xiàn)象等不足[11]。

由于辨識算法問題本質(zhì)是求“逆”問題，其求解過程和結(jié)果依賴于飛行試驗數(shù)據(jù)的信噪比，受飛行試驗輸入激勵類型、測量設(shè)備精度、系統(tǒng)干擾、飛行大氣環(huán)境等因素的影響顯著，特別是基于短時機動動作的試飛數(shù)據(jù)，增加了準確辨識頻域特性的難度。針對此問題，最普遍的解決方法是，根據(jù)飛行試驗目的，通過優(yōu)化機動動作，來降低系統(tǒng)辨識的難度[12-14]。比如，美國空軍和NASA特別建議，頻域辨識法的適宜機動為數(shù)分鐘的掃頻輸入信號[3,5,7]，德國DLR則推薦使用“3211”機動[5,10]，以得到較準確的辨識結(jié)果。但是，掃頻、“3211”機動等優(yōu)化機動雖然降低了系統(tǒng)辨識的難度，卻與試飛工程的易實現(xiàn)性、安全性和經(jīng)濟性相矛盾，特別是掃頻輸入可能超出飛行員操縱預期，增大試飛風險。對于大型飛機來說，尤其是為獲得橫航向重要的荷蘭滾模態(tài)，進行腳蹬掃頻輸入激勵往往是難以實現(xiàn)的[15-17]。而降低試飛周期和耗費，也是試飛工程高度關(guān)注的重要問題。根據(jù)美國各種新機試飛統(tǒng)計，通常僅原型機本身試飛耗費，就占整個型號研制經(jīng)費的10%～20%[17]。為此，在試飛工程中，期望采用脈沖、倍脈沖等較為簡單的短時機動作為激勵輸入信號[16-17]。但是，短時激勵信號頻譜范圍窄、信噪比低，即使CIFER等成熟的頻域辨識軟件對短時機動動作，往往也不能得到理想的動力學辨識結(jié)果[5,7]。

對于數(shù)秒時長的短時機動的系統(tǒng)辨識，2007年STI基于短時脈沖輸入試飛數(shù)據(jù)(記錄數(shù)據(jù)長度不足5 s)，研究了窄帶寬特征(NBS)的系辨識法，采用小波變換理論和窄帶巴特沃斯濾波技術(shù)辨識頻域響應特性，獲得了比傳統(tǒng)FFT辨識法更理想的結(jié)果[7,18-20]。與加窗的傅里葉變換相比，小波變換具有更好的時間和頻率分辨率，但也存在著理論復雜、小波基不唯一帶來的小波基選取困難[21]。

在有人飛機飛行品質(zhì)評價的工程應用中，根據(jù)駕駛員的頻率感知和操縱特性，主要在0.1～10.0 rad/s頻率范圍內(nèi)，根據(jù)有限頻率點的頻率響應特性，進行基于傳遞函數(shù)描述的響應模型的參數(shù)辨識并評價飛行品質(zhì)[22-23]。

所以，綜合上述試飛工程的安全性、經(jīng)濟性需求和有人駕駛飛機飛行品質(zhì)評價的工程應用特點，飛行動力學辨識算法的關(guān)鍵問題，就是如何設(shè)計在較短的時間內(nèi)盡量簡單的機動，更好地獲取所關(guān)心頻率范圍的響應特性問題。為此，本文針對短時機動飛行試驗，以倍脈沖機動為重點，探討基于譜估計進行時域-頻域的高效、穩(wěn)定轉(zhuǎn)換，并進行基于最小二乘法的等效擬配頻域系統(tǒng)辨識優(yōu)化方法的研究。

1 短時倍脈沖機動的數(shù)據(jù)特征

如前所述，有人駕駛飛機飛行品質(zhì)評價所關(guān)注的頻率范圍是基于駕駛員的感知靈敏度，一般在0.1～10.0 rad/s之間，具體與飛行模態(tài)的分布有關(guān)。表1給出了4種典型飛機的快變模態(tài)頻率，可見，隨著飛機尺寸量級的逐漸增加，其快變模態(tài)的頻率呈逐漸下降趨勢，對于有人駕駛飛機最大約為5.0 rad/s量級(輕小飛機)。因而，飛行品質(zhì)規(guī)范設(shè)置的頻率上限一般為10.0 rad/s[22-25]。而飛行品質(zhì)評價的頻率下限，主要考慮執(zhí)行任務(wù)時飛行員的駕駛負擔，一般設(shè)置為0.1 rad/s。

表1 典型飛機快變模態(tài)參數(shù)Table 1 Fast variable mode parameters of typical aircraft

從飛行品質(zhì)評價的方法本身，盡管感興趣的頻帶為0.1～10.0 rad/s，但是對于低頻的模態(tài)(沉浮、螺旋等)，主要基于時域方法，通過考察其發(fā)散程度予以評估。應用等效擬配進行頻域評估，重點針對高頻模態(tài)(包括各類增穩(wěn)形成的若干高頻模態(tài)的綜合)。飛行特性及其評價的物理本質(zhì)為基于短時機動的頻域辨識和品質(zhì)評價方法建立了基礎(chǔ)。

對于脈沖、倍脈沖、“211”等短時機動，倍脈沖的對稱性最好，有利于消除不對稱的穩(wěn)態(tài)誤差。當倍脈沖輸入的主要頻譜范圍與預先估計的系統(tǒng)模態(tài)頻率基本重合時，輸入激勵對系統(tǒng)模態(tài)頻率范圍有最強的激勵作用。若預估的系統(tǒng)快變模態(tài)頻率為ω，近似將倍脈沖看作是一個周期的正弦信號，則期望輸入的倍脈沖時間周期為

T=2π/ω

(1)

通過適當調(diào)整倍脈沖寬度，調(diào)整所希望激發(fā)的頻率范圍。

倍脈沖激勵信號的幅值，則要綜合考慮響應不可超出結(jié)構(gòu)模態(tài)邊界或激發(fā)非線性特性，又不能太弱以致信噪比太低，或未激發(fā)出真實的模態(tài)響應特性，直接影響辨識精度。通常，對于迎角和側(cè)滑角的響應峰峰值以3°～5° 為宜[14,18-19]。

典型的飛行試驗輸入輸出時間歷程如圖 1所示，t為時間，δlon為縱向操縱輸入，Nz為縱向過載，q為俯仰角速率。通常一個倍脈沖激勵，其有效激勵時間為1～5 s，由式(1)對應的頻率范圍約為1.2～6.3 rad/s，涵蓋了主要有人駕駛飛機的快變模態(tài)范圍(見表1)。在一個激勵完成后，若飛機基本處于平衡狀態(tài)，則可以繼續(xù)進行下一個輸入激勵的飛行試驗。

圖1 飛行試驗數(shù)據(jù)Fig.1 Flight test data

需要考慮的是，在倍脈沖激勵的有效時間之外(如激勵響應末端)，由于受到測量設(shè)備、大氣擾動、系統(tǒng)死區(qū)和間隙等因素影響，短時機動測量的試驗數(shù)據(jù)高頻噪聲可能很大，或者可能引入其他軸向的輸入響應[17,19]。因而需要進行機動動作篩選，去掉明顯的不協(xié)調(diào)數(shù)據(jù)。

2 短時機動時域-頻域轉(zhuǎn)換方法的改進

2.1 時域-頻域轉(zhuǎn)換基本理論

信號時域-頻域轉(zhuǎn)換可通過輸入x(t)、輸出y(t)時域響應數(shù)據(jù)的離散傅里葉變換，用離散的輸出和輸入幅值比和相位差數(shù)據(jù)點，來描述系統(tǒng)的輸入-輸出動態(tài)特性關(guān)系，在頻域內(nèi)建立非參數(shù)模型。

根據(jù)信號與系統(tǒng)理論，系統(tǒng)的頻域特性為[11]

(2)

式中：Y(ω)、X(ω)分別為輸出、輸入的傅里葉變換，ω為頻率；X*(ω)為X(ω)的共軛。式(2)將系統(tǒng)的頻域響應特性計算，轉(zhuǎn)變?yōu)檩斎胼敵龅幕スβ首V密度函數(shù)Pxy(ω)和輸入的自功率譜密度函數(shù)Pxx(ω)之比計算。

韋爾奇(Welch)平均周期圖法功率譜估計[11]，把長度為N的數(shù)據(jù)XN(n)分成L段，如圖2所示，每段的長度為M，段與段之間的重疊為M-K，第i個數(shù)據(jù)段經(jīng)加窗后為

i=0,1,…,L-1;n=0,1,…,M-1

(3)

式中：w(n)為窗函數(shù)；K為整數(shù)；L為分段數(shù)。它們滿足關(guān)系式

(L-1)K+M≤N

(4)

第i個數(shù)據(jù)段的周期圖為

(5)

由此得到平均后的功率譜為

(6)

Welch平均周期圖法功率譜估計，對于掃頻等長時間尺度的機動動作，通過合理的選擇窗函數(shù)類型、窗函數(shù)寬度及段間重疊，能獲得更好的頻率分辯率及方差性能。相反的，對于倍脈沖等短時機動來說，加窗引入了新的頻域辨識誤差，必須采用新的數(shù)據(jù)處理算法，消除Welch功率譜估計窗函數(shù)的邊緣衰減效應。

評價輸入和輸出之間因果性的相干函數(shù)r2定義為

(7)

文獻[5]指出，相干函數(shù)無跳變且r2≥0.6時，頻域特性辨識結(jié)果是可接受的。

相干函數(shù)的引入對頻域辨識很有意義。在非參數(shù)模型辨識階段，根據(jù)相干函數(shù)可自動實現(xiàn)頻域響應多窗口綜合；在參數(shù)化模型辨識階段，相干函數(shù)可以為模型等效擬配確定合適的頻率范圍，以強調(diào)高精度的頻率響應數(shù)據(jù)點。

圖2 Welch法加窗與重疊Fig.2 Windowing and overlap by Welch method

2.2 短時機動時域-頻域轉(zhuǎn)換的性能優(yōu)化

為了提高譜估計的頻域分辨率，獲得功率譜的漸進無偏估計，譜估計理論與算法仍處于不斷發(fā)展之中。Tischler和Remple提出了一種多窗頻譜優(yōu)化的方法，提高頻譜的計算精度，從而使得(偏)相干函數(shù)和頻率響應的計算更精確[5]。文獻[5]指出，小的窗口寬度，增加了窗口平均個數(shù)，有利于辨識高頻特性，但是對低頻特性辨識精度下降。

針對短時機動時域-頻域轉(zhuǎn)換的性能優(yōu)化需要，本文采用Tischler多窗口綜合方法，擴展多窗口綜合方式，提高頻域辨識的精度；為了消減窗函數(shù)對短時機動的邊緣衰減效應，探討了平移機動動作的數(shù)據(jù)預處理方式。

1) 多窗口綜合

多窗口綜合通過定義多個不同寬度的窗口，進行多次辨識，利用相干函數(shù)對頻域特性辨識結(jié)果進行自動加權(quán)綜合，避免選擇窗函數(shù)寬度對頻域特性辨識結(jié)果的影響，提高頻域特性辨識準確度。文獻[5]的復合分窗方法，先對功率譜密度加權(quán)平均，再計算頻率響應特性。本文探討不同的窗口綜合方法，即先計算多個窗口下頻域響應特性，再進行加權(quán)平均平滑，達到減小頻域辨識結(jié)果的方差的目的。

加權(quán)的依據(jù)是相干函數(shù)值，對依不同窗函數(shù)j辨識的高階幅值|Gj(ω)|、相位∠Gj(ω)，加權(quán)綜合公式為

(8)

(9)

對于加權(quán)后的頻域特性，相應的用加權(quán)相干函數(shù)進行度量，定義為

(10)

圖3 多窗口綜合后頻域響應對比Fig.3 Comparison of frequency-response with multi-windows composited

圖3示出32、16、8和4 s不同窗口寬度得到的頻域響應與多窗口綜合的相應對比示例?？梢姡訖?quán)綜合法保留了不同窗口寬度下有利的頻響特征，得到了各頻率下綜合優(yōu)化的辨識結(jié)果。

2) 窗函數(shù)邊緣衰減效應的消減

為減少頻譜能量泄漏，常采用哈明窗或海寧窗對飛行試驗數(shù)據(jù)進行截斷，獲得較為真實的頻譜[5,11]。數(shù)據(jù)長度為M的哈明窗表達式為

(11)

由于窗函數(shù)的非線性邊緣衰減，在單個窗口的始端和末端均會丟失相關(guān)數(shù)據(jù)信息。Welch功率譜估計，通過窗口重疊抵消邊緣縮減，避免引起頻域響應辨識結(jié)果的偏差，這對長時掃頻機動的飛行試驗數(shù)據(jù)來說，通常易得到理想的結(jié)果。對于激勵動作僅幾秒的倍脈沖短時機動，有效信號可能僅出現(xiàn)在首個窗口左側(cè)的非線性衰減邊緣，這時窗口重疊不能抵消窗函數(shù)的邊緣衰減效應，參見式(3)，示例如圖4(a)所示。如果短時機動居于窗函數(shù)中心附近，窗函數(shù)的非線性影響可顯著減弱，如圖4(b)所示。

圖4 窗函數(shù)對倍脈沖機動飛行試驗數(shù)據(jù)的衰減Fig.4 Window function tapering to doublet maneuver flight test data

下面以短時機動頻域響應辨識為例，進一步分析窗函數(shù)邊緣衰減效應的影響。某飛機縱向升降舵δe到過載Nz的傳遞函數(shù)為

其短周期模態(tài)頻率約1.34 rad/s，由式(1)，取倍脈沖機動時長4 s。為了對比分析，取哈明窗寬度為32 s。對比機動出現(xiàn)在首個窗口左側(cè)非線性衰減區(qū)域(約6～10 s)和中心區(qū)域(約14～18 s)兩組情況，頻域辨識結(jié)果與傳遞函數(shù)計算的頻域響應分別如圖5(a)和圖5(b)所示。

從辨識結(jié)果可知，對于倍脈沖激勵，當激勵靠近窗函數(shù)的中心，不存在嚴重的非線性邊緣衰減效應時，頻域辨識結(jié)果較準確(如圖5(b)所示)，否則，將引起較大的結(jié)果偏差(如圖5(a)所示)。對弱阻尼飛機狀態(tài)，邊緣衰減效應影響更大(有較大的諧振峰值誤差)。所以，針對倍脈沖等短時機動，窗函數(shù)的邊緣衰減效應必須避免。

圖5 邊緣衰減引起的頻域辨識結(jié)果誤差對比Fig.5 Comparison of frequency domain identification results error caused by windowing edge tapering effect

消減窗函數(shù)的邊緣衰減效應核心是控制機動出現(xiàn)在窗函數(shù)中的位置。對于飛機品質(zhì)評價頻域辨識，為獲得ωmin=2π/t≈0.1 rad/s的頻率分辨率，可將時域數(shù)據(jù)長度圓整化為t≈64 s?？紤]到快變模態(tài)不會低于0.2 rad/s，窗函數(shù)的最大窗口寬度不會超過32 s。因而辨識實施時，可以通過在短時機動前后補零、平移機動動作，或復制機動的方法，確保機動動作出現(xiàn)在窗函數(shù)中心。

3 基于短時機動的頻率自適應等效擬配

在電傳操縱飛機的飛行品質(zhì)評價中，基于頻域辨識法進行飛行特性的低階等效擬配是重要環(huán)節(jié)。亦即根據(jù)試飛機動，通過確立適當?shù)牡刃M配模型結(jié)構(gòu)，采用優(yōu)化算法進行頻域特性擬合，估計模型中的未知參數(shù)[3,5,22]，最終根據(jù)辨識的參數(shù)進行飛行品質(zhì)的定量評價。

3.1 擬配頻譜范圍與擬配節(jié)點的自適應加權(quán)

對于現(xiàn)代特性迥異的電傳操縱飛機，依據(jù)已有的飛行品質(zhì)規(guī)范，如MIL-STD-1797A[23]進行等效擬配，經(jīng)常存在擬配超出包線的情況。究其原因進行分析，有兩種主要情況：① 試驗激勵機動沒有充分激發(fā)0.1～10.0 rad/s的頻率范圍；② 被評飛機已經(jīng)不能用某種經(jīng)典的低階等效模型進行描述。對于第①種情況，最理想的方法是進行掃頻等試飛機動，但是這與高效、安全試飛相違背。而對于第②種情況，屬于等效評價方法的本質(zhì)問題，必須尋找新的品質(zhì)評價方法[26-27]。如第1節(jié)所分析，等效擬配評價主要針對快變模態(tài)，所以等效該頻段特性是關(guān)鍵。由此，本節(jié)尋找基于短時機動頻譜特征，適應性的確定等效擬配的頻率，獲得精確的快變模態(tài)參數(shù)，實現(xiàn)飛機快變模態(tài)品質(zhì)特性的準確、高效評價的目的。

在飛行試驗過程中，短時機動主要激發(fā)感興趣的飛機動態(tài)頻率范圍，且在這段有限的頻率范圍內(nèi)，系統(tǒng)動態(tài)特性是可準確辨識的。等效擬配時，應增加該段頻率范圍的權(quán)重，而減弱或忽略其他未充分激發(fā)頻段或辨識結(jié)果不可靠頻段的權(quán)重。而頻域特性自適應加權(quán)的依據(jù)是相干函數(shù)值(主要確定頻率特性辨識有效頻率范圍)和功率譜密度(確定基于短時機動激發(fā)的擬配節(jié)點的位置)。

為此，建立等效擬配過程中自適應的選取擬配頻率范圍和頻率點算法。

1) 設(shè)定等效擬配的頻率點個數(shù)Nω和擬配頻率范圍[ωmin,ωmax]。

(12)

(13)

(14)

(15)

6) 取ωP序列權(quán)重為WωP(WωP∈[0,1])，一般可取等值權(quán)重WωP=0.5)，ωr序列的權(quán)重為1-WωP，得到綜合加權(quán)的等效擬配頻率節(jié)點序列為

(16)

3.2 等效擬配結(jié)果的驗證

對某一飛機而言，在同一構(gòu)型、飛行狀態(tài)附近，工程上，往往進行多次輸入激勵，通過增加試驗樣本，獲得系統(tǒng)辨識結(jié)果的漸近無偏估計，理論上，其等效擬配的參數(shù)化模型結(jié)果應是一致的。這種一致性，表明了系統(tǒng)辨識結(jié)果具有試驗上的統(tǒng)計一致性，是系統(tǒng)辨識結(jié)果應達到的目標，也是評價辨識結(jié)果性能的一個標準[3,5,18,23]。

對于辨識模型失配的顯著性水平評估，引入曲線擬合重合度指標[28]：

(17)

其中：2-范數(shù)為

(18)

3.3 示例和分析

以某大型運輸類飛機(Ⅲ類)巡航飛行(B階段)短時機動縱向飛行試驗數(shù)據(jù)為例，考核前述優(yōu)化算法。這里，優(yōu)化算法進行了平移機動動作至窗函數(shù)中心的數(shù)據(jù)預處理、擬配頻譜范圍與擬配節(jié)點的自適應加權(quán)。

采用雙擬配模型結(jié)構(gòu)[23]，高階和低階等效的頻域、時域響應和擬配節(jié)點、失配度對比結(jié)果如圖6 所示，高度H=9 784 m，馬赫數(shù)Ma=0.32，迎角α=3.2°。圖6(a)、圖6(b)分別為縱向操縱輸入δlon到過載Nz和角速率q的常規(guī)方法高階、低階與優(yōu)化方法高階、低階的頻率響應對比。圖6(c)、圖6(d)分別為δlon到Nz和q的常規(guī)方法與優(yōu)化方法的等效擬配節(jié)點和失配度對比。圖6(e)為δlon輸入下，Nz和q的試驗數(shù)據(jù)、常規(guī)方法與優(yōu)化方法的時域響應對比。

由圖6(a)～圖6(d)可見，在關(guān)心的短周期模態(tài)頻率附近，優(yōu)化算法高低階頻域響應重合性明顯優(yōu)于常規(guī)方法，可獲得更準確的等效擬配模態(tài)頻率、阻尼和系統(tǒng)等效延遲時間。圖6(e)中的時域驗證對比，也表明優(yōu)化方法和飛行數(shù)據(jù)吻合，具有較強的真實飛行動力學特性等效能力。

優(yōu)化算法的成功，不僅源于針對短時機動消減窗函數(shù)邊緣衰減數(shù)據(jù)預處理方法，可獲得更準確的頻域響應特性，也與擬配頻率范圍與擬配節(jié)點的自適應加權(quán)算法，強調(diào)了譜密度能量和相干函數(shù)值較高的頻率點相關(guān)。如圖 7所示，分別給出了δlon和Nz、q的歸一化功率譜密度曲線PSDδlon、PSDNz、PSDq，取均值后的功率譜密度(PSD composited)，和依自適應算法選擇的等效擬配節(jié)點序列(Selected Matching Points)。由圖6(c)、圖6(d)和圖7可知，對于短時機動，擬配節(jié)點主要集中在譜密度和相干函數(shù)值較高的頻率。

概率密度的加權(quán)，使擬配頻率與試驗輸入激勵信號高度相關(guān)，當激勵為掃頻時，擬配頻率自適應的擴展到接近0.1～10.0 rad/s的整個頻率范圍，并近似均勻分布。既能兼顧掃頻寬頻長時試驗激勵的整體性，又能適應倍脈沖等窄頻短時試驗激勵的局部性。引入相干函數(shù)的加權(quán)，又考慮了表征輸入、輸出的線性相關(guān)性和頻域特性辨識結(jié)果的可信度。這種擬配頻率的自適應，在機動動作的頻率特性與低階等效擬配算法之間建立了緊密聯(lián)系，增強了頻域辨識法對不同輸入激勵信號的自適應性。

圖6 參數(shù)化動力學模型頻域辨識結(jié)果Fig.6 Parameterized dynamic model of frenquency domain identification results

采用優(yōu)化算法，對算例飛機5個不同的倍脈沖激勵低階等效擬配參數(shù)，其均值(Mean)、標準差(Std)如表2所示。

由表2可知，等效擬配失配度較小，表明在關(guān)心的頻率范圍內(nèi)，高低階頻域響應基本重合。不同機動的同一個擬配參數(shù)差異很小(如Std所示)，等效擬配結(jié)果具有較好的一致性。

圖7 功率譜密度曲線與擬配節(jié)點
Fig.7 Power spectral density curves and selected matching points

表2 縱向短周期等效擬配結(jié)果Table 2 Short-period parameters of LOES

ManeuverMismatch/dBDampingradioNaturalfrequency/(rad·s-1)Nume?ratorzero/sTimedelay/s149．90．4741．7120．7540．088228．90．4311．7240．7260．089319．90．5061．6360．6320．103417．10．6141．5780．5770．102563．10．4911．6230．5800．108Mean0．5031．6550．6540．098Std0．0550．0510．0670．007

圖8 動力學辨識結(jié)果的時域一致性驗證Fig.8 Time domain consistency verification of dynamic identification result

辨識結(jié)果的部分時域一致性驗證，如圖8所示，從高低階時域響應對比來看，除了低階等效模型不能復現(xiàn)飛行試驗數(shù)據(jù)中存在高階的擾動，在短周期模態(tài)頻率附近，高低階時域響應曲線重合很好。機動1～5與機動1高低階對比的重合度指標，如表 3所示。表 3中，不同的機動動作的重合度指標，說明了參數(shù)模型真實的反映了被辨識的系統(tǒng)特性，和辨識結(jié)果的可靠性、一致性。通過其他不同類型飛機的大量短時機動和少量掃頻飛行試驗數(shù)據(jù)飛行動力學模型辨識的工程應用，也驗證了擬配頻率自適應的動力學頻域辨識法，辨識結(jié)果準確、算法穩(wěn)定可靠，可滿足飛行品質(zhì)模態(tài)特性評價等工程應用需求。

表3動力學模型辨識對不同機動的一致性

Table3Comparisonofidentificationresult’sconsistencyofdynamicmodelfordifferentmaneuvers

ManeuverFlightstateH/kmMaα/(°)F(Nz)/%F(q)/%19．70．333．287．288．82vs19．80．323．287．789．43vs19．80．323．180．089．04vs19．70．323．678．588．15vs19．70．323．581．190．7

4 結(jié) 論

針對有人駕駛飛機飛行試驗高效、安全獲取飛行動力學特性和飛行品質(zhì)評價的問題，進行試飛參數(shù)辨識方法改進研究，得到：

1) 針對短時倍脈沖機動輸入作為激勵的有利試飛實施方式，以加權(quán)相干函數(shù)作為主要衡量指標，提出了窗函數(shù)邊緣衰減效應消減的數(shù)據(jù)預處理方法，并推廣Tischler多窗口綜合方法，取得了頻域特性辨識結(jié)果精度的明顯改進。

2) 從避免試飛機動的分散性可能導致飛行品質(zhì)辨識結(jié)果的分散性問題著手，為了提高評價結(jié)果的一致性及重合度，提出了相干函數(shù)與功率譜密度綜合加權(quán)的頻率自適應辨識方法，使得低階等效擬配與輸入激勵信號高度相關(guān)，能同時兼顧寬頻長時機動和窄頻短時機動頻域辨識法參數(shù)化模型辨識的自適應性。

上述算法改進，構(gòu)成了基于短時機動的飛行品質(zhì)評價的重要基礎(chǔ)，試驗證明能夠適應電傳操縱等具有更多、更豐富模態(tài)飛機的評價需求。

[1] JATEGAONKAR R V. Flight vehicle system identification: A time domain methodology[M]. Reston: AIAA, 2006: 9-14.

[2] HAMEL P G, JATEGAONKAR R V. Evolution of flight vehicle system identification[J]. Journal of Aircraft, 1996, 33(1): 9-28.

[3] MORELLI E A. Identification of low order equivalent system models from flight test data: NASA/TM-2000-210117[R].Washington, D.C.: NASA, 2000.

[4] KLEIN V, MORELLI E A. Aircraft system identification: Theory and practice[M]. Reston: AIAA, 2006: 1-4.

[5] TISCHLER M B, REMPLE R K. Aircraft and rotorcraft system identification: Engineering methods with flight test examples[M]. Reston: AIAA, 2006:1-35, 135-190.

[6] 蔡金獅, 汪清, 王文正, 等. 飛行器系統(tǒng)辨識學[M]. 北京: 國防出版社, 2003: 221-245.

CAI J S, WANG Q, WANG W Z, et al. Aircraft system identification[M]. Beijing: National Defense Industry Press, 2003:221-245 (in Chinese).

[7] KLYDE D H, BACHELDER E N, THOMPSON P M, et al. Flying qualities parameter identification using short duration flight test inputs[C]//AIAA Atmospheric Flight Mechanics Conference and Exhibit. Reston: AIAA, 2007： 1-13.

[8] WANG K C, ILIFF K W. Retrospective and recent examples of aircraft parameter identification at NASA Dryden Flight Research Center[J]. Journal of Aircraft, 2004, 41(4): 752-764.

[9] MORELLI E A, KLEIN V. Application of system identification to aircraft at NASA Langley Research Center[J]. Journal of Aircraft, 2005, 42(1): 12-25.

[10] JATEGAONKAR R V, FISCHENBERG D, GRUENHAGEN W. Aerodynamic modeling and system identification from flight data-recent applications at DLR[J]. Journal of Aircraft, 2004, 41(4): 681-691.

[11] 胡廣書. 數(shù)字信號處理理論、 算法與實現(xiàn)[M]. 北京: 清華大學出版社, 2012: 507-590.

HU G S. Theory, algorithm and implementation of digital signal processing[M]. Beijing: Tsinghua University Press, 2012: 507-590 (in Chinese).

[12] HAMEL P G, JATEGAONKAR R V. The role of system identification for flight vehicle applications-revisited[C]//Proceedings of RTO-SCI Symposium on “System Identification for Integrated Aircraft Development and Flight Testing”, RTO-MP. Neuilly-Sur-Seine： RTO/NATO， 1999: 2.1-2.12.

[13] NETO N S B, HEMERLY E M. Parameter estimation flight test maneuver optimization considering measurement colored residuals[C]//25th International Congress of the Aeronautical Sciences. Hamburg： ICAS， 2006： 1-10.

[14] MORELLI E A. Flight test validation of optimal input design and comparison to conventional inputs[C]//AIAA Atmospheric Flight Mechanics Conference. Reston: AIAA, 1997: 11-13.

[15] THEODORE C R, IVLER C M, TISCHLER M B, et al. System identification of large flexible transport aircraft[C]//AIAA Atmospheric Flight Mechanics Conference and Exhibit. Reston: AIAA, 2008: 6894.

[16] 李雅靜, 劉旭華. 頻域響應辨識優(yōu)化方法研究[J]. 飛行力學, 2014, 32(4): 342-346.

LI Y J， LIU X H. Optimization technique of frequency domain response identification[J]. Flight Dynamics, 2014, 32(4): 342-346 (in Chinese).

[17] 周自全. 飛行試驗工程[M]. 北京: 航空工業(yè)出版社， 2010: 4-10.

ZHOU Z Q. Flight test engineering[M]. Beijing: Aviation Industry Press, 2010: 4-10 (in Chinese).

[18] BACHELDER E, THOMPSON P, KLYDE D, et al. A new system identification method using short duration flight test inputs[C]//AIAA Atmospheric Flight Mechanics Conference. Reston: AIAA, 2011: 4059-4081.

[19] KLYDE D H, BACHELDER E N, MCMAHON R, et al. Use of in-flight simulation to create a flying qualities database: AIAA-2008-1624[R]. Reston: AIAA, 2008.

[20] BACHELDER E N, KLYDE D H, LEE D, et al. Dynamic detection of actuator frequency response[C]//43rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, Joint Propulsion Conferences. Reston: AIAA, 2007: 1-13.

[21] BOGGESS A, NARCOWICH F J. A first course in wavelets with Fourier analysis[M]. Hoboken： John Wiley & Sons, 2015: 234-268.

[22] HODGKINSON J. A history of low order equivalent systems for aircraft handling qualities analysis and design[C]//2003 AIAA Atmospheric Flight Mechanics Conference and Exhibit. Reston: AIAA, 2003: 11-14.

[23] STANDARD M. Flying qualities of piloted aircraft: MIL-STD-1797A[R]. Langley： Department of Defense, 2004: 183-280.

[24] DOROBANTU A, MURCH A, METTLER B, et al. System identification for small, low-cost, fixed-wing unmanned aircraft[J]. Journal of Aircraft, 2013, 50(4): 1117-1130.

[25] CARNDUFF S D. System identification of unmanned aerial vehicles[D]. Canfield： Canfield University, 2008.

[26] 譚文倩, 張曙光. 飛機橫航向等效擬配模型的對比研究[J].飛行力學, 2003, 21(4): 14-18.

TAN W Q, ZHANG S G. Research on the models of lateral and directional equivalent matching for an aircraft[J]. Flight Dynamics, 2003, 21(4): 14-18 (in Chinese).

[27] 張曙光, 方振平, 高立華. 縱向大迎角飛行品質(zhì)對參數(shù)的靈敏度分析[J]. 飛行力學， 1999， 17(2): 1-6.

ZHANG S G, FANG Z P, GAO L H. Analysis on sensitivity of high angles of attack longitudinal flying qualities to variation of parameters[J]. Flight Dynamics, 1999, 17(2): 1-6 (in Chinese).

[28] LENNART L. MATLAB R2012b system identification toolbox reference[M/OL]. Natick： MathWorks， 2012： 131-139 (2012-09)[2016-08-15]. www.mathworks.com.

(責任編輯： 鮑亞平)

An optimization method for dynamic identification based onshort-duration maneuvering flight test data

WANGBaoyin1,ZHANGShuguang1,*,JIAXiaopeng2

1.SchoolofTransportationsScienceandEngineering,BeihangUniversity,Beijing100083,China2.InstituteofAircraft,ChineseFlightTestEstablishment,Xi’an710089,China

How to generate the accurate response at the frequency range of interest from the maneuvering is one of the significant problems for dynamic identification. Although short-duration doublet is commonly used in flight test activities and is easy to implement with a good consideration of both security and economy, it is difficult to get the accurate estimation due to its limited spectrum and lower signal to noise ratio. Thus, a method focusing on improving estimation accuracy based on the short-duration maneuvering flight test data is developed in this paper. The main factors that affect the accuracy of the non-parametric model identification in the classical Welch method of spectrum estimate are analyzed. A novel data pre-processing method which can taper the window function edge is presented, and the multi-window composite technology is integrated to improve the identification accuracy. Regarding the limited frequency spectrum, an adaptive low-order equivalent matching method is developed to select the frequency range and node based on a weighted function of coherence function and power spectral density. This method make the low-order equivalent matching highly correlated with the maneuvering input signal, and improve the accuracy, consistency and adaptability of the model parameter identification. The dynamic identification optimization technique is applied to a large amount of short-duration maneuvering and several sweeping flight test data for various types of aircrafts. The results meet the accuracy requirement for flight quality evaluation application, while the algorithm is stable and reliable.

identification; flight test; short-duration maneuver; limited frequency spectrum; Welch method of spectrum estimate; adaptive; lower-order equivalent

2016-09-26；Revised2016-10-25；Accepted2016-11-09；Publishedonline2016-11-211439

2016-09-26；退修日期2016-10-25；錄用日期2016-11-09； < class="emphasis_bold">網(wǎng)絡(luò)出版時間

時間：2016-11-211439

www.cnki.net/kcms/detail/11.1929.V.20161121.1439.014.html

*

.E-mailgnahz@buaa.edu.cn

王保印， 張曙光， 賈曉鵬． 一種基于短時機動試飛數(shù)據(jù)的動力學辨識優(yōu)化方法J． 航空學報,2017,38(6):120815.WANGBY,ZHANGSG,JIAXP.Anoptimizationmethodfordynamicidentificationbasedonshort-durationmaneuveringflighttestdataJ.ActaAeronauticaetAstronauticaSinica,2017,38(6):120815.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0293

V212.1

A

1000-6893(2017)06-120815-12

URL：www.cnki.net/kcms/detail/11.1929.V.20161121.1439.014.html

*Correspondingauthor.E-mailgnahz@buaa.edu.cn