張 明，陸東亮，徐詩露，夏若平，何順帆

（1.武漢紡織大學電子與電氣工程學院，湖北武漢 430200；2.中南民族大學計算機科學學院，湖北武漢 430074）

0 引言

諧波治理一直都是電力系統(tǒng)中的關(guān)鍵問題。隨著電網(wǎng)中大量可再生能源的接入，加上非線性負載的廣泛使用，電網(wǎng)結(jié)構(gòu)異常復(fù)雜并由此產(chǎn)生了大量諧波[1-2]。因此，如何有效治理諧波成為了供用電雙方重點關(guān)注的問題，需要不斷提升電力諧波信號的檢測精度[3]。

對電力系統(tǒng)諧波進行快速準確的檢測，目前人們已經(jīng)提出了很多檢測諧波的方法[4]，從諧波檢測的狀態(tài)分類，主要可分為穩(wěn)態(tài)和動態(tài)諧波檢測，其中穩(wěn)態(tài)諧波比較常見也比較容易進行治理，而動態(tài)諧波僅出現(xiàn)于某種運行工況之下，檢測的難度較高[5-6]，故動態(tài)諧波的檢測算法是本文研究的重點。各種諧波檢測的文獻表明卡爾曼濾波（Kalman Filter，KF）及其擴展算法都可以檢測動態(tài)諧波，文獻[7]首次將局部集成變換的卡爾曼濾波器算法應(yīng)用于電力系統(tǒng)諧波檢測，能夠提高動態(tài)諧波的檢測精度，但是其要求處理的信號具有線性關(guān)系。文獻[8-9]提出了一種擴展卡爾曼濾波（Extended Kalman Filter，EKF）算法，這是一種局部方法，能夠為非線性系統(tǒng)提供次優(yōu)解，但是這種方法可能會因為模型的非線性化而導(dǎo)致結(jié)果出現(xiàn)偏差。文獻[10]在EKF 的基礎(chǔ)上引入了無跡卡爾曼濾波（Unscented Kalman Filter，UKF）算法，它減少了預(yù)測狀態(tài)的線性負擔，能夠?qū)⒕忍嵘炼A等級，但是在處理二階函數(shù)時仍然不能給出正確的二階時刻。文獻[11]提出了一種容積卡爾曼（Cubature Kalman Filter，CKF）算法來提高諧波的檢測精度，其利用了三階球面徑向規(guī)則，能夠解決高維狀態(tài)空間問題，且由于其采用的非線性模型，能夠克服線性化問題，其但是建立在噪聲矩陣不變的情況下。在動態(tài)諧波檢測中，不僅要考慮諧波信號的非線性化問題，還需要考慮高維問題以及因為系統(tǒng)噪聲引起的檢測誤差，因此，本文在CKF的基礎(chǔ)上通過引入噪聲估值的遺忘因子，提出了一種自適應(yīng)容積卡爾曼濾波（Adaptive Cubature Kalman Filter，ACKF）算法，仿真結(jié)果表明本文提出的算法在動態(tài)諧波檢測中具有更高的精度。

1 卡爾曼濾波算法原理

卡爾曼濾波是一種基于時域的濾波方法，將某一時刻測量的數(shù)據(jù)作為量測量來預(yù)測下一時刻的狀態(tài)量。設(shè)k時刻一組隨機的離散線性方程為[12-13]：

式中：Xk為系統(tǒng)的狀態(tài)向量；Ak-1為k-1 時刻的狀態(tài)轉(zhuǎn)移矩陣；Wk-1為k-1 時刻的系統(tǒng)噪聲向量；Zk為系統(tǒng)的量測向量；Hk為系統(tǒng)的量測矩陣；Vk為系統(tǒng)的量測噪聲向量。

假設(shè)系統(tǒng)中出現(xiàn)的噪聲都是高斯白噪聲，則其均值滿足為0，系統(tǒng)噪聲方差矩陣為Q，量測噪聲方差矩陣為R。

由式（1）可知，當前狀態(tài)的預(yù)測X′k-1是根據(jù)當前時刻測量的數(shù)據(jù)得來的，可以寫成關(guān)系式（2），而下一時刻的預(yù)測又是在當前預(yù)測的基礎(chǔ)上做出的，所以可以得到式（3）。

因為系統(tǒng)噪聲只對狀態(tài)向量產(chǎn)生影響，所以忽略E[Wk-1/Z1Z2…Zk-1]，E表示單位向量，由式（2）和式（3）就可以得到預(yù)測狀態(tài)向量：

而在實際應(yīng)用中，預(yù)測的狀態(tài)量與狀態(tài)量的真實值總存在一定的誤差，由式（1）可知，狀態(tài)量產(chǎn)生誤差也將會影響量測值Zk，所以就有如下推導(dǎo)：

式中：ΔX′k為預(yù)測狀態(tài)量與真實狀態(tài)量之間的誤差。

式中：Z′k為預(yù)測狀態(tài)向量出現(xiàn)誤差下的測量值；B為Zk與Z′k之間的殘差值。

由式（6）可以看出，通過處理殘差可以達到降低預(yù)測誤差的目的，引入濾波增益Kk，則可以得到修正后的狀態(tài)向量為：

在預(yù)測的過程中求得合適的濾波增益系數(shù)是本算法的關(guān)鍵。

卡爾曼濾波是一種最小均方誤差估計，在修正后的狀態(tài)向量X″k的基礎(chǔ)上繼續(xù)求得差值記為ΔX″k=Xk-X″k。

則結(jié)合式（1）、式（5）和式（7）可得：

于是有：

因為量測噪聲Vk不會對狀態(tài)向量產(chǎn)生影響，所以二者滿足不相關(guān)，就有[14]：

式中：Pk為系統(tǒng)的偏差矩陣，使其達到最小值則可以求得Kk的值：

然后將式（13）代入式（11）中可得系統(tǒng)的偏差矩陣為：

最后由狀態(tài)向量誤差推導(dǎo)出預(yù)測偏差矩陣為：

式中：Qk-1為k-1時刻系統(tǒng)噪聲的方差矩陣。

傳統(tǒng)的卡爾曼濾波算法能夠預(yù)測出動態(tài)諧波信號的相關(guān)信息，但它是建立在一種具有線性關(guān)系的基礎(chǔ)上，具有一定的局限性，需要進一步的分析研究。

2 基于自適應(yīng)容積卡爾曼濾波算法

傳統(tǒng)的卡爾曼濾波只能對線性關(guān)系的狀態(tài)方程進行處理，在實際中是無法運用的，因此需要在傳統(tǒng)方法的基礎(chǔ)上進行改進，針對此問題，提出了一種基于自適應(yīng)的容積卡爾曼濾波算法[15]。

2.1 容積卡爾曼濾波算法

對于高斯分布下的非線性濾波問題，實際上就是求解狀態(tài)量的后驗概率密度函數(shù)乘積的積分值，CKF 可以看作是一種利用三階球面徑向規(guī)則生成一組等權(quán)值的容積點來求解后經(jīng)驗期望的積分[16]。傳統(tǒng)的卡爾曼濾波的一般積分形式可以表示為：

式中：G(f)為待求的積分；Rn為n維的積分域；f(x)為非線性積分；x為狀態(tài)向量。

在容積卡爾曼濾波算法中需要將一般積分形式變換為球面徑向積分形式和三階球面徑向積分形式，經(jīng)過對式（16）的轉(zhuǎn)化，球面徑向積分形式可以表示為：

由式（17）可以發(fā)現(xiàn)，該變換是將狀態(tài)向量x的積分分解為半徑r與方向向量y相關(guān)的二重積分，式（17）中x=ry,r∈[0,+∞)，令yTy=1，因此可得xTx=r2，可以看出變換后積分中的球面單位為Un={y∈Rn|yTy=1}，σ(·)表示積分域的積分微元。

對式（17）進一步轉(zhuǎn)化可得：

式中：S(r)為球面積分，S(r)=∫Un f(ry)dσ(y)

在完成了球面徑向積分形式的轉(zhuǎn)化后，對式（18）和S(r)使用三階球面準則，可得mr×ms的球面徑向容積準則為：

式中：mr=1，ms=2n表示容積點的點數(shù)；{ai,ri} 為計算徑向積分容積點集合；{bj,yj} 為計算球面積分容積點集合。

為了更好地說明容積卡爾曼濾波算法在非線性動態(tài)諧波檢測上的應(yīng)用，設(shè)一個隨機的非線性離散方程為：

式中：xk為k時刻的狀態(tài)量；zk為k時刻的量測量；wk，vk分別為系統(tǒng)噪聲和量測噪聲，假設(shè)其為高斯白噪聲，則其方差也滿足高斯白噪聲，即wk～(0,Qk),vk～(0,Rk)，R，Q之間互不相關(guān)，且其均值為0。

容積卡爾曼濾波算法的實現(xiàn)可以分為初始化、狀態(tài)預(yù)測和狀態(tài)校正3 個步驟[17-19]。

2.1.1 初始化

在利用CKF 進行動態(tài)諧波狀態(tài)估計時需要給出前一時刻的量測量，本文給定k時刻的狀態(tài)量和估計偏差矩陣分別為xk和Pk|k。

當k=0 時刻時，假設(shè)x0～(xˉ0,P0),Q0,R0，則得濾波器的最優(yōu)初始化為：

2.1.2 狀態(tài)預(yù)測

狀態(tài)預(yù)測就是根據(jù)球面徑向規(guī)則生成的容積點集結(jié)合狀態(tài)方程來估計下一時刻的狀態(tài)量，并能夠計算狀態(tài)量的偏差矩陣，據(jù)此可以得到關(guān)于容積點的表達式如下：

式中：ξj為容積點，，其中[L]為n維空間的點集矩陣如式（23）所示；[L]j為第j個容積點；j取為權(quán)值；m為容積點的個數(shù)，在三階球面徑向規(guī)則下滿足m=2n。

根據(jù)預(yù)測的本質(zhì)，通過k時刻的估計偏差矩陣Pk|k計算諧波狀態(tài)的等權(quán)值容積點。

在得到了諧波的等權(quán)值容積點后，結(jié)合離散方程就可以估計出下一時刻系統(tǒng)的狀態(tài)向量，其表達式如下：

由容積點的非線性傳播和進一步的轉(zhuǎn)化可以計算得到諧波在k+1 時刻的估計偏差矩陣：

2.1.3 狀態(tài)校正

1）結(jié)合式（27）可計算諧波狀態(tài)向量測量的等權(quán)值容積點：

將其代入離散方程式（20）中就可以得到測量的量測向量。

再由式（26）就可以估計預(yù)測的諧波狀態(tài)測量值：

2）將估計得到的和實際測量得到的諧波狀態(tài)量對比，可以得到誤差系數(shù)，求得修正后的諧波狀態(tài)向量和諧波估計偏差矩陣。

則可得：

式中：Pvv,k+1為更新后的量測量估計偏差矩陣；Pvz,k+1|k為估計的交叉協(xié)方差矩陣。

最后由更新后的式（31）和式（32）計算得到濾波增益系數(shù)Kk+1，就可以得到k+1 時刻估計的諧波狀態(tài)向量均值和偏差矩陣Pk+1。

2.2 基于漸消記憶指數(shù)加權(quán)法的噪聲估值算法

容積卡爾曼濾波諧波動態(tài)估計可以得到很好的估計效果，但其是建立在電力系統(tǒng)中系統(tǒng)噪聲和量測噪聲都保持不變的前提下，而實際中的協(xié)方差矩陣Q和R都是在不斷變化的，這也就使得容積卡爾曼濾波只能停留在理想狀況下，不能捕捉到估計噪聲協(xié)方差所需的整體噪聲動態(tài)。因此，在容積卡爾曼濾波的基礎(chǔ)上加入了基于漸消記憶指數(shù)加權(quán)法的噪聲估值算法來得到一種自適應(yīng)容積卡爾曼濾波（ACKF）算法，以自適應(yīng)系統(tǒng)不斷變化的噪聲[20-22]。

基于漸消記憶指數(shù)加權(quán)法引入了噪聲估值遺忘因子α，這個因子決定了對過去的噪聲協(xié)方差估計和當前的噪聲瞬時估計的權(quán)重大小，并且可以提供噪聲協(xié)方差的過去和現(xiàn)在估計的加權(quán)平均值。而當遺忘因子為所有過去的估計提供恒定的權(quán)重時，可能無法跟蹤系統(tǒng)噪聲指標的變化[23]。因此，引入了衰落記憶的概念，通過不斷遺忘過去的陳舊數(shù)據(jù)而加大對現(xiàn)在估計生成的噪聲數(shù)據(jù)，以指數(shù)衰減的形式加權(quán)過去的估計，通過式（34）確定加權(quán)系數(shù)，加權(quán)系數(shù)由遺忘因子用于確定每個延遲的權(quán)重估計，過去的權(quán)重往往呈指數(shù)衰減從而產(chǎn)生衰減記憶。在動態(tài)諧波狀態(tài)估計中，由于新生成的噪聲數(shù)據(jù)與下一時刻的噪聲數(shù)據(jù)具有相似性，因此，漸消記憶指數(shù)加權(quán)法就是通過遺忘因子確定前后時刻的加權(quán)系數(shù)，由加權(quán)系數(shù)再結(jié)合前面的離散方程就可以得到新的系統(tǒng)噪聲和量測噪聲的協(xié)方差矩陣。

式中：dk為加權(quán)系數(shù)；α取0.85；Q′k+1為k+1 時刻生成的系統(tǒng)噪聲的協(xié)方差矩陣；R′k+1為k+1 時刻生成的量測噪聲的協(xié)方差矩陣，其中εk=Zk-h(xk)，表示殘差值。

在原始的算法中加入基于漸消記憶指數(shù)加權(quán)法的噪聲估值器，得到自適應(yīng)的容積卡爾曼濾波算法，其步驟為：

1）采用漸消記憶指數(shù)加權(quán)法生成系統(tǒng)噪聲和量測噪聲的時變協(xié)方差矩陣分別表示為Q′k+1和R′k+1；

2）用修正后的Q′k+1替換式（27）中的Qk+1得到k+1 時刻諧波狀態(tài)量的估計偏差矩陣，然后由式（28）和式（30）分別計算諧波量測量的等權(quán)容積點和諧波量測量，進而就可以求得在時變噪聲下的系統(tǒng)狀態(tài)向量和總的誤差偏差矩陣。

3 動態(tài)諧波檢測的仿真分析

3.1 仿真分析

為了驗證本文ACKF 算法用于動態(tài)諧波檢測的有效性，選擇KF，CKF 和ACKF 3 種算法，通過仿真分析比較這3 種算法進行動態(tài)諧波檢測的性能。

為了評估算法性能，引入評價指標幅值均方根誤差Ra和相角均方根誤差Rp分別衡量檢測諧波幅值和相角的準確度，其計算式分別為：

首先將ACKF 算法與CKF 算法進行比較驗證。仿真分析采用如表1 所示的動態(tài)諧波信號頻譜，其他仿真參數(shù)設(shè)置為采樣頻率4 000 Hz，采樣時間0.2 s，數(shù)據(jù)共800 個點，其中在0.1 s 時幅值增加1.5 倍，同時加入高斯噪聲，信噪比為50 dB。

表1 動態(tài)諧波信號頻譜表Table 1 Dynamic harmonic signal spectrum

CKF 算法用于動態(tài)諧波檢測的仿真結(jié)果如圖1～圖6 所示。

圖1 原始諧波信號與CKF算法跟蹤信號Fig.1 Original harmonic signal and tracking signal with CKF algorithm

圖1 為對原始動態(tài)諧波信號的動態(tài)跟蹤效果，圖2—圖6 展示了該算法實時估計出的各次諧波的幅值和相位。由圖1 可以看出CKF 算法可以追蹤到動態(tài)諧波信號，但在初始時刻出現(xiàn)了一些偏差。由圖2—圖6 中的波形可以得出CKF 算法對基波和各次諧波的幅值估計在初始位置存在誤差，需要進一步修正。

圖2 CKF算法估計的基波幅值與相角Fig.2 Fundamental amplitude and phase angle estimated by CKF algorithm

圖3 CKF算法估計的5次諧波幅值與相角Fig.3 Fifth harmonic amplitude and phase angle estimated by CKF algorithm

圖4 CKF算法估計的7次諧波幅值與相角Fig.4 Seventh harmonic amplitude and phase angle estimated by CKF algorithm

圖5 CKF算法估計的11次諧波幅值與相角Fig.5 11th harmonic amplitude and phase angle estimated by CKF algorithm

圖6 CKF算法估計的13次諧波幅值與相角Fig.6 13th harmonic amplitude and phase angle estimated by CKF algorithm

ACKF 算法用于動態(tài)諧波檢測的仿真結(jié)果如圖7—圖12 所示。其中圖7 為對原始動態(tài)諧波信號的動態(tài)跟蹤效果，圖8—圖12 展示了算法實時估計出的各次諧波的幅值和相位。

圖7 原始諧波信號與ACKF算法跟蹤信號Fig.7 Original harmonic signal and tracking signal with ACKF algorithm

圖8 ACKF算法估計的基波幅值與相角Fig.8 Fundamental amplitude and phase angle estimated by ACKF algorithm

