基于時變模態振型小波變換的結構損傷識別方法

王 盟，翁 順，余興勝，閆俊鋒，殷鵬程

(1.華中科技大學 土木工程與力學學院，武漢 430074；2.華中科技大學 控制結構湖北省重點實驗室，武漢 430074；3.中鐵第四勘察設計院集團有限公司，武漢 430063)

結構模態參數是大型土木工程結構優化設計、振動控制和損傷識別等研究領域的重要參數之一[1]。由于材料非線性、幾何非線性和狀態非線性等原因，土木工程結構本身是非線性的，其模態參數是時變的。結構損傷是一種非線性發展過程，也體現為時變的模態參數。結構損傷將導致結構模態參數的變化，如頻率、阻尼和振型。結構的振型包含了結構各個節點的信息，因此可以通過時變模態振型特征識別結構損傷位置和損傷程度。大多數傳統的模態參數識別方法都是基于時域或者頻域[2]，如快速傅里葉變換(fast Fourier transform，FFT)、時序分析法以及分區模態綜合法等。這些傳統方法雖然能識別出結構的模態參數[3]，但是也存在很多缺點：①頻域方法和時域方法都只能在頻域空間或時域空間中進行，而不能獲得信號在頻域空間和時域空間內的聯系，在實際應用中存在某些限制，例如，傅里葉變換只能得到信號中包含的頻率，但是卻不能得到這些頻率所對應的時域特征；②傳統方法大多數只適用于線性結構，然而實際上土木工程結構是非線性的，它的模態參數是時變的，傳統的識別方法難以展示時變特性。

目前時變模態參數識別方法主要有短時傅里葉變換(short time Fourier transform,STFT)、Hilbert-Huang變換(Hilbert-Huang transform,HHT)、Cohen類時頻分布、小波變換等[4]。STFT是使用一個固定的窗函數對信號進行傅里葉變換分析，通過移動窗函數來逐步分析信號，從而得到信號的時變特性[5]。劉小峰等[6]通過STFT方法，成功地對旋轉機械信號進行了精細時頻分析。STFT方法缺點在于其窗函數是固定的，因此它的分辨率是固定的，要改變分辨率必須更改窗函數。HHT方法主要包括經驗模態分解和Hilbert變換，這種方法的自適應性很強[7]。胡浩然等[8]結合奇異譜分析技術和HHT方法，準確地識別了故宮木結構古建筑的時變模態參數。HHT方法目前最大的問題是缺乏嚴謹完整的理論基礎。Cohen類時頻分布主要包括偽Wigner-Ville分布、平滑偽Wigner-Ville分布、Choi-Willams分布等[9]。郭奇等[10]結合二次聚合經驗模態分解和Wigner-Ville分布方法，改善了模態混疊問題。Cohen類時頻分布方法的缺點在于無法完全消除交叉項影響，從而造成虛假頻率成分產生。小波變換是一種窗口大小固定但形狀可變，時域和頻域分辨率都可改變的時頻局部化分析方法。它既克服了STFT分辨率固定的缺點，又能準確分析信號的真實頻率成分，并且有嚴謹完整的理論基礎，已經成為分析時變結構模態參數的重要方法[11]。

利用小波變換的多分辨率特性和時頻特性，可以得到瞬時頻率以及瞬時振型。Lardies[12]通過對結構信號進行連續小波變換，成功識別了線性時不變結構的頻率和振型。滕軍等[13]將復Morlet小波變換方法應用于大跨度空間結構，成功識別出了大跨度空間結構的頻率和振型。王超等[14]提出了一種基于動態規劃提取多成分信號小波脊和瞬時頻率的方法，通過找到小波系數模極大值，提取小波脊線，從而識別出時變結構的瞬時頻率。

將結構損傷當作一種非線性，它將引起結構模態振型的時變特征，因此通過時變模態振型可以識別結構的損傷狀況。閆桂榮等[15]總結了三大類結構損傷識別方法：基于模態域的方法、基于時間域的方法以及基于小波變換時頻域的方法。嚴平等[16]提出了一種基于單元模態應變能和小波變換的結構損傷識別方法，在單元模態應變能基礎上,利用小波變換系數的變化和分布情況構建單元模態應變能小波變換的結構損傷指標，準確識別出了結構的損傷位置以及損傷程度。目前，國內外學者已經研究了多種高精度、高效率的結構損傷識別方法，但是對于能同時識別結構損傷時間、位置以及程度的方法的研究依然匱乏。

