CN115032453A - 一种多频动态相量测量方法 - Google Patents

一种多频动态相量测量方法 Download PDF

Info

Publication number
CN115032453A
CN115032453A CN202210953016.9A CN202210953016A CN115032453A CN 115032453 A CN115032453 A CN 115032453A CN 202210953016 A CN202210953016 A CN 202210953016A CN 115032453 A CN115032453 A CN 115032453A
Authority
CN
China
Prior art keywords
frequency
signal
dynamic
phasor
equation
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.)
Pending
Application number
CN202210953016.9A
Other languages
English (en)
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.)
Sichuan University
Original Assignee
Sichuan University
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 Sichuan University filed Critical Sichuan University
Priority to CN202210953016.9A priority Critical patent/CN115032453A/zh
Publication of CN115032453A publication Critical patent/CN115032453A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R25/00—Arrangements for measuring phase angle between a voltage and a current or between voltages or currents
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
    • G01R23/16—Spectrum analysis; Fourier analysis
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y04—INFORMATION OR COMMUNICATION TECHNOLOGIES HAVING AN IMPACT ON OTHER TECHNOLOGY AREAS
    • Y04S—SYSTEMS INTEGRATING TECHNOLOGIES RELATED TO POWER NETWORK OPERATION, COMMUNICATION OR INFORMATION TECHNOLOGIES FOR IMPROVING THE ELECTRICAL POWER GENERATION, TRANSMISSION, DISTRIBUTION, MANAGEMENT OR USAGE, i.e. SMART GRIDS
    • Y04S10/00—Systems supporting electrical power generation, transmission or distribution
    • Y04S10/22—Flexible AC transmission systems [FACTS] or power factor or reactive power compensating or correcting units

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本发明公开了一种多频动态相量测量方法，包括以下步骤：建立多频动态信号的TFM模型和得出多频动态信号中DDC的近似表征；建立动态相量估计模型和重构动态相量方程；动态相量的设计步骤：建立基本性能测试场景、频率偏差测试场景、谐波震荡测试场景和抗干扰能力测试场景。本发明实现了足够的测量精度，在多频相位分析能力中减少了谐波间谐波分量的时变性及其对动态相量测量的影响，显著提高了测量精度。

Description

