随着分布式新能源和电力电子设备的大量接入，电力系统中谐波源呈现大幅增长，谐波污染问题逐步向复杂化、全网化发展，对电网的电能质量带来巨大影响[1]。谐波溯源是谐波治理的前提，其关键在于明晰区域内主导谐波源位置及其影响范围，对改善电能质量、维持电网稳定具有重要意义。

当前的谐波溯源问题可分为定性溯源和定量溯源。定性溯源主要根据谐波的功率方向展开研究，通过判断PCC点处有功功率或无功功率的流向来确定主导谐波源位置[2-5]。而定量溯源研究主要利用实测电流电压数据来量化谐波源在关注节点的谐波贡献度。对于仅有单个谐波源的电力网络，基于诺顿等效电路模型将网络划分为系统侧和用户侧，通过计算两侧的谐波阻抗即可得到谐波源对关注节点谐波畸变的影响程度[6-9]。对于含多谐波源的复杂网络，其关注母线的谐波畸变是所有谐波源交互影响、共同作用下的结果。Xu[10]和Hooman[11]等首次提出了多谐波责任划分量化模型及相应的基于最小二乘的求解方法，该模型通过定义各谐波源单独作用时产生的谐波电压在关注母线测得谐波电压上的投影来描述各谐波源在关注母线上的谐波责任，其求解关键在于找到背景谐波稳定且只有一个谐波源波动的有效数据集。然而，当谐波源数量增加时，这样的有效数据集将大幅减少。后续研究在此基础上提出了一系列改进的多谐波源责任求解方法。张安安等[12]提出了一种新颖的谐波阻抗估计算法，但因其计算量庞大，不利于推广。王杨等[13]提出一种多元回归算法用以估计量化谐波责任。该方法可最大限度降低背景谐波波动的影响，但要求各可疑谐波源在关注母线上产生的谐波电压与关注母线电压间的夹角不变。Wang[14]、徐方维[15]、Zhao[16]等利用复数独立分量分析法计算多谐波源谐波责任，该方法能够更好地抵御背景谐波干扰，但求解过程较为复杂，在需要评估的可疑谐波源过多时，运算量增大且计算准确度下降。此外，该方法要求谐波电压电流幅值相角已知，因此需要测量数据之间具有严格的时钟同步。

尽管谐波溯源问题在理论上取得了一定进展，但在谐波源高密度、分散化、全网化的背景下，上述方法的现实应用都面临以下3个困境：1）当系统可疑谐波源较多时，谐波责任量化算法的精度通常不理想，算法的计算负担较重；2）中国现有电能质量监测平台数据以谐波电压为主，仅有少量关键负荷配有谐波电流监测数据，且数据质量较差，难以满足谐波责任划分需求，具体在第1节中讨论；3）若安装额外的高精度电能质量监测装置或采用便携式电能质量监测装置，需明确主导谐波源位置。现有工程中大多通过经验或者负荷类型判断，误差较大。

面对上述工程难题，利用数据挖掘技术进行谐波污染分区，缩小主导谐波源范围，实现多谐波源的区域化定位，是一种有效且实际的解决思路。相关研究在近年取得一些进展，例如：于浩等[17]提出了一种基于时间序列模式匹配的电能质量污染区域化治理策略，通过动态时间弯曲的趋势–时间距离实现时间序列模式特征矩阵的模式匹配，并计算出节点间电能质量的长期耦合度，最后通过设置阈值划分电能质量治理区域。但在区域划分过程中同一个节点可能被划入多个区域，导致网络中交叉区域过多。并且随着网络节点的增多，该方法的分区难度也相应升高，分区结果更加模糊。石磊磊等[18]提出利用灰色关联理论对预处理后的等时维序列数组进行相关性分析，形成各节点电能质量序列的关联度矩阵，然后通过分层聚类将网络划分为多个治理区域。由于灰色关联度的大小易受参考序列等多种因素的影响[19]，使得节点间的关联度不具有唯一性，导致该方法的聚类输入不明确，影响网络的区域划分结果。此外，上述方法均是依靠人工设定分区阈值，不仅在实际工程中难以实现，而且可能由于阈值过高或过低造成分区结果准确性低、效果差，在大规模电力网络使用中具有较大的局限性。

