基于ESN延拓与HMM修正的端点效应处理及其应用

罗红玲，陈启卷，王卫玉，江 文，刘宛莹，席 慧，安宇晨

(武汉大学 水力机械过渡过程教育部重点实验室，武汉 430072)

自然界中许多信号，都是典型的非线性非平稳信号。传统的信号分析与处理方法大多以傅里叶变换为理论依据演变而来，傅里叶变换属于全局估计法，这类方法不仅难以精确描述频率随时间的变化，且自适应性较差[1,2]。将希尔伯特变换(Hilbert Transform, HT)应用到信号处理中[3]可以使分析分辨率达到观测信号的采样分辨率，但HT变换的对象必须是窄带信号，否则会产生与瞬时频率的物理意义相冲突的负频率[4]。基于此，1998年Huang等人提出了EMD方法[5]，自驱动地将信号按照不同频带划分，实现对信号进行平稳化处理，为非平稳信号进行HT奠定了基础，该过程被称为希尔伯特黄变换(Hilbert Huang Transform, HHT)。HHT具有优异的自适应性，且经过HHT得到的时频谱能够精确地反映信号的能量在时间和频率上的分布规律，因而在语音、生物医学信号、机械振动等领域得到了广泛的应用。

但Huang在提出HHT的同时指出，端点效应是限制其发展最为突出的问题之一。为此国内外学者提出了各类算法来抑制端点效应，在现有的算法中，延拓法得到普遍认可。根据原理的不同，可以将延拓法分为三大类：基于端部数据延拓、基于内部数据延拓、预测延拓。每一类算法虽然在不同场合都具有一定的效果，但都存在各自的局限性，目前难以找到一种具有普适性的方法，同时兼顾延拓准确性和运算效率。

对于故障诊断领域，有研究表明，通过分析振动信号来充分挖掘机组的故障信息，可对80%的潜在安全问题实现早期预警和防范[6-8]，尤其水电机组转速较低，骤发恶性事故概率极小，从故障征兆的出现发展到事故是一个缓变过程。EMD分解结果的准确性关系到振动信号的时频特征的有效提取。

预测延拓法较其他两类方法的准确性更高，因此考虑沿用预测延拓法以保留在其分解精度上的明显优势。神经网络预测[9]和支持向量机回归延拓[10]的应用较为广泛，两者都是将输入从低维空间映射到高维空间，在高维空间利用处理线性问题的理论和方法去处理问题，且各自仍存在局限性，并无明显优劣之分[11]。于是有研究者试图将各类方法结合运用，如支持向量机预测与镜像延拓结合[12]、支持向量机和窗函数结合[13]、B样条插值结合神经网络[14]和神经网络与镜像延拓结合[15]等，但这一类方法存在计算量偏重和复杂模型解决简单问题的弊端。

Jaeger在2001年提出的ESN[16,17]网络简单高效，克服了梯度消失等问题，其自带的反馈连接比支持向量机等核方法的前馈结构更适合时序相关的动态预测[18]。由此，本文借用ESN对原序列进行初步预测以保证延拓精度。然而，由于采用了基于数据点预测极值点的方式，使用的递归多步预测策略难以避免预测误差的累积。ESN模型产生的预测误差由两个部分组成，一是模型对序列变化趋势的预测能力(偏高、偏低或正常)，二是模型对随机波动的预测误差。因此，本文通过引入生成模型HMM[19,20]找到ESN模型状态序列的变化规律，对原始预测值进行调整，用来有效修正由于递归多步预测策略导致的误差累积，从而达到兼顾预测精度和运算效率的目的。

1 常见端点效应抑制算法及抑制效果评价

1.1 端点效应及常见抑制方法

端点效应是指在EMD分解过程中，利用3次样条函数求取信号上下包络线时，端点处由于缺乏极值点约束产生发散，分解结果出现失真的现象。根据原理的不同，可将数据延拓法分为以下4大类。

1.1.1 基于端部数据的延拓法

该方法根据信号端部有限个极值点呈现出的特征：①人为地在端点处增加极大极小值点，如边界局部特征尺度延拓法[21]、极值平均法[22]和根据斜率进行延拓的方法有SBM法[23]、改进的ISBM法[24]和RO-SBM法[25]以及极值平移法[26]、基于数据/极值联合对称延拓[27]等；②按照端点至首个极值之间的波形特征，通过延拓波形间接实现极值点的外延，如正弦波法[28]。这类方法虽然延拓简便、计算量小且实现了对端点处包络线的约束，但本质上都只是考虑端点处波动趋势，忽视了近端点处与信号内部规律的联系，缺乏对信号内部特征的真实反映，延拓精度难以保证。因此，适用于随机性较强的信号、长度较短的序列或对运算精度要求不是很高但对运算效率要求较高的场合。处于上述任意一种情况，不论选取哪种延拓方法，都难以实现信号的较准确预测，不可避免地引入误差，任意一种方法并不比其他方法更具适用性，因此可以在尽量保证分解效果的前提下采取最简单、快捷的办法。

1.1.2 基于内部数据延拓

这一类方法利用内部数据点的波动特征向外延拓：①直接在信号内部寻找相似波形进行延拓，如最大相关波形延拓[29]和自适应三角波匹配延拓[30]等；②利用极值符号[31]、互相关函数[32]等算法找出序列内与之最为匹配的部分向外延拓；③将整个序列进行镜像延拓[33]。虽然该方法考虑了延拓数据与内部波形之间的联系，但计算量大小由序列长度决定，且实测信号一般既有确定性成分(或确定性混沌成分)，又含有随机成分，难以保证在原序列中能匹配到相似波形，因此对规律性较强的信号抑制效果更好。

1.1.3 预测延拓法

预测延拓是通过数学方法或人工智能算法对数据序列进行预测，延拓部分数据与原序列连接处过渡平滑且预测准确性较高，适用于呈现出一定周期性和局部随机性较强的信号或对于分解精度要求较高的场合，但各类方法均存在其固有局限性。如ARMA模型[34]的延拓效果受制于模型复杂度和参数设置，不适用于随机性较强信号。最大Lyapunov指数预测[35]对具有混沌特征的信号更有效。而灰度预测模型[36]在形成预测公式时规定的已经条件是否足够合理有待探讨。支持向量机[37]参数与核函数的选择对预测结果影响较大。

1.2 端点效应抑制效果评价

通常，观察EMD分解图和希尔伯特时频谱两端是否发散可以直观展示抑制效果，为了避免定性观察所得结果主观性太强，下面列举出常用的两个评价指标来定量评估端点效应抑制效果。

1.2.1 正交系数(Index of Orthogonality, IO)

Huang通过大量实验指出，EMD方法具有正交性，即信号分解得到的各IMF分量之间具有不相关、互相正交的性质。检验正交性时，需要将残余分量作为一个额外的分量，因此原始信号x(t)可表示为：

(1)

式中：K表示IMF分量的个数(包括残余分量)；imfk(t)表示信号经过EMD分解得到的第k个IMF分量。

为了检验分解得到的各IMF分量间的正交性，需要将式(1)两边平方得到式(2)。如果分解结果完全正交，则式(2)右边的第二部分的交叉项应该为0，因此可以定义正交系数如式(3)。

(2)

(3)

式中：K和imfk(t)的含义同上；T表示信号长度。

由式(3)可以看出，正交系数越小，EMD端点效应抑制效果越好。

1.2.2 能量比θ

EMD分解具有完备性，即分解后各IMF分量的能量和理论上等于原信号能量。通过式(4)可计算长度为T数据序列x(t)的有效值RMS作为信号的能量值。

(4)

但由于端点效应的存在，端点处失真、虚假分量的产生都会使得EMD分解后IMF的能量和增加。因此可以通过计算分解前后的能量比θ来评估端点效应的影响程度[38]。

(5)

式中：RMSk表示信号经过EMD分解得到的第k个IMF的有效值；RMSoriginal为原始信号序列的有效值。

2 基于ESN-HMM的端点效应处理方法

ESN-HMM法主要分为3个步骤，首先根据原始数据序列构造样本集，利用训练好的ESN模型对数据点向外延拓。对比理论输出与实际模型输出得到检验误差序列，对误差序列进行HMM建模求出误差预测值。将初步预测值与误差预测值求和，即可实现对预测数据的误差校正，流程图如图1所示。

图1 基于ESN延拓与HMM修正的端点效应处理流程Fig.1 End point effect processing flow based on ESN-HMM

2.1 基于ESN的数据序列延拓

2.1.1 ESN的网络结构

ESN由输入层、储备池和输出层组成，结构如图2所示。储备池由大量的、稀疏的、随机连接的非线性节点组成，具有短期记忆功能。

图2 回声状态网络结构图Fig.2 The structure of echo state networks

设在n时刻的输入信号、输出信号和储备池状态分别为u(n)、y(n)、x(n)，Win、Wres、Wout、Wback分别代表输入层和储备池的连接权值、储备池内部的连接权值、储备池和输出层之间的连接权值以及反馈连接权值。储备池的状态更新方程和网络的输出计算由式(6)、式(7)给出。ESN预测模型的整体框架如图3所示。

x(n+1)=f[Winu(n+1)+Wresx(n)+Wbacky(n)]

(6)

y(n+1)=fout{Wout[u(n+1),x(n+1)]}

(7)

式中：f表示内部神经元的激活函数，一般为双曲正切函数；fout为输出单元的神经元激发函数，通常为对称型Sigmoid函数或线性函数。

图3 基于ESN的数据序列预测框架Fig.3 Echo state networks based framework for data sequence prediction

2.1.2 ESN网络的训练、预测及参数优选

在网络的训练过程中，Win、Wres和Wback都是随机生成且一旦初始化就保持不变，仅需对Wout进行训练，

而储备池到输出层是线性关系，因此存在全局唯一最优解且只需要通过求解线性回归问题即可得到。确定Wout的过程，主要包括以下4个步骤：

(1)网络初始化。假设有K个输入单元，L个输出单元。初始化ESN网络储备池的拓扑结构参数，主要包括储备池规模N、储备池内部连接权值谱半径SR。

(2)状态更新。在训练样本[u(n),d(n),n=1～l0+l]中选取长度为l0的初始化样本对储备池状态进行更新，用以去掉初始暂态的影响。

(3)状态收集。从l0+1时刻开始收集状态，经过l步长得到状态矩阵A(l,K+N)和由相应的模型实际输出y(n)构成的Y(l, 1)。

(4)权值计算。状态收集结束后，输出权值Wout通过最小化问题来求解，即最小化期望输出和模型实际输出的均方误差：

(8)

式中：d(n)为期望输出；l0表示用于更新状态的样本个数；l表示用于状态收集的样本个数。

应用训练好的ESN网络采用递归多步预测策略进行预测。

2.2 基于HMM的预测误差修正

2.2.1 基于HMM的误差预测模型

HMM描述由一个隐藏的马尔科夫链随机生成不可观测的状态序列I=(i1,i2,…,iT)(i表示不同时刻的状态，T表示离散时间序列的长度)，再由各状态生成一个观测序列O=(o1,o2,…,oT)(o表示不同时刻的观测值)的过程[34]。一个标准HMM，可用五元组来描述λ=(Q,V,π,A,B)。HMM模型的基本参数包括两组状态集合和3组概率：隐藏状态集合Q=(q1,q2,…,qM)，M表示隐藏状态个数，q表示不同的隐藏状态；观测状态集合V=(v1,v2,…,vN)，N表示观测状态个数，v表示不同的观测状态；初始状态概率向量表示为π=[πi]1×N,πi=P(i1-qi)、状态转移矩阵A=[aij]M×M,aij=P(it+1=qj|it=qi)以及观测概率矩阵B=[bj(k)]M×N,bj(k)=P(ot=vk|it=qj)，其中N表示观察状态个数)。对预测误差进行HMM建模包括以下几个步骤：

(1)确定状态序列。将ESN预测延拓方法得到的数据偏高、偏低或正常情况视为3种不同的隐藏状态，即M=3，根据预测误差的大小，按照式(9)可以确定状态序列I=(i1,i2,…,iT)。

