基于谱序列变换的高精度谐波参数估计算法

2020-10-31杨喜汪旭明3陈炳权张仁民王向明4雷可君

中南大学学报（自然科学版） 2020年9期

杨喜,2，汪旭明3，陈炳权，张仁民，王向明4，雷可君

(1.吉首大学信息科学与工程学院，湖南吉首，416000；2.东南大学信息科学与工程学院，江苏南京，210096；3.吉首大学物理与机电工程学院，湖南吉首，416000；4.上交所技术有限责任公司，上海，200120)

随着现代电力电子技术的快速发展，越来越多的非线性元器件被应用于电力系统[1]，在电网中产生大量谐波污染。谐波污染使得电能质量下降，从而严重影响电力系统的安全和经济运行。为了更好地进行谐波污染治理，必须对谐波参数进行高精度估计[2]。快速傅里叶变换(fast Fourier transform, FFT)是一种经典的谐波参数估计方法。然而，由于在实际运行过程中电网基波处于波动状态，因此，实际的采样过程通常处于非同步状态。直接采用FFT算法对谐波信号进行参数估计，其过程相当于对原始信号序列施加矩形窗。矩形窗实现方式虽然简单，但这种窗所引起的频谱泄露非常严重，其旁瓣峰值衰减仅为13 dB，由此导致FFT 算法对应的估计结果精度较低。为了减小频谱泄漏在谐波参数估计过程中带来的负面影响，采用在时域加窗处理和频域谱线插值处理相结合的方法来提高估计精度[3-13]。在前期的时域加窗处理过程中，常用的窗函数包括Hanning 窗、Blackman 窗、Nuttall 窗、Rife-Vincent窗和MSD窗。这些窗函数相对简单，但其旁瓣衰减速度有限而不利于抑制频谱泄露。为进一步加速窗函数旁瓣衰减，近年来，一些更加复杂的优化窗相继被提出，如乘法窗[11]、卷积窗[12-13]等。然而，当采用一些复杂的窗函数时，频谱泄露的抑制效果尽管有所提高，但另一方面也将不可避免地导致更高的计算复杂度而不利于信号参数的实时估计。通过加窗处理可以在一定程度上抑制长距离频谱泄漏的影响，但其对短距离范围内频谱泄露的抑制效果一般，因此，在随后的频域谱线插值处理过程中通过选择峰值谱线附近不同数量的谱线进行插值，以达到减小短距离频谱泄漏的影响，从而进一步提高参数估计精度。例如，文献[3，4，8]提出了加窗双谱线插值算法，文献[5，9]提出了加窗三谱线插值算法，文献[6，10]提出了加窗四谱线插值算法，文献[7]则提出了一种加窗六谱线插值算法。这些插值方法在一定程度上都可以提高谐波参数估计的精度。值得注意的是，上述基于不同数量谱线的插值方法有着各自的特点：对于利用较少谱线进行插值的方法，其插值公式推导和相关参数数值计算的复杂度较低，但其谐波参数估计的精度也相应较低；而对于利用较多谱线进行插值的方法，虽然参数估计的精度有所提高，但其计算复杂度也相应提高。相比较而言，文献[9]提出的基于Nuttall 窗的三峰插值谐波算法较好地兼顾了参数估计精度和计算复杂度。与此同时，考虑到FFT 方法在工程上简单易行，一些研究者提出了基于FFT的改进谐波参数分析方法，如文献[14-15]提出了基于FFT 序列进行变换处理的谐波参数估计算法，在此基础上，文献[16]通过在变换处理之前引入时域加窗处理环节以进一步减小频谱泄露，从而提高后继处理过程中谐波参数估计精度。整体而言，基于FFT 的改进算法虽然比前述基于加窗与插值处理相结合的算法其实现复杂度有所降低，但在实际应用过程中，谐波参数估计的精度仍有待进一步提高。本文在分析信号谱序列衰减特征的基础上，提出一种新的基于谱序列变换的高精度谐波参数估计算法。该算法通过对谱序列进行特定的加权变换处理，减少了频谱泄漏造成的谐波间干扰，极大地提高了谐波参数估计的精度。与此同时，由于该方法只需对信号FFT 变换后的谱序列直接进行简单的变换处理，因此，具有计算复杂度低的特点。

