考虑到半干旱区在全球变化和水文研究的重要意义以及遥感数据在该区域比较完备(全年有值天数达40%左右,远高于湿润地区)等因素,故选取半干旱区的老哈河流域作为研究区,其地理范围位于北纬41°—43°、东经117°—120°,基于遥感土壤湿度数据估算根系层土壤湿度时空分布。这里以兴隆坡水文站控制的集水区作为老哈河流域研究区,位于西辽河上游,流域地跨河北、辽宁两省和内蒙古自治区,海拔427～2 054 m,面积约为2.1万km²。

研究区采用的土壤湿度数据是1993—2006年的4—9月共84个月的农气站土壤湿度观测资料和CCI土壤湿度产品。农气站实测的SM数据来源于中国农作物气象数据网的中国农作物生长发育和农田土壤湿度旬值数据集,包括10 cm,20 cm,50 cm,70 cm和100 cm多个根区剖面不同深度土壤相对湿度的逐旬观测值,原始数据是相对土壤湿度(%)。本研究采用10 cm深度处的相对饱和值作为表层SM,而在10 cm以下根部剖面上相对饱和度的均值作为根系层SM。考虑到数据的有效性,选取1993—2006年的逐年4月、5月、6月、7月、8月、9月实测的逐旬数据,分别求取表层和根系层不同作物逐月各旬的均值作为月平均,再剔除个别缺测月份最终作为各农气站多年逐月实测的相对土壤湿度数据。由于SMAR模型采用的简化土壤水分平衡方程适用于干旱半干旱区,因此本文根据中国300～500 mm年等降雨量线选择了82个农气站点的数据(图1)。

图1 年降雨量300～500 mm农气站分布图

本文采用的NASA GLDAS V2.0数据集是GLDAS中NOAH陆面模式计算的蒸散发产品,数据来源:http:∥disc.sci.gsfc.nasa.gov/hydrology/data-holdings,数据格式是NetCDF,空间分辨率是0.25°,时间分辨率是月。

归一化植被指数NDVI数据采用的是GIMMS NDVI3 g V1数据集,来源于:https:∥ecocast.arc.nasa.gov/data/pub/gimms/3 g.v1/,原始数据格式为NetCDF,每年有两个nc4文件,每个nc4文件包含6个月NDVI数据,时间分辨率为15 d,空间分辨率为0.083°,共12景。本研究通过格式转换、最大值合成等预处理过程,使得每一个像元为该月最大NDVI值,从而获得研究区逐年逐月的NDVI数据集。

本文搜集农气站以及研究区逐像元的土壤属性数据来源于《中国土种志》^[17-19]。中国土种志以大区为单位,分别介绍了各大区主要土种在土壤分类系统中的归属、分布、面积、主要性状、典型剖面、生产性能和理化分析数据。本研究主要整理了82个农气站对应的典型剖面的土壤属性组成,包括粗砂(2～0.2 mm)、细砂(0.2～0.02 mm)、粉砂(0.02～0.002 mm)、黏粒(<0.002 mm)以及各深度对应的质地名称,将表层深度以下的各属性均值作为站点根系层土壤的各属性数据。

由于遥感获取的CCI土壤湿度产品和农气站实测的地表土壤湿度数据存在系统误差,导致这两组数据存在不一致性,因此为了使遥感数据与农气站实测数据更好地匹配,采用累积概率分布函数(CDF)^[16]重调整方法,基于实测土壤湿度数据对遥感数据进行标准化,使其具有相同的取值范围和累积概率分布特征。CDF重调整过程是以像元点为基础进行时间序列统计并调整。在实践中,农气站的分布数量有限,且不均匀,因此需要通过对站点数据进行空间插值获取所研究区域SM的整体时空分布状况。本文选择老哈河流域周围的赤峰等8个站点的1993—2006年4—9月的有效数据进行克里金空间插值^[16],以此获得老哈河流域表层实测SM的时空分布,再在此基础上利用CDF匹配技术将CCI土壤湿度产品按照与实测数据相同的时间序列进行重调整,建立研究区逐像元的CCI土壤湿度产品的改进模型。

CDF匹配技术^[20]是基于分段线性匹配,研究中,调整后的土壤湿度表达式如下:

CDF_site(x')=CDF_CCI(x)(1)

式中:CDF_site和CDF_CCI分别表示站点和CCI土壤湿度的累积分布函数; x表示调整之前的CCI土壤湿度。

Manfreda等^[1]假定土壤由两层组成,第一层是几厘米深的表层(合理的假设是卫星传感器观测深度可以代表5—10 cm的表层动态),而第二层位于第一层之下,延伸到相当于植被的生根深度。两层土壤最重要的水分交换是入渗,横向流动和毛细管上升等过程则可以忽略不计^[14]。Manfreda等^[1]定义了一个土壤水分平衡方程,它是由半干旱环境中一个简化的土壤水分平衡方程推导出来的,在这个方程中入渗期不是表示为降雨的函数,而是表层中土壤含水量的函数,分别用“1”和“2”区分上下层的变量和参数。只有当土壤水分超过田间持水量时,表层的水分通量才能被认为是显著的。他们提出了描述从顶层到下层的瞬时入渗通量的方程:

n₁Z_r1y(t)=n₁Z_r1y(s₁(t),t)=n₁Z_r1{(s₁(t)-s_c1),s₁(t)s_c1

0,s₁(t)<s_c1(2)

式中: y(t)[-]表示下层土壤饱和渗透分数; n₁[-]表示第一层土壤孔隙度; Z_r1表示第一层土层深度; s₁((θ₁)/(n₁))[-]表示第一层的相对饱和度(它是由第一层土壤体积含水量与土壤孔隙度的比值决定的); s_c1[-]表示第一层田间含水量相对饱和值。此公式意味着当相对饱和度高于田间含水量的任何值时,土壤的渗透性是无限大的,该模型不考虑下层的饱和效应。

第二层和深层土壤水分平衡受入渗和水土流失两个因素控制,水土流失主要又包括蒸发蒸腾损失和渗透。定义土壤第二层的“有效”土壤相对饱和度为:

x₂=(s₂-s_w2)/(1-s_w2)(3)

土壤贮水量为:

ω₀=(1-s_w2)n₂Z_r2(4)

则土壤水分平衡可以用式(5)来表述:

(1-s_w2)n₂Z_r2(dx₂(t))/(dt)=n₁Z_r1y(t)-V₂x₂(t)(5)

式中: s₂[-]表示第二层的相对饱和度; s_w2[-]表示凋萎点的相对饱和度; n₂[-]表示第二层的土壤孔隙度; Z_r2[L]表示第二层的土层深度; V₂[LT^-1]表示蒸散和渗透损失的土壤水分损失系数; x₂[-]表示第二层土壤的“有效”相对饱和度。此方程不考虑土壤水分高值时土壤流失函数的高度非线性,这种简化,再加上渗透不考虑饱和效应,意味该方法在潮湿环境中有一定的局限性,更适用于干旱半干旱环境。

上述方程可以用归一化系数和简化为:

a=(V₂)/((1-s_w2)n₂Z_r2)

b=(n₁Z_r1)/((1-s_w2)n₂Z_r2)(6)

“扩散系数”是土壤表层总可能储水量(n₁Z_r1)与总根区储水量((1-s_w2)n₂Z_r2)之比确定,“水分流失系数”控制根系层水分损失,是蒸散发ET与排水损失之和与总根区储水量之比确定,这些参数的取值可直接与两土层的深度和水土流失系数之比有关。因此,土壤水分平衡方程即为式(7):

(dx₂(t))/(dt)=by(t)-ax₂(t)(7)

假设相对饱和度x₂(t)初始条件等于0,可以推导出这种线性微分方程的解析解:

x₂(t)=∫^t₀be^a(ω-t)y(ω)dω(8)

在实际具体应用中,可能还需要离散的形式见公式(9):

x₂(t_j)=∑^j_i=0be^a(t_i-t_j)y(t_i)Δt(9)

扩展式中,假设Δt=(t_j-t_(j-1)),根据表层土壤湿度的时间序列可以推导出下层土壤湿度的公式(10):

x₂(t_j)=x₂(t_j-1)e^{-a(t_j-t_j-1)}+by(t_j)(t_j-t_j-1)(10)

也可以写成s₂的函数:

s₂(t_j)=[s_w2+(s₂(t_j-1)-s_w2)e^{-a(t_j-t_j-1)}]+[(1-s_w2)by(t_j)(t_j-t_j-1)](11)

式中:左半部分计算给定的时间步长中相对水分损失,水分损失函数随s₂线性减小,直至达到s_w2,此时水分损失为零,蒸散和排水损失假定在凋萎点值或以下可忽略不计,右半部分表示从表面到根区扩散的水分含量。Manfreda等^[1]提出的方法是一种代表了两个状态变量之间土壤水分分析关系式,需要确定的4个参数分别是:s_w2,a,b,s_c1,所有参数均可从土壤质地分数、土层深度和土壤水分损失等估计。

为了分别得到最优的模型参数组合,本文使用MATLAB的遗传算法工具箱GUI进行模型参数的优化估计^[14,21]。MATLAB遗传算法工具,即图形用户界面GUI,它使我们可以使用遗传算法而不用命令行的方式。遗传算法^[22]是一种基于自然选择、生物进化过程来求解问题的方法。在每一步中,遗传算法随机地从当前种群中选择若干个体作为父辈,并且使用它们产生下一代的子种群。在连续若干代后,种群朝着优化解的方向进化,最终得到最优解组合。本文以各站点表层SM数据模拟得到的根系层相对SM与农气站实际测量得到的根系层相对SM均值之间的RMSE最小为优化准则进行优化处理,假设每个参数的上、下边界分别相当于最大和最小可能值(例如,0≤s_c1≤1)。参数估计时,第二层土壤相对饱和度s₂初始值设置为各农气站实测的逐旬田间持水量相对饱和值的均值。应该注意的是,SMAR模型可能产生高于1的值,这些值需自动设置为1。

CDF重调整过程是以像元点为基础进行时间序列统计并调整,本文以研究区赤峰站所在像元进行示例分析,CDF重调整过程为:

(1)分别绘制农气站实测表层和遥感SM数据的CDF曲线,以0%,5%,10%,20%,30%,40%,50%,60%,70%,80%,90%,95%,100%这13个节点得到曲线上对应点的坐标,将曲线划分成12段,见图2A、B)所示。

(2)根据上一步得到各节点的CDF值建立CCI土壤湿度与农气站表层实测SM的分段线性回归曲线,见图2C所示。对每个部分进行分段线性回归,从而得到12组线性回归方程。

(3)利用得到的分段线性回归方程,进行CDF匹配调整,从而得到农气站实测土壤湿度、CCI土壤湿度和调整后的土壤湿度(图2D),图中实线表示农气站实测土壤湿度,粗虚线表示原始CCI土壤湿度,细虚线表示调整后的土壤湿度数据。

图2 实测土壤湿度的CDF曲线(A)、CCI土壤湿度的CDF曲线(B)、CCI与实测土壤湿度分段线性回归曲线(C)、CDF重调整前后比较(D)

以内蒙古赤峰站所在像元为例,已知1993—2006年逐年4—9月的站点实测土壤湿度值及CCI遥感土壤湿度,将实测数据转化成以体积含水量为单位,进行CDF重调整,从而获得调整后的土壤湿度数据时间序列曲线(图3),从图3可以看出,原始CCI SM整体上小于实测的SM,但是两套数据的趋势具有良好的一致性,调整后的数据普遍向上抬升,在数值和形态表现上都更加靠近基准的SM数据。对该站点所在像元CDF重调整前后的CCI SM与实测数据之间的误差进行统计分析(表1),可以看出CCI SM经重调整后各项误差指数都得到降低,尤其是MAE和RMSE减少12%以上,表明调整后的数据更加接近于实测的SM数据。综上,赤峰站所在像元的实测土壤数据经过CDF重调整后精度得到明显改善,整个研究区像元调整方法类似,误差也得到减小,精度得到提高。

表1 赤峰站所在像元CDF调整前后误差分析

本文首先根据中国农作物气象数据网的中国农作物生长发育和农田土壤湿度旬值数据集整理了分布在中国年降雨量线300～500 mm内的82个站点不同剖面深度实测相对SM,将农气站表层10 cm的相对SM的时间序列作为SMAR模型的输入数据,采用MATLAB的遗传算法优化由表层相对SM模拟作为输入得到的根系层相对SM与其实测的10 cm以下的根系层相对SM均值之间的RMSE,从而得到每个农气站SMAR模型最优的参数组合,为后续SMAR参数的统计分析奠定基础。该遗传算法能够在大多数站点位置找到每个参数的近似全局最小值。这里以老哈河流域周围的包括赤峰等6个站点为例运用SMAR模型估算站点尺度根系层的相对土壤湿度。

图3 赤峰站所在像元CDF匹配调整前后的土壤湿度时间序列曲线

图4和表2给出了SMAR模型在研究区及周边几个农气站14 a的初步应用,可以看到该方法能很好地捕捉下层土壤水分的动态变化,与实测根系层的相对SM趋势变化具有良好的一致性,并接近于每个站点的观察结果。RMSE和R分别表示农气站实测10 cm以下根系层的相对SM均值与SMAR模型模拟的根系层的相关关系。根据MATLAB遗传算法^[14]对站点的参数进行校准,在均方根误差方面表现较好,而在相关性方面各站点表现好坏存在差异,R最低的奈曼旗站为0.471,最高的阜新站为0.895。

表2中的参数值和误差值都与站点实测值存在密切关系,产生误差的一部分原因是SMAR模型估算根系层相对土壤湿度是具有时间序列的,但在实际运行时由于部分月份可能因为数据缺测而被剔除,造成前后时间间断,时间序列的中断也会对后续模型估算产生影响。从模型运行来看,根区田间持水量的初始值对于模型的预测影响是可以忽略不计的。另外各站点气候条件不同或区域的特定空间异质性也导致了它们SMAR模型参数值的不同。

表2 6个农气站点SMAR模型的参数组合

综上,由老哈河流域及周边实测站点的根系层土壤湿度估算可以看出,SMAR模型估算半干旱的站点根系层SM具有良好的效果,为区域根系层SM的估算奠定了基础。

考虑到SMAR模型参数的物理意义与土壤和植被的关系,因此可以将土壤性质、NDVI和ET^[3,6]用于预测SMAR参数(表3),所构建的回归关系与欧洲相关学者^[15]运用土壤性质作为变量建立的pedo—转换函数具有类似的预测能力。

本文先根据各农气站的质地编号结合《中国土种志》^[17-19]整理了各站点土壤不同剖面深度下的机械组成(包括粗砂、细砂、粉砂、黏土含量),剔除未知土壤类型的站点以及参数异常的站点后剩余38个站点用于多元回归分析。用这些农气站点土壤物理性质、NDVI和年平均ET数据(这里提取的是逐年4—9月的均值为年平均)作为变量,建立了估计SMAR参数的回归方程,并通过了SMAR参数与土壤性质之间(p<0.05)的双尾t检验。

站点根系层凋萎点水平与土壤性质的相关性最强,与根区细砂含量呈明显的负相关关系(双尾t检验,p<0.001)(图5A),从模型(表3)中也可以看出细砂含量对s_w2有明显的作用。年均ET是SMAR各参数的一个重要且具有统计学意义的变量,因为植物水分吸收和表面蒸发会影响近地表和根系层水分状况,与各参数都密切相关。其中扩散系数b与年平均ET呈显著的正线性相关(双尾t检验,p=0.008),这表明在b的回归中,应加入ET,通过近地表扩散率与表层田间持水量的联系,也应将其考虑到田间持水量s_c1的回归中。水分流失系数a反映的是蒸发蒸腾和排水损失之和,从图5C中可以看出与年平均ET呈正相关关系(双尾t检验,p=0.038),图5D可以看出与上层粉砂的含量也呈正相关(双尾t检验,p=0.04),建立的模型见表3。黏土含量和孔隙度与其他变量相比影响不大,在模型中也没有表现出来。NDVI能够检测植被的生长状态,可以反映植物冠层的背景影响,如土壤等,因此也可以用来估计SMAR模型的相关参数。从表3中可以看出扩散系数与NDVI呈负相关关系,这里的NDVI采取的是各站点所有时间序列的平均值。综上,表3所示的多元线性回归函数均能良好地估算SMAR模型的4个参数,并都通过了双尾t,p<0.05的显著性检验,这为区域逐像元模型的应用提供了支撑。

图4 6个站点实测的表层、根系层土壤湿度和SMAR模拟得到的土壤湿度时间序列曲线

表3 SMAR模型各参数的多元线性回归函数

由于农气站点数量与规模的限制以及遥感SM只能反映表层的土壤水分状况,因此本文将遥感SM产品与SMAR模型结合起来,进行流域空间内根系层土壤水分时空分布的估算。根据CCI遥感土壤湿度产品的空间分辨率0.25°×0.25°,提取覆盖老哈河流域的所有像元,共77个。根据《中国土种志》逐像元整理出对应的质地类型和土壤机械组成数据,并对应提取出各像元NDVI和年均ET数据,然后按照上文在农气站所使用的SMAR模型参数回归的方法进行参数回归。将重调整后的逐像元遥感SM作为SMAR模型的最终输入数据,再结合回归分析得到的模型参数,最终可以估算出覆盖流域范围的根系层相对土壤湿度的分布状况。以老哈河流域2006年4—9月为例,模拟出区域根系层相对土壤水分状况(图6),从2006年的模拟图中可以看出当遥感土壤湿度作为SMAR输入数据时,能够较好地模拟出根系层的相对土壤水分状态。

图5 SMAR模型参数与变量的相关关系散点图

图6 老哈河流域遥感数据作为输入模拟得到的2006年4—9月根系层土壤水分状况

为了进一步分析遥感SM数据作为SMAR模型输入的结果,分别以表层实测和重调整后的遥感数据作为输入,对覆盖老哈河流域的矩形栅格区域进行逐像元SMAR模拟,对两种结果进行相关关系和误差分析(图7)。可以看出,两者的相关系数总体较好,R平均值为0.692,经过SPSS双变量相关分析得到,97%的像元通过p<0.05的双尾t显著性检验。从4种误差分析来看,误差总体均小于0.1,其中平均绝对误差MAE表现最好,其余3种相差不大。总体来说,使用遥感数据与SMAR模型结合,能够较好地模拟出区域根系层的土壤湿度,为更大尺度预测根系层土壤水分提供了可行性参考,并为干旱预报和其他水文模拟业务提供了参考。

图7 老哈河流域逐像元SAMR模型估算的根系层相对土壤湿度相关分析