0 引言

生物多样性监测作为生态监测的重要组成部分, 在研究物种资源与生态环境, 发展资源科学、生态学等方面具有十分重要的意义。其中, 鸟类作为一种对栖息地改变和环境变化反映极为敏感的环境指示生物, 是生物多样性评估与监测、生态环境影响评价的重要指标^[1]。

自然环境下连续声学监测数据中潜在鸣声片段检测与自动分段是开展相关生态学研究工作的前提。对于长时间跨度的连续监测数据, 人工检测方法在检测效率、可行性等方面都有极大的局限性。借鉴语音活动性检测 (Voice Activity Detection, VAD) 的研究成果, 文献[2]采用基于短时能量和短时过零率的双阈值算法完成声音事件检测 (Acoustic Event Detection, AED), 但是在自然环境下由于风、雨、人造飞行物、农业活动等因素影响, 监测数据中环境噪声无法避免, 此时该方法的检测性能将显著下降。文献[3-4]比较了多种VAD方法, 结果表明高斯混合模型 (Gaussian Mixture Model, GMM) 方法在噪声环境下具有比较稳健的检测性能。文献[5]针对野外实测的凤头距翅麦鸡 (Vanellus chilensis) 鸣声数据进行实验, 结果也验证了GMM检测方法的有效性。

为了识别检测出的潜在鸣声片段对应物种, 文献[6-7]通过计算潜在鸣声片段与特定物种鸣叫信号模板之间的谱图互相关 (Spectrogram Cross Correlation, SPCC) 来实现自动物种识别。该方法虽然实现简单, 但是其判决阈值需要在实际应用时经过多次反复调整以避免阈值偏高时漏检样本显著增多以及阈值偏低时误判结果急剧增加。文献[8]针对固定长度的鸣声片段谱图采用角度径向变换 (Angular Radial Transform, ART) 提取了图像域特征, 然而该方法由于需要将谱图与ART基函数进行卷积操作因此对环境噪声较为敏感。文献[9]提取了基于Mel尺度的小波包分解子带倒谱系数 (Wavelet Packet decomposition Subband Cepstral Coefficient, WPSCC) 特征并结合支持向量机 (Support Vector Machine, SVM) 进行分类识别; 文献[10-11]分别基于灰度共生矩阵 (Gray Level Co-occurrence Matrix, GLCM) 以及和差统计法 (Sum and Difference Histograms, SDH) 提取了鸟鸣片段谱图对应的纹理特征向量, 并开展了识别实验。需要指出的是, 上述三种方法在实验验证时采用的数据样本均为手工截取的长度固定为2 s的片段, 并未涉及潜在鸣声片段的自动检测与提取。考虑到面向野外大时间跨度连续声监测应用时, 人工截取片段方式不具有可行性, 因此针对自动检测方法输出且持续时长变化的潜在鸣声片段, 上述方法的性能还需要进一步研究和验证。

认知科学领域的最新研究成果表明, 多个不同频带的谱模式 (spectrum pattern) 信息对于鸟鸣的感知识别尤为重要^[12]。受此启发, 本文提出了一种基于Mel子带参数化特征的鸟鸣自动识别方法。在采用GMM完成连续声监测数据的自动分段基础上, 基于自回归 (AutoRegressive, AR) 模型分别建模鸣声片段经过Mel带通滤波后各个子带输出的能量序列, 得到能够描述不同种类鸟鸣信号谱模式信息的参数化特征，然后利用SVM分类器进行了分类识别。基于野外自然环境11种鸟鸣信号的实验结果表明, 本文方法识别性能明显优于现有基于纹理特征的方法。

1 基于GMM的潜在鸣声片段检测

给定一条鸟鸣记录时长为T, 则该条记录可记为s(t), 0≤t≤T。经过分帧、加窗及短时傅里叶变换 (Short-Time Fourier Transform, STFT) 后相应功率谱图可以表示成如下矩阵形式:

$ \boldsymbol{SP}=[\boldsymbol{sp}\left( 1 \right),\boldsymbol{sp}\left( 2 \right),\cdots ,\boldsymbol{sp}\left( L \right)] $

