蚊虫是重要的病媒生物，它的生态习性是刺叮吸血，从而传播多种传染病。由于全球气候、环境、交通、城市化的变化以及耐药性等因素，蚊虫造成的人类疾病负担不断上升，同时对人居环境造成一定的影响^[1-7]。我们要因地制宜地开展蚊虫防制工作，就要对蚊虫空间分布及其影响因素有充分的了解和认识。本研究利用地理信息系统相关知识对南京蚊虫空间分布进行预测，并对影响蚊虫的土地利用驱动因素进行分析研究，为进一步防控蚊媒病和开展爱国卫生运动提供新的思路和方法。

1 材料和方法

1.1 材料

蚊虫监测数据：2018年3—11月，在南京除高淳以外的行政区设置固定的蚊虫监测点（图1）。采用国家标准GB/T23797⁃2009《病媒生物密度监测方法蚊虫》中推荐的的诱蚊灯法，参照同类研究设置城镇居民区、公园（含街心公园）、医院、农村民房和牲畜棚（猪圈、牛棚、羊圈、养殖场等）5类生境共96个监测点，每月开展2次监测收集，每个监测点设置1台诱蚊灯，光源离地1.5m，日落前1h接通电源，次日日出后1h拔出电源，取出捕蚊袋，对蚊虫进行冷冻处理或乙醚麻醉，鉴定计数，记录监测时间、地点 （经纬度），在江苏省疾病预防控制业务集成平台中的病媒生物监测网络直报平台进行数据汇总，其研究变量为蚊密度，计算公式为D=Nm/（Nl·T）^[8]，式中：D为蚊虫密度，单位为只（/灯·h），定义为每灯每小时捕获的蚊虫数量；Nm为蚊虫数量，单位为只；Nl为诱蚊灯的数量，单位为灯；T为诱蚊小时，单位为h。

遥感数据：基于2018年Landsat5遥感影像30m栅格土地利用数据（http：//ids.ceode.ac.cn/），根据影像光谱特征、结合野外实测资料，同时参照有关地理图件，对地物的几何形状、颜色特征、纹理特征和空间分布情况进行分析，建立包括耕地、林地、草地、水域、建设用地和未利用地在内的6个一级类型土地和包括林地、灌木林、疏木林、其他林地和高、中、低覆盖度草地等25个二级类型土地，具体分级分类和含义见表1，土地利用空间分布情况见图1。

1.2 方法

1.2.1 空间分布预测

克里金法是20世纪60年代法国Matheron在前人基础上总结并提出来的地统计方法，是在变异函数空间分析的基础上，对有限区域内的区域化变量取值进行无偏最优估计的一种方法，具有插值和估计的双重特点，它包括普通克里金、简单克里金、泛克里金和经验贝叶斯克里金等^[9-10]。经验贝叶斯克里金法可自动执行构建子集和模拟过程来设置内在参数，而不像其他克里金方法需要手动设置，其半变异函数参数使用受限最大似然法估计，而其他克里金方法通过已知的数据来计算半变异函数^[11-12]。

1.2.2 空间预测评价

对模型拟合评价方法交叉分析采用的是留一交叉验证法（leave⁃one⁃out cross validation，LOO⁃CV），是指假设有个N个采样点，对于第i个训练样本，将其取出，对剩下的N-1个样本进行训练，得到决策函数，并用其测试第i个训练样本，重复N次得到N个模型，用这N次交叉验证得到的预测误差来进行评价，其评价指标是预测误差的均值（MAE）、均方根误差 （RMSE）以及决定系数（R²）^[13]，计算公式如下：

$\begin{array}{c}MAE=\frac{1}{n}\sum _{i=1}^n\left|f_i-y_i\right|\\ RMSE=\sqrt[2]{\frac{1}{n}\sum _{i=1}^n{\left(f_i-y_i\right)}^2}\\ R^2=1-\frac{\sum _{i=1}^n{\left(f_i=y_i\right)}^2}{\sum _{i=1}^n{\left({\stackrel{-}{y}}_i-y_i\right)}^2}\end{array}$

式中y_i是实测值，f_i是预测值，${\stackrel{-}{y}}_i$是实测值的平均值。

1.2.3 驱动因素分析

地理探测器是中科院王劲峰等提出来的空间分析模型，其核心思想是如果某个自变量对某个因变量有影响，那么自变量的空间分布和因变量的空间分布应该趋于一致。那么探测X在多大程度解释 Y的空间分异可以用q值度量^[14-17]，其表达式为：

$\begin{array}{c}q=1-\frac{\sum _{h=1}^L{N}_h{\sigma }_h{}^2}{N{\sigma }^2}=1-\frac{SSW}{SST}\\ SSW=\sum _{h=1}^L{N}_h{\sigma }_h{}^2\\ SST=N{\sigma }^2\end{array}$

式中h=1，...，L 为变量 Y 或因子 X 的分层，即分类或分区；N_h和 N 分别是层 h 和全区的单元数；${\sigma }_h{}^2$和 σ² 分别是层h和全区Y值的方差。SSW和SST分别是层内方差之和和全区总方差。q 的取值为[0， 1]，值越大表明因素X对Y空间分布影响越大，反之则越小。

地理探测器主要分为因子探测器、风险探测器、交互探测器和生态探测器4个部分，本研究采用其中的因子探测器、交互探测器和生态探测器，而风险探测器主要探测等级或定性因素，不适用于本研究^[17]。交互探测器是识别不同风险因子之间的交互作用，即评估不同因子共同作用时是否会增加或减弱对因变量Y的解释力，分别计算两种因子对因变量的q值，并且计算交互时的q值，根据3个值的大小判断交互作用的类型。生态探测是比较不同因子对因变量的空间分布的影响是否有显著差异，用F值来衡量，其表达式为：

$\begin{array}{c}F=\frac{N_{x1}\left(N_{x2}-1\right){\mathrm{S}\mathrm{S}\mathrm{W}}_{x1}}{N_{x2}\left(N_{x1}-1\right){\mathrm{S}\mathrm{S}\mathrm{W}}_{x2}}\\ {\mathrm{S}\mathrm{S}\mathrm{W}}_{x1}=\sum _{h=1}^{L1}{\mathrm{N}}_h{\sigma }_h{}^2,{\mathrm{S}\mathrm{S}\mathrm{W}}_{x2}=\sum _{h=1}^{12}{\mathrm{N}}_h{\sigma }_h{}^2\end{array}$

式中 N_x1和 N_x2分别表示两个因子x1和x2的样本量；SSW_x1和SSW_x2分别表示由x1和x2形成的分层的层内方差之和；L1和L2分别表示变量x1和x2分层数目。

1.3 统计学方法

蚊虫数据和监测点2km范围内土地利用数据采用Excel 2007进行数据整理和汇总。空间分布及预测运用ArcGIS 10.4软件自带地统计学分析模块，交叉分析采用R软件4.0.3 ^[18] 进行计算。土地利用驱动因素分析采用地理探测器软件。

本研究首先计算各监测点的蚊虫密度，同时利用手持GPS仪记录各监测点的经纬度，基于国家基础地理信息中心（http：//www.ngcc.cn/ngcc/）的矢量地图，绘制南京地区地图以及标记监测点信息，利用经验贝叶斯克里金插值对蚊虫空间分布进行预测，并对预测结果进行交叉分析；然后根据蚊虫生态学知识，基于蚊虫监测点2km缓冲区范围内土地利用方式的构成^[19-21]，利用地理探测器模型进行潜在的土地利用驱动因素分析，探索影响蚊虫密度的土地利用驱动因素的相对重要性、是否存在差异及它们之间是否存在交互作用。P< 0.05为差异有统计学意义。

2 结果

2.1 蚊虫空间分布预测

蚊虫监测结果：2018年96个监测点共捕获蚊虫24 989只，所有监测点平均蚊密度为1.21只（/灯·h）。其中密度最低的是鼓楼区湖南路街道监测点，蚊密度为0.02只（/灯·h），密度最高的是六合区龙池街道监测点，蚊密度为8.10只（/灯·h）。

蚊虫分布预测：基于ArcGIS的Geostatistical Analyst模块，利用经验贝叶斯克里金技术进行插值 （图2）。根据预测结果，可以看出南京市蚊虫空间分布整体呈北高南低的趋势，预测的蚊虫密度高风险区域[密度值大于2只（/灯·h）]位于六合区，浦口区、江北新区、栖霞区处于中风险区域[1~2只（/灯·h）]，玄武区、秦淮区、建邺区、鼓楼区、雨花台区、江宁区、溧水区为低风险区域[小于1只（/灯·h）]。2018年高淳区开展蚊虫监测时监测点不固定，每月在不同的街道（镇）开展监测，全年所有监测点平均密度为0.14只（/灯·h），与之相邻的溧水区全年蚊虫监测点平均密度为0.30只（/灯·h），预测结果也提示高淳区处于低风险地区，预测结果基本符合实际情况。

交叉验证：根据实测值和预测值之间的预测误差，计算MAE=0.005 263 883，RMSE=0.741 404 1， R²=0.924 3，整体拟合效果较好。

2.2 土地利用驱动因素分析

土地利用情况：利用ArcGIS在蚊虫监测点设置2km的缓冲区，提取缓冲区内土地利用数据。96个监测点的土地利用按第一构成比排名，城镇用地最多（62个），其次为水田（22个）、旱地（9个）、有林地 （2个）、水库坑塘（1个）。

因子探测器：因子探测器q统计结果见表2，按q 统计大小排序为：水田、农村居民点、城镇用地、水库坑塘、旱地、河渠、有林地、疏林地、湖泊、其他林地、滩地、其他建设用地、裸岩石质地，其中水田、农村居民点、城镇用地、水库坑塘q值大于0.9，并且这4类因子P值均小于0.05。

交互探测器：大部分土地利用方式交互作用为非线性增强，但旱地与滩地、河渠与疏林地交互作用为双因子增强，水田与旱地交互作用为非线性减弱（表3）。

生态探测器：水田、旱地、水库坑塘分别与其他土地利用方式探测结果均不同（表3）。

