CN104020352A - 一种适用于m类pmu单元的同步相量测量方法 - Google Patents

Abstract

本发明公开了一种适用于M类PMU单元的同步相量测量方法，该方法包括：将针对相量因子的低通数字滤波器与离散傅立叶算法DFT相结合，计算得到消除频谱泄漏后的初始量测相量；利用二阶泰勒级数对所述动态相量进行拟合，获得DFT平均化效应造成的量测误差，并基于该误差对所述初始相量进行补偿，获得精确的动态量测相量；在相量上传处设置一针对M类相量量测单元PMU的低阶数字滤波器，将所述精确的量测相量作为输入，获得最终量测相量。通过采用本发明公开的方法，在输入为静态信号及动态信号时，均可准确且快速地进行相量量测，其相量量测精度可以满足IEEE C37.118.1的要求，并普遍可高出一个数量级。

Description

一种适用于M类PMU单元的同步相量测量方法

技术领域

本发明涉及同步相量测量技术领域，尤其涉及一种适用于M类PMU单元的同步相量测量方法。

背景技术

相量测量单元(Phasor Measurement Unit，PMU)的应用对电力系统的量测技术带来了革命性的变革。除了可提供同步相量，其高精度及高速的上传频率的优点，亦成为其在动态安全监控中广泛地作为相量数据源的原因。IEEE C37.118.1规定：M类PMU主要用于电力系统信号的同步量测及分析等应用，它不需要具有快速的响应速度，但要求精确相量量测的频带较宽，且必须具有消除带外频率分量的功能。因此，M类PMU能否在较宽的频带上获得精确的量测相量，且不受带外信号的影响，是基于M类PMU的应用能否对电力系统信号进行正确分析的关键。

离散傅立叶算法(Discrete Fourier transform，DFT)因为其可将额定频率分量从含有谐波分量的波形中提取出来以及其计算简单的特性，而在PMU中广泛地得到应用。但是，DFT是基于相量的各个参数在计算时间窗内不变，且频率为额定频率的假设基础上的。在电力系统动态过程中，这一假设是不成立的，相量的各个参数将随时间变化而变化。因此，由于DFT本身的缺陷将在相量量测中导致两个问题。第一，频率发生偏移以及动态信号输入都会造成频谱泄漏；第二，静态相量模型将导致DFT的平均化效应，从而在动态条件下产生较大误差。

2011年，发布了PMU标准IEEE C37.118.1。相比较于旧版PMU标准IEEE C37.118-2005，新的标准在量测测试范围上有了较大地升级。例如，增加了动态条件下的量测精度的规定，以及频率及频率变化率的量测精度要求。

随着IEEE以及中国PMU标准的发布及逐步完善，PMU动态条件下的量测精度受到越来越多的研究机构地重视，很多新的技术在算法中得到了应用。然而，目前并没有一个算法是基于动态相量模型，无法解决DFT因为动态输入信号而造成的频谱泄漏和平均化效应的问题。

发明内容

本发明的目的是提供一种适用于M类PMU单元的同步相量测量方法，该方法在输入为静态信号及动态信号时，均可准确且快速地进行相量量测，其相量量测精度可以满足IEEE C37.118.1的要求，并普遍可高出一个数量级。

本发明的目的是通过以下技术方案实现的：

一种适用于M类PMU单元的同步相量测量方法，该方法包括：

将针对相量因子的低通数字滤波器与离散傅立叶算法DFT相结合，计算得到消除频谱泄漏后的初始量测相量；

利用二阶泰勒级数对所述动态相量进行拟合，获得DFT平均化效应造成的量测误差，并基于该误差对所述初始相量进行补偿，获得精确的动态量测相量；

在相量上传处设置一针对M类相量量测单元PMU的低阶数字滤波器，将所述精确的量测相量作为输入，获得最终量测相量。

由上述本发明提供的技术方案可以看出，该测量方法无论是输入静态信号还是动态信号，都可以准确快速地进行相量量测，其相量量测精度不仅满足标准要求，并普遍高出标准一个数量级。该方法不仅可在频率偏移、电力系统振荡、失步以及故障时准确的量测相量，对整次谐波、非整次谐波及带外信号具有免疫功能。在该方法中采用了大量的查表技术，且并不需要较高的采样率，所以对现代的硬件来说计算负担并不重，可在实际装置中实现。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。

图1为本发明实施例一提供的一种适用于M类PMU单元的同步相量测量方法的流程图；

图2为本发明实施例一提供的幅值量测误差与频率的关系示意图；

图3为本发明实施例一提供的仿真频率扫面测试中不同基波频率下的TVE最大值示意图；

图4为本发明实施例一提供的仿真频率扫面测试中不同基波频率下的频率及频率变化率误差最大值示意图；

图5为本发明实施例一提供的仿真带外测试中不同调制频率下最大TVE示意图；

图6为本发明实施例一提供的仿真带外测试中不同调制频率下最大频率误差的示意图；

图7为本发明实施例一提供的仿真调制测试中不同调制频率下最大TVE示意图；

图8为本发明实施例一提供的仿真调制测试中不同调制频率下最大频率误差的示意图；

图9为本发明实施例一提供的仿真调制测试中不同调制频率下最大频率变化率误差的示意图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。

实施例一

图1为本发明实施例一提供的一种适用于M类PMU单元的同步相量测量方法的流程图。如图1所示，该方法主要包括如下步骤：

步骤11、将针对相量因子的低通数字滤波器与离散傅立叶算法DFT相结合，计算得到消除频谱泄漏后的初始量测相量。

本发明实施例中，所述针对相量因子的低通数字滤波器基于二阶泰勒级数的动态相量数学模型而获得；所述动态相量参数包括幅值、相角、频率以及频率变化率；上述参数在计算时间窗内都随时间变化而变化。

示例性的，首先将动态基频信号表示为：

其中，x(t)是信号的采样值，X_m(t)是相量幅值，f(t)是信号频率，是相量初相角，f₀是额定频率，Δf(t)是频率偏移；此时，信号的幅值与频率都是时间的函数。

将式(1)重新写为：

其中Re为实部。通常情况下认为参考坐标系以2πf₀的速度同步旋转，所以可忽略因此，式(1)可表示为这便是动态相量的数学模型公式：

频率变化率ROCOF为：

ROCOF＝d/dt[f(t)]＝Δf'(t)    (4)

其中，Δf'(t)为Δf(t)的导数。

因此，动态相量的所有参数都是时间的函数，为了在计算时间窗内逼近动态输入信号，利用二阶泰勒级数来模拟相量参数的非线性变化波形，如式(5-8)所示。

X_m(t)＝m₂t²+m₁t+m₀    (5)

f(t)＝p₂t²+p₁t+p₀    (7)

ROCOF(t)＝q₂t²+q₁t+q₀    (8)

其中，m₂＝d²X_m(t)/dt²|_t＝0，m₁＝dX_m(t)/dt|_t＝0，m₀＝X_m(0)，p₂＝d²f(t)/dt²|_t＝0，p₁＝df(t)/dt|_t＝0，p₀＝f(0)，q₂＝d²ROCOF(t)/dt²|_t＝0，q₁＝d ROCOF(t)/dt|_t＝0，q₀＝ROCOF(0)。

然后，基于上述动态相量的数学模型，分析了离散傅立叶算法(DFT)在动态相量输入条件下所产生的频谱泄露现象及对相量量测的影响，提出含有低通数字滤波器的DFT算法，以消除频谱泄露现象，并获得初始量测相量示例性的，上述与所述离散傅里叶算法DFT相结合的低通数字滤波器的截止频率可以为2.5Hz，时间窗长可以为动态相量输入信号的2个额定周期，初始相量计算频率可以为400Hz。

示例性的，现有技术中传统DFT算法公式为：

${\stackrel{\·}{X}}^{\′}=\frac{\sqrt{2}}{N}\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}x\left(k\right)e^{-\mathrm{jk}2\mathrm{\π}{\mathrm{\ω}}_0\mathrm{\Δt}}---\left(9\right)$

其中，为量测相量，N为采样个数，x(k)为第k个采样值，Δt为采样间隔。

可以看出，DFT的计算过程为将一个时间窗内的采样值与正交系数相乘所得到的乘积相加并取平均得到相量。现在将式(9)中采样值与正交系数的乘积用数列γ表示：

