CN110334087A - 一种基于改进鲁棒容积卡尔曼滤波的数据清洗方法 - Google Patents

Abstract

一种基于改进鲁棒容积卡尔曼滤波的数据清洗方法，属于数据处理技术领域。首先，整理原始量测数据，建立电力系统的动态模型。其次，建立容积卡尔曼滤波CKF基本方程，得到数据清洗要用到的前期准备数据。再次，对k+1时刻的电力系统量测装置采集到的原始量测数据z_k+1进行不良数据辨识、修正过程；对k+1时刻的原始量测数据z_k+1进行数据清洗。最后，计算得到下一时刻数据清洗的前期准备数据，并对下一时刻的电力系统量测数据进行清洗，得到数据清洗值的集合。本发明考虑量测数据之间的相关性，通过将突变数据辨识过程与数据相关性强弱程度辨识过程进行结合，对波动型数据和不良数据进行区分，并通过修正因子对不良数据完成清洗，提高数据清洗准确性。

Description

一种基于改进鲁棒容积卡尔曼滤波的数据清洗方法

技术领域

本发明属于数据处理技术领域，涉及电力系统数据清洗方法，特别涉及到一种基于改进鲁棒容积卡尔曼滤波的数据清洗方法。

背景技术

数据清洗，就是对缺失、错误数据等不良数据的辨识修正过程，即对“脏”数据的“清洗”过程，将“脏”数据补全或更正为更接近真实数据的可靠数据。由于电力系统量测误差及数据传输干扰等问题，通过量测装置最终收集得到的量测数据不可避免会存在“脏”数据。如果不对“脏”数据加以处理便直接使用，将会对后续分析带来巨大的误差。目前大多直接将超过量测标准差一定倍数的突变数据当作“脏”数据，但是随着风电等新能源并网趋势的发展，这种具有随机性、间歇性的新能源出力将导致电力系统功率产生波动，量测数据也随之成为变化不稳定的波动型数据。因此，目前电力系统中的突变数据不仅包括“脏”数据，还包含了大量量测正确的波动型数据，若不加区分的对所有突变数据进行数据清洗，可能反而会把正确数据处理成“脏”数据。

发明内容

本发明目的是在对量测数据的数据清洗过程中，提出一种基于改进鲁棒容积卡尔曼滤波的数据清洗方法，该方法考虑量测数据之间的相关性，通过将突变数据辨识过程与数据相关性强弱程度辨识过程进行结合，对波动型数据和不良数据进行区分，并通过修正因子对不良数据完成清洗，提高数据清洗准确性。

为了达到上述目的，本发明的技术方案是：

一种基于改进鲁棒容积卡尔曼滤波的数据清洗方法，包括以下步骤：

步骤1：整理原始量测数据。在某一时刻时，通过电力系统量测装置采集得到原始数据组成的向量为z，设定需要数据清洗的时刻个数为N，即得到需要清洗的总数据集合为{z₁，z₂，z₃，…，z_N}。以k代表时刻，则需要利用k时刻的数据对k+1时刻进行预测，从而对k+1时刻的数据z_k+1进行清洗。

步骤2：建立电力系统的动态模型，得到电力系统动态方程和量测方程如下所示：

x_k+1＝f(x_k)+q_k (1)

z_k+1＝h(x_k+1)+r_k+1 (2)

式中：k表示时刻；x_k和x_k+1分别是k时刻和k+1时刻的状态向量，包括电力系统电压幅值、相角等；z_k+1是k+1时刻的量测向量，包括电力系统功率、电压等；f(·)和h(·)分别是电力系统的状态转移函数和量测函数；q_k和r_k+1分别为k时刻的过程噪声和k+1时刻的量测噪声，满足均值为0、方差阵分别为Q和R的正态分布。

步骤3：设定起始时刻k＝1，则k+1＝2，表示利用第1个时刻的前期准备数据，对第2个时刻的数据z₂进行清洗。

步骤4：建立容积卡尔曼滤波(cubature Kalman filter，CKF)基本方程，并基于CKF基本方程，得到数据清洗要用到的前期准备数据。首先计算方程(3)～(7)得到电力系统的电压、相角等状态量在k+1时刻的预测值再计算方程(8)～(12)得到电力系统的功率等量测数据在k+1时刻的预测值然后计算方程(13)得到电力系统量测数据的预测误差方差阵P_vv,k+1；最后计算方程(14)，将k+1时刻量测装置采集的量测数据z_k+1减去预测值得到量测新息向量e_k+1。

所述的CKF主要包括预测和滤波两个基本过程：

1)预测

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

