CN112688325B - 基于二阶段改进itd算法的风电场次同步振荡监测方法 - Google Patents

Abstract

本发明公开了一种基于二阶段改进ITD算法的风电场次同步振荡监测方法，包括以下步骤：S1：采集风电场的电压电流信号，并进行预处理；S2：利用二阶段改进ITD算法实时识别预处理后电压电流信号中的次同步振荡参数；S3：基于次同步振荡参数和次同步振荡决策条件，实时监测风电场的次同步振荡。使用本发明中的监测方法可以实时监测电网中的次同步振荡信号，具有算法简单轻便、计算量小、检测速度快、动态性能强、抗噪声能力强且抗干扰能力强的优点。

Description

基于二阶段改进ITD算法的风电场次同步振荡监测方法

技术领域

本发明属于风电场次同步振荡监测技术领域，具体涉及一种基于二阶段改进ITD算法的风电场次同步振荡监测方法。

背景技术

随着可再生能源的快速发展，大规模风电系统接入电网，使得远距离和大规模风电外送成为常态。固定串补电容在提高系统的输电能力的同时也可能与双馈风机发生相互作用，引发次同步振荡。2011年，美国明尼苏达州风电系统的一条输电线路被切除，引发了9～13Hz的次同步振荡；2011年以来，中国河北沽源风电系统发生了多起频率在3～10Hz内波动的次同步振荡事故，致使大量风机脱网。

为了避免次同步震荡在发生后对电力系统造成危害，有必要实时监测电网中振荡信号。考虑大规模风电场运行时工况复杂，电力系统的动态变化引发次同步振荡具有随机性和强时变性，因此有必要选取一款动态跟踪性能好、抗噪声、抗干扰能力强且简单易于实现的算法，应用于电网振荡信号的监测。

目前，多种算法已被应用于次同步振荡的检测，其中主流的算法有以下几种：Prony算法可以快速得到振荡分量的频率、幅值和阻尼等模态参数，但是该算法对噪声敏感，在采样信号中存在的噪声会给算法的参数辨识结果带来误差；快速傅里叶算法FFT可以将采样的时域信号转换到频域，得到不同频率分量的信息，但是该算法依赖一个长窗来提高结果的分辨率，带来了较大的计算量并且牺牲了检测的实时性；时域信号分解算法有经验模态分解算法EMD和变分模态分解VMD算法，可以从采样信号中辨识出振荡的关键参数，但是这一类算法具有复杂的数学运算带来较大计算量。子空间类算法有ERA、MPM和TLS-ESPRIT等，这一类算法可以以较高分辨率提取出采样信号中的振荡分量，但是由于新型次超同步振荡存在的多频率分量可能性，算法的应用存在定阶困难的问题。由于算法无法自适应预先准确选择被测信号的模态数，算法的阶数与被测信号的模态数不匹配会使得振荡分量跟踪不准确。目前辅助算法定阶已有较多方法，这些辅助定阶方法很难发挥理想的作用。综上所述，现有振荡检测算法或多或少存在数据处理复杂、动态性能差和定阶困难等问题，难以在实际工程中开发和应用。

发明内容

本发明的目的是为了解决现有振荡检测算法存在数据处理复杂、动态性能差和定阶困难的问题，提出了一种基于二阶段改进ITD算法的风电场次同步振荡监测方法。

本发明的技术方案是：一种基于二阶段改进ITD算法的风电场次同步振荡监测方法包括以下步骤：

S1：采集风电场的电压电流信号，并进行预处理；

S2：利用二阶段改进ITD算法实时识别预处理后电压电流信号中的次同步振荡参数；

S3：基于次同步振荡参数和次同步振荡决策条件，实时监测风电场的次同步振荡。

本发明的有益效果是：

(1)使用本发明中的监测方法可以实时监测电网中的次同步振荡信号，具有算法简单轻便、计算量小、检测速度快、动态性能强、抗噪声能力强且抗干扰能力强的优点。

(2)本发明的监测方法考虑新型次超同步振荡存在的多频率分量可能性，无需预先估计被测信号的模态，使用方便。并且，本发明的监测方法可以直接嵌入DSP控制器，易于硬件实现且检测效果好。

(3)使用发明中所开发的次同步振荡监测方法，可以快速有效提取风电场关键节点的输出电压电流中的次同步分量，在次同步振荡发生后及时发出预警信号，可以为大规模风电场次同步振荡的检测、溯源和抑制提供有效的关键信息。

进一步地，步骤S1中，对电压电流信号进行预处理的方法为：采用8阶低通FIR滤波器进行滤波处理，其截止频率是120Hz。

上述进一步方案的有益效果是：在本发明中，信号预处理的目的是为了滤除采样电压电流信号中不需要的噪声和高次谐波。

进一步地，步骤S2包括以下子步骤：

S21：利用改进ITD算法计算预处理后电流信号的次同步分量频率f_sub和电压信号的基波分量频率f₀，完成一阶段频率估计；

