基于貝葉斯估計的故障電流過零點預測

劉云飛，周啟文，方太勛，2

（1.南京南瑞繼保電氣有限公司，南京 211102；2.南瑞集團有限公司（國網(wǎng)電力科學研究院有限公司），南京 211106）

0 引言

隨著新型電力系統(tǒng)的發(fā)展，電網(wǎng)故障電流水平持續(xù)攀升，導致斷路器開斷能力不足，嚴重影響電網(wǎng)安全穩(wěn)定運行，成為制約電力建設和地區(qū)經濟發(fā)展的瓶頸［1-2］。傳統(tǒng)故障電流限制措施從優(yōu)化電網(wǎng)結構、調整運行方式、安裝高阻抗設備等方面考慮，對電網(wǎng)運行的靈活性與經濟性都會帶來不同程度的負面影響。近年來，基于智能快速斷路器的故障電流抑制技術已發(fā)展成為一種新的針對故障電流超標問題的治理措施［3-5］，通過在故障后快速重構電網(wǎng)拓撲以增大系統(tǒng)阻抗，達到降低故障電流的目的。智能快速斷路器一般采用斥力機構驅動真空斷路器實現(xiàn)快速性的要求［6-8］，但目前商用真空斷路器的額定開斷電流最大僅有50 kA，同樣存在開斷能力不足的問題，在某些應用場合難以滿足應用需求。相控開斷技術通過預測故障電流過零點以控制斷路器的動作時刻，可有效控制燃弧時間，減小觸頭電侵蝕，提高斷路器開斷能力，同時對提升斷路器的可靠性和壽命具有十分重要的意義。

故障電流過零點預測算法的關鍵在于電流幅值、相位和衰減直流分量等狀態(tài)變量的快速估計，其中預測時間、預測誤差和算法復雜度是評價算法優(yōu)劣的重要指標。目前電力系統(tǒng)廣泛采用全波傅式算法或者半波傅式算法提取故障電流特征量，但在處理含衰減直流分量的電流波形時誤差較大。文獻［9］提出了基于安全點算法的相控開斷方法，即忽略故障電流衰減直流分量的變化，以安全點代替實際電流過零點，同樣存在預測誤差較大的問題。文獻［10-11］提出了基于RLS（遞推最小二乘法）算法的故障電流參數(shù)提取方法，即以均方誤差為代價函數(shù)擬合故障電流波形。文獻［12-15］提出了基于自適應神經網(wǎng)絡算法、WLMS（加權最小均方差）算法、自適應Prony 算法等的故障電流參數(shù)提取算法，這些算法通常均需要在線進行矩陣運算，計算量大且需存儲大量采樣數(shù)據(jù)，在邊緣計算中對處理器性能提出很高的要求。

本文研究了基于貝葉斯框架的卡爾曼濾波和粒子濾波在故障電流過零點預測中的應用，通過對故障電流衰減時間常數(shù)的不同處理策略建立故障電流模型，提出擴展卡爾曼濾波、標準卡爾曼濾波和粒子濾波3種算法及其實現(xiàn)，最后通過算例仿真對比分析了算法的預測結果。

1 故障電流過零點預測

1.1 貝葉斯估計

貝葉斯估計將狀態(tài)估計視為一個概率推理過程，即將估計問題轉化為利用貝葉斯公式求解后驗概率密度，包括預測和更新兩個過程，預測過程利用系統(tǒng)模型預測狀態(tài)的先驗概率密度，更新過程利用最新的量測值進行修正，得到后驗概率密度。

卡爾曼濾波是一種基于一階馬爾科夫模型的貝葉斯估計，其最大優(yōu)點在于計算量小，能夠利用前一時刻的狀態(tài)和量測值遞推計算當前時刻狀態(tài)的最優(yōu)估計。假設線性系統(tǒng)動態(tài)方程表示為：

