基于時序模型的顫振預測技術研究

劉立坤 張春蔚 王東森

摘要：本文主要研究了顫振邊界預測應用技術，針對湍流激勵下的結構響應數據，采用遺忘因子最小二乘法建立時序模型，結合Jury穩定性判據構造不同飛行速度下的穩定性參數，通過外推得到顫振臨界速度預測結果。以仿真試驗為例，討論了臨近邊界速度大小和響應數據時間長度等因素對預測精度的影響。最后結合顫振風洞試驗實測數據驗證了該方法的合理性和工程實用性。

關鍵詞：顫振邊界預測;時序模型;遺忘因子最小二乘法

中圖分類號：V217 文獻標識碼：A

顫振飛行試驗是新型或結構有重大改變的飛行器都必須進行的試飛項目之一，它使用真實飛行器在實際飛行條件下進行試驗，具有典型的高風險、高耗費和長周期的特點。顫振邊界預測技術是顫振試飛中確定顫振邊界的主要手段。由于顫振飛行試驗條件的復雜性，試驗數據的品質較低，顫振數學模型存在諸多不確定因素，顫振包線一般很難直接從試驗獲得，因此如何基于亞臨界響應測量與分析進而完成對包線的預估十分重要[1]。

國內在風洞顫振試驗中，常用的亞臨界顫振邊界預測方法主要包括阻尼外推法和穩定性參數法，但多數都需要通過模態參數辨識，獲取不同條件下結構的固有頻率和阻尼比。阻尼外推法通過直接外推臨界參數確定顫振邊界;穩定性參數法則由模態頻率和阻尼換算得到特征多項式的復數特征根，通過根和系數的關系構建勞斯判據或Jury判據穩定性參數，然后進行臨界參數外推確定顫振邊界。穩定性參數法在國內外已有廣泛的應用，常見的著作和論文主要有管德院士的《氣動彈性試驗》[2]、宇航出版社出版的《防空導彈結構與強度》[3]、中國空氣動力研究與發展中心郭洪濤等編寫的《顫振試驗數據處理技術研究》[4]等。然而，由于存在激勵力不足、測量點有限、試驗條件復雜及測試條件惡劣等問題，顫振試驗特別是飛行試驗測量的試驗數據品質較低，數據往往具有模態密集、有效樣本短和數據信噪比低等特征，這些問題的存在大大增加了模態參數識別的技術難度，進而影響到顫振邊界預測結果的準確性。

目前，國內基于時序模型的顫振邊界預測方法較為少見，本文通過直接識別時序模型多項式系數獲得穩定性參數，不必經過模態參數辨識，提供了一種基于時序模型的顫振邊界預測方法，以仿真試驗為例，討論臨近顫振邊界不同速度大小和響應數據時間長度等因素對預測精度的影響。最后結合顫振風洞試驗實測數據驗證了該方法的合理性和工程實用性。

1 方法介紹

1.1 遺忘因子最小二乘算法

遺忘因子最小二乘算法最初由Lennart Ljung于1999年提出[5]，廣泛應用于控制系統辨識領域，其思想是利用前一時刻的響應對下一時刻進行預測，根據當前的預測誤差設置適當的增益值，來影響預測的趨勢，經過迭代運算，使得最終的預測誤差達到最小，并以誤差最小時的參數作為系統辨識的結果。具體算法如下，假設時間序列{y（t）;t=1，2，…，N}是測量到的一組響應值，則有：式中：ε（t）是t時刻預測誤差Ψ班（t）是t時刻斜率，θ（t）是t時刻的參數估計值，y（t）是t時刻的響應值（t）是由t-1時刻響應得出的t時刻響應預測值，增益K（t）為當前預測誤差對參數估計的影響程度，λ是遺忘因子取值為1，通過不同時刻的循環迭代使得J（θ）最小，并以θ值作為系統辨識的結果[6～8]。

1.2 ARMA氣動彈性系統辨識[9～12]

假設{y（t）;t=1，2…，N}為N自由度氣動彈性系統經大氣湍流激勵的響應，將其表示成自回歸滑動平均（ARMA）模型的差分方程形式：式中：e（t）是零均值、方差桅σ²的白噪聲;變量z^-1為后移算子，z^-1y（t）=y（t-1），多項式系數A（z^-1）和B（z^-1）分別為：式中：n和m分別為自動回歸（AR）和滑動平均（MA）多項式階次。

由于白噪聲項e（t）無法測量，由過去時刻的數據y（t）得到的預測旬（t）為：

定義預測誤差為：

聯立式（8）和式（11）得：

將式（13）帶入遺傳因子最小二乘算法，經過迭代運算，即可得到ARMA多項式系數。

1.3 穩定裕度估計

假設經過遺忘因子算法辨識出ARMA多項式系數，將其表示成特征多項式形式：

對于離散系統使用Jury穩定性判據判斷系統的穩定性，對于k=1，3，…，n-1，定義k階方陣X_k、Y_k為：Jury判據要求：

Jury判據要求：

計算不同速度、不同模態結構的穩定性因子F^-（n-1），n為ARMA模型階次，繪制穩定性因子隨速度的變化曲線，并使用二次曲線擬合的方法進行外推，就能得到穩定性參數為零時對應的速度值，即為顫振臨界速度值。

2 仿真數據驗證分析

為驗證該方法在湍流激勵情況下的有效性，以得克薩斯A&M大學航空航天工程系的俯仰一滾轉氣動彈性風洞模型為算例[3]，主要研究了不同速度區間下的數據點數和響應數據時間長度對預測結果的影響。