式中：chol(·)表示Cholesky分解；P_k|k为k时刻状态估计误差方差阵；S_k|k为对P_k|k进行Cholesky分解得到的上三角矩阵；n为电力系统节点个数；I为n阶单位矩阵；ξ为根据Cubature准则生成的等权值矩阵；ξ_i表示矩阵ξ的第i列；为k时刻状态估计值；X_i,k|k为的Cubature点；为X_i,k|k的预测值；为状态预测值；Q_k为k时刻的过程噪声方差阵；P_k+1|k为状态预测误差方差阵。

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为状态量与量测量的协方差阵；W_k+1为滤波系数；为k+1时刻的电力系统状态估计值；P_k+1|k+1为k+1时刻状态估计误差方差阵。

步骤5：对k+1时刻的电力系统量测装置采集到的原始量测数据z_k+1进行不良数据辨识、修正过程，其流程图如图1所示。首先计算方程(20)～(27)，得到数据修正因子对角阵γ’_k+1。若γ’_k+1的对角元素等于1，表示该数据是正确数据，该数据对应的量测装置运行正常；若γ’_k+1对角元素大于1，表示该数据是不良数据，量测装置数据传输异常或者量测装置失效。将γ’_k+1代入方程(28)，计算得到修正后的电力系统量测总方差阵P_zz,k+1。

所述的不良数据辨识、修正过程如下：

5.1)在保证k+1时刻前的数据已不含不良数据的前提下，选取时间序列滑动窗口长度为T，将k+2-T时刻到k+1时刻的新息向量构成时间序列[e_k-T+2,e_k-T+3,...,e_k+1]，计算窗口内的平均新息协方差：

式中：P_e,k+1为平均新息协方差；T为时间序列滑动窗口长度；t代表时间序列的第t个向量，0≤t≤T，且t为整数；e_k+2-t为k+2-t时刻的新息向量。

当量测量突变时，由于e_k+1增大，必然使得平均新息协方差大于式(14)的观测噪声，即：

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

式中：P_vv,k+1为量测预测误差方差阵；R_k+1是量测噪声方差阵。

构造修正因子γ_k+1实现突变数据的辨识。γ_k+1可对R_k+1进行调整，其取值原则是使得上式左右两端相等，即：

求解得

则当可疑的突变量测量存在时，γ_k+1中的相应对角元素必将大于1。因此，取γ_k+1中大于1的对角元素对应的量测量，形成可疑数据集合D₁。D₁中的可疑数据即为突变数据，实际包含了正确的波动型数据和错误的不良数据。

5.2)选取时间序列滑窗长度为T，将k+2-T时刻到k+1时刻的量测向量构成量测量时间序列[z_k+2-T,z_k+3-T,...,z_k+1]，计算皮尔逊相关系数形成m阶相关矩阵，其对角元素为1，非对角元素计算如下：

其中：

式中：z_k+2-t为k+2-t时刻量测装置采集到的量测数据；为量测数据的平均值；r_ij为各量测量间的相关系数。

r_ij代表了第i个量测量与第j个量测量的相关性强弱，其范围为[-1,1]。|r_ij|越接近1，表明二者相关性越强；|r_ij|越接近0，表明二者相关性越弱。

计算第i个量测量与其他所有量测量的平均相关系数

式中：m为量测数据个数，亦为相关矩阵阶数。

若越小，则该量测量与其他量测量的总体相关性越弱，该量测量为不良数据的可能性越大。因此，为辨识可疑的不良数据，采用“末位淘汰”方法，选取与其他量相关性最弱的量测量形成可疑数据集合D₂。

采用定D₂元素个数的方法，其元素个数d₂的选取参考D₁的元素个数d₁，则可取d₂＝αd₁，α为可疑集长度参考倍数，可根据实际系统选取。

5.3)在得到上述可疑数据集合D₁和D₂之后，取二者交集得到不良数据集B＝D₁∩D₂，则B中的量测量才是真正的不良数据，其同时满足量测量发生突变和与整体量测数据相关性最弱两个条件。而B在D₁中的补集C_D1B实质为因新能源出力波动而产生的波动型数据，是正确的量测数据。

利用γ_k+1对不良数据进行修正。由于γ_k+1为非对角阵，为避免矩阵求逆运算发生奇异定义修正因子对角阵γ’_k+1，其对角元素取值为：

对式(15)进行如下修改以完成不良数据修正：

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

步骤6：对k+1时刻的原始量测数据z_k+1进行数据清洗。首先计算方程(16)～(18)，完成k+1时刻的电力系统状态估计过程，得到k+1时刻的状态估计值然后将代入方程(2)的量测函数，得到k+1时刻电力系统量测数据的清洗值

