CN108490257A - 一种基于频谱曲线拟合的短时窗间谐波测量方法 - Google Patents

Abstract

本发明公开了属于电力系统谐波测量技术领域的一种基于频谱曲线拟合的短时窗间谐波测量方法，包括：选择分析数据时长；使用汉宁窗抑制频谱泄露；根据测量精度需求，对加窗后的数据补零使得数据长度为2的整数次幂；对补零后的数据进行快速傅里叶变换，选取主瓣内的若干条谱线进行多项式拟合，求得多项式极值点对应的幅值、频率；由相频特性计算频率所对应的相角。本发明可以准确快速地进行间谐波测量，测量精度不仅满足标准要求，并高出标准一个数量级，在噪声较大的条件下依然保持较高的精度，计算量小，易于在微型嵌入式装置中实现。

Description

一种基于频谱曲线拟合的短时窗间谐波测量方法

技术领域

本发明属于电力系统谐波测量技术领域，尤其涉及一种基于频谱曲线拟合的短时窗间谐波测量方法。

背景技术

随着可再生能源的快速发展与电力电子装置的大量渗透，导致电网电气量中产生大量谐波、间谐波，引发电力系统的次同步振荡现象，严重时会导致火电机组轴系扭振保护动作、机组跳闸。近年来，新能源汇集地区次同步震荡现象频发，次同步震荡的监测与防治受到越来越多的研究机构地重视。相量测量单元(Phasor Measurement Unit，PMU)以其高精度及高速的上传频率，在电力系统动态安全监控中得到广泛应用。然而目前主流的商用PMU并不具备间谐波测量功能，只能依靠人为调取录波数据进行事后分析，无法对电力系统次同步震荡进行实时、准确地监测。

快速傅立叶变换(Fast Fourier transform，FFT)因为其算法简单，计算速度快，易于微机实现等优点，成为目前谐波测量中应用最广的方法。但是，FFT由于数据截断导致的频谱泄露和非整倍频采样造成的栅栏效应在短时窗的情况下会严重影响间谐波的测量精度。

发明内容

针对上述问题，本发明提出了一种基于频谱曲线拟合的短时窗间谐波测量方法，其特征在于，包括以下步骤：

步骤1：根据待测电压、电流信号的频率分布和动态变化情况选择分析数据时长T，在此基础上确定采样数据个数N；

步骤2：对N个采样数据加汉宁窗抑制频谱泄露；

步骤3：根据装置的计算性能和精度需求，对加窗后的采样数据补零；

步骤4：对补零后的采样数据进行快速傅里叶变换，得到离散的幅频特性与相频特性，选取幅频特性主瓣内若干条谱线进行多项式拟合，求得多项式极值点的幅值、频率；

步骤5：由相频特性计算步骤4得到的频率所对应的相角。

所述采样数据为电压与电流信号的采样值。

所述步骤1中，确定采样数据个数N的计算公式为：

N＝F_S*T+1

式中，Fs为采样率，N为采样点个数，T为分析数据时长。

所述步骤3对加窗后的采样数据补零，是在加窗后的数据末尾补若干个数值“0”，使补零后的数据长度N₀＝2^k，k为整数。

所述步骤4选取幅频特性主瓣内若干条谱线进行多项式拟合的计算方法为：

步骤1：选取幅频特性主瓣内的n条谱线，进行n-1次多项式拟合，得到多项式的各项系数[p₀ p₁ p₂ … p_n-1]^Τ，计算公式如下：

其中，[X₁ X₂ X₃ … X_n]^T为从左到右各条谱线的幅值，f₁为左侧第一条谱线对应的频率，Δf为补零后的频率分辨率，F_s为采样率，N₀为补零后的数据个数；

将式(1)简记为：

X＝BP (3)

其中，系数矩阵B为满秩方阵，可解得：

P＝B^-1X (4)

式中，X＝[X₁ X₂ X₃ … X_n]^T为各条谱线的幅值，P＝[p₀ p₁ p₂ … p_n-1]^Τ为多项式的各项系数；

步骤2：由多项式的各项系数[p₀ p₁ p₂ … p_n-1]^Τ求得多项式在区间(0,(n-1)Δf)内的极大值点(f_m,M)，得到间谐波测量值，包括幅值测量值M，频率测量值f₁+f_m。

所述步骤5中，由相频特性计算步骤4得到的频率所对应的相角，即间谐波相角，其计算公式为：

若则

若则

其中，f_l和f_r为相频特性中距离f最近的两个频率点，其对应的相角分别为和

本发明的有益效果在于：

(1)本发明可快速、准确地测量谐波/间谐波的频率与幅值，其测量精度不仅满足《新能源次同步振荡监测系统技术规范》的要求，并普遍高出标准一个数量级。

(2)本发明适用于微型PMU单元，适应性强，可根据测量的精度需求和装置的计算能力选择补零的数量，在噪声较大的条件下依然能够保持较高的精度，并且计算量不大，易于在微型嵌入式装置中实现。

附图说明

附图1为一种基于频谱曲线拟合的短时窗间谐波测量方法流程图；

附图2为矩形窗的幅频特性；

附图3为汉宁窗的幅频特性；

附图4为海明窗的幅频特性；

附图5为布莱克曼窗的幅频特性；

附图6为补零对栅栏效应的改善示意图；

附图7为对补零后的数据进行FFT变换后得到的离散谱线；

附图8为频谱曲线拟合取峰示意图；

附图9为新疆地区风电场主变高压侧电流频谱；

附图10为新疆地区风电场主变高压侧电流中35.7Hz频率分量的幅值变化情况；

附图11为对采样数据加窗后的信号；

附图12为采用补零法得到的信号；

具体实施方式

下面结合附图和实施例对本发明进行详细说明。

附图1为本发明提出的一种基于频谱曲线拟合的短时窗间谐波测量方法流程图。如图1所示，所述方法包括如下步骤：

步骤1：根据待测电压、电流信号的频率分布和动态变化情况选择分析数据时长T，在此基础上确定采样数据个数N；

步骤2：对N个采样数据加汉宁窗抑制频谱泄露；

步骤3：根据装置的计算性能和精度需求，对加窗后的采样数据补零；

步骤4：对补零后的采样数据进行快速傅里叶变换，得到离散的幅频特性与相频特性，选取幅频特性主瓣内若干条谱线进行多项式拟合，求得多项式极值点的幅值、频率；

步骤5：由相频特性计算步骤4得到的频率所对应的相角。

具体的，所述步骤1中根据待测电压、电流信号的频率分布和动态变化情况选择分析数据时长T，由于T的取值关系到动态响应时间和谱线的频率分辨率，如果选择太大则不利于观测信号的动态变化情况，如果选择太小则导致谱线过于稀疏，影响测量精度，因此需要根据采样数据的动态变化情况选择分析数据时长T，在确保信号变化不大的前提下，T的取值越大越好。在确定T的取值的基础上确定采样数据个数N，计算公式如下：

N＝F_S*T+1

式中，Fs为采样率，N为采样点个数，T为分析数据时长。

具体的，所述步骤2中对N个采样数据加汉宁窗抑制频谱泄露。频谱泄露会导致不同频率分量间发生旁瓣干涉，影响测量精度，故需对采样数据加窗抑制频谱泄露，窗函数的选取需考虑到频谱泄露程度以及频域主瓣的宽度，能量越集中的窗函数幅频特性的主瓣越宽，容易导致频率相近分量发生主瓣干涉，无法分辨。附图2-5分别为矩形窗、汉宁窗、海明窗和布莱克曼窗的幅频特性，对比图2-5可以看出，矩形窗频谱泄露严重，主瓣宽度较窄，海明窗与汉宁窗主瓣宽度均为矩形窗的两倍，海明窗的旁瓣小，衰减慢，汉宁窗的第一旁瓣较大，衰减快，布莱克曼窗基本没有泄露，主瓣宽度为矩形窗三倍，综合考虑，选择汉宁窗进行分析；

具体的，所述步骤3中，根据装置的计算性能和精度需求，对加窗后的采样数据补零。所述在加窗后的采样数据末尾补零是指在加窗后的数据末尾人为的补上一定数量的数值“0”，使得补零后的数据长度N₀＝2^k，k为整数。一方面可提高快速傅里叶变换的计算效率，另一方面，通过增加补零的数量，可减轻栅栏效应的影响，提高测量精度。附图6为补零对栅栏效应的改善示意图，由图6可以看出，通过增加补零的数量，可以在不改变主瓣宽度的情况下人为地提高频率分辨率，使得谱线更加密集，测量精度更高，但是，N₀越大意味着参与FFT的点数越多，运算量也会相应增加，因此需根据装置的计算性能和需要的测量精度进行选择。

具体的，所述步骤4中，对补“0”后的数据进行快速傅里叶变换(FFT)，得到如图7所示的离散谱线，选取频谱主瓣内的若干条谱线进行多项式拟合，求得多项式极值作为间谐波测量值，具体包括以下步骤：

步骤1：选取频谱主瓣内的n条谱线，进行n-1次多项式拟合，得到多项式的各项系数[p₀ p₁ p₃ … p_n-1]^Τ；

其中，[X₁ X₂ X₃ … X_n]^T为从左到右各条谱线的幅值，f₁为左侧第一条谱线对应的频率，Δf为补零后的频率分辨率，F_s为采样率，N₀为补零后的数据个数；

将式(1)简记为：

X＝BP (3)

系数矩阵B为满秩方阵，可解得：

P＝B^-1X (4)

式中，X＝[X₁ X₂ X₃ … X_n]^T为各条谱线的幅值，P＝[p₀ p₁ p₃ … p_n-1]^Τ为多项式的各项系数。

实际运行时可提前计算好矩阵B^-1，避免重复求逆运算，易于在微型嵌入式装置中实现。

步骤2：由多项式的各项系数[p₀ p₁ p₃ … p_n-1]^Τ求得多项式在区间(0,(n-1)Δf)内的极大值点(f_m,M)，得到幅值测量值M，频率测量值f＝f₁+f_m，如图8所示的频谱曲线拟合取峰示意图。

具体的，所述步骤5中，由相频特性求得间谐波相角的计算公式为：

其中，f_l和f_r为相频特性中距离f最近的两个频率点，其对应的相角为和若则若则

实施例1

为了对本发明进一步详细说明，下面以新疆地区风电场为例对电力系统进行谐波测量。

(1)根据采样数据的动态变化情况选择分析数据时长T，确定采样点个数N。如图9所示的新疆地区风电场主变高压侧电流频谱，其中含有35.7Hz和64.4Hz的间谐波，其中含35.7Hz频率分量的幅值变化情况如图10所示，由图10可以看出，幅值在1至2秒时间内变化约为10％，故时长T不宜超过1秒，本发明选取T＝1，假设装置的采样率Fs＝10000，通过计算得到采样值个数为N＝Fs*T+1＝10001。

(2)对上述计算得到的10001个采样点加窗抑制频谱泄露，对采样数据加汉宁窗，加窗后的信号如图11所示。

(3)根据装置的计算性能和精度需求，采用补零法提高频率分辨率；在加窗后的数据末尾补6383个“0”值，则补零后的采样点个数N₀＝16384＝2¹⁴，如图12所示。

(4)对补零后的数据进行快速傅里叶变换，得到离散频谱，选取频谱主瓣内的3条谱线，进行2次多项式的加权最小二乘拟合得到多项式的各项系数[p₀ p₁ p₂]，计算过程如下：

(5)由计算得到的多项式的各项系数[p₀ p₁ p₂]求得多项式在区间(0,1.22)内的极大值点(f_m,M)，得到幅值测量值M，频率测量值f₁+f_m，即间谐波测量值。

为了证实本方法的有效性，下面对上述测量方法进行仿真测试，包括基本测试、增加补零数量测试、增加谱线数量与多项式次数测试和噪声测试。

假设理想信号幅值为1，频率为34.567Hz，初始相角为0°，采样率Fs＝1200。

1、基本测试

测试参数：数据窗长T＝1,N＝1201，N₀＝2048，n＝3，测试结果如表1所示：

表1基本测试

测量量	理论值	测量值	测量误差
				频率(Hz)	34.567	34.5841	0.0172
幅值	1	1.000830	0.0830％
				相角(°)	0	-0.0057	-0.0057

2、增加补零数量测试

测试参数：数据窗长T＝1,N＝1201，n＝3，改变N₀的取值，测试结果如表2所示：

表2增加补零数量测试

	测量频率(Hz)	频率误差(Hz)	测量幅值	幅值误差(％)	测量相角(°)	相角误差(°)
							N0＝2048	34.5842	0.0172	1.000830	0.0830	-0.0057	-0.0057
N0＝4096	34.5755	0.0085	1.000831	0.0831	0.0028	0.0028
							N0＝9192	34.5708	0.0039	1.000825	0.0825	0.0014	0.0014

3、增加谱线数量与多项式次数测试

测试参数：数据窗长T＝1,N＝1201，N₀＝2048，改变n的取值，测试结果如表3所示：

表3增加谱线数量与多项式次数测试

4、噪声测试

测试参数：T＝1,N＝1201，含20dB噪声，测试结果如表4所示：

表4噪声测试

	频率(Hz)	频率误差(Hz)	幅值	幅值误差(％)	测量相角(°)	相角误差(°)
							N0＝2048，n＝3	34.5832	0.0162	0.999322	0.068	0.0440	0.0440
N0＝2048，n＝7，	34.5827	0.0157	0.999324	0.068	0.0424	0.0424
							N0＝4096，n＝3，	34.5749	0.0079	1.001704	0.170	0.6364	0.6364
N0＝4096，n＝7	34.5747	0.0077	1.001706	0.171	0.6154	0.6154

由表4可知，在含有20dB噪声情况下，测量精度依然能够满足要求。

测试结果表明，本发明提供的测量方法可以准确快速地进行间谐波测量，其测量精度不仅满足标准要求，并普遍高出标准一个数量级。并且该方法适应性强，可根据测量的精度需求和装置的计算能力选择补零的数量，在噪声较大的环境下依旧能够保持较高的测量精度。

此实施例仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (6)

1.一种基于频谱曲线拟合的短时窗间谐波测量方法，其特征在于，包括以下步骤：

步骤1：根据待测电压、电流信号的频率分布和动态变化情况选择分析数据时长T，在此基础上确定采样数据个数N；

步骤2：对N个采样数据加汉宁窗抑制频谱泄露；

步骤3：根据装置的计算性能和精度需求，对加窗后的采样数据补零；

步骤4：对补零后的采样数据进行快速傅里叶变换，得到离散的幅频特性与相频特性，选取幅频特性主瓣内若干条谱线进行多项式拟合，求得多项式极值点的幅值、频率；

步骤5：由相频特性计算步骤4得到的频率所对应的相角。

2.根据权利要求1所述的一种基于频谱曲线拟合的短时窗间谐波测量方法，其特征在于，所述采样数据为电压与电流信号的采样值。

3.根据权利要求1所述的一种基于频谱曲线拟合的短时窗间谐波测量方法，其特征在于，所述步骤1中，确定采样数据个数N的计算公式为：

N＝F_S*T+1

式中，Fs为采样率，N为采样点个数，T为分析数据时长。

4.根据权利要求1所述的一种基于频谱曲线拟合的短时窗间谐波测量方法，其特征在于，所述步骤3对加窗后的采样数据补零，是在加窗后的数据末尾补若干个数值“0”，使补零后的数据长度N₀＝2^k，k为整数。

5.根据权利要求1所述的一种基于频谱曲线拟合的短时窗间谐波测量方法，其特征在于，所述步骤4选取幅频特性主瓣内若干条谱线进行多项式拟合的计算方法为：

步骤1：选取幅频特性主瓣内的n条谱线，进行n-1次多项式拟合，得到多项式的各项系数[p₀ p₁ p₂ … p_n-1]^Τ，计算公式如下：

其中，[X₁ X₂ X₃ … X_n]^T为从左到右各条谱线的幅值，f₁为左侧第一条谱线对应的频率，Δf为补零后的频率分辨率，F_s为采样率，N₀为补零后的数据个数；

将式(1)简记为：

X＝BP (3)

其中，系数矩阵B为满秩方阵，可解得：

P＝B^-1X (4)

式中，X＝[X₁ X₂ X₃ … X_n]^T为各条谱线的幅值，P＝[p₀ p₁ p₂ … p_n-1]^Τ为多项式的各项系数；

步骤2：由多项式的各项系数[p₀ p₁ p₂ … p_n-1]^Τ求得多项式在区间(0,(n-1)Δf)内的极大值点(f_m,M)，得到间谐波测量值，包括幅值测量值M，频率测量值f₁+f_m。

6.根据权利要求1所述的一种基于频谱曲线拟合的短时窗间谐波测量方法，其特征在于，所述步骤5中，由相频特性计算步骤4得到的频率所对应的相角，即间谐波相角，其计算公式为：

若则

若则

其中，f_l和f_r为相频特性中距离f最近的两个频率点，其对应的相角分别为和