土壤养分累积迁移成为地表水和地下水的污染源，是全世界迫切解决的难题之一^[1]。农田面源污染中氮磷元素的迁移已导致我国大面积的湖泊和水库富营养化^[2]，受到了普遍关注，一些学者利用多元统计、地统计方法进行定性或定量研究^[3]。一般认为自然来源和人为来源是农田面源污染的两大来源，前者指在成土过程中地质高背景值带来的富集作用，后者主要指农业生产及农村生活对土壤的作用^[4-7]。蔬菜在居民生活中必不可少，我国已是世界上最大的蔬菜生产和消费国，蔬菜播种面积和产量均占世界的40%以上^[8]。根据中国国家统计局数据及公开发表文献，我国蔬菜面积已达到1 998.1万hm²，设施蔬菜面积370万hm^2[9]，蔬菜露天栽培面积1 628.1万hm²。我国蔬菜集约化种植区大量施用肥料，导致土壤养分流失，引起土壤酸化、盐渍化和地下水污染，使土壤生物活性降低、分解能力下降，造成了严重的经济及环境问题^[10]。农业景观空间格局是农业面源污染的主要影响因素之一^[11]，菜地面源污染风险空间分布格局识别与优化可为科学防范和治理农业面源污染提供决策依据^[12-14]。农业景观空间格局主要研究景观格局的时空变化特征^[15]，通常采用遥感、地理信息科学和景观格局分析方法，运用经典回归模型和空间滞后模型研究农业景观格局及其组分的变化，如傅伯杰^[16]对黄土区农业景观空间格局的研究，Liu等^[17]运用多元线性回归模型量化不同土地利用对流域中氮、磷浓度贡献的研究。目前国内外学者在农业景观研究领域已有相关研究，但主要局限于农业景观格局变化及其环境效应、评价方法及其应用等方面^[18]。现代农业要求农田集约化、田面平整化、田块规则化，这促使农田斑块均质，导致土壤生态系统自调节功能减弱，为探明生态风险区来源，需要可追溯污染源的空间分析方法^[19]。局部二元莫兰指数分析可以帮助检测空间异常值或热点，是一种有效的区域变量探索性空间分析方法^[20]。一些学者研究了我国菜地土壤养分的时空变化和环境质量变化^[10]。Valente等^[21]使用局部二元莫兰指数进行空间相关性分析，结果表明土壤电导率与土壤未利用磷呈强相关性，可作为精准农业管理的工具。Fu等^[22]研究发现，由于农民的集约化管理，土壤磷在农场周围和交通路线呈现高-高（High-High）空间分布。然而以往研究中缺乏时空尺度下菜地种植对土壤属性空间分布格局影响的研究，无法阐明研究区长时间尺度下菜地种植对土壤属性值变化的影响机制。因此，本研究以安徽省肥东县为实验区，通过面向对象分类方法获取2019年农业设施菜地空间分布位置，借助地面农业调查得到研究区全部菜地空间分布位置，结合采集的375个表层土壤样品，运用空间插值方法模拟土壤属性含量空间分布，运用莫兰指数对菜地空间密度分布是否影响土壤属性空间分布格局的问题进行研究，研究结果为农业土壤质量安全的管理与保护提供数据支撑，有助于基本农田的长期管理和保护。

1 材料和方法 1.1 研究区概况

肥东县位于安徽省中部，是巢湖流域面源污染优先控制区^[23]，面积22.12万hm²，耕地面积7.67万hm²，占总面积的35%，总人口110万，其中农业人口97万，是一个典型的农业大县。气候属北亚热带季风气候区，四季分明，雨量适中。菜地总面积6 159 hm²，其中露天菜地5 092.52 hm²，种植时间从1989年至2013年，设施菜地1 066.48 hm²，种植时间从2003年至2019年。

1.2 数据来源及技术路线 1.2.1 蔬菜种植类型获取

遥感影像来源于2019年下载的Google earth17级数据，主要根据设施的形态特征如长度、宽度、空间分布等，基于面向对象方法提取影像中的所有设施菜地。运用谷歌地球中的不同历史影像，对所有的设施菜地进行时空对比，将最早出现时间作为起始种植时间。设施菜地提取结果与2016年全县农业调查总菜地空间分布进行空间相交，分别得到不同种植年限露天菜地与设施蔬菜空间分布，如图 1所示。

1.2.2 土壤采样与分析

基于成土母质、土壤类型划分采样单元并简单随机布点采样，于2017年7月至12月采集了375个0~ 20 cm的表层土样，包含24个露天菜地采样点和10个设施菜地采样点，如图 2所示，其分析方法依据测土配方施肥技术规范（2011年修订版）执行。

1.3 地统计插值和核密度模拟

本文采用ARCMAP 10.2软件平台进行土壤属性的地统计学插值和菜地密度的模拟，地统计是依据协方差函数对随机过程/随机场进行空间建模和预测（插值）的一种回归算法^[24]。核密度估计是对离散柱状图连续替换，创建平滑曲线的估值方法^[25]，把空间配置和输入点数量计算的默认搜索半径作为带宽^[26]，避免了稀疏数据集中的空间异常值。

1.4 空间自相关分析

莫兰指数是计算空间自相关中最为经典的方法^[27]，一般分为全局型（Global spatial autocorrelation）和区域型（Local spatial autocorrelation）两种^[28]。全局型Moran′s I是基于统计学相关系数的共变量（Covariance）关系推算得来，共变量的大小程度即代表两组数的相关性大小。全局型Moran′s I值的功能在于描述某现象的整体分布聚集状况，5%显著水平下，Z（I）>1.96时，表示研究范围内某现象的分布有显著空间自相关性，Z（I） < -1.96时，表示研究范围内某现象的分布呈现负的空间自相关性。全局型的Moran′s I的公式如下：

(1)

式中：W_ij是研究范围内每一个空间单元i与j（j={1，2，3，…，n}）区空间单元的空间相邻权重矩阵，以1表示i与j相邻，而以0表示i与j不相邻；x_i为一组变量的每项值，x_j为另一组变量的每项值，x是其平均数。

局域型Moran′s I是空间自相关局部检验统计量，用于识别空间聚类和空间异常值的位置。计算方法如下：

(2)

各变量定义与公式（1）类似。依照公式（2）计算出的Moran′s I值介于-1到1之间，大于0为正相关，小于0为负相关。值越大空间分布相关性越大，即空间聚集分布明显，反之，值越小空间分布相关性越小，而当值趋于0时，即代表此时空间分布呈现随机分布的情形。采用GeoDa软件获取全局和局部莫兰指数空间自相关分类结果^[27]。

2 结果与分析 2.1 菜地土壤属性差异性分析 2.1.1 土壤属性特征分析

根据前期调研，研究区青菜有机肥用量为11.25 t·hm^-2，小白菜有机肥用量为7.5 t·hm^-2，用量最多的肥料品牌为河南莲鑫宝肥业有限公司的复合肥料（硫酸钾型，总养分≥14%，有机质≥45%）。土壤采样数据的分析结果如表 1所示，土壤有效磷的变异系数最高，速效钾、全氮、有机质及pH紧随其后，分别达到29.34%、17.73%、14.07%、12.35%、4.67%。

如表 2所示，距离城镇越近，有机质、全氮含量和pH相对越高，随着距离的增加整体呈下降趋势；有效磷含量和速效钾含量则随着距离的增加呈现先增加后下降的趋势，在2.0~2.5 km处出现最高值。

如表 3所示，从均值可知随着种植年限的增加，露天菜地土壤有机质、全氮、有效磷、速效钾、pH变化明显，种植1~20 a的露天菜地除速效钾外其余土壤属性均低于设施菜地。设施菜地种植时间越短土壤属性指标值越高，而露天菜地种植时间越短全氮、速效钾、速效磷的指标值越高，其他属性指标值规律不明显，露天菜地pH在种植时间为20~30 a时，指标值反而较高。从标准差可知露天菜地土壤属性中，有机质、全氮、速效钾离散程度随种植时间的增加而降低，pH离散程度随种植时间的增加而提高，有效磷离散程度规律不明显；设施菜地土壤属性离散程度随种植年限的增加而降低。

2.1.2 菜地核密度空间分布

由图 3可知，两种菜地都具有空间集聚特征，露天菜地集中区域分布较广，主要在乡镇和城市周边（图 3左图）。设施菜地比露天菜地更加分散，主要集中于肥东县西南区域（图 3右图）。

2.2 菜地土壤属性空间格局分析

如表 4所示，全局Moran′s I值表明各土壤属性均有显著空间自相关性（P < 0.01），全局Moran′s I值排序结果是速效钾>有机质>全氮>pH>有效磷。

二元局域性莫兰指数（Bivariate local Moran′s I）的结果表明土壤属性和两种菜地均有显著空间自相关性（P < 0.01），与土壤有机质、全氮及有效磷均呈空间正相关，与速效钾和pH呈空间负相关（表 5），除速效钾外，其余露天菜地的土壤属性相关性均大于设施菜地。

二元局域Moran′s I可以定位双因素空间自相关性，可确定菜地分布对周边区域土壤属性的影响类型，具体分为4种：（1）High（土壤属性值）-High（菜地核密度值）区域，代表菜地发展集中度高，土壤属性累积与菜地相关；（2）High（土壤属性值）-Low（菜地核密度值）区域，代表菜地发展集中度较低，土壤属性累积较为严重的区域，但与菜地种植无关；（3）Low（土壤属性值）-High（菜地核密度值）区域，代表菜地发展集中度高，无养分累积的区域；（4）Low（土壤属性值）-Low（菜地核密度值）区域，代表菜地产业发展集中度低，土壤属性累积较低的区域。如图 4和图 5所示，两种菜地对5种土壤属性中的空间分布影响相同，有机质和全氮在土壤中均呈现高高（High-High）空间正相关，有效磷和速效钾均呈现低高（Low - High）负相关，pH呈现高低（High-Low）负相关。

3 讨论 3.1 土壤属性特征分析

研究区土壤有效磷的变异系数最高，速效钾、全氮、有机质及pH紧随其后。土壤属性空间变异主要与人为生产或自然环境相关。因农户施肥习惯不同且磷在土壤中难以迁移造成其空间变异最高，钾因容易形成半径较小的离子且迁移性最强而空间变异位居其次，因土壤具有巨大的酸碱缓冲能力，因此土壤pH变异系数最低，而有机质与全氮迁移能力在速效钾与pH之间，这与前人研究结果相似^[29]。

各土壤属性与城镇距离的关系变化，源于土地利用变化，肥东县最早发展蔬菜产业的区域一般位于近郊，随着城市扩张，蔬菜产业种植中心发生位移，有效磷含量和速效钾含量在2~2.5 km处出现最高值。

露天菜地相比设施菜地土壤酸化明显，这与前人研究结果相似^[30]。露天蔬菜种植历史相对较长，部分土壤酸化程度高于设施菜地，pH相对较低，但种植时间最长的露天菜地因距离城市较近，随着城市扩建所用建筑材料的影响，pH值反而较高。设施菜地种植时间越短土壤属性越高，因近10 a内新增菜地多为土地流转后的农业合作社，为追求经济效益，其施肥量远超过露天菜地，二者相差4~10倍，各土壤属性除速效钾外均高于露天菜地。露天菜地种植时间越短全氮、速效钾值越高，其他属性规律不明显，这与农户长期形成的相同施肥习惯有关，但露天菜地地表容易产生径流，具有较高的面源污染风险。研究区露天菜地超20 a种植区域集中在龙塘乡周围，近10~20 a的遍布肥东县各乡镇周围，由于设施种植是露天种植利润的3~4倍^[31]，造成近10 a肥东县设施菜地兴起，露天菜地增长较慢。10 a以上设施菜地主要集中在牌坊回族满族乡内，是因为将设施菜地作为扶贫重点项目发展，而其他乡镇发展缓慢。近10 a设施菜地增长迅速，遍布研究区所有乡镇，但集中于其西南部，这与研究区南部大开发后人口集中有关，但到达六家畈镇后，因距巢湖太近，受到农业环保限制而戛然而止。

3.2 空间格局分析

局部Moran′s I指数对空间离群值敏感，可以对土壤中氮、磷等元素划定区域，可为污染源的识别与控制提供依据^[24]。一些研究表明该方法可得到各土壤养分元素含量分布的“高高”、“低低”聚集区和“低高”、“高低”区域的具体位置，可以分别从乡镇尺度^[32]、行政村尺度^[33]对土壤属性进行管理。

本文与木合塔尔·艾买提等^[34]的研究结果不同，后者结果中土壤pH和缓效钾的含量呈空间正相关，全氮在空间上不具有相关性，因本文以1 km的网格作为评价单元，而后者以行政村为评价单元，导致空间尺度差异，从而造成研究结果不同。本文中设施菜地和露天菜地种植模式不同，对土壤属性空间分布的影响具有一定差异性，但种植模式差异无法解释土壤属性所有的空间变异性。其他学者的研究也说明土壤属性空间分布除受到地形、土壤类型等自然因素的影响之外，还受灌溉能力和施肥方式等人为因素的影响，且人为因素越来越强，越来越复杂^[35]。

4 结论

（1）露天蔬菜种植历史相对较长，距今10~20 a种植面积增长较快，主要集中于各乡镇周边，种植时间越短其全氮、速效钾的平均值越高，种植时间最长区域的pH受到城市开发的影响反而最高。

（2）近10 a设施菜地因利润较高而种植面积增长迅速，种植时间越短土壤属性指标值越高，除速效钾外其余4种属性值均高于露天菜地。距离城镇越近，土壤有机质、全氮和pH的指标值越高，随着距离的增加呈下降趋势。

（3）两种菜地种植与土壤酸化、土壤属性累积相关性显著，与有机质和全氮呈现高高空间正相关，与有效磷和速效钾呈现低高空间负相关，与pH呈高低空间负相关。

（4）通过Moran′s I空间分析方法，可实现菜地对周边土壤属性空间格局影响的精准表达，实现菜地种植对不同区域影响差异的空间可视化，可以为进一步分析土壤属性扩散演化机制提供参考，实现菜地不同区域及尺度的土壤属性分区管理，有助于定义更严格的农业种植管理规范，实现基本农田的长期管理和保护。