遥感建模与应用综合实习实验报告学期2017-2018学年第二学期姓名学号指导教师闵爱莲1实验题目:土壤湿度遥感反演研究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图像、BASEIMAGE图像、研究区矢量边界图。4.实验步骤:4.1图像预处理以及NDVI制图。①大气校正。本实验选择黑暗像元法进行大气校正。打开TM200205,快速统计图像DN值,记录7个波段每个波段的黑暗像元。本实验选择DN值像元个数为十位数时对应的DN值作为黑暗像元DN值。所有像元均减去此黑暗像元的DN,达到校正目的。所依据原理为认为黑暗像元的辐射记录值为0,且认为图像各处大气影响一致,之所以不为0是由于大气程辐射的影响,减去此DN,即为整体消除大气程辐射的影响。②几何校正。加载并打开baseimage和TM200205图像,基于Map下的图2像到图像校正方法进行几何校正。以baseimage为基准图像,以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字体,输出如图:34.2热红外地温反演及TVDI制作(以下所有公式均在BandMath下输入)。①单位光谱范围内幅亮度值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的实时性强,却易受到周围地物的影响,因此常将两者综合起来研究干旱。4Price和Carlson等研究首先发现,如果研究区的植被覆盖度和土壤湿度值变化较大,则获得的植被指数NDVI和地表温度Ts所形成的散点图呈三角形关系,其中NDVI作为横坐标,Ts作为纵坐标。Moran等则分析NDVI和Ts的散点图为梯形关系。Sandholt等通过研究简化的NDVI-Ts三角形空间,提出了温度植被干旱指数TVDI,用以估算土壤含水状况。NDVI-Ts特征空间的示意图,体现了Ts与NDVI的关系。TVDI值为1是干边(Dryedge),代表土壤缺水;TVDI值为0则是湿边(Wetedge),具有最大的土壤蒸发蒸腾总量和无限的水分供应,反映了土壤水分的两个极端状态。Sandholt等[14]首先提出了温度植被干旱指数(TemperatureVegetationDrynessIndex,TVDI)的概念。其表达式为:TVDI=(Ts-Tsmin)-(Tsmax-Tsmin)本实验只统计NDVI对应下的Tsmin和Tsmax。NDVI只选择0的部分。得到干边方程为:Tsmax=-36.779*NDVI+323.55,R方值为0.95645Tamin=65.614*NDVI+268.99,R方值为0.9621代码:proTVDI;读取NDVIfn=pickfile(title='��NDVI�ļ�')openr,lun,fn,/get_lunndvi=fltarr(2205,1814)readu,lun,ndvitvscl,ndvi,order=1free_lun,lun;读取Temperaturfn1=pickfile(title='Temperature')openr,lun1,fn1,/get_lunwendu=fltarr(2205,1814)readu,lun1,wendutvscl,wendu,order=1free_lun,lun1;设置ndvi步长及组数a=max(ndvi)b=min(ndvi)step=0.01n=ceil((a-b)/step);mini=fltarr(n)maxi=fltarr(n);fori=0,n-1dobeginmini[i]=b+i*stepmaxi[i]=mini[i]+stependforminimum=fltarr(n)maximum=fltarr(n)fori=0,n-1dobegintemp1=where(ndvigtmini[i]andndviltmaxi[i])ifn_elements(temp1)eq1thenbegintemp2=06endifelsebegintemp2=wendu[temp1]endelseminimum[i]=min(temp2)maximum[i]=max(temp2)endforplot,mini,minimum,psy=2oplot,maxi,maximum,psym=2data=fltarr(3,n_elements(mini))data[0,*]=transpose(mini)data[1,*]=transpose(minimum)data[2,*]=transpose(maximum)openw,lun2,'d:/TVDI000.txt',/get_lunprintf,lun2,data,format='(f10,,f20,,f20)'free_lun,lun2END如图:⑤TVDI指数计算。输入公式为:(B1-(65.614*B2+268.99))/((-36.779*B2+323.55)-(65.614*B2+268.99))⑥云检测。本实验判定Ts290K的部分为云,因此做一个云腌膜文件,Datamin设为290,非云为1,云为0。Tsmin=65.614NDVI+268.99R²=0.9621Tsmax=-36.779NDVI+323.55R²=0.9564250.0275.0300.0325.0350.00.000.050.100.150.200.250.300.350.400.450.500.55TSNDVITs-NDVI散点图TminTmax线性(Tmin)线性(Tmax)7⑦水体检测。本实验判定MNDWI0.40的部分为水体,因此做一个水体腌膜文件,Datamax设为0.40,非水为1,水为0.⑧NDVI0检测。原理同上,去除建筑物。⑨得到最后的TVDI图像。输入公式:B1*B2*B3*B4,B1…B4依次对应上面⑥、⑦、⑧及TVDI图像。8⑩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.9⑥基于公式:-0.8225*B1+0.8375对TVDI图层进行土壤湿度反演。基于公式(B1eq0.8375)*0+(b1ne0.8375)*B1做云水体建筑物的掩膜。进行密度分割,得到结果。5.实验总结通过本次实验,认识了常用的遥感图像格式,掌握了遥感图像几何校正、研究区裁切、NDVI、TVDI、土壤湿度反演计算等流程,深入领会到了图像密度分割以及遥感制图技术方法,收获颇多,感谢闵爱莲老师对本次学生实验中的指导!