(1)

其中L表示记录中帧数目。第l帧 (1≤l≤L) 频域向量为：

$ \boldsymbol{sp}\left( l \right)={{[{{\left| S(0,l) \right|}^{2}},{{\left| S(1,l) \right|}^{2}},\cdots ,{{\left| S(\frac{N}{2}\;-1,l) \right|}^{2}}]}^{\rm T}} $

(2)

其中:T表示转置; S(k, l) 表示记录s(t) 在时频点 (k, l) 处的STFT结果, k(0≤k≤N/2-1) 为频率序号, N为每帧信号进行快速傅里叶变换 (Fast Fourier Transform, FFT) 的点数。考虑到记录中含有潜在鸣声的数据帧通常能量较大, 而环境噪声帧能量通常相对较小, 因此可以采用含有两个高斯分量的GMM建模记录中帧对数能量分布, 分别表征潜在鸣声帧以及环境噪声帧, 进而选取高似然率的数据帧组成候选声音事件, 完成潜在鸣声片段检测^[3-5]。在计算帧对数能量时, 由于很多鸟鸣信号相比语音而言频率较高, 且不同鸟类鸣声的频带范围也不同, 考虑到野外开展声学生态监测时通常需要考察多种不同鸟类, 因此普遍采用在一个较宽频带范围[f_L, f_H]内来分析数据帧。其中：f_L通常设为1 kHz用于部分滤除普遍存在的风噪以及人类活动影响；f_H则需要根据监测对象确定, 最高可达20 kHz^[5]。令N_L和N_H分别表示f_L和f_H对应的频率序号, 则第l帧的能量及其对数分别表示为:

$ e\left( l \right)=\sum\limits_{k={{N}_{L}}}^{{{N}_{H}}}{{{\left| S(k,l) \right|}^{2}}} $

(3)

$ le\left( l \right)={{\lg }}\left( e\left( l \right) \right) $

(4)

记录s(t) 对应的帧对数能量序列为:

$ \boldsymbol{le}={{[le\left( 1 \right),le\left( 2 \right),\cdots ,le\left( L \right)]}^{\rm T}} $

(5)

假设le中样本是由含有两个高斯分量的GMM生成, 则其概率密度函数可表示为:

$ \begin{align} & p\left( \boldsymbol{le};w,\mu ,\sigma \right)=\sum\limits_{m=1}^{2}{{{w}_{m}}N}\left( \boldsymbol{le}\left| {{\mu }_{m}},{{\sigma }_{m}} \right. \right) =\\ & \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_{m=1}^{2}{\frac{{{w}_{m}}}{{{\sigma }_{m}}\sqrt{2\pi }}}\exp \left( -\frac{{{\left( \boldsymbol{le}-{{\mu }_{m}} \right)}^{2}}}{2\sigma _{m}^{2}} \right) \\ \end{align} $

(6)

其中:对应于第m个高斯分量, μ_m为均值, σ_m²为方差。w_m表示高斯分量权重系数, 0≤w_m≤1, 且满足$\sum\limits_{m=1}^{2}{{{w}_{m}}}=1$。式 (6) 中参数集{w_m, μ_m, σ_m}, m=1, 2可以利用期望最大 (Expectation Maximization, EM) 算法获得相应的最大似然估计。

估计结果中的均值μ_m(m=1, 2) 依据大小可分别记为μ_H和μ_L, 二者各自对应的高斯分量可分别近似为潜在鸣声事件帧集合以及环境噪声帧集合的概率密度函数。此时对于le中的每一帧, 如果该帧属于潜在鸣声事件帧集合的后验概率大于属于环境噪声帧集合的后验概率, 则判决该帧归属某个潜在鸣声片段, 与该帧时间上连续且同样满足上述条件的其他帧也归属为该片段。利用上述检测方法输出的潜在鸣声片段集合记为D={AE₁, AE₂, …, AE_K}, 其中K为片段个数。

