基于声学指数的神农架国家公园声音多样性动态变化

罗 丽,侯亚男,杨敬元,于新文,高灵旺,欧阳萱,杨铭伦,高家军,郭安琪,刘昱坤

(1.中国林业科学研究院资源信息研究所,北京 100091;2.中国农业大学植物保护学院,北京 100083;3.成都东软学院计算机与软件学院,四川 成都 611844;4.神农架国家公园管理局科学研究院,湖北 神农架林区 442421;5.国家林业和草原局林业遥感与信息技术重点实验室,北京 100091)

生物多样性对人类的生存和发展具有十分重要的意义,2021年10月的《昆明宣言》[1]指出,过去10年的生物多样性战略计划取得了一定的进展,但这些进展不足以实现生物多样性目标[2]。因此,在全球气候变暖、自然景观日益破碎化的情况下,开发有效的大规模生物多样性监测方法是保护和研究生物多样性的重要环节[3]。

声音是生物之间交流的重要手段,包含着重要的生态信息,不仅可以作为物种识别的特征,也可以传达有关物种的行为、丰度、位置、体型等线索[4]。快速声学调查(rapid acoustic survey,RAS)[5]是近年来提出的以一种非侵入性的、具有成本效益的快速生物多样性评估方法,通过对大规模野外声音的采集和计算,分析和推断大时空尺度下的物种多样性。目前,生态声学领域对于声音数据的分析主要分为两个方向:一是运用或者改良现有的语音识别、机器学习等领域的技术达到识别、分类动物声音的目的[6-7],与传统生物多样性调查方法类似,此种研究方法所使用的数据往往是非连续的、片段式的;二是基于声景生态学(soundscape ecology)[8]方法,使用声学指数(acoustic indices)[9]或者根据实际的研究对象和数据设计新的声学指数来评估生物多样性及其动态变化情况,运用此方法的数据往往是长时间(>24 h,甚至是大 于1月)的连续音频数据。因此,基于声学指数的方法,能够反映生物多样性连续的变化特征。此外,声学指数通常被认为是快速评估生物多样性水平的生态指标[10-11]。声学指数是根据声音的振幅强度、时间/频率或振幅的复杂程度、声景的异质性来评估声音群落的强度、相对丰度、丰富度、异质性以及人为干扰对声学群落的影响[12],大致可以分为α声学指数和β声学指数两类,前者用于评估声景或动物群落的生物多样性和复杂性,后者用于评估不同动物群落或同一群落不同时间段内声音多样性的差异。目前大多数研究集中于探究声学特征与鸟类[13-14]、两栖动物[15-16]和昆虫[17-18]多样性的相关性,鲜有研究利用声音来探索生物多样性的整体变化特征。

2018—2019年中国林业科学研究院研究团队在神农架国家公园大龙潭部署了野外声音采集设备,基于声学指数对长时间序列的声音数据进行了声景分析和鸟类多样性检测[19-21],发现声学复杂性指数(acoustic complexity index, ACI)、生物声学指数(bioacoustic index, BI)和归一化声景差异指数(normalized difference soundscape index, NDSI)与该地区鸟类多样性的相关性较好,能反映出当地声景变化。但该研究仅局限在大龙潭小区域内,难以反映大区域尺度的生物多样性变化情况。本研究在此基础上,在更大的空间尺度上部署了相同的声音采集设备,利用物种鸣叫数量较多、活跃程度较大的夏季的声音数据,分析上述3个指数在神农架国家公园内的动态变化规律,一是为了继续探索声学指数对目标区域内动物鸣声多样性的响应,二是为了探究声学指数在神农架国家公园夏季的动态变化特征,旨在为当地的动物多样性调查和监测提供一种适用性的技术途径,同时为国内其他保护区开展生物多样性的量化研究提供参照。

1 材料与方法

1.1 研究区概况

研究区域位于神农架国家公园内。声音采集设备部署范围为110°00′～110°41E, 31°41′～31°53″N,覆盖面积约200 km2。研究区内多高山,气候属于北亚热带和暖温带的过渡区域,受山脉影响,气温垂直分带明显,年降水量800～2 500 mm,降水量随海拔升高而增加。

1.2 采样方法

1)采样点设置。相较于在大龙潭区域的研究[19-20],本研究共设置了9个采样点,覆盖的空间尺度更大,分布在区域内的坡底、坡中、坡顶、沟谷、溪流等不同生境类型内、采样点之间最小距离为2 150.51 m,最大距离为 38 601.83 m,采样点覆盖的海拔梯度为1 396～3 011 m,如表1所示。

表1 各采样点基本信息

2)数据采集及预处理。录音设备安装在距地面适宜高度的树干上,使用外部电池可以连续工作2个月,采用立体声录音,采样率为24 kHz,量化位数为16位的WAV格式自动保存在存储卡上。图1为9个采样点中2号观音洞采样点的采样设备及周围环境。采样设备连续工作,全天候不间断,声音数据每2 h存储1次,一天共获取12条2 h的声音数据。本次实验选取了2021年5月21日—7月17日样点的声音数据,在每个采样点收集了740条录音数据。其中,采样点7数据丢失,其余采样点数据完整,故采用8个采样点进行数据处理和分析,共5 920条,总时长为1 480 h,大小为952 G。

图1 2号观音洞采样点及其周围环境Fig. 1 Sampling sites of No. 2 and its surrounding environment

每个原始声音文件时长2 h,为避免数据冗余,使用Python脚本在音频每10 min截取开始前2 min生成小的声音文件,经过处理获得144个2 min的声音文件,以表示一整天的声音环境。本研究在每个采样点截取了7 488条(144×52 d)2 min的录音数据,共计59 904条(7 488×8个采样点),大小为718 GB。

使用Adobe Audition 2019软件对声音文件进行人工监听和频谱图扫描,确定采集到的声音主要分布在0～8 kHz范围内,其中(0,2] kHz的声音主要来自录音设备运行噪声,(2,6] kHz范围内的主要由鸟类鸣声构成,兼具少数其他动物声音,(6,8] kHz范围内的声音主要出现在夜间,来源于少数夜间活动的鸟类或昆虫。由于在野外采集不可避免会存在设备底噪,为了增强生物声音,对声音数据中设备底噪进行样本捕捉后,使用谱减法分别对每条声音数据进行降噪处理。

1.3 生态声学指数计算

指数计算均使用一个经过实验验证的python脚本实现(https://github.com/patricguyot/Acoustic_Indices),与以R语言中seewave/soundecology包的计算结果等效。前期研究发现,ACI、BI、NDSI与当地鸟类物种丰富度显著相关,因此,使用了ACI、BI和NDSI指数进行分析。其中,ACI、BI指数是基于同一个声谱图计算的,该声谱图计算为FFT的平方幅度,窗口和跳跃大小分别为512帧和256帧,指数的具体计算公式见文献[22-24]。

声学复杂度指数(ACI)[22]通过计算音频中强度的变异性来直接量化景观的复杂性。受背景噪声影响较小,可以用来反映复杂声音的频率和振幅的快速变化,数值越大,表明物种鸣叫数量越多、活跃程度越大。生物声学指数(BI)[23]用来衡量群落水平的动物声音多样性,计算频谱图中一定频率范围内、某一阈值以上的声强。BI用来反映物种群落的丰富度,数值越大,表明物种越多,低值表明很少或没有声学活动。归一化声景差异指数(NDSI)[24]在景观尺度上描述人类活动对声景的干扰,取值范围为[-1, 1],NDSI值高代表生物声在该生境中占主导地位,反之代表人类活动占主导地位。

1.4 数据分析

利用Mann-Kendall检验法对长时间序列的ACI、BI和NDSI指数进行了趋势分析和突变性检验。Mann-Kendall是一种时间尺度下的变化趋势性检验,属于一种非参数检验方法,可以统计评估变量随时间变化是否有单调上升或下降的趋势,单调上升(下降)意味着变量随时间序列增加(减少)。Uk为正序统计量,它是按时间序列X顺序X1,X2,…,Xn计算出的统计量序列,计算公式如下:

(1)

声学指数的时间序列作图与箱线图的制作在Origin 2019b中进行,声学指数的Mann-Kendall检验在Matlab 2016b中实现。

2 结果与分析

2.1 神农架夏季声学指数的动态变化特征

为了探索ACI、BI和NDSI指数在夏季对声学活动的动态变化响应,对2021年5月21日—7月17日采集的声学指数进行了连续性动态分析,为方便比较,各指数值均已进行Min-Max归一化处理,统一映射到[0, 1]区间上,结果如图2所示。

图2 各采样点夏季ACI、BI和NDSI指数的动态变化 Fig. 2 Changes of ACI, BI and NDSI at each sampling site in summer of Shengnongjia National Park

从图2可以看出,1号采样点的ACI指数值变化呈现多峰变化,BI和NDSI指数值变化趋势相似,生物声变化呈现先增加后下降的变化趋势。2号采样点ACI指数保持较小值且变化趋势不大,BI指数在7月初出现一个峰值,NDSI指数在7月初出现短暂的下降趋势,说明在这段时间内2～8 kHz出现的声学活动较少。3号采样点ACI与NDSI指数变化趋势相似,在7月7日ACI出现低的峰值而NDSI出现了高的峰值。9号采样点ACI指数5月至6月中旬变化趋势不大,6月末至7月呈下降趋势,BI、NDSI指数在5月至6月中旬呈上升趋势,6月末到7月呈下降趋势。

总的来说,ACI指数值在8个采样点变化的总体趋势不明显;BI指数值随时间变化大致呈现双峰变化,5月末开始上升,到6月中旬开始下降,7月初开始上升,7月中旬下降,说明采样点在2～8 kHz的频率范围内,声强变化呈现先上升后下降再上升的趋势;NDSI指数对生物声和人为活动声较敏感,变化幅度较大,大致呈现3峰/4峰变化趋势,5月21日至6月初、6月初至6月末、7月初至7月17日呈先上升后下降趋势。

2.2 声学指数的突变性检验及其在突变日期的日变化特征

图3 ACI、BI和NDSI指数的Mann-Kendall突变性检验Fig. 3 Mann-Kendall mutation tests for ACI, BI and NDSI

4号采样点ACI指数、BI指数和NDSI指数均发生了突变,且NDSI指数在该采样点发生了多次突变。5号采样点ACI指数在5月21—27日期间发生多次突变,在突变前显著上升,突变后显著下降,BI指数在6月8日15:00之前和7月12日0:00之后发生多次突变,NDSI指数发生突变的时间为6月16日2:00(X= 125.5),6号采样点BI指数与5号采样点BI指数变化相似。

综上所述,指数值的变化趋势较大,均有显著上升和显著下降的趋势,说明区域内动物鸣叫有显著性的变化,受温度的影响,动物的鸣叫行为会出现显著增加和显著减少的变化。

根据突变性检验结果,对声学指数的突变日期进行深入分析,探究造成突变的原因。ACI指数在4号、5号、8号采样点突变日期的日变化趋势如图4a所示,其中在5号采样点多次发生突变,选取了5月22、23、26日作日变化分析;BI指数与NDSI指数的突变日期如图4b、4c所示。

图4 ACI、BI和 NDSI指数在8个采样点突变日期的 日变化趋势Fig. 4 Diurnal variation trend of ACI, BI and NDSI at sampling points with mutation dates

结合突变性检验的结果可知,4号采样点ACI值在5:00出现最大值,也是突变发生的时刻,通过对原音频监听发现,由于天气变化,物种在当日鸣叫活动减少,突变性检验也有同样印证(Uk<0)。5号采样点在突变时刻受大雨天气影响,生物鸣叫活动受到雨声、风声的掩盖,开始显著减少。8号采样点发生突变当日天气晴朗,凌晨生物鸣叫活动增加。由图4b可知,BI指数多在夏季开始和结束时发生突变,3号采样点在突变日的鸣叫规律较不明显,结合音频分析,该采样点靠近交通道路,鸟类鸣叫容易受汽车鸣笛声影响。而4号采样点BI指数发生突变是因为受大雨天气影响,生物鸣叫行为减少。这种情况在5、6号采样点也有体现。由图4c可知,NDSI指数在各采样点突变日期均不同,1号采样点NDSI指数发生突变的原因主要是5月29日后,金丝猴研究基地内人工饲养活动增加,对声景产生影响。4号采样点在夏季开始和结束时多次发生突变,主要是由于恶劣天气使得生物声在声景中的占比减小,并且在突变后,NDSI指数值呈显著下降趋势,生物活动减少。5号采样点夜间蝉鸣声较明显,6号采样点晚间两栖动物鸣叫声较明显。

综上所述,引起突变的原因主要有两种:一是因为恶劣天气的影响,导致生物活跃程度减小,鸣叫行为减少;二是人为活动的影响,汽车鸣笛声、敲击声对声景造成较大的影响,这主要体现在1号采样点和3号采样点。此外,BI、NDSI指数受气候影响较大,恶劣气候极易影响指数的变化,ACI指数受气候影响较小,发生突变较少。

2.3 神农架声学指数夏季日动态变化特征

为了探究每个采样点在神农架大时空尺度内的日动态变化特征,根据一天所选时间点的数据对各采样点的3个指数分别进行平均,最后得出ACI、BI和NDSI指数在一天时间序列中的变化趋势,并绘制了8个采样点声学指数的日变化序列图,如图5所示。

图5 不同采样点的生态声学指数日动态变化特征Fig. 5 Diurnal dynamic changes of ecoacoustic indices at different sampling sites

从图5a可以看出,各采样点ACI指数变化趋势不一致。ACI指数最高值出现在6:30的4号采样点,最低值出现在0:10的5号采样点,大部分鸣叫动物夜间休息,早上逐渐苏醒并开始活动。2号采样点ACI值在夜间趋于平稳,黎明迅速上升,白天保持高值,傍晚下降,符合鸟类黎明/黄昏合唱的习性,说明在该采样点有高水平的鸟类活动。4号、9号采样点的日变化趋势有夜间指数值趋于平稳、黎明迅速上升、白天下降、傍晚出现暂时高值的特点,这是因为黎明大合唱消耗能量,物种会在合唱后进行觅食,从而减少了鸣叫行为[25]。3、4、9号采样点的ACI值变化幅度较大,鸣叫物种的数量较不稳定;而 1、6、8号采样点的ACI值在一天中变化幅度不大,声音的频率和振幅变化不大,物种数量及活跃程度在这些采样点较稳定。

从图5b可以看出,BI值在8个采样点的变化趋势相同,均符合夜间指数值较低、黎明迅速上升、白天保持高值、傍晚下降的特点,BI指数的日动态变化准确捕捉到了鸟类鸣唱活动。在日动态变化中,BI指数最高值出现在5:00的2号采样点,最低值出现在21:00的6号采样点。5号采样点夜间BI指数也保持较高的值,说明在该采样点,夜间物种鸣叫声音振幅高,频带间变化最小。

从图5c可以看出,NDSI值在2～9号采样点内日动态变化趋势相似,其中,2号采样点NDSI指数值较高,白天指数值接近1,该采样点海拔位置相对较高,人为活动较少,生物声在该声景中占主导地位。1号采样点的NDSI值变化没有明显的规律,该采样点在国家林业和草原局神农架金丝猴研究基地内,人为活动较频繁,NDSI指数对人类活动较敏感,因此变化幅度较大。5号采样点20:00至5:00间,NDSI值 <1,说明该采样点夜间声景受到人类干扰。

结合图5来看,在更大的区域内,一天中BI和NDSI指数值的变化趋势相同,5:00左右迅速上升,白天保持较高值,18:00左右开始下降,夜间趋于平稳。指数的变化符合鸟类黎明合唱和黄昏合唱的习性,说明在研究区域内声学指数的变化主要来自鸟类鸣叫。结合指数值的变化趋势可以推断,物种黎明合唱的时间为5:00左右,黄昏合唱的时间为18:00左右。5号采样点ACI值较低但BI值较高,这与该采样点螽斯类和蟋蟀类昆虫鸣声明显有关,重复且响亮的鸣声导致ACI值较低而BI值较高。

2.4 神农架声学指数随海拔梯度的变化特征

由于本研究采样点的海拔差异较大,为了对指数的空间差异进行研究,根据时间序列的日变化趋势,选取声学活动活跃程度较高的5:00和18:00的时刻,来探究声学指数与海拔之间的关系,并按照海拔梯度升序对各采样点进行排序分析,结果如图6所示。为了检测声学指数与海拔梯度之间的相关性,还计算了ACI、BI和NDSI指数与海拔的Spearman秩相关系数。

图6 神农架晨昏两个时刻(5:00和18:00)声学指数值随海拔梯度的变化Fig. 6 Acoustic index values as a function of altitude gradient at 5:00 and 18:00 in Shengnongjia National Park

从图6可以看出,ACI指数在5:00和18:00呈现相似的变化趋势,沿着海拔梯度呈现先上升后下降再上升,且在中海拔地区具有最高的声学指数值的变化规律,这与Colwell等[26]认为物种丰富度在海拔梯度分布呈现的驼峰形格局的规律一致。BI指数沿着海拔梯度均呈现先下降后上升再下降的变化规律,也是在中海拔地区具有最高值,但规律性不如ACI指数明显。NDSI指数随着海拔梯度的升高没有明显的变化规律。此外,从5:00和18:00的声学指数值看,5:00的声学指数值普遍高于18:00,说明黎明物种合唱时的生物多样性和活动性高于黄昏合唱,并且黎明合唱时是高值的采样点在黄昏合唱时同样也是高值。ACI、BI和NDSI指数在2号采样点,也就是2 375.36 m处出现最高值,说明在研究区域内海拔2 400 m左右的鸣叫物种丰富度较大。

计算可知,ACI指数与海拔相关性显著,相关性系数为0.762,而BI、NDSI和海拔相关系数分别为0.262、0.214,相关性不显著,因此,这3个指数中,仅ACI指数在本研究区域内可以较好地解释海拔梯度与声学丰富度之间的关系。

3 讨 论

生态声学生态位假说和声学适应假说解释了物种间的竞争会导致错位发声,动物发出的声音会根据所处的环境做出改变,因此,探究动物发声是评估物种多样性的有效工具之一,而声学指数是被证明与录音中生物声音多样性和丰富度密切相关。本研究基于生态声学的方法,对神农架动物鸣声连续时间序列数据进行了量化研究,分析了声学指数夏季时间序列的连续动态及日变化趋势,通过Mann-Kendall突变性检验的方法分析了指数变化特征。此外,还分析了声学指数随海拔梯度空间变化的特征。结果表明:①在所观测的夏季3个月时间内,ACI指数没有显著的变化规律,BI指数大致呈现双峰变化,NDSI指数呈现3峰/4峰变化;②Mann-Kendall突变性检验以及对突变日期的日变化分析发现,神农架国家公园内ACI指数在各采样点发生突变的次数较少,BI指数和NDSI指数突变发生次数较多,发生突变的主要原因是气候的变化和人为活动的干扰;③在大区域下,ACI指数没有明显的日变化趋势,BI指数和NDSI指数的日变化趋势符合动物黎明/黄昏合唱的习性;④ACI指数随海拔先上升后下降,在中海拔区域具有最大值,并且与海拔相关性显著,BI指数也在中海拔区域具有最大值,NDSI指数没有显著的变化趋势;物种黎明合唱时的生物多样性和活动性高于黄昏合唱。

RAS方法可超越物种水平,在景观水平上分析动物群落的声音特征,得到整体声景观多样性的定量描述,避免传统方法所必需的物种形态鉴定过程,既不需要鉴定物种,也不需要识别表型特征,而是估计群落所产生的复合声景的时间和频率的异质性,考虑了声源、环境、气候等因素影响声音多样性变化,能够快速对研究区域内的生物多样性做出评估。本研究的ACI指数连续动态变化没有表现出显著的变化规律,而BI指数呈双峰变化,NDSI指数呈3峰/4峰变化,结合Mann-Kendall突变性检验和突变日期日变化分析结果发现,对比3个指数,ACI指数对恶劣天气较不敏感,BI和NDSI指数较敏感。此外,BI和NDSI指数对人类干扰也较敏感。因此,ACI指数在少数采样点发生突变,而BI和NDSI指数几乎在每个采样点都发生了突变。结合突变日期的日变化分析来看,造成声学指数在夏季发生突变的原因主要有恶劣天气和人为干扰。游海林等[27]基于Mann-Kendall突变性检验的方法利用BI指数对鄱阳湖进行鸟类多样性监测发现,突变发生的主要原因是受人类活动导致,本研究指数发生突变的原因与此相似。此外,本研究连续性分析结果与Mann-Kendall突变性检验结果一致,证实了生物多样性特征及其活动习性会随着气候变化、人为活动增加等外部环境因素的变化而变化[4]。

在大区域下声学指数的日变化趋势不同于在小区域内所得出的结论[19],对比3个指数的日变化趋势来看,ACI指数不同于其他两个指数,没有表现出明显的日变化模式。分析认为,这可能是因为在大区域内,ACI指数与生物多样性的变化关系不强。在很多研究中,ACI指数都显示出与其他指数不同的日变化趋势,如Bradfer-Lawrence等[28]使用7个声学指数探索在巴拿马中部人类改良景观中的日变化和季节变化的结果一致,ACI指数是唯一没有表现出明显昼夜模式的指数。由于ACI指数可以有效地过滤掉一致的声音,比如重复的昆虫叫声[5],而这些声音是BI和NDSI指数产生昼夜差异的基础。ACI指数在生物多样性评估中表现也不一致,如Towsey等[29]发现ACI不能准确地代表(澳大利亚)布里斯班郊区的生物多样性,Doser 等[30]在研究密歇根州(美国)的生物多样性时得出了相同的结论。本研究结果与前人的研究结果一致。BI指数和NDSI指数准确反映了物种黎明/黄昏合唱,这是在温带森林中鸟类保卫繁殖区域、试图吸引伴侣或召集同伴的生态活动,说明在神农架国家公园内的生物声是由鸟类主导。

对于声学指数空间变化,有研究证明,样地跨越不同的景观或海拔梯度,对动物群落有强、弱或中等影响[23],如海拔会对一些区域生物多样性产生强烈的影响,海拔升高生物多样性呈线性下降[31],或者在中海拔具有最大的多样性[32],大多数声学指数支持动物声音复杂性随海拔升高而降低的观点,但是支持这一观点的研究主要是在热带和亚热带地区进行[33-34]。本研究是在温带森林中进行,结果显示ACI在中海拔区域具有最大值,并且与海拔相关性显著,BI指数也是在中海拔区域有最大值,但规律不如ACI指数明显,而NDSI指数没有明显的变化规律。因此,本研究中只有ACI指数能较好地支持声学指数沿海拔梯度有较显著空间变化特征的观点。由于缺乏神农架国家公园不同海拔的物种丰富度的数据,在今后研究中将会考虑物种丰富度与海拔的相关性,验证上述结论的科学性。

本研究结果说明BI、NDSI指数能较好地反映动物鸣声活动随时间变化的特征,而ACI指数能较好说明动物鸣声活动随海拔梯度空间变化的特征。因此,生态声学方法可以有效地帮助公园管理者对保护区内生物多样性进行快速地量化评估,为未来使用生态声学方法来监测保护区内的生物多样的变化提供技术参考。随着科学仪器硬件技术的发展,在未来研究中,将声学指数计算结合到前端声音传感器中进行实现,将大大提高生物多样性评估速度。