一种振荡信号模态参数的识别方法

本发明涉及振荡模态参数识别领域，尤其是涉及一种振荡信号模态参数的识别方法。
背景技术：
：随着电力技术的不断进步，大电网互联已逐步实现，加之快速高放大倍数励磁装置的广泛使用，由此带来的低频振荡问题时有发生。低频振荡通常出现在远距离、重负荷输电线路上，或者互联系统的弱联络线上，在采用快速响应高放大倍数励磁系统的条件下更容易出现。快速励磁调节器的时间常数大为减少，这有效地改善了电压调节特性，但其产生的附加阻尼为负值，抵消了系统本身所固有的正阻尼，使系统的总阻尼减少或成为负值；此外电网负荷过重和电网之间的弱互联也会使系统的阻尼降低，以至系统在扰动作用后的功率振荡长久不能平息。低频振荡的频率一般在0.2—2.5hz之间，低频振荡会引起联络线过流跳闸或系统与系统或机组与系统之间的失步而解列，严重威胁电力系统的稳定运行。及时准确地对低频振荡特征进行分析，在低频振荡可能对电网造成严重危害前发出预警信息，可以使电力部门采取相应措施、抑制电网严重低频振荡现象的发生，从而有效提高电网运行的稳定性和安全性。目前针对低频振荡问题所采用的分析方法主要有傅立叶变换、小波变换、kalman滤波法、矩阵束辨识法、prony算法、随机子空间识别(ssi)算法和希尔伯特—黄(hht)算法等。傅立叶变换无法分析出阻尼特性和局部特性，所以不适于非线性、非平稳信号。小波变换存在频率交叠和自适应基选取问题，只适合于瞬态和非平稳信号。kalman滤波法能消除噪声的影响，对不同输入信号的适应性较好，但计算精度和收敛速度受初始参数设置的影响很大。矩阵束辨识法能够准确估计系统的振荡模态，并具有较强的抗噪声能力，但若信号存在时变特性，该算法的计算误差较大，无法揭示振荡的动态特性。prony方法可以提取出振荡信号模式、相位角和阻尼等信息，但存在受噪声影响大、计算速度慢和定阶问题不确定等问题。hht算法是近年发展起来的一种新型的适于非平稳、非线性信号的分析方法，传统的hht算法受端点效应的影响，虽然能得到振荡模态的瞬时频率、瞬时幅值和衰减因子，难以达到较高的计算精度。技术实现要素：本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种振荡信号模态参数的识别方法。本发明的目的可以通过以下技术方案来实现：一种振荡信号模态参数的识别方法，包括以下步骤：s1：获取系统的原始信号，通过经验模态分解方法对原始信号进行处理，获取新信号；s2：用随机子空间识别法对s1获取的新信号进行处理，获取系统的频率和阻尼比；s3：采用prony方法对s1获取的新信号进行处理，获取系统的频率、振幅和相角；s4：基于频率相同法则对s2中获取的系统的频率和阻尼比、s3中获取的系统的频率、振幅和相角进行模态参数配对，获取模态参数。所述的经验模态分解方法即emd(empiricalmodedecomposition)方法。优选地，所述的步骤s1具体包括：s11：对系统的原始信号x(t)求导，提取x(t)的极值；s12：拟合x(t)的上包络线和下包络线，并得到包络线的极大值emax(t)和极小值emim(t)；s13：计算上包络线和下包络线的包络线均值m(t)，将原始信号x(t)减去包络线均值m(t)，得到一个差值信号d(t)；s14：判断d(t)是否满足本征模态函数条件，若是，将d(t)作为一个从x(t)分解出的一个本征模态函数imf分量，进入步骤s15，否则，根据公式：x(t)＝d(t)更新原始信号，返回步骤s11；s15：判断是否满足停止筛分条件，若是，得到原始信号的一组本征模态函数imf与一个残差r，进入步骤s16，否则，根据公式：x(t)＝x(t)-d(t)更新原始信号，返回步骤s11；s16：根据原始信号的一组本征模态函数imf与一个残差r获取新信号c(t)。优选地，所述的步骤s16中获取新信号c(t)的公式为：优选地，s14中所述的本证模态函数条件为：极值点和过零点的数目相等或至多差1，且分别连接局部极大值和局部极小值所形成的两条包络线的均值在任一点处为零。优选地，s15中满足下列条件中任意一条即满足停止筛分条件：相对误差小于误差阈值、或迭代次数大于最大迭代次数、生成的imf数量达到指定值、残差信号极值个数小于指定值。优选地，所述的s2具体包括：s21：构建系统离散状态空间模型：其中，xk为系统状态量，yk为测量输出量，a为状态矩阵，c为输出矩阵，wk为零均值过程噪声，vk为零均值测量噪声，并对新信号c(t)构建hankel矩阵，所述的hankel矩阵h为：s22：根据hankel矩阵中的参数yp、yf构建由输出协方差序列组成的toeplitz矩阵，所述的toeplitz矩阵t1/i为：其中，ri为输出协方差，s23：对toeplitz矩阵进行奇异值分解，获取状态矩阵a；s24：对状态矩阵a进行特征值分解，获取系统的特征值；s25：根据系统的特征值计算系统的频率和阻尼比：优选地，所述的步骤s23具体包括：s231：对toeplitz矩阵进行奇异值分解(svd)：其中，u、v分别为左奇异向量正交矩阵、右奇异向量正交矩阵，s为正奇异矩阵组成的对角阵；s1＝diag[σi]；σ1≥σ2≥…≥σn≥0，σi为奇异值，n为系统阶数；s232：根据公式：计算状态矩阵a。优选地，所述的步骤s24具体包括：s241：对状态矩阵进行特征值分解的公式为：a＝ψλψ-1其中，λ＝diag[μi]表示一个对角n阶矩阵，由离散时间复特征值μi组成，ψ是以特征向量为列向量组成的n阶矩阵；s242：获取系统的特征值λi的公式为：其中，δt为采样时间间隔。优选地，所述的步骤s3中的prony方法采用p个具有任意幅值、相位、频率和衰减因子的指数函数的线性组合对s1获取的新信号c(t)进行等间距采样，获取系统的频率、振幅和相角。优选地，所述的步骤s3具体包括：s31：选取离散时间模型其中，zm为，且zm＝exp[(αm+j2πfm)δt]，bm为，且bm＝amexp(jθm)，am、fm、θm、αm分别表示振荡的幅值、频率、初相和衰减因子，n为原始信号采样点数，δt为采样的时间间隔；s32：构建误差函数ε：并将转化为求解常系数线性差分方程：对系数aq进行最小二乘估计，可得一组线性矩阵方程：定义样本函数r(g,h)为：构建扩展矩阵re为：其中，pe为扩展后的阶数。s33：用奇异值分解和最小二乘法(svd-tls)确定扩展矩阵re的有效秩p，得prony算法的法方程形式：求解此法方程，即可得到系数a1,…,ap。s34：求特征多项式1+a1z-1+a2z-2+…+apz-p＝0的根z1,…,zp，并计算近似值s35：将离散时间模型转换为含未知参数b的线性方程：根据s34中获得的z1,…,zp求系数向量b1,…,bp：其中，b为系数向量b1,…,bp构成的矩阵，z为z1,…,zp构成的矩阵，h表示把原矩阵中每个元素求共轭再转置，即：zh＝(z*)t；s36：根据zm、bm计算信号的振幅、相角和频率。优选地，所述步骤s36中计算信号的振幅、相角和频率的公式为：优选地，所述的步骤s4具体包括：s41：获取s2求得的频率fs的维度ni与s3求得的频率fp的维度mi；s42：判断迭代次数ni是否小于等于ni，若是，则进入步骤s43，否则进入步骤s5；s43：提取fsi，计算fsi与fpm的绝对误差δf、相对误差δf，并获取绝对误差δf的最小值δfmin；s44：判断迭代次数ji是否小于等于最大迭代次数mi，若是，进入步骤s45，否则返回步骤s42；s45：依次判断fpm是否满足配对条件，若满足配对条件，进入步骤s46，否则返回步骤s44进行下一次迭代，所述的配对条件为相对误差δf小于误差允许值mpe且该绝对误差δf为最小值δfmin；s46：进行模态参数配对：fopd＝fsi,ξopd＝ξi,θopd＝θm,aopd＝am，其中，fo、ξo、θo、ao分别为将输出模态参数的频率、阻尼比、相位角和振幅，pd为成功配对的对数。与现有技术相比，本发明具有如下优点：(1)本发明分别利用prony算法、随机子空间算法对信号进行处理，并基于频率相同法则对两个方法获取的参数进行模态参数配对，prony算法对于信号频率、相角和振幅获取的精度较高，随机子空间识别法对信号的频率、阻尼比获取的精度较高，本发明基于频率相同法则进行模态参数配对后，将prony算法、随机子空间算法获取的精度高的参数作为最后匹配后的参数结果，取长补短，以获得振荡信号精确而完整的模态参数；(2)本发明基于经验模态分解方法对原始信号进行处理，在保留原始信号主要信息的前提下，降噪滤波以消除振荡信号中的高频噪声，克服了prony算法对噪声的敏感，避免了随机子空间算法在处理非线性、非平稳信号过程中所产生的虚假模态。附图说明图1为本发明方法流程框图；图2为经验模态分解方法流程图；图3为随机子空间算法流程图；图4为prony算法流程图；图5为参数配对流程图；图6为经验模态分解方法分解得到的imfs和残差r；图7为经验模态分解方法整形后的c(t)与原始信号x(t)；图8为各成分波形图；图9为重构后的模态信号与信号c(t)。具体实施方式下面结合附图和具体实施例对本发明进行详细说明。注意，以下的实施方式的说明只是实质上的例示，本发明并不意在对其适用物或其用途进行限定，且本发明并不限定于以下的实施方式。实施例一种振荡信号模态参数的识别方法，如图1所示，包括以下步骤：s1：获取系统的原始信号，通过经验模态分解方法对原始信号进行处理，获取新信号。本实施例中，获取的原始信号为一组低频振荡采样信号，采样频率100hz，采样时间0-20s，共计2000个采样点。如图2所示，步骤s1具体包括：s11：对系统的原始信号x(t)求导，提取x(t)的极值；s12：拟合x(t)的上包络线和下包络线，并得到包络线的极大值emax(t)和极小值emim(t)；本实施例中，采用三次样条插值函数拟合上包络线和下包络线。s13：计算上包络线和下包络线的包络线均值m(t)，将原始信号x(t)减去包络线均值m(t)，得到一个差值信号d(t)；s14：判断d(t)是否满足本征模态函数条件，若是，将d(t)作为一个从x(t)分解出的一个本征模态函数imf分量，进入步骤s15，否则，根据公式：x(t)＝d(t)更新原始信号，返回步骤s11；s14中所述的本证模态函数条件为：极值点和过零点的数目相等或至多差1，且分别连接局部极大值和局部极小值所形成的两条包络线的均值在任一点处为零；s15：判断是否满足停止筛分条件，若是，得到原始信号的一组本征模态函数imf与一个残差r，进入步骤s16，否则，根据公式：x(t)＝x(t)-d(t)更新原始信号，返回步骤s11；s15中满足下列条件中任意一条即满足停止筛分条件：相对误差小于误差阈值、或迭代次数大于最大迭代次数、生成的imf数量达到指定值、残差信号极值个数小于指定值。s16：根据原始信号的一组本征模态函数imf与一个残差r获取新信号c(t)。具体地，步骤s16中获取新信号c(t)的公式为：通过经验模态分解(emd)方法得到的一组imf与一个直流分量r，如图6所示；整形得到的新信号c(t)与原始信号x(t)波形几乎重合，如图7所示。由此可见，emd方法降噪滤波后仍保留了原始信息的重要成分。s2：用随机子空间识别法对s1获取的新信号进行处理，获取系统的频率和阻尼比；s3：采用prony方法对s1获取的新信号进行处理，获取系统的频率、振幅和相角。具体地，本实施例中，步骤s3中的prony方法采用p个具有任意幅值、相位、频率和衰减因子的指数函数的线性组合对s1获取的新信号c(t)进行等间距采样，获取系统的频率、振幅和相角。s4：基于频率相同法则对s2中获取的系统的频率和阻尼比、s3中获取的系统的频率、振幅和相角进行模态参数配对，获取模态参数。如图5所示，步骤s4具体包括：s41：获取s2求得的频率fs的维度ni与s3求得的频率fp的维度mi；s42：判断迭代次数ni是否小于等于ni，若是，则进入步骤s43，否则进入步骤s5；s43：提取fsi，计算fsi与fpm的绝对误差δf、相对误差δf，并获取绝对误差δf的最小值δfmin；s44：判断迭代次数ji是否小于等于最大迭代次数mi，若是，进入步骤s45，否则返回步骤s42；s45：依次判断fpm是否满足配对条件，若满足配对条件，进入步骤s46，否则返回步骤s44进行下一次迭代，所述的配对条件为相对误差δf小于误差允许值mpe且该绝对误差δf为最小值δfmin；s46：进行模态参数配对：fopd＝fsi,ξopd＝ξi,θopd＝θm,aopd＝am，其中，fo、ξo、θo、ao分别为将输出模态参数的频率、阻尼比、相位角和振幅，pd为成功配对的对数。本实施例中，完成模态匹配后，还包括步骤s5：s5：输出各信号成分的频率、振幅、阻尼比和相角。此外，本实施例中还输出各信号成分波形图；将提取出的各信号成分重构，并与c(t)对比。本实施例中，s2、s3的具体计算步骤分别为：如图3所示，s2具体包括：s21：构建系统离散状态空间模型：其中，xk为系统状态量，yk为测量输出量，a为状态矩阵，含信号的状态参数，c为输出矩阵，含所采集到的信号体现的信息，wk为零均值过程噪声，vk为零均值测量噪声，wk、vk为处理过程和建模误差引起的互不相关的噪声，测量输出量yk为c(t)。并对新信号c(t)构建hankel矩阵，所述的hankel矩阵h为：s22：根据hankel矩阵中的参数yp、yf构建由输出协方差序列组成的toeplitz矩阵，所述的toeplitz矩阵t1/i为：其中，ri为输出协方差，s23：对toeplitz矩阵进行奇异值分解，获取状态矩阵a；s231：对toeplitz矩阵进行奇异值分解(svd)：其中，u、v分别为左奇异向量正交矩阵、右奇异向量正交矩阵，s为正奇异矩阵组成的对角阵；s1＝diag[σi]；σ1≥σ2≥…≥σn≥0，σi为奇异值，n为系统阶数；s232：根据公式：计算状态矩阵a；s24：对状态矩阵a进行特征值分解，获取系统的特征值；s241：对状态矩阵进行特征值分解的公式为：a＝ψλψ-1其中，λ＝diag[μi]表示一个对角n阶矩阵，由离散时间复特征值μi组成，ψ是以特征向量为列向量组成的n阶矩阵；s242：获取系统的特征值λi的公式为：其中，δt为采样时间间隔。s25：根据系统的特征值计算系统的频率和阻尼比：本实施例中，步骤s2获取的频率和阻尼比见表1。表1.ssi求解得到的频率和阻尼比成分频率/hz阻尼比10.53900.010420.88550.042430.97950.089141.08010.013051.95420.0767如图4所示，步骤s3具体包括：s31：选取离散时间模型其中，zm为，且zm＝exp[(αm+j2πfm)δt]，bm为，且bm＝amexp(jθm)，am、fm、θm、αm分别表示振荡的幅值、频率、初相和衰减因子，n为原始信号采样点数，δt为采样的时间间隔；s32：构建误差函数ε：并将转化为求解常系数线性差分方程：对系数aq进行最小二乘估计，可得一组线性矩阵方程：定义样本函数r(g,h)为：构建扩展矩阵re为：其中，pe为扩展后的阶数。s33：用奇异值分解和最小二乘法(svd-tls)确定扩展矩阵re的有效秩p，得prony算法的法方程形式：求解此法方程，即可得到系数a1,…,ap。s34：求特征多项式1+a1z-1+a2z-2+…+apz-p＝0的根z1,…,zp，并计算近似值s35：将离散时间模型转换为含未知参数b的线性方程：根据s34中获得的z1,…,zp求系数向量b1,…,bp：其中，b为系数向量b1,…,bp构成的矩阵，z为z1,…,zp构成的矩阵，h表示把原矩阵中每个元素求共轭再转置，即：zh＝(z*)t；s36：根据zm、bm计算信号的振幅、相角和频率。s36中计算信号的振幅、相角和频率的公式为：本实施例中，步骤s2获取的98组结果振幅、相角、频率见表2。表2.prony求解得到的频率、振幅、相角本实施例中，模态参数配对结果见表3，根据模态参数配对结果，输出各信号成分的频率fo、振幅ao、阻尼比ξo和相角θo，将提取出的各信号成分重构，并与c(t)对比，如图9所示，由表3和图8可知，融合了emd、ssi、prony的模态参数识别法从原始信号x(t)中提取出了两个主要成分，且得到了相应的频率、振幅、阻尼比和相角。并且，由图9可见，本发明计算结果拟合出的曲线与原始信号曲线高度吻合，体现出良好的精度和鲁棒性。表3.模态参数配对结果成分振幅/um频率/hz阻尼比相角/rad10.05820.53900.01042.836920.35240.97950.08911.8578上述实施方式仅为例举，不表示对本发明范围的限定。这些实施方式还能以其它各种方式来实施，且能在不脱离本发明技术思想的范围内作各种省略、置换、变更。当前第1页12