楊玉坤，楊明玉

（1.吉林省電力勘測(cè)設(shè)計(jì)院，長(zhǎng)春 130022；2.華北電力大學(xué)電力系統(tǒng)保護(hù)與動(dòng)態(tài)安全監(jiān)控教育部重點(diǎn)實(shí)驗(yàn)室，保定 071003）

隨著電力電子裝置在電力系統(tǒng)中的廣泛應(yīng)用以及非線性負(fù)荷的日益增多，電力系統(tǒng)中的諧波污染越來越嚴(yán)重和復(fù)雜。諧波污染除了與基波成整數(shù)倍的諧波外，還存在許多非整數(shù)倍的間諧波。間諧波的存在不僅會(huì)對(duì)電力系統(tǒng)造成更大的危害，更增加了諧波檢測(cè)的難度。

目前常用的諧波檢測(cè)分析方法是基于傅里葉變換的算法，但是它們對(duì)間諧波的分析存在著頻譜泄漏和柵欄效應(yīng)現(xiàn)象［1］。近年來一些新的方法被用于諧波檢測(cè)當(dāng)中，其中現(xiàn)代譜估計(jì)方法中的Prony算法由于無需估計(jì)樣本自相關(guān)函數(shù)，只需求解兩組齊次線性方程和一個(gè)線性多項(xiàng)式便可以一次性辨識(shí)出信號(hào)中各分量的頻率、幅值、衰減因子和初相，并具有較高的辨識(shí)精度，該算法已在國(guó)內(nèi)外的諧波檢測(cè)研究中得到充分重視［2～4］。

以往Prony算法的諧波檢測(cè)研究中，都是依靠濾波或迭代的方法來提高Prony算法的辨識(shí)精度，但是它們都是以提升算法的計(jì)算量為代價(jià)的。而在線的諧波檢測(cè)要求具有很好的實(shí)時(shí)性，因此在保證精度的前提下，尋找一種能夠減小算法復(fù)雜度的方法具有十分重要的意義。另外，Prony算法的辨識(shí)精度與其參數(shù)選擇有很大的關(guān)系，而以往的研究中，沒有細(xì)致地討論過最優(yōu)參數(shù)的選取問題，對(duì)于如何識(shí)別由于較高模型階數(shù)引起的虛假頻率成分問題，也沒有一個(gè)可靠的方法。本文就針對(duì)該問題進(jìn)行研究。

1 Prony算法建模與求解的基本思路

Prony算法采用p個(gè)指數(shù)項(xiàng)的線性組合來擬合等間隔采樣數(shù)據(jù)（本文后面稱p為模型的階數(shù)）。設(shè)N個(gè)采樣數(shù)據(jù)為x（0），x（1），…，x（N－1），則有

定義

可由歐拉公式得

即只需求解bi和zi即可用p個(gè)具有任意幅值A(chǔ)i、頻率fi、相位θi和衰減因子αi的余弦分量來擬合采樣數(shù)據(jù)。式（3）是常系數(shù)線性差分方程

的齊次解其特征多項(xiàng)式為

并代入式（4）可得

Prony算法的求解即為求取ai使誤差的平方達(dá)到最小。將ai代入特征多項(xiàng)式（5）求得zi，再代入由式（1）組成的方程組求得bi，即

最后利用zi和bi的定義式求得p個(gè)余弦量的幅值A(chǔ)i、頻率fi、相位θi和衰減因子αi。

2 Prony算法求解方法的改進(jìn)

對(duì)于式（7）中ai的求解，通常的方法是構(gòu)造誤差ε（n）的平方和J（a）為代價(jià)函數(shù)，即

令

使其達(dá)到最小，最后得到方程組［5］為

本文提出的方法是直接令誤差ε（n）＝0，則由式（7）可得求解ai的方程組為

可簡(jiǎn)寫為Ya＝y(tǒng)

