改進(jìn)加窗頻譜峰值擬合算法及諧波分析應(yīng)用

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

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

1 引言

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

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

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

2 頻譜泄漏與三角自卷積窗函數(shù)

2.1 非同步采樣時(shí)的頻譜泄漏和柵欄效應(yīng)

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

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

信號在觀測周期Tw轉(zhuǎn)換為 N點(diǎn)離散序列相當(dāng)于被矩形窗wR(n)截?cái)?，其離散傅里葉變換為

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

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

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

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

式中，n∈[0,N?1]；m∈[0,M?1]；p為參與卷積的三角窗的個(gè)數(shù)，稱為窗的階數(shù)。為便于計(jì)算，卷積運(yùn)算后可對序列進(jìn)行補(bǔ)零，將p階三角自卷積窗的長度調(diào)整為N=pM。

根據(jù)卷積定理，函數(shù)在時(shí)域進(jìn)行卷積等效于頻域相乘，因此p階三角自卷積窗的頻率響應(yīng)為

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

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

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

表1 窗函數(shù)性能比較Tab.1 Comparison of window characteristics

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

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

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

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

不失一般性，設(shè)當(dāng)前測量的為第h項(xiàng)諧波（h≤H），為簡單起見，忽略其余各次諧波對第h項(xiàng)諧波的泄漏影響，此時(shí)，式（6）變?yōu)?/p>

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

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

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

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

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

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

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

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

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

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

3.3 精度可控的多項(xiàng)式擬合與諧波參數(shù)求解算法

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

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

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

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

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

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

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

由于第h項(xiàng)諧波對應(yīng)的真實(shí)頻率譜線位置應(yīng)為第h k0根，因此第h項(xiàng)諧波的頻率為

由式（10）可得，第h項(xiàng)諧波的幅值計(jì)算式為

同理，可得第h項(xiàng)諧波的初相角為

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

4 仿真與分析

4.1 復(fù)雜信號的諧波分析

采用文獻(xiàn)[12]給出的含 2～21次諧波的信號進(jìn)行分析

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

表2 復(fù)雜信號的基波及諧波參數(shù)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階三角自卷積窗進(jìn)行加權(quán)，利用本文提出的多項(xiàng)式擬合算法進(jìn)行離散頻譜校正后，N=512時(shí)頻率計(jì)算的相對百分比誤差為2.9×10?7%，N=1024時(shí)頻率計(jì)算的相對百分比誤差為-6.4×10?9%。當(dāng) N=1024時(shí)，幅值相對誤差EAh≤0.00094%，相位相對誤差Eφh＜0.012%，可實(shí)現(xiàn)復(fù)雜諧波信號參數(shù)的高準(zhǔn)確度分析。與Blackman窗、Blackman-Harris窗和 Nuttall窗相比，基于 4階三角自卷積窗的最小二乘多項(xiàng)式擬合頻譜校正算法具有更高的諧波分析準(zhǔn)確度，且計(jì)算簡單，易于實(shí)現(xiàn)。

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

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

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

圖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 基波頻率變動下的諧波分析

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

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

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

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

5 結(jié)論

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

[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]蔡濤,段善旭,劉方銳.基于實(shí)值MUSIC算法的電力諧波分析方法[J].電工技術(shù)學(xué)報(bào),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]龐浩,李東霞,俎云霄,等.應(yīng)用 FFT進(jìn)行電力系統(tǒng)諧波分析的改進(jìn)算法[J].中國電機(jī)工程學(xué)報(bào),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窗的高準(zhǔn)確度電力諧波相量計(jì)算方法[J].電工技術(shù)學(xué)報(bào),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].中國電機(jī)工程學(xué)報(bào),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]張介秋,梁昌洪,陳碩圃.基于卷積窗的電力系統(tǒng)諧波理論分析與算法[J].中國電機(jī)工程學(xué)報(bào),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ìn)算法[J].中國電機(jī)工程學(xué)報(bào),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].電工技術(shù)學(xué)報(bào),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.