基于Modis数据的北京市地表温度反演
评论区有下载原文和相关资料的链接，自己翻找即可。
操作平台
ENVI 5.5
ArcGIS 10.2
数据源
MODIS B1产品（包含1km 热红外波段）
数据来源
https://ladsweb.modaps.eosdis.nasa.gov/search/
研究区：北京市
研究时间：2019年9月
原理介绍
算法：劈窗算法
主要根据覃志豪研究成果进行，针对于陆地
T s = A 0 + A 1 T 31 − A 2 T 32 T_s=A_0+A_1T_{31}-A_2T_{32} Ts=A0+A1T31−A2T32
其中 T s {\ T}_s Ts 为地表温度，
A 0 、 A 1 、 A 2 A_0、A_1、A_2 A0、A1、A2 为参数，
T 31 、 T 32 T_{31}、T_{32} T31、T32 分别为 B 31 、 B 32 B31、B32 B31、B32 的亮度温度。
A 0 = a 31 D 32 ( 1 − C 31 − D 31 ) D 32 C 31 − D 31 C 32 − a 32 D 31 （ 1 − C 32 − D 32 ） D 32 C 31 − D 31 C 32 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 + D 31 D 32 C 31 − D 31 C 32 + b 31 D 32 ( 1 − C 31 − D 31 ) D 32 C 31 − D 31 C 32 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+D32C31−D31C32D31+D32C31−D31C32b31D32(1−C31−D31)
A 2 = D 31 D 32 C 31 − D 31 C 32 + b 32 D 31 ( 1 − C 31 − D 32 ) D 32 C 31 − D 31 C 32 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=D32C31−D31C32D31+D32C31−D31C32b32D31(1−C31−D32)
其中：
a 31 = − 64.60363 ， b 31 = 0.440817 a_{31}=-64.60363，b_{31}=0.440817 a31=−64.60363，b31=0.440817
a 32 = − 68.72575 ， b 32 = 0.47345 a_{32}=-68.72575，b_{32}=0.47345 a32=−68.72575，b32=0.47345
C i C_i Ci D i D_i Di 均为参数，运算如下
C i = ε i τ i C_i=\varepsilon_i\tau_i Ci=εiτi
ε i \varepsilon_i εi 为地表比辐射率， τ i \tau_i τi 为 i \ i i 热通道大气透过率
D i = ( 1 − τ i ) [ 1 + ( 1 − ε i τ i ) ] D_i=\left(1-\tau_i\right)[1+(1-\varepsilon_i\tau_i)] Di=(1−τi)[1+(1−εiτi)]
其中：
地表比辐射率计算
ε i = P v R v ε i v + ( 1 − P v ) R s ε i s + d ε i \varepsilon_i=P_vR_v\varepsilon_{iv}+(1-Pv)Rsεis+dεi εi=PvRvεiv+(1−Pv)Rsεis+dεi
P v P_v Pv 为地表植被覆盖度，可以通过归一化植被指数计算
R v R_v Rv 、 R s R_s Rs 分别为植被和裸土的辐射比率
其中：
ε 31 v = 0.98672 \varepsilon_{31v}=0.98672 ε31v=0.98672
ε 32 v = 0.98990 \varepsilon_{32v}=0.98990 ε32v=0.98990
ε 31 s = 0.96767 \varepsilon_{31s}=0.96767 ε31s=0.96767
ε 32 s = 0.97790 \varepsilon_{32s}=0.97790 ε32s=0.97790
R v = 0.92762 + 0.07033 P v R_v=0.92762+0.07033P_v Rv=0.92762+0.07033Pv
R s = 0.99782 + 0.08362 P v R_s=0.99782+0.08362P_v Rs=0.99782+0.08362Pv
P v = N D V I − N D V I S N D V I V − N D V I S P_v=\frac{NDVI-{NDVI}_S}{{NDVI}_V-{NDVI}_S} Pv=NDVIV−NDVISNDVI−NDVIS
其中：
N D V I V = 0.7 {NDVI}_V=0.7 NDVIV=0.7
N D V I S = 0.05 {NDVI}_S=0.05 NDVIS=0.05
d ε = 0.003796 min [ P v , ( 1 − P v ) ] \ d_\varepsilon=0.003796\min[P_{v\ },(1-P_v)] dε=0.003796min[Pv ,(1−Pv)]
N D V I = B 2 − B 1 B 2 + B 1 NDVI=\frac{B_2-B_1}{B_2+B_1} NDVI=B2+B1B2−B1
其中：
B i B_i Bi 为第 i 波段的反射率
大气透过率计算
W = ( α − ln r e f 19 r e f 2 ) / β W=(\alpha-\ln{\frac{{ref}_{19}}{{ref}_2}})/\beta W=(α−lnref2ref19)/β
其中：
W 为 水汽含量， r e f i {ref}_i refi 为 I 波段反射率
α = 0.02 \alpha=0.02 α=0.02
β = 0.651 \beta=0.651 β=0.651
y 31 = 5.72 − 4.69 e ( x / 42.05 ) y_{31}=5.72-4.69e^{(x/42.05)} y31=5.72−4.69e(x/42.05)
y 32 = − 1.25 + 2.28 e ( − x / 12.44 ) y_{32}=\ -1.25+2.28e^{(-x/12.44)} y32= −1.25+2.28e(−x/12.44)
亮温计算:
T 31 = K 31 , 2 ln ( 1 + K 31 , 1 R a d 31 ) T_{31}=\frac{K_{31,2}}{\ln{(1+\frac{K_{31,1}}{{Rad}_{31}})}} T31=ln(1+Rad31K31,1)K31,2
T 32 = K 32 , 2 ln ( 1 + K 32 , 1 R a d 32 ) T_{32}=\frac{K_{32,2}}{\ln{(1+\frac{K_{32,1}}{{Rad}_{32}})}} T32=ln(1+Rad32K32,1)K32,2
其中：
T i T_i Ti 为亮温
K 31 , 1 = 729.541636 K_{31,1}=729.541636 K31,1=729.541636
K 31 , 2 = 1304.413871 K_{31,2}=1304.413871 K31,2=1304.413871
K 32 , 1 = 474.684780 K_{32,1}=474.684780 K32,1=474.684780
K 32 , 2 = 1196.978785 K_{32,2}=1196.978785 K32,2=1196.978785
R a d i {Rad}_i Radi 为I 通道 辐射亮度
操作过程
技术流程图
具体操作
一、预处理（几何校正与辐射定标）
方法一：MCTK
选用ENVI提供的扩展工具MCTK，进行几何校正，几何校正后的结果同时也进行了定标。首先安装MCTK,然后即可在Toolbox中extensions中找到安装的扩展工具
打开工具，按照提示输入参数
方法二：手动定标
通过toolbox中的
工具查看，热红外数据集的scales和 offsets,并通过公式：
R a d i a n c e = s c a l e s × ( D N − o f f s e t s ) Radiance=scales\times(DN-offsets) Radiance=scales×(DN−offsets)
计算。
使用ENVI中的band math计算出结果
由于后续还需要用到NDVI，所以还需要对B1、B2进行定标。操作相同，不再赘述。
结果：MCTTK这种定标方式和手动定标结果有一定出入，所以暂时选择MCTK定标方式。
根据相关研究，做地表温度反演可以不用进行大气校正，所以就不进行大气校正了。
二、计算
- 计算亮温
计算T31、T32亮温
- 地表比辐射率计算
- NDVI计算
-
P v P_v Pv 计算
P v = b 1 g t 0.7 ∗ 1 + b 1 l t 0.05 ∗ 0 + b 1 g e 0.05 a n d b 1 l e 0.7 ∗ ( ( b 1 − 0.05 ) / ( 0.7 − 0.05 ) ) 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))
-
d ε d_\varepsilon dε 计算
d ε = ( b 1 e q 0 o r b 1 e q 1 ) ∗ 0 + ( b 1 g e 0 a n d b 1 l e 0.5 ) ∗ 0.003796 ∗ b 1 + ( b 1 g e 0.5 a n d b 1 l e 1 ) ∗ 0.0037968 ∗ 1 − b 1 + b 1 e q 0.5 ∗ 0.00189 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
- ε i \varepsilon_i εi 计算
b 1 ∗ ( 0.92762 + 0.07033 ∗ b 1 ) ∗ 0.98672 + ( 1 − b 1 ) ∗ ( 0.99782 + 0.08362 ∗ b 1 ) ∗ 0.96767 + b 2 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
- 大气透过率计算
-
W 水汽含量计算
-
T 31 − T 32 T31-T32 T31−T32 大气透过率计算
τ 31 = 5.72 − 4.69 e ( x / 42.05 ) \tau_{31}=5.72-4.69e^{(x/42.05)} τ31=5.72−4.69e(x/42.05)
-
τ 32 = − 1.25 − 2.28 e ( − x / 12.44 ) \tau_{32}=-1.25-2.28e^{(-x/12.44)} τ32=−1.25−2.28e(−x/12.44)
-
计算参数 C i D i C_i D_i CiDi
- C i C_i Ci 计算
C 31 = ε 31 × τ 31 C_{31}=\varepsilon_{31}\times\tau_{31} C31=ε31×τ31
C 32 = ε 32 × τ 32 C_{32}=\varepsilon_{32}\times\tau_{32} C32=ε32×τ32
-
D i D_i Di 计算
D i = 1 − τ i 1 + 1 − ε i τ i = ( 1 − τ i ) ( 2 − C i ) D_i=1-τi1+1-εiτi=(1-τi)(2-Ci) Di=1−τi1+1−εiτi=(1−τi)(2−Ci)D 31 = ( 1 − τ 31 ) ( 2 − C 31 ) D_{31}=(1-\tau_{31})(2-C_{31}) D31=(1−τ31)(2−C31)
D 32 = ( 1 − τ 32 ) ( 2 − C 32 ) D_{32}=(1-\tau_{32})(2-C_{32}) D32=(1−τ32)(2−C32)
-
计算参数 A 0 A 1 A 2 A_{0\ }\ \ A_1\ \ A_2 A0 A1 A2
- 计算 A 0 A_0 A0
A 0 = a 31 D 32 ( 1 − C 31 − D 31 ) D 32 C 31 − D 31 C 32 − a 32 D 31 （ 1 − C 32 − D 32 ） D 32 C 31 − D 31 C 32 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 A_1 A1
A 1 = 1 + D 31 D 32 C 31 − D 31 C 32 + b 31 D 32 ( 1 − C 31 − D 31 ) D 32 C 31 − D 31 C 32 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+D32C31−D31C32D31+D32C31−D31C32b31D32(1−C31−D31)
- 计算 A 2 A_2 A2
A 2 = D 31 D 32 C 31 − D 31 C 32 + b 32 D 31 ( 1 − C 31 − D 32 ) D 32 C 31 − D 31 C 32 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=D32C31−D31C32D31+D32C31−D31C32b32D31(1−C31−D32)
- 计算地表温度
T s = A 0 + A 1 T 31 − A 2 T 32 T_s=A_0+A_1T_{31}-A_2T_{32} Ts=A0+A1T31−A2T32
反演结果
- 由于中间过程步骤有一步单位需要换算，所以计算结果还需要 ÷ 10 \div10 ÷10
- ArcGIS 制图
在ARCGIS 中对反演结果进行可视化表达，加上等温线，显得更加直观。根据反演结果，发现北京市西南部温度比北部温度高，与现实情况，北京北部为怀柔、密云区，多山体和水域（密云水库）等，而北京市建成区在南部，所以呈现出这样的温度特征。说明反演结果较为准确。
![[MODIS数据处理#1]利用MRT工具预处理MODIS数据——以MOD16、MOD13为例](https://img-blog.csdnimg.cn/20190820215036284.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzM3OTQ4ODY2,size_16,color_FFFFFF,t_70)
