Prony算法在諧波、間諧波參數辨識中的應用

楊玉坤，楊明玉

（1.吉林省電力勘測設計院，長春 130022；2.華北電力大學電力系統保護與動態安全監控教育部重點實驗室，保定 071003）

隨著電力電子裝置在電力系統中的廣泛應用以及非線性負荷的日益增多，電力系統中的諧波污染越來越嚴重和復雜。諧波污染除了與基波成整數倍的諧波外，還存在許多非整數倍的間諧波。間諧波的存在不僅會對電力系統造成更大的危害，更增加了諧波檢測的難度。

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

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

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

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

定義

可由歐拉公式得

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

的齊次解其特征多項式為

并代入式（4）可得

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

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

2 Prony算法求解方法的改進

對于式（7）中ai的求解，通常的方法是構造誤差ε（n）的平方和J（a）為代價函數，即

令

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

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

可簡寫為Ya＝y

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

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

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

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

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

3 Prony算法的參數選取

3.1 模型階數的選取

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

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

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

3.2 采樣頻率的選取

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

本文提出一種自適應采樣頻率的方法：對于方程組（12），在矩陣理論中其最小二乘解為YTY的最小特征值對應的特征向量。設YTY滿足

式中：λi為特征值；xi和wi是λi的右特征向量和左特征向量，它們的二范數均為1。假如YTY受到了一個擾動變為YTY＋εF，其中ε為二范數為1的矩陣F的系數，則式（15）和（16）將變為

此時的最小特征值對應的特征向量設為xmin.index（ε），它和未經擾動時的最小特征值對應的特征向量xmin.index之間的關系［10］為

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

式（20）可以作為各個采樣頻率的評價指標。在進行Prony分析之前，先以較高的采樣頻率采集數據，再進行各倍率的抽取以得出各采樣頻率時的數據，分別用這些數據形成方程組（12），再利用式（20）進行靈敏度計算，最后用靈敏度最小時的方程組求取Prony算法中的ai。

3.3 時間窗的選取

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

4 真實頻率成分的提取方法

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

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

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

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

①反映了各頻率分量幅值和衰減的大小，因為由z的定義式可知，z的模值越大其衰減因子越小對應衰減越慢；

②所需的計算量十分小，因為zni的計算已在Prony算法的第二個線性方程組（8）中計算完畢，整個能量計算只需要少量的乘法和加法計算。

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

5 算例分析

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

圖2 原始波形Fig.2 Original waveform

1）最優采樣頻率的選取

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

表1 各采樣頻率所對應的靈敏度Tab.1 Sensitivities corresponding to each sampling frequency

由表1可知，采樣頻率在588 Hz時所對應的靈敏度最小，所以取此時的數據進行Prony運算。

2）Prony辨識結果

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

由此可知，本算法在信噪比為30 d B強噪聲環境中各參數的辨識誤差均很小，辨識結果與原信號十分接近。

表2 Prony算法辨識結果及誤差分析Tab.2 Identification results of Prony algorithm and error analysis

3）真實頻率成分的提取

根據本文第4節所提出的方法，當把各結果分量的能量從大到小排序，能量迅速減小時，可以認為后面的分量是虛假分量。

由圖4可知前3個分量是真實的頻率成分，后面的分量是由噪聲引起的，與原信號的情況相符。

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

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

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

用本文所提出的方法和原始方法分別對原信號（仍為30 dB的信噪比但噪聲不同）獨立運算100次，計算耗時（PC機主頻2.66GHz，480M內存）比較結果如表3所示，辨識精度比較結果如表4所示。在原始方法的計算過程中，還出現了幾次“矩陣奇異”的警告。由此可知，本文所提方法的計算耗時明顯優于原始方法，在辨識精度和穩定性方面也比原始方法好。

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

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

6 結論

（1）本文提出的Prony算法改進方案應用于電力系統諧波的檢測，無論是對工頻整數次還是非整數次諧波，都能一次性地辨識出信號中各成分的幅值、頻率、初相等參數，且具有較高的精度。

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

（3）本文提出了真實頻率成分的提取方法，解決了應用Prony算法進行諧波檢測時較高模型階數帶來的虛假成分問題，進一步提高了Prony算法的可靠性。

［1］ 楊洪耕，惠錦，侯鵬（Yang Honggeng，Hui Jin，Hou Peng）.電力系統諧波和間諧波檢測方法綜述（Detection methods of harmonics and inter－harmonics in power system）［J］.電力系統及其自動化學報（Proceedings of the CSU－EPSA），2010，22（2）：65－69.

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

［3］ 黃云江，謝維波（Huang Yunjiang，Xie Weibo）.基于迭代Prony算法的高精度諧波檢測（Precise harmonic analysis based iterative Prony method）［J］.電氣應用（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，France：2006.

［5］ 束洪春.電力工程信號處理應用［M］.北京：科學出版社，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］ 魏木生.廣義最小二乘問題的理論和計算［M］.北京：科學出版社，2006.

［8］ 王治平，鄭慧嬈，張莉（Wang Zhiping，Zheng Huirao，Zhang Li）.Hankel矩陣類的一種快速分解算法（A fast decomposition algorithm of Hankel matrix class）［J］.數學雜志（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.