CN106849131A - 一种基于四阶混合平均累积量与改进tls‑esprit算法的低频振荡模态辨识方法 - Google Patents

Abstract

本发明涉及一种基于四阶混合平均累积量与改进TLS‑ESPRIT算法的低频振荡模态辨识方法，利用四阶混合累积量抑制高斯色噪声，并对TLS‑ESPRIT算法进行改进，引入奇异值相对变化率定阶，利用测量数据构造矩阵进行电力系统低频振荡模态辨识，本发明使得检测出电力系统低频振荡信号不同振荡模态的振荡频率、衰减因子、振荡幅值和相位信息更加准确。

Description

一种基于四阶混合平均累积量与改进TLS-ESPRIT算法的低频 振荡模态辨识方法

技术领域

本发明涉及电力系统安全运行领域，特别是一种基于四阶混合平均累积量与改进TLS-ESPRIT算法的低频振荡模态辨识方法。

背景技术

随着电力技术的不断进步，大电网互联已经逐步实现，加之快速高放大倍数的励磁装置的广泛使用，导致低频振荡现象时有发生。电力系统低频功率振荡是发生在电力系统中的一种有功功率振荡，振荡频率通常为0.2Hz-2.5Hz。低频振荡问题影响电力系统的安全与稳定，因此对低频振荡模态进行有效的快速辨识，并研究相应的低频振荡抑制策略和方法，是保证电力系统安全可靠运行的关键。

目前，在电力系统低频振荡主导模态辨识及其特征提取方面，通常使用TLS-ESPRIT、快速傅立叶变换(FFT)、小波分析、Prony分析、希尔伯特—黄变换(HHT)。FFT方法有较好的准确性和鲁棒性，但该方法不能有效提取振荡的衰减特征；小波分析存在小波基难以选取的问题且难以实现快速在线辨识；Prony方法受噪声影响较大，在低性噪比的噪声条件下，辨识精度会存在很大的误差；HHT分析适合非线性、非平稳的信号特征参数的提取，缺点是是在提取信号过程中会出现端点效应和模态混叠，并且算法的复杂度高，不适用于在线辨识；基于空间旋转不变技术的ESPRIT方法，该方法具有优于Prony的抗干扰能力和运算效率。

电力系统中的谐波污染越来越严重，越来越复杂，从现场抽取的数据中包含着大量的噪声，干扰的存在会影响对信号辨识结果的准确性。目前，电力系统模态辨识一般都是通过自带的奇异值分解进行噪声的消除，这种情况对噪声的估计并不足，尤其对低信噪比情况下的噪声处理能力有限。

发明内容

有鉴于此，本发明的目的是提出一种基于四阶混合平均累积量与改进TLS-ESPRIT算法的低频振荡模态辨识方法，有利于更加准确地对电力系统低频振荡的模态进行辨识。

本发明采用以下方案实现：一种基于四阶混合平均累积量与改进TLS-ESPRIT算法的低频振荡模态辨识方法，具体包括以下步骤：

步骤S1：提取电力系统低频振荡信号Y(t)，电力系统低频振荡信号可以表示成如下所示：

式中：2P为信号模态数的两倍；A_m为幅值；σ_m为衰减因子；f_m为频率；θ_m为初相角；v(t)为高斯白噪声；

步骤S2：低频振荡信号Y(t)通过低通滤波器产生一定高度的高斯色噪声Y₁(t)；

步骤S3：对含高斯色噪声信号Y₁(t)去均值得到Y'(t)；

步骤S4：利用Y'(t)构造四阶混合平均累积量

式中：S₁＝max(0,-τ)；S₂＝min(N-1,N-1-τ)；τ≤N-1；Y'^*为Y'的共轭复数，N为采样点数；

步骤S5：利用构造Hankel矩阵Y：

式中，L＞2P；M＞2P；P为低频振荡的模态数，N为采样点数。

步骤S6：对Y进行奇异值分解，得到由奇异值矩阵所组成的对角矩阵，如下：

Y＝USV^T；

式中，U^HU＝I；V^HV＝I；U为L×L阶方阵；V为M×M阶方阵；H表示共轭转置；奇异值S为L×M阶对角阵，其中的奇异值为s，按照降序排列；V按照奇异值大小划分为信号子空间V_S和噪声子空间V_N；其中具体根据信号的强弱来划分。本发明提出利用奇异值相对变化率来进行定阶，取最大的奇异值相对变化率的下角标为系统模型的阶数，即信号子空间V_S。详细见步骤S71，S72，S73。

步骤S7：根据奇异值相对变化率法确定系统的阶数2P，V_S的列向量对应于矩阵Y的幅值最大的2P个奇异值的特征向量；

步骤S8：V_S删除第一行和第二行剩下矩阵分别为V₁、V₂，构造矩阵|V₁V₂|并对矩阵|V₁V₂|进行奇异值分解得到得到右特征向量向量对矩阵进行分块将分为四个2P×2P的矩阵：

步骤S9：计算的特征值λ_k(k＝1,2,…2P)，根据下列式子计算各分量的振荡频率f_k和衰减系数σ_k：

其中，T_s为采样周期；

步骤S10：通过最小二乘法求得信号的振荡幅值和初相角：

步骤S11：判断步骤S10得出的振荡幅值与初相角是否满足预设的拟合条件，若是，则输出低频振荡信息，得出电力系统低频振荡各模态的完整参数；否则，返回步骤S7。

进一步地，所述步骤S7具体包括以下步骤；

步骤S71：将得到的奇异值从大到小排列：

s₁≥s₂≥s_2P≥s_2P+1≈…≈s_L≈0；

步骤S72：定义奇异值相对变化率如下所示：

步骤S73：取RCRSV的最大值所对应的下标对应为系统模态的阶数2P。

进一步地，所述步骤S10具体包括以下步骤：

步骤S101：列出下列方程：

Y＝λc；

c＝|c₁,c₂,…,c_2P|；

步骤S102：用最小二乘法得到方程的解为c＝(λ^Tλ)-¹λ^TY，可得各个分量的幅值A_k和相位θ_k分别为：

与现有技术相比，本发明有以下有益效果：

1、本发明通过构造四阶混合平均累积量对电力系统高斯色噪声进行抑制，从而使辨识结果更加准确。

2、本发明提出的定阶方法在电力系统低频振荡辨识的过程中定阶计算量较小、速度快，受主观因素影响较小。

3、本发明方法可以快速、准确地辨识出电力系统低频振荡各模态参数，具有抗噪能力强、辨识精度高等优点。

附图说明

图1为本发明实施例中的方法流程示意图。

图2为本发明实施例中测试信号、含噪信号及辨识后拟合信号。

图3为本发明实施例中奇异值分布与奇异值相对变化率分布图。

图4为本发明实施例中对构建的测试信号进行辨识得到的参数表。

具体实施方式

下面结合附图及实施例对本发明做进一步说明。

如图1所示，本实施例对Y(t)＝e^-0.3tcos(2π×0.5t)+e^-0.6tcos(2π×1.1t)测试信号进行分析，对本发明的技术方案进行具体说明。

如图1所示，本实施例提供了一种基于四阶混合平均累积量与改进TLS-ESPRIT的低频振荡模态辨识方法，具体包括以下步骤：

步骤S1：提取电力系统低频振荡信号，电力系统低频振荡信号可以表示成如下所示：

式中：2P为假设的信号模态数的两倍；A_m为幅值；σ_m为衰减因子；f_m为频率；θ_m为初相角；v(t)为高斯白噪声。

步骤S2：低频振荡信号Y(t)通过低通滤波器产生一定高度的高斯色噪声Y₁(t)。

步骤S3：对含高斯色噪声信号Y₁(t)去均值得到Y'(t)。

步骤S4：利用Y'(t)构造四阶混合平均累积量

式中：S₁＝max(0,-τ)；S₂＝min(N-1,N-1-τ)；τ≤N-1；Y'^*为Y'的共轭复数。

