基于相关性分析的电网非同步监测数据场景谐波责任划分

编者按

“双碳”目标下，能源结构发生显著变化，分布式能源与非线性负荷大量接入电力系统，谐波源数量剧增，谐波污染问题越发严重和复杂。为有效治理谐波，国际上提出了谐波治理的“奖惩性方案”，而准确的谐波责任划分是该方案实施和有效谐波治理的重要前提。

《中国电力》2025年第1期刊发了陈仕龙等撰写的《基于相关性分析的电网非同步监测数据场景谐波责任划分》一文。文章提出一种综合考虑数据非同步性、场景划分和数据相关性的谐波责任划分方法。首先，利用分段聚合近似算法进行降噪预处理，而后利用形状动态时间规整算法（shape dynamic time warping，ShapeDTW）处理数据间的非同步性问题；其次，利用点排序识别聚类结构的聚类算法（ordering points to identify the clustering structure，OPTICS）进行划分场景，分场景讨论谐波责任；最后，采用大数据分析思想，利用相关性分析方法构建各场景谐波责任和总谐波责任指标，并将各场景时长占比考虑在内。通过仿真分析和电网实例分析对本文方法进行验证，其谐波责任结论具有准确性与合理性，可进一步进行工程应用验证。

摘要

针对传统谐波责任划分方法需采用专门同步设备监测数据，且需基于等值电路模型划分谐波责任，工程应用较为复杂等不足，采用现有谐波监测装置非同步测量数据，提出一种综合考虑了数据非同步性、场景划分和数据相关性的谐波责任划分方法。首先，对原始非同步监测数据集采用分段聚合近似算法进行降噪预处理，利用形状动态时间规整算法（shape dynamic time warping，ShapeDTW）实现数据匹配对齐；然后，利用点排序识别聚类结构的聚类算法（ordering points to identify the clustering structure，OPTICS）划分场景以处理电力系统中因负荷投切和无功补偿装置切换等情况导致的谐波责任变化；最后，基于相关性分析构建场景谐波责任和总谐波责任指标，在指标构建的过程中引入了场景时长占比这一因素以得到更加科学合理的总谐波责任值。通过仿真验证和电网实例验证，该方法能基于现有非同步性监测数据实现各用户合理时间尺度动态谐波责任划分，可为工程上的快速谐波责任划分提供一定的新思路和新方法。

01 监测数据非同步和谐波阻抗变化下的相关性分析

1.1  谐波监测数据的相关性分析

当电力系统中存在多个谐波源分散分布时，任一关注母线上的谐波电压畸变都是由所有谐波源注入谐波电流而引起的共同结果，其多谐波责任划分通常可等效为如图1所示的模型。

图1  多谐波源电力系统等效模型

Fig.1  Equivalent model of multi-harmonic source power system

以云南电网某变电站的谐波监测数据为例，分别对3条馈线和母线处采集24 h的7次谐波电流和7次谐波电压，其3条馈线的谐波电流和母线谐波电压数据变化趋势如图2所示。可以看出，母线处的谐波电压和各馈线谐波电流变化趋势存在一定的相关性，但不同馈线谐波电流与母线谐波电压的相关性存在明显差异，且不同时间段相关性存在变化。

图2  3馈线7次谐波电流与母线7次谐波电压变化趋势

Fig.2  Variation of 7 th harmonic current of three feeders and 7 th harmonic voltage of bus

以云南电网某变电站严格对齐且无明显谐波阻抗变化的监测数据为例，关注母线处的7次谐波电压幅值与所接入的馈线谐波源的7次谐波电流幅值之间的线性相关性散点图如图3所示，其斜率与馈线的转移谐波阻抗Z_sh有关，截距与非关注谐波源共同产生的谐波电压U_0h有关，多个谐波源存在多个线性相关性。

图3  谐波源线性相关性散点图

Fig.3  Linear correlation scatter plot of harmonic sources

从以上分析可看出，在数学上母线谐波电压近似由网络内所有谐波源线性组合而成，不同母线受各谐波源的影响存在差异。受同一种或多种谐波源影响的母线，其母线谐波电压数据与各馈线谐波电流存在一定的关联性特征，在时间序列上体现为数据间波动的相似性。

1.2  谐波监测数据的非同步性

电能质量监测装置监测点如图4所示，监测装置一般部署在10 kV母线公共连接点上，数据采样间隔为3 min，可获得母线谐波电压数据以及m条馈线的谐波电流数据。

图4  谐波监测数据采集示意

Fig.4  Harmonic monitoring data acquisition schematic

由于母线谐波电压采集和馈线谐波电流采集分属不同监测装置，而电能质量监测装置以本地时钟为参考基准进行数据采集，这就造成了不同监测点数据采集的非同步性。同时，非同步测量下现有监测装置难以获取谐波电压电流瞬时值。

而目前电网公司所使用的谐波监测装置输出的测量数据一般为监测周期内的统计值，如最大值、最小值、平均值和95%概率大值，母线谐波电压能反映电力系统及供电用户谐波综合作用的最终结果，因此常选取谐波电压统计值进行谐波责任分析。不同组统计数据间常常会在数据时间上出现错位和偏移，进一步加剧了谐波监测数据间的非同步性。

若以在正常工况下采集图3数据来源的变电站母线谐波电压监测数据和馈线谐波电流监测数据，以馈线谐波电流幅值为横坐标、母线谐波电压幅值为纵坐标的散点图如图5所示。可以看出，非同步采样下其监测数据不再线性分布，无法衡量出数据间的相关性程度。即使通过网络校时等方法来校正时钟，也很难实现数据间的完美同步。

图5  非同步采样下谐波监测数据散点图

Fig.5  Scatter plot of harmonic monitoring data by asynchronous sampling

1.3  含谐波阻抗变化的谐波监测数据

在实际运行的电力系统中，当谐波监测装置数据监测周期较长时，谐波阻抗可能因系统运行方式改变、负荷投切和无功补偿装置切换等情况发生改变，将对如何准确划分谐波源s应承担的谐波责任产生较大影响，不同时段谐波源s的应承担谐波责任将差别较大。

延长图3数据的监测周期，将该变电站含运行方式改变和负荷投切的时段纳入数据监测周期，其监测数据散点图如图6所示。可以看出，其数据间的相关性发生了显著变化，不能再以一个固定相关性系数区衡量该时段内谐波源应承担的谐波责任。

图6  谐波阻抗变化下谐波监测数据散点图

Fig.6  Scatter plot of harmonic monitoring data under harmonic impedance change

02 本文谐波责任划分方法

2.1  数据预处理

电能质量监测数据具有高噪声等特点，直接使用原始数据进行谐波责任划分会导致计算结果精确度不高，通过对谐波监测数据进行降噪可有效改善该问题。本文采用分段聚合近似算法（piecewise aggregation approximation，PAA）对谐波监测数据进行预处理。将谐波监测数据表示为时间序列v={v₁,v₂, ···, v_i, ···, v_m}，v_i表示第i个监测数据，m表示序列长度。

采用文献[21]中经典PAA算法对谐波监测数据进行降噪处理，即

式中：ω为时间窗口长度；为预处理后的第j个谐波监测数据；PAA降噪处理后的数据为长度为n。经PAA处理后，数据可留存原有信息，且数据噪声显著降低。

2.2  基于ShapeDTW算法的数据对齐

电能质量监测数据是一种典型的时序数据，不同监测点采集的数据序列存在局部位移。动态时间规整算法（dynamic time warping，DTW）能够对各时序数据点进行非同时刻映射，有效处理序列中的局部位移现象，衡量2个非时间序列之间的相似程度。在进行存在局部位移的曲线匹配时，DTW距离对应与传统欧式距离对应对比如图7所示，可见欧式距离在度量存在时移但局部相似的曲线时并不适用，而DTW距离可准确度量其对应关系并进行匹配。

图7  DTW距离对应与欧式距离对应对比

Fig.7  Comparison of DTW distance correspondence and Euclidean distance correspondence

通过求解最优匹配路径和对齐方式，可得到母线监测谐波电压序列和馈线监测谐波电流的距离矩阵，从而定量分析两者之间的相关性。设馈线谐波电流序列和母线谐波电压序列分别为x={x₁, x₂, ···, x_m}和y={y₁, y₂, ···, y_n}，序列长度分别为m和n。首先，构造一个m行n列的距离矩阵M，其中M[i,j]表示序列x的第i个数x_i与序列y的第j个数y_j的欧氏距离。其次，设累计距离矩阵为M_c，其第1行和第1列初值为

累计距离矩阵其余部分计算方法为

式中：2≤i≤m，2≤j≤n，i、j∈N。

最后，确定序列x和y的DTW距离。由式（3）可知，累计距离计算的过程相当于计算序列x和y的最优对齐方式，并将最优匹配方式下的累计距离计入矩阵末尾，即可得2组序列数据的整体最小累计距离D(x,y)= M_c[i,j]。其值代表2组序列在趋势特征和时间特征上的相似程度，值越大则相似度越高。

在实际工程应用中，利用DTW算法进行母线谐波电压和馈线谐波电流点对点数据匹配时，会出现图8中2点之间距离最小但忽略局部形状相似度的不合理匹配，导致A点与B′点匹配。但从图8中2组数据的形状相似度对比，可明显看出A点应与

图8  合理匹配与不合理匹配示意

Fig.8  Schematic of reasonable and unreasonable matching

为有效避免此种类似不合理匹配情况的出现，使具有相似局部形状的监测数据序列点趋于匹配，本文在DTW算法匹配的基础上，采用ShapeDTW算法将序列点周围的局部形状信息合并到动态规划匹配过程中，实现数据匹配对齐。ShapeDTW算法实现数据匹配对齐步骤如下。

1）本文进行的谐波责任划分是研究2序列x和y之间的责任关系，此处选择以序列x为基准序列，序列y为待对齐序列。以2序列每个元素为中心，截取长度分别为L_x(L_x≪m)和L_y(L_y≪n)的子序列，可得序列x的子序列矩阵X′(m×L_x维)和序列y的子序列矩阵Y′(n×L_y维)。

2）对于矩阵X′的每行数据序列x′(i)(i=1, 2, ···, m)，按如式（4）所示原则与矩阵Y′的每行数据序列y′(j)(j=1, 2, ···, n)进行最优匹配，设式（4）在处取得最小值，即表示序列x的元素x_i与序列y的元素y_j匹配。

3）保存基准序列x不变，将与x_i匹配的y_j分配至新序列y′中的第i个元素，即按以上规则对序列y重构，形成新序列y′。

2.3  基于OPTICS聚类的场景划分

采用本文方法将非同步的馈线谐波电流数据和母线谐波电压数据对齐后，可实现数据间的同步性。但针对1.3节中含谐波阻抗变化的谐波监测数据，如未考虑谐波阻抗变化对谐波责任划分的影响，其谐波责任划分结果的合理性与适用性将存疑，对此这里将不同谐波阻抗对应的运行场景划分为不同的簇，然后再对各场景下的谐波源进行谐波责任划分。

OPTICS聚类算法是一种基于密度的聚类算法，在DBSCAN聚类算法的基础上改进得到。DBSCAN算法存在易受输入领域参数(ε, β)的影响，当参数取值不同时，DBSCAN聚类结果也不同的情况。而OPTICS算法通过生成一个反映各样本点基于密度聚类结构的有序队列实现灵活聚类。从理论上来讲，OPTICS算法可对不同密度的数据进行聚类，获得任意可分形状的簇。

假设数据集为X={x₁, x₂, ···, x_i, ···, x_n}，邻域半径为ε，最小数目为β。核心对象定义为：如果c_a(x_i)≥β ，则x_i为X的核心对象点，其中c_a(x_i)表示点x_i邻域内所包含的元素数。直接密度可达定义为：若x_i属于x_j的邻域且x_j为核心对象点，则称x_i是x_j的直接密度可达点。核心距离定义为：x_i的核心距离为使x_i成为核心对象点的最小领域半径。可达距离定义为：x_j关于x_i的可达距离为x_i的核心距离和x_j与x_i之间欧几里得距离的距离最大值。以谐波阻抗为聚类依据的基于OPTICS聚类算法的场景划分步骤如下。

1）遍历样本集E中的元素，判断该元素是否为核心对象，是则归入到集合Ω中，否则继续判断下一个元素，直至遍历全部元素。

2）在集合Ω中任意选取一个未被处理的对象点p，将该点标为已处理，寻找该点所有直接密度可达点，并按可达距离大小，将所有直接密度可达点依次排序存入集合S中。

