基于TVDI的不同土地类型土壤湿度趋势研究

王汉文，原喜忠，雷胜友，郭长志，褚志成，陈虹宇，陈鹏辉

(长安大学 公路学院，陕西 西安 710064)

0 引 言

在地球系统中，土壤湿度状态在陆气能量交换中起到了重要的控制作用。传统依靠地面站点统计所计算得到的干旱湿度指数虽然可以较准确地表达附近地区的干湿情况，但是获取数据成本较高，并且往往数据的数量较少，分布存在较大的差异性，在时间序列方面，很难表达土壤湿度的动态变化。随着遥感卫星技术不断发展，并伴随诸多相关研究都表明其反演得到的湿度情况可以较为准确地反映实际情况，遥感以其宏观、动态、高效、受地面条件限制少等优点为开展区域尺度下大面积土壤湿度监测提供了有效手段[1-2]。在之前的研究中，大多数学者是通过建立植被指数(NDVI)与地表温度(Ts)的综合模型进行土壤湿度监测，比如温度植被指数(TVI)、植被供水指数(VSWI)、水分亏缺指数(WDI)、条件植被温度指数(VTCI)和温度植被干旱指数(TVDI)等[3]，其中VTCI与TVDI运用较为广泛，它们均与表层土壤湿度有较好相关性[4-5]，并且反演速度较快，因此本文采用其中之一的TVDI进行相关研究。最新发射的Landsat 8卫星比以前的Landsat系列卫星在扫描方式、波段设置、辐射分辨性能等方面都有较大改进[6]，并且由于其所测数据多为开源资料，所以它所提供的影像已成为现阶段多数土壤湿度研究的数据源。运用遥感影像进行湿度反演以及旱情动态监测的研究，之前的学者已经对诸多地区进行了探索，效果较为理想，但是针对具体地区不同土地类型一年内的湿度变化趋势，却鲜有人分析。为此，本文使用Landsat 8 OLI/TIRS产品数据，采用温度植被干旱指数(TVDI)对研究区一年内的各个月的土壤湿度进行反演，结合不同土地类型分析林地、草地等土地类型的年湿度变化趋势。

1 资料与方法

1.1 研究区概况

论文研究的矩形区域位于广东省北部的韶关市内113°36′36″～114°04′12″E，24°30′10″～24°58′12″N，含武江区与曲江区两大部分，研究区域总面积约2 927.8 km2，其影像位置见图1。

图1 研究区Landsat 8遥感影像位置

该区域属于中亚热带湿润型季风气候，春季阴雨连绵，秋季降水偏少，冬季寒冷，夏季偏热。年均降雨1 400～2 400 mm，年平均气温为18.8～21.6 ℃。研究区地形以丘陵为主，河谷盆地分布其中，海拔变化较大，为65～1 200 m。土地覆被主要有林地、草地以及耕地和市区的建筑等。韶关是全国重点林区，而且还拥有十分丰富的耕地资源，对其进行土壤湿度监测反演以及不同土地类型的湿度变化趋势研究，对了解区域内部不同植被的分布发育情况和保护森林、耕地资源，发展绿色生态有着极其重要的意义。

1.2 数据与方法

1.2.1 数据源与预处理

研究所使用的Landsat 8影像数据来源于地理空间数据平台，其详细信息可以参见表1。如果2017年内某月没有可使用的影像数据，也应选取接近年份对应月份的影像。将影像加载到软件ENVI 5.3中进行裁剪，几何校正、辐射定标以及FLAASH大气校正等预处理，相关精度均符合研究标准。将处理后的影像进行监督分类，分类后处理以及精度检验之后得到分辨率为30 m的土地类型分类图(图2)，5种土地类型为草地、耕地、林地、水体、建筑用地及裸地，再对这5种用地类型进行掩膜提取处理，并对TVDI适用的用地类型进行分区统计分析，计算出相应的平均值来表示这一地类的TVDI值。

表1 遥感数据信息

图2 土地类型分类图

1.2.2 归一化植被指数

归一化植被指数(NDVI)是反映土地覆盖植被状况的一种遥感指标。该指数在一定程度上能消除地形和群落结构的阴影、辐射干扰及太阳高度角和大气所带来的噪声[6]。其被定义为近红外波段的反射值与红光波段的反射值之差除以两者之和，公式为

(1)

式中：ρNIR为近红外波段反射率；ρRED为红光波段反射率。

NDVI的取值为-1～1，负值表示对可见光高反射的物质，接近0表示地物为岩石、裸地等[7]，正值表示存在植被覆盖，且值越大，覆盖度越大。