针对现有方法不足之处，本文深入研究了电能质量时间序列的数据特点，提出一种基于DTW[20]和DBA[21]的聚类算法kDBA++，对中国现有电能质量监测平台提供的谐波电压数据开展数据挖掘工作。该方法以kmeans++算法为骨架，通过DTW距离度量时间序列之间的形状相似性，并采用DBA算法求取聚类质心，进而完成聚类，因此称之为kDBA++。所提方法克服了阈值选取困难、分区结果不明确的缺点，能够适用于大规模电力网络的区域化分析。最后通过IEEE123节点配电系统仿真案例和某地区电网实测数据对所提方法进行验证和对比分析，分区结果证实了该算法在实际工程中的准确性，能够为谐波溯源提供有效的技术支撑。

电能质量在线监测系统已应用于全国27家省（市）公司，覆盖110 kV及以上电压等级的主网回路15.3万条，设立电网电压考核点2.1万个，装载设备256万台[22]。该系统主要由终端数据采集、通信网络传输、监测管理中心3部分组成，位于系统终端的电能质量监测装置以IEEE 1159.3—2019 电能质量数据交换格式（power quality data interchange format，PQDIF）[23]将采集到的基波和谐波数据通过以太网上传至监测管理中心，以实现电能质量的智能化统一管理。监测数据具有以下特点：

1）非实时性。电能质量监测装置每10个周期进行一次谐波分析计算，以3 min为间隔统计900组数据的最大值、平均值和最小值作为电能质量监测数据。在PQDIF文件中，包含的谐波相关数据见表1，一般仅要求测量19次谐波，但谐波次数最高可达50次。

2）非同步性：电能质量监测装置在分析数据的过程中以本地时钟为参考，这就造成了全网各点记录数据的非同步性，即使通过网络校时等手段来尽可能地校正时钟，各节点数据间仍难以完美同步。

电网公司通常从电压层面对全网进行考核，从而监测数据以电压为主，而电流数据较少。现有PQDIF所记录谐波数据虽然包括了谐波电流和电压，但大部分谐波电流数据为零，即无实测值。同时，非同步测量下通过现有电能质量监测平台难以获取谐波责任划分所需的电压电流瞬时同步数据，亟需寻求基于现有数据的谐波溯源方法。在数学上，节点的谐波电压近似由网络内所有谐波源线性组合而成，不同节点受各谐波源的影响不尽相同。因此，受同一种或多种谐波源影响的节点，其谐波电压数据存在一定的关联性特征，在时间序列上体现为节点间数据波动的相似性。但在谐波经线路逐级传导过程中，由于线路阻抗参数等因素的影响，节点谐波电压潜藏的污染特征量各不相同。

利用这一特征，本文通过数据波动相似性进行谐波污染分区研究，分析节点间谐波电压序列的变化规律，挖掘出潜在的耦合特性，对网络各节点进行区域化划分，使得区域内污染关联性强、区域间污染关联性弱。分区的意义：当某分区内发生谐波污染时，只需在该分区内进行谐波溯源，从而大大缩小了溯源范围，为进一步开展谐波源精确定位奠定基础。

本质上，谐波污染分区是一个聚类问题，其对象具有高维度、含噪声、时移性等特点，传统方法求解困难。为此，本文综合考虑了PAA、kmeans++、DTW、DBA各自特点，由此改进提出kDBA++算法，改进思路如图1所示。该算法各部分具体介绍见第2.1～2.5节。

电能质量监测数据常具有高噪声、高维度等特点[24]，直接应用原始数据开展分析将导致结果的精确性不足以及计算效率较低等问题。通过对监测数据进行降噪与压缩等预处理操作，可有效改善上述不足。本文采用PAA[25]对电能质量监测数据进行预处理。将监测数据按时间顺序表示，形成对应的时序v如下：

式中， ${{{v}}_i}$ 表示第i个监测数据，m表示该序列长度。

