一种基于逆瞬态分析法的非金属管道泄漏定位方法与流程

本发明涉及非金属管道泄漏检测技术领域，特别是涉及一种基于逆瞬态分析法的非金属管道泄漏定位方法。

背景技术：

随着现代工业化的发展，人民生活水平日益提高，管道在人们的生产生活中的地位也越来越高。非金属管材具有强度高、密度小、耐腐蚀性强、绝缘性优、使用寿命长等优点，因此非金属管道的使用成为研究的热点。但随着城市埋地管网越来越密集，尤其是输送天然气、成品油等管道建设项目越来越多，伴随着这些管道数量的增多和运行时间的增长，管道老化、泄漏问题日益凸显。在实际应用中存在管道接头处渗漏拉断、管体破裂、管道内部内衬管收缩变形等问题。因此对于非金属管道泄漏检测技术也日益兴起。

逆向瞬态分析(ita，inverse-transientanalysis)，最初由国外的pudar和liggett在1992年提出，他们指出逆向瞬态分析方法意味着流量，压力等系统状态是已知或者被测量得到，但如管道粗糙度，泄漏位置是未知的，利用列文伯格-马夸尔特算法(levenberg-marquard，lm)不断优化求解，使监测和计算压力水头差的平方和最小，从而减小误差，但是同时带来的问题就是计算速度大大减少。专利cn201810772672基于逆瞬态模型的城市非金属管道泄漏定位方法应用于输液管道的泄漏检测，应用的算法是全局搜索法pso和局部算法lm算法相结合的优化算法，来提高目标函数收敛性和计算精度。此方法的原理简单容易实现。但使用该专利方法实施计算速度缓慢，而且计算精度不够高。

本发明将蚁狮算法(alo)应用到逆瞬态检测管道泄漏的中求取收敛性最好的目标函数值摩阻系数λ，优于其他算法，能更好地提高计算速度，并提高了计算的精度。另外专利cn201810772672基于逆瞬态模型的城市非金属管道泄漏定位方法针对输液管道的泄漏检测，通过优化目标函数值jk和τk值代入蠕变函数方程，再将蠕变函数带入进moc特征方程中，从而得到h值。而本发明针对气体管道提出通过优化摩阻系数，来带入特征方程中，求得p值，因此本专利适用于城市燃气管道的泄漏检测。与其他泄漏检测方法相比，只需测量压力和流量参数即可，操作简单且快捷。

技术实现要素：

本发明所要解决的技术问题是：在现有的技术基础上，既保证目标函数的收敛速度和计算精度，提高了计算速度，本发明提供一种基于逆瞬态分析法的非金属管道泄漏定位方法。

本发明解决其技术问题所要采用的技术方案是：一种基于逆瞬态分析法的非金属管道泄漏定位方法，包括如下步骤：

s1：在实验室内搭建非金属管道模型，进行非金属管道未泄漏和泄漏两种状态的模拟试验，使用便携式无线油气管网健康诊断系统分别在未泄漏和泄漏两种状态下测量管道流量、温度及压力参数数据，找出泄漏与未泄漏的差异，总结非金属气体管道泄漏和未泄漏时的泄漏规律，从而得到泄漏管道和未泄漏管段的统计特性，其中未泄漏管段统计特性为该测量点压力不变，流量也不变，有泄漏管段的统计特性表现为当管道某点发生泄漏时泄漏点上游流量突然上升，泄漏点下游流量慢慢下降，但下降幅度不大，上下游流量差值逐渐增大，随着管道内气体趋于稳定时，流量差值慢慢接近稳定值。泄漏点上下游压力慢慢减小，渐渐趋于稳定。并获得实验管道测量到的各泄漏节点的测压头。

s2：根据s1得到的压力、流量和温度的变化规律，建立非金属管道逆瞬态模型：

s2.1拟用多阶段参数约束逆瞬态法中使用的目标函数定义为：

