鄂西北山区耕层土壤全氮空间变异及自相关分析

吴正祥,冯 莉

(1.南阳师范学院 乡村振兴研究院,河南 南阳 473061;2.河南省自然灾害遥感监测重点实验室,河南 南阳 473061)

氮素作为土壤组成的重要部分,直接影响土壤肥力,是植物生长发育过程中最重要的元素之一[1-2],含量过高或过低都会引起作物生长不良、水体环境污染等一系列农业问题[3]。土壤中氮的长期积累和循环不仅影响整个生态系统的平衡和稳定,而且对缓解全球温室效应也产生着深远影响[4]。土壤学、生态学、环境科学等不同领域对土壤养分供应、植物养分吸收、温室气体排放等的研究中均需测定土壤全氮(Soil total nitrogen,STN)含量[5]。全氮在土壤养分总量中占比不大,但在土壤微生态、养分循环以及农业可持续发展等方面却有着巨大作用[6],其含量与分布是评价土壤肥力和质量的关键指标[7]。掌握STN的空间变异规律和分布特征可以为精准施肥等农业管理措施、耕地质量的保护和提升提供理论依据[8]。目前,众多学者结合GIS和地统计学对STN空间分异以及STN时空分异进行了较为深入的研究,取得了丰硕的成果[9-11]。学者通过采用方差分析、线性回归分析等不同模型对STN的空间分布进行探究[12-13],相关研究尺度涉及农田[14-15]、流域[16-17]、行政区[18]等,揭示了地形、成土母质、土壤类型、土地利用方式及人为管理等在不同程度上对STN空间分异的影响。空间自相关分析是探求变量在空间内分布模式的另一种方法,应用范围较广,可用于人口分布、城镇家庭消费以及土壤理化性质的空间变异研究[19-21],和地统计学相比,它可以有效识别出土壤理化性质的空间聚类特征,并对空间相关的显著性进行检验,在衡量自相关时更为严格。将两种方法结合使用,可以对土壤养分的空间变异特征进行更深入的剖析,而目前相关研究成果较少。

STN含量不仅对土壤肥力和耕地质量有着重大的影响,而且还潜在地威胁着水质安全。目前关于该区域的相关研究成果较少,为此,本研究以十堰市为案例,综合应用地统计学方法及空间自相关方法,对该区域耕层STN的空间变异特征及规律进行探讨,并进一步判断STN的空间聚集区和空间孤立区的具体分布位置,为研究区确定合理的STN管理单元、实现农业生产土壤氮肥精准管理、确保南水北调中线水源地生态安全和经济社会可持续发展提供参考。

1 材料与方法

1.1 研究区概况

十堰市位于湖北省西北部,秦巴山区汉水谷地,是鄂、豫、陕、渝交界区域性中心城市,地理坐标介于东经109°29′～111°16′和北纬31°30′～33°16′之间。该区域属亚热带季风气候,区内多年平均气温15.4 ℃,年平均日照1650 h,年平均无霜期224 d。多年平均降水量为769.6 mm,降水集中分布在夏季;降水量表现为东南部多西北部少。该区域地势南北高,中间低,自西南向东北倾斜;地貌类型以山地为主。土壤类型以黄棕壤、石灰土和水稻土为主。成土母质主要有泥质岩类风化物、第四纪老冲积物、红砂岩类冲积物、碳酸盐类风化物、石英质岩风化物、河湖冲沉积物,还有少量紫色岩类风化物及结晶盐类风化物。该区域耕地具有明显的地域性和立体分布规律,耕地面积随高度的升高而减少。耕地利用以水田、旱地和水浇地为主。

1.2 数据来源与数据预处理

数据来源于2017年秋十堰各县市区耕地质量调查与评价。

对数据处理采用拉依达准则法(3倍的标准差),利用SPSS 24.0对采样点的STN原始数据进行特异值的识别与剔除,剔除异常值后总样点数为666个,本文相关内容的研究均采用剔除异常值后的数据。

