基于隨機(jī)子空間識別的模態(tài)參數(shù)不確定性量化方法

羅 杰, 康 杰, 孫嘉寶, 曾舒洪

(南京航空航天大學(xué) 航天學(xué)院,南京 211106)

模態(tài)參數(shù)識別是結(jié)構(gòu)健康監(jiān)測中常用的技術(shù),通過分析結(jié)構(gòu)的振動響應(yīng)信號,識別出結(jié)構(gòu)的特征參數(shù),如頻率、阻尼、模態(tài)振型等[1-3]。協(xié)方差驅(qū)動隨機(jī)子空間識別(covariance-driven stochastic subspace identification,SSI-COV)方法作為模態(tài)參數(shù)識別方法的一種,利用其強(qiáng)魯棒性、高精度等優(yōu)勢,被廣泛應(yīng)用于僅有輸出響應(yīng)數(shù)據(jù)的工作模態(tài)分析中,用于提取結(jié)構(gòu)的模態(tài)頻率、阻尼比和模態(tài)振型[4-5]。SSI-COV的計(jì)算模態(tài)參數(shù)過程可分為以下幾個步驟:首先,收集傳感器信號,將每一時(shí)刻的響應(yīng)按列排列,估計(jì)相關(guān)函數(shù)矩陣,并組成Hankel矩陣的形式;然后對Hankel矩陣進(jìn)行奇異值分解(singular value decomposition,SVD)運(yùn)算,根據(jù)分解的結(jié)果計(jì)算出狀態(tài)矩陣O,從而獲得系統(tǒng)矩陣A和輸出矩陣C;最后,對系統(tǒng)矩陣A進(jìn)行特征值分解運(yùn)算,并將分解的特征值、特征向量轉(zhuǎn)化成連續(xù)時(shí)間系統(tǒng)狀態(tài)矩陣的特征值和特征向量,根據(jù)連續(xù)系統(tǒng)的特征值和特征向量獲得結(jié)構(gòu)的模態(tài)頻率、振型和阻尼比等結(jié)構(gòu)模態(tài)參數(shù)。

由于振動測量噪聲、數(shù)據(jù)長度有限等因素的存在,SSI-COV方法在實(shí)際應(yīng)用中,模態(tài)參數(shù)識別結(jié)果會受到不確定性的影響,因此需要對模態(tài)參數(shù)的不確定度進(jìn)行量化[6]。獲得的不確定度可用于評估模態(tài)參數(shù)估計(jì)結(jié)果的質(zhì)量,并能夠剔除虛假模態(tài);還用于考慮結(jié)構(gòu)動力學(xué)參數(shù)不確定性的結(jié)構(gòu)設(shè)計(jì)、模型修正等問題的研究[7-8]。Reynders等[9-10]利用矩陣一階擾動分析方法量化了SSI-COV方法識別模態(tài)參數(shù)的不確定度,并給出了具體的公式推導(dǎo)過程。Carden等[11]利用SSI估計(jì)了大型民用基礎(chǔ)設(shè)施模態(tài)參數(shù)的置信區(qū)間。這些方法在計(jì)算各中間變量的擾動時(shí),均采用矩陣?yán)币约癒ronecker積等數(shù)學(xué)運(yùn)算來推導(dǎo)出模態(tài)參數(shù)的方差。這種計(jì)算方法會導(dǎo)致矩陣維度高、運(yùn)算效率低等問題。

針對該問題,本文提出了一種基于SSI-COV的模態(tài)參數(shù)不確定度高效計(jì)算方法。首先,計(jì)算振動響應(yīng)相關(guān)函數(shù)的方差,通過SVD,選取恰當(dāng)?shù)钠娈愔到財(cái)嚯A數(shù),由每階奇異向量組裝,獲得多組Hankel矩陣的擾動。其次,根據(jù)一階矩陣攝動理論,隱式計(jì)算SSI-COV算法中各中間變量的一階擾動,最終,由多組模態(tài)參數(shù)的擾動疊加計(jì)算出方差。

