基于GIS的环巢湖地区土壤养分空间变异特征研究

童童，梅帅，刘莹，常珺枫，梁华忠，范其龙，王强，马友华*

（1.安徽农业大学资源与环境学院，合肥 230036；2.合肥市农业经济技术服务管理总站，合肥 230000）

粮食是人类赖以生存的资本，随着人类的发展，粮食的需求量也在不断提高。施肥是保障粮食产量的重要手段，但盲目施肥导致养分流失严重[1]，造成河流湖泊等水体污染问题。根据土壤养分状况，合理配施有机无机肥、选择适宜的种植模式，能够有效提高肥料利用率，改善土壤理化性质，减少农业面源污染。土壤养分的空间变异由母质、地形、气象等结构性因素和施肥、种植习惯等人为活动的随机性因素共同影响[2]，研究土壤养分的空间变异特征对充分掌握养分分布规律及促进农业可持续发展具有重要的现实意义。

地统计学中的半变异函数法和Moran's I 指数法已被广泛地应用于空间自相关性的研究[3]。已有学者[4]通过在研究区域选取土壤样点并检测其养分指标，利用地统计学方法和GIS 技术探讨了土壤养分在研究区域的分布特征和空间异质性；陈文举等[5]则通过克里金插值法和单因素方差分析，研究了不同深度土壤pH 的空间分布特征，并了解了土壤酸化分层现象及酸化趋势；Chen 等[6]获取了不同年份的土壤样点数据，通过半方差函数法和克里金插值法，研究了土壤养分的时空变化特征。半变异函数的参数能够表明其空间自相关程度，而Moran's I 指数能够反映出研究对象的分布聚集状况。王强等[7]利用遥感影像数据提取出不同种植年限和种植模式下菜地的空间分布，并以1 km 的网格为评价单元计算局部Moran's I 指数，从而识别出菜地土壤养分空间聚类和空间异常值，最终结果表明菜地土壤养分受到种植时间与城镇距离因素的影响；陈清霞等[8]将半变异函数与Moran's I指数结合，分析出土壤pH 具有小范围、中等程度的空间自相关性，且其空间自相关性显著，体现在土壤pH呈斑块状的镶嵌分布。大多数针对土壤养分空间变异特征的研究是在小范围尺度或大范围尺度低密度采样的数据前提下，其研究结果的精度受限，且研究的土壤养分指标较为单一，少有养分之间全面的对比研究。

本研究以巢湖周边地区（巢湖市、肥东县、肥西县、庐江县、包河区）为研究范围，以同一时期高密度采样下土壤养分（pH、有机质、全氮、碱解氮、有效磷、速效钾、缓效钾）检测结果作为研究数据。基于ArcMap10.6软件，计算并选择最佳半变异函数参数拟合模型，分析其空间自相关程度；结合全局Moran's I指数反映土壤养分自相关性的大小及其是否显著[9]；在最佳拟合模型的基础上，用普通克里金插值法得到养分分布图，并计算局部Moran's I 指数制作LISA 聚类图。通过将半变异函数与全局/局部Moran's I 指数结合，从而得到更加全面有效的研究结果[10]。本研究结果可为环巢湖地区测土配方施肥、养分分区宏观调控和农业面源污染关键源区识别提供理论指导。

1 材料与方法

1.1 研究区概况

巢湖是我国五大淡水湖之一，水域面积780 km2，位于安徽省合肥市境内，属北亚热带湿润季风气候。研究区域位于30°57'~32°16'N，116°40'~117°58'E，占地面积共8 607 km2（包含巢湖水域），耕地面积共3 556.7 km2（巢湖市700.2 km2、肥东县1 054.2 km2、肥西县674.9 km2、庐江县1 105.4 km2、包河区22.0 km2），其中水田占比最大，为3 005.6 km2，占耕地面积的84.51%。研究区域为江淮丘陵地区，地形起伏变化明显，土壤母质主要为下蜀系黄土，土壤类型以水稻土为主。

1.2 样品采集与检测

基于土地利用方式、行政区划和耕地面积划分采样单元，并结合最新的GF 卫星影像数据和谷歌影像数据进行室内点位勘误和纠正，最后通过点位实地踏勘确定采样点位。研究区域共布点8 073 个（图1），其中巢湖市1 632 个、肥东县2 258 个、肥西县1 555个、庐江县2 295 个、包河区333 个。研究区域平均每0.465 km2耕地1个点位。

图1 环巢湖地区采样点分布图Figure 1 Distribution of sampling points around Chaohu Lake