$\mathrm{\γ}=\left\{{\mathrm{\γ}}_k\right|{\mathrm{\γ}}_k=x\left(k\right)e^{\mathrm{jk}2\mathrm{\π}{f}_0\mathrm{\Δt}}\},k=\mathrm{0,1},...,N-1---\left(10\right)$

分别将γ的实部与虚部提取出来，得到：

C_x＝Re(γ)＝XC    (11)

S_x＝Im(γ)＝XS    (12)

其中，X＝{diag[x(k)；n＝0,1,…,N-1]}，C^T＝[cos[2πf₀(0)Δt],cos[2πf₀(1)Δt],…,cos[2πf₀(N-1)Δt]]，S^T＝[sin[2πf₀(0)Δt],sin[2πf₀(1)Δt],…,sin[2πf₀(N-1)Δt]]。

C_x和S_x即为相量因子，其特性将直接关系到DFT的量测特性。将式(9)用相量因子可表示为：

${\stackrel{\·}{X}}^{\′}=\frac{\sqrt{2}}{N}\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{C}_x\left(k\right)-j{S}_x\left(k\right)---\left(13\right)$

若输入信号为静态的，且频率为额定值：

则C_x和S_x为：

$\begin{array}{c}{C}_x\left(k\right)=x\left(k\right)\cos \left(2\mathrm{\π}{f}_0\mathrm{k\Δt}\right)\\ =\frac{\sqrt{2}}{2}{X}_m\cos (4\mathrm{\πk}{f}_0\mathrm{\Δt}+{\mathrm{\φ}}_0)+\frac{\sqrt{2}}{2}{X}_m\cos \left({\mathrm{\φ}}_0\right)\end{array}---\left(15\right)$

$\begin{array}{c}{S}_x\left(k\right)=x\left(k\right)\sin \left(2{\mathrm{\π}}_0\mathrm{k\Δt}\right)\\ =\frac{\sqrt{2}}{2}{X}_m\sin (4\mathrm{k\π}{f}_0\mathrm{\Δt}+{\mathrm{\φ}}_0)-\frac{\sqrt{2}}{2}{X}_m\sin \left({\mathrm{\φ}}_0\right)\end{array}---\left(16\right)$

如式(15)和(16)所示，当k增加时，相量因子是由一个2f₀分量和一个直流分量组成的，对一个周期内的相量因子取平均，可把2f₀分量消除，从而获得真实相量，如式(17)所示，这就是整周期DFT算法。

上述DFT算法因为其可消除谐波并提取基波的优点而得到广泛地应用，但当频率偏移额定值时会发生频谱泄露，尤其当频率随时间变化而变化时，消除此频谱泄露较为困难，当输入动态信号时，亦会发生频谱泄露。

示例性的，将式(1)作为输入信号，则C_x和S_x变为：

$\begin{array}{c}{C}_x\left(k\right)=x\left(k\right)\cos \left({\mathrm{\ω}}_0\mathrm{k\Δt}\right)\\ =\frac{\sqrt{2}}{2}{X}_m\left(\mathrm{k\Δt}\right)\cos [2\mathrm{\π}{\∫}_0^{\mathrm{k\Δt}}f\left(t\right)\mathrm{dt}+2\mathrm{\πk}{f}_0\mathrm{\Δt}+{\mathrm{\φ}}_0]+\frac{\sqrt{2}}{2}\\ =C_1+C_2\end{array}{X}_m\left(\mathrm{k\Δt}\right)\cos (2\mathrm{\π}{\∫}_0^{\mathrm{k\Δt}}\mathrm{\Δf}\left(t\right)\mathrm{dt}+{\mathrm{\φ}}_0)---\left(18\right)$

$\begin{array}{c}{S}_x\left(k\right)=x\left(k\right)\sin \left({\mathrm{\ω}}_0\mathrm{k\Δt}\right)\\ =\frac{\sqrt{2}}{2}{X}_m\left(\mathrm{k\Δt}\right)\sin [2\mathrm{\π}{\∫}_0^{\mathrm{k\Δt}}f\left(t\right)\mathrm{dt}+2\mathrm{\πk}{f}_0\mathrm{\Δt}+{\mathrm{\φ}}_0]-\frac{\sqrt{2}}{2}\\ =S_1+S_2\end{array}{X}_m\left(\mathrm{k\Δt}\right)\sin (2\mathrm{\π}{\∫}_0^{\mathrm{k\Δt}}\mathrm{\Δf}\left(t\right)\mathrm{dt}+{\mathrm{\φ}}_0)---\left(19\right)$

相量真实值应如式(3)所示，由式(18)、(19)中的C₂和S₂组成。但是由于C₁和S₁的频率并不是f₀的整数倍，所以不能完全消除，因此便造成了频谱泄露。

如上述式(18)和(19)所示，导致频谱泄漏的频率分量的频率都在2f₀附近，而组成真实相量的频率分量的频率较低，IEEE C37.118.1对P类PMU的要求中，Δf和f₁都在在2Hz以下。所以，可利用低通数字滤波器在滤除相量因子中的高频分量的同时保留低频分量，因此本发明实施例提出带有低通数字滤波器的改进DFT算法，以消除频谱泄漏，从而得到初始量测相量其公式如式(20)所示：

$\stackrel{\widetilde{\·}}{X}=\frac{\sqrt{2}}{S}\underset{k=-\frac{N-1}{2}}{\overset{\frac{N-1}{2}}{\mathrm{\Σ}}}x(k+\frac{N-1}{2})h(k+\frac{N-1}{2})e^{-\mathrm{jk}{\mathrm{\ω}}_0\mathrm{\Δt}}---\left(20\right)$

其中，为初始量测相量，为低通滤波器系数， 为第个采样值。

值得注意的是，由于将时标打在了时间窗中部，所以k的初始值是–(N-1)/2。根据初始量测相量可进一步得到初始量测相量幅值与相角

进一步的，可以将式(24)重写为相量因子的形式：

可以看出，虽然频谱泄漏已被消除，但是量测相量的幅值取决X_m(t)，Δf以及滤波器系数。只有当X_m(t)为常数Δf(t)＝0时，才可获得准确的相量否则，幅值量测将会出现误差；但是，一旦数字滤波器确定后，幅值误差只与Δf有关。因此，将幅值误差看作Δf的二次函数进行拟合，可以得到两者的函数关系，如图2所示为幅值量测误差与频率的关系示意图 $\mathrm{Magnitude\; Error}=\frac{\left|\right|\stackrel{\·}{X}|-|\stackrel{\widetilde{\·}}{X}\left|\right|}{\left|\stackrel{\·}{X}\right|}\×100\%.$

步骤12、利用二阶泰勒级数对所述动态相量进行拟合，获得DFT平均化效应造成的量测误差，并基于该误差对所述初始相量进行补偿，获得精确的动态量测相量。

由于通过上述步骤11中加装数字滤波器的DFT，可以得到消除频谱泄漏后的初始相量，但是由于DFT仍是基于静态相量模型，所以当输入动态信号时，初始量测相量仍存在较大误差，无法满足标准要求。

示例性的，当电力系统发生系统振荡时，相量的幅值、相角、频率以及频率变化率都将发生非线性变化，下面以计算精确的动态量测相量幅值为例进行说明。

首先，将式(5)中的X_m(t)表达式带入式(21)中，此处为了简化推导过程，假设Δf(t)＝0，则化简之后的表达式为：

其中， $K_2=[{(-\frac{N-1}{2})}^2,{(-\frac{N-1}{2}+1)}^2,...,0,...,{\left(\frac{N-1}{2}\right)}^2],$ E＝[1,1…1]， $K_1=[-\frac{N-1}{2},-\frac{N-1}{2}+1,...,0,...,\frac{N-1}{2}],$ H＝[h(0),h(1),…,h(N-1)]^T；

由于数字滤波器的系数为对称结构，可以将H改写为： $H={[h\left(0\right),h\left(1\right),..,h\left(\frac{N-1}{2}\right),...,h\left(1\right),h\left(0\right)]}^T,$ 则有：

$\mathrm{EH}/s=\frac{h\left(0\right)+h\left(1\right)+...+h(N-1)}{h\left(0\right)+h\left(1\right)+...+h(N-1)}=1---\left(23\right)$

