我国农林业受到有害生物危害造成的生态和经济损失巨大, 掌握其空间分布格局及动态变化规律是进行防治的前提。随着地统计学(Geostatistics, GS)和地理信息系统(Geographic Information System, GIS)的发展和应用, 有害生物种群的时空动态变化研究有了新的思路和工具^[1]。地统计学与传统统计学不同, 不仅关注样本的数字特征以外, 还关注样本数据的空间分布特征, 侧重于区域化变量的空间结构分析、模拟、空间插值^[2]。地理信息系统具有强大的空间分析功能, 在查询、编辑、显示、以及海量数据处理方面具有优势, 可用于空间数据的管理、模拟与预测, 以及可视化展示等^[3—4]。将二者相结合应用于害虫种群的时空动态研究十分广泛, 比如, 季荣^[5]、娄国强^[3]、潘宏阳^[6]等人将GS和GIS分别应用于东亚飞蝗卵块、春尺蠖及松材线虫的时空分布研究中。

美国白蛾Hyphantria cunea(Drury)是世界性的检疫害虫, 不仅寄主广泛, 食性杂, 会危害林木、花卉、农作物、蔬菜等, 而且繁殖量大, 在我国有的地方一年2代, 有的地方甚至3代, 因此经常暴发成灾, 种群数量大时甚至会入户扰民, 影响人们的正常生活^[7]。自2003年在北京市平谷区首次监测到美国白蛾幼虫后, 美国白蛾在北京地区逐渐扩散, 对社会环境、经济、生态等各方面造成了极大的损失^[8]。本文利用GS和GIS, 对北京地区美国白蛾种群的时空分布特点及动态变化规律进行了分析, 明确了北京地区美国白蛾种群的分布情况, 并结合气候数据、地形数据、人口统计数据探究影响美国白蛾种群时空分布变化的因素。

1 材料与方法 1.1 研究区域概况

本研究地点为北京市, 位于东经115.7°—117.4°, 北纬39.4°—41.6°, 东部和南部是平原, 而西部和北部是山区。平原地区的年平均气温在12℃左右, 7月份温度最高, 平均气温26℃；1月份温度最低, 平均气温在-5℃左右。而海拔高的山区温度普遍比平原低, 年平均气温10℃左右, 最热月份和最冷月份均比平原地区低2—5℃, 属于暖温带半湿润半干旱季风气候。美国白蛾在北京地区一年有3代^[9], 危害极其严重。

1.2 数据来源

本研究中使用的美国白蛾监测数据从北京市林业保护站测报处获得, 主要利用性诱捕器(Trap)和诱虫灯(moth-killing lamp)在北京市16个区县的3000多个采样点诱捕美国白蛾。研究选取自2008—2016年的监测数据, 原始监测记录共38371条, 经过整合之后, 用于插值与分析的样本点共有14372个。气候数据来源于国家气象科学数据共享服务平台(http://data.cma.cn), 包括各月份的最高最低温度、平均温度、日照时数、平均相对湿度, 以及降水量。地形数据主要使用的是地理空间数据云网站(http://www.gscloud.cn)的数字高程模型数据。北京市常住人口数据来源于2014年北京统计年鉴。

1.3 研究方法 1.3.1 数据预处理

利用Geostatistical Analyst工具条中探索数据的Voronoi图^[10]功能, 对整理后美国白蛾的监测数据进行预处理, 去除观测中的局部异常值, 提高插值精度。异常值主要是由于在诱捕过程中, 诱捕器悬挂位置可能存在过密的情况, 导致相邻诱捕器的虫口差异较大；其次在人工统计过程中可能出现个别差错。处理异常值的方法是将Voronoi统计量选择“聚类”, 聚类结果通过不同的色带显示, 根据地理学第一定理可知, 空间上越相近的事务, 相关性越大^[11—12], 若局部值异常高或异常低, 则Voronoi图中该异常点所在的多边形颜色与周围不同, 该点则需要删除。

1.3.2 半变异函数分析

半变异函数分析是为了揭示区域化变量的空间依赖关系, 从而推测一定范围内该区域化变量的变化规律^[13—14]。半变异函数可以检验所采集的样本数据中是否存在空间自相关^[15], 若空间自相关弱或没有空间自相关则不能用克里金进行插值。其计算公式^[16]如下：

式中, γ (h) 是样点间距的半方差, 该值随着h的增加而增加, Z(x_i)和Z(x_i+h)表示在x_i和x_i+h位置上对美国白蛾诱捕的虫口数量, N(h)是间距为h的成对样本点总数^[17—18]。

半变异函数曲线如图 1所示：块金常数(C₀)能够反映半变异函数的空间异质性^[19—20]；基台值(Sill)是随着样本点的间距不断增大而趋于稳定的常数^[21]；变程(a)表示在一定观测尺度下, 空间相关性的作用范围, 在变程范围内, 样本点距离越近, 其相似性越大, 当距离大于变程, 样本则可以看作是完全独立的^[22—23]。常见的半变异函数模型有：指数模型、高斯模型、球状模型、线性有基台值模型、线性无基台值模型、对数模型等^[24—25]。通常, 残差平方和(RSS)越小, 决定系数(R²)越大, 模型精度越高^[26], 最佳模型的平均值标准化误差和平均误差接近于0。

1.3.3 空间模拟插值

克里金(Kriging)插值法是基于半变异函数理论分析的结果^[27], 无偏的估算一定范围内未知点取值的一种方法^[28]。该方法利用周围样本点的值确定待估计点的权重, 以估计未测量点的值, 充分考虑了样本点的大小和样本点与待估计点之间的空间自相关关系^[29]。克里金插值的结果准确度与样本量有关, 样本量越多, 结果越精确。本研究利用普通克里金插值法对北京地区美国白蛾种群的空间分布特征进行定量分析, 制作2008—2016年不同世代美国白蛾的空间分布地图, 直观地展示其分布情况, 进而为防治提供理论依据。

1.4 统计及分析软件

昆虫的发育受气候因素影响是综合性的, 且在不同年份之间表现出的差异较大。为了明确北京地区美国白蛾的成虫数量与气象因子之间的关系, 本研究使用R软件进行不同世代成虫数量与气象因子的Pearson相关性分析；然后利用Vegan包^[30]进行美国白蛾的成虫数量去趋势对应分析(Detrended Correspondence Analysis, DCA), 因为DCA最大梯度长度值小于3, 所以最终选用限制性排序法——冗余分析(Redundancy Analysis, RDA)^[31—33]。最后基于999个随机排列进行蒙特卡罗置换检验以确定每个气候因子在解释成虫数量变化中的相对重要性和显著性^[34—35]。

本研究使用SPSS 18.0对数据进行批量处理, 使用百度地图API获取美国白蛾诱捕点的坐标, 半变异函数分析与半变异理论模型拟合是在GS+ 7.0软件中完成, 数据预处理、地统计分析及结果制图在ArcGIS 10.2软件中完成。

2 结果与分析 2.1 北京地区美国白蛾种群的半变异函数分析

根据北京市2008—2016年各代美国白蛾成虫数量的半变异函数分析和理论模型拟合结果(表 1)可以看出, 各代美国白蛾成虫的种群分布均存在空间相关性。模型拟合程度较高, 块金值相较于基台值都比较小, 平均占比为35%, 最低为3.6%, 说明随机因素引起的空间变异小于空间相关性的影响。一般情况下, 球型和指数型半变异函数说明所研究的空间分布对象属聚集分布^[36], 表 1中理论模型结果显示, 球形模型和指数模型共有18个, 占比66.67%, 说明北京地区美国白蛾种群在多数情况下呈聚集分布。其他的理论模型有高斯模型、线性有基台值模型, 表明种群在一定程度上呈现中等程度聚集分布。

2.2 北京地区美国白蛾种群的空间模拟插值

利用普通克里金插值法, 结合变异函数结构拟合结果, 对2008—2016年北京地区美国白蛾种群的空间分布进行插值分析, 结果显示(表 2), 各世代的平均误差和平均标准化误差基本接近0, 说明测量值和预测值之间误差较低；均方根误差与平均值标准误差之间的差值用来衡量模型预测测量值的接近程度, 不同年份各世代的差值不同, 但都基本达到了要求的精度, 北京地区美国白蛾种群分布插值地图(图 2)与实际分布情况拟合度较好。

2.3 北京地区美国白蛾种群动态分析

克里金插值结果显示, 2008—2016年北京地区单年内美国白蛾种群动态的变化情况基本类似。各世代成虫的分布有多个聚集中心, 基本上都位于北京中心城区：东城区、西城区、朝阳区、海淀区、丰台区、顺义区等区县。其次, 美国白蛾种群的聚集中心数量越多的年份, 最大虫口密度越大, 尤其是2010年和2013年, 并且聚集中心的虫口密度与其它区域的虫口密度差值也较大, 反之, 聚集中心的数量减少, 与周边的种群密度差值变小, 说明美国白蛾的聚集程度与种群数量密切相关, 基本上呈正相关关系。

另外, 统计各世代种群数量发现, 北京地区美国白蛾的种群数量呈现有规律的波动变化(图 3), 单年内种群数量表现出下降的趋势, 第一代成虫与第二代成虫的种群接近, 但都低于当年越冬代成虫的数量, 到第二年越冬代成虫数量回升, 依此呈现周期性变化。由此可知, 今后的防治重点应在越冬代蛹羽化成虫之前就开始防治, 以控制美国白蛾越冬代成虫数量。

2.4 北京地区美国白蛾时空分布影响因素分析 2.4.1 气候因素

蒙特卡洛置换检验结果(表 3)表明, 气象因子对北京地区美国白蛾成虫数量的影响是显著的(P=0.017 < 0.05)。从RDA排序图(图 4)中可以看出, 前两个轴共解释了94.28%的数据, Tavg7、Tmax8S、RHU3与RDA1呈正相关关系, SSD2、Tmax4则与RDA1呈负相关关系。美国白蛾以滞育蛹越冬, 但温度升高会加快昆虫滞育期间的新陈代谢率, 所以2月份长时间的日照增加了美国白蛾越冬期间的能量消耗, 进而降低其越冬代的成活率^[37—38]；而美国白蛾解除滞育需经过一段较长时间的低温才能完成^[39—40], 并且喜欢潮湿的环境, 北京地区美国白蛾越冬代在4月初逐渐开始羽化成虫^[41], 4月份气温越高, 则相对湿度越低, 因此不利于解除滞育, 甚至会导致蛹干瘪死亡。7月份是第一代成虫羽化高峰期, 成虫羽化2d后即可产卵^[41], 由于北京地区夏季雨热同期, 所以高温潮湿的环境十分适合美国白蛾生长发育。

8月上旬为第二代幼虫野外网幕高峰期, 下旬越冬代幼虫开始出现^[41], 次年的越冬代种群数量与前一年第二代成虫的数量密切相关(r=0.785), 8月份最高温度与当年第二代成虫数量正相关, 进而影响次年美国白蛾种群数量。虽然Tmin9S的检验结果并不显著(P=0.118>0.05), 但其与RDA2的相关性极强(r=0.956)。对美国白蛾而言, 滞育诱导的敏感虫态是在幼虫期1—3龄, 诱导其越冬滞育的主要诱因之一就是低温^[42], 9月份的温度越低, 则幼虫会更敏感, 更早进行能量贮备, 成功越冬的机率更大, 导致次年越冬代成虫数量上升。

2.4.2 地形

从图 5可以看出, 北京地区美国白蛾的诱捕记录主要分布在城区中部和东南部的低海拔平原地区, 而在北部和西部的高海拔地区诱捕点比较少, 可以看出美国白蛾种群的空间分布受海拔高度的影响十分显著。有研究表明, 随着海拔的升高, 美国白蛾生存概率下降^[43], 并且, 美国白蛾的飞行能力弱, 成虫飞行距离一般为100m左右^[44]。因此, 在北京地区高海拔的山区美国白蛾种群数量很小。

2.4.3 人为干扰

北京地区美国白蛾的诱捕记录分布与常住人口密度高度一致。从图 6可以看出, 北京中心城区的常住人口显著高于周围郊区, 诱捕到美国白蛾的分布点则刚好集中在东城区、西城区、海淀区、朝阳区、丰台区, 而人口较少的郊区只零星分布着少量美国白蛾的诱捕记录, 如门头沟区、延庆区、怀柔区等区县。结合各年的插值结果地图可以看出, 人口密集的区域往往是美国白蛾种群密度最大的聚集中心, 虫口密度高出远郊区县几十倍, 说明人为干扰因素也会显著影响美国白蛾种群的分布。

3 讨论

相对于最早入侵的平谷区, 2008—2016年期间, 美国白蛾的聚集中心主要在人口密集的核心城区。克里金插值结果显示, 聚集中心区县的虫口密度远高于周围其他区县, 说明这些区域是美国白蛾的高度适生区^[9]。平谷区及其附近地区的美国白蛾虫口密度, 虽然低于聚集中心区县, 但仍高于其他非聚集中心的区县。作为美国白蛾在北京地区的首次入侵地, 平谷区可能一直受北京市以外其他地区美国白蛾种群的影响, 不断有美国白蛾传入, 从2012年的三张插值地图可以看出, 平谷区的美国白蛾种群密度显著高于其他地区。与平谷区交界的是河北省承德市兴隆县, 据报道该县经过连续九年的监控防治, 至2014年全县的美国白蛾发生面积和有害株率显著下降^[45], 从图 2的插值地图可以看出, 2013和2014年平谷区的虫口密度明显低于2012年, 因此推测该地区附近分布有一个较大的美国白蛾种群, 能够影响平谷区的美国白蛾种群数量, 进而影响着北京其他区县的美国白蛾种群数量。另外, 通州区和大兴区的美国白蛾种群也不可忽视, 在2015和2016年的插值地图中, 这两个区的美国白蛾种群密度较高, 同样地, 这两个地区也是河北省交界。因此, 在防治美国白蛾的过程中, 不仅要关注当地的危害情况, 也应加强与周边地区的联防联控。

本文运用地统计学方法(GS)和地理信息系统(GIS)对2008—2016年北京地区美国白蛾时空分布变化进行了分析, 研究结果与利用频次分布比较法^[46]、聚集度指标法^[47]等传统统计分析方法的研究结果一致, 均表明美国白蛾在空间上呈聚集分布的特征, 并且聚集程度随种群密度增大而增大。地统计学作为一种空间分析工具, 并不是独立存在的^{[5, 48]}, 将其与GIS结合应用于研究有害生物的空间分布动态格局, 为掌握有害生物扩散机制提供了有力支撑。本研究由于未获取到最近几年北京地区美国白蛾的诱捕数据, 无法与最新的时空分布比较, 同时, 收集的诱捕数据受人为经验影响, 加上统计美国白蛾数量时可能遗漏或重复, 导致插值结果存在一定误差，未来的研究可以从这几方面加以改进。

RDA分析结果显示, 有6年美国白蛾成虫数量变化与气象因子的关系变化是相似的, 而2010、2011、2014年则与其他几年的差异较大。根据各世代数量动态变化(图 3)显示, 2010年第一代和第二代成虫数量明显高于其他年份, 分析发现当年7月份平均温度和8月份最高温度比研究期间平均值高1.1℃、1℃, 充分说明了高温对美国白蛾第一代和第二代是十分有利的。2011年和2014年越冬代成虫数量显著高于其他年份, 这是由于前一年9月份最低温度比研究期间平均值低2℃左右, 可能增强了越冬代蛹滞育强度^[49—50], 并且前一年第二代成虫数量比相邻年份成虫数量高, 一定程度上增加了越冬代蛹的数量。

本文对北京地区美国白蛾成虫数量与气象因子的分析, 未加入降水量因素, 是因为降水量与相对湿度高度相关, 避免数据冗余。美国白蛾喜潮湿环境, 所以在防治过程中, 除了温度、湿度、日照时数, 也应当观测降水变化, 适当的调整防治措施。另外, 7月中旬至9月下旬, 世代重叠现象比较严重^[41], 应采取多种措施防治各龄期的美国白蛾。

4 结论

空间特征上, 美国白蛾种群呈现聚集分布, 在北京地区不同区县之间的分布存在较大的差异, 虫害发生较轻的区县虫口密度只有严重区县的几十分之一；整体上, 美国白蛾集中分布在海拔400m以下的地区, 东部和南部的平原地区较多, 西部和北部的山区少。

时间特征上, 北京地区的美国白蛾种群, 越冬代数量远大于第一代和第二代, 而第一代和第二代种群数量相近, 次年亦是如此, 呈周期性变化。

通过RDA分析气候因素对北京地区美国白蛾种群变化的影响, 结果表明：北京地区美国白蛾种群动态与3月份平均相对湿度、7月份平均气温、上一年8月份最高温度呈正相关关系；与2月份日照时数、4月份最高温度、上一年9月份最低温度呈负相关关系。