一种融合EWT的ISSA-VMD-KU的uwDAS振动信号去噪方法

本发明涉及uwdas振动信号去噪，特别是涉及一种融合ewt的issa-vmd-ku的uwdas振动信号去噪方法。

背景技术：

1、分布式光纤声波传感(distributed optical fiber acoustic sensing，简称das)系统可以实现振动和声场连续分布式探测，其中，基于菲索干涉的超弱光纤布拉格光栅(ultra weak fiber bragg grating,uwfbg)和增强相位敏感光时域反射(phase-sensitive optical ti me domain reflector，-otdr)技术的das系统，简称uwdas。uwfbg反射光的强度相比于瑞利后向散射光高2～3个数量级，大大提高了系统的信噪比、动态范围和灵敏度。uwdas具有大范围、高分辨率、高灵敏度、宽频带、抗电磁干扰等优点，广泛应用于管道安全监测、水声探测、地震检波、电力电缆监控等领域。

2、窄线宽激光器发出的连续光经过调制器调制为脉冲光，再经光放大器放大，并通过环形器进入超弱光纤光栅阵列，各光栅的反射光反射回环形器，随后进入迈克尔逊干涉仪，相邻光栅的反射光在3×3耦合器处干涉，形成三路相位差为120°的干涉光信号，经光电转换后，通过上位机解调相位变化，从而实现对uwfbg之间发生的扰动信号的监测。

3、uwdas受激光器、放大器、探测器等光学器件噪声的影响，还存在电噪声、或由于光缆过长、人为或环境等各种干扰造成的噪声，使采集的数据质量下降，对目标振动信号的处理和分析产生严重影响，因此对uwdas采集的信号进行去噪以提升信号质量、提高信噪比是该领域的研究热点。

4、当前去噪方法主要包括带通滤波法、低通滤波法、小波阈值法、经验小波变换(empirical wavelet transform，ewt)以及变分模态分解(variational modaldecomposition，vmd)去噪方法等。带通、低通滤波法适合平稳信号的处理，而uwdas信号具有非平稳特性；小波阈值法需预设小波基、小波分解层数以及合适的阈值，难以进行自适应滤波；ewt需要预设参数和分割频谱的方式，参数设置不合适会出现无用分量和频谱混叠问题；vmd算法可一定程度抑制模态混叠问题、可处理非平稳、非线性信号，被广泛应用于信号降噪领域。

5、现有技术的缺陷与不足：

6、(1)vmd算法可将任意信号分解为若干个本征模态函数(intrinsic modefunction，imf)和一个余项(residue，res)。uwdas采集的信号经过vmd分解后，得到从高频到低频不同中心频率的分量imf，噪声主要分布于高阶imf中。vmd分解效果受参数分解层数k以及惩罚因子alpha影响较大，当k值过大时，会导致分解结果出现混频现象，出现过多的分量，当k值较小时，又不利于捕获各个模态的正确中心频率，惩罚因子alpha控制数据保真度，不适当的惩罚因子不利于模态分量的保存。vmd算法的关键在于如何选取合适的分解层数k以及惩罚因子alpha进行后续有效分量imf的筛选。

7、(2)当前的imf筛选策略难以准确筛选有效的imf：相关系数筛选法适用于低信噪比的信号，对于复杂的uwdas信号，其中含噪种类多，相关系数无法正确筛选有效的imf；排列熵法的计算依赖于参数的选择，并且其基于统计特征，可能会舍弃imf中的一些重要信息。因此，合适的筛选方法十分重要。

8、(3)ewt对信号的局部特征提取能力弱，无法有效地处理非平稳信号；使用的基函数较少，不足以表示复杂的信号；需要预设参数和分割频谱的方式，参数设置不合适会出现无用分量和频谱混叠问题，仅使用ewt去噪效果并不理想。

