基于加窗插值和Prony算法的谐波与间谐波分析方法

技术领域

本发明属于电力系统信号分析领域，特别是一种基于加窗插值和Prony算法的谐波与间谐波分析方法。

背景技术

电能是由电力企业向用户提供的一种特殊商品，在理想情况下，公用电网是以恒定的频率、标准的正弦波形、规定的电压水平和对称的三相电压电流向用户供电，由电厂、电网和用电方共同保证电能质量。但在实际电网中，由于系统中的发电机、变压器、输电线路等设备的固有缺陷及非线性负荷的大量使用等原因，使得电力企业提供的电压、电流信号波形和理想波形相比存在着一定程度的偏差，当偏差较大时，会给电网带来严重危害。电压、频率和波形是电能质量的三个因素，如果它们任意一个不在限定的范围内，就会给电力系统和用户带来不同程度的危害。

近年来，随着电力电子设备和冲击性负荷的增多，使得电网中谐波及间谐波的含量大量增加。谐波的存在会降低变压器容量，增加电机损耗，干扰电气设备正常工作，严重时甚至会损毁设备。与传统谐波相比，间谐波可以是直流到高次谐波间的任意频率，因此间谐波不但具有谐波的危害，其特性还会导致电压闪变、继电器误动作及引起无源滤波器过载等多种危害。同时，随着分布式电源的兴起和智能电网的快速发展，间谐波的来源越来越多，对电力系统的安全运行造成严重的威胁，因此对间谐波的治理迫在眉睫。然而，治理间谐波的前提条件是能够准确检测出间谐波参数，只有分析出间谐波信号，才能对其进行抑制和治理，进而改善供电环境，确保供电设备正常运转。因此，准确快速的检测出间谐波，对于抑制间谐波进而改善提高电能质量具有重要的现实意义。

目前文献中一般使用FFT插值算法来分析谐波信号，当信号中不含频率相近的谐波和间谐波信号时，现有FFT方法及其改进算法的具有较高的精度和稳定性，但当待分析信号中含有频率相近的谐波和间谐波成分时，该方法的频率分辨率不能满足需求。而Prony算法可使用一组指数函数的线性组合来拟合待分析信号，在理论上具有无限小的频率分辨率，当电网信号为宽带信号且所含谐波和间谐波成分较多时，算法准确性会受到一定影响，因此将两种算法结合可以弥补两种算法的局限性，同时，能保证精确度的情况下，准确分析出频率相近的间谐波信号。

发明内容

本发明的目的在于提供一种结构简单、运算速度高、可靠性高的基于加窗插值和Prony算法的电力系统谐波与间谐波分析方法。

实现本发明目的的技术解决方案为：一种基于加窗插值和Prony算法的谐波与间谐波分析方法，包括以下步骤：

步骤1、分析比较基本窗函数的频域特性，选择性能最优的窗函数；

步骤2、对待测信号进行加窗得到加窗信号，对加窗信号进行FFT变换，得到信号的频谱；

步骤3、在频谱中搜索谱峰，仅搜索幅值大于设定阈值的谱线，计算每个谱峰两根相邻谱线之间的相位差；

步骤4、判断主瓣干扰是否存在：

若主瓣干扰不存在，则使用插值算法校正对应的谐波或间谐波成分的参数；

若主瓣干扰存在，则不处理当前主瓣，利用插值算法算出的参数构建一个新信号，用原信号减去新信号得到残差信号，用Prony算法分析残差信号，得到含谐波和间谐波的信号的参数值。

进一步地，步骤1所述的分析比较基本窗函数的频域特性，选择性能最优的窗函数，具体如下：

步骤1.1、hanning窗也称升余弦窗，长度为N的离散hanning窗时域表达式w

其中，n＝0,1,2,…,N-1，w

式中，W

步骤1.2、hamming窗也称改进升余弦窗，长度为N的离散hamming窗时域表达式w

其中，n＝0,1,2,…,N-1，w

式中，W

步骤1.3、长度为N的离散Blackman窗时域表达式w

其中，n＝0,1,2,…,N-1，w

式中，W

步骤1.4、长度为N的离散4项3阶Nuttall窗的时域表达式w

其中，n＝0,1,2,…,N-1，4项最低旁瓣Nuttall窗的DTFT为W

式中，W

步骤1.5、根据以下原则选取窗函数：1)能量尽可能集中在主瓣内；2)旁瓣峰值电平尽可能低且衰减快；3)窗函数时域、频域表达式尽可能简洁。

进一步地，步骤2所述的对待测信号进行加窗得到加窗信号，对加窗信号进行FFT变换，得到信号的频谱，具体如下：

步骤2.1、设定待分析信号为：

式中，I为谐波个数，i表示为第i次谐波，i＝1时，f

步骤2.2、对x(t)以采样频率f

式中，n＝0,1,…,N-1，N为采样点数；

