CN108870091B - 基于高低频混合检测的管道泄漏监测系统及方法 - Google Patents

基于高低频混合检测的管道泄漏监测系统及方法 Download PDF

Info

Publication number
CN108870091B
CN108870091B CN201810794586.1A CN201810794586A CN108870091B CN 108870091 B CN108870091 B CN 108870091B CN 201810794586 A CN201810794586 A CN 201810794586A CN 108870091 B CN108870091 B CN 108870091B
Authority
CN
China
Prior art keywords
frequency
signal
data
head end
low
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810794586.1A
Other languages
English (en)
Other versions
CN108870091A (zh
Inventor
马大中
张化光
于洋
冯健
刘金海
关勇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northeastern University China
Original Assignee
Northeastern University China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northeastern University China filed Critical Northeastern University China
Priority to CN201810794586.1A priority Critical patent/CN108870091B/zh
Publication of CN108870091A publication Critical patent/CN108870091A/zh
Application granted granted Critical
Publication of CN108870091B publication Critical patent/CN108870091B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F17—STORING OR DISTRIBUTING GASES OR LIQUIDS
    • F17D—PIPE-LINE SYSTEMS; PIPE-LINES
    • F17D5/00—Protection or supervision of installations
    • F17D5/02—Preventing, monitoring, or locating loss
    • F17D5/06—Preventing, monitoring, or locating loss using electric or acoustic means

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Acoustics & Sound (AREA)
  • Mechanical Engineering (AREA)
  • General Engineering & Computer Science (AREA)
  • Complex Calculations (AREA)

Abstract

本发明提出一种基于高低频混合检测的管道泄漏监测系统及方法，包括：首端超高频传感器、末端超高频传感器、首端加速度传感器、末端加速度传感器、首端下位机、末端下位机、首端上位机和末端上位机；还提供一种基于该系统的方法：采集高频低频混合历史数据并已知其工况分类，进行预处理，多信息域分析处理，特征提取，用协方差核切片逆回归算法处理，得到特征向量子空间，采集高频低频混合实时数据，按照历史数据的步骤得到特征向量子空间，用线性判别分析进行分类，并计算相似度，得到实时数据分类结果，若有泄漏计算泄漏位置，为排除故障提供可靠的数据依据。确保了泄漏的特征值的准确性，最大限度的减少误报，对小流量的泄漏有更准确的判断。

Description

基于高低频混合检测的管道泄漏监测系统及方法
技术领域
本发明涉及管道风险预测技术领域，尤其涉及一种基于高低频混合检测的管道泄漏检测系统及方法。
背景技术
管道运输在经济发展中的所扮演的角色越来与重要，如城市自来水管道、陆地原油管道、海底油气管道等，石油的运输大部分以成品油的形式在管道中运输。随着管网的逐年扩建，管道运输己经成为陆上油气运输的主要方式。但是管道的老化、锈蚀、突发性自然灾害及人为破坏等都会造成成品油管道的泄漏乃至破裂，如不及时发现并加以制止，不仅造成能源浪费、经济损失、环境污染，而且会危及人身安全，甚至造成灾难性事故。因此，对油气管道进行实时在线监测，对泄漏事故进行准确及时的报警，并准确估计出泄漏点的位置具有重要的意义。
现今的管道泄漏检测方法有很多种，如光纤测漏法、声波检测法、压力梯度法、负压波法等等。其中，负压波法是近年来国际上应用最广的管道泄漏检测方法，该方法具有反应时间短、可检测泄漏量范围广等特点，但是对于缓慢且流量较小的泄漏，由于其单位时间内压力变化缓慢，负压波法对其敏感度较低，容易产生漏报，而且由于管道输送系统的复杂工况调整，一些常见的操作如主输泵的启停、阀门的开关、调节阀开度的变化等都会引起负压波，且与泄漏引发的负压波具有很高的相似度，降低了负压波法的检测精度。
声波检测法泄漏检测对于缓慢且流量较小的泄漏有较好的监测精度，但是其对易将泄漏和人为或环境因素造成的声波异常相混淆，并且这种方法不能探测同时发生的多点泄漏。
发明内容
本发明要解决的技术问题是针对上述现有技术的不足，提供一种基于高低频混合检测的管道泄漏检测系统及方法，能更好地屏蔽噪声干扰，并能确保在信号源改变时，滤波后还原信号的准确性，最大限度的减少误报，对小流量的泄漏有更准确的判断。
一种基于高低频混合检测的管道泄漏监测系统，包括首端超高频传感器、首端加速度传感器、首端下位机、首端上位机、末端超高频传感器、末端加速度传感器、末端下位机和末端上位机；
所述首端超高频传感器和首端加速度传感器分别安装在待测试管道的首端，首端超高频传感器和首端加速度传感器分别通过电线与首端下位机相连接，首端下位机通过网线与首端上位机相连接；
所述末端超高频传感器和末端加速度传感器分别安装在待测试管道的末端，末端超高频传感器和末端加速度传感器分别通过电线与末端下位机相连接，末端下位机通过网线与末端上位机相连接；
所述首端超高频传感器和首端加速度传感器，与待测试管道中输送介质相接触，首端超高频传感器用于采集首端管道高频信号的实时数据，首端加速度传感器用于采集低频信号的实时数据；
所述首端下位机用于存储从首端超高频传感器和首端加速度传感器传递出的超高频信号和低频信号的混合实时数据，将该高低频混合实时数据进行模数转换且加上时间戳，传递给首端上位机；
所述首端上位机用于接收并保存由首端下位机传递出来的加上时间戳的高低频混合实时数据，并将该高低频混合实时数据保存到首端上位机数据库中，首端上位机通过编程来实现管道泄漏检测功能；
所述末端超高频传感器和末端加速度传感器，与待测试管道中输送介质相接触，末端超高频传感器用于采集末端管道高频信号的实时数据，末端加速度传感器用于采集低频信号的实时数据；
所述末端下位机用于存储从末端超高频传感器和末端加速度传感器传递出的道高频信号和低频信号的高低频混合实时数据，将该高低频混合实时数据进行模数转换且加上时间数据，传递给末端上位机；
所述末端上位机用于接收并保存由末端下位机传递出来的加上时间戳的高低频混合实时数据，并将该高低频混合实时数据保存到首端上位机数据库中，末端上位机通过编程来实现管道泄漏检测功能；
所述首端下位机和末端下位机均包括：A-CORE模块、B-CORE模块、电源模块、存储模块、GPS校时模块、通讯模块、调试模块、传感器滤波电路、AD转换模块；
所述电源模块、存储模块、GPS校时模块、通讯模块、调试模块分别与A-CORE模块相连接，其中，存储模块通过并行总线与A-CORE模块相连接，GPS校时模块通过USART或SPI通讯协议与A-CORE模块相连接，通讯模块通过以太网RS-485协议与A-CORE模块相连接；传感器滤波电路与AD转换模块相连接，AD转换模块与B-CORE模块相连接；A-CORE模块通过网线与B-CORE模块相连接；电源模块与A-CORE模块和B-CORE模块相连接；
所述A-CORE模块是主控模块，将B-CORE模块传递过来的数据进行存储，并加以时间戳转发给首端或者末端上位机；
所述B-CORE模块是A-CORE模块的从设备，接收超高频传感器和加速度传感器传递出的高频信号和低频信号的高低频混合实时数据，并通过传感器滤波电路完成对该实时数据的硬件滤波，并且通过AD转换模块完成模数信号的转换，然后将模数信号转换后的数据传递给A-CORE模块；
所述电源模块为A-CORE模块和B-CORE模块提供电源；
所述存储模块保存用GPS校时模块加时间戳后的从B-CORE模块传递给A-CORE的数据；
所述GPS校时模块用来给从B-CORE模块传递给A-CORE的数据加时间戳；
通讯模块具有以太网RS-485协议，用以与首端上位机或末端上位机进行通讯；
所述传感器滤波电路具有硬件滤波功能，为高低频混合实时数据进行第一次硬件滤波；
AD转换模块将高低频混合实时数据模拟信号转换成数字信号，此高低频混合实时数据为首端超高频传感器和首端加速度传感器或者末端超高频传感器和末端加速度传感器传递来的。
所述AD转换模块与B-CORE模块通过SPI通讯协议相连接；所述GPS校时模块通过USART或者SPI通讯协议相连接；所述A-CORE模块和B-CORE模块之间通过DMA通讯协议相连接。
首端加速度传感器与末端加速度传感器之间的距离不小于25Uq，U为低频信号的周期，q为低频信号在管道中的传播速度。
一种基于高低频混合检测的管道泄漏监测系统的方法，包括以下步骤：
步骤1：从首端和末端上位机数据库中采集高频低频混合历史数据，该历史数据包括以下情况下几种工况类别：正常工况信号、泄漏信号和工况调整信号，从首端超高频传感器传递给首端上位机的历史高频数据记作X1，从首端加速度传感器传递给首端上位机的历史低频数据记作Y1，从末端超高频传感器传递给末端上位机的历史高频数据记作X2，从末端加速度传感器传递给末端上位机的历史低频数据记作Y2；
步骤2：分别对首端和末端上位机采集的高频低频混合历史数据进行预处理，包括步骤2.1～步骤2.3：
步骤2.1：将采集的混合信号进行滤波去噪处理：按小波包算法展开，通过一个低通滤波器h和一个高通滤波器g进行滤波，分解得到低频信号d2n和高频信号d2n+1各一组，每一次分解就将上层的第n个频带进一步分解为下层的第2n和第2n+1两个子频带；
将首端和末端上位机混合历史数据X1+Y1和X2+Y2，作为小波包的输入信号，记作信号序列xy(t)，小波包算法描述如下：假设共扼滤波器h(n)满足
其中，k，l为小波变换距离，n为小波分解层数，δkl为参数，Z为负数空间；
令：
g(k)＝(-1)kh(1-k) (3)
则小波包分解的递推公式为：
其中，d0(t)＝xy(t)，dn(t)为小波分解的系数，n为小波分解的层数；
步骤2.2：将步骤2.1分解后的信号d2n和d2n+1进行小波包重构方法处理：选择一最佳小波包分解基，确定分解的层次n，对待处理信号进行小波包分解，最佳小波包分解基对可定义为具有最大信息量或最小熵的小波包分解基，根据熵的标准计算最佳小波包分解基；
小包基熵的判断公式其中，s代表信号,si表示信号在正交小波包基上的投影系数，判断E(s)达到最小，其中，E(si)是求si的期望值；
小波包重构方法：由{d2n(t)}与{d2n+1(t)}求dn+1(t)
将小波重构后的信号dn+1(t)，重新记作x(t)；
步骤2.3：采用Hilbert变换的包络解调方法对小波重构后的信号x(t)进行包络解调得到包络信号zh(t)，其中，Hilbert变换公式为公式(6)
其中，z(t)为Hilbert变换后的信号，x(t)为输入的重构后的高频或低频信号j为负数表示，τ为任意参数表示；
Hilbert解调公式为，公式(7)
其中，zh(t)为Hilbert解调后的信号，
步骤3：将步骤2得到的Hilbert解调后包络信号zh(t)进行多信息域分析处理，特征提取，包括步骤3.1～步骤3.4，其中，步骤3.1～步骤3.3不分先后顺序：
步骤3.1：将步骤2得到的Hilbert解调后包络信号zh(t)进行时域下的特征提取，包括步骤3.1.1～3.1.2：
步骤3.1.1：提取步骤2得到的Hilbert解调后包络信号zh(t)中低频信号无量纲时域指标，无量纲时域指标包括：峰值指标C、波形指标S、脉冲指标I、裕度指标L和峭度指标K；
峰值指标C:zp为信号zh(t)的最大瞬时值，具体取信号的时间范围可以根据需要进行人工设定，例如设定1s内，这样zp即为1s内信号的最大瞬时值，zrms为信号均方根值；
波形指标S:zrms为信号zh(t)均方根值,zav信号的绝对均值，
脉冲指标I:zp为信号的最大瞬时值,zav信号的绝对均值，
裕度指标L:zp为信号的最大瞬时值,zr信号的方根幅值，
峭度指标K:z为信号zh(t)，μz为信号zh(t)的均值，σz为信号zh(t)的方差；
步骤3.1.2：提取步骤2得到的Hilbert解调后包络信号zh(t)中高频信号无量纲时域指标，无量纲时域指标于步骤3.1.1中指标一致；
步骤3.2：将步骤2得到的Hilbert解调后包络信号zh(t)进行频域下的特征提取，包括步骤3.2.1～3.2.2：
步骤3.2.1：将步骤2得到的Hilbert解调后包络信号zh(t)中低频信号进行傅里叶变换x(ω)，提取频域下的特征参量，频域特征参量包括：衰减系数β、功率谱密度Sx、输入输出的功率谱密度Sxy、均方频率MSF、频率方差VF；
衰减系数β：ω为频率，c为传播速度，具体为其中，a为管道半径，h为管壁厚度；
功率谱密度Sx：z为信号zh(t)傅里叶变换后的信号z(ω)；
输入输出的功率谱密度Sxy:H(ω)为系统的频率响应，H(ω)＝e-jwx/c c为传播速度，y为输出响应，
均方频率MSF：S(ω)中为信号的功率谱，ω为频率；
频率方差VF:
步骤3.2.2：将步骤2得到的Hilbert解调后包络信号zh(t)中高频信号进行傅里叶变换，提取频域下的特征参量，频域特征参量于步骤3.2.1中指代一致；
步骤3.3：将步骤2得到的Hilbert解调后包络信号zh(t)进行时频域下的特征提取，包括步骤3.3.1～步骤3.3.4：
步骤3.3.1：定义故障特征向量：根据式(8)计算各频带信号的能量
设信号经小波包分解后,第j层第k个频带的Sjk所对应的能量为Ejk,则有
其中，p为数据长度；j为小波包分解层次，j根据需要自己定；另M＝2j-1,则k＝0,1,2,3,4…M,为分解频带的序号；zkm为Sjk的离散点的幅值，为了公式表达方便zh(t)，记作zkm，总能量E为各子频带的能量之和,即
当故障发生时,各频带内信号的能量有明显的变化,因此,可用信号分解后的各频带能量占总能量的百分比作为故障特征向量，因此，利用小波包理论提取到的时频域特征向量可表达为：
Wjk＝[Ej0,Ej0,…,EjM]T/E (10)
步骤3.3.2：用步骤2.1及步骤2.2的小波分解与重构的方法，将能量Ejk进行3层分解，提取8个频率成分信号特征，以S30表示d30的重构信号,S31表示d31的重构信号,其他依此类推；利用第三层所有子带的对应信号来表示总信号S为：
S＝S30+S31+S32+S33+S34+S35+S36+S37 (11)
步骤3.3.3：计算各频带重构信号的能量，设E3k为S3k(k＝0,1,2,…,7)对应的能量则有：
其中，zkm(k＝0,1,2,…,7；m＝1,2,…,N)表示S3K的离散点赋值；
步骤3.3.4：构造故障特征向量，特征向量:
T＝[E30,E31,E32,E33,E34,E35,E36,E37]T (13)
当E3k(k＝0,1,2,…,7)中有一个能量比其他能量值大3倍以上，此时,数据分析会有所不便，有必要对T进行归一化处理,令:
此时,特征表示为
W3k＝E3k/E (15)
通过上述特征向量提取方法,可得到各状态的维特征向量,建立特征向量与系统状态的对应关系；
步骤3.4：取步骤3.1～步骤3.3得到时域下的特征、频域下的特征和时频域下的特征组成故障特征向量：峰值指标C、波形指标S、脉冲指标I、裕度指标L和峭度指标K、衰减系数β、功率谱密度Sx、输入输出的胡功率谱密度Sxy、均方频率MSF、频率方差VF、W30、W31、W32、W33、W34、W35、W36、W37，共18个故障特征向量，组成特征矩阵XL；
步骤4：分别对首端和末端上位机采集的高频低频混合历史数据经过步骤2～步骤3后，得到特征矩阵XL,训练COIR-LDA判别，即协方差核切片逆回归算法-线性判别分析，COIR算法是由KSIR算法改进得到的，分别c类工况下的数据特征，KSIR算法即核切片逆回归算法，具体步骤包括步骤4.1～步骤4.8：
步骤4.1：由核函数定义的非线性变换，核函数取高斯核函数其中，exp就是e指数函数，||·||表示绝对距离，x为输入，x指代特征矩阵XL，xi为核函数中心，σ为函数的宽度参数，将18维测量特征矩阵XL映射至高维特征空间F，得λ为特征值；
步骤4.2：Y切成c份，分别为I1,…,Ic；Y＝f(b1x,b2x,b3x…bnx)，计算切片内I1,…,Ic的均值M；
其中，表示n*c的0/1指示矩阵,每一行只有一个元素是1，表示混合历史数据本属于某一切片，其余为0，其中，c表示工况的数量，n是样本数量，b1,b2…,bn是投影向量b，δh是个函数；
步骤4.3：计算高维特征空间F中每段的均值向量和加权协方差阵即：
V＝MPMT (17)
我们把矩阵V的前d个特征向量记成为切片内样本占全部样本的比例，l＝1,…,d，n是样本数量；
步骤4.4：中心子空间的方向可以通过获得，x输入的特征,也就是XL,n样本数量；
步骤4.5:利用核函数的内积运算Kij(xi,xj)＝<φ(xi),φ(xj)>求解特征值分解Vv＝λv，可以得到下面对偶方程：
得到特征值λ，其中n是样本数量，h表示工况的数量，α＝[α1，...αn]T；
步骤4.6：投影向量b从特征向量v得到，为了求解令其中β＝[β1,…,βn]T，并且两侧同乘以即可以得到非线性子空间β的闭式解：
步骤4.7:非线性中心子空间即可以通过d个从构成的基函数表示，如果给定一个新的样本x*，x*指的是XL，它的低维表示z*∈Rd可以通过将φ(x*)投影中心子空间得到，也就是说，z*的l个元素是：其中，k*＝[kx(x1,x*),…,kx(xn,x*))]T，Rd是前d个特征向量组成的R维空间；
步骤4.8:将剩余低频历史数据及高频历史数据经过步骤2进行预处理，步骤3特征提取，并进行步骤4，最终得到2c种高低频历史数据的非线性中心子空间z*,每个z*有对应的
步骤5：改进步骤4所述的(18)式，得到COIR算法，即协方差核切片逆回归算法，具体步骤包括步骤5.1～步骤5.5：
步骤5.1：对任意的函数f∈Rx，如果存在一个函数g∈Rx，使得对所有的y，有E[f(x)|y]＝g(y)，其中，E是求期望，var是求方差，那么：
步骤5.2：使用方差公式可以得到：
var[E[φ(x)|y]]＝var(φ(x))-E[var[φ(x)|y]] (21)
步骤5.3：通过步骤5.1和步骤5.2，步骤5.3的估计量可以写成
其中，
步骤5.4：COIR算法需要对上式进行特征值分解，给定观测样本再令φy＝[φ(y1),…,φ(yn)]，那么步骤5.4的估计量可以写成
步骤5.5：同样，特征值分解var[E[φ(x)|y]]v＝λv可以写成下式：
其中，ε为人工设定的参数。
步骤6：实时采集管道高频及低频数据，通过步骤2进行预处理，以秒为单位，送入步骤3进行时域、频域和时频域特征提取，再进入步骤4，通过训练得到的COIR-LDA判别模型，根据得到的特征向量与建立的四种工况情况的数据库进行比对，根据LDA判别函数得到实时数据对应于某类工况，具体包括步骤6.1～步骤6.2；
步骤6.1：对于n个步骤2得到的小波包重构后的信号的c类问题,假设ωi(i＝1,2,…,c)类模式的先验概率为p(ωi),显然有
步骤6.2：将步骤5的COIR算法作用于采集的数据X,提取和输出F的中心子空间，步骤4求出来的投影向量b的空间表示的样本数据Z,然后构造第i类别LDA的判别函数:
式中zt∈Rr为F的中心子空间中带判别的分类属性矢量，zi则为第i类训练样本的在F的中心子空间的平均分类属性矢量，为训练样本在F的中心子空间的类内离散度矩阵，其中Ww,z是矩阵xw xz的矩阵，xz为第i类特征样本xw是输入的样本，得到c种工况下对应的Li(zt)，建立管道泄漏检测的数据库。
步骤7：实时采集管道高频及低频数据，通过步骤2进行预处理，然后送入步骤3进行时域、频域和时频域特征提取，再进入步骤4，通过训练得到的COIR-LDA判别模型，根据得到的特征向量与建立的四种工况情况的数据库进行比对，计算两者之间的相似程度；
相似程度的计算公式如下所示：
其中，Sk为相似度，C(x)new为新采集的数据，为第k种类型的测试样本，为k类样本的一个集合；
若Si≥0.85，则当前数据为第k类数据；否则，结合现场状况的实际分析及人为诊断确定故障的类型，并将其收入管道泄漏检测的数据库中。
步骤8：根据步骤7，若实时数据判定为大泄漏信号或小泄漏信号时，则通过高频信号和低频信号混合计算泄漏点位置，对管道泄漏进行定位，否则返回步骤6，继续实时监测；
管道泄露进行定位的具体过程包括步骤8.1～步骤8.5：
步骤:8.1：根据首端和末端传感器的时间差值计算泄漏点，如下式所示，
式中，K1为泄漏点到首端超高频传感器的距离；S1为管线首端和末端超高频传感器之间的距离；τ1为高频信号传到首端和末端超高频传感器的时差；o1为高频信号在管线内的传播速度；
步骤8.2：在泄漏点的上游站及下游站各选取泄漏发生点前30秒的高频数据均值和泄漏点后15秒的高频数据均值，计算上下游站高频信号变化比，分别如下两式所示，
其中：δ1、δ2分别为首端和末端的高频信号变化比；X1+、X2+分别为泄漏发生点前30秒首端和末端的高频数据均值；X1-、X2-分别为泄漏发生点后15秒首端和末端的高频数据均值；
则管道泄漏时总高频信号变化比δp为：其中，μ为由现场管道长度和正常输送管道的大小确定的参数；
步骤8.3：对低频信号的特征值在首端即末端的时间差值，计算泄漏点，如下式所示，
式中，K2为泄漏点到首端加速度传感器的距离；S2为管线首端和末端加速度传感器之间的距离；τ2为声波传到首端和末端加速度传感器的时差；o2为低频信号在管线内的传播速度；
步骤8.4：在泄漏点的首端及末端各选取泄漏发生点前5秒的低频数据均值和泄漏点后3秒的低频数据均值，计算首端和末端低频信号变化比，分别如下两式所示，
其中，δ3、δ4分别为首端和末端的低频信号变化比；Y1+、Y2+分别为泄漏发生点前5秒首端和末端的低频信号数据均值；Y1-、Y2-分别为泄漏发生点后3秒首端和末端的低频数据均值；
则管道泄漏时总低频信号变化比δs为：
步骤8.5：根据泄漏点到上游站超高频传感器的最终距离为：
其中，T为低频信号的周期。
有益技术效果：
本发明使用的下位机采集数据精度更高、更稳定以及成本更低；采用多次滤波包括对高频和低频信号的单独硬件及软件滤波，更好的屏蔽了噪声的干扰；在数据预处理阶段，通过小波包法对信号进行分解与滤波，并通过Hilbert变换的包络解调方法对滤波去噪后信号进行包络解调得到包络信号，确保了得到能反应故障特征的微弱信号并且还原了信号的准确性；采用COIR-LDA判别分析的方法综合分析了高频和低频信号，确保了泄漏的特征值的准确性，最大限度的减少了误报，对小流量的泄漏有了更准确的判断；COIR算法与KSIR算法相比，COIR算法不需要对响应变量Y切片，可以处理多元变量的Y，并且COIR算法可以轻易地推广到半监督情况；综合利用高频和低频信号判断管道泄漏的类型，针对负压波法对小泄漏量和声波法对大泄漏量的判断不精确的问题，很好地解决了小流量泄漏精准判断问题；给出了判断泄漏后确定具体泄漏位置的方法，为排除故障提供可靠的数据依据。
附图说明
图1本发明实施例的基于高低频混合检测的管道泄漏监测系统；
图2本发明实施例的首端下位机或末端下位所包含模块示意图；
图3本发明实施例的高低频混合检测的管道泄漏监测系统整体流程图；
图4本发明实施例的高低频混合检测的管道泄漏监测系统预处理流程图；
图5本发明实施例的高低频混合检测的管道泄漏监测系统特征提取流程图；
图6本发明实施例的高低频混合检测的管道泄漏监测系统COIR处理流程图；
图7本发明实施例的高频大泄漏特征矩阵截图；
图8本发明实施例的低频正常工况特征矩阵截图；
图9本发明实施例的8种高低频历史数据计算得到的投影向量；
图10本发明实施例的线性判别分析分类结果；
图中：1-首端超高频传感器；2-首端加速度传感器；3-首端下位机；4-首端上位机；5-末端超高频传感器；6-末端加速度传感器；7-末端下位机；8-末端上位机；
具体实施方式
下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。
一种基于高低频混合检测的管道泄漏监测系统，包括首端超高频传感器1、首端加速度传感器2、首端下位机3、首端上位机4、末端超高频传感器5、末端加速度传感器6、末端下位机7和末端上位机8，如图1所示；
所述首端超高频传感器1和首端加速度传感器2分别安装在待测试管道的首端，首端超高频传感器1和首端加速度传感器2通过电线与首端下位机3相连接，首端下位机3通过网线与首端上位机4相连接；
所述末端超高频传感器5和末端加速度传感器6分别安装在待测试管道的末端，末端超高频传感器6和末端加速度传感器6通过电线与末端下位机7相连接，末端下位机7通过网线与末端上位机8相连接；
所述首端超高频传感器1和首端加速度传感器2，与待测试管道中输送介质相接触，首端超高频传感器1用于采集首端管道高频信号的实时数据，首端加速度传感器2用于采集低频信号的实时数据；
本实施例中，首端超高频传感器1和末端超高频传感器5均采用MPT202HT超高频传感器，其主要参数如下：
压力变送器的测量压力范围为0～6.6MPa；信号分辨率0.02％，准确度±0.08％，更新速率1kHz；输出信号为4～20mADC(二线制)，带负载能力不小于1000Ω，供电电源为24VDC(9-36VDC)；具有承受最大量程的2倍的过载能力；环境温度每变化50°F(28℃)的影响优于：±(0.02％量程上限+0.15％量程)；静压每变化1000psi(6.5MPa)的影响优于：±0.1％量程上限。
本实施例中，首端加速度传感器2和末端加速度传感器6采用AC192系列，具有以下的技术指标：
电压灵敏度：100mv/g；测量频率范围：0.2～2000Hz；量程：5g；线性度：≤1％；工作温度：-20～120℃；抗冲击：10g；输出方式：顶端M5；激励电流：2～10mA，激励电压：10-24VDC。
所述首端下位机3用于存储从首端超高频传感器1和首端加速度传感器2传递出的超高频信号和低频信号的混合实时数据，将该高低频混合实时数据进行模数转换且加上时间戳，传递给首端上位机4；
所述首端上位机4用于接收并保存由首端下位机3传递出来的加上时间戳的高低频混合实时数据，并将该高低频混合实时数据保存到首端上位机4数据库中，首端上位机4通过编程来实现管道泄漏检测功能；
所述末端超高频传感器5和末端加速度传感器6，与待测试管道中输送介质相接触，末端超高频传感器5用于采集末端管道高频信号的实时数据，末端加速度传感器6用于采用低频信号的实时数据；
所述末端下位机7用于存储从末端超高频传感器5和末端加速度传感器6传递出的道高频信号和低频信号的高低频混合实时数据，将该高低频混合实时数据进行模数转换且加上时间数据，传递给末端上位机8；
所述末端上位机8用于接收并保存由末端下位机7传递出来的加上时间戳的高低频混合实时数据，并将该高低频混合实时数据保存到首端上位机8数据库中，末端上位机8通过编程来实现管道泄漏检测功能；
所述首端下位机3和末端下位机7均包括：A-CORE模块、B-CORE模块、电源模块、存储模块、GPS校时模块、通讯模块、调试模块、传感器滤波电路、AD转换模块，如图2所示；
所述电源模块、存储模块、GPS校时模块、通讯模块、调试模块分别与A-CORE模块相连接，其中，存储模块通过并行总线与A-CORE模块相连接，GPS校时模块通过USART或SPI通讯协议与A-CORE模块相连接，通讯模块通过以太网RS-485协议与A-CORE模块相连接；传感器滤波电路与AD转换模块相连接，AD转换模块与B-CORE模块相连接；A-CORE模块通过网线与B-CORE模块相连接；电源模块与A-CORE模块和B-CORE模块相连接；
所述A-CORE模块是主控模块，将B-CORE模块传递过来的数据进行存储，并加以时间戳转发给首端或者末端上位机；本实例A-CORE模块具体型号为：STM32103ZE，B-CORE模块具体型号为：XC9536；
所述B-CORE模块是A-CORE模块的从设备，接收超高频传感器和加速度传感器传递出的高频信号和低频信号的高低频混合实时数据，并通过传感器滤波电路完成对该实时数据的硬件滤波，并且通过AD转换模块完成模数信号的转换，然后将模数信号转换后的数据传递给A-CORE模块；
所述电源模块为A-CORE模块和B-CORE模块提供电源，具体型号为220V转24V明纬型电源S-50-24以及WRB2405S-3WR2DC-DC电源芯片；
所述存储模块保存用GPS校时模块加时间戳后的从B-CORE模块传递给A-CORE的数据，具体为闪迪sd卡32g内存卡；
所述GPS校时模块用来给从B-CORE模块传递给A-CORE的数据加时间戳，具体为GPS北斗双模定位模块ATK-S1216；
通讯模块具有以太网RS-485协议，用以与首端上位机4或末端上位机8进行通讯，具体为W5500网络芯片；
所述传感器滤波电路具有硬件滤波功能，为高低频混合实时数据进行第一次硬件滤波，具体为MAX281AEPA直插DIP8有源滤波芯片；
AD转换模块将高低频混合实时数据模拟信号转换成数字信号，此高低频混合实时数据为首端超高频传感器1和首端加速度传感器2或者末端超高频传感器5和末端加速度传感器6传递来的，具体型号为Ad7606芯片；
其中，所述AD转换模块与B-CORE模块通过SPI通讯协议相连接；
所述GPS校时模块通过USART或者SPI通讯协议相连接；
所述A-CORE模块和B-CORE模块之间通过DMA通讯协议相连接；
首端加速度传感器2与末端加速度传感器6之间的距离不小于25Uq，U为低频信号的周期，本实施例中U为1秒，q为低频信号在管道中的传播速度，本实施例中q＝430m/s；
本实施例中首端下位机和末端下位机的技术指标的技术指标要求如下：
采样频率范围：0-2.5KHz；驱动信号：3.3V及5V；采样信号的通道数：至少4路；
采样精度：≥0.001；通信方式：以太网通信；供电方式：24V直流供电。
本实施例中，对低频信号采样频率选取为500Hz，对高频信号的采样频率甚至为500Hz。
本发明还提供一种基于高低频混合检测的管道泄漏监测系统的方法，采用上述的检测系统实现，该方法包括以下步骤，如图3所示：
步骤1：从首端和末端上位机数据库中采集高频低频混合历史数据，该历史数据包括四种情况下历史数据：正常工况信号、大泄漏信号、小泄露信号和工况调整信号，从首端超高频传感器传递给首端上位机的历史高频数据记作X1，从首端加速度传感器传递给首端上位机的历史低频数据记作Y1，从末端超高频传感器传递给末端上位机的历史高频数据记作X2，从末端加速度传感器传递给末端上位机的历史低频数据记作Y2，所有数据时间长度为1秒钟；
步骤2：分别对首端和末端上位机采集的高频低频混合历史数据进行预处理，包括步骤2.1～步骤2.3，如图4所示：
步骤2.1：将采集的混合信号按小波包算法展开，通过一个低通滤波器h和一个高通滤波器g进行滤波，分解得到低频信号d2n和高频信号d2n+1各一组，每一次分解就将上层的第n个频带进一步分解为下层的第2n和第2n+1两个子频带；
将首端和末端上位机混合历史数据X1+Y1和X2+Y2，作为小波包的输入信号，记作信号序列xy(t)，小波包算法描述如下：假设共扼滤波器h(n)满足
其中，k，l为小波变换距离，n为小波分解层数，δkl为参数，Z为负数空间；
令：
g(k)＝(-1)kh(1-k) (3)
则小波包分解的递推公式为：
其中，d0(t)＝xy(t)，dn(t)为小波分解的系数，n为小波分解的层数。
步骤2.2：将步骤2.1分解后的信号d2n和d2n+1进行小波包重构方法处理：选择一最佳小波包分解基，确定分解的层次n，对待处理信号进行小波包分解。最佳小波包分解基对可定义为具有最大信息量或最小熵的小波包分解基，根据熵的标准计算最佳小波包分解基；
小包基熵的判断公式其中，s代表信号,si表示信号在正交小波包基上的投影系数，判断E(s)达到最小，其中，E(si)是求si的期望值；
小波包重构方法：由{d2n(t)}与{d2n+1(t)}求dn+1(t)
将小波重构后的信号dn+1(t)，重新记作x(t)；
步骤2.3：采用Hilbert变换的包络解调方法对小波重构后的信号x(t)进行包络解调得到包络信号zh(t)，其中，Hilbert变换公式为公式(6)
其中，z(t)为Hilbert变换后的信号，x(t)为输入的重构后的高频或低频信号j为负数表示，τ为任意参数表示；
Hilbert解调公式为，公式(7)
其中，zh(t)为Hilbert解调后的信号，
步骤3：将步骤2得到的Hilbert解调后包络信号zh(t)进行多信息域分析处理，包括步骤3.1～步骤3.4，如图5所示：
步骤3.1：将步骤2得到的Hilbert解调后包络信号zh(t)进行时域下的特征提取，包括步骤3.1.1～3.1.2：
步骤3.1.1：提取步骤2得到的Hilbert解调后包络信号zh(t)中低频信号无量纲时域指标，无量纲时域指标包括：峰值指标C、波形指标S、脉冲指标I、裕度指标L和峭度指标K；
峰值指标C:zp为信号zh(t)的最大瞬时值，具体取信号的时间范围可以根据需要进行人工设定，例如设定1s内，这样zp即为1s内信号的最大瞬时值，zrms为信号均方根值；
波形指标S:zrms为信号zh(t)均方根值,zav信号的绝对均值，
脉冲指标I:zp为信号的最大瞬时值,zav信号的绝对均值，
裕度指标L:zp为信号的最大瞬时值,zr信号的方根幅值，
峭度指标K:z为信号zh(t)，μx为信号zh(t)的均值，σx为信号zh(t)的方差；
步骤3.1.2：提取步骤2得到的Hilbert解调后包络信号zh(t)中高频信号无量纲时域指标，无量纲时域指标于步骤3.1.1中指标一致；
步骤3.2：将步骤2得到的Hilbert解调后包络信号zh(t)进行频域下的特征提取，包括步骤3.2.1～3.2.2：
步骤3.2.1：将步骤2得到的Hilbert解调后包络信号zh(t)中低频信号进行傅里叶变换x(ω)，提取频域下的特征参量，频域特征参量包括：衰减系数β、功率谱密度Sx、输入输出的功率谱密度Sxy、均方频率MSF、频率方差VF；
衰减系数β：ω为频率，c为传播速度，具体为其中，a为管道半径，h为管壁厚度；
功率谱密度Sx：z为信号zh(t)傅里叶变换后的信号z(ω)；
输入输出的功率谱密度Sxy:H(ω)为系统的频率响应，H(ω)＝e-jwx/c c为传播速度，y为输出响应，
均方频率MSF：S(ω)中为信号的功率谱，ω为频率；
频率方差VF:
步骤3.2.2：将步骤2得到的Hilbert解调后包络信号zh(t)中高频信号进行傅里叶变换，提取频域下的特征参量，频域特征参量于步骤3.2.1中指代一致；
步骤3.3：将步骤2得到的Hilbert解调后包络信号zh(t)进行时频域下的特征提取，包括步骤3.3.1～步骤3.3.4：
步骤3.3.1：定义故障特征向量：根据式(8)计算各频带信号的能量
设信号经小波包分解后,第j层第k个频带的Sjk所对应的能量为Ejk,则有
其中，p为数据长度；j为小波包分解层次，j根据需要自己定；另M＝2j-1,则k＝0,1,2,3,4…M,为分解频带的序号；zkm为Sjk的离散点的幅值，为了公式表达方便zh(t)，记作zkm。总能量E为各子频带的能量之和,即
当故障发生时,各频带内信号的能量有明显的变化,因此,可用信号分解后的各频带能量占总能量的百分比作为故障特征向量。因此,利用小波包理论提取到的时频域特征向量可表达为：
Wjk＝[Ej0,Ej0,…,EjM]T/E (10)
步骤3.3.2：用步骤2.1及步骤2.2的小波分解与重构的方法，将能量Ejk进行3层分解，提取8个频率成分信号特征，以S30表示d30的重构信号,S31表示d31的重构信号,其他依此类推。利用第三层所有子带的对应信号来表示总信号S为：
S＝S30+S31+S32+S33+S34+S35+S36+S37 (11)
步骤3.3.3：计算各频带重构信号的能量。设E3k为S3k(k＝0,1,2,…,7)对应的能量则有：
其中，zkm(k＝0,1,2,…,7；m＝1,2,…,N)表示S3K的离散点赋值。
步骤3.3.4：构造故障特征向量。特征向量:
T＝[E30,E31,E32,E33,E34,E35,E36,E37]T (13)
当E3k(k＝0,1,2,…,7)中有一个能量比其他能量值大3倍以上，此时,数据分析会有所不便，有必要对T进行归一化处理,令:
此时,特征表示为
W3k＝E3k/E (15)
通过上述特征向量提取方法,可得到各状态的维特征向量,建立特征向量与系统状态的对应关系。
步骤3.4：取步骤3.1～步骤3.3得到时域下的特征、频域下的特征和时频域下的特征组成故障特征向量：峰值指标C、波形指标S、脉冲指标I、裕度指标L和峭度指标K、衰减系数β、功率谱密度Sx、输入输出的胡功率谱密度Sxy、均方频率MSF、频率方差VF、W30、W31、W32、W33、W34、W35、W36、W37，共18个故障特征向量，组成特征矩阵XL；
按照以上方法，得到高频大泄漏特征矩阵，高频小泄漏特征矩阵，高频正常工况特征矩阵，高频工况调整特征矩阵，低频大泄漏特征矩阵，低频小泄漏特征矩阵，低频正常工况特征矩阵，低频工况调整特征矩阵，其中，高频大泄漏特征矩阵如图7所示，低频正常工况特征矩阵如图8所示；
步骤4：分别对首端和末端上位机采集的高频低频混合历史数据经过步骤2～步骤3后，得到特征矩阵XL,训练COIR-LDA判别，即协方差核切片逆回归算法-线性判别分析，COIR算法是由KSIR算法改进得到的，分别c类工况下的数据特征，KSIR算法，即核切片逆回归算法具体步骤如下，如图6所示：
步骤4.1：由核函数定义的非线性变换，核函数取高斯核函数其中，exp就是e指数函数，||·||表示绝对距离，x为输入，x指代特征矩阵XL，xi为核函数中心，σ为函数的宽度参数，将18维测量特征矩阵XL映射至高维特征空间F，得λ为特征值；
步骤4.2：Y切成4份，分别为I1,…,I4；Y＝f(b1x,b2x,b3x…bnx)，计算切片内I1,…,I4的均值M；
其中，表示n*c的0/1指示矩阵，每一行只有一个元素是1，表示混合历史数据属于某一切片，其余为0，其中，c表示工况的数量c＝4，n是样本数量，b1,b2…,bn是投影向量b，δh是个函数；
步骤4.3：计算高维特征空间F中每段的均值向量和加权协方差阵即：
V＝MPMT (17)
我们把矩阵V的前d个特征向量记成为切片内样本占全部样本的比例，l＝1,…,d，n是样本数量；
步骤4.4：中心子空间的方向可以通过获得，x输入的特征,也就是XL,n样本数量；
步骤4.5:利用核函数的内积运算Kij(xi,xj)＝<φ(xi),φ(xj)>求解特征值分解Vv＝λv，可以得到下面对偶方程：
得到特征值λ，其中n是样本数量，h表示工况的数量，α＝[α1，...αn]T；
步骤4.6：投影向量b从特征向量v得到，为了求解令其中β＝[β1,…,βn]T，并且两侧同乘以即可以得到非线性子空间β的闭式解：
步骤4.7:非线性中心子空间即可以通过d个从构成的基函数表示。如果给定一个新的样本x*，x*指的是XL，它的低维表示z*∈Rd可以通过将φ(x*)投影中心子空间得到，也就是说，z*的l个元素是：其中，k*＝[kx(x1,x*),…,kx(xn,x*))]T，Rd是前d个特征向量组成的R维空间；
步骤4.8:将剩余低频历史数据及高频历史数据经过步骤2进行预处理，步骤3特征提取，并进行步骤4，最终得到8种高低频历史数据的非线性中心子空间z*,每个z*有对应的最后8种高低频历史数据计算得到的投影向量b如图9所示；
步骤5：改进步骤4所述的(18)式，得到COIR算法，即协方差核切片逆回归算法，具体步骤如下：
步骤5.1：对任意的函数f∈Rx，如果存在一个函数g∈Rx，使得对所有的y，有E[f(x)|y]＝g(y)，其中，E是求期望，var是求方差，那么：
步骤5.2：使用方差公式可以得到：
var[E[φ(x)|y]]＝var(φ(x))-E[var[φ(x)|y]] (21)
步骤5.3：通过步骤5.1和步骤5.2，步骤5.3的估计量可以写成
其中，
步骤5.4：COIR算法需要对上式进行特征值分解，给定观测样本再令φy＝[φ(y1),…,φ(yn)]，那么步骤5.4的估计量可以写成
步骤5.5：同样，特征值分解var[E[φ(x)|y]]v＝λv可以写成下式：
其中，ε为人工设定的参数；
步骤6：实时采集管道高频及低频数据，通过步骤2进行预处理，以秒为单位，送入步骤3进行时域、频域和时频域特征提取，再进入步骤4，通过训练得到的COIR-LDA判别模型，根据得到的特征向量与建立的四种工况情况的数据库进行比对，根据LDA判别函数得到实时数据对应于某类工况；
步骤6.1：对于n个步骤2得到的小波包重构后的信号的c类问题,假设ωi(i＝1,2,…,c)类模式的先验概率为p(ωi),显然有
步骤6.2：将步骤5的COIR算法作用于采集的数据X,提取和输出F的中心子空间，步骤4求出来的投影向量b的空间表示的样本数据Z，然后构造第i类别LDA的判别函数:
式中zt∈Rr为F的中心子空间中带判别的分类属性矢量，则为第i类训练样本的在F的中心子空间的平均分类属性矢量，为训练样本在F的中心子空间的类内离散度矩阵。其中Ww,z是矩阵xw xz的矩阵，xz为第i类特征样本xw是输入的样本，得到c种工况下对应的Li(zt)，建立管道泄漏检测的数据库；
LDA线性判别分析分类结果，如图10所示，其中，正方形代表大泄漏数据分类结果，星型代表小泄漏数据分类结果，圆形代表工况调整数据分类结果，三角形是正常情况数据分类结果；
步骤7：实时采集管道高频及低频数据，通过步骤2进行预处理，以秒为单位，送入步骤3进行时域、频域和时频域特征提取，再进入步骤4，通过训练得到的COIR-LDA判别模型，根据得到的特征向量与建立的四种工况情况的数据库进行比对，计算两者之间的相似程度；
相似程度的计算公式如下所示：
其中，Sk为相似度，C(x)new为新采集的数据，为第k种类型的测试样本，为k类样本的一个集合；
若Si≥0.85，则当前数据为第k类数据；否则，结合现场状况的实际分析及人为诊断确定故障的类型，并将其收入管道泄漏检测的数据库中。
步骤8：根据步骤:7，若实时数据判定为大泄漏信号或小泄漏信号时，则通过高频信号和低频信号混合计算泄漏点位置，对管道泄漏进行定位，否则返回步骤6，继续实时监测。
本实施例中当发生管道泄露时，进行定位的具体过程及计算结果为：
步骤8.1：根据首端和末端传感器的时间差值计算泄漏点，如下式所示，
式中，K1为泄漏点到首端超高频传感器的距离；S1为管线首端和末端超高频传感器之间的距离；τ1为高频信号传到首端和末端超高频传感器的时差；o1为高频信号在管线内的传播速度；
其中，S1为30km,τ1为2s,o1为1430m/s K1＝13570；
步骤8.2：在泄漏点的上游站及下游站各选取泄漏发生点前30秒的高频数据均值和泄漏点后15秒的高频数据均值，计算上下游站高频信号变化比，分别如下两式所示，
其中：δ1、δ2分别为首端和末端的高频信号变化比；X1+、X2+分别为泄漏发生点前30秒首端和末端的高频数据均值；X1-、X2-分别为泄漏发生点后15秒首端和末端的高频数据均值；
则管道泄漏时总高频信号变化比δp为：其中，μ为由现场管道长度和正常输送管道的大小确定的参数；
其中，X1+为5.2789Mpa，X2+为5.2808Mpa，X1-为1.5248Mpa，X2-为1.5215Mpa，
δ1＝3.4620，δ2＝3.4707；
则管道泄漏时总高频信号变化比δp为：其中，μ为由现场管道长度和正常输送管道的大小确定的参数；μ＝1，δp＝3.46635；
步骤8.3：对低频信号的特征值在首端即末端的时间差值，计算泄漏点。如下式所示，
式中，K2为泄漏点到首端加速度传感器的距离；S2为管线首端和末端加速度传感器之间的距离；τ2为声波传到首端和末端加速度传感器的时差；o2为低频信号在管线内的传播速度；
其中，S2为30km,τ2为2s,o2为430m/s，K1＝14570；
步骤8.4：在泄漏点的首端及末端各选取泄漏发生点前5秒的低频数据均值和泄漏点后3秒的低频数据均值，计算首端和末端低频信号变化比，分别如下两式所示，
其中，δ3、δ4分别为首端和末端的低频信号变化比；Y1+、Y2+分别为泄漏发生点前5秒首端和末端的低频信号数据均值；Y1-、Y2-分别为泄漏发生点后3秒首端和末端的低频数据均值；
则管道泄漏时总低频信号变化比δs为：
Y1+＝4.2657，Y2+＝4.2676，Y1-＝0.8849，Y2-＝0.8842；
δ3＝4.8205，δ4＝4.8265；
则管道泄漏时总低频信号变化比δs为δs＝4.8235；
步骤8.5：根据泄漏点到上游站超高频传感器的最终距离为：
其中，T为低频信号的周期。
其中，T为1，最后计算结果为：K＝16654m；
影响COIR-LDA模型分类正确率的因素有:F的EDR空间维数r、核函数的宽度系数σ以及各类别的先验概率P(ωi)等。由于本研究实例未提供各类别的先验概率信息,故一般简单取值为各类别样本个体数占训练样本容量的比例即P(ωi)＝Ph/I(h＝1,2,3，…，c)。F的中心子空间维数R,一般取值为R＝c-1,本专利研究的是4分类问题,即r＝3。核函数的宽度系数σ,可由经验试探法通过分类器对样本的最低分类错误个数确定。σ的不同取值与误判个数ne的关系为随着σ的减小,ne值呈现先下降、后又上升之趋势,在σ＝0.02,ne取得最小值,由此确定σ的优化取值。
具体应用结果表明：
采用COIR-LDA判别分析的方法综合分析了高频和低频信号，确保了泄漏的特征值的准确性，最大限度的减少了误报，对小流量的泄漏有了更准确的判断；并给出了判断泄漏后确定具体泄漏位置的方法，为排除故障提供可靠的数据依据。

