CN105468880B - 一种低频振荡快速衰减信号分量参数的提取方法 - Google Patents

Abstract

本发明公开了一种Prony法低频振荡快速衰减信号分量参数的提取方法，包括以下步骤：读入滤波后的测量数据；划分时段；在第二时段内分析慢速衰减信号分量参数；调整慢速衰减信号分量参数；在第一时段内分析快速衰减信号分量参数；分析结果输出。本发明把采样时段划分成两个时段，第二时段求低频振荡中慢速衰减分量参数，第一时段求快速衰减分量参数，可以有效提取快速衰减分量信息，同时也能提高慢速衰减分量参数的计算精度。本发明求得的低频振荡中各信号分量的频率、衰减因子、幅值和初相参数准确度高，可以更好反映低频振荡的特征，能为采取措施抑制低频振荡提供可靠依据。

Description

一种低频振荡快速衰减信号分量参数的提取方法

技术领域

本发明涉及一种电力系统的Prony法低频振荡的计算方法，特别是一种Prony法低频振荡快速衰减信号分量参数的提取方法。

背景技术

电力系统运行时，由于扰动会发生频率为0.1Hz～2.5Hz的低频振荡，严重影响电力系统的稳定性，危及电网及其相关设备的安全运行。因此必须严密监测可能出现的低频振荡现象，及时处理，防止造成破坏系统稳定性的后果。Prony分析是目前分析低频振荡最常用最有效的方法之一，通过对信号进行Prony分析，可以直接得到反映低频振荡的幅值、相位、频率及衰减因子等参数。

Prony算法用一系列的任意频率、衰减因子、幅值和初相的指数函数的线性组合来拟合一个函数，不用再通过频域响应来求解，其计算量大为减少。也就是说该数学模型可以由一组衰减的正弦分量组成，是一种使用线性方程组来求解非线性问题的分析方法。

假设输入信号x(n)有N个采样点，即输入信号为x(0),…,x(N-1)。对该输入信号建立信号模型如下：

式中，是x(n)的近似值，p为信号模型的阶数，p值与信号分量数有关，b_i和z_i的表达式为：

式中，A_i为幅值，θ_i为初相，α_i为衰减因子，f_i为频率，Δt为采样时间间隔。

为了求式(1)中z_i的值，由z_i构造如下的特征多项式

则z_i是特征多项式(4)构成的如下特征方程的根

式中，a₀＝1。

根据式(1)和式(5)，可以推导出满足递推的常系数线性差分方程，该差分方程为：

式(6)中是实际测量数据x(n)的近似值，它们之间存在误差e(n)，即

式(6)代入式(7)，得

对式(8)进行最小二乘估计，使得误差平方和最小，则导致一组非线性方程，求解困难。为了实现线性估计，令

则式(8)可写成

式(10)写成矩阵形式为

式(11)为特征方程系数求解方程组，它的方程个数为N–p，未知数个数为p，一般情况下方程个数多于未知数个数，可以用最小二乘法计算估计值。求解式(11)的线性最小二乘法称为扩展Prony法。

采用最小二乘法，使最小，得到最小二乘法的法方程

式中，样本函数r(i,j)为(下式中上标(*)表示共轭)

最小误差能量为

用高斯消去法解式(12)，就可以求出a＝[a₁ a₂ … a_p]^T。a的精确求解是Prony法的关键，求出的a代入特征方程(5)，求出z＝[z₁ z₂ … z_p]^T后，就可以计算出信号分量的频率和衰减因子。式(5)是一元高次方程，一般采用QR分解法求解。

信号分量的频率和衰减因子计算公式如下：

式中，arctan、ln、Im、Re分别为反正切函数、自然对数函数、取复数虚部函数、取复数实部函数。

计算出信号模型的z_i代入式(1)，得到关于未知数b的线性方程如下：

式中

b＝[b₁ b₂ … b_p]^T (18)

式(17)的Z是N×p维的范德蒙矩阵。由于z_i各不相同，范德蒙矩阵Z的各列线性独立，为满列秩的。求解式(16)计算出b后，就可以计算信号分量的幅值和初相了。式(16)也是方程个数多于未知数个数的方程组，仍然用最小二乘法计算估计值。

幅值和初相计算公式为