如图 1所示, 本文以一段野外实测的时长6 s的棕顶雀鹀鸣声数据为例来说明上述基于GMM的潜在鸣声片段检测与自动分段方法的有效性。其中：图 1(a)表示棕顶雀鹀鸣叫声的原始时频谱图；图 1(b)表示经过上述自动分段方法后提取的潜在鸣声片段的时频谱图。

2 Mel子带参数化特征提取

针对集合D={AE_i|1≤i≤K}中的各个候选片段, 首先利用符合人耳听觉特性的Mel滤波器组并行滤波, 然后针对各个子带输出的随时间变化的能量序列采用AR建模, 从而得到能够描述该事件的多个频带谱包络特性的参数化特征。

以第i个候选事件AE_i为例, 设其包含的帧数目为Q, 该事件相应的功率谱图矩阵可表示为SP _i=[sp _i(1), sp _i(2), …, sp _i(Q)], 其中每一帧的频域向量由式 (2) 定义。设置一个包含32个带通滤波器的Mel滤波器组, 其中第1个子带的下限Mel频率与第32个子带的上限Mel频率分别对应第1章中的物理频率f_L和f_H, 其转换关系由式 (7) 定义：

$ {{f}_{ \rm Mel}}=1127\times \ln \left( 1+{f}/{700}\; \right) $

(7)

Mel滤波器组中第j个带通滤波器离散化后在物理频率域的频率响应可以表示为:

$ \boldsymbol{h}\left( j \right)={{[{{h}_{j}}(0),{{h}_{j}}(1),\cdots ,{{h}_{j}}({N}/{2}\;-1)]}^{\rm T}};\\ \;\;\;\;\;\;\;\;\;j=1,2,\ldots ,32 $

(8)

需要注意的是, 式 (8) 中非零的频率响应仅存在于与[f_L, f_H]对应的[N_L, N_H]范围内。相应滤波器组的频率响应矩阵可表示为

$ \boldsymbol{H}={{[\boldsymbol{h}\left( 1 \right),\boldsymbol{h}\left( 2 \right),\cdots ,\boldsymbol{h}\left( 32 \right)]}^{\rm T}} $

(9)

候选事件AE_i经过并行滤波后的输出为：

$ \boldsymbol{Y}=\boldsymbol{H}\cdot \boldsymbol{S}{{\boldsymbol{P}}_{i}} $

(10)

令

$ \begin{align} & \boldsymbol{u}_{j}^{\rm T}=\underset{j}{\mathop{\rm row}}\,\left( \boldsymbol{Y} \right)-\rm {mean}\left( \underset{j}{\mathop{\rm row}}\,\left( \boldsymbol{Y} \right) \right) =\\ & \;\;\;\;\;\;\;[{{u}_{j}}\left( 1 \right),{{u}_{j}}\left( 2 \right),\cdots ,{{u}_{j}}\left( L \right)],j=1,2,\ldots ,32 \\ \end{align} $

(11)

其中: $\underset{j}{\mathop{\text{row}}}\,$ (·) 表示矩阵第j行, mean (·) 表示取均值。显然, u _j表示第j个Mel子带内信号能量随时间变化的信息, 而集合{ u ₁, u ₂, …, u ₃₂}则表征了事件AE_i在时频域上的谱模式信息。考虑到u _j具有近似谱包络特性, 因此AR模型是一种自然的合理选择。基于上述思路, u _j对应的M_j阶AR模型可表示为：

$ {{u}_{j}}\left( l \right)+a_{_{1}}^{(j)}{{u}_{j}}\left( l-1 \right)+\cdots +a_{_{{{M}_{j}}}}^{(j)}{{u}_{j}}\left( l-{{M}_{j}} \right)=z\left( l \right) $

(12)

其中:z(l) 为零均值白噪声序列, 模型参数a₁^(j), a₂^(j), …, a^(j)_Mj即为u_j对应的参数化特征。众所周知, 采用Levinson-Durbin算法可以快速高效地估计上述模型参数, 而由于不同Mel子带内信号能量随时间变化规律不同, 相应AR模型的阶次也不同, 本文基于常用的AIC (Akaike Information Criterion) 准则实现模型定阶。另外, 作者前期针对大量野外实测记录研究发现, 最高阶次设为10阶可以获得较为满意结果。综上所述, 用于表征一个潜在鸣声片段的参数化特征向量可以表示为：

