毕业论文-利用TVDI反演地表含水量方法研究

1、PAGE 第19页（共19页） 利用TVDI反演地表含水量方法研究作 者：指导教师：摘要：利用MODIS资料构建了地表温度(Ts)-增强植被指数(EVI)的特征空间,拟合了特征空间中的干、湿边方程,计算了温度植被干旱指数(TVDI),并推导出利用TVDI和干、湿边土壤水分计算土壤含水量的方程。在计算TVDI的过程中,为了减少高程的影响,利用数字高程模型(DEM)对Ts进行了订正;利用同期野外实测土壤湿度数据计算了干边上的土壤水分值,从而反演出平均土壤含水量。结果表明:TVDI方法能反演土壤表层水分,实测值与预测值之间平均绝对误差在15个百分点左右。高程校正后的TVDI能更好的反映土壤水分,与校

2、正前相比,平均绝对误差减少5个百分点,基本满足业务需要。使用野外与卫星同步采样的土壤湿度数据进行验证,发现TVDI指标与实测土壤湿度数据显著相关,能够较好地反映表层土壤湿度。关键词：TVDI；土壤水分；MODIS；地表温度;DEM土壤水分是监测土地退化的一个重要指标,是气候、水文、生态、农业等领域的主要参数,在地表与大气界面的水分和能量交换中起重要作用。遥感能够快速方便地获取大区域的地表信息,因此使用遥感监测土壤水分意义重大。20世纪80年代以来,以遥感监测地表辐射温度信息反演土壤水分的方法如应用于稀疏植被的热惯量法得到了广泛应用。在多数情况下,地面一般为植被不完全覆盖,因此遥感探测到的地表温

3、度必然受到植被覆盖度影响,较高的土壤背景温度会严重干扰旱情信息。为了消除土壤背景的影响,将地表温度和反映植被覆盖度的植被指数联合考虑将会是提高反演精度的有效途径。Price、Carlson等人发现植被指数(NDVI)与地表温度(Ts)之间存在负相关关系,而且如果研究区植被覆盖包含从裸土到全覆盖、土壤湿度从极干旱到极湿润的各种情况,以遥感资料获得NDVI和Ts为横、纵坐标的散点图呈三角形。在这个三角形特征空间中,Goetz研究发现Ts与NDVI的斜率与土壤湿度之间为一元线性关系1, Sandholt等进一步结合对该特征空间生态特征的解释,提出用温度植被旱情指数(TVDI,Temperature-

4、Vegetation Dryness Index)估算土壤表层水分状况,取得了良好的效果10。国内齐述华、王鹏新等利用TVDI进行了全国旱情监测研究,结果表明TVDI作为旱情评价指标是合理的11-12,姚春生等发现TVDI指标与实测表层土壤湿度数据显著相关,作为指示地表土壤湿度的指标有一定的合理性2。卿清涛等基于TVDI方法对四川省8月上旬的旱情进行了遥感监测评估3。国内学者主要探讨了TVDI作为旱情评价指标的合理性,并未深入利用TVDI反演土壤水分值。另外由于构成TVDI的一个重要因子是地表温度,高程因子对地表温度影响较大,对于地势起伏较大的地区应用TVDI方法时应考虑高程对TVDI的影响。

5、利用高程对TVDI进行校正方法也较少报道。本文拟采用MODIS卫星资料,利用TVDI方法反演土壤水分,并利用DEM数据对TVDI进行校正,结合野外同步实测土壤水分值利用TVDI方法反演地表土壤水分,并对高程校正前后的结果进行对比分析。1原理和方法1.1数据来源及预处理 选取江苏省仪征市作为研究区域。仪征市位于1190211922E,32143236N,总面积903 km2。地表温度和植被指数数据来源于美国航空航天局(NASA)陆地数据分发中心(LandProcesses DAAC)的MOD11A2和MOD13A2数据,对应时间为2004年5月0925日、9月20日10月6号及11月0117日的

6、合成数据。MOD11A2是8 d合成空间分辨率为1 km的陆地表面温度产品,包含白天、夜间地表温度、31和32波段通道发射率等资料。MOD13A2是16 d合成空间分辨率为1 km的植被指数产品,包含NDVI、EVI、红光、近红外、中红外、蓝光等波段反射率及其它辅助信息。为保证遥感数据时间同步性,将2个时相8 d合成的地表温度数据合成为1个时相的16 d合成地表温度数据。合成时,若2个时相都有值,取平均值作为合成后像元值;若只有1个时相有值,取其作为合成像元值。地面数据采用实测土壤湿度数据,距离间隔为5 km,网格点40个,采样点覆盖大面积区域,周围地物单一,保证采样数据与MODIS数据产品匹

7、配。采样使用GPS定位。采样时间为2004年05月1219日、09月2430日及11月0818日。采样时间持续7 d左右。采样深度为3545cm,采用烘干测重法测定土壤湿度。1.2反演模型的建立植被依据生态系统中水、热、气等状况,调控其内部与外部的物质、能量交换,植被覆盖度变化是地球内部(土壤母质、类型等)与外部作用(气温、降水等)的综合结果,是区域生态系统环境变化的重要指示。植被覆盖度作为植被的直观量化指标,很大程度上反映植被的基本情况,是研究水文、气象、生态等区域或全球性问题的基础数据。Sandholt等在研究土壤湿度时发现,Ts-ND-VI的特征空间中有很多等值线,于是提出了温度植被干旱

8、指数(Temperature-Vegetation Dryness In-dex,TVDI)的概念。TVDI由植被指数和地表温度计算得到,只依靠图像数据,其定义为:TVDI=Ts-TsminTsmax-Tsmin(1)其中:Tsmin表示最小地表温度,对应的是湿边;Ts是任意像元的地表温度;Tsmax=a+bNDVI为某一NDVI对应的最高温度,即干边;a、b是干边拟合方程的系数。在干边上TVDI=1,在湿边上TVDI=0。对于每个像元,利用NDVI确定Tsmax,根据T在Ts-NDVI梯形中的位置,计算TVDI。TVDI越大,土壤湿度越低,TVDI越小,土壤湿度越高。估计这些参数要求研究区域

9、的范围足够大,地表覆盖从裸土变化到比较稠密的植被覆盖,土壤表层含水量从萎蔫含水量变化到田间持水量。Moran等4在假设Ts-NDVI特征空间呈梯形的基础上,从理论上计算梯形4个顶点坐标的研究结果表明,在不同的植被覆盖度条件下,Ts-NDVI特征空间种最低温度(Tsmin)随植被覆盖度而不同。本文在将Ts-NDVI特征空间简化处理为三角形的同时,对Tsmax和Tsmin同时进行线性拟合,拟合的方程为:Tsmax=a1+b1*NDVI (2)Tsmin=a2+b2*NDVI (3)这样可以由下式计算温度植被干旱指数:TVDI=Ts-(a2+b2*NDVI)(a1+b1*NDVI)-(a2+b2*N

10、DVI)(4)其中:a1、b1是干边拟合方程的系数,a2、b2是湿边拟合方程的系数。2实验部分2.1温度植被干旱指数计算在NDVI和Ts构成的三角形特征空间中(图1),将不同植被指数条件下的最高下垫面温度(Tmax)相连,构成了三角形的干边,代表了该区域内的干旱上限,将其干旱指数定义为1。相应地将不同植被指数下的最低下垫面温度(Tmin)相连构成了三角形的湿边,代表了该区域内的最湿润区,将其干旱指数定义为0。基于此, Sandholt等5提出了温度植被干旱指数(TVDI)估测土壤表层水分状况。TVDI=Ts-TsminTsmax-Tsmin(1)Tsmax=EVIa+b;Tsmin=EVIc+

11、d (2)式中,Tsmax、Tsmin分别由植被指数与地表温度根据干边、湿边线性拟合获得,本研究中采用EVI植被指数(EnhancedVegetation Index), a、b、c、d分别是干边和湿边线性拟合方程的系数。在三角形区域内任一点的TVDI值介于0和1之间,TVDI值越大,对应的土壤湿度越低, TVDI值越小,对应的土壤湿度越高。图1Ts-NDVI特征空间示意图 2.2土壤水分计算NDVI-Ts的直线斜率可以反映土壤湿度状况。许多研究发现,在Ts-NDVI的特征空间中,存在一系列土壤湿度的等值线,这些等值线都可近似认为相交于干湿边的交点,并且该等值线近似为直线(图2),该直线的斜率

