運(yùn)用于電能質(zhì)量分析的卡爾曼濾波與仿真研究

1、運(yùn)用于電能質(zhì)量分析的卡爾曼濾波與仿真研究 姓名： 班級(jí)： 學(xué)號(hào)：運(yùn)用于電能質(zhì)量分析的卡爾曼濾波與仿真研究摘要：在現(xiàn)代電力系統(tǒng)中，由于新能源的接入，以及人們對(duì)電能質(zhì)量的要求越來(lái)越高，電能質(zhì)量分析變得越來(lái)越重要。這些變化使得系統(tǒng)狀態(tài)參數(shù)中往往混雜著噪聲，因此有必要采取適當(dāng)?shù)姆椒ǎ瑥碾S機(jī)干擾的觀測(cè)信號(hào)中提取有效的系統(tǒng)狀態(tài)參數(shù)?，F(xiàn)有運(yùn)用于電能質(zhì)量檢測(cè)的方法有很多，但由于人們對(duì)濾波過(guò)程中的的數(shù)值穩(wěn)定性、快速性、實(shí)用性和實(shí)用性提出越來(lái)越高的要求時(shí)，各種濾波方法的缺點(diǎn)越來(lái)越明顯。在此基礎(chǔ)上，卡爾曼濾波越來(lái)越受關(guān)注，并在此基礎(chǔ)上的各種改進(jìn)措施也越來(lái)越多。本文首先介紹卡爾曼濾波的基本原理，接著對(duì)運(yùn)用于電能質(zhì)量的

2、改進(jìn)的擴(kuò)展卡爾曼濾波以及無(wú)跡卡爾曼濾波進(jìn)行介紹。在比較改進(jìn)的兩種卡爾曼濾波器的利弊后，進(jìn)行仿真比較和驗(yàn)證，最后對(duì)卡爾曼濾波方法提出自己的一點(diǎn)理解與看法。關(guān)鍵詞：電能質(zhì)量，卡爾曼濾波，擴(kuò)展卡爾曼濾波，無(wú)跡卡爾曼濾波0 引言 電能質(zhì)量是指優(yōu)質(zhì)供電，合格的電能質(zhì)量是指供給敏感設(shè)備的電力的相應(yīng)的接地系統(tǒng)是適合該設(shè)備正常工作的。一個(gè)理想的三相交流電力系統(tǒng)應(yīng)該具有恒定的頻率(50Hz)和正弦的波形，A、B、C三相之間滿足相位相差120度的關(guān)系。而由于大規(guī)模電力電子設(shè)備、新能源和其他非線性元件在電力系統(tǒng)中的應(yīng)用，使得電能質(zhì)量變得越來(lái)越差，引起諸多問(wèn)題，例如：諧波，電壓驟降，電壓驟升，電壓波動(dòng)，電壓中斷，電壓

3、切痕，脈沖瞬變等。在以上列出的電能質(zhì)量問(wèn)題中，諧波是電力系統(tǒng)中最普遍存在的問(wèn)題。 美籍科學(xué)家卡爾曼在系統(tǒng)狀態(tài)空間模型的基礎(chǔ)上提出了著名的線性卡爾曼濾波器。這個(gè)系統(tǒng)可用包含正交狀態(tài)變量的微分方程模型來(lái)描述，這種濾波器是將過(guò)去的測(cè)量估計(jì)誤差合并到新的測(cè)量誤差中來(lái)估計(jì)將來(lái)的誤差?？柭鼮V波器存在的問(wèn)題是它的假設(shè)前提是線性系統(tǒng)，并在線性系統(tǒng)下設(shè)計(jì)的線性無(wú)偏、最小方差估計(jì)器。而電力系統(tǒng)是高度的非線性系統(tǒng)，故在此基礎(chǔ)上人們提出了許多改進(jìn)的卡爾曼濾波器方法，比如擴(kuò)展卡爾曼濾波器(EKF)、無(wú)跡卡爾曼濾波器(UKF)、中心差分卡爾曼濾波器(CDKF)、粒子濾波器(PF)。另外，也有人將卡爾曼濾波器與其他方法結(jié)

4、合起來(lái)，例如將采用相位修正小波變換即S變換和擴(kuò)展卡爾曼濾波相結(jié)合。 在電力系統(tǒng)中，卡爾曼也被應(yīng)用于很多方面。在電力系統(tǒng)短期負(fù)荷預(yù)測(cè)方面，采用最小二乘支持向量機(jī)(RLS-SVM)技術(shù)和卡爾曼濾波技術(shù)分別對(duì)節(jié)點(diǎn)有功分配因子和節(jié)點(diǎn)功率因數(shù)建立自適應(yīng)動(dòng)態(tài)預(yù)測(cè)模型，結(jié)果表明兩種算法滿足了系統(tǒng)運(yùn)算速度、魯棒性和預(yù)測(cè)能力的要求；在電力系統(tǒng)動(dòng)態(tài)狀態(tài)估計(jì)方面，人們提出了混合卡爾曼粒子濾波器(MKPF)，該濾波器在電力系統(tǒng)在受到擾動(dòng)之后，可以很好地收斂于真實(shí)值，和EKF和UKF相比具有更高的精度和穩(wěn)定性；在電能質(zhì)量分析方面，卡爾曼濾波尤其是擴(kuò)展卡爾曼濾波和無(wú)跡卡爾曼濾波在頻率跟蹤領(lǐng)域得到了廣泛的應(yīng)用，在傅里葉變換

5、方法、最小絕對(duì)值狀態(tài)估計(jì)、遺傳算法、小波分析、Teager 能量算子等方法中脫穎而出；在電力系統(tǒng)繼電保護(hù)方面，卡爾曼濾波具有響應(yīng)速度快、濾波效果好的特點(diǎn)，能夠滿足繼電保護(hù)快速跳閘和準(zhǔn)確動(dòng)作的要求，從而在電力系統(tǒng)繼電保護(hù)中呈現(xiàn)出了很好的應(yīng)用前景；在風(fēng)電場(chǎng)風(fēng)速預(yù)測(cè)中方面，通過(guò)卡爾曼濾波修正的風(fēng)速數(shù)據(jù)能夠很好地跟蹤實(shí)際風(fēng)速數(shù)據(jù)的變化趨勢(shì)，預(yù)測(cè)精度得到了明顯提高；在電機(jī)狀態(tài)和參數(shù)估計(jì)方面，卡爾曼濾波在同步電機(jī)尤其是永磁同步電機(jī)(PMSM)中發(fā)揮了極致作用。 本文著重從卡爾曼濾波原理、分類及其在電能質(zhì)量分析中的應(yīng)用進(jìn)展進(jìn)行系統(tǒng)的介紹和分析，并對(duì)卡爾曼濾波器存在的可能研究的問(wèn)題提出自己的看法。1 常規(guī)卡爾