基于時頻特征的結構損傷識別不僅可以識別結構損傷的位置和損傷程度，而且可以識別結構損傷發生的時間，特別適合于廣泛實施的實時健康監測系統。因此，本文提出了一種基于時變模態振型小波變換的結構損傷識別方法，利用小波變換先識別出結構的時變模態參數，然后對識別的時變振型差做小尺度的小波變換，通過小波系數的實部出現極大值的時間及其出現的位置識別出結構的損傷時間與損傷位置，通過極大值的大小判斷損傷程度，應用于懸臂梁結構以及昌贛客運專線贛江特大橋模型中，均成功識別出結構的損傷。

1 基本理論

小波變換是一種窗口大小固定但形狀可變，時域和頻域分辨率都可改變的時頻局部化分析方法。它可以對非穩態信號和非線性結構進行詳細的時頻分析和多分辨分析，已經成為分析時變結構模態參數的重要方法。

1.1 小波變換原理

設ψ(t)∈L2(R)(L2(R)為平方可積的實數空間，即能量有限的信號空間)為小波母函數，且ψ(t)滿足可容許條件

(1)

式中，ψ(ω)為小波母函數ψ(t)的傅里葉變換。將小波母函數通過伸縮和平移可以得到小波函數族

(2)

式中：a為尺度參數，可以控制小波母函數的伸縮；b為平移參數。

本文使用復Morlet小波作為小波母函數，其表達式為

(3)

式中：fb為帶寬參數；fc為小波的中心頻率，實際應用中，可以通過修改這兩個參數來達到精度的要求。復Morlet小波的頻域表達形式[17]為

(4)

對于任意信號x(t)，其連續小波變換可表示為

(5)

1.2 基于小波變換的瞬時模態參數識別方法

單自由度黏性阻尼系統的脈沖響應函數可以表示為

x(t)=Be-ζωntcos(ωdt+φ0)

(6)

xa(t)=x(t)+jH[x(t)]=Be-ζωntej(ωdt+φ0)

(7)

式中，H[x(t)]為x(t)的希爾伯特變換，將式(7)代入式(5)，并將Be-ζωnt在t=b附近泰勒展開得

(8)

(9)

多自由度黏性阻尼系統的脈沖響應函數可以表示為

(10)

同理，將式(10)代入式(5)，泰勒展開并忽略無窮小量可以得到

(11)

1.2.1 基于小波系數模極大值的瞬時頻率識別方法

對于單自由度黏性阻尼系統，設響應信號的采樣頻率為fs，代入式(6)可得

(12)

(13)

本文使用的是復Morlet小波，將復Morlet小波母函數代入式(13)，并對小波系數取模可得

(14)

由式(14)可知，在時刻k=b，當尺度參數a滿足aω′d=2πfc時，小波系數的模能取最大值。因此，在時刻k=b，只要找到能使小波系數模取最大值的尺度參數al，即可得到這一時刻的瞬時頻率為

(15)

因此，對于每個時刻，找到使小波系數模取最大值的尺度參數al，即可得到信號的瞬時頻率。

(16)

(17)

因此，在時刻k=b，只要找到能使小波系數模取第i個局部極大值的尺度參數ali，即可得到這一時刻的第i階瞬時頻率為

(18)

1.2.2 基于小波系數比值的瞬時振型識別方法

若一個結構有m個節點，假設：xp(t)為第p個節點的脈沖響應信號；xr(t)為第r個節點的脈沖響應信號。在尺度ali下對它們進行連續小波變換，由1.2.1節可知，當取aliω′di=2πfc時，多自由度系統就可解耦成只含第i階模態的系統。因此，Wψ,p(ali,b)和Wψ,r(ali,b)都只保留了各自節點的第i階模態振型。選取r節點為參照點，則在時刻t=bk，p節點相對于r節點的第i階振型χi,p,bk可表示為[18]