Claims (4)

1.一种基于高低频混合检测的管道泄漏监测系统，其特征在于，包括首端超高频传感器、首端加速度传感器、首端下位机、首端上位机、末端超高频传感器、末端加速度传感器、末端下位机和末端上位机；
所述首端超高频传感器和首端加速度传感器分别安装在待测试管道的首端，首端超高频传感器和首端加速度传感器分别通过电线与首端下位机相连接，首端下位机通过网线与首端上位机相连接；
所述末端超高频传感器和末端加速度传感器分别安装在待测试管道的末端，末端超高频传感器和末端加速度传感器分别通过电线与末端下位机相连接，末端下位机通过网线与末端上位机相连接；
所述首端超高频传感器和首端加速度传感器，与待测试管道中输送介质相接触，首端超高频传感器用于采集首端管道高频信号的实时数据，首端加速度传感器用于采集低频信号的实时数据；
所述首端下位机用于存储从首端超高频传感器和首端加速度传感器传递出的超高频信号和低频信号的混合实时数据，将该高低频混合实时数据进行模数转换且加上时间戳，传递给首端上位机；
所述首端上位机用于接收并保存由首端下位机传递出来的加上时间戳的高低频混合实时数据，并将该高低频混合实时数据保存到首端上位机数据库中，首端上位机通过编程来实现管道泄漏检测功能；
所述末端超高频传感器和末端加速度传感器，与待测试管道中输送介质相接触，末端超高频传感器用于采集末端管道高频信号的实时数据，末端加速度传感器用于采集低频信号的实时数据；
所述末端下位机用于存储从末端超高频传感器和末端加速度传感器传递出的道高频信号和低频信号的高低频混合实时数据，将该高低频混合实时数据进行模数转换且加上时间数据，传递给末端上位机；
所述末端上位机用于接收并保存由末端下位机传递出来的加上时间戳的高低频混合实时数据，并将该高低频混合实时数据保存到首端上位机数据库中，末端上位机通过编程来实现管道泄漏检测功能；
所述首端下位机和末端下位机均包括：A-CORE模块、B-CORE模块、电源模块、存储模块、GPS校时模块、通讯模块、调试模块、传感器滤波电路、AD转换模块；
所述电源模块、存储模块、GPS校时模块、通讯模块、调试模块分别与A-CORE模块相连接，其中，存储模块通过并行总线与A-CORE模块相连接，GPS校时模块通过USART或SPI通讯协议与A-CORE模块相连接，通讯模块通过以太网RS-485协议与A-CORE模块相连接；传感器滤波电路与AD转换模块相连接，AD转换模块与B-CORE模块相连接；A-CORE模块通过网线与B-CORE模块相连接；电源模块与A-CORE模块和B-CORE模块相连接；
所述A-CORE模块是主控模块，将B-CORE模块传递过来的数据进行存储，并加以时间戳转发给首端或者末端上位机；
所述B-CORE模块是A-CORE模块的从设备，接收超高频传感器和加速度传感器传递出的高频信号和低频信号的高低频混合实时数据，并通过传感器滤波电路完成对该实时数据的硬件滤波，并且通过AD转换模块完成模数信号的转换，然后将模数信号转换后的数据传递给A-CORE模块；
所述电源模块为A-CORE模块和B-CORE模块提供电源；
所述存储模块保存用GPS校时模块加时间戳后的从B-CORE模块传递给A-CORE的数据；
所述GPS校时模块用来给从B-CORE模块传递给A-CORE的数据加时间戳；
通讯模块具有以太网RS-485协议，用以与首端上位机或末端上位机进行通讯；
所述传感器滤波电路具有硬件滤波功能，为高低频混合实时数据进行第一次硬件滤波；
AD转换模块将高低频混合实时数据模拟信号转换成数字信号，此高低频混合实时数据为首端超高频传感器和首端加速度传感器或者末端超高频传感器和末端加速度传感器传递来的。
2.根据权利要求1所述的一种基于高低频混合检测的管道泄漏监测系统，其特征在于：
所述AD转换模块与B-CORE模块通过SPI通讯协议相连接；所述GPS校时模块通过USART或者SPI通讯协议相连接；所述A-CORE模块和B-CORE模块之间通过DMA通讯协议相连接。
3.根据权利要求1所述的一种基于高低频混合检测的管道泄漏监测系统，其特征在于：
首端加速度传感器与末端加速度传感器之间的距离不小于25Uq，U为低频信号的周期，q为低频信号在管道中的传播速度。
4.一种基于高低频混合检测的管道泄漏监测系统的方法，采用权利要求1所述的一种基于高低频混合检测的管道泄漏监测系统，其特征在于，包括以下步骤：
步骤1：从首端和末端上位机数据库中采集高频低频混合历史数据，该历史数据包括以下情况下几种工况类别：正常工况信号、泄漏信号和工况调整信号，从首端超高频传感器传递给首端上位机的历史高频数据记作X1，从首端加速度传感器传递给首端上位机的历史低频数据记作Y1，从末端超高频传感器传递给末端上位机的历史高频数据记作X2，从末端加速度传感器传递给末端上位机的历史低频数据记作Y2；
步骤2：分别对首端和末端上位机采集的高频低频混合历史数据进行预处理，包括步骤2.1～步骤2.3：
步骤2.1：将采集的混合信号进行滤波去噪处理：按小波包算法展开，通过一个低通滤波器h和一个高通滤波器g进行滤波，分解得到低频信号d2n和高频信号d2n+1各一组，每一次分解就将上层的第n个频带进一步分解为下层的第2n和第2n+1两个子频带；
将首端和末端上位机混合历史数据X1+Y1和X2+Y2，作为小波包的输入信号，记作信号序列xy(t)，小波包算法描述如下：假设共扼滤波器h(n)满足：
其中，k，l为小波变换距离，n为小波分解层数，δkl为参数，Z为负数空间；
令：
g(k)＝(-1)kh(1-k) (3)
则小波包分解的递推公式为：
其中，d0(t)＝xy(t)，dn(t)为小波分解的系数，n为小波分解的层数；
步骤2.2：将步骤2.1分解后的信号d2n和d2n+1进行小波包重构方法处理：选择一最佳小波包分解基，确定分解的层次n，对待处理信号进行小波包分解，最佳小波包分解基对可定义为具有最大信息量或最小熵的小波包分解基，根据熵的标准计算最佳小波包分解基；
小包基熵的判断公式其中，s代表信号,si表示信号在正交小波包基上的投影系数，判断E(s)达到最小，其中，E(si)是求si的期望值；
小波包重构方法：由{d2n(t)}与{d2n+1(t)}求dn+1(t)
将小波重构后的信号dn+1(t)，重新记作x(t)；
步骤2.3：采用Hilbert变换的包络解调方法对小波重构后的信号x(t)进行包络解调得到包络信号zh(t)，其中，Hilbert变换公式为公式(6)：
其中，z(t)为Hilbert变换后的信号，x(t)为输入的重构后的高频或低频信号j为负数表示，τ为任意参数表示；
Hilbert解调公式为，公式(7)：
其中，zh(t)为Hilbert解调后的信号，
步骤3：将步骤2得到的Hilbert解调后包络信号zh(t)进行多信息域分析处理，特征提取，包括步骤3.1～步骤3.4，其中，步骤3.1～步骤3.3不分先后顺序：
步骤3.1：将步骤2得到的Hilbert解调后包络信号zh(t)进行时域下的特征提取，包括步骤3.1.1～3.1.2：
步骤3.1.1：提取步骤2得到的Hilbert解调后包络信号zh(t)中低频信号无量纲时域指标，无量纲时域指标包括：峰值指标C、波形指标S、脉冲指标I、裕度指标L和峭度指标K；
峰值指标C:zp为信号zh(t)的最大瞬时值，具体取信号的时间范围可以根据需要进行人工设定，例如设定1s内，这样zp即为1s内信号的最大瞬时值，zrms为信号均方根值；
波形指标S:zrms为信号zh(t)均方根值,zav信号的绝对均值，脉冲指标I:zp为信号的最大瞬时值,zav信号的绝对均值，
裕度指标L:zp为信号的最大瞬时值,zr信号的方根幅值，
峭度指标K:z为信号zh(t)，μz为信号zh(t)的均值，σz为信号zh(t)的方差；
步骤3.1.2：提取步骤2得到的Hilbert解调后包络信号zh(t)中高频信号无量纲时域指标，无量纲时域指标于步骤3.1.1中指标一致；
步骤3.2：将步骤2得到的Hilbert解调后包络信号zh(t)进行频域下的特征提取，包括步骤3.2.1～3.2.2：
步骤3.2.1：将步骤2得到的Hilbert解调后包络信号zh(t)中低频信号进行傅里叶变换x(ω)，提取频域下的特征参量，频域特征参量包括：衰减系数β、功率谱密度Sx、输入输出的功率谱密度Sxy、均方频率MSF、频率方差VF；
衰减系数β：ω为频率，c为传播速度，具体为其中，a为管道半径，h为管壁厚度；
功率谱密度Sx：，z为信号zh(t)傅里叶变换后的信号z(ω)；
输入输出的功率谱密度Sxy:，H(ω)为系统的频率响应，H(ω)＝e-jwx/c ，c为传播速度，y为输出响应，
均方频率MSF：S(ω)中为信号的功率谱，ω为频率；
频率方差VF:
步骤3.2.2：将步骤2得到的Hilbert解调后包络信号zh(t)中高频信号进行傅里叶变换，提取频域下的特征参量，频域特征参量于步骤3.2.1中指代一致；
步骤3.3：将步骤2得到的Hilbert解调后包络信号zh(t)进行时频域下的特征提取，包括步骤3.3.1～步骤3.3.4：
步骤3.3.1：定义故障特征向量：根据式(8)计算各频带信号的能量，设信号经小波包分解后,第j层第k个频带的Sjk所对应的能量为Ejk,则有：
其中，p为数据长度；j为小波包分解层次，j根据需要自己定；另M＝2j-1,则k＝0,1,2,3,4…M,为分解频带的序号；zkm为Sjk的离散点的幅值，为了公式表达方便zh(t)，记作zkm，总能量E为各子频带的能量之和,即：
当故障发生时,各频带内信号的能量有明显的变化,因此,可用信号分解后的各频带能量占总能量的百分比作为故障特征向量，因此，利用小波包理论提取到的时频域特征向量可表达为：
Wjk＝[Ej0,Ej0,…,EjM]T/E (10)
步骤3.3.2：用步骤2.1及步骤2.2的小波分解与重构的方法，将能量Ejk进行3层分解，提取8个频率成分信号特征，以S30表示d30的重构信号,S31表示d31的重构信号,其他依此类推；利用第三层所有子带的对应信号来表示总信号S为：
S＝S30+S31+S32+S33+S34+S35+S36+S37 (11)
步骤3.3.3：计算各频带重构信号的能量，设E3k为S3k(k＝0,1,2,···,7)对应的能量则有：
其中，zkm(k＝0,1,2,···,7；m＝1,2,···,N)表示S3K的离散点赋值；
步骤3.3.4：构造故障特征向量，特征向量：
T＝[E30,E31,E32,E33,E34,E35,E36,E37]T (13)
当E3k(k＝0,1,2,···,7)中有一个能量比其他能量值大3倍以上，此时,数据分析会有所不便，有必要对T进行归一化处理,令：
此时,特征表示为：
W3k＝E3k/E (15)
通过上述特征向量提取方法,可得到各状态的维特征向量,建立特征向量与系统状态的对应关系；
步骤3.4：取步骤3.1～步骤3.3得到时域下的特征、频域下的特征和时频域下的特征组成故障特征向量：峰值指标C、波形指标S、脉冲指标I、裕度指标L和峭度指标K、衰减系数β、功率谱密度Sx、输入输出的胡功率谱密度Sxy、均方频率MSF、频率方差VF、W30、W31、W32、W33、W34、W35、W36、W37，共18个故障特征向量，组成特征矩阵XL；
步骤4：分别对首端和末端上位机采集的高频低频混合历史数据经过步骤2～步骤3后，得到特征矩阵XL,训练COIR-LDA判别，即协方差核切片逆回归算法-线性判别分析，COIR算法是由KSIR算法改进得到的，分别c类工况下的数据特征，KSIR算法即核切片逆回归算法，具体步骤包括步骤4.1～步骤4.8：
步骤4.1：由核函数定义的非线性变换，核函数取高斯核函数其中，exp就是e指数函数，||·||表示绝对距离，x为输入，x指代特征矩阵XL，xi为核函数中心，σ为函数的宽度参数，将18维测量特征矩阵XL映射至高维特征空间F，得λ为特征值；
步骤4.2：Y切成c份，分别为I1,…,Ic；Y＝f(b1x,b2x,b3x…bnx)，计算切片内I1,…,Ic的均值M；
其中，表示n*c的0/1指示矩阵,每一行只有一个元素是1，表示混合历史数据属于某一切片，其余为0，其中，c表示工况的数量，n是样本数量，b1,b2…,bn是投影向量b，δh是个函数；
步骤4.3：计算高维特征空间F中每段的均值向量和加权协方差阵即：
V＝MPMT (17)
我们把矩阵V的前d个特征向量记成为切片内样本占全部样本的比例，l＝1,…,d，n是样本数量；
步骤4.4：中心子空间的方向可以通过获得，x输入的特征,也就是XL,n样本数量；
步骤4.5:利用核函数的内积运算Kij(xi,xj)＝<φ(xi),φ(xj)>求解特征值分解Vv＝λv，可以得到下面对偶方程：
得到特征值λ，其中n是样本数量，h表示工况的数量，α＝[α1，…αn]T；
步骤4.6：投影向量b从特征向量v得到，为了求解令其中β＝[β1,…,βn]T，并且两侧同乘以即可以得到非线性子空间β的闭式解：
步骤4.7:非线性中心子空间即可以通过d个从构成的基函数表示，如果给定一个新的样本x*，x*指的是XL，它的低维表示z*∈Rd可以通过将φ(x*)投影中心子空间得到，也就是说，z*的l个元素是：其中，k*＝[kx(x1,x*),…,kx(xn,x*))]T，Rd是前d个特征向量组成的R维空间；
步骤4.8:将剩余低频历史数据及高频历史数据经过步骤2进行预处理，步骤3特征提取，并进行步骤4，最终得到2c种高低频历史数据的非线性中心子空间z*,每个z*有对应的
步骤5：改进步骤4所述的(18)式，得到COIR算法，即协方差核切片逆回归算法，具体步骤包括步骤5.1～步骤5.5：
步骤5.1：对任意的函数f∈Rx，如果存在一个函数g∈Rx，使得对所有的y，有E[f(x)|y]＝g(y)，其中，E是求期望，var是求方差，那么：
步骤5.2：使用方差公式可以得到：
var[E[φ(x)|y]]＝var(φ(x))-E[var[φ(x)|y]] (21)
步骤5.3：通过步骤5.1和步骤5.2，步骤5.3的估计量可以写成：
其中，
步骤5.4：COIR算法需要对上式进行特征值分解，给定观测样本再令φy＝[φ(y1),…,φ(yn)]，那么步骤5.4的估计量可以写成：
步骤5.5：同样，特征值分解var[E[φ(x)|y]]v＝λv可以写成下式：
其中，ε为人工设定的参数；
步骤6：实时采集管道高频及低频数据，通过步骤2进行预处理，以秒为单位，送入步骤3进行时域、频域和时频域特征提取，再进入步骤4，通过训练得到的COIR-LDA判别模型，根据得到的特征向量与建立的四种工况情况的数据库进行比对，根据LDA判别函数得到实时数据对应于某类工况，具体包括步骤6.1～步骤6.2；
步骤6.1：对于n个步骤2得到的小波包重构后的信号的c类问题,假设ωi(i＝1,2,…,c)类模式的先验概率为p(ωi),显然有
步骤6.2：将步骤5的COIR算法作用于采集的数据X,提取和输出F的中心子空间，步骤4求出来的投影向量b的空间表示的样本数据Z,然后构造第i类别LDA的判别函数：
式中zt∈Rr为F的中心子空间中带判别的分类属性矢量，则为第i类训练样本的在F的中心子空间的平均分类属性矢量，为训练样本在F的中心子空间的类内离散度矩阵，其中Ww,z是矩阵xw xz的矩阵，xz为第i类特征样本xw是输入的样本，得到c种工况下对应的Li(zt)，建立管道泄漏检测的数据库；
步骤7：实时采集管道高频及低频数据，通过步骤2进行预处理，然后送入步骤3进行时域、频域和时频域特征提取，再进入步骤4，通过训练得到的COIR-LDA判别模型，根据得到的特征向量与建立的四种工况情况的数据库进行比对，计算两者之间的相似程度；
相似程度的计算公式如下所示：
其中，Sk为相似度，C(x)new为新采集的数据，为第k种类型的测试样本，为k类样本的一个集合；
若Si≥0.85，则当前数据为第k类数据；否则，结合现场状况的实际分析及人为诊断确定故障的类型，并将其收入管道泄漏检测的数据库中；
步骤8：根据步骤7，若实时数据判定为大泄漏信号或小泄漏信号时，则通过高频信号和低频信号混合计算泄漏点位置，对管道泄漏进行定位，否则返回步骤6，继续实时监测；
管道泄露进行定位的具体过程包括步骤8.1～步骤8.5：
步骤:8.1：根据首端和末端传感器的时间差值计算泄漏点，如下式所示，
式中，K1为泄漏点到首端超高频传感器的距离；S1为管线首端和末端超高频传感器之间的距离；τ1为高频信号传到首端和末端超高频传感器的时差；o1为高频信号在管线内的传播速度；
步骤8.2：在泄漏点的上游站及下游站各选取泄漏发生点前30秒的高频数据均值和泄漏点后15秒的高频数据均值，计算上下游站高频信号变化比，分别如下两式所示，
其中：δ1、δ2分别为首端和末端的高频信号变化比；X1+、X2+分别为泄漏发生点前30秒首端和末端的高频数据均值；X1-、X2-分别为泄漏发生点后15秒首端和末端的高频数据均值；
则管道泄漏时总高频信号变化比δp为：其中，μ为由现场管道长度和正常输送管道的大小确定的参数；
步骤8.3：对低频信号的特征值在首端即末端的时间差值，计算泄漏点，如下式所示，
式中，K2为泄漏点到首端加速度传感器的距离；S2为管线首端和末端加速度传感器之间的距离；τ2为声波传到首端和末端加速度传感器的时差；o2为低频信号在管线内的传播速度；
步骤8.4：在泄漏点的首端及末端各选取泄漏发生点前5秒的低频数据均值和泄漏点后3秒的低频数据均值，计算首端和末端低频信号变化比，分别如下两式所示，
其中，δ3、δ4分别为首端和末端的低频信号变化比；Y1+、Y2+分别为泄漏发生点前5秒首端和末端的低频信号数据均值；Y1-、Y2-分别为泄漏发生点后3秒首端和末端的低频数据均值；
则管道泄漏时总低频信号变化比δs为：
步骤8.5：根据泄漏点到上游站超高频传感器的最终距离为：
其中，T为低频信号的周期。
CN201810794586.1A 2018-07-19 2018-07-19 基于高低频混合检测的管道泄漏监测系统及方法 Active CN108870091B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810794586.1A CN108870091B (zh) 2018-07-19 2018-07-19 基于高低频混合检测的管道泄漏监测系统及方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810794586.1A CN108870091B (zh) 2018-07-19 2018-07-19 基于高低频混合检测的管道泄漏监测系统及方法

