基于小波分析的绝缘子污闪泄漏电流信号稀疏表示方法

摘要

本发明公开了电力系统信号处理技术领域中的一种基于小波分析的绝缘子污闪泄漏电流信号稀疏表示方法。包括采用sym8正交小波基为稀疏基Ψ，构造小波稀疏变换正交矩阵

说明书

技术领域

本发明属于电力系统信号处理技术领域，尤其涉及一种基于小波分析的绝缘子污闪泄漏电流信号稀疏表示方法。

背景技术

在绝缘子污闪泄漏电流的实时监测中，污闪泄漏电流信号采样速率和处理速度要求越来越高，庞大的数据量给泄漏电流信号的采集传输和存储造成很大压力。压缩感知（CS）通过利用泄漏电流信号的可稀疏特性，在远小于Nyquist采样率的条件下，利用较少的测量值恢复原始信号，从而缓解了上述数据采集传输和存储压力。

基于CS思想，信号稀疏表示的稀疏分解过程中，信号通常分解在一组完备的正交基上，主要有傅立叶变换、Gabor变换、小波变换等。但是傅立叶变换只能在时、频域的某一域上对信号进行全局观察，还要求信号满足平稳特性。Gabor变换解决了对信号局部化问题，但是无法实现多分辨率分析。

本发明提出基于小波分析的绝缘子污闪泄漏电流信号稀疏表示方法，即采用小波变换算法在sym8小波基上对绝缘子泄漏电流的进行稀疏分解，并且通过采用正交匹配追踪法进行重构，从而实现了信号的稀疏表示。本发明不但降低了传感器采样频率的要求，而且在不影响信号质量的前提下提高了信号的采集和重构速度，为进行污闪泄漏电流的实时监测的实现打下了很好的基础。

发明内容

本发明的目的在于，提供一种基于小波分析的绝缘子污闪泄漏电流信号稀疏表示方法，在降低传感器采样频率和不影响信号质量的同时，可以提高泄漏电流信号的采集和重构速度，并且提高重构的准确性。

为了实现上述目的，本发明提出的技术方案是，一种基于小波分析的绝缘子污闪泄漏电流信号稀疏表示方法，其特征是所述方法包括：

步骤1：选择稀疏基并计算原始泄漏电流信号的稀疏度；

步骤2：构造高斯随机矩阵Φ，并计算测量向量y和传感矩阵Θ；

步骤3：利用测量向量y和传感矩阵Θ重构原始泄漏电流信号。

所述步骤1包括：

步骤:11：选择sym8正交小波基作为稀疏基，根据sym8正交小波基分别设计低通滤波器和高通滤波器；

步骤12：分别使用低通滤波器和高通滤波器对原始泄露电流信号进行从1到log₂(N)层的小波分解，得到各层的低通滤波器分解系数h和高通滤波器分解系数g，N为原始泄露电流信号的长度；

步骤13：构造小波稀疏变换正交矩阵

步骤14：根据公式计算原始泄漏电流信号的小波稀疏系数向量，其中x为原始泄漏电流信号在sym8正交小波基上的N×1维信号；

步骤15：将大于设定阈值ε的小波稀疏系数的数目作为原始泄漏电流信号的稀疏度K。

所述子步骤13包括：

子步骤1301：令j=log₂L+1，L为低通滤波器分解系数h的长度，令初始正交矩阵为单位矩阵I；

子步骤1302：利用第j层的低通滤波器分解系数提取原始泄露电流信号的低频特征系数，得到低频特征系数向量H_j；利用第j层的高通滤波器分解系数提取原始泄露电流信号的高频特征系数，得到高频特征系数向量G_j；

子步骤1303：令i=1；

子步骤1304：对低频特征系数向量H_j和高频特征系数向量G_j进行周期性延拓，即对低频特征系数向量H_j和高频特征系数向量G_j分别增加2^j-L列零向量，得到低频特征系数矩阵H_j′和高频特征系数矩阵G_j′；

子步骤1305：分别对低频特征系数矩阵H_j′和高频特征系数矩阵G_j′中的列向量进行单步圆周位移；

子步骤1306：判断i>2^j-1是否成立，如果i>2^j-1则执行步骤1307；否则，令i=i+1，返回步骤1304；

子步骤1307：对低频特征系数矩阵H_j′和高频特征系数矩阵G_j′进行正交化处理，得到正交矩阵

子步骤1308：判断j>log₂N是否成立，如果j>log₂N，则执行子步骤1309；否则令且令j=j+1，返回子步骤1302；

子步骤1309：此时结束j的迭代，得到最终的正交矩阵并令即为小波稀疏变换正交矩阵。

所述步骤2包括：

步骤21：根据原始泄漏电流信号的稀疏度K计算测量数M；M为距离Klog(N/K)最近的正整数；

步骤22：构造阶数为M×N的高斯随机矩阵Φ，且高斯随机矩阵Φ与稀疏基矩阵Ψ不相关；

步骤23：根据公式y=Φf计算测量向量y；其中，f为原始泄漏电流信号；

步骤24：根据公式Θ=ΦΨ计算传感矩阵Θ；其中，Ψ为稀疏基矩阵；

所述步骤3包括：

步骤31：分别设定冗余向量初始值γ₀、索引集合初始值Λ₀和初始迭代次数t；

其中，冗余向量初始值γ₀=y；

索引集合初始值Λ₀为空集；

初始迭代次数t=0；

步骤32：根据公式计算索引值δ_t，并令Λ_t=Λ_t-1∪{δ_t}；

其中，u_j为相关系数且u_j=|<γ_t,θ_j>|；

θ_j为传感矩阵Θ的第j列列向量；

N为原始泄露电流信号的长度；

步骤33：根据公式计算使残差y-Θx_i最小时的信号估计值其中，x_i为时的原始泄露电流信号；

步骤34：根据公式计算第t次迭代的冗余向量γ_t；

步骤35：判断t<2×K是否成立，如果t<2×K，则令t=t+1，返回步骤32；否则，执行子步骤36；

步骤36：将和小波稀疏变换正交矩阵的逆矩阵做逆小波变换得到重构的时域信号

本发明通过在sym8小波基上对绝缘子泄漏电流信号的进行稀疏分解，并且通过采用正交匹配追踪法进行重构，达到了信号的稀疏表示；并且，本发明不但降低了传感器采样频率的要求，而且在不影响信号质量的前提下提高了信号的采集和重构速度，同时保证了重构的准确性。

附图说明

图1是基于小波分析的绝缘子污闪泄漏电流信号稀疏表示方法流程图；

图2是实施例提供的绝缘子原始泄漏电流信号图；

图3是实施例提供的绝缘子原始泄漏电流信号的重构图。

具体实施方式

下面结合附图，对优选实施例作详细说明。应该强调的是，下述说明仅仅是示例性的，而不是为了限制本发明的范围及其应用。

实施例1

图1是基于小波分析的绝缘子污闪泄漏电流信号稀疏表示方法流程图。下面结合图1对本发明提供的方法的原理进行说明。如图1所示，本发明提供的基于小波分析的绝缘子污闪泄漏电流信号稀疏表示方法包括：

步骤1：选择稀疏基并计算原始泄漏电流信号的稀疏度。

压缩感知的原理是，对于长度为N的原始泄露电流信号f，用一组基向量Ψ=[ψ₁,ψ₂,..,ψ_N]的线性组合表示，即Ψ是N×N稀疏基矩阵。x是原始泄露电流信号f在Ψ域上的等价表示，当x仅有K（K<<N）个非零系数时，则信号是稀疏的，此时称Ψ为信号f的稀疏基，而此时的稀疏系数个数成为稀疏度K。

步骤1所述的过程包括如下步骤：

步骤11：选择sym8正交小波基作为稀疏基，根据sym8正交小波基分别设计低通滤波器和高通滤波器。

