一种基于相位调制的新型低频正弦频率计算方法

王越超，李 军

(广东电网有限责任公司电力科学研究院，广东广州510080)

0 引言

在系统频率测量和频率特性分析方面，各种算法不断产生［1－15］，在低频正弦频率测量上有多种算法［1－8］，如零交法［1］、基于幅值调制的测量算法［2］、基于带通滤波的算法［3］、基于小波变换算法［4］、基于人工智能的的算法［5］、基于离散傅里叶(DFT)变换算法［6］、基于一阶导数的算法［7］等。

高准确度的正弦频率测量技术是一些应用技术的基础，文献［1］指出了一种电力系统阻抗测量方法需要有精准的频率测量结果作为参考值。电力系统额定工频50 Hz［9］属于频率较低的正弦频率。但目前的频率测量技术在低频正弦频率测量方面普遍存在准确度不高、抗噪声包括抗谐波干扰性不强等问题。

零交法(zero-crossing algorithm)是低频正弦频率测量的基本方法［1］，如用于电力系统工频频率的测量，该方法原理是通过检测过程信号波形的过零点，利用1 个或几个周期过零点的时间间隔来计算出此段波形的频率值。但研究结果表明，在存在干扰的情况下，该方法测量出的频率值存在较大的误差［1］。

离散傅里叶(DFT)变换算法是用于正弦频率计算的基本数学方法［6］。但DFT 算法存在一定的局限性，其中信号的非整数周期截断所引起频谱泄漏问题是造成DFT 算法误差的主要内在原因［6］，由于频谱泄漏问题客观上难以避免，因此DFT 算法难以满足正弦频率的高准确度计算要求。

在现有高准确度低频正弦频率测量方法研究中［8］，一些方法在研究上没有充分考虑实际信号处理技术的局限性和信号构成的复杂性，如信号离散采样产生的数据量化背景噪声影响问题、信号中的分次谐波影响问题等。

1 新型正弦频率计算原理

新型正弦频率测量方法的基本原理是:对低频正弦信号序列进行相位调制，得到相位随信号频率变化、但相位变化方向相反的2 路相位调制序列;通过精密相位计算，得到2 路相位调制序列的精密相位;通过精密相位(差)计算，得到2 路相位调制序列的精密相位差;根据2 路相位调制序列的精密相位差，可得到信号序列的正弦频率，如图1所示。

图1 正弦频率测量原理示意图Fig.1 Schematic diagram of the sinusoidal frequency measurement

图1所示，正弦频率测量需要借助一个频率初测单元实现，频率初测单元作用是给出参考频率，并根据参考频率进行精密相位计算，允许频率初测单元存在± 0.25% 以内的相对误差。

2 相位调制原理

相位调制方式为:用信号序列和它的±1π 移相序列相减，得相位随信号频率变化、但相位变化方向相反的相位调制A 序列和相位调制B 序列。

首先令信号序列为

式中:Xi(n)为信号序列;A 为信号幅值，V;ω 为信号频率，rad/s;Tn为信号采样间隔，s;φ 为初相位，rad;Nπ为1π 移相序列长度，单位无量纲;N 为序列长度，单位无量纲。

根据参考频率计算1π 移相序列长度，为

式中:(int)代表取整数。由于Nπ整数化存在1 个采样间隔内的误差，需要根据实际1π 移相序列长度计算调相频率，为式(3):

式中:ωph为调相频率，rad/s。

信号频率与调相频率的频差为

式中:Ωph为调相频差，rad/s。

信号频率与参考频率的频差为

式中:Ω 为频差，rad/s。

相位调制A 序列为

式中:XPhA(n)为相位调制A 序列;kph为相位调制序列的幅值系数、单位无量纲。相位调制A 序列的调制相位为

式中:PHA为相位调制A 序列的调制相位，rad。

相位调制B 序列为式(8):

式中:XPhB(n)为相位调制B 序列。相位调制B 序列的调制相位为式(9):

式中:PHB为相位调制B 序列的调制相位，rad。

相位调制需要使用1 倍信号周期序列长度。

所述两种相位调制序列在频域的幅频特性完全相同，其中在参考频率100 πrad/s，得到的频域幅频特性，如图2所示。

图2 输出序列幅频特性示意图Fig.2 Schematic diagram of the amplitude frequency characteristics of the output sequence

图2所示，相位调制序列对偶次谐波有较好的抑制作用，对分次谐波也有较好的衰减作用。

3 精密相位计算原理

图1所示，相位计算准确度直接决定了新型频率测量方法的准确度。为了获得较高的相位计算准确度，文章提出了一种基于正交混频的精密相位计算方法，如图3所示。

图3所示，对输入序列进行正交混频，得到了互为正交的实频序列和虚频序列。

3.1 混频信号和混频干扰

图3所示，所谓的混频器实际上是乘法器，以实频混频信号为例，混频干扰分析如下:

设信号由基波、1/3、1/2、2，3 次谐波构成为例，为式(10):

图3 精密相位计算方法示意图Fig.3 Schematic diagram of the precise phase measurement

在信号基波频率等于参考频率，实频混频信号为

式中:XR(t)为实频混频信号;cos(φ)/2 为有用信号，其余的均为混频干扰信号。

以相位调制A 序列为例，在不考虑混频干扰前提下，实频序列和虚频序列为

式中:XRA(n)为实频序列;XIA(n)为虚频序列。

3.2 数字滤波

混频干扰严重影响相位计算准确度，是造成频率计算误差的主要内在原因，设计一种数字滤波对混频干扰进行抑制。从精密相位计算角度考虑，数字滤波必须能够完全滤除所述的混频干扰。数字滤波具体采用算术平均滤波算法，即对NT个连续离散值相加，然后取其算术平均值作为本次滤波值输出。为了方便计算数字滤波参数，在参考频率等于信号基波频率时，给出如下混频频率计算表1。

表1 混频频率计算表Tab.1 Calculation of mixing frequency

根据表1 计算结果，在数字滤波参数NT取值为2ω/3 频率单位周期序列长度时，可以对2ω/3和4ω/3 混频频率进行抑制。而NT取值为ω/2 频率单位周期序列长度时，可以对ω/2、3ω/2、1ω、2ω、3ω、4ω、5ω 混频频率进行抑制。因此，数字滤波由2 种参数的数字滤波所构成，为了提高对混频干扰的抑制性能，每种参数的数字滤波均由参数相同的三级数字滤波组成。

数字滤波为式(13):

式中:D1(n)为数字滤波1 输出序列;X(n)为数字滤波1 输入序列;NT1为滤波参数1、即连续离散值相加数量;D2(n)为数字滤波2 输出序列;NT2为滤波参数2、即连续离散值相加数量。2 级数字滤需要使用10.5 倍信号周期序列长度。

数字滤波用于实频序列，数字滤波输入序列为XRA(n)，数字滤波输出序列为XRAD(n);而用于虚频序列，数字滤波输入序列为XIA(n)，数字滤波输出序列为XIAD(n)。

实际根据参考频率计算滤波参数，在采样频率10 kHz，参考频率100 πrad/s，则NT1=300，NT2=400，计算得到数字滤波频域特性如图4所示。

图4 中，K(ω)为数字滤波频域幅频增益、单位dB，给出的频域特性对表1 给出的混频频率具有完全的抑制作用。

由于实际存在误差，包括参考频率误差，数字滤波参数误差。在参考频率变化时，数字滤波参数存在1 个采样间隔内的误差，通过改进算法，可将数字滤波参数误差控制在0.5 个采样间隔内。

图4 数字滤波频域特性示意图Fig.4 Schematic diagram of frequency characteristics of digital filter

在参考频率误差不大于±0.25%或者滤波参数误差在0.5 个采样间隔以内，图4所示的频域滤波特性对混频干扰仍然具有良好的抑制特性。其中在基波频率100 πrad/s，参考频率100.25 πrad/s，滤波参数NT1=299、NT2=399，得到结果如图5所示。垂直线为表1 给出的混频频率点，图5给出的最小抑制度为－210 dB，换算成衰减量约为3.16 ×10－11。

图5 混频频率抑制特性示意图Fig.5 Mixing frequency component suppression characteristic diagram