$\mathrm{EH}/s=\frac{h\left(0\right)+h\left(1\right)+...+h(N-1)}{h\left(0\right)+h\left(1\right)+...+h(N-1)}=1---\left(24\right)$

$K_{2}H/s=\frac{h\left(0\right){(-\frac{N-1}{2})}^2+h\left(1\right){(-\frac{N-1}{2}+1)}^2+...+h(N+1){\left(\frac{N-1}{2}\right)}^2}{h\left(0\right)+h\left(1\right)+...+h(N-1)}={\tilde{K}}_2---\left(25\right)$

其中，为K₂的加权平均值；将式(23)-(25)带入式(22)中，得到：

由于时标打在时间窗中部，所以相量的真实值为则幅值量测误差为：

$e_m=\left|\stackrel{\widetilde{\·}}{X}\right|-\left|\stackrel{\·}{X}\right|=m_2{\mathrm{\Δt}}^2{\tilde{K}}_2=\mathrm{\α}{m}_2---\left(27\right)$

其中，对于一个在装置中实施的算法，其采样间隔与滤波器系数都是固定的，因此，Δt与为固定值。由此可知，相量的幅值量测误差与m₂存在线性关系，且线性关系的斜率与采样间隔和数字滤波器系数有关。

为了补偿动态量测误差(幅值误差)，需要计算m₂，设相量M的公式为：

$\left[\begin{array}{c}{\tilde{X}}_m^0\\ {\tilde{X}}_m^1\\ {\tilde{X}}_m^2\\ .\\ .\\ .\\ {\tilde{X}}_m^{M-1}\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 1& {\mathrm{\Δt}}_c& {{\mathrm{\Δt}}_c}^2\\ 1& 2{\mathrm{\Δt}}_c& 2^2{{\mathrm{\Δt}}_c}^2\\ .& .& .\\ .& .& .\\ .& .& .\\ 1& (M-1){\mathrm{\Δt}}_c& {(M-1)}^1{{\mathrm{\Δt}}_c}^2\end{array}\right]\left[\begin{array}{c}{m}_0\\ m_1\\ m_2\end{array}\right]---\left(28\right)$

其中，为量测相量幅值，r＝0,1,2,...,M-1；Δt_c为计算间隔；上市用矩阵形式可表示为：

${\tilde{X}}_m=T_c{C}_m---\left(29\right)$

其中，T_c为系数矩阵；位置矩阵C_m可通过最小而乘法计算，表示为：

$C_m={\left[{T_c}^T{T}_c\right]}^{-1}{T_c}^T{\tilde{X}}_m---\left(30\right)$

同样，由于时标在时间窗中部，通过式(31)曲线拟合可得到

${\tilde{X}}_m^{\frac{M-1}{2}}={\left(\frac{M-1}{2}{\mathrm{\Δt}}_c\right)}^2{m}_2+\frac{M-1}{2}{\mathrm{\Δt}}_c{m}_1+m_0---\left(31\right)$

最后，根据图2对进行频率偏移静态补偿，得到再进一步根据式(27)，利用m₂对进行精度校准：

${\tilde{X}}_{\mathrm{mf}}={\tilde{X}}_{m1}-{\mathrm{\αm}}_2---\left(32\right)$

其中，补偿系数α可由仿真实验来确定。

另一方面，n₁、n₂也可以根据相角通过类似的方法计算得到，之后初始的频率和频率变化率可由式(33)、(34)计算得到：

$\tilde{f}=n_2\frac{M-1}{2}{\mathrm{\Δt}}_c/\mathrm{\π}+n_1/2\mathrm{\π}---\left(33\right)$

$\mathrm{ROCO}\tilde{F}=n_2/\mathrm{\π}---\left(32\right)$

对于随时间变化的f(t)、以及ROCOF(t)，它们的量测误差也与各自的泰勒表达式的二次系数呈线性关系，这通过仿真实验就可以得到验证。同样，可采用上述类似的方法进行计算，故不再赘述。

步骤13、在相量上传处设置一针对M类相量量测单元PMU的低阶数字滤波器，将所述精确的量测相量作为输入，获得最终量测相量。

