CN110412349B

CN110412349B - 基于插值dft的同步相量数据次同步振荡参数辨识法

Info

Publication number: CN110412349B
Application number: CN201910795655.5A
Authority: CN
Inventors: 杨晓梅; 王杨; 肖先勇; 张家宁; 郭林明
Original assignee: Sichuan University
Current assignee: Sichuan University
Priority date: 2019-08-27
Filing date: 2019-08-27
Publication date: 2021-07-23
Anticipated expiration: 2039-08-27
Also published as: CN110412349A

Abstract

基于插值DFT的同步相量数据次同步振荡参数辨识法，包括以下步骤：步骤一、获取实时高分辨率波形数据；步骤二、在主站获得同步相量X_r(m)；步骤三、在矩形窗内对同步相量X_r(m)进行离散傅里叶变换操作，得离散傅里叶变换谱F_r(k)；步骤四、确定次同步振荡幅度最大的频谱线k_m；步骤五、构造指标相量[k1，k2，k3，k4]，引入二阶差分光谱的频谱F_rs(k_i)，获得简化后的二阶差分比值R；步骤六、计算次同步振荡分量的频率f_s和次同步振荡分量的阻尼系数α_s的值；步骤七、对同步相量X_r(m)进行Hann窗口截取并执行离散傅里叶变换操作，反推得次同步振荡分量幅值A_s计算公式；通过本方法，能准确辨识次同步振荡分量的参数，利于缓解次同步振荡对系统稳定性和设备安全性的影响。

Description

基于插值DFT的同步相量数据次同步振荡参数辨识法

技术领域

本发明涉及辨识电网中次同步振荡参数从而制定电气控制策略的领域，具体是一种基于插值DFT的同步相量数据次同步振荡参数辨识法。

背景技术

次同步振荡(SSO)是电力系统机电侧相互作用引起的一种异常现象。近年来，随着可再生能源的快速增长和应用，全球范围发生了多次有关于次同步振荡的电力故障事件。在这些事故中，由于次同步振荡较为严重，对汽轮发电机转子轴系造成了严重的破坏，甚至削弱了电力系统的安全可靠性。因此，对次同步振荡进行监测和识别，对制定电气控制策略和实施次同步振荡监控具有重要意义。

发明内容

本发明目的在于解决由于电力系统机电侧相互作用引起的次同步振荡的参数监测不准确的问题，提供了一种基于插值DFT的同步相量数据次同步振荡参数辨识法，通过使用本方法，可以准确辨识次同步振荡的参数，有效制定电气控制策略，缓解次同步振荡对电力系统稳定性和设备安全性造成的影响。

本发明主要通过以下技术方案实现：

基于插值DFT的同步相量数据次同步振荡参数辨识法，对相量测量单元获取的相量数据进行频谱分析，从而获得次同步振荡的参数，包括以下步骤：

步骤一、采用广域测量系统内的相量测量单元获取实时高分辨率波形数据；

步骤二、在相量测量单元中对时域波形数据进行傅里叶变换得到同步相量数据，采用上报频率f_r，即以间距f_pr＝f_p/f_r对同步相量数据进行重新采样，在主站获得上报的同步相量X_r(m)，f_p为相量测量单元对时域波形数据的采样频率；

步骤三、在定窗宽为Nr的矩形窗内对同步相量X_r(m)进行离散傅里叶变换操作，得到离散傅里叶变换谱F_r(k)；

步骤四、用次同步振荡相量的最大幅值谱确定次同步振荡幅度最大的频谱线k_m；

步骤五、构造指标相量[k1，k2，k3，k4]，引入二阶差分光谱的频谱F_rs(k_i)，获得二阶差分的比值R，利用泰勒级数展开的一阶项对F_rs(k_i)进行重写，利用

且k_m≈L_rs，获得简化后的二阶差分比值R；

步骤六、获得次同步振荡分量的频率f_s和次同步振荡分量的阻尼系数α_s的计算公式，计算次同步振荡分量的频率f_s和次同步振荡分量的阻尼系数α_s的值；

