譚　博　鄭　華

西北工業(yè)大學,西安,710072

基于顫振試驗數(shù)據(jù)分析的矩陣束方法性能研究

譚博鄭華

西北工業(yè)大學,西安,710072

為將矩陣束引入顫振試驗數(shù)據(jù)處理的工程應用領域,基于隨機信號仿真,采用蒙特卡羅方法分析了該方法的數(shù)值性能,研究了樣本長度、信噪比及計算參數(shù)對計算性能的影響,并在飛機機翼氣彈模型的風洞顫振試驗中進行了驗證。與傳統(tǒng)頻域方法的比較分析表明,矩陣束方法性能良好,是一種可靠的模態(tài)參數(shù)估計方法。

矩陣束;蒙特卡羅方法;模態(tài)參數(shù)識別;顫振試驗數(shù)據(jù)分析

0　引言

飛機結構顫振是飛機研制過程中必須進行的實驗科目,具有風險高、周期長等特點，且其觀測信號具有有效樣本短、信噪比小、結構模態(tài)密集等特征，從而使得如何從實際觀測數(shù)據(jù)中提取結構模態(tài)參數(shù)成為該行業(yè)十分關注的問題。

矩陣束方法作為一種模態(tài)參數(shù)估計方法,在許多領域有了廣泛應用。但是，在引入飛機結構顫振試驗時,該方法的參數(shù)估計質量經(jīng)常會受到樣本長度、信噪比等的影響。為此，本文基于隨機信號仿真,應用蒙特卡羅方法對其相關的數(shù)值性能進行了分析研究。

1　矩陣束方法簡介

假設觀測到的系統(tǒng)響應信號可表示為M個模態(tài)的指數(shù)函數(shù)的線性組合：

(1)

其中,x(t)、n(t)、y(t)分別為系統(tǒng)響應、系統(tǒng)噪聲和含噪聲的實測信號,0≤t≤T;T為最大觀測時間。對第i個模態(tài),算子si=-αi+jωi可用于表示模態(tài)的頻率和阻尼比系數(shù),其中ωi為角頻率,αi為阻尼比系數(shù)。

式(1)的離散時間形式為

(2)

k=0,1,2,…,N

其中,zi為系統(tǒng)響應的極點,zi=exp(siTs);Ts為采樣周期;N為最大采樣點數(shù)。

由采集序列y(t)可以構造如下的Hankel矩陣：

(3)

式中,L為矩陣束參數(shù),通常取值位于N/4～N/3之間。

對矩陣Y進行奇異值分解,得到

Y=UDVT

(4)

式中,U為N-L階正交矩陣;V為L+1階正交矩陣;D為半正定的(N-L)×(L+1)階對角陣，其主對角線上的元素為奇異值(由大至小排列)。

在已知模態(tài)個數(shù)M的情況下,由D的前M個較大的非零奇異值形成新的矩陣：

(5)

Δ=diag(δ1,δ2,…,δM)

(6)

(7)

(8)

2　數(shù)值性能研究

2.1采樣點數(shù)影響

處理試驗數(shù)據(jù)時,依據(jù)不同的試驗環(huán)境、結果精度、速度要求以及測試設備條件等設置不同的采樣長度后，方能進行模態(tài)參數(shù)估計,而采樣長度往往會影響最終數(shù)據(jù)處理結果。因此本文首先通過仿真數(shù)據(jù)來研究樣本長度對矩陣束方法數(shù)值性能的影響。設置一個單模態(tài)自由度系統(tǒng),其阻尼比系數(shù)為0.03,振動頻率為11 Hz,激勵信號為隨機白噪聲，響應信號的采樣率為128,仿真信號的采樣長度從1 s逐秒增加至60 s，針對每一個采樣長度進行3000次運算，將計算所得模態(tài)的頻率和阻尼比進行統(tǒng)計，分別將其高斯分布中值作為最終結果,如圖1、圖2所示。圖1與圖2中，橫軸表示采樣點數(shù)，范圍為128～7500,對應的采樣時間為1～60 s。

圖1　計算阻尼比系數(shù)隨采樣點數(shù)變化曲線

圖2　計算頻率隨采樣點數(shù)變化曲線