3）若S为空集，则返回步骤2）；若S不为空，则选取集合S中可达距离最小的样本点q，标记为已处理，将该点存至有序列表M中，并判断点q是否为核心对象点，是则继续步骤4），否则返回步骤3）。

4）寻找点q所有直接密度可达点a_q(j)，若a_q(j)已存在M中则不做处理；否则判断S中是否已存在a_q(j)，存在则继续步骤5），不存在则跳到步骤6）。

5）若此时关于当前对象新的可达距离小于旧可达距离d_r(i)，则将其对应的可达距离替换为对S按可达距离重新排序，返回步骤3）。

6）插入点a_q(j)，对S按可达距离重新排序，返回步骤3）。

7）将样本集E中所有元素按步骤2）~6）处理完毕。

以处理顺序为横坐标，可达距离d_r(i)为纵坐标，生成有序队列图。根据有序队列图选取合适的邻域半径ε，若d_r(i)＜ε，则该可达距离有效，将其聚为一类，输出低谷数据，得到最终聚类划分结果。OPTICS聚类完成后，监测数据集被划分为不同场景的数据簇，可在各簇数据基础上利用相关性分析方法进行谐波责任划分。

03 谐波责任指标

为更有效地刻画出第1节中的馈线谐波电流和母线谐波电压的线性关系程度，本文采用相关性分析方法进行谐波责任划分。母线谐波电压和各馈线谐波电流变化趋势具有一定的相关性，不同馈线不同运行场景对应的相关性存在显著差异，需定量分析其相关性系数大小，进而根据相关性系数大小确定各谐波源的谐波责任。用变量之间的相关程度来描述变量间的线性关系，其相关性可以用系数r(x,y)为

式中：x_i和y_i分别为馈线谐波电流序列x的第i个元素和母线谐波电压序列y的第i个元素；分别为序列x的均值和序列y的均值。r(x,y)∈[−1,1]，当r(x,y)的值趋向于1时，代表2序列间具有良好的正相关性；当r(x,y)的值趋向于–1时，代表2序列间具有良好的负相关性；当r(x,y)的值趋向于0时，代表2序列间不存在相关性。

在同场景时段内，谐波责任保持相对稳定。若该监测周期内含k个同场景时段，对于第a个同场景时段(a=1, 2, ···,k)的谐波责任r_a为

式中：x(a)和y(a)分别为第a个同场景时段下的馈线谐波电流子序列和母线谐波电压子序列。

为更加直观地对比各馈线谐波责任的相对大小，并便于谐波责任奖惩方案的实施，本文对r_a进行归一化处理。在归一化过程中，以所有馈线谐波源在母线共同造成的结果为基准值1，使所有馈线谐波责任相加为1，其中负值代表该馈线并未产生谐波反而被迫吸收了谐波，为谐波责任划分中的受害者，其值为所有谐波源在该馈线共同作用的谐波。

将所有r_a中的正值处理为

式中：b为正值的个数；B_a(i)为第i条馈线在第a个场景下应承担责任的系数（假设该母线连接有N个谐波源馈线，i=1, 2, ···,N）。

归一化完成后，产生谐波的馈线应承担的责任为吸收谐波的馈线被迫承担的责任为其计算过程为

式中：为r_a中的负值，表示该馈线被迫承担的谐波责任；c为负值的个数。

综合考虑各场景下谐波责任对应的时间范畴，累计第i条馈线第1个至第k个同场景时段的谐波责任，馈线i在监测周期内的总谐波责任F_i为

式中：E_a(i)为第i条馈线在第a个场景归一化后的谐波责任值；T(a)为监测周期内第a个同场景时段的时长。此后，运用该归一化方法，对总谐波责任进行归一化处理，得到更加直观和便于比较且考虑时间范畴的各馈线总谐波责任。

04 算例分析

4.1  仿真分析

为验证本文方法的可行性和准确性，在Matlab/Simulink平台上搭建如图9所示的3馈线多谐波源系统，进行谐波责任划分仿真分析。

图9  多谐波源等效电路

Fig.9  Equivalent circuit of multi-harmonic sources