本文組織結(jié)構(gòu)如下:第1節(jié)概述SSI-COV方法原理和計(jì)算過程;第2節(jié)介紹本文提出的不確定度高效計(jì)算方法,并與已有方法相比,討論了計(jì)算量的不同;第3節(jié)利用桁架結(jié)構(gòu)數(shù)據(jù)對提出的方法進(jìn)行蒙特卡洛仿真(Monte Carlo simulation,MCS)驗(yàn)證,并對上一節(jié)的計(jì)算量對比分析具體化,第4節(jié)對本文進(jìn)行總結(jié)。

1 SSI-COV方法原理

1.1 離散時(shí)間狀態(tài)空間模型

具有N個自由度系統(tǒng),在外界環(huán)境激勵下的運(yùn)動微分方程可以表示為

(1)

式中:M為質(zhì)量矩陣;C為阻尼矩陣;K為剛度矩陣;q(t)為t時(shí)刻的位移列向量;f(t)為外部激勵向量。

在實(shí)際測試中,考慮到信號采集的時(shí)間離散性和測量時(shí)存在的噪聲干擾現(xiàn)象,對式(1)進(jìn)行采樣離散化處理和微分求解,得到的系統(tǒng)離散時(shí)間狀態(tài)空間模型為[12]

(2)

式中:yk和xk分別為第k個采樣時(shí)刻系統(tǒng)的輸出向量和狀態(tài)向量;wk和vk分別為環(huán)境干擾噪聲序列和傳感器自身測量誤差所帶來的噪聲序列;A為系統(tǒng)狀態(tài)矩陣;C為系統(tǒng)輸出矩陣。

此空間模型是由狀態(tài)方程和輸出方程組合,構(gòu)成一個系統(tǒng)完整的動態(tài)描述。其基本假設(shè)前提為

(1)噪聲項(xiàng)為均值為零的平穩(wěn)隨機(jī)過程,即

E[wk]=0, E[vk]=0

(3)

(2) 噪聲信號是均值為零的白噪聲序列,且與結(jié)構(gòu)真實(shí)狀態(tài)無關(guān),即

(4)

(3) 系統(tǒng)為線性時(shí)不變系統(tǒng)(既滿足疊加原理又 具有時(shí)不變特性),狀態(tài)序列為平穩(wěn)隨機(jī)過程,即

(5)

式中: E為數(shù)學(xué)期望;Σ為狀態(tài)協(xié)方差矩陣;k為時(shí)間。

1.2 模態(tài)參數(shù)識別

根據(jù)系統(tǒng)的輸出響應(yīng)求出相關(guān)函數(shù)矩陣,在實(shí)際測試中,通常只選取一段數(shù)據(jù)量參與計(jì)算,因此,輸出協(xié)方差矩陣可估計(jì)為[13]

(6)

式中:yk∈m×1;τ為相關(guān)函數(shù)時(shí)延;S為數(shù)據(jù)點(diǎn)個數(shù);k為第k個時(shí)刻且k=1,2,3,…,S;m為輸出通道數(shù)。

根據(jù)式(6)相關(guān)函數(shù)構(gòu)造Hankel矩陣為

(7)

式中,i和j分別為Hankel矩陣的列塊數(shù)和行塊數(shù)。對Hankel矩陣進(jìn)行SVD,得

(8)

式中:U為左奇異矩陣;Σ為奇異值矩陣;V為右奇異矩陣。

(9)

則系統(tǒng)矩陣A和輸出矩陣C可分別表示為

(10)

式中:O1:m為可觀矩陣O的前m行形成的新矩陣;O↑為去掉矩陣O后m行形成的新矩陣;O↓為去掉矩陣O前m行形成的新矩陣。

由矩陣A特征值分解得

(11)

式中:λ為特征值;ψ為右特征向量;β為左特征向量;H為取復(fù)共軛轉(zhuǎn)置。

則第d階結(jié)構(gòu)模態(tài)頻率fd、阻尼比ξd和模態(tài)振型φd可分別表示為

1.3 Hankel矩陣維度及模型階數(shù)確定

在對結(jié)構(gòu)進(jìn)行模態(tài)辨識之前,需要選取一系列模型參數(shù)。其中Hankel矩陣維度和模型階數(shù)的選擇對隨機(jī)子空間法的辨識結(jié)果影響十分顯著,因此針對不同的系統(tǒng)結(jié)構(gòu),選擇不同的參數(shù)具有重要意義[15]。Hankel矩陣維度可根據(jù)下列公式確定