由于事先不知道低频振荡中信号分量的个数，即式(11)中的p不确定，因此需要把信号模型的阶数取得大些(即取为p_e>p)，再通过奇异值分解法计算出方程组系数矩阵的有效秩r，即可得到信号模型的实际阶数p＝r。这样把式(11)改写为式(21)

式(21)的系数矩阵为扩展阶矩阵

扩展阶矩阵X_e为(N-p_e)×p_e阶矩阵，用奇异值分解法可以求出该扩展阶矩阵的有效秩r，即可得到信号模型的实际阶数p＝r。

如图1和图2所示，现有电力系统的Prony法低频振荡的分析方法，主要包括以下步骤：

A、读入滤波后的测量数据

测量数据中含有高频信号及噪声，需先对测量数据滤波去掉高频信号及噪声，才能进行Prony法低频振荡分析。

B、奇异值分解法信号模型实际阶数计算

根据测量数据形成扩展阶矩阵X_e，所取信号模型阶数p_e大于信号模型实际阶数p(p_e取满足的整数，符号为向下取整符号)，然后采用奇异值分解法对扩展阶矩阵X_e进行计算，求出扩展阶矩阵的所有奇异值，此奇异值为非负的，并按下列顺序排列：

σ₁₁≥σ₂₂≥…≥σ_hh≥0 (23)

式中

h＝min(N-p_e,p_e) (24)

求出所有奇异值并排顺后，按下式对奇异值归一化

σ_kk0＝σ_kk/σ₁₁,1≤k≤h (25)

选择一个较小的正数作为阈值，使σ_kk0大于此阈值的最大整数k取为扩展阶矩阵X_e的有效秩r，则信号模型的实际阶数p＝r。

C、特征方程系数计算

计算时，一般需要把信号模型的阶数m取得比信号模型的实际阶数p大一些，否则计算结果的精度很差。取方程未知数个数m大于信号模型实际阶数p(m取满足的整数，符号为向下取整符号)，重新列方程(11)得到新的特征方程系数求解方程组如下：

此方程组的方程数(N–m)大于未知数个数m，可以采用最小二乘法对未知数进行估计，计算出对应特征方程的系数。由于方程(26)是线性方程，形成相应法方程后采用高斯消去法求解。

D、频率和衰减因子计算

信号模型的阶数为m时的特征方程为：

式中，a₀＝1。

式(27)是一元高次方程，其解为实数或共轭复数，一般采用QR分解法求解，得到z后，就可以计算出信号分量的频率和衰减因子。

E、幅值和初相计算

对应式(26)，求解信号模型的系数b_i的公式(17)和(18)相应变为：

b＝[b₁ b₂ … b_m]^T (29)

把求出的z代入式(28),并根据式(16)计算出b后，就可以计算信号分量的幅值和初相。式(16)的方程数也大于未知数个数，仍然可以用最小二乘法对未知数进行估计。计算出信号分量按幅值由大到小排序，取前p个信号分量作为最终的信号分量结果。

计算出的信号分量为如下复指数函数，

一般情况下，复指数函数成对出现，即复指数函数和它的共轭值成对出现。函数x′_i(t)的共轭值为

信号分量x′_i(t)和信号分量x″_i(t)叠加为正弦分量如下：

由于计算特征方程系数矩阵时，未知数个数m大于信号模型实际阶数p，计算出信号分量比实际的信号分量多。这些多出的分量幅值很小，可以根据幅值剔除这些多余的分量。

以上传统分析方法求频率和衰减因子时，如果低频振荡中含有衰减较快的分量，分析就非常困难，分析得到的衰减较快信号分量的参数误差很大，甚至分析结果中不包含某些衰减较快的信号分量，丢失这些衰减较快信号分量的信息。

发明内容

为解决现有技术存在的上述问题，本发明要提出一种Prony法低频振荡快速衰减信号分量参数的提取方法，以便解决准确提取低频振荡中快速衰减分量的问题。

为了实现上述目的，本发明提出了采用Prony法分段分析方法来解决准确提取低频振荡中快速衰减分量的问题。

本发明的技术方案如下：一种Prony法低频振荡快速衰减信号分量参数的提取方法，包括以下步骤：

A、读入滤波后的测量数据

按照采样间隔Δt和采样时长T读入滤波后的测量数据，所述的采样时长T>10s，采样间隔Δt＝0.005s。

B、划分时段