以相位调制A 序列为例，在混频干扰得到完全抑制前提下，实频和虚频混频数字滤波序列为式(14):

式中:XRAD(n)为实频序列的数字滤波序列;XIAD(n)为虚频序列的数字滤波序列;K(Ω)为数字滤波在频差Ω 的增益，单位无量纲;β(Ω)为数字滤波在频差Ω 的移相，rad。

3.3 积分计算

以相位调制A 序列为例，对实频和虚频数字滤波序列进行积分计算，为式(15):

式中:RA为实频积分值;IA为虚频积分值。

积分计算使用0.5 倍信号周期序列长度是足够的。

3.4 相位计算

以相位调制A 序列为例，相位调制A 序列的相位计算为

式(16)给出相位计算范围在±π/2，当信号初相位超出±0.5π 时，可能会造成计算错误。因此对信号初相位有一定的要求，一般要求信号初相位不超出±0.45π 范围。通过改进相位计算方法，可将相位计算范围扩大为±π，这样允许信号初相位有更大的范围，可为±0.9π。式(16)的改进相位算法为式(17)

式中:PhAC为改进相位算法的计算相位、范围±π，＆代表与逻辑。

省略相位调制B 序列的相位计算过程，相位调制B 序列的相位计算为

4 精密相位差和正弦频率计算

根据式(17)和式(18)，相位差计算为

式(19)中，相位差计算得到了调制相位的相位差。

根据式(19)，频率计算为

5 仿真实验

进行了电力系统50 Hz 工频正弦频率和调相相位差计算仿真实验，仿真实验条件为:实验信号基波频率变化范围在45～55 Hz，信号的采样频率10 kHz，信号的离散数据量化位数24 bit，频率初测单元相对误差±0.25%，信号序列长度(尽量)取12 倍信号周期序列长度。

具体实验信号由基波，1/2、1/3、2、3、4、5 次谐波成分等构成:

在基波频率50 Hz、参考频率50.125 Hz，得到的计算结果，如表2所示。

表2 新型正弦频率测量方法实验结果Tab.2 Experimental results of novel sinusoidal frequency measurement

在信号基波频率45～55 Hz 变化、参考频率相对误差0.25%时，得到基波频率计算相对误差随信号基波频率变化的实验结果，如图6所示。图中，|Ferr(f)|为(随信号基波频率变化)基波频率计算相对误差绝对值、单位无量纲;f 为信号基波频率，Hz。分析图6 实验结果，基波频率相对误差表现出明显的随机性，产生原因主要是离散数据量化背景噪声引起的，也表明数字滤波能够对混频干扰频率进行深度抑制，滤波后的残余幅值已低于背景噪声水平。

图6 基波频率相对误差实验结果1Fig.6 The fundamental frequency of relative error of experimental results A

另外，为了考查新型频率测量方法的抗噪声干扰特性，进行了白噪声加扰实验，所谓白噪声是指功频谱密度为常数的随机噪声，通常用信噪比衡量信号的质量，表述式为

式中:S:N 为功率信噪比，dB，Xs(n)为信号序列;Es为信号序列在序列长度N 的方差;Xn(n)为白噪声序列，En为白噪声序列在序列长度N的方差。

在信号基波频率50 Hz 和信号窗口时间0.25 s，参考频率为50.125 Hz，得到频率相对误差随信噪比变化的实验结果，如图7所示。

图7 基波频率相对误差实验结果2Fig.7 The fundamental frequency of relative error of experimental results B

图中:|Ferr(S:N)|为(随信号信噪比变化)基波频率计算相对误差绝对值，单位无量纲;S:N 为信号信噪比，dB。图7 给出了白噪声干扰环境下的频率计算相对误差分布图，其中在信噪比40 dB时可实现10－6量级准确度的频率测量。

6 物理实验

进行了电力系统50 Hz 工频频率测量的物理实验，这里指采集实际的高准确度信号发生器或实际电力系统的信号进行频率计算。具体物理实验条件为:实验频率测量系统的频率基准采用准确度±1 ×10－8量级的恒温晶振，采集设备的数据量化位数为24 bit，采集设备的采样频率为10 kHz。

