含植被覆盖影响的石羊河流域土壤水分遥感估算及空间格局分析

杨志辉,赵 军,*,朱国锋,2,刘家良,郭文镈,黄治化,王炎强

1 西北师范大学地理与环境科学学院, 兰州 730070

2 中国科学院西北生态环境资源研究院,冰冻圈科学国家重点实验室, 兰州 730000

土壤水分(SM,soil moisture)是陆气系统中的变量[1- 2],是水文循环和水热平衡的组成部分[3],对地表植被生长也至关重要[4]。此外,SM对地表大气温度的演变[5]、降水的形成[6]和碳循环[7]有重要影响,对干旱事件的形成和发展也有重要的驱动作用[8]。因此土壤水分在气候、水文、生态和农业研究中起着重要的作用。

遥感技术的发展为大范围、及时、准确的SM遥感监测提供了可能[9],成为检测和监测SM的主要手段[10]。早期SM光学遥感估算方法主要有光谱反射率法、短波红外水分胁迫指数、植被状况指数和温度植被干燥指数等[11-13]。这些方法在估算SM时取得了较好的效果,但由于地表条件的差异,应用于不同研究区也存在一定的差异。表观热惯量法适用于裸露土壤,温度/植被指数法适用于植被覆盖度较大区域[14],但由于不同植被覆盖率的影响,降低了SM的估算精度。因此探索遥感估算SM的方法必须通过消除植被覆盖影响来实现。消除混合像素中植被反射率贡献的方法主要有：光谱混合分析法、光学植被覆盖度法和梯形特征空间法[15]。但在实际应用中光谱混合分析和光学植被覆盖度法计算处理复杂,估算条件较难满足。而梯形特征空间散点图是由地表温度/植被指数或反照率/植被覆盖度的散点图形成,对SM较为敏感且不依赖于地面观测数据,并可以与其他模型组合,所以梯形特征空间被广泛应用于估计植被蒸腾作用、土壤蒸发和土壤含水量[16]。

石羊河流域生态问题受广大学者关注[17]。该流域农业用水量比重较高,土壤水分的准确估算对农业灌溉用水的精确计算至关重要[18],因此,提高SM的估算精度、分析SM的空间格局是该流域目前亟待解决的问题之一。基于此本文引入散点图梯形特征空间来消除植被的影响,提高SM估算精度,进而分析空间格局,为土壤水分的有效利用提供依据。

1 数据与方法

1.1 研究区概况

石羊河流域位于甘肃河西走廊东部,祁连山北麓,介于101°22′—104°16′E,36°29′—39°27′N[19],属河西走廊三大内陆河流域之一,总面积4.16万km2(图1)。流域太阳辐射强,温差大,降水少,年降水量自南至北约为700—50mm、潜在蒸发量约为1200—2600mm,属于典型的大陆温带干旱气候,且气候特征具明显垂直分异。上游祁连山高寒半干旱湿润区,植被以高山和亚高山草甸为主,水土保持能力较强；中游走廊平原干旱区,是重要的灌溉农业区；下游温暖干旱区,主要为民勤绿洲和荒漠[20]。

图1 研究区位置概况Fig.1 Location of the study area

1.2 研究数据

本文采用空间分辨率500m,时间分辨率8d的MODIS-MYD09A1数据的1—7波段计算地表反照率,从2018年11月至2019年10月共46期数据,该数据能够满足生产活动对土壤湿度数据快速获取和及时更新的要求。由空间分辨率为250m,时间分辨率为16d的MODIS-MYD13A1 NDVI产品数据计算植被覆盖度(fractional vegetation cover,FVC),从2018年11月至2019年10月共24期。上述数据来源于NASA(https://ladsweb.modaps.eosdis.nasa.gov/)。

土地利用数据和土壤类型数据来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)。

于2019年8月下旬用土钻对石羊河流域0—10cm地表进行土壤采样,用手持GPS记录样点位置,样点分布如图1所示。将土壤样本装入铝盒中并现场称重,带回实验室放入105℃的烤箱中持续烘烤8h,然后将其放置在干燥器中冷却和称重,持续操作至两者差小于3mg,最后计算实际土壤含水量。实测SM数据用于构建土壤水分估算模型和精度验证。

1.3 研究方法

1.3.1地表反照率

采用Liang Shunlin在2001年提出的MODIS表面反照率模型来计算地表反照率[21],该模型在反演华北土壤水分时也得到了充分应用[22],公式如下：

a=0.16ρ1+0.29ρ2+0.243ρ3+0.116ρ4+0.112ρ5+0.08ρ7-0.0015

(1)

式中：a为全波段反照率,ρ1—ρ5、ρ7分别为MODIS相应波段的地物反射率。

1.3.2建立梯形特征空间

梯形由asd,avd,avw和asw四个点形成(图2)。其中asd和asw表示植被覆盖率为零的最大和最小地表反照率。点avd和点avw表示植被覆盖率为100%时最大和最小地表反照率。其中上边界称为“理论干燥边”,表示在给定的植被覆盖条件下土壤是非常干燥的。但在实际应用中,“理论干燥边”通常无法准确获得,取而代之的是“实际干边”,表示植被覆盖率变化时的实际表面反照率。“理论湿边”表示在一定植被覆盖条件下土壤水分是充足的。同样,在实际应用中使用“实际湿边”。位于“理论干边”和“理论湿边”之间的SM等值线便可以通过线性插值计算[23]。在SM等值线的基础上可以容易计算裸土反照率。本文根据像元二分模型和归一化植被指数计算研究区植被覆盖度[24]。

图2 表面反照率/植被覆盖梯形特征空间示意图Fig.2 Trapezoid space schematic diagram of surface albedo/vegetation coverage

1.3.3裸土反照率计算

研究区地表反照率主要取决于当地植被和土壤特征。因此,反照率是土壤和植被的加权总和[25],结合梯形特征空间得到裸土反照率计算公式,计算步骤详见文献[26]。

(2)

式中：asoil代表裸土反照率,a代表地表反照率。FVC表示植被覆盖度。kd和kw分别表示“实际干边”和“实际湿边”的斜率。ai,max和ai,min表示植被覆盖度较低时的实际地面最大和最小反照率。

1.3.4模型精度评估

计算SM估算值与实测土壤水分值的拟合度(R2),均方根误差(RMSE)、平均绝对百分比误差(MAPE)、不等式系数(TIC)、绝对误差(AE)和相对误差(RE)来评估土壤水分估算模型精度。此外,通过逐像元计算植被覆盖度和土壤水分的Pearson′s相关系数来评估植被覆盖对土壤水分空间分布的影响[27]。

1.3.5稳定性计算

标准差是衡量数据分布程度的方法,可用于评估时间序列上SM的稳定性[28]。公式如下：

(3)

1.3.6空间自相关

为分析整体尺度上SM在地理空间上是否存在集聚特征,采用全局Moran′sI(式(4))评估SM在空间上的自相关性。同时用局部空间相关 Moran′sI(式(5))系数进一步量化空间要素与周围要素之间的差异程度及空间分布特征[29-30]。

(4)

(5)

1.3.7地理探测器模型

地理探测器主要包含以下4种: 因子探测,探测每个影响因子对 SM空间分异的影响大小,q值越高则对SM的影响力越大。风险探测,探测因子对SM是否具有风险性,指示因子在不同等级内对SM分布的影响。交互作用探测,用于判断不同影响因子对SM 的交互作用。生态探测,用于比较两因子对SM空间分布的影响是否有显著差异。分析过程详见文献[31]。q值是衡量一个因子对因变量的解释力,表达式为：

(6)

2 结果与分析

2.1 生成梯形特征空间散点图

Toby Carlson的研究已表明：沙漠可视为地表裸露的干土,反照率几近于裸露的干土反照率；水体可视为地表裸露的湿土,反照率等于裸露的湿土反照率[32]。基于此并考虑了不同植被覆盖条件的地表反照率特征,在无覆盖、少覆盖、多覆盖和全覆盖四类植被覆盖度区域内划定样方(图3)。获得样方内对应格网点的植被覆盖率和地表反照率值,剔除异常值(FVC <0,或FVC> 1)后绘制反照率/植被覆盖梯形特征空间散点图(图4)。

图3 不同植被覆盖度下样方分布Fig.3 Sample distribution under different vegetation coverage

图4 植被覆盖度/表面反照率梯形特征空间散点图Fig.4 Trapezoidal feature space scatter plot of vegetation coverage / surface albedo

2.2 土壤湿度估算模型的建立

刘培君等研究表明SM与地表反照率或裸土反照率呈自然对数关系[33]。因此,根据实测SM值和对应样点的地表反照率或裸土反照率值通过回归分析建立估算模型和确定精度系数R2(表1),从R2可以看出,裸土反照率模型的准确性要明显高于地表反照率模型。

表1 地表反照率/裸土反照率的土壤水分回归模型Table 1 Sregression model of surface albedo / bare soil albedo

2.3 模型精度分析

实测数据与SM估算模型之间的误差分析(表2)表明：除01、02、11和12号点外,所有裸土反照率模型结果的AE和RE都小于地表反照率模型,RMSE、MAPE、TIC和R2值皆是裸土反照率模型小于地表反照率模型。此外,本文拟合度(R2=0.92)与同类研究(R2=0.57)和同流域研究(R2=0.79)比较精度都有所提高[34- 35]。结果表明裸土反照率模型能有效地消除植被覆盖对地表土壤光谱的影响,明显地提高了SM的估算精度。

表2 模型的精度分析与比较Table 2 Accuracy analysis and comparison of models

2.4 土壤湿度空间分布

参考WU Chengyong在青海对土壤水分的分级标准[36],结合研究区SM实际状况,对石羊河流域全年平均土壤水分分级(图5),结果显示流域内土壤水分整体较低,上游至下游土壤水分依次递减,且各级区域依次相邻分布。其中,SM大于25%的区域1254.53km2,占比最低为3.09%,主要分布在流域西南部祁连山山麓地带；SM为15%—25%的区域面积占比8.68%；这两级皆分布在海拔较高区域,中游仅红崖山水库区有小范围分布；SM为7%—15%的区域占比最高为33.91%,主要分布于凉州农灌区和民勤绿洲邻近水库部分；SM小于7%的区域则分布在民勤县荒漠区。

图5 土壤水分空间分布Fig.5 Spatial distribution of soil moisture

2.5 土壤水分稳定性

用2018年10月—2019年11月的46期SM数据,依据式(3)计算SM的稳定性(图6),反映土壤水分的年内变化。良好的稳定性代表SM几乎保持不变,受季节变化和外部影响较小；差的稳定性则变化较大,不利于SM的保持。图6表明：稳定性最好的区域主要分布于下游,因为该区域常年裸露且SM值低。稳定性较好的区域主要分布于中游和民勤绿洲,稳定性差的区域主要分布于上游,但占比最小,主要受林地植被覆盖随季节变化和积雪的影响。SM不同水平稳定性的面积占比如图6所示。

图6 土壤水分稳定性特征Fig.6 Characteristics of soil moisture stability

2.6 SM空间自相关特征

利用ArcGIS和Geoda软件分析研究区土壤水分的空间聚集特征,全局Moran′sI系数为0.88(Z-score=1852.94,P<0.01),表明石羊河流域SM存在强烈的全局空间自相关性； Moran′sI散点图如图7所示,样点主要分布在第一(高-高,正相关,占总散点41.95%)和第三象限(低-低,正相关,占总散点57.97%),表明研究区SM值呈现高-高和低-低聚集格局。

图7 土壤水分全局 Moran′s I 系数散点图Fig.7 Scatter plot of soil moisture global Moran′s I coefficient

此外,基于局域Moran′sI指数分析(P<0.05)得到局部空间自相关特征(图8)。SM聚类被分为四种类型：高-高(正空间自相关,占21.95%),主要分布于石羊河上游的肃南、天祝、武威、永昌等地；低-低(正空间自相关,占30.33%),主要集中分布于下游民勤荒漠区；低-高(负空间自相关,占0.04%)和高-低(负空间自相关,占0.01%)区域占比极低,且未出现聚集特征。中游区域则未表现出显著的空间自相关特征。

图8 土壤水分聚类分布图Fig.8 Soil moisture cluster distribution

2.7 空间分异分析

利用地理探测器模型对SM空间分异特征归因,按照不同因子对SM空间分异特征的解释能力排序:植被覆盖度(q=0.42,P<0.001)>土壤类型(q=0.37,P<0.001)>高程(q=0.22,P<0.001)>土地利用(q=0.09,P<0.001)。表明FVC是该地区SM空间分布的主导因素,土壤类型次之。风险探测分析(图9)和交互探测器分析(表3)表明：不同土地利用对SM的影响存在差异,其中林地SM均值最大,耕地次之,未利用地最小。不同土壤类型对SM的影响也存在较大差异。且SM往往不受单因子作用,而是多因子协同交互影响。表3得出,交互作用最强的因子组合为土壤类型和植被覆盖度(q=0.503),说明土壤类型的不同显著增加了FVC对 SM的解释力。因子组合皆产生了较高的q值,同时因子交互作用表现为较强的非线性增强交互作用,且两因子对SM空间分布的影响有显著差异。

图9 风险探测器结果Fig.9 Results of risk detector

表3 交互探测器和生态探测器结果Table 3 Results of interactive detector and ecological_detector

3 讨论

空间统计不同土壤类型的SM特征(表4)表明：石羊河流域内半淋溶土面积最小,主要分布在海拔较高的山谷和阴坡,其SM均值最高为24.28%。钙层土和干旱土依次与半淋溶土邻近分布,其SM均值也相对较高。漠土和初育土分布广泛,覆盖中下游广大区域,该土壤类型区内无河流分支,植被稀少,蒸散发强烈,故SM均值较低。半水成土主要分布在武威和民勤绿洲边缘区,处在初育土和人为土的过渡带,盐碱土零星分布在初育土中,SM值皆较低。而人为土分布在中游,主要由灌溉农田组成。高山土分布在流域内海拔最高的西北部高山区,该区域常年气温较低,林地植被郁闭度较高,土壤腐殖质层厚,蒸发量小,益于SM的保持。这种现象与其自身的理化性质、地形因素、局地气候和植被郁闭度等有关,这些因素对SM的有效保持和空间分布都产生了重要的影响[37]。