步骤2.3、用Nuttall窗w(n)对式(11)的信号x(n)进行加窗处理，得到加窗后信号的离散傅里叶变换为X(kΔf)：

式中，k为一谐波处的采样点值，Δf为频率分辨率，那么kΔf则为该谐波处的频率；用j表示虚数，i表示谐波次数；

步骤2.4、忽略负频点处频峰的旁瓣影响，则在正频点周边的离散谱线函数为：

式中，离散谱线频率分辨率Δf＝f

其中：

式中，a

设定信号的频点k

进一步地，步骤3所述的在频谱中搜索谱峰，仅搜索幅值大于设定阈值的谱线，计算每个谱峰两根相邻谱线之间的相位差，具体如下：

步骤3.1、由式(14)和式(15)得：

步骤3.2、得k

arg[X(k

arg[X(k

arg[X(k

由式(19)至式(21)知，k

进一步地，步骤4所述的判断主瓣干扰是否存在：若主瓣干扰不存在，则使用插值算法校正对应的谐波或间谐波成分的参数，具体如下：

步骤4.1、设置可动态调整的参量，通过增加一个临界信号来获得临界值；IEC61000-3-6规定当信号的幅值小于基波幅值的0.1％时，忽略此信号，故将残余信号的终止条件б

步骤4.2、判断无其他谐波谱线干扰，则用插值法进行计算，具体如下：

设

记式(22)为β＝h(α)，其反函数记为α＝h

f

步骤4.3、第i次谐波的幅值修正A

当N值大于阈值时，式(24)简化为：

A

利用多项式拟合逼近，得A

A

式中，m(α)为h(α)的拟合逼近式；

步骤4.4、由式(13)得相位修正公式为：

其中，Nuttall窗函数的三谱线插值7阶逼近式为：

α＝h(β)＝1.724339β-0.085167λ

m(α)＝1.724339+0.430783λ

进一步地，步骤4所述的若主瓣干扰存在，则不处理当前主瓣，利用插值算法算出的参数构建一个新信号，用原信号减去新信号得到残差信号，用Prony算法分析残差信号，得到含谐波和间谐波的信号的参数值，具体如下：

步骤5.1、采用p个具有任意幅值、相位、频率和衰减因子的指数函数的线性组合对振荡数据进行等间距采样，设N个原始信号为x(0),x(1),…,x(n),有数学模型：

式中，n＝0,1,…,N-1，A

步骤5.2、定义样本函数r(i,j)为：

式中，x(·)为原始信号，x*(·)为其共轭值；

利用样本函数r(i,j)构建矩阵为：

步骤5.3、线性参数估计看做求解方程：

R

式中，a

对R

0＝…σ

式中，σ

由于实际中存在噪声，且M>p，因此用噪声空间取代样本矩阵p-M维零空间；p-M为估计阶数与实际阶数的差值；

定义系统阶数估计w

式中，p为系统的估计阶数，由于w

步骤5.4、用零空间代替系统的噪声空间，得到样本矩阵R

A

式中：

其中，

步骤5.5、由式(36)求解预测参数组成的多项式得到z

步骤5.6、结合步骤4求出的不含主瓣干扰的谐波信号的频率、幅值、相角，计算含谐波和间谐波的信号的频率、幅值、相角。

本发明与现有技术相比，其显著优点在于：(1)分析比较基本常用窗函数的频域特性，选择性能较好的窗函数，对待测信号进行加窗得到加窗信号，对加窗信号进行FFT变换，得到信号的频谱，结构简单；(2)在频谱中搜索谱峰，计算每个谱峰附近两根谱线之间的相位差；若主瓣干扰不存在，则使用插值算法校正对应的谐波或间谐波成分的参数；若主瓣干扰存在，则暂不处理当前主瓣，运算速度快、可靠性高；(3)利用插值算法算出的参数构建一个新信号，用原信号减去新信号得到残差信号，用Prony算法分析残差信号，准确快速的得到频率相近的谐波和间谐波的参数，从而有效的分析信号中谐波与间谐波的参数，对于抑制间谐波进而改善提高电能质量具有重要的现实意义。

附图说明

图1是本发明加窗插值和Prony算法的谐波与间谐波分析方法的流程示意图。

图2是本发明实施例中部分常用窗函数的归一化频率特性图。

图3是本发明实施例中含频率相近的谐波与间谐波信号的幅频谱图。

具体实施方式

下面结合附图及具体实施例对本发明作进一步详细说明。

结合图1，本发明基于加窗插值和Prony算法的谐波与间谐波分析方法，包括以下步骤：

步骤1、分析比较基本常用窗函数的频域特性，选择性能较好的窗函数，结合图2，具体如下：

步骤1.1、hanning窗也称升余弦窗，长度为N的离散hanning窗时域表达式w

其中，n＝0,1,2,…,N-1，w

式中，W

步骤1.2、hamming窗也称改进升余弦窗，长度为N的离散hamming窗时域表达式w

其中，n＝0,1,2,…,N-1，w

式中，W

步骤1.3、长度为N的离散Blackman窗时域表达式w

其中，n＝0,1,2,…,N-1，w

式中，W

步骤1.4、长度为N的离散4项3阶Nuttall窗的时域表达式w

其中，n＝0,1,2,…,N-1，4项最低旁瓣Nuttall窗的DTFT为W

式中，W

步骤1.5、根据以下原则选取窗函数：1)能量尽可能集中在主瓣内；2)旁瓣峰值电平尽可能低且衰减快；3)窗函数时域、频域表达式尽可能简洁。

但窗函数无法同时具备主瓣窄、旁瓣峰值电平低且旁瓣衰减快的条件。因此，在采用FFT谐波检测算法时要综合考虑实时性和准确度两方面的要求。如果仅要求精确计算频率，而不考虑幅值、初相角等参数的准确度，则可选用主瓣宽度比较窄而便于分辨的矩形窗，例如测量物体的自振频率等；如果分析窄带信号，且有较强的干扰噪声，则应选用旁瓣幅度小的窗函数，如Hanning窗、三角窗等。本发明分析信号中谐波与间谐波的参数，综合考虑选择Nuttall窗。

步骤2、对待测信号进行加窗得到加窗信号，对加窗信号进行FFT变换，得到信号的频谱，具体如下：

步骤2.1、设定待分析信号为：

式中，I为谐波个数，i表示为第i次谐波，i＝1时，f

步骤2.2、对x(t)以采样频率f

式中，n＝0,1,…,N-1，N为采样点数；

步骤2.3、用Nuttall窗w(n)对式(11)的信号x(n)进行加窗处理，得到加窗后信号的离散傅里叶变换为X(kΔf)：

式中，k为某一谐波处的采样点值，Δf为频率分辨率，那么kΔf则为该谐波处的频率。用j表示虚数，i表示谐波次数。

步骤2.4、忽略负频点处频峰的旁瓣影响，则在正频点周边的离散谱线函数为：

式中，离散谱线频率分辨率Δf＝f

其中：

式中，a

设定信号的频点k

步骤3、在频谱中搜索谱峰，仅搜索幅值大于设定阈值的谱线，计算每个谱峰两根相邻谱线之间的相位差，具体如下：

步骤3.1、由式(14)和式(15)得：

步骤3.2、得k

arg[X(k

arg[X(k

arg[X(k

由式(19)至式(21)知，k

步骤4、判断主瓣干扰是否存在：

(1)若主瓣干扰不存在，则使用插值算法校正对应的谐波或间谐波成分的参数，具体如下：

步骤4.1、设置可动态调整的参量，通过增加一个临界信号来获得临界值；IEC61000-3-6规定当信号的幅值小于基波幅值的0.1％时，忽略此信号，故将残余信号的终止条件б

步骤4.2、判断无其他谐波谱线干扰，则用插值法进行计算，具体如下：

设

记式(22)为β＝h(α)，其反函数记为α＝h

f

步骤4.3、第i次谐波的幅值修正A

当N值较大时，式(24)可简化为：

A

利用多项式拟合逼近，可得A

A

式中，m(α)为h(α)的拟合逼近式。

步骤4.4、由式(13)得相位修正公式为：

其中，Nuttall窗函数的三谱线插值7阶逼近式为：

α＝h(β)＝1.724339β-0.085167λ

m(α)＝1.724339+0.430783λ

(2)若主瓣干扰存在，则不处理当前主瓣，利用插值算法算出的参数构建一个新信号，用原信号减去新信号得到残差信号，用Prony算法分析残差信号，得到含谐波和间谐波的信号的参数值，如图3所示，具体如下：

步骤5.1、采用p个具有任意幅值、相位、频率和衰减因子的指数函数的线性组合对振荡数据进行等间距采样，设N个原始信号为x(0),x(1),…,x(n)，有数学模型：

式中，n＝0,1,…,N-1，A

步骤5.2、定义样本函数r(i,j)为：

式中，x(·)为原始信号，x*(·)为其共轭值。

利用样本函数r(i,j)构建矩阵为：

步骤5.3、线性参数估计看做求解方程：

R

式中，a

对R

0＝…σ

式中，σ

由于实际中存在噪声，且M>p，因此用噪声空间取代样本矩阵p-M维零空间；p-M为估计阶数与实际阶数的差值。

定义系统阶数估计w

式中，p为系统的估计阶数，由于w

步骤5.4、用零空间代替系统的噪声空间，得到样本矩阵R

A

式中：

其中，

步骤5.5、由式(36)求解预测参数组成的多项式得到z

步骤5.6、结合步骤4求出的不含主瓣干扰的谐波信号的频率、幅值、相角，至此可以全部计算含谐波和间谐波的信号的频率、幅值、相角。

综上所述，本发明线路简单，硬件成本低，可靠性高，能够有效分析出信号中谐波与间谐波的参数。