9、(4)麻雀搜索算法(sparrow search algorithm，ssa)以最小包络熵为目标函数可用于优化vmd参数，但是容易陷入局部最优解，搜索能力较弱，寻优的k值总靠近设置寻优范围的上限或下限。因此，需要改进算法提高种群多样性及算法搜索能力。

技术实现思路

1、为解决上述技术问题，本发明提供一种融合ewt的issa-vmd-ku的uwdas振动信号去噪方法，相比传统的emd、vmd，该方法能够有效抑制模态混叠现象，分解出更有效的分量imf，提升了去噪效果。

2、本发明采取的技术方案为：

3、一种融合ewt的issa-vmd-ku的uwdas振动信号去噪方法，包括以下步骤：

4、步骤1：使用uwdas采集振动信号y(t)；

5、步骤2：对振动信号y(t)进行ewt预处理，得到m个mra模态分量；

6、步骤3：计算各阶mra与原始信号的皮尔森相关系数ρ，选取ρ＝[ρl,ρh]阈值范围的分量进行重构，得到重构后振动信号y′(t)；

7、步骤4：使用vmd分解重构后振动信号y′(t)，将最小包络熵作为适应度函数，使用issa优化vmd的分解层数k和惩罚因子alpha，直至达到最大迭代次数，输出最优参数组合[k_best，alpha_best]；

8、步骤5：使用最优参数[k_best，alpha_best]再对y'(t)进行vmd分解，得到k_best个imf与余项res；

9、步骤6：计算各阶imf的峭度值ku，选取ku＝[kul，kuh]阈值范围内的imf重构，得到去噪后的uwdas振动信号y″(t)。

10、所述步骤2包括以下步骤：

11、s2.1:对输入振动信号y(t)进行快速傅里叶变换，依据香农准则将信号频率规范化为[0,π]，并将其划分出n个连续部分，以相邻极大值点的中点ωn作为分割区间的边界，第n个频谱区间为[ωn-1,ωn],ω0＝0,ωn＝π；

12、s2.2:计算近似系数和细节系数，公式如下：

13、

14、

15、上式中，为近似系数；为细节系数；<y,φ1>表示输入的振动信号与经验小波尺度函数的内积，<y,ψn>表示输入的振动信号与经验小波函数的内积，φ1表示表示第1个经验小波尺度函数，ψn表示第n个经验小波函数；表示φ1傅里叶变换的复共轭；表示ψn傅里叶变换的复共轭；n＝1,2,3…n+1；表示y(t)的傅里叶变换；f-1[·]表示求傅里叶逆变换；和分别为经验小波尺度函数和经验小波函数，依据littlewood-paley和meyer的小波构造理论定义为：

16、

17、

18、上式中，τn＝γωn表示相邻两个频带中间的过渡区间，0<γ<1；ω为振动信号频率；β为关于x的函数：β(x)＝x4(35-84x+70x2-20x3)

19、s2.3:根据近似系数和细节系数，y(t)被分解成一个近似分量y0(t)和m＝n-1个细节分量ym(t)，m＝1,2,3…m：m表示第m个细节分量；m表示细节分量数。

20、

21、

22、上式中，φ1(t)表示第1个经验小波尺度函数；ψm(t)表示第m个经验小波函数，为近似系数；为细节系数

23、所述步骤3中，计算各阶mra即步骤2.3中的细节分量ym(t)与原始信号的皮尔森相关系数，公式如下：

24、

25、上式中：cov表示协方差，σ表示标准差；ρm表示皮尔森相关系数；σy(t)表示输入信号y(t)的标准差；表示细节分量ym(t)的标准差；cov()表示输入信号y(t)与细节分量ym(t)的协方差；

26、根据阈值ρ＝[ρl,ρh]筛选mra分量，ρl表示皮尔森相关系数阈值的下限，ρh表示皮尔森相关系数阈值的上限。重构公式如下：

27、

28、上式中，θ为根据阈值ρ＝[ρl,ρh]判断的有效mra的集合；[ym(t)]p表示属于有效mra集合的分量；y'(t)为ewt预处理后的uwdas振动信号。

