一种基于极值-留数的高背景噪声测试信号降噪方法研究

李 颖， 卢洪超， 周 琳， 陈文文， 齐聪山， 刘福顺

(1. 浙江科技学院 中德工程师学院， 杭州 310023； 2. 中国海洋大学 工程学院， 青岛 266100；3. 海洋石油工程(青岛)有限公司， 山东 青岛 266520)

对结构实测振动响应信号进行分析是获取结构动力特性的重要手段。然而，在土木工程、海洋工程领域，测试设备、环境等因素导致测试的结构振动信号中不可避免的包含各种噪声。信号中的噪声不仅使结构模态识别的精度降低，还会产生“噪声模态”[1]。当实测信号中包含可能淹没结构固有模态的高背景噪声时，结构模态参数识别十分困难。为此，在对工程结构进行模态参数识别、结构损伤诊断和健康检测时，为了保证分析结果的正确性，往往需要采取合理的技术手段降低测试噪声的影响[2]。

目前，已经有许多信号消噪技术应用于工程实践。Cadzow[3]提出了基于奇异值分解(Singular Value Decomposition, SVD)的低秩逼近算法，该方法运用结构响应数据构建Hankel矩阵，进行SVD分解，认为较小的奇异值对应噪声，将其置为零并重构Hankel矩阵，从而达到消噪的目的。该方法对信噪比(Signal-Noise Ratio, SNR)较小的信号消噪效果较好，当噪声能量较大时，噪声成分对应的奇异值可能比结构固有成分对应的奇异值更大，从而表现出消噪局限性；另外，当响应数据较大时，导致Hankel矩阵的维数较大，进行SVD分解会消耗大量计算机资源[4]，甚至无法完成计算。Shumway等[5]对时域平均法作了详细阐述，时域平均法用一段数据的平均值代替中心点的数值，进行平均的方式包括直接平均、多项式平均和周期回归平均等方法，并将各种方法应用于时间序列，对时间序列进行光滑。时域平均法具有一定的消噪效果，但各种平均方法的选取缺乏参照，没有对消噪效果进行定量讨论，另外运用该方法还需要时标信息，而且要有足够的数据量。Kim等[6]在频域内对频响函数(Frequency Response Function, FRF)运用小波分解进行消噪，并将该方法运用于2自由度质量-弹簧-阻尼系统，能够有效消除信号中添加的高斯白噪声，但小波消噪算法对小波核函数的选择有较大依赖性，在实际运用中受到一定局限。Wu等[7-8]对经验模态分解(Empirical Mode Decompositon, EMD)的研究表明，白噪声经EMD分解得到的固有模态函数(Intrinsic Mode Function, IMF)分量的统计特征满足正态分布，IMF的能量密度与对应的平均周期的积为常量，这些研究结果为EMD消噪提供了理论基础。对含噪声的信号进行EMD分解，得到频率从高到低的IMF，高频部分受噪声影响较大，低频部分受噪声影响较小，而信号的能量也主要集中在低频段，使用低频段的IMF进行重构完成消噪过程[9]。基于EMD分解的消噪方法对模拟的高信噪比信号具有较好的消噪效果，但该方法可能舍弃部分有用信号的能量，而且无法准确判定选取有效IMF的个数。另外，近年来还有学者运用EMD分解方法对非平稳信号的消噪进行消噪研究，并应用于转子振动信号的去噪，取得了较好效果[10]。

本文将从信号分解的角度出发，提出一种基于极值-留数的信号分解消噪方法，该方法提取代表信号中固有信息的极值和留数，并根据这些信息进行重构，从而达到消噪目的。与传统方法相比，本文方法避免反复进行SVD分解运算，具有较高的计算效率。另外，进行信号分解时，无需积分变换的核函数，便于实施。该方法能够解决实际工程中在遭遇高背景噪声时，其结构固有信息往往被当做环境噪声而剔除的问题。文中将首先应用一质量-弹簧-阻尼模型进行验证，然后应用实测海洋平台加速度响应数据进一步研究新方法在实际工程中的应用效果。

1 信号分解方式

1.1 傅里叶序列

连续函数y(t)的周期为T，其傅里叶分解为

(1)

(2)

对于实测的离散数字信号yk(k=0,2,…,N-1)，其傅里叶分解为

(3)

其中xn=ei(2πn/N)。

1.2 Prony序列

连续信号y(t)分解成复指数的形式为[11]

(4)

式中:p为分解阶次;λn为极值，由于y(t)为实数，λn必须为实数或为共轭复数对。令λn≡-αn+iωn，αn为阻尼系数，ωn为角频率。系数γn为留数，它与λn成对出现，或为实数，或为共轭复数。令γn≡Aneiθn，其中An为各个成分的幅值，θn为对应相位。

等间隔采样离散信号yk，其复指数分解为

(5)

式中:zn=eλnΔt。

2 基于极值-留数的信号降噪方法

对信号进行Prony分解，得到一系列极值和留数，每个极值和对应的留数表示信号中的一个成分。通过式(4)进行叠加，即可用极值和留数表示实测信号。当信号中包含噪声成分时，噪声成分也有对应的极值和留数，本文方法就是通过剔除噪声成分的极值和留数，从而达到消噪的目的。因此本文消噪方法的关键在于准确求解信号所对应的极值和留数。

为了求解式(4)中的极值和留数，引入高阶微分方程

(6)

(7)

式(4)为式(6)的通解，其中λn为式(6)对应特征方程式(7)的根。

对于实测的等间隔采样离散数据yk(k=0,1,…,N-1)，式(6)高阶微分方程变为差分方程

(8)

式中:bn为常系数。令bn=1，式(9)对应的特征方程为

(9)

通过以下步骤求解极值和留数

步骤1将式(8)写成矩阵形式，由于Y和y′已知，可求出向量b，即系数bk(k=0,1,…,p-1)

Yb=-y′

(10)

其中

(11)

(12)

(13)

步骤2将bk代入式(9)中，求得特征值zk(k=1,2,…,p)

步骤3将特征值zk代入式(5)中，即可求得对应的留数γk。

2.1 一阶差分方程求解方法

由于求解方程式(10)时，存在病态问题，很小的误差会使结果产生较大的偏差。为了解决病态问题，本文将高阶差分方程转化为一阶差分方程，即令

x1,k=yk,x2,k=yk+1,…,xp,k=yk+p-1

(14)

则

(15)

式(8)可以写成

xk+1=Fxk

(16)

其中

(17)

式(16)表明，其特征值即是方程式(7)的根。故可由响应数据yk构建Hankel矩阵

(18)

式中:ξ、η分别为Hankel矩阵的行数和列数，为了提高计算效率，ξ和η的值都取最接近于数据序列长度一半的整数。当k分别取0和1时，对H(0)进行奇异值分解，得

(19)

(20)

由式(21)得到系统矩阵A

(21)

对A进行特征值分析，得到特征值zn，则响应信号的极值λn=ln(zn)/Δt，将λn代入式(5)即可求得对应的留数γn。

2.2 噪声信号分离与剔除

海洋结构实测振动数据y(n)可以看作由两部分组成[12]

y(n)=yreal(n)+ynoise(n)

(22)

式中:yreal(n)为结构的真实响应信号;ynoise(n)为实测信号中的噪声信号。

对实测信号进行极值-留数分解，得到极值和留数

(23)

(24)

去除噪声对应的极值和留数后，得到实真响应信号的极值和留数

(25)

(26)

根据式(5)重构响应信号，得到真实信号

(27)

极值-留数分解消噪算法步骤如下：

步骤1对海洋结构实测振动信号y(n)进行傅里叶变换，得到响应信号的频谱；

步骤2根据频谱找到峰值对应的频率，并以此频率为中心频率设置频率窗口；

步骤3对实测信号y(n)进行极值-留数分解，取极值的绝对值，将复频率转化为实频率，选出步骤二频率窗口中的实频率对应的极值和留数；

步骤4通过式(5)，将步骤三得到极值和留数重构成时域信号，完成消噪过程。

3 数值算例

为了验证本文方法的有效性，构建一个6自由度质量-弹簧-阻尼模型，如图1所示。系统中各质量块的质量为mi=100 kg，刚度为ki=5×107N/m，阻尼系数为ci=50 N·s/m，各质量块的位移分别为xi，其中i=1,2,…,6。

图1 6自由度质量-弹簧-阻尼模型

根据上述模型的系统状态空间方程，利用MATLAB软件的impulse函数生成离散脉冲响应数据。脉冲响应的采样频率设为500 Hz，信号序列包含2 048个点。对于不同的激励点，整个6自由度系统会产生36个响应序列，为了便于分析，选取在质量m1处输入、m1处输出的脉冲响应数据进行研究，其时程如图2所示。

图2 分析数据的时程图

为了研究消噪效果，在信号中添加白噪声模拟实测信号中的噪声。使用SNR来描述噪声水平，SNR为信号功率Ps与噪声功率Pn之比，即SNR=10lg(Ps/Pn)。对分析数据分别添加不同信噪比的白噪声，模拟实测信号中的不同强度噪声。

对含信噪比为40 dB噪声的信号和不含噪声的信号进行傅里叶分析，得到信号的频谱，它们的时程和频谱对比如图3所示。

图3 不含噪声信号与含信噪比40 dB噪声的信号时程与频谱对比

Fig.3 The time domain and frequency spectrum comparison between clear signal and signal containing SNR 40 dB noise

根据频谱的峰值，分别设置频率范围(26.5，27.5)、(79，80)、(127，128)、(168，169)、(198.5，199.5)和(218，219)，6个峰值分别包含于6个频率范围之内。将信噪比为40 dB的信号代入式(18)构建Hankel矩阵，根据式(21)得到系统矩阵A，进而计算极值和留数。对极值取模获得实频率，根据实频率选取各频率范围内对应的极值和留数，代入式(5)中进行重构，得到消噪后信号。消噪后信号与不含噪声信号的时程和频谱对比,如图4所示。结果表明本文方法能有效消除信噪比为40 dB的噪声。

为了进一步验证本文方法消噪效果，将该方法应用于噪声水平较高的信噪比为10 dB的信号，消噪前后的时程和频谱如图5和图6所示。消噪前第六阶频率几乎淹没于噪声中，消噪后得到的频谱与不含噪声的频谱重合较好，表明本文方法对高噪声水平的信号也有较好的消噪效果。

图4 不含噪声信号与含信噪比40 dB噪声的信号消噪后的时程与频谱对比

Fig.4 The time domain and frequency spectrum comparison between clear signal and signals eliminated noise

图5 不含噪声信号与含信噪比10 dB噪声的信号时程与频谱对比

Fig.5 The time domain and frequency spectrum comparison between clear signal and signal containing SNR 10 dB noise

图6 不含噪声信号与含信噪比10 dB噪声的信号消噪后的时程与频谱对比

Fig.6 The time domain and frequency spectrum comparison between clear signal and signals eliminated noise

为了定量分析本文方法的消噪效果，引入噪声抑制率(Noise Suppression Ration, NSR)、信号失真率(Signal Distortion Ratio, SDR)[13]和决定系数(Coefficient of Determination)三个指标。

NSR=

(28)

(29)

式中：s(n)为不含噪声信号；f(n)为添加噪声后f′(n)为消噪后信号。

(30)

对于信噪比为40 dB的信号，本文消噪方法消噪后的噪声抑制率为0.797 8，信号失真率为0.002 008 4，决定系数为0.999 96。运用本文方法分别对含信噪比为30 dB、20 dB和10 dB噪声的信号进行消噪，其噪声抑制率和信号失真率列于表1中，消噪前后的决定系数如图7所示。对比文献[13]，其提出的消噪方法的噪声抑制率为0.872 7，信号失真率为0.174 4。本文方法消噪后得到不同噪声水平的噪声抑制率都约为0.8，信号失真率随着噪声的增大而逐渐增大，但总都维持在比较低的水平；消噪后的决定系数R2都大于0.995，随着噪声水平的提高，决定系数R2逐渐减小。三个指标的分析结果表明本文的消噪方法消噪效果较好。

表1 不同信噪比信号的噪声抑制率(NSR)和信号失真率(SDR)

4 实测数据应用

海洋结构物实测振动响应信号往往含有较高水平的背景噪声，尤其环境激励下其结构固有模态往往淹没于噪声之中，导致结构动力特性识别困难。为研究方法的有效性，选用冰击激励作用下的某海洋石油平台(见图8)实测数据进行研究。结构振动测试时，使用三向加速度传感器采集其加速度响应信号,传感器安装位置如图9所示。采样频率为200 Hz。选取主振动方向300～320 s的实测数据进行分析，其加速度响应时程如图10所示。对加速度响应数据进行傅里叶分析，得到频谱如图11所示。由于海洋平台实际环境激励情况比较复杂，激励能量较弱，得到频谱非常复杂，包含大量噪声成分。这些噪声使响应信号的进一步分析面临很大困难。

图7 决定系数随信噪比变化

图8 海洋石油平台

图9 加速度传感器安装位置

1994年，该平台所有单位对其进行了振动检测，应用传统的频谱分析技术，得到该石油平台的前2阶频率，分别为0.9 Hz、1.15 Hz与1.175 Hz。考虑到该平台至今仍在服役，其结构固有特性相对稳定。故将应用本文方法对图10所示信号进行消噪处理。设置频率窗口为0.8～1.3 Hz，并采用获取图5同样的技术流程，可以得到消噪后的信号，其频谱如图12 所示。从图12可知，实测原始信号中小于0.8 Hz、大于1.3 Hz的频率成分被成功剔除，其0.8～1.3 Hz范围内有2个明显的峰值，峰值处对应的频率分别约为0.85 Hz和0.95 Hz。为进一步验证消噪信号的正确性，对其应用ERA方法进行模态参数识别，其结果为0.916 19 Hz和1.036 7 Hz，相对1994年结果计算精度得到大幅提高。

图10 加速度响应信号

图11 频谱图

图12 消噪后频谱图

5 结 论

本文从信号分解的角度提出了一种基于极值-留数的信号分解消噪方法，解决了实际工程中在遭遇高背景噪声其结构固有信息往往被当作环境噪声而剔除的问题。选用的6自由度质量-弹簧-阻尼系统已经证实：当信噪比(SNR)分别为40 dB和10 dB时，消噪后信号与不含噪声信号的时程图和频谱图几乎重合，表明该方法能有效消除信号中强噪声的干扰。计算获得的不同信噪比下的噪声抑制率(NSR)和信号失真率(SDR)，NSR约为0.8，SDR维持在较低水平，当噪声能量变化较大时，使用本文方法消噪得到的NSR变化较小，SDR较大。不同信噪比下，消噪后信号与不含噪声信号的决定系数都大于0.995，表明消噪效果较好。为研究方法的有效性，选用了冰击激励作用下的某海洋石油平台实测数据进行研究，获得了预期频率区间的频率成分，基于此进行模态参数识别，可以大大提高模态参数识别精度，具有潜在的工程应用价值。