脈沖風洞天平短時振蕩測力數據穩態值提取的優化識別方法

王 鋒， 武 龍， 吳東升， 樂嘉陵

(中國空氣動力研究與發展中心 超高速所 高超聲速沖壓發動機技術重點實驗室， 四川 綿陽 621000)

激波風洞、炮風洞和燃燒風洞等脈沖式風洞是開展高超聲速飛行器測力試驗的重要設備。由于脈沖式風洞在極短時間內啟動，氣流以近似階躍形式作用于試驗模型，導致試驗系統結構振動，在幾十或幾百毫秒的有效試驗時間內，天平測力信號呈現大幅度振蕩，無法像連續式風洞那樣等衰減到穩定值再采集數據。因此，對于脈沖式風洞測力試驗，從振蕩信號中提取隱藏其中的穩態值是一個比較困難的問題。

對于這種短時振蕩測力信號，可以采取在試驗穩定段內取均值的方法來得到穩態測力值[1-2]。這種方法比較簡單，但理論上只有當測力信號是緩慢衰減的正弦波，可以清晰準確地獲取到整數個周期的數據，才能得到較準確的均值。然而，實際中的測力試驗系統大多是復雜的結構振動系統，特別是對于較大尺度的飛行器模型，本身低頻模態比較密集，在試驗時很多個模態被激發，測力信號形態復雜，難以判斷截取整數個周期，而且也幾乎不可能在所取時間段內，使參與振動的主要模態正好都取到整數個周期。這時，對信號取平均值作為穩態載荷測量值就有可能帶來較大誤差。當然，系統的基頻越高，穩態時段上的振動周期越多，誤差就會越小。但是，對于接近真實飛行器的大尺度試驗模型，質量很大(幾百千克，甚至超過一噸)，在保證天平靈敏度的前提下，試驗系統的基頻難以做到很高，在有效時間段內，測力信號中低頻成分的周期數可能只有3～5個甚至更少，這時平均法的誤差就可能比較大。因此，有必要采取新的方法從振動信號中提取可靠的穩態值，而不必對試驗系統的基頻提出過高要求。

1 問題的求解模型

根據脈沖式風洞的工作特點，其理想的工作過程可以描述為：在發出啟動信號后，在短暫時間內(幾毫秒至幾十毫秒)，風洞的噴管迅速建立起穩定的流場，并保持約幾百毫秒，然后流場快速減弱消失，如圖1所示。從模型的角度看，其在靜止狀態突然受到強大的瞬態氣動載荷的沖擊產生振動，然后持續受到穩定氣動載荷的作用，最后有個突然的卸載過程。可把此過程區分為加載段、保持段和卸載段。通常只關心保持段。

圖1 脈沖風洞理想工作過程 Fig.1 Ideal operation process of impulse tunnel

在加載段，試驗系統結構經歷復雜動態載荷作用下的受迫振動過程，其結果是為保持段提供了初始振動位移和速度。在保持段，試驗系統在前段提供的初始擾動下產生自由衰減振動，同時在穩態氣流作用下產生穩定的氣動載荷。測力的目的就是得到在保持段被振動信號掩蓋的穩態輸出數值，最后利用標定的天平方程計算載荷。至于卸載段則可以直接忽略。

從結構上，試驗系統一般由試驗模型、天平和支撐件組成。試驗中，系統僅允許微幅振動，因此可視為線性振動系統，可用模態疊加法來描述其振動運動，而每個模態就是一個單自由度振動系統。模型測力一般同時測量多個分量，因此天平有多路輸出，但數據處理時可逐個處理，求得每一路的穩態值，最后再一起代入天平方程求模型載荷。在此僅考慮一路輸出，基于前面的分析，根據單自由度系統的自由衰減振動響應和模態疊加原理[3]，在保持段，系統由初始擾動引起的理論輸出可寫為

(1)

式中：時間t以該段起點為零時刻，m是選取的系統振動模態數目，Ai、ξi、fi和φi分別是第i個模態的初振幅、阻尼比、頻率和初相位，A0是穩態輸出。這些參數均是未知的，但所關心的主要是A0，如果獲得了每個輸出通道的A0，利用天平方程就可計算出模型氣動載荷的各個分量。因此，數據處理的任務就是識別參數A0，當然其它參數也會同時獲得。

(2)

(3)

J中包含4m+1個未知參數，其顯式的表達式如下

(4)

下面討論如何求解該優化問題以獲得未知參數。