其中，of是目标函数，nms是测量节点数，nts是时间步数，是第i个测量站和第j个时间步长的测量压力值，pi，j(a)是根据校准模型计算得到的第i个测量站和第j个时间步长的压力值；a＝[a1，……an]^t是未知压力值，ak是第k段位置压力值，其中k＝1，2，3……n，n表示管段个数；在这种情况下，对压力信号的搜索通常以ak∈[amin，n，amax，n]，当n为1时则[amin，1，amax，1]由最大和最小可行参数值确定，[amin，n，amax，n]阶段中的搜索空间在阶段n-1中确定，并在下一节，下一阶段的搜索空间更新中解释；其中amin，n是压力信号最小可行参数值，amax，n是压力信号最大可行参数值。

对于m个不同独立逆瞬态模型运行的所有解决方案都保存在解决方案矩阵an中：

式中，ai，j，n为第j流域的第i次压力信号估算；n为管段数；m为独立逆瞬态模型运行数；n为阶段数。

本发明使用的目标函数以多阶段参数约束ita算法的形式，可以测量同一测量点的不同时间，以及不同测量点的不同时间的目标函数，全方面地得到所有情况下的目标函数，相比于在2010年1月力学学报第42卷第1期由李俊花等人发表的《一种新的长输管道泄漏点定位方法》使用的目标函数：来说更加全面精确，更适用于复杂的天然气管网泄漏检测。

s2.2泄漏与未泄漏管段分类

对于矩阵an给出的压力信号的数据集，对于给定的一组压力信号估计值：{a1，k，n....am，k，n}，m为独立逆瞬态模型运行数，表示在k管段，第n阶段，所有压力信号估计值作为一组数据集合，总结该组的统计特性，确定该组的统计特性是否和未泄漏管段的压力信号估计的统计特性相一致，如果第k管段的统计特性与未泄漏管段的特性一致，则该管段归类为未泄漏；如果不一致，则该管段归类为泄漏；

s2.3.1：根据所属的未泄漏管段集合c0，n和泄漏管段集合ca，n，为不同的区域分配不同的搜索空间间隔；泄漏管段的检测是参数估计的重点，因此，一个被划分为泄漏管段(ak，n，k∈ca，n)将保留原来的宽压力信号搜索区间{amin，1，amax，1}，以便估计策略仍然能够在全搜索空间内搜索；

s2.3.2：对一个划分为未泄漏管段(ak，n，k∈c0，n)分配一个较窄的压力信号搜索区间{amin，n+1，amax，n+1}；新的搜索空间边界amin，n+1和amax，n+1，用于分类为未泄漏的范围，由a0，n的百分位值确定：

amin，n+1＝p-thpercentileofa0，n(3)

amax，n+1＝q-thpercentileofa0，n(4)

whereq＞p；a0，n＝{ai，j，n：i∈c0，n，j＝1，...，m}(5)

a0，n是划分为未泄漏的所有管段的多个估计值的集合，因此a0，n的统计数据可以用来表示未泄漏管段的压力信号值所处的参数范围；考虑到在a0，n，p和q百分位数的估计中可能存在异常值，n被用来确定新的搜索空间界限。与原来的搜索空间{amin，1，amax，1}相比，新的搜索空间amin，n+1，amax，n+1的间隔更窄。数值验证部分将讨论百分位等级p和q的选择。

s2.4进行泄漏与未泄漏管段搜索空间更新的终止标准

当第n阶段的更新搜索空间内生成的m多个解都比第n-1阶段的先前最佳目标函数具有更大的目标函数时，该算法从第s2.2阶段的迭代终止。

s3：使用压力梯度法对泄漏管道进行初步泄漏定位，结合压力梯度法和区段分类法缩小定位管段范围。

假设气体管道进气口到漏失点的压降同漏失点到出气口处的压力降沿直线分布，通过公式可以得到泄漏点位置到进气口的距离。在计算那么管道进口到漏失点之间的压力降为：

漏失点到管道出口之间的压力降为：

联立上两式，得

其中，xl为漏失点到管道进气口的长度，m；l为管道长度值，m；hi为管道入口处的测压管水头，m；h0为管道出口处的测压管水头，m；qi为管道入口处体积流量，m³/h；q0管道出口处的体积流量，m³/h；

