基于指示克里格的青铜峡灌区土壤盐渍化风险与潜水埋深和矿化度关系研究

赵志鹏，黄小琴，方 磊

（宁夏回族自治区水文环境地质调查院，宁夏 银川 750011）

0 引 言

土壤是世界上最复杂多样的生态系统［1］。土壤盐渍化是可溶性盐分在土壤中积累致使土壤特性变化和质量下降的过程［2］。我国盐渍土面积达到3.69×107hm2，不仅是制约农业可持续发展的主要障碍，也可能造成经济损失、次生危害和生态环境脆弱［3，4］。在任何引水灌溉系统的干旱半干旱地区都会存在因盐渍化引起的土地退化问题［5］，与这些区域地势低平，排水不畅，地下水埋藏较浅关系很大［6］。地下水作为盐分传输、累积和排泄的主要载体，与土壤盐分关系密切，特别是潜水位埋深和矿化度对土壤盐分影响显著［7］。因此，在综合考虑土壤盐渍化与地下水关系的基础上，分析评价土壤盐渍化风险，对农业生产、水土资源利用和环境保护等方面均有重要意义。但在进行区域性土壤、地下水等特性的分析中经常会出现样本数据中存在比算术平均值或中位数要大（或小）得多的特异值，从而使样本偏离正态（或对数正态）分布，影响变异函数的稳健性。指示克立格是一种非参数估计方法，在处理有偏数据的明显优势，使其在生态学［8，9］、地质学［10-12］、环境科学［13-15］等领域被广泛应用。采用指示克里格方法开展土壤盐渍化风险评价的研究较多，周在明在环渤海低平原区的研究表明地下水矿化度大于2 g/L，地下水位埋深小于3 m 以及土壤盐分含量大于1 g/kg 的概率分布呈一致性［16］；姚荣江等的研究证明在黄河三角洲典型地区地下水埋深与土壤盐分的概率分布存在空间上的规律性与相似性［17］；李仙岳等利用指示克里格法比较了永济灌域春灌前和生育期不同阈值条件下土壤表层含盐量、地下水埋深和矿化度的概率分布，从概率空间分布的角度提出了不同时期防治土壤盐渍化的地下水临界埋深和矿化度［18］。相对其他土壤盐渍化发育地区，青铜峡灌区近年来土壤盐渍化的研究主要集中于利用多源遥感、实测高光谱数据开展盐分反演、评价及动态监测［19］，草甸湿地等特殊土地利用类型盐分特征分析［20］，灌区盐渍土空间分布特征［21］等方面，针对整个灌区尺度的盐渍化风险评价研究较少。本文以青铜峡灌区为研究对象，运用指示克里格法分析不同深度土壤盐分、潜水位埋深和矿化度的空间概率分布规律，分析土壤盐渍化与潜水的空间分布关系，以期为地下水资源调控、土地资源高效利用、盐渍土防治等提供参考。

1 研究区概况

青铜峡灌区南起青铜峡水利枢纽，北至石嘴山，西抵贺兰山，东邻鄂尔多斯台地西缘，可分为银北、银川和银南3 个片区（图1）。灌区多年平均降水量178 mm，年蒸发量接近年降水量的10 倍，日照时数2 870～3 080 h。灌区分布的土壤以灌於土、灰钙土为主。粮食作物以水稻、玉米、小麦为主，经济作物以蔬菜和瓜果类为主。青铜峡灌区主要由黄河冲积平原构成，海拔1 100～1 200 m，地势由南向北缓缓倾斜，坡降约1/4 000。灌区地下水类型为第四系松散岩类孔隙水，分为单一潜水区和多层结构区，前者位于灌区最南端，系黄河出青铜峡峡口所形成的冲积扇，从西南向东北，含水层岩性由卵砾石逐渐变为含砾粉细砂，含水层厚度变厚；多层结构区分布在单一潜水区以外地区，在大约250 m 深度以上的范围内，从上向下划分为潜水、第一承压、第二承压含水岩组，各含水岩组间通常具有相对较连续的弱透水层，其中潜水含水层厚度一般为20～60 m，含水岩组岩性表现出自灌区边缘向中心及由南向北逐渐变细、淤泥质含量增多的特点。引黄灌溉渠系渗漏补给和田间入渗补给占灌区地下水补给量的80%左右，潜水蒸发占地下水排泄量的47%左右［22］。1998-2017 年青铜峡灌区潜水位埋深下降了0.69 m，增加速率为0.038 m/a，年内地下水埋深呈双峰双谷特征［23］。2019 年潜水溶解性总固体为0.21～21.95 g/L。土壤盐渍化一直是制约青铜峡灌区农业发展的重要因素，而盐渍化的形成演化过程与平缓的地形、较浅的潜水位埋深、强烈的垂向蒸发排泄和局部地下水的高矿化度，以及引黄灌溉导致的盐分积累等因素密切相关。

图1 土壤采样点和潜水监测点分布图Fig.1 Distribution of soil sampling points and phreatic water monitoring points

2 数据与方法

2.1 数据来源

本研究使用了土壤可溶盐、潜水位埋深、潜水矿化度3类数据（图1）。其中:土壤可溶盐数据来自于2018 年4 月和2019 年4 月春灌前开展的2 次野外采集样品的测试，通过综合考虑灌区农业种植结构、土地利用类型和地下水状况等条件，尽可能均匀地布设了80 个采样点，分层采集了0～20、20～40、40～60 cm深度土壤可溶盐样品；各采样点使用GPS 定位，并对当地土地耕作、灌溉制度、农作物发育特征等信息进行调查；采集的土壤样品由宁夏回族自治区地质矿产中心实验室按照《土工试验方法标准（GB/T50123-1999）》测试，其中K++ Na+使用火焰光度法测定、SO2-4用比浊法测定、Ca2+和Mg2+使用EDTA 滴定法测定、HCO-3和CO2-3用双指示剂滴定法测试、Cl-用AgNO3滴定法测试、pH用pH计测定，采用离子加和法计算土壤全盐量。潜水位埋深数据来自于灌区内117眼国家级和省级监测井自动化监测数据，另外选用了218 个同时期潜水位调查点数据作为补充，数据时段均为2019年4月25-30日。潜水矿化度数据来自于灌区内103 眼国家级监测井水质测试报告，数据时段为2019 年4月。国家级监测点潜水位埋深和水质数据按照国家地下水监测工程要求的标准开展全流程质量控制［24］，补充的潜水位调查点采用水位测量仪测量地下水埋深，使用RTK 获取点的位置和高程信息。

2.2 研究方法

采用软件GS+9.0 确定不同阈值下土壤全盐量、潜水位埋深、潜水矿化度的变异函数模型，并将其模型参数输入Arc-GIS10.6 进行指示克里格插值，进而绘制概率分布图。步骤是:①根据研究需要确定土壤全盐量、潜水位埋深、潜水矿化度评价的阈值。②根据设定阈值利用指示函数对土壤全盐量、潜水位埋深和潜水矿化度数据进行二态指示变换，将得到的指示变换值（即1 或0）作为评价数据。③在GS+9.0 对获取的指示变换值进行空间变异分析，获取最佳变异函数模型及其参数。④将变异函数模型参数输入ArcGIS10.6 进行指示克里格插值，得到各指标满足相应阈值的空间概率分布图。经典统计及相应的图件制作采用SPSS25.0、GraphPad Prism 9.3.1软件进行。

3 结果与分析

3.1 统计特征值分析

表1 为统计特征值分析结果。0～20、20～40、40～60 cm 深度土壤全盐量平均值分别为12.27、4.85、4.75 g/kg，潜水位埋深、潜水矿化度的平均值分别为4.61 m 和2.55 g/L。变异系数Cv反映样点的离散程度，对于估计结果可以提供一些预警信息，Cv＜0.1为弱变异性、0.1≤Cv≤1为中等变异性，Cv＞1为强变异性并标识存在一些特别大的样本值［25］。青铜峡灌区不同深度土壤全盐量、潜水位埋深、潜水矿化度的Cv均大于1，都属于强变异，造成这种现象的原因首先与研究区范围较大而产生的地形地貌、地下水赋存条件、土壤及地层特征等结构性因素有关，另外也与灌排方式、种植结构、耕种方式等随机性因素的作用有关。单样本K-S 正态检验结果表明，均不服从正态分布或对数正态分布。从图2 可知，0～20、20～40、40～60 cm 土壤全盐量频率分布分别集中在小于上四分位数的10.35、4.69、4.98 g/kg 区间，潜水位埋深和潜水矿化度频率分布也集中在小于上四分位数的4.85 m 和2.44 g/L 的区间，由于都在局部地区存在高值，频率分布呈现明显的“高顶”现象。