以5次谐波为例，各时段的馈线与母线采样时间差以及系统侧和用户馈线侧仿真参数设置如表1所示，共采集组样本点数据。其中，系统侧5次谐波电压源U_S=50 V，用户侧3条馈线的5次谐波电流源分别为I_c1=7.1 A、I_c2=5.4 A和I_c3=3.5 A。为模拟电网中谐波源的扰动，分别向系统侧谐波电压源和3个馈线谐波电流源加入以初始值为中心，方差和的正态噪声扰动。

表1  仿真参数

Table 1  Simulation parameters

运用本文方法对仿真得到的样本数据进行预处理，时段1中馈线1的部分谐波电流测量数据预处理后效果对比如图10所示。

图10  数据预处理效果

Fig.10  Data preprocessing effect

从图10可见，预处理后数据分布比预处理前更趋于稳定，且离散点和奇异值减少，但其变化趋势和规律仍与原来一致，数据在保留了原有信息的同时噪声显著降低。

数据预处理后，由于仿真设置中母线谐波电压采集和馈线谐波电流采集之间存在非同步时间差，其时序数据之间并非存在一一对应关系。为此，以母线采集到的5次谐波电压为基准曲线，对预处理后的3条馈线时序5次谐波电流数据曲线采用本文方法实现数据匹配对齐。为清晰展示，以时段2馈线1匹配对齐的过程为例，截取部分数据进行分析验证。忽略曲线之间的幅值大小，保留曲线形状和变化趋势，将多条曲线置于同一坐标轴下进行对比分析，其对齐前后的对比如图11所示。

图11  部分数据的匹配对齐

Fig.11  Matching alignment of partial data

从图11中馈线5次谐波电流变化曲线和母线5次谐波电压变化曲线可见，在该场景下两者间具有强相关性，局部位移将导致两者间的相关性发生显著变化。若未对时序数据序列存在的局部位移进行有效处理，而直接基于相关性分析进行谐波责任划分其结果将失去真实性与可靠性。对齐后的馈线5次谐波电流曲线能与母线5次谐波电压曲线实现图11中A点与B点平行于y轴的一一对应，其馈线5次谐波电流曲线对齐前后同一点对应的时间差与仿真设置的条件一致。

将非同步采样下的时序数据序列匹配对齐后，不同场景下的谐波阻抗仿真参数具有显著变化，其变化也将导致母线谐波电压和馈线谐波电流的相关性发生改变。对此，采用OPTICS聚类算法进行聚类分析，以实现不同相关性的场景划分。图12中，横轴表示谐波数据的处理顺序，纵轴表示当前处理对象的可达距离。将可达距离d_r(i)与设置好的领域半径ε进行对比，若d_r(i)＜ε，则该对象可达距离有意义，对应的谐波数据被聚为一类。输出水平线ε与曲线d_r(i)相交处以下的低谷数据有4簇，对应仿真设置的4个场景。

图12  OPTICS聚类的有序队列

Fig.12  Ordered queue of OPTICS clustering

此后，运用本文基于相关性的谐波责任划分方法计算各场景下的谐波责任，并以同一场景下各馈线谐波源在母线造成的电压畸变为基准值进行归一化，各场景下馈线5次谐波责任的计算结果如表2所示。可以看出，同一馈线在不同场景下5次谐波责任具有明显差异，不同场景下各馈线间的5次谐波责任相对大小关系不同。在场景1中，馈线2被迫吸收了谐波，为谐波责任划分中的受害者，与母线共同承担了谐波源造成的后果，但随着场景的变化，该馈线由受害者转变为责任者。

表2  各场景下馈线的谐波责任

Table 2  Harmonic responsibility of the feeder under each scenario

从表1的仿真参数设置中可看出，不同场景对应的时长不同，表2的谐波责任为该场景下对应的责任，但各场景的谐波责任所对应和持续的时长不同，若按以往不考虑对应场景时长进行场景划分的方法把各场景谐波责任简单相加取平均值作为考虑场景划分后的总谐波责任，其结果难以信服。对此，运用本文考虑各场景时间范畴进行谐波责任划分的方法计算考虑场景时长占比后的各场景5次谐波责任如表3所示。

表3  考虑场景时长占比后的谐波责任

Table 3  Harmonic responsibility with consideration of the scenario duration percentage