對(duì)于方程組（12），由于一般取p＜N／2通常為等式個(gè)數(shù)大于未知數(shù)個(gè)數(shù)的超定方程組。對(duì)于求解超定方程組的一般方法有最小二乘法、奇異值分解SVD（singular value decomposition）和正交三角分解QR（quadrature right－triangle）。SVD算法和QR算法由于避免了矩陣的求逆運(yùn)算，所以較最小二乘法穩(wěn)定，且有相應(yīng)的成熟算法，求解速度非?？?。文獻(xiàn)［7］比較了SVD分解和QR分解算法的浮點(diǎn)運(yùn)算數(shù)量級(jí)，當(dāng)方程組（12）中的Y是一個(gè)階矩陣時(shí)，Golub－Reinsch SVD算法的數(shù)量級(jí)為4mn2＋8n3，而Householder QR算法的僅為2mn2－2n3／3。為了盡量縮短Prony算法在線諧波檢測(cè)的耗時(shí)，本文采用QR分解算法來求解上述超定方程組。

對(duì)矩陣Y進(jìn)行QR分解，就是把Y分解為一個(gè)正交矩陣Q和一個(gè)上三角矩陣R的乘積形式，即Y＝QR，代入最小二乘解的一般表達(dá)式中得

式（13）又可轉(zhuǎn)化為非常容易求解的上三角方程組

綜上所述，本文所提出的求解式（7）中ai的方法，在算法復(fù)雜度和穩(wěn)定性方面明顯優(yōu)于傳統(tǒng)方法。值得提出的是，對(duì)于方程組（8）中bi的求解，也是一個(gè)超定方程組，其原始解法是利用式（4）遞推出x∧（n），再用最小二乘法求解。應(yīng)用本文所提出的QR分解算法，先令x∧（n）＝x（n）再去求解，可以進(jìn)一步減小算法的計(jì)算量。

另外，求解ai的方程組（12）和bi的方程組（8）不是一般類型的方程組。如方程組（12）中的Y是各條主對(duì)角線上元素皆相同的“Toeplitz矩陣”；方程組（8）中矩陣的每一行依次是極點(diǎn)zi的0，1，2，…，N－1次冪，為“Vandermonde矩陣”。用QR分解這兩類矩陣時(shí)，都有相應(yīng)的快速算法［8，9］，這為進(jìn)一步提高Prony算法的計(jì)算效率提供了有力的數(shù)學(xué)支持。

3 Prony算法的參數(shù)選取

3.1 模型階數(shù)的選取

一般情況下，信號(hào)中頻率成分的個(gè)數(shù)是未知的且常常含有噪聲；選取較大的階數(shù)進(jìn)行Prony分析時(shí)，辨識(shí)結(jié)果中所增加的余弦分量可被用來擬合信號(hào)中的噪聲成分，有利于提高算法的性能。但表征噪聲的頻率分量會(huì)干擾真實(shí)成分的辨別與提取，過大的模型階數(shù)反而會(huì)降低算法的性能。大量實(shí)驗(yàn)結(jié)果表明，最佳階數(shù)范圍和總的采樣點(diǎn)個(gè)數(shù)密切相關(guān)，一般為（0.35～0.45）倍的采樣點(diǎn)個(gè)數(shù)N。

模型階數(shù)與采樣點(diǎn)數(shù)的關(guān)系對(duì)算法性能的影響可用“方程的時(shí)間跨度”的概念來解釋（如圖1）。觀察方程組（12），樣本矩陣Y的每一行都是由階數(shù)p個(gè)采樣點(diǎn)構(gòu)成，各行所使用的數(shù)據(jù)在序列上只相差1個(gè)數(shù)據(jù)點(diǎn)，相當(dāng)于一個(gè)子窗口在移動(dòng)。p越大，每一行所使用的采樣點(diǎn)數(shù)越多，由這p個(gè)采樣點(diǎn)覆蓋信號(hào)的時(shí)間越長(zhǎng)，越能表達(dá)信號(hào)的信息；若p過大，由于總的采樣點(diǎn)數(shù)一定，樣本矩陣Y的行數(shù)（即等式個(gè)數(shù)）將很少，超定方程組將趨近于適定方程組甚至欠定方程組，以致降低算法的性能。

圖1 方程的時(shí)間跨度Fig.1 Time span of the equation

3.2 采樣頻率的選取

