杨 晨，余一平，樊 陈，金 标，舒石泷

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

0 引言

随着电网规模的不断扩大、交直流系统互联以及高比例新能源与电力电子设备接入电网，低频振荡问题日益突出，制约了电能的远距离外送，严重影响着电力系统的安全稳定运行［1-2］。对低频振荡模式进行准确分析是抑制低频振荡的基础，这对保障电网安全稳定运行具有重要意义［3］。近年来，随着广域量测技术与宽频量测技术的发展，其为振荡模式在线辨识提供了技术支撑［4-5］。基于先进的量测技术能够实时获取系统动态信号，从而实现电网振荡模式的在线辨识。

实践证明，对于短路故障等扰动引发的电网振荡事件，采用Prony方法可以较为准确地从中获得系统模式信息。但是由于该类方法只能在明显振荡发生后产生告警且扰动的发生概率相对较小，因此其难以应用于日常监测电网的动态模式特征。电网日常运行中存在大量的小扰动，这些小扰动亦会激发电网的小幅动态变化，其动态轨迹中含有丰富的电网动态模式信息，若能基于日常小扰动响应识别电网存在的动态模式特征，这对运行人员进行振荡预警具有重要价值。具体而言，日常小扰动响应实测信号可分为2 种类型：一种是已知激励的系统响应信号，一般是由牵引负荷或是电容器投切等已知时间规律扰动引起的一定幅度的振荡；另一种是环境激励引发的系统响应信号，由系统运行过程中持续存在的负荷投切等随机性质小扰动引起。

常用的基于动态轨迹信号辨识模式信号方法有Prony 方法、希尔伯特-黄变换（Hilbert-Huang transform，HHT）方法等。Prony 方法可以直接提取信号中的振荡模态的幅值、初始相角、衰减因子和频率，计算简便，但该方法对噪声敏感且易产生虚假模态，为此文献［6-7］分别采用差分正交匹配追踪、形态滤波对Prony 方法进行改进；HHT 方法可以对非线性、非平稳信号进行分析，具备较好的时频分辨率，但该方法存在端点效应和模态混叠问题，文献［8］针对这些问题进行了改进，但仍存在虚假分量。从环境激励响应信号中提取低频振荡模式参数的方法有很多，文献［9-10］分别利用随机减量技术（random decrement technique，RDT）和变分模态分解（variational mode decomposition，VMD）从环境激励响应信号中提取自由振荡信号再进行振荡模式信息的辨识。该类方法的辨识效果取决于提取算法的准确度，且计算相对复杂。直接对环境激励响应信号进行分析的方法有随机子空间法、自回归滑动平均（auto-regressive moving average，ARMA）方法等。文献［11］提出了改进的随机子空间递推辨识方法，提高了抗噪性能，但该方法采用了奇异值分解使得计算速度受到限制。文献［12］采用滑动窗技术实现基于ARMA 方法的低频振荡模式参数追踪，但是滑动窗大小的选取对结果有影响，且更新结果需重新计算，计算量较大。文献［13］采用递推最小二乘法对ARMA 参数进行了追踪，但其白噪声估计过程计算较慢。文献［14］使用自回归（auto-regressive，AR）模型估计白噪声，提出用两段递推增广最小二乘方法求解ARMA参数，提高了参数估计速度，但是迭代收敛速度较慢。

本文根据高阶AR 模型与ARMA 模型等效原理［15］，采用递推最小二乘法拟合高阶AR 模型，并将它们等效转换为ARMA 模型，得到一种基于ARMA模型的两段最小二乘法，该方法不需要对噪声进行估计，迭代收敛速度较常规递推ARMA 辨识方法更快，辨识准确度更高，适用于日常小扰动响应中的低频振荡模式在线追踪。通过对给定ARMA 模型信号和新英格兰10 机39 节点系统仿真所得环境激励响应信号进行辨识，验证了本文所提方法的正确性。

