CN107329041B - 一种基于b样条双稳态去噪的配电网故障选线方法 - Google Patents

Abstract

本发明涉及一种基于B样条双稳态去噪的配电网故障选线方法，该方法包括如下步骤：(1)求取优化的B样条双稳态系统势函数及4阶龙格‑库塔算法的计算步长；(2)获取配电网各线路故障后1个工频周期的零序电流作为暂态零序电流；(3)将各线路的暂态零序电流导入B样条双稳态系统势函数并采用4阶龙格‑库塔算法求解得到各线路的特征暂态零序电流；(4)对各线路的特征暂态零序电流进行互相关系数计算并获得各线路的特征极性参数；(5)将特征极性参数最小的线路确定为故障线路。与现有技术相比，本发明可在强噪声背景下的微弱信号中实现准确选线。

Description

一种基于B样条双稳态去噪的配电网故障选线方法

技术领域

本发明涉及一种配电网故障选线方法，尤其是涉及一种基于B样条双稳态去噪的配电网故障选线方法。

背景技术

目前，故障选线方法主要利用小波变换、S变换、Prony算法、Hough变换等信号处理工具提取信号，然后，采用人工神经网络、支持向量机等建立选线判据。

小波变换具有良好的时域、频域局部化特性，能提取暂态信号在不同尺度的特征，但小波变换易受噪声影响。此外，不同的小波基函数将导致不同的暂态特征提取结果。S变换可确定出零序电流的主导特征频率，并通过比较主导特征频率上各馈线零序电流的S变换能量，以实现故障选线。Prony分段拟合故障T/4周期内的暂态零序电流信号，不仅有效避开电流互感器磁密饱和对采集信号的影响，而且在一定程度上提高了Prony整体拟合精度。Hough变换在参数空间不超过2维的情况下，有优异的表现，但若参数空间增大，其计算量急剧上升，同时耗费巨大的存储空间，耗时也随之猛增。

人工神经网络具有对周围环境自学习、自适应功能，可用于处理带噪声的、不完整的数据集，但存在局部最优，训练时间较长，可靠性有限等问题。支持向量机在解决小样本、非线性及高维模式识别问题中具有优势，但识别能力易受自身参数的影响。因此，对于故障选线，特别是强噪声背景下的微弱信号故障选线，还需进一步研究。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于B样条双稳态去噪的配电网故障选线方法。

本发明的目的可以通过以下技术方案来实现：

一种基于B样条双稳态去噪的配电网故障选线方法，该方法包括如下步骤：

(1)求取优化的B样条双稳态系统势函数及4阶龙格-库塔算法的计算步长；

(2)获取配电网各线路故障后1个工频周期的零序电流作为暂态零序电流；

(3)将各线路的暂态零序电流导入B样条双稳态系统势函数并采用4阶龙格-库塔算法求解得到各线路的特征暂态零序电流；

(4)对各线路的特征暂态零序电流进行互相关系数计算并获得各线路的特征极性参数；

(5)将特征极性参数最小的线路确定为故障线路。

所述的B样条双稳态系统势函数包括三种模式下的B样条双稳态系统势函数，分别为强故障模式、中等故障模式和弱故障模式，步骤(1)中求取每种故障模式下的优化的B样条双稳态系统势函数及4阶龙格-库塔算法的计算步长。

强故障模式：该模式下暂态零序电流的峰值大于等于50A，中等故障模式：该模式下暂态零序电流的峰值大于10A且小于50A，弱故障模式：该模式下暂态零序电流的峰值大于0A且小于等于10A。

所述的B样条双稳态系统势函数为：

a_G和b_G为G故障模式下的势函数参数，s(t)+Γ(t)为导入的暂态零序电流，x_G(t)为G故障模式下待求解的特征零序电流，t为采样时刻，G＝q、m、w，G＝q时表示强故障模式，G＝m时表示中等故障模式，G＝w时表示弱故障模式。

步骤(1)中采用遗传算法求取G故障模式下的优化的B样条双稳态系统势函数及4阶龙格-库塔算法的计算步长。

采用遗传算法求取B样条双稳态系统势函数及4阶龙格-库塔算法的计算步长具体为：

(11)设定G故障模式下B样条双稳态系统势函数中势函数参数a_G和b_G以及对应的4阶龙格-库塔算法数值计算步长h_G的范围，设置遗传算法的参数：包括最大迭代次数、种群大小、交叉概率为和变异概率；

(12)获取待优化种群，种群中第k个个体记作P_k＝[a_Gk,b_Gk,h_Gk]，k＝1,2……N，N为种群个体总数，将a_Gk,和b_Gk作为G故障模式下第k个双稳态系统的势函数参数a_G和b_G对应的值得到G故障模式下的N个双稳态系统，h_Gk作为第k个双稳态系统采用4阶龙格-库塔算法求解时的数值计算步长h_G对应的值；

