基于ETM数据地表温度反演实验单通道算法操作文档

基于ETM 数据煤田火区地表温度反演的研究实验操作步骤与流程
算法：单通道算法，其公式为
Τs =γ[ε−1(ψ1L sensor +ψ2)+ψ3]+δ （1）
γ={c 2L sensor T sensor 2[λ4
c 1L sensor +λ−1]}−1
（2） δ=−γL sensor +T sensor （3）
L sensor =L min (λ)+[L max (λ)−L min (λ)]Q DN Q max ⁄ （4）
T sensor =K 2ln (1+K 1L (λ)⁄)
⁄ （5） K 1=666.09(mW ∙cm −2∙sr −1∙um −1), K 2=1282.71K
ψ1=0.1471∙ω2−0.1558ω+1.1234 （6）
ψ2=−1.1836∙ω2−0.3761ω−0.5289 （7）
ψ3=−0.0455∙ω2+1.8719ω−0.3907 （8）
ω=0.177e +0.339 （9）
e =0.6108∗exp [17.27(Τ0−273)
237.3+Τ0−273]∗RH （10） 先来说明单通道算法公式（1）中γ，L sensor ，δ，ψ1，ψ2，ψ3等这些参数的计算过程，地表比辐射率ε的计算过程稍后在说明。

（1）对于ψ1，ψ2，ψ3的计算，只要查阅资料得知相对湿度RH ，与温度Τ0后，就可以算出大气中水蒸汽的含量ω，进而可以根据公式算出ψ1，ψ2，ψ3。

（2）对于L sensor 的计算，也就是辐射校正的过程，主要目的在于把影像中像元的灰度值转化成辐亮度L sensor ，公式（4）中的L min (λ)，L max (λ)，Q max 在影像头文件中可以找到，Q DN 就是所要进行校正的影像。

在ENVI 中的操作如下： Basic tools → band math ,然后点开出现如下左侧对话框：
对于ETM 数据热红外波段高增益就是L sensor =3.2+9.45∙Q DN 255⁄ 然后点ok 出现如下右侧对话框：
选择b1为需要校正的波段
指定输出文件夹点ok即可。

（3）辐射校正完了进行大气校正。

（特别注意：对于TM/ETM数据，大气校正的波段不能包含第六波段，大气校正之前需要把辐射校正完后的数据的BSQ格式转成BIL/BIP格式，完后还需要对影像各波段的波长中心值wavelength进行编辑。

）如果是单波段数据需要先进行波段叠加（layer stacking）
具体操作如下：
Basic tools→layer stacking,全部选中所有的波段，指定输出文件夹点ok即可。

完了进行格式转换，如下：
Basic tools→convert data (BSQ,BIL,BIP),选择波段叠加后的影像，
点击ok后出现：
指定输出文件夹点击ok即可。

对输入数据进行头文件编辑，主要是对波长wavelenth（即每一波段的波长中心值）和波长宽度fwhm（每一波段的波长范围）的编辑。

不是高光谱数据可以不对fwhm进行编辑。

（envi——file——Edit Envi Header）
完后可以进行大气校正（使用的是FLAASH）
envi→spectral→flaash,出现如下对话框：
该对话框分三部分，上面主要为输入输出文件夹的设定，中间部分包含影像中心坐标，传感器类型，卫星飞行时间，下面部分主要是大气模型与气溶胶模型的反演，对于多光谱数据可以不做光谱打磨（spectral polishing）和高级设置（advanced settings）具体不在此说明，各项参数设置完后点击apply即可。

（4）进行T sensor的计算，Basic tools→band math输入公式后如下左侧图，指定b1,如下右图，b1就是辐射定标后的波段L sensor
最后指定文件夹输出即可。

（5）计算γ，公式（2）中c1，c2，λ都是常数，再利用波段运算就可得到: Basic tools→band math,如下图左侧，指定b1与b2的波段如下右图：
其中b1为L sensor波段，b2为T sensor波段
指定输出文件夹就好。

（6）计算δ，公式（3）中涉及到的已经全部为已知量，再利用波段运算就好。

Basic tools→band math，如下图左侧，指定b1,b2和b3的波段如下右图:
其中b1为L sensor波段，b2为T sensor波段
B3为γ波段，指定输出文件夹就可以得到
δ波段。

下面主要介绍地表比辐射率ε的计算过程：
在文中我们主要通过归一化植被指数阈值法(NDVI THM)来确定比辐射率ε首先计算一副影像的NDVI，计算公式如下：
NDVI=ρ4−ρ3ρ4+ρ3
（1）在此之前需要对该景ETM影像的3,4波段进行表观反射率的计算（注意：不能用辐射校正后的数据，应该使用原始数据来计算表观反射率）具体操作如下：Basic Tools→Preprocessing→Calibration Utilities→Landsat Calibration后出现如下对话框，分别选择3波段与4波段，
点击ok之后，出现下面对话框：
传感器类型，飞行时间及太阳高度角都可以在头文件中找到。

如上红色椭圆处都设置好了后，点击Edit Calibration Parameters按钮，出现如下对话框：
点击ok，指定输出文件夹就完成了对3波段的表观反射率计算，4波段的操作同上。

（2）完后，运用波段运算进行NDVI的计算，如下左图，指定b1,b2波段后如下右图（b1为4波段，b2为3波段），
点击ok就可得到NDVI的影像如下：
在影像上点击右键，再点击Quick statistics,就有如下统计图：
（3）根据上图及NDVI所占的百分比，可以确定出NDVImin，NDVImax，在根据如下公式计算植被覆盖度Ρν：
Ρν=[
NDVI−NDVI min NDVI max−NDVI min
]
2
对于本文而言，NDVImin=−0.01，NDVImax=0.4，再用波段运算计算Ρν，如下左图所示，指定b1波段（b1为ndvi波段）如下右图，输出可得到Ρν波段。

（4）计算完Ρν后，可根据如下公式计算地表比辐射率ε：
ε=0.9625+0.061Ρv−0.0461Ρv2
利用波段运算计算地表比辐射率ε，如下左图，指定B1（B1就是上面的Ρν波段）后如下右图所示，输出可得到ε波段：
到此，单通道算法公式中所有的参数
都已经确定。

下面利用波段运算，按照单通道算法的公式，进行地表温度的反演：
本文中ω为0.4877，ψ1,ψ2和ψ3分别为1.0824，-0.9938，0.5114.波段运算如下图所示，
指定b1,b2,b3,b4,各个波段（b1为γ波段，b2为ε波段，b3为L sensor波段，b4为δ波段）。

如下图所示
最后的地表温度反演结果图如下：
最后进行密度分割，对于不同的温度区间着色，具体操作如下：
在主影像窗口点击overlay→density slice,选择结果影像图，清除默认区间，重新设置就好，并着色，如下图：
到此，地表温度反演的所有操作已经完成。