基于4种植被指数TVDI模型的三江平原土壤湿度反演

陈明星, 张玉虎

(首都师范大学 资源环境与旅游学院, 北京 100048)

土壤湿度是水文学、气象学以及农业科学研究领域的一个重要参数。大型农灌区土壤湿度状况的监测是农业用水管理以及农作物旱情预报的一个重要内容[1]，一直以来受到学者们的广泛关注[2-4]。遥感手段相对于传统的土壤湿度监测手段，能够快速获取大面积区域的土壤湿度状况，对于农业旱情监测具有较大的实用价值[5]。国内外学者对土壤湿度的遥感监测方法做了许多研究，提出了不同的模型和方法，总体上分为光学遥感和微波遥感两种方式，主要有：表观热惯量法[6-7]、温度植被干旱指数[8-9]、植被供水指数[10-11]、微波极化差异指数[12]、积分方程模型[13]、Dobson模型[14]等。微波遥感方式仅适合小尺度土壤湿度反演[15]，光学遥感手段中部分方法在实际应用中也具有一定的局限性，如表观热惯量方法仅适合于裸土地区土壤湿度的估算[16]。在众多的模型与方法中，通过LST-NDVI 特征空间计算TVDI反演土壤湿度的方法，由于较强的适用性在国内外得到了广泛的研究与关注。齐述华等[17]、曹雷等[18]研究表明TVDI可有效反演区域土壤水分，且精度较高，J Chen[19]和Li Z等[20]用TVDI方法分别研究了我国黄淮海平原和陕北半干旱黄土高原的土壤湿度状况，均得到了较好的效果；N.T.Sona等[21]探讨了利用TVDI反演的土壤湿度结果对农业干旱监测的适用性，验证结果显示该方法效果较好。NR Patel等[22]评估了TVDI监测土壤湿度状况的潜力，发现TVDI在植被覆盖稀疏时能够较好的监测土壤湿度的时间变化。当前大多数研究在TVDI模型计算中仅采用归一化植被指数，但是归一化植被指数对于不同时期植被覆盖差异的敏感性不同，影响了TVDI监测土壤湿度的准确性。

本文考虑对植被覆盖度敏感性不同的4种植被指数，以三江平原为研究区，计算不同植被指数下的TVDI，并用土壤相对湿度和降水量数据对TVDI监测土壤湿度状况的效果进行验证。对比不同植被指数下TVDI模型对土壤湿度状况的监测效果，确定反演土壤湿度效果最佳的植被指数，为三江平原地区土壤湿度的监测提供参考依据。

1 资料与方法

1.1 研究区概况

三江平原位于中国黑龙江省东北部，介于43°50′02″—48°24′41″N，129°11′49″—134°46′37″E，地理位置如图1所示。是由松花江、黑龙江和乌苏里江冲积形成的低平原，土地面积为1.09×105km2，地势总体特征是西南高东北低。本区属于温带季风气候，夏季炎热多雨，冬季寒冷干燥。1月份平均气温低于-18℃，7月平均气温21～22℃，年降水量500～650 mm，且集中在夏季，属于湿润、半湿润气候区。行政区域包括佳木斯、鸡西、鹤岗、双鸭山、七台河等地级市以及牡丹江市所属的穆棱县和哈尔滨所属的依兰县，共计23个县(市、区)。三江平原是我国重要的粮食产区，主要农作物为水稻、玉米、大豆等。

图1 研究区位置

1.2 原理和方法

温度植被干旱指数(Temperature Vegetation Dryness Index，TVDI)最早是由Sandholt等[23]提出，TVDI主要考虑了归一化植被指数(NDVI)和地表温度(LST)这两个描述土壤表层特征的重要参数[24]，以此构建NDVI--LST的特征空间，Sandholt等认为NDVI-LST的特征空间中存在一系列反映土壤湿度的等值线，据此提出了温度植被干旱指数(TVDI)的概念(图2)。其公式[25-26]为：

(1)

式中：LSTmin为NDVI相同值对应的最低地表温度，为NDVI-LST特征空间的湿边，LSTmax为NDVI相同值对应的最高地表温度，为NDVI-LST特征空间的干边。根据影像像元构造的特征空间，通过线性拟合获得特征空间中的干边和湿边方程，其公式[25]为：

LSTmax=α1+b1×NDVI

(2)

LSTmin=α2+b2×NDVI

(3)

式中：α1，b1，α2，b2分别是特征空间中拟合的干边和湿边方程的系数。TVDI值在0～1之间，在干边上TVDI为1，在湿边上TVDI为0。TVDI越趋向于0，表示土壤湿度越高；TVDI越趋向于1，表示土壤湿度越低。

图2 LST-NDVI特征空间

