CN103630742B - 一种动态信号参数的获取方法 - Google Patents

Abstract

本申请公开了一种动态信号参数的获取方法，包括：选取电网的动态采样信号序列组成自相关矩阵；确定自相关矩阵的有效秩和动态采样信号序列的频率分量数；建立AR模型，求解所述AR模型的模型参数；利用Prony算法，确定动态采样信号表达式及复序列，动态采样信号序列由所述复序列在满足平方误差最小条件下表示；将所述模型参数对应的特征多项式根带入所述复序列，求解所述动态采样信号序列的各项参数。本申请并不直接求解Prony算法中的参数，而是借助AR参数模型思路，将当前时刻信号看成由以前各时刻信号的线性组合形成，将非线性问题转换为线性估计问题，使得计算过程更加简单且计算结果更加精确。

Description

一种动态信号参数的获取方法

技术领域

本申请涉及电网谐波分析技术领域，更具体地说，涉及一种动态信号参数的获取方法。

背景技术

由于电力电子等非线性设备在电力系统中的广泛应用，不仅导致谐波和间谐波日益增多，而且还存在具有衰减振荡分量，严重影响了电力系统的安全运行。分析谐波、间谐波以及衰减振荡参数对电力系统有重要意义。

目前的谐波分析主要采用傅氏方法，将信号看成由一系列不衰减的正弦频率成分组成，因而无法给出动态信号中衰减振荡参数，同时在傅氏分析中的频谱泄漏和栅栏效应也会造成无法检测频率相近的接谐波的问题。自回归AR（Auto Regressive）参数谱估计方法通过建立参数模型逼近真实过程，较大提高了频率分辨率，可用于间谐波的频率分析中，但其无法给出谐波的幅值和相位参数。普罗尼Prony算法将动态信号看成由一系列具有任意幅值、相位、频率和衰减因子的衰减正弦分量组成，因而特别适合具有衰减振荡分量的非平稳过程研究。同时由于采用参数模型而克服了傅氏分析中频率分辨率受限于窗长的缺陷，因而还可用于间谐波检测中。但是，直接求解Prony算法中的幅值、相位频率和衰减因子参数将导致求解一个非线性最小二乘问题，其难度大且数值稳定性差。

因此，急需一种电网谐波分析中动态信号参数的获取方案，能快速准确的获取电网谐波中的动态信号参数。

发明内容

有鉴于此，本申请提供了一种动态信号参数的获取方法，用于快速准确的获取电网谐波中的动态信号参数。

为了实现上述目的，现提出的方案如下：

一种动态信号参数的获取方法，包括：

选取电网的动态采样信号序列，由所述动态采样信号序列组成自相关矩阵；

确定所述自相关矩阵的有效秩，根据所述有效秩确定所述动态采样信号序列的频率分量数；

建立AR模型，求解所述AR模型的模型参数；

利用Prony算法，将所述动态采样信号序列表示为一组衰减振荡的正弦分量；

确定所述动态采样信号序列的复序列，所述动态采样信号序列由所述复序列在满足平方误差最小条件下表示；

将所述模型参数对应的特征多项式根带入所述复序列，求解所述动态采样信号序列的各项参数，所述各项参数包括幅值、相位、衰减和频率。

优选地，所述自相关矩阵的阶数P_e满足以下公式：N/4<p_e<N/3，其中N为采样点的个数。

优选地，所述确定所述自相关矩阵的有效秩，根据所述有效秩确定所述动态采样信号序列的频率分量数具体为：

利用SVD方法对所述自相关矩阵进行分解：

将所述自相关矩阵分解为：R_e=USV^T，其中R_e代表所述自相关矩阵，U是p_e×p_e维正交矩阵，V是(p_e+1)×(p_e+1)维正交矩阵，S是p_e×(p_e+1)维非负对角阵；

