一种管道多点泄漏精确定位方法与流程

本发明涉及压力管道泄漏检测技术领域，更为具体来说，特别涉及一种基于改进ceemd和盲源分离的管道多点泄漏精确定位方法。

背景技术：

城市管道在长期使用过程中，因受内外腐蚀等多种因素影响而容易引发泄漏。尤其是早期发生的微小泄漏，由于其泄漏孔径小、泄漏流量小等特点，导致微小泄漏难于通过视觉、听觉和管道检测系统捕捉。实际工况中，管道的泄漏通常是多个泄漏点相互作用形成的，而多点泄漏往往具有信号干扰大，信号包含内容多等特点，导致多点泄漏信号难以提取，定位误差大，给实际检测带来诸多问题。

经验模态分解(emd)是一种针对非线性、非平稳信号的分析方法，这种算法将信号分解为一系列振荡函数,其强大的自适应特征使得emd广泛应用于信号分析领域。然而，emd本身存在模态混叠，其导致分解出的具有不同振幅特征的固有模态函数(imf)相互叠加；为此，在其基础上产生了集合经验模态分解(eemd)，它将白噪声加入待分析信号再进行emd分解，一定程度上缓解了模态混叠效应，但是eemd本身也存在一些问题，如信号分解过程中白噪声没有完全被消除等。针对上述问题，本文采用互补集合经验模态分解(ceemd)对泄漏信号进行分解，其有效缓解了eemd本身的缺陷，也大大提高了计算精度。但ceemd也并非完美无瑕,由于利用大量的总体集成平均去除添加噪声，因此其运行时间与传统的emd方法相比要长很多。

盲源分离技术(bbs)是目前管道故障诊断领域分离混杂信号的常用方法，管道故障诊断的实质就是在泄漏源信号的数目和源信号特征未知的情况下，将泄漏信号单独地分离出来。因此，将盲源分离技术引入管道故障诊断是十分适合的，且现已存在相当多的成功应用案例。然而，现有的盲源分离技术仍然存在一定的局限性：一方面，目前的盲源分离算法大多针对超定或者正定情况，即要求观测信号数目大于或等于源信号数目，对欠定情况(观测信号数目小于源信号数目)并不满足，从而造成分离出的源信号仍然存在一定的混叠，影响最终诊断结果。另一方面，真实的管道泄漏信号具有非平稳、非线性等特点，现有方法对处理这种信号的表现通常较差。

针对上述问题，本发明提出一种管道多点泄漏精确定位方法，旨在：①通过改进的ceemd方法减少分解原始信号所用时间，提高管道泄漏检测的工作效率。②通过盲源分离与信号稀疏性结合的方法解决实际应用中欠定情况下的管道故障诊断，提高管道泄漏检测的精度。

技术实现要素：

为克服现有技术的不足，本发明提供了一种管道多点泄漏精确定位方法，以实现对两点泄漏源的分离和精确定位。以此方法应用于管道多点泄漏源定位，进而实现泄漏多点源的精确定位。该定位方法包括步骤如下：

一种管道多点泄漏精确定位方法，包括以下步骤：

s1：在被测管道上下游两端分别设置一个声发射传感器，且两声发射传感器之间的距离为l，利用两端的声发射传感器采集管道泄漏信号，作为管道原始泄漏信号x(t)；

s2：利用改进的ceemd对管道原始泄漏信号x(t)进行分解处理，将管道原始泄漏信号x(t)分解为各阶imf[imf1,imf2,...,imfn]。其具体过程如下：

原始的ceemd具有较小重构误差，但需要通过多次集成平均去除加入噪声的影响，使得运行速度很慢。在使用原始的ceemd分解信号时，噪声信号、高频冲击信号会首先被分解出来，作为前几阶imf分量。当这些会引起模态混叠的分量被分离后，剩余的信号无需进行完整的ceemd分解过程。

因此，提出一种改进的ceemd方法，将原始的ceemd方法与峭度理论及自相关函数相结合。峭度值对冲击信号十分敏感，使用峭度值检验各阶imf中是否含有大的峰值，从而判定其是否为冲击信号；噪声信号的自相关函数具有快速衰减到零的特点，因此使用自相关函数极大值点的均方根值来判定各阶imf是否为噪声信号。

当冲击信号、噪声信号从原始信号中分离后，改用emd方法对剩余信号进行分解，旨在节省ceemd的分解时间，提高ceemd的计算效率。其具体操作步骤如下：

(1)在管道原始泄漏信号x(t)中加入一对符号相反的噪声，添加噪声的对数为n，于是得到2n个添加噪声后的信号，如下式：

