基于时序聚类及改进Pettitt 的系统谐波阻抗估计

杨 勇，李文涛，刘文飞，张旭军，谢映洲，苗 虹

（1.国网甘肃省电力公司电力科学研究院，兰州 730070；2.四川大学电气工程学院，成都 610065）

高比例新能源和高比例电力电子设备的接入使系统谐波含量增加、频次增高，从而引起谐波分量放大及谐波分量超限的问题。为了应对这一挑战，精确测量系统谐波阻抗十分重要。但在电力系统实际运行情况下，由于风力发电容易受到风速影响，光伏电源容易受到光照强度影响，以及电力电子器件开断等原因，使背景谐波电压发生波动。此外，电力系统运行方式发生变化、电力电子器件的投切或增减负荷等原因均有可能导致系统谐波阻抗发生变化[1]，而背景谐波波动和阻抗变化会使谐波电压与电流数据出现严重无序性，给谐波阻抗的准确估计带来困难。

目前，谐波阻抗测量方法分为“干预式”和“非干预式”两种。“干预式”通过改变系统网络拓扑结构或向系统注入谐波电流、通断支路和启停负荷等措施，来估算系统谐波阻抗，但注入的扰动可能对系统运行产生不利影响[2]。“非干预式”法是在不干扰系统正常运行情况下，利用负荷或系统本身的自然扰动及可测量参数等来估算谐波阻抗，这类方法因对系统正常运行不产生影响而被广泛采用[3]。现有方法分为波动量法[4]、线性回归法[5−8]、独立分量法[9]和协方差法[10]等。

线性回归法由于计算精度较高而被广泛使用，但该方法更多用于背景谐波相对平稳状态，当背景谐波波动较大或冲击性负荷较多时，由于数据的无序性恶化，该类方法是否依旧具有稳健性还待验证[11]。为了规避背景谐波波动影响，文献[1]利用相关系数筛选数据，找到背景谐波波动较小的数据段来估计谐波阻抗，但没有考虑谐波电压和电流幅值序列存在偏移而带来的计算误差；文献[5]提出利用互近似熵保留波形相似的数据段，从而规避了背景谐波大幅波动的影响，该方法提供了新的数据筛选思路，但对噪声比较敏感；文献[7]通过动态时间弯曲距离DTW（dynamic time warping）法对谐波电压与电流数据进行筛选，抑制背景谐波波动影响，解决了谐波序列偏移问题，但计算量过大及计算时间过长。

谐波阻抗突变会使阻抗估计结果实时性较差且误差较大，如何有效识别阻抗突变点进行分段计算成为解决问题的关键。文献[1]利用基于密度聚类算法OPTICS（ordering points to identify the cluster⁃ing structure）将不同谐波阻抗分为不同簇分别计算，但当阻抗变化不明显时，难以准确进行划分；文献[12]利用谐波电压与电流的比值近似代替谐波阻抗值，通过设定阈值来判断谐波阻抗是否发生变化，但仅限于背景谐波稳定的场合；文献[13]使用小波变换模极大值法对谐波阻抗突变点进行检验，能够有效识别突变点，但没有考虑到噪声等原因对结果的影响。

本文提出一种基于时序聚类及改进Pettitt的系统谐波阻抗估计方法。首先，运用时序聚类的互相关系数原理来减小背景谐波波动影响，使谐波电压和电流序列平移对齐并筛选出两者相关性强的数据。然后，为规避阻抗突变对估计结果的影响，提出基于二元分割的Pettitt 方法，对谐波阻抗突变点进行检验，以实现对谐波阻抗的实时性估计。最后，使用再生权最小二乘法对分段后的数据进行谐波阻抗估计，通过赋权重系数进一步减小异常值对结果的影响。仿真与实例分析验证了该方法比传统线性回归等方法更具稳健性。

1 线性回归法基本原理

1.1 谐波阻抗计算等效模型

由于电力系统是一个复杂的网络，难以直接用来分析计算，故常用诺顿等效电路将电力系统分为系统侧和用户侧，在等效电路的公共耦合点PCC（point of common coupling）测量谐波电压与电流数据来估计系统谐波阻抗，等效模型如图1所示。

