CN107560734A - 一种Landsat8遥感反演地表温度方法 - Google Patents

Abstract

本发明公开了一种Landsat8遥感反演地表温度方法，计算大气水汽含量w、计算Landsat8两个热红外波段对应的大气透过率τ₁₀和τ₁₁、计算Landsat8两个热红外波段对应的地表比辐射率ε₁₀和ε₁₁、反演计算地表温度T_S：完成步骤一至步骤三的基础上，将以上参数计算结果代入反演公式计算地表温度值。本发明通过不同的算法思路，使Landsat8遥感反演地表温度算法技术得到改进和提升，使得利用Landsat8热红外波段反演地表温度精度得到提高；通过简单而直接的逻辑关系，使得外专业使用者能够较好地理解和掌握该算法技术，使得该算法技术能够得到广泛应用，使其产生一定的综合效益和综合价值。

Description

一种Landsat8遥感反演地表温度方法

技术领域

本发明属于遥感对地观测领域，尤其涉及一种Landsat8遥感反演地表温度 方法。

背景技术

随着遥感对地观测技术的发展，通过卫星搭载热红外波段可以实现对地表 温度的提取和监测。90年代以来，经过多方研究，遥感反演地表温度的算法得 到广泛研究，截至目前，其大致形成了三种类型：单通道、双通道和多角度反 演算法。Qin等(2001)提出了适用于Landsat TM热红外遥感数据的地表温度单 窗反演算法，随后，Jiménez-&Sobrino(2003)和Jiménez-等(2009) 提出了一种通用的单通道温度反演算法。Li等(2004)、Ottlé等(1992)、Price 等(1983)等学者也深入研究了遥感单通道反演地表温度算法。同时，双通道(劈 窗)反演算法也在遥感地表温度反演中具有广泛应用和研究。单通道地表温度 遥感反演算法是通过一个热红外遥感波段提取地表温度，比起双通道反演算法需要多估计一个变量——大气有效温度，因此会比双通道反演算法具有更大的 误差，由于双通道反演算法利用了两个热红外波段的辐射信息，通过解算热红 外波段大气传输方程组，消去大气有效温度变量，使得地表温度反演结果比单 通道算法更加准确，另外多角度反演算法和多通道反演算法思想类似，但是由 于多角度热红外遥感数据并不普遍，使其没有双通道反演算法应用广泛。 Landsat8卫星2013年发射升空，其搭载两个新的热红外波段band10和band11， 对于Landsat8两个热红外波段，需要开发相应的地表温度反演算法。Jimenez-Munoz等(2014)、Rozenstein等(2014)、Du等(2014)在2014年研发 了适用于Landsat8遥感数据的地表温度双通道——劈窗反演算法，这些算法在 本质上都是对大气传输方程中的普朗克热辐射函数进行了统一的简化处理，并 且都采用了泰勒线性展开低阶部分作为算法模型主要部分。

Rozenstein等(2014)提出的Landat8遥感反演地表温度劈窗算法技术与本 发明最为接近，其主要技术是基于大气传输方程，具体如下式所示：

B_i(T_i)＝ε_iτ_iB_i(T_s)+(1-τ_i)(1+(1-ε_i)τ_i)B_i(T_a)

该算法技术通过对上式进行泰勒线性展开，并取第一阶展开部分化解上述计算公式，得出如下计算公式：

T_s＝A₀+A₁T₁₀-A₂T₁₁

其中：

A₀＝E₁a₁₀+E₂a₁₁

A₁＝1+A+E₁b₁₀

A₂＝A+E₂b₁₁

C_i＝ε_iτ_i(θ)

D_i＝[1-τ_i(θ)][1+(1-ε_i)τ_i(θ)]

A＝D₁₀/E₀

E₁＝D₁₁(1-C₁₀-D₁₀)/E₀

E₂＝D₁₀(1-C₁₁-D₁₁)/E₀

E₀＝D₁₁C₁₀-D₁₀C₁₁

该算法技术中需要估计大气透过率、地表比辐射率等参数，其中，大气透 过率的估计也是通过大气水汽含量w与大气透过率的关系间接求得，该算法技 术认为两者呈线性关系。

通过大气传输软件Modtran4.0模拟数据评估的该算法技术的地表温度反演 精度的均方根误差(RMSE)为0.93℃，最大绝对误差为1.9742℃。