χi,p,bk=Wψ,p(ali,bk)/Wψ,r(ali,bk)

(19)

因此，對每個節點信號都做同樣的小波變換，然后進行歸一化處理，即可得到結構在時刻t=bk的第i階瞬時振型

Φi,bk=[χi,1,bk,χi,2,bk,…,χi,r-1,bk,1,χi,r+1,bk,…,χi,m-1,bk,χi,m,bk]

(20)

同理，其他任意時刻的瞬時振型都可以由式(19)和式(20)得到，則結構在任意時刻的第i階振型為

(21)

式中：w為時間步數；m為節點數。因此，結構任意時刻的振型都可由式(21)得到。

1.3 基于時變振型差小波變換的損傷識別方法

結構損傷是一種非線性發展過程，將導致結構模態參數的時變特征。結構的振型包含了結構各個節點的信息，因此本文從結構時變振型入手，利用小波變換的奇異性檢測能力和局部放大特性[19]，對結構第一階時變振型差做小尺度的小波變換，通過變換后的結果來識別損傷。

對于一個含有m個節點的結構，假設結構損傷前的第一階理論振型為

Φ1,0=[χ1,1,0,χ1,2,0,…,χ1,r-1,0,1,χ1,r+1,0,…,χ1,m,0]

(22)

結構在任意時刻的振型可由式(19)以及式(20)計算得到，則結構在時刻t=bk的第一階振型差可以表示為

ΔΦ1,bk=[Δχ1,1,bk,Δχ1,2,bk,…,Δχ1,r-1,bk,0,Δχ1,r+1,bk,…,Δχ1,m,bk],

(23)

當結構發生損傷時，損傷位置的振型變化會大于其他位置的變化[20]。因此，利用本文復Morlet小波實部的奇異性檢測能力和局部放大特性，對結構的振型差做小尺度小波變換

(24)

若結構在這個時刻(t=bk)有損傷，則小波系數實部出現極大值的位置即為結構損傷位置。

同理，任意時刻的第一階振型差都可由式(23)得到，然后利用式(24)對振型差做小尺度小波變換，則小波系數的實部出現極大值的時間bk即為結構的損傷時間，出現極大值的位置即為結構損傷位置，極大值的大小反映了損傷的程度。

2 數值仿真

2.1 懸臂梁算例

通過一個1 m長懸臂梁數值模型，驗證時變模態參數小波變換方法用于損傷識別的有效性。懸臂梁由10.0個單元組成，其模型示意圖，如圖1所示。

圖1 懸臂梁結構模型Fig.1 The model of cantilever beam structure

結構在未損傷前，各個單元的抗彎剛度EI=4.73×103N·m2，抗拉剛度EA=1.42×106N，密度ρ=2 500 kg/m3。為了驗證本文方法的精確性，在懸臂梁上施加了Imperial Valley地震波，地震波采樣頻率為200 Hz，如圖2所示。

圖2 Imperial Valley地震波加速度曲線Fig.2 Acceleration curve of seismic wave in Imperial Valley

本文通過剛度改變模擬結構非線性損傷，有利于定量損傷程度，通過減小單元彈性模量E的值模擬單元剛度損失。當結構施加如圖2的地震波時，假設：當t=5 s時，結構的單元2發生了10%損傷；當t=10 s時，結構的單元3發生了5%損傷。因此，整個損傷過程可以分為3個階段，如表1所示。

表1 懸臂梁損傷狀況表Tab.1 Damage stage table of the cantilever beam

對于每個階段，運用結構動力學知識，由彈性模量以及密度等參數計算結構的剛度矩陣和質量矩陣，從而可以得到結構每個階段的第一階理論頻率為：f1=3.849 9 Hz，f2=3.797 0 Hz，f3=3.780 2 Hz。將節點11作為參照節點r，還可以得到每個階段結構的第一階理論振型，如圖3所示。由圖3可知：雖然結構在不同階段有不同的損傷，但理論振型變化很細微，僅憑結構在3個階段的振型變化難以識別結構損傷。

