基于小波包能量曲率差的古木結構損傷識別

王 鑫， 胡衛兵, 孟昭博

(1.西安建筑科技大學 土木工程學院, 西安 710055; 2聊城大學 建筑工程學院, 山東 聊城 252000)

我國是文明古國，古建筑是重要的歷史文化遺產組成部分，是中華民族的瑰寶。西安鐘樓位于東、西、南、北四條大街交匯處，是西安的標志性建筑，是國家重點保護文物。上世紀80年代地面交通引起的微幅振動使鐘樓二層地板上的展品發生移位，高臺基的墻體和木構架的橫梁裂縫不斷擴展。雖然鐘樓環道內側至高臺基邊緣的距離擴大到26.7 m，車輛引起的振動對鐘樓的影響有增無減，鐘樓二層地板有明顯的振感。西安是西部發展重心，現代工業迅速發展，人口不斷增加，地面交通流量不斷增大，為從根本上緩解地面交通壓力，西安市規劃建設6條地鐵線路，總長達251.8 km，鐘樓四側有四條隧道穿過，距離鐘樓僅15 m。對古木結構來說，交通環境振動是歷經數百年古建木結構壽命新威脅，也是正常情況下經受的最嚴重的威脅，長期的交通環境微振動能使木結構疲勞、榫卯松動、結構變形、壽命縮短。鐘樓一旦遭到破壞，不可再生，造成歷史文化遺產不可彌補的損失，這就要求結構健康監測與安全評估系統及時發現損傷并在第一時間預警。因此利用環境激勵下的結構振動響應對結構進行健康監測，已是結構健康監測中損傷檢測的新研究方向，同時對保護古建筑有著積極意義。

通常的損傷識別系統分為四個階段: ① 損傷是否存在; ② 損傷位置判定; ③ 損傷的嚴重程度; ④ 結構的實用性預測。損傷識別方法分為局部和整體兩種方法。目前有直觀判定、超聲波、電磁場及電渦流等檢測方法，需預先知道損傷的大體位置,對損傷部位進行儀器操作，整體損傷識別通過結構振動特性的變化評價結構的健康情況,只有將兩者結合起來,才能準確評價復雜結構的健康狀況[1-2]。傳統的模態參數頻率和振型不能準確識別損傷,大多數損傷評估方法通過對瞬態信號的傅里葉變換得到模態參數，而傅里葉變換的最大缺陷是對高頻模態分析不足。小波分析被譽為數學顯微鏡,本身具有放大、縮小和平移等功能,可通過檢查不同放大倍數下的變化來研究信號的特征,具有優良的時頻局部化特性，但有高頻段分辨率差的缺點。小波包分析是小波變換的擴充,能為信號提供一種更加精細的分析方法, 將小波分析沒有細分的高頻部分進一步分解, 具有任意的時-頻分辨率，根據分析信號的特征選擇相應的頻帶,使之與信號頻譜匹配。因此小波包分析在土木工程結構的健康監測與損傷診斷中具有非常廣闊的應用前景。

20世紀80年代開始，國內外學者開展了小波包分析在結構損傷診斷領域中的應用研究，取得了一定的研究成果。丁幼亮等[3]對Benchmark 鋼框架結構試驗數據和潤揚大橋懸索橋監測數據進行小波包能量譜損傷預警分析，在此基礎上詳細考察了不同小波函數和小波包分解層次的損傷預警效果。劉濤等[4]做了基于小波包能量譜的結構損傷預警方法試驗研究。鄧揚等[5]采用小波包分析進行了拉索損傷聲發射信號特征提取。范穎芳等[6]利用小波包對不同損傷情況下拱橋結構的觀測信號分析，確定了結構異常狀態的敏感特征。韓建剛等[7]提出了小波包變換的能量變化率指標對梁體進行了損傷識別的定位研究，并進行試驗驗證。余竹等[8]利用滄州子牙河新橋替換下的梁體進行兩種工況損傷模擬，用小波包能量曲率差法識別損傷，考察小波函數和分解層數對識別效果的影響。

本文以西安鐘樓為工程依托，提出小波包能量曲率差對古木結構進行損傷定位識別，把有限元分析得到梁上各節點的加速度響應信號進行小波包分解，計算該指標進行損傷定位，該指標對于損傷識別比較敏感,能準確判定古木結構損傷的具體位置。

1 小波包分析

小波包由一系列線性組合小波函數組成：

(1)

式中，i，j，k分別表示頻率因子、尺度因子和平移因子。

