基于TVDI的甘肃黑方台土壤湿度分布及影响因素分析

叶伟林,周自强,石三娥,康丽娟,王得楷

(1.甘肃省科学院地质自然灾害防治研究所,甘肃 兰州 730000;2.兰州大学资源环境学院,甘肃 兰州 730000;3.成都理工大学地质灾害防治与地质环境保护国家重点实验室,四川 成都 610059)

甘肃省黑方台属于大陆性干旱气候,降水量少而蒸发量大,土壤含水量低。但由于20世纪60年代人口迁入,对该区进行了大量的农业灌溉,导致干旱的区域土壤含水量增高。同时长期灌溉使得地下水位上升,台塬边黄土底部逐渐饱水软化,引起土体抗剪强度下降,导致黑方台稳定性降低,致使其边缘成为滑坡的易发区[1-2]。据统计70年代以来黑方台发生大小滑坡约150次,造成不同程度的人员伤亡和经济损失,引起了社会的广泛关注。此外,土壤中的水分是滑坡的主要诱因[3],它是地表层和大气能量交换过程中的重要参与者,决定着植被的蒸发以及光合作用,并且对陆地表面蒸发、水的运移以及碳循环等有很强的控制作用,对于研究预测区域干湿情况有着举足轻重的作用[4-5]。因此,监测黑方台土壤湿度分布及其变化对分析灌溉等因素对研究区土壤湿度的影响以及灌溉对滑坡的影响具有一定的意义。

土壤湿度的传统获取方法主要包括质量法、中子仪法在内的田间实测法[6],其采用单点进行实地监测,虽然可以直接准确地测量土壤剖面的含水量,但是采样点数目太多,需要大量的人力物力,不仅费时而且成本高,很难大范围、高效率地获取土壤水分。此外,区域在土壤、地形、植被覆盖上的空间差异使单点的代表性差,只能够反映出采样点局部区域的土壤湿度。因此,田间实测法在大范围监测土壤湿度的过程中存在很大的局限性[7-9]。遥感手段获取土壤湿度相对于传统的手段具有范围广、速度快、成本低、可多期监测等特点,分为微波遥感和光学遥感,包括热惯量法[7-9]、蒸散量法[10]、作物缺水指数法[10-12]、温度植被干旱指数(TVDI,temperature vegetation dryness index)法等[13-15]。TVDI通过构建地表温度(LST,land surface temperature)和归一化植被指数(NDVI,normalized difference vegetation index)的对应关系来反映土壤湿度状况。Price[16]于1990年发现LST与NDVI呈负相关关系,以它们为横纵轴构成的散点区域呈三角形,在此基础上,Sandholt等[17]通过总结前人研究,首次提出了TVDI来监测表层土壤水分状况,获得了较好的效果;赵杰鹏等[18]通过改进TVDI反演模型,较精确地估算了新疆地区的土壤水分状况;王纯枝等[19]利用MODIS产品数据NDVI和Ts构建了特征空间计算得到TVDI,经验证TVDI能准确的反映黄淮海平原土壤10～20 cm处的土壤水分状况;王行汉等[20]利用增强温度植被指数(ETVDI,enhanced temperature vegetation dryness index)和地表温度反演TVDI,克服了TVDI在高植被覆盖区植被指数饱和的缺陷,为南部高植被覆盖区的土壤水分监测提供了一种新的方法。

众多学者采用TVDI方法对不同区域土壤湿度进行了大量的研究,得到了很好的效果。黑方台为低植被区域,用TVDI来研究土壤湿度效果更好。基于此,利用2019年11月遥感数据反演TVDI分布状况,之后通过相同日期的采样数据建立TVDI与不同深度处土壤湿度的关系,分析TVDI与不同土壤深度关系的可靠性和准确性。在此基础上,分析黑方台土壤湿度空间分布以及动态变化,探讨灌溉等因素对黑方台土壤湿度的影响以及灌溉对滑坡的影响。

1 研究地区与研究方法

1.1 研究区概况

甘肃黑方台地处西北黄土高原,为黄河Ⅳ级阶地台塬,海拔1 700 m左右,由黑台和方台组成,属于温带大陆性气候,全年降水量少,昼夜温差大,日照时间长,四季分明。多年平均降水量为287.6 mm,年蒸发量达1 600 mm,最大降水量为431.9 mm,最小降水量为178.8 mm,年蒸发量约为年降水量的5.4倍[21-22]。20世纪60年代,由于刘家峡与盐锅峡水库的修建,大量的人员移居在此,为解决农业缺水问题,自1960年开始,大量黄河水被抽至台塬进行大水漫灌,长年灌溉使黄土底部形成了数十米厚的地下水位。灌溉水的入渗使得土体含水量增加,饱和度增大,引起土体抗剪强度下降,从而导致黑方台稳定性降低,容易发生滑坡[22]。

1.2 数据来源

(1) 遥感数据 选用2017—2019年间26景LandsatOLI/TIRS遥感影像(影像获取于https://www.usgs.gov/)。遥感影像预处理具体如下:首先利用Envi5.3所提供的辐射定标工具(radiometric calibration)对Band10、Band4、Band5进行辐射定标,将原始的DN值转化为辐射亮度值,后利用FLAASH模型进行大气校正,减少大气散射、吸收、反射引起的误差,得到真实的地表反照率,再对数据进行裁剪等处理。

(2) 野外采样数据 为了精确提取土壤湿度,在遥感影像获取当天进行了土壤数据采集,土壤湿度在该时间段内保持在稳定状态。采样点分布如图1所示,12个采样点分别分布在不同的土地利用类型和典型的滑坡体上。每个采样点分别在10 cm、20 cm、30 cm、40 cm、50 cm的埋深处采样,最终采集到5个深度的土壤质量水。

图1 研究区概况及采样点分布Fig.1 General data and distribution of sampling points of study area

(3)气温降水数据和其他数据 从中国气象数据网(https://data.cma.cn/)获取了永靖站点2017—2018年月平均气温、月降水、月日照时数、月相对湿度等气象数据;由于永靖站缺少月蒸发量数据,其经纬度为(103°18′E,35°58′N),而榆中县、民和县经纬度分别为(104°09′E,35°52′N)、(102°5′E,36°2′N),发现永靖县恰处于中间且离它们最近,因此利用榆中县和民和县的月蒸发量均值来代替永靖县;每日的灌溉量数据来自于永靖县盐锅峡镇水管所。

1.3 研究方法

(1) 归一化植被指数 基于地表植被对红光波段和近红外波段不同的吸收和反射能力,通过比值运算可以有效地反映土壤背景信息和地表植被覆盖状况。生长较茂盛的植被在近红外和红外波段内表现出强度的反射差异,而生长不好的植被反射差异较少[23]。植被的生长状况以及分布范围能被NDVI有效提取,因此NDVI被广泛应用于自然环境中植被信息的提取,其反演公式为

(1)

其中:ρ1、ρ2分别为Landsat8近红外和红外波段的反射率。

(2) 地表温度 应用Landsat热红外波段反演地表温度主要有3种方法:大气校正法、单窗算法、单通道算法。基于大气校正法,利用Landsat8 TIRS反演地表温度[24-25],反演公式为

Lλ=[εB(TS)+(1-ε)L↓]τ+L↑,

(2)

(3)

(4)

其中:Lλ是卫星传感器接收到的热红外辐射亮度值;L↑是大气向上的辐射亮度;L↓是大气向下的辐射亮度;ε是地表比辐射率;B(TS)是黑体的辐射亮度;TS是地表真实的温度;对于TIRS Band10,K1=774.89 W/(m2×μm×sr),K2=1 321.08 W/(m2×μm×sr)。

(3) 温度植被干旱指数 Price[16]发现,当研究区植被覆盖满足从裸地到完全植被覆盖时,植被指数和LST构成的散点图呈现出三角形。Sandholt 等[17]基于NDVI和LST建立了LSTi-NDVI三角形特征空间,提出了温度植被干旱指数的概念TVDI,并证明地表温度LSTi与NDVI存在明显的负相关关系,任意一个NDVI值对应唯一一组LSTmax(干边)和LSTmin(湿边),因而LSTi-NDVI特征空间的拟合线斜率可以反映区域土壤湿度的情况(见图2)。TVDI计算公式为

图2 LSTi-NDVI的特征空间Fig.2 The feature space of LSTi-NDVI

(5)

LSTmin=a1-b1NDVI,
LSTmax=a2-b2NDVI,

(6)

(7)

其中:LSTmax(干边)和LSTmin(湿边)分别表示当NDVI等于某一特定值时,地表温度的最大值和最小值,即特征空间的干边和湿边;a1、b1、a2、b2是干边、湿边拟合系数;LSTi表示任一像元地表温度。

(4) 偏差分析 偏差为每月TVDI值距多月平均TVDI的距离,可表示某一时段内TVDI偏离多月TVDI均值的程度,其中正值表示高于多年平均水平,反之则表示低于多年平均水平[26],计算公式为

Di=TVDIi-TVDIm,

(8)

其中:Di为第i月偏差值;TVDIi为第i月TVDI值;TVDIm为多月平均TVDI;i为月份。

(5) 线性倾向率 利用线性倾向率分析2017—2019年TVDI的月际变化趋势,若为正值则表示TVDI呈现增加趋势,否则表示下降趋势。基于像元的线性倾向率计算公式为

(9)

其中:Bslope为线性倾向值;TVDIi为第i月的TVDI值;i为月份变量(i=1,2,3,…,26),n=26。当Bslope>0时,表示随时间增加TVDI呈上升趋势,反之,TVDI呈下降趋势,Bslope的大小表示TVDI上升或下降的速率。

(6) 相关分析 为了研究黑方台土壤湿度状况与自然和人为影响因素的线性关系,采用卡尔·皮尔森(KarlPearson)在1880年提出的Pearson相关系数,该系数能检验出变量因子间相关性的强弱和方向,已广泛应用于各学科领域,其计算公式为

(10)

其中:rxy为相关系数,值域范围是(-1,1),值大于零表示正相关,值小于零表示负相关,其绝对值越大表示相关性越显著;xi为第i月的平均TVDI值;yi为当月相关变量的均值。

2 结果与分析

2.1 土壤湿度监测结果检验

国内学者比较了各温度植被干旱指数与土壤湿度之间的关系,结果表明TVDI相对于其他指数能更有效反映土壤湿度变化[27-28]。为了验证研究区TVDI与表层土壤湿度的相关性,选取12个野外采样点,得到土壤不同深度60个实测土壤湿度数据,将其分别与样点处TVDI进行线性回归(见图3)。发现10 cm土壤深度处(见图3(a)),P=0.002,R2=0.63;20 cm土壤深度处(见图3(b)),P=0,R2=0.79,相关性最好;30 cm土壤深度处(见图3(c)),P=0.001,R2=0.68,相关性次之;在40 cm(见图3(d))、50 cm(见图3(e))处随着土壤深度的加深,相关性呈减小趋势,P值分别为0.007、0.028。结果表明研究区TVDI和实测土壤湿度在20 cm、30 cm土壤深度处呈较强的负相关关系,而在40 cm、50 cm处随土壤深度逐渐加深,TVDI指数和土壤湿度的相关性减小。

图3 土壤湿度和TVDI线性拟合Fig.3 Linear fit of soil moisture and TVDI

2.2 土壤湿度空间分布特征

为了分析黑方台土壤湿度的空间分布状况,对2017—2019年26个月的TVDI均值进行统计并按自然断点法(natural breaks)将其分为5个等级(见图4(a)),1～5等级月均TVDI的范围值分别为0.05～0.32、0.32～0.50、0.50～0.61、0.61～0.74、0.74～0.99,并对各等级的面积比进行了统计分析(见图4(b)),其中等级越高代表土壤越干旱,土壤湿度越低。黑方台3等级占研究区面积最大,为32.44%,说明对研究区而言,整体土壤湿度较高,干旱程度较低。等级4和等级5的区域主要分布在黑方台台塬边周围,占总区域的35.79%,这些区域多为裸地,没有农业灌溉,加之干旱的气候导致该区域在20 cm、30 cm处土壤湿度非常低。台塬内部土壤湿度以2、3等级为主,占总区域的58.5%,且4等级呈零星状分布其中。由于靠近黄河,离水源近,台塬内部以农用水浇地为主,包括耕地、林地等,因此土壤湿度较台塬边缘高。综上说明黑方台土壤含水量呈现台塬内高,台塬边低的空间分布格局,这主要是受人类活动农业灌溉的影响。

图4 月均TVDI的空间分布和各等级面积占比Fig.4 The spatial distribution of monthly TVDI mean and the area ratio of each grade

2.3 土壤湿度空间变化特征

(1) 时间变化特征 从黑方台各月TVDI均值来看(见图5),干旱月份较少,土壤湿度较高。TVDI值呈现不规则的波动变化趋势,峰值出现在2018年7月和8月,分别为0.67和0.65,超过平均值0.57;低值区出现在2017年11月、12月,2019年1月为0.47,低于平均值,因为此时研究区有一定量的灌溉且蒸发量较少。各月TVDI均值在时间尺度上呈微小的下降趋势,表明土壤湿度在该时间尺度上呈微小的上升趋势。应用偏差分析得到2017年3月—2019年11月各月TVDI偏离多月平均水平的程度,TVDI的偏差值表现的特征为在每年3—8月值为正,表现为增加的趋势,而在每年的10月至次年1月值为负,呈减小的趋势。据黑方台灌溉调查资料显示,每年春末、秋初及整个夏季黑方台灌溉量增大,土壤湿度会随着灌溉量的增大而增大,但由于研究区随温度升高,蒸发量急剧增大,导致这些月份土壤湿度反而较小。而在冬季等气温较冷的月份,蒸发量比夏季小,加之有少量的灌溉,导致研究区台塬土壤湿度较高。

图5 2017年3月—2019年11月 TVDI 月变化Fig.5 Monthly variations of TVDI during March 2017 to November 2019

(2) 空间变化特征 为了进一步探讨研究区土壤湿度的空间分异特征,采用最小二乘法,逐像元对研究区年均TVDI进行回归分析,结果见图6。由图6看出,整个黑方台土壤湿度表现出明显的空间分异格局,其中2017年3月—2019年11月TVDI上升的区域主要分布在黑方台的西北区域,说明该时间段内黑方台的土壤湿度呈减小的趋势。在台塬内部发现,黑台中部的小片区域内,变化趋势小于零,说明在该时间段内TVDI越来越小,土壤湿度在变大。此外,台塬上的其他区域变化率不太显著,说明台塬上大部分区域在2017年3月—2019年11月土壤湿度没有发生明显变化。这主要是由于黑方台区域土壤湿度在长时间尺度上受人类定时灌溉活动的影响,其台塬上TVDI基本保持不变。

图6 2017—2019年TVDI线性倾向率Fig.6 The tendency rate of TVDI during 2017-2019

2.4 土壤湿度影响因素

黑方台年降水少而蒸发大,农业用地均为水浇地,鉴于黑方台以上环境条件,从自然因素和人为因素2个方面选取月平均气温(℃)、降水(mm)、日照时数(h)、蒸发量(mm)、相对湿度(%)和月灌溉量(m3)6个因子来分析其对土壤湿度的影响情况。

(1) 自然因素 通过对月平均气温、月降水量、月日照时数、月蒸发量、月相对湿度与月均TVDI进行Pearson统计(见表1),发现除相对湿度与TVDI呈不显著负相关外,其余因子均与TVDI呈正相关。从Pearson相关系数|R|来看,各影响因子与TVDI关系密切程度依次为月蒸发量>月平均气温>月降水量>月日照时数>月相对湿度。结果表明土壤湿度与月平均气温、月蒸发量、月降水量、月日照时数均在0.01水平上呈负相关,相关系数依次为-0.789、-0.733、-0.621、-0.557。黑方台2017—2019年气象数据见图7。由于研究区年降水量稀少且集中在夏季,与此同时研究区温度也急剧升高,导致夏季蒸发量增高而土壤湿度呈降低趋势。

表1 月均TVDI与各因素相关系数分析

图7 研究区2017年1月—2019年10月气象数据Fig.7 The meteorological data during January 2017 to October 2019

(2) 人为因素 黑方台属于大陆性干旱气候,年降水量较少,而蒸发量巨大,为了深入分析黑方台农用灌溉量对土壤含水量变化趋势的影响情况,获取2017—2019年该区域月灌溉量(见图8),通过计算每月的累计量来研究月均TVDI与灌溉量的Pearson关系,其相关系数为0.409,呈不显著正相关。研究表明在年际内随着月份的增大,农业需水量增多,其灌溉量也增多,灌溉水不仅满足地表水,随着浅层土壤水饱和而入渗成为地下水,造成该区域成为滑坡的高发区;与此同时大气温度逐渐升高,土壤温度也在升高,研究区蒸发量快速上升,导致在温度较高的季节随灌溉量的增多,TVDI值增大,而土壤中的含水率减小。

图8 研究区2017年3月—2019年11月灌溉量Fig.8 The irrigation volume during March 2017 to November 2019

3 讨论与结论

3.1 讨论

研究黑方台土壤湿度空间分布和时间变化发现,该区域年蒸发总量达到1 000 mm以上,成为继农业灌溉量后土壤湿度的第二大影响因素。由于研究区长期的农业灌溉引起地下水位上升,造成土壤抗剪强度下降,加之西高东低的地势导致研究区东南方向成为滑坡的高发区。许多学者对研究区灌溉引起的滑坡机理、过程等进行了大量的研究,张茂省等[29-30]研究了黑方台地表水的入渗、灌溉与地下水的演化关系;谷天峰等[31]研究了黑方台地下水位上升引起的黄土斜坡内部渗流场的变化特征,探讨了地下水位变化对斜坡稳定性的影响;亓星等[1]通过研究黑方台地下水位与台塬灌溉的响应规律,发现地下水具有133 d的滞后期;严冬冬[32]研究了灌溉引起的黑方台区域地下水的演化过程以及水位上升速度;朱立峰[33]对黑方台滑坡群密集发育的内在控制条件和外动力诱发因素进行了研究。

综上,这些研究都忽略了研究区气候因素——高蒸发量与滑坡的关系。上述研究表明,农业灌溉是黑方台土壤湿度上升的主要影响因素,而高蒸发量是土壤湿度降低的重要因素,说明用于农业灌溉的一部分水被蒸发损失。因此在考虑灌溉量对滑坡影响的研究中,蒸发量应作为不可忽视的一部分水分被考虑在黑方台滑坡影响因素内。在研究其他区域滑坡影响以及机理的过程中,研究区诸如蒸发量等气候因子是不容忽视的重要因素。

3.2 结论

研究利用多源遥感数据,构建指示生态干旱强弱的温度植被干旱指数,从而研究整个流域的土壤湿度变化特征,探讨其对气候变化的响应。主要结论如下:

(1) 从2017年3月—2019年11月TVDI月均值的空间分布来看,黑方台土壤湿度表现出台塬内高,边缘低的明显空间差异,这主要受研究区人为干扰的影响,台塬内多为可利用的水浇地,而台塬边缘由于经常发生滑坡,多为废弃的未利用地。TVDI月均值为0.57,说明对降水稀少的研究区而言,灌溉对土壤湿度影响很大。

(2) 从TVDI月均值的时间变化特征来看,年际土壤湿度存在减—增—减的变化态势,主要受蒸发量及灌溉量的影响,总体时间尺度上呈现土壤湿度增加的态势。

(3) 从黑方台土壤湿度的空间变化特征来看,台塬内部有小部分土壤湿度增大的区域,而大部分区域保持土壤湿度不变,此外台塬边缘土壤湿度也基本保持不变,主要受到人类活动的影响。土壤湿度减小的区域多分布在黑方台西北方向。

(4) 从土壤湿度的多个影响因子来看,月平均气温与其负相关关系最为显著,也进一步验证了研究区土壤湿度的大小主要受蒸发量和灌溉量影响,而空间分布格局受灌溉影响。