周 洲

(成都電業(yè)局，四川 成都 610021)

電網(wǎng)中大量非線性負(fù)荷的應(yīng)用導(dǎo)致諧波污染日益嚴(yán)重，對(duì)諧波成分的準(zhǔn)確分析有利于對(duì)電能質(zhì)量進(jìn)行合理評(píng)估。離散傅里葉變換(discrete Fourier transform，DFT)和快速傅里葉變換(fast Fourier transform，F(xiàn)FT)因適用于計(jì)算機(jī)計(jì)算而被廣泛應(yīng)用，但是電網(wǎng)頻率波動(dòng)會(huì)導(dǎo)致信號(hào)采樣的不同步，此時(shí)會(huì)出現(xiàn)影響測(cè)量精度的頻譜泄漏問(wèn)題［1］。

很多種方法已經(jīng)可以降低由于非同步采樣中頻譜泄漏引起的測(cè)量誤差。在頻域，可以選擇適當(dāng)?shù)拇昂瘮?shù)來(lái)降低旁瓣帶來(lái)的長(zhǎng)范圍泄漏，并通過(guò)對(duì)頻率、相位和幅值進(jìn)行插值的手段減少短范圍泄漏［2-10］。

在時(shí)域，可以在截?cái)鄶?shù)據(jù)后補(bǔ)零來(lái)降低柵欄效應(yīng)［1］，通過(guò)波峰檢測(cè)分析諧波參數(shù)，但此方法克服不了諧波間的頻譜干擾且計(jì)算量較大;時(shí)域同步插值通過(guò)對(duì)采樣序列進(jìn)行拉格朗日插值近似實(shí)現(xiàn)序列的同步采樣［11-14］，此方法需要知道信號(hào)的實(shí)際頻率。文獻(xiàn)［15］提出了一種基于相關(guān)計(jì)算的采樣窗口同步化方法，該方法無(wú)需知道信號(hào)的實(shí)際頻率，通過(guò)序列迭代計(jì)算使得采樣序列長(zhǎng)度自動(dòng)滿足同步化收斂條件。但是在實(shí)際應(yīng)用于穩(wěn)態(tài)諧波計(jì)算時(shí)，由于該方法采樣周期在迭代過(guò)程中保持恒定，且實(shí)際系統(tǒng)中電網(wǎng)頻率在一定范圍波動(dòng)，這樣會(huì)使得算法迭代時(shí)間長(zhǎng)，同步化序列長(zhǎng)，并且最終的同步化序列無(wú)法采用FFT算法。因此，提出了一種新的采樣周期自適應(yīng)方法用于諧波計(jì)算，通過(guò)計(jì)算固定長(zhǎng)度相鄰信號(hào)序列的相關(guān)系數(shù)并結(jié)合提出的二次修正方程，不斷迭代調(diào)整采樣周期，得到滿足FFT算法要求的同步化采樣序列。最后通過(guò)FFT計(jì)算得到各次諧波幅值相位。仿真算例證實(shí)了該方法的有效性。

1 基于序列相關(guān)計(jì)算的采樣窗口同步化方法［15］

連續(xù)信號(hào)等間距離散化后得到有限長(zhǎng)的離散時(shí)間序列，表示為信號(hào)向量 S={S0，S1，…SN-1}。兩個(gè)長(zhǎng)度相同的信號(hào)向量S1和S2的相關(guān)計(jì)算為這兩個(gè)向量的內(nèi)積。

其中，θ ?［0，π］，表示兩向量之間的夾角。

定義兩離散信號(hào)序列的相關(guān)系數(shù)如下。

Kp的取值范圍為［-1，1］，Kp值用來(lái)指示兩個(gè)離散信號(hào)之間的相似程度。Kp=1表示兩個(gè)離散信號(hào)相位一致;Kp=-1表示兩個(gè)離散信號(hào)相位相反。

基于相關(guān)系數(shù)的采樣窗口同步化計(jì)算步驟如下。

考慮一離散信號(hào)S

1)選擇p的初始值，確定信號(hào)Xp為迭代的初始序列，其中 Xp={S0，S1，…Sp-1}。

2)定義一個(gè)與Xp長(zhǎng)度相同的后繼信號(hào)Yp={Sp，Sp+1，…S2p-1}。

3)計(jì)算Xp與Yp的相關(guān)系數(shù)Kp。

4)視察Kp是否滿足收斂條件(Kp=1)，若不滿足，令p=p+1，重復(fù)步驟1)、2)和3)直到序列滿足收斂條件。