(13)给定G故障模式下无噪声的暂态零序电流s(t)，s(t)＝[s(1)、s(2)…s(q)…s(Q)]，s(q)为G故障模式下无噪声的暂态零序电流在第q个采样点的值，q＝1,2……Q，Q为采样点总个数；

(14)对s(t)分别添加噪声强度为D的高斯白噪声Γ(t)获取含噪声暂态零序电流i_zg(t)，用i_zg(t)分别替换步骤(12)中G故障模式下的N个双稳态系统的势函数中的s(t)+Γ(t)；

(15)采用4阶龙格-库塔算法分别对N个双稳态系统的势函数求解，得到第k个双稳态系统输出的中间特征电流

为第k个双稳态系统输出的中间特征电流在第q个采样点的值；

(16)给定B样条函数，对第k个双稳态系统输出的中间特征电流进行B样条曲线拟合得到特征暂态零序电流

为第k个双稳态系统的特征暂态零序电流在第q个采样点的值；

(17)求取和s(t)之间的互相关系数ρ^k，选取ρ^k中的最大值作为群体互相关系数ρ_max，其中：

(18)判断ρ_max是否大于设定值，若是则输出ρ_max所对应参数a_max，b_max和h_max为最优化参数，否则执行步骤(19)；

(19)判断是否到达最大迭代次数，若是则输出ρ_max所对应参数a_max，b_max和h_max为最优化参数，否则进行交叉变异得到新种群，并将该新种群作为待优化种群，返回步骤(12)。

步骤(16)中B样条函数为4次B样条函数。

所述的步骤(2)具体为：采用8kHz的采样频率记录故障时刻起1个工频周期内第i条线路的暂态零序电流x_i(t)，x_i(t)＝[x_i(1)、x_i(2)…x_i(q)…x_i(Q)]，x_i(q)为第i条线路的暂态零序电流在第q个采样点的值，i＝1,2……l，l为线路总条数，q＝1,2……Q，Q为采样点总个数。

步骤(3)具体为：

(31)将第i条线路的暂态零序电流在第q个采样点的值x_i(q)分别替换3种模式下的B样条双稳态系统势函数中的s(t)+Γ(t)，分别求解得到和为G模式下第i条线路的过程零序电流在第q个采样点的值，G＝q、m、w，G＝q时表示强故障模式，G＝m时表示中等故障模式，G＝w时表示弱故障模式，q＝1,2……Q；

(32)按照峰值最大原则确定第i条线路的特征暂态零序电流在第q个采样点的值则第i条线路的特征零序电流为

步骤(4)具体为：

(41)计算各线路的特征暂态零序电流求取互相关系数矩阵C：

互相关系数矩阵C中第i行第j列元素记作为ρ_ij，ρ_ij表示第i条线路和第j条线路特征暂态零序电流之间的互相关系数，i＝1,2……l，j＝1,2……l，l为线路总条数，ρ_ij具体为：

其中，为第i条线路在第q个采样点对应的特征暂态零序电流，为第j条线路在第q个采样点对应的特征暂态零序电流，q＝1,2……Q，Q为采样点总个数；

(42)计算第i条线路的综合互相关系数S_i：

(43)计算第i条线路的特征极性参数J_i：

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

(1)本发明通过遗传算法和B样条拟合方法实现双稳态系统势函数的势函数参数以及4阶龙格-库塔算法的计算步长，能使在强噪声背景下经双稳态系统提取的特征信号保持其幅值、频率种类和波形变化趋势，实现强噪声背景下的微弱信号故障选线，提高故障选线准确率；

(2)采用B样条拟合方法改善了双稳态系统跟随有效自身输出信号与输出信号的噪声含量之间的矛盾，提高了势函数参数以及4阶龙格-库塔算法的计算步长的准确率，进而提高故障选线准确率；

(3)本发明B样条双稳态系统势函数设置三种模式，进而得到三种模式下的B样条双稳态系统势函数，在故障选线时，由于无法预知故障模式，因此分别通过三种模式下的B样条双稳态系统势函数求解过程零序电流，然后通过峰值原则来确定特征暂态零序电流，进而通过特征暂态零序电流来进一步进行故障选线，此种方式避免了只采用一个固定的B样条双稳态系统势函数进行求解特征暂态零序电流导致的不准确而影响故障选线的准确率，本发明进一步提高了故障选线的准确率。

附图说明

图1为本发明基于B样条双稳态去噪的配电网故障选线方法的流程框图；