步骤七、对同步相量X_r(m)进行Hann窗口截取并执行离散傅里叶变换操作，用次同步振荡相量的最大幅值谱定位第k_m谱线，得到

并利用步骤六得到的f_s和α_s计算W_s(k_m)，反推得到次同步振荡分量幅值A_s的计算公式，计算次同步振荡分量幅值A_s的值。

目前，可再生能源的使用显著增加了次同步振荡现象发生的几率和严重性，次同步振荡会导致发电量的严重损失，并对系统稳定性和设备安全性构成巨大威胁，准确辨识次同步振荡的参数，对缓解次同步振荡对系统稳定性和设备安全性造成的影响至关重要。目前识别次同步振荡的方法分为时域法、频域法和时频域法，其中时域法利用电压/电流时间信号估计次同步振荡的参数，需要预先准确地确定模型阶数且信噪比高时，时域法才具有比较好的估计效果；时频域方法通过更为复杂的计算，利用希尔伯特-黄变换或变分模分解等方法将时间信号分解为多频分量再进行子序列识别，时频域法计算非常复杂；频域方法基于离散傅里叶变换(DFT)实现。在对次同步振荡参数的辨识过程中，始终需要高分辨率波形数据，本技术方案中采用广域测量系统(WAMS)的相量测量单元测量的同步相量数据，广域测量系统通过相量测量单元(PMU)的安装提供电力系统状态的在线估计，但是目前的相量测量单元只能测量频率为50Hz/60Hz的相量，因此，次同步振荡不能直接从同步相量测量，只能通过频谱泄漏分量表现，目前频谱泄漏分量中恢复次同步振荡信息的校正方法需要消耗大量的时间，同时不能识别次同步振荡的阻尼系数。本发明中，摒弃之前利用相量测量单元幅度数据来识别次同步振荡参数的方法，通过对相量测量单元的相量数据进行频谱分析来识别次同步振荡的参数，对上报频率为f_r，以间距为f_pr＝f_p/f_r对同步相量进行采样，获得同步相量X_r(m)，然后在定窗宽为Nr的矩形窗内对同步相量X_r(m)进行离散傅里叶变换操作，得到离散傅里叶变换谱F_r(k)，再利用次同步振荡相量的最大幅值谱确定次同步振荡幅度最大的频谱线k_m，再构造指标相量[k1，k2，k3，k4]，引入二阶差分光谱的频谱F_rs(k_i)，获得二阶差分的比值R，最后通过二阶差分的比值R获得次同步振荡分量的频率f_s和次同步振荡分量的阻尼系数α_s的计算公式，计算次同步振荡分量的频率f_s和次同步振荡分量的阻尼系数α_s的值，再对同步相量X_r(m)进行Hann窗口截取并执行离散傅里叶变换操作，反推得到次同步振荡分量幅值A_s的计算公式，计算次同步振荡分量幅值A_s的值，利用此方法，可以实时精准的估计次同步振荡分量的频率f_s、次同步振荡分量的阻尼系数α_s和次同步振荡分量的幅值A_s的值。

进一步的，所述步骤二获得的同步相量X_r(m)为：

X_r(m)＝[X_p(0)，…，X_p(mf_pr-1)，…，X_p((N_r-1)f_pr-1)]，m＝0，1，…，N_r-1。

进一步的，所述步骤三得到离散傅里叶变换谱F_r(k)的公式为：

其中F_r1是对应于基频分量的频谱，F_rs是对应次同步振荡分量的频谱。

进一步的，所述步骤五获得的二阶差分比值R公式为：

利用泰勒级数展开的一阶项对F_rs进行重写，利用

且k_m≈L_rs，

获得简化后的二阶差分比值R公式为：

其中

其中

进一步的，所述步骤六中获得的次同步振荡分量的频率f_s和次同步振荡分量的阻尼系数α_s的计算公式分别为：