物理实验结果表明，新型频率测量方法用于正弦频率测量具有较高的准确度，在45～55 Hz 频率范围内，在矩形窗口时间0.25 s 得到的正弦频率计算准确度优于±5.6 ×10－7，在窗口时间1.0 s 得到的正弦频率计算准确度优于±3.1 ×10－8，如图8所示。

图8 频率相对误差物理实验结果Fig.8 The relative error of the results of physical experiments frequency diagram

图中:Ferr(t)为(随过程时间变化)基波频率计算相对误差值，单位无量纲;t 为过程时间，s。另外，采集实际电力信号进行频率计算，同时与“零交法”频率测量进行对比，得到的结果，如图9所示。

图9 中F(t)为(随过程时间变化)基波频率计算值，Hz;t 为过程时间，s。图9所示，在20 s 时间内，信号频率呈缓慢变化趋势，采用新型正弦频率计算方法得到结果的波动幅度相对较小，而“零交法”频率测量结果的波动幅度相对较大，可见新型正弦频率计算方法相对“零交法”能够更真实的反映实际频率变化趋势。

图9 实际电力信号频率计算实验结果Fig.9 Calculating the results of experiments schematic diagram of the power frequency signal

7 结论

精密相位计算方法的本质是复数积分，文章换一个角度看问题，将复数积分看成是一种正交混频器，将频谱泄漏看成是一种混频干扰，指出了混频干扰是是造成正弦频率计算误差的主要内在原因。文章所设计的数字滤波，本质上是多种窗口特性的合成，对混频干扰具有深度的抑制作用，可实现较高准确度的低频正弦频率计算。通过数学计算、仿真试验和物理实验结果证明了所提出的新型频率测量方法的正确性和有效性。本文所提出的方法在电力科学研究、低频正弦频率的测量、低频率范围的精密测量仪器的研制等方面具有重要的用途和参考价值。

［1］肖遥，孟·让·柯洛德.电力系统频率测量误差成因分析［J］.电网技术，2002，26 (1):39－42.

［2］李军，王越超.一种基于幅值调制的新型电力系统正弦频率测量方法［J］.电工技术学报，2015，30(7):144－150.

［3］王兴国，黄少锋.基于复解析带通滤波器的固有频率自适应提取原理和方法［J］.电工技术学报，2009，24 (12):179－184.

［4］朱锋峰，任震，黄雯莹.基于小波变换修正幅值的电力系统频率偏移诊断方法［J］.电网技术，2004，28 (11):34－37.

［5］刘涤尘，夏利民，商志会.基于人工神经网络的电网频率测量方法［J］.电网技术，2000，24 (8):40－43.

［6］肖鲲，王莉娜，Shahzad M K.基于三线DFT 的航空电源频率实时检测算法［J］.电工技术学报，2012，27 (10):190－195.

［7］林达斌，江亚群，黄娟.基于一阶导数的电力系统频率测量新算法［J］.电测与仪表，2011，48(4):9－12.

［8］周卫平，吴正国，夏立.基波相位和频率的高精度检测及在有源电力滤波器中的应用［J］.中国电机工程学报，2004，24 (4):91－96.

［9］雷银照.我国供用电频率50Hz 的起源［J］.电工技术学报，2010，25 (3):20－26.

［10］王坚，陈明辉，文丽，等.评价电力系统频率质量的新指标［J］.华北电力大学学报，2010，37(3):7－13.

［11］星成武，马立新，孙大帅，等.基于频率跟踪的谐波分析［J］.华北电力大学学报，2010，37 (5):43－46.

［12］李军，万文军，张曦.一种基于阶跃响应的理想频率信号源及频域分析的研究［J］.动力工程学报，2012，32 (4):308－314.

［13］万文军，李军.一种基于LCR 发散振荡响应的控制系统频率特性辨识方法［J］.电机与控制学报，2014，18 (11):112－118.

［14］李军，万文军，刘志刚.一种阶跃函数在矩形时间窗口频域特性的分析方法［J］.电力自动化设备，2013，33 (11):111－116.

［15］李军，万文军，刘志刚，等.一种基于时域响应的控制系统频率特性分析方法［J］.中国电机工程学报，2012，32 (29):116－122.