曾林俊，肖輝，江維，鄧仕燊，楊俊琛

（1.長沙理工大學(xué)，長沙410114；2.國網(wǎng)株洲供電公司，湖南 株洲412000；3.國網(wǎng)福建龍海市供電公司，福建 龍海363100）

0 引 言

隨著我國電力工業(yè)的不斷發(fā)展，電力系統(tǒng)規(guī)模不斷擴大，已經(jīng)擴展成了大區(qū)域互聯(lián)電網(wǎng)。此種情況有利于改善輸配電的運行經(jīng)濟性和可靠性，但區(qū)域電網(wǎng)互聯(lián)有可能引發(fā)低頻振蕩現(xiàn)象［1］。在采用勵磁控制器的條件下，低頻振蕩發(fā)生的幾率不斷提高，如果無法很好的對低頻振蕩進行抑制，電力系統(tǒng)的穩(wěn)定性將遭到破壞。而對低頻振蕩的模式識別［2－3］，有助于了解系統(tǒng)的動態(tài)性能，為安全預(yù)警提供較好的分析結(jié)果。

現(xiàn)有的低頻振蕩模式分析方法，基于數(shù)值解的分析法如傳遞函數(shù)法，及基于實測信號的分析法，如ESPRIT、FFT、小波分析等大部分對線性系統(tǒng)有良好的辨識效果，傳遞函數(shù)法通過對系統(tǒng)輸入信號，再由反饋信號尋找到系統(tǒng)內(nèi)部的某種對應(yīng)關(guān)系，進而可判斷系統(tǒng)的穩(wěn)定狀況等，但當系統(tǒng)非線性程度較強時，傳遞函數(shù)法［4］所計算的變量關(guān)系將會十分復(fù)雜；ESPRIT算法［5-6］假定振蕩信號可表示為一系列指數(shù)函數(shù)的線性組合，未一起使用復(fù)觀測數(shù)據(jù)及其共軛數(shù)據(jù)進行參數(shù)估計，且在低信噪比下的性能較低，使得辨識精度不高，同時對非線性信號的分析有一定的限制；FFT［7-8］在分析含噪的非線性、非平穩(wěn)信號上顯得不夠理想，無時域分析功能，無法自適應(yīng)反映信號的時變特性，且頻率分辨率不高，難以滿足非平穩(wěn)信號分析的需要；小波分析法［9］可用于非線性系統(tǒng)，但小波基難以選擇且存在頻率分辨問題?，F(xiàn)有的分析方法對于非線性非高斯的復(fù)雜系統(tǒng)的去噪和參數(shù)估計方面處理效果不是很理想，有一定的誤差，然而電力系統(tǒng)大部分都是含噪的非高斯非線性復(fù)雜系統(tǒng)，現(xiàn)有方法無法很好地滿足非高斯非線性復(fù)雜系統(tǒng)的需要。

為克服現(xiàn)有方法的不足，本文提出一種識別低頻振蕩模式的新方法，即將前置粒子濾波去噪算法與時頻分析法-EEMD［10-12］相結(jié)合，利用粒子濾波對非高斯非線性系統(tǒng)處理的優(yōu)勢，對振蕩信號進行去噪，再對去噪后的信號使用EEMD進行模態(tài)分解，檢測出低頻振蕩模式。最后，通過算例仿真，驗證了本方法的正確性和可靠性。

1 粒子濾波算法

1.1 粒子濾波理論

粒子濾波［13-16］（Particle Filter，PF），又稱為序貴蒙特卡羅（Sequential Monte Carlo，SMC）方法，是一種基于SMC方法的貝葉斯濾波技術(shù)。其基本原理是尋找一組在狀態(tài)空間傳播的隨機粒子（樣本）描述系統(tǒng)的狀態(tài)，通過SMC法簡化貝葉斯估計中的積分運算，進而得出系統(tǒng)狀態(tài)的最小均方差估計。在粒子數(shù)接近無窮時能逼近服從任意概率分布的系統(tǒng)狀態(tài)。相對其它濾波技術(shù)而言，粒子濾波［17］不需要對系統(tǒng)狀態(tài)作任何先驗性假設(shè)，原理上能用于可用于狀態(tài)空間模型描述的任何系統(tǒng)。