j+i-1≤S

(13)

(14)

(15)

式中,fs和f0分別為采樣頻率和基頻。

2 模態(tài)參數(shù)不確定度量化

由于辨識的模態(tài)參數(shù)總是受到許多不同來源的統(tǒng)計(jì)不確定度的影響,如響應(yīng)測量噪聲與數(shù)據(jù)有限長度,因此,需要對已辨識的模態(tài)參數(shù)進(jìn)行不確定度量化。最開始計(jì)算不確定度的方法需顯式表示出Jacobian矩陣,導(dǎo)致矩陣運(yùn)算維度高、計(jì)算效率低。文獻(xiàn)[9]提出了一種無需顯式計(jì)算Jacobian矩陣的方法,本文采用該方法進(jìn)行模態(tài)參數(shù)的不確定度量化,基于一階矩陣攝動理論逐級計(jì)算誤差擾動,估計(jì)模態(tài)參數(shù)的方差。

2.1 Hankel矩陣的擾動估計(jì)

(16)

式中,vec(·)為矩陣?yán)边\(yùn)算。對cov(T)進(jìn)行SVD,則cov(T)可近似為

(17)

(18)

由矩陣T得到的擾動矩陣與矩陣T大小保持一致,所以將擾動向量轉(zhuǎn)換成jm×im矩陣,即為矩陣T的第ε階擾動,記作ΔTε,其中ε=1,2,3,…,n。

由于式(16)中Hankel矩陣T的方差矩陣維度為ijm2×ijm2,維度高,對計(jì)算機(jī)內(nèi)存的需求較大。觀察可知,Hankel矩陣T中存在很多重復(fù)的相關(guān)函數(shù)矩陣塊,利用該特點(diǎn)可以大幅度降低計(jì)算量。

提取其中一個相關(guān)函數(shù)矩陣塊,記為

(19)

(20)

對方差矩陣cov(R)進(jìn)行SVD,文獻(xiàn)[16]需要進(jìn)行nb次運(yùn)算,為了進(jìn)一步降低計(jì)算量,可選擇奇異值階數(shù)q進(jìn)行截?cái)?得到R矩陣的q階擾動為

(21)

(22)

將擾動向量Δδr重組成大小為m×(i+j-1)m的矩陣,再按式(7)中排列順序填入對應(yīng)的位置,即可獲得Hankel矩陣T的第r階擾動矩陣,記作ΔTr,其中r=1,2,…,q。

上述方法避免了直接計(jì)算Hankel矩陣的方差矩陣cov(T),并且進(jìn)行奇異值截?cái)?降低計(jì)算量的同時(shí)大幅節(jié)省了存儲方差矩陣的內(nèi)存空間。

2.2 模態(tài)參數(shù)的方差估計(jì)

首先,根據(jù)一階矩陣擾動傳播理論,由式(8)可知,Hankel矩陣將擾動傳遞給奇異值和奇異向量。

(23)

(24)

考慮T矩陣存在擾動ΔT,由式(24)可得

ΔTvc+TΔvc=Δσcuc+σcΔuc

(25)

(26)

(27)

(28)

將式(24)左右奇異向量分離變量,考慮矩陣T存在擾動ΔT,可得

(29)

將式(24)、式(28)代入式(29),得

(30)

式中:I為單位矩陣;T為矩陣轉(zhuǎn)置。

針對式(17)得到的每階擾動,均可計(jì)算出矩陣T奇異值、奇異向量的對應(yīng)擾動。由矩陣T第r階擾動引起的第l階奇異值、奇異向量的擾動為

(31)

則由矩陣T第r階擾動引起的奇異值、奇異向量矩陣的擾動可記為

(32)

考慮O矩陣存在擾動,由式(9)可得狀態(tài)矩陣O的擾動為[15]

(33)

考慮輸出矩陣C、系統(tǒng)矩陣A存在擾動,由式(10)可得[15]

ΔCr=S3ΔOr

(34)

(35)

式中:

S1=[I(j-1)m×(j-1)m0(j-1)m×m]

