黨存祿，陳國國

（蘭州理工大學(xué) 電氣工程與信息工程學(xué)院，蘭州 730050）

國家規(guī)范對電力系統(tǒng)的頻率有著嚴(yán)格的要求，一般情況下頻率的偏差值為±0.2 Hz，它反映的是系統(tǒng)有功功率的平衡[1]。風(fēng)力發(fā)電和光伏發(fā)電的裝機容量越來越大，但是在這些新型能源系統(tǒng)中，頻率是不穩(wěn)定的，是隨著時間不斷變化的[2]。例如：負(fù)載波動引起的有功功率不平衡，系統(tǒng)發(fā)生各種故障時，也會引起有功功率的不平衡。波動的頻率不僅會影響到其測量的準(zhǔn)確性，而且有可能會使低頻減載裝置誤動作或拒動[3]。

在新能源系統(tǒng)中，頻率的測量主要由嵌入式微機系統(tǒng)來實現(xiàn)[4]。隨著科學(xué)技術(shù)的發(fā)展，新能源系統(tǒng)中各種自動裝置對系統(tǒng)頻率的測量精度要求日趨嚴(yán)格[5]。當(dāng)前，離散傅里葉變換（DFT）、快速傅里葉變換（FFT）以及過零點測頻算法是嵌入式裝置中普遍使用的軟件測頻算法[6]。其中DFT 測頻算法當(dāng)采樣點數(shù)比較少時，會發(fā)生頻譜泄漏和柵欄效應(yīng)[7]，并且高次諧波對其影響很大，對以上的問題后期進(jìn)行了各種各樣的改進(jìn)算法，其中加窗函數(shù)修正的DFT 測頻算法[8]存在主要的問題是它的實時性達(dá)不到預(yù)計的要求；文獻(xiàn)[9]結(jié)合插值和迭代的優(yōu)點得到的DFT測頻算法計算復(fù)雜，反應(yīng)周期長，不利于應(yīng)用到嵌入式裝置中。

基于頻率的變化率，本文提出了時變頻率測量算法，解決了頻率隨時間變化時，現(xiàn)有的測頻方法精度不高、誤差大、計算復(fù)雜等問題。首先把原始的信號送入濾波器進(jìn)行濾波，濾除信號中所含的高次諧波、直流分量和噪聲；然后在每一段小區(qū)間內(nèi)對濾波后的信號采用過零點檢測算法[10]進(jìn)行頻率值的計算，將所得到的頻率值作為最小二乘非線性曲線擬合的迭代初值，在每一段小區(qū)間內(nèi)進(jìn)行正弦曲線擬合；最后對擬合后的曲線以調(diào)整采樣頻率的方式對其進(jìn)行二次采樣，將所得到的二次采樣數(shù)據(jù)進(jìn)行動態(tài)頻率測量算法計算頻率，解決現(xiàn)有的算法計算量大、測量精度不高和占用內(nèi)存大等問題。

1 頻率算法的基本原理

1.1 濾波器的設(shè)計

在實際運行的電力系統(tǒng)中，信號波形往往不是標(biāo)準(zhǔn)的正弦曲線，其中夾雜著高次諧波、直流分量以及各種噪聲的復(fù)雜曲線[11]；系統(tǒng)正常運行時，頻率在50 Hz 左右波動，其中的諧波頻率也在隨著工頻不斷變化[12]，這就給頻率的準(zhǔn)確測量帶來了很大的測量誤差，因此需要設(shè)計濾波器。

本文采用基于巴特沃斯型帶通濾波器的設(shè)計[13]來濾除原始信號中所含的噪聲和諧波，文中主要設(shè)計參數(shù)：通帶范圍為45～55 Hz，通帶上限為60 Hz，通帶下限為40 Hz，通帶最大衰減為2 dB，阻帶最小衰減為15 dB。帶通濾波器頻率響應(yīng)如圖1所示。

圖1 帶通濾波器頻率響應(yīng)Fig.1 Frequency response of bandwidth filter

分析圖1可得，當(dāng)信號經(jīng)過這款濾波器后，信號中的高次諧波分量和直流分量的衰減接近于0，最后只剩下50 Hz 附近的頻率可以通過濾波器。從而濾除掉原始信號中摻雜的直流分量和高次諧波分量，減小其對基波信號的干擾。

信號中不僅含有以上對信號的干擾項，而且含有各種隨機噪聲，帶通濾波器不可能把所有的干擾項全部濾除，經(jīng)過濾波器濾波后的信號中仍然含有一部分噪聲，這些未被濾除的噪聲很容易產(chǎn)生頻率泄漏，所以還要對信號的噪聲進(jìn)行處理。

1.2 基于列文伯格-馬夸爾特的最小二乘非線性曲線擬合原理

濾波后的信號中仍然含有部分的噪聲信號，這些噪聲信號很容易造成頻率的泄漏[14]，這就在很大程度上影響了現(xiàn)有測頻算法的測量精度。

本文運用基于列文伯格-馬夸爾特的最小二乘非線性曲線擬合[15]方法對濾波后的信號進(jìn)行正弦曲線擬合。其基本原理為在濾波后信號的每段區(qū)間Vi中把濾波后的信號擬合成標(biāo)準(zhǔn)的正弦曲線Yi，表示為

式中：Yi為第i 段的標(biāo)準(zhǔn)正弦擬合函數(shù)。

標(biāo)準(zhǔn)正弦擬合函數(shù)的各種參數(shù)，如相角、幅值和頻率的最優(yōu)值采用列文伯格-馬夸爾特方法來獲得，其中的初始值設(shè)為0.01。為了減小迭代次數(shù)，把迭代頻率的初始值設(shè)為過零檢測算法所得的頻率值，取0～π 中的任意值作為相角θ 的值。假如得到的參數(shù)值不在合理范圍值之內(nèi)，重新調(diào)整θ 的值對參數(shù)進(jìn)行優(yōu)化估計，直至參數(shù)值在合理范圍值之內(nèi)。最后對擬合后的正弦曲線進(jìn)行二次采樣，進(jìn)一步去除信號中的各種噪聲，從而提高頻率測量的算法精度。

1.3 頻率變化率與周期的關(guān)系

首先假設(shè)被測信號的幅值是Um，頻率為f，相角是φ（f，t），那么信號的表達(dá)式可以表示為

如果f 是隨時間變化的，那么f 的一階導(dǎo)數(shù)肯定不為零，則每一時刻的相角φ（f，t）可以表示為

式中：ω=2π f 是角頻率；dω/dt 是角加速度；φ0是零時刻（t=0）的初始相角；ω0=2π f0是零時刻的初始角速度，f0為初始頻率。以緊鄰3 個過零點處的相位角為研究對象，如圖2所示。

圖2 信號波形與過零點示意圖Fig.2 Signal waveform and zero crossing diagram

把過零點0 設(shè)置為起始點，則在t=0 處，相角φ0=0，頻率為f0。在過零點1 處，相角為2π，過零點2 處，相角為4π。

由式（3）得：

解得：

其中：

由式（5）和式（6）可以求出頻率的變化率d f/dt和過零點0 處的初始頻率f0。

根據(jù)以上所求得的結(jié)果可以計算出每一時刻的頻率。如圖2所示，假設(shè)過零點2 為最新的過零點，采樣點m 為當(dāng)前時刻的采樣點，過零點2 到當(dāng)前時刻之間的采樣值個數(shù)為m，因此當(dāng)前時刻的頻率為

式中：Δt 是過零點2 與采樣點1 之間的時間；Ts為相鄰兩個采樣點之間的時間，即采樣時間間隔。

1.4 計算周期T

1.4.1 Δt 的計算

如圖2，在采樣的時候信號的過零點經(jīng)常并不一定與采樣點重合，要準(zhǔn)確求出周期T，則必須首先要求出過零點與過零點后的第一個采樣點之間的時間Δt。

對于信號的真正過零點是沒有辦法直接求解出來的，但是可以通過其他的方法近似求解，比如，可以把零點附近的波形看作是直線，然后采用相似三角行的方法來求解[16]。

為了解決相似三角形方法解得的過零點值精度不高的問題，本文采用牛頓多項式[17]f（x）來逼近：

接下來就是計算式（8）的各項系數(shù)即a0，a1，a2，a3。

如圖3所示，依據(jù)式（8），從眾多的采樣點中取4 個采樣點，分別為過零點之后的第一個采樣點（tk，uk），過零點之前的三個采樣點（tk-1，uk-1）、（tk-2，uk-2）、（tk-3，uk-3）。根據(jù)這4 個采樣值可以建立向后差分表[18]，見表1。

