基于隨機子空間識別的模態參數不確定性量化方法

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

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

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

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

針對該問題,本文提出了一種基于SSI-COV的模態參數不確定度高效計算方法。首先,計算振動響應相關函數的方差,通過SVD,選取恰當的奇異值截斷階數,由每階奇異向量組裝,獲得多組Hankel矩陣的擾動。其次,根據一階矩陣攝動理論,隱式計算SSI-COV算法中各中間變量的一階擾動,最終,由多組模態參數的擾動疊加計算出方差。

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

1 SSI-COV方法原理

1.1 離散時間狀態空間模型

具有N個自由度系統,在外界環境激勵下的運動微分方程可以表示為

(1)

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

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

(2)

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

此空間模型是由狀態方程和輸出方程組合,構成一個系統完整的動態描述。其基本假設前提為

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

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

(3)

(2) 噪聲信號是均值為零的白噪聲序列,且與結構真實狀態無關,即

(4)

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

(5)

式中: E為數學期望;Σ為狀態協方差矩陣;k為時間。

1.2 模態參數識別

根據系統的輸出響應求出相關函數矩陣,在實際測試中,通常只選取一段數據量參與計算,因此,輸出協方差矩陣可估計為[13]

(6)

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

根據式(6)相關函數構造Hankel矩陣為

(7)

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

(8)

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

(9)

則系統矩陣A和輸出矩陣C可分別表示為

(10)

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

由矩陣A特征值分解得

(11)

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

則第d階結構模態頻率fd、阻尼比ξd和模態振型φd可分別表示為

1.3 Hankel矩陣維度及模型階數確定

在對結構進行模態辨識之前,需要選取一系列模型參數。其中Hankel矩陣維度和模型階數的選擇對隨機子空間法的辨識結果影響十分顯著,因此針對不同的系統結構,選擇不同的參數具有重要意義[15]。Hankel矩陣維度可根據下列公式確定

j+i-1≤S

(13)

(14)

(15)

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

2 模態參數不確定度量化

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

2.1 Hankel矩陣的擾動估計

(16)

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

(17)

(18)

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

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

提取其中一個相關函數矩陣塊,記為

(19)

(20)

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

(21)

(22)

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

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

2.2 模態參數的方差估計

首先,根據一階矩陣擾動傳播理論,由式(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為矩陣轉置。

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

(31)

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

(32)

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

(33)

考慮輸出矩陣C、系統矩陣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)

根據文獻[9],由第r階特征值擾動引起的第d階模態頻率、阻尼的擾動為

(38)

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

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

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

(39)

但是,通常結構模態振型為復數形式,需要進行歸一化處理。歸一化選擇對模態振型上的每一個點除以任意選擇的點φd,w進行歸一化,則模態振型的擾動則進一步轉化為[9]

(40)

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

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

(41)

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

(42)

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

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

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

(43)

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

3 桁架結構仿真算例驗證

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

圖1 桁架結構圖Fig.1 Truss structure

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

表2 模態阻尼比估計值的標準差Tab.2 Standard deviation of damping estimates

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

(44)

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

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

3.1 仿真結果

利用SSI-COV方法估計桁架結構模態參數,并根據本文所提方法計算模態參數標準差。將得到的模態參數標準差與MCS方法得到的模態參數標準差進行比較,從而驗證本文方法的有效性和可行性。

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

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

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

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

圖4 第一階模態頻率辨識結果Fig.4 First-order modal frequency identification results

圖5 第一階模態阻尼辨識結果Fig.5 First-order modal damping identification results

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

圖6 模態振型標準差Fig.6 Standard deviation of modal shape

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

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

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

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

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

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

3.2 不確定度量化的計算量分析

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

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

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

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

4 結 論

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

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