系統(tǒng)的狀態(tài)空間方程包含狀態(tài)方程和觀測方程，其模型描述如下：

狀態(tài)方程：

觀測方程：

式中xk∈Rnx是系統(tǒng)的狀態(tài)值；yk∈Rny是觀測值，f（·）、h（·）表示狀態(tài)方程和觀測方程；uk為過程噪聲；vk為觀測噪聲，其都服從高斯分布。模型中假設(shè)狀態(tài)轉(zhuǎn)移過程服從一階馬爾科夫模型，即時刻k的狀態(tài)xk只和前一時刻k－1的狀態(tài)xk－1有關(guān)，且模型中觀測值間相互獨立，即觀測值yk只與k時刻的狀態(tài)量xk有關(guān)。根據(jù)貝葉斯估計原理，在已知觀測量 y1：k的前提下，x1：k的后驗概率分布為：

其中，p（x0：k）表示先驗概率密度；p（y1：k｜x0：k）表示觀測量，為y1：k時的似然概率密度。

將上式（3）與Monte Carlo原理相結(jié)合，從后驗概率密度 p（x0：k｜y1：k）中抽取 N個相互獨立的樣本來求解函數(shù)估計值，如期望或方差。但是在實際很難直接通過后驗概率密度函數(shù)取樣，一般引入重要性密度函數(shù) q（x0：k｜y1：k）和其相應(yīng)的權(quán)值，權(quán)值計算式為：

為進一步解決計算量較大的問題，將權(quán)值按照下式進行歸一化處理：

從而根據(jù)以上各式，后驗概率密度可以近似的表示為：

其中，δ（）·表示狄拉克delta函數(shù)。

1.2 重要性函數(shù)及重采樣

PF算法存在的主要不足為粒子數(shù)匱乏，樣本權(quán)值退化，粒子數(shù)多樣性缺失?，F(xiàn)在對PF算法處理此不足最有效的方法為選擇重要性函數(shù)和使用重采樣方法［18-20］。

選取重要性函數(shù)的目的是要使重要性權(quán)重的方差最小，在本文中將狀態(tài)變量的轉(zhuǎn)移概率密度函數(shù)作為重要性概率密度函數(shù)，即：q（x0：k｜y1：）k＝p（xk｜xk－1）。

為了彌補權(quán)值退化的不足，Gordon等提出重采樣方法，使得PF算法又重新得到人們的重視，重采樣的思想即減少權(quán)重小的粒子數(shù)，增加權(quán)重較大的粒子數(shù)，此過程是對后驗概率密度函數(shù)再采樣N次，得到新的支撐點集合，重采樣操作得到的粒子獨立同分布，且采樣點被重新賦予的權(quán)值是相等的，然而重采樣方法也往往帶來負面的影響即樣本枯竭，因此我們需要考慮退化程度問題，當PF濾波退化到一定程度時才應(yīng)考慮重采樣。通過有效粒子數(shù)Neff來衡量粒子退化程度，按下式進行選?。?/p>

濾波前預(yù)先設(shè)定一個閾值門限Nthres，有效粒子數(shù)的值必須大于設(shè)定的閾值Nthres，否則重采樣；然后根據(jù)有效粒子數(shù)對應(yīng)的歸一化權(quán)值更新粒子集。

1.3 低頻振蕩信號去噪實現(xiàn)

通過獲取電力系統(tǒng)低頻振蕩時的信號，根據(jù)信號建立狀態(tài)方程和觀測方程，根據(jù)先驗概率密度p（x0：k）得到采樣粒子，對粒子集進行初始化及重要性采樣，通過重要性密度函數(shù)生成樣本，對樣本計算權(quán)值及歸一化權(quán)值，判定是否進行重采樣，迭代結(jié)束后得到粒子的最優(yōu)狀態(tài)估計，從而達到對原始信號去噪的結(jié)果。

在電力系統(tǒng)中，調(diào)度部門關(guān)注的是低頻振蕩的各個模式，為調(diào)度人員提供必要的信息，使得能夠?qū)Φ皖l振蕩進行一定的抑制。在對低頻振蕩模式進行辨識前，通過使用粒子濾波去除信號噪聲，還原出信號的特征，提高了接下來的模式辨識的精度。