步骤12：分别使用低通滤波器和高通滤波器对原始泄露电流信号进行从1到log₂(N)层的小波分解，得到各层的低通滤波器分解系数h和高通滤波器分解系数g，N为原始泄露电流信号的长度。

在给定小波基上求分解滤波器系数已经有成熟的实现方法，一般的数学计算软件（如MATLAB）中都提供该算法。利用MATLAB中的wfilters函数，即可计算出sym8基的低通滤波器分解系数h和高通滤波器分解系数g。

步骤13：构造小波稀疏变换正交矩阵其具体包括如下子步骤：

子步骤1301：令j=log₂L+1，L为低通滤波器分解系数h的长度。同时，令初始正交矩阵为单位矩阵I。

子步骤1302：利用第j层的低通滤波器分解系数提取原始泄露电流信号的低频特征系数，得到低频特征系数向量H_j；利用第j层的高通滤波器分解系数提取原始泄露电流信号的高频特征系数，得到高频特征系数向量G_j。

子步骤1303：令i=1。

子步骤1304：对低频特征系数向量H_j和高频特征系数向量G_j进行周期性延拓，即对低频特征系数向量H_j和高频特征系数向量G_j分别增加2^j-L列零向量，得到低频特征系数矩阵H_j′和高频特征系数矩阵G_j′。

子步骤1305：分别对低频特征系数矩阵H_j′和高频特征系数矩阵G_j′中的列向量进行单步圆周位移。单步圆周位移就是将矩阵中的列向量向右（或者向左）移动一列，而最右（最左）列移动到最左（或者最右）列上。

子步骤1306：判断i>2^j-1是否成立，如果i>2^j-1则执行步骤1307；否则，令i=i+1，返回步骤1304。

子步骤1307：对低频特征系数矩阵H_j′和高频特征系数矩阵G_j′进行正交化处理，得到正交矩阵

子步骤1308：判断j>log₂N是否成立，如果j>log₂N，则执行子步骤1309；否则令且令j=j+1，返回子步骤1302。

子步骤1309：此时结束j的迭代，得到最终的正交矩阵并令即为小波稀疏变换正交矩阵。

通过这种方式构造小波稀疏变换正交矩阵是为了将分解得到的各子带数据分别落到相互正交的L2(R)的子空间中，使各子带数据相关性减小，有利于精确重建原始信号。

步骤14：根据公式计算原始泄漏电流信号的小波稀疏系数向量，其中x为原始泄漏电流信号在sym8正交小波基上的N×1维信号。

步骤15：将大于设定阈值ε的小波稀疏系数的个数作为原始泄漏电流信号的稀疏度K。

阈值ε的设定与信号的特性和变化相关，本发明为了简单计算，通常采用全局阈值的方法，确定一个统一的阈值；同时，为了使统一的阈值可以根据不同的信号特征进行适应性选择，本发明采用了划分稀疏系数的区间然后取同一区间系数平均值为阈值的方法。

步骤2：构造高斯随机矩阵Φ，并计算测量向量y和传感矩阵Θ。

步骤2的实现过程具体包括：

步骤21：根据原始泄漏电流信号的稀疏度K计算测量数M，M为距离Klog(N/K)最近的正整数，采用如下公式计算：

$>M=\left(\begin{array}{c}\left[K\log (N/K)\right],K\log (N/K)-\left[K\log (N/K)\right]\le 0.5\\ \left[K\log (N/K)\right]+1,K\log (N/K)-\left[K\log (N/K)\right]>0.5\end{array}\right)$>

其中，[Klog(N/K)]是Klog(N/K)的取整运算，即当Klog(N/K)比其取整结果大0.5时，说明Klog(N/K)距离[Klog(N/K)]+1更近。反之，说明Klog(N/K)距离[Klog(N/K)]更近。

步骤22：构造阶数为M×N的高斯随机矩阵Φ，且高斯随机矩阵Φ与稀疏基矩阵Ψ不相关。

步骤23：根据公式y=Φf计算测量向量y；其中，f为原始泄漏电流信号。