在此也可以稱Kp為周期指示器。可以通過(guò)考察Kp與1的接近程度來(lái)估計(jì)采樣窗口的同步化程度。當(dāng)Kp=1時(shí)，表示序列Xp和Yp完全相同，采樣已經(jīng)實(shí)現(xiàn)了同步化;對(duì)序列Xp進(jìn)行DFT便可以精確得到信號(hào)各次諧波參數(shù)，此時(shí)頻譜泄漏的影響為零。Kp≠1時(shí)則繼續(xù)進(jìn)行迭代計(jì)算直到滿足收斂條件為止。在實(shí)際應(yīng)用中，若要求對(duì)各種信號(hào)都滿足收斂條件Kp=1，則會(huì)使得采樣序列非常長(zhǎng)。一般情況下根據(jù)精度要求以ε=|Kp-1|≤εset作為收斂條件。

該方法的優(yōu)點(diǎn)在于，不需要預(yù)先知道信號(hào)各成分的頻率，采樣窗口可以通過(guò)基于相關(guān)計(jì)算的迭代方法自適應(yīng)求得，從而實(shí)現(xiàn)信號(hào)的整周期截?cái)嗯c同步采樣。很顯然，此方法在信號(hào)非平穩(wěn)的情況下是不適用的。往往滿足收斂條件的同步化序列非常長(zhǎng)，在固定采樣周期的非同步采樣下尤為明顯。

2 基于相關(guān)系數(shù)的采樣周期自適應(yīng)

上節(jié)所述的信號(hào)周期自適應(yīng)方法可以在不知道信號(hào)頻率的情況下對(duì)其實(shí)現(xiàn)同步采樣，但是其存在下面兩個(gè)問(wèn)題:①要滿足收斂條件ε=| Kp-1|≤εset的采樣序列較長(zhǎng)，迭代運(yùn)算次數(shù)較多且DFT的運(yùn)算量較大;②通過(guò)迭代計(jì)算得到采樣序列的點(diǎn)數(shù)pmax無(wú)法保證滿足pmax=2M(M為正整數(shù))，以至在一般情況下對(duì)采樣序列不能使用FFT算法。

現(xiàn)考慮固定每周期采樣點(diǎn)數(shù)N=2M。初始采樣周期為T(mén)s0=Nf0(f0表示信號(hào)沒(méi)有頻偏時(shí)的頻率，對(duì)于50 Hz系統(tǒng)，f0=50 Hz)。

通過(guò)計(jì)算XN與YN(YN為XN的后繼N點(diǎn)序列)的相關(guān)系數(shù)，根據(jù)修正方程，不斷地修正信號(hào)的采樣周期，使得信號(hào)滿足收斂條件實(shí)現(xiàn)同步采樣。

計(jì)算中相關(guān)系數(shù)越是接近1，表示此時(shí)的XN與YN越接近，采樣周期的修正量就應(yīng)該越小;反之，相關(guān)系數(shù)越是遠(yuǎn)離1，表示XN與YN相似程度越差，采樣周期的修正量就應(yīng)該越大。這樣可以加快信號(hào)的收斂速度。

設(shè)第i+1次采樣周期修正為

ΔNi稱為采樣周期修正因子，i=0，1，2…。

根據(jù)上述要求，采用二次拋物線的修正方程為

其中，AN為迭代加速因子;N為采樣點(diǎn)數(shù);KNi表示第i-1次序列XNi與YNi的相關(guān)系數(shù)。

當(dāng)N=64時(shí)修正方程的函數(shù)圖形如圖1所示。

圖1 采樣周期修正因子ΔNi的修正方程函數(shù)圖

從圖1可以看出，相關(guān)系數(shù)越是接近1，修正量就越小，這樣使得收斂更快且保證信號(hào)序列由快到慢接近收斂點(diǎn)。在某種程度上來(lái)講，加速因子AN的取值越大，收斂速度越快。

基于相關(guān)系數(shù)的采樣周期自適應(yīng)具體步驟如下。

1)根據(jù)Ts0采樣得到2N點(diǎn)的信號(hào)序列，前N個(gè)點(diǎn)序列為XN0，后N個(gè)點(diǎn)序列為YN0，根據(jù)式(2)計(jì)算得到相關(guān)系數(shù)KN0。

2)根據(jù)修正方程ΔN0=AN和式(3)分別計(jì)算

