黃瑞，肖宇，劉謀海，單鉉昇，溫和

（1.湖南大學(xué) 電氣與信息工程學(xué)院，湖南 長沙 410082；2.國網(wǎng)湖南省電力有限公司，湖南 長沙 410004；3.湖南大學(xué)智能電氣量測(cè)與應(yīng)用技術(shù)湖南省重點(diǎn)實(shí)驗(yàn)室，湖南長沙 410082）

電網(wǎng)參數(shù)的準(zhǔn)確估計(jì)是電能計(jì)量、電能質(zhì)量檢測(cè)、繼電保護(hù)的重要基礎(chǔ).電力系統(tǒng)的干擾，例如故障，大負(fù)載的開/關(guān)，發(fā)電機(jī)的斷開，線路切換等，會(huì)導(dǎo)致發(fā)電機(jī)之間的轉(zhuǎn)子角發(fā)生振蕩從而導(dǎo)致電網(wǎng)產(chǎn)生低頻振蕩.如何在低頻振蕩情況下實(shí)現(xiàn)電網(wǎng)參數(shù)準(zhǔn)確快速測(cè)量，對(duì)維持電力系統(tǒng)穩(wěn)定運(yùn)行起著至關(guān)重要的作用[1].

近年來，很多學(xué)者與研究機(jī)構(gòu)對(duì)于電網(wǎng)動(dòng)態(tài)參數(shù)測(cè)量進(jìn)行了研究.其中常用的測(cè)量方法主要有：基于正弦信號(hào)模型的檢測(cè)算法；周期法及其改進(jìn)算法，主要包括過零檢測(cè)法、水平交點(diǎn)法、高次修正函數(shù)法和最小二乘多項(xiàng)式曲線擬合法等；隨機(jī)模型算法，主要包括最小二乘法[2]、最小絕對(duì)值近似法[3]、牛頓迭代算法和線性濾波算法等[4-5].基于周期信號(hào)模型的檢測(cè)算法應(yīng)用最為廣泛，主要包括離散傅里葉變換（Discrete Fourier transform，DFT）[6-7]和快速傅里葉變換類算法及其改進(jìn)算法.

理想情況下的電網(wǎng)信號(hào)為標(biāo)準(zhǔn)正弦信號(hào)，以正弦信號(hào)模型為基礎(chǔ)的電網(wǎng)頻率測(cè)量方法得到大量應(yīng)用，其中應(yīng)用最為廣泛的方法就是基于離散傅里葉變換的參數(shù)估計(jì)方法.在同步采樣情況下，僅用一個(gè)周波的采樣信息就可實(shí)現(xiàn)對(duì)參數(shù)的準(zhǔn)確估計(jì).但是在非同步采樣情況下，其精度受到頻譜泄漏和柵欄效應(yīng)的影響，需要采用加窗插值[8-9]等方法進(jìn)行優(yōu)化.

