现场测试参考值计算方法、系统、设备及存储介质

1.本发明涉及同步相量测量技术领域，尤其涉及一种现场测试参考值计算方法、系统、设备及存储介质。


背景技术：

2.同步相量量测装置(phasor measurement uni，pmu)因其快速性和同步性是监测电力系统动态过程的有效工具之一，国内目前已经安装了约4000台pmu，然而随着大规模新能源并网和电力电子装置的接入，这些pmu的数据质量是否能满足现在电力系统的动态监测有待研究。
3.近些年也发生过早期在现场安装的未经测试过的pmu，因高频谐波和频率偏差的影响进而导致低频振荡误报的事件，且现场信号也呈现出了更多的动态变化，因此有必要对现场的pmu进行现场测试与校准，以满足电网监测的正确性和可靠性。现有的测试校准系统主要分为两类，第一类是基于高精度信号源的测试系统，当信号源变为现场信号时，该类测试系统失去参考值，无法进行比对分析，因此不适用于现场测试。第二类是基于高精度校准器的测试系统，该类测试系统同时将现场信号输入给待测pmu和校准器，以校准器的量测值作为参考值，适用于为现场测试提供参考值，而校准器的参考值计算方法的精度需到达标准要求的10倍及以上。目前现有的校准器算法有频域算法和时域算法两种，频域算法主要是基于离散傅里叶变化的算法，该类算法计算量小，对稳态信号有较好的测量精度，然而对于动态信号，测量精度不满足校准要求。时域算法主要基于最小二乘法的算法，该类算法测量精度高，无论是稳态信号还是动态信号，都具有较好的测量精度，且满足校准要求。由于基于最小二乘法的方法计算复杂度较高，导致计算效率较低，因此，研究新的计算方法，在不失精度的前提下，减少计算量，提高计算效率，使算法可以在硬件中稳定运行，对pmu现场测试具有十分重要的作用，但现有技术中缺乏这样的解决方案。


技术实现要素：