把采样时长T分成两个时段，分段时刻为t_s，第一时段时长为t₁，t₁的最佳取值范围为0.5s≤t₁≤1.4s，剩余时段划为第二时段t₂，则第二时段时长t₂＝T-t₁。

C、在第二时段内分析慢速衰减信号分量参数

第二时段测量数据用来分析低频振荡中的慢速衰减信号分量参数，所述的信号分量参数为信号分量的频率、衰减因子、幅值和初相参数；采样间隔为Δt₂，Δt₂＝2k₂Δt,k₂＝5,6,…,14，从t_s时刻开始每间隔Δt₂从读入的测量数据中取一个值，形成第二时段的测量数据，采用扩展Prony法分析低频振荡中的慢速衰减信号分量参数。在第二时段内，以t_s时刻为扩展Prony法分析的0时刻，所取测量数据的个数N₂为

式中，T为原始测量数据的采样时长，t_s为分段时刻，Δt₂为第二时段测量数据的采样间隔，为向上取整符号。

D、调整慢速衰减信号分量参数

步骤C分析得到的慢速衰减信号分量参数是以分段时刻t_s为计时起点的信号分量参数，需要转换成以0时刻为计时起点的信号分量参数，各信号分量的频率和衰减因子不变，幅值和初相为

式中，A_i和θ_i分别是以t_s时刻为计时起点的第i个信号分量的幅值和初相，A_i0和θ_i0分别是以0时刻为计时起点的第i个信号分量的幅值和初相，f_i和α_i分别为第i个信号分量的频率和衰减因子。

根据式(34)得到的θ_i0需要调整为反余弦函数的主值，即使得θ_i0在[-π,π]范围内。

E、第一时段测量数据处理

步骤D得到以0时刻为计时起点的各信号分量的频率、衰减因子、幅值和初相参数后，从第一时段测量数据减去已经分析出的这些慢速衰减信号分量值，只保留快速衰减信号分量。

F、在第一时段内分析快速衰减信号分量参数

利用第一时段已处理过的测量数据分析低频振荡中的快速衰减信号分量参数，采样间隔为Δt₁，Δt₁＝k₁Δt,k₁＝1,2，从0s开始每间隔Δt₁从已处理过的测量数据中取一个值，形成第一时段的测量数据，采用扩展Prony法分析低频振荡中的快速衰减信号分量参数。第一时段所取测量数据的个数N₁为

G、分析结果输出

将第二时段得到的慢速衰减信号分量参数与第一时段得到的快速衰减信号分量参数合在一起，作为分段Prony法得到的所有信号分量参数。

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

1、本发明把采样时段划分成两个时段，第二时段求低频振荡中慢速衰减分量参数，第一时段求快速衰减分量参数，可以有效提取快速衰减分量信息，同时也能提高慢速衰减分量参数的计算精度。

2、本发明求得的低频振荡中各信号分量的频率、衰减因子、幅值和初相参数准确度高，可以更好反映低频振荡的特征，能为采取措施抑制低频振荡提供可靠依据。

附图说明

本发明共有附图3张。其中：

图1是Prony法低频振荡分析的流程图。

图2是扩展Prony法的流程图。

图3是本发明Prony法低频振荡分析的流程图。

具体实施方式

下面结合附图对本发明进行进一步地说明。采用图1所示的传统Prony法低频振荡分析的流程图、图2所示的扩展Prony法的流程图及图3所示的本发明Prony法低频振荡分析的流程图，对下列信号进行了分析。

x(t)＝10e^-0.0311tcos(2π×0.2473t+20°)+2e^-0.2652tcos(2π×0.42t+13°)+

10e^-0.0027tcos(2π×0.423t+110°)+1e^-0.2936tcos(2π×1.0349t+60°)+ (36)

6.3e^-55.1156tcos(2π×0.9t+60°)+5.8e^-45.8788tcos(2π×2.4t+10°)

上述信号包含6个正弦分量，等效为12个复指数分量，采用式(1)对上述信号建模时，信号模型的实际阶数p为12。算例信号的参数见表1，其中前4个信号分量为慢速衰减信号，后2个为快速衰减信号。

表1算例信号的参数

分量序号	频率f(Hz)	衰减因子α	幅值A	初相θ(°)
					1	0.24730	-0.03110	10.00000	20.00000
2	0.42000	-0.26520	2.00000	13.00000
					3	0.42300	-0.00270	10.00000	110.00000
