[发明专利]一种水声信号的粒子滤波降噪方法

权利要求书

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}{x}_{i_2}+...+\underset{i_1,i_2,\mathrm{...},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模型的基础上，选择一定的项来组成新的模型，组成模型的项的选择由遗传算法完成；

其中：为模型系数，i₁,i₂…i_M均从1到N；

选择一定的项组成新模型方法如下：

在建立实际水声时间序列的预测模型过程中，将预测精度作为模型好坏的评价标准，将选择一部分数据作为训练数据，另外一部分数据作为评价数据；即对于某一模型，使用训练数据，通过最小二乘计算出模型系数；在得到模型系数后，使用评价数据来计算模型的预测误差，将该误差作为模型好坏的评价准则；当误差满足要求时，这时的项就是合适的项，从而组成新的模型；

组成模型的项的选择通过交叉繁殖和变异实现；

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

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

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

$\begin{array}{cc}\underset{i=1}{\overset{M}{\mathrm{\Π}}}{x_i}^{n_i}& n_i\≥0\end{array}---\left(2\right)$

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

M表示Volterra模型阶数；N表示模型变量数，即x_i，…，x_N；

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

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

实际水声信号的状态方程和观测方程如下：

状态方程：

x_n＝f(x_n-1) (3)

观测方程：

y_n＝Hx_n+w_n (4)

式中，x_n＝(x_n,x_n-1,...,x_n-M+1)^T，

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

状态方程：n时刻的状态x_n由n-1时刻的状态x_n-¹状态噪声组成；

观测方程：n时刻的观测值y_n由n时刻的状态x_n和观测噪声w_n组成；

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

1)、初始化k＝0

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

2)、重要性采样

从重要性密度函数q(x_0:k|y_1:k)中采样产生粒子

3)、重要性权值

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

归一化权值

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

4)、重采样

从中根据重要性权值重新采样得到新的N个粒子并重新分配权值

5)、状态估计均值

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

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

相空间吸引子轨迹是相空间无穷多个点的集合，这些点对应于系统的混沌状态；

降噪前的信号为实测水声信号；

首先对实测水声信号建立状态方程和观测方程；其次利用粒子滤波对其进行滤波；最后得出滤波后的信号，即降噪后时序波形；

对降噪前和降噪后的信号分别画出其相空间吸引子轨迹，即横坐标为x(n)，纵坐标为x(n+τ)，其中τ为时间延迟；

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

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

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

Lyapunov指数定义为：

$\mathrm{\λ}=\underset{n\→\mathrm{\∞}}{\lim }\frac{1}{n}\underset{i=0}{\overset{n-1}{\Σ}}In|\frac{df\left(x\right)}{x}{|}_{x=x_i}---\left(9\right)$

关联维数D定义为：

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

K₂熵作为K熵的估计，定义为：