本发明属于声音信号处理技术领域，具体涉及一种基于高斯混合模型的宽吻海豚click信号检测方法。
背景技术：
时频分析技术是数字信号处理的研究内容之一，常用于信号处理、故障诊断和特征提取等领域，能够清楚地描述了信号频率随时间变化的关系。短时傅立叶变换基本思想是对信号加滑动时间窗(截断)，并对窗内信号做傅立叶变换，得到信号的频率域数据。加窗实质是用一个所谓的窗函数与原始的时域信号作乘积的过程，使得相乘后的信号似乎更好地满足傅立叶变换的周期性要求。
聚类是机器学习中“无监督学习”的一部分，常与数据挖掘、主成分分析和图像分割等学科关联研究。在“无监督学习”中，训练样本的标记信息是未知的，目标是通过对无标记训练样本的学习来揭示数据的内在性质，为进一步的数据分析提供基础。聚类是将物理或抽象对象的集合分成由类似的对象组成的多个类的过程。由聚类所生成的簇是一组数据对象的集合，这些对象与同一个簇中的对象彼此相似，与其他簇中的对象相异。传统的聚类分析算法有基于划分的方法、基于层次的方法、基于密度的方法、基于网格的方法和基于模型的方法。其中最常用的是基于划分的k-均值聚类方法和k-中心点聚类方法，还有基于模型的高斯混合模型聚类方法。k-均值聚类方法和k-中心点聚类方法的结果是每个数据点被划分到某一个聚类簇中，而高斯混合模型聚类方法的结果是这些数据点被划分到每个聚类簇的概率。因此，k-均值聚类方法和k-中心点聚类方法统称为硬聚类，高斯混合模型聚类方法被称为软聚类。
高斯混合模型就是用高斯概率密度函数精确地量化事物，将一个事物分解为若干的基于高斯概率密度函数形成的模型。简单来说，高斯混合模型就是若干个高斯概率密度函数的线性组合。如果一个图像的目标区域和背景区域相差比较大，且目标区域和背景区域在灰度上有一定的差距，则该图像的灰度直方图呈现双峰-谷形状，其中一个峰对应目标，另一个峰对应背景。而对于复杂图像而言，图像的直方图通常是多峰的，即对应多个高斯概率密度函数，就是高斯混合模型中的一个高斯分布。高斯混合模型聚类是通过对这些高斯分布的划分，将数据划分到不同的聚类簇。
二值化技术是数字图像处理领域最常用的技术之一，常用于图像预处理、边缘检测和图像分割等研究领域。图像二值化是指将图像灰度值在0到255范围内的像素点的灰度值设置为0或255，去掉中间的过渡像素，即整个图像呈现出非黑即白的视觉效果。图像二值化最常用的方法就是设置一个阈值t，用t将图像的像素点分为两部分，大于t的目标对象集和小于t的背景集。根据阈值选取方式的不同，图像二值化方法分为固定阈值方法和自适应阈值方法。顾名思义，固定阈值方法就是设置一个固定阈值不随图像的变化而改变，自适应阈值方法就是根据图像像素点的灰度分布自适应的计算一个最佳阈值。尽管到目前为止，二值化阈值的选取方法数以百计，但仍没有一个通用的方法对各种各样的图像都能达到最理像的效果。
技术实现要素：
本发明的目的在于提供一种基于高斯混合模型的宽吻海豚click信号检测方法，能够提取清楚完整的检测宽吻海豚click信号。
本发明采用的技术方案为，一种基于高斯混合模型的宽吻海豚click信号检测方法，具体按照以下步骤实施：
步骤1、采集海洋中的声音信号，利用时频分析技术将海洋中的声音信号的时间域信息与频率域信息联合分析，得到该信号的时频数据，并作出初始信号时频图，根据初始信号时频图得到能量强度；
步骤2、选取能量强度作为聚类的样本空间，使用高斯混合模型聚类方法对样本空间进行二聚类，得到两个聚类簇；
步骤3、根据两个聚类簇求自适应二值化的分割阈值，将步骤1中能量强度与分割阈值比较，根据比较结果显示二值化信号时频图，根据二值化信号时频图得到完整清晰的宽吻海豚click信号检测结果。
步骤1具体过程为：
步骤1.1、通过声音传感器采集海洋中的声音信号，使用matlab的audioread函数读入海洋中的声音信号；
步骤1.2、对读入海洋中的声音信号加滑动汉宁窗进行信号截断，得截断信号；
步骤1.3、对每个窗内截断信号进行傅里叶变换，得到关于时间和频率的函数；
步骤1.4、将根据关于时间和频率的函数作出初始信号时频图，该时频图以时间和频率为坐标轴，图上的每一点代表该时间该频率下海洋中的声音信号的能量强度。
步骤1.1海洋中的声音信号在matlab中以矩阵形式存在，该矩阵的行数代表读入该海洋中的声音信号时的采样样本点的个数，列数为该声音信号的音频通道的个数。
步骤1.2具体过程为：对读入海洋中的声音信号加汉宁窗，加窗后信号时域表达式＝原始信号的时域表达式*窗函数的时域表达式，汉宁窗每次移动其窗口大小的1/3，每次移动即为对读入海洋中的声音信号的截取；
其中，窗函数的时域表达式为：
步骤1.3傅里叶变换公式如下：
其中，f(w)表示频率为w时的能量强度，f(t)为截断信号。
步骤2具体过程为：
步骤2.1、以能量强度作为高斯混合模型聚类的训练集d＝{x1，x2，…，xm}，则该训练集d满足定义高斯混合模型定义：
其中，αi为混合系数，即组成高斯混合成分的每个高斯分布的比例，且满足
p(x|ui,∑i)满足高斯分布的定义：
步骤2.2、初始化高斯混合分布的模型参数{(αi,ui,∑i)|1≤i＜2}；
步骤2.3、设置循环终止条件：c＝1和循环终止条件：循环次数n；
步骤2.4、随机变量zj∈{1,2}表示生成训练集d的样本xj的高斯混合成分，其中j＝1,2,...,m，计算样本xj由第i个高斯分布生成的后验概率pm(zj＝i|xj)，将其记为γji(i＝1,2)，即
γji＝pm(zj＝i|xj)(1≤i＜2)(5)；
步骤2.5、对样本集d的随机变量的xj的后验概率pm(xj)求其极大似然估计，即最大化对数似然，如下：
步骤2.6、根据ll(d)计算新的均值向量ui'，新的协方差矩阵∑i'，新的混合系数αi'，将参数{(αi,ui,∑i)|1≤i≤2}更新为{(αi',ui',∑i')|1≤i≤2}；
步骤2.7、判是否满足循环终止条件；
若满足，则令聚类簇ci＝φ(1≤i≤2)；
若不满足，则循环次数c＝c+1，且跳到步骤2.4继续执行；
步骤2.8、针对随机变量xj，根据如下公式计算簇标记λj，将xj划入相应的聚类簇：cλj＝cλj∪{xj}；
步骤2.6具体过程为：
令
对混合系数αi而言，除了满足
两边同乘αi，对所有混合成分求和可知λ＝-m，则有：
步骤3具体过程为：
步骤3.1、根据高斯混合模型聚类方法得到该信号的两个聚类簇，对能量强度较小的聚类簇求其最大值max，对能量强度较大的聚类簇求其最小值min；
步骤3.2、对max和min求其平均值average，average作为自适应二值化的分割阈值；
步骤3.3、判断步骤1中得到的能量强度是否大于分割阈值：
若能量强度大于分割阈值，则该样本点的灰度值为255；
若能量强度不大于分割阈值，则该样本点的灰度值为0；
步骤3.4、显示二值化的信号时频图，根据二值化的信号时频图能够得到完整清晰的宽吻海豚click信号检测结果。
步骤3.1中能量强度较小的聚类簇为海洋噪声产生，能量强度较大的聚类簇为宽吻海豚的click信号。
本发明一种基于高斯混合模型的宽吻海豚click信号检测方法有益效果是：
1.本发明使用加窗傅里叶变换对采集的海洋信号进行时频分析，得到海洋声音信号的初始信号时频图，使得后续步骤实现简单，快捷；
2.本发明使用高斯混合模型的软聚类方法，对海洋声音信号进行二聚类，增加了检测的灵活性；
3.本发明的检测结果准确性高，完整性好，可以类比用于检测各种海豚的声音信号。
附图说明
图1是本发明一种基于高斯混合模型的宽吻海豚click信号检测方法流程图；
图2是本发明中时频分析方法流程图；
图3是本发明的高斯混合模型聚类方法流程图；
图4是本发明的自适应二值化方法流程图；
图5是本发明中海洋中的声音信号波形图；
图6是本发明中得到的初始信号时频图；
图7是本发明实施例中高斯混合模型的示例图；
图8是本发明实施例中检测结果示例图。
具体实施方式
下面结合附图和具体实施方式对本发明进行详细说明。
本发明一种基于高斯混合模型的宽吻海豚click信号检测方法，如图1所示，具体按照以下步骤实施：
步骤1、采集海洋中的声音信号，利用时频分析技术将海洋中的声音信号的时间域信息与频率域信息联合分析，得到该信号的时频数据，并作出初始信号时频图，根据初始信号时频图得到能量强度；
具体过程为：
步骤1.1、如图2所示，通过声音传感器采集海洋中的声音信号，使用matlab的audioread函数读入海洋中的声音信号；
海洋中的声音信号在matlab中以矩阵形式存在，该矩阵的行数代表读入该海洋中的声音信号时的采样样本点的个数，列数为该声音信号的音频通道的个数；
步骤1.2、对读入海洋中的声音信号加滑动汉宁窗进行信号截断，得截断信号；
对读入海洋中的声音信号加汉宁窗，加窗后信号时域表达式＝原始信号的时域表达式*窗函数的时域表达式，汉宁窗每次移动其窗口大小的1/3，每次移动即为对读入海洋中的声音信号的截取；
其中，窗函数的时域表达式为：
步骤1.3、对每个窗内截断信号进行傅里叶变换，得到关于时间和频率的函数；
傅里叶变换公式如下：
其中，f(w)表示频率为w时的能量强度，f(t)为截断信号；
步骤1.4、将根据关于时间和频率的函数作出初始信号时频图，该时频图以时间和频率为坐标轴，图上的每一点代表该时间该频率下海洋中的声音信号的能量强度。
步骤2、选取能量强度作为聚类的样本空间，使用高斯混合模型聚类方法对样本空间进行二聚类，得到两个聚类簇；
具体过程为：
步骤2.1、如图3所示，以能量强度作为高斯混合模型聚类的训练集d＝{x1，x2，…，xm}，则该训练集d满足定义高斯混合模型定义：
其中，αi为混合系数，即组成高斯混合成分的每个高斯分布的比例，且满足
p(x|ui,∑i)满足高斯分布的定义：
步骤2.2、初始化高斯混合分布的模型参数{(αi,ui,∑i)|1≤i＜2}；
步骤2.3、设置循环终止条件：c＝1和循环终止条件：循环次数n；
步骤2.4、随机变量zj∈{1,2}表示生成训练集d的样本xj的高斯混合成分，其中j＝1,2,...,m，计算样本xj由第i个高斯分布生成的后验概率pm(zj＝i|xj)，将其记为γji(i＝1,2)，即
γji＝pm(zj＝i|xj)(1≤i＜2)(5)；
步骤2.5、对样本集d的随机变量的xj的后验概率pm(xj)求其极大似然估计，即最大化对数似然，如下：
步骤2.6、根据ll(d)计算新的均值向量ui'，新的协方差矩阵∑i'，新的混合系数αi'，将参数{(αi,ui,∑i)|1≤i≤2}更新为{(αi',ui',∑i')|1≤i≤2}；
具体过程为：
令
对混合系数αi而言，除了满足
两边同乘αi，对所有混合成分求和可知λ＝-m，则有：
步骤2.7、判是否满足循环终止条件；
若满足，则令聚类簇ci＝φ(1≤i≤2)；
若不满足，则循环次数c＝c+1，且重新执行步骤2.4；
步骤2.8、针对随机变量xj，根据如下公式计算簇标记λj，将xj划入相应的聚类簇：cλj＝cλj∪{xj}；
步骤3、根据两个聚类簇求自适应二值化的分割阈值，将步骤1中能量强度与分割阈值比较，根据比较结果显示二值化信号时频图，根据二值化信号时频图得到完整清晰的宽吻海豚click信号检测结果；
具体过程为：
步骤3.1、如图4所示，根据高斯混合模型聚类方法得到该信号的两个聚类簇，对能量强度较小的聚类簇求其最大值max，对能量强度较大的聚类簇求其最小值min；
能量强度较小的聚类簇为海洋噪声产生，能量强度较大的聚类簇为宽吻海豚的click信号；
步骤3.2、对max和min求其平均值average，average作为自适应二值化的分割阈值；
步骤3.3、判断步骤1中得到的能量强度是否大于分割阈值：
若能量强度大于分割阈值，则该样本点的灰度值为255；
若能量强度不大于分割阈值，则该样本点的灰度值为0；
步骤3.4、显示二值化的信号时频图，根据二值化的信号时频图能够得到完整清晰的宽吻海豚click信号检测结果。
实施例
下面通过一个声音传感器采集到的声音信号作为实例来说明一种基于高斯混合模型聚类的宽吻海豚click信号的检测方法的实现过程。
步骤1.1、通过声音传感器采集海洋中的声音信号，使用matlab的audioread函数读入海洋中的声音信号；
海洋中的声音信号在matlab中以矩阵形式存在，该矩阵的行数代表读入该海洋中的声音信号时的采样样本点的个数，列数为该声音信号的音频通道的个数；
步骤1.2、对读入海洋中的声音信号加滑动汉宁窗进行信号截断，得截断信号；
对读入海洋中的声音信号加汉宁窗，加窗后信号时域表达式＝原始信号的时域表达式*窗函数的时域表达式，汉宁窗每次移动其窗口大小的1/3，每次移动即为对读入海洋中的声音信号的截取；
其中，窗函数的时域表达式为：
步骤1.3、对每个窗内截断信号进行傅里叶变换，得到关于时间和频率的函数；
傅里叶变换公式如下：
其中，f(w)表示频率为w时的能量强度，f(t)为截断信号；
步骤1.4、将根据关于时间和频率的函数作出初始信号时频图，该时频图以时间和频率为坐标轴，图上的每一点代表该时间该频率下海洋中的声音信号的能量强度。
读入海洋声音信号的波形图如图5所示，初始信号时频图如图6所示。
步骤2.1、以能量强度作为高斯混合模型聚类的训练集d＝{x1，x2，…，xm}，则该训练集d满足定义高斯混合模型定义：
其中，αi为混合系数，即组成高斯混合成分的每个高斯分布的比例，且满足
p(x|ui,∑i)满足高斯分布的定义：
步骤2.2、初始化高斯混合分布的模型参数{(αi,ui,∑i)|1≤i＜2}；
步骤2.3、设置循环终止条件：c＝1和循环终止条件：循环次数n；
步骤2.4、随机变量zj∈{1,2}表示生成训练集d的样本xj的高斯混合成分，其中j＝1,2,...,m，计算样本xj由第i个高斯分布生成的后验概率pm(zj＝i|xj)，将其记为γji(i＝1,2)，即
γji＝pm(zj＝i|xj)(1≤i＜2)(5)；
步骤2.5、对样本集d的随机变量的xj的后验概率pm(xj)求其极大似然估计，即最大化对数似然，如下：
步骤2.6、根据ll(d)计算新的均值向量ui'，新的协方差矩阵∑i'，新的混合系数αi'，将参数{(αi,ui,∑i)|1≤i≤2}更新为{(αi',ui',∑i')|1≤i≤2}；
具体过程为：
令
对混合系数αi而言，除了满足
两边同乘αi，对所有混合成分求和可知λ＝-m，则有：
步骤2.7、判是否满足循环终止条件；
若满足，则令聚类簇ci＝φ(1≤i≤2)；
若不满足，则循环次数c＝c+1，且重新执行步骤2.4；
步骤2.8、针对随机变量xj，根据如下公式计算簇标记λj，将xj划入相应的聚类簇：cλj＝cλj∪{xj}；
信号的时频数据划分为两个聚类簇，能量强度较高的宽吻海豚click信号聚类簇，和能量强度较弱的海洋噪声与混响聚类簇，有3个聚类中心的高斯混合模型的示例图如下图7所示。
步骤3.1、根据高斯混合模型聚类方法得到该信号的两个聚类簇，对能量强度较小的聚类簇求其最大值max，对能量强度较大的聚类簇求其最小值min；
能量强度较小的聚类簇为海洋噪声产生，能量强度较大的聚类簇为宽吻海豚的click信号；
步骤3.2、对max和min求其平均值average，average作为自适应二值化的分割阈值；
步骤3.3、判断步骤1中得到的能量强度是否大于分割阈值：
若能量强度大于分割阈值，则该样本点的灰度值为255；
若能量强度不大于分割阈值，则该样本点的灰度值为0；
步骤3.4、显示二值化的信号时频图，根据二值化的信号时频图能够得到完整清晰的宽吻海豚click信号检测结果。
由检测结果图8可以看出，下面的较短的横向线条是检测到宽吻海豚click信号。
通过上述方式，本发明一种基于高斯混合模型的宽吻海豚click信号检测方法，使用加窗傅里叶变换对采集的海洋信号进行时频分析，得到海洋声音信号的初始信号时频图，使得后续步骤实现简单，快捷；使用高斯混合模型的软聚类方法，对海洋声音信号进行二聚类，增加了检测的灵活性；检测结果准确性高，完整性好，可以类比用于检测各种海豚的声音信号。