基于宽线性最小方差无失真响应的非平衡电力系统频率估计方法与流程

本发明属于电力系统技术领域，具体涉及一种非平衡电力系统频率估计方法。

背景技术：

在电力系统中，常用的频率估计采用基于相位的方法，但此类方法不适用于发生大幅度的频率震荡的情况。

实际电力系统中由于存在谐波干扰、噪声和非平衡电压等无法预料的情况，在这些情况下对电力系统的频率进行快速而又准确的估计具有很大的实用意义。

由于标准单相技术的局限性，特别是当所选相遭受电压下降或瞬变时会严重影响频率估计的准确性。考虑线间电压时，由于在三相系统中存在六个不同的单相电压，难以选择最具代表性的单相信号来充分描述系统频率。

比较好的解决方案是设计一个同时考虑所有三相电压的框架，当任何相位遭受骤降，瞬变或谐波时可以对增强的鲁棒性作出统一估计。为此，clarke变换根据所有三相电压提供的信息构建了复信号，此变换使经典单相方法具有增强的鲁棒性，并且在复域中开发了许多已被证明比在真实域中的相应方法更可靠的解决方案。这些解决方案包括使用锁相环(pll)、最小二乘法、卡尔曼滤波和基于解调的方法。但是，在非平衡系统中这些解决方案往往不能满足实际需求。

技术实现要素：

发明目的：针对现有技术中存在的问题，本发明公开了一种基于宽线性最小方差无失真响应的非平衡电力系统频率估计方法，该方法更适用于非平衡系统。

技术方案：一种基于宽线性最小方差无失真响应的非平衡电力系统频率估计方法，包括如下步骤：

步骤1、建立三相电力系统的复电压信号：

其中va、vb、vc分别为三相电压信号的峰值电压，△t为采样间隔，等于1/fs，fs为采样频率，φ为初始相位，ω＝2πf0为电压信号的角频率，f0为系统频率；

步骤2、建立非平衡电力系统电压信号的wl-mvdr谱：

其中f＝[1,0]^t，vaug(k)＝[v^t(k),v^h(k)]^t为复电压信号v(k)的增广向量，为复电压信号v(k)的增广协方差矩阵，其中r为协方差矩阵，c为伪协方差矩阵，对求逆得到矩阵g和d，为激励矢量，s＝[1,e^-jω,…,e^-j(k-1)ω]^t，k为协方差矩阵r的维数；

符号t表示向量或矩阵的转置；h表示矩阵的共轭转置；符号*表示矩阵的伴随矩阵；符号-1表示矩阵的逆矩阵；

步骤3、求解使代价函数最小的角频率为系统角频率的估计值

所述代价函数j(k)为：

系统频率的估计值为：其中fs为采样频率。

优选地，步骤3中采用自适应算法计算系统角频率的估计值步骤如下：

(2.1)确定迭代初始值ω(0)，构建角频率ω更新模型：

ω(k+1)＝ω(k)-μ(coe1*sin(2ω(k))+coe2*cos(2ω(k))+coe3*sin(ω(k))+coe4*cos(ω(k)))

其中ω(k+1)表示第k次迭代估计的角频率，μ为迭代步长，

re(·)表示取括号内数值的实部，im(·)表示取括号内数值的虚部；

f(ω)＝s^hgs，g(ω)＝s^td*s，

(2.2)迭代计算步骤(2.1)中的ω(k+1)，检查本次迭代得到的角频率是否满足收敛条件，如满足，则ω(k+1)即为系统角频率的估计值

收敛条件为：|ω(k+1)-ω(k)|≤ε或迭代次数k达到预设迭代次数阈值l；ε是一个极小的正实数。

作为另一种优选，步骤3中采用直接求解法计算系统角频率的估计值步骤如下：

(3.1)建立关于角频率的一元四次方程：

其中

re(·)表示取括号内数值的实部，im(·)表示取括号内数值的虚部；

f(ω)＝s^hgs，g(ω)＝s^td^*s，

(3.2)求解步骤(3.1)中的一元四次方程，得到4个解，选择其中最接近系统角频率期望ω0＝2πf0/fs的作为系统角频率的估计值f0为系统频率期望；fs为采样频率。

上述计算中，矩阵g和d由增广协方差矩阵求逆计算得到，具体步骤为：

(4.1)根据复电压信号的观测值，建立数据阵v(k)：

其中m为用于计算的有效观测长度；n为总观测长度；

(4.2)增广协方差矩阵的逆可以由矩阵求逆引理得到：

