CN115906535B - 野值影响下基于事件触发机制的谐波检测方法 - Google Patents

野值影响下基于事件触发机制的谐波检测方法 Download PDF

Info

Publication number
CN115906535B
CN115906535B CN202310022658.1A CN202310022658A CN115906535B CN 115906535 B CN115906535 B CN 115906535B CN 202310022658 A CN202310022658 A CN 202310022658A CN 115906535 B CN115906535 B CN 115906535B
Authority
CN
China
Prior art keywords
value
measurement
harmonic
time
detection
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
CN202310022658.1A
Other languages
English (en)
Other versions
CN115906535A (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.)
Shandong University of Science and Technology
Original Assignee
Shandong University of Science and Technology
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 Shandong University of Science and Technology filed Critical Shandong University of Science and Technology
Priority to CN202310022658.1A priority Critical patent/CN115906535B/zh
Publication of CN115906535A publication Critical patent/CN115906535A/zh
Application granted granted Critical
Publication of CN115906535B publication Critical patent/CN115906535B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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/40—Arrangements for reducing harmonics

Landscapes

  • Geophysics And Detection Of Objects (AREA)
  • Measurement Of Current Or Voltage (AREA)

Abstract

本发明属于配电网技术领域，具体公开了一种野值影响下基于事件触发机制的谐波检测方法，以解决目前针对野值和通信资源受限影响下谐波检测精度不高的技术问题。本发明方法包括如下步骤：首先，针对含有直流衰减分量的谐波信号采用了一种动态检测模型，以提高谐波检测的精度；其次，为了降低系统性能对野值的敏感性并达到节省网络通信资源的目的，设计了一种融合野值检测的事件触发器；在此基础上，将事件触发机制引入卡尔曼滤波器，并充分考虑非触发误差对滤波性能的影响，以最小化误差协方差矩阵的上界为目标，获得滤波器的增益。本发明方法对于直流衰减分量及谐波的参数具有较高的检测精度。

Description

