基于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數據，遙感影像清晰無云，質量良好，便于解譯?；贓NVI5.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產品平均值差別較小?？傮w而言，植被指數混合模型法和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.

此外，針對像元尺度的地表發射率驗證，還需比較遙感計算和實地測量的地表發射率值，因此在后續工作中需要進一步拓展，以期更好地為提高反演地表發射率的精度提供線索。