當(dāng)電網(wǎng)信號(hào)處于動(dòng)態(tài)變化時(shí)，由于DFT 方法的局限性，相關(guān)動(dòng)態(tài)參數(shù)無法準(zhǔn)確估計(jì).低頻振蕩的場(chǎng)景下，電網(wǎng)中的電壓與電流信號(hào)不再是穩(wěn)態(tài)信號(hào)，其幅值會(huì)產(chǎn)生周期性的波動(dòng)，如不及時(shí)監(jiān)測(cè)與控制，將進(jìn)一步導(dǎo)致系統(tǒng)間產(chǎn)生功率震蕩，對(duì)傳輸線路和用電設(shè)備造成不良影響，破壞系統(tǒng)的安全穩(wěn)定性.為了實(shí)現(xiàn)對(duì)低頻振蕩信號(hào)及時(shí)準(zhǔn)確的測(cè)量，國內(nèi)外學(xué)者對(duì)相關(guān)算法進(jìn)行了研究，文獻(xiàn)[10]將Prony 方法與數(shù)字濾波相結(jié)合，相較于傳統(tǒng)Prony 方法提高了抗噪性.卡爾曼濾波算法對(duì)于動(dòng)態(tài)信號(hào)有著較好的跟蹤效果，文獻(xiàn)[11]提出了基于卡爾曼濾波的測(cè)量方法，有較好抗噪性和測(cè)量精度.迭代濾波算法[12-13]在頻率接近工頻的情況下能實(shí)現(xiàn)高精度的測(cè)量，但是在頻偏較大與幅值動(dòng)態(tài)變化的情況下誤差較大.小波變換[14-15]雖然適用于動(dòng)態(tài)信號(hào)分析，但是計(jì)算量大，實(shí)時(shí)性不足.文獻(xiàn)[16]提出了適用于含阻尼振蕩信號(hào)的插值方法，但該方法僅適用于事后離線分析，不適用于低頻振蕩在線監(jiān)測(cè).文獻(xiàn)[17]提出了基于同步向量測(cè)量數(shù)據(jù)的振蕩參數(shù)測(cè)量方法，但是每次運(yùn)算需要累計(jì)約2 s 的同步相量數(shù)據(jù).文獻(xiàn)[18]提出的泰勒傅里葉方法通過泰勒級(jí)數(shù)展開簡(jiǎn)化計(jì)算，在穩(wěn)態(tài)情況下精度高，諧波抑制能力強(qiáng)，在動(dòng)態(tài)情況下響應(yīng)速度快.文獻(xiàn)[19]通過將Prony 方法與泰勒傅里葉級(jí)數(shù)結(jié)合，提出了泰勒Prony 算法，實(shí)現(xiàn)了在低頻振蕩場(chǎng)景下的參數(shù)估計(jì).

以上提到的方法，均為利用單相信號(hào)實(shí)現(xiàn)參數(shù)估計(jì)，而在電力系統(tǒng)中通常要將三相作為一個(gè)整體進(jìn)行考量，考慮到三相系統(tǒng)的平衡特性，本文提出了一種基于復(fù)頻譜插值DFT 的參數(shù)估計(jì)方法.通過等幅值克拉克變換，將三相系統(tǒng)中的參數(shù)估計(jì)問題轉(zhuǎn)換為復(fù)指數(shù)信號(hào)的參數(shù)估計(jì)，本文通過引入指數(shù)衰減系數(shù)，構(gòu)建三相信號(hào)模型，并基于所構(gòu)建模型對(duì)各項(xiàng)估計(jì)參數(shù)進(jìn)行了推導(dǎo)，給出了各項(xiàng)參數(shù)的表達(dá)式.在穩(wěn)態(tài)信號(hào)和動(dòng)態(tài)信號(hào)的仿真驗(yàn)證中，通過與其他方法的對(duì)比，本文方法的準(zhǔn)確性得到驗(yàn)證.

1 信號(hào)模型

電力系統(tǒng)故障、線路切換、發(fā)電機(jī)斷開與連接以及斷開或接入大量負(fù)載會(huì)導(dǎo)致電力系統(tǒng)的波動(dòng)，進(jìn)而引發(fā)電壓或電流低頻振蕩.通常情況下，在低頻振蕩或?qū)ΨQ故障時(shí)，電力系統(tǒng)仍然保持對(duì)稱特性[20].

在文獻(xiàn)[21]和[22]中提出，指數(shù)衰減的正弦信號(hào)可以用來近似擬合低頻振蕩和對(duì)稱故障，本文將阻尼系數(shù)引入平衡的三相信號(hào)模型中，使得信號(hào)模型更加符合低頻振蕩下電網(wǎng)信號(hào)的特點(diǎn).考慮到低頻振蕩的對(duì)稱特性，將振蕩信號(hào)建模為一組對(duì)稱的指數(shù)衰減正弦信號(hào)，可以表示為：

式中：A，τ，φ 和ω 分別代表幅值，衰減系數(shù)，相位和角頻率.信號(hào)的采樣率為fs，Ts=1/fs.

2 基于復(fù)頻譜插值DFT 的參數(shù)估計(jì)方法

對(duì)于三相三線系統(tǒng)，其三相電流及以任意點(diǎn)為電壓參考點(diǎn)的三相電壓經(jīng)克拉克變換后得到的α、β兩相瞬時(shí)電壓uα、uβ和瞬時(shí)電流iα、iβ可完全表征此三相三線系統(tǒng).因此，基于前文給出的三相信號(hào)模型，做等幅值克拉克變換可得：