其中

Im(·)和Re(·)分别表示复变量的虚部和实部。

进一步的，在确定次同步振荡分量的频率f_s和次同步振荡分量的阻尼系数α_s之后，使用Hann窗，定义ω(m)为：

将采集的同步相量谱写为：

其中

和

分别为基频分量和次同步振荡分量的频谱；上式中的

重写为：

其中

其中W_s(k)为阻尼Hann窗的频谱，定义为：

其中，D为阻尼矩形窗的频谱；

由此次同步振荡振幅A_s可以由下式求得：

其中k_m是

峰值的编号。

使用矩形窗可以满足对f_s和α_s进行估计的性能要求，一旦确定α_s和f_s，可以在矩形窗内直接计算幅值A_s，但是此种方式会导致计算的幅值A_s存在较大误差。本发明中使用Hann窗对幅值A_s进行识别，通过Hann窗，定义ω(m)，然后谱写采集的同步相量谱

定义阻尼Hann窗的频谱W_s(k)，获得次同步振荡振幅A_s公式，利用公式计算次同步振荡振幅A_s的取值，采用Hann窗对幅值A_s进行精确的识别。

进一步的，所述阻尼矩形窗的频谱D由以下公式获得；

其中

为频谱中归一化后的阻尼系数。

本发明与现有技术相比，具有如下的优点和有益效果：

1、本发明中，摒弃之前利用相量测量单元获取的同步相量的幅值数据来识别次同步振荡参数的方法，通过对相量测量单元的复数量的相量数据进行频谱分析来识别次同步振荡的参数，对上报频率为f_r，以间距为f_pr＝f_p/f_r对同步相量进行采样，获得同步相量X_r(m)，然后在定窗宽为Nr的矩形窗内对同步相量X_r(m)进行离散傅里叶变换操作，得到离散傅里叶变换谱F_r(k)，再利用次同步振荡相量的最大幅值谱确定次同步振荡幅度最大的频谱线k_m，再构造指标相量[k1，k2，k3，k4]，引入二阶差分光谱的频谱F_rs(k_i)，获得二阶差分的比值R，最后通过二阶差分的比值R获得次同步振荡分量的频率f_s和次同步振荡分量的阻尼系数α_s的计算公式，利用此方法，可以实时精准的分析和估计次同步振荡分量的频率f_s和次同步振荡分量的阻尼系数α_s。

2、本发明中使用Hann窗对幅值A_s进行识别，通过Hann窗，定义ω(m)，然后谱写采集的同步相量谱

定义阻尼Hann窗的频谱W_s(k)，获得次同步振荡振幅A_s公式，利用公式计算次同步振荡振幅A_s的取值，采用Hann窗对可以对次同步振荡分量幅值A_s进行精确的识别。

附图说明

此处所说明的附图用来提供对本发明实施例的进一步理解，构成本申请的一部分，并不构成对本发明实施例的限定。在附图中：

图1为基于插值DFT的同步相量数据次同步振荡参数辨识法流程图；

图2为相量测量单元所采集的同步相量的频谱；

图3为使用矩形窗和Hann窗得到同步相量的频谱；

图4为示例性测试信号和相应同步相量的幅值图；

图5为图4fs在[10，45]Hz范围内变化时的估计误差图；

图6为对于不同的基频f₁，采集的具有相同的A₁，A_s，f_s，α_s和φ_s值的同步相量的频谱幅值图；

图7为使用标称和非标称频率情况下α_s函数的估计误差图

图8为A_s函数无测量噪声时的估计误差图；

图9为f_s在[10,45]Hz范围内变化时对于噪声瞬时信号的估计误差图；

图10为α_s在[0.01,1]范围内变化时对噪声瞬时信号的估计误差图；

图11为A_s函数在40dB信噪比测量噪声下的估计误差图；

图12为情况1的模拟结果图；

图13为情况2的模拟结果图；