從圖1、圖2可以明顯看出,樣本點過少會極大地影響矩陣束方法計算結果的精度,導致最終計算結果與設定真值相差過大。增加采樣長度可提高計算結果的精度，使其接近真實值。當采樣點數(shù)達到一定數(shù)量之后,單純地增加采樣點不能明顯改善矩陣束方法的數(shù)值性能,反而會延長計算時間，導致系統(tǒng)的運行時間間隔增加，影響方法的計算效率。由于在飛機結構顫振試驗中常有實時分析的要求，因此，采樣點數(shù)的設置應在計算效率和結果精度之間做出平衡。

2.2信噪比對數(shù)值性能的影響

試驗采集得到的信號總會受到噪聲的影響。因此模態(tài)參數(shù)估計方法的抗噪性能也是選擇方法時需考量的重點之一。

蒙特卡羅方法是一種通過對大量彼此獨立的試驗結果進行統(tǒng)計，得出統(tǒng)計對象數(shù)學特征的方法。由于實際中的試驗周期較長,無法進行大量的重復試驗，因此在本節(jié)中，將通過仿真信號對矩陣束方法在噪聲環(huán)境下的性能進行研究。設置一個雙模態(tài)自由度系統(tǒng),模態(tài)參數(shù)分別設置為:f1=11 Hz,f2=13 Hz,d1=0.05,d2=0.04,f與d分別表示模態(tài)參數(shù)中的頻率與阻尼,下標1、2表示模態(tài)1和模態(tài)2。將響應信號的信噪比依次置為無噪聲(信噪比無窮大)、20 dB、10 dB、6 dB、3 dB和0,使用矩陣束方法對加噪后的信號進行模態(tài)參數(shù)估計運算,對所得結果進行統(tǒng)計,表1所示為3000次計算的統(tǒng)計結果?？梢钥闯?，矩陣束方法估計的模態(tài)參數(shù)中,頻率受噪聲影響較小,隨著信噪比的下降沒有發(fā)生明顯變化，而估計所得阻尼隨著信噪比的下降有較大變化。

表1　不同信噪比下矩陣束方法的計算結果

2.3設置參數(shù)對數(shù)值性能的影響

從仿真實驗結果可以看出,在僅改變信噪比的情況下，矩陣束方法所得的阻尼估計結果會產(chǎn)生較大的波動。為削弱噪聲對矩陣束方法阻尼估計結果的影響，增強矩陣束方法的抗噪性能，可以適當調整方法的設置參數(shù)，即模態(tài)個數(shù)M。

模態(tài)個數(shù)M對應于運算時設置的極點個數(shù)。增加設置的模態(tài)個數(shù)M,會導致最終計算結果中出現(xiàn)實際系統(tǒng)中不存在的虛假模態(tài)，但可有效地削弱噪聲對真實模態(tài)的計算結果的影響。

下面以1個三模態(tài)自由度且有密集模態(tài)的系統(tǒng)為例，說明改變參數(shù)對矩陣束方法的數(shù)值性能的影響。系統(tǒng)的3個模態(tài)的參數(shù)設置分別為:f1=11 Hz,f2=12 Hz,f3=17 Hz,d1=0.05,d2=0.04,d3=0.01,其中,11 Hz、12 Hz的模態(tài)為系統(tǒng)的密集模態(tài)頻率。為系統(tǒng)的響應信號添加信噪比為6 dB的噪聲后得到采集信號。

(a)M=3

(b)M=5

(c)M=6圖3　矩陣束方法估計結果分布

圖3所示為統(tǒng)計得到的計算結果分布情況,圖中，橫軸表示計算所得模態(tài)參數(shù)的頻率,范圍為9～33 Hz,縱軸為阻尼比系數(shù)，范圍為0～0.16。圖3中的每個點代表通過矩陣束計算得到的一個模態(tài)。

圖3a所示為M=3的結果分布情況,可以看出,由于受到噪聲的干擾，系統(tǒng)的2個密集模態(tài)較難分辨,但是可以清晰地分辨出設計系統(tǒng)位于17 Hz的模態(tài)。圖3b所示為M=5的計算結果分布,可以看到當設置模態(tài)個數(shù)增加到5后,矩陣束方法對系統(tǒng)的2個密集模態(tài)有了較好的分辨能力,同時17 Hz的模態(tài)仍然被可以清晰地進行辨識,圖中右上方分布的計算結果點即為矩陣束方法計算結果中的虛假模態(tài)。將圖3a、圖3b進行比較后不難發(fā)現(xiàn)，圖3b所示的計算結果相較圖3a的結果更加接近設定系統(tǒng)的真實值，其分布也更加集中，這表明，增加模態(tài)個數(shù)不僅可以減小計算結果的誤差，而且可以減小結果的統(tǒng)計方差。在結構顫振試驗中，采集信號往往長度有限，處理數(shù)據(jù)時難以采用多次計算后進行統(tǒng)計的方式來減小誤差，在這樣的背景下，減小計算結果的統(tǒng)計方差可以增加結果的置信度。

圖3c所示為將模態(tài)個數(shù)再次增加后的計算結果分布。由圖3b、圖3c的比較可以看出，此時設置的模態(tài)個數(shù)為實際模態(tài)個數(shù)的2倍,盡管系統(tǒng)的2個密集模態(tài)仍可較為清晰地分辨出來，但是其結果分布明顯比圖3b所示結果更分散。這表示矩陣束方法的計算結果會在一個較大的范圍內波動，降低了計算結果的置信度，不利于確定系統(tǒng)的模態(tài)參數(shù)。因此若模態(tài)個數(shù)設置過大，反而會降低矩陣束方法的結果的估計精度，而且增加模態(tài)參數(shù)的個數(shù)會延長系統(tǒng)的運算間隔，降低計算效率，導致系統(tǒng)成本的增加，因此應避免設置模態(tài)數(shù)過大。

3　在顫振試驗數(shù)據(jù)分析中的應用

本文使用的試驗數(shù)據(jù)由某型飛機模型風洞試驗數(shù)據(jù)得來。經(jīng)零均值化處理后，選取顫振發(fā)生前的5個速度對應的5段數(shù)據(jù)作為計算依據(jù),分別應用矩陣束法與半功率帶寬法對數(shù)據(jù)進行分析，給出擬合曲線以及預測結果。本次使用的試驗數(shù)據(jù)中，將顫振發(fā)生時的風速設置為1，其余各速度依此進行歸一化處理。顫振發(fā)生時刻采集信號的時域波形及頻譜如圖4所示。

(a)時域波形

(b)頻譜圖圖4　顫振發(fā)生時的采集信號

由圖4可以看到,顫振發(fā)生時信號能量集中于8 Hz附近。依次使用矩陣束法和半功率帶寬法對顫振臨界點前的五段數(shù)據(jù)進行模態(tài)分析計算,對結果進行擬合后得到的曲線如圖5所示。

(a)矩陣束方法

(b)半功率帶寬法圖5　兩種計算方法的擬合結果

對比圖5a、圖5b可以看出，合理選擇算法參數(shù)后，矩陣束方法計算得到的阻尼結果趨勢較明顯，其最終估計的顫振速度為0.993，與真實值相差較小；半功率帶寬法的計算結果波動較大，擬合曲線不平滑，其最后估計的顫振速度為1.103，與真實值相差較大。產(chǎn)生這一結果的主要原因是，對于密集模態(tài)而言，半功率帶寬法受到自身算法的數(shù)學特性的限制，難以對模態(tài)之間的邊界進行清晰的區(qū)分,因此其對密集模態(tài)的分辨效果不如矩陣束法好，計算阻尼和真實值偏差較大。結構顫振試驗中，密集模態(tài)的分辨問題是一個經(jīng)常面對的問題,這也是半功率帶寬法在飛機結構顫振試驗中表現(xiàn)不佳，需要引入矩陣束法的原因之一。

4　結語

在保證系統(tǒng)運算效率的前提下，增加采樣點數(shù)可以提高矩陣束方法的計算精度；矩陣束方法用于模態(tài)參數(shù)估計時，其得到的頻率估計結果受到噪聲影響較小，且相對于阻尼比系數(shù)估計結果更加穩(wěn)定；在檢測信號含有噪聲的情況下，可以適當增大模態(tài)個數(shù)，用虛假模態(tài)來削弱噪聲的影響，以有效改善最終結果。

