基于Landsat 8數(shù)據(jù)反演地表發(fā)射率的幾種不同算法對(duì)比分析

王麗霞，孫津花，劉 招，張雙成，楊 耘

（1.長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院，陜西西安710054；2.國(guó)土資源部退化及未利用土地整治工程重點(diǎn)實(shí)驗(yàn)室，陜西西安710075；3.長(zhǎng)安大學(xué)地球科學(xué)與資源學(xué)院，陜西西安710054；4.長(zhǎng)安大學(xué) 環(huán)境科學(xué)與工程學(xué)院，陜西西安710054）

0 引 言

地表發(fā)射率是遙感數(shù)據(jù)獲取陸面溫度中的一個(gè)關(guān)鍵參數(shù)，它與地表組成成分、地表粗糙度、含水量等因素有關(guān)，普遍應(yīng)用于輻射傳輸過(guò)程和地氣系統(tǒng)的能量平衡［1］。熱紅外遙感傳感器所接收到的信號(hào)包括2部分，一部分是被測(cè)物體發(fā)射的熱輻射，另一部分是由于地物的發(fā)射率小于黑體而反射了環(huán)境輻射的冷發(fā)射［2］。因此，為了求得地物的真實(shí)溫度，還必須知道地物的發(fā)射率。一些學(xué)者進(jìn)一步研究發(fā)現(xiàn)，地表發(fā)射率有0.01的變化，估算出的地表溫度就會(huì)產(chǎn)生1～2℃的誤差［3-6］。因此，地表發(fā)射率的準(zhǔn)確獲取，對(duì)地表溫度反演精度的提高至關(guān)重要。

目前，利用RS反演區(qū)域地表發(fā)射率是獲取大范圍地表發(fā)射率的重要途徑［7］。張成才等采用Van和Owe提出的植被指數(shù)法，建立經(jīng)驗(yàn)公式并結(jié)合典型地表發(fā)射率值計(jì)算得到河南省白沙灌區(qū)的地表發(fā)射率［8］。林志斌采用三波段梯度差法簡(jiǎn)單估算了植被覆蓋度，并選取適合Landsat ETM+影像的植被指數(shù)混合模型算法計(jì)算得到了地表發(fā)射率［9］。唐伯惠等提出基于分類的地表發(fā)射率反演算法，由于地表狀況本身的復(fù)雜性以及衛(wèi)星觀測(cè)中的混合像元問(wèn)題，發(fā)射率反演精度通常達(dá)不到所需的0.01［10］?？紤]到NDVI閾值法的可操作性好及較高的精度，趙少華等應(yīng)用Sobrino等提出的NDVITEM算法獲取山西省地表發(fā)射率，并反演得到較高精度的地表溫度值［11］。覃志豪等對(duì)Sobrino等提出的NDVITEM算法進(jìn)行了改進(jìn)，考慮了不同地表形態(tài)下dε（粗糙表面的空腔效應(yīng)引起的發(fā)射率比例）取值以及城市像元發(fā)射率估算的問(wèn)題，并對(duì)山東省陵縣區(qū)進(jìn)行地表發(fā)輻射率反演，取得了較合理的地表溫度反演結(jié)果［12-13］。與之前的Landsat系列衛(wèi)星相比，Landsat8自2013年發(fā)射以來(lái)，在許多方面進(jìn)行了改進(jìn)，也被廣泛應(yīng)用于地表溫度的反演［14-16］。然而，僅針對(duì)Landsat 8數(shù)據(jù)進(jìn)行地表發(fā)射率反演算法的對(duì)比研究較為少見(jiàn)。

歸一化植被指數(shù)（Normalized Difference Vegetation Index，NDVI）是描述地表覆被生長(zhǎng)狀況的一個(gè)重要指標(biāo)［17-18］。由于地表熱紅外波段的發(fā)射率同NDVI之間具有非常高的相關(guān)性，目前基于NDVI獲取地表發(fā)射率的方法主要有4種。即：Van和Owe提出的植被指數(shù)法，Valor和Caselles提出的植被指數(shù)混合模型法，Sobrino等提出的NDVITEM算法，覃志豪等提出的NDVITEM改進(jìn)算法［19-21］。4種算法的表達(dá)式及基本參數(shù)設(shè)定都不相同。因此文中利用西安市的Landsat8數(shù)據(jù)，應(yīng)用ENVI，ArcGIS等軟件，分別用4種算法對(duì)研究區(qū)進(jìn)行了地表發(fā)射率反演，并結(jié)合同時(shí)期精度0.01的MODIS LSE產(chǎn)品，對(duì)比分析了4種算法的適用性和精度結(jié)果，以期獲取更精確的地表發(fā)射率反演算法。

1 數(shù)據(jù)來(lái)源及處理

1. 1 研究數(shù)據(jù)

西安市地處西北部關(guān)中盆地中部秦嶺北麓，地跨渭河南北兩岸，位于北緯 33°39′～34°45′，東經(jīng)107°40′～109°49′.下轄 11區(qū) 2縣并代管西咸新區(qū)，總面積 10 752 km2，2017年末戶籍人口905.68萬(wàn)。屬暖溫帶半濕潤(rùn)季風(fēng)氣候，平均海拔高424 m，1月份平均氣溫0.4℃，7月份平均氣溫26.6℃，年平均氣溫13.3℃，年平均降水量613.7 mm，年平均濕度69.6%［22-23］.研究數(shù)據(jù)采用2015年4月28日西安市的Landsat8數(shù)據(jù)，遙感影像清晰無(wú)云，質(zhì)量良好，便于解譯。基于ENVI5.3平臺(tái)對(duì)遙感影像進(jìn)行的前期處理有：影像鑲嵌和裁剪；OLI數(shù)據(jù)進(jìn)行輻射定標(biāo)；利用FLAASH模型對(duì)OLI數(shù)據(jù)進(jìn)行大氣校正。為了更好地比較不同算法、不同下墊面類型反演地表發(fā)射率的差異情況，文中使用分類回歸決策樹方法將地表分成植被、建筑物、裸土和水體4類。在ENVI 5.3光譜指數(shù)計(jì)算器中分別獲取NDVI（歸一化植被指數(shù)），NDBI（歸一化建筑指數(shù)），MNDWI（改進(jìn)型歸一化水體指數(shù)），結(jié)合人工判定閾值（根據(jù)10%置信度，NDVI＞0. 65為植被；NDBI＞-0. 05為裸土；MNDWI＞0. 17為水體）建立決策樹模型。地表分類結(jié)果如圖1所示。

圖1 西安市決策樹分類結(jié)果Fig.1 Xi’an decision tree classification results

1. 2 對(duì)比數(shù)據(jù)

文中使用了MODISMOD11A1地表發(fā)射率產(chǎn)品進(jìn)行Landsat 8反演的地表發(fā)射率結(jié)果精度驗(yàn)證。MOD11A1是經(jīng)過(guò)正弦投影后的全球每日地表發(fā)射率產(chǎn)品，該產(chǎn)品基于地表下墊面的分類思想，利用地表分類信息，通過(guò)模型計(jì)算實(shí)現(xiàn)了全球地表發(fā)射率產(chǎn)品的生成［9］。在晴朗天氣條件下，MOD11A1地表發(fā)射率產(chǎn)品精度可達(dá)到0.01.

2 地表發(fā)射率反演算法

2. 1 植被指數(shù)法

根據(jù)NDVI和地表發(fā)射率之間有高度的相關(guān)性為理論基礎(chǔ)，將博茨瓦納采集的發(fā)射率信息與NDVI相關(guān)聯(lián)，建立理論模型。模型如式（1）所示

式中，經(jīng)驗(yàn)公式是基于博茨瓦納熱帶草原總結(jié)的，且 NDVI的適用范圍為0. 157～0. 727.

2. 2 植被指數(shù)混合模型法

該方法的主要思想是：假定植被和裸土的地表發(fā)射率已知分別為0. 985和0. 96. 當(dāng)NDVI＜0. 1時(shí)，認(rèn)為是純裸土像元；當(dāng)NDVI＞0. 72時(shí)認(rèn)為是純植被像元；NDVI介于0. 1至0. 72之間時(shí)，認(rèn)為是混合像元，采用式（2）計(jì)算地表發(fā)射率

式中 ε為地表發(fā)射率；εV為純植被像元發(fā)射率；

式中 NDVI為歸一化植被指數(shù)；NDVImin=0. 1，NDVImax=0. 72，且當(dāng) NDVI＞0. 72時(shí)認(rèn)為地表覆被狀況良好，則植被覆蓋度賦值為1，當(dāng)NDVI＜0. 1時(shí)認(rèn)為地表覆被狀況較差，則植被覆蓋度賦值為0. 發(fā)射率比例dε計(jì)算公式如式（4）

式中，＜dε＞為地形幾何形狀因子，根據(jù)不同幾何分布可取均值0. 015.

2. 3 NDVI TEM算法

Sobrino等在前人研究的基礎(chǔ)上，對(duì)植被指數(shù)混合模型法進(jìn)行了改進(jìn)。其主要思想是：當(dāng)NDVI＜0. 2時(shí)，認(rèn)為地表覆被狀況較差，地表發(fā)射率為ε=0.98-0. 042ρ，式中ρ是紅外波段的地物反射率；NDVI＞0. 5時(shí)認(rèn)為是地表覆被狀況良好，地表發(fā)射率為常數(shù) 0. 99；NDVI在 0. 2至0. 5之間時(shí)認(rèn)為是混合區(qū)，通過(guò)式（5）計(jì)算地表發(fā)射率

