CN105759118A - 一种同步相量的量测方法 - Google Patents

Abstract

本发明提供一种同步相量的量测方法，所述量测方法利用同步相量动态信号模型，用Legendre多项式对信号进行拟合；计算拟合系数并获得瞬时同步相量。本发明提供的同步相量量测方法降低了矩阵条件数，可准确且快速地进行相量量测，提高相量量测精度。

Description

一种同步相量的量测方法

技术领域

本发明涉及同步相量测量技术领域，具体讲涉及一种同步相量的量测方法。

背景技术

相量测量单元(PhasorMeasurementUnit，PMU)除提供同步相量外，其高精度及高速的上传频率使其作为相量数据源在动态安全监控中广泛应用。相量测量单元的应用对电力系统的量测技术带来了革命性的变革。

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

2014年，IEEE发布了PMU国际标准IEEEC37.118.1a，相比较于旧版PMU标准IEEEC37.118-2005，新的标准修订了PMU在动态条件下的量测精度要求。国内也将于2015年颁布国家电网PMU标准《电力系统实时动态监测系统技术规范》，也全面规定了在电力系统静动态条件下PMU的量测精度要求。

随着IEEE以及中国PMU标准的发布及逐步完善，PMU动态条件下的量测精度受到越来越多的研究机构的重视，很多新的技术在算法中得到了应用。然而，大多数商用PMU仍采用的是DFT算法，而DFT算法基于静态相量模型，无法解决其在动态条件下的频谱泄漏和平均化效应的问题，动态量测精度无法满足新标准的要求。

为满足现有技术的发展需要，克服现有技术的缺陷，本发明提供一种高精度的同步相量的量测方法。

发明内容

本发明的目的是提供一种高精度的同步相量的量测方法，该量测方法中同步相量是一个动态数学模型，利用该动态模型，采用Legendre多项式对动态信号进行拟合，计算得到拟合系数之后，进一步获得精确的同步相量。Legendre多项式拟合方法属于线性最小二乘拟合算法，与传统Taylor-Fourie方法相比，降低了矩阵条件数，使得整个量测方法效率提高。该量测方法在输入为静态信号及动态信号时，均可准确且快速地进行相量量测，其相量量测精度可以满足国家电网企标《电力系统实时动态监测系统技术规范》。

一种同步相量的量测方法，其改进之处在于，所述量测方法包括以下步骤：

I、拟合动态信号；

II、计算拟合系数；

Ⅲ、确定瞬时同步相量。

进一步的，所述步骤I中，所述动态信号模型x(t)如下式所示：

其中，相量初相角；f₀：额定频率；Δf：频率偏移量；X_m(t)：相量幅值的时间函数；f(t)：信号频率的时间函数；

对式(1)所示的动态信号模型的拟合包括以下步骤：

(i)将式(1)变形为下式(2)：

其中，

(ii)式(2)所述的动态信号模型的展开式如下式所示：

(iii)用Legendre多项式展开式(3)中的权重函数和令则q(t)的k阶Legendre多项式展开式如下：

q(τ)＝q₀L₀(τ)+q₁L₁(τ)+…+q_kL_k(τ)+o(τ^k)(4)

其中，L_k(τ)：多项式的第k阶项；q_k：多项式第k阶项系数；o(τ^k)：多项式展开式的皮亚诺型余项；k为多项式展开的阶数，k＝0，1，2…n；

令则r(t)的k阶Legendre多项式展开式如下：

r(τ)＝r₀L₀(τ)+r₁L₁(τ)+…+r_kL_k(τ)+o(τ^k)(5)

其中，L_k(τ)：多项式的第k阶项；r_k：多项式展开式中第k阶项系数；o(τ^k)：多项式展开式的皮亚诺型余项；k为多项式展开的阶数，k＝0，1，2…n。

进一步的，时间坐标t∈[a,b]数据窗内的数据点的线性坐标变换如下式所示：

τ＝(2t-a-b)/(b-a)；(6)

将步骤(iii)中的式(4)和式(5)所示的展开式代入步骤(ii)的式(3)中，得如下式(7)所示的τ时间坐标下成立的动态信号模型x(τ)：

x(τ)＝[q₀L₀(τ)+q₁L₁(τ)+…+q_kL_k(τ)]cos(2πf_lτ)+[r₀L₀(τ)+r₁L₁(τ)+…+

r_kL_k(τ)]sin(2πf_lτ)＝

q₀L₀(τ)cos(2πf_lτ)+q₁L₁(τ)cos(2πf_lτ)+…+q_kL_k(τ)cos(2πf_lτ)+

r₀L₀(τ)sin(2πf_lτ)+r₁L₁(τ)sin(2πf_lτ)+…+r_kL_k(τ)sin(2πf_lτ)(7)