1 信号频谱分析

设单频率无限长余弦信号为

其中：A0，f0和φ0分别为信号的幅值、频率、相位。信号对应的抽样序列xd(n)为

其中：ω0= 2πf0/fs；fs为抽样频率。将xd(n)截断成1个长为N点的序列x(n)，这相当于将xd(n)乘上1个长为N点的窗函数wN(n)：

注意到对x(n)进行FFT 变换时相当于预先对x(n)加了1个矩形窗，此时，wN(n)可以表示为

原信号的复指数序列形式可以表示为

利用频域卷积定理，x(n)的离散时间傅里叶变换可表示为

式中：Xd(ejω)表示原序列xd(n)的离散时间傅里叶变换，

其中：δ(ω)为冲激函数。WN(ejω)表示wN(n)的离散时间傅里叶变换，

因此，原序列的频谱为

利用冲激函数的性质可以得到

式(11)为连续复谱函数，由于该复谱函数的正、负频率点具有对称性，且负频点对谐波参数估计的影响较小，所以，在分析过程中一般忽略负频率成分的影响，而只对信号的正频率部分加以分析和利用。在频域内取N点，以采样间隔Δf=fs/N对连续谱进行离散化。相应地，在式(11)中令ω= 2πn/N(n= 0,1,…,N- 1)，ω0= 2πk0/N即得到离散后的谱序列为

又由于N≫1，故谱序列可进一步近似为

2 基于谱序列变换的谐波分析

2.1 谱序列变换原理

由式(13)可知，当处于同步采样时仅当n=k0时，有X(n)=A0Nejφ0/2。此时，谱线的峰值恰好等于准确值，而在其他谱线上有X(n)= 0。电网中基波频率处于波动状态使得实际采样总是处于非同步状态，即k0为一个非整数，故设k0=k1+Δk，其中k1为k0的整数部分，Δk为其小数部分，故有

利用三角函数和差化积公式，式(14)可以表示为

由式(15)可知，当信号处于非同步采样时，频谱的峰值不等于准确值，而是在其他谱线上产生了分量，即产生了频谱泄漏，并且各次谐波之间将产生相互干扰，从而使得对谐波参数的准确估计十分困难，因此，如何降低频谱泄漏成为高精度谐波参数估计的关键。进一步令

则式(15)可以简洁地表示为

式(17)表明信号的幅度谱与||成反比。这意味着随着||增大，远离真实频率位置的谱线幅值只是以O()的速度衰减。由于这些谱线的幅值衰减速度不够快，频谱泄露将对谐波参数的估计造成很大干扰，极大地影响谐波参数估计的精度。为了减小频谱泄漏对谐波参数估计精度带来的负面影响，需要尽量加快偏离真实频率位置谱线幅值的衰减速度，为此，本文新设计出一个衰减速度为O(-11)的新序列：

需指出的是，尽管可以选择更多的谱线来构造衰减速度更快的新序列，但考虑到算法复杂度和估计精度的平衡，在式(18)所构造的新序列中只是选择了真实频率位置左右各5根谱线参与谐波参数的估计。注意到Y(n)可以看作1个有理多项式，因此，根据代数理论，新序列Y(n)可以进一步展开为如下表达式：

联立式(18)和式(19)可求得

由式(15)～(17)可得到

式中：i= 1,2,…,5。所以，Y(n)可以等价地表示为

可见新的谱序列Y(n)实际上是将真实频率附近的11 个点进行加权处理，由此频谱序列的衰减速度由原来的1/||变为衰减速度得到极大提高，使得离真实谱线较远处频点的幅值急剧减小，从而有效地抑制频谱泄漏。为验证本文算法对频谱泄漏的抑制效果，现对1个由基波和三次谐波组成的复合信号进行谱分析，信号为

其中：f0= 50 Hz；fs= 520.5Hz。对信号直接进行FFT处理和基于谱序列变换后得到的频谱特征分别如图1和图2所示(其中，频点为频域中离散谱线序号)。对比图1和图2可见：利用FFT算法对信号直接处理得到的谱序列频谱泄漏严重，而采用本文所提方法对谱序列进行变换后，新的谱序列中具有较大幅值的谱线主要集中在峰值谱线附近，减少了谐波之间的干扰，从而达到了进一步抑制频谱泄漏的目的，进而提高了谐波参数估计的精度。

图1 直接FFT频谱图Fig.1 Spectrum graph of direct FFT

2.2 谐波参数估计方法

非同步采样使得谐波真实频率的位置与峰值谱线的位置并不对应，而是位于相邻的2根谱线之间。设真实频率偏离峰值谱线的频率为Δk，若能求得Δk，则可准确计算真实谱线的位置，进而估计相关参数。为此，先在变换后的序列中搜寻幅值最大的2 根谱线Y(k1)和Y(k1+ 1)。由式(16)和式(18)得：

图2 基于谱序列变换后的频谱图Fig.2 Spectrum graph of spectral sequence transformation

设|Y(k1)|和|Y(k1+ 1)|的比值为ε0，则

可以求出偏移量为

结合式(16)和式(18)不难得到谐波幅值和相位的计算表达式：

其中：arg( )Y(k1)为Y(k1)的相位。经归纳，谐波参数估计的具体步骤如下。

1)对序列x(n)进行FFT 运算，得到谱序列X(n)；

2)由原谱序列X(n)根据式(22)计算得到新的谱序列Y(n)；

3)在变换后的序列Y(n)中搜索幅值最大的2根谱线Y(k1)和Y(k1+ 1)；

4)求出|Y(k1)|和|Y(k1+ 1)|的比值ε0；

5)通过ε0与Δk之间的关系式(27)计算偏移量Δk；

6)根据式(28)和(29)计算出幅值A0和相位φ0。

基于谱序列变换的谐波参数估计流程图如图3所示。

图3 谐波参数估计的流程图Fig.3 Flow chart of harmonic parameter estimation

上面以单频信号为例介绍了谐波参数估计的方法，该方法可以方便地推广到用于由基波和谐波组成的复合信号的参数估计。具体而言，在复合信号参数估计中，先计算出基波频率f1=(k1+Δk)·fs/N，再将第i次谐波的频域搜索范围设置为[i·f1- 5,i·f1+ 5]，并在该范围内搜索出幅值最大的2根谱线，根据式(24)～(27)即可依次计算出各次谐波的幅值和相位。

3 仿真与分析

为了验证本文所提算法的有效性，利用本文算法、文献[9]中算法和文献[16]中算法分别对单频信号和复合信号及相关的重要场景进行仿真实验对比。

3.1 单频信号场景

对单频信号作参数估计：

其中：A= 25V；f0= 35.5Hz；fs= 100 Hz；φ=1.5rad；采样点数N为64。

利用文献[9]中算法计算的幅值A和相位φ为：

利用文献[16]中算法计算的幅值A和相位φ为：

利用本文算法计算的幅值A和相位φ为：

由计算结果可知，利用文献[9]中算法估计幅值和相位的相对误差分别为4.80×10-9和9.34×10-8，利用文献[16]中算法估计幅值和相位的相对误差分别为8.43×10-9和4.33×10-8，而利用本文算法估计幅值和相位的相对误差分别为3.55×10-11和3.24×10-10。显然，对于单频信号而言，本文算法比文献[9]和文献[16]中所提算法具有更高的参数估计精度。

3.2 多频信号场景

设1个由基波和2～21次谐波组成的复合信号，其基波频率为f1= 50.2 Hz，表达式为

式中：A1为基波的幅值；Ai为各次谐波的幅值；φ1为基波的相位；φi为各次谐波的相位。在仿真过程中，具体参数设置如表1所示。

