技术领域
[0001] 本
发明涉及城市非金属管道泄漏检测技术领域，特别是涉及一种城市非金属管道泄漏定位方法。
背景技术
[0002] 因非金属管材具有强度高、
密度小、
腐蚀性强、绝缘性优、使用寿命长等优点，其逐步取代金属管道用于
能源、原材料等的管道运输。随着城市埋地管网越来越密集，伴随着这些管道数量的增多和运行时间的增长，非金属管道因接头渗漏、拉断，或遭第三方破坏等原因导致的管道泄漏隐患或事故日益增多。因此对管道故障的有效检测和及时处理显得尤为重要。近些年，国内外在管道泄漏检测技术以及泄漏定位等方面都做了大量的研究工作。其中基于逆瞬态法的管道泄漏检测方法研究发展尤为迅速。
[0003] 逆向瞬态分析(ITA，inverse-transient analysis)，最初由国外的Pudar和Liggett在1992年提出，这种方法通常用于配
水网络，是一种用于加压管道中
摩擦系数泄漏检测和校准的强大方法，通过这种方法，启动瞬态流动并在系统某处测量压
力，通过减小测量压力与实际压力的误差以达到校准泄漏参数，从而进行泄漏检测定位。为此目的，数值模型的建立以及各种优化技术的选择成为了ITA技术应用于管道泄漏检测的主要挑战。国外学者针对配水管网提出了液压瞬态求解器以提高模型
精度减小误差，针对ITA中的优化求解采取了许多
算法，2001年Vitkovsky采用的混合复杂演化(SCE)针对大孔泄漏(泄漏孔直径与管道直径的比值大于0.1)优化计算精度较好，却不适于小孔泄漏(泄漏孔直径与管道直径的比值等于小于0.1)，2010年B.S.Jung采用的粒子群算法(PSO)和
遗传算法(GA)全局搜索能力较强，但计算量大，对目标函数的收敛性较差，对摩擦系数的校准能力较低。
发明内容
[0004] 本发明所要解决的技术问题是：为了克服
现有技术中的不足，本发明提供一种城市非金属管道泄漏定位方法。本发明通过数值模拟管道泄漏的特征或
马尔科夫法对管道是否产生泄漏作出判断标准；针对输气管道采用了气体在管道中的连续性、运动方程从而描述气体在等温管道中瞬变流动，采用序列二次规划(SQP，sequence quadratic program)法来解决非线性约束问题从而实现最小化目标函数的目的，SQP方法具有很好的兼容性和快速收敛性，计算效率高，能够对小孔泄漏的泄漏参数精确校准，从而达到对输气管道泄漏的精确定位。
[0005] 本发明解决其技术问题所要采用的技术方案是：一种城市非金属管道泄漏定位方法，该种方法主要用于检测复杂的存在弯折，众多
阀门、仪表和管支的城市管道中较大流量压力下的泄漏，整个过程主要包括泄漏判断和泄漏定位两个步骤，其中，泄漏判断是指通过
软件数值模拟规律分析或马尔科夫法两种方法判断管道泄漏；泄漏定位是指通过逆瞬态的分析方法将假定泄漏点的计算压力与测量压力差值最小化获得符合现实情况的泄漏参数来判别泄漏
位置；只有准确判断出哪个管道存在泄漏，才能对泄漏位置进行定位，通过双重的筛选提高了检测效率。
[0006] 具体包括以下步骤：
[0007] S1：判断管道是否发生泄漏；
[0008] 在对管道泄漏状态进行筛选判断时，可以采用多种方法，本发明给出两种技术方案。
[0009] 技术方案一：采用Fluent软件模拟规律分析进行泄漏判断，具体包括以下步骤：
[0010] 用AnsysIcem对管道进行建模，使用Fluent软件模拟气体在非金属管道泄漏和无泄漏状态下的压力、流量参数分布，并通过实验进行验证，分别找出气体在非金属管道泄漏和无泄漏状态下的压力、流量参数随时间和空间变化的规律。对管道是否发生泄漏给出判断标准。
[0011] 步骤A1，对实验管道进行建模，采用结构化的四面体网格划分区域，并对泄漏口附近进行适当的网格加密处理；
[0012] 步骤A2，控制方程和算法选择：
[0013] 数值模型选用标准k-ε双方程模型，在模型中，用湍能k反映特征速度，用湍能耗散率ε反映特征长度尺度，利用Boussinesq假定简化而形成了k方程和ε方程，联合上述两个方程一道组成封闭的方程组
对流动进行模拟，具体为：
[0014] k方程：
[0015] ε方程：
[0016] 其中：xi，xj，xk为坐标分量；u为速度矢量；ρ为
流体密度；μ为粘性系数；μt为湍动
粘度；Gk、Gb为速度梯度、
浮力引发的湍
动能k的产生项；各常数C1ε＝1.44，C2ε＝1.92，C3ε＝0.09，σk＝1.0，σt＝1.3；Sk、Sε为用户定义源项；YM为可压
湍流中脉动扩张项；
[0017] 模型的数值计算算法选用SIMPLE算法，对于给定的压力，求解离散形式的动量方程，得出速度场，然后对给定的压力加以修正，反复
迭代，直到获得收敛的解；
[0018] 步骤A3，边界条件设置，包括对进口和出口的条件设置、流体和管道材质的物性参数设置；
[0019] 步骤A4，在给定的多组进口压力下，分别对管道泄漏和非泄漏状态的管道内压力及流量参数的变化进行模拟，找出压力和流量参数的变化规律，对管道是否发生泄漏给出判断标准，然后将实际测得的管道数据与该标准进行对比，确定泄漏管道和无泄漏管道；
[0020] 步骤A4.1，对数值模拟泄漏和非泄漏状态的管道内压力、流量流场进行分析，找出泄漏前后整个管道内压力、流量的变化规律；
[0021] 步骤A4.2，对实验管道泄漏和非泄漏状态下的压力、流量流场进行分析，找出泄漏前后整个管道内压力、流量的变化规律；
[0022] 步骤A4.3，将数值模拟和管道实验所得到的泄漏前后管道内压力、流量的变化规律进行对比分析，总结得到判断管道产生泄漏的具体标准。
[0023] 技术方案二：采用基于马尔科夫链的方法实现非金属管道泄漏的判断，具体包括以下步骤：
[0024] 步骤B1，利用安装在管线上的流量
传感器收集第一组X个流量变化率q1，q2，q3，q4，q5...qx，其中qi为流量的变化率＝ΔQn/Qn×100％，qi(i＝1，2，3，4，5...X)中，i为采集到的X个流量变化率的排序，ΔQn为传感器传输数据的时间点Tn的流量与时间点Tn-1的流量的变化量，Qn为时间段T内的实时流量，Tn为流量计收集传输流量数据的时间，n＝1，2，3，4，5…，然后设置四个变化状态，分别为第一状态：qi的值为0，表示管道为非泄漏状态；第二状态：qi的值为0-1(％)，表示管道为小型泄漏状态；第三状态：qi的值为1-3(％)，表示管道为泄漏扩大状态；第四状态：qi的值为3-100(％)，表示管道为大泄漏状态。
[0025] 设定条件：管道的泄漏为突变过程，泄漏状态的变化为渐变构成，管道的工作状态由于其他各类因素的影响可能会直接从运行状态变为泄漏状态中某一程度的泄漏，且管道的泄漏状态不可逆，但是泄漏程度可逆，根据以上条件，得到管道泄漏的邻接矩阵A为：
[0026]
[0027] 其中，当Aij＝1表示第i状态与第j状态之间存在转移关系，当Aij＝0表示第i状态与第j状态之间不存在转移关系，(i＝1，2，3，4；j＝1，2，3，4)。
[0028] 步骤B2，进行第一组流量变化率的状态分类，将各流量变化率分别归为这四个状态分类，以时间为序，按时间顺序统计在X个参数中，相邻两个参数选取的时间间隔，流量变化率的状态转移情况，形成状态转移概率矩阵P中的各个元素；
[0029] 步骤B2.1，得到处于第一状态的流量变化率数量为n1，由第一状态至第一状态转移的数量为n11，由第一状态至第二状态转移的数量为n12，由第一状态至第三状态转移的数量为n13，由第一状态至第四状态转移的数量为n14，统计这n1个流量变化率在下一时刻的状态转移情况：
[0030] 第一状态至第一状态的转移概率
[0031] 第一状态至第二状态的转移概率
[0032] 第一状态至第三状态的转移概率
[0033] 第一状态至第四状态的转移概率
[0034] 步骤B2.2，得到处于第二状态的流量变化率数量为n2，由第二状态至第二状态转移的数量为n22，由第二状态至第三状态转移的数量为n23，由第二状态至第四状态转移的数量为n24，统计这n2个流量变化率在下一时刻的状态转移情况：
[0035] 第二状态至第一状态的转移概率p21＝0；
[0036] 第二状态至第二状态的转移概率
[0037] 第二状态至第三状态的转移概率
[0038] 第二状态至第四状态的转移概率
[0039] 步骤B2.3，得到处于第三状态的流量变化率数量为n3，由第三状态至第二状态转移的数量为n32，由第三状态至第三状态转移的数量为n33，由第三状态至第四状态转移的数量为n34，统计这n3个流量变化率在下一时刻的状态转移情况：
[0040] 第三状态至第一状态的转移概率p31＝0；
[0041] 第三状态至第二状态的转移概率
[0042] 第三状态至第三状态的转移概率
[0043] 第三状态至第四状态的转移概率
[0044] 步骤B2.4，得到处于第四状态的流量变化率数量为n4，由第四状态至第二状态转移的数量为n42，由第四状态至第三状态转移的数量为n43，由第四状态至第四状态转移的数量为n44，统计这n4个流量变化率在下一时刻的状态转移情况：
[0045] 第四状态至第一状态的转移概率p41＝0；
[0046] 第四状态至第二状态的转移概率
[0047] 第四状态至第三状态的转移概率
[0048] 第四状态至第四状态的转移概率
[0049] 则该马尔科夫链中的状态转移概率矩阵P为：
[0050]
[0051] 且满足：n1+n2+n3+n4＝X。
[0052] 步骤B3，按照步骤B1的方法收集第二组X个流量变化率，并按照步骤B2的方法统计这X个流量变化率在每个变化状态内的数量，将此时实际流量变化率在每个区间的数量编为第一参数向量，设为 其中，α1表示第一参数向量，表示在第一参数向量下，X个流量变化率中，在第一状态下的流量变化率的数量；
表示在第一参数向量下，X个流量变化率中，在第二状态下的流量变化率的数量； 表示在第一参数向量下，X个流量变化率中，在第三状态下的流量变化率的数量； 表示在第一参数向量下，X个流量变化率中，在第四状态下的流量变化率的数量；
[0053] 由马尔科夫链计算，得到将来时间段的基于马尔可夫链的预测流量变化率的参数向量β1，其公式为：
[0054]
[0055] 由矩阵乘法：
[0056]
[0057]
[0058]
[0059]
[0060] 其中， 表示在预测得出的在第一状态下的流量变化率的数量， 表示在预测得出的在第二状态下的流量变化率的数量， 表示在预测得出的在第三状态下的流量变化率的数量， 表示在预测得出的在第四状态下的流量变化率的数量；
[0061] 在这X个流量参数中，以S表示为预测状态序列号，(S＝1，2，3，4)，选取含有流量参数呈现数量最多的状态，即选取 认作此时预测管道状态将呈现为第Smax状态。
[0062] 步骤B4，按照步骤A1的方法获取收集第三组X个流量变化率的值，并按照步骤B2的方法统计其在各个变化状态的匹配数量情况，设为在这X个流量参数中，以W表示为实际管道状态序列号(W＝1，2，3，4)，选取含有流量参数呈现数量最多的状态，即选取 认作管道实际处于第Wmax状态。
[0063] 步骤B5，比较预测状态 和实际状态 中状态序列号S和W的值：
[0064] ①若S＝W＝1，则判断管道处于正常运行状态；
[0065] ②若S＜W，则管道处于开始泄漏前期阶段，且判断其为第W状态；
[0066] ③若S＞W，则判断管道处于开始泄漏后期阶段，且判断其状态为第W状态；
[0067] ④若S＝W≠1，则判断其状态为泄漏稳定阶段，且判断其状态为第W状态。
[0068] 步骤B6，以步骤B3所得数据为第一组数据，步骤B4所得数据成为第二组数据，继续预测判定管道实时状态或泄漏情况，开始重新迭代数据。步骤B1-步骤B6的泄漏检测方法可作为利用次
声波泄漏定位的辅助手段，当利用此方法检测到泄漏状态时。
[0069] 筛选出存在泄漏的管道后，进入泄漏位置的定位步骤，以确定具体的泄漏位置，便于维护和维修。泄漏位置的定位具体包括以下步骤：
[0070] S2：管道泄漏实验：获取管道压力、流量数据
[0071] 具体方式为通过逐渐关闭阀门产生瞬态管流，阀门关闭时间持续足够长，以最小化不稳定和摩擦系数不确定性的影响，得到各个测量点不同时段的压力、流量数据。
[0072] S3：建立管道泄漏控制方程判断泄漏位置
[0073] 利用气体的连续性和运动方程求解出的非线性方程与泄漏孔泄漏量方程建立非金属管道气体泄漏状态的控制方程，将S2步骤测得监测点两端Δx处的压力PA、PB，流量MA、MB带入控制方程，得到关于监测点计算压力PL与流量M1、M2的方程，以此为依据，根据逆瞬态计算原理，将某一时间点监测点a计算压力PL与试验真实数据进行比较，并采用SQP算法对两者差值进行收敛达到最小值从而获得符合现实情况的泄漏参数；
[0074] 其中，气体连续性和运动方程为：
[0075]
[0076]
[0077] 其中：M为气体流速，kg·s-1；A为管道横截面积，m2；P为压力，Pa；ρ为气体密度，kg·m-3；D为管内径，m；λ为摩阻系数。
[0078] 采用特征线法对公式(1)和(2)进行求解，得到两组特征线方程，用C+，C-进行表示，具体如：
[0079]
[0080]
[0081]
[0082]
[0083] 其中，a为压力
波速，m/s，气体管道等温流动时，a为定值。
[0084] 公式(11)和(13)分别为C+，C-特征线AP，BP；公式(10)和(12)为满足各自特征线的相容性方程。将管道均分为N等份，步长为Δx，时间步长为Δt＝aΔx，分别沿着C+特征线AP和C-特征线BP对公式(10)和(12)进行积分，其中对相容性方程的第三项即摩阻项采用一阶近似，得到C+特征线和C-特征线上的两个非线性方程，具体如：
[0085]
[0086]
[0087] 管道泄漏视为小孔泄漏，则燃气泄漏公式取决于燃气在泄漏口处的流动速度，且燃气在泄漏口处的流动速度一般为亚音速流动，具体如：
[0088]
[0089] 其中：ML为泄漏孔泄漏流量，kg·s-1；Ae＝CA0为有效泄漏面积，m2；C为孔口系数，跟泄漏孔的形状有关；A0为泄漏孔面积，m2；Pd为管道起点压力，Pa；k为气体的绝热系数，无量-1 -1纲；R为气体常数，J·kg ·K ；T为气体的
温度，K。
[0090] 假设监测
节点为泄漏点，建立泄漏控制方程对监测节点的压力、流量数据进行分析，以泄漏孔处为边界处，流入泄漏孔处之前Δx距离内的流动特性满足公式(14)，流出泄漏孔处之后Δx距离内的流动特性满足公式(15)。设流经泄漏孔处之前的流体参数下标为1，流经泄漏孔处之后的流体参数下标为2。则建立的泄漏控制方程具体如：
[0091]
[0092]
[0093]
[0094] M1-M2＝ML (20)
[0095] P1＝P2＝PL (21)
[0096] 其中，P1、M1为流入泄漏点之前的一段距离的流体压力，Pa、流量，kg·s-1；P2、M2为流出泄漏点之后的一段距离的流体压力，Pa、流量，kg·s-1；PL为泄漏点处的压力，Pa。
[0097] S4：定义目标函数
[0098] 将监测点处的压力计算值与测量值之差的最小值作为控制目标，控制方程和边界条件作为限制条件，对确定泄漏参数进行反问题分析，定义具有最小二乘目标准则的目标函数，某个时间点距离监测点a两端Δx处测量点A、B的压力PA、PB和流量MA、MB由试验给出，单位时间步长Δt后监测点a两端的计算压力P1、P2(P1＝P2＝PL)和流量M1、M2是未知的，流量之差可以判断监测点a是否产生泄漏，根据S3泄漏控制方程将已知的数据带入并且表示出监测点a两端流量M1、M2关于压力PL的方程，通过将计算压力PL与真实压力的差值收敛到最小，得到符合真实情况的监测点流量数据，并将流量数据代入燃气泄漏量方程求出有效泄漏面积，从而确定监测点泄漏情况；
[0099] S4.1：定义最小二乘准则目标函数：
[0100]
[0101] 其中：E是目标函数；M为时间步数；Pi是计算出的压力，Pa；P′i是测得的压力，Pa。通过将目标函数E最小化为零，从而产生最符合实际情况的有效泄漏面积，并根据有效面积的值确定该节点是否泄漏。
[0102] S4.2：确定泄漏有效面积的限制范围：
[0103] 0≤Aei≤Aemax (23)
[0104] 其中：Aei为节点i处的有效泄漏面积；Aemax为泄漏面积的最大限制，确定为管道横截面积的合理比例。
[0105] S4.3：摩阻系数根据实验数据由摩阻计算公式计算给出，而输气管道的流动一般处于紊流区，具体为：
[0106] 伯拉修斯公式：
[0107]
[0108] 其中，雷诺系数计算公式采用经验公式：
[0109]
[0110]
[0111] 式中，入为摩阻系数，Re为雷诺系数，Qv为体积流量m3·s-1，D为管道内径，m，υ为运动粘度，m2·s-1，μ为动力粘度，Pa·s，ρ为气体密度，kg·m-3。
[0112] S5：算法优化
[0113] 混合复杂演化(SCE)针对大孔泄漏优化计算精度较好，却不适于小孔泄漏，粒子群算法(PSO)和遗传算法(GA)全局搜索能力较强，但计算量大，对目标函数的收敛性较差，对摩擦系数的校准能力较低。SQP算法针对求解非线性优化问题具有收敛性好、计算效率高、边界搜索能力强等优势，且可以解决小孔泄漏的泄漏参数校准精确问题。为了解决目标函数的非线性规划问题，本
专利采取SQP方法在S4设定监测点两端流量限制范围内根据S3所获得压力、流量数据将目标函数最小化确定有效泄漏面积。
[0114] S5.1：建立M文件，根据实验数据定义目标函数E(x)；
[0115] S5.2：建立M文件，定义约束条件中的非线性约束Ax≤b或Aeq·X＝Beq；
[0116] S5.3：确定迭代的初始值X0；
[0117] S5.4：确定变量的上下限VLB，VUB；
[0118] S5.5：建立主程序，非线性规划求解的函数为fmincon，运行求解。
[0119] S6：泄漏定位
[0120] S6.1：在S5最小化结束时，如果发现节点的泄漏面积为非零，则将其视为泄漏；
[0121] S6.2：由于逆瞬态法的定位方法主要是判断监测点处是否泄漏，所以定位精度与监测点的设置数量有密切关系，并不能给出传统方法的定位公式，所以如果泄漏点并不在监测点处，即两个相邻的监测点的泄漏面积都不为零，则认为泄漏点在两个相邻且泄漏面积都不为零的监测点之间，并在两个监测点之间布置等距的多个测量节点，将测得的数据和计算的数据代入目标函数重复S5，直至各节点距离小于管道检测长度的4％时终止重复步骤。
[0122] 本发明的有益效果是：本发明提供的一种城市非金属管道泄漏定位方法，与现有技术相比，该方法在对非金属管道定位的同时可以确定有效泄漏面积，通过定义最小二乘准则目标函数，将反问题应用于泄漏检测，采用序列二次规划(SQP)方法确定泄漏参数，保证了目标函数的收敛速度和计算精度，从而对每个节点的泄漏有效面积的检测更加准确，快速。
附图说明
[0123] 下面结合附图和
实施例对本发明作进一步说明。
[0124] 图1是本发明的流程示意图。
[0126] 图3是具有特定时间间隔的特征网格示意图。
[0127] 图4是泄漏网格划分示意图。
[0128] 图5是两种状态下模拟管道压力分布图。
[0129] 图6是两种状态下模拟管道速度矢量分布图。
[0131] 图中：1-上游流量传感器，2-上有
压力传感器，3-上游次声波传感器，4-泄漏阀，5-下游次声波传感器，6-下游压力传感器，7-下游流量传感器。
具体实施方式
[0132] 现在结合附图对本发明作详细的说明。此图为简化的示意图，仅以示意方式说明本发明的基本结构，因此其仅显示与本发明有关的构成。
[0133] 如图2所示，本发明的实验管段采用三段直管段依次连接形成U形结构，其中，X1-X2管段长3.2m，X2-X3长1.6m，X3-X4长3.2m，泄漏孔在X5处，距离管道进口2m处，泄漏孔径为1mm。其中，X1位于管道进口处，X4位于管道出口处，X6距离进口处2.8m处，且在管道进口端依次设置上游流量传感器1、上有压力传感器2和上游次声波传感器3，在管道出口端依次设置下游流量传感器7、下有压力传感器6和下游次声波传感器5，泄漏阀4设置在X3处。
[0134] 实施例一：本实施例中采用Fluent软件模拟规律分析进行泄漏判断。
[0135] 如图1所示，本发明的一种城市非金属管道泄漏定位方法，包括使用Fluent软件模拟非金属管道泄漏和非泄漏状态下的压力、流量参数分布，研究气体在非金属管道发生泄漏的特征，给出判断管道是否发生泄漏的判断标准，并将实验数据与标准相比较判断出泄漏管道和非泄漏管道，再结合逆瞬态法对发生泄漏的管道的泄漏点进行定位。
[0136] S1：判断管道是否发生泄漏。
[0137] 步骤A1，对实验的U形管道进行建模，尺寸设置与实验装置相同，两端管段长3.2m，中间管段1.6m，管径长0.0456m。泄漏孔位于距离管道首端2m处，泄漏孔径为1mm；网格划分采用结构化的四面体网格划分区域，并对泄漏口附近进行适当的网格加密处理，泄漏口处网格加密见图4。
[0138] 步骤A2，控制方程和算法选择：
[0139] 数值模型选用标准k-ε双方程模型，在模型中，用湍能k反映特征速度，用湍能耗散率ε反映特征长度尺度，利用Boussinesq假定简化而形成了k方程和ε方程，联合上述两个方程一道组成封闭的方程组对流动进行模拟，具体为：
[0140] k方程：
[0141] ε方程：
[0142] 其中：xi，xj，xk为坐标分量；u为速度矢量；ρ为流体密度；μ为粘性系数；μt为湍动粘度；Gk、Gb为速度梯度、浮力引发的湍动能k的产生项；各常数C1ε＝1.44，C2ε＝1.92，C3ε＝0.09，σk＝1.0，σt＝1.3；Sk、Sε为用户定义源项；YM为可压湍流中脉动扩张项；
[0143] 模型的数值计算算法选用SIMPLE算法，对于给定的压力，求解离散形式的动量方程，得出速度场，然后对给定的压力加以修正，反复迭代，直到获得收敛的解；
[0144] 步骤A3，边界条件设置
[0145] 入口边界条件设置：可压缩流体的入口边界条件应设为压力入口，本试验进行了三种入口压力下的数值模拟，分别设为0.3MPa、0.2MPa和0.1MPa；
[0146] 出口边界条件设置：出口边界条件设为压力出口，并且管道末端出口和泄漏孔的末端出口都设置为外界环境压力；
[0147] 流体介质设置：空气内介质为空气，空气的物理性质见表1；
[0148] 管道材质设置：管道的材质为PE管，
比热容为0.45cal/g℃，密度为0.956g/cm3，导热系数为0.5W/m·K。
[0149] 表1常温下空气的物理性质表
[0150]
[0151] 步骤A4，在给定的多组进口压力下，分别对管道泄漏和非泄漏状态的管道内压力及流量参数的变化进行模拟，找出压力和流量参数的变化规律，对管道是否发生泄漏给出判断标准，然后将实际测得的管道数据与该标准进行对比，确定泄漏管道和无泄漏管道；
[0152] 步骤A4.1，数值模拟管道压力、流速流场分布：
[0153] 见附图5中的(a)图，在0.3MPa的进口压力下未泄漏时管道内的压力值从上游到下游总体呈线性分布，并保持着一定范围内的压力梯度随管道长度增加而减小；见附图5中的(b)图在同样的进口压力下，管道产生泄漏时，管道内总体压力下降幅度较大，最为明显的是在泄漏孔附近的压力相比未泄漏时下降了大约4.5×104Pa，并且在泄漏孔前后的压力梯度明显减小；见附图5中的(b)、(c)和(d)图，对比不同进口压力
下管道泄漏状态的压力分布图，管道进口压力越大，管道内压力的变化梯度越大，在泄漏孔附近压力普遍减小并且泄漏孔前后管段的压力变化幅度减小。
[0154] 见附图6中的a图，0.3MPa未泄漏时管道内流速趋于平稳状态，整体的变化幅度很小；见附图6中的b图，0.3MPa泄漏时管道外壁的流速明显增加，这是因为管道泄漏时，管道内外的压差
加速了管道内气体的流动，从而导致管道外侧的流速与内侧的流速的差距变大；见附图6中的(b)、(c)和(d)图，对比不同进口压力下管道泄漏状态的速度矢量分布图，可知泄漏状态，管道进口压力越大，管道内整体流速越大，内外壁的流速差异越大。
[0155] 步骤A4.2，实验管道压力、流速流场分布：
[0156] 见附图7，正方形标记的折线表示上游流量，圆形标记的折线为下游流量，正三
角标记的折线为上游压力，倒三角标记的折线为下游压力，如附图7的(a)、(b)和(c)图所示，在发生泄漏时，泄漏孔前端的流量大幅增高，泄漏孔后端的流量小幅度减小，过一段时间，泄漏孔前段流量减小，泄漏孔两端流量保持不变，且进口压力越大，流量差越大；产生泄漏时，整个管道压力减小，上下游压力差别不明显。
[0157] 步骤A4.3，数值模拟和实验结果的比较分析
[0158] 实验对模拟城市
燃气管道的不同孔径泄漏前后管道状态变化的仿真和实验比较，可以得到总体相近的结论，但也有不同的结果。相近的结论：产生泄漏时，管道内整体压力降低，且压力梯度减小，管道内气体流速变快；不同的结论：实验中，泄漏孔前段管段流速加快，后段管段流速略微降低，模拟中整个管道内外壁流速相差较大，而泄漏孔前后的流速变化不明显。综合数值模拟和实验的分析，并结合相关文献，得到管道产生泄漏时参数变化的两个标准：(1)产生泄漏时，管道内整体压力降低，且压力梯度减小，管道内气体流速变快；(2)泄漏孔前段管段流速加快，后段管段流速略微降低，泄漏孔前后的流量差会在短时间内逐渐增大然后保持不变。
[0159] S2：管道泄漏实验：获取管道压力、流量数据
[0160] 通过逐渐关闭阀门产生瞬态管流，阀门关闭时间持续足够长，以最小化不稳定和摩擦系数不确定性的影响，得到测量点X1、X2、X3、X4不同时段的压力、流量数据以及距离测量点两端Δx即1.6m处的压力、流量数据，见表1、2；
[0161] 表10.3MPa泄漏各节点流量数据
[0162]
[0163]
[0164] 表20.3MPa泄漏各节点压力数据
[0165]
[0166] S3：建立管道泄漏控制方程。
[0167] 假设节点为泄漏点，利用气体的连续性和运动方程求解出的非线性方程与泄漏孔泄漏量方程建立非金属管道气体泄漏状态的控制方程，如：
[0168]
[0169]
[0170]
[0171] M1-M2＝ML (20)
[0172] P1＝P2＝PL (21)
[0173] 其中PA、MA为节点前Δx处的压力，Pa、流量，kg·s-1；PB、MB为节点前Δx处的压力，Pa、流量，kg·s-1；P1、M1为流入泄漏点之前的一段距离的流体压力，Pa、流量，kg·s-1；P2、M2-1为流出泄漏点之后的一段距离的流体压力，Pa、流量，kg·s ；PL为泄漏点处的压力，Pa；ML为泄漏孔泄漏流量，kg·s-1；Ae＝CA0为有效泄漏面积，m2；C为孔口系数，跟泄漏孔的形状有关；A0为泄漏孔面积，m2；Pd为管道起点压力，Pa；k为气体的绝热系数，无量纲；R为气体常数，J·kg-1·K-1；T为气体的温度，K；D为管内径，m；λ为摩阻系数；A为管道横截面积，m2。
[0174] S4：建立目标函数
[0175] S4.1：根据泄漏控制方程(17-21)表示出节点计算压力与真实压力的差值，如：
[0176]
[0177] S4.2：确定泄漏面积的有效范围，从而确定泄漏量的上限，本文取泄漏孔半径为5mm时的有效泄漏面积作为泄漏面积的最大限制。
[0178] 0≤Ae≤7.85×10-5m2 (23)
[0179] S4.3：确定管道内的摩阻系数。
[0180] 针对PE管道中气体不同的流动状态，对不同管材摩阻系数λ采用不同的计算公式，本文中采用PE管道计算公式是伯拉修斯公式。
[0181] 伯拉修斯公式：
[0182]
[0183] 其中，雷诺系数计算公式采用经验公式：
[0184]
[0185]
[0186] 式中，入为摩阻系数，Re为雷诺系数，Qv为体积流量m3/s，D为管道内径，m，v为运动粘度，m2/s，μ为动力粘度，Pa·s，ρ为气体密度，kg/m3。24摄氏度
时空气的动力粘度为1.83×10-5Pa·s，运动粘度为1.4364×10-5m2/s。将未泄漏时稳态的数据及参数带入式(24-26)，求出管道内摩阻系数大约为0.01。
[0187] S5：算法优化
[0188] S5.1：将实验数据中第14至15组的节点X2的压力、流量数据及X2前后两端1.6m处的压力流量数据和已知参数代入式(22)，得到目标函数：
[0189]
[0190] S5.2：定义非线性约束
[0191] M1-M2≥0 (28)
[0192] 0≤M1MB+Δt+MLmax (29)
[0193] MB+Δt≤M2MA+Δt (30)
[0194] MLmax是由泄漏量公式(3)代进最大限制的泄漏面积求出的最大泄漏量，为0.0258kg/s，得到：
[0195]
[0196]
[0197] S5.3：确定初始值X0＝[0.00343，0.00089]
[0198] S5.4：VLB＝[0，0.00089]，VUB＝[0.0293，0.00089]
[0199] S5.5：得出结果，M1＝0.0029；M2＝0.0009，即节点X2处前后两端计算出的流量差值为0.002kg/s，计算出泄漏面积6.09×10-7m2。按照步骤S5.1至5.5相应地计算出X1，X3，X4节点前后的流量差分别为0kg/s，0kg/s，0kg/s，并计算出泄漏面积分别为0m2，0m2，0m2。
[0200] S6：泄漏定位
[0201] S6.1，经过计算，节点X1，X3，X4前后的流量差为0，泄漏面积为0，故这三个节点都未产生泄漏；节点X2前后两端产生较明显的流量差，证明在X2节点前段1.6m至后段1.6m长的管段某个位置产生了泄漏，而X3节点计算结果正常，所以判断X3节点前段1.6m至后段1.6m长的管段没有产生泄漏，所以在X2节点前1.6m处至X2节点的管段再次进行节点的设置从而进行计算。
[0202] S6.2：继续划分节点X5，X6，见附图2，并根据步骤S5.1至5.4取Δx为0.4m时相应地计算出X5，X6节点前后的流量差，分别为0.002kg/s，0kg/s，计算出泄漏面积分别为6.09×10-7m2，0m2。证明在X5前后0.4m的管段存在泄漏，可依据需要再次设立节点继续计算，本实验X5节点为泄漏点，可以证明，该方法通过对监测节点压力、流量的设置计算，可以找取产生泄漏的管段，缺点是该方法过于依赖对节点的选取，并且模型的建立，参数的计算等方面都会对结果的准确性产生影响。
[0203] 实施例二：
[0204] 本实施例与实施例一的不同之处在于步骤S1中判断管道泄漏的方法不同，本实施例中采用基于马尔科夫链的方法实现非金属管道泄漏的判断。
[0205] S1：判断管道是否发生泄漏；基于马尔科夫链的非金属管道泄漏判断包括以下步骤：
[0206] 步骤B1，利用安装在管线上的流量传感器采集到第一组五个流量：2.78、2.79、2.82、3.07、3.7，单位：m3/h，根据公式qi＝ΔQi/Qi×100％计算得到第一组五个流量变化率q1＝0％，q2＝3.58％，q3＝1.06％，q4＝8.14％，q5＝17.02％。
[0207] 步骤B2，设置四个变化状态，第一状态：流量变化率qi的值为0，表示管道为非泄漏状态；第二状态：流量变化率qi的值为0-1(％)，表示管道为小型泄漏状态；第三状态：流量变化率qi的值为1-3(％)，表示管道为泄漏扩大状态；第四状态：流量变化率qi的值为3-100(％)，表示管道为大泄漏状态；
[0208] 设定条件：管道的泄漏为突变过程，泄漏状态的变化为渐变过程，管道的工作状态由于其他各类因素的影响可能会直接从运行状态变为泄漏状态中某一程度的泄漏，且管道的泄漏状态不可逆，但是泄漏程度可逆，根据以上条件，得到管道泄漏的邻接矩阵A为：
[0209]
[0210] 其中，Aij＝1(i＝1，2...4；j＝1，2...4)表示第i状态与第j状态之间存在转移关系，Aij＝0表示第i状态与第j状态之间不存在转移关系。
[0211] 将步骤1中获得第一组流量变化率的进行状态分类，将各流量变化率分别归为这四个状态分类，以时间为序，按时间顺序统计在这五个时间段中，流量变化率的状态转移情况，形成状态转移概率矩阵P中的各个元素。
[0212] 状态转移概率矩阵P的计算过程如下：
[0213] 步骤B2.1，根据分类得到处于第一状态的流量变化率数量为1，由第一状态至第一状态转移的数量为0，由第一状态至第二状态转移的数量为0，由第一状态至第三状态转移的数量为0，由第一状态至第四状态转移的数量为0，统计这5个流量变化率在下一时刻的状态转移情况：
[0214] 第一状态至第一状态的转移概率p11＝0；
[0215] 第一状态至第二状态的转移概率p12＝0；
[0216] 第一状态至第三状态的转移概率p13＝1
[0217] 第一状态至第四状态的转移概率p14＝0；
[0218] 步骤B2.2，得到处于第二状态的流量变化率数量为0，由第二状态至第二状态转移的数量为0，由第二状态至第三状态转移的数量为0，由第二状态至第四状态转移的数量为0，统计这0个流量变化率在下一时刻的状态转移情况：
[0219] 第二状态至第一状态的转移概率p21＝0；
[0220] 第二状态至第二状态的转移概率p22＝0；
[0221] 第二状态至第三状态的转移概率p23＝0；
[0222] 第二状态至第四状态的转移概率p24＝0；
[0223] 步骤B2.3，得到处于第三状态的流量变化率数量为1，由第三状态至第二状态转移的数量为0，由第三状态至第三状态转移的数量为0，由第三状态至第四状态转移的数量为1，统计这0个流量变化率在下一时刻的状态转移情况：
[0224] 第三状态至第一状态的转移概率p31＝0；
[0225] 第三状态至第二状态的转移概率p32＝0；
[0226] 第三状态至第三状态的转移概率p33＝0；
[0227] 第三状态至第四状态的转移概率p34＝1；
[0228] 步骤B2.4，得到处于第四状态的流量变化率数量为3，由第四状态至第二状态转移的数量为0，由第四状态至第三状态转移的数量为0，由第四状态至第四状态转移的数量为2，统计这0个流量变化率在下一时刻的状态转移情况：
[0229] 第四状态至第一状态的转移概率p41＝0；
[0230] 第四状态至第二状态的转移概率p42＝0；
[0231] 第四状态至第三状态的转移概率p43＝1/3
[0232] 第四状态至第四状态的转移概率p44＝2/3
[0233] 则该马尔科夫链中的状态转移概率矩阵P为：
[0234]
[0235] 步骤B3，收集第二组流量：4.79、6.62、7.16、7.66、8.22，单位：m3/h。得到第二组五个流量变化率：q1＝22.75q2＝27.64q3＝7.54％，q4＝6.52％，q5＝6.81％，同步骤2统计这五个流量变化率在每个区间内的数量，将此时实际流量变化率在每个区间的数量编为第一参数向量，设为α1＝[0 0 0 5]，其中α1表示第一参数向量，0表示在第一参数向量下，五个流量变化率中，在第一状态下的流量变化率的数量；0表示在第一参数向量下，五个流量变化率中，在第二状态下的流量变化率的数量；0表示在第一参数向量下，五个流量变化率中，在第三状态下的流量变化率的数量；5表示在第一参数向量下，五个流量变化率中，在第四状态下的流量变化率的数量；
[0236] 由马尔科夫链计算，得到第三个时间段的基于马尔可夫链的预测流量变化率的参数向量β1，其公式为：
[0237] β1＝α1×P，得到的β1＝[0 0 0 5]，
[0238] 由β1＝[0 0 0 5]，得到此时预测下一时刻处于第四状态的参数数量最多，含有5个参数，即 其中S＝4，则此时预测下一时刻状态为第四状态，即大泄漏状态。
[0239] 步骤B4，收集到第三组流量：8.66、8.66、9.01、9.22、9.42，单位：m3/h，得到第三组五个流量变化率的值q1＝5.08％，q2＝0％，q3＝3.88％，q4＝2.33％，q5＝2.12％，同步骤2统计其在各个区间的匹配数量情况，设为α2＝[1 021]，其中得到此时处于第三状态的参数数量最多，含有2个参数，即 其中W＝2，即此时实际状态为第三状态。
[0240] 步骤B5，比较 和 中S和W的数量值：
[0241] 由于S＝4，W＝2，S＞W，得到处于开始泄漏后期阶段，且判断其状态为第三状态。
[0242] 步骤B6，使第二组数据为第一组数据，第三组数据成为第二组数据，继续预测判定管道实时状态或泄漏情况。此时为泄漏第三状态，判定发生泄漏，触发报警程序，系统开始进行报警。
[0243] 本例对泄漏定位的计算，即步骤S2至S6同实施例一，不再重复赘述说明。
[0244] 总结：基于软件模拟的方法和马尔可夫法对管道的泄漏判断都有其优势和缺点，软件模拟的方法只能总结出泄漏和非泄漏两种状态流畅分布的规律，优点是对两种状态流场分布规律的分析，可以为泄漏的判断提供理论思路，给出理论依据，但在具体实践中，很难量化到具体数值，来判断管道是否发生微小泄漏；而马尔可夫法可以根据流量参数的捕捉对管道的泄漏情况进行判断，但是可能会由于其他因素对流量的影响从而产生误判，需要很多数据的
支撑来保证对泄漏判断的准确性。
[0245] 以上述依据本发明的理想实施例为启示，通过上述的说明内容，相关的工作人员完全可以在不偏离本发明的范围内，进行多样的变更以及
修改。本项发明的技术范围并不局限于
说明书上的内容，必须要根据
权利要求范围来确定其技术性范围。