1.3 数据来源及处理

本研究采用的遥感数据来自于美国USGS网站提供的MODIS Terra产品数据(https:∥lpdaac.usgs.gov/dataset_discovery/modis)，包括每8 d合成的空间分辨率为1 000 m的地表温度产品(MOD11A2)和空间分辨率为500 m的地表反射率数据(MOD09A1)。数据时相覆盖了2013年5—9月不同作物生长期的植被覆盖情况(东北地区植被覆盖度在一年内呈先上升后下降趋势，最高的时间段为7—8月份[27])。本次研究使用的气象数据和土壤湿度数据来源于中国气象数据网(http:∥data.cma.cn)，时间为2013年5月至9月。研究区共有12个土壤墒情站和22个气象站(图1)，选取12个与土壤墒情站对应的气象站。土壤墒情站提供了每旬的20 cm深度的土壤相对湿度数据，气象站提供了逐日降水量数据。遥感影像的选取尽可能保证与土壤相对数据时间上的对应，并覆盖5—9月的每个月份。

将MOD09A1，MOD11A2影像进行预处理后得到500 m分辨率的地表反射率和地表温度数据，反射率数据包括：红光波段(Red)和近红外(NIR)波段和蓝光波段(Blue)波段。使用ArcGIS软件分别计算得到归一化植被指数(NDVI)、增强型植被指数(EVI)、改进的修正土壤调整植被指数(MSAVI)以及比值植被指数(RVI)，在ENVI中采用基于IDL语言开发的TVDI计算工具，输入植被指数和地表温度数据获取构建特征空间的植被指数和地表温度值及计算出的TVDI数字影像。不同植被指数[25,28]的计算公式见表1。在预处理过程中，由于8月份的遥感影像云量较大，经去云处理后影像缺失范围较大，故不做分析。

表1 不同植被指数计算公式

注：ρNIR为近红光波段反射率，ρRED为红光波段反射率，ρBLUE为蓝光波段反射率。

2 结果与分析

2.1 不同植被指数构造的特征空间

TVDI在计算过程中一般不考虑水体存在的区域，利用植被指数值将这部分区域剔除。为了进一步降低误差，根据像元统计及目视解译，结合不同时期的植被指数值，选取不同的阈值(5月为0.04、6月为0.05、7月为0.1、9月为0.1)对临湖、临河或河中沙洲等距离水体较近，但由于植被指数计算误差存在少量水体或土壤水分饱和的裸地所在区域的少量像元进行剔除。在Excel软件中，以植被指数为横坐标，LST值为纵坐标，构建不同植被指数的LST-Ⅵ的特征空间，结果如图3。

根据图3，特征空间形状总体相似，地表温度的最大值均表现为随着植被指数值的增大，先有轻微上升之后逐渐下降，但总体上随植被指数增加呈减小趋势，上升区间主要处于植被指数小于0.2区间，这部分主要为裸地或有稀疏植被分布的地区，由于植被覆盖度较低，很难对区域内的地表温度产生调控作用。地表温度的最小值总体上呈现随植被指数的增加而增大的趋势，但是上升速率较小，地表温度最小值所在区域可能代表了湿地、径流量不大的水系等[26]，这些区域的地表温度在自然条件下处于比较稳定的水平。

图3 2013年5-9月LST-Ⅵ特征空间

同一时期，不同植被指数构造的特征空间存在一定区别，LST-NDVI构造的特征空间在不同时期相较于其他植被指数较差，除5月外，特征空间的形状与理论的三角形均存在差异。LST-RVI构造的特征空间在4种植被指数中效果最好，LST-EVI和LST-MSAVI构造的特征空间区别不大，效果均较好。特征空间总体上符合理论的三角形关系，与前人的相关研究结果一致[17-18]，构造的特征空间可以用来获取TVDI。

2.2 特征空间对应干、湿边的确定

在理想状态下，当地表覆盖情况满足从裸土到植被完全覆盖均匀变化，土壤表层含水量从凋萎系数到田间持水量均匀变化时，特征空间的干、湿边能拟合成直线状[29]。根据特征空间的构造结果，拟合不同植被指数的特征空间对应的干边和湿边方程，情况见表2。

从表2可以看出，特征空间的干边方程的拟合效果总体较好，R2大部分都在0.8以上；湿边方程的拟合效果总体较差，R2大部分低于0.5，三江平原存在大面积的湿地及森林，可能是由于此类覆盖物对地表温度的调控作用所导致的。植被指数对于干、湿边方程的拟合效果存在一定影响，总体上不同时期LST-RVI的干、湿边方程拟合效果最好，其次为LST-EVI和LST-MSAVI，LST-NDVI的拟合效果相对其他3种植被指数最差，干、湿边方程的R2均低于其他植被指数。总体而言，干、湿边方程的拟合结果能够满足TVDI计算的需要。

