CMIF方法中模態(tài)參數(shù)不確定性的計(jì)算

孫鑫暉，郝木明，李振濤，張令彌，王 彤

（1．中國石油大學(xué)（華東）化學(xué)工程學(xué)院，山東 青島 266580；2．南京航空航天大學(xué)振動工程研究所，江蘇 南京 210016）

引 言

模態(tài) 參 數(shù) 識 別 （Modal Parameter Identification，MPI）是試驗(yàn)?zāi)B(tài)分析（Experimental Modal Analysis，EMA）的重要內(nèi)容。由于實(shí)驗(yàn)過程中不可避免的受到環(huán)境噪聲、儀器測量誤差等不確定因素的影響，從而導(dǎo)致模態(tài)參數(shù)識別結(jié)果中同樣存在不確定性。從目前的模態(tài)參數(shù)識別方法來看，絕大多數(shù)方法將模態(tài)參數(shù)作為一個(gè)確定性參數(shù)進(jìn)行識別，并未考慮其不確定性的影響。模態(tài)參數(shù)的不確定性信息可以有多種應(yīng)用，例如用于區(qū)分結(jié)構(gòu)模態(tài)與計(jì)算模態(tài)［1］，用于模型修正中作為加權(quán)信息使用［2］，因此一種模態(tài)參數(shù)識別方法在獲得識別結(jié)果的同時(shí)，有必要給出其不確定性信息。當(dāng)測量噪聲滿足正態(tài)分布時(shí)，模態(tài)參數(shù)的不確定性可以通過其統(tǒng)計(jì)特性進(jìn)行描述（如方差）［3］。在模態(tài)參數(shù)不確定性的研究中，時(shí)域方法：Longman對ERA方法進(jìn)行攝動分析，確定出模態(tài)參數(shù)對測量誤差的靈敏度［4］；Paez使用Bootstrap方法，對模態(tài)參數(shù)統(tǒng)計(jì)特性進(jìn)行了估計(jì)［5］；Reynders，Carden討論了隨機(jī)子空間方法中模態(tài)參數(shù)的不確定性計(jì)算［3，6］。頻域方法：Guillaume分析了頻響函數(shù)有理分式模型中系數(shù)攝動對極點(diǎn)的影響［7］；Pintelon給出了運(yùn)行模態(tài)分析中模態(tài)參數(shù)不確定性的計(jì)算方法［8］；Troyer對模態(tài)參數(shù)置信區(qū)間的快速計(jì)算方法進(jìn)行了研究［9，10］。

本文以頻域中CMIF模態(tài)參數(shù)識別方法（Complex Modal Indicator Function）為研究對象［11］，通過Kronecker代數(shù)以及矩陣的靈敏度分析，詳細(xì)推導(dǎo)了由測試不確定性（頻響函數(shù)中的噪聲方差）獲得模態(tài)參數(shù)不確定性（模態(tài)參數(shù)方差）的計(jì)算公式。

1 CMIF模態(tài)參數(shù)識別方法理論基礎(chǔ)

CMIF方法最初作為一種模態(tài)指示函數(shù)來使用，用于指示模態(tài)的分布情況，進(jìn)而發(fā)展為一種模態(tài)參數(shù)識別方法，用于EMA情況下模態(tài)參數(shù)識別。經(jīng)過發(fā)展，該方法已經(jīng)推廣到運(yùn)行模態(tài)分析（Operational Modal Analysis，OMA），用于僅有響應(yīng)數(shù)據(jù)情況下的參數(shù)識別［12］。CMIF識別方法的思想是利用復(fù)模態(tài)指示函數(shù)曲線，對固有頻率處的頻響函數(shù)進(jìn)行奇異值分解，獲得增強(qiáng)頻響函數(shù)（Enhanced FRF）。增強(qiáng)頻響函數(shù)在固有頻率附近近似為單自由度系統(tǒng)，可運(yùn)用單自由度擬合方法進(jìn)行識別。由于其擬合頻段范圍很窄，屬于一種窄帶識別算法［13］，其識別過程可以按照以下兩步進(jìn)行。

Step1：計(jì)算固有頻率處的增強(qiáng)頻響函數(shù)

對頻響函數(shù)矩陣H（jω）∈CNo×Ni（No，Ni為輸出、輸入數(shù)目）在固有頻率附近處進(jìn)行奇異值分解，可以得到

Step2：對增強(qiáng)頻響函數(shù)進(jìn)行擬合獲得模態(tài)參數(shù)

式中λ為極點(diǎn)，待識別參數(shù)b1，b0，λ可以通過以下最小二乘方程得到