12、与土壤湿度间可以描述为一元线性函数关系。图2温度植被干旱指数(TVDI)原理示意图 图2中, a点的土壤湿度可以表示为:RSMa=a1+a2AC(3)式中,A、C标注在图2中。A/C为Ts-NDVI直线的斜率, a1、a2为方程系数。a点所在等值线上任一点的土壤湿度可以表示为:RSM=a1+a2ABtan(4)A=Ts-Tmin; B=Tmax-Tmin。tan=B/C,为一常数。进一步推导有:RSM=a1+a2Ts-TminTmax-Tmin=a1+a3TVDI (5)其中: a1是TVDI=0时的土壤湿度值,即Ts-NDVI特征空间中湿边上的土壤湿度, a3是TVDI=1时干边上的土壤湿度

13、减去湿边上的土壤湿度,进一步得到:RSM=MRSW-TVDI(RSMW-RSMD) (6)式中, RSM表示某一像元土壤相对含水量;RSMW是湿边对应的土壤相对含水量(最大);RSMD是干边对应的土壤相对含水量(最小)。用TVDI进行土壤湿度反演的影响因子分析用TVDI反演土壤湿度的原理,是基于水分蒸发和植被蒸腾对地表温度的降低作用。TVDI的值取决于tsNDVI空间的情况,对TVDI的影响,来源于对tsNDVI空间的影响。对tsNDVI空间影响因子的探讨有很多。Friedl等利用航空数据,研究了堪萨斯北部草原ts与NDVI的关系,发现ts与NDVI关系受到植被覆盖及土壤湿度情况的影响。Goe

14、tz研究了不同传感器获取的ts、NDVI之间的关系,发现NDVI/ts斜率在很大程度上取决于植被情况和土壤湿度,与传感器类型关系不大。Goetz的研究还发现,ts对土壤湿度的敏感性强于NDVI。Sandholt总结了对ts、NDVI关系影响因子的研究,认为ts、NDVI关系的影响因子有:植被覆盖、土壤水分蒸发蒸腾、表面热特性、净辐射、大气动力及表面粗糙度等。这些因子都影响地表温度,从而影响tsNDVI空间的情况6。Ge-oward用基于能量和水分平衡原理的SSIB模型和田间实测资料,研究各因子对ts、NDVI关系的影响情况,认为土壤湿度和入射辐射是决定ts、NDVI关系的重要因子,而风速和水汽

15、压等影响较小,可以忽略。土壤湿度是决定ts、NDVI关系的重要因素,从tsNDVI空间得到的TVDI,能从一定程度上反映土壤湿度的状况。ts、NDVI的关系不是由土壤湿度唯一确定,因而用TVDI反演土壤湿度,受到太阳入射辐射、植被覆盖、表面热特性、大气动力、水汽压等因子的影响。这些因子都影响地表的能量平衡,使平衡状态对应的ts值发生变化,从而改变tsNDVI空间的情况,影响TVDI值。下垫面与大气间湍流形式的热交换,使地表温度在一定程度上受到气温的影响:一般情况下,随高程增加,气温降低7。气温主要是通过大气与下垫面间热交换,即湍流的形式来影响地温。下垫面与大气间的热交换,与地表温度和气温的差值

16、相关。当气温变化时,地气温差发生改变,热交换能量发生变化,地表温度也随之变化。在进行热交换的过程中,地温影响到气温,反过来气温也影响地温。地温和气温之间是一种互相影响的关系。地表温度受气温的影响与高程相关。中国高程跨度近9 000 m,从西到东形成3级阶梯的地形,地表起伏明显,气温分布差异大,地表温度受气温影响显著。由于随高程增加,气温对地表温度的降温作用,TVDI与土壤湿度的关系也受到一定的影响。图3中的斜线是TVDI的等值线,代表相同的土壤湿度。越靠近干边的等值线,代表的土壤湿度越低;越靠近湿边的等值线,代表的土壤湿度越高。由于水分蒸发及植被蒸腾降低地表温度,在相同NDVI条件下,土壤湿度

17、高的像元对应ts值低,像元对应(ts,NDVI)点在空间中的分布更接近湿边,反演得到的土壤湿度高,如图中G、H点;但是由于气温对地表温度的影响,使得相同NDVI条件下高程大的像元对应的ts值有所降低,(ts,NDVI)点的分布也更接近湿边,高程大的地区反演得到的结果,也同样会表现出土壤湿度高的特性,如图中的E、F点。图3 高程和土壤湿度对ts-NDVI空间分布的影响从反演获得的全国TVDI分布图(图4)中可以更直接地看出,随高程增加,气温对地表温度的降低作用,给土壤湿度反演带来的影响。反演结果在青藏高原海拔高区土壤湿度明显偏高,甚至高于南部丘陵区。这种不合理的分布,是由于青藏高原区高程远远高于

18、南部地区,气温对地表温度的降低显著,因而虽然该区的土壤湿度小于南部,但是其区内像元对应的(ts,NDVI)值,在tsNDVI空间中的分布却更靠近湿边。图4 未经DEM订正的9月上旬全国TVDI分布图地表温度受到气温的影响,与高程相关,因而气温对土壤湿度反演过程的影响,可以用DEM来进行订正。先对ts用DEM进行订正,再与NDVI构建经过DEM订正的tsNDVI特征空间,并计算TVDI。这种方法可以减小气温对地表温度的影响,使得像元对应的(ts,NDVI)值,在tsNDVI空间中的分布情况,能更多地取决于土壤湿度。3全国土壤湿度反演的DEM订正3.1对TVDI反演全国土壤湿度进行DEM订正的方法

19、在进行土壤湿度反演时,先对地表温度进行订正,所用公式为td=ts+Ha (2)式中:td为订正后的地表温度;ts为订正前的地表温度;H为高程值;a是高程对土壤湿度反演过程的影响程度,也是地表温度受气温影响随高程增加的降低程度。用订正后的地表温度td与NDVI构成温度植被指数空间,对NDVI以0.01的步长求出不同NDVI值对应的最高、最低地表温度td max和td min,然后将td max和td min分别与NDVI进行线性拟合,获得干湿边方程td max=a1+b1NDVI, td min=a2+b2NDVI,根据方程求出任一像元NDVI值对应的最高、最低地表温度,带入式(1)中,计算出经

20、过DEM订正的温度植被干旱指数TVDI。3.2高程影响程度值的确定影响程度值表示的是地表温度受气温影响随高程增加的降低程度。一般认为,高程每升高100m,气温降低约0.6左右,而对受气温影响地表温度随高程变化情况的研究还不多。由于一般情况下,地表温度值和气温值相差不大,变化趋势也一致,认为地表温度受气温影响的变化情况,可用气温的变化作为参考。计算9月上旬全国ts与DEM 2幅图的相关系数发现,它们的相关系数R达到-0.58,呈显著负相关。随高程增加,气温对地表温度的降温作用,可以部分解释这种负相关关系,因而值应该在地表温度随高程增加降低程度的范围内。用9月上旬获取的地表温度与高程数据,分析ts

21、随高程增加的降低程度t（图5）。由于6 000 m已经在雪线高程以上,对6 000 m的情况不予考虑。对06 000m高程范围,以100 m为增量,计算各个高程区间地表温度th和高程的平均值Hh,按Hh升序排列,并计算相邻区间th和Hh的差值,得到t=(th1-th2)/(Hh2-Hh1)。为了能与气温变化情况对应,将t表示为以1/100为单位的形式。对t值与高程的关系进行分析,发现随高程的变化,t值波动较小,t值大致分布在-2/1004/100之间。高程与地表温度负相关,t出现负值可能是因为受到太阳辐射、地形、地表覆盖等的影响。由于值在ts随高程增加降低程度的范围内,的取值应该在-2/100