基于已知激励和环境激励的响应信号都可以实现电网低频振荡模式在线监测，但现有工作主要集中于这2 类信号各自辨识方法的原理与改进，因此本文在已知激励响应与环境激励响应基本原理有区别的基础上，提出了低频振荡在线模式信息辨识方案，系统地分析了已知激励响应信号与环境激励响应信号的不同，得出在确定激励位置、观测点选择和响应模式间对应关系时，已知激励响应的辨识效果更好的结论。实际工程中，经过系统模式分析和激励源调研选择合适的观测点后可以选取已知激励响应信号作为主要分析对象。

1 ARMA模型的两段最小二乘估计原理

1.1 ARMA模型

对于观测时序｛x(t)｝建立的ARMA（n，m）模型，当模型阶次给定时，模型形式为：

式中：φi为自回归参数；θi为滑动平均参数；n、m分别为自回归部分的阶数和滑动平均部分的阶数；a(t)为均值是0的高斯白噪声。

ARMA模型可视为一a(t)为输入、x(t)为输出的线性时不变系统，其传递函数为H(z)=D(z)/A(z)，其极点即为A(z)=0 的根λi和λ*i（i=1，2，…，k，其中k为极点对数），根据式（3）即可得到振荡模式频率和阻尼比。

式中：Δt为采样时间间隔；fi为振荡模式频率；ξi为阻尼比。

1.2 两段最小二乘法

根据式（2）所示模型的平稳性与可逆性假设［15］，可证明式（2）等价于无穷阶AR(∞)模型，即：

式中：β0=1，βj→0(j→∞)。若取n0充分大，则有近似的高阶AR(n0)模型为：

使用递推最小二乘法对上式进行求解可获得β的估计值为：

由式（2）和式（4）可得到近似关系如下：

比较式（7）两边z-1的系数，则存在如下的关系：

式中：di=0(i＞m)；βi=0(i＞n0)。

从而得到d和D的估计值，最后根据式（8）可计算得到a的估计值。

为验证本文所提方法的正确性，将该方法与文献［14］中的两段递推增广最小二乘法进行比较，给定如下ARMA模型：

式中：a1=0.8；a2=0.4；b1=0.4；a(t)具体表示均值为0、方差为1的高斯白噪声。

根据式（11）递推生成x(t)的5 000 个数据，使用2 种方法对3 个参数进行辨识，结果如图1 所示。由图可知，本文给出的两段最小二乘法结果收敛更快，辨识精度更高。该方法把ARMA模型近似为一个高阶AR 模型，根据模型等价原理采用最小二乘法直接根据时间序列｛x(t)｝辨识ARMA 模型参数，避免了文献［14］中两段递推增广最小二乘法中对噪声项的估计，参数辨识效果较好，因此本文选择使用该方法对环境激励响应信号进行辨识。

图1 2种方法估计结果对比Fig.1 Comparison of estimation results between two methods

2 已知激励与环境激励响应机理

2.1 已知激励响应机理

以附录A 图A1 所示的单机无穷大系统为例，假定发电机模型采用经典二阶模型（发电机暂态电动势E′恒定），发电机内电势为Eg∠δ；节点2 为负荷节点，节点电压为Vd∠θd；节点3 为无穷大节点，节点电压为Vs∠0°。若存在负荷扰动，则节点注入功率方程线性化后可以表示为［16］：

式中：ΔPe、ΔQe分别为发电机有功功率和无功功率的偏差；ΔPd、ΔQd分别为负荷有功功率和无功功率的偏差；Δδ、ΔEg分别为发电机内电势相角和幅值的偏差；Δθd、ΔVd分别为节点电压相角和幅值的偏差；H、N、M、L分别为Pe、Qe、Pd、Qd对δ、Eg、θd、Vd的偏导数，下标为g 对应发电机节点状态量，下标为d 对应负荷节点状态量。

在经典模型情况下，可以求得：

式中：KS为发电机同步转矩系数；KP和KQ分别为与有功负荷扰动和无功负荷扰动相关的系数，均为实数。

由文献［17］可知冲击性负荷可表示为阶跃负荷。根据式（13）可知，有功负荷和无功负荷阶跃变化都会使机组电磁功率存在冲击性的分量。为方便分析，假设仅存在有功阶跃扰动，即ΔPd=Fε(t)、ΔQd=0，其中ε(t)为单位阶跃函数，F为阶跃幅值，将其代入发电机线性化转子运动方程并忽略机械功率后可得：