S22：根据电流信号的次同步分量频率f_sub和电压信号的基波分量频率f₀，获取电压电流信号中每个模态的相量，完成二阶段次同步振荡的参数识别。

上述进一步方案的有益效果是：在本发明中，针对ITD算法在振荡监测应用中的缺点，本发明提出将算法结合最小二乘法做出二阶段实现的改进。

进一步地，步骤S21包括以下子步骤：

S211：根据预处理后电流信号中振荡分量的关键点计算过零点(t_z,L_z)；

S212：根据过零点(t_z,L_z)进行频率值更新，并根据相邻过零点之间的距离计算电流信号的次同步分量频率f_sub；

S213：根据预处理后电压信号计算电压信号的基波分量频率f₀，完成一阶段频率估计。

上述进一步方案的有益效果是：在本发明中，为了估计振荡分量的频率，由于振荡分量在电流信号中的振荡更明显，本发明采用电流信号提取，对频率的估计结果会更精确。改进后的ITD算法在第一阶段提供次同步分量和基波的频率。

进一步地，步骤S211中，电流信号中振荡分量的关键点的计算公式为：

其中，L_k+1表示振荡分量的关键点，X_k、X_k+1和X_k+2分别表示原始信号在给定区间内的三个连续极值点；τ_k、τ_k+1和τ_k+2分别表示X_k、X_k+1和X_k+2所对应的时刻；

步骤S211中，过零点(t_z,L_z)的计算方法为：根据满足L_k+1*L_k<0的两个相邻极值点确定过零点(t_z,L_z)，其中，过零点(t_z,L_z)中纵坐标L_z的计算公式为：

/>

其中，L_k表示振荡分量的关键点，X_z表示过零点对应的原始信号，

根据过零点(t_z,L_z)所对应的原始信号在两个相邻极值点之间的单调性，得到过零点(t_z,L_z)的横坐标t_z。

上述进一步方案的有益效果是：在本发明中，振荡分量的关键点由原始信号的极值点决定，对于连续满足L_k+1*L_k<0的两个点，必定存在一个过零点；由于原始信号在两个相邻极值点之间必定单调，由此可以得到X_z的横坐标，即过零点(t_z,L_z)的横坐标t_z。

进一步地，步骤S212中，电流信号的次同步分量频率f_sub的计算公式为：

f_sub＝1/(2dt₁*(t_z-t_z-1))

其中，dt₁表示数据的采样时间间隔，t_z表示电流信号在检测过程中最新更新的振荡分量的过零点的位置，t_z-1表示电流信号在检测过程中最新更新的振荡分量的相邻过零点的位置。

上述进一步方案的有益效果是：在本发明中，次同步分量的频率在低频振荡信号出现新的过零点时更新。

进一步地，步骤S213中，电压信号的基波分量频率f₀的计算公式为：

f₀＝1/(2dt₁*(tv_i-tv_i-1))

其中，tv_i表示电压信号最新更新的过零点位置，tv_i-1表示电压信号最新更新的相邻过零点位置，dt₁表示数据的采样时间间隔。

上述进一步方案的有益效果是：在本发明中，在计算基波频率时，由于在电能传输过程中电压受到干扰的影响更小，因此基波频率f₀根据采样的电压信号计算得到更准确。

进一步地，步骤S22包括以下子步骤：

S221：根据电流信号的次同步分量频率f_sub和电压信号的基波分量频率f₀，计算超同步分量频率f_sup并构建电压电流信号的特征值向量λ；

S222：根据特征值向量λ构建电压电流信号的范德蒙矩阵A；

S223：利用最小二乘法，根据范德蒙矩阵A得到电压电流信号中每个模态的相量，完成次同步振荡参数识别。

上述进一步方案的有益效果是：在本发明中，为了解决幅值辨识中出现的问题，在算法第二阶段结合了最小二乘法估计振荡分量的幅值参数。在本发明中，电压电流采样信号中使用10个点的数据在最小二乘法中进行参数估计，这使得次同步分量的辨识连续，同时，参数的辨识受到噪声的影响更小。

进一步地，步骤S221中，超同步分量频率f_sup的计算公式为：

f_sup＝2*f₀-f_sub

其中，f_sub表示电流信号的次同步分量频率，f₀表示电压信号的基波分量频率；

电压电流信号的特征值向量λ的计算公式为：

其中，ω_i＝2πf_i，λ_n,n＝1,2,...,2i表示对应于所有q个模态的q对共轭的特征值，T表示矩阵转置运算，q表示模态个数，σ_i表示所分析信号中第i个分量的阻尼，f_i表示所分析信号中第i个分量的频率，ω_i表示表示对应于频率f_i的角频率，j表示虚部单位；

步骤S222中，范德蒙矩阵A的表达式为：

其中，λ^p _q表示特征值向量中第i个元素的p次幂，p＝0,1,...,m-1，q＝1,2,...,2q，m表示电压电流号中利用的点的数量；

步骤S223中，电流信号中每个模态的相量C_I的计算公式为：

C_I＝(A^TA)^-1A^TY_I

其中，Y_I＝[p_I1,p_I2,p_I3,...,p_Im]^T表示在最小二乘法拟合的过程中所使用的电流信号的点组成的相量，p_Ii,Ii＝I1,I2,...,Im表示在最小二乘法拟合的过程中所使用的所有电流数据，

电压信号中每个模态的相量C_V的计算公式为：

C_V＝(A^TA)^-1A^TY_V

其中，Y_V＝[p_V1,p_V2,p_V3,...,p_Vm]^T表示在最小二乘法拟合的过程中所使用的电压信号的点组成的相量，p_Vi,Vi＝V1,V2,...,Vm表示在最小二乘法拟合的过程中所使用的所有电压数据。

上述进一步方案的有益效果是：在本发明中，由于σ_i＜＜ω_i，阻尼系数对特征值的贡献可以忽略不计，因此有

由此，范德蒙矩阵A∈C^m×2q可以由特征值λ构造而成。

进一步地，步骤S3中，根据次同步振荡决策条件判断风电场发生次同步振荡的方法为：若在60Hz系统中，振荡频率在5-55Hz范围内，则风电场发生次同步振荡，否则风电场未发生次同步振荡；

若基波幅值的10％大于预警阈值，则风电场发生次同步振荡，否则风电场未发生次同步振荡；

若最小检测时间大于设定时间阈值，则风电场发生次同步振荡，否则风电场未发生次同步振荡。

附图说明

图1为风电场次同步振荡监测方法的流程图；

图2为ITD估计幅值误差原理图；

图3为四种对比算法定三阶和改进ITD算法对二阶模型Y₁和Y₃的频率、幅值辨识结果图；

图4为四种算法定二阶和改进ITD算法对三阶模型Y₂的频率、幅值辨识结果图；

图5为五种算法在35dB噪声下的频率和幅值估计结果图；

图6为在基波频率随时间变化时五种算法对次同步分量幅值和频率的跟踪性能图；

图7为在次同步频率随时间变化时五种算法对次同步分量幅值和频率的跟踪性能图；

图8为在基波幅值阶跃变化时五种算法对次同步分量幅值和频率的跟踪性能图；

图9为在次、超同步分量幅值阶跃变化时五种算法对次同步分量幅值和频率的跟踪性能图；

图10为ERCOT地区风电场建模示意图；

图11为案例1中仿真电压电流波形图；

图12为案例1中各算法对次同步分量的频率和幅值的跟踪性能图；

图13为案例2中仿真电压电流波形图；

图14为案例2中各算法对次同步分量的频率和幅值的跟踪性能图；

图15为五种算法对河北沽源实测数据的频率、幅值的估计结果图；

图16为五种算法对新疆哈密实测数据的频率、幅值的估计结果图；

图17为硬件在环实验电压电流波形图；

图18为硬件在环实验中ITD算法对次同步分量的频率和幅值的估计结果图；

图19为ITD分解原理图。

具体实施方式

下面结合附图对本发明的实施例作进一步的说明。

在描述本发明的具体实施例之前，为使本发明的方案更加清楚完整，首先对本发明中出现的缩略语和关键术语定义进行说明：

ITD：Intrinsic Time-scale Decomposition，固有时间尺度分解算法；

ITD算法可以自适应分解信号成多个振荡分量。其提取原理是，通过线性插值拟合原始信号的包络线，利用信号极值点附近上下两条包络线确定信号内振荡分量的关键点，通过插值拟合关键点即可得到低频振荡信号，原始信号减去该低频信号可以得到高频振荡信号。ITD算法简单、运算速度快、时频分辨率高，而且该算法不需要任何预先的模态信息。但同时，该算法应用于次同步振荡监测时具有以下缺点：

(1)振荡的幅值由算法所提取的波形的峰值点决定，也就是说振荡的能量只能够在监测到遮挡波形峰值的时候才能进行评估，这使得振荡的监测不连续。

(2)所提取的振荡信号由原始信号的极值点通过线性插值确定，这在某些情况下可能会产生很大的误差。如图2所示，算法所拟合的峰值点在K点，而实际波形的峰值在P点。

(3)在发生次同步振荡时，超同步振荡有可能伴随发生，而其频率范围内的分量对50Hz或60Hz的基波波形的极值点影响较小，这导致ITD算法难以检测到超同步振荡的发生。

(4)次同步振荡的频率一般在10～40Hz之间，要监测一个完整的周期，至少需要存储100ms的数据，这会占用大量的硬件资源。

如图1所示，本发明提供了一种基于二阶段改进ITD算法的风电场次同步振荡监测方法，包括以下步骤：

S1：采集风电场的电压电流信号，并进行预处理；

S2：利用二阶段改进ITD算法实时识别预处理后电压电流信号中的次同步振荡参数；