在滿足采樣定理的前提下，適當(dāng)增大采樣頻率可以提高算法的精度。但在時(shí)間窗確定的情況下，增大采樣頻率會(huì)導(dǎo)致運(yùn)算量大幅度上升。一些文獻(xiàn)中提出采樣頻率應(yīng)取4～10倍信號(hào)中的最高頻率較為合適［5］。這種方法的可靠性較低，事實(shí)上實(shí)際需要分析的信號(hào)的頻率范圍是變化的，即使是某一固定的信號(hào)，在不同的噪聲環(huán)境中用相同的采樣頻率辨識(shí)，算法的性能也往往不同。

本文提出一種自適應(yīng)采樣頻率的方法：對(duì)于方程組（12），在矩陣?yán)碚撝衅渥钚《私鉃閅TY的最小特征值對(duì)應(yīng)的特征向量。設(shè)YTY滿足

式中：λi為特征值；xi和wi是λi的右特征向量和左特征向量，它們的二范數(shù)均為1。假如YTY受到了一個(gè)擾動(dòng)變?yōu)閅TY＋εF，其中ε為二范數(shù)為1的矩陣F的系數(shù)，則式（15）和（16）將變?yōu)?/p>

此時(shí)的最小特征值對(duì)應(yīng)的特征向量設(shè)為xmin.index（ε），它和未經(jīng)擾動(dòng)時(shí)的最小特征值對(duì)應(yīng)的特征向量xmin.index之間的關(guān)系［10］為

即YTY受到了一個(gè)擾動(dòng)后方程組的解和未經(jīng)擾動(dòng)時(shí)的解的差之間存在著一個(gè)最大值，在不同的采樣頻率時(shí)該最大值是不同的。那么找到了此最大值最小時(shí)對(duì)應(yīng)的采樣頻率，也就找到了受擾動(dòng)影響最小的采樣頻率，定義靈敏度公式為

式（20）可以作為各個(gè)采樣頻率的評(píng)價(jià)指標(biāo)。在進(jìn)行Prony分析之前，先以較高的采樣頻率采集數(shù)據(jù)，再進(jìn)行各倍率的抽取以得出各采樣頻率時(shí)的數(shù)據(jù)，分別用這些數(shù)據(jù)形成方程組（12），再利用式（20）進(jìn)行靈敏度計(jì)算，最后用靈敏度最小時(shí)的方程組求取Prony算法中的ai。

3.3 時(shí)間窗的選取

在噪聲水平和采樣頻率確定的情況下，時(shí)間窗越長(zhǎng)估計(jì)的精度越高，但過長(zhǎng)則可能無法辨識(shí)出衰減快的分量同時(shí)會(huì)增加計(jì)算量。大量的實(shí)驗(yàn)證明，時(shí)間窗長(zhǎng)度應(yīng)為信號(hào)中最低頻率分量的1～2個(gè)周期。考慮到諧波中存在頻率低于50 Hz的分量的情況，時(shí)間窗取為80 ms能夠滿足絕大多數(shù)諧波檢測(cè)的需要。

4 真實(shí)頻率成分的提取方法

根據(jù)上文所述，為了提高算法的精度應(yīng)選取較大的模型階數(shù)，但分析結(jié)果中會(huì)出現(xiàn)一些信號(hào)中原本不存在的、表征噪聲的頻率分量，它們干擾了真實(shí)成分的辨別與提取。本文提出一種基于能量的提取方法，這種方法可靠且計(jì)算量十分小。

由于虛假成分是被用來擬合信號(hào)中的噪聲成分，所以其幅值都不會(huì)太大，或者呈現(xiàn)快速衰減的規(guī)律。所定義的能量要同時(shí)反映這兩個(gè)參數(shù)，既要設(shè)法突出幅值大的或者（和）衰減慢的頻率分量，

式中：Ei為各頻率分量的能量；Ai為各頻率分量的幅值；為代表各頻率分量的極點(diǎn)的n次冪；N為采樣點(diǎn)數(shù)；p為模型階數(shù)。

本文定義該能量表達(dá)式的理由如下：

①反映了各頻率分量幅值和衰減的大小，因?yàn)橛蓏的定義式可知，z的模值越大其衰減因子越小對(duì)應(yīng)衰減越慢；

