地温重建方法笔记

后来看到一篇博客写的很好,分享一下: https://www.cnblogs.com/enviidl/p/16285384.html

构建地温计算方程

单通道 (single-channel)方法是一种被广泛使用的用于从遥感数据重建地表温度的方法,其核心是根据热辐射收支平衡原构建的辐射传输方程和则计算地表热辐射强度.

对于遥感星上在λ波长观测的辐射值Lλsensor构建基本辐射传输方程为:

Lλsensor = [εB(λ, Ts) + (1 − ε)Lλatm↓τλ + Lλatm↑]

其中Lλatm↓是大气中向下穿过的辐射值,Lλatm↑是大气向上穿过的辐射值,ε为地表发射的辐射值,Ts为待求地表温度值.而B(λ, Ts)是温度为Ts的黑体辐射值,根据普朗克辐射定律当λ单位为 μm时:

$$ B(\lambda,T_s) = \frac{c_1 \lambda^{-5} }{\exp(\frac{c_2}{\lambda T_s}) - 1} $$

其中定标常数c1 = 1.19104 × 108(w ⋅ μm4 ⋅ m − 2 ⋅ sr − 1), c2 = 1.43877 × 104(μm ⋅ K).

由于辐射值与地表温度被揭示线性相关,因此可对其使用泰勒展开并忽略高阶项简化计算:

$$ B(\lambda,T_s) = B(\lambda,T_{sensor}) + \frac{\partial B(\lambda,T_{sensor})}{\partial T} (T_s - T_{sensor}) $$

其中偏微分结果由下式给出:

$$ \frac{\partial B(\lambda,T_{sensor})}{\partial T} = \frac{\int (\partial B(\lambda,T_{sensor} / \partial T) f(\lambda) d \lambda)}{\int f(\lambda)d \lambda} $$

其中Tsensor为基于传感器结果得到的温度数据,f(λ)表示表示给定波段的频谱响应函数.由 (Wang, et al. 2016)提供的模拟结果可以近似地认为:

$$ \frac{\partial B(\lambda,T_{sensor})}{\partial T} = a T_{sensor} + b $$

其中a和b分别为传感器参数,对于LandSat8的热传感器a = 0.001190,b =  − 0.21298,确定系数R2 = 1.

联立得到 Ts = {(1 − ετ)Lsensor − [1 + (1 − ε)ψ]φ}(ετγ) − 1 + Tsensor

其中: $$ \begin{aligned} \psi&= \frac{L_{\lambda}^{atm \downarrow} \tau(\theta)}{L_{\lambda}^{atm \uparrow} } \\ \varphi&=L_{\lambda}^{atm \uparrow} \\ \gamma&=a T_{sensor} +b \end{aligned} $$

整理得到,三个辐射参数 (地表发射的辐射值ε,传感器辐射值Lsensor和传感器获得温度Tsensor)和三个大气参数 (τ,ψ, andφ)是计算地表温度Ts所需的参数.

获取计算参数

获取大气参数

大气参数受到众多大气成分的影响,其中大气水分含量 (water vapor)的变化最大且影响对热辐射吸收能力最强,而其他成分则相对稳定. 因此可对大气参数关于水汽含量建模:

$$ \begin{pmatrix} \tau \\ \psi \\ \varphi_b \\ \Delta \varphi_l \\ \Delta \varphi_t \end{pmatrix} = C \cdot \begin{pmatrix} w^2 \\ w \\ 1 \end{pmatrix} = \begin{pmatrix} c_{11} & c_{12} & c_{13} \\ c_{21} & c_{22} & c_{23} \\ c_{31} & c_{32} & c_{33} \\ c_{41} & c_{42} & c_{43} \\ c_{51} & c_{52} & c_{53} \\ \end{pmatrix} \cdot \begin{pmatrix} w^2 \\ w \\ 1 \end{pmatrix} $$

本文使用 (Wang, et al. 2016)基于MODTRAN 4模拟结果得到的Landsat8卫星传感器校正系数:

$$ C_{L8B10} = \begin{pmatrix} -0.0027 & -0.0978 & -0.9949 \\ 0.0404 & -0.4839 & -2.0436 \\ -0.0389 & 1.2263 & -0.4706 \\ 0.1709 & -0.9764 & 0.5466 \\ 0.0219 & -0.1080 & 0.0741 \end{pmatrix} $$

获取地表辐射发射值

许多地表辐射发射值计算方法被发展,其中以TES算法,NDVI阈值法和基于地表真值的方法被广泛应用.本文使用基于NDVI阈值的方法计算地表辐射发射值.

NDVI阈值法使用两个NDVI阈值NDVIs和NDVIv将土地分为裸地 (NDVI < NDVIs), 植被茂盛 (NDVIv < NDVI)和混合 (NDVIs < NDVI < NDVIv)三种类型.从而计算植被覆盖度为:

$$ P_v=(\frac{NDVI - NDVI_s}{NDVI_v - NDVI_s})^2 $$

并构建地表辐射发射值ελ计算方程:

$$ \varepsilon_\lambda = \begin{cases} a_\lambda + b \lambda \rho_{red} &, NDVI < NDVI_s \\ \varepsilon_{v \lambda} P_v + \varepsilon_{s \lambda} (1 - P_v) + C_\lambda &, NDVI_s \leq NDVI < NDVI_v \\ \varepsilon_{v \lambda} + C_\lambda &, NDVI_v \leq NDVI \end{cases} $$

其中ρred为红色波段的反射率值,系数a和b为查阅光谱数据得到.εvλ是植被的辐射发射值,εsλ是裸地的辐射发射值,Cλ一个考虑到表面粗糙度引起的空腔效应的项,由下式给出:

Cλ = (1 − εsλεvλ) ⋅ F′ ⋅ (1 − Pv)

其中F′是一个0到1之间的几何参数,为了简化计算,我们根据经验结果对植被茂密类型的区域取C = 0.005, 对于混合类型的区域按土地利用类型设置不同的F′值:

获取遥感辐射亮度值

对遥感影像进行辐射校正后即可获得星上辐亮度,并根据公式进行亮温换算:

$$ T = \frac{K_2}{\ln (\frac{1}{K_1 / L} + 1)} $$

其中K1,K2为辐射校正常数,查阅元数据信息可知对于Landsat8数据产品K1 = 774.8853,K2 = 1321.0789.