式中：Xk為n×1維的狀態(tài)向量；Zk為m×1維的量測向量；Φk/k-1、Γk/k-1、Hk是已知的系統(tǒng)結構參數(shù)，分別為n×n維的狀態(tài)一步轉移矩陣、n×l 維的系統(tǒng)噪聲分配矩陣和m×n維的量測矩陣；Wk為l×1維的過程噪聲向量；Vk為m×1 維的量測噪聲向量。兩者都是零均值的高斯白噪聲向量序列，且它們之間互不相關。

標準卡爾曼濾波過程可由以下方程進行遞推求解：

式中：λ為衰減記憶因子。當λ＞1 時，表示狀態(tài)估計過程更依賴于量測值，以解決由模型誤差帶來的估計發(fā)散問題。

粒子濾波通過非參數(shù)化的蒙特卡洛模擬方法來實現(xiàn)遞推貝葉斯估計，適用于任何能用狀態(tài)空間模型描述的非線性系統(tǒng)，精度可逼近最優(yōu)估計［16］。粒子濾波以有限個參數(shù)來近似后驗，后驗分布的樣本叫作粒子，一個粒子就是根據(jù)真實世界狀態(tài)在時刻t的一種可能假設。根據(jù)序貫重要性采樣原理，假設后驗概率為，在重要性分布中采樣獲得一組M個離散樣本點，則可以得到粒子的重要性權重：

經過歸一化得到，進而通過重采樣得到狀態(tài)估計值：

1.2 故障電流模型

電網(wǎng)故障電流主要由基波、諧波、衰減直流分量等成分組成，電流量測方程可以表示為：

式中：N為總諧波次數(shù)；In為n次諧波電流幅值；φn為n次諧波電流相位；I0為故障電流直流分量初值；τ為故障電流衰減時間常數(shù)；ω為電網(wǎng)頻率；v（t）為電流量測噪聲。

在故障期間可以認為電網(wǎng)頻率為已知量，即保持工頻不變，若選取In、φn、I0、τ分別作為未知狀態(tài)變量，并假設各狀態(tài)量在計算過程中保持不變，則狀態(tài)方程為線性方程，而量測方程為非線性方程。其中非線性項由直流分量的衰減時間常數(shù)τ引起，對于該非線性項的處理是故障電流狀態(tài)估計的難點。

采用卡爾曼濾波估計故障電流狀態(tài)變量，但卡爾曼濾波的局限性在于其只適用于線性系統(tǒng)，因此需要將非線性系統(tǒng)近似為線性系統(tǒng)，即擴展卡爾曼濾波，采用一階泰勒級數(shù)展開進行線性化處理。另一方面，根據(jù)電路原理，衰減時間常數(shù)τ主要由網(wǎng)架拓撲結構即電網(wǎng)等值電抗L和等值電阻R決定，對于某一運行方式下的區(qū)域電網(wǎng)，其故障電流衰減時間常數(shù)τ一般可認為是常數(shù)。因此可以借助數(shù)字孿生電網(wǎng)或者離線仿真等技術事先求解，此時量測方程即為線性方程，進而采用標準卡爾曼濾波即可估計故障電流狀態(tài)。同時，由于粒子濾波處理非線性問題的優(yōu)越性，采用粒子濾波進行故障電流狀態(tài)估計成為另外一種可行的方法，此時將不需要對系統(tǒng)動態(tài)方程作任何處理。文中對這3種算法分別進行分析討論。

1.3 狀態(tài)估計算法

1.3.1 擴展卡爾曼濾波算法

采用擴展卡爾曼濾波進行故障電流狀態(tài)估計，選取In、φn、I0、τ作為狀態(tài)變量，狀態(tài)轉移矩陣Φk/k-1為單位矩陣，采樣周期為Δ，根據(jù)故障電流模型和系統(tǒng)動態(tài)方程可以得到：

遞推估計出狀態(tài)變量Xk后，即可根據(jù)式（6）計算故障電流各次諧波幅值和相位、衰減直流分量，計算公式如下：

