基于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諧波檢測方法對比驗證本文所提檢測方法的高精確度和良好的抗噪性能。最后利用三相不可控負載的實測數據，通過和電能質量分析儀檢測結果分析對比，結果證明本文所提方法對諧波檢測準確率較高，為諧波檢測與分析提供了一種新的思路和方法。