图2为本实施例中暂态电流分布示意图；

图3为本实施例中单相接地零序暂态等值电路；

图4为本实施例中周期信号驱动下双稳态系统的R_out曲线图；

图5为本实施例中ATP仿真模型；

图6为本实施例案例A的相关电流图，图6(a)为案例A的无噪声暂态零序电流，图6(b)为案例A的含噪声暂态零序电流，图6(c)为案例A的暂态特征零序电流。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。

实施例

如图1所示，一种基于B样条双稳态去噪的配电网故障选线方法，该方法包括如下步骤：

(1)求取优化的B样条双稳态系统势函数及4阶龙格-库塔算法的计算步长；

(2)获取配电网各线路故障后1个工频周期的零序电流作为暂态零序电流；

(3)将各线路的暂态零序电流导入B样条双稳态系统势函数并采用4阶龙格-库塔算法求解得到各线路的特征暂态零序电流；

(4)对各线路的特征暂态零序电流进行互相关系数计算并获得各线路的特征极性参数；

(5)将特征极性参数最小的线路确定为故障线路。

B样条双稳态系统势函数包括三种模式下的B样条双稳态系统势函数，分别为强故障模式、中等故障模式和弱故障模式，步骤(1)中求取每种故障模式下的优化的B样条双稳态系统势函数及4阶龙格-库塔算法的计算步长。强故障模式：该模式下暂态零序电流的峰值大于等于50A，中等故障模式：该模式下暂态零序电流的峰值大于10A且小于50A，弱故障模式：该模式下暂态零序电流的峰值大于0A且小于等于10A。B样条双稳态系统势函数为：

W＝dx_G(t)/dt＝-dU(x_G)/dx_G+s(t)+Γ(t)， (1)

a_G和b_G为G故障模式下的势函数参数，s(t)+Γ(t)为导入的暂态零序电流，x_G(t)为G故障模式下待求解的特征零序电流，t为采样时刻，G＝q、m、w，G＝q时表示强故障模式，G＝m时表示中等故障模式，G＝w时表示弱故障模式。

步骤(1)中采用遗传算法求取G故障模式下的优化的B样条双稳态系统势函数及4阶龙格-库塔算法的计算步长。采用遗传算法求取B样条双稳态系统势函数及4阶龙格-库塔算法的计算步长具体为：

(11)设定G故障模式下B样条双稳态系统势函数中势函数参数a_G和b_G以及对应的4阶龙格-库塔算法数值计算步长h_G的范围，设置遗传算法的参数：包括最大迭代次数、种群大小、交叉概率为和变异概率，本实施例中a_G＝[-200,200]，b_G＝[0,100]，h_G＝[0.01,0.02]，最大迭代次数为100次，种群大小为100个，交叉概率为0.3，变异概率为0.1；

(12)获取待优化种群，种群中第k个个体记作P_k＝[a_Gk,b_Gk,h_Gk]，k＝1,2……N，N为种群个体总数，将a_Gk,和b_Gk作为G故障模式下第k个双稳态系统的势函数参数a_G和b_G对应的值得到G故障模式下的N个双稳态系统，h_Gk作为第k个双稳态系统采用4阶龙格-库塔算法求解时的数值计算步长h_G对应的值；

(13)给定G故障模式下无噪声的暂态零序电流s(t)，s(t)＝[s(1)、s(2)…s(q)…s(Q)]，s(q)为G故障模式下无噪声的暂态零序电流在第q个采样点的值，q＝1,2……Q，Q为采样点总个数；

(14)对s(t)分别添加噪声强度为D的高斯白噪声Γ(t)获取含噪声暂态零序电流i_zg(t)，用i_zg(t)分别替换步骤(12)中G故障模式下的N个双稳态系统的势函数中的s(t)+Γ(t)；

(15)采用4阶龙格-库塔算法分别对N个双稳态系统的势函数求解，得到第k个双稳态系统输出的中间特征电流

为第k个双稳态系统输出的中间特征电流在第q个采样点的值；

(16)给定B样条函数，对第k个双稳态系统输出的中间特征电流进行B样条曲线拟合得到特征暂态零序电流B样条函数为4次B样条函数：

为第k个双稳态系统的特征暂态零序电流在第q个采样点的值；

(17)求取和s(t)之间的互相关系数ρ^k，选取ρ^k中的最大值作为群体互相关系数ρ_max，其中：

(18)判断ρ_max是否大于设定值ρ，若是则输出ρ_max所对应参数a_max，b_max和h_max为最优化参数，否则执行步骤(19)，本实施例中ρ取值为0.85；