野值影响下基于事件触发机制的谐波检测方法
技术领域
本发明属于配电网技术领域，涉及一种野值影响下基于事件触发机制的谐波检测方法。
背景技术
近年来，具有波动性大、干扰性强的电力电子装置在电力系统中被广泛使用，分布式电源并退网现象频繁发生，由此引发的谐波问题已经引起学者们的广泛关注。谐波作为一种经典的扰动形式而广泛存在，严重影响着电网的稳定运行。因此，进行快速、准确、有效地谐波检测，对于提高配电网的电能质量具有的重要意义。目前已经提出了多种方法以解决谐波的检测问题。按照不同的角度可分为频域分析法、时频域分析法、参数估计法等。
其中，快速傅里叶变换(FFT)是谐波检测领域常用的一种频域分析法，但在应用过程中容易发生频谱泄露和栅栏效应，因此很多研究中提出了改进的FFT方法。时频域分析法可以同时在时域和频域对谐波信号进行检测，例如：现有文献曾提出首先利用经验小波变换对基波和谐波进行划分，然后利用Hilbert变换对谐波的参数进行检测，其最终结果容易受到小波基函数和分解层数的影响。卡尔曼滤波(Kalman filter, KF)作为一种参数估计方法，适用于平稳和非平稳过程，凭借算法简单、精度高等特性已成为故障预测、目标跟踪等领域最基本、最重要的方法之一。到目前为止，KF在谐波检测领域已经被广泛研究。
在实际应用过程中，KF的过程噪声协方差矩阵会对滤波性能造成直接影响。针对这一问题，现有技术文献曾提出一种动态追踪模型，与传统模型不同的是，该动态追踪模型重点关注了状态变量之间的相关性，从理论上推导出了过程噪声协方差矩阵，改进了传统模型中过程噪声协方差矩阵为单位阵的缺陷，进一步提高了检测精度。
另外，在电网的实际运行过程中，电压信号除了含有周期分量以外，还存在主要以指数形式衰减的非周期分量，这会对电能质量产生较大影响。因此考虑电网中直流衰减分量和谐波同时存在的情况，建立一种具有针对性的动态检测模型具有重要的实际工程意义。
众所周知，网络通信资源十分有限。大量量测数据的传输会造成传感器能量的浪费且极大地增加通信信道的压力。在已有的通信机制中，事件触发机制是解决这一问题的有效方案，通过设计合理的触发函数，该机制可以减少传输过程中的冗余数据，从而节省通信资源。
值得一提的是，在复杂环境干扰下PMU采集到的测量数据可能会产生突发性误差，即测量野值。考虑到其幅值较大的特点，野值与最新传输值之间的差值很有可能大于预先设定的事件触发阈值，使得野值通过事件触发机制传输到滤波器中，严重降低滤波性能。因此，综合考虑通信资源受限和测量野值的影响，设计一种传输机制，在节约通信资源的同时，能够降低野值对滤波性能的影响，是目前本领域技术人员面临的一个重要挑战。
发明内容
本发明的目的在于提出一种野值影响下基于事件触发机制的谐波检测方法，以提高野值和通信资源受限影响下谐波的检测精度。
本发明为了实现上述目的，采用如下技术方案：
野值影响下基于事件触发机制的谐波检测方法，包括如下步骤：
步骤1. 建立谐波动态检测状态空间模型和量测模型；
将通过传感器采样得到的电力谐波模拟信号转换为含有直流衰减分量的离散谐波信号；
根据含有直流衰减分量的离散谐波信号表达式，设置系统状态变量；给定噪声和野值的假设条件，并计算过程噪声协方差矩阵；建立谐波动态检测状态空间模型和量测模型；
步骤2. 建立融合野值检测的事件触发器，将传感器的量测值传输至数据处理中心；
具体为：通过构造触发函数以及野值检测函数，判断当前时刻量测值是否满足事件触发条件；若满足事件触发条件，则当前时刻量测值通过通信网络传输至数据处理中心；
在数据处理中心设置触发探测器以及滤波器；其中，触发探测器根据通信网络传输过来的数据判断是否发生了事件触发，滤波器用于量测值并估计谐波信号的幅值和相角；
若发生了事件触发，则触发探测器赋值为1，将包含有效信息的当前时刻量测值传输到滤波器；否则，触发探测器赋值为0，利用最新时刻传输值代替当前时刻量测值；
步骤3. 建立基于融合野值检测事件触发机制的卡尔曼滤波器；
建立卡尔曼滤波器结构；根据已建立的谐波动态检测状态空间模型，计算一步预测值和预测误差协方差矩阵，计算估计误差协方差矩阵及其上界；
通过最小化估计误差协方差矩阵及其上界的迹得到卡尔曼滤波器的增益参数；最后利用卡尔曼滤波器获得谐波状态的估计值，得到谐波信号的幅值和相角估计值。
本发明具有如下优点：
如上所述，本发明针对野值和通信资源受限影响下的谐波检测问题，提出了一种野值影响下基于事件触发机制的谐波检测方法。该方法首先针对含有直流衰减分量的谐波信号，采用了一种动态检测模型；然后根据测量野值的特征，利用脉冲函数对野值进行建模，设计了一种融合野值检测的事件触发机制，对量测数据进行筛选，从而只将有效信息传输到滤波器中，在保证检测精度的同时减小了通信信道的压力。在此基础上，本发明还考虑了非触发误差和相关噪声对滤波性能的影响，同时从触发和不触发两种角度来推导滤波器的增益。此外，本发明还给出了具体的仿真实验，通过利用仿真实验验证了本发明所提出方法的有效性。
附图说明
图1为本发明实施例中野值影响下基于事件触发机制的谐波检测方法的过程示意图。
图2为本发明实施例中建立谐波检测模型的示意图。
图3为本发明实施例中融合野值检测的事件触发机制的设计过程示意图。
图4为本发明实施例中基于融合野值检测事件触发机制的卡尔曼滤波器的过程示意图。
图 5为野值存在时含有直流衰减分量、谐波和噪声的电压信号波形图。
图6为KF算法和本发明方法（OD-ETM-KF）对x 1,k 进行检测的效果对比图。
图7为KF算法和本发明方法（OD-ETM-KF）对x 2,k 进行检测的效果对比图。
图8为KF算法和本发明方法（OD-ETM-KF）对A 1,k 进行检测的效果对比图。
图9为本发明方法（OD-ETM-KF）下基波幅值的检测效果图。
图10为本发明方法下相角检测效果图。
图11为本发明方法在三次谐波幅值的均方根误差示意图。
图12为本发明方法在三次谐波相角的均方根误差示意图。
图13为本发明方法在五次谐波幅值的均方根误差示意图。
图14为本发明方法在五次谐波相角的均方根误差示意图。
图15为本发明方法直流衰减分量的均方根误差示意图。
图16为本发明方法事件触发率及基波幅值和相角的均方根误差随触发阈值的变化图。
具体实施方式
下面结合附图以及具体实施方式对本发明作进一步详细说明：
本实施例述及了一种野值影响下基于事件触发机制的谐波检测方法，以解决目前野值和通信资源受限影响下谐波的检测精度低的技术问题。
基于上述技术问题，首先搭建野值影响下基于事件触发机制的谐波检测系统，其包括传感器、通信网络、数据处理中心等，在数据处理中心设置触发探测器以及滤波器。
本实施例中野值影响下基于事件触发机制的谐波检测方法的大致过程如下：
将传感器通过采样得到的电力谐波模拟信号转换为含有直流衰减分量的离散谐波信号。
然后对传感器配置融合野值检测的事件触发器，通过事件触发器判断是否向数据处理中心传输当前时刻的量测数据，如果满足触发条件则将数据通过通信网络进行传输。
另外，由于在数据处理中心设置了触发探测器，因此触发探测器能够根据通信网络传输过来的数据判断是否发生了事件触发，并采用零阶保持策略来更新量测数据。
最后将处理后的量测数据传输到滤波器中执行状态估计，最终得到谐波检测结果。
下面结合图1对本发明野值影响下基于事件触发机制的谐波检测方法进行详细说明。
野值影响下基于事件触发机制的谐波检测方法，包括如下步骤：
步骤1. 建立谐波动态检测状态空间模型和量测模型。
将通过传感器采样得到的电力谐波模拟信号转换为含有直流衰减分量的离散谐波信号。
根据含有直流衰减分量的谐波信号表达式，设置系统状态变量；计算过程噪声协方差矩阵；给定噪声和野值的假设条件；建立谐波动态检测状态空间模型和量测模型。
如图2所示，该步骤1具体为：
步骤1.1. 首先给出含有直流衰减分量的谐波信号表达式，如公式(1)所示。
y k =∑ n i=1 A i sin(iwkT+φ i )+A c e -akT (1)
其中，y k 表示理想情况下的谐波信号量测值，n表示谐波分量的最高阶数；
A i 和φ i 分别表示第i次谐波信号的幅值和相角，i=1,2,…,n，w为角频率，A c e -akT 为直流衰减分量，A c 和a为常数，k和T分别表示采样时刻和采样周期。
步骤1.2. 设置系统状态变量x k ，如公式(2)所示。
x k =[ x 1,k x 2,k … x 2n-1,k x 2n,k x c,k ] T (2)
其中，x 2i-1,k = A i sin(iwkT+φ i )，x 2i,k = A i iw cos(iwkT+φ i )，x c,k =A c e -akT 。
步骤1.3. 给定噪声和野值的假设条件，并计算过程噪声协方差矩阵。
设过程噪声ω k =[ω 1,k ,ω 2,k ,…,ω n,k ,ω c,k ]T和量测噪声v k ；其中，ω i,k 表示第i次谐波信号的过程噪声，ω c,k 表示直流衰减信号的过程噪声；ω i,k 以及ω c,k 的表达式如下所示。
Figure 133751DEST_PATH_IMAGE001
。
Figure 153660DEST_PATH_IMAGE002
。
其中，ω i (τ)表示均值为0，且方差为σ i 2=A i 2/2π的正弦信号的高斯噪声；ω c (τ)表示均值为0，且方差为σ c 2=A c 2/2π的直流衰减信号的高斯噪声；ω i , k 、ω c,k 和v k 满足不等式(3)。
||ω i,k || ≤ε i ，|ω c,k | ≤ε c ，|v k | ≤υ (3)
其中，ε i 、ε c 和υ为给定的正标量，分别对应表示ω i , k 、ω c,k 和v k 的上界。
得到
Figure 156251DEST_PATH_IMAGE003
，||v k || ≤υ。
ε是一个标量，表示||ω k ||的上界；过程噪声ω k 和测量噪声v k 相关，统计特性如下：
E{ω k }=0，E{v k }=0，E{x k ω k }=0，E{x k v k }=0；
Q k =E{ω k ω l T }= Q k δ k,l ，R k =E{v k v l T }= R k δ k,l ，S k =E{ω k v l T }= S k δ k,l 。
其中，x k 表示k时刻的系统状态，ω l 和v l 分别表示l时刻的过程噪声和量测噪声，δ k,l 为克罗内克函数，S k 为相关噪声协方差矩阵，R k 为测量噪声协方差矩阵。
Q k =diag{Q k (1,1),Q k (2,2),…,Q k (n,n),Q c.k }为过程噪声协方差矩阵。
Figure 945215DEST_PATH_IMAGE004
。
Figure 654021DEST_PATH_IMAGE005
。
其中，κ为具有可调范围的参数；定义m k 为间歇发生的野值，描述为：
Figure 161225DEST_PATH_IMAGE006
。
式中，m j,k 为第j个野值的幅值，t(j)表示第j个野值出现的时刻。
野值出现的时间间隔η和幅值m j,k 满足η＞q，||m j,k ||＞ζ。
其中，η=min{θ j } j≥1为野值出现的最小时间间隔，θ j =t(j+1)－t(j)，q和ζ为已知正标量。
步骤1.4. 考虑过程噪声、量测噪声以及野值的影响，建立谐波动态检测状态空间模型和量测模型，分别如公式(4)和公式(5)所示。
x k+1 = F x k + ω k (4)
y k+1= C x k+1+ v k+1+ m k+1 (5)
其中，x k+1和y k+1分别表示k+1时刻的系统状态和量测值，v k +1和m k+1分别表示k+1时刻的量测噪声和野值；F=diag{F 1,F 2,…,F n ,F c }为状态转移矩阵。
Figure 701928DEST_PATH_IMAGE007
，F c =e-aT ；C=[1 0 1 0 … 1 0 1]∈R1×(2n+1)为观测矩 阵。
步骤2. 建立融合野值检测的事件触发器，将传感器的量测值传输至数据处理中心。
通过构造触发函数以及野值检测函数，判断当前时刻量测值是否满足事件触发条件；若满足事件触发条件，则当前时刻量测值通过通信网络传输至滤波器。
触发探测器根据通信网络传输过来的数据判断是否发生了事件触发。
若发生了事件触发，则触发探测器赋值为1，包含有效信息的当前时刻量测值传输到滤波器；否则，触发探测器赋值为0，利用最新时刻传输值代替当前时刻量测值。
如图3所示，该步骤2具体为：
步骤2.1. 利用当前时刻量测y k 、最新时刻传输值y ks 以及预设触发阈值γ，设计触发函数h(k)，判断当前时刻量测是否为“不必要量测”；触发函数h(k)如公式(6)所示。
h(k)=(y k －y ks )(y k －y ks ) T －γ (6)
其中，ks表示第s个触发时刻，满足k0＜k1＜…＜ks＜…； 若h(k)<0，则当前时刻量测为“不必要量测”且不被传输；否则，进行下一步判断。
步骤2.2. 利用最近的q个时刻的量测y k k-q 以及噪声和野值的幅值上界，设计野值检测函数f(k)以及野值检测阈值ρ，初步判断野值。
首先，利用传递函数表示公式(4)和(5)中模型的输入输出关系，如公式(7)所示。
y(z)=G yω (z)ω(z)+ G yv (z)v(z)+ G ym (z)m(z) (7)
式中，y(z)、ω(z)、v(z)和m(z)分别对应为y k 、ω k 、v k 和m k 的z变换。
G yω (z)=C(zI－F)-1，G yv (z)=I，G ym (z)=I。
有理分式矩阵G yω (z)的左矩阵分式描述为：
G yω (z)=O L (z)-1 E L (z) (8)
其中，O L (z)=z q I+ z q-1 O L,1 +…+ z 0 O L,q ，E L (z)=z q-1 E L,1+ z q-2 E L,2 +…+ z 0 E L,q ；实常数矩阵O L,μ 和E L,μ 分别为多项式矩阵O L (z)和E L (z)每一项的系数，其中，μ=1,2,…,q。
令M L =[I O L,1 … O L,q ]，N L =[ E L,1 E L,2 … E L,q ]，y k k-q =[ y k y k-1 … y k-q ]。
然后，构造野值检测函数f(k)和野值检测阈值ρ，分别如公式(9)和(10)所示。
f(k)=|| M L y k k-q || (9)
ρ=|| N L || qε+|| M L ||( q+1)υ(10)
若f(k)<ρ，当前时刻量测值传输至滤波器；否则，进行下一步判断。
步骤2.3. 根据野值发生的间歇性特点，通过判断当前时刻量测是否满足k > t(j－1)－η且j≥1的条件，从而确定当前时刻量测是否为野值，具体为：
若k≤t(j－1)－η，则当前时刻量测值传输至滤波器；
若k>t(j－1)－η，则当前时刻量测值判断为野值，且不被传输。
步骤2.4. 触发探测器将探测结果传输至滤波器。
根据零阶保持策略，得到k时刻滤波器接收的量测值。
根据步骤2.1至步骤2.3，当满足触发函数h(k)<0时，则当前时刻量测为“不必要量测”，当同时满足f(k)>ρ和k>t(j－1)－η，当前时刻量测为“野值”。
在上述情况下，当前时刻量测不被传输，触发探测器参数d k =0，根据零阶保持策略，利用最新时刻传输值y ks 来代替；否则包含有效信息的当前时刻量测值传输到滤波器，d k =1。
因此，k时刻滤波器接收的量测值为z k ，z k =d k y k +(1－d k )y ks (11)。
步骤3. 将OD-ETM引入卡尔曼滤波器KF，设计基于融合野值检测事件触发机制的卡尔曼滤波器OD-ETM-KF。该步骤的大致过程如下：
建立卡尔曼滤波器结构。根据已建立的状态空间模型计算一步预测值和一步预测值协方差矩阵，计算估计误差协方差矩阵及其上界。
通过最小化估计误差协方差矩阵及其上界的迹得到卡尔曼滤波器的增益参数；最后利用卡尔曼滤波器获得谐波状态的估计值，得到谐波信号的幅值和相角估计值。
如图4所示，该步骤3具体为：
步骤3.1. 根据谐波检测模型和融合野值检测的事件触发机制，滤波器结构如下：
x * k+1|k =F x * k|k (12)
x * k+1|k+1= x * k+1|k +(1－d k )L k+1(z k+1 －y * k+1|k ) +d k K k+1(y k+1 －y * k+1|k ) (13)
其中，x * k+1|k 和x * k+1|k+1分别表示k+1时刻的预测值和估计值，x * k|k 表示k时刻的估计值，y * k+1|k 为k+1时刻的量测预测值，L k+1和K k+1分别为不触发和触发情况下的滤波器增益。
步骤3.2. 计算预测误差协方差P k+1|k 和估计误差协方差矩阵P k+1|k+1。
P k+1|k =FP k|k + Q k －d k (FK k S k T + S k K k T F T ) (14)
P k+1|k+1= d k [(I－K k+1 C) P k+1|k (I－K k+1 C) T + K k+1 R k+1 K k+1 T ]+(1－d k )[(1+β 1)×(I－ L k+1 C) P k+1|k (I－L k+1 C) T + (1+β 2) L k+1 R k+1 L k+1 T +(1+β 1 -1+β 2 -1) ×γL k+1 L k+1 T ] (15)
其中，P k|k 为k时刻的估计误差协方差矩阵，R k+1为k+1时刻的量测噪声协方差矩阵，K k 表示k时刻触发情况下的滤波器增益，β 1和β 2为给定的正标量；
估计误差协方差上界P * k+1|k+1为：P * k+1|k+1=(1+β 1) (I－L k+1 C) P k+1|k (I－L k+1 C) T +(1+β 2) L k+1 R k+1 L k+1 T +(1+β 1 -1+β 2 -1)γL k+1 L k+1 T (16)
步骤3.3. 通过最小化估计误差协方差矩阵或其上界的迹，计算滤波器增益。
当d k =1，K k+1= P k+1|k C T (C P k+1|k C T + R k+1)-1；
当d k =0，L k+1=(1+β 1) P k+1|k C T ((1+β 1) CP k+1|k C T +(1+β 2) R k+1+(1+β 1 -1+β 2 -1)γI) -1。
步骤3.4. 计算基波和第i次谐波的幅值A * i,k 和相角估计值φ * i,k 。
Figure 876558DEST_PATH_IMAGE008
。
其中，x * i,k|k 和x * i+1,k|k 分别表示第i维和第i+1维状态估计值。
令
Figure 477303DEST_PATH_IMAGE009
。
则相角估计值φ * i,k =arctan(ψ i,k /ψ i+1,k )。
此外，为了验证本发明方法的有效性，本发明还给出了如下实验：
由于配电网中偶次谐波及高次的奇次谐波含量很少，本发明主要对基波、三次谐波、五次谐波以及直流衰减分量进行检测，并采用MATLAB软件进行仿真和性能分析。
为了使检测结果更具一般性，本发明在进行50次蒙特卡罗实验的基础上进行性能对比分析，以均方根误差(root mean square error, RMSE)作为评价指标，均方根误差的定义如下：
Figure 471804DEST_PATH_IMAGE010
。
其中，s k 和s * k 分别为参数的真实值和估计值，sÎ{A,φ}，N为蒙特卡罗实验次数。
均方根误差RMSE越小代表检测精度越高。
为了衡量OD-ETM下的数据传输量，本发明采用触发率作为评价指标，触发率定义为：P ET =D ET / D y ，D ET 为满足事件触发条件的数据传输量，D y 为传感器测量的数据总量。
参数设置如下：设带有直流衰减分量的谐波电压信号为：
u k =220√2sin(wkT)+0.4´220√2sin(3wkT+π/3)+0.2´220√2sin(5wkT+π/6)+10e-akT +v k 。
直流衰减系数取典型值20，采样频率为104HZ，采样点为1000个，P 0|0=0.001I 7×7，κ取0.05，观测噪声设置为v k =H k ω k +ξ k ，其中H k =row{1}2n+1，ξ k 是方差为0.01的白噪声。
β 1和β 2分别设置为0.001和0.1，ε和υ分别为0.04和0.96，q为8，野值幅值设置为2.1ρ，野值发生的时间间隔设置为满足概率分布的随机序列：
prob{θ j =12}=0.3，prob{θ j =38}=0.4；prob{θ j =48}=0.1，prob{θ j =55}=0.2。
实验结果如图5至图16所示。
图5展示了野值存在时含有直流衰减分量、谐波和噪声的电压信号。
以基波信号为例，分别采用KF和OD-ETM-KF对其状态变量和参数进行检测，并将检测效果进行对比分析，触发阈值设置为60。
图6、图7、图8分别为采用两种算法对状态变量x 1,k 和x 2,k 以及A 1,k 检测的效果对比图。
由图6至图8可以看出，野值影响下采用KF的检测效果较差,甚至会引起滤波发散，尤其在对x 2,k 和A 1,k 进行检测时表现较为明显，主要是因为传统KF中没有考虑到野值对系统性能的影响，而采用本发明OD-ETM-KF的检测结果更加接近真实值。
图9和图10分别为野值影响下采用OD-ETM-KF对A 1,k 和φ 1,k 进行检测的效果图，看出采用OD-ETM-KF能够降低野值产生的影响，在野值出现时刻仍然保持较高的检测精度。
图11至图15以均方根误差为评价指标，来衡量采用OD-ETM-KF对直流衰减分量及谐波参数的检测性能。其中，图11至图14为分别对三次谐波、五次谐波的参数进行检测的均方根误差，图15为对直流衰减分量进行检测的均方根误差。
通过分析仿真实验结果得到：采用OD-ETM-KF在触发阈值设置为60时，数据传输量达到59.3%，从而保证只将有效信息传输到滤波器中，节省了通信资源。
当量测中存在野值时，采用OD-ETM-KF对对直流衰减分量及谐波参数进行检测的均方根误差都在允许范围内，因此，该方法还具有一定的抗野值性能。
图16为触发率及A 1,k 和φ 1,k 的均方根误差随触发阈值的变化情况。
从图中能够看出随着触发阈值的增大，A 1,k 和φ 1,k 的均方根误差略有增大，但触发率不断降低.在数据传输率为53.2%的时候A 1,k 的检测误差只有0.0082，仍能保持在允许范围内。
该触发机制能够将几乎所有的有效信息传输到滤波器中，保证检测精度，因此，根据实际需求合理选取触发阈值，在保证检测精度的前提下减小数据传输量，节省通信资源。
当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。

