CN102799892A - 一种mfcc水下目标特征提取和识别方法 - Google Patents

Abstract

一种MFCC水下目标特征提取和识别方法，包括步骤：1）获取数据序列x(n)；2）进行分帧处理得到x_i(n)；3）通过加窗操作得到y_i(n)；4）计算分帧加窗信号的单边功率谱密度p_i(l)；5）求解三角滤波器组的传递函数H_m(f)，通过三角滤波器组得到Q；6）进行对数变换得到E；7）分别进行t方向和f方向的Sobel算子和拉普拉斯算子3×3的模板块运算得到A、B和C；8）分别进行DCT变换得到特征集C_A,C_B,C_C，组合特征；9）通过基于EM算法的混合高斯模型的聚类分类器进行水下目标识别。本方法有利于提高水下目标的识别率。

Description

一种MFCC水下目标特征提取和识别方法

技术领域

本发明涉及水下目标的分类识别技术，尤其涉及一种新的MFCC水下目标特征提取和识别方法。

背景技术

海洋环境的复杂多变性使得目标特性很难定性描述，因此，水下目标的分类识别是水声信号处理的难题。

发明内容

技术问题：本发明所要解决的技术问题是：提供一种新的MFCC水下目标特征提取和识别方法，该方法利用人耳对声音的非线性辨别能力，通过图像锐化处理中的Sobel算子和拉普拉斯算子运算，获得基于图像处理的新的MFCC特征，并采用基于EM算法的混合高斯模型的聚类分类器，仿真结果具有高识别率。

技术方案：为解决上述技术问题，本发明的一种新的MFCC水下目标特征提取和识别方法：首先对水下目标信号（即舰船辐射噪声信号）分别进行分帧加窗的单边功率谱求取，求解基于Mel频率的三角滤波器的频域传递函数，将水下目标信号的功率谱与三角滤波器组相乘并进行对数变化，经过t方向和f方向的Sobel算子和拉普拉斯算子运算，最后进行DCT变化去除特征间的相关性，设计基于EM算法的混合高斯模型的聚类分类器，提高水下目标的识别率。

本发明的原理是利用了人耳对声音的非线性辨别能力，采用倒谱技术中的MFCC方法，将水下目标发声体和信道的频谱分开，从而提取发声体的特征，并利用图像锐化处理中的线性微分算子运算，提取有利的细微变化特征，设计基于EM算法的混合高斯模型的聚类分类器。

本发明方法的具体实现方式包括以下步骤：

第一步：获取舰船辐射噪声信号对应的数据序列x(n)：接收长度为N的实时采集数据或提取N点存储的现成数据作为待处理的数据序列x(n),n=0,1,…N-1，数据的采样率为f_s；

第二步：进行分帧处理得到x_i(n)：对第一步的数据序列x(n)进行分帧处理，得到每一帧数据序列x_i(n)，其中，

x_i(n)=x(n+id)           式（1）

在式（1）中，n=0,1,…L-1，i=0,1,…K-1，

表示向下取整运算；L表示x_i(n)的数据长度，d表示帧之间的重复数据数目，K表示帧的数目；

第三步：通过加窗操作得到y_i(n)：对每一帧数据序列x_i(n)(i＝0,1,…K-1)进行加窗处理，得到y_i(n)，其中，

y_i(n)=x_i(n)w(n)         式（2）

在式（2）中，i=0,1,…K-1,n=0,1,…L-1,w(n)=0.54-0.46cos(2πn/L)，；w(n)表示窗函数；

第四步：计算分帧加窗信号的单边功率谱密度p_i(l)：对数据序列y_i(n)求解单边功率谱密度p_i(l)，其中，

$p_i\left(l\right)={\left|\underset{n=0}{\overset{L-1}{\mathrm{\Σ}}}{y}_i\left(n\right)\exp \right\{-j2\mathrm{\πnl}/L\left\}\right|}^2$ 式（3）

在式（3）中，l=0,1,…L/2+1，||表示绝对值运算，j表示虚数单位，即i=0,1,2,…,K-1；

第五步：求解三角滤波器组的传递函数H_m(l)，通过三角滤波器组得到Q：求解三角滤波器组在频域的传递函数H_m(l)，从第四步的分帧的单边功率谱密度p_i(l)通过三角滤波器组得到Q，其中，

$Q(m,i)=\underset{l=0}{\overset{L/2+1}{\mathrm{\Σ}}}{H}_m\left(l\right)p_i\left(l\right)$ 式（4）

在式（4）中，H_m(l)表示三角滤波器的传递函数，l=0,1,…L/2+1；i=0,1,…K-1。

将范围在(0,f_s/2)的频率f均等分成L/2+1个频率点,按照式（5）计算Mel频率值

mel=2595log(1+f/700)        式（5）

将mel均等分成M+2个Mel频率值，mel的取值范围在(0,2595log(1+f_s/700))，依次对每个Mel频率值找f中对应最接近的频率值点，组成[f(0),f(2),…f(m),…f(M+1)]，按照式（6）计算H_m(f)，将H_m(f)中的f按照序号表示，得到H_m(l)，l=0,1,…L/2+1，

$H_m\left(f\right)=\left\{\begin{array}{c}0,f<f(m-1)\\ \frac{2(f-f(m-1))}{(f(m+1)-f(m-1))(f\left(m\right)-f(m-1))},f(m-1)\&le;f<f\left(m\right)\\ \frac{2(f(m+1)-f\left(m\right))}{(f(m+1)-f(m-1))(f(m+1)-f\left(m\right))},f\left(m\right)\&le;f<f(m+1)\\ 0,f\&GreaterEqual;f(m+1)\end{array}\right.$ 式（6）

在式（5）（6）中，m=1,2,…M，M表示三角滤波器的个数。

第六步：进行对数变换得到E：对Q进行对数变化得到E，其中，

E(m,i)=log(Q(m,i))           式（7）

在式（7）中，将E(m,i)小于1的置为1；m=1,2,…M，i＝0,1,…K-1。

第七步：分别进行t方向和f方向的Sobel算子和拉普拉斯算子3×3的模板块运算得到A，B和C：将E通过3×3的模板块运算得到A，B和C，其中，

$U(m,i)=\underset{T,F=-1}{\overset{1}{\mathrm{\&Sigma;}}}E(m+T,i+F)d_U(T,F)$ 式（8）

在式（8）中，U的第一行、最后一行、第一列和最后一列不做变换，U的其余行和列按照式（6）计算，d_U表示Sobel算子和拉普拉斯算子3×3的模板，U=A,B,C；m=1,2,…M，i=0,1,…K-1。

第八步：分别进行DCT变换得到特征集C_A,C_B,C_C，组合特征D：将第七步中得到的A，B和C进行离散余弦变换分别得到C_A,C_B,C_C，并对C_A,C_B,C_C按列取平均得到D_A,D_B,D_C，最后拼接成组合特征D，其中，

$C_U(n,i)=\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}U(m,i)\cos \{\mathrm{\&pi;}(m-0.5)n/M\}$ 式（9）

$D_U=\frac{1}{K}\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}{C}_U(n,i)$ 式（10）

在式（9）（10）中，n=1,2,…p，p表示每一帧数据的特征数目；U=A,B,C，i=0,1,…K-1。

第九步：通过基于EM算法的混合高斯模型的聚类分类器进行水下目标识别：按照第八步得到的特征数据集x={x¹,x²,...x^Num}作为基于EM算法的混合高斯模型的聚类分类器的输入，首先训练混合高斯模型的各参数，最后根据概率密度最高原则进行聚类，C={C¹,C²,…C^R},其中，混合高斯模型如式（11）（12）

$P=\underset{i=1}{\overset{R}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{P}_i\left(x\right|{\mathrm{\&upsi;}}_i)$ 式（11）

$P_i\left(x\right|{\mathrm{\&upsi;}}_i)={\left(2\mathrm{\&pi;}\right)}^{-\frac{d}{2}}\left|{{\mathrm{\&Sigma;}}_i}^{-\frac{1}{2}}\right|\exp (-\frac{1}{2}{(x-{\mathrm{\&mu;}}_i)}^T{\mathrm{\&Sigma;}}_i^{-1}(x-{\mathrm{\&mu;}}_i))$ 式（12）

在式（11）（12）中，μ_i表示第i-分支的期望，∑_i是第i-分支的协方差矩阵，R表示混合高斯模型的模型阶数。

设置初始值，依次按照式（13）（14）（15）（16）进行迭代，直到满足终止条件，停止迭代。

$w_{\mathrm{ij}}=\frac{\left|{\mathrm{\&Sigma;}}_i^{{\left(t\right)}^{-\frac{1}{2}}}\right|\exp (-\frac{1}{2}{(x^i-{\mathrm{\&mu;}}_j^{\left(t\right)})}^T{\mathrm{\&Sigma;}}_j^{\left(t\right)-1}(x^i-{\mathrm{\&mu;}}_j^{\left(t\right)}))}{\underset{k=1}{\overset{R}{\mathrm{\&Sigma;}}}\left|{\mathrm{\&Sigma;}}_k^{{\left(t\right)}^{-\frac{1}{2}}}\right|\exp (-\frac{1}{2}{(x^i-{\mathrm{\&mu;}}_k^{\left(t\right)})}^T{\mathrm{\&Sigma;}}_j^{\left(t\right)-1}(x^i-{\mathrm{\&mu;}}_k^{\left(t\right)}))}$ 式（13）

${\widehat{\mathrm{\&alpha;}}}_j=\frac{1}{\mathrm{Num}}\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}$ 式（14）

${\widehat{\mathrm{\&mu;}}}_j=\frac{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}{x}^i}{N{\widehat{\mathrm{\&alpha;}}}_j}$ 式（15）

${\widehat{\mathrm{\&Sigma;}}}_j=\frac{1}{\mathrm{Num}{\widehat{\mathrm{\&alpha;}}}_j}\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}{x}^i(x^i-{\widehat{\mathrm{\&mu;}}}_j){(x^i-{\widehat{\mathrm{\&mu;}}}_j)}^T$ 式（16）

在式（13）（14）（15）（16）中，j=1,2,…,R

有益效果：与现有技术相比，本发明具有以下有益效果：

1．将水下目标看成是发声体，海洋环境及其他的环境看成是信道，则人耳听到的水下目标声音使发声体和信道在时域上的卷积，频域上则是乘积，利用倒谱技术将发声体和信道分开，从而提取发声体相关的自身的物理特征，因此本发明具有明确的物理意义。

2．通过倒谱技术中的MFCC技术，利用人耳对声音的非线性辨别能力，将发声体和信道的频谱转化为Mel频率段的功率谱特征，有效细化了特征。又利用图像锐化处理，通过Sobel算子和拉普拉斯算计运算，以上两种算子是线性微分算子，微分即是求导变化，因此提取MFCC中的细微变化特征，本发明有利于提高水下目标的识别率。

附图说明

图1所示为本发明的流程图。

图2-1～2-3所示为实施例2中水下目标A经过t方向、f方向的Sobel算子及拉普拉斯算子运算的时频图；其中，

图2-1是t方向的Solbel算子时频图；

图2-2是f方向的Solbel算子时频图；

图2-3是拉普拉斯算子时频图；

图3-1～3-3所示为实施例2中水下目标B经过t方向、f方向的Sobel算子及拉普拉斯算子运算的时频图；其中，

图3-1是t方向的Solbel算子时频图；

图3-2是f方向的Solbel算子时频图；

图3-3是拉普拉斯算子时频图；

图4-1～4-3所示为实施例2中水下目标C经过t方向、f方向的Sobel算子及拉普拉斯算子运算的时频图；其中，

