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

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

模態(tài)參數(shù)識別是結(jié)構(gòu)健康監(jiān)測領(lǐng)域的關(guān)鍵問題之一。對于大型土木和水利工程結(jié)構(gòu)，通常采用工作模態(tài)分析(Operational modal analysis, OMA)技術(shù)進行模態(tài)識別[1]。不同于試驗?zāi)B(tài)分析，OMA 無需施加人工激勵，僅需環(huán)境激勵下結(jié)構(gòu)的振動響應(yīng)識別參數(shù)，能更真實反映運行狀態(tài)下的結(jié)構(gòu)特性。

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

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

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

1 基于解析二階盲辨識的模態(tài)參數(shù)識別

1.1 二階盲辨識算法

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

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

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

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

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

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

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

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

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

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

對于N自由度線性系統(tǒng)，由模態(tài)疊加原理可知，結(jié)構(gòu)振動響應(yīng)y(t)與模態(tài)響應(yīng)q(t)可以通過線性變換建立聯(lián)系，即：

式中， Ψ為模態(tài)向量組成的模態(tài)矩陣。

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

1.2 解析二階盲辨識算法

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

由于模態(tài)向量與模態(tài)響應(yīng)以共軛復數(shù)對的形式出現(xiàn)，則式(9)可表示為：

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

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

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

則由式(12)可得：

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

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

1.3 頻域參數(shù)擬合法估計模態(tài)參數(shù)

利用式(14)和式(15)的相似關(guān)系，可由解析SOBI 估計模態(tài)矩陣 Φ和與模態(tài)響應(yīng)q(t)相關(guān)聯(lián)的解析信號p?R(t)。由矩陣 Φ得到結(jié)構(gòu)振型，由可通過頻域參數(shù)擬合法估計固有頻率和阻尼比。

根據(jù)振動力學，N自由度系統(tǒng)的運動方程為：

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

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

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

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

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

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

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

2 基于模態(tài)特征指標和K 均值聚類的真假模態(tài)自動區(qū)分

2.1 模態(tài)特征指標

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

2.1.1 總包絡(luò)長度指標

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

2.1.2 譜熵指標

在頻域中，真實模態(tài)的自功率譜密度函數(shù)通常更穩(wěn)定，更有規(guī)律性?？梢詫⑵漕惐瘸筛怕拭芏群瘮?shù)。對于模態(tài)j，模態(tài)響應(yīng)曲線qj(t)的自功率譜密度函數(shù)表示為Pj(ω)，為滿足概率積分為1 的條件，可對其進行歸一化：

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

2.2 K 均值聚類區(qū)分真假模態(tài)

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

式中：S1和S2分別為真實和虛假模態(tài)構(gòu)成的簇；c1和c2分別為S1和S2對應(yīng)的聚類中心點；//·//2為L2 范數(shù)，表征歐式距離。算法步驟如下：

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

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

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

式中，nk為Sk中數(shù)據(jù)點的個數(shù)。

通過K 均值聚類可自動提取真實模態(tài)，依次采用1.3 節(jié)所述的頻域參數(shù)擬合法即可估計出各真實模態(tài)對應(yīng)的模態(tài)參數(shù)。

3 8 自由度質(zhì)量-彈簧-阻尼模型數(shù)值算例

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

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

圖1 模態(tài)特征指標的K 均值聚類結(jié)果Fig.1 The result of K-means clustering for modal metrices

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

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

表2 EFDD 識別的模態(tài)參數(shù)與理論值的比較Table 2 Comparisons between modal parameters identified by EFDD and theoretical values

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

式中：(·)H為共軛轉(zhuǎn)置；MAC 值越接近1，則振型越相關(guān)；反之越接近0，則越不相關(guān)。分別計算本文方法識別的各階振型與相應(yīng)理論值之間的MAC值，見表1。相應(yīng)地，EFDD 識別的各階振型與理論值之間MAC 見表2?？梢钥闯?，兩種方法識別的各階MAC 值均大于99.9%，具有較高精度。

4 龍羊峽大壩工程實例

4.1 工程概況

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

4.2 大壩模態(tài)參數(shù)識別

采用主壩2600 m 廊道中布置的5 個加速度傳感器共15 個監(jiān)測通道的監(jiān)測數(shù)據(jù)進行大壩模態(tài)識別，如圖2 所示。采用2016 年6 月6 日11 點至12 點的60 min 加速度響應(yīng)數(shù)據(jù)進行分析，如圖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。采用本文提出的方法進行模態(tài)識別，假設(shè)源信號數(shù)等于信號通道數(shù)15。對解析SOBI 識別的15 個模態(tài)響應(yīng)分別計算總包絡(luò)長度和譜熵指標，并采用K 均值聚類方法區(qū)分真假模態(tài)，如圖4 所示?？梢钥闯?，K 均值聚類可自動區(qū)分真假模態(tài)。根據(jù)K 均值聚類結(jié)果可得5 階真實模態(tài)對應(yīng)的模態(tài)矩陣及模態(tài)響應(yīng)相關(guān)聯(lián)的解析信號p?R(t)。由解析信號實部pR(t)計算得到各階模態(tài)響應(yīng)的自相關(guān)函數(shù)和自功率譜密度函數(shù)，分別如圖5 和圖6 所示。為簡化起見，各模態(tài)響應(yīng)自相關(guān)函數(shù)幅值歸一化到[-1, 1]區(qū)間?？梢钥闯?，各階模態(tài)時域和頻域均符合結(jié)構(gòu)響應(yīng)特征，解析SOBI 具有較好的模態(tài)分離效果。由模態(tài)矩陣得到大壩振型如圖7 所示，相應(yīng)描述見表2?？梢钥闯?，第1 階、第3 階和第4 階振型為對稱振型，第2 階和第5 階振型為反對稱振型。

圖4 模態(tài)特征指標的K 均值聚類結(jié)果Fig.4 The result of K-means clustering for modal metrices

圖5 各階模態(tài)響應(yīng)的自相關(guān)函數(shù)Fig.5 Autocorrelation functions of modal responses for each mode

圖6 各階模態(tài)響應(yīng)的自功率譜密度函數(shù)Fig.6 Auto power spectral density functions of modal responses for each mode

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

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

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

表4 本文方法和EFDD 方法識別的模態(tài)參數(shù)Table 4 Modal parameters identified by the proposed method and EFDD

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

圖8 隨機子空間的穩(wěn)定圖Fig.8 Stabilization diagram of SSI

5 結(jié)論

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

(1) 解析二階盲辨識具有較好的模態(tài)分離效果。分離出的各階模態(tài)響應(yīng)符合結(jié)構(gòu)響應(yīng)特征。

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

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

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