22、和4/100区间。图5 高程h对地表温度的影响程度通过相关性分析的方法来进一步确定值。取不同影响程度,得到全国土壤湿度的反演结果,通过探讨土壤湿度实测值与反演结果相关性的变化情况,来寻找最佳值。实测数据是从分布于全国各地的98个气象站点收集的同时期表层20 cm土壤相对含水量数据。虽然研究中所用的分辨率数据与气象站点数据在空间上难以对应,但考虑到气象观测点实际代表所在区域的整体气候特点,可以反映气象要素分布的区域性趋势,认为利用气象站点实测数据进行验证分析,是目前条件下可行的方法。对值分布的区间进行探讨,取最大值4/100进行订正,发现反演结果中高程的作用已经很明显,反演结果与实测值的相关性很

23、低,表明所取影响程度偏大,值应1/100时,反演结果与实测值相关性低,而当影响程度在01/100区间时,相关性高(图6)。在01/100区间,取步长为0.1/100进行循环搜索,寻找最佳影响程度值,发现相关系数在01/100区间呈先升后降的趋势,在a设为0.6/100时,相关系数值最大,TVDI与实测数据有最好的相关性,反演得到TVDI值能最好地反映土壤湿度的状况。这个值与气温随高程变化的值也相近,说明取影响程度为0.6/100是可信的。对取不同值反演得到的TVDI结果进行对比也发现,当取0.6/100时,反演得到土壤湿度结果最符合土壤湿度分布的一般规律,呈现出明显的南高北低,东高西低”的分布

24、特性,并且沙漠区缺水的特性也能很好地表现出来(图7)。图6 影响程度值与相关系数R的关系 图7 经过DEM订正的9月份上旬全国TVDI分布图根据相关性比较及反演结果等的分析,得到9月份高程对土壤湿度反演过程的影响程度值为0.6/100。3.3DEM订正前后的土壤湿度对比取为0.6/100,对9月上旬地表温度进行订正,获取订正后的土壤湿度,对订正前后获取的结果进行了进一步的分析。3.4反演土壤湿度结果的变化分析经过DEM订正后,青藏高原区TVDI值增加,而南方及东北地区TVDI值降低,全国土壤湿度的整个分布趋势更加合理(表1)。订正后获取的TVDI值,随实测土壤湿度的增加,呈现降低的趋势,而订正

25、前的TVDI值,与实测土壤湿度没有很好的对应关系。由影响程度的获取过程可知,订正后土壤实测值与TVDI的相关性增强。经过订正之后,TVDI与土壤湿度实测值的相关性R增加了0.1,而TVDI与实测值回归得到的线性方程R2增加了009,经过订正后的TVDI,能更好地反映土壤湿度的状况(图8)表1 订正前后TVDI值的对比 订正前 订正后图8 订正前后土壤湿度实测值与TVDI关系3.5土壤湿度的空间分布趋势分析分析进行DEM订正前后获取的TVDI分布图可以发现,经过DEM订正的TVDI分布图,能更好地表现土壤湿度空间分布的整体趋势。未经DEM订正的TVDI图,土壤湿度值最高值出现在西藏东部,并且西藏

26、、青海的土壤湿度高于南部,北部土壤湿度值也高于南部,这种分布规律不符合实际情况。经过DEM订正后的TVDI图,土壤湿度呈现明显的“南高北低,东高西低”的分布特性。新疆塔克拉玛干沙漠和古尔班通古特沙漠区、青海柴达木盆地及内蒙古巴丹吉林沙漠和腾格里沙漠干旱缺水的特性,能很好地表现出来;喜马拉雅山北麓雨影区特征明显;新疆天山北麓和塔里木盆地北部、宁夏及河套地区的绿洲特征明显;黑龙江三江平原及嫩江流域土壤湿度高,与站点实测值情况吻合;云南、四川、海南岛及广西、湖南省土壤湿度高。这些分布规律与中国气候资源地图集8中的10月份中国平均相对湿度图有很好的对应关系。9月上旬气象数据表明,这一时期降雨主要集中在

