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

譚　博　鄭　華

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

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

譚博鄭華

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

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

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

0　引言

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

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

1　矩陣束方法簡介

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

(1)

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

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

(2)

k=0,1,2,…,N

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

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

(3)

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

對矩陣Y進(jìn)行奇異值分解,得到

Y=UDVT

(4)

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

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

(5)

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

(6)

(7)

(8)

2　數(shù)值性能研究

2.1采樣點(diǎn)數(shù)影響

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

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

圖2　計(jì)算頻率隨采樣點(diǎn)數(shù)變化曲線

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

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

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

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

表1　不同信噪比下矩陣束方法的計(jì)算結(jié)果

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

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

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

下面以1個(gè)三模態(tài)自由度且有密集模態(tài)的系統(tǒng)為例，說明改變參數(shù)對矩陣束方法的數(shù)值性能的影響。系統(tǒng)的3個(gè)模態(tài)的參數(shù)設(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)的響應(yīng)信號添加信噪比為6 dB的噪聲后得到采集信號。

(a)M=3

(b)M=5

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

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

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

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

3　在顫振試驗(yàn)數(shù)據(jù)分析中的應(yīng)用

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

(a)時(shí)域波形

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

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

(a)矩陣束方法

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

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

4　結(jié)語

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

顫振試驗(yàn)數(shù)據(jù)的處理結(jié)果表明，在模態(tài)阻尼比估計(jì)中，矩陣束法計(jì)算結(jié)果比較理想，其擬合曲線較為平滑，而且與傳統(tǒng)的半功率帶寬法相比，其估計(jì)結(jié)果具有更高的數(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ì)[J].計(jì)算機(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]張敏，黃俐，李文雄，等.大型結(jié)構(gòu)模態(tài)參數(shù)識別研究[J]，建筑科學(xué)與工程學(xué)報(bào)，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]徐利,鄒傳云，陳民，等. 基于矩陣束算法的極點(diǎn)提取分析[J]. 通信技術(shù)，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)弱響應(yīng)且模態(tài)秘籍的參數(shù)識別方法[J]. 中國機(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]楊明華,奇異矩陣束的標(biāo)準(zhǔn)形與廣義逆[D]. 哈爾濱：哈爾濱工業(yè)大學(xué),2013.

[9]鄭敏,申凡,史東鋒,等.單獨(dú)利用響應(yīng)數(shù)據(jù)進(jìn)行模態(tài)分析[J]. 中國機(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]黃應(yīng)來,董大偉,閆兵.密集模態(tài)分離及其參數(shù)識別方法研究[J].機(jī)械強(qiáng)度,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

國家自然科學(xué)基金資助項(xiàng)目(11302175)

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

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