CN110048416A - S-g滤波和自适应mp算法的低频振荡模态辨识方法 - Google Patents

S-g滤波和自适应mp算法的低频振荡模态辨识方法 Download PDF

Info

Publication number
CN110048416A
CN110048416A CN201910376969.1A CN201910376969A CN110048416A CN 110048416 A CN110048416 A CN 110048416A CN 201910376969 A CN201910376969 A CN 201910376969A CN 110048416 A CN110048416 A CN 110048416A
Authority
CN
China
Prior art keywords
signal
matrix
rank
low
frequency oscillation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201910376969.1A
Other languages
English (en)
Other versions
CN110048416B (zh
Inventor
李昕
任长安
陈坚
陈利平
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hunan Institute of Technology
Original Assignee
Hunan Institute of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hunan Institute of Technology filed Critical Hunan Institute of Technology
Priority to CN201910376969.1A priority Critical patent/CN110048416B/zh
Publication of CN110048416A publication Critical patent/CN110048416A/zh
Application granted granted Critical
Publication of CN110048416B publication Critical patent/CN110048416B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J3/00—Circuit arrangements for ac mains or ac distribution networks
    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J3/00—Circuit arrangements for ac mains or ac distribution networks
    • H02J3/24—Arrangements for preventing or reducing oscillations of power in networks
    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J2203/00—Indexing scheme relating to details of circuit arrangements for AC mains or AC distribution networks
    • H02J2203/20—Simulating, e g planning, reliability check, modelling or computer assisted design [CAD]
    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J3/00—Circuit arrangements for ac mains or ac distribution networks
    • H02J3/002—Flicker reduction, e.g. compensation of flicker introduced by non-linear load

Landscapes

  • Engineering & Computer Science (AREA)
  • Power Engineering (AREA)
  • Complex Calculations (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

S‑G滤波和自适应MP算法的低频振荡模态辨识方法，包括以下步骤：S1：提取WAMS的低频振荡信号；S2：采用S‑G滤波器对低频振荡信号进行S‑G滤波预处理；S3：评估S‑G滤波器的消噪效果；S4：提取并保存经S‑G滤波预处理的低频振荡信号；S5：设置矩阵束参数L，构建矩阵H；S6：对H进行奇异值分解，用奇异值累计占比相邻增量比方法进行定阶；S7：构造矩阵束Y2‑λY1；S8：求解矩阵束Y2‑λY1的广义特征值λi；S9：求解频率fi和阻尼因子αi；S10：求解初始幅值Ai和初相位θi；S11：拟合出估计信号，判断拟合精度是否满足要求；S12：输出各模态参数值。本发明自适应性好、运算量小、定阶准确。

Description

S-G滤波和自适应MP算法的低频振荡模态辨识方法
技术领域
本发明涉及广域测量系统关于电网稳定性分析技术领域，具体是涉及一种S-G滤波和自适应MP算法的低频振荡模态辨识方法。
背景技术
国家电网规模日益庞大，同时风机、光伏、储能等新能源大规模并入电网，电网运行的规模和复杂性随之增加，相应电网阻尼弱而导致低频振荡问题的可能随之增加，威胁到电网的输电和稳定运行。快速、实时、准确辨识振荡模态并定位相应的振荡区域，进而加以抑制是电力系统稳定运行的关键技术。
随着电网互联规模的加大，较早的基于系统动态模型的特征值分析方法面临系统模型高阶的挑战和限制，同时受系统非线性问题影响而使分析产生较大误差，且计算复杂度增加，运算时间也无法满足在线辨识的需求。上世纪90年代起随着WAMS(Wide AreaMeasurement System，广域监测系统)的出现和迅速发展，该系统利用相位测量单元和全球定位系统时间戳确保测量的同步性，将数据回传电网高速处理中心完成实时数据侦测，现代低频振荡辨识方法便是基于这些来自电网各测量点提供的在线实时数据。Prony算法是较早应用于低频振荡模态辨识的方法，传统Prony算法可以获得各个模态的模式参数：频率、阻尼因子等，然而该方法对噪声敏感易致辨识结果失准，且定阶过程因需要人为阈值设定而缺少自适应性，从而导致虚假模态的产生或者真实主导模态的遗漏。希尔伯特-黄变换(Hilbert-Huang Transform，HHT)适用于分析非线性、非平稳信号，但缺乏完善的数学解释，而且在应用过程中存在模态混叠、端点效应等问题。MP(Matrix Pencil)算法可实现谐波恢复和幅值振荡余弦信号的参数估计，是电力系统辨识领域近年备受关注的代表性辨识工具。MP模态辨识的过程中，在系统设备及信号传输过程中的热噪声影响较大的情况下，将导致辨识的失准甚至错误；因此进行预处理尤其是消噪处理，减轻噪声对后续辨识的影响则显得非常必要。另一方面，在MP算法的分析过程中要进行SVD分解以及奇异值矩阵的准确定阶，那么如何减少或者避免凭主观人为经验设定阈值对准确定阶的影响，也是该算法在进行低频振荡模态辨识过程中亟待解决的关键问题。
综上可知，选择有效的消噪预处理方法进行低频振荡系统噪声的抑制，以及后续辨识过程中准确选定奇异值矩阵的阶数是电力系统低频振荡模态辨识中的关键问题，实现电力系统低频振荡实时、快速、准确的辨识，从而采取相应的抑制措施以保证电网稳定运行具有重要意义。
发明内容
本发明所要解决的技术问题是，克服上述背景技术的不足，提供一种S-G滤波和自适应MP算法的低频振荡模态辨识方法，在电力系统广域测量系统中低频振荡信号受噪声影响下，保证模态及其参数的实时、快速、准确辨识，保证电网采取合理有效的抑制措施，从而提高电网稳定运行能力。
本发明解决其技术问题采用的技术方案是，一种S-G滤波和自适应MP算法的低频振荡模态辨识方法，包括以下步骤：
步骤S1：提取WAMS的低频振荡信号；WAMS的低频振荡实测信号表示为y(n)＝x(n)+w(n)；x(n)是低频振荡理想信号，w(n)是噪声信号，n＝0,1,2,…,N-1，信号长度为N；
步骤S2：采用S-G滤波器对低频振荡信号进行S-G滤波预处理；经过S-G滤波预处理后的信号表示为r(n)是残余噪声；
步骤S3：评估S-G滤波器的消噪效果；
步骤S4：提取并保存经S-G滤波预处理的低频振荡信号n＝0,1,2,…,N-1；
步骤S5：设置矩阵束参数L，构建由预处理后的信号构造的Hankel矩阵H；
步骤S6：对H进行奇异值分解，即H＝UDVT，其中U为(N-L)×(N-L)阶正交矩阵，V为(L+1)×(L+1)阶正交矩阵，VT为V的共轭转置；D为(N-L)×(L+1)阶矩阵，其主对角元素为H的非负奇异值，用σi表示，i＝1,2,…,L+1；并用奇异值累计占比相邻增量比方法进行定阶；
步骤S7：构造矩阵束Y2-λY1；
步骤S8：求解矩阵束Y2-λY1的广义特征值λi；
步骤S9：求解低频振荡信号各模态的频率fi和阻尼因子αi；
步骤S10：求解低频振荡信号各模态的初始幅值Ai和初相位θi；
步骤S11：根据频率fi、阻尼因子αi、初始幅值Ai和初相位θi拟合出估计信号，判断估计信号的拟合精度是否满足要求；
步骤S12：输出各模态参数值：频率fi、阻尼因子αi、初始幅值Ai和初相位θi。
进一步，所述步骤S2中，还包括以下步骤：
步骤S21：S-G滤波器利用移动的窗口进行信号平滑，设置窗口大小为2S+1，多项式阶数为T，在平滑窗口沿各采样点n＝0,1,2,…,N-1移动的过程中，分别将各采样点设置为窗口的中心，由各采样点为中心前后延伸S个点组成一个局部点集，用多项式对该点集进行拟合，其中pT(i)是进行拟合的多项式，i表示窗口中的第i个采样点，i＝-S…，-1,0,1…,S，S约为平滑窗口大小的一半，T为拟合多项式阶数，且T<2S+1，ak是多项式pT(i)的系数；
步骤S22：求解窗口内多项式拟合的误差，用总体误差平方和描述为：
步骤S23：使E最小，即其中r＝0,1,2,…,T；导出令则有
步骤S24：利用采样点数据y(i)，分别取r＝0,1,2,…,T，解方程组求得系数a0-aT；将各系数回代多项式拟合出滤波后的函数值；逐点移动窗口，得出所有采样点滤波后的函数值，从而达到消噪的效果。
进一步，所述步骤S3中，采用信噪比指标评估S-G滤波器的消噪效果，所述信噪比指标表示为：要求指标SNR＞20dB。
进一步，所述步骤S4中，用P阶幅值指数变化的余弦函数模型来表示低频振荡信号该模型为：
其中n＝0,1,2,…,N-1，P是信号模态数，Ts是采样间隔，fi表示频率，αi表示阻尼因子，Ai表示初始幅值，θi表示初相位；扩展到复数域的讨论，该指数模型表示为：
其中zi是信号的极点且2P为复数域模态阶数，是相应实数域幅值指数变化的余弦函数模型中余弦函数个数的2倍。
进一步，所述步骤S5中，由预处理后的信号n＝0,1,2,…,N-1构造的Hankel矩阵H如下：
进一步，所述步骤S6中，还包括以下步骤：
步骤S61：计算奇异值对应的累计占比：矩阵H的SVD分解后的奇异值序列具有非递增特性，即σ1≥σ2≥…≥σJ≥σJ+1≥…≥σi≥…≥σL+1≈0，其中L是分解得到的奇异值总体个数，i分别代表奇异值的序号，i＝1,2，…,L+1，σi表示序号i对应的各个奇异值，J是i中某一序号，并预设为真实模态阶数，取前50个奇异值，即L+1＝50，计算从第1个到第50个奇异值对应的累计占比：其中PJ是第J个奇异值对应的累计占比。
步骤S62：构造奇异值累计占比相邻增量比：RJ＝ΔPJ+1/ΔPJ＝(PJ+1-PJ)/(PJ-PJ-1)，J≥2，其中RJ为第J个奇异值累计占比相邻增量比值，当J为真实模态阶数时，取i＝1,2,…,L，L+1＝50，由Ri组成的曲线将在i＝J处出现最低点，即min(Ri)＝RJ。
进一步，所述步骤S7中，构造矩阵束Y2-λY1的方法如下：
取D前2P列构造新矩阵D'，
其中D'为(N-L)×2P阶矩阵，其前2P行是对角线元素为σ1、σ2、…、σ2P构成的2P×2P阶对角阵，余下N-L-2P行由0元素填充，D'矩阵的构造降低了噪声的影响并简化了后续矩阵运算的复杂度。取矩阵V前2P列构成(L+1)×2P阶矩阵V'，删去V'最后一行形成L×2P阶矩阵V1，删去V'第一行构成L×2P阶矩阵V2，从而构造两个矩阵Y1、Y2：Y1＝UD'V1 T，Y2＝UD'V2 T；进而得出矩阵束Y2-λY1，具体表示为：Y2-λY1＝Z1B(Z0-λI)Z2，其中Z0是对角线元素为z1,z2,…,z2P构成的对角阵，B是对角线元素为B1,B2,…,B2P构成的对角阵， zi是信号的极点，Bi表示各个模态的初始幅值，i＝1,2,…,2P。
进一步，所述步骤S10中，Bi的求取方法如下:
根据r(n)作为残余噪声影响微弱忽略不计，因此可用来估计利用最小二乘法计算其中B＝[B1,B2,…,B2P]T，进而求解各个Bi。
进一步，所述步骤S11中，判断估计信号的拟合精度是否满足要求，包括以下步骤：
步骤S111：将估计信号与预处理信号作对比，计算拟合精度AFI；
步骤S112：设定拟合精度阈值，判断拟合精度AFI是否满足拟合精度阈值要求，若满足拟合精度阈值要求，转入步骤S12，若不满足拟合精度阈值要求，返回步骤S5。
进一步，所述步骤S111中，拟合精度其中，x和分别表示低频振荡预处理信号向量和经模态辨识后估计信号向量，||·||表示2范数，AFI要求大于10dB。
与现有技术相比，本发明的优点如下：
(1)本发明通过S-G滤波器技术，对电力系统低频振荡信号进行消噪预处理，有效抑制噪声对后续模态辨识的影响。
(2)本发明利用所提出的定阶方法实现自适应定阶并能快速计算，避免凭人为经验设定阈值带来的主观因素影响，从而实时、准确辨识低频振荡主导模态。
附图说明
图1是本发明实施例的流程图。
图2是本发明实施例中的低频振荡理想信号x(n)。
图3是本发明实施例中的低频振荡实测信号y(n)。
图4是本发明实施例中的预处理低频振荡信号
图5是本发明实施例中的奇异值累计占比相邻增量比定阶示意图。
图6是本发明的数值信号和辨识拟合信号。
具体实施方式
下面结合附图及具体实施例对本发明作进一步详细描述。
参照图1，本实施例包括以下步骤：
步骤S1：提取WAMS的低频振荡信号，设置信号y(n)＝x(n)+w(n)，其中y(n)是WAMS的低频振荡实测信号，x(n)是低频振荡理想信号，w(n)是噪声信号，n＝0,1,2,…,N-1，信号长度为N。本实施例中采用测试信号x(n)＝x1(n)+x2(n)+x3(n)，n＝0,1,2,…,N-1，其中x1(n)＝2·e-0.2·n·cos(2π·1.5·n)，x2(n)＝1.5·e-0.1·n·cos(2π·0.5·n+π/6)，x3(n)＝0.5·e-0.05·n·cos(2·π·0.2·n+π/3)进行分析。x(n)的信号图如图2所示；y(n)的信号图如图3所示。
步骤S2：采用S-G滤波器对低频振荡信号进行S-G滤波预处理；
低频振荡实测信号y(n)经过S-G滤波预处理后，得到其中是y(n)经S-G滤波预处理后的信号，r(n)是残余噪声。的信号图如图4所示。
在本发明实施例中，所述步骤S2中，还包括以下步骤：
步骤S21：S-G滤波器利用移动的窗口进行信号平滑，设置窗口大小为2S+1，多项式阶数为T，在平滑窗口沿各采样点n＝0,1,2,…,N-1移动的过程中，分别将各采样点设置为窗口的中心，由各采样点为中心前后延伸S个点组成一个局部点集，用多项式对该局部点集进行拟合，其中pT(i)是进行拟合的多项式，i表示窗口中的第i个采样点，i＝-S…，-1,0,1…,S，S约为平滑窗口大小的一半，T为拟合多项式阶数，且T<2S+1，ak是多项式pT(i)的系数。
步骤S22：求解窗口内多项式拟合的误差，用总体误差平方和描述为：
步骤S23：为使E最小，要求其中r＝0,1,2,…,T；可导出令 则有
步骤S24：利用采样点数据y(i)，分别取r＝0,1,2,…,T，解方程组可求得系数a0-aT；将各系数回代多项式拟合出滤波后的函数值；逐点移动窗口，得出所有采样点滤波后的函数值，从而达到消噪的效果。
步骤S3：采用信噪比(Signal Noise Ratio，SNR)指标评估S-G滤波器的消噪效果，判断经S-G滤波预处理后的信号是否满足指标要求，若满足指标要求，则转入步骤S4；若不满足指标要求，返回步骤S2；
在本发明实施例中，在所述步骤S3中，信噪比指标：
步骤S4：提取并保存经S-G滤波预处理的低频振荡信号n＝0,1,2,…,N-1；
所述步骤S4中，用P阶幅值指数变化的余弦函数模型来表示低频振荡信号该模型为：
其中n＝0,1,2,…,N-1，P是信号模态数，Ts是采样间隔，fi、αi、Ai和θi分别表示待辨识各模态参数：频率、阻尼因子、初始幅值和初相位；扩展到复数域的讨论，该指数模型也可以表示为：
其中zi是信号的极点且2P为复数域模态阶数，是相应实数域幅值指数变化的余弦函数模型中余弦函数个数的2倍。
步骤S5：设置矩阵束参数L，L通常取N/4～N/3；构建由预处理后的信号构造的Hankel矩阵H，由预处理后的信号n＝0,1,2,…,N-1构造的Hankel矩阵H如下：
步骤S6：对H进行奇异值分解(Signal Value Decomposition，SVD)，即H＝UDVT，其中U为(N-L)×(N-L)阶正交矩阵，V为(L+1)×(L+1)阶正交矩阵，VT为V的共轭转置；D为(N-L)×(L+1)阶矩阵，其主对角元素为H的非负奇异值，用σi表示，i＝1,2,…,L+1；并用奇异值累计占比相邻增量比方法进行定阶，定阶结果如图5所示。
本实施例步骤S6中，包括以下步骤：
步骤S61：计算奇异值对应的累计占比：矩阵H的SVD分解后的奇异值序列具有非递增特性，即σ1≥σ2≥…≥σJ≥σJ+1≥…≥σi≥…≥σL+1≈0，其中L是分解得到的奇异值总体个数，i分别代表奇异值的序号，i＝1,2，…,L+1，σi表示序号i对应的各个奇异值，J是i中某一序号，并预设为真实模态阶数，为简化计算，取前50个奇异值，即L+1＝50，之后的奇异值相较极小，故忽略或舍弃其余奇异值，计算从第1个到第50个奇异值对应的累计占比，例如：其中PJ是第J个奇异值对应的累计占比；
步骤S62：构造奇异值累计占比相邻增量比RJ＝ΔPJ+1/ΔPJ＝(PJ+1-PJ)/(PJ-PJ-1)，J≥2，其中RJ为第J个奇异值累计占比相邻增量比值，当J为真实模态阶数时，取i＝1,2,…,L，L+1＝50，由Ri组成的曲线将在i＝J处出现最低点，即min(Ri)＝RJ，从而达到确定真实模态阶数J的目的，即自适应定阶。
步骤S7：构造矩阵束Y2-λY1：取D前2P列构造新矩阵D'，
其中D'为(N-L)×2P阶矩阵，其前2P行是对角线元素为σ1、σ2、…、σ2P构成的2P×2P阶对角阵，余下N-L-2P行由0元素填充，D'矩阵的构造降低了噪声的影响并简化了后续矩阵运算的复杂度。取矩阵V前2P列构成(L+1)×2P阶矩阵V'，删去V'最后一行形成L×2P阶矩阵V1，删去V'第一行构成L×2P阶矩阵V2，从而构造两个矩阵Y1、Y2：Y1＝UD'V1 T，Y2＝UD'V2 T；进而得出矩阵束Y2-λY1，具体表示为：Y2-λY1＝Z1B(Z0-λI)Z2，其中Z0是对角线元素为z1,z2,…,z2P构成的对角阵，B是对角线元素为B1,B2,…,B2P构成的对角阵， zi是信号的极点，Bi表示各个模态的初始幅值，i＝1,2,…,2P。
步骤S8：求解矩阵束Y2-λY1的广义特征值λi，λi(i＝1，2，3，…，2P)，也就是求解Y1 +Y2的特征值，其中Y1 +为Y1的广义逆矩阵，从而解得特征值λi，也就是信号的极点zi，i＝1,2,…,2P。
步骤S9：求解低频振荡信号各模态的频率fi和阻尼因子αi：αi＝Re[(lnλi)/Ts]，
步骤S10：求解低频振荡信号各模态的初始幅值Ai和初相位θi：Ai＝2|Bi|，θi＝arctan[Im(Bi)/Re(Bi)]；
所述步骤S10中，Bi的求取方法如下:
根据r(n)作为残余噪声影响微弱忽略不计，因此可用来估计利用最小二乘法计算其中B＝[B1,B2,…,B2P]T，进而求解各个Bi。
步骤S11：根据上述求解的模态参数：频率fi、阻尼因子αi、初始幅值Ai和初相位θi拟合出估计信号，判断估计信号的拟合精度是否满足要求？若满足要求，转入步骤S12，若不满足要求，返回步骤S5；
所述步骤S11中，判断估计信号的拟合精度是否满足要求，包括以下步骤：
步骤S111：将估计信号与预处理信号作对比，计算拟合精度AFI；
步骤S112：设定拟合精度阈值，判断拟合精度AFI是否满足拟合精度阈值要求，若满足拟合精度阈值要求，转入步骤S12，若不满足拟合精度阈值要求，返回步骤S5。
所述步骤S111中，拟合精度其中，x和分别表示低频振荡预处理信号向量和经模态辨识后估计信号向量，||·||表示2范数，AFI要求大于10dB。
步骤S12：输出各模态参数值：频率fi、阻尼因子αi、初始幅值Ai和初相位θi。
较佳地，辨识后的各模态参数如表1所示，低频振荡预处理信号和辨识拟合信号如图6所示。
表1辨识后各模态参数表
本发明首先利用S-G滤波技术进行消噪预处理，然后采用自适应MP算法进行模态辨识。本发明可克服电网广域测量系统及传输过程中的热噪声(白噪声为主)对低频振荡信号的污染，确保后续低频振荡模态辨识的准确性。针对MP算法在模态辨识过程中阶数确定的问题，采取奇异值累计占比相邻增量比方法进行定阶，该方法自适应性好、运算量小、定阶准确。本发明提出的两阶段低频振荡模态辨识方法，具有较好的可行性和应用价值。
本领域的技术人员可以对本发明进行各种修改和变型，倘若这些修改和变型在本发明权利要求及其等同技术的范围之内，则这些修改和变型也在本发明的保护范围之内。
说明书中未详细描述的内容为本领域技术人员公知的现有技术。

Claims (10)

1.一种S-G滤波和自适应MP算法的低频振荡模态辨识方法，其特征在于，包括以下步骤：
步骤S1：提取WAMS的低频振荡信号；WAMS的低频振荡实测信号表示为y(n)＝x(n)+w(n)；x(n)是低频振荡理想信号，w(n)是噪声信号，n＝0,1,2,…,N-1，信号长度为N；
步骤S2：采用S-G滤波器对低频振荡信号进行S-G滤波预处理；经过S-G滤波预处理后的信号表示为r(n)是残余噪声；
步骤S3：评估S-G滤波器的消噪效果；
步骤S4：提取并保存经S-G滤波预处理的低频振荡信号n＝0,1,2,…,N-1；
步骤S5：设置矩阵束参数L，构建由预处理后的信号构造的Hankel矩阵H；
步骤S6：对H进行奇异值分解，即H＝UDVT，其中U为(N-L)×(N-L)阶正交矩阵，V为(L+1)×(L+1)阶正交矩阵，VT为V的共轭转置，D为(N-L)×(L+1)阶矩阵，其主对角元素为H的非负奇异值，用σi表示，i＝1,2,…,L+1；并用奇异值累计占比相邻增量比方法进行定阶；
步骤S7：构造矩阵束Y2-λY1；
步骤S8：求解矩阵束Y2-λY1的广义特征值λi；
步骤S9：求解低频振荡信号各模态的频率fi和阻尼因子αi；
步骤S10：求解低频振荡信号各模态的初始幅值Ai和初相位θi；
步骤S11：根据频率fi、阻尼因子αi、初始幅值Ai和初相位θi拟合出估计信号，判断估计信号的拟合精度是否满足要求；
步骤S12：输出各模态参数值：频率fi、阻尼因子αi、初始幅值Ai和初相位θi。
2.如权利要求1所述的S-G滤波和自适应MP算法的低频振荡模态辨识方法，其特征在于：所述步骤S2中，还包括以下步骤：
步骤S21：S-G滤波器利用移动的窗口进行信号平滑，设置窗口大小为2S+1，多项式阶数为T，在平滑窗口沿各采样点n＝0,1,2,…,N-1移动的过程中，分别将各采样点设置为窗口的中心，由各采样点为中心前后延伸S个点组成一个局部点集，用多项式对该点集进行拟合，其中pT(i)是进行拟合的多项式，i表示窗口中的第i个采样点，i＝-S…，-1,0,1…,S，S约为平滑窗口大小的一半，T为拟合多项式阶数，且T<2S+1，ak是多项式pT(i)的系数；
步骤S22：求解窗口内多项式拟合的误差，用总体误差平方和描述为：
步骤S23：使E最小，即其中r＝0,1,2,…,T；导出令则有
步骤S24：利用采样点数据y(i)，分别取r＝0,1,2,…,T，解方程组求得系数a0-aT；将各系数回代多项式拟合出滤波后的函数值；逐点移动窗口，得出所有采样点滤波后的函数值，从而达到消噪的效果。
3.如权利要求1所述的S-G滤波和自适应MP算法的低频振荡模态辨识方法，其特征在于：所述步骤S3中，采用信噪比指标评估S-G滤波器的消噪效果，所述信噪比指标表示为：要求指标SNR＞20dB。
4.如权利要求1所述的S-G滤波和自适应MP算法的低频振荡模态辨识方法，其特征在于：所述步骤S4中，用P阶幅值指数变化的余弦函数模型来表示低频振荡信号该模型为：
其中n＝0,1,2,…,N-1，P是信号模态数，Ts是采样间隔，fi表示频率，αi表示阻尼因子，Ai表示初始幅值，θi表示初相位；扩展到复数域的讨论，该指数模型表示为：
其中zi是信号的极点且2P为复数域模态阶数，是相应实数域幅值指数变化的余弦函数模型中余弦函数个数的2倍。
5.如权利要求1所述的S-G滤波和自适应MP算法的低频振荡模态辨识方法，其特征在于：所述步骤S5中，由预处理后的信号构造Hankel矩阵H如下：
6.如权利要求1所述的S-G滤波和自适应MP算法的低频振荡模态辨识方法，其特征在于：所述步骤S6中，还包括以下步骤：
步骤S61：计算奇异值对应的累计占比：矩阵H的SVD分解后的奇异值序列具有非递增特性，即σ1≥σ2≥…≥σJ≥σJ+1≥…≥σi≥…≥σL+1≈0，其中L是分解得到的奇异值总体个数，i分别代表奇异值的序号，i＝1,2，…,L+1，σi表示序号i对应的各个奇异值，J是i中某一序号，并预设为真实模态阶数，取前50个奇异值，即L+1＝50，计算从第1个到第50个奇异值对应的累计占比：其中PJ是第J个奇异值对应的累计占比；
步骤S62：构造奇异值累计占比相邻增量比：RJ＝ΔPJ+1/ΔPJ＝(PJ+1-PJ)/(PJ-PJ-1)，J≥2，其中RJ为第J个奇异值累计占比相邻增量比值，当J为真实模态阶数时，取i＝1,2,…,L，L+1＝50，由Ri组成的曲线将在i＝J处出现最低点，即min(Ri)＝RJ。
7.如权利要求1所述的S-G滤波和自适应MP算法的低频振荡模态辨识方法，其特征在于：所述步骤S7中，构造矩阵束Y2-λY1的方法如下：
取D前2P列构造新矩阵D'，
其中D'为(N-L)×2P阶矩阵，其前2P行是对角线元素为σ1、σ2、…、σ2P构成的2P×2P阶对角阵，余下N-L-2P行由0元素填充，D'矩阵的构造降低了噪声的影响并简化了后续矩阵运算的复杂度。取矩阵V前2P列构成(L+1)×2P阶矩阵V'，删去V'最后一行形成L×2P阶矩阵V1，删去V'第一行构成L×2P阶矩阵V2，从而构造两个矩阵Y1、Y2：Y1＝UD'V1 T，Y2＝UD'V2 T；进而得出矩阵束Y2-λY1，具体表示为：Y2-λY1＝Z1B(Z0-λI)Z2，其中Z0是对角线元素为z1,z2,…,z2P构成的对角阵，B是对角线元素为B1,B2,…,B2P构成的对角阵， zi是信号的极点，Bi是各个模态的初始幅值，i＝1,2,…,2P。
8.如权利要求1所述的S-G滤波和自适应MP算法的低频振荡模态辨识方法，其特征在于：所述步骤S10中，Bi的求取方法如下:
根据r(n)作为残余噪声影响微弱忽略不计，因此可用来估计利用最小二乘法计算其中B＝[B1,B2,…,B2P]T，进而求解各个Bi。
9.如权利要求1所述的S-G滤波和自适应MP算法的低频振荡模态辨识方法，其特征在于：所述步骤S11中，判断估计信号的拟合精度是否满足要求，包括以下步骤：
步骤S111：将估计信号与预处理信号作对比，计算拟合精度AFI；
步骤S112：设定拟合精度阈值，判断拟合精度AFI是否满足拟合精度阈值要求，若满足拟合精度阈值要求，转入步骤S12，若不满足拟合精度阈值要求，返回步骤S5。
10.如权利要求1所述的S-G滤波和自适应MP算法的低频振荡模态辨识方法，其特征在于：所述步骤S111中，拟合精度其中，x和分别表示低频振荡预处理信号向量和经模态辨识后估计信号向量，||.||表示2范数，AFI要求大于10dB。
CN201910376969.1A 2019-05-07 2019-05-07 S-g滤波和自适应mp算法的低频振荡模态辨识方法 Active CN110048416B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910376969.1A CN110048416B (zh) 2019-05-07 2019-05-07 S-g滤波和自适应mp算法的低频振荡模态辨识方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910376969.1A CN110048416B (zh) 2019-05-07 2019-05-07 S-g滤波和自适应mp算法的低频振荡模态辨识方法

Publications (2)

Publication Number Publication Date
CN110048416A true CN110048416A (zh) 2019-07-23
CN110048416B CN110048416B (zh) 2021-03-02

Family

ID=67281163

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910376969.1A Active CN110048416B (zh) 2019-05-07 2019-05-07 S-g滤波和自适应mp算法的低频振荡模态辨识方法

Country Status (1)

Country Link
CN (1) CN110048416B (zh)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110554683A (zh) * 2019-09-09 2019-12-10 北京航天自动控制研究所 一种周期性控制过程中多模态自适应动态激励添加方法
CN111046327A (zh) * 2019-12-18 2020-04-21 河海大学 适用于低频振荡与次同步振荡辨识的Prony分析方法
CN111310303A (zh) * 2020-01-17 2020-06-19 合肥工业大学 一种幅值指数衰减的正弦波参数识别方法
CN114372344A (zh) * 2021-12-02 2022-04-19 上海电力大学 一种次同步振荡阻尼特性影响因素的量化辨识方法
CN114462451A (zh) * 2022-01-24 2022-05-10 北京航空航天大学 用于机械故障诊断的特征模态分解方法
CN115219787A (zh) * 2022-09-20 2022-10-21 湖南大学 基于改进矩阵束的电网相量移动测量方法、系统及介质
CN117290640A (zh) * 2023-11-27 2023-12-26 天津领语未来智能科技有限公司 用于非线性信号处理的奇异谱谐波分解方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100015938A1 (en) * 2008-07-17 2010-01-21 Panasonic Corporation Rf receiver
CN104077480A (zh) * 2014-06-27 2014-10-01 福州大学 基于Matrix Pencil的电力系统低频振荡模态辨识方法
CN106451498A (zh) * 2016-11-28 2017-02-22 福州大学 一种基于改进广义形态滤波的低频振荡模态辨识方法
CN106505587A (zh) * 2016-10-19 2017-03-15 福州大学 基于广义形态滤波与改进mp算法的低频振荡模态辨识方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100015938A1 (en) * 2008-07-17 2010-01-21 Panasonic Corporation Rf receiver
CN104077480A (zh) * 2014-06-27 2014-10-01 福州大学 基于Matrix Pencil的电力系统低频振荡模态辨识方法
CN106505587A (zh) * 2016-10-19 2017-03-15 福州大学 基于广义形态滤波与改进mp算法的低频振荡模态辨识方法
CN106451498A (zh) * 2016-11-28 2017-02-22 福州大学 一种基于改进广义形态滤波的低频振荡模态辨识方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
金涛 等: "基于广义形态滤波与改进矩阵束的电力系统低频振荡模态辨识", 《电工技术学报》 *

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110554683A (zh) * 2019-09-09 2019-12-10 北京航天自动控制研究所 一种周期性控制过程中多模态自适应动态激励添加方法
CN110554683B (zh) * 2019-09-09 2020-12-18 北京航天自动控制研究所 一种周期性控制过程中多模态自适应动态激励添加方法
CN111046327A (zh) * 2019-12-18 2020-04-21 河海大学 适用于低频振荡与次同步振荡辨识的Prony分析方法
CN111310303A (zh) * 2020-01-17 2020-06-19 合肥工业大学 一种幅值指数衰减的正弦波参数识别方法
CN111310303B (zh) * 2020-01-17 2024-02-02 合肥工业大学 一种幅值指数衰减的正弦波参数识别方法
CN114372344A (zh) * 2021-12-02 2022-04-19 上海电力大学 一种次同步振荡阻尼特性影响因素的量化辨识方法
CN114372344B (zh) * 2021-12-02 2024-09-20 上海电力大学 一种次同步振荡阻尼特性影响因素的量化辨识方法
CN114462451A (zh) * 2022-01-24 2022-05-10 北京航空航天大学 用于机械故障诊断的特征模态分解方法
CN114462451B (zh) * 2022-01-24 2024-06-07 北京航空航天大学 用于机械故障诊断的特征模态分解方法
CN115219787A (zh) * 2022-09-20 2022-10-21 湖南大学 基于改进矩阵束的电网相量移动测量方法、系统及介质
CN117290640A (zh) * 2023-11-27 2023-12-26 天津领语未来智能科技有限公司 用于非线性信号处理的奇异谱谐波分解方法
CN117290640B (zh) * 2023-11-27 2024-01-26 天津领语未来智能科技有限公司 用于非线性信号处理的奇异谱谐波分解方法

Also Published As

Publication number Publication date
CN110048416B (zh) 2021-03-02

Similar Documents

Publication Publication Date Title
CN110048416A (zh) S-g滤波和自适应mp算法的低频振荡模态辨识方法
CN103368175B (zh) 电力系统动态稳定在线评估方法
CN106845010B (zh) 基于改进SVD降噪和Prony的低频振荡主导模式辨识方法
CN104993480B (zh) 基于递推随机子空间的电力系统低频振荡在线辨识方法
CN103944174B (zh) 基于互相关函数滤噪算法的低频振荡在线辨识方法
CN103401238B (zh) 一种基于总体测辨法的电力负荷建模方法
CN103198184B (zh) 一种电力系统低频振荡特征类噪声辨识方法
CN103995178A (zh) 一种基于时频聚集特性准则s变换的电压暂降检测方法
CN106646121B (zh) 一种配电网故障行波波头的辨识方法
CN107658881A (zh) 基于戴维南等值方法的电压稳定临界点判断方法
CN107370150B (zh) 基于同步相量量测的电力系统状态估计不良数据处理方法
CN107860469A (zh) 一种基于正交多项式拟合的变电站噪声预测方法
CN101216512A (zh) 一种非正弦周期信号实时高精度检测方法
CN106451498A (zh) 一种基于改进广义形态滤波的低频振荡模态辨识方法
CN109218073A (zh) 一种计及网络攻击和参数不确定性的动态状态估计方法
CN105137373A (zh) 一种指数信号的去噪方法
Kravtsov Randomness, determinateness, and predictability
CN102567630B (zh) 一种大跨桥梁结构风致振动响应的确定方法
CN110470236B (zh) 一种嵌入光纤光栅的柔性结构形变重构方法
CN115688288A (zh) 飞行器气动参数辨识方法、装置、计算机设备及存储介质
CN103245830B (zh) 一种结合ar谱估计与非线性优化的间谐波检测方法
CN117309079A (zh) 基于时差法的超声飞渡时间测量方法、装置、设备及介质
CN112710925A (zh) 基于改进vmd和s变换的高渗透率主动配电网故障测距方法
Thiébaux Extending estimation accuracy with anisotropic interpolation
Chitturi et al. Comparing performance of Prony analysis and matrix pencil method for monitoring power system oscillations

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant