基于时频分析和改进差分进化的多点泄漏定位方法

郎宪明， 朱永强， 李金娜， 曹江涛， 李 平， 宋华东

(1.辽宁石油化工大学 信息与控制工程学院，辽宁 抚顺 113001； 2.东北大学 信息科学学院，沈阳 110819； 3.沈阳仪表科学研究院有限公司，沈阳 110043)

管道是输送石油、天然气等流体的重要工具[1]。目前，有多种管道泄漏检测方法已被应用于管道泄漏监测。然而，现有的泄漏定位方法一般都是针对某一个泄漏点进行的，当流体管道局部管段出现压力过高、腐蚀、磨损、振动等，有时会出现2个或2个以上泄漏点。因此，针对单点泄漏的检测方法对管道进行多点泄漏定位往往失效[2]。

目前，基于声波的管泄漏检测方法主要是基于单点泄漏声波信号的智能分析[3-4]。利用Hilbert-Huang变换对管道首末站的声发射信号进行去噪处理，精确获得声波信号突变点时间，提高了泄漏定位的精度[5]。针对压力信号下降缓慢的缺点，提出了一种信号压缩方法，将压力信号转换为超声波声速信号，突出声波拐点的信息[6]。采用低通滤波与小波滤波相结合的方法，对管道首末站的压力信号进行实时高精度滤波，保证泄漏信息不被过滤掉[7]。由于小波去噪需要根据经验选择合适的小波基函数，而经验模态分解(empirical mode decomposition，EMD)不用凭经验选择基函数的优点，因此采用EMD对管道泄漏信号进行处理，提取管道泄漏信号的突变点信息，以减小外界环境噪声对泄漏点定位的影响[8]。变分模态分解(variational mode decomposition，VMD)是一种新的信号分解方法[9]，它可以降低EMD模态混叠对信号消噪的影响。

由于多数研究是关于管道的单点泄漏情况[10-15]，当管道同时发生多点泄漏，产生的瞬态泄漏声波信号混叠在一起，由于泄漏声波信号相互影响，以及外部噪声(如施工、管道阀门噪声等)叠加在管道泄漏声波信号上，从而影响瞬态声波信号在信号传播中的衰减规律，因此在混杂信号的情况下，多点泄漏信号不能有效地去噪。文献[16]采用扩展卡尔曼滤波方法对多点连续泄漏进行检测。文献[17]利用超声波声速的变化来进行多点连续泄漏进行定位。

当多点泄漏声波向管道首末站传播时，很难区分各个泄漏声波信号到达管道首末站的突变点时间。因此，基于单点泄漏的泄漏定位方法不适用于多点泄漏。然而，时频分析(time-frequency analysis，TFA)可以从时域和频域分析来分析非平稳的多点泄漏声波信号[18]，通过时频分析来确定多点泄漏的数量。

差分进化(differential evolution，DE)算法是一种新的群体智能优化算法[19]。它具有简单、易于实现、无梯度信息、参数少等优点，被广泛应用于工业过程参数优化中。然而，进化算法在进化后期容易陷入局部最优，收敛速度较慢。针对此问题，文献[20]引入了混合自适应局部搜索算法来提高DE算法的收敛速度。为了加快离散进化算法的收敛速度，文献[21]采用正交交替的方法来提高收敛速度。为了提高差分算法的收敛速度和精度，本文采用粒子群优化算法对差分进化算法进行改进，根据建立的多点泄漏定位目标函数，采用改进差分进化算法估计多点泄漏位置。

1 多点泄漏定位原理

1.1 改进变分模态分解方法

当VMD用于信号分解时，根据经验选择适当的模态分解个数m[22]。由于m决定了频率分辨率，如果选择m过小，则会出现信号欠分解，对泄漏声波信号变化不敏感，如果选择m太大，会产生伪分量，增加计算量。

将采集到的声波信号y(t)采用VMD分解，计算如下

(1)

式中，γ(t)是信号的重构误差。

能量函数的计算如下

Eγ=γ(t)2

(2)

Eγ能量越小，两个信号的相似度越高，如果误差能量为0，则说明两个信号相同。

由于VMD分解模态的数量对信号重构有影响，因此采用能量函数对VMD进行改进。改进变分模态分解(improved variational mode decomposing, IVMD)算法的计算步骤如下：

(1) VMD算法中的分解模态数量设定为m。

