基于非参数乘法回归模型的农业景观步甲空间分布模拟

丰思捷，贾 强，陈宝雄，张旭珠，王海锋，刘云慧,4*

(1. 中国农业大学 资源与环境学院，北京 100193; 2. 农业农村部农业生态与资源保护总站，北京 100125; 3. 河北工程大学 园林与生态工程学院，河北 邯郸 056038; 4. 中国农业大学 生物多样性与有机农业北京市重点实验室，北京 100193)

0 引 言

人类土地利用集约化导致的自然生境减少和景观结构改变显著地影响着生物多样性的空间分布，已经成为近代生物多样性丧失的重要原因[1-2]。维持农业景观中生物多样性是保障农业生产所需的诸如传花授粉、害虫控制、养分循环等多种生态服务功能的基础[3-5]。越来越多的研究显示，低集约化的农业景观比自然生境或保护区维持的物种更多，其原因在于低集约化的农业景观具有更高的生产力，可以为物种提供更多的资源[6-7]。然而，农业景观结构如何影响生物多样性，尤其是有益生物的空间分布，迄今为止尚未建立较为精确的空间分布及预测模型。

研究农业景观中物种多样性与景观结构的相关关系，通过构建物种空间精确分布模型，探讨以较小人力、物力了解物种个体(Abundance)、丰富度(Richness)以及多样性(Diversity)的空间分布及其随环境变化的状况，对于评估农业景观环境变化对物种多样性影响，制定合理的农业生物多样性保护战略，保障农业可持续发展具有重要意义。早期研究更多地关注景观结构，包括景观组成异质性(Composition Heterogeneity)和景观配置异质性(Configuration Heterogeneity)，以及与景观生物多样性之间的数量关系[8-11]。近年来，随着地理信息系统(Geographic Information System, GIS)技术以及统计分析中空间分布模型的不断完善[12-13]，涌现了多种物种分布模型(Species Distribution Models, SDMs)，即将物种的分布样本信息和对应的环境变量信息进行关联，得出物种的分布与环境变量之间的关系，并将这种关系应用于所研究的区域，对目标物种的分布进行估计的模型。目前常用的模型有生态位因子分析(Ecological Niche Factor Analysis, ENFA)、最大熵(Maximum Entropy, MaxEnt)、广义线性模型(Generalized Linear Models, GLMs)、广义加法模型(Generalized Additive Models, GAMs)、分类与回归树(Classification and Regression Tree, CART)、推动回归树(Boosted Regression Tree, BRT)、马氏距离(Mahalanobis Distance, MD)、边界函数(Border Function, BF)、随机森林(Random Forest, RF)和多元适应性回归样条(Multivariate Adaptive Regression Splines, MARS)等[14]。这些模型中大多使用简化的统计学模型，缺乏内部因素的相互作用，不能反映种群与栖息地之间的非线性、多因素关系。然而，非参数乘法回归(Non-parametric Multiplicative Regression, NPMR)模型是一个灵活可靠的在多维生态位空间中表示物种响应面的方法，其模拟结果与真实结果高度一致(判定系数R2=0.97)，模拟精度极高，现已广泛用于物种分布的监测和量化[15]。因此，将非参数乘法回归模型与GIS技术结合，对于构建物种分布与环境变量的关系，精确模拟物种空间分布意义重大。

步甲科是鞘翅目最大的科之一，约占昆虫种数的3%。由于其对环境变化敏感，且具有良好的分类学基础和标准化取样方法[16]，常被用作生境变化以及生物多样性研究的指示生物[17-18]。同时，大多数步甲物种在幼虫时期或者整个生命阶段为捕食性[19-20]，在农业生态系统的生物防治中具有重要作用。因此，了解步甲多样性分布格局及其对农业景观结构的响应模式，建立重要农业区域步甲的空间精确分布及预测模型，对农业景观天敌多样性保护，促进害虫生物防治及制定相应的景观管理策略均具有重要意义。

本文采用步甲作为指示生物，以华北平原典型农业景观为例，在调查农业景观中步甲群落多样性及其周围景观结构特征的基础上，定量分析景观结构与步甲物种多样性、步甲常见种分布之间的关系，构建物种-景观结构关系模型，并结合GIS技术建立研究区农业景观中步甲群落常见种空间分布预测模型，评估农业景观结构变化对物种多样性的影响，以期为华北平原农业景观步甲多样性保护和制定合理土地利用规划提供参考和技术支持。

1 材料与方法

1.1 研究区概况

研究区位于北京市密云区东南部(40°21′23″N～40°31′23″N，116°49′06″E～16°41′13″E)(图1)。该区属温带半湿润大陆性季风气候区，四季分明，温差较大，光照充足，盛行东北、西南风。年平均气温为10.8 ℃，年平均降水量为661.3 mm，无霜期平均为177 d。取样点主要分布于密云区西田各庄镇南部平原地区村庄，包括疃里、建新、大辛庄、沿村、水洼屯、恒河、渤海寨及西田各庄等8个村在内约35 km2的区域。根据2009年景观调查数据，农用地为该区域主要土地利用类型，面积为20.2 km2，占总面积的59.4%，其中玉米地、花生地、果园和半自然生境(林地、荒草地、防护林带)分别占研究区土地面积的41.2%、9.7%、6.8%和21.9%。

1.2 取样方法

选取研究区主要土地利用类型为玉米地40块(图1)，每块样地设置20 m×20 m的样方一个，采用五点法进行取样，陷阱分别置于样方中心及其与东、南、西、北4个方向对角线中部[图2(a)]。所用陷阱为上口径8.2 cm、下口径5.2 cm、高11.7 cm的塑料口杯，其上离地面5 cm处设置10 cm×10 cm的防雨罩[图2(b)]。每个陷阱中放置约占陷阱堤1/3的75%酒精作为防腐剂。取样时间为2009年5月至9月，每月中旬进行一次物种取样，每次布点6 d后收回[21]。取回陷阱时，每个样地的5个陷阱进行合并，用75%酒精保存回实验室分类鉴定到种。

图1 研究区域取样点分布和土地利用示意图

图2 步甲取样样方设计

1.3 景观分类

通过遥感图像处理软件ENVI分别将2006年8月21日和2008年8月31日的QuickBird遥感影像进行拼接并几何校正。步甲取样期间，结合遥感影像对当地土地利用情况进行田间调研，并在ArcGIS中数字化处理，绘制研究区土地利用图。

根据Walz的研究，与生物多样性关系紧密的景观要素包括半自然生境比例、斑块大小、斑块多样性、结构多样性、连通性以及地貌多样性[22]。而步甲的活动主要在50 m的范围内[23]。基于此，在ArcGIS9.2软件中，以取样点为中心，提取出边长为200 m的正方形景观，并在Fragstats3.3软件中计算斑块密度(Patch Density, PD)、香农多样性指数(Shannon Diversity Index, SHDI)、周长面积比(Perimeter-area Ratio, PARA)、聚集度(Contagion, CONTAG)、斑块丰度(Patch Richness, PR)、香农均匀度指数(Shannon Evenness Index, SHEI)[24]。同时，通过ArcGIS9.2软件自身的统计功能，计算半自然生境比例(Ratio of Semi-natural Habitat, RSNH)、道路比例(Ratio of Road, RR)以及半自然生境类型数(Number of Semi-natural Habitat Types, NSNHT)。最后，将上述参数作为景观结构指数。

1.4 数据分析

研究采用非参数乘法回归模型来拟合每个步甲物种和步甲整体群落对景观结构的响应。在非参数乘法回归模型中，自变量以乘法而非加法整合到模型中，即某个自变量的变化对因变量产生的效果也受到其他自变量的影响[15,25]。HyperNiche软件是非参数乘法回归模型拟合的高效实用工具，研究利用该软件进行非参数乘法回归模型拟合并输出物种响应曲线/面[26]；在此过程中，核函数用于平滑，而交叉检验用于选择拟合效果最好同时最简洁的模型。

1.4.1 模型构建

将5个月的取样数据合并，计算步甲群落物种数、活动密度(每个陷阱捕获的步甲数量)、Fisher’s α多样性指数以及香农多样性指数作为群落多样性指数。根据每个物种在各样点物种存在或不存在的状况分别整理为1或0二元数据。通过HyperNiche软件中局部均值(Local Mean, LM)法和局部最小二乘(Local Least Square, LLR)法(局部均值法使用Gaussian核函数，局部最小二乘法使用局部线性估计)，分别分析物种1或0二元数据和群落多样性指数与景观结构指数之间的关系。