式中 ε為地表發(fā)射率；εV為純植被像元發(fā)射率；εS純裸土像元發(fā)射率；PV為植被覆蓋度；dε為考慮了粗糙表面的空腔效應(yīng)引起的發(fā)射率比例。植被覆蓋度PV計(jì)算公式（6）εS為純裸土像元發(fā)射率；PV為植被覆蓋度；dε為考慮了粗糙表面的空腔效應(yīng)引起的發(fā)射率比例。植被覆蓋度（PV）計(jì)算公式如式（3）

式中 NDVI為歸一化植被指數(shù)；NDVImin=0. 2，NDVImax=0. 5，且當(dāng) NDVI＞0. 5時(shí)認(rèn)為地表覆被狀況良好，則植被覆蓋度賦值為1，當(dāng)NDVI＜0. 2時(shí)認(rèn)為地表覆被狀況較差，則植被覆蓋度賦值為0.

式中，F(xiàn)為值為0. 55的形態(tài)參數(shù)。

2. 4 NDVI TEM改進(jìn)算法

覃志豪等考慮了不同地表形態(tài)下dε的取值以及城市像元發(fā)射率估算的問(wèn)題。該方法認(rèn)為：對(duì)于水體像元，一般取 εw=0. 995來(lái)估計(jì)；PV≤0.5時(shí)，dε=0. 003 796PV；PV＞0. 5時(shí)，dε=0. 003 796（1-PV）. 發(fā)射率估算公式如下：自然表面像元發(fā)射率計(jì)算如下式

城鎮(zhèn)像元發(fā)射率計(jì)算如下式式中 εV，εS和εm分別為純植被、純裸土和純建筑表面像元的發(fā)射率，分別取 εV=0. 986，εS=0.972，εm=0. 970. RV，RS和 Rm分別為純植被、純裸土和純建筑表面像元的溫度比率，與植被覆蓋度有關(guān)，計(jì)算公式分別如下

3 反演結(jié)果與對(duì)比分析

3. 1 反演結(jié)果

基于以上4種算法，利用ENVI 5.3圖像處理軟件，對(duì)西安市影像進(jìn)行地表分類，然后計(jì)算NDVI及植被覆蓋度，最后反演得到地表發(fā)射率值。4種算法計(jì)算結(jié)果如圖2所示。

圖2 Landsat8 4種算法發(fā)射率反演結(jié)果Fig.2 LSE obtained from Lansat8 data using four LSE retrieval algorithms

3. 2 對(duì)比分析

3.2.1 4種算法反演結(jié)果與MODISLSE產(chǎn)品比較

4種算法不僅在反演模型上有差異，對(duì)PV為0和1時(shí)NDVI的設(shè)定也不同。這樣的差異隨著研究區(qū)地表覆被狀況的不同，可能會(huì)產(chǎn)生不同的影響。因此有必要對(duì)4種算法進(jìn)行比較，進(jìn)而得到適合本研究區(qū)的算法。

Landsat8熱紅外數(shù)據(jù)與MODIS發(fā)射率產(chǎn)品的空間分辨率不同。為了進(jìn)行便精確之間的比較，首先利用 MRT（MODIS Reprojection Tool）工具將MODIS LSE產(chǎn)品地圖投影轉(zhuǎn)換為UTM-WGS84，并進(jìn)行研究區(qū)范圍的裁剪，然后利用ENVI軟件采用最鄰近內(nèi)插法將Landsat8發(fā)射率反演結(jié)果重采樣到1 000m.文中采用最小值、最大值和平均值這些統(tǒng)計(jì)指標(biāo)，對(duì)反演結(jié)果進(jìn)行評(píng)價(jià)，統(tǒng)計(jì)結(jié)果見(jiàn)表1.

表1 Landsat 8各算法反演結(jié)果與MODIS發(fā)射率產(chǎn)品對(duì)比Table 1 Comparation between Landsat 8 retrieval results and MODIS LSE product

從西安區(qū)域的統(tǒng)計(jì)來(lái)看，利用4種算法反演得到的發(fā)射率結(jié)果，相對(duì)于MODIS LSE產(chǎn)品最小值差別分別為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.其中，植被指數(shù)混合模型法和NDVITEM改進(jìn)算法與MODIS LSE產(chǎn)品最小值差別較小，植被指數(shù)法和NDVITEM改進(jìn)算法與MODIS發(fā)射率產(chǎn)品最大值差別較小，植被指數(shù)混合模型法和NDVITEM改進(jìn)算法與MODIS LSE產(chǎn)品平均值差別較小。總體而言，植被指數(shù)混合模型法和NDVITEM改進(jìn)算法與MODIS LSE產(chǎn)品更接近，可信精度較高。

