CN104715146A - 一种水声信号的粒子滤波降噪方法 - Google Patents

Abstract

本发明提供一种水声信号的粒子滤波降噪方法，包括：步骤(a)、利用遗传算法对实测水声信号建立数学模型；步骤(b)、根据水声信号的遗传算法建模过程及原理，得出水声信号的状态方程；同时，将建模所得的数学模型转化为水声信号的观测方程；步骤(c)、利用粒子滤波对实测水声信号进行降噪，给出降噪前后的时序波形和相空间吸引子轨迹；步骤(d)、计算降噪前与降噪后的噪声强度，Lyapunov指数，关联维数和K₂熵；本发明利用粒子滤波对实测的水声信号进行降噪，给出降噪前后的信号时域波形及相空间吸引子轨迹；计算降噪前与降噪后的噪声强度，Lyapunov指数，关联维数，K₂熵等特征参数，为定量分析粒子滤波对实际水声信号的降噪效果提供依据。

Description

一种水声信号的粒子滤波降噪方法

技术领域

本发明属于水声信号处理技术领域，涉及一种水声信号的粒子滤波降噪方法。

背景技术

文献“非高斯背景条件下水声信号粒子滤波性能分析,火力与指挥控制,2014,Vol.39(4),p34-37,41”公开了一种非高斯背景条件下水声信号粒子滤波性能分析方法。该方法是利用假定的水声信号数学模型进行粒子滤波，但是需要说明的是，并未对实际水声信号进行粒子滤波降噪处理。

发明内容

本发明的目的在于克服上述现有方法的不足，提供一种水声信号的粒子滤波降噪方法，以四种不同类别的舰船辐射噪声作为实际水声信号，建立了基于遗传算法的实际水声信号的数学模型，提出了基于粒子滤波的水声信号降噪处理方法，并对实际舰船信号进行降噪处理。

为达到上述目的，本发明采用以下技术方案：

一种水声信号的粒子滤波降噪方法，包括以下步骤：

步骤(a)、利用遗传算法对实测水声信号建立数学模型；

步骤(b)、根据水声信号的遗传算法建模过程及原理，得出水声信号的状态方程；同时，将建模所得的数学模型转化为水声信号的观测方程；

步骤(c)、利用粒子滤波对实测水声信号进行降噪，给出降噪前后的时序波形和相空间吸引子轨迹，具体步骤如下：

1)、初始化k＝0

从初始分布p(x₀)中采样得到N个初始粒子

2)、重要性采样 

从重要性密度函数q(x_0:k|y_1:k)中采样产生粒子 $x_k^i~q\left(x_k^i\right|x_{0:k-1},y_{1:k}^i),i=\mathrm{1,2}...,N;$

3)、重要性权值 

当i＝2,…,N时， $w_k^i=\frac{p\left(y_k\right|x_k^i)p\left(x_k^i\right|x_{k-1}^i)}{q\left(x_k^i\right|x_{k-1}^i,y_k)}---\left(1\right)$

归一化权值 

$w_k^i=\frac{w_k^i}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^i}---\left(2\right)$

4)、重采样

从中根据重要性权值重新采样得到新的N个粒子并重新分配权值  $w_k^i=1/N;$

5)、状态估计均值

${\hat{x}}_k=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^i{x}_k^i---\left(3\right)$

步骤(d)、计算降噪前与降噪后的噪声强度，Lyapunov指数，关联维数和K₂熵。

进一步，所述步骤(a)建立数学模型具体包括以下步骤：

1)、使用Volterra模型作为预测模型，如下：

$\begin{array}{c}y=f(x_1,x_2,...,x_N)\\ =w_0+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_i{x}_i+\underset{i_1,i_2=1}{\overset{N}{\mathrm{\Σ}}}{w}_{i_1{i}_2}{x}_{i_1}{x}_{i_2}+...+\underset{i_1,i_2,...,i_M=1}{\overset{N}{\mathrm{\Σ}}}{w}_{i_1{i}_2...i_M}{x}_{i_1}{x}_{i_2}...x_{i_M}\end{array}---\left(4\right)$

对于N个变量的M阶的Volterra模型，一共有1+N+N^2+…+N^M项，在Volterra模型的基础上，选择一定的项来组成新的模型，组成模型的项的选择由 遗传算法完成；

2)、对每一项进行编码，对于模型的项采用以下通用公式表示：

$\underset{i=1}{\overset{M}{\mathrm{\Π}}}{x}_i^{n_i}{,n}_i\≥0---\left(5\right)$

式中，对于项用向量a＝(n₁,n₂,…,n_M)^T来表示；多个向量组成一个矩阵，表示一个模型，矩阵的每个列向量表示模型的一个项，矩阵的列数即为模型的项数；

3)、使用最小二乘法来计算模型系数，在得到模型系数后，使用评价数据来计算模型的预测误差，使用该误差作为模型好坏的评价；

交叉繁殖：将两个模型的项按照一定的规律分别排序，在两个序列中各随机选择一个交叉点，然后交换序列的后半部分，得到两个新的模型，作为交叉繁殖的后代；

变异：将随机生成模型的一个项随机替换掉模型原来的一个项，完成变异。

本发明的有益效果是：

本发明鉴于粒子滤波常用于非线性、非高斯动态系统最优估计，完全突破了Kalman滤波理论框架，提出了适用于滤除水声信号中噪声的新方法。利用粒子滤波对实测的水声信号进行降噪，给出降噪前后的信号时域波形及相空间吸引子轨迹；计算降噪前与降噪后的噪声强度，Lyapunov指数，关联维数，Kolmogorov熵等特征参数，定量分析粒子滤波对实际水声信号的降噪效果。

通过对含噪Lorenz模型和实测水声信号进行降噪处理，得到了较为满意的结果。结果，表明采用粒子滤波降噪方法对水声信号进行降噪处理是一种行之有效的方法，它在水声信号的混沌处理研究中具有重要作用。

附图说明

图1第一类水声信号的实测数据与一步预测数据；

a)实测数据与一步预测数据(2048点)；

b)实测数据与一步预测数据(300点)；

图2第二类水声信号的实测数据与一步预测数据；

a)实测数据与一步预测数据(2048点)；

b)实测数据与一步预测数据(300点)；

图3第三类水声信号的实测数据与一步预测数据；

a)实测数据与一步预测数据(2048点)；

b)实测数据与一步预测数据(300点)；

图4第四类水声信号的实测数据与一步预测数据；

a)实测数据与一步预测数据(2048点)；

b)实测数据与一步预测数据(300点)；

图5第一类水声信号粒子滤波降噪前后的时序波形和相空间吸引子轨迹；

a)降噪前的时序波形；b)降噪后的时序波形；c)降噪前的吸引子轨迹；d)降噪后的吸引子轨迹；

图6第二类水声信号粒子滤波降噪前后的时序波形和相空间吸引子轨迹；

a)降噪前的时序波形；b)降噪后的时序波形；c)降噪前的吸引子轨迹；d)降噪后的吸引子轨迹；

图7第三类水声信号粒子滤波降噪前后的时序波形和相空间吸引子轨迹；

a)降噪前的时序波形；b)降噪后的时序波形；c)降噪前的吸引子轨迹； d)降噪后的吸引子轨迹；

图8第四类水声信号粒子滤波降噪前后的时序波形和相空间吸引子轨迹；

a)降噪前的时序波形；b)降噪后的时序波形；c)降噪前的吸引子轨迹；d)降噪后的吸引子轨迹。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明方法完整的步骤：

(a)利用遗传算法对水声信号建模

选择一个基本模型框架，对每一个可能的模型进行遗传编码，建立模型评价标准，设计遗传算法操作，包括交叉、繁殖和变异，利用遗传算法建立了实测水声信号的预测模型。

利用遗传算法建立时间序列的预测模型，需要完成以下内容：

(1)选择一个基本模型框架，因为必须在一定的限制下，才能使用遗传算法完成模型的建立(实际上就是在基本模型框架基础上求解模型参数)。

(2)对每一个可能的模型进行遗传编码，就是用一串数据来表示一个模型。好的遗传编码策略可以很好的反应问题本身的规律，有助于遗传算法的求解。

(3)建立模型评价标准，对于每一个可能的模型，能够评价其好坏，遗传算法以此为依据搜索最好的模型。

(4)设计遗传操作，包括交叉繁殖，变异。遗传操作的设计要尽量考虑模型的特点，尽量好的保存优秀的模式。

在了解数据形成原理过程的情况下，可以根据数据的产生原理，建立准确 的预测模型(模型参数可以待解)。但是，要处理的是水声数据，并不了解水声数据形成的原理过程，不能建立准确的模型，这时一般使用可以通用的模型。这里，使用Volterra模型，如下

$\begin{array}{c}y=f(x_1,x_2,...,x_N)\\ =w_0+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_i{x}_i+\underset{i_1,i_2=1}{\overset{N}{\mathrm{\Σ}}}{w}_{i_1{i}_2}{x}_{i_1}{x}_{i_2}+...+\underset{i_1,i_2,...,i_M=1}{\overset{N}{\mathrm{\Σ}}}{w}_{i_1{i}_2...i_M}{x}_{i_1}{x}_{i_2}...x_{i_M}\end{array}---\left(1\right)$

对于N个变量的M阶Volterra模型，一共有1+N+N^2+…+N^M项。在这些项里面，并不是每一项都能反映时间序列的内部规律，一味的增加模型的项数，会导致模型预测性能下降。所以在Volterra模型的基础上，选择一定的项来组成新的模型，组成模型的项的选择由遗传算法完成。

为了用遗传算法求解模型(即选出合适的项作为模型的组成)，需要对每一项进行编码。对于模型的项，可以有以下的通用表示：

$\underset{i=1}{\overset{M}{\mathrm{\Π}}}{x}_i^{n_i},n_i\≥0---\left(2\right)$

式中， $\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{n}_i\≤M$

对于项可以用向量a＝(n₁,n₂,…,n_M)^T来表示。多个这样的向量组成一个矩阵，便表示一个模型，矩阵的每个列向量表示模型的一个项，矩阵的列数即为模型的项数。

对于时间序列的预测模型，一般使用预测精度作为模型好坏的评价标准。在本模型中，使用最小二乘来计算模型系数，所以模型的预测精度一定是随模型项数的增加而增加的，但实际上增加的多项式项并不一定反映数据的真实规律。为了解决这个问题，选择一部分数据作为训练数据，另外一部分数据最为 评价数据。对于某一模型，使用训练数据，通过最小二乘计算模型系数；在得到模型系数后，使用评价数据来计算模型的预测误差，使用该误差作为模型好坏的评价。

交叉繁殖：将两个模型的项按照一定的规律分别排序，在两个序列中各随机选择一个交叉点，然后交换序列的后半部分，这样便得到两个新的模型，作为交叉繁殖的后代。变异：将随机生成模型的一个项，然后随机替换掉模型原来的一个项，完成变异。

水声信号长度为2048点，每次预测6点，各点之间时间间隔为1即连续取，种群大小为60，繁衍代数为20，异化概率为0.03。下面给出四类水声信号的建模。

1)第一类水声信号建模

最优个体编码：

1 0 0 0 0 0 1 2 1 0 

0 1 0 0 0 1 2 0 1 1 

0 0 1 0 0 1 0 0 0 1 

0 0 0 1 0 0 0 1 1 1 

0 0 0 0 1 0 0 0 0 0 

第一行至第五行分别代表x_n-1，x_n-2，x_n-3，x_n-4，x_n-5；第一列至第十列共十列代表建模后数学表达式有十项组成；数值为1代表其本身，数值为2代表其本身的平方。例如：第七列为12000；表示建模后数学表达式的第七项为x_n-1x_n-2 ²。

通过最小二乘计算得到的系数为：

2.1291，-2.0046，1.2571，-0.5679，0.1492，-0.0025，-0.0603，-0.3920，0.4859， -0.0636。

通过最优个体编码和最小二乘计算得到的系数，其数学表达式为：

x_n＝2.1291x_n-1-2.0046x_n-2+1.2571x_n-3-0.5679x_n-4+0.1492x_n-5-0.0025x_n-2x_n-3

                                                                         (3) 

-0.0603x_n-1x_n-2 ²-0.3920x_n-1 ²x_n-4+0.4859x_n-1x_n-2x_n-4-0.0636x_n-2x_n-3x_n-4

第一类水声信号的时序波形及预测值如图1所示。

2)第二类水声信号建模

最优个体编码：

1 0 0 0 0 0 0 0 0 1 0 

0 1 0 0 1 1 0 1 0 1 0 

0 0 1 0 1 0 3 1 2 0 2 

0 0 0 1 0 1 0 1 1 0 0 

0 0 0 0 0 0 0 0 0 1 1 

第一行至第五行分别代表x_n-1，x_n-2，x_n-3，x_n-4，x_n-5；第一列至第十一列共十一列代表建模后数学表达式有十一项组成；数值为1代表其本身，数值为2代表其本身的平方，数值为3代表其本身的立方。例如：第七列为00300；表示建模后数学表达式的第七项为x_n-3 ³。

通过最小二乘计算得到的系数为：

1.8852，-0.6360，0.8201，-0.1708，0.0302，-0.0519，-0.3059，0.8315，-0.4886，-0.3684，0.1976

通过最优个体编码和最小二乘计算得到的系数，其数学表达式为：

x_n＝1.8852x_n-1-1.6360x_n-2+0.8201x_n-3-0.1708x_n-4+0.0302x_n-2x_n-3

-0.0519x_n-2x_n-4-0.3059x_n-3 ³+0.8315x_n-2x_n-3x_n-4-0.4886x_n-3 ²x_n-4        (4) 

-0.3684x_n-1x_n-2x_n-5+0.1976x_n-3 ²x_n-5

第二类水声信号的时序波形及预测值如图2所示。

3)第三类水声信号建模

最优个体编码：

1 0 0 0 0 0 0

0 1 0 0 1 1 1

0 0 1 0 0 1 0

0 0 0 1 0 0 1

0 0 0 0 0 0 0

0 0 0 0 1 1 1

第一行至第六行分别代表x_n-1，x_n-2，x_n-3，x_n-4，x_n-5，x_n-6；第一列至第七列共七列代表建模后数学表达式有七项组成；数值为1代表其本身。例如：第七列为010101；表示建模后数学表达式的第七项为x_n-2x_n-4x_n-6。

通过最小二乘计算得到的系数为：

1.7778，-1.6329，0.9791，-0.2537，0.0467，0.2955，-0.3004

通过最优个体编码和最小二乘计算得到的系数，其数学表达式为：

x_n＝1.7778x_n-1-1.6329x_n-2+0.9791x_n-3-0.2537x_n-4+0.0467x_n-2x_n-6

                                                             (5) 

+0.2955x_n-2x_n-3x_n-6-0.3004x_n-2x_n-4x_n-6

第三类水声信号的时序波形及预测值如图3所示。

4)第四类水声信号

最优个体编码：

1 0 0 0 0 0 0 1 0

0 1 0 0 0 0 1 0 0

0 0 1 0 0 0 1 1 2

0 0 0 1 0 0 1 0 0

0 0 0 0 1 0 0 1 1

0 0 0 0 0 1 0 0 0

第一行至第六行分别代表x_n-1，x_n-2，x_n-3，x_n-4，x_n-5，x_n-6；第一列至第九列共九列代表建模后数学表达式有九项组成；数值为1代表其本身，数值为2代表其本身的平方。例如：第七列为011100；表示建模后数学表达式的第七项为x_n-2x_n-3x_n-4。

通过最小二乘计算得到的系数为

1.5299，-0.5553，-0.0780，0.0520，0.0698，-0.0438，0.0848，-0.0419，-0.0304

通过最优个体编码和最小二乘计算得到的系数，其数学表达式为：

x_n＝1.5299x_n-1-0.5553x_n-2-0.078x_n-3+0.052x_n-4+0.0698x_n-5

                                                        (6) 

-0.0438x_n-6+0.0848x_n-2x_n-3x_n-4-0.0419x_n-1x_n-3x_n-5-0.0304x_n-3 ²x_n-5

第四类水声信号的时序波形及预测值如图4所示。

对实际水声信号建模的结果分析

均方根误差(RMSE：Root Mean Square Error)定义为：

$\mathrm{RMSE}=\sqrt{\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{(\hat{x}\left(n\right)-x\left(n\right))}^2}---\left(7\right)$

式中，x(n)表示n时刻的真实值，表示n时刻的预测值。

信误比(SER：Signal Error Ration)定义为：

$\mathrm{SER}=10\mathrm{lo}{g}_{10}\left[\frac{E{\left(x\left(n\right)\right)}^2}{E\left({(x-\left(n\right)-\hat{x}\left(n\right))}^2\right)}\right]---\left(8\right)$

水声信号预测误差与信误比表如表1所示。从表1的结果可以看出，采用遗传算法对实际水声信号进行建模是可行的，特别是当信号的时域波形起伏变化不大时，如第四类水声信号，它的建模效果更好，其均方根误差和信误比分别为0.0375和16.9dB。对于那些时域波形起伏变化较大的信号，如第三类水声 信号，其均方根误差和信误比分别为0.0884和9.3dB。

表1 水声信号预测误差与信误比

水声信号类别	均方根误差	信误比(dB)
			第一类水声信号	0.0425	13.3
第二类水声信号	0.0660	10.0
			第三类水声信号	0.0884	9.3
第四类水声信号	0.0375	16.9

通过遗传算法对四种水声信号分别进行建模，给出了相应的数学表达式，这将为水声信号的粒子滤波降噪奠定基础。

(b)实际水声信号的状态方程和观测方程

根据水声信号的遗传算法建模过程及原理，可以写出水声信号的状态方程；同时，建模所得的数学模型可以转化为水声信号的观测方程。实际水声信号的状态方程和观测方程如下：

状态方程： 

x_n＝f(x_n-1)  (9)

观测方程： 

y_n＝Hx_n+w_n  (10)

式中，x_n＝(x_n,x_n-1,...,x_n-M+1)^T， $x_n\left(m\right)=\left\{\begin{array}{cc}{x}_{n-1}(m-1)& m>=1\\ g\left(x_{n-1}\right)& m=0\end{array}\right.,$

H＝(1,0,...,0)_1×M，w_n＝randn，g(x_n-1)为水声信号建模的数学表达式。

例如，对于第一类信号，根据式(9)和式(10)，其状态方程和观测方程如下：

状态方程： 

$\begin{array}{c}\left[\begin{array}{c}{x}_n\left(0\right)\\ x_n\left(1\right)\\ x_n\left(2\right)\\ x_n\left(3\right)\\ x_n\left(4\right)\end{array}\right]=f\left(\begin{array}{c}\left[\begin{array}{c}{x}_{n-1}\left(0\right)\\ x_{n-1}\left(1\right)\\ x_{n-1}\left(2\right)\\ x_{n-1}\left(3\right)\\ x_{n-1}\left(4\right)\end{array}\right]\end{array}\right)\\ =\left[\begin{array}{c}2.1291X_{n-1}\left(0\right)-2.0046X_{n-1}\left(1\right)+1.2571X_{n-1}\left(2\right)-0.5679X_{n-1}\left(3\right)+\\ 0.1492X_{n-1}\left(4\right)-0.0025X_{n-1}\left(1\right)X_{n-1}\left(2\right)-0.0603x(n-1)X_{n-1}{\left(1\right)}^2-\\ 0.3920X_{n-1}{\left(0\right)}^2{X}_{n-1}\left(3\right)+0.4859X_{n-1}\left(0\right)X_{n-1}\left(1\right)X_{n-1}\left(3\right)-\\ 0.0636X_{n-1}\left(1\right)X_{n-1}\left(2\right)X_{n-1}\left(3\right)\\ X_{n-1}\left(0\right)\\ X_{n-1}\left(1\right)\\ X_{n-1}\left(2\right)\\ X_{n-1}\left(3\right)\end{array}\right]\end{array}$

观测方程： 

y_n＝x_n+w_n        (12) 

(c)利用粒子滤波对Lorenz信号进行降噪，给出降噪前后的时序波形；利用粒子滤波对实测水声信号进行降噪，给出降噪前后的时序波形和相空间吸引子轨迹。

1)初始化k＝0

从p(x₀)中采样粒子

2)重要性采样 

从重要性密度函数q(x_0:k|y_1:k)中采样产生粒子 $x_k^i~q\left(x_k^i\right|x_{0:k-1},y_{1:k}^i),i=\mathrm{1,2}...,N.$

3)重要性权值 

当i＝2,...,N时，

$w_k^i=\frac{p\left(y_k\right|x_k^i)p\left(x_k^i\right|x_{k-1}^i)}{q\left(x_k^i\right|x_{k-1}^i,y_k)}---\left(13\right)$

归一化权值 

$w_k^i=\frac{w_k^i}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^i}---\left(14\right)$

4)重采样

从中根据重要性权值重新采样得到新的N个粒子并重新分配权值  $w_k^i=1/N.$

5)状态估计均值 

${\hat{x}}_k=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^i{x}_k^i---\left(15\right)$

对于实测水声信号，为了验证粒子滤波的降噪性能，选取4种不同类别的舰船辐射噪声作为样本数据。各样本数据的采样率为20kHz，在降噪之前，各样本数据均先进行归一化处理。然后，利用粒子滤波对它们进行降噪处理。利用粒子滤波对它们进行降噪处理，处理结果分别如图5、图6、图7、图8，其中图a)、图b)分别为降噪前后的时域波形，为了便于比较，只画出了前200点的波形。图c)、图d)分别为降噪前后的相空间吸引子轨迹。从图中可以看出，降噪前后的信号时域波形及相空间吸引子轨迹都有明显的变化，粒子滤波对水声信号的噪声进行了较大程度的抑制，达到了降噪的目的。

(d)计算降噪前与降噪后的噪声强度，Lyapunov指数，关联维数，Kolmogorov熵等特征参数，并分析其降噪性能。

为了定量分析粒子滤波对实际水声信号的降噪效果，可以计算降噪前与降噪后的噪声强度，Lyapunov指数，关联维数，Kolmogorov熵等特征参数。

噪声强度： 

对于时间序列{x(n)}，n＝1，2，...，N，它的噪声强度(NI:Noise Intensity)可以近似表示为：

$\mathrm{NI}=\sqrt{\frac{1}{N-1}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{(x\left(n\right)-\stackrel{\‾}{x})}^2}---\left(13\right)$

式中：为时间序列的均值。

Lyapunov指数：

最大Lyapunov指数体现了相邻轨迹的平均分岔按指数规律增大的速度，它的定义如下所述：

假如是相空间中无限接近的两个点，它们之间的距离可以表示为 δ_Δn表示经过一段时间Δn后，起始于这两个点的两个轨迹之间的距离，即那么，最大Lyapunov指数λ由下式确定

δ_Δn≈δ₀e^λΔn          (14) 

式中，δ_Δn＜＜1，Δn＞＞1。

关联维数： 

若时间序列具有混沌特性，则其产生的奇怪吸引子具有分形结构，在几何上表现为分数维。分数维的定义有很多，关联维数是其中的一个。计算关联维数采用G-P算法.具体算法如下：

对于重构后相空间的点，选取适当的r(r＞0)，计算关联积分：

$C(m,r)=\frac{2}{M(M-1)}\underset{1\≤i\≤j\≤M}{\mathrm{\Σ}}\mathrm{\θ}(r-d_{\mathrm{ij}})---\left(15\right)$

当r增大时，关联积分C(m,r)与关联维数D有如下关系：

C(m,r)∝r^D           (16) 

因此，关联维数D定义为：

$D=\underset{r\→0}{\lim }\mathrm{InC}(m,r)/\mathrm{Inr}---\left(17\right)$

在计算时，用InC(m,r)-Inr图形的斜率估计D值。随着嵌入维数m的 增加，斜率趋于常数，这个常数就是关联维数。

Kolmogorov熵：

Kolmogoror熵(简称K熵)是刻划混沌系统的一个重要的量。在不同类型的动力学系统中，K熵的数值是不同的。K熵的数值可以用来区分规则运动、混沌运动和随机运动。在随机运动系统中，K熵是无界的；在规则运动系统中，K熵为零；在混沌运动系统中，K熵大于零，K熵越大，那么信息的损失速率越大,系统的混沌程度越大，或者说系统越复杂。

对于给定一个m维动力系统，把它的相空间划分成n个边长为r的m维立方体小盒子，并将时间轴划分成长度为τ的小段。x(t)是坐落在吸引子域中系统的一条轨线，p(i₁,i₂,…,i_n)表示x(τ)在盒子i₁中，x(2τ)在盒子i₂，…，x(nτ)在盒子i_n中的联合概率，Kolmogorov熵定义为：

$K=-\underset{\mathrm{\τ}\→0}{\lim }\underset{r\→0}{\lim }\underset{n\→\∞}{\lim }\frac{1}{\mathrm{n\τ}}\underset{(i_1,i_2,...,i_n)}{\mathrm{\Σ}}p(i_1,i_2,...i_n)\mathrm{Inp}(i_1,i_2,...,i_n)---\left(18\right)$

研究时间序列时，通常用K₂熵作为K熵的一个估计。K₂熵与关联积分C(r)存在如下关系：

$K_{2,m}\left(r\right)=\frac{1}{\mathrm{i\τ}}\mathrm{In}\frac{C_m\left(r\right)}{C_{m+i}\left(r\right)}---\left(19\right)$

对于给定的r，不断增加嵌入维数m，并求出相应的C_m(r)，由上式可得出K_2,m(r)，当其不随嵌入维数增加而改变时可用得到的K_2,m(r)作为K₂的估计值。

上述四类水声信号在降噪前与降噪后的噪声强度，Lyapunov指数，关联维数，K₂熵等特征参数值如表2所示。由表2可知，采用粒子滤波对水声信号进行降噪处理后，水声信号的上述特征参数同降噪前比较均有较明显的改善。主要体现在以下两方面，一是降噪后信号的噪声强度较显著降低，二是降噪后信号的关联维数，Lyapunov指数和K₂熵也同样获得较显著地减小。这一结果和理 论分析是一致的。

表2 四类水声信号粒子滤波降噪前后各特征参数比较

通过对含噪Lorenz模型和实测水声信号进行降噪处理，得到了较为满意的结果。结果，表明采用粒子滤波降噪方法对水声信号进行降噪处理是一种行之有效的方法，它在水声信号的混沌处理研究中具有重要作用。

Claims (2)

1.一种水声信号的粒子滤波降噪方法，其特征在于包括以下步骤：

步骤(a)、利用遗传算法对实测水声信号建立数学模型；

步骤(b)、根据水声信号的遗传算法建模过程及原理，得出水声信号的状态方程；同时，将建模所得的数学模型转化为水声信号的观测方程；

步骤(c)、利用粒子滤波对实测水声信号进行降噪，给出降噪前后的时序波形和相空间吸引子轨迹，具体步骤如下：

1)、初始化k＝0

从初始分布p(x₀)中采样得到N个初始粒子

2)、重要性采样

从重要性密度函数q(x_0:k|y_1:k)中采样产生粒子 $x_k^i~q\left(x_k^i\right|x_{0:k-1},y_{1:k}^i),i=\mathrm{1,2}...,N;$

3)、重要性权值

当i＝2,...,N时， $w_k^i=\frac{p\left(y_k\right|x_k^i)p\left(x_k^i\right|x_{k-1}^i)}{q\left(x_k^i\right|x_{k-1}^i,y_k)}---\left(1\right)$

归一化权值

$w_k^i=\frac{w_k^i}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^i}---\left(2\right)$

4)、重采样

从中根据重要性权值重新采样得到新的N个粒子并重新分配权值 $w_k^i=1/N;$

5)、状态估计均值

${\hat{x}}_k=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^i{x}_k^i---\left(3\right)$

步骤(d)、计算降噪前与降噪后的噪声强度，Lyapunov指数，关联维数和K₂熵。

2.根据权利要求1所述的水声信号的粒子滤波降噪方法，其特征在于：所述步骤(a)建立数学模型具体包括以下步骤：

1)、使用Volterra模型作为预测模型，如下：

$\begin{array}{c}y=f(x_1,x_2,...,x_N)\\ =w_0+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_i{x}_i+\underset{i_1,i_2=1}{\overset{N}{\mathrm{\Σ}}}{w}_{i_1{i}_2}{x}_{i_1{i}_2}+...+\underset{i_1,i_2,...,i_M=1}{\overset{N}{\mathrm{\Σ}}}{w}_{i_1{i}_2...i_M}{x}_{i_1}{x}_{i_2}...x_{i_M}\end{array}---\left(4\right)$

对于N个变量的M阶的Volterra模型，一共有1+N+N^2+…+N^M项，在Volterra模型的基础上，选择一定的项来组成新的模型，组成模型的项的选择由遗传算法完成；

2)、对每一项进行编码，对于模型的项采用以下通用公式表示：

$\underset{i=1}{\overset{M}{\mathrm{\Π}}}{x_i}^{n_i},n_i\≥0---\left(5\right)$

式中，对于项用向量a＝(n₁,n₂,…,n_M)^T来表示；多个向量组成一个矩阵，表示一个模型，矩阵的每个列向量表示模型的一个项，矩阵的列数即为模型的项数；

3)、使用最小二乘法来计算模型系数，在得到模型系数后，使用评价数据来计算模型的预测误差，使用该误差作为模型好坏的评价；

交叉繁殖：将两个模型的项按照一定的规律分别排序，在两个序列中各随机选择一个交叉点，然后交换序列的后半部分，得到两个新的模型，作为交叉繁殖的后代；

变异：将随机生成模型的一个项随机替换掉模型原来的一个项，完成变异。