基于一階自卷積Kaiser窗的分頻段諧波檢測DFT算法

謝業強，王鋼，曾德輝

(1.華南理工大學 電力學院，廣東 廣州 510641；2.廣州嘉緣電力科技有限公司，廣東 廣州 510610)

電力諧波檢測是電網諧波得到有效治理的前提。電力電子化電力系統的諧波特性已發生變化，諧波帶寬大幅拓寬以及噪聲干擾加劇給電力諧波研究提出了新的挑戰。目前許多先進的算法被用于諧波檢測，比如小波變換[1-3]、譜估計[4]、希爾伯特黃變換[5-6]以及神經網絡[7]。小波變換是時-頻算法，具有良好的時域特性，但當諧波成分過多時會造成嚴重的混疊現象，不適合弱諧波分量較多的寬頻信號檢測；譜估計算法理論上具有無限的分辨率，但該算法受噪聲干擾嚴重，計算量大且易出現虛假分量；希爾伯特黃變換具有良好的暫態特性，但存在嚴重的模態混疊現象，同樣不適合于寬頻信號檢測；神經網絡具有很強的自適應性，但檢測精度十分依賴于初始樣本。因此，綜合考慮帶寬信號所要求的處理效率和檢測精度，上述方法都不能滿足電力電子化條件下的電力諧波檢測要求。離散傅里葉變換(discrete Fourier transformation,DFT)[8-10]易于寫入嵌入式系統且具有計算量小、效率高、操作簡單等優點，然而DFT存在固有的性能缺陷，如頻譜泄漏和柵欄效應。

目前，針對DFT諧波檢測算法固有性能缺陷的研究，已取得一定成果。為了降低頻譜能量泄漏，最常見的策略是采用窗函數對信號進行加權，(比如三角窗[11]、平頂窗[12]、漢寧窗[13-14]、Blackman窗[15]、Nuttall窗[16]等經典窗函數)，但是經典窗函數的主瓣寬度是固定的，無法同時保證旁瓣峰值低和快速的旁瓣衰減特性，性能有待提高。與經典窗相比，Kaiser窗[17-18]能夠通過改變調節因子使主瓣能量與旁瓣能量之比近乎最大，且能夠自主選擇主瓣寬度和旁瓣高度間的比重[19]，因此信號加權更靈活。另外，通過對窗函數進行時域自乘運算[20]或卷積運算[21]可進一步提高窗函數抑制頻譜泄漏性能，但時域自乘運算會導致窗函數的主瓣寬度變寬[22]，降低頻率分辨率而不利于寬頻信號檢測。為了減少DFT柵欄效應引起的理論誤差，可借助校正措施給予修正，如譜線插值算法[23]。根據緊鄰頻域峰值點加權譜線的數量，可分為雙譜線、三譜線、四譜線插值算法[24]等，加權譜線數量越多校正效果越明顯，但加權譜線的增多也會導致計算量增加[25]。當前，校正算法須在頻率估計值已知的情況下，才可求解相位估計值，這會將頻率估計誤差累計到相位測量中，造成測相精度受到影響，在噪聲環境下更加顯著。文獻[26]提出一種名為全相位快速傅里葉變換算法(all phase fast Fourier transformation,APFFT)，其原理可簡單概括為在進行DFT算法前對原始數據進行全相位數據預處理，使得APFFT具有相位不變的優良特性，理論上可以借助精確的相位差值去校正頻率和幅值；但全相位數據預處理會導致柵欄效應加劇，使得頻率與幅值測量存在較大的理論誤差，尤其是弱諧波分量。

對此，為綜合提高寬頻諧波信號檢測的精度、穩定性和抗噪性能，在研究窗函數自卷積運算、譜線插值算法以及全相位數據預處理算法優缺點的基礎上，實現優勢互補，提出基于一階自卷積Kaiser窗的諧波檢測DFT算法。為降低寬頻信號的檢測時間，對寬頻信號進行分頻段處理，提出分頻段的諧波檢測DFT算法，有效地降低高頻段的檢測時間，在保證檢測精度下提高計算速度。所提算法通過頻段劃分能有效解決寬頻信號中高頻段諧波檢測計算時間冗余問題，以及通過自卷積運算實現全相位數據預處理算法和譜線插值算法的優勢互補，綜合增強寬頻諧波檢測算法性能。最后，構建包含50倍基頻多正弦復雜信號模型，并引入基波頻率波動和白噪聲干擾測試驗證所提算法對寬頻諧波檢測的可行性。

1 基于一階自卷積Kaiser窗的改進算法

加窗DFT算法經過譜線插值校正以后具有良好的頻率和幅值檢測精度，但存在頻率估計誤差蔓延問題，進而影響到相位檢測的精度。通過全相位數據預處理的DFT算法具有相位不變特性，能夠解決譜線插值DFT算法的頻率誤差蔓延問題；然而，原始數據經過全相位數據預處理算法后，數據長度變長即由原來的N個點變成2N-1個點，不能與DFT算法直接融合實現互補，而通過對窗函數進行自卷積運算能解決上述問題。故須研究Kaiser窗函數的頻譜泄漏抑制性能以及自卷積運算對Kaiser窗的影響。

1.1 Kaiser窗

Kaiser窗的時域表達式[27]

(1)

式中：n=0,1,…,N-1為采樣序列；N為窗函數長度；I0(β)為零階貝塞爾函數，其中β為調節因子，函數表達式為

(2)

長度為N的Kaiser窗頻域表達式

(3)

式中ω=2πk/N，k為譜線位置。

對Kaiser窗進行一階自卷積運算并歸一化得2N-1點的時域歸一化卷積窗

(4)

式中w(n)*w(n)為一階自卷積運算。

歸一化卷積窗的頻域表達式

(5)

式中上標*表示自卷積后的窗函數。

圖1給出長度N=25的Kaiser窗的歸一化對數頻譜，其中圖1(a)—(d)分別為β=0、4、8、11時，Kaiser窗的歸一化對數頻譜，其中“dB/oct”表示每2倍頻的分貝變化量。圖2展示一階自卷積Kaiser窗與Kaiser窗的歸一化對數頻譜對比。

圖1 Kaiser窗的歸一化對數頻譜Fig.1 Normalized logarithmic spectrum of Kaiser window

圖2 一階自卷積Kaiser窗與Kaiser窗的頻譜對比Fig.2 Comparisons of spectrum between 1st-order Kaiser self-convolution window and Kaiser window

選擇窗函數時需要考慮被檢測信號的特性。由于電力電子化電力系統的諧波帶寬大幅拓寬以及噪聲干擾嚴重，這不僅要求窗函數具有良好的旁瓣特性，而且還要求具有窄的主瓣寬度。但是對窗函數而言，主瓣寬度和旁瓣性能是固有的矛盾體，即頻率分辨率和頻譜泄漏抑制能力不能同時兼得。由圖1(a)—(d)可知，當Kaiser窗隨著β值的提升，旁瓣峰值與主瓣峰值差值越大且漸進衰減速度越快，具有更高的抑制頻譜泄漏能力；然而，隨著β值的增加，Kaiser窗函數的主瓣寬度越寬，這將導致頻率分辨率的減低，不利于寬頻諧波的檢測。由圖2可知，Kaiser窗在進行一階自卷積運算時，能夠在不改變主瓣寬度的情況下，進一步降低旁瓣峰值和加快衰減速度，從而提高抑制頻譜泄漏能力。此外，一階自卷積Kaiser窗(β=11)的旁瓣性能可以媲美調節因子β=20的Kaiser窗，且具有更窄的主瓣寬度，說明調節因子β=11的Kaiser窗經過一階自卷積改造后，在不影響頻率分辨率的情況下取得了更好的頻譜泄漏抑制能力，足以滿足目前的諧波檢測要求。因此，為平衡頻率分辨率和旁瓣性能，折中選擇調節因子β=11的Kaiser窗作為自卷積運算的母窗函數。

