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 基金項目 :甘肅省科技攻關項目 (2GS057-A52-008 作者簡介 :韓建平 (1970- , 男 , 教授 , 主要從事工程結構抗震、 減震控制及結構損傷識別研究 . E 2mail:j phanlut . cn文章編號 :1000-13

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

3、這兩種方法的可行性和實用性。關鍵詞 :H ilbert -Huang 變換 ; 隨機子空間識別 ; 經驗模態分解 ; 自然激勵技術 ; 模態參數識別 ; 振動臺試驗中圖分類號 :O329; T U311. 3 文獻標志碼 :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引言土木工程結構的損傷診斷和健康監測具有極其重要的現實意義 , 而振動信號的分析與處理是結構損傷 診斷和健康監測研究和實踐中的一個主要環節 , 同時也是難點所在 。 通過信號處理 , 識別結構的模態參數是 當前國內外研究的熱點問題之一 。傳統的信號處理方法主要是基于傅里葉變換 , 這是一種純頻域的分析方法 , 它用不同頻率的各復正弦分 量的疊加擬合原函數 。 傅里葉頻譜散布在頻率軸上 , 。 后來出現的小波變換 (W avelet Transfor m , 達到時頻局域化 分析的目的 , 它實質上是一種窗口可調的傅里葉變換 ,

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

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

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

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

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

15、條件 , 則將 h1作為新的原始數據 , 重復上述篩選過程 , 直到得到的 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 分量 , 代表原始數據序列中最高頻的組分 。將 r1作為原始數據重復以上過程 , 得到 x (t 的第 2個 I M F 分量 C2, 重復循環 n 次 , 得到信號 x (t 的 n 個 I M

16、F 分量 。當 r n 成為一個單 調函數不能再從中提取滿足 I M F 條件的分量時 , 循環結束 。因此 , 任何一個信號 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 構建 x (t 的解析信號 Y (t ,Y (t =x (t +i x (t =A (t e i (t (6 式中 , A (t 為瞬時幅值 ,

17、(t 為瞬時相位 , i =(-1 1/2。瞬時頻率 (t 為(t =d (t /d t (7 根據式 (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 時刻對應頻率 j 的瞬時幅值 , 因此 , 在時 -頻域內的幅值分布 , 記為 A (t, ; x ,A (t, ; x =6n j =1A j (t, j (9 式中 , A (t, ; x 稱為 x (t 的 H ilbert -Huang 譜 。因此 , 如果存

18、在 n 階本征模態函數 (I M F , 通過 H ilbert 變換可以得到任意 t 時刻的 n 個不同的頻率分 量 , H ilbert -Huang 譜 A (t, ; x 是信號 x (t 的幅值關于時間 t 和頻率 的分布 。2 隨機子空間識別線性振動系統的二階微分動力方程可以表示為確定性的連續狀態空間方程 6, 7,x (t =A c x (t +B c u (t y (t =C x (t +D u (t (11 在實際的信號采集過程中 , 經過離散采樣并引入隨機噪聲后 , 確定性的連續狀態空間方程變為離散的隨 機狀態空間方程 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 是離散狀態向量 , y k =y (k t 和 u k =u (k t 分別是離散的輸入和輸出向量 , A =e A c t 是離散狀態矩陣 , B =A-I A c -1Bc是輸入影響矩陣 , C 是輸出影響矩陣 , D 是直接傳輸矩陣 , wk是處理過程誤差 ,v k 是測量誤差 。 假定 w k 、 v k 為均值為零且互不相關的白噪聲 , 相關函數為 , 55第 1期 韓建平等 :基于 H ilbert -Huang 變換和隨機子空間識別的模態參數識別E w p v p (w T q v T q = Q S

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

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

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

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