基于譜序列變換的高精度諧波參數估計算法

楊喜,2，汪旭明3，陳炳權，張仁民，王向明4，雷可君

(1.吉首大學信息科學與工程學院，湖南吉首，416000；2.東南大學信息科學與工程學院，江蘇南京，210096；3.吉首大學物理與機電工程學院，湖南吉首，416000；4.上交所技術有限責任公司，上海，200120)

隨著現代電力電子技術的快速發展，越來越多的非線性元器件被應用于電力系統[1]，在電網中產生大量諧波污染。諧波污染使得電能質量下降，從而嚴重影響電力系統的安全和經濟運行。為了更好地進行諧波污染治理，必須對諧波參數進行高精度估計[2]。快速傅里葉變換(fast Fourier transform, FFT)是一種經典的諧波參數估計方法。然而，由于在實際運行過程中電網基波處于波動狀態，因此，實際的采樣過程通常處于非同步狀態。直接采用FFT算法對諧波信號進行參數估計，其過程相當于對原始信號序列施加矩形窗。矩形窗實現方式雖然簡單，但這種窗所引起的頻譜泄露非常嚴重，其旁瓣峰值衰減僅為13 dB，由此導致FFT 算法對應的估計結果精度較低。為了減小頻譜泄漏在諧波參數估計過程中帶來的負面影響，采用在時域加窗處理和頻域譜線插值處理相結合的方法來提高估計精度[3-13]。在前期的時域加窗處理過程中，常用的窗函數包括Hanning 窗、Blackman 窗、Nuttall 窗、Rife-Vincent窗和MSD窗。這些窗函數相對簡單，但其旁瓣衰減速度有限而不利于抑制頻譜泄露。為進一步加速窗函數旁瓣衰減，近年來，一些更加復雜的優化窗相繼被提出，如乘法窗[11]、卷積窗[12-13]等。然而，當采用一些復雜的窗函數時，頻譜泄露的抑制效果盡管有所提高，但另一方面也將不可避免地導致更高的計算復雜度而不利于信號參數的實時估計。通過加窗處理可以在一定程度上抑制長距離頻譜泄漏的影響，但其對短距離范圍內頻譜泄露的抑制效果一般，因此，在隨后的頻域譜線插值處理過程中通過選擇峰值譜線附近不同數量的譜線進行插值，以達到減小短距離頻譜泄漏的影響，從而進一步提高參數估計精度。例如，文獻[3，4，8]提出了加窗雙譜線插值算法，文獻[5，9]提出了加窗三譜線插值算法，文獻[6，10]提出了加窗四譜線插值算法，文獻[7]則提出了一種加窗六譜線插值算法。這些插值方法在一定程度上都可以提高諧波參數估計的精度。值得注意的是，上述基于不同數量譜線的插值方法有著各自的特點：對于利用較少譜線進行插值的方法，其插值公式推導和相關參數數值計算的復雜度較低，但其諧波參數估計的精度也相應較低；而對于利用較多譜線進行插值的方法，雖然參數估計的精度有所提高，但其計算復雜度也相應提高。相比較而言，文獻[9]提出的基于Nuttall 窗的三峰插值諧波算法較好地兼顧了參數估計精度和計算復雜度。與此同時，考慮到FFT 方法在工程上簡單易行，一些研究者提出了基于FFT的改進諧波參數分析方法，如文獻[14-15]提出了基于FFT 序列進行變換處理的諧波參數估計算法，在此基礎上，文獻[16]通過在變換處理之前引入時域加窗處理環節以進一步減小頻譜泄露，從而提高后繼處理過程中諧波參數估計精度。整體而言，基于FFT 的改進算法雖然比前述基于加窗與插值處理相結合的算法其實現復雜度有所降低，但在實際應用過程中，諧波參數估計的精度仍有待進一步提高。本文在分析信號譜序列衰減特征的基礎上，提出一種新的基于譜序列變換的高精度諧波參數估計算法。該算法通過對譜序列進行特定的加權變換處理，減少了頻譜泄漏造成的諧波間干擾，極大地提高了諧波參數估計的精度。與此同時，由于該方法只需對信號FFT 變換后的譜序列直接進行簡單的變換處理，因此，具有計算復雜度低的特點。

1 信號頻譜分析

設單頻率無限長余弦信號為

其中：A0，f0和φ0分別為信號的幅值、頻率、相位。信號對應的抽樣序列xd(n)為