3)根據(jù)步驟2)分別得到兩個(gè)采樣周期值，采樣并計(jì)算KN1和KN2。若KN1＞KN2，說(shuō)明信號(hào)頻率偏大，式(3)分母中的符號(hào)在迭代過(guò)程中定為“+”號(hào);若KN1＜KN2，說(shuō)明信號(hào)頻率偏小，定為“-”號(hào)。

4)根據(jù)步驟2)得到的采樣周期對(duì)信號(hào)繼續(xù)采樣，計(jì)算 KNi。

5)根據(jù)修正方程ΔNi=AN，和式(3)計(jì)算得到第i+1次的采樣周期。

6)根據(jù)新的采樣周期繼續(xù)采樣2N點(diǎn)的信號(hào)序列，計(jì)算KNi，判斷是否滿足收斂條件。滿足則結(jié)束計(jì)算，不滿足則返回步驟5)繼續(xù)迭代計(jì)算，直到滿足為止。

最后，對(duì)滿足收斂條件的最新采樣序列做FFT計(jì)算得到信號(hào)各次諧波的幅值和相對(duì)相位(由于采樣周期的不斷調(diào)整使得各次諧波相位發(fā)生變化，相位計(jì)算采用以基波相位作為零度基準(zhǔn)值)。

3 數(shù)字仿真計(jì)算

含有各次諧波的電力系統(tǒng)信號(hào)表示為

其中，m為正整數(shù);Am和φm分別表示m次諧波的幅值和相角，具體參數(shù)如表1所示。

表1 諧波參數(shù)設(shè)置

計(jì)算1:對(duì)兩種方法的計(jì)算量進(jìn)行比較。

首先，對(duì)信號(hào)進(jìn)行基于采樣窗口同步化的迭代計(jì)算［15］，采樣頻率為固定值3200 Hz，初始迭代序列Xp的p=60，收斂標(biāo)準(zhǔn)Kp＞0.99999。在不同頻率偏差的情況下得到收斂所需的迭代次數(shù)、加法和乘法運(yùn)算次數(shù)如表2所示。再對(duì)信號(hào)進(jìn)行基于采樣周期自適應(yīng)的迭代計(jì)算，此時(shí)，初始采樣頻率fs0=3200 Hz，收斂標(biāo)準(zhǔn) Kp＞0.9999999999。每周期采樣點(diǎn)數(shù)恒是N=64，通過(guò)迭代計(jì)算不斷地修正采樣周期最終使得采樣序列達(dá)到收斂條件，具體的運(yùn)算量亦如表2所示。

從表中看出，后者的迭代次數(shù)及運(yùn)算量要明顯小于前者。采用基于采樣周期自適應(yīng)的迭代計(jì)算時(shí)，隨著加速因子的增加可以有效地提高收斂速度。

計(jì)算2:對(duì)諧波計(jì)算精度進(jìn)行數(shù)字比較。

假設(shè)信號(hào)實(shí)際頻率f=49.5 Hz，對(duì)式(4)信號(hào)分別采用以下方法進(jìn)行諧波計(jì)算。方法1，直接采用FFT算法計(jì)算;方法2，采用基于信號(hào)周期自適應(yīng)的方法計(jì)算［15］;方法3，采用基于采樣周期自適應(yīng)的方法計(jì)算。方法1的采樣頻率fs=3200 Hz，窗口長(zhǎng)度為64點(diǎn)。方法2和方法3的采樣頻率、收斂條件等與計(jì)算1時(shí)相同。

計(jì)算結(jié)果如表3所示，從各次諧波的分布可以看出，信號(hào)中的偶次諧波含量要遠(yuǎn)比奇次諧波小。

現(xiàn)令

作為測(cè)量精度指示因子。ΔAm和Δφm為m次諧波幅值和相位的絕對(duì)誤差，由于基波和8次諧波存在零值，在此不參與α因子計(jì)算。最后單獨(dú)做絕對(duì)誤差比較。

表2 兩種方法運(yùn)算量的比較

根據(jù)表中數(shù)據(jù)計(jì)算，采用方法1時(shí)，奇次諧波的α因子為

偶次諧波α因子為

由于非同步采樣，尤其是對(duì)于含有量較小的偶次諧波而言測(cè)量誤差較大。

采用方法2時(shí)，奇次、偶次諧波精度指示因子分別為 α奇次=0.000862221和 α偶次=0.0014。其精度和FFT相比有著明顯地提高。

采用方法3時(shí)，根據(jù)最后一列數(shù)據(jù)，得到奇次、偶次諧波精度指示因子分別為 α奇次=0和α偶次=0.000016766。并且從表中可以看出收斂標(biāo)準(zhǔn)的提高導(dǎo)致測(cè)量精度的提高，但在相同加速因子的條件下會(huì)增加迭代次數(shù)。與方法2相比，增加后的運(yùn)算量仍甚小，并且精度要遠(yuǎn)高于方法2。

對(duì)于基波，方法3的測(cè)量誤差為零值。方法1、2都有一定程度的誤差存在。對(duì)于8次諧波，方法3計(jì)算得到的幅值誤差要遠(yuǎn)遠(yuǎn)小于其他方法，但是其相位誤差較大。

在采用方法2時(shí)，由于信號(hào)的頻率未知，導(dǎo)致滿足收斂標(biāo)準(zhǔn)的信號(hào)序列很長(zhǎng)，且不能保證滿足FFT的采樣點(diǎn)數(shù)要求(如表中的點(diǎn)數(shù)為1099點(diǎn))，這樣通過(guò)DFT計(jì)算諧波參數(shù)時(shí)計(jì)算量就會(huì)更加大。

4 應(yīng)用討論

4.1 加速因子的取值范圍

加速因子AN的大小決定了計(jì)算迭代的初始步長(zhǎng)。從這里的仿真結(jié)果可以看出，加速因子AN的取值越大，迭代計(jì)算的初始步長(zhǎng)越長(zhǎng)，算法收斂速度也就越快，諧波計(jì)算的實(shí)時(shí)性也就越高。但是，對(duì)于一個(gè)工頻周波信號(hào)，若盲目提高加速因子取值，可能會(huì)引發(fā)計(jì)算收斂序列變成若干個(gè)周期，從而影響計(jì)算實(shí)時(shí)性和DFT計(jì)算量。因此，加速因子AN應(yīng)根據(jù)實(shí)際信號(hào)的頻偏特性在一起范圍內(nèi)取值。這個(gè)取值范圍取決于信號(hào)頻率偏差范圍。

設(shè)信號(hào)的頻率偏差為Δf'，理想的采樣周期應(yīng)為

表3 3種方法諧波參數(shù)計(jì)算結(jié)果

根據(jù)式(8)可以計(jì)算得到

在采樣周期自適應(yīng)的迭代過(guò)程中，步長(zhǎng)ΔNi應(yīng)小于等于ΔN'。另由于迭代步長(zhǎng)是單調(diào)遞減的，因此第一次迭代計(jì)算的值應(yīng)為最大ΔNmax，必須滿足

得到AN的取值范圍為

根據(jù)電力系統(tǒng)的運(yùn)行規(guī)范和標(biāo)準(zhǔn)，電力系統(tǒng)的頻率往往被控制在50 Hz±0.5 Hz之間，因此信號(hào)的頻率偏差是有上限的，根據(jù)這個(gè)上限，可以確定加速因子AN的取值上限。取值為0.5 Hz。

4.2 收斂條件

計(jì)算結(jié)果的精度取決于用戶所定的計(jì)算收斂條件，即為:ε=| Kp-1|≤εset。若Kp取值越高，則計(jì)算誤差越小;反之，則計(jì)算誤差越大。但提高計(jì)算精度的代價(jià)是影響計(jì)算的實(shí)時(shí)性，因此，在實(shí)際應(yīng)用過(guò)程中，Kp的取值可以根據(jù)測(cè)量的精度要求、計(jì)算實(shí)時(shí)性和環(huán)境的噪聲水平等綜合考慮后取值，在先保證信號(hào)實(shí)時(shí)性的條件下，即采樣的時(shí)間必須小于某個(gè)數(shù)tmax的條件下，取相關(guān)系數(shù)最高的那組周波數(shù)據(jù)進(jìn)行DFT計(jì)算，得到相應(yīng)的諧波值。

5 結(jié)論

采用采樣周期自適應(yīng)的方法可以有效解決由于非同步采樣引起的頻譜泄露問(wèn)題。具體有以下幾個(gè)特點(diǎn)。

1)該方法無(wú)需知道電網(wǎng)信號(hào)的實(shí)際頻率，隨著電網(wǎng)頻率的偏移，迭代計(jì)算可以不斷適應(yīng)頻率的動(dòng)態(tài)偏移實(shí)現(xiàn)同步采樣。

