CN110412415B - 一种基于dft、多阶滤波和突变判据的同步相量计算方法 - Google Patents

一种基于dft、多阶滤波和突变判据的同步相量计算方法 Download PDF

Info

Publication number
CN110412415B
CN110412415B CN201910517446.4A CN201910517446A CN110412415B CN 110412415 B CN110412415 B CN 110412415B CN 201910517446 A CN201910517446 A CN 201910517446A CN 110412415 B CN110412415 B CN 110412415B
Authority
CN
China
Prior art keywords
signal
frequency
calculation
phase
measured
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910517446.4A
Other languages
English (en)
Other versions
CN110412415A (zh
Inventor
蔡玉朋
李敏
张武洋
刘爱民
廖于翔
卢盛阳
楚天丰
刘娟
史松杰
吴蒙
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan Zhongyuan Huadian Science & Technology Co ltd
State Grid Corp of China SGCC
State Grid Liaoning Electric Power Co Ltd
Electric Power Research Institute of State Grid Liaoning Electric Power Co Ltd
Original Assignee
Wuhan Zhongyuan Huadian Science & Technology Co ltd
State Grid Corp of China SGCC
State Grid Liaoning Electric Power Co Ltd
Electric Power Research Institute of State Grid Liaoning Electric Power Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Wuhan Zhongyuan Huadian Science & Technology Co ltd, State Grid Corp of China SGCC, State Grid Liaoning Electric Power Co Ltd, Electric Power Research Institute of State Grid Liaoning Electric Power Co Ltd filed Critical Wuhan Zhongyuan Huadian Science & Technology Co ltd
Priority to CN201910517446.4A priority Critical patent/CN110412415B/zh
Publication of CN110412415A publication Critical patent/CN110412415A/zh
Application granted granted Critical
Publication of CN110412415B publication Critical patent/CN110412415B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R31/00—Arrangements for testing electric properties; Arrangements for locating electric faults; Arrangements for electrical testing characterised by what is being tested not provided for elsewhere
    • G01R31/08—Locating faults in cables, transmission lines, or networks
    • G01R31/081—Locating faults in cables, transmission lines, or networks according to type of conductors
    • G01R31/086—Locating faults in cables, transmission lines, or networks according to type of conductors in power transmission or distribution networks, i.e. with interconnected conductors
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R31/00—Arrangements for testing electric properties; Arrangements for locating electric faults; Arrangements for electrical testing characterised by what is being tested not provided for elsewhere
    • G01R31/08—Locating faults in cables, transmission lines, or networks
    • G01R31/088—Aspects of digital computing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E40/00—Technologies for an efficient electrical power generation, transmission or distribution
    • Y02E40/70—Smart grids as climate change mitigation technology in the energy generation sector
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E60/00—Enabling technologies; Technologies with a potential or indirect contribution to GHG emissions mitigation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y04—INFORMATION OR COMMUNICATION TECHNOLOGIES HAVING AN IMPACT ON OTHER TECHNOLOGY AREAS
    • Y04S—SYSTEMS INTEGRATING TECHNOLOGIES RELATED TO POWER NETWORK OPERATION, COMMUNICATION OR INFORMATION TECHNOLOGIES FOR IMPROVING THE ELECTRICAL POWER GENERATION, TRANSMISSION, DISTRIBUTION, MANAGEMENT OR USAGE, i.e. SMART GRIDS
    • Y04S10/00—Systems supporting electrical power generation, transmission or distribution
    • Y04S10/22—Flexible AC transmission systems [FACTS] or power factor or reactive power compensating or correcting units

Abstract

本发明涉及电力系统同步相量的计算方法技术领域，尤其涉及一种基于DFT、多阶滤波和突变判据的同步相量计算方法，具体是一种基于DFT、采用多阶数滤波、突变判据计算的同步相量计算方法。包括以下步骤：初始计算、频率计算、幅值与相位的计算、频率变化率的计算、信号突变判断算法以及算法测试。本发明通过突变判据的计算将幅值突变、相位突变与其它几种稳态或准稳态情况区别开来处理。通过多阶滤波器的设计，提升计算的精度，同时又具备一定的响应速度。能够满足频率偏移、谐波、幅值突变、相位突变、频率斜坡变化、幅值相位调制等多种情况下的对于相量的高精度测量，可用于同步相量装置的核心测量算法或其它需要计算相量的场合。

Description

一种基于DFT、多阶滤波和突变判据的同步相量计算方法
技术领域
本发明涉及电力系统同步相量的计算方法技术领域，尤其涉及一种基于DFT、多阶滤波和突变判据的同步相量计算方法，具体是一种基于DFT、采用多阶数滤波、突变判据计算的同步相量计算方法。
背景技术
电力系统的不断发展使电网呈现出电压等级更高、传输距离更远、输送容量更大的特点。同时，超高压远距离交直流混合输电、广泛应用新技术和市场化运营机制等，也使系统的运行方式更加复杂，这对电力系统的安全性提出了新的挑战。为尽量避免电力事故的发生造成损失，研究电网的稳定性监测、异常状态判断、发展趋势预测、有效的保护控制策略等技术手段，就尤其重要。
基于全球定位系统(GPS)的同步相量测量装置(PMU)涉及到电力系统的监视、控制和保护等诸多领域，推动了电力系统的监视、控制和保护新方法和理论的发展，为电力系统的稳定控制和保护开辟一个新的领域，已成为电力系统中一项重要的量测技术。同步相量测量数据可用于预测系统状态、提高状态估计精度、加快潮流计算速度、输电线路的故障定位、发电机非线性励磁控制、阻尼系统振荡的SVC控制等，因此对于提高对电力系统动态稳定监测和故障时系统暂态过程的分析能力具有重要意义。
同步相量测量的基本功能是利用GPS信号对电网中母线或节点电压/电流同步测量，并进行分析计算，可以测量出本地频率、频率变化率、正序相量，或通过定制测量得到谐波、零序量、相电压、相电流等信息。对同步相量测量装置来说，同步相量测量算法是装置的核心，也是研究的重点，它主要指对于幅值、相角、频率、频率变化率的测量方法。相量测量算法的优劣，直接决定了同步相量装置的成本、性能和效用。
目前，常用的相量测量方法主要有过零检测法、定频采样的离散傅立叶算法(DFT)和变频采样的DFT方法。为了减小计算量或提升计算精度，各种改进型算法被设计出来，如递归式DFT、基于锁相环进行频率跟踪的DFT。此外，还有非DFT的相量计算方法，如神经网络、卡尔曼滤波、小波分析及混合算法等，不同算法各有优劣。
相量测量算法一直是研究中的重点与难点，这是因为，电网的频率不可能完全精确地运行在额定频率点，根据DFT的计算原理，在这样的情况下会出现频谱泄露现象，从而导致较大的计算误差。此外，电网电压和电流还可能包含各种整次、非整次谐波和其他暂态分量。由于故障等各种原因，电力系统可能产生各种暂态过程，如幅值突变、相位突变、频率斜坡变化、幅值和相位调制等，从而导致电压和电流波形中含有暂态分量。很多在稳态测量中精度较高的算法，对于暂态变化情况下的测量往往都表现不佳。因此，研究具有响应速度快、计算量小、误差小、适用范围广的相量测量算法，对于推动同步向量装置在电网的应用具有重要意义。
发明内容
针对上述现有技术中存在原不中之处，本发明提出了一种基于DFT、多阶滤波和突变判据的同步相量计算方法，其目的在于提出一种同步相量计算的新方法，该方法以DFT算法为基础，通过设计多阶滤波器的组合设计，在满足响应速度要求的前提下，尽量提升了计算精度。提出了波形突变判据，给出了波形突变的阈值大小，用以判断待测信号的稳定状态。若发生了波形突变，则通过算法对计算结果进行修正。对于稳态测量和暂态测量两种状态的区分，有效地提升了算法的适用范围和计算精度。
为了实现上述发明目的，本发明是通过以下方式实现的：
一种基于DFT、多阶滤波和突变判据的同步相量计算方法，包括以下步骤：
步骤1.初始计算：对同步采样数据采用DFT算法，计算出待测信号的幅值x1和相位φ1；
步骤2.频率计算：通过相邻时刻的相位差Δφ1计算出待测信号频率f1，通过三点滤波、定点数平均滤波、不定点数的平均滤波的多阶滤波器的组合，对计算得到的频率值进行滤波处理，能够得到待测信号的频率fc；在这个过程中，还会得到一个中间计算值f2，用于检测到波形发生突变时的相关计算；
步骤3.幅值与相位的计算：稳态情况，采用频率fc对同步采样数据进行再样，能够得到待测信号的幅值xc和相位φc；暂态情况，采用f2对同步数据进行再采样，得到相位p3作为相位的计算值位φc；
步骤4.频率变化率的计算：根据额定频率、频率偏移量、频率变化率和相位之间的关系，采用最小二乘法计算频率变化率f′；
步骤5.信号突变判断算法：根据波形突变会导致DFT计算结果发生突变的原理设计了突变判据J，通过仿真计算得到突变阈值大小Jmax；一旦J>Jmax，则认为波形发生了突变；此时，改变突变计数值C的大小；算法流程进入突变状态，并修正计算结果；每次循环修正后C值减1，当C等于0时且J<Jmax，认为待测信号再次趋于稳定，算法恢复至稳态测量状态；
步骤6.算法测试。
所述步骤1.初始计算，包括：
相量是电力系统稳态条件下定义的概念；一个纯正弦信号表示为：x(t)＝Xmcos(ωt+φ)，根据欧拉公式，表示为，x(t)＝Re{ej(ωt+φ)}，若忽略ejωt，则正弦信号用一个复数X表示，称之为信号的相量表示形式，即，X＝Xm‘ejφ＝Xm‘(cosφ+jsinφ)；
初始计算的算法令fs为采样率，对于额定频率为50Hz的正弦信号，一周期能够采样n个点；xk为k时刻的待测正弦信号的采样值；经过一段时间的采样，得到t＝-(k-1)，t＝-(k-2)，……，t＝-2，t＝-1，t＝0时刻n个采样值；采样数据以堆栈的方式存储于缓存数组W1中；其中，W1[0]存储了最新时刻的采样值，W1[1]存储了上一时刻的采样值，……以此类推；所有缓存数组皆以这种方式对某次计算得到的数值进行存储；采用DFT算法，通过计算得到待测信号的幅值Xm与相位φ；计算公式为：
Figure BDA0002095454990000031
上式中：xk表示在k时刻对待测信号进行采样，得到的信号瞬时值；a表示经过DFT计算得到的相量实部；b表示经过DFT计算得到的相量虚部；Xm表示经过DFT计算得到的信号峰值；φ'表示相量实部与虚部相比后得到的反正切角度值；
反正切函数arctan的值域范围为[-π/2,π/2]，而相量的角度范围为(-π,π]，故需要判断(a,b)所在象限，经过转换后得到待测信号的相角φ，相应计算公式为：
Figure BDA0002095454990000032
上式中：a表示经过DFT计算得到的相量实部；b表示经过DFT计算得到的相量虚部；φ为经过象限转换后得到的待测信号的相角；
相位与幅值的计算结果以堆栈方式存储于一定长度的缓存数组P1与X1中。
所述步骤2.频率计算，包括：
设k时刻、k+1时刻通过初始计算得到的相位φk、φk+1分别存储于缓存P1[1]和P1[0]中，待测信号的频率f1的计算公式为：
Figure BDA0002095454990000041
上式中：fs为对待测信号进行模数转换的采样频率；P为存储数据，括号内数字为数组下标；f1为用于频率计算的中间参数；
实际电力系统中，信号的频率会存在不同程度的频率偏移；当待测信号频率偏移额定频率时，DFT计算得到幅值和相位会含有近似的二次谐波；基于相位信息计算得到的频率f1，同样也会耦合二次谐波；除了计算误差，电网本身也存在着谐波，此外，信号还会受到各类噪声的影响；故计算得到的f1误差较大，不能直接作为待测信号的频率值；
为了提升频率计算的精度，需要进行滤波处理；滤波算法采用多阶结构，由不同的滤波算法组合而成；滤波算法包括三点滤波、递归式的定点数平均滤波、递归式的变点数的平均滤波；
对比三种滤波算法的优缺点，三点滤波算法计算量小、延时小，但是滤波效果一般；递归式定点数平均滤波算法计算量小，延时相对较长，滤波效果较好；递归式变点数平均滤波算法计算量小，延时相对较长，滤波效果最好，但是需要以获知待测信号大致频率为计算的前提条件；多个滤波器级联，能够构成多阶滤波算法；滤波效果若要越好，则滤波阶数也越高，同时滤波输出对于输入信号的变化的响应也就越慢；滤波算法的设计要兼顾精度与响应速度；
采用本发明算法，对THD谐波含量为0.2％、频率为46Hz的稳态正线信号进行频率测量,经过DFT计算得到相位后再直接计算得到f1由于待测信号频率为46Hz，偏移了50Hz的额定频率，故f1耦合了二次谐波，在40Hz与50Hz之间进行波动，误差较大；经过多阶滤波后，依次得到信号频率f3、f4、f5、f6，频率波动的幅度也依次减小；f3～f6对应的最大误差分别为400mHz、40mHz、5mH、0.06mHz；经过多阶滤波后得到的f6误差较其小，作为频率测量结果fc。
所述三点滤波算法的原理为，对相位相差为0°、60°、120°的三个量进行相加，能够在一定程度上滤除二次脉动谐波；对于待测信号额定频率为50Hz，采样率fs为1kHz的情况，若数组X连续存储了t＝0时刻、t＝-1时刻、t＝-2时刻、……、t＝-n时刻的计算结果，则相差0°、60°、120°的三个数值依次为：X(0)、X(3)、X(7)；对其求平均即能够滤除二次谐波，对应的计算公式为：
Figure BDA0002095454990000051
上式中：X为存储数据，括号内数字为数组下标；out为滤波计算后得到的输出结果。
所述递归式的定点数平均滤波原理为，对固定长度的数据进行平均化处理；对于长度为n的数组一维X进行平均滤波的计算公式为，
Figure BDA0002095454990000052
上式中：X为存储数据，括号内数字为数组下标；n为数组X的长度；out为滤波计算后得到的输出结果；
采用递归式的定点数平均滤波算法，有效减小滤波算法计算量；令数组sum为长度为n的数组X各成员数据之和，最新数据为xnew，最旧数据为xold，则平均滤波结果avg为，
Figure BDA0002095454990000053
上式中：sum为X数组所有数据之和；n为数组X的长度；avg为数组X的所有数据之和的平均值；
每次计算avg完成后，新数据xnew进栈，旧数据xold出栈剔除，以实现数据的不断刷新存储。
所述递归式的变点数的平均滤波，是指：当待测信号频率偏移额定频率时，频率计算结果耦合的谐波实际上只是近似于二次谐波，故采用固定点数的平均滤波效果仍然有限；采用频率跟踪的方式，进行变点数的平均滤波具有更好的滤波效果；每次滤波求取平均的点数n1的计算公式为，
Figure BDA0002095454990000061
上式中：fs为对待测信号进行模数转换的采样频率；f为通过初始计算得到的待测信号的频率，ceil函数为求不小于给定数的最小整数，n1为一采样周期内采样数据个数；
其中，ceil(x)函数的作用为，求不小于给定实数x的最小整数；fs为采样率，f为通过初始计算得到的待测信号的频率；与平均滤波算法类似，变点数滤波亦为递归形式。
所述步骤3.幅值与相位的计算，包括：
稳态情况，采用频率fc对同步采样数据进行再样，能够得到待测信号的幅值xc和相位φc；暂态情况，采用f2对同步数据进行再采样，得到相位p3作为相位的计算值位φc；
当待测信号偏移了额定频率时，还按照额定频率的DFT进行幅值和相位计算，是造成计算产生较大误差的原因；此时需要采用再采样算法，等价于DFT计算频率与待测信号频率相同，因而能够得到精度较高的幅值与相位计算结果；再采样算法原理具体如下:
设偏移了额定频率的输入信号为:
x(t)＝Xm cos(ωt+φ)
上式中：x(t)为待测信号t时刻瞬时值；Xm为待测信号峰值；ω为待测信号角频率；φ为待测信号初始相角；
按照额定频率的整N倍数频率Nω0进行采样，第k个采样点和k+1个采样值为：
Figure BDA0002095454990000062
上式中：xk为离散化后，k时刻对待测信号进行采样得到的信号瞬时值；xk+1为离散化后k+1时刻对待测信号进行采样得到的信号瞬时值；Xm为待测信号峰值；θ为k+1时刻与k时刻相比，待测信号的相位变化大小；φ为待测信号初始相角；
若按照与待测信号实际频率ω对应的整N倍数频率Nω进行再采样，则第n个点能够表示为：
xn＝Xm cos(kθ+γ+φ)
上式中：xn为采用与待测信号频率整倍数采样率，对待测信号进行采样得到的瞬时值；Xm为待测信号峰值；k为对时间进行离散化后的时刻值，φ为待测信号初始相角；γ为再采样后的信号与最近时刻采样信号之间的相位差大小；
若xn对应的采样时间介于xk与xk+1之间，推导xk、xk+1、xn之间的关系有：
Figure BDA0002095454990000071
上式中：xk为离散化后k时刻对待测信号进行采样得到的信号瞬时值，xk+1为离散化后k+1时刻对待测信号进行采样得到的信号瞬时值；θ为k+1时刻与k时刻间待测信号的相位变化；γ为再采样后的频率偏移量大小；
其中，f为再采样前经过额定频率DFT计算得到的待测信号频率，fs为采样频率，Δt为xn与xk采样点之间的时间差；θ＝2πf/fs，γ＝2πfΔt；W1经过再采样后得到新的采样序列W2，对W2再进行一次DFT计算，得到待测信号幅值x4、相位p2；x4存入缓存X4，对X4进行一次FIR低通滤波，得到的幅值计算结果Xc具有较高的精度；
当待测信号稳定时，信号突变计数C为0，此时采用p2作为相位计算结果；而当波形发生幅值突变或相位突变时，需要快速计算波形相位；此时采用响应较快的频率f2作为再采样算法的采样频率，随之计算得到的相位p3，为相位计算结果。
所述步骤4.频率变化率的计算，包括：
根据额定频率、频率偏移量、频率变化率和相位之间的关系，采用最小二乘法计算频率变化率f′；
其原理为，任意时刻t的角频率ω表示为：
ω(t)＝ω0+Δω+tω′
上式中：ω(t)为t时刻待测信号的瞬时角频率，ω0为额定角频率，Δω为偏移额定角频率的大小，ω′为角频率变化率；
相角φ是频率相对于时间的积分，有：
Figure BDA0002095454990000072
上式中：φ(t)为待测信号t时刻相位瞬时值，φ0为初始相位大小，ω0为额定角频率，Δω为偏移额定角频率的大小，ω′为角频率变化率；
可知φ(t)是关于t的二阶多项式函数；通过运算得到信号的相位角，将计算得到的n个相位角组成一维列向量：
[Φ]＝[φ0 φ1 φ2…φn-1]T
上式中：φx为第x次计算得到的待测信号相位大小；Φ为列向量，用于存储n个相位计算结果；
根据相位与频率的关系，列出方程：
Figure BDA0002095454990000081
上式中φx为第x次计算得到的待测信号相位大小，Φ为列向量，用于存储n个相位计算结果，Δt为两次采样间的时间间隔，ω0为额定角频率，Δω为偏移额定角频率的大小，ω’为角频率变化率；将上式简化为矩阵形式有[Φ]＝[B][A]；其中B为3列n行的矩阵，与采样间隔时间有关；A为1列3行的矩阵，称之为角频率特征参数矩阵；采用矩阵A、B用于简化表达式；
根据最小二乘法，计算得到矩阵A为：
Figure BDA0002095454990000082
上式中：B为简化后的矩阵，A为角频率特征参数矩阵，ω0为额定角频率，Δω为偏移额定角频率的大小，ω′为角频率变化率，Φ为列向量，用于存储n个相位计算结果；
求解A矩阵得到ω0、Δω、ω′，进一步计算频率变化率f′，有：
f′＝ω′/2π
上式中：f′为计算得到的待测信号频率变化率，ω′为待测信号角频率变化率；
当采样率恒定时，B矩阵就固定了，因此每次计算前不用再计算B，而是在初始化前程序时直接令：
G＝[BTB]-1B
上式中：G为n阶方矩，B为3列n行的矩阵，与采样间隔时间有关；
由于相位角φ的值域为[-π，π)，在计算之前，[Φ]需要进行单调化处理，以保证φk<φk+1(k＝0，1，2……，n-2)。
所述步骤5.信号突变判断算法，包括：
经过一轮计算后，得到信号的幅值x、相位φ、频率f、频率变化率r；由于在计算过程中采用了滤波算法，因此计算得到的结果相比原信号波形具有一定的延时；(a)为待测波形，圆圈标示的为采样值；待测信号为稳定的周期信号；(b)为算法计算波形，相比待测信号，计算值延时了k个单位；为了去除滤波掩饰的影响，在为数据打时标时，需要等价于将计算波形前移k个单位；前移后计算波形与待测波形相重合，计算误差较小；
然而，这样的波形前移的方式只适用于待测波形为稳态周期信号的情况，当波形发生突变时，波形前移会导致误差；(c)所示为在t＝i时刻发生突变的待测波形；由于滤波的效果，计算得到的波形(d)；若采取与稳态计算一样的算法，将计算波形提前，则(e)；待测波形与计算波形相比，额外多出的阴影部分，即为计算误差；此部分的误差还意味着，计算得到的波形突变时刻，与实际情况相比提前了k个时间长度；
为了能够区分稳态波形和突变波形两种状态、采取不同的数据处理方式以最小化计算误差，需要设计波形突变判据；根据DFT的计算原理，当信号发生突变时，若计算的时间窗了包含了突变数据，则DFT的计算结果也会发生相应的突变；故能够根据信号幅值的变化大小，作为信号发生突变的判据J，一旦J大于了某个阈值J，即认为待测信号发生了突变；突变判据J的计算公式为，
Figure BDA0002095454990000091
上式中：J为信号突变判据，X1为信号幅值计算存储数组，括号内数字为数组下标，abs为求绝对值函数；
公式中k的取值与采样率有关，一般要求X1(k)与X1(0)两者所对应的时间差要在8ms左右，这样计算得到的J大小合适；若k太大，则计算得到的信号突变时间点与实际信号突变时间点误差较大；若k太小，则J也较小，有可能无法判断信号突变是否发生；为了确定阈值Jt的大小，对各种情况下的J进行了测试，只有当发生幅值突变和相位突变的情况下，J大于1；故取阈值Jt＝1.0，当J>Jt时，则认为波形发生了突变；
当发生突变后，通过判断J的大小能够知道何时波形会再趋于稳定；但是由于算法结构中包含滤波器，故当突变后J小于一定数值时，计算结果可能仍未稳定；因此，需要采用突变计数C来记录波形突变的长度；波形未突变时，C＝0；当波形突变后，C的大小一般等于一周期采样点数n的7.5～7.8倍，令这个值为Ck；当C≠0时，每次完成采样计算后，C需要减1，当C再次等于0时，且J<1，意味着信号再次稳定；否则，C＝Ck，待测波形仍为不稳定的变化状态。
所述步骤6.算法测试，包括：频率偏移、谐波、幅值突变、相位突变、频率斜坡变化、幅值调制、相位调制及幅值相位同时调制；测试考察幅值、相位、频率、频率变化率、正序幅值以及正序相位。
本发明具有以优点及有益效果：
电力系统的同步相量测量属于高精度测量。电力系统可能出现的状况较为复杂。因而对于同步相量算法的测试内容较多，包含频率偏移、谐波、幅值突变、相位突变、频率斜坡变化、幅值相位调制等。同步相量算法往往不能兼顾，有可能对于其中一种或几种情况计算精度较高，而对于其它情况则计算精度较低。本发明提出一种同步相量的计算方法。通过突变判据的计算将幅值突变、相位突变与其它几种稳态或准稳态情况区别开来处理。通过多阶滤波器的设计，提升计算的精度，同时又具备一定的响应速度。本发明的算法能够满足频率偏移、谐波、幅值突变、相位突变、频率斜坡变化、幅值相位调制等多种情况下的对于相量的高精度测量，可用于同步相量装置的核心测量算法，或是其它需要计算相量的场合。
附图说明
图1是本发明的相量算法结构示意图；
图2是本发明的初始计算算法结构示意图；
图3是本发明的频率计算算法结构示意图；
图4、图5是本发明的多阶滤波波形图；
图6是本发明的幅值计算算法结构示意图；
图7是本发明的再采样算法原理示意图；
图8是本发明的频率变化率计算算法结构示意图；
图9是本发明的相位角的未处理示意图；
图10是本发明的相位角的单调化处理示意图；
图11是本发明的计算延时稳态信号示意图；
图12是本发明的计算延时突变信号示意图；
图13是本发明的种情况下突变判据Jmax数值大小对比图；
图14是本发明的波形突变判据J与突变计数C计算算法示意图。
具体实施方式
本发明提出了一种基于DFT、多阶滤波和突变判据的同步相量计算方法，如图1所示，图1是本发明的相量算法结构示意图。该算法主要由初始计算、波形突变计算、频率计算、幅值相位计算、频率变化率计算以及计算结果修正几部分组成。算法以DFT离散傅里叶变换为基础，设计了特殊的多阶滤波算法，以提升计算精度。通过计算波形突变判据J的数值大小，对计算结果进行修正，以解决正弦信号稳态情况与波形发生突变情况下，相量的高精度计算问题。下面详细介绍算法各部分原理。
一种基于DFT、多阶滤波和突变判据的同步相量计算方法，包括以下步骤：
步骤1.初始计算。
对同步采样数据采用DFT算法，计算出待测信号的幅值x1和相位φ1。
相量是电力系统稳态条件下定义的概念。一个纯正弦信号可以表示为：x(t)＝Xmcos(ωt+φ)，根据欧拉公式，可以表示为，x(t)＝Re{ej(ωt+φ)}，若忽略ejωt，则正弦信号可用一个复数X表示，称之为信号的相量表示形式，即，X＝Xm‘ejφ＝Xm‘(cosφ+jsinφ)。
初始计算的算法结构如图2所示。令fs为采样率，对于额定频率为50Hz的正弦信号，一周期可以采样n个点。xk为k时刻的待测正弦信号的采样值。经过一段时间的采样，得到t＝-(k-1)，t＝-(k-2)，……，t＝-2，t＝-1，t＝0时刻n个采样值。采样数据以堆栈的方式存储于缓存数组W1中。其中，W1[0]存储了最新时刻的采样值，W1[1]存储了上一时刻的采样值，……以此类推。若不加特殊说明，本发明中所有缓存数组皆以这种方式对某次计算得到的数值进行存储。采用DFT算法，通过计算可以得到待测信号的幅值Xm与相位φ。计算公式为：
Figure BDA0002095454990000121
上式中：xk表示在k时刻对待测信号进行采样，得到的信号瞬时值；a表示经过DFT计算得到的相量实部；b表示经过DFT计算得到的相量虚部；Xm表示经过DFT计算得到的信号峰值；φ'表示相量实部与虚部相比后得到的反正切角度值。
反正切函数arctan的值域范围为[-π/2,π/2]，而相量的角度范围为(-π,π]，故需要判断(a,b)所在象限，经过转换后得到待测信号的相角φ，相应计算公式为：
Figure BDA0002095454990000122
上式中：a表示经过DFT计算得到的相量实部；b表示经过DFT计算得到的相量虚部；φ为经过象限转换后得到的待测信号的相角。
相位与幅值的计算结果以堆栈方式存储于一定长度的缓存数组P1与X1中。
步骤2.频率计算。
通过相邻时刻的相位差Δφ1计算出待测信号频率f1，通过三点滤波、定点数平均滤波、不定点数的平均滤波的多阶滤波器的组合，对计算得到的频率值进行滤波处理，可以得到待测信号的频率fc。在这个过程中，还会得到一个中间计算值f2，用于检测到波形发生突变时的相关计算。具体如下：
设k时刻、k+1时刻通过初始计算得到的相位φk、φk+1分别存储于缓存P1[1]和P1[0]中，待测信号的频率f1的计算公式为，
Figure BDA0002095454990000131
上式中：fs为对待测信号进行模数转换的采样频率；P为存储数据，括号内数字为数组下标；f1为用于频率计算的中间参数。
实际电力系统中，信号的频率不会恰好是50Hz，会存在不同程度的频率偏移。当待测信号频率偏移额定频率时，DFT计算得到幅值和相位会含有近似的二次谐波。基于相位信息计算得到的频率f1，同样也会耦合二次谐波。除了计算误差，电网本身也存在着谐波，此外，信号还会受到各类噪声的影响。故计算得到的f1误差较大，不能直接作为待测信号的频率值。
为了提升频率计算的精度，需要进行滤波处理。滤波算法采用多阶结构，由不同的滤波算法组合而成。滤波算法包括三点滤波、递归式的定点数平均滤波、递归式的变点数的平均滤波。
三点滤波算法的原理为，对相位相差为0°、60°、120°的三个量进行相加，可以在一定程度上滤除二次脉动谐波。对于待测信号额定频率为50Hz，采样率fs为1kHz的情况，若数组X连续存储了t＝0时刻、t＝-1时刻、t＝-2时刻、……、t＝-n时刻的计算结果，则相差0°、60°、120°的三个数值依次为：X(0)、X(3)、X(7)。对其求平均即可以滤除二次谐波，对应的计算公式为：
Figure BDA0002095454990000132
上式中：X为存储数据，括号内数字为数组下标；out为滤波计算后得到的输出结果。
递归式的定点数平均滤波原理为，对固定长度的数据进行平均化处理。例如，对于长度为n的数组一维X进行平均滤波的计算公式为：
Figure BDA0002095454990000133
上式中：X为存储数据，括号内数字为数组下标；n为数组X的长度；out为滤波计算后得到的输出结果。
采用递归式的定点数平均滤波算法，可有效减小滤波算法计算量。令数组sum为长度为n的数组X各成员数据之和，最新数据为xnew，最旧数据为xold，则平均滤波结果avg为，
Figure BDA0002095454990000141
上式中：sum为X数组所有数据之和；n为数组X的长度；avg为数组X的所有数据之和的平均值。
每次计算avg完成后，新数据xnew进栈，旧数据xold出栈剔除，以实现数据的不断刷新存储。
当待测信号频率偏移额定频率时，频率计算结果耦合的谐波实际上只是近似于二次谐波，故采用固定点数的平均滤波效果仍然有限。采用频率跟踪的方式，进行变点数的平均滤波具有更好的滤波效果。每次滤波求取平均的点数n1的计算公式为，
Figure BDA0002095454990000142
上式中：fs为对待测信号进行模数转换的采样频率；f为通过初始计算得到的待测信号的频率，ceil函数为求不小于给定数的最小整数，n1为一采样周期内采样数据个数。
其中，ceil(x)函数的作用为，求不小于给定实数x的最小整数。fs为采样率，f为通过初始计算得到的待测信号的频率。与平均滤波算法类似，变点数滤波亦为递归形式。
对比三种滤波算法的优缺点。三点滤波算法计算量小、延时小，但是滤波效果一般；递归式定点数平均滤波算法计算量小，延时相对较长，滤波效果较好；递归式变点数平均滤波算法计算量小，延时相对较长，滤波效果最好，但是需要以获知待测信号大致频率为计算的前提条件。多个滤波器级联，可以构成多阶滤波算法。滤波效果若要越好，则滤波阶数也越高，同时滤波输出对于输入信号的变化的响应也就越慢。滤波算法的设计要兼顾精度与响应速度。最终采用的多阶滤波算法如图3所示，图3是本发明的频率计算算法结构示意图；
如图4、图5所示,图4、图5是本发明的多阶滤波波形图。为采用本发明算法，对THD谐波含量为0.2％、频率为46Hz的稳态正线信号进行频率测量的波形图。经过DFT计算得到相位后再直接计算得到f1。由图4可知，由于待测信号频率为46Hz，偏移了50Hz的额定频率，故f1耦合了二次谐波，在40Hz与50Hz之间进行波动，误差较大。经过多阶滤波后，依次得到信号频率f3、f4、f5、f6，频率波动的幅度也依次减小。f3～f6对应的最大误差分别为400mHz、40mHz、5mH、0.06mHz。经过多阶滤波后得到的f6误差较其小，作为频率测量结果fc。图3所示算法中的频率f2用于信号突变时的相关计算。
步骤3.幅值与相位计算。
稳态情况，采用频率fc对同步采样数据进行再样，可以得到待测信号的幅值xc和相位φc；暂态情况，采用f2对同步数据进行再采样，得到相位p3作为相位的计算值位φc。
如图6所示，图6是本发明的幅值计算算法结构示意图。当待测信号偏移了额定频率时，还按照额定频率的DFT进行幅值和相位计算，是造成计算产生较大误差的原因。此时需要采用再采样算法，等价于DFT计算频率与待测信号频率相同，因而可以得到精度较高的幅值与相位计算结果。再采样算法原理具体如下。
如图7所示，图7是本发明的再采样算法原理示意图
设偏移了额定频率的输入信号为：
x(t)＝Xm cos(ωt+φ)
上式中：x(t)为待测信号t时刻瞬时值；Xm为待测信号峰值；ω为待测信号角频率；φ为待测信号初始相角。
按照额定频率的整N倍数频率Nω0进行采样，第k个采样点和k+1个采样值为：
Figure BDA0002095454990000151
上式中：xk为离散化后，k时刻对待测信号进行采样得到的信号瞬时值；xk+1为离散化后k+1时刻对待测信号进行采样得到的信号瞬时值；Xm为待测信号峰值；θ为k+1时刻与k时刻相比，待测信号的相位变化大小；φ为待测信号初始相角。
若按照与待测信号实际频率ω对应的整N倍数频率Nω进行再采样，则第n个点可以表示为：
xn＝Xm cos(kθ+γ+φ)
上式中：xn为采用与待测信号频率整倍数采样率，对待测信号进行采样得到的瞬时值；Xm为待测信号峰值；k为对时间进行离散化后的时刻值，φ为待测信号初始相角；γ为再采样后的信号与最近时刻采样信号之间的相位差大小。
若xn对应的采样时间介于xk与xk+1之间，推导xk、xk+1、xn之间的关系有：
Figure BDA0002095454990000161
上式中：xk为离散化后k时刻对待测信号进行采样得到的信号瞬时值，xk+1为离散化后k+1时刻对待测信号进行采样得到的信号瞬时值；θ为k+1时刻与k时刻间待测信号的相位变化；γ为再采样后的频率偏移量大小。
其中，f为再采样前经过额定频率DFT计算得到的待测信号频率，fs为采样频率，Δt为xn与xk采样点之间的时间差。θ＝2πf/fs，γ＝2πfΔt。W1经过再采样后得到新的采样序列W2，对W2再进行一次DFT计算，得到待测信号幅值x4、相位p2。x4存入缓存X4，对X4进行一次FIR低通滤波，得到的幅值计算结果Xc具有较高的精度。
当待测信号稳定时，信号突变计数C为0，此时采用p2作为相位计算结果。而当波形发生幅值突变或相位突变时，需要快速计算波形相位。此时采用响应较快的频率f2作为再采样算法的采样频率，随之计算得到的相位p3，为相位计算结果。
步骤4.频率变化率的计算。
根据额定频率、频率偏移量、频率变化率和相位之间的关系，采用最小二乘法计算频率变化率f′。
如图8所示，图8是本发明的频率变化率计算算法结构示意图。
其原理为，任意时刻t的角频率ω可以表示为：
ω(t)＝ω0+Δω+tω′
上式中：ω(t)为t时刻待测信号的瞬时角频率，ω0为额定角频率，Δω为偏移额定角频率的大小，ω′为角频率变化率。
相角φ是频率相对于时间的积分，有：
Figure BDA0002095454990000162
上式中：φ(t)为待测信号t时刻相位瞬时值，φ0为初始相位大小，ω0为额定角频率，Δω为偏移额定角频率的大小，ω′为角频率变化率。
可知φ(t)是关于t的二阶多项式函数。通过运算可以得到信号的相位角，将计算得到的n个相位角组成一维列向量：
[Φ]＝[φ0 φ1 φ2…φn-1]T
上式中：φx为第x次计算得到的待测信号相位大小。Φ为列向量，用于存储n个相位计算结果。
根据相位与频率的关系，列出方程：
Figure BDA0002095454990000171
上式中φx为第x次计算得到的待测信号相位大小。Φ为列向量，用于存储n个相位计算结果。Δt为两次采样间的时间间隔。ω0为额定角频率，Δω为偏移额定角频率的大小，ω’为角频率变化率。将上式简化为矩阵形式有[Φ]＝[B][A]。其中B为3列n行的矩阵，与采样间隔时间有关；A为1列3行的矩阵，称之为角频率特征参数矩阵。采用矩阵A、B可用于简化表达式。
根据最小二乘法，计算得到矩阵A为：
Figure BDA0002095454990000172
上式中：B为时间矩阵，A为角频率特征参数矩阵。ω0为额定角频率，Δω为偏移额定角频率的大小，ω′为角频率变化率。Φ为列向量，用于存储n个相位计算结果。
求解A矩阵得到ω0、Δω、ω′，进一步计算频率变化率f′，有：
f′＝ω′/2π
上式中：f′为计算得到的待测信号频率变化率，ω′为计算得到的待测信号角频率变化率。
当采样率恒定时，B矩阵就固定了，因此每次计算前不用再计算B，而是在初始化前程序时直接令：
G＝[BTB]-1B
上式中：B为3列n行的矩阵，与采样间隔时间有关；G为n列3行的矩阵。
由于相位角φ的值域为[-π，π)，在计算之前，[Φ]需要进行单调化处理，以保证φk<φk+1(k＝0，1，2……，n-2)。对相角进行单调化的前后对比如9和图10所示，图9是本发明的相位角的未处理示意图，图10是本发明的相位角的单调化处理示意图。
步骤5.信号突变判断算法。
根据波形突变会导致DFT计算结果发生突变的原理设计了突变判据J，通过仿真计算得到突变阈值大小Jmax。一旦J>Jmax，则认为波形发生了突变。此时，改变突变计数值C的大小。算法流程进入突变状态，并修正计算结果。每次循环修正后C值减1，当C等于0时且J<Jmax，认为待测信号再次趋于稳定，算法恢复至稳态测量状态。具体如下：
经过一轮计算后，可以得到信号的幅值x、相位φ、频率f、频率变化率r。由于在计算过程中采用了滤波算法，因此计算得到的结果相比原信号波形具有一定的延时。如图11所示，图11是本发明的计算延时稳态信号示意图。图11中(a)为待测波形，圆圈标示的为采样值。待测信号为稳定的周期信号。图11中(b)为算法计算波形，相比待测信号，计算值延时了k个单位。为了去除滤波掩饰的影响，在为数据打时标时，需要等价于将计算波形前移k个单位。前移后计算波形与待测波形相重合，计算误差较小。
然而，这样的波形前移的方式只适用于待测波形为稳态周期信号的情况，当波形发生突变时，波形前移会导致误差。如图12所示，图12是本发明的计算延时突变信号示意图。图中(c)所示为在t＝i时刻发生突变的待测波形。由于滤波的效果，计算得到的波形如图中(d)所示为阻尼变化的波形图中为简单示意，去掉了阻力变化的波形细节。若采取与稳态计算一样的算法，将计算波形提前，则如图中(e)所示。待测波形与计算波形相比，额外多出的阴影部分，即为计算误差。此部分的误差还意味着，计算得到的波形突变时刻，与实际情况相比提前了k个时间长度。
为了能够区分稳态波形和突变波形两种状态、采取不同的数据处理方式以最小化计算误差，需要设计波形突变判据。根据DFT的计算原理，当信号发生突变时，若计算的时间窗了包含了突变数据，则DFT的计算结果也会发生相应的突变。故可以根据信号幅值的变化大小，作为信号发生突变的判据J，一旦J大于了某个阈值J，即认为待测信号发生了突变。突变判据J的计算公式为，
Figure BDA0002095454990000191
上式中：J为信号突变判据，X1为信号幅值计算存储数组；括号内数字为数组下标。abs为求绝对值函数。
公式中k的取值与采样率有关，一般要求X1(k)与X1(0)两者所对应的时间差要在8ms左右，这样计算得到的J大小合适。若k太大，则计算得到的信号突变时间点与实际信号突变时间点误差较大。若k太小，则J也较小，有可能无法判断信号突变是否发生。为了确定阈值Jt的大小，对各种情况下的J进行了测试，测试对比结果如图13所示，图13是本发明的种情况下突变判据Jmax数值大小对比图。由图13可知，只有当发生幅值突变和相位突变的情况下，J大于1。故可以取阈值Jt＝1.0，当J>Jt时，则认为波形发生了突变。
当发生突变后，通过判断J的大小可以知道何时波形会再趋于稳定。但是由于算法结构中包含滤波器，故当突变后J小于一定数值时，计算结果可能仍未稳定。因此，需要采用突变计数C来记录波形突变的长度。波形未突变时，C＝0；当波形突变后，C的大小一般等于一周期采样点数n的7.5～7.8倍，令这个值为Ck。当C≠0时，每次完成采样计算后，C需要减1，当C再次等于0时，且J<1，意味着信号再次稳定。否则，C＝Ck，待测波形仍为不稳定的变化状态。如图14所示，图14是本发明的波形突变判据J与突变计数C计算算法示意图。
步骤6.算法测试。
相量测试主要的内容有：频率偏移、谐波、幅值突变、相位突变、频率斜坡变化、幅值调制、相位调制、幅值相位同时调制等。测试主要考察幅值、相位、频率、频率变化率、正序幅值、正序相位这些参数的计算误差需要满足一定要求。表1汇总了对于本文相量设计方法的各项仿真测试结果，均满足同步相量测量的标准要求。
表1各项仿真测试结果一览表。
Figure BDA0002095454990000192
Figure BDA0002095454990000201

Claims (10)

1.一种基于DFT、多阶滤波和突变判据的同步相量计算方法，其特征在于：包括以下步骤：
步骤1.初始计算：对同步采样数据采用DFT算法，计算出待测信号的幅值x1和相位φ1；
步骤2.频率计算：通过相邻时刻的相位差Δφ1计算出待测信号频率f1，通过三点滤波、定点数平均滤波、不定点数的平均滤波的多阶滤波器的组合，对计算得到的频率值进行滤波处理，能够得到待测信号的频率fc；在这个过程中，还会得到一个中间计算值f2，用于检测到波形发生突变时的相关计算；
步骤3.幅值与相位的计算：稳态情况，采用频率fc对同步采样数据进行再样，能够得到待测信号的幅值xc和相位φc；暂态情况，采用f2对同步数据进行再采样，得到相位p3作为相位的计算值位φc；
步骤4.频率变化率的计算：根据额定频率、频率偏移量、频率变化率和相位之间的关系，采用最小二乘法计算频率变化率f′；
步骤5.信号突变判断算法：根据波形突变会导致DFT计算结果发生突变的原理设计了突变判据J，通过仿真计算得到突变阈值大小Jmax；一旦J>Jmax，则认为波形发生了突变；此时，改变突变计数值C的大小；算法流程进入突变状态，并修正计算结果；每次循环修正后C值减1，当C等于0时且J<Jmax，认为待测信号再次趋于稳定，算法恢复至稳态测量状态；
步骤6.算法测试。
2.根据权利要求1所述的一种基于DFT、多阶滤波和突变判据的同步相量计算方法，其特征在于：所述步骤1.初始计算，包括：
相量是电力系统稳态条件下定义的概念；一个纯正弦信号表示为：x(t)＝Xmcos(ωt+φ)，根据欧拉公式，表示为，x(t)＝Re{ej(ωt+φ)}，若忽略ejωt，则正弦信号用一个复数X表示，称之为信号的相量表示形式；
初始计算的算法令fs为采样率，对于额定频率为50Hz的正弦信号，一周期能够采样n个点；xk为k时刻的待测信号的采样值；经过一段时间的采样，得到t＝-(k-1)，t＝-(k-2)，……，t＝-2，t＝-1，t＝0时刻n个采样值；采样数据以堆栈的方式存储于缓存数组W1中；其中，W1[0]存储了最新时刻的采样值，W1[1]存储了上一时刻的采样值，……以此类推；所有缓存数组皆以这种方式对某次计算得到的数值进行存储；采用DFT算法，通过计算得到待测信号的幅值Xm与相位φ；计算公式为：
Figure FDA0003237072420000021
上式中：xk表示在k时刻对待测信号进行采样，得到的信号瞬时值；a表示经过DFT计算得到的相量实部；b表示经过DFT计算得到的相量虚部；Xm表示经过DFT计算得到的信号峰值；φ'表示相量实部与虚部相比后得到的反正切角度值；
反正切函数arctan的值域范围为[-π/2,π/2]，而相量的角度范围为(-π,π]，故需要判断(a,b)所在象限，经过转换后得到待测信号的相角φ，相应计算公式为：
Figure FDA0003237072420000022
上式中：a表示经过DFT计算得到的相量实部；b表示经过DFT计算得到的相量虚部；φ为经过象限转换后得到的待测信号的相角；
相位与幅值的计算结果以堆栈方式存储于一定长度的缓存数组P1与X1中。
3.根据权利要求1所述的一种基于DFT、多阶滤波和突变判据的同步相量计算方法，其特征在于：所述步骤2.频率计算，包括：
设k时刻、k+1时刻通过初始计算得到的相位φk、φk+1分别存储于缓存P[1]和P[0]中，待测信号的频率f1的计算公式为：
Figure FDA0003237072420000023
上式中：fs为对待测信号进行模数转换的采样频率；P为存储数据，括号内数字为数组下标；f1为用于频率计算的中间参数；
实际电力系统中，信号的频率会存在不同程度的频率偏移；当待测信号频率偏移额定频率时，DFT计算得到幅值和相位会含有近似的二次谐波；基于相位信息计算得到的频率f1，同样也会耦合二次谐波；除了计算误差，电网本身也存在着谐波，此外，信号还会受到各类噪声的影响；故计算得到的f1误差较大，不能直接作为待测信号的频率值；
为了提升频率计算的精度，需要进行滤波处理；滤波算法采用多阶结构，由不同的滤波算法组合而成；滤波算法包括三点滤波、定点数平均滤波、不定点数的平均滤波；
对比三种滤波算法的优缺点，三点滤波算法计算量小、延时小，但是滤波效果一般；递归式定点数平均滤波算法计算量小，延时相对较长，滤波效果较好；递归式变点数平均滤波算法计算量小，延时相对较长，滤波效果最好，但是需要以获知待测信号大致频率为计算的前提条件；多个滤波器级联，能够构成多阶滤波算法；滤波效果若要越好，则滤波阶数也越高，同时滤波输出对于输入信号的变化的响应也就越慢；滤波算法的设计要兼顾精度与响应速度；
采用该算法，对THD谐波含量为0.2％、频率为46Hz的稳态正弦信号进行频率测量,经过DFT计算得到相位后再直接计算得到f1由于待测信号频率为46Hz，偏移了50Hz的额定频率，故f1耦合了二次谐波，在40Hz与50Hz之间进行波动，误差较大；经过多阶滤波后，依次得到信号频率f3、f4、f5、f6，频率波动的幅度也依次减小；f3～f6对应的最大误差分别为400mHz、40mHz、5mH、0.06mHz；经过多阶滤波后得到的f6误差较其小，作为频率测量结果fc。
4.根据权利要求2所述的一种基于DFT、多阶滤波和突变判据的同步相量计算方法，其特征在于：所述三点滤波算法的原理为，对相位相差为0°、60°、120°的三个量进行相加，能够在一定程度上滤除二次脉动谐波；对于待测信号额定频率为50Hz，采样率fs为1kHz的情况，若数组X连续存储了t＝0时刻、t＝-1时刻、t＝-2时刻、……、t＝-n时刻的计算结果，则相差0°、60°、120°的三个数值依次为：X(0)、X(3)、X(7)；对其求平均即能够滤除二次谐波，对应的计算公式为：
Figure FDA0003237072420000031
上式中：X为存储数据，括号内数字为数组下标；out为滤波计算后得到的输出结果。
5.根据权利要求2所述的一种基于DFT、多阶滤波和突变判据的同步相量计算方法，其特征在于：所述定点数平均滤波原理为，对固定长度的数据进行平均化处理；对于长度为n的数组一维X进行平均滤波的计算公式为，
Figure FDA0003237072420000041
上式中：X为存储数据，括号内数字为数组下标；n为数组X的长度；out为滤波计算后得到的输出结果；
采用递归式的定点数平均滤波算法，有效减小滤波算法计算量；令数组sum为长度为n的数组X各成员数据之和，最新数据为xnew，最旧数据为xold，则平均滤波结果avg为，
Figure FDA0003237072420000042
上式中：sum为X数组所有数据之和；n为数组X的长度；avg为数组X的所有数据之和的平均值；
每次计算avg完成后，新数据xnew进栈，旧数据xold出栈剔除，以实现数据的不断刷新存储。
6.根据权利要求2所述的一种基于DFT、多阶滤波和突变判据的同步相量计算方法，其特征在于：所述不定点数的平均滤波，是指：当待测信号频率偏移额定频率时，频率计算结果耦合的谐波实际上只是近似于二次谐波，故采用固定点数的平均滤波效果仍然有限；采用频率跟踪的方式，进行变点数的平均滤波具有更好的滤波效果；每次滤波求取平均的点数n1的计算公式为，
Figure FDA0003237072420000043
上式中：fs为对待测信号进行模数转换的采样频率；f为通过初始计算得到的待测信号的频率，ceil函数为求不小于给定数的最小整数，n1为一采样周期内采样数据个数；
其中，ceil(x)函数的作用为，求不小于给定实数x的最小整数；fs为采样率，f为通过初始计算得到的待测信号的频率；与平均滤波算法类似，变点数滤波亦为递归形式。
7.根据权利要求1所述的一种基于DFT、多阶滤波和突变判据的同步相量计算方法，其特征在于：所述步骤3.幅值与相位的计算，包括：
稳态情况，采用频率fc对同步采样数据进行采样，能够得到待测信号的幅值xc和相位φc；暂态情况，采用f2对同步数据进行再采样，得到相位p3作为相位的计算相位φc；
当待测信号偏移了额定频率时，还按照额定频率的DFT进行幅值和相位计算，是造成计算产生较大误差的原因；此时需要采用再采样算法，等价于DFT计算频率与待测信号频率相同，因而能够得到精度较高的幅值与相位计算结果；再采样算法原理具体如下:
设偏移了额定频率的输入信号为:
x(t)＝Xmcos(ωt+φ)
上式中：x(t)为待测信号t时刻瞬时值；Xm为待测信号峰值；ω为待测信号角频率；φ为待测信号初始相角；
按照额定频率的整N倍数频率Nω0进行采样，第k个采样点和k+1个采样值为：
Figure FDA0003237072420000051
上式中：xk为离散化后，k时刻对待测信号进行采样得到的信号瞬时值；xk+1为离散化后k+1时刻对待测信号进行采样得到的信号瞬时值；Xm为待测信号峰值；θ为k+1时刻与k时刻相比，待测信号的相位变化大小；φ为待测信号初始相角；
若按照与待测信号实际频率ω对应的整N倍数频率Nω进行再采样，则第n个点能够表示为：
xn＝Xmcos(kθ+γ+φ)
上式中：xn为采用与待测信号频率整倍数采样率，对待测信号进行采样得到的瞬时值；Xm为待测信号峰值；k为对时间进行离散化后的时刻值，φ为待测信号初始相角；γ为再采样后的信号与最近时刻采样信号之间的相位差大小；
若xn对应的采样时间介于xk与xk+1之间，推导xk、xk+1、xn之间的关系有：
Figure FDA0003237072420000052
上式中：xk为离散化后k时刻对待测信号进行采样得到的信号瞬时值，xk+1为离散化后k+1时刻对待测信号进行采样得到的信号瞬时值；θ为k+1时刻与k时刻间待测信号的相位变化；γ为再采样后的频率偏移量大小；
其中，f为再采样前经过额定频率DFT计算得到的待测信号频率，fs为采样频率，Δt为xn与xk采样点之间的时间差；θ＝2πf/fs，γ＝2πfΔt；经过再采样后得到新的采样序列，对采样序列再进行一次DFT计算，得到待测信号幅值x4、相位p2；x4存入缓存X4，对X4进行一次FIR低通滤波，得到的幅值计算结果Xc具有较高的精度；
当待测信号稳定时，信号突变计数C为0，此时采用p2作为相位计算结果；而当波形发生幅值突变或相位突变时，需要快速计算波形相位；此时采用响应较快的频率f2作为再采样算法的采样频率，随之计算得到的相位p3，为相位计算结果。
8.根据权利要求1所述的一种基于DFT、多阶滤波和突变判据的同步相量计算方法，其特征在于：所述步骤4.频率变化率的计算，包括：
根据额定频率、频率偏移量、频率变化率和相位之间的关系，采用最小二乘法计算频率变化率f′；
其原理为，任意时刻t的角频率ω表示为：
ω(t)＝ω0+Δω+tω′
上式中：ω(t)为t时刻待测信号的角频率，ω0为额定角频率，Δω为偏移额定角频率的大小，ω′为角频率变化率；
相角φ是频率相对于时间的积分，有：
Figure FDA0003237072420000061
上式中：φ(t)为待测信号t时刻相位瞬时值，φ0为初始相位大小，ω0为额定角频率，Δω为频率偏移量的大小，ω′为角频率变化率；
可知φ(t)是关于t的二阶多项式函数；通过运算得到信号的相位角，将计算得到的n个相位角组成一维列向量：
[Φ]＝[φ0 φ1 φ2…φn-1]T
上式中：φx为第x次计算得到的待测信号相位大小；Φ为列向量，用于存储n个相位计算结果；
根据相位与频率的关系，列出方程：
Figure FDA0003237072420000071
上式中φx为第x次计算得到的待测信号相位大小，Φ为列向量，用于存储n个相位计算结果，Δt为两次采样间的时间间隔，ω0为额定角频率，Δω为频率偏移量的大小，ω’为角频率变化率；将上式简化为矩阵形式有[Φ]＝[B][A]；其中B为3列n行的矩阵，与采样间隔时间有关；A为1列3行的矩阵，称之为角频率特征参数矩阵；采用矩阵A、B用于简化表达式；
根据最小二乘法，计算得到矩阵A为：
Figure FDA0003237072420000072
上式中：B为简化后的矩阵，A为角频率特征参数矩阵，ω0为额定角频率，Δω为偏移额定角频率的大小，ω′为角频率变化率，Φ为列向量，用于存储n个相位计算结果；求解A矩阵得到ω0、Δω、ω′，进一步计算频率变化率f′，有：
f′＝ω′/2π
上式中：f′为计算得到的待测信号频率变化率，ω′为待测信号角频率变化率；
当采样率恒定时，B矩阵就固定了，因此每次计算前不用再计算B，而是在初始化前程序时直接令：
G＝[BTB]-1B
上式中：G为n阶方矩，B为3列n行的矩阵，与采样间隔时间有关；
由于相位角φ的值域为[-π，π)，在计算之前，[Φ]需要进行单调化处理，以保证φk<φk+1(k＝0，1，2……，n-2)。
9.根据权利要求1所述的一种基于DFT、多阶滤波和突变判据的同步相量计算方法，其特征在于：所述步骤5.信号突变判断算法，包括：
经过一轮计算后，得到信号的幅值x、相位φ、频率f、频率变化率r；由于在计算过程中采用了滤波算法，因此计算得到的结果相比原信号波形具有一定的延时；
为了能够区分稳态波形和突变波形两种状态、采取不同的数据处理方式以最小化计算误差，需要设计波形突变判据；根据DFT的计算原理，当信号发生突变时，若计算的时间窗了包含了突变数据，则DFT的计算结果也会发生相应的突变；故能够根据信号幅值的变化大小，作为信号发生突变的判据J，一旦J大于了某个阈值J，即认为待测信号发生了突变；突变判据J的计算公式为，
Figure FDA0003237072420000081
上式中：J为信号突变判据，X1为信号幅值计算存储数组，括号内数字为数组下标，abs为求绝对值函数；
公式中k的取值与采样率有关，一般要求X1(k)与X1(0)两者所对应的时间差要在8ms左右，这样计算得到的J大小合适；若k太大，则计算得到的信号突变时间点与实际信号突变时间点误差较大；若k太小，则J也较小，有可能无法判断信号突变是否发生；为了确定阈值Jmax的大小，对各种情况下的J进行了测试，只有当发生幅值突变和相位突变的情况下，J大于1；故取阈值Jmax＝1.0，当J>Jmax时，则认为波形发生了突变；
当发生突变后，通过判断J的大小能够知道何时波形会再趋于稳定；但是由于算法结构中包含滤波器，故当突变后J小于一定数值时，计算结果可能仍未稳定；因此，需要采用突变计数C来记录波形突变的长度；波形未突变时，C＝0；当波形突变后，C的大小一般等于一周期采样点数n的7.5～7.8倍，令这个值为Ck；当C≠0时，每次完成采样计算后，C需要减1，当C再次等于0时，且J<1，意味着信号再次稳定；否则，C＝Ck，待测波形仍为不稳定的变化状态。
10.根据权利要求1所述的一种基于DFT、多阶滤波和突变判据的同步相量计算方法，其特征在于：所述步骤6.算法测试，包括：频率偏移、谐波、幅值突变、相位突变、频率斜坡变化、幅值调制、相位调制及幅值相位同时调制；测试考察幅值、相位、频率、频率变化率、正序幅值以及正序相位。
CN201910517446.4A 2019-06-14 2019-06-14 一种基于dft、多阶滤波和突变判据的同步相量计算方法 Active CN110412415B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910517446.4A CN110412415B (zh) 2019-06-14 2019-06-14 一种基于dft、多阶滤波和突变判据的同步相量计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910517446.4A CN110412415B (zh) 2019-06-14 2019-06-14 一种基于dft、多阶滤波和突变判据的同步相量计算方法

Publications (2)

Publication Number Publication Date
CN110412415A CN110412415A (zh) 2019-11-05
CN110412415B true CN110412415B (zh) 2021-11-23

Family

ID=68359144

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910517446.4A Active CN110412415B (zh) 2019-06-14 2019-06-14 一种基于dft、多阶滤波和突变判据的同步相量计算方法

Country Status (1)

Country Link
CN (1) CN110412415B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111381285B (zh) * 2020-03-31 2022-06-10 唐山学院 多频电流过零方波信号处理方法及装置
CN113535812B (zh) * 2021-06-29 2024-01-30 浙江中控技术股份有限公司 工况稳态检测方法、过程优化方法
CN114665992B (zh) * 2022-03-22 2023-10-27 吉林省广播电视研究所(吉林省广播电视局科技信息中心) 基于平均极值测量信号指标的方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104635094A (zh) * 2015-03-02 2015-05-20 国电南瑞科技股份有限公司 一种提升pmu同步相量测量精度的方法
CN105116453A (zh) * 2015-08-14 2015-12-02 中国石油天然气股份有限公司 一种冻土带天然气水合物的瞬变电磁勘探方法及装置
CN107345984A (zh) * 2017-06-23 2017-11-14 华北电力大学 一种基于信号识别的自适应同步相量测量方法
CN108833058A (zh) * 2018-05-25 2018-11-16 国网上海市电力公司 一种广域测量系统通信过程动态数据压缩、解压的方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10101392B2 (en) * 2014-07-11 2018-10-16 Advantest Corporation Scan test multiplexing

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104635094A (zh) * 2015-03-02 2015-05-20 国电南瑞科技股份有限公司 一种提升pmu同步相量测量精度的方法
CN105116453A (zh) * 2015-08-14 2015-12-02 中国石油天然气股份有限公司 一种冻土带天然气水合物的瞬变电磁勘探方法及装置
CN107345984A (zh) * 2017-06-23 2017-11-14 华北电力大学 一种基于信号识别的自适应同步相量测量方法
CN108833058A (zh) * 2018-05-25 2018-11-16 国网上海市电力公司 一种广域测量系统通信过程动态数据压缩、解压的方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
PMU同步采样数据压缩算法研究;胡丽丽;《中国优秀硕士学位论文全文数据库 信息科技辑》;20160515;I140-259 *
基于同步相量的高密度量测数据压缩技术研究;蔡玉朋 等;《电测与仪表》;20200302;第57卷(第13期);第91-97页 *
瞬变电磁法数据滤波方法对比研究;张建锋 等;《工程勘察》;20150801;第43卷(第8期);第93-98页 *

Also Published As

Publication number Publication date
CN110412415A (zh) 2019-11-05

Similar Documents

Publication Publication Date Title
CN110412415B (zh) 一种基于dft、多阶滤波和突变判据的同步相量计算方法
Wang et al. A novel phase-locked loop based on frequency detector and initial phase angle detector
US8067932B2 (en) Advanced real-time grid monitoring system and method
US5721689A (en) System and method for phasor estimation and frequency tracking in digital protection systems
CN108896820B (zh) 一种适用于静止变频器启动的调相机启机保护相量计算方法
CN109412191B (zh) 一种用于高压直流输电系统的锁相方法、装置和设备
CN108614155B (zh) 一种加入汉明窗的同步相量测量方法及系统
CN103983847B (zh) 一种同步相量测量中基于rls的自适应频率跟踪测量方法
CN110535161A (zh) Lcl型储能变换器的有限控制集模型预测控制方法
CN107144734B (zh) 一种适用于pmu的配电网高精度相量测量方法
CN111817713B (zh) 可快速同步对称故障下电压相位的高压直流锁相环及方法
US5832413A (en) Generator protection system and method for phasor estimation and frequency tracking during frequency ramping
CN110954763A (zh) 一种基于谐波电流注入和谐波阻抗测量的微电网非破坏性孤岛检测方法
CN113805010A (zh) 一种配电网单相接地故障的研判方法及系统
CN107315103A (zh) 一种电力冲击负载检测方法
Ukil et al. Power systems frequency estimation using amplitude tracking square wave for low-end protective relays
JP5537095B2 (ja) ベクトル計測装置及びベクトル計測システム
CN108717141B (zh) 利用单相电压测量电气量频率的方法、测量系统
Destro et al. Implementation aspects of adaptive window moving average filter applied to PLLs—Comparative study
CN113013854B (zh) 基于动态过程中差动电流变化率的区内外故障预判方法
CN113258555B (zh) 一种柔性直流输电系统的谐波谐振抑制方法及系统
US20230314492A1 (en) Providing a frequency of an electrical quantity in an electrical power system
CN109298286A (zh) 一种判断电能质量暂降原因及设计暂降源定向算法的方法
WO2021082036A1 (zh) 一种电力系统频率测量方法、母线电压校正方法及装置
CN113672863A (zh) 一种发电机启机保护的简化相量计算方法及系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
EE01 Entry into force of recordation of patent licensing contract
EE01 Entry into force of recordation of patent licensing contract

Application publication date: 20191105

Assignee: ZHONGYU HIGH-TECH (WUHAN) INFORMATION TECHNOLOGY Co.,Ltd.

Assignor: STATE GRID LIAONING ELECTRIC POWER Research Institute

Contract record no.: X2023980034642

Denomination of invention: A Synchronous Phasor Calculation Method Based on DFT, Multi order Filtering, and Catastrophe Criterion

Granted publication date: 20211123

License type: Common License

Record date: 20230411

Application publication date: 20191105

Assignee: Beijing ANGLI Sifang Electrical Equipment Co.,Ltd.

Assignor: STATE GRID LIAONING ELECTRIC POWER Research Institute

Contract record no.: X2023980034637

Denomination of invention: A Synchronous Phasor Calculation Method Based on DFT, Multi order Filtering, and Catastrophe Criterion

Granted publication date: 20211123

License type: Common License

Record date: 20230411

Application publication date: 20191105

Assignee: NANJING NENGYUN ELECTRIC POWER TECHNOLOGY CO.,LTD.

Assignor: STATE GRID LIAONING ELECTRIC POWER Research Institute

Contract record no.: X2023980034640

Denomination of invention: A Synchronous Phasor Calculation Method Based on DFT, Multi order Filtering, and Catastrophe Criterion

Granted publication date: 20211123

License type: Common License

Record date: 20230411

Application publication date: 20191105

Assignee: Yilixin (Guangzhou) Electronic Technology Co.,Ltd.

Assignor: STATE GRID LIAONING ELECTRIC POWER Research Institute

Contract record no.: X2023980034641

Denomination of invention: A Synchronous Phasor Calculation Method Based on DFT, Multi order Filtering, and Catastrophe Criterion

Granted publication date: 20211123

License type: Common License

Record date: 20230411