取所述对角阵S的前p个奇异值构成的对角阵Σ_p作R_e的最佳逼近 ${\hat{R}}_e={\mathrm{U\&Sigma;}}_p{V}^T=U\left[\begin{array}{cc}{S}_p& 0\\ 0& 0\end{array}\right]V^T,$ 其中S_p=diag(σ₁,σ₂,...,σ_p)；

判断所述动态采样信号序列是否包含噪声；

若所述动态采样信号序列不包含噪声，则计算β_i=σ_i+1/σ_i,1≤i≤p_e-1，将β_i取得最大值时对应的i确定为有效秩P，将P/2的整数部分确定为所述频率分量数P′；

若所述动态采样信号序列包含噪声，则根据信噪比和β_i的局部最大值确定有效秩P，将P/2的整数部分确定为所述频率分量数P′。

优选地，所述建立AR模型的过程具体为：

将所述动态采样信号序列表示为：其中C为模型阶数，w(n)是零均值白噪声序列，a_k为C阶AR模型的模型参数。

优选地，所述求解所述AR模型的模型参数的过程具体为：

判断所述动态采样信号序列是否包含噪声；

若所述动态采样信号序列不包含噪声，取AR模型的阶数C为所述有效秩P；

若所述动态采样信号序列包含噪声，取AR模型的阶数C为所述自相关矩阵的阶数P_e；

使用协方差算法，求出所述模型参数a_k。

优选地，所述利用Prony算法，将所述动态采样信号序列表示为一组衰减振荡的正弦分量具体为：

将所述动态采样信号序列表示为：

式中T_s为采样周期，q为谐波次数。

优选地，所述确定所述动态采样信号序列的复序列具体为：

将所述复序列表示为：式中b_m=A_mexp(jθ_m)，z_m=exp[(α_m+j2πf_m)T_s]，A_m,θ_m,α_m,f_m对应为幅值、相位、衰减和频率参数。

优选地，所述平方误差最小条件具体为：

$\min [\mathrm{\&epsiv;}=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{|x\left(n\right)-\hat{x}\left(n\right)|}^2].$

优选地，所述将所述模型参数对应的特征多项式根带入所述复序列，求解所述动态采样信号序列的各项参数具体为：

由所述模型参数a_k构成特征多项式，求解其特征多项式的根z_k，z_k即对应为所述复序列表达式中的z_m；

将z_m带入所述复序列表达式中，利用最小二乘法确定参数b_m；