与本发明最为相近的算法技术的缺点主要有:

第一，该算法技术的反演精度稍低，原因在于其只取了大气传输方程泰勒 线性展开的低阶部分，其高阶部分影响没有得到充分考虑，另外，大气透过率 与大气水汽含量w的相关关系使用了简单的线性关系。

第二，该算法技术的逻辑结构比较复杂，对于外专业使用者而言较难理解， 在一定程度上影响到了其现实实用性。该缺点原因在于该算法技术涉及的间接 变量比较多，由此产生了较多的嵌套关系，使得外专业使用者理解具有一定困 难。

发明内容

本发明的目的在于提供一种实用的Landsat8遥感反演地表温度方法，旨在 解决现有类似算法技术的反演精度稍低和逻辑结构复杂的问题。通过对大气传 输方程中普朗克辐射函数采取不同简化方式，充分考虑了高阶部分，同时利用 三次多项式拟合大气透过率与大气水汽含量w的函数关系，在一定程度上提高 了Landsat8遥感地表温度反演精度。同时本发明的逻辑结构简约，非常有利于 外专业使用者理解和掌握。

本发明是这样实现的，其所述的Landsat8遥感反演地表温度方法包括以下 步骤：

步骤一：计算大气水汽含量w：利用同期MODIS遥感数据计算大气水汽含量 w，计算公式如下所示：

W＝((0.02-ln(b19/b2))/0.651)²

其中，b19和b2对应与MODIS数据中的第19波段和第2波段对应的辐射灰度 值；

步骤二：计算Landsat8两个热红外波段对应的大气透过率τ₁₀和τ₁₁：在步 骤一计算的基础上，将像元的大气水汽含量w代入如下公式计算大气透过率：

τ₁₀＝0.9570356-0.0277340w-0.0333734w²+0.0028800w³，

τ₁₁＝0.9456728-0.0857755w-0.0290912w²+0.0032169w³；

步骤三：计算Landsat8两个热红外波段对应的地表比辐射率ε₁₀和ε₁₁：首 先计算目标区域的值被归一化指数NDVI，然后按照以下公式计算Landsat8两个 热红外波段对应的地表比辐射率，

ε₁₀＝0.984p_v+0.964(1-p_v)+(1-0.964)(1-p_v)0.492，

ε₁₁＝0.980p_v+0.970(1-p_v)+(1-0.970)(1-p_v)0.490，

其中，p_v代表植被覆盖度，NDVI代表像元的植被归一化指数值，NDVI_min和 NDVI_max分别代表计算空间范围内植被归一化指数的最小值和最大值；

步骤四：反演计算地表温度T_S：完成步骤一至步骤三的基础上，将以上参数 计算结果代入以下公式计算地表温度值：

各参数具体计算如下：

A₁₀＝0.0006678ε₁₀τ₁₀，

A₁₁＝＝0.0006188ε₁₁τ₁₁，

B₁₀＝-0.2333226ε₁₀τ₁₀，

B₁₁＝-0.1990475ε₁₀τ1₀，

C₁₀＝0.1312942(1-τ₁₀)(1+(1-ε₁₀)τ₁₀)，

C₁₁＝0.1387986(1-τ₁₁)(1+(1-ε₁₁)τ₁₁)，

D₁₀＝21.1666266ε₁₀τ₁₀-26.7808503(1+(1-ε₁₀)τ₁₀)-B₁₀(T₁₀)，

D₁₁＝16.7224278ε₁₀τ₁₀-27.7043284(1+(1-ε₁₀)τ₁₀)-B₁₁(T₁₁)，

其中，B₁₀(T₁₀)和B₁₁(T₁₁)分别是Landsat8两个热红外波段对应亮温T₁₀和 T₁₁的普朗克热辐射值。

进一步，计算大气水汽含量w前，对遥感数据进行预处理、对MODIS影像 几何纠正、对Landsat8 OLI多光谱遥感数据进行大气校正。

本发明另一目的在于提供利用上述的Landsat8遥感反演地表温度方法的Landsat8遥感反演地表温度系统。

本发明通过不同的算法思路，使Landsat遥感反演地表温度技术得到改进 和提升，使得利用Landsat8热红外波段反演地表温度精度得到提高；通过简单 而直接的逻辑关系，使得外专业使用者能够较好地理解和掌握该技术，使得该 技术能够得到广泛应用，使其产生一定的综合效益和综合价值；本发明能够突 破传统遥感反演地表温度的劈窗算法技术的思路，为使用者提供一种实用而具 有较强的可操作性的Landsat8遥感反演地表温度技术。

