基于粒子濾波與EEMD的低頻振蕩模式識別方法研究*

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

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

0 引 言

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

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

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

1 粒子濾波算法

1.1 粒子濾波理論

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

系統的狀態空間方程包含狀態方程和觀測方程，其模型描述如下：

狀態方程：

觀測方程：

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

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

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

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

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

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

1.2 重要性函數及重采樣

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

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

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

濾波前預先設定一個閾值門限Nthres，有效粒子數的值必須大于設定的閾值Nthres，否則重采樣；然后根據有效粒子數對應的歸一化權值更新粒子集。

1.3 低頻振蕩信號去噪實現

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

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

2 基于EEMD的模式辨識

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

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

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

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

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

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

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

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

3 基于粒子濾波的EEMD分解

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

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

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

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

（4）對抽取的N個樣本利用式（4）計算相對應的權值wk；

（5）將步驟3得到的權值根據式（5）進行歸一化處理，得到歸一化權值

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

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

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

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

4 算例分析

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

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

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

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

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

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

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

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

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

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

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

表1 取均值后頻率與理想參數對比Tab.1 Comparison of the frequency after averages and the ideal parameter

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

圖5 Prony算法Fig.5 Prony algorithm

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

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

6 結束語

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

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