图1 诺顿等效电路Fig.1 Norton equivalent circuit

根据文献[12]中线性回归法基本原理，忽略相角影响，只取模值，则式（1）可表示为

1.2 谐波数据相关性筛选原理

背景谐波在电力系统实际运行情况下是一个随机波动过程，总会存在某些小时间段内没有波动或者波动幅度很小的情况。当背景谐波相对稳定且该时间段内的谐波阻抗不发生变化时，谐波电压幅值和谐波电流幅值呈现较强的正相关性。

图2 背景谐波波动情况Fig.2 Background harmonic fluctuations

2 时序聚类中互相关系数筛选原理

时序聚类互相关系数方法通过固定谐波电压幅值序列，平移谐波电流幅值序列，每移动一步计算两者距离测度SBD（shape based distance），当移动过程中两者SBD最小时两序列平移对齐，通过筛选对齐后SBD满足要求的子序列，可以减小背景谐波波动影响，其基本过程如下。

步骤2构建移动窗，设置移动窗的大小L和移动步长T。其中，移动窗大小L与步长T的最优参数通过多次实验获取。

（1）根据两子序列互相关关系，固定公共耦合点谐波电压幅值序列u͂的时间窗，平移另一条谐波电流幅值序列i͂的时间窗如图3所示。

图3 滑动过程Fig.3 Sliding process

式中，SBD 值范围为[0，2]。

步骤5滑动窗口向后移动一个距离单位T，重复步骤4，直至滑动窗口覆盖整个谐波数据序列。

本文利用互相关系数来计算谐波电压与电流幅值的相关性，SBD 越小两序列相关性越高。该方法有效解决谐波序列偏移问题，降低了背景谐波电压波动给系统谐波阻抗估计带来的误差。

3 基于改进Pettitt 方法阻抗突变检验

考虑到在实际情况中新型电力系统下系统侧谐波阻抗不再为定值，其值大小会因为无功补偿装置投切、系统运行方式改变、新发电设备的接入等原因而发生变化。因此在计算谐波阻抗时，应先判断其值是否发生变化，否则，不仅结果误差较大且计算得到的谐波阻抗实际意义并不明确。

3.1 谐波阻抗变化理想模型

当背景谐波波动较小时，谐波阻抗可近似估计为

系统谐波阻抗变化规律的简化如图4所示。

图4 阻抗突变等效模型Fig.4 Equivalence model of impedance mutation

图4 中，tn为谐波数据的采样总时间；t1,t2,…,tn−1为阻抗发生变化的时间点。由于在时间段(ti,ti+1)，i=1,2,…,n−1 内阻抗值不变，因此在该段时间内可建立线性回归模型求解系统谐波阻抗，而关键问题就是找到谐波阻抗突变点t1,t2,…,tn−1。

3.2 Pettitt 方法原理

Pettitt方法是一种对秩序列进行处理的非参数检验方法，其计算简便，可以明确突变时间点，并且能够较好的识别序列的突变点[15]。Pettitt方法的基本原理是假设将样本数据任意分为2 个子序列，2 个子序列均值相等，即，表明当随机变量序列x1,x2,…,xN被任意分割成x1,x2,…,xmp和xmp+1,xmp+2,…,xN两部分以后，如果两部分中随机变量的分布函数F1(x)≠F2(x)，则认为突变发生在mp 处，可定义为

式中：Ut,N为统计量，计算了第1个样本序列大于第2 个样本序列的次数；N为序列样本数；sgn(xt−xj)为阶跃函数，其可表示为

对应的显著性概率可表示为

一般认为当Pf＜0.01 时，满足显著性要求，即系统谐波阻抗在该点发生突变。

3.3 基于二元分割的Pettitt 方法

由于Pettitt突变检验方法只能识别长序列中最显著的一个突变点，但谐波阻抗序列中可能存在多个突变点，而且谐波阻抗序列较大时，该方法的计算量较大。本文提出引入二元分割法对Pettitt突变检验法进行改进，使其可以识别全部的突变点且减少计算量，具体步骤如下。

