某某机械设备有限公司欢迎您!

HHT变换讲义

时间:2020-03-16

  HHT变换讲义_理学_高等教育_教育专区。HHT 变换讲义 1.1 简介 传统的信号处理方法, 如傅立叶分析是一种纯频域的分析方法。它用频率不 同的各复正弦分量的叠加来拟合原函数,也即用 F ?? ? 来分辨 f ?? ? 。而 F ?? ?

  HHT 变换讲义 1.1 简介 传统的信号处理方法, 如傅立叶分析是一种纯频域的分析方法。它用频率不 同的各复正弦分量的叠加来拟合原函数,也即用 F ?? ? 来分辨 f ?? ? 。而 F ?? ? 在有 限频域上的信息不足以确定在任意小范围内的函数 f ?? ? ,特别是非平稳信号在 时间轴上的任何突变,其频谱将散布在整个频率轴上。而且,非平稳动态信号的 统计特性与时间有关, 对非平稳信号的处理需要进行时频分析,希望得到时域和 频域中非平稳信号的全貌和局域化结果。在傅立叶变换中,人们若想得到信号的 时域信息,就得不到频域信息。反之亦然。后来出现的小波(Wavelet)变换通 过一种可伸缩和平移的小波对信号变换,从而达到时频局域化分析的目的。但这 种变换实际上没有完全摆脱傅立叶变换的局限,它是一种窗口可调的傅立叶变 换,其窗内的信号必须是平稳的。另外,小波变换是非适应性的,小波基一旦选 定,在整个信号分析过程中就只能使用这一个小波基了。 HHT(Hilbert-Huang Transform)技术是(1998年由NASA的Norden E Huang 等提出的新的信号处理方法。 该方法适用于非线性非平稳的信号分析,被认为是 近年来对以傅立叶变换为基础的线性和稳态谱分析的一个重大突破。目前HHT 技术已用于地球物理学和生物医学等领域的研究,并取得了较好的结果。 存在的问题 尽管HHT技术在处理非线性、 非稳态信号方面有很大的优势, 但是这个方法 本身还是有许多的问题有待进一步研究。 正如Huang 在文章中指出的那样,对于这种新的信号处理方法,其基的完 备性还需要严密的证明。另外,在做Hilbert变换时出现的边界效应也需要更好的 方法来解决。但是,HHT技术中最严重,也是现今研究的最多的是EMD 分解中 的包络过程。从对EMD分解方法的介绍可以看出,包络线的构造影响着整个分 解的结果, 也决定了后面的Hilbert变换。Huang 采用的三次样条插值来拟和包络 线。在实际应用中,发现这样做会产生严重的边界效应,污染了原始数据。特别 是对短数据而言,这种影响使分析所得的结果失去了原有的意义。对此,Huang 等采用的是根据信号端点处的振幅和频率, 分别增加两组特征波的方式进行数据 延拓。Huang的这种延拓方法已经向NASA申请了专利。除此之外,还有人提出 了其它方法进行端点延拓。 比如国家海洋局的黄大吉等提出的镜像闭合延拓法和 极值延拓法。 其中镜像闭合延拓法是根据信号的分布特性,把镜子放在具有对称 性的极值位置, 通过镜像法把镜内信号映射成一个周期性的环形信号,不存在端 点,从根本上避免了端点效应。而极值延拓法具有镜像闭合延拓法相当的效果, 而不增加信号序列本身的长度, 计算较快,尤其在处理非对称波形信号时比较优 越。还有青岛海洋大学的邓拥军等,利用神经网络分析方法来延拓数据端。这些 方法对边界效应的抑制都有一定的作用,但是也都需要更深一步的研究。 1.2 瞬时频率 频率是个极其重要的物理量,定义为信号周期倒数,其物理含义显而易见。 对于正弦信号,它的频率为恒值。但是对于大部分信号,它的频率是随时间变化 的函数, 故提出瞬时频率概念。 瞬时频率即表征信号在局部时间点上瞬态频率特 性,整个持续期上的瞬时频率反映了信号频率的时变规律。 对于随机时间序列 X(t),对其进行 Hilbert 变换[39],可以得到 Y(t)如下: Y ?t ? ? 1 X ?? ? t ?? ? PV ( ? ? ?? d? ) (1-1) 1 其中,PV 为柯西主值(Cauchy principal value)。该式表示 Y(t)是 X(t)与 的 ?? 卷积。通过这个定义,X(t)和 Y(t)组成了一个共轭复数对,于是可以得到一个解 析信号 Z(t)如下 Z ? t ? ? X ? t ? ? iY ? t ? ? a ? t ? e i? ? t ? (1-2) 其中 ? a ? t ? ? ? X 2 ? t ? ? Y 2 ? t ? ?1 2 ? ? ? ? ? ? Y ?t ? ? ?? ? t ? ? arctan ? ? X ?t ? ? ? ? ? ? ? (1-3) 从理论上讲, 虚部的定义方法有很多种。但是 Hilbert 变换为其提供了一个 唯一的虚部值,这就使得其结果成为一个解析函数。得到了相位,就可以得到瞬 时频率,因为瞬时频率就是相位导数。 ?? d? ? t ? dt (1-4) 1 的卷积,因此它强 ?? 调了 X(t)的局部特性。公式(1-2)的极坐标表达式更进一步澄清了这个表达式的 局部特性:它的幅值和相位的最佳局部匹配随三角函数变化。 从本质上说,公式(1-1)将 Hilbert 变换定义为 X(t)与 1.3 固有模态函数 IMF(Intrinsic Mode Function) 由瞬时频率的物理意义可知,并不是任意的信号都能用瞬时频率来讨论。只 有当信号满足只包括一种振动模式,而没有复杂叠加波的情况时才行。实际上, 定义一个有意义的瞬时频率的必要条件就是要求函数关于局部零平均值对称, 并 且零交叉点和极值点数量相同。基于此种原因,提出了固有模态函数的概念。固 有模态函数满足以下两个条件: (1)整个数据范围内,极值点和过零点的数量相等或者相差一个; (2)在任意点处,所有极大值点形成的包络线和所有极小值点形成的包络线 的平均值为零。 第一个条件是显而易见的, 它类似于平稳过程中传统的稳定且满足高斯分布 的窄带信号条件。 第二个条件把传统的全局条件调整到局部情况。只有满足了这 个条件, 得到的瞬时频率才不会因为不对称波形的存在而引起不规则波动。所以 这一点是得到正确瞬时频率的必要条件。而这一点是必须的,因为这样瞬时频率 就可以不包含由于不对称波形造成的波动。 为了使用瞬时频率定义, 必须要把随机数据归结为 IMF 组件,这样才可以为 每个 IMF 组件定义瞬时频率。 为了将数据归结为所需的 IMF 组件,接下来引入经 验模式分解方法。 1.4 经验模式分解 EMD(Empirical Mode Decomposition) 经验模态分解方法的大体思路是利用时间序列上下包络的平均值确定 “瞬时 平衡位置” ,进而提取固有模态函数[39]。这种方法基于如下假设: (1)信号至少有两个极点—-个极大值和一个极小值; (2)信号特征时间尺度是由极值间的时间间隔来确定的; (3)如果数据没有极值而仅有拐点,可以通过微分、分解、再积分的方法 获得 IMF。 在此假设基础上,Huang 等人进一步指出:可以用经验模态分解方法将信号 的固有模态筛选出来。 经验模式分解过程就是个筛选过程, 实现振动模式的提取。 该方法的基本思路是用波动上、下包络的平均值去确定“瞬时平衡位置” ,进而 提取出固有模态函数。 上、 下包络线是由三次样条函数对极大值点和极小值点进 行拟合得到的。 经验模式分解过程的基本过程可概括如下: (1)寻找信号x(t)所有局部极大值和局部极小值,为更好保留原序列的特性, 局部极大值定义为时间序列中的某个时刻的值, 它只要满足既大于前一时刻的值 也大于后一时刻的值即可。 局部极小值的提取同理,即该时刻的值满足既小于前 一时刻的值也小于后一时刻的值。使用三次样条函数进行拟合,获得上包络线 xmax(t)和下包络线)计算上、下包络线的均值m(t)=[xmax(t)+xmin(t)]/2; (3)用原信号 x(t)减去均值 m(t),得到第一个组件 h(t)=x(t)-m(t);由于 原始序列的差异,组件 h(t)不一定就是一个 IMF,如果 h(t)不满足固有模态函数 两个条件,就把 h(t)当成原始信号,重复(1)-(3),直到满足条件为止,这时满足 固有模态函数条件的 h(t)作为一个 IMF,令 I1(t)=h(t),至此第一个 IMF 已经成功 的提取了。由于剩余的 r(t)=x(t)-I1(t)仍然包含具有更长周期组件的信息,因此可 以把它看成新的信号,重复上述过程,依次得到第二个 I2(t),第三个 I3(t),?, 当 r(t)满足单调序列或常值序列条件时,终止筛选过程,可以认为完成了提取固 有模态函数的任务, 最后的 r(t)称为余项,它是原始信号的趋势项。由此可得 x(t) 的表达式 x ? t ? ? ? I i ? t ? ? r ? t ? ,即原始序列是由 n 个 IMF 与一个趋势项组成。如 i ?1 n 上所述, 整个过程就像筛选过程,根据时间特性把固有模态函数从信号中提取出 来。 1.5 端点延拓 HHT 方法的分析质量很大程度上取决于 EMD 分解的质量,而在应用 EMD 方法 时的一个非常棘手的问题是, 由于信号两端不可能同时处于极大值和极小值,这 样, “筛”过程中构成上下包络的三次样条函数在数据序列的两端就会出现发散 现象(如图 1.1(a)所示,原始信号是中心频率为 1KHz、5 个波峰、以 100 个数 据描绘一周期的窄带信号) ,并且这种发散的结果会随着“筛”过程的不断进行 逐渐向内“污染”整个数据序列而使所得结果严重失线(b)所示) 。 另外,在进行 Hilbert 变换时,信号的两端也会出现严重的端点效应。对于 一个较长的数据序列来讲, 可以根据极值点的情况不断抛弃两端的数据来保证所 得到的包络的失真度达到最小。 但对于一个数据点数少的序列来讲,这样的操作 就变得完全不可行。因此,必须对信号或其极值向外进行延拓,以确保包络线抵 达端点。 端点延拓的目的是确保上、下包络都与端点相交,以便有与每一个信号点 相对应的局部平均值。而上、下包络是由极大值和极小值连结而成的,因此只要 对极大值和极小值进行延拓, 而不必对信号本身进行延拓。极大值和极小值是相 间分布的,同时考虑到样条插值的要求,所以只要在信号左、右两端分别延拓两 个极大值和两个极小值即可。由于端点以外没有信号,任何延拓都是人为的,通 过对不同形式的信号的多次试验, 发现一种端点延拓的方法以端点的一个特征波 为依据进行延拓具有较好的效果,具体做法如下[40]: 原始信号 上下包络线 (a) 拟合结果 (b)EMD 分解 图 3.1 未作端点延拓的分析效果 原始信号 上下包络线 (a) 拟合结果 (b)EMD 分解 图 3.2 以端点的一个特征波为依据进行延拓的分析效果 设离散信号: t ? [t (1), t (2), ???, t (n)] ? [t1, t2 , ???, tn ] X (t ) ? [ x(t1 ), x(t2 ), ???, x(tn )] ? [ x1 , x2 , ???, xn ] (1-5) (1-6) 其采样步长为 ?t , X (t ) 有 M 个极大值和 N 个极小值, 对应的序列下标 I m , I n ) ( 、 时间( Tm , Tn )和函数值( U , V )记为: I m ? [ I m (1), I m (2),?, I m (M )] I n ? [ I n (1), I n (2),?, I n ( N )] Tm (i) ? tIm , Tn (i ) ? tIn , U (i) ? xIm , i ? 1,?, M V (i ) ? xI n , i ? 1,?, N (1-7) (1-8) (1-9) (1-10) 具体作法如下: 1.左端 信号左端第一个特征波包含的信号点数为 k1 ? I m (2) ? I m (1), 当I m (1) ? I n (1) ? k1 ? ? I n (2) ? I n (1), 当I m (1) ? I n (1) ?2 I (1) ? I (1) , 当M =N =1 n ? m 向外延拓的两个极值的位置( Tm , Tn )和数值( U , V )为: Tm (0) ? Tm (1) ? k1?t , Tm (?1) ? Tm (1) ? 2k1?t , Tn (0) ? Tn (1) ? k1?t Tn (?1) ? Tn (1) ? 2k1?t U (0) ? U (1) U (?1) ? U (1) V (0) ? V (1) V (?1) ? V (1) (1-11) (1-12) (1-13) (1-14) (1-15) 2.右端 信号右端第一个特征波包含的信号点数为 k2 ? I m ( M ) ? I m ( M ? 1), ? k2 ? ? I n ( N ) ? I n ( N ? 1), ?2 I ( M ) ? I ( N ) , n ? m Tm ( M ? 1) ? Tm ( M ) ? k2 ?t , Tn ( N ? 1) ? Tn ( N ) ? k2 ?t , Tn ( N ? 2) ? Tn ( N ) ? 2k2 ?t , 当I m ( M ) ? I n ( N ) 当I m ( M ) ? I n ( N ) 当M =N =1 (1-16) 向外延拓的两个极值的位置( Tm , Tn )和数值( U , V )为: U ( M ? 1) ? U ( M ) (1-17) (1-18) (1-19) (1-20) Tm ( M ? 2) ? Tm ( M ) ? 2k2 ?t , U ( M ? 2) ? U ( M ) V ( N ? 1) ? V ( N ) V ( N ? 2) ? V ( N ) 3.当端点的数值比近端点的第一个极大值大或极小值小时, 要进行特殊的处 理,以避免信号落到包络线 , U (0) ? x1 , V (0) ? x1 , 当x1 ? U (1) 当x1 ? V (1) 当xn ? V ( N ) (1-21) (1-22) (1-23) (1-24) Tm ( M ? 1) ? tn , U ( M ? 1) ? xn , 当xn ? U ( M ) Tn ( N ? 1) ? tn , V ( N ? 1) ? xn , 由图 1.2 可见, 以端点的一个特征波为依据进行延拓,分别在信号两端增加 两个极大值和两个极小值, 从而使原始信号被延长的包络线所限制,有效地抑制 了端点效应,使得 EMD 分解得到合理的各个 IMF 模态。 另外, 国家海洋局的黄大吉等还提出了镜像闭合延拓法,即根据信号的分布 特性, 把镜子放在具有对称性的极值位置,通过镜像法把境内信号映射成一个周 期性的环形信号。这样做就不存在端点,从根本上避免了端点效应[40]。青岛海洋 大学的邓拥军等人利用神经网络分析的方法来进行端点延拓, 在数据的两端各得 到两个附加的极大值点和两个附加的极小值点,由此有效的抑制了端点效应[41]。 1.6 EMD 结束准则 根据理论分析IMF必须满足两个条件,但是在实际筛选过程中,严格满足这 两个条件的信号几乎不存在,因此如果以这两个条件为依据判定IMF,可能根本 就得不到结果或者要以冗长的程序执行时间为代价。 因此筛选过程要小心进行, 为了保证IMF组件的调幅和调频都具有物理意义, 同时考虑到程序可行性,必须制定一个结束筛选的准则。习惯上,可以通过标准 偏差 SD 来完成,它可以由连续两个筛选结果得到[18]: ? ? h1? k ?1? ? t ? ? h1k ? t ? SD ? ? ? h12? k ?1? ? t ? t ?0 ? ? T ? ? 2 ? ? ? ? ? (1-25) h1(k-1)(t)和h1k(t)分别表示两个连续的筛选结果,一般来说,SD值越小,所得的IMF 分量的线性和稳定性就越好,典型的SD值为0.2-0.3,这样对于不同筛选过程就 有严格的SD限制了。但是在筛选过程中发现,SD的这个范围过于死板,对很多 信号而言,以0.2-0.3为限制的 SD 值并不科学,因此在实际筛选过程中,可以以 0.2-0.3为参考值,根据实际情况进行适当调整。另外,如果不是进行信号精确 处理,可以考虑控制筛选次数来决定筛选结果,这样不仅保障了程序的有效性, 也控制了程序的执行时间,在具体问题中取得了很好的效果。 1.7 Hilbert 谱 完成了经验模式分解过程就得到了所有可提取的固有模态函数,在 Hilbert 变换的基础上,只要根据公式(1-1)-(1-4)计算瞬时频率就可以了。对固有模态 函数进行 Hilbert 变换后,可以用下面方式表示信号 X ? t ? ? ? a j ? t ? exp i ? ? j ? t ? dt j ?1 n ? ? (1-26) 在这儿没有考虑余项 rn,因为它只是单调函数或常量。尽管 Hilbert 变换可以把 单调函数看成是振动的一部分, 但是其余项中的能量很小,一般情况可以不用考 虑。 由希尔伯特变换得出的振幅和频率都是时间的函数, 如果用三维图形表达幅 值、 频率和时间之间的关系, 或者把振幅用灰度的形式显示在频率-时间平面上, 就可以得到 Hilbert 谱 H ( w, t ) 。 如果把 H ( w, t ) 对时间积分,就得到希尔伯特边际谱 h( w) : h(w) ? ? H (w, t )dt 0 T (1-27) 边际谱提供了对每个频率的总振幅的量测, 表达了整个时间长度内累积的振 幅。另外,作为希尔伯特边际谱的附加结果,可以得到如式(1-28)定义的希尔 伯特瞬时能量: IE (t ) ? ? H 2 ( w, t )dw w (1-28) 瞬时能量提供了信号能量随时间的变换情况。事实上,如果振幅的平方对时 间积分,可以得到希尔伯特能量谱: ES (w) ? ? H 2 (w, t )dt 0 T (1-29) 希尔伯特能量谱提供了对于每个频率的能量的量测, 表达了每个频率在整个 时间长度内所累积的能量。 1.8 HHT 方法的软件实现 HHT变换的整个过程包括EMD与Hilbert变换,整个过程是通过matlab软件编 写的M文件实现的,其程序流程图如图1.3所示。 开始 输入信号 x(t) r(t)=x(t),n=0 x(t) 求出 x(t)的所有极值点 拟和上下包络线 x(t)=h(t) x(t)=h(t) 计算包络线的平均值 m(t) h(t)=x(t)-m(t) h 是否满足 IMF 条件 是 否 n=n+1,Cn (t)=h(t), r (t)=r (t)-Cn (t) 否 Cn(t) 是 否 满 足 EMD 中止条件 是 对 Cn(t)进行 Hilbert 变换 构造解析函数 求瞬时频率 边际谱,希尔波特能量 谱瞬时能量谱等 结束 图 1.3 经验模式分解与 Hilbert 变换流程图 1.9 HHT 方法的优越性 1、HHT方法是一种全新的信号分析方法,能够描绘出信号的时频谱图、边际 谱,能量谱等,是一种更具有适应性的时频局域化分析方法[28]。而传统的信号处 理方法, 如傅立叶变换是一种纯时域的分析方法,它用频率从零到无穷大的各复 正弦分量的叠加来拟合原函数f(t)在每个时刻的值,也即用F(ω )在有限频域上 的信息不足以确定在很小范围内的函数f(t), 特别是非平稳信号在时间轴上的任 何突变,其频谱将散布在整个频率轴上[39]。所以,这种分析方法适用于确定性的 平稳信号,而在非线性、非平稳过程的处理上,傅立叶变换将只作为一种数学变 化手段。而且,非平稳信号的统计特性与时间有关,对非平稳信号的处理需要进 行时频分析, 希望得到时域和频域中非平稳信号的全貌和局域化结果。但在傅立 叶变换中,若想得到信号的时域信息,就得不到频域信息,反之亦然。后来出现 的小波变换通过一种可伸缩和平移的基小波对信号变换, 从而达到时域局域化分 析的目的, 但这种变化实际上并没有完全摆脱傅立叶变换的局限,它是一种窗口 可调的傅立叶变换,其窗内的信号必须是平稳的[12]。 2、HHT局部性能良好而且是自适应的,对平稳信号和非平稳信号都能进行分 析。它没有固定的先验基底,分解完全基于数据本身进行。固有模态函数是基于 序列数据的时间特征尺度得出的不同的数据有不同的固有模态函数, 每个固有模 态函数可以认为是信号中固有的一个模态,所以通过Hilbert变换得到的瞬时频 率具有清晰的物理意义, 能够表达信号的局部特征。而傅立叶变换是以余弦函数 为基底进行信号分解的, 这样难免产生许多实际上并不存在的虚假分量。小波变 换是选定适当的小波基,再对信号进行分解的,且小波基一但选定,在整个信号 分析过程中就只能用这个基。另外,小波变化的解释也不直观,用小波变化时, 为了确定某一信号的局部变化, 即使这一局部变化只是发生在低频范围内,也必 须从高频范围内开始去寻找这一结果,频率越高,小波变换得局域化性能就越好 [28] 。 3、式(1-26)中信号幅值和频率都是时间函数,但同样的信号如果用傅立叶 变换表示的话,可得 X ?t ? ? ? a je j ?1 n i? j t , (1-30) 其中, aj 和 ?j 是常量。公式(1-26)和(1-30)的差别在于,IMF表示了广义的 傅立叶扩展, 可变的幅值和瞬时频率不仅强化了信号信息还使之适用于非平稳信 号。通过IMF可以清楚的区分调幅和调频,这样就打破了傅立叶变换中固有幅值 和固有频率的限制,允许分解出的IMF的幅值随时间变化,使信号分析更加灵活 方便。同时,由于Hilbert变化通过微分法来定义瞬时频率,因此不需要大量的 信号点来定义振动。 4、Hilbert谱把各IMF分量的幅值以灰度的形式表示在频率—时间图上,其 中幅值以点的灰度表示:点越亮,幅值越大;点越暗,幅值越小。Hilbert谱三 维灰度图的形式对各分量的瞬时频率和瞬时振幅都进行了比较确切的刻画[42]。 比 传统的频谱图更直观清晰。 5、很多信号不是围绕水平位置振动的,这就需要从测量数据中去掉这个背 景场,使之成为水平直线]。EMD方法就可以有效地提取或取出数据 序列的均值,消除序列的趋势项,把复杂的数据分解成若干线性、平稳的模态, 且不改变原数据的物理特性。即式(1-31) 。而小波变换和傅立叶变换,它们的 基底都是以水平轴为平衡位置的波动现象。 ? I (t ) ? x(t ) ? r (t ) i ?1 i n (1-31)