交流信号的干扰识别方法

Interference Identification Method of AC Signal

摘要: 微机保护中通常使用全周傅氏算法对离散的交流采样数据进行分析，但是交流信号中存在干扰数据时，会对相量分析产生影响，得到错误的系统量测结果。本文提出了一种仅利用原始交流采样序列生成若干个校验序列的方法，以此为基础提出了交流信号中干扰数据的识别方法，通过仿真计算的结果验证了判据的可行性和有效性。该方法还可以应用于电气量的测量中。

Abstract: Full-cycle Fourier algorithm is usually used to analyze discrete AC sampling data in microcomputer protection. When there is interference data in AC signal, it will affect phasor analysis and get wrong system measurement results. In this paper, a method of generating several check sequences using only the original AC sampling sequence is proposed, and based on this, a method of identifying the interference data in AC signal is proposed. The feasibility and validity of the criterion are verified by simulation results. This method can also be applied to the measurement of electrical volume.

1. 引言

目前在电力系统发生故障时，绝大多数计算机继电保护是用数学算法通过故障信号的基波相量来对故障进行分析和判断 [1] [2] [3] [4] [5]。在现在的电力系统继电保护应用中，通常采用全周傅氏算法对离散的交流采样数据进行计算处理，求得故障信号的基波有效值和相位，得到系统的量测数据，然后故障判断以实现各种继电保护的功能 [6] - [13]。但由于受到电磁扰动、传输畸变等影响，交流采样数据中可能会存在干扰数据 [14] [15]，此时全周傅式算法计算得到的基波有效值和相角与实际值之间存在较大误差，导致保护误动或拒动 [16]。此外，对其它测量的应用也会产生影响。

文献 [17] [18] [19] [20] [21] 指出，微机保护装置硬件方面抗干扰措施主要有电源滤波、屏蔽、隔离、接地等技术。然而装置结构、使用环境的不同，干扰程度也会不一样，硬件措施无论如何完善，也很难彻底消除干扰数据。如果干扰数据不幸被采集到，保护在软件上应进行数据的甄别。由于被干扰的是原始信号，难以通过量测结果直接甄别，而基于间断点识别、滤波分析的检测方案会引入较大的计算量和时延环节，不利于交流信号的实时分析 [22] [23]，为此本文在全周傅式算法原理的基础上，先针对单个干扰数据提出了一种仅利用原始交流采样序列生成若干个校验序列的方法，并提出了参考的识别门槛，通过仿真对其有效可行性进行判断。此方法还可以推广到若干个干扰数据的识别和判断。

2. 大数对全周傅氏算法的影响

对于采样得到的离散的电压、电流原始数据，微机保护通常采用全周傅氏算法计算工频量的有效值和相位，并通过相量计算得到保护判据。但是一次和二次系统均可能出现数据丢失、受电磁干扰等情况，从而产生很大异常数值的采样值(该异常采样值可简称为“大数”)，受大数的影响，通过傅氏算法无法得到正确的电量信息。大数的存在会对数据处理、整定计算产生影响，严重情况下甚至会导致保护装备产生拒动与误动。

假设此处的输入信号为： $x\left(t\right)=A\sin \left(\omega t+\phi \right)$，采样点数量N，大数 $x^{\prime }\left(p\right)=kx\left(p\right)$。当出现极端情况，AD满量程(正的最大或负的最大)输出时，可以用电流信号 $I^{\prime }=30\sqrt{2}I$ 近似替代，其中，I为额定电流；30为电流最大值的倍数。例如，当 $A=1$， $\omega =100\pi \left(\omega =2\pi f,f=50\text{Hz}\right)$，采样点数为30，干扰发生在第10个采样点时，有 $x^{\prime }\left(10\right)=30\sqrt{2}x\left(10\right)$，应用全周傅氏算法会发现此时计算值最大误差可达247%。

实际工程中，当线路发生故障时，保护系统需要快速准确地获取系统电量信息以构成保护判据。如果受到了严重的电磁扰动，相对误差将远远大于微机保护的误差允许值，对微机保护造成误动或拒动的影响。由于受干扰的是被测的原始信号，难以对受大数影响的计算结果进行直接判断。而基于间断点识别和滤波的检测方式会引入较大的计算量以及时延环节，不利于保护的快速动作。

本文主要利用截断误差的原理设定判据的门槛值，介绍一种利用电气量固有的特征实现对单个大数快速有效的识别。本文以单个干扰数据为例进行推导说明，此方法还可以推广到多个干扰数据上。

3. 单个大数识别算法

3.1. 基本思路

为清晰展示算法的思路，先用单一的工频信号进行阐述，假设输入信号 $x\left(t\right)=A\sin \left(\omega t+\phi \right)$，每周采样点数量为N，采样间隔为T_S，那么微机保护可以得到 $N+1$ 个采样值的原始采样数据集 $Q=\left\{x\left(0\right),x\left(1\right),x\left(2\right),\cdots ,x\left(N\right)\right\}$。在此条件下，关键点是：将原始采样数据集Q按照数据序列的奇偶顺序再得到两个附加的校验数据集 $R=\left\{x\left(1\right),x\left(3\right),\cdots ,x\left(N-1\right)\right\}$ 和 $S=\left\{x\left(0\right),x\left(2\right),x\left(4\right),\cdots ,x\left(N\right)\right\}$，简称为奇序列R，偶序列S，其中，奇、偶序列相当于采样间隔为2T_S。三个序列的示意图如图1所示。

Figure 1. Sampling data set, calibration data set schematic diagram

图1. 原始采样数据集、校验数据集示意图

利用上述三个采样序列集，可以实现：1) 正常无大数影响时，Q、R、S三个序列的计算结果误差较小且基本一致；2) 当采样值存在一个大数影响时，在Q、R、S三个序列中包含异常大数的计算结果误差较大，序列中不包含异常大数的计算结果误差较小，特征是三者的结果不一致。因此，利用(1)、(2)的特征差异，就可以实现识别单个干扰大数据。

在MATLAB中验证这一想法：以输入信号 $x\left(t\right)=3\sin \left(100\pi t+0.1\right)$ 为例，设采样点数量为 $N=30$，构建三个序列，运算结果如表1所示。

Table 1. Results of full-cycle Fourier algorithm for primitive, odd and even sequences

表1. 原、奇、偶序列的全周傅氏算法运算结果

3.2. 判据的门槛值

3.2.1. 傅氏算法截断误差分析

在全周傅氏算法中，通过使用梯形积分公式近似处理傅里叶级数中的系数积分项。梯形积分公式原理是将被积区间分为若干个小区间，利用梯形来近似计算曲线下方面积，实现数值积分的计算。利用梯形积分公式近似计算积分时将产生截断误差，截断误差受信号类型以及采样点数量的影响。数值分析知识指出，对于 $f\left(x\right)$，利用梯形积分计算其在区间 $\left[a,b\right]$ 上的定积分时所产生的截断误差满足式(1)。

$\left|E_{Tra}\right|=\left|{\displaystyle {\int }_a^{b}f\left(x\right)-Tra\left(N\right)}\right|\le \frac{b-a}{12}{\left(\frac{b-a}{N}\right)}^{2}M$ (1)

式中： $E_{Tra}$ 代表截断误差； $Tra\left(N\right)$ 代表采样点数量为N时梯形公式计算结果；M代表 $\left|f^{″}\left(x\right)\right|$ 在区间 $\left[a,b\right]$ 的最大值。

傅氏算法的基本思路来源于傅里叶级数，根据傅里叶级数的原理，基波的正、余弦项系数 $a_1$ 、 $b_1$ 可以表示为：

$a_1=\frac{2}{T}{\displaystyle {\int }_0^{T}x\left(t\right)\cos \left(\omega t\right)\text{d}t}$ (2)

$b_1=\frac{2}{T}{\displaystyle {\int }_0^{T}x\left(t\right)\sin \left(\omega t\right)\text{d}t}$ (3)

根据 $a_1$ 、 $b_1$ 可以计算出基波的有效值 $X_1$ 为：

$X_1=\sqrt{\frac{a_1^2+b_1^2}{2}}$ (4)

利用微机进行处理时，式(2)、(3)可以用积分梯形公式求得：

$a_{1N}=\frac{1}{N}\left[2{\displaystyle {\sum }_{k=1}^N{x}_k\sin \left(k\frac{2\pi }{N}\right)}\right]$ (5)

$b_{1N}=\frac{1}{N}\left[x_0+2{\displaystyle {\sum }_{k=1}^{N-1}{x}_k\cos \left(k\frac{2\pi }{N}\right)+x_N}\right]$ (6)

利用微机进行处理时将引入截断误差，为进一步分析截断误差对全周傅氏算法的影响，设输入信号 $x\left(t\right)=A\sin \left(\omega t+\phi \right)$，则可得到采样信号：

$x_k=A\sin \left(\omega t_k+\phi \right)=A\sin \left(\frac{2\pi }{T}\frac{k}{N}T+\phi \right)=A\sin \left(\frac{2\pi }{N}k+\phi \right)$ (7)

将 $x_k$ 代入式(5)，并对其进行化简，可以得到：

$\begin{array}{c}{a}_{1N}=\frac{1}{N}\left[2{\displaystyle {\sum }_{k=1}^N{x}_k\sin \left(k\frac{2\pi }{N}\right)}\right]\\ =\frac{A}{N}\left[{\displaystyle {\sum }_{k=1}^{N}2\sin \left(k\frac{2\pi }{N}+\phi \right)\sin \left(k\frac{2\pi }{N}\right)}\right]\\ =\frac{A}{N}\left[{\displaystyle {\sum }_{k=1}^N\cos \left(\phi \right)-\cos \left(k\frac{4\pi }{N}+\phi \right)}\right]\\ =A\cos \left(\phi \right)-\frac{A}{2\pi }{\displaystyle {\sum }_{k=1}^N\frac{2\pi }{N}}\cos \left(2\frac{2\pi }{N}k+\phi \right)\end{array}$ (8)

为了对式(8)的取值范围进行分析，构造 $g\left(\theta \right)=\cos \left(2\theta +\phi \right)$， $g\left(\theta \right)$ 在区间 $\left[0,2\pi \right]$ 上积分可表示为：

${\displaystyle {\int }_0^{2\pi }\cos \left(2\theta +\phi \right)\text{d}\theta }=\underset{N\to \infty }{\lim }{\displaystyle {\sum }_{k=1}^N\frac{2\pi }{N}\cos \left(2\frac{2\pi }{N}k+\phi \right)}=0$ (9)

利用梯形公式计算时，引入的截断误差：

$\begin{array}{c}\left|E_g\right|=\left|{\displaystyle {\int }_0^{2\pi }\cos \left(2\theta +\phi \right)\text{d}\theta }-{\displaystyle {\sum }_{k=1}^N\frac{2\pi }{N}}\cos \left(2\frac{2\pi }{N}k+\phi \right)\right|\\ =\left|-{\displaystyle {\sum }_{k=1}^N\frac{2\pi }{N}}\cos \left(2\frac{2\pi }{N}k+\phi \right)\right|\\ \le \frac{2\pi }{12}{\left(\frac{2\pi }{N}\right)}^{2}4=\frac{2\pi }{3}{\left(\frac{2\pi }{N}\right)}^2\end{array}$ (10)

将式(10)代入式(8)中化简，可得到基波的正弦项系数 $a_1$ 的取值范围：

$A\cos \left(\phi \right)-\frac{A}{3}{\left(\frac{2\pi }{N}\right)}^2\le a_1\le A\cos \left(\phi \right)+\frac{A}{3}{\left(\frac{2\pi }{N}\right)}^2$ (11)

同理，对于余弦项系数 $b_1$ 的取值范围有：

$A\sin \left(\phi \right)-\frac{A}{3}{\left(\frac{2\pi }{N}\right)}^2\le b_1\le A\sin \left(\phi \right)+\frac{A}{3}{\left(\frac{2\pi }{N}\right)}^2$ (12)

本节分析截断误差对全周傅氏算法影响的目的主要在于得出识别算法判据的形式为后续的数值拟合做准备，而不是得出具体的公式。经验证，无论初相角 $\phi$ 多大多小， $a_1$ 、 $b_1$ 为正或为负，都不会影响下文得出的判据形式。因此，为了简便，仅以 $a_1$ 、 $b_1$ 都为正来分析。

利用式(4)可以计算出基波有效值范围满足：

$\{\begin{array}{l}{X}_1=\sqrt{\frac{a_1^2+b_1^2}{2}}\le \frac{1}{\sqrt{2}}\sqrt{{\left(A\cos \left(\phi \right)+\frac{A}{3}{\left(\frac{2\pi }{N}\right)}^2\right)}^2+{\left(A\sin \left(\phi \right)+\frac{A}{3}{\left(\frac{2\pi }{N}\right)}^2\right)}^2}\\ X_1=\sqrt{\frac{a_1^2+b_1^2}{2}}\ge \frac{1}{\sqrt{2}}\sqrt{{\left(A\cos \left(\phi \right)-\frac{A}{3}{\left(\frac{2\pi }{N}\right)}^2\right)}^2+{\left(A\sin \left(\phi \right)-\frac{A}{3}{\left(\frac{2\pi }{N}\right)}^2\right)}^2}\end{array}$ (13)

对式(13)进行整理，此时 $a_1$ 、 $b_1$ 都为正，因此必定可得 $\sin \left(\phi +\frac{\pi }{4}\right)>0$，另外，为了下文拟合表达式更加简洁实用略去 ${\left(\frac{A}{3}{\left(\frac{2\pi }{N}\right)}^2\right)}^2$ 项(误差影响不会超过2%，且在后续的验证中可知影响不大)。当信号采样点数量为N，信号初相角为 $\phi$ 时，考虑截断误差影响时基波有效值满足：

$\frac{A}{\sqrt{2}}{\left(1-\frac{2\sqrt{2}}{3}{\left(\frac{2\pi }{N}\right)}^2\sin \left(\phi +\frac{\pi }{4}\right)\right)}^{\frac{1}{2}}\le X_1\le \frac{A}{\sqrt{2}}{\left(1+\frac{2\sqrt{2}}{3}{\left(\frac{2\pi }{N}\right)}^2\sin \left(\phi +\frac{\pi }{4}\right)\right)}^{\frac{1}{2}}$ (14)

式(14)中 $X_1$ 的取值范围为后续计算门槛值中的相对误差提供理论基础。

3.2.2. 校验集相对误差分析

为消除未知待求的幅值影响，以下分析采样数据集Q与校验数据集R的结果关系，数据集S类同。

将采样数据集Q与校验数据集R代入全周傅氏算法中，计算得到信号的有效值与相角 $\left(X_Q,{\phi }_Q\right)$， $\left(X_R,{\phi }_R\right)$。根据3.1中式(14)，可知其基波有效值满足：

$\frac{A}{\sqrt{2}}{\left(1-\frac{2\sqrt{2}}{3}{\left(\frac{2\pi }{N}\right)}^2\sin \left(\phi +\frac{\pi }{4}\right)\right)}^{\frac{1}{2}}\le X_Q\le \frac{A}{\sqrt{2}}{\left(1+\frac{2\sqrt{2}}{3}{\left(\frac{2\pi }{N}\right)}^2\sin \left(\phi +\frac{\pi }{4}\right)\right)}^{\frac{1}{2}}$ (15)

$\frac{A}{\sqrt{2}}{\left(1-\frac{2\sqrt{2}}{3}{\left(\frac{4\pi }{N}\right)}^2\sin \left(\phi +\frac{\pi }{4}\right)\right)}^{\frac{1}{2}}\le X_R\le \frac{A}{\sqrt{2}}{\left(1+\frac{2\sqrt{2}}{3}{\left(\frac{4\pi }{N}\right)}^2\sin \left(\phi +\frac{\pi }{4}\right)\right)}^{\frac{1}{2}}$ (16)

因此，校验数据集与采样数据集计算结果之间的相对误差 ${\delta }_R$ 满足：

$\frac{X_{R.\min }-X_{Q.\max }}{X_{Q.\max }}\le {\delta }_R=\frac{X_R-X_Q}{X_Q}\le \frac{X_{R.\max }-X_{Q.\min }}{X_{Q.\min }}$ (17)

将式(15)、(16)代入式(17)，可得 ${\delta }_R$ 范围：

$\{\begin{array}{l}{\delta }_R\le {\delta }_{R.\max }=\frac{\sqrt{1+4h\left(N,b\right)}-\sqrt{1-h\left(N,b\right)}}{\sqrt{1-h\left(N,b\right)}}\\ {\delta }_R\ge {\delta }_{R.\min }=\frac{\sqrt{1-4h\left(N,b\right)}-\sqrt{1+h\left(N,b\right)}}{\sqrt{1+h\left(N,b\right)}}\end{array}$ (18)

式中， $h\left(N,b\right)=\frac{2\sqrt{2}}{3}{\left(\frac{2\pi }{N}\right)}^2\sin \left(\phi +\frac{\pi }{4}\right)$。

虽然傅氏算法的截断误差受输入信号的幅值影响，但利用采样数据集与校验数据集来计算相对误差后消除了信号幅值的影响，同时得到了输入信号不受大数影响情况下相对误差 ${\delta }_R$ 的取值区间。将相对误差 ${\delta }_R$ 的区间边界做等价无穷小代换并略去高阶无穷小后，可以近似认为其区间边界满足：

$\{\begin{array}{l}{\delta }_{R.\max }\left(N,\phi \right)=\frac{c_1\sin \left(c_0\phi +c_2\right)+c_3}{N}\\ {\delta }_{R.\min }\left(N,\phi \right)=\frac{c_4\sin \left(c_0\phi +c_5\right)+c_6}{N}\end{array}$ (19)

虽然傅氏算法的截断误差受输入信号的幅值影响，但利用采样数据集与校验数据集来计算相对误差后消除了信号幅值的影响，同时得到了输入信号不受大数影响情况下相对误差 ${\delta }_R$ 的取值区间。将相对误差 ${\delta }_R$ 的区间边界做等价无穷小代换并略去高阶无穷小后，可以近似认为其区间边界满足：

式中， $\phi$ 表示采样正弦信号的初相角；N表示采样点数； $c_0~c_6$ 为拟合系数(下文推导具体数值，且分析表明与 $\phi$ 、N无关)。

由此可以得出结论，当采样数据中不存在大数时，利用全周傅氏变换计算校验数据集R(或S)的有效值 $X_R$ (或 $X_S$ )与采样数据集Q的有效值 $X_Q$， $X_R$ (或 $X_S$ )与 $X_Q$ 之间的相对误差 $\delta$ 落在 $\left[{\delta }_{R.\min },{\delta }_{R.\max }\right]$ (或 $\left[{\delta }_{S.\min },{\delta }_{S.\max }\right]$ )区间，反之则不落在区间内。

4. 验证

利用实际数据搭建计算模型，改变输入信号的初相角 $\phi$ 与采样点数N计算相对误差，绘制出 ${\delta }_R$ (由于 ${\delta }_R\cdot N$ 是确定数值，同时为了降低图像维度使结果清晰明了)的分布图如图2所示。

Figure 2. ${\delta }_R\cdot N$ distribution without single large number

图2. 不存在大数时 ${\delta }_R\cdot N$ 分布图

由 ${\delta }_R$ 分布的上下限根据最小二乘法可以求得 $c_0~c_6$ 的具体数值，得到区间边界的表达式如式(20)所示。经校验，该表达式下最小二乘拟合的 $R^2\ge 0.995$。改变数据进行重复仿真校验，利用最小二乘法求得的上下限表达式始终为相对误差 ${\delta }_R$ 分布的包络线。

$\{\begin{array}{l}{\delta }_{R.\max }\left(N,\phi \right)=\frac{0.4975\sin \left(2\phi +1.412\right)-0.027}{N}\\ {\delta }_{R.\min }\left(N,\phi \right)=\frac{0.4798\sin \left(2\phi +1.325\right)-0.081}{N}\end{array}$ (20)

同理可以求出，校验数据集S与采样数据集Q间相对误差 ${\delta }_S$ 的取值区间 $\left[{\delta }_{S.\min },{\delta }_{S.\max }\right]$ 满足：

$\{\begin{array}{l}{\delta }_{S.\max }\left(N,\phi \right)=\frac{0.4975\sin \left(2\phi +\frac{2\pi }{N}+1.412\right)-0.027}{N}\\ {\delta }_{S.\min }\left(N,\phi \right)=\frac{0.4798\sin \left(2\phi +\frac{2\pi }{N}+1.325\right)-0.081}{N}\end{array}$ (21)

考虑到短路电流中包含非周期分量与谐波分量的影响，在实际应用时可以将判别阈值适度放大，为 $\left[{\delta }_{\min }\left(1-\Delta d\right),{\delta }_{\max }\left(1+\Delta d\right)\right]$， $\Delta d$ 为区间裕度，此处取 $\Delta d=5\%$，限于篇幅本文不再赘述。

在这个基础上利用实际数据再次搭建计算模型，在输入信号中混入了单个大数 $x^{\prime }\left(3\right)=1.2x\left(3\right)$，改变输入信号的初相角 $\phi$ 与采样点数N计算相对误差 $\delta$，绘制出 $\delta \cdot N$ (N是确定数值，同时降低图像维度)的分布图以及 $\delta \cdot N$ 取值区间 $\left[{\delta }_{\min }\cdot N,{\delta }_{\max }\cdot N\right]$ 如图3所示。

Figure 3. ${\delta }_R\cdot N$ distribution diagram with single large number ( $x^{\prime }\left(3\right)=1.2x(\; 3\; ))$

图3. 存在单个大数( $x^{\prime }\left(3\right)=1.2x\left(3\right)$ )时 ${\delta }_R\cdot N$ 分布图

由图3可以看出，除了在系数过零点附近(过零点处的大数影响可以忽略不计)，利用本文提出的方法均可实现单个大数的识别。

5. 结语

针对交流采样数据中可能会存在干扰数据导致保护误动或拒动的问题，本文提出了一种通过利用原始交流采样序列生成若干个校验序列来识别的方法，并提出了参考的识别判据门槛。

本文以单个干扰数据为例，将送入微机的采样间隔为T_S的采样数据集按序列的奇偶性分为两个采样间隔为2T_S的校验数据集，分别利用全周傅氏算法计算采样结果，再通过对全周傅氏算法的截断误差进行分析，得到了单个大数存在的判据门槛值，最后分析采样结果实现对单个大数的识别。针对非单个干扰数据，可以采用如图1类似的方法再构建出采样间隔为3T_S (甚至4T_S)的校验序列。

仿真计算结果验证了该方法的可行性和有效性。该方法避免了使用间断点识别和高频滤波算法，简化计算的同时减少了时延环节，具有形式简单、精准度高的特点，可以实现“实时计算、实时决策”。