魏 国,张 蓓,撒文彬,孙金玮

(哈尔滨工业大学电气工程学院,哈尔滨 150001)

电力系统无功功率是电力系统中一个非常重要的技术指标。控制无功的基础是无功功率测量的实时性、有效性和准确性。基于C.Budeanu对无功功率的定义[1],目前主要的无功功率测量方法有移相法[2]、平均功率瞬时检测法[3]、基于快速傅里叶变换FFT(fast Fourier transform)的方法[4],还有基于H ilbert变换[5,6]和基于H ilbert数字滤波器[7～10]的方法等。其中Hilbert变换法较其他方法精度更高,能更精确地满足无功功率的定义,是应用较广泛的一种方法。

传统离散时间Hilbert变换是基于离散傅里叶变换DFT(discrete Fourier transform)和离散傅里叶反变换IDFT(inverse discrete Fourier transform)实现的。根据这种方法,在实际应用中只需离线求出一个Hilbert变换矩阵H,然后用采样序列与H相乘即可得到移相90°后的信号。经实际验证,该算法具有较高的精度,在信号的采样点数少的情况下,测量较为方便,计算时间也能基本满足实际工程需求。但进一步测试发现,若采样点数增多,H矩阵的维数迅速增大,信号的Hilbert变换需要消耗大量时间,无功测量的实时性和有效性已经不能满足实际要求。经分析,若采用加窗FFT和IFFT方法实现信号Hilbert变换,可有效缩短计算时间,提高无功测量的实时性,其变换效果也优于传统方法。

1 Hilber t变换原理

理想的H ilbert变换是幅频特性为1、负频率成分进行90°相移、正频成分进行-90°相移的线性变换,其频域传递函数为

设有实信号f(t),其傅氏变换为F(ω)。令R1(ω)、I1(ω)分别表示ω＞0时的频谱实部与虚部,有:

令f′(t)=f(t)+j(t),则f′(t)的傅里叶变换为

由式(4)可知,f′(t)的正频分量为f(t)的两倍,负频分量为0。实际应用中,可先求出信号f(t)的傅氏变换F(ω),然后根据式(4)求得F′(ω),进而得到F′(ω)对应的时域信号f′(t)。f′(t)的虚部即为经过H ilbert变换后的信号(t)。

2 传统的离散时间Hilbert变换方法

以采样频率f s对电压信号采样m个周期,每周期采样M点,得N点序列u(n)。令u(n)基波频率真值为f1,所含最高谐波次数为h,传统的离散时间H ilbert变换方法如下。

1)对u(n)进行DFT运算,得到信号各次谐波分量。

式中U(0),U(m),…,U(mh)等效为序列u(n)的离散时间傅里叶变换DTFT(discrete time Fourier transform)在[0,f s]区间上以m×为间隔的h+1点采样,分别表示信号的直流分量、基波、…、h次谐波分量。

2)对各次谐波分量加权,然后进行IDFT运算。

3)综合前两步,序列u(n)的Hilbert变换过程等效于H矩阵与u(n)乘积的过程,即u^(n)=H u(n),其中H为N×N阶实矩阵。

传统方法简单直观,易于软件实现,在采样点数较少时不失为一种好方法。但信号采样点数增多后,H矩阵维数变大导致计算量快速增大,降低了无功测量的实时性。下文对传统方法进行了改进。

3 加窗FFT和IFFT实现Hilbert变换

3.1 用FFT替代传统方法中DFT运算

对于长度为N的序列x(n),其N点DFT为

令序列长度为N=2M,对x(n)定义另一种形式的DFT运算:

则有

考虑到

且N=2M,式(10)变为:

若x(n)长度为N=m×M,参照上述推导有

就等同于

容易看出,Xs(k)序列等效为x(n)的DTFT在区间[0上以m×为间隔的M点采样。

参照上述推导,对于式(5)中的电压序列u(n),若将其m个周期上的采样点对应累加到一个周期并做M点FFT,所得M点序列Us(k)即等效为u(n)的DTFT在[0,f s]上以f为间隔的M点采样。其中,显然,U s(k)与式(5)中U(mk)的意义完全相同。由此,式(5)可变为:

求取U(0)～U(mh)的过程简化为:(1)将N点序列u(n)累加到一个周期,得M点序列us(n);(2)对u s(n)做M点FFT,得U s(k),k=0,1,…,M-1。Us(0),Us(1),…,Us(h)即为U(0),U(m),…,U(mh),分别表示信号的直流分量、基波、…、h次谐波分量。

3.2 用IFFT替代传统方法中IDFT运算

在式(6)中,令:

对B做如下处理。

1)在B的每两个元素间补m-1个零,在最后一个元素后补N-1-mh个零,得:

