（五）Landat_5 TM 遥感影像计算NDVI、MNDWI、NDBI以及地表温度反演

流程概况图

一、植被覆盖指数(NDVI)

所谓植被指数，就是利用多波段遥感图像的可见光波段以及近红外波段相组合，形成能够反映突出植被分布状况和强度的指数。植被指数的求解方法不唯一，定义的植被指数已经多达50多种。通过红波段以及近红外波段，可以有效的消除大气对植被指数的影响。
公式：
NDVI=（band4-band3）/（band4+band3）=（NIR-RED）/（NIR+RED）

%计算植被覆盖指数NDVI
NIR=double(image{4});
RED=double(image{3});
NDVI=(NIR-RED)./(NIR+RED);
imshow(NDVI);
c_map=zeros(100,3);
c_map(1:12,:)=repmat([0.5,0.25,0],12,1);
c_map(13:24,:)=repmat([0.35,0.25,0],12,1);
c_map(25:36,:)=repmat([0.1,0.25,0],12,1);
c_map(37:48,:)=repmat([0,0.2,0],12,1);
c_map(49:61,:)=repmat([0,0.4,0],13,1);
c_map(62:74,:)=repmat([0,0.65,0],13,1);
c_map(75:87,:)=repmat([0,0.8,0],13,1);
c_map(88:100,:)=repmat([0.5,1,0],13,1);
colormap(c_map);
colorbar;
title('植被指数---NDVI((NIR-RED) / (NIR+RED))');

二、归一化水体指数(MNDWI)

该归一化水体指数MNDWI是通过NDWI修改之后以增强开放水域要素，通过绿波段和中红外波段相组合，从而使整幅图的水域均衡显示。我之前正在做实验时，用到另外两个归一化水体，一个提取不充分，另一个冗余度太高，只有此归一化指数，在提取显示方面具有均衡性的特点。
公式：
MNDWI=（band2-band3）/（band2+band3）=（GREEN-MIR）/（GREEN+MIR）

%计算水体归一化指数NDWI
GREEN=double(image{2});
MIR=double(image{5});
MNDWI=(GREEN-MIR)./(GREEN+MIR);
MNDWI=adapthisteq(MNDWI.*100);
imshow(MNDWI);
title('归一化水体指数---MNDWI((GREEN-MIR)/(GREEN+MIR))');

三、归一化建筑指数

归一化建筑指数是利用近红外和中红外波段结合，从而突出显示建筑区域，NDBI图像颜色越深，代表值越大，表明建筑用地比例越高，建筑密度越高，通过目视解译，不断调整NDBI阈值，直到找到合适的阈值，提取建设用地，调整阈值后获得不透水面分布 。
公式：
NDBI=（band5-band4）/（band5+band4）=（MIR-NIR）/（MIR+NIR）

%计算归一化建筑指数
NDBI=(MIR-NIR)./(MIR+NIR);
NDBI=adapthisteq(NDBI.*100);
imshow(NDBI);
title('归一化建筑指数---NDBI((MIR-NIR)/(MIR+NIR))');

四、地表温度反演

地表温度反演不同于其他指数那样，直接利用波段进行结合计算得到，在反演地表温度之前，需要先计算植被覆盖率、地表比辐射率、黑体辐射亮度等值，我计算的这片区域，温度基本都在26~28摄氏度之间，高温主要聚集在城市裸地区域，低温主要聚集在河流植被区域。
公式：
T=（K2/ ln（K1/Dt）+1）-273
K1和K2的值可以在头文件中获取：

在计算Dt时，大气在热红外的透过率、大气向下辐射亮度、大气向上辐射亮度需要去官网进行查询。

%>>>>>>>>>>>>>>>>>>>>>>>>>>计算植被覆盖率PV>>>>>>>>>>>>>>
w1=NDVI;
w1(w1<0.03)=0;
w1(w1>0.35)=1;

w2=w1;
w2(w2==1)=0;
PV=w1+w2.*((NDVI-0.03)/(0.35-0.03));
%>>>>>>>>>>>>>>>>>>>>>>>>>计算地表比辐射率x>>>>>>>>>>>>>>
x=0.004.*PV+0.986;
%>>>>>>>>>>>>>>>>>>>>>>>>>计算黑体辐射亮度>>>>>>>>>>>>>>>
ca=0.74;%大气在热红外的透过率
cb=3.45;%大气向下辐射亮度
cc=2.07;%大气向上辐射亮度
b2=double(fs_image{6});%热红外波段的辐射标定值
Dt=(b2-cc-ca.*(1-x).*cb)./(ca.*x);
%>>>>>>>>>>>>>>>>>>>>>>>>>>>计算地表温度>>>>>>>>>>>>>>>>>>>
K1=table2array(head_data(173,"L1_METADATA_FILE"));
K2=table2array(head_data(174,"L1_METADATA_FILE"));
T=K2./log((K1./Dt)+1)-273;
%剔除异常值
T(T>40)=0;
%为了方便出图，上面的T为真实地表温度值
TT=(T-27).*0.1;
imshow(TT);
colormap("hot")
title('地表温度反演值---(K2/ln(K1./DT)+1)-273');