中国指挥与控制学会会刊 
相位敏感型光时域反射计(Φ-OTDR, Phase Sensitive Optical Time-Domain Reflectometer)技术因其可以探测振动信号而发展形成了分布式光纤声波传感技术(DAS, Distributed Fiber Acoustic Sensing)[1]。Φ-OTDR传感器利用相干光源的脉冲感知光纤全段发生的振动。当光纤某一位置发生扰动时,后向瑞利散射光会产生相位变化,进而影响接收端测量的后向散射光强度。这种原理可用于长距离分布式振动测量。
自1993年H.F. Taylor首次提出Φ-OTDR以来,该技术就因分布式特性受到广泛关注[2]。在地质勘测、土木工程、电力设施和周界安防等领域得到了广泛应用[3]。然而,由于Φ-OTDR自身噪声和环境噪声的影响,采集数据呈现非平稳非高斯的噪声特征[4]。为此,需要一种能够在低信噪比条件下降噪的方法。研究人员提出了各种信号处理方法来对Φ-OTDR信号进行去噪,主要包括短时傅里叶变换、小波阈值去噪和经验模态分解[5]。短时傅里叶变换和小波去噪方法可以对Φ-OTDR信号进行去噪,但噪声信号和有用信号通常不能在单一维度上有效区分,有用信号在去噪过程中容易丢失。短时傅里叶变换只能进行频域滤波。小波去噪方法需要手动选择小波函数且分解尺度适应性程度不高[6]。 经验模态分解(EMD, Empirical Mode Decomposition)具有较高的适应性,已在Φ-OTDR信号的降噪中得到应用[7]。然而,EMD方法存在模态混叠现象。因此,对于Φ-OTDR信号的降噪,需要进一步改进现有的方法,以提高降噪效果。
针对以上问题,本文提出一种基于样本熵和小波阈值去噪的改进自适应噪声完全集成经验模态分解的降噪方法。首先,使用ICEEMDAN分解输入信号,通过计算每个固有模态函数(IMF,Intrinsic Mode Functions)的样本熵(SampEn,Sample Entropy)来判断含噪的IMF。其次,利用小波阈值去噪对含噪IMF进行去噪,进一步提取含噪IMF中包含的有用信号。最后,将降噪IMF与剩余的IMF结合并重建,获得去噪结果。通过实测信号验证了所提方法对信号的降噪效果。与小波阈值去噪(WTD,wavelet threshold denoising)方法、谱减法(SS,Spectral Subtraction)、EMD方法、EMD-WTD 方法和ICEEMDAN方法相比,本文所提方法的性能更好。
1 算法原理
1.1 Φ-OTDR传感原理
Φ-OTDR中空间分辨率与激光脉冲宽度之间的关系可以表示为
Δz=
其中,c代表真空中的标准光速,nRef代表所使用光纤的折射率,Tp代表激光器发射激光的脉冲宽度,Δz代表空间分辨率。脉冲宽度会影响空间分辨率和动态范围,较窄线宽的激光脉冲意味着较低的空间分辨率。然而,较窄线宽的激光脉冲携带较低的光功率,从而导致较低的动态范围。本文采用基于自适应噪声、相关函数和小波阈值去噪的改进完全集成经验模态分解来提高信噪比和动态范围,这在窄线宽激光脉冲下更为显著。通常情况下,在许多现有的Φ-OTDR配置下,信噪比的提高要么是以增加复杂度为代价,要么是以减小测量动态范围为代价。
1.2 Φ-OTDR噪声分析
Φ-OTDR中的两个主要噪声源是光学噪声和自然噪声,其中,光学噪声包括相干噪声和ASE噪声[8]。光电转换及采集过程也会带来噪声,但是对信号的影响很小,不需要针对其进行特殊处理。相干噪声的产生是因为使用了相干长度大于传感距离的脉冲光源。当使用非相干光源作为传统Φ-OTDR的信号源时,接收端后向散射信号的噪声将导致信噪比与信号通带内的空间分辨率与传感距离的比值成正比。而Φ-OTDR传感器中,光源的相干长度可以大于测量距离。脉冲持续时间内的后向散射信号会与在同一波段存在的“0”能级的后向散射信号进行相干混合。因此,在接收机上会有一个额外的跳变,它会产生噪声。
ASE噪声源是光纤放大器和半导体光放大器中的一种噪声源,产生于激活粒子从激发态返回基态并放大光信号时的随机非相干自发辐射。该放大器自发辐射噪声,即ASE噪声,伴随光信号产生,对系统性能产生重要影响。
Φ-OTDR光纤传感监测技术以其高灵敏度而闻名,能够探测到各种扰动产生的振动信号。这种高灵敏度也意味着在实际的传感信号采集中,可能会受到环境振动所引起的干扰的影响。这些干扰对光纤传感振动信号的质量产生一定的影响。自然噪声源自于这种高灵敏度下放大的周围环境噪声。
1.3 ICEEMDAN原理
1998年由Huang等人提出了EMD,其在面对分析非线性非平稳信号的问题时具有显著的优势[9]。然而,EMD算法也存在模式混叠的问题。为克服模态混叠,Wu和Huang于2009年提出了集成经验模态分解(EEMD)[10],该方法利用滤波器组的行为和白噪声频谱均匀分布的统计特性,有效抑制由于间歇性高频分量等因素引起的模态混叠现象。为解决EEMD噪声残留的问题,2011年Torres等人提出了具有自适应噪声的完整集成经验模态分解(CEEMDAN)。在分解的初始阶段,CEEMDAN存在产生虚假模式的倾向。为了消除CEEMDAN的虚假模式,Colominas等人提出了ICEEMDAN[11],这一算法能够获得更具物理意义的模式组件。ICEEMDAN方法选择了一种特殊的原则引入白噪声。该方法将白噪声进行EMD分解,选择其中第K个IMF分量加入分解过程中。
根据Colominas等人提出的算法改进思路,ICEEMDAN的具体分解步骤如下。其中,<·>是平均值,x是原始信号,x(i)是加入白噪声后的信号,ω(i) 是均值和单位方差为零的白噪声,M(·) 是信号的局部平均值, Ek(·) 是EMD分解后的第k个模态分量,βk是第k次调整时添加噪声的系数。ICEEMDAN的具体分解过程如下:
步骤1:构造x(i)=x+β0E1(ω(i)),计算x(i)的局部均值M(x(i))以获得第一次分解的残差r1=<M(x(i))>。
步骤2:减去第一个残差,得到从原始信号处获得的第一个分量IMF1=x-r1。
步骤3:构造第二个残差r2=<M(r1+β1E2(ω(i)))>,计算第二个分量IMF2=r1-r2=r1-<M(r1+β1E2(ω(i)))>。
步骤4:对于k=3,4,…,K,减去第k个残差以获得第k个分量的值IMFk=rk-1-rk。
步骤5:重复步骤4,直到获得所有分量。
选择常量βk-1来调整添加噪声和残差的信噪比。当k=1时,β0=ε0std(x)/std(E1(ω(i)));当k>2时,βk=ε0std(rk),其中,ε0是输入信号和加入的首个噪声所需信噪比的倒数,std(·)是标准差计算。
1.4 样本熵
样本熵是由Richman等学者提出的一种统计量,与近似熵有所区别,它不考虑自身匹配的计数。样本熵通过衡量信号中新模式产生的概率来评估时间序列的复杂性[12]。当新模式产生的概率较高时,序列的复杂性也相应增加。
对于由N个数据组成的时间序列 =x(1),x(2),…,x(N),样本熵的计算方法如下:
步骤1:按照序号组成的一组维数为m的向量序列,Xm(1),Xm(2),…,Xm ,其中Xm(i)= ,1≤i≤N-m+1。这些向量代表从第i点开始的m个连续x的值。
步骤2:定义向量Xm(i)和Xm(j)之间的距离d[Xm(i),Xm(j)]为两者对应元素中最大差值的绝对值。即
d[Xm(i),Xm(j)]=maxk=0,…,m-1
步骤3:对于向量Xm(i),统计向量Xm(i)与向量Xm(j)之间距离小于r的j(1≤j≤N-m,j≠i)的个数,记为Bi。对于1≤i≤N-m,定义:
(r)= Bi
步骤4:定义B(m)(r)为
B(m)(r)= (r)
步骤5:增加维数到m+1,计算Xm+1(i)与Xm+1(j)(1≤j≤N-m,j≠i)距离小于等于r的个数,记为Ai。1≤i≤N-m,定义:
(r)= Ai
步骤6:定义A(m)(r)为
A(m)(r)= (r)
B(m)(r)是两个序列在相似容限下r匹配m个点的概率,是两个序列匹配m+1个点的概率。样本熵的定义为
SampEn =-ln
1.5 小波阈值去噪
小波阈值去噪由Donoho等学者提出[13]。该方法通过小波变换对信号进行分解和重构,结合阈值处理的策略有效地去除噪声成分。小波阈值去噪的过程如下:
步骤1:选择合适的小波函数和分解尺度来分解输入信号。
步骤2:在小波阈值去噪中,需要选择适当的阈值函数和阈值。若小波系数低于阈值,则将其归属于噪声信号,并将其设为零;若小波系数不低于阈值,则认定其为有效信号的成分。这一策略有助于有效区分噪声和有用信号,并对小波系数进行相应的调整。
步骤3:去噪信号由阈值小波系数的逆小波变换得到。阈值的选择对算法效果至关重要。阈值的选取确定小波系数应该被保留或丢弃,决定了去噪效果。
硬阈值函数会保留大于或等于阈值的小波系数,并将小于阈值的小波系数设为零。软阈值函数也将小于阈值的小波系数设置为零,但软阈值函数会将不小于阈值的小波系数缩小。硬阈值函数和软阈值函数分别可以表示为:
=
=
式中,WJ,K是处理前的小波系数; 是处理后的小波系数;λ是设置的阈值;Sgn(·)是符号函数;J是最大分解尺度,j=1,2,…,J;K是小波系数的长度,k=1,2,…,nj,nj=N/2J-j+1;N是信号的长度。
在实际使用过程中,硬阈值函数的不连续性可能会导致去噪信号出现跳点。软阈值函数的去噪效果更加连续和平滑。针对噪声IMF分量,出现跳点对信号的影响是很大的。
1.6 ICEEMDAN-样本熵-WTD方法
考虑ICEEMDAN、样本熵、WTD的优点,本文提出了ICEEMDAN-样本熵-WTD方法,具体步骤如下:
步骤1:输入信号由ICEEDAN分解为一系列IMF。
步骤2:计算每个IMF分量的样本熵值。如果大于或等于阈值p,那么它被认为是噪声IMF;否则,它被认为是IMF的有用信号。当阈值p设置为0.1,可以有效区分噪声IMF和有用信号IMF。
步骤3:使用 WTD 方法对噪声 IMF 进行降噪。
步骤4:将去噪后的噪声IMF与剩余的有用信号IMF相结合以进行重建。
由于原始信号中含有部分中高频分量,原有的ICEEMDAN算法的去噪模式会直接舍弃显著含噪的IMF分量,这一行为在降低噪声的同时也带来了有用信号的损失。本文方法采用改进原有的ICEEMDAN算法,通过样本熵判据的方法判断含噪IMF分量,将原本被舍弃的IMF分量进行小波阈值去噪,从而提取其中的有用信号,在提高信噪比的同时减少有用信号的损失。
2 实验与结果分析
2.1 实验设置
实验使用的Φ-OTDR结构如图1 所示。本实验采用具有窄线宽的相干光源,其中,心波长为1 550 nm,线宽为3 kHz。通过应用声光调制器对光信号进行调制,产生脉冲光。掺铒光放大器对脉冲光进行放大,注入环形器。然后,光信号通过单模待测传感光纤传播,传感光纤产生的后向瑞利散射光经由环形器的另一接口输出,进入密集波分复用器。经过耦合后光信号被光电探测器检测,然后通过具有250 MS/s采样率的高速采集卡进行处理,最终输出到主机终端进行数据处理。
实验中,将用于音频播放的音响放置于传感光纤2 km处,通过主机控制采集卡采集并获得实验数据。对实验数据进行预处理,裁剪信号使得采集信号和原始信号对齐,对两个信号进行去直流和归一化,结果如图2 所示。
使用信噪比(SNR, Signal to Noise Ratio)、均方误差(MSE, Mean Squared Error)、波形互相关系数(CC, Cross Correlation)定量比较六种降噪算法的性能。SNR可以反映降噪算法的去噪效果,MSE和CC可以反映去噪后信号和真实信号之间的接近程度。信噪比、均方误差和波形互相关系数的计算公式如下:
SNR(dB)=10log10 =10log10
MSE= (Yi- )2
Rxy(m)= x(k)y(k-m)
其中,Yi是原始信号幅值, 是去噪信号的幅值,x(k)是原始信号序列,y(k-m)是去噪信号序列,m是两个序列的移位长度,N是信号序列的长度。
2.2 实验结果分析
由于ICEEMDAN是一种自适应分解算法,因此,无须针对不同信号提前设置分解次数。这一优点可以减少参数设置不当对去噪效果造成的负面影响。与常见的用于Φ-OTDR去噪的奇异值分解和小波阈值去噪相比,本文算法具有显著的优势,无须设置分解层数使其可以针对不同类型的信号进行降噪处理。设定噪声方差为0.2,添加噪声1次,最大分解层数为500。ICEEMDAN对实测信号的分解结果如图3 所示。
对得到的每个IMF分量计算样本熵,重构维数为1,重构阈值大小为0.25倍的信号标准差,将求得的样本熵与设定的阈值0.1进行比较,如表1 所示。依据样本熵阈值判据,IMF1和IMF2的样本熵超过了设定阈值,说明IMF1和IMF2中存在噪声信号。
表1 各IMF分量样本熵Tab.1 Sample entropy of each IMF component |
| 信号 | 样本熵 |
|---|---|
| IMF1 | 0.446 |
| IMF2 | 0.118 |
| IMF3 | 0.095 |
| IMF4 | 0.017 |
| IMF5 | 0.001 |
| IMF6 | 0 |
| IMF7 | 0 |
| 设定阈值 | 0.1 |
从图4 可以看出,IMF1和IMF2中存在中高频的噪声频谱,印证了样本熵阈值判据的有效性。根据阈值选择,选择IMF1和IMF2进行小波软阈值去噪处理,然后将它们与其他分量进行合并,重新构建去噪信号。
根据所提方法的步骤对采集的实际信号进行去噪,并将所提方法与小波阈值去噪方法、谱减法、EMD法、EMD-WTD法和ICEEMDAN法进行比较。小波阈值法采用软阈值去噪,选择db4小波作为分解的小波基函数,分解尺度为5。保留更多信息。谱减法相比于其他信号多截取前0.5 s的背景噪声信号,用以模拟稳态噪声分布。EMD方法和EMD-WTD方法的模态数与ICEEMDAN方法模态一致,判断所有EMD类算法得到的IMF分量样本熵,高于阈值的分量选择去除,余下的IMF分量参与信号的重建。通过本文所提降噪方法去噪后的信号结果如图5 所示。
从图5 可以看出,本文采用的算法对信号能起到一定的还原效果,但在稳态部分仍存在噪声的轮廓,呈现明显的波动。
表2 各算法实测信号去噪后的SNR、MSE、CC值Tab.2 SNR, MSE, and CC values of measured signals denoised by each algorithm |
| 算法 | SNR | MSE | CC |
|---|---|---|---|
| 原始数据 | 5.24 | 0.022 0 | 0.883 |
| 小波阈值 | 7.65 | 0.008 3 | 0.912 |
| 谱减法 | 5.94 | 0.014 8 | 0.897 |
| EMD | 7.64 | 0.104 7 | 0.924 |
| EMD+WTD | 8.85 | 0.096 2 | 0.931 |
| ICEEMDAN | 9.21 | 0.008 3 | 0.928 |
| 本文算法 | 10.58 | 0.007 2 | 0.940 |
以上分析表明,本文提出的算法不仅具有优异的降噪性能,而且在去噪信号的重建过程中保持了有用信号的完整性。但是对硬件系统造成的稳态波形变化,降噪算法很难将其还原至原始状态。在实验中计算了各方法的计算时间,虽然计算时间变化趋势和信噪比相反,但在实际使用过程中,可以通过提升计算机的性能和并行计算来减少运行时间。Φ-OTDR在实际应用中通常监测较长时间,采样间隔的时间较大。因此,所提方法对于实际应用是可行的。
3 结束语
本文提出了用于Φ-OTDR信号降噪的ICEEMDAN-SE-WTD方法。通过在实测信号上实验验证了所提方法在提高信号质量、增大信噪比、保持信号原始信息方面具有显著的优势。与常见的Φ-OTDR信号降噪算法相比,本文所提算法各项参数均有领先。相比于原始采集信号,信噪比提高了5.34 dB,均方误差减少了0.0148,为原来的32.7%,波形互相关系数提高了5.7%,与本文实验的其他降噪算法相比,改进方法具有显著的信噪比提升,可以有效提高Φ-OTDR的性能,是一种实用的降噪方法,能够提高分布式光纤声波传感系统的测量效果,为后续声波信号的分析和分类提供了条件。