一种风电次同步振荡的关键参数辨识方法与流程

[0001]
本发明属于电力系统技术领域，更为具体地讲，涉及一种风电次同步振荡的关键参数辨识方法。


背景技术：

[0002]
近年来，随着以风力发电为代表的可再生能源的日益普及，新型次同步振荡(sub-synchronous oscillation，sso)频发，与汽轮发电机组轴系扭振所引起的传统sso仅限于电力系统的局部区域不同，可再生能源电力系统中发生的sso涉及多种组件，并且会在更大范围内传播。在已报告的事故中，发现sso具有极强的破坏性，从而极大地威胁了设备的安全性和系统的稳定运行。因此，sso的监测和辨识对于次同步振荡的抑制至关重要。
[0003]
现有辨识sso参数的方法可分为时域法、频域法和时频域法。其中，时域法主要利用瞬时电压或电流测量数据进行振荡辨识，如prony、旋转不变技术估计信号参数和粒子群优化，但这些方法需要模型的先验知识，预设参数对结果影响大。时频域方法例如，ibrahim时域模态辨识法和变分模态分解可以提取每个模态的参数，但这些方法会产生较大计算负担，且边缘趋势效应严重，易造成分解不准确。频域法主要基于傅立叶变换(discrete fourier transform,dft)实现，并通过改善dft性能来提高振荡参数辨识精度。为降低由栅栏效应和频谱泄露引起的误差，现有文献采用基于加窗函数的插值方法，常用的窗函数由hamming窗、高斯窗和rife-vincent窗等，其中rife-vincent窗与其他经典窗函数相比有良好的主瓣和旁瓣的特性，可获得较好的频谱泄露补偿效果。但是，这些方法所使用的瞬时数据仅存储在故障记录器(digital fault recorder，dfr)中，不易收集且数据不同步，主要用于事后的离线分析，无法实现sso的在线监测。
[0004]
同步相量测量数据更适用于在线监测，用于获取这些数据的相量测量单元(phasor measurement unit，pmu)已广泛应用于电力系统的动态监控中，现已记录多起sso事故。pmu使用固定的频率和窗口长度进行采样，其相量测量方法主要包括dft及其改进方法，如通过从pmu数据中检索sso信息的方法，但仅获得sso频率，而其他参数如幅值、阻尼系数和相位等辨识困难。近年来研究者提出了一些新技术来解决此问题，如提出通过pmu数据监测间谐波的方法，或通过校正同步相量的频谱泄漏来辨识sso的频率和幅度，但这两种方法没有考虑振荡信号的阻尼系数，会使sso幅值辨识结果不准确。最新研究提出了一种使用同步相量中的频谱泄漏信息来辨识关键sso参数的方法，但当基频和sso频率接近时，参数辨识的误差很大，且该方法在噪声的条件下鲁棒性较低。


技术实现要素：