由于一階泰勒級數(shù)展開引入截斷誤差，需要對非線性項I0和τ的估計結果進行修正，即根據(jù)擬合的隨計算窗口時間變化的校正系數(shù)補償線性化帶來的估計誤差。

1.3.2 標準卡爾曼濾波算法

采用標準卡爾曼濾波進行故障電流狀態(tài)估計，選取In、φn、I0作為狀態(tài)變量，將衰減時間常數(shù)τ作為已知常數(shù)，此時狀態(tài)估計結果的準確性將直接取決于衰減時間常數(shù)τ。同樣，根據(jù)故障電流模型和系統(tǒng)動態(tài)方程可以得到線性方程如下：

卡爾曼算法涉及矩陣運算，然而，由卡爾曼狀態(tài)估計方程可知，在初始條件確定的情況下，矩陣Pk/k-1、Kk、Pk的計算僅由系統(tǒng)參數(shù)Φk/k-1、Hk、Qk、Rk決定，與量測值Zk無關。因此可以采用離線、在線相結合的計算策略，事先對卡爾曼增益Kk和量測矩陣Hk進行離線計算，計算結果作為宏參數(shù)保存在內存中，當系統(tǒng)實時運行時直接調用，不同通道的狀態(tài)計算過程可以復用參數(shù)。進一步地，由于狀態(tài)轉移矩陣Φk/k-1為單位矩陣，卡爾曼狀態(tài)估計方程的在線實時計算過程可以簡化為僅對的遞推求解。式（9）給出了電流狀態(tài)量離散計算過程。

式中：下標表示第k個計算時刻；上標表示對應向量的第n個元素；Zk為第k個時刻電流采樣值。

可以看到，電流單個狀態(tài)量的計算平均僅需兩次加法和兩次乘法，計算量遠小于其他故障特征提取算法。目前電力系統(tǒng)控制保護裝置采樣及程序執(zhí)行周期一般為833 μs（常規(guī)保護裝置）或者100 μs（故障限流控制裝置），CPU 主頻通?？梢赃_到幾百至上千兆赫，因此上述卡爾曼算法對控制保護裝置負載率不會產生明顯影響。

卡爾曼狀態(tài)估計中P0、Q、R等參數(shù)影響狀態(tài)估計性能。協(xié)方差矩陣P0為系統(tǒng)對初始狀態(tài)X0的置信度，由于通常X0無法準確獲取，P0取值應盡量大，以加快狀態(tài)估計過程的收斂速度，本文取P0=κ×diag（a1，a2，…，an），其中κ為常數(shù)。協(xié)方差矩陣Q表示模型誤差，例如模型的線性化、離散化誤差等；協(xié)方差矩陣R為量測誤差，與傳感器的特性相關。這兩個參數(shù)通過影響卡爾曼增益K的值，進而影響預測值和量測值的權重。在不知道P0、Q、R的準確先驗信息的情況下，應適當增大Q的取值，以增大對實時量測值的利用權重，進而根據(jù)準確度、動態(tài)性能等要求進行調整。

1.3.3 粒子濾波算法

由于標準粒子濾波在低信噪比條件下可能出現(xiàn)估計結果發(fā)散的情況，本文針對狀態(tài)方程為線性、量測方程為非線性的情況采用粒子濾波混合卡爾曼濾波的方法求解，即先采用標準粒子濾波對狀態(tài)變量進行預估計，由于狀態(tài)方程為線性方程，進而可以通過卡爾曼濾波進一步估計。

當觀測噪聲較大時，后驗概率存在誤差Ek，因此可以將作為量測值，構建如下的線性系統(tǒng)：

對此線性系統(tǒng)采用標準卡爾曼濾波再次進行狀態(tài)估計，得到最終估計值。

在故障電流狀態(tài)估計時，選取In、φn、I0、τ作為狀態(tài)變量，此時無需對故障電流量測方程進行線性化處理，只需要選取一組狀態(tài)量的樣本點進行遞推求解。

假設粒子數(shù)為M，狀態(tài)變量個數(shù)為2N+2，生成初始粒子、生成高斯分布隨機數(shù)、計算先驗概率和重采樣的計算量分別為A1、A2、A3、A4，則粒子濾波的時間復雜度可以表示為：

可以看到，與卡爾曼濾波相比，粒子濾波的計算量大大增加，且主要取決于粒子數(shù)。粒子濾波算法通常要求采用高性能DSP（數(shù)字信號處理器）或FPGA（現(xiàn)場可編程門陣列）來實現(xiàn)，并需要進行算法的優(yōu)化設計，對控制保護裝置的計算能力要求較高。

計算出故障電流模型參數(shù)In、φn、I0、τ后，可由式（5）中的故障電流模型直接預測未來一段時間內的電流值，然后根據(jù)二分法搜索故障電流第一、第二等過零點時刻，根據(jù)斷路器分閘時間和燃弧時間計算相控開斷所需等待時間。所述相控開斷過程如圖1所示。圖1中：T0為故障發(fā)生時刻；T1為檢測出故障、狀態(tài)估計開始時刻；T2為狀態(tài)估計結束、預測故障電流時刻；T3為發(fā)出分閘指令時刻；T4為過零開斷時刻。相應地，t1為故障檢測時間；t2為狀態(tài)估計時間；t3為等待時間；t4為斷路器分閘和燃弧時間。相控開斷技術通過過零點預測算法將開關設備的燃弧時間轉化為控制保護裝置的等待時間，因此不會增加繼電保護或者限流設備的總體動作時間，對系統(tǒng)不會產生負面影響。

圖1 故障電流過零點預測及相控時序

2 算例分析

假設故障電流量測值周期分量幅值為20 kA，相角為π/3，直流分量幅值為10 kA，衰減時間常數(shù)為45 ms，狀態(tài)估計時長為40 ms。同時由于實際應用中采樣、處理等環(huán)節(jié)引入的噪聲，本文在電流波形中加入信噪比為40 dB的高斯白噪聲，采樣頻率為10 kHz，分別采用上節(jié)所述擴展卡爾曼濾波、標準卡爾曼濾波、粒子濾波進行狀態(tài)估計，同時與常規(guī)基于最小二乘法的故障電流特征參數(shù)提取方法進行比較，其中最小二乘法的故障電流模型與標準卡爾曼濾波模型保持一致，即滿足式（8）。為減小計算量、加快狀態(tài)估計收斂速度，本文故障電流模型中的周期分量僅考慮基波。同時，粒子濾波算法中，根據(jù)初始的量測值和量測方程開環(huán)計算出狀態(tài)量的初始估計值，并在初始估計值附近按照20%計算誤差的均勻分布產生初始粒子，粒子數(shù)取200。標準卡爾曼濾波算法中，衰減時間常數(shù)按照+10%的計算誤差（即49.5 ms）進行估計。

圖2 為幾種算法關于In、φn、I0、τ這4 個狀態(tài)變量的估計結果。隨著時序數(shù)據(jù)的遞推求解，4種算法最終均收斂于真實值，其中擴展卡爾曼濾波由于衰減時間常數(shù)收斂較慢，在15 ms左右時基本得到可用的狀態(tài)估計值，此時4個狀態(tài)量的估計誤差分別為0.12%、2.15%、1.38%、0.93%，標準卡爾曼濾波和粒子濾波收斂時間分別為6 ms和10 ms 左右，對應誤差分別為0.41%、2.26%、0.34%、10%和0.24%、1.33%、2.70%、0.46%，最小二乘法收斂過程與標準卡爾曼濾波基本相同，但誤差比標準卡爾曼濾波要大，分別為0.74%、2.86%、1.44%、10%，這是由算法的特性所決定的，即最小二乘法作為卡爾曼濾波的特例，不考慮狀態(tài)轉移方程中的過程噪聲，因此對于相同的線性模型，常以卡爾曼濾波作為狀態(tài)的最優(yōu)估計。當然，也可以采用最小二乘法求解式（6）中的故障電流模型，此時計算結果除誤差較大外與擴展卡爾曼濾波算法也基本類似，此處不再贅述。