式中：ω0=2πf0，f0为系统基准频率；TJ为发电机惯性时间常数；KD为发电机阻尼系数。

式（15）等号右侧第1 项为由初始状态引起的自由振荡分量；第2 项为由阶跃扰动引起的自由振荡分量，其幅值与初始条件无关，与负荷扰动的幅值有关；第3 项为由阶跃扰动造成的稳态分量。可见在冲击性负荷等已知激励下发电机的时域响应形式为自由振荡。

2.2 环境激励响应机理

电网在日常运行过程中，负荷随机波动等环境激励会给系统带来小扰动。在环境激励下，发电机的转子运动方程可表示为：

式中：W(t)为随机激励；σ0为随机激励强度。

由式（16）可以看出发电机为单自由度系统，根据随机振荡理论［18］，系统的频率响应函数H(ω)可表示为：

式中：C为常数。

若已知系统的频率响应函数，其响应的功率谱Su(ω)即可由激励的功率谱SW(ω)求得，即：

可见在环境激励响应信号的功率谱中存在振荡频率为ωn的振荡模式。因此，可以通过电力系统在环境激励下的响应信号中提取振荡模式参数。

3 2种类型信号的比较

基于上述对已知激励响应和环境激励响应机理的介绍，初步比较这2种信号的不同之处有如下3个方面。

1）激励大小与响应形式不同。已知激励一般为冲击性负荷扰动或电容器投切，激励的具体形式已知，扰动幅值相对较大，响应信号也较强，响应形式一般为负阻尼振荡；而环境激励扰动幅值小，响应信号为类噪声信号，其幅值较小，从中难以观察到振荡现象。

2）信号成分不同。已知激励响应信号主要包含振荡模式信息，而环境激励响应信号除了包含振荡模式信息外混杂着噪声。

3）数据量不同。牵引负荷和变电站无功投切等冲击性负荷扰动时间相对固定，每天各个时段都会发生，且地理分布广，数据量较多；环境激励的信号时刻存在，数据量庞大。因此冲击性负荷这一类已知激励的响应信号和环境激励的类噪声信号都可用于振荡模式信息的在线监测。

根据上述比较可知，当量测信号为已知激励响应信号时可使用Prony方法进行辨识，量测信号为环境激励响应信号时可采用第1 节提出的基于两段最小二乘估计的ARMA方法进行辨识。低频振荡在线模式信息辨识方案流程见附录A图A2。

4 算例分析

为了对已知激励响应与环境激励响应进行比较，本文使用基于MATLAB 的PSAT（Power System Analysis Toolbox）搭建新英格兰10机39节点系统进行仿真测试，其结构见附录A图A3。

通过对10 机39 节点系统中冲击性负荷扰动获得的已知激励响应信号与负荷随机扰动获得的环境激励响应信号进行功率谱分析，并分别使用Prony方法和基于两段最小二乘估计的ARMA方法进行模式辨识，将辨识结果与特征值分析结果进行比较，从而得到用于辨识电网模式特征的2种信号的优缺点。

将系统在平衡点处线性化进行小干扰稳定分析，由分析结果可知10 机39 节点系统中共存在9 个低频振荡模式，其中阻尼比相对较弱的7 个振荡模式如附录A表A1所示。

4.1 冲击性负荷扰动

根据文献［17］中冲击性负荷的特性，在仿真系统的节点设置负荷阶跃变化以模拟实际电网中的冲击性负荷给电网带来的扰动。

实际电网中相量量测单元（phasor measurement unit，PMU）和宽频量测可以测量多种机械量和电气量，如发电机功角、转速、电压、电磁功率和无功功率等，上述几种量测量都包含振荡模式信息，但实际监测中几种量测信号的辨识结果会存在一定差异。为选择出合适的分析对象，首先在10机39节点系统中节点25 处分别设置30 MW 的有功冲击性负荷扰动和30 Mvar 的无功冲击性负荷扰动，对7 号机组功角δ7、转速ω7、机端电压V7、机组电磁功率Pe7和无功功率Q7的响应轨迹进行功率谱分析，除功角外各个量测量都转换为标幺值进行计算，得到的功率谱如图2所示。