3.2.2 4種算法之間的比較與分析

為了進(jìn)一步對(duì)比4種反演算法在區(qū)域尺度的差異，文中對(duì)利用各種算法反演的2015年4月28日西安市Landsat 8數(shù)據(jù)發(fā)射率結(jié)果進(jìn)行比較分析，圖3是研究區(qū)內(nèi)4種算法反演結(jié)果兩兩之間的差值統(tǒng)計(jì)圖。

圖3 4種方法反演結(jié)果兩兩之間的差值統(tǒng)計(jì)Fig.3 Statistical analysis of the differences between the inversion results of fourmethods

從圖3可以看出，植被指數(shù)混合模型法和NDVITEM改進(jìn)算法反演結(jié)果差異最小，差值范圍控制在0.04之內(nèi)，且大部分像元發(fā)射率差值靠近0值。NDVITEM算法與NDVITEM改進(jìn)算法反演結(jié)果差異最大。主要原因是：兩者在NDVI最大值和最小值的設(shè)定上有所不同，NDVITEM算法中，NDVImin=0. 2，NDVImax=0. 5. NDVITEM改進(jìn)算法中，NDVImin=0.05，NDVImax=0. 7，且 NDVITEM改進(jìn)算法對(duì)水體像元直接賦給0. 995的發(fā)射率值，而NDVITEM算法在計(jì)算NDVI小于0. 2的像元時(shí)，利用紅外波段的反射率計(jì)算其發(fā)射率值。

為揭示不同算法、不同地表覆被類型反演的地表發(fā)射率的差異情況，文中還基于上文決策樹分類所提取的植被，城鎮(zhèn)，水體，裸土4種地表類型，對(duì)不同算法不同地表覆蓋類型的發(fā)射率平均值兩兩之間的差值進(jìn)行分析見(jiàn)表2.

從表2可以看出，在分類尺度上，對(duì)于植被區(qū)域，植被指數(shù)法與NDVITEM算法差異最大達(dá)到了-0. 009 4，植被指數(shù)混合模型法與 NDVITEM改進(jìn)算法差異最小，僅0. 001 2，對(duì)于城鎮(zhèn)區(qū)域（未考慮植被指數(shù)法），植被指數(shù)混合模型法與NDVITEM算法差異最大，達(dá)到了0. 011 7，NDVITEM算法與NDVITEM改進(jìn)算法差異最小，達(dá)到了0. 002 9，對(duì)于水體區(qū)域（未考慮植被指數(shù)法），NDVITEM算法與NDVITEM改進(jìn)算法差異最大，達(dá)到了0. 033 3，植被指數(shù)混合模型法與NDVITEM改進(jìn)算法差異最小，達(dá)到了-0.009 9，對(duì)于裸土區(qū)域（未考慮植被指數(shù)法），NDVITEM算法與NDVITEM改進(jìn)算法差異最大，達(dá)到了0. 012 0，植被指數(shù)混合模型法與NDVITEM算法差異最小，達(dá)到了0.005 9.總體而言，地表類型為植被的情況下，算法之間的差異最小，其次為城鎮(zhèn)和裸土，地表類型為水體的情況下，算法之間的差異最大。

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

4 結(jié) 論

1）從像元尺度上，植被指數(shù)混合模型法與NDVITEM改進(jìn)算法總體精度相對(duì)較高，植被指數(shù)法與NDVITEM算法總體精度次之；

2）從下墊面類型尺度上，對(duì)于植被區(qū)域，4種算法之間的差異最小，其次為城鎮(zhèn)和裸土。對(duì)于水體區(qū)域，4種算法之間的差異最大；

3）從算法的適用性而言，植被指數(shù)法適用性最弱，該模型只能用來(lái)反演NDVI值范圍為0. 157～0. 727的自然表面，且反演精度較低。植被指數(shù)混合模型法與NDVITEM改進(jìn)算法適用性較好，且算法之間的差異最小。這與兩者模型中NDVI最大值和最小值的設(shè)定上有關(guān)，植被指數(shù)混合模型法中，NDVImin=0. 1，NDVImax=0. 72；NDVITEM改進(jìn)算法中，NDVImin=0. 05，NDVImax=0.7.

此外，針對(duì)像元尺度的地表發(fā)射率驗(yàn)證，還需比較遙感計(jì)算和實(shí)地測(cè)量的地表發(fā)射率值，因此在后續(xù)工作中需要進(jìn)一步拓展，以期更好地為提高反演地表發(fā)射率的精度提供線索。