其中，f_l：t时间坐标下的额定频率f₀映射到τ时间坐标下的对应频率。

进一步的，步骤II中所述拟合系数的计算包括：步骤(I)、将采样点的时间坐标和采样值代入式(7)，得式(8)所示的τ时间坐标下的动态信号模型x(τ_i)的表示式：

x(τ_i)＝q₀L₀(τ_i)cos(2πf_lτ_i)+q₁L₁(τ_i)cos(2πf_lτ_i)+…+q_kL_k(τ_i)cos(2πf_lτ_i)+

r₀L₀(τ_i)sin(2πf_lτ_i)+r₁L₁(τ_i)sin(2πf_lτ_i)+…+r_kL_k(τ_i)sin(2πf_lτ_i)(8)

其中，τ_i表示各采样点的时间坐标，i＝0,1,2…n；x(τ_i)为采样点对应的采样值；

步骤(II)、由式(8)得下式(9)所示的矩阵形式：

$x\left(\mathrm{\τ}\right)=L\left[\begin{array}{c}q\\ r\end{array}\right]---\left(9\right)$

其中，x(τ)＝[x(τ₀),x(τ₁),…,x(τ_N-1)]^T：动态信号模型各采样点对应的采样值构成的矩阵；x(τ_N-1)：采样点τ_N-1处的采样值；N＝1,2…n+1；q＝[q₀，q₁，…，q_N-1]^T与r＝[r₀，r₁，…，r_N-1]^T为拟合系数构成的矩阵，q_N-1和r_N-1为第N-1个时间坐标点对应的系数；[●]^T表示转置矩阵；辅助运算矩阵

$L=\left[\begin{array}{cc}{L}_0\left({\mathrm{\τ}}_0\right)\cos \left(2{\mathrm{\πf}}_l{\mathrm{\τ}}_0\right)...L_k\left({\mathrm{\τ}}_0\right)\cos \left(2{\mathrm{\πf}}_l{\mathrm{\τ}}_0\right),& L_0\left({\mathrm{\τ}}_0\right)\sin \left(2{\mathrm{\πf}}_l{\mathrm{\τ}}_0\right)...L_k\left({\mathrm{\τ}}_0\right)\sin \left(2{\mathrm{\πf}}_l{\mathrm{\τ}}_0\right)\\ L_0\left({\mathrm{\τ}}_1\right)\cos \left(2{\mathrm{\πf}}_l{\mathrm{\τ}}_1\right)...L_k\left({\mathrm{\τ}}_1\right)\cos \left(2{\mathrm{\πf}}_l{\mathrm{\τ}}_1\right),& L_0\left({\mathrm{\τ}}_1\right)\sin \left(2{\mathrm{\πf}}_l{\mathrm{\τ}}_1\right)...L_k\left({\mathrm{\τ}}_2\right)\sin \left(2{\mathrm{\πf}}_l{\mathrm{\τ}}_1\right)\\ .& .\\ .& .\\ .& .\\ L_0\left({\mathrm{\τ}}_{N-1}\right)\cos \left(2{\mathrm{\πf}}_l{\mathrm{\τ}}_{N-1}\right)...L_k\left({\mathrm{\τ}}_{N-1}\right)\cos \left(2{\mathrm{\πf}}_l{\mathrm{\τ}}_{N-1}\right),& L_0\left({\mathrm{\τ}}_{N-1}\right)\sin \left(2{\mathrm{\πf}}_l{\mathrm{\τ}}_{N-1}\right)...L_k\left({\mathrm{\τ}}_{N-1}\right)\sin \left(2{\mathrm{\πf}}_l{\mathrm{\τ}}_{N-1}\right)\end{array}\right]$

进一步的，所述步骤Ⅲ中，根据步骤II中得出的拟合系数矩阵计算τ_i时间点的动态同步相量包括下述步骤：

(1)将t_i时间点变换成下式所示的τ_i时间点：

τ_i＝(2t_i-a-b)/(b-a)；(10)

(2)按下式计算动态同步相量幅值X_m(τ_i)：

$X_m\left(t_i\right)=X_m\left({\mathrm{\τ}}_i\right)=\sqrt{q{\left({\mathrm{\τ}}_i\right)}^2+r{\left({\mathrm{\τ}}_i\right)}^2};---\left(11\right)$

(3)按下式计算动态同步相量初相角

其中，q(τ_i)和r(τ_i)为时间窗t∈[a,b]中任意时间点τ_i处的拟合系数；符号angle(a+jb)表示求取复数(a+jb)的角度；表示时间坐标在t_i处的动态同步相量相角；X_m(t_i)表示时间坐标在t_i处的动态同步相量幅值。

进一步的，用式(12)得出的同步相量相角计算下式所示瞬时频率f：

$f=\frac{1}{2\mathrm{\π}}\frac{d\mathrm{\φ}}{dt}\≈\frac{1}{2\mathrm{\π}}\frac{\mathrm{\Δ}\mathrm{\φ}}{\mathrm{\Δ}t}---\left(13\right)$

其中，为相量初相角偏移量；Δt为时间增量；

瞬时频率变化率R_f按下式计算：

$R_f=\frac{df}{dt}\≈\frac{\mathrm{\Δ}f}{\mathrm{\Δ}t}---\left(14\right)$

其中，Δf为频率偏移量。

与最接近的现有技术比，本发明具有以下有益效果：

1、本发明提供的同步相量的量测方法无论是输入静态信号还是动态信号，都可以准确快速地进行相量量测，其相量量测精度能够满足标准要求。

2、本发明提供的同步相量的量测方法可在频率偏移、电力系统振荡及失步时准确的量测相量。

附图说明

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

图1为本发明实施例一提供的仿真稳态50Hz测试中幅值误差示意图；

图2为本发明实施例一提供的仿真稳态50Hz测试中相角误差示意图；

图3为本发明实施例一提供的仿真稳态50Hz测试中频率误差示意图；

图4为本发明实施例一提供的仿真稳态50Hz测试中频率变化率误差示意图。

具体实施方式

以下将结合具体实施例和说明书附图对技术方案进行清楚、完整的描述。

本发明的实施例中，所述基于Legendre多项式拟合二段同步相量的量测方法基于动态信号数学模型获得。

电力系统中的动态信号模型表示为：

x(t)＝X_mcos[2πf(t)t+φ₀](1)

一个更为通用的电力系统瞬时波形的定义为：

$\begin{array}{c}x\left(t\right)=X_m\left(t\right)\cos \mathrm{\φ}\left(t\right)\\ =X_m\left(t\right)\cos [2\mathrm{\π}\∫f\left(t\right)\mathrm{dt}+{\mathrm{\φ}}_0]\end{array}---\left(2\right)$

其中，x(t)是信号的采样值，X_m(t)表示相量幅值的时间函数，f(t)表示信号频率的时间函数，φ₀是相量初相角。

式(2)中，波形的瞬时值表示为瞬时幅值与瞬时相角余弦的乘积，同步相量的相角测量基准为标准时钟信号，即同步相角φ(t)是相对于世界协调时间(UTC)同步的额定频率余弦信号的相角。

若将瞬时频率表示为f(t)＝f₀+g(t)，其中f₀是额定频率，g(t)表示实际瞬时频率与系统额定频率之差，则(2)式可表示为：

x(t)＝X_m(t)cos(2πf₀t+2π∫g(t)dt+φ₀)(3)

对应式(3)的瞬时同步相量定义为：

$\begin{array}{c}X\left(t\right)=X_m\left(t\right)/\sqrt{2}\exp j\mathrm{\φ}\left(t\right)\\ =X_m\left(t\right)/\sqrt{2}\exp j(2\mathrm{\π}\∫gdt+{\mathrm{\φ}}_0)\end{array}---\left(4\right)$

由上述分析可见，同步相量的定义是一个动态的概念，即理论上每个时间点处对应各自的动态同步相量。

本发明的技术方案是基于Legendre多项式拟合求取动态同步相量，以下式的Legendre方程表示式为例对所述Legendre多项式及其相关的重要性质进行说明：

$(1-t^2)\frac{d^{2}y}{{\mathrm{dt}}^2}-2t\frac{dy}{dt}+n(n+1)y=0---\left(5\right)$

Legendre多项式是Legendre方程的幂级数解，可表示为：

$L_n\left(t\right)=\frac{1}{2^{n}n!}\frac{d^n}{{\mathrm{dt}}^n}{(t^2-1)}^n:---\left(6\right)$

L_n(t)：Legendre方程的第n级幂级数解

Legendre多项式具有如下性质：

1、Legendre多项式的递推关系 $\begin{array}{c}{L}_0\left(t\right)=0;L_1\left(t\right)=t;\\ (n+1)L_{n+1}\left(t\right)=(2n+1)t{L}_n\left(t\right)-n{L}_{n-1}\left(t\right),n\≥2,t\∈[-\mathrm{1,1}]\end{array}---\left(7\right)$

式(7)中，下标n表示Legendre多项式阶数。

2、Legendre多项式的正交性

Legendre多项式在t∈[-1,1]区间上体现出正交性，即：

${\∫}_{-1}^1{L}_m\left(t\right)L_n\left(t\right)dt=0,m\≠n---\left(8\right)$

其中，下标m和n表示Legendre多项式的阶数。

上述式(8)积分可以通过线性变换转换到任意区间t∈[a,b]上的积分。

$\mathrm{\τ}=\frac{2t-a-b}{b-a}\⇔t=\frac{1}{2}\[(b-a)\mathrm{\τ}+a+b\]---\left(9\right)$

从而，对于任意函数f(t)，可得：

${\∫}_a^{b}f\left(t\right)dt={\∫}_{-1}^{1}f\left(\frac{(b-a)\mathrm{\τ}+a+b}{2}\right)\left(\frac{(b-a)}{2}\right)d\mathrm{\τ}---\left(10\right)$

利用式(8)、(9)和(10)，可以将Legendre多项式的正交性应用到任意区间上。

3、Legendre多项式的奇偶性

L_n(-t)＝(-1)ⁿL_n(t),t∈[-1,1](11)

利用Legendre多项式的奇偶性和三角函数的奇偶性，可以得到如下结论：

${\∫}_{-1}^1{L}_m\left(t\right)cos\left(\mathrm{\ω}t\right)L_n\left(t\right)cos\left(\mathrm{\ω}t\right)dt=0,\frac{m+n+1}{2}=k,k\∈N$

${\∫}_{-1}^1{L}_m\left(t\right)cos\left(\mathrm{\ω}t\right)L_n\left(t\right)sin\left(\mathrm{\ω}t\right)dt=0,\frac{m+n}{2}=k,k\∈N---\left(12\right)$

${\∫}_{-1}^1{L}_m\left(t\right)sin\left(\mathrm{\ω}t\right)L_n\left(t\right)sin\left(\mathrm{\ω}t\right)dt=0,\frac{m+n+1}{2}=k,k\∈N$

4、Legendre多项式的离散形式

由于实际物理系统中处理的信号多为离散信号，因此将上述Legendre多项式的性质转化为离散形式，便于实际应用。

由定积分的定义，可以将连续定义域中的积分转化为离散定义域中的求和形式：

$\begin{array}{c}{\∫}_{-1}^1{L}_m\left(t\right)L_n\left(t\right)dt=\underset{\mathrm{\δ}t\→0}{\lim }{\mathrm{\Σ}}_{i=0}^{N-1}{L}_m\left(t_i\right)L_n\left(t_i\right)\mathrm{\δ}t\\ \≈\frac{2}{N}{\mathrm{\Σ}}_{i=0}^{N-1}{L}_m\left(t_i\right)L_n\left(t_i\right)\end{array}---\left(13\right)$

式中，N表示对积分区间[-1,1]进行分段的段数，且δt·N＝2，δt为积分区间[-1,1]分段后的区间长度。因此Legendre多项式的正交性可以表示为：

${\&Integral;}_{-1}^1{L}_m\left(t\right)L_n\left(t\right)dt=0\&DoubleLeftRightArrow;<L_m\left(t\right),L_n\left(t\right)>=0,m\&NotEqual;n---\left(14\right)$

其中，表示向量α(t)和向量β(t)的内积。

将式(12)离散化：

${\mathrm{\&Sigma;}}_{i=0}^{N-1}{L}_m\left(t_i\right)cos\left({\mathrm{\&omega;t}}_i\right)L_n\left(t_i\right)cos\left({\mathrm{\&omega;t}}_i\right)=0,\frac{m+n+1}{2}=k,k\&Element;N$

${\mathrm{\&Sigma;}}_{i=0}^{N-1}{L}_m\left(t_i\right)cos\left({\mathrm{\&omega;t}}_i\right)L_n\left(t_i\right)sin\left({\mathrm{\&omega;t}}_i\right)=0,\frac{m+n}{2}=k,k\&Element;N---\left(15\right)$

${\mathrm{\&Sigma;}}_{i=0}^{N-1}{L}_m\left(t_i\right)sin\left({\mathrm{\&omega;t}}_i\right)L_n\left(t_i\right)sin\left({\mathrm{\&omega;t}}_i\right)=0,\frac{m+n+1}{2}=k,k\&Element;N$

将式(15)表示为向量内积的形式：

$<L_m\left(t\right)cos\left(\mathrm{\&omega;}t\right),L_n\left(t\right)cos\left(\mathrm{\&omega;}t\right)>=0,\frac{m+n+1}{2}=k,k\&Element;N$

$<L_m\left(t\right)cos\left(\mathrm{\&omega;}t\right),L_n\left(t\right)sin\left(\mathrm{\&omega;}t\right)>=0,\frac{m+n}{2}=k,k\&Element;N---\left(16\right)$

$<L_m\left(t\right)sin\left(\mathrm{\&omega;}t\right),L_n\left(t\right)sin\left(\mathrm{\&omega;}t\right)>=0,\frac{m+n+1}{2}=k,k\&Element;N$

即向量L_n(t)cos(ωt)和L_n(t)sin(ωt)的正交性由(m+n)的奇偶性决定。

下面对基于Legendre多项式拟合的同步相量的量测方法进行具体说明。

