光电互感器中多数字源采样率转换的快速滤波方法

摘要

本发明公开了一种光电互感器中多数字源采样率转换的快速滤波方法，通过将内插和抽取进行级联的方式实现任意分数倍采样频率的转换，以巴特沃斯模拟滤波器为基础，根据不同数字源采样频率计算出模拟滤波器截止频率，以实现抗混叠和抗镜像效果；然后根据滤波器通带最大平坦和过渡带最窄的原则，并采样双线性变换，得到z域低通滤波器方程参数；最后采用FRR方法进行滤波和反相设计，最终得到最大平坦和零相位失真的输出信号序列。本发明提供的方法，在光电互感器差动保护中，可实现较低阶数的快速滤波，以得到最大的输出平坦和零相位；简单实用，有利于保护的快速动作，能够减少频率混叠造成的误差，实现了采样频率的归一化。

说明书

技术领域

本发明涉及一种光电互感器中多数字源采样率转换的快速滤波方法，基于零相位和最大平坦算法对信号进行处理，属于无限冲激响应（IIR）滤波器技术。

背景技术

光电互感器（OCT）是电子式电流互感器（Electronic Current Transformer，ECT）的一种。近年来，电子式电流互感器的研究已经取得了巨大的进展，同时随着数字化变电站的推广，光电互感器的应用更是进入到了实用化阶段。ECT/OCT往往会集成在气体绝缘组合开关（GIS）内部或安装在断路器和变压器套管上，由各自的一次设备制造商供货，因此，当OCT由不同的一次设备供应商供货时，可能会出现各个OCT的采样数据输出频率不一致的情况。而基于采样值差动保护原理的保护装置通常需要获得瞬时差电流值，要求与被保护一次设备相连的所有支路电流采样值的采样频率必须一致，因此当不同数字源的采样频率不一致时，必须通过重采样将采样频率转换成一致，即实现采样频率的归一化，以确保保护装置的正确判别故障区域和动作。

传统的低通滤波器采用有限冲激响应（FIR）滤波器。FIR滤波器获得零响应需要附加的延时，延时为(N-1)/2个采样间隔，其中N为滤波器长度，采样间隔时间为1/(Uf_m)（U为插值量f_m为转换前采样频率）。滤波器长度取得越长，越接近理想幅频特性，从而输出精度越高，但同时输出延时也越长，从而大大延长了保护的动作时间，不利于保护的快速动作。在选取滤波器长度时，应在满足保护精度要求的前提下尽可能减小滤波长度，从而降低对保护动作速度的影响。

FIR滤波器的常用设计方法有窗函数法，频率采样法及等纹波逼近法。窗函数法、频率采样法设计简单，容易应用，但它们对设计指标只能进行有限控制；基于切比雪夫原理的等纹波设计方法可对理想滤波器做最佳逼近，因而可获得较好的通带和阻带特性；但是切比雪夫通带的纹波特性，影响了输出的幅值特性。

现代的多数字源采样处理过程中，通常采用三次样条插值对信号进行高频转低频处理；由于三次样条插值中混合的较高频率数据，造成频谱混乱，产生了较大的差动电流，引起了装置的误动。

发明内容

发明目的：为了克服现有技术中存在的不足，本发明提供一种光电互感器中多数字源采样率转换的快速滤波方法，是一种零相位最大平坦滤波方法，该方法基于巴特沃斯模拟滤波方法，在最大平坦性上采用FRR（forward filter-reverse filter-reverse output）方法，以得到精确零相位失真的输出序列，提高保护装置快速准确动作。

技术方案：为实现上述目的，本发明采用的技术方案为：

光电互感器中多数字源采样率转换的快速滤波方法，通过将内插和抽取进行级联的方式实现任意分数倍采样频率的转换，将内插中的抗镜像滤波器和抽取中的抗混叠滤波器合并为一个低通滤波器；然后根据滤波器通带最大平坦和过渡带最窄的原则，求解得到低通滤波器方程系数向量；最后采用FRR方法进行滤波和反相设计，最终得到最大平坦和零相位失真的输出信号序列。具体设计时，以巴特沃斯滤波器为原型，采用双线性变换法，设计出最大平坦无限脉冲响应数字滤波特性；然后根据FRR方法，先将输出序列按顺序滤波，然后将所得结果逆转后反向通过滤波器，再将所得结果逆转后输出，以得到精确的零相位失真的输出序列；该方法具体包括如下步骤：

（1）插值采样比U/D计算：根据归一化原理，对两个数字源光电互感器的信号采样值进行频率归一化，根据两数字源采样频率比值，得到插值U和采样D；

（2）截止频率设计：通过将内插和抽取进行级联的方式实现任意分数倍采样频率的转换，将内插中的抗镜像滤波和抽取中的抗混叠滤波合并在一个IIR数字低通滤波器内；根据计算得到的插值U和采样D，设计IIR数字低通滤波器的截止频率，所述截止频率同时满足抗镜像滤波和抗混叠滤波的频率要求；

（3）幅频响应通带的最大平坦性设计：根据低通滤波器通带最大平坦和过渡带最窄的原则，先设计一个等效模拟滤波器，然后根据双线性变换将所述等效模拟滤波器映射为IIR数字低通滤波器，设计参数包括截止频率和放大倍数；具体方法为：

设计的等效模拟滤波器为一个兼具抗镜像滤波和抗混叠滤波的n阶巴特沃斯（Butterworth）低通滤波器，所述n阶巴特沃斯低通滤波器的传递函数H(s)为：

$>H\left(s\right)=\frac{a}{b_n{s}^n+b_{n-1}{s}^{n-1}+···+b_{1}s+b_0}$>

双线性变换公式为：

$>S=\frac{2}{T_S}\frac{(1-Z^{-1})}{(1+Z^{-1})}$>

其中，T_S为周期系数；

将双线性变换公式带入H(s)得到Z域的5阶巴特沃斯滤波器传递函数H(z)为：

$>H\left(z\right)=\frac{a_0+a_1{z}^{-1}+a_2{z}^{-2}+a_3{z}^{-3}+a_4{z}^{-4}+a_5{z}^{-5}}{b_0+b_1{z}^{-1}+b_2{z}^{-2}+b_3{z}^{-3}+b_4{z}^{-4}+b_5{z}^{-5}}$>

所述Z域的5阶巴特沃斯滤波器传递函数即作为IIR数字低通滤波器的传递函数；

（4）零相位设计：根据FRR（正通反向滤波再反向）方法，使数字源光电互感器的信号先顺序通过IIR数字低通滤波器，然后反相再通过IIR数字低通滤波器，最后将输出反向得到滤波结果，即精确零相位失真序列。

本发明中，等效模拟滤波器以巴特沃斯滤波器低通滤波器为原型，设计出了最大平坦无限脉冲响应数字滤波特性，实现了较低阶数较好的选频；在由等效模拟滤波器向数字低通滤波器的转换过程中，采用了双线性变换，该方法不会产生频率特性的混叠失真，能够克服模拟频率和非线性关系造成的幅频特性失真。

下面就本发明的设计过程进行说明。

为了避免在抽取过程中的频谱混叠和差值过程中的高频镜像，在实现抽取和差值的过程中必须设置IIR数字低通滤波器对输入信号进行滤波处理；首先对频率为f_s的输入信号x(n)进行差值，每相邻两个点之间插入(U-1)个零，然后使插值后的信号通过抗镜像滤波器，即可得到信号u(n)；对信号u(n)进行抗混叠滤波处理，滤除信号中的高频分量，再进行抽取处理，每隔D个点抽取一个点，得到输出信号y(n)，输出信号y(n)的频率为f_y=Uf_s/D。

设置h₁(n)为抗混叠滤波器，h₂(n)为抗镜像滤波器，其中h₁(n)和h₂(n)的输入信号频率均为f=Uf_s，故可以将抗混叠滤波器和抗镜像滤波器等效设计为一个IIR数字低通滤波器h(n)，所述IIR数字低通滤波器h(n)的理想特性为：

其中，ω为信号角频率。

为了实现差电流幅值输出的准确性，需要对IIR数字低通滤波器进行最大平坦性设计；考虑到采样值差动保护原理对采样点幅值要求较高的特性，应尽可能减少在滤波过程中带来的幅值误差，所以采用最大平坦通带以减小在频率转换过程中造成的幅值误差；而模拟滤波器中的巴特沃斯滤波器具有最大平坦幅度特性，因此以巴特沃斯低通滤波器为原型，采用双线性变换法可以设计出最大平坦IIR数字低通滤波器H(z)；5阶的最大平坦IIR数字低通滤波器系统函数为：

$>H\left(z\right)=\frac{a_0+a_1{z}^{-1}+a_2{z}^{-2}+a_3{z}^{-3}+a_4{z}^{-4}+a_5{z}^{-5}}{b_0+b_1{z}^{-1}+b_2{z}^{-2}+b_3{z}^{-3}+b_4{z}^{-4}+b_5{z}^{-5}}$>

完成最大平坦设计后，为了保证不失真，需要进行零相位数字滤波设计；零相位数字滤波的实现方法可以采用FRR方法，其主要过程为：先将输入序列按顺序滤波（Forward Filter），然后将所得结果逆转（last in first out，LIFO）后反向通过滤波器（Reverse Filter），再将所得结果逆转后输出（Reverse Output），即得精确零相位失真的输出序列；其原理分析如下：

y₁(n)=x(n)*h(n)

y₂(n)=y₁(M-1-n)

y₃(n)=y₂(n)*h(n)

y(n)=y₃(M-1-n)

其中，x(n)为输入信号序列；M为输入信号序列x(n)的信号序列长度；h(n)为IIR数字低通滤波器的冲激响应序列；y(n)为第二次滤波后的逆转序列即FRR的输出序列。

FRR响应的频域表示为：

Y₁(e^jω)=X(e^jω)H(e^jω)

Y₂(e^jω)=e^-jω(N-1)Y₁(e^-jω)

Y₃(e^jω)=Y₂(e^jω)H(e^jω)

Y(e^jω)=e^-jω(N-1)Y₃(e^jω)

可得：Y(e^jω)=X(e^jω)|H(e^jω)|²。

输出序列Y(e^jω)和输入序列X(e^jω)之间不存在相移，零相位滤波在理论上可以实现精确的零相位滤波；通过对过去时刻的数据x(n_P-M+1)与当前时刻的数据x(n_P)之间的M个数据进行滤波处理，可得到信号数据x(n_P-M+1)的滤波修正值；其中n_P为x(n)序列中任意值，n_P≥M-1。

在充分考虑M对延时与稳定性的影响后，选取M=7，采样频率为4kHz时，缓存延时为：t_d=M/f_s=1.75ms，1.75ms的缓存延时不足1/4个工频周期（5ms），对于保护的动作时间来讲可以满足其响应要求。

有益效果：本发明提供的光电互感器中多数字源采样率转换的快速滤波方法，在光电互感器差动保护中，可以实现较低阶数的快速滤波，以得到最大的输出平坦和零相位；本发明相对简单实用，有利于保护的快速动作，能够减少频率混叠造成的误差，实现了采样频率的归一化，一定程度上降低了输出延迟，实现了快速滤波，有利于保护的快速动作。

附图说明

图1为本发明的流程框图；

图2为信号采样率转换示意图；

图3为N=3/4/5时巴特沃斯滤波器幅频响应示意图；

图4为N=5时巴特沃斯滤波器阶跃响应示意图；

图5为PSRC EMTP模型图；

图6(a)为OCT1采样电流I_1H图；

图6(b)为OCT2采样电流I_2L图；

图7(a)为I_1H序列频谱图；

图7(b)为I_2L序列频谱图；

图8(a)为OCT转换后的电流I_1L图；

图8(b)为差电流I_dL图；

图9(a)为OCT1转换后的电流I′_1L图；

图9(b)为差电流I′_dL图；

图10(a)为OCT2提高采样率后电流I_2H图；

图10(b)为差电流I_dH图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

光电互感器中多数字源采样率转换的快速滤波方法，首先以巴特沃斯滤波器为原型，采用双线性变换法，设计出最大平坦无限脉冲响应数字滤波特性；然后根据FRR方法，先将输出序列按顺序滤波，然后将所得结果逆转后反向通过滤波器，再将所得结果逆转后输出，以得到精确的零相位失真的输出序列。

结合图1，就本发明的实现步骤给予说明。

（1）插值采样比U/D计算：根据归一化原理，对两个数字源光电互感器的信号采样值进行频率归一化，根据两数字源采样频率比值，得到插值U和采样D。

（2）截止频率设计：通过将内插和抽取进行级联的方式实现任意分数倍采样频率的转换，将内插中的抗镜像滤波和抽取中的抗混叠滤波合并在一个IIR数字低通滤波器内；根据计算得到的插值U和采样D，设计IIR数字低通滤波器的截止频率，所述截止频率同时满足抗镜像滤波和抗混叠滤波的频率要求。

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

组合后的IIR数字低通滤波器的理想低通频率响应特性为：

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

（3）幅频响应通带的最大平坦性设计：根据低通滤波器通带最大平坦和过渡带最窄的原则，先设计一个等效模拟滤波器，然后根据双线性变换将所述等效模拟滤波器映射为IIR数字滤低通波器，设计参数包括截止频率和放大倍数。

考虑到采样值差动保护原理对采样点幅值要求较高的特性，应尽可能减少在滤波过程中带来的幅值误差，所以采用最大平坦型通带以减小在频率转换过程中造成的幅值误差。模拟滤波器中巴特沃斯滤波器具有具有最大平坦幅度特性，因此借助巴特沃斯模拟滤波器为原型，采用双线性变换法可以设计出最大平坦IIR数字低通滤波器H(z)。

设计的等效模拟滤波器为一个兼具抗镜像滤波和抗混叠滤波的n阶巴特沃斯低通滤波器，所述n阶巴特沃斯低通滤波器的传递函数H(s)为：

$>H\left(s\right)=\frac{a}{b_n{s}^n+b_{n-1}{s}^{n-1}+···+b_{1}s+b_0}$>

双线性变换公式为：

$>S=\frac{2}{T_S}\frac{(1-Z^{-1})}{(1+Z^{-1})}$>

其中，T_S为周期系数。

从图3以N表示Z域的5阶巴特沃斯滤波器的阶数，当N=3/4/5时的幅频响应曲线可以看出，当N越大时，滤波器的过渡带越窄，滤波效果越好，但滤波时的计算量和所需存储的数据量也越大；当N=5时，过渡带带宽为600Hz，可以满足设计需要；本案的N=5的最大平坦IIR数字低通滤波器的系统函数为：

$>\left(\begin{array}{c}H\left(z\right)=\frac{0.1084+0.5419z^{-1}+{1.0837z}^{-2}+{1.0837z}^{-3}+{0.5419z}^{-4}+{0.1084z}^{-5}}{1+{0.9854z}^{-1}+{0.9739z}^{-2}+{0.3864z}^{-3}+{0.1112z}^{-4}+{0.0113z}^{-5}}\\ =0.10837·\frac{1+z^{-1}}{1+{0.15838z}^{-1}}·\frac{1+{2z}^{-1}+z^{-2}}{1+{0.34929z}^{-1}+{0.13031z}^{-2}}·\frac{1+{2z}^{-1}+z^{-2}}{1+0.47655z^{-1}+{0.54572z}^{-2}}\end{array}\right)$>

（4）零相位设计：根据FRR方法，使数字源光电互感器的信号先顺序通过IIR数字低通滤波器，然后再反相通过IIR数字低通滤波器，最后将输出反向得到滤波结果，即精确零相位失真序列。

零相位数字滤波的实现方法可以采用FRR方法，其主要过程为：先将输入序列按顺序滤波（Forward Filter），然后将所得结果逆转（last in first out，LIFO）后反向通过滤波器（Reverse Filter），再将所得结果逆转后输出（Reverse Output），即得精确零相位失真的输出序列；其原理分析如下：

y₁(n)=x(n)*h(n)

y₂(n)=y₁(M-1-n)

y₃(n)=y₂(n)*h(n)

y(n)=y₃(M-1-n)

其中，x(n)为输入信号序列；M为输入信号序列x(n)的信号序列长度；h(n)为IIR数字低通滤波器的冲激响应序列；y(n)为第二次滤波后的逆转序列即FRR的输出序列。

FRR响应的频域表示为：

Y₁(e^jω)=X(e^jω)H(e^jω)

Y₂(e^jω)=e^-jω(N-1)Y₁(e^-jω)

Y₃(e^jω)=Y₂(e^jω)H(e^jω)

Y(e^jω)=e^-jω(N-1)Y₃(e^jω)

可得：Y(e^jω)=X(e^jω)|H(e^jω)|²。

输出序列Y(e^jω)和输入序列X(e^jω)之间不存在相移，零相位滤波在理论上可以实现精确的零相位滤波；通过对过去时刻的数据x(n_P-M+1)与当前时刻的数据x(n_P)之间的M个数据进行滤波处理，可得到信号数据x(n_P-M+1)的滤波修正值；其中n_P为x(n)序列中任意值，n_P≥M-1。

该FRR方法在实现时需要缓存M个采用周期的数据，其缓存延时t_d由信号序列长度M与信号采用频率f_s决定，序列长度M越小，缓存延时越短。同时从图4观察到所设计的IIR数字低通滤波器的阶跃响应至少需要6个点才能够稳定，因此对于进入IIR数字低通滤波器的信号序列，其长度M越到，受到滤波器边缘效应的影响越小，滤波后输出的最终值越稳定。

在充分考虑M对延时与稳定性的影响后，选取M=7，采样频率为4kHz时，缓存延时为：t_d=M/f_s=1.75ms，1.75ms的缓存延时不足1/4个工频周期（5ms），对于保护的动作时间来讲可以满足其响应要求。

下面结合实施例对本发明做出进一步的说明。

仿真模型采用IEEE PES保护专委会PSRC（Power System Relaying Committee）推荐的EMTP参考电网模型，如图5所示。该模型是一个三机系统，包含四条输电线路，其中有一对是双回线，另外还包含了一个T型接线，其中所有输电线路均采用分布参数模型。

选取线路L4作为测试线路，将OCT1和OCT2装设在线路L4的首段和末端。根据IEC60044-8标准对ECT输出数据采样频率的规定，选取OCT1输出数据采样频率为4kHz，OCT2输出数据采样频率为2.4kHz。

在t=0s时母线4区外K点发生单相接地故障，OCT1、OCT2在采样前经本案设计的IIR数字低通滤波器进行抗混叠滤波，以A相为例进行分析。图6(a)和图6(b)为故障发生时t=0到t=20ms时OCT1和OCT2的采样数据I_1H、I_2L及其对应的采样图。由于OCT1采样频率高于OCT2，从频谱图，图7(a)和图7(b)可以清晰的看到，OCT1的频谱能量在1200～2000Hz正是OCT2所不含的。如果直接进行简单的抽取和插值操作，而不经过滤波处理则会造成较大的差电流。

因此通过数字算法将不同采样频率的OCT采样数据转化为同一采样频率的同时，需要进行滤波算法处理。滤波器的长度对算法的执行时间起着至关重要的作用，合适的选择滤波器的长度是平衡滤波效果和保护动作时间之间的重要环节。

