基于CEEMD与改进Prony的谐波检测方法*

吴新忠，邢强，周涛，成江洋

（中国矿业大学信息与电气工程学院，江苏徐州221008）

0 引 言

近年来，随着国民经济的快速发展，大大推进了电网建设的步伐。但由于越来越多电力电子设备的投入，在电网中产生了大量的谐波与间谐波。由于谐波对电力系统和用电设备的安全、稳定和可靠运行有着重要影响，因此对谐波快速和准确的检测，为谐波污染和治理提供了前提和保证［1-2］。

目前，谐波检测作为电能质量研究的重心和出发点，现有的谐波检测和分析方法主要包括基于傅里叶分析法［3］、小波变换法［4］、粒子群算法［5］、S变换［6］算法等等。傅里叶变换检测精度高、功能多、实现简单且使用方便，是电网谐波检测应用最为普遍的一种方法，但存在频率检测分辨率低、在非同步采用条件下容易造成频谱泄露和栅栏效应的缺点。小波变换比到FFT具有优良的时频特性，采用小波检测法可以有效的检测出谐波与间谐波分量，但小波算法对噪声比较敏感，频率分辨率低且运算量较大。而基于粒子群算法和S变换算法则存在着计算量大、实时性不强的缺点。

文献［7］将Prony分析法应用到谐波检测中，Prony方法只需求解两组齐次线性方程和一个线性多项式可以简单有效地检测出谐波和间谐波的幅值、频率和相位。但传统Prony方法易受噪声影响，且模型阶数的选取和辨识准确率具有密切的关系。阶数选取过小，拟合程度达不到理想效果，而阶数选取较高，虽然可以增加拟合度但其计算量也随之增大。

而在处理非线性非平稳谐波信号时，HHT方法具有自适应性和抗干扰能力［8］。该方法首先把原始信号通过经验模态分解（Empirical Mode Decomposition，EMD）分解为不同频率的固有模态函数（Intrinsic Mode Function，IMF），然后根据IMF的瞬时变化特征设定阈值进行时空滤波和平稳化处理，有效地去除高频噪声分量和虚假分量，最后对分解的固有模态函数进行希尔伯特变换（Hilbert-Transformation），检测分析出谐波信号的瞬时特征参数。但EMD分解方法存在一定的模态混叠和端点效应等缺点，对瞬时幅值和频率的精确提取造成影响，且直接进行Hilbert变换检测幅值和相位结果往往不准确。

针对上述问题，文章结合EMD分解和Prony方法在谐波检测时的各自特点，提出使用完全经验模态分解（Complete Ensemble Empirical Mode Decomposition，CEEMD）与改进的Prony算法相结合实现对谐波分量进行辨识的新方法。CEEMD是由Torres［9］等人在EMD和总体平均经验模态分解（Ensemble Empirical Mode Decomposition，EEMD）［10］基础上改进提出的，对原始待分解信号加入正、负成对形式的辅助噪声，然后再对叠加信号求其平均值，得到更为彻底的分解结果，因此EMD存在的模态混叠现象能够得到进一步缓解，减少了虚假分量。因此，首先采用CEEMD将谐波信号进行分解得到各个模态函数IMF，然后利用能量门限法去除虚假分量和高频噪声分量提取出真实的谐波分量。再改进Prony算法系统阶数估计的过程，满足复杂谐波信号辨识精度要求。最后将真实的谐波成分进行改进Prony变换，分析出信号的幅值、频率和相位参数。

1 基本原理

1.1 CEEMD原理

CEEMD和EEMD都是通过加入辅助噪声达到分解目的。其中，CEEMD第一个IMF的求取方法与EEMD相同。定义运算符号Ej（·）为EMD分解的第j个固有模态函数，wi（t）是单位方差为零均值的高斯白噪声，εk是每个阶段设定的信噪比系数，如果x（t）是原始信号，则CEEMD算法的具体步骤如下：

（1）首先，对原始信号 x（t）加入高斯白噪声wi（t），对目标信号 x（t）＋ε0wi（t）进行 n次 EMD分解，获得第一个固有模态函数c1。

