石羊河下游民勤盆地土壤盐分空间分异特征

刘 新， 郝媛媛， 花立民

（甘肃农业大学草业学院，草业生态系统教育部重点实验室，国家林业草原高寒草地鼠害防控工程技术研究中心，甘肃 兰州 730070）

土壤盐渍化是自然环境和人类活动共同作用影响下的土壤退化现象，也是全球性的生态环境问题之一，不仅制约着农业经济的发展，也危害了生态环境的稳定性［1］。土壤盐分一般在表土层（0～30 cm）分布较为广泛［2］，其量化后可以表征盐渍化的进程，同时盐分在不同土层、方向的分布与移动导致了盐渍化的纷繁多变［3-4］。而在盐分运移过程中受母质、环境和生物等因素影响而具有空间异质性，并通过空间上斑块之间的差异来体现该性质，进而呈现出不同斑块之间的（自）相关性［5］。通过获取空间异质性与其自相关性有利于土壤资源的科学利用和明确盐渍化斑块对环境的影响。因此，通过量化土壤盐分来获取区域化的盐分数据，可以在不同区域中展示研究区内盐渍化的空间分布，并对土壤盐渍化动态监测具有重要意义。同时，细化表土层后，研究发现盐分受到结构因素和随机因素的综合影响，从而导致土壤盐分累积或散失量的差异［6-7］。通过对比不同土层之间的盐分变化可以掌握不同土层中盐分的差异，有助于更清楚地了解土壤盐渍化进程［8］。综上所述，结合土壤盐分量化和土层之间的空间分异特征可以更深入地理解盐渍化现象。

民勤盆地作为干旱区荒漠-绿洲生态系统的典型代表，风沙、盐渍和干旱使境内土壤盐渍化较严重［9］。现有研究多从地下水、地表水的水体矿化度出发［10］，到土壤盐分离子的变化［11］，以及不同植被群落的土壤盐分与养分的空间变化［12-13］，而缺乏对民勤盆地土壤盐分在空间上的系统分异研究。本文以土壤盐分量化为基础，辅以表土层不同层次之间的含量和区域差异，以此分析民勤盆地盐渍化程度分布变化与空间分异特征，以期为民勤盆地的土壤盐渍化变化提供基础数据，为生态恢复及土壤退化治理提供科学依据。

1 材料与方法

1.1 研究区概况

民勤盆地（38°24′～39°15′N,102°40′～103°55′E）位于甘肃省武威市民勤县，处于石羊河流域下游且由腾格里沙漠和巴丹吉林沙漠三面（除西面）包围（图1）。属温带大陆性干旱气候区，四季分明，干燥少雨多风，年均温7～8 ℃，风沙频发，年均风速2.7 m·s-1，降 水量110 mm 左右，年平均蒸发量2483 mm。作为流域的尾闾段，沙漠侵袭和咸水灌溉导致土壤沙化、盐化加剧。土壤类型主要是灰棕漠土、沙土和盐化土。植被以典型的沙生、旱生和盐生植物为主，白刺（Nitrariaspp.）、梭梭（Haloxylon ammodendron）和柽柳（Tamarix chinensis）是其主要的优势灌木。

图1 研究区地理位置及采样点分布Fig.1 Location and distribution of sampling points in the study area

1.2 土壤盐分数据采集及测定

选取人为因素影响较小且植被分布相对均匀的区域，每个采样点以南北方向间隔100 m 布设10 m×10 m 样方，3 个记为重复。样方内南北对角线取3 个点的土样（深度为0～10 cm、10～20 cm 和20～30 cm）并分层混匀，四分法装袋后带回实验室备用。2020 年（79 个样点）和2021 年（51 个样点）共计130个采样点，经自然风干、拣去杂物、研磨过筛后取标准土样用于盐分测定。利用样品电导率及与烘干残渣法联合获取盐分含量［14］。参考《中国盐渍土》和《土壤农化分析》盐渍化类型划分指标对研究区土壤盐渍化进行分类分级（<2 g·kg-1非盐化土；2～3 g·kg-1轻度盐渍土；3～5 g·kg-1中度盐渍土；5～10 g·kg-1强度盐渍土；>10 g·kg-1盐渍土）［15-17］，并评价盐渍化程度。

1.3 研究方法

1.3.1 正态性检验及方差分析 Kolmogorov-Smirnov检验（K-S 检验）是一种常用的拟合优度检验方法［18］，可根据样本数据推断其来自的总体是否服从某一特定理论分布，在5%置信水平下进行土壤盐分数据的正态分布检验，并根据P值大小得出结论［19］。单因素方差分析一般用于两个及以上样本平均数差别的显著性检验，用于研究某一个类型的自变量与另一个数值型的因变量之间的差异（5%）［11］。

1.3.2 空间自相关性 利用Moran’sI指数展现土壤盐分的空间依赖性（自相关性），计算公式为：

式中：n为需要估计的样点数目；xi和xj分别为变量x在样点i和j处的观测值；xˉ为x的平均值；Wij为相邻权重。值域为［-1,1］，当I=0 时，表示变量在此空间分布上无相关性，且呈随机分布；当I值趋于1时，变量之间的空间分布关系呈现完全正相关和聚集状态。反之，I值趋于-1 时，变量呈完全负相关且分散分布。

1.3.3 空间异质性 利用半变异函数的拟合参数解析土壤盐分分布的空间异质性特征，并从随机性和结构性变异程度两个方面进行解释。选取块金值（Nugget，C0）、基台值（Sill，Ci+C）、偏基台值（Partial Sill，C）和变程（Range，A0）4 个主要参数。C0+C表示变量的总变异性，其中C0表示随机因素引起的空间异质性，C表示结构（系统）因素引起的；二者之比C0/(C0+C），即基底效应展现由空间自相关的部分引发的空间变异程度；A0表示变量具有空间相关性的阈值大小，即在阈值范围内该变量具有空间相关性，超出范围则独立无相关性［9,20］。

1.3.4 Kriging插值 Kriging插值是地统计学中评估空间分布的一种方法，包含多种算法，常用的有普通Kriging、简单Kriging 和泛Kriging［21］。本研究采用普通Kriging法对土壤盐分进行插值。公式［22］为：

式中：S0为待估算盐分含量的样点；Si为已知样点i的盐分含量；Wi为样点i的权重值；x为参与计算的已知盐分含量样点的数量。

1.3.5 模型精度检验 为了保证预测值的无偏性，采用平均误差（Mean Error，ME）、平均标准误差（Mean Standardized Errors，MSE）、标准均方根预测误差（Root-Mean-Square Standardized Errors，RMSSE）、均方根预测误差（Root-Mean-Square Prediction Error，RMSE）、平均标准误差（Average Standard Error，ASE）和决定系数（R2）6 个指标对Kriging 插值结果进行精度检验。

1.4 数据分析

运用GS+9.0 进行半变异函数模型筛选和Moran’sI指数分析，Kriging 插值以及空间分布图使用ArcGIS 10.3 建模和绘制，利用Excel 2016 进行统计分析。

2 结果与分析

2.1 统计分析及正态性检验

研究区内土壤表层土上层（0～10 cm）、中层（10～20 cm）和下层（20～30 cm）的盐分含量呈逐渐递减的趋势，即随土层变深盐分含量均值逐渐减少（5.44 g·kg-1→5.23 g·kg-1→5.07 g·kg-1），且盐分变化范 围 缩 小（1.19～32.08 g·kg-1→1.24～19.05 g·kg-1→1.43～16.29 g·kg-1）。同时，3个土层中土壤盐分含量的中位数均小于其均值，且标准差为3.50±0.38 g·kg-1，说明50%以上的样点盐分含量在均值水平以下且离散程度较大。K-S 正态性检验（95%置信区间）结果表明：0～10 cm、10～20 cm 和20～30 cm 3 个土层土壤盐分经对数变换后的P值分别为0.2902、0.7952 和0.5170，均＞0.05，呈偏态（偏度0.2076～0.8057）正态分布（图2）。方差分析结果表明，不同土层盐分无显著差异（P=0.164），因此可用地统计学法探索其空间分异特征。

图2 样点盐分含量正态分布Fig.2 Normal distribution of salt content in samples

2.2 空间自相关性

3 个土层土壤盐分的Moran’sI指数均呈“∽”型波动趋势（图3），即随着距离的增加，Moran’sI指数在0 值上下波动，且正空间自相关性大于负空间自相关性；值域分别为［-0.151, 0.302］、［-0.185,0.392］和［-0.259, 0.552］，范围随土层加深而变大，空间自相关性和空间差异性均随土层加深而变大；除0～10 cm 负的空间自相关性（随距离增大而增大）外，其他土层的正（负）空间自相关性均随距离的增大而减小。

图3 土壤各层盐分Moran’s I指数Fig.3 Moran’s I index of soil salinity in each layer

2.3 空间异质性

通过对线性、球状、指数和高斯4种变异函数模型模拟的土壤盐分结果进行反复比较，筛选出最优半变异函数模型及其插值相关参数（表1）。3 个土层均以指数模型为最优（R2随土层加深而变大（0.62→0.69→0.81）），且RSS（残差）均是0.0002，说明3个土层的模拟结果均可以有效表征各土层土壤盐分的空间异质性。C0（块金值）在0.008～0.060，是C的13.80%～66.67%，说明在模拟过程中随机因素得到了很好的控制，较为稳定而不影响模拟结果。C0/(C0+C）（基底效应）随土层加深由弱到强烈再到中等，说明盐分含量在不同土层中的空间异质性受随机性和结构性因素的影响不同。A0（变程）的变化与C0/(C0+C）一致。其中，0～10 cm 土层因受雨雪、风沙、践踏、垦殖等随机性因素的影响较大，而结构性因素的影响相对较弱，以致盐分含量变化较其余两层剧烈，空间异质性最弱且A0最大。

表1 土壤各层盐分的变异函数理论模型及相关参数Tab.1 Variogram theoretical model and related parameters of soil salinity in each layer

2.4 空间分布格局

2.4.1 插值及精度验证 对插值结果（图4）进行精度验证（表2）可知，3个土层ME（平均误差）（-0.004～0.006）与MSE（平均标准误差）（-0.018～0.017）均趋近于0，筛选的变异模型满足对土壤盐分预测的准确性与有效性；ASE（平均标准误差）与RMSE（均方根预测误差）差值的绝对值均≤0.03，对变异程度的预测控制较稳定可信；RMSSE（标准均方根预测误差）均趋近于1，对变异程度的预测在合理范围内。因此，参数设置合理，插值精度均符合要求，且10～20 cm＞20～30 cm＞0～10 cm。

表2 不同土层土壤盐分插值模型验证结果Tab.2 Verification results of interpolation model of soil salinity in different layers

图4 不同土层的土壤盐分空间分布Fig.4 Spatial distribution of soil salinity in different soil layers

2.4.2 空间分布特征 3 个土层土壤盐分在空间分布上存在一定差异性（图4），非盐化土和4种不同程度盐渍土在3 个土层中均有分布。水平方向上，盐渍化程度均自西南向东北逐渐增加，与石羊河流向一致；垂直方向的分布规律由简单到复杂，从聚集到分散。研究区98%以上区域为不同程度的盐渍土（表3），以中度盐渍土为主，面积均> 70.00%且20～30 cm >0～10 cm >10～20 cm；强度盐渍土次之，主要分布在腾格里沙漠边缘（研究区东北-北-南沿线边缘）区域，面积随土层加深而减小（由0～10 cm的23.43%减少到20～30 cm 的15.49%）且逐渐趋于斑块化；轻度盐渍土位列第三，主要分布于研究区西南毗邻武威盆地区域，随着土层加深研究区内部也逐渐出现斑块化分布区域，面积变化规律与强度盐渍土恰好相反（3.12%→12.20%→14.15%）；盐土面积最小（<0.40%）。非盐渍土主要分布在研究区西南部的盆地边缘，20～30 cm 土层占1.27%，面积最大；0～10 cm 土层次之，为20～30 cm 土层的1/2；10～20 cm土层最小，仅占0.03%。

表3 不同土层的土壤盐分面积占比Tab.3 Proportion of soil salt area in different soil layers at different times

3 讨论

3.1 采样点质量与盐分研究精度的关系

盐土的盐分组成与母岩的类型和成分有密切的联系，为了获取更具有代表性、典型性的样品，根据研究区的盐分状况、植被类型等确定优质的采样点［23］。采样点的确定需要结合实际，综合运用多种评价准则来分析获取。从采样点数量来看，以黄河三角洲为例，通过研究土壤盐分空间变异的变程指标确定出每1000 km2需采集107 个样点［24］，而不同地点有着不同的样点分布数量标准；从采样点属性来看，不同的间距与协变量（生境、地表覆盖、利用方式等）的选择对区域盐分含量估算有影响［25］；从处理方式来看，干燥与否、过筛大小、测定方法选取等不同的处理都会对盐分含量拟合模型精度有影响［23］；从结果获取来看，拟合（数学归纳、统计建模等）与反演（卫星、雷达、无人机以及影像的光谱选取）都会给结果带来误差［26］。因此，在样品选取与采集过程中需要综合考量方法、协变量、处理、结果展现等的影响，才能使得研究更加准确并符合实际。

3.2 沙化对盐渍化程度空间异质性的影响

不同土层不同斑块表现出空间异质性［27］。盐渍化程度分布基本呈现灌区向四周加深的趋势，恰与两大沙漠侵袭民勤盆地的方向相反，土壤沙化程度可能对盐渍化程度有影响［28］。插值结果表明，靠近腾格里沙漠区域的盐渍化程度高于靠近巴丹吉林沙漠区域（图4），实测原位数据结果亦是如此（表4），原因在于巴丹吉林沙漠的侵袭程度要高于腾格里沙漠［29-30］，在基于出入沙地、人员配置等因素考量设置腾格里沙漠边缘区域的采样点较多且分布广泛，而布设在巴丹吉林沙漠边缘区域的采样点相对较少且集中分布在沙漠边缘两端（图5），也会影响两大沙漠边缘区域土壤盐分的空间异质性。其潜在分布与实际情况之间复杂的变化与运移方式有待进一步研究［31］。

表4 沙漠边缘区域统计信息对比Tab.4 Comparison of statistical information of desert edge area

图5 沙漠边缘采样区域的对比Fig.5 Comparison of sampling areas at desert edge

3.3 民勤盆地盐渍化程度的变化趋势探究

民勤盆地是基于石羊河流域的农业垦殖区［32］。20世纪70年代以前，自然条件下的蒸发积盐是该地区土壤盐渍化的主要原因［33-34］；70年代之后，大量地下水开采灌溉引起的土壤盐分人为再分配加剧了该地区的次生盐渍化［35］。陈丽娟等［34］和顾梦鹤等［35］的研究发现，青土湖地区在石羊河流域综合治理生态恢复工程［36］中，初期（2010年）与中期（2016年）以强度盐渍土为主［34-35］，而本研究发现，工程结束后（2020—2021年）中度盐渍已达研究区的73.55%，说明近15 a来，研究区盐渍化程度有所改善，由强度逐渐向中度转换，间接表明了石羊河流域综合治理生态恢复工程对民勤盆地土壤盐渍化的恢复也具有一定效果。

3.4 空间异质性变化的探究

诸多学者对表层土壤的盐渍化程度变化提出了不同的思考。不同环境下，生物与非生物因素的互作关系是影响盐分分布或聚集多寡的关键因素［37-38］。民勤地区因所处环境限制，植被所延伸的丰富度［39］和生物量［40］等因素导致了不同土层变化程度的差异性；地下水［41］和土壤湿度［42］等水分条件又会从非生物的角度抑制或促进这些变化。不同的土层划分标准（间隔不同）也会对土壤盐分空间异质性的表达造成影响，综合多个土层的土壤剖面研究［43］与仅针对表土层甚至其他单一土层的结果不同且存在一定差距［44］。误差与样本量对于试验结果的影响是造成空间异质性差异的又一影响因素，包括方法选择、人为误差、数据累积和迭代程度以及系统变化等［45］。

4 结论

通过筛选最优变异模型获取研究区土壤盐分的空间分布格局，表现不同土层盐分的空间差异，结论如下：

（1）3 个土层土壤盐分含量和变化范围均随土层加深而变小，且均呈偏正态分布；Moran’sI指数均呈“∽”型波动趋势，值域范围随土层加深而变大，空间自相关性均随距离的增大而减小。

（2）土壤盐分的变异函数均以指数模型最佳，模拟结果均可有效表征各土层土壤盐分的空间异质性且随机性因素控制较为稳定；空间异质性随土层加深由弱到强烈再到中等，0～10 cm 土层因受雨雪、风沙、践踏、垦殖等随机性因素的影响较大，空间异质性最弱且变程最大。

（3）经交叉验证插值精度（10～20 cm＞20～30 cm＞0～10 cm）均符合要求。插值结果表明，3 个土层土壤盐分均存在空间分异，水平方向（同一土层）上，盐渍化程度均自西南向东北逐渐增加，与石羊河走向一致；垂直方向（不同土层）上，分布规律由简单到复杂，从聚集状态（0～10 cm）逐渐趋于斑块化（10～20 cm）和斑点化（20～30 cm）。

（4）非盐化土和4 种不同程度盐渍土在0～30 cm 土层中均有分布，以中度盐渍土为主（面积均＞70.00%）；强度盐渍土次之，主要分布在腾格里沙漠边缘（研究区东北-北-南沿线边缘）区域，面积随土层加深而减小。