$ \boldsymbol{v}={{[\boldsymbol{v}_{1}^{\rm T},\boldsymbol{v}_{2}^{\rm T},\cdots ,\boldsymbol{v}_{32}^{\rm T}]}^{\rm T}} $

(13)

$ {{\boldsymbol{v}}_{j}}={{[a_{_{1}}^{(j)},a_{_{2}}^{(j)},\cdots ,a_{_{M}}^{(j)}]}^{\rm T}},j=1,2,\ldots ,32 $

(14)

其中:M=10。当某个Mel子带输出序列u_j对应模型阶次M_j小于M时, 对应的v _j最后M-M_j个系数设置为零。

需要指出的是, 当前已有相关研究采用Mel系数特征, 通常处理过程是对每一帧数据提取Mel频率倒谱系数 (Mel Frequency Cepstrum Coefficient, MFCC) 特征, 然后借助隐马尔可夫模型 (Hidden Markov Model, HMM) 等分类器进行识别, 属于帧级的特征提取与识别^[13]。而本文方法是针对一个鸟鸣片段提取事件级特征, 在提取过程中采用Mel带通滤波器组用于近似模拟人耳听觉特性, 然后对每个鸟鸣事件的时频特性进行AR建模, 与已有研究明显不同。

3 基于SVM的分类识别

SVM是基于统计学理论的机器学习方法, 通常用于进行模式识别、分类以及回归分析^[14]。其基本思想是寻找一个满足分类要求的超平面, 使得训练集中的样本点与分类面距离尽可能远。SVM在解决小样本、非线性及高维度的模式识别问题中有独特的优势, 已广泛应用于图像识别、文本分类、语音识别等领域。本文在提取式 (13) 所示的参数化特征后, 采用基于“one-versus-one”策略^[15]构造的多分类SVM实现潜在鸣声片段的自动物种分类。

4 实验与分析 4.1 实验设计

实验采用的11种鸟鸣数据均来自Xeno-canto鸟声数据库^[16], 相关物种明细如表 1所示。该数据库内含的鸟鸣数据均在野外自然复杂声学环境下录取, 每条记录时长数十秒至数十分钟, 且对应物种均已标记。为了便于后续信号处理, 每条记录均统一转换为采样率32 kHz, 采样精度为16 bit。实验中谱图计算所用参数为帧长10 ms (320个采样点), 帧移5 ms, 海明 (Hamming) 窗加权, 512点FFT。在第1章所述自动分段方法处理后经过人工检查确定真实鸟鸣片段共2 762个。基于每个片段的起止时间信息, 帧移减小为2 ms重新计算该片段对应谱图以便获得更高的时间分辨率, 进而提取相应的Mel子带参数化特征 (分类实验中每个片段对应特征向量称为一个样本), 并采用SVM分类器进行物种识别实验。实验采用LIBSVM工具箱^[17], 选用径向基核函数 (Radial Basis Function, RBF), 并采用自动寻优方式设置惩罚因子c和核参数g, 最终对应参数为c=8, g=0.062 5。

需要指出的是, 鸟类鸣声 (vocalization) 包括鸣叫 (call) 和鸣唱 (song), 其中鸣唱相对较复杂, 一般由多个音节 (syllable) 组成, 并在一次鸣唱中多次重复出现。上述检测出的每个鸟鸣片段即对应于一次鸣叫或一个音节。通过考察野外实测数据发现, 同种鸟的鸣声片段持续时长变化较为剧烈, 如长嘴沼泽鹪鹩检测出的鸣声片段中最短持续时长为25 ms, 而最长可达595 ms; 歌带鹀鸣声片段中最短持续时长为50 ms, 最长则达400 ms。而文献[8-11]方法都是针对人工截取的固定长度的鸣声片段, 这对于野外长时间连续监测数据分析任务而言, 上述方法的应用将受到极大限制。