2)在B后补上M-h-1个零,得:

写出N点IFFT的标准系数矩阵A1。

易得:

可见,式(6)的运算过程等效于对B1的N点IFFT所得时域序列b1(n)取虚部的过程,即

由于B1中部分元素为0,b1(n)可写为

易得:

参照式(22)、(24)和(26),式(6)的运算最终等效为:1)对B2做M点IFFT运算并将所得M点时域序列除以m,得序列b1(n)(n=0,1,…,M-1);2)根据采样周期数m延拓b1(n)(n=0,1,…,M-1),得b1(n)(n=0,1,…,N-1);3)对N点长度的b1(n)序列取其虚部即得到电压移相信号(n)。

3.3 窗函数应用于Hilbert变换

传统的H ilbert变换方法需先求出信号各次谐波参数。由于电网信号的频率通常会在额定频率附近波动,同步采样不可能实现。对应式(14)有:

采样的非同步性导致信号FFT分析时各谐波间出现频谱泄漏。同时FFT算法本身又会引入栅栏效应,这两因素都将降低谐波参数的测量精度。对应式(15)中,U(0)～U(mh)的值不再准确。

若选择合适的窗函数则可有效地抑制频谱泄漏,同时,栅栏效应可以通过基于窗函数的插值算法加以克服。但插值算法要求精确求解信号采样同步误差,其计算量通常较大,消耗时间多,影响无功功率测量的实时性。实际应用中,只需根据信号的采样周期数,选择合适的窗函数进行运算,所得谐波参数测量结果就足以满足精度要求。目前比较常用的窗函数为组合余弦窗函数,如Hanning窗、Blackman窗、Blackman-H arris窗、矩形自卷积窗[6]等。

设长度为N的窗序列为:

在H ilbert变换中,对信号u(n)加窗得序列:

综上所述,只需对长度为m×M的周期离散序列加窗,然后进行1次M点的FFT和1次M点的IFFT运算即可准确实现其H ilbert变换。与传统的离散Hilbert变换方法相比,改进方法有效简化了计算,并且对信号进行了预加窗处理,改善了变换的效果。称上述算法为基于加窗FFT和IFFT的H ilbert变换方法。

4 计算量比较及验证

令1次实数乘法所需运算时间为a,1次实数加法运算时间为b,1次复数乘法运算时间为a1,1次复数加法运算时间为b1。忽略加窗所需要的计算时间。

但进一步分析发现,式(7)中的H矩阵具有一定的对称特性,具体来说,H可表示为

式中,H′是H的子矩阵,为M阶方阵。

若将电压序列u(n)写成如下形式:

其中,uMi表示u(n)中的第i个周期的采样点构成的序列,为M长度的列向量。

可见,使用传统方法对u(n)移相的过程最终可简化成使用H′矩阵与u s(n)做乘积然后将所得时域序列周期延拓的过程。实际应用中,只需离线求出H′矩阵,即可实现信号的Hilbert变换,这样可以提高计算速度、节省大量的存储空间。

考虑到应用传统方法时可以做如上调整,其计算时间缩短为

下文提到的传统算法都是指调整后的传统算法。

改进算法不需要离线计算,但需对信号进行1次M点FFT和1次M点IFFT,其计算时间约为

1次复乘相当于4次实乘和2次实加,1次复加相当于2次实加,则

式(35)变为

参考(34)、(37),传统算法与改进算法计算时间之差为

若只考虑乘法运算:

可见,ΔT只与单周期采样点数M有关。

用ΔT/T′o表示改进算法相对于传统方法节省计算量的相对值,表1给出了其与M之间的关系(单位:1)。

表1 ΔT/-M关系表Tab.1 Relationship betweenΔT/ and M

表1 ΔT/-M关系表Tab.1 Relationship betweenΔT/ and M

642560 a 4096 a 62.512812800 a 16384 a 78.125657344 a 65536 a 87.5512243712 a 262144 a 93.010241007616 a 1048576 a 96.120484104192 a 4194304 a 97.9

由表1知,随着M的增加,ΔT/T′o的值快速增大,改进算法的速度优势更加明显。

改进算法已应用于“三相数字式多功能测控电表”的无功功率测量中,该电表使用自带12位A/D转换器的C8051F系列单片机。实验表明,改进算法相比传统算法,无功功率计算速度快出很多且测量精度亦有所提高。单周期采样128点,采样2周波时,传统算法实现Hilbert变换所需时间为89 m s,而改进算法只需要42 ms,计算时间减少了约52.8%。虽然根据理论推导,在M=128时,改进算法节省时间量应为78.1%,但上述推导只考虑了FFT蝶形运算所需时间,而程序具体实现中需对采样信号加窗,并进行码位倒置和寄存器赋值等一系列额外运算,所以说该实验结果与理论推导保持一致。

5 结语

本文对基于DFT和IDFT实现H ilbert变换的传统算法进行改进,提出了一种基于加窗FFT和IFFT的H ilbert变换快速算法。该算法只需对加窗信号做1次M点FFT运算和1次M点IFFT运算(M为单周波采样点数)即可实现采样序列的准确90°相移。理论分析表明,采用此算法能有效缩短信号移相计算时间,提高无功测量的精度。改进算法已经应用于工程实践,在应用中提高了程序运行的效率和无功测量的实时性、准确性。

[1] 郑小平(Zheng Xiaoping).关于无功功率的定义及其计算方法(The definition and arithmetic for the reactive pow er)[J].电测与仪表(ElectricalM easurement&Instrumentation),2006,43(6):1-4,16.

[2] 陈国通,吴杰康,张宏亮(Chen Guotong,Wu Jiekang,Zhang Hongliang).无功功率和电能的移相算法(Phase-shift algorithms for reactive electric power and energy in pow er system s)[J].电力学报(Journal of Electric Pow er),2007,22(4):428-431.

[3] 薛蕙,杨仁刚(Xue Hui,Yang Rengang).改进的瞬时无功和谐波电流检测理论(A novel detection theory of instaneous reactive and harmonic current)[J].电力系统及其自动化学报(Proceedings of the CSUEPSA),2002,14(2):8-11.

[4] 邱海锋,周浩(Q iu H aifeng,Zhou H ao).电力系统无功测量方法综述(Summary on reactive pow er measurement in pow er system)[J].电测与仪表(Electrical M easurement&Instrumentation),2007,44(1):5-9.

[5] 郑常宝,王群京,方斌,等(Zheng Changbao,W ang Qun jing,Fang Bin,et a l).用小波变换和希尔伯特变换测量无功功率(Reactive power measurement by wavelet transform and H ilbert transform)[J].系统仿真学报(Journalof System Simulation),2005,17(4):822-824.

[6] 黄纯,江亚群(Huang Chun,Jiang Yaqun).电功率微机测量新算法(Novel algorithm for measuring active and reactive pow er)[J].中国电机工程学报(Proceedings o f the CSEE),2008,28(4):54-58.

[7] 俎云霄,庞浩,李东霞,等(Zu Yunxiao,Pang H ao,Li Dongxia,et a l).一种基于H ilbert数字滤波的无功功率测量方法(A method of reactive power measurement based on H ilbert digital filtering)[J].电力系统自动化(Automation of Electric Pow er System s),2003,27(16):50-52,70.

[8] 庞浩,王赞基,陈建业,等(Pang Hao,Wang Zan ji,Chen Jianye,etal).基于2对H ilbert移相滤波器的无功功率测量方法(Method of reactive pow er measurementbased on tw o pairs of H ilbert digital filters)[J].电力系统自动化(Automation of Electric Pow er Systems),2006,30(18):45-48.

[9] AnsariRashid.IIR discrete-time H ilbert transformers[J].IEEE Trans on Acoustics,Speech,and Signal Processing,1987,35(8):1116-1119.

[10]Ansari Rashid.Ellip tic filter design for a class of generalized halfband filters[J].IEEE T rans on A-coustics,Speech,and Signal Processing,1985,33(5):1146-1150.