基于Modis数据的地表温度反演

基于Modis数据的北京市地表温度反演

图片压缩有点看不清 ，留言可发word原文。

操作平台

ENVI 5.5
ArcGIS 10.2

数据源

MODIS B1产品（包含1km 热红外波段）
数据来源
https://ladsweb.modaps.eosdis.nasa.gov/search/
研究区：北京市
研究时间：2019年9月

原理介绍

算法：劈窗算法
主要根据覃志豪研究成果进行，针对于陆地
$T_s=A_0+A_1T_{31}-A_2T_{32}$Ts​=A0​+A1​T31​−A2​T32​

其中 ${\ T}_s$ Ts​ 为地表温度，
$A_0、A_1、A_2$A0​、A1​、A2​ 为参数，
$T_{31}、T_{32}$T31​、T32​ 分别为 $B31、B32$B31、B32 的亮度温度。

$A_0=\frac{a31D32(1-C31-D31)}{D32C31-D31C32}-\frac{a32D31（1-C32-D32）}{D32C31-D31C32}$A0​=D32C31−D31C32a31D32(1−C31−D31)​−D32C31−D31C32a32D31（1−C32−D32）​

$A_1=1+\frac{D_{31}}{D_{32}C_{31}-D_{31}C_{32}}+\frac{b_{31}D_{32}(1-C_{31}-D_{31})}{D_{32}C_{31}-D_{31}C_{32}}$A1​=1+D32​C31​−D31​C32​D31​​+D32​C31​−D31​C32​b31​D32​(1−C31​−D31​)​

$A_2=\frac{D_{31}}{D_{32}C_{31}-D_{31}C_{32}}+\frac{{b_{32}D}_{31}(1-C_{31}-D_{32})}{D_{32}C_{31}-D_{31}C_{32}}$A2​=D32​C31​−D31​C32​D31​​+D32​C31​−D31​C32​b32​D31​(1−C31​−D32​)​
其中：
$a_{31}=-64.60363，b_{31}=0.440817$a31​=−64.60363，b31​=0.440817
$a_{32}=-68.72575，b_{32}=0.47345$a32​=−68.72575，b32​=0.47345

$C_i$Ci​ $D_i$Di​ 均为参数，运算如下

$C_i=\varepsilon_i\tau_i$Ci​=εi​τi​

$\varepsilon_i$εi​ 为地表比辐射率，$\tau_i$τi​ 为 $\ i$ i 热通道大气透过率

$D_i=\left(1-\tau_i\right)[1+(1-\varepsilon_i\tau_i)]$Di​=(1−τi​)[1+(1−εi​τi​)]

其中：

地表比辐射率计算

$\varepsilon_i=P_vR_v\varepsilon_{iv}+(1-Pv)Rsεis+dεi$εi​=Pv​Rv​εiv​+(1−Pv)Rsεis+dεi

$P_v$Pv​ 为地表植被覆盖度，可以通过归一化植被指数计算
$R_v$Rv​ 、$R_s$Rs​ 分别为植被和裸土的辐射比率

其中：
$\varepsilon_{31v}=0.98672$ε31v​=0.98672
$\varepsilon_{32v}=0.98990$ε32v​=0.98990
$\varepsilon_{31s}=0.96767$ε31s​=0.96767
$\varepsilon_{32s}=0.97790$ε32s​=0.97790
$R_v=0.92762+0.07033P_v$Rv​=0.92762+0.07033Pv​
$R_s=0.99782+0.08362P_v$Rs​=0.99782+0.08362Pv​
$P_v=\frac{NDVI-{NDVI}_S}{{NDVI}_V-{NDVI}_S}$Pv​=NDVIV​−NDVIS​NDVI−NDVIS​​
其中：
${NDVI}_V=0.7$NDVIV​=0.7
${NDVI}_S=0.05$NDVIS​=0.05
$\ d_\varepsilon=0.003796\min[P_{v\ },(1-P_v)]$ dε​=0.003796min[Pv ​,(1−Pv​)]

$NDVI=\frac{B_2-B_1}{B_2+B_1}$NDVI=B2​+B1​B2​−B1​​
其中：
$B_i$Bi​ 为第 i 波段的反射率

大气透过率计算

$W=(\alpha-\ln{\frac{{ref}_{19}}{{ref}_2}})/\beta$W=(α−lnref2​ref19​​)/β
其中：
W 为 水汽含量，${ref}_i$refi​ 为 I 波段反射率
$\alpha=0.02$α=0.02
$\beta=0.651$β=0.651

$y_{31}=5.72-4.69e^{(x/42.05)}$y31​=5.72−4.69e(x/42.05)
$y_{32}=\ -1.25+2.28e^{(-x/12.44)}$y32​= −1.25+2.28e(−x/12.44)

亮温计算:

$T_{31}=\frac{K_{31,2}}{\ln{(1+\frac{K_{31,1}}{{Rad}_{31}})}}$T31​=ln(1+Rad31​K31,1​​)K31,2​​

$T_{32}=\frac{K_{32,2}}{\ln{(1+\frac{K_{32,1}}{{Rad}_{32}})}}$T32​=ln(1+Rad32​K32,1​​)K32,2​​

其中：

$T_i$Ti​ 为亮温
$K_{31,1}=729.541636$K31,1​=729.541636
$K_{31,2}=1304.413871$K31,2​=1304.413871
$K_{32,1}=474.684780$K32,1​=474.684780
$K_{32,2}=1196.978785$K32,2​=1196.978785