2)迭代過(guò)程中通過(guò)對(duì)采樣周期進(jìn)行修正，并引入修正加速因子可以提高信號(hào)序列收斂速度，即使在頻率偏差較大時(shí)也能即時(shí)使序列很快收斂。

3)在迭代過(guò)程中采樣點(diǎn)數(shù)恒定，每次迭代的計(jì)算量相同;而采樣周期固定時(shí)，每次迭代的計(jì)算量隨著迭代次數(shù)的增加而增加。當(dāng)?shù)螖?shù)較多時(shí)計(jì)算量很大。

4)采樣點(diǎn)數(shù)滿足N=2M，因此可以采用FFT計(jì)算，提高運(yùn)算速度。

5)實(shí)際分析時(shí)也可用時(shí)域多項(xiàng)式插值手段對(duì)采樣信號(hào)進(jìn)行采樣周期調(diào)整的重采樣。

［1］Ignacio，Santamaria-Caballero，Carlos J，Pantaleon-Prieto，Jesus Ibanez，et al.Improved Procedures for Estimating Amplitudes and Phases of Harmonics with Application to Vibration Analysis［J］.IEEE Transactions on Instrumentation and Mesurement，1998，47(1):209-2144.

［2］Dusan Agrez.Weighted Multipoint Interpolated DFT to Improve Amplitude Estimation of Multifrequency Signal［J］.IEEE Transactions on Instrumentation and Mesurement，2002，51(2):287-292.

［3］D.Belega，D.Dallet.Frequency Estimation via Weighted Multipoint Interpolated DFT［J］.IET Science，Measurement＆ Technology，2008，2(1):1-8.

［4］潘文，錢(qián)俞壽，周鶚.基于加窗差值FFT的電力諧波測(cè)量理論(I)窗函數(shù)研究［J］.電工技術(shù)學(xué)報(bào)，1994，9(1):50-54.

［5］潘文，錢(qián)俞壽，周鶚.基于加窗差值FFT的電力諧波測(cè)量理論(II)雙插值 FFT理論［J］.電工技術(shù)學(xué)報(bào)，1994，9(2):53-56.

［6］張伏生，耿中行，葛耀中.電力系統(tǒng)諧波分析的高精度FFT 算法［J］.中國(guó)電機(jī)工程學(xué)報(bào)，1999，19(3):63-66.

［7］孫紅偉，李梅，袁健華，等.用于電力系統(tǒng)諧波分析的加窗FFT 算法研究［J］. 高電壓技術(shù)，2004，30(8):52-55.

［8］錢(qián)昊，趙榮祥.基于插值FFT算法的間諧波分析［J］.中國(guó)電機(jī)工程學(xué)報(bào)，2005，25(21):87-91.

［9］趙文春，馬偉明，胡安.電機(jī)測(cè)試中諧波分析的高精度FFT算法［J］.中國(guó)電機(jī)工程學(xué)報(bào)2001，21(12):83-87.

［10］柴旭崢，文習(xí)山，關(guān)根志，等.一種高精度的電力系統(tǒng)諧波分析算法［J］.中國(guó)電機(jī)工程學(xué)報(bào)，2003，23(9):67-70.

［11］陳隆道，錢(qián)照明，張圣訓(xùn).周期域分析中非同步取樣數(shù)據(jù)的同步化［J］.電子學(xué)報(bào)，2001，29(7):1-4.

［12］曾澤昊，余有靈，許維勝.一種減小泄漏的同步化算法［J］.電測(cè)與儀表，2005(479):12-14.

［13］曾澤昊，余有靈，許維勝.基于插值同步算法的頻譜泄漏分析與仿真［J］.系統(tǒng)仿真技術(shù)，2005，1(3):148-152.

［14］惠錦，楊洪耕.改進(jìn)的時(shí)域同步插值算法在基于TDA諧波和間諧波測(cè)量中的應(yīng)用［J］.電工技術(shù)學(xué)報(bào)，2008，23(7):124-129.

［15］T.X.Zhu.Exact Harmonics/Imterharmonics Calculation Using Adaptive Window Width［J］.IEEE Transactions on Power Delivery，2007，47(1):2279-2288.

［16］鄭恩讓?zhuān)瑮顫?rùn)賢，高森.關(guān)于電力系統(tǒng)FFT諧波檢測(cè)存在問(wèn)題的研究［J］. 繼電器，2006，34(18):52-57.

［17］布賴母，柳群譯.快速傅立葉變換［M］.上海:上海科學(xué)技術(shù)出版社，1979.