CN107425906A

CN107425906A - 面向地下管网安全监测的分布式光纤传感信号处理方法

Info

Publication number: CN107425906A
Application number: CN201710613247.4A
Authority: CN
Inventors: 吴慧娟; 张伟; 刘香荣; 饶云江
Original assignee: University of Electronic Science and Technology of China
Current assignee: University of Electronic Science and Technology of China
Priority date: 2017-07-25
Filing date: 2017-07-25
Publication date: 2017-12-01
Anticipated expiration: 2037-07-25
Also published as: CN107425906B; US10620038B2; US20180080812A1

Abstract

本发明公开了一种面向地下管网安全监测的分布式光纤传感信号处理方法，属于基础设施安全监测领域，旨在提高现有分布式光纤声音/振动传感系统在复杂应用条件下的智能检测与识别能力。本发明利用分布式光纤声音/振动传感系统沿探测光缆拾取整条线路的声音或振动信号；对整个监测范围内各空间点的声音或振动信号分别自定义提取基于信号相对量的短时特征和长时特征；基于采集信号特征和各自不同的背景噪声等先验知识构建和训练各空间点的贝叶斯识别分类网络；在线检测时，各空间点分别基于现场采集信号进行异常检测，检测到异常则提取相应特征，输入相对应的贝叶斯识别分类网络进行识别和分类。

Description

面向地下管网安全监测的分布式光纤传感信号处理方法

技术领域

一种面向地下管网安全监测的分布式光纤传感信号处理方法，用于地下管网及长输管道安全监测，属于地下管网、长输管道、通信光缆、电力线缆、周界及结构安全监测领域。

背景技术

城市地下管网如水管、气管、热力管道和油管等是城市的血管和命脉，它的安全运行是人民生命财产安全的重要保证。管道泄漏事故频发，青岛“11.22”输油管道泄漏爆燃，大连输油管线爆裂火灾、台湾高雄丙烷泄漏等事故就是惨痛的教训，分析其主要事故原因均是由外力破坏油气管道引起。另外，据统计美国39.6％，欧洲49.7％，中国58.7％的泄漏事故是由外力破坏或干扰引起，远大于由腐蚀老化造成的影响。因此，利用物联网技术对地下管网进行实时在线监测，对外力破坏等隐患事件进行及时预测和报警，防患于未然，是当前管网安全管理亟需要解决的重要问题。

相敏光时域反射技术(Φ-OTDR)作为分布式光纤传感技术的代表，利用光纤感测沿线环境中声波、振动等物理量的空间分布和时间变化信息，该技术长距离多点定位能力强，同时感测灵敏度高，光纤中无功能器件，寿命长，单端探测，工程施工和维护简便，因此是实现大范围环境安全监测的一种重要技术手段，在管道安全、光电缆安全、轨道安全、土木结构健康监测及周界安防等应用领域中发挥重要作用。Φ-OTDR技术按其传感解调方式主要分为两大类：1)基于直接或外差探测的强度解调方式，得到光电流信号与振动源产生的相位变化为非线性关系，I～cos(Δφ(t))，对事件响应敏感度低，常用于分布式振动传感；2)基于非平衡Mach-Zehnder、Micelson干涉等或其他方法进行线性解调，得到光电流信号与振动源产生的相位变化为线性关系，I～Δφ(t)，灵敏度高，可用于分布式声音和振动传感。

现有基于线性解调的分布式光纤声音/振动传感系统在实际应用中，因光纤感测的信号随距离衰减导致整条线路各点位采集的信号幅度分布不均匀，大范围监测干扰环境复杂，特别是不同段或不同点位埋设环境千差万别，例如，有些段埋在火车道附近，有些段埋在公路附近，有些段则埋在工厂或河流附近，不同点位感知的干扰源不同。目前仅利用某一点位的信号特征或识别网络进行整条线路不同干扰源条件下所有点位的信号检测识别，容易造成算法自适应能力差、系统误报率偏高的问题，因此本发明专利提出基于不同点位传感信号与干扰源的先验知识分别构建各空间点的识别网络，以及信号相对量的短时、长时特征提取方法，提高分布式光纤声音/振动传感系统大范围监测时，在全线各点位信号幅度分布不均匀和各区段埋设环境不一致时声音和振动目标源的分布式智能辨识能力。

发明内容

本发明的目的在于：解决现有分布式光纤声音/振动传感系统监测中，因整条光纤感测的信号随距离衰减导致各点位采集的信号幅度分布不均匀，以及大范围监测环境复杂，特别是不同段和不同点位埋设环境千差万别，造成算法不具备自适应能力、系统误报率偏高的问题；提供了一种面向地下管网安全监测的分布式光纤传感信号处理方法。

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

一种面向地下管网安全监测的分布式光纤传感信号处理方法，其特征在于，包括如下步骤：

步骤1、利用监测范围内布设的探测光缆，拾取探测光缆沿线上各空间点的声音或振动信号；

步骤2、采集整个监测范围内各空间点的声音或振动信号进行基于信号相对量自定义的短时特征提取；

步骤3、各空间点根据短时特征检测是否有异常事件，对有异常事件的信号进行长时特征累积计算，累积计算后进行离散化处理，再输入对应该空间点训练好的贝叶斯识别分类网络进行识别和分类。

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

步骤3.1、给定判断阈值，根据短时特征对各空间点声音或振动信号进行异常判断；

步骤3.2、若检测出某空间点发生异常事件，对该空间点从发生异常事件的时刻开始进行长时特征累积计算，累积计算后进行离散化处理；

步骤3.3、将离散处理后的结果输入对应该空间点训练好的贝叶斯识别分类网络进行识别和分类。

进一步，所述步骤3中，各空间点的贝叶斯识别分类网络构建的具体步骤如下：

步骤31、累积各空间点的声音或振动信号并进行预处理，预处理后得到的长时间段信号构建典型事件数据训练集；

步骤32、对典型事件数据训练集进行基于信号相对量自定义的短时特征提取，提取后累积计算长时特征，将典型事件数据训练集转化为对应的长时特征训练集；

步骤33、基于长时特征训练集中各长时特征的概率分布信息，结合各空间点背景环境的先验概率信息，构建贝叶斯识别分类网络拓扑结构并进行参数学习，得到训练好的贝叶斯识别分类网络。

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

步骤311、基于周期性触发脉冲采集各空间点的声音或振动信号，第k个触发脉冲周期内采集的整个监测范围的声音或振动信号为：

X_k＝{x_ki(i＝1,2,…,N)}＝[x_k1,x_k2,…,x_kN] (1)；

式中，N为整个监测范围内的数据采集长度；

步骤312、将每个光脉冲触发周期内采集的声音或振动信号沿时间轴进行累积，构建一个空间N维，时间T维的时空信号响应矩阵为：

XX＝{x_ki(k＝1,2,…,T；i＝1,2,…,N)} (2)；

步骤313、空间N维，时间T维的时空信号响应矩阵中，第i个空间点累积得到的短时间段信号为：

w_i＝{x_ki,k＝1,2,…,M} (3)；

步骤314、各空间点在时间轴上连续累积m个短时段信号构成一个长时间段信号为：

W_i＝{x_ki(k＝1,2,…,L)} (4)；

L＝m×M m≥1，

式中，L为长时间段信号长度,M为短时间段信号长度；

步骤315、累积空间长度为L的长时间段信号，按实际发生的事件类型分别给长时间段信号贴上事件类型标签，并添加到数据库中构建典型事件数据训练集。

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

步骤321、对短时间段信号提取4个短时特征，分别为：主冲击强度特征A₁、次冲击强度特征A₂、主次冲击幅值比例特征A₃和小波包能量谱特征A₄；

步骤322、在短时特征的基础上累积计算得到6个长时特征，分别为：长时主冲击强度特征B₁，长时次冲击强度特征B₂，长时幅值比例特征B₃，长时小波包能量谱特征B₄，长时环境背景特征B₅和长时严格冲击个数特征B₆；

步骤323、根据长时特征和其对应的事件类型标签，将典型事件数据训练集转换为对应的长时特征训练集。

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

主冲击强度特征A₁的提取方法：