圖3 懸臂梁第一階理論振型Fig.3 The first theoretical mode shape of the cantilever beam

施加地震波后，結構每個節點的動力響應信號可以由Newmark法模擬得到。本文選取位移響應信號分析，節點11前20 s的原始位移響應，如圖4所示。

圖4 節點11前20 s的位移信號Fig.4 The displacement signal for the first 20 s of node 11

2.1.1 瞬時頻率與瞬時振型的識別

考慮到模擬損傷程度的大小和實際工程中噪聲的影響，本文在原始位移響應中添加了1%的高斯白噪聲，分別探究無噪聲和有噪聲環境下本文算法的有效性。使用復Morlet小波作為母小波，設置小波中心頻率fc=3 Hz。將兩種情況下節點11的位移信號作為x(t)代入連續小波變換。由1.2.1節可知，只要找到能使小波系數模取第i個局部極大值的尺度參數ali，即可得到這一時刻的第i階瞬時頻率。因此，通過提取每個時刻滿足條件的尺度參數al1，代入式(18)，即可得到結構的第一階瞬時頻率，如圖5所示。

由圖5可知，無噪聲識別結果比有噪聲識別結果得到的曲線更加光滑。兩種情況下，頻率在t=5 s以及t=10 s時都發生了突變，說明結構在這兩個時刻發生了損傷。時變頻率分為3個階段，0～5 s對應表1的階段1，5～10 s對應表1的階段2，10～20 s對應表1的階段3。將每個階段頻率的平均值作為識別值，與理論值的對比如表2所示。

圖5 懸臂梁第一階瞬時頻率Fig.5 First order instantaneous frequency of the cantilever beam

表2 懸臂梁頻率對比表Tab.2 Frequency comparison table of the cantilever beam Hz

由表2可知，兩種情況下，各個階段第一階頻率的識別值與理論值很相近。進一步地，通過計算發現，兩種情況下，各個階段識別值與理論值的相對誤差絕對值都在2%以下，說明本文小波變換方法能高精度地提取結構的瞬時頻率，且具有一定的抗噪性。

在無噪聲和有噪聲這兩種情況下，分別對每個節點的位移信號做同樣的小波變換，并選取節點11為參照節點r，運用式(19)和式(20)即可計算結構任意時刻的第一階瞬時振型，如圖6、圖7所示。

圖6 無噪聲下懸臂梁第一階瞬時振型Fig.6 First order instantaneous mode shape of the cantilever beam without noise

由圖6、圖7可知，瞬時振型曲線與理論振型相似，并且通過計算，發現每個階段識別的瞬時振型與對應階段的理論振型的模態置信度(modal assurance criterion,MAC)值接近1，這說明通過小波變換方法能高精度地提取結構的瞬時振型。

圖7 有噪聲下懸臂梁第一階瞬時振型Fig.7 First order instantaneous mode shape of the cantilever beam under noise

2.1.2 損傷識別

圖8 無噪聲下第一階段損傷識別Fig.8 Damage identification of the first stage without noise

圖9 無噪聲下第二階段損傷識別Fig.9 Damage identification of the second stage without noise

圖10 無噪聲下第三階段損傷識別Fig.10 Damage identification of the third stage without noise

圖11 無噪聲下全過程損傷識別Fig.11 Damage identification of the whole stage without noise

圖12 有噪聲下全過程損傷識別Fig.12 Damage identification of the whole stage under noise

由圖8可知，沒有極大值出現，說明在階段1結構無損傷；由圖9可知，在t=5 s之后，結構在節點3處(對應單元2)有極大值，說明在階段2結構的單元2發生了損傷；由圖10可知，在t=10 s之后，結構在節點3(對應單元2)和節點4處(對應單元3)均出現了極大值，說明階段3結構的單元2、單元3都發生了損傷；由圖11和圖12可知，有噪聲和無噪聲環境下均能準確識別出損傷位置；由圖8～圖12可知，結構在t=5 s時，單元2發生損傷，在t=10 s時，單元2、單元3都發生損傷，識別的結果與表1中設置的剛度損失工況一樣，說明本文方法能夠精確的識別出結構的損傷時間和損傷位置，且具有一定的抗噪性。