一种多频动态相量测量方法
技术领域
本发明涉及电子信息领域，特别是涉及一种多频动态相量测量方法。
背景技术
随着电网的监测技术的不断进步，以及相量测量单元(PMUs)等量测设备在电网中的投入使用，电力系统的动态监控成为了电网可靠运行和控制的重要保证。其中，谐波、间谐波信息可以用于保护和监控电力系统应用，如高阻断故障识别、智能孤岛检测、谐波状态估计等。电网谐波及间谐波也会破坏电力系统的安全、稳定和经济运行。另一方面，电力电子设备等非线性负荷的持续投入导致大量谐波注入电网，电压波形严重“畸变”，电能质量下降。因此，快速准确地检测电网中的谐波和间谐波分量具有重要的工程实用价值。
由于谐波和间谐波广泛存在于输配电网络中，越来越多的相量测量单元试图兼备基本相量测量和谐波相量测量功能。但二者独特的约束条件和成本效益要求亟需开发一种兼顾基本相量测量和谐波相量测量的新算法。目前，有大量关于相量测量的文献，而精度一定程度上取决于算法中测量波形的建模方式。
为了应对配电网全新测量环境的挑战，各国学者提出了很多改进型或者全新的PMU测量算法，主要可以分为两类：基于模型的相量测量方法和基于离散傅里叶变换（DFT）的方法，第一种方法将动态相量近似用不同数学模型展开，其中包括如泰勒傅里叶变换、递归最小二乘法(RLS)方法、自回归移动平均(ARMA)方法以及基于sinc插值函数的谐波相量测量算法。
现阶段动态谐波分析方法中，为了解决间谐波对估计器精度的干扰，减少相量解的计算负担，近年来机器学习领域的研究将传统相量估计方法与压缩感知算法相结合，可以使用较少的数据点精确恢复特定信号，同时保持相对较短的观测间隔并抵抗间谐波干扰，避免了大量冗余数据对运算速度的影响，有学者研究了基于压缩感知（CS）的超分辨率离散傅里叶变换（DFT）分析，其中较短的采样序列可用于PMU和谐波分析的应用。然而，它只能用于静态相量模型。
发明内容
基于此，针对上述技术问题，本发明提出了一种多频动态相量测量方法，旨在精确测量全TFM波形模型中的谐波和相间分量。
本发明的技术方案是：一种多频动态相量测量方法，包括以下步骤：
a、DDC分量的动态相量估计步骤：电力系统的多频动态信号模型由振幅和相位时变的正弦分量和表示：
Figure 94823DEST_PATH_IMAGE001
(1)
其中，其中
Figure 302076DEST_PATH_IMAGE002
为信号中的DDC分量，
Figure 371663DEST_PATH_IMAGE003
为基波及谐波信号，
Figure 162901DEST_PATH_IMAGE004
和
Figure 976136DEST_PATH_IMAGE005
分别为DDC分量的幅值和时间常数，
Figure 536431DEST_PATH_IMAGE006
和
Figure 635974DEST_PATH_IMAGE007
分别为谐波相量的幅值和时间常数，
Figure 413044DEST_PATH_IMAGE008
是包含实际电力系统频率
Figure 498811DEST_PATH_IMAGE009
，其谐波倍数和可能的间谐波频率的通用频率集合，然后建立多频动态信号的TFM模型和得出多频动态信号中DDC的近似表征；
b、动态相量的估计步骤：设定
Figure 444771DEST_PATH_IMAGE010
为原始信号，
Figure 856160DEST_PATH_IMAGE011
为目标信号，计算其
Figure 621991DEST_PATH_IMAGE012
个线性测量值来获得观测数据，用以下公式表示:
Figure 901663DEST_PATH_IMAGE013
(14)
其中，
Figure 577495DEST_PATH_IMAGE014
为
Figure 520305DEST_PATH_IMAGE015
矩阵，
Figure 914377DEST_PATH_IMAGE016
为观测值，矩阵
Figure 466581DEST_PATH_IMAGE014
从
Figure 387133DEST_PATH_IMAGE017
映射到
Figure 874746DEST_PATH_IMAGE018
，它表示因为
Figure 880748DEST_PATH_IMAGE019
而产生的降维，假设信号
Figure 377589DEST_PATH_IMAGE020
是可通过正交变换
Figure 668760DEST_PATH_IMAGE021
进行压缩且
Figure 327274DEST_PATH_IMAGE022
，然后将公式(1)重写为：
Figure 554993DEST_PATH_IMAGE023
(15)
其中，
Figure 589945DEST_PATH_IMAGE024
是原始信号
Figure 485089DEST_PATH_IMAGE020
的变换系数，采用高斯随机矩阵
Figure 580084DEST_PATH_IMAGE014
和离散小波逆变换
Figure 295099DEST_PATH_IMAGE025
，且
Figure 133742DEST_PATH_IMAGE026
的组合矩阵也是高斯的并且概率满足RIP条件，在
Figure 384857DEST_PATH_IMAGE027
足够高的情况下，重建系数
Figure 509808DEST_PATH_IMAGE024
，对原信号
Figure 587485DEST_PATH_IMAGE020
进行逆变换
Figure 88874DEST_PATH_IMAGE028
即可恢复原信号
Figure 568397DEST_PATH_IMAGE020
，如果确定原信号
Figure 129828DEST_PATH_IMAGE020
，就通过公式计算谐波的频率
Figure 694802DEST_PATH_IMAGE029
、振幅
Figure 498416DEST_PATH_IMAGE030
以及频率变化率
Figure 832446DEST_PATH_IMAGE031
；然后建立动态相量估计模型和重构动态相量方程；
c、动态相量的设计步骤：建立基本性能测试场景、频率偏差测试场景、谐波震荡测试场景和抗干扰能力测试场景。
步骤a中，“建立多频动态信号的TFM模型”的步骤如下：
将公式(1)中
Figure 564778DEST_PATH_IMAGE032
的每个分量与一个动态相量相关联，定义为：
Figure 351469DEST_PATH_IMAGE033
Figure 460239DEST_PATH_IMAGE034
(2)
将公式（1）中的
Figure 648775DEST_PATH_IMAGE035
方程改写为：
Figure 552009DEST_PATH_IMAGE036
(3)
基于电力系统的大惯性特征，将用泰勒级数近似表示为：
Figure 91574DEST_PATH_IMAGE037
(4)
其中，
Figure 239921DEST_PATH_IMAGE038
为
Figure 548543DEST_PATH_IMAGE039
的
Figure 357099DEST_PATH_IMAGE040
阶导数，
Figure 508595DEST_PATH_IMAGE040
为泰勒展开阶数，K取3，即：
Figure 568954DEST_PATH_IMAGE041
(5)
所以，
Figure 856716DEST_PATH_IMAGE042
的离散表达式为:
Figure 977119DEST_PATH_IMAGE043
(6)
其中，T为采样间隔，∗表示共轭计算符；
将公式(2)代入公式(6)得信号的泰勒展开表达式，再以采样频率
Figure 854726DEST_PATH_IMAGE044
进行采样.设
Figure 718777DEST_PATH_IMAGE045
是一个长度为
Figure 126625DEST_PATH_IMAGE046
的有限样本序列，设
Figure 417929DEST_PATH_IMAGE046
为偶数且
Figure 12858DEST_PATH_IMAGE047
，采样间隔
Figure 805233DEST_PATH_IMAGE048
，因此，动态相量计算的时间参考位于样本记录中的
Figure 942954DEST_PATH_IMAGE049
，离散化后的信号表达式
Figure 31258DEST_PATH_IMAGE050
为：
Figure 254429DEST_PATH_IMAGE051
(7)
其中，
Figure 584916DEST_PATH_IMAGE014
是原式中
Figure 577143DEST_PATH_IMAGE052
的系数，表示大小为
Figure 600462DEST_PATH_IMAGE053
的矩阵，该矩阵为指数项的泰勒傅里叶基矩阵，为防止与式(2)中同类字母相量混淆，另设
Figure 45350DEST_PATH_IMAGE054
，作为长度为
Figure 913949DEST_PATH_IMAGE055
的列向量，当谐波
Figure 26262DEST_PATH_IMAGE056
时
Figure 453438DEST_PATH_IMAGE057
的集合，∗表示共轭计算符，即
Figure 651202DEST_PATH_IMAGE058
。
步骤a中，“得出多频动态信号中DDC的近似表征”的步骤如下：
将较窄时间窗内
Figure 57912DEST_PATH_IMAGE059
的衰减直流分量近似表示成一个动态、频率较低的余弦分量，即:
Figure 149365DEST_PATH_IMAGE060
（8）
其中，
Figure 655433DEST_PATH_IMAGE061
和
Figure 933967DEST_PATH_IMAGE062
分别为模型DDC分量的振幅和初始相位，
Figure 285314DEST_PATH_IMAGE063
为算法观测时间窗长度，
Figure 732738DEST_PATH_IMAGE064
为一个较低的频率，当DDC时间常数很大时，DDC分量倾向于直流分量，在此条件下，
Figure 144128DEST_PATH_IMAGE065
Hz；
接下来，基于频域采样定理，将DDC频率余弦分量对应的动态相量记为
Figure 175538DEST_PATH_IMAGE066
，对幅值有限的时域信号
Figure 799417DEST_PATH_IMAGE067
，基于频域采样定理对其进行参数化建模，具体表示为:
Figure 865462DEST_PATH_IMAGE068
（9）
其中，
Figure 447753DEST_PATH_IMAGE069
为相量
Figure 700880DEST_PATH_IMAGE070
在频率
Figure 11339DEST_PATH_IMAGE071
处的频率采样值，
Figure 541678DEST_PATH_IMAGE072
表示一个向下取整算子，
Figure 419504DEST_PATH_IMAGE073
是频域采样间隔，
Figure 35293DEST_PATH_IMAGE074
代表用来进行
Figure 656767DEST_PATH_IMAGE070
参数化建模的频域采样点数，为提高上式模型的准确度，要求
Figure 307191DEST_PATH_IMAGE075
，对
Figure 90340DEST_PATH_IMAGE070
进行建模后，实现动态DDC余弦分量和动态基波分量的近似表征；
将
Figure 459004DEST_PATH_IMAGE076
设为一个奇数，使得
Figure 120055DEST_PATH_IMAGE077
时刻在观测窗口的中心，公式(8)中拟合的DDC分量的离散形式用公式(10)表示：
Figure 749619DEST_PATH_IMAGE078
(10)
其中，
Figure 579035DEST_PATH_IMAGE079
为公式(8)右侧表达式中余弦信号的
Figure 559629DEST_PATH_IMAGE080
个样本点的列向量，式中
Figure 398272DEST_PATH_IMAGE081
为一个包含了信号
Figure 616764DEST_PATH_IMAGE082
的
Figure 882660DEST_PATH_IMAGE080
个采样点的列向量，
Figure 583507DEST_PATH_IMAGE083
是一个包含了
Figure 225841DEST_PATH_IMAGE084
的行向量，
Figure 564418DEST_PATH_IMAGE085
为一个每列包含了
Figure 125849DEST_PATH_IMAGE080
个
Figure 425244DEST_PATH_IMAGE086
采样点的矩阵,为共轭算子；
公式(10)的最小二乘技术提供了最好的参数，因为它获得了
Figure 730323DEST_PATH_IMAGE083
和二阶泰勒逼近之间的最小误差.此时是满足下列约束条件下的最优解：
Figure 64352DEST_PATH_IMAGE087
(11)
其中，
Figure 298150DEST_PATH_IMAGE088
表示欧几里得范数，再引入拉格朗日算子和埃尔米特运算符进行求导，相量的系数向量计算为：
Figure 350420DEST_PATH_IMAGE089
(12)
式中H表示埃尔米特算子，根据公式(9)和公式(12)，得到DDC分量的估计值为:
Figure 928031DEST_PATH_IMAGE090
(13)。
步骤b中，“建立动态相量估计模型”的步骤如下：
将测量结果还原为原始信号，则重建信号的目标函数改为：
Figure 506780DEST_PATH_IMAGE091
(16)
其中，
Figure 285380DEST_PATH_IMAGE092
是任意向量
Figure 949580DEST_PATH_IMAGE093
的
Figure 471828DEST_PATH_IMAGE094
范数，其
Figure 409478DEST_PATH_IMAGE094
范数计算定义为
Figure 358980DEST_PATH_IMAGE095
；
公式(16)的正则化引入辅助信号
Figure 979317DEST_PATH_IMAGE096
，得到一个新的目标函数
Figure 898731DEST_PATH_IMAGE097
：
Figure 61859DEST_PATH_IMAGE098
(17)
Figure 306896DEST_PATH_IMAGE099
(18)
Figure 821054DEST_PATH_IMAGE100
(19)
其中，
Figure 45624DEST_PATH_IMAGE101
是梯度算子，
Figure 328838DEST_PATH_IMAGE102
是是垂直于
Figure 10355DEST_PATH_IMAGE096
的梯度的单位向量，
Figure 746230DEST_PATH_IMAGE103
是信号
Figure 600922DEST_PATH_IMAGE104
，所以
Figure 830653DEST_PATH_IMAGE105
是
Figure 558437DEST_PATH_IMAGE106
的梯度，
Figure 906242DEST_PATH_IMAGE107
,
Figure 112096DEST_PATH_IMAGE108
是垂直方向和水平方向的偏导数；
新的正则化是参考信号特征的方向向量
Figure 228956DEST_PATH_IMAGE102
和目标信号的梯度向量
Figure 127642DEST_PATH_IMAGE109
的点积，如果
Figure 697164DEST_PATH_IMAGE096
和
Figure 706708DEST_PATH_IMAGE106
边缘方向夹脚太小，所以优化正则化表达式
Figure 445119DEST_PATH_IMAGE110
得到：
Figure 249127DEST_PATH_IMAGE111
(20)
其中，
Figure 571524DEST_PATH_IMAGE112
是
Figure 853601DEST_PATH_IMAGE093
方向上的差算子,
Figure 945054DEST_PATH_IMAGE113
是
Figure 919963DEST_PATH_IMAGE114
方向上的差算子；
新的目标函数(18)有两个正则项:一个是系数
Figure 729656DEST_PATH_IMAGE115
的
Figure 81003DEST_PATH_IMAGE116
正则化，另一个是受参考信号
Figure 277496DEST_PATH_IMAGE117
约束的二次正则化，
Figure 813520DEST_PATH_IMAGE118
包含目标信号
Figure 720296DEST_PATH_IMAGE106
和参考信号
Figure 734388DEST_PATH_IMAGE096
的边缘方向，
Figure 675799DEST_PATH_IMAGE119
指向新正则化的方向，
Figure 851566DEST_PATH_IMAGE120
控制新正则化的强度，如果
Figure 245638DEST_PATH_IMAGE096
的梯度较小，而
Figure 564886DEST_PATH_IMAGE102
方向上的
Figure 360804DEST_PATH_IMAGE121
投影较大，则执行正则化惩罚，向量
Figure 238630DEST_PATH_IMAGE122
控制平滑的方向，
Figure 119998DEST_PATH_IMAGE123
控制平滑的强度.混合正则化提高压缩感知重构的性能。
步骤b中，“重构动态相量方程”的步骤如下：
引入一个辅助变量，多次迭代求解布雷格曼模型来减小信号结构损失，新的目标函数中，基于分裂的准则，引入辅助变量
Figure 475893DEST_PATH_IMAGE124
、
Figure 126317DEST_PATH_IMAGE125
，对于这两个约束，
Figure 175045DEST_PATH_IMAGE126
被新的正则化参数
Figure 543709DEST_PATH_IMAGE127
和
Figure 936251DEST_PATH_IMAGE128
取代，用来控制二次惩罚函数项.并利用二次惩罚项转化为一个无约束的最优化分裂模型：
Figure 565816DEST_PATH_IMAGE129
(21)
基于布雷格曼分割迭代，计算
Figure 660811DEST_PATH_IMAGE115
的目标函数，迭代过程被分解为几个步骤，信号重建通过求解下述优化问题来实现：
Figure 641405DEST_PATH_IMAGE130
(22)
Figure 214469DEST_PATH_IMAGE131
(23)
Figure 698540DEST_PATH_IMAGE132
(24)
在公式(22)式的迭代过程中,
Figure 964436DEST_PATH_IMAGE133
单调递减,直至为0，其次,在公式(22)的迭代求解中,只要满足
Figure 668212DEST_PATH_IMAGE134
，
Figure 310546DEST_PATH_IMAGE135
单调趋近真实解
Figure 914703DEST_PATH_IMAGE136
；
在从公式(22)到(24)的每个迭代周期中，公式(22)是一个可微优化问题，公式(24)显式求解并直接更新，利用定义
Figure 351500DEST_PATH_IMAGE137
函数来对变量
Figure 41107DEST_PATH_IMAGE138
进行更新，求解式(14)，即：
Figure 955974DEST_PATH_IMAGE140
(25)
其中，
Figure 680216DEST_PATH_IMAGE141
函数定义为：
Figure 287915DEST_PATH_IMAGE143
(26)
在公式(22)中，
Figure 703634DEST_PATH_IMAGE144
，
Figure 546825DEST_PATH_IMAGE145
,
Figure 940DEST_PATH_IMAGE146
,
Figure 904174DEST_PATH_IMAGE147
和
Figure 443740DEST_PATH_IMAGE148
都视为常数矩阵和向量，
Figure 90622DEST_PATH_IMAGE135
是可变的，最小化(22)，在第
Figure 399243DEST_PATH_IMAGE149
次迭代中，将对
Figure 709264DEST_PATH_IMAGE150
的一阶导数设为0，求出
Figure 470547DEST_PATH_IMAGE151
相对于
Figure 921120DEST_PATH_IMAGE135
的一阶导数，然后得到：
Figure 84248DEST_PATH_IMAGE153
(27)
基于公式(27)，定义线性正则化算子
Figure 329284DEST_PATH_IMAGE154
为：
Figure 843442DEST_PATH_IMAGE156
(28)
根据公式(28)中
Figure 832127DEST_PATH_IMAGE154
的定义和函数(18)的一阶导数得到：
Figure 849761DEST_PATH_IMAGE158
(29)
通过将公式(20)用
Figure 764234DEST_PATH_IMAGE159
表示，得到了
Figure 765688DEST_PATH_IMAGE159
的一个封闭解：
Figure 26906DEST_PATH_IMAGE161
(30)
其中，
Figure 289260DEST_PATH_IMAGE162
是一个单位矩阵，公式(21)是线性的，求解
Figure 17044DEST_PATH_IMAGE159
，将公式(23)、(24)、(30)循环合并更新，当正则化参数
Figure 99270DEST_PATH_IMAGE163
时，计算系数
Figure 570702DEST_PATH_IMAGE135
，成为基本优化；
解得系数
Figure 189028DEST_PATH_IMAGE164
后通过
Figure 822134DEST_PATH_IMAGE165
可得到重构相量
Figure 657235DEST_PATH_IMAGE166
（即
Figure 666780DEST_PATH_IMAGE167
），其每一行对应于每个频率分量在不同时间的相量，引入插值因子F，频率分辨率为
Figure 638147DEST_PATH_IMAGE168
，重构相量频率表示为
Figure 707734DEST_PATH_IMAGE169
，
Figure 764551DEST_PATH_IMAGE170
为频率索引，确定
Figure 577787DEST_PATH_IMAGE166
，所求的谐波的频率
Figure 896336DEST_PATH_IMAGE171
、振幅
Figure 136824DEST_PATH_IMAGE172
以及频率变化率
Figure 946518DEST_PATH_IMAGE173
通过公式(31-33)进行计算：
Figure 32285DEST_PATH_IMAGE175
(31)
Figure 978244DEST_PATH_IMAGE177
(32)
Figure 124055DEST_PATH_IMAGE179
(33)
其中，
Figure 155465DEST_PATH_IMAGE180
为
Figure 310503DEST_PATH_IMAGE166
的第列，
Figure 346854DEST_PATH_IMAGE181
，
Figure 53779DEST_PATH_IMAGE182
和
Figure 447851DEST_PATH_IMAGE183
分别为
Figure 265635DEST_PATH_IMAGE180
的零阶导数，一阶导数和二阶导数。
步骤c中，“建立基本性能测试场景”的步骤如下：
设定采样窗长为5个周期，基频带宽为1Hz，构造式(34)所示含基波及动态分量的信号模型：
Figure 61552DEST_PATH_IMAGE185
(33)
公式（33）中，f1为基波频率，此处设为50Hz，
Figure 939378DEST_PATH_IMAGE186
、
Figure 555167DEST_PATH_IMAGE187
分别代表基波和各次谐波相角，取
Figure 675177DEST_PATH_IMAGE188
范围中的任意值。低频段谐波次数h取值为2-13，采样频率设置为5kHz；
Figure 60022DEST_PATH_IMAGE189
,
Figure 108749DEST_PATH_IMAGE190
均初始化为全零矩阵，正则化参数
Figure 477414DEST_PATH_IMAGE191
，设定
Figure 637000DEST_PATH_IMAGE192
，
Figure 985DEST_PATH_IMAGE193
为参数间平衡系数，在
Figure 95980DEST_PATH_IMAGE194
内算法稳定性较高，k为迭代次数，在参数
Figure 312460DEST_PATH_IMAGE193
、k变化情况下，对重构效果及算法运行时间进行分析。
步骤c中，“建立频率偏差测试场景”的步骤如下：
基频带宽设置为1Hz，采用
Figure 151103DEST_PATH_IMAGE195
= 1.15Hz，获得每个谐波频率附近良好的通带和阻带性能，采样频率设置为10kHz，多频动态信号如下所示:
Figure 900753DEST_PATH_IMAGE197
(34)
式中，
Figure 901070DEST_PATH_IMAGE198
的变化范围为49.5-50.5 Hz，阶跃为0.2Hz。
步骤c中，“建立谐波震荡测试场景”的步骤如下：
建立多频动态信号如下：
Figure 103381DEST_PATH_IMAGE200
(35)
采样率设定为5kHz，采样周期长度为5个周期，其中
Figure 480136DEST_PATH_IMAGE201
为调制频率，设置为5Hz。
步骤c中，“建立抗干扰能力测试场景”的步骤如下：
信号引入信噪比为60dB的高斯白噪声，具体的动态信号如下：
Figure 84293DEST_PATH_IMAGE203
(36)
其中
Figure 521090DEST_PATH_IMAGE204
为间谐波频率，采样率设定为5kHz，采样周期长度为5个周期。
本发明的有益效果是：
a、首先，基于泰勒-傅里叶多频变换建立可以同时估计谐波分量和间谐波分量的动态信号模型；其次，基于动态信号分布的规律性和稀疏性，引入辅助信号将相量问题转化为压缩感知问题，最后，将动态相量重构问题转化为分裂布雷格曼迭代对混合正则化算法的优化问题，对信号进行重构，得到动态相量估计；
b、本发明还考虑了高阶导数，允许同时估计相量，这使得相量估计更加准确和高效，因此，实现了足够的测量精度，在多频相位分析能力中减少了谐波间谐波分量的时变性及其对动态相量测量的影响，显著提高了测量精度。
附图说明
图1为重构效果及运行时间随参数变化情况图；
图2为第2-13次谐波的最大TVEs、FEs和Fes图；
图3为第2-13次谐波的参数估计结果图；
图4[a]为间谐波和噪声干扰条件下各方法的估计结果图；
图4[b]为间谐波和噪声干扰条件下各方法的估计结果图；
图4[c]为间谐波和噪声干扰条件下各方法的估计结果图；
图5为本发明将动态谐波频域分布的规律性作为动态相量恢复的优化目标的算法流程图。
具体实施方式
下面结合附图对本发明的实施例进行详细说明。
实施例：
1、考虑DDC分量的动态相量估计
对含有基波、动态谐波、间谐波和DDC分量的信号来说，电力系统的多频动态信号模型
Figure 715092DEST_PATH_IMAGE205
可以由振幅和相位时变的正弦分量和表示：
Figure 895538DEST_PATH_IMAGE206
(1)
其中，其中
Figure 619780DEST_PATH_IMAGE207
为信号中的DDC分量，
Figure 227479DEST_PATH_IMAGE208
为基波及谐波信号，
Figure 138803DEST_PATH_IMAGE209
和
Figure 122940DEST_PATH_IMAGE210
分别为DDC分量的幅值和时间常数，
Figure 436109DEST_PATH_IMAGE211
和
Figure 840808DEST_PATH_IMAGE212
分别为谐波相量的幅值和时间常数，
Figure 114795DEST_PATH_IMAGE213
是包含实际电力系统频率
Figure 761677DEST_PATH_IMAGE214
，其谐波倍数和可能的间谐波频率的通用频率集合。
1.1、多频动态信号的TFM模型
(1)中
Figure 804719DEST_PATH_IMAGE215
的每个分量可以与一个动态相量相关联，定义为 ：
Figure 878854DEST_PATH_IMAGE216
(2)
所以（1）中的
Figure 905716DEST_PATH_IMAGE217
方程可改写为：
Figure 90710DEST_PATH_IMAGE218
(3)
基于电力系统的大惯性特征，可将用泰勒级数近似表示为：
Figure 519417DEST_PATH_IMAGE219
(4)
其中，
Figure 997409DEST_PATH_IMAGE220
为
Figure 511567DEST_PATH_IMAGE221
的
Figure 234673DEST_PATH_IMAGE222
阶导数，
Figure 642520DEST_PATH_IMAGE223
为泰勒展开阶数，综合考虑模型精确度和算法运算量等因素，本发明中K取3，即：
Figure 933824DEST_PATH_IMAGE224
(5)
所以，
Figure 794333DEST_PATH_IMAGE225
的离散表达式为:
Figure 462075DEST_PATH_IMAGE226
(6)
其中，T为采样间隔，∗表示共轭计算符。
将公式(2)代入(6)可得信号的泰勒展开表达式，再以采样频率
Figure 225894DEST_PATH_IMAGE227
进行采样,设
Figure 422520DEST_PATH_IMAGE228
是一个长度为
Figure 770325DEST_PATH_IMAGE229
的有限样本序列，设
Figure 241757DEST_PATH_IMAGE229
为偶数且
Figure 358618DEST_PATH_IMAGE230
,采样间隔
Figure 850779DEST_PATH_IMAGE231
,因此，动态相量计算的时间参考位于样本记录中的
Figure 561246DEST_PATH_IMAGE232
，离散化后的信号表达式
Figure 24037DEST_PATH_IMAGE233
为：
Figure 870771DEST_PATH_IMAGE234
(7)
其中，
Figure 64992DEST_PATH_IMAGE235
是原式中
Figure 997176DEST_PATH_IMAGE236
的系数，表示大小为
Figure 935045DEST_PATH_IMAGE237
的矩阵，该矩阵为指数项的泰勒傅里叶基矩阵,为防止与式(2)中同类字母相量混淆，另设
Figure 636284DEST_PATH_IMAGE238
，作为长度为
Figure 1407DEST_PATH_IMAGE239
的列向量，描述了当谐波
Figure 420887DEST_PATH_IMAGE240
时
Figure 132753DEST_PATH_IMAGE241
的集合，∗表示共轭计算符,即
Figure 954078DEST_PATH_IMAGE242
。
1.2、多频动态信号中DDC的近似表征
电力系统信号中往往含有DDC分量，而在常用的谐波相量检测方法中，DDC分量往往因为其含量小，检测困难而被忽略,然而，一旦DDC分量发生偏置，它会严重干扰原始信号的低阶谐波分量,此时，DDC分量影响较大，不能被忽略。
一般来说，在狭窄的时间窗口中，可以将较窄时间窗内
Figure 490102DEST_PATH_IMAGE243
的衰减直流分量近似表示成一个动态、频率较低的余弦分量[23]，即:
Figure 396878DEST_PATH_IMAGE244
（8）
其中
Figure 410970DEST_PATH_IMAGE245
和
Figure 352382DEST_PATH_IMAGE246
分别为模型DDC分量的振幅和初始相位，
Figure 793727DEST_PATH_IMAGE247
为算法观测时间窗长度，
Figure 810968DEST_PATH_IMAGE248
为一个较低的频率，当DDC时间常数很大时，DDC分量倾向于直流分量,在此条件下，
Figure 504118DEST_PATH_IMAGE249
Hz。
接下来，基于频域采样定理，将DDC频率余弦分量对应的动态相量记为
Figure 424669DEST_PATH_IMAGE250
,对幅值有限的时域信号
Figure 912283DEST_PATH_IMAGE251
，可基于频域采样定理对其进行参数化建模[24]，具体表示为:
Figure 918285DEST_PATH_IMAGE252
（9）
其中，
Figure 415125DEST_PATH_IMAGE253
为相量
Figure 190183DEST_PATH_IMAGE254
在频率
Figure 114277DEST_PATH_IMAGE255
处的频率采样值，
Figure 843461DEST_PATH_IMAGE256
表示一个向下取整算子，
Figure 878413DEST_PATH_IMAGE257
是频域采样间隔;
Figure 507977DEST_PATH_IMAGE258
代表用来进行
Figure 602972DEST_PATH_IMAGE259
参数化建模的频域采样点数,为提高上式模型的准确度，一般要求
Figure 317987DEST_PATH_IMAGE260
，对
Figure 156630DEST_PATH_IMAGE261
进行建模后，可以实现动态DDC余弦分量和动态基波分量的近似表征。
本发明中，将
Figure 906280DEST_PATH_IMAGE262
设为一个奇数，使得
Figure 172177DEST_PATH_IMAGE263
时刻在观测窗口的中心,式(8)中拟合的DDC分量的离散形式可用(10)表示。
Figure 878882DEST_PATH_IMAGE264
(10)
其中，
Figure 255637DEST_PATH_IMAGE265
为(8)右侧表达式中余弦信号的
Figure 594215DEST_PATH_IMAGE266
个样本点的列向量,式中
Figure 155646DEST_PATH_IMAGE265
为一个包含了信号
Figure 986199DEST_PATH_IMAGE267
的
Figure 25699DEST_PATH_IMAGE266
个采样点的列向量，
Figure 359728DEST_PATH_IMAGE268
是一个包含了
Figure 593526DEST_PATH_IMAGE269
的行向量，
Figure 645795DEST_PATH_IMAGE270
为一个每列包含了
Figure 754566DEST_PATH_IMAGE266
个
Figure 943102DEST_PATH_IMAGE271
采样点的矩阵,
Figure 580756DEST_PATH_IMAGE272
为共轭算子。
表达式(10)的最小二乘技术提供了最好的参数，因为它获得了
Figure 385901DEST_PATH_IMAGE268
和二阶泰勒逼近之间的最小误差,此时是满足下列约束条件下的最优解：
Figure 767204DEST_PATH_IMAGE273
(11)
其中，
Figure 75826DEST_PATH_IMAGE274
表示欧几里得范数，再引入拉格朗日算子和埃尔米特运算符进行求导，相量的系数向量可以计算为：
Figure 648496DEST_PATH_IMAGE275
(12)
式中H表示埃尔米特算子,如此，根据式(9)和式(12)，就可以得到DDC分量的估计值为:
Figure 675358DEST_PATH_IMAGE276
(13)
2、动态相量的估计
Figure 860352DEST_PATH_IMAGE277
为原始信号，为了区别于参考信号，
Figure 23480DEST_PATH_IMAGE278
也被称为目标信号,我们可以通过计算其
Figure 534095DEST_PATH_IMAGE279
个线性测量值来获得观测数据,为使公式更加美观简洁，这里用代替(7)中
Figure 782674DEST_PATH_IMAGE280
，用以下公式表示:
Figure 505780DEST_PATH_IMAGE281
(14)
而
Figure 415092DEST_PATH_IMAGE282
为
Figure 706396DEST_PATH_IMAGE283
矩阵，
Figure 566905DEST_PATH_IMAGE284
为观测值,矩阵
Figure 969067DEST_PATH_IMAGE282
从
Figure 231421DEST_PATH_IMAGE285
映射到
Figure 959206DEST_PATH_IMAGE286
，它表示因为
Figure 41431DEST_PATH_IMAGE287
而产生的降维,假设信号
Figure 247285DEST_PATH_IMAGE288
是可通过正交变换
Figure 856821DEST_PATH_IMAGE289
(例如,离散小波变换)进行压缩且
Figure 755507DEST_PATH_IMAGE290
,然后可以将式(1)重写为：
Figure 590608DEST_PATH_IMAGE291
(15)
Figure 334573DEST_PATH_IMAGE292
是原始信号
Figure 571519DEST_PATH_IMAGE293
的变换系数，可采用高斯随机矩阵
Figure 375527DEST_PATH_IMAGE282
和离散小波逆变换
Figure 697924DEST_PATH_IMAGE294
，且
Figure 245580DEST_PATH_IMAGE295
的组合矩阵也是高斯的并且概率满足RIP条件,因此，在
Figure 838498DEST_PATH_IMAGE296
足够高的情况下，可以很好地重建系数
Figure 78986DEST_PATH_IMAGE292
，对原信号进行逆变换
Figure 623100DEST_PATH_IMAGE297
即可恢复原信号
Figure 974447DEST_PATH_IMAGE293
,而一旦确定原信号
Figure 920406DEST_PATH_IMAGE293
，就可以通过公式计算谐波的频率
Figure 331796DEST_PATH_IMAGE298
、振幅
Figure 363206DEST_PATH_IMAGE299
以及频率变化率
Figure 252664DEST_PATH_IMAGE300
。
2.1、动态相量估计模型
电能质量监测数据往往会被噪声和尖峰异常值污染,异常值具有稀疏性，则利用L1范数进行描述,而为了获得具有合理的精度和对谐波与简谐波的鲁棒性的估计，经常依赖于原始信号的先验信息,然后将
Figure 551665DEST_PATH_IMAGE301
的估计表示为一个基于正则化的优化问题，该问题包含了这一先验信息，即为将测量结果还原为原始信号，则重建信号的目标函数改为：
Figure 258590DEST_PATH_IMAGE303
(16)
其中，
Figure 652662DEST_PATH_IMAGE304
是任意向量
Figure 204866DEST_PATH_IMAGE305
的
Figure 735205DEST_PATH_IMAGE306
范数，其
Figure 613031DEST_PATH_IMAGE307
范数计算定义为
Figure 494399DEST_PATH_IMAGE308
,
目前该领域有大量研究有助于(3)中目标函数的优化，其中大部分利用了变换域系数
Figure 617338DEST_PATH_IMAGE292
的稀疏性,由于谐波间谐波信号的特殊特性，(16)的正则化可以引入辅助信号
Figure 267763DEST_PATH_IMAGE309
，得到一个新的目标函数
Figure 50911DEST_PATH_IMAGE310
：
Figure 278630DEST_PATH_IMAGE311
(17)
Figure 313582DEST_PATH_IMAGE312
(18)
Figure 208726DEST_PATH_IMAGE313
(19)
其中，
Figure 38141DEST_PATH_IMAGE314
是梯度算子，
Figure 523130DEST_PATH_IMAGE315
是是垂直于
Figure 96194DEST_PATH_IMAGE316
的梯度的单位向量,
Figure 845844DEST_PATH_IMAGE317
是信号
Figure 111740DEST_PATH_IMAGE293
，所以
Figure 48472DEST_PATH_IMAGE318
是
Figure 690806DEST_PATH_IMAGE293
的梯度,
Figure 294963DEST_PATH_IMAGE319
,
Figure 731761DEST_PATH_IMAGE320
是垂直方向和水平方向的偏导数。
新的正则化是参考信号特征的方向向量
Figure 922833DEST_PATH_IMAGE315
和目标信号的梯度向量
Figure 103278DEST_PATH_IMAGE321
的点积,如果
Figure 561942DEST_PATH_IMAGE316
和
Figure 169640DEST_PATH_IMAGE293
边缘方向夹脚太小，所以优化正则化表达式
Figure 346544DEST_PATH_IMAGE322
得到：
Figure 65101DEST_PATH_IMAGE323
(20)
其中，
Figure 643850DEST_PATH_IMAGE324
是
Figure 422450DEST_PATH_IMAGE325
方向上的差算子,
Figure 319606DEST_PATH_IMAGE326
是
Figure 107433DEST_PATH_IMAGE327
方向上的差算子。
新的目标函数(18)有两个正则项:一个是系数
Figure 603005DEST_PATH_IMAGE328
的
Figure 411561DEST_PATH_IMAGE329
正则化，另一个是受参考信号
Figure 438423DEST_PATH_IMAGE316
约束的二次正则化,
Figure 124882DEST_PATH_IMAGE330
包含目标信号
Figure 288010DEST_PATH_IMAGE331
和参考信号
Figure 673992DEST_PATH_IMAGE316
的边缘方向,
Figure 47204DEST_PATH_IMAGE315
指向新正则化的方向，
Figure 35889DEST_PATH_IMAGE332
控制新正则化的强度,如果
Figure 53523DEST_PATH_IMAGE316
的梯度较小，而
Figure 735041DEST_PATH_IMAGE315
方向上的
Figure 470915DEST_PATH_IMAGE333
投影较大，则需要执行正则化惩罚,这也可以理解为，向量
Figure 513825DEST_PATH_IMAGE315
控制平滑的方向，
Figure 385966DEST_PATH_IMAGE334
控制平滑的强度,混合正则化可以提高压缩感知重构的性能。
2.2、重构动态相量方程
分裂布雷格曼方法的基本思想是引入一个辅助变量，替换原目标泛函中较难处理部分，从而简化问题的求解,本发明在全变分的基础上,通过多次迭代求解布雷格曼模型来减小信号结构损失的问题,新的目标函数中，基于分裂的准则，引入辅助变量
Figure 238384DEST_PATH_IMAGE335
,
Figure 461555DEST_PATH_IMAGE336
；对于这两个约束，
Figure 792042DEST_PATH_IMAGE337
被新的正则化参数
Figure 784269DEST_PATH_IMAGE338
和
Figure 807589DEST_PATH_IMAGE339
取代，用来控制二次惩罚函数项,并利用二次惩罚项转化为一个无约束的最优化分裂模型：
Figure 612996DEST_PATH_IMAGE340
(21)
基于分裂布雷格曼迭代，可以有效地计算
Figure 622540DEST_PATH_IMAGE328
的目标函数,迭代过程被分解为几个步骤，信号重建可以通过求解下述优化问题来实现:
Figure 859487DEST_PATH_IMAGE342
(22)
Figure 663495DEST_PATH_IMAGE344
(23)
Figure 985892DEST_PATH_IMAGE346
(24)
与传统的求解等式约束优化问题的惩罚函数与延拓方法比，布雷格曼迭代具有收敛速度快与数值稳定性好的优点,非常适用于求解
Figure 533548DEST_PATH_IMAGE347
范数正则化问题，也符合信号重建的求解问题。
其中有两个重要的性质：首先,在公式(22)的迭代过程中,
Figure 625000DEST_PATH_IMAGE348
单调递减,直至为0.其次,在(22)式的迭代求解中,只要满足
Figure 599910DEST_PATH_IMAGE349
，
Figure 908138DEST_PATH_IMAGE328
单调趋近真实解
Figure 259485DEST_PATH_IMAGE350
因此,布雷格曼方法具有稳定的收敛性。
在从(22)到(24)的每个迭代周期中，(22)是一个可微优化问题,式(24)显式求解并直接更新,利用定义
Figure 939865DEST_PATH_IMAGE351
函数来对变量
Figure 351255DEST_PATH_IMAGE352
进行更新，可以有效地求解式(14)，即：
Figure 382665DEST_PATH_IMAGE354
(25)
其中，
Figure 272123DEST_PATH_IMAGE351
函数定义为：
Figure 338168DEST_PATH_IMAGE356
(26)
在公式(22)中，
Figure 654880DEST_PATH_IMAGE357
，
Figure 675051DEST_PATH_IMAGE358
,
Figure 368200DEST_PATH_IMAGE359
,
Figure 23173DEST_PATH_IMAGE360
和
Figure 900999DEST_PATH_IMAGE361
都可以被视为常数矩阵和向量,只有
Figure 782367DEST_PATH_IMAGE328
是可变的,为了最小化(22)，在第
Figure 138262DEST_PATH_IMAGE362
次迭代中，将对
Figure 788686DEST_PATH_IMAGE363
的一阶导数设为0.然而，在此之前，我们首先求出
Figure 341808DEST_PATH_IMAGE364
相对于
Figure 710473DEST_PATH_IMAGE328
的一阶导数，然后我们得到：
Figure 870059DEST_PATH_IMAGE366
(27)
基于(27)，我们定义线性正则化算子
Figure 374989DEST_PATH_IMAGE367
为：
Figure 594618DEST_PATH_IMAGE369
(28)
根据(28)中
Figure 450579DEST_PATH_IMAGE367
的定义和函数(18)的一阶导数得到：
Figure 882697DEST_PATH_IMAGE371
(29)
通过将(20)用
Figure 507713DEST_PATH_IMAGE363
表示，得到了
Figure 665287DEST_PATH_IMAGE363
的一个封闭解：
Figure 477386DEST_PATH_IMAGE373
(30)
其中，
Figure 244353DEST_PATH_IMAGE374
是一个单位矩阵,由于(21)是线性的，可以有效地求解
Figure 723876DEST_PATH_IMAGE363
,在本发明算法的迭代中，将(23)、(24)、(30)循环合并更新,当正则化参数时，该算法只能计算系数
Figure 285308DEST_PATH_IMAGE328
，成为基本优化。
解得系数
Figure 850281DEST_PATH_IMAGE328
后通过
Figure 889781DEST_PATH_IMAGE375
可得到重构相量
Figure 489390DEST_PATH_IMAGE376
（即
Figure 720258DEST_PATH_IMAGE377
），其每一行对应于每个频率分量在不同时间的相量,为得到更精细的频域结果，引入了插值因子F,频率分辨率为
Figure 772527DEST_PATH_IMAGE378
，重构相量频率可近似表示为
Figure 350139DEST_PATH_IMAGE379
，
Figure 928888DEST_PATH_IMAGE380
为频率索引,而一旦确定
Figure 973067DEST_PATH_IMAGE376
，所求的谐波的频率
Figure 371688DEST_PATH_IMAGE381
、振幅
Figure 893936DEST_PATH_IMAGE382
以及频率变化率
Figure 828656DEST_PATH_IMAGE383
可通过式(31-33)进行计算。
Figure 778158DEST_PATH_IMAGE385
(31)
Figure 664074DEST_PATH_IMAGE387
(32)
Figure 990013DEST_PATH_IMAGE389
(33)
其中，
Figure 277775DEST_PATH_IMAGE390
为
Figure 398178DEST_PATH_IMAGE391
的第
Figure 36969DEST_PATH_IMAGE392
列，
Figure 635441DEST_PATH_IMAGE393
，
Figure 270385DEST_PATH_IMAGE394
和
Figure DEST_PATH_IMAGE395
分别为
Figure DEST_PATH_IMAGE396
的零阶导数，一阶导数和二阶导数。
综上所述，本发明将动态谐波频域分布的规律性作为动态相量恢复的优化目标，应用基于布雷格曼分割迭代正则化模型算法来重构含有各次谐波泰勒展开系数的动态相量，算法流程如图5所示。
3、动态相量的估计
根据IEEE测量标准，动态相量测量误差主要包括相量总误差(total vectorerror，TVE)、频率误差(frequency error，FE)和频率变化率误差(ROCOF error，RFE)。本节中，将在不同的测试条件下比较本发明提出的估计算法与文献中对比算法的精度，并进行讨论。其中对比算法包括：基于分裂布雷格曼迭代的压缩感知方法、泰勒傅里叶变换(TFT)、压缩感知泰勒傅里叶多频分析方法(CSTFM)、IpDFT和O样条设计的FIR滤波器(OFF)算法。测试场景包括基本性能测试、频率偏差测试、谐波震荡测试以及测试。
本发明在不同频率偏差和谐波失真、高斯白噪声和光噪声干扰的稳态条件下，对各种谐波信号进行仿真试验。此外，设计了频率渐变和动态调制试验来验证其动态性能。5种对比算法都使用相同的矩形观测窗口，对比算法的相量导数和谐波RFE估计滤波器也基于本发明算法相同方案设计。
3.1、基本性能测试
为了验证本发明算法在信号频率偏移时的有效性，假设采样窗长为5个周期，基频带宽为1Hz，构造式(34)所示含基波及动态分量的信号模型：
Figure DEST_PATH_IMAGE397
(34)
式中，f1为基波频率，此处设为50Hz，
Figure DEST_PATH_IMAGE398
、
Figure DEST_PATH_IMAGE399
分别代表基波和各次谐波相角，取
Figure DEST_PATH_IMAGE400
范围中的任意值。低频段谐波次数h取值为2-13，采样频率设置为5 kHz。
在运用本发明算法时，
Figure DEST_PATH_IMAGE401
,
Figure DEST_PATH_IMAGE402
均初始化为全零矩阵，正则化参数。又有
Figure DEST_PATH_IMAGE403
。
Figure DEST_PATH_IMAGE404
为参数间平衡系数，在
Figure DEST_PATH_IMAGE405
内算法稳定性较高，k为迭代次数。本发明根据相似数据点大多分布在样本点的邻近位置，将k-最近邻聚类算法搜索窗口设置以当前样本点为中心。在参数
Figure DEST_PATH_IMAGE406
、k变化情况下，对重构效果及算法运行时间进行分析，结果如图1所示。
由图1可知，随着迭代次数k的增大，总相量误差逐渐减小，但是算法运行时间却不断增加，当k达到800后，改变迭代次数对精度影响较小且此时重构趋于稳定。当平衡参数
Figure 561689DEST_PATH_IMAGE406
在0.2到0.3之间变化时，计算精度先增后减，在0.25附件达到峰值。因此，适当减少迭代次数数量k，选择适当的平衡参数能够有效提高重构效果。综合重构性能和测量精度两方面考虑，本发明应用算法中参数k=800，平衡参数
Figure 297564DEST_PATH_IMAGE406
=0.25。
选择OFF、CSTFM、TFT以及IpDFT作为对比算法，本发明与对比算法的总相量误差估计、频率和频率变化率误差估计结果如下表所示，各算法估计精度比较：
Figure DEST_PATH_IMAGE407
本发明与对比算法的总相量（TVE）误差估计结果表
Figure DEST_PATH_IMAGE408
本发明与对比算法的总相量频率（FE）误差估计结果表
Figure DEST_PATH_IMAGE409
本发明与对比算法的总相量频率变化率（RFE）误差估计结果表
由上表可知，SBCS算法TVE、FE、RFE指标的最大值分别为0.591%、0.057Hz、2.181Hz/s。根据IEEE标准TVE，FE和RFE值为1.5%，0.06Hz和2.3Hz/s，可以看出本发明方法完全可以满足IEEE的测量标准。且本发明算法SBCS的TVE，RFE，FE值都小于其他算法，即本发明算法SBCS方法对于含ddc分量的动态信号具有较好的检测能力。
OFF算法误差估计接近本发明算法，但其TVE、FE、RFE指标仍然不满足IEEE测量标准。原因是该方法由于空间步长重构过程中存在较大的噪声，导致频率误差较大。至于TFT，CSTFM和IpDFT方法，它们都是通过二阶泰勒模型来重构ddc分量，其泰勒模型固有的展开阶数会在重构过程中存在一定的误差，故精度不高。而CSTFM方法在高次谐波附近的精度仅次于本发明的SBCS方法，但在低次谐波附近的精度接近TFT方法，TVE、FE、RFE指标的最大值分别为7.418%、1.736Hz、25.211Hz/s，难以满足IEEE的测量标准。IpDFT和TFT方法最高总相量误差为10.456%和8.872% FE与RFE测量结果精度也不够理想。动态条件下，傅里叶变换模型不能跟踪观察窗口中发生的相量变化，从而导致不正确的相量评估。
3.2、频率偏差测试
当电力系统中发生失步时，电压信号的频率将发生连续变化，此时PMU的相量、频率、频率变化率的量测精度对失步解列控制尤为重要。IEEE规定电力系统频率绝对偏差应始终小于0.5Hz。采用
Figure DEST_PATH_IMAGE410
=1.15Hz，以获得每个谐波频率附近良好的通带和阻带性能。如前所述，本发明特别考虑13阶以下的谐波。采样频率设置为10kHz。每次测试的结果是经过1000次重复运行，其中基相和谐相均为随机数分布。因此，为评估测试算法频率偏差条件下带来的影响，本节仍以上述的五种算法作为对比算法，具体的多频动态信号如下所示:
Figure DEST_PATH_IMAGE412
(35)
图2是第2-13次谐波的最大TVEs、FEs和Fes图。式中，
Figure DEST_PATH_IMAGE413
的变化范围为49.5-50.5Hz，阶跃为0.2 Hz。在这种情况下，我们可以看到SBCS算法在谐波相量、频率和ROCOF估计方面总是比其他四种算法更准确。这是因为基于动态相量高阶导数的模型有助于获得更好的通带和阻带性能，特别是对于高阶谐波。在此条件下，SBCS的最大TVE、FE和RFE分别为0.30%、0.025 Hz和0.2 Hz/s。在IEEE标准中，它们的阈值分别为1.5%、0.01 Hz和0.4 Hz/s。因此，SBCS都满足估计要求。当信号波形受到频率偏移影响时，所提出的估计算法能表现出更高的性能。
图2中可以看出，在2-13次谐波估计方面，其他方法均不满足IEEE测量标准。IpDFT方法的TVE、FE和RFE相比其他方法误差较大，可认为其估计精度最低，对频率偏移调制的敏感性相当高。OFF，TFT和CSTFM方法的估计精度显然要好于IpDFT方法。其中，CSTFM方法是依据动态模型来估计相量测量，误差值受频率变化影响较IpDFT更小。由于TFT方法是根据泰勒二阶模型来估计相量，所以测量精确度不高，且受频率偏移影响较大。OFF算法利用O样条作为采样算子，得到最佳的泰勒-傅立叶系数。它能够实现在谐波频率上的调制，测量精度接近IEEE标准。
3.3、谐波震荡测试
本节测试中，具体的多频动态信号如下：
Figure DEST_PATH_IMAGE414
(36)
本节测试中采样率设定为5kHz，采样周期长度为5个周期。其中
Figure DEST_PATH_IMAGE415
为调制频率，设置为5Hz。其他参数与上章节的值相同。
本发明算法估计结果图中仅展示2-13次谐波的参数估计结果。估计结果如图3所示。谐波振荡对SBCS算法在低阶谐波(如第2-7次)的影响结果较大。然而，本发明算法在高次谐波参数估计方面比其他方法更精确，特别是在谐波频率和ROCOF估计方面。在此条件下，SBCS的TVE、FE和RFE整体最大值分别为2.27%、0.05Hz和2.14Hz/s。CSTFM方法为6.42%，1.89Hz和7.82 Hz/s，TFT方法为8.9%，2.2Hz和18.74 Hz/s，OFF方法为3.84%，0.25Hz和4.25Hz/s，IpDFT方法为25.47%，3.44Hz和21.29Hz/s。以上数据表明，本发明算法的TVE，FE和RFE值最小，可见本发明提出的方法在低频振荡的场景下有更高、更稳定的估计精度。此外，在4.3节测试条件下，IEEE标准中对应的阈值分别为3.5%、0.08 Hz和2.5 Hz/s，其他四个对比算法都不满足要求。TFT方法在本节谐波震荡情况下，相邻谐波间会发生严重干扰，进而影响到TFT方法的估计精度。OFF提供了一种最优的振荡数据压缩算法，误差由样条的阶数和埃尔米特样条的阶数控制。所以误差影响较小，接近测量标准。CSTFM方法引入了TF多频（TFM）模型来描述动态相量，但模型频率无法准确选择。这种不正确的信号模型将导致较大的误差。IpDFT方法在本节谐波振荡情况下，产生严重频谱泄漏和互谐波干扰，不能很好地适应谐波频率，误差也较大。
3.4、抗干扰能力测试
通常的情况下，电力系统信号中都会含有一定的间谐波和噪声，这严重影响到谐波相量的估计。本节测试中，信号引入信噪比为60dB的高斯白噪声。具体的动态信号如下：
Figure DEST_PATH_IMAGE416
(37)
其中
Figure DEST_PATH_IMAGE417
为间谐波频率。本节测试中采样率设定为5kHz，采样周期长度为5个周期。具体的仿真结果如图4[a]、图4[b]、图4[c]所示。图4[a]、图4[b]、图4[c]是间谐波和噪声干扰条件下各方法的估计结果图，图4[a]Max,TVE，图4[b]Max,FE，图4[c]Max,RFE。
从图4[a]、图4[b]、图4[c]中可以看出本发明算法SBCS方法的TVE、RFE、FE值都小于其他算法。这可认为本发明算法在动态条件下估计精度最高。尽管在此过程中出现了局部的噪声，但重建效果几乎是完美。本发明方法的最大TVE，FE和RFE值分别为3.52%，0.44Hz和2.35Hz/s。OFF方法的最大TVE，FE和RFE值分别为4.84%，1.26Hz和5.1Hz/s。TFT方法的最大TVE，FE和RFE值分别为9.97%，2.65Hz和14.93Hz/s。IpDFT方法的最大TVE，FE和RFE值分别为11.36%，2.65Hz和14.98Hz/s。CSTFM方法的最大TVE，FE和RFE值分别为8.42%，1.32Hz和6.88Hz/s。根据本节测试条件下的IEEE标准，其最大TVE，FE和RFE值的阈值为3.4%，0.45Hz，和2.5Hz/s。本发明算法完全满足标准要求，OFF算法具有带限信号的最佳空间采样器，为振荡提供了一个最佳的数据压缩算法，以样条的程度递增，提供了一个强大的最优状态估计器，有效抑制了噪声和间谐波干扰。TFT方法在本节噪声干扰情况下，相邻谐波间会发生严重干扰，进而使TFT方法的估计误差变大。IpDFT算法不能有效抑制位于基波附近的间谐波分量对基波频率估计带来的影响，并且受基频偏移的影响。本发明算法SBCS在多频相位分析能力中减少了间谐波分量的时变性及其噪声对动态相量测量的影响，显著提高了测量精度。
4、结论
本发明提出了基于CS-TFM的动态谐波相量估计算法，用于估计动态谐波相量、频率和频率变化率。通过分析计算和仿真测试，验证了本发明算法的优越性，并得到了一些结论。
在各类仿真和实例测试中，即便将系统处于动态条件、间谐波和噪声干扰状态下，本发明算法的性能是值得接受的，并且对于动态谐波相量、频率和频率变化率大部分时间内满足对于谐波相量估计的IEEE标准。
本发明旨在精确测量全TFM波形模型中的谐波和相间分量，首先，基于泰勒-傅里叶多频变换建立可以同时估计谐波分量和间谐波分量的动态信号模型，其次，基于动态信号分布的规律性和稀疏性，引入辅助信号将相量问题转化为压缩感知问题，最后，将动态相量重构问题转化为分裂布雷格曼迭代对混合正则化算法的优化问题，对信号进行重构，得到动态相量估计，该方法中还考虑了高阶导数，允许同时估计相量，这使得相量估计更加准确和高效，因此，实现了足够的测量精度，在多频相位分析能力中减少了谐波间谐波分量的时变性及其对动态相量测量的影响，显著提高了测量精度。
以上所述实施例仅表达了本发明的具体实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。

Claims (9)

1.一种多频动态相量测量方法，其特征在于，包括以下步骤：
a、DDC分量的动态相量估计步骤：电力系统的多频动态信号模型
Figure 30390DEST_PATH_IMAGE001
由振幅和相位时变的正弦分量和表示：
Figure 535321DEST_PATH_IMAGE002
（1）
其中，其中
Figure 364737DEST_PATH_IMAGE003
为信号中的DDC分量，
Figure 689539DEST_PATH_IMAGE004
为基波及谐波信号，
Figure 26717DEST_PATH_IMAGE005
和
Figure 386154DEST_PATH_IMAGE006
分别为DDC分量的幅值和时间常数，
Figure 386471DEST_PATH_IMAGE007
和
Figure 198569DEST_PATH_IMAGE008
分别为谐波相量的幅值和时间常数，
Figure 575324DEST_PATH_IMAGE009
是包含实际电力系统频率
Figure 54847DEST_PATH_IMAGE010
，其谐波倍数和可能的间谐波频率的通用频率集合，然后建立多频动态信号的TFM模型和得出多频动态信号中DDC的近似表征；
b、动态相量的估计步骤：设定
Figure 730460DEST_PATH_IMAGE011
为原始信号，
Figure 295433DEST_PATH_IMAGE012
为目标信号，计算其
Figure 210300DEST_PATH_IMAGE013
个线性测量值来获得观测数据，用以下公式表示:
Figure 544329DEST_PATH_IMAGE014
（14）
其中，
Figure 152028DEST_PATH_IMAGE015
为
Figure 938718DEST_PATH_IMAGE016
矩阵，
Figure 922855DEST_PATH_IMAGE017
为观测值，矩阵
Figure 609926DEST_PATH_IMAGE018
从
Figure 122947DEST_PATH_IMAGE019
映射到
Figure 396933DEST_PATH_IMAGE020
，它表示因为
Figure 919182DEST_PATH_IMAGE021
而产生的降维，假设信号
Figure 227803DEST_PATH_IMAGE022
是可通过正交变换
Figure 911725DEST_PATH_IMAGE023
进行压缩且
Figure 437122DEST_PATH_IMAGE024
，然后将公式(1)重写为：
Figure 497482DEST_PATH_IMAGE025
(15)
其中，
Figure 395031DEST_PATH_IMAGE026
是原始信号
Figure 781013DEST_PATH_IMAGE027
的变换系数，采用高斯随机矩阵
Figure 29592DEST_PATH_IMAGE028
和离散小波逆变换
Figure 893643DEST_PATH_IMAGE029
，且
Figure 911277DEST_PATH_IMAGE030
的组合矩阵也是高斯的并且概率满足RIP条件，在
Figure 701116DEST_PATH_IMAGE031
足够高的情况下，重建系数
Figure 436991DEST_PATH_IMAGE032
，对原信号
Figure 104733DEST_PATH_IMAGE033
进行逆变换
Figure 976874DEST_PATH_IMAGE034
即可恢复原信号
Figure 173500DEST_PATH_IMAGE035
，如果确定原信号
Figure 131092DEST_PATH_IMAGE036
，就通过公式计算谐波的频率
Figure 336945DEST_PATH_IMAGE037
、振幅
Figure 821848DEST_PATH_IMAGE038
以及频率变化率
Figure 720534DEST_PATH_IMAGE039
；然后建立动态相量估计模型和重构动态相量方程；
c、动态相量的设计步骤：建立基本性能测试场景、频率偏差测试场景、谐波震荡测试场景和抗干扰能力测试场景。
2.根据权利要求1所述的一种多频动态相量测量方法，其特征在于，步骤a中，“建立多频动态信号的TFM模型”的步骤如下：
将公式(1)中
Figure 165422DEST_PATH_IMAGE040
的每个分量与一个动态相量相关联，定义为：
Figure 909387DEST_PATH_IMAGE041
Figure 21699DEST_PATH_IMAGE042
(2)
将公式（1）中的
Figure 560128DEST_PATH_IMAGE043
方程改写为：
Figure 757891DEST_PATH_IMAGE044
(3)
基于电力系统的大惯性特征，将用泰勒级数近似表示为：
Figure 538503DEST_PATH_IMAGE045
(4)
其中，
Figure 239743DEST_PATH_IMAGE046
为
Figure 949073DEST_PATH_IMAGE047
的
Figure 368553DEST_PATH_IMAGE048
阶导数，
Figure 454320DEST_PATH_IMAGE049
为泰勒展开阶数，K取3，即：
Figure 508602DEST_PATH_IMAGE050
(5)
所以，
Figure 919992DEST_PATH_IMAGE051
的离散表达式为:
Figure 826768DEST_PATH_IMAGE052
(6)
其中，T为采样间隔，∗表示共轭计算符；
将公式(2)代入公式(6)得信号的泰勒展开表达式，再以采样频率
Figure 450647DEST_PATH_IMAGE053
进行采样.设
Figure 392058DEST_PATH_IMAGE054
是一个长度为
Figure 974349DEST_PATH_IMAGE055
的有限样本序列，设
Figure 368422DEST_PATH_IMAGE056
为偶数且
Figure 294527DEST_PATH_IMAGE057
，采样间隔
Figure 824866DEST_PATH_IMAGE058
，因此，动态相量计算的时间参考位于样本记录中的
Figure 578058DEST_PATH_IMAGE059
，离散化后的信号表达式
Figure 459426DEST_PATH_IMAGE060
为：
Figure 956267DEST_PATH_IMAGE061
Figure 341112DEST_PATH_IMAGE062
(7)
其中，
Figure 999626DEST_PATH_IMAGE063
是原式中
Figure 872685DEST_PATH_IMAGE064
的系数，表示大小为
Figure 642058DEST_PATH_IMAGE065
的矩阵，该矩阵为指数项的泰勒傅里叶基矩阵，为防止与式(2)中同类字母相量混淆，另设
Figure 412568DEST_PATH_IMAGE066
，作为长度为
Figure 241984DEST_PATH_IMAGE067
的列向量，当谐波
Figure 97944DEST_PATH_IMAGE068
时
Figure 671008DEST_PATH_IMAGE069
的集合，∗表示共轭计算符，即
Figure 30445DEST_PATH_IMAGE070
。
3.根据权利要求2所述的一种多频动态相量测量方法，其特征在于，步X骤a中，“得出多频动态信号中DDC的近似表征”的步骤如下：
将较窄时间窗内
Figure 794876DEST_PATH_IMAGE071
的衰减直流分量近似表示成一个动态、频率较低的余弦分量，即:
Figure 872554DEST_PATH_IMAGE072
（8）
其中，
Figure 249309DEST_PATH_IMAGE073
和
Figure 728831DEST_PATH_IMAGE074
分别为模型DDC分量的振幅和初始相位，
Figure 165629DEST_PATH_IMAGE075
为算法观测时间窗长度，
Figure 730603DEST_PATH_IMAGE076
为一个较低的频率，当DDC时间常数很大时，DDC分量倾向于直流分量，在此条件下，
Figure 911048DEST_PATH_IMAGE077
Hz；
接下来，基于频域采样定理，将DDC频率余弦分量对应的动态相量记为
Figure 245078DEST_PATH_IMAGE078
，对幅值有限的时域信号
Figure 85732DEST_PATH_IMAGE079
，基于频域采样定理对其进行参数化建模，具体表示为:
Figure 138002DEST_PATH_IMAGE080
（9）
其中，
Figure 856559DEST_PATH_IMAGE081
为相量
Figure 310674DEST_PATH_IMAGE082
在频率
Figure 89275DEST_PATH_IMAGE083
处的频率采样值，
Figure 363261DEST_PATH_IMAGE084
表示一个向下取整算子，
Figure 885509DEST_PATH_IMAGE085
是频域采样间隔，
Figure 692666DEST_PATH_IMAGE086
代表用来进行
Figure 642167DEST_PATH_IMAGE087
参数化建模的频域采样点数，为提高上式模型的准确度，要求
Figure 137871DEST_PATH_IMAGE088
，对
Figure 463810DEST_PATH_IMAGE089
进行建模后，实现动态DDC余弦分量和动态基波分量的近似表征；
将
Figure 626938DEST_PATH_IMAGE090
设为一个奇数，使得
Figure 747341DEST_PATH_IMAGE091
时刻在观测窗口的中心，公式(8)中拟合的DDC分量的离散形式用公式(10)表示：
Figure 995920DEST_PATH_IMAGE092
(10)
其中，
Figure 376084DEST_PATH_IMAGE093
为公式(8)右侧表达式中余弦信号的
Figure 393718DEST_PATH_IMAGE094
个样本点的列向量，式中
Figure 950601DEST_PATH_IMAGE095
为一个包含了信号
Figure 420897DEST_PATH_IMAGE096
的
Figure 88639DEST_PATH_IMAGE097
个采样点的列向量，
Figure 695201DEST_PATH_IMAGE098
是一个包含了
Figure 688564DEST_PATH_IMAGE099
的行向量，
Figure 144691DEST_PATH_IMAGE100
为一个每列包含了
Figure 350545DEST_PATH_IMAGE101
个
Figure 77192DEST_PATH_IMAGE102
采样点的矩阵,
Figure 975878DEST_PATH_IMAGE103
为共轭算子；
公式(10)的最小二乘技术提供了最好的参数，因为它获得了
Figure 420766DEST_PATH_IMAGE104
和二阶泰勒逼近之间的最小误差.此时是满足下列约束条件下的最优解：
Figure 430310DEST_PATH_IMAGE105
(11)
其中，
Figure 542623DEST_PATH_IMAGE106
表示欧几里得范数，再引入拉格朗日算子和埃尔米特运算符进行求导，相量的系数向量计算为：
Figure 845166DEST_PATH_IMAGE107
(12)
式中H表示埃尔米特算子，根据公式(9)和公式(12)，得到DDC分量的估计值为:
Figure 42929DEST_PATH_IMAGE108
(13)。
4.根据权利要求1所述的一种多频动态相量测量方法，其特征在于，步骤b中，“建立动态相量估计模型”的步骤如下：
将测量结果还原为原始信号，则重建信号的目标函数改为：
Figure 590585DEST_PATH_IMAGE109
(16)
其中，
Figure 291825DEST_PATH_IMAGE110
是任意向量
Figure 266734DEST_PATH_IMAGE111
的
Figure 951793DEST_PATH_IMAGE112
范数，其
Figure 303140DEST_PATH_IMAGE113
范数计算定义为
Figure 124466DEST_PATH_IMAGE114
；
公式(16)的正则化引入辅助信号
Figure 768811DEST_PATH_IMAGE115
，得到一个新的目标函数
Figure 675587DEST_PATH_IMAGE116
：
Figure 299467DEST_PATH_IMAGE117
(17)
Figure 240878DEST_PATH_IMAGE118
(18)
Figure 292011DEST_PATH_IMAGE119
(19)
其中，
Figure 686083DEST_PATH_IMAGE120
是梯度算子，
Figure 618048DEST_PATH_IMAGE121
是是垂直于
Figure 679545DEST_PATH_IMAGE122
的梯度的单位向量，
Figure 167158DEST_PATH_IMAGE123
是信号
Figure 48526DEST_PATH_IMAGE124
，所以
Figure 279787DEST_PATH_IMAGE125
是
Figure 930211DEST_PATH_IMAGE126
的梯度，
Figure 854305DEST_PATH_IMAGE127
,
Figure 222970DEST_PATH_IMAGE128
是垂直方向和水平方向的偏导数；
新的正则化是参考信号特征的方向向量
Figure 490878DEST_PATH_IMAGE129
和目标信号的梯度向量
Figure 995808DEST_PATH_IMAGE130
的点积，如果
Figure 90803DEST_PATH_IMAGE131
和
Figure 681185DEST_PATH_IMAGE132
边缘方向夹脚太小，所以优化正则化表达式
Figure 254248DEST_PATH_IMAGE133
得到：
Figure 879265DEST_PATH_IMAGE134
(20)
其中，
Figure 378117DEST_PATH_IMAGE135
是
Figure 455794DEST_PATH_IMAGE136
方向上的差算子,
Figure 832549DEST_PATH_IMAGE137
是
Figure 46493DEST_PATH_IMAGE138
方向上的差算子；
新的目标函数(18)有两个正则项:一个是系数
Figure 483290DEST_PATH_IMAGE139
的
Figure 48264DEST_PATH_IMAGE140
正则化，另一个是受参考信号
Figure 697551DEST_PATH_IMAGE141
约束的二次正则化，
Figure 795695DEST_PATH_IMAGE142
包含目标信号
Figure 137814DEST_PATH_IMAGE143
和参考信号
Figure 190084DEST_PATH_IMAGE144
的边缘方向，
Figure 908641DEST_PATH_IMAGE145
指向新正则化的方向，
Figure 362756DEST_PATH_IMAGE146
控制新正则化的强度，如果
Figure 875777DEST_PATH_IMAGE147
的梯度较小，而
Figure 415343DEST_PATH_IMAGE148
方向上的
Figure 164688DEST_PATH_IMAGE149
投影较大，则执行正则化惩罚，向量
Figure 473309DEST_PATH_IMAGE150
控制平滑的方向，
Figure 422811DEST_PATH_IMAGE151
控制平滑的强度.混合正则化提高压缩感知重构的性能。
5.根据权利要求4所述的一种多频动态相量测量方法，其特征在于，步骤b中，“重构动态相量方程”的步骤如下：
引入一个辅助变量，多次迭代求解布雷格曼模型来减小信号结构损失，新的目标函数中，基于分裂的准则，引入辅助变量
Figure 184094DEST_PATH_IMAGE152
、
Figure 244453DEST_PATH_IMAGE153
，对于这两个约束，
Figure 407582DEST_PATH_IMAGE154
被新的正则化参数
Figure 527984DEST_PATH_IMAGE155
和
Figure 540677DEST_PATH_IMAGE156
取代，用来控制二次惩罚函数项.并利用二次惩罚项转化为一个无约束的最优化分裂模型：
Figure 404728DEST_PATH_IMAGE157
Figure 156784DEST_PATH_IMAGE158
(21)
基于布雷格曼分割迭代，计算
Figure 713667DEST_PATH_IMAGE159
的目标函数，迭代过程被分解为几个步骤，信号重建通过求解下述优化问题来实现：
Figure 449542DEST_PATH_IMAGE160
(22)
Figure 851704DEST_PATH_IMAGE161
(23)
Figure 723845DEST_PATH_IMAGE162
(24)
在公式(22)式的迭代过程中,
Figure 950165DEST_PATH_IMAGE163
单调递减,直至为0，其次,在公式(22)的迭代求解中,只要满足
Figure 907757DEST_PATH_IMAGE164
，
Figure 113610DEST_PATH_IMAGE165
单调趋近真实解
Figure 105837DEST_PATH_IMAGE166
；
在从公式(22)到(24)的每个迭代周期中，公式(22)是一个可微优化问题，公式(24)显式求解并直接更新，利用定义
Figure 738943DEST_PATH_IMAGE167
函数来对变量
Figure 183831DEST_PATH_IMAGE168
进行更新，求解式(14)，即：
Figure 193376DEST_PATH_IMAGE169
(25)
其中，
Figure 538644DEST_PATH_IMAGE170
函数定义为：
Figure 608231DEST_PATH_IMAGE171
(26)
在公式(22)中，
Figure 274836DEST_PATH_IMAGE172
，
Figure 822492DEST_PATH_IMAGE173
,
Figure 789311DEST_PATH_IMAGE174
,
Figure 764220DEST_PATH_IMAGE175
和
Figure 449279DEST_PATH_IMAGE176
都视为常数矩阵和向量，
Figure 39442DEST_PATH_IMAGE177
是可变的，最小化(22)，在第
Figure 860767DEST_PATH_IMAGE178
次迭代中，将对
Figure 6578DEST_PATH_IMAGE179
的一阶导数设为0，求出
Figure 913354DEST_PATH_IMAGE180
相对于
Figure 68392DEST_PATH_IMAGE181
的一阶导数，然后得到：
Figure 213065DEST_PATH_IMAGE182
(27)
基于公式(27)，定义线性正则化算子
Figure 293891DEST_PATH_IMAGE183
为：
Figure 687964DEST_PATH_IMAGE184
(28)
根据公式(28)中
Figure 115534DEST_PATH_IMAGE185
的定义和函数(18)的一阶导数得到：
Figure 911452DEST_PATH_IMAGE186
(29)
通过将公式(20)用
Figure 664644DEST_PATH_IMAGE187
表示，得到了
Figure 14854DEST_PATH_IMAGE188
的一个封闭解：
Figure 511694DEST_PATH_IMAGE189
(30)
其中，
Figure 660654DEST_PATH_IMAGE190
是一个单位矩阵，公式(21)是线性的，求解
Figure 319168DEST_PATH_IMAGE191
，将公式(23)、(24)、(30)循环合并更新，当正则化参数
Figure 687832DEST_PATH_IMAGE192
时，计算系数
Figure 722785DEST_PATH_IMAGE193
，成为基本优化；
解得系数
Figure 227715DEST_PATH_IMAGE194
后通过
Figure 57131DEST_PATH_IMAGE195
可得到重构相量
Figure 647512DEST_PATH_IMAGE196
（即
Figure 719111DEST_PATH_IMAGE197
），其每一行对应于每个频率分量在不同时间的相量，引入插值因子F，频率分辨率为
Figure 344128DEST_PATH_IMAGE198
，重构相量频率表示为
Figure 610024DEST_PATH_IMAGE199
，
Figure 687701DEST_PATH_IMAGE200
为频率索引，确定
Figure 798877DEST_PATH_IMAGE201
，所求的谐波的频率
Figure 12820DEST_PATH_IMAGE202
、振幅
Figure 449618DEST_PATH_IMAGE203
以及频率变化率
Figure 530705DEST_PATH_IMAGE204
通过公式(31-33)进行计算：
Figure 711150DEST_PATH_IMAGE205
(31)
Figure 310759DEST_PATH_IMAGE206
(32)
Figure 918458DEST_PATH_IMAGE207
(33)
其中，
Figure 705148DEST_PATH_IMAGE208
为
Figure 423706DEST_PATH_IMAGE209
的第
Figure 877821DEST_PATH_IMAGE210
列，
Figure 889377DEST_PATH_IMAGE211
，
Figure 163363DEST_PATH_IMAGE212
和
Figure 685611DEST_PATH_IMAGE213
分别为
Figure 728654DEST_PATH_IMAGE214
的零阶导数，一阶导数和二阶导数。
6.根据权利要求1所述的一种多频动态相量测量方法，其特征在于，步骤c中，“建立基本性能测试场景”的步骤如下：
设定采样窗长为5个周期，基频带宽为1Hz，构造式(34)所示含基波及动态分量的信号模型：
Figure 678155DEST_PATH_IMAGE215
(34)
公式（34）中，f1为基波频率，此处设为50Hz，
Figure 705017DEST_PATH_IMAGE216
、
Figure 765377DEST_PATH_IMAGE217
分别代表基波和各次谐波相角，取
Figure 692620DEST_PATH_IMAGE218
范围中的任意值，低频段谐波次数h取值为2-13，采样频率设置为5kHz；
Figure 813022DEST_PATH_IMAGE219
,
Figure 796022DEST_PATH_IMAGE220
均初始化为全零矩阵，正则化参数
Figure 660073DEST_PATH_IMAGE221
，设定
Figure 943286DEST_PATH_IMAGE222
，
Figure 969011DEST_PATH_IMAGE223
为参数间平衡系数，在
Figure 970465DEST_PATH_IMAGE224
内算法稳定性较高，k为迭代次数，在参数
Figure 871163DEST_PATH_IMAGE225
、k变化情况下，对重构效果及算法运行时间进行分析。
7.根据权利要求6所述的一种多频动态相量测量方法，其特征在于，步骤c中，“建立频率偏差测试场景”的步骤如下：
基频带宽设置为1Hz，采用
Figure 8883DEST_PATH_IMAGE226
= 1.15Hz，获得每个谐波频率附近良好的通带和阻带性能，采样频率设置为10kHz，多频动态信号如下所示:
Figure 205509DEST_PATH_IMAGE227
(35)
公式（35）中，
Figure 163101DEST_PATH_IMAGE228
的变化范围为49.5-50.5Hz，阶跃为0.2Hz。
8.根据权利要求7所述的一种多频动态相量测量方法，其特征在于，步骤c中，“建立谐波震荡测试场景”的步骤如下：
建立多频动态信号如下：
Figure 900113DEST_PATH_IMAGE229
(36)
采样率设定为5kHz，采样周期长度为5个周期，其中
Figure 626760DEST_PATH_IMAGE230
为调制频率，设置为5Hz。
9.根据权利要求8所述的一种多频动态相量测量方法，其特征在于，步骤c中，“建立抗干扰能力测试场景”的步骤如下：
信号引入信噪比为60dB的高斯白噪声，具体的动态信号如下：
Figure 259867DEST_PATH_IMAGE231
(37)
其中
Figure 474729DEST_PATH_IMAGE232
为间谐波频率，采样率设定为5kHz，采样周期长度为5个周期。
CN202210953016.9A 2022-08-08 2022-08-08 一种多频动态相量测量方法 Pending CN115032453A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210953016.9A CN115032453A (zh) 2022-08-08 2022-08-08 一种多频动态相量测量方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210953016.9A CN115032453A (zh) 2022-08-08 2022-08-08 一种多频动态相量测量方法