（2）然后，计算得到一阶残差 r1（t）。

（3）对 r1（t）＋ε1E1（wi（t））继续分解，i＝1，2...n，将它的第一个固有模态函数分量作为CEEMD的c2。

（4）以此类推计算第 k个剩余残量 rk，k＝2，3...K。

（5）再继续分解 rk（t）＋εkEk（wi（t）），i＝1，2...n，将它的第一个固有模态函数分量作为CEEMD的 c（k＋1）。

（6）继续分解直到剩余残差满足结束条件（残差的极值最多不超过两个），否则返回到步骤（4）～（6）中进行计算，最终获得的剩余残差R（t），其中k是IMF的总数。

式（7）表明CEEMD能够完全分解，并且可以获得精确的重构信号。

1.2 能量门限法

在CEEMD分解过程中，可以将信号分解为噪声分量、真实谐波分量和虚假分量，为了将这三种信号分量检测区分出来，可以采用能量门限法［11］分辨分量成分。能量门限法的核心思想：每个待分解的信号都由类型不同分量组成，而不同类型的信号分量能量等级一般不同，经过CEEMD分解得到的IMF也应该有相同等级的能量Mi，相对较大的主导IMF可判断为谐波分量，较小的IMF可以认为是虚假分量和噪声分量。设：

进行归一化可知：

对ei设置阈值，去除噪声分量和虚假分量，从而得到真实的谐波分量。

1.3 改进Prony方法原理

Prony方法可以简单有效的提取出信号的特征参数，通过降阶模型拟合系统的实测数据，因此在数字信号处理领域中有着广泛应用。但经典的Prony算法在系统阶数辨识和线性预测参数求解效果［12］并不理想，影响算法对非平稳信号拟合的精度。虽然通过提高Prony算法的阶数可以提高检测精度，但阶数越高往往会增加运算的计算量。因此，为了提高Prony算法检测的精度和拟合效果，本文对经典Prony方法进行改进并应用在谐波检测分析中。

Prony算法假定的数据模型是一组P个具有任意幅值、相位、频率和衰减因子的指数函数：设观测数据为x（n），运用Prony方法拟合离散时间函数。

改进的Prony方法主要步骤如下：

（1）构造样本矩阵

由输入信号序列计算样本函数r（i，j），并构造样本矩阵Rc，定义样本函数：

（2）系统阶数估计改进

线性参数估计可以看作求解方程组（12）的过程。

对样本矩阵Rc进行奇异值分解，可得到它的奇异值分布如式（13）所示。

实际中由于存在噪声，使样本矩阵中p-M维零空间被噪声空间取代。定义w（i）为系统阶数估计参数：

式中p为系统的估计阶数，由于w（i）是单调递增的，当i值从1向p递增时，w（i）的值会向1逼近，而信号空间的奇异值大于噪声空间的奇异值，因此当取到某值使得 w（i）大于限值（一般取 λ＝0.995）时，可认为此时的即为系统的实际阶数M。

（3）线性预测参数优化

确定系统的实际阶数M后，将系统的噪声空间以零空间替代，即得到样本矩阵Rc的最优近似矩阵。去除噪声空间影响后，式（12）中的参数矩阵Q＝［a1a2…ap］只有M个独立参数，则可构造p＋1-M维的方程组。

S（M）是由酉矩阵子向量组成的，则必存在逆矩阵S－（M），可得求取预测参数的解的表达式为：

（4）特征参数求解

（a）由式（18）求取预测参数a1a2…ap的估计值后，进一步代入式（19）中通过多项式求根得zi；

（b）然后将 zi代入式（10）中求得 bi；

（c）最后将 zi和 bi代入式（20）计算幅值 Ai、频率fi和初相位θi，完成谐波分量特征参数提取：

2 谐波检测方法流程

传统的EMD分解与Prony算法相结合［13］的分析方法提取特征参数，EMD直接分解模态混叠现象较为严重且拟合效果差无法获得精确的检测结果。因此本文采用CEEMD与改进Prony相结合的新方法对谐波信号进行分析，其方法基本步骤为：

（1）对原始含噪信号进行分解时，用CEEMD取代EMD进行分解，将信号从高频频率分量到低频频率分量依次抽取得到一系列固有模态函数ck（t）；

（2）对（1）中所得ck（t）代入式（8）和式（9）得到各分量能量，设定门限阈值ei，运用能量门限法去除虚假分量和噪声分量得到真实的谐波分量；

（3）对经典的Prony算法实际阶数和线性预测参数的求解过程进行综合改进，使其应用在谐波检测分析中，提高算法的辨识精度和拟合效果；

（4）最后，把真实的谐波分量进行改进Prony变换，估计出信号的频率、幅值和相位。

3 仿真信号分析

在实际电网中，电网信号既包含谐波和间谐波成分还有随机噪声扰动信号，为了进行仿真验证，通过公（21）构造真实的电网谐波信号：

式中N＝4；f1＝25；f2＝150；f3＝155；f4＝250；A1＝8.7；A2＝14.2；A3＝4.5；A4＝1.5；θ1＝π/3；θ2＝π/3；θ3＝π/4；θ4＝π/5；n（t）为随机加入的高斯白噪声。

采用Matlab在信噪比为20 dB的随机噪声环境下生成谐波信号，如图1所示。其中采样频率为8 000 Hz，采样长度为0.25 s，共采集2 000个数据点。

图1 原始谐波信号Fig.1 Original harmonic signal

3.1 去噪性能分析

为了提高本文方法检测时的精度和拟合效果，减少噪声对Prony方法特征提取的干扰，对原始信号应采取滤波消噪处理。因此，首先对原始信号分别进行EMD和CEEMD分解，其中CEEMD添加幅值为0.15的高斯白噪声，进行120次分解运算，分解结果如图2（a）和2（b）所示。

对比图2（a）和2（b）可以看出，在噪声影响的情况下，EMD分解可以得到8条IMF，CEEMD为9条。虽然EMD能将信号分解为不同频率的模态函数，但模态混叠现象严重，C4～C5模态分量由于混叠了不同频率成分存在较大扰动。而由于CEEMD加入的辅助噪声采用正、负成对的形式，分解效果最为彻底，受噪声干扰较小，能够有效的分解出各个分量，所以选取CEEMD对信号进行分解有效的解决了模态混叠现象。

图2 信号分解结果对比Fig.2 Results contrast of signal decomposition

为了去除噪声分量和虚假分量，提取真实的谐波成分，本文设置能量门限阈值ei＝0.1，根据式（8）和式（9）计算CEEMD分解出的各 Ci的能量 Mi和ei，若大于0.1则判断为谐波成分，小于0.1则为虚假分量和噪声分量。由表1可知，C3～C6为谐波分量，其余为噪声分量和虚假分量，予以舍弃。

CEEMD与能量门限法相结合既可以去除虚假分量，又可以起到平稳消噪的效果。为了定量检验CEEMD的去噪效果，将真实的谐波分量进行叠加重构［14］，与原始不含噪的理想信号对比。为了从宏观上衡量去噪重构的效果，选择重构信噪比（SNR）和均方误差百分值（MSE）评价指标，信噪比越大，幅值误差越小，说明去噪效果越好；微观衡量标准选择能量恢复系数（ERP）和波形畸变率（JBL），能量恢复系数越大，波形畸变率越小，说明波形畸变的越小去噪效果越好。

表1 固有模态函数能量值Tab.1 Energy values of intrinsic mode functions

并采用小波阈值去噪算法和基于基于3δ准则的EMD的阈值去噪算法与本文方法对比。对于小波算法，选用的母小波是db2，分解层数为6层。表2给出了不同去噪算法在四个评价指标的对比。

表2 不同算法的评价指标值Tab.2 Evaluation indicator values for different algorithms

对于本文去噪来说，虽然SNR比EMD的阈值去噪低1.13%，但是MSE小47.54%，ERP高6.69%，JBL低36.9%。综合评价可知，本文提出的方法要优于基于EMD的阈值去噪算法，而本文方法和db2小波去噪算法进行比较，它的各项指标都要优于db2。

3.2 检测性能分析

其后，为了验证改进方法前后的拟合效果，现使用原始Prony方法和改进的Prony方法对真实的谐波分量C3～C6进行拟合，实验结果如图3（a）和3（b）所示。可见改进Prony方法比到原始Prony方法对信号分量的拟合度更高，改善效果更为明显。

图3定性比较了Prony方法改进前后的拟合波形，为了更客观的定量评价拟合效果，图4（a）和图4（b）给出了两种方法均方误差对比图，原始Prony算法均方误差达到103.5，而改进Prony方法均方误差为34.6，可见改进Prony方法在误差意义上优于传统Prony方法。

图3 Prony方法拟合效果图Fig.3 Fitting effect diagram of Pronymethod

图4 Prony方法均方误差对比Fig.4 Mean square error contrast of Pronymethod

最后通过式（20）求得各谐波分量的幅值、频率和相位，完成特征参数的检测工作。表3～表5给出了用本文方法在信噪比为20 dB，实际阶数M＝8时，检测到谐波的特征参数，同时为了比较，分别用文献［3］中的FFT方法和文献［12］中的Prony方法对信号进行检测。从表中结果可以看出，FFT和Prony方法幅值检测误差为2.88%和1.06%；频率误差分别为3.30%和2.03%；相位误差分别为4.32%和2.73%；而本文方法的平均误差为0.27%、1.16%和1.51%。显然本文方法对谐波辨识精度要优于FFT与Prony方法，且对幅值检测精度最高，对频率和相位检测精度稍逊，对谐波的检测精度要高于间谐波。

表3 幅值检测结果Tab.3 Result of amplitudes detection

表4 频率检测结果Tab.4 Result of frequencies detection

表5 相位检测结果Tab.5 Result of phases detection

3.3 抗噪性能分析

抗噪性能是检验算法稳定性的重要标准，为了验证本文提出算法的抗噪性能，在不同的信噪比条件下，对频率为155 Hz谐波分量进行检测分析，根据均方根误差（RMSE）式（22）做了30组对比试验得到误差曲线如图5所示。

图5表明在不同的信噪比条件下，本文所提方法均方误差较低，明显优于其他2种方法。随着信噪比的提高，均方误差降低明显，能够准确的检测谐波频率分量，可见本文所提方法具有良好的抗噪性能和稳定性。

4 实测数据验证

为了更好的说明本文提出的谐波检测方法，在仿真试验的基础上，再通过对谐波源为三相不可控负载的实测数据进行验证。

图5 各检测方法不同信噪比下的均方根误差Fig.5 RMSE of three detection methods in SNRs

带负载三相不可控整流器参数设置：电容为2 200μF，电感 1.5 mH，负载电阻 8Ω。利用Fluke435电能质量检测仪对整流器负载电流进行测试，Fluke采样频率为10 kHz，检测谐波电流波形、频谱含量如图6（a）、图6（b）所示。取L1相做研究，并将实测数据导入MATLAB通过本文所提方法对实测电流数据进行分析，检测出谐波成分与含量如表6所示，通过和实际仪器对比分析，验证本文方法对谐波辨识精度较高。

图6 负载实验检测结果Fig.6 Test results of load experiment

表6 负载实验结果对比Tab.6 Comparison of load experiment results

5 结束语

文章在EMD和Prony算法的基础上，对传统的Prony方法进行改进，并与CEEMD相结合引入到谐波检测领域中。通过对Prony方法实际阶数和线性预测参数求解过程进行改进，实现对谐波分量频率、幅值和相位特征参数的有效提取。以MATLAB仿真实验为平台，并与FFT和Prony谐波检测方法对比验证本文所提检测方法的高精确度和良好的抗噪性能。最后利用三相不可控负载的实测数据，通过和电能质量分析仪检测结果分析对比，结果证明本文所提方法对谐波检测准确率较高，为谐波检测与分析提供了一种新的思路和方法。