吳新忠，邢強(qiáng)，周濤，成江洋

（中國(guó)礦業(yè)大學(xué)信息與電氣工程學(xué)院，江蘇徐州221008）

0 引 言

近年來(lái)，隨著國(guó)民經(jīng)濟(jì)的快速發(fā)展，大大推進(jìn)了電網(wǎng)建設(shè)的步伐。但由于越來(lái)越多電力電子設(shè)備的投入，在電網(wǎng)中產(chǎn)生了大量的諧波與間諧波。由于諧波對(duì)電力系統(tǒng)和用電設(shè)備的安全、穩(wěn)定和可靠運(yùn)行有著重要影響，因此對(duì)諧波快速和準(zhǔn)確的檢測(cè)，為諧波污染和治理提供了前提和保證［1-2］。

目前，諧波檢測(cè)作為電能質(zhì)量研究的重心和出發(fā)點(diǎn)，現(xiàn)有的諧波檢測(cè)和分析方法主要包括基于傅里葉分析法［3］、小波變換法［4］、粒子群算法［5］、S變換［6］算法等等。傅里葉變換檢測(cè)精度高、功能多、實(shí)現(xiàn)簡(jiǎn)單且使用方便，是電網(wǎng)諧波檢測(cè)應(yīng)用最為普遍的一種方法，但存在頻率檢測(cè)分辨率低、在非同步采用條件下容易造成頻譜泄露和柵欄效應(yīng)的缺點(diǎn)。小波變換比到FFT具有優(yōu)良的時(shí)頻特性，采用小波檢測(cè)法可以有效的檢測(cè)出諧波與間諧波分量，但小波算法對(duì)噪聲比較敏感，頻率分辨率低且運(yùn)算量較大。而基于粒子群算法和S變換算法則存在著計(jì)算量大、實(shí)時(shí)性不強(qiáng)的缺點(diǎn)。

文獻(xiàn)［7］將Prony分析法應(yīng)用到諧波檢測(cè)中，Prony方法只需求解兩組齊次線性方程和一個(gè)線性多項(xiàng)式可以簡(jiǎn)單有效地檢測(cè)出諧波和間諧波的幅值、頻率和相位。但傳統(tǒng)Prony方法易受噪聲影響，且模型階數(shù)的選取和辨識(shí)準(zhǔn)確率具有密切的關(guān)系。階數(shù)選取過(guò)小，擬合程度達(dá)不到理想效果，而階數(shù)選取較高，雖然可以增加擬合度但其計(jì)算量也隨之增大。

而在處理非線性非平穩(wěn)諧波信號(hào)時(shí)，HHT方法具有自適應(yīng)性和抗干擾能力［8］。該方法首先把原始信號(hào)通過(guò)經(jīng)驗(yàn)?zāi)B(tài)分解（Empirical Mode Decomposition，EMD）分解為不同頻率的固有模態(tài)函數(shù)（Intrinsic Mode Function，IMF），然后根據(jù)IMF的瞬時(shí)變化特征設(shè)定閾值進(jìn)行時(shí)空濾波和平穩(wěn)化處理，有效地去除高頻噪聲分量和虛假分量，最后對(duì)分解的固有模態(tài)函數(shù)進(jìn)行希爾伯特變換（Hilbert-Transformation），檢測(cè)分析出諧波信號(hào)的瞬時(shí)特征參數(shù)。但EMD分解方法存在一定的模態(tài)混疊和端點(diǎn)效應(yīng)等缺點(diǎn)，對(duì)瞬時(shí)幅值和頻率的精確提取造成影響，且直接進(jìn)行Hilbert變換檢測(cè)幅值和相位結(jié)果往往不準(zhǔn)確。

針對(duì)上述問(wèn)題，文章結(jié)合EMD分解和Prony方法在諧波檢測(cè)時(shí)的各自特點(diǎn)，提出使用完全經(jīng)驗(yàn)?zāi)B(tài)分解（Complete Ensemble Empirical Mode Decomposition，CEEMD）與改進(jìn)的Prony算法相結(jié)合實(shí)現(xiàn)對(duì)諧波分量進(jìn)行辨識(shí)的新方法。CEEMD是由Torres［9］等人在EMD和總體平均經(jīng)驗(yàn)?zāi)B(tài)分解（Ensemble Empirical Mode Decomposition，EEMD）［10］基礎(chǔ)上改進(jìn)提出的，對(duì)原始待分解信號(hào)加入正、負(fù)成對(duì)形式的輔助噪聲，然后再對(duì)疊加信號(hào)求其平均值，得到更為徹底的分解結(jié)果，因此EMD存在的模態(tài)混疊現(xiàn)象能夠得到進(jìn)一步緩解，減少了虛假分量。因此，首先采用CEEMD將諧波信號(hào)進(jìn)行分解得到各個(gè)模態(tài)函數(shù)IMF，然后利用能量門限法去除虛假分量和高頻噪聲分量提取出真實(shí)的諧波分量。再改進(jìn)Prony算法系統(tǒng)階數(shù)估計(jì)的過(guò)程，滿足復(fù)雜諧波信號(hào)辨識(shí)精度要求。最后將真實(shí)的諧波成分進(jìn)行改進(jìn)Prony變換，分析出信號(hào)的幅值、頻率和相位參數(shù)。

1 基本原理

1.1 CEEMD原理

CEEMD和EEMD都是通過(guò)加入輔助噪聲達(dá)到分解目的。其中，CEEMD第一個(gè)IMF的求取方法與EEMD相同。定義運(yùn)算符號(hào)Ej（·）為EMD分解的第j個(gè)固有模態(tài)函數(shù)，wi（t）是單位方差為零均值的高斯白噪聲，εk是每個(gè)階段設(shè)定的信噪比系數(shù)，如果x（t）是原始信號(hào)，則CEEMD算法的具體步驟如下：

（1）首先，對(duì)原始信號(hào) x（t）加入高斯白噪聲wi（t），對(duì)目標(biāo)信號(hào) x（t）＋ε0wi（t）進(jìn)行 n次 EMD分解，獲得第一個(gè)固有模態(tài)函數(shù)c1。

（2）然后，計(jì)算得到一階殘差 r1（t）。

（3）對(duì) r1（t）＋ε1E1（wi（t））繼續(xù)分解，i＝1，2...n，將它的第一個(gè)固有模態(tài)函數(shù)分量作為CEEMD的c2。

（4）以此類推計(jì)算第 k個(gè)剩余殘量 rk，k＝2，3...K。

（5）再繼續(xù)分解 rk（t）＋εkEk（wi（t）），i＝1，2...n，將它的第一個(gè)固有模態(tài)函數(shù)分量作為CEEMD的 c（k＋1）。

（6）繼續(xù)分解直到剩余殘差滿足結(jié)束條件（殘差的極值最多不超過(guò)兩個(gè)），否則返回到步驟（4）～（6）中進(jìn)行計(jì)算，最終獲得的剩余殘差R（t），其中k是IMF的總數(shù)。

式（7）表明CEEMD能夠完全分解，并且可以獲得精確的重構(gòu)信號(hào)。

1.2 能量門限法

在CEEMD分解過(guò)程中，可以將信號(hào)分解為噪聲分量、真實(shí)諧波分量和虛假分量，為了將這三種信號(hào)分量檢測(cè)區(qū)分出來(lái)，可以采用能量門限法［11］分辨分量成分。能量門限法的核心思想：每個(gè)待分解的信號(hào)都由類型不同分量組成，而不同類型的信號(hào)分量能量等級(jí)一般不同，經(jīng)過(guò)CEEMD分解得到的IMF也應(yīng)該有相同等級(jí)的能量Mi，相對(duì)較大的主導(dǎo)IMF可判斷為諧波分量，較小的IMF可以認(rèn)為是虛假分量和噪聲分量。設(shè)：

進(jìn)行歸一化可知：