S3：基于次同步振荡参数和次同步振荡决策条件，实时监测风电场的次同步振荡。

在本发明实施例中，如图1所示，步骤S1中，对电压电流信号进行预处理的方法为：采用8阶低通FIR滤波器进行滤波处理，其截止频率是120Hz。

在本发明中，信号预处理的目的是为了滤除采样电压电流信号中不需要的噪声和高次谐波。

在本发明实施例中，如图1所示，步骤S2包括以下子步骤：

S21：利用改进ITD算法计算预处理后电流信号的次同步分量频率f_sub和电压信号的基波分量频率f₀，完成一阶段频率估计；

S22：根据电流信号的次同步分量频率f_sub和电压信号的基波分量频率f₀，获取电压电流信号中每个模态的相量，完成二阶段次同步振荡的参数识别。

在本发明中，针对ITD算法在振荡监测应用中的缺点，本发明提出将算法结合最小二乘法做出二阶段实现的改进。

在本发明实施例中，如图1所示，步骤S21包括以下子步骤：

S211：根据预处理后电流信号中振荡分量的关键点计算过零点(t_z,L_z)；

S212：根据过零点(t_z,L_z)进行频率值更新，并根据相邻过零点之间的距离计算电流信号的次同步分量频率f_sub；

S213：根据预处理后电压信号计算电压信号的基波分量频率f₀，完成一阶段频率估计。

在本发明中，为了估计振荡分量的频率，由于振荡分量在电流信号中的振荡更明显，本发明采用电流信号提取，对频率的估计结果会更精确。改进后的ITD算法在第一阶段提供次同步分量和基波的频率。

在本发明实施例中，如图19所示，步骤S211中，电流信号中振荡分量的关键点的计算公式为：

其中，L_k+1表示振荡分量的关键点，X_k、X_k+1和X_k+2分别表示原始信号在给定区间内的三个连续极值点；τ_k、τ_k+1和τ_k+2分别表示X_k、X_k+1和X_k+2所对应的时刻；

步骤S211中，过零点(t_z,L_z)的计算方法为：根据满足L_k+1*L_k<0的两个相邻极值点确定过零点(t_z,L_z)，其中，过零点(t_z,L_z)中纵坐标L_z的计算公式为：

其中，L_k表示振荡分量的关键点，X_z表示过零点对应的原始信号，

根据过零点(t_z,L_z)所对应的原始信号在两个相邻极值点之间的单调性，得到过零点(t_z,L_z)的横坐标t_z。

在本发明中，振荡分量的关键点由原始信号的极值点决定，对于连续满足L_k+1*L_k<0的两个点，必定存在一个过零点；由于原始信号在两个相邻极值点之间必定单调，由此可以得到X_z的横坐标，即过零点(t_z,L_z)的横坐标t_z。

在本发明实施例中，如图1所示，步骤S212中，电流信号的次同步分量频率f_sub的计算公式为：

f_sub＝1/(2dt₁*(t_z-t_z-1))

其中，dt₁表示数据的采样时间间隔，t_z表示电流信号在检测过程中最新更新的振荡分量的过零点的位置，t_z-1表示电流信号在检测过程中最新更新的振荡分量的相邻过零点的位置。

在本发明中，次同步分量的频率在低频振荡信号出现新的过零点时更新。

在本发明实施例中，如图1所示，步骤S213中，电压信号的基波分量频率f₀的计算公式为：

f₀＝1/(2dt₁*(tv_i-tv_i-1))

其中，tv_i表示电压信号最新更新的过零点位置，tv_i-1表示电压信号最新更新的相邻过零点位置，dt₁表示数据的采样时间间隔。

在本发明中，在计算基波频率时，由于在电能传输过程中电压受到干扰的影响更小，因此基波频率f₀根据采样的电压信号计算得到更准确。

在本发明实施例中，如图1所示，步骤S22包括以下子步骤：

S221：根据电流信号的次同步分量频率f_sub和电压信号的基波分量频率f₀，计算超同步分量频率f_sup并构建电压电流信号的特征值向量λ；

S222：根据特征值向量λ构建电压电流信号的范德蒙矩阵A；

S223：利用最小二乘法，根据范德蒙矩阵A得到电压电流信号中每个模态的相量，完成次同步振荡参数识别。

在本发明中，为了解决幅值辨识中出现的问题，在算法第二阶段结合了最小二乘法估计振荡分量的幅值参数。在本发明中，电压电流采样信号中使用10个点的数据在最小二乘法中进行参数估计，这使得次同步分量的辨识连续，同时，参数的辨识受到噪声的影响更小。

在本发明实施例中，如图1所示，步骤S221中，超同步分量频率f_sup的计算公式为：

f_sup＝2*f₀-f_sub

其中，f_sub表示电流信号的次同步分量频率，f₀表示电压信号的基波分量频率；

电压电流信号的特征值向量λ的计算公式为：