图2 冲击性负荷扰动下各量测信号的功率谱Fig.2 Power spectrum of measurement signals under impulsive load disturbance

由图2 可知：有功冲击性负荷扰动下各个量测量的频率成分相同，其中电磁功率轨迹的功率谱幅值最大；无功冲击性负荷扰动下各个量测量的频率成分存在差别，电磁功率的功率谱包含的振荡模式多一些且幅值最大。因此在冲击性负荷扰动下低频振荡在线辨识的分析对象可以选用机组电磁功率。

然后分析冲击性负荷扰动响应信号作为低频振荡在线监测数据的适用性，先对不同位置的已知激励扰动响应信号进行成分分析，分别设置节点12、节点23、节点25处负荷在10 s时增加30 MW 有功冲击性负荷扰动，同样地在上述几个节点处的负荷设置增加30 Mvar 无功冲击性负荷扰动，将3 号机组、7号机组和8 号机组的电磁功率作为分析对象。图3为冲击性负荷扰动下7 号机组的电磁功率（标幺值）的响应轨迹，其余2 台机组电磁功率响应轨迹见附录A图A4。

图3 冲击性负荷扰动下电磁功率响应曲线Fig.3 Electromagnetic power response curves under impulsive load disturbance

由图3 和图A4 可见，冲击性有功负荷扰动和无功负荷扰动的机组响应都使机组发生了负阻尼振荡，其中与扰动源电气距离近的机组电磁功率响应振幅较大，且有功负荷扰动的振幅要大于无功负荷扰动的振幅。虽然扰动点不同，但电磁功率的时域响应波形都具有相同的特征，电磁功率在振荡结束后稳定在一个与扰动前不同的值，这是由于负荷功率阶跃变化引起的增量，与式（15）中等号右侧的第3项相对应，其中无功负荷扰动引起的增量较小。

对图3 中7 号机组在不同位置扰动下的电磁功率响应轨迹进行功率谱分析，结果如图4 所示。由图可知，冲击性负荷扰动响应都包含了一个区域间振荡模式（将该模式记作模式1）与多个本地机组间振荡模式（如模式2 等），但不同节点处的有功和无功扰动引起本地机组间振荡模式不同，这说明扰动源位置与扰动类型的不同会激发不同的本地机组间振荡模式。

图4 不同位置冲击性负荷扰动下电磁功率功率谱Fig.4 Power spectrum of electromagnetic power under impulsive load disturbance at different positions

对图3 中节点23 处冲击性有功负荷扰动下7 号机组电磁功率响应信号进行Prony分析，辨识结果如表1 所示。由表可见，频率和阻尼比与附录A 表A1中小干扰稳定分析结果之间的误差都较小，辨识结果较为准确，能够满足在线辨识的需求，因此冲击性负荷扰动响应信号可以用于振荡模式分析。

表1 冲击性负荷扰动响应信号分析结果Table 1 Analysis result of impulsive load disturbance response signals

然后分析在同一位置冲击性负荷扰动下响应信号的成分，对节点12 处负荷设置的有功功率和无功功率阶跃扰动产生的响应进行分析，将距离扰动源由小到大的3 号机组、6 号机组与8 号机组的电磁功率作为分析对象。3 台机组电磁功率响应轨迹的功率谱见附录A 图A5。由图3和图A5可见，不同机组激发了不同的本地机组间振荡模式，这说明与扰动源不同电气距离的观测点会包含不同的本地机组间振荡模式。由于冲击性无功负荷扰动响应的幅值比冲击性有功负荷扰动响应的幅值小，其功率谱幅值也相对小一些。

在区域电网中，已知激励扰动响应信号都包含区域间振荡模式成分，因此能用于监测区域间振荡模式，但由于量测信号观测点的不同，监测到的本地机组间振荡模式不同，所以本地机组间振荡模式需要选择合适的观测点来进行监测。