利用Minitab 17.1软件对STN含量进行数据变换处理,选取最优的正态分布拟合系数λ,使用SPSS 24.0进行描述性统计分析及相关分析;其中相关分析用于揭示不同因素与STN的相关性是否显著;在ArcGIS 10.7中对样点进行坐标转换和投影(Beijing_1954_3_Degree_GK_CM_111E),利用GS+9.0软件进行半方差分析、空间自相关分析及理论模型的拟合,普通克里格插值在ArcGIS 10.7中进行,通过GeoDa对STN空间自相关图进行绘制。

1.3 研究方法

1.3.1 地统计分析法

方差函数是地统计学的理论基础,在本研究中用来反映区域化变量SOM在研究区内的空间变异与相关程度。计算公式如下:

(1)

式中:r(h)为半方差函数;N(h)是以h为间距的所有观测点的成对数目。h为两个样本点的分隔距离;Z(xi)和Z(x+h)分别为区域化变量z(x)在xi和xi+h处的实测值。

1.3.2 空间自相关分析方法

空间自相关是描述事物间相互关系的方法。它通过事物之间在地理空间的相互关系来判断事物之间的相互关系(集聚或分散),在地理学上应用广泛。包括全局空间自相关和局部空间自相关,空间自相关分析主要用莫兰指数(Moran’sI)来反映空间要素在整个研究区域范围内的空间聚集程度,Moran’sI计算公式如下:

(2)

对于全局Moran’sI,一般用标准化统计量阈值Z(I)进行显著性检验,其一般公式为:

(3)

式(3)中Z(I)为标准化统计量;E(I)为期望值;Var(I)为理论方差。当|Z|>1.96时,表明该属性在空间上存在显著的空间自相关性,否则,不存在空间自相关。

全局自相关忽视了区域内部要素间的相互关系,因此用局部空间自相关(Local indicators of spatial association,LISA)的方法来分析每个空间要素属性与其临近空间要素属性的关联程度。Moran’sI的LISA计算公式为:

(4)

各变量含义与公式(3)中相同,全局空间自相关分析主要由Moran’sI反映采样区域STN在研究区空间分布情况,Moran’sI的范围在[-1,1]之间。若Moran’sI>0且通过显著性检验,则说明研究区域存在空间集聚,表示空间数据属性值高的样点相互邻近(记为H-H),或者是空间数据属性低的样点相互邻近(记为L-L);反之,当Moran’sI小于0,则说明研究区域存在空间负相关,表示空间数据属性值高的样点与空间数据属性低的样点相互邻近(记为H-L),或者是空间数据属性低的样点与空间数据属性值高的样点相互邻近(记为L-H)。

2 结果与分析

2.1 基本统计学特征

2.1.1 STN的描述性统计

根据第二次土壤普查中STN分级方法,研究区STN含量分为6级(表1)。表1显示,1到6级的样本数分别为7.70%,13.55%,28.53%,21.54%,16.83%和11.84%。研究区STN含量集中分布在中等、缺乏及较缺乏三个级别,占总数的63.62%。研究区STN含量范围为0.124～3.720 g/kg,均值是1.118 g/kg,表明研究区耕层STN处于中等水平,变异系数为49.374%,属于中等程度的空间变异性。

表1 STN值分级及各级分布

全部701个原始样点的STN数据剔除异常值,获取666个有效样点,剔除异常值后标准差从0.552 g/kg降低到0.453 g/kg,进行box-cox变换后,标准差降低到0.246 g/kg,可以初步判断出全氮分布离散程度有所降低;原始样本的偏度系数为0.957,峰度系数为1.025,说明STN数据右偏,剔除异常值后,偏度系数从0.957减少到0.446,进行box-cox变换后,偏度系数减少到0.033,接近0,全氮数据更接近正态分布。

2.1.2 STN与环境因子的相关分析

本研究中,运用SPSS 24.0对十堰市STN与地形因子(高程、坡度、水平曲率、剖面曲率)、气候因子(年均温、年降水量)、主要土壤因子(土壤pH值、有机质)、距交通道路距离、距河流距离进行相关分析,结果见表2。如表2显示,十堰市STN分别与高程及有机质呈显著正相关,与土壤pH值及年均温呈显著负相关,说明这四个因素是影响研究区STN空间分布的重要因子。其中,高程通过影响土壤中水热资源的再分配与土壤生态系统的物质循环过程与强度来影响土壤氮素分布;微生物是驱动土壤氮循环的引擎,而微生物的代谢活性受到温度的制约;土壤pH值主要是通过影响土壤微生物的活性进而影响不同区域土壤有机质及氮素水平;STN含量的高低取决于有机质的积累和分解。由于结构性因素和随机因素的影响,十堰市不同区域土壤pH值及有机质空间分布存在差异,进而对STN的空间分异产生影响。

表2 STN与环境因子的相关系数

2.2 STN空间变异分析

2.2.1 各向同性结构特征

运用GS+9.0软件对研究区STN的各向同性半方差函数进行拟合,拟合结果显示,指数模型的决定系数(R2)最大而残差(RSS)最小;根据决定系数(R2)最大,残差平方和(RSS)最小,并优先考虑RSS的原则[22],得到研究区STN各向同性条件下最优半方差函数模型(图1)及相关拟合参数(表3)。半方差理论模型为指数模型时,决定系数(R2)为0.972,残差(RSS)为3.908E-05,拟合程度较高。研究区STN半方差模型的块金值为0.0360,说明在当前采样尺度下仍然存在由试验误差、耕作及施肥等随机因素引起不同程度的变异。土壤养分的空间变异性受结构性因素和随机性因素的共同作用。研究区STN的块金系数值为0.348,属于中等程度的空间自相关,略大于0.25,表明十堰市STN空间变异受结构性要素影响较大,受随机要素影响较小。变程为321 600 m,远大于实际采样间距,说明各采样点均在空间变异范围内,具有代表性。

表3 STN半方差函数模型及参数

图2 STN空间分布图

图3 STN各向同性(a)和各向异性(b)条件下全局自相关系数

2.2.2 各向异性结构特征

研究区地形比较复杂,不同地貌类型下的土地利用差异较大。为此,本文对STN的各向异性进行深入分析,绘制0°(N-S)、45°(NE-WS)、90°(E-W)和135°(SE-WN)共4个方向的半方差函数图(图1)。各向异性分析表明:STN含量在间距小于20 km时,各方向变化整体上呈现各向同性;当间距在20～80 km时,各方向各向异性均呈增大趋势,其中0°(N-S)方向和45°(NE-WS)方向变异程度较大,90°(E-W)方向和135°(SE-WN)方向变异程度较小,但4个方向变异整体差异较小;当间距大于80 km时,90°(E-W)方向变异呈减小趋势,其他3个方向各向异性进一步增大,其中0°(N-S)方向变化最为剧烈。综合STN含量的各向异性分析,当间隔距离增大时,0°(N-S)、45°(NE-WS)和135°(SE-WN)方向上变异均表现出随间隔距离的增大而增大的趋势,其中0°(N-S)方向变化较为剧烈,90°(E-W)方向上变异则表现出随间隔距离的增大先增大后变小的确实,整体上变化较小,处于稳定状态。这可能与研究区的地形差异有关。由于0°(N-S)和45°(N-S)方向自南向北穿越大巴山、武当山、汉江谷底,秦岭及其他山间谷地,地势起伏较大,土地利用类型跨度大,而90°(E-W)和135°(SE-WN)方向上处于汉江、堵河及其支流的河谷地带,地势起伏较小,土地利用类型跨度较小。综上所述,在0°(N-S)、45°(NE-WS)和135°(SE-WN)方向上,变异程度随着采样粒度的增加逐渐增大,而90°(E-W)方向上随采样粒度的增大,变异程度先增大后减小,该方向上整体变化不大。

2.3 STN的空间分布

将GS+9.0软件模拟克里格插值最优参数输入ArcGis 10.7中,通过普通克里格插值,绘制STN的空间分布图(图2)。由图2可知,十堰市STN取值范围为0.12～2.96 g/kg。研究区高值区主要分布在竹溪县西南部、竹山县南部和房县南部;低值区主要分布在三个区域:郧西县中东部、竹山县中部和丹江口市内丹江口水库周边区域,整体上呈现从南向北逐渐递减的趋势。研究区STN含量空间特征与研究区地形变化具有较强的一致性:高值区分布海拔较高区域,低值区主要分布在汉江及其支流堵河河谷等海拔较低区域,STN含量呈现随高程降低而不断减小的趋势。其中0°(N-S)方向和45°(NE-WS)方向地势起伏比较大,局部气候、耕地利用、农作物种植类型及种植制度复杂,导致STN含量空间变异比较剧烈,呈高低相间分布,而90°(E-W)方向和135°(SE-WN)方向地形以河谷平原为主,地势及地表景观变化小,多为耕作区,土地利用方式相对简单,种植制度比较单一,STN含量空间变化微弱。从以上分析可以看出,研究区各向半方差函数模型表现出不同的特征,各向异性受地形及由其产生的地表环境差异的影响较大。

2.4 STN全局空间自相关特征分析

本研究用GeoDa对STN数据进行空间权重的计算和空间自相关性分析,结果表明,Z(I)值为18.21,表明研究区STN存在强烈的全局空间自相关性。

为进一步描述STN含量的空间自相关随样点距离的变化关系,根据其空间自相关分析结果,利用GS+9.0软件绘制STN各向同性(图3(a))及各向异性(图3(b))条件下空间自相关Moran’sI指数随间隔距离的变化图。由图3(a)可以看出,在各向同性下,STN的自相关性随着间隔距离的增加呈现较大变化。在0～80 000 m范围内,STN含量的Moran’sI指数均大于0,说明在结构上呈现正相关,存在空间集聚,随着间隔距离的增加,指数逐渐变小,说明其空间自相关性逐渐减弱,其中在60 000～80 000 m范围内,Moran’sI指数保持在0.00～0.05之间趋于恒定,自相关性较弱。在80 000～105 000 m范围内,STN含量的Moran’sI指数小于0,说明在结构上呈现负相关,存在空间孤立,随着间隔距离的增大,指数绝对值逐渐增大,说明其空间自相关性逐渐增强;在间隔距离>95 000 m后,STN空间自相关性表现为平稳状态。从地统计学角度来说,近距离的样点具有相同的表生作用、气候条件、地形条件及成土母质等,这些因素呈斑块状或连续性分布,呈现出较强的空间相关性,反映出的空间结构性较好,随着样点间距的增加,这些因素表现出较为明显的差异,另外受随机因素的影响,其空间结构性会受到影响[17]。

图2 STN空间分布图

图3 STN各向同性(a)和各向异性(b)条件下全局自相关系数

通过对空间自相关各向异性分析,可反映不同方向上STN的空间聚集和空间孤立特征。由图3(b)可以看出,在各向异性条件下,研究区STN在不同方向均表现出不同的空间自相关特性。在0°(N-S)方向上,STN含量的Moran’sI指数在间隔距离0～65 000 m范围内大于0,STN含量表现为正相关,在此范围内存在H-H和L-L空间聚集,随着距离的增加,正相关呈现逐渐减弱的趋势;当间隔距离在65 000～105 000 m范围内时,STN含量的Moran’sI指数小于0,且随着距离的增加,其绝对值越来越大,STN含量呈现H-L和L-H的空间孤立现象,在局部上变化比较明显;整体说明在该方向上>65 000 m表现为不同性质的斑块,从0°(N-S)方向上地形地貌、高程及成土母质均有所差异。在45°(NE-WS)方向上,Moran’sI指数在间隔距离为0～62 000 m的范围以正相关占主导,在此范围内存在H-H和L-L空间聚集;当间隔距离在62 000～105 000 m范围内,Moran’sI指数小于0,其绝对值随间隔距离增加逐渐增大,存在H-L孤立现象。在90°(E-W)方向上,Moran’sI指数在间隔距离0～10 500 m范围内均大于0,随着间隔距离的增加,Moran’sI指数呈现缓慢减小的趋势,在间隔距离>60 000 m时,Moran’sI指数在0.10附近趋于恒定,说明在此方向上不存在负相关,STN含量在此方向上存在H-H和L-L空间聚集;究其原因,可能受研究区地形地势影响较大,整个研究区地势南北高,中间低,汉江自西向东穿越全境。90°(E-W)方向上地形以河谷和山间盆地为主,成土母质和土壤类型走向受地形影响较大,该区域的土地利用与管理相对较单一,在结构因素和随机因素差异不大的情况下,STN的空间自相关程度就更加显著;在斑块性质相同的条件下,STN含量的空间结构性在此方向上得到了较好的反映。在135°(SE-WN)方向上,STN含量的Moran’sI指数在间隔距离0～105 000 m范围内以正相关为主,负相关不显著,主要呈现H-H和L-L的空间聚集特点;其成因与90°(E-W)方向大致相似,研究区内众多河流基本上呈现西北东南或东南西北流向,河谷面积较大,但受到中部武当山的影响,其土地利用和管理存在一定的变化,空间自相关没有90°(E-W)方向表现的显著,但结构性在此方向上也得到了较好的反映。总的来说,研究区内各向同性和各向异性自相关特征的特征,与区域内部不同尺度和不同方向上的地形地貌特征以及受其影响造成的土地利用方式及人为管理的差异性有着密切关系,其具体的影响过程有待进一步验证。

图4 STN局部Moran’s I系数散点图

2.5 STN局部空间自相关特征分析

为揭示研究区不同样点STN空间自相关分布格局,运用GeoDa软件对研究区耕层STN含量进行局部空间自相关分析,绘制了研究区STN的局部Moran’sI指数散点图(图4)。

结果显示,研究区STN局部Moran’sI指数为0.3576,说明不同土壤样点STN空间分布存在显著的空间集聚。H-H和L-L集中分布在第一、三象限,即STN高值区(或低值区)被STN高值区(或低值区)包围;第二、四象限为L-H和H-L的分布模式。即STN高值区(或低值区)被STN低值区(或高值区)包围。总体来说,十堰市耕层STN主要呈现H-H和L-L并存的格局,其中分布在第二、四象限样点较少,在高低值混合分布情形中,低值零星分布于高值区和高值零星分布于低值区的情况大致相同。

基于局部Moran’sI指数散点图,对局部Moran’sI指数进行显著性检验,并在研究区范围内进行可视化分析(图5)。图5揭示了在95%的置信水平下,研究区STN含量在空间分布上呈现的特性。STN局部Moran’sI指数共有482个样点未通过显著性检验(P=0.05),达到显著性水平的以L-L型最多,有82个点,主要分布在郧阳县和丹江口市;该类型分布区域主要在汉江河谷和丹江口水库分布区域,海拔较低,气温较高;其次是H-H型,共有69个点,主要分布在竹溪县和房县境内,位于秦岭山区,地势较高,气温较低。环境因子相关分析结果表明(表2),研究区气温与STN呈显著负相关,而降水与STN相关性并不显著,说明和降水相比,气温对STN的影响更大。一般而言,同一时期,气温较低有利于STN的积累[21]。由于研究区地形条件复杂,气候的垂直地带性得以体现,气温和降水等条件直接影响STN的迁移和转化,决定着母质分化和成土过程的方向和强度;同时影响着耕作制度和耕作方式,海拔较低区域耕种条件好,农作物产量高,另外,耕作制度和海拔较高处存在差别,作物从土壤中带走更多的土壤养分,因此海拔较高处,STN含量较高,海拔较低处则相反。

L-H和H-L的样点较少,分别为19个和14个;其中L-H主要分布在郧阳县和丹江口市,H-L主要分布在竹溪县和房县境内。L-H和H-L样点的分布体现了局部地区STN分布的空间异向性。除了地形、气候、母质等的影响外,研究区耕地地块破碎,当以小农经营为主,农户农田管理差异较大,从而导致STN含量出现L-H和H-L的空间结构。对于L-H型样点,周边被高值环绕,农户可以借鉴周边的农田管理措施,促使其向H-H型转化;而H-L型样点STN含量高于周边,其原因可做进一步分析,以期为周边区域培肥土壤提供参考。

通过图2和图5的对比可以看出,由普通克里格插值法预测的STN含量的空间分布情况与局部空间相关性分布情况结果基本一致,说明在此研究尺度下,十堰市STN在半变异函数计算范围和空间自相关距离一致时含量呈现显著的空间异质性和空间相关性。

图5 STN局部空间自相关特征分布图

3 结论

