改進加窗頻譜峰值擬合算法及諧波分析應用

溫 和 滕召勝 黎福海 王 永 胡曉光

（1.湖南大學電氣與信息工程學院 長沙 4100822.北京航天航空大學自動化工程學院 北京 100191）

1 引言

電力系統中，諧波的準確檢測能夠為諧波源定位、諧波治理、電能質量管理等提供科學依據，對電力系統的安全與穩定運行具有重要的意義[1-3]?？焖俑道锶~變換（FFT）因其易于實現，是目前最主要的電力諧波分析方法，但非同步采樣時，FFT固有的頻譜泄漏和柵欄效應影響了諧波參數（幅值、相位、頻率）估計的準確性[4-6]。

國內外學者研究并提出了基于Hanning窗[7-8]、Bl ackman-Harris窗[9]、 Rife-Vincent窗[10]、 Nutt all窗[11-12]、矩形卷積窗[13-14]、三角自卷積窗[15-16]（Triangular Self-Convolution Windows，TSCW）等的加權電力諧波分析方法，在一定程度上提高了諧波檢測的準確度[17-18]。但是采用上述窗函數進行諧波參數分析時，仍存在離散頻譜辨識困難、修正公式求解復雜、計算準確度低等問題[19-20]，特別是對弱信號提取或對復雜信號進行分析時仍存在較大誤差。

在對非同步采樣和非整周期截斷時的信號頻譜進行分析的基礎上，本文結合三角自卷積窗的頻譜特性，采用最小二乘法（Least Square Method，LSM）逼近諧波參數分析的多項式，建立簡便、易行且可根據計算精度要求進行調節的諧波幅值、頻率、初相角參數計算式，該方法具有計算量小、實時性高、易于嵌入式系統實現的優點，較好地解決了非同步采樣時離散頻譜校正中存在的計算準確度與實時性的矛盾。仿真結果表明：在非同步采樣和非整周期截斷的條件下，本文算法能有效消除各次諧波相互干擾、抑制白噪聲和基波頻率波動對諧波參數分析的影響，提高諧波參數分析準確度。

2 頻譜泄漏與三角自卷積窗函數

2.1 非同步采樣時的頻譜泄漏和柵欄效應

考慮一個周期為T、包含H次諧波的離散信號

式中，h為諧波次數；f0=1 T代表基波頻率；Ah、?h分別代表第h次諧波的幅值和初相角；為采樣頻率； n =0,1,…, N?1。

信號在觀測周期Tw轉換為 N點離散序列相當于被矩形窗wR(n)截斷，其離散傅里葉變換為

同步采樣時，觀測周期Tw為信號周期T的整數倍。而非同步采樣時，離散頻譜中信號各次諧波頻率與頻率分辨率之間不再為對應的整數倍關系，即出現所謂的柵欄效應。信號各次頻率分量不再是單一譜線，而是分布于整個頻率軸，即能量不再集中，并向臨近譜線泄漏，即所謂的頻譜泄漏現象。

非同步采樣情況下直接采用FFT進行諧波參數分析，即直接利用頻譜峰值計算諧波參數時，頻譜泄漏和柵欄效應引起的誤差較大。此外，非同步采樣時，若信號包含有多個頻率分量，則其各頻率分量會產生泄漏疊加干擾（矢量之和），影響諧波參數檢測的準確度。

2.2 三角自卷積窗與頻譜泄漏抑制

由式（2）可知，頻率 f0處的頻譜形狀等于信號截斷時所加窗函數的窗譜形狀。因此，通過選擇窗函數減少頻譜泄漏，降低諧波分量之間的相互干擾。三角自卷積窗wT(n) 具有優良旁瓣性能，可由若干個三角窗wg(m) 進行自卷積運算獲得[23]，即

式中，n∈[0,N?1]；m∈[0,M?1]；p為參與卷積的三角窗的個數，稱為窗的階數。為便于計算，卷積運算后可對序列進行補零，將p階三角自卷積窗的長度調整為N=pM。

根據卷積定理，函數在時域進行卷積等效于頻域相乘，因此p階三角自卷積窗的頻率響應為

圖1給出了三角窗（M=128）進行1～4階自卷積得到1～4階三角自卷積窗的幅頻響應曲線。圖1中標注了各階三角自卷積窗的旁瓣峰值電平和旁瓣衰減速率。

圖1 p階三角自卷積窗幅頻響應Fig.1 Frequency responses of the p-order TSCW

表1為4階和8階三角自卷積窗與常用窗函數在主瓣寬度、旁瓣峰值電平和旁瓣衰減速率等方面的比較。由圖1和表1可見，三角自卷積窗的旁瓣電平與衰減速率隨卷積階數的升高而加強；與現有經典窗函數相比，三角自卷積窗具有較好的旁瓣峰值電平和旁瓣衰減速率，能有效減少頻譜泄漏引入的誤差。

表1 窗函數性能比較Tab.1 Comparison of window characteristics

3 基于LSM的離散頻譜峰值擬合算法

3.1 基于TSCW的離散頻譜峰值校正原理

對式（1）所示信號 x(n) 加三角自卷積窗后，其加窗離散傅里葉變換為

式中，k =0,1,L,N?1；WT(*)是三角自卷積窗的離散頻譜函數。

不失一般性，設當前測量的為第h項諧波（h≤H），為簡單起見，忽略其余各次諧波對第h項諧波的泄漏影響，此時，式（6）變為

設第h項諧波的頻率hf0在離散頻譜中對應的位置為hk0處，即 hf0=hk0Δf。參見圖2，非同步采樣時，hk0為非整數，即第h項諧波頻率hf0與離散譜線不重合，且位于離散頻譜幅值最大譜線k1和次大譜線k2之間（k1≤hk0≤k2=k1+1）。

參見圖 2，令這兩條峰值譜線的幅值分別為y1和y2，即

式（10）即為第h項諧波離散頻譜的峰值譜線（第k1與k2根譜線）與真實頻率譜線（第hk0根）的相對位置參數λ的約束條件。因此，理論上，可通過對離散頻譜峰值進行校正，得到參數λ，從而實現諧波參數求解。

圖2 非同步采樣時離散頻譜的柵欄效應Fig.2 Picket fence effect of the spectral with non-synchronized sampling

3.2 基于LSM的離散頻譜峰值擬合算法

式（10）為復雜的有理式，若直接求其反函數，計算量較大。此外，式（10）建立的基礎是忽略諧波間的相互泄漏影響，僅考慮非同步采樣引起的頻譜泄漏，即所描述的頻譜峰值與相對位置參數λ之間的關系為一種近似關系。因此，本文考慮采用基于LSM的多項式擬合方法獲得求解λ的計算式。

根據式（10），取L組待擬合試驗數據(ηi, λi)（i=0,1,L,L?1），由于待擬合試驗數據本身存在偏差，因此不要求擬合多項式λ=S(η)經過所有已知待擬合數據點(ηi,λi)，而僅要求在給定點ηi上的誤差γi=S(ηi)?λi按誤差平方和最小，即

顯然，式（12）為 a0,a1,…,aK的多元函數，因此多項式擬合過程即為求式（12）極值的過程。由多元函數求極值的必要條件可得

式中，j =0,1,…,K?1。

由式（13）建立了關于a0,a1,L,aK的線性方程組，解出系數ak，從而可得基于 LSM的頻譜峰值擬合多項式 SK(η)。

3.3 精度可控的多項式擬合與諧波參數求解算法

由于 hk0=k1+λ，可知λ取值范圍為[0，1]，在該范圍內取一組λ 值，可由式（10）得到對應的一組η值，即建立待擬合試驗數據(ηi, λi)。根據式（13），可求解擬合多項式系數 a0,a1,…,aK。

工程實踐中，可以根據測量精度的要求確定基于LSM的擬合多項式的次數。為確保諧波參數計算實時性，擬合多項式一般不超過7次。根據式（13），得到的基于4階三角自卷積窗的頻譜峰值最小二乘擬合5次多項式為

式中，參數 λ的擬合誤差為 9.1587×10?5。

基于4階三角自卷積窗的頻譜峰值最小二乘擬合6次多項式為

式中，參數λ的擬合誤差為 1.7016×10?5。

基于4階三角自卷積窗的頻譜峰值最小二乘擬合7次多項式為

式中，參數λ的擬合誤差為3.3323×10?6。

由于第h項諧波對應的真實頻率譜線位置應為第h k0根，因此第h項諧波的頻率為

由式（10）可得，第h項諧波的幅值計算式為

同理，可得第h項諧波的初相角為

信號中基波和各次諧波分量參數的計算可按上述過程逐次計算。

4 仿真與分析

4.1 復雜信號的諧波分析

采用文獻[12]給出的含 2～21次諧波的信號進行分析

式中，基波頻率f1=50.5Hz；采樣頻率fs=2500Hz；基波和各次諧波的幅值Ah和初相位h?在表2中給出。仿真實驗結果由表3和表4給出，其中a E-b代表a×10?b。其中 Blackman、Blackman-Harris、Nuttall窗采用的數據長度為N=512；4階三角自卷積窗分別長度為N=512和N=1024（最小二乘多項式擬合次數均為7次）。

表2 復雜信號的基波及諧波參數Tab.2 Parameters of the simulated harmonic signal

表3 幅值與頻率相對百分比誤差Tab.3 Percentage relative errors in calculating amplitude and frequency（%）

表4 初相位測量相對百分比誤差Tab.4 Percentage relative errors in calculating phase（%）