Publications (2)

Publication Number Publication Date
CN108870091A CN108870091A (zh) 2018-11-23
CN108870091B true CN108870091B (zh) 2019-09-10

Family

ID=64303476

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810794586.1A Active CN108870091B (zh) 2018-07-19 2018-07-19 基于高低频混合检测的管道泄漏监测系统及方法

Country Status (1)

Country Link
CN (1) CN108870091B (zh)

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109798449B (zh) * 2018-12-20 2021-08-10 国能大渡河沙坪发电有限公司 基于机器视觉单元神经网络的供水系统巡检方法及系统
CN109922082B (zh) * 2019-04-10 2021-09-21 杭州数梦工场科技有限公司 流量异常的检测方法及装置和计算机可读存储介质
CN110532635B (zh) * 2019-08-05 2023-03-28 上海第二工业大学 一种基于时域的管道泄漏检测算法
CN112817291B (zh) * 2019-11-15 2022-03-08 中国科学院沈阳自动化研究所 基于混合特性评价和子空间分解的层次故障监测方法
CN111062127B (zh) * 2019-12-16 2023-08-04 辽宁石油化工大学 管道漏点的检测方法及装置、存储介质、终端
CN113048402A (zh) * 2019-12-27 2021-06-29 柯思科技股份有限公司 中低压管线监控系统及其方法
CN111695465B (zh) * 2020-06-01 2024-03-05 浙江英集动力科技有限公司 基于压力波模式识别的管网故障诊断与定位方法及系统
CN112413413B (zh) * 2020-11-20 2022-03-01 成都艾斯皮埃尔科技有限公司 结合深度学习和多次测量技术的管道泄漏监测及定位方法
CN116501623B (zh) * 2023-04-13 2024-09-17 四川九洲电器集团有限责任公司 一种用量子计算方法实现的自动测试系统及方法
CN116557793B (zh) * 2023-07-10 2023-12-05 中建安装集团有限公司 一种融合压力传感和温度传感的供热管道运行状态监测系统及方法
CN117053129B (zh) * 2023-10-12 2023-12-15 山西瑞赛科环保科技有限公司 一种液氯的安全风险预警平台、方法、电子设备以及介质