步骤1设由式（9）近似计算得到的谐波阻抗序列为x(t),t=1,2,…,N，为减小计算量，将谐波阻抗序列均分为s个子序列，每个子序列包含q个数据，本文取第一个子序列进行分析

用Pettitt 方法对谐波阻抗序列x[t1∶tq]进行突变检验，当子序列存在多个突变点时，记检测出的突变点位置为mp，由于检验出的突变点mp 可能会发生偏移，使检验结果不准确。突变点偏移如图5所示。

图5 突变点偏移Fig.5 Mutation point offset

步骤2为了判断检测出的突变点是否发生偏移，从突变点mp 进行分段，分别对分段后的子序列再进行Pettitt 检验，若分段后的两子序列检测不到新的突变点，则mp 为真突变点；若检验出新的突变点，则证明检验出的突变点mp 发生偏移，将检验出的新突变点记为i、j。检验为真突变点过程如图6所示。

图6 检验为真突变点Fig.6 Tested positive for mutation point

步骤3重复上述过程对两子序列进行Pettitt检验，直到检验不到突变点，则说明该子序列检测的为真突变点，检验结果准确无误。

同理对mp 段序列均用此方法进行计算求解突变点，并将所有检测得到的真突变点汇总为集合Mcp并排序，则此时集合Mcp作为谐波阻抗真突变点检测的最终结果。

4 基于再生权最小二乘法的谐波阻抗估计

4.1 再生权最小二乘法

根据突变点分段后，需要采取精准高效的方法对系统谐波阻抗进行估计，本文提出使用再生权最小二乘法SBWLS（self−born weighted least squares）对谐波阻抗进行估计。SBWLS 是一种基于观测值残差条件方程给予的信息来构建观测值权重的稳健回归估计方法，与传统最小二乘法不同，它不利用最小二乘法计算得到的观测值残差来构建观测值的权重，当数据存在多个异常值时，SBWLS 相较于一般的稳健回归方法，在消除或削弱异常值对参数估计的影响方面表现更为卓越[16]。

把PCC处谐波电流作为自变量，谐波电压作为因变量可列出线性回归观测方程为

误差方程可表示为

式中：V为观测值Y的残差；g为误差方程常数项，g=−(d+βX0−Y)。

式中：P为观测值Y的权阵，其元素为pj；r为多余观测值的数量，r=z−t,r＞0，z为观测值的数量。

将式（15）划分为两个部分，即

式中：βt为t×t阶的满秩矩阵，对系数矩阵β进行线性变换来定义；t为未知数的数量。

将式（19）变形可得

将式（21）代入到式（20）中可得再生权函数式为

由系统谐波阻抗线性关系式可得

式中：Ui为第i个谐波电压真实值；Vi为第i个谐波电压残差。

将式（29）矩阵化

在谐波阻抗线性回归模型中，有z组谐波数据观测值，其中未知量个数t=2，即谐波阻抗Zs与背景谐波电压Us，该回归模型的误差方程为式（15）。

由SBWLS 基本原理，对误差方程中的系数矩阵β进行线性变换并确定组成线性无关组的2 个误差方程式，则该模型误差方程式的分块矩阵可表示为

由式（27）可得再生权函数式为

以式（33）求得的观测值的再生权作为下一次迭代的观测值的权重，通过最小二乘法进行迭代计算，迭代的终止条件为相邻两次观测值残差小于0.1，则此时根据式（21）解得的背景谐波电压Us和系统谐波阻抗Zs逼近真实值。

4.2 谐波阻抗估计方法流程

根据本文方法原理，可得谐波阻抗估计方法流程如图7所示。首先，对采集到的谐波数据序列进行相关性检验，运用时序聚类中互相关系数原理来减小背景谐波波动影响，使谐波电压和电流序列平移对齐并筛选出两者相关性强的数据；然后，使用本文提出的基于二元分割的Pettitt方法对阻抗突变点进行检验；最后，使用SBWLS 对分段后的系统谐波阻抗进行计算。

图7 谐波阻抗估计流程Fig.7 Harmonic impedance estimation process

5 仿真与实例分析

5.1 仿真分析

为了验证本文方法的有效性，在Matlab上搭建如图8 所示的3 馈线电力系统模型进行仿真分析。以5次谐波为例，利用4种方法（方法1为波动量法、方法2 为线性回归法、方法3 为独立矢量法、方法4为本文方法）进行对照，估计系统侧谐波阻抗值。

图8 3 馈线系统Fig.8 Three-feeder system

为了模拟系统谐波阻抗突变，每4 s 改变系统谐波阻抗值，并使谐波源的幅值和相角在一定范围内波动；仿真共16 s，每0.004 s 测量一次谐波电压和电流值，仿真模型中谐波源与基波阻抗参数设置如表1所示。

表1 3 馈线系统仿真参数Tab.1 Simulation parameters of three-feeder system

仿真一共采集了4 000组PCC处的5次谐波电压和电流数据，如图9所示。其中，5次谐波电压幅值数据如图9（a）所示，电流幅值数据如图9（b）所示。

图9 谐波电压、电流幅值数据Fig.9 Data of harmonic voltage and current amplitude

利用互相关系数进行数据筛选，设置滑动窗口长度为40、滑动步长为5，计算相应窗口内谐波电压与谐波电流的SBD 值。从793组距离值中筛选出SBD值小于0.1的数据，再以谐波电流幅值为横坐标，以谐波电压幅值为纵坐标作散点图如图10所示。

图10 筛选前后谐波电压与电流散点图Fig.10 Scatter plot of harmonic voltage and current before and after screening

由图10 可知，未经互相关系数筛选前，谐波电压与电流幅值存在某种线性关系，但是受到背景谐波波动影响，这种线性关系并不明显；筛选后两者出现清晰的线性关系，故说明了所提方法能够有效消除背景谐波波动影响。筛选后谐波数据呈现4条直线，由谐波阻抗回归方程可知，直线斜率对应着4种不同谐波阻抗值，但需要从谐波数据中区分出这4条直线，然后分段进行计算。

本文通过式（9）计算出谐波阻抗近似值，接着使用基于二元分割的Pettitt 方法检验突变点，检验结果如图11所示。在图11中，UT为统计量，由式（10）得到，从而计算2个子序列中元素大小比较的总体情况，即，当UT 为正值时表示第1个子序列元素总体上大于第2个子序列，其为负值则正好相反；K为极值点，代表最显著的突变点，由式（12）得到。在采样点877、1 200、1 726、2 585 处，UT 曲线出现极值点K，由式（13）计算这些极值点的显著性水平Pf，可得利用采样点877、1 726、2 585 计算的Pf≈0，远小于给定的显著性水平Pf＜0.01，即说明在该极值点处谐波阻抗发生突变。这与谐波阻抗近似值中阻抗突变点一致，且3个极值点分别对应所设置的4 s、8 s、12 s 突变时间；而利用采样点1 200计算的Pf≈2，不满足给定显著性水平要求，不能说明该点谐波阻抗发生突变，检验结果与预设相同，故该方法能够有效检验出谐波阻抗突变点。此外，在4 s、8 s突变点的UT曲线先减小后增大，说明谐波阻抗在该点发生突变增大，而12 s处正好相反，说明谐波阻抗在极值点处发生突变减小，故该方法还可以预测阻抗突变趋势。

图11 3 馈线系统阻抗突变点检验结果Fig.11 Test results of three-feeder system impedance mutation point

由图11可以看出，曲线存在不光滑的区域，可能也会导致极值点K的出现，但本文所提Pettitt 方法只检测一段序列中最显著突变点，图中不光滑区域相比检测到的极值点来说是不显著的，而且可以通过多次实验得到给定显著性水平Pf的最优值，并设置分段最小序列阈值来避免不光滑区域带来的检测误差，以达到提高谐波阻抗突变检测精度的目的。

在检测出系统谐波阻抗突变时间点后，对筛选后的谐波电压、电流数据从突变点分段，再分别对这几段数据采用SBWLS 进行系统谐波阻抗估计，4种方法估计结果如表2所示。

表2 系统谐波阻抗估计结果Tab.2 Result of system harmonic impedance estimation

由表2 中可以看出，本文方法估计得到的系统谐波阻抗值与理论值的误差很小，且相对于其他方法精度更高，故本文方法在3馈线电力系统谐波阻抗计算上十分有效。

5.2 实例分析

为了进一步验证本文方法在电力系统实际运行情况下的有效性，选择使用标准的IEEE 14 节点系统来模拟电力系统复杂的运行环境从而进行仿真分析，如图12所示。

图12 IEEE 14 节点系统Fig.12 IEEE 14-bus system

IEEE 14节点模拟系统包括2台发电机、3台调相机、14 条母线节点、15 条输电线路及3 台变压器。各元件和负荷参数均按照IEEE 14节点的标准参数设定，设置系统额定电压为500 kV，估计特定谐波源作用时系统的谐波阻抗。

将母线11 作为研究母线，HL 则是研究母线负荷侧谐波源；背景谐波源负荷HS由母线6接入。本文对5次谐波进行分析，系统将母线11上的负荷看作用户侧负荷，除去用户侧负荷的其余部分均看作系统侧，通过仿真计算HL 谐波源负荷作用情况下该系统的等效谐波阻抗。

为了模拟无功装置投切、增减负荷等原因而导致的系统谐波阻抗发生变化，仿真实验分为3 部分：①测量系统正常运行情况下的谐波阻抗；②测量切除母线6上的同步调相机情况下的谐波阻抗；③测量减小母线6 谐波源负荷HS 阻抗情况下的谐波阻抗。仿真共48 s，每一部分16 s，每0.032 s通过傅里叶变换得到母线11 处的谐波电压和谐波电流幅值数据，共1 500个采样点。

同样利用互相关系数进行数据筛选，设置滑动窗口长度为40，滑动窗口从第1 个采样点开始，滑动一次，窗口向后移动5 个采样点，计算相应窗口内谐波电压与谐波电流的SBD 值。从293 组距离值中筛选出SBD 值小于0.1 的数据，再使用基于二元分割的Pettitt 方法进行阻抗突变检验，检验结果如图13所示。

图13 IEEE14 节点系统阻抗突变点检测结果Fig.13 Test result of IEEE 14-bus system impedance mutation point

由图13 可知，在采样点395、791 处UT 曲线出现极值K，由式（13）计算极值点的显著性水平Pf，其结果均小于给定的显著性水平，即Pf＜0.01，说明该极值点处谐波阻抗发生突变。这与谐波阻抗近似值中阻抗突变点一致，且2个极值点分别对应所设置的16 s、32 s 突变时间，故该方法在IEEE 14节点检验系统谐波阻抗突变点效果良好。检测出突变点后进行分段，使用SBWLS 估计每一段的谐波阻抗值，估计结果如表3所示。

表3 IEEE 14 节点系统谐波阻抗估计结果Tab.3 Result of IEEE 14-bus system harmonic impedance estimation

为了得到IEEE 14 节点系统谐波阻抗的理论值，将背景谐波源负荷HS电流置零，同时去除用户侧谐波源负荷HL 幅值和相角的扰动，即谐波源负荷HL单独稳定作用情况下测量得到的系统谐波阻抗为其理论值。由表3可以看出，估计得到的系统谐波阻抗值与理论值相差不大，这一仿真结果验证了本文方法的可靠性与有效性。

6 结 论

本文提出的基于时序聚类及改进Pettitt的系统谐波阻抗估计方法，结论如下。

（1）背景谐波大幅波动会使谐波电压与电流出现严重无序性，给阻抗估计结果带来误差，本文提出的时序聚类的互相关系数方法能够解决谐波电压和电流幅值序列偏移问题，并通过两者相关性筛选出背景谐波相对稳定的时段，有效地减小谐波阻抗估计的误差。

（2）本文提出的基于二元分割的Pettitt 方法能够对谐波阻抗发生突变的时间点进行有效识别，实现对系统谐波阻抗实时性计算。

（3）本文方法先筛选出背景谐波相对稳定的数据，然后使用SBWLS 估计谐波阻抗。仿真和实例分析结果验证了该方法与线性回归等传统方法相比更加准确，当异常值较多时更具稳健性。