在三相系統(tǒng)中，經(jīng)等幅值克拉克變換得到的yα和yβ可以通過如下方式重構(gòu)信號(hào)，可得到含有正序分量和負(fù)序分量的復(fù)信號(hào)ycomples：

式中，下標(biāo)+表示正序分量的參數(shù)，下標(biāo)-表示負(fù)序分量的參數(shù).當(dāng)三項(xiàng)系統(tǒng)為對(duì)稱系統(tǒng)時(shí)，負(fù)序分量為0，公式（3）可簡(jiǎn)化為：

由此，可以得到一個(gè)新的復(fù)信號(hào)序列，且根據(jù)克拉克變換的性質(zhì)，復(fù)信號(hào)序列與三相系統(tǒng)有著相同的頻率、幅值和相位.因此，三相系統(tǒng)的參數(shù)估計(jì)問題轉(zhuǎn)化成了對(duì)復(fù)信號(hào)ycomplex的參數(shù)估計(jì).接下來將對(duì)A，τ，φ 和ω 的求解進(jìn)行推導(dǎo).

使用長度為N 的矩形窗對(duì)復(fù)信號(hào)ycomplex（n）加權(quán)，并對(duì)序列做N 點(diǎn)離散傅里葉變換可得：

由公式（6）可知，λ 和ρ 通過這兩個(gè)未知量包含了幅值、相位和頻率的參數(shù)，因此，通過求解λ 和ρ便能得出所有的目標(biāo)參數(shù).將公式（6）改寫為矩陣形式可得：

DFT 序列Y（k）可由未知參數(shù)λ 和ρ 的線性組合表示，在N>2 的情況下，有兩個(gè)未知參數(shù)和N 個(gè)可用方程，那么未知參數(shù)λ 和ρ 的值可通過解線性方程的方法求得，選擇幅值最大和次大的兩根譜線，其序號(hào)分別記為k1和k2，可構(gòu)建線性方程組：

由此求解可得到未知參數(shù)λ 和ρ，根據(jù)λ 和ρ與幅值頻率相位的關(guān)系可得：

在三相不平衡的情況下，由于存在負(fù)序分量，需要考慮負(fù)頻成分的干擾，為實(shí)現(xiàn)對(duì)正序分量參數(shù)的準(zhǔn)確估計(jì)，可將公式（6）擴(kuò)展為如公式（10）所示：

需要求解的未知數(shù)增加為4 個(gè)，因此需通過在代表正序分量和負(fù)序分量的兩處峰值附近選擇最大和次大共計(jì)四根譜線構(gòu)建矩陣如式（11）所示，其中等式左側(cè)矩陣代表未知量矩陣.

根據(jù)式（11）以及λ+、λ-和ρ+、ρ-的關(guān)系，可分別求得λ+、λ-、ρ+、ρ-.將λ+、ρ+代入公式（9），即可得到正序分量的相關(guān)參數(shù).

通常情況下，幅值最大和次大兩根譜線的選擇十分簡(jiǎn)單，但是在信號(hào)頻率十分接近工頻時(shí)，峰值譜線左右兩側(cè)的譜線幅值很小，次大幅值譜線的選擇會(huì)因?yàn)樵肼暤挠绊憣?dǎo)致誤判，從而影響估計(jì)精度.

將信號(hào)的實(shí)際角頻率記為ω，最大幅值譜線位置記為k1，則

由公式（6）和（10）可得

則任意譜線與Y（k1）存在如下關(guān)系式

記k2=k1+1，k3=k1-1，在δ 趨近于0 時(shí)，由公式（15）可得：

由公式（16）（17）可知，在δ 趨近于0 時(shí)，最大幅值譜線兩側(cè)的譜線幅值十分接近，在噪聲的影響下，直接地比較幅值大小容易造成次大譜線選擇的誤判.為了降低噪聲對(duì)譜線篩選造成的影響，本文采取文獻(xiàn)[23]中的方式進(jìn)行次大譜線的選擇：

其中Re 表示取實(shí)部，

當(dāng)Δ--Δ+＞0 時(shí)，k2=k1+1，Δ--Δ+＜0 時(shí)，k2=k1-1.

基于復(fù)頻譜插值DFT 的參數(shù)估計(jì)方法的流程可總結(jié)為如下幾個(gè)步驟：

