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