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

回到主页

操作平台

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} \]

其中 \({\ T}_s\) 为地表温度，
\(A_0、A_1、A_2\) 为参数，
\(T_{31}、T_{32}\) 分别为 \(B31、B32\) 的亮度温度。

\[ A_0=\frac{a31D32(1-C31-D31)}{D32C31-D31C32}-\frac{a32D31（1-C32-D32）}{D32C31-D31C32} \]
\[ 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}} \]
\[ 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}} \]

其中：
\(a_{31}=-64.60363，b_{31}=0.440817\)
\(a_{32}=-68.72575，b_{32}=0.47345\)

\(C_i\) \(D_i\) 均为参数，运算如下

\[ C_i=\varepsilon_i\tau_i \]

\(\varepsilon_i\) 为地表比辐射率，\(\tau_i\) 为 \(\ i\) 热通道大气透过率

\[ D_i=\left(1-\tau_i\right)[1+(1-\varepsilon_i\tau_i)] \]

其中：

地表比辐射率计算

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

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

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

\[ NDVI=\frac{B_2-B_1}{B_2+B_1} \]

其中：
\(B_i\) 为第 i 波段的反射率

大气透过率计算

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

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

亮温计算:

\[ T_{31}=\frac{K_{31,2}}{\ln{(1+\frac{K_{31,1}}{{Rad}_{31}})}} \]
\[ T_{32}=\frac{K_{32,2}}{\ln{(1+\frac{K_{32,1}}{{Rad}_{32}})}} \]

其中：

\(T_i\) 为亮温
\(K_{31,1}=729.541636\)
\(K_{31,2}=1304.413871\)
\(K_{32,1}=474.684780\)
\(K_{32,2}=1196.978785\)

\({Rad}_i\) 为I 通道 辐射亮度

操作过程

技术流程图

具体操作

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

方法一：MCTK

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

方法二：手动定标

通过toolbox中的工具查看，热红外数据集的scales和 offsets,并通过公式：
\(Radiance=scales\times(DN-offsets)\)
计算。
使用ENVI中的band math计算出结果 由于后续还需要用到NDVI，所以还需要对B1、B2进行定标。操作相同，不再赘述。
结果：MCTTK这种定标方式和手动定标结果有一定出入，所以暂时选择MCTK定标方式。 根据相关研究，做地表温度反演可以不用进行大气校正，所以就不进行大气校正了。

二、计算

  1. 计算亮温 计算T31、T32亮温

  2. 地表比辐射率计算

    1. NDVI计算

    1. \(P_v\) 计算

    \[ 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)) \]

    1. \(d_\varepsilon\) 计算

    \[ 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 \]

    1. \(\varepsilon_i\) 计算

    \[ b1*(0.92762+0.07033*b1)*0.98672+(1-b1)*(0.99782+0.08362*b1)*0.96767+b2 \]

  3. 大气透过率计算

    1. W 水汽含量计算

    2. \(T31-T32\) 大气透过率计算

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

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

  4. 计算参数 \(C_i D_i\)

    1. \(C_i\) 计算

    \[ C_{31}=\varepsilon_{31}\times\tau_{31} \]

    \[ C_{32}=\varepsilon_{32}\times\tau_{32} \]

    1. \(D_i\) 计算

    \[ D_i=1-τi1+1-εiτi=(1-τi)(2-Ci) \]
    \[ D_{31}=(1-\tau_{31})(2-C_{31}) \]

    \[ D_{32}=(1-\tau_{32})(2-C_{32}) \]

  5. 计算参数 \(A_{0\ }\ \ A_1\ \ A_2\)

    1. 计算 \(A_0\)

    \[ A_0=\frac{a31D32(1-C31-D31)}{D32C31-D31C32}-\frac{a32D31（1-C32-D32）}{D32C31-D31C32} \]

    1. 计算 \(A_1\)

    \[ 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}} \]

    1. 计算 \(A_2\)

    \[ 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}} \]

  6. 计算地表温度

\[ T_s=A_0+A_1T_{31}-A_2T_{32} \]

反演结果

  1. 由于中间过程步骤有一步单位需要换算，所以计算结果还需要 \(\div10\)

  2. ArcGIS 制图 在ARCGIS 中对反演结果进行可视化表达，加上等温线，显得更加直观。根据反演结果，发现北京市西南部温度比北部温度高，与现实情况，北京北部为怀柔、密云区，多山体和水域（密云水库）等，而北京市建成区在南部，所以呈现出这样的温度特征。说明反演结果较为准确。

本文章使用limfx的vsocde插件快速发布