基于ZoomFFT 和多谱线插值算法的谐波分析方法*

2024-01-15陈乐柱史琛筠

电气工程学报 2023年4期

朱坤陈乐柱许胜史琛筠

(1.安徽工业大学电气与信息工程学院马鞍山 243000；2.泰州市电能变换与控制工程技术研究中心泰州 225300)

1 引言

随着电力电子行业的发展和非线性负荷设备的增加，电网谐波污染进一步恶化，从而影响电力系统安全平稳运行。为了能够有效治理电网中存在的电能质量问题，首先需要对电网的电力参数进行检测，为保证对电力系统的有效滤波或补偿，需要保证检测参数的精准性，同时为保证对电力系统的及时治理，还需要保证参数的实时性。

目前，基于快速傅里叶变换(Fast Fourier transform，FFT)的算法因具备计算速度快等优点被广泛运用，但仍存在频谱泄漏、栅栏效应和频谱混叠等问题[1-4]。

理想状态下时域下的信号序列通常认为是无限长，而实际上在分析时域信号时都是截取某一段信号时序，需要加窗函数，相当于频域卷积，并且通常情况下在进行时域截取信号时，为非整数信号周期，在分析的频谱图上就会存在旁瓣。为了减少频谱泄漏，需要乘以加权的窗函数，例如Hanning 窗、Hamming 窗、Kaiser 窗、Blackman 窗、Nuttall 窗等。实际在对信号进行处理时，时域和频域都是离散的，连续的信号不会造成频率分量丢失，但离散信号时序经过FFT 从时域到频域会造成频谱中频率分量丢失，导致频谱不完整。

结合工程应用，采集的电信号中会存在与整数次谐波相近频率的间谐波含量，采用传统的快速傅里叶变换很难将其区分，此时需要提高频率分辨率，可以采用提高采样点数N或者减小采样频率的方法。其中也有学者提出频谱细化算法和加窗插值算法对频谱进行分析，例如线性调频 Z 变换(或CZT)[5]、小波变换[6]、复调制频谱细化(或ZFFT)[7-8]、Kaiser 窗双谱线插值算法[9]、Kaiser 窗三谱线插值算法[10]和Nuttall 窗三谱线插值[11]等。这些算法都能够提高频率分辨率，但仍有其局限性，例如CZT 算法不适用于多频电信号且相对密集的有效分析，且实时性相对较差；Kaiser 窗和Nuttall 窗的多谱线插值算法测量精度不能满足所有条件。

针对电网信号测量分析中需要数据精度的同时也需要满足其实时性的要求，本文结合复调制频谱细化(ZFFT)算法和基于Blackman-Harris 窗四谱线插值算法对测量信号进行分析，并结合主瓣干扰判定条件对算法进行处理，最后通过仿真验证改进型方法的精准度和实时性。

2 基于Blackman-Harris 窗四谱线插值校正算法

针对基于Blackman-Harris 窗双谱线插值算法[12]和基于Blackman-Harris 窗三谱线插值校正算法[13]存在的误差，本文提出基于Blackman-Harris 窗四谱线插值算法来减少校正误差，提高测量精确性。

2.1 Blackman-Harris 窗

加窗函数能够在发生频谱泄漏时有效改善并降低频谱泄漏带来的影响，同时提高频谱分辨率。通常窗函数具有如下几种特征参数：主瓣宽度、旁瓣衰减速率、最高旁瓣高度、幅值失真度等。因此，在选择窗函数时要结合实际测量电信号和分析窗函数的旁瓣特性。

表1 所示为部分窗函数的特性[14]，根据窗函数的特性和电信号的性质来选择所需的窗函数。

表1 部分窗函数特性

由表1 可知，Blackman-Harris 窗函数的旁瓣峰值电平相比较于其他窗函数较低(-92 dB)，旁瓣衰减速率为6 dB/oct。衰减速率越大，旁瓣峰值电平越小，对频谱泄漏的治理能力越好。综合考虑实际情况，本文选取Blackman-Harris 窗函数。

Blackman-Harris 窗的本质是一个4 项系数的加权余弦窗，其函数的时域表达式为

式中，K为余弦窗的项数；N为窗函数的数据长度，n=1,2,3,…,N-1。取K=3，得到一般表达式为

式中，a1=0.358 75，a2=0.488 29，a3=0.141 28，a4=0.011 68。

2.2 四谱线插值算法

为了方便理解和计算，对电信号x(t)采取只含单一次谐波分量的信号进行分析(式(3))，其中以采样频率Fs对信号x(t)进行采样得到的离散时间信号如式(4)所示

式中，f0表示单一次谐波的频率；A表示单一次谐波的幅值；φ0表示单一次谐波的相位[15]。

现对信号x(n)进行加窗函数w(n)处理，即得到式(5)

代入式(4)，对式(5)进行离散傅里叶变换，得到表达式为

忽略负频点的旁瓣影响，最终得到简化表达式为

式中，W[·]为窗函数的离散傅里叶变换形式；Δf为频率分辨率，Δf=Fs/N，kr=f0/Δf，其中kr为峰值频点，由于是离散采样，所以通常kr都是非整数。

对式(2)中Blackman-Harris 窗的时域表达式进行离散傅里叶变换，得到表达式为

式中，b1=0.358 75，b2=0.224 145，b3=0.070 64，b4=0.005 8。WR(·)为矩形窗的频谱幅度函数表达式。式(8)中W[·]的参数根据现有条件可以转换为[16]

WR(·)的表达式为

由式(8)可知ω的表达式(10)，将其代入式(12)进行离散采样

因式(9)中含有WR(·)形式，为利于后面计算所以写成WR(ω±2πm/N)的形式。

通常情况下对信号采样很难做到同步采样，所以f0表示的峰值频率通常不在离散的频谱点上，如图1 所示，在峰值频点kr两侧分别有两条最靠近kr频点的谱线，分别是kr1、kr2、kr3、kr4。分析单谱线、双谱线和三谱线算法可以知道，靠近峰值频点kr的谱线理论上拥有更多的频谱信息，所以要尽可能多地分析谱线信息。

图1 信号频谱图

引入变量ρ，ρ代表位于峰值频点kr两边的幅值之比，表达式为

将式(8)、(10)一并代入式(13)得到表达式为

从式(14)可以看出，ρ可以表示成γ的函数表达式，ρ=y(γ)，其中反函数可以表示为γ=y-1(ρ)，可以根据ρ求出γ，利用Matlab 中的polyfit 多项式曲线拟合函数拟合出表达式，在调用polyfit 函数时需要选择拟合多项式的系数n，n取值过大会增加运算成本，通过仿真遍历可知，当n取值大于7 后，多项式中阶数大于7 的项对计算结果影响已经很小，且会增加运算时间，所以通过polyfit 多项式曲线拟合函数拟合的通用表达式为

式中，c1、c3、c5、c2n+1分别代表每个奇数次项(2n+1)前面的系数。

其中信号频率的校正公式为

由条件可知测量信号的初始相位修正表达式为

测量信号的幅值是根据图1 中的四条谱线通过加权推算得来，且因为kr2、kr3是最靠近峰值频点的两条谱线，所以需要乘以更大的权值，这四条谱线的权值分别取1、3、3、1。由此得到信号幅值修正表达式为

其中，当频谱分析中N取值较大时，此时可以表达为

g(γ)经过拟合后的表达式可表示为

式中，b0、b2、b2n分别代表偶数次项(2n)前面的系数。

将式(9)和式(12)代入式(14)，利用Matlab 中的曲线拟合函数polyfit()进行求解，为保证数据的精度，在选择拟合次数时不能过小，此次拟合选取的次数为7，在γ取值范围内取值代入，通过函数polyfit(ρ，γ，7)拟合求得Blackman-Harris 窗四谱线插值算法中γ为

如图2 所示，图2a 为通过拟合得到的曲线图，图2b 为真实数值与拟合曲线数值的差值。

图2 拟合曲线和差值图

将式(9)、式(12)和式(18)结合，利用曲线拟合函数polyfit(γ，g(γ)，6)求得Blackman-Harris 窗四谱线插值算法中g(γ)为

将求得的式(21)、(22)代入式(16)、(17)、(19)，即可得到基于Blackman-Harris 窗函数四谱线插值算法的频率、相位和幅值。

3 主瓣干扰判定

3.1 频率相近分量频谱分析

一般测量得到的电信号中不仅含有整数次谐波，还含有间谐波，且当间谐波与整数次谐波相近时就会发生频谱混叠，导致的结果就是不能区分相近的频率分量。

设待测信号的表达式为

式中，A1、A2表示两个频率分量的幅值；f1、f2表示信号中含有的两个频率分量；φ1、φ2表示两个频率分量的初相位。相近频率频谱图如图3 所示。

图3 相近频率频谱图

由图3 可以看出，当两个频率相差比较近时，在频谱图中只能看到一个波峰，也就是发生了频谱混叠，只存在一个主瓣，且发生了主瓣干扰，使得被干扰频率的谱线幅值也发生了失真现象。如果采用加窗插值算法对测量信号进行分析处理，得到的数据准确度会产生偏差，且不能反映出真实的频谱情况。

3.2 主瓣干扰分析

现就需要对测量的电信号进行判定，判定的依据为是否发生了主瓣间的干扰。一般可以认为，当没有发生主瓣干扰时，且不受信号噪声影响时，主瓣内谱线间的相位和幅值有关联，主瓣内相邻谱线间的相位差绝对值等于π，所以在判定主瓣干扰时可以作为依据。

余弦组合窗是目前运用最多的一类窗函数，以此为例，其时域表达式为式(1)，可知其频谱表达式为

式中，Wo(·)为矩形窗的离散傅里叶变换形式。

将式(25)代入式(24)化简可得

将式(26)写成

设加窗信号经过离散傅里叶变换为X(k)，其中在测量第i项谐波时，若不考虑负频率和周边谐波的泄漏影响可知

设在主瓣内，k1和k2是相邻的两根谱线，k2=k1+1，可知

可知两条谱线的相位差为

由此可知在不受信号噪声及主瓣干扰的情况下，主瓣内相邻谱线间的相位差为π。

引入参数Δθ，根据表达式

式中，δ1、δ2分别是靠近峰值频点最近的两个谱线相位。按照实际工程应用来说，Δθ通常都不等于0，所以需要选定一个恰当的值£来权衡是否发生了主瓣干扰，£值的大小直接决定了测量时长和精准度。如果£的取值过大，就会导致发生了主瓣干扰但是没有判别出来的情况，从而导致漏判；若£的取值过小，将会导致过于灵敏，使得计算时长增加[17]。本文结合实际工程应用，对不同频率差信号分量进行仿真遍历，在引入较少运算量的同时对谐波分量进行有效分离检测，最终选取£=8×10-6。

4 复调制频谱细化算法

当测量电信号的某一频段上判定了主瓣干扰，说明在这一频段上存在不止一个频率分量，且它们的频率值相差很小，已经发生了频谱混叠的现象，这时就需要利用频率分辨率更高的处理算法来区分在这一频段上存在的频率分量[18]。本文选取的是复调制频谱细化算法ZoomFFT(或ZFFT)。

ZFFT 算法能够实现在某一频段上对频谱进行放大处理，使其在这一频段上频率分辨率提高，放大D倍[19]。对于频率分辨率Δf=Fs/N，根据公式可知想要提高频率分辨率，可以改变Fs和N这两个值。

(1) 保持N值不变，减小采样频率Fs的值，使得Δf减小，也就是频率分辨率提高，但是根据奈奎斯特采样定理可知，如果采样频率Fs小于两倍信号中的最高频率，测量的信号频谱将会发生频谱混叠现象，即要保证Fs≥2fmax。

(2) 保持采样频率Fs不变，增大采样点数N，也能提高频率分辨率，但是随着N值增大，带来的是计算量的增大，从而导致计算时长的增加，这对一些实时性要求高的场合就不太适用，且会占用大量内存。

ZoomFFT 算法实现的流程图如图4 所示。

图4 ZFFT 算法流程图

ZFFT 算法实现的过程总共分为7 个部分，其具体步骤如下所示。

(1) 滤波。根据奈奎斯特采样定理可知，FFT的频率分析范围有限，2/Fs为最高的分析频率，所以需要滤除其频率以上的频率，防止在分析频谱时出现混叠。

(2) 离散采样。对时域信号x(t)进行离散化处理，得到x(n)，采样长度为放大倍数D与采样点数的乘积。

(3) 频移。经过处理后的离散信号x(n)乘以旋转因子exp(-j2πfn)，目的是为了让所需分析的频段中心频率f移至频率零点，得到新的信号。

(4) 低通滤波。目的是为了滤除所需频谱细化频段以外的信号分量，只保留所需细化的频段。

(5) 重采样。对经过低通滤波器之后保留的频段进行重新采样，采样时每个点的间隔为D，采样频率由之前的Fs变为Fs/D，由此可见，采样频率降低了D倍。

(6) FFT 分析。对重采样后的离散信号进行FFT分析，分析的点数仍然为N，但是采样频率由之前的Fs变为Fs/D，所以频率分辨率也会提升D倍。

(7) 频率调整。由步骤(3)可知信号经过频移，所以需要对最终得到的频谱信号进行调整。

综上，通过ZFFT 算法处理过后的频谱在以频率f为中心的频段上，将频率分辨率由之前的Δf=Fs/N降为Δf=Fs/ND，也就是将频率分辨率提升了D倍。放大的频段范围是ΔF=f±Fs/2D。

观察ZFFT 算法流程步骤可知，ZFFT 算法只能针对想要观察的目标频段进行放大处理，不能够对全频域内进行放大。结合如今的电力系统测量分析可知，大多分析为50 次谐波以内，所以一般分析0～2.5 kHz 谐波含量，若只对存在主瓣干扰的频段采用频谱细化分析，对只含有整数次谐波的频段采用加窗插值算法分析，这将会大大减少计算时间，不仅提高了谐波检测精度，同时也提高了实时性。放大倍数的选择也决定了ZFFT 算法的有效性，如果放大倍数过小，则不能够对存在主瓣干扰的频率实现有效区分，若放大倍数过大，则会导致采样的长度过长。

综上所述，对测量电能质量提出了改进方法，流程如图5 所示。

图5 改进型测量方法流程图

根据图5 可知，结合Blackman-Harris 窗插值算法和ZFFT 算法，利用条件主瓣干扰判定选用算法，弥补了ZFFT 算法不能全频域对频谱进行细化，基于Blackman-Harris 窗四谱线插值不能够精确区分频率相近分量的缺点。

5 仿真验证与分析

为了验证所提改进型方法的准确度和实时性，这里对整数次谐波和间谐波进行了频谱分析。主要的电力参数设置如表2 所示。

表2 电力参数

结合实际电力系统工程分析，此次试验分析的谐波范围为0～50 次，相应频谱范围为0～2.5 kHz，最高的频率为2.5 kHz。因此，根据奈奎斯特采样定理选取采样频率为5 kHz，采样点数N=1 024，ZoomFFT 选取的放大倍数D=10。

此外，在表2 中，参数中除了含有基波和其他整数次谐波含量，还设置了间谐波且频率相近。

图6 为改进型方法仿真局部频谱图。根据仿真结果实测数据得到参数结果如表3 所示。

图6 改进型方法仿真频谱图

表3 改进型方法仿真结果数据

同时分别对传统的 FFT 谐波检测方法和Blackman-Harris 窗插值算法进行仿真验证对比，图7为基于Blackman-Harris 窗四谱线插值算法仿真频谱图，图8 为FFT 算法仿真频谱图。信号仿真设置参数如表2 所示。

图7 基于Blackman-Harris 窗四谱线插值算法仿真频谱图

图8 传统FFT 算法仿真频谱图

对比图6、图7、图8 和表3、表4 可知，信号分量X2在基于Blackman-Harris 窗四谱线插值算法中无法被检测出，不能将所有的间谐波区分开来，同样条件下，采用本文所提方法能够将所有间谐波全部区分并取得精确的参数值。

表4 基于Blackman-Harris 窗四谱线插值算法仿真结果数据

通过对比不同谐波检测分析的方法可知，对电信号加不同的窗函数和采用不同的插值算法得到的时间成本都不相同[20]，其中ZoomFFT 算法在提高频率分辨时是以分析更多的采样点数得到的。如果只使用ZFFT 算法对信号进行分析，将会耗费大量的时间和内存单元。通过仿真对比，不同加窗插值算法的平均计算时长如表5 所示。由表5 可知，ZFFT算法精确度高，但计算时间过长，不利于测量要求的实时性，且会占据大量内存。因此，本文采用的方法兼容了ZFFT 和加Blackman-Harris 窗插值算法的优点。