Landsat8热红外波段地表比辐射率估算方法改进

何 好， 张 运,3， 汤 志， 黄 静， 丁 玥

(1.安徽师范大学 地理与旅游学院，安徽 芜湖 241000；2.资源环境与地理信息工程安徽工程技术研究中心，安徽 芜湖 241000；3.自然灾害过程与防控安徽省省级重点实验室，安徽 芜湖 241000)

地表比辐射率是遥感方法反演地表温度的重要参数，作为地物本身内在的一种物理属性，受到地表地物类型，土壤含水量，植被覆盖度与结构等因素的影响。在地表温度反演中，地表比辐射率0.01的误差会导致反演得到的地表温度误差达到1℃以上[1-2]。因此，对地表比辐射率进行准确计算对地表温度反演尤其重要。目前，遥感数据反演地表比辐射率方法主要分为两大类：第一类是直接反演，运用多个热红外波段数据计算地表比辐射率，例如Gillespie等提出的ASTER温度地表比辐射率分离算法[3]、Becker与Li提出的光谱指数方法[4]以及Watson提出的双温法[5]等；第二类是利用归一化植被指数(NDVI)间接估算地表比辐射率，例如Van和Owe提出的植被指数法[6]、Valor和Caselles提出的基于植被指数混合模型的方法[7]、Sabrino等提出的NDVI阈值法[8]与覃志豪等提出的NDVI阈值改进方法[9]。第一类算法需要在8～12μm之间至少3个以上的热红外波段[10]，目前大多数星载热红外传感器都不满足该条件。第二类算法在目前运用相对较为广泛，其中，Van和Owe植被指数法中的经验公式只考虑了自然地表，但在实际应用中，由于热红外波段必然存在纯裸土与纯植被等纯像元，故Van和Owe的植被指数法很难满足目前反演地表比辐射率的精度要求。而Valor和Caselle提出的基于植被指数混合模型的方法，考虑了纯植被像元与纯裸土像元，使地表比辐射率的估算使用范围更加宽广，但是将水体归入纯裸土像元，导致所得到的水体地表比辐射率数值低于正常值。Sabrino等提出的NDVI阈值法除了没有单独考虑水体地表比辐射率之外，其在纯植被像元与纯裸土像元的划分时所设置的NDVI阈值随着遥感影像空间分辨率的提高，不完全适用于Landsat8 热红外波段。覃志豪等提出的NDVI阈值改进方法，对NDVI阈值范围更加精确，但是其将水体提取出来后直接进行赋值0.995，导致所得到的水体地表比辐射率数值上比实际的偏大，且其提出的公式系数适用于TM6波段，该公式是否可以适用于Landsat8 热红外波段，并没有得到验证。目前所提出的地表比辐射率的估算大多还是停留在针对TM/ETM波段得到的公式，其中利用的经验系数与参数也是针对TM/ETM数据。随着Landsat8卫星2013年的升空，由于其具有高分辨率、具有两个热红外波段等特点，越来越多的研究者们，利用Landsat8数据作为自己的数据源进行地表温度的反演。为了便于Landsat8 热红外波段提取地表温度，有必要对Landsat8 热红外波段的地表比辐射率估算进行更深入的探讨。

本文首先介绍植被指数混合模型，在此基础上结合ASTER波谱库数据提出新的估算方法得到Landsat8热红外波段的地表比辐射率，从简单实用的角度考虑，不考虑温度与土壤水分对地表比辐射率的影响。运用MODIS LSE产品、LST产品与地面气象站点实测气温数据对反演结果进行验证，证明本文的计算方法得到的地表比辐射率具有较高精度，可以用于Landsat8热红外波段地表比辐射率的估算,也可以用于地表温度反演研究。

1 方法介绍

1.1 植被指数混合模型

在植被指数法的基础上，Valor和Caselle将表达式由自然像元扩展到地表像元[11]，通过裸土与植被的地表比辐射率以及植被与土壤组合比例，运用NDVI计算出植被覆盖度，之后估算得到地表比辐射率。该方法根据不同地表类型NDVI值不同的特点将地表类型分为裸土、自然地表与植被三种。采用公式(1)计算自然地表像元地表比福射率:

ε=εvPV+εS(1-PV)+dε

(1)

式中，ε为地表比辐射率，εv为纯植被像元比辐射率，εs为裸土像元比辐射率，Pv为植被覆盖度，dε为地形几何形状。dε主要与植被覆盖度、植被距地高度以及植被间距离有关，在地形相对平坦时，可以忽略不计，在地形起伏较大时，dε采用简化公式[12]计算：

dε=4PV(1-PV)

(2)

式中，为地形几何均值，可以根据不同种植被距地高度与它们之间的距离合理取值求得。

1.2 利用ASTER波谱库对模型进行改进

地表比辐射率主要取决于地表物质结构和遥感器的热红外波段区间[13}，会随着热红外波段波长的变化而变化[14]，故估算Landsat8地表比辐射率时，必须将Landsat8的第10波段与第11波段地表比辐射率运用相关参数分别进行计算。从Landsat8遥感影像上观察，可以将地表结构分为5类，分别是水体、纯裸土、城镇、自然表面和纯植被。水面像元、纯裸土像元与纯植被像元统称为纯像元，城镇像元与自然表面称为混合像元，其中城镇像元可以看成由建筑物表面与植被组成，自然表面主要是由土壤、植被等组成。对于地表比辐射率的计算，自然表面与城镇像元是应该关注的重点。

1.2.1 典型地物比辐射率 由于Landsat8 TIRS两个热波段的波长范围不同，在计算Landsat8 TIRS热红外波段地表比辐射率过程中，需要分别利用ASTER波谱库中不同波长范围、不同地表类型与Landsat8 TIRS热红外波段的波谱响应函数得到典型地表比辐射率[15-17]。Landsat8 TIRS典型地表比辐射率数值如表1。水体的比辐射率值：第10波段水体比辐射率ε10w=0.99383，第11波段水体比辐射率ε11w=0.99254。植被的比辐射率相对来说较高，一般认为在0.98～0.99之间[15-17]，利用ASTER波谱库中的植被比辐射率波谱，结合ASTER光谱响应函数求取针叶林、阔叶林和草地的比辐射率，求其三者平均值作为其植被比辐射率，得到第10波段植被比辐射率ε10v=0.98672，第11波段植被比辐射率ε11v=0.98990，由于本次研究不考虑地形、土壤类型与土壤含水量的影响，取土壤平均比辐射率，第10波段裸土比辐射率为ε10s=0.96767,第11波段裸土比辐射率ε11s=0.97790,对于建筑物表面而言，第10波段比辐射率ε10m=0.964885,第11波段比辐射率ε11m=0.975115。

表1 Landsat8热红外波段典型地物地表比辐射率

1.2.2 纯像元地表比辐射率 在地表比辐射率计算时，直接对纯像元地表比辐射率进行赋值计算，其中水体比辐射率分别为ε10w=0.99383、ε11w=0.99254；裸土比辐射率分别为ε10s=0.96767、ε11s=0.97790；植被比辐射率分别为ε10m=0.98672、ε11m=0.98990。

1.2.3 自然表面比辐射率的估计方法 基于植被指数混合模型，不考虑地形起伏与植被结构的影响，假设每一个Landsat8 TIRS波段的自然表面都可以分为植被与裸土两种地物类型的比例组合，则自然表面像元的地表比辐射率可表示为

ε=εvPv+εs(1-Pv)

(3)

式中，εv，εs分别为植被比辐射率与裸土比辐射率，Pv为Landsat8 TIRS每个自然表面像元的植被覆盖度。结合表1数据，可以得到两个热红外波段自然地表比辐射率公式分别为

ε10=0.01905Pv+0.96767

(4)

ε11=0.012Pv+0.97790

(5)

1.2.4 城镇地表比辐射率的估计方法 实际上，Landsat8热红外波段中的城镇像元可以分成建筑物表面与植被的比例组合，基于植被指数混合模型，城镇像元地表比辐射率可以表示为