在MATLAB中進行數值仿真試驗，其動力學模型見式（18），具體取值見參考文獻[6]。該模型為二自由度系統，模型階次真值為4，顫振速度V_flutter為12.7m/s，對應的速度頻率和速度阻尼曲線如圖1和圖2所示。

采用零均值高斯白噪聲模擬真實大氣湍流激勵形式進行激勵，首先分析了響應信號時長為60s、無噪聲情況下不同速度區間數據的預測結果，見表1。對于新型或重大改型飛機，速度在0.8V_d～0.98V_d之間須經過顫振飛行試驗驗證，（V_d為最大設計俯沖速度，根據CCAR25部相關條款V_d≤V_flutter/1.15），表1給出了8種典型速度區間的預測結果。其中，fz-fit為時序模型預測結果，g-fit為阻尼外推法預測結果。圖3和圖4給出了兩種情況的顫振速度預測結果圖[7-9]。

分析結果表明，低速試驗點時，阻尼較為分散，阻尼外推法預測誤差較大，時序模型預測結果比阻尼外推法精度高且較為保守。隨著速度逐漸接近臨界顫振速度，兩種方法預測結果趨于一致并且都較為準確。

在顫振飛行試驗中，為了得到有足夠置信度的信息，大氣湍流激勵方法常要求飛行員持續保持飛機穩定平飛狀態一定時間，采集足夠長的結構響應信號，這大大增加了試驗的成本和風險。選取速度區間0.8V_d～0.98V_d，針對不同時間長度的響應信號數據進行了分析，表2給出了具體分析結果。

分析結果表明，隨著響應信號數據樣本的減少，基于時序模型的預測方法預測精度逐漸降低。但即使是5s的短時數據樣本情況，該方法仍然具有一定的精度。

3 顫振飛行試驗數據處理中的應用

某型飛機風洞試驗，采用大氣湍流激勵，利用布置在平尾、垂尾等位置的加速度傳感器測量結構響應。在速度201～381m/s的部分試驗點進行試驗，圖5～圖7給出了部分速度下左平尾結構響應自功率譜曲線，在速度381m/s時，頻域兩模態頻率歸一、響應幅值明顯放大，速度381m/s已非常接近顫振速度。考慮到試驗的安全性，在該速度點終止試驗。分別使用基于時序模型的預測方法和阻尼外推法進行顫振速度預測[10～12]。

圖8給出了部分試驗點的顫振速度預測結果，fz為穩定性因子值，g為結構阻尼比，fz-fit為穩定性因子外推擬合曲線，g-fit為阻尼外推擬合曲線。圖中，阻尼外推法顫振速度預測值約為411m/s，基于時序模型的顫振預測方法預測值約為403m/s，兩種方法預測結果一致性較好。

4 結論

采用基于時序模型的顫振邊界預測方法，對湍流激勵下的氣動彈性模型進行顫振邊界預測，并與傳統阻尼外推法進行了對比分析。仿真試驗結果表明，基于時序模型的顫振預測方法在接近顫振臨界速度時與阻尼外推法預測結果一致，并且由于穩定性因子隨速度成近似線性變化的關系，在低速試驗點時其預測精度較阻尼外推法有一定的優勢。

將該方法應用于某型飛機風洞試驗數據，對平尾進行顫振速度預測，結果表明，該方法能夠得出工程適用的顫振速度預測值。此外，與阻尼外推法相比，該方法對于湍流激勵得到的短時響應信號也有一定的適應能力，在實際飛行試驗中，能夠有效保障科研試飛的安全，降低飛行試驗的風險，具有一定的工程使用價值。

參考文獻

[1]張偉偉，鐘華壽，肖華，等.顫振飛行試驗的邊界預測方法回顧與展望[J].航空學報，2015，36（5）：1367-1384.

[2]管德.氣動彈性試驗[M].北京：北京航空航天大學出版社，1986.

[3]許椿蔭.防空導彈結構與強度[M].北京：中國宇航出版社，2006.

[4]郭洪濤.顫振試驗數據處理技術研究[D].綿陽：中國空氣動力研究與發展中心，2005.

[5]Jung L.System identification：Theory for the user[M].Beijing：Tsinghua University Press，2002.

[6]Baldelli D H，Zeng J，Lind R，et al.Flutter-prediction toolfor flight-test-based aeroelastic parameter-varying models[J].Journal of Guidance Control&Dynamics，2009，32（1）：158-171.

[7]Matsuzaki Y Flutter boundary prediction based on Jury，s stabil-ity criterion：a review[C]//Aiaa/asme/asce/ahs/asc Structures，Structural Dynamics and Materials Conference，2011.

[8]Matsuzaki Y.Flutter boundary prediction of digitalized smartmultimode systems using steady-state responses[C]//Aiaa/asme/asce/ahs/asc Structures，Structural Dynamics and MaterialsConference，2006.

[9]Torii H，Matsuzaki Y Flutter margin evaluation for three-modediscrete-time systems[J].Journal ofAircraft，2013，38（1）：42-47.

[10]AlAA.Flutter prediction of multimode systems from their tur-bulent responses by Jurys criterion[C]//Aiaa/asme/asce/ahs/asc Structures，Structural Dynamics and Materials Conference，2013.

[11]Matsuzakir Y Flutter boundary prediction of an adaptive smartwing during process of adaptation[R].Proc Spie，2005：5764.

[12]Matsuzaki Y，Torii H.Flutter boundary prediction of a smartwing in process of structural adaptation[C]//Aiaa/asme/asce/ahs/asc Structures，Structural Dynamics and Materials Conference，2007.