Publications (1)

Publication Number Publication Date
CN115032453A true CN115032453A (zh) 2022-09-09

Family

ID=83129942

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210953016.9A Pending CN115032453A (zh) 2022-08-08 2022-08-08 一种多频动态相量测量方法

Country Status (1)

Country Link
CN (1) CN115032453A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117948872A (zh) * 2024-03-27 2024-04-30 常州同惠电子股份有限公司 一种基于lm法的电感式位置传感器校准方法及校准系统

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2009300104A (ja) * 2008-06-10 2009-12-24 Mitsubishi Electric Corp 送電線事故点標定装置および送電線事故点標定方法
CN103049923A (zh) * 2012-12-10 2013-04-17 深圳先进技术研究院 磁共振快速成像的方法
EP2713172A1 (en) * 2012-09-26 2014-04-02 Rheinisch-Westfälisch-Technische Hochschule Aachen Measurement apparatus for electricity distribution grids
CN108132013A (zh) * 2017-12-22 2018-06-08 交通运输部科学研究院 一种基于组合差分gnss的桥梁形变多频动态分析方法
CN114781196A (zh) * 2022-06-20 2022-07-22 四川大学 一种基于稀疏采集模型的谐波检测方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2009300104A (ja) * 2008-06-10 2009-12-24 Mitsubishi Electric Corp 送電線事故点標定装置および送電線事故点標定方法
EP2713172A1 (en) * 2012-09-26 2014-04-02 Rheinisch-Westfälisch-Technische Hochschule Aachen Measurement apparatus for electricity distribution grids
CN103049923A (zh) * 2012-12-10 2013-04-17 深圳先进技术研究院 磁共振快速成像的方法
CN108132013A (zh) * 2017-12-22 2018-06-08 交通运输部科学研究院 一种基于组合差分gnss的桥梁形变多频动态分析方法
CN114781196A (zh) * 2022-06-20 2022-07-22 四川大学 一种基于稀疏采集模型的谐波检测方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
AOBING CHI: "A Bregman-Split-Based Compressive Sensing Method for Dynamic Harmonic Estimation", 《ENTROPY》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117948872A (zh) * 2024-03-27 2024-04-30 常州同惠电子股份有限公司 一种基于lm法的电感式位置传感器校准方法及校准系统
CN117948872B (zh) * 2024-03-27 2024-06-04 常州同惠电子股份有限公司 一种基于lm法的电感式位置传感器校准方法及校准系统

Similar Documents

Publication Publication Date Title
Guo et al. Simplified algorithms for removal of the effect of exponentially decaying DC-offset on the Fourier algorithm
Zhan et al. Dynamic single-phase synchronized phase and frequency estimation at the distribution level
Zygarlicki et al. A reduced Prony's method in power-quality analysis—parameters selection
CN108037361B (zh) 一种基于滑动窗dft的高精度谐波参数估计方法
CN107102255B (zh) 单一adc采集通道动态特性测试方法
CN109409762A (zh) 基于改进数据关联分析的谐波责任划分方法
CN114781196A (zh) 一种基于稀疏采集模型的谐波检测方法
Lima et al. A novel method for power quality multiple disturbance decomposition based on Independent Component Analysis
CN109669072B (zh) 一种配电网的自适应同步相量量测方法
Abbas A new fast algorithm to estimate real-time phasors using adaptive signal processing
CN115032453A (zh) 一种多频动态相量测量方法
Petrović et al. Computational effective modified Newton–Raphson algorithm for power harmonics parameters estimation
Liu A wavelet based model for on-line tracking of power system harmonics using Kalman filtering
Zamora-Mendez et al. Model-based synchrophasor estimation by exploiting the eigensystem realization approach
Betancourt et al. A spatio-temporal processing Padé approach for visualizing harmonic distortion propagation on electrical networks
Tao et al. A robust parametric method for power harmonic estimation based on M-estimators
CN110112757B (zh) 基于sure小波消噪和改进hht的低频振荡分析方法
Shetty et al. Performance analysis of data-driven techniques for detection and identification of low frequency oscillations in multimachine power system
Yi et al. Instantaneous harmonic decomposition technique for three‐phase current based on multiple reference coordinates
CN110647720A (zh) 一种嵌入式平台下非平稳信号电能计量的方法
CN110320400A (zh) 准同步采样和改进能量算子的电压闪变包络参数提取方法
CN115219787A (zh) 基于改进矩阵束的电网相量移动测量方法、系统及介质
Ray et al. Improved recursive newton type algorithm based power system frequency estimation
Sundaram et al. Assessment of power quality disturbances in the distribution system using Kalman filter and fuzzy expert system
Xiong et al. A Multiple-Frequency-Taylor-Model based estimator for dynamic synchrophasor considering decaying DC component

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
RJ01 Rejection of invention patent application after publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20220909