1.4.2 模型选择与评价

研究引入logB来指示二元数据响应模型拟合结果的质量。当logB>2.00时，自变量和因变量之间的关系被认为是“非常确定”，模型精度较高。对于定量数据，引入xR2(CrossRSquare)来检验模型拟合的质量[26]。选出合适的模型之后对模型进行“微调”(Tune)，该过程为每个变量调整公差范围，即1%变量值的范围。这些更精细的调试在原模型中4%公差值的范围内发生。当一个模型的平均邻域小于8时，就不对其进行建立模型推理下限的分析[27]。

1.5 步甲空间分布预测

在ArcGIS9.2软件中，将研究区划分为200 m×200 m的若干格网，按照第1.3节中所述的方法分别计算每个方格的景观结构指数，作为预测矩阵输入到HyperNiche软件中，根据非参数乘法回归模型拟合后挑选出来的模型，对每个方格的物种出现概率以及群落多样性指数进行预测。最后，将每个物种的出现概率通过ArcGIS9.2软件出图展示，同时将受景观结构显著影响的物种出现概率进行累加，得到研究区步甲常见物种出现概率分布图。

2 结果分析

2.1 步甲群落物种组成

2009年5月至9月对北京市密云区西田各庄南部农田区步甲进行取样调查，共捕获步甲6 286头，分属10个属21个种(表1)。其中，草原婪步甲、蠋步甲、后斑青步甲、通缘步甲和淡鞘婪步甲为研究区步甲优势种，分别占捕获步甲总个体数的65.19%、10.37%、7.92%、7.13%和5.46%；此外，黄斑青步甲占捕获步甲总个体数的1.75%，为该地区的步甲常见种；其余16种步甲个体数量占比均在1%以下，合计占捕获步甲总个体数的2.18%，为该地区稀有种。总体来看，研究区域物种群落结构往往出现极个别物种个体数占绝对优势，而多于一半的物种捕获量较为稀少的状况。

表1 研究区域玉米地步甲群落的物种组成

2.2 步甲物种和群落对景观结构的响应

2.2.1 物种水平上步甲对景观结构的响应

去除在所有取样点中都出现的物种以及出现次数过少的物种之后，将物种的1或0二元数据作为因变量，景观结构指数作为自变量，以非参数乘法回归模型中局部均值法进行拟合，得到剩余14个物种-景观结构模型。结果显示：有6个步甲物种的出现概率与景观结构之间关系“非常确定”，即满足logB>2.00(表2)。这几个步甲物种包括3个优势种(后斑青步甲、蠋步甲、淡鞘婪步甲)、1个常见种(黄斑青步甲)以及2个稀有种(单齿婪步甲、单齿蝼步甲)。具体表现为：在200 m × 200 m景观尺度下，随景观中斑块密度的增加，黄斑青步甲出现概率显著降低(p<0.05)，而淡鞘婪步甲出现概率显著增加(p<0.05)；随着景观中半自然生境类型数的增加，黄斑青步甲和后斑青步甲出现概率显著降低(p<0.05)，而单齿婪步甲出现概率显著增加(p<0.05)；随着景观中半自然生境比例的增加，单齿蝼步甲出现概率显著增加(p<0.05)。此外，由于非参数乘法回归模型拟合中自变量之间会相互影响，蠋步甲出现概率受景观结构的影响较为复杂。

表2 步甲物种存在对景观结构的响应

2.2.2 群落水平上步甲多样性对景观结构的响应

以群落多样性指数作为因变量，景观结构指数作为自变量，以非参数乘法回归中局部最小二乘法进行拟合。结果表明：在200 m×200 m景观尺度上，景观结构对步甲群落物种数、活动密度、Fisher’s α多样性指数及香农多样性指数均无显著影响(p>0.05)(表3)。

表3 步甲群落多样性对景观结构的响应

2.3 步甲物种出现概率预测

根据“微调”之后的6个物种模型，对研究区域黄斑青步甲、后斑青步甲、蠋步甲、淡鞘婪步甲、单齿婪步甲和单齿蝼步甲在各个景观栅格中的出现概率进行预测，并在ArcGIS9.2软件中进行展示。结果显示：黄斑青步甲、后斑青步甲、蠋步甲、淡鞘婪步甲较为普遍地出现在研究区域中，而在研究区域西北部出现较为明显的空白；单齿婪步甲受景观结构影响较大，仅在特定区域才有较高出现频率；而单齿蝼步甲较为均匀地分布在整个研究区域(图3)。

图3 物种分布的概率预测

将受景观影响显著的6个物种预测结果叠加，即：物种在某位点的出现概率若大于0.5，则记作“1”，反之记作“0”；然后将该位点6个物种出现的数据相加，得到这6个对景观变化敏感物种的空间分布梯度[图4(a)]。取出物种分布概率梯度大于4的样点，得到景观变化敏感物种多样性的空间分布[图4(b)]。结果显示：就模型中6个受景观结构显著影响的物种而言，出现概率比较高的地区位于半自然植被和农业用地之间的狭长地带，或者位于农业用地中半自然植被附近，而单纯半自然植被比例较高的地区、无半自然植被覆盖的农用地中或建筑用地中，步甲出现概率均较低[图4(b)]。因此，若要对研究区域步甲物种进行保护，需要重点保护这些农田中的半自然植被，同时也需要合理配置农田中半自然植被在空间上的分布。

图4 6个对景观结构敏感物种的空间分布梯度

3 讨 论

根据本研究中非参数乘法回归建模的结果来看，200 m×200 m尺度下，步甲群落多样性指标与景观结构指标之间的关系并不显著。对物种存在/不存在与景观结构关系的模拟显示，在捕获的23个物种中，9个物种分布范围过于广泛或者狭窄，8个物种分布对景观结构的响应并不显著，仅有6个物种与景观结构指标建立了“非常确定”的关系。

非参数乘法回归模型适用于模拟物种分布的概率且模拟精度较高。在适应性方面，Safaei等利用逻辑回归(Logistic Regression, LR)、非参数乘法回归和生态位因子分析等模型研究伊朗中部黄芪(AstragalusverusOlivier)的潜在分布情况，发现虽然生态位因子分析模型精度最高，但其强调的是环境适宜性值，非参数乘法回归和逻辑回归模型则关注的是物种的发生概率，因而更适用于物种模拟[27]。此外，相比逻辑回归模型，非参数乘法回归模型较简单，且提供了物种对各种环境梯度的响应曲线，从而量化了物种存在的最佳条件。在精度方面，Lintz等通过检验模拟面的二元和连续响应的预测误差来比较分类与回归树、随机森林和非参数乘法回归模型的预测精度，发现与其他模型相比，非参数乘法回归模型显示了最小的可变性和最大的准确性[28]。本次研究利用非参数乘法回归模型模拟步甲物种对农田景观结构的响应，显示logB>2.00，因此，利用非参数乘法回归模型可以很好地模拟和预测步甲在景观中的空间分布变化。

欧洲研究显示，步甲对景观结构的响应可能因物种的特征不同而不同，例如广布型(Generalist)步甲的物种丰度和多度显示与景观结构显著相关，但是专性分布(Specialist)的步甲不受景观结构的影响[29]；此外，物种迁移能力、体长等的不同也可能使其对景观变化的敏感程度不同[11,30]。本次研究用于建模的14个物种中，仅有6个物种极显著地受到景观结构的影响，可能原因有3种。第一，未建立显著相关关系的8个步甲物种多数为植食性或兼食性物种(如婪步甲属、暗步甲属)，属于较低营养级步甲物种；而与景观结构建立显著相关关系的物种多数为捕食性物种(如青步甲属、蠋步甲)，属于较高营养级步甲。早期研究证明，位于食物网较高营养级的物种对于周围景观结构的改变更加敏感，而位于较低营养级的物种往往不受周围景观结构的显著影响[11,31]。第二，本次研究取样只考虑了玉米地。一方面，取到样本较多的物种可能是已很好地适应农田生境及其管理方式的物种，因而这些物种较少受周围景观结构的影响；另一方面，对于某些更倾向于在自然、半自然生境中生存的物种(如暗步甲属、麻步甲等)，并不能单纯依据农田取样结果建立适宜于整个景观的物种空间分布模型。第三，不同物种对周围景观结构响应的尺度不同[32]，因此，取样点200 m范围的景观尺度可能因为物种之间的差异化响应而对群落水平的多样性不表现显著影响。此外，在群落水平上，步甲多样性不受景观结构影响也可能是由于不同物种对于景观结构指标的响应方式区别很大甚至相反，从而表现出群落水平步甲的多样性不受景观结构指标的影响。

在受景观结构指标显著影响的6个物种中，由模型可知影响这些步甲物种分布的景观指数包括景观斑块密度、半自然生境类型数、道路比例、香农多样性指数、半自然生境比例以及聚集度。然而，不同物种对于周围景观结构的响应不同。在200 m×200 m景观尺度下，半自然生境类型数增加，可能意味着景观中小斑块数量增加，景观破碎化严重，从而使景观中捕食性步甲(如黄斑青步甲和后斑青步甲)的出现概率显著降低；而在此景观尺度上，斑块密度的增加同样意味着景观破碎化程度的增加，因此，黄斑青步甲的出现概率会随之降低。此外，半自然生境类型增加、景观中斑块密度增加，可以为以草本科种子作为补充食源的兼食性步甲(如淡鞘婪步甲、单齿婪步甲)提供更加丰富的替代食源，从而增加其在景观中的出现概率。

在欧洲，半自然生境在农业景观生物多样性保护和生态系统服务维持中的重要功能和作用成为研究和关注的焦点[33-35]。一般结论认为，农业景观中的半自然生境(包括农田边界、灌木带、林地、水塘、沟渠和休耕地等)的保护和建立，可以为农业景观中绝大部分物种提供生存所需的资源，如食物源、物种源、避难所、繁育场所和迁移廊道等，有利于农田景观物种多样性的维持和保护[36-39]，具有较高自然、半自然生境比例的景观中，农业景观具有较高天敌昆虫多样性和较少害虫爆发[40-41]。本次研究同样发现，增加景观中半自然生境的比例，可以显著增加景观中捕食性步甲(单齿蝼步甲)的出现概率。

本次研究的重要性在于，通过建立随机取样点的物种-景观数量关系模型，在ArcGIS9.2软件中对周围整个景观中步甲的出现概率做出了空间分布预测。从发现的6个与景观结构有明确关系的物种多样性分布梯度来看，这类物种出现概率较高的区域主要位于半自然生境与农用地之间或者农用地中靠近半自然生境附近，而单纯的半自然生境中(如西北部大片的半自然生境)或周围半自然生境较少的大片农用地中步甲的分布梯度均不高[图4(a)、(b)]。经过长期的人为干扰和集约化生产，现有景观中的物种可能主要为已经适应目前集约化生产的“农田物种”；农田较高的生产力会为其提供较为丰富的食物来源，因而会在农田中而非在半自然生境中具有较高的多样性，但是半自然生境较农田所受干扰较小，在为农田生物提供替代的食物来源、避难所、繁殖地、越冬场所等方面起到了重要作用[42-44]；步甲在农田和半自然生境间不断的迁入、迁出，因而可能通过生境间溢出(Cross-habitat Spillover)效应使得半自然生境与农用地之间或者农用地中靠近半自然生境附近具有较高的多样性[45-46]。因此，增加农田地块间半自然生境对增加步甲群落的多样性将具有积极意义。然而，农田管理过程中，应尽量避免景观过于破碎化而造成的农田生物移动受到阻碍，在增加农业景观半自然生境比例的同时，合理优化景观配置异质性，从而达到增加农业景观有益生物多样性保护的目的。

4 结 语

(1)利用非参数乘法回归模型可以很好地模拟步甲物种对农田景观结构的响应，模拟和预测步甲在景观中的空间分布变化。

(2)200 m×200 m尺度下，尽管步甲群落多样性不受景观结构的显著影响，但是存在部分步甲物种对景观结构响应显著，但不同物种对景观结构的响应不同。

(3)增加农田地块间半自然生境，并合理改善景观配置异质性将有利于对景观结构变化敏感步甲物种的维持。