一种基于时域Hermite插值的准同步谐波分析装置及方法

一种基于时域hermite插值的准同步谐波分析装置及方法
技术领域
1.本发明涉及谐波分析技术领域，具体为一种基于时域hermite插值的准同步谐波分析装置及方法。


背景技术：

2.随着大规模可再生能源的接入及负荷侧的再电气化过程，大量的特性各异的电源、负荷、储能等装备以电力电子为接口接入现有电力系统，使电力系统向着高比例可再生能源和高比例电力电子设备趋势快速发展。由此带来的谐波畸变会对系统中的设备造成严重的危害，对电力系统的安全稳定构成了极大的威胁。因此，准确地检测电力系统信号中地谐波成分具有重要意义。
3.目前，关于谐波的检测方法主要有三种。第一种是基于离散傅里叶变换(discrete fourier transform,dft)的方法。在同步采样的条件下，dft具有较高的精度。但是，电力系统的基频可能会因负载变化而产生波动以至于同步采样无法实现，导致dft在检测谐波时不可避免地会遇到频谱泄漏和栅栏效应等问题。第二种是现代谱估计方法，这些方法从另一个角度检测谐波分量，从根本上避免了dft固有的缺陷，具有较高的检测精度。然而这些方法获得精确结果的前提是准确估计信号中的频率分量个数。并且这些方法对噪声敏感，计算量大。因此，它的应用非常有限。第三种方法是现代类方法。这种方法可以在一定程度上消除频谱泄漏问题，但收敛速度和稳定性较差，应用也受到一定的限制，需要进一步研究和改进。在以上三种方法中，基于dft的谐波检测方法是最成熟且常用的，尽管它在非同步采样条件下存在频谱泄露和栅栏现象。
4.现有技术的缺点如下：
5.1)基于dft的传统方法，在检测谐波时不可避免地会遇到频谱泄漏和栅栏效应等问题。尽管通过hanning窗，nuttall窗等加窗插值算法能一定程度上改善频谱泄漏和栅栏效应，但这些窗难以同时满足主瓣能量集中与旁瓣衰减快等特点，并且非稳态下的加窗插值算法仍然存在修正公式求解复杂，计算精度难以提高等缺点。
6.2)现代谱估计方法，在检测谐波时需要一定的先验信息，其获得精确谐波分析结果的前提是已知信号中所含谐波分量的个数，并且这些方法对噪声敏感，计算量大。因此，它的应用非常有限。
7.3)现代类方法，在谐波分析时的收敛速度和稳定性较差，应用也受到一定的限制，需要进一步研究和改进。


技术实现要素：

8.针对上述问题，本发明的目的在于提供一种基于时域hermite插值的准同步谐波分析装置及方法，通过构建准同步采样序列以获得全周期截断信号解决频谱泄露和栅栏现象，有效提高了谐波分析的准确性。技术方案如下：
9.一种基于时域hermite插值的准同步谐波分析方法，包括以下步骤：
10.步骤1：利用电力系统采样信号构造轨迹矩阵并进行奇异值分解；
11.步骤2：将奇异值分奇偶位置进行重新排序，计算奇异值均值序列，并将奇异值均值序列中各元素及其序号归一化；再计算差值曲线，确定有效阶次，并根据有效阶次重构信号，完成去噪处理，得到去噪信号；
12.步骤3：基于变分模态分解提取去噪信号中的基波分量；
13.步骤4：利用穿越时间后的采样点和穿越时间前的采样点，可获得估计的穿越时间，从而得到基波周期；
14.步骤5：基于数值微分计算各个采样点的导数；
15.步骤6：分别计算各个子区间上的三次hermite插值多项式，并计算准同步采样时刻，将准同步采样时刻带入对应的三次hermite插值多项式获得准同步采样数据；
16.步骤7：对准同步采样数据应用dft进行谐波分析。
17.进一步的，所述步骤1中采样信号为电压信号时，表示为：
18.uo(kts)＝u
p
(kts)+n(kts)
ꢀꢀ
(1)
19.式中，u
p
(kts)表示去噪信号序列；n(kts)表示噪声序列；k为序列索引，ts为采样周期；
20.考虑一个长度为n的原始电压采样序列uo＝[u
o1 u
o2
ꢀ…ꢀuon
]，将其分解为k＝n-l+1个长度为l的向量，表示为：
[0021]hi
＝[u
oi u
o(i+1)
ꢀ…ꢀuo(i+l-1)
]
t
,1≤i≤k
ꢀꢀ
(2)
[0022]
将这些向量组成轨迹矩阵h：
[0023][0024]
令矩阵s＝hh
t
，计算矩阵s的l个特征值λ1≥λ2≥
…
≥λ
l
≥0，以及其所对应的标准正交向量u1,u2,
…
,u
l
；
[0025]
令d＝rank(h)，以及则h表示为：
[0026]
h＝h1+h2+
…
+hdꢀꢀ
(4)
[0027]
式中，1≤j≤d，是轨迹矩阵h的奇异值，uj和vj是对应于的左特征向量和右特征向量；d为轨迹矩阵h的秩。
[0028]
更进一步的，所述步骤2具体为：
[0029]
步骤2.1：将奇异值分奇偶位置进行重新排序，并舍弃掉奇数序列最后一个奇异值：
[0030][0031]
式中，σ
odd
表示奇数序列，σ
even
表示偶数序列；
[0032]
步骤2.2：计算奇偶序列的均值，得到奇异值均值序列：
[0033][0034]
令x为奇异值均值序列σm中各元素的序号，即x＝[1 2
…
end]；然后将x和σm归一化：
[0035][0036]
式中，和x
ai
分别为第ai个归一化的序号和第ai个未归一化的序号；和分别为第ai个归一化的奇异值均值和第ai个未归一化的奇异值均值；
[0037]
步骤2.3：将归一化曲线的“肘部”转化为“膝部”[0038][0039]
式中，为转换后的第ai个归一化的奇异值均值，{σn}为归一化的奇异值均值序列；
[0040]
曲线上离σs＝xn最远的点即为曲线的膝点，根据相似三角形的原理，将曲线上的点到σs＝xn的距离转换为该点横纵坐标的差值；则计算差值曲线：
[0041][0042]
式中，为第ai个差值；
[0043]
步骤2.4：假设差值曲线中的最大值是{σ
di
}中的第m个点，则有效阶次r＝2
×
(m-1)；
[0044]
步骤2.5：取式(4)中的前r个分量相加，得到矩阵h
i1
；
[0045]
假设h
i1
中的元素为h
i1(ki,kj)
，令l
*
＝min(l,k),k
*
＝max(l,k),n＝l+k-1，则按式(10)重构得到去噪信号u
p
(kts)：
[0046][0047]
式中，l
*
为向量长度l和向量数量k中较小的数值；k
*
为向量长度l和向量数量k中较大的数值；u
pk
为为去噪信号序列u
p
(kts)中的第k个元素。
[0048]
更进一步的，所述步骤3具体为：
[0049]
步骤3.1：将原始的去噪信号u
p
(kts)分解为z个有限带宽的模态函数，并采用希尔伯特变换对每个模态函数uz(t)进行信号解析，得到单边频谱；
[0050]
步骤3.2：在每个模态函数uz(t)中加入指数项对中心频率ωz进行修正，并将每个模态的频谱调制至基频带宽；
[0051]
步骤3.3：通过计算解调信号的梯度二范数估计每个模态函数uz(t)的基频带宽，构建带约束条件的变分模型：
[0052][0053]
式中，{uz}为模态集合；{ωz}为中心频率集合；为偏导运算符；δ(t)为单位脉冲函数；u
p
为去噪信号；
[0054]
步骤3.4：将问题转化为非约束条件的变分模型，简化计算，引入lagrange算子λ(t)和二次惩罚项α；得如下表达式：
[0055][0056]
式中，u
p
(t)为输入的去噪信号；uz(t)为模态函数；
[0057]
步骤3.5：采用交替方向乘子法求解式(12)，更新直至获得变分模型的最优解；
[0058][0059][0060]
式中，为当前剩余量的维纳滤波；为u
p
(t)的傅里叶变换；为ui(t)的傅里叶变换；为λ(t)的傅里叶变换；为当前模态函数功率谱的重心；
[0061]
步骤3.6：对最终获得的进行傅里叶逆变换得到分解出的z个分量，其中频率最接近50hz的即为基波分量uf(kts)。
[0062]
更进一步的，所述步骤4具体为：
[0063]
利用穿越时间后的3个采样点(kts，uf(kts))，((k+1)ts，uf((k+1)ts))，((k+2)ts，uf((k+2)ts))和穿越时间前的2个采样点((k-1)ts，uf((k-1)ts))，((k-2)ts，uf((k-2)ts))，通过下式获得估计的穿越时间，从而得到基波周期：
[0064][0065]
式中，z为穿越时间所对应的去噪信号的值；f(
·
)表示相应采样点的均差；
[0066]
将通过式(15)获得的时间序列表示为tc＝[t
c1 t
c2
ꢀ…ꢀ
t
c(end)
]，因此，第aj个基波周期的时间则表示为：
[0067]
t
pi
＝t
c(i+1)-t
ci
ꢀꢀ
(16)
[0068]
式中，和分别为第bj+1个和第bj个穿越时间；
[0069]
则每个基波周期的持续时间构成的序列为t
p
＝[t
p1 t
p2
ꢀ…ꢀ
t
p(end)
]，据此确定hermite插值重采样时刻。
[0070]
更进一步的，所述步骤5具体为：
[0071]
步骤5.1：对于由n个数据构成的去噪信号u
p
(kts)，将其对应的采样时刻序列tn＝[t
s 2tsꢀ…ꢀ
nts]看作将区间[t
s nts]等分成n-1份；由于：
[0072][0073]
式中，d(t)代表u
p
(t)的导数；
[0074]
步骤5.2：对式(17)的右边利用simpson公式，得到：
[0075][0076]
步骤5.3：对于式(18)，记mk为d(kts)的近似值，如果已知边界条件m1＝d(ts)和mn＝d(nts)，则有：
[0077][0078]
对于边界条件m1＝d(ts)和mn＝d(nts)，通过非中点的三点公式计算；
[0079]
步骤5.4：通过点(ts,u
p1
),(2ts,u
p2
),(3ts,u
p3
)的二次插值多项式表示为：
[0080][0081]
式中，u
p1
、u
p2
、u
p3
分别表示去噪信号u
p
(kts)中的第1个、第2个和第3个元素；
[0082]
则(ts,u
p1
)处导数值为：
[0083][0084]
将上式化简得：
[0085][0086]
式中，δk为kts处的差商，由下式计算：
[0087][0088]
同理：
[0089][0090]
将(22)和(24)代入(19)即求得所有采样点的导数。
[0091]
更进一步的，所述步骤6中，每个子区间kts≤t≤(k+1)ts上的插值多项式pk(t)通过下式计算：
[0092][0093]
令则第mi个准同步采样时刻通过下式计算：
[0094][0095]
式中，nr为期望的每周期采样点数；为第q个基波周期的时间；为第q个穿越时间。
[0096]
一种基于时域hermite插值的准同步谐波分析装置，包括：
[0097]
轨迹矩阵构造模块：利用电力系统采样信号构造轨迹矩阵并进行奇异值分解；
[0098]
有效阶次重构模块：将奇异值分奇偶位置进行重新排序，计算奇异值均值序列，并将奇异值均值序列中各元素及其序号归一化；再计算差值曲线，确定有效阶次，并根据有效阶次重构信号，完成去噪处理，得到去噪信号；
[0099]
基波分量提取模块：基于变分模态分解提取去噪信号中的基波分量；
[0100]
基波周期计算模块：利用穿越时间后的采样点和穿越时间前的采样点，可获得估计的穿越时间，从而得到基波周期；
[0101]
采样点导数计算模块：基于数值微分计算各个采样点的导数；
[0102]
准同步采样数据计算模块：分别计算各个子区间上的三次hermite插值多项式，并计算准同步采样时刻，将准同步采样时刻带入对应的三次hermite插值多项式获得准同步采样数据；
[0103]
谐波分析模块：对准同步采样数据应用dft进行谐波分析。
[0104]
本发明的有益效果是：本发明提出了一种基于时域hermite插值的准同步谐波分析算法，该方法以电力系统采样数据为基础，将改进的奇异谱分析方法应用于采用数据进行去噪，通过改进奇异谱分析方法提高了去噪效果，增强了谐波分析的抗噪性；采用vmd和牛顿插值获得准确的基波周期，通过数值微分方法获得采样点对应的导数，从而获得分段三次hermite插值多项式，最后带入准同步采样时刻获得准同步采样数据，达到消除dft运算时的频谱泄露和栅栏现象的目的。并且，该方法以电力系统采样数据为基础，并不需要先验信息。
附图说明
[0105]
图1为本发明基于时域hermite插值的准同步谐波分析方法的流程图。
具体实施方式
[0106]
下面结合附图和具体实施例对本发明做进一步详细说明。
[0107]
电力系统的谐波畸变会对系统中的设备造成严重的危害，对电力系统的安全稳定构成了极大的威胁，为了准确检测电力系统信号中的谐波，本发明提出一种基于时域hermite插值的准同步谐波分析方法，基本流程图如图1所示，具体步骤如下：
[0108]
s1：利用电力系统采样信号构造轨迹矩阵h并进行奇异值分解。
[0109]
电力系统中以等间隔时间采样的离散时间样本对应一个一维时间序列。以电压信
号为例，可以表示为：
[0110]
uo(kts)＝u
p
(kts)+n(kts)
ꢀꢀ
(27)
[0111]
考虑一个长度为n的原始电压采样序列uo＝[u
o1 u
o2
ꢀ…ꢀuon
]，将其分解为k＝n-l+1个长度为l的向量，l的选择一般为n/2，如果n为奇数的话，l的值可取(n+1)/2。
[0112]hi
＝[u
oi u
o(i+1)
ꢀ…ꢀuo(i+l-1)
]
t
,1≤i≤k
ꢀꢀ
(28)
[0113]
将这些向量组成轨迹矩阵h：
[0114][0115]
令s＝hh
t
,计算s的l个特征值λ1≥λ2≥
…
≥λ
l
≥0和其所对应的标准正交向量u1,u2,
…
,u
l
.
[0116]
令d＝rank(h)以及则h可以表示为：
[0117]
h＝h1+h2+
…
+hdꢀꢀ
(30)
[0118]
式中，1≤j≤d，是轨迹矩阵h的奇异值，uj和vj是对应于的左特征向量和右特征向量；d为轨迹矩阵h的秩。
[0119]
s2：确定有效阶次r，并根据r重构信号，完成去噪处理。
[0120]
首先，将奇异值分奇偶位置进行重新排序。如果奇异值数目为奇数，重新排序后会导致奇数序列比偶数序列多一个，此时为了方便后面计算，可舍弃掉最后一个奇异值：
[0121][0122]
计算奇偶序列的均值，得到奇异值均值序列：
[0123][0124]
令x为奇异值均值序列σm中各元素的序号，即x＝[1 2
ꢀ…ꢀ
end]。然后将x和σm归一化：
[0125][0126]
式中，和x
ai
分别为第ai个归一化的序号和第ai个未归一化的序号；和分别为第ai个归一化的奇异值均值和第ai个未归一化的奇异值均值。
[0127]
接着，将归一化曲线的“肘部”转化为“膝部”.
[0128][0129]
式中，为转换后的第ai个归一化的奇异值均值，{σn}为归一化的奇异值均值序列、
[0130]
由于归一化曲线一定经过(0，0)和(1，1)，因此曲线上离σs＝xn最远的点即为曲线的膝点。根据相似三角形的原理，曲线上的点到σs＝xn的距离可以转换为该点横纵坐标的差值。因此，接着计算差值曲线：
[0131][0132]
式中，为第ai个差值。
[0133]
假设差值曲线中的最大值是中的第m个点，则有效阶次r＝2
×
(m-1)。
[0134]
取式(4)中的前r个分量相加，得到矩阵h
i1
。假设h
i1
中的元素为h
i1(ki,kj)
，令l
*
＝min(l,k),k
*
＝max(l,k),n＝l+k-1，则可按式(10)重构得到去噪信号u
p
(kts)
[0135][0136]
s3：基于vmd提取去噪信号u
p
中的基波分量uf。
[0137]
vmd(variational mode decomposition变分模态分解)是一种信号分解方法，该方法在获取分解分量的过程中通过迭代搜寻变分模型最优解来确定每个分量的频率中心和带宽，从而能够自适应地实现信号各分量的有效分离。该方法的具体过程如下：
[0138]
(1)将原始的去噪信号u
p
(kts)分解为z个有限带宽模态函数，并采用希尔伯特变换对每个模态函数uz(t)进行信号解析，得到单边频谱。
[0139]
(2)每个模态函数uz(t)中加入指数项对中心频率ωz进行修正，并将每个模态的频谱调制至基频带宽。
[0140]
(3)通过计算解调信号的梯度二范数估计每个模态函数uz(t)的基频带宽，构建带约束条件的变分模型如式(11)。
[0141][0142]
式中，{uz}为模态集合；{ωz}为中心频率集合；为偏导运算符；δ(t)为单位脉冲函数；u
p
为去噪信号
[0143]
(4)为了将问题转化为非约束条件的变分模型，进而简化计算，引入lagrange算子λ(t)和二次惩罚项α；可得如下表达式：
[0144][0145]
式中，u
p
(t)为输入的去噪信号；uz(t)为模态函数
[0146]
(5)采用交替方向乘子法求解式(12)，更新直至获得变分模型的最优解。
[0147][0148][0149]
式中，为当前剩余量的维纳滤波；为u
p
(t)的傅里叶变换；为ui(t)的傅里叶变换；为λ(t)的傅里叶变换；为当前模态函数功率谱的重心。
[0150]
对最终获得的进行傅里叶逆变换即可得到分解出的z个分量，其中频率最接近50hz的就是基波分量uf(kts)。
[0151]
s4：确定基波周期。
[0152]
对于周期信号，如果满足一定条件的点在一个周期内只出现一次，则相邻点之间的时间间隔可视为信号的周期。例如正弦曲线上与y＝z相交且斜率为负的点一周期只出现一次。如果能准确的找到这些点，就能确定准确的基波周期。但是由于信号采样的限制，很难直接确定曲线穿过阈值z的时间，但是该时刻前后的采样点我们可以很容易的获得。利用这些采样点，使用插值算法构造这一小段时间内信号的解析表达式，即可获得近似的穿越时间。
[0153]
在估计基波周期过程中，牛顿四次插值多项式的误差很小，可以忽略不计。因此，利用穿越时间后的3个采样点(kts,uf(kts)),((k+1)ts,uf((k+1)ts)),((k+2)ts,uf((k+2)ts))和穿越时间前的2个采样点((k-1)ts,uf((k-1)ts)),((k-2)ts,uf((k-2)ts))，通过下式即可获得估计的穿越时间，从而得到基波周期。
[0154][0155]
通过式(15)获得的时间序列可以表示为tc＝[t
c1 t
c2
ꢀ…ꢀ
t
c(end)
]，因此，第aj个基波周期的时间则可表示为：
[0156][0157]
所以，t
p
＝[t
p1 t
p2
ꢀ…ꢀ
t
p(end)
]是每个基波周期的持续时间构成的序列，据此则可确定hermite插值重采样时刻。
[0158]
s5：基于数值微分计算各个采样点的导数。
[0159]
对于由n个数据构成的去噪信号u
p
(kts)而言，其对应的采样时刻序列tn＝[t
s 2tsꢀ…ꢀ
nts]，可以看成将区间[t
s nts]等分成n-1份。由于：
[0160][0161]
式中，d(t)代表u
p
(t)的导数。
[0162]
对式(17)的右边利用simpson公式，可得：
[0163][0164]
对于式(18)，记mk为d(kts)的近似值，如果已知边界条件m1＝d(ts)和mn＝d(nts)，则有：
[0165][0166]
对于边界条件m1＝d(ts)和mn＝d(nts)，可通过非中点的三点公式计算。
[0167]
通过点(ts,u
p1
),(2ts,u
p2
),(3ts,u
p3
)的二次插值多项式可表示为：
[0168][0169]
则(ts,u
p1
)处导数值为：
[0170][0171]
将上式化简可得：
[0172][0173]
式中，δk为kts处的差商，由下式计算：
[0174][0175]
同理：
[0176][0177]
将(22)和(24)代入(19)即可求得所有采样点的导数。
[0178]
s6：基于三次hermite插值多项式，计算准同步采样数据。
[0179]
计算出各点的导数后，就可以分别计算各个子区间上的三次hermite插值多项式了。每个子区间kts≤t≤(k+1)ts上的插值多项式pk(t)可通过下式计算:
[0180][0181]
令则第mi个准同步采样时刻通过下式计算：
[0182][0183]
式中，式中，nr为期望的每周期采样点数；为第q个基波周期的时间；为第q个穿越时间。
[0184]
将准同步采样时刻带入对应的三次hermite插值多项式即可获得准同步采样数据。
[0185]
s7：对准同步采样数据应用dft进行谐波分析。
[0186]
综上，通过改进奇异谱分析方法提高了去噪效果，增强了谐波分析的抗噪性，最终得到的准同步采样数据克服了进行dft运算时的频谱泄露和栅栏现象。