基于kaiser窗的电力系统谐波和间谐波参数估计

基于Kaiser窗的插值FFT谐波与间谐波的参数估计摘要：随着精密仪器在电网中的应用日益增多，由间谐波引起的如感应电动机的振动与噪声、电压闪变等问题也日益突出。间谐波的精确检测是解决这些问题的前提，因此研究间谐波的检测方法具有重要的现实意义。基于凯泽窗优良的窗函数特性，提出了加凯泽窗的分段逼近参数估计方法，分析了该方法的可行性和优点，在时域与频域对凯泽窗和其他窗函数的间谐波检测性能进行了详细比较，并通过仿真验证了凯泽窗的设计灵活性，以及该方法的正确性和有效性。关键词：间谐波；加窗插值；频谱泄露；栅栏效应；凯泽窗TheDetectionandRealizationofPower-harmonicsand

interharmonicsbasedonLabviewABSTRACT：Withthewidespreadofsophisticatedequipmentsinpowersystem,problemsthatareassociatedwithinterharmonics,suchasthenoiseandvibrationofinduction-motor,voltageflickerandsoon,hasbeenincreasinglyprominent.Todetectinterharmonicsprecislyistheprerequisiteofsolvingsuchproblems,therefore,theresearchoninterharmonicsdetectionisofpracticalsignificance.BasedontheexcellentwindowfunctionperformanceoftheKaiserwindow,themethodofsegmentationapproximationisproposedtoestimatetheparametersofinterharmonics,andthefeasibilityandmeritsofthismethodisanalyzed.Inaddition,thecomparisonofinterharmonicsdetectionperformancesbetweentheKaiserwindowandotherwindowedinterpolationalgorithmsinboththetimeandfrequencydomainaregiventoveritytheflexibilityoftheKaiserwindowmethod.Thevalidityofthemethodisprovedthroughsimulationresults,andtheirprecisionininterharmonicdetectionarealsoanalysizedaccordingtothecorrespondingnationalstandard.Keywords:Interharmonics;WindowedInterpolationalgorithm;frequencyspectrumleakage;picketfenceeffects;Kaiserwindow0引言随着经济和社会的发展，各种基于电力电子技术的波动性负荷如变频调速设备、整流供电装置、循环变流器件等的应用日益广泛，这些设备运行时电流中包含大量的间谐波。这些间谐波流入电力系统后会引起白炽灯管或CTR射线管等显示设备闪烁，干扰电网中的低频控制信号，引起音频设备和感应电动机噪声增大，

甚至可能导致滤波器因过流跳闸和滤波失败等问题。 目前间谐波已成为继谐波之后，国际上公认的公用电网污染的主要公害之一⑴。在对间谐波的研究过程中，实现精确可靠的间谐波检测是对间谐波产生机理和传播规律进行理论分析的前提，也能为研究其危害并制定相应抑制措施提供依据，因此研究间谐波的检测方法并实现间谐波的可靠检测具有重要的理论和应用价值。1凯泽窗与基于凯泽窗的间谐波检测算法1.1凯泽窗（Kaiser）人们研究窗函数主要是它可以有效地抑制频谱泄漏，本论文中采用的加窗插值FFT方法，窗函数的选择非常重要。在间谐波检测的频谱分析时，最好的是窗函数主瓣窄、旁瓣低且跌落快，但对于同一窗函数，以上两个要求又是互相矛盾的。因为要是增加主瓣的宽度，旁瓣就会降低，反之，若想主瓣能变得又高又窄，旁瓣就会增高。在选择窗的时候，应根据间谐波特征和研究目的来选择。目前，已经有20多种窗函数，在电力系统谐波与间谐波检测中常用的窗函数有矩形窗、海宁窗（Hanning）、汉明窗（Hamming和布来克曼窗（Blackman）。数字信号处理领域较早就提出了凯泽窗（Kaiser），它通过改变参数可以达到不同的性能，正由于其优良的窗函数特性，广泛应用于高通、低通、带通、带阻等各种滤波器的设计。凯泽窗（Kaiser）是利用贝塞尔函数来逼近需要的理想窗，其时域函数形式如式3-1w(n)二也；1-[1-2形式如式3-1w(n)二也；1-[1-2n/(N-1)]2)