根据IEEE C37.118.1的规定，M类PMU应能消除带外信号，其带外信号的频带由PMU上传频率(F_s)的奈奎斯特频率决定。

由式(18)和(19)可看出，到f₁大于F_s的奈奎斯特频率时，可通过将低通滤波器的截止频率设置到f₁以下的方法来实现。为了满足IEEE C37.118.1中带外测试的规定，低通滤波器需要有足够长时间窗。由于本发明实施例中的方法用二阶泰勒展开式来模拟一个时间窗内的动态信号，因此需要准确量测的调制信号的调制周期需至少是时间窗的两倍。IEEE C37.118.1中调制测试中的最大调制频率为5Hz，也就是说滤波器的时间窗最长为100ms。但是，对于上传频率较低的PMU来说，100ms的低通数字滤波器是不够的。

因此，对于M类PMU本发明在动态相量校准后使用了第二个低通数字滤波器。其时间窗长取决于是否能够满足带外测试，但同时也受带外测试中动态响应时间的限制。经该数字滤波器后得到最终量测相量，N’为第二个低通数字滤波器的时间窗长。

为了进一步说明本发明，下面再以具体的实例对上述测量方法进行仿真测试，具体来说：

IEEE C37.118.1中规定了静动态测试以及相量量测最大误差，以完整全面地模拟电力系统中静动态过程。本发明提出的算法根据标准进行了包括频率偏移、谐波影响、振荡、失步及故障条件下的仿真。算法的仿真误差与标准规定的误差进行了对比，结果表明发明的算法在静动态条件下都可很好的满足标准要求。仿真结果所用算法的额定频率为50Hz，上传频率为50Hz，采样频率为1200Hz，N＝48，M＝7，N’＝80。算法中如式(30)C_m等参数可提前计算，采用相量估计过程中采用查表方法。整个相量、频率及频率变化率的计算过程中共有465次乘法，403次加法。这对现代的处理器来说计算负担并不重。

1、频率扫描测试

不同电力系统的运行模式将导致电力系统信号的实际频率偏移其额定值。此外，故障亦会导致较大的频率偏移。此处列出了本发明实施例的方法在不同程度的频率偏移条件下的量测精度。输入信号的频率以1Hz的步长从45Hz变为55Hz，且在每个频率点保持稳态。IEEE C37.118.1中的频率测试中规定最大的综合矢量误差(Total vector error,TVE)为1％，最大频率误差(Frequency error,FE)为0.005Hz，最大频率变化率误差(ROCOF error,RFE)为0.01Hz/s。相应的本发明方法的测试结果如图3和图4所示，可以看出本发明方法的量测精度远高于标准要求。

2、谐波影响测试

如表1所示，本次测试采用了在49.5Hz、50Hz及50.5Hz下的10％的二次、三次、八次及十三次谐波的电力系统信号来验证发明算法的消除谐波的能力，标准要求的误差要求(std.)也在表1中列出。

表1  谐波影响下的相量量测误差统计

表1列出了算法在包含不同谐波次数的信号下的量测误差。可以看出，偶数次谐波对算法的影响比奇数次谐波要大。但是，所有的误差都远在标准误差极限以下。所以发明的方法对谐波影响有较好的免疫能力。

3、带外信号测试

当信号的变化频率大于Fs的奈奎斯特频率时，会发生频率混叠。带外信号测试的目的是验证PMU单元的带外信号抗混叠滤波器的滤波效果。测试的方法是在基波信号的基础上叠加一个10％的带外正弦信号。由于F_s＝50Hz，所以叠加的信号频率分别为10Hz、20Hz、25Hz、75Hz、90Hz和100Hz。IEEE标准中规定带外测试中TVE和频率的最大误差1.3％和0.01Hz。算法的测试结果如图5和图6所示，可以看出越靠近带通边缘的信号对量测精度影响越大。但是测试结果表明发明的方法可以有效地抑制带外信号对量测精度的影响。

4、调制信号测试

调制信号测试用来模拟电力系统发生振荡时信号的幅值和相角的波动波形。通常来说，电网节点的正序电压信号的幅值与相角同时发生振荡，且振荡角度相差180°。在测试过程中，信号的幅值与相角以正弦波形发生变化。一般误差在振荡的波峰或波谷等非线性程度最大的地方最大。