2 參數識別算法

式(3)是一個非線性優化問題，可以通過適當的優化算法將未知參數辨識出來。目標函數比較復雜而且參數眾多，為了能夠得到全局最小點并提高求解效率，確定較好的參數初始值十分重要。

2.1 參數初值確定方法

(5)

取Y向量的前半部分，在其幅頻曲線上取前m個最大峰值點，這些點對應的物理頻率、幅值和相位分別作為m個模態的頻率fi、幅值Ai和φi的初值

(6)

式中：fs為信號采樣頻率；k是峰值點對應的離散頻率點序號；YR、YI分別是Y(k)的實部和虛部；sgn(·)是符號函數。

關于模態數量m的確定：取Y的幅值的最大值作為參考量，將各峰值點的幅值與該量的比值作為選擇的依據，設定一個合適的閾值cr(比如0.1)，m就取比值大于cr的峰值點的個數。

這些初值的精度與信號本身以及DFT的頻率分辨率有關，而頻率分辨率由信號時長決定，為

(7)

式中：Δt是采樣周期；n為所取測量數據點數。雖然對于確定參數的大致初值來講，并不需要很精確，但是，如果試驗系統的前m階模態中存在兩個模態的頻率差小于Δf，則對應振動信號在DFT中將產生主瓣混疊，在幅頻曲線上體現為一個峰值，無法將其直接區分，這時對m的確定就有問題。例如，信號時長為100 ms時，僅為10 Hz，此種情況很有可能發生。因此，當數據的時長較短時，為了準確的確定系統重要的模態數目，需要在數據長度有限的情況下，提高DFT的頻率分辨率。對于無阻尼或弱阻尼信號，Sacchi等給出的高分辨率DFT方法[4-5]能提高短時數據DFT的頻率分辨率，可用來解決該問題。

因為實際結構阻尼比通常很小，ξi的初值就取零。

2.2 優化求解算法

將所有待辨識參數列在一起組成參數向量，令

(8)

解問題式(3)即求滿足

(9)

的參數θ。這是一個非線性方程組，在此采用擬牛頓法算法族[6]中的BFGS算法來求解。擬牛頓法避免了在迭代過程中求F(θ)關于θ的導數矩陣，也就是J(θ)關于θ的Hessian矩陣及其逆矩陣，在保持較快收斂速率的同時降低計算量。BFGS算法的迭代流程如下：

(0)準備初值θ0和一個(4m+1)×(4m+1)的正定對稱矩陣B(通常取單位矩陣)，設定收斂判斷參數εr，計算向量f0=F(θ0)；

(1)在迭代的第i步，計算fi=F(θi)，令

si=-Bifi

(10)

作為下一步的搜索方向；

(2)沿向量si執行一維搜索，求J(θi+αsi)的極小值，得到

(11)

更新參數

θi+1=θi+α*si

(12)

(3)如果‖F(θi+1)‖<εr，取θ*=θi+1，終止迭代，否則繼續下一步；

(4)更新矩陣Bi至Bi+1

(13)

其中

di=θi+1-θi，gi=F(θi+1)-F(θi)

(5)執行下一迭代步：i→i+1，轉到第(1)步。

3 應用算例

本節用4個算例來對方法的可行性進行檢驗。第一個算例是人為構造的仿真信號，用來檢驗辨識算法的準確性和魯棒性，其余3個算例均針對實際風洞試驗的測量信號，檢驗算法的實用性。

3.1 數值仿真算例

表1 信號參數真值Tab.1 True values of the parameters

本例中，已知一個由多模態振動疊加的信號，并考慮適當的噪聲，然后，由此信號辨識出該信號的直流分量和各模態的振動參數。假設信號取為

(14)

式中：ε取標準差為0.2的高斯白噪聲，其它參數任給合理的數值，如表1所示，其中第1、2模態的頻率間隔僅0.6 Hz。信號時長取200 ms，采樣頻率取2 kHz，分別用普通DFT和高分辨率DFT對信號進行分析，其幅頻曲線如圖2所示，因信號時長太短，普通DFT不能區分第1、2模態，而誤認信號只含兩個模態，高分辨率DFT可以分辨出3個模態。分別用兩種DFT結果給出的信號參數初值進行辨識計算，結果列于表 2，可見，盡管高分辨率DFT給出了3個模態的較好初值，但受噪聲影響，優化算法并不能準確收斂到三個模態的參數值。但對于我們所關心的直流項，用兩模態和三模態辨識結果精度都足夠高。用識別參數計算得到的信號與原始信號的對比如圖3所示，兩模態與三模態的重建曲線幾乎都與原始無噪聲信號重合。因此，當存在接近重頻的信號時，優化算法會收斂到一個合適的“假模態”，使得重建信號能很好擬合原始信號，而信號直流項的辨識精度并不受明顯影響。

現考察算法對噪聲的魯棒性。改變信號的噪聲水平，每個水平下隨機生成50次信號，用兩模態進行辨識，計算直流項辨識結果的相對誤差絕對值的平均值，其隨噪聲水平的變化如圖4所示，誤差隨噪聲水平近似線性增長，即使在噪聲標準差取0.5時，誤差也僅有2.5%，而直接對原始無噪信號取平均得到的值則為5.21，其相對誤差為4.2%。

圖2 普通DFT與高分辨率DFT的幅頻曲線 Fig.2 Amplitude-frequency plot of general DFT and high resolution DFT

參數i兩模態12三模態123Ai3.21751.00400.02763.23390.9852ci0.62400.13990.83770.66000.4127ωi22.59437.37717.14322.59237.346?i0.65060.1193-1.60370.64970.1175A04.99434.9974

圖3 重建信號與原始信號的對比 Fig.2 Comparison between reconstructed and original signals

圖4 直流項的辨識誤差隨信號噪聲水平的變化 Fig.4 Error of estimated DC value versus noise level of the signal

3.2 簡單模型單分量測力試驗數據

在脈沖燃燒風洞針對鈍頭錐模型進行了阻力測量試驗[7]，因為鈍頭錐模型是接近實心的模型，同時支撐系統剛度也很高，而天平僅設計了一個應變測量通道，所以整個測力系統非常接近單自由度振動系統。圖5給出了試驗測得的風洞燃燒室總壓曲線、來流皮托壓力曲線和天平輸出信號。天平信號在約0.55 s后的突變源于風洞擴壓段反射氣流的作用，是無效信號。總壓和皮托壓曲線反應了此風洞的典型運行特征。可見，在所謂穩定段(約0.1～0.55 s)，來流其實不很平穩，皮托壓一直有明顯波動。因此，在這一段，模型其實依然是受迫振動，只是由于載荷波動較小，系統振動主要由初始的沖擊引起。圖5中的天平信號可以定性地反映這點，其形態是帶有畸變的、振幅逐漸衰減的正弦曲線。

圖5 鈍頭錐吹風試驗結果 Fig.5 Test results of a blunt cone model

取0.1～0.5 s的信號作為穩態段信號，辨識得到穩態值為29.744 mV，直接取平均為29.27 mV，二者差別不太大，因為信號的周期較多。用辨識參數重建該段信號，與原始信號對比于圖6中。因為實際上該段載荷并不穩定，所以重建信號也不可能像上一例中那樣很好地吻合原信號，其偏差主要源于載荷波動引起的振動響應。

圖6 天平原始信號、重建信號與穩態值 Fig.6 The original and reconstructed balance signals and the estimated steady-state value

因為載荷不穩定，所以取不同時間段上的數據辨識得到的穩態值應該也會不同。為驗證這一點，只取時長為0.1 s的數據，從第0.1 s開始，依次后移0.05 s，直到末端到達第0.5 s，共辨識得到7個穩態值，如圖7所示，各段得到的穩態值大小是不同的，其趨勢大致與皮托壓的變化趨勢一致，說明識別方法能夠跟蹤到一定時段上模型平均載荷的變化。

圖7 天平信號分段辨識的穩態值 Fig.7 Piecewise estimated steady-state values of the balance signal

3.3 某升力體飛行器的測力試驗數據

現在以某大尺度飛行器模型在脈沖燃燒風洞的測力試驗數據來分析。僅取其升力通道的信號，如圖8中的實線所示，盡管可以看出信號包含一優勢頻率，但顯然不是單頻的簡諧振動信號。從風洞燃燒室總壓曲線觀察，大約0.8～1.1 s的時段，風洞進入穩定的試驗狀態。取0.825～1.1s時段的數據來辨識該段的穩態值，結果為3.148 mV，在圖8中以水平點劃線示出，而平均法的結果則為2.68 mV。同時辨識得到該段信號的3個主要頻率成分，分別是20.521 Hz、41.646 Hz和72.317 Hz，其幅值分別為7.483 mV、22.917 mV和2.313 mV。用辨識參數重建的信號如圖8中的點線所示，基本與原信號重合，說明參數的辨識是準確的。

由于此信號的穩態值與信號振幅相比較小，因此直接取平均的方法容易引起較大的相對誤差。

圖8 某飛行器測力信號及辨識結果 Fig.8 The original balance signal and the estimated results of a aircraft model

3.4 JF12風洞尖錐模型測力數據

JF12激波風洞[8]由于有效試驗時間稍短，約100 ms，則模型的振動周期數較少。以該風洞上10°尖錐模型三分量測力試驗數據為例，取頻率最低的俯仰通道的數據，如圖9中的實線所示。根據皮托管信號，在約0.13～0.23 s上是比較穩定的時段，取該段數據進行辨識，得到穩態值為-0.231 4 mV，以點劃線示于圖9之中，而直接平均的結果為-0.152 1 mV，由于時間短、周期數少，與辨識結果的差別較大。根據辨識參數重建的信號用點線繪于圖9中，與原信號吻合良好，間接說明辨識結果的可信性。另外，天平輸出信號中雖含有大量毛刺，但并未對參數辨識產生明顯的影響，表明方法有較好的魯棒性。

圖9 尖錐模型測力信號及辨識結果 Fig.9 The original balance signal and the estimated results of a cone model

4 結 論

本文將脈沖風洞平穩段的測力天平信號看作是試驗系統在初始條件下引起的圍繞未知穩態值的多模態振動信號的疊加，利用優化算法辨識信號主要成分的頻率、幅值、阻尼、相位等未知參數，其中零頻信號的幅值即為所要提取的天平穩態輸出值。這種方法比直接取數據均值的方法更客觀和準確，特別是對于風洞穩定段時間短、信號周期數較少的情況，優勢更加明顯。從另一方面講，采用此方法后，可適當降低對試驗系統基頻的要求，即適當降低天平測力元件的剛度，提高天平靈敏度，從而有助于提高測力精度。

但實際中，即使在平穩段，模型載荷也并非恒值，作為進一步的深化研究，可以將該段載荷也進行合理的參數化近似建模，與振動信號參數同時進行辨識，以得到更接近實際的時變載荷。

致謝

感謝中國空氣動力研究與發展中心吸氣式高超聲速技術研究中心賀偉研究員和中國科學院力學研究所劉云峰副研究員提供試驗數據。

[ 1 ] 賀偉, 童澤潤, 李宏斌. 單模塊超燃發動機推力測量天平研制[J]. 航空動力學報, 2010,25(10): 2285-2289.

HE Wei, TONG Zerun, LI Hongbin. Investigation of thrust balance for the single module scramjet[J]. Journal of Aerospace Power, 2010, 25(10): 2285-2289.

[ 2 ] 劉云峰, 汪運鵬, 苑朝凱, 等. JF12激波風洞10°尖錐標模氣動力特性試驗研究[C]∥ 第十七屆全國激波與激波管學術會議, 成都, 2016.

[ 3 ] THORBY D. Structural dynamics and vibration in practice[M]. Oxford, UK: Elsevier, 2008.

[ 4 ] SACCHI M D, ULRYCHlrych T J, WALKER C J. Interpolation and extrapolation using a high-resolution discrete fourier transform[J]. IEEE Transactions on Signal Processing, 1998, 46(1): 31-38.

[ 5 ] SACCHI M D, ULRYCH T J. Estimation of the discrete Fourier transform, a linear inversion approach[J]. Geophysics, 1996,61(4): 1128-1136.

[ 6 ] RAO S S. Engineering optimization: theory and practice[M]. 4th Ed. Hoboken, New Jersey: John Wiley & Sons, 2009.

[ 7 ] 王鋒, 任虎, 周正, 等. 載荷辨識方法用于脈沖風洞模型阻力測量研究[J]. 振動與沖擊, 2015, 34(24): 202-208.

WANG Feng, REN Hu, ZHOU Zheng, et al. Drag force measurement in impulse facilities by identification method[J]. Journal of Vibration and Shock, 2015, 34(24): 202-208.

[ 8 ] 俞鴻儒. 大幅度延長激波風洞試驗時間[J]. 中國科學：物理學 力學 天文學, 2015,45(9): 94-701.

YU Hongru. A big increase in shock tunnel test times[J]. Scientia Sinica Physica, Mechanica & Astronomica, 2015, 45(9): 94-701