求解各项参数的表达式为： $\left\{\begin{array}{c}{A}_m=\left|b_m\right|\\ {\mathrm{\&theta;}}_m={\tan }^{-1}[\mathrm{Im}\left(b_m\right)/\mathrm{Re}\left(b_m\right)]\\ {\mathrm{\&alpha;}}_m=\ln \left|z_m\right|/T_s\\ f_m={\tan }^{-1}[\mathrm{Im}\left(z_m\right)/\mathrm{Re}\left(z_m\right)]/2\mathrm{\&pi;}{T}_s\end{array}\right..$

优选地，在所述求解所述动态采样信号序列的各项参数之后还包括：

根据求解的结果，判断频点数是否等于所述频率分量数P′，若是则结束，否则选取幅值较大的前P′个分量。

从上述的技术方案可以看出，本申请公开的电网动态信号参数的获取方法，首先确定动态信号的频率分量数，再利用AR方法确定动态信号的模型参数，最后使用Prony算法来求解动态信号中的频率、幅值、相位和衰减参数。本申请并不直接求解Prony算法中的参数，而是借助AR参数模型思路，将当前时刻信号看成由以前各时刻信号的线性组合形成，将非线性问题转换为线性估计问题，使得计算过程更加简单且计算结果更加精确。

附图说明

为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。

图1为本申请实施例公开的一种动态信号参数的获取方法流程图；

图2为本申请实施例公开的一种确定动态信号的频率分量数的方法流程图；

图3为本申请实施例公开的一种确定动态信号频率分量数及动态信号AR模型参数的方法流程图；

图4为本申请实施例公开的另一种动态信号参数的获取方法流程图；

图5为本申请实施例公开的又一种动态信号参数的获取方法流程图。

具体实施方式

下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其它实施例，都属于本申请保护的范围。

参见图1

图1为本申请实施例公开的一种动态信号参数的获取方法流程图。

如图1所示，该方法包括：

步骤101：选取电网的动态采样信号序列，由所述动态采样信号序列组成自相关矩阵；

具体地，选取待分析的采样信号序列x(n)，其采样点数为N，选定模型阶数为P_e，且满足N/4<p_e<N/3，模型阶数P_e在范围内任意选取整数。自相关矩阵R_e表示为：

$R_e=\left[\begin{array}{cccc}r\left(\mathrm{1,0}\right)& r\left(\mathrm{1,1}\right)& ...& r(1,p_e)\\ r\left(\mathrm{2,0}\right)& r\left(\mathrm{2,1}\right)& ...& r(2,p_e)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ r(p_e,0)& r(p_e,1)& & r(p_e,p_e)\end{array}\right]---\left(1\right)$

各元素r(i,j)定义为：

$r(i,j)=\underset{n=p_e}{\overset{N-1}{\mathrm{\&Sigma;}}}x(n-j)x(n-i),i,j=\mathrm{0,1},...,p_e---\left(2\right)$

步骤102：确定所述自相关矩阵的有效秩，根据所述有效秩确定所述动态采样信号序列的频率分量数；

具体地，计算上述公式（1）矩阵的有效秩P，然后根据有效秩即可确定动态信号的频率分量数。

步骤103：建立AR模型，求解所述AR模型的模型参数；

具体地，AR模型假定信号x(n)是由零均值白噪声序列w(n)激励一个全极点的线性时不变离散时间系统得到的，即：

其中C为模型阶数，w(n)是零均值白噪声序列，a_k为C阶AR模型的模型参数。求解AR模型的模型参数。

步骤104：利用Prony算法，将所述动态采样信号序列表示为一组衰减振荡的正弦分量；

具体地，将所述动态采样信号序列表示为：

式中T_s为采样周期，q为谐波次数。

步骤105：确定所述动态采样信号序列的复序列，所述动态采样信号序列由所述复序列在满足平方误差最小条件下表示；

步骤106：将所述模型参数对应的特征多项式根带入所述复序列，求解所述动态采样信号序列的各项参数，所述各项参数包括幅值、相位、衰减和频率。

本实施例公开的电网动态信号参数的获取方法，首先确定动态信号的频率分量数，再利用AR方法确定动态信号的模型参数，最后使用Prony算法来求解动态信号中的频率、幅值、相位和衰减参数。本申请并不直接求解Prony算法中的参数，而是借助AR参数模型思路，将当前时刻信号看成由以前各时刻信号的线性组合形成，将非线性问题转换为线性估计问题，使得计算过程更加简单且计算结果更加精确。

实施例二

本实施例中，我们详细介绍确定动态信号频率分量数的过程。

通过实施例一我们已经确定了自相关矩阵R_e，接下来可以应用SVD算法确定矩阵R_e的有效秩P，进而通过有效秩P来确定动态信号频率分量数，具体地是将自相关矩阵R_e分解为：

R_e=USV^T         （3）

其中R_e代表所述自相关矩阵，U是p_e×p_e维正交矩阵，V是(p_e+1)×(p_e+1)维正交矩阵，S是p_e×(p_e+1)维非负对角阵，其对角线上元素σ_kk即为矩阵R_e的奇异值，且满足可见，矩阵R_e较大的奇异值集中在对角阵S的前段，因此可以取其前P个奇异值构成的对角阵Σ_p作R_e的最佳逼近

${\hat{R}}_e={\mathrm{U\&Sigma;}}_p{V}^T=U\left[\begin{array}{cc}{S}_p& 0\\ 0& 0\end{array}\right]V^T,$ 其中S_p=diag(σ₁,σ₂,...,σ_p)   （4）

有效秩P及频率分量数的确定过程：

判断动态采样信号序列是否包含噪声，若信号x(n)不含噪声，则S的前P个奇异值显著大于后面的奇异值，可以计算β_i=σ_i+1/σ_i,1≤i≤p_e-1，将β_i取得最大值时对应的i确定为有效秩P，信号的频率分量数P′为P/2的整数部分；若信号x(n)包含噪声，则可以根据信噪比并结合β_i的局部最大值确定有效秩P，信号的频率分量数P′为P/2的整数部分。

参见图2，图2为本申请实施例公开的一种确定动态信号的频率分量数的方法流程图。

上述过程概括起来如下：

步骤201：接收动态信号，并组成自相关矩阵；

步骤202：利用SVD分解所述自相关矩阵；

步骤203：判断动态信号是否包含噪声；

步骤204：若判断结果为否，则计算β_i=σ_i+1/σ_i，将β_i最大时对应的i确定为有效秩P，将P/2的整数部分确定为频率分量数；

步骤205：若判断结果为是，则根据信噪比和β_i的局部最大值确定有效秩P,将P/2的整数部分确定为频率分量数。

通过上面的过程，即可确定电网动态信号的频率分量数。

且采用SVD方法即使在短采样周期也有很高的频率分辨率，可准确确定动态信号频率分量数，可较好分别出信号的间谐波分量，同时也克服了AR模型阶数选取困难问题。

实施例三

本实施例中我们详细介绍确定动态信号模型参数的过程。

建立AR模型，AR模型假定信号x(n)是由零均值白噪声序列w(n)激励一个全极点的线性时不变离散时间系统得到的，即：

$x\left(n\right)=-\underset{k=1}{\overset{C}{\mathrm{\&Sigma;}}}{a}_{k}x(n-k)+w\left(n\right)---\left(5\right)$

上式中，C为模型阶数，a_k为C阶AR模型的模型参数。

经过上一实施例确定的有效秩P，对于不含噪声的信号，取AR模型阶数为P阶；而对含噪声的信号，AR模型阶数需大幅度提高，可取为P_e阶，通过协方差算法可给出模型参数a_k为{1,a₁,a₂...a_p}或则对应的分别为AR(P)模型或AR（P_e）模型。

参见图3，图3为本申请实施例公开的一种确定动态信号频率分量数及动态信号AR模型参数的方法流程图。

在上一实施例的步骤205之后，进一步增加：

步骤206：取AR(P)模型计算a_k；

步骤207：取AR（P_e）模型计算a_k。

实施例四

本实施例中我们详细介绍动态信号的参数确定过程。

利用Prony算法把信号x(n)看成由一组衰减振荡的正弦分量组成，即：

$x\left(n\right)=\underset{i=1}{\overset{q}{\mathrm{\&Sigma;}}}{A}_i{e}^{{\mathrm{\&alpha;}}_{i}n{T}_s}\cos (2\mathrm{\&pi;}{f}_{i}n{T}_s+{\mathrm{\&theta;}}_i)---\left(6\right)$

式中，T_s为采样周期，q为谐波次数。

动态信号x(n)可以其复序列在满足平方误差最小条件下近似表示，复序列表示为：

$\hat{x}\left(n\right)=\underset{m=1}{\overset{2q}{\mathrm{\&Sigma;}}}{b}_m{z}_m^n,n=\mathrm{0,1},...,N-1---\left(7\right)$

式中，b_m=A_mexp(jθ_m)，z_m=exp[(α_m+j2πf_m)T_s]，A_m,θ_m,α_m,f_m分别对应为幅值、相位、衰减和频率参数。

平方误差最小条件为：

$\min =[\mathrm{\&epsiv;}=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{|x\left(n\right)-\hat{x}\left(n\right)|}^2]---\left(8\right)$

由的表达式可知，其为常系数线性差分方程的齐次解形式，结合实施例三中的公式（5）x(n)的差分表示形式，可知实施例三中求得的AR模型参数a_k即对应为公式（7）差分方程的系数，进而由模型参数a_k构成的特征多项式根z_k即对应为复序列表达式中的z_m。再将求得的z_m带入表达式并应用最小二乘法可确定参数b_m，最终可给出A_m,θ_m,α_m,f_m的计算公式为：

$\left\{\begin{array}{c}{A}_m=\left|b_m\right|\\ {\mathrm{\&theta;}}_m={\tan }^{-1}[\mathrm{Im}\left(b_m\right)/\mathrm{Re}\left(b_m\right)]\\ {\mathrm{\&alpha;}}_m=\ln \left|z_m\right|/T_s\\ f_m={\tan }^{-1}[\mathrm{Im}\left(z_m\right)/\mathrm{Re}\left(z_m\right)]/2\mathrm{\&pi;}{T}_s\end{array}\right.---\left(9\right)$

图4为本申请实施例公开的另一种动态信号参数的获取方法流程图。

在上一实施例的基础上，本实施例进一步增加了：

步骤208：利用Prony算法确定动态信号表达式x(n)及复序列表达式x

步骤209：计算模型参数a_k对应特征多项式的根z_k即为复序列中的z_m；

步骤210：应用最小二乘法确定复序列中的b_m；

步骤211：由z_m和b_m确定动态信号的幅值、相位、衰减和频率参数。

通过将AR方法和Prony算法结合，利用AR方法求出z_m，再利用Prony算法确定幅值、相位、衰减和频率参数，克服了AR方法仅能给出频率信息的局限，也避免了直接求解Prony模型需要解决一个非线性最小二乘问题。

实施例五

参见图5，图5为本申请实施例公开的又一种动态信号参数的获取方法流程图。

由于在确定模型参数a_k时出现了两种情况，即有噪声和无噪声的情况，在有噪声时选取的AR模型阶数P_e，而P_e显著大于频率分量数P′即P/2，所以在最后计算出的参数中，对应于有噪声的情况下，频点数肯定大于P′，所以我们增加了频点数判断过程，即步骤212：判断频点数是否等于频率分量数P′；步骤213：若判断结果为否，则选取幅值较大的前P′个分量；若判断结果为是，则结束流程。这样，即可确定出P′个参数。

实施例六

本实施例中我们将对比传统的利用Prony算法来获取动态信号的参数和本申请的获取方法。

算例一：

我们选取电网动态信号模型为：

x(t)=3cos(2π×25t+π/5)+150cos(2π×50t+π/4)+20cos(2π×150t+π/6)+2cos(2π×180t+π/3)+15cos(2π×250t+π/8)

选取采用频率Fs=2000Hz，采样时间0.04s，采样点数80个，分别考虑不加噪声、加40dB的噪声情况下传统Prony方法与本申请的方法的计算结果。参见下表1：

表1

算例二：

我们选取含间谐波和衰减分量的动态电网信号模型为;

x(t)=150e^-0.4πtcos(2πf₁t+π/3)+10e^-0.6πtcos(2πf₂t+π/4)+2e^-0.2πtcos(2πf₃t+π/5)

其中f1=50Hz，f2=148Hz，f3=245Hz，取采样频率Fs=2000Hz，采样时间0.1s，采样点数200个，分别考虑不加噪声、加40dB的噪声情况下传统Prony法与本申请方法的计算结果。参见下表2：

表2

从表1和表2的对比可见，在谐波、间谐波和衰减分量的检测中，采用本申请的参数获取方法能够大幅提高计算精度，而且对噪声有较好的适应性。

最后，还需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本申请。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本申请的精神或范围的情况下，在其它实施例中实现。因此，本申请将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (9)

1.一种动态信号参数的获取方法，其特征在于，包括：

选取电网的动态采样信号序列，由所述动态采样信号序列组成自相关矩阵；

确定所述自相关矩阵的有效秩，根据所述有效秩确定所述动态采样信号序列的频率分量数；

建立AR模型，求解所述AR模型的模型参数；

利用Prony算法，将所述动态采样信号序列表示为一组衰减振荡的正弦分量；

确定所述动态采样信号序列的复序列，所述动态采样信号序列由所述复序列在满足平方误差最小条件下表示；

将所述模型参数对应的特征多项式根带入所述复序列，求解所述动态采样信号序列的各项参数，所述各项参数包括幅值、相位、衰减和频率；

所述求解所述AR模型的模型参数的过程具体为：

判断所述动态采样信号序列是否包含噪声；

若所述动态采样信号序列不包含噪声，取AR模型的阶数C为所述有效秩P；

若所述动态采样信号序列包含噪声，取AR模型的阶数C为所述自相关矩阵的阶数P_e；

使用协方差算法，求出所述模型参数a_k。

2.根据权利要求1所述的方法，其特征在于，所述自相关矩阵的阶数P_e满足以下公式：N/4<p_e<N/3，其中N为采样点的个数。

3.根据权利要求2所述的方法，其特征在于，所述确定所述自相关矩阵的有效秩，根据所述有效秩确定所述动态采样信号序列的频率分量数具体为：

利用SVD方法对所述自相关矩阵进行分解：

将所述自相关矩阵分解为：R_e＝USV^T，其中R_e代表所述自相关矩阵，U是p_e×p_e维正交矩阵，V是(p_e+1)×(p_e+1)维正交矩阵，S是p_e×(p_e+1)维非负对角阵；

取所述对角阵S的前p个奇异值构成的对角阵Σ_p作R_e的最佳逼近 其中S_p＝diag(σ₁,σ₂,...,σ_p)；

判断所述动态采样信号序列是否包含噪声；

若所述动态采样信号序列不包含噪声，则计算β_i＝σ_i+1/σ_i,1≤i≤p_e-1，将β_i取得最大值时对应的i确定为有效秩P，将P/2的整数部分确定为所述频率分量数P′；

若所述动态采样信号序列包含噪声，则根据信噪比和β_i的局部最大值确定有效秩P，将P/2的整数部分确定为所述频率分量数P′。

4.根据权利要求3所述的方法，其特征在于，所述建立AR模型的过程具体为：

将所述动态采样信号序列表示为：其中C为模型阶数，w(n)是零均值白噪声序列，a_k为C阶AR模型的模型参数。

5.根据权利要求1所述的方法，其特征在于，所述利用Prony算法，将所述动态采样信号序列表示为一组衰减振荡的正弦分量具体为：

将所述动态采样信号序列表示为：

式中T_s为采样周期，q为谐波次数。

6.根据权利要求5所述的方法，其特征在于，所述确定所述动态采样信号序列的复序列具体为：

将所述复序列表示为：,n＝0,1,…,N-1，式中b_m＝A_mexp(j_θm)，z_m=exp[(α_m+j2πf_m)T_s], A_m,θ_m,α_m,f_m对应为幅值、相位、衰减和频率参数。

7.根据权利要求6所述的方法，其特征在于，所述平方误差最小条件具体为：

8.根据权利要求7所述的方法，其特征在于，所述将所述模型参数对应的特征多项式根带入所述复序列，求解所述动态采样信号序列的各项参数具体为：

由所述模型参数a_k构成特征多项式，求解其特征多项式的根z_k，z_k即对应为所述复序列表达式中的z_m；

将z_m带入所述复序列表达式中，利用最小二乘法确定参数b_m；

求解各项参数的表达式为：

9.根据权利要求8所述的方法，其特征在于，在所述求解所述动态采样信号序列的各项参数之后还包括：

根据求解的结果，判断频点数是否等于所述频率分量数P′，若是则结束，否则选取幅值较大的前P′个分量。