6、曼濾波 根據(jù)經(jīng)典物理學(xué)可知，在沒(méi)有外部干擾時(shí)，一個(gè)系統(tǒng)未來(lái)的狀態(tài)可根據(jù)目前的狀態(tài)從已知的運(yùn)動(dòng)方程中確定出來(lái)。但任何實(shí)際的系統(tǒng)總是存在外部干擾，人們對(duì)其運(yùn)動(dòng)方程的描述不可能很準(zhǔn)確?？梢哉J(rèn)為：任何實(shí)際物理系統(tǒng)的行為都由2部分組成，一部分是根據(jù)已知運(yùn)動(dòng)方程正確地預(yù)測(cè)出來(lái)，另一部分是均方為零的隨機(jī)分量。 在分析期望響應(yīng)存在的情況下的線性最優(yōu)濾波器，我們由Wiener濾波器。在期望響應(yīng)未知，線性最優(yōu)濾波器即是Kalman濾波器。Kalman濾波理論是Wiener濾波理論的發(fā)展，它最早用于隨機(jī)過(guò)程的參數(shù)估計(jì)，后來(lái)很快在各種最優(yōu)濾波的最優(yōu)控制問(wèn)題中得到了極其廣泛的應(yīng)用。Kalman濾波器具有以下特點(diǎn)：(1)

7、其數(shù)學(xué)公式用狀態(tài)空間概念描述；(2)它的解是遞推計(jì)算的，級(jí)與Wiener濾波器不同，Kalman濾波器是一種自適應(yīng)濾波器。下面對(duì)其進(jìn)行簡(jiǎn)單介紹。 假設(shè)線性離散方程為 (1) (2)式中：為狀態(tài)向量；為量測(cè)向量；為系統(tǒng)噪聲向量；為量測(cè)噪聲向量；為狀態(tài)轉(zhuǎn)移矩陣；為量測(cè)轉(zhuǎn)移矩陣。假設(shè)系統(tǒng)噪聲和量測(cè)噪聲是互不相關(guān)的零均值高斯白噪聲，方差陣分別為、，為了要遞推地在每次獲得觀測(cè)量之后估計(jì)狀態(tài)量。定義狀態(tài)量的一步預(yù)測(cè)為，其他類推，則卡爾曼濾波遞推方程如下。狀態(tài)1步預(yù)測(cè)為 （3）1步預(yù)測(cè)誤差方差陣為 (4)狀態(tài)估計(jì)為 (5)估計(jì)誤差方差陣為 (6)濾波增益矩陣為 (7)式中為單位陣。式(3)(7)就是隨機(jī)線性

8、離散系統(tǒng)卡爾曼濾波的基本方程。只要給定初值和及、，根據(jù)時(shí)刻的觀測(cè)值就可以遞推計(jì)算得時(shí)刻的狀態(tài)估計(jì)（=1,2，···）?？柭鼮V波有如下特點(diǎn)： 由于Kalman濾波算法將被估計(jì)的信號(hào)看作在白噪聲作用下的一個(gè)隨機(jī)線性系統(tǒng)的輸出,并且其輸入輸出關(guān)系是由狀態(tài)方程和輸出方程在時(shí)間域內(nèi)給出的,因此這種濾波方法僅適用于平穩(wěn)序列的濾波,而且特別適用于非平穩(wěn)或平穩(wěn)馬爾可夫序列或高斯一馬爾可夫序列的濾波,因此其應(yīng)用范圍是十分廣泛的。 由于濾波的基本方程是時(shí)間域內(nèi)的遞推形式,其計(jì)算過(guò)程是一個(gè)不斷的“預(yù)測(cè)一修正”過(guò)程,在求解時(shí)不要求存儲(chǔ)大量數(shù)據(jù),并且一旦觀測(cè)到了新的數(shù)據(jù),隨時(shí)可以算得新的

9、濾波值,因此這種濾波方法非常便于實(shí)時(shí)處理和計(jì)算機(jī)實(shí)現(xiàn)。 由于濾波器的增益矩陣與觀測(cè)無(wú)關(guān),因此它可預(yù)先離線算出,從而可以減少實(shí)時(shí)在線計(jì)算量;在求濾波器增益矩陣,時(shí),要求一個(gè)矩陣的逆,即要計(jì)算，它的階數(shù)只取決于觀測(cè)方程的維數(shù)m,而m通常是很小的。這樣,上面的求逆運(yùn)算是比較方便的;另外,在求解濾波器增益的過(guò)程中, 隨時(shí)可以算得濾波器的精度指標(biāo),其對(duì)角線上的元素就是濾波誤差向量各分量的方差。 增益矩陣與初始方差陣、系統(tǒng)噪聲方差陣、以及觀測(cè)噪聲方差陣之間具有如下關(guān)系：·，和同乘一個(gè)相同的標(biāo)量時(shí)，值不變。·當(dāng)增大時(shí), 就變小,這在直觀上是很容易理解的，因?yàn)槿绻^測(cè)噪聲增大,那么濾波增益

10、就應(yīng)取小一些（因?yàn)檫@時(shí)的新信息里的誤差比較大），以減弱觀測(cè)噪聲對(duì)濾波值的影響·如果變小, 變小,或兩者都變小，則變小，也變小，從而變小。這也是很自然的。因?yàn)樽冃?表示初始估計(jì)較好, 變小,表示系統(tǒng)噪聲變小,于是增益矩陣也應(yīng)小些以便給于較小的修正。綜上所述,可以簡(jiǎn)單地說(shuō),增益矩陣與成正比,而與成反比。通過(guò)仿真附錄一中的卡爾曼濾波實(shí)驗(yàn)程序，我們可得如下圖形:. 和的圖： . 真實(shí)值與估計(jì)值的比較圖 . 各自的誤差圖 .通過(guò)用卡爾曼濾波器的狀態(tài)誤差協(xié)方差矩陣畫(huà)出的和：分析：中的方差是中的誤差平方后取均值，是均方誤差。誤差直接由真實(shí)值減去估計(jì)值，有正有負(fù)，而均方誤差沒(méi)有這個(gè)缺陷，更能綜合的表

11、示濾波的效果。由圖(2)可知，從真實(shí)值與估計(jì)值的比較可以看出，卡爾曼濾波器還是能夠有效的跟蹤真實(shí)值，但是由圖(3)明顯看出，在跟蹤的過(guò)程中，誤差在零值上下擺動(dòng)，存在不穩(wěn)定的問(wèn)題。真實(shí)值與估計(jì)值之間的誤差幅值一直處于抖動(dòng)狀態(tài)，由此，我們可以得出，運(yùn)用卡爾曼濾波，雖然能夠?qū)崿F(xiàn)運(yùn)用較少的信息量實(shí)現(xiàn)平穩(wěn)與非平穩(wěn)過(guò)程的追蹤，但也有當(dāng)運(yùn)動(dòng)目標(biāo)長(zhǎng)時(shí)間被遮擋時(shí)會(huì)存在目標(biāo)跟蹤丟失的情況 。2. 擴(kuò)展卡爾曼濾波由于在實(shí)際中廣泛存在的是非線性狀態(tài)空間模型，使得常規(guī)卡爾曼濾波在電能質(zhì)量分析中的應(yīng)用存在困難，于是便出現(xiàn)了諸多針對(duì)非線性模型的次優(yōu)方法，其中應(yīng)用最廣泛的是擴(kuò)展卡爾曼濾波(extended Kalman fi

12、ltering，EKF)。EKF 是將非線性系統(tǒng)線性化，與線性卡爾曼濾波公式完全類似。其主要思想是對(duì)非線性函數(shù)的泰勒展開(kāi)式進(jìn)行截?cái)?，?shí)現(xiàn)非線性函數(shù)的線性化。根據(jù)泰勒展開(kāi)式進(jìn)行的是1階還是2階截取，EKF主要分為1階EKF(first order EKF)和2階EKF(second order EKF)。電能質(zhì)量分析中最常用的是1階EKF，原理簡(jiǎn)述如下。假如非線性系統(tǒng)可表示為 （8） （9）式中： 為系統(tǒng)狀態(tài)向量；為系統(tǒng)量測(cè)向量；和 是關(guān)于狀態(tài)的非線性函數(shù)； 和 均是均值為零的高斯白噪聲。式(8)、(9)分別是狀態(tài)方程和量測(cè)方程。為了使卡爾曼濾波應(yīng)用到非線性系統(tǒng)中，非線性系統(tǒng)必須在指定位置進(jìn)行泰

13、勒展開(kāi)，實(shí)現(xiàn)線性化。推導(dǎo)過(guò)程如下：利用泰勒公式，分別在和處對(duì)狀態(tài)方程和觀測(cè)方程進(jìn)行1階泰勒展開(kāi)，可得 （10） （11）假設(shè) （12） （13） （14） （15）則式(10)、(11)可改寫(xiě)成與常規(guī)卡爾曼方程相似的形式： （16） （17）1階EKF遞推方程組與常規(guī)卡爾曼濾波遞推方程組在形式上相同，不同的是：中的 和 被 1階EKF中的Jacobian矩陣和代替，并且預(yù)測(cè)平均值和預(yù)測(cè)的冗余在EKF中也分別計(jì)算，其遞推方程與卡爾曼濾波相同。在電能質(zhì)量分析中，、 矩陣的設(shè)計(jì)略有不同。1991 年，Beides H M 和Heydt G T 提出用擴(kuò)展卡爾曼濾波獲得電力系統(tǒng)諧波的動(dòng)態(tài)狀態(tài)估計(jì)，經(jīng)過(guò)

14、實(shí)驗(yàn)室仿真和實(shí)測(cè)試驗(yàn)證明擴(kuò)展卡爾曼濾波能動(dòng)態(tài)地追蹤諧波內(nèi)容和時(shí)間。1993 年，Kamwa也將 EKF引入電力系統(tǒng)電能質(zhì)量分析中，用于測(cè)量閃變。基于卡爾曼濾波的實(shí)時(shí)監(jiān)控電壓閃變的算法，建立的模型是具有線性時(shí)不變過(guò)渡矩陣的系統(tǒng)電壓隨機(jī)狀態(tài)空間，用 EKF來(lái)測(cè)量帶噪聲的單相電壓。一種測(cè)量 50/60Hz低頻信號(hào)調(diào)制的卡爾曼濾波方法，用于確定信號(hào)和隨機(jī)信號(hào)的調(diào)制，產(chǎn)生的參數(shù)將電壓閃變控制在一個(gè)可接受的范圍內(nèi)，這種方法運(yùn)用了擴(kuò)展卡爾曼濾波模型，經(jīng)過(guò)典型事例的試驗(yàn)證明了其測(cè)量電壓閃變的準(zhǔn)確性和可靠性，并具有預(yù)測(cè)電壓閃變的功能。一種基于 13狀態(tài)系統(tǒng)的電壓暫降的擴(kuò)展卡爾曼濾波檢測(cè)和分析方法，實(shí)驗(yàn)室仿真和實(shí)

15、際電壓暫降仿真表明，這種方法效果明顯優(yōu)于IEC 電能質(zhì)量標(biāo)準(zhǔn)的RMS方法。雖然擴(kuò)展卡爾曼濾波有很好的發(fā)展前景，但它在實(shí)際應(yīng)用中存在明顯的缺陷：一是線性化有可能產(chǎn)生極不穩(wěn)定的濾波；二是EKF需要計(jì)算Jacobian矩陣的導(dǎo)數(shù)，實(shí)現(xiàn)起來(lái)較為復(fù)雜，而對(duì)于一些不可微的情況，EKF可能失效。在模型非線性較強(qiáng)以及系統(tǒng)噪聲非高斯時(shí)，估計(jì)的精度嚴(yán)重降低，甚至?xí)斐蔀V波器發(fā)散。近幾年，更多的人針對(duì)EKF對(duì)諧波的跟蹤和檢測(cè)提出了更好的改進(jìn)方法，其中最重要的是擴(kuò)展復(fù)卡爾曼濾波器或復(fù)數(shù)型擴(kuò)展卡爾曼濾波(extended complex Kalman filter，ECKF)。一種用于估計(jì)混有噪聲的信號(hào)頻率的ECKF。

16、由于復(fù)數(shù)濾波器公式需要多傳感器來(lái)測(cè)量三相信號(hào)，通過(guò)重置由輸出誤差幅度決定的協(xié)方差矩陣Q，使設(shè)計(jì)的模型仿真不僅考慮了頻率、幅值和相位的測(cè)量，也考慮到了實(shí)際復(fù)雜的網(wǎng)絡(luò)情況，實(shí)驗(yàn)結(jié)果表明其性能優(yōu)于普通的擴(kuò)展卡爾曼濾波器。一種 PEKF(the proposed variant of the extended Kalman filtering)用于測(cè)量電力系統(tǒng)諧波，并在數(shù)學(xué)上證明了其穩(wěn)定性，與 ADALINE (adaptive linear combiner)、ANF(adaptive notch filter)、MFT(multi frequency tracker)以及現(xiàn)有的常規(guī)EKF(ordi

17、nary extended Kalman filter，OEKF)和ECKF 相比，試驗(yàn)表明其具有更好的性能和更小的計(jì)算量，更適合實(shí)時(shí)實(shí)現(xiàn)。Huang 提出一種基于ECKF 的魯棒性算法來(lái)提高電力系統(tǒng)的頻率估計(jì)效果，利用在壓縮數(shù)據(jù)中產(chǎn)生的不良數(shù)據(jù)來(lái)提高頻率估計(jì)效率，并與ECKF算法進(jìn)行比較，證明其提高了系統(tǒng)執(zhí)行效率。3. 無(wú)跡卡爾曼濾波 為了更精確地?cái)M合非線性函數(shù)，Julier 提出了(unscented transformation，UT)和無(wú)跡卡爾曼濾波(unscented Kalman filtering，UKF)。通過(guò)結(jié)合無(wú)跡變換和無(wú)跡卡爾曼濾波來(lái)實(shí)現(xiàn)非線性系統(tǒng)的狀態(tài)估計(jì)，是近年來(lái)用于

18、解決該問(wèn)題的一種新的熱點(diǎn)方法。它通過(guò)一組精確選擇的 Sigma 點(diǎn)來(lái)匹配隨機(jī)量的統(tǒng)計(jì)特性，UKF 沒(méi)有涉及非線性映射函數(shù)的Jacobian矩陣計(jì)算問(wèn)題，從而使算法的實(shí)現(xiàn)EKF更為容易, 在保持相當(dāng)運(yùn)算量的同時(shí)，具有更高的估計(jì)精度和更廣泛的適用范圍。傳統(tǒng)的線性化方法是對(duì)非線性映射本身做某種線性近似，然后再應(yīng)用線性估計(jì)的各種方法。而Julier S. J.等人提出的無(wú)跡變換則是基于用有限的參數(shù)來(lái)近似隨機(jī)量的概率統(tǒng)計(jì)特性要比近似任意的非線性映射函數(shù)更為容易的思想。無(wú)跡變換的基本步驟可概括為：關(guān)于 的Sigma點(diǎn)集的產(chǎn)生不確定性的非線性變換與傳遞關(guān)于的統(tǒng)計(jì)特性的推算。無(wú)跡變換 Sigma點(diǎn)集的選取方式

19、不同，會(huì)產(chǎn)生很多種變換的演變形式，其目的主要是進(jìn)一步提高變換的精度，增強(qiáng)算法的穩(wěn)定性和減小運(yùn)算量等。在應(yīng)用UKF時(shí)首先要對(duì)狀態(tài)量進(jìn)行擴(kuò)展，也就是將模型噪聲也作為狀態(tài)量的一部分，相應(yīng)地，無(wú)跡變換中用到的Sigma點(diǎn)也需要擴(kuò)展，具體表示如下。擴(kuò)展?fàn)顟B(tài)方程的初始值： （16） （17）式中：為模型初始狀態(tài)變量；和 分別為擴(kuò)展?fàn)顟B(tài)向量的均值和協(xié)方差陣。Sigma點(diǎn)集的創(chuàng)建通過(guò)下式實(shí)現(xiàn)： （18）式中：為預(yù)測(cè)空間維數(shù)；表示矩陣平方根的第 個(gè)行向量或列向量，而矩陣平方根的常見(jiàn)求法是采用Cholesky 分解；為點(diǎn)處的協(xié)方差；，決定Sigma取的點(diǎn)數(shù)，是由和參數(shù)決定的函數(shù)，為控制周?chē)母唠A非線性值的參數(shù)，是

20、介于 0.00011之間的一個(gè)常數(shù)，是次要的比例參數(shù)，通常設(shè)置為 0 或3n ，以確保Sigma點(diǎn)分布的峭度與高斯分布的峭度一致。Sigma點(diǎn)向量通過(guò)狀態(tài)方程的非線性影射得到： （19）擴(kuò)展?fàn)顟B(tài)量的1步預(yù)測(cè)為 （20）擴(kuò)展?fàn)顟B(tài)量1步預(yù)測(cè)的協(xié)方差陣為 （21）然后計(jì)算量測(cè)空間。量測(cè)空間的Sigma點(diǎn)集的創(chuàng)建通過(guò)下式實(shí)現(xiàn)： （22）1 步預(yù)測(cè)的擴(kuò)展?fàn)顟B(tài)Sigma點(diǎn)向量經(jīng)過(guò)觀測(cè)方程的非線性映射得到： （23）觀測(cè)量的1步預(yù)測(cè)為 （24）式中：；為合并高階狀態(tài)分布的先驗(yàn)知識(shí)，高斯分布的最佳選擇是2。 觀測(cè)量1步預(yù)測(cè)的協(xié)方差陣為 （25）1 步預(yù)測(cè)的互協(xié)方差陣為 （26）最后得到量測(cè)更新過(guò)程。更新誤差協(xié)

21、方差： （27）計(jì)算濾波增益矩陣： （28）更新估計(jì)： （29）UKF在保持相當(dāng)運(yùn)算量的同時(shí)具有更高的估計(jì)精度和更為廣泛的適用范圍，它在國(guó)內(nèi)的相關(guān)研究起步較晚，但發(fā)展很快?？刹榈墓_(kāi)資料主要集中于最近的23年內(nèi)，而且在電力系統(tǒng)尤其是電能質(zhì)量方面的研究成果比較少。電能質(zhì)量幅度和頻率估計(jì)方面比較了無(wú)跡卡爾曼濾波和擴(kuò)展卡爾曼濾波的差異。仿真試驗(yàn)表明，在相等的噪聲條件下，無(wú)跡卡爾曼濾波與擴(kuò)展卡爾曼濾波的計(jì)算精度相似，但無(wú)跡卡爾曼濾波的計(jì)算量更小。羅諶持、張明應(yīng)用無(wú)跡變換，將復(fù)數(shù)型Sigma點(diǎn)卡爾曼濾波(complex Sigma point Kalman filtering，CSPKF)中的擴(kuò)展復(fù)卡爾

22、曼濾波算法進(jìn)行比較，通過(guò)變換，首先將三相電壓信號(hào)轉(zhuǎn)換成復(fù)電壓信號(hào)，再利用CSPKF算法對(duì)發(fā)生諧波畸變和隨機(jī)噪聲干擾的電力系統(tǒng)電壓信號(hào)的頻率進(jìn)行動(dòng)態(tài)估計(jì)和跟蹤的過(guò)程進(jìn)行改進(jìn)。算法仿真表明，CSPKF 算法具有優(yōu)異的動(dòng)態(tài)跟蹤性能，迅速跟蹤頻率和幅值變化的同時(shí)又保持了較低的跟蹤誤差。4. 應(yīng)用現(xiàn)狀隨著卡爾曼濾波的快速發(fā)展，除了上述3種形式及其改進(jìn)算法以外，還出現(xiàn)了其他一些形式，例如集合變換卡爾曼濾波(the ensemble transform Kalman filtering，ETKF)也在現(xiàn)代工程應(yīng)用中發(fā)揮了越來(lái)越重要的作用。就近十幾年的發(fā)展而言，卡爾曼濾波在電能質(zhì)量中的主流發(fā)展主要集中于對(duì) E

23、KF 算法進(jìn)行改進(jìn)、UKF 的改進(jìn)算法研究以及與其他方法相結(jié)合進(jìn)行電能質(zhì)量分析。方程中一些初始量(如噪聲協(xié)方差矩陣 Q、濾波增益矩陣 K 等)的選擇關(guān)系到動(dòng)態(tài)追蹤特性。一種自適應(yīng)卡爾曼濾波，該方法的狀態(tài)模型是線性的，頻率作為狀態(tài)變量與諧波母線電壓獨(dú)立。它利用統(tǒng)計(jì)規(guī)則在2個(gè)基本噪聲協(xié)方差矩陣Q 模式間進(jìn)行切換，取代了直接優(yōu)化的Q模式來(lái)進(jìn)行動(dòng)態(tài)諧波狀態(tài)估計(jì)。采用這種方法可以重置卡爾曼濾波增益K，從而更好地在暫態(tài)過(guò)程中實(shí)現(xiàn)快速追蹤，并在新西蘭220k V電網(wǎng)上進(jìn)行仿真，證明了方法的可行性。一種卡爾曼濾波的自校正方法來(lái)測(cè)量諧波，并與噪聲協(xié)方差矩陣 Q 為定值零的方法進(jìn)行比較，證明了其追蹤信號(hào)波動(dòng)的優(yōu)越

24、性。 除了單獨(dú)采用卡爾曼濾波進(jìn)行電能質(zhì)量分析外，也有人提出將其與其他方法相結(jié)合來(lái)檢測(cè)電能質(zhì)量擾動(dòng)。一種新的檢測(cè)、定位、分類短期擾動(dòng)的方法，采用相位修正小波變換即S 變換和擴(kuò)展卡爾曼濾波相結(jié)合。S 變換有很好的時(shí)頻特性，提供適合于電力系統(tǒng)自動(dòng)識(shí)別的檢測(cè)、定位和視覺(jué)識(shí)別，而EKF提供頻繁出現(xiàn)的頻率短期擾動(dòng)的自動(dòng)分類和測(cè)量。實(shí)驗(yàn)證明，結(jié)合ST和EKF可以完整地分類和測(cè)量電能質(zhì)量短期擾動(dòng)，在有噪聲的情況下具有很高的精度。小波與卡爾曼濾波相結(jié)合的方法來(lái)自動(dòng)檢測(cè)和分析電壓暫降。小波用來(lái)檢測(cè)電壓波形和更好地估算電壓驟降的時(shí)間相關(guān)的參數(shù)，而卡爾曼濾波在這一過(guò)程中用來(lái)估計(jì)電壓幅值和相位角，這種方法可用來(lái)處理實(shí)時(shí)

25、系統(tǒng)，并且已通過(guò)計(jì)算機(jī)仿真和實(shí)驗(yàn)對(duì)低壓配電系統(tǒng)的測(cè)試，證明了其測(cè)量電壓暫降的完整性和準(zhǔn)確性。Kumar A、Das B 和Sharma J 提出了魯棒擴(kuò)展卡爾曼濾波(robust extended Kalman filtering，REKF)技術(shù)來(lái)進(jìn)行動(dòng)態(tài)諧波狀態(tài)估計(jì)，人工神經(jīng)網(wǎng)絡(luò)和改進(jìn)的EKF相結(jié)合的方法，人工神經(jīng)網(wǎng)絡(luò)用于測(cè)量諧波儀的位置和提供虛擬測(cè)量(pseudo-measurement)的粗略估計(jì)，REKF 用于進(jìn)行動(dòng)態(tài)諧波狀態(tài)估計(jì)。文中通過(guò)正常負(fù)荷 情況、負(fù)荷突然變化和存在不良數(shù)據(jù)的IEEE-14 系統(tǒng)測(cè)量這3組實(shí)驗(yàn)，將REKF與EKF進(jìn)行比較，證明REKF的性能優(yōu)于EKF。狀態(tài)估計(jì)器