Claims (4)

1.野值影响下基于事件触发机制的谐波检测方法，其特征在于，包括如下步骤：
步骤1. 建立谐波动态检测状态空间模型和量测模型；
将通过传感器采样得到的电力谐波模拟信号转换为含有直流衰减分量的离散谐波信号；
根据含有直流衰减分量的离散谐波信号表达式，设置系统状态变量；给定噪声和野值的假设条件，并计算过程噪声协方差矩阵；建立谐波动态检测状态空间模型和量测模型；
步骤2. 建立融合野值检测的事件触发器，将传感器的量测值传输至数据处理中心；
具体为：通过构造触发函数以及野值检测函数，判断当前时刻量测值是否满足事件触发条件；若满足事件触发条件，则当前时刻量测值通过通信网络传输至数据处理中心；
在数据处理中心设置触发探测器以及滤波器；其中，触发探测器根据通信网络传输过来的数据判断是否发生了事件触发，滤波器用于量测值并估计谐波信号的幅值和相角；
若发生了事件触发，则触发探测器赋值为1，将包含有效信息的当前时刻量测值传输到滤波器；否则，触发探测器赋值为0，利用最新时刻传输值代替当前时刻量测值；
步骤3. 建立基于融合野值检测事件触发机制的卡尔曼滤波器；
建立卡尔曼滤波器结构；根据已建立的谐波动态检测状态空间模型，计算一步预测值和预测误差协方差矩阵，计算估计误差协方差矩阵及其上界；
通过最小化估计误差协方差矩阵及其上界的迹得到卡尔曼滤波器的增益参数；最后利用卡尔曼滤波器获得谐波状态的估计值，得到谐波信号的幅值和相角估计值。
2.根据权利要求1所述的野值影响下基于事件触发机制的谐波检测方法，其特征在于，
所述步骤1具体为：
步骤1.1. 首先给出含有直流衰减分量的谐波信号表达式，如公式(1)所示；
y k =∑ n i=1 A i sin(iwkT+φ i )+A c e -akT (1)
其中，y k 表示理想情况下的谐波信号量测值，n表示谐波分量的最高阶数；
A i 和φ i 分别表示第i次谐波信号的幅值和相角，i=1,2,…,n，w为角频率，A c e -akT 为直流衰减分量，A c 和a为常数，k和T分别表示采样时刻和采样周期；
步骤1.2. 设置系统状态变量x k ，如公式(2)所示；
x k =[ x 1,k x 2,k … x 2n-1,k x 2n,k x c,k ] T (2)
其中，x 2i-1,k = A i sin(iwkT+φ i )，x 2i,k = A i iw cos(iwkT+φ i )，x c,k =A c e -akT ；
步骤1.3. 给定噪声和野值的假设条件，并计算过程噪声协方差矩阵；
设过程噪声ω k =[ω 1,k ,ω 2,k ,…,ω n,k ,ω c,k ]T和量测噪声v k ；其中，ω i,k 表示第i次谐波信号的过程噪声，ω c,k 表示直流衰减信号的过程噪声；ω i,k 以及ω c,k 的表达式如下所示；
Figure 389843DEST_PATH_IMAGE001
；
Figure 585945DEST_PATH_IMAGE002
；
其中，ω i (τ)表示均值为0，且方差为σ i 2=A i 2/2π的正弦信号的高斯噪声；ω c (τ)表示均值为0，且方差为σ c 2=A c 2/2π的直流衰减信号的高斯噪声；ω i , k 、ω c,k 和v k 满足不等式(3)；
||ω i,k || ≤ε i ，|ω c,k | ≤ε c ，|v k | ≤υ (3)
其中，ε i 、ε c 和υ为给定的正标量，分别对应表示ω i , k 、ω c,k 和v k 的上界；
得到
Figure 631261DEST_PATH_IMAGE003
，||v k || ≤υ；
ε是一个标量，表示||ω k ||的上界；过程噪声ω k 和测量噪声v k 相关，统计特性如下：
E{ω k }=0，E{v k }=0，E{x k ω k }=0，E{x k v k }=0；
Q k =E{ω k ω l T }= Q k δ k,l ，R k =E{v k v l T }= R k δ k,l ，S k =E{ω k v l T }= S k δ k,l ；
其中，x k 表示k时刻的系统状态，ω l 和v l 分别表示l时刻的过程噪声和量测噪声，δ k,l 为克罗内克函数，S k 为相关噪声协方差矩阵，R k 为测量噪声协方差矩阵；
Q k =diag{Q k (1,1),Q k (2,2),…,Q k (n,n),Q c.k }为过程噪声协方差矩阵；
Figure 557629DEST_PATH_IMAGE004
；
Figure 637580DEST_PATH_IMAGE005
；
其中，κ为具有可调范围的参数；定义m k 为间歇发生的野值，描述为：
Figure 725622DEST_PATH_IMAGE006
；
式中，m j,k 为第j个野值的幅值，t(j)表示第j个野值出现的时刻；
野值出现的时间间隔η和幅值m j,k 满足η＞q，||m j,k ||＞ζ；
其中，η=min{θ j } j≥1为野值出现的最小时间间隔，θ j =t(j+1)－t(j)，q和ζ为已知正标量；
步骤1.4. 考虑过程噪声、量测噪声以及野值的影响，建立谐波动态检测状态空间模型和量测模型，分别如公式(4)和公式(5)所示；
x k+1 = F x k + ω k (4)
y k+1= C x k+1+ v k+1+ m k+1 (5)
其中，x k+1和y k+1分别表示k+1时刻的系统状态和量测值，v k +1和m k+1分别表示k+1时刻的量测噪声和野值；F=diag{F 1,F 2,…,F n ,F c }为状态转移矩阵；
Figure 523814DEST_PATH_IMAGE007
，F c =e-aT ；C=[1 0 1 0 … 1 0 1]∈R1×(2n+1)为观测矩阵。
3.根据权利要求2所述的野值影响下基于事件触发机制的谐波检测方法，其特征在于，
所述步骤2具体为：
步骤2.1. 利用当前时刻量测y k 、最新时刻传输值y ks 以及预设触发阈值γ，设计触发函数h(k)，判断当前时刻量测是否为“不必要量测”；触发函数h(k)如公式(6)所示；
h(k)=(y k －y ks )(y k －y ks ) T －γ(6)
其中，ks表示第s个触发时刻，满足k0＜k1＜…＜ks＜…； 若h(k)<0，则当前时刻量测为“不必要量测”且不被传输；否则，进行下一步判断；
步骤2.2. 利用最近的q个时刻的量测y k k-q 以及噪声和野值的幅值上界，设计野值检测函数f(k)以及野值检测阈值ρ，初步判断野值；
首先，利用传递函数表示公式(4)和(5)中模型的输入输出关系，如公式(7)所示；
y(z)=G yω (z)ω(z)+ G yv (z)v(z)+ G ym (z)m(z) (7)
式中，y(z)、ω(z)、v(z)和m(z)分别对应为y k 、ω k 、v k 和m k 的z变换；
G yω (z)=C(zI－F)-1，G yv (z)=I，G ym (z)=I；
有理分式矩阵G yω (z)的左矩阵分式描述为：
G yω (z)=O L (z)-1 E L (z) (8)
其中，O L (z)=z q I+ z q-1 O L,1 +…+ z 0 O L,q ，E L (z)=z q-1 E L,1+ z q-2 E L,2 +…+ z 0 E L,q ；实常数矩阵O L,μ 和E L,μ 分别为多项式矩阵O L (z)和E L (z)每一项的系数，其中，μ=1,2,…,q；
令M L =[I O L,1 … O L,q ]，N L =[ E L,1 E L,2 … E L,q ]，y k k-q =[ y k y k-1 … y k-q ]；
然后，构造野值检测函数f(k)和野值检测阈值ρ，分别如公式(9)和(10)所示；
f(k)=|| M L y k k-q || (9)
ρ=|| N L || qε+|| M L ||( q+1)υ(10)
若f(k)<ρ，当前时刻量测值传输至滤波器；否则，进行下一步判断；
步骤2.3. 根据野值发生的间歇性特点，通过判断当前时刻量测是否满足k > t(j－1)－η且j≥1的条件，从而确定当前时刻量测是否为野值，具体为：
若k≤t(j－1)－η，则当前时刻量测值传输至滤波器；
若k>t(j－1)－η，则当前时刻量测值判断为野值，且不被传输；
步骤2.4. 触发探测器将探测结果传输至滤波器；
根据零阶保持策略，得到k时刻滤波器接收的量测值；
根据步骤2.1至步骤2.3，当满足触发函数h(k)<0时，则当前时刻量测为“不必要量测”，当同时满足f(k)>ρ和k>t(j－1)－η，当前时刻量测为“野值”；
在上述情况下，当前时刻量测不被传输，触发探测器参数d k =0，根据零阶保持策略，利用最新时刻传输值y ks 来代替；否则包含有效信息的当前时刻量测值传输到滤波器，d k =1；
因此，k时刻滤波器接收的量测值为z k ，z k =d k y k +(1－d k )y ks (11)。
4.根据权利要求3所述的野值影响下基于事件触发机制的谐波检测方法，其特征在于，
所述步骤3具体为：
步骤3.1. 根据谐波检测模型和融合野值检测的事件触发机制，滤波器结构如下：
x * k+1|k =F x * k|k (12)
x * k+1|k+1= x * k+1|k +(1－d k )L k+1(z k+1 －y * k+1|k ) +d k K k+1(y k+1 －y * k+1|k ) (13)
其中，x * k+1|k 和x * k+1|k+1分别表示k+1时刻的预测值和估计值，x * k|k 表示k时刻的估计值，y * k+1|k 为k+1时刻的量测预测值，L k+1和K k+1分别为不触发和触发情况下的滤波器增益；
步骤3.2. 计算预测误差协方差P k+1|k 和估计误差协方差矩阵P k+1|k+1；
P k+1|k =FP k|k + Q k －d k (FK k S k T + S k K k T F T ) (14)
P k+1|k+1= d k [(I－K k+1 C) P k+1|k (I－K k+1 C) T + K k+1 R k+1 K k+1 T ]+(1－d k )[(1+β 1)×(I－L k+ 1 C) P k+1|k (I－L k+1 C) T + (1+β 2) L k+1 R k+1 L k+1 T +(1+β 1 -1+β 2 -1) ×γL k+1 L k+1 T ] (15)
其中，P k|k 为k时刻的估计误差协方差矩阵，R k+1为k+1时刻的量测噪声协方差矩阵，K k 表示k时刻触发情况下的滤波器增益，β 1和β 2为给定的正标量；
估计误差协方差上界P * k+1|k+1为：P * k+1|k+1=(1+β 1) (I－L k+1 C) P k+1|k (I－L k+1 C) T +(1+β 2)L k+1 R k+1 L k+1 T +(1+β 1 -1+β 2 -1)γL k+1 L k+1 T (16)
步骤3.3. 通过最小化估计误差协方差矩阵或其上界的迹，计算滤波器增益；
当d k =1，K k+1= P k+1|k C T (C P k+1|k C T + R k+1)-1；
当d k =0，L k+1=(1+β 1) P k+1|k C T ((1+β 1) CP k+1|k C T +(1+β 2) R k+1+(1+β 1 -1+β 2 -1)γI) -1；
步骤3.4. 计算基波和第i次谐波的幅值A * i,k 和相角估计值φ * i,k ；
Figure 739026DEST_PATH_IMAGE008
；
其中，x * i,k|k 和x * i+1,k|k 分别表示第i维和第i+1维状态估计值；
令
Figure 673484DEST_PATH_IMAGE009
；
则相角估计值φ * i,k =arctan(ψ i,k /ψ i+1,k )。
CN202310022658.1A 2023-01-06 2023-01-06 野值影响下基于事件触发机制的谐波检测方法 Active CN115906535B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310022658.1A CN115906535B (zh) 2023-01-06 2023-01-06 野值影响下基于事件触发机制的谐波检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310022658.1A CN115906535B (zh) 2023-01-06 2023-01-06 野值影响下基于事件触发机制的谐波检测方法

