时空克里金评估河套灌区土壤盐分时空格局_孙贯芳.pdf

第40卷第2期2023年2月Vol.40No.2Feb.2023干 旱 区 研 究ARIDZONERESEARCHhttp:/DOI：10.13866/j.azr.2023.02.03时空克里金评估河套灌区土壤盐分时空格局孙贯芳1，高照良1，朱焱2，杨金忠2，屈忠义3（1.西北农林科技大学水土保持研究所，陕西 杨凌712100；2.武汉大学水资源与水电工程科学国家重点实验室，湖北 武汉430072；3.内蒙古农业大学水利与土木建筑工程学院，内蒙古 呼和浩特010018）摘要：区域土壤盐分时空变异性大，采用经典统计和地统计方法无法准确判断取样时间不规则、空间位置不一致的土壤盐分的时空变化趋势。本文以内蒙古河套灌区隆胜研究区68个监测点01.8 m土壤剖面4582个土壤盐分数据为基础，利用时空地统计方法分析区域土壤盐分时空变化特征，比较时空克里金较传统空间克里金插值的精度提升效果，并验证时空地统计方法在监测点减少50%情况下预测区域盐分时空动态的能力。结果表明：（1）该研究区土壤盐分空间变异系数的变化范围是0.431.14，为中强变异，00.6 m根系层生育期积盐、非生育期脱盐，0.61.8 m土壤剖面生育期脱盐、非生育期积盐、农田土壤盐分有明显的季节性规律。（2）和度量模型能较好拟合盐分时空经验半方差，各层土壤盐分预测值和观测值间的均方根误差RMSE均小于0.21 dSm-1，较传统空间克里金的RMSE小0.020.09 dSm-1。（3）采用该方法在减少50%监测点情况下确定的土壤盐渍化分布与所有取样点确定的结果一致性较高，00.6 m和0.61.2 m土壤盐分面积间的相对误差MRE分别为-13.20%和-8.35%，RMSE为466.67hm2和494.43 hm2，决定系数R2为0.79和0.72。时空克里金同时利用了土壤盐分时间和空间上的更多信息，实现了稀疏盐分监测点数据集土壤盐分时空动态的精确估计，可极大提高区域土壤盐分时空格局监测的效率。关键词：土壤盐分；时空克里金；河套灌区；时空格局；和度量模型内蒙古河套灌区是中国最大的一首制自流灌区，盐碱地占全部灌溉耕地（57104hm2）的50%以上1，土壤盐碱化是该地区灌溉农业发展面临的主要障碍2。精确评估和预测土壤盐分时空分布，实现农田盐分的精准管理和调控对提高农业生产力意义重大。然而，河套灌区土壤及地下水状况复杂，加之强烈的人类活动的影响，土壤盐分时空变异性强，实现其时空动态预测十分困难。学者们一直致力于用经典统计学和地统计学来研究土壤盐分的变异，并取得了一些成果3-4。在内蒙古河套灌区，徐英5基于非饱和带土壤垂向水分和盐分的时空变异性，用条件模拟理论再现了沙壕试验区季节性冻融土水盐的波动过程；陈亚新等6研究了稳健统计学的水盐空间变差函数逼近方法；刘全明7将指示克立格与神经网络技术融合，用于土壤水盐空间变异性评价；史海滨等8用地统计学方法对比分析了盐渍化灌区节水改造前后土壤盐分时空变化规律；王瑞萍等9利用经典统计学和地统计学方法研究了河套灌区乌拉特灌域春季土壤盐碱化空间分布特征。土壤盐分是具有时空动态变化属性的物理量10，目前对土壤盐分时空变化的评估大多仅考虑空间关系3-4,11。土壤盐分的时间变化趋势一般采用多个时期土壤盐分统计特征进行对比分析，该方法不适用于取样时间不规则且取样时间内空间位置不一致的数据集，尤其不适用于时间和空间上的稀疏数据集12。时空克里金是地统计学向时空地统计学的延伸，其方差同时是时间和空间的函数，可为参数估计和预测提供更多的信息，许多研究已经证实时空克里金优于空间克里金12-14。到目前为止，时空克里金已经成功应用于环境科学15-16、气象学17-18、水文学19和土壤科学20。在时空克里金预测土壤盐收稿日期：2022-08-12；修订日期：2022-10-12基金项目：国家重点研发计划项目（2021YFD1900805）；中央高校基本科研业务费专项资金（2452021081）；内蒙古自治区科技成果转化专项资金（2021CG0022）。作者简介：孙贯芳（1989-），男，助理研究员，主要研究方向为农业水土资源与环境.E-mail:通讯作者：朱焱.E-mail:182193页2期孙贯芳等：时空克里金评估河套灌区土壤盐分时空格局分时空动态方面，Douaik等21研究了匈牙利GreatHungarian平原一块25 hm2田块00.4 m土层的土壤盐分时空变化；Gasch等22利用和度量模型建立了美国华盛顿普尔曼附近一块37 hm2农田3维土壤电导率动态模型。以上2个土壤盐分时空克里金的研究均是田间尺度，区域尺度盐分变异性增强，时空克里金方法在区域土壤盐分时空建模与预测中的应用效果有待深入研究。本文以内蒙古河套灌区隆胜研究区为典型研究区，基于20172018年研究区内68个监测点01.8 m土壤剖面4582个土壤盐分数据，利用时空地统计方法分析区域土壤盐分时空变化特征，比较时空克里金较传统空间克里金插值的精度提升效果，并验证时空地统计方法在监测点大幅减少情况下获取区域盐分时空动态的能力，研究成果可为土壤盐碱化时空动态建模与预测提供理论和科学依据。1材料与方法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。研究区地质构造为河湖相交替沉积形成的湖相和河相沉积层，01.0 m土壤质地以粉砂质壤土、粉土及砂质壤土为主，水平方向土壤质地分布极不均匀，空间变异强烈。02.5 m垂直方向上土体构型较为复杂，多有黏土、细砂夹层、砂黏土和粉砂土互层的土体结构23。1.2 土壤盐分采样及测试20172018年在研究区内均匀布置了68个农田土壤水盐观测点（图1），每个观测点距离900 m左右。分别于2017年5月初、2017年9月末、2018年5月初、2018年9月末对01.8 m深度范围内每0.2 m一层的土壤水盐进行了 4 次观测，分别标记为Y1705、Y1709、Y1805和Y1809。春季部分观测点地下水埋深浅，取到地下水面处停止取样。每个观测点用土钻取2个重复（孔），共获取4582个盐分有效样本。土壤盐分用电导率仪(上海雷磁DDSJ-308F)测试水土比为5:1的土壤浸提液的电导率。文中所用盐分值为2个重复的均值。图1中的32个长期监测点是根据改进时间稳定性方法确定的，将其用以验证时空地统计方法在监测点减少情况下获取区域盐分时空动态的能力，具体确定方法见文献 24。1.3 时空地统计时空克里金方法可用来评估土壤盐分的时空动态。通常普通时空克里金有以下 4 个步骤14：（1）时空经验半方差函数的计算；（2）用理论时空半方差拟合经验半方差；（3）用时空克里金做时空预测；（4）用留一交叉验证方法评估预测结果。各图1 隆胜研究区及盐分取样点位置图Fig.1 Longsheng study area and soil salinity samplinglocations18340卷干旱区研究步骤的详细计算方法如下：用时空经验半方差函数(hS,hT)描述土壤参数时空分布的相关属性13：(hS,hT)=12N(hS,hT)i=1N(hS,hT)Z(si+hS,ti+hT)-Z(si,ti)2（1）式中：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)，三者关系表示为：(hS,hT)=S(hS)+T(hT)+ST(hST)（2）式中：S(hS)，T(hT)和ST(hST)是各自独立的，通常由空间理论半方差函数模型如球状模型、高斯模型、指数模型来表示。在联合方差ST(hST)中常引入几何各项异性比来计算时空滞后距离hST：hST=h2S+(hT)2（3）该研究区土壤盐分的3个独立方差S(hS)，T(hT)和ST(hST)均用球状模型来拟合。球状模型的定义如下24-25：i(hi)=0,hi=0C0i+Ci3hi2ai-h3i2a3i,0ai(i=S/T/ST)（4）式中：C0为块金值；C为拱高；C0+C为基台值；a为相关长度(变程)；i表示空间、时间和时空变量。参照空间相关度的概念，用各部分独立半方差的块金值之和与基台值之和的比值(STD)来评价时空相关性。STD0.75表示时空相关性较弱，其余表示中等强度时空相关性26。克里金(Kriging)插值是对区域变量取值进行最优线性无偏估计的方法。时空克里金插值类似于空间克里金，一个栅格中的土壤盐分值Z*(s0，t0)由邻近的取样点盐分值(si，ti)进行估算27。Z*()s0,t0=i=1niZ()si,ti（5）式中：i为已知观测点对估值点的加权系数。根据无偏估计和方差最小两项要求可确定加权系数28，经转化可由以下矩阵方程求得权重系数：11(hS,hT)12(hS,hT)1n(hS,hT)121(hS,hT)22(hS,hT)2n(hS,hT)1n1(hS,hT)n2(hS,hT)nn(hS,hT)1111012n=10(hS,hT)20(hS,hT)n0(hS,hT)1（6）式中：ij(hS，hT)为空间滞后距离为hS和时间滞后距离为hT的两点(si，ti)、(sj，tj)的时空半方差值；为极小化处理时的拉格朗日乘数；n为搜索域内的观测数。留一交叉验证用来评估时空克里金插值效果。交叉验证时，验证点的观测值被移出数据集，用剩余的数据来预测该时空点的土壤盐分值12。然后，比较土壤盐分观测值和预测值的拟合优度来评估时空克里金模型。拟合优度评价指标用平均相对误差MRE、均方根误差RMSE和决定系数R2，具体计算方法见文献 24。本研究基于Matlab环境，采用遗传算法估计变异函数参数29，并完成试验半方差散点图、理论模型拟合曲面图的绘制；时空插值结果分布图在ArcGIS 10.3中完成。土壤盐分的经典统计分析利用Excel 2019软件完成。2结果与分析2.1 土壤盐分季节性变化特征由表1可知，不同土层土壤盐分的空间均值在0.250.35 dSm-1之间变化，平均处于非盐碱化水平1。土壤盐分的峰度变化范围是0.7118.52，偏度变化范围为1.054.18，数据右侧拖尾严重，所有取样时间各层土壤盐分均是非正态分布的。土壤盐分空间变异系数的变化范围是0.431.14，为中强变异。00.4 m越接近表层，土壤盐分受降雨、灌溉、蒸发等上边界作用越强烈，盐分随时间波动越剧烈。0.40.6 m土壤盐分均值不同观测期变化范围为0.310.32 dSm-1，变化幅度不超过5%，是盐分最为稳定的土层。2017年、2018年生育初期及生育末期00.6 m根系层土壤盐分分别为0.32 dSm-1、0.34dSm-1、0.29 dSm-1、0.30 dSm-1，根系层土壤生育期处于积盐状态。以2017年生育期末和2018年生育1842期孙贯芳等：时空克里金评估河套灌区土壤盐分时空格局期初根系层土壤盐分为例进行分析，秋浇和冻融作用使根系层土壤盐分降低了14.7%。这与孙贯芳等30在河套灌区秋浇灌溉洗盐试验的结果基本一致：秋浇灌黄河水180 mm后，次年春播前01.0 m土