一种基于主瓣宽度的多谱线插值谐波分析方法及系统与流程

本发明涉及谐波分析领域，具体来讲涉及一种基于主瓣宽度的多谱线插值谐波分析方法及系统。
背景技术：
：电网中谐波的问题日益严重，不仅恶化电能质量，对电网的安全稳定和经济运行也造成较大影响。因此，对电网中谐波参数进行高精度分析与测量，对于提高电能质量，减少谐波危害，维护电网安全稳定、高效运行是十分必要的。电力系统谐波分析最常用的算法是采用fft(fastfouriertransform，快速傅里叶变换)分析算法。但由于电网频率的不固定，fft很难做到同步采样和整周期截断，所带来的频谱泄露和栅栏效应可能会导致谐波的频率、幅值和相位等参数测量不准，无法满足测量要求。加窗插值fft算法被广泛用来解决这一问题，具体做法是：运用各种窗函数对信号进行截断，然后结合谱线插值fft进行谐波分析，从而提高测量精度。常用窗函数例如汉宁(hanning)窗函数、布莱克曼(blackman)窗函数、布莱克曼汉斯(blackman-harris)窗函数、纳托尔(nuttall)窗函数、莱夫文森特(rife-vincent)窗函数以及各种组合窗函数。随着复杂组合窗函数的出现，主瓣宽度越来越宽，主瓣内的所含谱线数目越来越多，常规的双插值和三插值会出现丢失有用信息的情况，导致谐波分析精度降低。技术实现要素：针对现有技术中存在的缺陷，本发明的目的在于提供一种基于主瓣宽度的多谱线插值谐波分析方法及系统，有效提高谐波分析的精度。为达到以上目的，本发明提供一种基于主瓣宽度的多谱线插值谐波分析方法，包括步骤：s1.上位机接收到互感器采集到的离散的电网信号x(n)，n为采样点的序数，n为自然数；采用离散余弦窗函数w(n)对电网信号x(n)进行加窗截断，得到加窗信号xw(n)：xw(n)＝x(n)w(n)式(1)对式(1)的加窗信号xw(n)进行fft变换后，忽略负频率点处旁瓣的影响，得到加窗fft频谱x(k)：其中w(·)为窗函数的频谱，k为正整数，x(k)表示第k次谐波的频谱，ak为第k次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，k0为真实频谱的谱线位置；s2.整数次谐波之间间隔为50hz，非同步采样带来的泄露频谱线在谐波间隔中最多有h1根，h1＝int(50/δf)，其中δf为离散频率间隔；主瓣内的谱线根数h2为h2＝int(δω/2πδf)，多谱线的根数h确定为h＝min(h1,h2)；设峰值频率点左右多条的谱线分别为k1＜…＜kh＜k＜kh+1＜…＜kh，对应的谱线幅值分别为y1，…，yh，yh+1，…，yh；记α＝k-kh-0.5，由于0≤k-kh≤1，则-0.5≤α≤0.5，另记其中α和β均为系数，ci为加权系数，i＝1…h，ci为yi的权值，令ci＝c(h-i+1)，根据离散余弦窗函数w(n)和频域主瓣表达式，加余弦窗时得到s4.根据式(4)和式(5)计算多谱线插值算法的修正式；s5.计算基波参数，包括基波频率f0和基波幅值a1，根据式(4)得出基波的相位s6.确定基波频率f0后，在范围(kf0-5,kf0+5)内重复步骤s2～s5，直到所有谐波参数计算完毕。在上述技术方案的基础上，离散余弦窗函数w(n)的表达式为其中，n为采样点数，n为正整数，n＝0,1,2…n-1；σ表示求和；m为窗函数的累加次数，m＝0,1,2…m-1；m为窗函数项数，m为正整数；bm为窗函数系数。在上述技术方案的基础上，所述窗函数系数bm满足约束条件在上述技术方案的基础上，对式(1)的加窗信号进行fft变换后，得到加窗fft频谱其中，第一次谐波为基波，f0为基波频率，且δf＝fs/n。在上述技术方案的基础上，所述s4中，根据式(4)和式(5)得到ci为加权系数，i＝1…h，ci为yi和ym-i的权值，当n大于1000时，式(7)简化为β＝g(α)，其反函数为α＝g-1(β)；采用多项逼近方法计算奇函数β＝g-1(α)，α≈p11×β+p13×β3+…+p1pβp，其中，p11，p13，…，p1p为多项式逼近的奇次项系数，由式(4)可知信号幅值为得到通过逼近方法，当n大于1000时，窗函数系数为实系数，式(10)可表示为u(·)为偶函数，逼近多项式不含奇次项，则多谱线修正逼近多项式为u(α)＝(p20+p22α2+…+p2dαd)式(11)其中，p20，p22，…，p2d为多项式逼近的偶次项系数。在上述技术方案的基础上，所述s5中，f0＝k·δf＝(k2+α+0.5)δf，在上述技术方案的基础上，所述电网信号包括电流信号、电压信号，所述真实频谱的谱线位置k0为小数。本发明还提供一种基于主瓣宽度的多谱线插值谐波分析系统，设置于上位机，包括信号预处理单元、谱线数目确定单元、权值确定单元、修正公式计算单元、基波参数计算单元、谐波参数确定单元；所述信号预处理单元，用于对上位机中来自互感器的离散的电网信号x(n)进行预处理；采用离散余弦窗函数w(n)对电网信号x(n)进行加窗截断，得到加窗信号xw(n)xw(n)＝x(n)w(n)式(1)对式(1)的加窗信号xw(n)进行fft变换后，忽略负频率点处旁瓣的影响，得到加窗fft频谱x(k)：其中w(·)为窗函数的频谱，k为正整数，x(k)表示第k次谐波的频谱，ak为第k次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，k0为真实频谱的谱线位置；所述谱线数目确定单元，用于确定多根谱线；整数次谐波之间间隔为50hz，非同步采样带来的泄露频谱线在谐波间隔中最多有h1根，h1＝int(50/δf)，其中δf为离散频率间隔；主瓣内的谱线根数h2为h2＝int(δω/2πδf)，多谱线的根数h确定为h＝min(h1,h2)；权值确定单元，用于计算多谱线所对应得权值大小；设峰值频率点左右多条的谱线分别为k1＜…＜kh＜k＜kh+1＜…＜kh，对应的谱线幅值分别为y1，…，yh，yh+1，…，yh；记α＝k-kh-0.5，由于0≤k-kh≤1，则-0.5≤α≤0.5，另记其中α和β均为系数，ci为加权系数，i＝1…h，ci为yi的权值，令ci＝c(h-i+1)，根据离散余弦窗函数w(n)和频域主瓣表达式，加余弦窗时得到所述修正公式计算单元，用于根据式(4)和式(5)计算多谱线插值算法的修正公式；所述基波参数计算单元，用于计算基波参数，包括基波频率f0和基波幅值a1，以及基波的相位所述谐波参数确定单元，用于确定基波频率f0后，在范围(kf0-5,kf0+5)内重复步骤s2～s5，直到所有计算得到所有谐波参数。在上述技术方案的基础上，系统还包括误差分析单元，在窗函数主瓣宽度不断增加，参与运算的谱线数目也相应增加的情况下，进行多谱线插值计算的误差分析。在上述技术方案的基础上，所述电网信号包括电流信号、电压信号，所述真实频谱的谱线位置k0为小数。本发明的有益效果在于：通过从电力系统的电网信号加窗后的频域表达式入手，根据所加窗函数的主瓣宽度确定参与运算的谱线数目，进而推导出各个谱线权值的大小。并在此基础上进行谐波分析，计算谐波相关参数。通过算例表明，本发明利用主瓣内所有谱线进行插值分析，提高了精度，实现了谐波参数的高精度测量，具有较好的使用价值。附图说明图1为本发明基于主瓣宽度的多谱线插值谐波分析方法流程图；图2为本发明实施例中谐波信号的波形图；图3为1阶hanning窗函数的幅频特性图；图4为2阶hanning窗函数的幅频特性图；图5为3阶hanning窗函数的幅频特性图；图6为4阶hanning窗函数的幅频特性图；图7为本发明基于主瓣宽度的多谱线插值谐波分析系统示意图。具体实施方式以下结合附图及实施例对本发明作进一步详细说明。如图1所示，本发明基于主瓣宽度的多谱线插值谐波分析方法，包括如下步骤：s1.信号预处理。具体的，上位机接收互感器采集到的离散的电网信号x(n)，n为采样点的序数，n为自然数；电网信号包括电流信号、电压信号。采用离散余弦窗函数w(n)对电网信号x(n)进行加窗截断，得到加窗信号xw(n)：xw(n)＝x(n)w(n)式(1)离散余弦窗函数w(n)的表达式为：其中，n为采样点数，n为正整数，n＝0,1,2...n-1；σ表示求和；m为窗函数的累加次数，m＝0,1,2...m-1；m为窗函数项数，m为正整数；bm为窗函数系数。窗函数系数bm满足约束条件如下：对式(1)的加窗信号xw(n)进行fft变换后，得到加窗fft频谱：其中，w(·)为窗函数的频谱，k为正整数，x(k)表示第k次谐波的频谱，ak为第k次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，第一次谐波为基波，fs为采样频率，f0为基波频率，δf为离散频率间隔，且δf＝fs/n。令：k0＝f0/δf，k0为真实频谱的谱线位置，真实频谱的谱线位置k0为小数，忽略负频率点处旁瓣的影响，式(3)变为：s2.根据主瓣宽度确定谱线数目。具体的，对于整数次谐波之间间隔为50hz，那么非同步采样带来的泄露频谱线在谐波间隔中最多有h1根，h1＝int(50/δf)。若所加窗函数的主瓣宽度为δωrad/s，主瓣内的谱线根数h2为h2＝int(δω/2πδf)，多谱线的根数h确定为h＝min(h1,h2)。s3.确定权值。具体的，令对信号的非同步采样或非整周期数据截断时，由于栅栏效应，信号的峰值频率f0＝k0·δf很难刚好位于离散谱线的频点上，即k0一般不为整数。设峰值频率点左右多条的谱线分别为k1＜…＜kh＜k＜kh+1＜…＜kh，对应的谱线幅值分别为y1，…，yh，yh+1，…，yh。记α＝k-kh-0.5，由于0≤k-kh≤1，则-0.5≤α≤0.5，另记：其中α和β均为系数，ci为加权系数，i＝1…h，ci为yi的权值。考虑到计算的方便，令ci＝c(h-i+1)，余弦窗函数的表达式如式(2)所示时，根据频域主瓣表达式，加余弦窗时候，可以求得：s4.计算多普线插值算法的修正公式。具体的，根据式(4)和式(5)得到其中，ci为加权系数，i＝1…h，ci为yi和ym-i的权值。当n大于1000时，式(7)可以简化为β＝g(α)，其反函数为α＝g-1(β)；由于所采用的余弦窗系数均为实系数，其频率响应是偶对称的，因而g(·)和g-1(·)均为奇函数。采用多项逼近方法计算奇函数β＝g-1(α)，表达式为：α≈p11×β+p13×β3+…+p1pβp式(8)其中，p11，p13，…，p1p为多项式逼近的奇次项系数，由式(4)可知信号幅值ai为：可以得到：通过类似于式(10)的逼近方法，当n大于1000时，窗函数系数为实系数，式(10)可表示为：u(·)为偶函数，逼近多项式不含奇次项，则多谱线修正逼近多项式为：u(α)＝(p20+p22α2+…+p2dαd)式(11)其中，p20，p22，…，p2d为多项式逼近的偶次项系数。s5.计算基波参数。具体的，包括计算基波频率f0、基波幅值a1：f0＝k·δf＝(k2+α+0.5)δf式(12)根据式(4)得出基波的相位s6.确定谐波参数。具体的，确定基波频率f0后，在范围(kf0-5,kf0+5)内重复步骤s2～s5，直到所有谐波参数计算完毕。上述步骤s1～s6完成后，可以增加一个误差分析的步骤，当窗函数主瓣宽度不断增加，参与运算的谱线数目也相应增加，进行多谱线插值算法的误差分析，确保高精度。下面通过一个具体实施例来进行详细说明。如图2所示，是一个选定电网信号谐波的波形，基于主瓣宽度的多谱线插值谐波分析方法包括：a1.信号预处理。具体的，将互感器采集到的离散电网信号x(n)，信号模型参见表1所示。为了验证所提算法的精度，进行10次谐波仿真分析。信号模型为：其中，基波频率f0为50.5hz，采样频率fs为5120hz，采样点数n为1024。表1、电网信号的谐波参数自乘法窗函数的定义为：若干个相同窗函数进行相乘所构成的新的窗函数。乘法窗函数的公式：其中，w(n)为基本窗函数，p为参加乘法的窗函数个数，称为乘法窗函数的阶数，p＝1时即为原窗函数。乘法窗函数的主瓣特性和旁瓣特性同基本窗函数选取和乘法窗的阶数相关。根据原窗函数的不同和乘法窗阶数的不同可以构造多种乘法窗函数。本实施例选取经典的hanning窗函数作为原窗函数构造新的乘法窗函数。加窗fft变换后得到信号频谱：a2.确定谱线数目。由于频域分辨率δf＝(2π)/n，p阶hanning乘法窗函数的主瓣宽度(mb)为(4(p+1)π)/n，主瓣内谱线根数为2(p+1)。hanning自乘法窗函数的频谱特性如表2所示。表2、hanning自乘法窗函数频谱特性a3.计算相关谱线系数，参见表3。表3、谱线的加权系数窗类型c1c2c3c4c51阶0.632910002阶0.41400.82721003阶0.30060.63290.9004104阶0.23440.50500.75710.93561由于4阶hanning乘法窗函数主主瓣内有10根谱线，两次谐波间隔中最多也为10根。a4.根据步骤s4的推导，表2中的4种窗函数的修正公式α＝g-1(β)，u(α)分别如下：(1)1阶hanning自乘法窗函数：α＝1.1301β-0.1836β3+0.0811β5u(α)＝0.5355+0.1758α2+0.1043α4如图3所示，为1阶hanning自乘法窗函数幅频特性。(2)2阶hanning自乘法窗函数：α＝1.5438β-0.3134β3+0.1618β5u(α)＝0.7543+0.2249α2+0.1545α4如图4所示，为2阶hanning自乘法窗函数幅频特性。(3)3阶hanning自乘法窗函数：α＝1.8496β-0.4045β3+0.2181β5u(α)＝0.5745+0.1716α2+0.1053α4如图5所示，为3阶hanning自乘法窗函数的幅频特性。(4)4阶hanning自乘法窗函数：α＝2.1895β-0.5473β3+0.2937β5u(α)＝0.4668+0.1854α2+0.1512α4如图6所示，为4阶hanning自乘法窗函数的幅频特性。a5.计算基波参数，公式如下：f0＝k·δf＝(k2+α+0.5)δf；a6.确定基波频率f0后，在范围(kf0-5,kf0+5)内重复步骤a2～a5，直到所有谐波参数计算完毕。a7.进行误差分析。仿真实验的误差分析比较结果参见表3～表6所示，其中，dai表示基波和各次谐波幅值测量值的相对误差；df0表示基波频率测量值的相对误差；表示基波和各次谐波初始相位测量值的相对误差，均用百分比表示。表4、hanning各阶乘法窗的基波频率相对误差阶次1234df0-1.97e-71.43e-82.06e-9-2.55e-10表5、hanning窗各阶幅值相对误差比较项1阶2阶3阶4阶da1-1.23e-5-2.37e-6-7.65e-7-7.35e-8da2-7.96e-5-6.25e-5-1.96e-5-7.53e-6da3-4.53e-6-2.12e-6-6.49e-7-3.66e-7da46.59e-54.63e-62.51e-6-3.66e-7da52.58e-5-3.69e-51.38e-69.87e-7da67.75e-53.45e-62.42e-69.57e-7da7-2.61e-53.80e-61.17e-6-7.85e-7da8-7.37e-5-4.25e-6-2.75e-69.94e-7da9-5.55e-5-3.52e-6-1.84e-6-1.00e-6da10-1.10e-4-3.53e-5-1.13e-5-1.09e-5表6、hanning窗各阶相位相对误差从表4～表6数据可以看出，本发明实施例所推导的多谱线插值算法具有良好的精度。如图7所示，本发明基于主瓣宽度的多谱线插值谐波分析系统，包括信号预处理单元、谱线数目确定单元、权值确定单元、修正公式计算单元、基波参数计算单元和谐波参数确定单元。所述信号预处理单元，用于对上位机中来自互感器的离散的电网信号x(n)进行预处理；n为采样点的序数，n为自然数；电网信号包括电流信号、电压信号。采用离散余弦窗函数w(n)对电网信号x(n)进行加窗截断，得到加窗信号xw(n)：xw(n)＝x(n)w(n)式(1)离散余弦窗函数w(n)的表达式为：其中，n为采样点数，n为正整数，n＝0,1,2...n-1；σ表示求和；m为窗函数的累加次数，m＝0,1,2...m-1；m为窗函数项数，m为正整数；bm为窗函数系数。窗函数系数bm满足约束条件如下：对式(1)的加窗信号xw(n)进行fft变换后，得到加窗fft频谱：其中，w(·)为窗函数的频谱，k为正整数，x(k)表示第k次谐波的频谱，ak为第k次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，第一次谐波为基波，fs为采样频率，f0为基波频率，δf为离散频率间隔，且δf＝fs/n。令：k0＝f0/δf，k0为真实频谱的谱线位置，真实频谱的谱线位置k0为小数，忽略负频率点处旁瓣的影响，式(3)变为：所述谱线数目确定单元，用于确定多根谱线。具体的，对于整数次谐波之间间隔为50hz，那么非同步采样带来的泄露频谱线在谐波间隔中最多有h1根，h1＝int(50/δf)。若所加窗函数的主瓣宽度为δωrad/s，主瓣内的谱线根数h2为h2＝int(δω/2πδf)，多谱线的根数h确定为h＝min(h1,h2)。所述权值确定单元，用于计算多谱线所对应得权值大小；具体的，令对信号的非同步采样或非整周期数据截断时，由于栅栏效应，信号的峰值频率f0＝k0·δf很难刚好位于离散谱线的频点上，即k0一般不为整数。设峰值频率点左右多条的谱线分别为k1＜…＜kh＜k＜kh+1＜…＜kh，对应的谱线幅值分别为y1，…，yh，yh+1，…，yh。记α＝k-kh-0.5，由于0≤k-kh≤1，则-0.5≤α≤0.5，另记：其中α和β均为系数，ci为加权系数，i＝1…h，ci为yi的权值。考虑到计算的方便，令ci＝c(h-i+1)，余弦窗函数的表达式如式(2)所示时，根据频域主瓣表达式，加余弦窗时候，可以求得：所述修正公式计算单元，用于计算多谱线插值算法的修正公式。具体的，根据式(4)和式(5)得到其中，ci为加权系数，i＝1…h，ci为yi和ym-i的权值。当n大于1000时，式(7)可以简化为β＝g(α)，其反函数为α＝g-1(β)；由于所采用的余弦窗系数均为实系数，其频率响应是偶对称的，因而g(·)和g-1(·)均为奇函数。采用多项逼近方法计算奇函数β＝g-1(α)，表达式为：α≈p11×β+p13×β3+…+p1pβp式(8)其中，p11，p13，…，p1p为多项式逼近的奇次项系数，由式(4)可知信号幅值ai为：可以得到：通过类似于式(10)的逼近方法，当n大于1000时，窗函数系数为实系数，式(10)可表示为：u(·)为偶函数，逼近多项式不含奇次项，则多谱线修正逼近多项式为：u(α)＝(p20+p22α2+…+p2dαd)式(11)其中，p20，p22，…，p2d为多项式逼近的偶次项系数。所述基波参数计算单元，用于计算基波参数。具体的，包括计算基波频率f0、基波幅值a1：f0＝k·δf＝(k2+α+0.5)δf式(12)根据式(4)得出基波的相位所述谐波参数确定单元，用于确定谐波参数。具体的，确定基波频率f0后，在范围(kf0-5,kf0+5)内重复步骤s2～s5，直到所有谐波参数计算完毕。本发明还可以设置一个误差分析单元，在窗函数主瓣宽度不断增加，参与运算的谱线数目也相应增加的情况下，进行多谱线插值计算的误差分析，验证本发明分析的精确度。本发明不局限于上述实施方式，对于本
技术领域：
的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围之内。本说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。当前第1页12