對(duì)ei設(shè)置閾值，去除噪聲分量和虛假分量，從而得到真實(shí)的諧波分量。

1.3 改進(jìn)Prony方法原理

Prony方法可以簡(jiǎn)單有效的提取出信號(hào)的特征參數(shù)，通過(guò)降階模型擬合系統(tǒng)的實(shí)測(cè)數(shù)據(jù)，因此在數(shù)字信號(hào)處理領(lǐng)域中有著廣泛應(yīng)用。但經(jīng)典的Prony算法在系統(tǒng)階數(shù)辨識(shí)和線性預(yù)測(cè)參數(shù)求解效果［12］并不理想，影響算法對(duì)非平穩(wěn)信號(hào)擬合的精度。雖然通過(guò)提高Prony算法的階數(shù)可以提高檢測(cè)精度，但階數(shù)越高往往會(huì)增加運(yùn)算的計(jì)算量。因此，為了提高Prony算法檢測(cè)的精度和擬合效果，本文對(duì)經(jīng)典Prony方法進(jìn)行改進(jìn)并應(yīng)用在諧波檢測(cè)分析中。

Prony算法假定的數(shù)據(jù)模型是一組P個(gè)具有任意幅值、相位、頻率和衰減因子的指數(shù)函數(shù)：設(shè)觀測(cè)數(shù)據(jù)為x（n），運(yùn)用Prony方法擬合離散時(shí)間函數(shù)。

改進(jìn)的Prony方法主要步驟如下：

（1）構(gòu)造樣本矩陣

由輸入信號(hào)序列計(jì)算樣本函數(shù)r（i，j），并構(gòu)造樣本矩陣Rc，定義樣本函數(shù)：

（2）系統(tǒng)階數(shù)估計(jì)改進(jìn)

線性參數(shù)估計(jì)可以看作求解方程組（12）的過(guò)程。

對(duì)樣本矩陣Rc進(jìn)行奇異值分解，可得到它的奇異值分布如式（13）所示。

實(shí)際中由于存在噪聲，使樣本矩陣中p-M維零空間被噪聲空間取代。定義w（i）為系統(tǒng)階數(shù)估計(jì)參數(shù)：

式中p為系統(tǒng)的估計(jì)階數(shù)，由于w（i）是單調(diào)遞增的，當(dāng)i值從1向p遞增時(shí)，w（i）的值會(huì)向1逼近，而信號(hào)空間的奇異值大于噪聲空間的奇異值，因此當(dāng)取到某值使得 w（i）大于限值（一般取 λ＝0.995）時(shí)，可認(rèn)為此時(shí)的即為系統(tǒng)的實(shí)際階數(shù)M。

（3）線性預(yù)測(cè)參數(shù)優(yōu)化

確定系統(tǒng)的實(shí)際階數(shù)M后，將系統(tǒng)的噪聲空間以零空間替代，即得到樣本矩陣Rc的最優(yōu)近似矩陣。去除噪聲空間影響后，式（12）中的參數(shù)矩陣Q＝［a1a2…ap］只有M個(gè)獨(dú)立參數(shù)，則可構(gòu)造p＋1-M維的方程組。

S（M）是由酉矩陣子向量組成的，則必存在逆矩陣S－（M），可得求取預(yù)測(cè)參數(shù)的解的表達(dá)式為：

（4）特征參數(shù)求解

（a）由式（18）求取預(yù)測(cè)參數(shù)a1a2…ap的估計(jì)值后，進(jìn)一步代入式（19）中通過(guò)多項(xiàng)式求根得zi；

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

（c）最后將 zi和 bi代入式（20）計(jì)算幅值 Ai、頻率fi和初相位θi，完成諧波分量特征參數(shù)提取：

2 諧波檢測(cè)方法流程

傳統(tǒng)的EMD分解與Prony算法相結(jié)合［13］的分析方法提取特征參數(shù)，EMD直接分解模態(tài)混疊現(xiàn)象較為嚴(yán)重且擬合效果差無(wú)法獲得精確的檢測(cè)結(jié)果。因此本文采用CEEMD與改進(jìn)Prony相結(jié)合的新方法對(duì)諧波信號(hào)進(jìn)行分析，其方法基本步驟為：

（1）對(duì)原始含噪信號(hào)進(jìn)行分解時(shí)，用CEEMD取代EMD進(jìn)行分解，將信號(hào)從高頻頻率分量到低頻頻率分量依次抽取得到一系列固有模態(tài)函數(shù)ck（t）；

（2）對(duì)（1）中所得ck（t）代入式（8）和式（9）得到各分量能量，設(shè)定門限閾值ei，運(yùn)用能量門限法去除虛假分量和噪聲分量得到真實(shí)的諧波分量；

（3）對(duì)經(jīng)典的Prony算法實(shí)際階數(shù)和線性預(yù)測(cè)參數(shù)的求解過(guò)程進(jìn)行綜合改進(jìn)，使其應(yīng)用在諧波檢測(cè)分析中，提高算法的辨識(shí)精度和擬合效果；

（4）最后，把真實(shí)的諧波分量進(jìn)行改進(jìn)Prony變換，估計(jì)出信號(hào)的頻率、幅值和相位。

3 仿真信號(hào)分析

在實(shí)際電網(wǎng)中，電網(wǎng)信號(hào)既包含諧波和間諧波成分還有隨機(jī)噪聲擾動(dòng)信號(hào)，為了進(jìn)行仿真驗(yàn)證，通過(guò)公（21）構(gòu)造真實(shí)的電網(wǎng)諧波信號(hào)：

式中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）為隨機(jī)加入的高斯白噪聲。

采用Matlab在信噪比為20 dB的隨機(jī)噪聲環(huán)境下生成諧波信號(hào)，如圖1所示。其中采樣頻率為8 000 Hz，采樣長(zhǎng)度為0.25 s，共采集2 000個(gè)數(shù)據(jù)點(diǎn)。

圖1 原始諧波信號(hào)Fig.1 Original harmonic signal

3.1 去噪性能分析

為了提高本文方法檢測(cè)時(shí)的精度和擬合效果，減少噪聲對(duì)Prony方法特征提取的干擾，對(duì)原始信號(hào)應(yīng)采取濾波消噪處理。因此，首先對(duì)原始信號(hào)分別進(jìn)行EMD和CEEMD分解，其中CEEMD添加幅值為0.15的高斯白噪聲，進(jìn)行120次分解運(yùn)算，分解結(jié)果如圖2（a）和2（b）所示。

對(duì)比圖2（a）和2（b）可以看出，在噪聲影響的情況下，EMD分解可以得到8條IMF，CEEMD為9條。雖然EMD能將信號(hào)分解為不同頻率的模態(tài)函數(shù)，但模態(tài)混疊現(xiàn)象嚴(yán)重，C4～C5模態(tài)分量由于混疊了不同頻率成分存在較大擾動(dòng)。而由于CEEMD加入的輔助噪聲采用正、負(fù)成對(duì)的形式，分解效果最為徹底，受噪聲干擾較小，能夠有效的分解出各個(gè)分量，所以選取CEEMD對(duì)信號(hào)進(jìn)行分解有效的解決了模態(tài)混疊現(xiàn)象。

圖2 信號(hào)分解結(jié)果對(duì)比Fig.2 Results contrast of signal decomposition

為了去除噪聲分量和虛假分量，提取真實(shí)的諧波成分，本文設(shè)置能量門限閾值ei＝0.1，根據(jù)式（8）和式（9）計(jì)算CEEMD分解出的各 Ci的能量 Mi和ei，若大于0.1則判斷為諧波成分，小于0.1則為虛假分量和噪聲分量。由表1可知，C3～C6為諧波分量，其余為噪聲分量和虛假分量，予以舍棄。

CEEMD與能量門限法相結(jié)合既可以去除虛假分量，又可以起到平穩(wěn)消噪的效果。為了定量檢驗(yàn)CEEMD的去噪效果，將真實(shí)的諧波分量進(jìn)行疊加重構(gòu)［14］，與原始不含噪的理想信號(hào)對(duì)比。為了從宏觀上衡量去噪重構(gòu)的效果，選擇重構(gòu)信噪比（SNR）和均方誤差百分值（MSE）評(píng)價(jià)指標(biāo)，信噪比越大，幅值誤差越小，說(shuō)明去噪效果越好；微觀衡量標(biāo)準(zhǔn)選擇能量恢復(fù)系數(shù)（ERP）和波形畸變率（JBL），能量恢復(fù)系數(shù)越大，波形畸變率越小，說(shuō)明波形畸變的越小去噪效果越好。

表1 固有模態(tài)函數(shù)能量值Tab.1 Energy values of intrinsic mode functions

并采用小波閾值去噪算法和基于基于3δ準(zhǔn)則的EMD的閾值去噪算法與本文方法對(duì)比。對(duì)于小波算法，選用的母小波是db2，分解層數(shù)為6層。表2給出了不同去噪算法在四個(gè)評(píng)價(jià)指標(biāo)的對(duì)比。

表2 不同算法的評(píng)價(jià)指標(biāo)值Tab.2 Evaluation indicator values for different algorithms

對(duì)于本文去噪來(lái)說(shuō)，雖然SNR比EMD的閾值去噪低1.13%，但是MSE小47.54%，ERP高6.69%，JBL低36.9%。綜合評(píng)價(jià)可知，本文提出的方法要優(yōu)于基于EMD的閾值去噪算法，而本文方法和db2小波去噪算法進(jìn)行比較，它的各項(xiàng)指標(biāo)都要優(yōu)于db2。

3.2 檢測(cè)性能分析

其后，為了驗(yàn)證改進(jìn)方法前后的擬合效果，現(xiàn)使用原始Prony方法和改進(jìn)的Prony方法對(duì)真實(shí)的諧波分量C3～C6進(jìn)行擬合，實(shí)驗(yàn)結(jié)果如圖3（a）和3（b）所示?？梢?jiàn)改進(jìn)Prony方法比到原始Prony方法對(duì)信號(hào)分量的擬合度更高，改善效果更為明顯。

圖3定性比較了Prony方法改進(jìn)前后的擬合波形，為了更客觀的定量評(píng)價(jià)擬合效果，圖4（a）和圖4（b）給出了兩種方法均方誤差對(duì)比圖，原始Prony算法均方誤差達(dá)到103.5，而改進(jìn)Prony方法均方誤差為34.6，可見(jiàn)改進(jìn)Prony方法在誤差意義上優(yōu)于傳統(tǒng)Prony方法。

圖3 Prony方法擬合效果圖Fig.3 Fitting effect diagram of Pronymethod

圖4 Prony方法均方誤差對(duì)比Fig.4 Mean square error contrast of Pronymethod

最后通過(guò)式（20）求得各諧波分量的幅值、頻率和相位，完成特征參數(shù)的檢測(cè)工作。表3～表5給出了用本文方法在信噪比為20 dB，實(shí)際階數(shù)M＝8時(shí)，檢測(cè)到諧波的特征參數(shù)，同時(shí)為了比較，分別用文獻(xiàn)［3］中的FFT方法和文獻(xiàn)［12］中的Prony方法對(duì)信號(hào)進(jìn)行檢測(cè)。從表中結(jié)果可以看出，F(xiàn)FT和Prony方法幅值檢測(cè)誤差為2.88%和1.06%；頻率誤差分別為3.30%和2.03%；相位誤差分別為4.32%和2.73%；而本文方法的平均誤差為0.27%、1.16%和1.51%。顯然本文方法對(duì)諧波辨識(shí)精度要優(yōu)于FFT與Prony方法，且對(duì)幅值檢測(cè)精度最高，對(duì)頻率和相位檢測(cè)精度稍遜，對(duì)諧波的檢測(cè)精度要高于間諧波。