图14为情况3的模拟结果图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施例和附图，对本发明作进一步的详细说明，本发明的示意性实施方式及其说明仅用于解释本发明，并不作为对本发明的限定。

实施例1：

如图1所示，基于插值DFT的同步相量数据次同步振荡参数辨识法，对相量测量单元获取的相量数据进行频谱分析，从而获得次同步振荡的参数，包括以下步骤：

且k_m≈L_rs，获得简化后的二阶差分比值R；

优选的，所述步骤二获得的同步相量X_r(m)为：

优选的，所述步骤三得到离散傅里叶变换谱F_r(k)的公式为：

优选的，所述步骤五获得的二阶差分比值R公式为：

利用泰勒级数展开的一阶项对F_rs进行重写，利用

且k_m≈L_rs，

获得简化后的二阶差分比值R公式为：

其中

其中

优选的，所述步骤六中获得的次同步振荡分量的频率f_s和次同步振荡分量的阻尼系数α_s的计算公式分别为：

其中

Im(·)和Re(·)分别表示复变量的虚部和实部。

优选的，在确定次同步振荡分量的频率f_s和次同步振荡分量的阻尼系数α_s之后，使用Hann窗，定义ω(m)为：

将采集的同步相量谱写为：

其中

和

分别为基频分量和次同步振荡分量的频谱；上式中的

重写为：

其中

其中W_s(k)为阻尼Hann窗的频谱，定义为：

其中，D为阻尼矩形窗的频谱；

由此次同步振荡振幅A_s可以由下式求得：

其中k_m是

峰值的编号。

优选的，所述阻尼矩形窗的频谱D由以下公式获得；

其中

为频谱中归一化后的阻尼系数。

本实施例中对存在次同步振荡的同步相量进行了分析，在不失一般性的前提下，基频分量和次同步分量组成的瞬时信号x(t)建模为：

其中A₁、f₁和φ₁分别是基频分量的振幅，频率和初始相位，f_s和φ_s分别是次同步振荡的频率和初始相位，α_s是次同步振荡的阻尼系数。假设相量测量单元对采样频率为f_p固定的信号进行采样，则第n次采样时(1)中的x(t)可表示为：

其中f_p1＝f₁/f_p是基频分量相对于相量测量单元采样频率f_p的归一化频率，α_ps＝α_s/f_p是次同步振荡相对于f_p的归一化阻尼系数，f_ps＝f_s/f_p是次同步振荡分量相对于f_p的归一化频率。

通过在公式(2)上配合矩形窗应用离散傅里叶变换来获得同步相量，矩形窗口的一个单周期长度为N_p，其中N_p＝f_p/f₀，f₀为电力系统正常运行的标称频率，可取50Hz或60Hz。因此，归一化基频f_p1可以改写为：

其中L₁＝L₀+σ₀是在频谱区间表示的归一化基频，σ₀是频间偏移位置。当f₁等于归一化频率f₀时(L₁＝L₀＝1,σ₀＝0)，信号中基频分量与第二个谱线相关。否则，L₁不是整数且σ₀≠0，此时说明存在频谱泄露。

同样归一化后次同步振荡的α_ps和f_ps分量可以被表示为：

其中

和L_s分别为在频谱中归一化后的阻尼和频率，其中f_s<f₀且L_s<1。

根据已定义的频谱，可以得到第r个滑动时间窗处的次同步振荡谱如下：

其中r为样本数，r＝1,2,...,k＝0,1,...,N_p-1，k为谱数，

公式(5)中右侧第一项为基频分量的频谱，将其记为X₁。根据欧拉方程和级数求和法，X₁的表达式可改写为如公式(6)所示：

其中

其中X₁ ⁺和X₁ ^-分别表示基频分量的正谱和负谱，(.)^*表示共轭算子。

公式(5)右侧第二项为次同步振荡分量的频谱，记为X_s。同样，X_s可以简化为：

其中

其中X_s ⁺和X_s ^-分别表示次同步振荡分量的正谱和负谱。