理論上可以通過更高階的自卷積運算提高窗函數的旁瓣性能，但高階的自卷積運算會導致運算數據的增加，比如二階自卷積運算的數據將由原來的N點增加至3N-2個點，而且高階的自卷積運算同樣會導致窗函數主瓣寬度的增加[28]；因此，綜合考慮頻譜泄漏的抑制能力、頻率分辨能力以及全相位數據預處理算法的相位不變特性，選擇一階自卷積Kaiser窗(β=11)對諧波信號進行處理。

1.2 基于一階自卷積Kaiser窗的DFT算法

對已知諧波原始信號s(n)均勻采樣2N-1個點所得到的離散時間信號

(6)

式中：n=0,1,…,2N-2為采樣序列；H為諧波項數；Ah為第h次諧波的幅值；φh為第h次諧波的初始相位；fs為采樣頻率；f1、A1、φ1分別為基波的頻率、幅值、初相。

為抑制頻譜泄漏，對s(n)加歸一化卷積窗w*(n)，并忽略負頻點-hf1處信號的旁瓣干擾，得加窗后的DFT表達式

(7)

式中hf1=khΔf為第h次諧波的頻率，Δf=fs/(2N-1)。

為簡化推導，設第h次諧波為檢測量，忽略頻譜泄漏，故式(7)簡化為

(8)

考慮到加權譜線數量的增多會導致計算量的增加，且已采用具有良好旁瓣性能的卷積窗函數加權抑制頻譜泄漏，峰值譜線附近譜線幅值較小，采用過多加權譜線校正后精度提高不明顯，故引入雙譜線插值算法校正頻率和幅值偏差，具體求解參照文獻[29]。

依據雙譜線插值算法可推導出第h次諧波的頻率估計表達式

(9)

第h次諧波的幅值修正公式

(10)

式(9)、(10)中：L1=|S′(khaΔf)|、L2=|S′(khbΔf)|分別為靠近被測頻點附近的最大譜線kha和次最大譜線khb的幅值，滿足kha≤kh≤khb=kha+1；α=kh-kha-0.5，α∈[-0.5,0.5]。

為解決譜線插值校正算法的頻率誤差蔓延到相位檢測問題，引入全相位數據預處理算法，具體計算過程參照文獻[26]，則式(8)經過全相位數據預處理算法后得到的APFFT頻譜為

(11)

故第h次諧波的初相位

(12)

綜上，式(9)、(10)和(12)構成了高精度穩定的諧波參數估計表達式，可求得各次諧波的頻率、幅值和相位。

2 分頻段諧波檢測DFT算法

針對電力電子化電力系統諧波帶寬情況，若各頻段諧波的檢測均采用同樣的采樣頻率和采樣窗寬數據，勢必增加諧波檢測的計算負擔，故提出分頻段諧波檢測思想，即針對不同頻段的諧波采用不同的采樣頻率和采樣窗寬對信號進行一階自卷積Kaiser窗的DFT算法分析。

本文采用的一階自卷積Kaiser窗的分頻段DFT算法的分析計算過程如下：

a)對原始諧波信號s(n)進行Nf個頻段劃分，并設置各頻段采樣頻率fsi和采樣窗寬dti，其中i∈Nf，i=1時為劃分的最低頻段。

b)利用一階自卷積Kaiser窗的DFT算法，分析計算最低頻段(i=1)的各次諧波分量的頻率f′ij、幅值A′ij和相位φ′ij，其中j為各頻段的諧波次數；

c)依據已分析頻段的諧波分量信息，進行信號的時域重構：

(13)

d)求取下一頻段(i+1頻段)待檢測信號的殘余分量

Δsi(n)=s(n)-s′(n).

(14)

e)對Δsi(n)進行一階自卷積Kaiser窗的DFT計算，求得相應頻段內的諧波分量的頻率f′ij、幅值A′ij和相位φ′ij。

f)重復過程c)—e)，直至滿足fmax∈fim，fim為第i頻段的理論最高頻率，fmax=max(f′ij)為第i頻段檢測結果中的最高頻率。

一階自卷積Kaiser窗的分頻段DFT算法流程如圖3所示。

圖3 分頻段DFT算法流程Fig.3 Segmented DFT algorithm flow

