基于Landsat 8和MODIS Terra卫星数据的多种单通道反演算法对比

赵 阳，亓进磊，黎 慧

阜阳师范大学历史文化与旅游学院，安徽阜阳,236037

陆地表面温度LST(Land Surface Temperature，地表温度)作为常见的特征物理量，在地表与大气能量交换过程中扮演着重要角色，在城市热环境、地表辐射能量平衡、全球气候变化和资源环境监测等领域都有着重要的研究价值[1-3]。较为传统的LST获取的方法需要耗费大量的人力物力，多是以点的数据形式对附近温度进行推演，如果想要得到较为准确的大范围LST数据，需要区域内大量数据观测点的支持。近几十年间，遥感卫星影像的分辨率和精度都得到了大幅度的提升，应用领域也在不断拓宽。热红外遥感探测技术的发展更是使得大区域温度的获取和研究区域时空变化成为了可能，通过热红外遥感获取地表热红外波谱的辐射能量信息，根据地表不同物体反射率的差异，使用相关的遥感影像处理软件可以快速方便的获取其热力学温度，在LST反演方面发挥着重要作用[4-5]。遥感影像一般是通过无量纲的数字量化值(Digital Number，DN)记录的，不同地物反映到遥感卫星传感器上的DN值有所不同。在对遥感影像数据做定量化研究时，可以通过辐射定标实现DN值与辐射亮度值、比辐射率和亮度温度等物理量的转化。基于热红外遥感的LST反演与应用取得了一系列的研究成果，并取得了很好的温度反演效果。

目前LST的遥感反演方法有单通道算法、多通道算法、多角度算法、多时相算法和高光谱反演算法等[1]。其中较为常用的是单通道算法和多通道算法，国内外的众多学者也相继提出一系列的改进型算法，并在验证过程中取得了很理想的温度反演效果。本文主要使用Landsat 8数据与MODIS卫星数据，而美国地质调查局USGS(United States Geological Survey)指出Landsat 8 TIRS传感器存在一定的问题，如第11波段会出现定标不稳定等，在对其进行定量研究时，会出现影像数据与真实温度的差异，不建议使用(USGS，2014)。而分裂窗温度反演算法过程中需要分别对两个热红外影像进行辐射定标和大气校正，如果使用第10、11波段进行分裂窗算法，其反演结果误差较大。因此，本文只利用Landsat 8第10波段热红外波段进行多种单通道温度反演，以确保温度反演结果的准确性。

本文主要使用四种算法，即大气校正法[6-8]、普适性单通道算法[9]、覃志豪单窗算法[10]和胡德勇TIRS10_SC单通道算法[3]。通过对这四种单通道反演算法的研究对比，得到较为符合本研究区的温度反演算法，为区域内的城市热岛、温度时序变化研究等提供依据。工作流程如图1所示。

1 研究区和数据源

1.1 研究区

阜阳市位于淮北平原腹地，安徽省西北部，地跨东经115°37′58″～115°59′46″、北纬32°46′17.6″～33°02′32.5″，下辖三区四县一市，总面积9 775 km2，2016年统计年鉴显示阜阳市户籍总人口1 061.55万人，常住人口799.1万人，是安徽省人口最多的地级市。由于Landsat 8 遥感影像数据过境时间与图幅镶嵌裁剪等出现的误差，为了确保反演数据的精度与准确性，本文只以阜阳市市区为研究区域，研究区内地形平坦，地面平均标高海拔一般为28～31 m，主要的土地利用类型为建设用地、耕地、裸地等。

1.2 数据源

自从Landsat系列卫星上空后，在获取免费方便、数据持续更新时间长、中等分辨率和数据信息成熟等多方面优势下，使其在全球的教育教学和科研等方面取得重要地位。MODIS为低分辨率卫星，凭借其数据获取周期性短、获取时间固定等优势，在全球大气监测、大范围宏观温度反演等方面有着天然的优势。

1.2.1 遥感数据及预处理

本文共使用了两种遥感影像数据，具有相近过境时间的Landsat 8遥感卫星数据和MODIS遥感卫星数据(如表1)。使用的Landsat 8数据为2018年9月8日的Landsat8 OLI_RIRS C1 Level-1数据产品(经过几何校正的Landsat 8一级数据产品)，图幅编号为LC81220372018251LGN00，数据条带号为122，行编号为37，时间为当地时间10:48，选择此数据的主要因素是无云层覆盖有利于后期地表反射率的反演数据的准确性。Landsat 8上携带两个传感器，在本次研究中主要使用OLI的第3、第4、第5、第6波段数据和TIRS第10波段TIR1数据。MODIS数据主要采用的是与Landsat 8遥感数据时刻接近的MODIS Terra卫星产品数据，主要使用MOD021KM产品数据进行几何校正和大气水汽含量以及大气透过率的估算，数据时间为当地时间11:55，当天晴朗无云且与Landsat 8数据时间相隔67分钟，由于当天天气晴朗无云，可以认为这两个数据的水汽含量相似，为下文的大气透过率的反演精度提供了保障。

表1 遥感影像数据相关参数

本文遥感数据预处理主要是通过ENVI、ArcMap等软件实现的，其中对Landsat 8的OLI传感器波段进行辐射定标和FLAASH大气校正，得到与实际地表相同的地表反射信息，然后根据校正好的遥感影像信息通过分类回归决策树CART(Classification And Regression Tree)对地表典型地物进行分类，最后根据分类结果对地表比辐射率进行估算；对第10波段TIR1数据进行辐射定标，然后通过普朗克函数得出其亮度温度值。值得注意的是由于使用的是经过几何校正的Landsat 8一级数据产品，第10波段TIR1数据已经被重采样为30 m的空间分辨率(与OLI波段分辨率保持一致)，这里不需要对其进行几何校正。MODIS数据主要是通过第2波段与第19波段来估算研究区的大气水汽含量，在进行估算前需要对数据进行配准重采样处理，使其与Landsat 8数据具有相同的分辨率和投影坐标信息。

1.2.2 实测数据