②所需的計(jì)算量十分小，因?yàn)閦ni的計(jì)算已在Prony算法的第二個(gè)線性方程組（8）中計(jì)算完畢，整個(gè)能量計(jì)算只需要少量的乘法和加法計(jì)算。

真實(shí)頻率成分和噪聲成分在本文所定義的能量上具有明顯的差異：真實(shí)頻率成分的能量要遠(yuǎn)大于噪聲成分的能量。將各成分按能量從大到小排序，當(dāng)能量迅速地減小時(shí)可認(rèn)為后面的是虛假頻率分量。又要考慮計(jì)算量的大小以適合在線計(jì)算快速性需要。本文定義各頻率分量能量的表達(dá)式為

5 算例分析

為了闡述本文所提出的Prony算法完整的求解過程，說明其在電力系統(tǒng)諧波參數(shù)辨識(shí)中應(yīng)用的可行性，本文采用文獻(xiàn)［3］中沒有安裝補(bǔ)償裝置的電弧爐電流波形來進(jìn)行分析。信號(hào)是由基波（50 Hz，100 A，相位30°）、高次諧波（125 Hz，75 A，相位－60°）和間諧波（25 Hz，65 A，相位90°）組成，另外還加有信噪比為30 dB的高斯白噪聲，如圖2所示。

圖2 原始波形Fig.2 Original waveform

1）最優(yōu)采樣頻率的選取

根據(jù)本文3.2節(jié)所提出的方法，先以較高的采樣頻率10000 Hz采集數(shù)據(jù)80 ms得到800個(gè)數(shù)據(jù)點(diǎn)，再分別以每10、12、17、20、25個(gè)點(diǎn)進(jìn)行各倍率的抽取得到采樣頻率分別為1000 Hz、833 Hz、588 Hz、500 Hz、400 Hz的數(shù)據(jù)，分別用這些數(shù)據(jù)形成方程組（12）（模型階數(shù)均取為0.4倍的數(shù)據(jù)點(diǎn)數(shù)），再利用式（20）進(jìn)行靈敏度計(jì)算，得出結(jié)果見表1。

表1 各采樣頻率所對(duì)應(yīng)的靈敏度Tab.1 Sensitivities corresponding to each sampling frequency

由表1可知，采樣頻率在588 Hz時(shí)所對(duì)應(yīng)的靈敏度最小，所以取此時(shí)的數(shù)據(jù)進(jìn)行Prony運(yùn)算。

2）Prony辨識(shí)結(jié)果

應(yīng)用本文所提算法得到的Prony辨識(shí)結(jié)果及誤差分析如表2所示（已按能量從大到小排序，僅列出11項(xiàng)中的前8項(xiàng)），原信號(hào)波形與辨識(shí)結(jié)果擬合出的信號(hào)波形的誤差曲線見圖3。

由此可知，本算法在信噪比為30 d B強(qiáng)噪聲環(huán)境中各參數(shù)的辨識(shí)誤差均很小，辨識(shí)結(jié)果與原信號(hào)十分接近。

表2 Prony算法辨識(shí)結(jié)果及誤差分析Tab.2 Identification results of Prony algorithm and error analysis

3）真實(shí)頻率成分的提取

根據(jù)本文第4節(jié)所提出的方法，當(dāng)把各結(jié)果分量的能量從大到小排序，能量迅速減小時(shí)，可以認(rèn)為后面的分量是虛假分量。

由圖4可知前3個(gè)分量是真實(shí)的頻率成分，后面的分量是由噪聲引起的，與原信號(hào)的情況相符。

圖3 原始波形與Prony算法擬合波形的誤差曲線Fig.3 Error curve of the original waveform and the fitting wave of the Prony algorithm

圖4 Prony辨識(shí)結(jié)果各分量的能量大小變化Fig.4 Energy variation of the Prony identification result components

（4）與原始Prony算法的比較

用本文所提出的方法和原始方法分別對(duì)原信號(hào)（仍為30 dB的信噪比但噪聲不同）獨(dú)立運(yùn)算100次，計(jì)算耗時(shí)（PC機(jī)主頻2.66GHz，480M內(nèi)存）比較結(jié)果如表3所示，辨識(shí)精度比較結(jié)果如表4所示。在原始方法的計(jì)算過程中，還出現(xiàn)了幾次“矩陣奇異”的警告。由此可知，本文所提方法的計(jì)算耗時(shí)明顯優(yōu)于原始方法，在辨識(shí)精度和穩(wěn)定性方面也比原始方法好。

表3 本文方法和原始方法的耗時(shí)比較結(jié)果Tab.3 Time－consuming comparison of the proposed and original method s

表4 本文方法和原始方法的辨識(shí)精度比較結(jié)果Tab.4 Identification accuracy comparison of the proposed and original method

6 結(jié)論

（1）本文提出的Prony算法改進(jìn)方案應(yīng)用于電力系統(tǒng)諧波的檢測(cè)，無論是對(duì)工頻整數(shù)次還是非整數(shù)次諧波，都能一次性地辨識(shí)出信號(hào)中各成分的幅值、頻率、初相等參數(shù)，且具有較高的精度。

（2）本文提出了減少算法復(fù)雜度的方法，以及在線自適應(yīng)采樣頻率選取方案、模型階數(shù)和時(shí)間窗的選取原則，縮短了算法的計(jì)算耗時(shí)，提高了辨識(shí)精度。

（3）本文提出了真實(shí)頻率成分的提取方法，解決了應(yīng)用Prony算法進(jìn)行諧波檢測(cè)時(shí)較高模型階數(shù)帶來的虛假成分問題，進(jìn)一步提高了Prony算法的可靠性。

［1］ 楊洪耕，惠錦，侯鵬（Yang Honggeng，Hui Jin，Hou Peng）.電力系統(tǒng)諧波和間諧波檢測(cè)方法綜述（Detection methods of harmonics and inter－h(huán)armonics in power system）［J］.電力系統(tǒng)及其自動(dòng)化學(xué)報(bào)（Proceedings of the CSU－EPSA），2010，22（2）：65－69.

［2］ 丁屹峰，程浩忠，呂干云，等（Ding Yifeng，Cheng Haozhong，LüGanyun，et al）.基于Prony算法的諧波和間諧波頻譜估計(jì)（Spectrum estimation of harmonics and interharmonics based on Prony algorithm）［J］.電工技術(shù)學(xué)報(bào)（Transactions of China Electrotechnical Society），2005，20（10）：94－97.

［3］ 黃云江，謝維波（Huang Yunjiang，Xie Weibo）.基于迭代Prony算法的高精度諧波檢測(cè)（Precise harmonic analysis based iterative Prony method）［J］.電氣應(yīng)用（Electrotechnical Application），2007，26（4）：97－100.

［4］ Costa F F，Cardoso A J M.Harmonic and interharmonic identification based on improved Prony＇s meth－od［C］∥32nd Annual Conference on IEEE Industrial Electronics，Paris，F(xiàn)rance：2006.

［5］ 束洪春.電力工程信號(hào)處理應(yīng)用［M］.北京：科學(xué)出版社，2009.

［6］ Kamel N，Sankaran P，Venkatesh B.Fault－current predetermination using time－limited CT secondary side measurements by the covariance－Prony method［J］.IEE Proceedings：Generation，Transmission and Distribution，2004，151（6）：735－739.

［7］ 魏木生.廣義最小二乘問題的理論和計(jì)算［M］.北京：科學(xué)出版社，2006.

［8］ 王治平，鄭慧嬈，張莉（Wang Zhiping，Zheng Huirao，Zhang Li）.Hankel矩陣類的一種快速分解算法（A fast decomposition algorithm of Hankel matrix class）［J］.數(shù)學(xué)雜志（Journal of Mathematics），1998，18（S1）：145－147.

［9］ Demeure Cedric J，Scharf Louis L.Fast least squares

solution of Vandermonde systems of equations［C］∥IEEE International Conference on Acoustics，Speech

and Signal Processing，Glasgow，Scotland：1989.［10］Lee Joon－Ho，Kim Hyo－Tae.Selecting sampling interval of transient response for the improved Prony method［J］.IEEE Trans on Antennas and Propagation，2003，51（1）：74－77.