基于Landsat 8数据反演地表发射率的几种不同算法对比分析

王丽霞，孙津花，刘 招，张双成，杨 耘

（1.长安大学地质工程与测绘学院，陕西西安710054；2.国土资源部退化及未利用土地整治工程重点实验室，陕西西安710075；3.长安大学地球科学与资源学院，陕西西安710054；4.长安大学 环境科学与工程学院，陕西西安710054）

0 引 言

地表发射率是遥感数据获取陆面温度中的一个关键参数，它与地表组成成分、地表粗糙度、含水量等因素有关，普遍应用于辐射传输过程和地气系统的能量平衡［1］。热红外遥感传感器所接收到的信号包括2部分，一部分是被测物体发射的热辐射，另一部分是由于地物的发射率小于黑体而反射了环境辐射的冷发射［2］。因此，为了求得地物的真实温度，还必须知道地物的发射率。一些学者进一步研究发现，地表发射率有0.01的变化，估算出的地表温度就会产生1～2℃的误差［3-6］。因此，地表发射率的准确获取，对地表温度反演精度的提高至关重要。

目前，利用RS反演区域地表发射率是获取大范围地表发射率的重要途径［7］。张成才等采用Van和Owe提出的植被指数法，建立经验公式并结合典型地表发射率值计算得到河南省白沙灌区的地表发射率［8］。林志斌采用三波段梯度差法简单估算了植被覆盖度，并选取适合Landsat ETM+影像的植被指数混合模型算法计算得到了地表发射率［9］。唐伯惠等提出基于分类的地表发射率反演算法，由于地表状况本身的复杂性以及卫星观测中的混合像元问题，发射率反演精度通常达不到所需的0.01［10］。考虑到NDVI阈值法的可操作性好及较高的精度，赵少华等应用Sobrino等提出的NDVITEM算法获取山西省地表发射率，并反演得到较高精度的地表温度值［11］。覃志豪等对Sobrino等提出的NDVITEM算法进行了改进，考虑了不同地表形态下dε（粗糙表面的空腔效应引起的发射率比例）取值以及城市像元发射率估算的问题，并对山东省陵县区进行地表发辐射率反演，取得了较合理的地表温度反演结果［12-13］。与之前的Landsat系列卫星相比，Landsat8自2013年发射以来，在许多方面进行了改进，也被广泛应用于地表温度的反演［14-16］。然而，仅针对Landsat 8数据进行地表发射率反演算法的对比研究较为少见。

归一化植被指数（Normalized Difference Vegetation Index，NDVI）是描述地表覆被生长状况的一个重要指标［17-18］。由于地表热红外波段的发射率同NDVI之间具有非常高的相关性，目前基于NDVI获取地表发射率的方法主要有4种。即：Van和Owe提出的植被指数法，Valor和Caselles提出的植被指数混合模型法，Sobrino等提出的NDVITEM算法，覃志豪等提出的NDVITEM改进算法［19-21］。4种算法的表达式及基本参数设定都不相同。因此文中利用西安市的Landsat8数据，应用ENVI，ArcGIS等软件，分别用4种算法对研究区进行了地表发射率反演，并结合同时期精度0.01的MODIS LSE产品，对比分析了4种算法的适用性和精度结果，以期获取更精确的地表发射率反演算法。

1 数据来源及处理

1. 1 研究数据

西安市地处西北部关中盆地中部秦岭北麓，地跨渭河南北两岸，位于北纬 33°39′～34°45′，东经107°40′～109°49′.下辖 11区 2县并代管西咸新区，总面积 10 752 km2，2017年末户籍人口905.68万。属暖温带半湿润季风气候，平均海拔高424 m，1月份平均气温0.4℃，7月份平均气温26.6℃，年平均气温13.3℃，年平均降水量613.7 mm，年平均湿度69.6%［22-23］.研究数据采用2015年4月28日西安市的Landsat8数据，遥感影像清晰无云，质量良好，便于解译。基于ENVI5.3平台对遥感影像进行的前期处理有：影像镶嵌和裁剪；OLI数据进行辐射定标；利用FLAASH模型对OLI数据进行大气校正。为了更好地比较不同算法、不同下垫面类型反演地表发射率的差异情况，文中使用分类回归决策树方法将地表分成植被、建筑物、裸土和水体4类。在ENVI 5.3光谱指数计算器中分别获取NDVI（归一化植被指数），NDBI（归一化建筑指数），MNDWI（改进型归一化水体指数），结合人工判定阈值（根据10%置信度，NDVI＞0. 65为植被；NDBI＞-0. 05为裸土；MNDWI＞0. 17为水体）建立决策树模型。地表分类结果如图1所示。

