基于大气校正法的Landsat8TIRS反演地表温度
分享
-
2016-06-09
热红外遥感（Infrared Remote Sensing）是指传感器工作波段限于红外波段范围之内的遥感。即利用星载或机载传感器收集、记录地物的热红外信息，并利用这种热红外信息来识别地物和反演地表参数如温度、湿度和热惯量等。目前有很多的卫星携带了热红外传感器，包括ASTER、AVHRR、MODIS、TM/ETM+/ TIRS等。 目前，地表温度反演算法主要有以下三种：大气校正法（也称为辐射传输方程：Radiative Transfer Equation——RTE）、单通道算法和分裂窗算法。
本实例是基于大气校正法，利用Landsat8 TIRS反演地表温度。
基本原理：首先估计大气对地表热辐射的影响, 然后把这部分大气影响从卫星传感器所观测到的热辐射总量中减去, 从而得到地表热辐射强度, 再把这一热辐射强度转化为相应的地表温度。
具体实现为：卫星传感器接收到的热红外辐射亮度值Lλ由三部分组成：大气向上辐射亮度L↑，地面的真实辐射亮度经过大气层之后到达卫星传感器的能量；大气向下辐射到达地面后反射的能量。卫星传感器接收到的热红外辐射亮度值Lλ的表达式可写为（辐射传输方程）：
Lλ = [εB(TS) + (1-ε)L↓]τ + L↑ （1.1）
式中，ε为地表比辐射率，TS为地表真实温度(K)，B(TS)为黑体热辐射亮度，τ为大气在热红外波段的透过率。则温度为T的黑体在热红外波段的辐射亮度B(TS)为：
B(TS) = [Lλ - L↑- τ(1-ε)L↓]/τε （1.2）
Ts可以用普朗克公式的函数获取。
TS = K2/ln(K1/ B(TS)+ 1) （1.3）
对于TM，K1 =607.76 W/(m2*µm*sr)，K2 =1260.56K。
对于ETM+，K1=666.09 W/(m2*µm*sr)，K2 =1282.71K。
对于TIRS Band10，K1= 774.89 W/(m2*µm*sr)，K2 = 1321.08K。
从上可知此类算法需要2个参数：大气剖面参数和地表比辐射率。大气剖面参数 在NASA提供的网站(http://atmcorr.gsfc.nasa.gov/)中，输入成影时间以及中心经纬度可以获取大气剖面参数。适用于只有一个热红外波段的数据，如Landsat TM /ETM+/TIRS数据。
主要内容就是使用BandMath工具计算公式（1.2）和公式（1.3），处理流程如下图所示。
注：OLI大气校正是可选项，根据前人研究成果，大气校正处理对结果影响不是很大。
图1.1 基于大气校正法的TIRS反演流程图
1、图像辐射定标
(1)在主界面中，选择File→Open，在文件选择对话框中选择"LC81230322013276LGN00_MTL.txt"文件，ENVI自动按照波长分为五个数据集：多光谱数据（1-7波段），全色波段数据（8波段），卷云波段数据（9波段），热红外数据（10，11波段）和质量波段数据（12波段）。
（2）在Toolbox工具箱中，选择Radiometric Correction/Radiometric Calibration。在File Selection对话框中，选择数据LC81230322013132LGN02_MTL_Thermal，单击Spectral Subset选择Thermal Infrared1（10.9），打开Radiometric Calibration面板。
（3）在Radiometric Calibration面板中，设置以下参数：
（4）选择输出路径和文件名1-LC81230322_band10_rad.dat，单击OK按钮执行定标处理。得到Band10辐射亮度图像。
2、 地表比辐射率计算
TIRS的Band10热红外波段与TM/ETM+6热红外波段具有近似的波谱范围，本例采用TM/ETM+6相同的地表比辐射率计算方法。使用Sobrino提出的NDVI阈值法计算地表比辐射率。
ε=0.004Pv+0.986 （1.4）
其中， Pv是植被覆盖度，用以下公式计算：
Pv = [(NDVI- NDVISoil)／(NDVIVeg - NDVISoil)] （1.5）
其中，NDVI为归一化植被指数，NDVISoil为完全是裸土或无植被覆盖区域的NDVI值，NDVIVeg则代表完全被植被所覆盖的像元的NDVI值，即纯植被像元的NDVI值。取经验值NDVIVeg = 0.70和NDVISoil = 0.05，即当某个像元的NDVI大于0.70时，Pv取值为1；当NDVI小于0.05，Pv取值为0。
注：这里采用简化的植被覆盖度计算模型，感兴趣的可以使用更加精确的植被覆盖度计算模型。
（1）在Toolbox工具箱中，双击Spectral/Vegetation/NDVI工具，在文件输入对话框中，选择Landsat8 OLI多光谱图像。
图1.2 NDVI文件输入对话框
提示：覃志豪提出使用原始的DN值图像计算NDVI对反演结果影响不大。
（2）在NDVI Calculaton parameters对话框中，自动识别NDVI计算波段：Red：4，Near IR：5。
图1.3 NDVI对话框
（3）选择输出文件名和路径。
（4）在Toobox中，选择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。
计算得到植被覆盖度图像。
（5）在Toobox中，选择Band Ratio/Band Math，输入表达式：
0.004*b1+0.986
其中，b1：植被覆盖度图像。
计算得到地表比辐射率图像。
提示：为了得到更精确的地表比辐射率数据，可使用覃志豪等提出的先将地表分成水体、自然表面和城镇区，分别针对三种地表类型计算地表比辐射率：
3、黑体辐射亮度与地表温度计算
在NASA公布的网站查询（http://atmcorr.gsfc.nasa.gov），输入成影时间：2013-10-03 02:55和中心经纬度（Lat：40.32899857，Lon：116.70610046），以及其他相应的参数，得到大气剖面信息为：
提示：由于缺少地表相关参数（气压、温度、相对湿度等信息），得到的结果是基于模型计算的结果。
（1）依据公式（1.2），在Toolbox工具箱中，双击Band Ratio/Band Math工具，输入表达式：
(b2-0.75-0.9*(1-b1)*1.29)/(0.9*b1)
其中，b1：地表比辐射率图像
b2：Band10辐射亮度图像
计算得到同温度下的黑体辐射亮度图像。
（2）依据公式（1.3），在Toobox中，双击Band Ratio/Band Math工具，输入表达式：
(1321.08)/alog(774.89/b1+1)-273
其中，b1：同温度下的黑体辐射亮度图像
得到单位为摄氏度的地表温度图像。
提示：公式（1.3中），TIRS Band10的K1和K2 是从*_MTL.txt元数据文件中获取。
（3）在图层管理器（Layer Manager）中的地表温度图像图层，右键选择 Raster Color Slices。将温度划分为四个区间：
提示：缺少同步温度测量数据用于验证反演结果，查询2013年10月3号北京市最低气温为10°，最高气温为22°。本示例反演结果大部分在这个区间内，反演结果有一定的参考价值。
文章来源：http://blog.sina.com.cn/s/blog_764b1e9d0102wa8s.html
本实例是基于大气校正法，利用Landsat8 TIRS反演地表温度。
基本原理：首先估计大气对地表热辐射的影响, 然后把这部分大气影响从卫星传感器所观测到的热辐射总量中减去, 从而得到地表热辐射强度, 再把这一热辐射强度转化为相应的地表温度。
具体实现为：卫星传感器接收到的热红外辐射亮度值Lλ由三部分组成：大气向上辐射亮度L↑，地面的真实辐射亮度经过大气层之后到达卫星传感器的能量；大气向下辐射到达地面后反射的能量。卫星传感器接收到的热红外辐射亮度值Lλ的表达式可写为（辐射传输方程）：
Lλ = [εB(TS) + (1-ε)L↓]τ + L↑ （1.1）
式中，ε为地表比辐射率，TS为地表真实温度(K)，B(TS)为黑体热辐射亮度，τ为大气在热红外波段的透过率。则温度为T的黑体在热红外波段的辐射亮度B(TS)为：
B(TS) = [Lλ - L↑- τ(1-ε)L↓]/τε （1.2）
Ts可以用普朗克公式的函数获取。
TS = K2/ln(K1/ B(TS)+ 1) （1.3）
对于TM，K1 =607.76 W/(m2*µm*sr)，K2 =1260.56K。
对于ETM+，K1=666.09 W/(m2*µm*sr)，K2 =1282.71K。
对于TIRS Band10，K1= 774.89 W/(m2*µm*sr)，K2 = 1321.08K。
从上可知此类算法需要2个参数：大气剖面参数和地表比辐射率。大气剖面参数 在NASA提供的网站(http://atmcorr.gsfc.nasa.gov/)中，输入成影时间以及中心经纬度可以获取大气剖面参数。适用于只有一个热红外波段的数据，如Landsat TM /ETM+/TIRS数据。
主要内容就是使用BandMath工具计算公式（1.2）和公式（1.3），处理流程如下图所示。
注：OLI大气校正是可选项，根据前人研究成果，大气校正处理对结果影响不是很大。
图1.1 基于大气校正法的TIRS反演流程图
1、图像辐射定标
(1)在主界面中，选择File→Open，在文件选择对话框中选择"LC81230322013276LGN00_MTL.txt"文件，ENVI自动按照波长分为五个数据集：多光谱数据（1-7波段），全色波段数据（8波段），卷云波段数据（9波段），热红外数据（10，11波段）和质量波段数据（12波段）。
（2）在Toolbox工具箱中，选择Radiometric Correction/Radiometric Calibration。在File Selection对话框中，选择数据LC81230322013132LGN02_MTL_Thermal，单击Spectral Subset选择Thermal Infrared1（10.9），打开Radiometric Calibration面板。
（3）在Radiometric Calibration面板中，设置以下参数：
- 定标类型（Calibration Type）：辐射亮度值（radiance）。
- 其他选择默认参数。
（4）选择输出路径和文件名1-LC81230322_band10_rad.dat，单击OK按钮执行定标处理。得到Band10辐射亮度图像。
2、 地表比辐射率计算
TIRS的Band10热红外波段与TM/ETM+6热红外波段具有近似的波谱范围，本例采用TM/ETM+6相同的地表比辐射率计算方法。使用Sobrino提出的NDVI阈值法计算地表比辐射率。
ε=0.004Pv+0.986 （1.4）
其中， Pv是植被覆盖度，用以下公式计算：
Pv = [(NDVI- NDVISoil)／(NDVIVeg - NDVISoil)] （1.5）
其中，NDVI为归一化植被指数，NDVISoil为完全是裸土或无植被覆盖区域的NDVI值，NDVIVeg则代表完全被植被所覆盖的像元的NDVI值，即纯植被像元的NDVI值。取经验值NDVIVeg = 0.70和NDVISoil = 0.05，即当某个像元的NDVI大于0.70时，Pv取值为1；当NDVI小于0.05，Pv取值为0。
注：这里采用简化的植被覆盖度计算模型，感兴趣的可以使用更加精确的植被覆盖度计算模型。
（1）在Toolbox工具箱中，双击Spectral/Vegetation/NDVI工具，在文件输入对话框中，选择Landsat8 OLI多光谱图像。
图1.2 NDVI文件输入对话框
提示：覃志豪提出使用原始的DN值图像计算NDVI对反演结果影响不大。
（2）在NDVI Calculaton parameters对话框中，自动识别NDVI计算波段：Red：4，Near IR：5。
图1.3 NDVI对话框
（3）选择输出文件名和路径。
（4）在Toobox中，选择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。
计算得到植被覆盖度图像。
（5）在Toobox中，选择Band Ratio/Band Math，输入表达式：
0.004*b1+0.986
其中，b1：植被覆盖度图像。
计算得到地表比辐射率图像。
提示：为了得到更精确的地表比辐射率数据，可使用覃志豪等提出的先将地表分成水体、自然表面和城镇区，分别针对三种地表类型计算地表比辐射率：
- 水体像元比辐射率：0.995
- 自然表面像元比辐射率：εsurface = 0.9625 + 0.0614Pv - 0.0461Pv2
- 城镇区像元比辐射率：εbuilding = 0.9589 + 0.086Pv - 0.0671Pv2
3、黑体辐射亮度与地表温度计算
在NASA公布的网站查询（http://atmcorr.gsfc.nasa.gov），输入成影时间：2013-10-03 02:55和中心经纬度（Lat：40.32899857，Lon：116.70610046），以及其他相应的参数，得到大气剖面信息为：
- 大气在热红外波段的透过率τ：0.90
- 大气向上辐射亮度L↑：0.75 W/(m2·sr·μm)
- 大气向下辐射亮辐射亮度L↓：1.29W/(m2·sr·μm)
提示：由于缺少地表相关参数（气压、温度、相对湿度等信息），得到的结果是基于模型计算的结果。
（1）依据公式（1.2），在Toolbox工具箱中，双击Band Ratio/Band Math工具，输入表达式：
(b2-0.75-0.9*(1-b1)*1.29)/(0.9*b1)
其中，b1：地表比辐射率图像
b2：Band10辐射亮度图像
计算得到同温度下的黑体辐射亮度图像。
（2）依据公式（1.3），在Toobox中，双击Band Ratio/Band Math工具，输入表达式：
(1321.08)/alog(774.89/b1+1)-273
其中，b1：同温度下的黑体辐射亮度图像
得到单位为摄氏度的地表温度图像。
提示：公式（1.3中），TIRS Band10的K1和K2 是从*_MTL.txt元数据文件中获取。
（3）在图层管理器（Layer Manager）中的地表温度图像图层，右键选择 Raster Color Slices。将温度划分为四个区间：
- 25℃以上
- 22℃至25℃
- 20℃至22℃
- 15℃至20℃
- 低于15℃ （4）分别浏览几个温度区间的空间分布范围。 （5）统计反演结果得出81%区域的温度集中在15~22°区间。
提示：缺少同步温度测量数据用于验证反演结果，查询2013年10月3号北京市最低气温为10°，最高气温为22°。本示例反演结果大部分在这个区间内，反演结果有一定的参考价值。
文章来源：http://blog.sina.com.cn/s/blog_764b1e9d0102wa8s.html