其中：ω0= 2πf0/fs；fs為抽樣頻率。將xd(n)截斷成1個長為N點的序列x(n)，這相當于將xd(n)乘上1個長為N點的窗函數wN(n)：

注意到對x(n)進行FFT 變換時相當于預先對x(n)加了1個矩形窗，此時，wN(n)可以表示為

原信號的復指數序列形式可以表示為

利用頻域卷積定理，x(n)的離散時間傅里葉變換可表示為

式中：Xd(ejω)表示原序列xd(n)的離散時間傅里葉變換，

其中：δ(ω)為沖激函數。WN(ejω)表示wN(n)的離散時間傅里葉變換，

因此，原序列的頻譜為

利用沖激函數的性質可以得到

式(11)為連續復譜函數，由于該復譜函數的正、負頻率點具有對稱性，且負頻點對諧波參數估計的影響較小，所以，在分析過程中一般忽略負頻率成分的影響，而只對信號的正頻率部分加以分析和利用。在頻域內取N點，以采樣間隔Δf=fs/N對連續譜進行離散化。相應地，在式(11)中令ω= 2πn/N(n= 0,1,…,N- 1)，ω0= 2πk0/N即得到離散后的譜序列為

又由于N?1，故譜序列可進一步近似為

2 基于譜序列變換的諧波分析

2.1 譜序列變換原理

由式(13)可知，當處于同步采樣時僅當n=k0時，有X(n)=A0Nejφ0/2。此時，譜線的峰值恰好等于準確值，而在其他譜線上有X(n)= 0。電網中基波頻率處于波動狀態使得實際采樣總是處于非同步狀態，即k0為一個非整數，故設k0=k1+Δk，其中k1為k0的整數部分，Δk為其小數部分，故有

利用三角函數和差化積公式，式(14)可以表示為

由式(15)可知，當信號處于非同步采樣時，頻譜的峰值不等于準確值，而是在其他譜線上產生了分量，即產生了頻譜泄漏，并且各次諧波之間將產生相互干擾，從而使得對諧波參數的準確估計十分困難，因此，如何降低頻譜泄漏成為高精度諧波參數估計的關鍵。進一步令

則式(15)可以簡潔地表示為

式(17)表明信號的幅度譜與||成反比。這意味著隨著||增大，遠離真實頻率位置的譜線幅值只是以O()的速度衰減。由于這些譜線的幅值衰減速度不夠快，頻譜泄露將對諧波參數的估計造成很大干擾，極大地影響諧波參數估計的精度。為了減小頻譜泄漏對諧波參數估計精度帶來的負面影響，需要盡量加快偏離真實頻率位置譜線幅值的衰減速度，為此，本文新設計出一個衰減速度為O(-11)的新序列：

需指出的是，盡管可以選擇更多的譜線來構造衰減速度更快的新序列，但考慮到算法復雜度和估計精度的平衡，在式(18)所構造的新序列中只是選擇了真實頻率位置左右各5根譜線參與諧波參數的估計。注意到Y(n)可以看作1個有理多項式，因此，根據代數理論，新序列Y(n)可以進一步展開為如下表達式：

聯立式(18)和式(19)可求得

由式(15)～(17)可得到

式中：i= 1,2,…,5。所以，Y(n)可以等價地表示為

可見新的譜序列Y(n)實際上是將真實頻率附近的11 個點進行加權處理，由此頻譜序列的衰減速度由原來的1/||變為衰減速度得到極大提高，使得離真實譜線較遠處頻點的幅值急劇減小，從而有效地抑制頻譜泄漏。為驗證本文算法對頻譜泄漏的抑制效果，現對1個由基波和三次諧波組成的復合信號進行譜分析，信號為

其中：f0= 50 Hz；fs= 520.5Hz。對信號直接進行FFT處理和基于譜序列變換后得到的頻譜特征分別如圖1和圖2所示(其中，頻點為頻域中離散譜線序號)。對比圖1和圖2可見：利用FFT算法對信號直接處理得到的譜序列頻譜泄漏嚴重，而采用本文所提方法對譜序列進行變換后，新的譜序列中具有較大幅值的譜線主要集中在峰值譜線附近，減少了諧波之間的干擾，從而達到了進一步抑制頻譜泄漏的目的，進而提高了諧波參數估計的精度。

圖1 直接FFT頻譜圖Fig.1 Spectrum graph of direct FFT

2.2 諧波參數估計方法

