CN115508618B - 一种基于时域Hermite插值的准同步谐波分析装置及方法 - Google Patents
一种基于时域Hermite插值的准同步谐波分析装置及方法 Download PDFInfo
- Publication number
- CN115508618B CN115508618B CN202211253712.5A CN202211253712A CN115508618B CN 115508618 B CN115508618 B CN 115508618B CN 202211253712 A CN202211253712 A CN 202211253712A CN 115508618 B CN115508618 B CN 115508618B
- Authority
- CN
- China
- Prior art keywords
- sequence
- sampling
- quasi
- signal
- synchronous
- 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
Links
Classifications
-
- 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/13—Differential equations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
- G06F17/141—Discrete Fourier transforms
-
- 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
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
- Y02E40/00—Technologies for an efficient electrical power generation, transmission or distribution
- Y02E40/40—Arrangements for reducing harmonics
Abstract
本发明公开了一种基于时域Hermite插值的准同步谐波分析装置及方法,属于谐波分析技术领域,轨迹矩阵构造模块利用电力系统采样信号构造轨迹矩阵并进行奇异值分解;有效阶次重构模块确定有效阶次,并根据有效阶次重构信号,完成去噪处理;基波分量提取模块基于VMD提取去噪信号中的基波分量;基波周期计算模块确定基波周期;采样点导数计算模块基于数值微分计算各个采样点的导数;准同步采样数据计算模块基于三次Hermite插值多项式,计算准同步采样数据;谐波分析模块对准同步采样数据应用DFT进行谐波分析。本发明通过改进奇异谱分析方法提高了去噪效果,增强了谐波分析的抗噪性,最终得到的准同步采样数据克服了进行DFT运算时的频谱泄露和栅栏现象。
Description
技术领域
本发明涉及谐波分析技术领域,具体为一种基于时域Hermite插值的准同步谐波分析装置及方法。
背景技术
随着大规模可再生能源的接入及负荷侧的再电气化过程,大量的特性各异的电源、负荷、储能等装备以电力电子为接口接入现有电力系统,使电力系统向着高比例可再生能源和高比例电力电子设备趋势快速发展。由此带来的谐波畸变会对系统中的设备造成严重的危害,对电力系统的安全稳定构成了极大的威胁。因此,准确地检测电力系统信号中地谐波成分具有重要意义。
目前,关于谐波的检测方法主要有三种。第一种是基于离散傅里叶变换(DiscreteFourier Transform,DFT)的方法。在同步采样的条件下,DFT具有较高的精度。但是,电力系统的基频可能会因负载变化而产生波动以至于同步采样无法实现,导致DFT在检测谐波时不可避免地会遇到频谱泄漏和栅栏效应等问题。第二种是现代谱估计方法,这些方法从另一个角度检测谐波分量,从根本上避免了DFT固有的缺陷,具有较高的检测精度。然而这些方法获得精确结果的前提是准确估计信号中的频率分量个数。并且这些方法对噪声敏感,计算量大。因此,它的应用非常有限。第三种方法是现代类方法。这种方法可以在一定程度上消除频谱泄漏问题,但收敛速度和稳定性较差,应用也受到一定的限制,需要进一步研究和改进。在以上三种方法中,基于DFT的谐波检测方法是最成熟且常用的,尽管它在非同步采样条件下存在频谱泄露和栅栏现象。
现有技术的缺点如下:
1)基于DFT的传统方法,在检测谐波时不可避免地会遇到频谱泄漏和栅栏效应等问题。尽管通过Hanning窗,Nuttall窗等加窗插值算法能一定程度上改善频谱泄漏和栅栏效应,但这些窗难以同时满足主瓣能量集中与旁瓣衰减快等特点,并且非稳态下的加窗插值算法仍然存在修正公式求解复杂,计算精度难以提高等缺点。
2)现代谱估计方法,在检测谐波时需要一定的先验信息,其获得精确谐波分析结果的前提是已知信号中所含谐波分量的个数,并且这些方法对噪声敏感,计算量大。因此,它的应用非常有限。
3)现代类方法,在谐波分析时的收敛速度和稳定性较差,应用也受到一定的限制,需要进一步研究和改进。
发明内容
针对上述问题,本发明的目的在于提供一种基于时域Hermite插值的准同步谐波分析装置及方法,通过构建准同步采样序列以获得全周期截断信号解决频谱泄露和栅栏现象,有效提高了谐波分析的准确性。技术方案如下:
一种基于时域Hermite插值的准同步谐波分析方法,包括以下步骤:
步骤1:利用电力系统采样信号构造轨迹矩阵并进行奇异值分解;
步骤2:将奇异值分奇偶位置进行重新排序,计算奇异值均值序列,并将奇异值均值序列中各元素及其序号归一化;再计算差值曲线,确定有效阶次,并根据有效阶次重构信号,完成去噪处理,得到去噪信号;
步骤3:基于变分模态分解提取去噪信号中的基波分量;
步骤4:利用穿越时间后的采样点和穿越时间前的采样点,可获得估计的穿越时间,从而得到基波周期;
步骤5:基于数值微分计算各个采样点的导数;
步骤6:分别计算各个子区间上的三次Hermite插值多项式,并计算准同步采样时刻,将准同步采样时刻带入对应的三次Hermite插值多项式获得准同步采样数据;
步骤7:对准同步采样数据应用DFT进行谐波分析。
进一步的,所述步骤1中采样信号为电压信号时,表示为:
Uo(kTs)=Up(kTs)+N(kTs) (1)
式中,Up(kTs)表示去噪信号序列;N(kTs)表示噪声序列;k为序列索引,Ts为采样周期;
考虑一个长度为N的原始电压采样序列Uo=[uo1 uo2…uoN],将其分解为K=N-L+1个长度为L的向量,表示为:
hi=[uoi uo(i+1) … uo(i+L-1)]T,1≤i≤K (2)
将这些向量组成轨迹矩阵H:
令矩阵S=HHT,计算矩阵S的L个特征值λ1≥λ2≥…≥λL≥0,以及其所对应的标准正交向量U1,U2,…,UL;
令d=rank(H),以及则H表示为:
H=H1+H2+…+Hd (4)
式中,1≤l≤d,/>是轨迹矩阵H的奇异值,Ul和Vl是对应于/>的左特征向量和右特征向量;d为轨迹矩阵H的秩。
更进一步的,所述步骤2具体为:
步骤2.1:将奇异值分奇偶位置进行重新排序,并舍弃掉奇数序列最后一个奇异值:
式中,σodd表示奇数序列,σeven表示偶数序列;
步骤2.2:计算奇偶序列的均值,得到奇异值均值序列:
令x为奇异值均值序列σm中各元素的序号,即x=[1 2 … end];然后将x和σm归一化:
式中,和xai分别为第ai个归一化的序号和第ai个未归一化的序号;/>和/>分别为第ai个归一化的奇异值均值和第ai个未归一化的奇异值均值;
步骤2.3:将归一化曲线的“肘部”转化为“膝部”
式中,为转换后的第ai个归一化的奇异值均值,{σn}为归一化的奇异值均值序列;
曲线上离σs=xn最远的点即为曲线的膝点,根据相似三角形的原理,将曲线上的点到σs=xn的距离转换为该点横纵坐标的差值;则计算差值曲线:
式中,为第ai个差值;
步骤2.4:假设差值曲线中的最大值是{σdi}中的第M个点,则有效阶次r=2×(M-1);
步骤2.5:取式(4)中的前r个分量相加,得到矩阵HI1;
假设HI1中的元素为hI1(ki,kj),令L*=min(L,K),K*=max(L,K),N=L+K-1,则按式(10)重构得到去噪信号Up(kTs):
式中,L*为向量长度L和向量数量K中较小的数值;K*为向量长度L和向量数量K中较大的数值;upk为为去噪信号序列Up(kTs)中的第k个元素。
更进一步的,所述步骤3具体为:
步骤3.1:将原始的去噪信号Up(kTs)分解为z个有限带宽的模态函数,并采用希尔伯特变换对每个模态函数uz(t)进行信号解析,得到单边频谱;
步骤3.2:在每个模态函数uz(t)中加入指数项对中心频率ωz进行修正,并将每个模态的频谱调制至基频带宽;
步骤3.3:通过计算解调信号的梯度二范数估计每个模态函数uz(t)的基频带宽,构建带约束条件的变分模型:
式中,{uz}为模态集合;{ωz}为中心频率集合;为偏导运算符;δ(t)为单位脉冲函数;Up为去噪信号;
步骤3.4:将问题转化为非约束条件的变分模型,简化计算,引入Lagrange算子λ(t)和二次惩罚项α;得如下表达式:
式中,Up(t)为输入的去噪信号;uz(t)为模态函数;
步骤3.5:采用交替方向乘子法求解式(12),更新直至获得变分模型的最优解;
式中,为当前剩余量/>的维纳滤波;/>为Up(t)的傅里叶变换;/>为ui(t)的傅里叶变换;/>为λ(t)的傅里叶变换;/>为当前模态函数功率谱的重心;
步骤3.6:对最终获得的进行傅里叶逆变换得到分解出的z个分量,其中频率最接近50Hz的即为基波分量Uf(kTs)。
更进一步的,所述步骤4具体为:
利用穿越时间后的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)),通过下式获得估计的穿越时间,从而得到基波周期:
式中,Z为穿越时间所对应的去噪信号的值;f(·)表示相应采样点的均差;
将通过式(15)获得的时间序列表示为Tc=[tc1 tc2 … tc(end)],因此,第aj个基波周期的时间则表示为:
tpi=tc(i+1)-tci (16)
式中,和/>分别为第bj+1个和第bj个穿越时间;
则每个基波周期的持续时间构成的序列为Tp=[tp1 tp2 … tp(end)],据此确定Hermite插值重采样时刻。
更进一步的,所述步骤5具体为:
步骤5.1:对于由N个数据构成的去噪信号Up(kTs),将其对应的采样时刻序列TN=[Ts 2Ts…NTs]看作将区间[Ts NTs]等分成N-1份;由于:
式中,D(t)代表Up(t)的导数;
步骤5.2:对式(17)的右边利用Simpson公式,得到:
步骤5.3:对于式(18),记mk为D(kTs)的近似值,如果已知边界条件m1=D(Ts)和mN=D(NTs),则有:
对于边界条件m1=D(Ts)和mN=D(NTs),通过非中点的三点公式计算;
步骤5.4:通过点(Ts,up1),(2Ts,up2),(3Ts,up3)的二次插值多项式表示为:
式中,up1、up2、up3分别表示去噪信号Up(kTs)中的第1个、第2个和第3个元素;
则(Ts,up1)处导数值为:
将上式化简得:
式中,δk为kTs处的差商,由下式计算:
同理:
将(22)和(24)代入(19)即求得所有采样点的导数。
更进一步的,所述步骤6中,每个子区间kTs≤t≤(k+1)Ts上的插值多项式pk(t)通过下式计算:
令则第mi个准同步采样时刻通过下式计算:
式中,Nr为期望的每周期采样点数;为第q个基波周期的时间;/>为第q个穿越时间。
一种基于时域Hermite插值的准同步谐波分析装置,包括:
轨迹矩阵构造模块:利用电力系统采样信号构造轨迹矩阵并进行奇异值分解;
有效阶次重构模块:将奇异值分奇偶位置进行重新排序,计算奇异值均值序列,并将奇异值均值序列中各元素及其序号归一化;再计算差值曲线,确定有效阶次,并根据有效阶次重构信号,完成去噪处理,得到去噪信号;
基波分量提取模块:基于变分模态分解提取去噪信号中的基波分量;
基波周期计算模块:利用穿越时间后的采样点和穿越时间前的采样点,可获得估计的穿越时间,从而得到基波周期;
采样点导数计算模块:基于数值微分计算各个采样点的导数;
准同步采样数据计算模块:分别计算各个子区间上的三次Hermite插值多项式,并计算准同步采样时刻,将准同步采样时刻带入对应的三次Hermite插值多项式获得准同步采样数据;
谐波分析模块:对准同步采样数据应用DFT进行谐波分析。
本发明的有益效果是:本发明提出了一种基于时域Hermite插值的准同步谐波分析算法,该方法以电力系统采样数据为基础,将改进的奇异谱分析方法应用于采用数据进行去噪,通过改进奇异谱分析方法提高了去噪效果,增强了谐波分析的抗噪性;采用VMD和牛顿插值获得准确的基波周期,通过数值微分方法获得采样点对应的导数,从而获得分段三次Hermite插值多项式,最后带入准同步采样时刻获得准同步采样数据,达到消除DFT运算时的频谱泄露和栅栏现象的目的。并且,该方法以电力系统采样数据为基础,并不需要先验信息。
附图说明
图1为本发明基于时域Hermite插值的准同步谐波分析方法的流程图。
具体实施方式
下面结合附图和具体实施例对本发明做进一步详细说明。
电力系统的谐波畸变会对系统中的设备造成严重的危害,对电力系统的安全稳定构成了极大的威胁,为了准确检测电力系统信号中的谐波,本发明提出一种基于时域Hermite插值的准同步谐波分析方法,基本流程图如图1所示,具体步骤如下:
S1:利用电力系统采样信号构造轨迹矩阵H并进行奇异值分解。
电力系统中以等间隔时间采样的离散时间样本对应一个一维时间序列。以电压信号为例,可以表示为:
Uo(kTs)=Up(kTs)+N(kTs) (1)
考虑一个长度为N的原始电压采样序列Uo=[uo1 uo2…uoN],将其分解为K=N-L+1个长度为L的向量,L的选择一般为N/2,如果N为奇数的话,L的值可取(N+1)/2。
hi=[uoi uo(i+1) … uo(i+L-1)]T,1≤i≤K (2)
将这些向量组成轨迹矩阵H:
令S=HHT,计算S的L个特征值λ1≥λ2≥…≥λL≥0和其所对应的标准正交向量U1,U2,…,UL.
令d=rank(H)以及则H可以表示为:
H=H1+H2+…+Hd (4)
式中,1≤l≤d,/>是轨迹矩阵H的奇异值,Ul和Vl是对应于/>的左特征向量和右特征向量;d为轨迹矩阵H的秩。
S2:确定有效阶次r,并根据r重构信号,完成去噪处理。
首先,将奇异值分奇偶位置进行重新排序。如果奇异值数目为奇数,重新排序后会导致奇数序列比偶数序列多一个,此时为了方便后面计算,可舍弃掉最后一个奇异值:
计算奇偶序列的均值,得到奇异值均值序列:
令x为奇异值均值序列σm中各元素的序号,即x=[1 2 … end]。然后将x和σm归一化:
式中,和xai分别为第ai个归一化的序号和第ai个未归一化的序号;/>和/>分别为第ai个归一化的奇异值均值和第ai个未归一化的奇异值均值。
接着,将归一化曲线的“肘部”转化为“膝部”.
式中,为转换后的第ai个归一化的奇异值均值,{σn}为归一化的奇异值均值序列、
由于归一化曲线一定经过(0,0)和(1,1),因此曲线上离σs=xn最远的点即为曲线的膝点。根据相似三角形的原理,曲线上的点到σs=xn的距离可以转换为该点横纵坐标的差值。因此,接着计算差值曲线:
式中,为第ai个差值。
假设差值曲线中的最大值是中的第M个点,则有效阶次r=2×(M-1)。
取式(4)中的前r个分量相加,得到矩阵HI1。假设HI1中的元素为hI1(ki,kj),令L*=min(L,K),K*=max(L,K),N=L+K-1,则可按式(10)重构得到去噪信号Up(kTs)
S3:基于VMD提取去噪信号Up中的基波分量Uf。
VMD(Variational mode decomposition变分模态分解)是一种信号分解方法,该方法在获取分解分量的过程中通过迭代搜寻变分模型最优解来确定每个分量的频率中心和带宽,从而能够自适应地实现信号各分量的有效分离。该方法的具体过程如下:
(1)将原始的去噪信号Up(kTs)分解为z个有限带宽模态函数,并采用希尔伯特变换对每个模态函数uz(t)进行信号解析,得到单边频谱。
(2)每个模态函数uz(t)中加入指数项对中心频率ωz进行修正,并将每个模态的频谱调制至基频带宽。
(3)通过计算解调信号的梯度二范数估计每个模态函数uz(t)的基频带宽,构建带约束条件的变分模型如式(11)。
式中,{uz}为模态集合;{ωz}为中心频率集合;为偏导运算符;δ(t)为单位脉冲函数;Up为去噪信号
(4)为了将问题转化为非约束条件的变分模型,进而简化计算,引入Lagrange算子λ(t)和二次惩罚项α;可得如下表达式:
式中,Up(t)为输入的去噪信号;uz(t)为模态函数
(5)采用交替方向乘子法求解式(12),更新直至获得变分模型的最优解。
式中,为当前剩余量/>的维纳滤波;/>为Up(t)的傅里叶变换;/>为ui(t)的傅里叶变换;/>为λ(t)的傅里叶变换;/>为当前模态函数功率谱的重心。
对最终获得的进行傅里叶逆变换即可得到分解出的z个分量,其中频率最接近50Hz的就是基波分量Uf(kTs)。
S4:确定基波周期。
对于周期信号,如果满足一定条件的点在一个周期内只出现一次,则相邻点之间的时间间隔可视为信号的周期。例如正弦曲线上与y=Z相交且斜率为负的点一周期只出现一次。如果能准确的找到这些点,就能确定准确的基波周期。但是由于信号采样的限制,很难直接确定曲线穿过阈值Z的时间,但是该时刻前后的采样点我们可以很容易的获得。利用这些采样点,使用插值算法构造这一小段时间内信号的解析表达式,即可获得近似的穿越时间。
在估计基波周期过程中,牛顿四次插值多项式的误差很小,可以忽略不计。因此,利用穿越时间后的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)),通过下式即可获得估计的穿越时间,从而得到基波周期。
通过式(15)获得的时间序列可以表示为Tc=[tc1 tc2 … tc(end)],因此,第aj个基波周期的时间则可表示为:
所以,Tp=[tp1 tp2…tp(end)]是每个基波周期的持续时间构成的序列,据此则可确定Hermite插值重采样时刻。
S5:基于数值微分计算各个采样点的导数。
对于由N个数据构成的去噪信号Up(kTs)而言,其对应的采样时刻序列TN=[Ts2Ts…NTs],可以看成将区间[Ts NTs]等分成N-1份。由于:
式中,D(t)代表Up(t)的导数。
对式(17)的右边利用Simpson公式,可得:
对于式(18),记mk为D(kTs)的近似值,如果已知边界条件m1=D(Ts)和mN=D(NTs),则有:
对于边界条件m1=D(Ts)和mN=D(NTs),可通过非中点的三点公式计算。
通过点(Ts,up1),(2Ts,up2),(3Ts,up3)的二次插值多项式可表示为:
则(Ts,up1)处导数值为:
将上式化简可得:
式中,δk为kTs处的差商,由下式计算:
同理:
将(22)和(24)代入(19)即可求得所有采样点的导数。
S6:基于三次Hermite插值多项式,计算准同步采样数据。
计算出各点的导数后,就可以分别计算各个子区间上的三次Hermite插值多项式了。每个子区间kTs≤t≤(k+1)Ts上的插值多项式pk(t)可通过下式计算:
令则第mi个准同步采样时刻通过下式计算:
式中,式中,Nr为期望的每周期采样点数;为第q个基波周期的时间;/>为第q个穿越时间。
将准同步采样时刻带入对应的三次Hermite插值多项式即可获得准同步采样数据。
S7:对准同步采样数据应用DFT进行谐波分析。
综上,通过改进奇异谱分析方法提高了去噪效果,增强了谐波分析的抗噪性,最终得到的准同步采样数据克服了进行DFT运算时的频谱泄露和栅栏现象。
Claims (4)
1.一种基于时域Hermite插值的准同步谐波分析方法,其特征在于,包括以下步骤:
步骤1:利用电力系统采样信号构造轨迹矩阵并进行奇异值分解;
步骤2:将奇异值分奇偶位置进行重新排序,计算奇异值均值序列,并将奇异值均值序列中各元素及其序号归一化;再计算差值曲线,确定有效阶次,并根据有效阶次重构信号,完成去噪处理,得到去噪信号;
步骤3:基于变分模态分解提取去噪信号中的基波分量;
步骤4:利用穿越时间后的采样点和穿越时间前的采样点,可获得估计的穿越时间,从而得到基波周期;
步骤5:基于数值微分计算各个采样点的导数;
步骤6:分别计算各个子区间上的三次Hermite插值多项式,并计算准同步采样时刻,将准同步采样时刻带入对应的三次Hermite插值多项式获得准同步采样数据;
步骤7:对准同步采样数据应用DFT进行谐波分析;
所述步骤2具体为:
步骤2.1:将奇异值分奇偶位置进行重新排序;当轨迹矩阵H的秩d为奇数时,舍弃掉奇数序列最后一个奇异值:
式中,σodd表示奇数序列,σeven表示偶数序列;当轨迹矩阵H的秩d为偶数时,end=d;
当轨迹矩阵H的秩d为奇数时,end=d-1;
步骤2.2:计算奇偶序列的均值,得到奇异值均值序列:
令x为奇异值均值序列σm中各元素的序号,即x=[1,2…end];然后将x和σm归一化:
式中,和xai分别为第ai个归一化的序号和第ai个未归一化的序号;/>和/>分别为第ai个归一化的奇异值均值和第ai个未归一化的奇异值均值;
步骤2.3:将归一化曲线的“肘部”转化为“膝部”
式中,为转换后的第ai个归一化的奇异值均值,{σn}为归一化的奇异值均值序列;曲线上离σs=xn最远的点即为曲线的膝点,根据相似三角形的原理,将曲线上的点到σs=xn的距离转换为该点横纵坐标的差值;则计算差值曲线:
式中,为第ai个差值;
步骤2.4:假设差值曲线中的最大值是中的第M个点,则有效阶次r=2×(M-1);
步骤2.5:取式(4)中的前r个分量相加,得到矩阵HI1;
假设HI1中的元素为hI1(ki,kj),令L*=min(L,K),K*=max(L,K),N=L+K-1,则按式(10)重构得到去噪信号Up(kTs):
式中,L*为向量长度L和向量数量K中较小的数值;K*为向量长度L和向量数量K中较大的数值;upk为去噪信号序列Up(kTs)中的第k个元素;
所述步骤3具体为:
步骤3.1:将原始的去噪信号Up(kTs)分解为z个有限带宽的模态函数,并采用希尔伯特变换对每个模态函数uz(t)进行信号解析,得到单边频谱;
步骤3.2:在每个模态函数uz(t)中加入指数项对中心频率ωz进行修正,并将每个模态的频谱调制至基频带宽;
步骤3.3:通过计算解调信号的梯度二范数估计每个模态函数uz(t)的基频带宽,构建带约束条件的变分模型:
式中,{uz}为模态集合;{ωz}为中心频率集合;为偏导运算符;δ(t)为单位脉冲函数;Up为去噪信号;
步骤3.4:将问题转化为非约束条件的变分模型,简化计算,引入Lagrange算子λ(t)和二次惩罚项α;得如下表达式:
式中,Up(t)为输入的去噪信号;uz(t)为模态函数;
步骤3.5:采用交替方向乘子法求解式(12),更新直至获得变分模型的最优解;
式中,为当前剩余量/>的维纳滤波;/>为Up(t)的傅里叶变换;为ui(t)的傅里叶变换;/>为λ(t)的傅里叶变换;/>为当前模态函数功率谱的重心;
步骤3.6:对最终获得的进行傅里叶逆变换得到分解出的z个分量,其中频率最接近50Hz的即为基波分量Uf(kTs);
所述步骤4具体为:
利用穿越时间后的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)),通过下式获得估计的穿越时间,从而得到基波周期:
式中,Z为穿越时间所对应的去噪信号的值;f(·)表示相应采样点的均差;
将通过式(15)获得的时间序列表示为Tc=[tc1 tc2 … tc(end)],因此,第aj个基波周期的时间则表示为:
式中,和/>分别为第bj+1个和第bj个穿越时间;
则每个基波周期的持续时间构成的序列为Tp=[tp1 tp2 … tp(end)],据此确定Hermite插值重采样时刻;
所述步骤5具体为:
步骤5.1:对于由N个数据构成的去噪信号Up(kTs),将其对应的采样时刻序列TN=[Ts2Ts…NTs]看作将区间[Ts NTs]等分成N-1份;由于:
式中,D(t)代表Up(t)的导数;
步骤5.2:对式(17)的右边利用Simpson公式,得到:
步骤5.3:对于式(18),记mk为D(kTs)的近似值,如果已知边界条件m1=D(Ts)和mN=D(NTs),则有:
对于边界条件m1=D(Ts)和mN=D(NTs),通过非中点的三点公式计算;
步骤5.4:通过点(Ts,up1),(2Ts,up2),(3Ts,up3)的二次插值多项式表示为:
式中,up1、up2、up3分别表示去噪信号Up(kTs)中的第1个、第2个和第3个元素;
则(Ts,up1)处导数值为:
将上式化简得:
式中,δk为kTs处的差商,由下式计算:
同理:
将(22)和(24)代入(19)即求得所有采样点的导数。
2.根据权利要求1所述的基于时域Hermite插值的准同步谐波分析方法,其特征在于,所述步骤1中采样信号为电压信号时,表示为:
Uo(kTs)=Up(kTs)+N(kTs) (1)
式中,Up(kTs)表示去噪信号序列;N(kTs)表示噪声序列;k为序列索引,Ts为采样周期;考虑一个长度为N的原始电压采样序列Uo=[uo1 uo2 … uoN],将其分解为K=N-L+1个长度为L的向量,表示为:
hi=[uoi uo(i+1) … uo(i+L-1)]T,1≤i≤K (2)
将这些向量组成轨迹矩阵H:
令矩阵S=HHT,计算矩阵S的L个特征值λ1≥λ2≥…≥λL≥0,以及其所对应的标准正交向量U1,U2,…,UL;
令d=rank(H),以及则H表示为:
H=H1+H2+…+Hd (4)
式中,1≤l≤d,/>是轨迹矩阵H的奇异值,Ul和Vl是对应于/>的左特征向量和右特征向量;d为轨迹矩阵H的秩。
3.根据权利要求1所述的基于时域Hermite插值的准同步谐波分析方法,其特征在于,所述步骤6中,每个子区间kTs≤t≤(k+1)Ts上的插值多项式pk(t)通过下式计算:
令则第mi个准同步采样时刻通过下式计算:
式中,Nr为期望的每周期采样点数;为第q个基波周期的时间;/>为第q个穿越时间。
4.一种采用权利要求1所述基于时域Hermite插值的准同步谐波分析方法的准同步谐波分析装置,其特征在于,包括:
轨迹矩阵构造模块:利用电力系统采样信号构造轨迹矩阵并进行奇异值分解;
有效阶次重构模块:将奇异值分奇偶位置进行重新排序,计算奇异值均值序列,并将奇异值均值序列中各元素及其序号归一化;再计算差值曲线,确定有效阶次,并根据有效阶次重构信号,完成去噪处理,得到去噪信号;
基波分量提取模块:基于变分模态分解提取去噪信号中的基波分量;
基波周期计算模块:利用穿越时间后的采样点和穿越时间前的采样点,可获得估计的穿越时间,从而得到基波周期;
采样点导数计算模块:基于数值微分计算各个采样点的导数;
准同步采样数据计算模块:分别计算各个子区间上的三次Hermite插值多项式,并计算准同步采样时刻,将准同步采样时刻带入对应的三次Hermite插值多项式获得准同步采样数据;
谐波分析模块:对准同步采样数据应用DFT进行谐波分析。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202211253712.5A CN115508618B (zh) | 2022-10-13 | 2022-10-13 | 一种基于时域Hermite插值的准同步谐波分析装置及方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202211253712.5A CN115508618B (zh) | 2022-10-13 | 2022-10-13 | 一种基于时域Hermite插值的准同步谐波分析装置及方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN115508618A CN115508618A (zh) | 2022-12-23 |
CN115508618B true CN115508618B (zh) | 2023-10-03 |
Family
ID=84510578
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202211253712.5A Active CN115508618B (zh) | 2022-10-13 | 2022-10-13 | 一种基于时域Hermite插值的准同步谐波分析装置及方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN115508618B (zh) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116501131B (zh) * | 2023-06-27 | 2023-12-05 | 北京智芯微电子科技有限公司 | 频率偏差的确定方法、实时时钟的补偿方法及其系统 |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO1999046731A1 (en) * | 1998-03-13 | 1999-09-16 | The University Of Houston System | Methods for performing daf data filtering and padding |
CN105606892A (zh) * | 2015-12-17 | 2016-05-25 | 西安理工大学 | 一种基于sst变换的电网谐波与间谐波分析方法 |
CN108061821A (zh) * | 2016-11-05 | 2018-05-22 | 南京理工大学 | 一种改进的双馈风力发电系统谐波检测方法 |
CN109145476A (zh) * | 2018-08-31 | 2019-01-04 | 重庆水利电力职业技术学院 | 用于电力系统信号处理的时域自适应分段复指数级数法 |
CN113032716A (zh) * | 2019-12-24 | 2021-06-25 | 南京理工大学 | 基于加窗插值和Prony算法的谐波与间谐波分析方法 |
CN113702861A (zh) * | 2021-08-27 | 2021-11-26 | 广东省科学院电子电器研究所 | 一种基于大数据分析的电源故障位置预测方法及装置 |
Family Cites Families (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2004034231A2 (en) * | 2002-10-11 | 2004-04-22 | Flint Hills Scientific, L.L.C. | Intrinsic timescale decomposition, filtering, and automated analysis of signals of arbitrary origin or timescale |
TWI376507B (en) * | 2008-06-06 | 2012-11-11 | Univ Ishou | Method of optimizing spectrum of power analyzer |
-
2022
- 2022-10-13 CN CN202211253712.5A patent/CN115508618B/zh active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO1999046731A1 (en) * | 1998-03-13 | 1999-09-16 | The University Of Houston System | Methods for performing daf data filtering and padding |
CN105606892A (zh) * | 2015-12-17 | 2016-05-25 | 西安理工大学 | 一种基于sst变换的电网谐波与间谐波分析方法 |
CN108061821A (zh) * | 2016-11-05 | 2018-05-22 | 南京理工大学 | 一种改进的双馈风力发电系统谐波检测方法 |
CN109145476A (zh) * | 2018-08-31 | 2019-01-04 | 重庆水利电力职业技术学院 | 用于电力系统信号处理的时域自适应分段复指数级数法 |
CN113032716A (zh) * | 2019-12-24 | 2021-06-25 | 南京理工大学 | 基于加窗插值和Prony算法的谐波与间谐波分析方法 |
CN113702861A (zh) * | 2021-08-27 | 2021-11-26 | 广东省科学院电子电器研究所 | 一种基于大数据分析的电源故障位置预测方法及装置 |
Non-Patent Citations (6)
Title |
---|
Harmonic and Inter-harmonic Signal Analysis Based on Hilbert-Huang Transform;Fu Wei.etc;2011 Second International Conference on Digital Manufacturing & Automation;全文 * |
基于Hermite插值同步化算法的电网谐波与无功测量方法;章梦哲;王军;赵书娟;许珉;;电力系统保护与控制(12);全文 * |
基于三次样条插值的算术傅里叶谐波分析方法;武婕;刘开培;乐健;陈宜皇;;电测与仪表(01);全文 * |
基于奇异值总体最小二乘法的间谐波估计算法;沈睿佼;杨洪耕;吴昊;;电网技术(23);全文 * |
慕昆;马小雨;何国锋.基于DFT及谐波群泄露能量最小的电力谐波/间谐波分析方法.现代电子技术.2017,(01),全文. * |
熊杰锋;王柏林;孙艳.电力系统间谐波和谐波分析的海宁窗插值算法.自动化仪表.2010,(04),全文. * |
Also Published As
Publication number | Publication date |
---|---|
CN115508618A (zh) | 2022-12-23 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Ferdi | Computation of fractional order derivative and integral via power series expansion and signal modelling | |
CN115508618B (zh) | 一种基于时域Hermite插值的准同步谐波分析装置及方法 | |
CN107192878A (zh) | 一种基于压缩感知的电力系统谐波检测方法及装置 | |
US4658367A (en) | Noise corrected pole and zero analyzer | |
CN103744106A (zh) | 一种基于高斯滤波成形多道脉冲幅度分析装置 | |
CN111046791A (zh) | 一种基于含可变因子广义s变换的电流信号滤波去噪方法 | |
WO2006117191A2 (en) | Use of digital filters in a detector and method for obtaining energy spectra of nuclear radiation | |
US4654808A (en) | Noise corrected pole and zero analyzer | |
Vese et al. | Reduced non-convex functional approximations for image restoration & segmentation | |
Stefan et al. | Sparsity enforcing edge detection method for blurred and noisy Fourier data | |
Beylkin et al. | Nonlinear inversion of a band-limited Fourier transform | |
Bhanu et al. | On the computation of the complex cepstrum | |
Flickner et al. | Periodic quasi-orthogonal spline bases and applications to least-squares curve fitting of digital images | |
CN113655534B (zh) | 基于多线性奇异值张量分解核磁共振fid信号噪声抑制方法 | |
JP2599670B2 (ja) | 伝達関数の極・ゼロ点決定方法 | |
Mouffak et al. | Noncausal forward/backward two-pass IIR digital filters in real time | |
US20140174946A1 (en) | System and method for estimating impedance in electrochemical impedance spectroscopy | |
CN113191317B (zh) | 一种基于极点构造低通滤波器的信号包络提取方法和装置 | |
US4654809A (en) | Noise corrected pole and zero analyzer | |
Demigny et al. | Fast recursive implementation of the Gaussian filter | |
Gueguen et al. | Factorial linear modelling, algorithms and applications | |
CN117713505B (zh) | 一种开关电源高边电流检测方法及系统 | |
Hanke | An ϵ-free a posteriori stopping rule for certain iterative regularization methods | |
Bittanti et al. | Poles identification of an analog filter for nuclear spectroscopy via subspace-based techniques | |
Mi et al. | An adaptive algorithm for identification of Hammerstein models in frequency domain |
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 |