S2=[0(j-1)m×mI(j-1)m×(j-1)m]

S3=[Im×m0m×(j-1)m]

考慮特征值、特征向量存在擾動,由式(11)可得

(36)

(37)

根據(jù)文獻(xiàn)[9],由第r階特征值擾動引起的第d階模態(tài)頻率、阻尼的擾動為

(38)

式中:μd=ln(λd)/Δt;Im(·)為取復(fù)數(shù)的虛部。

同理,由第r階特征向量的擾動引起的第d個模態(tài)振型的擾動為[9]

Δφd,r=ΔCrαd+CΔαd,r

(39)

但是,通常結(jié)構(gòu)模態(tài)振型為復(fù)數(shù)形式,需要進(jìn)行歸一化處理。歸一化選擇對模態(tài)振型上的每一個點(diǎn)除以任意選擇的點(diǎn)φd,w進(jìn)行歸一化,則模態(tài)振型的擾動則進(jìn)一步轉(zhuǎn)化為[9]

(40)

式中,S4∈1×m為在w處為1其他處均為0的行向量。

根據(jù)式(38)、式(40)可得由第r階矩陣T的擾動引起的第d階頻率、阻尼、模態(tài)振型的方差為

(41)

借助于式(41)計(jì)算出的模態(tài)參數(shù)方差,式(21)中的擾動截?cái)嚯A數(shù)q可根據(jù)如下準(zhǔn)則確定

(42)

式中,r為截?cái)嚯A數(shù)。滿足該準(zhǔn)則的最小r值即為最終確定的擾動截?cái)嚯A數(shù)q。

在已有方法文獻(xiàn)[14]中,通常計(jì)算T矩陣的擾動,先將數(shù)據(jù)響應(yīng)分段分別計(jì)算矩陣T,利用每一段矩陣T和所有矩陣T的均值相減直接獲得其擾動。而在第2.1節(jié)和本節(jié)的方法中,先計(jì)算相關(guān)函數(shù)的方差矩陣,進(jìn)行SVD,并選擇適當(dāng)?shù)钠娈愔到財(cái)嚯A數(shù),可有效降低擾動矩陣計(jì)算的次數(shù),最后由奇異值、奇異向量計(jì)算出的擾動矩陣來組裝T矩陣的擾動。因此,本文在計(jì)算T矩陣的擾動以及最終模態(tài)參數(shù)的擾動時(shí)更加高效。另外,在文獻(xiàn)[16]中,擾動傳播表達(dá)式的推導(dǎo)均利用Kronecker積和矩陣?yán)边\(yùn)算,將擾動寫在方程最右側(cè),例如式(43)寫為

Δσc=(vc?uc)Tvec(ΔT)

(43)

由于上式中包含Kronecker積和矩陣?yán)边\(yùn)算,vc?uc的維度為1×ijm2,需要ijm2次乘法計(jì)算,Δσc總的計(jì)算量包含2ijm2次乘法和(ijm2-1)次加法;而式(28)中矩陣維度更小,總的計(jì)算量為ijm2+jm乘法和(ijm2-1)次加法,計(jì)算量更小。

3 桁架結(jié)構(gòu)仿真算例驗(yàn)證

本節(jié)采用文獻(xiàn)[17]中的桁架結(jié)構(gòu)模型驗(yàn)證所提方法的有效性。桁架結(jié)構(gòu)模型如圖1所示,彈性模量為6.98×1010Pa,材料密度為2 770 kg/m3,并在節(jié)點(diǎn)①、節(jié)點(diǎn)②、節(jié)點(diǎn)③、節(jié)點(diǎn)④上加入454 kg的額外集中質(zhì)量,Ai為第i根桿件橫截面積,其阻尼矩陣與質(zhì)量矩陣成正比且第一模態(tài)阻尼等于1%,表1和表2給出了桁架結(jié)構(gòu)真實(shí)模態(tài)頻率和阻尼比。

圖1 桁架結(jié)構(gòu)圖Fig.1 Truss structure

表1 模態(tài)頻率估計(jì)值的標(biāo)準(zhǔn)差Tab.1 Standard deviation of frequency estimates