1.2.3 地表温度

地表温度是反映地球表面能量状态的一个重要指标，在研究区干湿情况反演中，其作用是要与归一化植被指数一同构建Ts-NDVI特征空间[8]。本研究中使用辐射传导方程法(大气校正法)来计算地表温度，其原理是利用卫星拍摄影像的热辐射强度减去大气对地表的热辐射强度，最终得到地表的热辐射强度，由于不同物体表面的比辐射率不同，通过确定不同地表的比辐射率，之后根据普朗克公式即可反算出地表的温度[9]。辐射传输方程的表达式为

Lλ=[ε·B(Ts)+(1-ε)L↓]·τ+L↑， (2)

式中：Lλ为Landsat 8的热红外辐射亮度值；ε为地表比辐射率；B(Ts)为黑体热辐射强度，是地表温度Ts的函数；τ为大气在热红外波段的透过率；L↑,L↓分别为大气上行辐射强度与大气下行热辐射强度。τ,L↑,L↓这3个大气剖面参数可通过NASA提供的网站查询得到。公式变换得到B(Ts)表达式，即：

B(Ts)=[Lλ-L↑-τ·(1-ε)L↓]/(τε)， (3)

将其代入普朗克公式，得到地表温度值

(4)

式中：K1，K2为热红外波段的定标常数，K1=774.89 W/(m2·μm·sr)，K2=1 321.08 K。

1.2.4 温度植被干旱指数

温度植被干旱指数是利用Ts-NDVI特征空间提取的水分胁迫指标，它是估算陆面表层土壤湿度的一种重要方法[3，10]。对于某一区域而言，随着植被覆盖变化，植物蒸腾所转化为的潜热与显热能量就处于动态平衡之中，进而表现为地面温度的动态变化，又由于土壤的湿度状态与地面温度存在关系，于是就能间接反演出土壤的湿度状态。在之前的研究中，S.M.Moran等[11]发现Ts-NDVI特征空间呈梯形，而I Sandholt等[10]则在此基础上将其简化为三角形，并提出了温度植被干旱指数TVDI(Temperature Vegetation Dryness Index)[12]的概念，计算式为

(5)

式中：TS为影像中任意像元对应的地表温度；TSmax，TSmin分别为影像中同一NDVI值对应的最高与最低地表温度。

在简化的特征空间中，对干边(TSmax)、湿边(TSmin)进行线性拟合，即

TSmin=a1+b1·NDVI,

(6)

TSmax=a2+b2·NDVI,

(7)

式中：a1，b1为TVDI湿边方程的系数；a2，b2为TVDI干边方程的系数。将式(6)～(7)代入式(5)得

(8)

TVDI取值为[0，1]，其中湿边的TVDI值为0，干边的TVDI值为1，取值越接近0表示土壤湿度接近正常湿润土壤，越接近1，表示土壤越干旱缺水[13-15]，故土壤湿度与TVDI取值呈负相关关系。以TVDI作为土壤湿度分级标准，将土壤湿度分为5个层次[16]，如表2所示。

表2 TVDI土壤湿度分级标准

2 结果与分析

2.1 NDVI结果分析

通过逐月计算得到研究区每月的NDVI空间分布图，共12幅。均采用密度分级得到4个等级[14]，以10月份的NDVI空间分布图为例(图3)。经与土地类型分类图比对得到：NDVI值小于0的区域多为水体区域，经像素统计后得到水域面积约占总面积的1.34%；NDVI值在0～0.3之间多为建筑用地、裸地以及稀疏草地等低植被覆盖的区域，面积约占总面积的4.11%，其中0～0.2基本为建筑用地及裸地[5]；NDVI值在0.3～0.6之间主要是长势较差的耕地、草地等较高植被覆盖的区域，面积约占总面积的12.90%；对于NDVI值大于0.6的高植被覆盖区域，主要包括林地以及长势较好的耕地、草地，其约占总面积的81.65%，其中大于0.9的部分只占总区域面积的1.56%。NDVI空间分布效果基本可以较好地反映土地类型的基本情况。

2.2 Ts-NDVI特征空间及干湿边方程

2.2.1 特征空间中NDVI范围的确定

NDVI<0对应的地表类型为水体(可以认定湿度为100%)，在特征空间构建时不予讨论[4]。根据相关研究[5，17-18]发现，当植物覆盖度较小时，NDVI值很难指示该地区的植物生物量，当植物覆盖度太大时，由于植物覆盖密集，NDVI值增加延缓而呈现饱和状态，对植物覆被检测的灵敏度将降低，无法准确显示植被覆盖情况，故根据NDVI像元分布直方图及多次拟合尝试，对于0≤NDVI<0.2及NDVI>0.9的值予以删除。综上所述，构建特征空间的NDVI取值为0.2～0.9。该区间的选择也与前人的研究[5，18]接近。

2.2.2 特征空间构建及干湿边方程拟合

将反演得到的TS以及NDVI数据导入ArcGIS软件中，并提取出选定区间内某一NDVI及其所对应的最大和最小地表温度TSmax，TSmin。以一定步长的NDVI为横坐标(自变量x)，反演得到的地表温度为纵坐标(应变量y)，构建了一年内各个月份的Ts-NDVI特征空间，并拟合出了各个时相的干湿边方程，如图4所示。

图3 10月份NDVI等级分布

从图4可以看出，不同时相的Ts-NDVI散点图基本都具有相似的形状，呈三角形或梯形，与前人研究[19-21]结论相同。总体而言，地表温度最大值TSmax随NDVI增大呈减小的趋势，地表温度最小值TSmin，除几个月因删除了NDVI大于0.9后的散点图导致有减小的趋势外，其他多数均有增大的趋势，并且地表温度的最大与最小值的差值逐渐缩小。

根据Ts-NDVI特征空间，拟合得到各个时相的干湿边方程，见表3，其中干边的斜率均小于0，湿边斜率大多大于0，这也间接反映了地表温度随NDVI值变化的趋势情况。表3中干边方程的比例系数的绝对值均大于湿边方程的比例系数的绝对值，说明对该地区而言，干边所拟合的地表温度对植被指数的变化更加敏感。从拟合效果而言，干边的拟合效果较湿边整体较好，R2的平均值达到了0.837。

2.3 TVDI等级图分布

根据表3中的干湿边方程，计算得到不同时相的TVDI值，以TVDI值作为土壤湿度分级的标准，结合表2将其分为5个层次，由于在构建特征空间时未考虑NDVI<0.2及NDVI>0.9的部分，为减少土壤湿度分布误差，故再叠加水体、建筑用地及裸地并将NDVI>0.9的区域归为无数据区，进而得到本研究区年内各个月的土壤湿度分布(图5)。

2.3.1 土壤湿度空间分布特征

对于研究区域的土壤湿度状况，整体呈现出西侧较为干旱，东侧较为湿润。结合土地类型分类图可以看出，东侧大多为林地或者草地，地表的黏性土层也相对西侧较厚，故储存水分的能力相对较高，相应地，植被覆盖指数图也能反映这一方面；另一方面，东侧的地势较高且多为山地丘陵，反演得到的地表温度较低，植物的蒸发蒸腾能力较低，进而地表的土壤湿度相应较高。局部区域可以看出，建筑用地及裸地周围的土壤湿度明显较低，林地密集的区域土壤湿度则相对较高，草地和耕地之间的土壤湿度比较与时间有一定的关系，主要是存在人为因素导致耕地土壤湿度在某些特定时间段内发生较大变化。而对于山地丘陵区域，通过提取阴阳坡位置处的土壤湿度发现，相同高程阴坡的TVDI值要比阳坡低0.2左右，即湿度大约相差一个层次。

2.3.2 土壤湿度时间分布特征

土壤湿度变化在一年内大体呈现出先降低再升高的循环，湿度最低即干旱程度最严重(极干旱)主要集中在4—6月，通过ArcGIS相关的分类统计工具，计算得到了4—6月极干旱区的面积比例分别为43.67%，30.62%，14.65%，其他月份极干旱区则均小于10%。研究区西侧大片区域表现为极干旱的类型，结合土地类型、作物生长、气候变化等因素可知，西侧多为耕地及草地。农作物一般4月播种，生长期多为4，5两月，而草地的快速生长期一般也在4月，对水分需求都较大，加之该月份内区域降水较少并且地表温度也处于上升趋势，所以导致了大面积极干旱的情况。为了直观了解一年内研究区的土壤湿度整体情况，按照计算的TVDI值，将表2中的极湿润(0

图4 年内各月的TS-NDVI特征空间

表3 年内各月Ts-NDVI特征空间干湿边拟合方程

图5 年内各月份的TVDI土壤湿度等级分布

图6 一年内湿度变化整体趋势Fig.6 Overall trend of humidity change of the year

2.4 不同用地类型的土壤湿度趋势变化

结合之前的土地类型分类图2以及一年内逐月的TVDI土壤湿度等级分布图5掩膜提取出各个用地类型，并对林地、耕地、草地每月的湿度情况进行分区统计分析，计算得到TVDI的月平均值，进而绘制出这3种不同用地类型的土壤湿度趋势如图7所示。

从图7可以看出，几种用地类型的土壤湿度情况在一年的变化整体表现为：3月，TVDI值迅速增大，土壤湿度减小；4—8月湿度基本维持不变；9月，TVDI值又减小，相应土壤湿度则呈现增大趋势；之后的10月至来年2月各个用地类型的TVDI值又基本保持恒定。这一规律也基本符合该地区植被在生长周期内与土壤水分的相互作用关系。局部区域，各个用地类型土壤湿度的情况又可以分为两大阶段：第一阶段是3—10月，土壤湿度情况表现为林地>草地>耕地；第二阶段是一年内的其余时间，表现为林地>耕地>草地。两个阶段中湿度变化主要发生在耕地与草地之间。对于耕地草地之间的变化，经查阅资料得知，与该地区农作物播种后因生长需要吸收较多土壤水分以及在该时间段内降雨量大大减小有密切关系。对于林地，由于其植被覆盖相对于其他类型一直较高，地表土壤多为黏性土且性质稳定，保水能力较强，故TVDI值相对较低，湿度一直比其他用地类型高。通过计算一年内林地、草地、耕地TVDI值的变化幅度值(TVDImax-TVDImin)分别为0.226，0.241和0.324，得到3种用地类型的土壤湿度变化幅度情况为耕地>草地>林地。其中林地土壤湿度变化幅度最小，土壤水分的保持能力最强，耕地由于受到较大的人为因素影响，加之保水能力较差，土壤湿度变化幅度最大。

图7 不同用地类型的土壤湿度趋势

综上，几种用地类型在一年内土壤湿度情况为林地>耕地或草地；土壤湿度变化幅度情况为耕地>草地>林地。几种用地类型在同月中TVDI值差别不是很大，变化趋势基本相同，且均符合大区域范围下的土壤湿度变化情况，因此具有一定的参考性。

3 讨 论

本研究讨论了一年中每月的湿度情况，相对于其他研究而言，研究时期较长，更具有代表性和动态特性，但由于研究区内Landsat 8一个月中可用的影像数据十分有限，用某期的影像进行土壤湿度反演，进而用于代表该月的土壤湿度情况，其中可能存在一些偶然性因素导致较大的土壤湿度误差。

由于NDVI在低值区易受到土壤背景因素的影响，在高值区对高覆盖植被反应不敏感，故会导致TVDI的计算存在一定误差[18]。虽然本文中采取排除的方法提高准确性，但处理过程过于麻烦，故可以考虑用EVI[22-23]及SAVI[24-25]取代NDVI来构建特征空间。

在使用TVDI进行不同土地类型土壤湿度研究时受到了诸多条件限制，今后可以考虑使用其他相关性较好的湿度指数，例如VTCI等，用以验证TVDI的准确性并对比研究湿度指数对不同地类的适用性。

研究中对于用地类型的分类数量不多，且精度不高。例如未对林地中的二级分类进行详细区分，今后可以考虑选取地类丰富的区域，进行更加详细的分类与深入研究。

4 结 论

本文利用Landsat 8相关数据对广东省北部韶关市部分区域进行了土地类型分类及土壤湿度监测反演研究，得到了一年内各个月具有代表性的归一化植被指数NDVI、地表温度TS和温度植被干旱指数TVDI等数值，并结合土地类型分类，得到了以下结论。

(1)对于研究区的土壤湿度空间分布，整体呈现西侧干旱，东侧湿润的情况，局部区域与人为因素、气候条件以及地形地貌存在一定联系。

(2)对于研究区土壤湿度时间分布，在一年内大体呈现先降低至基本恒定，再升高并保持稳定的循环。湿度最低主要集中在4—6月。结合TVDI将湿度情况分为潮湿、干燥两大类，4—9月，干燥区面积大于潮湿区面积，一年中其他月份，潮湿区面积大于干燥区面积。

(3)基于TVDI，几种用地类型在一年内土壤湿度情况，3—10月，林地>草地>耕地，其余时间为林地>耕地>草地，其中耕地与草地的湿度大小情况与研究区内人为因素及气候关系密切。通过计算TVDImax-TVDImin，得到一年内几种用地类型的土壤湿度变化幅度情况为耕地>草地>林地。