(19)判断是否到达最大迭代次数，若是则输出ρ_max所对应参数a_max，b_max和h_max为最优化参数，否则进行交叉变异得到新种群，并将该新种群作为待优化种群，返回步骤(12)。

本发明中可获取三类故障模式所对应的B样条双稳态系统的优化参数：1)强故障模式下：a_q＝-182.1387，b_q＝1.5362，h_q＝0.0106；2)中等故障模式下：a_m＝-147.5044，b_my＝7.3695，h_my＝0.0106；3)弱故障模式下：a_w＝-9.9026，b_w＝0.0005，h_w＝0.020。

步骤(2)具体为：采用8kHz的采样频率记录故障时刻起1个工频周期内第i条线路的暂态零序电流x_i(t)，x_i(t)＝[x_i(1)、x_i(2)…x_i(q)…x_i(Q)]，x_i(q)为第i条线路的暂态零序电流在第q个采样点的值，i＝1,2……l，l为线路总条数，q＝1,2……Q，Q为采样点总个数。

步骤(3)具体为：

(31)将第i条线路的暂态零序电流在第q个采样点的值x_i(q)分别替换3种模式下的B样条双稳态系统势函数中的s(t)+Γ(t)，分别求解得到和为G模式下第i条线路的过程零序电流在第q个采样点的值，G＝q、m、w，G＝q时表示强故障模式，G＝m时表示中等故障模式，G＝w时表示弱故障模式，q＝1,2……Q；

(32)按照峰值最大原则确定第i条线路的特征暂态零序电流在第q个采样点的值则第i条线路的特征零序电流为

步骤(4)具体为：

(41)计算各线路的特征暂态零序电流求取互相关系数矩阵C：

互相关系数矩阵C中第i行第j列元素记作为ρ_ij，ρ_ij表示第i条线路和第j条线路特征暂态零序电流之间的互相关系数，i＝1,2……l，j＝1,2……l，l为线路总条数，ρ_ij具体为：

其中，为第i条线路在第q个采样点对应的特征暂态零序电流，为第j条线路在第q个采样点对应的特征暂态零序电流，q＝1,2……Q，Q为采样点总个数；

(42)计算第i条线路的综合互相关系数S_i：

(43)计算第i条线路的特征极性参数J_i：

具体地，本发明的工作原理如下：

1、小电流接地故障特征分析

小电流接地系统发生单相接地故障时，其电流分布示意如图2所示。当采用消弧线圈以后，单相接地时的电流分布将发生重大的变化。图2中，在电源的中性点接入了消弧线圈，当线路2上A相接地以后，电容电压的大小和分布与不接消弧线圈时是一样的，区别之处是在接地点又增加了一个电感分量的电流因此，从接地点流回的总电流为：

式中：为全系统的对地电容电流，为消弧线圈的电流，若用L表示它的电感，则：

由于和的相位相差180°，因此将因消弧线圈的补偿而减小。且当消弧线圈运行于过补偿方式时，接地点残余电流呈感性，流经故障线路电容性无功功率的方向是由母线流向线路，与非故障线路相同，故障信息比较微弱。

利用单相接地零序暂态等值电路进行暂态零序电流分析，如图3所示。其中：C₀为线路零序电容，L₀为线路零序等值电感，R_g为接地点的过渡电阻，R_p和L_p分别为消弧线圈的等效电阻和电感，e(t)为零序电压。

在补偿电网中发生故障的瞬间，由图3可得流过故障点的暂态零序电流i_0.t为：

式中：i_0L.t为暂态零序电流中的电感电流分量，i_0C.t为暂态零序电流的电容电流分量，I_Lm和I_Cm分别为电感电流和电容电流的初值(I_Cm＝U_phmωC，I_Lm＝U_phm/ωL)，U_phm为相电压的幅值，ω为工频角频率，ω_f和δ分别为暂态零序电流容性分量的振荡角频率和衰减系数，τ_L为电感电流的衰减时间常数，为接地时故障线路相电压的初始相位。

由以上分析可知，当单相接地故障发生后，在故障点便有衰减很快的暂态电容电流和衰减很慢的暂态电感电流流过。暂态接地电流的幅值和频率均主要由暂态电容电流所确定。线路越长时，自振荡频率越低，暂态电容电流的自由振荡分量的幅值也会降低，同时，自由振荡的持续时间一般也会减少至半个工频周波左右。再者配电网的结构本身就是一个小电流接地系统，实际运行中又由于噪声干扰，信噪比小，进而使故障线路的暂态零序电流与非故障线路的差别不大，含有明显故障信息的信号微弱，给配电网准确选线带来了一定难度。

2、龙格-库塔方法

为了避免计算高阶导数，龙格-库塔方法利用f(x,y)在某些点处的值的线性组合，构造一类计算公式，使其按泰勒级数展开后，与初值问题的解的泰勒展开式比较，存在尽可能多的项完全相同，从而保证算式有较高的精度。这种方法间接利用了泰勒展开的思想，避免了计算高阶导数的困难。

一般的龙格-库塔方法的形式为：

其中α₂,α₃,…α_n，c₁,c₂,…c_n，μ₂₁,μ₃₁,μ₃₂,…μ_n1,μ_n2,…μ_n,n-1,均为待定参数，g为计算步长，f()表示待求解函数，选取上述待定参数这些参数的是要求式(14)中第1式右端在(x_n,y_n)处作泰勒展开式，并按g的幂次从低到高的排列式与微分方程解的泰勒展开式有尽可能多的项重合，也就是要求符合式(15)：

χ₁＝f_n,χ₂＝f_n',χ₃＝f_n”，…， (15)

这里f_n,f_n',f_n”，…，表示y'(x_n)＝f(x_n,y_n)，y”(x_n)，y”'(x_n),…。通常把式(14)称为N级龙格-库塔方法，简记为N级龙格-库塔方法。更高阶的龙格-库塔方法由于计算量较大，一般不采用。本发明采用4阶龙格-库塔方法，其计算式如式(16)所示：

3、双稳态系统理论

弱暂态零序电流的随机共振不同于周期信号，常用的随机共振双稳态系统可表示为：

W＝dx/dt＝-dU(x)/dx+s(t)+Γ(t)， (17)

式中：U(x)为势函数，U(x)＝-a(x(t))²/2+b(x(t))⁴/4，a和b为势函数参数，s(t)为输入信号，t为时间，x(t)为输出信号，Γ(t)的强度为D。

当s(t)为周期信号Asin(2πf₀t)、b＝1时，在绝热近似条件下(即噪声强度D、信号幅值A和频率f₀都远小于1，可求得系统的输出信噪比R_out为：

当A＝0.2、a＝1时改变D值和当A和D依次为0.2、0.15时改变a值，可分别得到R_out随D与a变化的曲线，如图4所示。图4曲线的形状类似于力学中的单峰共振曲线，且在适当的噪声强度D或势函数参数下，加入噪声会使输入端的信噪比降低，但能使其输出端的信噪比R_out增加，改变势函数参数还可进一步提高R_out，这就是随机共振现象。

但是，当引入非周期信号s(t)的作用后，势函数将受到其调制，也即为：

U″(x)＝-a+3bx²， (20)

进而得到双稳态系统单势阱近似模型为：

x＝-U″(x_s)(x-x_s)+Γ(t)， (21)

式中：x_s为满足-ax+bx³-s(t)＝0的一个稳态解，x即为x(t)。

式(21)所对应的系统输出的概率密度函数ρ(x,t)满足Fokker–Planck方程，即：

式(22)中的最小非零特征值U″(x_s)决定了系统响应的快慢，对系统的瞬时行为起决定作用。

当系统响应速度足够快时，可假设输入恒定s(t)＝s₀，也即系统稳态输出可看成一个确定成分x₀和随机成分Ψ(t)的组合，系统输出可进一步表示为：

s_out≈s₀+U″(x₀)ψ(t)， (23)

s_out的方差为：

D₁(s_out)＝U″(x_s)D， (24)

由式(24)可得输出方差D₁(s_out)与输入噪声方差D(Γ(t))的比值为：

D₁(s_out)/D(Γ(t))＝U″(x_s)h/2， (25)

式中：h为信号采样间隔。

由式(25)和式(22)对比可知，增大U″(x_s)可使双稳态系统跟随信号变化，即信号失真减小。但增大U″(x_s)将导致输出方差与输入噪声方差的比值增加，即增加信号中的噪声成分。而U″(x_s)的大小又由势函数参数a和b所决定，因此，势函数参数对于非周期信号的检测性能起着关键作用。

4、遗传算法基本原理

弱暂态零序电流的随机共振不同于周期信号，它对系统参数的合理选取更为敏感。此外，在变尺度思想条件下，四阶龙格—库塔算法数值计算步长h的选择也十分关键。据已有文献可知，遗传算法本质上是一种启发式的随机搜索算法，由遗传算法得出的结果每次都不尽相同，再加上自变量在给定的约束条件下进行了无缝编码，所以从理论上讲，遗传算法总有很多机会得到全局最优结果。因此，本发明采用遗传算法对双稳态系统表达式中势函数的参数a、b以及四阶龙格—库塔算法的数值计算步长h进行优化，使得变尺度双稳态系统能够更好地检测强噪声背景下的弱暂态零序电流。

遗传算法的概念是由Michigan大学的J.Holland教授在1962年提出的。遗传算法是模拟达尔文生物进化论的自然选择和孟德尔遗传学机理的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法。它把自然界“优胜劣汰，适者生存”的生物进化原理引入优化参数形成的编码串联群体中，按照所选择的适应度函数并通过遗传中的选择、交叉和变异对个体进行筛选，使适应度值好的个体被保留，适应度差的个体被淘汰，新的群体既继承了上一代的信息，又优于上一代。这样反复循环，直至满足条件。遗传算法基本的操作分为：

选择操作：选择操作是指从旧群体中以一定概率选择个体到新群体中，个体被选中的概率跟适应度值有关，个体适应度值越好，被选中的概率越大。遗传算法选择操作有轮盘赌法、锦标赛法等多种方法，本发明选择轮盘赌法，即基于适应度比例的选择策略，每个个体j的选择概率P_j：

o_j＝k₁S_j， (26)

式中：S_j为个体j的适应度值，本发明中令互相关系数ρ^k为S_j，k₁为系数，N为种群个体数目，j表示个体编号。

交叉操作：交叉操作是指从个体中选择两个个体，通过两个染色体的交换组合，来产生新的优秀个体。交叉过程为从群体中任选两个染色体，随机选择一点或多点染色体位置进行交换。由于个体采用实数编码，所以交叉操作方法采用实数交叉法。

变异操作：变异操作是指从群体中任选一个个体，选择染色体中的一点进行变异以产生更优秀的个体。选择第j个个体的第m个基因r_jm进行变异。

5、B样条拟合原理

工程中应用的拟合曲线有2类：1)生成的曲线通过所有给定点，如抛物样条曲线和三次参数样条曲线；2)不一定通过给定的点，而是比较好地接近已给定的点，代表性曲线为Bezier曲线与B样条曲线。Bezier曲线的形状通过一组特征多边折线的顶点确定，且只有第1个顶点和最后一个顶点在曲线上，其余顶点则用于定义曲线的阶数、导数及形状，第1条和最后一条边为曲线在两端点处的切线方向。因此Bezier曲线在应用中存在一定不足，特征多边形顶点个数确定后，曲线的阶次无法更改，导致Bezier曲线缺乏灵活，当顶点数过多时，曲线的阶次较高，特征多边形对曲线形状的控制能力降低。另外，所定义的曲线在定义域内均受到全部顶点的影响，使得曲线不能进行局部修改。B样条曲线是在保留Bezier曲线优点的同时，克服其由于整体表示带来的不能局部修改的缺点，及解决在描述复杂形状时带来的连接问题下提出的灵活曲线。