步骤S5：利用构造Hankel矩阵Y，式子如下所示：

式中，L＞2P；M＞2P；P为低频振荡的模态数；N为采样点数。

步骤S6：对Y进行奇异值分解，得到由奇异值矩阵所组成的对角矩阵，如下：Y＝USV^T。

式中，H表示共轭转置；U^HU＝I；V^HV＝I；U为L×L阶方阵；V为M×M阶方阵；奇异值S为L×M阶对角阵，按照降序排列；V按照奇异值大小划分为信号子空间V_S和噪声子空间V_N。；其中具体根据信号的强弱来划分。本发明提出利用奇异值相对变化率来进行定阶，取最大的奇异值相对变化率的下角标为系统模型的阶数，即信号子空间V_S。详细见步骤S71，S72，S73。

步骤S7：根据奇异值相对变化率法确定系统的阶数2P，V_S的列向量对应于矩阵Y的幅值最大的2P个奇异值的特征向量。

步骤S8：V_S删除第一行和第二行剩下矩阵分别为V₁、V₂，构造矩阵|V₁V₂|并对矩阵|V₁V₂|进行奇异值分解得到得到右特征向量向量对矩阵进行分块将分为四个2P×2P的矩阵为

步骤S9：计算的特征值λ_k(k＝1,2,…2P)，根据下列式子计算各分量的振荡频率和衰减系数，T_s为采样周期。

步骤S10：进一步通过最小二乘法求得振荡幅值和初相角，过程如下所示：

Y＝λc

c＝|c₁,c₂,…,c_2P|

用最小二乘法得到方程的解为c＝(λ^Tλ)-¹λ^TY，可得各个分量的幅值和相位分别为：

步骤S11：得出电力系统低频振荡各模态的完整参数。

较佳的，辨识后的模态参数如图4所示，测试信号、含噪信号及辨识后的拟合信号如图2所示。

在本实施例中，所述步骤S7具体包括以下步骤：

步骤S71：将得到的奇异值从大到小排列：

s₁≥s₂≥s_2P≥s_2P+1≈…≈s_L≈0；

步骤S72：定义奇异值相对变化率如下所示：

步骤S73：取RCRSV的最大值所对应的下标对应为系统模态的阶数2P。

其中，奇异值分布与奇异值相对变化率分布图如图3所示。从图中可以看出系统的模态阶数为4。

综上所述，本实施例可以有效的抑制高斯色噪声对振荡模态辨识过程的影响；定阶计算量较小、速度快，受主观因素影响较小；可以快速、准确地辨识出电力系统低频振荡各模态参数。

以上所述仅为本发明的较佳实施例，凡依本发明申请专利范围所做的均等变化与修饰，皆应属本发明的涵盖范围。

Claims (3)

1.一种基于四阶混合平均累积量与改进TLS-ESPRIT算法的低频振荡模态辨识方法，其特征在于：包括以下步骤：

步骤S1：提取电力系统低频振荡信号Y(t)，电力系统低频振荡信号可以表示成如下所示：

$Y\left(t\right)=x\left(t\right)+v\left(t\right)=\underset{m=1}{\overset{2P}{\Σ}}{A}_m{e}^{{\mathrm{\σ}}_{m}t}cos(2{\mathrm{\πf}}_{m}t+{\mathrm{\θ}}_m)+v\left(t\right)$

式中：2P为信号模态数的两倍；A_m为幅值；σ_m为衰减因子；f_m为频率；θ_m为初相角；v(t)为高斯白噪声；

步骤S2：低频振荡信号Y(t)通过低通滤波器产生一定高度的高斯色噪声Y₁(t)；

步骤S3：对含高斯色噪声信号Y₁(t)去均值得到Y'(t)；

步骤S4：利用Y'(t)构造四阶混合平均累积量

${\stackrel{\‾}{C}}_{4Y}\left(\mathrm{\τ}\right)=-\frac{2}{N}\underset{t=S_1}{\overset{S_2}{\Σ}}{Y}^{\′}\left(t\right)Y^{\′}(t+\mathrm{\τ}){\left(Y^{\′*}\right(t+\mathrm{\τ}\left)\right)}^2;$