图1 西安市决策树分类结果Fig.1 Xi’an decision tree classification results

1. 2 对比数据

文中使用了MODISMOD11A1地表发射率产品进行Landsat 8反演的地表发射率结果精度验证。MOD11A1是经过正弦投影后的全球每日地表发射率产品，该产品基于地表下垫面的分类思想，利用地表分类信息，通过模型计算实现了全球地表发射率产品的生成［9］。在晴朗天气条件下，MOD11A1地表发射率产品精度可达到0.01.

2 地表发射率反演算法

2. 1 植被指数法

根据NDVI和地表发射率之间有高度的相关性为理论基础，将博茨瓦纳采集的发射率信息与NDVI相关联，建立理论模型。模型如式（1）所示

式中，经验公式是基于博茨瓦纳热带草原总结的，且 NDVI的适用范围为0. 157～0. 727.

2. 2 植被指数混合模型法

该方法的主要思想是：假定植被和裸土的地表发射率已知分别为0. 985和0. 96. 当NDVI＜0. 1时，认为是纯裸土像元；当NDVI＞0. 72时认为是纯植被像元；NDVI介于0. 1至0. 72之间时，认为是混合像元，采用式（2）计算地表发射率

式中 ε为地表发射率；εV为纯植被像元发射率；

式中 NDVI为归一化植被指数；NDVImin=0. 1，NDVImax=0. 72，且当 NDVI＞0. 72时认为地表覆被状况良好，则植被覆盖度赋值为1，当NDVI＜0. 1时认为地表覆被状况较差，则植被覆盖度赋值为0. 发射率比例dε计算公式如式（4）

式中，＜dε＞为地形几何形状因子，根据不同几何分布可取均值0. 015.

2. 3 NDVI TEM算法

Sobrino等在前人研究的基础上，对植被指数混合模型法进行了改进。其主要思想是：当NDVI＜0. 2时，认为地表覆被状况较差，地表发射率为ε=0.98-0. 042ρ，式中ρ是红外波段的地物反射率；NDVI＞0. 5时认为是地表覆被状况良好，地表发射率为常数 0. 99；NDVI在 0. 2至0. 5之间时认为是混合区，通过式（5）计算地表发射率

式中 ε为地表发射率；εV为纯植被像元发射率；εS纯裸土像元发射率；PV为植被覆盖度；dε为考虑了粗糙表面的空腔效应引起的发射率比例。植被覆盖度PV计算公式（6）εS为纯裸土像元发射率；PV为植被覆盖度；dε为考虑了粗糙表面的空腔效应引起的发射率比例。植被覆盖度（PV）计算公式如式（3）

式中 NDVI为归一化植被指数；NDVImin=0. 2，NDVImax=0. 5，且当 NDVI＞0. 5时认为地表覆被状况良好，则植被覆盖度赋值为1，当NDVI＜0. 2时认为地表覆被状况较差，则植被覆盖度赋值为0.

式中，F为值为0. 55的形态参数。

2. 4 NDVI TEM改进算法

覃志豪等考虑了不同地表形态下dε的取值以及城市像元发射率估算的问题。该方法认为：对于水体像元，一般取 εw=0. 995来估计；PV≤0.5时，dε=0. 003 796PV；PV＞0. 5时，dε=0. 003 796（1-PV）. 发射率估算公式如下：自然表面像元发射率计算如下式

城镇像元发射率计算如下式式中 εV，εS和εm分别为纯植被、纯裸土和纯建筑表面像元的发射率，分别取 εV=0. 986，εS=0.972，εm=0. 970. RV，RS和 Rm分别为纯植被、纯裸土和纯建筑表面像元的温度比率，与植被覆盖度有关，计算公式分别如下

3 反演结果与对比分析

3. 1 反演结果

基于以上4种算法，利用ENVI 5.3图像处理软件，对西安市影像进行地表分类，然后计算NDVI及植被覆盖度，最后反演得到地表发射率值。4种算法计算结果如图2所示。

图2 Landsat8 4种算法发射率反演结果Fig.2 LSE obtained from Lansat8 data using four LSE retrieval algorithms

3. 2 对比分析

3.2.1 4种算法反演结果与MODISLSE产品比较

4种算法不仅在反演模型上有差异，对PV为0和1时NDVI的设定也不同。这样的差异随着研究区地表覆被状况的不同，可能会产生不同的影响。因此有必要对4种算法进行比较，进而得到适合本研究区的算法。

Landsat8热红外数据与MODIS发射率产品的空间分辨率不同。为了进行便精确之间的比较，首先利用 MRT（MODIS Reprojection Tool）工具将MODIS LSE产品地图投影转换为UTM-WGS84，并进行研究区范围的裁剪，然后利用ENVI软件采用最邻近内插法将Landsat8发射率反演结果重采样到1 000m.文中采用最小值、最大值和平均值这些统计指标，对反演结果进行评价，统计结果见表1.

