基于解析二階盲辨識的結構工作模態分析及真假模態區分

李澗鳴，包騰飛,3，周喜武，高瑾瑾，顧 昊

(1.河海大學水文水資源與水利工程科學國家重點實驗室，南京 210098；2.河海大學水利水電學院，南京 210098；3.三峽大學水利與環境學院，宜昌 443002；4.江蘇省水利工程科技咨詢股份有限公司，南京 210029；5.水利部交通運輸部國家能源局南京水利科學研究院，南京 210029)

模態參數識別是結構健康監測領域的關鍵問題之一。對于大型土木和水利工程結構，通常采用工作模態分析(Operational modal analysis, OMA)技術進行模態識別[1]。不同于試驗模態分析，OMA 無需施加人工激勵，僅需環境激勵下結構的振動響應識別參數，能更真實反映運行狀態下的結構特性。

為實現結構服役安全狀態的長期、實時和在線監測，需開發具備自動識別能力的OMA 方法。現有研究大多基于隨機子空間方法(Stochastic subspace identification, SSI)和穩定圖進行自動分析[2-3]。由于SSI 系統階次冗余，且結構在運行中易受到非平穩激勵和噪聲干擾，會不可避免地引入虛假模態，往往需要借助專家經驗才能從穩定圖中分辨出真實模態[4-6]。為了減少人為因素，實現自動識別，國內外學者致力于研究通過聚類、指標閾值和智能算法等方法自動區分虛假模態和真實模態。

MAGALHAES 等[7]利用譜系聚類對穩定極點分組，將超過一定數量極點的分組自動歸為真實模態。REYNDERS 等[8]基于譜系聚類和K 均值聚類提出了一種三階段聚類自動識別方法，并通過諸多模態指標來判斷真實模態。鄭沛娟等[9]引入圖論聚類通過二次聚類法識別結構真實模態。張永祥等[10]構造了2 個不同維度的分塊Hankel 矩陣并進行同階極點匹配得到更清晰的穩定圖，進而通過譜系聚類進行自動識別。祝青鑫等[11]綜合應用了主成分分析、K 均值聚類和譜系聚類提出了一種模態參數自動識別方法。雖然通過聚類等方法可以一定程度實現SSI 自動識別，但并不能根本解決穩定圖難以分辨的困難，尤其對于激勵相對較弱的結構，穩定圖中極點往往分布混亂，聚類結果很難保證準確可靠。

鑒于上述原因，本文另辟蹊徑，提出了一種基于解析二階盲辨識的模態參數識別方法，并基于模態指標和K 均值聚類提出了相應的真假模態區分方法，以期實現模態參數自動識別。該方法避免了引入穩定圖造成的模態選擇困難，能夠根據結構振動響應有效分離出各階模態響應，通過構造反映時、頻域模態響應特征的模態特征指標區分真實模態和虛假模態，進而利用K 均值聚類自動篩選出真實模態，并采用頻域參數擬合法估計模態參數。通過8 自由度數值算例和混凝土重力拱壩工程實例驗證本文方法的準確性和優越性。

1 基于解析二階盲辨識的模態參數識別

1.1 二階盲辨識算法

二階盲辨識(Second-order blind identification,SOBI)[12]是一種原本用于信號處理的盲源分離技術，PONCELET 等[13]成功將該技術引入結構模態分析中，并證明了其適用性。

盲源分離模型假定實測信號和未知源信號滿足線性混合關系，即對于n通道混合實測信號x(t)=[x1(t),x2(t),···,xn(t)]T，m組獨立的未知源信號s(t)=[s1(t),s2(t),···,sm(t)]T，滿足如下關系：

式中，A稱為混合矩陣。BSS 的目標是根據實測信號x(t)同時估計混合矩陣A和源信號s(t)，滿足：

式中，B為A的(廣義)逆矩陣，稱為分離矩陣。

SOBI 的基本思想是尋求使得實測信號的時延協方差矩陣近似對角化的源信號。定義時延協方差矩陣R(τ)=E{x(t)x(t-τ)T}，其中，τ為時間延遲系數，1≤τ≤p。采用近似聯合對角化[12]技術求得矩陣J，使得JR(τ)JT近似為對角矩陣。則問題為：

式中，off(·)為矩陣的非對角元。

實際計算中，需進行白化處理。引入線性變換：

式中，W為白化矩陣，可由主成分分析確定。對協方差矩陣R(0)=E{x(t)x(t)T}進行特征分解：

式中，E和D分別為R(0)的正交特征向量矩陣和特征值矩陣。則W計算如下：

容易證明，E{z(t)z(t)T}為單位矩陣。白化后的時延協方差矩陣為R(τ)=E{z(t)z(t-τ)T}，通過Jacobi算法可求出正交矩陣J，滿足式(3)。則混合矩陣A和源信號s(t)為：

對于N自由度線性系統，由模態疊加原理可知，結構振動響應y(t)與模態響應q(t)可以通過線性變換建立聯系，即：

式中， Ψ為模態向量組成的模態矩陣。

可以看出，式(9)與式(1)具有相似的形式，模態矩陣 Ψ對應于混合矩陣A，模態響應q(t)對應于源信號s(t)。因此，可以通過SOBI 算法得到Ψ和q(t)，進而識別模態參數。

1.2 解析二階盲辨識算法

由結構動力學可知，式(9)中只有當阻尼滿足經典阻尼假定(如比例阻尼)的情況下 Ψ才為N階實矩陣，故SOBI 理論上僅適用于阻尼簡化的實模態情形。在振動測試中往往需考慮更一般的復模態情形，此時Ψ ∈CN×2N，q(t)∈C2N×1。因此，有必要將SOBI 推廣以處理復模態問題。

由于模態向量與模態響應以共軛復數對的形式出現，則式(9)可表示為：

式中， (·)?為復共軛。

與p(t)實部pR(t)對應的解析信號表示為：

與式(11)類似，x(t)對應的解析信號為：

則由式(12)可得：

相應地，可將式(1)推廣為如下形式：

則1.1 節中的SOBI 可以自然推廣為復數形式的解析SOBI。即根據式(15)，由白化處理后的y?(t)得到時延協方差矩陣，再由近似聯合對角化求出酉矩陣Jc使得時延協方差矩陣近似對角化。類似式(7)和式(8)得到混合矩陣Ac與解析源信號s?(t)。相關過程與1.1 節完全類似，故不再贅述。

1.3 頻域參數擬合法估計模態參數

利用式(14)和式(15)的相似關系，可由解析SOBI 估計模態矩陣 Φ和與模態響應q(t)相關聯的解析信號p?R(t)。由矩陣 Φ得到結構振型，由可通過頻域參數擬合法估計固有頻率和阻尼比。

根據振動力學，N自由度系統的運動方程為：

式中：u(t)為位移向量；f(t)為外部激勵向量；M、C、K分別為質量矩陣、阻尼矩陣和剛度矩陣。系統的頻響函數H(iω)具有如下形式：

式中：ωnr第r階自振圓頻率； ζr為阻尼比。

在頻響函數的第r階共振峰附近，忽略共軛項，則式(17)可作以下近似：

假設激勵為白噪聲，則pR(t)的自功率譜密度函數P(ω)可作為H(ω)的估計。由式(19)可得：

P(ω)譜峰附近取l個頻率點ω1,ω2,···,ωl，則：

由式(21)可通過最小二乘法估計λr，則由式(18)

可得固有頻率fr和阻尼比 ζr為：

2 基于模態特征指標和K 均值聚類的真假模態自動區分

2.1 模態特征指標

由于結構被充分激勵的模態數目難以事先確定，在解析SOBI 分析中假定源信號的數目與實測信號維數相同，計算中會引入多余的虛假模態。另外，實測信號中還存在系統和測量噪聲、非平穩激勵以及濾波等引入的虛假模態。為實現模態參數自動識別，需要對結構真實模態和虛假模態進行自動區分，最后通過頻域參數擬合法識別真實模態參數。通過比較模態表現出的不同模態響應特性，分別從時域和頻域角度出發，提出了總包絡長度指標和譜熵指標兩個模態特征指標以區分真假模態。

2.1.1 總包絡長度指標

在時域中，真實模態的模態響應通常具有更短的包絡線[15]。根據解析信號理論，對于模態j，模態響應曲線qj(t)的包絡線可表示為。由勾股定理， dt時段內的包絡線長度近似為dl=則總包絡線長度為：

2.1.2 譜熵指標

在頻域中，真實模態的自功率譜密度函數通常更穩定，更有規律性。可以將其類比成概率密度函數。對于模態j，模態響應曲線qj(t)的自功率譜密度函數表示為Pj(ω)，為滿足概率積分為1 的條件，可對其進行歸一化：

在信息論中可通過信息熵度量信源的不確定性。則類似定義譜熵描述功率譜密度的不確定性：

2.2 K 均值聚類區分真假模態

對于解析SOBI 計算的各模態響應，分別計算總包絡長度和譜熵兩個模態特征指標，并歸一化到[0,1]區間以均衡權重。對于模態j，其總包絡長度和譜熵指標值構成二維數據點。本文采用聚類數為2 的K 均值聚類算法(K-means Clustering)自動將模態特征指標構成的二維數據點分成兩個簇，分別代表真實模態和虛假模態構成的集合。K 均值聚類算法的目標是尋求使得各簇中所有數據點與最近的聚類中心的距離之和最小的最優聚類方案，可表示如下：

式中：S1和S2分別為真實和虛假模態構成的簇；c1和c2分別為S1和S2對應的聚類中心點；//·//2為L2 范數，表征歐式距離。算法步驟如下：

步驟1：對于各數據點組成的集合{p1,p2,···,pm}，隨機選擇2 個初始聚類中心點c1、c2。

步驟2：分別計算各數據點pj到各聚類中心ck的歐式距離//pj-ck//2，尋找與各數據點距離最小的聚類中心，并將數據點分配到對應的簇中。

步驟3：采用平均的方法重新計算分配完成后的新簇的聚類中心，計算公式為：

式中，nk為Sk中數據點的個數。

通過K 均值聚類可自動提取真實模態，依次采用1.3 節所述的頻域參數擬合法即可估計出各真實模態對應的模態參數。

3 8 自由度質量-彈簧-阻尼模型數值算例

為驗證本文提出方法的有效性，并檢驗其自動區分真假模態的能力，采用MATLAB 平臺建立了一個8 自由度質量-彈簧-阻尼模型，其質量、剛度和阻尼矩陣分別為：

模型考慮非經典阻尼，阻尼矩陣可由滿足瑞利阻尼的部分添加非經典項形成。其中瑞利阻尼部分的第1 階和第8 階模態的阻尼比設定為2%。對模型的各自由度施加高斯白噪聲以模擬環境激勵。采用50 Hz 采樣頻率進行信號采集，各自由度采集60 min 時長的加速度響應。對各加速度信號均添加信噪比為10 dB 的高斯白噪聲以模擬真實測試環境中的測量噪聲。結合模型的固有頻率特點，對信號進行重采樣，奈奎斯特頻率為6.25 Hz。經過重采樣，上述加速度響應信號中僅包含模型前4 階模態信息。采用解析SOBI 進行模態識別時，假設源信號數等于信號通道數，因此會引入4 階虛假模態。如2.1 節所述，對解析SOBI 識別的8 個模態響應分別計算總包絡長度和譜熵指標，并采用K 均值聚類方法區分真假模態，如圖1 所示。可以看出，K 均值聚類算法成功將各模態分成兩個簇，其中4 階真實模態(“+”表示)構成的簇比4 階虛假模態(“×”表示)構成的簇的聚類中心更靠近坐標原點，可實現真假模態自動區分。

圖1 模態特征指標的K 均值聚類結果Fig.1 The result of K-means clustering for modal metrices

根據K 均值聚類結果可得4 階真實模態對應的模態矩陣及模態響應相關聯的解析信號。由頻域參數擬合法識別各階固有頻率和阻尼比，見表1。表1 中同時列出了相應理論值，并計算相對誤差。固有頻率和阻尼比的相對誤差絕對值分別小于0.04%和6%。作為對比，采用增強頻域分解法(Enhance frequency domain decomposition, EFDD)[16]對固有頻率和阻尼比進行識別，見表2。可以看出，本文方法對固有頻率和阻尼比的識別結果均具有較高精度，其中固有頻率的精度與EFDD 相近，而阻尼比的精度比EFDD 更高。

表1 本文方法識別的模態參數與理論值的比較Table 1 Comparisons between modal parameters identified by the proposed method and theoretical values

表2 EFDD 識別的模態參數與理論值的比較Table 2 Comparisons between modal parameters identified by EFDD and theoretical values

由模態矩陣得到結構振型，可通過模態置信準則(Modal assurance criterion, MAC)評價振型識別結果。第r階和第s階模態的模態向量vr和vs之間的MAC 值定義為：

式中：(·)H為共軛轉置；MAC 值越接近1，則振型越相關；反之越接近0，則越不相關。分別計算本文方法識別的各階振型與相應理論值之間的MAC值，見表1。相應地，EFDD 識別的各階振型與理論值之間MAC 見表2。可以看出，兩種方法識別的各階MAC 值均大于99.9%，具有較高精度。

4 龍羊峽大壩工程實例

4.1 工程概況

龍羊峽水利樞紐位于黃河上游，由混凝土重力拱壩主壩、兩岸重力墩、兩岸重力副壩、泄水建筑物、引水建筑物和水電站廠房等組成。擋水前緣總長度1226 m，其中主壩長396 m。最大壩高178 m，最大底寬80 m。建基標高2432 m，水庫正常蓄水位2600 m。工程處于高地震烈度區，基本烈度為VII 度，大壩按VIII 度地震設防。大壩于1986 年布置了強震監測系統，并于2009 年進行了技術改造，實現了震動實時監測能力。監測系統采用TDA-23 型三向加速度傳感器和TDA324 FA/CA 數字采集器，共布置13 個測點，采樣頻率為200 Hz。

4.2 大壩模態參數識別

采用主壩2600 m 廊道中布置的5 個加速度傳感器共15 個監測通道的監測數據進行大壩模態識別，如圖2 所示。采用2016 年6 月6 日11 點至12 點的60 min 加速度響應數據進行分析，如圖3 所示。