mi⁺(t)＝x(t)+ni⁺(t)，(i＝1,2,3......n)

mi^-(t)＝x(t)+ni^-(t)，(i＝1,2,3......n)(1)

式中，mi(t)为第i个添加噪声后的信号，x(t)为原始信号，ni(t)为添加的第i个噪声信号，“+”和“-”表示符号相反。

(2)利用emd算法中筛分过程，对步骤(1)得到的2n个添加噪声后的信号mi⁺(t)，mi^-(t)进行分解，得到2n个一阶imf分量imfi,⁺1和imfi,^-1。原始泄漏信号x(t)第一个imf分量imf1为2n个一阶imf分量的集成平均值，即：

之后对剩余函数ri⁺(t)及ri^-(t)进行更新：

ri⁺(t)＝mi⁺(t)-imfi,⁺1

ri^-(t)＝mi^-(t)-imfi,^-1(3)

(3)将剩余函数ri⁺(t)以及ri^-(t)作为新的分解对象重复步骤(2)，得到第二个imf分量imf2。当第p个imf分量imfp既不是冲击性分量也不是噪声分量时，停止迭代。

(4)从原始泄漏信号x(t)中减去前三个步骤中得到的p-1个imf分量后，改用emd方法分解剩余信号，最终可以得到：

式中，imfi为步骤(1)～(3)获得的imf分量；imf′i为步骤(4)获得的imf分量，rn(t)为步骤(4)获得的余量。

s3：对步骤s2分解得到的各阶imf分别求其包络信号c(t)＝[c1,c2,...,cn]。

包络信号主要与故障信息有关，对各阶imf求其包络，可使各阶imf的物理意义更加明确，有利于噪声的进一步消除。传统的信号包络提取是通过翻转和低通滤波来实现的，希尔伯特变换(hilberttransform,ht)提供了一个非常好的替代选择。原始泄漏信号x(t)的希尔伯特变换定义为：

具体地，对ceemd分解之后的imf分量其形式为：

对应的包络信号为：

其中，ci表示第i个imf分量对应的包络信号，i＝1,2,3......n，将所有imf分量的包络信号算出后，则c(t)＝[c1,c2,...,cn]。

s4：对步骤s3得到的包络信号c(t)进行样本熵计算，并根据所求样本熵值的大小对各阶imf进行优选，从而得到观测信号h(t)。

样本熵是在近似熵的基础上发展来的，近似熵是时间序列的复杂度的一种表示，利用这一特点，可将其作为判断故障的特征参数，但近似熵存在自匹配的问题，且其计算严重依赖于数据的长度。为了克服这一缺点，研究人员对近似熵算法进行改进，从而形成样本熵。样本熵不仅减少了对序列长度的依赖，同时克服了近似熵自身匹配的缺点。

对原始泄漏信号序列x(t)＝{x(1),x(2),...,x(n)}，其样本熵的计算过程如下：

(1)由原始泄漏信号序列构成一个m维向量：

xm(i)＝{x(i),x(i+1),...,x(i+m-1)},1≤i≤n-m+1(8)

(2)定义xm(i)和xm(j)之间的距离：

(3)给定容限r，统计每一个xm(i)对应的d[xm(i),xm(j)]≤r的数目，记为ai。把ai与n-m+1的比值记为：

(4)求出bi^m(r)的平均值b^m(r)：

(5)用同样的方法求出b^m+1(r)，则该时间序列的样本熵定义为：

已知ceemd分解后的各阶imf是由噪声信号分量、有效信号分量和冗余分量组成的，且求其包络信号并未影响其组成结构。因此，通过分别计算其包络信号样本熵实现imf的优选。样本熵值越大，对应其非稳定时间序列的复杂程度越高，即片段中包含的泄漏信号越多；而噪声分量的包络信号表现为周期震荡函数，其复杂程度并不高，即样本熵值较低。因此，最终的观测信号为优选的样本熵值较大的imf之和，即：

h(t)＝imfa+imfb+...+imfn(13)

式中，imfa，imfb，imfn为优选的imf。

s5：将s4得到的观测信号重构、升维，利用盲源分离和信号稀疏性表示相结合的方法，得到两点泄漏分离后的有效泄漏信号y(t)；

经典的盲源分离算法大多针对超定或者正定情况(观测信号数目大于或等于源信号数目)，而实际工程应用中更多的为欠定情况(观测信号数目小于源信号数目)，经典盲源分离并不适用，从而造成分离出的源信号仍然存在混叠，影响最终诊断结果。本发明针对实际欠定情况，将盲源分离与信号稀疏性原理进行结合以获得多个相互独立的源信号。

首先将升维后的观测信号h(t)作为盲源分离的输入分量，则：

h(t)＝as(t)+n(t)(14)

式中，h(t)＝[h1(t),h2(t),...,hm(t)]^t是m维观测信号，s(t)＝[s1(t),s2(t),...,sn(t)]^t是n维源信号，a是一个m×n维混合矩阵，其中m＜n，n(t)＝[n1(t),n2(t),...,nm(t)]^t是m维附加噪声信号。

不考虑噪声影响的，上式可写为如下矩阵形式：

改写为向量形式，即：

h(t)＝a1s1(t)+a2s2(t)+...+ansn(t)(16)

由于源信号不具有稀疏性或稀疏性不强，通过傅里叶变换将其变换到时域中，使其在时域中具有稀疏性。满足稀疏特性后，将只有某个源信号的取值很大，其他源信号取值很小或者接近于零。假设对于某一采样时刻t，其他源信号幅值皆很小或为零，只有源信号si(t)取值非零，则式上式可以写成：

h(t)＝aisi(t),i＝1,2,...,n(17)

通过以上步骤，每次可以从观测信号中分离出一个源信号，经过多次分离可将所需有效泄漏信号从混合信号中分离出来。

s6：对s5得到的有效泄漏信号进行互相关计算，确定泄漏信号的时差δt；

互相关时延计算是计算两信号相关性的常用方法，利用互相关函数可对两路有效泄漏信号到达上、下游传感器的时差进行计算，其公式为：

式中，y1和y2为两路有效泄漏信号，m和n均为整数。

于是，两路有效泄漏信号的时差δt可通过下式计算：

δt＝±m×g(19)

式中，m为得到的互相关函数中最大相关峰所对应的m值，g为信号的采样周期，正、负号表示y(1)早于或晚于y(2)。

s7：通过实地测量和查阅工程手册确定泄漏信号在管道中的传播速度v；

s8：针对金属管道两点泄漏的定位采用时差定位公式进行计算；

式中，x为被检测管道泄漏源位置，即泄漏点到上游声发射传感器的距离；l为两声发射传感器之间的距离，δt为泄漏信号到达两传感器的时间差，v为泄漏信号在管道中的传播速度；由计算得到的泄漏信号到达两传感器的时间差δt和泄漏信号在管道中的传播速度v，结合时差定位公式最终计算得到的泄漏定位结果。

本发明可以达到以下有益效果：

1、利用改进的ceemd分解，可以大幅度减少运算时间，从而提高管道泄漏检测效率。

2、利用样本熵分析，可同时分析“低频成分”和“高频成分”的熵信息，避免了遗漏在“高频成分”里的信息，能够提取更加丰富的故障信息。

3、利用盲源分离和信号稀疏性表示相结合的方法，可解决更符合实际应用的信号欠定混叠问题，更准确地获得非平稳振动信号的时频变化特征，实现真实信号的特征提取。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1是本发明最佳实施例的流程示意图；

图2是实验管道多点泄漏系统图；

图3是上游ceemd处理结果图；

图4是上游各阶imf对应的包络信号图；

图5是上游各阶imf样本熵值结果汇总图；

图6是上游观测信号图；

图7是下游观测信号图；

图8是上游有效泄漏信号图；

图9是下游有效泄漏信号图；

图10是互相关时延计算结果图。

具体实施方式

现在结合附图对本发明作详细的说明。此图为简化的示意图，仅以示意方式说明本发明的基本结构，因此其仅显示与本发明有关的构成。

如图1所示，本发明的一种管道多点泄漏精确定位方法，包括以下步骤：

s1：在被测管道两端分别设置一个声发射传感器，且两声发射传感器之间的距离为l，利用两端的声发射传感器采集原始泄漏信号x(t)；

实验管道如图2，管道长度为50m，管道材质为钢材，管道规格为dn150，介质为压缩空气，管道内介质处于流动状态；上游声发射传感器1放置于零点位置，距离上游声发射传感器16m和30m处分别为两个泄漏孔径为2.0mm的泄漏孔；距离上游声发射传感器40m的地方安装固定下游声发射传感器2，即两传感器间距40m。利用声发射传感器对泄漏信号进行采集，管道压力0.3mpa，实验随机提取30个泄漏信号数据作为后续处理的原始泄漏信号x(t)，如表1所示。

表10.3mpa压力下金属管道两点泄漏原始数据

从表中可以看出，泄漏管道声发射粗定位存在较大的误差，需要进一步分析处理；本实施例随机取第1个数据进行具体计算。

s2：利用改进的ceemd对原始泄漏信号x(t)进行分解处理，将原始泄漏信号x(t)分解为各阶imf[imf1,imf2,...,imfn]；

