cankao基于Hilbert_Huang变换和随机子空间识

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对式

24、 (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. 3 E MD 和 SSI 结合识别模态参数利用 E MD 对每个测点的振动响应数据进行

25、处理 , 得到每个测点的前 n 阶模态响应 x j (t , 且是平稳的随 机信号 。 当所有测点的振动响应数据都经过处理后 , 则可以运用 SSI 方法进行模态参数识别 。由于此时识别所用的数据是平稳的单阶模态响应信号 , 根据 SSI 方法的特点 , 既可以更好地发挥该方法 的优势 , 也能克服其不足 。 计算得到的稳定图具有明显的稳定轴 , 各阶模态之间的影响大大减少 , 因此模态 参数的识别结果更准确 。65地 震 工 程 与 工 程 振 动 第 30卷4 基于振动台试验数据的模态参数识别利用本文提出的上述两种方法分别对一 12层钢筋混凝土框架模型振动台试验的加速度记录数据进行 了分析

26、处理 , 识别了结构的模态参数 。该模型为单跨 12层钢筋混凝土框架 , 几何相似关系为 1/10,原型及模型结构 、 测点布置及试验方案详 见文献 10。 选取 El Centr o 波激励下顶层测点 A 7的 x 方向加速度响应数据为研究对象 。图 1所示为原始加速度时程 , 图 2为原始信号的傅里叶谱 。 对原始信号进行 E , 得到结构的第 一阶模态响应 x 1(t , 对 x 1(t 进行 H ilbert 变换 , 得到第 11t 3。 选取测点 A 5为参考点 , 运用 NExT 方法对 x 1(t 进行处理 , x 1(x 1(t 得到 ln A 1(t 的与时间 t 的关系以

27、及由最小二乘法对 ln A 1(t , 可以识别出结构的第 2、 3阶模态参数 。 图 6为运用 E MD 和 SSI 此外 , 借助大型通用有限元软件 ANSYS, 对该模型结构也进行了有限元分析 。75第 1期 韩建平等 :基于 H ilbert -Huang 变换和随机子空间识别的模态参数识别58 地 震 工 程 与 工 程 振 动 第 30 卷 模态参数的识别结果以及有限元分析结果如表 1 所示 。图 7 为利用 SSI识别得到的结构 X 方向前三阶 振型与有限元分析结果的对比 。 表 1 12 层 RC框架模型的模态参数识别结果 Table 1 Identified modal pa

28、ram eters of 12 2storey RC frame model 模态阶数 1 2 3 FFT HHT + NExT EMD + SSI 频率 /Hz 4. 1758 14. 5203 27. 6171 频率 /Hz 4. 1678 14. 5622 27. 9869 阻尼比 / % 6. 80 6. 37 4. 89 频率 /Hz 4. 0847 14. 8668 27. 4989 阻尼比 / % 6. 04 6. 19 4. 01 有限元分析 频率 /Hz 3. 8028 14. 4281 27. 8172 图 7 SSI振型识别结果和有限元分析结果的比较 Fig . 7 Co

29、mparison of first 3 mode shapes betw een SSI and finite element analysis 表 1 结果表明 ,用 FFT、 HHT +NExT、EMD + SSI这三种方法识别得到的结构前三阶固有频率基本一致 , 本文两种方法识别的阻尼比基本吻合 。与该模型结构的有限元分析结果相比 ,第 1 阶固有频率偏差较大 ,其 余各阶频率比较接近 。图 7 结果表明 ,识别的振型与有限元分析的振型吻合较好 ,特别是第 3 阶振型的识别 结果 ,说明该方法对识别难度较大的振型具有较强的识别能力 。 5 结语 HHT方法是一种自适应性的时频分析方法 ,

30、 不仅非常适合于处理非平稳及非线性数据 , 而且具有很好 的时频局域化分析的能力 ,同时 HHT也提供了一种很好的识别结构阻尼的方法 。从理论上讲 ,通过 EMD 方 法得到的平稳信号可以很好地运用 NExT进行处理 ,且得到的自由衰减信号不会含有其它模态分量 ,因此通 过 H ilbert变换以及模态分析的基本理论可以识别出结构的单阶模态阻尼 。这种方法理论上直观明了 ,实现 过程简单方便 。另外 ,通过 EMD 方法 ,可以得到每个测点的单阶模态响应 ,而且是平稳的随机信号 , 再结合 SSI方法进行模态参数识别 ,计算得到的稳定图具有明显的稳定轴 , 各阶模态之间的影响大大减小 , 因此

31、模 态参数的识别结果更准确 。 12 层钢筋混凝土框架模型的模态参数识别结果表明 ,本文提出的两种方法的识别结果同目前广泛应用 的 FFT方法的识别结果非常接近 ; 与有限元分析结果相比 ,固有频率基本一致 ,振型拟合良好 。由此表明本 文提出的这两种方法是可行的 、 实用的 。 1994-2010 China Academic Journal Electronic Publishing House. All rights reserved. 第 1 期 韩建平等 : 基于 H ilbert - Huang变换和随机子空间识别的模态参数识别 59 需要注意的是 ,若结构是时不变系统 ,则较难体

32、现出 HHT方法时频局域化分析的优越性 ; 用 HHT方法 识别模态参数时需要根据结构响应的傅里叶谱确定结构固有频率的大致范围 ,也还没有完全摆脱傅里叶变 换的影响 。另外 ,环境激励下的结构模态参数识别方法 (如本文选用的 SSI和 NExT一般均假设激励为白噪 声信号 ,而本文选用地震动激励下模型结构的加速度响应是违背这一假设的 , 这对识别结果的精准性具有 一定的影响 。 参考文献 : 1 L i Q S, W u J R. Tim e frequency analysis of typhoon effects on a 79 - storey tall building J . Jou

33、rnal ofW ind Engineering and Industrial Aero2 dynam ics, 2007, 95 ( 12 : 1648 - 1666. 2 韩建平 , 王飞行 , 李 慧 . 基于振动台试验的模态参数识别算法比较研究 J . 华中科技大学学报 : 城市科学版 , 2008, 25 ( 3 : 57 - 60. 3 Huang N E, Shen Z, Long S R, et al . The emp irical mode decomposition and H ilbert spectrum for nonlinear and nonstationary tim e series analysis J . Proceedings of the Royal Society A: Mathematical , Physical and Engineer