基于指示Kriging的土壤盐渍化风险与地下水环境分析

李仙岳 崔佳琪 史海滨 孙亚楠 马红雨 贾 飚

(1.内蒙古农业大学水利与土木建筑工程学院， 呼和浩特 010018；2.内蒙古自治区水利科学研究院， 呼和浩特 010020)

0 引言

土壤盐渍化是全球化的生态环境问题，是制约农业可持续发展的主要障碍[1-3]。据统计，全球盐渍化土壤面积达9.5亿hm2，占陆地面积的7.26%。我国盐渍化土壤分布广泛，面积约为0.37亿hm2，占全国可利用土地面积的4.88%，其中耕地盐渍化面积达到0.09亿hm2，占全国耕地面积的6.62%，主要集中在西北、华北、东北及沿海地区，其中陕、甘、宁、青、蒙、新的盐渍土面积占全国盐渍化土壤面积的69.03%[4]。内蒙古河套灌区是西北地区典型的盐渍化灌区，轻度、中度和重度盐渍化耕地分别占总面积的29.8%、17.2%和9.2%[5]，因此，土壤盐渍化防治已成为该区农业和生态可持续发展的主要工作。地下水是盐分传输、累积和排泄的主要载体，地下水矿化度和地下水埋深是影响土壤盐渍化程度的最主要因素[6-9]，直接影响土壤盐分的变化，即盐随水移是土壤中盐分迁移的主要方式，水分的运移带动了盐分在土壤中的移动和累积，当地下水埋深小于临界深度或地下水矿化度超出一定值时，地下水中的盐分随上升毛管水不断迁移到作物根层和地表，此时土壤发生盐渍化的风险增大[10-12]。因此，研究土壤盐渍化与地下水环境的关系、控制合理的地下水埋深是防治土壤盐渍化的重要途径[13]。

目前，关于土壤盐渍化与地下水环境关系的研究较多，如常用土壤毛管水上升高度法[14]和野外调查统计法[15]研究特定观测点或典型小区发生土壤盐渍化时的地下水临界深度。为研究大区域范围内土壤盐渍化和地下水环境的关系，WANG等[16]应用经典统计学、地质统计学和GIS分析了北疆内陆河流域绿洲土壤盐渍化与地下水埋深和土地利用的关系，发现灌溉区土壤盐渍化累计面积大于非灌溉区，农业灌溉用水导致浅埋地下水上升，从而间接导致土壤盐渍化；周在明等[6]基于指示Kriging法对环渤海低平原区表层土壤、地下水埋深和地下水矿化度进行了空间变异分析，发现地下水矿化度大于2 g/L、地下水埋深小于3 m及土壤含盐量大于1 g/kg时的空间概率分布呈一致性。上述研究主要针对新疆内陆河流域和环渤海低平原区，而河套灌区属于沿黄流域，地下水环境受黄河来水和配水的影响非常大。徐英等[17]基于指示Kriging法分析了河套灌区解放闸灌域土壤盐渍化与地下水埋深之间的关系，发现区域尺度上土壤表层发生中度和轻度盐渍化时地下水位临界埋深分别为2.0 m和2.5 m，但仅考虑了地下水埋深，并未考虑地下水矿化度。以上研究均基于一个时期进行讨论，且仅定性比较二者的概率分布，而河套灌区由于受到气候和灌溉等因素的影响，春灌前和生育期发生土壤盐渍化时的地下水临界埋深和矿化度差异明显。故考虑灌水和非灌水期以及地下水埋深和矿化度等多因素过程，探索河套灌区地下水环境与土壤盐渍化的关系，对进一步揭示不同时期河套灌区土壤盐渍化与地下水环境的定量关系具有重要意义。

本文以河套灌区永济灌域为研究区，运用指示Kriging法分析春灌前和生育期的土壤表层盐分与地下水埋深、矿化度的空间概率分布规律，评价二者的空间分布关系并计算匹配率，确定不同时期灌域尺度上满足土壤盐渍化防治的地下水埋深和矿化度临界值，以期为盐渍化灌区土壤盐渍化的防治提供参考。

1 材料与方法

1.1 研究区概况

永济灌域(40°36′～41°13′N，107°13′～107°42′E)位于内蒙古河套灌区中游，现灌溉面积约1.122×105hm2，属温带大陆性气候，干旱少雨，日照时数长，蒸发强烈，年均气温8.9℃，年均降雨量135.2 mm，且多发生于生育期(4—9月)，年均蒸发量2 327.3 mm，属于典型的盐渍化灌域。灌域内耕地多以砂壤土和壤砂土为主，灌溉作物包括粮食作物和经济作物，粮食作物以夏玉米和春小麦为主，经济作物包括向日葵、蔬菜、瓜果等，是典型的灌溉农业区。灌域以引黄河水灌溉为主，每年约在4月中下旬开灌，以夏灌(4—9月)和秋浇(10—11月)为主，到11月中下旬秋浇结束，排泄方式以强烈的蒸腾蒸发为主，土壤盐渍化成为制约灌域农业发展的主要因素。随着节水改造工程的实施，灌域地下水埋深呈下降趋势，地下水呈矿化趋势，多年(1998—2017年)年均地下水埋深和矿化度分别为2.14 m和3.35 g/L，截至到2011年，永济干渠、永刚分干渠、隆盛示范区东济支渠和隆盛示范区西济支渠共衬砌53.35 km，干渠和分干渠共整治9.72 km，干沟、分干沟和支沟共整治608.758 km，建筑物配套改造127座，节水效果明显，引水量由1998—1999年的10.10亿m3降至2000—2017年的8.65亿m3，用水量减少，灌溉效率得到极大的提高。

1.2 样品采集与处理

在研究区共布置51眼地下水观测井，每眼井周围均布置盐分采样点，土壤盐分采样点同时考虑覆盖不同盐渍化程度的土壤，土壤盐分采样点和地下水观测井布置如图1所示，共有68个土壤盐分采样点。考虑到地下水对土壤返盐的滞后效应，于2019年4月中下旬(春灌前)和7月末(生育期)采集0～30 cm范围内土壤(每个采样点3个重复)，在实验室内自然风干、研磨、过2 mm筛后，按1∶5土水质量比配制浸提液，测定土壤含盐量[18]；考虑到地下水对土壤返盐影响具有一定滞后效应[17]，于2019年3月末(春灌前)和7月初(生育期)用皮尺与测绳测定(每个观测井重复观测3次)地下水埋深并采集水样，应用电导率法计算地下水矿化度，用以探讨春灌前和生育期土壤表层(0～30 cm)含盐量与地下水环境(地下水埋深和地下水矿化度)之间的关系。

1.3 数据处理方法

指示 Kriging法(Indicator Kriging)是一种非参数估计方法[19-23]，不仅可以削弱有偏分布，还可以抑制特异值对变异函数稳健性的影响，可用于估计满足给定阈值的指示变量条件概率，绘制相应的概率空间分布图，或称为风险分布图。本文采用地统计学软件 GS+7.0 确定不同阈值下土壤表层含盐量、地下水埋深和地下水矿化度的变异函数模型，并将其模型参数输入 ArcGIS 10.4 进行指示Kriging插值估计，绘制不同阈值条件下土壤表层盐分、地下水埋深和地下水矿化度的概率分布图，并应用GIS中的重叠面积提取模块计算各指标概率区面积比例构成及不同阈值条件下土壤表层含盐量与地下水埋深、矿化度概率空间分布图的匹配率。

2 结果与分析

2.1 统计特征值分析

由表1可知，研究区春灌前和生育期土壤表层含盐量均值分别为3.95、3.62 g/kg，变异系数分别为133%和158%，属于中度盐渍土，地下水埋深均值分别为2.75、2.22 m，变异系数分别为37%和56%，地下水矿化度均值分别为3.43、4.00 g/L，变异系数分别为110%和115%，属于微咸水。从不同指标来看，地下水埋深属于中等变异性，地下水矿化度和土壤表层含盐量属于强变异性，可能由于研究区范围较大，地势差异、土地耕种方式差异和人为排灌等因素的共同作用，导致土壤表层含盐量和地下水矿化度变异性较大。从不同时期来看，受耕作措施、灌溉制度和降雨等人为因素和气象因素的影响，生育期各指标变异性比春灌前较大。

表1 土壤表层含盐量、地下水埋深和地下水矿化度统计特征值

2.2 指示变异函数模型构建

根据轻度(含盐量2.0 g/kg)和中度(含盐量3.0 g/kg)盐渍化土壤的划分标准确定土壤表层含盐量阈值。结合类似灌域相关研究[12,14]和永济灌域实际地下水埋深情况，选取1.8、2.2、2.6、3.0 m作为地下水埋深的阈值。根据地下水矿化度的划分标准[24-26]，分别选取2.0、2.5、3.0 g/L作为地下水矿化度的阈值。本文中，当土壤表层含盐量大于等于阈值(2.0、3.0 g/kg)和地下水矿化度大于等于阈值(2.0、2.5、3.0 g/L)时，指示变换值为 1，否则为 0，当地下水埋深小于等于阈值(1.8、2.2、2.6、3.0 m)时，指示变换值为1，否则为0。

对指示变换后的土壤表层含盐量、地下水埋深和矿化度进行半变异函数分析，由表2可知，各变量指示变异函数均可用球状模型拟合，变程在6.47～15.51 km之间；不同阈值下各区域变量块金值较小，说明在本研究尺度上由采样误差等引起的变异不大[6]。块金值与基台值百分比均在25%～75%之间，呈中等程度的空间自相关性，即其变异是结构性因素(气候、地形、水文地质条件等非人为因素)和随机性因素(灌溉制度、地下水利用等人为活动)共同作用的结果。将变异函数模型参数输入GIS 10.4进行指示Kriging插值，绘制各变量满足相应阈值的概率空间分布图并计算不同概率的面积占比及不同变量之间有相同概率的面积匹配率，颜色由浅色到深色表示大于等于(小于等于)阈值的概率越大(图2、3、4)，以概率50%为分界线，大于50%为高概率区，反之为低概率区。

表2 研究区土壤表层含盐量、地下水埋深和矿化度指示变异函数理论模型

2.3 土壤表层盐分与地下水环境的概率空间分布特征

图2a和图2b为春灌前土壤易发生轻度盐渍化以上(含盐量大于等于2.0 g/kg)和中度盐渍化以上(含盐量大于等于3.0 g/kg)的空间概率分布图，结合表3可知，研究区49%的区域为土壤发生轻度盐渍化的高风险区(概率在50%以上)且主要分布在研究区的东南和西北部，对比图1可知，高风险区多分布在排水干沟和总排水干渠附近；23%的区域为发生中度盐渍化的高风险区，且均包含在轻度盐渍化高风险区范围内；生育期较春灌前(图3c、3d)，土壤发生轻度盐渍化和中度盐渍化的高风险区缩至32%和12%，缩小了17、11个百分点，研究区西北和东南部仍是高风险区集中地带。图3a～3d为春灌前不同阈值地下水埋深的概率分布，图中颜色越深，表明该处地下水埋深小于等于该阈值的概率越大，当阈值为1.8 m时，仅21%的区域为高概率区且主要分布在北部和东南端部；当阈值为2.2 m时，高概率区在原有基础上扩大，35%的区域为高概率区；当阈值为2.6 m时，西部地区由低概率区过渡为高概率区，高概率区面积比例增至63%；当阈值为3.0 m时，仅24%的区域为低概率区且主要集中分布于北端部、南部及西部边缘地区。生育期不同阈值下的地下水埋深高概率区呈扩大趋势(图3e～3h)，较春灌前分别扩至44%、59%、79%和88%，其中阈值为1.8 m和2.2 m时，面积扩大了23、24个百分点，阈值为2.6 m和3.0 m时扩大了16、12个百分点，当阈值为3.0 m时，整个研究区中北部均为高概率区，小阈值的高概率区包含在大阈值的高概率区，其分布区域随阈值的增大而逐渐扩大。图4a～4c为春灌前不同阈值下地下水矿化度的概率分布，图中颜色越深，表明该处地下水矿化度大于等于该阈值的概率越大，地下水矿化度阈值为2.0、2.5、3.0 g/L时高概率区面积比例分别为47%、35%和17%，且中部和北部地区是高概率集中分布区，低概率区主要集中于研究区中偏北和西南腹部，大阈值的高概率区包含在小阈值的高概率区，其分布区域随阈值的增大逐渐缩小；而生育期，不同阈值下的地下水矿化度高概率区分别扩至58%、43%和30%，较春灌前分别扩大了11、8、13个百分点，低概率区仍主要集中于中偏北和西南腹部。

表3 研究区不同阈值条件下各指标概率区面积比例构成

春灌前与生育期3个指标的空间概率分布变化主要由于春季返盐激烈，土壤表层盐分积聚量大[27]，而生育期多次灌水(夏灌到秋灌间一般灌6次水)使得土壤中部分盐分淋洗到地下水中[28]，地下水埋深变浅，同时，强烈的蒸发使地下水各离子组分浓缩，溶质浓度变大[29]，所以进入7月(生育期内的灌溉后期)后，浅埋地下水高概率区扩大，土壤盐渍化高风险区缩小，地下水矿化高风险区扩大。

2.4 土壤表层盐分与地下水环境的概率空间分布关系

对比图2和图3可知，春灌前地下水埋深小于等于2.6 m时的概率空间分布与土壤表层含盐量大于等于2.0 g/kg的概率空间分布有较大相似性，地下水埋深小于等于2.6 m的高概率区几乎覆盖轻度盐渍化的高风险区，且二者匹配率最高，为69%(图5a)；地下水埋深小于等于2.2 m时的概率空间分布与土壤表层含盐量大于等于3.0 g/kg的概率空间分布具有较高的相似度，地下水埋深小于等于2.2 m的高概率区基本覆盖中度盐渍化的高风险区，二者匹配率高达74%。生育期，地下水埋深小于等于2.2 m时的概率空间分布与土壤表层含盐量大于等于2.0 g/kg的概率空间分布匹配率最高，为62%；地下水埋深小于等于1.8 m时的概率空间分布与土壤表层含盐量大于等于3.0 g/kg的概率空间分布具有较高的相似度，匹配率为63%。因此，可判定春灌前土壤表层发生轻度盐渍化的地下水埋深临界值为2.6 m，发生中度盐渍化的地下水埋深临界值为2.2 m；生育期，土壤表层发生轻度盐渍化的地下水埋深临界值为2.2 m，发生中度盐渍化的地下水埋深临界值为1.8 m。

对比图2和图4可知，春灌前地下水矿化度大于等于2.0 g/L时的概率空间分布，与土壤表层含盐量大于等于2.0 g/kg(轻度盐渍化)的概率空间分布匹配率最高，为61%(图5b)；地下水矿化度大于等于2.5 g/L时的概率空间分布与土壤表层含盐量大于等于3.0 g/kg(中度盐渍化)的概率空间分布具有较高的相似度，二者匹配率为60%。而生育期，地下水矿化度大于等于2.5 g/L时的概率空间分布与土壤表层含盐量大于等于2.0 g/kg的概率空间分布匹配率最高，为59%；地下水矿化度大于等于3.0 g/L时的概率空间分布与土壤表层含盐量大于等于3.0 g/kg的概率空间分布具有较高的相似度，地下水矿化度大于等于3.0 g/L的高概率区基本覆盖中度盐渍化的高风险区，二者匹配率高达64%。可见，春灌前土壤发生轻度盐渍化时的临界地下水矿化度为2.0 g/L，发生中度盐渍化时的临界地下水矿化度为2.5 g/L；生育期土壤发生轻度盐渍化时的临界地下水矿化度为2.5 g/L，发生中度盐渍化时的临界地下水矿化度为3.0 g/L。

对比春灌前和生育期各指标之间的匹配率可知，春灌前各指标匹配率优于生育期，因为生育期土壤盐渍化受人为因素(种植结构、灌溉制度、耕作措施等)影响更大，导致生育期地下水环境对土壤盐渍化的影响不如春灌前明显；同时，土壤表层含盐量与地下水埋深的概率空间分布匹配率明显优于土壤表层含盐量与地下水矿化度的匹配率，地下水埋深对于土壤盐渍化的影响明显大于地下水矿化度，其原因为地下水中的盐分主要通过上升毛管水不断迁移到作物根层和地表，因此，埋深较浅时，地下水中的盐分更易导致土壤盐渍化。

3 讨论

当地下水埋深较浅时，地下水中的盐分易于随上升毛管水迁移到作物根层，蒸发后盐分表聚，易于产生次生盐渍化，特别是当地下水中的矿化度较大时，越易发生次生盐渍化。通过河套灌区地下水与土壤盐渍化关系的研究表明(表4)，河套灌区地下水埋深和土壤盐渍化密切相关，临界地下水埋深范围在1.5～2.7 m之间。目前更多的研究主要通过小区试验建立土壤和地下水临界埋深的关系，孔繁瑞等[30]和张义强等[31]利用测坑法和田间试验得出生育期适宜作物(葵花和玉米)生长的地下水埋深，分别为2.0 m左右和2.0～2.5 m，本研究在生育期防治轻度和中度盐渍化土壤时的临界地下水埋深分别为1.8 m和2.2 m，结果较相近；杨会峰等[32]通过数值模拟对全年不同时期的临界地下水埋深进行了研究，在春灌前临界地下水埋深较深，生育期较浅，总体趋势与本研究一致，且本研究的临界地下水埋深数值在该研究范围内。这些小区尺度的研究基本没有考虑不同区域土壤质地和灌溉、排水、降雨的差异，故从区域尺度进行探索防治土壤盐渍化的临界地下水埋深更具有实际意义。徐英等[23]在区域尺度上通过指示Kriging法得出解放闸灌域防治轻度和中度盐渍化土壤的临界地下水埋深分别为2.0 m和2.5 m，与本文研究的结果相似，但略低0.1～0.2 m，这主要由于解放闸灌域轻度以上盐渍化土占耕地的面积小于永济灌域[5]，因此，土壤更易发生盐渍化，另外该研究没有考虑地下水矿化度的影响。故总体上本研究得到的临界地下水埋深与目前相关研究结论相似，且把地下水埋深和矿化度进行综合考虑，可见应用指示Kriging法研究灌域或区域尺度上地下水临界埋深和矿化度是可行的。

表4 河套灌区防治土壤盐渍化/适宜作物生长的临界地下水埋深和矿化度

春灌前较生育期发生土壤盐渍化时的临界埋深大[32]，临界矿化度小，主要由于春灌前土壤处于融冻期，返盐激烈且不能及时脱盐[33]，盐分表聚加重，土壤表层含盐量大，更易发生土壤盐渍化，因此当地下水埋深小于等于2.2 m或地下水矿化度大于等于2.5 g/L时，土壤易发生中度盐渍化。而生育期，多次引黄灌溉的同时，降雨也多集中于此时期[29]，盐随水走，土壤表层盐分部分被淋洗至地下水中，土壤表层盐分变少，便不易发生土壤盐渍化，即发生盐渍化时的埋深临界值减小，矿化度临界值增大。同时，本文中春灌前土壤盐渍化与地下水环境概率分布的匹配率大于生育期，主要由于研究区春灌前降雨稀少，蒸发量也较小，无引黄水量，地下水受到较少的外界因素(气象因素和人为因素)干扰，而降雨、引水和强烈的蒸腾蒸发多集中于生育期[28]，这个时期种植结构及一系列耕作措施也会使土壤盐分变异性增大，因此，春灌前地下水环境对土壤盐渍化的影响更强烈。

4 结论

(1)河套灌区永济灌域春灌前和生育期地下水埋深的变异系数分别为37%和56%，属中等变异性，土壤表层含盐量和地下水矿化度变异系数分别为133%、158%和110%、115%，属强变异性，且生育期各指标变异系数相比春灌前较大；土壤表层含盐量、地下水埋深和地下水矿化度的指示变异函数理论模型均为球状模型，且不同阈值条件下，各变量呈中等强度的空间自相关性。

(2)春灌前，永济灌域49%的区域为土壤盐渍化(轻度盐渍化)的高风险区，35%的区域为浅埋地下水(埋深小于等于2.2 m)高概率区，17%的区域为地下水矿化(矿化度大于等于3.0 g/L)的高风险区；生育期，32%的区域为土壤盐渍化(轻度盐渍化)的高风险区，59%的区域为浅埋地下水(埋深小于等于2.2 m)高概率区，30%的区域为地下水矿化(矿化度大于等于3.0 g/L)的高风险区。春灌前较生育期土壤表层盐渍化高风险区扩大，浅埋地下水高概率区缩小，地下水矿化高风险区缩小。

(3)春灌前，永济灌域土壤表层发生轻度、中度盐渍化时的地下水临界埋深分别为2.6、2.2 m，发生轻度、中度盐渍化时的地下水临界矿化度分别为2.0、2.5 g/L；生育期，土壤表层发生轻度、中度盐渍化时的地下水临界埋深分别为2.2、1.8 m，发生轻度、中度盐渍化时的地下水临界矿化度分别为2.5、3.0 g/L，春灌前更易发生土壤盐渍化，研究区北部、东南部和中部小部分区域为地下水小于临界埋深且大于临界矿化度的高概率区，是土壤返盐的高风险区，建议进一步完善排水系统，降低局部盐渍化风险区的地下水埋深和矿化度。

(4)春灌前，土壤盐分受外界因素(气象因素和人为因素)影响小，且较生育期土壤表层含盐量、地下水埋深和矿化度变异性较小，地下水环境对土壤盐渍化的影响更强烈。