CN109557367A - 一种高频分辨率谐波和间谐波Prony方法及装置 - Google Patents
一种高频分辨率谐波和间谐波Prony方法及装置 Download PDFInfo
- Publication number
- CN109557367A CN109557367A CN201811236990.3A CN201811236990A CN109557367A CN 109557367 A CN109557367 A CN 109557367A CN 201811236990 A CN201811236990 A CN 201811236990A CN 109557367 A CN109557367 A CN 109557367A
- Authority
- CN
- China
- Prior art keywords
- prony
- coefficient
- prony coefficient
- power system
- matrix
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 158
- 239000011159 matrix material Substances 0.000 claims abstract description 69
- 239000013598 vector Substances 0.000 claims abstract description 56
- 238000013016 damping Methods 0.000 claims abstract description 21
- 238000005070 sampling Methods 0.000 claims abstract description 20
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 19
- 230000006870 function Effects 0.000 claims description 7
- 238000004891 communication Methods 0.000 claims description 6
- 238000012545 processing Methods 0.000 claims description 6
- 230000001052 transient effect Effects 0.000 claims description 4
- 238000004590 computer program Methods 0.000 claims description 2
- 238000004458 analytical method Methods 0.000 abstract description 22
- 238000010586 diagram Methods 0.000 description 27
- 238000005516 engineering process Methods 0.000 description 9
- 238000002474 experimental method Methods 0.000 description 6
- 238000007792 addition Methods 0.000 description 3
- 230000003321 amplification Effects 0.000 description 3
- 230000006378 damage Effects 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 238000003199 nucleic acid amplification method Methods 0.000 description 3
- 238000001228 spectrum Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 230000005611 electricity Effects 0.000 description 2
- 230000017105 transposition Effects 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000002411 adverse Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 239000003990 capacitor Substances 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000010438 heat treatment Methods 0.000 description 1
- 230000006698 induction Effects 0.000 description 1
- 239000004615 ingredient Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
- G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
- G01R23/16—Spectrum analysis; Fourier analysis
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
- G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
- G01R23/02—Arrangements for measuring frequency, e.g. pulse repetition rate; Arrangements for measuring period of current or voltage
-
- H—ELECTRICITY
- H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
- H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
- H02J3/00—Circuit arrangements for ac mains or ac distribution networks
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Engineering & Computer Science (AREA)
- Power Engineering (AREA)
- Complex Calculations (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明实施例提供一种高频分辨率谐波和间谐波Prony方法及装置，该方法包括：根据采样的电力系统信号，构造自相关矩阵，将所述自相关矩阵进行特征分解，得到对应的特征向量表示的所述自相关矩阵；基于子空间方法和第一Prony系数的对称性，根据所述特征向量和所述电力系统信号中的正弦波数量，得到第一Prony系数；根据所述第一Prony系数得到所述电力系统信号的每个分量的阻尼因子和频率。本发明实施例提供的一种高频分辨率谐波和间谐波Prony方法及装置，利用Prony系数的对称性和子空间方法的优点，大大降低了计算负担，提高了电力系统高频分辨率谐波和间谐波分析的精度。
Description
技术领域
本发明涉及电力电子领域，尤其涉及一种高频分辨率谐波和间谐 波Prony方法及装置。
背景技术
电力电子系统和工业中时变非线性负载的使用越来越多，导致了 严重的谐波和间谐波失真。电力系统中的谐波和间谐波分量可能导致 额外的功率损耗，设备发热和损坏。它们还会干扰通信电路，导致电 网中的谐振以及保护和控制设备的异常操作。准确分析谐波和间谐波 对于有效防止不良影响和全面了解电能质量问题至关重要。
目前，已经提出了许多技术来分析电力系统中的谐波和间谐波。 在这些分析方法中，快速傅里叶变换(FFT)是一种强大的工具。计算 效率高和简单的求解过程是FFT的优点。然而，基于周期性采样提出了 DFT，并且FFT技术直接应用于频谱分析将导致间谐波和基频偏差出现 不准确的情况。加窗插值DFT(WIDFT)算法可以显着减少由DFT的 频谱泄漏和栅栏引起的不准确性，并且不会显着增加计算负担。但是， DFT和WIDFT均由固定的频率分辨率来表征，这些分辨率由数据窗口 长度决定。如果存在接近频率的谐波和间谐波，则DFT和WIDFT不能 分离这些成分，并且可能给谐波/间谐波分析提供高度不准确的结果。
近年来，一些先进的用于电力系统谐波/间谐波分析的高频分辨率 方法已经开发出来。Prony方法是一种不以固定频率分辨率为特征的方 法。然而，Prony系数的计算通常涉及高计算负担，并且采样信号需要 预滤波或其他技术来提供平滑信号并消除噪声。
另一种强大的高频分辨率方法是子空间方法。子空间方法通过已 知协方差函数的背景噪声中的随机正弦信号之和来处理模拟的随机信 号。协方差矩阵的特征向量被分为两个正交组：跨越信号空间的特征 向量和跨越噪声空间的特征向量。最重要的技术之一是min-norm方法。 Min-norm方法使用一个向量进行频率估计。这个矢量属于噪声子空间，具有最小欧几里德范数，其第一个元素等于1。另一个基于子空间理论 的重要技术是通过旋转不变技术估计信号参数，即ESPRIT方法。 ESPRIT使用信号子空间的旋转不变性来估计信号的参数。子空间方法 的特征不在于固定的频率分辨率。然而，子空间方法需要构造和分解 通常涉及高计算负担的自相关矩阵，并且需要事先估计估计次数。
发明内容
本发明实施例为克服上述技术缺陷，提供一种高频分辨率谐波和 间谐波Prony方法及装置。
第一方面，本发明实施例提供一种高频分辨率谐波和间谐波Prony 方法，包括：
根据采样的电力系统信号，构造自相关矩阵，将所述自相关矩阵 进行特征分解，得到对应的特征向量表示的所述自相关矩阵；
基于子空间方法和第一Prony系数的对称性，根据所述特征向量 和所述电力系统信号中的正弦波数量，得到第一Prony系数；
根据所述第一Prony系数得到所述电力系统信号的每个分量的阻 尼因子和频率。
第二方面，本发明实施例提供一种高频分辨率谐波和间谐波Prony 装置，包括：
特征分解模块，用于根据采样的电力系统信号，构造自相关矩阵， 将所述自相关矩阵进行特征分解，得到对应的特征向量表示的所述自 相关矩阵；
Prony系数求解模块，用于基于子空间方法和第一Prony系数的对 称性，根据所述特征向量和所述电力系统信号中的正弦波数量，得到 第一Prony系数；
处理模块，用于根据所述第一Prony系数得到所述电力系统信号 的每个分量的阻尼因子和频率。
第三方面，本发明实施例提供一种电子设备，包括存储器和处 理器，所述处理器和所述存储器通过总线完成相互间的通信；所述 存储器存储有可被所述处理器执行的程序指令，所述处理器调用所 述程序指令能够执行如第一方面所述的方法。
第四方面，本发明实施例提供一种非暂态计算机可读存储介质， 其上存储有计算机程序，该计算机程序被处理器执行时实现如第一方 面所述高频分辨率谐波和间谐波Prony方法。
本发明实施例提供的一种高频分辨率谐波和间谐波Prony方法及 装置，利用Prony系数的对称性和子空间方法的优点，大大降低了计 算负担，提高了电力系统高频分辨率谐波和间谐波分析的精度。
附图说明
图1为本发明实施例提供的一种高频分辨率谐波和间谐波Prony 方法的流程示意图；
图2为本发明实施例提供的示例一的FPE示意图；
图3为本发明实施例提供的示例二的FPE示意图；
图4为本发明实施例提供的示例三的FPE示意图；
图5为本发明实施例提供的示例三的FPE放大视图；
图6为本发明实施例提供的示例四的FPE示意图；
图7为本发明实施例提供的示例四的FPE放大视图；
图8为本发明实施例提供的示例五的FPE示意图；
图9为本发明实施例提供的示例五的FPE放大视图；
图10为本发明实施例提供的示例六的FPE示意图；
图11为本发明实施例提供的实验室单相整流器/逆变器电路的结 构图；
图12为本发明实施例提供的电路整流模式操作期间的采样电流 图；
图13为本发明实施例提供的正弦分量的数量设定为M＝25时使用 改进的Prony方法的分析结果图；
图14为本发明实施例提供的使用图13所示的分析结果的重建波 形图；
图15为本发明实施例提供的一种高频分辨率谐波和间谐波Prony 装置的结构示意图；
图16为本发明实施例提供的一种电子设备的实体结构示意图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合本发 明实施例中的附图，对本发明中的技术方案进行清楚、完整地描述， 显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施 例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性 劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
图1为本发明实施例提供的一种高频分辨率谐波和间谐波Prony 方法的流程示意图，如图1所示，该方法包括：
步骤11，根据采样的电力系统信号，构造自相关矩阵，将所述自 相关矩阵进行特征分解，得到对应的特征向量表示的所述自相关矩阵；
步骤12，基于子空间方法和第一Prony系数的对称性，根据所述 特征向量和所述电力系统信号中的正弦波数量，得到第一Prony系数；
步骤13，根据所述第一Prony系数得到所述电力系统信号的每个 分量的阻尼因子和频率。
考虑采样的电力系统信号x(n)，x(n)由正弦分量和添加的高斯白 噪声(AWGN)w(n)组成：
其中Ts是采样时间，w(n)是增加的噪声，M是正弦波的数量。Ak， ak，ωk＝2πfk和是第k个分量的幅度，阻尼因子，角速度和初始相 位。Prony方法包括首先求解下面的线性方程组，找出阻尼因子和频率：
其中n＝2M,2M+1,…,N-1。方程(2)构成2M未知数的线性方 程系统，即a(m)系数。
一旦a(m)系数已知，每个分量的阻尼因子和频率可以通过以下函 数F(z)的根来计算：
可以通过求解将这些未知与采样数据相关联的第二组线性方程来 计算每个指数的幅度和相位角，如下所述：
ZA＝X (4)
其中X＝[x(0),x(1),…,x(N-1)]H,
zi(i＝1,2,...,2M)是(3)中函数F(z)的根， 符号H表示共轭转置操作。
另一种强大的谐波和间谐波分析方法是子空间方法。子空间方法 将样本投影到信号和噪声子空间中。考虑如(1)所示的采样电力系统信 号x(n)。L×L自相关矩阵是一个L×L的厄米特矩阵，Rx＝Rx H,Rx可 以用它的本征分解来表示，如下
其中Λ＝diag(λ1,λ2,…λL)包含Rx的特征值的降序， V＝[v1,v2,…vL]是相应的特征向量的矩阵。特征向量可以分为两个正 交组：信号特征向量的L×2M矩阵和噪声特征向量的L×(L-2M)矩阵 Vn＝[v2M+1,v2M+2,…vL]。显然，Vs与Vn正交。
Vs⊥Vn (6)
基于分离信号和噪声子空间中的数据的重要技术是最小范数方 法。min-norm方法使用一个向量进行频率估计。该向量属于噪声子空 间，具有最小欧几里德范数，其第一个元素等于1。另一项重要技术是 通过旋转不变技术(ESPRIT)估计信号参数。ESPRIT方法使用信号 子空间的旋转不变性特征。
本发明实施例提供的一种高频分辨率谐波和间谐波Prony方法， 结合子空间方法，在现有的Prony方法上进行改进，首先根据采样的 电力系统信号x(n)，构造(2M+1)×(2M+1)自相关矩阵Rx，将所述自 相关矩阵进行特征分解，用对应的特征值和特征向量表示所述自相关 矩阵。
所述得到对应的特征向量表示的所述自相关矩阵，具体包括：
其中，Rx为(2M+1)×(2M+1)的所述自相关矩阵， V＝[v1,v2,…,v2M+1]为所述特征向量的矩阵，Λ＝diag(λ1,λ2,…,λ2M+1)是 由Rx的特征值组成的向量，λ1,λ2,…,λ2M+1为Rx的特征值的降序排列， M为所述正弦波数量。
其次，基于子空间方法和第一Prony系数a(m)的对称性，根据特 征向量V＝[v1,v2,…,v2M+1]和采样的电力系统信号中的正弦波数量M， 得到第一Prony系数a(m)，所述第一Prony系数a(m)用于求解所述电 力系统信号中的每个分量的阻尼因子和频率。根据第一Prony系数a(m) 的对称性，能够将(2M+1)×1的特征向量降低为(M+1)×1，降低了计 算量。
最后，根据所述第一Prony系数a(m)得到所述每个分量的阻尼因 子和频率，从而实现高频分辨率谐波和间谐波的分析。
间谐波的特点是放大电压闪变和音频干扰，造成感应电动机振动 及异常。对于由电容、电感和电阻构成的无源滤波器电路，间谐波可 能会被放大，严重时会使滤波器因谐波过载而不能正常运行，甚至造 成损坏。间谐波的影响和危害等同整数次谐波电压的影响和危害。
本发明实施例提供的一种高频分辨率谐波和间谐波Prony方法， 利用Prony系数的对称性和子空间方法的优点，大大降低了计算负担， 提高了电力系统高频分辨率谐波和间谐波分析的精度，从而减小了间 谐波对于电路的影响和危害。
在上述实施例的基础上，所述根据所述特征向量和所述电力系统 信号中的正弦波数量，得到第一Prony系数，具体包括：
基于子空间方法和所述特征向量，表示所述第一Prony系数；
基于所述第一Prony系数的对称性，根据第一Prony系数表示第二 Prony系数；
根据所述特征向量和所述正弦波数量，求解得到所述第二Prony 系数；
根据所述第二Prony系数得到所述第一Prony系数。
首先，基于子空间方法和特征向量，表示第一Prony系数，具体 包括：
根据x(n)＝s(n)+w(n)将所述特征向量的矩阵V＝[v1,v2,…,v2M+1] 分为包含信号特征向量的(2M+1)×(2M)矩阵Vs和噪声特征向量的 (2M+1)×1矩阵Vn，其中，x(n)为所述电力系统信号，由正弦分量s(n) 和高斯白噪声w(n)组成，Vs＝[v1,v2,…,v2M]，Vn＝[v2M+1]；
根据Vs与Vn正交，得到其中 n＝2M+1,2M+2,...,N-1；
结合线性方程组和得 到所述第一Prony系数的表达式a(m)＝v2M+1(2M-m)，其中，a(m)为 所述第一Prony系数。
考虑采样的电力系统信号x(n)，x(n)由正弦分量和添加的高斯白噪 声(AWGN)w(n)组成：
其中Ts是采样时间，w(n)是增加的噪声，M是正弦波的数量。Ak， ak，ωk＝2πfk和是第k个分量的幅度，阻尼因子，角速度和初始相 位。构造(2M+1)×(2M+1)自相关矩阵Rx。并且Rx可以用其特征分解 表示如下：
其中Λ＝diag(λ1,λ2,…λ2M+1)包含降序Rx的特征值， V＝[v1,v2,…v2M+1]是相应特征向量的矩阵。特征向量可以分为两个正交 组：信号特征向量的(2M+1)×(2M)矩阵Vs＝[v1,v2,…v2M]和噪声特征 向量的(2M+1)×1矩阵Vn＝[v2M+1]。注意电力系统信号x(n)中的信号 s(n)必须在信号子空间中，并且特征向量v2M+1必须在噪声子空间中。 然后，可以基于Vs⊥Vn获得以下关系。
其中n＝2M+1,2M+2,...,N-1。
等式(8)示出了特征向量v2M+1可用于表示(2)中的a(m)如下：
a(m)＝v2M+1(2M-m) (9)
其次，基于所述第一Prony系数a(m)的对称性，根据第一Prony 系数a(m)表示第二Prony系数b(m)，具体包括：
根据第一Prony系数的对称性a(m)，即：
a(m)＝a(2M-m)，
得到第二Prony系数b(m)的表达式：
其中，b(m)为所述第二Prony系数，M为正弦波数量。
最后，根据特征向量和正弦波数量M，求解得到第二Prony系数 b(m)，具体包括：
根据线性方程组和第一Prony系数a(m)的对称 性a(m)＝a(2M-m)，得到：
结合信号特征向量和噪声特征向量的正交性Vs⊥Vn，根据得到：
b(m)＝vM+1(m)。
第一Prony系数a(m)具有对称性，即：
a(m)＝a(2M-m)，m＝0,1,...,M，
因此，(2)式可以改写为：
这里
方程(10)是M未知的线性方程组，即b(m)系数。为了获得(10)中的 b(m)，采样的信号被转换为
y(n)＝x(n)+x(2M-n)
＝[s(n)+s(2M-n)]+[w(n)+w(2M-n)](11)
＝s'(n)+w'(n)
y(n)的(M+1)×(M+1)自相关矩阵Ry可以用它的本征分解来表 示如下：
其中Λ＝diag(λ1,λ2,…λM+1)包含Ry的特征值的降序， V＝[v1,v2,…vM+1]是相应特征向量的矩阵。
特征向量可以分为两个正交分组：信号特征向量的(M+1)×M矩 阵Vs＝[v1,v2,…vM]和噪声特征向量的(M+1)×1矩阵Vn＝[vM+1]。
(8)式中的信号s'(n)＝s(n)+s(2M-n)必 须在信号子空间中，并且特征向量vM+1必须在噪声子空间中。然后， 可以基于(6)式Vs⊥Vn获得以下关系：
等式(2)，(10)和(13)表明，特征向量vM+1可用于表示(10)中的b(m)如 下：
b(m)＝vM+1(m) (14)
组合(10)和(14)，然后使用如下简单关系获得a(m)：
获取到a(m)系数后，通过函数的根计算电力系 统信号中的每个分量的阻尼因子和频率，然后根据下式计算每个分量 的幅度和相位角：
ZA＝X。
本发明实施例提供的一种高频分辨率谐波和间谐波Prony方法， 利用Prony系数的对称性和子空间方法的优点，大大降低了计算负担， 提高了电力系统高频分辨率谐波和间谐波分析的精度。
在上述实施例的基础上，所述根据所述第一Prony系数得到所述 电力系统信号的每个分量的阻尼因子和频率，具体包括：
根据函数F(z)的根计算得到所述电力系统信号的每个分量的阻尼 因子和频率，其中，式中a(m)为所述第一Prony 系数，M为所述正弦波数量。
求出第一Prony系数a(m)系数后，每个分量的阻尼因子和频率可 以通过以下函数F(z)的根来计算：
可以通过求解将这些未知与采样数据相关联的第二组线性方程来 计算每个指数的幅度和相位角，如下所述。
ZA＝X (4)
其中X＝[x(0),x(1),…,x(N-1)]H,
zi(i＝1,2,...,2M)是(3)中函数F(z)的根， 符号H表示共轭转置操作。
本发明实施例提供的一种高频分辨率谐波和间谐波Prony方法， 利用Prony系数的对称性和子空间方法的优点，大大降低了计算负担， 提高了电力系统高频分辨率谐波和间谐波分析的精度。
在上述实施例的基础上，当处理Prony方法和子空间方法的应用 时，最重要的问题之一是找到最合适的正弦曲线(指数)或估计次数 的数量，以符合结果准确性和计算工作。如果次数设置得太低，会获 得高度平滑的光谱，并且可能无法将近似频率的元件分开；另一方面， 如果次数选择得太高，就有可能在频谱中引入虚假的低电平峰值。
已经提出了一些标准来解决这个问题，最终预测误差(FPE)就是 其中之一。使用平方预测误差之和的期望来定义FPE。在本发明实施 例中，FPE用于确定正弦曲线的数量。FPE定义如下：
其中第一个因素是所考虑的模型保真度的度量，并且随着M的 增加趋于减小。第二个项随M增加并构成一个惩罚项，旨 在防止使用太多的参数。
考虑到正弦波数量的确定，进行了一些实验来分析传统Prony方 法和提出的Prony方法之间的不同特征。在本发明实施例的实验中， 采样频率设置为3200Hz，标称基频设置为50Hz。对于改进的Prony 和传统的Prony，信号矩阵的大小分别设置为(M+1)×64，2M×64。因此， 只要1/30s的数据窗口足以获得两种方法的测量结果。
由于基于奈奎斯特采样定理，谐波的最大阶数为31，因此正弦波 的数量被设置在[1，31]的范围内。MATLAB用于软件实现。为了更好 地理解与应用本发明提出的一种高频分辨率谐波和间谐波Prony方法， 本发明进行以下示例，且本发明不仅局限于以下示例。
示例一：样本信号的模型如下：x(t)＝sin(2π50.1t)。采样信号中只 有一个正弦波，基频的偏差为0.1Hz。
首先，将估计次数正确地设置为M＝1，传统的Prony方法和本发 明实施例提供的改进的Prony方法都可以给出准确的分析。然而，如 果没有正确设置估计次数，例如，M＝2,3,...,31，传统的Prony方法将 给出大误差的分析，而本发明实施例提供的改进的Prony方法仍然可 以在大多数时间给出准确的分析。
图2为本发明实施例提供的示例一的FPE示意图，其中图2中的 (a)图为本发明实施例提供的改进的Prony方法的FPE示意图，M＝30， 图2中的(b)图为传统的Prony方法的FPE示意图，M＝21，如图2所 示，表明：对于传统的Prony方法，除M＝1外，FPE的值很大。然而， 对于改进的Prony方法，除了M＝16,23,28之外，FPE的值小于1.5e-9。 这表明传统的Prony只能给出M＝1的准确估计，但是，改进的Prony 方法可以给出大范围的M的准确估计，实验结果与这个结论一致。
进一步分析表明，对于M＝16,23,28，仍然可以使用改进的Prony 方法准确地估计频率。然而，在使用(4)计算幅度时出现奇异矩阵， 并且不能给出准确的幅度估计。
示例二：样本信号的模型如下：
这里基频频率偏差为0.1Hz，信号中有20％的二次，三次和四次 谐波。
采样信号中有四个正弦曲线，实验表明，即使将估计次数正确设 置为M＝4，传统的Prony方法也不能给出准确的估计。然而，当M≥4 时，改进的Prony方法可以给出所有分量的准确估计。
图3为本发明实施例提供的示例二的FPE示意图，其中图3中的 (a)图为本发明实施例提供的改进的Prony方法的FPE示意图，M＝30， 图3中的(b)图为传统的Prony方法的FPE示意图，M＝21，如图3所 示，表明：对于传统的Prony方法，FPE的值很大，最小值为4.6145。 然而，对于改进的Prony方法，对于M＞4，FPE的值小于7e-8。这表 明，无论M是什么值，传统的Prony都无法给出准确的估计。然而， 改进的Prony可以在很宽的M范围内给出准确的估计，实验结果与该 结论一致。
示例三：采样信号的模型如下：
x(t)＝[1+0.1sin(2π1t)]sin(2π50t)，
其中信号中存在1Hz幅度调制。
采样信号可表示为：
x(t)＝sin(2π50t)-0.05cos(2π49t)+0.05cos(2π51t)，
在采样信号中有三个具有接近频率的正弦分量为49Hz，50Hz和 51Hz，用短数据窗口将这三个组件分开为40ms是很困难的。
计算FPE的值。对于改进的Prony方法，FPE的最小值在M＝30 时获得，FPE＝1.4e-17。对于传统的Prony方法，FPE的最小值在M＝1 时获得，FPE＝1.3e-4。
图4为本发明实施例提供的示例三的FPE示意图，其中图4中的 (a)图为本发明实施例提供的改进的Prony方法的FPE示意图，图4中 的(b)图为传统的Prony方法的FPE示意图，如图4所示，其中对于改 进的Prony，M＝30，对于传统的Prony方法，M＝1。图5为本发明实 施例提供的示例三的FPE放大视图，其中图5中的(a)图为图4中的(a) 图对应的FPE放大视图，图5中的(b)图为图4中的(b)图对应的FPE放 大视图。图4-5表明：改进的Prony方法可以准确地分离三个正弦分量， 然而，传统的Prony方法不能分离三个正弦分量。频率为49Hz和51Hz 的分量与基波分量合并，使估计的基波幅度高于实际值。
示例四：采样信号的模型如下：
测试信号包含基波分量和第2，第3，第4，第5，第6，第7谐波， 以及频率为135Hz的间谐波，w(t)是白噪声，信噪比(SNR)是50dB。
对于改进的Prony方法，FPE的最小值在M＝31时获得，FPE＝0.013。 对于传统的Prony方法，FPE的最小值在M＝24时获得，FPE＝7.380。
图6为本发明实施例提供的示例四的FPE示意图，其中图6中的(a) 图为本发明实施例提供的改进的Prony方法的FPE示意图，图6中的(b) 图为传统的Prony方法的FPE示意图，如图6所示，其中对于改进的 Prony，M＝31，对于传统的Prony方法，M＝24。图7为本发明实施例 提供的示例四的FPE放大视图，其中图7中的(a)图为图6中的(a)图对应 的FPE放大视图，图7中的(b)图为图6中的(b)图对应的FPE放大视图，图 6-7表明：改进的Prony方法可以在嘈杂的环境中准确地分离出8个正弦 分量。然而，传统的Prony方法无法给出准确的估计，频率为135Hz和 150Hz的分量合并在一起，二次谐波的估计幅度远高于实际值。
将所提出的改进的Prony方法与一些最近提出的高频分辨率谐波/ 间谐波分析方法进行比较，包括：传统的Prony方法，Min-norm方法 和Esprit方法。在我们的实验中，使用上述四种方法找到正弦分量的 频率，然后通过求解如(4)所示的线性方程组来计算每个分量的幅度 和相位角。MATLAB用于软件实现，采样频率设置为3200Hz。
为了便于比较，对于改进的Prony方法，传统的Prony方法， Min-norm方法和ESPRIT方法，信号矩阵的大小设置为(M+1)×64， 2M×64，64×64，64×64。因此，长达1/30s的数据窗口足以获得所有上 述方法的测量结果。
示例五：采样信号的模型如下：
其中，f1＝50Hz，采样信号包含20％的二次，三次，四次谐波和 一个频率为187Hz的间谐波，w(t)为白噪声，信噪比为60dB。
对于改进的Prony方法，传统的Prony方法，Min-norm方法和 ESPRIT方法，分别在M＝31，M＝28，M＝5，M＝30时获得FPE的 最小值。因此，相应地设置四种方法的估计次数。
图8为本发明实施例提供的示例五的FPE示意图，其中图8中的 (a)图为本发明实施例提供的改进的Prony方法的FPE示意图，M＝31， (b)图为传统的Prony方法的FPE示意图，M＝28，(c)图为Min-norm方 法的FPE示意图，M＝5，(d)图为ESPRIT方法的FPE示意图，M＝30， 图9为本发明实施例提供的示例五的FPE放大视图，其中图9中的(a) 图为图8中的(a)图对应的FPE放大视图，图9中的(b)图为图8中的(b) 图对应的FPE放大视图。图8-9显示改进的Prony方法可以精确估计 基波分量和所有谐波/间谐波分量，它提供了所有上述方法中的最佳估 计。
在本实验中，对于改进的Prony方法，Min-norm方法和ESPRIT 方法，自相关矩阵的大小分别设置为31×31,64×64,64×64。需要自相关 矩阵的构造和特征分解，其计算复杂度分别为O(313)，O(643)和O(643)。 传统的Prony方法不需要分解自相关矩阵。然而，它需要2M×2M矩 阵的逆操作，其计算复杂度为O(563)。因此，与其他提到的方法相比， 改进的Prony方法的计算复杂度已经大大降低。近频信号在嘈杂的环 境中通常难以区分具有近频率的信号，该实验用于比较在嘈杂环境中 区分具有近频率的分量的方法的能力。
示例六：采样信号的模型如下：
x(t)＝sin(2πf1t)+0.2sin(2πfit)+w(t)，
其中f1＝50Hz，fi＝55Hz，w(t)是白噪声，SNR＝60dB。
对于改进的Prony，传统的Prony，Min-norm和ESPRIT方法，在 M＝30，M＝21，M＝2，M＝31时获得FPE的最小值，因此，相应地 设置四种方法的估计次数。
图10为本发明实施例提供的示例六的FPE示意图，其中图10中 的(a)图为本发明实施例提供的改进的Prony方法的FPE示意图， M＝30，(b)图为传统的Prony方法的FPE示意图，M＝21，(c)图为 Min-norm方法的FPE示意图，M＝2，(d)图为ESPRIT方法的FPE示 意图，M＝31，如图10所示，图10显示改进的Prony方法在所有提到 的方法中给出了最佳估计。
图11为本发明实施例提供的实验室单相整流器/逆变器电路的结 构图，如图11所示，直流侧是48V直流电池，交流侧连接到交流负载 和电网。该电路有两种工作条件：逆变器模式和整流器模式。在逆变 器模式下，电池中的能量转换为电网。在整流器模式中，电网能量转 换为DC侧，或者为DC电池充电，在图11中表示为iinv的电流是采 样电流。图12为本发明实施例提供的电路整流模式操作期间的采样电 流图，如图12所示，采样频率为3000Hz。
计算在[1，29]范围内的M的FPE值，并且对于改进的Prony方法， 在M＝25时获得FPE的最小值。图13为本发明实施例提供的正弦分量 的数量设定为M＝25时使用改进的Prony方法的分析结果图，图14为 本发明实施例提供的使用图13所示的分析结果的重建波形图，图14 表示重建波形精确地接近原始采样波形，这表明图13的分析结果是高 度准确的。
本发明实施例提供的一种高频分辨率谐波和间谐波Prony方法， 利用Prony系数的对称性和子空间方法的优点，大大降低了计算负担， 提高了电力系统高频分辨率谐波和间谐波分析的精度。
图15为本发明实施例提供的一种高频分辨率谐波和间谐波Prony 装置的结构示意图，如图15所示，包括特征分解模块151、Prony系 数求解模块152和处理模块153，其中：
特征分解模块151，用于根据采样的电力系统信号，构造自相关矩 阵，将所述自相关矩阵进行特征分解，得到对应的特征向量表示的所 述自相关矩阵；
Prony系数求解模块152，用于基于子空间方法和第一Prony系数 的对称性，根据所述特征向量和所述电力系统信号中的正弦波数量， 得到第一Prony系数；
处理模块153，用于根据所述第一Prony系数得到所述电力系统信 号的每个分量的阻尼因子和频率。
本发明实施例提供的一种高频分辨率谐波和间谐波Prony装置， 结合子空间方法，在现有的Prony方法上进行改进，首先特征分解模 块151根据采样的电力系统信号x(n)，构造(2M+1)×(2M+1)自相关 矩阵Rx，将所述自相关矩阵进行特征分解，用对应的特征值和特征向 量表示所述自相关矩阵。
所述用对应的特征值和特征向量表示所述自相关矩阵，具体包括：
其中，Rx为(2M+1)×(2M+1)的所述自相关矩阵， V＝[v1,v2,…,v2M+1]为所述特征向量的矩阵，Λ＝diag(λ1,λ2,…,λ2M+1)包 含降序Rx的特征值，M为所述正弦波数量。
其次，Prony系数求解模块152基于子空间方法和第一Prony系数 a(m)的对称性，根据特征向量V＝[v1,v2,…,v2M+1]和采样的电力系统信 号中的正弦波数量M，得到第一Prony系数a(m)，所述第一Prony系 数a(m)用于求解所述电力系统信号中的每个分量的阻尼因子和频率。 根据第一Prony系数a(m)的对称性，能够将(2M+1)×1的特征向量降 低为(M+1)×1，降低了计算量。
最后，处理模块153根据所述第一Prony系数a(m)得到电力系统 信号的每个分量的阻尼因子和频率，从而实现高频分辨率谐波和间谐 波的分析。
本发明实施例提供的装置是用于执行上述方法实施例的，具体流 程请参加上述方法实施例，此处不再赘述。
本发明实施例提供的一种高频分辨率谐波和间谐波Prony装置， 利用Prony系数的对称性和子空间方法的优点，大大降低了计算负担， 提高了电力系统高频分辨率谐波和间谐波分析的精度。
图16为本发明实施例提供的一种电子设备的实体结构示意图，如 图16所示，该电子设备可以包括：处理器(processor)161、通信接口 (Communications Interface)162、存储器(memory)163和总线164，其中， 处理器161，通信接口162，存储器163通过总线164完成相互间的通 信。总线164可以用于电子设备与传感器之间的信息传输。处理器161 可以调用存储器163中的逻辑指令，以执行如下方法：根据采样的电 力系统信号，构造自相关矩阵，将所述自相关矩阵进行特征分解，得 到对应的特征向量表示的所述自相关矩阵；基于子空间方法和第一 Prony系数的对称性，根据所述特征向量和所述电力系统信号中的正弦波数量，得到第一Prony系数；根据所述第一Prony系数得到所述电力 系统信号的每个分量的阻尼因子和频率。
此外，上述的存储器163中的逻辑指令可以通过软件功能单元的 形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可 读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说 对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的 形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干 指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网 络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前 述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟 或者光盘等各种可以存储程序代码的介质。
本发明实施例提供一种非暂态计算机可读存储介质，该非暂态计 算机可读存储介质存储计算机指令，该计算机指令使计算机执行上述 实施例所提供的一种高频分辨率谐波和间谐波Prony方法，例如包括： 根据采样的电力系统信号，构造自相关矩阵，将所述自相关矩阵进行 特征分解，得到对应的特征向量表示的所述自相关矩阵；基于子空间 方法和第一Prony系数的对称性，根据所述特征向量和所述电力系统 信号中的正弦波数量，得到第一Prony系数；根据所述第一Prony系数 得到所述电力系统信号的每个分量的阻尼因子和频率。
以上所述仅为本发明的优选实施例，并不用于限制本发明。本发 明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的 修改或补充，但并不会偏离本发明的精神或者超越所附权利要求书定 义的范围。
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而 非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领 域的技术人员应当理解：其依然可以对前述各实施例所记载的技术方 案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或 者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的 精神和范围。
Claims (10)
1.一种高频分辨率谐波和间谐波Prony方法，其特征在于，包括：
根据采样的电力系统信号，构造自相关矩阵，将所述自相关矩阵进行特征分解，得到对应的特征向量表示的所述自相关矩阵；
基于子空间方法和第一Prony系数的对称性，根据所述特征向量和所述电力系统信号中的正弦波数量，得到第一Prony系数；
根据所述第一Prony系数得到所述电力系统信号的每个分量的阻尼因子和频率。
2.根据权利要求1所述的方法，其特征在于，所述得到对应的特征向量表示的所述自相关矩阵，具体包括：
其中，Rx为(2M+1)×(2M+1)的所述自相关矩阵，V＝[v1,v2,…,v2M+1]为所述特征向量的矩阵，Λ＝diag(λ1,λ2,…,λ2M+1)是由Rx的特征值组成的向量，λ1,λ2,…,λ2M+1为Rx的特征值的降序排列，M为所述正弦波数量。
3.根据权利要求2所述的方法，其特征在于，所述根据所述特征向量和所述电力系统信号中的正弦波数量，得到第一Prony系数，具体包括：
基于子空间方法和所述特征向量，表示所述第一Prony系数；
基于所述第一Prony系数的对称性，根据第一Prony系数表示第二Prony系数；
根据所述特征向量和所述正弦波数量，求解得到所述第二Prony系数；
根据所述第二Prony系数得到所述第一Prony系数。
4.根据权利要求3所述的方法，其特征在于，所述基于子空间方法和所述特征向量，表示所述第一Prony系数，具体包括：
根据x(n)＝s(n)+w(n)将所述特征向量的矩阵V＝[v1,v2,…,v2M+1]分为包含信号特征向量的(2M+1)×(2M)矩阵Vs和噪声特征向量的(2M+1)×1矩阵Vn，其中，x(n)为所述电力系统信号，由正弦分量s(n)和高斯白噪声w(n)组成，Vs＝[v1,v2,…,v2M]，Vn＝[v2M+1]；
根据Vs与Vn正交，得到其中n＝2M+1,2M+2,...,N-1；
结合线性方程组和得到所述第一Prony系数的表达式a(m)＝v2M+1(2M-m)，其中，a(m)为所述第一Prony系数。
5.根据权利要求4所述的方法，其特征在于，所述基于所述第一Prony系数的对称性，根据第一Prony系数表示第二Prony系数，具体包括：
根据所述第一Prony系数的对称性，a(m)＝a(2M-m)，得到所述第二Prony系数的表达式其中，b(m)为所述第二Prony系数。
6.根据权利要求5所述的方法，其特征在于，所述根据所述特征向量和所述正弦波数量，求解得到所述第二Prony系数，具体包括：
根据所述线性方程组和所述第一Prony系数的对称性得到
结合信号特征向量和噪声特征向量的正交性，根据所述得到b(m)＝vM+1(m)。
7.根据权利要求1-6任一项所述的方法，其特征在于，所述根据所述第一Prony系数得到所述电力系统信号的每个分量的阻尼因子和频率，具体包括：
根据函数F(z)的根计算得到所述电力系统信号的每个分量的阻尼因子和频率，其中，式中a(m)为所述第一Prony系数，M为所述正弦波数量。
8.一种高频分辨率谐波和间谐波Prony装置，其特征在于，包括：
特征分解模块，用于根据采样的电力系统信号，构造自相关矩阵，将所述自相关矩阵进行特征分解，得到对应的特征向量表示的所述自相关矩阵；
Prony系数求解模块，用于基于子空间方法和第一Prony系数的对称性，根据所述特征向量和所述电力系统信号中的正弦波数量，得到第一Prony系数；
处理模块，用于根据所述第一Prony系数得到所述电力系统信号的每个分量的阻尼因子和频率。
9.一种电子设备，其特征在于，包括存储器和处理器，所述处理器和所述存储器通过总线完成相互间的通信；所述存储器存储有可被所述处理器执行的程序指令，所述处理器调用所述程序指令能够执行如权利要求1至7任一所述的方法。
10.一种非暂态计算机可读存储介质，其上存储有计算机程序，其特征在于，该计算机程序被处理器执行时实现如权利要求1至7任一项所述高频分辨率谐波和间谐波Prony方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811236990.3A CN109557367B (zh) | 2018-10-23 | 2018-10-23 | 一种高频分辨率谐波和间谐波Prony方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811236990.3A CN109557367B (zh) | 2018-10-23 | 2018-10-23 | 一种高频分辨率谐波和间谐波Prony方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109557367A true CN109557367A (zh) | 2019-04-02 |
| CN109557367B CN109557367B (zh) | 2020-09-08 |
Family
ID=65865074
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811236990.3A Expired - Fee Related CN109557367B (zh) | 2018-10-23 | 2018-10-23 | 一种高频分辨率谐波和间谐波Prony方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109557367B (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110830133A (zh) * | 2019-12-23 | 2020-02-21 | 华中科技大学 | 一种基于多阶信道预测方法、预测系统及应用 |
| CN113032716A (zh) * | 2019-12-24 | 2021-06-25 | 南京理工大学 | 基于加窗插值和Prony算法的谐波与间谐波分析方法 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101368987A (zh) * | 2008-09-27 | 2009-02-18 | 武汉大学 | 一种电力系统稳态谐波和/或间谐波测量方法 |
| CN101566649A (zh) * | 2009-05-27 | 2009-10-28 | 重庆大学 | 一种电力系统间谐波检测方法 |
| CN103630742A (zh) * | 2013-12-16 | 2014-03-12 | 国家电网公司 | 一种动态信号参数的获取方法 |
| CN104181389A (zh) * | 2014-07-02 | 2014-12-03 | 中国农业大学 | 电力系统中相量测量方法 |
| CN105606893A (zh) * | 2016-01-26 | 2016-05-25 | 江苏科技大学 | 基于空间平滑修正music的电力间谐波检测方法 |
| WO2017043143A1 (ja) * | 2015-09-10 | 2017-03-16 | 株式会社日立製作所 | 分散電源発電量推定装置および方法 |
| CN107545253A (zh) * | 2017-09-11 | 2018-01-05 | 国网江西省电力公司电力科学研究院 | 一种电力系统低频振荡监测方法 |
-
2018
- 2018-10-23 CN CN201811236990.3A patent/CN109557367B/zh not_active Expired - Fee Related
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101368987A (zh) * | 2008-09-27 | 2009-02-18 | 武汉大学 | 一种电力系统稳态谐波和/或间谐波测量方法 |
| CN101566649A (zh) * | 2009-05-27 | 2009-10-28 | 重庆大学 | 一种电力系统间谐波检测方法 |
| CN103630742A (zh) * | 2013-12-16 | 2014-03-12 | 国家电网公司 | 一种动态信号参数的获取方法 |
| CN104181389A (zh) * | 2014-07-02 | 2014-12-03 | 中国农业大学 | 电力系统中相量测量方法 |
| WO2017043143A1 (ja) * | 2015-09-10 | 2017-03-16 | 株式会社日立製作所 | 分散電源発電量推定装置および方法 |
| CN105606893A (zh) * | 2016-01-26 | 2016-05-25 | 江苏科技大学 | 基于空间平滑修正music的电力间谐波检测方法 |
| CN107545253A (zh) * | 2017-09-11 | 2018-01-05 | 国网江西省电力公司电力科学研究院 | 一种电力系统低频振荡监测方法 |
Non-Patent Citations (1)
| Title |
|---|
| 张经纬 等: "基于四阶累积量的多信号分类法间谐波检测研究", 《继电器》 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110830133A (zh) * | 2019-12-23 | 2020-02-21 | 华中科技大学 | 一种基于多阶信道预测方法、预测系统及应用 |
| CN113032716A (zh) * | 2019-12-24 | 2021-06-25 | 南京理工大学 | 基于加窗插值和Prony算法的谐波与间谐波分析方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109557367B (zh) | 2020-09-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Su et al. | Power harmonic and interharmonic detection method in renewable power based on Nuttall double‐window all‐phase FFT algorithm | |
| Zygarlicki et al. | A reduced Prony's method in power-quality analysis—parameters selection | |
| Yao et al. | Fast S-transform for time-varying voltage flicker analysis | |
| CN107102255B (zh) | 单一adc采集通道动态特性测试方法 | |
| CN110648088B (zh) | 一种基于鸟群算法与svm的电能质量扰动源判断方法 | |
| Jain et al. | An adaptive time-efficient technique for harmonic estimation of nonstationary signals | |
| CN106483374A (zh) | 一种基于Nuttall双窗全相位FFT的谐波间谐波检测方法 | |
| CN115480203B (zh) | 一种电流互感器误差状态在线定量评估方法及系统 | |
| Belega et al. | Accuracy of the normalized frequency estimation of a discrete-time sine-wave by the energy-based method | |
| Tse | Practical application of wavelet to power quality analysis | |
| Jaiswal et al. | FDST‐based PQ event detection and energy metering implementation on FPGA‐in‐the‐loop and NI‐LabVIEW | |
| Molina‐Moreno et al. | Time domain harmonic state estimation in unbalanced power networks based on optimal number of meters and the principle of half‐wave symmetry | |
| CN109557367A (zh) | 一种高频分辨率谐波和间谐波Prony方法及装置 | |
| Petrović et al. | Computational effective modified Newton–Raphson algorithm for power harmonics parameters estimation | |
| Zolfaghari et al. | Evaluation of windowed ESPRIT virtual instrument for estimating Power Quality Indices | |
| CN107315103A (zh) | 一种电力冲击负载检测方法 | |
| CN111353415B (zh) | 一种脉冲响应中谐波成分的检测方法 | |
| Afrandideh et al. | Two modified DFT‐based algorithms for fundamental phasor estimation | |
| Lin | Separation of adjacent interharmonics using maximum energy retrieving algorithm | |
| Alaca et al. | CNN-Based Phase Fault Classification in Real and Simulated Power Systems Data | |
| Juhlin et al. | Fast gridless estimation of damped modes | |
| Lin | Identification of interharmonics using disperse energy distribution algorithm for flicker troubleshooting | |
| Ding et al. | Improved sparse component analysis for multi-point harmonic contribution evaluation under incomplete measurements | |
| CN114184838A (zh) | 基于sn互卷积窗的电力系统谐波检测方法、系统及介质 | |
| Jianze et al. | Time-varying transient harmonics measurement based on wavelet transform |
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 |
Granted publication date: 20200908 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |