ENVI大气校正方法反演Landsat 7地表温度
本文介绍基于ENVI软件，实现对Landsat 7遥感影像加以大气校正方法的地表温度反演操作。
1 图像前期处理与本文理论部分
更新：基于GEE的地表温度Landsat反演可以看文章单窗算法的地表温度反演：谷歌地球引擎GEE代码，自动批量操作，处理更快。
本文操作部分将直接由植被覆盖度计算展开；而对于一个完整的地表温度反演计算过程，在求解植被覆盖度这一步骤之前仍有很多数据准备、预处理等工作。为了更好理解整个实验过程，将我们未进行的步骤梳理如下。其中，具体的前期操作方法大家可以参考文章ENVI实现QUAC、简化黑暗像元、FLAASH方法的遥感影像大气校正。
首先获取需要的遥感图像数据，并对其进行包括上述文章内容在内的预处理步骤——数据导入、辐射定标、几何校正、大气校正、图像拼接与裁剪等。
其中，辐射定标需要分两步骤，即对可见光波段数据（如上述文章中的1、2、3、4、5、7六个波段）与热红外波段数据（如上述文章中的6波段）分别进行辐射定标。
其次，需要计算NDVI（即Normalized Difference Vegetation Index，归一化植被指数，而非植被覆盖度）。NDVI是指一幅遥感影像中，近红外波段的反射值与红光波段的反射值之差比上这二者之和；其可以用来检测植被生长状态、植被覆盖度，还可以消除部分辐射误差等。NDVI的具体取值范围限制在-1到1之间，其负值表示地面覆盖为云、水、雪等，即对可见光具有高反射；0值表示地面覆盖有岩石或裸土等，从而使得NIR（Near Infrared，近红外波段）和R（Red，红光波段）近似相等；其正值则表示地面有植被覆盖，且植被覆盖度越高，其数值越高。目前，在一些网站（如NASA官方网站）具有NDVI成品数据，可供我们直接下载、利用；而通过初始遥感影像中的近红外波段数据和红光波段数据，我们可以直接利用前述定义公式，即
对其加以计算。计算NDVI时需要注意，所选用的遥感影像不能具有过多的云干扰。
再次，需要视实际情况对计算得出的NDVI结果图像进行重采样。这是由于，在本文中需要多次利用“Band Math”工具对图像数据进行计算，而这一工具要求输入的数据在图像分辨率（即像元大小）及行、列数等方面完全一致。同时，我们即将使用的热红外数据（即Landsat ETM+第六波段数据），其分辨率为60米；而计算得到的NDVI数据图像分辨率为30米。因此，我们需要对分辨率精度更高的NDVI数据图像加以重采样处理，使得二者分辨率一致。
重采样功能可以通过ENVI软件中选择“Basic Tools”→“Resize Data (Spatial/Spectral)”加以实现。
目前主要的地表温度单波段反演算法包括大气校正法（又名辐射传输方程法，Radiative Transfer Equation，RTE）、单通道算法和单窗算法。本文我们使用大气校正方法。大气校正法的基本原理为：估计得到大气对地表热辐射的影响，然后将这一部分大气影响由卫星传感器所接收到的热辐射总量中减去，得到地表热辐射强度；最后将地表热辐射强度转化为对应的地表温度即可。
卫星传感器接收到的热红外辐射亮度值Lλ由三部分组成，其分别为：大气向上辐射亮度L↑、地表物体的真实辐射亮度在经过大气层后到达卫星传感器的能量、大气向下辐射亮度在经过地面反射后的能量。因此，结合上述这一理论过程，可以用辐射传输方程来表示卫星传感器接收到的热红外辐射亮度值：
其中，ε为地表比辐射率，T_s为地表真实温度，B(T_s )为地表在T_s这一真实温度下的黑体热辐射亮度，τ为大气的热红外波段透过率。地表比辐射率ε又称为发射率，而根据基尔霍夫定律，发射率与吸收率相等，则(1-ε)可以表示反射率。因此，[εB(T_s )τ]即为地表物体的真实辐射亮度在经过大气层后到达卫星传感器的能量，而[(1-ε)L↓τ]则表示大气向下辐射亮度在经过地面反射后的能量。
物体的比辐射率是物体向外辐射电磁波的能力表征，是指在同一温度下地表发射的辐射量与一黑体发射的辐射量的比值。其不仅依赖于地表物体的组成，而且与物体的表面状态（如表面粗糙度）及物理性质（如介电常数、含水量等）有关，并随着所测定的波长和观测角度等因素变化。对地表比辐射率的精确定量测量难度较大，因此本文依据经验法对地表比辐射率加以估计：
其中，F_C为植被覆盖度。由此可知，我们需要同时计算出地表植被覆盖度，用以确定地表比辐射率计算公式并参与计算。本文我们采取混合像元分解法求解植被覆盖度。同上述地表比辐射率计算公式较为一致，我们依然将地表分为水体、植被与建筑三个部分；其中，依据NDVI数值对这一类别加以具体区分。NDVI小于0时，认为地物为水体，植被覆盖度为0；NDVI大于0.7时，认为地物为植被，植被覆盖度为1；NDVI取值在[0,0.7]时，认为地物处于水体与植被之间，植被覆盖度依据公式计算。针对不同的地物，计算植被覆盖度：
最后，完成上述全部计算并依据辐射传输方程求得B(T_s )后，可以依据普朗克公式反函数求出地表真实温度。公式为：
所得地表真实温度单位为开尔文（K），我们需要将其转换为常见的摄氏度。温度单位的转换可以在得到地表真实温度图像后单独计算，亦可以直接在上式中计算。
2 实际操作
2.1 植被覆盖度计算
依据上述分析，首先我们需要借助60米分辨率的研究区NDVI数据图像，计算研究区的植被覆盖度。
（1） 打开ENVI Classic 5.3（64-bit）软件，选择“File”→“Open Image File”，在弹出的文件选择窗口中选择“TM-NDVI-60m.img”文件；点击“打开”。
（2） 选择“Basic Tools”→“Band Math”，在弹出的公式创建窗口中输入本次实验的第一个公式，即植被覆盖度公式。输入公式完成后，点击下方“Add to List”按钮，即可将公式存入待选择区内。为了减少后期不必要的工作量，可以每次编辑完成一个公式后点击“Save”按钮，以将待选择区内的公式保存。
（3） 保存公式完成后，点击左下角“OK”按钮，即可开始公式的计算。在弹出的公式变量文件选择窗口中，将这一公式的变量“B1”选择为我们刚刚添加的文件“TM-NDVI-60m.img”。随后，配置输出文件地址等信息。
（4） 配置完成后，点击左下角“OK”按钮，即可开始公式的运行。运行结束，将所得到的研究区植被覆盖度结果图像导入ENVI软件中，显示如下。
2.2 地表比辐射率计算
如前所述，我们完成植被覆盖度的计算其实是为地表比辐射率计算做准备。同样是将地物分为水体、城镇与自然表面三个类别，依据经验法，运用与植被覆盖度形式类似的分段计算公式，为不同类别地物赋予不同的地表比辐射率计算公式。
考虑到在计算植被覆盖度时使用的公式中已包含“B1”这个变量波段名称，为了避免不同公式之间的变量相互混淆，因此在一开始我在每一个公式中都使用了不同的变量名称——如在此处地表比辐射率计算公式中使用“B2”“B3”等。但通过后期的实验过程发现，其实每一次公式的变量名称即使一致也不会对实验造成明显的影响。
（1） 选择“Basic Tools”→“Band Math”，在弹出的公式创建窗口中输入本次实验的第二个公式，即地表比辐射率公式。输入公式完成后，点击下方“Add to List”按钮，即可将公式存入待选择区内。为了减少后期不必要的工作量，可以在编辑完成这一公式后点击“Save”按钮，并将最近一次保存的公式文件覆盖，以将待选择区内的两条公式统一保存。
（2） 保存公式完成后，点击左下角“OK”按钮，即可开始公式的计算。在弹出的公式变量文件选择窗口中，将这一公式的变量“B2”依然选择为我们最开始添加的图像文件“TM-NDVI-60m.img”，并将这一公式的变量“B3”选择为通过上述步骤获得的植被覆盖度结果图像文件。随后，配置输出文件地址等信息。
（3） 配置完成后，点击左下角“OK”按钮，即可开始公式的运行。运行结束，将所得到的研究区地表比辐射率结果图像导入ENVI软件中，显示如下。
2.3 相同温度下黑体辐射亮度值计算
正本文第一部分分析所得，已知研究区域地表比辐射率与热红外波段亮度，我们便可以计算相同温度下黑体辐射亮度值。
（1） 在ENVI Classic 5.3（64-bit）软件中选择“File”→“Open Image File”，在弹出的文件选择窗口中选择“TM6-rad-subset-jz-xiangfan.img”文件；点击“打开”。
（2） 选择“Basic Tools”→“Band Math”，在弹出的公式创建窗口中输入本次实验的第三个公式，即黑体辐射亮度值公式。输入公式完成后，点击下方“Add to List”按钮，即可将公式存入待选择区内。为了减少后期不必要的工作量，可以在编辑完成这一公式后点击“Save”按钮，并将最近一次保存的公式文件覆盖，以将待选择区内的三条公式统一保存。
（3） 保存公式完成后，点击左下角“OK”按钮，即可开始公式的计算。在弹出的公式变量文件选择窗口中，将这一公式的变量“B4”选择为我们刚刚添加的图像文件“TM6-rad-subset-jz-xiangfan.img”，并将这一公式的变量“B5”选择为通过上述步骤获得的地表比辐射率结果图像文件。随后，配置输出文件地址等信息。
（4） 配置完成后，点击左下角“OK”按钮，即可开始公式的运行。运行结束，将所得到的研究区相同温度下黑体热红外波段辐射亮度值计算结果图像导入ENVI软件中，显示如下。
2.4 地表真实温度计算
由前述分析可知，通过上述步骤得到的黑体热红外波段辐射亮度并不是地表实际温度，我们依然需要通过普朗克公式反函数实现二者之间的转换关系。
与此同时，这一步骤得到的地表实际温度结果单位为开尔文（K），并不是我们平日里经常使用的摄氏度（℃）。因此，我们还需要实现温度单位的转换。
由开尔文温度转换为摄氏度只需要在原温度基础之上减去273.15即可，较为简单，没有必要单独转换。因此我选择直接在这一步骤将温度单位的转换完成。
（1） 选择“Basic Tools”→“Band Math”，在弹出的公式创建窗口中输入本次实验的第四个公式，即地表真实温度公式。输入公式完成后，点击下方“Add to List”按钮，即可将公式存入待选择区内。为了减少后期不必要的工作量，可以在编辑完成这一公式后点击“Save”按钮，并将最近一次保存的公式文件覆盖，以将待选择区内的四条公式统一保存。
（2） 保存公式完成后，点击左下角“OK”按钮，即可开始公式的计算。在弹出的公式变量文件选择窗口中，将这一公式的变量“B6”选择为我们刚刚计算获得的黑体热红外波段辐射亮度结果图像文件，随后，配置输出文件地址等信息。
（3） 配置完成后，点击左下角“OK”按钮，即可开始公式的运行。运行结束，将所得到的研究区地表真实温度计算结果图像导入ENVI软件中，显示如下。
（4） 在得到的结果图像任意位置处右键，选择“Cursor Location/Value”或“Quick Stats”选项，均可以查看图像各个像元的像素信息。在这里，我使用“Cursor Location/Value”查看图像像元信息，得到如下所示的结果。
（5） 可以看到，其中像元像素值（即“Data”值）出现了“302.902771”这一数值。而这一步骤所得到的结果为单位为摄氏度的地表实际温度数值，不可能出现三百多的数据。因此，说明很可能在前面操作部分出现错误。
（6） 通过返回检查发现，在其中计算地表比辐射率时，所输入的公式出现了输入错误。将这一错误纠正后，重新生成地表实际温度图像，并对其进行统计检查。
（7） 可以看到，经过修改后的地表真实温度图像数据符合实际情况，可以认为错误已被排除。
2.5 图像导出
由于在实验后期需要制作专题地图，以将温度分为不同等级并比较不同地物的温度特性，因此需要将ENVI中得到的地表实际温度结果图像另存为“.tif”格式，从而方便利用ArcGIS系列软件对其加以进一步分级、美化与出图处理。
（1） 首先，我在ENVI图像窗口中进行图像的另存为操作。选择“File”→“Save Image As”→“Image File”，在弹出的保存配置窗口中选择文件格式为“TIFF/GeoTIFF”，并配置好结果图像文件保存路径、保存文件名等。
（2） 得到保存的结果图像文件后，将其添加进入AcrMap 10.2软件中。此时发现，得到的图层文件在像素上已被拉伸至0-255范围。通过AcrMap的“识别”功能，发现像元像素已全部被拉伸，即其原值均已成为上述0-255范围的数据。因此，认为图像应该不会再通过相关操作恢复原有的温度数值。随后，尝试用“设置栅格属性”等工具加以调整，均以失败告终。
（3） 通过查阅网络资源，看到有人指出这种情况是保存了“图像”而不是“影像”；在这一提醒下，尝试直接利用ENVI软件左上角工具栏中的“File”→“Save File As”→“TIFF/GeoTIFF”对图像加以保存。在弹出的保存配置窗口中选择文件，并配置好结果图像文件保存路径、保存文件名等。
（4） 将第二次保存的“.tif”格式图像结果文件导入AcrMap 10.2软件中，看到结果数值为正常状态。
（5） 得到正常结果后，对两种不同的保存结果加以思考。本文所使用的软件为ENVI Classic 5.3（64-bit）版本，而使用非Classic版本同样发现这一问题。回顾在图像界面保存图片的过程，进一步发现，这样的保存应当是将图像格式直接设置成为了8bit整形存储的格式，从而丢失其原有的像素数据信息与意义——即其只保存了这幅图像的“外貌”，运用0-255之间并无实际意义的数据来表示图像的灰度；但其原有图像，即我们希望得到的图像，其数据、灰度等级应是温度的表示。若将这张原有由不同温度划分出不同灰度等级的图像改变为由0-255之间数字划分出不同灰度等级的图像，其自然失去了原有的温度意义。
（6） 将正常图像导入AcrMap 10.2软件中后，同时发现其影像的黑边无法去除——运用“识别”功能可以看到，黑边部分并无数据，其原本就已均为“NoData”状态。针对这一现象，可能是不同软件在保存、读取数据文件时出现的常见错误，无需在意，后期对温度重新划分等级后即可消失。
（7） 最终导入图像及其信息如下图所示。其中，有两个结果相近的图层是因为在一开始地表实际温度转换过程中，我减去的数值为“273”而并非“273.15”。
3 专题地图制作
导入AcrMap 10.2软件中的图像为像元像素由22.3997至47.0569，若不对其重新划分温度区间处理，将会导致所成彩色图像各个颜色十分离散，不利于观察、分析。
（1） 导入后的图像自身不具有统计信息，无法在“已分类”模块中对其各个数值划分区间。因此，需要首先使用工具对其数据信息加以统计。点击“Data Management Tools.tbx”→“栅格”→“栅格属性”→“计算统计数据”。
（2） 在弹出的配置窗口中选择需要统计数据的图层，其它项目不需要做调整。
（3） 点击“确定”，即可对图层图像完成统计数据的计算。
（4） 在图层名称处右键，选择“属性”→“符号系统”→“已分类”，打开“分类”按钮，将分类方法选择为“手动”，并依据要求设置四个“中断值”——30，35，39和图像像素的最大值47.0569。点击确定即可保存间断点设置。
（5） 回到图层属性窗口后，对四个范围所表示的内容加以标注，以方便后期制图。同时，针对温度专题图的特色，调节合适的配色方案。原本我准备使用红色系作为配色，但发现出图后整体区分效果并不是很好。因此选择另一种对比相对较强的配色方案。
4 不同地物地表温度对比分析
为了更加清晰地对比不同地物地表温度的差异，借助腾讯地图中卫星地图模块，以襄阳市汉江为对照，对比研究区域卫星地图与专题地图。得出结论如下：
（1） 水体温度明显低于周围城市、乡村等区域。这一特征可谓温度专题地图中最为明显的要素之一——由专题地图可以清晰识别出一条低温带由襄阳市西部进入，由西北向东南方向延申；并在城市中部偏西南方向形成冲击岛，随后继续向南行进。此外，在冲击岛北部，由唐白河形成的低温带同样较为明显；其清晰可辨，甚至可以沿低温带追溯至该河流的发源处。城市东南处的秦咀水库同样如此。由水体的温度色带可以看出，在成像时刻，其温度在30℃以下。由于水体的温度和光照有密切关系，尤其是汉江、唐白河这一类地上河；而本次实验所用图像的成像时间为上午10:30左右，当天日照时间还不长，因此水温较低。
（2） 城镇温度明显高于周围区域。由温度专题地图可以看出，主要的高温聚集区域均位于汉江、唐白河沿岸；结合卫星地图可以看出，这些大多数均为襄阳市主要城市聚集中心，建筑、人口密集，经济发达。由城镇温度色带可以看出，在成像时刻，其温度在39℃以上。在另一方面，对于一些建筑物分布较为零散的郊区、农村等，温度相对没有城市中心那么高；但其温度普遍在35-39℃区间内，较之周围区域温度还是较高。由于城镇人口多，汽车、生活燃烧等排放较多，加之工业生产频繁，因此温度较高。
（3） 农田、耕作用地等温度较高。35℃以上区域主要分布在汉江、唐白河沿岸的另一原因是这些耕作用地和城镇一样，多分布于地势平坦、土壤肥沃、灌溉方便的河流两岸。由专题地图可以看出，多数耕作用地的温度在35-39℃这一区间；当然，同样亦有不小面积的农田温度在30-35℃区间。因此，可以认为耕作用地温度处于35℃左右。其温度较之城镇低，是由于乡村地广人稀，人口密度显著小于城市地区，因此产生的热量较少；而其温度较高于森林、山区则是因为地形平坦、人为耕作调节土温等。
（4） 森林、山区等温度明显较低。在温度专题地图中，有三处主要的低温聚集中心，分别为襄阳市西南部位、西部汉江北岸地区、中部汉江东部地区。其中最明显、面积最大的低温聚集中心——城市西南方的低温处为山区，而其他两处区域则多为森林。可以看到，这些低温地区的温度都在30℃以下。由于森林中树木的遮荫、吸热作用，使得森林地表温度明显低于其他地区；而对于山区而言，由于气压低，空气稀薄，大气保温效果较差；同时高山表面距离地面较远，无法吸收地面热量，使得其表面温度较低。与此同时，二者同样具有人口密度很低的特点。