相量测量单元在第2个谱线(即k＝1)处采集的同步相量信号对应于在第r个时间窗处的频谱X(r，1)。为了简化分析，当f_p>>f₀时(或当L_p>>1)时，负谱图像(X₁ ^-+X_s ^-)对(X₁ ⁺+X_s ⁺)的影响可以被忽略，因此将同步相量X_p(r)简化为：

其中

从公式(10)可以看出，X_p不仅包含基频分量构成的相量，还包含次同步振荡构成的相量。本文的目的是从Xp中估计次同步振荡的参数。现有技术中利用同步相量的幅值|X_p(r)|来实现这一目的。

虽然|X_p(r)|中含有A_s的信息，但由于存在A₁的影响且需要进行非线性运算，很难直接从公式(13)中计算A_s。因此，需要进行大量的离线实验才能找到A_s、A₁与|X_p(r)|之间的数值关系。事实上，相量X_p本身提供了比|X_p(r)|更详细的次同步振荡信息。因此，本实施例中基于公式(10)中的相量X_p识别f_s、A_S和φ_s等参数。

本实施例中由相量测量单元提供的同步相量以确定的上报频率传输到主站，对于50Hz系统，通常采用50Hz或100Hz的频率。假设上报频率为f_r，对同步相量X_p(r)以f_pr＝f_p/f_r的间距重新采样，如公式(14)所示：

X_r(m)＝[X_p(0),…，X_p(mf_pr-1),…，X_p((N_r-1)f_pr-1)]

(m＝0，1，…，N_r-1) (14)

对相量测量单元上报的同步相量X_r(m)执行离散傅里叶变换运算，由于进行重新采样，公式(10)中的变量r应被替换为mf_pr，可得：

rf_p1＝mf_r1，rf_ps＝mf_rs，rα_ps＝mα_rs (15)

其中

上式中的各变量分别为对应f_r的归一化频率和阻尼系数。综上，对相量测量单元提供的同步相量X_r(m)在定窗宽为N_r的矩形窗内进行的离散傅里叶变换谱可以确定为：

其中F_r1是对应于基频分量的频谱，F_rs是对应次同步振荡分量的频谱。类比于公式(4)，公式(16)中的α_rs和f_rs可以在频谱中重写为：

其中，当f_s与f_r不一致时，L_rs不为整数，则F_rs为：

其中D为阻尼矩形窗的频谱，由公式(20)求得：

考虑到L_rs可能不是整数，应该更严格地表示F_rs如公式(21)所示。

式中LX_r1(k)为基分量的光谱泄漏量。当f_s接近f₁时，LX_r1(k)对F_rs的影响增大。解决这一问题的有效方法之一是使用更宽的窗口长度N_r对相量进行离散傅里叶变换分析。这是因为较大的N_r提高了离散傅里叶变换的频率分辨率，使得一些非相干频率f₁能位于整数谱线上。根据发明者的测试，2秒的数据足以最小化LX_r1(k)的影响。此外由于f_s与f_r也不一致，F_rs出现了如图2所示的栅栏效应。次同步振荡频谱的谱线(图中点线)不位于整数频线上。换言之，次同步振荡的实际谱线位于k_m-1和k_m或k_m和k_m+1的频谱线之间，其中第k_m为次同步振荡幅度最大的频谱线。

为了精准确定f_s和α_s，本实施例中引入了二阶差分光谱的F_rs(k_i)。定义二阶差分的比值R为

其中谱线数向量[k1,k2,k3,k4]由公式(23)给出

公式(19)中F_rs的分母部分可以用泰勒级数展开的一阶项重写。考虑到

且k_m≈L_rs，第k_m个谱线处的频谱F_rs可简化为：

其中

据此，公式(22)中F_rs(k₁)、F_rs(k₂)和F_rs(k₃)的二阶差分可确定为：