29、所述步骤4包括以下步骤：

30、s4.1:初始化参数，优化vmd参数个数设置为2，并设置其范围，分解层数k设置为[kl，kh]，kl表示设置分解层数k范围的下限，kh表示设置分解层数k范围的上限；惩罚因子alpha设置为[al，ah]，al表示设置惩罚因子alpha范围的下限，ah表示表示设置惩罚因子alpha范围的上限；设置最大迭代寻优次数为itermax，设置种群数量为np。

31、s4.2:初始化麻雀种群，设置发现者比例pdrate、预警者比例sdrate、安全阈值st；

32、s4.3:对y′(t)进行vmd分解，计算初始种群的适应度值即imf的最小包络熵；

33、s4.3.1:初始化本征模态函数各个模态分量的中心频率以及拉格朗日乘子

34、s4.3.2:更新ωk以及公式如下：

35、

36、

37、

38、上式中，n′为迭代次数，n′初始化为1；为第n′+1次迭代的模态分量；为y′(t)的频域表示；为第n次迭代第i个模态分量；为第n′次迭代的拉格朗日乘子，拉格朗日乘子初始值为0；为第n′+1次迭代的中心频率；第n′+1次迭代的拉格朗日乘子；τ′为拉格朗日乘子的更新参数；k表示模态分量数；i表示第i个模态分量；α表示二阶惩罚因子；ωn′k表示第n′次迭代模态分量对应的中心频率；表示第n′+1次迭代模态分量对应的中心频率。

39、s4.3.3:判断是否达到如下收敛条件：

40、

41、上式中：表示第n′+1次迭代的模态分量，表示第n′次迭代的模态分量。∈为收敛容差，表示2范数的平方，若达到收敛条件则输出当前k个模态分量以及对应中心频率，否则重复步骤4.3.2。

42、s4.3.4:计算初始种群的适应度值，首先计算包络信号：

43、

44、上式中，aj(s)表示第j个分量通过希尔伯特变换后得到的包络信号；j＝1,2,3…k，s＝1,2,…,s，s为采样点数；imfj(s)为vmd分解后的imf分量；h[imfj(s)]表示imfj(s)的希尔伯特变换。

45、对包络信号归一化：

46、

47、上式中，pj(s)是第j个分量包络的归一化形式。

48、计算包络熵ep：

49、

50、上式中，ep(j)为第j个分量包络熵；lgpj(s)表示对第j个分量分量包络归一化后求其对数。

51、以包络熵的最小值minep为适应度函数:

52、minep＝min{ep(1) ep(2) … ep(k)}

53、ep(1) ep(2) … ep(k)分别表示第1个到第k个分量的包络熵。

54、s4.4:更新发现者位置、追随者位置和警觉者位置，根据步骤4.3中的公式计算此时的适应度值。

55、s4.4.1:在每次迭代过程中，在发现者位置更新公式中引入正弦余弦思想，同时，引入一个非线性正弦学习因子，用来提升局部探索能力。改进后的发现者位置公式为：

56、ω＝ωmin+(ωmax-ωmin)·sin(tπ/itermax)

57、

58、上式中，r1、r2分别为[0,2π]、[0,2]的随机数；学习因子取ωmin＝0.3、ωmax＝1。t′表示当前迭代的次数；分别表示迭代t′+1和t′次时第i′个麻雀第j′个维度的值(i′＝1,2…np,j′＝1,2)，xbest表示全局最优位置，r2(r2∈[0,1])和st(st∈[0.5,1.0])分别表示预警值和安全阈值。当r2<st，表示麻雀种群周围没有捕食者，此时发现者进入广域搜索模式；当rk≥st，表示一些麻雀已经发现捕食者，此时麻雀种群迅速飞往安全区域。

59、s4.4.2:在每次迭代过程中，追随者位置更新公式如下：

60、

