CN110454687B - 一种基于改进vmd的管道多点泄漏定位方法 - Google Patents

Abstract

本发明涉及一种基于改进VMD的管道多点泄漏定位方法，包括：采集管道原始泄漏信号；对所述原始泄漏信号进行总体局域均值分解，得到若干PF分量；计算各PF分量的相关系数，根据所述相关系数筛选出所需PF分量，并根据筛选出来的PF分量进行信号重构，确定变分模态分解的k值；对重构信号进行变分模态分解，得到若干IMF分量，计算各IMF分量的多尺度熵值，并根据各IMF分量的多尺度熵值筛选IMF分量；对筛选出来的IMF分量进行信号重构，通过对盲源分离后的各泄漏信号进行互相关定位计算完成管道泄漏定位。本发明能够对管道泄漏信号进行有效提取，消除了原始泄漏信号中的低相关成分和噪声的影响，使得最终定位的结果更精确。

Description

一种基于改进VMD的管道多点泄漏定位方法

技术领域

本发明涉及一种管道泄漏定位方法，尤其涉及一种基于改进VMD的管道多点泄漏定位方法。

背景技术

城市管道已成为现代城市发展不可或缺的工具，随着其规模的不断扩大，由于设备自然老化、气候环境以及人为破坏等影响，管道故障事件呈上升趋势，尤其是燃气管道一旦泄漏，很容易引起火灾、爆炸、中毒、环境污染等恶性事故。因此如何有效识别管道初期的微小泄漏，准确找到泄漏点位置，具有良好的经济价值和社会意义。

近年来，随着计算机技术的发展，管道泄漏检测技术正向软硬件结合的方向发展，各种新的管道泄漏检测方法仍是各国的前沿研究方向，而非平稳信号的分析和研究一直是研究的热点，如小波分解、经验模态分解(EMD)、局域均值分解(LMD)等，但它们也都有各自的不足。比如，小波分解不能进行自适应选择，EMD、LMD是自适应的时频分析方法，但却都存在模态混叠、过包络和欠包络等问题。变分模态分解(Variational Mode Decomposition，VMD)是Dragomiretskiy等人于2014年提出的一种新型自适应故障诊断方法。与EMD、LMD方法相比，VMD方法具有理论基础牢固、收敛速度快、分解结果的特征明显等优点。VMD作为一种分解算法，与EMD、LMD方法相似，可以将故障信号根据高低频分解为若干本征模态函数(Intrinsic Mode Functions，IMF)。但是VMD在分解前，需要确定分解层数k。分解层数如果过多，分解后的分量会出现断层现象，分解层数过少的话分解后的分量的泄漏信息变少，影响定位结果。但目前有关VMD相关研究中，如：《基于VMD分量相对熵分析的管道泄漏定位方法》，其分解层数K的确定多凭人为经验确定，可能导致分解不准确、定位误差过大等问题。所以在使用VMD方法前，需要对信号进行一定的降噪处理，来选取最优的参数K。

本发明正是在这种背景下产生，在对原始泄漏信号进行变分模态分解(VMD)之前通过分析总体局域均值分解(ELMD)的结果来优化选取VMD分解的个数，确定最优参数，分解过程中添加白噪声减少了重构误差，将相关系数较高的分量保留，提高了分解效率，通过基于独立分量分析的盲源分离方法进行信号估计并与原始泄漏信号进行互相关计算，最终完成管道多点泄漏的精确定位。

发明内容

本发明的目的在于克服现有管道多点泄漏定位方法中信号特征提取困难、定位不准确等问题，提供了一种基于改进VMD的管道多点泄漏定位方法，通过计算ELMD分解结果各分量的相关系数，分析选取VMD方法中的分解层数，更为有效地提取泄漏信号，提高泄漏定位的准确性，增强方法的实用性和可靠性。

本发明所采取的技术方案为：一种基于改进VMD的管道泄漏定位方法，该方法的步骤如下，

S1：在被测管道上设置声发射传感器，采集管道原始泄漏信号；

S2:对原始泄漏信号进行总体局域均值分解(ELMD)，得到一系列乘积函数(PF分量)；

S3：计算各PF分量的相关系数，筛选出有效且包含原信号信息较多的PF分量，并根据筛选出来的PF分量进行信号重构，并根据PF分量的相关系数计算分析，确定VMD分解的最优k值；

