1. 中国科学院大学资源与环境学院, 北京 101408;
2. 中国科学院怀柔生态环境综合观测研究站, 北京 101408
2016年05月03日 收稿; 2016年05月19日 收修改稿
基金项目: 国家自然科学基金（31270512）资助
通信作者: 张娜, E-mail:zhangna@ucas.ac.cn
摘要: 探讨亚洲小车蝗的空间分布，以及在不同尺度上影响其分布的生境因子。基于在内蒙古锡林郭勒盟镶黄旗敖包音高勒调查和实测的亚洲小车蝗密度与生境条件数据，获得研究区蝗虫生境适宜性评分值（S值）的空间分布图。对整个研究区的S值和相关生境变量分别进行Morans' I空间自相关分析、融合Moran's I空间自相关尺度图的尺度方差分析、小波分析、或空隙度分析，以识别其等级结构和特征尺度。结果表明，亚洲小车蝗在 < 0.3 km、 < 0.7~1.0 km和 < 1.9 km这3个尺度域上出现斑块化分布。这是不同尺度上相关生境因子共同作用的结果。植被类型和植被盖度对蝗虫在这3个尺度域上的分布均有影响；坡度、坡向和地形湿度指数对蝗虫在中小尺度域上的分布有影响；高程对蝗虫在中大尺度域上的分布有影响；土壤沙粒含量、pH、有机质含量和植物叶片含氮量对蝗虫在中尺度域上的分布有影响。
关键词: 蝗虫密度生境适宜性特征尺度尺度域空间格局分析
Multi-scale spatial distributions of Oedaleus decorus asiaticus
LI Yao1, ZHANG Na1,2
1. College of Resources and Environment, University of Chinese Academy of Sciences, Beijing 101408, China;
2. Huairou Eco-Environmental Observatory, Chinese Academy of Sciences, Beijing 101408, China
Abstract: This study discusses the spatial distributions of Oedaleus decorus asiaticus (O.d. asiaticus) and the habitat factors which are responsible for the spatial distributions of O. d. asiaticus at different scales. Using the field survey data of O. d. asiaticus density and habitat conditions in Aobaoyingaole village, Xianghuangqi County of Xilingol League, Inner Mongolia, China, we obtained the spatial distribution of grasshopper habitat suitability value (S) within the study area. Then we carried out the spatial autocorrelation analysis of Moran's I, scale variance analysis coupled with Moran's I scalogram, wavelet analysis, or lacunarity analysis to detect the spatial patterns of S value and other related habitat variables. The results showed that O. d. asiaticus exhibited patchy distributions within the three scale domains: < 0.3 km, < 0.7~1.0 km, and < 1.9 km, which might result from the combined effects of many habitat factors at different scales. Vegetation type and vegetation coverage were responsible for the distributions of grasshopper individuals within all the three scale domains; slope, aspect, and topographic wetness index were responsible for the distributions within the moderate and small scale domains; elevation was responsible for the distributions within the large and moderate scale domains; and soil sand content, soil pH value, soil organic matter content, and leaf nitrogen content were responsible for the distribution within the moderate scale domain.
Key words: grasshopper densityhabitat suitabilitycharacteristic scalescale domainspatial pattern analysis
斑块化分布广泛存在于自然界的大多数动植物物种中[1-3]，其中，构成某个复合种群的亚种群斑块是最为典型的例子。蝗虫在草地中的分布不是随机的，也不是均匀的，可能在一些区域很少发生，而在另一些区域却频繁发生，也即，蝗虫的空间分布呈现斑块化[4]。我们在内蒙古锡林郭勒盟草原野外样地的调查中发现，亚洲小车蝗 (Oedaleus decorus asiaticus) 常常聚集在一个小斑块中，斑块内部可假设其发生随机分布，但斑块之间的蝗虫密度可能差异很大，从而呈现明显的斑块化分布[5-6]。另外，这种斑块化分布常常是多尺度的：在几千平方公里的区域尺度上，蝗虫可以分布在适宜其生存的不同景观中；在一个几十平方公里的景观中，可以分布在不同的适宜斑块中；而在几公顷的亚种群斑块内部，还有可能出现几百平方米的更小斑块。因此，有必要获悉自然界蝗虫发生的多尺度斑块化格局，这对揭示格局形成的原因十分关键。然而，目前尚未见到有关报道。
在不考虑人类干扰情况下，蝗虫的发生尺度是蝗虫个体感知相关生境因子 (如地形、土壤和植被) 并产生反应的尺度，而某时某刻所见的蝗虫的空间分布可能是不同尺度上相关生境因子共同作用的结果。国内外****已对蝗虫密度与生境因子之间的关系做了大量研究[7-9]，其中包括一些空间研究[8-9]。然而，目前的研究仍停留于单尺度的探讨上，尚未发现蝗虫在不同尺度上的发生及其与相关生境因子之间关系的报道。究其原因，与生态学中对尺度问题的重视程度不足有关，大多数实验生态学仍立足于几百平方米的单尺度样地。
通过观察蝗虫发生过程和相关生境变量在不同空间尺度上的动态，可以识别它们的等级结构和特征尺度。每个等级水平上过程或变量发生的尺度范围被称为尺度域，而特征尺度是两个相邻尺度域之间表示过渡或转折的某个 (或某些) 尺度。在某个等级水平上 (或尺度域内)，蝗虫发生过程及其影响因素的空间格局特征变化很小或随着尺度的变化单调变化；而在不同等级水平 (或尺度域) 之间，由于各种原因，它们的空间格局特征可能出现突然的转折[10]。
目前，识别等级结构和特征尺度的多尺度空间格局分析方法主要包括空间统计学方法、景观指数法和分维分析法。空间统计学方法主要包括空隙度分析、半方差分析、空间自相关分析、尺度方差分析和小波分析。这类方法的最大优势在于其本身是多尺度的，可根据相应指标随尺度变化的趋势和转折，检测过程或变量的空间格局特征，包括等级水平及其特征尺度、各尺度域大小、平均斑块大小、空间自相关性、某种格局出现的尺度及相对重要尺度等。除空隙度分析可用于类型数据外，其他方法多适用于数值型数据[10]。蔡博峰和于嵘[11]以中国三北防护林为例，对比半方差分析、尺度方差分析、小波分析和空隙度分析方法的特点和优劣，指出目前的尺度分析方法各有利弊。为消除用单一方法进行尺度识别造成的误差，本研究采用多种空间统计学方法识别蝗虫发生和相关生境变量的等级结构和特征尺度，以相互印证和补充。
本文着重探究以下问题：蝗虫的空间分布是否呈现多尺度斑块化格局？若是，那么蝗虫在不同尺度上的空间分布究竟与哪些生境因子有关？拟将多尺度思想引入蝗虫研究，并通过多尺度空间格局分析使蝗虫发生与生境因子关系的研究更加深入。
1 研究方法1.1 研究区概况研究区位于内蒙古锡林郭勒盟西南端镶黄旗的敖包音高勒 (42°12′11″~42°14′42″N，114°07′34″~114°11′02″E)(图 1)，总面积1.35×103 hm2，平均高程1 341 m。该区属温带大陆性气候，气候干燥，降水量少，蒸发强烈，且日照充足。年日照时数3 031.6 h，年均气温3.1 ℃，年降水量267.9 mm，降水主要集中在6、7和8月。研究区的主要植被类型是以克氏针茅 (Stipa krylovii)+糙隐子草 (Cleistogenes squarrosa) 为优势的温带丛生禾草草原。还有以大针茅 (S. grandis)+冷蒿 (Artemisia frigida) 为主的温带丛生禾草和矮半灌木草原，以及以大针茅+羊草 (Leymus chinensis)+猪毛菜 (Salsola collina) 为主的温带丛生禾草、杂类草草甸草原。地带性土壤为栗钙土，隐域性土壤包括风沙土、草甸土、盐碱土和石质土[5, 12]。
Fig. 1
 | Download: JPG larger image |
图 1 研究区的位置和野外调查样点 (2012年) Fig. 1 Location of the study area and field patches sampled in 2012 | |
研究区内主要发生的草地蝗虫有亚洲小车蝗 (O. d. asiaticus)、毛足棒角蝗 (Dasyhippus barbipes)、白边迦蝗 (Bryodema luctuosum) 和宽须蚁蝗 (Myrmeleotettix palpalis)。本研究将亚洲小车蝗作为研究对象。亚洲小车蝗为中期种，孵化期一般在6月上中旬，发生高峰期在7月上中旬。亚洲小车蝗为寡食性昆虫，主要以禾本科植物为主，在食料缺乏时，也取食莎草科 (Cyperaceae) 和鸢尾科 (Iridaceae) 植物[5, 13]。
1.2 野外调查和数据获取我们将整个研究区划分为203个斑块。每个斑块内部的生境条件近似。斑块的平均面积为6.63 hm2。在2012年7月上旬蝗虫开始迁移之前，用GPS测定每个斑块的大小、位置和空间布局，并绘制斑块图 (图 1c)。之后，调查每个斑块内部的蝗虫密度和生境状况，并假设这些信息在一个斑块内部的不同位置上均相同。
影响亚洲小车蝗空间分布的生境因子主要包括地形、植被和土壤等方面[8-9, 14]。地形因子包括高程、坡度和坡向。高程信息 (图 2(a)) 来自国际科学数据服务平台 (http://datamirror.csdb.cn/) 提供的数字地形图 (DEM)。通过ArcGIS 10.0表面分析模块得到坡向 (图 2(b)) 和坡度 (图 2(c)) 信息。另外，鉴于地表土壤水分与地形密切相关，本研究用DEM计算地形湿度指数 (topographic wetness index，TWI)(图 2(d))，并用TWI代替表层土壤含水量。
Fig. 2
 | Download: JPG larger image |
图 2 研究区生境因子的空间分布图 Fig. 2 Spatial distribution map of all the habitat factors in study area | |
土壤因子包括含砂量 (SSC)、有机质含量 (SOM) 和pH值。我们在研究区采集389个表层0~5 cm的土样，在实验室用比重计法测定含砂量，用容量法测定有机质，用电位法测定pH值。最后用克里金插值法得到整个研究区的信息 (图 2(h)，2(i) 和2(j))。
对植被类型进行调查，得到研究区的植被类型图，最终合并为3类，其优势种或次优势种分别为：禾本科植物、非禾本科+禾本科植物、非禾本科植物 (图 2(f))。根据植物组成的差异，我们选取131个样点，在每个样点中选取1个有代表性的1 m×1 m样方，调查样方内的所有植物种类及其盖度；剪取所有植物的叶片样品，在实验室用元素分析仪测得叶片的含氮量；用各植物种类的盖度加权平均得到该样点单位面积的叶片含氮量 (LN)，并用克里金插值法得到整个研究区的LN (图 2(e))。另外，用2012年7月中旬获取的Landsat影像反演得到研究区的植被盖度信息 (图 2(g))。
所有生境因子数据的空间分辨率均为30 m×30 m。尽管亚洲小车蝗密度与生境因子的关联程度会因某些生境因子 (如植被盖度) 的年际变化而变化，但是本研究关注的是二者的空间关系，而非时间动态。为此，仅取一年 (2012年) 的调查数据进行分析。
1.3 亚洲小车蝗生境适宜性评价假设蝗虫密度与生境适宜性之间有高度的对应关系，
可用蝗虫生境适宜性表示蝗虫潜在发生的可能性，此为亚洲小车蝗虫卵孵化的生境适宜性。在进行生境适宜性评价之前需对所选取的生境因子进行分类。先利用ARCGIS 10.0软件中的自然间断法对各生境因子进行初始分类，再结合野外调查结果对初始分类进行修正，从而使分类后的各个子类别能更好地体现蝗虫密度的变异。自然间断法是一种非监督分类方法，可使各类之间的差异最大化，而类内的差异最小。
对研究区调查得到的蝗虫密度 (GD) 也进行分类：GD < 3 ind·m-2(第1类，很低)，3≤GD < 5 ind·m-2(第2类，低)，5≤GD < 10 ind·m-2(第3类，较高)，GD≥10 ind·m-2(第4类，很高)。
使用Shen等[8]提出的方法进行草地蝗虫生境适宜性评价。该方法既考虑单因子对蝗虫发生的影响，也考虑各因子间两两交互作用的影响。但由于本研究中各因子间的两两交互作用不显著，故忽略交互作用，仅考虑各单因子的作用。
| $S(k) = \sum\limits_{i = 1}^N {{P_i}{Q_{ik}}}, $ | (1) |
| $S = \sum\limits_{k = 1}^4 {S(k){R_k}} .$ | (2) |
对所选取的167个野外调查样点的数据进行统计分析，确定各生境因子的不同类别在各密度等级上的隶属度 (表 1)。由表 1可见，研究区2012年蝗虫密度整体偏小，等级为“很高”的调查点很少。
Table 1
 表 1 各类别生境因子对不同蝗虫密度等级的隶属度Table 1 Membership of grasshopper density rank for each category of habitat factors
| 表 1 各类别生境因子对不同蝗虫密度等级的隶属度Table 1 Membership of grasshopper density rank for each category of habitat factors | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
利用Excel地理探测器 (Excel-GeogDetector) 对选取的生境因子进行因子探测分析，获得单个因子对蝗虫密度的影响。Excel-GeogDetector是王劲峰等[15]与Hu等[16]提出的一种空间分析方法。最初用于分析相关空间因子对某种疾病发生的影响小。将疾病的流行率或发生概率等作为健康风险指标H，计算H在影响因子D的不同类别上的总方差 (Vard) 与H在整个研究区的总体方差 (σH2)。通过比较Vard与σH2的相对大小确定D对H的贡献率。当σH2较大，且Vard相对较小时，意味着D的分类可以很好地解释H的变异，也表明D对H有较大影响。一般用PD (power determinant) 值表示这种影响力
| ${\rm{PD}} = 1-\frac{{{\rm{Va}}{{\rm{r}}_{d}}}}{{\sigma _H^2}} = 1-\frac{{\sum\limits_{i = 1}^L {({n_{D, i}}\sigma _{{H_{D, i}}}^2)} }}{{n\sigma _H^2}}, $ | (3) |
| $n = \sum\limits_{i = 1}^L {{n_{D, i}}} .$ | (4) |
基于生境变量的空间数据 (图 2)，用单个生境因子对蝗虫密度的PD值及生境因子各类别对不同蝗虫等级的隶属度，计算每个空间单元的S值，获得整个研究区S值的空间分布。
1.4 草地蝗虫生境适宜性和相关生境因子的多尺度空间格局分析采用Moran’s I空间自相关分析、尺度方差分析和小波方差分析识别蝗虫生境适宜性及各相关生境因子的等级结构和特征尺度。
1.4.1 Moran’s I空间自相关分析变量的自相关程度会随着样点间距的变化而变化。可用某个变量的Moran’s I系数值在不同样点间距上的空间自相关图识别该变量的空间自相关性及特征尺度。在空间自相关图中，一个变量的特征尺度在其Moran’s I系数由显著正值变到显著负值的间距转折处，即空间自相关曲线与零轴相交处[10]。
本研究使用ROOKCASE软件[17]计算草地蝗虫生境适宜性评分值 (S值) 和各生境变量在不同间距上的Moran’s I系数值。选取的最小间距为30 m，共选取100个间距，绘制Moran’s I空间自相关图。用系列Bonferroni校正法对所有间距上的Moran’s I系数值进行显著性检验[10]。
1.4.2 尺度方差分析尺度方差分析是Moellering和Tobler[18]提出的一种空间等级分析方法。最初被用来确定一个已知巢式等级系统各层次上相对变异的大小；现已被用来检测空间连续型数值变量的多尺度和等级结构特征，以确定那些在功能上值得研究的等级层次或空间尺度[10, 18]。以尺度方差与等级层次 (或空间尺度) 作图，可得到尺度方差图。通常情况下，尺度方差达到峰值时所对应的尺度水平是所研究格局斑块性表现突出的尺度，也指示不同等级层次之间的特征尺度。同时，通过比较尺度方差各个峰值的大小，也可以判别在某个或某些尺度水平上格局或过程的变化对整个系统变异性贡献的大小[10]。
鉴于尺度方差分析自身存在的缺陷，本研究采用Zhang和Zhang[19]提出的一种融合Moran’s I空间自相关尺度图的尺度方差分析方法。该法可解决尺度方差分析本身无法验证尺度方差峰值在统计上的显著性的问题，同时也可以验证划区方案的合理性。如果在尺度方差出现峰值的所有尺度水平上，Moran’s I尺度图均出现转折 (即空间自相关特征发生性质变化，如由显著正相关到显著负相关或由显著正相关到不相关的变化)，那么这些尺度水平上的尺度方差峰值在统计上显著，这些尺度水平对应着特征尺度，同时所使用的划区方案也是适宜的。如果不存在这种吻合，则表明所使用的划区方案不适宜，需选用其他划区方案。
在本研究中，一个空间单元的大小为30 m×30 m。对S值与各生境变量均使用4种划区方案 (表 2) 进行重采样。最小尺度水平为1×1个空间单元 (30 m×30 m)，最大为64×64个空间单元 (1 920 m×1 920 m)。计算每个划区方案不同尺度水平下的尺度方差 (或平方和SS)，并以尺度水平为横轴、尺度方差 (或SS) 的对数为纵轴绘制尺度方差图。
Table 2
表 2 划区和尺度变换方案 (引自Zhang和Zhang[19])Table 2 Four zoning and scaling systems (ZSS)
| 表 2 划区和尺度变换方案 (引自Zhang和Zhang[19])Table 2 Four zoning and scaling systems (ZSS) | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
1.4.3 小波方差分析小波方差是某个空间尺度上沿样带每个位置上的小波变换值的平方的平均值，是尺度的函数，与位置无关，因此，在识别多尺度格局上更加直观和方便[10, 20]。某个尺度上出现较高的小波方差值说明具有大量该尺度上的低强度格局或一些该尺度上的高强度格局[10]，由此可通过小波方差峰值出现的尺度检测格局主要发生的特征尺度，并指示这些尺度的相对重要性。
不同的小波函数对格局识别的结果不同。我们选用Morlet复小波，因其为连续小波，波形与原始数据近似，且识别特征尺度的能力强于常用的Haar小波和Mexican Hat小波[21]。
在研究区内沿着南北方向选取一条具有代表性的样带。使用MATLAB 2014 a软件小波分析模块中的Morlet复小波函数对该样带上的S值和各生境变量进行小波变换，并计算不同尺度的小波方差，绘制小波方差图。为消除边缘效应对小波分析结果的影响，最终参与小波方差计算的样带长度小于初始进行Morlet复小波变换的样带长度。
1.4.4 空隙度分析空隙度可度量一个几何体的空隙偏离均匀分布的程度。可通过空隙度随滑动框大小 (尺度) 变化的对数曲线图来识别景观中某个类型变量的多尺度格局。一般来说，该类型的特征尺度位于空隙度对数曲线的转折处。这种转折可能表现为由缓慢的下降变为迅速下降，也可能表现为由下降的直线转变为曲线或由下降的曲线转变为直线[10]。
本研究使用APACK 2.23软件[22]计算植被类型和坡向这两个类型变量当滑动框取不同大小时的空隙度，并基于空隙度-滑动框大小的对数曲线图进行多尺度格局分析，识别特征尺度。滑动框大小在1~240个空间单元之间变化，每个单元的大小为30 m×30 m。
2 结果与分析2.1 草地蝗虫生境适宜性评价结果及其验证计算167个样点对应的S值。并将S值分为3类，分别对应的蝗虫密度范围为0~5、5~10和10 ind·m-2以上。
对预留的未参与隶属度和PD值计算的35个验证点 (图 1)，计算相应的S值，并划入不同的等级。结果发现，25个样点 (占总验证点的71.4%) 的生境适宜性等级与实测蝗虫密度完全吻合。
我们使用蝗虫四龄之前的调查数据进行生境适宜性评价。这时蝗虫尚未发生迁移，某地蝗虫密度与该地上年的产卵情况密切相关。如果上年蝗虫未在该地产卵，则即使该地的生境非常适宜蝗虫生存，也几无蝗虫，这种生境适宜性与蝗虫密度的不匹配会降低模拟精度。为消除这种可能影响，我们于2012年7月下旬到8月初蝗虫达到四龄之后，对研究区的蝗虫密度和生境条件进行了第二次调查。四龄之后蝗虫可能发生迁移，某地的适宜生境有可能吸引其他地方的蝗虫迁入，造成该地蝗虫密度上升。因此，可用第二次调查的蝗虫密度数据对第一次调查的数据进行校正。分析后发现，在生境适宜性等级与第一次实测蝗虫密度不吻合的10个样点中，有2个点的生境适宜性等级与第二次实测蝗虫密度吻合。由此，校正后，共有27个样点 (占总验证点的77.1%) 的生境适宜性等级与实测密度完全吻合。以上验证结果说明，S值与蝗虫密度有很好的对应关系，本研究所用的生境适宜性评价方法是合理可靠的。
在此基础上，我们获得了整个研究区的S值 (图 3(a)) 及其分类结果 (图 3(b))。
Fig. 3
 | Download: JPG larger image |
图 3 草地蝗虫生境适宜性 (a) 和分类等级 (b) Fig. 3 Habitat suitability for grasshoppers and classifications | |
2.2 草地蝗虫生境适宜性的多尺度空间格局草地蝗虫生境适宜性评分值 (S值) 的空间自相关分析结果表明，Moran’s I系数值随着间距的增大而单调降低，经历由正到零到负的变化 (图 4(a))。说明S值在整个景观中呈梯度变化，在较小间距上呈现显著正相关，在较大间距上呈现显著负相关；在由正到负的转折点 (1.98 km) 附近，S值不具有空间自相关性，而是具有随机分布特征。该转折点对应着一个特征尺度，指示此处S值的空间格局发生性质上的变化。
Fig. 4
 | Download: JPG larger image |
| (a) 中，实心三角形表示用系列Bonferroni法校正后呈现显著空间自相关，空心三角形表示呈现空间不相关。(c)~(f) 中，空心菱形表示尺度方差在该处出现峰值；空心正方形表示该处对应的Moran’s I系数值发生显著变化，空间格局由显著正相关变为不相关。 图 4 草地蝗虫生境适宜性评分值的多尺度空间格局分析 Fig. 4 Muti-scale spatial pattern analysis for habitat suitability value of grasshoppers | |
S值的小波方差分析结果表明，小波方差在尺度为1.89 km处达到最大峰值，在尺度为1.04和0.72 km处出现第2和第3个峰值，在尺度为0.33 km处出现第3个峰值 (图 4(b))。因此，从小波方差分析可识别出S值的3个特征尺度：1.89、0.72~1.04和0.33 km。
尺度方差分析中，当使用第1种划区方案 (表 2) 时，S值的尺度方差在16×32尺度水平 (0.48 km×0.96 km) 和32×64尺度水平 (0.96 km×1.92 km) 处出现峰值。但Moran’s I系数值在这2个尺度水平上并未出现显著变化，说明空间自相关特征并未发生性质变化 (图 4(c))。因此，不能确定这2个尺度水平为S值的特征尺度。当使用第3划区方案 (表 2) 时，S值的尺度方差在4×8尺度水平 (0.12 km×0.24 km) 和64×64尺度水平 (1.92 km×1.92 km) 处出现峰值。但Moran’s I系数值在0.12 km×0.24 km的尺度水平上未出现显著变化，说明空间自相关特征并未发生性质变化 (图 4(e))。因此，不能确定这2个尺度水平为S值的特征尺度。当使用第2和第4种划区方案 (表 2) 时，S值的尺度方差在64×64尺度水平 (1.92 km×1.92 km) 处出现峰值，且此处的Moran’s I系数值也出现显著变化 (图 4(d)和4(f))，因此，1.92 km为S值的特征尺度。
综上所述，亚洲小车蝗生境适宜性空间分布的特征尺度约为：0.3、0.7~1.0和1.9 km。
2.3 生境变量的多尺度空间格局2.3.1 数值型变量的多尺度空间格局因篇幅所限，仅以植被盖度为例来说明对生境变量多尺度格局的识别结果。植被盖度的空间自相关分析结果表明，Moran’s I随着间距的增大而单调递减 (图 5(a))，并在间距为1.92 km处与零轴相交。说明植被盖度在整个景观中呈梯度变化，其由正相关向负相关转折的一个特征尺度为1.92 km。植被盖度的小波方差分析结果表明，小波方差在尺度为1.83和0.48 km处分别出现大小峰值 (图 5(b))，即植被盖度的特征尺度为0.48和1.83 km。在尺度方差分析中，当使用第4种划区方案时，植被盖度的尺度方差在64×32尺度水平 (1.92 km×0.96 km) 处出现峰值 (图 5(c))，且此处的Moran’s I系数值也发生明显变化，空间格局由显著正相关变为不相关。因此，植被盖度的一个特征尺度约为1.36 km。综上，植被盖度空间分布的特征尺度约为：0.48、1.36和1.83~1.92 km。
Fig. 5
 | Download: JPG larger image |
图 5 植被盖度的多尺度空间格局分析 Fig. 5 Muti-scale spatial pattern analysis for vegetation coverage | |
对其他7个数值型生境变量也做同样的分析 (表 3)。
Table 3
表 3 生境变量的特征尺度Table 3 Characteristic scale of habitat variables
| 表 3 生境变量的特征尺度Table 3 Characteristic scale of habitat variables |
2.3.2 类型变量的多尺度空间格局对禾本科植物，当滑动框小于210 m时，空隙度-滑动框对数曲线下降非常缓慢；当滑动框大于210 m时，对数曲线近直线下降；当滑动框大于1.98 km时，呈负对数曲线变化，直至达到零值 (图 6(a))。说明禾本科植物在小尺度上呈明显的聚集分布，在中尺度上具有分形结构，而在大尺度上呈随机分布。曲线出现转折的滑动框大小对应禾本科植物空间分布的特征尺度。对非禾本科+禾本科植物和非禾本科植物，空隙度-滑动框对数曲线的变化趋势与禾本科植物近似，只是曲线出现转折的滑动框大小不同 (图 6(a))，对应的空间分布的特征尺度为0.21和1.2 km。
Fig. 6
 | Download: JPG larger image |
图 6 植被类型和坡向的空隙度与滑动框大小的对数图 Fig. 6 Log plots of lacunarity vs. gliding box size for vegetation type and aspect | |
对半阴坡和半阳坡坡面，当滑动框小于420 m左右时，空隙度-滑动框对数曲线下降较为缓慢；当滑动框大于420 m时，对数曲线近直线下降；当滑动框大于1.20~1.59 km时，呈负对数曲线变化，直至达到零值 (图 6(b))。指示半阴坡和半阳坡坡面在小尺度上呈一定的聚集分布，在中尺度上具有分形结构，而在大尺度上呈随机分布。与半阴坡和半阳坡坡面不同，对阳坡和阴坡坡面，对数曲线只在滑动框为1.20 km时出现转折 (图 6(b))。指示阳坡和阴坡坡面在中小尺度上具有分形结构，而在大尺度上呈随机分布。曲线出现转折的滑动框大小对应坡向空间分布的特征尺度。
3 讨论与结论对生境适宜性进行多尺度空间格局分析的结果表明，亚洲小车蝗在 < 0.3 km、 < 0.7~1.0 km和 < 1.9 km这3个尺度域上出现斑块化分布。在不同尺度上发生变化的生境因子可能对此有不同的贡献。坡度、坡向、TWI、植被类型、植被盖度在小尺度 (小于0.5 km范围内) 上呈现很强的聚集性 (表 3)，与蝗虫在小尺度上的分布相对应，暗示这些生境因子对蝗虫小尺度聚集的可能影响。高程、坡度、坡向、TWI、植被类型、植被盖度、单位面积的LN、SSC、SOM和土壤pH在中尺度 (0.6~1.5 km范围内) 上呈现很强的聚集性 (表 3)，与蝗虫在中尺度上的分布相对应，暗示这些生境因子对蝗虫中尺度分布的可能影响。高程、植被类型和植被盖度在较大尺度 (1.8~2.2 km范围内) 上呈现很强的聚集性 (表 3)，与蝗虫在大尺度上的分布相对应，暗示这些生境因子对蝗虫大尺度聚集的可能影响。
研究区禾本科植物的斑块面积总体较大 (图 2(f))，平均可达1.08 km2，故可在大尺度上影响蝗虫分布。而非禾本科植物+禾本科植物和非禾本科植物的平均斑块面积较小，镶嵌于禾本科植被类型中，仅能在中小尺度上影响蝗虫分布。研究区植被盖度跨度很大 (9.20%~84.30%)，具有较高的空间异质性，且这种异质性在小、中和大尺度上均有明显表现 (图 2(g))。植被盖度的这种多尺度变化会对蝗虫分布产生影响。本研究LN是用样点上不同植物种类的盖度 (相当于叶生物量的比重) 和叶片含氮量计算的单位面积上的值，因此，这里的LN既与植被类型有关，也与植被盖度有关，其空间分布也受到它们的共同影响。面积占绝对优势的禾本科植物 (克氏针茅、大针茅、糙隐子草、羊草) 本身的LN差异很小，平均为2.0%~2.7%，决定了LN在小尺度上的变化不明显。另外，无法识别小尺度上的变化，也可能与取样数量 (131个) 较少有关。但因植被盖度在中小尺度上变化明显，LN在中尺度 (约1.32 km) 上仍表现出较大的变化。
研究区坡度的异质性较强 (图 2(c))。即使是相邻不同方向的坡面也可能会有明显差异，使其在300~500 m尺度上呈现较大的变化。最西部和东南部地势较高，有坡度较大的坡面，而其他地方地势起伏较为平缓，地貌上的这种变化决定了坡度在中尺度 (1.0~1.5 km) 上的变化。研究区有很多小丘陵，致使坡向在局域范围内变化很大 (图 2(b))，因为在自然界有面南坡就必定有与之相应的其他坡向的相邻坡面。同时，研究区也有一些大丘陵，在较大范围内坡向趋于一致 (图 2(b))。坡度和坡向的这种空间变化可影响蝗虫在中小尺度上的变化。TWI是单位等高线上的汇流面积与坡度之比，与坡度类似，其分布可影响蝗虫在中小尺度上的变化。
研究区土壤以栗钙土为主，土壤理化性质的量值总体变化范围不大，空间分布也主要在中尺度上发生变异，而在局域范围内的变化很小。SOM在西部整体较高，而在东部整体较低 (图 2(h))；与之近乎相反，SSC和pH值在西部和西南部较低，而在北部、东北部和东南部较高 (图 2(j))。在这些不同方位，SOM、SSC和pH均在约1.5 km的范围内有很强的聚集性，形成中尺度上的热点区。土壤理化性质的这种中尺度变化会对蝗虫分布产生影响。
研究区高程大多为1 300~1 365 m (占86%)，空间异质性较低，在北部整体较低，差异较小；而在南部整体较高，差异也较小 (图 2(a))。也即，高程仅在中大尺度 (大于1.36 km) 上，才会发生较大变化，进而影响蝗虫分布。在小尺度 (小于1.44 km) 上高程的变化很小，可忽略对蝗虫分布的影响。之前认为高程是影响亚洲小车蝗空间分布的最重要生境因子的研究结果[5, 23]是基于大中尺度数据得来的。
由此可见，研究区亚洲小车蝗的空间分布是不同尺度上相关生境因子共同作用的结果。
本研究没有考虑影响蝗虫的人为因素，如飞机灭蝗和放牧等。一般来说，在一段时间内，放牧是局域性活动，其对蝗虫分布的影响集中在中小尺度上；而飞机灭蝗可在大尺度上影响蝗虫分布。在今后的研究中，应该考虑这些人为因素的影响。另外，蝗虫发生和生境变量可能还存在一个更小的尺度域，但受所用方法及最小空间单元的限制，尚无法获取。本研究没有考虑时间因素，鉴于蝗虫的发生会有年际变化，在今后的研究中还需对不同年份蝗虫的空间分布进行探讨。
在本研究使用的方法中，空间自相关分析法的格局识别率最高，可以保证至少识别出一个特征尺度。用小波方差分析会出现无法识别的现象，可能与样带选取有关；但该法对多尺度格局的检测最有效，这与以往的研究结果一致[11, 21]。在进行小波分析时，通过小波转换将原数据分解成具有不同频率 (空间尺度) 的多个组分，从而可将空间结构定量为尺度和样带位置的函数，由此可反映原数据的多尺度格局特征。尺度方差分析通常只能识别出较大的特征尺度。空隙度分析法可以对类型数据进行特征尺度的识别，但当空隙度-滑动框大小的对数曲线转折不明显时，识别特征尺度并不容易[11]。鉴于这些方法各有优劣，在进行特征尺度识别时，最好同时使用，以相互印证，提高识别的准确性。然而，使用哪些方法并没有统一的规则，在具体研究中，综合运用哪几种方法能够达到最好的识别效果，并且省时省力，是一个非常值得探讨的问题，且相关方法尚需不断地发展和完善。
?感谢中国农业科学院植物保护研究所张泽华研究员和吴惠惠，以及内蒙古锡林郭勒盟镶黄旗草原工作站工作人员对本研究给予的大力支持。
参考文献
| [1] | Nee S, May R M. Dynamics of metapopulations: habitat destruction and competitive coexistence[J].Journal of Animal Ecology, 1992, 61(1):37–40.DOI:10.2307/5506 |
| [2] | 刘任涛, 毕润成, 闫桂琴. 复合种群生态学研究现状与展望[J].山西师范大学学报 (自然科学版), 2006, 20(3):56–60. |
| [3] | Tilman D, May R M, Lehman C L, et al. Habitat destruction and the extinction debt[J].Nature, 1994, 371(6492):65–66.DOI:10.1038/371065a0 |
| [4] | 赵成章, 李丽丽, 王大为, 等. 黑河上游天然草地蝗虫空间异质性与分布格局[J].生态学报, 2012, 32(13):4166–4172. |
| [5] | 张红艳, 张娜, 陈晓燕. 内蒙古镶黄旗草地蝗虫潜在发生可能性评价[J].应用生态学报, 2012, 23(1):223–234. |
| [6] | Zhang N, Jing Y C, Liu C Y, et al. A cellular automaton model for grasshopper population dynamics in Inner Mongolia steppe habitats[J].Ecological Modelling, 2016, 329:5–17.DOI:10.1016/j.ecolmodel.2016.03.002 |
| [7] | 张军霞, 赵成章, 殷翠琴, 等. 黑河上游天然草地亚洲小车蝗蝗蝻与成虫多度分布与地形关系的GAM分析[J].昆虫学报, 2012, 55(12):1368–1375. |
| [8] | Shen J, Zhang N, Gexigeduren, et al. Construction of a GeogDetector-based model system to indicate the potential occurrence of grasshoppers in Inner Mongolia steppe habitats[J].Bulletin of Entomological Research, 2015, 105:335–346.DOI:10.1017/S0007485315000152 |
| [9] | Zhang N, Zhang H Y, He B, et al. Spatiotemporal heterogeneity of the potential occurrence of Oedaleus decorus asiaticus in Inner Mongolia steppe habitats[J].Journal of Arid Environments, 2015, 116:33–43.DOI:10.1016/j.jaridenv.2015.01.019 |
| [10] | 张娜. 景观生态学[M].北京: 科学出版社, 2014. |
| [11] | 蔡博峰, 于嵘. 景观生态学中的尺度分析方法[J].生态学报, 2008, 28(5):2279–2287. |
| [12] | 杨婷婷, 吴新宏, 姚国征, 等. 基于TM影像的镶黄旗草原沙化治理工程生态效益研究:以2004年封沙育林工程为例[J].水土保持研究, 2009, 16(1):204–207. |
| [13] | 陈永林. 中国主要蝗虫及蝗灾的生态学治理[M].北京: 科学出版社, 2007. |
| [14] | 查高德, 赵成章, 张军霞, 等. 黑河上游天然草地亚洲小车蝗雌雄虫的空间格局[J].生态学杂志, 2013, 32(11):3022–3028. |
| [15] | 王劲峰, 廖一兰, 刘鑫. 空间数据分析教程[M].北京: 科学出版社, 2010. |
| [16] | Hu Y, Wang J F, Li X.H, et al. Geographical detector-based risk assessment of the under-five mortality in the 2008 Wenchuan earthquake, China[J].PLoS ONE, 2011, 6(6):e21427.DOI:10.1371/journal.pone.0021427 |
| [17] | Sawada M. ROOKCASE: an excel 97/2000 visual basic (VB) add-in for exploring global and local spatial autocorrelation[J].Bulletin of the Ecological Society of America, 1999, 80(4):231–234.DOI:10.1890/0012-9623(1999)080[0231:TT]2.0.CO;2 |
| [18] | Moellering H, Tobler W. Geographical variances[J].Geographical Analysis, 1972, 4(1):34–50. |
| [19] | Zhang N, Zhang H. Scale variance analysis coupled with Moran's I scalogram to identify hierarchy and characteristic scale[J].International Journal of Geographical Information Science, 2011, 25(9):1525–1543.DOI:10.1080/13658816.2010.532134 |
| [20] | Bradshaw G A, Spies T A. Characterizing canopy gap structure in forests using wavelet analysis[J].Journal of Ecology, 1992:205–215. |
| [21] | 郁文, 刘茂松, 徐驰, 等. 南京市城市景观的特征尺度[J].生态学报, 2007, 27(4):1480–1488. |
| [22] | University of Wisconsin-Madison. APACK. Users guide and application: USA[EB/OL].(2012-04-20)[2016-04-20].http://forestandwildlifeecology.wisc.edu/staticsites/mladenofflab/Projects/Apack/index.htm. |
| [23] | 倪绍祥, 蒋建军, 巩爱歧, 等. 环青海湖地区草地生境的蝗虫潜在发生可能性评价[J].生态学报, 2002, 22(3):285–290. |