${Rad}_i$Radi​ 为I 通道 辐射亮度

操作过程

技术流程图

具体操作

一、预处理（几何校正与辐射定标）

方法一：MCTK

选用ENVI提供的扩展工具MCTK，进行几何校正，几何校正后的结果同时也进行了定标。首先安装MCTK,然后即可在Toolbox中extensions中找到安装的扩展工具

打开工具，按照提示输入参数

方法二：手动定标

通过toolbox中的工具查看，热红外数据集的scales和 offsets,并通过公式：
$Radiance=scales\times(DN-offsets)$Radiance=scales×(DN−offsets)
计算。
使用ENVI中的band math计算出结果
由于后续还需要用到NDVI，所以还需要对B1、B2进行定标。操作相同，不再赘述。

结果：MCTTK这种定标方式和手动定标结果有一定出入，所以暂时选择MCTK定标方式。
根据相关研究，做地表温度反演可以不用进行大气校正，所以就不进行大气校正了。

二、计算

1. 计算亮温
   计算T31、T32亮温
   
2. 地表比辐射率计算
   1. NDVI计算

2. $P_v$Pv​ 计算
   $P_v=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))$Pv​=b1gt0.7∗1+b1lt0.05∗0+b1ge0.05andb1le0.7∗((b1−0.05)/(0.7−0.05))
   

3. $d_\varepsilon$dε​ 计算
   $d_\varepsilon=\left(b_1\ eq\ 0\ or\ b_1\ eq\ 1\right)\ast0+\left(b_1\ ge\ 0\ and\ b_1\ le\ 0.5\right)\ast0.003796\ast b_1+\left(b_1\ ge\ 0.5\ and\ b_1\ le\ 1\right)\ast0.0037968\ast1-b1+b1eq 0.5*0.00189$dε​=(b1​ eq 0 or b1​ eq 1)∗0+(b1​ ge 0 and b1​ le 0.5)∗0.003796∗b1​+(b1​ ge 0.5 and b1​ le 1)∗0.0037968∗1−b1+b1eq0.5∗0.00189

4. $\varepsilon_i$εi​ 计算
   $b1*(0.92762+0.07033*b1)*0.98672+(1-b1)*(0.99782+0.08362*b1)*0.96767+b2$b1∗(0.92762+0.07033∗b1)∗0.98672+(1−b1)∗(0.99782+0.08362∗b1)∗0.96767+b2

3. 大气透过率计算

   1. W 水汽含量计算
      

   2. $T31-T32$T31−T32 大气透过率计算

      $\tau_{31}=5.72-4.69e^{(x/42.05)}$τ31​=5.72−4.69e(x/42.05)

$\tau_{32}=-1.25-2.28e^{(-x/12.44)}$τ32​=−1.25−2.28e(−x/12.44)

4. 计算参数 $C_i D_i$Ci​Di​

   1. $C_i$Ci​ 计算

   $C_{31}=\varepsilon_{31}\times\tau_{31}$C31​=ε31​×τ31​

$C_{32}=\varepsilon_{32}\times\tau_{32}$C32​=ε32​×τ32​

2. $D_i$Di​ 计算
   $D_i=1-τi1+1-εiτi=(1-τi)(2-Ci)$Di​=1−τi1+1−εiτi=(1−τi)(2−Ci)

   $D_{31}=(1-\tau_{31})(2-C_{31})$D31​=(1−τ31​)(2−C31​)
   
   $D_{32}=(1-\tau_{32})(2-C_{32})$D32​=(1−τ32​)(2−C32​)

4. 计算参数 $A_{0\ }\ \ A_1\ \ A_2$A0 ​  A1​  A2​

   1. 计算 $A_0$A0​

   $A_0=\frac{a31D32(1-C31-D31)}{D32C31-D31C32}-\frac{a32D31（1-C32-D32）}{D32C31-D31C32}$A0​=D32C31−D31C32a31D32(1−C31−D31)​−D32C31−D31C32a32D31（1−C32−D32）​
   

2. 计算 $A_1$A1​

$A_1=1+\frac{D_{31}}{D_{32}C_{31}-D_{31}C_{32}}+\frac{b_{31}D_{32}(1-C_{31}-D_{31})}{D_{32}C_{31}-D_{31}C_{32}}$A1​=1+D32​C31​−D31​C32​D31​​+D32​C31​−D31​C32​b31​D32​(1−C31​−D31​)​

3. 计算 $A_2$A2​

$A_2=\frac{D_{31}}{D_{32}C_{31}-D_{31}C_{32}}+\frac{{b_{32}D}_{31}(1-C_{31}-D_{32})}{D_{32}C_{31}-D_{31}C_{32}}$A2​=D32​C31​−D31​C32​D31​​+D32​C31​−D31​C32​b32​D31​(1−C31​−D32​)​

6. 计算地表温度
   $T_s=A_0+A_1T_{31}-A_2T_{32}$Ts​=A0​+A1​T31​−A2​T32​
   

反演结果

1. 由于中间过程步骤有一步单位需要换算，所以计算结果还需要 $\div10$÷10
2. ArcGIS 制图
   在ARCGIS 中对反演结果进行可视化表达，加上等温线，显得更加直观。根据反演结果，发现北京市西南部温度比北部温度高，与现实情况，北京北部为怀柔、密云区，多山体和水域（密云水库）等，而北京市建成区在南部，所以呈现出这样的温度特征。说明反演结果较为准确。