步骤7：计算方程(19)，得到k+1时刻状态估计误差方差阵P_k+1|k+1，以用作下一时刻数据清洗的前期准备数据。判断k的大小，若k≤N-1，则使k＝k+1，循环步骤4至步骤6，开始对下一时刻的电力系统量测数据进行清洗；若k>N-1，则数据清洗过程结束，得到数据清洗值的集合

本发明的效果和益处是：本发明提供了一种基于改进鲁棒容积卡尔曼滤波的数据清洗方法，该方法考虑了量测数据之间的相关性，通过将原鲁棒容积卡尔曼滤波的突变数据辨识过程与数据相关性强弱程度辨识过程进行结合，对波动型数据和不良数据进行区分，并通过修正因子对不良数据完成清洗，提高了数据清洗的准确性，将“脏”数据变为可靠数据。

附图说明

图1是上述步骤5的数据辨识修正过程。

图2是IEEE14节点系统拓扑图。

图3是风机出力的实际曲线与预测曲线。

图4是负荷系数曲线。

图5是场景(1)下，原始CKF、鲁棒CKF、改进鲁棒CKF三种方法的误差对比曲线。

图6是场景(2)下，原始CKF、鲁棒CKF、改进鲁棒CKF三种方法的误差对比曲线。

具体实施方式

以下以改进的IEEE14节点系统为例，拓扑图如图2所示，将2号节点的发电机处理为风力发电机。量测装置为数据采集与监控系统和同步相量测量单元，其中后者的配置节点为2、6、9。风机实际出力与预测功率曲线如图3所示，负荷曲线如图4所示。可疑集长度参考倍数α取2。

IEEE14节点系统的支路参数如表1所示，表中数据均为标幺值。

表1 IEEE14节点系统支路参数

结合技术方案和附图详细叙述本发明的具体实施方式。其具体步骤为：

步骤1：整理原始数据。设定需要数据清洗的时刻个数N＝96；根据风机出力与负荷数据，通过96次潮流计算，并将潮流计算值叠加服从均值为0的高斯分布误差，构成需要清洗的总数据集合为{z₁，z₂，z₃，…，z₉₆}。其中数据采集与监控系统量测误差的标准差为0.02，同步相量测量单元电压幅值量测误差的标准差为0.005，同步相量测量单元相角量测误差的标准差为0.002。

步骤2：建立电力系统动态状态方程为：

其中：

式中：J_k为潮流计算的雅可比矩阵；x_k和x_k+1分别为k时刻和k+1时刻的电力系统状态值；P_L,k和Q_L,k是k时刻的负荷有功和无功；和分别是k+1时刻的负荷预测有功和无功；ΔP_L和ΔQ_L分别为负荷有功和无功的预测偏差量；ΔP_G为发电机预测有功出力变化量，分为不可调度单元出力变化量ΔP_Gu与可调度单元出力变化量ΔP_Gd两类。

对于风电等不可调度的发电单元，其出力变化量为：

式中：P_Gu,k和分别为不可调度发电单元在k时刻的有功出力和k+1时刻的预测有功出力。

对于可调度的发电单元，其出力变化量为：

式中：n为节点个数；n_u为不可调度发电单元个数；β为各可调度发电机的负荷分配系数，没有可调度发电机的节点β(i)＝0；

步骤3：设定起始时刻k＝1，则k+1＝2，表示利用第1个时刻的前期准备数据，对第2个时刻的数据z₂进行清洗。

步骤4：基于CKF基本方程，得到数据清洗要用到的前期准备数据。首先计算方程(3)～(7)得到电力系统的电压、相角等状态量在k+1时刻的预测值计算方程(8)～(12)得到电力系统的功率等量测数据在k+1时刻的预测值然后分别计算方程(13)得到电力系统量测数据的预测误差方差阵P_vv,k+1；最后计算方程(14)，将k+1时刻量测装置测得的量测数据z_k+1减去其预测值得到新息向量e_k+1。

步骤5：对k+1时刻的原始量测数据z_k+1进行不良数据辨识、修正过程。首先计算方程(20)～(27)，取可疑集长度参考倍数α＝2，得到数据修正因子对角阵γ’_k+1；再将γ’_k+1代入方程(28)，计算得到修正后的电力系统量测总方差阵P_zz,k+1。

步骤6：对k+1时刻的原始量测数据进行数据清洗。首先计算方程(16)～(18)，完成k+1时刻的电力系统状态估计过程，得到k+1时刻的状态估计值然后将代入方程(2)的量测函数，得到k+1时刻电力系统量测数据清洗值