其中摩阻公式使用伯拉修斯公式：

雷诺系数计算公式采用经验公式：

式中，λ为摩阻系数，re为雷诺系数，qv为体积流量m³/s，d为管道内径m，v为运动粘度m²/s。

对各初步泄漏点位置设置上下游节点，分别得到上下游节点以及初步泄漏点处各时间点的压力、流量数据。

s4：在matlab的环境下使用蚁狮算法(alo)通过对目标函数进行优化，对管道的摩阻系数进行优化，在减少了目标函数的收敛时间，提高了计算精度的基础上，提高了计算速度，得到收敛性最佳的摩阻系数λ。

对于现在已提出的全局搜索法pso算法和局部算法lm算法相结合的优化算法，在实施时发现计算缓慢，以及计算精度不够，且文中没有明确定位误差相关数据。本发明在其基础上，提出使用蚁狮算法(alo)对摩阻系数进行算法优化得到收敛性好的目标函数，可以大大提高计算精度和减少收敛时间。

本发明提出使用蚁狮算法(alo)对摩阻系数进行算法优化得到收敛性好的目标函数，alo算法比pso算法相比较全局寻优能力好、收敛速度快和易于实现等方面，alo算法具有更好的全局能力和收敛速度，与lm算法相比，可以大大提高计算的速度以及保证计算的精度。

s5：将收敛性最好的摩阻系数λ代入气体特征方程，从而得出测压点计算压力pi，j值。

如s5所述的气体特征方程由联列封闭管道中一维瞬态流动的运动方程和连续性方程进行有限差分，忽略偏微分方程中的惯性项，其中对摩阻项采用二阶近似，得到c⁺、c^-特征线具体为：

其中pi-1，j-1为测压点的前一位置，前一时间点的测量压力值，mpa；pi+1，j-1为测压点后一位置，后一时间点测量压力值，mpa；pij为测压点计算压力值，mpa；mi-1，j-1为测压点的前一位置，前一时间点的流速，kg/s；mi+1，j-1为测压点后一位置，后一时间点的流速，kg/s；mij为测压点流速，kg/s；g为重力引起的加速度；a为管的横截面积，m²；a0为泄漏面积，m²；λ为摩阻系数；d为管内径，m；b为压力波速，m/s，气体管道等温流动时，b为定值。

s6：将实验测量到的测压头与计算出的测压头共同代入逆瞬态模型，得到目标函数值of。

s7：在实验所得多处的of值中，挑选最小值的位置节点，为泄漏节点。

s8：将最终得到的节点位置所求的各参数值代入压力梯度法中，得到计算泄漏点位置x′leak。

本发明的有益效果是：本发明提供的一种基于逆瞬态分析法的非金属管道泄漏定位方法，与现有技术相比，该方法只需测量管道压力、流量的参数即可，易实现。在现有的技术基础上，提高计算的速度，且保证目标函数的收敛时间和计算精度，实现了计算速度的提高。

附图说明

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

图1是本发明的流程示意图。

图2是多阶段参数约束ita算法流程图。

图3是实验管道示意图。

图4是实验管道实验节点布置图。

图5是具有特定时间间隔的特征网格示意图。

图中：1、上游涡街流量计，2、上游压力传感器，3、上游温度传感器，4、下游涡街流量计，5、下游压力传感器，6、下游温度传感器，7、上游泄漏阀，8、主管球阀，9、主管泄漏阀，10、支管球阀，11、支管泄漏阀。

具体实施方式

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

如图1所示，本发明的一种基于逆瞬态分析法的非金属管道泄漏定位方法，包括以下步骤：

s1：在实验室内搭建如图3所示全长7.9米的非金属管道模型，管道类型为地面管道，传输介质为空气，确定上游泄漏阀7模拟泄漏点i，下游泄漏阀9模拟泄漏点ii，分别距离进气口2.50m，4.15m。在管道上、下游分别安装上游涡街流量计1和下游涡街流量计4，上游压力传感器2和下游压力传感器5，上游温度传感器3和下游温度传感器6，分别测量管道上、下游管道的流量、压力和温度数据。上游流量计1距进气口1.35m，下游流量计4距进气口7.25m，上游压力传感器2距进气口2.05m，下游压力传感器5距进气口5.90m。进行非金属管道未泄漏和泄漏模拟试验，使用便携式无线油气管网健康诊断系统进行在未泄漏和泄漏两种状态的实验，并找出泄漏状态和未泄漏状态下的规律，突出泄漏状态下流量和压力的规律，并获得如图4所示实验管道各节点的压力流量值。

