基于奇异值分解的信号消噪技术

1、基于奇异值分解的信号消噪技术摘 要模态参数识别是从结构不同位置的动力响应信号中提取出结构的模态参数，即：从动力测试响应信号数据中确定结构的模态参数(模态振型、固有频率和阻尼比)。每一个结构都有其固有的模态参数，并且如果结构动力特性发生变化了，那么结构的模态参数也将发生相应的变化。显见，结构的模态参数识别是非常重要的，为诊断结构健康状况提供了依据。基于输出的模态参数识别方法利用的信息主要是系统的自由振动信号，要获得自由振动信号首先需获得结构的响应信号。由于环境激励的不充分和噪声等干扰因素的存在，导致信号测试信号不能直接用于参数辨识，需要对信号进行消噪处理。即从大量背景噪声中提取出可用于模态参数辨

2、识的有用信号成分,剔除干扰因素，提取有用信息。此时，信号消噪技术研究变得尤为重要。 本文采用了一种将Hankel矩阵和奇异值分解相结合的消噪方法。该方法首先对测量信号构造的Hankel矩阵进行奇异值分解，再利用测量信号快速傅立叶变换结果中主频率的个数来确定有效秩阶次，接着通过消噪信号的信噪比和均方差大小确定重构矩阵结构，最后通过反对角线平均法得到消噪后的信号数据。通过数值仿真，对不同信号进行定秩和消噪，从结果可以知道这种方法具有较好的消噪效果。关键词：信号消噪；奇异值分解；快速傅立叶变换；信噪比；均方差 A Method for Noise Reduction Based on Singula

3、r Value DecompositionAbstract Accurate estimate of the modal parameters of an offshore structure is crucial to many practical engineering issues, such as finite element (FE) model updating and validation, damage detection, etc. Modal parameter identification method uses the the response signal of st

4、ructure ,but actual response signal often contains a lot of noise, which will affect the accuracy of signal recognition. The test signal de-noising processing is an important step in signal processing. Using Singular Value Decomposition（SVD）of constructed Hankel matrix by measured signal is an effec

5、tive method for eliminating the random noise. The key is to choose the rank of the Hankel matrix and determine the structure of the reconstruction matrix. In this paper, it is using the number of the main frequency in the result of using signal fast Fourier transform to determine the rank of the Han

6、kel matrix, and through SNR(Signal to Noise Ratio) and MSE(Mean Square Error) to determine reconstruction matrix structure.Simulation and experiment validated this method. The results shows that the number of rank is double of the main frequency, and the best lines of reconstruction matrix is half o

7、f the length of the signal data. You can easy to choose the rank of the matrix and get a better noise elimination result.Keywords:Signal de-noising; Singular value decomposition; Fast Fourier transform; Signal to noise ratio; Mean square error目 录1 引言12 SVD分解消噪理论52.1 Hankel矩阵52.2 SVD分解的基本理论62.3 对测量信号

8、进行SVD分解63 有效秩阶次和重构矩阵结构的确定83.1 有效秩阶次的确定83.2 重构矩阵结构的确定144 消噪后的信号重构165 数值仿真175.1 Matlab仿真结果分析175.2 Matlab程序236 结束语26参考文献27基于奇异值分解的信号消噪技术1 引言 随着社会的发展，人类社会对石油的需求日益提高，海上采油区域不断扩大，有越来越多的海洋平台建造并投入使用，而这些海洋平台结构在复杂的服役环境中将受到设计载荷的作用以及各种突发性外在因素的影响而面临结构的损伤积累的问题，从而使结构的安全受到威胁。大型工程结构一旦出现事故，所带来的命和财产损失将是巨大的，对社会的影响更是深远和难

9、以估量的。特别对于海洋平台而言，其结构复杂，造价昂贵，一旦发生事故，不仅会对海洋环境造成很大的污染，还会带来不可估量的经济损失和人员伤亡，造成不好的社会政治影响。众所周知，海洋平台结构长期服役在恶劣的海洋环境中，并受到各种载荷的交互作用，如风载荷、海流、波浪载荷、冰载荷等，有时还要遭到地震、台风、海啸、船碰撞等意外打击，结构本身还要遭受环境腐蚀、海洋生物附着、海底冲刷等影响的作用。在这些恶劣的环境载荷长期作用下，再加上设计或使用的不当，结构容易产生各种形式的损伤，使结构的承载能力下降，严重的还会导致平台失效嘲。在国内外海洋开发工程中，曾发生过多起灾难性海洋平台事故，造成了巨大的人员伤亡、经济损

10、失以及不良的社会影响。 随着石油开采向海洋发展，海洋平台的数量成倍增加，合适的设计方法确保结构能够抵抗住不可预测的载荷造成的损伤，但是损伤在海洋平台结构的服役期间是不可避免的，确保人的生命安全和减少财产损失的唯一方法是诊断出结构的损伤，并能及时进行修复。由此可见，提高海洋平台结构及设备的可靠性，确保海洋作业安全的问题日益突出，新平台的质量评价、旧平台的残余寿命估计和在役平台的结构安全保证将成为日益突出的问题，海洋平台结构的健康监测与损伤诊断已成为刻不容缓的重要课题。 当前的结构损伤检测的方法很多，除了人工目测外，还有超声波、磁场法、放射法、热力场等局部检测方法。然而，较弱的视觉观测条件以及损伤

11、部位有可能被生长的海洋生物覆盖着，所以利用这类局部损伤检测技术方法对海洋平台结构进行损伤诊断是不可靠的。此外，这些技术要求结构的损伤区域是己知作为先决条件，要求配备特殊额外的测试设备和专业人员，因此，这些方法的检测成本较昂贵。与上述方法相比较，基于振动测试的结构健康监测技术是相对简单、成本较低的，被公认为是较有发展前景的全局性损伤诊断方法。这种方法的基本原理是：损伤将导致结构的系统刚度和阻尼矩阵发生改变，因而导致结构的动力特性参数(如结构的频响函数，模态参数等)的变化。换言之，结构动力特性参数能够作为结构损伤诊断的指标。这类方法最突出的优点是整个损伤诊断操作过程不会影响结构的正常工作。 在土木

12、工程领域，海洋平台结构、桥梁和大坝等工程结构被视为“系统”，而“识别”则意味着从振动测试数据中识别出结构的动力特性参数(模态参数)。对结构物而言，模态参数是结构的“指纹”，它是一系列独特的数据，能够反映结构本身的固有动力特性。每一个结构都有其固有的模态参数，并且如果结构动力特性发生变化了，那么结构的“指纹”也将发生相应的变化。因此，模态分析是结构动态设计以及设备故障诊断的重要方法。基于输出的模态参数识别方法利用的信息主要是系统的自由振动信号，要获得自由振动信号首先需获得结构的响应信号。由于环境激励的不充分和噪声等干扰因素的存在，响应信号中包含的有用信息十分微弱，尤其对海洋平台结构而言，其服役环

13、境恶劣而复杂，现场测试的信号中包含了较多的噪声成分，因此不能直接用于参数辨识，这是导致参数识别方法难以奏效的主要原因。因此，需要对信号进行消噪处理，从大量背景噪声中提取出可用于模态参数辨识的有用信号成分,剔除响应信号中的干扰因素，提取有用信息成为关键。此时，信号消噪技术研究变得尤为重要。 目前，有大量文献对信号降噪技术进行了研究，提出了多种降噪方法，如时域平均法1、小波降噪技术2、频域特征抽取技术3、自适应滤波技术4等。然而，各种方法在实际应用中都有各自的局限性，时域平均法需要有足够的数据量，并且在使用过程中必须有时标信息的支持；小波降噪和自适应滤波技术很大程度上依赖于滤波器性能；频域特征抽取

14、技术过度依赖于信号的幅值、频率、相位信息，计算起来很不方便，而且对多谱勒等变频信号而言，无法成功降噪。 近年也发展了一些其它方法，其中基于奇异值分解(SVD)的降噪技术因其计算方法简单易用引起了国内外相当一部分专家学者的重视。奇异值分解技术在声学、智能控制、电子学、信号处理等领域得到了广泛的应用。Fort等人5利用矩阵SVD方法和一个新的标准(dynamic mean evaluation，DME)确定模型阶次，再利用AR模型进行谱估计，将其应用到benchmark和Doppler信号分析中；Sanliturk等人6将Hankel矩阵和SVD算法结合，从复杂的噪声信号中获得较高精确度的频响函数

15、；Vrabie等人7在SVD算法的基础上引入独立分量分析的概念，可以在传感器相互干扰很大的情况下，在低通子空间很好的分离原波形，并应用于垂直地震剖面的分析；马寨璞8等人在卡尔曼滤波的基础上，提出了利用矩阵的奇异值将数据矩阵进行SVD分解的新的简化方法，得到能够描述原状态向量的新的较少维数向量的有效秩，并和HAMSON模式结合，利用渤海区域的SST进行实验得到验证。Eckart和Young9提出了截断奇异值分解技术，即只保留系统阶数范围内的几个最大奇异值,范围外的奇异值均设为零。该理论成为低秩逼近问题的核心,为后来的降噪技术发展提供了新的思路。近年来,低秩逼近问题己发展到结构低秩逼近,矩阵从一般

16、形式发展到特殊形式，如Hankel、Toeplitz矩阵等。已有很多学者提出结构完全最小二乘法。Aoki和Yue最早做了对结构完全最小二乘法问题的研究工作，在多年后的文献中才首次出现但这一问题的提法。Cadzowt,Bresler和Macovski提出了另一种求解方法，最终这些方法被证明只满足次优的L2优化准则,尽管如此,它们由于简单易用而被广泛采用。基于SVD分解的消噪技术,它是从矩阵的角度出发,将包含信号特征的矩阵分解到一系列奇异值和奇异值矢量对应的子空间中，近些年的研究对SVD算法进行了有益的应用和改进。该方法有两个关键: 1：如何确定分解后重构的有效秩阶次; 2：如何确定重构矩阵的行列

17、数。 针对有效秩阶次的选择， 常用的方法是试凑法和阈值法，均依赖于经验， 缺乏依据。目前,对此问题,已经出现了一些研究方法。如稳定图法，该方法通过在频谱图上标示出满足一定条件的稳定极点，并被认为是系统的真实极点。但这种方法不能完全排除噪声模态，特别是随着模型阶数的升高，一些拟合模型的拟合模态往往容易趋于稳定，用稳定图很难完全并正确确定模型阶次。有人提出了基于奇异值分解的模型定阶与降噪技术，该类方法一般将奇异值由大到小按降序排列，并将奇异值以最大值归一化，通过画出奇异值归一化曲线，在曲线上找到突降的位置，该处对应的奇异值个数即为模型阶次，也即信号中包含模态数目的两倍。对受噪声影响的数据，奇异值曲

18、线突降不明显，而是趋向于一条水平渐近线。一般认为奇异值曲线开始变为水平的点对应模型的阶次。此外还有朱启兵10提出了基于结构风险最小化原则的奇异值分解降噪方法，该方法依据统计学习理论， 把有效秩阶次的选择看作是一个学习过程，利用结构风险最小化原则来代替传统的经验风险最小化，从而自动得到奇异值分解降噪中矩阵的有效秩。王维11提出了基于非监督动态聚类算法来确定矩阵有效重构阶次的方法，该方法利用含噪声信号的奇异谱图中表征噪声的噪声平台平缓和集中的特性，通过向谱图纵轴投影，应用动态聚类合理确定噪声平台的边界，进而有效地确定奇异值分解降噪中矩阵的有效重构阶次。康春玉12采用主分量分析的方法，根据奇异值的大

19、小来确定有效秩的阶次。孙鑫晖13提出了通过奇异熵增量确定降噪阶次的方法，信号的奇异熵是信息熵的一种改进形式，能够反映信号包含信息量的多少。当信号受到宽频带噪声干扰，信号的奇异熵随着阶次升高一直增大。在较低阶次时，在信号与噪声的共同作用下，奇异熵增长速度较快。当达到一定阶次后，奇异熵增长速度放缓。这时信号的有效特征信息量已经趋于饱和，之后的增量是由于噪声所致。因此,选择奇异熵达到饱和阶次作为重构阶次，能够保留信号信息同时去除噪声。 针对重构矩阵的行列数的选择问题，普遍采取的方法是根据具体信号选择不同行列数进行试凑，这种试凑法需要大量的计算，并且严重依赖使用者的信号分析经验。Kanjilal14提

20、出了通过奇异值比谱来确定行列数的方法，但是该方法在信号由多个周期分量组成或噪声较严重的情况下效果并不明显。赵学智15提出通过分析所有分解分量信号所含信息量的变化趋势来确定合理的矩阵结构。上述方法在实际应用中取得了较好的效果，但也存在着一定的局限性。本文将Hankel矩阵与SVD分解相结合，首先对测量信号构造的Hankel矩阵进行奇异值分解，得到测量信号的奇异值（包括有用信号和噪声信号的奇异值），再利用测量信号快速傅立叶变换结果中主频率的个数来确定有效秩阶次，剔除噪声信号的奇异值，接着通过消噪信号的信噪比和均方差大小确定重构矩阵结构，得到有用信号的矩阵，最后通过反对角线平均法得到消噪后的信号数据

21、。2 SVD分解消噪理论2.1 Hankel矩阵 在信息与信号处理领域，要想将有用信号不失真地变换和处理几乎是不可能的。误差的来源有：信息传输处理时，信道或设备不理想造成的；传输处理过程中串入的一些其它信号也能引起误差。噪声是使信号产生失真的误差源。来自外部的噪声也称为干扰。从功率谱的角度看，如果一个随机过程的功率谱密度是常数，无论是什么分布，都称为白噪声。白噪声的频率成分非常丰富。假设有一个载有信息的、以时间t为变量的信号，由于在传输过程中受到了某种加性噪声或干扰的污染，致使观察到的信号发生变化，其可以表示为： 对于一个测得的信号，其中为有用信号，为噪声信号，基于相空间重构理论16，可以由其

22、构造Hankel矩阵 消噪的目的就是从A矩阵中求出的最佳逼近矩阵，再求出。 2.2 SVD分解的基本理论 设是一个秩为r的维矩阵，则的奇异值分解是指，存在矩阵和以及，使得： 其中，U，V分别为，维正交矩阵，0为零元素矩阵。 为维对角阵，其对角线元素为矩阵H的非零奇异值，且以非增顺序排列，即 ，且有。 矩阵的秩为r，从式（2）中除去H的零奇异值，得到奇异值分解的精简形式 其中：，分别为U，V的第个行向量。2.3 对测量信号进行SVD分解 根据式（2）对构造的Hankel矩阵（即式（1）进行SVD分解： 都为维矩阵，为维正交矩阵，为维正交矩阵，为维对角矩阵，其对角元素分别为的奇异值，同时还可看出。

23、上述奇异值分解式根据式（4）还可以进一步写成 由上式可以看出，如果信号中没有噪声或信噪比特别高，则矩阵是奇异的， 即的数目，的大小和系统有关。如果有噪声或信噪比不高，则矩阵是非奇异的，即的数目。因源信号是由有用信号和噪声信号共同组成，则矩阵也是由有用信号和噪声信号共同组成的矩阵，那么矩阵的奇异值可以反映信号和噪声信号能量集中的情况。前r个较大的奇异值将主要反映有用信号，较小的奇异值则主要反映噪声信号，把这部分反映噪声的奇异值置零就可以去除信号中的噪声。再利用奇异值分解的逆过程得到矩阵，即矩阵的秩为r的最佳逼近矩阵。相对于其噪声已被大大压缩。将中对应的元素相加平均，就可以得到降噪后的信号 。3

24、有效秩阶次和重构矩阵结构的确定3.1 有效秩阶次的确定 由实测经验可知，当实际测量环境很好时，测量得到的信号一般为光滑曲线；而当受到外界的随机噪声干扰时，测量得到的信号中就会含有大量的“毛刺”。利用噪声污染信号构造Hankel 矩阵进行奇异值分解降噪，就是对含噪信号进行逼近、剔除毛刺的过程。经过对大量仿真结果的研究，发现去噪结果中剩余噪声对信号的干扰影响能够通过信号中“毛刺”的数量及大小来判断，这与实测经验相吻合。当选取的奇异值数目较小时，大部分噪声都被剔除掉了，但也损失了大量有用信号。随着奇异值数目的逐渐增大，有用信号信息趋于完整，但噪声的干扰也会逐渐增加。当奇异值数目达到某一值时，降噪后的

25、信号既保留了大部分有用信息，也剔除掉了大部分噪声，这就是要选取的最佳奇异值数目。由第个非零奇异值重构得到重构信号分量， 分别对这些分量信号进行傅里叶变换后发现，分量的频率成分均为源信号的频率成分组成，而由较大的奇异值重构得到的分量信号其频率成分与源信号中主频率相对应。显然， 当有用信号未被噪声完全淹没时，源信号中主频率是有用信号的频率。因此可以通过奇异值与有用信号频率之间的某种对应关系来确定有效秩的阶次，从而达到降噪的效果。 给定一个源信号，其中：为有用信号，为强度为1的高斯白噪声。将和快速傅立叶变换FFT的结果如图1和图2所示，比较两张图，发现它们都包含两个主频率成分，图2不含其他频率，由此

26、可知这两个主频率为有用信号的频率成分，其他均为噪声的频率成分。图1 源信号的时域图和功率谱密度图图2 有用信号的时域图和功率谱密度图分别对源信号和有用信号利用式(1)构造相应的Hankel矩阵，并对其进行奇异值分解。随着矩阵的行数变化，其奇异值的个数分布如图3和图4所示：图3 源信号奇异值个数随矩阵行数变化情况图4 有用信号奇异值个数随矩阵行数变化情况 可以看到：经过对无噪声信号和受噪声污染信号的研究，发现一般由无噪声理想信号构造的Hankel矩阵的大部分奇异值为零。根据非零奇异值的数目，可以很容易地判断矩阵的有效阶次。因噪声具有随机和不相关的特点，因此，由受随机噪声污染信号构造的Hankel

27、 矩阵呈列满秩或行满秩状态（取决于行和列哪个维数更小）。具体分析如下：1：当行数或列数小于主频个数的两倍时，无论是有用信号还是源信号，其非零奇异值的个数等于矩阵行数和列数之中的最小值。2：当行数大于主频个数的两倍时，有用信号的非零奇异值一直是4个，且不随矩阵行数的变化而变化；当行数小于列数时，源信号的非零奇异值等于行数的值，当行数大于列数时，源信号的非零奇异值等于列数的值。由此可知，源信号经过奇异值分解等到的对角矩阵的对角元素都不为零。 再对源信号和有用信号的矩阵经过奇异值分解之后得到的对角矩阵的对角元素大小进行分析，如图5和图6所示：图5 源信号奇异值大小随矩阵行数变化情况因上图点数较多，不

28、能确定大奇异值个数，随选择较短行数以便观察。图51 矩阵行数在1-20之间时，奇异值大小情况图52 矩阵行数在70-120之间时，奇异值大小情况图6 有用信号奇异值大小随矩阵行数变化情况对比图5和图6中的奇异值曲线可知，每个奇异值受到噪声干扰的程度不同：较大的奇异值受噪声干扰的影响较小，而较小的奇异值受噪声干扰的影响则较大。由于由受噪声污染信号构造的Hankel 矩阵的每个奇异值都是由“有用信号”和“噪声信号”两部分组成，实际选取奇异值数目就是一个对有用信号和噪声信号进行取舍的过程。从图6中可以看出，随着行数的增加，奇异值的个数恒为4，是源信号快速傅立叶变换中主频个数的2倍；而图5中较大奇异值

29、的个数也一直为4，说明其主要反映有用信号的信息，将这些奇异值称为大奇异值，而图5中其他奇异值相对较小且分布比较集中，说明其反映出了噪声的特点。所以，在消噪过程中，选取的奇异值个数越接近有用信号的奇异值个数，就越是接近不受噪声污染的原始信号，当奇异值数目过小时，就会漏掉有用信号，奇异值过大时，降噪后的信号中就会含有大量的噪声。由此可知：有效秩的阶次与源信号的主频个数存在一个确定的关系， 即有效秩的阶次为源信号快速傅里叶变换后主频的个数的两倍。利用这种倍数关系可以确定有效秩的阶次，取得最佳的奇异值数目，不仅可以过滤掉大量的噪声，而且还能最大程度地保留原始信号的信息。3.2 重构矩阵结构的确定确定有

30、效秩的阶次之后，就可以提取出信号对角矩阵中的有用奇异值，即剔除了测试信号中的噪声奇异值，接着要重构有用信号的矩阵，这需要确定重构矩阵的结构，即矩阵的行列数，因为不同行列数会导致不同的消噪效果。本文根据降噪效果的好坏来确定最佳的矩阵结构。 降噪效果一般用信号的均方误差( mean square error，简称MSE)和信噪比(signal to noise ratio，简称SNR) 来衡量: MSE 越小，SNR 越大，降噪效果越好。MSE 与SNR 的定义形式如下 其中：为含噪声信号的第k个数据点；为无噪声信号的第k个数据点；N为信号长度。 选取不同行数L(L>5)对信号(N = 20

31、0)，噪声强度为1dB，进行奇异值分解与重构，降噪信号的SNR与MSE随L变化如图7所示。可以看出，当L增加到一定程度时，降噪信号的SNR基本稳定在22dB左右，仅有小的波动，此时MSE的变化也呈现出这种规律。由结果可知，当L=102时，降噪信号的MSE取得最小值，SNR取得最大值，此时降噪效果最好。图7 均方误差和信噪比随行数变化 进行分析后发现，SNR与MSE分别取最大和最小值时的L值不一定相等， 但最佳L值基本出现在处的一个领域内，并且L在该邻域内取值时，降噪效果较好且差异较小，均能满足要求。因此，重构矩阵的结构可以根据N来确定，实验应用中不妨取( 当N不是偶数时，舍弃最后一个数据点，不

32、影响最终结果)。4 消噪后的信号重构 基于SVD分解的消噪技术,它是从矩阵的角度出发,将包含信号特征的矩阵分解到一系列奇异值和奇异值矢量对应的子空间中，通快速过傅立叶变换结果中主频率个数确定有用信号的奇异值个数以保留矩阵的有用信号奇异值，并将其余奇异值置零以剔除信号噪声的奇异值。再根据有用信号奇异值重构有用信号的矩阵，因为不同行列数会导致不同的消噪效果，根据均方误差和信噪比可以知道当行数为信号数据长度的一半时（即矩阵为方阵），降噪效果最好。奇异值个数和重构矩阵行数确定后，就可以依据此得到消噪后的信号矩阵，但矩阵并不等于由真实信号构成的Hankel矩阵A，不是严格的Hankel矩阵。通过观察式（

33、1），可以看出真实信号与矩阵A的各元素之间存在如下关系： 其中，。也就是对矩阵A的反对角线求平均值即可得到真实信号在每一时刻的值，。根据这一思路，同样对矩阵的反对角线求平均值，从而得到经过消噪后，信号在每一时刻的估计值，即： 其中，。5 数值仿真在信息与信号处理领域，将有用信号不失真地变换和处理是不可能的，因为在信息传输处理时，信道或设备的不理想会造成误差，或者在传输处理过程中会串入一些其它信号即噪声。本文将Hankel矩阵与SVD分解相结合，首先对测量信号构造的Hankel矩阵进行奇异值分解，因源信号是由有用信号和噪声信号共同组成，则矩阵也是由有用信号和噪声信号共同组成的矩阵，那么矩阵的奇异

34、值可以反映信号和噪声信号能量集中的情况。前r个较大的奇异值将主要反映有用信号，较小的奇异值则主要反映噪声信号，把这部分反映噪声的奇异值置零就可以去除信号中的噪声。本文利用测量信号快速傅立叶变换结果中主频率的个数来确定有效秩阶次，剔除噪声信号的奇异值。因为重构矩阵行数的不同会影响降噪效果，通过消噪信号的信噪比和均方差大小比较可发现，SNR与MSE分别取最大和最小值时的矩阵行数L值不一定相等，但最佳L值基本出现在处的一个领域内，N为信号数据长度，并且L在该邻域内取值时，降噪效果较好且差异较小。最后根据重构矩阵和信号方程之间的关系，通过反对角线平均法得到消噪后的信号数据。5.1 Matlab仿真结果

35、分析 基于上述分析，对于一个含噪声的测试信号，其降噪的基本步骤如下:(1) 取信号数据长度的一半作为重构矩阵的行数，根据式(1)构造Hankel矩阵并进行奇异值分解;(2) 对信号进行快速傅里叶变换，确定主频个数n，以2n作为有效秩的阶次;(3) 用前2n个奇异值根据式(4)进行重构，得到重构矩阵，将中对应的元素根据式（10）相加后平均就可得到降噪后的信号。 分别用不同频率成分的信号对该方法进行验证：信号1：，高斯白噪声强度为1dB，10dB，15dB；信号2：,高斯白噪声强度为10dB；信号3：，高斯白噪声强度为10dB；信号4：，高斯白噪声强度为10dB。信号1，2，3取数据长度为200，

36、可将其构成矩阵的行数设为100；信号4取数据长度为400，可将其构成矩阵的行数设为200。表1 信号1消噪前后SNR和MSE的变化 为了考察这种降噪方法在不同噪声水平下的表现,分别进行了对含1dB, 10dB, 15dB高斯白噪声信号1的消噪处理, 发现在不同信噪比下, 通过这种方法得到的消噪波形都能较好地保留目标信号的波形特征,分别对各个信号进行噪声消除，得到结果如下：信号1噪声强度/dB源信号信噪比SNR消噪后信号信噪比SNR源信号均方误差MSE消噪后信号均方误差MSE136.63237.4870.0558960.0490791011.43114.4680.866270.65448157.

37、88239.69091.8171.6677图8 信号1-噪声1dB图9 信号1-噪声10dB图10 信号1-15dB 从表1的可以看出，降噪后信号的信噪比SNR都提高了，均方误差MSE都明显降低了，这说明本文给出的降噪方法是有效的。再看图8到图10中的波形，可以看出，随着噪声强度的增加，源信号波形失真越严重，功率谱密度图中噪声频率的个数也逐渐增多，并且逐渐淹没主频率。当噪声强度不至于淹没有用信号时，降噪后的波形与原波形吻合较好，证明了这种方法的降噪效果。图11 信号2-10dB图12 信号3-10dB图13 信号4-10dB 从傅立叶变换结果可以看出信号2，3，4的主频个数分别为2，3，4，那

38、么可以确定重构矩阵的有效秩阶次分别为4，6，8，由此可以剔除信号中的噪声信号。因为信号2和3的数据长度为200，信号4的数据长度为400，根据上文分析得出的当重构矩阵的行数为数据长度的一半时信噪比最大的理论，信号2和3重构矩阵行数为100，信号4的重构矩阵的行数为200。因为当信号受到噪声干扰时，测量得到的信号中含有大量的“毛刺”，利用噪声污染信号构造Hankel 矩阵进行奇异值分解降噪，就是对含噪信号进行逼近、剔除毛刺的过程，所以可以根据去噪结果中“毛刺”的数量及大小来判断消噪效果。观察各图波形，发现消噪后信号波形基本都是光滑曲线，几乎没有“毛刺”，这表明大部分噪声都被剔除掉了。并且发现消噪

39、后信号波形与有用信号波形吻合较好，证明了这种消噪方法的实用性。实验结果表明, 利用源信号主频个数来确定有效秩的阶次以及取信号数据长度的一半确定重构矩阵的行数的方法可以得到较好的降噪效果，该方法可以获得较高的信噪比,同时也较好的保留了原信号的特征波形，证明该方法是十分有效的。5.2 Matlab程序clear;Fs=1;%采样频率n=400;%采样个数 L=n/2;%矩阵最佳行数，数据长度一半 A=1:L;%第一列 B=L:n;%最后一行 c=hankel(A,B);%构造hankel矩阵 h=sin(0.01*pi*c)+2*cos(0.03*pi*c)+cos(0.05*pi*c)-2*si

40、n(0.07*pi*c);%有用信号 y=wgn(size(c,1),size(c,2),10);%加噪声 h=y+h;%源信号 for i=1:n p=max(1,i-L+1); q=min(n-L+1,i); he=0; for j=p:q he=he+(h(i-j+1,j); end ca=q-p+1; x1(i)=he/ca; end figure(2); subplot(311); i=1:n; plot(i,x1,'r');%输出源信号波形 grid; title('消噪前'); axis(0,n,-6,6);window=boxcar(length

41、(x1);%矩形窗nfft=n;%采样点数Pxx,f=periodogram(x1,window,nfft,Fs);%直接法求功率谱密度subplot(312);plot(f,10*log10(Pxx);title('功率谱密度图');xlabel('频率Hz');ylabel('功率谱密度');axis(0,0.5,0,30);%画出源信号功率谱密度图grid; U,S,V=svds(h,8);%奇异值分解，确定有效秩阶次 h2=U*S*V'%重构有用信号的矩阵 for i=1:n p=max(1,i-L+1); q=min(n-L+1

42、,i); he=0; for j=p:q he=he+(h2(i-j+1,j); end ca=q-p+1; x2(i)=he/ca; end subplot(313); i=1:n; plot(i,x2);%输出消噪后信号波形 hold on; grid; x0=sin(0.01*pi*i)+2*cos(0.03*pi*i)+cos(0.05*pi*i)-2*sin(0.07*pi*i); plot(i,x0,'-r');%输出有用信号波形 title('消噪后'); axis(0,n,-6,6); 6 结束语 本文将Hankel矩阵与奇异值分解法相结合，根据

43、测量信号中有用奇异值是不变的，将测量信号构造的Hankel矩阵进行奇异值分解，得到一个对角矩阵，对角矩阵中包含的元素即测量信号中的奇异值（包括有用信号和噪声信号的奇异值）。因为测量信号快速傅立叶变换结果中主频率的个数的2倍即有效秩阶次，可以据此剔除噪声信号的奇异值。接着通过消噪信号的信噪比和均方差大小确定重构矩阵结构，利用奇异值分解的逆过程即可得到有用信号的矩阵，最后通过反对角线平均法得到消噪后的信号数据。仿真结果表明，对于不同频率，不同噪声强度的信号，该方法可以剔除大部分噪声，获得较高的信噪比,同时也能较好的保留原信号的特征波形，证明该方法是十分有效的。参考文献1G.J.JanacekPractical time seriesLondon：Arnold，20012YY Kim，JC Hong，NYLeeFrequency response function estimation via a robustwavelet de-noising methodJ Sound Vib2001，244(4)：635-6493LLPresti，GOlm