CN110826017A - 一种基于参数优化Hankel矩阵和奇异值分解的信号去噪方法 - Google Patents

Abstract

一种基于参数优化Hankel矩阵和奇异值分解的信号去噪方法，包括：选择矩阵的列数，对于不同的列数，利用含噪信号分别构造Hankel矩阵；对Hankel矩阵分别进行奇异值分解，计算奇异值的四阶中心矩均方根，取最大值对应的列数为最优列数；对最优列数的Hankel矩阵，计算其奇异值贡献率的差分谱，利用后向门限比较法选择有效信号奇异值；利用有效信号奇异值进行重构，得到去噪后的Hankel矩阵，并利用相空间法对去噪后的Hankel矩阵进行重构，得到去噪后的信号。本发明可较好衡量有效信号奇异值和噪声奇异值之间的分布距离，有利于选择最优列数；可较好衡量奇异值的波动，有利于准确选择有效信号奇异值，计算简单，可有效提高信号去噪效果，适于雷达和通信信号的分析与处理。

Description

一种基于参数优化Hankel矩阵和奇异值分解的信号去噪方法

技术领域

本发明涉及属于数字信号处理领域，更具体地说，涉及一种基于参数优化Hankel矩阵和奇异值分解的信号去噪方法。

背景技术

含噪信号的噪声抑制是信号处理中的一个重要步骤，对于后续的信号检测与分析具有重要意义，是信号处理的研究热点之一。

常用的信号去噪方法有傅立叶变换法、小波阈值法、经验模态分解法和子空间投影法等。傅立叶变换法主要通过设计频域滤波器实现噪声抑制，它只能反映信号的整体特征, 不适用于频率随时间变化的非平稳信号。小波阈值法中小波基函数、分解层数以及阈值的选择都依赖于主观经验，缺乏自适应性。经验模模态分解(Empirical ModeDecomposition， EMD)法具有较好自适应性，但是存在模态混叠的问题。集合经验模态分解(Ensemble Empirical Mode Decomposition,EEMD)是一种改进的EMD方法，它可在一定程度上抑制模态混叠问题，但不能完全消除模态混叠的影响。子空间投影法是一种非参数化分析方法，它将含噪信号分解成合适的子空间分量，选择含有效信号的子空间分量进行重构以抑制噪声。子空间投影法包括主分量分析(Princinal Componet Analysis，PCA)和独立分量分析 (Independent Component Analysis,ICA)两种。ICA方法要求各个源信号之间瞬时统计独立，且至多只能有一个源信号为高斯分布，主要应用于盲源分离的场合。PCA方法的核心是利用特定的正交矩阵对信号矩阵进行正交变换，得到相互正交的对角主成分矩阵，具有较好的数值稳健性和自适应性，应用范围较广。奇异值分解(Singular ValueDecomposition,SVD) 是PCA的一种主要实现方法，奇异值对应的信号分量称为主分量，SVD去噪的基本原理就是通过选择有效信号的奇异值重构来获得去噪后的信号。

对于一维信号，利用奇异值分解法去噪的首要问题是如何构造合适的矩阵。目前常见的矩阵包括Toeplitz矩阵、周期矩阵和Hankel矩阵。其中，基于Hankel矩阵的奇异值分解和小波分析有类似的效果，在信号去噪中得到了广泛应用。

对于Hankel矩阵，列数是其唯一的参数，因此矩阵构造的核心是如何选择最优的列数。目前列数的选取主要有两种方法。一种方法是选取矩阵列数为信号长度的一半，该方法只适合于连续信号，对于间歇性信号效果较差。另一种方法是根据奇异值比谱来确定，但是该方法只利用前两个奇异值信息，仅适用于单周期分量的信号，不适用于多个周期分量的信号。对于奇异值分解，其核心是如何选择有效信号奇异值。目前常用的方法有差分谱法、奇异值曲线法和奇异值曲谱率法。差分谱法通过判断突变点来选择有效信号奇异值，主要依赖经验。奇异值曲线法需要噪声的先验知识，适应性较差。奇异值曲谱率法计算复杂，容易受噪声强度的影响。因此，如何以较小的计算复杂度选择最优矩阵列数和有效信号奇异值，对于提高基于Hankel矩阵和奇异值分解法的去噪性能具有重要意义。

发明内容

为了解决基于Hankel矩阵的奇异值分解方法在去噪时存在的问题，本发明提供了一种基于参数优化的Hankel矩阵和奇异值分解的去噪方法。该方法利用奇异值的四阶中心矩均方根对Hankel矩阵的列数进行优化，可使不同信号奇异值分布之间距离最大，有利于后续的判断；利用奇异值贡献率的差分谱选择有效信号奇异值，稳定性好，其计算简单，可有效提高信号的去噪性能。

本发明解决其技术问题，所采用的基于参数优化Hankel矩阵和奇异值分解的信号去噪方法，包括以下步骤：

(1)输入含噪信号；

(2)选择矩阵列数，利用含噪声信号构造Hankel矩阵；

(3)对利用含噪声信号构造的Hankel矩阵进行奇异值分解；

(4)计算矩阵奇异值的四阶中心矩均方根，

(5)改变矩阵列数，重复步骤(2)－(4)，取四阶中心矩均方根最大值对应的列数为最优列数；

(6)对于最优列数所对应的利用含噪声信号构造的Hankel矩阵，计算其奇异值贡献率的差分谱，利用后向门限比较法选择有效信号奇异值；

(7)对有效信号奇异值采用奇异值分解逆运算进行重构，得到去噪后的Hankel矩阵；

(8)利用相空间法对去噪后的Hankel矩阵进行重构，得到去噪后的信号。

进一步地，在本发明的基于参数优化Hankel矩阵和奇异值分解的信号去噪方法中，所述步骤(2)中利用含噪声信号构造Hankel矩阵，具体如下：

对于长度为N的含噪信号X＝[x(1),x(2),…,x(N)]，选择矩阵列数n，构造Hankel矩阵如公式(1)所示：

其中，列数取值范围为

floor表示向下取整，m＝N-n+1且m≥n，则A⁽ⁿ⁾∈R^m×n。

进一步地，在本发明的基于参数优化Hankel矩阵和奇异值分解的信号去噪方法中，所述步骤(3)中对利用含噪声信号构造的Hankel矩阵进行奇异值分解，如下公式(2)所示：

A⁽ⁿ⁾＝UΣV^T (2)

其中，A⁽ⁿ⁾为利用含噪声信号构造的Hankel矩阵，U∈R^m×m和V∈R^n×n为正交矩阵，T表示转置；Σ∈R^m×n为非负对角阵，

是矩阵A⁽ⁿ⁾的奇异值，且r＝min(m,n)＝n。

进一步地，在本发明的基于参数优化Hankel矩阵和奇异值分解的信号去噪方法中，所述步骤(4)中，奇异值四阶中心矩均方根的计算公式如下公式(3)：

其中，为利用含噪声信号构造的Hankel矩阵A⁽ⁿ⁾的第i个奇异值，

为奇异值的均值。

进一步地，在本发明的基于参数优化Hankel矩阵和奇异值分解的信号去噪方法中，所述步骤(5)中，最优列数的计算如下公式(4)：

进一步地，在本发明的基于参数优化Hankel矩阵和奇异值分解的信号去噪方法中，步骤(6)中，选择有效信号奇异值的具体步骤如下：

(1)对于最优列数n_opt，计算奇异值贡献率，如下公式(5)所示：

其中，

为矩阵

的第k个奇异值；

(2)计算η_k的差分谱，如下公式(6)所示：

dη_k＝η_k-η_k+1,k＝1,2,…,n_opt-1 (6)

(3)选取门限T,采用后向比较法确定有效信号奇异值的阶次，如下公式(7)所示：

km＝max[k|dη_k＞T],k＝n_opt-1,n_opt-2,…,1 (7)

其中，max[]表示取最大值，此时选择前km个奇异值为有效信号奇异值。

进一步地，在本发明的基于参数优化Hankel矩阵和奇异值分解的信号去噪方法中，步骤(7)中，采用奇异值分解逆运算进行矩阵重构，重构公式如下公式(8)所示：

其中，

为去噪后的Hankel矩阵，u_i、v_i的定义如下：

其中

m_opt＝N-n_opt+1， n_opt表示最优列数，且满足下述公式(9)：

进一步地，在本发明的基于参数优化Hankel矩阵和奇异值分解的信号去噪方法中，步骤(8)中，对去噪后的Hankel矩阵进行重构进行相空间重构，重构公式如下所示：

其中，为去噪后的Hankel矩阵，X_s为进行相空间重构后信号，即去噪后的信号， m_opt＝N-n_opt+1，n_opt表示最优列数。

实施本发明的基于参数优化Hankel矩阵和奇异值分解的信号去噪方法，具有以下有益效果：

1.本发明利用奇异值的四阶中心矩均方根来选择Hankel矩阵的最优列数，可使得最优列数下有效信号奇异值和噪声奇异值之间的分布距离最大，有利于更好的分离有效信号和噪声；

2.本发明利用奇异值贡献率可衡量不同奇异值所包含的信息量，其差分谱可反映该信息量的突变情况，有利于后续的判断；

3.本发明利用后向门限比较法选择来有效信号奇异值，克服了传统差分谱只取最大值的局限，可准确选择有效信号奇异值，提高去噪效果。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1为本发明一实施例的实现流程图；

图2(a)为原始理想信号1，图2(b)为含噪信号1；

图3为不同列数下含噪信号1构造的Hankel矩阵的奇异值四阶中心矩均方根；

图4(a)－图4(d)为列数n分别取100，200，350和500时，含噪信号1的去噪结果；

图5为不同列数下含噪信号1去噪后的信噪比；

图6(a)为原始理想信号2，图6(b)为含噪信号2；

图7为不同列数下含噪信号2构造的Hankel矩阵的奇异值四阶中心矩均方根；

图8(a)－图8(d)为列数n分别取50，150，300和500时，含噪信号2的去噪结果；

图9为不同列数下含噪信号2去噪后的信噪比。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

如图1所示，本实施例所述的基于参数优化Hankel矩阵和奇异值分解的信号去噪方法，包括以下步骤：

1.输入含噪信号；

2.选择矩阵列数，利用含噪声信号构造Hankel矩阵；

3.对利用含噪声信号构造的Hankel矩阵进行奇异值分解；

4.计算矩阵奇异值的四阶中心矩均方根，

5.改变矩阵列数，重复步骤2－4，取四阶中心矩均方根最大值对应的列数为最优列数；

6.对于最优列数所对应的利用含噪声信号构造的Hankel矩阵，计算其奇异值贡献率的差分谱，利用后向门限比较法选择有效信号奇异值；

7.对有效信号奇异值采用奇异值分解逆运算进行重构，得到去噪后的Hankel矩阵；

8.利用相空间法对去噪后的Hankel矩阵进行重构，得到去噪后的信号。

以图2(a)与图2(b)、图6(a)与图6(b)所示的2组原始理想信号和含噪信号为例，来说明本发明实施例的信号去噪方法，包括以下步骤：

1.输入含噪信号

(1)含噪信号1

原始理想信号1为雷达和通信中常用的多分量连续波信号。该信号包含两个频率和幅度不同的正弦波，如图2(a)所示，表达式如下公式(11)：

s₁(t)＝2sin(2πf₁t)+sin(2πf₂t) (11)

其中，f₁为10KHz，f₂为25KHz，信号采样率为2000KHz，采样点数为1000，对应的时间长度为0.5ms。

含噪信号1图2(b)所示，表达式如下公式(12)：

x₁(t)＝s₁(t)+n₁(t) (12)

其中，n₁(t)为高斯白噪声，信噪比为0dB。

(2)含噪信号2

原始理想信号2为雷达中常用的脉冲波信号，脉冲波形为Ricker子波。该信号包含两个幅度不同的脉冲波，如图6(a)所示，表达式如下公式(13)：

其中，信号采样率为1GHz，采样点数为1000，对应的时间长度为1us，两个脉冲的中心频率f₀均为20MHz，分别位于0.15us和0.75us处。

含噪信号2如图6(b)所示，表达式如下公式(14)：

x₂(t)＝s₂(t)+n₂(t) (14)

其中，n₂(t)为高斯白噪声，信噪比为0dB。

2.选择矩阵列数n，利用式(1)对两组含噪声信号x₁(t)和x₂(t)构造Hankel矩阵，记为

和

3.利用式(2)对矩阵和进行奇异值分解，得到两组奇异值，记为

和

4.利用式(3)计算两组奇异值的四阶中心矩均方根，得到

和

5.重复步骤2－4，其中列数取值范围为50≤n≤500，步进为50，这样得到一系列的

和

分别如图3和图7所示，取其中最大值对应的列数为最优列数，得到n_1opt＝350，n_2opt＝150；

6.对于最优列数n_1opt和n_2opt，利用式(5)和式(6)分别计算奇异值贡献率的差分谱，选取门限T＝0.01，利用式(7)确定有效信号奇异值的阶次km₁＝4，km₂＝10；

7.将km₁和km₂代入式(8)，得到去噪后的Hankel矩阵

和

8.利用式(10)对

和

进行相空间重构，得到去噪后的信号xd₁(t)和xd₂(t)，分别如图4(c)和图8(b)所示。

为了对比去噪效果，图4(a)－图4(d)和图8(a)－图8(d)中分别给出了x₁(t)和x₂(t)在不同列数下的去噪结果。

此外，为了量化评估去噪效果，图5和图9分别给出了x₁(t)和x₂(t)在不同列数下去噪后的信噪比(Signal to Noise Ratio,SNR)，信噪比的表达式如下：

其中，N为信号长度，s_i为原始理想信号，xd_i为去噪后的信号。

从图4(a)－图4(d)和图5中可看出，对于含噪信号1，当n＝100时，重构信号含有较多的噪声，此时SNR为11.53dB；当n＝200时，重构信号在一些点存在少量噪声和波形畸变，此时SNR为17.57dB；当n＝350时，重构信号在比较光滑，波形细节保留较好，此时SNR 为20.95dB；当n＝500时，重构信号仍然存在少量噪声和波形畸变，此时SNR为19.52dB。

从图8(a)－图8(d)和图9中可看出，对于含噪信号2，当n＝50时，重构信号含有较多的噪声，此时SNR为5.15dB；当n＝150时，重构信号中噪声得到较好抑制，此时SNR为6.70dB；当n＝300时，重构信号较好的抑制了噪声，但同时抑制了第2个Ricker子波有效信号，此时SNR为5.76dB；当n＝500时，重构信号在抑制噪声的同时也去除了第2个Ricker 子波有效信号，此时SNR为3.68dB。

上述2组信号的实验结果表明，本发明可准确选择最优矩阵列数和有效信号奇异值，具有较好的自适应性和稳定性，可改善基于Hankel矩阵和奇异值分解方法的去噪效果。

上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

Claims (8)

1.一种基于参数优化Hankel矩阵和奇异值分解的信号去噪方法，其特征在于，包括以下步骤：

(1)输入含噪信号；

(2)选择矩阵列数，利用含噪声信号构造Hankel矩阵；

(3)对利用含噪声信号构造的Hankel矩阵进行奇异值分解；

(4)计算矩阵奇异值的四阶中心矩均方根，

(5)改变矩阵列数，重复步骤(2)－(4)，取四阶中心矩均方根最大值对应的列数为最优列数；

(6)对于最优列数所对应的利用含噪声信号构造的Hankel矩阵，计算其奇异值贡献率的差分谱，利用后向门限比较法选择有效信号奇异值；

(7)对有效信号奇异值采用奇异值分解逆运算进行重构，得到去噪后的Hankel矩阵；

(8)利用相空间法对去噪后的Hankel矩阵进行重构，得到去噪后的信号。

2.根据权利要求1所述的一种基于参数优化Hankel矩阵和奇异值分解的信号去噪方法，其特征在于，所述步骤(2)中利用含噪声信号构造Hankel矩阵，具体如下：

对于长度为N的含噪信号X＝[x(1),x(2),…,x(N)]，选择矩阵列数n，构造Hankel矩阵如下：

其中，列数取值范围为floor表示向下取整，m＝N-n+1且m≥n，则A⁽ⁿ⁾∈R^m×n。

3.根据权利要求2所述的一种基于参数优化Hankel矩阵和奇异值分解的信号去噪方法，其特征在于，所述步骤(3)中对利用含噪声信号构造的Hankel矩阵进行奇异值分解，如下所示：

A⁽ⁿ⁾＝UΣV^T

其中，A⁽ⁿ⁾为利用含噪声信号构造的Hankel矩阵，U∈R^m×m和V∈R^n×n为正交矩阵，T表示转置；Σ∈R^m×n为非负对角阵，

是矩阵A⁽ⁿ⁾的奇异值，且r＝min(m,n)＝n。

4.根据权利要求2所述的一种基于参数优化Hankel矩阵和奇异值分解的信号去噪方法，其特征在于，所述步骤(4)中，奇异值四阶中心矩均方根的计算公式如下：

其中，

为利用含噪声信号构造的Hankel矩阵A⁽ⁿ⁾的第i个奇异值，

为奇异值的均值。

5.根据权利要求4所述的一种基于参数优化Hankel矩阵和奇异值分解的信号去噪方法，其特征在于，所述步骤(5)中，最优列数的计算如下：

6.根据权利要求3所述的一种基于参数优化Hankel矩阵和奇异值分解的信号去噪方法，其特征在于，步骤(6)中，选择有效信号奇异值的具体步骤如下：

(1)对于最优列数n_opt，计算奇异值贡献率，如下所示：

其中，

为矩阵

的第k个奇异值；

(2)计算η_k的差分谱，如下所示：

dη_k＝η_k-η_k+1,k＝1,2,…,n_opt-1；

(3)选取门限T,采用后向比较法确定有效信号奇异值的阶次，如下所示：

km＝max[k|dη_k＞T],k＝n_opt-1,n_opt-2,…,1；

其中，max[]表示取最大值，此时选择前km个奇异值为有效信号奇异值。

7.根据权利要求6所述的一种基于参数优化Hankel矩阵和奇异值分解的信号去噪方法，其特征在于，步骤(7)中，采用奇异值分解逆运算进行矩阵重构，重构公式如下所示：

其中，

为去噪后的Hankel矩阵，u_i、v_i的定义如下：

其中

m_opt＝N-n_opt+1，n_opt表示最优列数，且满足：

8.根据权利要求1所述的一种基于参数优化Hankel矩阵和奇异值分解的信号去噪方法，其特征在于，步骤(8)中，对去噪后的Hankel矩阵进行重构进行相空间重构，重构公式如下所示：

其中，

为去噪后的Hankel矩阵，X_s为进行相空间重构后信号，即去噪后的信号，m_opt＝N-n_opt+1，n_opt表示最优列数。

Images