1）設(shè)定采樣參數(shù)，同步采集三相電壓信號(hào)ya、yb、yc.

2）通過克拉克變換，將采集的三相電壓信號(hào)序列轉(zhuǎn)換為復(fù)信號(hào)序列ycomplex.

3）以長度為N 的矩形窗截取復(fù)信號(hào)序列ycomplex并做N 點(diǎn)DFT 得到DFT 序列Y（k）.

4）在三相平衡情況下，篩選最大和次大譜線，序號(hào)記為k1和k2，并根據(jù)公式（8）求解未知參數(shù)λ 和ρ.在三相不平衡的情況下，在代表正序分量和負(fù)序分量的兩處峰值附近選擇最大和次大共計(jì)四根譜線，并根據(jù)公式（11）求解未知參數(shù)λ+/-和ρ+/-.

5）根據(jù)得到的λ 和ρ 以及公式（9）求解幅值、頻率和相位參數(shù).

3 仿真分析

3.1 白噪聲影響下的仿真分析

為檢驗(yàn)本文所提出方法的性能，本節(jié)將在穩(wěn)態(tài)信號(hào)與動(dòng)態(tài)信號(hào)兩種場(chǎng)景下下進(jìn)行仿真驗(yàn)證.在穩(wěn)態(tài)信號(hào)的場(chǎng)景下，通過在不同信噪比（SNR）的條件下，估計(jì)穩(wěn)態(tài)信號(hào)的參數(shù)，并將本文方法與其他方法的估計(jì)精度進(jìn)行對(duì)比.在動(dòng)態(tài)信號(hào)場(chǎng)景下，分別在幅值階躍、相位階躍以及頻率斜升三種情況下對(duì)算法進(jìn)行測(cè)試，對(duì)比本文方法與其他方法的動(dòng)態(tài)跟蹤性能.為模擬三相不平衡的情況，所有仿真試驗(yàn)中均在模擬三相信號(hào)中加入負(fù)序分量，幅值為正序分量的30%.參與對(duì)比的方法有兩類，一類是加窗插值法，加窗類型分別為最大旁瓣衰減窗和漢寧窗[24]；另一類為基于Prony 的參數(shù)估計(jì)方法，分別為傳統(tǒng)Prony方法以及泰勒Prony 方法[19]。

穩(wěn)態(tài)情況下仿真參數(shù)設(shè)置為：三相信號(hào)的基頻f=49.5 Hz，采樣率fs=6 kHz，幅值A(chǔ)=1，初相位φ=0.1 rad.采樣窗長設(shè)置為N=256，SNR 變化范圍為30 dB 到90 dB，相位頻率和幅值估計(jì)的均方根誤差如圖1 所示.

圖1 相位、頻率與幅值估計(jì)誤差隨信噪比變化圖Fig.1 RMSE of phase，frequency and amplitude estimation results versus SNR

在白噪聲的影響下，隨著信噪比的增加，所有方法的相位，頻率和幅值估計(jì)結(jié)果的均方根誤差均呈現(xiàn)下降趨勢(shì).不論是在信噪比較低的情況下，還是在高信噪比的條件下，本文方法的估計(jì)精度均高于其余方法.基于Prony 的算法在信噪比較低的情況下估計(jì)精度較低，但是隨著信噪比的提高，估計(jì)精度有所提高，并在信噪比高于70 dB 的情況下優(yōu)于加窗插值算法.

動(dòng)態(tài)情況下，幅值階躍設(shè)置為100%，相位階躍設(shè)置為π/4，頻率斜升速率設(shè)置為1 Hz/s.初始頻率均設(shè)置為f=50 Hz.采樣率fs=6 kHz，采樣窗長設(shè)置為N=256，SNR 設(shè)置為50 dB.

幅值階躍情況下的幅值估計(jì)結(jié)果如圖2 所示.所有方法均在1 個(gè)周波左右的時(shí)間后實(shí)現(xiàn)了對(duì)幅值的準(zhǔn)確跟蹤.其中，加窗插值方法延遲略小于一個(gè)周波的時(shí)間，且波動(dòng)最小.本文方法在階躍瞬間估計(jì)結(jié)果出現(xiàn)超調(diào)，約為階躍值的30%.基于Prony 的估計(jì)方法在階躍瞬間的估計(jì)結(jié)果出現(xiàn)較大幅度的波動(dòng).

