CN101291055B - 一种输电线路故障行波初始波头到达时刻精确标定方法 - Google Patents

Abstract

本发明涉及一种利用希尔伯特-黄变换(HHT)的输电线路故障行波初始波头到达时刻精确标定方法。本方法为：对故障行波数据做HHT变换，首先进行EMD分解，其IMF₁分量集中体现了原始数据的高频信息，而且完全自适应于原始数据的时间尺度，以IMF₁分量为考察对象，求其一阶差分曲线并作Hilbert变换求得其瞬时频率曲线，结合两者的优势寻找极值点并与原始数据相比对获得行波波头到达时刻的准确标定，误差被控制在一个采样点之间。在强噪环境下，选取IMF₂或IMF₃，依上方法，同样可以达到行波波头到达时刻的准确标定。其他方数学方法均不能达到这样的精度。原理分析和仿真数据、工程数据验证表明，该方法对输电线路故障行波初始波头到达时刻标定准确。

Description

一种输电线路故障行波初始波头到达时刻精确标定方法

技术领域

本发明涉及电力系统继电保护技术领域，具体地说是一种利用希尔伯特-黄变换(HHT)的输电线路故障行波初始波头到达时刻精确标定方法。

背景技术

与传统的基于工频量的故障测距技术相比而言，基于行波的测距技术^[1-6]测距精度高、不易受系统运行方式、过渡电阻、CT饱和、线路分布电容的影响，因而成为继电保护领域热门的研究课题之一。其中行波波头到达时刻的精确定位是提高测距精度的关键技术之一。小波分析^[4，7，8]、数学形态学^[9-10]以及小波-数学形态学相结合^[12]的方法在行波波头到达时刻定位中的应用比较成熟。但是小波分析方法存在小波基的选择问题，在实际应用环境中选择怎样的小波基对分析结果影响很大，同时小波基的选择却还没有成熟的理论指导。数学形态学同样存在类似的问题，结构元素的选取没有理论指导。在行波波头检测的工程应用中，三次B样条小波被认为是具有良好信号奇异性检测能力的小波基，数学形态学也构造了独特的具有不同原点的扁平结构元素并结合多分辨形态梯度达到信号奇异性的良好检测效果，但是这些方法在反应波头位置时并不能精确到采样点，依然不能精确定位波头时刻。

HHT变换^[13-16]是一种独特的完全自适应的时频分析方法，它把EMD分解与经典的Hilbert变换相结合，直接从信号本身的时间尺度反应信号的时频本质。在信号奇异性检测方面具有独特的优势。本发明将HHT变换这一特性成功应用于故障行波波头到达时刻精确标定，仿真数据验证和实际工程数据验证表明了本方法的正确性、可行性。

参考文献：

[1]葛耀中，徐丙垠，陈平.利用暂态行波测距的研究[J].西安交大学报，1995，29(3)：70-75

[2]徐丙垠，葛耀中，朱锡贵.利用暂态电流行波的输电线路故障测距技术[C].第五届全国继电保护学术会议论文集：125-132

[3]徐丙垠.利用暂态行波的输电线路故障测距技术[D].西安交通大学博士学位论文，1991

[4]覃剑.小波变换应用于输电线路行波故障测距的研究[D].电力科学院，1998

[5]李友军，王俊生，郑玉平等.几种行波测距算法的比较.电力系统自动化[J]，2001，25(7)：36-39

[6]陈平，葛耀中，索南加乐等.基于故障开断暂态行波信息的输电线路故障测距研究[J].中国电机工程学报，2000，20(8)：56～60

[7]Mallat，S.and Hwang，W.L.Singularity detection and processing with wavelets[J].IEEE Trans on information Theory，Vol.38 1992.

[8][美]崔锦泰著，程正兴译.小波分析导论[M].西安交通大学出版社，1994.

[9]邹力，刘沛，赵青春.级联形态梯度变换及其在继电保护中的应用[J].中国电机工程学报，2004，24(12)：113-118

[10]林湘宁，刘沛，高艳.基于故障暂态和数学形态学的超高速线路方向保护[J].中国电机工程学报，2005，25(4)：13-18.

[11]吴青华，张东江.形态滤波技术及其在继电保护中的应用[J].电力系统自动化，2003，27(7)：45-49.

[12]林湘宁，刘沛，刘世明等.电力系统超高速保护的形态学-小波综合滤波算法[J].中国电机工程学报，2002，22(9)：19-24

[13]N.E.Huang，Z.Shen，S.R.Long，M.L.Wu，.H.Shih，Q.Zheng，N.C.Yen，C.C.Tungand H.H.Liu，″The empirical mode decomposition and Hilbert spectrum for nonlinear andnon-stationary time series analysis”[J]，Proc Roy Sec.London A，Vol.454，pp.903-995，1998.

[14]Norden E.Huang，Zheng Shen，Steven R.Long.″A new view of nonlinear waterwaves：The Hilbert Spectrum″[J]，Annu.Rev.FluidMech，vol.31，pp.417-457，1999.

[15]谭善文，多分辨希尔伯特一黄(Hilbert-Huang)变换方法的研究[D].重庆大学，2001.

[16]钟佑明.希尔伯特一黄变换局瞬信号分析理论的研究[D].重庆大学，2002.

发明内容

本发明的目的在于克服现有故障行波初始波头到达时刻判断技术的不足，提供一种利用.希尔伯特一黄变换(HHT)变换的输电线路故障行波初始波头到达时刻精确标定方法。它能发现EMD对信号时间尺度的自适应能力能够精确标定故障行波初始波头到达时刻，从而实现高精度双端故障行波测距。

本发明的理论基础如下：

1.希尔伯特(Hilbert)变换：

对于任一连续的时间信号x(t)，可得到它的Hilbert变换Y(t)为：

$Y\left(t\right)=\frac{1}{\mathrm{\π}}{\∫}_{-\∞}^{+\∞}\frac{X\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}\mathrm{d\τ}---\left(1\right)$

其反变换为：

$X\left(t\right)=\frac{1}{\mathrm{\π}}{\∫}_{-\∞}^{+\∞}\frac{Y\left(\mathrm{\τ}\right)}{\mathrm{\τ}-t}\mathrm{d\τ}---\left(2\right)$

得到解析信号：

Z(t)＝X(t)+iY(t)＝A(t)e^iθ(t)(3)

式中：A(t)为瞬时振幅；θ(t)为相位，其中：

A(t)＝[X(t)²+Y(t)²]^1/2(4)

$\mathrm{\θ}\left(t\right)=\arctan \frac{Y\left(t\right)}{X\left(t\right)}---\left(5\right)$

瞬时频率按下式计算：

$f\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{\mathrm{d\θ}\left(t\right)}{\mathrm{dt}}---\left(6\right)$

上式定义瞬时频率要求θ(t)必须满足一定的条件，即θ(t)必须是本征模态函数。

2.固有模态分量(IMF)：

Hilbert变换有一定的适用范围，只有当信号的表达式满足一定条件时才能讨论瞬时参数。要得到有实际物理意义的瞬时参数，所考察的信号应满足下面两个条件：

(1)在整个数据中，极值的个数与穿过零点的个数相等或至多相差1；

(2)信号上任意一点，由局部极大值点确定的包络线和局部极小值确定的包络线的均值为零，即信号关于时间轴局部对称。

满足上述条件的信号称为固有模态信号，相应的函数称为固有模态函数，并认为只有对IMF分量求出的瞬时参数才有实际的物理意义。

3.固有模态分解(EMD)过程：

对于一个时间序列x(t)，其经验模态分解过程如下：

(1).确定原始信号x(t)的所有极大值点和极小值点；

(2).采用样条函数求出x(t)的上、下包络线，并计算均值m(t)；

(3).作差h(t)＝x(t)-m(t)；

(4).将h(t)视为新的x(t)重复第3步，得到一个新的h(t)；

(5).用前后两个h(t)作标准差sd＝[h_k-1(t)-h_k(t)]²/h² _k-1(t)；

(6).判断sd是否小于0.25，满足则提取一个IMF分量imf(t)，否则重复第3步和第4步；

(7).作差r(t)＝x(t)-imf(t)；

(8).将r(t)视为新的x(t)重复第7步，得一个新的r(t)；

(9).作前后两个r(t)的标准差；

(10).判断上步的值是否小于0.25，满足则提取第二个IMF分量，否则重复第7步第8步；

(11).重复7到10步，求得所有IMF分量，不能提取的为残量；

$x\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{im}{f}_i\left(t\right)+r\left(t\right)---\left(10\right)$

式中r(t)为残余函数，它是一单调函数。

4.固有模态分解(EMD)中止条件：

EMD分解包含两个循环过程，求取IMF分量的过程称为筛选过程，那么在什么情况下停止IMF的筛选、以及什么时候中止分解成为EMD分解过程的两个关键问题。这两个终止条件分别被称为分量终止条件和分解终止条件。

(1).分量终止条件

筛选过程就是求取一个IMF的过程，其基本方法是从原始信号中不断找出极值。按照分解步骤不断进行筛选，直到满足一定的条件为止。这个筛选过程的目的就是不断的减小信号的不对称性，使波形趋向于关于零均值线对称。从而满足IMF分量的基本特征，能够通过Hilbert变换来计算瞬时频率。为了保证分解得到的IMF分量具有足够的调频调幅的物理意义，筛选的循环过程不能过多，太多的循环次数会使IMF分量过度平滑，失去原有意义；而过少的循环次数则使得得到的IMF分量不能满足IMF分量的基本特征，也就无法获取有意义的瞬时频率。因此对于分量终止条件，研究者提出了两种终止准则：类柯西收敛准则和简单终止准则。