Publications (2)

Publication Number Publication Date
CN115906535A CN115906535A (zh) 2023-04-04
CN115906535B true CN115906535B (zh) 2023-05-23

Family

ID=85733639

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310022658.1A Active CN115906535B (zh) 2023-01-06 2023-01-06 野值影响下基于事件触发机制的谐波检测方法

Country Status (1)

Country Link
CN (1) CN115906535B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116527060B (zh) * 2023-05-29 2024-01-05 北京理工大学 基于事件触发采样的信息压缩与异常检测方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104485669A (zh) * 2014-12-31 2015-04-01 广州航海学院 一种挖泥船电网谐波监测方法与系统
CN105425039A (zh) * 2015-12-29 2016-03-23 南京因泰莱电器股份有限公司 基于自适应卡尔曼滤波的谐波检测方法
WO2018014602A1 (zh) * 2016-07-19 2018-01-25 东南大学 适于高维gnss/ins深耦合的容积卡尔曼滤波方法
CN109871591A (zh) * 2019-01-24 2019-06-11 武汉大学 一种igbt功率模块在线估算结温的方法

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1600947A3 (en) * 2004-05-26 2005-12-21 Honda Research Institute Europe GmbH Subtractive cancellation of harmonic noise
US7323673B1 (en) * 2005-03-16 2008-01-29 Apache Technologies, Inc. Modulated laser light detector with discrete fourier transform algorithm
CN106707348A (zh) * 2016-12-09 2017-05-24 济南赛英立德电子科技有限公司 基于时频分析的微波多普勒静态人体探测方法及探测器
US9835433B1 (en) * 2017-05-09 2017-12-05 Tesa Sa Touch trigger probe
CN109239403B (zh) * 2018-10-17 2020-10-27 西北工业大学 一种基于时间测量的单器件虚拟加速度计及其实现方法
CN111245099B (zh) * 2020-03-19 2021-03-12 山东科技大学 基于事件触发传输机制和混合量测的配电网状态估计方法
CN112180899B (zh) * 2020-09-30 2021-08-24 山东科技大学 一种间歇异常测量检测下系统的状态估计方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104485669A (zh) * 2014-12-31 2015-04-01 广州航海学院 一种挖泥船电网谐波监测方法与系统
CN105425039A (zh) * 2015-12-29 2016-03-23 南京因泰莱电器股份有限公司 基于自适应卡尔曼滤波的谐波检测方法
WO2018014602A1 (zh) * 2016-07-19 2018-01-25 东南大学 适于高维gnss/ins深耦合的容积卡尔曼滤波方法
CN109871591A (zh) * 2019-01-24 2019-06-11 武汉大学 一种igbt功率模块在线估算结温的方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
SOC estimation of EV batteries based on Temperature-Compensation model and Robust Extended Kalman Filtering;Xingzhen bai;IEEE;全文 *

