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

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

0 引 言

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

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

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

1 粒子濾波算法

1.1 粒子濾波理論

粒子濾波［13-16］（Particle Filter，PF），又稱為序貴蒙特卡羅（Sequential Monte Carlo，SMC）方法，是一種基于SMC方法的貝葉斯濾波技術(shù)。其基本原理是尋找一組在狀態(tài)空間傳播的隨機(jī)粒子（樣本）描述系統(tǒng)的狀態(tài)，通過SMC法簡化貝葉斯估計(jì)中的積分運(yùn)算，進(jìn)而得出系統(tǒng)狀態(tài)的最小均方差估計(jì)。在粒子數(shù)接近無窮時(shí)能逼近服從任意概率分布的系統(tǒng)狀態(tài)。相對其它濾波技術(shù)而言，粒子濾波［17］不需要對系統(tǒng)狀態(tài)作任何先驗(yàn)性假設(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)移過程服從一階馬爾科夫模型，即時(shí)刻k的狀態(tài)xk只和前一時(shí)刻k－1的狀態(tài)xk－1有關(guān)，且模型中觀測值間相互獨(dú)立，即觀測值yk只與k時(shí)刻的狀態(tài)量xk有關(guān)。根據(jù)貝葉斯估計(jì)原理，在已知觀測量 y1：k的前提下，x1：k的后驗(yàn)概率分布為：

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

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

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

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

其中，δ（）·表示狄拉克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）。

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

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

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

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

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

2 基于EEMD的模式辨識

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

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

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

2）將xi（）t分別進(jìn)行EMD分解，得到各個(gè)IMF分量記為 Cij（t），與一個(gè)余項(xiàng)記為 rin（t），則有：

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

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

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

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

3 基于粒子濾波的EEMD分解

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

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

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

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

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

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

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

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

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

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

4 算例分析

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

該信號包含三個(gè)振蕩模式，頻率分別為0.7 Hz，1.5 Hz，0.5 Hz，該信號包含信噪比SNR為20 dB的高斯白噪聲ω（n），利用粒子濾波作為前置濾波單元，對上述信號進(jìn)行去噪處理，首先分別建立粒子濾波所需的狀態(tài)方程和觀測方程，將原始信號s（t）作為狀態(tài)方程，而觀測方程則由原始信號加上隨機(jī)噪聲得到如式（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），時(shí)間迭代步長為50，粒子數(shù) N＝1 000，定義有效區(qū)間為［5，95］，為了防止邊界效應(yīng)的影響，在下述操作過程中將區(qū)間擴(kuò)大為［0，100］，根據(jù)以上參數(shù)進(jìn)行粒子濾波去噪處理，其仿真結(jié)果如圖2所示。

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

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

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

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

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

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

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

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

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

為進(jìn)一步測試本文方法的有效性和正確性，采用Prony分析法對上述信號進(jìn)行模態(tài)辨識處理。取采樣頻率為fc＝50 Hz，采樣點(diǎn)數(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算法辨識非線性、帶噪聲的信號時(shí)，對噪聲及非線性情況比較敏感，辨識出的結(jié)果誤差較大，無法準(zhǔn)確的辨識出振蕩模式，而本文方法在此種情況下，適用性更好，更能較好的辨識出在帶噪聲的非線性信號的振蕩模式。但在本方法中，經(jīng)過粒子濾波及EEMD分解后，得到的瞬時(shí)頻率存在邊界效應(yīng)，存在誤差，導(dǎo)致計(jì)算精度存在不足，在后續(xù)研究中，將把減小邊界效應(yīng)的影響作為研究的方向。

6 結(jié)束語

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

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