表3得到的谐波责任相对大小关系代表该馈线在该场景下对监测全周期的总谐波责任贡献度的相对大小。结合表2和表3，以馈线1为例进行分析。馈线1在场景1时谐波责任最大，但在综合考虑各场景相对时长占比和谐波责任大小后，场景3下的谐波责任对总谐波责任结果贡献最大。馈线1在场景1时为短时局部谐波责任最大，场景3虽谐波责任和时长占比都并非最大，但在考虑谐波责任大小和持续时间后，该场景下的谐波责任对馈线1的总谐波责任贡献度最大。

为验证本文方法考虑监测数据非同步性、谐波阻抗变化和场景持续时长的必要性，将本文方法计算得到的总谐波责任与未进行上述考虑（在对比过程中，仅以其中一种作为自变量）的结果对比分析，结果如表4所示。

表4  总谐波责任对比

Table 4  Total harmonic responsibility comparison

从表4可以看出，如未对局部位移造成的谐波数据序列非同步性进行有效处理，而直接基于相关性分析进行谐波责任划分，其计算结果已不再具有工程实际价值。未考虑谐波阻抗变化不进行场景划分的计算结果与进行场景划分并考虑场景时间范畴的计算结果存在一定差异，但其数值大小整体相对接近。这从侧面验证了本文方法的准确性，进行场景划分后谐波责任计算结果的精度有较大提升，能有效应对谐波阻抗变化对谐波责任计算带来的影响。考虑场景划分后各场景持续时长占比的谐波责任划分方法可改进以往只进行场景划分（默认各场景时长相同）带来的不足，有效处理监测全时段内谐波责任的短时剧烈波动，使谐波责任计算结果更加精准。

4.2  实例分析

为验证本文方法在实际工程运用中的有效性，以云南电网公司某110 kV变电站10 kV母线及该母线所包含的3条馈线为例进行谐波责任划分。电能质量监测装置监测点如图4所示，监测周期为该变电站谐波污染较为严重的7天。由于数据体量较大，为清晰展示，综合考虑该变电站监测周期内各次谐波电压含有率、各馈线间不同次谐波电流含有量差距和各次背景谐波波动情况后，本文以7次谐波为例进行谐波责任划分分析。

按本文方法步骤进行7次谐波责任划分，结果如表5所示。表5中各场景谐波责任为未考虑时长占比的馈线谐波责任值，全时段谐波责任为考虑各场景时长占比后的监测全时段总谐波责任值。

表5  实例分析的谐波责任划分

Table 5  Harmonic responsibility division in case study

从表5可以看出，3条馈线的7次谐波责任在监测周期内存在明显波动，3条馈线在不同场景下谐波责任有显著差异且变化趋势不相同。虽谐波责任存在变化，但在各场景中馈线1的谐波责任均为最大，与全时段总谐波责任结论一致。以馈线1为例，虽该馈线在场景2时谐波责任最大，但持续时间不长，综合考虑各场景时长占比后，场景3下的谐波责任对该馈线总谐波责任值贡献最大。若未考虑各场景时长占比，其总谐波责任值将受短时剧烈波动的影响而远离于整体值。馈线1应承担7次谐波主要责任，馈线3承担次要责任，馈线2为谐波责任划分中的受害者。

实际生产运行中很难切断某条馈线来验证本文方法的合理性，经调研和查阅该变电站资料，对3条馈线所接负荷类型和大小分析，与馈线1和馈线3主要产生谐波的分析结果相符。根据该变电站的运行资料，本文方法划分的场景个数及其大致时长占比与该变电站的运行场景情况基本一致。综合上述分析，本文方法所得谐波责任划分结果与电网实际运行情况较为符合。

05 结语

本文提出了一种在考虑电网现有谐波监测数据非同步性和谐波阻抗变化后，基于相关性分析进行谐波责任划分的方法，并通过仿真分析和实例分析验证了该方法的可行性与合理性。

该方法为实际工程应用中谐波责任划分问题提供了新思路和新方法，但在长时间尺度内面临复杂工况时的适用性和可靠性还有待探究。下一步将利用谐波责任划分结果进一步制定合理性经济性指标，并对在线应用进行研究。

注：本文内容呈现略有调整，如需要请查看原文。