表2 模態(tài)阻尼比估計(jì)值的標(biāo)準(zhǔn)差Tab.2 Standard deviation of damping estimates

本算例中采用三種類型的激勵(即白噪聲、諧波激勵和非白噪聲激勵),其中諧波激勵的頻率為12 Hz;非白噪聲激勵由白噪聲經(jīng)過單自由度系統(tǒng)過濾得到,濾波函數(shù)為:

(44)

式中,s為Laplace變換算子,單自由度系統(tǒng)的阻尼比ζ0=0.5%,固有頻率ω0=40 Hz。

本文在節(jié)點(diǎn)①的x、y方向同時(shí)作用諧波激勵與非白噪聲激勵,且沿節(jié)點(diǎn)①的x、y方向的激勵完全相關(guān),節(jié)點(diǎn)②、節(jié)點(diǎn)③、節(jié)點(diǎn)④的x、y方向同時(shí)作用不相關(guān)高斯白噪聲激勵。使用Newmark-β方法計(jì)算桁架結(jié)構(gòu)的位移響應(yīng),響應(yīng)信號所加信噪比為20 dB,時(shí)間步長為1/1 024 s,總時(shí)長320 s;隨后進(jìn)行128 Hz重采樣,因此位移響應(yīng)數(shù)據(jù)點(diǎn)數(shù)為40 960。

3.1 仿真結(jié)果

利用SSI-COV方法估計(jì)桁架結(jié)構(gòu)模態(tài)參數(shù),并根據(jù)本文所提方法計(jì)算模態(tài)參數(shù)標(biāo)準(zhǔn)差。將得到的模態(tài)參數(shù)標(biāo)準(zhǔn)差與MCS方法得到的模態(tài)參數(shù)標(biāo)準(zhǔn)差進(jìn)行比較,從而驗(yàn)證本文方法的有效性和可行性。

圖2 Hankel矩陣奇異值變化曲線Fig.2 Singular value curve of Hankel matrix

對于相關(guān)函數(shù)方差矩陣奇異值截?cái)嚯A數(shù)q的選取,本文提取了第4個模態(tài)的頻率方差進(jìn)行計(jì)算,其截?cái)嚯A數(shù)為20時(shí)對應(yīng)的頻率方差為4.699 66×10-05,截?cái)嚯A數(shù)為21時(shí)對應(yīng)的頻率方差為4.707 37×10-05,經(jīng)式(42)計(jì)算,當(dāng)r=20時(shí)的截?cái)嘞禂?shù)為0.001 6,同理,當(dāng)r=18時(shí)計(jì)算出的截?cái)嘞禂?shù)為0.12,r=22時(shí)計(jì)算出的截?cái)嘞禂?shù)為0.000 7,故確定出有效截?cái)嚯A數(shù)q為20。如圖3所示,顯示的四條曲線表示系統(tǒng)前四階模態(tài)頻率的方差,橫軸表示相關(guān)函數(shù)方差矩陣進(jìn)行SVD所截?cái)嗟碾A數(shù),縱軸表示估計(jì)的模態(tài)頻率方差。不難發(fā)現(xiàn),隨著截?cái)嚯A數(shù)的增加,約20階以后,估計(jì)的模態(tài)頻率方差開始趨于穩(wěn)定。所以,本文選擇q=20作為相關(guān)函數(shù)方差矩陣的截?cái)嚯A數(shù),不僅能確保模態(tài)參數(shù)不發(fā)生較大偏差,還降低了同原來方法一半的計(jì)算量,大幅提高了原有算法的計(jì)算效率。

圖3 奇異值截?cái)嚯A數(shù)對模態(tài)頻率方差的影響Fig.3 The effect of singular value truncation order on modal frequency variance

確定固定參數(shù)后,進(jìn)行500次MCS,通過SSI-COV方法估計(jì)桁架模態(tài)參數(shù)。估計(jì)的第一階模態(tài)頻率、阻尼比如圖4和圖5所示。結(jié)果顯示,所有離散的樣本點(diǎn)都在數(shù)值1上下波動,辨識值與真實(shí)值接近,并未出現(xiàn)異常的辨識結(jié)果,說明仿真結(jié)果是可靠的。