式中：S₁＝max(0,-τ)；S₂＝min(N-1,N-1-τ)；τ≤N-1；Y'^*为Y'的共轭复数，N为采样点数；

步骤S5：利用构造Hankel矩阵Y：

式中，L＞2P；M＞2P；P为低频振荡的模态数，N为采样点数。

步骤S6：对Y进行奇异值分解，得到由奇异值矩阵所组成的对角矩阵，如下：

Y＝USV^T；

式中，U^HU＝I；V^HV＝I；U为L×L阶方阵；V为M×M阶方阵；H表示共轭转置；奇异值S为L×M阶对角阵，其中的奇异值为s，按照降序排列；V按照奇异值大小划分为信号子空间V_S和噪声子空间V_N；

步骤S7：根据奇异值相对变化率法确定系统的阶数2P，V_S的列向量对应于矩阵Y的幅值最大的2P个奇异值的特征向量；

步骤S8：V_S删除第一行和第二行剩下矩阵分别为V₁、V₂，构造矩阵|V₁V₂|并对矩阵|V₁V₂|进行奇异值分解得到得到右特征向量向量对矩阵进行分块将分为四个2P×2P的矩阵：

$\stackrel{\‾}{R}=\left|\begin{array}{cc}{R}_{11}& R_{12}\\ R_{21}& R_{22}\end{array}\right|;$

步骤S9：计算的特征值λ_k(k＝1,2,…2P)，根据下列式子计算各分量的振荡频率f_k和衰减系数σ_k：

$\left\{\begin{array}{c}{f}_k=\frac{{\mathrm{arg\λ}}_k}{2{\mathrm{\πT}}_s}\\ {\mathrm{\σ}}_k=-\frac{ln\left|{\mathrm{\λ}}_k\right|}{T_s}\end{array}\right.;$

其中，T_s为采样周期；

步骤S10：通过最小二乘法求得信号的振荡幅值和初相角：

步骤S11：判断步骤S10得出的振荡幅值与初相角是否满足预设的拟合条件，若是，则输出低频振荡信息，得出电力系统低频振荡各模态的完整参数；否则，返回步骤S7。

2.根据权利要求1所述的一种基于四阶混合平均累积量与改进TLS-ESPRIT算法的低频振荡模态辨识方法，其特征在于：所述步骤S7具体包括以下步骤；

步骤S71：将得到的奇异值从大到小排列：

s₁≥s₂≥s_2P≥s_2P+1≈…≈s_L≈0；

步骤S72：定义奇异值相对变化率如下所示：

${\mathrm{RCRSV}}_i=\frac{s_i-s_{i+1}}{s_{i+1}},i=1,2,...L;$

步骤S73：取RCRSV的最大值所对应的下标对应为系统模态的阶数2P。

3.根据权利要求1所述的一种基于四阶混合平均累积量与改进TLS-ESPRIT算法的低频振荡模态辨识方法，其特征在于：所述步骤S10具体包括以下步骤：

步骤S101：列出下列方程：

Y＝λc；

c＝|c₁,c₂,…,c_2P|；

$\mathrm{\λ}=\left(\begin{array}{cccc}1& 1& \mathrm{...}& 1\\ {\mathrm{\λ}}_1& {\mathrm{\λ}}_2& \mathrm{...}& {\mathrm{\λ}}_{2P}\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ {\mathrm{\λ}}_1^{N-1}& {\mathrm{\λ}}_2^{N-1}& \mathrm{...}& {\mathrm{\λ}}_{2P}^{N-1}\end{array}\right);$

步骤S102：用最小二乘法得到方程的解为c＝(λ^Tλ)-¹λ^TY，可得各个分量的幅值A_k和相位θ_k分别为：

$\left\{\begin{array}{c}{A}_k=2\left|c_k\right|\\ {\mathrm{\θ}}_k=\arg \left(c_k\right)\end{array}\right..$