基于宽频量测间谐波潮流计算的次／超同步振荡溯源方法

金 标，余一平，樊 陈，易文飞，鞠 平

（1. 河海大学能源与电气学院，江苏南京 211100；2. 中国电力科学研究院有限公司（南京），江苏南京 210003；3. 国网江苏省电力有限公司，江苏南京 211103）

0 引言

我国提出“力争2030 年前实现碳达峰，2060 年前实现碳中和”的战略目标。这要求我们构建以新能源为主体的新型电力系统，届时电力系统将呈现高比例可再生能源接入、高比例电力电子设备的“双高”特性［1］。而随着电力电子设备和可再生能源在电网中的比例逐渐升高，国内外频繁发生次／超同步振荡事故，其中典型的如2015 年以来我国新疆哈密直驱风机发生20～80 Hz 的振荡导致火电机组发生轴系扭振而跳闸［2］。新能源发电产生的次／超同步振荡传播可能引发邻近火电机组的轴系扭振，从而对电力系统的安全稳定运行造成更严重的威胁。文献［3⁃5］定义次同步频率下的功率，对次同步振荡的传播和分布进行分析，并利用次同步频率下功率的特性进行溯源，实现了次同步控制作用机理扰动源的辨识；文献［6］通过将实时监测数据和聚合阻抗模型相结合实现基于次同步控制相互作用机理的次同步振荡扰动源的快速定位与切除；文献［7⁃9］将暂态能量流推广至次同步强迫振荡扰动源定位及阻尼评估，建立了暂态能量流与发电机组电功率以及阻尼之间的联系，并实现了次同步强迫振荡的振荡源定位。如何快速辨识新能源发电引发的次／超同步振荡传播路径和分析其对周围邻近火电机组的影响，这对于未来新型电力系统的安全稳定运行具有重要意义。文献［3⁃9］均为基于量测的次同步振荡定位方法，依赖量测的部署，而当前针对宽频振荡的量测配置较少，限制了次同步振荡定位方法的实现。目前针对次同步振荡的研究仍不完善，不同机理下的判据适用性还有待验证。

宽频量测技术的快速发展为研究新能源发电次／超同步振荡现象提供了间谐波振荡相量数据［10⁃11］。文献［12］通过同步相量测量单元PMU（Phasor Measurement Unit）的测量数据还原算法实现对间谐波振荡的相量测量。文献［13］基于自适应滤波算法提出了宽频电磁振荡多模态特征的相量检测方法，并实践开发了多模态振荡检测装置，解决了由多间谐波频率混叠导致的电磁振荡不准确问题。文献［14］针对电力电子化电网电气量的宽频特征，设计了多功能宽频量测装置。宽频量测可以获得新能源发电次／超同步振荡过程中间谐波电压电流的分布，有利于调度运行人员识别次／超同步振荡传播路径，并进行溯源定位［15］。但目前电网中针对次／超同步振荡相关的间谐波频率电压电流分量的宽频量测设备配置不足，只能在新能源发电并网节点和部分特定重要厂站、节点加以装设。由于测量装置改造和数据处理量大等因素，对电网中所有节点都装设间谐波频率分量的量测设备在短期内是不现实的，这导致了新能源发电次／超同步振荡发生时，电网中绝大多数节点的间谐波电压电流分量分布或振荡传播路径，电网运行调度人员是未知的。如果能够通过类似潮流计算的方式获取电网各节点间谐波频率下的电压电流分布，即通过将新能源并网点振荡宽频量测与间谐波潮流计算相结合获取以全网各节点的间谐波电压电流分布为特征的振荡传播信息，对于运行人员进行振荡分析控制具有重要实用价值。而新能源发电次／超同步振荡引发的间谐波潮流分布与谐波潮流分布具有相似性，以往的研究针对基波和整倍频的谐波潮流计算已经非常成熟。文献［16］提出了经典谐波潮流的详细建模和计算步骤，先给出非线性电力电子谐波源的计算公式，然后对基波与谐波单独进行求解并联立满足节点状态约束条件的计算公式。