ε=εvPv+εm(1-Pv)

(6)

式中，εv，εm分别为植被比辐射率与建筑物表面比辐射率，Pv为Landsat8 TIRS每个城镇像元的植被覆盖度。结合表1数据可以得到城镇像元第10与第11波段地表比辐射率分别为

ε10=0.021835Pv+0.964885

(7)

ε11=0.014785Pv+0.975115

(8)

1.3 地物类型分类方法

首先利用Landsat8陆地成像仪中红光波段与近红外波段计算得到影像NDVI值。根据不同NDVI值所对应地物类型不尽相同对地物进行分类，由于水体、云与雪的NDVI数值为负，可以根据其特点将水体提取出来。对于纯植被、纯裸土、自然表面与城镇像元的提取则可以通过植被覆盖度进行计算。运用NDVI计算植被覆盖度，植被覆盖度公式如下：

Pv=(NDVI-NDVIs)/(NDVIv-NDVIs)

(9)

式中，NDVIv，NDVIs分别是研究区植被与裸土的NDVI值。若是遥感影像位置位于山区或者是影像范围内有明显的植被茂密区，则取该植被区域的平均NDVI作为NDVIv值；同样的，若有明显的裸露土壤区域，则取该裸露土壤区域的平均NDVI值作为NDVIs值。当NDVI大于NDVIv时，表明地表类型为纯植被；当NDVI小于NDVIs且大于等于0时，表明地表类型为裸土；当NDVI介于NDVIs与NDVIv两者中间时，表明地表类型为城镇或自然表面。通过公式(4)、(5)与公式(7)、(8)分别计算自然像元与城镇像元的地表比辐射率。在大多数情况下，植被茂密区的NDVI数值都在0.7以上，有时达到0.8甚至更高，裸土的NDVI值的范围一般为0.03～0.08。因此，虽然不同研究区的植被与裸土都有各自的NDVI，NDVIv与NDVIs也有一定大小的差异。但若是研究区上没有明显的植被茂密区与裸土地区，一般可以用NDVIv=0.7与NDVIs=0.05来进行植被覆盖度的计算[18]。

图1 研究区2018年4月10日Landsat8影像

2 研究区与数据介绍

2.1 研究区域

为了对上述提出的计算方法进行验证，选择包含上述5种地表类型区域为研究区。研究区选取安徽省合肥市，合肥市地处中纬度地带，属亚热带季风性湿润气候，年降水量适中，地势相对平坦，平均海拔约为25m，该研究区主要地表类型为城镇，植被，水体，裸土，自然地表。Landsat8影像过境GMT时间为2018年4月10日02时42分。

2.2 数据来源

2.2.1 Landsat8 数据 Landsat8 数据来自美国地质勘探局(USGS)，过境日期为2018年4月10日，云量为0.03%。其中，Landsat8 OLI传感器包含多光谱波段，空间分辨率为30m(其中全色波段空间分辨率为15m)，用于计算NDVI与植被覆盖度，进而计算Landsat8热红外波段地表比辐射率。Landsat8 TIRS热红外波段数据，空间分辨率为100m，用来反演地表温度。

