本发明属于电力系统谐波信号检测技术领域，尤其涉及一种基于多项式拟合及非干扰区域划分的间谐波检测方法。
背景技术：
随着用户对电能质量要求的提高，如何精确地检测并过滤谐波显得尤为重要。目前，有源滤波技术能够较好地治理谐波问题。该技术中一个重要环节是采用离散傅里叶变换(dft)处理信号，从而分析出谐波相关的信息。
近年来，各种高级窗函数以及各种复杂的插值公式被引入傅里叶变换的线路中，常规谐波的检测精度已能够达到较为理想水平。“频谱泄漏”以及“栅栏效应”所造成的误差均能被有效抑制。然而，这些工具在面对间谐波所造成的“主瓣干扰”时，其检测进度仍有较大的改善空间。存在的问题主要有两方面：
(1):主瓣干扰会对常规谐波的谱线造成巨大影响。这使得插值过程中所选用的谱线无法准确描述该谐波的信息，导致所计算出的结果与理论值之间存在较大偏差。
(2):目前许多方法的思路，都是从混合的谱线中分离出两种谐波谱线。这类“分离主瓣式”的方法未必能够精准、完美地分离两种谐波。无论是分离的量不足，还是分离的量过多，都会对检测结果的精度造成影响。同时，由于谐波检测对结果的精度要求极高，此类方法未必能够充分满足其精度要求。
考虑以上因素，本发明专利公开一种新的间谐波检测方法。该方法针对间谐波所造成的主瓣干扰问题，采用多项式拟合以及划分非干扰区域等手段，提高间谐波的参数检测精度。
技术实现要素：
为解决间谐波与相邻正弦波频率相近所造成的主瓣干扰现象，本发明公开了一种基于多项式拟合及非干扰区域划分的间谐波检测方法，该方法能够有效分辨频域内的主瓣干扰现象，且能够分别检测出相互干扰的两个正弦波的参数，具有较好的间谐波检测能力。
本发明采取的技术方案为：
基于多项式拟合及非干扰区域划分的间谐波检测方法，包括以下步骤：
步骤1：通过电流互感器，将电网中的电流转化为二次侧可供仪器测量的信号，再通过采样获得一定长度的离散信号。
步骤2：构造与离散信号长度一致的唯主瓣类窗函数。并对该窗的主瓣进行多项式拟合，得到主瓣频率与响度之间的关系。
步骤3；将步骤2所构造出的窗函数与离散信号相乘，并对乘积进行离散傅里叶变换(dft)，获得信号的频谱信息。
步骤4；搜索频域上各峰值，通过类信噪比公式计算出峰值附近谱线相较于峰值的响度，所得响度图类似“小包”状。
步骤5；以一定数值的响度作为阈值截断各小包，通过所截断的宽度判断该频域内是否存在主瓣干扰。
步骤6；对于存在主瓣干扰的小包，将其一侧的阈值以下、-300db以上的谱线按对称轴翻转至另一侧。并以对称的结果为依据，划分出非干扰区间。
步骤7；采用非干扰区间中的谱线建立插值公式，计算出偏移量δ，继而求解出谐波对应的频率、幅值、相角。
所述步骤1中，所述离散信号为电流信号、或者电压信号。
所述步骤2中，唯主瓣类窗特征为：仅主瓣的峰值高于-300db，所有旁瓣峰值均低于-300db。允许零星波动的频点数值高于-300db，但波动值不得高于-280db。
所述步骤2中，多项式拟合为：以各谱线数值相对于k＝0hz谱线数值的响度作为自变量；以该谱线对应的频率作为因变量。
所述步骤2中，多项式拟合以各谱线的数值相对于谱线数值的响度作为自变量，以谱线对应的频率作为因变量，通过最小二乘法计算出其表达式。其结果形似：
其中：
h为拟合项数，ch为每项相应的系数；
h为最大项数，该参数的选择可视检测精度而定，一般7至10次即可；
x为主瓣上谱线的响度，单位为db；
r(x)为响度为x的谱线对应的频率，单位为hz。
所述步骤3中，离散傅里叶变换的所选参数与窗函数的阶数有关，以8阶blackman-harris自卷积窗为例，所采用的dft公式为：
其中，x(n)为离散信号；w(n)为离散的窗函数。
所述步骤4中，类信噪比公式为：g＝20·lg(s/n)。其中，s为各谱线数值，n为局部最大值。
所述步骤5中，阈值为：其数值比旁瓣的响度高80db左右，且能够保证所有谐波的旁瓣在波动状态下，均无法超过该数值。
所述步骤5中，判断存在主瓣干扰的依据为：以阈值截断单个正弦波频域小包所得的宽度为基础，向上增加1至2根谱线作为裕度。若实际截断宽度大于该数值，则存在主瓣干扰；反之则不存在主瓣干扰。
所述步骤6中，按对称轴翻转为：通过步骤2建立的关系式求解出对称轴。将小包单侧的谱线：其数值在阈值以下、-300db以上，以对称轴为基准进行对称。
所述“按对称轴翻转”分为两个步骤：
step1对称轴的构造部分，其特征为：以小包某一侧的谱线(其数值在阈值以下、-300db以上)数值作为基础，通过多项式拟合计算出对应的“宽度”。然后将谱线的频率向小包内侧偏移相应宽度。平均各偏移结果，将均值所对应的谱线作为对称轴。
step2谱线对称部分，其特征为：将小包某一侧的谱线(其数值在阈值以下、-300db以上)按照对称轴，对称至另一侧。
所述步骤6中，边界划分方法为：将对称结果最外侧的谱线，以及结果同侧的第一根响度低于-300db的谱线作为边界。
所述步骤7中，建立插值公式为：所采用谱线位于非干扰区域中，且位于偏向小包内部的一侧。
本发明一种基于多项式拟合及非干扰区域划分的间谐波检测方法，有益效果如下：
(1)、能够辨明主瓣干扰：
本发明所公布的算法通过截断频域上的小包，来判断是否存在主瓣干扰。该方法以不受主瓣干扰的正弦波的截断宽度为依据，计算出划分阈值。该阈值能够定量比较实际的截断宽度，从而准确判断是否存在主瓣干扰。
(2)、对频谱泄露抗性强：
通常而论，唯主瓣类窗在频域中的旁瓣数值普遍较低。这一特性有利于减少频谱泄漏，提高谐波参数检测精度。
(3)、检测精度高：
相比与常规的“分离主瓣式”间谐波检测方法，本发明所公布的方法通过划分非干扰区域，避开了主瓣干扰的影响范围，同时也不会产生由于分离不够充分或分离过于激进带来的计算误差。因此该方法具有较高检测精度。
附图说明
图1为本发明检测方法流程图。
图2(a)为三种唯主瓣类窗函数的整体频谱图。
图2(b)为三种唯主瓣类窗函数旁瓣处的局部放大图。
图3为以一定数值的响度作为阈值截断频域中的小包示意图。
图4为通过多项式拟合构造的对称轴、以及谱线对称的过程图。
图5为非干扰区域内的谱线与理论频点之间的关系图。
具体实施方式
基于多项式拟合及非干扰区域划分的间谐波检测方法，包括以下步骤：
步骤1；获取含有间谐波的离散信号。
此处定义该信号如式(1)所示：
其中：
n为信号总长度，各采样点n＝0,1,2,...,n-1。
fz、az、
ts为相邻两个采样点之间的时间间隔，若采样频率为fs，则ts＝1/fs。
步骤2；构造与离散信号一致的唯主瓣类窗函数。
此处以8阶blackman-harris自卷积窗为例，其表达式为：
w(n)＝[(b*b)*(b*b)]*[(b*b)*(b*b)](2)
其中符号b指长度为n/8单个blackman-harris窗，表达式见公式(3)所示。
其中：
m为单个窗总长度，有m＝n/8；
aq(q＝0,1,2,3)为blackman-harris窗的各项系数：a0＝0.35875、a1＝0.48829、a2＝0.14128、a3＝0.01168。
将8个窗函数按照式(2)的次序卷积后，即可得到最终的8阶卷积窗w(n)。对其进行离散傅里叶变换，如式(4)：
其中，w为窗函数信号；m为乘积长度的八分之一；k为频率。
此处有两点值得注意：
1.在卷积窗函数的过程中，每次卷积完成后须在结果末尾增补一个零点，以保证经过8次卷积后，所得结果长度正好等于n。
2.上述dft公式中，指数函数中的分母由常规的n调整为m。其目的在于方便离散卷积定理在插值过程中的运用。
获得频域内的谱线后，采用如式(5)所示的类信噪比公式处理各谱线：
g＝20·lg(s/n)(5)
其中，s为各谱线的数值，n为k＝0所对应谱线(峰值谱线)的数值。计算所得结果量纲与响度一致，单位应为“分贝(db)”。
通过式(5)可得到各谱线相对于峰值的响度图。以各响度g作为自变量，所对应谱线的频率作为因变量，拟合出多项式如式(6)所示：
其中：
h为多项式次数，ch为对应的系数；
h为最大次数，该参数的选择可视检测精度的要求而定，一般7至10次即可；
x为主瓣上谱线的响度(单位为db)；
r(x)为响度为x的谱线对应的频率(单位为hz)。
步骤3；将窗函数与信号相乘，获得其频谱信息。
由于窗函数与离散信号具有相同的长度n，因此二者可以相乘。对其乘积进行离散傅里叶变换如式(7)，则可得到其频域信息：
同样考虑到离散卷积定理的使用，式(7)中的分母与式(4)的一致。
值得注意的是，后续的阈值截断部分对频谱分辨率要求较高。故此处频点的采集应尽量致密。拟定相邻谱线之间的间隔为d＝1/64。
步骤4；搜索频域上各峰值，计算出峰值附近谱线相较于峰值的响度。
首先，找到频域中各有效峰值。
然后通过类信噪比公式，同式(5)，计算其附近谱线的数值相较于该数值的响度。所得响度图类似“小包”状。对于采用类信噪比公式处理的谱线而言，其覆盖范围不必过大。频率范围能够完全覆盖整个小包，且最左、最右范围内的谱线趋于平整即可。
步骤5；拟定阈值，截断各小包，通过所截断宽度判断是否存在主瓣干扰。
观察8阶blackman-harris窗的频谱图可知，其旁瓣的响度普遍在-300db以下。实际检测过程中，基波的幅值通常比谐波幅值高出近百倍，折算为响度，约为40db。故阈值至少应拟定为-260db，以防止基波的旁瓣对谐波的主瓣造成干扰。
另一方面，由于窗函数的旁瓣较为波动，阈值还应包括一定裕度。考虑以上信息，此案例中的阈值最终拟定为snr＝-220db。以该数值为依据截断频域中的小包(计算该频域时频点步长取d＝1/64)，在不存在主瓣干扰的情况下，所得截断宽度为172根谱线至173根谱线不等。因此，判断存在主瓣干扰的依据为：
若某小包被snr＝-220db截断后宽度大于或等于174，则该频点附近存在主瓣干扰。
若某小包被snr＝-220db截断后宽度小于174，则该频点附近不存在主瓣干扰。
步骤6；对于存在主瓣干扰的小包，将部分谱线对称并划分非干扰区域。
对于存在主瓣干扰的小包，首先搜索其一侧位于snr＝-220db与snr＝-300db之间的所有谱线。以右侧的谱线为例，所得频点集合记为：
pi＝{(xi,yi)|-220>yi>-300}(8)
其中，pi为点集，下标为i的点所对应的坐标为(xi,yi)。要求该点集中所有点的y坐标均在-220db至-300db之间。
将各点的响度(yi)作为自变量带入多项式(6)，即可算出理论主瓣峰值距离这些谱线的距离(ri)。然后，将各谱线的频率向小包内侧偏移等同于该距离的长度，如式(9)所示：
mi＝xi-ri(9)
每一个属于集合p内的频点(xi,yi)都能够计算出对应的mi值。平均这一数值，即可得到理论峰值所在的频率，如式(10)。
以此频率作为对称轴，将集合中所有频点对称，即可得到对称的结果。该对称结果位于小包左侧，记为式(11)：
其中，qi为点集，下标为i的点所对应的坐标为(x′i,y′i)；该点集中所有点由点集pi中的点变换而来。
以该集合内最外侧(也是数值最小的谱线)的谱线x＝min(x′i)和原本位于小包左侧、响度低于-300db的第一根谱线为依据，划分出非干扰区域。
从这一划分方案不难看出，处于该区域内的谱线具有两方面特点：
(1)位于谱线x＝min(x′i)左侧；这一特点表明，右侧正弦波的主瓣对该谱线的影响已经衰减至-300db以下。因此基本可认为该部分谱线不受主瓣干扰。
(2)所处位置的响度高于-300db；这一特点表明，旁瓣对谱线的干扰不高，以该区域内的谱线构筑插值公式，计算出的结果具有较高精度。
步骤7；建立插值公式，计算偏移量δ，并求解出谐波参数。
该步骤为常规谐波检测过程中的必要步骤，在此以双谱线公式作简单说明。取位于非干扰区域内、最偏向小包内侧的两根谱线(从左至右记为y1、y2)建立双谱线插值公式如式(12)：
然后根据理论频点与这两根谱线频率之间的关系，计算出偏移量δ。最后，按照式(13)计算出各次谐波的频率、幅值以及相角：
其中，k1为谱线y1对应的频率；arg[·]表示括号内相量的角度；az、fz、
至此，两个相互干扰的正弦波中，左侧正弦波所对应的参数已经求解出。再以小包左侧谱线为例，重复步骤6及步骤7，即可求解出右侧正弦波对应的参数。
图2(a)为三种唯主瓣类窗函数的整体频谱图。观察图2(a)发现，唯主瓣类窗在频域中仅主瓣峰值位于-300db以上。因此距离该主瓣较远的部分基本不受其影响。这一特征有利于分析谐波的影响范围。
图2(b)为三种唯主瓣类窗函数旁瓣处的局部放大图；从图2(b)中可知，唯主瓣类窗的旁瓣数值普遍较低，因此这类窗也能够较好地控制频谱泄露所造成的误差。
图3为以一定数值的响度作为阈值截断频域中的小包示意图。观察图3发现，当频域中出现主瓣干扰现象时，一个显著的特征就是，存在主瓣干扰的小包的截断宽度，远比不存在主瓣干扰的小包宽。因此，以-220db作为阈值对小包进行截断，并以其截断处的宽度作为依据判断是否存在主瓣干扰现象这一方法十分合理。
图4为通过多项式拟合构造的对称轴、以及谱线对称的过程图。图4详细地展示了以对称轴为基准进行对称，且以对称结果划分非干扰区域的过程。从图中可看出，谱线的对称结果与谐波本身的轮廓较为一致。这说明以对称结果作为依据划分非干扰区域的方案十分合理，该算法的思路能够治理主瓣干扰问题。
图5为非干扰区域内的谱线与理论频点之间的关系图。从图5中可知，在建立插值公式时，宜采用位于非干扰区域内、偏向小包内侧的谱线。同时应注意到，须用左侧的谱线计算左侧相应的谐波；以右侧谱线计算右侧相应的谐波。这一对应关系不可混用。