电子式互感器中采样率转换的低延迟滤波器设计方法

摘要

本发明公开了一种电子式互感器中采样率转换的低延迟滤波器设计方法。该方法首先采用内插器和抽取器的级联来实现任意分数倍采样频率的转换，将内插器中的抗镜像滤波器和抽取器中的抗混叠滤波器合并为一个低通滤波器；然后采用均方误差最小化准则来求解该滤波器的系数向量，以滤波器通带幅值与阻带幅值为约束条件，以均方误差最小化为优化目标；最后将基于约束最小二乘法设计的滤波器系数向量的求解过程转化为求解一个正定二次规划问题，即可直接求解得到滤波器系数向量。该方法解决了直接线性卷积滤波过程带来较大输出延迟，滤波器阶数越高群延迟也越大的问题，极大地提高了保护动作的快速性。

说明书

技术领域

本发明涉及一种互感器中采样率转换的低延迟滤波器设计方法，基于约束最小二乘 法设计，属于FIR滤波器技术领域。

背景技术

在智能变电站中，微机保护系统过程层采用电子式互感器，间隔层二次保护装置 与过程层的连接通过网络接口实现，保护装置从通信接口中接收由不同数字源A/D采 样后的数字量，输入信号由传统的单一模拟信号变成了多数字源采样信号的混合输 入。基于差动保护原理的保护装置通常需要获得瞬时差电流值，要求与被保护一次设 备相连接的所有支路的电流采样值的采样频率必须一致。当各支路的电子式电流互感 器输出数据的采样频率不完全一致时，须将采样频率进行归一化。

FIR滤波器主要采用非递归结构，在Z平面上不存在极点，因此在继电保护算法 上广泛采用的是FIR数字滤波算法。但FIR滤波器的选频特性较差，实现与IIR滤波 器同等幅频响应特性，其阶数通常比IIR高数倍，反映在滤波器输出上，对应着群延 迟较大，输出响应速度较慢，这是FIR滤波器的主要问题所在。

传统的FRR滤波方法是先将输入序列按顺序滤波，然后将所得结果逆转后反向通 过滤波器，再将所得结果逆转后输出，即得到精确零相位失真的输出序列。该方案在 滤波环节进行了两次滤波，即在时域进行两次卷积计算，而每一次卷积滤波都会带来 群延迟，两次滤波其群延迟增加一倍，使得在采用高阶FIR滤波器时，输出严重滞后 于输入序列。尽管输出和输入在相位上不存任何相位上的失真，但引入了较大的输出 延迟，不利于保护的快速动作。

因此，发明一种性能更为优越、应用范围更为广泛的互感器采样率转换的新算法 成为亟需解决的课题。

发明内容

发明目的：针对上述问题，本发明提供了一种互感器中采样率转换的低延迟滤波 器设计方法，致力于解决数字低通FIR滤波器的群延迟较大的问题，实现快速重采样。

技术方案：为达到上述目的，本发明采取的技术方案为：

一种电子式互感器中采样率转换的低延迟滤波器设计方法，包括如下步骤：

1）采用内插器和抽取器的级联来实现任意分数倍采样频率的转换，将内插器中的 抗镜像滤波器和抽取器中的抗混叠滤波器合并，建立一个低通滤波器模型；

2）采用均方误差最小化准则来求解该滤波器的系数向量，以滤波器通带幅值与阻 带幅值为约束条件，以均方误差最小化为优化目标；

3）将通带幅值约束条件与阻带幅值约束条件转化为约束关系函数；

4）将基于约束最小二乘法设计的滤波器系数向量求解过程转化为求解一个正定二 次规划问题，即可直接进行求解得到滤波器系数向量。

步骤2中采用均方误差最小化准则来求解滤波器系数向量b，其加权误差函数为：

$E\left(x\right)=\frac{1}{2}\underset{-\pi }{\overset{\pi }{\int }}W\left(\omega \right){|H_D\left(\omega \right)-H\left(\omega \right)|}^2\mathrm{d\omega }$

式中x=b，H_D(ω)为理想幅频响应，H(ω)为幅频响应，W(ω)为误差加权函数。

步骤3中将通带幅值约束条件转化为约束关系函数的方法为：

1）将通带幅值约束条件||H(ω)|-1|≤δ_p，化简为 1-δ_p≤|H(ω)|≤1+δ_p，δ_p表示通带允许的偏差，P表示通带内的一系列点集；

2）将H(ω)经过K次迭代得到H^(k)(ω)，其通带群延迟为K，H^(k)(ω)与H(ω)的关 系如下：

H(ω)≈e^-jKω|H^(k)(ω)|

$\left|H^{\left(k\right)}\left(\omega \right)\right|\approx c_k^T\left(\omega \right)b_k$

式中c_k(ω)=[cos(Kω),cos((K-1)ω),...cos((K-n)ω)]^T。

通带群延迟的误差满足以下条件

|τ_g(ω)-K|≤δ_g，ω∈P={ω_i^P,i=1,...n_P}

式中，τ_g(ω)为通带群延迟的误差，δ_g通带群延迟允许的偏差。

3）由以上式子得通带幅值约束关系式如下

$A_p^{\left(k\right)}{b}_k\le q_p^{\left(k\right)}$

式中

$A_p^{\left(k\right)}={\left(\begin{array}{c}{c}_k^T\left({\omega }_1^{\left(p\right)}\right)\\ ...\\ c_k^T\left({\omega }_{n_p}^{\left(p\right)}\right)\\ -c_k^T\left({\omega }_1^{\left(p\right)}\right)\\ ...\\ -c_k^T\left({\omega }_{n_p}^{\left(p\right)}\right)\end{array}\right)}_{2n_p\times (n+1)},$$q_p^{\left(k\right)}=\left(\begin{array}{c}1+{\delta }_p\\ ...\\ 1+{\delta }_p\\ {\delta }_p-1\\ ...\\ {\delta }_p-1\end{array}\right);$

步骤3中将阻带幅值约束条件转化为约束关系函数的方法为：

1）定义阻带幅值约束公式|H(ω)|≤δ_s，ω∈S={ω_i^S,i=1,...n_S}中

s(ω)=[0,sin(ω),...,sin(nω)]^T

c(ω)=[1,cos(ω),...,cos(nω)]^T

式中δ_s表示阻带允许的偏差，S表示阻带内的一系列点集；

2）令

$\left|c^T\left(\omega \right)b_k\right|\le \frac{{\delta }_s}{2}$

$\left|s^T\left(\omega \right)b_k\right|\le \frac{{\delta }_s}{2}$

则

|H^(k)(ω)|≤δ_s

3）将幅值约束简化为如下形式：

$A_s^{\left(k\right)}{b}_k\le q_s^{\left(k\right)}$

式中

$A_s^{\left(k\right)}={\left(\begin{array}{c}C\left(\omega \right)\\ -C\left(\omega \right)\\ S\left(\omega \right)\\ -S\left(\omega \right)\end{array}\right)}_{4n_s\times (n+1)},$$q_s^{\left(k\right)}=\frac{{\delta }_s}{2}{\left(\begin{array}{c}1\\ ...\\ 1\end{array}\right)}_{4n_s\times 1}$

$C\left(\omega \right)={\left(\begin{array}{c}{c}^T\left({\omega }_1^{\left(s\right)}\right)\\ ...\\ c^T\left({\omega }_{n_s}^{\left(s\right)}\right)\end{array}\right)}_{n_s\times (n+1)},$$S\left(\omega \right)={\left(\begin{array}{c}{s}^T\left({\omega }_1^{\left(s\right)}\right)\\ ...\\ s^T\left({\omega }_{n_s}^{\left(s\right)}\right)\end{array}\right)}_{n_s\times (n+1)}$

步骤4中将约束最小二乘法设计转化为解一个正定二次规划问题后进行求解的方法 为：

1）将基于约束最小二乘法设计滤波器的目标函数：

转化为解一个正定二次规划问题：

$\underset{b_k}{\min }\frac{1}{2}{b}_k^T{H}_k{b}_k+b_k^T{p}_k$

式中ω∈[0,π]={ω_i,i=1,...M}，T代表转置，H_k为(n+1)×(n+1)阶正定Hessian矩阵，其中：

2）将通带和阻带约束统一表示为如下形式

A_kb_k≤q_k

式中$A_k=\left(\begin{array}{c}{A}_p^{\left(k\right)}\\ A_s^{\left(k\right)}\end{array}\right),$$q_k=\left(\begin{array}{c}{q}_p^{\left(k\right)}\\ q_s^{\left(k\right)}\end{array}\right)$

3）定义二次规划的数学模型为

$\underset{x}{\min }\frac{1}{2}\mathrm{Hx}+f^{T}x$

其中约束条件为

Ax≤b

A_eqx=b_eq

lb≤x≤ub

式中，H为二次型矩阵，A、Aeq分别为不等式和等式约束的系数矩阵，其余为向 量；

4）利用Matlab调用quadprog函数命令进行直接求解。

本发明的有益效果：本发明提出的上述互感器中采样率转换的低延迟滤波器设计方 法，极大地提高了保护动作的快速性，同时在区外故障时转换后的误差电流能够满足保 护精度的要求，不会造成区外故障时的误动作，有效地实现了测量环节和智能保护装置 之间的无缝连接。该方法可有效弥补传统的采样频率频率转换方法的缺陷，即直接线性 卷积滤波过程带来了较大的输出延迟，以及滤波器阶数越高，群延迟也越大。

附图说明

图1本发明的流程框图；

图2低通滤波器分别实现抽取和插值示意图；

图3按有理因子U/D的采样频率转换方法示意图；

图4低通滤波器的频谱示意图；

图5500kV系统仿真模型图；

图6算例系统接线图；

图7a低延迟FIR滤波器的幅频响应图；

图7b低延迟FIR滤波器的相频响应图；

图7c低延迟FIR滤波器的群延迟图；

图7d低延迟FIR滤波器的幅值误差图；

图8母线N区外故障时两侧电流采样值图；

图9ECT2转换后电流采样数据图；

图10差动保护动作图；

图11区外故障时差电流图；

图12方法2转换后电流的频域分析图。

具体实施方式

以下结合附图和实施例对本发明作进一步说明。

如图1所示，一种电子式互感器中采样率转换的低延迟滤波器设计方法，包括如下 步骤：

1）在采样频率转换过程中，变换因子是任意的有理数U/D，可以采用抽取和插值 的级联来实现任意分数倍采样频率的转换，实现过程如图2所示。由于内插器中的抗镜 像滤波器和抽取器中的抗混叠滤波器均按相同的采样频率U_fs操作，则可将两者合并成 一个低通滤波器，如图3所示。

组合滤波器h(n)的理想低通频率响应特性为：

$H\left(e^{\mathrm{j\omega }}\right)=\left(\begin{array}{cc}U& 0\le \left|\omega \right|\le \min (\frac{\pi }{U},\frac{\pi }{D})\\ 0& \min (\frac{\pi }{U},\frac{\pi }{D})\le \left|\omega \right|\le \pi \end{array}\right)$

式中，π/U和π/D分别是抗镜像滤波器和抗混叠滤波器的截止频率，组合滤波器的 截止频率应取二者中的最小值。

2）采用均方误差最小化准则来求解滤波器系数向量，列出加权误差函数最小化的 约束条件，即通带幅值约束条件与阻带幅值约束条件。

图4给出了低通滤波器的幅频响应。其中，ω_p和ω_s分别表示通带和阻带的归一化 截止频率，δ_p和δ_s分别表示通带和阻带允许的偏差。

采用均方误差最小化准则来求解滤波器系数向量，其加权误差函数为

$E\left(x\right)=\frac{1}{2}\underset{-\pi }{\overset{\pi }{\int }}W\left(\omega \right){|H_D\left(\omega \right)-H\left(\omega \right)|}^2\mathrm{d\omega }$

式中，x=b，H_D(ω)为理想幅频响应，其中W(ω)为误差加权函数，它在不同 的频带中取值可以不同，逼近精度高的频带，W(ω)取值大；逼近精度低的频带，W(ω) 取值小。加权误差函数最小化的约束条件如下：

通带幅值约束条件：

$\left|\right|H\left(\omega \right)|-1|\le {\delta }_p,{\omega \in P=\{{\omega }_i}^P,i=1,...n_P\}$

式中，P表示通带内的一系列点集；

阻带幅值约束条件：

|H(ω)|≤δ_s，ω∈S={ω_i^S,i=1,...n_S}

式中，S表示阻带内的一系列点集。

3）通带幅值约束条件可以化简为

1-δ_p≤|H(ω)|≤1+δ_p

考虑到通带幅值约束条件，经过K次迭代得到H^(k)(ω)的通带群延迟为K，与H(ω) 的关系如下

H(ω)≈e^-jKω|H^(k)(ω)|

$\left|H^{\left(k\right)}\left(\omega \right)\right|\approx c_k^T\left(\omega \right)b_k$

式中

c_k(ω)=[cos(Kω),cos((K-1)ω),...cos((K-n)ω)]^T

必须考虑通带群延迟的误差满足以下条件

|τ_g(ω)-K|≤δ_g，ω∈P={ω_i^P,i=1,...n_P}

由K与H(ω)的关系式得通带幅值约束关系式如下