由于y能恢复原始信号的充要条件是满足约束等距准则（RIP），而RIP的等价条件稀疏变换基Ψ和观测阵Φ不相关，而本文中稀疏变换基Ψ已经固定，所以需要观测阵满足RIP。有理论证明高斯随机矩阵正好满足RIP，则设观测矩阵Φ为N×M维高斯随机矩阵，代入表达式y=ΦΨx得到线性测量，即测量向量y。

步骤24：根据公式Θ=ΦΨ计算传感矩阵Θ；其中，Ψ为稀疏基矩阵。

步骤3：利用测量向量y和传感矩阵Θ重构原始泄漏电流信号。

重构原始泄露电流信号的过程又包括：

步骤31：分别设定冗余向量初始值γ₀、索引集合初始值Λ₀和初始迭代次数t。其中，冗余向量初始值γ₀=y，索引集合初始值Λ₀为空集，初始迭代次数t=0；

步骤32：根据公式计算索引值δ_t，并令Λ_t=Λ_t-1∪{δ_t}。其中，u_j为相关系数且u_j=|<γ_t,θ_j>|，即冗余向量与传感矩阵相应列向量的内积的绝对值，θ_j为传感矩阵Θ的第j列列向量，N为原始泄露电流信号的长度。

步骤33：根据公式采用最小二乘法计算使残差y-Θx_i最小时的信号估计值其中，x_i为时的原始泄露电流信号。

步骤34：根据公式计算第t次迭代的冗余向量γ_t。

步骤35：判断t<2×K是否成立，如果t<2×K，则令t=t+1，返回步骤32；否则，执行子步骤36。

步骤36：将和小波稀疏变换正交矩阵的逆矩阵做逆小波变换得到重构的时域信号

本发明在上述重构步骤中，将迭代次数设为固定的稀疏度K的两倍，减少了为了求解迭代次数的时间，并且能平衡了信号精确度和重构时间之间的矛盾，简化了正交匹配追踪算法（OMP）。

实施例2

下面以实际的绝缘子污闪泄漏电流信号为例，说明本发明的实施过程。本发明采用的是在盐度(ESDD)为0.4、灰度(NSDD)为2，电压为110KV条件下用10kHz的采样频率对标准型绝缘子表面的污闪泄漏电流信号，采集其中的2048个点进行案例分析，采样数据集为：

[0.863;0.924;1.006;1.22;1.644;2.234;2.854;3.369;3.721;4.01;4.49;...;-0.1]

单位为毫安（mA），用MATLAB画出的信号图如图2所示。

步骤1：选择正交小波基为稀疏基Ψ，构造小波稀疏变换正交矩阵并求出原始泄漏电流信号的稀疏度K。

本发明选择sym8正交小波基，且当将长度为N的原始泄露电流信号f，用一组基向量Ψ=[ψ₁,ψ₂,..,ψ_N]的线性组合表示，即Ψ是N×N矩阵，x是信号f在Ψ域上的等价表示，当x仅有K（K<<N）个非零系数时，则信号是稀疏的，此时称Ψ为信号f的稀疏基，而此时的稀疏系数个数成为稀疏度K。在本例中N=2048。

步骤11：选择sym8正交小波基作为稀疏基，根据sym8正交小波基分别设计低通滤波器和高通滤波器。

步骤12：分别使用低通滤波器和高通滤波器对原始泄露电流信号进行从1到log₂(N)层的小波分解，得到各层的低通滤波器分解系数h和高通滤波器分解系数g，其中N为原始泄露电流信号的长度，log2（N）=log2（2048）=11。

保留小数点后四位后分别得到：

h=[-0.0034,-0.0005,0.0317,0.0076,-0.1433,-0.0613,0.4814,0.7772,0.3644,-0.0519.-0.0272,0.0491,0.0038,-0.0150.-0.0003,0.0019]，

g=[-0.0019,-0.0003,0.0150,0.0038,-0.0491,-0.0272,0.0519,0.3644,-0.7772,0.4814,0.0613,-0.1433,-0.0076,0.0317,0.0005,-0.0034]。

其中sym8小波基为具有8阶消失矩的symlet小波，由尺度方程$>\phi \left(x\right)=\sqrt{2}\underset{n\in Z}{\Sigma }h\left(n\right)\phi (2x-n),$>小波基构造方程$>\Psi \left(x\right)=\sqrt{2}\underset{n\in Z}{\Sigma }{g}_n\phi (2x-n)$>以及关系式g(n)=±(-1)^1-nh(1-n)可以得出小波基函数Ψ(x)。h(n)可以由下面的方程构造：

$>\left(\begin{array}{c}W\left(z\right)=U\left(z\right)\overline{U(\frac{1}{z})}={\left[{\cos }^2\left(\frac{w}{2}\right)\right]}^{N}P\left[{\sin }^2\left(\frac{w}{2}\right)\right]\\ P\left(y\right)=\underset{k=0}{\overset{N-1}{\Sigma }}{C}_k^{N-1+K}{y}^k\end{array}\right)$>

其中z=e^iw，ω是角频率。当U(z)的根的模的平方大于1时构造的是symlet小波。

步骤13：构造小波稀疏变换正交矩阵

子步骤1301：令j=log₂L+1，L为低通滤波器分解系数h的长度，L=16，令初始正交矩阵为单位矩阵I。

子步骤1302：利用第j层的低通滤波器分解系数提取原始泄露电流信号的低频特征系数，得到低频特征系数向量H_j；利用第j层的高通滤波器分解系数提取原始泄露电流信号的高频特征系数，得到高频特征系数向量G_j。

子步骤1303：令i=1。

子步骤1304：对低频特征系数向量H_j和高频特征系数向量G_j进行周期性延拓，即对低频特征系数向量H_j和高频特征系数向量G_j分别增加2^j-L列零向量，得到低频特征系数矩阵H_j′和高频特征系数矩阵G_j′。

子步骤1305：分别对低频特征系数矩阵H_j′和高频特征系数矩阵G_j′中的列向量进行单步圆周位移。

子步骤1306：判断i>2^j-1是否成立，如果i>2^j-1则执行步骤1307；否则，令i=i+1，返回步骤1304。

子步骤1307：对低频特征系数矩阵H_j′和高频特征系数矩阵G_j′进行正交化处理，得到正交矩阵

子步骤1308：判断j>log₂N是否成立，如果j>log₂N，则执行子步骤1309；否则令且令j=j+1，返回子步骤1302。

子步骤1309：此时结束j的迭代，得到最终的正交矩阵并令即为小波稀疏变换正交矩阵。

因为矩阵太过庞大，故不在此列出，可以使用某些编程工具计算，比如MATLAB，结果可以直接代入下一步骤。

步骤14：根据公式计算原始泄漏电流信号的小波稀疏系数向量，其中x为原始泄漏电流信号在sym8正交小波基上的2048×1维信号，则XWAVE为2048×1维：

XWAVE=[0.0177;0.1378;0.1823;0.2183;0.0045;0.2838;0.1165;0.0354;0.0243;0.1461;…;6,2695e-4;3.2671e-5;4.6706e-5;3.9453e-7;2.4082e-6]，数据太多，只在此列出XWAVE的前十个数和后五个数。

步骤15：设定小波稀疏系数的阈值为ε，将大于阈值ε的系数个数设为稀疏度K，即为原始泄漏电流信号的的稀疏度。

阈值的设定与信号的特性和变化相关，本发明为了简单计算，通常采用全局阈值的方法，确定一个统一的阈值；同时，为了使统一的阈值可以根据不同的信号特征进行适应性选择，本发明采用了划分稀疏系数的区间然后取同一区间系数平均值为阈值的方法。本例中将信号的稀疏系数XWAVE划分为20个区间，并计算各个区间的的系数平均值和系数的个数，将系数个数最大的区间的系数平均值作为阈值ε，将大于阈值ε的系数个数设为稀疏度K，本例中，将包含了1950个点的区间作为恢复系数，阈值为0.0248，稀疏度K为146。

