CN110375206A - 基于glt及opt方法的供热管网运行状态监测方法 - Google Patents

Abstract

基于GLT及OPT方法的供热管网运行状态监测方法，步骤如下：1)将监测终端安装于城市供热管网的薄弱环节；2)安装监测终端时使安装距离＜200m或保证GPRS通信信号；使用GPRS上报数据；当GPRS信号受影响，采用RF射频通信将数据发送到相邻的终端由GPRS信号良好的终端发送到上位机；3)运行上位机，在上位机进行第一轮数据接收完成后，及上位机以及接收到管网中全部终端发送的数据后，为每一个终端标定位置，为监测系统的搭建做准备；4)监测系统投入运行后，由上位机通过GLT方法或者OPT方法将温度数据以及电导率数据作为输入；5)泄漏报警检测完成，将数据加入数据库，单次接收结束；该方法优于普通的阈值设定法，提高了方法的准确性。

Description

基于GLT及OPT方法的供热管网运行状态监测方法

技术领域

本发明属于城市供热管网泄漏监测技术领域，具体涉及基于GLT及OPT方法的供热管网运行状态监测方法。

背景技术

随着城市现代化水平的提高，冬季城市供热由过去的分散供热逐渐转变为集中供热。该方式是北方城镇的主要供热方式。集中供热具有污染低、能耗小、效率高等优点。但由于集中供热传输距离远，管道网络错综复杂等原因，如果管网被损坏发生泄漏，轻则影响热能生产和生活取暖，重则致使供热中断事故，不仅导致巨大的能源浪费，而且会造成严重的经济损失，最重要的是在部分地区，冬季严寒，供热事故将严重影响人们的正常生活。

据统计，近几年1400起供热事故中，有80％是由于热网补偿器泄漏导致的。但是热网补偿器泄漏过程极为隐蔽，人工巡检的方式不现实，需要耗费大量的人力物力。因此需要使用一种有效的泄漏监测手段，能够让工作人员在办公室甚至在办公室外就能随时随地监测管网状态，及时高效地发现异常状态而且定位到泄漏位置，方便检修。

解决管网泄漏报警的方法在理论上有很多，例如：神经网络以及模拟退火等方法，虽然能够从理论上解决问题，但是方法必须需要大量历史数据，对于实际工程问题，最终的决策方法要求简便实用，能够在微控制器中使用，切实解决实际问题，但现有的简单的阈值法缺乏理论依据，阈值的选择没有一个可靠的理论支撑，所以研究的实时泄漏报警方法及GLT泄漏报警方法以及OPT泄漏报警方法显得尤为重要。

发明内容

为克服上述现有技术的不足，本发明的目的在于提供基于GLT及OPT方法的供热管网运行状态监测方法，基于GLT方法的泄漏报警方法以及基于OPT方法的泄漏报警方法，以解决现有复杂方法的局限性，实用性不好的问题，以及简单方法没有理论支持的问题，提高泄漏报警决策的实用性和可靠性。

为实现上述目的，本发明采用的技术方案是，基于GLT及OPT方法的供热管网运行状态监测方法，包括以下步骤：

步骤1，将监测终端安装于城市供热管网的薄弱环节，一个监测终端包括四路温度数据、两路电导率数据，同时支持一路压力数据以及一路流量数据；其中四路温度数据以及两路电导率数据在使用GLT方法或者OPT方法时必须使用；对于检测同一数据类型的传感器，使其安装位置为检测同一目标而位于不同位置维度，电导率传感器的安装类似；

步骤2，安装监测终端时使安装距离＜200m或保证GPRS通信信号；在GPRS信号良好时，使用GPRS上报数据；由于安装环境问题，当GPRS信号受影响，采用RF射频通信将数据发送到相邻的终端由GPRS信号良好的终端发送到上位机；

步骤3，监测终端安装完成后，运行上位机，在上位机进行第一轮数据接收完成后，及上位机以及接收到管网中全部终端发送的数据后，在上位机运行地图展示，为每一个终端标定位置，为监测系统的搭建做准备；

步骤4，监测系统投入运行后，每次接收到数据后，都由上位机通过GLT方法或者OPT方法将温度数据以及电导率数据作为输入，GLT方法的泄漏检测函数为：

式中，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵，σ²为噪声方差，p为输入Z的线性变换，p'为p的转置，

p＝V'Z (2)

OPT的泄漏检测函数为：

ρ_k ^*为k位置传感器检测值相应的最优奇偶残差，|ρ_k ^*|为其绝对值，T为通过正态分布以及给定的虚警率确定的阈值，如果输出大于阈值，则发出泄漏报警；除此之外，两个方法还能对每一个传感器的泄漏进行判断，得出最有可能发生泄漏的传感器位置；

步骤5，泄漏报警检测完成，将数据加入数据库，单次接收结束。

步骤1所述的GLT方法，将m个余度温度传感器或电导率传感器的量测方程为：

Z＝HX+f+ξ (4)