圖3 牛頓插值多項式示意圖Fig.3 Schematic diagram of Newton interpolation polynomial

表1 牛頓向后差分表Tab.1 Newton backward difference table

其中，

則，式（8）又可表示為

x 與Δt 之間的關(guān)系式為

從而求得式（8）的各項系數(shù)如下：

令f（x）=0，求解牛頓三次插值多項式過零點：

利用牛頓迭代方法求取f（x）=0 的根：

為了減少迭代的次數(shù)，首先采用相似三角形的方法求出迭代初值x0，見圖4所示，由相似三角形原理可得：

圖4 計算過零點時間方法示意圖Fig.4 Schematic diagram of method for calculating zero crossing time

把Δt=-xTs代入式（16），可得：

所以迭代初值可以取為

當(dāng)?shù)冗_(dá)到所要求的精度時可以得到x 的值，最后計算出Δt=-xTs的值。

1.4.2 周期T 的計算

假設(shè)兩個緊鄰過零點之間的采樣點數(shù)為N，那么周期T 的計算公式為

式中：Δt1和Δt2是過零點與過零點之后第一個采樣點之間的時間值。

2 算法的實現(xiàn)過程

（1）對過零點的判斷。當(dāng)uk*uk-1<0 時在uk與uk-1之間必有一個零點。

（2）按式（9）、式（12）計算擬合函數(shù)f（x）的系數(shù)a0～a3。

（3）求解過零點。首先依據(jù)式（18）算出迭代初值，再根據(jù)式（8）、式（13）、式（14）計算出過零點的近似解x 的值。

（4）當(dāng)前周期T 的計算。利用式（11）計算出Δt的值，再根據(jù)式（19）得到當(dāng)前周期T。

（5）最后計算采樣點處的頻率。采用式（7）計算頻率f。

3 算法的仿真

接下來分別對標(biāo)準(zhǔn)正弦波信號、含高次諧波和噪聲、頻率波動情況下，對本文提出的算法進(jìn)行MATLAB 仿真，將仿真所得的結(jié)果與文獻(xiàn)[19]中的DFT 算法和文獻(xiàn)[20]中的過零點算法進(jìn)行比較，以驗證本算法的優(yōu)越性。

3.1 標(biāo)準(zhǔn)正弦波

設(shè)標(biāo)準(zhǔn)正弦信號的數(shù)學(xué)模型為u（t）＝Umsin（2π ft+φ），頻率f 取45～54 之間，φ 角隨機選取。測量數(shù)據(jù)如表2所示。

表2 標(biāo)準(zhǔn)正弦信號下的頻率誤差Tab.2 Frequency error under standard sinusoidal signal

3.2 含有噪聲和諧波

設(shè)含有噪聲和諧波信號的數(shù)學(xué)模型為

頻率f 取45～54 之間，φ 角隨機選取，ξ 為50 dB的白噪聲。測量數(shù)據(jù)如表3所示。

表3 含有噪聲和諧波信號下的頻率誤差Tab.3 Frequency error with noise and harmonic signal

3.3 頻率發(fā)生波動

對發(fā)電機啟動過程的頻率進(jìn)行模擬，取開始的頻率為35 Hz，初始相角為-10°，頻率變化率10，那么頻率和輸入電壓信號模型為

式中：ω0=2π f0；ω=2π f。測量數(shù)據(jù)如表4所示。

表4 頻率發(fā)生波動時的頻率誤差Tab.4 Frequency error in case of frequency fluctuation

4 結(jié)語

在頻率發(fā)生波動的條件下，以緊鄰兩個過零點之間的相位為依據(jù)，推算出了每一時刻頻率的數(shù)學(xué)表達(dá)式，從而提出了一種動態(tài)頻率測量算法。本算法原理簡單，實驗仿真結(jié)果表明，本文算法在標(biāo)準(zhǔn)正弦信號下的頻率測量精度達(dá)到10-6Hz 級，在含有噪聲和諧波信號下的頻率測量精度可達(dá)10-3Hz 級，在頻率波動下的頻率測量精度達(dá)到10-4Hz 級。

本文算法測量頻率的精度高、測量范圍廣、計算速度快、能有效地抵抗噪聲的干擾、在算法的準(zhǔn)確性和穩(wěn)定性方面有所提高，驗證了本算法的可行性和有效性，解決了頻率波動條件下，頻率測量精度不高、計算量大、占用嵌入式內(nèi)存大等問題。