其中，ω_i＝2πf_i，λ_n,n＝1,2,...,2i表示对应于所有q个模态的q对共轭的特征值，T表示矩阵转置运算，q表示模态个数，σ_i表示所分析信号中第i个分量的阻尼，f_i表示所分析信号中第i个分量的频率，ω_i表示表示对应于频率f_i的角频率，j表示虚部单位；

步骤S222中，范德蒙矩阵A的表达式为：

其中，λ^p _q表示特征值向量中第i个元素的p次幂，p＝0,1,...,m-1，q＝1,2,...,2q，m表示电压电流号中利用的点的数量；

步骤S223中，电流信号中每个模态的相量C_I的计算公式为：

C_I＝(A^TA)^-1A^TY_I

其中，Y_I＝[p_I1,p_I2,p_I3,...,p_Im]^T表示在最小二乘法拟合的过程中所使用的电流信号的点组成的相量，p_Ii,Ii＝I1,I2,...,Im表示在最小二乘法拟合的过程中所使用的所有电流数据，

电压信号中每个模态的相量C_V的计算公式为：

C_V＝(A^TA)^-1A^TY_V

其中，Y_V＝[p_V1,p_V2,p_V3,...,p_Vm]^T表示在最小二乘法拟合的过程中所使用的电压信号的点组成的相量，p_Vi,Vi＝V1,V2,...,Vm表示在最小二乘法拟合的过程中所使用的所有电压数据。

在本发明中，由于σ_i＜＜ω_i，阻尼系数对特征值的贡献可以忽略不计，因此有

由此，范德蒙矩阵A∈C^m×2q可以由特征值λ构造而成。

在本发明实施例中，如图1所示，步骤S3中，根据次同步振荡决策条件判断风电场发生次同步振荡的方法为：若在60Hz系统中，振荡频率在5-55Hz范围内，则风电场发生次同步振荡，否则风电场未发生次同步振荡；

若基波幅值的10％大于预警阈值，则风电场发生次同步振荡，否则风电场未发生次同步振荡；

若最小检测时间大于设定时间阈值，则风电场发生次同步振荡，否则风电场未发生次同步振荡。

在本发明实施例中，次同步振荡的决策由以下三个约束条件决定：

(1)在60Hz系统中，次同步振荡的频率一般在5～55Hz范围内，否则不被定义为次同步振荡，不会触发振荡监测装置的预警决策；

(2)次同步振荡幅值的大小代表了振荡的能量大小，其阈值的选取必须要在振荡监测的可靠性和预警速度之间得到平衡。在本发明中，选取基波幅值的10％作为预警阈值。在其他的保护系统中，可以根据保护需求选取合理阈值；

(3)为了避免由于意外干扰造成的误操作，需要当装置持续检测到振荡一段时间后才发出预警。在本发明中，最小检测时间选取50ms。

下面结合具体实施例对本发明的性能进行验证，其中所有算法都采用相同的1000Hz的采样率和40ms长的滑动窗口。

1.基于合成信号的测试

合成信号Y₀由次同步分量、基波分量和超同步分量建模而成：

其中，

表示基波分量的幅值和初始相位，/>

表示次同步分量的幅值和初始相位，/>

表示超同步分量的幅值和初始相位，f₀＝60Hz，A_sub＝A_sup＝A_s*10％。考虑振荡能量的逐渐增强，次、超同步分量的阻尼系数选为0.2。

1)算法的模态选择测试

考虑到大规模风电场引起的次同步振荡的随机性，有超同步振荡伴随发生的可能性。因此四种对比算法需要谨慎选择阶数，而所提出的改进ITD算法只需要在第二阶段固定q＝3即可。构造以下信号并进行测试：

Y₁＝100cos(2π*60t+10°)+10e^0.2tcos(2π*13t+20°)

Y₂＝100cos(2π*60t+10°)+10e^0.2tcos(2π*13t+20°)+10e^0.2tcos(2π*(120-13)t+30°)

如图3(a)-图3(c)所示，所有算法采用三阶时，测试二阶模型Y₁的频率和幅值的跟踪结果，可见所有算法都可以准确跟踪模型。但是当次同步频率随时间发生轻微变化时，其二阶模型为Y₃，由图3(b)-图3(d)所示，四种对比算法均不能准确跟踪振荡。可见错误的模态选择使得算法的动态跟踪能力变差。

Y₃＝100cos(2π*60t+10°)+10e^0.2tcos(2π*(13+t)t+20°)

如图4所示，当所有算法采用二阶时，测试三阶模型Y₂的跟踪性能。可见，除了改进ITD外，四种对比算法都不能准确跟踪振荡参数。

由此可见，在改进ITD算法中引入超同步分量的辨识，使得算法能够同时在二阶和三阶模型中有较好的性能。在之后的验证中，所有的信号采用三阶模型，而四种对比算法也相应取三阶。

2)抗噪声性能测试：

采样信号的噪声信噪比在40dB左右，因此，在合成信号Y4中加入35dB的白噪声以验证各算法的抗噪声能力，构造信号Y₄：