27、四川东北部、甘肃东南部、陕西南部、山西南部、河南、河北南部及江苏北部区,订正后获取的结果,这些地区土壤湿度明显偏高。反演得到的结果与降雨数据在空间分布上有很好的对应关系。对订正前后获取的TVDI分布图,用中国水文分区图9进行了叠加分析,发现订正后的TVDI分布图,能更好地体现水文分区的分布情况。4结果分析 将两个时相8 d合成的地表温度数据合成为一个时相的16 d合成的地表温度数据,并提取01范围内相同Fr下的不同像元对应的所有温度中的最大、最小陆地表面温度,Fr的精度为0.01,获取3个月的每16 d的Ts/Fr特征空间(图9)。由图9可以看出, 3个月的特征空间的干边和湿边形状相似。随着F

28、r的增大,陆地表面温度的最大值减小,最小值升高,且与Fr呈近似线性关系。图9 仪征市2004年5月、9月、11月的Ts/Fr特征空间表2 Ts/Fr特征空间的干边和湿边方程选取Fr在0.100.80线性趋势区域进行拟合,表2为拟合的干边和湿边方程,各方程均通过0.01置信度的t检验。反演土壤相对湿度,将实测数据随机分为建模和检验数据,其中建模数据样本数26个,占总样本数2/3,检验数据14个,约占1/3。使用40 cm深度建模数据的土壤湿度与Ts/Fr斜率值进行拟合,求取RSMD,即干边最小土壤湿度方程。拟合方程为:RSMD=0.667+23.847。(10)依据公式(8),计算得到仪征地区的

29、土壤相对湿度(图10)。图10 仪征市2004年5月、9月、11月土壤湿度反演分布图从图10可以看出,仪征市5月、9月和11月土壤湿度值有明显变化, 5月土壤湿度值最大, 11月土壤湿度值最小。由于受到江淮流域初夏梅雨天气影响,仪征市5月的土壤湿度值明显大于其它两月,同时5月土壤湿度值分布不均,南部及北部较高,中部区域较低,这与仪征南部属于长江冲积平原区,北部为谷底平原,河流汇集,灌溉便利有关,而中部属于缓岗丘陵区,由于地势原因,雨水不易存留,造成土壤湿度较小。11月,全市土壤湿度小且分布较均匀,这与11月仪征市降水减少有很大关系。表2为土壤相对湿度的反演与实测数据的统计分析。可以看出,实测与

30、反演平均值非常接近,尤其是9月和11月,两者差值的绝对值不到1。相对误差较小,平均值为7.64%;最小为9月,4.01%。均方根误差与相对误差类似。由此可见,本研究所采用的遥感方法反演出的土壤湿度精度较高,可服务于实际应用。表3 土壤湿度的统计量值 % 5结论与讨论本文推导了利用TVDI和干、湿边上的土壤湿度值计算土壤水分的方程，利用该方程反演了仪征市土壤表层的相对湿度。验证结果表明反演值和实测值之间的平均绝对误差为15个百分点左右。采用TVDI方法反演土壤水分是基于Ts-EVI的直线斜率与土壤水分之间为一元线性关系,这种关系被许多研究证明,但是否是适合反演土壤水分的最佳模型有待进一步分析。对

31、TVDI进行简单的高程校正,从一定程度上减少了由于高程增加气温降低,从而地表温度降低对TVDI指数的影响。但是地表温度、气温和高程之间的关系受高程、地表覆盖等许多因素影响,对TVDI指数进行更加精确的高程校正还需要进一步研究。参 考 文 献1.Goetz SJ.Multi-sensor analysis ofNDVI, surface temperature and biophysical variables at amixed grassland siteJ. International Journal ofRemote Sensing, 1997, 18(1): 71-94.2.姚春生,张增祥,汪潇.使用温度植被干旱指数法(TVDI)反演新疆土壤湿度 J.遥感技术与应用, 2004, 19(6):473-478.3.卿清涛,侯美亭,张顺谦.基于NOAA/AVHRR资料的四川省干旱监测方法初探J.中国农业气象, 2008, 29