非同步采樣使得諧波真實頻率的位置與峰值譜線的位置并不對應，而是位于相鄰的2根譜線之間。設真實頻率偏離峰值譜線的頻率為Δk，若能求得Δk，則可準確計算真實譜線的位置，進而估計相關參數。為此，先在變換后的序列中搜尋幅值最大的2 根譜線Y(k1)和Y(k1+ 1)。由式(16)和式(18)得：

圖2 基于譜序列變換后的頻譜圖Fig.2 Spectrum graph of spectral sequence transformation

設|Y(k1)|和|Y(k1+ 1)|的比值為ε0，則

可以求出偏移量為

結合式(16)和式(18)不難得到諧波幅值和相位的計算表達式：

其中：arg( )Y(k1)為Y(k1)的相位。經歸納，諧波參數估計的具體步驟如下。

1)對序列x(n)進行FFT 運算，得到譜序列X(n)；

2)由原譜序列X(n)根據式(22)計算得到新的譜序列Y(n)；

3)在變換后的序列Y(n)中搜索幅值最大的2根譜線Y(k1)和Y(k1+ 1)；

4)求出|Y(k1)|和|Y(k1+ 1)|的比值ε0；

5)通過ε0與Δk之間的關系式(27)計算偏移量Δk；

6)根據式(28)和(29)計算出幅值A0和相位φ0。

基于譜序列變換的諧波參數估計流程圖如圖3所示。

圖3 諧波參數估計的流程圖Fig.3 Flow chart of harmonic parameter estimation

上面以單頻信號為例介紹了諧波參數估計的方法，該方法可以方便地推廣到用于由基波和諧波組成的復合信號的參數估計。具體而言，在復合信號參數估計中，先計算出基波頻率f1=(k1+Δk)·fs/N，再將第i次諧波的頻域搜索范圍設置為[i·f1- 5,i·f1+ 5]，并在該范圍內搜索出幅值最大的2根譜線，根據式(24)～(27)即可依次計算出各次諧波的幅值和相位。

3 仿真與分析

為了驗證本文所提算法的有效性，利用本文算法、文獻[9]中算法和文獻[16]中算法分別對單頻信號和復合信號及相關的重要場景進行仿真實驗對比。

3.1 單頻信號場景

對單頻信號作參數估計：

其中：A= 25V；f0= 35.5Hz；fs= 100 Hz；φ=1.5rad；采樣點數N為64。

利用文獻[9]中算法計算的幅值A和相位φ為：

利用文獻[16]中算法計算的幅值A和相位φ為：

利用本文算法計算的幅值A和相位φ為：

由計算結果可知，利用文獻[9]中算法估計幅值和相位的相對誤差分別為4.80×10-9和9.34×10-8，利用文獻[16]中算法估計幅值和相位的相對誤差分別為8.43×10-9和4.33×10-8，而利用本文算法估計幅值和相位的相對誤差分別為3.55×10-11和3.24×10-10。顯然，對于單頻信號而言，本文算法比文獻[9]和文獻[16]中所提算法具有更高的參數估計精度。

3.2 多頻信號場景

設1個由基波和2～21次諧波組成的復合信號，其基波頻率為f1= 50.2 Hz，表達式為

式中：A1為基波的幅值；Ai為各次諧波的幅值；φ1為基波的相位；φi為各次諧波的相位。在仿真過程中，具體參數設置如表1所示。

設置采樣頻率fs= 2 300 Hz，采樣點數N=1 024，對信號分別采用本文算法、文獻[9]中算法和文獻[16]中算法進行參數估計，這3 種算法的幅值和相位估計相對誤差對比結果如圖4所示。從圖4(a)可見：對于基波的幅值估計相對誤差，本文算法所得結果比文獻[9]中算法所得結果高近7個數量級，比文獻[16]中算法結果高3 個以上數量級；對于基波以外的諧波，本文算法得到的估計相對誤差范圍為7.39×10-11～3.03×10-14，文獻[9]中算法結果相對誤差范圍為8.74×10-8～4.82×10-10，文獻[16]中算法結果相對誤差范圍為7.18×10-9～6.10×10-11。從圖4(b)可見：對于基波的相位估計相對誤差，本文算法所得結果比文獻[9]中算法結果高將近4個數量級，比文獻[16]中算法結果高近3 個數量級；進一步，對于基波以外的各次諧波，本文所提算法的相對誤差范圍為1.28×10-9～8.14×10-13，而文獻[9]中算法結果的相對誤差范圍為2.29×10-6～2.08×10-9，文獻[16]中算法結果的相對誤差范圍則為6.21×10-7～6.34×10-10。總體來看，對于基波和諧波組成的復合信號，本文算法在幅值和相位參數的估計結果精度與文獻[9]和文獻[16]中算法結果精度相比顯著提高。