Y₄＝300cos(2π*60t+10°)+30cos(2π*13t+20°)+30cos(2π*(120-13)t+30°)

如图5所示，五种算法在35dB噪声下的频率、幅值估计图，图5(a)-图5(c)中，窗长40ms，图5((b)-图5((d)作为对照，窗长100ms。如表1所示，为五种算法在35dB噪声下的辨识误差和方差。

表1

算法	Prony	ERA	MPM	TLS-ESPRIT	ITD
						40点频率平均误差	0.3208	0.3608	0.2678	0.3722	0.1591
40点频率方差	0.3593	0.2677	0.2774	0.2753	0.0443
						40点幅值平均误差	1.522	1.6295	1.3239	1.6163	0.4084
40点幅值方差	12.1185	2.8772	6.3141	2.8212	0.3100

如图5(a)-图5(c)所示，展示了在40ms窗长下对振荡分量频率和幅值的辨识，表1展示了其结果与实际值的平均误差和方差。与四种算法相比，ITD算法在噪声环境中具有更好的鲁棒性。另外，在图5(b)-图5(d)中五种算法采用100ms的窗长作为对比。可见，更长的窗口可以提高算法的抗噪声能力，使得其结果更准确，但是由于其结果是窗长的平均值，算法的动态性能变差。

3)动态性能测试

在大规模风电系统中发生次同步震荡时，运行工况复杂多变，具有随机性与时变性，因此，次同步振荡的特性也并不是一成不变的。在振荡过程中常常出现振荡分量频率、幅值等参数的变化。选用动态性能好的算法，有助于对振荡参数有效动态跟踪和提高振荡监测装置的稳定性和可靠性。

A频率偏移

考虑基波频率随时间变化时，建立数学模型Y₅：

Y₅＝100cos(2π*(60+t)t+10°)+10e^0.2tcos(2π*13t+20°)+10e^0.2tcos(2π*(120-13)t+30°)

如图6所示，当基波频率变化时，除prony之外的四种算法都可以跟踪次同步分量的频率和幅值，受到基频变化的影响，这四种算法的跟踪都存在轻微的振荡。

当次同步分量频率随时间变化时，建立数学模型Y₆：

Y₆＝100cos(2π*60t+10°)+10e^0.2tcos(2π*(13+t)t+20°)+10e^0.2tcos(2π*(120-13-t)t+30°)

如图7所示，除了prony以外所有算法都具有较好的性能。

B幅值阶跃

考虑基波幅值发生阶跃，建立数学模型Y₇：

Y₇＝100[1-0.1ε(t-0.02)]cos(2π*60t+10°)+10e^0.2tcos(2π*13t+20°)+10e^0.2tcos(2π*(120-13)t+30°)

同样，考虑次、超同步分量发生幅值阶跃时，建立数学模型：

Y₈＝100cos(2π*60t+10°)+10e^0.2t[1-0.1ε(t-0.02)]cos(2π*13t+20°)+10e^0.2t[1-0.1ε(t-0.02)]cos(2π*(120-13)t+30°)

如图8所示，为在基波幅值阶跃变化时五种算法对次同步分量幅值和频率的跟踪性能图；如图9所示，在次、超同步分量幅值阶跃变化时五种算法对次同步分量幅值和频率的跟踪性能图；如表2所示，为基波、次超同步分量幅值阶跃变化时算法估计方差表。

表2

算法	Prony	ERA	MPM	TLS-ESPRIT	ITD
						F-Y<sub>7</sub>	409.043	3.893	11.914	12.058	0.004
A-Y<sub>7</sub>	1515.906	118.487	20.182	120.843	0.582
						F-Y<sub>8</sub>	101.875	0.140	0.312	0.225	0.000
A-Y<sub>8</sub>	646.986	0.847	1.508	0.965	0.304

其中，F-Y₇和A-Y₇分别代表了五种算法对模型Y₇的频率和幅值的估计方差，F-Y₈和A-Y₈代表了对Y₈的估计方差。

由图8和图9，各算法对次同步分量频率和幅值的跟踪结果，与四种传统算法相对比，ITD在基波和次同步分量发生阶跃变化时，具有更强的动态跟踪能力。通过表2可见，ITD算法的跟踪具有更小的方差。

2.基于仿真信号的测试

在MATLAB/Simscape电力系统电磁暂态仿真软件中建立ERCOT地区次同步振荡测试模型并验证算法的性能，如图10所示，该系统包含5个基于DFIG的风力发电场，在线路5和线路6之间放置一个串联电容器。采样在系统发生次同步振荡时的电压电流波形供算法离线测试。在这一节的测试中，关闭最小检测时间模块以观察各个算法对振荡的跟踪性能。

在案例1中，t＝0.5s时电容的串补度由25％提高到35％而引发次同步振荡。其电压电流波形和五种算法的跟踪性能如图11和12所示，表3展示了各算法检测到振荡的时间。可见所有的算法都可以准确检测并跟踪次同步振荡，而由表3，改进ITD算法具有最快的检测速度。

表3

在案例2中，t＝1s在5和8之间的线路发生三相短路故障，100ms之后故障被切除，随后发生次同步振荡，如图13所示。与案例1相比，各算法面临在故障阶段的振荡监测。可见ITD算法所辨识的结果更加接近参考值且在振荡频率落在相应范围内后才发出预警，而四种对比算法在故障切除之前给出了错误的检测结果，并且在次同步振荡开始之前就给出了预警决策，如表4所示各算法检测到次同步振荡的时间。如图14所示，案例2中各算法对次同步分量的频率和幅值的跟踪性能

表4

3.基于实测数据的测试

为了更进一步测试算法的性能，利用河北沽源和新疆哈密的次同步振荡事件中的录波数据对五种算法进行对比测试。如图15(a)-图15(c)、图16(a)-图16(c)。参考值由一个高阶的低通滤波器滤除基波后的结果表示。另外，四种对比算法选取效果最好的定阶进行比较。

由表5，与四种对比算法相比较，ITD算法具有更轻微的波动，因而具有更小的方差。通过图15(b)-图15(d)，图16(b)-图16(d)，其中图15(a)和图15(c)采用40ms窗长，图15(b)和图15(d)采用100ms窗长，图16(a)和图16(c)采用40ms窗长，图16(b)和图16(d)采用100ms窗长，各算法采用更长的100ms的窗作为前者的对比，可见在100ms的窗下受到的影响更小。因此，ITD算法在40ms的短窗中运行可以在保持精度的同时，获得更好的动态性能。表5为五种算法对沽源和哈密实测数据的频率、幅值的估计方差结果。

表5

其中，GY-40-F、GY-40-A、GY-100-F和GY-100-A分别表示了在40ms窗口和100ms窗口下测试沽源实测数据的频率和幅值的方差，HM同理代表了算法在哈密实测数据下的方差。其中，100ms窗长作为40ms窗长的对照实验。

4.硬件在环实验测试

最后，利用硬件在环实验对所提出的二阶段改进ITD算法进行测试。在搭建的硬件在环实验平台中，在实时仿真软件Typhoon中建模图10中的ERCOT风电场模型，并通过设备“Typhoon HIL 402”进行实时仿真。仿真结果中的电压电流信号通过杜邦线传输进嵌入了ITD算法的DSP控制器中。信号的采样率采用1000Hz，算法采用40ms的窗长。由于在相同条件下四种对比算法均不能在DSP中开发实现，因此，利用采样信号离线滤除基波后的结果与算法实时分析的信号对比验证。在本节中，启用最小检测时间模块来测试ITD算法对振荡的实时监测能力。

如图17所示，仿真时长1.5s，当t1＝0.549s时，接入串联电容以触发次同步振荡。如图18所示，在振荡发生18ms后t＝0.567s其能量逐渐增大并满足频率和幅值的阈值条件，设备检测到振荡信号。当最小检测时间50ms满足后，设备在t3＝0.617s发出预警信号。因此，所开发的算法可以在振荡发生的初始阶段快速检测振荡，并识别振荡参数。综上，该算法简单，易于硬件实现，便于实时监测。

本发明的工作原理及过程为：首先采集风电场的电压电流信号，并进行预处理；再利用改进ITD算法实时监测预处理后电压电流信号中的次同步振荡分量，完成风电场次同步振荡监测。

本发明的有益效果为：

(1)使用本发明中的监测方法可以实时监测电网中的次同步振荡信号，具有算法简单轻便、计算量小、检测速度快、动态性能强、抗噪声能力强且抗干扰能力强的优点。

(2)本发明的监测方法考虑新型次超同步振荡存在的多频率分量可能性，无需预先估计被测信号的模态，使用方便。并且，本发明的监测方法可以直接嵌入DSP控制器，易于硬件实现且检测效果好。

(3)使用发明中所开发的次同步振荡监测方法，可以快速有效提取风电场关键节点的输出电压电流中的次同步分量，在次同步振荡发生后及时发出预警信号，可以为大规模风电场次同步振荡的检测、溯源和抑制提供有效的关键信息。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (3)

1.一种基于二阶段改进ITD算法的风电场次同步振荡监测方法，其特征在于，包括以下步骤：

S1：采集风电场的电压电流信号，并进行预处理；

S2：利用二阶段改进ITD算法实时识别预处理后电压电流信号中的次同步振荡参数；

S3：基于次同步振荡参数和次同步振荡决策条件，实时监测风电场的次同步振荡；

所述步骤S2包括以下子步骤：

S21：利用改进ITD算法计算预处理后电流信号的次同步分量频率f_sub和电压信号的基波分量频率f₀，完成一阶段频率估计；

S22：根据电流信号的次同步分量频率f_sub和电压信号的基波分量频率f₀，获取电压电流信号中每个模态的相量，完成二阶段次同步振荡的参数识别；

所述步骤S21包括以下子步骤：

S211：根据预处理后电流信号中振荡分量的关键点计算过零点(t_z,L_z)；

S212：根据过零点(t_z,L_z)进行频率值更新，并根据相邻过零点之间的距离计算电流信号的次同步分量频率f_sub；

S213：根据预处理后电压信号计算电压信号的基波分量频率f₀，完成一阶段频率估计；

所述步骤S211中，电流信号中振荡分量的关键点的计算公式为：

其中，L_k+1表示振荡分量的关键点，X_k、X_k+1和X_k+2分别表示原始信号在给定区间内的三个连续极值点；τ_k、τ_k+1和τ_k+2分别表示X_k、X_k+1和X_k+2所对应的时刻；

所述步骤S211中，过零点(t_z,L_z)的计算方法为：根据满足L_k+1*L_k<0的两个相邻极值点确定过零点(t_z,L_z)，其中，过零点(t_z,L_z)中纵坐标L_z的计算公式为：

其中，L_k表示振荡分量的关键点，X_z表示过零点对应的原始信号，

根据过零点(t_z,L_z)所对应的原始信号在两个相邻极值点之间的单调性，得到过零点(t_z,L_z)的横坐标t_z；

所述步骤S212中，电流信号的次同步分量频率f_sub的计算公式为：

f_sub＝1/(2dt₁*(t_z-t_z-1))

其中，dt₁表示数据的采样时间间隔，t_z表示电流信号在检测过程中最新更新的振荡分量的过零点的位置，t_z-1表示电流信号在检测过程中最新更新的振荡分量的相邻过零点的位置；

所述步骤S213中，电压信号的基波分量频率f₀的计算公式为：

f₀＝1/(2dt₁*(tv_i-tv_i-1))

其中，tv_i表示电压信号最新更新的过零点位置，tv_i-1表示电压信号最新更新的相邻过零点位置，dt₁表示数据的采样时间间隔；

所述步骤S22包括以下子步骤：

S221：根据电流信号的次同步分量频率f_sub和电压信号的基波分量频率f₀，计算超同步分量频率f_sup并构建电压电流信号的特征值向量λ；

S222：根据特征值向量λ构建电压电流信号的范德蒙矩阵A；

S223：利用最小二乘法，根据范德蒙矩阵A得到电压电流信号中每个模态的相量，完成次同步振荡参数识别；

所述步骤S221中，超同步分量频率f_sup的计算公式为：

f_sup＝2*f₀-f_sub

其中，f_sub表示电流信号的次同步分量频率，f₀表示电压信号的基波分量频率；

电压电流信号的特征值向量λ的计算公式为：

其中，ω_i＝2πf_i，λ_n,n＝1,2,...,2i表示对应于所有q个模态的q对共轭的特征值，T表示矩阵转置运算，q表示模态个数，σ_i表示所分析信号中第i个分量的阻尼，f_i表示所分析信号中第i个分量的频率，ω_i表示对应于频率f_i的角频率，j表示虚部单位；

所述步骤S222中，范德蒙矩阵A的表达式为：

其中，λ^p _q表示特征值向量中第i个元素的p次幂，p＝0,1,…,m-1，q＝1,2,…,2q，m表示电压电流号中利用的点的数量；

所述步骤S223中，电流信号中每个模态的相量C_I的计算公式为：

C_I＝(A^TA)^-1A^TY_I

其中，Y_I＝[p_I1,p_I2,p_I3,...,p_Im]^T表示在最小二乘法拟合的过程中所使用的电流信号的点组成的相量，p_Ii,Ii＝I1,I2,...,Im表示在最小二乘法拟合的过程中所使用的所有电流数据，

电压信号中每个模态的相量C_V的计算公式为：

C_V＝(A^TA)^-1A^TY_V

其中，Y_V＝[p_V1,p_V2,p_V3,...,p_Vm]^T表示在最小二乘法拟合的过程中所使用的电压信号的点组成的相量，p_Vi,Vi＝V1,V2,...,Vm表示在最小二乘法拟合的过程中所使用的所有电压数据。

2.根据权利要求1所述的基于二阶段改进ITD算法的风电场次同步振荡监测方法，其特征在于，所述步骤S1中，对电压电流信号进行预处理的方法为：采用8阶低通FIR滤波器进行滤波处理，其截止频率是120Hz。

3.根据权利要求1所述的基于二阶段改进ITD算法的风电场次同步振荡监测方法，其特征在于，所述步骤S3中，根据次同步振荡决策条件判断风电场发生次同步振荡的方法为：若在60Hz系统中，振荡频率在5-55Hz范围内，则风电场发生次同步振荡，否则风电场未发生次同步振荡；

若基波幅值的10％大于预警阈值，则风电场发生次同步振荡，否则风电场未发生次同步振荡；

若最小检测时间大于设定时间阈值，则风电场发生次同步振荡，否则风电场未发生次同步振荡。

Images