在当前能源变革背景下，大量分布式能源与可控负荷接入电网^[1-3]，电力系统的“电力电子化”程度愈加显著。电力电子技术的广泛应用有效地提升了电能的转换和传输效率，但同时也向电力系统引入了大量谐波^[4]。由谐波导致的电能质量问题威胁电网的安全、可靠和稳定运行，给电网的优质供电带来挑战。

针对谐波污染亟需治理的现状，国际上提出了谐波治理的“奖惩性方案”^[5]，而谐波源定位与责任划分是经济奖惩措施有效实施的前提^[6-11]。谐波源辨识作为谐波污染治理的基础，其研究开展相对较早^[12-13]。谐波源责任划分方法主要有线性回归法^[14-15]，波动量法^[16]，概率类方法^[17]和盲源分离法^[18-19]等。当前已开展的研究工作多集中于评估公共耦合点(point of common coupling, PCC)系统侧和用户侧的谐波责任。然而，实际电力系统中PCC处的谐波畸变通常由多个谐波源共同作用产生^[20]，除了计算单个谐波源对PCC的谐波贡献，还应考虑谐波源之间的相互影响。特别是在现代电力系统中，随着可再生能源和电动汽车的广泛接入，还需考虑不同谐波源的时变特性^[21]。由于在实际系统中，运行状态以及设备、负荷的时变性，包括准确获取系统参数，都存在一定难度，而直接基于实测历史数据评估谐波责任具有较好的应用前景^[22-23]。

本文提出了基于数据统计相关性分析的多谐波源责任评估方法。基于PCC点谐波畸变电压与谐波源电流之间的相关规律，提出利用偏相关分析法和多元线性回归分析谐波责任，选择出主要谐波源，基于偏最小二乘回归法建立了评价多谐波源系统各谐波源的责任指标；为选择满足分析要求的实测数据段，提出利用负荷水平分割法将谐波超标节点的谐波电压和谐波电流数据分段；结合统计相关性分析衡量谐波源的谐波电流和超标节点的谐波电压之间的统计规律，计算出谐波指标，确定谐波源的谐波责任定量大小。

1 谐波数据相关性分析 1.1 谐波电压与谐波电流的线性相关性

系统中PCC节点的谐波电压幅值与谐波源的谐波电流幅值之间存在线性相关性^[24]。以某系统的实测数据为例，谐波电压超标节点X的5次谐波电压幅值( $V_5^X$ )与接入的谐波源负荷G的5次谐波电流幅值( ${I_{5A}}$ )之间的线性相关性如图1所示。

根据图1中谐波电压和电流的关系，利用线性回归法得到：

$V_5^X = 0.26I{}_{5A} + 4.82$

(1)

利用式(1)的回归系数对谐波责任进行评估，从而获得适用于各谐波源相互独立的系统^[11]。但在实际电力系统中，往往存在多个分散式谐波源^[16]，且谐波源的谐波电流间存在相关关系，利用式(1)评估各谐波电流对PCC谐波电压的贡献大小存在难度。

1.2 谐波源偏相关性分析

在实际系统中，谐波源数目的增多带来复杂的谐波数据选择问题。针对某一特定的谐波问题，选择出贡献较大的3~4个谐波源，即可有针对性的进行谐波治理，提高评估精度。基于前述谐波电压和谐波电流数据之间的线性相关性，提出基于偏相关分析^[18]的数据预处理。

假设电力系统中含有n个谐波源，即有n个相关变量x₁，x₂，…，x_n，取m组观测数据，其n–2级偏相关系数的计算过程如下。首先计算直线相关系数 ${r_{ij}}$ ：

$ {r_{ij}} \!=\! \frac{{\displaystyle\sum {({x_i} - {{\bar x}_i})} ({x_j} - {{\bar x}_j})}}{{\sqrt {{{\displaystyle\sum {\!(\!{x_i} - {{\bar x}_i}\!)\!} }^2} \!\cdot\! {{\displaystyle\sum {\!(\!{x_j} - {{\bar x}_j}\!)\!} }^2}} }}\left( {i \!=\! 1,2, \cdots, n;\;j \!=\! 1,2, \cdots, n} \right) $

由简单相关系数 ${r_{ij}}$ 组成相关系数矩阵R，然后求R的逆矩阵C：

${{C}} = {{{R}}^{ - 1}} = \left[ {\begin{array}{*{20}{c}} {{c_{11}}}&{{c_{12}}}& \cdots &{{c_{1n}}} \\ {{c_{21}}}&{{c_{22}}}& \cdots &{{c_{2n}}} \\ \vdots & \vdots & & \vdots \\ {{c_{n1}}}&{{c_{n2}}}& \cdots &{{c_{nn}}} \end{array}} \right]$

则相关变量x₁，x₂，…，x_n之间的偏相关系数为

$ {p_{ij}} = \frac{{ - {c_{ij}}}}{{\sqrt {{c_{ii}}{c_{jj}}} }}\;\;\left( {i = 1,2, \cdots, m;j = 1,2, \cdots, m;i \ne j} \right) $

(2)

式(2)中的偏相关系数越大，说明变量间相关性越强，谐波源间的耦合和相互干扰越强。利用偏相关分析，对实测数据进行预处理，找出对观测节点谐波电压影响较大的主要谐波源。基于对PCC点和谐波源、谐波源与谐波源间的谐波电压和谐波电流数据的相关性分析，可进一步利用偏最小二乘回归划分谐波责任。

2 基于数据统计相关性的谐波责任划分

假设系统中节点X(观测节点)的谐波问题由n个谐波源共同作用产生。各谐波源的h次谐波电流分别表示为 ${{{I}}_{h1}},{{{I}}_{h2}}, \cdots ,{{{I}}_{hi}}, \cdots ,{{{I}}_{hn}}$ (h为谐波次数，i为谐波源的馈线编号)，节点X处的h次谐波电压为 ${{V}}_h^X$ ：

${{V}}_h^X = \underbrace {{{Z}}_{hi}^X{{{I}}_{hi}}}_{{{V}}_{hi}^X} + \underbrace {{{Z}}_{h1}^X{{{I}}_{h1}} + {{Z}}_{h2}^X{{{I}}_{h2}} + \cdots + {{Z}}_{hn}^X{{{I}}_{hn}} + {{{V}}_{h0}}}_{{{E}}_h^X}$

式中： ${{Z}}_{hi}^X$ 为观测节点X与谐波源节点i之间的h次谐波阻抗， ${{{V}}_{h0}}$ 为背景谐波。

X点的h次谐波电压由2部分组成： ${{V}}_{hi}^X$ 来自谐波源i， ${{E}}_h^X$ 来自其他谐波源。其中衡量谐波源i在X点处产生的谐波责任指标，由h次谐波电压 ${{V}}_{hi}^X$ 在X点的h次谐波总畸变电压 ${{V}}_h^X$ 上的投影表示，如图2所示，其中 $\alpha $ 为 ${{V}}_{hi}^X$ 和 ${{V}}_h^X$ 的夹角。

若量化单个谐波源的责任贡献，谐波源i对X点的h次谐波电压的谐波贡献率 $H_i^X$ 为

$H_i^X = \frac{{\left| {{{V}}_{hi}^X} \right|}}{{\left| {{{V}}_h^X} \right|}}\cos \alpha \times 100\text{%} $

(3)

式(3)可结合谐波潮流法计算得出，作为谐波贡献划分的依据。整理式(1)得到：

$ \begin{array}{l} \left| {{{V}}_h^X} \right| = \underbrace {\left| {{{Z}}_{h1}^X} \right|\left| {{{{I}}_{hi}}} \right|\cos {\alpha _1}}_{{{V}}_{h1}^X} + \underbrace {\left| {{{Z}}_{h2}^X} \right|\left| {{{{I}}_{h2}}} \right|\cos {\alpha _2}}_{{{V}}_{h2}^X} + \cdots \\ \qquad \underbrace {\left| {{{Z}}_{hn}^X} \right|\left| {{{{I}}_{hn}}} \right|\cos {\alpha _n}}_{{{V}}_{hn}^X} + \underbrace {\left| {{{{V}}_{h0}}} \right|\cos {\alpha _0}}_{{{V}}_{h0}^X} \end{array} $

(4)

式中： ${\alpha _i}\left( {i = 1,2, \cdots ,n} \right)$ 为 ${{V}}_{hi}^X$ 和 ${{V}}_h^X$ 的夹角， ${\alpha _0}$ 为背景谐波 ${{V}}_{h0}^X$ 和 ${{V}}_h^X$ 的夹角。

各相量之间的关系如图3所示。

式(4)为多元线性方程，自变量由各谐波源的对应次数谐波电流矩阵表示，因变量由观测节点的对应次数谐波电压矩阵表示。利用偏最小二乘法回归估计得到式(5)中的系数 ${\lambda _1},{\lambda _2}, \cdots ,{\lambda _n}$ 和常数项 ${\lambda _0}$ 。将式(4)中的各成分简化表示为

$\left\{ \begin{array}{l} {\lambda _0} = \left| {{{{V}}_{h0}}} \right|\cos {\alpha _0} \\ {\lambda _1} = \left| {{{Z}}_{h1}^X} \right|\cos {\alpha _1} \\ {\lambda _2} = \left| {{{Z}}_{h2}^X} \right|\cos {\alpha _2} \\ \begin{array}{*{20}{c}} {}&{}&{} \end{array} \vdots \\ {\lambda _n} = \left| {{{Z}}_{hn}^X} \right|\cos {\alpha _n} \\ \end{array} \right.$

(5)

对一段时间内测量的一组数据进行回归分析，自变量矩阵X和因变量矩阵Y分别为

${{X}} = \left[ {\begin{array}{*{20}{c}} {\left| {{{I}}_{h1}^X({t_1})} \right|}&{\left| {{{I}}_{h2}^X({t_1})} \right|}& \cdots &{\left| {{{I}}_{hn}^X({t_1})} \right|}\\ {\left| {{{I}}_{h1}^X({t_2})} \right|}&{\left| {{{I}}_{h2}^X({t_2})} \right|}& \cdots &{\left| {{{I}}_{hn}^X({t_2})} \right|}\\ \vdots & \vdots & & \vdots \\ {\left| {{{I}}_{h1}^X({t_m})} \right|}&{\left| {{{I}}_{h2}^X({t_m})} \right|}& \cdots &{\left| {{{I}}_{hn}^X({t_m})} \right|} \end{array}} \right]$

${{Y}} = {\left[ {{{\left| {{{V}}_h^X({t_1})} \right|}^2}\quad{{\left| {{{V}}_h^X({t_2})} \right|}^2}\quad...\quad{{\left| {{{V}}_h^X({t_m})} \right|}^2}} \right]^{\rm{T}}}$

式中： $\left| {{{V}}_h^X({t_k})} \right|$ 为各个采样时刻节点X的谐波电压的幅值， $\left| {{{I}}_{hi}^X({t_k})} \right|$ (k=1,2,…,m)为节点i的h次谐波电流的幅值。

回归得到各谐波电流幅值的系数估计值，分别为 $\;{\beta _1},{\beta _2}, \cdots ,{\beta _n}$ 和 $\;{\beta _0}$ ，估计误差 $\xi $ 。即谐波源i(i=1,2,…,n)对X点谐波电压的贡献 $\left| {{{V}}_{hi}^X} \right|$ 可以通过谐波电流幅值 $\left| {{{I}}_{hi}^X} \right|$ 与其系数 $\;{\beta _i}$ 的乘积表示。因此，回归计算后的X点的谐波电压表示为

$\left| {{{V}}_h^X} \right| = {\beta _0} + {\beta _1}\left| {{{I}}_{h1}^X} \right| + {\beta _2}\left| {{{I}}_{h2}^X} \right| + \cdots + {\beta _n}\left| {{{I}}_{hn}^X} \right| + \xi $

谐波源i对节点X的谐波责任为

$F_i^X = \frac{{{\beta _i}\left| {{{{I}}_{hi}}} \right|}}{{\left| {{{V}}_h^X} \right|}} \times 100\text{%} $

(6)

3 统计数据段的选择

在前述谐波责任评估指标计算中，假设随 $\left| {{{V}}_{hi}^X} \right|$ 的增加， $\left| {{{E}}_h^X} \right|$ 不变。若能找到一个时间段，保证在该时段内仅负荷i变化，而其他谐波源基本保持恒定，则在这个时间段内，即可认为 $\left| {{{E}}_h^X} \right|$ 为常量。考虑到实际系统中，满足要求的时间段可能会比较短，而且随谐波源个数的增加，满足要求的时间段在长度和数量上都会明显减少，因此提出基于负荷水平等级分割的数据选择方法，将负荷变化范围分成不同的水平等级。例如，如果谐波源i的谐波电流波动在0~10%，则可认为该数据段内的所有数据处于同一负荷水平，即认为负荷i基本保持不变。实际应用时可根据精度需要确定每一负荷水平的波动范围。

负荷水平等级分割法将所选择出的数据按等级存储。以电力系统中含有3个谐波源为例，将谐波源分别编号为A、B、C，图4给出了该情况下的数据存储结构。谐波源A、B、C的负荷水平都处于第一等级的所有时间段内的数据存储在单元格D₁₁₁内。同理，谐波源B的负荷水平在第四等级，A和C的负荷水平处于第一等级的所有时间段的数据存储于D₁₄₁内。因此图4中的数据存储结构的每一列数据均可用来进行相关性分析。例如，图4中灰色单元格内的数据可用来评估谐波源A的谐波责任。在这些数据中，谐波源B和C基本是保持不变的，只有谐波源A负荷在变动。对谐波源B、C的任意负荷水平等级，都存在一列数据，可用来评估谐波源A的谐波责任。实际应用时，选择数据最多的一列进行评估。

4 仿真分析 4.1 谐波责任评估流程

要评估多谐波源系统中的各谐波源贡献，首先要获取谐波数据，流程概括如下。

1)计算基波潮流。

将所评估的谐波源视作有功功率P和无功功率Q。已知的P、Q负荷，利用牛–拉法计算基波潮流，得到各负荷的基波电流 ${I_{{\rm{rated}}}}$ 。

2)计算谐波电流。

基于恒流源法，选择要研究的谐波次数，根据式(7)推算谐波源在该次谐波的注入电流：

$\left\{ \begin{array}{l} {I_h} = {I_{{\rm{rated}}}}\dfrac{{{I_{h{\rm{ - spectrum}}}}}}{{{I_{1{\rm{ - spectrum}}}}}} \\ {\theta _h} = {\theta _{h{\rm{ - spectrum}}}} + h({\theta _1} - {\theta _{{\rm{1 - spectrum}}}}) \\ \end{array} \right.$

(7)

式中： ${I_{h{\rm{ - spectrum}}}}$ 为谐波电流典型频谱^[1]中h次谐波电流的幅值(以基波幅值 ${I_{1{\rm{ - spectrum}}}}$ 为参考点)， ${\theta _{h{\rm{ - spectrum}}}}$ 为h次谐波电流的相角(以基波相角 ${\theta _{1{\rm{ - spectrum}}}}$ 为参考点)。

3)计算节点导纳矩阵。

基于电力系统网络的支路连接情况和各元件的谐波参数，计算各次谐波对应的节点导纳矩阵 ${{{Y}}_h}$ 。

4)计算谐波潮流。

基于网络节点导纳阵，计算系统的各次谐波潮流，得到观测节点的谐波电压。

5)计算各谐波电压。

根据式(4)把各次谐波电压表示成各谐波源谐波电流的组合。

完成数据采集后，对数据进行预处理，利用负荷水平等级分割法选择出满足要求的数据段。利用偏相关分析判断谐波源间的耦合关系，选定主要谐波源；利用谐波责任指标评估各谐波源的谐波责任，评估整体流程如图5所示。

4.2 算例分析

为验证本文所提算法的正确性，选用IEEE 37节点系统在Matlab/Simulink中仿真和计算。IEEE 37节点网络拓扑图如图6所示。该系统有3个发电机，1个平衡节点，1个PV节点及7个PQ节点。接入分布式电源，研究这些节点对其他节点的谐波影响，谐波源接入位置及其容量如表1所示。设置所有PQ节点的有功、无功在±110%额定范围内随机变化。仿真进行了5 000次计算，采集了该系统中10个节点的谐波电流和谐波电压，生成了谐波责任评估所用的基础数据。在此基础上对谐波源在观测节点的谐波责任进行了评估。

取节点713为观测节点，以5、7、11、13次谐波为例，研究分布式电源和负荷谐波源对该节点的谐波影响。利用负荷水平等级分割法对谐波数据进行预处理，并选取了单元格D₁₁₁内的数据段进行偏相关分析。将各测量节点的5次谐波数据与节点713的5次谐波数据的直线相关系数 ${r_{ij}}$ 和偏相关系数 ${p_{ij}}$ 计算结果整理如表2所示。

由表2数据可知，节点705、735、741与节点713的5次谐波存在高度相关性，因此可以确定这3个节点接入了对观测节点5次谐波电压影响很大的谐波源。该判断与表1中谐波源的实际接入情况一致。同理，对各测量节点与观测节点713的7、11、13次谐波数据进行偏相关性分析。

利用式(6)的谐波责任指标，分别针对不同次数的谐波，通过偏最小二乘法计算各谐波源的责任。各谐波源谐波责任的评估值与准确值的对比结果如表3—表5所示。例如，谐波源705在11次谐波下对节点713的谐波责任为负，表明在该次谐波下，谐波源705削弱了节点713的电压畸变。根据表3—表5中数据对比显示，节点735接入了系统中的主要谐波源。

对比本文所提方法与准确谐波责任值，绘制各谐波源的谐波责任对比如图7所示。从仿真结果可得，该方法对各谐波源的谐波责任评估值均能较好地与准确值吻合，能有效评估各谐波源的谐波责任。同时，图7反映了在各次谐波下，所提方法对主要谐波源的谐波责任评估结果准确度更高。

5 实测数据分析

实测数据来自某变电站35 kV母线，该变电站所供负荷含有较多的谐波源。所测母线接有5条专供馈线。同步采样母线电压及各条馈线的电流，每隔3 s记录6个周波数据，每个周波记录128个数据点。连续采集当日24 h内数据，并利用傅里叶分析得到采集数据的各次谐波值。以污染较为严重的3次谐波为例，先利用负荷水平等级分割法预处理，选取实测数据中适合分析的数据段，然后利用偏相关分析法分析3个主要谐波源A、B、C，分别对应于化工厂专供线、电气化铁路专供线和飞机场专供线。在某段选取分析的数据段内，由1 200个连续采样点采集观测母线的谐波电压和3条主要谐波源馈线的谐波电流数据，绘制如图8所示波形。该数据段内，由偏最小二乘法回归评估结果如表6所示。

在观察谐波责任动态变化时，设置统计数据段选择时的负荷波动范围为6%。先将24 h的谐波数据归类到如图4所示的存储结构中，然后在各数据段内计算偏相关系数，综合定位谐波主要来源。以2 h为时间尺度，即每2 h进行一次各条馈线对母线的谐波责任评估。绘制24 h内3条馈线的5次谐波责任变化趋势如图9所示。由图9结果显示，在5次谐波下，B线(化工厂专供线)对母线的电压畸变贡献最大，其次是A线(电气化铁路专供线)，而C线(飞机场专供线)对母线的电压畸变有抑制作用。

通过实际结果推断，在复杂的多谐波源电网系统中，谐波源有可能吸收谐波功率。同时，针对时变特性进行分析可知：A线和B线在凌晨对于5次谐波的贡献是全天中最大的，正午其次，即每日出现2次大的波动；而C线则在下午4时呈现爆发式的谐波贡献，其余时间的波动极小。同一谐波源在一天内不同时刻的谐波责任的波动，与负荷自身的功率及开启时间的长短等因素有关。评估结果与系统的实际运行状况吻合，可为谐波治理的责任划分提供参考依据。

6 结论

1)本文采用基于统计规律的相关性分析方法，从实测数据入手，分析了各谐波源电流与谐波超标节点谐波电压之间的映射关系，并通过偏相关分析识别出主要谐波源，进而评估多谐波源的谐波责任，解决了因多谐波源相关性导致谐波责任评估不准确的问题。

2)本文考虑了谐波的时变和波动特性，采用负荷水平等级分割法将谐波数据归类储存，针对数据段逐个进行偏相关分析，通过偏相关系数识别主要谐波源；通过改变时间尺度，研究了各谐波源的谐波责任时变特性。

3)本文考虑分布式电源接入电网场景下的谐波问题，提出了基于数据统计相关性分析的多谐波源责任评估方法，可评估含多个分散式谐波源系统的时变谐波责任。

仿真分析和实测结果表明，本文所提方法在含多谐波源的电网系统中能有效地评估各谐波源的谐波责任，且对主要谐波源有较好的评估精度。方法直观简易、工程实践性强，可为建立公平有效的谐波奖惩机制和标准来限制和治理谐波污染提供参考。