将动态信号模型x(t)表示为：

x(t)＝X_m(t)cos(2π∫fdt+φ₀)＝X_m(t)cos(2πf₀t+2π∫Δfdt+φ₀)(17)

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

记则式(17)可重新写为：

x(t)＝X_m(t)cos[2πf₀t+φ(t)](18)

对应t时刻的瞬时同步相量为：

$X\left(t\right)=\frac{X_m\left(t\right)}{\sqrt{2}}\exp \&lsqb;j\mathrm{\&phi;}\left(t\right)\&rsqb;---\left(19\right)$

对式(18)的电力系统瞬时信号模型x(t)展开得到：

x(t)＝X_m(t)cosφ(t)cos(2πf₀t)-X_m(t)sinφ(t)sin(2πf₀t)(20)

式(20)中的电力系统瞬时信号由基频信号cos(2πf₀t)和sin(2πf₀t)经加权后合成得到，其中权重函数X_m(t)cosφ(t)和X_m(t)sinφ(t)分别看作余弦信号和正弦信号的包络线。

令则利用k阶Legendre多项式将q(t)展开，并进行坐标变换如下：

q(τ)＝q₀L₀(τ)+q₁L₁(τ)+…+q_kL_k(τ)+o(τ^k)(21)

其中L_k(τ)：多项式的第k阶展开式；q_k：多项式第k阶展开式的系数；o(τ^k)：多项式展开式的皮亚诺型余项；k＝0，1，2…n；

令则利用k阶Legendre多项式将r(t)展开，并进行坐标变换如下：

r(τ)＝r₀L₀(τ)+r₁L₁(τ)+…+r_kL_k(τ)+o(τ^k)(22)

其中，L_k(τ)：多项式的第k阶展开式；r_k：多项式第k阶展开式的系数；o(τ^k)：多项式展开式的皮亚诺型余项；k＝0，1，2…n。

由于Legendre多项式只在自变量区间[-1,1]上有定义，故对数据窗内数据点的时间坐标t∈[a,b]进行线性坐标变换，即令τ＝(2t-a-b)/(b-a)。同时，为了保证时间坐标和电压采样值坐标的对应关系，应保证坐标变换前后波形形状、周期数和采样点数不变。

式(21)和(22)中的o(τ^k)为多项式展开的皮亚诺型余项，即拟合误差是τ^k项的高阶无穷小。

式(20)在τ时间坐标下仍成立，如下式所示：

$\begin{array}{c}x\left(\mathrm{\&tau;}\right)=X_m\left(\mathrm{\&tau;}\right)\cos \mathrm{\&phi;}\left(\mathrm{\&tau;}\right)\cos \left(2\mathrm{\&pi;}{f}_l\mathrm{\&tau;}\right)-X_m\left(\mathrm{\&tau;}\right)\sin \mathrm{\&phi;}\left(\mathrm{\&tau;}\right)\sin \left(2\mathrm{\&pi;}{f}_l\mathrm{\&tau;}\right)\\ =q\left(\mathrm{\&tau;}\right)\cos \left(2\mathrm{\&pi;}{f}_l\mathrm{\&tau;}\right)+r\left(\mathrm{\&tau;}\right)\sin \left(2\mathrm{\&pi;}{f}_l\mathrm{\&tau;}\right)\end{array}---\left(23\right)$

其中，f_l：t时间坐标下的额定频率f₀映射到τ时间坐标下的对应频率。

将式(21)和(22)代入式(23)，得到在[-1,1]区间上成立的以下表示式：

$\begin{array}{c}x\left(\mathrm{\&tau;}\right)=[q_0{L}_0\left(\mathrm{\&tau;}\right)+q_1{L}_1\left(\mathrm{\&tau;}\right)+...+q_k{L}_k\left(\mathrm{\&tau;}\right)]\cos \left(2\mathrm{\&pi;}{f}_l\mathrm{\&tau;}\right)\\ +[r_0{L}_0\left(\mathrm{\&tau;}\right)+r_1{L}_1\left(\mathrm{\&tau;}\right)+...+r_k{L}_k\left(\mathrm{\&tau;}\right)]\sin \left(2\mathrm{\&pi;}{f}_l\mathrm{\&tau;}\right)\\ =q_0{L}_0\left(\mathrm{\&tau;}\right)\cos \left(2\mathrm{\&pi;}{f}_l\mathrm{\&tau;}\right)+q_1{L}_1\left(\mathrm{\&tau;}\right)\cos \left(2\mathrm{\&pi;}{f}_l\mathrm{\&tau;}\right)+...+q_k{L}_k\left(\mathrm{\&tau;}\right)\cos \left(2\mathrm{\&pi;}{f}_l\mathrm{\&tau;}\right)\\ +r_0{L}_0\left(\mathrm{\&tau;}\right)\sin \left(2\mathrm{\&pi;}{f}_l\mathrm{\&tau;}\right)+r_1{L}_1\left(\mathrm{\&tau;}\right)\cos \left(2\mathrm{\&pi;}{f}_l\mathrm{\&tau;}\right)+...+r_k{L}_k\left(\mathrm{\&tau;}\right)\sin \left(2\mathrm{\&pi;}{f}_l\mathrm{\&tau;}\right)\end{array}---\left(24\right)$