S4：对重构信号进行变分模态分解(VMD)，得到一系列IMF分量，通过计算各IMF分量的多尺度熵(MSE)，筛选IMF分量；

S5：对筛选出来的IMF分量进行信号重构，得到观测信号，并利用独立分量分析的盲源分离方法分离出各泄漏点的信号(估计泄漏信号是盲源分离后的估计信号，采用盲源分离方法对原始泄漏信号进行分离，通过计算估计泄漏信号与原始信号的相关系数来区分是哪个泄漏点的信号，区分之后得到各泄漏点的信号，)，通过互相关定位算法进行管道泄漏定位。

进一步地，步骤S2中对原始泄漏信号x(t)进行ELMD分解，具体过程为：

S2.1：确定总体平均次数M及加入的白噪声等级，在信号x(t)中加入确定等级的白噪声n_m(t)，则泄漏信号表示为：

x_m(t)＝x(t)+n_m(t) (1)

其中，m为分解次数，t是时间，n_m(t)为白噪声信号，x_m(t)为加入白噪声后的混合泄漏信号；

S2.2：对混合泄漏信号x_m(t)进行LMD分解，再做总体平均得到最终PF分离，用公式(2)表示为：

其中，ε_n，m(n＝1，2，...，.N)，N为正整数；ε_n，m为第m次分解得到的第n个初始PF分量，将N个初始PF分量M次分解的均值y_n(n＝1，2，...，N)作为最终PF分量，作为后面步骤S3所使用的PF分量。

进一步地，步骤S3中分析ELMD分解结果并确定VMD参数k值的具体过程为：

传统方法如专利《基于VMD分量相对熵分析的管道泄漏定位方法》中，采用观察中心频率的方法确定K。通过观察不同K值对应的各个模态中心频率，当出现了中心频率相近的两个模态分量时，就认定上一层的模态层数为K的取值。这种方法通过人为观察模态中心频率来选取K值，会造成有效泄漏信号分解、提取结果不够准确，以致带来VMD分解后重构误差较大。因此本发明采用皮尔逊积矩相关系数计算各个PF分量y_n同原始信号x(t)的相关系数r。相关系数r的取值范围为[-1，1]，相关系数的绝对值越大，则两个样本的线性相关程度越高。峭度是无量纲参数，用来描述信号波形的尖峰度，对冲击信号较为敏感。泄漏信号中主要成分为冲击成分，其它成分则认为是干扰噪声信号，而峭度指标值为3，即认定峭度值小于3的信号为干扰信号。而在管道泄漏中，峭度值越大，说明信号的泄漏信息越多。因此，通过求各个PF分量同原始信号之间的相关系数，由相关系数定义可知，相关系数绝对值大于0.5的即为强相关分量，因此选取相关系数绝对值大于0.5的PF分量，对其计算峭度值，筛选峭度值小于3的PF分量，以此降低采集管道泄漏信号时非相关量的干扰；分析筛选出来的PF分量，确定VMD分解中的分量个数(k的取值)；然后对筛选后的PF分量进行信号重构，得到重构信号x′(t)。与王秀芳等人的专利《一种天然气管道泄漏定位方法和装置》不同，本发明通过对各PF分量求相关系数，根据相关系数筛选绝对值大于0.5的PF分量，然后通过计算分量的峭度值来判断是否含有噪声信号，以此确定VMD分解层数K的取值,K是分量个数，也是分解层数.

进一步地，步骤S4中对重构信号x′(t)进行VMD分解并筛选IMF分量的具体过程为：

S4.1：对重构信号x′(t)，利用VMD分解成k个固有模态函数分量，如下式所示：

其中，u_k为第k个固有模态函数分量，k＝1，2，...，n

而VMD算法实质是对变分问题的求解，其过程如下：

对于每一个IMF分量u_k(t)，采用Hilbert变换可以得到其解析信号u′_k(t)为：

其中δ(t)是狄里克莱函数，j是虚数符号，对各模态解析信号预估一个中心频率ω_k并进行混合，将每个模态的频谱调制到相应的基频带，构建新的解析信号，记为u″_k(t)如下所示：

