本发明属于配电网技术领域，具体公开了一种野值影响下基于事件触发机制的谐波检测方法，以解决目前针对野值和通信资源受限影响下谐波检测精度不高的技术问题。本发明方法包括如下步骤：首先，针对含有直流衰减分量的谐波信号采用了一种动态检测模型，以提高谐波检测的精度；其次，为了降低系统性能对野值的敏感性并达到节省网络通信资源的目的，设计了一种融合野值检测的事件触发器；在此基础上，将事件触发机制引入卡尔曼滤波器，并充分考虑非触发误差对滤波性能的影响，以最小化误差协方差矩阵的上界为目标，获得滤波器的增益。本发明方法对于直流衰减分量及谐波的参数具有较高的检测精度。
1.野值影响下基于事件触发机制的谐波检测方法，其特征在于，包括如下步骤：步骤1. 建立谐波动态检测状态空间模型和量测模型；
将通过传感器采样得到的电力谐波模拟信号转换为含有直流衰减分量的离散谐波信号；
根据含有直流衰减分量的离散谐波信号表达式，设置系统状态变量；给定噪声和野值的假设条件，并计算过程噪声协方差矩阵；建立谐波动态检测状态空间模型和量测模型；
步骤2. 建立融合野值检测的事件触发器，将传感器的量测值传输至数据处理中心；
具体为：通过构造触发函数以及野值检测函数，判断当前时刻量测值是否满足事件触发条件；若满足事件触发条件，则当前时刻量测值通过通信网络传输至数据处理中心；
在数据处理中心设置触发探测器以及滤波器；其中，触发探测器根据通信网络传输过来的数据判断是否发生了事件触发，滤波器用于量测值并估计谐波信号的幅值和相角；
若发生了事件触发，则触发探测器赋值为1，将包含有效信息的当前时刻量测值传输到滤波器；否则，触发探测器赋值为0，利用最新时刻传输值代替当前时刻量测值；
步骤3. 建立基于融合野值检测事件触发机制的卡尔曼滤波器；
建立卡尔曼滤波器结构；根据已建立的谐波动态检测状态空间模型，计算一步预测值和预测误差协方差矩阵，计算估计误差协方差矩阵及其上界；
通过最小化估计误差协方差矩阵及其上界的迹得到卡尔曼滤波器的增益参数；最后利用卡尔曼滤波器获得谐波状态的估计值，得到谐波信号的幅值和相角估计值。
2.根据权利要求1所述的野值影响下基于事件触发机制的谐波检测方法，其特征在于，所述步骤1具体为：步骤1.1. 首先给出含有直流衰减分量的谐波信号表达式，如公式(1)所示；
yk=∑i=1Aisin(iwkT+φi)+Ace (1)其中，yk表示理想情况下的谐波信号量测值，n表示谐波分量的最高阶数；
Ai和φi分别表示第i次谐波信号的幅值和相角，i=1,2,…,n，w为角频率，Ace 为直流衰减分量，Ac和a为常数，k和T分别表示采样时刻和采样周期；
步骤1.2. 设置系统状态变量xk，如公式(2)所示；
xk=[ x1,k x2,k… x2n‑1,k x2n,k xc,k] (2)‑akT
其中，x2i‑1,k = Aisin(iwkT+φi)，x2i,k = Aiiw cos(iwkT+φi)，xc,k =Ace ；
步骤1.3. 给定噪声和野值的假设条件，并计算过程噪声协方差矩阵；
设过程噪声ωk=[ω1,k ,ω2,k ,…,ωn,k,ωc,k] 和量测噪声vk；其中，ωi,k表示第i次谐波信号的过程噪声，ωc,k表示直流衰减信号的过程噪声；ωi,k以及ωc,k的表达式如下所示；
其中，ωi(τ)表示均值为0，且方差为σi =Ai /2π的正弦信号的高斯噪声；ωc(τ)表示均
值为0，且方差为σc =Ac /2π的直流衰减信号的高斯噪声；ωi,k、ωc,k和vk满足不等式(3)；
||ωi,k|| ≤εi，|ωc,k| ≤εc，|vk| ≤υ (3)其中，εi、εc和υ为给定的正标量，分别对应表示ωi,k、ωc,k和vk的上界；
ε是一个标量，表示||ωk||的上界；过程噪声ωk和测量噪声vk相关，统计特性如下：E{ωk }=0，E{vk}=0，E{xkωk }=0，E{xkvk }=0；
Qk=E{ωkωl }= Qkδk,l，Rk=E{vkvl }= Rkδk,l，Sk=E{ωkvl }= Skδk,l；
其中，xk表示k时刻的系统状态，ωl和vl分别表示l时刻的过程噪声和量测噪声，δk,l为克罗内克函数，Sk为相关噪声协方差矩阵，Rk为测量噪声协方差矩阵；
Qk=diag{Qk ,Qk ,…,Qk ,Qc.k}为过程噪声协方差矩阵；
其中，κ为具有可调范围的参数；定义mk为间歇发生的野值，描述为：；
式中，mj,k为第j个野值的幅值，t(j)表示第j个野值出现的时刻；
野值出现的时间间隔η和幅值mj,k满足η＞q，||mj,k||＞ζ；
其中，η=min{θj}j≥1为野值出现的最小时间间隔，θj=t(j+1)－t(j)，q和ζ为已知正标量；
步骤1.4. 考虑过程噪声、量测噪声以及野值的影响，建立谐波动态检测状态空间模型和量测模型，分别如公式(4)和公式(5)所示；
xk+1 = F xk+ ωk (4)yk+1= C xk+1+ vk+1+ mk+1 (5)其中，xk+1和yk+1分别表示k+1时刻的系统状态和量测值，vk +1和mk+1分别表示k+1时刻的量测噪声和野值；F=diag{F1,F2,…,Fn,Fc}为状态转移矩阵；
，Fc=e ；C=[1 0 1 0 … 1 0 1]∈R 为观测矩阵。
3.根据权利要求2所述的野值影响下基于事件触发机制的谐波检测方法，其特征在于，所述步骤2具体为：步骤2.1. 利用当前时刻量测yk、最新时刻传输值yks以及预设触发阈值γ，设计触发函数h(k)，判断当前时刻量测是否为“不必要量测”；触发函数h(k)如公式(6)所示；
h(k)=(yk－yks)(yk－yks) －γ(6)其中，ks表示第s个触发时刻，满足k0＜k1＜…＜ks＜…； 若h(k)<0，则当前时刻量测为“不必要量测”且不被传输；否则，进行下一步判断；
步骤2.2. 利用最近的q个时刻的量测yk‑q以及噪声和野值的幅值上界，设计野值检测函数f(k)以及野值检测阈值ρ，初步判断野值；
首先，利用传递函数表示公式(4)和(5)中模型的输入输出关系，如公式(7)所示；
y(z)=Gyω(z)ω(z)+ Gyv(z)v(z)+ Gym(z)m(z) (7)式中，y(z)、ω(z)、v(z)和m(z)分别对应为yk、ωk、vk和mk的z变换；
Gyω(z)=C(zI－F) ，Gyv(z)=I，Gym(z)=I；
Gyω(z)=OL(z) EL(z) (8)q q‑1 0 q‑1 q‑2 0
其中，OL(z)=zI+ z OL,1+…+ zOL,q，EL(z)=z EL,1+ z EL,2+…+ z EL,q；实常数矩阵OL,μ和EL,μ分别为多项式矩阵OL(z)和EL(z)每一项的系数，其中，μ=1,2,…,q；
令ML=[I OL,1 … OL,q]，NL=[ EL,1 EL,2 … EL,q]，yk‑q=[ yk yk‑1 … yk‑q]；
然后，构造野值检测函数f(k)和野值检测阈值ρ，分别如公式(9)和(10)所示；
f(k)=|| MLyk‑q || (9)ρ=|| NL || qε+|| ML ||( q+1)υ(10)若f(k)<ρ，当前时刻量测值传输至滤波器；否则，进行下一步判断；
步骤2.3. 根据野值发生的间歇性特点，通过判断当前时刻量测是否满足k > t(j－
1)－η且j≥1的条件，从而确定当前时刻量测是否为野值，具体为：若k≤t(j－1)－η，则当前时刻量测值传输至滤波器；
若k>t(j－1)－η，则当前时刻量测值判断为野值，且不被传输；
根据步骤2.1至步骤2.3，当满足触发函数h(k)<0时，则当前时刻量测为“不必要量测”，当同时满足f(k)>ρ和k>t(j－1)－η，当前时刻量测为“野值”；
在上述情况下，当前时刻量测不被传输，触发探测器参数dk=0，根据零阶保持策略，利用最新时刻传输值yks来代替；否则包含有效信息的当前时刻量测值传输到滤波器，dk=1；
因此，k时刻滤波器接收的量测值为zk，zk=dkyk+(1－dk)yks (11)。
4.根据权利要求3所述的野值影响下基于事件触发机制的谐波检测方法，其特征在于，所述步骤3具体为：步骤3.1. 根据谐波检测模型和融合野值检测的事件触发机制，滤波器结构如下：* *
xk+1|k =F xk|k (12)* * * *
xk+1|k+1= xk+1|k+(1－dk)Lk+1(zk+1－yk+1|k) +dkKk+1(yk+1－yk+1|k) (13)* * *其中，x k+1|k和x k+1|k+1分别表示k+1时刻的预测值和估计值，x k|k表示k时刻的估计值，*yk+1|k为k+1时刻的量测预测值，Lk+1和Kk+1分别为不触发和触发情况下的滤波器增益；
步骤3.2. 计算预测误差协方差Pk+1|k和估计误差协方差矩阵Pk+1|k+1；
Pk+1|k=FPk|k+ Qk－dk(FKkSk + SkKkF) (14)T T
Pk+1|k+1= dk[(I－Kk+1C) Pk+1|k(I－Kk+1C) + Kk+1 Rk+1 Kk+1 ]+(1－dk)[(1+β1)×(I－Lk+ T T ‑1 ‑1 T
1C) Pk+1|k(I－Lk+1C) + (1+β2) Lk+1 Rk+1 Lk+1+(1+β1 +β2 ) ×γLk+1 Lk+1] (15)其中，Pk|k为k时刻的估计误差协方差矩阵，Rk+1为k+1时刻的量测噪声协方差矩阵，Kk表示k时刻触发情况下的滤波器增益，β1和β2为给定的正标量；
估计误差协方差上界Pk+1|k+1为：P k+1|k+1=(1+β1) (I－Lk+1C) Pk+1|k(I－Lk+1C) +(1+β2) T ‑1 ‑1 T Lk+1 Rk+1 Lk+1+(1+β1 +β2 )γLk+1 Lk+1 (16)步骤3.3. 通过最小化估计误差协方差矩阵或其上界的迹，计算滤波器增益；
当dk=1，Kk+1= Pk+1|kC (C Pk+1|kC + Rk+1) ；
当dk=0，Lk+1=(1+β1) Pk+1|kC ((1+β1) CPk+1|kC +(1+β2) Rk+1+(1+β1 +β2 )γI) ；
步骤3.4. 计算基波和第i次谐波的幅值Ai,k和相角估计值φ i,k；
其中，xi,k|k和xi+1,k|k分别表示第i维和第i+1维状态估计值；
则相角估计值φ i,k=arctan(ψi,k/ψi+1,k)。
野值影响下基于事件触发机制的谐波检测方法
技术领域
[0001] 本发明属于配电网技术领域，涉及一种野值影响下基于事件触发机制的谐波检测方法。
背景技术
[0002] 近年来，具有波动性大、干扰性强的电力电子装置在电力系统中被广泛使用，分布式电源并退网现象频繁发生，由此引发的谐波问题已经引起学者们的广泛关注。谐波作为一种经典的扰动形式而广泛存在，严重影响着电网的稳定运行。因此，进行快速、准确、有效地谐波检测，对于提高配电网的电能质量具有的重要意义。目前已经提出了多种方法以解决谐波的检测问题。按照不同的角度可分为频域分析法、时频域分析法、参数估计法等。
[0003] 其中，快速傅里叶变换(FFT)是谐波检测领域常用的一种频域分析法，但在应用过程中容易发生频谱泄露和栅栏效应，因此很多研究中提出了改进的FFT方法。时频域分析法可以同时在时域和频域对谐波信号进行检测，例如：现有文献曾提出首先利用经验小波变换对基波和谐波进行划分，然后利用Hilbert变换对谐波的参数进行检测，其最终结果容易受到小波基函数和分解层数的影响。卡尔曼滤波(Kalman filter, KF)作为一种参数估计方法，适用于平稳和非平稳过程，凭借算法简单、精度高等特性已成为故障预测、目标跟踪等领域最基本、最重要的方法之一。到目前为止，KF在谐波检测领域已经被广泛研究。
[0004] 在实际应用过程中，KF的过程噪声协方差矩阵会对滤波性能造成直接影响。针对这一问题，现有技术文献曾提出一种动态追踪模型，与传统模型不同的是，该动态追踪模型重点关注了状态变量之间的相关性，从理论上推导出了过程噪声协方差矩阵，改进了传统模型中过程噪声协方差矩阵为单位阵的缺陷，进一步提高了检测精度。
[0005] 另外，在电网的实际运行过程中，电压信号除了含有周期分量以外，还存在主要以指数形式衰减的非周期分量，这会对电能质量产生较大影响。因此考虑电网中直流衰减分量和谐波同时存在的情况，建立一种具有针对性的动态检测模型具有重要的实际工程意义。
[0006] 众所周知，网络通信资源十分有限。大量量测数据的传输会造成传感器能量的浪费且极大地增加通信信道的压力。在已有的通信机制中，事件触发机制是解决这一问题的有效方案，通过设计合理的触发函数，该机制可以减少传输过程中的冗余数据，从而节省通信资源。
[0007] 值得一提的是，在复杂环境干扰下PMU采集到的测量数据可能会产生突发性误差，即测量野值。考虑到其幅值较大的特点，野值与最新传输值之间的差值很有可能大于预先设定的事件触发阈值，使得野值通过事件触发机制传输到滤波器中，严重降低滤波性能。因此，综合考虑通信资源受限和测量野值的影响，设计一种传输机制，在节约通信资源的同时，能够降低野值对滤波性能的影响，是目前本领域技术人员面临的一个重要挑战。
发明内容
[0008] 本发明的目的在于提出一种野值影响下基于事件触发机制的谐波检测方法，以提高野值和通信资源受限影响下谐波的检测精度。
[0009] 本发明为了实现上述目的，采用如下技术方案：
[0010] 野值影响下基于事件触发机制的谐波检测方法，包括如下步骤：
[0011] 步骤1. 建立谐波动态检测状态空间模型和量测模型；
[0012] 将通过传感器采样得到的电力谐波模拟信号转换为含有直流衰减分量的离散谐波信号；
[0013] 根据含有直流衰减分量的离散谐波信号表达式，设置系统状态变量；给定噪声和野值的假设条件，并计算过程噪声协方差矩阵；建立谐波动态检测状态空间模型和量测模型；
[0014] 步骤2. 建立融合野值检测的事件触发器，将传感器的量测值传输至数据处理中心；
[0015] 具体为：通过构造触发函数以及野值检测函数，判断当前时刻量测值是否满足事件触发条件；若满足事件触发条件，则当前时刻量测值通过通信网络传输至数据处理中心；
[0016] 在数据处理中心设置触发探测器以及滤波器；其中，触发探测器根据通信网络传输过来的数据判断是否发生了事件触发，滤波器用于量测值并估计谐波信号的幅值和相角；
[0017] 若发生了事件触发，则触发探测器赋值为1，将包含有效信息的当前时刻量测值传输到滤波器；否则，触发探测器赋值为0，利用最新时刻传输值代替当前时刻量测值；
[0018] 步骤3. 建立基于融合野值检测事件触发机制的卡尔曼滤波器；
[0019] 建立卡尔曼滤波器结构；根据已建立的谐波动态检测状态空间模型，计算一步预测值和预测误差协方差矩阵，计算估计误差协方差矩阵及其上界；
[0020] 通过最小化估计误差协方差矩阵及其上界的迹得到卡尔曼滤波器的增益参数；最后利用卡尔曼滤波器获得谐波状态的估计值，得到谐波信号的幅值和相角估计值。
[0021] 本发明具有如下优点：
[0022] 如上所述，本发明针对野值和通信资源受限影响下的谐波检测问题，提出了一种野值影响下基于事件触发机制的谐波检测方法。该方法首先针对含有直流衰减分量的谐波信号，采用了一种动态检测模型；然后根据测量野值的特征，利用脉冲函数对野值进行建模，设计了一种融合野值检测的事件触发机制，对量测数据进行筛选，从而只将有效信息传输到滤波器中，在保证检测精度的同时减小了通信信道的压力。在此基础上，本发明还考虑了非触发误差和相关噪声对滤波性能的影响，同时从触发和不触发两种角度来推导滤波器的增益。此外，本发明还给出了具体的仿真实验，通过利用仿真实验验证了本发明所提出方法的有效性。
附图说明
[0023] 图1为本发明实施例中野值影响下基于事件触发机制的谐波检测方法的过程示意图。
[0024] 图2为本发明实施例中建立谐波检测模型的示意图。
[0025] 图3为本发明实施例中融合野值检测的事件触发机制的设计过程示意图。
[0026] 图4为本发明实施例中基于融合野值检测事件触发机制的卡尔曼滤波器的过程示意图。
[0027] 图 5为野值存在时含有直流衰减分量、谐波和噪声的电压信号波形图。
[0028] 图6为KF算法和本发明方法（OD‑ETM‑KF）对x1,k进行检测的效果对比图。
[0029] 图7为KF算法和本发明方法（OD‑ETM‑KF）对x2,k进行检测的效果对比图。
[0030] 图8为KF算法和本发明方法（OD‑ETM‑KF）对A1,k进行检测的效果对比图。
[0031] 图9为本发明方法（OD‑ETM‑KF）下基波幅值的检测效果图。
[0032] 图10为本发明方法下相角检测效果图。
[0033] 图11为本发明方法在三次谐波幅值的均方根误差示意图。
[0034] 图12为本发明方法在三次谐波相角的均方根误差示意图。
[0035] 图13为本发明方法在五次谐波幅值的均方根误差示意图。
[0036] 图14为本发明方法在五次谐波相角的均方根误差示意图。
[0037] 图15为本发明方法直流衰减分量的均方根误差示意图。
[0038] 图16为本发明方法事件触发率及基波幅值和相角的均方根误差随触发阈值的变化图。
具体实施方式
[0039] 下面结合附图以及具体实施方式对本发明作进一步详细说明：
[0040] 本实施例述及了一种野值影响下基于事件触发机制的谐波检测方法，以解决目前野值和通信资源受限影响下谐波的检测精度低的技术问题。
[0041] 基于上述技术问题，首先搭建野值影响下基于事件触发机制的谐波检测系统，其包括传感器、通信网络、数据处理中心等，在数据处理中心设置触发探测器以及滤波器。
[0042] 本实施例中野值影响下基于事件触发机制的谐波检测方法的大致过程如下：
[0043] 将传感器通过采样得到的电力谐波模拟信号转换为含有直流衰减分量的离散谐波信号。
[0044] 然后对传感器配置融合野值检测的事件触发器，通过事件触发器判断是否向数据处理中心传输当前时刻的量测数据，如果满足触发条件则将数据通过通信网络进行传输。
[0045] 另外，由于在数据处理中心设置了触发探测器，因此触发探测器能够根据通信网络传输过来的数据判断是否发生了事件触发，并采用零阶保持策略来更新量测数据。
[0046] 最后将处理后的量测数据传输到滤波器中执行状态估计，最终得到谐波检测结果。
[0047] 下面结合图1对本发明野值影响下基于事件触发机制的谐波检测方法进行详细说明。
[0048] 野值影响下基于事件触发机制的谐波检测方法，包括如下步骤：
[0049] 步骤1. 建立谐波动态检测状态空间模型和量测模型。
[0050] 将通过传感器采样得到的电力谐波模拟信号转换为含有直流衰减分量的离散谐波信号。
[0051] 根据含有直流衰减分量的谐波信号表达式，设置系统状态变量；计算过程噪声协方差矩阵；给定噪声和野值的假设条件；建立谐波动态检测状态空间模型和量测模型。
[0052] 如图2所示，该步骤1具体为：
[0053] 步骤1.1. 首先给出含有直流衰减分量的谐波信号表达式，如公式(1)所示。
[0054] yk=∑ni=1Aisin(iwkT+φi)+Ace‑akT (1)
[0055] 其中，yk表示理想情况下的谐波信号量测值，n表示谐波分量的最高阶数；
[0056] Ai和φi分别表示第i次谐波信号的幅值和相角，i=1,2,…,n，w为角频率，Ace‑akT为直流衰减分量，Ac和a为常数，k和T分别表示采样时刻和采样周期。
[0057] 步骤1.2. 设置系统状态变量xk，如公式(2)所示。
[0058] xk=[ x1,k x2,k… x2n‑1,k x2n,k xc,k] T (2)
[0059] 其中，x2i‑1,k = Aisin(iwkT+φi)，x2i,k = Aiiw cos(iwkT+φi)，xc,k =Ace‑akT。
[0060] 步骤1.3. 给定噪声和野值的假设条件，并计算过程噪声协方差矩阵。
[0061] 设过程噪声ωk=[ω1,k ,ω2,k ,…,ωn,k,ωc,k]T和量测噪声vk；其中，ωi,k表示第i次谐波信号的过程噪声，ωc,k表示直流衰减信号的过程噪声；ωi,k以及ωc,k的表达式如下所示。
[0062] 。
[0063] 。
[0064] 其中，ωi(τ)表示均值为0，且方差为σi2=Ai2/2π的正弦信号的高斯噪声；ωc(τ)表2 2
示均值为0，且方差为σc =Ac /2π的直流衰减信号的高斯噪声；ωi,k、ωc,k和vk满足不等式(3)。
[0065] ||ωi,k|| ≤εi，|ωc,k| ≤εc，|vk| ≤υ (3)
[0066] 其中，εi、εc和υ为给定的正标量，分别对应表示ωi,k、ωc,k和vk的上界。
[0067] 得到 ，||vk|| ≤υ。
[0068] ε是一个标量，表示||ωk||的上界；过程噪声ωk和测量噪声vk相关，统计特性如下：
[0069] E{ωk }=0，E{vk}=0，E{xkωk }=0，E{xkvk }=0；
[0070] Qk=E{ωkωlT }= Qkδk,l，Rk=E{vkvlT }= Rkδk,l，Sk=E{ωkvlT }= Skδk,l。
[0071] 其中，xk表示k时刻的系统状态，ωl和vl分别表示l时刻的过程噪声和量测噪声，δk,l为克罗内克函数，Sk为相关噪声协方差矩阵，Rk为测量噪声协方差矩阵。
[0072] Qk=diag{Qk(1,1),Qk(2,2),…,Qk(n,n),Qc.k}为过程噪声协方差矩阵。
[0073] 。
[0074] 。
[0075] 其中，κ为具有可调范围的参数；定义mk为间歇发生的野值，描述为：
[0076] 。
[0077] 式中，mj,k为第j个野值的幅值，t(j)表示第j个野值出现的时刻。
[0078] 野值出现的时间间隔η和幅值mj,k满足η＞q，||mj,k||＞ζ。
[0079] 其中，η=min{θj}j≥1为野值出现的最小时间间隔，θj=t(j+1)－t(j)，q和ζ为已知正标量。
[0080] 步骤1.4. 考虑过程噪声、量测噪声以及野值的影响，建立谐波动态检测状态空间模型和量测模型，分别如公式(4)和公式(5)所示。
[0081] xk+1 = F xk+ ωk (4)
[0082] yk+1= C xk+1+ vk+1+ mk+1 (5)
[0083] 其中，xk+1和yk+1分别表示k+1时刻的系统状态和量测值，vk +1和mk+1分别表示k+1时刻的量测噪声和野值；F=diag{F1,F2,…,Fn,Fc}为状态转移矩阵。
[0084] ，Fc=e‑aT；C=[1 0 1 0 … 1 0 1]∈R1×(2n+1)为观测矩阵。
[0085] 步骤2. 建立融合野值检测的事件触发器，将传感器的量测值传输至数据处理中心。
[0086] 通过构造触发函数以及野值检测函数，判断当前时刻量测值是否满足事件触发条件；若满足事件触发条件，则当前时刻量测值通过通信网络传输至滤波器。
[0087] 触发探测器根据通信网络传输过来的数据判断是否发生了事件触发。
[0088] 若发生了事件触发，则触发探测器赋值为1，包含有效信息的当前时刻量测值传输到滤波器；否则，触发探测器赋值为0，利用最新时刻传输值代替当前时刻量测值。
[0089] 如图3所示，该步骤2具体为：
[0090] 步骤2.1. 利用当前时刻量测yk、最新时刻传输值yks以及预设触发阈值γ，设计触发函数h(k)，判断当前时刻量测是否为“不必要量测”；触发函数h(k)如公式(6)所示。
[0091] h(k)=(yk－yks)(yk－yks)T－γ (6)
[0092] 其中，ks表示第s个触发时刻，满足k0＜k1＜…＜ks＜…； 若h(k)<0，则当前时刻量测为“不必要量测”且不被传输；否则，进行下一步判断。
[0093] 步骤2.2. 利用最近的q个时刻的量测ykk‑q以及噪声和野值的幅值上界，设计野值检测函数f(k)以及野值检测阈值ρ，初步判断野值。
[0094] 首先，利用传递函数表示公式(4)和(5)中模型的输入输出关系，如公式(7)所示。
[0095] y(z)=Gyω(z)ω(z)+ Gyv(z)v(z)+ Gym(z)m(z) (7)
[0096] 式中，y(z)、ω(z)、v(z)和m(z)分别对应为yk、ωk、vk和mk的z变换。
[0097] Gyω(z)=C(zI－F)‑1，Gyv(z)=I，Gym(z)=I。
[0098] 有理分式矩阵Gyω(z)的左矩阵分式描述为：
[0099] Gyω(z)=OL(z)‑1 EL(z) (8)
[0100] 其中，OL(z)=zqI+ zq‑1OL,1+…+ z0OL,q，EL(z)=zq‑1 EL,1+ zq‑2 EL,2+…+ z0EL,q；实常数矩阵OL,μ和EL,μ分别为多项式矩阵OL(z)和EL(z)每一项的系数，其中，μ=1,2,…,q。
[0101] 令ML=[I OL,1 … OL,q]，NL=[ EL,1 EL,2 … EL,q]，ykk‑q=[ yk yk‑1 … yk‑q]。
[0102] 然后，构造野值检测函数f(k)和野值检测阈值ρ，分别如公式(9)和(10)所示。
[0103] f(k)=|| MLykk‑q || (9)
[0104] ρ=|| NL || qε+|| ML ||( q+1)υ(10)
[0105] 若f(k)<ρ，当前时刻量测值传输至滤波器；否则，进行下一步判断。
[0106] 步骤2.3. 根据野值发生的间歇性特点，通过判断当前时刻量测是否满足k > t(j－1)－η且j≥1的条件，从而确定当前时刻量测是否为野值，具体为：
[0107] 若k≤t(j－1)－η，则当前时刻量测值传输至滤波器；
[0108] 若k>t(j－1)－η，则当前时刻量测值判断为野值，且不被传输。
[0109] 步骤2.4. 触发探测器将探测结果传输至滤波器。
[0110] 根据零阶保持策略，得到k时刻滤波器接收的量测值。
[0111] 根据步骤2.1至步骤2.3，当满足触发函数h(k)<0时，则当前时刻量测为“不必要量测”，当同时满足f(k)>ρ和k>t(j－1)－η，当前时刻量测为“野值”。
[0112] 在上述情况下，当前时刻量测不被传输，触发探测器参数dk=0，根据零阶保持策略，利用最新时刻传输值yks来代替；否则包含有效信息的当前时刻量测值传输到滤波器，dk=1。
[0113] 因此，k时刻滤波器接收的量测值为zk，zk=dkyk+(1－dk)yks (11)。
[0114] 步骤3. 将OD‑ETM引入卡尔曼滤波器KF，设计基于融合野值检测事件触发机制的卡尔曼滤波器OD‑ETM‑KF。该步骤的大致过程如下：
[0115] 建立卡尔曼滤波器结构。根据已建立的状态空间模型计算一步预测值和一步预测值协方差矩阵，计算估计误差协方差矩阵及其上界。
[0116] 通过最小化估计误差协方差矩阵及其上界的迹得到卡尔曼滤波器的增益参数；最后利用卡尔曼滤波器获得谐波状态的估计值，得到谐波信号的幅值和相角估计值。
[0117] 如图4所示，该步骤3具体为：
[0118] 步骤3.1. 根据谐波检测模型和融合野值检测的事件触发机制，滤波器结构如下：
[0119] x*k+1|k =F x*k|k (12)
[0120] x*k+1|k+1= x*k+1|k+(1－dk)Lk+1(zk+1－y*k+1|k) +dkKk+1(yk+1－y*k+1|k) (13)[0121] 其中，x*k+1|k和x*k+1|k+1分别表示k+1时刻的预测值和估计值，x*k|k表示k时刻的估计*值，yk+1|k为k+1时刻的量测预测值，Lk+1和Kk+1分别为不触发和触发情况下的滤波器增益。
[0122] 步骤3.2. 计算预测误差协方差Pk+1|k和估计误差协方差矩阵Pk+1|k+1。
[0123] Pk+1|k=FPk|k+ Qk－dk(FKkSkT+ SkKkTFT) (14)
[0124] Pk+1|k+1= dk[(I－Kk+1C) Pk+1|k(I－Kk+1C)T+ Kk+1 Rk+1 Kk+1T]+(1－dk)[(1+β1)×(I－T T ‑1 ‑1 T
Lk+1C) Pk+1|k(I－Lk+1C) + (1+β2) Lk+1 Rk+1 Lk+1+(1+β1 +β2 ) ×γLk+1 Lk+1] (15)[0125] 其中，Pk|k为k时刻的估计误差协方差矩阵，Rk+1为k+1时刻的量测噪声协方差矩阵，Kk表示k时刻触发情况下的滤波器增益，β1和β2为给定的正标量；
[0126] 估计误差协方差上界P*k+1|k+1为：P*k+1|k+1=(1+β1) (I－Lk+1C) Pk+1|k(I－Lk+1C) T+(1T ‑1 ‑1 T +β2) Lk+1 Rk+1 Lk+1+(1+β1 +β2 )γLk+1 Lk+1 (16)
[0127] 步骤3.3. 通过最小化估计误差协方差矩阵或其上界的迹，计算滤波器增益。
[0128] 当dk=1，Kk+1= Pk+1|kCT(C Pk+1|kCT+ Rk+1)‑1；
[0129] 当dk=0，Lk+1=(1+β1) Pk+1|kCT((1+β1) CPk+1|kCT+(1+β2) Rk+1+(1+β1‑1+β2‑1)γI) ‑1。
[0130] 步骤3.4. 计算基波和第i次谐波的幅值A*i,k和相角估计值φ*i,k。
[0131] 。
[0132] 其中，x*i,k|k和x*i+1,k|k分别表示第i维和第i+1维状态估计值。
[0133] 令 。
[0134] 则相角估计值φ*i,k=arctan(ψi,k/ψi+1,k)。
[0135] 此外，为了验证本发明方法的有效性，本发明还给出了如下实验：
[0136] 由于配电网中偶次谐波及高次的奇次谐波含量很少，本发明主要对基波、三次谐波、五次谐波以及直流衰减分量进行检测，并采用MATLAB软件进行仿真和性能分析。
[0137] 为了使检测结果更具一般性，本发明在进行50次蒙特卡罗实验的基础上进行性能对比分析，以均方根误差(root mean square error, RMSE)作为评价指标，均方根误差的定义如下：
[0138] 。
[0139] 其中，sk和s*k分别为参数的真实值和估计值，sÎ{A,φ}，N为蒙特卡罗实验次数。
[0140] 均方根误差RMSE越小代表检测精度越高。
[0141] 为了衡量OD‑ETM下的数据传输量，本发明采用触发率作为评价指标，触发率定义为：PET=DET / Dy，DET为满足事件触发条件的数据传输量，Dy为传感器测量的数据总量。
[0142] 参数设置如下：设带有直流衰减分量的谐波电压信号为：
[0143] uk=220√2sin(wkT)+0.4´220√2sin(3wkT+π/3)+0.2´220√2sin(5wkT+π/6)+10e‑akT+vk。
[0144] 直流衰减系数取典型值20，采样频率为104HZ，采样点为1000个，P0|0=0.001I7×7，κ取0.05，观测噪声设置为vk=Hkωk+ξk，其中Hk=row{1}2n+1，ξk是方差为0.01的白噪声。
[0145] β1和β2分别设置为0.001和0.1，ε和υ分别为0.04和0.96，q为8，野值幅值设置为2.1ρ，野值发生的时间间隔设置为满足概率分布的随机序列：
[0146] prob{θj=12}=0.3，prob{θj=38}=0.4；prob{θj=48}=0.1，prob{θj=55}=0.2。
[0147] 实验结果如图5至图16所示。
[0148] 图5展示了野值存在时含有直流衰减分量、谐波和噪声的电压信号。
[0149] 以基波信号为例，分别采用KF和OD‑ETM‑KF对其状态变量和参数进行检测，并将检测效果进行对比分析，触发阈值设置为60。
[0150] 图6、图7、图8分别为采用两种算法对状态变量x1,k和x2,k以及A1,k检测的效果对比图。
[0151] 由图6至图8可以看出，野值影响下采用KF的检测效果较差,甚至会引起滤波发散，尤其在对x2,k和A1,k进行检测时表现较为明显，主要是因为传统KF中没有考虑到野值对系统性能的影响，而采用本发明OD‑ETM‑KF的检测结果更加接近真实值。
[0152] 图9和图10分别为野值影响下采用OD‑ETM‑KF对A1,k和φ1,k进行检测的效果图，看出采用OD‑ETM‑KF能够降低野值产生的影响，在野值出现时刻仍然保持较高的检测精度。
[0153] 图11至图15以均方根误差为评价指标，来衡量采用OD‑ETM‑KF对直流衰减分量及谐波参数的检测性能。其中，图11至图14为分别对三次谐波、五次谐波的参数进行检测的均方根误差，图15为对直流衰减分量进行检测的均方根误差。
[0154] 通过分析仿真实验结果得到：采用OD‑ETM‑KF在触发阈值设置为60时，数据传输量达到59.3%，从而保证只将有效信息传输到滤波器中，节省了通信资源。
[0155] 当量测中存在野值时，采用OD‑ETM‑KF对对直流衰减分量及谐波参数进行检测的均方根误差都在允许范围内，因此，该方法还具有一定的抗野值性能。
[0156] 图16为触发率及A1,k和φ1,k的均方根误差随触发阈值的变化情况。
[0157] 从图中能够看出随着触发阈值的增大，A1,k和φ1,k的均方根误差略有增大，但触发率不断降低.在数据传输率为53.2%的时候A1,k的检测误差只有0.0082，仍能保持在允许范围内。
[0158] 该触发机制能够将几乎所有的有效信息传输到滤波器中，保证检测精度，因此，根据实际需求合理选取触发阈值，在保证检测精度的前提下减小数据传输量，节省通信资源。
[0159] 当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。
