ENVI5.3.1基于Landsat 8影像进行辐射定标和大气校正
–
一、为什么要进行辐射定标和大气校正？
（1）辐射定标：辐射定标是将遥感传感器记录的无量纲DN值转换为具有实际物理意义的表现辐射亮度或者表观反射率的过程。
目的:消除传感器本身的误差，确定传感器入口处的准确辐射值。
（2）大气校正：将辐射亮度或者表面反射率转换为地表实际反射率。
目的：消除大气散射、吸收、反射引起的误差。
二、详细步骤
1. 数据获取
地理空间数据云（https://www.gscloud.cn/#page1/1）
获取了Landsat-8卫星三个时间三个地区的数据：
例如某个数据集文件名为：LC08_L1TP_130037_20210212_20210212_01_RT，LC08表示卫星名称为Landsat-8。L1TP（Level 1 Precision Terrain (Corrected)）表示L1级数据，TP表示该数据已经进行地形校正和几何校正。130037表示获取的数据地区编号在WRS-2参考系统中的条带号为130，行编号为37。第一个20210212表示数据获取时间为2021年2月12日，第二个20210212表示为数据处理时间为2021年2月12日。
2、数据预处理
2.1 辐射定标
2.1.1 多光谱波段
打开Landsat 8头文件，显示未处理的真彩色图像，如图：
图1 原始图像
打开Radiometric Correction中的Radiometric Calibration工具，选择多光谱数据。定标类型选择辐射亮度值（Radiance），输出格式为BIL，因为后续FLAASH大气校正的输入数据类型为BIL。输出数据类型为Float，系数为0.1（或者直接点击Apply FLAASH Settings）。确定输出路径后点击确定。这里输出数据名为OIL1_Radiometric_result.dat。如图：
定标结果验证：显示辐射定标结果图像，选择Display>Profiles>Spectral查看波谱曲线，看到定标后的数值主要集中在0-10范围内，单位是µW/(cm2 * sr * nm)如图。
2.1.2 热红外波段（获得Band10辐射亮度图像）
打开Radiometric Correction中的Radiometric Calibration工具。在File Selection对话框中，选择数据
Lco8_L1Tr_130037_20210212_20210212_01_KT_MTL_Thermal，单击Spectral Subset选择Thermal Infrared1（10.9），即第10波段，打开Radiometric Calibration面板，点击OK，获得辐射亮度图像，如图：
landsat8在这一步处理完后数值一般在5-11之间，温度越高的地方可能就越高。所以由于热岛效应，城市的值会更高。选择夏季月份的图像比较好。
2.1.3 辐射亮度温度计算（仅单窗算法需要）
公式：
K1、K2是卫星发射前的预设常量，具体值如下图所示：
Ti是辐射亮度温度（辐射亮温），Lλ是上面获取的辐射亮度（就是值大概在5-11之间的那幅图）。
输入公式：1321.08/alog(774.89/b1+1)
b1：Band10辐射亮度图像
2.2 大气校正
快速大气校正工具（简称QUAC）自动从图像上收集不同物质的波谱信息，获取经验值完成高光谱和多光谱的快速大气校正，它得到结果的精度近似FLAASH或者其他基于辐射传输模型的+/-15%。
打开Radiometric Correction——Atmospheric Correction Module——Quick Atmospheric Correction工具，导入经过辐射定标的数据OIL1_Radiometric_result.dat,点击OK。
再选好输出位置和名字。如图：
2.3 图像镶嵌
打开已经过辐射定标和大气校正的三幅影像，如图：
打开工具箱中的Mosaicking——Seamless Mosaic，单击加号按钮，把需要镶嵌的两幅图像选中。如果影像背景有黑色，可以在Data Ignore value中设置为0从而忽略黑色背景。
在Color Correction选项卡中选择Histogram Matching选项，Overlap Area Only是仅统计重叠区直方图进行匹配，Entire Scene是统计整幅图像直方图进行匹配。
随后在Seamlines中选择自动生成接边线Auto Generate Seamlines。
回到Main选项，在Feathering Distance列表中设置羽化距离，单位是像素。如果边缘或者接边线附近图像青光现象较为明显，建议设置大一点。同时在Seamlines/Feathering——Feathering选项卡内可以选择是边缘羽化Edge Feathering还是接边线羽化Seamline Feathering。这里选择接边线羽化，羽化距离500个像素。
最后在Export选项卡中导出镶嵌结果，选好输出位置和文件名，这里命名为OIL_mosaic_result.dat，重采样方法选择双三次卷积Cubic Convolution。点击finish。
因为本鼠鼠做的是温度反演，所以需要分别对OIL波段和TIRS波段进行图像镶嵌，但是之前没有对TIRS波段进行大气校正，因为只有2个波段吧，无法校正，如图：
2.4 图像裁剪
2.4.1 矢量数据下载
推荐全国地理信息资源目录服务系统
( https://www.webmap.cn/main.do?method=index)
点击“数据下载”，进入如下界面。
点击1：100万数据库（根据自己的需求选择），进入如下界面，这里有数据库的详细介绍。
然后选择你要下载的区域，想要哪个区域，就点一下，将其加入成果车。将所有你想下载的区域全部加入成果车以后就可以去结算了。点击成果车进去。按照提示将该填的信息填好，然后提交订单即可。提交完成后，点击全部数据，便出现了下载，然后点击下载，就ok啦！有的数据可能会等几天的时间才能申请下来，有的数据当场就可以下载。
2.4.2 ArcGIS扣出自己需要的shp文件
这里推荐这个up主的视频，里面有详细的介绍
ArcGIS快速抠出shp文件局部区域
2.4.3 详细步骤
将若尔盖行政区shp文件导入ENVI，并加载到视图中。
选择Regions of Interest——Subset Data from ROIs工具，在弹出的窗口中选择卫星影像数据。
点击OK后设置参数，Select Input ROIs选择若尔盖行政区shp文件，是否掩膜多边形外的像元Mask pixels outside of ROI?选择是，同时掩膜裁剪背景值设置为0，设置输出路径名后点击OK即可。
裁剪结果如下：
3.温度反演
3.1基于单窗算法
3.1.1 公式
其中 ε 表示地表比辐射率，t为当天大气的透射率；
a、b为常量，a=-67.355351，b=0.458606；
C、D为中间变量；
Ta为大气平均作用温度（单位：K）；
Tb可以用普朗克公式的反函数获取。
大气平均作用温度 Ta与近地面气温 T0（一般为2m）存在以下的线性关系：（注：Ta、T0的单位均为华氏度 K，而不是摄氏度，且影像用的是格林尼治时间，历史天气检索时用北京时间）
格林尼治时间+8=北京时间：（http://www.timebie.com/cn/greenwichmeanbeijing.php）
其中，Ta是大气平均作用温度，T0 是遥感图像获取时当地的气温。
Tb是辐射亮度温度（详见2.1.3）
3.1.2 计算
3.1.2.1NDVI计算
在Toolbox中输入NDVI，选择裁切好的多光谱OIL波段数据，点击OK。
打开NDVI calculator parameters对话框，在NDBI band 栏的red和near IR分别输入4、5（OLI的第4、5波段对应红和近红外波段），设置好输出位置和名称，计算得到NDVI。
3.1.2.2植被覆盖度数据计算
Pv为植被覆盖度，通常使用公式Pv = [(NDVI- NDVISoil)／(NDVIVeg - NDVISoil)]来计算，其中，NDVI为归一化植被指数，NDVISoil为完全是裸土或无植被覆盖区域的NDVI值，NDVIVeg则代表完全被植被所覆盖的像元的NDVI值，即纯植被像元的NDVI值。取经验值NDVIVeg = 0.70和NDVISoil = 0.05，即当某个像元的NDVI大于0.70时，Pv取值为1；当NDVI小于0.05，Pv取值为0。
在Toolbox里面选择Band math输入公式：
(b1 gt 0.7)*1+(b1 lt 0.05)0+(b1 ge 0.05 and b1 le 0.7)((b1-0.05)/(0.7-0.05))
b1即为NDVI，计算之后可以得到植被覆盖率图像
3.1.2.3地表比辐射率计算
通用比辐射率计算公式 0.004*b1+0.986，但是一般不采用
更准确的计算可利用ndvi将影像分为三部分，水体，自然表面，城镇，分别计算比辐射率做法为：
将ndvi值小于等于0的比辐射率赋值为0.995，大于0小于0.7的比辐射率赋值为(0.9589+0.086b2-0.0671b2b2)，大于等于0.7的比辐射率赋值为(0.9625+0.0614b2-0.0461*b2^2)
公式为：
(b1 le 0)0.995+(b1 gt 0 and b1 lt 0.7)(0.9589+0.086b2-0.0671b2^2)+(b1 ge 0.7)(0.9625+0.0614b2-0.0461*b2^2)
b1：ndvi
b2 植被覆盖度
计算之后可以得到地表比辐射率数据
3.1.2.3 大气透射率
在NASA官网(http://atmcorr.gsfc.nasa.gov/)中输入成影时间以及中心经纬度，则会提供上式中所需要的参数。得到：
大气在热红外波段的透过率τ为0.67，
大气向上辐射亮度L↑为2.97 W/(m2·sr·μm)，
大气向下辐射亮辐射亮度L↓为4.69W/(m2·sr·μm)。
3.14 地表温度计算
第一步，计算C: b1b2
b1：地表比辐射率； b2：大气透射率
第二步，计算D: (1-b1)*(1+(1-b1)b2)
b1：地表比辐射率； b2：大气透射率
第三步，计算LST：
(-67.355351(1-b1-b2)+b3(0.458606*(1-b1-b2)+b1+b2)-b2*11.1377265)/b1
【当温度介于0-70°C 时，a=-67.3554，b=0.4586】
b1:指的是C; b2：指的是D;
b3：指的是辐射亮温Tb；b4：指的是大气平均作用温度Ta
关于大气平均作用温度计算的Ta：
选择当地气温值。依据：我们平时生活中所了解到的气温，是在离地面1.5米高的百叶箱中的温度表上测得，具有较好的代表性。
3.2 基于大气校正法
3.2.1 原理
实现原理：卫星传感器接收到的热红外辐射亮度值Lλ由三部分组成：大气向上辐射亮度L↑，地面的真实辐射亮度经过大气层之后到达卫星传感器的能量；大气向下辐射到达地面后反射的能量。卫星传感器接收到的热红外辐射亮度值Lλ的表达式可写为（辐射传输方程）：
Lλ = [εB(TS) + (1-ε)L↓] τ + L↑
式中，ε为地表比辐射率，TS为地表真实温度(K)，B(TS)为黑体热辐射亮度，τ为大气在热红外波段的透过率。则温度为T的黑体在热红外波段的辐射亮度B(TS)为：
B(TS) = [Lλ - L↑- τ(1-ε)L↓]/τ ε
Ts可以用普朗克公式的函数获取。
TS = K2/ln(K1/ B(TS)+ 1)
对于TM，K1 =607.76 W/(m2µmsr)，K2 =1260.56K。
对于ETM+，K1=666.09 W/(m2µmsr)，K2 =1282.71K。
对于TIRS Band10，K1= 774.89 W/(m2µmsr)，K2 = 1321.08K。
3.2.2 计算NDVI
跟上面操作一样
3.2.3 计算植被覆盖度数据
跟上面操作一样
3.2.4 计算地表比辐射率
跟上面操作一样
3.2.5 辐射亮度计算
查询大气剖面数据(http://atmcorr.gsfc.nasa.gov/)，输入相关参数可得到：
大气在热红外波段的透过率τ：0.60
大气向上辐射亮度L↑：3.22 W/(m2·sr·μm)
大气向下辐射亮辐射亮度L↓：5.03W/(m2·sr·μm)
公式：B(TS) = [Lλ - L↑- τ(1-ε)L↓]/τ ε
在bandmath中输入表达式：(b2-3.22-0.60*(1-b1)5.03)/(0.60*b1)
其中b2指的是最开始求的第10波段的辐射亮度值
b1指的是刚求的地表比辐射率
计算之后可以得到黑体辐射亮度影像
3.2.6 地表温度计算
公式：TS = K2/ln(K1/ B(TS)+ 1)
对于TM，K1 =607.76 W/(m2µmsr)，K2 =1260.56K。
对于ETM+，K1=666.09 W/(m2µmsr)，K2 =1282.71K。
对于TIRS Band10，K1= 774.89 W/(m2µmsr)，K2 = 1321.08K。
在bandmath中输入表达式：(1321.08)/alog(774.89/b1 + 1) - 273
（盗用一下别人的图图），本鼠鼠暂时还打开不了大气校正参数计算的网站。
其中b1指的是黑体辐射亮度
计算之后即可得到地表温度数据了
4 结果浏览与输出
在图层管理器（Layer Manager）中的地表温度图像图层，右键选择 Raster Color Slices。将温度划分为五个区间：
30℃以上；25℃至30℃；20℃至25℃；15℃至20℃；低于15℃。
分别浏览几个温度区间的空间分布范围。
得出统计反演结果。
参考文章：
记 / ENVI5.3温度反演-单窗算法-操作要点（小白尔）文章来源：https://www.toymoban.com/news/detail-438706.html
ENVI5.3土地覆盖分类、地表温度反演详细操作流程记录**
使用大气校正法对landsat-8tirs地表温度进行反演**
ENVI5.3.1使用Landsat 8影像进行辐射定标和大气校正实例操**作
ENVI5.3.1使用Landsat 8影像进行图像镶嵌和图像裁剪实例操作
ENVI5.3.1使用Landsat 8影像进行NDVI计算实例操作
ENVI53 辐射校正、大气校正、影像裁剪超详细教程
全国地理信息资源服务系统行政边界矢量数据下载教程文章来源地址https://www.toymoban.com/news/detail-438706.html
到了这里，关于ENVI5.3.1Landsat 8影像基于单窗算法和辐射传输方程进行地表温度反演的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