其中X，θ，Y的表達(dá)式為：

式中Nf為擬合頻段內(nèi)包含的譜線數(shù)。固有頻率、阻尼比可以通過式（6）計(jì)算

2 CMIF方法中模態(tài)參數(shù)不確定性的計(jì)算

圖1 模態(tài)參數(shù)方差的計(jì)算流程Fig．1 The calculation flow of modal parameter variance

根據(jù)CMIF方法的識別過程，模態(tài)參數(shù)方差的計(jì)算流程按照圖1分4步完成。

2．1 頻響函數(shù)中噪聲方差的估計(jì)

噪聲的方差信息可以在頻響函數(shù)估計(jì)過程中獲得（如H1，H2估計(jì)），對于頻響函數(shù)矩陣H（jω）中的每一頻響函數(shù)Ho，i（jω）（o＝1，2，…，No，i＝1，2，…，Ni），方差的估計(jì)方法如下［14］

式中 var表示方差，E表示期望，M為頻響函數(shù)估計(jì)過程中的平均次數(shù)，γo，i（jω）為相干函數(shù)。

2．2 增強(qiáng)頻響函數(shù)協(xié)方差的計(jì)算

通過線性代數(shù)知識可知［15］，對于矩陣A，B，C存在如下關(guān)系

為了以后的推導(dǎo)過程中表示方便，這里引入2個(gè)下標(biāo)xre和xRe，二者的定義［8］

式中x可以為標(biāo)量，也可以為矩陣；Re（x），Im（x）表示x的實(shí)部與虛部。容易證明以上兩個(gè)符號滿足如下的關(guān)系

其中cov（［vec（H（jω））］re）∈R2NoNi×2NoNi為頻響函數(shù)矩陣按列展開后實(shí)部、虛部的協(xié)方差矩陣。

2．3 極點(diǎn)協(xié)方差的計(jì)算

對式（4）兩端左乘XH得到如下正規(guī)方程

式中M＝XHX，N＝XHY，對式（17）進(jìn)行一階靈敏度分析可以得到

將式（16）中的高階項(xiàng)略去，得到

式中

將式（18）代入式（17）得到Δθ的表達(dá)式

根據(jù)式（5）可知λ位于列向量θ的第3行，因此令P（jωk）為列向量（jωk－λ）M－1X（jωk）H的第3行，則Δλ可以表示為

將式（12）應(yīng)用于式（20），（Δλ）re可以表示為

之所以要將Δλ表示（Δλ）re的形式，原因是式（24）中計(jì)算模態(tài)參數(shù)方差時(shí)需要使用λre的協(xié)方差矩陣，而不是λ的方差。由式（21）得到λre的協(xié)方差為

假設(shè)頻響函數(shù)中噪聲在不同頻率分量之間相互獨(dú)立［9］，式（22）可以簡化為如下形式

2．4 模態(tài)參數(shù)方差的計(jì)算

由于CMIF方法在s域中識別，s域中極點(diǎn)的cov（λre）與模態(tài)參數(shù)方差之間關(guān)系可通過極點(diǎn)的靈敏度分析得到，計(jì)算公式如下，詳細(xì)推導(dǎo)見文獻(xiàn)［8］。

3 仿真算例

圖2 5自由度仿真算例Fig．2 Five DOF simulation case

采用圖2中的一個(gè)5自由度數(shù)值算例進(jìn)行仿真檢驗(yàn)。首先由表1數(shù)據(jù)計(jì)算出模態(tài)參數(shù)，如表2所示。在各階模態(tài)中添加1%阻尼比，利用模態(tài)疊加法得到理論頻響函數(shù)，頻率范圍0～3Hz，譜線數(shù)為1 024。在理論頻響函數(shù)的基礎(chǔ)上生成200個(gè)頻響函數(shù)樣本，在每個(gè)樣本中添加5%的高斯白噪聲，圖3為一個(gè)頻響函數(shù)樣本以及噪聲的標(biāo)準(zhǔn)差。通過Monte Carlo模擬（簡稱MC）對每個(gè)樣本進(jìn)行識別，得到200組模態(tài)參數(shù)，統(tǒng)計(jì)出模態(tài)參數(shù)的方差，并以此作為真值。然后由本文方法計(jì)算出模態(tài)參數(shù)的方差，將此結(jié)果與MC結(jié)果作比較。

表1 模型的物理參數(shù)Tab．1 Physical parameters

表2 模態(tài)參數(shù)的理論值Tab．2 Theoretical modal parameters