(9)

式中：it表示t时刻的隐藏状态；cre为预测误差与真实数据的比值。

(2)求取观测序列。利用K-means算法将预测误差分为偏低、略偏低、接近真实数据、略偏高、偏高5类，即N=5，cl、csl、csh、ch为5个类别的边界阈值，c1、c2、c3、c4和c5为5类中心。通过边界阈值可以形成新的观测值序列O=(o1,o2,…,oT)，如式(10)所示。

(10)

式中：ot表示t时刻的观测值。

(3)HMM参数估计。以I和O作为hmmestimate概率估计函数的输入，求得状态转移矩阵A和观测概率矩阵B。

(4)误差预测。预测t+1时刻的误差e，需要先通过式(11)求出t时刻模型所处状态概率分布pi，然后再通过式(12)对观测值的离散化阈值H进行还原。

pi=SiTi

(11)

e=piH

(12)

式中：Si为t时刻的状态；Ti为在状态i向其他状态转移的概率分布。

2.2.2 预测误差修正后处理

利用HMM得到的预测误差e对ESN模型得到的初步预测值Y进行修正，最终预测值Pre为：

Pre=Y+e

(13)

3 算例分析与验证

3.1 仿真信号对比分析

为了检验ESN-HMM端点效应抑制方法的有效性，现以式(14)生成的仿真信号为例进行分析。仿真信号由3个正弦信号组成，其时域波形如图4所示。

x(t)=3 sin(2 πf1t)+5 sin(2 πf2t)+

5 sin(2 πf3t)，t∈[0,4]

(14)

式中：f1、f2和f3表示3个正弦信号的频率值，分别为6、2.5和0.5 Hz。

图4 仿真信号时域波形Fig.4 Time domain waveform diagram of the simulated signal

当仿真信号x(t)未加任何边界处理经EMD分解，分解结果如图5所示。从图5可以看出，由于端点效应的不良影响，分解结果中各模态分量两端严重“飞翼”，IMF1和IMF2截去两端形变部分仍能识别出其对应频率成分6 Hz和2.5 Hz。但0.5

图5 无边界条件的EMD分解结果和时频谱Fig.5 EMD and time-frequency distribution without any boundary condition

Hz的低频分量由于极值点个数较少且时间间隔较大，波形失真使其失去了其原有的物理意义，且在时频谱上显示和其他成分两端产生了交叠。根据EMD分解过程可知，特征模态函数的筛分是依据高频分量到低频分量的次序进行，随着迭代的深入，由于分解误差会逐渐累积，加上低频分量处插值精度的降低，导致端点效应更为明显。

基于ESN延拓后进行EMD分解的效果和时频谱如图6所示。由图6可知，经ESN延拓的EMD分解结果端点处的“飞翼”现象相较于无延拓，得到了极大的改善。IMF1、IMF2和IMF3分别对应原信号中的6Hz、2.5Hz和0.5Hz三个频率分量，这说明数据延拓添加极值点能够有效地抑制包络线的发散，保证了分解结果的有效性。但分解结果中多出了幅值非常小的虚假分量IMF4，这是因为极值点的预测过程采用的是基于数据点递归预测的策略，迭代产生的误差会随着预测点数的增加而进一步累积，降低了分解结果的准确性。由此，利用HMM模型对ESN延拓的数据进行误差修正，再进行EMD分解得到如图7所示的抑制端点效应后的分解结果和时频谱图，可以看出，应用了ESN-HMM法后不仅能够有效抑制端点效应，HMM误差修正模型的加入使得所添加的极值点更为符合原信号特点，有效抑制了虚假分量的产生，从而能够准确分解出原信号中3个分量的幅值、频率特征。

图6 基于ESN延拓的EMD分解结果和时频谱Fig.6 EMD and time-frequency distribution of ESN

为了进一步检验ESN-HMM方法，从现有算法中选取了典型的斜率法和正弦波法(属于端部数据延拓法)、镜像延拓和自适应波形匹配法(属于内部数据延拓法)进行定量比较。表1列举了各类方法各项指标的具体数值，包括正交系数IO、能量比θ、耗时t和IMF分量个数IMFS。

