Signal feature extraction method based on MEEMD-HHT for distributed optical fiber vibration sensing system

分布式光纤振动传感系统通过检测振动引起的光纤中光信号的相位以及强度变化实现扰动事件的实时、定位监测，已广泛应用于周界安防、非法施工监测、机械或结构异常分析等多个领域^[1]。随着应用需求的不断细化，在实现扰动事件定位监测的同时，进一步准确判别扰动事件的类型是目前该领域的研究热点。而在模式识别算法实现过程中，信号特征的提取是至关重要的一个步骤，其自身的性能及其对信号特征分类算法复杂度的影响，直接影响分布式光纤振动传感系统模式识别结果的实时性、可靠性和准确性。

目前，国内外研究人员提出了多种分布式光纤振动传感系统特征提取算法，如过零率^[2]、短时平均能量^[3]及全相滤波器^[4]等时域特征提取方法，该类方法较适用于强度检测的直接探测结构，且最终识别准确率较低；在频域特征提取方面，有快速傅里叶变换、Mel频率倒谱系数^[5]等方法，由于该类方法为了在数学上拟合原始数据的非平稳波形不得不引入大量高频的伪谐波分量，导致其对低频能量的低估，对自然环境中的复杂振动信号的特征提取误差较大，获得的特征向量较少。经验模态分解方法和Hilbert变换相结合是一种提取信号瞬时频率和瞬时能量的有效方法，2016年，王亮^[6]采用经验模态分解（Empirical Mode Decomposition, EMD）方法对光纤振动传感系统进行信号特征的提取，针对不同的振动信号进行分解，并将Hilbert变换后求得振动信号的边际谱作为扰动信号的特征；朱程辉等^[7]采用总体经验模态分解（Ensemble Empirical Mode Decompostion, EEMD）的方法对振动信号进行特征提取，并结合其他特征提取方法增加了特征向量种类及精度。2017年，蒋立辉等^[8]采用互补经验模态分解（Complementary Ensemble Empirical Mode Decompostion, CEEMD）方法对获得的固有模态函数（Intrinsic Mode Function, IMF）进行奇异值分解，构建特征向量。2019年，陈勇等人^[9]在EMD算法中引入两次极值波延拓，将端点效应限制在延拓区域中，采用改进后的EMD算法分解信号，并利用分量的相关性提取有效IMF分量；曲洪权等人^[10]采用EEMD方法对入侵信号进行分解，计算各固有模态分量的能量占比，通过方差分析进一步筛选重要的特征分量。2020年，Wang等人^[11]对原始信号和差分信号分别进行EMD分解得到11维能量特征向量；Fouda等人^[12]采用EMD分解和Hilbert变换进一步得到信号的时频熵特征和高、低频段重心频率特征。然而，上述方法在特征提取过程中存在不同程度的模态混叠、伪分量多、重构误差大、计算耗时等问题，增加了特征分类算法复杂度，降低了系统模式识别的准确率或识别效率。

文中采用一种改进的经验模态分解方法（Modified Ensemble Empirical Mode Decompostion，MEEMD）应用于分布式光纤振动传感系统的信号特征提取，可有效解决分解过程中模式混叠以及伪分量问题，引入排列熵的评价机制，优化随机噪声迭代次数，简化了算法，提升了信号特征提取的准确率和提取效率。首先，详细介绍了MEEMD方法的原理，对其数值模型及实验中用于特征提取评价的Hilbert谱和边际谱计算方法进行推导；随后，通过仿真分析，采用正交性和完备性作为评价指标，验证了MEEMD方法相对于EMD、CEEMD方法可有效抑制模态混叠、伪分量问题，具有最小的重构误差，同时耗时相对较少。最后，搭建了实验平台，采用信号发生器模拟单频和混频信号驱动PZT分别作用于传感光纤，采用三种经验模态方法对系统解调的扰动点相位随时间变化曲线进行分解，并采用Hilbert谱和边际谱对结果进行了分析。

目前，EMD方法对振动信号经验模态分解获得的IMF分量中存在较多的伪分量，且有模态混叠现象，影响振动信号特征提取的准确性。CEEMD方法是对EMD方法和EEMD方法^[13]的改进，将一对互为相反数的正负白噪声作为辅助噪声加入源信号中，以消除原来EEMD 方法分解后重构信号中残留的多余的辅助白噪声，同时减少分解时所需的迭代次数，然而白噪声幅度由人为设定，直接影响上述问题解决效果。因此，对经验模态分解方法进行进一步改进，在分解过程中引入基于排列熵的时间序列随机性评价机制，减少了EEMD和CEEMD方法中冗余的集成平均，减小了计算量和重构误差，使IMF分量更具代表性，下面对其原理进行详细介绍。

首先，对相位解调后持续时间为T的非平稳随机扰动信号时间序列加入均值为零的白噪声信号${n_p}(t)$和$ - {n_p}(t)$，得到：

式中：${a_p}$表示第$p$次添加噪声信号的幅值，$p = 1, $$ 2,\cdots,{N_e}$，${N_e}$为加入噪声的总次数。对$x_p^ + (t)$和$x_p^ - (t)$分别进行EMD分解，得到第一组IMF分量序列集合$\left\{ {I_{1p}^ + (t)} \right\}$和$\left\{ {I_{1p}^ - (t)} \right\}$，将其进行累加、平均后，得到第一组IMF分量${I_1}(t)$：

随后采用排列熵（Permutation Entropy, PE）检查与${I_1}(t)$相对应的数字序列${I_1}(i)$是否为异常信号，对其进行延时空间重构^[14]，$i = 1,2,\cdots,k,\cdots,N$，$N$为数字序列总长度，得到如下序列：

式中：$\tau $为延迟时间；$m$为嵌入维数；将每个${R_1}(i) = $$ \{ {I_1}(i),{I_1}(i + \tau ),\cdots,{I_1}(i + (m - 1)\tau )\} $按照数值由小到大进行排列，得到一组新的序列：

若存在数值大小相等的情况，则按照$j$值的大小进行排列，所以任意一个${R_1}(i)$都可以得到一组符号排列${S_1}(g) = \{ {j_{i1}},{j_{i2}},\cdots,{j_{im}}\} $，m个不同的符号共有$m!$种不同的符号组合$g = 1,2,\cdots,k,k \leqslant m!$，${S_1}(g)$只是其中一种排列方式，通过计算可获得每一种符号序列出现的概率${P_{\rm{1}}},{P_{\rm{2}}},\cdots,{P_k}$以及序列的排列熵值：

经标准归一化处理，可以获得振动信号离散序列的随机程度的大小${H_p}$，${H_p}$越大表明时间序列的随机性越强，存在随机噪声等异常现象。需重复上述步骤优化噪声叠加次数，直至${I_1}(t)$不被判定为异常信号，对于其他IMF分量处理方式相同并将其从原始信号中逐个分离，剩余的为残余分量。

针对监测环境振动信号特点设定合适的阈值，将求得的排列熵值与设定的阈值相比较，将超出阈值的IMF分量累加求和，并将其从原始信号中分离，达到去除噪声以及环境因素的干扰的效果，将剩余数据进行EMD分解，最终得到振动信号的各个IMF分量，极大程度上降低了模态混叠现象和伪分量现象。

将分解的各IMF分量进行Hilbert变换，获得与之对应的的幅值和相位信息^[15]。

式中：$I(t)$表示IMF分量的时间序列；$y(t)$表示对应的解析信号。构造的解析函数${\textit{z}}(t)$为：

其中，相位函数：

而相位函数的导数即为瞬时频率：

采用Hilbert谱$H(w,t)$表达时间、频率和幅值之间的关系，对$H(w,t)$进行时间积分，得到Hilbert边际谱$h(w)$：

将Hilbert谱、Hilbert边际谱作为振动信号的特征向量，从而区分出不同类型的振动信息。

实验装置如图1所示。窄线宽、单频相干光源（Coherent Laser，中心波长：1550.12 nm）发出的光首先通过分光比为10∶90的光纤耦合器（Coupler），90%的光作为探测光通过声光调制器（Acousto Optic Modulator, AOM）调制，掺饵光纤放大器（Erbium Doped Fiber Amplifier, EDFA）放大后由光纤环形器（Circulator）注入到传感光纤中，10%的光作为本振光，传感光纤中的瑞利散射光^[16]为${E_{\rm{R}}}(t){{{\rm{exp}}i(}}(\omega + \Delta \omega )t + \varphi (t))$，其中，${E_{\rm{R}}}(t)$为散射光的振幅，$\omega $为光源的角频率，$\Delta \omega $为AOM调制产生的频移量，$\varphi (t)$为光纤折射率改变引起的相位变化；后向瑞利散射光与本振光在2×2耦合器（3 dB coupler）处进行光学混频，由平衡光电探测器（Balanced Photoelectric Detector，BPD）进行光电转换，保留其交流信号，${E_{{\rm{LO}}}}(t)$为本振光振幅，交流部分对散射光的振幅、相位以及偏振状态$\theta $敏感，通过数据采集卡（DAQ）后进行傅里叶积分的数字域处理^[16]，获得振幅和相位信息：

式中：$\Delta T$为调制后脉冲光的宽度，它决定了系统的空间分辨率，通过相位解调算法^[17]可以获得传感光纤的在振动位置处的相位$\varphi \left( t \right)$，相位信息准确还原了真实的振动情况，因此以时间-相位数据作为振动信号特征提取的原始信号。

Figure 1.  Experimental device diagram of distributed optical vibration sensing system

在MEEMD方法实现过程中，时间序列重构延迟时间$\tau $、嵌入维数m以及熵阈值$\eta $的设定对算法的执行效果有着至关重要的影响。

延迟时间$\tau $通过互信息法来判断$I(t)$和$I(t + \tau )$ 之间的相关性，将第一个IMF分量${I_1}(t)$作为原始数据序列R，Q作为${I_1}(t)$的时间延迟$\tau $数据序列${I_1}(t + \tau )$，两者互信息计算公式是：

式中：${P_r}({r_i})$、${P_q}({q_j})$为边缘分布的概率；${P_{rq}}({r_i},{q_j})$为当${{R}} = {r_i},Q = {q_j}$时重构的联合分布概率，延迟时间$\tau $是当互信息函数第一次达到极小值时对应的时间，文中将延迟时间设置为1。

嵌入维数m通过伪近邻法进行确定，逐渐增加m值，计算${I_1}(t)$的虚假最近邻点的比例，直到伪近邻点不再随着m的增加而减少时，m为最小嵌入维数。通过对数据特点的分析以及实验经验总结^{[14, 18]}，文中将嵌入维数设置为6。

为确定排列熵阈值，仿真生成典型的有用信号和噪声信号，采样频率设置为1 kHz，采样时间为2 s，分别计算各自的排列熵值，其中非异常信号包括高频谐波信号${x}_{1}\left(t\right)$、低频谐波信号${x_2}\left( t \right)$、调幅调频信号${x_3}\left( t \right)$、混频信号${x_4}\left( t \right)$、脉冲信号${x_5}\left( t \right)$和分段信号${x_6}\left( t \right)$；异常信号包括间歇噪声${x_7}\left( t \right)$和高斯白噪声${x_8}\left( t \right)$，可分别以公式（14）中表示：

八种信号的排列熵值如表1所示。由表1可知，非异常信号中，高频谐波信号${x_1}\left( t \right)$的排列熵最高，为0.4610；而异常信号中，间歇噪声${x_7}\left( t \right)$的排列熵值最低，为0.6163。参考文献[14]建议排列熵阈值取0.55~0.6较为合适，结合仿真信号排列熵结果考虑，文中将排列熵${H_p}$的阈值设置为0.55。

采用频率为80 Hz的高频周期信号${x_1}$、频率为10 Hz的低频周期信号${x_2}$、调幅信号${x_3}$以及功率为0.2的高斯白噪声信号$n\left( t \right)$进行信号仿真，以模拟解调后的振动信号：

公式（15）中，t=2 s，采样频率为1 kHz，合成的仿真振动信号及其各成分的时域波形图如图2所示。

Figure 2.  Vibration simulation signal and the composition of the time-domain waveform

对上述仿真信号分别进行EMD、CEEMD以及MEEMD分解，分别如图3~5所示。IMFi表示第几个分量，RS表示剩余残余的趋势项。

Figure 3.  EMD decomposition results of simulation signal

Figure 4.  CEEMD decomposition results of simulation signal

Figure 5.  MEEMD decomposition results of simulation signal

由图3~5可以看出，由于存在噪声的干扰，EMD分解存在模态混叠，且虚假分量较多；CEEMD通过添加噪声的方式，对模态混叠有一定程度的抑制作用，但是存在较多的伪分量，耗时最长；而MEEMD方法克服了上述方法的缺点，分解后的分量与实际的信号成分非常吻合，且强度损耗最小。三种方法的重构误差图如图6所示。

Figure 6.  Reconstruction error of three decomposition methods

可见，MEEMD方法对无序的噪声信号有较强的抑制作用，与EMD和CEEMD方法相比，重构误差最小，对振动信号特征提取有一定的优势。与算法相关的各项指标如表2所示。其中完备性以重构误差的均方差（Mean Squared Error，MSE）表示，正交性（Index of Orthogonality，IO）以公式（17）表示^[19]，IO越接近0，表明各模态分量${C_j}\left( t \right)$、${C_k}\left( t \right)$越各自独立，模态之间的混叠越小。

由表2中的各个参数可知，由于EMD方法没有添加噪声信号抑制模态混叠，运算时间最短，但信号重构误差最大，分解准确性最差；MEEMD方法时间相对CEEMD方法时间较短，IMF个数少，几乎没有伪分量的存在，完备性和正交性最优。

采用如图1所示的分布式光纤振动传感系统，传感光纤长度为6000 m，在4985 m处接入一个压电陶瓷PZT模拟振动源，分别采用已知的、不同频率的单频和混频的信号驱动PZT，以验证MEEMD-HHT方法用于振动信号特征提取的有效性。振动事件的频率范围一般集中在4~512 Hz之间^[20]，实验中选取100 Hz和200 Hz两组单频信号以及50、200、350 Hz和150、250、450 Hz两组混频信号，分别驱动PZT以模拟振动源。振动信号经相位解调后的时间-相位变化如图7(a)~(d)所示。

Figure 7.  Phase-time change curves of vibration signals of different frequencies demodulated by the system

对解调所得上述四种时间-相位信号分别进行EMD、CEEMD和MEEMD分解，并对各自获得的IMF分量分别进行Hilbert变换，得到的相对应的Hilbert谱。其中100 Hz单频信号分别经三种分解方法的Hilbert谱图如图8(a)~(c)所示。由图8（a）可知，EMD-HHT方法得到的Hilbert谱图在100 Hz处频率随时间变化不稳定，特别是在0.02~0.06 s之间极易受其他率分量的影响，频率成分复杂，使其出现模态混叠问题，且含有较多的伪分量；由图8（b）可知，CEEMD-HHT方法得到的Hilbert谱图在100 Hz频率处得到了明显改善，但是在100 Hz附近位置处仍存在伪分量问题；由图8（c）可知，MEEMD-HHT方法得到的Hilbert谱图频率随时间稳定分布在100 Hz处，有效避免了EMD和CEEMD方法出现的模态混叠和伪分量问题。

Figure 8.  Hilbert spectrum after three kinds of empirical mode decomposition and Hilbert transforn of vibration signal whose frequency is 100 Hz

对于200 Hz单频信号，同样经三种分解方法获得对应的Hilbert谱图，如图9(a)~(c)所示。由图9（a）可知，EMD-HHT方法得到的Hilbert谱图在200 Hz频率处随时间变化极不稳定，特别是在0.02 s、0.04 s、0.06 s、0.08 s和0.1 s附近处频率变化很大，且含有多个频率成分；由图9（b）可知，CEEMD-HHT方法得到的Hilbert谱图效果明显改善，但仍在200 Hz附近存在起伏变化；由图9（c）可知，MEEMD-HHT方法得到的Hilbert谱图频率稳定分布在200 Hz位置处。

Figure 9.  Hilbert spectrum after three kinds of empirical mode decomposition and Hilbert transforn of vibration signal whose frequency is 200 Hz

同样地，经相位解调后的50、200、350 Hz混频时间-相位信号，分别求其EMD、CEEMD和MEEMD方法的Hilbert谱图，如图10(a)~(c)所示，其中①、②、③分别代表三种分解方法提取的50、200、350 Hz频率特征。由图10（a）可知，EMD-HHT方法提取的三种频率特征极其杂乱，在高频350 Hz处尤为明显，在200 Hz处也存在较大的起伏，在0~0.01 s、0.04~0.05 s处未提取到50 Hz的特征；由图10（b）可知，CEEMD-HHT方法得到的Hilbert谱图在高频350 Hz处杂乱程度得到了改善，但是仍存在较大的起伏，不能正确提取频率特征；由图10（c）可知，MEEMD-HHT方法得到的Hilbert谱图可稳定、清晰地提取的三种频率特征，且无伪分量产生。

Figure 10.  Hilbert spectrum after three kinds of empirical mode decomposition and Hilbert transforn of vibration signal whose frequency is 50 Hz，200 Hz，350 Hz

150、250、450 Hz的混频时间-相位信号，经三种分解方法得到的Hilbert谱图分别如图11(a)~(c)所示。其中①、②、③分别代表三种分解方法提取的150 、250、450 Hz频率特征。由图11（a）可知，EMD-HHT提取的150 Hz频率特征随时间存在较大起伏，而250 Hz的特征并未提取到，对于450 Hz频率特征随时间变化更为杂乱，模态混叠现象严重；由图11（b）可知，CEEMD-HHT方法虽然改善了150 Hz特征的起伏程度，但是仍然没有提取到250 Hz特征，且高频450 Hz特征随时间变化程度较大；由图11（c）可知，MEEMD-HHT方法可以准确、稳定地获取150、250、450 Hz三种频率特征。