步骤2：构造高斯随机矩阵Φ，并计算测量向量y和传感矩阵Θ。

步骤21：根据原始泄漏电流信号的稀疏度K计算测量数M，M为距离Klog(N/K)最近的正整数。

步骤22：构造阶数为M×N的高斯随机矩阵Φ，且高斯随机矩阵Φ满足与小波稀疏变换正交矩阵无关性，将信号f根据公式y=Φf进行线性投影，则信号由N维降到M维，其中y为M维测量阵，M为测量数,M≥cKlog(N/K)。本实施例计算得M=386。

步骤23：根据公式y=Φf计算测量向量y；其中，f为原始泄漏电流信号。

本例可以直接采用MATLAB中的随机矩阵函数randn(M,N)=randn(386,2048)构造386×2048高斯随机矩阵，则根据y=ΦΨx计算得到386×1的测量向量：

y=[-2.8229e+3;1.4013e+3;4.7197e+3;223.7540;5.3618e+3;4.7203+3;1.5410e+3;-1.9145e+3;-43.6646;-1.4499e+3;…;-3.1721e+3;1.3935e+3;270.0417;-3.8352e+3;72.4483]，这里只列出y的前十个数和后五个数。

步骤24：根据公式Θ=ΦΨ计算传感矩阵Θ；其中，Ψ为稀疏基矩阵。

步骤3：利用测量向量y和传感矩阵Θ重构原始泄漏电流信号。

重构原始泄露电流信号的过程又包括：

步骤31：分别设定冗余向量初始值γ₀、索引集合初始值Λ₀和初始迭代次数t。其中，冗余向量初始值γ₀=y，索引集合初始值Λ₀为空集，初始迭代次数t=0；

步骤32：根据公式计算索引值δ_t，并令Λ_t=Λ_t-1∪{δ_t}。其中，u_j为相关系数且u_j=|<γ_t,θ_j>|，即冗余向量与传感矩阵相应列向量的内积的绝对值，θ_j为传感矩阵Θ的第j列列向量，N为原始泄露电流信号的长度。

步骤33：根据公式采用最小二乘法计算使残差y-Θx_i最小时的信号估计值其中，x_i为时的原始泄露电流信号。

步骤34：根据公式计算第t次迭代的冗余向量γ_t。

本实施例中，计算得到的最终冗余向量γ_t为386×1的向量：γ_t=[-0.5379;0.0614;0.2896;0.2158;0.3833;…;1.6716;-2.2200;-2.0355;-3.1984;1.4061]。

步骤35：由于本实施例中，稀疏度K=146，因此判断t<2×146是否成立，如果t<2×146，则令t=t+1，返回步骤32；否则，执行子步骤36。

将迭代次数设为固定的稀疏度K的两倍，减少了求解迭代次数的时间，并且能平衡了信号精确度和重构时间之间的矛盾，简化了正交匹配追踪算法（OMP）。

步骤36：将和小波稀疏变换正交矩阵的逆矩阵做逆小波变换得到重构的时域信号。

最后得到的最大相关系数u_t所在的索引位置pos的信号估计则为待重构的谱域向量r_y，本例计算得到1×2048维的待重构谱域向量为：

r_y=[1.8551,142.4470,-219.7036,179.3667,0,-283.7255,…,0,50.7893,0,0,0,8.2076]。

将待重构的谱域向量和小波稀疏变换正交矩阵的逆矩阵做逆小波变换得到重构的2048×1维时域信号为：

$>\hat{f}=[22.0010;13.1960;20.5435;19.6215;10.8332;...;-6.6397;1.4932;23.0058;-29.3130;39.8138].$>

用MATLAB画出来得到的信号重构图如图3所示,可以看出达到不错的重构效果。

本发明提出的方法能够在较低采样率下，对污闪泄漏电流信号进行稀疏表示，达到迅速而准确的采样和重构，为进行污闪泄漏电流的实时监测的实现打下了很好的基础。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。