计算以上解调信号的梯度的平方L²范数，估计出各模态分量的带宽。建立的约束变分模型Y为：

其中，{u_k}＝{u₁，u₂，…u_k}表示分解得到的k个IMF分量，{ω_k}＝[ω₁，ω₂，…ω_k}表示各分量的中心频率，s.t.是约束条件的英文缩写。

为求解上述约束变分模型Y，引入二次惩罚因子α和拉格朗日乘子λ(t),构造扩展拉格朗日表达式如下：

利用乘子交替方向算法不断更新各IMF及其中心频率：

λⁿ⁺¹，其中，u_k是分解得到的k个IMF分量，ω_k是u_k相对应的中心频率、λ是拉格朗日乘子，n是正整数。直到满足以下迭代停止条件：

其中，判别精度θ＞0，最终求得式(7)的鞍点即为最优解，得到k个{u_k}。

S4.2：，通过计算多尺度熵对k个IMF分量进行筛选，多尺度熵值越大，表明其分量包含的泄漏信息越多。设IMF分量u_k的原始时间序列为{x₁，x₂，...，x_n}，长度为N，设定其嵌入维数w，相似容限z，对其进行粗粒化变换，得到新的时间序列：

其中，k＝1，2，...，s/τ，μ，k为正整数且μ≠k，s为离散时间序列长度，τ为尺度因子，原始序列被分割成τ段且每段长为s/τ的粗粒序列，当τ＝1时，新的时间序列就是原始序列；通过维数w和相似容限z，构建时间序列的w维向量，记为x_μ(w)＝{x_μ，x_μ+1，…x_μ+w-1}，其中x_μ为时间序列中某一维度的向量，μ＝(k-1)τ+1，x_μ是原始序列里的值。

对每个长度为s/τ的粗粒时间序列分别求样本熵，得到其样本熵值，并绘制成尺度因子的函数，按公式(10)得到多尺度熵值：

其中，MSE为多尺度熵，SampE为样本熵，w是嵌入维数，k是正整数，τ是尺度因子，z是相似容限，s是时间序列长度，u_k是IMF分量。

S4.5：选取在各尺度下最优的多尺度熵MSE_ｎ，最终计算得到各IMF分量的多尺度熵值记为：MSE₁，MSE₂，MSE₃…MSE_n。与《基于PSO-VMD算法的管道泄漏检测装置及检测方法》专利相比，本发明通过ELMD对VMD分解层数K值进行了优化，并且根据多尺度熵值对分解后的IMF分量进行各个维度上的分析，相比信息熵更能提取到有效泄漏信息。而为了进行有效性分析，再和基于单一VMD分解的方法对比，进行验证。

进一步地，步骤S5中通过基于独立分量分析的盲源分离得到估计的泄漏信号并与原始泄漏信号进行互相关定位的具体过程为：

S5.1：通过上述过程得到各IMF分量的多尺度熵值MSE_ｎ，进行筛选。每个IMF分量表示多个尺度下含有不同程度的泄漏信息的模态分量，提取筛选后每个IMF分量的时域、频谱进行重构，得到观测信号Nimf(t)＝[b₁(t)，b₂(t)…b_n(t)]，该观测信号是信号的数学表现形式，与时间相关的函数，通过MATLAB对IMF分量的时域、频域进行重构得到；构成的多维观测信号很好解决了盲源欠定问题，即当传感器数目小于信号源数目；