公式(26)的推导假设Lx_r1(k_i)(i＝1,2,3)彼此之间存在较小的差异。由于f_s通常远离f₁，这个假设是有效的，因此使用二阶差分可以有效地减少Lx_r1对ω_rs的影响。同样，公式(22)中F_rs(k₂)、F_rs(k₃)和F_rs(k₄)的二阶差分为：

将公式(26)、(27)代入公式(22)，则比值R由下式计算：

因此，ω_rs可以通过计算公式(28)得到，R由公式(22)给出，F_rs直接从被提取的次同步振荡相量的离散傅里叶变化谱中获得。最后基于公式(18)和公式(25)，f_s频率和阻尼系数α_s可以由公式(29)确定。

其中

式中，Im(·)和Re(·)分别表示复变量的虚部和实部。

最后根据公式(30)和公式(29)，计算α_s和f_s。

一旦确定α_s和f_s，可以方便地根据公式(26)和(27)直接计算幅值A_s，然而这样的解决方案会导致很大的误差。原因是公式(26)和公式(27)中多个α_rs二阶差分虽然对f_s的估计影响不大，但对A_s的估计有较大的影响，为了解决这一问题，本文中利用公式(19)中的零阶差分计算A_s，即使用公式(19)中的F_rs(k_m)计算A_s，其中F_rs(k_m)对应于次同步振荡相量的最大幅度谱。

其中C_s、ω_rs分别由公式(12)和公式(25)给出。

如前所述，当f₁与采样频率f_p不相干时，f₁基相量的频谱泄漏会影响公式(31)中F_rs(k_m)的大小，如图3(a)所示。在本实施例中，使用了一个2秒的数据窗口，使频谱泄露对频率计算的影响最小化。然而，这种数据截取方法并不足以计算振幅的精确性。为了解决这一问题，进一步采用Hann窗来减小次同步振荡相量上基相量的频谱泄漏现象，如图3(b)所示。理论上，f_s的估计也可以使用Hann窗进行。然而，Hann窗截取后次同步振荡相量谱的表达式过于复杂，(公式(34)-(35)中的

)，使得计算公式(22)中的R较繁琐，其次使用矩形窗可以满足对f_s和α_s进行估计的性能要求，因此Hann窗仅用来进行A_s的识别，在图3中，其中A₁＝100，f₁＝49.7Hz，A_s＝10，f_s＝44.8Hz，α_s＝0.125，f_p＝10KHz，fr＝100Hz，(a)矩形窗，(b)Hann窗。

通过使用Hann窗，将ω(m)定义为：

所采集的同步相量谱写为：

其中

和

分别为基频分量和次同步振荡分量的频谱。公式(33)中的

重写为：

其中W_s(k)为阻尼Hann窗的频谱，定义为：

其中，D由公式(20)算得。

f_rs和

分别通过计算公式(29)和公式(30)得到，至此除A_s之外所有参数均已知，因此A_s可以由下式求得：

其中k_m是

峰值的编号。

本实施例中f_s为次同步振荡分量的实际频率；

α_s为次同步振荡分量的实际阻尼系数；

A_s为次同步振荡分量的实际幅值；

为次同步振荡分量频率的估计值；

为次同步振荡分量阻尼系数的估计值；

为次同步振荡分量幅值的估计值；

f₁为基频分量频率；

A₁为基频分量幅值。

为了验证本方法的可行性，本实施例中使用合成数据和模拟的次同步振荡数据进行验证研究：

本发明在合成信号中的应用，

合成信号建模与公式(1)相同，其参数定义如下：A1固定为100，同时考虑标称频率(f₁＝50Hz)和非标称频率(f₁＝49.5、49.7和50.3Hz)。fs变化范围为[10，45]HZ，频率变换间隔为0.5Hz，A_s变化的范围为[5％、50％]×A₁，变化间隔为5％×A₁，α_s变化的范围[0.01，1]。为了验证算法参数估计的有效性，100条满足上述条件且初相以π/10的间隔从-π变化到π的信号被用作实验数据。