圖2 主壩2600 m 廊道加速度傳感器布置圖Fig.2 Positions of accelerometers in the 2600 m dam gallery

圖3 大壩60 min 加速度時程曲線Fig.3 Dam acceleration time series in 60 min

對信號進行重采樣，奈奎斯特頻率為10 Hz。采用本文提出的方法進行模態識別，假設源信號數等于信號通道數15。對解析SOBI 識別的15 個模態響應分別計算總包絡長度和譜熵指標，并采用K 均值聚類方法區分真假模態，如圖4 所示。可以看出，K 均值聚類可自動區分真假模態。根據K 均值聚類結果可得5 階真實模態對應的模態矩陣及模態響應相關聯的解析信號p?R(t)。由解析信號實部pR(t)計算得到各階模態響應的自相關函數和自功率譜密度函數，分別如圖5 和圖6 所示。為簡化起見，各模態響應自相關函數幅值歸一化到[-1, 1]區間。可以看出，各階模態時域和頻域均符合結構響應特征，解析SOBI 具有較好的模態分離效果。由模態矩陣得到大壩振型如圖7 所示，相應描述見表2。可以看出，第1 階、第3 階和第4 階振型為對稱振型，第2 階和第5 階振型為反對稱振型。

圖4 模態特征指標的K 均值聚類結果Fig.4 The result of K-means clustering for modal metrices

圖5 各階模態響應的自相關函數Fig.5 Autocorrelation functions of modal responses for each mode

圖6 各階模態響應的自功率譜密度函數Fig.6 Auto power spectral density functions of modal responses for each mode

圖7 本文方法的振型識別結果Fig.7 Mode shapes results identified by the proposed method

根據頻域參數擬合法，由pR(t)計算得到自功率譜密度函數P(ω)，各譜峰附近取7 個頻率點ω1,ω2,···,ω7，并相應得到P(ω1),P(ω2),···,P(ω7)，通過式(21)求得各階模態復極點 λr，見表3。最后由式(22)和式(23)得到各階固有頻率和阻尼比，見表4。作為對比，采用EFDD 對固有頻率和阻尼比進行識別，見表4。由于第3 階和第4 階模態固有頻率較接近，EFDD 難以識別第4 階模態，其余模態參數與本文方法識別的結果相近。計算本文方法和EFDD 識別的各階振型之間的MAC 值，見表4。可以看出，除了第3 階和第4 階密集模態外，其余模態振型的MAC 值均大于97%。可見，采用本文方法能有效識別結構的模態參數，尤其對傳統方法較難識別的密集模態具有更好的識別效果。

表3 頻域參數擬合法計算表Table 3 Calculation table of frequency domain parameter fitting method

表4 本文方法和EFDD 方法識別的模態參數Table 4 Modal parameters identified by the proposed method and EFDD

由SSI 方法得到穩定圖如圖8 所示。穩定圖中的極點尤其是第3 階和第4 階密集模態附近分布不清晰，較難分辨真實和虛假模態，更難以進行自動識別。而采用本文方法避免了引入穩定圖，各階模態尤其是密集模態分離效果較好，較易區分真假模態，識別過程清晰，更適合進行模態參數自動識別。

圖8 隨機子空間的穩定圖Fig.8 Stabilization diagram of SSI

5 結論

本文基于解析二階盲辨識提出了一種結構模態參數識別方法，并相應提出了真假模態自動區分方法。通過解析二階盲辨識分離模態響應，然后構造時、頻域模態特征指標表征模態特性，并利用K 均值聚類區分真假模態，采用頻域參數擬合法估計模態參數。通過數值算例和龍羊峽混凝土大壩工程實例驗證了方法的有效性，得到以下結論：

(1) 解析二階盲辨識具有較好的模態分離效果。分離出的各階模態響應符合結構響應特征。

(2) 本文提出的總包絡長度和譜熵兩個模態特征指標可以有效區分真實和虛假模態，通過K 均值聚類可以自動篩選出真實模態。

(3) 本文方法對模態參數的識別結果具有較高精度。數值算例中識別的固有頻率和阻尼比與理論值的相對誤差絕對值分別小于0.04%和6%。識別的各階振型與理論值之間的MAC 值均大于99.9%。本文方法識別的頻率和振型結果精度與EFDD 接近，而阻尼比結果精度比EFDD 更高。工程實例中識別的大壩固有頻率和阻尼比除第4 階模態外均與EFDD 方法識別的結果相近，除第3 階、第4階模態外，本文方法與EFDD 識別振型的MAC 值均大于97%。

(4) 本文方法對傳統方法較難識別的密集模態具有更好的識別效果。且可以克服由于引入穩定圖造成的極點分布不清晰，模態難以分辨的問題。更適合進行模態參數自動識別。