CN106202977A - 一种基于盲源分离算法的低频振荡模式分析方法 - Google Patents

Abstract

本发明公开了一种基于盲源分离算法的低频振荡模式分析方法，包括以下步骤：采集电力系统的低频振荡信号作为输入，利用Takens嵌入定理将该单通道时间序列构造成多通道观测信号矩阵；利用盲源分离算法对该矩阵进行处理，分解得到不同的低频振荡模式；结合时频分析方法，求出不同的单模式信号的瞬时幅值和频率；对瞬时频率取平均值得到平均频率，对瞬时幅值的对数随时间的变化曲线进行最小二次拟合得到平均衰减系数。本发明方法将Takens嵌入定理、盲源分离、时频分析的理论结合在一起，不仅能够准确地从单通道时间序列中分离出不同的振荡模式，并且能够跟踪不同模式的瞬时频率和幅值的变化情况，还具有一定的抗噪声能力。

Description

一种基于盲源分离算法的低频振荡模式分析方法

技术领域

本发明涉及一种电力系统低频振荡信号模式的分析方法，尤其是一种基于盲源分离算法的低频振荡模式分析方法，属于电力系统稳定性分析技术领域。

背景技术

在日渐自由化的电力市场背景下，伴随着电力供应商越来越激烈的竞争、越来越大的市场压力、发电机储存电量的减少、长距离输电的需要以及传输电力线受到的物理限制等因素，电力系统的运行点越来越容易接近稳定极限，从而引发小干扰稳定性问题。低频振荡是小干扰稳定性问题中最为突出的一种，它的频率范围在0.1-2.5Hz之间。如果不对低频振荡进行有效的监测和控制，它将破坏系统的运行稳定，甚至引发区域断电或大面积停电。然而，有效的监测和控制基于对低频振荡信号进行的准确模式分析，这一分析可以帮助系统运行人员了解系统的运行情况，同时可以为电力系统稳定器参数的设置提供至关重要的信息。

如今已有许多用于低频振荡模式分析的方法。其中应用最为广泛的方法是Prony算法和HHT算法。然而，这两种方法在一定程度上都存在不足之处。例如，Prony算法结果的准确性容易受到噪声的影响，其模型阶数难以确定且不能够动态跟踪不同模式的瞬时幅值和频率的变化情况；HHT算法由于包含经验模态分解(Empirical mode decomposition,EMD)步骤，其分解结果在低频振荡模式相互接近或采样频率不够高的情况下，可能生成虚假模式或造成分解模式混叠。

发明内容

本发明的目的是针对上述现有技术的不足，提供一种基于盲源分离算法的低频振荡模式分析方法。

本发明的目的可以通过如下技术方案实现：

一种基于盲源分离算法的低频振荡模式分析方法，所述方法包括以下步骤：

1)通过装设在电力系统中的测量元件采集电力系统的低频振荡信号作为单通道时间序列；

2)根据Takens嵌入定理利用步骤1)所述的单通道时间序列构造多通道观测信号矩阵；

3)利用盲源分离算法对步骤2)所述的多通道观测信号矩阵进行分解，得到分离信号矩阵，所述分离信号矩阵包括单通道时间序列中所包含的低频振荡信号的不同模式；

4)结合时频分析理论，求出步骤3)所述的低频振荡信号的不同模式的瞬时幅值和瞬时频率；

5)对步骤4)所述的瞬时频率进行平均求解得到平均频率，对步骤4)所述的瞬时幅值的对数相对时间的变化曲线进行最小二次拟合求得平均衰减系数。

优选的，步骤1)中所述电力系统的低频振荡信号一般为任意两台发电机的相对角速度、传输线上的有功功率、节点电压或节点间相对角偏差中的一种。

优选地，所述任意两台发电机的相对角速度，通过下式进行计算：

ω_mn＝ω_m-ω_n

其中，ω_m和ω_n分别表示第m台发电机与第n台发电机的角速度，ω_mn为第m台发电机与第n台发电机的相对角速度，ω_mn即作为单通道时间序列。

优选地，步骤2)中所述的多通道观测信号矩阵的构造过程如下：

$X\left(t\right)=\left[\begin{array}{c}x\left(t\right)\\ x(t+d)\\ \·\\ \·\\ \·\\ x\left(t\right(k-1)d)\\ \·\\ \·\\ \·\\ x(t+(M-1)d)\end{array}\right]=\left[\begin{array}{c}{X}_1\\ X_2\\ \·\\ \·\\ \·\\ X_k\\ \·\\ \·\\ \·\\ X_M\end{array}\right]$