测试中幅值调制深度为10％，相角调制深度为0.1rad，调制频率从0.1Hz到5Hz。最大TVE规定为3％，最大频率误差为0.3Hz，最大频率变化率误差为30Hz/s。本发明的方法的测试结果如图7、图8和图9所示。可以看出，量测误差随着调制频率的增大而增大。这是由于随着调制频率越大，信号在一个时间窗内的变化也越来越快。不过相比于误差要求，可证明发明的算法具有较为精确的动态信号跟踪能力。

5、频率斜坡测试

频率斜坡测试用来模拟电力系统失步过程，不同于频率扫描测试，其基频从45Hz以1Hz/s连续变化至55Hz。通过测试可知，发明的算法可在持续变化的频率下准确的量测相量、频率、频率变化率，如表2所示。

表2  频率斜坡测试误差统计

6、阶跃测试

当电力系统中发生故障或者换线操作时，电压与电流波形的幅值和相角会发生突变。对M类PMU来说，其动态相应速度要求没有P类PMU高。但是，适当的动态跟踪速度对其在电力系统动态安全监控也至关重要。

该测试中，输入信号会分别发生10％的幅值阶跃和10°。发明的算法的相应时间如表III所示。可以看出，算法可满足相应时间的标准要求。

表3  阶跃测试响应时间

本发明实施例提供的测量方法无论是输入静态信号还是动态信号，都可以准确快速地进行相量量测，其相量量测精度不仅满足标准要求，并普遍高出标准一个数量级。该方法不仅可在频率偏移、电力系统振荡、失步以及故障时准确的量测相量，对整次谐波、非整次谐波及带外信号具有免疫功能。在该方法中采用了大量的查表技术，且并不需要较高的采样率，所以对现代的硬件来说计算负担并不重，可在实际装置中实现。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (5)

1.一种适用于M类PMU单元的同步相量测量方法，其特征在于，该方法包括：

将针对相量因子的低通数字滤波器与离散傅立叶算法DFT相结合，计算得到消除频谱泄漏后的初始量测相量；

利用二阶泰勒级数对所述动态相量进行拟合，获得DFT平均化效应造成的量测误差，并基于该误差对所述初始相量进行补偿，获得精确的动态量测相量；

在相量上传处设置一针对M类相量量测单元PMU的低阶数字滤波器，将所述精确的量测相量作为输入，获得最终量测相量。

2.根据权利要求1所述的方法，其特征在于，所述针对相量因子的低通数字滤波器基于二阶泰勒级数的动态相量数学模型而获得；

所述动态相量参数包括幅值、相角、频率以及频率变化率。

3.根据权利要求2所述的方法，其特征在于，所述动态相量的数学模型公式为：

频率变化率ROCOF为：ROCOF＝d/dt[f(t)]＝Δf'(t)；

其中，x(t)是信号的采样值，X_m(t)是相量幅值，f(t)是信号频率，是相量初相角，Δf(t)是频率偏移，Δf’(t)为Δf(t)的导数；

利用二阶泰勒级数来模拟相量参数的非线性变化波形，表示为：

X_m(t)＝m₂t²+m₁t+m₀；

f(t)＝p₂t²+p₁t+p₀；

ROCOF(t)＝q₂t²+q₁t+q₀；

其中，m₂＝d²X_m(t)/dt²|_t＝0，m₁＝dX_m(t)/dt|_t＝0，m₀＝X_m(0)， p₂＝d²f(t)/dt²|_t＝0，p₁＝df(t)/dt|_t＝0，p₀＝f(0)，q₂＝d²ROCOF(t)/dt²|_t＝0，q₁＝d ROCOF(t)/dt|_t＝0，q₀＝ROCOF(0)。

4.根据权利要求1或2或3所述的方法，其特征在于，所述计算得到消除频谱泄漏后的初始量测相量的公式为：

$\stackrel{\widetilde{\·}}{X}=\frac{\sqrt{2}}{S}\underset{k=-\frac{N-1}{2}}{\overset{\frac{N-1}{2}}{\mathrm{\Σ}}}x(k+\frac{N-1}{2})h(k+\frac{N-1}{2})e^{-\mathrm{jk}{\mathrm{\ω}}_0\mathrm{\Δt}};$