B样条曲线表达式如下：

式中，0≤δ≤1，k＝0,1,2,...,β，α＝0,1,2,...,M-β，M为控点个数，β为样条曲线的阶次，整条曲线由M-β段B样条曲线平滑连接而成，每段曲线由β+1个控制点生成，P_α+k是第(α+k)个控制点的坐标，F_k,β为β次B样条基函数，P_α,β为第α段曲线上任意一点坐标。

在得到一组有序电流数据点Q(t)(t＝1,2,...,N)，N为采样点总个数，若要想点Q(t)在B样条曲线上，也即点Q(t)要满足方程(29)，于是有：

将式(30)写成矩阵形式有：

Q＝FP， (31)

其中，Q为包含N个数据点的N×β矩阵，F是N×M的B样条基函数系数矩阵，P是包含M个未知控制点的M×β矩阵。

进一步可以推出高斯正交方程组：

F^TFP＝F^TQ， (32)

由于F是满秩的，所以F^TF是M阶对称可逆矩阵，方程存在唯一解。得到最小二乘拟合的B样条曲线的控制点，求解方程如下：

P＝(F^TF)^-1F^TQ， (33)

可以证明，在最小二乘法原理意义下式(33)的解是最优解，即解出的控制点形成的B样条曲线是最佳拟合曲线。

为确保拟合精度和降低计算量，本发明选择4次B样条曲线进行电流曲线的拟合，样条基函数如下：

本实施例利用ATP做单相接地仿真试验，仿真模型如图5所示，其中线路l₁和l₂为架空线，线路长度分别为13.5km、24km，线路l₃为缆-线混合线路，其中电缆线长度为5km，架空线长度为12km，线路l₄为电缆线，长度10km。模型具体电气参数如下：