小波函數ψi遞推關系式為:

(2)

(3)

式中:ψ表示小波母函數，h(k)和g(k)為與尺度函數及小波母函數相關的積分鏡像濾波器系數。

對于任意信號的第j階和第j+1階水平小波包分解遞推關系為:

(4)

(5)

(6)

其中，H和G分別為h(k)和g(k)構成的濾波算子，

(7)

(8)

經過j水平的小波包分解后，初始信號f(t)為:

(9)

(10)

小波包系數為:

(11)

小波包系數滿足正交條件:

(12)

小波包信號能量為:

(13)

將式(10)代入式(13)，并利用式(12)得到：

(14)

(15)

小波包組分能量對信號的變化十分敏感，可用于結構的損傷識別。

2 古木結構的損傷識別

本文把古木結構的損傷識別分為兩部分:① 對古木結構梁上各節點的加速度響應信號進行小波包分解；② 計算小波包能量曲率差進行結構的損傷定位識別，其中包括小波函數和小波包分解層次的選擇。

2.1 合理的選擇小波函數

從消失矩和支撐長度考慮，選用Daubechies為小波函數，簡記為dbN(N為階次)。N越大，Daubechies小波的消失矩越高，時域的分辨率越好；但同時Daubechies小波支撐長度越寬，小波的時域局域性越差。因此應合理確定Daubechies小波階次N。

對結構動力響應f(N,k)進行第i層小波包分解,fij表示第i層分解節點(i,j)的結構響應,每個頻帶內結構響應fij能量[9]:

(16)

則結構動力響應f(N,k)第i分解層的小波包能量譜向量Ei:

(17)

為衡量小波函數好壞，定義i分解層各頻帶能量系數系列{Eij}的代價函數M{Eij}。小波包能量譜中各頻帶能量系數Eij的時頻集中程度由代價函數M{Eij}反映。采用lp范數熵為代價函數，在同一小波包分解層上，計算不同小波函數的代價函數值并比較，確定較適合的Daubechies小波階次N。通常不同的階次計算小波函數的代價函數值越小越好,lp范數熵(1≤p≤2)定義為[9]：

(18)

2.2 合理的選擇小波包分解層數

在工程應用中，對結構動力響應進行小波包分解，計算每一分解層次上的小波包能量譜的代價函數，從代價函數和計算時間考慮確定適當的小波包分解層次，通常小波包能量譜代價函數值越小，計算機計算過程耗時越少，小波包分解層次越好。類似小波函數階次的選擇方法，采用lp范數熵為代價函數，確定合適的小波包分解層數，lp范數熵(1≤p≤2)定義[9]：

(19)

2.3 不等間距曲率求解方法[8]

在實際工程中，曲率一般由變量的二階差分(斜率的變化率)得到。

不等間距情況下曲率求解式為:

(20)

式中:分子為節點左右兩段曲線斜率差，分母為節點左右兩端斜率差間距。若節點等間距，hi-1=hi+1，則:

式(21)為二階差分法求解等間距曲率公式。 將式(21)中y換成小波包能量譜，則為小波包能量曲率。將完好狀態與損傷狀態各節點的小波包能量曲率進行插值，得到損傷狀態的小波包能量曲率差為:

(22)

3 算例

3.1 古木結構的有限元模擬

本文以西安鐘樓為工程依托，選取其中一榀框架進行分析。通過Ansys有限元軟件對環境激勵下古木結構進行損傷模擬分析，以鐘樓為參考，選取木框架計算參數，選取木梁長4 m,木柱高6 m,梁截面尺寸為300×700 mm2，柱截面直徑500 mm，用beam188梁單元模擬木柱、木梁，用combin14單元模擬榫卯節點，榫卯連接的彎曲剛度為[10]：1×1010kN·m/rad，木材的彈性模量取1×1010N/m2，泊松比為0.25，密度為410 kg/m3，采用Rayleigh定義的粘性比例阻尼。柱子擱置在有凹槽的柱礎上，不能完全限制柱子轉動，在荷載引起的微幅振動下不計其線位移，柱與基礎的連接簡化成固定鉸支座的力學模型符合實際情況[11-12]，建立古木結構的有限元模型如圖1所示。

圖1 古木結構的有限元模型(a)和節點詳圖(b)

西安鐘樓處于地面交通和地鐵運行的復雜交通環境下，地面交通振動通過高臺基傳播引起鐘樓振動，運行的地鐵對軌道產生的沖擊作用產生振動，通過隧道結構傳到周圍地層，并經過地層向周圍傳播激勵鐘樓產生振動，因而在該木框架的柱底節點1處沿x軸正方向施加隨機激勵荷載來模擬環境激勵對鐘樓的影響[13]，隨機激勵荷載的時程曲線及頻譜如圖2所示，獲得結構的加速度荷載時程，在此基礎上運用Matlab程序計算了小波包能量譜。

古木結構的損傷程度通過折減損傷單元的彈性模量來實現，其中10%、18%、20%分別指損傷單元的彈性模量減少10%、18%、20%[9]，如表1所示。

對損傷工況1、2進行分析，得出完好結構和損傷工況1、2梁跨中第31節點的豎向加速度時程曲線如圖3所示。從圖3看出各損傷工況的信號有細微差別，但很難判斷古木結構的損傷情況。因此下面采用小波包能量曲率差對古木結構進行損傷識別。

圖2 激勵荷載的時程曲線(a)和頻譜曲線(b)

表1 古木結構損傷工況

3.2 選擇計算參數

3.2.1 小波函數的選擇

梁跨中撓度較大，是易出現損傷部位，因而對完好結構梁跨中第31節點的豎向加速度響應，選擇不同階次的Daubechies進行小波包分解，分解層次取4，計算lp范數熵代價函數值如表2所示。從表2看出，當小波階次為20時，lp范數熵為5 592.46，其值相對其他小波階次最小，因此損傷識別小波函數選擇Daubechies20。

圖3 第31節點加速度響應

表2 分解層次為4時不同db小波的代價函數值

3.2.2 小波包分解層數的選擇

采用Daubechies20對梁跨中第31節點完好狀態下的豎向加速度響應進行小波包分解，分解層次取1～8，計算lp范數熵的代價函數值，并記錄計算機計算耗費的時間如表3所示。從表3看出，當小波包分解層次為4時代價函數值為5 592.46,計算機計算時間為0.109 s，代價函數和計算時間均相對較小，因此損傷識別小波包分解層次取4。

表3 Daubechies 20 不同分解層次的代價函數值和計算時間

注:計算機的CPU為Intel (R) Core(TM) i5 M2430 2.40 GHz。

3.3 損傷識別

3.3.1 自振頻率

為了研究古木結構不同損傷程度的損傷識別指標,列出完好結構、損傷工況1、損傷工況2的自振頻率f0、f1、f2見表4所示。

從表4可以看出古木結構的損傷對自振頻率的影響非常小,損傷工況1與完好結構的自振頻率最大誤差僅為0.63%，損傷工況2與完好結構的自振頻率最大誤差僅為1.39%，看來利用結構自振頻率的變化來發現古木結構的損傷十分困難,因此本文提出了基于小波包能量曲率差的損傷識別指標。

表4 古木結構的自振頻率

3.3.2 小波包能量曲率差

(1) 損傷工況1的小波包能量曲率差

對古木結構損傷工況1損傷前后梁上27到35節點的豎向加速度響應進行小波包分解，選擇小波函數Db20，分解層數為4層，得到16個小波包系數和能量值，分析前8個能量的小波包能量曲率差如圖4所示。

從圖4(b)、(f)、(i) 可以看出，損傷工況1的損傷發生在30到32節點之間的損傷單元52、53位置，這正好是梁假定的損傷單元52、53所在位置，與損傷工況1假定的損傷位置完全吻合，可以判定在此位置發生了損傷，說明小波包能量曲率變化量可以用于古木結構的損傷定位。

由于篇幅有限，損傷工況2損傷前后的小波包能量曲率差就不一一列舉了。

(2)損傷工況1、2的小波包能量曲率差(前8個分量疊加)。

從圖4看出損傷工況1的小波包能量曲率差的前8個能量圖中只有部分圖形能確定結構的損傷位置，在此基礎上將損傷工況1、2的小波包能量曲率差前8個分量進行疊加來分析損傷定位的效果如圖5所示。

圖4 損傷工況1的小波包能量曲率差

由圖5(a)、(b)看出，將損傷工況1、2損傷前后的小波包能量曲率差的前8個分量進行疊加，在節點31處發生突變，其值最大，損傷位置最明顯，隨著節點離跨中越遠，損傷指標值越小，到梁兩端的27和35節點時，損傷指標值達到最小，說明該指標對損傷位置較敏感，小波包能量曲率差的前8個分量疊加比小波包能量曲率差更能識別古木結構的損傷位置，它可以作為損傷識別指標進行古木結構的準確定位，并且損傷程度越大，該指標值越大。