其中X(t)表示多通道观测信号矩阵，x(t)表示单通道时间序列，d表示时间延迟，X_k表示第k个通道观测信号，k＝1,2，…，M，M为x(t)的傅里叶变换频谱中的主成分个数的两倍。

优选地，步骤3)中所述的利用盲源分离算法对多通道观测信号矩阵进行分解采用的算法是二阶盲辨识(Second-order blind identification,SOBI)，其处理过程如下：

一、对多通道观测信号矩阵X(t)进行预白化处理，得到预白化矩阵Z(t)，其中Z(t)＝W_mX(t)，W_m为白化矩阵；

二、通过联合近似对角化技术(Joint approximate diagonalization,JAD)将预白化矩阵Z(t)的一组时延协方差矩阵R_z(τ)(τ∈{τ_j|j＝1,2,…,p})同时最大程度对角化，求出最大程度对角化的矩阵Ψ。其中，保证对角化程度最大的方法是令非对角线元素的和最小，即：

$min\underset{j=1}{\overset{p}{\Σ}}off\left({\mathrm{\Ψ}}^T{R}_z\right({\mathrm{\τ}}_j\left)\mathrm{\Ψ}\right)$

其中off表示矩阵Ψ^TR_z(τ_j)Ψ的非对角线元素，τ_j表示时延；

三、根据矩阵Ψ、白化矩阵W_m以及多通道观测信号矩阵X(t)求解混矩阵W和分离信号矩阵Y(t)：

W＝Ψ^TW_m

Y(t)＝WX(t)

分离信号矩阵Y(t)包含单通道时间序列中所包含的低频振荡信号的不同模式。

优选地，步骤4)中所述的低频振荡信号的不同模式y_i(t)的瞬时幅值A_i(t)和瞬时频率f_i(t)的计算方法包括Hilbert变换法、傅里叶变换法、小波变换法和基于Teager能量算子的方法。

优选地，所述利用Hilbert变换法计算低频振荡信号的不同模式y_i(t)的瞬时幅值A_i(t)和瞬时频率f_i(t)的公式如下：

$H\left(y_i\right(t\left)\right)=\frac{1}{\mathrm{\π}}\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}\frac{y_i\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}d\mathrm{\τ}$

$A_i\left(t\right)=\sqrt{y_i{\left(t\right)}^2+H{\left(y_i\right(t\left)\right)}^2}$

$f_i\left(t\right)=\frac{d\left(arctan\right(H\left(y_i\right(t))/y_i\left(t\right)))}{2\mathrm{\π}dt}$

其中H()表示Hilbert变换，i＝1,2,…N，N为源信号个数。

优选地，步骤5)中所述的求解平均频率f_iav和平均衰减系数σ_i的方法如下：

$f_{iav}=\underset{k=1}{\overset{K}{\Σ}}{f}_i\[k\]$

${\mathrm{\σ}}_i=\frac{d\left(ln\right(A_i\left(t\right)\left)\right)}{dt}$

其中，f_iav表示第i个低频振荡模式的频率平均值，K表示瞬时频率f_i在对应的数据窗口内的采样点总数，平均衰减系数σ_i可以利用最小二乘法线性拟合ln(A_i(t))-t曲线求得。

本发明与现有技术相比，具有如下优点和有益效果：

1、本发明利用Takens嵌入定理，根据单通道低频振荡时间序列建立多通道观测信号矩阵，应用盲源分离算法对其进行处理，实现不同振荡模式的分解，再利用时频分析方法求解不同模式的瞬时幅值和频率，能够准确地从单通道时间序列中分离出不同的振荡模式。

2、本发明利用Takens嵌入定理根据单通道时间序列，建立多通道观测信号矩阵供盲源分离算法分解使用，由于只需要在一个低频振荡信号上做处理，降低了数据处理量和多变量数据采样的同步性问题。

3、本发明利用盲源分离算法处理多通道观测信号矩阵，不仅可以准确地分离出不同的模式，还具有噪声抵抗性。

4、本发明利用如Hilbert变换这种信号的时频分析方法来处理盲源分离的分解结果，可以跟踪不同模式对应的瞬时幅值和频率的变化情况。

附图说明

图1为一种基于盲源分离算法的低频振荡模式分析方法的流程示意图；

图2为经典两区四机电力系统仿真模型图；

图3为图2所示的发电机G2和发电机G1在15s-25s间的相对角速度变化曲线ω₂₁(t)；