表4 石羊河流域主要土壤类型的土壤水分特征Table 4 Soil moisture characteristics of the main soil types in the Shiyang River Basin

为探讨不同土地利用类型的土壤水分特征,依据土地利用图(图10)空间统计各土地利用的SM状况(表5)。草地大部分处于中等以下水平,其中较干区域占29.79%,中等湿度区域占53.48%。耕地SM大部分处于中等水平,其中中等湿度区域占69.92%。结合图6,发现草地与耕地的SM稳定性较好。林地SM整体较高,中等湿度及以上部分占80.97%。未利用地大部分的SM较低,其中干土壤44.7%,较干土壤49.63%。不同的土地利用类型对SM有重要影响,林地在涵养土壤水分方面优于耕地和草地[38]。

图10 土地利用类型Fig.10 Types of land use

表5 不同土地利用的不同土壤水分百分比Table 5 Different soil moisture percentages for different land uses

土壤水分和植被覆盖度空间相关性分析(图11)表明：SM和FVC具有较强的相关性,植被覆盖度对土壤水分产生了重要的影响,凉州区和民勤县低海拔区表现为负向强相关(R=-0.57,P<0.01),其中草地占66.51%,耕地占25.96%。主要因为该区域降水稀少、蒸散发强烈,浅根系植被水分吸收以地表土壤水为主,不利于土壤水的保持,但耕地长期耕作土壤质地松软且有灌溉水源,因此耕地占比较小。FVC与SM存在负相关(R=-0.37,P<0.01)区域则主要分散在下游,该区域靠近巴丹吉林沙漠和腾格里沙漠,沙漠化效应影响较大。结论与柳菲等研究结论基本一致[39]。上游高海拔区域SM与FVC则为正向强相关(R=0.49,P<0.01),其中林地占80.16%,海拔较高区域的FVC与SM为正相关(R=0.13,P<0.01),其中林地占38.25%,因此林地对SM表现出了正向影响。

图11 土壤水分与植被覆盖相关性Fig.11 Correlation between soil moisture and vegetation cover

李浙华等[40]研究显示,土壤水分通常受到气候、水文、地形、植被覆盖度以及人类活动等因素的影响。但在小尺度上气象和水文要素基本一致,地形、植被和土地利用的变化则成为影响SM保持和变化的重要因素。植被覆盖度高的区域地表蒸散发较小,有利于SM保持。同时,不同高程也有影响作用,上游高海拔区域温度较低且多为林地,使地表蒸散发更小,且植被多为高大灌木,根系发达,水分吸收多源于深层土壤。而下游植被则多为草本和农田,冠层密度低且根系较浅,主要吸收地表土壤水。从而不同的高程和土地利用影响了植被对土壤水分的影响程度。

限于作者水平,高山地区冰川积雪对土壤水分遥感估算的影响有待于后续深入研究,以期进一步提高土壤水分的遥感估算精度。

4 结论

(1)通过地表反照率和植被覆盖度构建梯形特征空间计算裸土反照率,进而估算地表土壤水分,可明显地消除植被覆盖的影响,提高遥感土壤水分估算精度,为流域尺度的SM计算提供了方法参考。该模型适用于石羊河流域,结果符合研究区实际情况。

(2)石羊河流域土壤水分表现出西南高东北低的态势,上游至下游SM值逐渐降低,SM值大于15%的区域面积为11.77%,主要分布在上游高海拔区域,SM值7%—15%的区域面积为33.91%,主要分布在武威绿洲和民勤绿洲邻近水库部分,SM值小于7%的区域面积为54.32%,主要分布在下游民勤境内的荒漠区。

(3)石羊河流域SM空间自相关特征显著,高-高聚集主要分布于石羊河上游林地,占研究区21.95%,低-低聚集主要集中分布于下游民勤的荒漠,占研究区30.33%。上游林地区SM呈高高聚集,下游荒漠区SM呈低低聚集。土壤水分和植被覆盖度存在显著相关性且不同土地利用影响程度存在差异,上游林地与SM呈显著正相关(R=0.49,P<0.01),中游草地与SM呈显著负相关(R=-0.57,P<0.01)。

(4)石羊河流域植被覆盖因素对SM的空间差异解释力最强,各因子对SM空间分异的解释能力存在显著差异,且因子间交互作用均增强了其对SM空间分异的解释力。

(5)石羊河流域土壤水分稳定性整体较好,中下游皆表现为稳定。不同土地利用类型的SM差异较大,未利用地大部分SM小于7%,大部分草地和耕地SM处于中等水平,SM值为7%—15%,林地SM值大于25%,水平整体较高。