表3中，Ef1表示基波頻率的測量值相對于真值的誤差。由表3和表4可見，采用4階三角自卷積窗進行加權，利用本文提出的多項式擬合算法進行離散頻譜校正后，N=512時頻率計算的相對百分比誤差為2.9×10?7%，N=1024時頻率計算的相對百分比誤差為-6.4×10?9%。當 N=1024時，幅值相對誤差EAh≤0.00094%，相位相對誤差Eφh＜0.012%，可實現復雜諧波信號參數的高準確度分析。與Blackman窗、Blackman-Harris窗和 Nuttall窗相比，基于 4階三角自卷積窗的最小二乘多項式擬合頻譜校正算法具有更高的諧波分析準確度，且計算簡單，易于實現。

4.2 白噪聲影響下的諧波分析

在白噪聲存在的情況下，對式（20）所示的信號添加白噪聲，信噪比（SNR）的范圍為[10dB，100dB]，變化步長為10dB。分別采用Blackman窗、Blackman-Harris窗、Nuttall窗和本文算法（4階三角自卷積窗N=512和N=1024）進行諧波參數估計。以幅值微弱的第21次諧波分析為例，圖3～圖5分別給出了白噪聲存在情況下幅值、初相角、頻率的絕對誤差曲線。

參見圖3～圖5，白噪聲對相角估計準確度存在一定的影響。當信噪比＜50dB時，采用 Blackman窗、Blackman-Harris窗和四項三階Nuttall窗均存在較大誤差；當信噪比＞50dB時，采用Blackman窗、Rife-Vincent（I）窗和 Nuttall窗獲得的幅值、初相角、頻率分析準確度有所提高，采用本文算法的準確度最高。由仿真結果可見，采用本文算法可以有效減少白噪聲對諧波分析準確度的影響，如在信噪比水平為[50dB，100dB]范圍時，采用本文算法可實現諧波參數的準確分析。

圖3 白噪聲情況下第21次諧波幅值分析絕對誤差Fig.3 Absolute errors of amplitude with white noise of the 21st harmonic

圖4 白噪聲情況下第21次諧波初相角分析絕對誤差Fig.4 Absolute errors of phase with white noise of the 21st harmonic

圖5 白噪聲情況下第21次諧波頻率分析絕對誤差Fig.5 Absolute errors of frequency with white noise of the 21st harmonic

4.3 基波頻率變動下的諧波分析

實際應用中，電網基波頻率并非恒定，往往存在小幅波動。設基波頻率波動范圍為[49.5Hz，50.5Hz]范圍，步長為 0.1Hz，采用本文算法（4階三角自卷積窗N=512）對式（20）所示的信號進行諧波參數分析。圖6和圖7分別給出了基波頻率波動情況下基波和各次諧波幅值和初相角的絕對誤差分布情況。

圖6 基波頻率波動時基波與諧波幅值分析絕對誤差Fig.6 Absolute errors of amplitude with frequency changing of fundamental and harmonic components

圖7 基波頻率波動時基波與諧波初相角分析絕對誤差Fig.7 Absolute errors of phase with frequency changing of fundamental and harmonic components

由圖6和圖7可知，當基波頻率發生波動時，諧波參數分析的準確度受到一定程度的影響。圖 6和圖7中，第2次諧波參數分析的準確度受基波頻率波動影響最大，主要原因是第2次諧波幅值較弱，且臨近幅值最大的基波分量，受頻譜泄漏和柵欄效應影響較大。仿真結果表明，在基波頻率波動情況下，采用本文算法實現較為準確的諧波參數分析，基波和各次諧波幅值分析絕對誤差＜0.00005V，初相角分析絕對誤差＜0.075°，可滿足實際測量要求。

5 結論

針對非同步采樣時離散頻譜校正中存在的計算準確度與實時性的矛盾，本文利用三角自卷積窗所具有良好的頻譜泄漏抑制特性，結合三角自卷積窗的頻譜函數，采用最小二乘法建立諧波參數求解的多項式，構建了計算簡便、易于實現的諧波幅值、頻率、初相角參數計算式。在實際使用中，可根據不同精度要求選擇擬合多項式，因此該方法具有可調節、實時性高、易于嵌入式系統實現的優點。對存在白噪聲和基波頻率波動等情況的諧波參數仿真分析結果表明，本文算法能有效克服白噪聲、基波頻率波動的影響，實現高準確度的諧波參數分析，具有一定的實用價值和參考意義。

[1]Chang G W,Chen C I,Liu Y J,et al.Measuring power system harmonics and interharmonics by an improved fast Fourier transform-based algorithm[J].IET Generation,Transmission & Distribution,2008,2(2): 193-201.