$A_p^{\left(k\right)}{b}_k\le q_p^{\left(k\right)}$

式中

$A_p^{\left(k\right)}={\left(\begin{array}{c}{c}_k^T\left({\omega }_1^{\left(p\right)}\right)\\ ...\\ c_k^T\left({\omega }_{n_p}^{\left(p\right)}\right)\\ -c_k^T\left({\omega }_1^{\left(p\right)}\right)\\ ...\\ -c_k^T\left({\omega }_{n_p}^{\left(p\right)}\right)\end{array}\right)}_{2n_p\times (n+1)},$$q_p^{\left(k\right)}=\left(\begin{array}{c}1+{\delta }_p\\ ...\\ 1+{\delta }_p\\ {\delta }_p-1\\ ...\\ {\delta }_p-1\end{array}\right)$

当k→∞时，H^(k)(ω)→H(ω)，此时||H(ω)|-1|≤δ_p与等波纹设计一致。

4）阻带幅值误差约束如下式

|H(ω)|≤δ_s，ω_s≤ω≤π

定义

s(ω)=[0,sin(ω),...,sin(nω)]^T

c(ω)=[1,cos(ω),...,cos(nω)]^T

如果

$\left|c^T\left(\omega \right)b_k\right|\le \frac{{\delta }_s}{2}$

$\left|s^T\left(\omega \right)b_k\right|\le \frac{{\delta }_s}{2}$

则满足阻带幅值约束条件

|H^(k)(ω)|≤δ_s

因此幅值约束可以写成如下形式

$A_s^{\left(k\right)}{b}_k\le q_s^{\left(k\right)}$

式中

$A_s^{\left(k\right)}={\left(\begin{array}{c}C\left(\omega \right)\\ -C\left(\omega \right)\\ S\left(\omega \right)\\ -S\left(\omega \right)\end{array}\right)}_{4n_s\times (n+1)},$$q_s^{\left(k\right)}=\frac{{\delta }_s}{2}{\left(\begin{array}{c}1\\ ...\\ 1\end{array}\right)}_{4n_s\times 1}$

$C\left(\omega \right)={\left(\begin{array}{c}{c}^T\left({\omega }_1^{\left(s\right)}\right)\\ ...\\ c^T\left({\omega }_{n_s}^{\left(s\right)}\right)\end{array}\right)}_{n_s\times (n+1)},$$S\left(\omega \right)={\left(\begin{array}{c}{s}^T\left({\omega }_1^{\left(s\right)}\right)\\ ...\\ s^T\left({\omega }_{n_s}^{\left(s\right)}\right)\end{array}\right)}_{n_s\times (n+1)}$

当k→∞时，H^(k)(ω)→H(ω)，当满足以上两式时，有|H(ω)|≤|c^T(ω)b_k|+ |s^T(ω)b_k|≤δ_s，即满足阻带幅值约束条件。

5）将约束最小二乘法设计转化为解一个正定二次规划问题，并将通带和阻带约束 统一表示。

基于约束最小二乘法设计滤波器的目标函数为

式中ω∈[0,π]={ω_i,i=1,...M}。

约束最小二乘法设计可转化为解一个正定二次规划问题：

$\underset{b_k}{\min }\frac{1}{2}{b}_k^T{H}_k{b}_k+b_k^T{p}_k$

式中H_k为(n+1)×(n+1)阶正定Hessian矩阵，其中：

通带和阻带约束可统一表示为如下形式

A_kb_k≤q_k

式中$A_k=\left(\begin{array}{c}{A}_p^{\left(k\right)}\\ A_s^{\left(k\right)}\end{array}\right),$$q_k=\left(\begin{array}{c}{q}_p^{\left(k\right)}\\ q_s^{\left(k\right)}\end{array}\right)$

6）利用MATLAB，根据二次规划的数学模型和约束条件进行直接求解。

采用约束最小二乘法设计的滤波器，其求解过程转化为求解一个正定二次规划问 题，二次规划的数学模型为

$\underset{x}{\min }\frac{1}{2}\mathrm{Hx}+f^{T}x$

其中约束条件为

Ax≤b

A_eqx=b_eq

lb≤x≤ub

式中，H为二次型矩阵，A、Aeq分别为不等式和等式约束的系数矩阵，其余为向量。

MATLAB提供了相关求解二次规划的函数，可以调用quadprog函数命令进行直接 求解。

实施例：

对500kV两圈变压器进行各种情况仿真分析，系统模型如图5所示，线路长度为 330km，变压器采用Yn/D-11接线，变比为500/13.8kV，转角方式为减零序的相电流差 动方式，采用星形侧向三角形侧转角。

IEC60044-8标准规定ECT的额定采样频率可取为20倍、48倍、80倍工频频率中 的任意一种，为使各条支路电流采样值采样频率一致，可以确定相应的抽取因子和插值 因子，其结构如图6所示。因此所设计的低通滤波器可能工作在两个频率下，即4kHz 和12kHz。

ECT1和ECT2为分别装设在变压器首端和末端的电子式电流互感器，为变压器差动 保护装置提供采样数据。设ECT1输出数据采样频率为1kHz，ECT2输出数据采样频率 为4kHz，采样频率转化算法在差动保护装置中进行，以高采样频率的采样数据向低采 样频率转化为例，设计FIR数字低通滤波器，其性能指标为：F_s=4kHz，f_p=450，f_s=550， 即通带和阻带的截止频率分别为0.225π和0.275π，阶数N=80，加权函数在通带和阻带 均取值为1，ω∈{kπ/1000；k=0,1,···,1000}，期望群延迟为25，通带纹波不大于1dB，阻带 衰减不小于25dB，迭代停止条件为tol=0.0001。

图7a-7d为该FIR数字低通滤波器的幅频响应、相频响应、群延迟和时域的单位脉 冲响应图，由图7(c)群延迟曲线可以看出，在通带内群延迟基本保持在25附近波动， 最大群延迟误差为2.716，峰值瞬时误差小于3.683%，其中基波幅值误差小于0.52%。 表1为所设计的FIR数字低通滤波器系数。

表1滤波器系数h(n)

由于以4kHz采样频率向1kHz采样频率转换，可确定插值因子U=1，抽取因子D=4， FIR滤波器采用上述基于约束最小二乘法设计的系数，并采用切比雪夫等波纹逼近法设 计的阶数N=80的低通滤波器，与其进行对比。采用切比雪夫等波纹逼近法设计的低通 滤波器，其群延迟在通带和阻带内均为常数，即延迟了40个采样点。分别对采样信号 进行采样频率转换，方法1为U=1倍插值、等波纹设计滤波算法和D=4倍抽取三者级 联的采样频率转换方法；方法2为U=1倍插值、基于约束最小二乘法设计的低延迟滤波 算法和D=4倍抽取三者级联的采样频率转换方法。

在t=0.04s时高压侧母线N区外F₁点发生三相短路故障，互感器ECT1、ECT2在 采样前经低通滤波器进行抗混叠滤波。以A相为例进行分析，图8给出了从ECT1和 ECT2的A相采样数据中恢复的连续波形I_1H，I_2L。

图9给出了高采样率的ECT2采样数据经过转角后，分别采用方法1和方法2后 采样电流输出值I_1L，可见方法1的转换方法输出滞后于方法2的输出，延迟约为 3.75ms。因此，采用等波纹设计滤波算法的普通滤波输出滞后于基于约束最小二乘法 设计的低延迟滤波算法，制约了保护动作的快速性，而低延迟滤波算法在完成滤波同 时不带来附加的输出延迟。

变压器保护方案为采样值差动保护判据，采用双折线采样值差动判据，转换后采 样频率为1kHz，因此选取数据窗长R=8，判别点数S=6，即连续8次判别中有6次及 以上则输出动作信号1，否则输出0，其动作信号如图10所示。方法1必须等待到正 确的滤波输出后方能进入差动保护算法，其动作时间明显大于方法2的动作时间。滤 波器阶数越高，其动作延迟时间也越大。

对转换后的电流采样数据进行误差分析，图11为采样频率转换后的计算的差电 流I_dL，I_dL=I_1L-I_2L，其最大瞬时值为1.029A，通过合理的整定，能够保证在区外故障 时差动保护的正确判别。

从图12中可以看出方法2的转换算法在小于500Hz频率段误差较小，很好地滤 除了I_1H中的高频分量，并且较完整地保留了I_1H中的低频分量，工频基波的幅值误差 相当于实际值的-60dB。计算方法2的有效值相对误差与DFT基波幅值相对误差，由 表2数据可见，由于频谱吻合良好，误差较低。

表2有效值相对误差与DFT基波幅值相对误差

t₀/ms 有效值误差(%) 基波幅值误差(%) 10 0.2184 0.1028 20 -0.1082 -0.0876 30 0.0323 0.0192 40 0.0795 0.0574 50 -0.0532 -0.0511