随着现代技术的发展，传统的气象收集手段也发生着改变，气象现代化不仅可以节省大量的人力物力让大众获取更加实时准确的气象信息，还能在现代大数据的支持下让数据预测变得触手可及。至2015年，阜阳市全市共建成自动气象站157个，气象观测网络基本覆盖辖区内乡镇，通过无线传感器网络可以实现对地表气象数据全天候连续不断采集，主要的气象信息包括气温、风向、降水量等，这些气象站点对阜阳市的气象准确的监测与预测起着重大的作用[11]。2017年阜阳市气象局阜阳自动气象站数据显示平台(http://fy.ahqx.gov.cn/amos/table.php)的开通，为地面实时实测气象数据的方便获取提供了可能。

本次研究所使用的Landsat 8遥感卫星过境时刻为当地上午10:48分，而阜阳市自动气象站上传频率为每十分钟一次，自动气象站实测数据基本上可满足本次研究需求，将以研究区内的22个可用实测点数据与四种算法最终结果做对比分析。

2 四种单通道地表温度反演算法

2.1 大气校正法

辐射传输方程(大气校正法)是根据辐射传输方程模型提出来的，是整个LST反演算法理论体系的基础，作为最早的温度反演算法之一，其具有很强的适用性，可以广泛地应用于各个热红外数据。

辐射传输方程法在使用的过程中需要使用相关的大气模型软件(LOWTRAN、MODTRAN或6S等)来模拟实时的大气数据(包括不同高度的气温、气压、CO2含量等)，根据不同物理要素对地表热辐射的不同影响，模拟得出数据相关路径上的大气透过率、大气下行辐射和大气下行辐射等大气剖面参数，然后再从遥感卫星传感器接收到的辐射总量中减去受到这三部分影响的大气影像，得到真实地表辐射值，如果在地表辐射率已知的情况下，通过普朗克函数求逆可以获取真实的地表温度值，相关的公式如下：

Lλ=[εB(TS)+(1-ε)L↓]τ+L↑

B(TS)=[Lλ-L↑-τ(1-ε)L↓]/(τ·ε)

TS=K2/ln[K1/B(TS)+1]

根据辐射传输方程，卫星传感器接收到的能量值Lλ(通过定标的热红外辐射亮度值)与ε地表反射率、B(TS)同温度下黑体热辐射亮度值、L↓大气向下辐射亮度、L↑大气向上辐射亮度和τ大气在热红外波段的透过率(大气透过率)等有关。地表的真实温度Ts(单位K)可以根据普朗克公式求逆得出。K1、K2为两个常数值，对于Landsat 8数据而言，可以直接通过其元数据文件直接查到其对应的K1、K2值，第10波段TIR1中K1=774.885 3 W/m2·μm·sr，K2=1 321.078 9 K。

2.2 覃志豪单窗算法

基于辐射传输方程的LST反演算法计算过程较为复杂且参数具有很强的不确定性，需要遥感卫星过境时地区较为准确的大气轮廓信息数据，而大气轮廓信息的获取复杂且成本较高，一般只是使用全球标准大气轮廓信息代替，导致LST反演产生的误差较大。

单窗算法是覃志豪等[10]根据地表热辐射传导方程，在一系列的假设条件下，推导出一种适用于Landsat TM热红外波段影像简单可行、反演精度较高的LST反演算法。该算法把大气和地表的影响因素直接考虑到反演公式当中，根据大气透过率、地表比辐射率和大气平均温度三个参数进行LST的推算，根据之间相对应的关系，将单窗算法表达为大气透过率和大气平均作用温度的函数。公式如下：

TS=[a(1-C-D)+(b(1-C-D)+C+D)Tsen-DTa]/C

Ts为地表的真实温度(单位K),C和D是地表反射率(ε)和大气透过率(τ)所构成的中间量，其中C=ε·τ，D=(1-τ)[1+(1-ε)τ]，Tsen是卫星传感器所探测到的像元亮度温度(单位K)。a和b是根据回归分析得到的两个常量，在不同的温度变化范围下，表现的系数值和相关系数也有所不同。国内外众多的学者[12]根据大气辐射传输软件LOWTRAN模拟Landsat 8不同温度范围内TIR1(第10波段)的反演回归系数得出：温度范围越小所对应的相关性越强。根据阜阳市夏季实际情况，选择a10=-62.806，b10=0.434，相关系数的平方(r2)为0.999 2。大气平均作用温度Ta与近地表温度T0(一般距地面2 m处)之间呈一定的函数关系，根据热红外波段大气平均作用温度估算方程，在中纬度夏季大气模式下，Ta=16.011 0+0.926 21T0。

2.3 Jiménez-Muoz普适性单通道算法

Ts=γ[(ψ1Lsen+ψ2)/ε+ψ3]+δ

γ≈Tsen2/(bλLsen)

δ≈Tsen-Tsen2/bλ

Ts为地表真实温度(单位K)；Lsen是传感器探测的辐射亮度值；ε地表反射率；Tsen是卫星传感器所探测到的像元亮度温度(单位K)，λ是遥感影像热红外波段的中心波长；bλ为常数，在TIRS10波段为1 324；ψ1、ψ2、ψ3是大气功能参数，与大气水汽含量有关，在大气参数(τ、L↓、L↑)已知的情况下，大气功能参数可用式计算：

ψ1=1/τ,ψ2=-L↓-L↑/τ,ψ3=L↓

一般情况下，实时的大气参数很难获取，Jiménez-Muoz等人选取GAPRI数据库中的4 838条大气轮廓数据，通过大气辐射传输软件MODTRAN4.0模拟计算大气剖面参数，得出相应的系数矩阵，推出Landsat 8中三个大气函数与大气水汽含量之间的二次多项式。具体方程式如下，其中ω为大气水汽含量：

ψ1=0.040 19ω2+0.029 36ω+1.015 23

ψ2=-0.383 33ω2-1.502 94ω+0.203 24

ψ3=0.009 18ω2+1.360 72ω-0.275 14

2.4 TIRS10_SC算法

TIRS10_SC算法是胡德勇等人[3]综合TIRS 10特性和热辐射传输方程，建模地表温度和亮度、大气平均作用温度、大气透过率和地表发射率等参数之间的关系，提出的针对Landsat 8 TIRS第10波段的单窗算法。

φ1=ε·τ

φ2=(1-τ)[1+(1-ε)τ]

T10为TIRS 10的亮温值(单位K)；K2可以通过元数据文件查到，对于Landsat 8，K2=1 321.078 9 K；Ta为大气平均作用温度，通过表根据当地气象资料可以得到。

3 单窗算法参数的量化

3.1 亮度温度的推算

当某一波段下一个物体的辐射亮度与绝对黑体的辐射亮度值相同时，则该黑体的物理温度就称为该波段下该物体的亮度温度(简称“亮温”)，所以亮温只是衡量物理温度的一个指标，并不能代替物体真实的物理温度值。亮温计算只要分为两个步骤，首先要将卫星获取的影像像元DN值转化为与之相对应的大气层上界热辐射强度(Ii)，然后根据Planck函数利用热辐射强度推算出该物体的亮温值，公式如下：

Ti=Ki2/ln(1+Ki1/Ii)

Ti为物体的亮度温度，i表示热红外遥感第几波段，Ki1、Ki2分别代表两个常数值，对于Landsat 8第10波段TIR1数据而言，K10 1=774.885 3 W/m2·μm·sr，K10 2=1 321.078 9 K，可用通过ENVI软件中的定标工具(Radiometric Calibration)直接将遥感影像的DN值定标为亮度温度值。

3.2 大气水汽含量反演

水汽是影响大气透过率的重要因素，在LST反演过程中，大气水汽含量估算的准确度直接影响到最后研究结果的精度。由于实时的大气剖面数据获取困难，水汽含量一般直接通过MOTRAN、6S等大气模型软件用标准大气进行模拟。由于Landsat 8数据水汽反演的不成熟与其局限性，本文采用了与Landsat8卫星数据过境时间相近的MOD05_L2数据进行大气水汽含量的反演。

MODIS包含36个波段，其中17、18和19波段为大气吸收波段，2和5波段为大气窗口波段。Kaufman等[13]指出使用MODIS数据通道比值法进行大气水汽含量的反演，可以部分有效地去除地表反射率随波长变换对大气透过率产生的影响，提高反演精度。根据毛克彪研究结果[14-15]，选取MOD05_L2第2波段和第19波段数据进行通道比值法来反演大气水汽含量，具体公式如下：

ω=[α-ln(ρ19/ρ2)/β]2

其中，ω表示大气水汽含量(g/cm2)；ρ2、ρ19分别表示MODIS第2波段和第19波段的地表反射率；α和β是两个校正系数，根据阜阳市的实际情况，本文选择混合型地表的参数，其中α=0.02，β=0.651(Kaufman等[13])。根据公式与相关系数计算得出阜阳市的大气水汽含量影像，如图2。

图2 大气水汽含量影像

3.3 大气透过率反演

Rozenstein等[12]通过LOWTRAN软件模拟分析Landsat 8 热红外波段在不同水汽范围内大气透过率与水汽含量之间的关系，得出在不同大气模式下Landsat 8热红外波段大气透过率估算方程，如表2所示。

表2 Landsat 8 0.5-3.0g/cm2水汽范围内大气透过率与水汽含量的关系

本研究区域范围在北纬32°～33°，属于中纬度地区，且研究时间位于夏季，所以采用中纬度夏季的大气模式，得出第10波段的大气透过率与大气水汽含量之间的关系：τ10=-0.113 4ω+1.033 5。根据大气水汽含量得出大气大气透过率影像，如图3。

图3 大气透过率影像

3.4 地表反射率的估算

在LST反演过程中，地表反射率的估算是反演的重点，主要有差值法、独立温度光谱指数法(TISI)和NDVI门槛值(NDVITHM)等方法。对于Landsat 系列影像图像，地表反射率主要取决于地表的物质结构和遥感的波段区间。地球表面不同区域的地表结构虽然很复杂，但从卫星像元的尺度来看，可以大体视作由4 种类型构成:水面、建筑、裸土和植被。建筑包括城市和农庄，主要由道路、各种建筑和房屋组成;植被包括林地和农田等。使用回归决策树CART方法，通过各种归一化的指数(NDVI、MNDWI和NDBI等)阈值对地表进行分类，如图4所示。

图4 决策树模型分类流程图

为了得到更精确的地表反射率数据，宋挺等[12]根据不同典型地物的比辐射率特征以及参照ASTER提供的常用地表比辐射率光谱库(http://speclib.jpl.nasa.gov)设定landsat 8第10通道比辐射率值如下：ε10w=0.996 83，ε10v=0.986 72，ε10s=0.967 67，ε10m=0.964 885，ε10w、ε10v、ε10s和ε10m分别表示第10波段下典型水体、植被、裸土和建筑的地表比辐射率(图5)。根据覃志豪等[15]提出的地表比辐射率的估算方法，将数据像元分为自然表面像元、城镇像元等分类进行运算，算出不同混合像元的地表比辐射率，具体的公式如下：

图5 分类影像

自然表面像元：εi=PvRvεiv+(1-Pv)Rsεis+dε

城镇像元：εi=PvRvεiv+(1-Pv)Rmεim+dε

式中，Pv为植被覆盖率，Pv=(NDVI-NDVIs)/(NDVIv-NDVIs)，其中NDVIv与NDVIs分别为植被和裸土的NDVi值，在没有明显的完全植被或裸土像元时，本文用NDVIv=0.7和NDVIs=0.05来近似估算植被覆盖度(图6)。在地表高低相差较大的情况下，可取dε=0。由于热辐射相互作用在植被与裸土各占一半时达到最大，提出如下经验公式来估算dε：

图6 地表比辐射率影像

当Pv<0.5时，dε=0.003 8Pv;当Pv>0.5时，dε=0.003 8(1-Pv)

当Pv<0.5时，dε最大，dε=0.001 9

Rv、Rm、Rs分别表示植被、建筑物和裸土的温度比率，覃志豪等根据主要地表类型的温度差异进行模拟，得出自然表面比辐射率的估算经验公式：

Rv=0.933 2+0.058 5Pv

Rm=0.988 6+0.128 7Pv

Rs=0.990 2+0.106 8Pv

4 地表温度反演结果对比

反演算法精度可以通过获取地面实测值或其他标准值进行验证，本文利用阜阳市自动气象站的实测数据来验证4种算法反演地表温度的精度。逐一计算各个参数之后，进行4种算法的地表温度反演，反演结果影像(取置信区间99.9%以去除差异值)如图7-10所示。

根据2018年9月8日4种算法温度反演结果与地面实测温度数据对比(如表3所示)可以看出，大气校正法、Qin_SC、J&M_SC10和TIRS10_SC的温度反演结果与地面实测温度相比的绝对误差的平均值分别为0.95 ℃、2.04 ℃、0.94 ℃和2.09 ℃。

表3 2018年9月8日4种算法温度反演结果与地面实测温度数据对比

如表4所示，大气校正法与J&M_SC10算法反演结果相近且误差相对较小，两种算法的最大绝对误差都为2.58 ℃，大气校正法与J&M_SC10算法绝对误差的最小值分别为0.05 ℃和0.06 ℃；而Qin_SC算法与TIRS10_SC算法的误差较大，Qin_SC算法的最大绝对误差为4.61 ℃，TIRS10_SC算法的最大绝对误差为4.78 ℃。通过对比不难发现，大气校正法和J&M_SC10算法均比较适合阜阳市地区的地表温度反演，反演结果与实际温度相差较小。四种算法反演温度与温度站实测数据不尽相同，尤其是Qin_SC算法与TIRS10_SC算法明显高于实测地表温度，由图也可以看出，Qin_SC算法与TIRS10_SC算法比大气校正法与J&M_SC10算法明显高2 ℃左右，分析原因主要有：(1)相关算法参数的不同。大气校正法与J&M_SC10算法所考虑的参数大致相同，反演结果相近，TIRS10_SC算法是辐射传输方程与Qin_SC算法的改进型，所以TIRS10_SC算法与Qin_SC算法考虑的因素参数较为相近，从而导致了四种算法反演温度两两相近的结果；(2)数据时间的差异。实测数据为当地时间10:50，Landsat 8影像数据为当地时间10:48，MODIS影像数据为当地时间11:55，时间相差所导致的反演精度变化在所难免；(3)算法计算参数的不确定性。在参数的计算过程中，很多的经验公式只是提供典型的夏冬季节，没有更多更加精确的时间经验公式选择，从而加大了各项参数在估算过程中的误差导致最后算法LST反演结果的误差增大。

表4 4种算法最大值、最小值和平均值对比

5 结 论

本文使用阜阳市2018年9月8日阜阳市地区的Landsat 8遥感数据和MODIS数据，通过大气校正法、覃志豪单窗算法(Qin_SC)、Jiménez-Muoz 普适性单通道算法(J&M_SC10)和胡德勇单窗算法(TIRS10_SC)四种算法进行LST反演对比。通过四种反演算法结果与实测数据对比，可以明显看出：大气校正法与J&M_SC10算法LST反演结果较为理想，绝对误差的最小值分别为0.05 ℃和0.06 ℃，与地面实测温度相比的绝对误差的平均值分别为0.95 ℃、0.94 ℃，比较适用于阜阳市地区；而Qin_SC算法与TIRS10_SC算法LST反演结果误差较大，与地面实测温度相比的绝对误差的平均值均大于2 ℃，与地面实测数据有较大的出入，在用于阜阳市地区时需要对其进行改进。

本文只是使用四种算法的计算公式做对比，其中所有的参数计算都是基于其中的经验公式，没有使用相关的大气模型软件(LOWTRAN、MODTRAN或6S等)对阜阳市这一特定地区进行大气参数的模拟，公式中所使用的相关的参数值也无法与阜阳市的具体实际情况相匹配，导致研究有些许误差，这是以后改进的重点。