(2) 经VMD信号分解后，获得了一定数量的本征模态函数(intrinsic mode functions, IMFs)。

(3) 计算不同数量分解本征模态函数的误差能量函数。

(4) 绘制误差能量图，首先计算出各误差能量函数的平均值阈值。一般情况下，如果误差能量值远小于阈值，称为有效本征模态函数。有效本征模态函数包含泄漏声波信号的突变点时刻、信号幅值变化等信息。

因此，泄漏声波信号经IVMD选择有效本征模态函数对信号进行消噪。

1.2 多点泄漏定位算法

当管道出现多点泄漏同时发生时，多点泄漏引起的声波幅值变化混合在一起[23-24]，由于多点泄漏信号之间的相互作用，在管道首末站获得多点泄漏声波的叠加信号，多点泄漏声波传播和泄漏衰减的振幅描述为

(3)

(4)

式中：xi是i个泄漏点产生的泄漏声波信号原始振幅;a是衰减因子;xa和xb分别是管道首末站获取的多点泄漏声波信号;li是第i个泄漏点;L是管道的长度。

当出现单点泄漏时，根据式(3)和式(4)，按照首末站获得的泄漏声波信号幅值变化可计算出单点泄漏位置。

根据多点泄漏声波信号幅值变化的特点，假设每个泄漏点产生的声波信号原始振幅大小一样，根据式(3)和式(4)，则多点泄漏的定位计算如下

(5)

根据管道首末站的获得的多点泄漏声波信号振幅变化，根据式(5)可以计算出多点泄漏位置。然而，同时发生的多点泄漏引起的振幅变化是重叠的，无法区分多点泄漏引起的振幅变化。因此，根据已知多点泄漏数量的情况下，根据式(5)建立多点泄漏定位函数，然后通过最小化管道首末站的多点泄漏声波信号振幅误差来估计多点泄漏位置。

2 多点泄漏定位方法

2.1 时频分析

采用Wigner分布对管道多点泄漏声波信号振幅的变化进行分析[25-27]。

Sx(t,ω)=ASPW(x)

(6)

式中：x是多点泄漏信号的振幅变化;Sx(t,ω)为多点泄漏声波振幅变化的时频分布。

Wigner的分布描述如下

ASPW(t,ω)=

(7)

式中，g(t)和h(t)分别是多点泄漏声波信号时域和频域的光滑函数。

多点泄漏信号的振幅变化时频分布具有不同的分布曲线，可以根据分布曲线确定泄漏点的数量。

2.2 差分进化算法

DE算法是一种随机并行的搜索算法[28]。它是从随机的初始个体开始，按照一定的迭代规则。根据每个个体的适应度值，保留较好的个体，引导搜索过程达到最优解。

对于具有n个变量的全局优化问题。全局优化问题可转化为最小值问题，求解如下函数

(8)

式中：D是问题空间解的维数；aj和bj分别是χj的上限和下限。

DE算法主要包括变异操作、交叉操作和选择操作。

(1) 突变操作

(9)

式中：t是进化次数；χr1,χr2和χr3是随机选取的三个单独个体；ui是修正因子；F是突变因子。

(2) 交叉操作

为了增加种群多样性，采用了交叉。

(10)

式中：j=1,2,…,D，D是维度；RC是交叉率，RC∈[0,1]。

(3) 选择操作

由式(10)从试验向量和原始向量中选出较好的个体。

(11)

式中，f(χ)是χ的目标函数。

2.3 改进差分进化(improved differential evolution，IDE)算法

在优化过程中，DE存在收敛速度慢、早熟等缺点，严重影响了算法的性能。因此，采用了一种改进的差分进化算法。在DE算法中引入PSO原理，用PSO搜索整个种群，使其快速跳出局部最优，避免早熟现象[29]。

粒子群优化的个体变异操作如下

(12)

其中，χbest是优势个体。

IDE的流程图如图1所示。

图1 IDE流程图Fig.1 Flow chart of IDE

2.4 多点泄漏定位方法

多点泄漏定位方法的流程图如图2所示。从图2中，可以使用IVMD对泄漏声波信号进行去噪。传统的泄漏定位方法很难区分多点泄漏个数，多点泄漏声波信号经IVMD去噪后，然后采用TFA算法来估计泄漏点数量。然后，根据式(5)建立的多点泄漏定位目标函数，采用IDE估计多点泄漏的位置。

图2 多点泄漏定位的流程图Fig.2 Flow chart of multiple leaks localization method

