cankao基于Hilbert_Huang變換和隨機(jī)子空間識

1、第 30卷 第 1期 2010年 2月 地 震 工 程 與 工 程 振 動 JOURNAL OF E ARTHQUAKE E NGI N EER I N G AND E NGI N EER I N G V I B RATI O N Vol . 30No . 1 Feb . 2010 收稿日期 :2009-01-16; 修訂日期 :2010-01-13 基金項(xiàng)目 :甘肅省科技攻關(guān)項(xiàng)目 (2GS057-A52-008 作者簡介 :韓建平 (1970- , 男 , 教授 , 主要從事工程結(jié)構(gòu)抗震、 減震控制及結(jié)構(gòu)損傷識別研究 . E 2mail:j phanlut . cn文章編號 :1000-13

2、01(2010 01-0053-07基于 Hilbert -Huang 變換和隨機(jī)子空間識別的模態(tài)參數(shù)識別韓建平 1, 2, 李達(dá)文 1, 1(11蘭州理工大學(xué) , 21同濟(jì)大學(xué) , 摘 要 :基于 H ilbert 2Huang 法。 方法一是基于 ilbert , 通過經(jīng)驗(yàn)?zāi)B(tài)分解和 H ilbert 變換提取信號的; 方法二是基于經(jīng)驗(yàn)?zāi)? 通過經(jīng)驗(yàn)?zāi)B(tài)分解對信號進(jìn)行預(yù)處理 , 進(jìn)而運(yùn)用隨機(jī)子空間識別方法利用這兩種方法 , 通過對一 12層鋼筋混凝土框架模型振動臺試驗(yàn)測點(diǎn)加速度記錄的處理 , 識別了該模型結(jié)構(gòu)的模態(tài)參數(shù)。 識別結(jié)果與傳統(tǒng)的基于傅里葉變換的識別結(jié)果及有限元分析結(jié)果的對比驗(yàn)證了

3、這兩種方法的可行性和實(shí)用性。關(guān)鍵詞 :H ilbert -Huang 變換 ; 隨機(jī)子空間識別 ; 經(jīng)驗(yàn)?zāi)B(tài)分解 ; 自然激勵技術(shù) ; 模態(tài)參數(shù)識別 ; 振動臺試驗(yàn)中圖分類號 :O329; T U311. 3 文獻(xiàn)標(biāo)志碼 :AM oda l param eter i den ti f i ca ti on ba sed on H ilbert 2Huangtran sform and stocha sti c subspace i den ti f i ca ti onHAN J ianp ing 1, 2, L IDa wen 1,WANG Feixing 1(1. I nstitute o

4、f Earthquake Pr otecti on and D isasterM itigati on, Lanzhou University of Technol ogy, Lanzhou 730050, China;2. State Key Laborat ory of D isaster Reducti on in Civil Engineering, Tongji University, Shanghai 200092, China Abstract:T wo app r oaches are p r oposed for modal para meter identificati o

5、n of civil engineering structures based on H ilbert 2Huang transf or m (HHT and st ochastic subs pace identificati on (SSI . The first app r oach is based on HHT and natural excitati on technique (NExT . E mp irical mode decompositi on (E MD and H ilbert transfor m (HT are used t o extract the insta

6、ntaneous characteristics of the original signal . Then, NExT and basic modal analysis theory are app lied t o identify modal para meters . The second app r oach is based on E MD and SSI . Each single modal re 2s ponse is obtained thr ough p r ocessing the original signal by E MD, then the modal para

7、 meters are identified by SSI . Finally, the original signals fr om the shaking table test of a 122st orey reinf orced concrete fra me model are p r ocessed and modal para meters are identified by these app r oaches, res pectively .I dentificati on results and comparis on with the results fr om trad

8、iti onal fast Fourier transfor m (FFT and finite ele ment analysis indicate that the p r oposed ap 2p r oaches are feasible and p ractical .Key words:H ilbert 2Huang transf or m; st ochastic subs pace identificati on; e mp irical mode decompositi on; natural excitati on technique; modal para meter i

9、dentificati on; shaking table test引言土木工程結(jié)構(gòu)的損傷診斷和健康監(jiān)測具有極其重要的現(xiàn)實(shí)意義 , 而振動信號的分析與處理是結(jié)構(gòu)損傷 診斷和健康監(jiān)測研究和實(shí)踐中的一個主要環(huán)節(jié) , 同時(shí)也是難點(diǎn)所在 。 通過信號處理 , 識別結(jié)構(gòu)的模態(tài)參數(shù)是 當(dāng)前國內(nèi)外研究的熱點(diǎn)問題之一 。傳統(tǒng)的信號處理方法主要是基于傅里葉變換 , 這是一種純頻域的分析方法 , 它用不同頻率的各復(fù)正弦分 量的疊加擬合原函數(shù) 。 傅里葉頻譜散布在頻率軸上 , 。 后來出現(xiàn)的小波變換 (W avelet Transfor m , 達(dá)到時(shí)頻局域化 分析的目的 , 它實(shí)質(zhì)上是一種窗口可調(diào)的傅里葉變換 ,

10、性的 , 1, (例如峰值 拾取法 、 頻域分解法等 2H ilbert -Huang 變換 (H Huang 等提出的時(shí)間序列信號處理的新方 法 , 目前己在海洋 、 地震 3。 HHT 由經(jīng)驗(yàn)?zāi)?態(tài)分解 (E mp on, E MD 及 H ilbert 變換 (H ilbert Transf or m 兩部分組成 , 其核心是 E MD 。 E MD , 它依據(jù)數(shù)據(jù)本身的時(shí)間尺度特征進(jìn)行分解 , 是 自適應(yīng)的 , 因此更適合于處理非線性非平穩(wěn)數(shù)據(jù) 。本文提出通過 E MD 和 H ilbert 變換提取結(jié)構(gòu)振動信號的 瞬時(shí)特性 , 在此基礎(chǔ)上利用自然激勵技術(shù) (Natural excit

11、ati on technique, NET 和模態(tài)分析的基本理論識別結(jié) 構(gòu)的模態(tài)參數(shù) 。 該方法不僅能夠準(zhǔn)確地識別結(jié)構(gòu)的固有頻率 , 還提供了一種較好的識別結(jié)構(gòu)阻尼的方法 。 隨機(jī)子空間識別 (St ochastic subs pace ldentificati on, SSI 是目前較為先進(jìn)的環(huán)境激勵條件下結(jié)構(gòu)模態(tài)參 數(shù)識別的時(shí)域方法 , 其數(shù)學(xué)模型為狀態(tài)空間方程 , 通過求解狀態(tài)空間方程的系統(tǒng)矩陣和輸出矩陣得到結(jié)構(gòu)的 模態(tài)參數(shù) 。 由于該方法假定輸入為白噪聲 , 而且計(jì)算得到的穩(wěn)定圖包含了結(jié)構(gòu)各階模態(tài)信息 , 給準(zhǔn)確識別模 態(tài)參數(shù)帶來了困難 。 因此 , 本文提出首先利用 E MD 處理振

12、動數(shù)據(jù) , 使其只包含某一階的模態(tài)信息 , 而且是平 穩(wěn)的隨機(jī)信號 , 然后利用 SSI 方法識別模態(tài)參數(shù) , 以得到較為理想的識別結(jié)果 。1 Hilbert -Huang 變換利用 H ilbert 變換識別結(jié)構(gòu)的模態(tài)參數(shù)是經(jīng)典的動力學(xué)反問題 , 但 H ilbert 變換只能提取信號的主分量 , 而其它分量則被忽略或被處理成一個畸變系數(shù) , 因此對非線性多自由度體系則不適用 。 而 E MD 方法能將非 線性非平穩(wěn)的振動信號分解成一系列的本征模態(tài)函數(shù) (I ntrinsic mode functi ons, I M Fs , Huang 對 E MD 的定 義使分解得到的 I M F 可以

13、很好地進(jìn)行 H ilbert 變換 3, 4。HHT 方法主要分為兩步 :第一步是通過 E MD 對信號進(jìn)行預(yù)處理 , 得到一系列的 I M Fs, 每階 I M F 都很好 地滿足 H ilbert 變換的條件 ; 第二步是對分解得到的 I M Fs 進(jìn)行 H ilbert 變換并構(gòu)建解析信號 , 求得各自的瞬時(shí) 頻率 , 并畫出時(shí)頻圖 , 進(jìn)而根據(jù)模態(tài)分析的基本理論 , 識別結(jié)構(gòu)的模態(tài)參數(shù) 。1. 1 E MDE MD 是將原始信號相鄰峰值點(diǎn)間的時(shí)延定義為時(shí)間尺度 , 據(jù)此把信號分解成包含不同時(shí)間尺度的若干 階 I M F 和一個余量的和 5。 設(shè)一原始數(shù)據(jù)序列為 x (t , 篩選步驟如

14、下 :首先 , 找出數(shù)據(jù)序列的所有局部極大值和極小值 , 并分別用三次樣條函數(shù)進(jìn)行插值 , 得到原序列的上 、 下 包絡(luò)線 , 然后對上 、 下包絡(luò)線上的每個時(shí)刻的值取平均 , 得到瞬時(shí)平均值 m1, 用原數(shù)據(jù)序列 x (t 減去 m 1, 得 到 :h 1=x (t -m 1(1 每一個 I M F 必須滿足兩個條件 :整個數(shù)據(jù)段內(nèi) , 極值點(diǎn)的個數(shù)和過零點(diǎn)的個數(shù)必須相等或相差最多不 能超過 1個 ; 任何一點(diǎn) , 由局部極大值點(diǎn)形成的包絡(luò)線和由局部極小值點(diǎn)形成的包絡(luò)線的平均值為 0。在實(shí) 際運(yùn)用時(shí) , 其平均值的絕對值小于某一個很小的特定值即可 。若 h 1 不滿足以上 I M F 的兩個

15、條件 , 則將 h1作為新的原始數(shù)據(jù) , 重復(fù)上述篩選過程 , 直到得到的 h 1k滿足I M F 的條件 , 從而得到第一個 I M F 分量 C 1,45地 震 工 程 與 工 程 振 動 第 30卷h 11=h 1-m 11C 1=h 1k =h 1(k -1 -m 1k (2 從原始信號中分離出 C1, 得到r 1=x (t -C 1(3 C 1 即為信號 x (t 的第 1個 I M F 分量 , 代表原始數(shù)據(jù)序列中最高頻的組分 。將 r1作為原始數(shù)據(jù)重復(fù)以上過程 , 得到 x (t 的第 2個 I M F 分量 C2, 重復(fù)循環(huán) n 次 , 得到信號 x (t 的 n 個 I M

16、F 分量 。當(dāng) r n 成為一個單 調(diào)函數(shù)不能再從中提取滿足 I M F 條件的分量時(shí) , 循環(huán)結(jié)束 。因此 , 任何一個信號 x (t 都可以分解為 n ,x (t =6n j =1C j (t r t (4式中 , 分量 C1, C 2, , C n , 余量 r n 則表示信號 x (t 的中心趨 勢 。1. 2 H ilbert 變換x t , 其 H ilbert 變換記為 x (t ,x (t =HT x (t =+- (t -d (5 構(gòu)建 x (t 的解析信號 Y (t ,Y (t =x (t +i x (t =A (t e i (t (6 式中 , A (t 為瞬時(shí)幅值 ,

17、(t 為瞬時(shí)相位 , i =(-1 1/2。瞬時(shí)頻率 (t 為(t =d (t /d t (7 根據(jù)式 (4 (7 , 得到以下解析方程 ,Y (t =6n j =1A j (t e i j (t d t =6n j =1A j (t, j e i j (t d t (8式中 , Aj(t, j 為第 j 階 I M F 在 t 時(shí)刻對應(yīng)頻率 j 的瞬時(shí)幅值 , 因此 , 在時(shí) -頻域內(nèi)的幅值分布 , 記為 A (t, ; x ,A (t, ; x =6n j =1A j (t, j (9 式中 , A (t, ; x 稱為 x (t 的 H ilbert -Huang 譜 。因此 , 如果存

18、在 n 階本征模態(tài)函數(shù) (I M F , 通過 H ilbert 變換可以得到任意 t 時(shí)刻的 n 個不同的頻率分 量 , H ilbert -Huang 譜 A (t, ; x 是信號 x (t 的幅值關(guān)于時(shí)間 t 和頻率 的分布 。2 隨機(jī)子空間識別線性振動系統(tǒng)的二階微分動力方程可以表示為確定性的連續(xù)狀態(tài)空間方程 6, 7,x (t =A c x (t +B c u (t y (t =C x (t +D u (t (11 在實(shí)際的信號采集過程中 , 經(jīng)過離散采樣并引入隨機(jī)噪聲后 , 確定性的連續(xù)狀態(tài)空間方程變?yōu)殡x散的隨 機(jī)狀態(tài)空間方程 6,x k +1=A x k +B u k +w ky

19、 k =C x k +D u k +v k(12式中 , xk=x (k t 是離散狀態(tài)向量 , y k =y (k t 和 u k =u (k t 分別是離散的輸入和輸出向量 , A =e A c t 是離散狀態(tài)矩陣 , B =A-I A c -1Bc是輸入影響矩陣 , C 是輸出影響矩陣 , D 是直接傳輸矩陣 , wk是處理過程誤差 ,v k 是測量誤差 。 假定 w k 、 v k 為均值為零且互不相關(guān)的白噪聲 , 相關(guān)函數(shù)為 , 55第 1期 韓建平等 :基于 H ilbert -Huang 變換和隨機(jī)子空間識別的模態(tài)參數(shù)識別E w p v p (w T q v T q = Q S

20、S T R p q (13式中 , E 是數(shù)學(xué)期望 , pq 是 Kronrcker delta 函數(shù) 。 環(huán)境激勵是不可測量的隨機(jī)激勵 , 且強(qiáng)度基本和噪聲影響相似 , 較難將兩者區(qū)分清楚 , 因此 , 可以將式 (12 中的輸入項(xiàng) u k 和噪聲項(xiàng) w k 、 v k 合并得到離散隨機(jī)子空間方法的基本方程 ,x k +1=A x k +w ky k =C x k +v v (14 基于式 (14 , 可以有不同的方法實(shí)現(xiàn)振動系統(tǒng)的識別 , 而 SSI 。 SSI , 并 以此來預(yù)測未來 。 SSI 采用有效的數(shù)學(xué)工具如矩陣的 QR (value decompositi on, S VD 以

21、及最小二乘等來識別系統(tǒng)狀態(tài)矩陣 , 3 3. 1 HHT 假設(shè) x (t 。 通過 x (t 的傅里葉譜 , 可以得到結(jié)構(gòu)前 n 階固 有頻率 (1, 2, 3, , n 的初步估計(jì) , 如第 1階固有頻率 1L 11H 。 要得到第 j 階模態(tài)響應(yīng) , 則選取合適的帶通濾波頻率 j L 、j H , 并將 x (t 通過該帶通濾波器 , 依次類推 , 得到 n 個時(shí)間序列信號 , 記為 x j (t (j =1, 2, 3, , n 。 分別對每一個時(shí)間序列信號 x j (t 進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解 , 得到的第 1階 I M F 一般就是對應(yīng)的結(jié) 構(gòu)第 j 階模態(tài)響應(yīng) , 記為 x j (t 。

22、因此 , 由環(huán)境激勵得到的結(jié)構(gòu)加速度響應(yīng)被分解成結(jié)構(gòu)的每一階模態(tài)響應(yīng) , 而且響應(yīng)是平穩(wěn)的隨機(jī)信號 , 適合進(jìn)行 H ilbert 變換 。對每個 x j (t 進(jìn)行 H ilbert 變換并構(gòu)建解析方程 Y j (t , 根據(jù)公式 (6 、(7 可以得到第 j 階模態(tài)的瞬時(shí)頻 率 j (t , 即第 j 階模態(tài)響應(yīng)的瞬時(shí)頻率在不同 t 時(shí)刻的分布 , 對 j (t 取平均值 , 則可得到該階模態(tài)頻率 j 。3. 2 HHT 方法同 NExT 結(jié)合識別阻尼比由于得到的每階模態(tài)響應(yīng) x j (t 是平穩(wěn)的隨機(jī)信號 , 同樣適合運(yùn)用 NExT 或者隨機(jī)減量技術(shù)得到對應(yīng)的 自由振動模態(tài)響應(yīng) , 記為

23、 x j (t 。根據(jù)結(jié)構(gòu)模態(tài)分析的基本理論 , 結(jié)構(gòu)的自由振動模態(tài)響應(yīng)可以表示成解析方程x j (t =R j e -j j t cos (d t -j (15 式中 , j 是第 j 階固有頻率 , j 為第 j 階阻尼比 , d y 為第 j 階有阻尼固有頻率 , R j 是與模態(tài)初始條件和阻尼比 有關(guān)的量 。 對式 (15 中的 x j (t 進(jìn)行 H ilbert 變換并構(gòu)建解析信號 Y j (t ,Y j (t =A j (t e i j (t =R j e -j j t e i (dj-j (16可以得到幅值 A j (t 為 A j (t =R j e -j j t (17對式 (17 兩邊取自然對數(shù) ln A j (t =-j j t +ln R j(18 式 (18 表明 , ln A j (t 與時(shí)間 t 是線性關(guān)系 , 其圖形是斜率為 -j j 的直線 。由于 j 可以用前述方法求得 , 則阻尼比 j 也可以確定 。需要指出的是 , 當(dāng)為非小阻尼情況時(shí) , ln A j (t 和時(shí)間 t 之間不是嚴(yán)格的線性關(guān)系 , ln A j (t 的值將在一條 直線附近振蕩 8, 9, 可以運(yùn)用最小二乘擬合法將 ln A j (t -t 圖形擬合成一條直線并進(jìn)而求得阻尼比 j 。3