[發(fā)明專利]一種基于改進魯棒容積卡爾曼濾波的數(shù)據(jù)清洗方法

有效

• 申請?zhí)枺?/td>: 201910478867.0
• 申請日: 2019-06-04
• 公開（公告）號: CN110334087B
• 公開（公告）日: 2020-12-11
• 發(fā)明（設計）人: 孫輝;胡姝博;孫越峰;高正男;彭飛翔;周瑋
• 申請（專利權）人: 大連理工大學
• 主分類號: G06F16/215
• 分類號: G06F16/215
• 代理公司: 大連理工大學專利中心 21200
• 代理人: 李曉亮;潘迅

1.一種基于改進魯棒容積卡爾曼濾波的數(shù)據(jù)清洗方法，其特征在于，包括以下步驟：

步驟1：在某一時刻時，通過電力系統(tǒng)量測裝置采集得到原始數(shù)據(jù)組成的向量為z，設定需要數(shù)據(jù)清洗的時刻個數(shù)為N，得到需要清洗的總數(shù)據(jù)集合為{z₁，z₂，z₃，…，z_N}；以k代表時刻，利用k時刻的數(shù)據(jù)對k+1時刻進行預測，進而對k+1時刻的數(shù)據(jù)z_k+1進行清洗；

步驟2：建立電力系統(tǒng)的動態(tài)模型，得到電力系統(tǒng)動態(tài)方程和量測方程；

步驟3：設定起始時刻k＝1，則k+1＝2，表示利用第1個時刻的前期準備數(shù)據(jù)，對第2個時刻的數(shù)據(jù)z₂進行清洗；

步驟4：建立容積卡爾曼濾波CKF基本方程，并基于CKF基本方程，得到數(shù)據(jù)清洗要用到的前期準備數(shù)據(jù)；首先計算方程(3)～(7)得到電力系統(tǒng)的狀態(tài)量在k+1時刻的預測值再計算方程(8)～(12)得到電力系統(tǒng)的量測數(shù)據(jù)在k+1時刻的預測值然后計算方程(13)得到電力系統(tǒng)量測數(shù)據(jù)的預測誤差方差陣P_vv,k+1；最后計算方程(14)，將k+1時刻量測裝置采集的量測數(shù)據(jù)z_k+1減去預測值得到量測新息向量e_k+1；

所述的CKF主要包括預測和濾波兩個基本過程：

1)預測

S_k|k＝chol(P_k|k) (3)

式中：chol(·)表示Cholesky分解；P_k|k為k時刻狀態(tài)估計誤差方差陣；S_k|k為對P_k|k進行Cholesky分解得到的上三角矩陣；n為電力系統(tǒng)節(jié)點個數(shù)；I為n階單位矩陣；ξ為根據(jù)Cubature準則生成的等權值矩陣；ξ_i表示矩陣ξ的第i列；為k時刻狀態(tài)估計值；X_i,k|k為的Cubature點；為X_i,k|k的預測值；為狀態(tài)預測值；Q_k為k時刻的過程噪聲方差陣；P_k+1|k為狀態(tài)預測誤差方差陣；f(·)是電力系統(tǒng)的狀態(tài)轉移函數(shù)；

2)濾波

S_k+1|k＝chol(P_k+1|k) (9)

Z_i,k+1|k＝h(X_i,k+1|k) (11)

P_zz,k+1＝P_vv,k+1+R_k+1 (15)

式中：S_k+1|k為對P_k+1|k進行Cholesky分解得到的上三角矩陣；X_i,k+1|k為的Cubature點；Z_i,k+1|k為量測預測值的Cubature點；為量測預測值；P_vv,k+1為量測預測誤差方差陣；e_k+1為新息向量；R_k+1是量測噪聲方差陣；P_zz,k+1為量測總方差；P_xz,k+1為狀態(tài)量與量測量的協(xié)方差陣；W_k+1為濾波系數(shù)；為k+1時刻的電力系統(tǒng)狀態(tài)估計值；P_k+1|k+1為k+1時刻狀態(tài)估計誤差方差陣；h(·)是電力系統(tǒng)的量測函數(shù)；

步驟5：對k+1時刻的電力系統(tǒng)量測裝置采集到的原始量測數(shù)據(jù)z_k+1進行不良數(shù)據(jù)辨識、修正過程；首先計算方程(20)～(27)，得到數(shù)據(jù)修正因子對角陣γ’_k+1；若γ’_k+1的對角元素等于1，表示該數(shù)據(jù)是正確數(shù)據(jù)，該數(shù)據(jù)對應的量測裝置運行正常；若γ’_k+1對角元素大于1，表示該數(shù)據(jù)是不良數(shù)據(jù)，量測裝置數(shù)據(jù)傳輸異常或者量測裝置失效；將γ’_k+1代入方程(28)，計算得到修正后的電力系統(tǒng)量測總方差陣P_zz,k+1；

