利用热红外遥感影像提取地表温度的草垛
一、简介及数据来源
地表温度作为地球环境分析的重要指标，而遥感技术作为现代重要的观测手段，使得基于遥感图像的地表温度反演的研究越来越多。主要的地表温度反演方法有：大气校正法，单窗算法，单通道法等等。本文以2013年花溪区的遥感影像为基础，介绍用辐射传输方程法对地表温度进行反演。
以Landsat8 LOI_TIRS卫星数字产品，下载花溪大致在图像中的遥感影像为基础数据，用“花溪区行政边界裁剪”得到预处理的图像。
二、技术流程
第一步：数据预处理
热红外数据使用的是Landsat8遥感卫星中的第十波段，已经做了传感器定标、几何校正、工程区裁剪，详细流程如下：
1.在主界面中，选择【File】打开元数据“LC81270422013272LGN00_MTL.txt”文件，ENVI会自动按照波长分为5个数据集：多光谱数据（1-7波段）、全色波段（8波段）、卷云波段数据（9波段）、热红外数据（10、11波段）和质量数据（12波段）。
注意事项：热红外数据有10和11两个波段，我们用第10波段来做辐射亮度值，用来后期的计算。
2.在Toolbox工具箱中，选择【Radiometric Correction】工具，运用其中的【Radiometric Calibration】。在【File Selection】对话框中，选择数据“LC81270422013272LGN00_MTL_Thermal”，单击“Spectral Subest”选择“Thermal Infrared1”，打开Radiometric Calibration面板。如图1所示：
图 1
3.在Radiometric Calibration面板中，设置一下参数：
l定标类型（Calibration Type）：辐射亮度值（radiance）。
l其他选择默认的参数。如图2所示：
图 2
4.选择输出路径和文件名称，单击OK按钮，执行定标处理。得到Band10辐射亮度图像。
5.工程区裁剪
（1）打开矢量文件“花溪区行政边界.shp”，并显示在视窗中，选择Linear2%拉伸显示。
（2）在Toolbox工具箱中，双击【Regions of Interest】/【Subest from ROIs】；在弹出的Select Input File Subset via ROIs面板中选择“fushedingbiao10.dat”，单击OK按钮。
（3）在Spatial Subset via ROI Parameters面板中，设置一下参数：
lSelect Iuput ROIs（在ROI列表中）：选择矢量文件“花溪区行政边界”。
lMask pixels outside of ROI？（是否掩膜多边形外的像元）：选择Yes。
lMask Background Value（裁剪背景值）：0。
（4）选择输出路径及文件名，单击OK按钮，裁剪图像。
如图3所示：
图 3
第二步：地表比辐射率的计算
物体的比辐射率是物体向外辐射电磁波的能力表征。它不仅依赖于地表物体的组成，而且与物体的表面状态(表面粗糙度等)及物理性质(介电常数、含水量等)有关，并随着所测定的波长和观测角度等因素有关。在大尺度上对比辐射率精确测量的难度很大，目前只是基于某些假设获得比辐射率的相对值，本文主要根据可见光和近红外光谱信息来估计比辐射率。
(一)植被覆盖度计算
1.计算植被覆盖度Pv采用的是混合像元分解法，将整景影像的地类大致分为水体、植被和建筑，具体的计算公式如下：
ε=0.004Pv+0.986
Pv= (NDVI- NDVIS)／(NDVIV - NDVIS)
式中，NDVI为归一化植被指数；NDVIS为完全被裸土或无植被覆盖区域的NDVI值，NDVIV则代表完全被植被所覆盖的像元的NDVI值，即纯植被像元的NDVI值。取经验值NDVIV=0.70和NDVIS=0.05，即当某个像元的NDVI大于0.70时，Pv取值为1；当NDVI小于0.05，Pv取值为0。
（1）在Toolbox工具箱中，双击【Spetral】/【Vegetation】中的NDVI工具，在文件输入对话框中选择Landsat8OLI大气校正结果。如图4所示：
图 4
注：覃志豪等（2004）提出使用原始DN值图像计算NDVI工具对反演结果影响不大。
（2）在NDVI Calculation parameters对话框中，选择NDVI计算波段：Red：4；Near IR ：5。
注：Landsat8 OLI数据Red是4，NIR为5，如果换为TM数据则Red为3，NIR为4。
（3）选择输出文件名和路径。如图5所示;
图 5
（4）在Toolbox中，选择【Band Ratio】/【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图像，计算得到植被覆盖度图像。如图6所示：
图 6
（5）在Toolbox中，选择【Band Ratio】/【Band Math】，输入表达式：
0.004*b1+0.986
其中b1为植被覆盖度图像，计算得到地表比辐射率图像。如图7所示：
图 7
注：为了得到更精准的地表比辐射率数据，可以使用覃志豪等（2004）提出的先将地表分成水体、自然表面和城镇区，分别针对3种地表类型计算地表比辐射率：
l水体像元比辐射率：0.995
l自然表面像元比辐射率：εsurface = 0.9625 + 0.0614PV - 0.0461PV2、εbuilding = 0.9589 + 0.086PV - 0.0671PV2
式中，εsurface和εbuilding分别代表自然表面像元和城镇像元的比辐射率。
第三步：计算相同温度下黑体的辐射亮度值
辐射传输方程法，又称大气校正法，其基本思路为：首先利用与卫星过空时间同步的大气数据来估计大气对地表热辐射的影响。然后把这部分大气影响从卫星高度上传感器所观测到的热辐射总量中减去。从而得到地表热辐射强度．再把这一热辐射强度转化为相应的地表温度.
卫星传感器接收到的热红外辐射亮度值Lλ由三部分组成：大气向上辐射亮度L↑，地面的真实辐射亮度经过大气层之后到达卫星传感器的能量；大气向下辐射到达地面后反射的能量。卫星传感器接收到的热红外辐射亮度值的表达式可写为（辐射传输方程）：
Lλ = [ε·B(TS) + (1-ε)L↓]·τ + L↑
这里，ε为地表辐射率，TS为地表真实温度，B(TS)为普朗克定律推到得到的黑体在TS的热辐射亮度，τ为大气在热红外波段的透过率。则温度为T的黑体在热红外波段的辐射亮度B(TS)为：
B(TS) = [Lλ - L↑- τ·(1-ε)L↓]/τ·ε
在NASA公布的网站查询，此时需要通过https://atmcorr.gsfc.nasa.gov/网站获取大气剖面参数，输入的信息可以在数据云网站上查找。
注：输入的中心纬度；成像时间等信息不能从地理空间数据云中找到该遥感影像获取，则需从下载解压之后的元数据中的txt文件中获取得到，最后得到的大气剖面参数为（如图8）：
图 8
l大气在热红外波段的通过率为：0.89
l大气向上辐射亮度：0.92 (m2·sr·μm)
l大气向下辐射亮度：1.55 (m2·sr·μm)
提示：由于缺少地表相关参数（气压、温度和相对湿度等信息），得到的结果是基于模型计算的。
（1）在Toolbox中，选择【Band Ratio】/【Band Math】，输入表达式：
(b2-0.92-0.89*(1-b1)*1.55)/(0.89*b1)
式中，b1为地表比辐射率图像；
b2为Band10辐射亮度图像。计算得到相同温度下的黑体辐射亮度图像。如图9所示：
图 9
第四步：反演地表温度
在获取温度为TS的黑体在热红外波段的辐射亮度后，根据普朗克公式的反函数，求得地表真实温度TS：
TS = K2/ln(K1/ B(TS)+ 1)
对于OLI数据：K1 =774.89W/(m2·sr·μm)，K2 =1321.08K
对于ETM+，K1 =666.09W/(m2·sr·μm)，K2 =1282.71K
对于TM数据，K1 =607.76W/(m2·sr·μm)，K2 =1260.56K
注：此数据的来源是根据元图像的txt文件中获取。
（1）在Toolbox中，选择【Band Ratio】/【Band Math】，输入表达式：
(1321.08)/alog(774.89/b1+1)-273
式中：b1为相同温度下的黑体辐射亮度值图像，得到地表温度图像。如图10所示：
图 10
第五步:结果浏览与输出
在Display中显示温度值，是一个灰度的单波段图像。
（1）在图层管理器（Layer Manager）中的地表温度图像图层上，右键选择“Raster Color Slices”，将温度划分为五个区间：<17℃、17-21℃、21-26℃、26-31℃和>31℃
五个等级。如图11所示：
（2）分别浏览各个温度区间的空间分布范围。
（3）统计反演结果得出80.7%区域的温度集中在21-26℃之间。
注意事项：
1.在用元文件数据计算Band10的辐射亮度值的时候，计算得到的数据是没有裁剪过的，如果在用Band Math进行计算的时候，会发现无法利用以及辐射定标获取的Band10辐射亮度值带入公式中进行计算，此时需要把图像按照“花溪区行政边界”进行裁剪，方可进行后续操作。
2.在运用Band Math时，输入公式，需要用计算器能识别的语言才能进行计算，一定不能马虎大意，哪怕只错一个小数点的位置，都会导致公式的错误而不能进行下一步的操作，这里一般检查的就是输入进去的公式是否有误。
3.注意得到的每一个结果，保存的路径要用英文名称，并且路径不能太长。否则会导致错误发生，而不能进行下一步的操作。
4.最后得到的温度值，输出的数据单位即摄氏度，大致是在我们生活的温度区间范围内，如果出现几千、几百或几万的数据，就是错误的。查看结果的值在【Layer Manager】中，单击图层右键，打开Quick Stats即可查看范围。
5.缺少同步数据温度测量数据用于验证反演结果，查询2013年9月29日贵阳市花溪区的气温，根据做出来的图像对比，才能证明图像有一定的参考价值。