表3 幅值檢測(cè)結(jié)果Tab.3 Result of amplitudes detection

表4 頻率檢測(cè)結(jié)果Tab.4 Result of frequencies detection

表5 相位檢測(cè)結(jié)果Tab.5 Result of phases detection

3.3 抗噪性能分析

抗噪性能是檢驗(yàn)算法穩(wěn)定性的重要標(biāo)準(zhǔn)，為了驗(yàn)證本文提出算法的抗噪性能，在不同的信噪比條件下，對(duì)頻率為155 Hz諧波分量進(jìn)行檢測(cè)分析，根據(jù)均方根誤差（RMSE）式（22）做了30組對(duì)比試驗(yàn)得到誤差曲線如圖5所示。

圖5表明在不同的信噪比條件下，本文所提方法均方誤差較低，明顯優(yōu)于其他2種方法。隨著信噪比的提高，均方誤差降低明顯，能夠準(zhǔn)確的檢測(cè)諧波頻率分量，可見(jiàn)本文所提方法具有良好的抗噪性能和穩(wěn)定性。

4 實(shí)測(cè)數(shù)據(jù)驗(yàn)證

為了更好的說(shuō)明本文提出的諧波檢測(cè)方法，在仿真試驗(yàn)的基礎(chǔ)上，再通過(guò)對(duì)諧波源為三相不可控負(fù)載的實(shí)測(cè)數(shù)據(jù)進(jìn)行驗(yàn)證。

圖5 各檢測(cè)方法不同信噪比下的均方根誤差Fig.5 RMSE of three detection methods in SNRs

帶負(fù)載三相不可控整流器參數(shù)設(shè)置：電容為2 200μF，電感 1.5 mH，負(fù)載電阻 8Ω。利用Fluke435電能質(zhì)量檢測(cè)儀對(duì)整流器負(fù)載電流進(jìn)行測(cè)試，F(xiàn)luke采樣頻率為10 kHz，檢測(cè)諧波電流波形、頻譜含量如圖6（a）、圖6（b）所示。取L1相做研究，并將實(shí)測(cè)數(shù)據(jù)導(dǎo)入MATLAB通過(guò)本文所提方法對(duì)實(shí)測(cè)電流數(shù)據(jù)進(jìn)行分析，檢測(cè)出諧波成分與含量如表6所示，通過(guò)和實(shí)際儀器對(duì)比分析，驗(yàn)證本文方法對(duì)諧波辨識(shí)精度較高。

圖6 負(fù)載實(shí)驗(yàn)檢測(cè)結(jié)果Fig.6 Test results of load experiment

表6 負(fù)載實(shí)驗(yàn)結(jié)果對(duì)比Tab.6 Comparison of load experiment results

5 結(jié)束語(yǔ)

文章在EMD和Prony算法的基礎(chǔ)上，對(duì)傳統(tǒng)的Prony方法進(jìn)行改進(jìn)，并與CEEMD相結(jié)合引入到諧波檢測(cè)領(lǐng)域中。通過(guò)對(duì)Prony方法實(shí)際階數(shù)和線性預(yù)測(cè)參數(shù)求解過(guò)程進(jìn)行改進(jìn)，實(shí)現(xiàn)對(duì)諧波分量頻率、幅值和相位特征參數(shù)的有效提取。以MATLAB仿真實(shí)驗(yàn)為平臺(tái)，并與FFT和Prony諧波檢測(cè)方法對(duì)比驗(yàn)證本文所提檢測(cè)方法的高精確度和良好的抗噪性能。最后利用三相不可控負(fù)載的實(shí)測(cè)數(shù)據(jù)，通過(guò)和電能質(zhì)量分析儀檢測(cè)結(jié)果分析對(duì)比，結(jié)果證明本文所提方法對(duì)諧波檢測(cè)準(zhǔn)確率較高，為諧波檢測(cè)與分析提供了一種新的思路和方法。