所述的不良數(shù)據(jù)辨識、修正過程如下：

5.1)在保證k+1時刻前的數(shù)據(jù)已不含不良數(shù)據(jù)的前提下，選取時間序列滑動窗口長度為T，將k+2-T時刻到k+1時刻的新息向量構成時間序列[e_k-T+2,e_k-T+3,...,e_k+1]，計算窗口內的平均新息協(xié)方差：

式中：P_e,k+1為平均新息協(xié)方差；T為時間序列滑動窗口長度；t代表時間序列的第t個向量，0≤t≤T，且t為整數(shù)；e_k+2-t為k+2-t時刻的新息向量；

當量測量突變時，由于e_k+1增大，必然使得平均新息協(xié)方差大于式(14)的觀測噪聲，即：

P_e,k+1＞P_vv,k+1+R_k+1 (21)

式中：P_vv,k+1為量測預測誤差方差陣；R_k+1是量測噪聲方差陣；

構造修正因子γ_k+1實現(xiàn)突變數(shù)據(jù)的辨識；γ_k+1可對R_k+1進行調整，其取值原則是使得上式左右兩端相等，求解得：

則當可疑的突變量測量存在時，γ_k+1中的相應對角元素必將大于1；因此，取γ_k+1中大于1的對角元素對應的量測量，形成可疑數(shù)據(jù)集合D₁；D₁中的可疑數(shù)據(jù)即為突變數(shù)據(jù)，實際包含正確的波動型數(shù)據(jù)和錯誤的不良數(shù)據(jù)；

5.2)選取時間序列滑窗長度為T，將k+2-T時刻到k+1時刻的量測向量構成量測量時間序列[z_k+2-T,z_k+3-T,...,z_k+1]，計算皮爾遜相關系數(shù)形成m階相關矩陣，其對角元素為1，非對角元素計算如下：

其中：

式中：z_k+2-t為k+2-t時刻量測裝置采集到的量測數(shù)據(jù)；為量測數(shù)據(jù)的平均值；r_ij為各量測量間的相關系數(shù)；

r_ij代表第i個量測量與第j個量測量的相關性強弱；

計算第i個量測量與其他所有量測量的平均相關系數(shù)

式中：m為量測數(shù)據(jù)個數(shù)，亦為相關矩陣階數(shù)；

為辨識可疑的不良數(shù)據(jù)，采用“末位淘汰”方法，選取與其他量相關性最弱的量測量形成可疑數(shù)據(jù)集合D₂；

采用定D₂元素個數(shù)的方法，其元素個數(shù)d₂的選取參考D₁的元素個數(shù)d₁，則可取d₂＝αd₁，α為可疑集長度參考倍數(shù)，可根據(jù)實際系統(tǒng)選?。?/p>

5.3)在得到上述可疑數(shù)據(jù)集合D₁和D₂之后，取二者交集得到不良數(shù)據(jù)集B＝D₁∩D₂，則B中的量測量才是真正的不良數(shù)據(jù)，其同時滿足量測量發(fā)生突變和與整體量測數(shù)據(jù)相關性最弱兩個條件；而B在D₁中的補集C_D1B實質為因新能源出力波動而產生的波動型數(shù)據(jù)，是正確的量測數(shù)據(jù)；

利用γ_k+1對不良數(shù)據(jù)進行修正；由于γ_k+1為非對角陣，為避免矩陣求逆運算發(fā)生奇異定義修正因子對角陣γ’_k+1，其對角元素取值為：

對式(15)進行如下修改以完成不良數(shù)據(jù)修正：

P_zz,k+1＝P_vv,k+1+γ'_k+1R_k+1 (28)

步驟6：對k+1時刻的原始量測數(shù)據(jù)z_k+1進行數(shù)據(jù)清洗；首先計算方程(16)～(18)，完成k+1時刻的電力系統(tǒng)狀態(tài)估計過程，得到k+1時刻的狀態(tài)估計值然后將代入方程(2)的量測函數(shù)，得到k+1時刻電力系統(tǒng)量測數(shù)據(jù)的清洗值

步驟7：計算方程(19)，得到k+1時刻狀態(tài)估計誤差方差陣P_k+1|k+1，以用作下一時刻數(shù)據(jù)清洗的前期準備數(shù)據(jù)；判斷k的大小，若k≤N-1，則使k＝k+1，循環(huán)步驟4至步驟6，開始對下一時刻的電力系統(tǒng)量測數(shù)據(jù)進行清洗；若kN-1，則數(shù)據(jù)清洗過程結束，得到數(shù)據(jù)清洗值的集合