3 讨论

在以往的蚊虫监测研究中，一般采取在区域内开展长期连续系统的监测，选择城镇居民区、农村居民区、牲畜棚等特殊生境，利用这些监测点数据的算数平均数来代表整个区域的蚊虫密度，这种研究方法可操作性强，方便快捷，可以快速评估区域内蚊虫侵害及蚊媒病传播的风险，但是对于区域内部蚊虫风险的不平衡性和影响蚊虫分布的环境因素无法量化分析，制约了蚊虫监测研究向深度和广度扩展的能力。因此本研究首次在南京地区开展利用空间插值方法来预测蚊虫空间分布，进而对南京蚊虫分布的风险等级有了新的认识。ArcGIS集成的空间分析方法已经广泛应用在公共卫生监测、突发疫情处置、新发传染病的分析等各方面^[22-24]，本研究采用ArcGIS中的经验贝叶斯克里格方法，操作比较简单，需要手动设置的区域较少，预测标准误差比其他克里金插值方法更小，可准确预测一般程度上不稳定的数据，对于小型数据集，也比其他克里金法更准确^[25-26]。交叉分析也显示拟合效果较好，能直观展示蚊虫空间分布现状，对于开展爱国卫生运动和政府购买有害生物防制服务具有一定的指导意义。

交互探测器：前标^*，表示双因子增强；前标^△，表示非线性减弱；前标^#，表示独立；前标无，表示非线性增强；生态探测器：后标^§，表示两两之间方差有统计学差异；后标无，表示两两之间方差无统计学差异。

通过地理探测器的因子探测研究，可以看到影响蚊虫密度的驱动因子分别是水田、农村居民点、城镇用地和水库坑塘。其中水田、水库坑塘是静态水体，蚊虫生活史4个环节中的3个环节都离不开水，水田和水库坑塘一般也认为是按蚊和库蚊孳生繁殖的生境^[27]。本研究分析得出，影响蚊虫密度的第一驱动因子是水田，刘美德等^[28] 通过GIS分析云南自然村中蚊虫群落与环境因子关系的研究中也发现，水稻田面积对于样点蚊虫群落最有影响，且水稻田面积与蚊虫群落的种类数量、个体数、多样性和优势度均有统计学关系。而农村居民点和城镇用地均是蚊虫活动的人居环境，可能存在大量的人工积水容器，是伊蚊孳生繁殖的生境^[27]，同时在农村居民点和城镇用地存在蚊虫刺叮吸血的动物目标和人群目标^[27]。

通过交互作用探测分析，旱地与滩地为双因子增强，可能因为滩地存在大量积水，但易受水位影响，蚊虫孳生环境不稳定，旱地缺少水域，依靠外部水源供给，因此双方交互为双因子增强^[29]。河渠与疏林地交互作用也为双因子增强，可能因为河渠水为流动性水域，疏林地存水功能较为薄弱，河渠水引入疏林地后，对于蚊虫栖息繁殖均有利^[30]。水田与旱地交互作用为非线性减弱，水田一般是蚊虫栖息繁殖的重要场所，旱地部分有积水区域也会存在蚊虫孳生，水田交互旱地后可能减少了作为蚊虫孳生繁殖生境的能力，另外大量水进入旱地后，破坏了原有的蚊虫孳生环境，因此双方交互为非线性减弱^[31]。

通过对本区域的蚊虫密度的土地利用驱动因素分析，建议在城乡规划以及人居环境建设时考虑土地利用模式对蚊虫的作用，例如区域内开展蚊虫防制等爱国卫生运动时，不仅要考虑物理防制与化学防制相结合，更需要着眼于常态化手段和精准性措施，从生态、环境、资源等多角度出发，对影响蚊虫密度的土地利用驱动因素进行有针对性的“治本清源”，如水田采用间歇灌溉、适时烤田、合理密植等^[32]，水库坑塘营建蚊虫天敌栖息地、净化水体，合理选择和配置植物等^[30]，对农村居民点和城镇开展翻盘倒罐和健康教育等^[33]。

本研究不足之处：第一，是缺少高淳区的数据，因此在插值时影响了整个南京地区的蚊虫空间分布预测。根据2018年监测情况，因高淳区的监测点不固定，每月在不同的街道（镇）开展监测，故没有纳入插值模型中，但计算高淳区全年所有监测点平均密度为0.14只（/灯·h），与之相邻的溧水区全年蚊虫监测点平均密度为0.30只（/灯·h），预测结果也提示高淳区处于低风险地区，预测结果基本符合实际情况，对高淳处于低风险的预测是相对科学准确的。第二，本研究只探索蚊虫的土地利用驱动因素影响，属于流行病学的生态学研究。第三，由于时间原因，只对1年的蚊虫密度和土地利用进行研究。未来可以在研究中加入时间趋势对蚊虫生态学进行时空加权回归研究^[34]，另外可以细化不同类型（如田塘型、缓流型、坑洼型、地下型、容器型）蚊虫的土地利用模式，从而丰富蚊虫环境生态学研究。

参考文献