2.2.2 MODIS数据 MOD11A1数据来自于NASA官网(http://earthdata.nasa.gov/)，空间分辨率1000m,包含有地表温度数据MODIS LST与MODIS31，32波段地表辐射率数据MOD11A1_LSE，过境日期为2018年4月10日，GMT时间为02时35分，与上述的Landsat8影像过境时间基本吻合，值得注意的是，MODIS LST数据是按照16位图像进行显示的，在进行地表温度数据验证时需要将DN值乘0.02得到绝对温度。而MOD11A1_LSE数据按照8位图像进行显示，在计算第31,32波段地表比辐射率需将DN值除以256以得到地表比辐射率，并以此来验证本文计算得到的地表比辐射率。

2.2.3 气象站点气温数据 运用本文提出的方法计算出研究区地表比辐射率后，再反演得到地表温度数据，在反演出地表温度后，利用安徽省合肥市6个气象站点气温数据(来自中国气象数据网)进行地表温度验证，以此来验证本文所采用的方法得到的地表比辐射率精度。

3 研究结果与分析

为了评估得到的地表比辐射率的精度，采用两种方法对计算结果进行验证。一是将本文估算得到的LSE与MODIS LSE产品数据进行验证；二是利用Landsat8数据运用本文得到的地表比辐射率数据反演地表温度，然后运用地面气象站点实际气温数据对反演得到的地面温度数据进行精度验证。

3.1 与MODIS LSE产品的对比验证

因为Landsat8热红外的两个波段光谱范围分别与MODIS的31、32波段光谱范围非常接近，具体参数见表2，地表比辐射率又与光谱范围有关，而MODIS LSE产品数据精度较高，故运用MODIS LSE产品对文中计算出的地表比辐射率结果进行对比验证。首先Landsat8 热红外波段空间分辨率与MODIS LSE产品的空间分辨率不同，为了进行更加精确的比较，我们将反演所得到的Landsat8的地表比辐射率数据重采样到与MODIS LSE产品的空间分辨率相同，图2分别为2018年4月10日Landsat8 TIRS第10波段与第11波段利用本文的方法估算出的地表比辐射率影像。从像元整体观察可以发现，Landsat8第10波段与第11波段地表比辐射率基本趋势相同。之后运用结果中最小值、最大值、平均值与标准差进行精度评估，结果如表3。通过表3可以发现，Landsat8热红外波段的LSE与MODIS LSE产品的标准差、平均值、最大值上相差较小，在最小值上Landsat8 TIRS 11波段与MODIS 32波段地表比辐射率相差较多，最小值相差达到0.02，在地表比辐射率平均值上第11波段大于第10波段，通过表1可以基本确定本文计算得到的地表辐射率是符合ASTER光谱库数据的。

(a)第10波段地表比辐射率 (b)第11波段地表比辐射率

表2 Landsat8与MODIS热红外波段参数

为了进一步验证文中计算得到的地表比辐射率的精度，采用平均误差Bias、均方根误差Rmse与相关系数R来评价本文计算得出的地表比辐射率与MODIS LSE地表比辐射率之间的差异，分为Landsat8 第10波段与MODIS第31波段、Landsat8 第11波段与MODIS第32波段两组进行统计对比分析。结果如表3所示。从表3结果中可以得出，第10波段与第11波段同对应MODIS波段的平均误差均小于0.004，均方根误差均小于0.005，相关系数均大于0.8，为显著相关。对比结果表明本文方法计算的地表比辐射率与MODIS LSE产品具有较好的一致性。第10波段光谱范围为10.6μm～11.19μm,第11波段光谱范围为11.5μm～12.51μm，而地表比辐射率在10～12um区间内较为平稳，在此区间以外则变化幅度较大，故第11波段数据相关参数稍大于第10波段。

表3 Landsat8与MODIS热红外波段参数对比

单独从以上指标并不能完全说明本文所反演的地表比辐射率精度，为了验证本文方法的稳定性，利用计算出的Landsat8热红外波段地表比辐射率重采样后与MODIS LSE产品进行差值运算，并把差值分为三种类型：0～0.005，精度较高；0.005～0.01，精度一般；>0.01，精度较低。差值运算如图3。

(a)第10波段与MODIS 31差值图 (b)第11波段与MODIS 32差值图

图中灰色区域是反演精度相对较高的地区，浅绿色是反演精度一般的地区，红色是反演精度较低的地区。从图3观察，Landsat8两热红外波段反演精度相对较高，第10波段最大偏差为0.013，第11波段反演最大偏差为0.026。为了更加直观观察反演精度，我们通过计算像元所占比例来判断反演精度，结果如表4。

可以发现，Landsat8热红外波段利用本文方法计算得到的地表比辐射率精度较高区域在80%以上，一般区域第10波段与第11精度分表为12%与15%，较差区域所占面积分别为4%、5%。总体来说本文估算得到的地表比辐射率精度较高，可以满足目前研究要求。由表3，表4可知，本文估算地表比辐射率方法与MODIS LSE产品具有较高的一致性，精度较高，可适用于Landsat8热红外波段地表比辐射率的估算。

表4 差值结果比例

3.2 基于地表温度的验证

为了进一步验证文中估算地表比辐射率方法的准确性，运用本文估算出的地表比辐射率反演得到研究区的地表温度，并利用地面气象站点气温数据与MODIS地表温度产品(MOD11A1)对反演得到的地表温度数据进行对比验证。目前利用Landsat8 热红外数据反演地表温度的方法有辐射传输方程法[19](RTE)、普适性单通道算法(SC)、基于影像算法(IB)、覃志豪单窗算法(MW)与劈窗算法(SW)。由于因为SC、IB,MW与SW均需要运用相对较多的大气参数，为了防止其他参数对本文计算结果产生影响，故选用参数较少的辐射传输方程法进行地表温度反演。第10波段与第11波段反演结果经重采样后结果如图4。

(a)第10波段地表温度 (b)第11波段地表温度

利用地表比辐射率反演得到的地表温度采用最小值、最大值与平均值对反演结果进行评价，结果如表5，可以知道Landsat8第10波段与第11波段反演结果相对于MODIS LST产品最小值分别差0.9℃、0.02℃，最大值差别分别是-2.85℃、-1.35℃，平均值差别分别为-0.52℃、-0.58℃。从统计指标上来看，Landsat8 热红外波段利用本文计算得到的地表比辐射率反演的地表温度可信度相对较高。

表5 Landsat8地表温度与MODIS LST产品对比统计表

验证文中由地表比辐射率计算出的地表温度精度不可仅仅从统计指标进行验证，需进一步运用实测温度数据进行验证，利用合肥市6个地面气象站点气温数据对反演得到的地表温度进行进一步的验证。结果表明，MODIS LST产品、反演得到的Landsat8 第10波段和第11波段与地面气象站点相比的绝对误差分别为1.35℃、0.82℃与0.84℃。其中MODIS温度产品与气温数据绝对误差相对较大，而Landsat8热红外波段反演出的地表温度同气温数据绝对误差相对较小，这是由于MODIS温度产品为北京时间10时35分过境而Landsat8在北京时间10时52分过境，且当日合肥市无降水，天气晴朗，可知，MODIS LST产品绝对误差较大是符合当日气候特点的。Landsat8卫星过境日期为2018年4月10日，由地-气温差规律[20]可知，地表温度需大于气象站点实测气温数据，而MODIS LST产品与反演得到的地表温度数据平均值均大于气象站点气温平均值，这一点也可以证明地表温度适用性较好。

通过表5、表6对比验证，可以说明利用本文方法计算的地表比辐射率反演得到的地表温度精度与适用性都较好，也可以间接验证该方法求出的地表比辐射率结果较为准确，本文提出的地表比辐射率计算方法可用于Landsat8热红外波段估算地表比辐射率，进而可用于输入参数对Landsat8地表温度进行反演。

表6 气象站点数据、反演得到的地表温度与MODIS温度产品数据对比

4 结语

针对Landsat8 热红外波段地表比辐射率估算方法较复杂的问题，本文基于植被指数混合模型，提出一种适合Landsat8热红外波段并结合ASTER光谱库数据计算地表比辐射率的方法。通过MODIS 地表比LSE产品、MODISLST产品与地面气象站实测气温数据对本文计算出的地表比辐射率进行对比验证后得出该方法精度较高，可以满足目前对Landsat8 热红外波段估算地表比辐射率的精度要求，利用本文计算得到的地表比辐射率反演出的地表温度也可以满足地表温度反演的精度要求。在今后对Landsat8进行地表比辐射率与地表温度反演时，可以更加简单精确得到实验结果。但只对合肥市进行了验证与分析，在其他地方还需要进行更广泛的对比验证，此外研究区地形起伏相对较小,在此后的研究中还需要进一步考虑地形起伏等因素对地表比辐射率的影响。