步骤7：计算方程(19)，得到为k+1时刻状态估计误差方差阵P_k+1|k+1。若k≤95，使k＝k+1，循环步骤4～步骤6，开始对下一时刻的电力系统量测数据进行清洗；若k>95，则数据清洗过程结束，得到量测数据清洗完成后的数据清洗值集合

为证明相比较于原始CKF(没有步骤5的方法)、鲁棒CKF(没有步骤5.2和步骤5.3的方法)，本方法(以下简述为改进鲁棒CKF)提高了数据清洗的准确性，采用数据清洗值的平均绝对误差作为性能指标进行对比，性能指标为：

式中：为k+1时刻量测数据清洗值；为k+1时刻量测数据真实值；m为量测数据个数。

数据清洗在两种场景下进行分别为：

场景(1)量测数据中不含不良数据时，原始CKF、鲁棒CKF、改进鲁棒CKF三种方法的误差曲线如图5所示。

场景(2)量测数据中含不良数据时，在32、40、48时刻分别添加个数不等、偏离程度不同的不良数据，模拟数据传输异常；在64-80时段，将两条支路的功率量测数据分别持续置零和扩大一倍，模拟量测装置异常。三种方法误差曲线如图6所示。

从图5可以直观地看出三种方法数据清洗值偏移真实值的程度。鲁棒CKF因会对量测正确的风电功率波动型数据进行错误清洗，因此其偏移量测真实值的程度最高；改进鲁棒CKF大大提高了数据辨识的准确度，避免数据错误清洗，其数据清洗性能与原始CKF基本相同。

从图6可以看到，虽然在不含不良数据时，原始CKF的数据清洗误差最小，但是在含不量数据时清洗误差急剧增大。改进鲁棒CKF无论原始数据是否含有不良数据，其数据清洗性能均优于原鲁棒CKF，清洗误差最小；同时，在含有不良数据时，改进鲁棒CKF的清洗误差同样大大低于原始CKF。

以上所述实施例仅表达了本发明的实施方式，但并不能因此而理解为对本发明专利的范围的限制，应当指出，对于本领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些均属于本发明的保护范围。

Claims (3)

1.一种基于改进鲁棒容积卡尔曼滤波的数据清洗方法，其特征在于，包括以下步骤：

步骤1：在某一时刻时，通过电力系统量测装置采集得到原始数据组成的向量为z，设定需要数据清洗的时刻个数为N，得到需要清洗的总数据集合为{z₁，z₂，z₃，…，z_N}；以k代表时刻，利用k时刻的数据对k+1时刻进行预测，进而对k+1时刻的数据z_k+1进行清洗；

步骤2：建立电力系统的动态模型，得到电力系统动态方程和量测方程；

步骤3：设定起始时刻k＝1，则k+1＝2，表示利用第1个时刻的前期准备数据，对第2个时刻的数据z₂进行清洗；

步骤4：建立容积卡尔曼滤波CKF基本方程，并基于CKF基本方程，得到数据清洗要用到的前期准备数据；首先计算方程(3)～(7)得到电力系统的状态量在k+1时刻的预测值再计算方程(8)～(12)得到电力系统的量测数据在k+1时刻的预测值然后计算方程(13)得到电力系统量测数据的预测误差方差阵P_vv,k+1；最后计算方程(14)，将k+1时刻量测装置采集的量测数据z_k+1减去预测值得到量测新息向量e_k+1；

所述的CKF主要包括预测和滤波两个基本过程：

1)预测

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

式中：chol(·)表示Cholesky分解；P_k|k为k时刻状态估计误差方差阵；S_k|k为对P_k|k进行Cholesky分解得到的上三角矩阵；n为电力系统节点个数；I为n阶单位矩阵；ξ为根据Cubature准则生成的等权值矩阵；ξ_i表示矩阵ξ的第i列；为k时刻状态估计值；X_i,k|k为的Cubature点；为X_i,k|k的预测值；为状态预测值；Q_k为k时刻的过程噪声方差阵；P_k+1|k为状态预测误差方差阵；

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为状态量与量测量的协方差阵；W_k+1为滤波系数；为k+1时刻的电力系统状态估计值；P_k+1|k+1为k+1时刻状态估计误差方差阵；

步骤5：对k+1时刻的电力系统量测装置采集到的原始量测数据z_k+1进行不良数据辨识、修正过程；首先计算方程(20)～(27)，得到数据修正因子对角阵γ’_k+1；若γ’_k+1的对角元素等于1，表示该数据是正确数据，该数据对应的量测装置运行正常；若γ’_k+1对角元素大于1，表示该数据是不良数据，量测装置数据传输异常或者量测装置失效；将γ’_k+1代入方程(28)，计算得到修正后的电力系统量测总方差阵P_zz,k+1；

所述的不良数据辨识、修正过程如下：

5.1)在保证k+1时刻前的数据已不含不良数据的前提下，选取时间序列滑动窗口长度为T，将k+2-T时刻到k+1时刻的新息向量构成时间序列[e_k-T+2,e_k-T+3,...,e_k+1]，计算窗口内的平均新息协方差：

式中：P_e,k+1为平均新息协方差；T为时间序列滑动窗口长度；t代表时间序列的第t个向量，0≤t≤T，且t为整数；e_k+2-t为k+2-t时刻的新息向量；

当量测量突变时，由于e_k+1增大，必然使得平均新息协方差大于式(14)的观测噪声，即：

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

式中：P_vv,k+1为量测预测误差方差阵；R_k+1是量测噪声方差阵；

构造修正因子γ_k+1实现突变数据的辨识；γ_k+1可对R_k+1进行调整，其取值原则是使得上式左右两端相等，求解得：

则当可疑的突变量测量存在时，γ_k+1中的相应对角元素必将大于1；因此，取γ_k+1中大于1的对角元素对应的量测量，形成可疑数据集合D₁；D₁中的可疑数据即为突变数据，实际包含正确的波动型数据和错误的不良数据；

5.2)选取时间序列滑窗长度为T，将k+2-T时刻到k+1时刻的量测向量构成量测量时间序列[z_k+2-T,z_k+3-T,...,z_k+1]，计算皮尔逊相关系数形成m阶相关矩阵，其对角元素为1，非对角元素计算如下：

其中：

式中：z_k+2-t为k+2-t时刻量测装置采集到的量测数据；为量测数据的平均值；r_ij为各量测量间的相关系数；

r_ij代表第i个量测量与第j个量测量的相关性强弱；

计算第i个量测量与其他所有量测量的平均相关系数

式中：m为量测数据个数，亦为相关矩阵阶数；

为辨识可疑的不良数据，采用“末位淘汰”方法，选取与其他量相关性最弱的量测量形成可疑数据集合D₂；

采用定D₂元素个数的方法，其元素个数d₂的选取参考D₁的元素个数d₁，则可取d₂＝αd₁，α为可疑集长度参考倍数，可根据实际系统选取；

5.3)在得到上述可疑数据集合D₁和D₂之后，取二者交集得到不良数据集B＝D₁∩D₂，则B中的量测量才是真正的不良数据，其同时满足量测量发生突变和与整体量测数据相关性最弱两个条件；而B在D₁中的补集C_D1B实质为因新能源出力波动而产生的波动型数据，是正确的量测数据；

利用γ_k+1对不良数据进行修正；由于γ_k+1为非对角阵，为避免矩阵求逆运算发生奇异定义修正因子对角阵γ’_k+1，其对角元素取值为：

对式(15)进行如下修改以完成不良数据修正：

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

步骤6：对k+1时刻的原始量测数据z_k+1进行数据清洗；首先计算方程(16)～(18)，完成k+1时刻的电力系统状态估计过程，得到k+1时刻的状态估计值然后将代入方程(2)的量测函数，得到k+1时刻电力系统量测数据的清洗值

步骤7：计算方程(19)，得到k+1时刻状态估计误差方差阵P_k+1|k+1，以用作下一时刻数据清洗的前期准备数据；判断k的大小，若k≤N-1，则使k＝k+1，循环步骤4至步骤6，开始对下一时刻的电力系统量测数据进行清洗；若k>N-1，则数据清洗过程结束，得到数据清洗值的集合

2.根据权利要求1所述的一种基于改进鲁棒容积卡尔曼滤波的数据清洗方法，其特征在于，所述的步骤2中电力系统动态方程和量测方程，如下所示：

x_k+1＝f(x_k)+q_k (1)

z_k+1＝h(x_k+1)+r_k+1 (2)

式中：k表示时刻；x_k和x_k+1分别是k时刻和k+1时刻的状态向量；z_k+1是k+1时刻的量测向量；f(·)和h(·)分别是电力系统的状态转移函数和量测函数；q_k和r_k+1分别为k时刻的过程噪声和k+1时刻的量测噪声。

3.根据权利要求1或2所述的一种基于改进鲁棒容积卡尔曼滤波的数据清洗方法，其特征在于，所述的步骤(24)中r_ij范围为[-1,1]。