图4-1是t方向的Solbel算子时频图；

图4-2是f方向的Solbel算子时频图；

图4-3是拉普拉斯算子时频图。

具体实施方式

下面结合附图，对本发明的技术方案进行详细的说明。

如图1所示，本发明的一新的MFCC水下目标特征提取和识别方法，包括以下步骤：

第一步：获取数据序列x(n)：接收长度为N的实时采集数据或提取N点存储的现成数据作为待处理的数据序列x(n),n=0,1,…N-1，数据的采样率为f_s=22.5KHz；

第二步：进行分帧处理得到x_i(n)：对第一步的数据序列x(n)进行分帧处理，得到每一帧数据序列x_i(n)，其中，

x_i(n)=x(n+id)        式（1）

在式（1）中，n=0,1,…L-1，i＝0,1,…K-1，d和L的取值分别为256和512，根据

公式计算帧的数目；

第三步：通过加窗操作得到y_i(n)：对每一帧数据序列x_i(n)(i＝0,1,…K-1)进行加窗处理，得到y_i(n)，其中，

y_i(n)=x_i(n)w(n)      式（2）

在式（2）中，i＝0,1,…K-1,n=0,1,…L-1,其中w(n)为海明窗，窗函数为w(n)=0.54-0.46cos(2πn/L)，；

第四步：计算分帧加窗信号的单边功率谱密度p_i(l)：对数据序列y_i(n)求解单边功率谱密度p_i(l)，其中，

$p_i\left(l\right)={\left|\underset{n=0}{\overset{L-1}{\mathrm{\&Sigma;}}}{y}_i\left(n\right)\exp \right\{-j2\mathrm{\&pi;nl}/L\left\}\right|}^2$ 式（3）

在式（3）中，l=0,1,…L/2+1，||表示绝对值运算，j表示虚数单位，即

i=0,1,2,…,K-1；

在第四步中，离散傅里叶变换可用快速傅里叶变换。

第五步：求解三角滤波器组的传递函数H_m(l)，通过三角滤波器组得到Q：求解三角滤波器组在频域的传递函数H_m(l)，从第四步的分帧的单边功率谱密度p_i(l)通过三角滤波器组得到Q，其中，

$Q(m,i)=\underset{l=0}{\overset{L/2+1}{\mathrm{\&Sigma;}}}{H}_m\left(l\right)p_i\left(l\right)$ 式（4）

在式（4）中，H_m(l)表示三角滤波器的传递函数，l=0,1,…L/2+1；i=0,1,…K-1。

将范围在(0,f_s/2)的频率f均等分成L/2+1个频率点,按照式（5）计算Mel频率值

mel=2595log(1+f/700)        式（5）

将mel均等分成M+2个Mel频率值，mel的取值范围在(0,2595log(1+f_s/700))，依次对每个Mel频率值找f中对应最接近的频率值点，组成[f(0),f(2),…f(m),…f(M+1)]，按照式（6）计算H_m(f)，将H_m(f)中的f按照序号表示，得到H_m(l)，l=0,1,…L/2+1，

$H_m\left(f\right)=\left\{\begin{array}{c}0,f<f(m-1)\\ \frac{2(f-f(m-1))}{(f(m+1)-f(m-1))(f\left(m\right)-f(m-1))},f(m-1)\&le;f<f\left(m\right)\\ \frac{2(f(m+1)-f\left(m\right))}{(f(m+1)-f(m-1))(f(m+1)-f\left(m\right))},f\left(m\right)\&le;f<f(m+1)\\ 0,f\&GreaterEqual;f(m+1)\end{array}\right.$ 式（6）

在式（5）（6）中，m=1,2,…M，M表示三角滤波器的个数，取为36。

第六步：进行对数变换得到E：对Q进行对数变化得到E，其中，

E(m,i)=log(Q(m,i))          式（7）

在式（7）中，将E(m,i)小于1的置为1；m=1,2,…M，i＝0,1,…K-1。

第七步：分别进行t方向和f方向的Sobel算子和拉普拉斯算子3×3的模板块运算得到A，B和C：将E通过3×3的模板块运算得到A，B和C，其中，

$U(m,i)=\underset{T,F=-1}{\overset{1}{\mathrm{\&Sigma;}}}E(m+T,i+F)d_U(T,F)$ 式（8）

在式（8）中，U的第一行、最后一行、第一列和最后一列不做变换，U的其余行和列按照式（6）计算，d_U表示Sobel算子和拉普拉斯算子3×3的模板，U＝A，B，C；m＝1，2，…M，i＝0，1，…K-1。

其中t方向和f方向的Sobel算子和拉普拉斯算子3×3模板为：

$\left(\begin{array}{ccc}1& 2& 1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right)$ $\left(\begin{array}{ccc}1& 0& -1\\ 2& 0& -1\\ 1& 0& -1\end{array}\right)$ $\left(\begin{array}{ccc}0& -1& 0\\ -1& 4& -1\\ 0& -1& 0\end{array}\right)$

X方向Sobel算子    Y方向Sobel算子    拉普拉斯算子

第八步：分别进行DCT变换得到特征集C_A，C_B，C_C，组合特征D：将第七步中得到的A，B和C进行离散余弦变换分别得到C_A，C_B，C_C，并对C_A，C_B，C_C按列取平均得到D_A，D_B，D_C，最后拼接成组合特征D，其中，

$C_U(n,i)=\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}U(m,i)\cos \{\mathrm{\&pi;}(m-0.5)n/M\}$ 式(9)

$D_U=\frac{1}{K}\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}{C}_U(n,i)$ 式(10)

在式(9)(10)中，n＝1,2,…p，p表示每一帧数据的特征数目；U＝A，B，C，i＝0，1，…K-1。

第九步：通过基于EM算法的混合高斯模型的聚类分类器进行水下目标识别：按照第八步得到的特征数据集x＝{x¹，x²，...x^Num}，作为基于EM算法的混合高斯模型的聚类分类器的输入，首先训练混合高斯模型的各参数，最后根据概率密度最高原则进行聚类，C＝{C¹，C²，…C^R},其中，混合高斯模型如式(11)(12)

$P=\underset{i=1}{\overset{R}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{P}_i\left(x\right|{\mathrm{\&upsi;}}_i)$ 式(11)

$P_i\left(x\right|{\mathrm{\&upsi;}}_i)={\left(2\mathrm{\&pi;}\right)}^{-\frac{d}{2}}\left|{{\mathrm{\&Sigma;}}_i}^{-\frac{1}{2}}\right|\exp (-\frac{1}{2}{(x-{\mathrm{\&mu;}}_i)}^T{\mathrm{\&Sigma;}}_i^{-1}(x-{\mathrm{\&mu;}}_i))$ 式(12)

在式(11)(12)中，μ_i表示第i-分支的期望，∑_i是第i-分支的协方差矩阵，混合高斯模型的模型阶数为3。

μ_i和∑_i的初始值设置采用K-means方法，依次按照式(13)(14)(15)(16)进行迭代，直到满足前后两次迭代的参数误差平方和不超过0.01，停止迭代。

$w_{\mathrm{ij}}=\frac{\left|{\mathrm{\&Sigma;}}_i^{{\left(t\right)}^{-\frac{1}{2}}}\right|\exp (-\frac{1}{2}{(x^i-{\mathrm{\&mu;}}_j^{\left(t\right)})}^T{\mathrm{\&Sigma;}}_j^{\left(t\right)-1}(x^i-{\mathrm{\&mu;}}_j^{\left(t\right)}))}{\underset{k=1}{\overset{R}{\mathrm{\&Sigma;}}}\left|{\mathrm{\&Sigma;}}_k^{{\left(t\right)}^{-\frac{1}{2}}}\right|\exp (-\frac{1}{2}{(x^i-{\mathrm{\&mu;}}_k^{\left(t\right)})}^T{\mathrm{\&Sigma;}}_j^{\left(t\right)-1}(x^i-{\mathrm{\&mu;}}_k^{\left(t\right)}))}$ 式（13）

${\widehat{\mathrm{\&alpha;}}}_j=\frac{1}{\mathrm{Num}}\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}$ 式（14）

${\widehat{\mathrm{\&mu;}}}_j=\frac{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}{x}^i}{N{\widehat{\mathrm{\&alpha;}}}_j}$ 式（15）

${\widehat{\mathrm{\&Sigma;}}}_j=\frac{1}{\mathrm{Num}{\widehat{\mathrm{\&alpha;}}}_j}\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}{x}^i(x^i-{\widehat{\mathrm{\&mu;}}}_j){(x^i-{\widehat{\mathrm{\&mu;}}}_j)}^T$ 式（16）

在式（13）（14）（15）（16）中，j=1,2,3。

基于EM算法的混合高斯模型的聚类算法属于现有技术，例如《通信技术》2010年，11(43)的第150页至第152页中公开的内容。

具体来说：

实施例1

仿真模型：线谱用周期信号产生，环境噪声用零均值的高斯白噪声产生，

f_i分别取30Hz，47Hz，65Hz，89Hz，110Hz，247Hz，580Hz，信噪比

其中σ²为高斯白噪声的方差，采样频率f_s为10KHz，每类的样本点数为4000，每类的样本数为100，信噪比取为-14dB,得到线谱信号和线谱信号加噪声的两类目标的识别率见表1。

表1两类目标信号不同信噪比时的识别率

实施例2

三类水下目标A，B，C，采样频率为22.5KHz，每类样本的数目分别为48，在第一步分帧处理中，每一帧的长度为512，帧间的重复数据数目为256。三角滤波器组的数目为36，进行DCT变换的特征维数为24，则最后组合的特征维数为66。基于EM算法的混合高斯模型的聚类分类器，均值和方差的初始值设置采用K-means方法设置，分类器停止迭代的条件为前后两次的参数误差平方和不超过0.01。其中图2，3，4分别为水下目标A，B，C，经过三角滤波器组后通过线性微分算子（t方向和f方向的Sobel算子和拉普拉斯算子运算后的时频图）。三类水下目标的识别率见表2。

表2三类水下目标的识别率

从实施例1的结果可以看出，本发明特征提取和识别方法具有一定的抗噪能力，从实施例2中的结果可以看出，本发明特征提取和识别方法具有高识别率。

Claims (5)

1.一种MFCC水下目标特征提取和识别方法，其特征在于，步骤包括：

首先，对水下目标信号即舰船辐射噪声信号分别进行分帧加窗的单边功率谱求取；

接下来，求解基于Mel频率的三角滤波器的频域传递函数；

然后，将水下目标信号的功率谱与三角滤波器组相乘并进行对数变化，经过t方向和f方向的Sobel算子和拉普拉斯算子运算；

最后，进行DCT变化去除特征间的相关性，设计基于EM算法的混合高斯模型的聚类分类器。

2.按照权利要求1所述的一种新的MFCC水下目标特征提取和识别方法，其特征在于，本方法具体包括以下步骤：

第一步：获取舰船辐射噪声信号对应的数据序列x(n)：

接收长度为N的实时采集数据，或提取N点存储的现成数据作为待处理的数据序列x(n),n=0,1,…N-1，数据的采样率为f_s；

第二步：进行分帧处理得到x_i(n)：

对第一步的数据序列x(n)进行分帧处理，得到每一帧数据序列x_i(n)，其中，

x_i(n)=x(n+id)        式（1）

在式（1）中，n=0,1,…L-1，i=0,1,…K-1，

表示向下取整运算；L表示x_i(n)的数据长度，d表示帧之间的重复数据数目，K表示帧的数目；

第三步：通过加窗操作得到y_i(n)：

对每一帧数据序列x_i(n)(i＝0,1,…K-1)进行加窗处理，得到y_i(n)，其中，

y_i(n)=x_i(n)w(n)      式（2）

在式（2）中，i=0,1,…K-1,n=0,1,…L-1,w(n)=0.54-0.46cos(2πn/L)，；w(n)表示窗函数；

第四步：计算分帧加窗信号的单边功率谱密度p_i(l)：

对数据序列y_i(n)求解单边功率谱密度p_i(l)，其中，

$p_i\left(l\right)={\left|\underset{n=0}{\overset{L-1}{\mathrm{\&Sigma;}}}{y}_i\left(n\right)\exp \right\{-j2\mathrm{\&pi;nl}/l\left\}\right|}^2$ 式（3）

在式（3）中，l=0,1,…L/2+1，||表示绝对值运算，j表示虚数单位，即

i=0,1，2,…,K-1；

第五步：求解三角滤波器组的传递函数H_m(l)，通过三角滤波器组得到Q，其中，Q表示运算的结果，无特指含义：

求解三角滤波器组在频域的传递函数H_m(l)，从第四步的分帧的单边功率谱密度p_i(l)通过三角滤波器组得到Q，其中，

$Q(m,i)=\underset{l=0}{\overset{L/2+1}{\mathrm{\&Sigma;}}}{H}_m\left(l\right)p_i\left(l\right)$ 式（4）

在式（4）中，H_m(l)表示三角滤波器的传递函数，l=0,1,…L/2+1；i=0,1，…K-1；

将范围在(0,f_s/2)的频率f均等分成L/2+1个频率点,按照式（5）计算Mel频率值mel=2595log(1+f/700)        式（5）

将mel均等分成M+2个Mel频率值，mel的取值范围在(0,2595log(1+f_s/700))，依次对每个Mel频率值找f中对应最接近的频率值点，组成[f(0),f(2),…f(m),…f(M+1)]，按照式（6）计算H_m(f)，将H_m(f)中的f按照序号表示，得到H_m(l)，l=0,1,…L/2+1，

$H_m\left(f\right)=\left\{\begin{array}{c}0,f<f(m-1)\\ \frac{2(f-f(m-1))}{(f(m+1)-f(m-1))(f\left(m\right)-f(m-1))},f(m-1)\&le;f<f\left(m\right)\\ \frac{2(f(m+1)-f\left(m\right))}{(f(m+1)-f(m-1))(f(m+1)-f\left(m\right))},f\left(m\right)\&le;f<f(m+1)\\ 0,f\&GreaterEqual;f(m+1)\end{array}\right.$ 式（6）

在式（5）（6）中，m=1,2,…M，M表示三角滤波器的个数；

第六步：进行对数变换得到E：

对Q进行对数变化得到E，其中，

E(m,i)=log(Q(m,i))        式（7）

在式（7）中，将E(m,i)小于1的置为1；m=1,2,…M，i＝0,1,…K-1；

第七步：分别进行t方向和f方向的Sobel算子和拉普拉斯算子3×3的模板块运算得到A，B和C，其中A，B和C是运算得到的结果，无特指含义：

将E通过3×3的模板块运算得到A，B和C，其中，

$U(m,i)=\underset{T,F=-1}{\overset{1}{\mathrm{\&Sigma;}}}E(m+T,i+F)d_U(T,F)$ 式（8）

在式（8）中，U的第一行、最后一行、第一列和最后一列不做变换，U的其余行和列按照式（6）计算，d_U表示Sobel算子和拉普拉斯算子3×3的模板，U=A,B,C；m=1,2,…M，i=0,1,…K-1；

第八步：分别进行DCT变换得到特征集C_A,C_B,C_C，组合特征D：

将第七步中得到的A，B和C进行离散余弦变换分别得到C_A,C_B,C_C，并对C_A,C_B,C_C按列取平均得到D_A,D_B,D_C，最后拼接成组合特征D，其中，

$C_U(n,i)=\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}U(m,i)\cos \{\mathrm{\&pi;}(m-0.5)n/M\}$ 式（9）

$D_U=\frac{1}{K}\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}{C}_U(n,i)$ 式（10）

在式（9）（10）中，n=1,2,…p，p表示每一帧数据的特征数目；U=A,B,C，i=0,1,…K-1；

第九步：通过基于EM算法的混合高斯模型的聚类分类器进行水下目标识别。

3.按照权利要求1所述的一种MFCC水下目标特征提取和识别方法，其特征在于，所述第九步中，按照第八步得到的特征数据集x={x¹,x²,...x^Num}作为基于EM算法的混合高斯模型的聚类分类器的输入；

首先训练混合高斯模型的各参数，然后根据概率密度最高原则进行聚类，C={C¹,C²,…C^R},其中，混合高斯模型如式（11）、（12）

$P=\underset{i=1}{\overset{R}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{P}_i\left(x\right|{\mathrm{\&upsi;}}_i)$ 式（11）

$P_i\left(x\right|{\mathrm{\&upsi;}}_i)={\left(2\mathrm{\&pi;}\right)}^{-\frac{d}{2}}\left|{{\mathrm{\&Sigma;}}_i}^{-\frac{1}{2}}\right|\exp (-\frac{1}{2}{(x-{\mathrm{\&mu;}}_i)}^T{\mathrm{\&Sigma;}}_i^{-1}(x-{\mathrm{\&mu;}}_i))$ 式（12）

在式（11）（12）中，μ_i表示第i-分支的期望，∑_i是第i-分支的协方差矩阵，R表示混合高斯模型的模型阶数；

设置初始值，依次按照式（13）（14）（15）（16）进行迭代，直到满足终止条件，停止迭代。

$w_{\mathrm{ij}}=\frac{\left|{\mathrm{\&Sigma;}}_i^{{\left(t\right)}^{-\frac{1}{2}}}\right|\exp (-\frac{1}{2}{(x^i-{\mathrm{\&mu;}}_j^{\left(t\right)})}^T{\mathrm{\&Sigma;}}_j^{\left(t\right)-1}(x^i-{\mathrm{\&mu;}}_j^{\left(t\right)}))}{\underset{k=1}{\overset{R}{\mathrm{\&Sigma;}}}\left|{\mathrm{\&Sigma;}}_k^{{\left(t\right)}^{-\frac{1}{2}}}\right|\exp (-\frac{1}{2}{(x^i-{\mathrm{\&mu;}}_k^{\left(t\right)})}^T{\mathrm{\&Sigma;}}_j^{\left(t\right)-1}(x^i-{\mathrm{\&mu;}}_k^{\left(t\right)}))}$ 式（13）

${\widehat{\mathrm{\&alpha;}}}_j=\frac{1}{\mathrm{Num}}\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}$ 式（14）

${\widehat{\mathrm{\&mu;}}}_j=\frac{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}{x}^i}{N{\widehat{\mathrm{\&alpha;}}}_j}$ 式（15）

${\widehat{\mathrm{\&Sigma;}}}_j=\frac{1}{\mathrm{Num}{\widehat{\mathrm{\&alpha;}}}_j}\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}{x}^i(x^i-{\widehat{\mathrm{\&mu;}}}_j){(x^i-{\widehat{\mathrm{\&mu;}}}_j)}^T$ 式（16）

在式（13）（14）（15）（16）中，j=1,2,…,R。

4.按照权利要求1所述的一种MFCC水下目标特征提取和识别方法，其特征在于，所述的第七步中，Sobel算子和拉普拉斯算子3×3的模板为

$\left(\begin{array}{ccc}(m-1,i-1)& (m-1,i)& (m-1,i+1)\\ (m,i-1)& (m,i)& (m,i+1)\\ (m+1,i-1)& (m+1,i)& (m+1,i+1)\end{array}\right).$

5.按照权利要求1所述的一种MFCC水下目标特征提取和识别计方法，其特征在于，所述的第九步中，初始值设置采用K-means设置方法。

Images