基于经验小波变换-噪声辅助分析的桥梁信号降噪方法

2022-11-21罗烨钶陈永高李升才

振动与冲击 2022年21期

罗烨钶，陈永高，李升才

(1.浙江工业职业技术学院建筑工程学院，浙江绍兴 312000； 2.华侨大学土木工程学院，福建泉州 362000)

桥梁健康监测系统通过监测结构关键部位的荷载与环境、结构响应以实现结构状态的评估，从而辅助桥梁运维、确保交通运输安全运行。实际工程中的桥梁运营环境极为复杂，监测系统往往会受到噪声干扰，使得能够反映结构特征信息的监测信号淹没于噪声中，给桥梁结构状态评估带来极大影响。因此，在实桥监测信号分析前需通过合理的方法对采集信号开展降噪处理，最大程度地降低噪声干扰，为桥梁结构状态评估提供合理可靠的数据基础。

目前，桥梁信号降噪中应用较多的方法有经验模态分解(empirical mode decomposition，EMD)、小波变换(wavelet transform，WT)和奇异值分解(singular value decomposition，SVD)。EMD和WT通过对信号进行分解、筛选、重构实现降噪，SVD通过选定合适的奇异值阶次进行信号重构实现降噪。路华丽[1]基于EMD滤波及切比雪夫滤波对某单塔自锚式悬索桥变形监测数据进行了有效降噪；颜潮勇[2]基于WT提出了分层、分段降噪方法，并通过苏通长江公路大桥GPS(global positioning system)监测数据进行了验证；熊春宝等[3]将完备集合经验模态分解(complete ensemble empirical mode decomposition，CEEMDAN)与WT结合，对某大跨度斜拉桥主梁竖向位移数据进行了降噪；董林鹭等[4]提出了一种基于局部均值分解和 SVD 的降噪方法，通过 SVD 对局部均值分解法分解筛选后的分量信号进行降噪，实现了对微震信号中高频噪声的去除；叶锡钧[5]提出了基于GA-SVD(genetic algorithm and singular value decomposition)的自适应降噪算法，并以某斜拉桥和某连续梁拱桥的振动信号为研究对象进行了验证。

桥梁结构体型质量巨大、运营环境复杂，其监测信号具有信噪比小及细节成份丰富的特点。在桥梁振动监测信号的降噪中，EMD易产生模态混叠，端点效应严重，且需要依托外部滤波算法实现降噪，WT需要考虑桥梁结构特征对信号的影响来选择小波基及阈值函数，SVD的奇异值阶次确定大多依赖人工判别，3种方法体系均具有一定的局限性。

经验小波变换(empirical wavelet transform，EWT)是由法国学者Gilles于2013年提出的一种新的时频域信号分析方法[6]。EWT能够通过频谱分割算法实现频带划分，并基于Meyer小波基构造正交小波滤波器组，拥有完备理论基础同时兼备较高的计算效率，在机械故障诊断领域得到了成熟应用[7]。相比EMD,WT,SVD等算法，EWT可通过频谱分割算法实现频带精确划分，并基于Meyer小波基进行本征模态函数(intrinsic mode function，IMF)滤波，针对桥梁振动监测信号具有较好的适应度。本文将EWT引入桥梁监测信号降噪中，并结合工程实际应用需求，提出一种基于噪声辅助分析理论的桥梁结构响应自适应降噪方法，以实现对复杂环境激励下桥梁监测信号的降噪。先通过仿真信号对所提方法应用效果进行测试，后以某斜拉桥加速度监测信号为例，验证了所提方法在实桥监测信号降噪分析中的有效性。

1 EWT原理及降噪分析

1.1 EWT算法

EWT步骤如下。

图1 EWT频谱分割示意图Fig.1 EWT spectrum segmentation diagram

步骤2基于频带边界坐标构造滤波器组。根据Little wood-Paley理论构建小波的方法，EWT中Meyer小波基的尺度函数φn(t)及小波函数Ψn(t)在频域内表达详见式(1)及式(2)。

(1)

(2)

τn=γωn

(3)

(4)