式中,X∈Rⁿ是n个待量测的状态变量，这些待测变量可以是土壤的温度及电导率，Z∈R^m是m个传感器的量测值，方法要求m≥n；H是传感器配置的测量矩阵，也可称为安装矩阵，主要用来描述从待测量到传感器数据的过程；f为泄漏向量，ξ为噪声向量；若待测量为温度一个量，即n＝1；有四路温度检测，即m＝4；显然温度传感器直接测量温度的大小，即：

H＝[1 1 1 1]' (5)

f是泄漏向量，描述泄漏发生的位置，ξ用来描述传感器在测量时产生的误差，使用m维零均值、协方差为σ²的高斯白噪声序列描述，即E(ξ)＝0,E(ξξ^T)＝σ²I_m，I_m为单位矩阵，由方差定义：

式中，X表示变量，此处为温度传感器的输出温度，μ为总体均值，这里使用温度实际值，若传感器测量误差为1℃，则易得(X-μ)²＝1，因此可得σ²＝1。

步骤1所述的GLT方法，是通过提出有泄漏状态和无泄漏状态两种假设，采用基于奇偶空间的广义似然比检验方法来进行泄漏的报警；对于泄漏检测，使用对泄漏解耦的奇偶残差实现泄漏检测函数，考虑量测值Z的线性变换：

p＝V'Z＝V'HX+V'f+V'ξ (7)

式中，f为泄漏向量，ξ为噪声向量，H为安装矩阵，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵；

让p对泄漏位置解耦，需要V满足能够分离出传感器测量的目标值与泄漏位置，即如下约束条件：

V'H＝0 (8)

带入式(8)，此时待观测的温度由于乘了一个0阵即式(7)化简为:

p＝V'f+V'ξ (9)

从式(9)可知，p独立于待测状态X而仅与噪声或可能的泄漏故障有关，实际上，当无泄漏时，p仅是噪声的函数，即:

p＝V'ξ (10)

接下来规定泄漏报警决策的假设检验：