(4.3)迭代计算步骤(4.2)中的检查相邻两次迭代得到的是否满足收敛条件，如满足，则即为增广协方差矩阵的逆矩阵；

(4.4)根据可得到协方差矩阵g和伪协方差矩阵d。

为了提高频率估计的精度，计算出矩阵g和d后，对g和d进行系数补偿：

g＝c1g′，d＝c2d′

其中c1和c2为补偿系数，g′和d′为补偿前的矩阵，g和d为补偿后的矩阵。

优选地，c1＝2.859，c2＝1。

有益效果：与现有的技术相比，本发明具有以下优点：1.充分利用了三项电压完整的二阶信息，增强了频率估计的鲁棒性。2.与传统的线性自适应估计相比，该方法更适用于非平衡系统并且给出了无偏的频率估计。3.该模型对于三相电压振幅随时间的推移和高次谐波存在的变化也不敏感。4.充分考虑输入信号与输出信号的噪声，鲁棒抗噪性能好。5.算法稳定，计算复杂度低，收敛性能，稳态收敛精度均明显高于同类其他算法。

附图说明

图1为typeb型电压暂降的非平衡系统下基于mvdr算法与基于wlmvdr算法的频率估计对比图；

图2为步长与迭代的wl-mvdr算法估计均方误差之间的关系；

图3为typeb型电压暂降的非平衡系统下，信噪比snr从20db到50db进行100次独立重复试验，基于mvdr算法与基于wlmvdr算法的频率估计偏差对比图；

图4typeb型电压暂降的非平衡系统下，信噪比snr从20db到50db进行100次独立重复试验，mvdr算法与wlmvdr算法的方差以及克拉美劳下界；

图5为当电压幅度调制情况下两种算法的频率估计；

图6为当电压频率调制情况下两种算法的频率估计；

图7为当添加谐波分量造成不平衡情况的频率估计；

图8为当va(k)变化造成真实数据不平衡情况的频率估计图；8(a)为真实数据相位a电压变化的三相电压图；8(b)为真实数据相位a电压变化时频率估计图。

具体实施方式

下面结合附图和具体实施方式，进一步阐明本发明。

一种基于宽线性最小方差无失真响应的非平衡电力系统频率估计方法，包括如下步骤：

步骤1、建立三相电力系统的复电压信号；

对三相电力系统中的电压信号进行采样，得到三相电力系统电压模型，可表示成下面的离散形式：

va(k)＝vacos(ωkδt+φ)z(1)

其中，va、vb、vc为峰值电压，△t为采样间隔，等于1/fs，这里fs为采样频率，φ为初始相位，ω＝2πf0为电压信号的角频率。

将三相电压通过克拉克变换变换到二相静止坐标系上，即与时间有关的三相电压可以通过正交克拉克变换得到一个零序列v0和两个相互垂直的分量vα、vβ，即：

其中，是用来确保变换中系统功率不变。

在没有噪声的平衡系统条件下，va、vb、vc是相等的，即：

其中幅度值是一个常数，vα(k)和vβ(k)是位置随时间变化的速率正比于系统频率的一个点的纵坐标。在实际中，通常只用vα和vβ部分来组成复电压v(k)，即：

这里ω0＝ωδt，为归一化角频率，vα(k)和vβ(k)表示某一时刻电压的正交表示。对于给定的采样频率fs，v(k)的概率密度函数是旋转不变的，因为对于任意正的θ，v和ve^θ有着相同的分布。这意味着v(k)是一个二阶圆信号，vα(k)和vβ(k)的能量相等，所以协方差矩阵r＝e[v(k)v(k)^h]能充分描述其二阶统计特性，并且伪协方差矩阵c＝e[v(k)v(k)^t]＝0。

当三相电力系统偏离正常状态，例如存在不同程度的凹陷或瞬变，va、vb、vc变得不相等，则复电压信号为：

其中va、vb、vc分别为三相电压信号的峰值电压，△t为采样间隔，等于1/fs，fs为采样频率，φ为初始相位，ω＝2πf0为电压信号的角频率，f0为系统频率。

步骤2、建立非平衡电力系统电压信号的wl-mvdr谱；

将复电压信号v(k)分别通过两个p阶复数fir带通滤波器后再合并为y(n)，则复电压信号的功率为y(n)的功率除以滤波器带宽。设两个带通滤波器的系数为：gi＝[gi(0),gi(1),…gi(p)]^h，i∈{1,2}；则：

其中，vaug(k)＝[v^t(k),v^h(k)]^t，

在约束条件

g^hs1＝1g^hs2＝0(9)

下，使得输出功率最小的g的解为：

其中，s＝[s1,s2]，s1＝[s^t,0^t]^t，s2＝[0^t,s^h]^t，

s＝[1,e^-jω,…,e^-j(k-1)ω]^t(11)

为增广协方差矩阵。

其中r为协方差矩阵，c为伪协方差矩阵。k为协方差矩阵r的维数。

对于非平衡系统，当v(k)由式(7)出时，可得：

同理可得：

±[v(k-1)v^*(k-1)]＝|a|²+|b|²(17)

所以，

同理，

可知协方差矩阵r为共轭对称阵，矩阵c为对称阵。

将代入输出功率p，可得：

其中，f＝[1,0]^t,

由数学矩阵求逆规律可知，可表示为:

且矩阵g为共轭对称阵，d为对称阵。代入输出功率p表达式可得wl-mvdr谱为：

其中f＝[1,0]^t，vaug(t)＝[v^t(k),vh(k)]^t为复电压信号v(k)的增广向量，为复电压信号v(k)的增广协方差矩阵，其中r为协方差矩阵，c为伪协方差矩阵，对求逆得到s为激励矢量，s＝[1,e^-jω,…,e^-j(k-1)ω]^t，k为协方差矩阵r的维数；

步骤3、求解使代价函数最小的角频率为系统角频率的估计值

wl-mvdr谱在期望频率f0处表现出一个峰值。

当s的基础频率ω和输入信号中所包含的频率分量(由协方差矩阵表示)一样的时候，wl-mvdr估计算法返回该频率分量处的能量。然而，在有噪声的情况下，假设每个频率分量处的能量都大于噪声的能量，pwlmvdr的倒数在每一个频率分量处有一个极小值。理论上，每一个极小值所在的位置都和实际存在的频率分量相对应，其中，绝对值最小的极小值所在的位置为基频分量。当只存在基频和白噪声的时候，协方差矩阵r的维数k＝2。

wl-mvdr频率估计基于最小化代价函数j(k)：

可以采用自适应算法或直接求解法计算系统角频率的估计值

采用自适应算法计算系统角频率的估计值的步骤如下：

(2.1)确定迭代初始值ω(0)，构建角频率ω更新模型；

使用lms(least-mean-square，最小均方)梯度：

其中μ为迭代步长，用来权衡收敛速度和估计精度。设：

f(ω)＝s^h5sg(ω)＝s^td^*s

于是，

所以：

而

所以，

其中，

角频率ω更新模型为：

ω(k+1)＝ω(k)-μ(coe1*sin(2ω(k))+coe2*cos(2ω(k))+coe3*sin(ω(k))+coe4*cos(ω(k)))(35)

其中ω(k+1)表示第k次迭代估计的角频率，re(·)表示取括号内数值的实部，im(·)表示取括号内数值的虚部；

(2.2)迭代计算步骤(2.1)中的ω(k+1)，检查本次迭代得到的角频率是否满足收敛条件，如满足，则ω(k+1)即为系统角频率的估计值

收敛条件为：|ω(k+1)-ω(k)|≤ε或迭代次数k达到预设迭代次数阈值l；ε是一个极小的正实数。

采用直接求解法计算系统角频率的估计值步骤如下：

(3.1)建立关于角频率的一元四次方程；

在稳定的状态下，式(35)收敛，梯度趋近于0，即：

coe1*sin(2ω)+coe2*cos(2ω)+coe3*sin(ω)+coe4*cos(ω)≈0(36)

由万能公式：

以及二倍角公式：

sin(2ω)＝2sin(ω)cos(ω)cos(2ω)＝2cos²(ω)-1

式(35)可以化为：

其中coe1～coe4分别为式(31)～式(34)取ω＝ω0时的值，即：

(3.2)求解步骤(3.1)中关于的一元四次方程，得到的解x后，利用即可得到系统角频率的估计值

由于一元四次方程有4个解，选择使系统角频率的估计值最接近系统角频率期望ω0＝2πf0/fs的作为最终的系统角频率的估计值f0为系统频率期望；fs为采样频率。

计算出后，系统频率的估计值为：其中fs为采样频率。

上述计算中，矩阵g和d由增广协方差矩阵求逆计算得到，具体步骤为：

(4.1)根据复电压信号的观测值，建立数据阵v(k)：

其中m为用于计算的有效观测长度；n为总观测长度；

(4.2)增广协方差矩阵的逆可以由矩阵求逆引理得到：

(4.3)迭代计算步骤(4.2)中的检查相邻两次迭代得到的是否满足收敛条件，如满足，则即为增广协方差矩阵的逆矩阵；

(4.4)根据可得到协方差矩阵g和伪协方差矩阵d。

由于对矩阵求逆的计算过程会影响到结果的精度，以γ＝0.7的typeb型非平衡电压系统为例，采样频率为1000hz，将计算出来的与相乘，结果为：

计算结果与单位矩阵之间存在较大的误差。为了解决这个问题，一种办法是寻求更加精确的求逆算法，另一种办法是对现有求逆算法进行补偿，使得求逆后的矩阵与原矩阵的乘积尽可能地接近于单位阵i，即：

令g＝c1g′，d＝c2d′，其中c1和c2为补偿系数，g′和d′为补偿前的矩阵，g和d为补偿后的矩阵。用补偿后的矩阵g和d进行计算，即：

r(k)g+c(k)d^*＝i(40)

c^*(k)g+r^*(k)d^*＝o(41)

r(k)d+c(k)5^*＝o(42)

c^*(k)d+r^*(k)g^*＝i(43)

经过计算，c1＝2.859，c2＝1时计算结果的精度较高。

实施例1：

将估计频率初始化为49.5hz，系统实际频率为f0＝50hz，迭代的wl-mvdr算法的迭代步长设为μ＝0.1，估计结果如图1所示。其中f-mvdr表示用mvdr算法进行的频率估计，mvdr表示的是没有迭代的mvdr频率估计算法。由图1可知，迭代的wl-mvdr算法都能精确地估计系统频率，而mvdr算法有一定的偏差。

接下来考虑在有噪声的情况下，两种算法的估计偏差和方差的比较。

首先考察步长μ对迭代的wl-mvdr算法的影响。步长在自适应算法中对于最小化代价函数有着很重要的作用，步长越小，算法越稳定，收敛得越慢；反之，步长越大，算法稳态估计误差越大，收敛得越快。图2表明在信噪比snr＝50db电压为typeb型下陷的情况下，不同的步长μ对迭代的wl-mvdr估计算法性能的影响，此处性能用均方误差(mse)表示。由图2可知，步长越小，均方误差越小，估计越精确。

接下来，通过与估计方差的理论克拉美-劳(crlb)下界作比较，分析这两种算法的估计方差。图3和图4表示的是不同的信噪比情况下，μ＝0.1时，100次独立重复仿真后几种算法的估计偏差和方差的比较。其中f-mvdr表示用mvdr算法进行的频率估计，f-wlmvdr表示用wl-mvdr算法进行的频率估计。结果表明，在信噪比很大的情况下，迭代的wl-mvdr算法表现出了无偏性，而mvdr算法有一定的偏差。这与理论相符。在方差分析图中，wl-mvdr算法比i-mvdr算法更接近于crlb下界，估计性能较好。

然后讨论当typeb型下陷的电力系统中调幅的情况下的算法性能，(无噪声)即：

va(k)＝1+0.05sin(2πkδt)

va(k)＝1+0.1sin(2πkδt)

va(k)＝1+0.15sin(2πkδt)

仿真结果如图5，该图表明，在一点调幅的情况下wl-mvdr表现出了很小的估计误差，而mvdr算法振荡相对比较大。

snr＝50db的条件下，假设系统频率在某一时刻突然改变，考察集中算法的跟踪性能及反应。即设：

图6表明两种算法对于频率变化的跟踪情况，由该图可以看出，迭代的wl-mvdr算法与f0最为接近，跟踪效果最好，反应最快；i-mvdr算法估计误差较大。

最后考虑高阶谐波对集中估计算法的性能影响。仍然取typeb型电压下陷，snr＝50db加上20％的3阶谐波、10％的5阶谐波以及10％的7阶谐波，加入了谐波之后的电压和电压跟踪结果如图7，可以清楚地看出迭代的wl-mvdr算法在高阶谐波的干扰下估计更加稳定精确。

实施例2：

将该方案用于真实环境中进行频率估计。实际中，三相电压由一个110/20/10kv的变压站记录。由abb公司生产的rel531数字线路距离保护终端被安装在站内，用于监测三个“相-地”电压的变化。该设备被用来为记录每当相电压值下降到其正常值的90％以上超过20ms的情况。在1khz下对系统频率为50hz的测量的三个“相—地”电压进行采用，如图8(a)所示，在t＝0.06s附近，相位va与地短路，电压下降到其正常值的44％。同时，相位vb和vc电压分别骤升36％，50％，得到的非圆度γ＝0.0448。

wl-mvdr算法与mvdr算法频率估计能力分别如图8(b)所示，wl-mvdr算法在正常操作条件下提供了准确的响应，然而，mvdr算法处理不平衡的能力较差。