基于HHT的泄流結(jié)構(gòu)損傷在線監(jiān)測方法研究

練繼建 ，李成業(yè)，劉 昉 ，馬 斌

(天津大學(xué)水利工程仿真與安全國家重點實驗室，天津 300072)

1 前言

在水利工程中，泄流結(jié)構(gòu)的工作條件極其復(fù)雜，由于水流的強烈紊動，產(chǎn)生的脈動壓力作用在結(jié)構(gòu)物上，極有可能造成結(jié)構(gòu)物的強烈振動，甚至導(dǎo)致結(jié)構(gòu)的破壞。為此，必須尋找有效的監(jiān)測方法，以快速、準(zhǔn)確地發(fā)現(xiàn)泄流結(jié)構(gòu)中存在的隱患，有效預(yù)報其安全狀況。現(xiàn)有的損傷檢測方法大都事先假設(shè)已經(jīng)獲得了比較詳盡、準(zhǔn)確的有限元模型信息，根據(jù)實測結(jié)構(gòu)的模態(tài)參數(shù)進(jìn)行損傷識別，這往往會給實際問題的解決帶來困難?；谛盘柼幚砑夹g(shù)的損傷識別是當(dāng)前結(jié)構(gòu)健康監(jiān)測領(lǐng)域的一種新方法，該方法直接對結(jié)構(gòu)振動信號進(jìn)行分析，從中提取完備信息并以此對結(jié)構(gòu)損傷進(jìn)行識別，具有較好的適用性和敏感性，并有助于克服由噪聲、模態(tài)參數(shù)識別精度、損傷指標(biāo)敏感性等因素對損傷結(jié)果的影響。泄流結(jié)構(gòu)損傷過程具有強隨機性和耦合性特點，加之環(huán)境背景噪聲的干擾作用，使得振動響應(yīng)中往往含有非線性成分、奇異點及不規(guī)則的突變部分，應(yīng)用傳統(tǒng)的基于傅立葉變換的信號處理方法存在一定的局限性。

1998年，N.E.Huang等人提出了 Hilbert-Huang 變換 (HHT)[1，2]，該方法由經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition，EMD)和希爾伯特變換(Hilbert Transform，HT)組成。通過經(jīng)驗?zāi)B(tài)分解方法對結(jié)構(gòu)的振動響應(yīng)信號進(jìn)行預(yù)處理，可以把非平穩(wěn)的原始信號分解成一系列具有不同特征尺度的平穩(wěn)信號，把原信號的能量分配到各個固有模態(tài)函數(shù)分量上，而且這種分解是自適應(yīng)的，因而更能反映損傷的本質(zhì)信息。目前，HHT法已被廣泛應(yīng)用于結(jié)構(gòu)健康監(jiān)測與損傷識別領(lǐng)域[3～6]。

筆者等提出了一種基于HHT的泄流結(jié)構(gòu)損傷在線監(jiān)測方法，利用信號能量譜構(gòu)造能量分布向量，進(jìn)而定義結(jié)構(gòu)異常指標(biāo)，通過計算相鄰時間段內(nèi)的異常指標(biāo)，實現(xiàn)結(jié)構(gòu)損傷的在線監(jiān)測和識別。

2 HHT算法的基本原理

N.E.Huang等人定義滿足下面兩個條件的函數(shù)為固有模態(tài)函數(shù):a.整個數(shù)據(jù)序列中，極值點的數(shù)量與過零點數(shù)量相等或至多相差1;b.信號關(guān)于時間軸局部對稱。

實測泄流結(jié)構(gòu)振動信號通常不滿足IMF條件，N.E.Huang等人提出經(jīng)驗?zāi)B(tài)分解法。對任意原始信號x(t)，EMD方法的計算過程是:找出x(t)所有的極大值點，應(yīng)用三次樣條函數(shù)擬合成原始信號的上包絡(luò)線;找出x(t)所有的極小值點，應(yīng)用三次樣條函數(shù)擬合成原始信號的下包絡(luò)線;計算上下包絡(luò)線的均值m1(t)，將原始信號x(t)減去m1(t)得到一個新的信號h1(t)。通常情況下，h1(t)并不滿足一個IMF分量的條件，此時將h1(t)看成原信號，重復(fù)上述步驟，得到第k次篩選的數(shù)據(jù)h1k(t)=h1(k-1)(t)-m1k(t)。此時需要有一個終止準(zhǔn)則來判斷h1k(t)是不是一個IMF分量，否則過多地重復(fù)處理過程，會導(dǎo)致IMF分量失去意義。這個準(zhǔn)則可以用兩個連續(xù)的處理結(jié)果之間的標(biāo)準(zhǔn)差Sd表示，即:

實踐證明，當(dāng)Sd的值在0.2～0.3時，既能保證基本模態(tài)函數(shù)的線性和穩(wěn)定性，又能使其具有相應(yīng)的物理意義。當(dāng)h1k(t)滿足Sd值的要求，則h1k(t)為原信號的第一階IMF分量，記為c1(t)，將原始信號x(t)減去c1(t)，重復(fù)上述過程，直到最后一個殘余信號rn(t)不可再分解為止。最后，原始信號x(t)可以表示成:

文獻(xiàn)[7]應(yīng)用EMD算法分析爆破振動信號，驗證了該方法的完備性和正交性。應(yīng)用EMD分解可得到信號x(t)的多個IMF分量的組合，對ci(t)作Hilbert變換:

式(3)中，PV代表柯西主值(cauchy principal value)，因此構(gòu)造解析信號z(t):

在相位函數(shù)θi(t)的基礎(chǔ)上定義瞬時頻率為:

瞬時頻率是時間t的函數(shù)，它表示某一時刻信號能量在頻域分布的量度。

3 基于HHT的結(jié)構(gòu)損傷在線監(jiān)測指標(biāo)

3.1 信號的能量分布向量

對原始信號的每個IMF分量進(jìn)行Hilbert變換后，可得Hilbert譜的表達(dá)形式:

式(6)中，i=1，2，…，n表示 IMF 分量的個數(shù)，省略了殘余分量rn(t);Re表示取實部。如果Hω，()t對時間積分，就可以得到Hilbert邊際譜:

邊際譜表達(dá)了每個頻率在全局上的幅度，它代表了統(tǒng)計意義上的全部累加幅度。如果振幅的平方對時間積分，就可以得到Hilbert能量譜:

Hilbert能量譜表達(dá)了每個頻率在整個時間長度內(nèi)所積累的能量，將信號的頻率區(qū)間[ωmin，ωmax]分成n個頻帶 ωi(i=1，2，3，…，n)，則各個頻帶上的能量可表示為E(ωi)，進(jìn)而定義能量分布向量ηi為:

3.2 結(jié)構(gòu)損傷異常指標(biāo)

結(jié)構(gòu)發(fā)生損傷引起固有頻率的變化，其振動響應(yīng)能量在頻域以及各階模態(tài)上的分布也會發(fā)生改變，且改變量和損傷位置及程度有關(guān)。利用信號Hilbert能量譜構(gòu)造結(jié)構(gòu)能量分布向量ηi，再進(jìn)一步用它來構(gòu)造結(jié)構(gòu)異常指數(shù)，對損傷是否發(fā)生進(jìn)行識別。設(shè)、分別為無損傷和未知情況下的結(jié)構(gòu)能量分布向量，i=1，2，3，…，n。定義和的互相關(guān)系數(shù)為 Φ()，構(gòu)造損傷指標(biāo)N[8]:f

當(dāng)Nf等于0或者非常小，表明響應(yīng)中模態(tài)成分沒有發(fā)生改變或僅發(fā)生微小變化，無損傷發(fā)生;當(dāng)結(jié)構(gòu)出現(xiàn)損傷，能量在各階模態(tài)上的分布會發(fā)生明顯改變，模態(tài)成分匹配度降低，Nf將遠(yuǎn)大于0。

結(jié)構(gòu)出現(xiàn)損傷，能量在頻域上的分布必然發(fā)生變化，導(dǎo)致能量分布向量ηi峰值所對應(yīng)的頻率發(fā)生變化，構(gòu)建損傷指標(biāo)Rf:

4 水電站導(dǎo)墻結(jié)構(gòu)在線監(jiān)測數(shù)值仿真試驗

文章以某工程導(dǎo)墻結(jié)構(gòu)為例，導(dǎo)墻截面尺寸如圖1所示。實際工程導(dǎo)墻順?biāo)鞣较蛎?4 m設(shè)一條結(jié)構(gòu)縫，因此取24 m長的導(dǎo)墻段建立有限元模型，模型中地基模擬深度為54 m，兩側(cè)寬度各取30 m，地基采用全約束，將與水流方向垂直的兩個導(dǎo)墻面加上法向約束，沿高程自上而下依次布置1#～6#測點，有限元模型及測點布置如圖2所示?；炷敛牧厦芏热? 400 kg/m3，動彈性模量取35 GPa，泊松比取0.167，在計算時將地基材料屬性與導(dǎo)墻取為一致。

圖1 導(dǎo)墻斷面圖(單位:m)Fig.1 The section of the guide wall(unit:m)

