ENVI土壤湿度遥感反演与干旱灾害监测

遥感建模与应用综合实习
实验报告
学期2017-2018学年第二学期
姓名
学号
指导教师闵爱莲
实验题目:土壤湿度遥感反演研究
1.实验目的
1.1熟悉遥感图像大气校正、几何校正的流程，掌握NDVI的计算及密度分割的过程。

1.2 学会利用TM热红外波段反演地表温度，了解TVDI的原理，制图输出黄骅市TVDI图像。

1.3 掌握基本的统计学评价指标（R2值、相关系数、线性拟合模型精度评价等），建立干旱指数TVDI和土壤相对湿度之间的关联模型。

2.实验要求
2.1 大气校正，几何校正后NDVI计算并制图输出。

2.2地表温度反演；Ts-NDVI散点图制作；TVDI指数计算；除云、除水、除建筑物，保留植被区域的TVDI。

2.3建立TVDI和土壤相对湿度之间的关联模型。

模型精度评价。

土壤相对湿度制图。

3.实验数据
TM200205图像、BASE IMAGE图像、研究区矢量边界图。

4.实验步骤:
4.1图像预处理以及NDVI制图。

①大气校正。

本实验选择黑暗像元法进行大气校正。

打开TM200205,快速统计图像DN值，记录7个波段每个波段的黑暗像元。

本实验选择DN值像元个数为十位数时对应的DN值作为黑暗像元DN值。

所有像元均减去此黑暗像元的DN,达到校正目的。

所依据原理为认为黑暗像元的辐射记录值为0，且认为图像各处大气影响一致，之所以不为0是由于大气程辐射的影响，减去此DN，即为整体消除大气程辐射的影响。

②几何校正。

加载并打开base image和TM200205图像，基于Map下的图
像到图像校正方法进行几何校正。

以base image为基准图像，以TM200205为待校正图像。

本实验选择20个地面同名地物点作为控制点，总的均方根误差为0.8257个像元，选择的原则是：控制点要分布均匀，图像边缘部分要多选取。

选择大桥桥头、河流拐弯处、道路交叉口、机场等易于区分且不易随时间变化太大的地物。

几何位置变换公式为多项式变换，二次多项式。

重采样方法为线性插法。

③裁剪图像。

基于ENVI菜单栏下的Basic 工具里的裁剪数据功能进行，空间子集选择边界图层。

裁剪之后，用边界图层.evf文件数据进行腌膜，得到研究区数据，其他地方为0值。

④NDVI计算。

基于公式（float（b1）-float（b2））/（b1+b2）或者Transform 功能下的NDVI计算工具直接进行计算。

⑤NDVI密度分割。

将NDVI从-0.6800~0.4208全部分为8类，外加一个0~0，即背景值，赋值为white
⑥制图输出。

在Annotation下设置图名、图例、指北针、比例尺。

需要指出的是，本实验在图例方面只加8个（不加0~0的图例），所有汉字均采用ENVI 字体中的161~180里的KaiTi，所有英文全部采用Roman3字体，输出如图：
4.2热红外地温反演及TVDI制作（以下所有公式均在Band Math下输入）。

①单位光谱围幅亮度值Rb。

输入公式：
(B1*(1.896-0.1534)/255+0.1534)/1.239
②地温计算，用传感器温度表征地表温度Ts。

输入公式：
1260.56/alog(60.766/B1+1)
③图像MNDWI计算。

输入公式：
（float（b1）-float（b2））/（b1+b2）
式中，B1为绿光波段，B2为短波红外波段（第五波段）。

④TVDI指数计算原理。

植被指数NDVI 是利用植被对太阳辐射各光谱波段的吸收和反射作用，得出反映植被生长状况的信息，因此可作为干旱监测的一种指标，但其监测作物干旱状态具有一定的滞后性；而地表温度Ts的实时性强，却易受到周围地物的影响，因此常将两者综合起来研究干旱。

Price 和Carlson 等研究首先发现，如果研究区的植被覆盖度和土壤湿度值变化较大，则获得的植被指数NDVI 和地表温度Ts所形成的散点图呈三角形关系，其中NDVI 作为横坐标，Ts作为纵坐标。

Moran 等则分析NDVI 和Ts的散点图为梯形关系。

Sandholt等通过研究简化的NDVI-Ts 三角形空间，提出了温度植被干旱指数TVDI，用以估算土壤含水状况。

NDVI-Ts特征空间的示意图,体现了Ts与NDVI的关系。

TVDI值为1是干边( Dry edge), 代表土壤缺水；TVDI 值为0则是湿边( Wet edge)，具有最大的土壤蒸发蒸腾总量和无限的水分供应，反映了土壤水分的两个极端状态。