根据实验数据总结得到未泄漏管道压力、流量的瞬变规律和管道泄漏时压力、流量的瞬变规律。未泄漏管道压力、流量的瞬变规律为：当未发生泄漏的管道产生瞬变时，若末端阀门工作状况发生改变，管道下游先于上游管道产生瞬变流。且不论末端阀门突然开启或者突然关闭，下游流量都大于上游流量，当由瞬变状态慢慢趋于稳定时，管道上下游流量差逐渐减小直到稳定值。管道泄漏时压力、流量的瞬变规律为：当管道某点发生泄漏时泄漏点上游流量突然上升，泄漏点下游流量慢慢下降，但下降幅度不大，上下游流量差值逐渐增大，随着管道内气体趋于稳定时，流量差值慢慢接近稳定值。泄漏点上下游压力慢慢减小，渐渐趋于稳定。

以下泄漏点ii发生泄漏时测压点的实验数据，控制管道初始压力为0.3mpa，泄漏点泄漏阀门的开度为30°。

表1泄漏点ii上下游的流量、压力值

s2：根据s1得到的压力、流量和温度的变化规律，建立如图2所示的非金属管道逆瞬态模型：

s2.1拟用ita方法中使用的目标函数定义为：

其中of是目标函数，nms是测量节点数，nts是时间步数，是第i个测量站和第j个时间步长的测量压力值，pi，j(a)是根据校准模型计算得到的第i个测量站和第j个时间步长的压力值；a＝[a1，......an]^t是未知压力值，ak是第k段位置压力值，其中k＝1，2，3……n，n表示管段个数；在这种情况下，对压力信号的搜索通常以ak∈[amin，n，amax，n]，当n为1时则[amin，1，amax，1]由最大和最小可行参数值确定，[amin，n，amax，n]阶段中的搜索空间在阶段n-1中确定，并在下一节，下一阶段的搜索空间更新中解释；其中amin，n是压力信号最小可行参数值，amax，n是压力信号最大可行参数值。

s2.2泄漏与未泄漏管段分类

如图4所示为实验管道各节点图，设置x1、x2、x3、x4为管道的节点，x1为进气口处，x2距离进气口3.24m，x3距离进气口4.84m处，x4为出气口处，距离进气口7.90m。将各节点间管道确定为一个管段，如图4，分别为x1-x2管段、x2-x3管段、x3-x4管段，共三个管段，分别测量各节点的压力、流量值。

表2x1、x2、x3、x4点的流量、压力值

对管道进行区段分类，其中x1-x2管段、x3-x4管段的统计特性符合未泄漏管段的统计特性，判断为未泄漏管段，x2-x3管段的统计特性符合泄漏管段的统计特性，判断为泄漏管段。

s3：使用压力梯度法对泄漏管道进行初步泄漏定位。

假设气体管道进气口到漏失点的压降同漏失点到出气口处的压力降沿直线分布，通过公式可以得到泄漏点位置到进气口的距离。在计算那么管道进口到漏失点之间的压力降为：

漏失点到管道出口之间的压力降为：

联立上两式，得

其中，xl为漏失点到管道进气口的长度，m；l为管道长度值，m；hi为管道入口处的测压管水头，m；h0为管道出口处的测压管水头，m；qi为管道入口处体积流量，m³/h；q0管道出口处的体积流量，m³/h；

其中摩阻公式使用伯拉修斯公式：

雷诺系数计算公式采用经验公式：

式中，λ为摩阻系数，re为雷诺系数，qv为体积流量m³/s，d为管道内径m，v为运动粘度m²/s