顫振試驗數(shù)據(jù)的處理結果表明，在模態(tài)阻尼比估計中，矩陣束法計算結果比較理想，其擬合曲線較為平滑，而且與傳統(tǒng)的半功率帶寬法相比，其估計結果具有更高的數(shù)值精度。

[1]Hua Yingbo,Sarkar Tapan K.Matrix Pencil Method for Estimation Parameters of Exponentially Damped/undamped Sinusoids in Noise[J].IEEE Transactions on Acoustics Speech and Signal Processing,1990,38(5):814-824.

[2]Thomas A J,Chard J,John E,et al.Defining a Bearing Replacement Strategy Using Monte Carlo Methods[J].International Journal of Quality & Reliability Management,2011,28(2):155-167.

[3]Potts D,Tasche M.Parameter Estimation for Nonincreasing Exponential Sums by Prony-like Methods[C]//17th Conference of the International Linear Algebra Society.Braunschweig,Germany,2013:1024-1039.

[4]朱瑞可,李興源,王渝紅,等.基于矩陣束算法的諧波和間諧波參數(shù)估計[J].計算機仿真,2012,40(3):388-391.

Zhu Ruike,Li Xingyuan,Wang Yuhong,et al.Parameter Identification of Harmonics and Interharmonics Based on Matrix Pencil Algorithm[J].East China Electric Power,2012,40(3):388-391.

[5]張敏，黃俐，李文雄，等.大型結構模態(tài)參數(shù)識別研究[J]，建筑科學與工程學報，2013，30(2)：49-54，75.

Zhang Min，Huang Li，Li Wenxiong，et al. Research on Modal Parameter Identification on Large-scale Structure[J]. Journal of Architecture and Civil Engineering，2013，30(2)：49-54，75.

[6]徐利,鄒傳云，陳民，等. 基于矩陣束算法的極點提取分析[J]. 通信技術，2012，45(6)：58-60.

Xu Li,Zou Chuanyun,Chen Min,et al.Analysis on Pole Extraction Based on Matrix Pencil Algorithm [J].Communications Technology，2012，45(6)：58-60.

[7]賈天嬌，岳林. 一種模態(tài)弱響應且模態(tài)秘籍的參數(shù)識別方法[J]. 中國機械工程，2012，23(11)：1313-1317.

Jia Tianjiao，Yue Lin. A Parameter Identification Method for Weak Modal Response with Close Modes[J]. China Mechanical Engineering，2012，23(11)：1313-1317.

[8]楊明華,奇異矩陣束的標準形與廣義逆[D]. 哈爾濱：哈爾濱工業(yè)大學,2013.

[9]鄭敏,申凡,史東鋒,等.單獨利用響應數(shù)據(jù)進行模態(tài)分析[J]. 中國機械工程,2006,17(4):405-409.

Zheng Min,Shen Fan,Shi Dongfeng,et al.Modal Analysis from Response-only[J].China Mechanical Engineering,2006,17(4)：405-409.

[10]黃應來,董大偉,閆兵.密集模態(tài)分離及其參數(shù)識別方法研究[J].機械強度,2009,31(1):8-13.

Huang Yinglai,Dong Dawei,Yan Bing.Study on Closely Spaced Modes Decomposition and Modal Parameter Identification[J].Journal of Mechanical Strength,2009,31(1):8-13.

(編輯張洋)

Research on Numerical Performance of Matrix Pencil Based on Flutter Test Analysis

Tan BoZheng Hua

Northwestern Polytechnical University,Xi’an,710072

In order to apply MP algorithm into flutter test data processing,the characteristics were studied through Monte Carlo method,the effects of signal length,signal-noise ratio and parameters were mainly concerned.Then MP algorithm was tested with real collected data from physical experiments.The results show that MP is a valuable method to identify modal parameters for flutter test analyzses compared with traditional ways of data processing.

matrix pencil(MP);Monte Carlo method;modal parameter identification;flutter test

2013-10-10

國家自然科學基金資助項目(11302175)

TP391.9DOI：10.3969/j.issn.1004-132X.2015.02.019

譚博,男,1987年生。西北工業(yè)大學動力與能源學院博士研究生。主要研究方向為人機環(huán)境與工程。發(fā)表論文7篇。鄭華,男,1983年生。西北工業(yè)大學動力與能源學院博士后研究人員。