采用文献[25]中经典PAA算法对监测数据进行压缩降噪，如下所示：

式中， $ \omega $ 为时间窗口宽度， $ v{'_j} $ 为预处理后的第j个监测数据。PAA压缩降噪处理后的数据记为 $ {\boldsymbol{v}}' $ ，如式（3）所示，其长度为n，有 $ n \lt m $ 。

基于PAA算法，某地3次谐波电压幅值时序数据压缩效果如图2所示（时间窗口宽度 $ \omega $ =4），由图2可见，经PAA压缩后，数据可保留原有信息，同时噪声大大降低，且时序维度仅为压缩前1/4，极大降低了后续计算负担。

在各聚类算法中，kmeans算法因其聚类效果优良、计算逻辑简单，得到广泛应用[26]。在kmeans算法基础上，kmeans++算法改善了初始质心的选择策略，从而避免迭代过程中出现质心粘黏，达到提高算法的聚类准确率与计算效率的目标[27]。其初始质心的选择策略如下：

将电能质量监测时序数据作为不同样本，先随机选择一个样本作为初始质心，再计算各样本与该质心的距离，则样本 $ {\boldsymbol{v}}{'_i} $ 被选为下一个初始质心的概率P(i)为：

式中， $ d(i) $ 表示样本i与初始质心的距离。

根据各 $ P(i) $ 生成下一个初始质心，当存在多个初始质心时， $ d(i) $ 表示样本i与最近初始质心的距离。通过式（4）反复迭代，直至选出K个初始质心（K值可通过实际需求指定），从而完成初始迭代。

以上述初始质心为基础，进行kmeans迭代，从而得到最终聚类方案。kmeans迭代过程可通过式（5）和（6）实现，详细流程见文献[28]。

式中， $ {{\boldsymbol{v}}_i} $ 为第i个样本，j*代表离样本点 $ {\boldsymbol{v}}{'_i} $ 最近的质心编号， ${{\boldsymbol{C}}_{{j^*}}}$ 代表包含质心 $ {\bar {\boldsymbol{v}}_{{j^*}}} $ 簇，Dist(a,b)代表向量a和b之间的距离，K代表质心个数。

然而，非同步采样下不同监测节点的电能质量时序存在时移现象，使得传统kmeans++所采用的欧式距离度量函数无法量化数据间的相似性，导致该算法难以直接应用于谐波污染分区。为此，本文以DTW作为距离度量函数，对kmeans++算法进行改进。

DTW 距离度量函数并非计算同时刻数据的差异，而是对时序间数据点进行非同时刻映射，此时各时序点间对应关系如图3所示。由图3可见，相比于传统欧式距离，该对应方式可准确度量存在时移但局部相似的时序间距离。由于该对应关系使不同时刻时序点得以匹配，故称其实现了时间轴“错位”[29]。

为计算长度分别为 $ {n_1} $ 和 $ {n_2} $ 的时序 $ {\boldsymbol{v}}{'_1} $ 和 $ {\boldsymbol{v}}{'_2} $ 的DTW距离，首先需构建 $ {n_1} \times {n_2} $ 维的矩阵 ${{\boldsymbol{D}}_{{\rm{base}}}}$ ，矩阵中元素 $ {d_{{\text{base}},ij}} $ 如式（7）所示：

式中， $ {\boldsymbol{v}}{'_{1,j}} $ 和 $ {\boldsymbol{v}}{'_{2,j}} $ 分别代表时序 $ {{\boldsymbol{v}}'_1} $ 、 $ {{\boldsymbol{v}}'_2} $ 的第i、j个数据， ${d_{{\rm{base}}(i,j)}}$ 表示 $ {\boldsymbol{v}}{'_{1,j}} $ 和 $ {\boldsymbol{v}}{'_{2,j}} $ 之间的特征距离。

从 ${d_{{\rm{base}}}}(1,1)$ 到 ${{\boldsymbol{D}}_{{\rm{base}}}}$ 上任意点都对应一个最小累加特征距离，遍历 ${{\boldsymbol{D}}_{{\rm{base}}}}$ 上所有点，可得最小累加特征距离矩阵 $ {\boldsymbol{Q}} $ ，具体计算过程如下[29]：

首先，根据式（8）对矩阵 $ {\boldsymbol{Q}} $ 中元素q进行初始化：

对于 $ {\boldsymbol{Q}} $ 中剩余元素有：

式（8）～（9）中，i、j均为整数，且有 $ i \in [1,{n_1}] $ ， $ j \in [1,{n_2}] $ 。

在通过上述方法建立的矩阵 $ {\boldsymbol{Q}} $ 中，元素 $ q( {n_1} + 1, {n_2} + 1) $ 即为时序 $ {\boldsymbol{v}}{'_1} $ 、 $ {\boldsymbol{v}}{'_2} $ 间的DTW距离。

然而，将DTW引入kmeans++算法后，质心的求解将成为一个新的难题。通常，质心位于离簇内所有样本距离平方之和的最小的位置，因此DTW距离下质心定义如式（10）所示：

式中，h为待求向量， $ {\boldsymbol{v}}{'_i} $ 为第i个时序向量， ${\rm{DTW}}({\boldsymbol{v}}{'_i},{\boldsymbol{h}})$ 为向量 $ {\boldsymbol{v}}{'_i} $ 与h之间的DTW距离， $ \overline {\boldsymbol{T}} $ 代表簇 $[{\boldsymbol{v}}{'_1}, {\boldsymbol{v}}{'_2}, \cdots ,{\boldsymbol{v}}{'_{{p_1}}}]$ 在DTW距离下的质心。

依此定义，DTW距离不同于欧氏距离，其质心求解问题具有非线性、离散化等特点[30]，传统方法求解精度不佳，为此，本文采用DBA克服这一难题[21]。

DBA的主要思想是，利用DTW找到两序列间的最小距离，并将序列坐标更新为与质心相关联的坐标，从而将复杂的DTW距离度量转换为欧式距离形式。DTW质心的求解过程如下[21]：

式中， $\overline T(l)$ 为式（10）中质心 $\overline {\boldsymbol{T}}$ 的第l个元素， ${\boldsymbol{V}}'$ 为质心 $\overline {\boldsymbol{T}}$ 所在簇，SS为质心 $\overline {\boldsymbol{T}}$ 与簇内其他时序向量的距离平方和， ${{{M}}_i}$ 为更新后的序列坐标空间， $ e' $ 为更新后的时间序列，L为时间序列长度。

鉴于式（11）无最大值，因此可以通过梯度下降法求得其极小值，即：

式（12）化简得式（13）：

对（13）求解得：

按照式（11）～（14）所述方法即可得到DTW距离下质心 $\overline {\boldsymbol{T}}$ ，据此可得到DBA算法流程如下：

1）随机设置簇内任一向量 $ {\boldsymbol{y}}' $ 迭代向量，其余向量为参考向量。设t为参考向量长度，并设置t个集合 $ {g_1},{g_2},\cdots,{g_t} $ ，且均初始化为空集；

2）设 ${{\boldsymbol{D}}_{{{{\rm{base}},yr}}}}$ 为 $ {\boldsymbol{y}}' $ 与某一参考向量 $ {\boldsymbol{v}}{'_r} $ 对应特征距离矩阵， $ {\boldsymbol{W}} $ 为 ${{\boldsymbol{D}}_{{m{{\rm{base}},yr}}}}$ 从始点到末点最小累加特征距离路径的坐标矩阵，每一列对应路径上某一点的坐标。设k为 $ {\boldsymbol{W}} $ 列数，则矩阵 $ {\boldsymbol{W}} $ 维度为 $ 2 \times k $ 。并初始化迭代次数x=0；