类柯西收敛准则：

$\mathrm{SD}=\frac{\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{|h_{1,k-1}\left(t\right)-h_{1,k}\left(t\right)|}^2}{\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{\left|h_{1,k-1}\left(t\right)\right|}^2}---\left(14\right)$

介于0.2与0.3之间时，筛选过程终止。

事实上，EMD分解的目的在于分解出IMF分量，而IMF分量具备两个条件：极值点数和过零点数相等或者至多相差一个，局部对称。而在IMF分量的两个条件中局部对称是主要的，在一般情况下，只要信号满足局部对称性，根据Hilbert计算出来的瞬时频率也就具有了实际物理意义。因此，定义如下分量终止准则：

mean[h_i，k]≤ε(15)

其中mean[h_i，k(t)]≤ε表示h_i，k(t)的均值线，ε为一个预先设定好的足够小的数值。

此外，可以通过附加设置一个最大循环次数来终止当前的筛选过程，强制终止，提高算法的鲁棒性。

(2).分解终止条件

EMD分解得到的是IMF分量，所以分解的第二个循环必须有终止条件，决定EMD分解何时结束。按照前文的分解流程，当前一个IMF分量被提取出来后，总会留下一个新的残余量r_i，Huang等人认为当r_i成为一个常量，或者是一个单调函数时，已经不能再分解出IMF分量，分解过程终止。同时，如果信号本身具有某种趋势的话，那么最后剩下的r_i就可以看作信号的趋势。

5.希尔伯特-黄变换(HHT)多分辨分析：

从EMD的过程来看，EMD分解是先将信号的高频分量IMF₁先提取出来，再依次从高频到低频将各个IMF分量(IMF₂、IMF₃...IMF_n)提取出来，以f_imax表示IMF₁的最高频率，那么

f_imax＞f_(i+1)max，i＝1，2，...，n(16)

各分量的频率分辨率的公式为：

Δf_i＝Δf_i+1/N    (17)

其中N为信号采样点数，因此有：

Δf_i＞Δf_i+1(18)

由式(15)可知，每个IMF分量的频率分辨率式是不一样的，包含高频成分的频率分辨率较低，包含低频成分的频率分辨率较高，而且含有不同频率成分的各个IMF分量和它们对应的Hilbert谱中的最高频率及频率分辨率是自适应得到的。

本发明的实现步骤如下：

(1)对原始行波数据进行数据截断，选取故障前后1-5ms的数据；

(2)依据本发明原理部分说明3，对截断数据进行EMD分解；

(3)根据EMD分解原理，按式(16)f_imax＞f_(i+1)max，i＝1，2，...，n可知IMF₁包含原始数据的最高频信息，因此，选取IMF₁分量作为研究对象；

(4)对IMF₁做一阶差分，并且按式(1) $Y\left(t\right)=\frac{1}{\mathrm{\π}}{\∫}_{-\∞}^{+\∞}\frac{X\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}\mathrm{d\τ}$ 进行Hilbert变换，按式(5) $\mathrm{\θ}\left(t\right)=\arctan \frac{Y\left(t\right)}{X\left(t\right)}$ 和(6) $f\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{\mathrm{d\θ}\left(t\right)}{\mathrm{dt}}$ 求取其瞬时频率；

(5)寻找一阶差分曲线极值点和瞬时频率极大值点，一阶差分的极大值点对应行波波头上升最快点、极小值点对应行波波头下降最快点，即准确标定了行波波头到达时刻和极性，结合瞬时频率实现对突变幅值较小的波头的准确识别。

(6)在强噪干扰环境下，选取IMF₂或者IMF₃重复步骤(4)、(5)。

本发明与现有技术相比较具有的优点如下：

(1)EMD分解过程完全从信号本身的时间尺度出发，以纯时域的计算方法进行，不但很好的适应了信号本身的特点，而且计算速度快；

(2)正是由于EMD分解过程的自适应性，在第一个IMF分量中完全表现了行波波头的突变特征，没有任何由于滤波而产生的波头信息损失；

(3)对所考察的IMF分量进行一阶差分，其极值点明确标定了行波波头到达时刻并体现了波头极性，在精度上达到了其他方法所不可能达到的精度；

(4)结合对所考察的IMF分量的瞬时频率极大值点的判断，弥补了差分法对突变幅度较小的波头检测不灵敏的缺点，达到了在更复杂情况下行波波头的精确定位。

附图说明

图1HHT变换EMD分解过程；

图2本发明对故障行波仿真数据的HHT分析；

图3图2行波波头下降沿局部放大；

图4图2行波波头上升沿局部放大；

图5应用本发明对云南电网曲靖供电局虹桥变虹沾II回2008年5月29日20时43分25秒故障C相故障行波数据的EMD分解结果，包括前3个IMF分量；

图6应用本发明对图5的EMD分解结果IMF₁分量做一阶差分所得结果；

图7图6局部放大；

图8应用本发明对云南电网曲靖供电局沾益变沾虹II回2008年5月29日20时43分25秒故障C相故障行波数据的EMD分解结果，包括前3个IMF分量；

图9应用本发明对图8的EMD分解结果IMF₃分量做一阶差分的局部放大。

具体实施方式

实施例1：

此例为云南电网曲靖供电局虹桥变虹沾II回2008年5月29日20时43分25秒故障C相故障行波电流分析。

(1)在线路发生故障时，故障行波分析装置启动，进行快速、高采样率(1MHz)录波；

(2)对行波数据文件进行截断，见图5，选取初始波头前后2ms数据时窗；

(3)对截断数据按图1所示流程进行EMD分解，获得前几个IMF分量，EMD分解结果见图5(IMF1)、(IMF2)、(IMF3)曲线；

(4)选取IMF₁为考察对象，

按式(1)

$Y\left(t\right)=\frac{1}{\mathrm{\π}}{\∫}_{-\∞}^{+\∞}\frac{X\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}\mathrm{d\τ}---\left(1\right)$

对其做Hilbert变换，

按式(5)、(6)

$\mathrm{\θ}\left(t\right)=\arctan \frac{Y\left(t\right)}{X\left(t\right)}---\left(5\right)$

$f\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{\mathrm{d\θ}\left(t\right)}{\mathrm{dt}}---\left(6\right)$

求取其瞬时频率，并求其一阶差分，见图6(IMF)、(DIFF)曲线所示；

(5)见图7，可得一阶差分曲线极大值为1018点(1.018ms)，即精确定位故障行波初始波头到达时刻为1.018ms并标定其波头极性为正；

(6)在步骤(1)所选数据包含较强噪声干扰的情况下，为了避免噪声干扰，选取IMF₂或者IMF₃，重复步骤(4)、(5)。

实施例2：

此例为故障行波仿真数据HHT分析，包含不同极性的行波波头。

(1)在线路发生故障时，故障行波分析装置启动，进行快速、高采样率(1MHz)录波；

(2)对行波数据文件进行截断，见图2，选取初始波头前后1.5ms数据时窗；

(3)对截断数据按图1所示流程进行EMD分解，获得前几个IMF分量；

(4)选取IMF1为考察对象，按式(1)对其做Hilbert变换，按式(5)、(6)求取其瞬时频率，并对其求其一阶差分，见图2(IMF1分量瞬时频率)、(IMF1分量一阶差分)曲线所示；

(5)见图3，可得一阶差分曲线极小值点为514点(0.514ms)，即精确定位故障行波波头到达时刻为0.514ms，并标定其极性为负。

见图4，可得一阶差分曲线极大值点为774点(0.774ms)，即精确定位故障行波波头到达时刻为0.774ms，并标定其极性为正。

Claims (1)

1.一种利用HHT的输电线路故障行波初始波头到达时刻精确标定方法，其特征在于：利用HHT变换核心算法EMD分解分离行波信号高频部分，对其高频IMF分量求一阶差分曲线并作Hilbert变换求得其瞬时频率曲线，结合两者的特点精确定位故障行波初始波头到达时刻和极性信息，并按以下步骤进行：

(1)对原始行波数据进行数据截断，选取故障前后1-5ms的数据；

(2)对截断数据进行EMD分解；

(3)选取IMF1分量作为研究对象；

(4)对IMF1做一阶差分，并且按式

$Y\left(t\right)=\frac{1}{\mathrm{\π}}{\∫}_{-\∞}^{+\∞}\frac{X\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}\mathrm{d\τ}$

对IMF1进行Hilbert变换，按式

$\mathrm{\θ}\left(t\right)=\arctan \frac{Y\left(t\right)}{X\left(t\right)}$

$f\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{\mathrm{d\θ}\left(t\right)}{\mathrm{dt}}$

求取其瞬时频率；

(5)寻找IMF1分量一阶差分曲线极值点和瞬时频率极大值点，一阶差分的极大值点对应行波波头上升最快点、极小值点对应行波波头下降最快点，即准确标定了行波波头到达时刻和极性，结合瞬时频率实现对突变幅值较小的波头的准确识别。

Images