圖2 幅值估計(jì)結(jié)果圖Fig.2 Estimation result of amplitude and phase

相位階躍情況下的相位估計(jì)結(jié)果如圖3 所示.所有方法均在1 個(gè)周波左右的時(shí)間后實(shí)現(xiàn)了對(duì)相位的準(zhǔn)確跟蹤.本文方法在階躍瞬間超調(diào)最小，其余方法的結(jié)果均出現(xiàn)了較大的波動(dòng).

圖3 相位估計(jì)結(jié)果圖Fig.3 Estimation result of amplitude and phase

在頻率斜升情況下的頻率估計(jì)結(jié)果如圖4 所示，本文方法的頻率估計(jì)結(jié)果波動(dòng)最小，全程實(shí)現(xiàn)對(duì)頻率的準(zhǔn)確跟蹤，其他方法的估計(jì)結(jié)果均存在波動(dòng)，其中基于Prony 方法的結(jié)果波動(dòng)最大.

圖4 頻率斜升估計(jì)結(jié)果圖Fig.4 Estimation result of frequecy ramp

3.2 低頻振蕩情況下的仿真分析

在低頻振蕩的情況下，電網(wǎng)信號(hào)的各項(xiàng)參數(shù)都處于動(dòng)態(tài)變化，本節(jié)將通過幅度相位調(diào)制信號(hào)模擬典型的低頻振蕩信號(hào)，對(duì)本文所提出方法的性能進(jìn)行評(píng)估，參與對(duì)比的方法為Prony 與泰勒Prony 的參數(shù)估計(jì)方法.

基于公式（1）的穩(wěn)態(tài)三相系統(tǒng)模型，對(duì)幅度與相位進(jìn)行調(diào)制，調(diào)制后的幅值與相位的表達(dá)式如下所示：

其中A0=1，φ0=0.5，A1=0.1，φ1=0.05，fa=fφ=5，τ1=0.5，τ2=0.4.

基頻頻率f=50 Hz，采樣頻率fs=5 kHz，信號(hào)的信噪比SNR 設(shè)置為60 dB，考慮到在低頻振蕩的情況下測(cè)量結(jié)果的實(shí)時(shí)性，參數(shù)估計(jì)的延遲最好不要超過一個(gè)周波，因此采樣窗長設(shè)置為N=128.仿真過程采取滑窗的形式，每次滑動(dòng)1 個(gè)點(diǎn)，估計(jì)各個(gè)時(shí)刻幅值和相位參數(shù)，算法估計(jì)精度通過綜合矢量誤差（TVE）來評(píng)估.

綜合矢量誤差（TVE）結(jié)合了幅值誤差和相角誤差，是對(duì)相量誤差的綜合衡量，TVE 能夠更全面反映出動(dòng)態(tài)相量估計(jì)的準(zhǔn)確度，本文方法與TFT 方法的TVE 如圖5 所示，圖中反映了在50 個(gè)工頻周期數(shù)內(nèi)的最大綜合向量誤差，可以看到本文的TVE 均保持在0.5%以下，Prony 方法的TVE 最大可達(dá)到1.5%，泰勒Prony 方法的TVE 介于本文方法與Prony 方法之間.以此可見本文提出的方法在低頻振蕩的場(chǎng)景下有更高、更穩(wěn)定的估計(jì)精度.

圖5 綜合矢量誤差Fig.5 Total vector error

3.3 三機(jī)九節(jié)點(diǎn)系統(tǒng)中的仿真分析

為了驗(yàn)證該方法在低頻振蕩條件下的參數(shù)估計(jì)性能，本文在IEEE 標(biāo)準(zhǔn)的三機(jī)九節(jié)點(diǎn)系統(tǒng)進(jìn)行仿真模擬，三機(jī)九節(jié)點(diǎn)系統(tǒng)如圖6 所示.三機(jī)九節(jié)點(diǎn)系統(tǒng)采用PSCAD 軟件進(jìn)行仿真.仿真條件設(shè)置如下：繼電保護(hù)裝置R 位于7 號(hào)和5 號(hào)總線之間的線路中.為了在系統(tǒng)中產(chǎn)生低頻振蕩，在總線7 和總線8 之間的線路中引入三相接地故障，故障在t=1 s 時(shí)開始，并在0.2 s 后通過斷開位于該線路兩端的斷路器來清除故障，此時(shí)該系統(tǒng)產(chǎn)生低頻振蕩，繼電保護(hù)裝置R 處可以觀測(cè)到低頻振蕩信號(hào).