由于電力電子化電力系統諧波檢測帶寬的要求已經擴展到50倍基頻，且實際電網中以低頻諧波成分為主[30]，將50次諧波劃分為高低2個頻段；以1 000 Hz作為高低2個頻段的分頻點，故劃分的低頻段為[10 Hz,1 000 Hz]，高頻段為[950 Hz,2 500 Hz]。采樣頻率取fs=20 kHz、采樣時間取1 s，并通過間隔數據二次采樣形式降低各頻段的采樣頻率。低高頻段的劃分及采樣頻率、采樣點數的取值見表1。

表1 頻段劃分Tab.1 Frequency band segmentation

3 仿真分析

3.1 寬頻信號諧波參數估計

為檢驗所提基于一階自卷積Kaiser窗的分頻段DFT算法對寬頻信號諧波參數估計的有效性，構建50倍基頻的多正弦基波和諧波的原始信號，其離散時間信號

n=0,1,…,2N-2.

(15)

式中：基波頻率f1為50.2Hz；各頻段的采樣頻率和采樣點數見表1，給定仿真參數見表2。

表2 給定仿真參數Tab.2 Given simulation parameters

仿真對比3種雙譜線插值算法，分別為基于Kaiser窗(β=11) (算法1)、基于一階自卷積Kaiser窗(β=11)(算法2)以及分頻段一階自卷積Kaiser窗 (算法3)。為確保仿真對比的有效性，算法1和算法2的采樣點數與算法3一致，總共為2×(2N-1)=3 198個點，采樣頻率和原始數據的采樣頻率fs一樣為20 kHz。各次諧波頻率、幅值和相位相對誤差分別如圖4—圖6所示。

由圖4—圖5可知，經過一階自卷積運算的算法2和算法3的頻率和幅值的檢測精度較算法1高，這說明相同條件下通過自卷積運算能夠提高Kaiser窗函數的旁瓣性能，進而提高抑制頻譜泄漏能力。由圖6可知，經過全相位數據預處理的算法2和算法3的相位估計精度總體達到10-7，而算法1的相位估計精度較低，這說明全相位數據預處理算法能夠解決譜線插值算法的頻率誤差累計到相位檢測誤差中的問題。另外由圖4—圖6可知，基于分頻段的算法3總體精度不及算法2，但仍具有較高的檢測精度，這說明在算法2的基礎上對原始寬頻信號進行分頻段的處理是可行的，但會產生一定的理論誤差。

圖4 頻率估計值相對誤差Fig.4 Relative errors of frequencies

圖5 幅值估計值相對誤差Fig.5 Relative errors of amplitudes

圖6 相位估計值相對誤差Fig.6 Relative errors of phases

為了分析所提分頻段的諧波檢測思路對減輕計算負擔的有效性，將算法3的計算時間與算法1和算法2進行比較，計算機處理器為Intel(R)i5-4210M，主頻為2.6 GHz。算法1、算法2和算法3在MATLAB軟件中的計算時間分別為0.039 s、0.095 s、0.079 s。可見，算法2由于包括一階自卷積運算和全相位數據預處理運算，計算時間有所增加。另外，算法3是在算法2的基礎上對寬頻信號進行分頻段處理，計算負擔有所下降。

3.2 基頻波動

當采樣頻率fs恰好是基波頻率f1的整數倍時，實現同步采樣，此時諧波檢測算法具有低的頻譜泄漏和柵欄效應；但是電網基波頻率時刻波動，基波頻率波動會導致非同步采樣，從而加劇頻譜泄漏和柵欄效應，尤其在發生故障時，對諧波參數估計會造成較大影響。對3種算法抗基波波動的性能進行仿真實驗，設基波頻率以0.1 Hz的步差在49.5～50.5 Hz區間內變化[31]，采樣頻率及其他參數同3.1節。據此獲得的基波頻率估計值相對誤差如圖7所示，各次諧波的頻率、幅值及相位的相對誤差分布分別如圖8—圖16所示。

圖7 基波頻率估計值相對誤差Fig.7 Relative errors of fundamental frequencies

圖8 算法1頻率估計值相對誤差Fig.8 Relative errors of frequencies of algorithm 1