采样时间为2021 年9 月至11 月，采样方法按照《测土配方施肥技术规程》（NY/T 2911—2016），取表层（0~30 cm）样5 点混合，多余样品用四分法舍弃，最终取样量在1~2 kg。对所有土壤样品检测pH、有机质、全氮、碱解氮、有效磷、速效钾和缓效钾，具体指标检测规程和量纲如表1所示。

表1 土壤养分指标检测规程Table 1 Methods for determination of soil nutrients

1.3 研究方法

1.3.1 数据预处理

为确保数据分析的准确性，需要对检测数据进行剔除异常值和正态分布检验，由于该研究样本量较大，剔除异常值采用3 倍标准差法（平均值±3 倍标准差）[11]，异常值用3 倍标准差范围内的最大最小值替代，从而保证样点数不变，并采用频率直方图的偏度和峰度来检验数据是否符合正态分布。

1.3.2 地统计分析

地统计方法是基于半变异函数来研究土壤元素空间变异特征的常用方法，能够揭示随机性和结构性因素对区域化变量的影响程度[12]。半变异函数公式如下：

式中：γ（h）为空间距离为h 时的半方差；N（h）为空间距离为h 的所有样点的点对数；Z（xi+h）和Z（xi）分别为空间位置在xi+h和xi时的实测值。

半变异函数的参数中，块金值（C0）代表人为活动和采样检测误差等因素引起的随机变异，基台值（C0+C）为总空间变异程度，块基比[C0/（C0+C）]为随机变异在总空间变异中的占比，块基比的大小能够衡量空间自相关程度，其变程为空间自相关的理论距离。

1.3.3 空间自相关分析

空间自相关为空间上与相邻区域之间观测值的相互依赖性[13]，全局Moran's I 指数（IN）能够反映整个研究区域的空间自相关性大小，其计算公式如下：

式中：N 为样本数；xi和xj分别代表空间上在i处和j处的实测值；Wij代表空间权重；xˉ为样本平均值。

全局Moran's I 指数在-1 到1 之间，大于0 表示正相关，小于0 表示负相关，其绝对值越大表明其自相关性越明显，全局Moran's I 指数等于0 则表示不相关，空间分布为随机分布。通过Z得分可以检验其显著性，在0.01置信水平下，当|Z|≥2.58时表示为极显著空间自相关。

1.3.4 空间聚类分析

空间聚类分析是基于局部Moran's I指数（Ii）来反映局部空间的自相关性，其计算公式如下：

局部Moran's I 最终将研究区域划分为5 种情况，即①不显著型：无显著空间相关性，呈随机分布；②高-高型：高值聚集区域；③高-低型：低值围绕高值区域；④低-高型：高值围绕低值区域；⑤低-低型：低值聚集区域。通过LISA 聚类图直接反映出局部区域土壤属性的关联性[14]。

2 结果与分析

2.1 环巢湖地区土壤养分含量统计特征

对采样点土壤检测结果剔除异常值后的描述性统计结果及正态分布检验结果如表2 所示。土壤pH、有机质、全氮、碱解氮、有效磷、速效钾和缓效钾的均值分别为6.07、19.97 g·kg-1、1.16 g·kg-1、112.95 mg·kg-1、18.21 mg·kg-1、130.88 mg·kg-1和378.58 mg·kg-1，按照安徽省耕地质量监测指标分级标准，研究区土壤pH 为弱酸性，有机质、全氮、有效磷和速效钾为中等水平，缓效钾较低。根据变异系数小于10%为弱变异、10%~90%为中等变异、大于90%为强变异的划分标准[15]，有效磷为强度变异，其变异系数达到了175.16%，其他土壤养分指标都为中等变异，其中pH变异系数最小，为13.60%。正态分布检验结果显示pH、有机质、全氮和碱解氮符合正态分布，有效磷、速效钾和缓效钾经对数转换后符合正态分布。

表2 环巢湖地区土壤养分描述性统计特征Table 2 Descriptive statistical characteristics of soil nutrients around Chaohu Lake

2.2 环巢湖地区土壤养分空间分布特征

2.2.1 地统计拟合模型的选择

计算并选择不同拟合模型下的半变异函数参数，是进行地统计分析和克里金插值的前提，表3 是不同拟合模型下的交叉验证预测误差。按照预测误差的平均值和标准平均值接近0、均方根最小、标准均方根接近1、平均标准误差接近均方根[11]，选择各土壤养分的最佳拟合模型。结果表明土壤pH、全氮、有效磷、速效钾和缓效钾的最佳拟合模型为指数模型，有机质和碱解氮的最佳拟合模型为高斯模型。

2.2.2 土壤养分空间分布格局

在最佳拟合模型的半变异函数模型及参数的基础上，使用ArcMap10.6 软件通过普通克里金插值法制作环巢湖地区土壤各养分分布图（图2）。综合7项土壤养分来看，环巢湖地区土壤养分都存在明显的南北差异，具有不同程度的破碎化分布，空间变异特征明显。表现在土壤pH 西南部偏酸性，东北部总体偏碱性，但在巢湖水面至东北方向有一条pH 低值条带状分布。土壤有机质、全氮和碱解氮的总体分布较为一致，由东南向西北方向其养分含量逐渐下降，并且在最北部有一块明显低值的三角形区域。有效磷小范围分布破碎化较为明显，整体上西部和北部地区偏低，但在西部和北部中间有一片显著高值区域。速效钾和缓效钾都呈现为中北部显著高于其他区域，从巢湖水面向东北方向有一条明显的高低值分界线。

图2 环巢湖地区土壤养分空间分布Figure 2 Spatial distribution map of soil nutrients around Chaohu Lake

2.3 环巢湖地区土壤养分空间变异特征

2.3.1 土壤养分地统计分析

在地统计分析中，步长大小的选择十分重要，步长过大可能会掩盖短程的自相关性，步长过小会导致落在计算范围内的点对太少，从而影响研究结果[16]。本研究通过计算采样点与最相邻点之间的平均距离（451.04 m），采用略低于该值的步长（450 m）和系统默认的12 个步长数，从而确保所有的步长都会有若干个点对。环巢湖地区土壤养分拟合模型半变异函数参数如表4 所示，其中块基比代表了随机因素引起的空间变异在总变异中的占比，用来表示空间自相关程度，块基比小于25%表示强烈的空间自相关性，25%~75%为中等的空间自相关性，大于75%则空间自相关性较弱[17]。环巢湖地区土壤pH 和缓效钾为中等空间自相关性，表明空间变异由随机性因素（施肥、种植习惯等人为活动）和结构性因素（母质、地形等自然条件）共同影响。其他养分指标（有机质、全氮、碱解氮、有效磷、速效钾）空间自相关性均较弱，表明其变异主要受随机性因素影响。变程代表区域变量的空间变异范围，7 项养分指标的变程依次为碱解氮>有效磷>全氮>有机质>缓效钾>速效钾>pH，碱解氮的变程为4 978 m，说明在4 978 m 的范围内其空间自相关性较弱。pH 的变程最小（1 968 m），说明环境因素在小范围上影响着pH 的空间异质性，较其他土壤养分空间分布更加破碎化[18]。

表4 环巢湖地区土壤养分半变异函数参数与全局Moran′s I指数Table 4 Semivariogram parameters and global Moran′s I index for soil nutrients of around Chaohu Lake

2.3.2 土壤养分空间自相关分析

全局Moran's I 指数反映了区域变量在整体上的聚集状态，其绝对值越大，说明空间自相关性越明显。环巢湖地区土壤各养分均为极显著空间正相关（Mo⁃ran's I>0，P<0.01，Z>2.58），其全局Moran's I指数依次为有效磷>速效钾>pH>缓效钾>全氮>有机质>碱解氮，有效磷的全局Moran's I值为0.711，说明采样点有效磷含量与邻近采样点接近，其空间聚类更为明显。有机质、全氮和碱解氮的全局Moran's I 值较小，表明其在空间上的分布更为分散。

2.3.3 土壤养分空间聚类分析

基于普通克里金插值结果，以行政村为评价单元，将行政村内插值结果的平均值空间连接到评价单元，并进行局部Moran's I指数计算，形成LISA 聚类图（图3）。pH高值聚集区域主要分布在北部和东部，且伴有少部分的低-高型异常区域，低值聚集区域主要分布在南部。有机质、全氮和碱解氮的聚类特征较为分散，其低值聚集区域都主要分布在中北部，高值聚集区域较为破碎化，在南部和中部均有不同程度的分布，同时高-低型和低-高型异常区域在高、低值聚集区域周边都有零散分布。有效磷的空间聚集特征最为明显，低值聚集区域在北部和西部，北部与西部中间有片明显的高值聚集区域，且异常值区域较少。速效钾和缓效钾的高值聚集区域分布在中北部，低值聚集区域分布在中部和南部，且在高低值聚集区周边伴随少量的异常值区域。

图3 环巢湖地区土壤养分LISA聚类分布图Figure 3 LISA map of soil nutrients around Chaohu Lake

3 讨论

3.1 养分变异特征分析与分区管理建议

3.1.1 土壤pH

环巢湖地区土壤pH 的块基比（73.92%）相比较其他养分较小，其pH 的变程（1 968 m）和变异系数（13.60%）在各养分中最低，这与董佳琦等[19]的研究结果相似，表明pH的空间连续性较小，在小范围内具有中等程度的空间自相关性，在土壤pH 分布图上表现为小块斑状聚集分布。研究区整体上自东北至西南方向土壤逐渐偏酸性，这可能与其有机质南高北低的大致分布有关[20]，腐殖质分解的腐植酸会降低土壤pH[21]，针对偏酸性地区应选择适宜的种植模式、保持水土、减少单质化肥的投入[22-23]，并且施用石灰能够调节土壤pH 从而减少重金属的危害[24]。地形地貌能够通过影响成土母质、土壤类型、土地利用方式及水热等条件，从而影响土壤属性，研究区自巢湖水面至东北方向有一条带状分布的土壤pH 低值聚集区，此处为大别山余脉，其山地分布与土壤pH 低值聚集区呈现高度重合，这主要是由于山地地区土壤淋溶作用较强，盐基离子易被氢离子取代成为盐基不饱和土壤，从而导致土壤偏酸性[25]，对于山地酸性土壤地区可适当选择种植果树与烟叶等，并增施有机肥来提高土壤肥力。

3.1.2 土壤有机质、全氮和碱解氮

土壤有机质、全氮和碱解氮的相关性较高，空间分布特征较为相似，其块基比（93.99%、91.39%、90.92%）和变程（3 874、4 467、4 978 m）较高，而全局Moran's I 指数（0.153、0.155、0.107）较小，表明在较大范围内人为活动等随机性因素主导其空间变异，且空间分布较为分散，在其养分分布图中呈现为小块斑状零散分布，这是由于土壤有机质受基施有机肥的影响较大[26]，小农户经营的现状致使种植模式难以统一，从而导致了在一定范围内土壤有机质、全氮和碱解氮分布破碎化。土壤有机质、全氮和碱解氮的变异系数（37.89%、36.91%、44.97%）相对较小，说明整体离散程度较低，具有破碎化严重又向均质化发展的分布特征。因此，在测土配方工作中应适当提高样点密度，从而更精确地掌握该养分分布特征，促进化肥减量增效。在土壤养分指标中，有机质和全氮含量对作物产量影响较大，其与土地利用方式和有机肥的施用相关性较高[27]，针对有机质和全氮偏低的地区实行有机无机配施、作物轮作[28]、秸秆还田和间套种植模式[27]，能够有效改善土壤理化性质并提高有机质与有机碳等含量。

3.1.3 土壤有效磷

有效磷的块基比（85.68%）和变程（4 878 m）与有机质（和全氮、碱解氮）较为接近，说明有效磷在大范围内的强度变异主要由人为活动等随机性因素导致[26]。其变异系数（175.16%）和全局Moran's I 指数（0.711）显著高于其他养分，且空间自相关性极为显著（Z>2.58），表现为高低值聚类特征明显，整体离散程度也最大，综合说明了有效磷具有高度积累的现象。结合LISA 图了解到有效磷在西部与北部有明显的大片低值聚集区域，但在两地中间却有一片高值聚集区域，这可能是中间区域为合肥市城区，人为活动较为集中，导致有效磷富集严重。土壤中的磷主要以吸附态和磷酸盐为主，其移动性较弱，极易被固定，当季利用率较低[29]。针对有效磷集中偏高的市区，应当寻找导致其含量普遍偏高的具体原因，采取“源头控制-过程拦截-末端治理”的防控思路并结合政策调控，从而避免造成农业面源污染[30]；而对有效磷偏低的地区，搭配有机肥与磷肥配施能够提高磷的有效性，并且施用磷肥能够明显提高有效磷含量[31]，同时有机质也在一定程度上能够影响有效磷的含量[32]。

3.1.4 土壤速效钾和缓效钾

速效钾和缓效钾在空间分布上保持着高度统一，但速效钾的块基比（81.77%）明显高于缓效钾（63.20%），说明速效钾的空间变异由人为活动等随机性因素的影响相对较大，而缓效钾更多的由成土母质等结构性因素影响，这是由于土壤中的速效钾为交换性钾和水溶性钾，而缓效钾为黏土矿物固定的钾，因此与成土母质的关系更为密切[33]。同时，速效钾的全局Moran's I 指数（0.397）也高于缓效钾（0.213），体现为速效钾的空间自相关性更为明显，更容易呈现聚集性分布，这是由于当作物吸收而导致速效钾不足时，缓效钾会迅速转化为速效钾，从而维持速效钾的含量。与土壤pH 分布特征类似，速效钾和缓效钾的含量在南部丘陵地区与东北方向的条带状山地区偏低，这是由于海拔较高的低山丘陵区受水流侵蚀作用较强，较土壤氮、磷元素，钾元素更容易流失[13]。作物秸秆中含有大量钾素，对缺钾地区促进秸秆规范化还田，注重施用钾肥能够有效提高土壤钾含量。

3.2 研究创新点与改进思路

本研究结果中，土壤养分空间自相关性总体上较弱，其空间变异更多的是由人为活动等随机性因素导致，且主导了土壤养分变异朝均质化方向发展[19]。在地统计分析结果中，土壤pH更容易受到母质、地形等结构性因素的影响，这与其他地统计学的研究结论一致[8，34]，并且土壤pH、速效钾和缓效钾的分布与地形海拔相关性较高，有许多学者已经证实了这一点[13，35]；有效磷的变异系数较其他养分更大，这与牛文鹏等[36]在多养分中的对比结果一致，有效磷的全局Moran's I 指数最大也证实了其更容易富集的特点。有机质、全氮和碱解氮的块基比最大，这与姜霓雯等[13]的研究中有机质与全氮块基比最大的结果类似，说明了其在土壤养分中更容易受到人为活动等随机性因素影响。从研究数据上看，本研究数据与其他相关研究偏差较大，但养分之间差异的总体规律较为相似，这可能与本研究区布设样点数较多有关。这是由于区域变量的空间变异特征结果随着采样点尺度变化而变化[16]，当采样点的尺度明显大于母质、地形等自然因素分布的尺度时，地统计分析结果更偏向于由结构性因素主导其空间变异[36]，当采样点尺度较小时结果则相反，在部分高密度采样的研究结果中已经证实了这一点[37]。

本研究基于环巢湖地区大范围尺度高密度随机布点采样的检测数据，使用半变异函数参数的块基比和变程分别说明其养分空间自相关程度和作用范围，用全局Moran's I 指数反映其空间自相关性的大小及显著性，并将研究结果与养分统计特征的变异系数呼应，最后通过克里金插值和局部Moran's I 指数制作养分分布图和LISA 聚类图，系统地呈现了研究区多种土壤养分的空间变异特征及内在联系。虽然利用地统计学和Moran's I 指数研究土壤养分空间变异特征的研究较多，但同时使用两种方法交叉进行结果论证的研究却较少，同时本研究涉及到的土壤养分较为全面，进一步对比分析更能为其他土壤养分空间变异相关研究提供可靠的数据参考。本研究详细介绍了环巢湖地区土壤养分的空间变异特征，但未深入探讨养分与种植模式、土壤质地、坡度坡向、海拔和植被类型等多种驱动力之间的关系，随着遥感技术的发展，下一步研究可以通过遥感影像提取或预测这些影响因素，从而更好地了解土壤养分空间分布的关键因素以及发展趋势。

4 结论

（1）环巢湖地区土壤pH、有机质、全氮、碱解氮、有效磷、速效钾和缓效钾的均值分别为6.07、19.97 g·kg-1、1.16 g·kg-1、112.95 mg·kg-1、18.21 mg·kg-1、130.88 mg·kg-1和378.58 mg·kg-1，除有效磷为强度变异外，其他养分均为中等程度变异。

（2）环巢湖地区土壤pH、全氮、有效磷、速效钾和缓效钾的半变异函数最佳拟合模型为指数模型，有机质和碱解氮的最佳拟合模型为高斯模型。pH 和速效钾为中等空间自相关性，其他养分空间自相关性均较弱，表明土壤养分空间变异主要受人为活动等随机性因素影响。

（3）环巢湖地区土壤各养分均为极显著空间正相关，其空间相关性大小依次为有效磷>速效钾>pH>缓效钾>全氮>有机质>碱解氮，且有效磷的空间自相关性显著大于其他养分。

（4）环巢湖地区土壤除有效磷外的其他养分空间分布较为分散，且不同养分空间分布规律不同，总体上看，南北方向各养分含量逐渐上升或下降，其养分高、低值聚集区域也有明显的南北分布。而有效磷聚集特征明显，表现为西部至北部两端低值聚集，中间区域高值聚集。