4.本发明的目的是提供一种现场测试参考值计算方法、系统、设备及存储介质，可以在 pmu现场测试校准时，提供准确的参考值，保障pmu现场测试校准的精度。
5.本发明的目的是通过以下技术方案实现的：
6.一种现场测试参考值计算方法，包括：
7.建立拟合模型对现场信号进行拟合，并构建出具有正交性的拟合参数矩阵；
8.分析不同求解方法的复杂度，以及不同窗长和拟合阶数对测量精度和计算时间的影响，选取复杂度最低的求解方法以及最优窗长和拟合阶数，对所述具有正交性的拟合参数矩阵进行分解，迭代求解出待定拟合的参数矩阵；
9.利用获得的所述待定拟合的参数矩阵计算现场测试参考值。
10.一种现场测试参考值计算系统，包括：
11.信号拟合与拟合参数矩阵构建单元，用于建立拟合模型对现场信号进行拟合，并
构建出具有正交性的拟合参数矩阵；
12.求解方法及最优参数选取与迭代求解单元，用于分析不同求解方法的复杂度，以及不同窗长和拟合阶数对测量精度和计算时间的影响，选取复杂度最低的求解方法以及最优窗长和拟合阶数，对所述具有正交性的拟合参数矩阵进行分解，迭代求解出待定拟合的参数矩阵；
13.现场测试参考值计算单元，用于利用获得的所述待定拟合的参数矩阵计算现场测试参考值。
14.一种处理设备，包括：一个或多个处理器；存储器，用于存储一个或多个程序；
15.其中，当所述一个或多个程序被所述一个或多个处理器执行时，使得所述一个或多个处理器实现前述的方法。
16.一种可读存储介质，存储有计算机程序，当计算机程序被处理器执行时实现前述的方法。
17.由上述本发明提供的技术方案可以看出，可以准确计算参考值，精度高于标准要求 10倍以上，从而在pmu现场测试校准时提供准确的参考值；并且，计算效率较高，可大幅减轻硬件计算负担。
附图说明
18.为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。
19.图1为本发明实施例提供的一种现场测试参考值计算方法的流程图；
20.图2为本发明实施例提供的t坐标下波形图；
21.图3为本发明实施例提供的τ坐标下波形图；
22.图4为本发明实施例提供的不同窗长和拟合阶数下的幅值误差图；
23.图5为本发明实施例提供的不同窗长和拟合阶数下的相角误差图；
24.图6为本发明实施例提供的各求解方法复杂度对比图；
25.图7为本发明实施例提供的传统单项式拟合矩阵零元素分布图；
26.图8为本发明实施例提供的legendre多项式拟合矩阵零元素分布图；
27.图9为本发明实施例提供的相同窗长不同拟合阶数下的计算时间对比图；
28.图10为本发明实施例提供的一种现场测试参考值计算系统的示意图；
29.图11为本发明实施例提供的一种处理设备的示意图。
具体实施方式
30.下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。
31.首先对本文中可能使用的术语进行如下说明：
32.术语“包括”、“包含”、“含有”、“具有”或其它类似语义的描述，应被解释为非排它性的包括。例如：包括某技术特征要素(如原料、组分、成分、载体、剂型、材料、尺寸、零件、部件、机构、装置、步骤、工序、方法、反应条件、加工条件、参数、算法、信号、数据、产品或制品等)，应被解释为不仅包括明确列出的某技术特征要素，还可以包括未明确列出的本领域公知的其它技术特征要素。
33.下面对本发明所提供的一种现场测试参考值计算方法、系统、设备及存储介质进行详细描述。本发明实施例中未作详细描述的内容属于本领域专业技术人员公知的现有技术。本发明实施例中未注明具体条件者，按照本领域常规条件或制造商建议的条件进行。
34.实施例一
35.本发明实施例提供一种现场测试参考值计算方法，如图1所示，该方法主要包括如下步骤：
36.步骤1、建立拟合模型对现场信号进行拟合，并构建出具有正交性的拟合参数矩阵。
37.本发明实施例中，建立基于legendre多项式拟合模型对现场信号进行拟合，表示为：
[0038][0039]
其中，干扰信号已被滤除，x(t)为基波信号，li(t)为i阶legendre多项式，ai和bi均为 legendre多项式系数，n为拟合阶数，i＝0,1,...,n，f0为基波频率，t为时间。
[0040]
将基于legendre多项式拟合模型对现场信号进行拟合的表达式转换为矩阵形式：
[0041][0042]
其中，x(t)为多个基波信号x(t)构成的列向量，它表示时间窗内采样序列，h为待定拟合的参数矩阵；q＝[a0,a1…an
]
t
，r＝[b0,b1…bn
]
t
，正体t为转置符号；时间t为列向量n为窗内采样点的数量，fs为采样频率，为窗长；
[0043]
p为n
×
(2n 2)阶的具有正交性的拟合参数矩阵，表示为：
[0044][0045]
其中，c＝cos(2πf
l
τ)，s＝sin(2πf
l
τ)，定义c与s主要是为了简化表达式，不同角标的τ为不同的时间t进行坐标变换后的值，f
l
为基波频率f0进行坐标变换后的值，表示为：
[0046]
τ＝(2t-a-b)/(b-a)
[0047]fl
＝f0(b-a)/2＝t/2t0[0048]
其中，a与b分别为时间窗的首末两端(即分别为0与)，斜体t为窗长(即)，t0为基频周期。
[0049]
上式转换公式中a与b为固定值，不同的时间t带入第一个式子获得不同角标的τ，以τ0为例，带入中第一个值进行计算，其余角标的τ也是如此计算。
[0050]
如图2～图3所示，分别为t坐标下波形图，τ坐标下波形图，通过上述坐标变换，可以保证变换前后坐标一一对应的关系不会变化。
[0051]
由于legendre多项式中t∈[-1,1]上体现出正交性即：
[0052][0053]
其中，m＝0,1,...,n,s＝0,1,...,n；
[0054]
该正交性质结合拟合参数矩阵p，可知前述p的矩阵式子中：
[0055][0056][0057][0058]
其中，w表示正整数，《a,b》表示两个向量的内积，a、b分别代表上述三个等式左侧括号内的两个向量，τ表示p中所有不同角标的τ构成的向量。
[0059]
其正交性表明，拟合参数矩阵p中的列向量两两正交，即拟合参数矩阵p具有部分正交性。正是由于该性质，为后续计算效率的提高提供了基础。
[0060]
步骤2、分析不同求解方法的复杂度，以及不同窗长和拟合阶数对测量精度和计算时间的影响，选取复杂度最低的求解方法以及最优窗长和拟合阶数，对所述具有正交性的拟合参数矩阵进行分解，迭代求解出待定拟合的参数矩阵。
[0061]
本发明实施例中，通过测试不同分析窗长和不同拟合阶数组合下的测试精度与计算时间；选出满足计算时间要求下测试精度最高时的窗长和拟合阶数，作为最优窗长和拟合阶数。主要原理如下：
[0062]
对于稳态信号，选取较短的窗长和较低的阶数即可准确测量，但由于现场信号呈现复杂的动态变化，因此需考虑窗长和拟合阶数对计算精度的影响。如图4和图5所示，以调制信号为例，分别分析了窗长和拟合阶数对计算精度的影响，两个图中，t
l
表示窗长，以秒为单元；可以看出当窗长为3个周波和4个周波时，随着窗长和拟合阶数的增加，测量误差整体随之减小，当窗长为2个周波时，测量误差并不会随着拟合阶数增加而减小，原因是较短的窗长不能反映信号的动态变化。通过仔细对比可发现，窗长为3个周波，拟合阶数为8阶时，测量误差最小，且计算时间也较少。当然，此处所提供的具体窗长数值与拟合阶数数值仅为举例，并非构成限制，在实际应用中，用户可以根据实际需求或者经验或者测试结果调整相关数值。
[0063]
本发明实时中，在选择具体求解方式时，考虑了不同求解方法的复杂度，为了简化，以n
×
n阶矩阵p为例可得到同样的结果，例如：
[0064]
最小二乘法复杂度为n3 2nn2 nn阶；
[0065]
svd分解复杂度为20nn2 2n2 n2阶；
[0066]
qr分解复杂度为2nn2 2n2 n2阶；
[0067]
cholesky分解复杂度为1/3n3 (n 3)n2 nn阶。
[0068]
通过对比不同求解方法的复杂度，优选使用cholesky分解迭代求解待定拟合的参数矩阵，主要步骤如下：
[0069]
1)由于cholesky分解只能对方阵进行分解，因此，对矩形形式的x(t)≈ph进行转换，获得如下等式：
[0070]
p
t
x(t)＝p
t
ph
[0071]
2)设置坐标变换后基波频率f
l
的初值f
l0
，使具有正交性的拟合参数矩阵p为常系数矩阵，并使用cholesky分解方法对上述等式中的p
t
p部分进行分解，转换为两个互为转置的下三角矩阵的乘积，得到：
[0072]
p
t
x(t)＝ll
th[0073]
其中，l为下三角矩阵，它是n阶可逆矩阵。
[0074]
本发明实施例中，由于是在离散时间采样，因此拟合参数矩阵p是满秩矩阵，所以矩阵p
t
p具有正定性，由正定矩阵的性质可知，其可以进行cholesky分解，分解成两个互为转置的下三角矩阵的乘积，分解出的，即p
t
p＝ll
t
，其中l为下三角矩阵。再由各个方法的复杂度可知，有两个可变参数，分别为拟合阶数和采样点个数，因此设置不同变量可得如图6所示的复杂度对比，图6中，从下至少四个部分依次对应cholesky分解、最小二乘法、qr分解、svd分解的复杂度，可见出基于cholesky分解的求解方法计算复杂度最低且低于传统的最小二乘法。
[0075]
同时cholesky分解本身会占用一定的计算时间，其是由lu分解简化而来，而lu分解的本质是高斯消元的过程，因此矩阵零元素含量越多，高斯消元的过程越快，cholesky 分解的过程也越快，由步骤1中legendre多项式构成的正交拟合矩阵可知，拟合参数矩阵 p由于正交性，其零元素的含量远小于传统单项式构成的拟合矩阵，零元素具体分布如图 7和图8所示，证明了理论分析的可行性。
[0076]
3)推导求解出待定拟合的参数矩阵h：
[0077]
h＝(l
t
)-1
l-1
p
t
x(t)。
[0078]
基于上述求解步骤，当已知频率f
l
，即可用上述方法求解出h，因此设坐标变换后基波频率f
l
初值为f
l0
，并带入具有正交性的拟合参数矩阵p中，作为其中的f
l
。以基频 50hz为例，窗长3个周波为例，频率初值f
l0
为1.5，进而可通过迭代计算的形式，直到满足要求的精度或超出最大迭代次数为止。需要说明的是，此处所涉及的具体数值仅为举例，并非构成限制，在实际应用中，用户可根据需求或者经验设定相应的频率初值大小。
[0079]
步骤3、利用获得的所述待定拟合的参数矩阵计算现场测试参考值。
[0080]
示例性的，场测试参考值的类型可以包括：幅值、相角与频率。
[0081]
得到最终的矩阵h后，可求解出时间窗内任意时刻的同步相量参考值，经仿真和实验验证后，把时标打在窗首，测量结果最准确，对应的τ0＝-1，则同步相量的幅值相角分别为：
[0082][0083][0084]
频率迭代修正公式为：
[0085][0086]
其中，t表示迭代次数，t＝1时，f
l(0)
为设置的坐标变换后基波频率的初值(也即前文定义的f
l0
)，r'(-1)与p'(-1)均表示求导数，δf
l(t-1)
为第t-1次迭代的频率修正量，a(τ0)为幅值，为相角，为相角，a0～an，b0～bn可由矩阵h获得。
[0087]
最后一次迭代获得频率f
lt
，由于f
lt
为变换后的频率，所以求解频率f要对f
lt
作反变换，变换公式为：
[0088][0089]
为了体现本发明在不失精度的前提下，大大提高了计算效率，对比了传统最小二乘法 (twls)与本发明上述方法(cal)在相同窗长不同拟合阶数下的计算时间，如图9所示，其中计算时间为计算一次相量所用的时间。从图9可以看出传统最小二乘法在拟合阶数过高时没有实际值，原因是矩阵求逆的计算量大，硬件计算负担大，导致计算失败。而所提方法在高阶时仍能保持良好的计算速度。
[0090]
表1与表2为pmu标准规定的各种测试类型的测量误差，表1与表2中cal为本发明上述方法，std为pmu标准规定的误差；ae表示幅值误差，pe表示相角误差，fe表示频率误差，rfe表示rocof误差。
[0091][0092]
表1最大幅值和相角误差
[0093][0094]
表2最大频率和rocof误差
[0095]
从以上两个表中所示结果可见，本发明计算的参考值均满足标准要求的10倍以上，能够保障pmu现场测试校准的精度。
[0096]
实施例二
[0097]
本发明还提供一种现场测试参考值计算系统，其主要基于前述实施例一提供的方法实现，如图10所示，该系统主要包括：
[0098]
信号拟合与拟合参数矩阵构建单元，用于建立拟合模型对现场信号进行拟合，并
构建出具有正交性的拟合参数矩阵；
[0099]
求解方法及最优参数选取与迭代求解单元，用于分析不同求解方法的复杂度，以及不同窗长和拟合阶数对测量精度和计算时间的影响，选取复杂度最低的求解方法以及最优窗长和拟合阶数，对所述具有正交性的拟合参数矩阵进行分解，迭代求解出待定拟合的参数矩阵；
[0100]
现场测试参考值计算单元，用于利用获得的所述待定拟合的参数矩阵计算现场测试参考值。
[0101]
所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，仅以上述各功能模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能模块完成，即将系统的内部结构划分成不同的功能模块，以完成以上描述的全部或者部分功能。
[0102]
以上各单元所涉及的相关原理在之前的实施例一中已经做了详细的说明，故不再赘述。
[0103]
实施例三
[0104]
本发明还提供一种处理设备，如图11所示，其主要包括：一个或多个处理器；存储器，用于存储一个或多个程序；其中，当所述一个或多个程序被所述一个或多个处理器执行时，使得所述一个或多个处理器实现前述实施例提供的方法。
[0105]
进一步的，所述处理设备还包括至少一个输入设备与至少一个输出设备；在所述处理设备中，处理器、存储器、输入设备、输出设备之间通过总线连接。
[0106]
本发明实施例中，所述存储器、输入设备与输出设备的具体类型不做限定；例如：
[0107]
输入设备可以为触摸屏、图像采集设备、物理按键或者鼠标等；
[0108]
输出设备可以为显示终端；
[0109]
存储器可以为随机存取存储器(random access memory，ram)，也可为非不稳定的存储器(non-volatile memory)，例如磁盘存储器。
[0110]
实施例四
[0111]
本发明还提供一种可读存储介质，存储有计算机程序，当计算机程序被处理器执行时实现前述实施例提供的方法。
[0112]
本发明实施例中可读存储介质作为计算机可读存储介质，可以设置于前述处理设备中，例如，作为处理设备中的存储器。此外，所述可读存储介质也可以是u盘、移动硬盘、只读存储器(read-only memory，rom)、磁碟或者光盘等各种可以存储程序代码的介质。
[0113]
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。