本文在结合宽频量测和谐波潮流算法的基础上，提出了一种基于新能源并网点宽频量测间谐波潮流计算的次／超同步溯源方法。所提间谐波潮流计算方法利用新能源并网点的较少量测设备可以获取全网次同步振荡传播路径，弥补量测不足，并基于间谐波潮流给出新能源发电系统振荡引发邻近火电机组轴系扭振时的振荡甄别定位。首先给出了间谐波潮流计算的网络方程和电力元件模型，提出了间谐波潮流计算的详细步骤；其次当新能源场站发生次／超同步振荡时，基于部分新能源并网点宽频量测间谐波潮流方法计算间谐波电压电流在全网中的分布，进而获得次／超同步振荡在电网中的传播路径；然后调度运行人员通过新能源预想振荡情况下计算全网间谐波潮流，与部分节点宽频实测量的比较分析，可以实现新能源发电次／超同步振荡引发邻近火电机组轴系扭振溯源定位；最后通过四机两区域系统和IEEE 10 机39 节点系统算例验证了所提方法的准确性，为调度运行人员从电网角度进行次／超同步振荡分析溯源提供一种新的手段。

1 间谐波潮流分析方法

1.1 次／超同步振荡特征分析

新能源发电发生次／超同步振荡时，电网中各节点电压和各支路电流中存在次／超同步频率电压电流分量。此时不考虑迅速发散的暂态过程，对持续振荡阶段进行分析［17⁃18］，如图1 所示，图中电流为标幺值。当振荡波形呈现出周期性的平稳振荡时，对其进行傅里叶变换可以解耦出工频与间谐波频率分量，即解耦为基频分量和间谐波频率分量的正弦波。

图1 次／超同步振荡波形图Fig.1 Sub-／super-synchronous oscillation waveforms

因此，发生次／超同步振荡时节点电压和支路电流可以写成不同频率瞬时量的叠加形式，具体如下：

式中：u(t)、i(t)分别为节点电压、支路电流；u0、us1、us2分别为基频、频率为s1、频率为s2电压分量的幅值；i0、is1、is2分别为基频、频率为s1、频率为s2电流分量的幅值；ϕ0、ϕs1、ϕs2分别为基频、频率为s1、频率为s2电压分量的相位；φ0、φs1、φs2分别为基频、频率为s1、频率为s2电流分量的相位；ω0、ωs1、ωs2分别为基频、频率为s1、频率为s2的角频率。

综上所述，在次／超同步振荡平稳阶段时，节点电压、支路电流瞬时值可以表示为不同频率下的瞬时分量叠加，基频以外的瞬时周期性分量，称之为间谐波分量，而间谐波分量的电压电流相量形式及其功率分布，可以表述为间谐波潮流。

1.2 间谐波潮流网络方程

次／超同步振荡平稳阶段中，新能源发电发生次／超同步振荡时，电网中各节点电压和各支路电流中存在间谐波电压电流分量。电网中间谐波电压电流分布对掌握次／超同步振荡的传播路径至关重要。本文将通过间谐波潮流的方式将电压电流各状态量在频域上进行解耦计算，获得振荡频率分量下间谐波电压电流分布。

对于有n个节点的电力系统网络，各个频率分量始终满足基尔霍夫电流定律（KCL）条件，即满足如下电流平衡方程：

式中：In为n个节点流入或流出的某一频率的间谐波电流矩阵，Yn为对应频率下n个节点导纳矩阵，Vn为对应频率下n个节点间谐波电压矩阵。

对于一个大电网而言，已知新能源场站基于宽频量测获得的次／超同步振荡间谐波注入电流，所有接地节点的间谐波电压为0。因此，电流平衡方程可以写成如下形式：

式中：IB、II和IG分别为系统新能源并网边界节点、内部节点和接地节点的注入电流；UB、UI和UG分别为系统边界节点、内部节点和接地节点的电压；YBB、YII、YGG分别为系统边界节点、内部节点和接地节点的自导纳；YBI、YIB和YIG、YGI分别为边界节点与内部节点和内部节点与接地节点之间的互导纳。

为分析方便，本文研究单个新能源场站发生次同步振荡的情况，即电网中只存在单个振荡源、单个边界节点，不涉及次同步振荡引发其他电力系统元件振荡产生新的振荡源情况。若有多个新能源并网点，则仅需增加模型中的边界节点个数，在数学上建立的潮流网络方程没有特殊的改变。