借助matlab编制相应代码，对原始泄漏信号x(t)进行ceemd处理，将原始泄漏信号x(t)分解为各阶imf的形式，其结果如图3所示，上游原始泄漏信号被分解为8阶imf分量。

s3：对s2分解得到的各阶imf分别求其包络信号c(i)＝[c1,c2,...,cn]；

借助matlab编制相应代码及公式，对各阶imf进行希尔伯特变换，其计算结果如图4所示，上游各阶imf分别对应变换为包络信号的形式。

s4：对s3得到的包络信号c(i)进行样本熵计算，并根据所求样本熵值的大小对各阶imf进行优选，从而得到观测信号h(t)；

借助matlab编制相应代码及公式，计算上游原始泄漏信号1～8阶imf对应的包络信号的样本熵，并进行汇总统计，计算结果如图5所示。根据样本熵值大小对各阶imf进行优选：样本熵值越大，对应其非稳定时间序列的复杂程度越高，即片段中包含的泄漏信号越多；而噪声分量的包络信号表现为周期震荡函数，其复杂程度并不高，即样本熵值较低。从图5中可清晰看出，前5阶imf对应的包络信号有明显的样本熵值，其中第1阶的样本熵值最大，而6～8阶的样本熵值几乎为零，可将其忽略不计。因此选择前6阶imf进行重构，最终获得上游观测信号如图6所示。同样地，对下游原始泄漏信号进行s2～s4步骤所述的计算处理，得到下游观测信号如图7所示。

s5：将s4得到的观测信号h(t)进行重构、升维，利用盲源分离和信号稀疏性表示相结合的方法，得到两点泄漏分离后的有效泄漏信号y(t)；

利用matlab中的盲源分离工具箱分别对上、下游观测信号h(t)进行时域分析，得到两点泄漏分离后的有效泄漏信号y(t)，如图8和图9所示。

s6：对s5得到的有效泄漏信号y(t)进行互相关计算，确定泄漏信号的时差δt；

借助matlab编制相应代码，对有效泄漏信号y(t)进行互相关时延计算，确定有效泄漏信号到达上、下游传感器的时间差δt，如图10所示。从图中可以看出，两路信号的互相关函数rxy在横坐标取值-7.36时达到最大值，因此两路信号的时延值为7.36毫秒，即有效泄漏信号到达上、下游传感器的时间差δt＝0.00736s。

s7：通过实地测量和查阅工程手册确定泄漏信号在管道中的传播速度v；

管道泄漏声速v不仅受管道材料的影响，还受到不同介质、工况的影响，且使用不同方法检测到的波速都有所差别。目前，国内外尚未形成统一的泄漏声波速度计算方法，各研究人员检测计算得到的波速值也不一致。根据沈功田《声发射检测技术及应用》(科学出版社，2015年版)：压缩空气介质所产生的声发射波在钢制管道中的传播波速为850～1050m/s，并通过大量实验对典型管道泄漏声波速度进行了测量，其中，对65m长的φ159×60mm钢制管道的测量恰好与本实施例工况一致，最终确定其声速平均值为950m/s。

s8：根据两声发射传感器之间的距离为l、泄漏信号到达两传感器的时差δt、泄漏信号在管道中的传播速度v，利用时差定位公式对泄漏点进行精确定位；

式中，x为被检测管道泄漏源位置，即泄漏点到上游声发射传感器的距离；l为两声发射传感器之间的距离，δt为泄漏信号到达两个传感器的时间差，v为泄漏信号在管道中的传播速度。

将s6，s7步骤得到的时间差和声速代入上式即可计算得到泄漏定位结果，即x＝(40-0.00736×950)/2＝16.50m。

将本发明信号处理定位结果与未处理粗定位信号计算结果对比，得到表2。

表2信号处理结果前后对比

由以上计算过程可知，粗定位的相对误差为23.25％，经过本发明处理方法处理后，泄漏定位的相对误差降低至3.13％。可见，本发明方法可使得泄漏源定位的精确度得到有效提高。

另外，对于同一信号，先后采用本发明中改进的ceemd与传统的ceemd进行分解计算，分别计算10次，然后求其平均时间并进行对比，得到表3。

表3两种ceemd运行时间对比

由此可见，改进的ceemd可大幅度减少运算时间，从而提高管道泄漏检测效率。

以上述依据本发明的理想实施例为启示，通过上述的说明内容，相关的工作人员完全可以在不偏离本发明的范围内，进行多样的变更以及修改。本项发明的技术范围并不局限于说明书上的内容，必须要根据权利要求范围来确定其技术性范围。