3 现场试验

在中国石油某个成品油管线进行了多次多点泄漏试验。管道长度为20 000 m，管道内径为100 mm。管道上游压力为2 MPa，管道下游压力为0.5 MPa，流量为200 L/min，密度是830 kg/m3，环境温度是25 ℃，通过声传感器(PCB 106B)获取管道首末站多点泄漏的声波信号。声波数据的采集试验装置原理图如图3所示。管道试验系统由成品油管道和多点泄漏定位系统组成。多点泄漏定位系统主要由PCB106b、工控机和数据采集器NI CDAQ-9181组成。泄漏定位方法的关键是采集管道首末站多点泄漏声波信号的时间同步，本试验时间同步采用GPS卫星时钟信号校时，管道首末站的泄漏声波信号采样频率是1 000 Hz。为了模拟多点泄漏，3个泄漏位置分别距离管道首站14 200 m、15 500 m和15 800 m处，根据式(5)，模拟3点泄漏的泄漏阀门孔径分别为15 mm、15 mm和15 mm。实际管道运行过程中，超过3点同时泄漏会给管道运行造成非常严重的事故。另外，3点同时泄漏是目前研究热点问题，本文只针对管道3点同时泄漏情况进行相关试验。

图3 试验装置示意图Fig.3 Schematic diagram of test apparatus

在距离首站14 200 m位置，模拟泄漏孔径为15 mm，泄漏持续时间为2 s，管道首末站分别获得的单点泄漏的声波信号变化如图4(a)所示。在距离首站14 200 m和15 500 m的2个泄漏位置，模拟泄漏孔径分别是15 mm和15 mm，泄漏持续时间为2 s，管道首末站分别获得的2点泄漏的声波信号变化如图4(a)所示。另外，在距离首站14 200 m、15 500 m和15 800 m的3个泄漏位置，模拟泄漏孔径分别是15 mm、15 mm和15 mm，泄漏持续时间为2 s，管道首末站分别获得的3点泄漏的声波信号变化如图4(c)所示。

(a) 单点泄漏

(b) 2点泄漏

(c) 3点泄漏图4 获取的多点泄漏声波信号Fig.4 Acquired acoustic wave signal under multiple leaks

从图4可以看出，随着泄漏点数量的增加，声压信号幅值逐渐增大。管道末站获得单点泄漏的声压信号幅值波动范围为0～2.2，2点泄漏时的声压信号幅值波动范围为0～6.2，3点泄漏的声压信号幅值波动范围为0～7.8。

根据PCB106b传感器获得的泄漏声波信号，经VMD进行5层本征模态分解，5个IMFs可以较好地消除被测信号的噪声。根据声波信号的5个IMFs，计算误差能量函数并绘制误差能量图，本征模态函数的误差能量图如图5所示。

从图5可以看出，当本征模态的数目为2个时，误差能量函数值为0.92，这些本征模式包含更有效的泄漏声波信号突变点、声波信号幅值变化等信息。因此，我们可以设定阈值为0.93，并且在本试验中选择了2个本征模态函数。泄漏声波信号经VMD的5层本征模态分解和经IVMD的2个本征模态函数的声波信号消噪对比曲线如图6所示。

图5 本征模态函数的误差能量函数值Fig.5 Error energy diagram of intrinsic mode

图6 经VMD和IVMD消噪后的声波信号对比曲线Fig.6 Comparison curve of acoustic signals de-noised with VMD and IVMD

从图6中可以看出，IVMD的2个IMFs可以较好地恢复泄漏声波信号的变化，并且相比VMD的计算量要小。

根据误差能量函数计算，IVMD选择2个IMFs对多点泄漏声波信号进行信号消噪处理。通过IVMD对单点泄漏、2点泄漏、3点泄漏三种情况下的泄漏声波信号进行处理。根据管道首末站传感器PCB106b获得的泄漏声波信号进行IVMD信号消噪处理，经IVMD和EMD对单点泄漏、2点泄漏、3点泄漏声波信号的去噪对比曲线如图7所示。

(a) 单点泄漏

(b) 2点泄漏

(c) 3点泄漏图7 经IVMD和EMD消噪后的多点泄漏声波信号Fig.7 Multiple leaks signals de-noised with IVMD and EMD

通过对单点泄漏、2点泄漏、3点泄漏信号消噪后进行TFA分析，其时频分析曲线如图8所示。