实验数据的采样频率为10KHz，对N_p＝200(一个周期)的瞬时时间数据进行DFT算法。提取的同步相量最终以100Hz的上报频率传输到主站。如图4所示为示例性测试信号和相应同步相量的幅值，图4测试信号示例(a)时域中的正弦信号(f₁＝49.7Hz，A₁＝30A，f_s＝10.4Hz，A_s＝10A，α₁＝0，α_s＝0.25，φ_s＝π/6，f_p＝10KHz)，(b)PMU以100Hz上报率所采集相量的幅值可达到10％。将本发明提出的方法每次应用于200个同步相量。为了评估计算结果的准确性，选取相对误差的最大值作估计误差(EE)为，定义为：

其中P与

分别为被估计参数的真实值和估计值。

1)f_s的影响：当A₁＝100，A_s＝10且α_s＝0.25时，不同f_s时对三个参数的估计误差如图5所示。由图可知，当f₁＝50Hz和49.5Hz时，估计误差小于f₁＝49.7和50.3Hz时的估计误差。原因是当Nr＝200个同步相量用于执行离散傅里叶变换时(例如频率分辨率为0.5Hz)，此时f₁＝49.5Hz和f₁＝50Hz的同步相量的频谱位于整数频谱线中。这意味着f₁＝49.5Hz的频谱泄漏影响小于其对f₁＝49.7和50.3Hz的影响，如图6所示。此外，图5表明当f_s越接近f₁时，由于基本相量的频谱泄漏，估计误差越大。实验证明，当f_s在[10，40]Hz范围内时，三个参数的估计误差均低于1％。图5为在[10，45]Hz范围内变化时的估计误差图，无测量噪声下，次同步振荡分量的频率fs在[10,45]Hz范围内变化时，估计得到的参数值与真实值之间的最大估计误差，其中A₁＝100，A_s＝10且α_s＝0.25，(a)为频率误差，(b)为阻尼误差，(c)为幅度误差；图6为对于不同的基频f₁，所采集具有相同的A₁，A_s，f_s，α_s和φ_s值的同步相量的频谱幅值图，(a)f₁＝50.0Hz，(b)f₁＝49.5Hz，(c)f₁＝49.7Hz，(d)f₁＝50.3Hz。

2)α_s的影响：当f_s＝30.25Hz且A_s＝0.1×A₁时，不同α_s下三个参数的估计误差如图7所示。当f_s＝30.25Hz时，对次同步振荡相量进行非相干采样，如图7(a)所示，不同的α_s对标称或非标称频率下的f_s和A_s的估计结果影响很小，而α_s的估计误差随着α_s的增加而逐渐减小。原因是较大的α_s有助于突出次同步振荡相量的频谱(参考公式(10)中的X_s ⁺(r，1))从而减小了对α_s的估计误差；图7为使用标称和非标称频率情况下α_s函数的估计误差图，无测量噪声下，次同步振荡分量的α_s在[0.01,1]范围内变化时，估计得到的参数值与真实值之间的最大估计误差，f_s＝30.25Hz，A_s＝0.1×A₁＝10，φ_s在[-π，π]变化。(a)频率误差(b)阻尼误差(c)幅度误差。

3)A_s的影响：当f_s＝10.75Hz，α_s＝0.25时，不同A_s下三个参数的估计误差如图8所示。从图中可以看出，即使在A_s很小的情况下，f_s和α_s的估计结果依然准确。如果A_s小到0.05×A₁，则A_s的估计误差可达到10％。图8为A_s函数无测量噪声时的估计误差图，次同步振荡分量的A_s在[0.05,0.5]A₁范围内变化时，估计得到的参数值与真实值之间的最大估计误差，其中A₁＝100，f_s＝10.75Hz，α_s＝0.25，φ_s在-π和π之间变化。(a)频率误差(b)阻尼误差(c)幅度误差。

