孙贯芳， 高照良， 朱 焱， 杨金忠， 屈忠义
（1.西北农林科技大学水土保持研究所，陕西 杨凌 712100；
2.武汉大学水资源与水电工程科学国家重点实验室，湖北 武汉 430072；
3.内蒙古农业大学水利与土木建筑工程学院，内蒙古 呼和浩特 010018）
内蒙古河套灌区是中国最大的一首制自流灌区，盐碱地占全部灌溉耕地（57×104hm2）的50%以上［1］，土壤盐碱化是该地区灌溉农业发展面临的主要障碍［2］。精确评估和预测土壤盐分时空分布，实现农田盐分的精准管理和调控对提高农业生产力意义重大。然而，河套灌区土壤及地下水状况复杂，加之强烈的人类活动的影响，土壤盐分时空变异性强，实现其时空动态预测十分困难。学者们一直致力于用经典统计学和地统计学来研究土壤盐分的变异，并取得了一些成果［3-4］。在内蒙古河套灌区，徐英［5］基于非饱和带土壤垂向水分和盐分的时空变异性，用条件模拟理论再现了沙壕试验区季节性冻融土水盐的波动过程；
陈亚新等［6］研究了稳健统计学的水盐空间变差函数逼近方法；
刘全明［7］将指示克立格与神经网络技术融合，用于土壤水盐空间变异性评价；
史海滨等［8］用地统计学方法对比分析了盐渍化灌区节水改造前后土壤盐分时空变化规律；
王瑞萍等［9］利用经典统计学和地统计学方法研究了河套灌区乌拉特灌域春季土壤盐碱化空间分布特征。土壤盐分是具有时空动态变化属性的物理量［10］，目前对土壤盐分时空变化的评估大多仅考虑空间关系［3-4,11］。土壤盐分的时间变化趋势一般采用多个时期土壤盐分统计特征进行对比分析，该方法不适用于取样时间不规则且取样时间内空间位置不一致的数据集，尤其不适用于时间和空间上的稀疏数据集［12］。
时空克里金是地统计学向时空地统计学的延伸，其方差同时是时间和空间的函数，可为参数估计和预测提供更多的信息，许多研究已经证实时空克里金优于空间克里金［12-14］。到目前为止，时空克里金已经成功应用于环境科学［15-16］、气象学［17-18］、水文学［19］和土壤科学［20］。在时空克里金预测土壤盐分时空动态方面，Douaik 等［21］研究了匈牙利Great Hungarian平原一块25 hm2田块0～0.4 m土层的土壤盐分时空变化；
Gasch 等［22］利用和度量模型建立了美国华盛顿普尔曼附近一块37 hm2农田3维土壤电导率动态模型。以上2个土壤盐分时空克里金的研究均是田间尺度，区域尺度盐分变异性增强，时空克里金方法在区域土壤盐分时空建模与预测中的应用效果有待深入研究。
本文以内蒙古河套灌区隆胜研究区为典型研究区，基于2017—2018 年研究区内68 个监测点0～1.8 m 土壤剖面4582 个土壤盐分数据，利用时空地统计方法分析区域土壤盐分时空变化特征，比较时空克里金较传统空间克里金插值的精度提升效果，并验证时空地统计方法在监测点大幅减少情况下获取区域盐分时空动态的能力，研究成果可为土壤盐碱化时空动态建模与预测提供理论和科学依据。
1.1 研究区概况
河套灌区深处内陆，属于中温带干旱半干旱大陆性气候，冬季严寒少雪，夏季高温干旱。隆胜研究区是本研究的主要监测区，位于内蒙古河套灌区永济灌域，西边为永济干渠，北为永刚分干沟，东至东济支渠，南临永刚分干渠，西南东北长约15.5 km，西北东南宽约8.0 km，总土地面积8219.75 hm2（图1）。隆胜研究区每年11月中旬土壤开始封冻，次年5 月上旬融通。邻近的临河气象站年均降雨量148.8 mm，年均蒸发量（20 cm 蒸发皿）2327 mm。2017年、2018年降雨量分别为100.5 mm、176.2 mm，生育期降雨量分别为53.1 mm、156.6 mm。
研究区地质构造为河湖相交替沉积形成的湖相和河相沉积层，0～1.0 m 土壤质地以粉砂质壤土、粉土及砂质壤土为主，水平方向土壤质地分布极不均匀，空间变异强烈。0～2.5 m垂直方向上土体构型较为复杂，多有黏土、细砂夹层、砂黏土和粉砂土互层的土体结构［23］。
1.2 土壤盐分采样及测试
2017—2018 年在研究区内均匀布置了68 个农田土壤水盐观测点（图1），每个观测点距离900 m左右。分别于2017年5月初、2017年9月末、2018年5月初、2018 年9 月末对0～1.8 m 深度范围内每0.2 m一层的土壤水盐进行了4 次观测，分别标记为Y1705、Y1709、Y1805和Y1809。春季部分观测点地下水埋深浅，取到地下水面处停止取样。每个观测点用土钻取2个重复（孔），共获取4582个盐分有效样本。土壤盐分用电导率仪(上海雷磁DDSJ-308F)测试水土比为5:1的土壤浸提液的电导率。文中所用盐分值为2个重复的均值。
图1 隆胜研究区及盐分取样点位置图Fig.1 Longsheng study area and soil salinity sampling locations
图1 中的32 个长期监测点是根据改进时间稳定性方法确定的，将其用以验证时空地统计方法在监测点减少情况下获取区域盐分时空动态的能力，具体确定方法见文献［24］。
1.3 时空地统计
时空克里金方法可用来评估土壤盐分的时空动态。通常普通时空克里金有以下4 个步骤［14］：（1）时空经验半方差函数的计算；
（2）用理论时空半方差拟合经验半方差；
（3）用时空克里金做时空预测；
（4）用留一交叉验证方法评估预测结果。各步骤的详细计算方法如下：
用时空经验半方差函数γ(hS,hT)描述土壤参数时空分布的相关属性［13］：
式中：hS为两样本点空间分隔距离或滞后距离(m)；
hT为两样本点时间滞后距离(d)；
N(hS，hT)是时空间距为(hS，hT)时所有观测样点的对数；
Z(si，ti)和Z(si+hS，ti+hT)为土壤盐分在时空位置(si，ti)和(si+hS，ti+hT)的值。
用和度量模型拟合时空经验半方差γ(hS，hT)［10,15］，该模型将γ(hS，hT)分为空间方差γS(hS)、时间方差γT(hT)和联合方差γST(hST)，三者关系表示为：
式中：γS(hS)，γT(hT)和γST(hST)是各自独立的，通常由空间理论半方差函数模型如球状模型、高斯模型、指数模型来表示。在联合方差γST(hST)中常引入几何各项异性比α来计算时空滞后距离hST：
该研究区土壤盐分的3 个独立方差γS(hS)，γT(hT)和γST(hST)均用球状模型来拟合。球状模型的定义如下［24-25］：
式中：C0为块金值；
C为拱高；
C0+C为基台值；
a为相关长度(变程)；
i表示空间、时间和时空变量。参照空间相关度的概念，用各部分独立半方差的块金值之和与基台值之和的比值(STD)来评价时空相关性。STD＜0.25 表示时空相关性很强烈，STD＞0.75表示时空相关性较弱，其余表示中等强度时空相关性［26］。
克里金(Kriging)插值是对区域变量取值进行最优线性无偏估计的方法。时空克里金插值类似于空间克里金，一个栅格中的土壤盐分值Z*(s0，t0)由邻近的取样点盐分值(si，ti)进行估算［27］。
式中：λi为已知观测点对估值点的加权系数。根据无偏估计和方差最小两项要求可确定加权系数［28］，经转化可由以下矩阵方程求得权重系数：
式中：γij(hS，hT)为空间滞后距离为hS和时间滞后距离为hT的两点(si，ti)、(sj，tj)的时空半方差值；
μ为极小化处理时的拉格朗日乘数；
n为搜索域内的观测数。
留一交叉验证用来评估时空克里金插值效果。交叉验证时，验证点的观测值被移出数据集，用剩余的数据来预测该时空点的土壤盐分值［12］。然后，比较土壤盐分观测值和预测值的拟合优度来评估时空克里金模型。拟合优度评价指标用平均相对误差MRE、均方根误差RMSE 和决定系数R2，具体计算方法见文献［24］。
本研究基于Matlab 环境，采用遗传算法估计变异函数参数［29］，并完成试验半方差散点图、理论模型拟合曲面图的绘制；
时空插值结果分布图在Arc-GIS 10.3 中完成。土壤盐分的经典统计分析利用Excel 2019软件完成。
2.1 土壤盐分季节性变化特征
由表1 可知，不同土层土壤盐分的空间均值在0.25～0.35 dS·m-1之间变化，平均处于非盐碱化水平［1］。土壤盐分的峰度变化范围是0.71～18.52，偏度变化范围为1.05～4.18，数据右侧拖尾严重，所有取样时间各层土壤盐分均是非正态分布的。土壤盐分空间变异系数的变化范围是0.43～1.14，为中强变异。0～0.4 m越接近表层，土壤盐分受降雨、灌溉、蒸发等上边界作用越强烈，盐分随时间波动越剧烈。0.4～0.6 m 土壤盐分均值不同观测期变化范围为0.31～0.32 dS·m-1，变化幅度不超过5%，是盐分最为稳定的土层。2017年、2018年生育初期及生育末期0～0.6 m根系层土壤盐分分别为0.32 dS·m-1、0.34 dS·m-1、0.29 dS·m-1、0.30 dS·m-1，根系层土壤生育期处于积盐状态。以2017年生育期末和2018年生育期初根系层土壤盐分为例进行分析，秋浇和冻融作用使根系层土壤盐分降低了14.7%。这与孙贯芳等［30］在河套灌区秋浇灌溉洗盐试验的结果基本一致：秋浇灌黄河水180 mm 后，次年春播前0～1.0 m土壤盐分平均下降10.86%～26.14%。从整个剖面来看，5月初土壤盐分0.6～1.0 m土层盐分最大，两端较小，秋浇及冻融作用使盐分在根系层底部以下剖面聚集，生育期蒸发蒸腾作用又使部分盐分向上运动，导致0.6 m 以下各层土壤盐分生育初期明显大于生育末期。忽略2018年较大降雨对0～0.4 m土层盐分的淋洗作用，整个土壤剖面9 月底盐分分布自上到下依次减小，属于典型的灌溉蒸发型土壤盐分剖面。
表1 4次取样时间土壤盐分统计特征Tab.1 Statistical characteristics of soil salinity at four sampling times
2.2 土壤盐分时空半方差
用对数变换后的0～1.8 m 土壤盐分数据计算经验时空半方差，结果如图2。所有半方差在变程范围内均随着时间滞后距离和空间滞后距离的增大而增大，说明土壤盐分的相似性随着时间和空间距离的增大而减小。土壤盐分的时空半方差值逐渐达到一个常数值，对数转换后的土壤盐分数据在空间上是均匀的，并具有时间稳定性［12,14，31］。用和度量模型拟合经验时空半方差，并对所有的时空位置重复进行交叉验证，基于交叉验证结果选择的时空半方差模型如表2 所示。通常很难区分空间、时间和联合部分的块金值［15］，表中参数C0是3 部分块金值和，不同土层土壤盐分块金值的变化范围为0.0028～0.0186，基台值变化范围为0.0371～0.0618。0.4～0.6 m土层的时空相关度STD等于0.32，属于中等相关，其余的时空相关度变化范围为0.06～0.18，均小于0.25，时空相关性强烈。
表2 土壤盐分时空半方差函数模型及交叉验证结果Tab.2 The semivariogram models and parameters of soil salinity
图2 各层经验(蓝色圆圈)和拟合的土壤盐分时空半方差曲面Fig.2 Empirical(blue circle)and fitted spatio-temporal semivariogram surfaces of the soil salinity in each layer
2.3 土壤盐分时空分布的交叉验证与预测
土壤盐分时空克里金交叉验证结果如图3 所示，各层土壤盐分预测值和观测值间的MRE变化范围为3.38%～5.59%。0～0.2 m 土壤盐分的RMSE 是0.21 dS·m-1，其他层土壤盐分RMSE均小于0.15 dS·m-1。各层土壤盐分预测值和观测值间的决定系数R2在0.50～0.66 间变化。该土壤盐分交叉验证的R2结果远好于Gasch等［22］报道的土壤盐分时空克里金交叉验证结果0.18，基本上接近于Douaik 等［21］报道的土壤盐分时空克里金验证结果0.76。由图3 可知，除少数土壤盐分值较大而预测结果偏小的情况外，大部分预测的土壤盐分与实测数据吻合良好，和度量模型能够合理表征土壤盐分的时空结构，并取得精确的预测结果。
图3 土壤盐分时空克里金交叉验证结果Fig.3 Soil salinity cross-validation results for the spatio-temporal Kriging
此外，传统空间克里金方法计算的各时期半方差模型见文献［24］。土壤盐分的空间结构受自然和人类活动的影响，4 个取样时间有较大的差别。以0.4～0.6 m 土 层 为 例，Y1705、Y1709、Y1805、Y1809，4个时期变程分别为1610 m、1680 m、1801 m和6180 m，若取样时间刚好是Y1809，将会错误的认为该区域0.4～0.6 m 土壤盐分的变程是6000 m 左右，远超实际变程1600 m，用Y1809 取得的结果认识区域土壤盐分空间结构特征，将存在严重的偏差，仅由1 次取样获取土壤盐分的空间结构是不可靠的。因此，需多次监测综合分析确定区域土壤盐分的空间结构特征。时空克里金能够同时利用不同监测时间的盐分信息，是解决这一问题的有力手段。图4展示了传统空间克里金和时空克里金方法计算的各时期各土层间土壤盐分交叉验证的RMSE值，时空克里金的RMSE较传统空间克里金小0.02～0.09 dS·m-1，显著提高了土壤盐分预测精度。
图4 传统空间克里金和时空克里金土壤盐分交叉验证结果比较Fig.4 The cross-validation results of traditional spatial Kriging and spatio-temporal Kriging
基于半方差函数分析，用时空克里金插值方法评估了不同取样时间各层土壤盐分的空间分布。土壤盐分空间分布具有时间稳定性，不同时间盐分空间分布格局基本一致［24］。图5展示了Y1805各层土壤盐分的空间分布。同一时间不同层间土壤盐分分布具有相似性，各层土壤盐分高低值分布位置较为接近，特别是相邻层间盐分斑块位置较为一致，仅斑块的大小和范围有所差别，说明土壤盐分垂向分布较为连续，各层土壤盐分含量密切相关。0～0.6 m 各层土壤盐分高低相间，呈插花式分布，总体表现为研究区东部及东北部盐分含量高，而西北部和南部盐分含量较低。0.6～1.8 m 各层土壤盐分在研究区中部和北部沿东北-西南方向呈条带状分布，自西向东表现为低-高-低-高的盐分分布格局。该研究区土壤盐分空间分布格局受地势、土质、地下水埋深及矿化度等多种因素的影响，研究区东部排水沟经过位置处盐分含量反而较其他位置明显大，主要是由于该区域地势低，地下水埋深浅矿化度高，且土壤多为粉土和黏土［23］，排水沟年久失修未能有效发挥作用所致。
图5 0～1.8 m各层土壤盐分空间分布(Y1805)Fig.5 Soil salinity maps of different soil layers within the depth 0-1.8 m in Y1805
2.4 时空地统计获取稀疏监测点盐分时空动态的能力
0～0.6 m根系层和0.6～1.2 m土层时空半方差模型及参数见表2，留一交叉验证结果如表3。长期监测点确定的0～0.6 m和0.6～1.2 m土壤盐分预测值与观测值的MRE 分别为2.52%和2.54%，RMSE 为0.12 dS·m-1和0.09 dS·m-1，R2为0.73 和0.70。长期监测点确定的根系层土壤盐分空间分布格局(图6a1～6d1)和所有取样点确定的空间分布图(图6a2～6d2)基本一致。如图7a、图7c 所示，由土壤盐分长期监测点确定的不同盐碱化水平的面积与所有取样点确定的面积吻合较好，两者间的MRE 为-13.20%，RMSE 为466.67 hm2，R2为0.79，说明由土壤盐分长期监测点确定根系层土壤盐分的空间分布是可行的。同样地，由土壤盐分长期监测点确定的0.6～1.2 m 土壤层盐分分布(图6e2～6h2)和所有取样点确定的盐分分布图(图6e1～6h1)也高度一致。两者确定的相同盐碱化水平的面积统计结果也较为接近(图7b、图7d)，面积间的MRE 为-8.35%，RMSE 为494.43 hm2，R2为0.72，说明土壤盐分长期监测点可以用来确定0.6～1.2 m 土壤层盐分的空间分布。进一步比较时空克里金与普通克里金获取稀疏监测点盐分时空动态的能力，两者土壤盐分长期监测点确定的不同盐碱化水平的面积与所有取样点面积间的评价指标结果见表4，时空克里金的结果明显优于普通克里金，说明在数据监测点较为稀疏的情况下，相较于普通克里金，时空克里金也提高了稀疏监测点土壤盐分时空动态的预测精度。
表4 时空克里金和普通克里金获取稀疏监测点土壤盐分时空动态能力的比较Tab.4 Comparison of ability of spatio-temporal Kriging and ordinary Kriging to obtain soil salinity dynamics under sparse monitoring locations
图6 由土壤盐分长期监测点(LSML)和所有监测点(ASML)确定的0～0.6 m根系层和0.6～1.2 m土壤层盐分空间分布比较Fig.6 The spatial distribution of root zone(within the depth of 0-0.6 m)and 0.6-1.2 m layer soil salinity determined by all soil salinity monitoring locations(ASML)and long-term soil salinity monitoring locations(LSML)at the 4 sampling times
图7 由土壤盐分长期监测点(LSML)和所有取样点(ASML)确定的不同土壤盐碱化面积的比较Fig.7 The area of different soil salinity determined by all soil salinity monitoring locations(ASML)and long-term soil salinity monitoring locations(LSML)
表3 由所有盐分监测点(ASML)和长期监测点(LSML)确定的土壤盐分时空克里金验证结果Tab.3 Cross-validation results for 0-0.6 m and 0.6-1.2 m layer soil salinity with all soil salinity monitoring locations(ASML)and long-term soil salinity monitoring locations(LSML)
综上所述，由土壤盐分长期监测点确定的土壤盐分空间分布状况与所有取样点较为一致，利用时空克里金在土壤盐分监测点数据减少为原来一半情况下仍可以精确估计土壤盐分时空动态，以较少的采样数据实现了土壤盐分时空变化特征的精确评估，极大提高区域土壤盐分监测的效率。
本文以内蒙古河套灌区隆胜研究区为典型研究区，基于2017—2018 年研究区内68 个监测点0～1.8 m土壤盐分数据，分析了土壤盐分的季节性变化特征，利用时空地统计方法研究了区域土壤盐分时空变化特征及时空克里金的插值精度，并进一步验证时空地统计方法在监测点不足原来一半情况下获取区域盐分时空动态的能力。主要结论有：
（1）该研究区土壤盐分有明显的季节性规律，0～0.6 m 根系层生育期积盐，非生育期脱盐；
0.6～1.8 m土壤剖面生育期脱盐，非生育期积盐。
（2）和度量模型能较好拟合盐分时空经验半方差，时空克里金对土壤盐分交叉验证的RMSE 较传统空间克里金小0.02～0.09 dS·m-1，时空克里金同时利用了土壤盐分时间和空间上的更多信息，较普通克里金显著提高土壤盐分的预测精度。
（3）同一时间不同层间土壤盐分分布具有相似性，相邻土层间盐分斑块位置较为一致，仅斑块的大小和范围有所差别。0～0.6 m 各层土壤盐分高低相间，呈插花式分布，0.6～1.8 m 各层土壤盐分在研究区中部和北部沿东北-西南方向呈条带状分布，自西向东表现为低-高-低-高的盐分分布格局。
（4）以不足原来所有观测点数量一半的稀疏监测点土壤盐分数据，结合时空克里金计算的土壤盐分空间分布状况与所有取样点结果较为一致，实现了土壤盐分时空动态的精确估计，极大提高区域土壤盐分监测的效率。
猜你喜欢克里盐分方差我可以咬一口吗？知识窗(2023年2期)2023-03-05概率与统计(2)——离散型随机变量的期望与方差中学生数理化(高中版.高考数学)(2021年3期)2021-06-09你今天真好看风流一代·经典文摘(2019年12期)2019-09-10方差越小越好?中学生数理化·七年级数学人教版(2019年6期)2019-06-25计算方差用哪个公式中学生数理化·七年级数学人教版(2019年6期)2019-06-25你今天真好看读者(2018年24期)2018-12-04方差生活秀初中生世界·九年级(2017年10期)2017-11-08要借你个肩膀吗？知识窗(2017年3期)2017-03-09长期膜下滴灌棉田根系层盐分累积效应模拟中国农业文摘-农业工程(2016年5期)2016-04-12摄影欣赏读者·校园版(2016年6期)2016-03-07