為驗證指標(biāo)Nf和Rf的穩(wěn)定性，根據(jù)泄流結(jié)構(gòu)激勵形式的特點，設(shè)計激勵位置相同但大小不同的白噪聲點激勵、激勵位置和大小均不同的白噪聲點激勵以及激勵位置和大小均不同的白噪聲面激勵三種工況，如表1中工況1～工況3，激勵值及其改變量均由模型試驗[9]確定。由模型試驗和有限元計算結(jié)果可知，導(dǎo)墻在截面A、高程260 m處的動應(yīng)力值最大，如圖3所示，該處出現(xiàn)裂縫的可能性也最大，為驗證指標(biāo)Nf和Rf對不同損傷程度的敏感性，文章設(shè)計了截面A處、裂縫深度分別為導(dǎo)墻厚度1/8、1/4和1/2三種程度的損傷，見表1工況4～工況6。

圖2 導(dǎo)墻有限元模型及測點布置圖Fig.2 The finite element model and measured points locations

圖3 模型試驗最不利工況下導(dǎo)墻動應(yīng)力Fig.3 Guide wall dynamic stress for the most unfavorable condition of test

表1 數(shù)值試驗工況表Table1 The conditions of numerical simulation test

續(xù)表

4.1 無損傷情況下各指標(biāo)穩(wěn)定性分析

對無損工況1～工況3條件下各測點的位移響應(yīng)進(jìn)行EMD分解，應(yīng)用式(6)～(8)計算信號的能量譜，應(yīng)用式(9)即可得到能量分布向量。通過計算分析，各工況下所有測點激勵變化(t=25 s)前后ηi的峰值位置以及在各頻率段上的分布規(guī)律基本相同，由于篇幅限制，文中僅列出1工況下3#、4#測點激勵變化前后的ηi分布圖，如圖4所示。根據(jù)式(10)可計算出各測點不同工況下的Nf值，見表2。從各測點的能量譜圖可以看出，信號能量主要在0～25 Hz范圍內(nèi)，找出t時刻前后能量分布向量ηi各峰值所對應(yīng)的頻率和，由式(11)計算出各測點不同工況下的Rf值，見表2。

圖4 工況1改變激勵前后3#、4#測點的能量分布向量圖Fig.4 Energy distribution vector of point 3#，4#before and after excitation changed for condition 1

表2 工況1～工況3不同測點的Nf和Rf指標(biāo)Table 2 The Nfand Rfof every measured point from condition 1 to condition 3

從表2中可以看出，工況1條件下各測點的Nf值均接近于0，而工況2、工況3激勵位置和大小都發(fā)生改變的情況下，各測點的Nf指標(biāo)均有所增加，但由模型試驗確定的最不利荷載條件下(工況3)，各測點指標(biāo)Nf的最大值為0.089 2，因此本例定義超過此最大值10%的數(shù)值作為損傷預(yù)警指數(shù)，即，當(dāng)Nf＞0.1時，即認(rèn)為結(jié)構(gòu)出現(xiàn)異常。與Nf值受激勵位置和激勵條件變化的影響不同，各工況下各測點信號的Rf指標(biāo)均為0，不因外界條件的變化而變化，只與結(jié)構(gòu)自身的特性相關(guān)，當(dāng)Rf＞0時，即認(rèn)為結(jié)構(gòu)出現(xiàn)異常。在實際泄流結(jié)構(gòu)損傷監(jiān)測過程中，可綜合分析Nf和Rf兩項指標(biāo)值，作為判斷結(jié)構(gòu)出現(xiàn)損傷的依據(jù)。

4.2 泄流結(jié)構(gòu)損傷的判定

4.2.1 在線監(jiān)測信號IMF分量的奇異點和瞬時頻率

文獻(xiàn)[10]和文獻(xiàn)[11]指出結(jié)構(gòu)剛度發(fā)生突變，其振動響應(yīng)EMD結(jié)果中的高階固有模態(tài)函數(shù)在時域上將出現(xiàn)奇異點。泄流結(jié)構(gòu)出現(xiàn)損傷必然導(dǎo)致剛度的降低，結(jié)構(gòu)在線監(jiān)測信號表現(xiàn)出一個小的不連續(xù)性，信號出現(xiàn)了奇異點，直接觀察振動信號很難發(fā)現(xiàn)這一不連續(xù)點，而對結(jié)構(gòu)響應(yīng)進(jìn)行EMD分解所得的高頻IMF信號可以檢測出該信號的奇異點及其發(fā)生的時刻。對4、5、6三個損傷工況下各測點響應(yīng)進(jìn)行EMD分解，分解結(jié)果的第一階IMF分量在25 s時刻均出現(xiàn)了明顯的突變，任取工況4和工況6測點2的響應(yīng)進(jìn)行分析，如圖5所示，表明t=25 s時刻結(jié)構(gòu)剛度可能發(fā)生突變。

圖5 工況4和工況6測點2響應(yīng)EMD分解的第一階IMF分量Fig.5 The first IMF component of point 2 for condition 4 and condition 6 respectively

外界載荷和工況的突然變化也會使EMD的分解結(jié)果出現(xiàn)奇異點，因此信號出現(xiàn)奇異點不足以判斷結(jié)構(gòu)出現(xiàn)損傷。信號IMF分量的瞬時頻率是時間的函數(shù)，反映結(jié)構(gòu)固有頻率隨時間的變化規(guī)律，對于健康的結(jié)構(gòu)其瞬時頻率是不變的，瞬時頻率的衰減反映了結(jié)構(gòu)已受到某種損傷。通過分析，損傷工況下各測點IMF分量的瞬時頻率在損傷前后均有一定的下降。任取工況6測點2分解后第四階IMF分量的瞬時頻率，如圖6所示。從圖6中可以看出，結(jié)構(gòu)響應(yīng)的瞬時頻率在第25 s時突然下降，瞬時頻率的平均值也由4.12 Hz下降到2.77 Hz，這表明此時刻結(jié)構(gòu)發(fā)生某種突然的損傷。

圖6 工況6測點2響應(yīng)EMD分解第四階IMF分量的瞬時頻率Fig.6 The instantaneous frequency of the fourth IMF component from point 2 for condition 6

4.2.2 異常指標(biāo)對損傷的敏感性分析

計算工況4～工況6各測點的能量分布向量ηi，通過分析可知，各工況下所有測點損傷前后ηi值的峰值位置及其在各頻率段上的分布均發(fā)生了明顯變化。由于篇幅限制，文中僅列出5工況下3#、4#測點損傷前后的ηi分布圖，如圖7所示。

按照3.2節(jié)中的方法，分別計算工況4～工況6各測點的Nf和Rf指標(biāo)，見表3。從表3可以看出，損傷特征指標(biāo)Nf和Rf對各種不同損傷模式均較為敏感，可以判斷結(jié)構(gòu)損傷的發(fā)生。當(dāng)結(jié)構(gòu)損傷位置一定時，隨著損傷程度的增加，各測點指標(biāo)值相應(yīng)變大，說明根據(jù)損傷指標(biāo)值的大小可以定性地判斷結(jié)構(gòu)的損傷程度。但是相同工況下，各測點的指標(biāo)值沒有明顯的規(guī)律性，根據(jù)能量定義的特征指標(biāo)Nf和Rf無法判斷損傷發(fā)生的位置，這是因為大型水工結(jié)構(gòu)不同部位的振動情況復(fù)雜，對損傷的敏感程度也不相同。

表3 工況4～工況6不同測點的Nf和Rf指標(biāo)Table 3 The Nfand Rfof every measured point from condition 4 to condition 6

4.3 指標(biāo)抗噪性研究

信號采集過程中不可避免地受到噪聲的影響，為了分析噪聲對損傷指標(biāo)的影響，文章在各工況下1#測點響應(yīng)信號中分別加入標(biāo)準(zhǔn)差為10%、20%、30%、40%、50%、60%、70%、80%倍信號標(biāo)準(zhǔn)差的白噪聲，定義為噪聲水平1～8。EMD分解本身就是一個由高頻到低頻的濾波過程，由文獻(xiàn)[12]可知，白噪聲EMD分解后能量主要集中在第一階分量(約占50% ～70%)，應(yīng)用信號能量譜計算損傷指標(biāo)時，除去第一階IMF分量以及最后幾階低頻干擾分量，計算結(jié)果見表4。

圖7 工況5損傷前后3#、4#測點的能量分布向量圖Fig.7 Energy distribution vector of point 3#，4#before and after damage appearance for condition 5

從表4可以看出，在各種噪聲水平下，無損傷工況的Nf指標(biāo)均小于預(yù)先設(shè)定的預(yù)警指數(shù)0.1，與表2中不含噪聲情況下的Nf值相比沒有明顯的規(guī)律性;而損傷工況下，隨著噪聲水平的不斷提高，信號Nf指標(biāo)逐漸減小，但減小的最大幅度(25.52%)遠(yuǎn)低于噪聲的增幅。信號各工況下的Rf指標(biāo)受噪聲的影響較小，基本與無噪情況下的計算結(jié)果相同。因此，在一定噪聲水平的情況下，文中方法仍可以判斷結(jié)構(gòu)的損傷。

表4 各工況1#測點含噪信號的Nf和Rf指標(biāo)Table 4 The Nfand Rfof point 1#responses disturbed by different noise levels for every condition

基于HHT能量分布指標(biāo)的泄流結(jié)構(gòu)在線監(jiān)測方法無需知道結(jié)構(gòu)的模態(tài)參數(shù)等其他特征信息，也無需激勵信息，只需對任一位置處響應(yīng)進(jìn)行監(jiān)測即可實現(xiàn)，能夠滿足結(jié)構(gòu)健康監(jiān)測的動態(tài)、實時性要求。但對于損傷位置和損傷程度的進(jìn)一步精確辨識，還需進(jìn)一步研究。

5 結(jié)語

文章利用HHT法在描述信號瞬時特性和信號處理方面獨特的優(yōu)勢，提出了一種基于HHT的泄流結(jié)構(gòu)損傷在線監(jiān)測方法。首先應(yīng)用HHT法計算響應(yīng)Hilbert能量譜，構(gòu)造能量分布向量ηi，根據(jù)結(jié)構(gòu)發(fā)生損傷時其振動能量在頻域和各階模態(tài)上的分布也會隨之改變的特點，定義結(jié)構(gòu)異常指標(biāo)Nf和Rf，通過計算相鄰時間段內(nèi)的異常指標(biāo)，實現(xiàn)結(jié)構(gòu)損傷在線監(jiān)測。

應(yīng)用文中方法對一導(dǎo)墻結(jié)構(gòu)進(jìn)行數(shù)值仿真試驗，模擬了不同激勵情況下的三種無損傷工況以及相同損傷位置不同損傷程度下的三種有損工況。分別進(jìn)行了無損傷情況下各指標(biāo)的穩(wěn)定性分析、損傷情況下各指標(biāo)的敏感性分析以及各工況下指標(biāo)的抗噪性分析，結(jié)果表明，文中的異常指標(biāo)能夠較精確地判斷結(jié)構(gòu)的損傷，同時具有一定的抗噪能力，其中指標(biāo)Rf的穩(wěn)定性、敏感性以及抗噪性更為突出。但是應(yīng)用文中方法對損傷位置和損傷程度進(jìn)行精確辨識，還需要進(jìn)一步的深入研究。

[1]Huang N E，Shen Z.The empirical mode decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Proceedings of the Royal Society of London，Series A.1998，454:903-995.

[2]Huang N E，Steven R L.A new view of nonlinear water waves:the Hilbert Spectrum[J].Annual Review of Fluid Mesh，1999，31(2):417-457.

[3]Huang N E，Huang K.基于希爾伯特-黃變換的鐵路橋梁結(jié)構(gòu)健康監(jiān)測[J].中國鐵道科學(xué)，2006，27(l):l-7.

[4]舒忠平.基于希爾波特-黃變換的結(jié)構(gòu)損傷檢測方法研究[D].西安:西北工業(yè)大學(xué)，2006.

[5]石志曉，李 昕，周 晶.損傷檢測的經(jīng)驗?zāi)B(tài)分解法[J].大連理工大學(xué)學(xué)報，2005，45(3):401-404.

[6]Qu W Z，Chen M X，Zeng Y L，et al.Structural damage detection using the Hilbert-Huang transforms method[C]//Li A，Sih G C，Nied H F，eds.Proceedings of International Conference on Health Monitoring of Structure，Materials and Environment.Nanjing:Southeast University Press，2007:293-298.

[7]張義平，李夕兵，左宇軍，等.爆破振動信號的HHT分析與應(yīng)用[M].北京:冶金工業(yè)出版社，2008.

[8]李傳習(xí)，陳富強.基于HHT的結(jié)構(gòu)損傷特征量與異常診斷[J].長沙理工大學(xué)學(xué)報(自然科學(xué)版)，2007，4(3):29-33.

[9]練繼建.向家壩水電站消力池底板與導(dǎo)(隔)墻結(jié)構(gòu)水彈性模型試驗研究[R].天津:天津大學(xué)水利水電工程系，2010.

[10]Xu Y L，Chen J.Structural damage detection using empirical mode decomposition:experimental investigation[J].Journal of Engineering Mechanics，2004，130(11):1279-1288.

[11]Yang J N，Lei Y，Lin S，et al.Hilbert-Huang based approach for structural damage detection[J].Journal of Engineering Mechanics，2004，130(1):85-95.

[12]Zhao Huawu，Norden E H.A study of the characteristics of white noise using the empirical mode decomposition method[J].The Royal Society，Proceedings Royal Society of London，A(2004)460:1597-1611.