表1 Landsat 8各算法反演结果与MODIS发射率产品对比Table 1 Comparation between Landsat 8 retrieval results and MODIS LSE product

从西安区域的统计来看，利用4种算法反演得到的发射率结果，相对于MODIS LSE产品最小值差别分别为0.035 6，0.002 0，0.010 5，0.0003，最大值差别分别为 0.000 6，0.005 1，0.004 1，0.001 8，平均值差别分别为 0.012 2，0.000 5，0.005 2，0.001 4.其中，植被指数混合模型法和NDVITEM改进算法与MODIS LSE产品最小值差别较小，植被指数法和NDVITEM改进算法与MODIS发射率产品最大值差别较小，植被指数混合模型法和NDVITEM改进算法与MODIS LSE产品平均值差别较小。总体而言，植被指数混合模型法和NDVITEM改进算法与MODIS LSE产品更接近，可信精度较高。

3.2.2 4种算法之间的比较与分析

为了进一步对比4种反演算法在区域尺度的差异，文中对利用各种算法反演的2015年4月28日西安市Landsat 8数据发射率结果进行比较分析，图3是研究区内4种算法反演结果两两之间的差值统计图。

图3 4种方法反演结果两两之间的差值统计Fig.3 Statistical analysis of the differences between the inversion results of fourmethods

从图3可以看出，植被指数混合模型法和NDVITEM改进算法反演结果差异最小，差值范围控制在0.04之内，且大部分像元发射率差值靠近0值。NDVITEM算法与NDVITEM改进算法反演结果差异最大。主要原因是：两者在NDVI最大值和最小值的设定上有所不同，NDVITEM算法中，NDVImin=0. 2，NDVImax=0. 5. NDVITEM改进算法中，NDVImin=0.05，NDVImax=0. 7，且 NDVITEM改进算法对水体像元直接赋给0. 995的发射率值，而NDVITEM算法在计算NDVI小于0. 2的像元时，利用红外波段的反射率计算其发射率值。

为揭示不同算法、不同地表覆被类型反演的地表发射率的差异情况，文中还基于上文决策树分类所提取的植被，城镇，水体，裸土4种地表类型，对不同算法不同地表覆盖类型的发射率平均值两两之间的差值进行分析见表2.

从表2可以看出，在分类尺度上，对于植被区域，植被指数法与NDVITEM算法差异最大达到了-0. 009 4，植被指数混合模型法与 NDVITEM改进算法差异最小，仅0. 001 2，对于城镇区域（未考虑植被指数法），植被指数混合模型法与NDVITEM算法差异最大，达到了0. 011 7，NDVITEM算法与NDVITEM改进算法差异最小，达到了0. 002 9，对于水体区域（未考虑植被指数法），NDVITEM算法与NDVITEM改进算法差异最大，达到了0. 033 3，植被指数混合模型法与NDVITEM改进算法差异最小，达到了-0.009 9，对于裸土区域（未考虑植被指数法），NDVITEM算法与NDVITEM改进算法差异最大，达到了0. 012 0，植被指数混合模型法与NDVITEM算法差异最小，达到了0.005 9.总体而言，地表类型为植被的情况下，算法之间的差异最小，其次为城镇和裸土，地表类型为水体的情况下，算法之间的差异最大。

表2 4种方法反演的不同植被类型的LSE平均值两两之间差值Tab. 2 Differences of average LSE values for different land coverage between pairs of four LSE retrieval results

4 结 论

1）从像元尺度上，植被指数混合模型法与NDVITEM改进算法总体精度相对较高，植被指数法与NDVITEM算法总体精度次之；

2）从下垫面类型尺度上，对于植被区域，4种算法之间的差异最小，其次为城镇和裸土。对于水体区域，4种算法之间的差异最大；

3）从算法的适用性而言，植被指数法适用性最弱，该模型只能用来反演NDVI值范围为0. 157～0. 727的自然表面，且反演精度较低。植被指数混合模型法与NDVITEM改进算法适用性较好，且算法之间的差异最小。这与两者模型中NDVI最大值和最小值的设定上有关，植被指数混合模型法中，NDVImin=0. 1，NDVImax=0. 72；NDVITEM改进算法中，NDVImin=0. 05，NDVImax=0.7.

此外，针对像元尺度的地表发射率验证，还需比较遥感计算和实地测量的地表发射率值，因此在后续工作中需要进一步拓展，以期更好地为提高反演地表发射率的精度提供线索。