表2 不同植被指数的干边和湿边方程

2.3 TVDI的精度验证

2.3.1 TVDI与土壤相对湿度的回归分析 根据三江平原地区12个土壤墒情站点测量的土壤相对湿度数据，以及由站点坐标原位提取的TVDI值，在SPSS中利用最小二乘法进行回归分析。考虑到设置的观测站点能够代表相应地区的气候特点，所以利用站点观测数据进行验证是可行的[17]。回归分析的结果见图4。结果表明，4个时期不同植被指数计算的TVDI与土壤相对湿度均表现随着TVDI值的上升，土壤湿度均呈下降趋势，即二者之间为负相关关系，TVDI值越高，土壤湿度越低。

根据回归分析结果，在植被覆盖度较低的5月份，4种植被指数计算的TVDI与土壤相对湿度回归分析的决定系数R2均在0.52左右，总体上相差不大，且均达到0.01的显著水平。6月份植被覆盖度有所增加，该时期R2均在0.2左右，植被指数的差异对于回归分析影响较小。7月份植被覆盖度较高，该时期计算的TVDI_EVI和TVDI_MSAVI与土壤相对湿度的回归分析结果相较于TVDI_NDVI和TVDI_RVI有显著优势，R2明显高于TVDI_NDVI和TVDI_RVI，且达到0.05的显著水平，TVDI_NDVI和TVDI_RVI的回归分析结果未达显著水平。相关研究也表明EVI和MSAVI 对高植被覆盖区域更为敏感，由此计算的TVDI比其他植被指数更能反映高植被覆盖时的土壤湿度状况[26,30]。9月植被覆盖度处于较低水平，决定系数R2均在0.55左右，植被指数的差异对于回归分析结果影响较小，该时间段TVDI反演的土壤湿度和实测的土壤相对湿度相关性均较高，且达到0.01的显著水平。

对比4种植被指数计算的TVDI反演的土壤湿度和实测的土壤相对湿度的回归分析结果。在低植被覆盖的5月和9月，决定系数R2均显著高于植被覆盖度较高的6月和7月的分析结果，低植被覆盖时的反演效果好于高植被覆盖时。不同植被指数在5月、6月和9月份时R2相差较小，反演效果差异较小；在7月份，增强型植被指数(EVI)和修正土壤调节植被指数(MSAVI)的R2高于另外两种植被指数，反演效果有所提升。

2.3.2 降水量与TVDI的相关分析 土壤湿度在一定程度上受降水量影响[31]，研究中采用降水量探讨其与TVDI和土壤相对湿度的关系，进一步确定TVDI方法监测土壤湿度状况的准确性。主要是将降水量与土壤相对湿度和TVDI做相关分析，结果见表3。

根据相关分析结果，7月份的降水量与土壤相对湿度呈显著正相关，即降水量增加，土壤相对湿度会增大，其他时间段均为不显著正相关。综合表明降水量对土壤相对湿度存在一定影响。

不同植被指数计算的TVDI和降水量的相关分析结果中，仅5月份的NDVI和EVI计算的TVDI通过了显著性检验，其余均未通过显著性检验。在5月、6月和9月期间，植被指数的不同对相关分析的结果并未造成明显的差异，同一时期不同植被指数的相关系数值均较为接近；在7月份，基于EVI和MSAVI计算的TVDI值与降水量的相关系数明显高于另外两种植被指数。总体情况与前文TVDI和土壤相对湿度的回归分析结果一致，均表现为除7月外，其他月份不同植被指数计算TVDI验证效果差异较小。降水量与不同植被指数计算的TVDI值的相关系数均为负值，即降雨量越大，TVDI值越小；结合降雨量与土壤湿度的相关分析结果(降雨量越大，土壤湿度越高)，表明TVDI值越小，土壤湿度越大，与前面分析结果一致。

表3 2013年5-9月降水量分别与土壤相对湿度及TVDI相关系数

注：*表示在0.05水平上显著相关。

图4 2013年5-9月土壤相对湿度与不同植被指数TVDI的拟合效果比较

2.4 三江平原土壤湿度状况时空分布