2 基于EEMD的模式辨識

對去噪后的信號通過EEMD［21-24］方法進行模式辨識?；舅枷霝橥ㄟ^疊加高斯白噪聲后的多次EMD分解。向原始信號加入頻率均勻的白噪聲背景時，不同時間尺度的信號區(qū)域都能自動映射到和背景白噪聲相關(guān)的合適的尺度上。通過不相關(guān)隨機序列的統(tǒng)計均值為0，對多次加入噪聲后的總體取均值后，噪聲將會被消除，總體均值被認為是信號本身，達到消除模態(tài)混疊問題。算法主要步驟為：

1）給原始信號x（t）中多次加入等長度隨機正態(tài)分布的白噪聲ni（t）（白噪聲標準差為原始信號標準差的0.1～0.4倍）即：

式中xi（）t為第i次加入高斯白噪聲的合成信號。所加噪聲的幅度會直接影響EEMD的分解效果，如果幅度過小，噪聲將失去補充尺度的作用，無法構(gòu)成不同頻率上的連續(xù)性。

2）將xi（）t分別進行EMD分解，得到各個IMF分量記為 Cij（t），與一個余項記為 rin（t），則有：

3）重復(fù)步驟1和步驟2操作N次。將以上各個IMF進行總體取均值運算，得到最終的IMF分量組為：

式中Cj（）t表示EEMD分解后得到最終的第j個IMF分量。N表示加入高斯白噪聲的總體個數(shù)，總體個數(shù)越大，對應(yīng)的白噪聲的IMF和將更接近于0，得到的結(jié)果也更趨近于真實值。

將分解出的各階IMF分量Cj（）t進行Hilbert變換可得：

可以得到瞬時相位 θi（t）、角頻率 ωi（t）、頻率fi（t）如下：

3 基于粒子濾波的EEMD分解

本文所述方法，將粒子濾波作為前置單元，去除非線性、非高斯低頻振蕩信號的噪聲，基于上兩節(jié)的理論，算法流程圖如圖1所示，本文方法的主要步驟如下：

（1）根據(jù)低頻振蕩信號建立粒子濾波的狀態(tài)方程和觀測方程，如式（1）和式（2）；

（2）初始化；k＝0，從先驗概率密度 p（x0：k）采樣得到采樣粒子；

（3）重要性采樣：在采樣粒子中依據(jù)重要性密度函數(shù) q（x0：k｜y1：k）中隨機抽取N個有限的樣本；

（4）對抽取的N個樣本利用式（4）計算相對應(yīng)的權(quán)值wk；

（5）將步驟3得到的權(quán)值根據(jù)式（5）進行歸一化處理，得到歸一化權(quán)值

（6）利用式（7）判定重采樣：有效粒子數(shù)的值必須大于設(shè)定的閾值Nthres，否則重采樣，繼續(xù)返回步驟（2）繼續(xù)執(zhí)行；

（8）對去噪后的信號使用EEMD進行分解，得到IMF；

（9）對 IMF進行 Hilbert變換，根據(jù)式（11）～式（14）計算得出各個IMF的瞬時頻率，從而辨識出低頻振蕩模式。

圖1 算法流程圖Fig.1 Algorithm flow chart

4 算例分析

為測試本方法對低頻振蕩模式辨識的有效性，取以下信號（單位為MW）：

該信號包含三個振蕩模式，頻率分別為0.7 Hz，1.5 Hz，0.5 Hz，該信號包含信噪比SNR為20 dB的高斯白噪聲ω（n），利用粒子濾波作為前置濾波單元，對上述信號進行去噪處理，首先分別建立粒子濾波所需的狀態(tài)方程和觀測方程，將原始信號s（t）作為狀態(tài)方程，而觀測方程則由原始信號加上隨機噪聲得到如式（16）、式（17）所示。