表1 基波和各次諧波的幅值和相位Table1 Amplitude and phase of the fundamental and harmonics

圖4 不同算法下諧波幅值和相位估計相對誤差對比Fig.4 Comparison of harmonic amplitude and phase estimation relative error under different algorithms

在采樣點數N為512，1 024 和2 048 時，幅值、相位估計相對誤差仿真結果見圖5。從圖5可以看出：當采樣點數由512 增加到1 024 時，幅值估計相對誤差至少提高了2個數量級，對4次和11次諧波則提高了近5個數量級；相位估計相對誤差至少提高了3個數量級，對7次和14次諧波則提高了近4 個數量級；當采樣點數由512 增加到2 048時，幅值估計相對誤差至少提高了3個數量級，對2次和15次諧波則提高了近6個數量級；同時，相位估計相對誤差至少提高了3個數量級，對5次和9 次諧波則提高了近7 個數量級。顯然，采樣點數的增加可以有效地提高算法的估計精度。

圖5 采樣點數N不同時諧波幅值和相位估計相對誤差對比Fig.5 Comparison of harmonic amplitude and phase estimation relative error when sampling points are different

表2 不同基波頻率下幅值估計相對誤差Table2 Amplitude estimation relative error at different fundamental frequencies

表3 不同基波頻率相位估計相對誤差Table3 Phase estimation relative error at different fundamental frequencies

3.3 基波頻率變化場景

電網實際運行過程中的基波頻率并不是固定不變的，通常在50 Hz左右波動，因此，考察算法在基波頻率變化情況下參數估計性能的穩定性非常必要。GB/T 15945—2008 規定電力系統基波頻率偏差最大范圍為±0.5 Hz，故在仿真過程中基波頻率設定在49.5～50.5 Hz[17-18]。信號時域表達式如式(31)所示，仿真過程中設置采樣頻率fs=2 300 Hz，采樣點數N=1 024。表2和表3所示仿真結果表明：在基波頻率變化時，本文所提算法仍然表現出優良的估計性能和良好的穩定性。

3.4 間諧波場景分析

電網在實際運行過程中可能會產生間諧波，其存在會造成電壓的閃變并對電網產生危害，因此，在處理諧波過程中，對間諧波參數的高精度估計非常必要[19-20]。為了驗證本文算法對間諧波參數估計的有效性，在仿真過程中采用1個由7次間諧波組成的復合信號(其基波頻率f1= 50.1Hz)：

具體參數設置如表4所示。另外，在仿真過程中設置采樣頻率fs=2 300 Hz，采樣點數N=1 024。采用文獻[9]中算法、文獻[16]中算法和本文算法分別對上述間諧波信號的幅值和相位參數進行估計，得到各次間諧波的幅值、相位估計的相對誤差如表5和表6所示。從仿真結果可以看出：文獻[9]和文獻[16]中算法在一定程度上提高了間諧波參數估計的精度，但本文算法與文獻[9]和[16]中算法相比，估計精度顯著提高。

表4 各次間諧波的幅值和相位Table4 Amplitude and phase of each inter-harmonic

表5 采用不同算法的時間諧波幅值估計相對誤差Table5 Amplitude estimation relative error of inter-harmonics under different algorithms

表6 采用不同算法的時間諧波相位估計相對誤差Table6 Phase estimation relative error of inter-harmonics under different algorithms

4 結論

1）在對諧波信號FFT 譜序列的衰減特征進行分析的基礎上，通過對原譜序列實施特定加權變換，加快非真實頻率幅值的衰減速度，從而達到有效抑制頻譜泄漏的目的，并在此基礎上給出了一種諧波幅值和相位的高精度估計方法。

2）該算法能有效抑制頻譜泄漏。與經典的加窗插值算法和FFT 改進算法所得結果相比，該算法估計精度顯著提高，且該算法在基波頻率變化及間諧波場景下均表現出優良的估計性能。

3）與經典的加窗插值算法相比，本文所提出的新算法無需復雜的加窗處理環節，也無需復雜的插值處理操作，實現過程更簡單，計算復雜度低。