由圖7可知，基波頻率發生波動時，3種算法檢測基波都有較高的精度，但是算法2和算法3比算法1的檢測精度高，說明一階自卷積運算提高了窗函數旁瓣性能。算法3是在算法2的基礎上對原始寬頻信號進行分頻段處理，進而對諧波信號進行檢測，由圖7可見算法3的檢測精度略微低于算法2，但高于算法1，表明分頻段處理是可行的且對基波波動的干擾影響不大。

圖9 算法2頻率估計值相對誤差Fig.9 Relative errors of frequencies of algorithm 2

圖10 算法3頻率估計值相對誤差Fig.10 Relative errors of frequencies of algorithm 3

圖11 算法1幅值估計值相對誤差Fig.11 Relative errors of amplitudes of algorithm 1

圖12 算法2幅值估計值相對誤差Fig.12 Relative errors of amplitudes of algorithm 2

由圖8—圖13知，在基波頻率波動的情況下，算法2和算法3的諧波檢測精度較算法1更高。盡管算法3檢測精度較算法2略微降低，但仍然具有良好的抗基波波動干擾性能。

圖13 算法3幅值估計值相對誤差Fig.13 Relative errors of amplitudes of algorithm 3

由圖14—圖16可見，算法1在基波頻率波動時，相位檢測的性能明顯低于其他2種算法，而經過全相位數據預處理算法改進得到的算法2和算法3具有良好的檢測精度；算法3引入分頻段產生的細微理論誤差，使其相位檢測精度略低于算法2。

圖14 算法1相位估計值相對誤差Fig.14 Relative errors of phases of algorithm 1

圖15 算法2相位估計值相對誤差Fig.15 Relative errors of phases of algorithm 2

圖16 算法3相位估計值相對誤差Fig.16 Relative errors of phases of algorithm 3

3.3 白噪聲干擾

實際電網中，噪聲不可避免，且以高斯白噪聲為主，主網信噪比通常為60 ～80 dB，低壓配電網通常為60 dB，最小可達到45 dB[32]。隨著電力系統電力電子化，噪聲愈加復雜，將影響諧波參數的測量；因此，有必要進行更低信噪比的仿真實驗以驗證算法的抗噪性能，故給定信噪比為30 dB的高斯白噪聲環境，對比3種算法的抗噪性能，仿真參數同3.1節，各次諧波頻率、幅值和相位估計值的相對誤差如圖17—19圖所示。

由圖17—圖19可見，在信噪比為30dB的強白噪聲干擾下，經一階自卷積運算的算法2較算法1具有較高的檢測精度，尤其是相位精度；算法3是在算法2的基礎上進行分頻段處理，在噪聲干擾下精度有所下降，但是下降不明顯，仍具有較高的穩定性。

圖17 噪聲背景下頻率相對誤差Fig.17 Relative frequency errors in noisy environment

圖18 噪聲背景下幅值相對誤差Fig.18 Relative errors of amplitudes in noisy environment

圖19 噪聲背景下相位相對誤差Fig.19 Relative errors of phases in noisy environment

4 結論

基于電力電子化電力系統寬頻諧波特性背景，提出基于一階自卷積Kaiser窗的分頻段諧波檢測DFT算法，結論如下：

a)一階自卷積運算能夠在不改變Kaiser窗主瓣寬度的情況下進一步降低旁瓣峰值和加快旁瓣衰減速度。

b)一階自卷積運算會導致運算數據由原來的N個點增加至2N-1個點，但能實現譜線插值算法和全相位數據預處理算法的優勢互補，可以在不明顯增加計算負擔的情況下提供較好的頻率、幅值和相位的檢測精度。

c)分頻段諧波檢測思路對寬頻信號的檢測具有可行性，能在一定程度上降低計算負擔，但會因為引入理論誤差而導致檢測精度略微降低。

d)所提分頻段算法適合寬頻信號檢測，能有效克服基波頻率波動和噪聲干擾，諧波參數檢測精度高且具有較高的穩定性和較低的計算負擔。

本文未能實現Kaiser窗的調整因子β自適應調節，另外沒有給出通用的頻段選擇原則和自動的實現方式。