一种城市气体管道微小泄漏诊断方法

1.本发明涉及城市埋地管道泄漏检测技术领域，具体涉及一种城市气体管道微小泄漏诊断方法。


背景技术：

2.随着城市化进程的加快，管道作为城市基础设施也迅速发展，对输送城市必需能源、提高人民生活水平质量做出了巨大贡献。然而随着管道使用年限的增长、施工缺陷和腐蚀以及人为破坏，城市管道泄漏常常发生，进而引发火灾、爆炸等恶性事故，造成人员及财产损失。再者，与长输油气管道相比，城市管道的管径更细、输送压力更低，其早期泄漏信号更微弱，泄漏检测与定位更加困难。因此，有关城市管道微小泄漏诊断的研究是一个非常有实用价值和社会意义的课题。
3.近年来，由于次声波检测技术的高效性与便捷性，在管道泄漏检测领域被大量应用，积累了一定的理论知识和技术成果，基于次声波的管道泄漏检测与定位技术也日趋成熟。当管道发生泄漏时，管道内外产生压力差，管内介质由泄漏口喷射出。管内流体与管壁相互作用而诱发不同频率的振动，导致空气振动产生不同频率的声波信号，并向管道两端传播，在传播过程中高频信号迅速衰减，当声波信号传递到管道两侧传感器时，主要成分为次声波。次声波检测技术是借助次声波检测系统对管道中的次声波信号进行采集、记录、分析并对声源的强度、位置、发生条件等性质进行评价的技术。其检测原理是：泄漏孔处声波信号会沿着管道上下游传播。在泄漏点的上、下游安装次声波传感器，用于接收来自泄漏点的次声波信号，由于上、下游传感器与泄漏点距离不同，从而在检测同一信号突变时会产生时间差，根据此时间差就能够确定泄漏点的具体位置。但，目前还存在城市管道微小泄漏(微小泄漏为泄漏孔径在1mm～5mm之间)不易察觉、定位误差较大的问题。


技术实现要素：