无泄漏状态H₀：E(p)＝0,E(pp')＝σ²V'V；

有泄漏状态H₁：E(p)＝μ≠0,E((p-μ)(p-μ)')＝σ²V'V；其中μ＝V'f，

由于p是ξ的线性函数，因此p服从与ξ一致的m-n维的正态分布，在无泄漏时似然函数为式(11)：

式中，ρ(p|H₀)表示对无泄漏状态H₀构建的似然函数，K为与奇偶向量与噪声方差相关的比例系数，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵，σ²为噪声方差，p为输入Z的线性变换，p'为p的转置；

在无泄漏时似然函数为式(12)

式中，ρ(p|H₁)表示对有泄漏状态H₁构建的似然函数，其中μ＝V'f，其余参数与式(11)一致。

在式(11)与式(12)中，比例系数K为：

联立式(11)与式(12)构造对数似然比函数为：

为了求出最大值可对式(14)关于μ求导得出极大似然函数值为：

所以，可以用p搭建泄漏报警函数

则通过p服从的正态分布可知FD_GLT～χ²(m-n)，

于是泄漏报警的准则为：FD_GLT≥T时，说明泄漏事故发生，触发报警；否则，判决无泄漏；

阈值T的选择由式(16)推得可知泄漏报警函数符合自由度为m-n的卡方分布即服从中心χ²(m-n)分布，该分布的概率密度函数为

假设给定一个阈值T，则可以通过式(17)计算出虚警率P_F:

显然若要求给定的报警准确率为1-α，显然α为虚警率，则只需查χ²分布表，得1-α分位数就可求出需要使用的阈值

步骤1所述的GLT方法，还可以定位出最有可能发生泄漏的位置，即对在膨胀节一圈安装的温度传感器，哪一个传感器所在的位置最有可能发生了泄漏；对于泄漏定位，提出泄漏定位函数，该函数是假设各个传感器位置发生泄漏时，使对应的奇偶残差的似然函数取得最大时的统计量；

泄漏位置定位包括m＝4个假设检验，假设H_i为第i个温度传感器处发生了泄漏，i＝1、2、...、m，：E(p)＝μ_i≠0，E((p-μ_i)'(p-μ_i))＝σ²其中μ_i＝V'e_if_i，e_i为第i行值，i＝1、2、...、m，为1，其余值为0的m维列向量，这里为4维列向量；记V'e_i＝V^*为奇偶矩阵V的第i个行向量的转置向量，

则相应于第i个温度传感器检测数据的对数似然函数为：

式中，表示相应于第i个温度传感器检测数据的对数似然函数，K为与奇偶向量与噪声方差相关的比例系数，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵，V_i ^*为奇偶矩阵V的第i行构成向量的转置向量，σ²为噪声方差，p为输入Z的线性变换，p'为p的转置，f_i为第i个温度传感器的泄漏向量，且μ_i＝V'e_if_i，对式(18)关于f_i求导求出对于f_i的极大似然估计:

极大似然估计值为：

式中，表示相应于第i个温度传感器检测数据的对数似然函数的极大似然估计值，K为与奇偶向量与噪声方差相关的比例系数，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵，V_i ^*为奇偶矩阵V的第i行构成向量的转置向量，σ²为噪声方差，p为输入Z的线性变换，p'为p的转置，如果设计时令V'V＝I_m-n，I_m-n为m-n维单位矩阵，由上式(20)可确定针对各个位置上的温度传感器的泄漏判断函数为：

式中，FI_GLT(i)表示第i个温度传感器的泄漏判断函数值，K为与奇偶向量与噪声方差相关的比例系数，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵，V_i ^*为奇偶矩阵V的第i行构成向量的转置向量，σ²为噪声方差，p为输入Z的线性变换，p'为p的转置，在判断泄漏时FI_GLT(i)越大，说明在该位置发生泄漏的可能性越大，对应于所得的极大似然估计值就越大；

如果第k个温度传感器检测的温度值经过计算的泄漏位置定位函数值FI_GLT(k)是所有检测到的温度计算出的函数值中的最值，m为传感器个数，则:

则说明第k个温度传感器所在的位置最有可能已发生了泄漏。

步骤1所述的OPT方法将泄漏检测与定位同时进行，

对于式(1)，增加D表示泄漏输入矩阵，F噪声输入矩阵则：

Z＝HX+Df+Fξ (23)

对第i个温度传感器所在的位置发生泄漏事故，选取最优奇偶向量v_i，带入式(23)得m＝4个奇偶方程：

r_i＝v_i'Z＝v_i'Df+v_i'Fξ,i＝1,2,…,m (24)

式中，r_i为与第i个温度传感器所在的位置发生泄漏事故这一事件对应的残差；系统发生泄漏时，可以通过下式判断最有可能发生泄漏的位置，式(25)的情况说明最有可能发生泄漏事故的位置在温度传感器i的位置上；

|r_i|＞|r_j|,j＝1,2,…,m,j≠i (25)

为了设计r，建立性能指标函数：

其中，v′_iDe_i表示对第i个温度传感器故障的敏感度，v′_iDe_j表示对第j个温度传感器故障的敏感度，e_i为m阶单位矩阵I_m的第i个列向量，v_i为设计的向量，||v′_iF||表示对此处使用的温度传感器测量误差的敏感程度；

使S_i最大，是为了让函数对在温度传感器i位置上发生的泄漏事故敏感，对其它位置上的泄漏事故以及测量时产生的误差不敏感，即在第i个温度传感器位置上发生泄漏事故时使|r_i|尽量大，

考虑到：

带入式(26)的的极值问题整理成求：

其中，u＝V′De_i，M_B＝V′(FF′+DD′-De_ie_i′D′)V，为了求得式(26)的极值解，由柯西—士瓦兹不等式定理及其相关推论得：

式(26)在出取得最大值，于是：

其中，u＝V′De_i，M_B＝V′(FF′+DD′-De_ie_i′D′)V，

最后，将v_i标准化，相应的最优奇偶残差为：

标准化以后的奇偶残差ρ_i～N(0,1)，则

最后说明T的选择，对给定的要求的准确率1-α，显然α为虚警率，通过正态分布求得阈值T＝u_1-α/2。

步骤1所述的GLT方法以及OPT方法中最优奇偶阵需要满足以下条件：

式中，V阵可以采用Potter方法求取，具体如下：

Step1，通过给定的H阵，通过w＝I-Hinv(H'H)H'计算得到w阵；式中H表示安装矩阵，H'为其转置，I为m维单位矩阵；

Step2，根据求得V阵首元素，

v_1j＝w_1j/v₁₁,j＝2,...,m；

v_ij＝0,i＝2,3,...,m-n,j＝1,2,...,i；i＞j；更新第一行除首元素外全部元素以及剩余下三角全部元素为0；式中，w_ij表示w阵第i行第j列的元素，v_ij表示V阵第i行第j列的元素；

Step3，从i＝2开始更新对角元素，w_ii表示w阵第i行第i列的元素，v_ii表示V阵第i行第i列的元素，

然后执行Step4；

Step4，w_ij表示w阵第i行第j列的元素，v_ij表示V阵第i行第j列的元素，利用由Step3使用的i更新该行没有计算的其余元素；

Step5，i增加到下一行，调到Step3继续执行，指导更新得到完整的V阵。

本发明的有益效果是：

本发明的关键方法包括基于奇偶空间的广义似然比检验方法(GLT)的泄漏报警方法以及基于奇偶空间的最优奇偶向量检验方法(OPT)的泄漏报警方法。系统由多个监测终端向上位机发送监测数据，由上位机通过GLT或者OPT方法进行数据处理，对可能出现的泄漏情况进行报警。该系统使用GPRS+RF通讯方式提高了终端与上位机通信的可靠性，采用GLT以及OPT方法在降低传感器测量噪声的影响下，通过判断检测函数的输出值大小确定泄漏事故的发送与否，阈值是以设定的误警率确定的。该方法优于普通的阈值设定法，即通过简单的经验值确定异常正常的界限的方法，提高了方法的准确性，而且方法的虚警率是可以事先由方法给定的，在理论上取用了在特定噪声下给定虚警率的检测模型。

可以广泛用于城市供热管网泄漏监测与报警，而且对整个城市的管网运行状态进行有效的监测，包括系统自运行开始的全部历史数据，通过多种功能方便地查询数据。本发明采用的基于GLT方法的泄漏报警方法以及基于OPT方法的泄漏报警方法，简化了复杂方法，而且实现了任意给定误报率的泄漏报警方法，有效减少了泄漏诊断中的漏报、误报现象。监测终端的通讯采用GPRS+RF通讯方式提高了终端与上位机通信的可靠性。

附图说明

图1为本发明系统总体原理图。

图2为本发明的流程图。

图3为本发明具体实施时模拟的测试数据。

图4(a)为GLT方法泄漏报警函数输出图。

图4(b)为OPT方法泄漏报警函数输出图。

图4(c)为GLT方法泄漏传感器的定位结果图。

图4(d)为OPT方法泄漏传感器的定位结果图

具体实施方式

下面结合附图和具体实施例对本发明作进一步详细阐述，但本发明并不限于这些实施方式。

图1为本发明系统总体原理图，按照该原理图搭建系统，图中有三部分，监测终端、数据处理中心以及APP。系统有多个ID互不相同的监测终端安装于城市供热管网的薄弱环节，有数据处理中心通过GPRS以及RF射频与数据处理中心通讯，由数据处理中心通过泄漏报警方法处理后发出相应警报；APP通过向数据处理中心发出请求可远程获取数据处理中心的历史数据实现对管网系统的远程监测。对本发明的基于GLT方法以及OPT方法的供热管网运行状态监测系统，流程图如图2所示，具体按照以下步骤实施：

步骤1，将监测终端安装于城市供热管网的薄弱环节，一个监测终端包括四路温度数据、两路电导率数据，同时支持一路压力数据以及一路流量数据；其中四路温度数据以及两路电导率数据在使用GLT方法或者OPT方法时必须使用；对于检测同一数据类型的传感器，使其安装位置为检测同一目标而位于不同位置维度，电导率传感器的安装类似；

步骤2，安装监测终端时使安装距离＜200m或保证GPRS通信信号；在GPRS信号良好时，使用GPRS上报数据；由于安装环境问题，当GPRS信号受影响，采用RF射频通信将数据发送到相邻的终端由GPRS信号良好的终端发送到上位机；

步骤3，监测终端安装完成后，运行上位机，在上位机进行第一轮数据接收完成后，及上位机以及接收到管网中全部终端发送的数据后，在上位机运行地图展示，为每一个终端标定位置，为监测系统的搭建做准备；

步骤4，监测系统投入运行后，每次接收到数据后，都由上位机通过GLT方法或者OPT方法将温度数据以及电导率数据作为输入，GLT方法的泄漏检测函数为：

式中，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵，σ²为噪声方差，p为输入Z的线性变换，p'为p的转置，

p＝V'Z (2)

OPT的泄漏检测函数为：

ρ_k*为k位置传感器检测值相应的最优奇偶残差，|ρ_k ^*|为其绝对值，T为通过正态分布以及给定的虚警率确定的阈值，如果输出大于阈值，则发出泄漏报警；除此之外，两个方法还能对每一个传感器的泄漏进行判断，得出最有可能发生泄漏的传感器位置；

步骤5，泄漏报警检测完成，将数据加入数据库，单次接收结束。

步骤1所述的GLT方法，将m个余度温度传感器或电导率传感器的量测方程为：

Z＝HX+f+ξ (4)

式中,X∈Rⁿ是n个待量测的状态变量，这些待测变量可以是土壤的温度及电导率，Z∈R^m是m个传感器的量测值，方法要求m≥n；H是传感器配置的测量矩阵，也可称为安装矩阵，主要用来描述从待测量到传感器数据的过程；f为泄漏向量，ξ为噪声向量；若待测量为温度一个量，即n＝1；有四路温度检测，即m＝4；显然温度传感器直接测量温度的大小，即：

H＝[1 1 1 1]' (5)

f是泄漏向量，描述泄漏发生的位置，ξ用来描述传感器在测量时产生的误差，使用m维零均值、协方差为σ²的高斯白噪声序列描述，即E(ξ)＝0,E(ξξ^T)＝σ²I_m，I_m为单位矩阵，由方差定义：

式中，X表示变量，此处为温度传感器的输出温度，μ为总体均值，这里使用温度实际值，若传感器测量误差为1℃，则易得(X-μ)²＝1，因此可得σ²＝1。

步骤1所述的GLT方法，是通过提出有泄漏状态和无泄漏状态两种假设，采用基于奇偶空间的广义似然比检验方法来进行泄漏的报警；对于泄漏检测，使用对泄漏解耦的奇偶残差实现泄漏检测函数，考虑量测值Z的线性变换：

p＝V'Z＝V'HX+V'f+V'ξ (7)

式中，f为泄漏向量，ξ为噪声向量，H为安装矩阵，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵；

让p对泄漏位置解耦，需要V满足能够分离出传感器测量的目标值与泄漏位置，即如下约束条件：

V'H＝0 (8)

带入式(8)，此时待观测的温度由于乘了一个0阵即式(7)化简为:

p＝V'f+V'ξ (9)

从式(9)可知，p独立于待测状态X而仅与噪声或可能的泄漏故障有关，实际上，当无泄漏时，p仅是噪声的函数，即:

p＝V'ξ (10)

接下来规定泄漏报警决策的假设检验：

无泄漏状态H₀：E(p)＝0,E(pp')＝σ²V'V；

有泄漏状态H₁：E(p)＝μ≠0,E((p-μ)(p-μ)')＝σ²V'V；其中μ＝V'f，

由于p是ξ的线性函数，因此p服从与ξ一致的m-n维的正态分布，在无泄漏时似然函数为式(11)：

式中，ρ(p|H₀)表示对无泄漏状态H₀构建的似然函数，K为与奇偶向量与噪声方差相关的比例系数，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵，σ²为噪声方差，p为输入Z的线性变换，p'为p的转置；

在无泄漏时似然函数为式(12)

式中，ρ(p|H₁)表示对有泄漏状态H₁构建的似然函数，其中μ＝V'f，其余参数与式(11)一致。

在式(11)与式(12)中，比例系数K为：

联立式(11)与式(12)构造对数似然比函数为：

为了求出最大值可对式(14)关于μ求导得出极大似然函数值为：

所以，可以用p搭建泄漏报警函数

则通过p服从的正态分布可知FD_GLT～χ²(m-n)，

于是泄漏报警的准则为：FD_GLT≥T时，说明泄漏事故发生，触发报警；否则，判决无泄漏；

阈值T的选择由式(16)推得可知泄漏报警函数符合自由度为m-n的卡方分布即服从中心χ²(m-n)分布，该分布的概率密度函数为

假设给定一个阈值T，则可以通过式(17)计算出虚警率P_F:

显然若要求给定的报警准确率为1-α，显然α为虚警率，则只需查χ²分布表，得1-α分位数就可求出需要使用的阈值

步骤1所述的GLT方法，还可以定位出最有可能发生泄漏的位置，即对在膨胀节一圈安装的温度传感器，哪一个传感器所在的位置最有可能发生了泄漏；对于泄漏定位，提出泄漏定位函数，该函数是假设各个传感器位置发生泄漏时，使对应的奇偶残差的似然函数取得最大时的统计量；

泄漏位置定位包括m＝4个假设检验，假设H_i为第i个温度传感器处发生了泄漏，i＝1、2、...、m，：E(p)＝μ_i≠0，E((p-μ_i)'(p-μ_i))＝σ²其中μ_i＝V'e_if_i，e_i为第i行值，i＝1、2、...、m，为1，其余值为0的m维列向量，这里为4维列向量；记V'e_i＝V^*为奇偶矩阵V的第i个行向量的转置向量，

则相应于第i个温度传感器检测数据的对数似然函数为：

式中，表示相应于第i个温度传感器检测数据的对数似然函数，K为与奇偶向量与噪声方差相关的比例系数，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵，V_i ^*为奇偶矩阵V的第i行构成向量的转置向量，σ²为噪声方差，p为输入Z的线性变换，p'为p的转置，f_i为第i个温度传感器的泄漏向量，且μ_i＝V'e_if_i，对式(18)关于f_i求导求出对于f_i的极大似然估计:

极大似然估计值为：

式中，表示相应于第i个温度传感器检测数据的对数似然函数的极大似然估计值，K为与奇偶向量与噪声方差相关的比例系数，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵，V_i ^*为奇偶矩阵V的第i行构成向量的转置向量，σ²为噪声方差，p为输入Z的线性变换，p'为p的转置，如果设计时令V'V＝I_m-n，I_m-n为m-n维单位矩阵，由上式(20)可确定针对各个位置上的温度传感器的泄漏判断函数为：

式中，FI_GLT(i)表示第i个温度传感器的泄漏判断函数值，K为与奇偶向量与噪声方差相关的比例系数，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵，V_i*为奇偶矩阵V的第i行构成向量的转置向量，σ²为噪声方差，p为输入Z的线性变换，p'为p的转置。在判断泄漏时FI_GLT(i)越大，说明在该位置发生泄漏的可能性越大，对应于所得的极大似然估计值就越大；

如果第k个温度传感器检测的温度值经过计算的泄漏位置定位函数值FI_GLT(k)是所有检测到的温度计算出的函数值中的最值，m为传感器个数，则:

则说明第k个温度传感器所在的位置最有可能已发生了泄漏。

步骤1所述的OPT方法将泄漏检测与定位同时进行，

对于式(1)，增加D表示泄漏输入矩阵，F噪声输入矩阵则：

Z＝HX+Df+Fξ (23)

对第i个温度传感器所在的位置发生泄漏事故，选取最优奇偶向量v_i，带入式(23)得m＝4个奇偶方程：

r_i＝v_i'Z＝v_i'Df+v_i'Fξ,i＝1,2,…,m (24)

式中，r_i为与第i个温度传感器所在的位置发生泄漏事故这一事件对应的残差；系统发生泄漏时，可以通过下式判断最有可能发生泄漏的位置，式(25)的情况说明最有可能发生泄漏事故的位置在温度传感器i的位置上；

|r_i|＞|r_j|,j＝1,2,…,m,j≠i (25)

为了设计r，建立性能指标函数：

其中，v_i′De_i表示对第i个温度传感器故障的敏感度，v_i′De_j表示对第j个温度传感器故障的敏感度，e_i为m阶单位矩阵I_m的第i个列向量，v_i为设计的向量，||v_i′F||表示对此处使用的温度传感器测量误差的敏感程度；

使S_i最大，是为了让函数对在温度传感器i位置上发生的泄漏事故敏感，对其它位置上的泄漏事故以及测量时产生的误差不敏感，即在第i个温度传感器位置上发生泄漏事故时使|r_i|尽量大，

考虑到：

带入式(26)的的极值问题整理成求：

其中，u＝V′De_i，M_B＝V′(FF′+DD′-De_ie_i′D′)V，为了求得式(26)的极值解，由柯西—士瓦兹不等式定理及其相关推论得：

式(26)在出取得最大值，于是：

其中，u＝V′De_i，M_B＝V′(FF′+DD′-De_ie_i′D′)V，

最后，将v_i标准化，相应的最优奇偶残差为：

ρ_i ^*＝v_i′＝v_i′Fε,i＝1,2,...m (31)

标准化以后的奇偶残差ρ_i～N(0,1)，则：

最后说明T的选择，对给定的要求的准确率1-α，显然α为虚警率，通过正态分布求得阈值T＝u_1-α/2。

步骤1所述的GLT方法以及OPT方法中最优奇偶阵需要满足以下条件：

式中，V阵可以采用Potter方法求取，具体如下：

Step1，通过给定的H阵，通过w＝I-Hinv(H'H)H'计算得到w阵；式中H表示安装矩阵，H'为其转置，I为m维单位矩阵；

Step2，根据求得V阵首元素，

v_1j＝w_1j/v₁₁,j＝2,...,m；

v_ij＝0,i＝2,3,...,m-n,j＝1,2,...,i；i＞j；更新第一行除首元素外全部元素以及剩余下三角全部元素为0；式中，w_ij表示w阵第i行第j列的元素，v_ij表示V阵第i行第j列的元素；

Step3，从i＝2开始更新对角元素，w_ii表示w阵第i行第i列的元素，v_ii表示V阵第i行第i列的元素，

然后执行Step4；

Step4，w_ij表示w阵第i行第j列的元素，v_ij表示V阵第i行第j列的元素，利用由Step3使用的i更新该行没有计算的其余元素；

Step5，i增加到下一行，调到Step3继续执行，指导更新得到完整的V阵。

以四路温度传感器为例，可计算得其安装矩阵为：

以四路温度传感器为例，并采用95％的准确率测试，假定由安装以及传感器造成的总噪声方差为1。首先设计模拟1000组数据，数据初始值为50℃，为了配合测量误差，为其注入均值为0方差为1的高斯白噪声序列模拟测量误差。考虑两种类型的泄漏故障，阶跃的突发性泄漏事故以及斜坡的缓慢的泄漏事故。为模拟这两种泄漏事故，对温度1在1到100个数据上增加20℃的阶跃泄漏，对温度2在101到200个数据上增加20℃的阶跃泄漏，对温度3在201到300个数据上增加步长为0.1℃的缓慢的泄漏事故，对温度4在501到700个数据上增加步长为0.1℃的缓慢的泄漏事故。四个温度传感器的模拟数据见图3。

分别采用GLT方法以及OPT方法，方法的判断结果见图4(a)～(d)，其中图4(a)为GLT方法泄漏报警函数输出图，图4(b)为OPT方法泄漏报警函数输出图，图4(c)为GLT方法泄漏传感器的定位结果图，图4(d)为OPT方法泄漏传感器的定位结果图。从图4(a)与图4(b)可看出两个方法的泄漏报警原理并且验证了方法的正确性，从图4(c)与图4(d)中可以清楚的看到两个方法都能很好的对泄漏做出报警并且同时给出最有可能发生泄漏的传感器位置。

Claims (6)

1.基于GLT及OPT方法的供热管网运行状态监测方法，其特征在于，包括以下步骤：

步骤1，将监测终端安装于城市供热管网的薄弱环节，一个监测终端包括四路温度数据、两路电导率数据，同时支持一路压力数据以及一路流量数据；其中四路温度数据以及两路电导率数据在使用GLT方法或者OPT方法时必须使用；对于检测同一数据类型的传感器，使其安装位置为检测同一目标而位于不同位置维度，电导率传感器的安装类似；

步骤2，安装监测终端时使安装距离＜200m或保证GPRS通信信号；在GPRS信号良好时，使用GPRS上报数据；由于安装环境问题，当GPRS信号受影响，采用RF射频通信将数据发送到相邻的终端由GPRS信号良好的终端发送到上位机；

步骤3，监测终端安装完成后，运行上位机，在上位机进行第一轮数据接收完成后，及上位机以及接收到管网中全部终端发送的数据后，在上位机运行地图展示，为每一个终端标定位置，为监测系统的搭建做准备；

步骤4，监测系统投入运行后，每次接收到数据后，都由上位机通过GLT方法或者OPT方法将温度数据以及电导率数据作为输入，GLT方法的泄漏检测函数为：

式中，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵，σ²为噪声方差，p为输入Z的线性变换，p'为p的转置，

p＝V'Z (2)

OPT的泄漏检测函数为：

ρ_k ^*为k位置传感器检测值相应的最优奇偶残差，|ρ_k ^*|为其绝对值，T为通过正态分布以及给定的虚警率确定的阈值，如果输出大于阈值，则发出泄漏报警；除此之外，两个方法还能对每一个传感器的泄漏进行判断，得出最有可能发生泄漏的传感器位置；

步骤5，泄漏报警检测完成，将数据加入数据库，单次接收结束。

2.根据权利1要求的基于GLT方法以及OPT方法的供热管网运行状态监测系统，其特征在于，步骤1所述的GLT方法，将m个余度温度传感器或电导率传感器的量测方程为：

Z＝HX+f+ξ (4)

式中,X∈Rⁿ是n个待量测的状态变量，这些待测变量可以是土壤的温度及电导率等。Z∈R^m是m个传感器的量测值，方法要求m≥n；H是传感器配置的测量矩阵，也可称为安装矩阵，主要用来描述从待测量到传感器数据的过程；f为泄漏向量，ξ为噪声向量；若待测量为温度一个量，即n＝1；有四路温度检测，即m＝4；显然温度传感器直接测量温度的大小，即：

H＝[1 1 1 1]' (5)

f是泄漏向量，描述泄漏发生的位置，ξ用来描述传感器在测量时产生的误差，使用m维零均值、协方差为σ²的高斯白噪声序列描述，即E(ξ)＝0,E(ξξ^T)＝σ²I_m，I_m为单位矩阵，由方差定义：

式中，X表示变量，此处为温度传感器的输出温度，μ为总体均值，这里使用温度实际值，若传感器测量误差为1℃，则易得(X-μ)²＝1，因此可得σ²＝1。

3.根据权利1要求的基于GLT方法以及OPT方法的供热管网运行状态监测系统，其特征在于，步骤1所述的GLT方法，是通过提出有泄漏状态和无泄漏状态两种假设，采用基于奇偶空间的广义似然比检验方法来进行泄漏的报警；对于泄漏检测，使用对泄漏解耦的奇偶残差实现泄漏检测函数，考虑量测值Z的线性变换：

p＝V'Z＝V'HX+V'f+V'ξ (7)

式中，f为泄漏向量，ξ为噪声向量，H为安装矩阵，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵；

让p对泄漏位置解耦，需要V满足能够分离出传感器测量的目标值与泄漏位置，即如下约束条件：

V'H＝0 (8)

带入式(8)，此时待观测的温度由于乘了一个0阵即式(7)化简为:

p＝V'f+V'ξ (9)

从式(9)可知，p独立于待测状态X而仅与噪声或可能的泄漏故障有关，实际上，当无泄漏时，p仅是噪声的函数，即:

p＝V'ξ (10)

接下来规定泄漏报警决策的假设检验：

无泄漏状态H₀：E(p)＝0,E(pp')＝σ²V'V；

有泄漏状态H₁：E(p)＝μ≠0,E((p-μ)(p-μ)')＝σ²V'V；其中μ＝V'f，

由于p是ξ的线性函数，因此p服从与ξ一致的m-n维的正态分布，在无泄漏时似然函数为式(11)：

式中，ρ(p|H₀)表示对无泄漏状态H₀构建的似然函数，K为与奇偶向量与噪声方差相关的比例系数，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵，σ²为噪声方差，p为输入Z的线性变换，p'为p的转置；

在无泄漏时似然函数为式(12)

式中，ρ(p|H₁)表示对有泄漏状态H₁构建的似然函数，其中μ＝V'f，其余参数与式(11)一致；

在式(11)与式(12)中，比例系数K为：

联立式(11)与式(12)构造对数似然比函数为：

为了求出最大值可对式(14)关于μ求导得出极大似然函数值为：

所以，可以用p搭建泄漏报警函数：

则通过p服从的正态分布可知FD_GLT～χ²(m-n)，

于是泄漏报警的准则为：FD_GLT≥T时，说明泄漏事故发生，触发报警；否则，判决无泄漏；

阈值T的选择由式(16)推得可知泄漏报警函数符合自由度为m-n的卡方分布即服从中心χ²(m-n)分布，该分布的概率密度函数为

假设给定一个阈值T，则可以通过式(17)计算出虚警率P_F:

或

显然若要求给定的报警准确率为1-α，显然α为虚警率，则只需查χ²分布表，得1-α分位数就可求出需要使用的阈值

4.根据权利1要求的基于GLT方法以及OPT方法的供热管网运行状态监测系统，其特征在于，步骤1所述的GLT方法，可定位出最有可能发生泄漏的位置，即对在膨胀节一圈安装的温度传感器，哪一个传感器所在的位置最有可能发生了泄漏；对于泄漏定位，提出泄漏定位函数，该函数是假设各个传感器位置发生泄漏时，使对应的奇偶残差的似然函数取得最大时的统计量；

泄漏位置定位包括m＝4个假设检验，假设H_i为第i个温度传感器处发生了泄漏，i＝1、2、...、m，：E(p)＝μ_i≠0，E((p-μ_i)'(p-μ_i))＝σ²其中μ_i＝V'e_if_i，e_i为第i行值，i＝1、2、...、m，为1，其余值为0的m维列向量，这里为4维列向量；记V'e_i＝V^*为奇偶矩阵V的第i个行向量的转置向量。

则相应于第i个温度传感器检测数据的对数似然函数为：

式中，表示相应于第i个温度传感器检测数据的对数似然函数，K为与奇偶向量与噪声方差相关的比例系数，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵，V_i ^*为奇偶矩阵V的第i行构成向量的转置向量，σ²为噪声方差，p为输入Z的线性变换，p'为p的转置，f_i为第i个温度传感器的泄漏向量，且μ_i＝V'e_if_i,对式(18)关于f_i求导求出对于f_i的极大似然估计:

极大似然估计值为：

式中，表示相应于第i个温度传感器检测数据的对数似然函数的极大似然估计值，K为与奇偶向量与噪声方差相关的比例系数，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵，V_i ^*为奇偶矩阵V的第i行构成向量的转置向量，σ²为噪声方差，p为输入Z的线性变换，p'为p的转置；如果设计时令V'V＝I_m-n，I_m-n为m-n维单位矩阵，由上式(20)可确定针对各个位置上的温度传感器的泄漏判断函数为：

式中，FI_GLT(i)表示第i个温度传感器的泄漏判断函数值，K为与奇偶向量与噪声方差相关的比例系数，V表示根据系统安装设计的奇偶矩阵，V'为其转置矩阵，V_i ^*为奇偶矩阵V的第i行构成向量的转置向量，σ²为噪声方差，p为输入Z的线性变换，p'为p的转置，在判断泄漏时FI_GLT(i)越大，说明在该位置发生泄漏的可能性越大，对应于所得的极大似然估计值就越大；

如果第k个温度传感器检测的温度值经过计算的泄漏位置定位函数值FI_GLT(k)是所有检测到的温度计算出的函数值中的最值，m为传感器个数，则:

则说明第k个温度传感器所在的位置最有可能已发生了泄漏。

5.根据权利1要求的基于GLT方法以及OPT方法的供热管网运行状态监测系统，其特征在于，步骤1所述的OPT方法将泄漏检测与定位同时进行，

对于式(1)，增加D表示泄漏输入矩阵，F噪声输入矩阵则

Z＝HX+Df+Fξ (23)

对第i个温度传感器所在的位置发生泄漏事故，选取最优奇偶向量v_i，带入式(23)得m＝4个奇偶方程：

r_i＝v_i'Z＝v_i'Df+v_i'Fξ,i＝1,2,…,m (24)

式中，r_i为与第i个温度传感器所在的位置发生泄漏事故这一事件对应的残差；系统发生泄漏时，可以通过下式判断最有可能发生泄漏的位置，式(25)的情况说明最有可能发生泄漏事故的位置在温度传感器i的位置上；

|r_i|＞|r_j|,j＝1,2,…,m,j≠i (25)

为了设计r，建立性能指标函数：

其中，v_i′De_i表示对第i个温度传感器故障的敏感度。v_i′De_j表示对第j个温度传感器故障的敏感度，e_i为m阶单位矩阵I_m的第i个列向量，v_i为设计的向量，||v_i′F||表示对此处使用的温度传感器测量误差的敏感程度；

使S_i最大，是为了让函数对在温度传感器i位置上发生的泄漏事故敏感，对其它位置上的泄漏事故以及测量时产生的误差不敏感，即在第i个温度传感器位置上发生泄漏事故时使|r_i|尽量大，

考虑到：

带入式(26)的的极值问题整理成求：

其中，u＝V′De_i，M_B＝V′(FF′+DD′-De_ie_i′D′)V，为了求得式(26)的极值解，由柯西—士瓦兹不等式定理及其相关推论得：

式(26)在出取得最大值，于是：

其中，u＝V′De_i，M_B＝V′(FF′+DD′-De_ie_i′D′)V，

最后，将v_i标准化，相应的最优奇偶残差为：

标准化以后的奇偶残差ρ_i～N(0,1)，则

最后说明T的选择，对给定的要求的准确率1-α，显然α为虚警率，通过正态分布求得阈值T＝u_1-α/2。

6.根据权利1要求的基于GLT方法以及OPT方法的供热管网运行状态监测系统，其特征在于，步骤1所述的GLT方法以及OPT方法中最优奇偶阵需要满足以下条件：

式中，V阵可以采用Potter方法求取，具体如下：

Step1，通过给定的H阵，通过w＝I-Hinv(H'H)H'计算得到w阵；式中H表示安装矩阵，H'为其转置，I为m维单位矩阵

Step2，根据求得V阵首元素，

v_1j＝w_1j/v₁₁,j＝2,...,m；

v_ij＝0,i＝2,3,...,m-n,j＝1,2,...,i；i＞j；更新第一行除首元素外全部元素以及剩余下三角全部元素为0；式中，w_ij表示w阵第i行第j列的元素，v_ij表示V阵第i行第j列的元素；

Step3，从i＝2开始更新对角元素，w_ii表示w阵第i行第i列的元素，v_ii表示V阵第i行第i列的元素，

然后执行Step4；

Step4，w_ij表示w阵第i行第j列的元素，v_ij表示V阵第i行第j列的元素，利用由Step3使用的i更新该行没有计算的其余元素；

Step5，i增加到下一行，调到Step3继续执行，指导更新得到完整的V阵。