61、上式中，是发现者所在的最佳位置；表示当前全局最差的位置；a表示一个1×2的矩阵，其中每个元素分配1或-1；a+＝at(aat)-1，at为a的转置，a+为a的moore-penrose逆矩阵；l为一个1×2的矩阵，其中内部元素都为1；当i′>np/2，表示适应度值较差的第i′个搜索者最有可能挨饿。

62、s4.4.3:意识到危险的麻雀即警觉者在种群中的初始位置是随机产生的，其位置更新公式如下：

63、

64、上式中，表示当前全局最优的位置；β′作为一个步长控制参数，是一个均值为0，方差为1的服从正态分布的随机数；fi′是当前的适应度值；fg和fw是当前全局最佳和最差适应度值；ε是最小的常数，以避免零除法误差；当fi′>fg时，表示警觉者处于种群的边缘；当fi′＝fg时，表明位于种群中间的麻雀意识到了危险，需要向其他麻雀靠近；k′∈[-1,1]表示麻雀移动的方向，也是步长控制系数。进行步长因子动态调整增大算法中位置更新的随机性，公式如下：

65、

66、

67、上式中，rand为[0,1]的随机数。

68、s4.5:t′加1，重复步骤4.4，直至t′达到最大迭代次数itermax，输出全局最小适应度值对应的最优参数组合[k_best，alpha_best]。

69、所述步骤5中，使用步骤4输出的最优参数组合[k_best，alpha_best]再对y'(t)进行vmd分解，得到k_best个imf与一个余项res。

70、vmd分解公式在s4.3中已经提到，上述用最优的参数再进行一次vmd分解。

71、所述步骤6中，计算各阶imf的峭度值ku，选取ku＝[kul，kuh]阈值范围内的im f重构，得到去噪后的uwdas振动信号y″(t)。

72、峭度值ku具体计算公式如下：

73、

74、上式中，ss＝1,2,…,s，imf(ss)是使用最优参数进行vmd分解后imf分量；μ是信号的平均值；rms是均方根值，ku是峭度值，ku反应了振动信号的冲击性。

75、根据ku＝[kul，kuh]阈值范围，kul表示设置峭度值ku阈值范围的下限，kuh表示设置峭度值ku阈值范围的上限；选择有效分量imf后重构uwdas振动信号，

76、峭度值ku根据信号特征确定，具体重构公式如下：

77、

78、上式中，为根据峭度值判断的有效imf的阶数的集合；y”(t)为去噪后重构的uwdas振动信号；imfg表示属于有效imf集合的分量。

79、本发明一种融合ewt的issa-vmd-ku的uwdas振动信号去噪方法，技术效果如下：

80、1)本发明使用ewt对uwdas振动信号进行预处理。ewt与vmd的结合可以更准确的分解uwdas振动信号，捕捉不同频率范围的振动成分，选择性的去除其中的噪声，达到更好的去噪效果。

81、2)本发明使用issa优化vmd的分解层数k以及惩罚因子alpha，相比优化前的vmd算法，可分解出更合适的层数，有效抑制模态混叠，可以提高分解精度。在issa中加入非正弦因子改进发现者位置更新，加入动态因子调整改进警觉位置更新，增强算法的搜索能力，避免陷入局部最优解；

82、3)峭度值是一种统计量，对冲击信号特别敏感，可以反映振动信号的冲击特性。本发明使用峭度值筛选imf，作为区分噪声和有用信号的依据，提取有用信号成分、减少噪声，并且具有较好的自适应性，以筛选有效imf。

83、4)本发明步骤3中，利用皮尔森相关系数可以更好准确地评估分量与参考信号之间的相似性，更有利于保留的有效信号成分。与ewt结合使用，可以更加灵活的处理和分析信号；

84、5)本发明步骤4中，将最小包络熵作为适应度函数可以有效地引导优化过程，使得优化的vmd能够更精确地提取信号的主要成分。结合改进的ssa算法来优化vmd的方法具有较好的信号处理和分析性能；