根据冲击性负荷扰动响应信号的分析结果可知，冲击性负荷响应信号包含多个振荡模式，可以反映系统中真实的振荡频率、阻尼信息，在实际工程中如牵引列车和变电站电容器投切的冲击性负荷扰动时间固定，响应数据较多，可以达到与振荡事件信号分析同样的效果，能够有效实现在固定时段对系统振荡模式特征进行在线监测，但由于冲击性负荷扰动程度比短路故障扰动小得多，实际电网中并非所有节点都能检测到已知激励响应信号，因此需要根据日常电网激励特征选择激励点和根据模式可观测性和节点响应强度选择观测点。

4.2 随机扰动

由文献［19］可知实际电网负荷随机扰动呈高斯分布，为模拟实际电网中小幅度随机扰动得到环境激励响应信号，分别向负荷节点注入随机扰动信号，该信号为高斯白噪声信号。

与4.1节类似，为选择出随机扰动下合适的分析对象，在10 机39 节点系统中向节点12 处的负荷注入标准差为0.1 MW 的高斯白噪声，采集7 号机组功角、转速、机端电压、机组电磁功率和无功功率，将除功角外的各个量测量转换为标幺值后计算得到如图5 所示的功率谱。可见在随机扰动下电磁功率的功率谱幅值较大，即电磁功率响应信号较其他量测信号明显，且包含的振荡模式也较多，因此在随机扰动下将电磁功率作为低频振荡信息辨识分析对象。

图5 随机负荷扰动下各量测信号的功率谱Fig.5 Power spectrum of measurement signals under random load disturbance

然后比较环境激励响应信号与已知激励响应信号包含的信息，分析3 号机组的电磁功率Pe3与6 号机组的电磁功率Pe6响应信号的功率谱，并与已知激励响应信号功率谱进行对比，如图6所示。

图6 环境激励响应信号与已知激励响应信号功率谱对比Fig.6 Comparison of power spectrum between environmental excitation response signals and known excitation response signals

由图6 可知，2 种类型信号的频率成分相同，已知激励响应信号功率谱与环境激励响应信号功率谱的包络线形状类似，但由于环境激励响应幅值较小，其功率谱幅值比已知激励响应信号功率谱的幅值要小一些。环境激励响应信号的功率谱中振荡模式对应频率的幅值较大，因此根据环境激励响应信号能获取系统振荡模式信息，但环境激励响应信号功率谱中振荡频率附近的频率对应幅值也较大，使得环境激励下6 号机组电磁功率响应信号中包含的2 个本地机组间振荡模式区分得不太明显。环境激励响应信号功率谱中振荡频率范围外的其他频率成分幅值相当，类似于高斯白噪声，会对振荡模式辨识造成干扰。

为分析环境激励响应信号作为低频振荡在线监测数据的适用性，先分析随机扰动源位置对响应信号成分的影响，分别向节点23和节点25处的负荷注入标准差为0.1 MW 的高斯白噪声，分析3 号机组电磁功率，3 号机组在节点23 处的负荷随机扰动下的响应情况见附录A图A6。

相较4.1节中的已知激励响应信号而言，环境激励响应信号幅值比已知激励响应小很多，且时域波形较为杂乱。对节点23 和节点25 处随机扰动的电磁功率响应信号进行功率谱分析，其功率谱见附录A 图A7。可见不同节点的负荷扰动激发的电磁功率响应信号的频率成分存在不同，节点23和节点25处的随机扰动都激发了振荡模式1，但激发的本地机组间振荡模式不同，因此与已知激励响应信号类似，本地机组间振荡模式需要选择合适的观测点来进行监测。

然后向节点21、节点23、节点24 处的负荷注入高斯白噪声进行仿真，分析电网中观测点和随机负荷扰动源距离与响应信号间的关系，选取与随机扰动源的电气距离较近的6 号机组和电气距离较远的4 号机组的电磁功率作为分析对象，2 台机组电磁功率的功率谱如图7 所示，图中Pe4为4 号机组的电磁功率。

图7 电磁功率功率谱Fig.7 Power spectrum of electromagnetic power

由图7中2台机组的电磁功率响应信号可知，环境激励都激发出了区域间振荡模式1，由于6号机组与扰动源电气距离小，6 号机组电磁功率功率谱幅值比节点3电压功率谱的幅值大，因此6号机组电磁功率响应信号中噪声幅值也较大。6 号机组电磁功率响应信号中存在的明显噪声分量会给振荡模式辨识带来一定的误差。