其中，为初始量测相量，为低通滤波器系数，

5.根据权利要求4所述的方法，其特征在于，所述获得精确的动态量测相量包括：获得精确的动态量测相量幅值、相角、频率以及频率变化率；

其中，获得精确的动态量测相量幅值的步骤包括：

将初始量测相量的计算公式改写为相量因子的形式：

将相量幅值X_m(t)＝m₂t²+m₁t+m₀带入上式，且设Δf(t)＝0，经过化简后得到式子：

其中， $K_2=[{(-\frac{N-1}{2})}^2,{(-\frac{N-1}{2}+1)}^2,...,0,...,{\left(\frac{N-1}{2}\right)}^2],$ E＝[1,1…1]， $K_1=[-\frac{N-1}{2},-\frac{N-1}{2}+1,...,0,...,\frac{N-1}{2}],$ H＝[h(0),h(1),…,h(N-1)]^T；

数字滤波器的系数为对称结构，将H改写为： $H={[h\left(0\right),h\left(1\right),..,h\left(\frac{N-1}{2}\right),...,h\left(1\right),h\left(0\right)]}^T,$ 则：

$\mathrm{EH}/s=\frac{h\left(0\right)+h\left(1\right)+...+h(N-1)}{h\left(0\right)+h\left(1\right)+...+h(N-1)}=1;$

$\mathrm{EH}/s=\frac{h\left(0\right)+h\left(1\right)+...+h(N-1)}{h\left(0\right)+h\left(1\right)+...+h(N-1)}=1;$

$K_{2}H/s=\frac{h\left(0\right){(-\frac{N-1}{2})}^2+h\left(1\right){(-\frac{N-1}{2}+1)}^2+...+h(N+1){\left(\frac{N-1}{2}\right)}^2}{h\left(0\right)+h\left(1\right)+...+h(N-1)}={\tilde{K}}_2;$

其中，为K₂的加权平均值；

将上述三个表达式带入化简后的表达式，得到：

则幅值量测误差为： $e_m=\left|\stackrel{\widetilde{\·}}{X}\right|-\left|\stackrel{\·}{X}\right|=m_2{\mathrm{\Δt}}^2{\tilde{K}}_2=\mathrm{\α}{m}_2;$

其中，为相量真实值： $\mathrm{\α}={\mathrm{\Δt}}^2{\tilde{K}}_2,$ Δt与为固定值；

m₂则根据相量M来计算，相量M表示为：

$\left[\begin{array}{c}{\tilde{X}}_m^0\\ {\tilde{X}}_m^1\\ {\tilde{X}}_m^2\\ .\\ .\\ .\\ {\tilde{X}}_m^{M-1}\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 1& {\mathrm{\Δt}}_c& {{\mathrm{\Δt}}_c}^2\\ 1& 2{\mathrm{\Δt}}_c& 2^2{{\mathrm{\Δt}}_c}^2\\ .& .& .\\ .& .& .\\ .& .& .\\ 1& (M-1){\mathrm{\Δt}}_c& {(M-1)}^1{{\mathrm{\Δt}}_c}^2\end{array}\right]\left[\begin{array}{c}{m}_0\\ m_1\\ m_2\end{array}\right];$

其中，为量测相量幅值，r＝0,1,2,...,M-1；Δt_c为计算间隔；

将上式转换为矩阵形式，则有：其中，T_c为系数矩阵；矩阵 $C_m={\left[{T_c}^T{T}_c\right]}^{-1}{T_c}^T{\tilde{X}}_m;$

再利用曲线拟合计算其公式为：

${\tilde{X}}_m^{\frac{M-1}{2}}={\left(\frac{M-1}{2}{\mathrm{\Δt}}_c\right)}^2{m}_2+\frac{M-1}{2}{\mathrm{\Δt}}_c{m}_1+m_0;$

对所述进行静态补偿，获得再进行精度校准，获得精确的动态量测相量幅值

${\tilde{X}}_{\mathrm{mf}}={\tilde{X}}_{m1}-{\mathrm{\αm}}_2.$