S5.2：通过对观测信号Nimf(t)进行分离，计算观测信号的相关矩阵R_x＝E[Nimf(t)Nimf^H(t)]，根据特征值估计泄漏信号数目μ。并对观测信号Nimf(t中心化，并作白化处理，最终得到盲源分离后管道的μ个估计信号Y_μ(t)；再对估计信号和原始泄漏信号进行相关系数计算，确定估计泄漏信号分别与哪个原始泄漏信号一一对应；

S5.3：通过互相关计算上下游盲源分离后的估计泄漏信号Y_μ(t)和Y′_μ(t)的时间延迟Δt，在压力、孔径和介质等参数确定的情况下，得到传播速度v，依据时差定位法即可确定管道泄漏点的位置：

式中，x为估计泄漏定位值，即泄漏点到上游压力传感器的距离，l为上下游的传感器距离，v是泄漏声发射信号在管道中的传播速度；Δt是时间延迟；Y_μ(t)和Y′_μ(t)分别是盲源分离后上游得到的估计泄漏信号和下游得到的估计泄漏信号。

本发明的有益效果是：本发明提供一种基于改进VMD的管道多点泄漏定位方法，该方法通过声发射仪采集管道泄漏信号，利用总体局域均值分解和变分模态分解信号进行降噪处理。通过相关系数分析ELMD分解后的结果，优化选取VMD方法中的分解层数K，再通过计算IMF分量的多尺度熵最终筛选出含有大量泄漏信息的IMF分量，通过基于独立分量分析的盲源分离对IMF分量进行重构、降噪，最后用互相关对盲源分离后得到的估计泄漏信号进行定位。本发明能够对管道泄漏信号进行有效提取，消除了原始泄漏信号中的低相关成分和噪声的影响，使得最终定位的结果更精确。

附图说明

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

图2是管道泄漏定位布置示意图；

图3是上游原始泄漏信号的时频图；

图4是ELMD分解结果图；

图5是VMD分解结果图；

图6是盲源分离后的估计泄漏信号时域图。

具体实施方式

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

本发明的一种基于改进VMD的管道多点泄漏定位方法，如图1为本发明的具体流程图。结合了总体局域均值分解算法、变分模态分解算法、多尺度熵和基于独立分量分析的盲源分离等方法，具体步骤如下：

1、实验参数如下：管道长50m，管道初始压力为0.3MPa，两个传感器#1#2间距为42m，分别放置于管道上游与下游，泄漏点L1、L2、L3分别位于距上游传感器8m，18m，30m处，泄漏孔径为1mm。实验管道布置示意如图2所示。声发射仪采用美国PAC公司生产的6006型4通道数字式检测仪，声发射传感器采用低频传感器，型号为DP3I，响应频率为20～100KHz，谐振频率为30KHz，采样率为1Msps,内置放大器增益为40dB。利用声发射仪采集管道的泄漏信号，得到的上下游泄漏信号分别记为x₁(t)、x₂(t)。以上游传感器采集得到的原始泄漏声发射信号x₁(t)为例，图3是x₁(t)的时域、频域分布。

2、对泄漏信号x₁(t)进行总体局域均值分解(ELMD)，根据经验取平均次数为100次，以减弱添加的噪声对真实PF的影响；而若添加白噪声的幅值过小起不到改变极值点分布的作用；如果幅值太大固然能减弱所添加噪声影响，但也会使分解过慢。因此设置添加的白噪声的幅值为原信号的标准差的0.15倍；对所有残留噪声的各阶PF分量分别做总体平均运算，最终得到ELMD分解结果如图4所示。

从图4中可以看出，原始泄漏信号x₁(t)经过ELMD分解得到了8个PF分量和1个残余分量u_i，并且可以看出其分量受到噪声干扰较大，不能很好地提取出原始泄漏信号的波形，而且PF₇—PF₈分量幅值有较大的变化，因此ELMD方法具有一定的消噪能力，但是在分解时会发生模态混叠现象。

3、根据相关系数准则计算各PF分量的相关系数进行PF分量的筛选。相关系数是反映变量之间的相关关系密切程度。采用皮尔逊积矩相关系数计算各个PF分量同原始信号x(t)的相关系数r和峭度值。如表1所示：

表1各PF分量的相关系数和峭度值

根据表1得到的结果，PF₁-PF₇分量的相关系数值大于0.5，而PF₈相关系数小于0.5，由相关系数定义可知，相关系数大于0.5即为强相关分量，因此，舍弃PF₈分量。然后PF₇的峭度值小于3，因此确定VMD分解中分量的个数为6，即k＝6。为了去除重构信号中的低相关成分，选择PF₁-PF₆进行重构，得到重构信号记为x′₁(t)。

4、对重构信号x′₁(t)进行VMD分解，其中k＝6，α＝3000，分解结果如图5所示。从图5中可以看出，经过ELMD降噪后，信号经VMD分解得到6个IMF分量和一个残余分量，为了进一步选择IMF分量，计算0.3Mpa、0.25Mpa和0.2MPA下各IMF分量的多尺度熵值，结果见表2。

表2不同压力下各IMF分量的多尺度熵值

从表2中可以看出IMF₁、IMF₂、IMF₃、IMF₄的多尺度熵值较大，IMF₅和IMF₆分量的多尺度熵值较小。因此，信号经VMD分解后得到的IMF分量中，IMF₁、IMF₂、IMF₃、IMF₄包含泄漏信号大量的特征信息，其他分量几乎不含有效泄漏信息。故本文选取IMF₁、IMF₂、IMF₃、IMF₄作为主要分量进行重构。这样可以大幅度的减少背景噪声、残余分量对泄漏信号的影响。

5、经过上述过程选择IMF₁—IMF₄分量，提取每个IMF分量的时域、频谱进行重构，得到观测信号Nimf(t)＝[b₁(t)，b₂(t)…b_n(t)]；构成的多维观测信号很好解决了盲源欠定问题，计算相关矩阵R_x的特征值，然后对观测信号进行源数目估计。IMF₁-IMF₄分量对应的特征值分别为1.4393、1.0558、1.0265、0.497。前3个特征值较大，因此估计泄漏信号数目为3。利用独立分量分析的盲源分离方法得到管道上游的3个估计信号Y_μ(t)如图6所示，同理可以得到下游的3个估计信号Y′_μ(t)。在实际多源泄漏定位中，由于经盲源分离后分离出的信号均为泄漏信号，无法从分离后的时域图区别哪个估计信号对应哪个泄漏源，因此计算上游三个分离信号与下游三个分离信号两两间的相关系数，结果见表3。具有相关性的信号相关系数较大，以此来区别信号是否来自同一泄漏源，获取对应的上下游泄漏信号。

表3上下游估计泄漏信号相关系数表

由表3可以看出，Y₁和Y′₁、Y₂和Y′₃、Y₃和Y′₂分别为同一泄漏源。因此对Y₁和Y′₁进行互相关分析得到其互相关函数的时频谱峰值对应的频率为27.1945kHz，峰值对应的延时为0.01749s，根据经验值，声发射信号在水载作用下的波速为1500m/s，以此将延时与声速带入时差定位公式可以确定上游传感器到泄漏位置的距离7.88m，而实际距离为8m，因此定位绝对误差为0.12m，相对误差为1.5％，以此类推可计算出剩余两个泄漏点的位置以及相对误差，为了验证其有效性，与直接互相关和基于VMD分解的管道泄漏定位进行对比，结果如表4所示。

表4 3种方法定位结果比较

根据表4中3种方法定位结果和定位误差可以得到，直接互相关结果的平均误差为7.45％，VMD方法定位平均误差为4.3％，本发明提出的方法最小定位误差仅为1.28％。从表中可以很明显地看出直接互相关定位结果的误差远远大于另外两种，因为VMD方法能够自适应的分解信号，去除一定的干扰，本发明所提出的方法使原始泄漏信号信号经过改进VMD分解后提取了更多的泄漏信息，消除了VMD分解层数K值对分解结果的影响，减少了VMD分解后重构的误差，并且能够实现管道多点泄漏的精确定位，以上数据结果证明了所述方法的有效性。以上述依据本发明的理想实施例为启示，通过上述的说明内容，可以在不偏离本发明的范围内，进行多样的变更以及修改。本项发明的技术范围并不局限于说明书上的内容，必须要根据权利要求范围来确定其技术性范围。

Claims (6)

1.一种基于改进VMD的管道多点泄漏定位方法，其特征在于：包括如下步骤

S1采集管道原始泄漏信号；

S2对所述原始泄漏信号进行总体局域均值分解，得到若干PF分量；

S3计算各PF分量的相关系数，根据所述相关系数筛选出所需PF分量，并根据筛选出来的PF分量进行信号重构，确定变分模态分解的k值；

S4对重构信号进行变分模态分解，得到若干IMF分量，计算各IMF分量的多尺度熵值，并根据各IMF分量的多尺度熵值筛选IMF分量；

S5对筛选出来的IMF分量进行信号重构，得到观测信号，并利用独立分量分析的盲源分离方法对所述观测信号进行分离处理，得到估计泄漏信号，根据所述估计泄漏信号通过互相关定位算法进行管道泄漏定位；

步骤S4具体为

S4.1：对重构信号x′(t)进行VMD分解得到k个IMF分量；

S4.2：通过计算多尺度熵进行IMF分量的筛选；

所述多尺度熵的获取方法为：

(1)设IMF分量u_k的原始时间序列为{x₁，x₂，...，x_N}，长度为N，设定其嵌入维数w，相似容限z，对其进行粗粒化变换，得到新的时间序列u_k(τ)：

其中，k＝1，2，...，s/τ，μ，k为正整数且μ≠k，s为离散时间序列长度，τ为尺度因子，原始序列被分割成τ段且每段长为s/τ的粗粒序列，当τ＝1时，新的时间序列就是原始序列；

(2)对每个长度为N/τ的粗粒时间序列分别求样本熵，得到其样本熵值，并绘制成尺度因子的函数，按公式(5)得到多尺度熵值MSE(u_k，τ，w，z)：

式中SampE{u_k ^τ，w，z，s}为样本熵，w是嵌入维数；

(3)根据所述多尺度熵值的大小来选取最优的IMF分量，最终得到各IMF分量的多尺度熵值记为：MSE₁，MSE₂，MSE₃…MSE_n。

2.根据权利要求1所述的基于改进VMD的管道多点泄漏定位方法，其特征在于：步骤S2中总体局域均值分解方法为

S2.1：在原始泄漏信号x(t)中加入设定等级的白噪声n_m(t)得到混合泄露信号x_m(t)，所述泄漏信号表示为：

x_m(t)＝x(t)+n_m(t) (1)

其中，m为分解次数，t是时间，n_m(t)为白噪声信号，x_m(t)为加入白噪声后的混合泄漏信号；

S2.2：对混合泄漏信号x_m(t)进行LMD分解，再做总体平均得到最终PF分量，记为：

其中，ε_n，m(n＝1，2，...，.N)，N为正整数；ε_n，m为第m次分解得到的第n个初始PF分量，将N个初始PF分量M次分解的均值y_n(n＝1，2，...，N)作为最终PF分量。

3.根据权利要求1所述的基于改进VMD的管道多点泄漏定位方法，其特征在于：步骤S3具体为

S3.1采用皮尔逊积矩相关系数计算各PF分量y_n同原始信号x(t)的相关系数r；

S3.2根据相关系数r选取相关系数绝对值大于a的PF分量，然后对其计算峭度值，舍弃峭度值小于b的PF分量得到筛选后的PF分量；

S3.3最后确定VMD分解中的分量个数k，然后对筛选后的PF分量进行信号重构，得到重构信号x′(t)。

4.根据权利要求3所述的基于改进VMD的管道多点泄漏定位方法，其特征在于：a取值0.5，b取值3。

5.根据权利要求1所述的基于改进VMD的管道多点泄漏定位方法，其特征在于：通过互相关定位算法计算上下游盲源分离后的估计泄漏信号Y_μ(t)和Y′_μ(t)的时间延迟Δt，依据式(11)确定管道泄漏点的位置，

式中，x为估计泄漏定位值，l为上下游的传感器距离，v是泄漏声发射信号在管道中的传播速度；Δt为时间延迟。

6.根据权利要求1所述的基于改进VMD的管道多点泄漏定位方法，其特征在于：步骤S5具体为

S5.1：提取筛选后每个IMF分量的时域、频谱进行重构，得到观测信号

Nimf(t)＝[b₁(t)，b₂(t)…b_n(t)]；

S5.2：通过对观测信号Nimf(t)进行分离，计算观测信号的相关矩阵R_x＝E[Nimf(t)Nimf^H(t)]，根据特征值估计泄漏信号数目μ；并对观测信号Nimf(t)中心化，并作白化处理，最终得到盲源分离后管道的μ个估计泄漏信号Y_μ(t)；再对估计信号和原始泄漏信号进行相关系数计算，确定估计信号所对应的各泄漏源；

S5.3：通过互相关计算上下游盲源分离后的估计泄漏信号Y_μ(t)和Y′_μ(t)的时间延迟Δt，依据时差定位法确定管道泄漏点的位置：

式中，x为估计泄漏定位值，l为上下游的传感器距离，v是泄漏声发射信号在管道中的传播速度；Δt是时间延迟。

Images