本发明涉及电流互感器噪声滤波领域，尤其是一种基于时间序列的光纤电流互感器随机噪声实时滤波方法。
背景技术：
电子式互感器是智能变电站中一种重要的过程层设备，为智能变电站间隔层设备提供电流、电压信息，实现电力系统继电保护、电能计量、故障录波、状态监测等功能。光纤电流互感器实现了电流信号的实时测量和监测，为电力系统继电保护、电能计量等设备提供输入信号。
光纤电流互感器技术已经日趋成熟，但仍存在输出噪声和随机误差等参数值偏高的问题，严重影响了光纤电流互感器的工程应用。为降低光纤电流互感器输出信号的噪声和随机误差，可采用滤波的方法。有效滤波的前提是能够精确建立光纤电流互感器的输出模型，尤其对于卡尔曼滤波，光纤电流互感器的输出模型精确程度直接关系着滤波效果的好坏。目前，在针对光纤电流互感器输出信号的模型建立，可采用固定的模型，如AR(2)、ARMA(2，1)，然后加以滤波，但不能忽略光纤电流互感器输出序列的模型统计分析以及数据分析过程输出序列的非平稳性、非随机和非正态等特性，对输出序列进行独立、平稳、正态、零均值以及趋势项处理，而直接采用ARMA时间序列模型分析从基础上不符合时间序列适用平稳序列这一特性，因此直接建立的模型误差较大。卡尔曼滤波是一种利用线性系统状态方程，通过系统输入输出观测数据，对系统状态进行最优估计的算法。在测量方差已知的情况下能够从一系列存在测量噪声的数据中，估计动态系统的状态，在求出模型参数和噪声统计特性的基础上，经过卡尔曼滤波器，可达到良好的消噪效果，现已在通信、导航、制导与控制等多领域得到了较好的应用。
技术实现要素：
本发明所要解决的技术问题在于，提供一种基于时间序列的光纤电流互感器随机噪声实时滤波方法，可克服现有光纤电流互感器随机噪声滤波方法的不足。
为解决上述技术问题，本发明提供一种基于时间序列的光纤电流互感器随机噪声实时滤波方法，包括如下步骤：
(1)采集光纤电流互感器输出的某相电流数据{x0(t)}；
(2)电流序列平均滤波处理；对有效序列进行平均滤波处理，计算每十个数据的均值作为新的序列{x1(t)}；
(3)建立时间序列模型；
(a)平稳性检验；采用单位根检验法则，不满足平稳性要求则进行差分处理，以获取平稳的电流序列{xn}；
(b)正态性检验；采用偏峰态检验法，电流序列{xn}具有以下4个表示其总体概率密度函数的参数如下：
当g1和g2约等于0时，可以认为电流序列{xn}为正态时序；
(c)周期性检验；周期性检验用来识别光纤电流互感器输出数据中是否包含有随机量以外的周期性分量，这样在解释数据分析的结果时可以避免出现错误；周期性检验的方法是直接考察从输出数据中得到的概率密度函数或自相关函数或功率谱密度函数的图形；
(d)建立时间序列模型；依次选取p和q值计算各模型的AIC值，然后选择最小的AIC值，确定时间序列模型的阶次，即p和q值，采用最小二乘法拟合出模型参数，写出时间序列数学模型；
(4)写出与确定的时间序列模型相对应的卡尔曼滤波方程的状态空间模型；
状态方程：Xk＝AXk+BVk；
输出方程：Yk＝CXk+Wk；
其中，Vk和Wk的统计特性为：
系统的状态方程为过程噪声为Vk＝[rk，0]T；
对于AR(p)模型，
Vk＝[rk，0]T；
对于ARMA(p，q)模型，
Vk＝[rk，...，rk-q]T；
卡尔曼滤波处理，采用以下卡尔曼滤波器对光纤电流互感器的时间序列输出信号进行滤波处理
Kk＝Pk，k-1CT[CPk，k-1CT+Rk]-1
pk，k-1＝APk，k-1AT+BQk-1，kBT
Pk，k＝[I-KkC]Pk，k-1
式中，为滤波状态的进一步估计，为k时刻滤波器的状态，Kk为k时刻滤波器的增益矩阵，R为系统量测噪声误差，Q为系统过程噪声方差，Pk，k为滤波器误差协方差矩阵，为k时刻滤波器的输出。
优选的，步骤(2)中采用平均值滤波的方法对原始数据进行处理。
优选的，步骤(3)中采用AIC信息准则来确定时间序列模型的阶次，计算各个模型的AIC值，取AIC值最小的模型为使用模型；采用快速算法RLS进行AR模型的参数估计，采用长自归白噪声估计法建立ARMA模型；采用检验模型残差是否为白噪声来确定模型的可用性。
优选的，步骤(3)中建立时间序列模型包括以下步骤：
(a)确定模型阶次；光纤电流互感器的建模采用AR模型和ARMA模型：AR模型指任一时刻k上的数值yk表示过去p个时刻上数值的线性组合加上k时刻的白噪声，表示为：
yk＝a1yk-1+…+apyk-p+rk
式中，常整数p为模型的阶次，常系数a1,…,ap为模型参数，{rk}为均值为0、方差为σ2的白噪声，p阶模型记为AR(p)；
ARMA模型是在AR模型的基础上减去过去q个时刻上白噪声的线性组合，表示为：
yk＝a1yk-1+…+apyk-p+rk-θ1rk-1-θ2rk-2-…-θqrk-q
式中，常整数p、q为模型的阶次，常系数a1,…,ap,θ1,…θq为模型参数，模型记为ARMA(p,q)；
一般首先对数据序列进行AR(p)建模，若找不到适用模型，再进行ARMA(p,q)建模；估计模型参数之前，需要确定模型的阶次；采用AIC信息准则来确定模型阶次；AIC信息准则的简化公式为：
式中，p、q为ARMA模型阶次；n为序列中数据个数；为噪声a(t)的方差；建立AR或ARMA时间序列模型，计算各个模型的AIC值，取AIC值最小的模型为使用模型，确定模型阶次；
(b)AR模型参数估计；采用快速算法RLS进行AR模型的参数估计；基于已知观测值(yk,yk-1,…,y0,…,y1-p)求估计值和式可写成如下向量形式：
式中，T表示矩阵的转置，αT＝a1,…,ap；
定义AR(p)模型参数α的估计公式如下：
初值和是利用少量观测数据(y1,…,y20)通过以下两式求得：
其中，
(c)ARMA模型参数估计；采用长自回归白噪声估计法建立ARMA模型，主要步骤为：
①建立长自回归模型AR(pN)；阶次取lgN的适当倍数，即pN＝(lgN)1+δ，选δ为一个0到1之间的正数；AR(pN)的自回归参数为
由线性最小二乘估计得到
其中，YpN＝(xpN+1,...,xN)T
②求长回归模型残差，检验其独立性；用步骤①所得即样本值x＝(x1,…,xN)T，计算残差检验的独立性；若不独立，则增大pN，再重新进行①、②两步，否则进行下一步；
③估计ARMA模型参数，联合白噪声估计值t＝pN+1,...,N，N和样本值xt，t＝1,...,N，按线性最小二乘估计ARMA(p,q)模型的各参数值：
其中，Ypq＝(xpN+1,xpN+2,...,xN)T，
上式中在步骤②由长自回归估计得到，所以该算法不涉及非线性求解问题，只用到线性最小二乘估计；
(d)模型适用性检验；检验残差是否为白噪声，如果模型残差为白噪声，则认为模型可用；反之，不可用；检验是否为白噪声通过分析自相关函数图和偏自相关函数图来确定。
优选的，步骤(4)中采用卡尔曼滤波的方法对建立的时间序列模型进行滤波处理。
本发明的有益效果为：通过平均值滤波方法的处理，提高了数据的有效性；经过平稳性、周期性和正态性检验，保证了建模的可行性；AR模型和ARMA模型的建立、阶次选择及模型检验，确保了数据模型的普适性；卡尔曼滤波算法处理光纤电流互感器的输出数据，有效的滤除随机噪声，提高了测量的精确性。
附图说明
图1是本发明的方法流程图。
具体实施方式
如图1所示，一种基于时间序列的光纤电流互感器随机噪声实时滤波方法，包括如下步骤：
步骤一：采集光纤电流互感器输出的某相电流数据{x0(t)}；
步骤二：电流序列平均滤波处理
对有效序列进行平均滤波处理，计算每十个数据的均值作为新的序列{x1(t)}；
步骤三：建立时间序列模型
(1)平稳性检验
采用单位根检验法则，不满足平稳性要求则进行差分处理，以获取平稳的电流序列{xn}；
(2)正态性检验
采用偏峰态检验法，电流序列{xn}具有以下4个表示其总体概率密度函数的参数如下，
均值：
方差：
标准偏度系数：
标准峰度系数：
当g1和g2约等于0时，可以认为电流序列{xn}为正态时序。
(3)周期性检验
周期性检验用来识别光纤电流互感器输出数据中是否包含有随机量以外的周期性分量，这样在解释数据分析的结果时可以避免出现错误。周期性检验的方法是直接考察从输出数据中得到的概率密度函数或自相关函数或功率谱密度函数的图形。
(4)建立时间序列模型，依次选取p和q值计算各模型的AIC值，然后选择最小的AIC值，确定时间序列模型的阶次，即p和q值，采用最小二乘法拟合出模型参数，写出时间序列数学模型；
a)确定模型阶次
光纤电流互感器的建模采用AR模型和ARMA模型：AR模型指任一时刻k上的数值yk表示过去p个时刻上数值的线性组合加上k时刻的白噪声，表示为：
yk＝a1yk-1+…+apyk-p+rk
式中，常整数p为模型的阶次，常系数a1，…，ap为模型参数，{rk}为均值为0、方差为σ2的白噪声，p阶模型记为AR(p)；
ARMA模型是在AR模型的基础上减去过去q个时刻上白噪声的线性组合，表示为：
yk＝a1yk-1+…+apyk-p+rk-θ1rk-1-θ2rk-2-…-θqrk-q
式中，常整数p、q为模型的阶次，常系数a1，…，ap，θ1，…θq为模型参数，模型记为ARMA(p，q)；
一般首先对数据序列进行AR(p)建模，若找不到适用模型，再进行ARMA(p,q)建模；估计模型参数之前，需要确定模型的阶次；采用AIC信息准则来确定模型阶次；AIC信息准则的简化公式为：
式中，p、q为ARMA模型阶次；n为序列中数据个数；为噪声a(t)的方差。建立AR或ARMA时间序列模型，计算各个模型的AIC值，取AIC值最小的模型为使用模型，确定模型阶次；
b)AR模型参数估计
采用快速算法RLS进行AR模型的参数估计；基于已知观测值(yk,yk-1,…,y0,…,y1-p)求估计值和式可写成如下向量形式：
式中，T表示矩阵的转置，αT＝a1,…,ap；
定义AR(p)模型参数α的估计公式如下：
初值和P0是利用少量观测数据(y1,…,y20)通过以下两式求得：
其中，
c)ARMA模型参数估计
采用长自回归白噪声估计法建立ARMA模型，主要步骤为：
①建立长自回归模型AR(pN)；阶次取lgN的适当倍数，即pN＝(lgN)1+δ，选δ为一个0到1之间的正数；AR(pN)的自回归参数为
由线性最小二乘估计得到
其中，YpN＝(xpN+1,...,xN)T
②求长回归模型残差，检验其独立性；用步骤①所得即样本值x＝(x1,…,xN)T，计算残差检验的独立性；若不独立，则增大pN，再重新进行①、②两步，否则进行下一步；
③估计ARMA模型参数，联合白噪声估计值t＝pN+1,...,N，N和样本值xt，t＝1,...,N，按线性最小二乘估计ARMA(p,q)模型的各参数值：
其中，Ypq＝(xpN+1,xpN+2,...,xN)T，
上式中在步骤②由长自回归估计得到，所以该算法不涉及非线性求解问题，只用到线性最小二乘估计。
d)模型适用性检验
检验残差是否为白噪声，如果模型残差为白噪声，则认为模型可用；反之，不可用。检验是否为白噪声通过分析自相关函数图和偏自相关函数图来确定。
步骤四：写出与确定的时间序列模型相对应的卡尔曼滤波方程的状态空间模型
状态方程：Xk＝AXk+BVk
输出方程：Yk＝CXk+Wk
其中，Vk和Wk的统计特性为：
系统的状态方程为过程噪声为Vk＝[rk,0]T；
Vk＝[rk，0]T
对于ARMA(p,q)模型，
Vk＝[rk,…,rk-q]T；
卡尔曼滤波处理，采用以下卡尔曼滤波器对光纤电流互感器的时间序列输出信号进行滤波处理
Kk＝Pk,k-1CT[CPk,k-1CT+Rk]-1
Pk,k-1＝APk,k-1AT+BQk-1,kBT
Pk,k＝[I-KkC]Pk,k-1
式中，为滤波状态的进一步估计，为k时刻滤波器的状态，Kk为k时刻滤波器的增益矩阵，R为系统量测噪声误差，Q为系统过程噪声方差，Pk,k为滤波器误差协方差矩阵，为k时刻滤波器的输出。
尽管本发明就优选实施方式进行了示意和描述，但本领域的技术人员应当理解，只要不超出本发明的权利要求所限定的范围，可以对本发明进行各种变化和修改。