圖4 第一階模態(tài)頻率辨識結(jié)果Fig.4 First-order modal frequency identification results

圖5 第一階模態(tài)阻尼辨識結(jié)果Fig.5 First-order modal damping identification results

其次,將所提方法計(jì)算的模態(tài)參數(shù)標(biāo)準(zhǔn)差與MCS得到的樣本統(tǒng)計(jì)量進(jìn)行比較。在每次仿真中,通過第2.2節(jié)中提出的方法估計(jì)模態(tài)頻率、阻尼比和模態(tài)振型的標(biāo)準(zhǔn)差,隨后計(jì)算500次估計(jì)標(biāo)準(zhǔn)差的均值。將此標(biāo)準(zhǔn)差均值與500次MCS得到的樣本標(biāo)準(zhǔn)差進(jìn)行比較,模態(tài)頻率和阻尼比結(jié)果如表1和表2所示,模態(tài)振型對比結(jié)果如圖6所示。

圖6 模態(tài)振型標(biāo)準(zhǔn)差Fig.6 Standard deviation of modal shape

結(jié)果顯示,SSI-COV方法能夠完整地估計(jì)8階結(jié)構(gòu)模態(tài),識別結(jié)果中頻率為12 Hz與40 Hz的兩階模態(tài)是由諧波激勵和有色噪聲激勵引起的,屬于虛假模態(tài)。由此發(fā)現(xiàn),SSI-COV方法在模態(tài)參數(shù)識別過程中無法剔除虛假模態(tài)。但是,由計(jì)算獲得的模態(tài)參數(shù)不確定度結(jié)果表明,對于由非白噪聲激勵或諧波激勵引起的已識別的虛假模態(tài),其頻率和阻尼比的不確定度顯著大于結(jié)構(gòu)模態(tài)的不確定度,可以利用該特性來完成對虛假模態(tài)的剔除。且所提方法計(jì)算的模態(tài)參數(shù)標(biāo)準(zhǔn)差總體小于MCS中模態(tài)參數(shù)樣本計(jì)算的標(biāo)準(zhǔn)差。

對于模態(tài)振型的標(biāo)準(zhǔn)差,選擇了第1、第2階的結(jié)構(gòu)模態(tài)振型不確定度和虛假模態(tài)(諧波激勵、非白噪聲激勵)的模態(tài)振型不確定度結(jié)果進(jìn)行比較。如圖6所示,實(shí)線表示由本文方法計(jì)算的不確定度結(jié)果,虛線則由MCS模態(tài)參數(shù)樣本計(jì)算的不確定度結(jié)果。通過觀察,由本文方法計(jì)算的模態(tài)振型的不確定度同樣比仿真結(jié)果偏大,且虛假模態(tài)計(jì)算的振型標(biāo)準(zhǔn)差顯著大于結(jié)構(gòu)模態(tài)的振型標(biāo)準(zhǔn)差,由振型標(biāo)準(zhǔn)差的結(jié)果,可進(jìn)一步作為判斷虛假模態(tài)的依據(jù)。

下面通過討論不同Hankel維度對所辨識的模態(tài)參數(shù)不確定度的影響。文中討論了Hankel矩陣行分塊數(shù)為30、60、90等三種不同情況對辨識出的模態(tài)參數(shù)不確定度的影響,除矩陣行塊數(shù)j以外其余的固定參數(shù)設(shè)置保持不變。

模態(tài)頻率、阻尼的不確定度如圖7和圖8所示,觀察可知,隨著Hankel矩陣維度的變化,由諧波激勵產(chǎn)生的模態(tài)頻率、阻尼的不確定度變化幅度較大且沒有規(guī)律,該特點(diǎn)也可作為判斷虛假模態(tài)依據(jù)。由非白噪聲產(chǎn)生的模態(tài)頻率不確定度總體上與結(jié)構(gòu)真實(shí)模態(tài)頻率的不確定度變化一致,隨著Hankel矩陣維度增加而增大,但非白噪聲激勵引起的模態(tài)頻率對應(yīng)的不確定度顯著大于結(jié)構(gòu)模態(tài)頻率的不確定度。而對于阻尼模態(tài)頻率的不確定度來說,結(jié)構(gòu)模態(tài)的不確定度隨Hankel矩陣維度變化很小,虛假模態(tài)的不確定度變化較為明顯。

圖7 不同Hankel矩陣維度對模態(tài)頻率不確定度的影響Fig.7 Influence of different Hankel matrix dimensions on modal frequency uncertainty

圖8 不同Hankel矩陣維度對模態(tài)阻尼不確定度的影響Fig.8 Influence of different Hankel matrix dimensions on modal damping uncertainty

3.2 不確定度量化的計(jì)算量分析

下面就本文方法和傳統(tǒng)方法的計(jì)算量作具體對比,分析在實(shí)際計(jì)算過程中,兩者計(jì)算量上存在的差距,從而驗(yàn)證本文所提方法的有效性。以下討論的固定參數(shù)同仿真分析數(shù)值一致。

傳統(tǒng)方法中,例如文獻(xiàn)[12]的方法,計(jì)算T矩陣擾動,先將數(shù)據(jù)響應(yīng)分段分別計(jì)算矩陣T,利用每一段的矩陣T和所有矩陣T的均值相減直接獲得其擾動。根據(jù)響應(yīng)數(shù)據(jù)分段次數(shù)為40,所以還需計(jì)算40次T的擾動,整個模態(tài)參數(shù)的擾動計(jì)算也隨之計(jì)算40次。在文獻(xiàn)[11]中計(jì)算各中間變量的不確定度時(shí),如式(43)計(jì)算奇異值的擾動,vc?uc的維度為1×38 400,需要38 400次乘法計(jì)算,Δσc總的計(jì)算量包含76 800次乘法和38 399次加法。

在本文所提方法中,利用式(19)、式(20)先求相關(guān)函數(shù)矩陣塊R的方差,進(jìn)行SVD,根據(jù)截?cái)鄿?zhǔn)則在截?cái)嚯A數(shù)為20時(shí)進(jìn)行截?cái)?。而已有文獻(xiàn)[15]中數(shù)據(jù)分段數(shù)nb為40,所以擾動T原本需要40次運(yùn)算現(xiàn)縮減至20次。同理,各中間變量的擾動計(jì)算次數(shù)也隨之減少一半,整個模態(tài)參數(shù)不確定度計(jì)算效率提高50%。另外,繞過矩陣?yán)币约癒ronecker積等數(shù)學(xué)運(yùn)算,利用式(28)中計(jì)算Δσc,總的計(jì)算量為38 640次乘法和38 399次加法,與上述過程相比,減少了一半的乘法,提高了計(jì)算效率。同理,在計(jì)算各中間變量的每一階擾動時(shí),本文所提方法跟傳統(tǒng)方法相比,均減小了一定的計(jì)算量。

所以,在整個模態(tài)參數(shù)不確定度計(jì)算過程中,本文提出的方法大幅地提高了原有算法的計(jì)算效率、降低了對內(nèi)存的需求。

4 結(jié) 論

本文提出了一種基于SSI-COV模態(tài)參數(shù)不確定度高效計(jì)算方法,通過相關(guān)函數(shù)矩陣的方差進(jìn)行SVD,根據(jù)奇異值截?cái)嚯A數(shù)準(zhǔn)則,選擇合適的截?cái)嚯A數(shù)q,可顯著降低計(jì)算量,然后由各階擾動向量組裝成q組Hankel矩陣的擾動。根據(jù)一階矩陣擾動傳播理論和SSI算法,繞過傳統(tǒng)的Kronecker積和矩陣?yán)钡葦?shù)學(xué)運(yùn)算,逐級計(jì)算各中間變量的擾動,從而獲得多組模態(tài)參數(shù)的擾動,通過擾動平方求和計(jì)算模態(tài)參數(shù)的方差。

最后應(yīng)用本文方法對桁架結(jié)構(gòu)算例仿真,通過比較結(jié)構(gòu)模態(tài)和虛假模態(tài)的方差,證明了不確定度計(jì)算可用來剔除模態(tài)參數(shù)識別中的虛假模態(tài)。與已有方法比較,本文所提方法能夠大幅地提高了計(jì)算效率、降低了對內(nèi)存的需求。