图4为图3所示的相对角速度变化曲线ω₂₁(t)的盲源分离结果中的第一个振荡模式y₁(t)图；

图5为与图4所对应的瞬时幅值A₁(t)和瞬时频率f₁(t)的变化曲线；

图6为图3所示的相对角速度变化曲线ω₂₁(t)的盲源分离结果中的第二个振荡模式y₂(t)图；

图7为与图6所对应的瞬时幅值A₂(t)和瞬时频率f₂(t)的变化曲线。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例1：

一种基于盲源分离算法的低频振荡模式分析方法，所述方法流程图如图1所示，包括以下步骤：

1)通过装设在电力系统中的测量元件采集电力系统的低频振荡信号作为单通道时间序列；

2)根据Takens嵌入定理利用步骤1)所述的单通道时间序列构造多通道观测信号矩阵；

3)利用盲源分离算法对步骤2)所述的多通道观测信号矩阵进行分解，得到近似的源信号矩阵，所述源信号矩阵包括单通道时间序列中所包含的低频振荡信号的不同模式；

4)结合时频分析理论，求出步骤3)所述的低频振荡信号的不同模式的瞬时幅值和瞬时频率；

5)对步骤4)所述的瞬时频率进行平均求解得到平均频率，对步骤4)所述的瞬时幅值的对数相对时间的变化曲线进行最小二次拟合求得平均衰减系数。

优选的，步骤1)中所述电力系统的低频振荡信号一般为任意两台发电机的相对角速度、传输线上的有功功率、节点电压或节点间相对角偏差中的一种。

优选地，所述任意两台发电机的相对角速度，通过下式进行计算：

ω_mn＝ω_m-ω_n

其中，ω_m和ω_n分别表示第m台发电机与第n台发电机的角速度，ω_mn为第m台发电机与第n台发电机的相对角速度，ω_mn即作为单通道时间序列。

优选地，步骤2)中所述的多通道观测信号矩阵的构造过程如下：

$X\left(t\right)=\left[\begin{array}{c}x\left(t\right)\\ x(t+d)\\ \·\\ \·\\ \·\\ x\left(t\right(k-1)d)\\ \·\\ \·\\ \·\\ x(t+(M-1)d)\end{array}\right]=\left[\begin{array}{c}{X}_1\\ X_2\\ \·\\ \·\\ \·\\ X_k\\ \·\\ \·\\ \·\\ X_M\end{array}\right]$

其中X(t)表示多通道观测信号矩阵，x(t)表示单通道时间序列，d表示时间延迟，X_k表示第k个通道观测信号，k＝1,2，…，M，M为X(t)的傅里叶变换频谱中的主成分个数的两倍。

该多通道观测信号矩阵的特点是：仅根据单通道时间序列构造多通道观测信号矩阵，只要合理选择时间延迟d和矩阵维数M，就可以包含复杂动态系统的全部信息，从而有利于分离出所有不同的低频振荡模式。

优选地，步骤3)中所述的利用盲源分离算法对多通道观测信号矩阵进行分解采用的算法是二阶盲辨识(Second-order blind identification,SOBI)，其处理过程如下：

一、对多通道观测信号矩阵X(t)进行预白化处理，得到预白化矩阵Z(t)，其中Z(t)＝W_mX(t)，W_m为白化矩阵；

二、通过联合近似对角化技术(Joint approximate diagonalization,JAD)将预白化矩阵Z(t)的一组时延协方差矩阵R_z(τ)(τ∈{τ_j|j＝1,2,…,p})同时最大程度对角化，求出最大程度对角化的矩阵Ψ。其中，保证对角化程度最大的方法是令非对角线元素的和最小，即：

$min\underset{j=1}{\overset{p}{\Σ}}off\left({\mathrm{\Ψ}}^T{R}_z\right({\mathrm{\τ}}_j\left)\mathrm{\Ψ}\right)$

其中off表示矩阵Ψ^TR_z(τ_j)Ψ的非对角线元素，τ_j表示时延；

三、根据矩阵Ψ和白化矩阵W_m求解混矩阵W和分离信号矩阵Y(t)：

W＝Ψ^TW_m

Y(t)＝Wx(t)

分离信号矩阵Y(t)包含单通道时间序列中所包含的低频振荡信号的不同模式；

优选地，步骤4)中所述的低频振荡信号的不同模式y_i(t)的瞬时幅值A_i(t)和瞬时频率f_i(t)的计算方法包括Hilbert变换法、傅里叶变换法、小波变换法和基于Teager能量算子的方法。