圖6 三機(jī)九節(jié)點(diǎn)系統(tǒng)Fig.6 IEEE 9 bus system

繼電保護(hù)裝置的采樣率設(shè)置為5 kHz，采樣窗長設(shè)置為N ＝128，約為0.8 個(gè)工頻周期.本節(jié)將對(duì)繼電保護(hù)裝置R 處所采集到的三相電壓信號(hào)進(jìn)行分析，所有數(shù)據(jù)在MATLAB 中處理，信噪比SNR=40 dB，其中a 相電壓信號(hào)如圖7 所示.

圖7 a 相電壓Fig.7 Voltage signal of phase a

幅值估計(jì)結(jié)果如圖8 所示，通過對(duì)比可以看出，本文方法在故障發(fā)生前、故障發(fā)生時(shí)、以及故障清除后的三種情況下都能夠準(zhǔn)確地實(shí)現(xiàn)幅值的估計(jì)，而且在故障發(fā)生時(shí)能夠?qū)﹄妷后E降實(shí)現(xiàn)快速準(zhǔn)確的跟蹤，在故障發(fā)生后約0.02 s 后，本文方法實(shí)現(xiàn)了對(duì)幅值的正確跟蹤，期間與真實(shí)值的偏差不超過10%，相較于對(duì)比方法，本文方法在故障發(fā)生階段的估計(jì)結(jié)果波動(dòng)最小.而Prony 方法與泰勒Prony 方法更易受噪聲影響，在故障發(fā)生后，估計(jì)結(jié)果存在較大的波動(dòng)，偏離真實(shí)值最大超過40%以上，無法準(zhǔn)確跟蹤故障期間電壓的變化.

圖8 幅值估計(jì)結(jié)果比較（b 為a 中A 區(qū)域的放大圖）Fig.8 .Comparison of amplitude measurement results（area A is magnified in b）

4 結(jié)論

本文提出了一種基于復(fù)頻譜插值DFT 的動(dòng)態(tài)振蕩信號(hào)測(cè)量方法.通過等幅值克拉克變換，將三相系統(tǒng)中的參數(shù)估計(jì)問題轉(zhuǎn)換為復(fù)指數(shù)信號(hào)的參數(shù)估計(jì)，通過引入指數(shù)衰減系數(shù)，構(gòu)建三相信號(hào)模型，并基于所構(gòu)建模型對(duì)各項(xiàng)估計(jì)參數(shù)進(jìn)行了推導(dǎo)，給出了各項(xiàng)參數(shù)的表達(dá)式.通過理論仿真與電力系統(tǒng)模型仿真驗(yàn)證，結(jié)果均表明本算法具有良好的性能，主要體現(xiàn)在以下方面：

1）本文所提出的方法可在采樣點(diǎn)數(shù)較少的情況下實(shí)現(xiàn)對(duì)頻率幅值的準(zhǔn)確測(cè)量，且能在信噪比較低的場(chǎng)景下提供較為準(zhǔn)確的估計(jì)結(jié)果.

2）在處理動(dòng)態(tài)信號(hào)時(shí)，仿真實(shí)驗(yàn)的結(jié)果表明，本文所提出的方法在抗噪性能和動(dòng)態(tài)參數(shù)估計(jì)性能上都具有一定優(yōu)勢(shì)，可以在噪聲干擾的情況下對(duì)低頻振蕩狀態(tài)下的三相信號(hào)參數(shù)進(jìn)行較為準(zhǔn)確的估計(jì).

3）本文方法最多只需要四根DFT 譜線的參數(shù)即構(gòu)建矩陣實(shí)現(xiàn)對(duì)參數(shù)的估計(jì)，運(yùn)算量主要為DFT 運(yùn)算，而且只需采集一個(gè)周波左右長度的數(shù)據(jù)即可，現(xiàn)有DSP 的運(yùn)算能力足以滿足本文方法的需求，因此本文方法可移植至嵌入式平臺(tái)進(jìn)行應(yīng)用.