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

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

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

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

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

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

1 SSI-COV方法原理

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

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

(1)

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

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

(2)

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

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

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

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

(3)

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

(4)

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

(5)

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

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

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

(6)

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

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

(7)

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

(8)

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

(9)

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

(10)

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

由矩陣A特征值分解得

(11)

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

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

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

在對結(jié)構(gòu)進行模態(tài)辨識之前,需要選取一系列模型參數(shù)。其中Hankel矩陣維度和模型階數(shù)的選擇對隨機子空間法的辨識結(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)計不確定度的影響,如響應測量噪聲與數(shù)據(jù)有限長度,因此,需要對已辨識的模態(tài)參數(shù)進行不確定度量化。最開始計算不確定度的方法需顯式表示出Jacobian矩陣,導致矩陣運算維度高、計算效率低。文獻[9]提出了一種無需顯式計算Jacobian矩陣的方法,本文采用該方法進行模態(tài)參數(shù)的不確定度量化,基于一階矩陣攝動理論逐級計算誤差擾動,估計模態(tài)參數(shù)的方差。

2.1 Hankel矩陣的擾動估計

(16)

式中,vec(·)為矩陣拉直運算。對cov(T)進行SVD,則cov(T)可近似為

(17)

(18)

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

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

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

(19)

(20)

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

(21)

(22)

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

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

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

首先,根據(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)得到的每階擾動,均可計算出矩陣T奇異值、奇異向量的對應擾動。由矩陣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ù)文獻[9],由第r階特征值擾動引起的第d階模態(tài)頻率、阻尼的擾動為

(38)

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

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

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

(39)

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

(40)

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

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

(41)

借助于式(41)計算出的模態(tài)參數(shù)方差,式(21)中的擾動截斷階數(shù)q可根據(jù)如下準則確定

(42)

式中,r為截斷階數(shù)。滿足該準則的最小r值即為最終確定的擾動截斷階數(shù)q。

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

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

(43)

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

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

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

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

表1 模態(tài)頻率估計值的標準差Tab.1 Standard deviation of frequency estimates

表2 模態(tài)阻尼比估計值的標準差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é)點①的x、y方向同時作用諧波激勵與非白噪聲激勵,且沿節(jié)點①的x、y方向的激勵完全相關(guān),節(jié)點②、節(jié)點③、節(jié)點④的x、y方向同時作用不相關(guān)高斯白噪聲激勵。使用Newmark-β方法計算桁架結(jié)構(gòu)的位移響應,響應信號所加信噪比為20 dB,時間步長為1/1 024 s,總時長320 s;隨后進行128 Hz重采樣,因此位移響應數(shù)據(jù)點數(shù)為40 960。

3.1 仿真結(jié)果

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

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

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

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

確定固定參數(shù)后,進行500次MCS,通過SSI-COV方法估計桁架模態(tài)參數(shù)。估計的第一階模態(tài)頻率、阻尼比如圖4和圖5所示。結(jié)果顯示,所有離散的樣本點都在數(shù)值1上下波動,辨識值與真實值接近,并未出現(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

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

圖6 模態(tài)振型標準差Fig.6 Standard deviation of modal shape

結(jié)果顯示,SSI-COV方法能夠完整地估計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)。但是,由計算獲得的模態(tài)參數(shù)不確定度結(jié)果表明,對于由非白噪聲激勵或諧波激勵引起的已識別的虛假模態(tài),其頻率和阻尼比的不確定度顯著大于結(jié)構(gòu)模態(tài)的不確定度,可以利用該特性來完成對虛假模態(tài)的剔除。且所提方法計算的模態(tài)參數(shù)標準差總體小于MCS中模態(tài)參數(shù)樣本計算的標準差。

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

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

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

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

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

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

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

4 結(jié) 論

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

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