3.4 噪聲對損傷識別的影響

在結構的健康監測過程中,由傳感器采集的數字信號難免要受到外界噪聲的干擾，對于信號中存在的噪聲,一般假定為高斯白噪聲,在原有信號基礎上疊加服從正態分布均值為0的量作為測量噪聲,通過信號的信噪比(SNR)來衡量信號的噪聲水平,信噪比定義為:

(23)

式中：AS為信號x(n)的均方根，AN為噪聲y(n)的均方根。

圖5 損傷工況1、2的小波包能量曲率差疊加

下面以損傷工況1為例來分析在信噪比SNR=10、20、30、40、50情況下對損傷定位效果的影響。在有限元分析得到的梁上各節點的豎向加速度信號中分別加入不同分貝的高斯白噪聲來研究噪聲對損傷識別的影響。仍選用db20小波函數進行小波包分解,分解層數為4，不同噪聲水平下的小波包能量曲率差的前8個分量疊加如圖6所示。

從圖6看出測試數據含有白噪聲對信號的高頻部分影響較大,故在含噪信號的小波包分解中應選用低頻概貌信號和低頻細節信號作為損傷識別的判別依據。當信噪比SNR小于或者等于20 db時, 該損傷指標受噪聲影響較大,對損傷定位效果影響較大,已不具備損傷定位的能力。當信噪比SNR等于30 db時，受噪聲影響較小，已能基本進行損傷定位了。當信噪比SNR大于或者等于40 db時, 損傷定位能力已不受噪聲影響，該損傷指標對損傷定位效果與無噪聲信號基本相當,說明隨著信噪比的提高, 該損傷指標對損傷識別的敏感性逐漸增加，損傷定位效果越好，該損傷指標具有一定的抗噪聲干擾能力。當信噪比較小時,受噪聲的影響較大,因此需對含噪聲信號進行消噪處理,盡可能還原為原始信號,才能保留損傷信息，以便對古木結構的損傷進行準確定位。

圖6 不同噪聲水平下小波包能量曲率差疊加

3.5 損傷程度的判定

由圖5看出, 對于無噪聲信號在同一損傷位置不同損傷程度時，損傷指標柱狀圖基本相似,只是數值大小有差別。因此設想對于同一損傷位置不同損傷程度，若能找到損傷程度和損傷指標之間的函數關系，繪出其關系曲線，就能由該關系曲線對損傷程度進行判斷了。

損傷程度的判定方法:首先判定結構是否存在損傷,若存在損傷，再確定損傷的具體位置,然后針對該損傷位置,對不同損傷工況進行數值模擬,得到損傷位置上不同損傷工況下的損傷指標,再運用matlab進行數值擬合,繪出損傷指標與損傷程度之間的關系曲線,由該關系曲線判定該損傷位置的損傷程度。

仍假設古木結構梁跨中出現損傷,對梁跨中損傷程度10%、15%、20%、25%、30%、35%、40%、45%、50%各損傷工況進行有限元分析,得到梁上各節點的豎向加速度信號，計算梁跨中31節點的損傷指標—小波包能量曲率差的前8個分量疊加如表5所示，采用matlab進行數值擬合，找到損傷位置31節點處損傷程度和損傷指標之間的函數關系，繪制其函數關系曲線如圖7所示。

表5 損傷指標

圖7 損傷程度和損傷指標之間的關系曲線

從圖7得出損傷指標和損傷程度之間的函數關系式為:y=0.062 36x2+0.947 9x+10.13，對古木結構損傷前后進行有限元模擬，得到梁上各節點的豎向加速度信號進行小波包分解，求出損傷指標,由該損傷指標就能在圖7中找到相應的損傷程度了。

圖8 損傷工況3的損傷指標

下面驗算其適用性：假設梁跨中損傷程度為18%,對損傷工況3進行數值分析,得到梁上各節點的豎向加速度響應信號進行小波包分解,計算并繪制損傷指標如圖8所示，得到31節點的損傷指標y= 48.53,由圖7可以逆推出x=18.35%,其相對誤差為1.9%,可見其誤差非常小，因此該函數關系式能對古木結構的損傷程度進行較準確的識別。

4 結 論

本文以西安鐘樓為工程依托，對隨機激勵作用下的古木結構梁上各節點的加速度響應信號進行小波包分解,提出了小波包能量曲率差損傷識別指標,通過此指標進行古木結構的損傷定位，得出結論：