线路：架空线路正序参数R₁＝0.17Ω/km，L₁＝1.2mH/km，C₁＝9.697nF/km；零序参数R₀＝0.23Ω/km，L₀＝5.48mH/km，C₀＝6nF/km。电缆线路正序参数R₁₁＝0.193Ω/km，L₁₁＝0.442mH/km，C₁₁＝143nF/km；零序参数R₀₀＝1.93Ω/km，L₀₀＝5.48mH/km，C₀₀＝143nF/km。变压器：110/10.5kV；高压侧单相中性点线圈电阻0.40Ω，电感12.2Ω；低压侧单相线圈电阻0.006Ω，电感0.183Ω；励磁电流0.672A，励磁磁通202.2Wb，磁路电阻400kΩ。负荷：一律采用三角形接法，Z_L＝400+j20Ω。消弧线圈：在消弧线圈接地系统仿真时，消弧线圈电感为L_N＝1281.9mH。

其中，消弧线圈的电阻值取电抗值的10％，经计算为40.2517Ω。仿真模型采样频率f＝10⁵Hz，仿真时长0.06s，故障发生时刻设定为0.02s。

为验证本发明一种利用改进振子系统的故障选线方法的正确性，分别在以下故障情况下进行验证：不同接地电阻值及故障初相角、不同信噪比、不同线路故障。

①不同接地电阻值及故障初相角情况下：以线路l₁在不同接地电阻值及故障初相角故障时为例，具体计算结果如表1所示：

以案例A为例，来说明B样条双稳态系统去噪的有效性，相关的电流波形见图6。案例A:故障线路为线路l₁，故障距离为距离母线5km，故障初相角为90°，接地电阻为1800Ω。

在图6(c)与图6(b)相比较，图6(c)的噪声含量明显少于图6(c)。以互相关系数来说，线路l₃的互相关系数矩阵中的C₃为例，图6(a)的C₃为[-0.242,0.950,1.000,0.963]，图6(b)的C₃为[-0.219,0.379,1.000,0.834]，图6(c)的C₃为[-0.316,0.836,1.000,0.979]。由图6(a)和图6(b)的互相关系数矩阵C₃可知，噪声使特征暂态零序电流之间的互相关系数减小，也就是为什么相关理论会在选线中失效的一个原因。在图6(c)中，图6(c)的与图6(a)的基本一致。因此，本发明所提的B样条双稳态系统去噪的有效性。

在表1中，线路l₁的J₁是最小的，所以线路1被判定为故障线路，这与实际故障一致。表1中的数据表明了本发明所提选线方法的正确性。

表1不同故障电阻下的选线结果

②不同信噪比情况下：当接地电阻为500Ω，故障线路为l₃，故障距离为5km，故障初相角为0°，信噪比依次为-1db,5db,20db和50db时，按照本发明所提方法，具体计算结果如表2所示，由表2可知本发明的有效性。

表2不同信噪比下的选线结果

③不同线路故障情况下：本发明中线路2为架空线路，线路3为缆线混合线路，线路4为电缆线路。当这三条线路发生不同的接地电阻故障时，其选线结果见表3。表3表明了本发明所提的选线方法不受线路类型的影响。

表3不同故障线路的选线结果

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于B样条双稳态去噪的配电网故障选线方法，其特征在于，该方法包括如下步骤：

(1)求取优化的B样条双稳态系统势函数及4阶龙格-库塔算法的计算步长；

(2)获取配电网各线路故障后1个工频周期的零序电流作为暂态零序电流；

(3)将各线路的暂态零序电流导入B样条双稳态系统势函数并采用4阶龙格-库塔算法求解得到各线路的特征暂态零序电流；

(4)对各线路的特征暂态零序电流进行互相关系数计算并获得各线路的特征极性参数；

(5)将特征极性参数最小的线路确定为故障线路；

所述的B样条双稳态系统势函数包括三种模式下的B样条双稳态系统势函数，分别为强故障模式、中等故障模式和弱故障模式，步骤(1)中求取每种故障模式下的优化的B样条双稳态系统势函数及4阶龙格-库塔算法的计算步长；

所述的B样条双稳态系统势函数为：

W＝dx_G(t)/dt＝-dU(x_G)/dx_G+s(t)+Γ(t)，

a_G和b_G为G故障模式下的势函数参数，s(t)+Γ(t)为导入的暂态零序电流，x_G(t)为G故障模式下待求解的特征零序电流，t为采样时刻，G＝q、m、w，G＝q时表示强故障模式，G＝m时表示中等故障模式，G＝w时表示弱故障模式；

所述的步骤(2)具体为：采用8kHz的采样频率记录故障时刻起1个工频周期内第i条线路的暂态零序电流x_i(t)，x_i(t)＝[x_i(1)、x_i(2)…x_i(q)…x_i(Q)]，x_i(q)为第i条线路的暂态零序电流在第q个采样点的值，i＝1,2……l，l为线路总条数，q＝1,2……Q，Q为采样点总个数；