将各采样点的时间坐标和采样值代入(24)式中：

x(τ_i)＝q₀L₀(τ_i)cos(2πf_lτ_i)+q₁L₁(τ_i)cos(2πf_lτ_i)+...+q_kL_k(τ_i)cos(2πf_lτ_i)

+r₀L₀(τ_i)sin(2πf_lτ_i)+r₁L₁(τ_i)cos(2πf_lτ)+...+r_kL_k(τ_i)sin(2πf_lτ_i)

(25)

整理为矩阵形式为：

$x\left(\mathrm{\&tau;}\right)=L\left[\begin{array}{c}q\\ r\end{array}\right]---\left(26\right)$

其中，x(τ)＝[x(τ₀),x(τ₁),…,x(τ_N-1)]^T，q＝[q₀,q₁…,q_N-1]^T，r＝[r₀,r₁…,r_N-1]^T，

$L=\left[\begin{array}{cc}{L}_0\left({\mathrm{\&tau;}}_0\right)\cos \left(2{\mathrm{\&pi;f}}_l{\mathrm{\&tau;}}_0\right)...L_k\left({\mathrm{\&tau;}}_0\right)\cos \left(2{\mathrm{\&pi;f}}_l{\mathrm{\&tau;}}_0\right),& L_0\left({\mathrm{\&tau;}}_0\right)\sin \left(2{\mathrm{\&pi;f}}_l{\mathrm{\&tau;}}_0\right)...L_k\left({\mathrm{\&tau;}}_0\right)\sin \left(2{\mathrm{\&pi;f}}_l{\mathrm{\&tau;}}_0\right)\\ L_0\left({\mathrm{\&tau;}}_1\right)\cos \left(2{\mathrm{\&pi;f}}_l{\mathrm{\&tau;}}_1\right)...L_k\left({\mathrm{\&tau;}}_1\right)\cos \left(2{\mathrm{\&pi;f}}_l{\mathrm{\&tau;}}_1\right),& L_0\left({\mathrm{\&tau;}}_1\right)\sin \left(2{\mathrm{\&pi;f}}_l{\mathrm{\&tau;}}_1\right)...L_k\left({\mathrm{\&tau;}}_2\right)\sin \left(2{\mathrm{\&pi;f}}_l{\mathrm{\&tau;}}_1\right)\\ .& .\\ .& .\\ .& .\\ L_0\left({\mathrm{\&tau;}}_{N-1}\right)\cos \left(2{\mathrm{\&pi;f}}_l{\mathrm{\&tau;}}_{N-1}\right)...L_k\left({\mathrm{\&tau;}}_{N-1}\right)\cos \left(2{\mathrm{\&pi;f}}_l{\mathrm{\&tau;}}_{N-1}\right),& L_0\left({\mathrm{\&tau;}}_{N-1}\right)\sin \left(2{\mathrm{\&pi;f}}_l{\mathrm{\&tau;}}_{N-1}\right)...L_k\left({\mathrm{\&tau;}}_{N-1}\right)\sin \left(2{\mathrm{\&pi;f}}_l{\mathrm{\&tau;}}_{N-1}\right)\end{array}\right]---\left(27\right)$

矩阵L的列向量具有部分正交性，矩阵L条件数较小，奇异值分解的效率更高，从而使得整个量测方法效率提高。

在计算得到拟合系数q＝[q₀,q₁…,q_N-1]^T和r＝[r₀,r₁…,r_N-1]^T后，利用式(21)和式(22)进一步得到时间窗t∈[a,b]中任意时间点t_i处的q(τ_i)、r(τ_i)，而后利用q(τ_i)和r(τ_i)计算τ_i时间点处所对应的动态同步相量，如下式(28)所示：

(1)由下式将τ_i时间点变换成t_i时间点：

$t_i=\frac{(b-a){\mathrm{\&tau;}}_i+a+b}{2};---(28-1)$

(2)由下式计算动态同步相量幅值X_m(τ_i)：

$X_m\left(t_i\right)=X_m\left({\mathrm{\&tau;}}_i\right)=\sqrt{q{\left({\mathrm{\&tau;}}_i\right)}^2+r{\left({\mathrm{\&tau;}}_i\right)}^2};---(28-2)$

(3)由下式计算动态同步相量初相角

符号angle(a+jb)表示求取复数(a+jb)的角度。

由于同步相量的相角范围是[-π,π]，因此这里不能利用arctan(a+jb)求取相量相角。

式(28)计算出的同步相量的时标也是打在窗头t₀，对应τ₀＝-1处最准确，即：

$\begin{array}{c}{X}_m\left(t_0\right)=X_m\left({\mathrm{\&tau;}}_0\right)=\sqrt{q{(-1)}^2+r{(-1)}^2}\\ \mathrm{\&phi;}\left(t_0\right)=\mathrm{\&phi;}\left({\mathrm{\&tau;}}_0\right)=\mathrm{angle}[q(-1)-\mathrm{jr}(-1)]\end{array}---\left(29\right)$

在求取同步相量相角之后，利用相角计算值的差分形式求取瞬时频率；

$f=\frac{1}{2\mathrm{\&pi;}}\frac{d\mathrm{\&phi;}}{\mathrm{\&omega;}}\&ap;\frac{1}{2\mathrm{\&pi;}}\frac{\mathrm{\&Delta;}\mathrm{\&phi;}}{\mathrm{\&Delta;}t}---\left(30\right)$

其中，为相量初相角偏移量；Δt为时间增量；

利用频率计算值进一步差分求取瞬时频率变化率。

$R_f=\frac{df}{dt}\&ap;\frac{\mathrm{\&Delta;}f}{\mathrm{\&Delta;}t}---\left(31\right)$

其中，Δf为频率偏移量。

实施例一

对上述量测方法进行仿真测试，具体来说，对8阶Legendre多项式拟合算法的测量精度进行仿真分析及对比，其中，仿真中数据窗长取为两个工频周波，即40ms；采样率为5000Hz。仿真结果如表1中所示：

表1基于Legendre多项式拟合的同步相量的测量精度

下面以稳态50Hz为例具体说明。

输入信号频率为50Hz的正弦信号，信号初相角为60°，采样频率为5000Hz，上传频率为500Hz，采用数据窗长为两个工频周波，即40ms，仿真时间长度为1s，仿真结果如图1、图2、图3和图4所示。

由图1仿真测试的幅值误差示意图中得到最大幅值误差为4.3×10^-8％；由图2的相角误差示意图得到最大相角误差为2.4×10-8°；由图3的频率误差示意图得到最大频率误差为4.1×10^-8Hz；由图4的频率变化率误差示意图得到最大频率变化率误差为1.7×10^{- 5}Hz/s。根据《电力系统实时动态监测系统技术规范》规定，最大幅值误差为0.2％，最大相角误差为0.2°，最大频率误差为0.002Hz，最大频率变化率误差为0.01Hz/s，由上述数据对比，并结合附图和同步相量的测量精度表数据，可知，本发明方法的量测精度远高于PMU标准要求。

本发明实施例提供的量测方法无论是输入静态信号还是动态信号，都可以准确快速地进行相量量测，其相量量测精度能够满足标准要求。该方法可在频率偏移、电力系统振荡及失步时准确地量测相量，具有较高的计算效率，对现代的硬件来说计算负担并不重，可在实际装置中实现。

最后应当说明的是：以上实施例仅用以说明本申请的技术方案而非对其保护范围的限制，尽管参照上述实施例对本申请进行了详细的说明，所属领域的普通技术人员应当理解：本领域技术人员阅读本申请后依然可以对申请的具体实施方式进行种种变更、修改或者等同替换，但这些变更、修改或者等同替换，均在申请待批的权利要求保护范围之内。

Claims (6)

1.一种同步相量的量测方法，其特征在于，所述量测方法包括以下步骤：

I、拟合动态信号；

II、计算拟合系数；

Ⅲ、确定瞬时同步相量。

2.根据权利要求1所述的量测方法，其特征在于，所述步骤I中，所述动态信号模型x(t)如下式所示：

其中，相量初相角；f₀：额定频率；Δf：频率偏移量；X_m(t)：相量幅值的时间函数；f(t)：信号频率的时间函数；

对式(1)所示的动态信号模型的拟合包括以下步骤：

(i)将式(1)变形为下式(2)：

其中，

(ii)式(2)所述的动态信号模型的展开式如下式所示：

(iii)用Legendre多项式展开式(3)中的权重函数和

令则q(t)的k阶Legendre多项式展开式如下：

q(τ)＝q₀L₀(τ)+q₁L₁(τ)+…+q_kL_k(τ)+o(τ^k)(4)

其中，L_k(τ)：多项式的第k阶项；q_k：多项式第k阶项系数；o(τ^k)：多项式展开式的皮亚诺型余项；k为多项式展开的阶数，k＝0，1，2…n；

令则r(t)的k阶Legendre多项式展开式如下：

r(τ)＝r₀L₀(τ)+r₁L₁(τ)+…+r_kL_k(τ)+o(τ^k)(5)