3）将 $ {\boldsymbol{W}} $ 中第(k–x)列元素记为[i,j]T（上标T为转置），执行 $ {g_i} = {g_i} \cup {\boldsymbol{v}}{'_r}(j) $ ，x=x+1；

4）重复步骤3），当 $ x = k $ 时，停止迭代，并更换参考向量；

5）重复步骤2）～4），直至遍历簇中所有参考向量；

6）对每个 $ g $ 集合里元素求均值，得到 $ \bar g $ ，并构成向量 $ {\boldsymbol{f}}=[{\bar g_1},{\bar g_2}, \cdots ,{\bar g_t}] $ ，随后将所有集合g重置为空集；

7）以f为新的迭代向量，并重复步骤2）～6），直至两次得到的f差距小于某一阈值，即认为f为DTW距离下的质心。

对于具有时移特性的时间序列，采用欧氏距离法和DBA求解得质心如图4所示，图4中，78条曲线代表78个节点谐波电压幅值时序数据。算数平均法求得的质心波形与各样本差异较大且严重失真。相比之下，DBA求得的质心与样本更相似且波形未畸变。从而可知，相比于传统方法，DBA可有效求解时移序列的质心。

综上，本文基于海量具有时移特性的电能质量监测数据，采用kDBA++算法开展谐波污染分区工作：在传统kmeans++中引入DTW算法改进距离度量函数，并通过DBA算法求解DTW距离下的质心，从而完成监测数据与相应节点的聚类，且每一聚类对应一个谐波污染分区。然而，在实际工程应用中，往往缺乏确定K值的先验知识，从而需要引入聚类效果评估指标，以寻找最优聚类数目K。

常用的聚类评价指标有I指标、COP指标、Silhouette系数。其中：I指标过于依赖用户参数设置，评价结果具有主观性[31]；COP指标在数据规模增大时，时间复杂度急剧上升，不适合评价大规模电能质量数据聚类结果[31]。相比之下，Silhouett系数不依赖用户主观参数选取，时间复杂度不高，因此本文采用该指标评价聚类结果，其计算过程如下[32]：

1）对于任意样本 $ {\boldsymbol{v}}{'_i} $ ，计算其和同簇内剩余样本距离的均值，记作 $ a\left( i \right) $ ，用于量化簇内的凝聚度。

式中， $ {\boldsymbol{v}}{'_j} $ 为簇内除 ${\boldsymbol{v}}{'_i} $ 的其他向量，p2为同簇样本数量。

2）对于同一样本 $ {\boldsymbol{v}}{'_i} $ ，计算其和最近簇中所有元素距离的均值，记作 $ b\left( i \right) $ ，用于量化簇间的相异度。

式中， $ {\boldsymbol{u}}{'_j} $ 为最近簇中第j个样本，q为该簇所含样本个数。

3）样本 $ {\boldsymbol{v}}{'_i} $ 的Silhouette系数为：

4）计算所有样本点的Silhouette系数，其均值即为当前聚类结果的Silhouette系数，记为指标S。

Silhouette系数的取值在[–1,1]之间，为保障聚类效果，在实际应用中应使聚类结果的Silhouette系数大于0.8。

考虑到实际工程中难以通过先验知识确定谐波分区个数，本文利用Silhouette系数自适应选择K值。将K值从2开始逐步增大至Kmax，在每个K值下进行聚类，并采用Silhouette系数评价各次聚类结果，使Silhouette系数最大的K值即为最优聚类数目。Kmax通常根据式(18)确定[33]。

式中， $ \left\lceil . \right\rceil $ 表示向下取整，c为样本总量。

采用Silhouette系数和kDBA++算法开展谐波污染分区工作的整体流程如图5所示，具体步骤如下：

1）采用PAA算法预处理每个节点的电能质量时序数据，完成降维降噪工作。

2）将K值从2开始逐步增大至Kmax。

3）以K为聚类数目，通过kDBA++算法对数据进行聚类。

4）采用Silhouette系数评价各K值下的聚类结果。

5）找出使Silhouette系数最大的K值，相应的聚类结果即为谐波污染最优分区。

通过上述步骤，对谐波污染进行区域化划分，当某一分区谐波污染严重时，仅需在该分区内开展谐波溯源，即可为进一步精确定位主导谐波源奠定基础。

仿真以图6所示IEEE123系统为参考模型，该系统包含123个节点、2台变压器、4个无功补偿电容器。在系统中随机选取8个节点[84,96,114,6,20,24,39,9]注入3次谐波电流。其中：前4个节点注入较大谐波，作为主导谐波源节点；后4个节点注入较小谐波。各谐波源分布如图6所示。

主导谐波源节点注入电流的最大值为0.56 A，考虑谐波源的投切，谐波电流趋势存在一定突变。后4个节点注入谐波幅值为均匀分布的随机变量，作非主导谐波源，设置其均值为0.05 A，极差为0.1 A。上述谐波源时序曲线如图7所示。其中，左边4条曲线分别对应4个主导谐波源节点，右边4条曲线对应4个非主导谐波源节点。

仿真数据颗粒度与实际电能质量监测数据相同（1样本/3 min），模拟时长72 h，共计数据1 440组，如图8所示，图8中132条曲线对应系统内132个节点。下面将从聚类结果的波形相似度及Silhouette系数两方面对比分析各算法的谐波污染分区效果。

采用kDBA++算法所得谐波污染分区如图9所示，可见，系统被分为A、B、C、D 4个区域，且各区域只含有一个主导谐波源。聚类结果的Silhouette系数评价为0.823 5，满足实际工程要求。

采用kDBA++算法求得各区域节点谐波电压幅值波形，如图10所示。

图10中132条曲线对应系统内132个节点。各区域内节点谐波电压波形相似，不同区域间波形差异大，符合聚类原则。

文献[16]所提灰色关联度层次聚类法分区结果如图11所示。

图11中132条曲线对应系统内132个节点，虽然系统也被分为A、B、C、D 等4个区域，但A、B区域内节点的谐波电压波形混乱，聚类结果的Silhouette系数为0.454 4，低于kDBA++算法。

采用文献[17]中所提出的阈值聚类法求得聚类数目共34类。其中具有代表性的分类如图12所示。

图12 中4幅子图中共计245条曲线，分别代表245个节点。该数目远超真实样本数量（132），表明结果中一点多区现象严重。且由图12可见，每一分区内均波形杂乱。同时，Silhouette系数仅为0.232 9，难以满足实际工程需求，因此可认为该算法实际效果欠佳。

采用经典聚类方法kmeans++分区结果如图13所示。图13中132条曲线对应系统内132个节点，由于其难以度量时移时序间距离，因此所得聚类结果不理想，虽然B区域内时序曲线变化趋势相似，但 A、C、D区内部波形杂乱，差异较大，且该算法聚类结果Silhouette系数为0.537 3，难以满足工程需求。因此可认为kmeans++算法在处理时移时序聚类问题上效果欠佳。

为进一步验证kDBA++算法求得谐波污染分区的正确性，以分区B为例，采用文献[34]所提谐波贡献度计算方法，分析8个谐波源对B区中11个非谐波源节点的谐波电压贡献度如图14所示。图14中，每一柱形均代表了不同谐波源节点对B区内部非谐波源节点的谐波电压贡献度。可见，B区内主导节点96对该区所有非谐波源节点的谐波贡献度最大，相比之下，其他谐波源贡献度明显较小，从而再次论证了分区结果的准确性。其他3区情况与B区类似。因此，当网络中某一母线3次谐波电压含量超标时，只需在该母线所在区域内搜索可疑谐波源，就可大大缩小主导谐波源追溯范围。

值得注意的是，谐波污染分区的意义在于，同一分区内的谐波畸变主要由该分区内的谐波源导致，而受其他分区谐波源影响较小。从而本质上，同一谐波污染分区中各节点间谐波转移阻抗较大，易受该分区内谐波源影响，不同分区之间各节点谐波转移阻抗较小，因此谐波源对非同一分区的节点谐波畸变影响较小。图9 C区中之所以包含大量节点，正是因为这些节点之间的谐波转移阻抗较大，从而被聚为一类。为进一步说明不同分区中各节点之间转移阻抗大小关系，将C、B两分区之间各节点转移阻抗呈现出来，见表2。由表2可见，C区、B区内各节点之间谐波转移阻抗较大，而这两个分区之间各节点谐波转移阻抗较小，再次论证了分区结果的正确性。

考虑到实际电网中存在其他频次谐波，故基于5次谐波污染仿真数据，再次进行验证。该频次谐波污染同样由[84,96,114,6,20, 24,39,9]8个谐波源引起，前4个谐波源电流幅值较大，起主导作用，后4个电流幅值较小，以模拟背景谐波。5次谐波电压数据共计1 440组，如图15所示。图15中132条曲线对应系统内132个节点。

针对该数据，kDBA++聚类效果如图16所示。图16中132条曲线对应系统中132个节点。

同样地，kDBA++将谐波污染分为4个区域，区域内谐波电压关联度较强，区域间较弱，聚类效果Silhouette系数为0.884 8，满足工程需求（Silhouette系数应大于0.8）。且每一区域均包含有一个主导谐波源。本次仿真验证充分表明在多次谐波污染叠加的情况下，kDBA++算法仍然能开展谐波污染分区任务。

各方法计算耗时为秒级，满足工程需求，因此本文未对耗时问题做深入对比。

采用中国某地区220 kV输电网的65个节点实测数据对所提方法做进一步验证。以图17所示5次谐波电压为例，65条曲线对应65个不同节点，其中：地区A数据36组，编号为1～36，地区B数据29组，编号为37～65。调研发现，两地区电气耦合很弱，且每个地区均有一个大型工业负荷作为主导谐波源，该先验知识可用于算例结果的验证。

通过kDBA++和Silhouette系数分析，最优聚类数目为K=2，聚类效果如图18所示。图18中65条曲线代表系统内65个节点。此时Silhouette系数S=0.913 2，满足工程应用要求。由图18可见，数据被清晰地分为了A、B两个区域，且各区域内电压波形相似，区域间波形差异较大。此外，A分区内节点编号为1～36，B分区内节点编号为37～65，与地理位置关系完全一致，因此可认为kDBA++完成了谐波污染分区的任务。

基于轮廓系数和波形将kDBA++的聚类结果与文献[16]中灰色关联度层次聚类法和文献[17]中阈值聚类法进行对比。首先基于灰色关联度层次聚类法对实测数据进行聚类分析，数据被分为A、B两区，如图19所示，65条曲线对应系统内65个节点。可见，尽管A区中的电压波形相似，但B区波形杂乱，聚类精度低，且该方法聚类结果Silhouette系数仅有0.657 5，低于工程需求（Silhouette系数应大于0.8）。

随后采用文献[17]中基于形状相似性的阈值聚类法开展谐波污染分区，其结果共计14类，远超据前述案例背景介绍中先验知识所估计的分区数量（2类）。从这14类结果中选取具有代表性的2类结果呈现如图20所示，图20（a）、（b）中共计88条曲线分别对应88个节点，曲线数量远大于实际节点数目，一点多区现象严重。且由图20可见，该方法将不同趋势的曲线归为同一类，聚类结果不可靠。此外，聚类结果的Silhouette系数为0.110 2， 低于工程应用要求。因此，可认为该方法难以有效实现谐波污染分区。

本文提出了一种基于kDBA++聚类算法的谐波污染分区溯源策略并利用仿真和实测数据进行了详细验证。针对中国电能质量数据非同步特点，以kmeans++算法为骨架，通过DTW距离度量时间序列之间的形状相似性，并采用DBA算法求取聚类质心，提出了新型kDBA++聚类算法。进一步结合Silhouette系数实现聚类效果评价并确定最优分区数目，提升工程应用可靠性。最后，采用IEEE123系统的仿真数据和中国某地实测数据对kDBA++和现有方法进行了比较，结果表明，理论上kDBA++无论是准确度还是在抗非主导谐波源干扰能力都具有一定优越性，更适用于复杂大电网下的谐波分区溯源工作。

所提方法已实现将主导谐波源搜索范围缩小到区域级，如何通过海量电能质量监测数据并结合精细化聚类技术进一步将搜索范围缩小到母线级，最终锁定主导谐波源，仍需进一步研究。