圖3 一個(gè)FRF樣本及噪聲標(biāo)準(zhǔn)差（5%噪聲）Fig．3 A sample FRF and noise standard deviation（5%noise）

圖4為5%噪聲下一個(gè)頻響函數(shù)樣本的識別結(jié)果。圖中曲線從上向下依次為頻響函數(shù)奇異值分解得到的奇異值曲線，用于指示模態(tài)。第一條奇異值曲線（最上面）的每個(gè)峰值對應(yīng)一階模態(tài)，如果存在密集模態(tài)，則第一、第二兩條奇異值曲線會同時(shí)出現(xiàn)峰值。采用MC模擬對200組頻響函數(shù)樣本進(jìn)行參數(shù)識別，結(jié)果如圖5所示。

圖4 一個(gè)頻響函數(shù)樣本的CMIF方法識別結(jié)果Fig．4 CMIF curve and identification results

圖5 Monte Carlo模擬Fig．5 Monte Carlo simulation

表3為5%噪聲情況下，頻率、阻尼方差的MC模擬結(jié)果與本文的計(jì)算結(jié)果，其中頻率標(biāo)準(zhǔn)差最大誤差為6．6%，阻尼標(biāo)準(zhǔn)差最大誤差為8．3%。另外在頻響函數(shù)中添加3%，10%的高斯白噪聲，這3種情況的MC模擬結(jié)果與計(jì)算結(jié)果如圖6所示，隨著噪聲量級的增加，模態(tài)參數(shù)的方差也在增加。在本算例中，各階模態(tài)隨著階次增加方差也在增大，原因是各階模態(tài)的強(qiáng)度依次減弱（圖4中各階模態(tài)的峰值高度），因此受到同樣量級噪聲干擾的情況下，弱模態(tài)的抗噪聲干擾能力低，導(dǎo)致其模態(tài)參數(shù)方差增加。

按照文獻(xiàn)［9］計(jì)算了5%噪聲情況下ployLSCF算法中模態(tài)參數(shù)的方差，其中識別為10，結(jié)果如表4所示。由于采用兩種不同的模態(tài)參數(shù)方法，其識別結(jié)果的方差除了與噪聲大小有關(guān)之外，還與識別算法自身有關(guān)，因此二者之間數(shù)值并不相同，但是結(jié)果的數(shù)量級以及各階模態(tài)方差的相對大小基本一致。

表3 計(jì)算結(jié)果與MC模擬結(jié)果（CMIF）Tab．3 Estimation and MC simulation results（CMIF）

表4 計(jì)算結(jié)果與MC模擬結(jié)果（polyLSCF）Tab．4 Estimation and MC simulation results（polyLSCF）

圖6 不同量級噪聲下計(jì)算結(jié)果與MC模擬結(jié)果Fig．6 Estimation and MC results for different noise level

4 實(shí)測算例

采用圖7中T型構(gòu)件進(jìn)行實(shí)驗(yàn)驗(yàn)證。實(shí)驗(yàn)由錘擊法完成，數(shù)據(jù)通過法國OROS公司的OR34動態(tài)分析儀采集。由分析儀所獲得數(shù)據(jù)為UFF58格式［16］，該數(shù)據(jù)格式廣泛應(yīng)用于結(jié)構(gòu)動力學(xué)分析中。本文中的程序在MATLAB2007下完成，為了讀取UFF格式文件，需要參照文檔進(jìn)行數(shù)據(jù)轉(zhuǎn)換［16］。

圖7 實(shí)驗(yàn)裝置與測量儀器Fig．7 Experimental setup and instruments

如圖8所示，測試重復(fù)進(jìn)行40次，每一次敲擊可以估計(jì)一組頻響函數(shù)；圖9為一個(gè)頻響函數(shù)樣本及其噪聲分布。噪聲的方差通過40組頻響函數(shù)統(tǒng)計(jì)得到，模態(tài)參數(shù)方差的計(jì)算過程與仿真算例相同。

圖8 一個(gè)測點(diǎn)多組激勵與響應(yīng)Fig．8 Multi－group input and output

圖9 實(shí)測頻響函數(shù)及其噪聲標(biāo)準(zhǔn)差Fig．9 Measurement FRF and noise standard deviation

圖10 T型構(gòu)件一個(gè)FRF樣本的識別結(jié)果Fig．10 Identification result of T shape structure

圖11 計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的比較Fig．11 Estimation and experiment results

圖10為一個(gè)頻響函數(shù)樣本的識別結(jié)果，在分析范圍內(nèi)包含7階模態(tài)。模態(tài)參數(shù)方差的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果如圖11所示，二者吻合程度相對于仿真算例有所降低，但也能夠接受。誤差主要來自兩方面：（1）頻響函數(shù)矩陣中的所有頻響函數(shù)互不相關(guān)假設(shè)不能完全滿足，因此計(jì)算噪聲的協(xié)方差矩陣有一定誤差；（2）式（23）中假設(shè)頻響函數(shù)中噪聲在不同頻率分量之間相互獨(dú)立，由圖9可知，噪聲與頻響函數(shù)之間具有一定的相關(guān)性。

5 結(jié) 論

（1）利用矩陣的靈敏度分析與Kronecker代數(shù)理論，通過測量頻響函數(shù)中噪聲的方差信息，給出了CMIF識別算法中模態(tài)參數(shù)不確定性的計(jì)算方法，在識別結(jié)果的基礎(chǔ)上增加了統(tǒng)計(jì)信息。

（2）仿真結(jié)果表明，推導(dǎo)過程中采用的一階靈敏度分析方法已經(jīng)具有較高的精度。實(shí)測算例中由于一些假設(shè)不能完全滿足，計(jì)算精度相對于仿真結(jié)果有所降低。

（3）在工程應(yīng)用中，如果原始數(shù)據(jù)只有頻響函數(shù)，沒有噪聲的方差信息，這時(shí)模態(tài)參數(shù)方差的計(jì)算將受到限制。

［1］ Verboven P，Cauberghe B．User－assisting tools for a fast frequency－domain modal parameter estimation method［J］．MSSP，2004，18（4）：759—780．

［2］ Cauberghe B．Applied frequency－domain system identification in the field of experimental and operational modal analysis［D］．Belgium：Vrije Universiteit Brussel，2004．

［3］ Reynders E，Pintelon R，Roecka G．Uncertainty bounds on modal parameters obtained from stochastic subspace identification［J］．Mechanical Systems and Signal Processing，2008，22（4）：948—969．

［4］ Longman R W，Juang J N．Variance Based Confidence Criterion for ERA Identification Modal Parameters［R］．AAS 87－454：581—601．

［5］ Paez T L，Hunter N F．Statistcal series：part25：fundamental concepts of the bootstrap for statistical analysis of mechanical systems［J］．Experimental Techniques，1998，22（3）：3 523．

［6］ Carden E P，Mita A．Challenges in developing confidence intervals on modal parameters estimated for large civil infrastructure with stochastic subspace identification［J］．Structural Control and Health Monitoring，2009，23（1）：217—225．

［7］ Guillaume P，Schoukens J，Pintelon R．Sensitivity of roots to errors in the coefficients of polynomials ob－tained by frequency domain estimation methods［J］．IEEE Transactions on Instrumentation and Measurement，1989，38（6）：1 050—1 056．

［8］ Pintelon R，Guillaume P，Schoukens J．Uncertainty calculation in（operational）modal analysis［J］．MSSP，2007，21（6）：2 359—2 373．

［9］ Troyer T，Guillaume P，Steenackers G．Fast variance calculation of polyreference least－squares frequency－domain estimates［J］．MSSP，2009，23（5）：1 423—1 433．

［10］Troyer T，Guillaume P，Pintelon R，et al．Fast calculation of confidence intervals on parameter estimates of least－squares frequency－domain estimators［J］．MSSP，2009，23（2）：261—273．

［11］Shih C Y，Tsuei Y G．Complex mode indication function and its applications to spatial domain parameter estimation［J］．MSSP，1988，2（4）：367—377．

［12］Zhang L M，Wang T，Tamura Y．A frequency－spatial domain decomposition（FSDD）method for operational modal analysis［J］．MSSP，2010，24（5）：1 227—1 239．

［13］Zhang L M，Sun X H，Wang T．Narrow－band，selectband vs broadband modal identification：their features and comparisons［A］．Proc of the 26th International Modal Analysis Conference［C］．USA，F(xiàn)ebruary 2009．

［14］貝達(dá)特，皮爾索，著．相關(guān)分析和譜分析的工程應(yīng)用［M］．令福根，譯．北京：國防工業(yè)出版社，1983．

Bendat J S，Piersol A G．Engineering Applications of Correlation and Spectral Analysis［M］．Beijing：National Defence Industry Press，1983（in Chinese）．

［15］戴華．矩陣論［M］．北京：科學(xué)出版社，2001．Dai H．Matrix Theory［M］．Beijing：Science Press，2001．

［16］UC－SDRL．Universal File Formats for Modal Analysis Testing［EB／OL］．http：／／www．sdrl．uc．edu／universal－file－formats－for－modal－analysis－testing－1，2008．