Also Published As

Publication number Publication date
CN115906535A (zh) 2023-04-04

Similar Documents

Publication Publication Date Title
Peng et al. Polynomial chirplet transform with application to instantaneous frequency estimation
CN106845010B (zh) 基于改进SVD降噪和Prony的低频振荡主导模式辨识方法
CN115906535B (zh) 野值影响下基于事件触发机制的谐波检测方法
CN106597408B (zh) 基于时频分析和瞬时频率曲线拟合的高阶pps信号参数估计方法
CN104360251B (zh) 一种变压器局部放电的超声波信号时延估计方法
CN110224394B (zh) 适用于非平稳功率振荡信号特征提取的傅里叶分解算法
Ma et al. Centralized fusion estimators for multi-sensor systems with multiplicative noises and missing measurements
CN104833852A (zh) 一种基于人工神经网络的电力系统谐波信号估计测量方法
CN104215833B (zh) 电力系统频率测量方法及装置
CN113158907A (zh) 基于小波和混沌理论的微弱舰船辐射特征信号检测方法
CN112670990A (zh) 基于MEEMD-Prony联合算法的电力系统低频振荡特征参数提取方法
CN109067433B (zh) 适用于智能电能表的低压电力线载波通信噪声抑制方法
Ibrahim Gabr et al. Hybrid detection algorithm for online faulty sensors identification in wireless sensor networks
Kolosok et al. Wavelet analysis of PMU measurements for identification of cyber attacks on TCMS
CN112014811B (zh) 一种雷达载波频率的精细估计方法
CN111579868B (zh) 一种高次谐波的测量方法及装置
Overbey et al. Analysis of local state space models for feature extraction in structural health monitoring
Sarmah et al. Influence of system size on spatiotemporal dynamics of a model for plastic instability: Projecting low-dimensional and extensive chaos
Gao et al. Health monitoring of controller area network in hybrid excavator based on the message response time
Xie et al. Fault Diagnosis of Industrial Robots Based on Phase Difference Correction Method
Romulus A comparison between instantaneous frequency estimation methods of frequency modulated signals covered with gaussian noise
CN111505420B (zh) 一种线路避雷器状态的在线监测与诊断方法及系统
Dardanelli et al. Model-based Kalman filtering approaches for frequency tracking
CN116735957B (zh) 计及主瓣重叠干扰的近频谐波与间谐波测量方法及系统
CN117309079B (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