将一个短时间段信号w_i作为处理对象，对该短时间段信号的幅度值进行直方图统计，将短时间段信号w_i的幅度取值区间[min(w_i),max(w_i)]平均分成10个连续的区间，其中，max(w_i)表示w_i的最大值，min(w_i)表示w_i的最小值，每一个区间记为s_l(l＝1,2,…,10)，即：

s_l＝[(l-1)×averlength+min(w_i),l×averlength+min(w_i)],(l＝1,2,…,10) (5)；

在公式(5)中，averlength为每一个区间的宽度，即：

分别计算短时间段信号w_i在相应的10个区间采样点的频数分布直方图及每个区间频数百分比percent，即：

式(7)中，n_l表示短时间段信号所有采样点的幅值在区间s_l的采样点个数，M为短时间段信号长度，假设percent_l(l∈[1,10])在l＝q处取得最大值，q为直方图最高处所对应的幅度区间，那么短时主冲击强度特征A₁为：

次冲击强度特征A₂的提取方法：

假设短时间段信号w_i＝{x_ki,k＝1,2,…,M}，在k＝t处取得绝对值最大值，t为短时间段信号绝对值最大时对应的时刻，将w_i在k＝t附近的M/4个采样点去除，得到剩余的3M/4个采样点构成时间序列w′_i，M为短时间段信号长度，即：

最后，以处理后的w′_i为对象，运用与提取主冲击强度特征A₁相同的办法，来提取次冲击强度特征A₂；

主次冲击幅值比例特征A₃的提取方法：

式(10)中，max(|w′_i|)表示时间序列信号w′_i的最大幅度绝对值，max(|w_i|)表示时间序列信号w_i的最大幅度绝对值；

小波包能量谱特征A₄的提取方法：

利用db6小波包对短时间段信号w_i进行三层小波包分解，将分解后的各子频带按照频率由低到高顺序，记为s₀～s₇，最后计算各频带的小波包系数能量，构成小波包系数能量特征向量E_wp＝[E₀,E₁,…,E₇]，可得：

式(11)中，E₂，E₃代表中低频子频段的能量，E₀，E₁代表低频子频段的能量。

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

步骤3221、初始化：令变量t＝1，t为短时间段信号绝对值最大时对应的时刻；令集合将长时间段信号看作m个短时间段信号w_i,t(t＝1,2,…,m)，w_i,t表示第t个短时间段信号；

步骤3222、将m个短时间段信号的主冲击强度特征A₁、次冲击强度特征A₂、主次冲击幅值比例特征A₃和小波包能量谱特征A₄分别记为a1、a2、a3、a4，若a1>thr，表示将该短时间段信号判断为异常事件，那么将a1、a2、a3和a4分别添加到集合T1、T2、T3、T4中；若a1<thr，表示将该短时间段信号判断为正常事件，则将a1放进集合T5中；若a1>thr且a2<thr1，a3<thr2，满足严格冲击条件，则将a1放进集合T6中，thr、thr1、thr2分别为检测冲击强度的经验阈值，对集合T1、T2、T3、T4、T5中的短时特征值分别求平均值，依次可以得到长时主冲击强度特征B₁，长时次冲击强度特征B₂，长时幅值比例特征B₃，长时小波包能量谱特征B₄，长时环境背景特征B₅的值为：

B₁＝average(T1),B₂＝average(T2),B₃＝average(T3),B₄＝average(T4) (12)；

B₅＝average(T5) (13)；

计算长时严格冲击个数特征B₆为：

B₆＝length(T6) (14)；

式(14)表示集合T6中数据个数，即该长时间段内满足短时严格冲击条件的个数。

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

步骤331、对长时特征训练集进行离散化处理，构建离散长时特征训练集；

步骤332、根据各空间点的离散长时特征训练集，获得所有长时特征属性节点与类节点的互信息，以及在类节点条件下的特征属性节点间的条件互信息，构建贝叶斯网络拓扑结构；

步骤333、构建各空间点的贝叶斯网络拓扑结构，基于构建好的贝叶斯网络拓扑结构进行参数学习，得到训练好的贝叶斯识别分类网络。

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

步骤3321、计算长时特征属性节点与类节点的互信息：

计算所有长时特征属性节点B_j与类节点C的互信息，即

式(15)中，S_b表示长时特征属性节点B_j(j＝1,2,…,6)所有可能取离散值集合，p(b'_j,c)表示B_j＝b'_j,C＝c的联合概率，p(b′_j)表示B_j＝b'_j的概率；p(b'_j,c)和p(b′_j)的值由离散长时特征训练集统计得到，计算方法为：p(b'_j,c)＝离散长时特征训练集中B_j＝b'_j,C＝c的记录条数/离散长时特征训练集中记录的总条数；p(b′_j)＝离散长时特征训练集中B_j＝b'_j的记录条数/离散长时特征训练集中记录的总条数；

步骤3322、计算在类节点条件下的特征属性节点间的条件互信息：

根据长时特征训练集，计算两两特征属性节点B_i,B_j(i≠j,且i＝1,2,…6,j＝1,2,…,6)在类节点C条件下的条件互信息，

式(16)中，b′_i表示长时特征属性节点B_i的值，n表示B_i可能取值的最大值，b'_j表示B_j的值，h表示B_j的值的最大值，c表示类节点C的值，p(b′_i,b'_j,c)表示B_i＝b′_i,B_j＝b'_j,C＝c的联合概率，p(b′_i,b′_j|c)、p(b′_i|c)以及p(b'_j|c)是相对应的条件概率；计算方法为：p(b′_i,b'_j,c)＝离散长时特征训练集中B_i＝b′_i,B_j＝b'_j,C＝c的记录条数/离散长时特征训练集中记录的总条数，p(b′_i,b′_j|c)＝离散长时特征训练集中B_i＝b′_i,B_j＝b'_j,C＝c的记录条数/离散长时特征训练集中C＝c的记录条数，p(b′_i|c)和p(b'_j|c)计算方法与p(b′_i,b′_j|c)的计算方法类似；

步骤3223、初始化：令集合集合S2＝{B₁,B₂,…,B₆}；在S2中选择一个长时特征属性节点B_j(j＝1,2,…,6)，使I(B_j；C)的值最大，添加一条弧线，由类节点C指向B_j，即类节点C是B_j的一个父节点；

步骤3224、在S2中选择一个长时特征属性节点B_j，使I(B_j；C)的值最大，添加一条弧线，由C指向B_j；令k＝min(|S1|,2)，在S1中选择前k个长时特征属性节点，使I(B_i；B_j|C)最大，为网络添加k条弧线，且由B_i指向B_j，将长时特征属性节点B_j添加到S1中，且将B_j从S2中删除，直到集合S1中包含所有的长时特征属性节点为止；

确定贝叶斯网络拓扑结构时，先将B₆从集合S2中删除，并将其添加到集合S1中；集合S2中剩余的节点为{B₁,B₂,B₃,B₄,B₅}，此时，集合S2中，B₂与类节点C的互信息最大，k＝min(|S1|,2)＝1，在S1中选择k个I(B_i；B₂|C)的节点作为B₂的父节点，即将B₆作为B₂的父节点，将B₂从集合S2中删除，并将B₂添加到集合S1中；这时，S2＝{B₁,B₃,B₄,B₅}，S1＝{B₂,B₆}，集合S2中，B₅与类节点C的互信息最大，k＝2，在S1中选择k个I(B_i；B₅|C)的节点作为B₅的父节点，并将B₅从集合S2中删除，并将B₅添加到集合S1中；如此循环，直到S2为空集合为止。

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

由后验概率定义式，分别计算离散长时特征为当前值时该事件属于真实异常事件的概率PS₁和属于环境干扰事件的概率PS₂：

式(17)中，p(c＝1)表示事件类别是真实异常事件的概率，是一个先验概率值；p(b′₁,…,b′₆|c＝1)是一个条件概率，表示在类别为真实异常事件的条件下，长时特征为当前值的概率；p(b′₁,…,b'₆)表示长时特征为当前值的概率；p(b′_i|pa_i,c＝1)表示在长时特征属性节点B_i的父节点离散值为pa_i且类别为真实异常事件的条件下，长时特征属性节点B_i的离散值为b′_i的概率，式(18)中各符号的物理意义与式(17)类似，由于PS₁和PS₂实际计算过程中p(b′₁,…,b′₆)的值是一样的，且p(b′₁,…,b′₆|c＝1)等价于p(b′_i|pa_i,c＝1)是训练过程得到的参数，PS₁和PS₂分别等效为PS'₁和PS'₂：

依据计算得到的后验概率值确定该事件的所属类别，如果某空间点的后验概率计算得到PS'₁>PS'₂，则判断该空间点发生的异常事件以PS'₁概率判断为真实异常事件，分类结果class＝1；否则以PS'₂概率判断为环境干扰事件，分类结果class＝2。

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

一、本发明中考虑了各空间点的环境背景差异建立分布式贝叶斯识别分类网络，真正实现各空间点差异化的分布式识别和分类，提高整个系统在大范围监测时全线各点位信号幅度分布不均匀和各区段埋设环境不一致等复杂条件下的正确检测率，并降低系统误报率；

二、本发明是一种典型全分布式光纤传感系统，单端探测，监测范围可达几十公里，大范围区域内监测成本低；光纤感知前端无源，不受天气、气候及光线等环境因素影响，感知效率高，寿命比电类传感网络长；

三、本发明可直接利用现有沿管道铺设的通信光缆空余纤芯进行长距离单端探测，施工方便，维护简单；

四、本发明的检测结果可由光时域反射原理确定其位置，该系统还具备多点异常的同时检测、识别和定位能力。

附图说明

图1为本发明的系统结构及方法流程示意图；

图2为本发明中基于贝叶斯识别分类网络的分布式光纤传感信号处理流程；

图3为本发明中分布式光纤传感信号累积及预处理：(a)分布式光纤传感信号累积及预处理流程；(b)累积及预处理得到的异常事件时空响应图；

图4为本发明中三种典型事件主冲击强度特征A₁值的概率分布：(a)正常事件；(b)真实异常事件；(c)环境干扰事件；

图5为本发明中某空间点长时间段信号的长时特征B₁、B₂、B₃、B₄、B₅的概率分布：(a)B₁的概率分布；(b)B₂的概率分布；(c)B₃的概率分布；(d)B₄的概率分布；(e)B₅的概率分布；

图6为本发明中某空间点长时间段信号的长时特征B₆的概率分布；

图7为某空间点各特征属性节点与类节点间的互信息；

图8为某空间点各特征属性节点在类节点下的条件互信息；

图9为本发明中某空间点构建的贝叶斯识别分类网络。

具体实施方式

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

实施例一

面向地下管网安全监测的分布式光纤传感信号处理方法的系统结构及其工作原理，如图1所示，基于线性解调Φ-OTDR技术原理的分布式光纤声音/振动传感系统，硬件主要由三个部分组成，探测光缆、光信号解调设备、信号处理主机。探测光缆通常采用普通单模通信光纤，一般沿地下管道、输电线缆、城镇道路埋地铺设，也可直接利用沿管道或道路铺设的通信光缆空余纤芯。光信号解调设备是该系统的核心，其内部组成器件主要包括光学器件和电学器件两类。由超窄线宽激光器产生一路连续相干光信号，经声光或电光调制器调制成光脉冲信号，光脉冲信号由掺铒光纤放大器(EDFA)集中放大，放大后的光脉冲信号依次经隔离器、环形器的1端口、2端口注入探测光缆；光脉冲信号沿光缆传输过程产生瑞利散射，其后向瑞利散射光信号沿光缆返回，由环形器的2端口、3端口接收，经光学滤波器滤除噪声信号后经第一耦合器耦合后，再将其注入到一个非平衡马赫曾德尔或迈克尔逊干涉仪等，具体由解调方法确定，经由3*3的第二耦合器输出三路相位差为120度的外界扰动引入的相位变化信息，即可获得声波及振动在光纤上的作用信号，解调出的光信号由光电探测器转换成电信号，再由波形发生卡控制的同步触发模数转换器进行信号同步采集，最后数字电信号通过网络等接口实时传输给信号处理主机。信号处理主机为普通电脑主机(PC)或FPGA/DSP嵌入式主板，用于光纤探测信号的分析、处理，通过特定信号处理算法得到引起声波、振动等的事件信息，并由光时域反射原理确定其位置，以及对感测事件进行智能识别和分类。

面向地下管网安全监测的分布式光纤传感信号处理方法：探测光纤本身相当于沿光纤均匀分布多个传感节点，每个传感节点对应于系统空间点一个采集节点，由这些分布式空间传感节点或采集节点协同工作拾取整条线路上的声音或振动信号，并对逐点采集的声音、振动信号分别进行可分辨特征提取，即本发明中特别定义的基于信号相对量的短时特征和长时特征提取，结合各空间点不同背景环境的事件概率等先验知识，构建分布式贝叶斯识别分类网络，最后，基于构建好的贝叶斯识别分类网络对各空间点感知信号提取的短时特征和长时特征进行异常事件的分布式检测，然后进行识别和分类，其中，采集整个监测范围内各空间点的声音或振动信号进行短时特征提取的方法与构建贝叶斯识别分类网络的方法中的短时特征提取方法相同。

实施例二

在实施例一的基础上，基于分布式贝叶斯识别分类网络的分布式光纤传感信号处理流程，如图2所示，主要分为两个部分：

第一部分，基于训练数据库构建各空间点的贝叶斯识别分类网络：首先累积各个空间的声音/振动传感信号并进行预处理，预处理后得到的长时间段信号构建典型事件数据训练集；对典型事件数据训练集进行基于信号相对量的短时特征提取，并通过提取的短时特征累积计算长时特征，离散化处理后基于离散长时特征训练集中各长时特征的概率分布信息，结合空间各点背景环境的先验概率信息，构建各自的贝叶斯识别分类网络拓扑并进行参数学习，得到训练好的分布式贝叶斯识别分类网络。

第二部分，基于现场事件信号进行在线检测、然后识别和分类：给定判断阈值，基于信号相对量提取的短时特征对各空间的声音或振动信号进行异常判断，若检测出某空间点发生异常事件，对该空间从发生异常事件的短时特征开始进行长时特征累积计算，累积计算后进行离散化处理，离散化处理后输入构建好的分布式贝叶斯识别分类网络对现场事件信号进行在线识别和分类。

实施例三

在实施例二的基础上，各空间点的光纤声音/振动传感信号累积及预处理过程，如图3(a)所示：光信号解调设备以一个触发脉冲作为一个信号采集周期，该时间单元内采集的信号是一条沿光缆分布的呈指数衰减的OTDR轨迹，设第k个触发脉冲周期内采集的整个监测范围的信号为：

X_k＝{x_ki(i＝1,2,…,N)}＝[x_k1,x_k2,…,x_kN] (1)；

其中，i表示空间的序号，N为整个监测范围内的数据空间采集点数，可以看作为N个空间传感节点。随着光脉冲周期触发，沿着时间轴将每个光脉冲触发周期内采集的原始OTDR轨迹进行累积，连续累积T条采集的原始OTDR轨迹，构建得到一个空间N维，时间T维的时空信号响应矩阵为：

XX＝{x_ki(k＝1,2,…,T；i＝1,2,…,N)} (2)；

T为第T个触发脉冲周期。如图3(b)所示，为时空信号响应矩阵XX的图形显示，即分布式光纤传感信号累积、预处理得到的异常事件时空响应图。本发明实施例中，对典型事件数据训练集中的每一条记录进行基于信号相对量的短时特征提取和长时特征累积计算，首先要进行短时间段信号和长时间段信号累积，短时间段信号和长时间段信号累积时间分别为1s和30s。设第i个空间或采样节点累积得到的短时间段信号定义为，

w_i＝{x_ki,k＝1,2,…,M} (3)；

本发明实施例中，M为时间轴上1s累积的数据点数，由于周期触发频率为508Hz，1s内短时间段信号数据点M为508。各空间点在时间轴上连续累积m个短时信号可构成一个长时间段信号，

W_i＝{x_ki(k＝1,2,…,L)} (4)；

L＝m×M m≥1，

式中，L为长时间段信号长度，m为正整数。本发明实施例中，长时间段信号累积长度为30s，故m取30。

实施例四

在实施例三的基础上，基于各空间点累积的典型事件长时间段信号，分别构建其典型事件数据训练集，具体操作过程如下：累积各空间点长度为L的长时间段信号，按实际发生的事件类型分别贴上事件类型标签。地下管网安全监测过程中，各空间点的典型事件类型通常有：没有异常事件的平稳环境噪声，人为挖掘，容易误判的交通干扰和工厂干扰等。本实施例中，将平稳环境噪声称为正常事件，事件类别标签设为0；将人为挖掘、交通和工厂干扰等事件统称为异常事件，其中人为挖掘为真实异常事件，事件类别标签设为1，而交通和工厂干扰则属于环境干扰事件，事件类别标签设为2。根据事件类型标签分别添加到数据库中，完成三类典型事件数据训练集构建。

实施例五

在实施例四的基础上，对典型事件数据训练集的每条数据记录进行本发明中基于信号相对量自定义的短时特征提取和长时特征累积计算，将典型事件数据训练集转化为对应的长时特征训练集，具体方法如下：

1、短时间段信号的短时特征提取方法

为避免大范围监测不同位置处光纤感知信号幅度不均匀对检测、识别的影响，本发明中提取的均是相对量特征。对于1s短时间段信号提取4个基于信号相对量自定义的短时特征分别为：主冲击强度特征A₁、次冲击强度特征A₂、主次冲击幅值比例特征A₃和小波包能量谱特征A₄。其中主冲击强度特征A₁，提取该段信号的最大冲击特性；次冲击强度特征A₂，提取除去最大冲击部分信号后剩余信号的冲击特性；主次冲击幅值比例特征A₃，提取最大冲击部分与去除最大冲击后剩余信号冲击幅度比；小波包能量谱特征A₄，提取冲击信号的重要频段能量谱特征。

(1)主冲击强度特征A₁的提取方法

鉴于正常事件通常属于平稳的随机噪声信号，而典型的异常事件信号通常伴随较大的振动冲击。本发明提取的第一个短时特征为主冲击强度特征A₁。提取短时特征时，将一个短时间段信号w_i作为处理对象，首先对该短时间段信号的幅度值进行直方图统计，将短时间段信号w_i的幅度取值区间[min(w_i),max(w_i)]平均分成10个连续的区间，其中，max(w_i)表示w_i的最大值，min(w_i)表示w_i的最小值，每一个区间记为s_l(l＝1,2,…,10)，即：

在公式(5)中，averlength为每一个区间的宽度，即：

公式(7)中，n_l表示短时间段信号所有采样点的幅值在区间s_l的采样点个数，M为短时间段信号长度，假设percent_l(l∈[1,10])在l＝q处取得最大值，q为直方图最高处所对应的幅度区间，那么短时主冲击强度特征A₁为：

对于平稳的环境噪声即正常事件的短时间段信号而言，其横向幅值分布区间较小，且每一个区间的采样点个数即直方图高度相差不明显，此时，计算得到的A₁值很小；而对于异常事件的短时间段信号而言，其横向幅值分布区间会异常增大，同时绝大部分采样点的幅值绝对值较小，所以大部分的采样点集中在中间的两个较小幅值分布区间，此时，计算得到的A₁值较大。因此将A₁作为第一个短时可分辨特征，即主冲击强度特征进行提取操作。如图4所示，三类典型事件的主冲击强度特征A₁值的概率分布存在显著差异，图4(a)中正常事件信号的A₁值偏小，分布中心在0.4附近，而图4(b)中真实异常事件和图4(c)中环境干扰事件信号的A₁值均偏大，分布中心均在0.9附近。

(2)次冲击强度特征A₂的提取方法

如图4所示，真实异常事件和环境干扰事件的两类异常事件短时间段信号，通常都具有较大的主冲击强度，因此继续提取第二个短时特征，次冲击强度特征A₂。经分析，人为挖掘的真实异常事件短时间段信号一般只有一个强度较大冲击，将该冲击去除后，剩余短时间段信号接近于平稳的环境噪声信号；而交通干扰等环境干扰事件短时间段信号，由于路面不平或石子影响常会伴有多个冲击强度不同的冲击，将其最大的冲击去除后，剩余信号仍然包含一个以上的冲击信号。次冲击强度特征A₂的提取方法如下：假设短时间段信号w_i＝{x_ki,k＝1,2,…,M}，在k＝t处取得绝对值最大值，t为短时间段信号绝对值最大时对应的时刻，将w_i在k＝t附近的M/4个采样点去除，得到剩余的3M/4个采样点构成时间序列w_i'，M为短时间段信号长度，即：

最后，以处理后的w_i'为对象，运用与提取主冲击强度特征A₁相同的办法，如式(5)、(6)、(7)、(8)所示，来提取次冲击强度特征A₂。通常，人为挖掘事件的短时间段信号的次冲击强度特征A₂的值较小；而交通和工厂干扰事件的短时间段信号次冲击强度特征A₂值较大。

(3)主次冲击幅值比例特征A₃的提取方法

由于去除最大冲击后的真实异常事件短时间段信号类似平稳的环境噪声，与最大冲击的幅值相比，其幅值较小；而去除最大冲击后的交通和工厂干扰事件的短时间段信号却有多个冲击强度不同的冲击，与最大冲击的幅值相比，其幅值较大。因此，作为第三个短时特征，主次冲击幅值比例特征A₃的计算方法：

公式(10)中，max(|w′_i|)表示时间序列信号w′_i的幅度最大绝对值，max(|w_i|)表示时间序列信号w_i的幅度最大绝对值。

(4)小波包能量谱特征A₄的提取方法

通常，在真实异常事件的瞬时冲击中，中低频信息比较丰富；而交通或者工厂等环境干扰事件的短时间段信号瞬时冲击中，低频信息较为丰富。基于此，利用db6小波包对短时间段信号w_i进行三层小波包分解，将分解后的各子频带按照频率由低到高顺序，记为s₀～s₇，最后计算各频带的小波包系数能量，构成小波包系数能量特征向量E_wp＝[E₀,E₁,…,E₇]。由此得到小波包能量谱特征A₄：

式(11)中，E₂，E₃代表中低频子频段的能量，E₀，E₁代表低频子频段的能量，真实异常事件的短时间段信号中低频信息丰富，所以A₄值偏大，而交通或工厂干扰的环境干扰事件短时信号的低频信息丰富，所以A₄值一般偏小。

2、长时间段信号的长时特征累积计算方法

短时间段信号是对事件的局部观察，基于短时特征进行事件分类识别准确性不高，因此提取短时特征后，继续累积长时间段信号及其特征对事件进行长时持续观察，基于更丰富的事件信息进行更为准确的识别。本发明中，提取的自定义长时间段信号的长时特征共计6个，具体包括：长时主冲击强度特征B₁，长时次冲击强度特征B₂，长时幅值比例特征B₃，长时小波包能量谱特征B₄，长时环境背景特征B₅和长时严格冲击个数特征B₆，均是以提取的自定义短时特征为基础进行累积计算得到。其中，长时主冲击强度特征B₁，表示所有含有异常事件的短时间段信号的主冲击强度的平均值；长时次冲击强度特征B₂，表示所有含有异常事件的短时间段信号的次冲击强度的平均值；长时幅值比例特征B₃，表示所有含有异常事件的短时间段信号的幅值比例特征平均值；长时小波包能量谱特征B₄，表示所有含有异常事件的短时间段信号的小波包能量谱特征的平均值；长时环境背景特征B₅，表示长时间段信号的背景起伏程度；长时严格冲击个数特征B₆，则表示长时间段信号内检测到的严格冲击个数。长时间段信号的长时特征无需重新提取，基于该长时间段信号已经得到的m个短时特征累积或平均计算得到，具体步骤为：

(a)初始化

令变量t＝1；令集合将该长时间段信号看作m个短时间段信号w_i,t(t＝1,2,…,m)，w_i,t表示第t个短时间段信号。

(b)重复以下操作，直到t＝m为止

将某一数据记录短时间段信号的短时特征A₁、A₂、A₃、A₄的值分别记为a1、a2、a3、a4。根据长时间段信号观察窗口内m个短时间段信号的四类短时特征值大小分别放置在不同的集合T1、T2、T3、T4、T5、T6中，为提取6个长时特征做准备。例如，若a1>thr，表示将该短时信号判断为异常事件，那么将a1、a2、a3和a4分别添加到集合T1、T2、T3、T4中；若a1<thr，表示将该短时间段信号判断为正常事件，则将a1放进集合T5中；若a1>thr且a2<thr1，a3<thr2，满足严格冲击条件，则将a1放进集合T6中。thr、thr1、thr2分别为判断冲击强弱的经验阈值。当m个短时间段信号的短时特征值全部放入集合T1、T2、T3、T4、T5、T6中，则开始计算该长时间段信号的6个长时特征值。其中，对集合T1、T2、T3、T4中的短时间段特征值分别求平均值，依次可以得到B₁、B₂、B₃、B₄的值：

集合T5中包含的是所有正常事件短时间段信号的a1。由于人为挖掘等真实异常事件是间断出现，其附近的短时间段信号较为平稳，a1值较小，而交通车辆、工厂干扰等环境干扰事件较为复杂，其附近平稳的短时间段信号也含有较小的波动，a1值较大，最后，对集合T5中所有的a1值求平均值，可得到长时间段信号的环境背景波动情况。因此，长时环境背景特征B₅计算为：

B₅＝average(T5) (13)；

由于环境干扰事件信号中的冲击与真实异常事件信号的冲击有明显差异，为尽量减少环境干扰事件引入的误报，需要对环境干扰事件产生的冲击信号进行严格筛选和排除，因此增加长时严格冲击个数特征B₆作为异常事件判断的依据之一。由单一经验阈值thr简单判断得到的冲击实际包含两类，伪冲击和严格冲击；继续增加两个经验阈值thr1、thr2判断筛选，得到严格冲击，作为排除伪冲击信号和检测识别真实异常事件的依据。具体操作为，当短时主冲击强度特征A₁的值a1>thr时，当且仅当短时次冲击强度特征A₂的值a2<thr1，且短时幅值比特征A₃的值a3<thr2的时候，才判断该短时间段信号w_i,t是严格冲击信号。其中，经验阈值thr1和阈值thr2的选取准则依据事件数据训练库的训练信号确定。本实施例中，对于人为挖掘引起的真实异常事件信号，通过阈值thr、thr1、thr2处理后，将以70％以上的概率被判断为严格冲击信号；而对于交通干扰和工厂干扰等环境异常事件信号，通过阈值thr、thr1、thr2处理后，仅有10％以下的概率被判断为严格冲击冲击信号。最后，长时严格冲击个数特征B₆计算为：

B₆＝length(T6) (14)；

3、构建典型事件数据训练集相对应的长时特征训练集

按照上述特征提取方法，对典型事件数据训练集中的每一条数据记录进行短时特征提取和长时特征累积计算，加上该条数据记录的事件类型标签，最终将得到该记录所对应的长时特征向量，记为Train_Feature＝[class,b1,b2,b3,b4,b5,b6]，其中class表示该条记录的事件类型标签，b1,b2,b3,b4,b5,b6分别为长时特征B₁，B₂，B₃，B₄，B₅，B₆的值。最后，将整个典型事件数据训练集转化为所对应的长时特征训练集TrainingFeature。

实施例六

在实施例五的基础上，对典型事件数据训练集中的长时特征训练集进行离散化处理，具体步骤如下：

1、长时特征属性值的离散化

由于特征集{B₁，B₂，B₃，B₄，B₅，B₆}中前5个特征属性值是连续的，而贝叶斯识别分类网络中的每一个属性节点的值是离散值，所以需要对长时特征的特征属性值进行离散化处理，具体方法：以某空间点为例，由该空间点的连续长时特征训练集，统计长时特征的概率分布，然后依据概率分布进行区间划分和离散化处理。

(1)B₁特征属性值的离散化处理

基于连续长时特征训练集计算得到某空间点长时间段信号的长时特征B₁的概率分布，如图5(a)所示，在区间[0.80,0.91]，环境干扰事件的概率密度大；在区间[0.91,0.92],两类异常事件的概率密度相差不大；在区间[0.92,1.0]，真实异常事件的概率密度较大。因此将区间(0,0.8),[0.80,0.91]、[0.91,0.92]、[0.92,1.0]上的B₁的长时特征值分别置为0、1、2、3。

(2)B₂特征属性值离散化处理

基于连续长时特征训练集计算得到某空间点长时间段信号的长时特征B₂的概率分布，如图5(b)所示，在区间[0,0.63)两类事件的概率密度相差不大；在区间[0.63,0.68)上，真实异常事件的概率密度大；在区间[0.68,1.0]上,环境干扰事件的概率密度较大。因此将区间[0,0.63)、[0.63,0.68)、[0.68,1.0]上的B₂属性特征值分别置为0、1、2。

(3)B₃特征属性值的离散化处理

基于连续长时特征训练集计算得到某空间点长时间段信号的长时特征B₃的概率分布，如图5(c)所示，在区间[0,0.32)上，真实异常事件的概率密度大；在区间[0.32,0.39)两类事件的概率密度相差不大；在区间[0.39,1.0]上，环境干扰事件的概率密度较大。因此将区间[0,0.32)、[0.32,0.39)、[0.39,1.0]上的B₃属性特征值分别置为0、1、2。

(4)B₄特征属性值的离散化处理

基于连续长时特征训练集计算得到某空间点长时间段信号的长时特征B₄的概率分布，如图5(d)所示,在区间上[0,0.6)，环境干扰事件的概率密度较大；在区间[0.6,0.8)两类事件的概率密度相差不大；在区间[0.8,3.0]上,环境干扰事件的概率密度较小，真实异常事件的概率密度较大。因此将区间[0,0.6)、[0.6,0.8)、[0.8,3.0]、(3.0,∞)的B₄属性特征值分别置为0、1、2、3。

(5)B₅特征属性值的离散化处理

基于连续长时特征训练集计算得到某空间点长时间段信号的长时特征B₅的概率分布，如图5(e)所示,在区间[0,0.59)上，真实异常事件的概率密度大；在区间[0.59,0.67)两类事件的概率密度相差不大；在区间[0.67,1.0]上，环境干扰事件的概率密度较大。因此发明实施例中将区间[0,0.59)、[0.59,0.67)、[0.67,1.0]上的B₅属性特征值分别置为0、1、2。

长时特征B₆是长时间段内对m个短时间段信号严格冲击的计数值，本身是离散的，无需再进行离散化处理。某空间长时间段信号的长时特征B₆的概率分布，如图6所示，可以直接用以构建贝叶斯识别分类网络。

2、构建离散的长时特征训练集

根据上述特征属性值离散化方法，对连续长时特征训练集中每一条数据记录所对应的特征值离散化处理，将实施例五得到的长时特征训练集TrainingFeature转化为离散长时特征训练集TrainingFeature'，其中离散长时特征训练集中的特征矢量记为：Train_Feature'＝[class,b1',b2',b3',b4',b5',b6']，b1'，b2'，b3'，b4'，b5'，b6'分别为长时特征B₁，B₂，B₃，B₄，B₅，B₆离散化处理后对应的特征值。

实施例七

在实施例六的基础上，构建各空间点(对应各个空间)的分布式贝叶斯识别分类网络，具体包含该空间点的贝叶斯网络拓扑结构学习及网络参数学习。

以某空间点的贝叶斯识别分类网络构建为例，依据该点的离散长时特征训练集获得所有长时特征属性节点与类节点的互信息，以及在类节点条件下的特征属性节点间的条件互信息，进行贝叶斯网络拓扑结构学习；基于构建好的贝叶斯网络拓扑结构进行参数学习，即在网络拓扑结构基础上，计算每一个特征属性节点的条件概率，具体步骤如下：

1、依据互信息和条件互信息进行贝叶斯网络拓扑结构学习

(1)计算长时特征属性节点与类节点的互信息

计算所有长时特征属性节点B_j与类节点C的互信息，即

式(15)中，S_b表示长时特征属性节点B_j(j＝1,2,…,6)所有可能取离散值集合，p(b'_j,c)表示B_j＝b'_j,C＝c的联合概率，p(b′_j)表示B_j＝b'_j的概率；p(b′_j,c)和p(b′_j)的值由离散长时特征训练集统计得到，计算方法为：p(b'_j,c)＝离散长时特征训练集中B_j＝b'_j,C＝c的记录条数/记录离散长时特征训练集的总条数；p(b′_j)＝离散长时特征训练集中B_j＝b'_j的记录条数/离散长时特征训练集中记录的总条数。本实施例中，根据上述方法计算所有离散长时特征属性节点B_j与类节点C的互信息I(B_j；C)，计算结果如图7所示。

(2)计算在类节点条件下特征属性节点间的条件互信息

根据离散长时特征训练集，计算两两特征属性节点B_i,B_j(i≠j,且i＝1,2,…6,j＝1,2,…,6)在类节点C条件下的条件互信息。

式(16)中，b′_i表示B_i的值，n表示B_i可能取值的最大值，b'_j表示B_j的值，h表示B_j的值的最大值，c表示类节点C的值，p(b′_i,b'_j,c)表示B_i＝b′_i,B_j＝b'_j,C＝c的联合概率，p(b′_i,b′_j|c),p(b′_i|c)以及p(b'_j|c)是相对应的条件概率；计算方法为：p(b′_i,b'_j,c)＝离散长时特征训练集中B_i＝b′_i,B_j＝b'_j,C＝c的记录条数/离散长时特征训练集中记录的总条数，p(b′_i,b′_j|c)＝离散长时特征训练集中B_i＝b′_i,B_j＝b'_j,C＝c的记录条数/离散长时特征训练集中C＝c的记录条数，p(b′_i|c)和p(b'_j|c)计算方法与p(b′_i,b′_j|c)的计算方法类似。本实施例中，根据上述方法计算在类节点条件下特征属性节点间的条件互信息I(B_i；B_j|C)，计算结果如图8所示。

(3)构建贝叶斯网络拓扑结构

根据长时特征属性节点与类节点的互信息，以及在类节点条件下特征属性节点间的条件互信息，确定贝叶斯网络拓扑结构，具体步骤如下：

(a)初始化

令集合集合S2＝{B₁,B₂,…,B₆}；在S2中选择一个特征属性节点B_j(j＝1,2,…,6)，使I(B_j；C)的值最大，添加一条弧线，由类节点C指向B_j，即类节点C是B_j的一个父节点；

(b)重复以下操作，确定贝叶斯网络拓扑结构

在S2中选择一个长时特征属性节点B_j，使I(B_j；C)的值最大，添加一条弧线，由C指向B_j；令k＝min(|S1|,2)，在S1中选择前k个长时特征属性节点，使I(B_i；B_j|C)最大，为网络添加k条弧线，且由B_i指向B_j，将长时特征属性节点B_j添加到S1中，且将B_j从S2中删除。直到集合S1中包含所有的长时特征属性节点为止。

确定贝叶斯网络拓扑结构时，由于B₆与类节点C的互信息最大，所以先将B₆从集合S2中删除，并将其添加到集合S1中；集合S2中剩余的节点为{B₁,B₂,B₃,B₄,B₅}，此时，集合S2中，B₂与类节点C节点的互信息最大，k＝min(|S1|,2)＝1，在S1中选择k个I(B_i；B₂|C)的节点作为B₂的父节点，即将B₆作为B₂的父节点，将B₂从集合S2中删除，并将其添加到集合S1中；这时，S2＝{B₁,B₃,B₄,B₅}，S1＝{B₂,B₆}，集合S2中，B₅与类节点C的互信息最大，k＝2，在S1中选择k个I(B_i；B₅|C)的节点作为B₅的父节点，并将B₅从集合S2中删除，并将其添加到集合S1中；如此循环，直到为空集合为止。本实施例中，依据图7和图8构建的贝叶斯网络拓扑结构，如图9所示。

2、贝叶斯网络参数学习

在构建的贝叶斯网络拓扑结构上进行网络参数学习，也就是计算贝叶斯网络拓扑结构每一个特征属性节点的条件概率p(b′_i|Pa_i,c)(i＝1,2,…,6)，b′_i表示长时特征属性节点B_i的值，Pa_i是由上述训练过程训练得到的属性节点B_i的父节点，计算方法为：p(b′_i|Pa_i,c)＝离散长时特征训练集中B_i＝b′_i,Pa_i＝pa_i,C＝c的记录条数/离散长时特征训练集中Pa_i＝pa_i,C＝c的记录条数，Pa_i为B_i的父节点Pa_i的值。

以上是以某空间点为例进行贝叶斯识别分类网络的构建和训练，以此类推，可以对光纤沿线所有空间点的贝叶斯识别分类网络分别进行分布式构建和训练。

实施例八

在实施例七的基础上，依据各空间点构建和训练好的分布式贝叶斯识别分类网络进行异常事件在线检测与识别过程：首先基于短时特征对各空间点的信号进行异常判断，若检测出某空间点发生异常事件，对该点进行长时特征累积计算和离散化处理，然后将离散化处理结果输入到构建好的分布式贝叶斯识别分类网络中，对该点异常事件进行识别和分类，即具体辨别该异常事件是人为挖掘引起的真实异常事件，还是交通或工厂等典型干扰引起的环境干扰事件。具体检测过程：对所有各空间点的采集信号沿时间轴进行短时(1s)信号累积，以1s短时间段信号为单元提取短时主冲击强度特征A₁进行异常事件检测。如果短时特征A₁的值小于经验阈值thr，则判断该短时间段信号属于正常事件，并遍历到下一空间点对其短时间段信号进行类似的短时特征提取和异常事件检测；若A₁超过经验阈值thr则判断该空间点发生异常事件。

实施例九

在实施例八的基础上，基于各空间点训练好的贝叶斯识别分类网络对现场检测到的异常事件进行在线识别和分类的具体过程：检测到空间某空间点发生异常事件后，从此刻开始对空间该点信号进行长时(30s)累积，当某一空间点的异常信号长时累积时间到，则类似于训练阶段累积计算相应的长时特征B₁，B₂，B₃，B₄，B₅，B₆，并对其进行离散化处理，输入到针对该点训练好的分布式贝叶斯识别分类网络中，对该异常事件进行识别和分类，输出该异常事件的判断类型，如真实异常事件、环境干扰事件等。具体步骤如下：

由后验概率定义式，分别计算长时特征为当前值时该事件属于真实异常事件的概率PS₁和属于环境干扰事件的概率PS₂：

式(17)中，p(c＝1)表示事件类别是真实异常事件的概率，是一个先验概率值；p(b′₁,…,b′₆|c＝1)是一个条件概率，表示在类别为真实异常事件的条件下，长时特征为当前值的概率；p(b′₁,…,b′₆)表示长时特征为当前值的概率；p(b′_i|pa_i,c＝1)表示在长时特征属性节点B_i的父节点离散值为pa_i且类别为真实异常事件的条件下，长时特征属性节点B_i的离散值为b_i'的概率。式(18)中各符号的物理意义与式(17)类似，式(17)为真实异常事件后验概率，式(18)为环境干扰事件的后验概率，由于PS₁和PS₂实际计算过程中p(b′₁,…,b'₆)的值是一样的，且p(b′₁,…,b′₆|c＝1)等价于p(b′_i|pa_i,c＝1)是训练过程得到的参数，PS₁和PS₂分别等效为PS'₁和PS'₂：

依据计算得到的后验概率值确定该事件的所属类别。例如，如果某空间点的后验概率计算得到PS'₁>PS'₂，则判断该点发生的异常事件以PS'₁概率判断为真实异常事件，分类结果class＝1；否则以PS'₂概率判断为环境干扰事件，分类结果class＝2。

本发明专利公开了一种面向地下管网安全监测的分布式光纤传感信号处理方法，不仅适用于基于相敏光时域反射技术的分布式光纤声波/振动传感系统，还适用于基于M-Z、Micelson、Sagnac干涉等其他分布式光纤声音/振动传感系统,以及基于光纤光栅、光纤珐玻阵列的准分布式光纤声波/振动传感系统；该方法中的短时、长时特征提取和离散化具体方法可根据实际应用选取。本发明基于各空间点感知的真实异常和环境干扰信号构建分布式识别分类网络的方法，可直接应用于地下管网如油气管道、供水管道、热力管道等的安全监测应用中，也适用于通信光缆、电力线缆、周界安防、结构健康监测等分布式光纤传感应用领域。以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims

1.一种面向地下管网安全监测的分布式光纤传感信号处理方法，其特征在于，包括如下步骤：

2.根据权利要求1所述的一种面向地下管网安全监测的分布式光纤传感信号处理方法，其特征在于，所述步骤3的具体步骤如下：

3.根据权利要求1、2所述的一种面向地下管网安全监测的分布式光纤传感信号处理方法，其特征在于，所述步骤3中，各空间点的贝叶斯识别分类网络构建的具体步骤如下：

4.根据权利要求3所述的一种面向地下管网安全监测的分布式光纤传感信号处理方法，其特征在于，所述步骤31的具体步骤如下：

X_k＝{x_ki(i＝1,2,…,N)}＝[x_k1,x_k2,…,x_kN] (1)；

式中，N为整个监测范围内的数据采集长度；

XX＝{x_ki(k＝1,2,…,T；i＝1,2,…,N)} (2)；

w_i＝{x_ki,k＝1,2,…,M} (3)；

W_i＝{x_ki(k＝1,2,…,L)} (4)；

L＝m×M m≥1，

式中，L为长时间段信号长度,M为短时间段信号长度；

5.根据权利要求3所述的一种面向地下管网安全监测的分布式光纤传感信号处理方法，其特征在于，所述步骤32的具体步骤如下：

6.根据权利要求5所述的一种面向地下管网安全监测的分布式光纤传感信号处理方法，其特征在于，所述步骤321的具体步骤如下：

主冲击强度特征A₁的提取方法：

在公式(5)中，averlength为每一个区间的宽度，即：

<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>p</mi> <mi>e</mi> <mi>r</mi> <mi>c</mi> <mi>e</mi> <mi>n</mi> <mi>t</mi> <mo>=</mo> <mo>{</mo> <msub> <mi>percent</mi> <mi>l</mi> </msub> <mo>,</mo> <mi>l</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mn>...</mn> <mo>,</mo> <mn>10</mn> <mo>}</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>percent</mi> <mi>l</mi> </msub> <mo>=</mo> <mfrac> <msub> <mi>n</mi> <mi>l</mi> </msub> <mi>M</mi> </mfrac> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

<mrow> <msub> <mi>A</mi> <mn>1</mn> </msub> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>p</mi> <mi>q</mi> </msub> <mo>+</mo> <msub> <mi>p</mi> <mrow> <mi>q</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> </mrow> </mtd> <mtd> <mrow> <mo>(</mo> <mi>q</mi> <mo>&Element;</mo> <mo>&lsqb;</mo> <mn>2</mn> <mo>,</mo> <mn>9</mn> <mo>&rsqb;</mo> <mo>,</mo> <msub> <mi>p</mi> <mrow> <mi>q</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>&le;</mo> <msub> <mi>p</mi> <mrow> <mi>q</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>P</mi> <mi>q</mi> </msub> <mo>+</mo> <msub> <mi>P</mi> <mrow> <mi>q</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> </mrow> </mtd> <mtd> <mrow> <mo>(</mo> <mi>q</mi> <mo>&Element;</mo> <mo>&lsqb;</mo> <mn>2</mn> <mo>,</mo> <mn>9</mn> <mo>&rsqb;</mo> <mo>,</mo> <msub> <mi>P</mi> <mrow> <mi>q</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>></mo> <msub> <mi>P</mi> <mrow> <mi>q</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>p</mi> <mi>q</mi> </msub> <mo>+</mo> <msub> <mi>p</mi> <mrow> <mi>q</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> </mrow> </mtd> <mtd> <mrow> <mo>(</mo> <mi>q</mi> <mo>=</mo> <mn>1</mn> <mo>)</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>p</mi> <mi>q</mi> </msub> <mo>+</mo> <msub> <mi>p</mi> <mrow> <mi>q</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> </mrow> </mtd> <mtd> <mrow> <mo>(</mo> <mi>q</mi> <mo>=</mo> <mn>10</mn> <mo>)</mo> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

次冲击强度特征A₂的提取方法：

假设短时间段信号w_i＝{x_ki,k＝1,2,…,M}，在k＝t处取得绝对值最大值，t为短时间段信号绝对值最大时对应的时刻，将w_i在k＝t附近的M/4个采样点去除，得到剩余的3M/4个采样点构成时间序列w'_i，M为短时间段信号长度，即：

<mrow> <msubsup> <mi>w</mi> <mi>i</mi> <mo>&prime;</mo> </msubsup> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mi>i</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mn>...</mn> <mo>,</mo> <mi>t</mi> <mo>-</mo> <mfrac> <mi>M</mi> <mn>8</mn> </mfrac> <mo>,</mo> <mi>t</mi> <mo>+</mo> <mfrac> <mi>M</mi> <mn>8</mn> </mfrac> <mo>,</mo> <mi>t</mi> <mo>+</mo> <mfrac> <mi>M</mi> <mn>8</mn> </mfrac> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mn>...</mn> <mo>,</mo> <mi>M</mi> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mfrac> <mi>M</mi> <mn>8</mn> </mfrac> <mo><</mo> <mi>t</mi> <mo><</mo> <mfrac> <mrow> <mn>7</mn> <mi>M</mi> </mrow> <mn>8</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mi>i</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>=</mo> <mfrac> <mi>M</mi> <mn>4</mn> </mfrac> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mfrac> <mi>M</mi> <mn>4</mn> </mfrac> <mo>+</mo> <mn>2</mn> <mo>,</mo> <mn>...</mn> <mo>,</mo> <mi>M</mi> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mi>t</mi> <mo>&le;</mo> <mfrac> <mi>M</mi> <mn>8</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mi>i</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mn>...</mn> <mo>,</mo> <mfrac> <mrow> <mn>3</mn> <mi>M</mi> </mrow> <mn>4</mn> </mfrac> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mi>t</mi> <mo>&GreaterEqual;</mo> <mfrac> <mrow> <mn>7</mn> <mi>M</mi> </mrow> <mn>8</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

最后，以处理后的w'_i为对象，运用与提取主冲击强度特征A₁相同的办法，来提取次冲击强度特征A₂；

主次冲击幅值比例特征A₃的提取方法：

<mrow> <msub> <mi>A</mi> <mn>3</mn> </msub> <mo>=</mo> <mfrac> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> <mrow> <mo>(</mo> <mo>|</mo> <msubsup> <mi>w</mi> <mi>i</mi> <mo>&prime;</mo> </msubsup> <mo>|</mo> <mo>)</mo> </mrow> </mrow> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> <mrow> <mo>(</mo> <mo>|</mo> <msub> <mi>w</mi> <mi>i</mi> </msub> <mo>|</mo> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

式(10)中，max(|w'_i|)表示时间序列信号w'_i的最大幅度绝对值，max(|w_i|)表示时间序列信号w_i的最大幅度绝对值；

小波包能量谱特征A₄的提取方法：

7.根据权利要求5所述的一种面向地下管网安全监测的分布式光纤传感信号处理方法，其特征在于，所述步骤322的具体步骤如下：

B₅＝average(T5) (13)；

计算长时严格冲击个数特征B₆为：

B₆＝length(T6) (14)；

8.根据权利要求3所述的一种面向地下管网安全监测的分布式光纤传感信号处理方法，其特征在于，所述步骤33的具体步骤如下：

9.根据权利要求7所述的一种面向地下管网安全监测的分布式光纤传感信号处理方法，其特征在于，所述步骤332的具体步骤如下：

步骤3321、计算长时特征属性节点与类节点的互信息：

计算所有长时特征属性节点B_j与类节点C的互信息，即

<mrow> <mi>I</mi> <mrow> <mo>(</mo> <msub> <mi>B</mi> <mi>j</mi> </msub> <mo>;</mo> <mi>C</mi> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mo>&Sigma;</mo> <mrow> <mi>b</mi> <mo>&Element;</mo> <msub> <mi>S</mi> <mi>b</mi> </msub> </mrow> </munder> <munder> <mo>&Sigma;</mo> <mrow> <mi>c</mi> <mo>&Element;</mo> <mo>{</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>}</mo> </mrow> </munder> <mi>p</mi> <mrow> <mo>(</mo> <msubsup> <mi>b</mi> <mi>j</mi> <mo>&prime;</mo> </msubsup> <mo>,</mo> <mi>c</mi> <mo>)</mo> </mrow> <mi>log</mi> <mfrac> <mrow> <mi>p</mi> <mrow> <mo>(</mo> <msubsup> <mi>b</mi> <mi>j</mi> <mo>&prime;</mo> </msubsup> <mo>,</mo> <mi>c</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>p</mi> <mrow> <mo>(</mo> <msubsup> <mi>b</mi> <mi>j</mi> <mo>&prime;</mo> </msubsup> <mo>)</mo> </mrow> <mi>p</mi> <mrow> <mo>(</mo> <mi>c</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>,</mo> <mrow> <mo>(</mo> <mi>j</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mn>...</mn> <mo>,</mo> <mn>6</mn> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>15</mn> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

式(15)中，S_b表示长时特征属性节点B_j(j＝1,2,…,6)所有可能取离散值集合，p(b'_j,c)表示B_j＝b'_j,C＝c的联合概率，p(b'_j)表示B_j＝b'_j的概率；p(b'_j,c)和p(b'_j)的值由离散长时特征训练集统计得到，计算方法为：p(b'_j,c)＝离散长时特征训练集中B_j＝b'_j,C＝c的记录条数/离散长时特征训练集中记录的总条数；p(b'_j)＝离散长时特征训练集中B_j＝b'_j的记录条数/离散长时特征训练集中记录的总条数；

式(16)中，b'_i表示长时特征属性节点B_i的值，n表示B_i可能取值的最大值，b'_j表示B_j的值，h表示B_j的值的最大值，c表示类节点C的值，p(b'_i,b'_j,c)表示B_i＝b'_i,B_j＝b'_j,C＝c的联合概率，p(b'_i,b'_j|c)、p(b'_i|c)以及p(b'_j|c)是相对应的条件概率；计算方法为：p(b'_i,b'_j,c)＝离散长时特征训练集中B_i＝b'_i,B_j＝b'_j,C＝c的记录条数/离散长时特征训练集中记录的总条数，p(b'_i,b'_j|c)＝离散长时特征训练集中B_i＝b'_i,B_j＝b'_j,C＝c的记录条数/离散长时特征训练集中C＝c的记录条数，p(b'_i|c)和p(b'_j|c)计算方法与p(b'_i,b'_j|c)的计算方法类似；

10.根据权利要求2所述的一种面向地下管网安全监测的分布式光纤传感信号处理方法，其特征在于，所述步骤3.3具体如下：

<mrow> <msub> <mi>PS</mi> <mn>1</mn> </msub> <mo>=</mo> <mfrac> <mrow> <mi>p</mi> <mrow> <mo>(</mo> <mi>c</mi> <mo>=</mo> <mn>1</mn> <mo>)</mo> </mrow> <mi>p</mi> <mrow> <mo>(</mo> <msubsup> <mi>b</mi> <mn>1</mn> <mo>&prime;</mo> </msubsup> <mo>,</mo> <mn>...</mn> <mo>,</mo> <msubsup> <mi>b</mi> <mn>6</mn> <mo>&prime;</mo> </msubsup> <mo>|</mo> <mi>c</mi> <mo>=</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> <mrow> <mi>p</mi> <mrow> <mo>(</mo> <msubsup> <mi>b</mi> <mn>1</mn> <mo>&prime;</mo> </msubsup> <mo>,</mo> <mn>...</mn> <mo>,</mo> <msubsup> <mi>b</mi> <mn>6</mn> <mo>&prime;</mo> </msubsup> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>&Proportional;</mo> <mi>p</mi> <mrow> <mo>(</mo> <mi>c</mi> <mo>=</mo> <mn>1</mn> <mo>)</mo> </mrow> <munderover> <mo>&Pi;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>6</mn> </munderover> <mi>p</mi> <mrow> <mo>(</mo> <msubsup> <mi>b</mi> <mi>i</mi> <mo>&prime;</mo> </msubsup> <mo>|</mo> <msub> <mi>pa</mi> <mi>i</mi> </msub> <mo>,</mo> <mi>c</mi> <mo>=</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>17</mn> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

<mrow> <msub> <mi>PS</mi> <mn>2</mn> </msub> <mo>=</mo> <mfrac> <mrow> <mi>p</mi> <mrow> <mo>(</mo> <mi>c</mi> <mo>=</mo> <mn>2</mn> <mo>)</mo> </mrow> <mi>p</mi> <mrow> <mo>(</mo> <msubsup> <mi>b</mi> <mn>1</mn> <mo>&prime;</mo> </msubsup> <mo>,</mo> <mn>...</mn> <mo>,</mo> <msubsup> <mi>b</mi> <mn>6</mn> <mo>&prime;</mo> </msubsup> <mo>|</mo> <mi>c</mi> <mo>=</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow> <mrow> <mi>p</mi> <mrow> <mo>(</mo> <msubsup> <mi>b</mi> <mn>1</mn> <mo>&prime;</mo> </msubsup> <mo>,</mo> <mn>...</mn> <mo>,</mo> <msubsup> <mi>b</mi> <mn>6</mn> <mo>&prime;</mo> </msubsup> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>&Proportional;</mo> <mi>p</mi> <mrow> <mo>(</mo> <mi>c</mi> <mo>=</mo> <mn>2</mn> <mo>)</mo> </mrow> <munderover> <mo>&Pi;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>6</mn> </munderover> <mi>p</mi> <mrow> <mo>(</mo> <msubsup> <mi>b</mi> <mi>i</mi> <mo>&prime;</mo> </msubsup> <mo>|</mo> <msub> <mi>pa</mi> <mi>i</mi> </msub> <mo>,</mo> <mi>c</mi> <mo>=</mo> <mn>2</mn> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>18</mn> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

式(17)中，p(c＝1)表示事件类别是真实异常事件的概率，是一个先验概率值；p(b'₁,…,b'₆|c＝1)是一个条件概率，表示在类别为真实异常事件的条件下，长时特征为当前值的概率；p(b'₁,…,b'₆)表示长时特征为当前值的概率；p(b'_i|pa_i,c＝1)表示在长时特征属性节点B_i的父节点离散值为pa_i且类别为真实异常事件的条件下，长时特征属性节点B_i的离散值为b'_i的概率，式(18)中各符号的物理意义与式(17)类似，由于PS₁和PS₂实际计算过程中p(b'₁,…,b'₆)的值是一样的，且p(b'₁,…,b'₆|c＝1)等价于p(b'_i|pa_i,c＝1)是训练过程得到的参数，PS₁和PS₂分别等效为PS'₁和PS'₂：

<mrow> <msubsup> <mi>PS</mi> <mn>1</mn> <mo>&prime;</mo> </msubsup> <mo>=</mo> <mi>p</mi> <mrow> <mo>(</mo> <mi>c</mi> <mo>=</mo> <mn>1</mn> <mo>)</mo> </mrow> <munderover> <mo>&Pi;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>6</mn> </munderover> <mi>p</mi> <mrow> <mo>(</mo> <msubsup> <mi>b</mi> <mi>i</mi> <mo>&prime;</mo> </msubsup> <mo>|</mo> <msub> <mi>pa</mi> <mi>i</mi> </msub> <mo>,</mo> <mi>c</mi> <mo>=</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>19</mn> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

<mrow> <msubsup> <mi>PS</mi> <mn>2</mn> <mo>&prime;</mo> </msubsup> <mo>=</mo> <mi>p</mi> <mrow> <mo>(</mo> <mi>c</mi> <mo>=</mo> <mn>2</mn> <mo>)</mo> </mrow> <munderover> <mo>&Pi;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>6</mn> </munderover> <mi>p</mi> <mrow> <mo>(</mo> <msubsup> <mi>b</mi> <mi>i</mi> <mo>&prime;</mo> </msubsup> <mo>|</mo> <msub> <mi>pa</mi> <mi>i</mi> </msub> <mo>,</mo> <mi>c</mi> <mo>=</mo> <mn>2</mn> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>20</mn> <mo>)</mo> </mrow> <mo>;</mo> </mrow>