由表1可知，在階段2結構發生了單損傷，在階段3結構發生了多損傷。在階段2時，為了進一步判斷損傷程度，依次設置單元2的損傷程度為10%，15%，20%，…，30%，用本文的方法依次提取每種損傷程度下振型差小波系數的實部相對極大值，如圖13所示。

圖13 不同損傷程度下單元2小波系數實部相對極大值Fig.13 The relative maximum value of the real part of the wavelet coefficients of unit 2 in different damage degree

由圖13可知，實部相對極大值隨著損傷程度的增大而增大，即可判斷：單損傷下，振型差小波系數的實部相對極大值越大，單元的損傷程度越大。

階段3結構發生了多損傷，為了進一步判斷單元3的損傷程度，保持單元2的損傷不變(10%)，依次設置單元3的損傷程度為5%，10%，15%，…，30%，用本文的方法依次提取每種損傷程度下單元3振型差小波系數的實部相對極大值，如圖14所示。

由圖14可知，實部相對極大值隨著損傷程度的增大而增大，即可判斷：多損傷下，在其他單元不變的情況下，所求單元振型差小波系數實部相對極大值越大，該單元的損傷程度越大。

圖14 不同損傷程度下單元3小波系數實部相對極大值Fig.14 The relative maximum value of the real part of the wavelet coefficients of unit 3 in different damage degree

綜上所述，本文方法能夠精確地識別出懸臂梁結構的損傷時間、損傷位置和損傷程度，且具有一定的抗噪性。

2.2 昌贛客運專線贛江特大橋算例

昌贛客運專線贛江特大橋(以下簡稱贛江橋)位于江西省贛州市，其主橋橋位位于章水、貢江兩江匯合口下游1.9 km處，位于既有贛江公路大橋上游1.1 km處。主橋結構采用(35+40+60+300+60+40+35)m混合梁斜拉橋，半漂浮體系，全梁長572.1 m(含梁縫)。索塔全高分別為124.5 m，127.6 m，橋面以上塔高88 m，為主跨的1/3.409。斜拉索采用空間雙索面體系，扇形布置，全橋共48對斜拉索。本文通過贛江橋實測響應數據不斷進行有限元模型修正，得到了實際橋梁的等效模型，其有限元模型如圖15所示。

圖15 贛江橋模型Fig.15 The model of Ganjiang River Bridge

橋梁的主跨由108個節點(節點81～節點188)以及107個單元(單元79～單元185)組成。根據美國國家公路與運輸協會標準(AASHTO)對于船舶碰撞荷載提出的經驗公式，船橋正撞時的設計船舶撞擊力按式(25)計算

(25)

式中：P為撞擊力；DWT為載質量噸位；v為船舶的速度。現實中，船橋碰撞時，相撞船舶周圍水的質量對荷載也有影響，可以將這部分影響以船體附加質量的形式加以修正，即

DWT=m1+m2=(1+k)m1

(26)

式中：m1為船舶載質量噸位；m2為附加水質量；k為附加因子，一般取0.02～0.07，本文取0.05。

利用ANSYS軟件模擬速度為6 m/s，質量為850 t的船撞向贛江橋的節點123(對應單元120)，由式(25)、式(26)可得沖擊力為22 MN，使用ANSYS軟件完成了船橋碰撞瞬態分析，并通過有限元模型修正方法獲取了單元等效剛度損失為30%。為了便于定量損傷程度，本文通過設定碰撞后單元彈性模量折減30%證明本文方法的精度。通過設置標準差為0.01的白噪聲激勵信號模擬環境激勵，當t=5 s時，通過22 MN荷載模擬碰撞沖擊力，采樣頻率為50 Hz，由Newmark法模擬得到各個節點的動力響應，節點123碰撞前后10 s的位移響應信號，如圖16所示。

圖16 節點123的位移信號Fig.16 The displacement signal of node 123

