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

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

Info

Publication number
CN110048416B
CN110048416B CN201910376969.1A CN201910376969A CN110048416B CN 110048416 B CN110048416 B CN 110048416B CN 201910376969 A CN201910376969 A CN 201910376969A CN 110048416 B CN110048416 B CN 110048416B
Authority
CN
China
Prior art keywords
matrix
signal
frequency oscillation
low
filtering
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.)
Active
Application number
CN201910376969.1A
Other languages
English (en)
Other versions
CN110048416A (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

Images

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)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
  • Complex Calculations (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滤波预处理后的信号表示为
Figure GDA0002895354650000021
r(n)是残余噪声；
步骤S3：评估S-G滤波器的消噪效果；
步骤S4：提取并保存经S-G滤波预处理的低频振荡信号
Figure GDA0002895354650000031
n＝0,1,2,…,N-1；
步骤S5：设置矩阵束参数L，构建由预处理后的信号
Figure GDA0002895354650000032
构造的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个点组成一个局部点集，用多项式
Figure GDA0002895354650000033
对该点集进行拟合，其中pT(i)是进行拟合的多项式，i表示窗口中的第i个采样点，i＝-S…，-1,0,1…,S，S约为平滑窗口大小的一半，T为拟合多项式阶数，且T<2S+1，ak是多项式pT(i)的系数；
步骤S22：求解窗口内多项式拟合的误差，用总体误差平方和描述为：
Figure GDA0002895354650000041
步骤S23：使E最小，即
Figure GDA0002895354650000042
其中r＝0,1,2,…,T；
Figure GDA0002895354650000043
导出
Figure GDA0002895354650000044
令
Figure GDA0002895354650000045
则有
Figure GDA0002895354650000046
步骤S24：利用采样点数据y(i)，分别取r＝0,1,2,…,T，解方程组求得系数a0-aT；将各系数回代多项式
Figure GDA0002895354650000047
拟合出滤波后的函数值；逐点移动窗口，得出所有采样点滤波后的函数值，从而达到消噪的效果。
进一步，所述步骤S3中，采用信噪比指标评估S-G滤波器的消噪效果，所述信噪比指标表示为：
Figure GDA0002895354650000048
要求指标SNR＞20dB。
进一步，所述步骤S4中，用P阶幅值指数变化的余弦函数模型来表示低频振荡信号
Figure GDA0002895354650000049
该模型为：
Figure GDA00028953546500000410
其中n＝0,1,2,…,N-1，P是信号模态数，Ts是采样间隔，fi表示频率，αi表示阻尼因子，Ai表示初始幅值，θi表示初相位；扩展到复数域的讨论，该指数模型表示为：
Figure GDA0002895354650000051
其中
Figure GDA0002895354650000052
zi是信号的极点且
Figure GDA0002895354650000053
2P为复数域模态阶数，是相应实数域幅值指数变化的余弦函数模型中余弦函数个数的2倍。
进一步，所述步骤S5中，由预处理后的信号
Figure GDA0002895354650000054
构造的Hankel矩阵H如下：
Figure GDA0002895354650000055
进一步，所述步骤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个奇异值对应的累计占比：
Figure GDA0002895354650000056
其中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'，
Figure GDA0002895354650000061
其中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构成的对角阵，
Figure GDA0002895354650000062
Figure GDA0002895354650000063
zi是信号的极点，Bi表示各个模态的初始幅值，
Figure GDA0002895354650000064
Figure GDA0002895354650000065
进一步，所述步骤S7中，Bi的求取方法如下:
根据
Figure GDA0002895354650000066
r(n)作为残余噪声影响微弱忽略不计，因此可用
Figure GDA0002895354650000067
来估计
Figure GDA0002895354650000068
利用最小二乘法计算
Figure GDA0002895354650000069
其中
Figure GDA0002895354650000071
B＝[B1,B2,…,B2P]T，
Figure GDA0002895354650000072
进而求解各个Bi。
进一步，所述步骤S11中，判断估计信号的拟合精度是否满足要求，包括以下步骤：
步骤S111：将估计信号与预处理信号作对比，计算拟合精度AFI；
步骤S112：设定拟合精度阈值，判断拟合精度AFI是否满足拟合精度阈值要求，若满足拟合精度阈值要求，转入步骤S12，若不满足拟合精度阈值要求，返回步骤S5。
进一步，所述步骤S111中，拟合精度
Figure GDA0002895354650000073
其中，x和
Figure GDA0002895354650000074
分别表示低频振荡预处理信号向量和经模态辨识后估计信号向量，||.||表示2范数，AFI要求大于10dB。
与现有技术相比，本发明的优点如下：
(1)本发明通过S-G滤波器技术，对电力系统低频振荡信号进行消噪预处理，有效抑制噪声对后续模态辨识的影响。
(2)本发明利用所提出的定阶方法实现自适应定阶并能快速计算，避免凭人为经验设定阈值带来的主观因素影响，从而实时、准确辨识低频振荡主导模态。
附图说明
图1是本发明实施例的流程图。
图2是本发明实施例中的低频振荡理想信号x(n)。
图3是本发明实施例中的低频振荡实测信号y(n)。
图4是本发明实施例中的预处理低频振荡信号
Figure GDA0002895354650000082
图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滤波预处理后，得到
Figure GDA0002895354650000083
其中
Figure GDA0002895354650000085
是y(n)经S-G滤波预处理后的信号，r(n)是残余噪声。
Figure GDA0002895354650000084
的信号图如图4所示。
在本发明实施例中，所述步骤S2中，还包括以下步骤：
步骤S21：S-G滤波器利用移动的窗口进行信号平滑，设置窗口大小为2S+1，多项式阶数为T，在平滑窗口沿各采样点n＝0,1,2,…,N-1移动的过程中，分别将各采样点设置为窗口的中心，由各采样点为中心前后延伸S个点组成一个局部点集，用多项式
Figure GDA0002895354650000081
对该局部点集进行拟合，其中pT(i)是进行拟合的多项式，i表示窗口中的第i个采样点，i＝-S…，-1,0,1…,S，S约为平滑窗口大小的一半，T为拟合多项式阶数，且T<2S+1，ak是多项式pT(i)的系数。
步骤S22：求解窗口内多项式拟合的误差，用总体误差平方和描述为：
Figure GDA0002895354650000091
步骤S23：为使E最小，要求
Figure GDA0002895354650000092
其中r＝0,1,2,…,T；
Figure GDA0002895354650000093
可导出
Figure GDA0002895354650000094
令
Figure GDA0002895354650000095
Figure GDA0002895354650000096
则有
Figure GDA0002895354650000097
步骤S24：利用采样点数据y(i)，分别取r＝0,1,2,…,T，解方程组可求得系数a0-aT；将各系数回代多项式
Figure GDA0002895354650000098
拟合出滤波后的函数值；逐点移动窗口，得出所有采样点滤波后的函数值，从而达到消噪的效果。
步骤S3：采用信噪比(Signal Noise Ratio，SNR)指标评估S-G滤波器的消噪效果，判断经S-G滤波预处理后的信号
Figure GDA0002895354650000099
是否满足指标要求，若满足指标要求，则转入步骤S4；若不满足指标要求，返回步骤S2；
在本发明实施例中，在所述步骤S3中，信噪比指标：
Figure GDA00028953546500000910
要求指标SNR＞20dB。
步骤S4：提取并保存经S-G滤波预处理的低频振荡信号
Figure GDA00028953546500000911
n＝0,1,2,…,N-1；
所述步骤S4中，用P阶幅值指数变化的余弦函数模型来表示低频振荡信号
Figure GDA0002895354650000101
该模型为：
Figure GDA0002895354650000102
其中n＝0,1,2,…,N-1，P是信号模态数，Ts是采样间隔，fi、αi、Ai和θi分别表示待辨识各模态参数：频率、阻尼因子、初始幅值和初相位；扩展到复数域的讨论，该指数模型也可以表示为：
Figure GDA0002895354650000103
其中
Figure GDA0002895354650000104
zi是信号的极点且
Figure GDA0002895354650000105
2P为复数域模态阶数，是相应实数域幅值指数变化的余弦函数模型中余弦函数个数的2倍。
步骤S5：设置矩阵束参数L，L通常取N/4～N/3；构建由预处理后的信号
Figure GDA0002895354650000106
构造的Hankel矩阵H，由预处理后的信号
Figure GDA0002895354650000107
构造的Hankel矩阵H如下：
Figure GDA0002895354650000108
步骤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个奇异值对应的累计占比，例如：
Figure GDA0002895354650000111
其中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'，
Figure GDA0002895354650000112
其中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构成的对角阵，
Figure GDA0002895354650000121
Figure GDA0002895354650000122
zi是信号的极点，Bi表示各个模态的初始幅值，
Figure GDA0002895354650000123
Figure GDA0002895354650000124
步骤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]，
Figure GDA0002895354650000125
步骤S10：求解低频振荡信号各模态的初始幅值Ai和初相位θi：Ai＝2|Bi|，θi＝arctan[Im(Bi)/Re(Bi)]；
所述步骤S10中，Bi的求取方法如下:
根据
Figure GDA0002895354650000126
r(n)作为残余噪声影响微弱忽略不计，因此可用
Figure GDA0002895354650000127
来估计
Figure GDA0002895354650000128
利用最小二乘法计算
Figure GDA0002895354650000129
其中
Figure GDA00028953546500001210
B＝[B1,B2,…,B2P]T，
Figure GDA00028953546500001211
进而求解各个Bi。
步骤S11：根据上述求解的模态参数：频率fi、阻尼因子αi、初始幅值Ai和初相位θi拟合出估计信号，判断估计信号的拟合精度是否满足要求？若满足要求，转入步骤S12，若不满足要求，返回步骤S5；
所述步骤S11中，判断估计信号的拟合精度是否满足要求，包括以下步骤：
步骤S111：将估计信号与预处理信号作对比，计算拟合精度AFI；
步骤S112：设定拟合精度阈值，判断拟合精度AFI是否满足拟合精度阈值要求，若满足拟合精度阈值要求，转入步骤S12，若不满足拟合精度阈值要求，返回步骤S5。
所述步骤S111中，拟合精度
Figure GDA0002895354650000131
其中，x和
Figure GDA0002895354650000132
分别表示低频振荡预处理信号向量和经模态辨识后估计信号向量，||.||表示2范数，AFI要求大于10dB。
步骤S12：输出各模态参数值：频率fi、阻尼因子αi、初始幅值Ai和初相位θi。
较佳地，辨识后的各模态参数如表1所示，低频振荡预处理信号和辨识拟合信号如图6所示。
表1辨识后各模态参数表
Figure GDA0002895354650000133
本发明首先利用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滤波预处理后的信号表示为
Figure FDA0002895354640000011
r(n)是残余噪声；
步骤S3：评估S-G滤波器的消噪效果；
步骤S4：提取并保存经S-G滤波预处理的低频振荡信号
Figure FDA0002895354640000012
Figure FDA0002895354640000013
步骤S5：设置矩阵束参数L，构建由预处理后的信号
Figure FDA0002895354640000014
构造的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个点组成一个局部点集，用多项式
Figure FDA0002895354640000021
对该点集进行拟合，其中pT(i)是进行拟合的多项式，i表示窗口中的第i个采样点，i＝-S…，-1,0,1…,S，S约为平滑窗口大小的一半，T为拟合多项式阶数，且T<2S+1，ak是多项式pT(i)的系数；
步骤S22：求解窗口内多项式拟合的误差，用总体误差平方和描述为：
Figure FDA0002895354640000022
步骤S23：使E最小，即
Figure FDA0002895354640000023
其中r＝0,1,2,…,T；
Figure FDA0002895354640000024
导出
Figure FDA0002895354640000025
令
Figure FDA0002895354640000026
则有
Figure FDA0002895354640000027
步骤S24：利用采样点数据y(i)，分别取r＝0,1,2,…,T，解方程组求得系数a0-aT；将各系数回代多项式
Figure FDA0002895354640000028
拟合出滤波后的函数值；逐点移动窗口，得出所有采样点滤波后的函数值，从而达到消噪的效果。
3.如权利要求1所述的S-G滤波和自适应MP算法的低频振荡模态辨识方法，其特征在于：所述步骤S3中，采用信噪比指标评估S-G滤波器的消噪效果，所述信噪比指标表示为：
Figure FDA0002895354640000031
要求指标SNR＞20dB。
4.如权利要求1所述的S-G滤波和自适应MP算法的低频振荡模态辨识方法，其特征在于：所述步骤S4中，用P阶幅值指数变化的余弦函数模型来表示低频振荡信号
Figure FDA0002895354640000032
该模型为：
Figure FDA0002895354640000033
其中n＝0,1,2,…,N-1，P是信号模态数，Ts是采样间隔，fi表示频率，αi表示阻尼因子，Ai表示初始幅值，θi表示初相位；扩展到复数域的讨论，该指数模型表示为：
Figure FDA0002895354640000034
其中
Figure FDA0002895354640000035
zi是信号的极点且
Figure FDA0002895354640000036
2P为复数域模态阶数，是相应实数域幅值指数变化的余弦函数模型中余弦函数个数的2倍。
5.如权利要求1所述的S-G滤波和自适应MP算法的低频振荡模态辨识方法，其特征在于：所述步骤S5中，由预处理后的信号
Figure FDA0002895354640000037
构造Hankel矩阵H如下：
Figure FDA0002895354640000041
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个奇异值对应的累计占比：
Figure FDA0002895354640000042
其中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'，
Figure FDA0002895354640000051
其中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构成的对角阵，
Figure FDA0002895354640000052
Figure FDA0002895354640000053
zi是信号的极点，Bi是各个模态的初始幅值，
Figure FDA0002895354640000054
Figure FDA0002895354640000055
8.如权利要求7所述的S-G滤波和自适应MP算法的低频振荡模态辨识方法，其特征在于：所述步骤S7中，Bi的求取方法如下:
根据
Figure FDA0002895354640000056
r(n)作为残余噪声影响微弱忽略不计，因此可用
Figure FDA0002895354640000057
来估计
Figure FDA0002895354640000058
利用最小二乘法计算
Figure FDA0002895354640000059
其中
Figure FDA0002895354640000061
B＝[B1,B2,…,B2P]T，
Figure FDA0002895354640000062
进而求解各个Bi。
9.如权利要求1所述的S-G滤波和自适应MP算法的低频振荡模态辨识方法，其特征在于：所述步骤S11中，判断估计信号的拟合精度是否满足要求，包括以下步骤：
步骤S111：将估计信号与预处理信号作对比，计算拟合精度AFI；
步骤S112：设定拟合精度阈值，判断拟合精度AFI是否满足拟合精度阈值要求，若满足拟合精度阈值要求，转入步骤S12，若不满足拟合精度阈值要求，返回步骤S5。
10.如权利要求9所述的S-G滤波和自适应MP算法的低频振荡模态辨识方法，其特征在于：所述步骤S111中，拟合精度
Figure FDA0002895354640000063
其中，x和
Figure FDA0002895354640000064
分别表示低频振荡预处理信号向量和经模态辨识后估计信号向量，||.||表示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 CN110048416A (zh) 2019-07-23
CN110048416B true 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)

Families Citing this family (5)

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

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP5067294B2 (ja) * 2008-07-17 2012-11-07 パナソニック株式会社 電波受信装置
CN104077480A (zh) * 2014-06-27 2014-10-01 福州大学 基于Matrix Pencil的电力系统低频振荡模态辨识方法
CN106505587B (zh) * 2016-10-19 2019-01-18 福州大学 基于广义形态滤波与改进mp算法的低频振荡模态辨识方法
CN106451498B (zh) * 2016-11-28 2019-01-18 福州大学 一种基于改进广义形态滤波的低频振荡模态辨识方法

Also Published As

Publication number Publication date
CN110048416A (zh) 2019-07-23

Similar Documents

Publication Publication Date Title
CN110048416B (zh) S-g滤波和自适应mp算法的低频振荡模态辨识方法
CN105699804B (zh) 一种配电网大数据故障检测与定位方法
Grant et al. Comparison of matrix pencil and prony methods for power system modal analysis of noisy signals
CN106845010B (zh) 基于改进SVD降噪和Prony的低频振荡主导模式辨识方法
CN105374367B (zh) 异常帧检测方法和装置
CN104375976B (zh) 基于张量正则分解的欠定盲源分离中的混合矩阵识别方法
CN110967599A (zh) 一种电能质量扰动检测与定位算法
Angrisani et al. Ultrasonic time-of-flight estimation through unscented Kalman filter
CN103439688A (zh) 一种用于分布式麦克风阵列的声源定位系统及定位方法
Belega et al. Accuracy of the normalized frequency estimation of a discrete-time sine-wave by the energy-based method
CN114218778A (zh) 一种用于声爆试验数据的分析方法及装置
CN108281961B (zh) 一种自适应鲁棒扩展卡尔曼的参数辨识方法
CN106972949B (zh) 一种基于自适应补偿技术的分数阶网络系统状态估计方法
Wahlberg et al. Identification of Wiener systems with process noise is a nonlinear errors-in-variables problem
CN112688324B (zh) 基于FastICA与TLS-ESPRIT的电力系统低频振荡模态辨识方法
Angrisani et al. New proposal for uncertainty evaluation in indirect measurements
So et al. Efficient frequency estimation of a single real tone based on principal singular value decomposition
Qing-Fang et al. Small-time scale network traffic prediction based on a local support vector machine regression model
Mehra et al. Modes preserving wavelet based multi-scale PCA algorithm for compression of smart grid data
CN116561541A (zh) 随机动载荷识别模型、模型的训练方法以及识别方法
Orosz et al. Analysis of resonator-based harmonic estimation in the case of data loss
Xianmin A new method with high confidence for validation of computer simulation models of flight systems
CN104868876B (zh) 一种针对过程噪声协方差矩阵Q未知情况下的Kalman滤波方法
CN111342478B (zh) 一种基于最优变量投影的电力系统动态稳定评估方法
CN110032758B (zh) 计算电信号的能量的方法、装置和计算机存储介质

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