【Envi】基于单窗算法的地表温度反演实验操作记录

比前言还前

2021年7月14日：虽然这一年一直在做遥感的工作，但是由于研究尺度和效率的问题我主要使用python和gee，所以这篇可能是有关Envi的唯一一篇博客。这篇博客写于刚学习Envi的时候，其内容是唯一一次通过Envi完成的工作（很基础的工作），所以我对于Envi的理解相当浅薄，所以这篇博客内容存在着一些问题。评论区一些小伙伴们已经指出了一些问题，有一些我实在是不懂得解决或者解决不了（或者没看到）。所以不要盲目复制工作，特别是在做研究的时候。复制这篇博客的工作的时候，最好还是看一遍评论区。

前言

前天刚下的Envi，高强度学习了两天摸索出了点门道，不一定准确。本文谨作为个人记录，如有问题欢迎指正。
写这篇博客的时候，我是边做边写的，想到啥写啥，可能会比较混乱一点。后面有空的话我会整理一下。

参考博客

在实践操作这方面，相关博客给的启发更大一点；而文献大部分只参考了公式或者是经验数值，太多了太杂了就不一一写出来了。
记 / ENVI5.3温度反演-单窗算法-操作要点
地表温度热红外遥感反演－以单窗算法为例
地表温度反演操作总结

技术流程

1.数据预处理

1.1数据获取

首先在地理空间数据云下载landsat8的系列影像。
因为我的研究区是长春，覆盖范围是三幅行列号分别是119029、118029、1118030，其中118029是主图。

1.2辐射亮度温度

1.2.1热红外波段辐射定标与亮度

经过一系列摸索，我喜欢先做热红外波段的辐射定标，因为可以在一定程度上遇见最终的温度反演输出结果。
首先导入图像，进行辐射定标，波段选取热红外Thermal。

然后下面的spectral subset只勾选上面那个（10.9000），然后点击确定。

进入下一个页面，上面的部分保持默认，只需要填output filename就行了。我的命名规则是行列号+rhwfsdb（热红外辐射定标）（很土我知道）。然后确定，就开始允许了。

辐射定标之后它会自动加载到图层里，完事就得到一张热红外波段定标后的图层（如下）。
可以通过鼠标值先查看一下结果是否靠谱。
landsat8在这一步处理完后数值一般在5-11之间，温度越高的地方可能就越高。
所以由于热岛效应，城市的值会更高。
但是我之前选了个3月初的长春的遥感图，反演出来的结果不太理想，那幅图的结果显示的温度主要呈现出西高东低，城市内的温度不是特别显著，我后来思考了一下，可能三月初这个时间段正在化雪，所以就导致温度不太准确。
吸取教训，这次我选择的是九月底的遥感图像。

处理完第一幅图后，对其他两幅图也进行相同的操作。
把三幅图放在一起看的效果长这样（下图）。为了方便观察，导入事先准备好的长春市shp。嗐，本来两幅图就行了，偏偏左边多出了个边边，还我要用三幅图。

再然后可以准备开始拼接图像了。但是在拼接之前，经过我多次实验发现，它这个黑色的框框会影响拼接（可能只有这三幅图会），所以我决定用ROI把它裁一下，特别是119029这幅图，我只用到它的一小块，所以更需要裁一下，而且还能提高拼接速度（虽然没提高几秒）。
首先裁下118029这幅图，右键要裁的图层选ROI调出窗口。为了避免裁多了，我还是把长春市的shp放着作参考。

用ROI分别裁剪完三幅图片，我这边是这个样子的（下图）。

然后可以正式开始拼接了。搜索seamless mosaic，导入全部图层。

他会选定一幅作为reference，另外两幅根据reference进行调整（好像是这样）。而且这个reference好像是默认选定你打开到envi的第一幅图，使用seamless mosaic后好像更改不了（如果可以更改的话请教教我）（我这段话是不是太多“好像“了）。

color correction栏里，勾选histogram matching。

然后再点进seamlines里的auto generate seamlines，它就开始生成缝合线了（原谅我这里截图不了）。

其他保持默认。点进export，选择导出路径，就可以开始导出了。我把它命名为rhwfsdb_mosaiced。

完事之后图长这个样子（下图）。我顺便把长春shp导进来参考一下，完美！

然后为了裁出长春的形状，用感兴趣区裁剪，命名为CC_rhwfsdb。

裁剪结果（下图）。

然后我把他配色变换一下，已经能看得出有一点像了（比我之前做的那副三月份的效果好多了）。
虽然已经能看出点东西，但是这还是最简单的一步。

1.2.2辐射亮度温度计算

需要用到这个公式，Ti是辐射亮度温度（辐射亮温），Lλ是上面获取的辐射亮度（就是值大概在5-11之间的那幅图）。

关于这个K1和K2，查了一堆参考文献，基本能达成共识。K1=774.89，K2=1321.08。
然后就可以开始我们的波段运算了。
输入公式：1321.08/alog(774.89/b1+1)

下一个页面它让你选b1是啥，选择我们刚刚得到的CC_rhwfsdb，就是还是上面那图。命名为CC_fslw。
计算完成，得出一幅值为300左右的图，这个就是它的辐射亮温了，单位是开尔文。

1.3地表比辐射率

这一部分主要获取植被覆盖度图片。具体的流程是先对三幅多光谱的遥感图像分别做辐射定标，然后进行拼图，再然后计算ndvi，再再然后利用ndvi计算植被覆盖度。

1.3.1多光谱辐射定标