26、即使在存在不良數(shù)據(jù)時(shí)也能實(shí)時(shí)估計(jì)電力系統(tǒng)諧波狀態(tài)，從IEEE-14節(jié)點(diǎn)系統(tǒng)仿真結(jié)果來(lái)看，估計(jì)結(jié)果與真值相近。近年來(lái)，無(wú)跡變換和無(wú)跡卡爾曼濾波的研究取得了一些進(jìn)展。無(wú)跡變換主要有：對(duì)稱 Sigma點(diǎn)集無(wú)跡變換(symmetric Sigma set UT)、單形 Sigma點(diǎn)集無(wú)跡變換(simplex Sigma set UT)、變尺度無(wú)跡變換(scaled UT)、高階無(wú)跡變換(high order UT)等?？梢愿鶕?jù)具體的場(chǎng)合提出特定分布的無(wú)跡變換的衍變形式，當(dāng)然，通常只能根據(jù)特定的場(chǎng)合達(dá)到其中的部分目標(biāo)。也可引入代價(jià)函數(shù)或懲罰函數(shù)來(lái)盡可能滿足某種條件而不要求一定能夠達(dá)到?？紤]到提高算法的數(shù)

27、值穩(wěn)定性和減小算法運(yùn)算量，Merwe 和Wan借用平方根濾波(square-root filtering)，提出應(yīng)用QR分解和Cholesky分解更新的方法直接非線性傳遞更新協(xié)方差陣的平方根，命名為平方根UKF(square-root UKF，SR-UKF)。5.實(shí)驗(yàn)驗(yàn)證與分析 通過(guò)附錄二的程序，我們可以得到如下圖形：(5) EKF、UKF估計(jì)值與真實(shí)值比較 由圖可以看出，擴(kuò)展卡爾曼濾波器和無(wú)跡卡爾曼濾波器都能跟蹤真實(shí)值，且都存在一定的跟蹤誤差，從兩者的整體趨勢(shì)來(lái)看，無(wú)跡卡爾曼除了在40步長(zhǎng)等少數(shù)點(diǎn)誤差較大外，其余的跟蹤點(diǎn)都比較準(zhǔn)確，所以無(wú)跡卡爾曼具有比擴(kuò)展卡爾曼跟優(yōu)良的跟蹤性能，與前文理論部

28、分符合。(6) EKF估計(jì)值與真實(shí)值及95%置信區(qū)間比較 在此，我們單獨(dú)將EKF估計(jì)值與真實(shí)值和95%置信區(qū)間相比較，由此可以看出，雖然EKF大致能跟蹤實(shí)際真實(shí)值，但明顯有很大一部分?jǐn)?shù)據(jù)沒(méi)有跟蹤到，存在一定誤差。這里，雖然將誤差放寬到95%，只有少部分點(diǎn)在誤差區(qū)間，還有很大一部分點(diǎn)在置信區(qū)間外，說(shuō)明，在卡爾曼濾波器基礎(chǔ)上改進(jìn)的EKF濾波器，仍然存在誤差較大問(wèn)題，雖然已經(jīng)能運(yùn)用在電能質(zhì)量這種非線性系統(tǒng)中。(7) UKF估計(jì)值與真實(shí)值及95%置信區(qū)間比較 圖(7)中，我們可以看出，相比于圖(6)，UKF有明顯的改善。UKF不僅能跟蹤非線性系統(tǒng)，而且在此圖中，只有幾個(gè)點(diǎn)在95%置信區(qū)間外，并且都偏離

29、置信區(qū)間曲線不遠(yuǎn)。由此，我們可以得出結(jié)論，UKF比EKF有著跟好的濾波性能。(8) EKF與UKF估計(jì)值的誤差比較 經(jīng)過(guò)圖(6)和圖(7)的比較，我們能夠得出初步的結(jié)論，這里，在圖(8)，我們將EKF和UKF的估計(jì)值誤差擺放在一起，我們可以更好地得出我們的結(jié)論：UKF比EKF在電能質(zhì)量分析中擁有更好的性能。6. 結(jié)論 卡爾曼濾波是由卡爾曼在 20世紀(jì)60年代提出的基于貝葉斯濾波原理，用一個(gè)狀態(tài)方程和一個(gè)量測(cè)方程來(lái)完整地描述線性動(dòng)態(tài)過(guò)程的方法。擴(kuò)展卡爾曼濾波算法中傳統(tǒng)的解決方法是借助 Taylor 級(jí)數(shù)展開(kāi)，其后也有很多改進(jìn)算法。而無(wú)跡變換和無(wú)跡卡爾曼濾波是通過(guò)一組精確選擇的 Sigma 點(diǎn)來(lái)匹

30、配隨機(jī)量的統(tǒng)計(jì)特性，易于實(shí)現(xiàn)，目前雖取得了一些理論成果，但工程應(yīng)用較少。因此，卡爾曼濾波的理論研究和分析方法仍有發(fā)展空間。 1) 算法的運(yùn)算量和估算精度的平衡。影響算法運(yùn)算量的一個(gè)重要因素是 Sigma點(diǎn)集中 Sigma點(diǎn)的數(shù)目，而估算精度是UT 和UKF的主要優(yōu)勢(shì)，如何平衡最合適的Sigma點(diǎn)數(shù)目很重要。2) 算法中參數(shù)的智能選取問(wèn)題。UT 中涉及到的很多智能參數(shù)很有靈活性，其選取影響到程序的穩(wěn)定性，甚至導(dǎo)致程序崩潰。因此，智能參數(shù)的選取顯得至關(guān)重要。3) 算法的實(shí)際應(yīng)用問(wèn)題。針對(duì)電能質(zhì)量中信號(hào)處理的一些特點(diǎn)，可以針對(duì) UT 和 UKF 進(jìn)行一些改進(jìn)，而如何改進(jìn)也顯得十分重要。4) 保證卡爾

31、曼濾波算法最優(yōu)性的一個(gè)先決條件是建立準(zhǔn)確的系統(tǒng)方程和量測(cè)方程，這要求清晰了解系統(tǒng)的模型。偏離理想的系統(tǒng)模型或者動(dòng)態(tài)模型誤差在濾波時(shí)未模型化，這都必然會(huì)對(duì)濾波結(jié)果造成影響。另外，如何檢測(cè)并修復(fù)模型誤差，以使模型最佳。這都仍有待進(jìn)一步的研究。5) 保證卡爾曼濾波算法最優(yōu)性的另一個(gè)先決條件是正確描述系統(tǒng)噪聲和量測(cè)噪聲的統(tǒng)計(jì)特性。卡爾曼算法中要求兩噪聲都是白噪聲且互不相關(guān)，但在實(shí)際應(yīng)用中，系統(tǒng)噪聲和量測(cè)噪聲極有可能不是白噪聲或二者相關(guān)。即有色噪聲或者相關(guān)噪聲條件下的卡爾曼濾波算法是可以繼續(xù)加深研究的熱點(diǎn)。附錄一 卡爾曼濾波改進(jìn)程序%卡爾曼濾波實(shí)驗(yàn)程序clc;y1=3.29691969,3.387365

32、15,7.02830641,9.71212521,11.42018315,15.97870583,22.06934285,28.30212781,30.44683831,38.75875595; %觀測(cè)值y1(k) y2=2.10134294,0.47540797,3.17688898,2.49811140,2.91992424,6.17307616,5.42519274,3.05365741,5.98051141,4.51016361; %觀測(cè)值y2(k) p0=1,0;0,1;p=p0; %均方誤差陣賦初值A(chǔ)k=1,1;0,1; %轉(zhuǎn)移矩陣Qk=1,0;0,1; %系統(tǒng)噪聲矩陣Ck=1,0

33、;0,1; %量測(cè)矩陣Rk=1,0;0,2; %測(cè)量噪聲矩陣x0=0,0'xk=x0; %狀態(tài)矩陣賦初值for k=1:10 Pk=Ak*p*Ak'+Qk; %濾波方程3 Hk=Pk*Ck'*inv(Ck*Pk*Ck'+Rk); %濾波方程2 yk=y1(k);y2(k); %觀測(cè)值xk=Ak*xk+Hk*(yk-Ck*Ak*xk); %濾波方程1 x1(k)=xk(1); x2(k)=xk(2); %記錄估計(jì)值p=(eye(2)-Hk*Ck)*Pk; %濾波方程4 pk(:,k)=p(1,1),p(2,2)' %記錄狀態(tài)誤差協(xié)方差矩陣end figur

34、e %畫(huà)圖表示狀態(tài)矢量的估計(jì)值subplot(2,1,1) i=1:10; plot(i,x1(i),'k') h=legend('x1(k)的估計(jì)值') set(h,'interpreter','none') subplot(2,1,2) i=1:10; plot(i,x2(i),'k') h=legend('x2(k)的估計(jì)值') set(h,'interpreter','none') X1=0,1.65428714,3.50300702,5.997852924,

35、9.15040740,12.50873910,16.92192594,21.34483352,25.89335144,31.54135330,36.93605670; %由模擬得到的實(shí)際狀態(tài)值X1(k) X2=0,1.65428714,1.84871988,2.47552222,3.17187816,3.35833170,4.41318684,4.42290758,4.54851792,5.64800186,5.394470340; %由模擬得到的實(shí)際狀態(tài)值X2(k) figure %在同一幅圖中畫(huà)出狀態(tài)矢量的估計(jì)值與真實(shí)值subplot(2,1,1) i=1:10; plot(i,x1(i)

36、,'k',i,X1(i+1),'b') h=legend('x1(k)的估計(jì)值','x1(k)的真實(shí)值') set(h,'interpreter','none') subplot(2,1,2) i=1:10; plot(i,x2(i),'k',i,X2(i+1),'b') h=legend('x2(k)的估計(jì)值','x2(k)的真實(shí)值') set(h,'interpreter','none') for i