网络中的节点依据边界条件的不同可以划分为图2 所示的接地节点、内部节点和边界节点。接地节点的间谐波频率分量电压为0；内部节点的间谐波频率分量注入电流为0；边界节点的次同步频率分量电流可以通过宽频量测辨识获取。

图2 系统节点分类图Fig.2 Diagram of system bus classification

因此，式（3）可以重新组合为如下形式：

不妨假设边界节点和内部节点共m个，编号为1—m。利用高斯-赛德尔迭代法求解，UB&I可由式（5）计算得到。

由式（5）可计算得到UI，又已知接地电压为0，即UG=0，可由式（6）计算得到IG。

由式（4）—（6）可知，在网络参数已知的系统中，仅需在新能源并网点设置宽频量测设备，发生新能源场站次／超同步振荡时便可计算得到间谐波电压电流在全网的传播与分布。

实际网络中，获取完全准确的网络间谐波模型是比较困难的，尤其是电网中复杂的动态元件、非线性的电力电子设备以及负荷等效模型。此时可以通过在未知结构参数的设备与模型出口侧加设宽频量测设备通过辨识来解决，例如在输电网与配电网联络线处或在电力电子出口侧设置宽频量测设备来获得负荷和电力电子设备间谐波频率下的等效阻抗。网络中宽频量测点的布置一般需要考虑新能源并网点、部分发电机组出口和负荷接入电网点。新能源并网点的等效建模困难且机理复杂需要设置量测点；发电机组和负荷可基于现有模型采用扫频法建立间谐波频率等效模型，当条件具备时可以安装宽频量测设备建立更精准的等效模型。与传统基频潮流相比，由于同步发电机等效为阻抗，间谐波潮流的量测节点只需要新能源并网点，远少于传统基频潮流。同理，与谐波潮流相比，间谐波潮流的量测节点少于统一解法和交替迭代解法的谐波潮流。

1.3 间谐波潮流元件模型

系统中线路、负荷、变压器等模型均与系统中流过的电流的频率相关，因此需要对原系统中基频分量下的模型进行修正［19］。重新构建谐波频率下的阻抗等值网络，得到新的节点导纳矩阵。

1）输电线路。

π型等值电路阻抗计算公式如下：

式中：zl为线路阻抗；yl为线路对地导纳；r、x、g和b分别为基频下线路电阻、电抗、电导和电纳；f和fn分别为间谐波频率和基频频率。

线路的电阻部分不受流过线路的电流频率变化影响，仅电抗部分与电流频率相关，需要根据频率进行修正。

2）负荷模型。

当仅考虑静态负荷模型时，恒阻抗负荷计算公式如下：

式中：R、Xc和Xl分别为负荷等效电阻、容抗和感抗；Vn为负荷额定电压；Pn、Qc和Ql分别为负荷有功功率、容性无功功率和感性无功功率。

对于异步电动机负荷模型，不考虑转子电阻，化简后的异步电动机负荷阻抗计算公式如下：

式中：ZM为异步电动机负荷阻抗；Rs为定子电阻；Ls、Lr、Lsr分别为定子电感、转子电感和互感；fs为异步电动机转差频率。

3）变压器。

双绕组变压器阻抗计算公式如下：

式中：kT为变压器高低侧变比的比值；下标j、i分别为变压器高、低压侧节点编号；Uj、Ui分别为变压器高、低压侧电压；Ij、Ii分别为变压器高、低压侧电流；Zb为变压器等效阻抗。

需要注意的是，式（10）中的Zb同样需要进行基于频率的阻抗变换。具体过程与输电线路的变换保持一致。对于三绕组变压器，简化为3 个双绕组变压器进行处理。

4）同步发电机。

对于恒励磁的同步发电机在非共振间谐波频率下的模型，忽略阻抗简化后的同步发电机阻抗计算公式如下：

式中：ZG为同步发电机的等效阻抗；X″为同步发电机次暂态电抗。

对于同步发电机建模，需要了解间谐波电流在同步发电机中引发的电磁变化。系统发生次／超同步振荡，发电机端口流入间谐波频率为f的电流，该电流在发电机中产生时变的磁场，即气隙中存在频率为fn-f的dq轴磁通量，在绕组上产生频率为fn-f的阻尼电流。由于流入电机的间谐波电流恒定，导致阻尼电流也持续振荡，所以应该选择建立计及阻尼绕组的次暂态方程。非共振频率下转速约为1 p.u.，恒励磁条件下励磁电压的间谐波频率分量为0。综上，联立求解可得到同步发电机的阻抗计算公式如式（11）所示。

1.4 间谐波潮流计算步骤

间谐波潮流计算框图如附录A 图A1 所示，详细步骤如下。

1）基于新能源并网点宽频量测设备，获取并网点振荡电流、电压瞬时数据，通过辨识方法对量测数据进行解耦，本文采用Prony方法得到并网点注入的间谐波频率电流分量；若存在多个间谐波频率分量，则选取主导次同步振荡频率分量进行间谐波潮流计算，得到该频率下间谐波潮流分布。

2）对基频网络元件进行校正，形成间谐波频率f下的节点导纳矩阵Yn。

3）设置初始值（将接地电压设置为0，设置边界节点电流或电压）。

4）设置收敛指标，利用式（3）—（5）对间谐波潮流迭代计算。

5）打印、输出结果。

2 次同步振荡溯源分析

2.1 振荡传播路径获取

新能源并网发生次／超同步振荡时，利用并网点宽频量测信息获得并网点间谐波电流，基于间谐波潮流计算方法，获取全网所有节点次／超同步频率分量电流电压信息。间谐波电流流过输电线路，间谐波电压将产生一个线路压降，由间谐波分量电流流向判定线路上间谐波分量传播的方向；流过网络节点，间谐波分量仍然满足电流守恒，可以由KCL判定电流的分流与汇聚；由线路间谐波电流大小及传播方向和节点上间谐波电流流入流出方向获取次／超同步振荡传播路径。网络线路与节点图见附录A图A2。传播路径计算如式（12）和式（13）所示。

式中：Ui和Uj分别为同一线路上节点i和节点j的电压，Z为该线路上节点i和节点j间的等效阻抗，I为流过该线路的电流；Ii为流入某一节点i的不同支路电流。

2.2 邻近火电机组被激发轴系扭振溯源

新能源发电通过风火打捆向外输电时，当风电发生次同步振荡且振荡频率与邻近火电机组扭振自然频率互补，形成共振，将激发火电机组轴系扭振。此时火电机组的间谐波频率等效阻抗将发生突变，与非共振频率下火电机组的等效阻抗有较大差异。调度运行人员选择火电机组非扭振情况下的等效阻抗进行预想新能源并网次同步振荡下间谐波潮流分析，若预想振荡间谐波潮流计算结果与宽频量测实测结果存在显著差别，可能是由于火电机组发生轴系扭振造成间谐波电流真实分布发生改变。由于存在共振机组等效阻抗的差别，因此体现在预想振荡间谐波潮流中为发生扭振机组节点的间谐波电压电流与宽频量测实际值的误差较大，远大于其他节点，且不在一个数量级上。因此可以通过计算各个火电机组预想计算值与实测量测值的相对误差，存在显著差异的机组即为轴系扭振机组。这一现象可作为新能源并网次同步振荡引发邻近火电机组轴系扭振溯源定位的依据。

3 间谐波潮流振荡传播路径分析

3.1 四机两区域系统

为验证间谐波潮流计算的准确性，通过时域仿真结果与间谐波潮流计算结果进行对比分析。基于MATLAB／Simulink软件搭建时域仿真模型，本节选用四机两区域系统进行验证。系统模型由4台同步发电机、对称的2个区域构成，2条并联220 km输电线路连接2 个区域，并在节点N9处加设模拟新能源发电次／超同步振荡源，向系统注入次同步频率分量电流，其间谐波系统阻抗网络见附录A 图A3，具体网络参数见附录B表B1。

新能源并网振荡设置三相20 Hz 正弦波恒定电流，幅值为5 kA。通过宽频量测得到新能源并网点向系统注入间谐波电流5 kA，利用间谐波潮流方法进行计算，求解得到间谐波电压电流的分布情况。

由宽频量测设备获取新能源并网点出口电压电流信息，其频域如图3（a）所示，主要监测到60 Hz的工频分量与20 Hz 的间谐波频率分量。选择1.0～1.2 s的振荡窗口，由Prony辨识工频分量和间谐波频率分量在1 s时间断面上的幅值和相位，并由辨识结果计算得到拟合波形，结果如图3（b）所示，拟合波形与量测波形保持一致，可见次同步振荡电压电流瞬时值可以解耦为基频分量和间谐波分量。通过间谐波潮流计算，对其他电网节点电压、支路电流进行评估，其中支路N5-N6的电流波形如图3（c）所示，可见量测波形与间谐波潮流计算后重新叠加基频波形保持一致。综上所述，在新能源并网发生次同步振荡的情况下，可以由相量法的间谐波潮流计算来获取间谐波电压电流在网络中的分布。在具体振荡发生时，振荡频率由宽频量测辨识获取，因此间谐波潮流表征振荡仅需确定该频率分量的振荡幅值、振荡相位和相位对应的时刻即可推导出在时域上连续的振荡波形。为简便表达，下文只展示振荡差异的幅值和相位。

图3 间谐波潮流计算结果Fig.3 Calculative results of interharmonic power flow

具体地，间谐波潮流计算结果与时域仿真结果对比见表1和表2。

表1 四机两区域系统间谐波电压分布结果Table 1 Results of interharmonic voltage distribution of four-machine two-area system

表2 四机两区域系统间谐波电流分布结果Table 2 Results of interharmonic current distribution of four-machine two-area system

在四机两区域系统中，由时域仿真结果与间谐波潮流计算结果对比可见，间谐波潮流计算结果更准确。节点间谐波电压电流计算结果无论是在幅值还是在相位上都有较好的精确度。

通过间谐波潮流计算结果可以分析次同步振荡在系统中传播与分布路径。由计算结果可知，随着与次同步振荡源的电气距离增大，次同步频率分量电压逐渐减小。在单个振荡源情况下，由于网络呈现树干式结构，电压单调下降。由间谐波潮流计算电流结果可知次同步振荡在网络各个方向传播分布的情况。与间谐波振荡源电气距离最近的GG4分流得到的振荡电流最大，区域2 节点振荡电流同比要高于区域1对称节点振荡电流，区域2受次同步振荡影响更大。

3.2 IEEE 10机39节点系统

本节再选择IEEE 10机39节点系统验证间谐波潮流计算的准确性，IEEE 10 机39 节点系统间谐波阻抗网络图见附录B 图B1，具体网络参数见附录B表B2。该系统包括10 台发电机、39 个节点、12 台变压器、34 条线路，代表新英格兰的一个345 kV 电力网络，其中机组G1为系统与外网连接的等值机组。本文选择在电力网络中间位置母线16 处，加设模拟新能源并网次／超同步振荡源，向系统注入次同步频率分量电流。

新能源并网点设置三相20 Hz 正弦波电流，幅值为500 A。IEEE 10 机39 节点系统较大，结构复杂，支干路较多，对电流分流作用较大，单条线路上基频电流较小，在IEEE 10 机39 节点系统内观察间谐波在大电网传播分布的详细情况。

IEEE 10 机39 节点系统所有发电机参数未知。计及新能源并网点节点16 注入电流500 A，系统共存在11 个边界节点，需要在发电机出口、负荷与电网联络线及新能源并网点设置宽频量测设备。同四机两区域系统一致，首先通过新能源并网点辨识间谐波振荡频率，然后由Prony辨识出边界节点间谐波分量信息，最后通过间谐波潮流计算求解间谐波电压电流分布，其电压时域仿真结果与间谐波潮流计算结果如图4所示。

图4 IEEE 10机39节点系统间谐波潮流电压分布Fig.4 Voltage distribution of interharmonic power flow of IEEE 10-machine 39-bus system

由图4可见，IEEE 10机39节点系统的间谐波潮流计算结果与时域仿真结果两者相近，大致准确。仿真得到的支路电流结果见附录B表B2。

节点12 和节点31 间谐波潮流计算结果与时域仿真结果存在较大误差。这是由于时域仿真得到的量测数据在Prony辨识中解耦效果不好，导致仿真数据失真与间谐波潮流计算结果相差较大。总体而言间谐波潮流计算结果与时域仿真结果一致，验证了间谐波潮流计算的准确性。

在单振荡源情形下，通过间谐波潮流获取网络电压信息和电流信息，监测得到在IEEE 10 机39 节点系统中间谐波潮流传播图见附录B图B2［10］。

4 间谐波潮流火电扭振溯源

4.1 风火打捆系统扭振

本节构建三质量块的火电机组模型，由多质量块扭振自然频率计算公式得到自然扭振频率f1=32.407 Hz和f2=24.719 Hz。

通过MATLAB／Simulink 在单机无穷大系统中搭建上述火电机组模型，通过频率扫描法监测机组端口的电压与电流数据。将量测数据通过Prony 辨识、解耦获取其中的间谐波电压与间谐波电流成分，然后通过计算获得火电机组的等效阻抗。采用频率扫描法得到的火电机组等效阻抗如图5所示。

图5 采用频率扫描法得到的火电机组等效阻抗Fig.5 Equivalent impedance of thermal power unit with frequency scanning method

由图5 可知，随着频率升高，发电机的等效阻抗幅值逐渐增大，相位保持在90°，与式（11）所示同步发电机的等效阻抗模型一致。在火电机组的扭振自然频率互补f0-f1、f0-f2处，由于电机发生轴系扭振，转子角速度偏移较大，不满足非共振频率下转速约为1 p.u.，式（11）不再适用，其中f0为基频。此时，共振频率互补导致发电机的等效阻抗幅值和相位都发生突变，即火电机组轴系扭振导致机组等效阻抗突变，使间谐波潮流结果存在较大偏差。在本文的2个仿真案例中，输电线路、变压器、恒功率负荷、电动机负荷等都为无源元件，其元件电抗对频率响应为线性的，不产生突变。因此，网络中间谐波阻抗突变是由火电机组轴系扭振所导致的。

4.2 火电机组扭振溯源定位

调度运行人员日常可以通过间谐波潮流计算获得预想新能源并网次同步振荡情况下全网的振荡传播分布。本节继续选用IEEE 10机39节点系统进行仿真分析，运行人员通过间谐波潮流计算已经掌握电网中预想未发生火电机组轴系扭振情况下的节点间谐波电压分布。而当发电机G8由于新能源场站次同步振荡传播激发轴系扭振时，其间谐波等效阻抗发生突变，附近节点的间谐波电压电流分布将发生变化，与之前预想未发生扭振情况下潮流结果进行对比，可求得各个节点间谐波电压的相对偏差。机组端口电压计算值与实际值对比如表3所示。

表3 机组端口电压计算值与实际值对比Table 3 Comparison of calculative and actual values of generator port voltage

火电机组发生扭振会引起间谐波等效阻抗变化，从而影响间谐波潮流分布。由表3 可得，机组G8的相对误差最大，偏差高达-37.551%，远大于其他节点，具有良好的可辨识性。而机组G8正是设定发生扭振的机组，辨识正确。实际应用中当电网中部分节点配置宽频量测装置时，如果实测得到的间谐波电压电流分布与预想未发生扭振情况下的间谐波电压电流分布存在较大偏差，则可能引发了火电机组扭振。

5 结论

本文针对电网中出现的次／超同步振荡现象，提出了扩展至工频及其整倍频外的间谐波潮流计算方法，给出了完整可知网络系统下间谐波潮流模型。对电力系统中的电力元件模型进行了校正，给出了间谐波频率下电力系统元件的等效阻抗模型。详细介绍了间谐波潮流计算方法的步骤、适用条件。在四机两区域系统以及IEEE 10机39节点系统中与时域仿真方法对比验证所提间谐波潮流计算方法的有效性和准确性。最后给出了间谐波潮流分布在新能源次同步振荡引发邻近火电机组轴系扭振识别中的应用。

本文所提基于新能源并网点宽频量测的间谐波潮流计算方法，逻辑清晰、计算简单便捷。利用所提方法能够有效地分析新能源次／超同步振荡时间谐波电压电流在电网中的传播和分布，进而获得振荡传播路径，并识别出风电次同步振荡激发邻近火电机组轴系扭振位置。

附录见本刊网络版（http：//www.epae.cn）。