设置采样频率fs= 2 300 Hz，采样点数N=1 024，对信号分别采用本文算法、文献[9]中算法和文献[16]中算法进行参数估计，这3 种算法的幅值和相位估计相对误差对比结果如图4所示。从图4(a)可见：对于基波的幅值估计相对误差，本文算法所得结果比文献[9]中算法所得结果高近7个数量级，比文献[16]中算法结果高3 个以上数量级；对于基波以外的谐波，本文算法得到的估计相对误差范围为7.39×10-11～3.03×10-14，文献[9]中算法结果相对误差范围为8.74×10-8～4.82×10-10，文献[16]中算法结果相对误差范围为7.18×10-9～6.10×10-11。从图4(b)可见：对于基波的相位估计相对误差，本文算法所得结果比文献[9]中算法结果高将近4个数量级，比文献[16]中算法结果高近3 个数量级；进一步，对于基波以外的各次谐波，本文所提算法的相对误差范围为1.28×10-9～8.14×10-13，而文献[9]中算法结果的相对误差范围为2.29×10-6～2.08×10-9，文献[16]中算法结果的相对误差范围则为6.21×10-7～6.34×10-10。总体来看，对于基波和谐波组成的复合信号，本文算法在幅值和相位参数的估计结果精度与文献[9]和文献[16]中算法结果精度相比显著提高。

表1 基波和各次谐波的幅值和相位Table1 Amplitude and phase of the fundamental and harmonics

图4 不同算法下谐波幅值和相位估计相对误差对比Fig.4 Comparison of harmonic amplitude and phase estimation relative error under different algorithms

在采样点数N为512，1 024 和2 048 时，幅值、相位估计相对误差仿真结果见图5。从图5可以看出：当采样点数由512 增加到1 024 时，幅值估计相对误差至少提高了2个数量级，对4次和11次谐波则提高了近5个数量级；相位估计相对误差至少提高了3个数量级，对7次和14次谐波则提高了近4 个数量级；当采样点数由512 增加到2 048时，幅值估计相对误差至少提高了3个数量级，对2次和15次谐波则提高了近6个数量级；同时，相位估计相对误差至少提高了3个数量级，对5次和9 次谐波则提高了近7 个数量级。显然，采样点数的增加可以有效地提高算法的估计精度。

图5 采样点数N不同时谐波幅值和相位估计相对误差对比Fig.5 Comparison of harmonic amplitude and phase estimation relative error when sampling points are different

表2 不同基波频率下幅值估计相对误差Table2 Amplitude estimation relative error at different fundamental frequencies

表3 不同基波频率相位估计相对误差Table3 Phase estimation relative error at different fundamental frequencies

3.3 基波频率变化场景

电网实际运行过程中的基波频率并不是固定不变的，通常在50 Hz左右波动，因此，考察算法在基波频率变化情况下参数估计性能的稳定性非常必要。GB/T 15945—2008 规定电力系统基波频率偏差最大范围为±0.5 Hz，故在仿真过程中基波频率设定在49.5～50.5 Hz[17-18]。信号时域表达式如式(31)所示，仿真过程中设置采样频率fs=2 300 Hz，采样点数N=1 024。表2和表3所示仿真结果表明：在基波频率变化时，本文所提算法仍然表现出优良的估计性能和良好的稳定性。

3.4 间谐波场景分析

电网在实际运行过程中可能会产生间谐波，其存在会造成电压的闪变并对电网产生危害，因此，在处理谐波过程中，对间谐波参数的高精度估计非常必要[19-20]。为了验证本文算法对间谐波参数估计的有效性，在仿真过程中采用1个由7次间谐波组成的复合信号(其基波频率f1= 50.1Hz)：

具体参数设置如表4所示。另外，在仿真过程中设置采样频率fs=2 300 Hz，采样点数N=1 024。采用文献[9]中算法、文献[16]中算法和本文算法分别对上述间谐波信号的幅值和相位参数进行估计，得到各次间谐波的幅值、相位估计的相对误差如表5和表6所示。从仿真结果可以看出：文献[9]和文献[16]中算法在一定程度上提高了间谐波参数估计的精度，但本文算法与文献[9]和[16]中算法相比，估计精度显著提高。