本发明属于电力系统动态频率测量技术，具体涉及一种最小二乘拟合动态频率测量方法。
背景技术：
频率是电能质量的重要参数之一，能够反映电力系统瞬时运行状态的变化。从电力供应网络的全局来评价，频率是电网有功功率平衡状态的指示量。测量频率有助于监测电力系统电能发、输、配、用的全过程。从电能质量的角度来分析，频率是衡量电能质量的重要标准之一。频率偏高或者偏低都会影响设备的运行效率，甚至会对设备造成不可逆转的损害。以电力系统安全稳定运行为出发点和归宿点，频率是电力系统继电保护设备的重要控制参数之一。当电力系统处于非正常运行状态时，频率监测有助于辨识电力系统运行状态，并为继电保护装置、低频减载装置的动作提供判断依据。总的来说，实时、精确地频率测量有助于电网安全、稳定、经济运行。
目前，国内外大量学者已经在数字化测量技术的基础上提出了很多软件测频理论和算法。与硬件测量法相比，软件测频法不依赖于硬件电路，研究的重点主要在数学方法上。因而软件测频法灵活多变，适用场景广泛，能够匹配多样化的电网运行状态。一般而言，频率测量包括三个步骤：信号预处理、频率测量、结果再处理。其中，信号预处理主要是滤波环节和筛选数据，可以利用硬件电路或者滤波辅助算法达到预处理的目标。结果再处理同样为辅助算法，能够检验异常值或者通过辅助算法优化测量结果，使测量结果达到工程应用的要求。频率测量算法是频率测量步骤的核心部分。经典的测频算法主要有过零点检测法、函数解析法、函数逼近法、离散傅里叶变换算法及其改进算法、小波分析法等等。
上述算法都是在基本观测模型的基础上，不断改进电网信号的数学模型，目的在于逼近真实的物理信号。同时，上述算法在适用范围、测量精度、响应时间、算法复杂度、谐波抑制和噪声抗干扰等方面各有特点。但是，上述算法只对稳定频率测量有较好的效果，而对动态频率的测量缺乏支撑。
技术实现要素：
发明目的：本发明的目的在于解决现有技术中存在的不足，提供一种最小二乘拟合动态频率测量方法，改进传统的最小二乘测频法，不仅可用于稳定频率的测量，而且在频率动态变化特别是频率线性变化的状态下，能够实现频率的实时测量；还能直接测量出频率变化率，增强对频率变化速率的监测。
技术方案：本发明所述的一种最小二乘拟合动态频率测量方法，依次包括以下步骤：
一种最小二乘拟合动态频率测量方法，其特征在于：依次包括以下步骤：
(1)设置移动数据窗的时间长度t0和采样时间间隔△t，将单个移动数据窗的电网连续信号通过a/d转换采样，并将第一个采样点的时间参数置零，从而形成该数据窗内的等时间间隔离散信号；
(2)利用fir数字滤波器对待测离散信号进行工频基波提取；
(3)从工频基波提取后的离散信号中提取频率偏差和频率变化量偏差，对该数据窗内的采样点进行二元函数泰勒展开，构建量测状态模型；该模型的矩阵方程为：
[a][y]＝[x]
[a]为由采样时间间隔△t、采样点序数n和预估频率f0组成的常系数矩阵，[y]为含有待求频率参数组成的10×1的矩阵，[x]为n个采样点组成的n×1的矩阵，
(4)利用预估的频率值，离线算出常系数矩阵[a]；
(5)应用最小二乘法，按最小平方误差原理进行电网信号曲线拟合，求解量测状态矩阵方程，得到[y]；
(6)应用公式
应用公式
其中，y1，y2，y3，y4，y5，y6依次分别为[y]中的第1至6列；
(7)该数据窗的最后一个采样点的动态频率可以根据该点的采样时间、频率偏差量和频率变化率偏差量求得，即公式f＝f0+△f+n△k△t；f0为预估频率，f为实际频率；
(8)当n<m，时间向前推进一个采样时间间隔，数据窗也随之剔除第1个采样点和加入新的采样点，返回步骤(3)，形成新的量测状态矩阵，进行最新采样点动态频率测量流程；n是指采样点序数，m是指采样点的数量；
(9)当n≥m，结束测量流程。
进一步的，实现实时性与精确度的统一，所述步骤(1)中的移动数据窗的时间长度t0取值范围是0.01-0.04s，t0不宜过长，会加大算法的复杂度。
进一步的，所述步骤(2)中的fir数字滤波器采用通带频率范围在40～60hz的8阶fir数字滤波器，该滤波器的差分方程表达式为
z(n)＝0.02712x(n)+0.09165x(n-1)+0.17275x(n-2)+
0.23402x(n-3)+0.23402x(n-4)+0.17275x(n-5)
+0.09165x(n-6)+0.02712x(n-7)
式中，x(n)为第n个采样点的数值。
进一步的，所述步骤(3)中，基波状态下的量测状态模型为：
(3.1)在电力系统出现大量有功缺额的情况下，电网的频率变化是一个非常复杂的动态过程，若将电力系统视为等值的单机系统，频率特性可以用数学公式表达为：
tf为系统频率变化过程中的时间常数，一般在4s～6s间变化；△f为频率偏差，满足公式△f＝f-f0，f为实际频率，f0为基波频率；△p为该单机系统的有功缺额；kl为负荷的频率调节系数；又因为
在非常短的时间内，将上式进行线性化处理，令频率变化率为k，则有f＝f'+kt，从而短时间内的频率动态变化问题可以转化为频率线性变化问题；电网信号函数可以表示为
x(t)为电网波形的单相电压或者电流，a1为基波峰值，f‘为待测量的基波频率，k为频率变化率，
(3.2)电网信号函数是关于频率变化率和频率的二元函数，由于频率关于时间的函数为正弦函数的隐函数，自变量f‘和k很难直接通过正弦函数解出，然后通过该二元函数的泰勒展开公式提取自变量△f和△k，从而便于建立自变量f‘和k的线性方程，单个采样点的线性方程可以表示为：
n为采样点的序号，△t为采样时间间隔，将数据窗内所有采样点的方程联立起来，分离未知量和常系数量，就构成了量测状态矩阵方程[a][y]＝[x]，[x]为单个数据窗所有采样点数值组成的矩阵；
常系数矩阵[a]中的第n行可以表示为
未知矩阵[y]中的参数为
(3.3)常系数矩阵[a]是由采样时间间隔△t、采样点序数n和预估频率f0构成，第一个数据窗的初始预估频率采用三点测频法进行估计，频率公式为：
其中，x(k-1)、x(k)和x(k+1)为相邻三个采样点的数值，第二个数据窗的初始预估频率采用第一个数据窗所测量的频率值，以此类推，每个数据窗的常系数矩阵可以离线求出，从而减少测量方法的复杂度。
进一步的，所述步骤(5)中，利用最小二乘法进行曲线拟合，通过最小化误差的平方和寻找矩阵方程的最优解，即使用公式[y]＝{[a]t[a]}-1[a]t[x]确定含有频率和频率变化率的未知参数矩阵[y]。
有益效果：发明不仅可以用于稳定频率的测量，而且在频率动态变化，特别是频率线性变化的状态下，能够实现频率的实时测量。同时，该方法还能直接测量出频率变化率，增强对频率变化速率的监测。
附图说明
图1为本发明的流程示意图。
具体实施方式
下面对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。
本发明利用二元函数泰勒展开提取电网信号模型的频率和频率变化率参数，在获取等时间间隔的离散信号后建立量测矩阵方程。然后，通过最小二乘法解决线性方程组的拟合问题，求取矩阵方程的参数。为了避免采样数据的递增对计算量带来的影响，设定固定时间的数据窗进行参数计算，在保证精度的同时，能够提高算法的实时性和对突变信号的敏感性。
实施例1：
如图1所示，本实施例步骤具体如下：
1、设置数据窗时间t＝0.02s和采样频率fs＝1600hz，对电网信号进行采样。每个数据窗的采样点为24个，起始采样时间置零。
2、利用fir数字滤波器对数据窗的信号进行处理，得到较为纯净的基波信号并形成信号矩阵[x]。
3、采用二元函数泰勒级数从待测离散信号中提取频率偏差和频率变化量偏差，对该数据窗内的采样点进行二元函数泰勒展开。单个采样点展开公式为：
式中，n为采样点的序列号，
4、应用最小二乘法拟合量测状态矩阵方程，求解未知参数矩阵[y]。求解公式为:[y]＝{[a]t[a]}-1[a]t[x]。
5、由于未知参数矩阵[y]包含了频率偏差量和频率变化率偏差量，利用公式分别求取频率偏差量和频率变化率偏差量：
式中，y1，y2，y3，y4，y5，y6分别为[y]中的第1、2、3、4、5、6列。再根据式f＝f0+△f+0.02*△k就可以得到该数据窗内第24个采样点的频率。
6、当n<n，时间向前推进一个采样时间间隔，数据窗也随之剔除第n-23个采样点和加入第n+1个采样点，返回步骤(2)，进行该数据窗的频率测量。n为所需测量频率的最后一个采样点的序号。
7、当n≥n，结束流程。
至此，完成了电网信号动态频率的实时测量和电网信号频率变化率的跟踪。
从上述实施例可以看出，本发明不仅可以适用于稳定频率的测量，而且在频率动态变化，特别是频率线性变化的状态下，能够实现频率的实时测量；本发明以多元泰勒级数为基础，建立电网信号测量矩阵方程，并且采用最小二乘法解决方程线性拟合问题，在跟踪频率的同时，还能对频率变化率进行监测，具有精确度高，实时性好和抗干扰性强等特点。经验证，本发明的误差在10-3hz以内，测量速度略大于0.02s。此外，该算法还具有一定的抗干扰性和鲁棒性。