[0005]
本发明的目的在于克服现有技术的不足，提供一种风电次同步振荡的关键参数辨识方法，通过添加阻尼rv-m窗函数的dft插值处理，完成风电次同步振荡的频率、阻尼系数、幅值和相位的辨识。
[0006]
为实现上述发明目的，本发明一种风电次同步振荡的关键参数辨识方法，其特征
在于，包括以下步骤：
[0007]
(1)、在风电系统的输出端提取包含基频分量和风电次同步振荡分量的瞬时电流或电压信号x(t)；
[0008][0009]
其中，a1,f1,φ1和a
s
,f
s
,φ
s
分别表示基频和风电次同步振荡分量的振幅、频率和相位，α
s
是风电次同步振荡的阻尼系数，下标s表示风电次同步振荡；
[0010]
(2)、相量测量单元以固定的采样频率fp对信号x(t)进行采样，并以上报频率f
r
向主站传输，主站得到采样后的离散信号x(n)；
[0011][0012]
其中，n为采样个数，f
p1
＝f1/f
p
，f
ps
＝f
s
/f
p
，f
p1
、f
ps
分别是基频和风电次同步振荡分量的归一化频率，α
ps
＝α
s
/f
p
是阻尼系数的归一化形式，f
pr
＝f
p
/f
r
是对同步相量进行重新采样的间隔；
[0013]
(3)、对离散信号x(n)加余弦组合窗rv-m，再进行傅立叶变换dft，其中，x(n)加窗后的dft变换表示为：
[0014][0015]
其中，w(n)为余弦组合窗，k＝0,1,
…
,n
p-1，n
p
为窗长；
[0016]
(4)、计算dft变换后第m个滑动时间窗的频谱；
[0017][0018]
其中，m＝1,2,
…
,m，m为窗函数的阶数；f
p1
＝f1/f
p
＝l1/n
p
，l1为归一化基频在频谱中的位置，f
ps
＝f
s
/f
p
＝l
s
/n
p
，l
s
为归一化次同步振荡频率在频谱中的位置；
[0019]
(5)、提取频谱x
w
(m,k)的基频x1(m,k)和风电次同步振荡分量频谱x
s
(m,k)；
[0020][0021][0022]
其中，(
·
)
*
是复数共轭符号，w1(
·
)代表rv-m窗函数的离散傅里叶变换形式；w
s
(
·
)代表阻尼rv-m窗函数的离散傅里叶变换形式，
[0023]
(6)、辨识风电次同步振荡的关键参数
[0024]
(6.1)、将频谱x
s
(m,k)简化为：
[0025][0026]
其中，ξ
x1
(k)表示来自基频分量的频谱泄露值；
[0027]
(6.2)、对简化后的x
s
(m,k)的进行三点差值；
[0028]
(6.2.1)、计算频谱x
s
(m,k)中的比率r1和r2；
[0029][0030]
其中，k2是频谱最高点位置，v(k2)是对应位置的幅值，k3和k1是频谱第二高和第三高点位置，v(k3)与v(k1)是对应第二和第三高位置的幅值；
[0031]
(6.2.2)、计算振荡频率f
s
；
[0032][0033][0034]
其中，d为频率校正因子，re(
·
)为复数的实部；
[0035]
(6.2.3)、计算阻尼系数α
s
；
[0036]
判断频率校正因子d是否等于0.5，如果d≠0.5，则按照如下公式计算阻尼系数α
s
；
[0037][0038]
如果d＝0.5，则此时辨识不再准确，需要在信号x(t)的末尾补一个零，然后再返回步骤(2)；
[0039]
(6.2.4)、计算幅值a
s
；
[0040][0041]
其中，x
s
(m,k2)是风电次同步振荡的频谱最大幅值，w
s
(k2)为频谱最高位置对应的阻尼rv-m窗函数的dft幅值；
[0042]
(6.2.5)、计算相位φ
s
；
[0043][0044]
至此风电次同步振荡的关键参数f
s
、α
s
、a
s
和φ
s
辨识完成。
[0045]
本发明的发明目的是这样实现的：
[0046]
本发明一种风电次同步振荡的关键参数辨识方法，首先在待辨识的风电系统的输出端提取包含基频分量和风电次同步振荡分量的瞬时电流或电压信号，然后相量测量单元以固定的采样频率对提取信号进行采样，并以上报频率向主站传输，主站得到采样后的离散信号，最后采用添加阻尼rv-m窗函数和dft三点插值的方法对风电次同步振荡的关键参数进行辨识。
[0047]
同时，本发明一种风电次同步振荡的关键参数辨识方法还具有以下有益效果：
[0048]
(1)、本发明利用系统基频分量和sso分量在频域中呈线性叠加的特征对sso关键参数辨识，无需振荡或电力系统模型的先验信息，因此能够适用电力系统的复杂运行条件。
[0049]
(2)、针对传统dft算法的两个主要系统误差源，即频谱泄漏和栅栏效应，本发明采用加阻尼rv-m窗函数和dft三点插值的方法对sso参数进行辨识，仿真结果表明本方法能有效抑制误差使辨识准确。
[0050]
(3)、本发明在非标称、噪声和动态等多种条件下的仿真结果表明，本发明参数辨识误差小于1％，其中频率辨识误差最小，小于10-4
％，具有高辨识精度；另外动态条件下低阶窗函数稳定性优于高阶。
[0051]
(4)、本发明方法是对同步相量数据使用加rv-m窗的dft插值法，计算复杂度低，能与pmu设备相结合，适用于sso参数的在线辨识。
附图说明
[0052]
图1是本发明一种风电次同步振荡的关键参数辨识方法流程图；
[0053]
图2是发生sso时的同步相量频谱；
[0054]
图3是示例信号示意图；
[0055]
图4是f
s
在10-45hz变化时的辨识误差示意图；
[0056]
图5是α
s
在0.01-1变化时的辨识误差示意图；
[0057]
图6是a
s
在10-50变化时的辨识误差示意图；
[0058]
图7是噪声条件下1阶rv-m窗的辨识误差示意图；
[0059]
图8是噪声条件下4阶rv-m窗的辨识误差示意图；
[0060]
图9是1阶rv-m窗、基频斜升、相位阶跃和幅值阶跃变化测试辨识误差示意图；
[0061]
图10是4阶rv-m窗、基频斜升、相位阶跃、幅值阶跃变化测试辨识误差示意图；
[0062]
图11是1-4阶rv-m窗φ
s
在-π到π变化时的相位辨识误差示意图。
具体实施方式
[0063]
下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。
[0064]
实施例
[0065]
图1是本发明一种风电次同步振荡的关键参数辨识方法流程图。
[0066]
在本实施例中，如图1所示，本发明一种风电次同步振荡的关键参数辨识方法，包括以下步骤：
[0067]
s1、在待处理的风电系统的输出端提取包含基频分量和风电次同步振荡分量的瞬时电流或电压信号x(t)，在大多数sso发生的情况下只存在一组sso分量，不失一般性，可将电流或电压信号建模为：
[0068][0069]
其中，a1,f1,φ1和a
s
,f
s
,φ
s
分别表示基频和风电次同步振荡分量的振幅、频率和相位，α
s
是风电次同步振荡的阻尼系数，下标s表示风电次同步振荡；
[0070]
s2、相量测量单元以固定的采样频率f
p
对信号x(t)进行采样，并以上报频率f
r
向主站传输，主站得到采样后的离散信号x(n)；
[0071][0072]
其中，n为采样个数，f
p1
＝f1/f
p
，f
ps
＝f
s
/f
p
，f
p1
、f
ps
分别是基频和风电次同步振荡分量的归一化频率，α
ps
＝α
s
/f
p
是阻尼系数的归一化形式，f
pr
＝f
p
/f
r
是对同步相量进行重新采样的间隔；
[0073]
s3、传统的dft方法存在的系统误差分析，因此本发明通过添加rv-m窗和dft变换两种措施来减少频谱泄露和栅栏效应引起的误差，也就是先对离散信号x(n)加余弦组合窗rv-m，再进行dft变换，其中，x(n)加窗后的dft变换表示为：
[0074][0075]
其中，w(n)为余弦组合窗，k＝0,1,
…
,n
p-1，n
p
为窗长；
[0076]
在本实施例中，rv-m窗是一种典型的余弦组合窗，w(n)的离散形式为：
[0077][0078]
其中，m为窗函数的阶数；
[0079]
s4、计算dft变换后第m个滑动时间窗的频谱；
[0080][0081]
其中，m＝1,2,
…
,m，m为窗函数的阶数；f
p1
＝f1/f
p
＝l1/n
p
，l1为归一化基频在频谱中的位置，f
ps
＝f
s
/f
p
＝l
s
/n
p
，l
s
为归一化次同步振荡频率在频谱中的位置；
[0082]
s5、提取频谱x
w
(m,k)的基频x1(m,k)和风电次同步振荡分量频谱x
s
(m,k)；
[0083][0084][0085]
其中，(
·
)
*
是复数共轭符号，w1(
·
)代表rv-m窗函数的离散傅里叶变换形式；w
s
(
·
)代表阻尼rv-m窗函数的离散傅里叶变换形式，
[0086]
s6、辨识风电次同步振荡的关键参数
[0087]
s6.1、将频谱x
s
(m,k)简化为：
[0088]
[0089]
其中，ξ
x1
(k)表示来自基频分量的频谱泄露值；
[0090]
s6.2、对简化后的x
s
(m,k)的进行三点差值；
[0091]
s6.2.1、在本实施例中，传统插值是采用实数波形数据不同，而本发明进行改进，以处理复数同步相量数据为基础，基于加rv-m窗插值dft的方法，通过计算频谱x
s
(m,k)中的比率r1和r2来精确辨识f
s
、α
s
和a
s
；其中，
[0092][0093]
其中，k2是频谱最高点位置，v(k2)是对应位置的幅值，k3和k1是频谱第二高和第三高点位置，v(k3)与v(k1)是对应第二和第三高位置的幅值，关系如图2所示；
[0094]
s6.2.2、计算振荡频率f
s
；
[0095][0096][0097]
其中，d为频率校正因子，re(
·
)为复数的实部；
[0098]
s6.2.3、计算阻尼系数α
s
；
[0099]
本发明讨论电力系统中需要被关注的sso，即阻尼系数大于0的情况，而频率校正因子的大小决定阻尼系数的准确度，因此我们先判断频率校正因子d是否等于0.5，如果d≠0.5，则按照如下公式计算阻尼系数α
s
；
[0100][0101]
如果d＝0.5，则此时辨识不再准确，需要在信号x(t)的末尾补一个零，通过补零后改变了原信号x(t)，然后再返回步骤s2，重复上述过程；
[0102]
s6.2.4、计算幅值a
s
；
[0103][0104]
其中，x
s
(m,k2)是风电次同步振荡的频谱最大幅值，w
s
(k2)为频谱最高位置对应的阻尼rv-m窗函数的dft幅值；
[0105]
s6.2.5、计算相位φ
s
；
[0106][0107]
至此风电次同步振荡的关键参数f
s
、α
s
、a
s
和φ
s
辨识完成。
[0108]
实验验证
[0109]
为了便于本领域人员理解发明的技术内容以及技术性能，下面结合附图对本发明内容和性能进一步阐释。
[0110]
建立含有sso的合成仿真信号，利用传输至主站的同步相量数据进行sso参数辨识。目前，50hz系统的传输速率有两种类型：f
r
＝f0和f
r
＝2f0。当f
r
＝f0时，由获得的同步相量不能辨识f
s
＞f0/2的sso。因此，本发明将pmu的传输速率设置为100hz，对发生sso时引发
的超同步谐波不在本发明的研究范围内。将pmu的采样频率f
p
设定为10khz，根据测试，窗长n
p
设置为200可弱化频谱泄露的影响以得到相对精确的辨识结果，这意味着相量数据的间隔频率为100hz，频率分辨率为0.5hz，每次处理200个同步相量，窗长对辨识结果的影响也在后文中进行了分析。在风电系统的输出端提取包含基频分量和风电次同步振荡分量的瞬时电流或电压信号x(t)，其中，f1＝50hz、f
s
＝44.6hz、a1＝100、a
s
＝40、α
s
＝0.15、φ1＝π/9及φ
s
＝-π/18，如图3所示，其中，(a)是瞬时波形，(b)是同步相量幅值波形，(c)是加矩形窗的同步相量频谱，(d)是加阻尼rv-m窗的同步相量频谱，通过对比可见阻尼rv-m窗的加入可更好抑制频谱泄露以精确辨识f1和f
s
接近的情况。
[0111]
引入辨识偏差ies(identification errors)以评估所提出方法的准确性，其计算公式为：
[0112][0113]
其中，p和分别为参数的真实值和估计值。为了验证该方法在不同条件下的性能，下面分析振荡参数对辨识结果的影响，由计算公式发现相位参数不影响幅值、频率、衰减系数参数的辨识，因此单独讨论φ
s
和rv-m窗阶数变化对相位辨识的影响。
[0114]
1)、f
s
变化对辨识精度的影响；
[0115]
建立合成仿真信号参数设定如下：a1＝100、a
s
＝10、φ1＝π/9、φ
s
＝π/35、α
s
＝0.15。通常规定电力系统的频率波动幅度不超过基频的2％，即50
±
1hz，因此设置基频为标称和非标称(49.1、49.5、50、50.3、50.5和50.9hz)，并设置振荡频率f
s
在10-45hz之间，步长为0.3hz；rv-m窗函数设置为1-4阶分别进行测试。
[0116]
各个参数辨识结果如图4所示，其中，(a)是频率辨识结果；(b)是阻尼系数辨识结果；(c)是幅值辨识结果，m＝1，2，3，4；当f1＝49.5、50和50.5hz时比其他频率情况下辨识误差低；当f
s
是0.5hz整数倍时，辨识误差为极小值，并在25hz时达到最小。这表明在这频率分辨率整数倍位置上时辨识结果更加准确；当f
s
接近f1时，辨识误差因频谱混叠随之增大，显然随着rv-m窗函数阶数的增加，辨识精度逐渐增大，受matlab计算精度限制，超过10-16
％结果为0，因此在3和4阶中频率和幅值某些点辨识结果显示为空白。总之，采用本发明提出的方法，三种参数的辨识误差均低于1％，且频率辨识误差最低，可以达到10-10％
，精度高于ieee标准中最大总矢量误差(max total vector error，max tve)为1％的要求。
[0117]
2)、α
s
变化对辨识精度的影响；
[0118]
建立合成仿真信号中参数设定如下：a1＝100、a
s
＝10、φ1＝π/9、φ
s
＝π/35。为防止采样频率与次同步频率相干，f
s
设置为27.25hz。在不同阶rv-m窗，标称频率和非标称频率的f1情况下，α
s
在0.01-1范围内变化，步长为0.01。各个参数辨识结果如图5所示，其中，(a)是频率辨识结果；(b)是阻尼系数辨识结果；(c)是幅值辨识结果，m＝1，2，3，4；从图5可以看出，α
s
的变化对频率和幅值的辨识影响较小，但会影响阻尼系数的估算。这是由于在α
s
变大时，sso频谱对于阻尼系数辨识变得更明显。与第一种情况相同，由于频谱泄漏，在非标称频率的情况下，这三个参数的辨识偏差较大，而增加rv-m窗阶数会减小辨识误差。
[0119]
3)、a
s
变化对辨识精度的影响
[0120]
合成仿真信号中参数设定如下:a1＝100、α
s
＝0.25、f
s
＝27.25、φ1＝π/9、φ
s
＝π/35。a
s
变化范围为10-50，步长为5，各个参数辨识结果如图6所示，其中，(a)是频率辨识结
果；(b)是阻尼系数辨识结果；(c)是幅值辨识结果，m＝1，2，3，4。三种参数辨识结果有着相同的趋势，当基频为49.5、50和50.5hz时，幅值的变化对辨识结果影响较小。而基频为其他频率时，辨识误差会随着a
s
的增大而减小，这是由于a
s
增大使x
s
在频谱中更为明显，从而使辨识更为容易，结果则更加精确。总之，辨识误差在1％以内，精度高于标准要求。
[0121]
4)、噪声对辨识结果的影响
[0122]
pmu数据信噪比(snr)通常为45db左右，snr定义为：
[0123][0124]
其中，σ2是高斯白噪声的方差。为了验证所提出算法的抗噪能力，分别对1)、2)和3)三种仿真案例中加入30db和40db高斯白噪声，验证本发明提出的算法，为便于观察和篇幅简洁，将m设置为1和4阶。辨识结果如图7和8所示，分别为m＝1和m＝4时噪声在30db和40db下f
s
、α
s
和a
s
对辨识结果的影响，合成模拟信号参数与前文设置相同，其中，(a)
–
(c)snr＝30db，(d)
–
(f)snr＝40db。辨识结果总结如下：1)在噪声情况下，f
s
比α
s
、a
s
辨识更加准确。当m＝1、f
s
在10-45hz时，三种参数辨识误差分别约为10-4
％、1％和10-2
％；2)由于引入的rv-m窗具有良好的抑制频谱泄漏性能，即使f
s
与f1非常接近，参数辨识误差波动趋势较小；3)随着α
s
和a
s
增加或rv-m窗阶数提高，辨识误差呈下降趋势，原因可解释为：当α
s
或a
s
增大时，sso在频谱中更加明显从而便于辨识；4)通常sso事件中α
s
＞0.1且a
s
＞0.15a1，因此所提方法能适用于大多数情况。当α
s
＜0.1或a
s
＜15(a
s
＜0.15a1)时，辨识结果误差较大，原因是此时sso分量较小，参数辨识更易受噪声影响，可通过提升rv-m窗阶数来减小误差，如图8采用4阶rv-m窗，a
s
在10-20变化时α
s
的辨识误差由10％降低至1％左右；α
s
在0.01-0.1变化时，幅值辨识误差由1％降低至0.1％左右。
[0125]
5)、动态条件下方法鲁棒性验证
[0126]
根据ieee标准，需要在动态条件(频率斜升变化，幅度阶跃变化和相位阶跃变化)对辨识算法进行验证。在动态条件下，信号模型可以改写为如下形式：
[0127][0128]
其中，u1(t)是阶跃函数，k
a
和k
x
分别是相位和幅值的调制系数。k
a
从-0.1到0.1变化，步长为0.02；k
x
从-π/18到π/18变化，步长为π/90。当k
a
变化时，k
x
设为0，反之亦然，两种情况下f1＝49.1、49.5、50、50.3、50.5和50.9hz，f
s
＝27.25hz。同样设置rv窗阶数为1阶和4阶。
[0129]
三个参数的辨识误差如图9和10所示，其中，(a)是频率辨识误差，(b)是阻尼系辨识误差数，(c是)幅值辨识误差；在频率斜升变化下，基频f1在48-52hz范围内以0.5hz/s的步长线性变化。为验证不同sso频率下的动态辨识结果，测试了10-45hz区间内五个sso频率。在该条件且窗函数不同阶数下，当基频动态变化时，三种参数辨识误差均在1％以内。随着次同步频率接近基本频率，辨识精度相应提高。此外，当f1＝49.5、50、50.5hz，即位于频率分辨率的整数位置时，辨识结果在所有情况下都最为准确。在幅度阶跃变化和相位阶跃变化下，两个动态条件对低阶窗函数的辨识结果影响较小，而高阶窗函数稳定性较差。
[0130]
根据ieee c37.118.1标准的定义，动态响应时间是将阶跃变化应用于输入之前和之后两次稳态测量之间转换的时间。为此，在动态测试中的频率、幅值和相位阶跃变化下进行了多组测试，三种动态条件下的平均动态响应时间分别为160ms、156ms和192ms。实际振
荡会持续一段时间(通常为几秒钟)，不会立即损坏设备，上百毫秒的动态响应时间和辨识误差能够满足实际应用的要求。
[0131]
6)、相位变化对参数辨识的影响
[0132]
为验证相位辨识精度，设置合成仿真信号的参数如下：f
s
＝27.25hz、f1＝49.1、49.5、50、50.3、50.5和50.9hz，a1＝100、a
s
＝100、α
s
＝0.25、φ1＝0，φ
s
变化范围是-π到π，步长为π/20。如图11所示，其中(a)～(d)分别对应1-4阶rv-m窗φ
s
在-π到π变化时的相位辨识误差结果；受matlab计算精度限制，低于10-16
％误差显示为0，相位辨识误差随着φ
s
的变化出现波动现象，随着rv-m窗阶数的提高而降低，且rv-m窗阶数越高，波动越大；另外在f1＝49.5、50和50.5hz时辨识误差相近。可见φ
s
变化会影响辨识精度，所提方法对相位的辨识误差在10-3
％以内，rv-m窗阶数的提高会提升辨识精度。
[0133]
尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。