(1) 在無噪聲干擾下，該損傷指標對于古木結構的損傷定位比較敏感,可以準確判定古木結構的損傷位置，并且損傷指標隨損傷程度的加大而增大。

(2) 該損傷指標在高斯白噪聲干擾下,當信噪比SNR大于或者等于40 db時，損傷識別能力已不受噪聲影響，能對古木結構的損傷進行準確定位,說明該損傷指標具有一定的抗噪聲干擾能力,在實際工程應用中能夠取得較好的定位識別效果。

(3) 得出了損傷指標和損傷程度之間的函數關系式，用其進行損傷程度的判斷并驗算其適用性，為研究環境激勵下西安鐘樓的損傷預警提供了理論依據。

參 考 文 獻

[1]Coifman R R, Wickerhauser M V. Entropy-based algorithms for best basis selection [J]. IEEE Trans. Inf. Theory, 1992, 38:713-718.

[2]Doebling S W, Farrar C R, Prime M B. A summary review of vibration-based damage identification methods [J]. Shock Vib.Dig, 1998, 30(2): 91-105.

[3]丁幼亮,李愛群,鄧揚. 面向結構損傷預警的小波包能量譜識別參數[J].東南大學學報(自然科學版), 2011,41(4):824-828.

DING You-liang, LI Ai-qun, DENG Yang. Parameters for identification of wavelet packet energy spectrum for structural damage alarming[J]. Journal of Southeast University (Natural Science Edition), 2011,41(4):824-828.

[4]劉濤,李愛群,丁幼亮,等. 基于小波包能量譜的結構損傷預警方法試驗研究[J]. 振 動 與 沖 擊，2009,28(4)：4-9.

LIU Tao, LI Ai-qun, DING You-liang, et al. Experimental study on structural damage alarming method based on wavelet packet energy spectrum[J]. Journal of Vibration and Shock,2009,28(4)：4-9.

[5]鄧揚,丁幼亮,李愛群.基于小波包分析的拉索損傷聲發射信號特征提取[J].振 動 與 沖 擊，2010，29(6)：154-158.

DENG Yang, DING You-liang, LI Ai-qun. Feature extraction of acoustic emission signals for cable damage based on wavelet packet analysis[J]. Journal of Vibration and Shock, 2010，29(6)：154-158.

[6]范穎芳，胡志強，周晶,等. 基于小波包分析的肋拱橋結構損傷狀態研究[J].工 程 力 學，2008,25(5)：182-188.

FAN Ying-fang, HU Zhi-qiang, ZHOU Jing, et al. Study on damage state of ribbed arch bridge using wavelet packet analysis [J]. Engineering Mechanics, 2008,25(5)：182-188.

[7]HAN Jian-Gang, REN Wei-Xin, SUN Zeng-Shou. Wavelet packet based damage identification of beam structures [J]. International Journal of Solids and Structures, 2005,42: 6610-6627.

[8]余竹，夏禾, Goicolea J M,等. 基于小波包能量曲率差法的橋梁損傷識別試驗研究[J]. 振動與沖擊,2013,32(5)：20-25.

YU Zhu，XIA He，Goicolea J M, et al. Experimental study on bridge damage identification based on wavelet packet energy curvature difference method [J]. Journal of Vibration and Shock, 2013,32(5)：20-25.

[9]李愛群,丁幼亮.工程結構損傷預警理論及其應用[M].北京:科學出版社, 2007.

[10]孟昭博.西安鐘樓的交通振動響應分析及評估[D].西安：西安建筑科技大學，2009.

[11]孟昭博, 袁俊, 吳敏哲,等. 古建筑高臺基對地震反應的影響[J]. 西安建筑科技大學學報(自然科學版), 2008,40(6):835-840.

MENG Zhao-bo,YUAN Jun,WU Min-zhe, et al. The influence of high-station base of ancient building on its seismic responses [J]. J.Xi’an Univ. of Arch. & Tech(Natural Science Edition), 2008, 40(6):835-840.

[12]趙均海，俞茂宏，高大峰,等.中國古代木結構的彈塑性有限元分析[J].西安建筑科技大學學報，1999,31(2):131-133.

ZHAO Jun-hai, YU Mao-hong, GAO Da-feng, et al. FEM analysis on the elasto-plasticity of ancient wooden structure [J]. J.Xi’an Univ. of Arch. &Tech, 1999, 31(2):131-133.

[13]韓廣森.城市軌道交通微幅振動對古建筑的影響[D].西安：西安建筑科技大學，2011.