首先进行辐射定标。
老操作，不过这次要选的是戴MultiSpectral的，然后点击确定。

在下一个页面里，output interleave换成BIL，然后点击apply flaash settings把scale factor换成0.10，在output filename设置输出，我这里的命名规则是行列号+fsdb。最后点击确定，等待处理完成。

按照以上的操作把全部三幅图的多光谱都辐射定标处理完成后，放在一起大概长这样（下图）。

然后开始拼图操作，和前文的操作一样；拼图完后，用长春shp把研究区域裁下来，这里就不细说了。最后输出的命名为CC_mosaiced。得到下面这样式的图。

1.3.2NDVI计算

把前文多光谱辐射定标+拼接+裁剪后得到的图进行ndvi处理。直接在工具栏里搜索ndvi，有集成好的工具直接可以用了。
直接选定需要处理的图层，其他默认，点OK。

下一个页面，也是其他默认，然后选择输出，我直接命名为ndvi。静待它ndvi计算完成。

ndvi计算结束后得到这样一幅图，这个就是归一化植被指数图，调一下色，把值高的调成绿色看看。和现实的比起来，还挺像的。

为了确保ndvi数值准确（在-1到1之间），这里我进入统计-空间统计检查一下数值。

点ok后到这个页面，勾选histograms更清楚。然后再点ok。

结果都分布在-1到1之间，效果拔群。

然后查看下面的表格，记录累计百分比在5%和95%的ndvi数值。（注意，这里图里选的是0.5%是错误的，感谢评论区提醒）
我这边分别是-0.488609和0.450226。
这个值用于计算植被覆盖度。

1.3.3植被覆盖度计算

根据这个公式：
*(b1 lt NDVI MIN)0+(b1 gt NDVI MAX)1+(b1 ge NDVI MIN and b1 le NDVI MAX)((b1-NDVI MIN)/ (NDVI MAX-NDVI MIN))
把上面ndvi获取的两个数值（-0.488609和0.450226）分别替换NDVI MIN和NDVI MAX，得到：
*(b1 lt -0.488609)0+(b1 gt 0.450226)1+(b1 ge -0.488609 and b1 le 0.450226)((b1+0.488609)/ (0.450226+0.488609))
打开波段计算，把上面这个公式放进去，点OK。

b1选之前得到的ndvi，输出命名为fvc，点OK。

得到这样一幅植被覆盖度图，并且给他调了下色（下图），这幅图的值在0-1之间，可以参考上面统计ndvi的方法检查一下。

1.3.4地表比辐射率计算

获取了植被覆盖度图，就可以计算地表比辐射率了。
利用这个公式：0.004*b1+0.986
点开波段计算，b1选择植被覆盖度图fvc，输出我命名为emiss。

然后得到这样一幅图（下图），这个就是地表比辐射率了。

2.地表温度反演

前面这个数据预处理好像有点长……不管了……
接下来没剩几步了，直接开始地表温度反演操作。
需要用到以下公式

有两个中间参数，分别是C和D。

2.1中间参数C

首先我们计算中间参数C，C是由地表比辐射率和大气透射率相乘得来。
地表比辐射率就是我们做的emiss图层，那大气透射率在哪呢？
在NASA官网中输入成影时间以及中心经纬度，则会提供大气透射率。（不过这个网站需要邮箱，QQ邮箱好像不能用……）
我根据我主图（行列号118029）查询大气透射率。
输入以下必要的信息。

等一会会返回这样一个网页，我们只需要记录下τ=0.73。

好了，计算中间参数C的参数我们都获取到了。
导入地表比辐射率，点开波段计算，公式为：b1*0.73，计算中间参数C。
输出命名为C。

得到下面这样一幅图。

2.2中间参数D

D的计算也是只要用到地表比辐射率和大气透射率。
同样是用波段计算，公式为：*(1-0.73)*(1+(1-b1)0.73)
输出命名为D。

处理完成后，得到下面这样一幅图。

2.3地表温度LST

万事俱备只欠东风。
根据上面的LST公式，我们可以写出这样的公式：
(-62.735657(1-b1-b2)+b3(0.434036*(1-b1-b2)+b1+b2)+b2Ta)/b1**
（当温度介于0-70°C 时，a=-62.735657，b=0.434036，这个是经验系数，不同的人可能经验系数不一样，这块我没具体地了解）
上式b1是C，b2是D，b3是辐射亮温，Ta是大气平均作用温度。
前三个温度我们都获取了，就差一个大气平均作用温度。
这里我用这个公式，Ta是大气平均作用温度，T0是近地面温度。

通过查询历史天气作为近地面温度计算得大气平均作用温度。因为我主图是凌晨2：21拍摄的，所以接近当日的最低温，查询得当日最低温为5摄氏度，在这里把它变换成开氏温度为278.15。所以T0=278.15，然后Ta=273.63353。
所有参数都确定了，补全我们的公式：（当然最后别忘减去273.15转化为摄氏度）
((-62.735657(1-b1-b2)+b3(0.434036(1-b1-b2)+b1+b2)-b2*273.63353)/b1)-273.15**
然后丢进波段运算进行计算，命名为LST_ra

然后得到这样一幅图，调个色什么的。温度大概在30度左右，说实话有点偏高了，经验系数ab那块可以再参考一下。

如果要生成可以导入ArcGIS的格式的话，可以在输出的末尾加一个.tif。
下图是导入ArcGIS里的效果。如果要使用的话，还需要掩膜一下。

大概就这样了。

-----------------------分割线（以下是乞讨内容）-----------------------