β(ω)为满足式(5)的任意函数

(5)

步骤3通过小波基函数滤波得到分量信号。根据WT的方法将EWT的细节系数和近似系数分别定为式(6)和式(7)。细节系数和近似系数分别代表了信号的高频成分和低频成分。

(6)

(7)

通过分解所得信号分量为

(8)

(9)

(10)

1.2 EWT降噪分析

随着EWT的进一步应用，已有诸多学者尝试将EWT用于机械领域的信号降噪中[8]。EWT实现降噪主要有3类途径：① 提升频谱分割精度，使不同频率成分的特征分量彼此之间得到分离，避免模态混叠，同时使噪声分量与特征分量分离，避免经Meyer小波基滤波所得IMF分量被噪声污染[9]；② 提高特征IMF分量的识别精度，通过构造对特征或噪声分量具有一定敏感度的筛分指标，进而判定IMF中的噪声分量或无效分量，对其筛分后重构达到降噪目的[10]；③ 对IMF分量进行滤波，如基于WT等方法对IMF分量滤波重构实现降噪[11]。

根据工程应用需求不同，现有基于EWT的降噪案例大多基于上述3类途径进行组合，但受制于EWT算法的特点，上述3种途径均具有一定短板。如：途径①的效果极度依赖频谱分割方法的分割精度，若频谱分割不当，会严重影响后续的滤波质量，目前大多在频谱分割前对信号进行谱分析或构造拟合曲线，进而突出信号中的特征谱峰、抑制噪声谱峰，之后对已平滑的频谱曲线进行“伪”频谱分割，得到的频带坐标又返至原始信号进行“真”频谱分割，而此种方法取决于平铺曲线的平滑效果，不同工程场景中无法统一适用；途径②的指标建立在人为选择的基础上，与EWT中固定的频谱分割算法不同，指标选择很大程度依赖研究人员的专业判断能力，不同的指标对信号的敏感程度不一，如何基于信号特征进行自适应指标构建尚无定论；途径③将EWT滤波效果依赖于其他算法，构建基于EWT和所选算法的组合降噪框架，该方法虽能够提升最终的降噪效果，但增加了需要调整的参数，在提高应用难度的同时也降低了算法的自适应性。

2 改进EWT降噪算法

2.1 染噪信号的构造

为抑制EMD分解时的模态混叠和端点效应、降低信号噪声对分解过程的干扰，Huang基于噪声辅助分析理论(noise-assisted data analysis，NADA)对EMD进行了改进。NADA理论主要体现为：向原始信号添加高斯白噪声，基于EMD获得每次加噪后分解所得IMF分量，再对各组IMF分量集成平均，进而抑制IMF分量中噪声的影响，重构后即可实现原始信号降噪。

根据Huang所提噪声辅助分析的思想，向原始信号x(t)中添加标准差比例为σ的高斯白噪声，根据文献[12]及EWT算法特点，将σ拟为0.2，将集成次数考虑为迭代次数Ne拟为100。当σ和Nei(i=1,2,…,100)确定时，构造染噪信号为

Xi(t)=X(t)+noise·σ

(11)

式中：noise为高斯白噪声;Xi(t)为σ下第Nei次集成所得染噪信号。

2.2 有效频带的选取

实桥测试信号在环境激励下的傅里叶谱往往会被噪声“波浪”影响，这种影响体现在整个频域内。随着激励信息的复杂化及剧烈化，“波浪”会越来越高，直至“淹没”代表结构特征的幅值。在频域内，噪声分布较为均匀且随机，代表结构固有信息的特征分量分布则较为稳定，且在大多工程案例中，特征分量的频域幅值往往比噪声分量高很多。

根据特征分量能量高且稳定、噪声能量低且随机的分布规律，本节提出一种基于能量准则的频带划分方法。首先，基于EWT中的Scalespace-Otsu方法对Xi(t)频带进行划分，得到频带坐标集合W={w1,w2,…,wj}(j∈N+)，根据Gilles所定规则，将该W标准化于[0,π]内，即w1=0,wj=π。然后，提取j-1个频带内的主峰幅值，得到主峰幅值集合A={a1,a2,…,aj}，此处为避免求取频带能量时受到大幅值噪声影响，故以主峰幅值为判断依据。提取最大主峰幅值的频带区间[wm,wm+1](1≤m≤j-1)，利用Meyer小波滤波器组滤波得到FIM,m。定义1阶染噪信号为

(12)

按照以上分析流程，在对大多数特征分量所在区间进行定位并滤波得到IMF分量后，残量频谱形态主要由两类构成，即低幅值含噪特征分量与噪声。针对低幅值含噪特征分量，将以误差指标进行判断；针对噪声分量，理论上当频域范围内没有明显的高幅值谱峰时，证明特征分量已近似提取完毕，此时无需再进行EWT滤波。然而，噪声的幅值特征是低且随机的，当频域范围内皆为噪声分量时，频谱会呈现出整体幅值较低、所有极大值点的幅值较为接近的多峰形态，本文将其定义为“伪染噪信号”。

2.3 分解终止阈值分析

针对“伪染噪信号”的判断标准，本节提出一种基于概率密度函数(probability density function，PDF)曲线的自适应判定方法。原始含噪信号傅里叶谱中，含噪特征信号对应的高幅值区域占比较少，噪声信号对应的低幅值“波浪”区域占比较多，通过对比可知，两者的分布规律具有较大的差异性，对两者频数直方图的PDF曲线进行峰值位置估计，若峰值位置达到某临界状态，则可判定当前信号属于完全噪声信号。以某含有4个特征分量及1个噪声信号的含噪信号为例，X1为原始含噪信号，X2～X5分别为从中去除1～4个特征分量的余量信号。X1～X5的傅里叶谱及对应的PDF曲线如图2所示。

图2 X1～X5的频域统计规律Fig.2 Frequency domain statistics of X1～X5

由图2可知，随着信号特征频点减少，其PDF曲线峰值将逐渐右移，当信号不含特征频点时，PDF曲线的峰值位置将会稳定于某范围(见图2(b)X5)。经过大量试验，确定该范围中心稳定于全局分布区间的33%左右，即以33%为临界阈值。

针对EWT的算法特点，当傅里叶谱某一频带经Meyer小波基滤波后，其残量的傅里叶谱中该区间的幅值曲线将趋近于0，同时呈现出平滑特征，这种与噪声规律不同的分布形式会使PDF曲线估计产生偏差，因此制定规则如下：

(1) 在单次迭代过程中，第2次及后续分割前需将之前分割得到的能量最大的频带剔除，剔除后方可进行PDF曲线估计。

(2) 在单次迭代过程中，当某次分割得到的能量最大区间[wn,wn+1]的频带边界之一wn+1在上一步骤的能量最大区间[wm,wm+1]内时，将wn+1替换为wm。

2.4 IMF分量集成

假设经过p次循环分解后PDF曲线到达临界阈值，此时可得到p个频带区间及p个IMF分量，通过p个频带区间又可得p个主峰幅值和2p个频带边界坐标，同时基于随机减量法又可得p个IMF分量的频率值。基于噪声辅助分析思想，对各IMF分量设定集成规则：首先，将各IMF分量按照识别的频率大小从低频到高频排序；其次，为避免集成过程中不同频段IMF的混叠，定义集成系数Th为

Th=

(13)

第2次迭代结束后可对两组IMF分量进行第1次集成，小于Th的两个IMF分量叠加后取平均，对应的两组特征信息(w,f,a)经同样平均处理后作为新IMF分量的特征信息，大于Th的两个IMF分量及对应的特征信息则予以保留，第1、第2次的集成结果将与第3次迭代结果进行第2次集成，后续迭代以此类推。

2.5 有效IMF分量筛分

由于本方法在集成过程中会保留大于集成阈值Th的IMF分量，故当100次迭代结束后，必然会存在几乎没有被集成平均或平均次数极少的IMF分量。基于噪声辅助分析理论可知，该部分IMF分量中仍存在较强的噪声影响，故应筛分后剔除。针对有效IMF分量的筛分，本节提出一种筛分系数，将能量密度与平均周期乘积ET和JS(Jensen-Shannon)散度结合，通过信号自身噪声特性及与原始信号间的相关性来进行判别。

Wu等[13]基于能量密度与平均周期的乘积ET对噪声IMF分量进行了有效判定。两者定义见式(14)和式(15)

(14)

(15)

式中：E为能量密度;T为平均周期;N为IMF分量的数据长度;Nzc为IMF的过零点个数。当IMF分量特性越接近特征响应，ET值就越大，反之越接近噪声信号则ET值越小。

JS散度是在KL(Kullback-Leibler)散度基础上改进的一种衡量相似性的指标，能够避免KL散度在对两种不同的分布计算时产生的不对称性，可较为合理地反映两种分布的相似程度[14]。以分布P和分布Q为例，两者JS散度的定义见式(16)

(16)

式中,KL(·)为KL散度。JS散度越大，即证明P和Q的分布越接近，在信号分析中，若某IMF分量与原始信号的JS散度越大，即代表该信号蕴含的特征信息较多，反之则越少。

由于ET指标和JS散度具有不同的分布特征，故通过标准化方法对两者处理后进行线性平均，形成筛分系数Ts。经过大量试验发现，当Ts>0.7时，可判定该IMF分量为特征分量，反之则予以筛除。对最终集成所得各IMF分量筛分后再进行特征信号的重构，即可得到最终的降噪信号。本文所提的基于噪声辅助分析的改进EWT方法降噪具体流程,如图3所示。

图3 改进EWT方法降噪流程图Fig.3 Process of improved EWT noise reduction

3 仿真信号试验

3.1 降噪评价指标

为对改进EWT方法的降噪效果进行定量分析，引入信噪比(signal-noise ratio，SNR)、均方根误差(root mean square error，RMSE)、噪声抑制率(noise suppression ration，NSR)等3个指标，其定义如下

(17)

(18)

(19)

式中：s(i)为初始纯净信号;f′(i)为染噪信号;f(i)为对染噪信号降噪后的信号;N为采样点数。SNR和NSR随降噪效果的提升而增大，RMSE反之。

3.2 仿真信号构造

构造正弦信号x1及x4、余弦信号x2、调幅调频信号x3、自由衰减信号x5，并添加模拟高斯环境噪声的噪声信号noise构成含噪仿真合成信号，见式(20)

(20)

将采样频率拟为10 Hz，采样时长2 min，为体现各信号细节，选择0～30 s的时域波形进行展示，如图4所示。合成信号加噪前后傅里叶谱如图5和图6所示。

(a)

(b)

(c)

(d)

(e)

(f)图4 信号构成图Fig.4 Signal composition diagram

(a) 时程曲线

(b) 傅里叶谱图5 纯净模拟信号Fig.5 Pure simulation signal

(a) 时程曲线

(b) 傅里叶谱图6 加噪后模拟信号Fig.6 Simulation signal after adding noise

3.3 改进EWT分析及评价

当迭代次数i=1时，染噪信号经改进EWT方法的分析流程如图7和图8所示，循环编码为Ci(i=1，…，N)。6次循环过程对应的峰值位置比例,如表1所示。

表1 前两组IMF分量的集成系数Tab.1 Integration coefficients of the first two IMF

(a) C1-分割结果

(b) C2-分割结果

(d) C4-分割结果

(e) C5-分割结果

(f) C6-分割结果图7 染噪信号首次迭代分割规律Fig.7 The first iterative segmentation rule

(a) C1-PDF曲线估计

(b) C2-PDF曲线估计

(d) C4-PDF曲线估计

(e) C5-PDF曲线估计

(f) C6-PDF曲线估计图8 残量概率密度曲线Fig.8 Probability density curves of residual

当循环次数C为6时，此时PDF曲线峰值位置比例为0.363，超出了33%的阈值限界，故以前5次循环所得IMF分量为有效IMF分量。此外，C3的位置比例较C1,C2偏小，PDF曲线的峰值位置比例从C4开始恢复正常规律，究其原因，可能是C2中的[b8,b9]频带受噪声影响较大，故在滤除该区间后傅里叶谱的整体信噪比有所增强，故对应的峰值位置比例变小。针对实际工程而言，当循环次数达到一定程度时，残量频域形态必然会达到接近于噪声的临界状态，基于EWT算法自身特性，无论频带受到噪声的影响强弱与否，改进EWT方法均能够做出准确的临界状态判定。

第1次迭代共得到5阶IMF分量，按同样步骤完成第2次迭代，根据集成规则对两次迭代过程所得IMF分量进行集成平均，两组IMF分量对应的集成系数如表2所示。由表2可知，第1次迭代所得5阶IMF分量可分别与第2次迭代的IMF21，IMF22，IMF23，IMF25，IMF26等分量进行集成，两组数据之间的集成系数均大于0.8，满足集成条件。

表2 前两组IMF分量的集成系数Tab.2 Integration coefficient of the first two IMF

按相同流程迭代至第100次，此时共得到12个IMF分量，按照2.5节所述流程，计算各IMF分量的能量密度与平均周期乘积ET以及JS散度，对两者进行标准化后再通过线性平均形成筛分系数，如表3所示。对满足筛分条件的各IMF分量重构得到降噪信号。

表3 最终迭代所得各IMF分量筛分指标Tab.3 Screening indexes of IMF obtained by final iteration

为量化本文所提改进EWT方法的降噪效果，同时对比该方法的应用性能，分别采用熊春宝等研究中的CEEMDAN-WT降噪、董林鹭等研究中的SVD降噪、改进EWT降噪3种方法对模拟含噪信号进行降噪。其中：CEEMDAN-WT按熊春宝等的研究确定IMF筛分阈值为0.5，选用db6小波基、软阈值；SVD按董林鹭等的研究采用加权能量贡献率确定奇异值有效阶数，阈值为0.1%。

CEEMDAN-WT和改进EWT分解所得IMF分量，如图9所示。由图9可知，改进EWT的IMF分量要比前者更少，且各分量形式与原始分量信号接近，表明经过改进EWT处理后各特征IMF分量可有效保留并分离，且具有较好的分解精度，CEEMDAN-WT中IMF1～IMF4与原始分量虽较为接近，但未能有效反映原始分量的周期规律，且在分析过程中易混入无效IMF分量。对各方法重构后信号进行对比，时程图、重构误差、傅里叶谱，如图10～图12所示。由图10及图11可知，SVD与改进EWT的重构信号与纯净合成信号相近，CEEMDAN-WT效果则相对较差，且前两者的重构误差均比后者要低，改进EWT的最大重构误差为1.358‰，相比SVD的最大重构误差4.026‰更低。观察图12可知，CEEMDAN-WT重构信号傅里叶谱与原始含噪信号傅里叶谱接近，其中高频部分降噪效果较低频部分更明显，但傅里叶谱全局降噪效果较差。SVD与改进EWT的降噪效果更好，其中，SVD在降噪过程中过滤了x1特征频率，x3虽然能观察到频率主峰，但其调幅调频特性并未在重构信号傅里叶谱中体现。改进EWT相比SVD降噪效果更好，x3的调幅调频特征得到保留，在降噪过程中x4受到x5影响较小，未出现SVD降噪后出现的x4特征频率幅值减小的情况。经综合分析，改进EWT的全局降噪效果相比CEEMDAN-WT及SVD更好，且相比后两者能够更加完整地保留分量信号的频域特征。

图9 模拟信号各方法分解结果Fig.9 Decomposition results of simulation signals

(a) CEEMDAN-WT

(b) SVD

(a) CEEMDAN-WT

(b) SVD

计算重构信号的降噪指标，如表4所示。经前述可知，SNR与NSR越大、RMSE越小，表明降噪效果越好。由表4可知，经改进EWT降噪后重构信号的SNR，NSR，RMSE分别为25.254 dB,1.653×10-3,0.874。观察可得，SVD的三项降噪指标与改进EWT在数量级保持一致，但仍有一定差距，如SVD的RMSE为4.847×10-3，接近改进EWT对应指标的3倍，而SNR与NSR则较为接近。CEEMDAN-WT的降噪效果最差，其SNR为9.254 dB，是3种方法中的最小值，接近SVD的1/2、改进EWT的1/3，对应的NSR保持同样的数值规律。此外，由于CEEMDAN-WT未能对信号中的高频噪声分量进行有效降噪，导致其RMSE值达到2.986。

表4 重构信号降噪指标Tab.4 Noise reduction index of reconstructed signal

经综合对比，改进EWT相比CEEMDAN-WT、SVD具有更好的降噪效果，该方法能够在有效提取特征分量的同时保持较低重构误差。

4 工程验证

为验证改进EWT方法在实桥监测数据分析中具有良好的降噪效果，选用某斜拉桥为研究对象，基于Hilbert-Huang谱对降噪效果进行评价。

4.1 工程概况

该斜拉桥主桥由464 m中跨和两侧对称布置的216.5 m边跨组成，主桥全长897 m，属双塔双索面半漂浮体系。大桥跨度大、高程相对较高、车流量多，风荷载和车辆荷载导致桥面振感明显，该桥加速度传感器布置形式，如图13所示。

图13 加速度传感器布置图(m)Fig.13 Acceleration sensor layout (m)

本文选取该桥中跨跨中截面的Z-6#竖向加速度信号为研究对象，采集所用传感器为941B竖向拾振器，采样频率为50 Hz，Z-6#加速度时程曲线及傅里叶谱详，如图14所示。由图14可知，跨中截面特征频率主要集中在[5,15]Hz，但受环境噪声影响，特征频率区间被噪声淹没，出现毛刺现象，无法辨别各特征频率信息。

(a) 时程图

(b) 傅里叶谱图14 Z-6#监测数据Fig.14 Z-6# monitoring data

4.2 降噪效果验证

分别采用CEEMDAN-WT、SVD、改进EWT方法对Z-6#测点的采集信号进行降噪，CEEMDAN-WT与改进EWT方法分解所得IMF分量，如图15所示;对应的时频曲线,如图16所示;3种降噪方法重构信号的傅里叶谱，如图17所示。

图15 实桥信号各方法分解结果Fig.15 Decomposition results of bridge signals

图16 改进EWT分解后的 Hilbert-Huang 谱Fig.16 Hilbert-Huang spectrum of Bridge signal

(a) CEEMDAN-WT

(b) SVD

由图15可知，在时域内，CEEMDAN-WT分解得到8个IMF分量，经分析可知IMF1～IMF4为高于15 Hz高频分量，IMF5,IMF6代表了[5,15]Hz的特征分量，IMF7,IMF8为接近5Hz的低频分量，改进EWT分解得到6个IMF分量，经分析得IMF1～IMF6均处于[5,15]Hz。在分解数量上改进EWT更具优势，且通过与原始信号的傅里叶谱对比，CEEMDAN-WT仅有IMF5,IMF6为有效信号分量，而改进EWT的各IMF分量均处于特征频带区间，表明改进EWT能够准确定位信号特征分量所处频带范围。由图15可知，在时频域内，改进EWT所得IMF分量的时频曲线更加清晰和连续，反映出改进EWT方法具备良好的降噪性能。通过图16中各时频曲线对应的频率区间可知，改进EWT的IMF分量特征频率主要集中在[5,15]Hz，与图1中原始信号傅里叶谱基本规律一致，即经过频带提取、集成平均、分量筛选后，改进EWT能够在保持分解数量的优势的同时保留信号特征信息。由图17可知，在频域内，SVD与改进EWT均能够有效抑制环境噪声对信号中特征分量的影响，在保留特征分量频率信息的同时能够避免无效噪声对特征频带的毛刺干扰，CEEMDAN-WT的降噪能力则相对较差。通过观察可知，改进EWT对6.3 Hz主频附近的噪声进行了剔除，而SVD在该位置未能有效降噪，在整体范围内改进EWT比SVD具备更好的降噪效果。