CN109145476B - 用于电力系统信号处理的时域自适应分段复指数级数法 - Google Patents
用于电力系统信号处理的时域自适应分段复指数级数法 Download PDFInfo
- Publication number
- CN109145476B CN109145476B CN201811015179.2A CN201811015179A CN109145476B CN 109145476 B CN109145476 B CN 109145476B CN 201811015179 A CN201811015179 A CN 201811015179A CN 109145476 B CN109145476 B CN 109145476B
- Authority
- CN
- China
- Prior art keywords
- power system
- signal
- formula
- sub
- complex exponential
- 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
- 238000000034 method Methods 0.000 title claims abstract description 94
- 238000012545 processing Methods 0.000 title claims abstract description 27
- 238000005070 sampling Methods 0.000 claims abstract description 37
- 230000003044 adaptive effect Effects 0.000 claims description 25
- 239000011159 matrix material Substances 0.000 claims description 13
- 238000004364 calculation method Methods 0.000 claims description 12
- 238000000354 decomposition reaction Methods 0.000 claims description 8
- 239000013598 vector Substances 0.000 claims description 8
- 238000004422 calculation algorithm Methods 0.000 claims description 5
- 230000002950 deficient Effects 0.000 claims description 2
- 230000001575 pathological effect Effects 0.000 claims description 2
- 238000006467 substitution reaction Methods 0.000 claims description 2
- 230000017105 transposition Effects 0.000 claims description 2
- 230000011218 segmentation Effects 0.000 abstract description 23
- 238000003672 processing method Methods 0.000 abstract description 7
- 230000008569 process Effects 0.000 abstract description 3
- 238000004088 simulation Methods 0.000 description 29
- 230000000694 effects Effects 0.000 description 4
- 238000002474 experimental method Methods 0.000 description 4
- 241001123248 Arma Species 0.000 description 3
- 230000005540 biological transmission Effects 0.000 description 3
- 238000004445 quantitative analysis Methods 0.000 description 3
- 230000009471 action Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 238000011166 aliquoting Methods 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 238000003780 insertion Methods 0.000 description 1
- 230000037431 insertion Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 230000002688 persistence Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 239000013589 supplement Substances 0.000 description 1
- 230000002459 sustained effect Effects 0.000 description 1
- 230000001052 transient effect Effects 0.000 description 1
- 238000004804 winding Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- 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/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/06—Energy or water supply
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Business, Economics & Management (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Health & Medical Sciences (AREA)
- Economics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Human Resources & Organizations (AREA)
- Public Health (AREA)
- Primary Health Care (AREA)
- Strategic Management (AREA)
- Tourism & Hospitality (AREA)
- General Business, Economics & Management (AREA)
- General Health & Medical Sciences (AREA)
- Computing Systems (AREA)
- Water Supply & Treatment (AREA)
- Marketing (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Algebra (AREA)
- Computer Hardware Design (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Complex Calculations (AREA)
Abstract
本发明公开了用于电力系统信号处理的时域自适应分段复指数级数法,涉及电力系统和现代信号处理技术领域,解决了现有电力系统信号处理方法不能直接计算出电力系统信号分量的幅值、频率、相位和衰减因子4个特性参数及不能处理带不连续非平稳分量的电力系统信号的问题,其技术方案要点是:首先,采用复指数级数对电力系统信号的时域采样序列进行拟合;然后,以方均根误差最小为目标自适应调整所截取采样序列的长度,实现采样序列的分段;最后,计算各子段上的4个特性参数。本发明不仅具有能同时获取电力系统信号分量的4个特性参数及处理带不连续非平稳分量的电力系统信号的优点,还可用于电力系统信号重构,在电力系统领域具有广阔的应用前景。
Description
技术领域
本发明属于电力系统和现代信号处理技术领域,更具体地说,其涉及用于电力系统信号处理的时域自适应分段复指数级数法。
背景技术
电力系统的信号包括电压、电流、频率、有功功率、无功功率和视在功率等。电力系统信号处理是电力系统分析的基础,高质量的电力系统信号处理是实现电力系统安全、稳定、经济、优质运行的前提条件。电力系统信号处理的本质是用数学方法将电力系统信号所含的各种分量分离开来,并采用特定的数学工具予以描述,从而提取或揭示重要的信息用以指导决策,包括时域和频域两类方法。
目前广泛采用的电力系统信号处理方法有离散傅里叶变换(discrete Fouriertransform,DFT)、小波变换(wavelet transform,WT)和自回归滑动平均建模(autoregressive moving average,ARMA)等。DFT及其各种变化形式(如傅里叶级数和快速傅里叶变换)仅适于处理平稳信号,当采样序列的时间长度不能准确地为各信号分量周期的整数倍或当信号含有随时间衰减或增长的直流分量、周期分量及非平稳噪声分量时,计算结果将受到极大影响。WT是时频处理方法,适于处理非平稳信号,但小波基函数的选取带有经验性,且分量的时频表征方式不能直观反映电力系统的特性,实际应用不方便。ARMA模型适宜描述平稳随机过程,故对具有非平稳性的电力系统信号进行建模缺乏严格的理论支持,且ARMA模型仅能给出信号分量的频率和幅值,在应用中具有局限性。
综上,为满足电力系统各分支专业的应用需求,要求信号处理方法应能同时获取信号分量的幅值、频率、相位和衰减因子4个特性参数,且须在时域内进行表征。因此,如何设计出能同时获取信号分量的幅值、频率、相位和衰减因子4个特性参数并能处理带不连续非平稳分量的电力系统信号的信号处理方法是目前迫切需要解决的问题。
发明内容
针对现有信号处理方法存在的问题,本发明提供了用于电力系统信号处理的时域自适应分段复指数级数法,不仅具有能同时获取电力系统信号分量的幅值、频率、相位和衰减因子4个特性参数及处理带不连续非平稳分量的电力系统信号的优点,还可用于电力系统信号的重构,即根据需要对已有电力系统信号进行插值或延拓,在电力系统领域具有广阔的应用前景。
本发明的上述技术目的是通过以下技术方案实现的:首先,采用复指数级数对电力系统信号的时域采样序列进行拟合;然后,以方均根误差(root mean square error,RMSE)最小为目标自适应地调整所截取的采样序列的长度实现采样序列的分段;最后,计算给出各子段上电力系统信号分量的幅值、频率、相位和衰减因子4个特性参数。
进一步,所述用于电力系统信号处理的时域自适应分段复指数级数法包括:
对电力系统的时域连续信号x(t)按采样频率fs等距采样后,得到N个时域采样数据x(nTs),简记为时域采样序列xn,n=0,1,···,N–1;采用p个复指数组成级数对各xn进行线性表示,即拟合:
式中,0<p≤N/2;Ts=1/fs为采样间距/周期(单位s);j为虚数单位;Ai、θi、αi和fi分别为第i个复指数的幅值(单位视电力系统信号的具体物理意义而定)、初相角(单位°)、衰减因子(单位s–1)和频率(单位Hz);设:
bi为复留数,si为复频率。
问题归结为已知xn,求解未知参数bi和zi。方法是视xn为p阶齐次前向差分方程:
的解。该齐次差分方程的通解形式为xn=czn(c为待定常数),将其代入式(3)可得一特征方程:
约去常数c,两边同乘z–n,式(4)变为:
求解式(5)可得p个特征根zi,并得到式(3)的齐次解:
式中,ci为待定常数,由边界条件确定。令ci=bi,则式(6)就与式(1)相同。
综上,复指数级数法的基本步骤如下:先由式(3)解出ak,k=0,1,···,p–1;再由式(5)解出zi;最后由式(1)或式(6)解出bi。
由式(3)求解ak时,为使用全部N个采样数据以最大限度地利用xn所携带的信息,采用xn的后向排序序列进行后向预测,即将xn依次向后延时(平移)一个时间单位(序号)后代入式(3),得一方程组,其矩阵形式为:
简记为:
Xa=xA (8)
其中,X为(N–p)×p矩阵,a和xA均为(N–p)×1列向量。式(8)为超定方程,须采用最小二乘法(least squares,LS)求解,得到LS解为:
aLS=(XTX)–1XTxA (9)
式中,上标T表示矩阵转置运算,XTX常严重秩亏和病态,无法正常求逆,需用Moore-Penrose的广义逆(伪逆)替代计算(方法1),或依据LS与奇异值分解(singularvaluedecomposition,SVD)的关系进行数值估计(方法2)。方法1以Frobenius范数||(XTX)(XTX)+–I||F最小化为目标,其中I为p×p单位矩阵,(XTX)+为XTX的伪逆,用以替代(XTX)–1;方法2以2范数||Xa–xA||2最小化为目标,与LS的目标等价。方法2比方法1在理论上更严谨,而实验也表明方法2的精度比方法1高3~4个数量级,故本发明采用方法2计算aLS,步骤为:先求X的SVD,得到X=UΣVT,其中U和V分别为(N–p)×(N–p)和p×p正交矩阵,Σ为(N–p)×p矩阵(其对角元素为X的p个奇异值σi,其余元素为0);再将U和V改写为列向量的形式U=[u1,u2,···,uN–p]和V=[v1,v2,···,vp];最后由下式计算aLS:
将求得的ak代入式(5)求解zi。为克服传统迭代法求解的低效性及其无法求取重根或复根的缺点,本发明提出一快速求解方法。设p×p友矩阵:
则A的p阶特征多项式为:
式中,I为p×p单位矩阵;符号|·|表示行列式计算;z为A的特征值。将Bp按第1行展开,得:
Bp=zBp–1+(–1)1+p(–1)p–1a0=zBp–1+a0 (13)
用相同的方法求取Bp–1,并对式(13)进行递推可得:
Bp=zBp–1+a0=z(zBp–2+a1)+a0=z2(zBp–3+a2)+a1z+a0=···
=zp+ap–1zp–1+···+a2z2+a1z+a0 (14)
式(14)的系数取负后即为A第p列的元素,故可由式(14)的系数ak按式(11)写出A;令式(14)等于0,即得到A的特征方程,且该方程与式(5)相同,故式(5)的根就是A的特征值。因此,对A进行特征值分解(eigenvalue decomposition,EVD),所得特征值就是方程(5)的根zi。该方法经一次EVD即可求出方程(5)的所有根,故求解效率和速度比传统迭代法高很多。另外,SVD可由两步EVD来实现,故信号处理程序仅需编制一个EVD子程序,即可实现SVD和方程(5)的求解。
求得zi后,将其与xn一起代入式(1)得一方程组,其矩阵形式为:
简记为:
Zb=xB (16)
Z为Vandermonde矩阵。可用求解式(8)的方法求解式(16),故不再赘述。
最后,由求得的bi和zi根据式(2)求出各信号分量的4个特性参数:
Ai=|bi|,θi=argbi,αi=ln|zi|/Ts,fi=argzi/(2πTs) (17)
式中arg[·]表示求复数的幅角主值,有argbi、arg zi∈(–π,π]。求幅角主值时须判断幅角所处的象限,以argbi为例,设其实部和虚部为Re(bi)和Im(bi),有:
式中,arctan[·]∈(–π/2,π/2)。注意,按式(17)计算得到的是复指数级数的参数,并非电力系统常用的实数型正、余弦信号的参数,还需通过Euler公式将2个复指数相加合成为1个实余弦信号,故除衰减因子αi外,实余弦信号的幅值、相角和频率分别为2Ai、|θi|和|fi|(对于实余弦信号分量,θi和fi是正负对称出现的)。
经实验发现,当复指数的个数p取最大值N/2时,可完全保证复指数级数对采样序列xn的拟合性,不会出现算法失稳的问题。但考虑到N可能为奇数,应用中应取p=int_[N/2],符号int_表示向负向取整,比如int_[5.95]=5。然而,经实验还发现,p=int_[N/2]虽能保证拟合性,但可能出现过拟合现象,即拟合数据在采样序列的尾部出现轻微振荡,使拟合误差有所增加。而改取:
p=int_[N/2]–1 (19)
则可在保证拟合性的同时避免出现过拟合现象。
进一步,所述用于电力系统信号处理的时域自适应分段复指数级数法的自适应分段策略具体包括:
因为复指数是关于时间t的连续函数,故要求被拟合的电力系统信号在采样序列占据的整个时间段内各个分量持续存在且特性参数不会发生突变。但部分电力系统信号,比如由数字故障录波器记录的输电线路自动重合闸的断续信号,可能不满足该要求。因此,本发明提出自适应分段策略,用以自动确定不连续信号分量的起点和终点,在起点和终点间的各个子段上分别采用复指数级数法实现拟合,从而保证复指数级数法的适用性和精确性。
设ns和ne分别为子段起点和终点的序号,nΔ为搜索步长(子段每次增加的点数),ef为子段的RMSE,Lmin、NΔm和Em分别为预先指定的子段的最小长度、最大搜索步长和最大RMSE。RMSE的计算公式如下:
式中,以符号∧标识拟合值;子段复指数的个数按式(19)确定,ps=int_[(ne–ns+1)/2]–1;注意,n的起始值不是ns,而是与式(1)中相同的0,相当于将xn的0起点右移到了当前子段的起点ns,也就是将每个子段均视为一由0起始的新序列。
进一步,所述用于电力系统信号处理的时域自适应分段复指数级数法的算法流程如图1,具体步骤为:
开始时,取子段长度为Lmin,可能因数据点少使ef≥Em;此时,保持nΔ=1不变并执行分支④,避免nΔ>1使当前子段的ne跨越分界点延伸至下一子段,继续使ef≥Em,陷入死循环。在循环执行若干次分支④后,找到使ef首次小于Em的ne,接着在满足ef<Em的条件下,取nΔ=NΔm,经循环执行分支②搜索子段的最大长度(即终点)。当再次出现ef≥Em时,转而执行分支③,以确定出精确的分界点,方法是将nΔ从NΔm逐步减半,采用图中的函数int45[nΔ/2]实现(符号int45表示4舍5入取整,比如int45[3.4]=3、int45[3.5]=4)。当ef<Em再次满足时,程序又转回执行分支②,因此时nΔ已改变故不再置其为NΔm。分支②和③交替执行,直到nΔ被减为0,表明精确分界点已找到,则计算并存储当前子段的特性参数,完成后将nΔ重置为1并调整ns和ne经分支①进行下一子段的搜索。
附图说明
图1为本发明提供的用于电力系统信号处理的时域自适应分段复指数级数法的算法流程;
图2为本发明实施例1提供的电力系统仿真信号及其拟合信号的波形;
图3为本发明实施例1提供的5种信号插值方法的插值结果及理想插值结果;
图4为本发明实施例1提供的采用时域自适应分段复指数级数法对电力系统仿真信号同时进行延拓和插值后的结果;
图5为本发明实施例2提供的真实电力系统故障信号及其拟合信号的波形;
图6为本发明实施例2提供的采用时域自适应分段复指数级数法对真实电力系统信号同时进行延拓和插值后的结果。
具体实施方式
以下结合附图1~6和实施例,对本发明作进一步说明。实施例1将时域自适应分段复指数级数法应用于一个由仿真产生的已知各分量特性参数且带有不连续非平稳分量的电力系统信号,以验证其有效性;实施例2将时域自适应分段复指数级数法应用于一个各分量特性参数未知的真实电力系统信号,以验证其实用性。
实施例1:包括电力系统仿真信号的拟合及其特性参数求取、电力系统仿真信号的插值及延拓3个子例。
1、电力系统仿真信号的拟合及其特性参数求取。
电力系统仿真信号,以下简称仿真信号,其离散表达式为:
式中,N=600,Ts=1ms,fs=1kHz;各信号分量xi(nTs)的表达式为:
具体参数如表1,其包含4个间谐波、1个3次谐波和1个衰减直流分量。
表1 6个信号分量的参数
序号i | Ai | αi/s–1 | fi/Hz | θi/(°) |
1 | 3.0 | –2.0 | 12.0 | 30.0 |
2 | 2.0 | –1.0 | 10.0 | 20.0 |
3 | 1.0 | –0.5 | 5.0 | 10.0 |
4 | 0.7 | –3.0 | 8.0 | 15.0 |
5 | 3.5 | –5.0 | 250.0 | 160.0 |
6 | 1.5 | –1.5 | 0.0 | 0.0 |
为模拟突然接入的及非持续的信号分量,设置:
式中,序号乘以Ts即得到相应的时间值;x1(t)从0s接入,持续到0.122s消失;x5(t)在0~0.122s内不存在,从0.123s接入,持续到0.299s消失;x6(t)在0~0.299s内不存在,从0.3s接入,一直持续到仿真信号的终止时刻;其余3个信号分量均由0s接入,一直持续到仿真信号的终止时刻。6个信号分量按式(22)合成的仿真信号波形如图2的(a)。
首先,暂不加入自适应分段策略,直接将复指数级数法用于对仿真信号进行整体的拟合。在式(21)中取ne–ns=N–1,将计算得到的拟合波形叠加绘制在仿真信号上,如图2的(b)。可见,虽然拟合数据点在整体上很接近仿真信号,但在某些局部(如图中的虚线圈部分),可看到拟合数据点偏离了仿真信号,此时p=299,已达到复指数个数允许的最大值300,但ef≈1.731×103,表明拟合精度仍非常低;若减小p,则ef更大,甚至出现无法拟合的情况。此时计算出的特性参数误差过大,故不再列出。因此,如果不加入自适应分段策略,复指数级数法在处理不连续信号分量时会存在缺陷。
接下来,用时域自适应分段复指数级数法对仿真信号进行拟合,选取实验参数Lmin=20、NΔm=20、Em=1×10–6。同样,将由式(21)计算得到的各子段的拟合信号依次叠加绘制在仿真信号上,如图2的(c)。图中的虚竖线为各子段的分界线,可见,时域自适应分段复指数级数法将仿真信号自动分为了4个子段,拟合数据点全部都在仿真信号上,从视觉上已无法看出差别。为评判分段效果和拟合精度,在表2中列出了各子段的详细信息及特性参数的计算值。可见,时域自适应分段复指数级数法准确地寻找到了式(24)设定的序号122、299和300,且各子段的ef数量级均在10–20左右,说明拟合精度非常高。
最后,通过对比表2和表1,进一步确认计算得出的特性参数的精度。先将表2子段1的4组特性参数与表1的前4行参数对比(根据式(24)的设定子段1仅包含x1(t)~x4(t)),可见两者相互对应,仅排列顺序相反,表明时域自适应分段复指数级数法准确地找到了子段1应包含的4个信号分量,且定量分析显示,计算得出的特性参数与其真实值间的相对误差在10–7量级。再将表2子段2的4组特性参数与表1的2~5行参数对比(根据式(24)的设定子段2仅包含x2(t)~x5(t)),可见两者的αi、fi相差很小,而Ai、θi却相差甚远。
表2各子段的详细信息及特性参数计算值
注意,这并不是计算有误,而是由时间偏移导致的“正确结果”:子段2各信号分量的Ai、θi实际上是x2(t)~x5(t)的Ai、θi经0.122s后的结果,比如,x4(t)的幅值初值为0.7,初相角为15.0°,在αi、fi的作用下经0.122s后,分别变为0.7e–3.0×0.122和arg[15.0°+2π×8.0×0.122×180°/π],具体数值见表3子段2数据的第1行;同样处理后,将x3(t)、x2(t)和x5(t)的Ai、θi经0.122s后的数值列于表3子段2数据的2~4行,这时再次将表3子段2的数据与表1的2~5行参数对比,可见两者的数值一致了,定量分析显示两者的相对误差也在10–7量级。后续的第3、第4子段也有类似情况,可仿照上述方法处理,其各信号分量经相应时间偏移后的Ai、θi的数值也列于表3中。总之,若已知一信号分量的接入时
表3考虑时间偏移量后各子段Ai、θi的计算值
间,欲使其在各子段上的分析结果一致,在计算Ai、θi时应计及相应的时间偏移量,这是在具体实施时要特别注意的。
2、电力系统仿真信号的插值。
信号插值是多采样率系统统一采样频率的主要手段,用于将低采样频率序列向高采样频率序列转换;此外,对一采样周期较大的采样序列,有时需获知其采样点间某些特定时刻的采样点,此时也需对原采样序列进行插值,以补出那些未被记录下来的采样点。
此处算例的目标是通过插值将仿真信号的采样频率由fs=1kHz提高为fs1=10kHz,采样周期由Ts=1ms缩小为Ts1=0.1ms。这需要在仿真信号的每2个点间插入9个点,即将1个Ts等分为10份,每1份为1个Ts1。因仿真信号各分量对fs、fs1均满足采样定理,故插值不会造成频谱混叠,无需增加滤波环节。
采用时域自适应分段复指数级数法实现信号插值的方法为:对于仿真信号的数据点,将其直接拷贝至相应的时间点上;而对于仿真信号数据点间的点,则由式(21)的如下变形形式产生:
上式的实质是在仿真信号的2个整数序号间插入9个浮点数序号,比如0,0.1,···,0.9,1和1,1.1,···,1.9,2等,而与这些浮点数序号对应的数据点和即为0~1和1~2间应插入的数据点。由式(1)可知,指数n+0.1k在计算时会与Ts相乘形成对应的时间(n+0.1k)Ts,这正是各整数序号所对应时间之间更细的时间,故即为这些更细时间上的数据。
以下将时域自适应分段复指数级数法的信号插值性能与多项式插值、填0滤波插值这两类目前使用最广的传统信号插值方法进行对比,以验证其有效性。图3依次示出了以上信号插值方法的插值结果和理想的插值结果,图中以实线表示fs时的仿真信号,以点线表示fs1时的插值信号。说明:①因理论证明高次插值不可取,且实验也表明4次及以上多项式的插值结果并不比3次多项式好,故图中仅给出3例低次多项式的插值结果;②默认以fs1时的仿真波形作为理想的插值结果即插入数据的真实值,由此计算各信号插值方法插值结果的RMSE(ef),ef也于各子图中示出;③因插值后数据总量增加10倍,插入数据点过于密集,从图上无法看出差异,故仅截取了0.115~0.135s的一段波形(因其包含0.122s这一分界点,很具代表性)进行横向展宽,以较清晰地揭示插值的细节。
图3的(a)为1次Lagrange多项式的插值结果。因两插值结点(仿真信号的数据点)间用直线连接,故插入数据点在仿真信号两点间的直线上排列齐整,看似很完美,然而ef却很大,对比图3的(f)可知,这是因为除插值结点外,多数插入的数据点均非fs1时波形上的数据点。
图3的(b)为3次Hermite多项式的插值结果。Hermite方法在Lagrange方法的基础上增加了插值结点处的导数约束,比Lagrange方法多使用了一些数据信息,故插值结果与理想插值结果稍接近了一些,ef也比图3的(a)小,但仍很大。
图3的(c)为3次样条函数的插值结果。样条插值的波形不仅光滑,且其ef较之前两种方法进一步减小,波形已比较接近图3的(f)的理想波形,但主要缺陷在于各子段的末端会因样条插值固有的末端效应而出现摆动(比如第1子段的末端)。
图3的(d)为填0滤波的插值结果。该方法是在拟插入数据的位置上填上0,用0替代其数值,然后经抗混叠滤波器滤除折叠频率以上因填0作用复制产生的频率分量,插值效果与滤波器密切相关。经实验发现,长度10、截止数字域频率0.5(对应模拟域频率0.5fs1/(2π)≈795.77Hz)的数字滤波器可使ef基本为最小,且使插值精度又较图3的(c)提高了一些。然而,这是在理想插值结果已知的情况下通过多次尝试设计的数字滤波器,在实际应用中不具通用性,难以保证数字滤波器参数对所有信号都是最优的;其次,数字滤波器也存在末端效应。
图3(e)为时域自适应分段复指数级数法的插值结果。其插值精度较之图3的(d)有大幅提高,ef的量级降为10–7。
小结:造成4种传统插值方法插值误差大的主要原因在于其均不具备信号信息的提取和利用能力;而时域自适应分段复指数级数法能充分提取和利用信号信息,自适应地、目标明确地实现插值。
3、电力系统仿真信号的延拓。
由一段已知信号推知其终止时刻后或起始时刻前的波形,称为信号的正向延拓(也称为预测)或负向延拓(也称为反演)。若电力系统信号经时域自适应分段复指数级数法实现了分段,则正向延拓波形与最后一个子段的相关性最强,负向延拓波形与第1子段的相关性最强,故分别按最后一个子段和第1个子段的复指数级数拟合式实现正向和负向延拓。因延拓和插值同时进行的实用价值更大,仅延拓、不插值的情况可视为其简化形式,故以仿真信号延拓、插值同时进行的情况为例。
此处算例的目标是对仿真信号分别进行正、负向延拓,各延拓0.1s时长,并同时进行插值。设NF和NB分别为正、负向延拓的点数,因0.1s/Ts1=1000,故正、负向需各延拓NF=1000点和NB=1000点。
先进行负向延拓。因仿真信号的起始时刻为0,延拓信号的时间值为负值,故各数据点的序号不仅为浮点数,且为负数。负向延拓的公式为:
式中,参数为计算得到的第1子段的参数;注意,计算式(26)前,应先提取第1子段的主导信号分量(ps取主导信号分量的个数),去除虚假分量,用主导分量的参数进行计算,因的指数为负数,相当于做除法,若某些虚假分量的/>过小,将使计算出的/>有很大误差。
接下来进行正向延拓。延拓数据点的序号为正浮点数,无负向延拓的问题,故无需提取主导分量,但为保证波形光滑,须从最后一个子段的最后一个数据点开始添加延拓数据。设最后一个数据点的序号为nlast(相当于最后一个子段的ns,由第4子段的信息知nlast=146),正向延拓的公式为:
经正、负向延拓并插值后的仿真信号如图4的点线。图中,实线为由式(22)~(24)按Ts1在–0.1~0.7s上计算得到的精确波形,可见,–0.1~0s的负向延拓波形和0.6~0.7s的正向延拓波形与精确波形重合得很好,定量分析表明其ef均在10–7量级。
实施例2:包括真实电力系统信号的拟合及其特性参数求取、插值和延拓2个子例。
1、真实电力系统信号的拟合及其特性参数求取。
图5的(a)为一由数字故障录波器记录的电力系统中某条输电线路发生短路故障前后的实际故障电流信号,为便于处理对数值进行了归一化处理即取标么值。
首先,暂不加入自适应分段策略,直接将复指数级数法用于对实际故障电流信号进行整体的拟合,结果如图5的(b)。可见,在0~0.25s拟合信号与实际故障电流信号重合较好,但到末端时拟合信号出现了发散,相应的ef很大,求取特性参数已无意义。这再次说明,如果不加入自适应分段策略,复指数级数法在处理含有不连续信号分量的电力系统信号时会存在缺陷。
接下来,用时域自适应分段复指数级数法对实际故障电流信号进行拟合,选取实验参数Lmin=20、NΔm=20、Em=5×10–4,适当加大Em是因为实际故障电流信号成分复杂,寻求过高的拟合精度是不现实的,若Em过小,则判据ef<Em可能永远无法得到满足。拟合结果如图5的(c),图中用虚竖线示出了各子段的分界线。因实际故障电流信号所含噪声分量较多,故仅在表4中列出各子段的ef及主导分量的特性参数。时域自适应分段复指数级数法将实际故障电流信号分为了6个子段:1~3子段的基波的幅值基本相同且较小,据此可判定此3子段为故障前的正常运行电流;到第4子段,基波的幅值迅速增大为前3段的约5倍,据此可判定此子段为短路暂态电流;到第5子段,基波的幅值较之第4子段减小很多,可判定继电保护在动作切除故障输电线路;到第6子段,基波缺失了,仅剩下1个伴随整个故障过程的3次谐波,这是由电力系统的非线性元件产生的1个环绕于电力变压器三角形绕组的固有环流。可见,正是因为实际故障电流信号中基波的非持续性导致了不加入自适应分段策略时算法的失效。注意,表中前5个子段的RMSE均满足ef<Em,而第6子段的ef>Em,这是因为其为最后一段,无再多的数据使ef继续下降,只能将其独立为一段,但其ef还是可接受的。本子例还可视为是时域自适应分段复指数级数法在电力系统故障诊断方面的一个应用案例。
表4故障信号各子段的ef及若干主导分量的特性参数
2、真实电力系统信号的插值和延拓。
本子例的目标是将实际故障电流信号的采样频率从原先的1kHz经插值后提高到10kHz,且在正、负向延拓各0.1s(即5个工频周波)。将式(26)和(27)用于实际故障电流信号,经插值、延拓后的实际故障电流信号如图6,注意,因实际故障电流信号在采样前会经过抗混叠滤波,而时域自适应分段复指数级数法不会无端地产生频率分量,故插值、延拓不会造成频谱混叠。图中,0s虚竖线左侧为负向延拓信号,0.3s虚竖线右侧为正向延拓信号。可见,正向延拓信号有逐渐增大且出现振荡的趋势,这表明与经典预测方法一样,随着预测时长的增加,预测的精度和可信度都将下降,因为当前的已知信息绝不能完全取代未来的信息;同样,对负向延拓(反演)也如此。不过,若对延拓时长加以限制,比如仅延拓0.1s左右,则精度和可信度均可得到保证。
以上所述仅为本发明的实施例,并不用以限制本发明,凡在本发明精神和原则内所作的任何修改、等同替换和改进等,均应包含在本发明的保护范围内。
Claims (5)
1.用于电力系统信号处理的时域自适应分段复指数级数法,其特征在于,采用复指数级数对电力系统信号的时域采样序列进行拟合:
对电力系统的时域连续信号x(t)按采样频率fs等距采样后,得到N个时域采样数据x(nTs),简记为时域采样序列xn,n=0,1,…,N–1;采用p个复指数组成级数对各xn进行线性表示,即拟合:
式中,0<p≤N/2;Ts=1/fs为采样间距/周期,单位为s;j为虚数单位;Ai为第i个复指数的幅值,其单位视电力系统信号的具体物理意义而定;θi为初相角,单位为°;αi为衰减因子,单位为s–1;fi为频率,单位为Hz;设:
已知xn求解未知参数bi和zi,设xn为p阶齐次前向差分方程的解:
该齐次前向差分方程的特征方程为:
复指数级数法的基本步骤如下:先由式(3)解出ak,k=0,1,…,p–1;再由式(4)解出zi;最后由式(1)解出bi;
由式(3)求解ak时,采用xn的后向排序序列进行后向预测,即将xn依次向后延时、平移一个时间单位、序号后代入式(3),得一方程组,其矩阵形式为:
简记为:
Xa=xA (6)
其中,X为(N–p)×p矩阵,a和xA均为(N–p)×1列向量;式(6)为超定方程,采用最小二乘法(least squares,LS)求解,得到LS解为:
aLS=(XTX)–1XTxA (7)
式中,上标T表示矩阵转置运算;XTX常严重秩亏和病态,无法正常求逆;依据LS与奇异值分解(singular value decomposition,SVD)的关系对aLS进行数值估计,该方法以2范数||Xa–xA||2最小化为目标,与LS的目标等价,步骤为:
先求X的SVD,得到X=UΣVT,其中U和V分别为(N–p)×(N–p)和p×p正交矩阵,Σ为(N–p)×p矩阵,其对角元素为X的p个奇异值σi,其余元素为0;再将U和V改写为列向量的形式U=[u1,u2,…,uN–p]和V=[v1,v2,…,vp];最后由下式计算aLS:
将求得的ak代入式(4)求解zi,所采用的快速求解方法为:首先,用ak组成p×p友矩阵:
然后,对A进行特征值分解(eigenvalue decomposition,EVD),所得特征值就是方程(4)的根zi;
求得zi后,将其与xn一起代入式(1)得一方程组,其矩阵形式为:
简记为:
Zb=xB (11)
用求解式(6)的方法求解式(11);
最后,由求得的bi和zi根据式(2)求出各信号分量的4个特性参数:
Ai=|bi|,θi=argbi,αi=ln|zi|/Ts,fi=arg zi/(2πTs) (12)
式中arg[·]表示求复数的幅角主值,有arg bi、arg zi∈(–π,π];求幅角主值时须判断幅角所处的象限,以arg bi为例,设其实部和虚部为Re(bi)和Im(bi),有:
式中,arctan[·]∈(–π/2,π/2);按式(12)计算得到的是复指数级数的参数,通过Euler公式将2个复指数相加合成为1个实余弦信号,除衰减因子αi外,实余弦信号的幅值、相角和频率分别为2Ai、|θi|和|fi|;
复指数的个数p按下式确定:
p=int–[N/2]–1 (14)
式中,符号int–表示向负向取整。
2.如权利要求1所述的用于电力系统信号处理的时域自适应分段复指数级数法,其特征在于,以方均根误差最小为目标自适应地调整所截取的采样序列的长度,实现采样序列的分段,并计算给出各子段上电力系统信号分量的幅值、频率、相位和衰减因子4个特性参数:
设ns和ne分别为子段起点和终点的序号,n△为搜索步长,即子段每次增加的点数,ef为子段的RMSE,Lmin、N△m和Em分别为预先指定的子段的最小长度、最大搜索步长和最大RMSE,RMSE的计算公式如下:
式中,以符号∧标识拟合值;子段复指数的个数按式(14)确定,ps=int–[(ne–ns+1)/2]–1;其中n的起始值不是ns,而是与式(1)中相同的0。
3.如权利要求2所述的用于电力系统信号处理的时域自适应分段复指数级数法,其特征在于,所述时域自适应分段复指数级数法的算法流程如下:
开始时,取子段长度为Lmin,当数据点少ef≥Em时,保持n△=1不变并执行分支④;
在循环执行若干次分支④后,找到使ef首次小于Em的ne,在满足ef<Em的条件下,取n△=N△m,经循环执行分支②搜索子段的最大长度,即子段的终点;
当再次出现ef≥Em时,转而执行分支③,采用函数int45[n△/2]将n△从N△m逐步减半,确定出精确的分界点,此处符号int45表示4舍5入取整;
当ef<Em再次满足时,又转回执行分支②,但n△不再置为N△m;
分支②和③交替执行,当n△被减为0,表明精确分界点已找到,则计算并存储当前子段的特性参数,完成后将n△重置为1并调整ns和ne经分支①进行下一子段的搜索。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811015179.2A CN109145476B (zh) | 2018-08-31 | 2018-08-31 | 用于电力系统信号处理的时域自适应分段复指数级数法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811015179.2A CN109145476B (zh) | 2018-08-31 | 2018-08-31 | 用于电力系统信号处理的时域自适应分段复指数级数法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109145476A CN109145476A (zh) | 2019-01-04 |
CN109145476B true CN109145476B (zh) | 2023-05-26 |
Family
ID=64826047
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201811015179.2A Active CN109145476B (zh) | 2018-08-31 | 2018-08-31 | 用于电力系统信号处理的时域自适应分段复指数级数法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109145476B (zh) |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110888365A (zh) * | 2019-11-30 | 2020-03-17 | 国网辽宁省电力有限公司锦州供电公司 | 电网系统非同步采样基波数据同步化方法 |
CN111310303B (zh) * | 2020-01-17 | 2024-02-02 | 合肥工业大学 | 一种幅值指数衰减的正弦波参数识别方法 |
CN112379216A (zh) * | 2020-11-05 | 2021-02-19 | 山东大学 | 输电线故障测距基波分量提取中量测误差抑制方法及系统 |
CN115508618B (zh) * | 2022-10-13 | 2023-10-03 | 四川大学 | 一种基于时域Hermite插值的准同步谐波分析装置及方法 |
Family Cites Families (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
GB2514595B (en) * | 2013-05-30 | 2017-10-18 | Imp Innovations Ltd | Method and apparatus for estimating frequency domain representation of signals |
CN103630742B (zh) * | 2013-12-16 | 2015-09-30 | 国家电网公司 | 一种动态信号参数的获取方法 |
CN108152651A (zh) * | 2017-12-27 | 2018-06-12 | 重庆水利电力职业技术学院 | 基于gmapm和som-lvq-ann的输电线路故障综合识别方法 |
-
2018
- 2018-08-31 CN CN201811015179.2A patent/CN109145476B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN109145476A (zh) | 2019-01-04 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109145476B (zh) | 用于电力系统信号处理的时域自适应分段复指数级数法 | |
Yu | A discrete Fourier transform-based adaptive mimic phasor estimator for distance relaying applications | |
Hassanpour et al. | Time domain signal enhancement based on an optimized singular vector denoising algorithm | |
Damera-Venkata et al. | Design of optimal minimum-phase digital FIR filters using discrete Hilbert transforms | |
Simon et al. | An efficient nonlinear least square multisine fitting algorithm | |
CN109301842B (zh) | 一种基于负阻尼贡献的风电场次同步振荡切除方法 | |
CN109684937B (zh) | 一种基于fft及数学形态法的信号去噪方法及装置 | |
Bi et al. | Sampling rate conversion in the frequency domain [dsp tips and tricks] | |
CN111948454B (zh) | 一种抗衰减直流分量干扰的同步相量测量方法及系统 | |
CN109783770A (zh) | 频率分析方法及装置 | |
CN111650422B (zh) | 高压直流系统的同步触发方法、系统及存储介质 | |
CN109358327A (zh) | 一种非线性调频信号的迭代生成方法、终端和存储介质 | |
CN113267676A (zh) | 格策尔算法下的频谱确定方法、系统、设备及存储介质 | |
Das et al. | Application of Wavelet transform in denoising synchronising signal in line synchronised power electronics converters | |
US20190360470A1 (en) | Estimating a characteristic of a wind turbine electrical signal | |
CN111835038B (zh) | 预同步方法、系统及终端设备 | |
CN114460363A (zh) | 交流信号同步采样方法、装置和交流采样设备 | |
CN110082642B (zh) | 基于全相位微分滤波的电网工况故障时刻检测方法及装置 | |
CN114184838A (zh) | 基于sn互卷积窗的电力系统谐波检测方法、系统及介质 | |
Karamta et al. | Implementation of extended Kalman filter based dynamic state estimation on SMIB system incorporating UPFC dynamics | |
CN111308170A (zh) | 电机相电流的采样方法及装置 | |
Treichler | γ-LMS and its use in a noise-compensating adaptive spectral analysis technique | |
CN108631315A (zh) | 基于泰勒级数展开的重复控制分数延迟滤波器设计方法 | |
CN117074778B (zh) | 基于负载场景的谐波提取方法、装置和计算机设备 | |
CN108767821B (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 |