4.针对现有城市管道微小泄漏(微小泄漏为泄漏孔径在1mm～5mm之间的泄漏)不易察觉、定位误差较大的问题不足问题，本发明的目的在于提出一种城市气体管道微小泄漏诊断方法，旨在解决城市管道(包括金属管道和非金属管道)微小泄漏的诊断问题。
5.一种城市气体管道微小泄漏诊断方法，包括以下步骤为：
6.s1：结合城市管道的实际运行情况，搭建埋地气体管道模拟试验系统(本申请同时适用于金属管道和非金属管道的泄漏诊断，且处理过程一致，因此本申请中的金属管道和非金属管道统称为气体管道)，分别采集埋地气体管道模拟试验系统在不同压力下管道未泄漏和微小泄漏的管内次声波的时域信号，并对所述时域信号进行频谱转换得到模拟试验系统的次声波频谱，总结出城市气体管道微小泄漏的频谱特征规律，进而得到城市气体管道微小泄漏频谱判断准则。
7.频谱转换具体过程为：
8.a、将采集到的次声波时域信号x(t)转换为频域谱密度函数x(f)：
[0009][0010]
式中：x(f)为频域谱密度函数；ω为模拟角频率，ω＝2πf，f表示频率；e为自然对数底；j为虚部单位；t表示时刻。
[0011]
b、计算频域中谱密度函数x(f)的幅值谱估计函数
[0012]
将谱密度函数x(f)分为m个等长数据段，将谱密度函数x(f)的m点观察数据x
m
(f)视为有限信号；然后对x
m
(f)进行傅里叶变换，将其变换为x
m
(k)；最后取x
m
(k)幅值的平方并除以数据段数m，作为谱密度函数x(f)的幅值谱估计即：
[0013][0014]
c、根据幅值谱估计函数绘制次声波信号频谱：
[0015]
利用matlab软件的绘图功能，绘制幅值谱估计函数的函数图像，即得到与之对应的次声波信号频谱图。
[0016]
进一步地，城市气体管道微小泄漏频谱判断准则为城市气体管道未泄漏次声波频谱信号满足只有一个幅值高峰的特性，且相应频率不大于5hz；而城市气体管道微小泄漏时，次声波频谱具有至少2个幅值高峰波动特性，且相应频率范围在5hz～20hz内。
[0017]
s2：对待测城市管道进行次声波时域信号采集，将采集的待测城市管道的次声波时域信号进行频谱转换，得到待测城市管道的次声波频谱，对照所述城市气体管道微小泄漏频谱判断准则判定测试管道是否发生泄漏。
[0018]
s3：若管道泄漏，综合考虑城市燃气管道实际情况，建立基于次声波检测技术的管道泄漏定位模型。
[0019]
管道泄漏定位模型为：假设s(t)，n1(t)，n2(t)互不相关，
[0020][0021]
式中：l1为泄漏位置到上游次声波传感器的距离，单位m；l为上游次声波传感器和下游次声波传感器之间的距离，单位m；τ0为时间延迟，单位s；v为次声波信号的波速，单位m/s；s(t)为未知的原始泄漏源信号；n1(t)和n2(t)为表示环境噪声的零均值随机噪声信号，t为时刻。
[0022]
s4：考虑噪声信号对次声波的干扰，采用相干函数提取频谱中连续峰值对应的频带为特征频带，若存在至少2个连续峰值，则取最大值所在频带为特征频带，再使用带通滤波器提取特征频带中的泄漏信号。
[0023]
进一步地，提取泄漏信号的具体方法如下：
[0024]
设管道上、下游次声波传声器接收到的次声波时域信号分别为x1(t)和x2(t)，则：
[0025][0026]
式中：x1(t)和x2(t)分别是上、下游次声波传感器输出的时域信号；s(t)为未知的原始泄漏源信号；n1(t)和n2(t)为表示环境噪声的零均值随机噪声信号；a是声波衰减系数；τ0为次声波传递至管道上、下游两端的时间延迟，单位s；t为时刻；
[0027]
x1(t)和x2(t)信号的互相关函数由下式定义：
[0028][0029]
式中：为输出的时域信号之间的互相关函数；e[
·
]为期望算子；τ为时间延迟；
[0030]
对公式(5)进行傅里叶变换，得到信号x1(t)和x2(t)的互谱公式：
[0031][0032]
式中：为两输出的时域信号的互谱函数；s
ss
(f)为两个传感器接收到的泄漏源信号的互谱函数，j为虚部单位，f为频率，为两输出信号的相位差：
[0033][0034]
由式(7)可知：在信号源的频带中，若两输出的时域信号的时间延迟相同，则随频率变化而线性变化，斜率为2πτ0；
[0035]
将输出的时域信号x1(t)和x2(t)分为m个等长数据段，给每个数据段加一个余弦锥形窗口，以抑制频谱估计中的旁瓣泄漏；再对每段数据进行傅里叶变换，设i为段数，则x
1i
(t)和x
2i
(t)经过傅里叶变换为x
1i
(f)和则第i段数据的互谱为计算两个信号的互谱函数为：
[0036]
利用互谱函数表示单边互谱
[0037][0038]
式中：为互谱密度函数，为正交谱密度函数，两者分别为的实部和虚部；对应的相位差表示为：
[0039][0040]
通过welch平均周期图法得到两个输出的时域信号的单边自谱g
i
和g
π
，两个输出的时域信号x1(t)和x2(t)的相干函数表示为：
[0041][0042]
泄漏输出信号的自谱、源信号自谱和噪声自谱之间的关系是：
[0043][0044][0045]
式中：g
ss
(f)为源信号自谱，为噪声自谱；
[0046]
由式(6)和(10)可将相干函数转换为：
[0047][0048]
式中：和分别表示信号x1(t)和x2(t)的信噪比的倒数，假设两组
输出的时域信号在频率f下的源信号的信噪比分别为：snr1和snr2，将snr1和snr2代入式(13)，相干函数可表示为：
[0049][0050]
若假设x1(t)中没有噪声信号，则snr1无穷大。
[0051][0052]
转换得：根据相干函数，因此，频率为f的源信号的最小信噪比为snr2；同理，当x2(t)中没有噪声信号时，频率为f的源信号的最小信噪比为snr1；因此证得：频率为f的源信号的最小信噪比为：
[0053][0054]
由公式(16)能够知道：相干函数在频域测量两个信号的相关程度，一个频带内的相干函数越大，该频带内源信号的信噪比越高，因此，根据其幅度来提取源信号的频带，从而得到输出的时域信号的特征频带；
[0055]
使用带通滤波器提取特征频带中的泄漏信号，以特征频带的上下限为带通滤波器的上、下截止频率f0和f1，求出带通滤波器的频率带宽δf＝f1‑
f0，中心频率带通滤波器g(f)为：
[0056][0057]
s5：采用广义互时频法计算提取的上、下游次声波输出泄漏信号的时间延迟。
[0058]
进一步地，采用广义互时频法精确计算提取的上、下游次声波输出泄漏信号的时间延迟的具体方法为：
[0059]
通过计算次声波信号的互功率谱，再结合带通滤波器和互谱函数求解修正加权函数ψ
12
(f)，并在频域范围内加权，加强信号中的谱分量，提升信号的信噪比，得到广义互时频函数；再利用广义互时频函数将频域反变换到时域，根据信号相关性原理计算次声波信号的时差，得到次声波信号的时间延迟。
[0060]
进一步地，计算提取的上、下游次声波输出泄漏信号的时间延迟具体步骤为：
[0061]
管道上、下游次声波传声器接收到的次声波时域信号x1(t)、x2(t)的互功率谱为：
[0062][0063]
式中：x(f)，s(f)，n(f)分别为时域信号x(t)，s(t)，n(t)的频域表示，j为虚数单位，τ1为泄漏信号从泄漏位置传递至上游次声波传感器的时间，τ2为泄漏信号从泄漏位置传递至下游次声波传感器的时间；
[0064]
由于x(t)，s(t)，n(t)之间彼此不相关，则式(18)简化为：
[0065][0066]
结合带通滤波器g(f)和互谱函数求解修正加权函数ψ
12
(f)：
[0067][0068]
对互功率谱函数进行加权，x1(t)、x2(t)的广义互时频函数表示为：
[0069][0070]
根据式(21)绘制广义互时频分析图，根据图中相关函数最大值对应的采样点数和信号采样频率求得次声波信号传递至管道上、下游两端的时间延迟τ0；
[0071]
τ0＝nf
s
ꢀꢀꢀ
(22)
[0072]
式中：n为最大幅值点对应的采样时间，μs；f
s
为信号采样频率，hz。
[0073]
s6：采集管道泄漏状态ts内各位置处压力的平均值作为管道压力值p，并根据管道压力值p求得次声波信号的波速。
[0074]
s6.1：在城市气体管道上设置若干压力传感器，再通过压力监测平台调取采集管道泄漏状态ts内各测压位置压力的平均值作为管道压力值，设时间步长为1s，位置步长为2.5m，即：
[0075][0076]
式中：t为管道泄漏时间，s；t为时刻，s；l为上游次声波传感器和下游次声波传感器之间的距离，m；l0为管道压力传感器位置，m；[
·
]为取整符号。
[0077]
s6.2：在城市管道实际运行过程中，管内次声波传播速度并非理想状态下的声速，尽管两者在数值上相差不大，但由此进行管道泄漏定位而产生的误差却不可忽视。因此，考虑到管内介质压缩系数、介质密度、管道大小等因素，管道泄漏次声波信号波速为：
[0078][0079]
式中：δ为介质压缩系数，无量纲；ρ为介质密度，kg/m3；e为管道弹性模量，pa；d为管道直径，mm；e为管道壁厚，mm；c为修正系数，无量纲。
[0080]
其中，气体压缩系数δ的计算公式为：
[0081][0082]
式中：z为压缩因子，p为管内压力，pa。对于中低压干燥气体，压缩因子z＝1。
[0083]
s7：最后将精确的时间延迟和声波速度代入定位公式对泄漏点精确定位。
[0084]
本发明的有益效果为：本发明的诊断方法精确度高、简便易行，通过不同工况的次声波频谱分析试验总结了城市管道微小泄漏的频谱特征规律，建立了城市管道微小泄漏频谱判断准则，其泄漏判断效率高、误判率低；采用相干函数确定频谱中的特征频带，再使用带通滤波器提取特征频带中的泄漏信号，此步骤能有效地滤除环境噪声；采用广义互时频法进一步精确计算了次声波信号的时间延迟，并将时间延迟代入管道泄漏计算模型，实现了管道泄漏精确定位。
附图说明
[0085]
下面结合附图和具体实施方式来详细说明本发明：
[0086]
图1是本发明的流程图；
[0087]
图2是试验管道及传感器仪表布置图；
[0088]
图3是试验管道未泄漏情境下的次声波频谱图；
[0089]
图4是试验管道泄漏情境下的次声波频谱图；
[0090]
图5是实施例1的泄漏检测时的次声波频谱图；
[0091]
图6是实施例1的上、下游泄漏次声波信号的相干函数图；
[0092]
图7是实施例1的广义互时频分析图；
[0093]
图8是实施例1的管道压力分布图；
[0094]
图9是实施例1的上、下游泄漏信号时域波形图。
[0095]
图10是实施例2的泄漏检测时的次声波频谱图；
[0096]
图11是实施例2的上、下游泄漏次声波信号的相干函数图；
[0097]
图12是实施例2的广义互时频分析图；
[0098]
图13是实施例2的管道压力分布图；
[0099]
图14是实施例2的上、下游泄漏信号时域波形图。
[0100]
图15是实施例3的泄漏检测时的次声波频谱图；
[0101]
图16是实施例3的上、下游泄漏次声波信号的相干函数图；
[0102]
图17是实施例3的广义互时频分析图；
[0103]
图18是实施例3的管道压力分布图；
[0104]
图19是实施例3的上、下游泄漏信号时域波形图。
[0105]
图中附图标记的含义为：1
‑
入口球阀，2
‑
上游次声波传感器，3
‑
第一压力传感器，4
‑
第二压力传感器，5
‑
第三压力传感器，6
‑
第四压力传感器，7
‑
泄漏孔，8
‑
第五压力传感器，9
‑
第六压力传感器，10
‑
下游次声波传感器，11
‑
出口球阀。
具体实施方式
[0106]
为使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施方式，进一步阐述本发明。
[0107]
实施例1：
[0108]
如图1所示，本发明的管道微小泄漏诊断方法，包括以下步骤：
[0109]
s1：本发明结合城市管道的实际运行情况，建立了由气体管道系统、次声波采集系统和压力监测系统三部分构成的试验管道系统。在实验室内搭建如图2所示的气体管道系统模型，管道为u型pe管，管道型号为φ63
×
1.5mm，总长度为12.8m；在管道的左侧直线段，自管道进口端依次设置有入口球阀1、上游次声波传感器2、第一压力传感器3、第二压力传感器4；在管道的右侧直线段，自管道出口端依次设置有出口球阀11、下游次声波传感器10、第六压力传感器9、第五压力传感器8、泄漏孔7；在管道的弯线段设置依次设置有第三压力传感器5和第四压力传感器6；6个压力传感器在管道上分别为相邻两个压力传感器相邻2.5米布置。
[0110]
泄漏孔7距上游次声波传感器2的9.05m处，采用可调节式泄漏孔7模拟城市管道发
生微小泄漏，采用空压机输气模拟城市气体管道的运行过程。在管道的两端分别安装上、下游次声波传感器，并通过电缆连接到次声波采集仪，再由网线传至pc端的分析软件中。在管道上每隔两米布置压力传感器，通过物联网技术，将采集的压力数据通过网关上传至gprs模块(物联网卡)继而上传至pc端的压力监测平台。
[0111]
s1.1：分别在0.1mpa、0.2mpa和0.3mpa管道压力下，采集多组管道未泄漏和1～5mm微小泄漏的次声波时域信号x(t)。
[0112]
s1.2：对采集所得次声波时域信号进行频谱转换，得到与之对应的频谱，以不同压力下管道未泄漏与泄漏各一组数据为例，如图3和图4所示。通过大量的对比、分析，总结得到管道未泄漏和微小泄漏次声波信号的频谱特征规律，最终建立起城市管道微小泄漏频谱判断准则如下：城市管道未泄漏次声波频谱信号满足只有一个幅值高峰的特性，且相应频率不大于5hz；而微小泄漏时，次声波频谱具有至少2个幅值高峰波动特性，且相应频率范围在5hz～20hz内。
[0113]
s2：将试验管道的运行压力调节至0.2mpa，待管道平稳运行后打开泄漏孔7模拟管道2mm微小泄漏，采集管道次声波时域信号，并对其进行频谱转换，得到城市管道的次声波信号频谱如图5所示。由图5可知，该次声波信号频谱的幅值波动主要集中在5hz和8.5hz处，此外在13.5hz和24hz处也存在较弱的幅值波动，即该频谱具有多峰值波动特性且其所处频率范围属于5hz～20hz。对照s1.2中的频谱判断准则，判断管道此时处于微小泄漏。
[0114]
s3：综合考虑城市燃气管道实际情况，建立基于次声波检测技术的管道泄漏定位模型。
[0115][0116]
s4：采用相干函数提取频谱中的特征频带如图6所示。由图6可知函数在14.5hz处取最大值，取其所在频带12～16hz为特征频带。再使用带通滤波器提取特征频带中的泄漏信号，以特征频带的上下限为带通滤波器的上、下截止频率，f0＝12hz和f1＝16hz，中心频率f
c
＝14hz。
[0117]
s5：广义互时频法精确计算次声波信号的时间延迟，如图7所示。
[0118]
由图7可知，上、下游次声波泄漏信号的相关函数在0～200hz频率范围内取得最大值，且相关函数最大值对应的采样时间n＝100μs，信号采样频率f
s
＝187根据该采样点数精确计算出次声波信号传递至管道上、下游两端的时差：τ0＝187
×
10
‑4＝0.0187s。
[0119]
s6：采集管道泄漏状态ts内各测压位置压力的平均值作为管道压力值，并根据管道压力值求得次声波信号的波速。
[0120]
s6.1：通过压力监测平台调取管道上各压力传感器采集的压力数据，如图8所示。然后取管道采集状态30s内各测压位置处压力的平均值作为管道压力值，时间步长为1s，位置步长为2.5m，即：
[0121][0122]
s6.2：由式(24)可知，且已知壁厚e为1.5mm，空气密度ρ为1.29kg/m3，弹性模量e为2.9gpa，泄漏孔7直径d为2mm，修正系数c为0.91，δ为气体压缩系数与管内压力有关，联立式
(24)和(25)并代入参数值，得到次声波信号波速：
[0123][0124]
s7：最后将精确的时间延迟和次声波速度代入定位公式对泄漏点精确定位。
[0125][0126]
计算定位误差为：
[0127][0128]
与未经处理的上、下游泄漏信号粗定位结果进行对比。上、下游次声波泄漏信号的时域波形图如图9所示。
[0129]
由图9可知，上、下游泄漏信号时域波形中出现幅值尖峰的时间分别为t
′1＝0.0279，t
′2＝0.0492，则次声波信号时差τ
′
＝t
′2‑
t
′1＝0.0213s，因此根据次声波定位原理得到粗定位结果：l
′1＝(12.8+0.0213
×
340)/2＝10.02m。
[0130]
计算定位误差为：
[0131][0132]
计算结果表明，本发明方法的定位精度比粗定位精度提高了3.29％。
[0133]
实施例2：
[0134]
本实施例与实施例1的区别在于：待管道平稳运行后打开泄漏孔7以模拟管道1mm微小泄漏。进行验证本发明的诊断方法，具体如下：
[0135]
采集管道次声波时域信号，并进行频谱转换如图10所示。由图10可知，该次声波信号频谱的幅值波动主要集中在10.5hz和18.5hz处，即该频谱具有多峰值波动特性且其所处频率范围属于5hz～20hz。对照城市管道微小泄漏频谱判断准则，判断管道此时处于微小泄漏。
[0136]
采用相干函数提取频谱中的特征频带如图11所示。由图11可知函数在13.5hz处取最大值，取其所在频带12～15hz为特征频带。再使用带通滤波器提取特征频带中的泄漏信号，以特征频带的上下限为带通滤波器的上下截止频率，f0＝12hz和f1＝15hz，中心频率f
c
＝13.5hz。
[0137]
为进一步精确计算次声波信号时差，对管道上、下游泄漏信号x1(t)，x2(t)进行广义互时频分析，如图12所示。
[0138]
由图12可知，上、下游次声波泄漏信号的相关函数在0～200hz频率范围内取得最大值，且相关函数最大值对应的采样时间n＝527μs，信号采样频率f
s
＝34hz，根据该采样点数精确计算出次声波信号传递至管道上、下游两端的时差：τ0＝527
×
34
×
10
‑6＝0.0179s。
[0139]
通过压力监测平台调取管道上各压力传感器采集的压力数据，如图13所示。然后取管道采集状态30s内各测压位置处压力的平均值作为管道压力值，即：
[0140][0141]
将管道压力值代入次声波速度计算公式得：
[0142][0143]
将精确的时间延迟和次声波速度代入定位公式对泄漏点精确定位。
[0144][0145]
计算定位误差为：
[0146][0147]
与未经处理的上、下游泄漏信号粗定位结果进行对比。上、下游次声波泄漏信号的时域波形图如图14所示。
[0148]
由图14可知，上、下游泄漏信号时域波形中出现幅值尖峰的时间分别为t
′1＝0.0315s，t
′2＝0.0565s，则次声波信号时差τ
′
＝t
′2‑
t
′1＝0.025s，因此根据次声波定位原理得到粗定位结果：l
′1＝(12.8+0.025
×
340)/2＝10.65m。
[0149]
计算定位误差为：
[0150][0151]
计算结果表明，本发明方法的定位精度比粗定位精度提高了8％。
[0152]
实施例3：
[0153]
本实施例与实施例2的区别在于：待管道平稳运行后打开泄漏孔7以模拟管道5mm微小泄漏，进行验证本发明的诊断方法，具体如下：
[0154]
采集管道次声波时域信号，并进行频谱转换如图15所示。由图15可知，该次声波信号频谱的幅值波动主要集中在6hz和18hz处，即该频谱具有多峰值波动特性且其所处频率范围属于5hz～20hz。对照城市管道微小泄漏频谱判断准则，判断管道此时处于微小泄漏。
[0155]
采用相干函数提取频谱中的特征频带如图16所示。由图11可知函数在11.7hz处取最大值，取其所在频带11.5～15hz为特征频带。再使用带通滤波器提取特征频带中的泄漏信号，以特征频带的上下限为带通滤波器的上下截止频率，f0＝11.5hz和f1＝15hz，中心频率f
c
＝13.25hz。
[0156]
为进一步精确计算次声波信号时差，对管道上、下游泄漏信号x1(t)，x2(t)进行广义互时频分析，如图17所示。
[0157]
由图17可知，上、下游次声波泄漏信号的相关函数在0～200hz频率范围内取得最大值，且相关函数最大值对应的采样时间n＝840μs，信号采样频率f
s
＝27hz，根据该采样点数精确计算出次声波信号传递至管道上、下游两端的时差：τ0＝840
×
27
×
10
‑6＝0.0227s。
[0158]
通过压力监测平台调取管道上各压力传感器采集的压力数据，如图18所示。然后取管道采集状态30s内各测压位置处压力的平均值作为管道压力值，即：
[0159][0160]
将管道压力值代入次声波速度计算公式得：
[0161][0162]
将精确的时间延迟和次声波速度代入定位公式对泄漏点精确定位。
[0163][0164]
计算定位误差为：
[0165][0166]
与未经处理的上、下游泄漏信号粗定位结果进行对比。上、下游次声波泄漏信号的时域波形图如图19所示。
[0167]
由图19可知，上、下游泄漏信号时域波形中出现幅值尖峰的时间分别为t
′1＝0.039s，t
′2＝0.058s，则次声波信号时差τ
′
＝t
′2‑
t
′1＝0.019s，因此根据次声波定位原理得到粗定位结果：l
′1＝(12.8+0.019
×
340)/2＝9.63m。
[0168]
计算定位误差为：
[0169][0170]
计算结果表明，本发明方法的定位精度比粗定位精度提高了3.23％。
[0171]
实施案例1～3表明：本发明方法能较大程度地提高管道泄漏定位精度，可实现管道泄漏的精确定位。
[0172]
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。