Sandholt等[14] 首先提出了温度植被干旱指数( Temperature Vegetation Dryness Index, TVDI )的概念。

其表达式为:TVDI=(Ts-Tsmin)-(Tsmax-Tsmin)
本实验只统计NDVI对应下的Tsmin和Tsmax。

NDVI只选择>0的部分。

得到干边方程为：
Tsmax= -36.779*NDVI+323.55,R方值为0.9564
Tamin= 65.614*NDVI+268.99，R方值为0.9621
代码：pro TVDI
;读取NDVI
fn=pickfile(title='��NDVI�ļ�')
openr,lun,fn,/get_lun
ndvi=fltarr(2205,1814)
readu,lun,ndvi
tvscl,ndvi,order=1
free_lun,lun
;读取Temperatu rֵ
fn1=pickfile(title='Temperature')
openr,lun1,fn1,/get_lun
wendu=fltarr(2205,1814)
readu,lun1,wendu
tvscl,wendu,order=1
free_lun,lun1
;设置ndvi步长及组数
a=max(ndvi)
b=min(ndvi)
step=0.01
n=ceil((a-b)/step)
;
mini=fltarr(n)
maxi=fltarr(n)
;ֵ
for i=0,n-1do begin
mini[i]=b+i*step
maxi[i]=mini[i]+step
endfor
minimum=fltarr(n)
maximum=fltarr(n)
for i=0,n-1do begin
temp1=where(ndvi gt mini[i] and ndvi lt maxi[i])
if n_elements(temp1) eq1then begin
temp2=0
endif else begin
temp2=wendu[temp1]
endelse
minimum[i]=min(temp2)
maximum[i]=max(temp2)
endfor
plot,mini,minimum,psy=2
oplot,maxi,maximum,psym=2
data=fltarr(3,n_elements(mini))
data[0,*]=transpose(mini)
data[1,*]=transpose(minimum)
data[2,*]=transpose(maximum)
openw,lun2,'d:/TVDI000.txt',/get_lun
printf,lun2,data,format='(f10," ",f20," ",f20)'
free_lun,lun2
END
如图：
⑤TVDI指数计算。

输入公式为：
(B1-(65.614*B2+268.99))/((-36.779*B2+323.55)-(65.614*B2+268.99))
⑥云检测。

本实验判定Ts<290K的部分为云，因此做一个云腌膜文件，Data min设为290，非云为1，云为0。

⑦水体检测。

本实验判定MNDWI>0.40的部分为水体，因此做一个水体腌膜文件，Data max设为0.40，非水为1，水为0.
⑧NDVI<0检测。

原理同上，去除建筑物。

⑨得到最后的TVDI图像。

输入公式：
B1*B2*B3*B4，B1…B4依次对应上面⑥、⑦、⑧及TVDI图像。

⑩TVDI密度分割，制图输出，步骤同4.1⑥。

如图：
4.3土壤湿度遥感反演。

①土壤站点矢量图生成。

在Arcmap中加载X/Ydata，选择实测土壤湿度文件.txt，字段要对应好，此处不详细介绍，坐标系选择WGS-1984,，生成图层，导出shapfile文件。

②ENVI打开矢量图层，生成evf文件。

加载生成的evf文件到TVDI窗口图像中，在“可以使用的矢量列表中”选择导出evf图层到ROI，属性文件为TVDI.
③获取土壤湿度对应经纬度下的TVDI值。

在ENVI下启动ROI工具，导入生成的ROI文件。

导出ROI文件到ASCII文件，属性文件依旧选择TVDI.
④将生成的文本文档和实测土壤湿度文本文档用EXCEL打开，按空格顺序打开。

按照经度升序排列，选取前70对数据建立模型，拟合模型得到线性模型为y=-0.8225*x+0.8375 ,R方值为0.6857，X为TVDI，Y为土壤湿度。

70对数据相关系数为-0.82805，查表知，其通过了0.001的显著性检验，置信度为99.9%，可以使用。

⑤用剩余的20对数据作为模型检验的数据。

代入上述公式，得到遥感反演土壤湿度数据，计算得MSE为0.243569，RMSE为0.091851.
⑥基于公式：-0.8225*B1+0.8375 对TVDI图层进行土壤湿度反演。

基于公式(B1 eq 0.8375)*0+(b1 ne 0.8375)*B1做云水体建筑物的掩膜。

进行密度分割，得到结果。

5.实验总结
通过本次实验，认识了常用的遥感图像格式，掌握了遥感图像几何校正、研究区裁切、NDVI、TVDI 、土壤湿度反演计算等流程，深入领会到了图像密度分割以及遥感制图技术方法，收获颇多，感闵爱莲老师对本次学生实验中的指导！。