使用基于两段最小二乘估计的ARMA 方法对2台机组的电磁功率响应信号中的模式1 进行辨识，算法迭代收敛后20 min数据的振荡模式辨识结果如图8 所示。对2 台机组电磁功率的辨识结果求均值得到如表2所示的最终的辨识结果。

图8 振荡模式的辨识结果对比Fig.8 Comparison of identification results for oscillation mode

表2 环境激励响应信号分析结果Table 2 Analysis result of environmental excitation response signal

由图8 和表2 可知，通过2 台机组的响应信号都能较为准确地辨识出区域间振荡模式的频率，但由于6 号机组距离随机扰动源近，信号中噪声的成分多，阻尼比辨识误差相对大一些。对比表1和表2可知，直接对环境激励响应信号进行分析的结果相较已知激励响应信号的辨识结果差一些。

由式（18）可知，仿真信号中的噪声成分是由随机激励产生的，实际工程应用中还会存在量测噪声和数据丢失的情况，使得辨识结果出现偏差。在6号机组电磁功率响应信号中增加高斯白噪声并去除少量数据点以模拟现实中出现的量测噪声和数据丢失的情况，将该信号记作受扰信号，分析受扰信号时丢失的数据采用线性插值的方式进行填补，在已知振荡频率后对受扰信号进行简单低通滤波处理，原信号与滤波前后受扰信号的振荡模式1 辨识结果对比如图9所示，对辨识结果求均值后得到如表3所示的最终辨识结果。

由图9和表3可知，当数据出现量测噪声和数据丢失时会使得振荡模式信息的辨识结果出现一定的偏差，经滤波处理后信号中噪声成分减少，阻尼比的辨识结果有所改善。本文给出的两段最小二乘法辨识结果较为准确，能满足日常监测要求，辨识结果与信号受干扰程度有关。因此对系统振荡模式在线辨识时，环境激励响应信号需要预先获取该网络的振荡模式并进行滤波等处理。

表3 环境激励响应信号分析结果Table 3 Analysis result of environmental excitation response signal

综上所述，相较于已知激励响应信号而言，环境激励响应信号包含的振荡模式信息相对较少，为监测到更多的振荡模式需要更多的量测信号；环境激励信号中包含一些噪声，噪声成分的多少与观测点和扰动源电气距离有关，若使用环境激励响应进行振荡模式辨识时需要对数据做预处理；电网中大多数节点都能够监测到环境激励响应，因此其观测点部署的难度较已知激励响应信号观测点部署低一些。

5 结论

基于PMU 或宽频量测信号进行振荡模式辨识是进行电力系统振荡模式在线监测的有效方式，本文给出了一种适用于日常小扰动响应振荡模式信息辨识的两段最小二乘的估计方法，能够准确有效地从系统响应数据中提取振荡模式信息。对日常小扰动响应的已知激励响应和随机环境激励响应进行了机理介绍，并采用功率谱与辨识方法进行综合分析比较，用于电网振荡模式特征在线监测的2 种响应信号优缺点总结如下。

1）2 种类型信号的激励形式不同，因此响应形式也存在差异，已知激励响应信号一般为负阻尼振荡，随机环境激励响应信号为类噪声信号。已知激励响应信号幅值较随机环境激励响应信号大，但其数据量相对较少。

2）已知激励响应信号和随机环境激励响应信号包含的频率成分类似，都能监测到振荡模式信息，但环境激励响应信号中振荡模式信号相对微弱且噪声成分较多，随机环境激励响应信号观测点部署相对简单。

3）随机环境激励响应信号易受噪声干扰且幅值较小，其辨识的难度与复杂程度相对较大，直接进行分析会存在一定误差，辨识效果与信号受干扰程度有关。

因此在进行系统模式分析和激励源调研分析，预先确定激励位置、观测点选择以及响应模式间的对应关系的情况下，已知激励响应的效果更好。在这种情况下可以将已知激励响应辨识作为日常系统振荡模式特征在线识别的主要手段，在无已知激励响应信号的情况下使用随机环境激励响应信号辨识作为补充。

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