圖2 4種算法估計的狀態(tài)變量時序值

基于本文所提算法分別以3個時刻的狀態(tài)估計值代入故障電流模型方程中預測未來一段時間內的電流波形，如圖3所示。可以看到，3種算法預測的故障電流波形均與原始電流基本重合。

圖3 貝葉斯估計算法預測的故障電流

在標準卡爾曼濾波算法中，由于衰減時間常數(shù)事先計算，其計算誤差影響狀態(tài)估計結果。本文為討論衰減時間常數(shù)的計算誤差對故障電流過零點預測結果的影響，分別選取了時間常數(shù)為45 ms、60 ms、90 ms、120 ms 和150 ms 時不同時間常數(shù)計算值引起的過零點預測誤差，如圖4 所示，其中誤差取未來兩個過零點中的誤差最大值，每個點取10 次計算平均值?？梢钥吹?，總體預測誤差均在2 ms 以內，取決于計算值偏離真實值的程度，考慮50%的時間常數(shù)計算誤差時，過零點預測誤差基本在0.5 ms 以內，且隨著時間常數(shù)的增大，衰減特征更加明顯，從而過零點預測誤差也更小。

圖4 不同衰減時間常數(shù)下的過零點預測誤差

考慮5%的3次諧波和10%的5次諧波，同樣采用上述只計及基波的策略進行估計，估計結果如圖5 所示?？梢钥吹剑鳡顟B(tài)變量最終均會收斂，但收斂時間增加，標準卡爾曼濾波收斂時間約10 ms，粒子濾波雖然快速收斂，但估計誤差也有所增加，10 ms 時最大誤差達到5%。擴展卡爾曼濾波對衰減時間常數(shù)的估計值在20 ms之后才逐漸穩(wěn)定，考慮故障電流模型中的諧波分量后再次使用擴展卡爾曼濾波算法進行估計，結果如圖5中的擴展卡爾曼濾波2曲線，在20 ms基本能夠準確估計衰減時間常數(shù)在內的各個狀態(tài)?？梢钥吹?，考慮故障電流模型中的諧波分量后可以提高狀態(tài)估計精度，但也增加了算法復雜度，同時由于特征諧波次數(shù)與實際工況密切相關，需要根據(jù)實際應用折中取舍。

圖5 加入諧波后貝葉斯算法估計的狀態(tài)變量時序值

綜上所述，標準卡爾曼濾波算法由于提前計算衰減時間常數(shù)這一非線性項，對其他狀態(tài)量的在線估計速度較快；擴展卡爾曼濾波算法則適用于衰減時間常數(shù)無法提前計算且對收斂時間要求不高的場合；粒子濾波由于計算量大，在控制保護裝置中的實現(xiàn)存在挑戰(zhàn)。由于目前故障限流裝置動作時間一般要求20 ms以內，常規(guī)保護裝置出口時間一般在20 ms以上，因此以本文方法預測時間能夠滿足相關控制保護裝置動作時間的要求。

3 結語

本文研究了基于貝葉斯框架的故障電流過零點預測算法，通過對非線性項的不同處理策略提出了擴展卡爾曼濾波、標準卡爾曼濾波和粒子濾波3 種算法及其實現(xiàn)。粒子濾波雖然整體性能較優(yōu)，但存在計算量大的問題，在邊緣計算中難以實現(xiàn)。而通過事先計算衰減時間常數(shù)進而采用標準卡爾曼濾波進行狀態(tài)估計的方法可以快速預測故障電流過零點，且預測誤差相對較小，具有工程應用價值。對于含有低頻諧波的故障電流，在故障電流模型中考慮相應的特征次諧波，可以減小估計誤差，但同時也將增加狀態(tài)估計的計算量，需要在實際應用中綜合考慮。基于貝葉斯估計的故障電流過零點預測算法在控制保護裝置中的實現(xiàn)和驗證工作有待進一步開展。