圖9 ACKF算法估計的5次諧波幅值與相角Fig.9 Fifth harmonic amplitude and phase angle estimated by ACKF algorithm

圖10 ACKF算法估計的7次諧波幅值與相角Fig.10 Seventh harmonic amplitude and phase angle estimated by ACKF algorithm

圖11 ACKF算法估計的11次諧波幅值與相角Fig.11 11th harmonic amplitude and phase angle estimated by ACKF algorithm

圖12 ACKF算法估計的13次諧波幅值與相角Fig.12 13th harmonic amplitude and phase angle estimated by ACKF algorithm

由圖7 中的跟蹤波形可以看出，ACKF 算法的準確度更高，在初始時刻的波形抖動很小。由圖8—圖12 可以看出ACKF 算法在幅值和相位上能準確估計出真實值。

最后將KF 算法同樣應(yīng)用于對基波、5 次諧波、7次諧波、11 次諧波和13 次諧波的動態(tài)諧波信號檢測中，比較結(jié)果見表2。從表2 中可以看出，在有噪聲干擾下，ACKF 算法的動態(tài)諧波檢測精度在總體上要高于CKF 算法與KF 算法，能夠很好地跟蹤動態(tài)諧波信號，同時也能很精確地估計出各階次諧波的幅值和相位。

表2 3種算法的Ra，Rp值比較Table 2 Comparison of Ra and Rp among three algorithms

3.2 在有源濾波器諧波檢測中的應(yīng)用

有源濾波器（Active Power Filter，APF）[24-25]主要用于抑制電力系統(tǒng)的動態(tài)諧波，其高效工作的基礎(chǔ)是實時、準確地檢測諧波，傳統(tǒng)的諧波檢測方法如傅立葉變換（Fast Fourier Transform，F(xiàn)FT），雖然有廣泛應(yīng)用，但一般會有1 個周期以上的時間延遲，同時存在檢測精度不高等問題。這里嘗試將ACKF 算法應(yīng)用到APF 中，其系統(tǒng)原理框圖如圖13 所示。

圖13 有源濾波器系統(tǒng)結(jié)構(gòu)框圖Fig.13 Structure diagram of active filter system

基于Simulink 搭建400 V 有源濾波器仿真系統(tǒng)。采樣頻率設(shè)置為4 000 Hz，仿真時間為0.2 s，采集電流數(shù)據(jù)共800 個點，其中在0.1 s 時將負荷增加1.2 倍，同時加入高斯噪聲，信噪比為50 dB，其仿真結(jié)果如圖14—圖17 所示。

圖14 ACKF算法估計的5次諧波幅值與相角Fig.14 Fifth harmonic amplitude and phase angle estimated by ACKF algorithm in active filter simulation

圖15 ACKF算法估計的7 次諧波幅值與相角Fig.15 Seventh harmonic amplitude and phase angle estimated by ACKF algorithm in active filter simulation

圖16 ACKF算法估計的11次諧波幅值與相角Fig.16 11th harmonic amplitude and phase angle estimated by ACKF algorithm in active filter simulation

圖17 ACKF算法估計的13次諧波幅值與相角Fig.17 13th harmonic amplitude and phase angle estimated by ACKF algorithm in active filter simulation

圖14—圖17 展示了ACKF 算法實時估計出的5 次，7 次，11 次，13 次諧波的幅值和相位，由仿真結(jié)果可以看出，在有噪聲干擾下，ACKF 算法應(yīng)用于有源濾波器諧波電流檢測，實時性更好、精度更高。

4 結(jié)語

以電網(wǎng)某些運行工況下出現(xiàn)動態(tài)諧波信號為切入點，介紹了卡爾曼濾波動態(tài)諧波檢測算法的基本原理，針對傳統(tǒng)卡爾曼濾波只能處理線性關(guān)系下的動態(tài)諧波信號等問題，采用了一種容積卡爾曼濾波算法，在描述了容積卡爾曼濾波的原理后詳細推導(dǎo)了其實現(xiàn)的3 個步驟，然后又針對容積卡爾曼濾波算法存在的缺點引入了基于漸消記憶指數(shù)加權(quán)法的噪聲估值算法，將兩者結(jié)合起來構(gòu)成了自適應(yīng)容積卡爾曼濾波（ACKF）算法，最后通過仿真驗證了本文提出的ACKF 算法在動態(tài)諧波檢測中的有效性。