由表1可知，相较于无边界延拓，经过延拓再分解的信号的正交系数和能量比都有大幅减小，说明延拓算法对于端点效应均能起到抑制作用。其中，斜率法和正弦波延拓法由于算法简便，因而延拓耗时最短，但其延拓依据是信号端部数据的特征，割裂了与信号序列内部的联系，因而准确性略低于另外两类方法。由于原始信号波形具有一定的对称性和规律性，因此基于内部数据延拓一类算法的分解效果也比较好。基于ESN延拓与HMM修正的预测延拓法在分解准确性上表现最优，运算效率虽略低但是仍能满足实时要求。

图7 基于ESN-HMM的预测延拓的EMD分解结果和时频谱Fig.7 EMD and time-frequency distribution of ESN-HMM method

表1 各类方法端点效应抑制效果比较Tab.1 Comparison of inhibition effects of various methods

3.2 实测压力脉动信号分析

为了进一步验证基于ESN延拓与HMM修正的端点效应抑制方法的有效性和工程实用性，以某水电站200 MW混流式水轮机在部分负荷工况下的进口压力脉动信号为例，分析其幅值超标的原因。

机组额定转速为150 r/min，额定转频为2.5 Hz，机组旋转一周传感器采集256个数据，采样频率640 Hz。由于端点效应对短数据序列、极值点个数较少的情况的不良影响更明显，因此仅截取2 048个数据点进行分析，信号的原始波形如图8所示。

图8 实测尾水管进口压力脉动波形Fig.8 Measured pressure fluctuation waveform of draft tube inlet

对压力脉动信号直接进行EMD分解得到的时频谱如图9所示，各模态分量两端“飞翼”严重，分解结果完全失真，失去了原有的物理意义，无法识别出相应的时频分布。

图9 无边界延拓的时频谱Fig.9 Time-frequency distribution without any boundary condition

同样地，对各类延拓算法的分解结果进行对比，各类方法的抑制效果详见表2。

表2 各类方法端点效应抑制效果比较Tab.2 Comparison of inhibition effects of various methods

正弦波延拓和镜像延拓效果的分解结果如图10、图11所示。在所有的延拓方法中，正弦波延拓效果最差，各模态分量两端发散严重，无法获取原信号中的频率成分。镜像延拓法得到的时频谱可以辨认出信号中含有频率值小于转频2.5 Hz的低频分量，但各模态分量的频率谱线之间多处交叠，尤其端点处更为严重，因此难以准确地提取出原信号的时频特征。

图10 基于正弦波延拓的时频谱Fig.10 Time-frequency distribution of Sine wave extension

图11 基于镜像延拓的时频谱Fig.11 Time-frequency distribution of mirror extension

基于ESN预测延拓的时频谱如图12所示，可以看出，低频分量在端点处的交叉重叠程度较小，能够分辨出其中一个低频分量的谱线轮廓，但由于延拓精度不够，虚假分量仍存在。基于ESN延拓与HMM修正后得到的时频图谱线清晰，虚假分量减少，两端振荡微弱，低频脉动分量的幅值特征明显，如图13所示，可以看出此时低频分量占据主导地位。而国内外工程实践说明，在部分负荷工况下，转轮后通常会产生螺旋涡带，脉动频率为机组转频的0.15～0.3倍[39]，与信号分析结果吻合，验证了此时机组尾水管中压力脉动幅值较大主要是受到了低频涡带的影响。

图12 基于ESN预测延拓的时频谱Fig.12 Time-frequency distribution of ESN

图13 基于ESN延拓与HMM修正的预测延拓的时频谱Fig.13 Time-frequency distribution of ENS-HMM

4 结 语

针对水电机组故障诊断对EMD端点效应抑制算法的延拓准确性和运算效率的要求，提出了一种基于ESN延拓与HMM修正的预测延拓法，首先利用ENS预测延拓保证了分解精度，其次运用误差校正的思想减小递归预测导致的误差累积，并通过与典型方法对仿真信号和实测压力脉动信号的对比分析，进一步验证了ESN-HMM方法的有效性。

□