采用降低频率实现频率归一化

根据MRDSP的基本理论，将采样率为4KHz的OCT1采样数据转化为2.4KHz，取内插因子U=3，抽取因子D=5，实现分数倍频率转换。本案所设计的IIR数字低通滤波器的阶数N=5，信号序列长度M=7，滤波器的输出延迟t_d约为1.8ms。

OCT1转换后的电流I_1L如图8(a)所示，图8(b)给出了差电流I_dL=I_1L-I_2L，其最大瞬时差电流值为66.68A，远小于差动保护的启动门槛。可见，线路区外故障时，采样频率转换算法计算误差不会影响保护的动作特性。为了进行对比分析，在这里给出将OCT1的采样数据I_1H进行三次样条函数插值得到采样频率为2.4KHz的I′_1L，图9(a)给出了OCT1转换后的电流I′_1L，图9(b)为其对应的差电流I′_dL，从图中可以看出最大差电流值达到343.86A，为瞬时最大峰值电流的11.95%。如果直接进入保护装置，则可能会引起保护的误动。根据香农采样定理，其中差电流形成的主要原因是在对OCT1的采样数据I_1H进行三次样条插值来降低采样频率过程中，导致OCT1采样数据中1200Hz～2000Hz的频谱映射到低频段，造成频谱混叠，这是形成较大差电流的主要原因。由此可见，在将高采样率数据转换为低采样率数据的同时必须滤除高采样率数据中的高频分量以避免频谱混叠。

采用提高采样频率实现频率归一化

为了改善对信号序列的内插，文中使用最简单且最近似的方法，就是将上采样操作中的所插入的那U-1个零值分别用当前的x(n)值进行替换。这样的功能能够通过一个数字采样-保持（SH元件）来实现，即对其输入端的脉冲周期T_i内所接收的输入值，以一个提高的频率f_o=Uf_i=U/T_i，在输出端共输出U次。具有这种功能的SH元件被表示为一个长度为N=U的抗镜像滤波器，其冲激响应为：

$>g\left(n\right)=1,\forall n=\mathrm{0,1},···,N=U-1$>

相应的振幅响应为：

$>G_o\left(\Omega \right)=\frac{\sin (\mathrm{N\Omega }/2)}{\sin (\Omega /2)}=\frac{\sin (\mathrm{U\Omega }/2)}{\sin (\Omega /2)}$>

其中，Ω为信号角频率

虽然该方法简单，但内插的质量很低，文中对其进行改善；将I=2个这样的元件级联，对于内插因子U=5的级联结构的冲激响应为：

$>g_s\left(n\right)=\frac{1}{5}·g\left(n\right)*g\left(n\right)$>

振幅相应为：

$>{G_o}^{\left(2\right)}\left(\Omega \right)=5{\left[\frac{\sin (5\Omega /2)}{5·\sin (\Omega /2)}\right]}^2$>

式中，级联系统的尺度与子系统的个数I无关，它总是等于内插因子U=5。文中采用I=2的级联系统，因为此时内插质量与单个的SH元件相比，具有明显的改善，而其必要的附加能耗较小。这种级联系统会在原信号的相邻采样值之间，进行具有线性特性的内插，采用此内插算法，结合多采样率转换方法，可以更加有效的实现频率转换，由图10(b)仿真的差电流可以看到其最大的差电流仅为I_dH=17.59A，输出延迟为2.417ms，效果较好。

通过以上仿真可以看到，采用降低采样频率和提高采样频率实现采样频率归一化，均需要滤除高采样率采样数据中的高频分量，避免频率混叠。仿真验证了本案所设计的滤波方法能够以较低的阶数获得较好滤波效果优良特性，同时实现了零相位滤波，本案所设计的滤波算法的输出延迟较小，有利于保护的快速动作。在提高采样率实现采样频率归一化的过程中，文中提出具有线性特性的内插算法，仿真发现其差点流值较小，但输出延迟有所增加。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。