式中的xk與yk分別表示系統(tǒng)的狀態(tài)量和觀測量，并且假設(shè)初始狀態(tài)值為x0＝0.1，初始狀態(tài)分布p（ x0） ～N（0，1），觀測噪聲 rk～N（0，5），時間迭代步長為50，粒子數(shù) N＝1 000，定義有效區(qū)間為［5，95］，為了防止邊界效應(yīng)的影響，在下述操作過程中將區(qū)間擴大為［0，100］，根據(jù)以上參數(shù)進行粒子濾波去噪處理，其仿真結(jié)果如圖2所示。

圖2 粒子濾波去噪結(jié)果Fig.2 Particle filtering de－noising results

根據(jù)圖2可以看出，經(jīng)過粒子濾波去噪后的信號基本可以與原信號很好的吻合，誤差較小，基本保留了原始信號的特征，去噪效果較好。因此，使用粒子濾波對信號進行去噪處理是有效的、精確的。

為進一步辨識出信號的低頻振蕩模式，對去噪后的信號使用改進的EMD方法EEMD來進行分解，分解得到3個IMF1～IMF3結(jié)果如圖3（a）～圖3（c）所示。

圖3 EEMD分解后的各IMF分量Fig.3 IMFs by EEMD decomposition

對得到的IMF1～IMF3進行提取瞬時頻率參數(shù)操作，提取的IMF1～IMF3的瞬時頻率結(jié)果如圖4（a）～圖4（c）所示。

經(jīng)過上述操作后得到的瞬時頻率存在邊界效應(yīng)，為了更好的對低頻振蕩模式進行識別，減少邊界效應(yīng)對最終模式辨識的結(jié)果影響，對得到的瞬時頻率f1～f3在有效區(qū)間［5，95］內(nèi)進行求均值處理，各瞬時頻率的均值見圖4（a）～圖4（c），結(jié)果如下所示。

圖4 IMF分量得到的瞬時頻率Fig.4 Instantaneous frequency of the IMF component

由表1可知，原信號中包含了三個振蕩模式，存在兩個局部振蕩模式和一個區(qū)域振蕩模式，這與文獻［25］中所用方法得到的結(jié)果一致。因此，本文所述方法可以很好的去除噪聲，辨識出低頻振蕩模式，并且誤差較小。

表1 取均值后頻率與理想?yún)?shù)對比Tab.1 Comparison of the frequency after averages and the ideal parameter

為進一步測試本文方法的有效性和正確性，采用Prony分析法對上述信號進行模態(tài)辨識處理。取采樣頻率為fc＝50 Hz，采樣點數(shù)N＝500，模型階數(shù)p＝18，得到Prony分析的模式辨識結(jié)果如下所示。

圖5 Prony算法Fig.5 Prony algorithm

表2 Prony方法辨識結(jié)果（頻率單位：Hz）Tab.2 Prony method identification results（frequency unit：Hz）

由圖5及表2可知，根據(jù)Prony算法辨識非線性、帶噪聲的信號時，對噪聲及非線性情況比較敏感，辨識出的結(jié)果誤差較大，無法準確的辨識出振蕩模式，而本文方法在此種情況下，適用性更好，更能較好的辨識出在帶噪聲的非線性信號的振蕩模式。但在本方法中，經(jīng)過粒子濾波及EEMD分解后，得到的瞬時頻率存在邊界效應(yīng)，存在誤差，導(dǎo)致計算精度存在不足，在后續(xù)研究中，將把減小邊界效應(yīng)的影響作為研究的方向。

6 結(jié)束語

本文提出一種基于粒子濾波與EEMD的低頻振蕩模式辨識的新方法，在對粒子濾波算法進行詳細的分析與探討之后，利用粒子濾波作為前置濾波單元，對信號進行去噪處理后，再使用EEMD對去噪后的信號進行分解得出低頻振蕩模式。從算例結(jié)果可以得出本文所提方法能夠很好的還原初始信號的特征，證明了本方法的準確性及精確性。

PF算法是一種非線性統(tǒng)計算法，能夠適用于非線性非高斯系統(tǒng)的去噪、狀態(tài)估計等，在濾波精度和穩(wěn)定度上都具有較大的優(yōu)勢，具有良好的應(yīng)用前景，粒子濾波算法計算效率的提高和有效處理邊界效應(yīng)是今后研究的方向。