Figure 11.  Hilbert spectrum after three kinds of empirical mode decomposition and Hilbert transforn of vibration signal whose frequency is 150 Hz，250 Hz，450 Hz

由以上结果可以看出，无论是对单频信号还是混频信号的分析，相对于其他两种方法，MEEMD方法所对应的Hilbert谱中振动信号频率所对应的能量最为集中，无模态混叠和伪分量问题，有用特征最为清晰。若在Hilbert图某一频率上存在着能量，意味着振动信号一定含有该频率分量，而该振动出现的具体时刻已在Hilbert谱中给出，将上述Hilbert谱在0.1 s时间上进行积分，使之从时间-频率-幅值三者间的关系转变为频率-幅值两者间的关系，即为边际谱，它可以更精确地反映信号的实际频率成分。四种信号求得的边际谱如图12(a)~(d)所示。

Figure 12.  Marginal spectra of three decomposition methods on different frequencies

由图12（a）可知，EMD方法求得的边际谱峰值频率点为99 Hz，存在特征提取偏差，且在0~50 Hz范围内仍存在伪分量频率成分；CEEMD方法求得的边际谱虽然得到了正确的100 Hz频率特征，但是在100 Hz附近也存在伪分量频率；而MEEMD方法能够准确提取到100 Hz频率特征，且幅值最高，分布更为集中。

由图12（b）可知，EMD方法求得的边际谱峰值频率点为201 Hz，与正确频率相差1 Hz，且在0~200 Hz范围内存在虚假的频率成分；CEEMD方法提取到的频率为198 Hz，同样存在误差；而MEEMD方法准确提取到200 Hz的频率特征。

由图12（c）可知，与EMD、CEEMD方法相比，MEEMD方法准确提取到50、200、350 Hz的频率特征，且分布高度集中。

由图12（d）可知，MEEMD方法同样可以准确、集中地提取到150、250 、450 Hz的三种频率特征，相对于EMD和CEEMD方法，频率混叠和伪分量成分占比非常小，且有用的特征频率幅值较高，获得的振动信号的特征量最为准确。

为进一步显示MEEMD方法的准确性，分别选取10、50、100、150、200、300、400、500 Hz信号各50组单频信号，依次采用三种方法进行频率特征提取的测试和实验，频率检测在±0.05 Hz范围内视为正确，检测结果如表3所示。

由表3可以看出，三种方法均有被误判的情况发生，而MEEMD对单一频率的振动信号有较好的提取准确性，随着样本的数量增加，准确率会更高，适用于分布式光纤振动传感系统。此外，对于20组包含三种不同频率的混频信号进行特征提取测试，频率检测在±2 Hz范围内视为正确，对比结果如表4所示。

由表4可以看出，MEEMD特征提取方法在获得相对较少的特征向量个数的情况下，保证对混频信号较高的识别准确率，且Hilbert变换提取特征向量的过程中耗时最短，随着频率成分的增加，该方法相对于其他经验模态分解方法的各性能指标的优势将更加明显，能够满足分布式光纤振动传感系统模式识别过程中对于特征提取准确率、实时性的要求。

文中对MEEMD-HHT方法原理及其执行过程中关键参数的选取进行了详细的理论推导；通过仿真模拟，验证了MEEMD方法相对于EMD、CEEMD方法在完备性、正交性、算法耗时等关键性能上的优势；以单频和混频信号驱动PZT，模拟外界的振动信号，对分布式光纤振动传感系统解调的相位-时间变化信号采用三种经验模态分解方法进行分解，结合Hilbert变换获得Hilbert谱和边际谱，通过对比，MEEMD-HHT获得的振动信号的特征最为接近真实的仿真信号，对于单频振动信号MEEMD的平均频率识别准确率最高，达99.2%，耗时居中；对于混频振动信号MEEMD方法的平均频率识别准确率达98.1%，相对于CEEMD方法准确率提高7%，相对于EMD方法提高15.6%，耗时最短，可有效抑制模态混叠和伪分量问题，对光纤后向散射中的随机噪声序列以及环境干扰有抑制作用，提高了实际应用中的振动信号特征提取效果，对分布式光纤振动传感系统的模式识别功能的准确、高效实现具有重大意义。在今后的研究工作中，可通过对实验样本数据的大量积累对该方法性能指标与信号频率成分的相关性给出定性评价，并对方法执行过程中时间序列重构的嵌入维数参数作出自适应判定。