整個船橋碰撞過程分為兩個階段，如表3所示。當t<5 s時，結構為未損狀態；當t≥5 s時，船橋碰撞，碰撞處發生了損傷。通過數值模型計算可以得到橋梁每個階段的第一階理論頻率為：f1=0.414 1 Hz，f2=0.412 9 Hz。將節點132作為參照節點r(振型為“1”的節點)，進行歸一化處理，還可以得到橋梁碰撞前后第一階理論振型，如圖17所示。

表3 “船撞橋”損傷狀況表Tab.3 Damage stage table of “ship-bridge collision”

圖17 橋梁碰撞前后第一階理論振型Fig.17 The first theoretical mode shape before and after collision

由圖17可知，雖然橋梁在碰撞后有損傷，但碰撞前后理論振型變化很細微，僅憑橋梁在碰撞前后振型變化難以識別橋梁損傷。

2.2.1 瞬時頻率與瞬時振型的識別

使用復Morlet小波作為母小波，設置小波中心頻率fc=5 Hz。將節點123的位移信號作為x(t)代入連續小波變換。找到每個時刻滿足極大值條件的尺度參數al1，代入式(18)，即可得到結構的第一階瞬時頻率，如圖18所示。

圖18 贛江橋第一階瞬時頻率Fig.18 First order instantaneous frequency of Ganjiang River Bridge

由圖18可知，頻率在t=5 s時發生了突變，說明結構在這個時刻發生了損傷。時變頻率分為兩個階段，0～5 s對應表3的階段1，5～10 s對應表3的階段2，將每個階段頻率的平均值作為識別值，與理論值的對比如表4所示。

表4 贛江橋頻率對比表Tab.4 Frequency comparison table of Ganjiang River Bridge

由表4可知，第一階頻率的識別值與理論值相對誤差絕對值非常小，各個階段的誤差都在2%以下。充分說明了小波變換方法識別瞬時頻率的精確性。

對每個節點的位移信號做同樣的小波變換，并選取節點132為參照節點r(振型為“1”的節點)，運用式(19)和式(20)可計算贛江橋任意時刻的第一階瞬時振型，如圖19所示。由圖19可知，瞬時振型曲線與理論振型相似，并且通過計算，發現每個階段識別的瞬時振型與對應階段的理論振型的MAC值接近于1，這說明通過本文小波變換方法能高精度地提取贛江橋的瞬時振型。

圖19 贛江橋第一階瞬時振型Fig.19 First order instantaneous mode shape of Ganjiang River Bridge

2.2.2 損傷識別

圖20 碰撞前損傷識別Fig.20 Damage identification before collision

圖21 碰撞后損傷識別Fig.21 Damage identification after collision

圖22 全過程損傷識別Fig.22 Damage identification of the whole stage

由圖20可知，當t=0～5 s時，沒有極大值出現，說明在這一階段橋梁無損傷；由圖21可知，當t=5～10 s時，橋梁在節點123處(對應單元120)有極大值，說明在這一階段橋梁的單元120發生了損傷。因此，可以判斷當t=5 s時，發生了“船撞橋”事件，且撞擊的位置為單元120，識別的結果與損傷狀況吻合，說明本文方法能夠精確地識別出船橋碰撞的時間和位置。

3 結 論

(1)本文提出了一種基于時變模態振型小波變換的結構損傷識別方法。通過提取當小波系數模取極大值時的尺度參數，識別出結構的瞬時頻率。然后選取特定尺度參數，對結構各個節點的信號進行連續小波變換，從而將信號解耦到只含某階模態，通過小波系數的比值進行歸一化處理，從而得到結構的瞬時振型。最后，對結構第一階時變振型差做小尺度的小波變換，通過小波系數的實部出現極大值的時間及其出現的位置識別出結構的損傷時間與損傷位置，通過極大值的大小判斷損傷程度。

(2)將該方法應用于一個地震波作用下的懸臂梁結構中，結果表明，該方法不僅能準確地識別出結構單損傷、多損傷發生的時間，還能準確地識別出損傷的位置以及損傷程度。

(3)將該方法應用于一個沖擊荷載作用下的“船撞橋”事件模擬。結果表明，該方法不僅能夠精確地識別出船橋碰撞的時間，并且能夠準確地識別出碰撞位置。