其中，L_k(τ)：多项式的第k阶项；r_k：多项式展开式中第k阶项系数；o(τ^k)：多项式展开式的皮亚诺型余项；k为多项式展开的阶数，k＝0，1，2…n。

3.根据权利要求2所述的量测方法，其特征在于，时间坐标t∈[a,b]数据窗内的数据点的线性坐标变换如下式所示：

τ＝(2t-a-b)/(b-a)；(6)

将步骤(iii)中的式(4)和式(5)所示的展开式代入步骤(ii)的式(3)中，得如下式(7)所示的τ时间坐标下成立的动态信号模型x(τ)：

x(τ)＝[q₀L₀(τ)+q₁L₁(τ)+…+q_kL_k(τ)]cos(2πf_lτ)+[r₀L₀(τ)+r₁L₁(τ)+…+

r_kL_k(τ)]sin(2πf_lτ)＝

q₀L₀(τ)cos(2πf_lτ)+q₁L₁(τ)cos(2πf_lτ)+…+q_kL_k(τ)cos(2πf_lτ)+

r₀L₀(τ)sin(2πf_lτ)+r₁L₁(τ)sin(2πf_lτ)+…+r_kL_k(τ)sin(2πf_lτ)(7)

其中，f_l：t时间坐标下的额定频率f₀映射到τ时间坐标下的对应频率。

4.根据权利要求1所述的量测方法，其特征在于，步骤II中所述拟合系数的计算包括：步骤(I)、将采样点的时间坐标和采样值代入式(7)，得式(8)所示的τ时间坐标下的动态信号模型x(τ_i)的表示式：

x(τ_i)＝q₀L₀(τ_i)cos(2πf_lτ_i)+q₁L₁(τ_i)cos(2πf_lτ_i)+…+q_kL_k(τ_i)cos(2πf_lτ_i)+

r₀L₀(τ_i)sin(2πf_lτ_i)+r₁L₁(τ_i)sin(2πf_lτ_i)+…+r_kL_k(τ_i)sin(2πf_lτ_i)(8)

其中，τ_i表示各采样点的时间坐标，i＝0,1,2…n；x(τ_i)为采样点对应的采样值；

步骤(II)、由式(8)得下式(9)所示的矩阵形式：

$x\left(\mathrm{\&tau;}\right)=L\left[\begin{array}{c}q\\ r\end{array}\right]---\left(9\right)$

其中，x(τ)＝[x(τ₀),x(τ₁),…,x(τ_N-1)]^T：动态信号模型各采样点对应的采样值构成的矩阵；x(τ_N-1)：采样点τ_N-1处的采样值；N＝1,2…n+1；q＝[q₀，q₁，…，q_N-1]^T与r＝[r₀，r₁，…，r_N-1]^T为拟合系数构成的矩阵，q_N-1和r_N-1为第N-1个时间坐标点对应的系数；[●]^T表示转置矩阵；辅助运算矩阵

5.根据权利要求1所述的量测方法，其特征在于，所述步骤Ⅲ中，根据步骤II中得出的拟合系数矩阵计算τ_i时间点的动态同步相量包括下述步骤：

(1)将t_i时间点变换成下式所示的τ_i时间点：

τ_i＝(2t_i-a-b)/(b-a)；(10)

(2)按下式计算动态同步相量幅值X_m(τ_i)：

$X_m\left(t_i\right)=X_m\left({\mathrm{\&tau;}}_i\right)=\sqrt{q{\left({\mathrm{\&tau;}}_i\right)}^2+r{\left({\mathrm{\&tau;}}_i\right)}^2};---\left(11\right)$

(3)按下式计算动态同步相量初相角

其中，q(τ_i)和r(τ_i)为时间窗t∈[a,b]中任意时间点τ_i处的拟合系数；符号angle(a+jb)表示求取复数(a+jb)的角度；表示时间坐标在t_i处的动态同步相量相角；X_m(t_i)表示时间坐标在t_i处的动态同步相量幅值。

6.根据权利要求5所述的量测方法，其特征在于，用式(12)得出的同步相量相角计算下式所示瞬时频率f：

$f=\frac{1}{2\mathrm{\&pi;}}\frac{d\mathrm{\&phi;}}{dt}\&ap;\frac{1}{2\mathrm{\&pi;}}\frac{\mathrm{\&Delta;}\mathrm{\&phi;}}{\mathrm{\&Delta;}t}---\left(13\right)$

其中，为相量初相角偏移量；Δt为时间增量；

瞬时频率变化率R_f按下式计算：

$R_f=\frac{df}{dt}\&ap;\frac{\mathrm{\&Delta;}f}{\mathrm{\&Delta;}t}---\left(14\right)$

其中，Δf为频率偏移量。