管道内径为0.0456m，24°时空气的动力粘度1.83×10^-5pa.s运动粘度为1.43640722135×10^-5m²/s

将代入中去，再将ki、k0值代入公式：

中去，得到初始泄漏点到管道进口的距离，计算误差值，结果如下表3所示。

表3压力梯度法初步定位

由结果可知泄漏前以及刚开始泄漏时，用压力梯度法所得到的初始泄漏定位数据是极大的，远超过实际实验管道，应该舍去，因此舍去1-6组数据，而7、8组数据因为所得数值超过实际管道长度，因此也舍去，上表第9组数据开始为有效数据。通过误差值的计算，发现误差范围从(3.96％～45.29％)，可以看出直接使用压力梯度法，误差范围很大，很难确定真正的误差范围，在之后的逆瞬态法处理时有很大的麻烦。

因此泄漏与未泄漏的管道分类，已知泄漏点在x2和x3节点之间，因此算得初始泄漏点距离进气口距离应在3.24m到4.84m之间。满足该条件的只有4.6722m、3.8370m、3.5278m，对应的误差值为6.6100％、3.9615％、7.8758％。该方法极大的方便后文的数据处理。

提取出该三组数据绘制如下表格4

表4压力梯度法结合区段分类初步定位

s4：在matlab的环境下使用蚁狮算法(alo)通过对目标函数进行优化，对管道的摩阻系数进行优化，在减少了目标函数的收敛时间，提高了计算精度的基础上，提高了计算速度，得到如表5收敛性最佳的水力摩阻系数λ。

表5优化后的摩阻系数

s5：将收敛性最好的摩阻系数λ，代入气体特征方程，从而得出测压点计算压力pi，j值。

如s5所述的气体特征方程由联列封闭管道中一维瞬态流动的运动方程和连续性方程进行有限差分，利用边界条件法对数据进行处理，如图5所示具有特定时间间隔的特征网格示意图，忽略偏微分方程中的惯性项，其中对摩阻项采用二阶近似，得到c⁺、c^-特征线具体为：

其中pi-1，j-1为测压点的前一位置，前一时间点的测量压力值，mpa，pi+1，j-1为测压点后一位置，后一时间点测量压力值，mpa，pij为测压点计算压力值，mpa，mi-1，j-1为测压点的前一位置，前一时间点的流速，kg/s，mi+1，j-1为测压点后一位置，后一时间点的流速，kg/s，mij为测压点流速，kg/s；g为重力引起的加速度；a为管的横截面积，m²，a0为泄漏面积，m²；λ为摩阻系数；d为管内径，m；b为压力波速，m/s，气体管道等温流动时，b为定值。

压力波速约等于空气中的声速：

时间步长δt＝0.001s，δx＝0.34m，直径d＝0.0456m

截面积公式：

将步骤s3所得三个位置设为初步泄漏节点，在各泄漏节点位置安装压力表和流量计，测得各测压节点测得的压力、流量值以及各测压节点其所在位置的前一位置，前一时刻，后一位置后一时刻的压力、流量值。得到如表6-表8所示，将各数据代入特征线方程中进行计算。

表64.6722m所在位置的前一位置，前一时刻，后一位置后一时刻的泄漏量和压力

表73.8370m所在位置的前一位置，前一时刻，后一位置后一时刻的泄漏量和压力

表83.5278m所在位置的前一位置，前一时刻，后一位置后一时刻的泄漏量和压力

s6：将实验测量到的测压头与计算出的测压头共同代入逆瞬态模型，得到如下表9所示的目标函数值of。

表9目标函数值0f

由计算可得目标函数分别为0.088804％、0.014884％、0.138384％，其中最小值为0.014884％。

s7：在实验所得多处的of值中，挑选最小值的位置节点，即3.8370m处。

s8：将最终得到的节点位置所求的各参数值代入压力梯度法中，得到最后精确定位结果如表10所示。

表10最终定位结果及误差

由表10可知本专利提出的一种基于逆瞬态分析法的非金属管道泄漏定位方法可以进行气体非金属管道的泄漏检测及定位，定位误差最小为2.09％。

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