根据以上分析结果，MSAVI和EVI计算的TVDI在分别在5月和7月回归分析R2最高，且特征空间拟合效果较好；NDVI计算的TVDI在6月和9月回归分析R2最高，但特征空间拟合效果较差，EVI和RVI的R2分别仅次于NDVI且特征空间拟合效果较好。综合考虑，在不同的时间选取相应的植被指数计算TVDI(5月为MSAVI，6月为EVI，7月为EVI，9月为RVI)并绘制土壤湿度等级图，分析不同时间三江平原的土壤湿度状况。根据TVDI值将研究区的土壤湿度状况划分为7个等级[32]，分别为：极度湿润(0～0.1)、中度湿润(0.1～0.2)、轻度湿润(0.2～0.4)、正常(0.4～0.6)、轻度干旱(0.6～0.8)、中度干旱(0.8～0.9)、极度干旱(0.9～1)。三江平原土壤湿度分布情况见图5。

图5 不同时期土壤湿度分布

由TVDI反映的三江平原地区的土壤湿度分布状况可知，5月全区受旱范围比较广泛，主要分布在三江平原的西部，受旱地区多为轻度干旱，中部的集贤县、西南部鸡东县以及密山市西部地区出现了中度干旱，同时南部地区存在小范围的极度干旱。东部地区多为轻度湿润或正常情况，土壤湿度状况分区明显，呈现显著的从东北到西南由湿润向干旱变化的特征。6月全区大部分区域土壤湿度状况正常，中西部及西南部区域存在轻度干旱，中部存在零散的中度干旱，范围较小，主要集中在富锦市、宝清县和集贤县。7月全区大部分地区为正常，中部及南部地区存在零散轻度干旱，中部的富锦市部分地区存在中度干旱，在东部饶河县的西部地区发生了中度干旱，范围较小，全区总体上土壤湿度状况良好。9月份全区受旱范围广泛，除三江平原东北部以及西北部地区外，全区其他地方均发生不同程度的干旱，干旱等级主要为轻度干旱，西部的勃利县以及中部的宝清县部分区域为中度干旱，研究区的土壤湿度状况整体呈现偏干旱的情势。

对比4个时期的土壤湿度分布图，在三江平原东部及东北部，土壤湿度一直保持为正常或轻度湿润。相关研究[33]表明三江平原东部及东北部地物类型主要为水田或沼泽等湿地及森林，可能是森林具有的对地表温度的调控作用[34]，间接影响了土壤湿度。湿地区域由于土壤水分长期处于过饱和状态，遥感手段反演的土壤湿度状况显示为正常或轻度湿润。

3 讨论与结论

3.1 讨 论

针对单一植被指数对不同植被覆盖状况下敏感性不同的问题，本文对比分析了不同植被指数下TVDI对土壤湿度状况的监测效果，得到了合理的结论，TVDI能够适用于三江平原的土壤湿度状况监测，与前人研究[9,18,24,29]结论一致。在不同时期选用何时的植被指数计算TVDI能够更好的反演土壤湿度状况。但是在研究中仍然存在不足之处有待提高：(1) 本研究在不同植被覆盖状况下选取不同的植被指数计算TVDI，虽然能够在一定程度上提高土壤湿度状况反演的准确性，但是土壤湿度状况受多种因素的影响，下一步可以考虑建立综合多种影响因素的土壤湿度状况监测模型，进一步提高土壤湿度状况监测的准确性。(2) 考虑到遥感影像和土壤湿度数据的时间对应及云量问题，本研究选用了四期遥感影像进行研究，对于不同植被覆盖状况为定性研究，下一步可以考虑实测土壤湿度，选用多期遥感影像进行分析，对于植被覆盖状况进行定量化，从而提高TVDI模型在时间上的适用性。(3) 基于MODIS遥感数据能够监测大尺度土壤湿度状况，本研究中也取得了较好的结果，但是MODIS影响容易受云和雨雪天气影响，且分辨率较低，后期可以考虑利用雷达影像在小区域做更为细致的研究。

3.2 结 论

(1) TVDI方法能够用来反演三江平原地区土壤湿度状况,基于不同植被指数的TVDI与土壤湿度基本上呈负相关关系,即TVDI值越大,土壤湿度越低。(2) 不同植被指数计算的TVDI在5月、6月、9月这3个时期,回归分析R2数值相近(5月均在0.52左右、6月均在0.20左右、7月均在0.55左右),即不同植被指数的TVDI在这3个时期对土壤湿度的反演效果差异较小,4种植被指数均适合用来反演这3个时间段的土壤湿度;7月份的EVI和MSAVI计算的TVDI的R2(均在0.36左右)明显高于植被指数NDVI和RVI的R2(均在0.15左右),表明EVI和MSAVI计算的TVDI在7月份对土壤湿度的反演效果好于其他两种植被指数, EVI和MSAVI更适合用来反演该时期的土壤湿度。(3) 三江平原2013年5—9月大部分地区湿度状况属于正常或轻度干旱,干旱现象主要出现在中部和西南部;东部及东北部地区由于森林和湿地的覆盖,在不同时期均保持在正常和轻度湿润状况。