表1 青铜峡灌区土壤全盐量、潜水位埋深和潜水矿化度统计特征值Tab.1 Descriptive statistics of total soil salt、groundwater depth and TDS in Qingtongxia Irrigation Area

图2 青铜峡灌区土壤全盐量、潜水位埋深和潜水矿化度小提琴图Fig.2 Violin plot of total soil salt， groundwater depth and TDS in Qingtongxia Irrigation Area

3.2 指示变异函数模型构建

3.2.1 阈值的选择

综合考虑青铜峡灌区地形地貌、水文地质条件以及土壤盐渍化程度等多方面因素确定各变量阈值。其中根据中国盐渍土划分标准［2］和相关研究成果［26］，选取土壤全盐量1.0、2.0、4.0、10.0 g/kg 分别作为土壤轻度、中度、重度盐化及盐土化的阈值；参考Zhang 等人的研究［27］，选取潜水矿化度为1.0、2.0、2.5、3.0 g/L 为水质咸化阈值；参考金晓媚［28］、邓丽［29］、周在明［16］等人研究，选取1.5、2.0、2.5、3.0 m作为潜水蒸发作用的临界水位埋深。二态指示变换规则是，土壤全盐量和潜水矿化度大于相应阈值为1、小于相应阈值为0；潜水位埋深为大于相应阈值为0、小于相应阈值为1。

3.2.2 指示变异函数模型

指示克立格法通过指定阈值对数据进行二态指示变换，可以削弱偏态分布和特异值对变异函数稳健性产生的影响。表2是各变量不同阈值指示变换后的指示变异函数模型和相关参数。可以看出:青铜峡灌区不同深度、不同阈值条件下土壤全盐量的指示半方差函数符合球状模型、指数模型、高斯模型等不同类型。潜水位埋深的指示半方差函数符合指数模型。潜水矿化度阈值为1.0 和2.0 g/L 时指示半方差函数符合指数模型，阈值为2.5 和3.0 g/L 时指示半方差函数符合高斯模型。从决定系数和残差看，理论模型均显示出较好的拟合程度。

表2 不同阈值下土壤全盐量、潜水位埋深和潜水矿化度的指示变异函数理论模型Tab.2 Indicator semivariogram models of soil salinity， phreatic water depth and phreatic TDS for different threshold values

3.2.3 不同变量的空间异质性分析

块金值C0由测量误差和空间变异共同作用引起，反映区域化变量的随机异质性程度；基台值Sill表示结构性和随机性变异构成的总变异，表示变量空间异质性的强弱，C0/Sill为＜25%、25%～75%、＞75%对应空间相关性较强、中等和很弱，在整个尺度上具有恒定变异时其值将接近100%。变程是不同变量空间异质性的尺度函数，用来衡量最大变异程度的空间距离，在变程之内具有空间自相关性［30］。区域化变量的空间变异是由结构性、随机性因素共同引起的，其中气候、地形地貌、土壤母质、地质与水文地质条件等非人为的结构性因素是变量具有空间连续性的原因，而灌溉制度、种植结构、耕作措施、地下水开发等随机因素则会弱化其空间自相关性。根据表2，各模型块金值C0都较小，说明在本研究尺度上由于采样误差、短距离的变异、随机和固有变异引起土壤全盐量、潜水位埋深和潜水矿化度的变异都不大。土壤全盐量C0/Sill的不同反映出在不同深度由空间自相关引起的空间异质性所占的比例有差异，总体表现出随土层深度增加呈倒U 型特点:0～20 cm 阈值为1.0 g/kg，20～40 cm 阈值为1.0 g/kg 和4.0 g/kg 和40～60 cm 阈值为4.0 g/kg 时，C0/Sill均小于25%，表现出较强的空间自相关，即空间变异主要由结构型因素作用；而20～40 和40～60 cm 阈值为10.0 g/kg 时，C0/Sill为100%，说明在整个研究尺度具有恒定变异；其余条件下C0/Sil处于25%～75%之间，表现出中等的空间自相关，也即空间变异由结构型因素和随机性因素共同作用。潜水位埋深阈值为最小的1.5 m 和最大的3.0 m 时，C0/Sill分别为5.82%和12.59%，表现出较强的空间自相关，空间变异主要由水文地质条件等结构型因素作用，而阈值为2.0 和2.5 m 时，C0/Sill处于25%～75%之间，表现出中等的空间自相关，说明人类活动等随机性因素对空间变异的影响程度增加。研究区潜水矿化度的C0/Sill在4种阈值条件下均＜25%，表现出较强的空间自相关，即其空间变异主要由结构性因素控制的，具体表现在由于含水层岩性、透水能力、水力坡度、水循环能力等的不同，造成灌区单一潜水区和多层结构区潜水位埋深和TDS 特征的明显区别，即单一浅水区潜水埋藏较深而TDS 多小于0.5 g/L；而在多层结构区，潜水的水位埋深和水化学成分在空间上的变化都较大，永宁县以南地下水水力坡度大、径流条件好、水交替作用强烈，大部分为低矿化淡水，而在永宁县以北随着地下水径流条件的逐渐变差和水位埋深的逐渐变浅，地下水质逐渐变差。另外，从表2 知，潜水位埋深和潜水矿化度的变程分别在8.70～45.30 和12.99～19.20 km，自相关范围均不大（约为青铜峡灌区最大直线范围的1/10 到1/3）。不同深度土壤全盐量的变程总体随着阈值的变大而增大，比如0～20 cm 阈值为1.0 和2.0 g/kg 时自相关范围较小，而阈值为4.0 和10.0 g/kg 时，自相关距离已超过青铜峡灌区最大直线范围的1/2，反映出当阈值较小时表层土壤全盐量受自然因素和人类活动的影响较大，而随着阈值的增大土壤自身性状对其盐分特征的影响变大。

3.3 土壤全盐量与潜水环境的概率空间分布特征

将变异函数模型参数输人ArcGIS10.6 进行指示克里格插值，绘制各指标满足相应阈值的概率空间分布图，并定义概率大于40%为高概率区、反之为低概率区。图3（a）～（d）、（e）～（h）和（i）～（l）分别是不同深度土壤处于轻度（＞1 g/kg）、中度（＞2 g/kg）、重度盐化（＞4 g/kg）和盐土化（＞10 g/kg）的概率分布图，表3是不同阈值条件下各指标高低概率区面积统计情况。可以看出，春灌前青铜峡灌区0～20 cm 深度土壤发生轻度、中度、重度盐化和盐土化的高概率区占灌区总面积的比例为92.38%、73.71%、51.83%、24.61%；20～40 cm深度土壤发生轻度、中度、重度盐化和盐土化的高概率区占灌区总面积的比例为92.15%、50.21%、6.87%、0.53%，40～60 cm深度土壤发生轻度、中度、重度盐化和盐土化的高概率区占灌区总面积的比例为89.16%、43.06%、6.09%、0.46%；即随着阈值增大，各深度高概率区明显减少，这与灌区盐化表现出的表聚型特点相符［21］。对比图1 发现，当土壤深度为0～20 cm，阈值为2 g/kg 时，银北片区的大部分、银川片区的全部均为高概率区，其面积约占灌区总面积的73.71%；当阈值为4 g/kg 时，高概率区面积减少至3 174.25 km2，主要分布在银川片区全部区域和银北片区的局部；当阈值为10 g/kg 时，除银川片区中部和银北片区石嘴山南部的小范围外约1 507.33 km2区域外，其余均演化为低概率区。当土壤深度为20～40和40～60 cm，发生中度盐化（阈值为2 g/kg）的高概率区分别占灌区总面积的50.21%和43.06%，主要分布在银北片区东北部、西南部和银川片区的北部；发生重度盐化（阈值为4 g/kg）的高概率区面积分别减少至420.74 和372.84 km2，主要集中在银北片区的西南部；发生盐土化（阈值为10 g/kg）的高概率区面积进一步减少至32.48 和28.45 km2，占比也降至0.50%左右。总的来看，不同深度大阈值的高概率区均包含在上一级小阈值的高概率区，并在空间分布上具有明显的相似性，尤其是深部土壤重度盐化和盐土化高风险区主要集中银北片区的西南部一带。

图3 青铜峡灌区不同阈值下不同深度土壤全盐量概率空间分布Fig.3 Probability maps of soil salinity at different depths under different thresholds in Qingtongxia Irrigation Area

图4（a）～（d）为不同阈值潜水位埋深的概率分布，对比图1发现，不同阈值概率分布图的共同特点是:银北片区是高概率集中分布区，低概率区主要集中在银川片区中部和银南片区。当潜水位埋深阈值为1.5 m 时，研究区潜水埋深概率分布区间主要在0～0.2 的低概率之间，仅13.42%的区域为高概率区且主要分布在银北片区的东北部和银南片区东南部局部小区域；当阈值为2.0 m 和2.5 m 时时，高概率区面积在原有基础上分别扩大至约1 852.48、2 564.33 km2，占比增加至30.25%和41.87%，增加部分均主要来源于银北片区；随着阈值的增加为3.0 m 时，高概率风险区继续增加，仅约34.82%的区域为低概率区且主要集中分布于银川片区中的城市建设区周围和银南片区的单一潜水区，其余地区基本为高概率区。分析不同阈值下潜水位埋深概率图，发现小阈值的高概率区包含在大阈值的高概率区，其分布区域随阈值的增大而扩散。

图4 青铜峡灌区不同阈值下潜水位埋深概率空间分布Fig.4 Probability spatial distribution of phreatic water depth under different thresholds in Qingtongxia Irrigation Area

图5（a）～（d）为不同阈值下潜水矿化度的概率分布，当阈值为1.0、2.0、2.5、3.0 g/L 时高概率区面积分别是3 621.55、3 381.99、2 067.43、1 957.72 km2，占比分别为59.13%、55.22%、33.76%和31.96%。结合图1发现，不同阈值下银川片区和银北片区都是高概率集中分布区，低概率区主要集中于银南片区，同时表现出大阈值的高概率区包含在小阈值的高概率区，其分布区域随阈值的增大逐渐缩小的变化特点。

图5 青铜峡灌区不同阈值下潜水矿化度概率空间分布Fig.5 Probability spatial distribution of phreatic TDS under different thresholds in Qingtongxia Irrigation Area

3.4 土壤全盐量与潜水环境的概率空间分布关系

分析图3、图4 和图5 发现3 个不同深度土壤发生轻度盐化（阈值为1 g/kg）高风险区基本覆盖整个灌区，这与由于本研究土壤可溶盐样品采集于4月份的春灌前，此时土壤处于融冻期，返盐激烈且不能及时脱盐导致盐分积累加重的事实有关；另外3个深度发生盐土化（阈值为10 g/kg）的范围分布很小。即分析不同深度土壤发生轻度盐化和盐土化的概率空间分布，与潜水埋深和矿化度的关系意义不大。以下重点对比深度为20～40 cm 时土壤发生中度（阈值为2 g/kg）、重度盐化（阈值为4 g/kg）的概率空间分布与潜水位埋深和矿化度之间的关系。使用ArcGIS10.6 的Intersect 功能计算不同阈值土壤全盐量高概率区与潜水位埋深、矿化度高概率区相交的部分，即同时满足土壤全盐量相应阈值和潜水环境相应阈值高概率分布的区域，并统计匹配率发现（表4），青铜峡灌区潜水位埋深越浅，表层土壤含盐量越高，反之愈轻；而潜水矿化度愈大，土壤盐渍化愈严重，反之也愈轻。这与王金哲等在西北干旱区［31］和樊自立等在塔里木河流域［32］的研究结论相同。可能的原因是潜水埋深越浅，地下水中的盐分越易通过潜水蒸发由地下水带至土壤耕层而表聚；特别是当潜水矿化度较大时，越易发生次生盐渍化。因此可以认为存在临界埋深和矿化度，即当潜水位埋深小于某临界值时、矿化度大于某临界值时，土壤表层发生盐渍化的风险就大。对比图3 和图4 发现，潜水埋深2.0 m 及以上阈值的概率空间分布与土壤全盐量2.0 g/kg和4.0 g/kg的概率空间分布有较大相似性，且随着潜水埋深的增大匹配率逐渐增大，并且明显看出潜水埋深2.0 m 是匹配率突增的临界值，可以初步判定发生中度和重度盐化的潜水埋深临界值为 2.0 m，与黄权中在河套灌区的研究结果相同［33］。再对比图3 和图5 可知，土壤全盐量阈值为2.0 g/kg 的高概率分区范围与潜水矿化度1.0 g/L 的相似性最大，匹配率达到74.22%；与潜水矿化度2.0 g/L 的相似性也较大，匹配率达到69.34%；随着潜水矿化度阈值增加至2.5 g/L，匹配率大幅降低至39.02%；土壤全盐量阈值为4.0 g/kg 的高概率分区范围与潜水矿化度2.0 g/L 的相似性最大，匹配率达到59.45%；与潜水矿化度1.0、2.5 和3.0 g/L 匹配率都仅为30%左右；即可以初步判定发生中度和重度盐化的潜水矿化度临界阈值为2.0 g/L。以上分析表明青铜峡灌区土壤盐分的空间变异性与潜水位埋深和矿化度的空间分布密切相关，但同时也发现，土壤全盐量的高概率空间分布并不与潜水位埋深和矿化度的高概率空间分布完全相同，这是因为土壤盐分特征还受到的气象、地形地貌、土地利用方式等因素的共同影响［21］，而农业生产活动的频繁程度、作物类型、施肥措施和灌溉水量及时间也是影响土壤盐分的重要因素［25，34］。因此在青铜峡灌区影响农田土壤盐渍化的一个重要因素是灌区潜水埋深和矿化度。因此，有必要进一步健全灌排系统，将大水漫灌用田间畦灌、滴灌、喷灌等节水型灌溉措施来替代，同时结合其他土壤改良方法，多渠道联合控制土壤盐渍化。

表4 20～40cm土壤全盐量与潜水环境匹配率统计表Tab.4 Statistical table of matching rate between total salt content of 20～40 cm soil and phreatic environment

4 结 论

根据对青铜峡灌区2018、2019两年4月份春灌前0～20、20～40 和40～60 cm 深度土壤全盐量、潜水水位埋深、潜水矿化度进行空间变异性分析表明。

（1）青铜峡灌区土壤全盐量、潜水位埋深、潜水矿化度的变异系数均大于1，属于强变异；频率分布呈现明显的“高顶”现象。

（2）青铜峡灌区不同深度不同阈值土壤全盐量的指示半方差函数符合球状模型、指数模型、高斯模型等不同类型。潜水位埋深的指示半方差函数符合指数模型。潜水矿化度阈值为1.0 和2.0 g/L 时指示半方差函数符合指数模型，阈值为2.5 和3.0 g/L 时指示半方差函数符合高斯模型。不同阈值条件下的土壤全盐量和潜水位埋深呈中等—较强的空间自相关，潜水矿化度表现出较强的空间自相关。

（3）青铜峡灌区0～20 cm 深度土壤发生轻度、中度、重度盐化和盐土化的高概率区占灌区的92.38%、73.71%、51.83%、24.61%；20～40 cm 深度比例为92.15%、50.21%、6.87%、0.53%；40～60 cm 深度比例为89.16%、43.06%、6.09%、0.46%。潜水位阈 值 为1.5、2.0、2.5、3.0 m 时 高 概 率 区 占 灌 区 的13.42%、30.25%、41.87%、65.18%。潜水矿化度阈值为1.0、2.0、2.5、3.0 g/L 时高概率区面积占灌区的59.13%、55.22%、33.76% 和31.96%。随着阈值增加，高风险区面积均明显减少。

（4）青铜峡灌区土壤发生中度和重度盐渍化时的潜水临界埋深为2.0 m，潜水临界矿化度为2.0 g/L。