4	1.03490	-0.29360	1.00000	60.00000
					5	0.90000	-55.11560	6.30000	60.00000
6	2.40000	-45.87880	5.80000	10.00000

对上述信号每0.005s计算一个值，共计算4000个值，用来模拟测量数据。则此测量数据的采样间隔为0.005s，时长为20s，测量数据个数为4000。

传统的Prony法和本发明方法计算结果见表2，计算时，传统的Prony法的测量数据的采样间隔取0.05s，信号模型阶数m取36；本发明方法的分段时刻为1.0s，第1时段测量数据的采样间隔取0.01s，信号模型阶数m取16，第2时段测量数据的采样间隔取0.1s，信号模型阶数m取36。

表2两种方法的计算结果比较

由表2参照表1中各信号分量的实际参数值可见，对于低频振荡中的慢速衰减信号分量(分量1～分量4)，传统方法和本发明算法都能很好地辨识，所得正弦分量误差很小；对于快速衰减信号分量(分量5和分量6)，本发明算法能够较好地辨识，所得正弦分量误差较小，传统方法则未辨识出分量5，虽然能辨识出分量6，但误差很大。

本算法可以采用任何一种编程语言和编程环境实现，如C语言、C++、FORTRAN、Delphi等。开发环境可以采用Visual C++、Borland C++Builder、Visual FORTRAN等。

Claims (1)

1.一种Prony法低频振荡快速衰减信号分量参数的提取方法，其特征在于：包括以下步骤：

A、读入滤波后的测量数据

按照采样间隔Δt和采样时长T读入滤波后的测量数据，所述的采样时长T>10s，采样间隔Δt＝0.005s；

B、划分时段

把采样时长T分成两个时段，分段时刻为t_s，第一时段时长为t₁，t₁的最佳取值范围为0.5s≤t₁≤1.4s，剩余时段划为第二时段t₂，则第二时段时长t₂＝T-t₁；

C、在第二时段内分析慢速衰减信号分量参数

第二时段测量数据用来分析低频振荡中的慢速衰减信号分量参数，所述的信号分量参数为信号分量的频率、衰减因子、幅值和初相参数；采样间隔为Δt₂，Δt₂＝2k₂Δt,k₂＝5,6,…,14，从t_s时刻开始每间隔Δt₂从读入的测量数据中取一个值，形成第二时段的测量数据，采用扩展Prony法分析低频振荡中的慢速衰减信号分量参数；在第二时段内，以t_s时刻为扩展Prony法分析的0时刻，所取测量数据的个数N₂为

式中，T为原始测量数据的采样时长，t_s为分段时刻，Δt₂为第二时段测量数据的采样间隔；为向上取整符号；

D、调整慢速衰减信号分量参数

步骤C分析得到的慢速衰减信号分量参数是以分段时刻t_s为计时起点的信号分量参数，需要转换成以0时刻为计时起点的信号分量参数，各信号分量的频率和衰减因子不变，幅值和初相为

式中，A_i和θ_i分别是以t_s时刻为计时起点的第i个信号分量的幅值和初相，A_i0和θ_i0分别是以0时刻为计时起点的第i个信号分量的幅值和初相，f_i和α_i分别为第i个信号分量的频率和衰减因子；p为低频振荡中慢速衰减信号分量的个数；

根据式(2)得到的θ_i0需要调整为反余弦函数的主值，即使得θ_i0在[-π,π]范围内；

E、第一时段测量数据处理

步骤D得到以0时刻为计时起点的各信号分量的频率、衰减因子、幅值和初相参数后，从第一时段测量数据减去已经分析出的慢速衰减信号分量值，只保留快速衰减信号分量；

F、在第一时段内分析快速衰减信号分量参数

利用第一时段已处理过的测量数据分析低频振荡中的快速衰减信号分量参数，采样间隔为Δt₁，Δt₁＝k₁Δt,k₁＝1,2，从0s开始每间隔Δt₁从已处理过的测量数据中取一个值，形成第一时段的测量数据，采用扩展Prony法分析低频振荡中的快速衰减信号分量参数；第一时段所取测量数据的个数N₁为

G、分析结果输出

将第二时段得到的慢速衰减信号分量参数与第一时段得到的快速衰减信号分量参数合在一起，作为分段Prony法得到的所有信号分量参数。