3.1 根据第二次土壤普查中STN分级标准,研究区STN含量集中分布在中等、缺乏及较缺乏三个级别,呈现中等程度的空间变异性。空间模型采用指数模型较好,研究区STN变程为321 600 m,步长为7050 m,块金系数为0.348,略大于0.25,说明区域内存在中等强度的空间自相关性,结构性因素对研究区STN空间变异起主导作用。说明土地利用方式及生产方式等随机因素的影响处于次要位置。各向异性变异结构中0°(N-S)方向变异程度占主导作用。

3.2 插值结果表明:高值区主要分布在竹溪县和房县南部等海拔较高区域;低值区主要分布在海拔较低区域:郧西县中南部、竹山县中部和丹江口市内丹江口水库周边区域及竹溪县和房县交界处,整体上呈现从南向北逐渐递减的空间分布格局。

3.3 研究区耕层STN的Moran’sI指数为0.3576,标准化Z(I)值为18.21, STN存在强烈的全局空间自相关性,耕层STN含量值趋于空间集聚。在全局空间自相关分析中,90°(E-W)方向上空间自相关性较好,STN含量的空间结构性在此方向上得到了较好的反映。

3.4 绘制局部Moran’sI指数散点图可直接显示耕层STN空间集聚与空间孤立的具体位置,其中耕层STN空间分布以L-L型和H-H型为主,其中L-L型主要分布在郧阳县和丹江口市,H-H型主要分布在竹溪县和房县境内。L-H和H-L的样点较少,其中L-H主要分布在郧阳县和丹江口市,H-L主要分布在竹溪县和房县境内。

4 讨论

耕层STN的空间分布特征受地形地貌、气候因子、成土母质等结构性因素及土地利用、施肥、耕种制度及人为管理等随机性因素的影响。其空间分异规律是成土母质及外源长期输入的结果。掌握STN的空间特征,划定不同的农业生产区域,为STN的精准管理及农业生产投入提供依据,是实现耕地土壤可持续利用的前提[23-24]。研究区地形地貌复杂,成土母质各异,同时该地为地形、气候过渡带,水热状况垂直差异明显,强烈影响着STN的空间分布。全思懋等[14]通过对江苏省农田STN含量变化及其影响因子进行研究发现,气候对STN含量变化的影响并不显著,与本研究结果存在一定的差异,这可能与该区域气候因子的地带性差异较小有关。刘合满等[25]对亚热带红壤丘陵典型区STN的空间变异特征进行探索,认为高程决定了区域的气温、湿度和植被条件等,一方面,不同高程常年气温和降水存在差异,导致不同高程区域土壤有机氮分解速率和雨水淋失产生差别;另一方面,高程对土地利用方式影响较大,作物种植类型、耕作制度、生产方式及人为管理的变化在一定程度上可以造成STN的变化。王勤等[26]探究了合肥老城区绿地土壤pH和氮磷的空间变异特征,发现土壤pH对土壤氮素转化有显著影响,较高的土壤pH值导致土壤硝化作用增强,可能会增加氮素损失。赵业婷等[27]通过对渭北台塬区耕地土壤有机质与全氮空间特征研究,认为氮素主要以有机质形态存在,STN含量的消长取决于各地区有机质的积累和分解作用的相对强度。本研究结合两种方法,能够更精确掌握研究区STN的空间分布格局,并能够准确识别STN的空间集聚和空间孤立区,从整体上发现耕地土壤养分退化和重建过程中存在的问题,并采取针对性措施进行保护与培育,因地制宜地挖掘区域最大的农业生产潜力,提高粮食综合生产能力。同时,还有助于对农业面源氮污染重点控制区进行识别,为划定农业生产中环境保护的重点区域提供依据。

本研究未对影响其空间分异的结构性因素和随机因素做深入分析,另外,不同施肥条件及大气沉降等因素对区域耕层STN空间异质性的影响,受数据获取所限,本文未做深入研究;由于篇幅所限,未能采用多种方法进行STN空间分布进行预测及对不同方法的预测结果进行精度评价。未来可以通过对研究区内部不同尺度不同时期采样数据进行时空分析,扩大数据搜集范围,采用多种方法对STN的时空变化规律、局域性格局动态变化特点及其成因进行探讨将是笔者下一步的工作。