优选地，所述利用Hilbert变换法计算低频振荡信号的不同模式y_i(t)的瞬时幅值A_i(t)和瞬时频率f_i(t)的公式如下：

$H\left(y_i\right(t\left)\right)=\frac{1}{\mathrm{\π}}\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}\frac{y_i\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}d\mathrm{\τ}$

$A_i\left(t\right)=\sqrt{y_i{\left(t\right)}^2+H{\left(y_i\right(t\left)\right)}^2}$

$f_i\left(t\right)=\frac{d\left(arctan\right(H\left(y_i\right(t))/y_i\left(t\right)))}{2\mathrm{\π}dt}$

其中H()表示Hilbert变换，i＝1,2,…N，N为源信号个数。

优选地，步骤5)中所述的求解平均频率f_iav和平均衰减系数σ_i的方法如下：

$f_{iav}=\underset{k=1}{\overset{K}{\Σ}}{f}_i\[k\]$

${\mathrm{\σ}}_i=\frac{d\left(\ln \right(A_i\left(t\right)\left)\right)}{dt}$

其中，f_iav表示第i个低频振荡模式的频率平均值，K表示瞬时频率f_i在对应的数据窗口内的采样点总数，平均衰减系数σ_i可以利用最小二乘法线性拟合ln(A_i(t))-t曲线求得。

本实施例的电力系统仿真模型如图2所示，该电力系统的所有参数都按照经典四机电力系统设置，为诱发该系统产生低频振荡信号，在时间t＝15s时，令图2所示的节点7处的负荷增加1％，该系统的频率和采样频率分别为60Hz和100Hz，数据窗的长度为10s，即1000个采样点，在步骤2)中构造多通道观测信号矩阵时，时间延迟d取10，矩阵维数M取4，取图2所示的发电机2和发电机1之间的相对角速度ω₂₁作为低频振荡信号进行研究，ω₂₁在15s-25s内随时间的变化曲线如图3所示，从图3中可以看出，在15s之后，即在节点7处的负荷增加之后，该电力系统仿真模型发生了低频振荡。

处理信号发现：矩阵Y(t)中包含两个源信号y₁(t)和y₂(t)，分别如图4和图6所示。从图中可以看出，分解得到的两个源信号是指数型衰减的正弦函数，分别代表两种低频振荡模式。其对应的瞬时幅值和频率随时间的变化情况分别如图5和图7所示。由图5可以求出y₁(t)的平均频率为0.6071Hz，平均衰减系数为-0.0450；由图7可以求出y₂(t)的平均频率为1.0028Hz，平均衰减系数为-0.4859。这两种模式的参数计算结果，与Prony算法分解得到的两种模式的参数十分接近(0.6100Hz，-0.0480；1.000Hz，-0.4600)。但比起Prony算法，本发明的方法具有动态跟踪参数变化的优点。虽然HHT算法也具有动态跟踪参数变化的特点，但是其分解准确度受到经验模态分解的影响，当低频振荡模式频率相互接近或采样频率不够高的情况下，经验模态分解的分解结果可能出现虚假模式或模式混叠情况。在本实施例中，用HHT处理ω₂₁(t)得到的分解结果就偏离Prony算法的结果(0.5469Hz，-0.1488；0.9571Hz，-0.1002)。所以相比于Prony算法和HHT算法，本发明方法分别在动态参数追踪性和模式分解准确性上更胜一筹。

以上所述，仅为本发明专利较佳的实施例，但本发明专利的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明专利所公开的范围内，根据本发明专利的技术方案及其发明专利构思加以等同替换或改变，都属于本发明专利的保护范围。

Claims (8)

1.一种基于盲源分离算法的低频振荡模式分析方法，其特征在于：所述方法包括以下步骤：

1)通过装设在电力系统中的测量元件采集电力系统的低频振荡信号作为单通道时间序列；

2)根据Takens嵌入定理利用步骤1)所述的单通道时间序列构造多通道观测信号矩阵；

3)利用盲源分离算法对步骤2)所述的多通道观测信号矩阵进行分解，得到分离信号矩阵，所述分离信号矩阵包括单通道时间序列中所包含的低频振荡信号的不同模式；

4)结合时频分析理论，求出步骤3)所述的低频振荡信号的不同模式的瞬时幅值和瞬时频率；

5)对步骤4)所述的瞬时频率进行平均求解得到平均频率，对步骤4)所述的瞬时幅值的对数相对时间的变化曲线进行最小二次拟合求得平均衰减系数。