步骤(3)具体为：

(31)将第i条线路的暂态零序电流在第q个采样点的值x_i(q)分别替换3种模式下的B样条双稳态系统势函数中的s(t)+Γ(t)，分别求解得到和为G模式下第i条线路的过程零序电流在第q个采样点的值，G＝q、m、w，G＝q时表示强故障模式，G＝m时表示中等故障模式，G＝w时表示弱故障模式，q＝1,2……Q；

(32)按照峰值最大原则确定第i条线路的特征暂态零序电流在第q个采样点的值则第i条线路的特征零序电流为

2.根据权利要求1所述的一种基于B样条双稳态去噪的配电网故障选线方法，其特征在于，强故障模式：该模式下暂态零序电流的峰值大于等于50A，中等故障模式：该模式下暂态零序电流的峰值大于10A且小于50A，弱故障模式：该模式下暂态零序电流的峰值大于0A且小于等于10A。

3.根据权利要求1所述的一种基于B样条双稳态去噪的配电网故障选线方法，其特征在于，步骤(1)中采用遗传算法求取G故障模式下的优化的B样条双稳态系统势函数及4阶龙格-库塔算法的计算步长。

4.根据权利要求3所述的一种基于B样条双稳态去噪的配电网故障选线方法，其特征在于，采用遗传算法求取B样条双稳态系统势函数及4阶龙格-库塔算法的计算步长具体为：

(11)设定G故障模式下B样条双稳态系统势函数中势函数参数a_G和b_G以及对应的4阶龙格-库塔算法数值计算步长h_G的范围，设置遗传算法的参数：包括最大迭代次数、种群大小、交叉概率为和变异概率；

(12)获取待优化种群，种群中第k个个体记作P_k＝[a_Gk,b_Gk,h_Gk]，k＝1,2……N，N为种群个体总数，将a_Gk,和b_Gk作为G故障模式下第k个双稳态系统的势函数参数a_G和b_G对应的值得到G故障模式下的N个双稳态系统，h_Gk作为第k个双稳态系统采用4阶龙格-库塔算法求解时的数值计算步长h_G对应的值；

(13)给定G故障模式下无噪声的暂态零序电流s(t)，s(t)＝[s(1)、s(2)…s(q)…s(Q)]，s(q)为G故障模式下无噪声的暂态零序电流在第q个采样点的值，q＝1,2……Q，Q为采样点总个数；

(14)对s(t)分别添加噪声强度为D的高斯白噪声Γ(t)获取含噪声暂态零序电流i_zg(t)，用i_zg(t)分别替换步骤(12)中G故障模式下的N个双稳态系统的势函数中的s(t)+Γ(t)；

(15)采用4阶龙格-库塔算法分别对N个双稳态系统的势函数求解，得到第k个双稳态系统输出的中间特征电流

为第k个双稳态系统输出的中间特征电流在第q个采样点的值；

(16)给定B样条函数，对第k个双稳态系统输出的中间特征电流进行B样条曲线拟合得到特征暂态零序电流

为第k个双稳态系统的特征暂态零序电流在第q个采样点的值；

(17)求取和s(t)之间的互相关系数ρ^k，选取ρ^k中的最大值作为群体互相关系数ρ_max，其中：

(18)判断ρ_max是否大于设定值，若是则输出ρ_max所对应参数a_max，b_max和h_max为最优化参数，否则执行步骤(19)；

(19)判断是否到达最大迭代次数，若是则输出ρ_max所对应参数a_max，b_max和h_max为最优化参数，否则进行交叉变异得到新种群，并将该新种群作为待优化种群，返回步骤(12)。

5.根据权利要求4所述的一种基于B样条双稳态去噪的配电网故障选线方法，其特征在于，步骤(16)中B样条函数为4次B样条函数。

6.根据权利要求1所述的一种基于B样条双稳态去噪的配电网故障选线方法，其特征在于，步骤(4)具体为：

(41)计算各线路的特征暂态零序电流求取互相关系数矩阵C：

互相关系数矩阵C中第i行第j列元素记作为ρ_ij，ρ_ij表示第i条线路和第j条线路特征暂态零序电流之间的互相关系数，i＝1,2……l，j＝1,2……l，l为线路总条数，ρ_ij具体为：

其中，为第i条线路在第q个采样点对应的特征暂态零序电流，为第j条线路在第q个采样点对应的特征暂态零序电流，q＝1,2……Q，Q为采样点总个数；

(42)计算第i条线路的综合互相关系数S_i：

(43)计算第i条线路的特征极性参数J_i：