[2]Hagh M T,Taghizadeh H,Razi K.Harmonic minimization in multilevel inverters using modified species-based particle swarm optimization[J].IEEE Transactions on Power Electric,2010,24(10): 2259-2267.

[3]蔡濤,段善旭,劉方銳.基于實值MUSIC算法的電力諧波分析方法[J].電工技術學報,2009,24(12):149-155.Cai Tao,Duan Shanxu,Liu Fangrui.Power harmonic analysis based on real-valued spectral MUSIC algorithm[J].Transactions of China Electrotechnical Society,2009,24(12): 149-155.

[4]Yang Rengang,Xue Hui.A novel algorithm for accurate frequency measurement using transformed consecutive points of DFT[J].IEEE Transactions on Power Systems,2008,23(3): 1057-1062.

[5]Ferrero A,Ottoboni R.High-accuracy Fourier analysis based on synchronous sampling techniques[J].IEEE Transactions on Instrumentation and Measurement,1992,41(6): 780-786.

[6]龐浩,李東霞,俎云霄,等.應用 FFT進行電力系統諧波分析的改進算法[J].中國電機工程學報,2003,23(6): 50-54.Pang Hao,Li Dongxia,Zu Yunxiao,et a1.An improved algorithm for harmonic analysis of power system using FFT technique[J].Proceedings of the CSEE,2003,23(6): 50-54.

[7]Chen K F,Li Y F.Combining the Hanning windowed interpolated FFT in both directions[J].Computer Physics Communications,2008,178(12): 924-928.

[8]Agrez D.Dynamics of frequency estimation in the frequency domain[J]. IEEE Transactions on Instrumentation and Measurement,2007,56(6):2111-2118.

[9]Harris F J.On the use of windows for harmonic analysis with the discrete Fourier transform[J].Proceedings of IEEE,1978,66(1): 51-83.

[10]曾博,滕召勝,高云鵬,等.基于Rife-Vincent窗的高準確度電力諧波相量計算方法[J].電工技術學報,2009,24(8): 154-159.Zeng Bo,Teng Zhaosheng,Gao Yunpeng,et al.An accurate approach for power harmonic phasor calculation based on rife-vincent window[J].Transactions of China Electrotechnical Society,2009,24(8): 154-159.

[11]Nuttall A H.Some windows with very good sidelobe behavior [J].IEEE Transactions on Acoustics Speech Signal Processing,1981,29(1): 84-91.

[12]卿柏元,滕召勝,高云鵬,等.基于 Nuttall窗雙峰插值FFT的高精度電力諧波分析方法[J].中國電機工程學報,2008,28(2): 48-52.Qing Baiyuan,Teng Zhaosheng,Gao Yunpeng,et a1.An approach for electrical harmonic analysis based on Nuttall window double-spectrum-line interpolation FFT[J].Proceedings of the CSEE,2008,28(2): 48-52.

[13]張介秋,梁昌洪,陳碩圃.基于卷積窗的電力系統諧波理論分析與算法[J].中國電機工程學報,2004,24(11): 48-52.Zhang Jieqiu,Liang Changhong,Chen Yanpu.Power system harmonic theory analysis and algorithm based on convolution windows[J].Proceedings of the CSEE,2004,24(11): 48-52.

[14]黃純,江亞群.諧波分析的加窗插值改進算法[J].中國電機工程學報,2005,25(15): 26-32.Huang Chun,Jiang Yaqun.Improved window and interpolation algorithm for analysis of power system harmonics[J].Proceedings of the CSEE,2005,25(15):26-32.

[15]溫和,滕召勝,王一,等.基于三角自卷積窗的介損角高精度測量算法[J].電工技術學報,2009,24(2): 164-169.Wen He,Teng Zhaosheng,Wang Yi,et al.High accuracy dielectric loss angle measurement algorithm based on triangular self-convolution window[J].Transactions of China Electrotechnical Society,2009,24(2): 164-169.

[16]Wen H,Teng Z S,Guo S Y,et al.Triangular self-convolution window with desirable sidelobe behaviors for harmonic analysis of power system[J].IEEE Transactions on Instrumentation and Measurement,2010,59(3): 543-552.

[17]Courses E,Surveys T.Dynamics of frequency estimation in the frequency domain[J].IEEE Transactions on Instrumentation and Measurement,2007,56(6): 2111-2118.

[18]Antoni J,Schoukens J.A comprehensive study of the bias and variance of frequency-response-function measurements: optimal window selection and overlapping strategies[J].Automatica,2007,43(10):1723-1736.

[19]Agre D.Interpolation in the frequency domain to improve phase measurement[J].Measurement,2008,41(2): 151-159.

[20]Belega D,Dallet D.Amplitude estimation by a multipoint interpolated DFT approach [J].IEEE Transactions on Instrumentation and Measurment,2009,58(5): 1316-1323.