CN107590516B - 基于光纤传感数据挖掘的输气管道泄漏检测识别方法 - Google Patents

Abstract

本发明公开了一种基光纤传感数据挖掘的输气管道泄漏检测识别方法，利用分布式光纤声波/振动传感器拾取沿管道传播的泄漏声波/振动信号，对各空间点采集的泄漏声波/振动信号进行Mel倒谱和AR模型特征提取与特征选择，通过改进的特征规则挖掘方法和正反样本挖掘，建立所选择特征属性与泄漏事件类型的关联规则，以此对输气管道泄漏事件进行实时在线智能检测、识别与分类，解决复杂噪声环境下输气管道在线泄漏检测问题。本发明不仅能检测识别单纯的管道泄漏信号，而且能够将混有不同干扰源的泄漏信号检测并识别出来。

Description

基于光纤传感数据挖掘的输气管道泄漏检测识别方法

技术领域

基于光纤传感数据挖掘的输气管道泄漏检测识别方法，用于输气管道泄漏检测，涉及管道安全、光纤声波/振动信号处理及数据挖掘技术领域。

背景技术

天然气输送管道在生产制作过程中存在一些小缺陷，在长时间的输送荷载以及输送介质的高温、腐蚀环境作用下，管道腐蚀会加快、加重，最终导致管道破损泄漏事故发生，造成巨大的直接经济损失，同时诱发一系列次生灾害，如火灾和环境污染等，对人民生命安全构成威胁。因此，天然气管道输送过程中的安全问题，特别是天然气管道在线实时泄漏检测，一直是安全生产需要考虑的重中之重。

目前输气管道的泄漏监测方法有许多种，根据所监测的管内气体流动参数不同，可分为质量/体积平衡法、应用统计法、负压波法、瞬态模型法、分布式光纤法等。其中基于分布式光纤传感的输气管道泄漏检测方法因其具有灵敏度好、定位精度高、检测时间短、适应能力强、系统寿命长等优点而受到广泛关注和应用。当输气管道的某一点破裂产生泄漏时，将打破输气管道中正常压力平衡，管道系统内流体弹性能释放，引起瞬间振荡而产生声波。该泄漏声波由气体泄漏的激动能量产生，为持续发射的信号，频率可达几百kHz，主要沿管壁进行传播。利用沿管道安装的光纤可以监听并采集这种声波信号，通过对声波信号进行特征提取及识别，判断管道是否发生泄漏。

但是基于分布式光纤声波/振动传感的输气管道泄漏检测方法，通常由于泄漏信号弱、淹没在复杂的环境噪声中，而且不同管段埋设环境不同，环境噪声复杂多变，造成实际泄漏检测困难。此外，管道参数如管材、管径和管内压力不同，对泄漏声/振动信号频谱等特征影响较大。因此，在实际应用中，输气管道泄漏的检测与识别仍是一个复杂的应用难题。

发明内容

本发明的目的在于：解决现有分布式光纤声波/振动传感的输气管道泄漏检测方法，因泄漏信号弱，淹没在了复杂多变的环境噪声中，造成泄漏检测困难的问题；提供了一种基于分布式光纤声波/振动传感及数据挖掘的输气管道泄漏检测识别方法。

本发明采用的技术方案如下：

一种基于光纤传感数据挖掘的输气管道泄漏检测识别方法，其特征在于，如下步骤：

步骤1、沿管道铺设的探测光缆，感知和采集各空间点上的声波/振动信号；

步骤2、基于构建好的典型事件数据库构建关联规则分类器；

步骤3、对构建好的关联规则分类器对声波/振动信号进行在线识别和分类，若出现新的未知事件，更新典型事件数据库。

进一步，所述步骤2的具体步骤如下：

步骤2.1、基于输气管道实际发生的典型事件类型，分别构建不同类型事件的样本数据库；

步骤2.2、对样本数据库中的每条数据记录，进行信号特征提取和特征选择；

步骤2.3、对提取和选择后得到的优化特征矩阵及不同类型事件标签进行二值化处理，并将二值化处理后的结果进行拼接，得到二值化训练集；

步骤2.4、基于二值化特征训练集进行关联分析并构建关联规则分类器。

进一步，所述步骤2.1的具体步骤如下：

步骤2.11、对某空间点内采集的时间序列进行分帧处理，将每一帧时间信号作为一条数据记录放在样本数据库中；

步骤2.12、根据实际发生的典型事件类型对每条数据记录分别进行标注，构建不同类型事件的样本数据库。

进一步，所述步骤2.2的具体步骤如下：

步骤2.21、对每条数据记录提取N₁阶Mel频率倒谱系数；

步骤2.22、对每条数据记录提取N₂阶AR模型系数；

步骤2.23、将步骤2.21和步骤2.22提取的特征值合并进行特征选择。

进一步，所述步骤2.23的具体步骤如下：

将步骤2.21提取的C＝{C(i),i＝1,2,…N₁}和步骤2.22提取的A＝{A(j),j＝ 1,2,…,N₂}特征值合并，构成每条数据记录的N^*维特征向量{F(l),l＝1,2,…,N^*}＝ {C(i),A(j),(i＝1,2,…,N₁,j＝1,2,…,N₂)}，对样本数据库中的Cnum条数据记录，都分别进行以上特征值提取的操作，得到典型泄漏事件训练数据的特征矩阵F＝{F(c)(l),(c＝ 1,2,…,Cnum,l＝1,2,…,N^*)}，特征矩阵F中的行代表数据记录序数，列代表特征维数；对特征矩阵F＝{F(c)(l),(c＝1,2,…,Cnum,l＝1,2,…,N^*)}中的每一个特征列F_l＝{F_l(c),C＝ 1,2,…,Cnum}，分别按照：(F_l(c)-min(F_l))/(max(F_l)-min(F_l))进行最大/最小规范化，得到N^*维取值在0-1范围内的特征列F_l′＝{F_l′(c),C＝1,2,…,Cnum}，然后求其方差

其中Fsev是该列特征的平均值

对每一个特征列的方差值FT_l(l＝1,2,…,N^*)从大到小排序，选取其中方差较大的前N个特征，0<N≤N^*，所选特征在原特征矩阵F中的对应列的标号{n₁,n₂,…n_N}∈{1,2,…,N^*}，删除原特征矩阵F中未选择的特征列，得到一个Cnum×N维的优化特征矩阵 F_s＝{F_s(c)(n),(c＝1,2,…,Cnum,n＝1,2,…,N)}。

进一步，所述步骤2.3的具体步骤如下：

步骤2.31、基于模糊C均值聚类算法对选取的特征矩阵进行二值化处理，具体如下：

通过模糊C均值聚类算法将提取和选择后的特征矩阵 F_s＝{F_s(c)(n),(c＝1,2,…,Cnum,n＝1,2,…,N)}转换成二值化特征矩阵 FA＝{FA(c)(n),(c＝1,2,…,Cnum,n＝1,2,…,N)}，具体操作如下：将某一列的特征值 F_s(:,n)＝{F_s(c)(n),(c＝1,2,…,Cnum)}通过模糊C均值聚类算法分为两堆，每一堆各有一个中心点，较大的中心点的值为big(n)，较小中心点的值为sma(n)，与当前特征值F_s(c)(n)比较大小，若|F_s(c)(n)-big(n)|<|F_s(c)(n)-sma(n)|，则令FA(c)(n)＝1，反之则令 FA(c)(n)＝0；另设一个与FA取值相反的二值化特征矩阵

后续关联规则挖掘时将综合两组二值化特征矩阵的分析结果；将二值化特征矩阵的每一列，作为后续数据挖掘事务的一项，记为：feature(n),n＝1,2,…,N；另外，保存特征矩阵F_s的N对特征聚类中心点{Cen(n)＝{big(n),sma(n)},n＝1,2,…,N}，作为后续在线检测和识别时的特征二值化参考依据；

步骤2.32、对不同类型事件标签进行二值化处理，具体如下：

将V类事件分别表示为一个V位的二值化数值，若当前数据记录c为第β类，则第β位为1，其余都为0；以此类推，Cnum条样本记录的事件类型标注将转化为一个Cnum×V的二值矩阵，最终得到二值化事件标签矩阵L_Cnum×V；

将二值化事件标签矩阵的每一列，也作为后续数据挖掘的一项，记为：label(β),β＝ 1,2,…,V；

步骤2.33、拼接特征矩阵与事件标签矩阵，具体如下：

将记为feature(n),n＝1,2,…,N的特征列与记为label(β),β＝1,2,…,V的事件标签列拼接，对已有的两个取值相反的二值化特征矩阵FA_Cnum×N和FB_Cnum×N，分别与二值化事件标签矩阵L_Cnum×V对应拼接，得到两个包含特征列和事件标签的二值化训练集：FA′_Cnum×(N+V),FB′_Cnum×(N+V)。

进一步，所述步骤2.4的具体步骤如下：输入二值化特征训练集： FA′_Cnum×(N+V),FB′_Cnum×(N+V)，通过Apriori算法挖掘提取feature→label的关联规则，为每一条规则计算以下四个规则评价指标：支持度s(feature→label)＝σ(feature∪label)/σ(样本)，σ为样本计数值、置信度c(feature→label)＝σ(feature∪label)/ σ(feature)、Kulc度量k(feature→label)＝(c(feature→label)+(label→feature))/ 2和不平衡度IR(feature→label)＝|σ(feature)-σ(label)|/|σ(feature)+σ(label)- σ(feature∪label)|，设置相应的四个阈值：最小支持度mins、最小置信度minc、最小Kulc 度量mink、最大不平衡度maxIR，其中，mins的值依据实际样本数量和样本组成设定；minc 是评价规则的重要程度；mink的值与规则的正反两个置信度有关；不平衡度IR通常与Kulc 度量一起使用，按照设定的阈值修剪不满足阈值要求的关联规则；

对正反两组二值化训练集输入，将得到两组有关特征参数与事件类别feature与label的重要关联规则作为输出，提取这些规则的前件，即所提取的关联规则中feature的集合，用 RA(x)表示从二值化特征训练集A′_Cnum×(N+V)挖掘的每条规则前件，用RB(y)表示从二值化特征训练集B′_Cnum×(N+V)挖掘的每条规则前件，则按不同的事件类别β分别形成规则集，得到各类事件的规则集R＝{RA(β)(x),RB(β)(y),(β＝1,2,…,V；x＝1,2,…,X_β；y＝1,2,…,Y_β)}， X_β，Y_β分别为事件β所对应的两组规则前件个数，保存此规则集R完成关联规则分类器的构建，将其作为测试信号分类的参照特征和判决条件。

进一步，所述步骤3的具体步骤如下：

对采集到的某空间点时序测试声波/振动信号按q毫秒分帧，每一帧信号XX^*＝{XX^*(m),m＝1,2,…,M}即为一条测试数据记录，包含

个采样点，其中，有U 点的重叠区域，M取整数；

对当前数据记录提取N₁阶Mel频率倒谱系数{C(i),i＝1,2,…,N₁}，和N₂阶AR模型系数 {A(j),j＝1,2,…,N₂}，组成N^*维特征向量F^*＝{F^*(l),l＝1,2,…,N^*}＝{C(i),A(j),(i＝ 1,2,…,N₁,j＝1,2,…,N₂)}，特征选择时，按照保存的N个特征标号{n₁,n₂,…n_N}，从特征向量F^*中提取这N维系数，得到N维特征向量F_s ^*＝{F_s ^*(n),n＝1,2,…,N}；

采用二值化方法和保存的N对聚类中心点{Cen(n)＝{big(n),sma(n)},n＝1,2,…,N}，对特征向量F_s ^*＝{F_s ^*(n),n＝1,2,…,N}的每一个值计算min{|F_s ^*(n)-big(n)|,|F_s ^*(n)- sma(n)|}，如果距较大中心点更近，则令FA^*(n)＝1，反之则令FA^*(n)＝0，得到测试数据记录的特征向量F_s ^*的二值化特征向量FA^*＝{FA^*(n),n＝1,2,…,N}；

通过比较二值化特征向量FA^*与关联规则分类器中的规则集R满足的关系来进行训练事件类型的识别与分类，具体如下：设FA^*中所有为1的序号集合为D_A，所有为0的序号集合为D_B，则FA^*中为1和0的序号集D＝{D_A＝{ε},D_B＝{γ},(FA^*(ε)＝1,FA^*(γ)＝0)}，对比分类器关联规则分类器中的规则集 R＝{RA(β)(x),RB(β)(y),(β＝1,2,…,V；x＝1,2,…,X_β；y＝1,2,…,Y_β)}，设RA_β＝ {RA(β)(x),x＝1,2,…,X_β}，RB_β＝{RB(β)(y),y＝1,2,…,Y_β}表示针对事件类型β的两组规则集，分别计算序号集D满足集合R中各类事件关联规则集的数量占比Count(β)：

依据上式统计结果，比较各类事件关联规则集的满足数量占比，计算最大值 maxC＝max{Count(β),β＝1,2,…,V}，当最大值满足阈值，即maxC>thr时，thr为事件判定需满足规则占比的最小阈值，该测试信号所属的事件类别为Count(β)为最大时对应的那个事件类别标号，记为β_k；当最大值不满足阈值，则令β_k＝0，得到β_k的值后，即可进行事件类型判定；

若出现新的未知事件，作为待定事件类型β^*，继续补充事件样本数据库，作为后期更新关联规则分类器的依据。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

1、本发明利用分布式光纤声波/振动传感器拾取沿管道传播的泄漏声波/振动信号，对各空间点采集的泄漏声波/振动信号进行Mel倒谱和AR模型特征提取与特征选择，通过关联分析方法挖掘和建立所选择特征属性与泄漏事件类型的关联规则，以此基于分布式光纤声波/振动传感及数据挖掘对输气管道泄漏事件进行实时在线智能检测、识别与分类。

2、本发明不仅能检测识别单纯的管道泄漏信号，而且能够将混有不同干扰源的泄漏信号检测并识别出来，抗干扰能力比较强，特别适合复杂噪声环境中的泄漏声波检测识别。

3、本发明在进行关联规则挖掘之前，提出一种特征选择方法，用于提取和选择可分辨性较大的特征项，去除冗余特征项，大大提升关联规则挖掘的效率。

4、本发明中采用一种改进的Apriori算法，基于目标事件类别只挖掘包含某类事件项的频繁项集，得到特征项集与事件项的关联规则，以此有效提升关联规则的挖掘效率；并且基于Kulc度量和IR不平衡度修剪关联规则，可以有效提升基于关联规则构建的分类器的分类效果和效率。

5、本发明中采用一种改进的特征规则挖掘方法，对信号的特征矩阵二值化后，获得正反两组特征布尔矩阵，利用正反样本并行去挖掘特征与事件的关联规则，使得挖掘规则更加全面，可以有效的提升基于关联规则所构建的分类器的分类效果。

附图说明

图1是本发明采集的输气管道泄漏声信号时空矩阵；

图2是本发明的输气管道泄漏检测识别方法流程图；

图3是本发明特征选择后得到的待识别5类事件的9维特征图；

图4是本发明中对不同类型事件标签进行二值化处理的示意图；

图5是本发明管道泄漏声波幅值衰减随距离的变化规律示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

实施例一

基于分布式光纤声波/振动传感(DAS)的输气管道在线监测系统，系统硬件主要由三个部分组成，探测光缆、光信号解调设备、信号处理主机。探测光缆通常采用普通单模通信光纤或声波、振动增敏型传感光缆，一般沿输气管道内壁或外壁铺设。光信号解调设备是该系统的核心，其内部组成器件主要包括光学和电学器件两类。在解调设备中超窄线宽激光器发出连续光信号经声光/电光调制器调制成光脉冲，依次经掺铒光纤放大器(EDFA)、分布式拉曼放大器、隔离器、环形器注入到探测光缆，沿光缆传输过程中产生的后向瑞利散射光信号沿光缆返回，并由环形器接收，滤波后经马赫曾德尔干涉仪(MZI)或其他干涉仪等，解调输出外界声波或扰动带来的相位变化信息，即可获得管道泄漏或其他干扰事件产生的声波或振动信号，解调出的信号经网口等形式接口实时传输给信号处理主机。信号处理主机为普通电脑主机或FPGA/DSP等嵌入式主机，用于对获得的声波或振动传感信号进行实时分析、处理，对管道泄漏或其他干扰等异常事件进行检测、识别和分类，并对真实泄漏点进行预警、定位。

基于分布式光纤声波/振动传感系统采集的输气管道泄漏声信号时空矩阵如图1所示，图中横坐标为空间轴，表示管道沿线的数据采集空间点，两空间点的距离为ΔL；纵坐标为时间轴，两采样点采样间隔ΔT＝1/f_s，f_s为时间轴上脉冲触发频率，即时间采样频率。在进行泄漏在线监测时，首先针对各空间点的时间序列进行光纤声波/振动信号的识别和分类，然后基于识别出的泄漏信号空间点的分布情况及声波传播规律进行管道泄漏点定位。

实施例二

在实施例一的基础上，各空间点泄漏信号的识别和分类是本发明的关键，具体流程框架如图2所示，分为两部分：第一部分，基于输气管道的典型事件数据库进行关联规则挖掘和分类器训练，即构建关联规则分类器；第二部分，利用训练好的关联规则分类器对管道泄漏进行在线识别和分类。

实施例三

在实施例二的基础上，基于输气管道实际发生的典型事件类型，分别构建不同类型事件的样本数据库，具体操作过程如下：对某空间点t秒内采集的时间序列进行分帧处理，将每一帧时间信号作为一条数据记录放在样本数据库中。例如，对k秒内采集的时间序列

t是采样时间，为正整数，f_s为时间采样频率，单位为Hz。每隔q毫秒分为一帧，q一般取20～80，则每一帧包含

(取整)个采样点： XX＝{XX(m),m＝1,2,…,M}。随时间滑动依次选取不同时间段帧长为M的时间序列。为保证数据的连续性，一段设置U个采样点的重叠区域，U一般取M的

或

因此，对t秒时间序列进行分帧后得到的数据帧数为Ψ＝t×1000÷q(取整)，相当于t秒时间序列经过分帧后得到Ψ条数据记录。

根据采集信号的事件类型对每条数据记录分别进行标注：纯管道泄漏事件，标注为1；泄漏加走路声干扰，标注为2；说话声干扰，标注为3；走路声干扰，标注为4；正常无事件，标注为5。将所有类型事件数据样本添加到样本数据库中，完成样本数据库的构建，设样本数据库中不同类型事件的数据样本记录总数为Cnum。

实施例四

在实施例三的基础上，对样本数据库中的每条数据记录，记录序号为c(c＝1,2,…,Cnum)，进行信号特征提取和特征选择，具体方法如下：

1、提取N₁阶Mel频率倒谱系数

对每条数据记录，即每一帧声波/振动信号XX＝{XX(m),m＝1,2,…,M}，M为一帧的数据样本点数，乘上汉明窗：

一般情况下a取0.46，得到：

XX′＝XX×W (1)

再对XX′进行M点快速傅里叶变换并求模的平方，得到声波/振动信号在频谱上的能量分布E＝{E(m),m＝1,2,…M}，频谱数据的频率间隔为

f_s为时间采样频率。

定义一组有K个滤波器的Mel尺度滤波器组

K通常取 22～26，采用的滤波器为三角滤波器，中心频率设为{f(k),k＝1,2,…,K}，各中心频率f(k)之间的间隔随着k的减小而缩小，随着k的增大而增宽，三角形滤波器的频率响应定义如下：

将声波/振动信号在频谱上的能量分布E作为Mel滤波器组的输入，得到输出：

M_1,K＝E_1,M×H_M,K (3)

式(3)中输出的M_1,K＝{M(k),k＝1,2,…,K}即为声波/振动信号的Mel频谱。

定义离散余弦变换(DCT)系数为：

k＝1,2,…,K,i＝ 1,2,…N₁}，对M取对数，并通过离散余弦变换(DCT)求倒谱，得到：

将式(4)中得到的N₁阶Mel频率倒谱系数

作为要提取的第一组信号特征，记为C＝{C(i),i＝1,2,…N₁}，N₁一般为12～16。

2、提取N₂阶AR模型系数

AR模型是现代功率谱估计方法之一，其假设研究过程是由一个输入序列激励一个线性系统G(z)的输出。具体为：对每条数据记录，即每一帧声波/振动信号XX＝{XX(m),m＝1,2,…,M}，计算一个假设参数模型

通过这个假设参数模型来估计线性系统传递函数

中参数数组A，其中z是系统传递函数G的自变量，是一个复数，N₂是模型的阶数。利用Levinson-Durbin递推算法可以求解AR模型中的参数数组A，具体为：设A_j(θ)为AR模型在j阶次时的第θ个系数(θ＝1,2,…,j)，ρ_j为j阶系统白噪声的方差，设声波/振动信号第j阶系统的输入自相关函数为R_XXj＝ E[XX(m)XX(m+j)]，则基于声波/振动信号自相关系数，首先可以得到对一阶AR模型的初始值：

A₁(1)＝-R_XX(1)/R_XX(0) (5)

ρ₁＝R_XX(0){1-A₁ ²(1)} (6)

式(5)中的R_XX(0)，R_XX(1)分别为声波/振动信号时延为0和1时的自相关系数。由此，AR模型的高阶参数则依据Levinson-Durbin递推算法分别推导得到：

式(7)中，j是系统递推中的当前阶次，A_j(j)是AR模型在阶次时的最后一个系数，A_j-1(θ) 是AR模型在j-1阶次时的第θ个系数(θ＝1,2,…,j-1)。由式(7)中的A_j-1(θ),A_j(j)可以推导出模型在j阶次时的前j-1个系数{A_j(θ),θ＝1,2,…,j-1}，以及i阶系统白噪声方差ρ_j：

A_j(θ)＝A_j-1(θ)+A_j(j)A_j-1(j-θ) (8)

ρ_j＝ρ_j-1[1-A_j ²(j)] (9)

对N₂阶的AR模型，当前阶次有j＝1,2,…N₂，依次利用式(7)、(8)、(9)得到不同阶次的参数。在递推过程中，若式(7)中求得|A_j(j)|＝1时，则提前停止递推；否则就继续递推，直至求解到最高阶的最后一个系数。

将求解出的N₂阶AR模型的N₂个参数

作为提取的第二组信号特征，记为数组A＝{A(j),j＝1,2,…,N₂}，N₂一般为3～6。

将步骤1和步骤2提取的特征值合并，得到N^*维(N^*＝N₁+N₂)特征向量，构成每条数据记录的N^*维特征向量{F(l),l＝1,2,…,N^*}＝{C(i),A(j),(i＝1,2,…,N₁,j＝1,2,…,N₂)}。对样本数据库中的Cnum条数据记录，都分别进行以上特征提取的操作，得到典型泄漏事件训练数据的特征矩阵F＝{Fc(l),(c＝1,2,…,Cnum,l＝1,2,…,N^*)}，特征矩阵F中的行代表数据记录序数，列代表特征维数,C代表的是Mel系数，小c是特征矩阵的索引。

3、特征选择

基于步骤1和步骤2中得到的特征矩阵F维数较高，为了避免特征冗余，提炼可分辨度大的特征项，提升特征规则挖掘效率，需要进行特征选择来提取最优特征列，具体步骤如下：对特征矩阵F＝{F(c)(l),(c＝1,2,…,Cnum,l＝1,2,…,N^*)}中的每一个特征列 F_l＝{F_l(c),c＝1,2,…,Cnum}，分别按照：(F_l(c)-min(F_l))/(max(F_l)-min(F_l))进行最大 /最小规范化，得到N^*维取值在0-1范围内的特征列F_l′＝F_l′(c),c＝1,2,…,Cnum}，然后求其方差

其中Fsev是该列特征的平均值

对每一个特征列的方差值FT_l(l＝1,2,…,N^*)从大到小排序，选取其中方差较大的前N个特征，0<N≤N^*，建议为6～10，保存所选特征在原特征矩阵F中的对应列的标号{n₁,n₂,…n_N}∈{1,2,…,N^*}，删除原特征矩阵F中未选择的特征列，得到一个 Cnum×N维的优化特征矩阵F_s＝{F_s(c)(n),(c＝1,2,…,Cnum,n＝1,2,…,N)}。

在本实施例中，提取了12维MFCC，5维AR模型系数，对这些特征量进行特征选择后，最终得到的9维特征对5类事件的区分效果如图3所示。由图3可见，利用本发明方法选择的9维特征对待识别5类事件的可分辨度比较大。

实施例五

在实施例四的基础上，对提取和选择后的优化特征矩阵及不同类型事件标签进行二值化及其他数据挖掘预处理，具体步骤如下：

1、基于模糊C均值聚类算法(FCM)对选取的特征矩阵进行二值化处理

通过模糊C均值聚类算法(FCM)将提取和选择后的特征矩阵F_s＝{F_s(c)(n),(c＝1,2,…,Cnum,n＝1,2,…,N)}转换成二值化特征矩阵FA＝{FA(c)(n),(c＝1,2,…,Cnum,n＝ 1,2,…,N)}，具体操作如下：将某一列的特征值F_s(:,n)＝{F_s(c)(n),(c＝1,2,…,Cnum)}通过模糊C均值聚类算法(FCM)分为两堆，每一堆各有一个中心点，较大的中心点的值为big(n)，较小中心点的值为sma(n)，与当前特征值F_s(c)(n)比较大小，若|F_s(c)(n)-big(n)|< |F_s(c)(n)-sma(n)|，则令FA(c)(n)＝1，反之则令FA(c)(n)＝0。为了避免数据挖掘时忽略某些事件类型中二值化为0的特征值的作用，利用正反样本并行挖掘关联规则，提升分类效果，另设一个与FA取值相反的二值化特征矩阵

后续关联规则挖掘时将综合两组二值化特征矩阵的分析结果。将二值化特征矩阵的每一列，作为后续数据挖掘事务的一项，记为：feature(n),n＝1,2,…,N。另外，保存特征矩阵F_s的N对特征聚类中心点{Cen(n)＝{big(n),sma(n)},n＝1,2,…,N}，作为后续在线检测和识别时的特征二值化参考依据。

2、对不同类型事件标签进行二值化处理

本发明实施例中样本数据库共标注了5类事件，分别为：纯泄漏事件，标注为1；泄漏加走路声干扰，标注为2；说话声干扰，标注为3；走路声干扰，标注为4；正常无事件，标注为5。将这5类事件分别表示为一个5位的二值化数值，若当前数据记录c为第β类，则第β位为1，其余都为0。例如β＝3时，二值化表示如图4所示。

以此类推，Cnum条样本记录的事件类型标注将转化为一个Cnum×5的二值矩阵，最终得到二值化事件标签矩阵L_Cnum×5。

将二值化事件标签矩阵的每一列，也作为后续数据挖掘的一项，记为：label(β),β＝ 1,2,…,5。

3、拼接特征矩阵与事件标签矩阵

将记为feature(n),n＝1,2,…,N的特征列与记为label(β),β＝1,2,…,5的事件标签列拼接，对已有的两个取值相反的二值化特征矩阵FA_Cnum×N和FB_Cnum×N，分别与二值化事件标签矩阵L_Cnum×5对应拼接，得到两个包含特征列和事件标签得二值化训练集：FA′_Cnum×(N+5),FB′_Cnum×(N+5)。

实施例六

在实施例五的基础上，基于二值化特征训练集进行关联分析并构建关联规则分类器，具体方法如下：输入二值化特征训练集：FA′_Cnum×(N+5),FB′_Cnum×(N+5)，基于改进的Apriori算法，只挖掘包含某类事件项的频繁项集，提取feature→label的关联规则。为每一条规则计算以下四个规则评价指标：支持度s(feature→label)＝σ(feature∪label)/ σ(样本)，σ为样本计数值、置信度c(feature→label)＝σ(feature∪label)/ σ(feature)、Kulc度量k(feature→label)＝(c(feature→label)+c(label→ feature))/2和不平衡度IR(feature→label)＝|σ(feature)-σ(label)|/|σ(feature)+ σ(label)-σ(feature∪label)|，设置相应的四个阈值：最小支持度mins、最小置信度minc、最小Kulc度量mink、最大不平衡度maxIR。一般来说，mins和minc是关联分析中最常用的基本评价指标，mins的值依据实际样本数量和样本组成设定；minc是评价规则的重要程度，推荐设置为0.5以上；另外引入了两个关联规则修剪指标，不平衡度IR通常与Kulc度量，最小Kulc度量mink的值与规则的正反两个置信度有关，推荐设置为略大于minc的值；不平衡度IR通常与Kulc度量一起使用，一般认为最大不平衡度maxIR为0.5以下的规则是比较平衡的。按照设定的阈值修剪不满足阈值要求的关联规则，修剪后得到规则的是排除了冗余信息、可用于分类的重要关联规则。

对正反两组二值化训练集作为输入进行挖掘，将得到两组有关特征参数与事件类别 (feature与label)的重要关联规则作为输出。提取这些规则的前件，即所提取的关联规则中feature的集合，用RA(x)表示从二值化特征训练集FA′_Cnum×(N+5)挖掘的每条规则前件，用RB(y)表示从二值化特征训练集FB′_Cnum×(N+5)挖掘的每条规则前件，则按不同的事件类别β分别形成规则集合，得到各类事件的规则集R＝{RA(β)(x),RB(β)(y),(β＝1,2,…,5；x＝1,2,…,X_β；y＝1,2,…,Y_β)}，X_β，Y_β分别为事件β所对应的两组规则前件个数。保存此规则集R完成关联规则分类器的构建，将其作为测试信号分类的参照特征和判决条件。

本实施例中，针对样本数据库中的5类事件，得到了5类事件的规则集R，例如规则集中第2类事件的正反两组关联规则集合可以表示为：

RA_β＝2＝{{4},{3,4},{4,6},{4,7},{4,9},{4,6,7},{4,6,9},{4,7,9}}

RB_β＝2＝{{2},{2,5},{1,2},{1,2,5},{8},{5,8}}

关联规则集中的每一个数字集合表示当前事件挖掘的每一条关联规则，每一条关联规则包含的是特征维数的集合。对关联规则按照阈值修剪前后，规则数量对比如表1所示，大大提高了后面规则挖掘的效率。

事件类别	修剪前规则数量	修剪后规则数量
			1：纯管道泄漏事件	78	47
2：泄漏加走路声干扰	82	13
			3：说话声干扰	106	25
4：走路声干扰	90	0
			5：正常无事件	70	4

表中，修剪前规则是指挖掘过程中，只设置最小支持度mins和最小置信度minc阈值后得到的关联规则；修剪后规则是指再设置最小Kulc度量mink和最大不平衡度maxIR阈值后得到的关联规则。

实施例七

在实施例六的基础上，基于关联规则分类器的泄漏信号在线检测与识别过程，具体方法如下：

参照实施例三中的样本信号时间分帧方法，对采集到的某空间点时序测试声波/振动信号按q毫秒分帧，每一帧信号XX^*＝{XX^*(m),m＝1,2,…,M}即为一条测试数据记录，包含

个采样点，有U点的重叠区域，M取整数。

参照实施例四中的声波/振动信号特征提取和选择方法，对当前数据记录提取N₁阶Mel 频率倒谱系数{C(i),i＝1,2,…,N₁}，和N₂阶AR模型系数{A(j),j＝1,2,…,N₂}，组成N^*维特征向量F^*＝{F^*(l),l＝1,2,…,N^*}＝{C(i),A(j),(i＝1,2,…,N₁,j＝1,2,…,N₂)}。特征选择时，按照保存的N个特征标号{n₁,n₂,…n_N}，从特征向量F^*中提取这N维系数，得到N维特征向量F_s ^*＝{F_s ^*(n),n＝1,2,…,N}。

参照实施例五中的二值化处理方法，以及保存的N对聚类中心点 {Cen(n)＝{big(n),sma(n)},n＝1,2,…,N}，对特征向量F_s ^*＝{F_s ^*(n),n＝1,2,…,N}的每一个值计算min{|F_s ^*(n)-big(n)|,|F_s ^*(n)-sma(n)|}，如果距较大中心点更近，则令 FA^*(n)＝1，反之则令FA^*(n)＝0，得到测试数据记录的特征向量F_s ^*的二值化特征向量FA^*＝{FA^*(n),n＝1,2,…,N}。

通过比较二值化特征向量FA^*与关联规则分类器中的规则集R满足的关系来进行训练事件类型的识别与分类，具体如下：设FA^*中所有为1的序号集合为D_A，所有为0的序号集合为D_B，则FA^*中为1和0的序号集D＝{D_A＝{ε},D_B＝{γ},(FA^*(ε)＝1,FA^*(γ)＝0)}，对比分类器关联规则分类器中的规则集 R＝{RA(β)(x),RB(β)(y),(β＝1,2,…,5；x＝1,2,…,X_β；y＝1,2,…,Y_β)}，设RA_β＝ {RA(β)(x),x＝1,2,…,X_β}，RB_β＝{RB(β)(y),y＝1,2,…,Y_β}表示针对事件类型β的两组规则集，分别计算序号集D满足规则集R中各类事件关联规则集的数量占比Count(β)：

依据上式统计结果，比较各类事件关联规则集的满足数量占比，计算最大值 maxC＝max{Count(β),β＝1,2,…,5}，当最大值满足阈值，即maxC>thr时，thr为事件判定需满足规则占比的最小阈值，推荐设置为0.1～0.3左右，该测试声波/振动信号所属的事件类别为Count(β)为最大时对应的那个事件类别标号，记为β_k；当最大值不满足阈值，则令β_k＝0。得到β_k的值后，即可进行事件类型判定：如果β_k＝1，则判定为纯泄漏事件，如果β_k＝2，则判定为泄漏加走路声干扰，如果β_k＝3/4/5，则判定为此刻无泄漏事件发生，管道运行一切正常，如果β_k＝0，则判定为未知事件，即当前分类器无法识别的事件。

基于以上方法进行试验，试验中对每一类事件取3000条样本记录，5类事件共15000 条样本记录，采用十折交叉验证的方法，即将各类事件样本均分成10份，依次取其中1份作为测试集，剩余9份作为训练集，一共试验10次，得到此输气管道泄漏识别方法的平均识别效果如表2所示：

表2中精确率(precision)表示的是判定为某类的样本中有多少是真正为该类样本，召回率(recall)表示的是某类样本中有多少被正确判定了。例如对某类事件β，设β类的样本中被正确判定为β的数量为tp，β类的样本中被错误判定为其他事件的数量为fp,不是β类的其他事件样本中被判定为β的数量为fn，不是β类的其他事件样本中被判定为除β类事件以外的其他事件的数量为tn，则

综合评价指标(f-measure) 是精确率(precision)和召回率(recall)的加权调和平均：

当参数a＝1时，就是最常见的f1-measure了。精确率(precision)和召回率(recall)指标有的时候是矛盾的，综合考虑这两个度量，可以更有力的说明分类效果。由表2的识别结果，可以看出该发明方法对于识别噪声环境下的泄漏信号的有效性。

若出现新的未知事件，作为待定事件类型β^*，继续补充事件样本数据库，作为后期更新关联规则分类器的依据。更新关联规则分类器时，利用实施例四、五、六中的信号处理和关联分析方法挖掘事件关联规则，基于挖掘出的关联规则测试β^*类事件的识别率，若识别率满足要求则可以在关联规则分类器中添加新类别β^*对应的规则集合

后期基于更新后的分类规则集R^*实现新添加事件的检测与识别。

实施例八

在施例七的基础上，对所有空间点泄漏信号分别进行识别分类后，基于泄漏信号空间点的分布情况及声波传播规律进行管道泄漏点定位，具体方法为：如图5所示，由声波传播规律可知，输气管道发生泄漏时，泄漏声波信号从泄漏点沿管径向管道两端传播，且传播规律呈指数衰减：

P＝P_A×e^-ux (11)

式(11)中P_A为泄漏点信号幅值，P为到泄漏点距离为x的空间点的信号幅值，衰减系数 u＝σα，其中α为阻尼吸收系数，σ为修正因子，对输气管道一般为0.5-1.5。实际应用中，管道上下游阻尼系数不同，根据公式(11)可知通过两个空间点的信号幅值即可求解一组衰减系数，通过上下游衰减系数和两个空间点幅值可以计算出当前空间点距泄漏点的距离，即利用公式(11)通过四个空间点可以实现泄漏点定位。

设管道沿线此时共有W个已被关联规则分类器判定为泄漏的连续空间点，各空间点间距为ΔL，在该段包含W个泄漏点的泄漏区间内(W>4)，找出原始信号峰值处的那个空间点，设其幅值为P₀，在该点上游取最近的两个空间点，设其幅值分别为P_-2,P_-1，下游取最近的两个空间点，设其幅值分别为P₁,P₂；设上游的衰减系数为u^-，下游的衰减系数为u⁺；设泄漏点产生的声波幅值为P_A，P_A到P₁的距离为X，则P_A到P_-1的距离为2ΔL-X，P_A到P₂的距离为ΔL+X，P_A到P_-2的距离为3ΔL-X。将这4个空间点的相关参数代入公式(11)可以得到一组泄漏声波传播方程：

由(12)/(13)、(14)/(15)可分别解出上下游衰减系数：

再由(12)/(14)可得：

ln(P₁/P_-1)＝u^-(2ΔL-X)-u⁺X (18)

将(16)、(17)代入(18)，可解得

Location就是泄漏点到空间点P₁的距离，即为当前方法确定的泄漏点位置。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种基于光纤传感数据挖掘的输气管道泄漏检测识别方法，其特征在于，如下步骤：

步骤1、沿管道铺设的探测光缆，感知和采集各空间点上的声波/振动信号；

步骤2、基于构建好的典型事件数据库构建关联规则分类器；

步骤3、对构建好的关联规则分类器对声波/振动信号进行在线识别和分类，若出现新的未知事件，更新典型事件数据库；

所述步骤2的具体步骤如下：

步骤2.1、基于输气管道实际发生的典型事件类型，分别构建不同类型事件的样本数据库；

步骤2.2、对样本数据库中的每条数据记录，进行信号特征提取和特征选择；

步骤2.3、对提取和选择后得到的优化特征矩阵及不同类型事件标签进行二值化处理，并将二值化处理后的结果进行拼接，得到二值化训练集；

步骤2.4、基于二值化特征训练集进行关联分析并构建关联规则分类器；

所述步骤2.1的具体步骤如下：

步骤2.11、对某空间点内采集的时间序列进行分帧处理，将每一帧时间信号作为一条数据记录放在样本数据库中；

步骤2.12、根据实际发生的典型事件类型对每条数据记录分别进行标注，构建不同类型事件的样本数据库；

所述步骤2.2的具体步骤如下：

步骤2.21、对每条数据记录提取N₁阶Mel频率倒谱系数；

步骤2.22、对每条数据记录提取N₂阶AR模型系数；

步骤2.23、将步骤2.21和步骤2.22提取的特征值合并进行特征选择；

所述步骤2.23的具体步骤如下：

将步骤2.21提取的C＝{C(i),i＝1,2,…N₁}和步骤2.22提取的A＝{A(j),j＝1,2,…,N₂}特征值合并，构成每条数据记录的N^*维特征向量{F(l),l＝1,2,…,N^*}＝{C(i),A(j),(i＝1,2,…,N₁,j＝1,2,…,N₂)}，对样本数据库中的Cnum条数据记录，都分别进行以上特征值提取的操作，得到典型泄漏事件训练数据的特征矩阵F＝{F(c)(l),(c＝1,2,…,Cnum,l＝1,2,…,N^*)}，特征矩阵F中的行代表数据记录序数，列代表特征维数；

对特征矩阵F＝{F(c)(l),(c＝1,2,…,Cnum,l＝1,2,…,N^*)}中的每一个特征列F_l＝{F_l(c),c＝1,2,…,Cnum}，分别按照：(F_l(c)-min(F_l))/(max(F_l)-min(F_l))进行最大/最小规范化，得到N^*维取值在0-1范围内的特征列F_l′＝{F_l′(c),c＝1,2,…,Cnum}，然后求其方差

其中Fsev是该列特征的平均值

对每一个特征列的方差值FT_l(l＝1,2,…,N^*)从大到小排序，选取其中方差较大的前N个特征，0<N≤N^*，所选特征在原特征矩阵F中的对应列的标号{n₁,n₂,…n_N}∈{1,2,…,N^*}，删除原特征矩阵F中未选择的特征列，得到一个Cnum×N维的优化特征矩阵F_s＝{F_s(c)(n),(c＝1,2,…,Cnum,n＝1,2,…,N)}；

所述步骤2.3的具体步骤如下：

步骤2.31、基于模糊C均值聚类算法对选取的特征矩阵进行二值化处理，具体如下：

通过模糊C均值聚类算法将提取和选择后的特征矩阵F_s＝{F_s(c)(n),(c＝1,2,…,Cnum,n＝1,2,…,N)}转换成二值化特征矩阵FA＝{FA(c)(n),(c＝1,2,…,Cnum,n＝1,2,…,N)}，具体操作如下：将某一列的特征值F_s(:,n)＝{F_s(c)(n),(c＝1,2,…,Cnum)}通过模糊C均值聚类算法分为两堆，每一堆各有一个中心点，较大的中心点的值为big(n)，较小中心点的值为sma(n)，与当前特征值F_s(c)(n)比较大小，若|F_s(c)(n)-big(n)|<|F_s(c)(n)-sma(n)|，则令FA(c)(n)＝1，反之则令FA(c)(n)＝0；另设一个与FA取值相反的二值化特征矩阵

后续关联规则挖掘时将综合两组二值化特征矩阵的分析结果；将二值化特征矩阵的每一列，作为后续数据挖掘事务的一项，记为：feature(n),n＝1,2,…,N；另外，保存特征矩阵F_s的N对特征聚类中心点{Cen(n)＝{big(n),sma(n)},n＝1,2,…,N}，作为后续在线检测和识别时的特征二值化参考依据；

步骤2.32、对不同类型事件标签进行二值化处理，具体如下：

将V类事件分别表示为一个V位的二值化数值，若当前数据记录c为第β类，则第β位为1，其余都为0；以此类推，Cnum条样本记录的事件类型标注将转化为一个Cnum×V的二值矩阵，最终得到二值化事件标签矩阵L_Cnum×V；

将二值化事件标签矩阵的每一列，也作为后续数据挖掘的一项，记为：label(β),β＝1,2,…,V；

步骤2.33、拼接特征矩阵与事件标签矩阵，具体如下：

将记为feature(n),n＝1,2,…,N的特征列与记为label(β),β＝1,2,…,V的事件标签列拼接，对已有的两个取值相反的二值化特征矩阵FA_Cnum×N和FB_Cnum×N，分别与二值化事件标签矩阵L_Cnum×V对应拼接，得到两个包含特征列和事件标签的二值化训练集：FA′_Cnum×(N+V),FB′_Cnum×(N+V)。

2.根据权利要求1所述的一种基于光纤传感数据挖掘的输气管道泄漏检测识别方法，其特征在于，所述步骤2.4的具体步骤如下：输入二值化特征训练集：FA′_Cnum×(N+V),FB′_Cnum×(N+V)，通过Apriori算法挖掘提取feature→label的关联规则，为每一条规则计算以下四个规则评价指标：支持度s(feature→label)＝σ(feature∪label)/σ(样本)，σ为样本计数值、置信度c(feature→label)＝σ(feature∪label)/σ(feature)、Kulc度量k(feature→label)＝(c(feature→label)+c(label→feature))/2和不平衡度IR(feature→label)＝|σ(feature)-σ(label)|/|σ(feature)+σ(label)-σ(feature∪label)|，设置相应的四个阈值：最小支持度mins、最小置信度minc、最小Kulc度量mink、最大不平衡度maxIR，其中，mins的值依据实际样本数量和样本组成设定；minc是评价规则的重要程度；mink的值与规则的正反两个置信度有关；不平衡度IR与Kulc度量一起使用，按照设定的阈值修剪不满足阈值要求的关联规则；

对正反两组二值化训练集输入，将得到两组有关特征参数与事件类别feature与label的重要关联规则作为输出，提取这些规则的前件，即所提取的关联规则中feature的集合，用RA(x)表示从二值化特征训练集A′_Cnum×(N+V)挖掘的每条规则前件，用RB(y)表示从二值化特征训练集B′_Cnum×(N+V)挖掘的每条规则前件，则按不同的事件类别β分别形成规则集，得到各类事件的规则集R＝{RA(β)(x),RB(β)(y),(β＝1,2,…,V；x＝1,2,…,X_β；y＝1,2,…,Y_β)}，X_β，Y_β分别为事件β所对应的两组规则前件个数，保存此规则集R完成关联规则分类器的构建，将其作为测试信号分类的参照特征和判决条件。

3.根据权利要求2所述的一种基于光纤传感数据挖掘的输气管道泄漏检测识别方法，其特征在于，所述步骤3的具体步骤如下：

对采集到的某空间点时序测试声波/振动信号按q毫秒分帧，每一帧信号XX^*＝{XX^*(m),m＝1,2,…,M}即为一条测试数据记录，包含

个采样点，其中，有U点的重叠区域，M取整数；

对当前数据记录提取N₁阶Mel频率倒谱系数{C(i),i＝1,2,…,N₁}，和N₂阶AR模型系数{A(j),j＝1,2,…,N₂}，组成N^*维特征向量F^*＝{F^*(l),l＝1,2,…,N^*}＝{C(i),A(j),(i＝1,2,…,N₁,j＝1,2,…,N₂)}，特征选择时，按照保存的N个特征标号{n₁,n₂,…n_N}，从特征向量F^*中提取这N维系数，得到N维特征向量F_s ^*＝{F_s ^*(n),n＝1,2,…,N}；

采用二值化方法和保存的N对聚类中心点{Cen(n)＝{big(n),sma(n)},n＝1,2,…,N}，对特征向量F_s ^*＝{F_s ^*(n),n＝1,2,…,N}的每一个值计算min{|F_s ^*(n)-big(n)|,|F_s ^*(n)-sma(n)|}，如果距较大中心点更近，则令FA^*(n)＝1，反之则令FA^*(n)＝0，得到测试数据记录的特征向量F_s ^*的二值化特征向量FA^*＝{FA^*(n),n＝1,2,…,N}；

通过比较二值化特征向量FA^*与关联规则分类器中的规则集R满足的关系来进行训练事件类型的识别与分类，具体如下：设FA^*中所有为1的序号集合为D_A，所有为0的序号集合为D_B，则FA^*中为1和0的序号集D＝{D_A＝{ε},D_B＝{γ},(FA^*(ε)＝1,FA^*(γ)＝0)}，对比分类器关联规则分类器中的规则集R＝{RA(β)(x),RB(β)(y),(β＝1,2,…,V；x＝1,2,…,X_β；y＝1,2,…,Y_β)}，设RA_β＝{RA(β)(x),x＝1,2,…,X_β}，RB_β＝{RB(β)(y),y＝1,2,…,Y_β}表示针对事件类型β的两组规则集，分别计算序号集D满足集合R中各类事件关联规则集的数量占比Count(β)：

依据上式统计结果，比较各类事件关联规则集的满足数量占比，计算最大值maxC＝max{Count(β),β＝1,2,…,V}，当最大值满足阈值，即maxC>thr时，thr为事件判定需满足规则占比的最小阈值，该测试信号所属的事件类别为Count(β)为最大时对应的那个事件类别标号，记为β_k；当最大值不满足阈值，则令β_k＝0，得到β_k的值后，即可进行事件类型判定；

若出现新的未知事件，作为待定事件类型β^*，继续补充事件样本数据库，作为后期更新关联规则分类器的依据。

Images