基于鲁棒容积卡尔曼滤波器的发电机动态状态估计

毕天姝 陈 亮 薛安成 杨奇逊（新能源电力系统国家重点实验室（华北电力大学） 北京 102206）



基于鲁棒容积卡尔曼滤波器的发电机动态状态估计

毕天姝 陈 亮 薛安成 杨奇逊
（新能源电力系统国家重点实验室（华北电力大学） 北京 102206）

摘要同步相量测量单元（PMU）能够对电力系统动态过程中发电机功角进行直接量测。然而，坏数据有可能导致状态估计准确度下降甚至失效。提出了一种基于鲁棒性容积卡尔曼滤波（CKF）的机电暂态过程发电机动态状态估计方法。在CKF中构造时变多维观测噪声尺度因子，根据量测新息对PMU量测误差进行调整，使得量测量能够对状态量预报值进行准确修正。给出了时变多维观测噪声尺度因子的具体构造方法。针对滤波增益求逆发生奇异的问题，提出解决方案，对鲁棒CKF动态状态估计过程进行说明。仿真结果表明该方法能够有效抑制量测坏数据对动态状态估计的影响。

关键词：机电暂态 动态状态估计 容积卡尔曼滤波 鲁棒性

国家重点基础研究发展计划（973计划）（2012CB215206），国家自然科学基金（51222703），高等学校博士学科点专项科研基金（20120036110009）和“111”计划（B08013）资助项目。

Dynamic State Estimator for Synchronous Machines Based on Robust Cubature Kalman Filter

Bi Tianshu Chen Liang Xue Ancheng Yang Qixun
（State Key Laboratory of Alternate Electrical Power System with Renewable Energy Sources North China Electric Power University Beijing 102206 China）

Abstract Phasor measurement unit (PMU) can measure the rotor angle of synchronous machine in power system dynamic process. However, the bad data may decrease the accuracy of state estimations, even lead to the failure of the estimator. Based on the robust cubature Kalman filter (CKF), a novel dynamic state estimator for synchronous machine in the electromechanical transient process is proposed. A time-varying multi-dimensional scale factor is introduced into CKF. The PMU measurement covariance can be adjusted according to the innovation. As a result, the PMU measurements will correct the state predictions precisely. The formulation of the scale factor is clarified, and the method for dealing with the problem of the gain matrix singularity is addressed. The detailed process of dynamic state estimation based on robust CKF is given. The simulation results show that the method can prevent the influence of bad data on the precision of the dynamic state estimation.

Keywords：Electromechanical transient process, dynamic state estimator, cubature Kalman filter, robustness

0 引言

同步相量测量单元（Phasor Measurement Unit， PMU）的出现[1-3]为电力系统暂态稳定分析与控制提供了新的技术手段[4-9]。然而，在干扰、测量设备故障和同步信号丢失等情况下，往往导致坏数据的出现。坏数据可能使得PMU在应用过程中出现错误的分析结果和控制策略。状态估计能够剔除量测量中存在的坏数据，因此，研究基于PMU的电力系统机电暂态过程中动态状态估计至关重要。

近年来，基于PMU的机电暂态过程发电机动态状态估计问题已经引起学者广泛关注[10-14]。文献[10]建立了动态状态估计模型，并采用卡尔曼滤波进行求解。由于发电机动态方程为非线性函数，文献[11,12]采用扩展卡尔曼滤波（Extented Kalman Filter，EKF）求解动态状态估计模型。EKF的线性化过程导致截断误差过大，因此，文献[13,14]提出了基于无迹卡尔曼滤波（Unscented Kalman Filter，UKF）的动态状态估计。UKF需要确定参数值，灵活性差，应用不便。针对该问题，加拿大学者S. Haykin 于2009年提出容积卡尔曼滤波（Cubature Kalman Filter，CKF）[15,16]算法。然而，无论UKF还是CKF，当量测量存在坏数据时，估计精度都会在一定程度上受到影响，甚至导致估计结果严重偏离实际值，动态状态估计失败。尽管也有学者提出抗差状态估计方法[17]和坏数据的检测与辨识方法[18]，但是这些方法都针对最小二乘法静态状态估计。文献[19]提出一种基于EKF的鲁棒动态状态估计方法，但EKF截断误差过大的问题依然存在。