37、=1:10 %計(jì)算x(k)的誤差e1(i)=X1(i+1)-x1(i); e2(i)=X2(i+1)-x2(i); end figure %畫(huà)出誤差圖subplot(2,1,1) i=1:10; plot(i,e1(i),'r') h=legend('x1(k)的誤差') set(h,'interpreter','none') subplot(2,1,2) i=1:10; plot(i,e2(i),'r') h=legend('x2(k)的誤差') set(h,'interpreter

38、9;,'none') figure %通過(guò)用爾曼濾波器的狀態(tài)誤差協(xié)方差矩陣畫(huà)出E1(k/k)2和E2(k/k)2i=1:10; subplot(2,1,1) plot(i,pk(1,i),'r') h= legend('由狀態(tài)誤差協(xié)方差矩陣得到的E1(k/k)2')set(h,'Interpreter','none') subplot(2,1,2) plot(i,pk(2,i),'r') h= legend('由狀態(tài)誤差協(xié)方差矩陣得到的E2(k/k)2')set(h,'Int

39、erpreter','none') 附錄二% EKF UKF 的兩個(gè)算法的仿真與比較clear;% tic;x = 0.1; % 初始狀態(tài) x_estimate = 1;%狀態(tài)的估計(jì)e_x_estimate = x_estimate; %EKF的初始估計(jì)u_x_estimate = x_estimate; %UKF的初始估計(jì)Q = 10;%input('請(qǐng)輸入過(guò)程噪聲方差Q的值: '); % 過(guò)程狀態(tài)協(xié)方差 R = 1;%input('請(qǐng)輸入測(cè)量噪聲方差R的值: '); % 測(cè)量噪聲協(xié)方差 P =5;%初始估計(jì)方差e_P = P; %EK

40、F方差u_P = P;%UKF方差tf = 50; % 模擬長(zhǎng)度 x_array = x;%真實(shí)值數(shù)組e_x_estimate_array = e_x_estimate;%EKF最優(yōu)估計(jì)值數(shù)組u_x_estimate_array = u_x_estimate;%UKF最優(yōu)估計(jì)值數(shù)組u_k = 1; %微調(diào)參數(shù)u_symmetry_number = 4; % 對(duì)稱的點(diǎn)的個(gè)數(shù)u_total_number = 2 * u_symmetry_number + 1; %總的采樣點(diǎn)的個(gè)數(shù)linear = 0.5;close all;for k = 1 : tf % 模擬系統(tǒng) x = linear * x

41、+ (25 * x / (1 + x2) + 8 * cos(1.2*(k-1) + sqrt(Q) * randn; %狀態(tài)值 y = (x2 / 20) + sqrt(R) * randn; %觀測(cè)值 %擴(kuò)展卡爾曼濾波器 %進(jìn)行估計(jì) 第一階段的估計(jì) e_x_estimate_1 = linear * e_x_estimate + 25 * e_x_estimate /(1+e_x_estimate2) + 8 * cos(1.2*(k-1); e_y_estimate = (e_x_estimate_1)2/20; %這是根據(jù)k=1時(shí)估計(jì)值為1得到的觀測(cè)值；只是這個(gè)由我估計(jì)得到的 第24行

42、的y也是觀測(cè)值 不過(guò)是由加了噪聲的真實(shí)值得到的 %相關(guān)矩陣 e_A = linear + 25 * (1-e_x_estimate2)/(1+e_x_estimate2)2);%傳遞矩陣 e_H = e_x_estimate_1/10; %觀測(cè)矩陣 %估計(jì)的誤差 e_p_estimate = e_A * e_P * e_A' + Q; %擴(kuò)展卡爾曼增益 e_K = e_p_estimate * e_H'/(e_H * e_p_estimate * e_H' + R); %進(jìn)行估計(jì)值的更新 第二階段 e_x_estimate_2 = e_x_estimate_1 + e_

43、K * (y - e_y_estimate); %更新后的估計(jì)值的誤差 e_p_estimate_update = e_p_estimate - e_K * e_H * e_p_estimate; %進(jìn)入下一次迭代的參數(shù)變化 e_P = e_p_estimate_update; e_x_estimate = e_x_estimate_2; %采樣點(diǎn)的選取 存在x(i) u_x_par = u_x_estimate; for i = 2 : (u_symmetry_number+1) u_x_par(i,:) = u_x_estimate + sqrt(u_symmetry_number+u_k

44、) * u_P); end for i = (u_symmetry_number+2) : u_total_number u_x_par(i,:) = u_x_estimate - sqrt(u_symmetry_number+u_k) * u_P); end %計(jì)算權(quán)值 u_w_1 = u_k/(u_symmetry_number+u_k); u_w_N1 = 1/(2 * (u_symmetry_number+u_k); %把這些粒子通過(guò)傳遞方程 得到下一個(gè)狀態(tài) for i = 1: u_total_number u_x_par(i) = 0.5 * u_x_par(i) + 25 * u

45、_x_par(i)/(1+u_x_par(i)2) + 8 * cos(1.2*(k-1); end %傳遞后的均值和方差 u_x_next = u_w_1 * u_x_par(1); for i = 2 : u_total_number u_x_next = u_x_next + u_w_N1 * u_x_par(i); end u_p_next = Q + u_w_1 * (u_x_par(1)-u_x_next) * (u_x_par(1)-u_x_next)' for i = 2 : u_total_number u_p_next = u_p_next + u_w_N1 * (

46、u_x_par(i)-u_x_next) * (u_x_par(i)-u_x_next)' end % %對(duì)傳遞后的均值和方差進(jìn)行采樣 產(chǎn)生粒子 存在y(i)% u_y_2obser(1) = u_x_next;% for i = 2 : (u_symmetry_number+1)% u_y_2obser(i,:) = u_x_next + sqrt(u_symmetry_number+k) * u_p_next);% end% for i = (u_symmetry_number + 2) : u_total_number% u_y_2obser(i,:) = u_x_next -

47、sqrt(u_symmetry_number+u_k) * u_p_next);% end%另外存在y_2obser(i) 中； for i = 1 :u_total_number u_y_2obser(i,:) = u_x_par(i); end %通過(guò)觀測(cè)方程 得到一系列的粒子 for i = 1: u_total_number u_y_2obser(i) = u_y_2obser(i)2/20; end %通過(guò)觀測(cè)方程后的均值 y_obse u_y_obse = u_w_1 * u_y_2obser(1); for i = 2 : u_total_number u_y_obse = u_

48、y_obse + u_w_N1 * u_y_2obser(i); end %Pzz測(cè)量方差矩陣 u_pzz = R + u_w_1 * (u_y_2obser(1)-u_y_obse)*(u_y_2obser(1)-u_y_obse)' for i = 2 : u_total_number u_pzz = u_pzz + u_w_N1 * (u_y_2obser(i) - u_y_obse)*(u_y_2obser(i) - u_y_obse)' end %Pxz狀態(tài)向量與測(cè)量值的協(xié)方差矩陣 u_pxz = u_w_1 * (u_x_par(1) - u_x_next)* (u

49、_y_2obser(1)-u_y_obse)' for i = 2 : u_total_number u_pxz = u_pxz + u_w_N1 * (u_x_par(i) - u_x_next) * (u_y_2obser(i)- u_y_obse)' end %卡爾曼增益 u_K = u_pxz/u_pzz; %估計(jì)量的更新 u_x_next_optimal = u_x_next + u_K * (y - u_y_obse);%第一步的估計(jì)值 + 修正值； u_x_estimate = u_x_next_optimal; %方差的更新 u_p_next_update =

50、u_p_next - u_K * u_pzz * u_K' u_P = u_p_next_update; %進(jìn)行畫(huà)圖程序 x_array = x_array,x; e_x_estimate_array = e_x_estimate_array,e_x_estimate; u_x_estimate_array = u_x_estimate_array,u_x_estimate; e_error(k,:) = abs(x_array(k)-e_x_estimate_array(k); u_error(k,:) = abs(x_array(k)-u_x_estimate_array(k); end t = 0 : tf; f