从图8可以看出，单点泄漏、2点泄漏、3点泄漏泄漏引起的频率分量变化明显，可以通过TFA有效地区分。因此，可以通过对多点泄漏声波信号的TFA分析，估计出泄漏点的数量。

通过TFA方法估计泄漏点的数量，根据式(5)建立多点泄漏定位目标函数，采用IDE对多点泄漏定位目标函数进行求解，估计出多点泄漏的位置。在DE算法中，突变因子F=0.5、迭代次数t=100和维数D，其中维数D由TFA估计出泄漏点的数量决定，交叉率RC=0.9和种群数量为100。另外，将IDE的参数与DE的数值一样，粒子群数量是40，维数是20，学习因子是2，惯性权重是1.2。

当管道同时发生3点泄漏时，IDE和DE的适应度函数收敛性能曲线如图9所示。

图9 IDE和DE的适应度函数收敛曲线Fig.9 Fitness function convergence performance of IDE and DE

从图9可以看出，当IDE算法的迭代次数为75时，适应度函数最终收敛到0.35左右，而DE算法的适应度函数为72，但适应度函数最终收敛到6，因此，IDE算法在估计泄漏点位置求解速度和求解精度方面相比DE算法更具优势。由于IDE算法通过引入PSO算法对其进行了改进，使得IDE算法在求解多点泄漏定位函数时跳出局部极值，收敛到全局最优解。因此，IDE算法求解多点泄漏定位函数。

管道首末站采用PCB106b传感器获得多点泄漏声波信号数据，信号采样率为1 000 Hz。建立管道单点泄漏、2点泄漏、3点泄漏共3种工况数据，每种工况进行30组试验。采用IVMD方法对单点泄漏、2点泄漏、3点泄漏的泄漏声波信号进行降噪处理。利用2个本征模态函数分别对单点泄漏、2点泄漏、3点泄漏的泄漏声波信号进行重构。同时，EMD算法选取6～12个IMF分别对单点泄漏、2点泄漏、3点泄漏的泄漏声波信号进行消噪。然后，利用TFA和IDE对采集到的多次泄漏声数据进行处理。不同方法的泄漏试验定位结果如表1所示。由于多点泄漏是一个小概率事件，并且超过3点同时泄漏会引起流体管道出现重大事故发生。目前，3点泄漏是研究的热点问题，在现场试验中，不考虑管道出现超过3点同时泄漏的情况。

从表1可以看出，多点泄漏的互相关分析中只有一个峰值。因此，无法计算2点或3点泄漏的泄漏位置。TFA可以用来确定多个泄漏点的数量，然后使用IDE来估计多点泄漏的位置。在距离管道首站14 200 m，15 500 m和15 800 m位置同时发生3点泄漏，经过IVMD对泄漏声波信号进行消噪，然后将消噪后的信号进行TFA计算，估计出泄漏点数量，并根据式(5)，建立多点泄漏定位目标函数，采用IDE算法求解此目标函数。经TFA和IDE算法估计的3点泄漏位置如图10所示。

表1 不同模型的多点泄漏定位结果对比

图10 经TFA和IDE估计的3点泄漏位置Fig.10 Three leaks localization by TFA and IDE

从图10可以看出，用IDE算法迭代次数为14，估计的第1点泄漏定位值为14 230 m，定位误差为30 m；第2点泄漏的定位估计值为15 530 m，IDE算法迭代次数为10次，定位误差为30 m；第3点漏定位估计值为15 818 m，IDE算法迭代次数为16次，定位误差为18 m。因此，该方法可以估计多点泄漏的位置。由于管道实际运行中，管道首末站的次声波传感器位置及模拟多点泄漏位置不能随意调整。因此，在接下来的研究中会根据实际工况调整传感器及泄漏点的位置，体现多点泄漏定位方法的实用性。

4 结 论

本文提出一种基于时频分析和改进的差分进化算法的多点泄漏定位方法。通过IVMD对泄漏声波信号进行消噪处理，将消噪信号进行TFA分析，估计多点泄漏的泄漏点数量，根据估计的泄漏点数量及泄漏定位目标函数，采用IDE算法估计多点泄漏位置。在某成品油管线进行了多点泄漏定位试验，泄漏点定位误差最小是18 m。试验结果表明，该方法可以在3点泄漏同时发生的情况下，有效地估计3点泄漏的位置。在后续的研究中，课题组将对管道3点以上同时泄漏点的定位进行研究。