本文提出一种基于鲁棒CKF的发电机动态状态估计方法。将时变多维观测噪声尺度因子引入到CKF动态状态估计中，根据量测新息对量测误差进行在线调整，使其更加逼近真实噪声。再利用调整后的误差计算滤波增益，使其能够在量测量含有坏数据的情况下对状态量预报值进行准确修正，得到准确的状态量估计值。文章给出了时变多维观测噪声尺度因子的具体构造方法，针对滤波增益求逆发生奇异的问题，提出解决方案，对基于鲁棒CKF的动态状态估计过程进行了详细说明。将提出的鲁棒CKF动态状态估计应用于某实际大区域电网。结果表明，当量测量出现连续多点坏数据时，该方法依然能够得到准确的估计结果，且满足实时性的要求。

1 容积卡尔曼滤波

CKF是基于高斯假设的贝叶斯估计基本框架，针对贝叶斯估计中需要求解非线性状态方程与高斯概率密度乘积积分的问题，采用三阶球面-径向规则生成一组等权值的Cubature点对该积分进行逼近。以下对CKF的基本过程进行说明。

针对非线性动态系统

式中，x、u和z分别为状态向量、控制向量和量测向量；F和H分别为非线性状态方程和量测方程；v和w分别为过程噪声和量测噪声，满足均值为0、误差方差阵分别为Q和R的正态分布；下标k表示时刻。

CKF滤波算法包括两个基本过程为预报步和滤波步。预报步根据球面-径向规则生成状态量的Cubature点，从而利用状态方程对状态量进行一步预报，并计算预报误差方差；滤波步根据量测值对状态量预报值进行修正，从而得到更加准确的状态量估计值。具体步骤如下：

1）预报步

式中，Pk|k和分别为k时刻状态量估计误差方差阵和估计值；Sk|k为Pk|k的平方根矩阵；Xi,k|k为的 Cubature点；i=1,2,…,2n，n为状态量维数；[1]∈Rn，[1]表示点集

2）滤波步

式中，Sk+1|k为Pk+1|k的二次方根矩阵；Xi,k+1|k为的Cubature点；Zi,k+1|k为量测量预报值的Cubature点；为量测量预报值；Pvv,k+1为量测量预报误差方差阵；Pzz,k+1|k为量测量总方差；Pxz,k+1|k为交叉协方差阵；Wk+1为滤波系数。

由上述CKF基本步骤可见，CKF无需选择任何参数，应用灵活。然而，如果量测量出现坏数据，式（12）中量测误差方差阵R就会与实际量测误差不符，造成估计结果不准确。

2 多维观测噪声尺度因子

针对量测出现坏数据导致量测误差方差与实际值不符的问题，将一个时变多维观测噪声尺度因子引入CKF。由式（15）可见，CKF滤波算法需要根据新息向量对状态量的预报值进行修正，进而得到状态量估计值。修正程度由新息ek+1和滤波增益Wk+1共同决定。当实际量测噪声与给定的量测误差方差阵R相符时，ek+1和Wk+1能够对预报值进行正确修正，CKF能够得到准确的估计结果。然而，当量测量出现坏数据时，新息向量ek+1中坏数据对应的元素突然增大，而Wk+1并未随之进行调整，对状态量预报值的修正不准确，导致估计结果准确度下降。

由式（12）和式（14）可见，Wk+1是量测噪声方差阵R的函数，为了使Wk+1能够对状态量预报值进行正确修正，可以构造时变多维观测噪声尺度因子γk+1，对R进行在线调整。当量测量出现坏数据时，γk+1的值随之变化，进而调整Wk+1使其能够对状态量预报值进行准确修正。对式（12）进行如下修改

式中，γk+1为m维时变多维观测噪声尺度因子，m为量测量维数。由此，式（14）改为

式中，等号右边求逆部分为观测噪声方差，为使Wk+1能够准确修正状态量预报值，需要满足

当量测量没有坏数据时，γk+1为m维单位矩阵。当量测量出现坏数据时，新息ek+1会增大，γk+1不再为单位矩阵，具体取值情况需要另行计算。

利用开窗估计法计算新息实时协方差

式中，Pe,k+1为新息协方差；M为开窗估计法的窗长。当量测出现坏数据时，必然有Pe,k+1＞Pvv,k+1+Rk+1，则需要计算γk+1，进而对量测误差方差进行调整。γk+1的取值原则是使得新息方差与观测噪声相等，即

由式（21）求解得到γk+1

使用开窗估计法计算得到的新息协方差可能会导致式（22）计算的时变多维观测噪声尺度因子不是对角阵，从而使得式（18）等号右边矩阵求逆发生奇异。针对该问题，定义对角阵

式中，γk+1,ii为γk+1的第i个对角元素。由于为对角阵且对角元素均大于0，因此，将作为时变多维观测噪声尺度因子计算式（18）的滤波增益就可以保证求逆不发生奇异。则滤波增益计算公式为

由式（25）可见，滤波增益中引入时变多维观测噪声尺度因子，可以根据量测新息对量测误差方差矩阵进行调整。当量测量出现坏数据后，在观测噪声尺度因子的作用下，滤波增益减小，进而降低坏数据对估计结果的影响。

3 基于鲁棒CKF的发电机动态状态估计

系统动态过程中拓扑改变，难以实时获取，系统状态量发生突变，无法预测。因此，利用PMU实测的发电机出口电压相量实现动态方程与外部网络的解耦，选择不会发生突变且满足动态方程的发电机状态量作为动态状态估计的状态量，发电机动态状态估计的状态方程就为发电机动态方程。由于现有控制设备难以对次暂态过程做出及时反应，因此，忽略次暂态，动态方程由六阶降为四阶。发电机动态状态估计状态方程为

式中，δ 和ω 分别为发电机转子绝对功角和电角速度标幺值；TJ为发电机惯性时间常数；Tm为机械转矩；U和φ 分别为发电机出口电压相量的幅值和相角；和分别为q轴和d轴瞬变电抗；和分别为q轴和d轴瞬态电动势；D为阻尼系数；和分别为q轴和d轴开路瞬变时间常数；Ef为定子励磁电动势；Xq和Xd分别为q轴和d轴同步电抗。

PMU能够对发电机功角和角速度进行直接测量。利用GPS可以得到电磁功率的同步量测值。因此，测量量包括发电机功角、电角速度和输出电磁功率。动态状态估计的量测方程为式中，δz、ωz和Pe分别为转子绝对功角、电角速度和输出电磁功率的PMU量测值。对照式（1）中状态方程有和Ef在机电暂态过程中假设恒定，U和φ 采用PMU量测值。

利用导数的基本概念将式（26）的状态方程进行离散化，变为式（1）中差分方程的形式，就可以利用鲁棒CKF对发电机动态状态估计模型进行求解。

下面给出基于鲁棒CKF的发电机动态状态估计具体步骤。

1）预报步

（1）由k-1时刻估计结果得到k时刻发电机状

（3）根据式（4）中差分形式的发电机动态方程对发电机状态量的每一个Cubature点进行变换，得

（4）根据式（5）对所有发电机状态量Cubature点的预报值进行加权求和，进而得到状态量预报值。利用球面-径向规则生成的Cubature点权重为1/m，m=2n。动态状态估计状态量维数为4，则Cubature点权值为1/8。

（5）根据式（6）计算预报误差协方差阵Pk+1|k。

2）滤波步

（2）根据式（8）生成发电机状态量预报值的Cubature点Xi,k+1|k。

（3）根据式（9）对每个状态量预报值Cubature点进行变换，得到PMU量测量预报值的Cubature 点Zi,k+1|k。

（4）根据式（10）对所有量测量预报值的Cubature点进行加权求和，进而得到量测量预报值根据球面-径向规则，Cubature点的权重为1/8。

（5）得到k+1时刻PMU量测值zk+1= (δzk+1ωzk+1Pek+1)T后，计算新息向量再分别根据式（11）和式（22）计算PMU量测量预报误差方差阵Pvv,k+1和时变多维观测噪声尺度因子γk+1，再根据式（23）和式（24）计算对角阵。

（6）根据式（13）计算发电机状态量预报值和PMU量测量预报值之间的互协方差矩阵Pxz,k+1|k，再根据式（25）计算卡尔曼滤波增益Wk+1。

（7）利用k+1时刻的新息ek+1，通过卡尔曼滤波增益Wk+1对发电机状态量预报值进行滤波，

即⁀根据式（15）计算得到k+1时刻状态量估计值

xk+1|k+1。

（8）计算发电机状态量估计误差方差阵Pk+1|k+1得到k+1时刻状态量估计值和估计误差方差阵Pk+1|k+1后，令k=k+1，并返回预报步，进行下一次迭代。

4 仿真分析

将本文所提出的基于鲁棒CKF的发电机动态状态估计应用于某实际大区域电网。该电网中含有96个厂站节点、19台机组和187条线路。所研究的发电机参数取值情况如下：=0.196，=0.637，Xq=1.836，Xd=1.836，=7.5，=1.2，D=2， TJ=6.059 5。在t=0时某线路发生三相金属性短路故障，两周波后故障消失。利用BPA仿真软件计算结果作为真值，在真值基础上叠加高斯白噪声作为量测值，在t=4s时对功角量测设置12个连续坏数据。仿真时间为10s，仿真步长Δt为1周波，即0.02s。

分别采用本文提出的鲁棒CKF滤波算法以及CKF滤波算法对发电机动态状态估计进行求解。基于这两种方法的发电机状态量估计结果如图1所示。可见，在t=4s时，由于发电机功角量测存在坏数据，在基于CKF的发电机动态状态估计中，量测误差方差矩阵与功角的实际量测误差不一致，从而使得滤波增益无法对功角预报值进行准确修正，最终导致功角估计结果偏离真值。由于滤波增益还包含了互协方差矩阵Pxz,k+1|k，因此，功角量测的坏数据还会对转子电角速度和发电机内电动势预报值的修正产生影响，从而导致估计值产生较大偏差。

图1 发电机动态状态估计结果Fig.1 The results of dynamic state estimators for synchronous machines

然而，基于鲁棒CKF的发电机动态状态估计在功角量测存在坏数据的情况下仍然能够得到较为准确的估计结果。这是因为当坏数据出现时，鲁棒CKF会根据量测新息计算时变多维观测噪声尺度因子，进而调整滤波增益，使其能够对发电机状态量的预报值进行准确修正。

图2 发电机功角的观测噪声尺度因子Fig.2 The scale factor of power angles

图2给出了发电机功角量测的时变观测噪声尺度因子。可见，在t=4s时，功角量测出现连续多点坏数据，时变观测噪声尺度因子的值也会相应增大。这是由于当坏数据出现后，坏数据点对应的新息向量中元素值会突然增大，根据新息计算得到的观测噪声尺度因子就会增大，进而修正量测误差，降低坏数据权值，避免坏数据对状态估计结果的影响。

为了获取量化的估计指标，定义滤波系数η 以及指标ε 为

表1 动态状态估计指标Tab.1 The indexes of the dynamic state estimation

将基于CKF和鲁棒CKF两种方法的发电机动态状态估计分别用于实际大区域电网和IEEE 9节点测试系统。表1给出了量测量不存在坏数据情况下，两种方法的估计指标。可见，这两种发电机动态状态估计方法的滤波系数均小于1，说明两种方法都能够有效滤除量测量中存在的随机误差。除了IEEE 9节点测试系统中发电机G1的d轴暂态电动势和发电机G2的电角速度以外，基于鲁棒CKF的动态状态指标均小于CKF动态状态估计。这是由于鲁棒动态状态估计中的时变多维观测噪声尺度因子在正常量测条件下取值基本为1，这与容积卡尔曼滤波的估计过程基本一致。此外，CKF和鲁棒CKF的单步计算时间分别为0.206ms和0.272ms。尽管鲁棒CKF需要计算时变多维观测噪声尺度因子，导致计算效率略低于CKF，但是仍然能够满足实时性需求，并且估计准确度明显提升。

5 结论

针对PMU量测量中坏数据导致动态状态估计结果偏离真实值的问题，提出了一种基于鲁棒CKF的抗坏数据动态状态估计方法。将时变多维观测噪声尺度因子引入到CKF滤波算法中。采用某实际大区域电网进行仿真的结果表明，提出的鲁棒CKF动态状态估计方法在量测量存在坏数据的情况下依然能够得到准确的估计结果。此外，虽然鲁棒CKF在计算效率上略低于CKF，但是依然能够满足实时性需求。

参考文献

[1] Phadke A G. Synchronized phasor measurements in power systems[J]. IEEE Computer Applications in Power, 1993, 6(2): 10-15.

[2] Burnett R O, Phadke A G, et al. Synchronized phasor measurements of a power system event[J]. IEEE Transactions on Power System, 1994, 9(3): 1643-1650.

[3] Phadke A G. Synchronized phasor measurements-a historical overview[C]//IEEE/PES Transmission and Distribution Conference and Exhibition: Asia Pacific, 2002, 1: 476-479.

[4] 秦晓辉, 毕天姝, 杨奇逊, 等. 基于WAMS动态轨迹的电力系统功角失稳判据[J]. 电力系统自动化, 2008, 32(23): 18-22.

Qin Xiaohui, Bi Tianshu, Yang Qixun, et al. Power system transient stability assessment based on WAMS dynamic trajectories[J]. Automation of Electric Power Systems, 2008, 32(23): 18-22.

[5] 徐慧明, 毕天姝, 黄少锋, 等. 基于WAMS的潮流转移识别算法[J]. 电力系统自动化, 2006, 30(14): 14-19.

Xu Huiming, Bi Tianshu, Huang Shaofeng, et al. WAMS based flow transfer identification algrithm[J]. Automation of Electric Power Systems, 2006, 30(14): 14-19.

[6] 杨东俊, 丁坚勇, 邵汉桥, 等. 基于WAMS的负阻尼低频振荡与强迫功率振荡的特征判别[J]. 电力系统自动化, 2013, 37(13): 57-62.

Yang Dongjun, Ding Jianyong, Shao Hanqiao, et al. WAMS based characteristic discrimination of negative damping low-frequency oscillation and forced power oscillation[J]. Automation of Electric Power Systems, 2013, 37(13): 57-62.

[7] 宋洪磊, 吴俊勇, 郝亮亮, 等. 基于WAMS和改进拉普拉斯特征映射的同调机群在线识别[J]. 电网技术, 2013, 37(8): 2157-2164.

Song Honglei, Wu Junyong, Hao Liangliang, et al. On-line identification of coherent generator based on WAMS and improved Laplacian eigenmap algorithm[J]. Power System Technology, 2013, 37(8): 2157-2164.

[8] 马世英, 刘道伟, 吴萌, 等. 基于WAMS及机组对的电网暂态稳定态势在线量化评估方法[J]. 电网技术, 2013, 37(5): 1323-1328.

Ma Shiying, Liu Daowei, Wu Meng, et al. An online quantitative evaluation method of power system transient stability situation based on WAMS and generator pair[J]. Power System Technology, 2013, 37(5): 1323-1328.

[9] 丁涛, 董柏峰, 顾伟, 等. 基于PMU 的电压稳定动态线性化指标优化切负荷算法[J]. 电力系统保护与控制, 2013, 41(9): 27-33.

Ding Tao, Dong Baifeng, Gu Wei, et al. Optimization of a load shedding scheme using dynamic voltage stability linearized index based on PMU[J]. Power System Protection and Control, 2013, 41(9): 27-33.

[10] 秦晓辉, 毕天姝, 杨奇逊. 基于WAMS的电力系统机电暂态过程动态状态估计[J]. 中国电机工程学报, 2008, 28(7): 19-25.

Qin Xiaohui, Bi Tianshu, Yang Qixun. Dynamic state estimator based on WAMS during power system transient process[J]. Proceedings of the CSEE, 2008, 28(7): 19-25.

[11] Huang Zhenyu, Schneider K, Nieplocha J. Feasibility studies of applying kalman filter techniques to power system dynamic state estimation[C]//International Conference on Power Engineering, 2007: 376-382.

[12] Ghahremani E, Kamwa I. Dynamic state estimation in power system by applying the extended Kalman filter with unknown inputs to phasor measurements[J]. IEEE Transations on Power Systems, 2011, 26(4): 2556-2565.

[13] Ghahremani E, Kamwa I. Online state estimation of a synchronous generator using unscented Kalman filter from phasor measurements units[J]. IEEE Transactions on Energy Conversion, 2011, 26(4): 1099-1108.

[14] Wang Shaobu, Gao Wenzhong, Meliopoulos A P S. An alternative method for power system dynamic state estimation based on unscented transform[J]. IEEE Transactions on Power System, 2012, 27(2): 942-950.

[15] Arasaratnam I, Haykin S. Cubature Kalman filters[J]. IEEE Transactions on Automatic Control, 2009, 54(6): 1254-1269.

[16] Arasaratnam I, Haykin S, Hurd T R. Cubature Kalman filter for continuous-discrete systems: theory and simulations[J]. IEEE Transactions on Signal Processing, 2010, 58(10): 4977-4993.

[17] 巨云涛, 吴文传, 张伯明. 支持大规模电流量测的配网抗差状态估计方法[J]. 中国电机工程学报, 2011, 31(19): 82-89.

Ju Yuntao, Wu Wenchuan, Zhang Boming. A new method for distributed state estimation accommodating current measurements[J]. Proceedings of the CSEE, 2011, 31(19): 82-89.

[18] 李碧君, 薛禹胜, 顾锦汶, 等. 基于快速分解正交变换状态估计算法的坏数据检测与辨识[J]. 电力系统自动化, 1999, 23(20): 1-4.

Li Bijun, Xue Yusheng, Gu Jinwen, et al. Detection and identification of bad data in state estimation using fast decoupled orthogonal transformation[J]. Automation of Electric Power Systems, 1999, 23(20): 1-4.

[19] Shih K R, Huang S J. Application of a robust algorithm for dynamic state estimation of a power system[J]. IEEE Transaction on Power System, 2002, 17(1): 141-147.

毕天姝 女，1973年生，教授，博士生导师，研究方向为电力系统保护与控制、WAMS应用及故障诊断。

E-mail: tsbi@ncepu.edu.cn（通信作者）

陈 亮 男，1984年生，博士研究生，研究方向为WAMS高级应用。E-mail: ch.lg@163.com

作者简介

收稿日期2014-01-12 改稿日期 2014-04-09

中图分类号：TM711