实验中, 采用和差统计法提取纹理特征^[11], 然后采用SVM分类器进行分类识别, 并与本文方法进行性能对比。两种方法的物种识别实验均进行了50次, 每次实验采用无放回方式从每类中随机抽取60%样本作为总体训练集, 余下每类40%样本共同组成总体测试集, 采用查准率precision、查全率recall和F1度量作为性能评价指标^[18]。定义分别如下：

$ precision=\frac{TP}{TP+FP} $

(15)

$ recall=\frac{TP}{TP+FN} $

(16)

$ F1=\frac{2\times precision\times recall}{precision\text{+}recall} $

(17)

其中: TP和FP分别表示测试集中该类所有正例被正确和错误分类的数目; 而FN则指测试集中该类所有负例被错误分类的数目。

4.2 结果分析

根据上述实验设计进行50次物种识别实验后, 两种方法的各类平均查准率、查全率以及F1度量结果分别如图 2~4所示。其中, 为简洁起见, 图中11种鸟依据表 1所示分别简记为C-C、M-M、C-P、G-T、S-P、S-A、A-H、C-B、B-C、H-M和P-C。

从图 2~3可以看出, 本文方法对每一类鸟鸣声识别的查准率和查全率均明显优于基于纹理特征的方法。其中, 本文方法对于各类鸟鸣的查准率结果中最高为99%, 最低为89%, 而纹理特征方法相应分别为89%和58%;查全率结果中本文方法最高为98%, 最低为92%, 纹理特征方法相应分别为90%和49%。尤其是对于每一类样本而言, 本文方法性能均明显优于纹理特征方法。上述结果也在图 4所示的查准率与查全率的综合性能评估指标F1度量中得到直观体现, 可以看出本文方法对于各类鸟鸣识别综合性能较为均匀, 且均不低于90%, 而纹理特征方法则差异度比较明显, 最高可达90%, 而最低则至仅55%, 表明其适用鸟鸣范围存在一定局限性。考虑到样本中有多类鸟鸣片段的持续时长变化较为剧烈, 上述实验结果也说明本文所提参数化特征对于不同时长的鸣声片段具有较好的普适性。此外, 针对实验中采用的野外实测记录数据, 本文方法并未事先进行去噪处理即可获得良好的分类结果; 而基于纹理特征的识别方法则还须通过音频增强算法抑制噪声影响^[10-11]。综上所述, 本文方法更适用于野外自然环境下连续监测数据的自动分析。

需要指出的是, 不同Mel子带数量会影响鸣声片段谱模式特性的建模效果。为了进一步研究Mel子带数量对本文方法性能的影响, 在不同Mel带通滤波器数量情况下分别进行50次物种识别实验, 每次实验中测试集和训练集产生方式与前述实验相同, 得到的各类平均F1度量结果如图 5所示。可以看出, 对于24和32个子带情况, 所有种类识别的综合评价指标F1均不低于90%, 而对于40个子带情况, 存在3类的综合评价指标低于90%。总体而言, 不同Mel子带数目对于本文方法性能的影响有限, 但是对比结果表明过多的子带数目对于不同物种鸟鸣信号而言并不总是利于识别性能的提升。

5 结语

本文提出了一种基于Mel子带参数化特征的鸟鸣自动识别方法。该方法首先利用基于GMM的AED处理过程提取潜在鸟鸣片段, 然后针对每个片段各个Mel子带输出能量序列进行AR建模, 提取能够描述不同种类鸟鸣信号时频谱模式特性的参数化特征, 并采用SVM进行分类识别。基于野外鸟鸣数据的实验结果表明, 本文方法对不同类型的鸟类物种具有良好的识别性能, 且明显优于现有纹理特征方法。但是需要注意的是, 考虑到自然环境下采集的鸣声数据中不可避免地会出现并发多声源情况, 即有多只和/或多类物种同时鸣叫或鸣唱, 此时直接使用本文提出的自动分段及识别方法并不合适, 原因在于此时检测出的潜在鸣声片段会是上述多个不同鸣声信号的合成结果, 相应提取的参数化特征显然不能正确表征相应物种的鸣声。因此未来的研究工作可以考虑利用麦克风阵列进行多声源分离后再使用本文提出的鸟鸣自动识别方法。