l°(B)n=0,1,2N-1式3-1式中，1。「）为零阶第一类修正的贝赛尔函数，可用式3-2的级数表示。COb(X)=COb(X)=1 '[n=±(x/2)n

n!]2式3-2凯泽窗（Kaiser）的幅度公式为:（N」）/2W（）=w（0） 2'w（n）cos‘n 式3-3n二由公式可知，凯泽窗的B值与贝赛尔函数级数项数n是两个独立的参数，但都会对凯泽窗造成影响，为进一步的了解和研究凯泽窗，有必要分别对其分别介绍。1.1.1贝赛尔函数的项数（n）确定凯泽窗（Kaiser）参数B的前提下，贝塞尔函数的项数（n）设定对凯泽窗的影响就变得非常重要了。贝塞尔函数的项数（n）对凯泽窗的影响本质上就是利用项数的增大来近似逼近理想的凯泽窗（Kaiser）。n的取值决定了需要的凯泽窗的精确程度，通常用 15~25有限项去近似表达这个无穷级数[11][12][13]。matlab自带的kaiser函数已经满足本文两种加窗插值算法的精度要求， 下图是matlab自带kaiser函数与项数设置为43的凯泽窗函数时域比较图（B固定为20，N取1024）。图3-1自带kaiser函数与项数设置为43的凯泽窗函数时域比较利用matlab的科学研究中完全可以采用其自带的凯泽窗函数“ kaiser（N，beta）”，但这并不妨碍对贝塞尔函数的项数（n）的详细研究，因为间谐波检测的工业应用千差万别，特定情况下可能对准确度很敏感，某些情况下又可能对计算量很在乎。图3-2为不同的项数（n）所对应的凯泽窗（Kaiser）时域比较示意图，以及部分n值与matlab自带kaiser函数进行幅度特性比较图（采样点数为64，B固定为20）。

图3-2时域与频域比较贝塞尔函数项数（n）的数值在低于15的时候，凯泽窗的时域变动比较明显。超过15之后一直到70，图形几乎重叠，图形上已经无法分辨各个曲线，实际上从表3-1的详细数据上也能看出同样的变化规律：当项数大于 15之后自定义的凯泽窗（Kaiser）变化很微小。Kaiser）详细数据（采样点数为64，B为20)表3-1不同项数的凯泽窗（i贝塞尔函数项数泄露系数（%旁瓣相对衰减（dB）主瓣-3dB带宽（rad）matlab自带0-154.90.074219n=700-154.90.074219n=350-154.90.074219n=150-99.80.074219n=100-67.30.070313n=60-48.20.058594n=50-43.40.054688n值为15以下如n等于10、6、5的时候，其凯泽窗主瓣、旁瓣变动较大，而当n值为15以上即n值为35、70以及matlab自带的kaiser函数时，其凯泽窗主瓣、旁瓣几乎不变。如果对时实性与计算量的要求不是很高的情况下，为了实现程序的简洁与契合度，算法的仿真可直接使用matlab内部kaiser函数，当然如有其他的设计要求，可灵活设置该参数。1.1.2B参数（beta）上节3.1.1已经详细的阐述了贝赛尔函数级数的项数变化对凯泽窗（Kaiser）的影响，本节中将详细阐述B值大小对凯泽窗（Kaiser）的影响（确定贝塞尔函数级数的项数）。随着B增大，主瓣加宽，旁瓣幅度减少［12］。图3-3给出了不同1值所对应的凯泽窗（Kaiser）时域和频率图形（采样点取 1024,利用matlab自带贝赛尔函数）凯湮奮（Kji&ar）京时域图那来样点n利用matlab自带贝赛尔函数）凯湮奮（Kji&ar）京时域图那来样点n图3-3凯泽窗（Kaiser）时域和频域图形时域图中，从上往下,1分别为0，1.231，2.341，3.440，4.394，5.440，分别摘取1为0、4.394、6.451，7.321，分别摘取1为0、4.394、8.500、10.399来进行幅度特性对比， 明显看出1为一个可自由选择的凯泽窗函数形状参数，能调节主瓣和旁瓣的宽度［12］。利用MATLA的求得不同的1值的凯泽窗（Kaiser）详细数据（采样点数定为1024），总结出表3-2表3-2不同B值凯泽窗（Kaiser）的详细数据（采样点数为 1024）Kaiser窗函数泄露系数（%旁瓣相对衰减（dB）主瓣-3dB带宽（rad）B=09.15-13.30.001709B=1.2315.48-15.40.001709B=2.3411.65-20.10.0019531B=3.4400.33-26.40.0021973B=4.3940.07-32.60.0021973B=5.4400.01-39.80.0024414B=6.4510-47.10.0026855B=7.3210-53.50.0029297B=8.5000-62.40.0031738B=9.4120-69.50.0031738B=10.3990-77.30.003418根据式3-1凯泽窗(Kaiser)的公式，当一：=0时，Kaiser窗和矩形窗函数是一样的图形。当1=3.440时，Kaiser窗旁瓣相对衰减速率为-26.4dB，继续增加］值，旁瓣峰值电平将继续下降，渐进衰减速率继续增加。1值越大，其主瓣宽度也越大，频谱的旁瓣也随着变小。利用matlab和本文分段逼近的加窗插值算法(见415本文3.4.2)，对信号f(t)=247sin(2nt一+13°)进行间谐波检测(采样频率5120为5120Hz采样点数为1024)，得到结果：表3-3不同卩值检测幅值表测量值(B=18)实际值测量值(B=5.44)测量值(B=6.121)测量值(B=7.865)频率41.54040.000340.137941.4017幅值247276.2669274.2383267.0376247.5005初始角1366.972866.960261.990316.5354曲5.44、6.121、7.865分别可以近似取代汉明窗、海宁窗、布来克曼窗(见3.2)，结果证明在单一间谐波的情况下凯泽窗可以调节 竝以提高主瓣的能量比重，达到比其他常用窗更好的测试精度。同样的采样频率和采样点数对频率50Hz、幅值200V、初始相位200的基波和频率151Hz、幅值0.6V、初始相位10°的间谐波的叠加电压信号进行幅值检测，见表3-3。实际幅值表3-4不同卩值检测幅值表测量幅值(B=5)测量幅值(B=10)测量幅值(B=15)测量幅值(B=20)基波(50Hz)200v200.1118v200.0706v200.0518v199.9592V间谐波(151HZ0.6v0.6540v0.6357v0.6103v0.6000V当频率差比较大的时候(即主瓣影响可以忽略)，两信号分量的幅值的变化只与两者旁瓣影响有关，旁瓣跌落厉害则幅值测试准确，旁瓣跌落很小则幅值检测相对就不那么准确，故上表随着B值的增加，旁瓣衰减变大，幅值测试准确度明显提高。间谐波检测平台完全可以利用凯泽窗可变动主瓣与旁瓣的特点来提高在特定情况下的测试精度。2基于凯泽窗的分段逼近算法实际电网中，电力信号x(t)是含有基波分量、谐波分量、间谐波分量的复杂信号。以fs采样频率对它进行均匀采样离散化，得到x(t)的离散时间信号

x(n)八sin(2二x(n)八sin(2二n7zkfl

fS-;:k)式3-13式中M――间谐波与谐波的总项数；Ak――k项间谐波/谐波的幅值；fi――基波的频率；fs――采样频率；\――k项间谐波/谐波的初相位；Zk――k项间谐波/谐波相比fi的倍数。特别要指出：由于信号中含有基波非整数倍的间谐波分量，Zk不一定是正整数。设凯泽窗(Kaiser)的窗函数为w(n),利用w(n)对信号x(n)进行加权截断，得到加窗后的离散信号xw(n)=x(n)w(n)，再经过离散傅里叶变换(DFT)得到的表达式:Xw(h：Xw(h：f)「AkejkW产y2j. fs)_e」kW(2：(hfZkf1))fs式3-14式3-14中W(x)为凯泽窗(Kaiser)的频谱函数表达式:W(x)二式3-15Nsinh(「：2-(Nx/2)W(x)二式3-151。(：)「2-(Nx/2)2凯泽窗(Kaiser)窗函数的优良特性使负频点yh的旁瓣对Xw(hf)的影响非常微弱，可以忽略，当然这是在我们已经选定参数的凯泽窗( Kaiser)来截断的前提下，如果采用其他的窗函数或其他方式截断使频率分辨点离得过近， 则忽略负频点的影响将带来较大误差。利用凯泽窗(Kaiser)对信号进行离散抽样后，忽略负频点处频峰的旁瓣影响，表达式简化为：Xw(hf)心Xw(hf)心22二(hf-ZkfJfs式3-16其中・讦为离散抽样间隔fs/N，再设Zkf,=hk•讦，由于Zk不一定是正整数,故这里hk也不一定是正整数。转换式3-16为：Xw(h'f)讣Xw(h'f)讣kW(竿严式3-17只要推出检测其中一项间谐波参数的方法，检测其他项间谐波参数的方法类似，即它们都有一般性。设需要检测的是第i项间谐波，简化计算量和复杂度，忽略掉基波，其他谐波、间谐波对其的频谱泄露影响，此时，表达式为：Xw(hf)-hXw(hf)-hi)2j式3-18需要注意的是上面的估算是在一定的前提条件下成立的， 其不仅假定正负频率分量之间是无干扰的，而且假定各频率间也是无干扰的。但是由于主瓣能量较大，当间谐波频率离基波较近时，其不仅会受到基波旁瓣频谱的泄露影响，而且当主瓣宽度大于二者频率间隔的时候，间谐波也将受到基波主瓣的影响，甚至会出现频谱淹没现象导致无法检测出来，这也是加窗插值算法检测间谐波的技术难题。假定i项间谐波的峰值频点左右分别有幅值最大和幅值次最大的第 h1和h2条谱线。实际情况中h1，h2谱线对应的幅值不一定刚好是最大和次最大， 还可能是次最大与最大的情况，不过两种情况进行参数估计的原理和方法是一致的， 为了简化理解，我们假定hi，h2的幅值分别为最大和次最大。其中hi：：：h「：h^hi1，这两条谱线幅值分别是Y=Xw(hQf)|和丫2=Xw(h^f)o由于函数为超越函数，无法由确定两点的谱值数，得到准确的实际频率分辨点，只能利用插值估计方法来逼近其真实值，故引入参数：二h-0.5，数值：的范围在卜0.5,0.5］，分段逼近算法引入的是参数Xr，其表达式：式3-31式中丫1、丫2分别是我们假定求参数的i项间谐波峰值频点左右最大h1谱线的幅值与次最大h2谱线的幅值。可以根据丫1、Y2的定义以及与式3-20同样的变量替换，可将式3-31转换为：2二-' 0.5W[式3-32X式3-322二_0.5W[由此可见，上式为一个关于S的函数frCJ且满足Xr二frC)o不同的S值

对应不同的fr（J，我们给定一组n个S值，假定为乩3,&……Sn,得到n个幅值比fr（"）,fr（、・2）,fr（、3）……fr（n）。根据fr（、i）和、i我们可以构造分段线性函数，如图3-11所示。其中和fr5,fr52是已知的测试点值。图3-11图3-11分段线性函数这样当我们已知第i项间谐波的双峰幅值比仆*时，便可以按照上图中的线性关系即直线斜率相等原理，得到表达式：式3-33§h2—国h1 _ §hi一§式3-33fr(「h2)—frGh1)「frChi)—讥廿)如此可用苛、十、frCh2）、fr（F）、fr、；hi来表示「hi'f._f(h;-f(h1'f._f(h;-f(h1)h1fr('hi)-frGh1)1式3-34然后将-hi代入matlab函数polyfit推导出系数后的g（S,再根据频率、幅值、初始相位的修正公式求出各个参数。式3-19I一丫2-丫1

丫2丫式3-19将式3-18代入，并用：变量替换可推得公式:

式3-20W[2n(」+0.5)/N]|—|W[2ji(」式3-20W[2二(—、.0.5)/N] |W[2「:(」-0.5)/N]引入参数厶使-=，同时为了让凯泽窗（Kaiser）的频谱函数满足［0,N-1］范围，将频域信号平移N-1］范围，将频域信号平移宁，转换后的频谱函数表达式:式3-21sinh(,M)2-[(N-1)x/2] 2)—j字式3-211。(二),e)21。(二),e)2-[(N-1)x/2]空丄，和上文的Z和hk情况类似，h不一定是整数，将其代入式3-21N得:W(2nh] N—1sinh(J(M)2-[(N-1)Jih/N]2)jj罟右N1。(得:W(2nh] N—1sinh(J(M)2-[(N-1)Jih/N]2)jj罟右N1。(二).(二)-[(n-1)二h/N]2式3-22如h"工」'「0.5，此时-0.5乞1，且N很大，N-1可近似为则式3-22的表达式转换为:N-1sinh(r2-(-、0.5)2)怙(二)二「2-(一、 0.5)2式3-23如h--「0.5，此时--0.5叮，同上可将3-22的表达式转换为:w(2l^052)NN-1sinh(r2 -0.5)2)1。(二)二i2-、-0.5)2式3-24将式3-23和式3-24代入式3-20，并且由于［和••都是常数，则可得到J二f（「J，其中f（「•）是关于的多项式。微型计算机不能处理连续的模拟信号，只能处理离散的数字信号，这就决定了在实际应用中「•并不是一个连续的值，它所对应的函数值f「）也一样不是一个连续的值。本论文作者将步长参数在程序中设置为“step”并将其定义为0.001，即在■■:的取值范围［-0.5Q.5］均匀取1001个离散点对应1001个离散的」值，然后推求一个解析函数f'（"）*（」）使其通过或近似通过有限序列的资料点），通过多项式函数逼近的方法求得此拟合函数，即所谓的曲线拟合（curvefitting）。由于曲线拟合（curvefitting）要选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系，所以必须要确定下自变量的阶数，根据实际应用需要和减少计算量的综合考虑，本论文作者选择5阶取得了较好的效果。MATLAB中polyfit函数很方便实现上述功能，具体程序为“P=polyfit(mu,delta,5)”其中mu表示丄,delta表示：，数字5代表多项式最高阶为5阶，P为曲线拟合后每阶常系数值。h1和h2两根谱线加权平均得到第i项间谐波/谐波的幅值修正公式：A2W严穿'A2W严穿'+A⑷严1"）]N■|W［"（h2—hi）】+W［2沢（h1—hi）］NW严3+0.5）］+27r（」_0.5）］NN2(丫1丫2)式3-25式中A、A分别是谱线h2、hi的幅值，当N取值很大时，公式还可转换为:式3-26式3-25与式3-26相等，整理可得:v(、Jv(、J二W【2心十°5)]|十W［2JI（4°.5）］NN2N式3-27同上利用MATLAB平台中曲线拟合polyfit函数进行多项