Family Cites Families (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007163255A (ja) * 2005-12-13 2007-06-28 Mitsubishi Materials Corp 漏水検出方法及び漏水検出システム並びにｒｆｉｄタグ
CN102023186B (zh) * 2010-12-29 2013-07-31 钢铁研究总院 电磁超声探头以及使用该电磁超声探头检测管道的方法
CN102563362B (zh) * 2011-12-31 2013-08-28 杭州哲达科技股份有限公司 压缩空气系统管网泄漏智能检测方法及系统
CN102797979B (zh) * 2012-08-29 2014-07-02 上海海事大学 地下管道泄漏点的探测装置及方法
CN105889763B (zh) * 2014-10-10 2017-12-15 保定市金迪科技开发有限公司 一种管道泄漏的检测设备及其检测方法
CN104930353B (zh) * 2015-07-10 2017-08-15 浙江大学 一种城市供水管网的调压减漏评估系统
CN105760839A (zh) * 2016-02-22 2016-07-13 重庆大学 基于多特征流形学习与支持向量机的轴承故障诊断方法
CN106523928B (zh) * 2016-11-24 2018-08-03 东北大学 基于声波实时数据二级筛选的管道泄漏检测方法
CN107063582A (zh) * 2016-12-30 2017-08-18 天津市誉航润铭科技发展有限公司 一种输水管道泄露监测传感器
CN207316488U (zh) * 2017-10-30 2018-05-04 成都川盛塑胶有限公司 一种用于管道输送系统漏液或漏气的远程无线自动监测器
CN107835261A (zh) * 2017-12-15 2018-03-23 佛山三维二次方科技有限公司 基于物联网的油气管道远程数据采集系统

Also Published As

Publication number Publication date
CN108870091A (zh) 2018-11-23

Similar Documents

Publication Publication Date Title
CN108870091B (zh) 基于高低频混合检测的管道泄漏监测系统及方法
CN106338385B (zh) 一种基于奇异谱分解的旋转机械故障诊断方法
Du et al. Wavelet leaders multifractal features based fault diagnosis of rotating mechanism
Yan et al. Improved Hilbert–Huang transform based weak signal detection methodology and its application on incipient fault diagnosis and ECG signal analysis
CN110454687A (zh) 一种基于改进vmd的管道多点泄漏定位方法
CN104198184A (zh) 基于第二代小波变换与bp神经网络的轴承故障的诊断方法
CN101201386B (zh) 一种模拟集成电路参数型故障的定位方法
CN104655425A (zh) 基于稀疏表示和大间隔分布学习的轴承故障分类诊断方法
CN103020478A (zh) 一种海洋水色遥感产品真实性检验的方法
CN103115789A (zh) 金属结构损伤剩余寿命的第二代小波支持向量机评估方法
CN105488520A (zh) 基于多分辨奇异谱熵和svm的泄漏声发射信号识别方法
CN112213687B (zh) 基于伪异常点辨识的关口电能表数据异常检测方法及系统
CN114152825A (zh) 变压器的故障诊断方法、装置和变压器的故障诊断系统
CN113325357A (zh) 基于输出时间序列差分的电压互感器误差评估方法及系统
Qian et al. A new health indicator for rolling bearings based on impulsiveness and periodicity of signals
CN103335840A (zh) 一种矿用钻机变速箱故障智能诊断方法
Hai et al. Rolling bearing fault feature extraction using non-convex periodic group sparse method
CN109815940A (zh) 小波包能量谱法损伤识别方法
Wang et al. Detrended Fluctuation Analysis and Hough Transform Based Self‐Adaptation Double‐Scale Feature Extraction of Gear Vibration Signals
Bernardi et al. Characterization of the Ethanol-Water Blend by Acoustic Signature Analysis in Ultrasonic Signals
Jiao et al. A gas regulator fault detecting method based on acoustic emission technology
CN116772122A (zh) 一种天然气管道泄漏故障诊断方法、系统、设备及介质
Zhao et al. Pipeline leak fault feature extraction based on wavelet packet analysis and application
Ogundare et al. Review of fault detection techniques for health monitoring of helicopter gearbox
Dong et al. EMD Method for Minimizing the Effect of Seasonal Trends in Detrended Cross‐Correlation Analysis

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant