CN110837001B - 一种电力系统中谐波和间谐波的分析方法与装置 - Google Patents

一种电力系统中谐波和间谐波的分析方法与装置 Download PDF

Info

Publication number
CN110837001B
CN110837001B CN201911108332.0A CN201911108332A CN110837001B CN 110837001 B CN110837001 B CN 110837001B CN 201911108332 A CN201911108332 A CN 201911108332A CN 110837001 B CN110837001 B CN 110837001B
Authority
CN
China
Prior art keywords
inter
harmonic waves
harmonic
harmonics
frequency
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201911108332.0A
Other languages
English (en)
Other versions
CN110837001A (zh
Inventor
蔡晓峰
张鸿博
窦小磊
王佳欣
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Henan Institute of Engineering
Original Assignee
Henan Institute of Engineering
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Henan Institute of Engineering filed Critical Henan Institute of Engineering
Priority to CN201911108332.0A priority Critical patent/CN110837001B/zh
Publication of CN110837001A publication Critical patent/CN110837001A/zh
Application granted granted Critical
Publication of CN110837001B publication Critical patent/CN110837001B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01RMEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R23/00Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
    • G01R23/16Spectrum analysis; Fourier analysis
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/14Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141Discrete Fourier transforms

Landscapes

  • Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Discrete Mathematics (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Complex Calculations (AREA)

Abstract

本发明涉及一种电力系统中谐波和间谐波的分析方法与装置,属于电网谐波分析领域,方法包括获取电力系统信号,采用Root‑MUSIC算法求出信号的谐波和间谐波频率,根据谐波和间谐波中相邻两个频率之间的频率间隔和设定的间隔阀值对谐波和间谐波进行分组;对于加窗处理后的信号,利用离散傅里叶变换进行谱分析,得到加窗处理后信号的频谱;当分组中谐波和间谐波数量M≥2时,选取与分组中各谐波和间谐波距离最近的M条谱线,通过构建谱线方程组求解谐波和间谐波的幅值和相位。本发明估计出谐波和间谐波的频率后,根据频率间隔对(间)谐波进行分组求解,运算量较小,能有效克服频率相近的间谐波间以及谐波和间谐波之间的主瓣干扰,从而获得较高精度的幅值和相位。

Description

一种电力系统中谐波和间谐波的分析方法与装置
技术领域
本发明属于电网谐波分析领域,具体涉及一种电力系统中谐波和间谐波的分析方法与装置。
背景技术
近年来随着新能源发电和电力电子技术在电力系统中的广泛应用,电网中出现大量谐波和间谐波,而于此同时,各种精密的、对电能质量敏感的用电设备也在不断普及,使得谐波和间谐波问题更加突出,准确进行电网(间)谐波分析是电能质量监测、谐波治理、电能计量等技术的前提,下面对现有技术中谐波和间谐波分析的几种常用方法进行介绍:
FFT(Fast Fourier Transform,快速傅里叶变换)是目前(间)谐波分析的常用方法,但在非同步采样时会出现频谱泄漏和栅栏效应问题。加窗插值FFT算法可在一定程度上抑制非同步采样引起的频谱泄漏和栅栏效应,但当谐波与间谐波或不同频率的间谐波之间频率相近时,严重的频谱干扰影响加窗插值FFT算法对谐波、间谐波参数的精确估计。虽然采用上述优良的窗函数可以抑制频谱泄露,但一般性能优良的窗函数主瓣也较宽,容易造成频率相近的间谐波或谐波与间谐波之间的主瓣重叠,此时通过常规的插值方式进行间谐波的分析已无法避免主瓣间的干扰。
多重信号分类法(Multiple Signal Classification,MUSIC)是一种估计信号空间参数的现代谱估计方法,其基本思想是对由采样数据构成的自相关矩阵进行特征分解,从而求得信号子空间和噪声子空间,然后利用信号子空间和噪声子空间的正交性来估计信号的频率,它具有较高的分辨率、估计精度及稳定性。Root-MUSIC算法是对MUSIC改进,以对多项式求根计算代替在频域轴上的谱峰搜索,减少计算量,同时提高了精度。但是MUSIC和Root-MUSIC算法都只能够检测信号中间谐波的频率,无法估计间谐波的幅值和相位。作者李新在期刊《电力自动化设备》中发表了论文《基于实值Root-MUSIC和Prony算法的间谐波参数估计》,该方法先用Root-MUSIC算法估计出谐波间谐波的频率,然后利用Prony算法估计出幅值和相角,但Prony算法易受噪声影响,且只能将所有谐波间谐波一并求出,而不能有选择地求解特定谐波或间谐波信号,导致谐波间谐波含量较多、采样点数较大时算法运算量大。
发明内容
本发明的目的是提供一种电力系统中谐波和间谐波的分析方法与装置,用于解决现有技术对谐波和间谐波参数估计的精度低且计算量大的问题。
基于上述目的,一种电力系统中谐波和间谐波的分析方法的技术方案如下:
获取电力系统的信号,采用Root-MUSIC算法求出信号的谐波和间谐波频率,根据谐波和间谐波中相邻两个频率之间的频率间隔和设定的间隔阀值,将谐波和间谐波进行分组;
采用窗函数对所述信号进行加窗处理,利用离散傅里叶变换对加窗处理后的信号进行谱分析,得到加窗处理后信号的频谱;
当所述分组中的谐波和间谐波数量为M,且M≥2时,选取与所述分组中各谐波和间谐波距离最近的M条谱线,通过构建以下谱线方程组并求解,获得谐波和间谐波的幅值和相位;
X=WA
方程组中,矩阵X、W、A是复数向量或矩阵,其中X为所述M条谱线对应的复数值构成的向量,W为对应M条谱线的频谱函数值构成的矩阵,
Figure BDA0002271972690000021
m=1,2,…,M,其中A和θ分别为分组中待求谐波和间谐波的幅值和相位。
基于上述目的,一种电力系统中谐波和间谐波的处理装置的技术方案如下:
包括处理器,用于执行指令以实现上述电力系统中谐波和间谐波的分析方法。
上述两个技术方案的有益效果是:
本发明先利用Root-MUSIC算法估计出谐波和间谐波的频率,然后利用谐波和间谐波的频率分布密集程度,根据谐波和间谐波中相邻两个频率之间的频率间隔和设定的间隔阀值对(间)谐波进行分组求解,每组中矩阵W的阶数均很小,因此能够很快获得谐波和间谐波的幅值和相位。与现有技术求解大型方程组相比,本发明每一组分组后的矩阵阶数远小于现有技术中求解时矩阵的阶数,运算量较小;且能有效克服频率相近的间谐波之间以及谐波和间谐波之间的主瓣干扰,从而获得较高精度的幅值和相位。另外,与现有技术相比,本发明通过Root-MUSIC算法结合构造谱线方程组的方法,在谐波和间谐波幅值的估计上具有更好的抗噪声性能。
对于只有一个谐波或只有一个间谐波的分组,由于已经采用窗函数进行加窗处理,远距离的频率成分产生的频谱泄露被大幅抑制,因此直接采用如下谐波和间谐波的幅值和相位计算式:
A=2Xw(k1r)|/|W(2π(k1r-k1)/N)
θ=arg[2jXw(k1r)/W(2π(k1r-k1)/N)]
式中,A为待求谐波或间谐波的幅值,Xw(k1r)为距离原始信号角频率ω1最近的第k1r条谱线的值,k1为角频率ω1对应的谱线坐标,W为第k1r条谱线的频谱函数值,N为对加窗处理后的信号进行利用离散傅里叶变换的抽样点数;θ为待求谐波或间谐波的相位,arg表示求复数的相位。
当所述窗函数采用组合余弦窗函数,如Nuttall窗时,所述设定的间隔阀值是根据窗函数对应的主瓣宽度确定的。
进一步,若分组内的谐波和/或间谐波的频率相近,那么,与分组中各谐波和间谐波距离最近的M条谱线可能为同一条。基于此,当与所述分组中各谐波和间谐波距离最近的M条谱线为同一条谱线K时,剩余的M-1条谱线从谱线K的相邻谱线中选择。
附图说明
图1是本发明方法实施例中对谐波和间谐波进行分组的示意图;
图2是本发明方法实施例中信号的快速傅里叶变换的幅度谱图;
图3是本发明方法实施例中的谐波和间谐波分析方法流程图。
具体实施方式
下面结合附图对本发明的具体实施方式作进一步的说明。
方法实施例:
本实施例提出一种电力系统中的谐波和间谐波分析方法,如图3所示,包括以下步骤:
一、先用Root-MUSIC算法(求根多重信号分类算法)估计出信号的(间)谐波频率,根据(间)谐波频率之间的频率间隔和设定的间隔阀值,对(间)谐波进行分组。
具体的,设电力系统的谐波和间谐波采样信号(电力信号)表示为:
Figure BDA0002271972690000041
其中,x(n)为采样信号;n为采样点数;Ts为采样周期;p为谐波和间谐波的个数;fi为第i个谐波或间谐波分量的频率;Ai、θi分别为第i个谐波或间谐波分量的幅值和初始相位;v(n)为零均值且方差为σ2的高斯白噪声。
截取一段长度为N的信号采样序列,构造数据阵列X:
Figure BDA0002271972690000042
首先求矩阵X的协方差矩阵Rx,并对矩阵Rx进行特征值分解,得到从小到大排列的m个特征值为:λ1≈λ2≈…≈λm-2pm-2p+1<…<λm。其中λ1…λm-2p这m-2p个特征值明显小于λm-2p+1…λm这2p个特征值,则可判断信号中谐波和间谐波的个数为p。较大的2p个特征值对应的特征向量张成信号子空间为US,较小的m-2p个特征值对应的特征向量张成噪声子空间为UN
为了利用噪声特征向量提取空间参数信息,构造Root-MUSIC型多项式:
Figure BDA0002271972690000043
其中,f(z)为多项式,p(z)=[1,z,…,zm-1]T,z为复变量,上标H表示共轭转置,UN为m-2p个特征值对应的特征向量张成噪声子空间。根据该式可得出,该式的零点即对应着信号的频率,但式并不完全为z的多项式,还存在z的共轭的幂次项,因此使方程f(z)=0的求根过程变得复杂,基于此可将上式修正为:
Figure BDA0002271972690000044
该式即为Root-MUSIC算法的多项式,此多项式的最高次数为2(m-1),即有(m-1)对根,该(m-1)对根中有p对根正好分布在单位圆上,这p对根就对应于采集信号的p个频率,但在实际应用中由于存在多方面的误差,p对根不一定正好位于单位圆上,此时只需求出p对最接近于单位圆的根即可。设最靠近单位圆的根为z,则对应采样信号(该信号为正弦信号)的频率为f=arg(z/(2πTs)),其中arg表示求复数的相位,Ts表示采样时间间隔。
根据欧拉公式,一个实正弦可以表示成两个复正弦的和或差,因此应用Root-MUSIC算法求得采样信号的频率为正负成对的形式,取其中正频率作为谐波和间谐波的频率估计值即可。
Root-MUSIC算法能够检测信号中谐波和间谐波的各个频率,为了计算谐波与间谐波的幅值和相位,对谐波和间谐波进行分组,即在频率间隔在间隔阀值6Δf以上的相邻两个频率成分之间划一个分界线,把谐波和间谐波分成若干组。
之所以采用6Δf(Δf为离散傅立叶变换对应的谱线间隔,即频率分辨率)作为划分宽度(即设定的间隔阀值),是因为6Δf对应为下面步骤二中加窗FFT所加6项余弦窗函数的主瓣宽度,如果两种频率成分相隔6Δf以上,说明这两种频率成分均已处于对方的窗函数主瓣之外,由于6项余弦窗优良的性能,主瓣以外的泄露频谱被大大抑制,因此可以忽略不同分组之间的泄露干扰。举例来说,假设一信号中含有f1-f6的6种频率成分,如图1所示,根据该图中的频率间隔关系,则f1-f3分为一组,f4-f5分为一组,f6为一组。
二、然后利用加窗FFT对采集信号进行谱分析。即采用窗函数对信号进行加窗处理,利用离散傅里叶变换对加窗处理后的信号进行谱分析,得到加窗处理后信号的频谱。
具体的,对于频率为f1、幅值为A、初相位为θ的单一频率信号x(t),经过采样率为fs的模数变换后得到离散信号x(n):
x(n)=Asin(2πf1n/fs+θ)=Asin(ω1n+θ) (5)
其中,ω1=2πf1/fs,为信号x(n)对应的数字角频率。
定义长度为N的组合余弦窗的时域形式为w(n),其连续频谱为W(2πf),则采用窗函数进行加窗处理后的信号记为xw(n),若忽略xw(n)负频点-ω1处频峰的旁瓣影响,则xw(n)在正频点ω1附近的连续频谱函数Xw(ω)可以表达为:
Figure BDA0002271972690000051
其中,j为复数符号;对上式进行N点离散抽样,即可得到Xw(ω)的N点离散傅立叶变换的表达式Xw(kΔω)为:
Figure BDA0002271972690000052
其中,Δω为离散频率间隔,N为数据截断长度,数字频率间隔Δω=2π/N,对应的模拟频率间隔为Δf=fs/N,k为谱线坐标,k=0,1,2,…,N-1,k1为频率ω1对应的谱线坐标,k1=ω1/Δω=f1/Δf=Nf1/fs
由于x(n)的频率ω1很难正好位于离散谱线的频点上,因此,k1一般不是整数;假设距离ω1最近的谱线为第k1r条,则有:
Figure BDA0002271972690000061
Ae=2jXw(k1r)/W(2π(k1r-k1)/N) (9)
其中,A为待求谐波或间谐波的幅值,Xw(k1r)为距离原始信号角频率ω1最近的第k1r条谱线的值,k1为角频率ω1对应的谱线坐标,W为第k1r条谱线的频谱函数,N为对加窗处理后的信号进行利用离散傅里叶变换的抽样点数;θ为待求谐波或间谐波的相位.
结合上面两个公式,得到信号幅值A为:
A=2|Xw(k1r)|/|W(2π(k1r-k1)/N)| (10)
信号初相位θ为:
θ=arg[2jXw(k1r)/W(2π(k1r-k1)/N)] (11)
其中,||表示求复数的模,arg表示求复数的相位。
为了较好地抑制非同步采样下的频谱泄露,在进行离散傅立叶变换之前,采用窗函数对原始信号在时域进行了加窗处理,本实施例采用组合余弦窗,N点组合余弦窗的计算式为:
Figure BDA0002271972690000062
其中i=0,1,2,...,N-1,D为窗函数的项数,ad为窗函数的系数。
N点组合余弦窗频谱W(ω)可表示为:
Figure BDA0002271972690000063
其中:
Figure BDA0002271972690000071
为了便于分析,令ω=2πδ/N化简可得:
Figure BDA0002271972690000072
Figure BDA0002271972690000073
根据式(6)和式(15),令δ=k-k1,可得加窗FFT进行谱分析所得离散谱线可表示为:
Figure BDA0002271972690000074
本实施例采用的余弦窗为6项余弦窗,其对应的窗函数的系数分别为:
a0=0.24609375,a1=0.41015625,a2=0.234375,a3=0.087890625,a4=0.01953125,a5=0.001953125。
表1对比分析了6项余弦窗与常用窗函数的旁瓣特性参数,可以看出6项余弦窗具有旁瓣峰值小且衰减速度快的良好性能。
表1
Figure BDA0002271972690000075
三、对于谐波和间谐波密集的分组,通过构建谱线方程组并求解获得谐波和间谐波参数,即幅值和相位。具体过程如下:
当分组中含有多个频率成分时,采用加窗FFT进行谱分析的结果是单个频率成分频谱的叠加,例如x(n)为3个频率成分的叠加:
x(n)=A1sin(ω1n+θ1)+A2sin(ω2n+θ2)+A3sin(ω3n+θ3) (17)
则根据式(7)有:
Figure BDA0002271972690000081
式中k1=ω1/Δω,k2=ω2/Δω,k3=ω3/Δω,找出h1、h2、h3三条谱线结合上式(18)可列写三个方程构成一个方程组,写成矩阵形式为:
Figure BDA0002271972690000082
式中,h1,h2,h3表示选取的三条谱线,且三条谱线的位置坐标依次为h1、h2、h3,X(h1),X(h2),X(h3)分别对应三条谱线的值,A1,A2,A3为该分组中待求的三个幅值,θ123为该分组中待求的三个相位。
式(19)中,W是一个矩阵,矩阵中每个元素都是由N点组合余弦窗频谱
Figure BDA0002271972690000083
即上面式(15)计算得来的,只是每个元素对应的δ值不同。式(19)简记为:X=WA,其中X=2j×[X(k1) X(k2) X(k3)]T
上式中向量X中的元素可由加窗FFT得到,矩阵W由于Root-MUSIC算法已经求得信号中包含的频率值也可以计算得到,因此仅向量A是未知量,于是可求得A=W-1X,W-1为W的逆阵。A中元素的模对应各频率成分的幅值,相角对应各频率成分的初相角。
以上以x(n)中含有3个频率成分为例进行了说明,当x(n)中含有任意M个频率成分(M为大于等于2的整数)时,即:
x(n)=A1sin(ω1n+θ1)+A2sin(ω2n+θ2)…AMsin(ωMn+θM) (20)
则根据式(7)有:
Figure BDA0002271972690000084
式中k1=ω1/Δω,k2=ω2/Δω,…,,kM=ωM/Δω,找出h1、h2、…,,hM共M条谱线结合上式(20)可列写M个方程构成一个方程组,写成矩阵形式为:
Figure BDA0002271972690000085
式中,h1,h2,…,hM表示选取的M条谱线,且M条谱线的位置坐标依次为h1、h2、…,hM,X(h1),X(h2),…,X(hM)分别对应M条谱线的值,A1,A2,…,AM为该分组中待求的M个幅值,θ12,…,θM为该分组中待求的M个相位。
需要说明的是,本实施例中的分组是至少两个频率成分的谐波和/或间谐波的分组,分组中既可以只有谐波,也可以只有间谐波,还可以谐波和间谐波均有,因此每个分组中的频谱阶数与分为一组的谐波和间谐波数量相同,至少为两阶(即2×2矩阵)。
本实施例中,对于只有一种频率成分的谐波和间谐波组,由于已经采用性能优良的窗函数进行加窗处理,远距离的频率成分产生的频谱泄露被大幅抑制,因此可直接根据式(10)、式(11)计算出谐波和间谐波参数,而对于含有两种及以上频率成分的分组,,则采用上述构建谱线方程组的方法(即步骤三中的方法)求解。
本实施例中,在谱线的选择上,首先选择距离分组内各次(间)谐波频率最近的谱线,因为时域信号经过傅立叶变换后信号频率点附近的谱线幅值大,抗噪声能力强。如果间谐波频率相近,距离它们最近的谱线可能是同一条谱线,如谱线K,那么剩余的谱线就从谱线K的相邻谱线中选择。这样做,不必对每种频率成分都采用解方程组的方式求解,而是仅在需要时进行,同时采用分组求解的方式也避免了求解大型方程组,相对现有技术的算法(如Prony算法)而言降低了运算量。而Prony算法只能将所有谐波和间谐波一并求出,需要分析的谐波和间谐波次数较多、采样点数较大,算法运算量大,相对于本发明的方法不能有选择地对谐波和间谐波分组求解。
此外,本实施例中的方程组是在频域构造和求解的,由于采用了加6项余弦窗处理,6项余弦窗具有优良的频谱泄露抑制能力,距离较远的其他频率成分的频谱泄露被大大抑制,为获得高精度的间谐波估计结果提供了条件。
下面将本方案与现有技术中的Root-MUSIC+Prony算法的运算量相比,二者均需要Root-MUSIC算法求解谐波间谐波的频率,因此频率估计上二者的运算量相同,但在幅值和相位的估计上,二者运算量存在很大差异,Prony算法需要根据频率估计结果构造范德蒙矩阵,然后根据最小二乘法求解,即:
a=(ZHZ)-1ZHx
其中,a为待求信号幅值和相位构成的向量,x为采样点构成的向量,Z为范德蒙矩阵,H表示共轭转置,其阶数p等于信号中谐波间谐波数量,由于1个实正弦可以表示成2个复正弦的和或差,因此对于实信号,Z的阶数p是信号中谐波间谐波数量的2倍。。根据电网谐波测量标准,通常需要测量25次谐波,再加上测量间谐波的需要,Z的阶数p不宜小于50阶,因此涉及到一个50阶矩阵的求逆运算,另外还涉及多次矩阵乘法运算,具体为:
ZHZ是p×N和N×p矩阵相乘(可得p×p矩阵,矩阵每个元素均需N次乘法,共p×p×N次乘法),ZHZ的逆阵再乘ZH是p×p的矩阵和p×N矩阵相乘(可得p×N矩阵,而矩阵每个元素均需p次乘法,共p×N×p次乘法),最后(ZHZ)-1ZH和向量x相乘是p×N矩阵和一个N×1向量相乘(p×N次乘法),N为采样点数,最终,Prony算法估计谐波间谐波的幅值和相位共需2×p×p×N+p×N次乘法和一个p×p阶矩阵求逆。
而本发明中的方法需要一次快速傅里叶变换运算的N/2×log2N次乘法,假设每周期采样102点(理论上可分析51次谐波),共采样1024点,快速傅里叶变换仅需5×N次乘法,远小于Prony法的2×p×p×N+p×N(p>50)次乘法。另外,本发明的方法根据频率间隔对谐波和间谐波进行分组求解,构造的谱线方程组规模一般较小,求解运算量远小于p×p(p>50)阶矩阵求逆的运算量。因此,本发明的方法在谐波和间谐波的幅值与相位的估计上的运算量远小于Prony算法。
下面按照本发明的方法进行仿真如下分析:
仿真信号包括谐波成分和间谐波成分,表达式如下所示:
Figure BDA0002271972690000101
式中,p为谐波个数,Ai为i次谐波幅值,f1为基波频率,θi为i次谐波初相位;q为间谐波个数,fj为j号间谐波的频率,Aj为j号间谐波的幅值,θj为j号间谐波初相位。
根据电网谐波测量标准要求,需要测量40次谐波,但40次谐波的参数太多,由篇幅所限,下面仅考虑一些低次谐波和间谐波,各频率成分参数如下表所示:
表2
Figure BDA0002271972690000102
Figure BDA0002271972690000111
采样频率设为5120Hz,采样10个额定工频周期(0.2s)共得1024点。
信号的快速傅里叶变换的幅度谱如图2所示,由于53Hz间谐波距离基波很近,幅值又远小于基波,造成二者主瓣重叠,且53Hz间谐波被基波主瓣严重干扰,已经无法从幅度谱上观察出来。同时,由于53Hz间谐波的存在,对基波频率的测量也产生了很大影响,按照作者徐艳春在期刊《电力系统保护与控制》中论文《基于六项余弦窗四谱线插值FFT的高精度谐波检测算法》提供的加窗插值原理,计算基波频率为50.2562Hz,误差达到0.0562Hz,由于在快速傅里叶变换的谐波分析中通常以i×f0作为高次谐波的频率,因此基波频率的偏差无疑会给高次谐波的幅值、相位校正带来更大误差,因此,对于存在间谐波尤其是基波附近存在间谐波的信号,按照插值FFT算法进行谐波和间谐波分析效果不理想。而66Hz间谐波在幅度谱上虽然还有所显现,但也存在被基波主瓣干扰的问题,不能利用插值FFT算法获取高精度的参数估计结果。
本实施例中Root-MUSIC算法的频率估计结果如表3所示,表中的误差为相对误差,aE±n表示a×10±n(下同),从表3中可以看出,Root-MUSIC算法的频率分辨率高,能相当准确地分析出谐波间谐波的频率,为接下来分析间谐波的幅值和相位奠定了基础;不过阵元数M直接影响到Root-MUSIC算法精度,较大的M值能带来更高的精度,但需要的运算时间长。通过仿真发现,在采样频率为5120Hz、采样1024点的情况下,取M=400,可在保证精度的同时获得较快的运算速度,而且有较好的抗噪能力,因此在Root-MUSIC算法中阵元数M取400。
无噪情况下Root-MUSIC算法频率估计结果如下表所示:
表3
Figure BDA0002271972690000112
Figure BDA0002271972690000121
为对比分析谐波间谐波的幅值、相位测量精度,采用作者李新在期刊《电力自动化设备》中论文《基于实值Root-MUSIC和Prony算法的间谐波参数估计》提供的Root-MUSIC+Prony算法(称算法1)和本发明的方法(称算法2)分别进行了谐波和间谐波的仿真分析,需要说明的是,上述论文已仿真对比了几种常用的间谐波分析算法,验证了Root-MUSIC+Prony算法具有较高精度,因此这里不再考虑其他方法,仿真结果如表4所示,无噪情况下的幅值、相位估计结果如下表所示:
表4
Figure BDA0002271972690000122
表4表明,在无噪情况下,算法1和算法2对于幅值和相位的估计均达到了很高的精度,且算法2对于大部分谐波和间谐波的幅值测量精度更高,虽然算法2对2次、4次谐波的幅值测量精度略低于算法1,但差别并不大;在相位测量精度上,算法1的精度优于算法2,容易看出两种算法都有精度高于对方的成分,也都有精度不如对方的成分,相互之间的优势不甚明显。
实际采集的电力信号中往往含有噪声,在含义噪声情况下对信号进行仿真的结果如表5、表6所示,其中表5为信噪比为60dB时算法1和算法2的仿真结果,表6为信噪比为40dB时算法1和算法2的仿真结果。
表5
Figure BDA0002271972690000131
表6
Figure BDA0002271972690000132
由表5、表6表明,有噪声情况下两种算法精度均有所下降,且随着信噪比的降低,算法精度也随之降低(由于噪声条件下测量误差有一定的随机性,因此表中误差均取20次测量的平均值)。进一步分析表5、表6不难发现,算法2在信噪比为60dB和40dB情况下对于谐波和间谐波幅值的测量精度高于算法1,表现出更好的抗噪声性能。但同时从表5、表6中两种算法相位误差的对比来看(表5、表6中加粗字体表示的项为算法1精度优于算法2,非加粗字体表示的项为算法2精度优于算法1的成分),两种算法各都有精度高于对方的成分,也有精度不如对方的成分,相互之间的优势不明显。不过在电力谐波、间谐波的分析中,往往更加关注谐波和间谐波的幅值大小,在这些应用场合,采用算法2有利于获得更精确的幅值测量结果。
综上,本发明的方法具有以下优点:
1、本发明提出的基于Root-MUSIC算法与谱线方程组的谐波和间谐波分析方法,先利用Root-MUSIC算法估计出谐波和间谐波的频率,然后根据谐波和间谐波之间的频率间隔对谐波和间谐波进行分组求解,避免了求解大型方程组,运算量较小。
2、本实施例中的谱线方程组是在加6项余弦窗FFT的基础上建立的,由于6项余弦窗具有优良的频谱泄露抑制能力,距离较远的其他频率成分的频谱泄露被大大抑制,为获得高精度的谐波和间谐波估计结果提供了条件。
3、对于间谐波密集的分组,本发明通过构建谱线方程组并进行求解,从而获得间谐波参数,能有效克服频率相近的间谐波之间的主瓣干扰,获得较高的分析精度。
4、本发明的谐波和间谐波分析方法,与现有Root-MUSIC+Prony算法的间谐波分析方法相比,本发明的谐波和间谐波分析方法在谐波和间谐波幅值的估计上表现出更好的抗噪声性能。
装置实施例:
本实施例提出一种电力系统中谐波和间谐波的处理装置,包括处理器,用于执行指令以实现以下步骤:
获取电力系统的信号,采用Root-MUSIC算法求出信号的谐波和间谐波频率,根据谐波和间谐波中相邻两个频率之间设定的频率间隔,将谐波和间谐波进行分组;
采用窗函数对信号进行加窗处理,利用离散傅里叶变换对加窗处理后的信号进行谱分析,得到加窗处理后信号的频谱;
当分组中的谐波和间谐波数量为M,且M≥2时,选取与分组中各谐波和间谐波距离最近的M条谱线,通过构建以下谱线方程组并求解,获得谐波和间谐波的幅值和相位;
X=WA
方程组中,矩阵X、W、A均是复数矩阵,其中X为选取M条谱线的频谱函数值,W为对应M条谱线的频谱函数值构成的矩阵,
Figure BDA0002271972690000151
m=1,2,…,M,其中A和θ分别为分组中待求谐波和间谐波的幅值和相位。
上述实施例中所指的处理装置,实际上是基于本发明方法流程的一种计算机解决方案,即一种软件构架,上述装置即为与方法流程相对应的处理进程。由于对上述方法已经在方法实施例中介绍的足够清楚完整,故不再详细进行描述。
另外,本实施例中的处理器既可以是计算机,也可以是微处理器,如ARM等,还可以是可编程芯片,如FPGA、DSP等。
以上所述仅为本发明的优选实施例,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。例如,本实施例中也可以采用旁瓣衰减更加迅速的窗函数,对应地,在间谐波分组时根据窗函数对应的主瓣宽度调整分组需要的频率间隔。
因此,凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的权利要求范围之内。

Claims (8)

1.一种电力系统中谐波和间谐波的分析方法,其特征在于,包括以下步骤:
获取电力系统的信号,采用Root-MUSIC算法求出信号的谐波和间谐波频率,根据谐波和间谐波中相邻两个频率之间的频率间隔和设定的间隔阀值,将谐波和间谐波进行分组;
采用窗函数对所述信号进行加窗处理,利用离散傅里叶变换对加窗处理后的信号进行谱分析,得到加窗处理后信号的频谱;
当所述分组中的谐波和间谐波数量为M,且M≥2时,选取与所述分组中各谐波和间谐波距离最近的M条谱线,通过构建以下谱线方程组并求解,获得谐波和间谐波的幅值和相位;
X=WA
方程组中,矩阵X、W、A是复数向量或矩阵,其中X为所述M条谱线对应的复数值构成的向量,W为对应M条谱线的频谱函数值构成的矩阵,
Figure FDA0002271972680000011
其中A和θ分别为分组中待求谐波和间谐波的幅值和相位。
2.根据权利要求1所述的电力系统中谐波和间谐波的分析方法,其特征在于,对于分组中只有单一频率的谐波或间谐波,其幅值和相位计算式如下:
A=2|Xw(k1r)|/|W(2π(k1r-k1)/N)|
θ=arg[2jXw(k1r)/W(2π(k1r-k1)/N)]
式中,A为待求谐波或间谐波的幅值,Xw(k1r)为距离原始信号角频率ω1最近的第k1r条谱线的值,k1为角频率ω1对应的谱线坐标,W为第k1r条谱线的频谱函数值,N为对加窗处理后的信号进行利用离散傅里叶变换的抽样点数;θ为待求谐波或间谐波的相位,arg表示求复数的相位。
3.根据权利要求1所述的电力系统中谐波和间谐波的分析方法,其特征在于,当所述窗函数为组合余弦窗函数时,所述设定的间隔阀值是根据窗函数对应的主瓣宽度确定的。
4.根据权利要求1所述的电力系统中谐波和间谐波的分析方法,其特征在于,当与所述分组中各谐波和间谐波距离最近的M条谱线为同一条谱线K时,剩余的M-1条谱线从谱线K的相邻谱线中选择。
5.一种电力系统中谐波和间谐波的处理装置,其特征在于,包括处理器,用于执行指令以实现以下步骤:
获取电力系统的信号,采用Root-MUSIC算法求出信号的谐波和间谐波频率,根据谐波和间谐波中相邻两个频率之间的频率间隔和设定的间隔阀值,将谐波和间谐波进行分组;
采用窗函数对所述信号进行加窗处理,利用离散傅里叶变换对加窗处理后的信号进行谱分析,得到加窗处理后信号的频谱;
当所述分组中的谐波和间谐波数量为M,且M≥2时,选取与所述分组中各谐波和间谐波距离最近的M条谱线,通过构建以下谱线方程组并求解,获得谐波和间谐波的幅值和相位;
X=WA
方程组中,矩阵X、W、A是复数向量或矩阵,其中X为所述M条谱线对应的复数值构成的向量,W为对应M条谱线的频谱函数值构成的矩阵,
Figure FDA0002271972680000021
其中A和θ分别为分组中待求谐波和间谐波的幅值和相位。
6.根据权利要求5所述的电力系统中谐波和间谐波的处理装置,其特征在于,对于分组中只有单一频率的谐波或间谐波,其幅值和相位计算式如下:
A=2|Xw(k1r)|/|W(2π(k1r-k1)/N)|
θ=arg[2jXw(k1r)/W(2π(k1r-k1)/N)]
式中,A为待求谐波或间谐波的幅值,Xw(k1r)为距离原始信号角频率ω1最近的第k1r条谱线的值,k1为角频率ω1对应的谱线坐标,W为第k1r条谱线的频谱函数值,N为对加窗处理后的信号进行利用离散傅里叶变换的抽样点数;θ为待求谐波或间谐波的相位,arg表示求复数的相位。
7.根据权利要求5所述的电力系统中谐波和间谐波的处理装置,其特征在于,当所述窗函数为组合余弦窗函数时,所述设定的间隔阀值是根据窗函数对应的主瓣宽度确定的。
8.根据权利要求5所述的电力系统中谐波和间谐波的处理装置,其特征在于,当与所述分组中各谐波和间谐波距离最近的M条谱线为同一条谱线K时,剩余的M-1条谱线从谱线K的相邻谱线中选择。
CN201911108332.0A 2019-11-13 2019-11-13 一种电力系统中谐波和间谐波的分析方法与装置 Expired - Fee Related CN110837001B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911108332.0A CN110837001B (zh) 2019-11-13 2019-11-13 一种电力系统中谐波和间谐波的分析方法与装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911108332.0A CN110837001B (zh) 2019-11-13 2019-11-13 一种电力系统中谐波和间谐波的分析方法与装置

Publications (2)

Publication Number Publication Date
CN110837001A CN110837001A (zh) 2020-02-25
CN110837001B true CN110837001B (zh) 2021-10-01

Family

ID=69576264

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911108332.0A Expired - Fee Related CN110837001B (zh) 2019-11-13 2019-11-13 一种电力系统中谐波和间谐波的分析方法与装置

Country Status (1)

Country Link
CN (1) CN110837001B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2022016300A1 (zh) * 2020-07-23 2022-01-27 刘保国 一种有限复数信号测量系统与高精度分解方法
CN112485522B (zh) * 2020-12-09 2023-05-16 国网四川省电力公司电力科学研究院 基于电能数据感知的平顶窗函数同步相量测量方法及装置
CN112781723B (zh) * 2021-01-27 2023-09-12 南京微动智测信息技术有限公司 一种基于频谱方差的谐波成分检测方法
CN113189398A (zh) * 2021-04-29 2021-07-30 云南电网有限责任公司电力科学研究院 一种置零点频域加窗的高阶谐波分析方法及装置
CN113341226B (zh) * 2021-06-21 2022-04-29 合肥美的暖通设备有限公司 谐波检测方法、装置、变频器及存储介质
CN113468474B (zh) * 2021-09-06 2021-11-19 南京易司拓电力科技股份有限公司 基于根Mini-Norm的电网频率估计方法
CN115015633B (zh) * 2022-06-24 2024-08-09 南京航空航天大学 一种电力系统中的谐波与间谐波的频率估计方法
CN116735957B (zh) * 2023-06-07 2024-02-27 四川大学 计及主瓣重叠干扰的近频谐波与间谐波测量方法及系统

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH09127164A (ja) * 1995-10-27 1997-05-16 Hioki Ee Corp 高調波電流解析装置
CN101113995A (zh) * 2007-08-29 2008-01-30 湖南大学 基于Nuttall窗双峰插值FFT的基波与谐波检测方法
CN201352236Y (zh) * 2008-12-02 2009-11-25 湖南海兴电器有限责任公司 基于Kaiser窗双谱线插值FFT的谐波电能表
CN106771594A (zh) * 2016-12-08 2017-05-31 清华大学 一种电力系统的次/超同步谐波检测方法
CN108133704A (zh) * 2018-02-22 2018-06-08 成都启英泰伦科技有限公司 一种声源锁定系统
CN110174553A (zh) * 2019-06-27 2019-08-27 河北工业大学 一种基于解析模态分解的密集频率谐波/间谐波检测方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH09127164A (ja) * 1995-10-27 1997-05-16 Hioki Ee Corp 高調波電流解析装置
CN101113995A (zh) * 2007-08-29 2008-01-30 湖南大学 基于Nuttall窗双峰插值FFT的基波与谐波检测方法
CN201352236Y (zh) * 2008-12-02 2009-11-25 湖南海兴电器有限责任公司 基于Kaiser窗双谱线插值FFT的谐波电能表
CN106771594A (zh) * 2016-12-08 2017-05-31 清华大学 一种电力系统的次/超同步谐波检测方法
CN108133704A (zh) * 2018-02-22 2018-06-08 成都启英泰伦科技有限公司 一种声源锁定系统
CN110174553A (zh) * 2019-06-27 2019-08-27 河北工业大学 一种基于解析模态分解的密集频率谐波/间谐波检测方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
《一种新的谐波、间谐波参数估计算法》;林达斌;《电网技术》;20120630;第170-174页 *

Also Published As

Publication number Publication date
CN110837001A (zh) 2020-02-25

Similar Documents

Publication Publication Date Title
CN110837001B (zh) 一种电力系统中谐波和间谐波的分析方法与装置
Fan et al. Frequency estimator of sinusoid based on interpolation of three DFT spectral lines
CN104897960B (zh) 基于加窗四谱线插值fft的谐波快速分析方法及系统
Wen et al. Novel three-point interpolation DFT method for frequency measurement of sine-wave
CN106483374A (zh) 一种基于Nuttall双窗全相位FFT的谐波间谐波检测方法
CN101113995A (zh) 基于Nuttall窗双峰插值FFT的基波与谐波检测方法
CN106018956B (zh) 一种加窗谱线插值的电力系统频率计算方法
Borkowski et al. Comparison of sine-wave frequency estimation methods in respect of speed and accuracy for a few observed cycles distorted by noise and harmonics
Belega et al. Impact of harmonics on the interpolated DFT frequency estimator
CN111856401A (zh) 一种基于互谱相位拟合的时延估计方法
Ma et al. Harmonic and interharmonic analysis of mixed dense frequency signals
Fan et al. Frequency estimator of sinusoid by interpolated DFT method based on maximum sidelobe decay windows
Li et al. Frequency estimation based on symmetric discrete Fourier transform
CN105372492B (zh) 基于三条dft复数谱线的信号频率测量方法
Becejac et al. Impact of PMU data errors on modal extraction using matrix pencil method
CN112035790A (zh) 井间定位信号频率估计方法
CN114487589A (zh) 电网宽频信号自适应测量方法、装置及系统
Ruan et al. Improved Prony method for high-frequency-resolution harmonic and interharmonic analysis
CN109030942B (zh) 谐相角分析方法
CN114184838A (zh) 基于sn互卷积窗的电力系统谐波检测方法、系统及介质
CN113129912A (zh) 一种单音信号的检测方法
Liu et al. An approach to power system harmonic analysis based on triple-line interpolation discrete Fourier transform
Fan et al. A new frequency estimator of sinusoid based on interpolated FFT
CN116735957B (zh) 计及主瓣重叠干扰的近频谐波与间谐波测量方法及系统
Lv et al. Adaptive algorithm based on FFT for frequency estimation

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20211001