本发明Landsat8卫星遥感数据具有信噪比高、现势性好和获取成本低等优 势，通过简单易用的地表温度反演技术，可以大大促进热环境改善、生态环境 演变等方面的研究水平，同时也可以实现Landsat8遥感监测地表温度业务化、 流程化和工具化，带来广泛的社会、经济和生态效益。

附图说明

图1是本发明实施例提供的Landsat8遥感反演地表温度方法流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例， 对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以 解释本发明，并不用于限定本发明。

以太原市城区为实例区域，选取覆盖实例区域的无云晴朗天气条件下夏季 2013年6月27日的Landsat8OLI多光谱和TIR热红外遥感数据以及对应同日 相近时相的MODIS数据做为地表温度反演数据源。遥感数据投影为通用横轴墨 卡托投影(UTM)北49区，投影基准为WGS-84椭球。

如图1所示，本发明实施例提供的Landsat8遥感反演地表温度方法，步骤 如下：

S101：遥感数据进行预处理。对MODIS影像几何纠正，对Landsat8OLI多 光谱遥感数据进行大气校正。

S102：大气透过率参数计算。预处理完成后，利用MODIS数据产品MOD021KM 反演大气水汽含量W₁₀、w₁₁，再根据大气透过率与水汽含量的三次拟合公式计 算大气透过率τ₁₀和τ₁₁。

S103：地表比辐射率计算。根据Landsat8OLI多光谱数据计算归一化值被 指数NDVI和植被覆盖度p_v，然后代入根据相关公式计算地表比辐射率ε₁₀和ε₁₁。

S104：计算Landsat8热红外波段星上亮温T₁₀和T₁₁,具体计算公式可查阅Landsat8官方技术文档，也可由遥感处理软件ENVI5.1直接计算而得。

S105：反演计算地表温度。在大气透过率、地表比辐射率以及星上亮温都 计算出以后，根据地表温度反演公式计算实例区域地表温度。

本发明是根据其对地表温度反演的误差影响不同，采取了不同的数学处理， 热红外波段亮温热辐射函数B_i(T_i)通过原始函数公式直接计算，地表温度热辐射 函数B_i(T_S)二次函数拟合处理，大气有效温度热辐射函数B_i(T_a)线性拟合处理， 然后通过解算方程组，即可得出地表温度T_S的计算公式，该过程逻辑简单，并没 有嵌套的数学结构，所用到的数学知识简单易懂，因此大大降低了使用门槛。

本发明Landsat8卫星遥感数据具有信噪比高、现势性好和获取成本低等优 势，通过简单易用的地表温度反演技术，可以大大促进热环境改善、生态环境 演变等方面的研究水平，同时也可以实现Landsat8遥感监测地表温度业务化、 流程化和工具化，带来广泛的社会、经济和生态效益。

本发明技术同时兼具有精度高、可操作性强和数学逻辑简单等特点，在逻 辑上满足奥卡姆剃刀原则，即“如无必要，勿增实体”的“简单有效原理”。为 了与现有相近技术比较，同样通过大气传输软件Modtran4.0模拟数据评估本发 明的地表温度反演精度，结果表明本发明反演的地表温度结果的均方根误差 (RMSE)为0.51℃，最大绝对误差为0.9933℃，均低于现有最接近的算法技术 精度，表1列出了具体的精度验证结果。除此之外，还利用实测数据对本发明 技术反演的地表温度进行了验证，验证结果表明均方根误差(RMSE)为0.70℃， 最大绝对误差为1.1℃，也均小于现有接近技术反演精度。

另外，本发明技术以遥感大气传输方程为基础，如下式所示：

B_i(T_i)＝ε_iτ_iB_i(T_S)+(1-τ_i)(1+(1-ε_i)τ_i)B_i(T_a)

所有劈窗算法的其主要技术核心都是对上式中的普朗克热辐射函数的简化 处理，

表1 Landsat8遥感反演地表温度实用技术的精度验证结果表

现有技术采取的思路都是对大气传输方程中的普朗克热辐射函数进行统一 处理，而本发明是根据其对地表温度反演的误差影响不同，采取了不同的数学 处理，热红外波段亮温热辐射函数B_i(T_i)通过原始函数公式直接计算，地表温度 热辐射函数B_i(T_S)二次函数拟合处理，大气有效温度热辐射函数B_i(T_a)线性拟合 处理，然后通过解算方程组，即可得出地表温度T_S的计算公式，该过程逻辑简单， 并没有嵌套的数学结构，所用到的数学知识简单易懂，因此大大降低了使用门 槛。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发 明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明 的保护范围之内。

Claims (6)

1.一种Landsat8遥感反演地表温度方法，其特征在于，所述Landsat8遥感反演地表温度方法包括以下步骤：

步骤一：计算大气水汽含量w：利用同期MODIS遥感数据计算大气水汽含量w，计算公式如下所示：

w＝((0.02-ln(b19/b2))/0.651)²

其中，b19和b2对应与MODIS数据中的第19波段和第2波段对应的辐射灰度值；

步骤二：计算Landsat8两个热红外波段对应的大气透过率τ₁₀和τ₁₁：在步骤一计算的基础上，将像元的大气水汽含量w计算大气透过率；

步骤三：计算Landsat8两个热红外波段对应的地表比辐射率ε₁₀和ε₁₁：首先计算目标区域的值被归一化指数NDVI；然后计算Landsat8两个热红外波段对应的地表比辐射率，

步骤四：反演计算地表温度T_S：完成步骤一至步骤三的基础上，将以上参数计算结果代入以下公式计算地表温度值：

上式中各参数是由具体由计算出的大气透过率τ₁₀和τ₁₁、地表比辐射率ε₁₀和ε₁₁计算而得。

2.如权利要求1所述的Landsat8遥感反演地表温度方法，其特征在于，计算大气水汽含量w前，对遥感数据进行预处理、对MODIS影像几何纠正、对Landsat8OLI多光谱遥感数据进行大气校正。

3.如权利要求1所述的Landsat8遥感反演地表温度方法，其特征在于，步骤二中，将像元的大气水汽含量w代入如下公式计算大气透过率：

τ₁₀＝0.9570356-0.0277340w-0.0333734w²+0.0028800w³，

τ₁₁＝0.9456728-0.0857755w-0.0290912w²+0.0032169w³。

4.如权利要求1所述的Landsat8遥感反演地表温度方法，其特征在于，步骤三中，计算Landsat8两个热红外波段对应的地表比辐射率ε₁₀和ε₁₁：

ε₁₀＝0.984p_v+0.964(1-p_v)+(1-0.964)(1-p_v)0.492，

ε₁₁＝0.980p_v+0.970(1-p_v)+(1-0.970)(1-p_v)0.490；

首先计算目标区域的值被归一化指数NDVI，计算公式如下：

然后按照以下公式计算Landsat8两个热红外波段对应的地表比辐射率，

其中，p_v代表植被覆盖度，NDVI代表像元的植被归一化指数值，NDVI_min和NDVI_max分别代表计算空间范围内植被归一化指数的最小值和最大值。

5.如权利要求1所述的Landsat8遥感反演地表温度方法，其特征在于，步骤四中，所述

式中各参数是由具体由计算出的大气透过率τ₁₀和τ₁₁、地表比辐射率ε₁₀和ε₁₁计算而得，具体计算公式如下所示：

A₁₀＝0.0006678ε₁₀τ₁₀，

A₁₁＝0.0006188ε₁₁τ₁₁，

B₁₀＝-0.2333226ε₁₀τ₁₀，

B₁₁＝-0.1990475ε₁₀τ₁₀，

C₁₀＝0.1312942(1-τ₁₀)(1+(1-ε₁₀)τ₁₀)，

C₁₁＝0.1387986(1-τ₁₁)(1+(1-ε₁₁)τ₁₁)，

D₁₀＝211666266ε₁₀τ₁₀-26.7808503(1+(1-ε₁₀)τ₁₀)-B₁₀(T₁₀)，

D₁₁＝16.7224278ε₁₀τ₁₀-27.7043284(1+(1-ε₁₀)τ₁₀)-B₁₁(T₁₁)，

其中，B₁₀(T₁₀)和B₁₁(T₁₁)分别是Landsat8两个热红外波段对应亮温T₁₀和T₁₁的普朗克热辐射值。

6.一种利用权利要求1所述的Landsat8遥感反演地表温度方法的Landsat8遥感反演地表温度系统。