4)噪声的影响：为了研究该方法在噪声条件下的性能，在测试信号中加入具有40dB信噪比(SNR)的零均值高斯噪声后重新验证算法效果，定义SNR为

σ²是噪声的方差。结果如图9-图11所示。可以看出，在有噪声的情况下所有参数的估计仍然具有可接受的准确度。然而，当α_s或A_s非常小时，它们的估计准确度可能容易受到噪声的影响。图9为f_s在[10,45]Hz范围内变化时对于噪声瞬时信号的估计误差图，其中A₁＝100，A_s＝10,α_s＝0.25,φ_s在[-π，π]变化。(a)频率误差(b)阻尼误差(c)幅度误差；图10为α_s在[0.01，1]范围内变化时对噪声瞬时信号的估计误差，其中f_s＝30.25Hz，A_s＝0.1×A₁＝10，φ_s在[-π，π]变化，(a)频率误差(b)阻尼误差(c)幅度误差；图11为A_s函数在40dB信噪比测量噪声下的估计误差图，其中A₁＝100，f_s＝10.75Hz，α_s＝0.25，φ_s在[-π，π]变化。(a)频率误差(b)阻尼误差(c)幅度误差。

5)Hann窗口对A_s估计的重要性：表I显示了使用矩形窗口同步相量估算A_s时的结果。通过比较，可以得出结论，比起使用矩形窗口，使用Hann窗口进行A_s估计会大大减小估计误差。

表1为Hann和矩形窗口估计A_s的最大误差对比表，A1＝100，A_s＝10，α_s＝0.25。

表1 Hann和矩形窗口估计A_s的最大误差对比表

本发明在模拟数据中的应用

为了通过模拟的次同步振荡数据进一步测试本技术方案。在基于Matlab/Simulink平台的风电场系统里建模。选取感应电机的六阶模型，用双质量传动系模型表示发电机轴。表2中提供了关键系统参数，风速设定为9米/秒。

表2系统参数

系统电压	220KV
		等效系统电抗	19.98Ω
导线电感，L	0.3H
		导线电阻，R	5.3Ω
串联电容，C(20％补偿)	110μF
		变压器电压	220kV/690V

当线路的串联补偿电平从线路电感的20％增加到30％时，次同步振荡事件在t＝1s时启动。该事件是由传输线的串联电容器与系统其余部分的电感之间的电谐振引起的。共振频率(s_ω)下的DFIG滑差为负。这种负滑差使DFIG转子的等效电阻也为负；当该负电阻的幅度大于系统其余部分的电阻(DFIG定子和网络电阻)时，谐振变得不稳定。图12为情况1的模拟结果，如图12(a)所示为通过线路的A相电流的波形，图12(b)示出了相应的A相电流的PMU数据。

选取2s到4s之间的相量测量数据作为验证本文方法的输入数据。获得的次同步振荡信息如表3所示。将结果与以波形数据作为输入的Prony算法计算的结果进行比较，从表3中可以看出，两种方法得到的结果是一致的，表明本文所提出的方法可以准确地从相量测量单元数据中提取次同步振荡信息。

表3不同补偿水平和风速下SSO参数估计结果

本实施例中模拟了另外两种情况(案例2和案例3)以进一步验证。在这两种情况下，风速为10米/秒。线路的补偿水平分别从30％增加到35％和40％。如图13和图14所示，谐振在情况2中衰减，但在情况3中迅速增长，表明系统的补偿水平对系统稳定性起着重要作用。验证结果如表3所示。同样，通过本文算法估计的次同步振荡信息与通过Prony算法获得的次同步振荡信息高度一致。图13为情况2的模拟结果，其中(a)来自A相的电流波形，(b)A相电流的PMU数据。图14为情况3的模拟结果，其中(a)来自A相的电流波形，(b)A相电流的PMU数据。

综上所述，本实施例中验证了与现有技术相比，该方法更加有效和准确。仿真结果表明，即使在非标称和噪声条件下，所提出的方法也可以达到可接受的精度。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。