2.根据权利要求1所述的一种基于盲源分离算法的低频振荡模式分析方法，其特征在于：步骤1)中所述电力系统的低频振荡信号为任意两台发电机的相对角速度、传输线上的有功功率、节点电压或节点间相对角偏差中的一种。

3.根据权利要求2所述的一种基于盲源分离算法的低频振荡模式分析方法，其特征在于：所述任意两台发电机的相对角速度，通过下式进行计算：

ω_mn＝ω_m-ω_n

其中，ω_m和ω_n分别表示第m台发电机与第n台发电机的角速度，ω_mn为第m台发电机与第n台发电机的相对角速度，ω_mn即作为单通道时间序列。

4.根据权利要求1所述的一种基于盲源分离算法的低频振荡模式分析方法，其特征在于：步骤2)中所述的多通道观测信号矩阵的构造过程如下：

$X\left(t\right)=\left[\begin{array}{c}x\left(t\right)\\ x(t+d)\\ .\\ .\\ .\\ x(t+(k-1\left)d\right)\\ .\\ .\\ .\\ x(t+(M-1\left)d\right)\end{array}\right]=\left[\begin{array}{c}{X}_1\\ X_2\\ .\\ .\\ .\\ X_k\\ .\\ .\\ .\\ X_M\end{array}\right]$

其中X(t)表示多通道观测信号矩阵，x(t)表示单通道时间序列，d表示时间延迟，X_k表示第k个通道观测信号，k＝1,2，…，M，M为x(t)的傅里叶变换频谱中的主成分个数的两倍。

5.根据权利要求1所述的一种基于盲源分离算法的低频振荡模式分析方法，其特征在于：步骤3)中所述的利用盲源分离算法对多通道观测信号矩阵进行分解采用的算法是二阶盲辨识，其处理过程如下：

一、对多通道观测信号矩阵X(t)进行预白化处理，得到预白化矩阵Z(t)，其中Z(t)＝W_mX(t)，W_m为白化矩阵；

二、通过联合近似对角化技术将预白化矩阵Z(t)的一组时延协方差矩阵R_z(τ)(τ∈{τ_j|j＝1,2,…,p})同时最大程度对角化，求出最大程度对角化的矩阵Ψ，其中，保证对角化程度最大的方法是令非对角线元素的和最小，即：

$min\underset{i=1}{\overset{p}{\Σ}}off\left({\mathrm{\Ψ}}^T{R}_z\right({\mathrm{\τ}}_j\left)\mathrm{\Ψ}\right)$

其中off表示矩阵Ψ^TR_z(τ_j)Ψ的非对角线元素，τ_j表示时延。

6.根据权利要求1所述的一种基于盲源分离算法的低频振荡模式分析方法，其特征在于：步骤4)中所述的低频振荡信号的不同模式y_i(t)的瞬时幅值A_i(t)和瞬时频率f_i(t)的计算方法包括Hilbert变换法、傅里叶变换法、小波变换法和基于Teager能量算子的方法。

7.根据权利要求6所述的一种基于盲源分离算法的低频振荡模式分析方法，其特征在于：所述利用Hilbert变换法计算低频振荡信号的不同模式y_i(t)的瞬时幅值A_i(t)和瞬时频率f_i(t)的公式如下：

$H\left(y_i\right(t\left)\right)=\frac{1}{\mathrm{\π}}\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}\frac{y_i\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}d\mathrm{\τ}$

$A_i\left(t\right)=\sqrt{y_i{\left(t\right)}^2+H{\left(y_i\right(t\left)\right)}^2}$

$f_i\left(t\right)=\frac{d\left(\arctan \right(H\left(y_i\left(t\right)\right)/y_i\left(t\right)\left)\right)}{2\mathrm{\π}dt}$

其中H()表示Hilbert变换，i＝1,2,…N，N为源信号个数。

8.根据权利要求1所述的一种基于盲源分离算法的低频振荡模式分析方法，其特征在于：步骤5)中所述的求解平均频率f_iav和平均衰减系数σ_i的方法如下：

$f_{iav}=\underset{k=1}{\overset{K}{\Σ}}{f}_i\[k\]$

${\mathrm{\σ}}_i=\frac{d\left(ln\right(A_i\left(t\right)\left)\right)}{dt}$

其中，f_iav表示第i个低频振荡模式的频率平均值，K表示瞬时频率f_i在对应的数据窗口内的采样点总数，平均衰减系数σ_i可以利用最小二乘法线性拟合ln(A_i(t))-t曲线求得。