Landsat8熱紅外波段地表比輻射率估算方法改進(jìn)

何 好， 張 運,3， 湯 志， 黃 靜， 丁 玥

(1.安徽師范大學(xué) 地理與旅游學(xué)院，安徽 蕪湖 241000；2.資源環(huán)境與地理信息工程安徽工程技術(shù)研究中心，安徽 蕪湖 241000；3.自然災(zāi)害過程與防控安徽省省級重點實驗室，安徽 蕪湖 241000)

地表比輻射率是遙感方法反演地表溫度的重要參數(shù)，作為地物本身內(nèi)在的一種物理屬性，受到地表地物類型，土壤含水量，植被覆蓋度與結(jié)構(gòu)等因素的影響。在地表溫度反演中，地表比輻射率0.01的誤差會導(dǎo)致反演得到的地表溫度誤差達(dá)到1℃以上[1-2]。因此，對地表比輻射率進(jìn)行準(zhǔn)確計算對地表溫度反演尤其重要。目前，遙感數(shù)據(jù)反演地表比輻射率方法主要分為兩大類：第一類是直接反演，運用多個熱紅外波段數(shù)據(jù)計算地表比輻射率，例如Gillespie等提出的ASTER溫度地表比輻射率分離算法[3]、Becker與Li提出的光譜指數(shù)方法[4]以及Watson提出的雙溫法[5]等；第二類是利用歸一化植被指數(shù)(NDVI)間接估算地表比輻射率，例如Van和Owe提出的植被指數(shù)法[6]、Valor和Caselles提出的基于植被指數(shù)混合模型的方法[7]、Sabrino等提出的NDVI閾值法[8]與覃志豪等提出的NDVI閾值改進(jìn)方法[9]。第一類算法需要在8～12μm之間至少3個以上的熱紅外波段[10]，目前大多數(shù)星載熱紅外傳感器都不滿足該條件。第二類算法在目前運用相對較為廣泛，其中，Van和Owe植被指數(shù)法中的經(jīng)驗公式只考慮了自然地表，但在實際應(yīng)用中，由于熱紅外波段必然存在純裸土與純植被等純像元，故Van和Owe的植被指數(shù)法很難滿足目前反演地表比輻射率的精度要求。而Valor和Caselle提出的基于植被指數(shù)混合模型的方法，考慮了純植被像元與純裸土像元，使地表比輻射率的估算使用范圍更加寬廣，但是將水體歸入純裸土像元，導(dǎo)致所得到的水體地表比輻射率數(shù)值低于正常值。Sabrino等提出的NDVI閾值法除了沒有單獨考慮水體地表比輻射率之外，其在純植被像元與純裸土像元的劃分時所設(shè)置的NDVI閾值隨著遙感影像空間分辨率的提高，不完全適用于Landsat8 熱紅外波段。覃志豪等提出的NDVI閾值改進(jìn)方法，對NDVI閾值范圍更加精確，但是其將水體提取出來后直接進(jìn)行賦值0.995，導(dǎo)致所得到的水體地表比輻射率數(shù)值上比實際的偏大，且其提出的公式系數(shù)適用于TM6波段，該公式是否可以適用于Landsat8 熱紅外波段，并沒有得到驗證。目前所提出的地表比輻射率的估算大多還是停留在針對TM/ETM波段得到的公式，其中利用的經(jīng)驗系數(shù)與參數(shù)也是針對TM/ETM數(shù)據(jù)。隨著Landsat8衛(wèi)星2013年的升空，由于其具有高分辨率、具有兩個熱紅外波段等特點，越來越多的研究者們，利用Landsat8數(shù)據(jù)作為自己的數(shù)據(jù)源進(jìn)行地表溫度的反演。為了便于Landsat8 熱紅外波段提取地表溫度，有必要對Landsat8 熱紅外波段的地表比輻射率估算進(jìn)行更深入的探討。

本文首先介紹植被指數(shù)混合模型，在此基礎(chǔ)上結(jié)合ASTER波譜庫數(shù)據(jù)提出新的估算方法得到Landsat8熱紅外波段的地表比輻射率，從簡單實用的角度考慮，不考慮溫度與土壤水分對地表比輻射率的影響。運用MODIS LSE產(chǎn)品、LST產(chǎn)品與地面氣象站點實測氣溫數(shù)據(jù)對反演結(jié)果進(jìn)行驗證，證明本文的計算方法得到的地表比輻射率具有較高精度，可以用于Landsat8熱紅外波段地表比輻射率的估算,也可以用于地表溫度反演研究。

1 方法介紹

1.1 植被指數(shù)混合模型

在植被指數(shù)法的基礎(chǔ)上，Valor和Caselle將表達(dá)式由自然像元擴展到地表像元[11]，通過裸土與植被的地表比輻射率以及植被與土壤組合比例，運用NDVI計算出植被覆蓋度，之后估算得到地表比輻射率。該方法根據(jù)不同地表類型NDVI值不同的特點將地表類型分為裸土、自然地表與植被三種。采用公式(1)計算自然地表像元地表比福射率:

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

(1)

式中，ε為地表比輻射率，εv為純植被像元比輻射率，εs為裸土像元比輻射率，Pv為植被覆蓋度，dε為地形幾何形狀。dε主要與植被覆蓋度、植被距地高度以及植被間距離有關(guān)，在地形相對平坦時，可以忽略不計，在地形起伏較大時，dε采用簡化公式[12]計算：

dε=4PV(1-PV)

(2)

式中，為地形幾何均值，可以根據(jù)不同種植被距地高度與它們之間的距離合理取值求得。

1.2 利用ASTER波譜庫對模型進(jìn)行改進(jìn)

地表比輻射率主要取決于地表物質(zhì)結(jié)構(gòu)和遙感器的熱紅外波段區(qū)間[13}，會隨著熱紅外波段波長的變化而變化[14]，故估算Landsat8地表比輻射率時，必須將Landsat8的第10波段與第11波段地表比輻射率運用相關(guān)參數(shù)分別進(jìn)行計算。從Landsat8遙感影像上觀察，可以將地表結(jié)構(gòu)分為5類，分別是水體、純裸土、城鎮(zhèn)、自然表面和純植被。水面像元、純裸土像元與純植被像元統(tǒng)稱為純像元，城鎮(zhèn)像元與自然表面稱為混合像元，其中城鎮(zhèn)像元可以看成由建筑物表面與植被組成，自然表面主要是由土壤、植被等組成。對于地表比輻射率的計算，自然表面與城鎮(zhèn)像元是應(yīng)該關(guān)注的重點。

1.2.1 典型地物比輻射率 由于Landsat8 TIRS兩個熱波段的波長范圍不同，在計算Landsat8 TIRS熱紅外波段地表比輻射率過程中，需要分別利用ASTER波譜庫中不同波長范圍、不同地表類型與Landsat8 TIRS熱紅外波段的波譜響應(yīng)函數(shù)得到典型地表比輻射率[15-17]。Landsat8 TIRS典型地表比輻射率數(shù)值如表1。水體的比輻射率值：第10波段水體比輻射率ε10w=0.99383，第11波段水體比輻射率ε11w=0.99254。植被的比輻射率相對來說較高，一般認(rèn)為在0.98～0.99之間[15-17]，利用ASTER波譜庫中的植被比輻射率波譜，結(jié)合ASTER光譜響應(yīng)函數(shù)求取針葉林、闊葉林和草地的比輻射率，求其三者平均值作為其植被比輻射率，得到第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 純像元地表比輻射率 在地表比輻射率計算時，直接對純像元地表比輻射率進(jìn)行賦值計算，其中水體比輻射率分別為ε10w=0.99383、ε11w=0.99254；裸土比輻射率分別為ε10s=0.96767、ε11s=0.97790；植被比輻射率分別為ε10m=0.98672、ε11m=0.98990。

1.2.3 自然表面比輻射率的估計方法 基于植被指數(shù)混合模型，不考慮地形起伏與植被結(jié)構(gòu)的影響，假設(shè)每一個Landsat8 TIRS波段的自然表面都可以分為植被與裸土兩種地物類型的比例組合，則自然表面像元的地表比輻射率可表示為

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

(3)

式中，εv，εs分別為植被比輻射率與裸土比輻射率，Pv為Landsat8 TIRS每個自然表面像元的植被覆蓋度。結(jié)合表1數(shù)據(jù)，可以得到兩個熱紅外波段自然地表比輻射率公式分別為

ε10=0.01905Pv+0.96767

(4)

ε11=0.012Pv+0.97790

(5)

1.2.4 城鎮(zhèn)地表比輻射率的估計方法 實際上，Landsat8熱紅外波段中的城鎮(zhèn)像元可以分成建筑物表面與植被的比例組合，基于植被指數(shù)混合模型，城鎮(zhèn)像元地表比輻射率可以表示為

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

(6)

式中，εv，εm分別為植被比輻射率與建筑物表面比輻射率，Pv為Landsat8 TIRS每個城鎮(zhèn)像元的植被覆蓋度。結(jié)合表1數(shù)據(jù)可以得到城鎮(zhèn)像元第10與第11波段地表比輻射率分別為

ε10=0.021835Pv+0.964885

(7)

ε11=0.014785Pv+0.975115

(8)

1.3 地物類型分類方法

首先利用Landsat8陸地成像儀中紅光波段與近紅外波段計算得到影像NDVI值。根據(jù)不同NDVI值所對應(yīng)地物類型不盡相同對地物進(jìn)行分類，由于水體、云與雪的NDVI數(shù)值為負(fù)，可以根據(jù)其特點將水體提取出來。對于純植被、純裸土、自然表面與城鎮(zhèn)像元的提取則可以通過植被覆蓋度進(jìn)行計算。運用NDVI計算植被覆蓋度，植被覆蓋度公式如下：

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

(9)

式中，NDVIv，NDVIs分別是研究區(qū)植被與裸土的NDVI值。若是遙感影像位置位于山區(qū)或者是影像范圍內(nèi)有明顯的植被茂密區(qū)，則取該植被區(qū)域的平均NDVI作為NDVIv值；同樣的，若有明顯的裸露土壤區(qū)域，則取該裸露土壤區(qū)域的平均NDVI值作為NDVIs值。當(dāng)NDVI大于NDVIv時，表明地表類型為純植被；當(dāng)NDVI小于NDVIs且大于等于0時，表明地表類型為裸土；當(dāng)NDVI介于NDVIs與NDVIv兩者中間時，表明地表類型為城鎮(zhèn)或自然表面。通過公式(4)、(5)與公式(7)、(8)分別計算自然像元與城鎮(zhèn)像元的地表比輻射率。在大多數(shù)情況下，植被茂密區(qū)的NDVI數(shù)值都在0.7以上，有時達(dá)到0.8甚至更高，裸土的NDVI值的范圍一般為0.03～0.08。因此，雖然不同研究區(qū)的植被與裸土都有各自的NDVI，NDVIv與NDVIs也有一定大小的差異。但若是研究區(qū)上沒有明顯的植被茂密區(qū)與裸土地區(qū)，一般可以用NDVIv=0.7與NDVIs=0.05來進(jìn)行植被覆蓋度的計算[18]。

圖1 研究區(qū)2018年4月10日Landsat8影像

2 研究區(qū)與數(shù)據(jù)介紹

2.1 研究區(qū)域

為了對上述提出的計算方法進(jìn)行驗證，選擇包含上述5種地表類型區(qū)域為研究區(qū)。研究區(qū)選取安徽省合肥市，合肥市地處中緯度地帶，屬亞熱帶季風(fēng)性濕潤氣候，年降水量適中，地勢相對平坦，平均海拔約為25m，該研究區(qū)主要地表類型為城鎮(zhèn)，植被，水體，裸土，自然地表。Landsat8影像過境GMT時間為2018年4月10日02時42分。

2.2 數(shù)據(jù)來源

2.2.1 Landsat8 數(shù)據(jù) Landsat8 數(shù)據(jù)來自美國地質(zhì)勘探局(USGS)，過境日期為2018年4月10日，云量為0.03%。其中，Landsat8 OLI傳感器包含多光譜波段，空間分辨率為30m(其中全色波段空間分辨率為15m)，用于計算NDVI與植被覆蓋度，進(jìn)而計算Landsat8熱紅外波段地表比輻射率。Landsat8 TIRS熱紅外波段數(shù)據(jù)，空間分辨率為100m，用來反演地表溫度。

2.2.2 MODIS數(shù)據(jù) MOD11A1數(shù)據(jù)來自于NASA官網(wǎng)(http://earthdata.nasa.gov/)，空間分辨率1000m,包含有地表溫度數(shù)據(jù)MODIS LST與MODIS31，32波段地表輻射率數(shù)據(jù)MOD11A1_LSE，過境日期為2018年4月10日，GMT時間為02時35分，與上述的Landsat8影像過境時間基本吻合，值得注意的是，MODIS LST數(shù)據(jù)是按照16位圖像進(jìn)行顯示的，在進(jìn)行地表溫度數(shù)據(jù)驗證時需要將DN值乘0.02得到絕對溫度。而MOD11A1_LSE數(shù)據(jù)按照8位圖像進(jìn)行顯示，在計算第31,32波段地表比輻射率需將DN值除以256以得到地表比輻射率，并以此來驗證本文計算得到的地表比輻射率。

2.2.3 氣象站點氣溫數(shù)據(jù) 運用本文提出的方法計算出研究區(qū)地表比輻射率后，再反演得到地表溫度數(shù)據(jù)，在反演出地表溫度后，利用安徽省合肥市6個氣象站點氣溫數(shù)據(jù)(來自中國氣象數(shù)據(jù)網(wǎng))進(jìn)行地表溫度驗證，以此來驗證本文所采用的方法得到的地表比輻射率精度。

3 研究結(jié)果與分析

為了評估得到的地表比輻射率的精度，采用兩種方法對計算結(jié)果進(jìn)行驗證。一是將本文估算得到的LSE與MODIS LSE產(chǎn)品數(shù)據(jù)進(jìn)行驗證；二是利用Landsat8數(shù)據(jù)運用本文得到的地表比輻射率數(shù)據(jù)反演地表溫度，然后運用地面氣象站點實際氣溫數(shù)據(jù)對反演得到的地面溫度數(shù)據(jù)進(jìn)行精度驗證。

3.1 與MODIS LSE產(chǎn)品的對比驗證

因為Landsat8熱紅外的兩個波段光譜范圍分別與MODIS的31、32波段光譜范圍非常接近，具體參數(shù)見表2，地表比輻射率又與光譜范圍有關(guān)，而MODIS LSE產(chǎn)品數(shù)據(jù)精度較高，故運用MODIS LSE產(chǎn)品對文中計算出的地表比輻射率結(jié)果進(jìn)行對比驗證。首先Landsat8 熱紅外波段空間分辨率與MODIS LSE產(chǎn)品的空間分辨率不同，為了進(jìn)行更加精確的比較，我們將反演所得到的Landsat8的地表比輻射率數(shù)據(jù)重采樣到與MODIS LSE產(chǎn)品的空間分辨率相同，圖2分別為2018年4月10日Landsat8 TIRS第10波段與第11波段利用本文的方法估算出的地表比輻射率影像。從像元整體觀察可以發(fā)現(xiàn)，Landsat8第10波段與第11波段地表比輻射率基本趨勢相同。之后運用結(jié)果中最小值、最大值、平均值與標(biāo)準(zhǔn)差進(jìn)行精度評估，結(jié)果如表3。通過表3可以發(fā)現(xiàn)，Landsat8熱紅外波段的LSE與MODIS LSE產(chǎn)品的標(biāo)準(zhǔn)差、平均值、最大值上相差較小，在最小值上Landsat8 TIRS 11波段與MODIS 32波段地表比輻射率相差較多，最小值相差達(dá)到0.02，在地表比輻射率平均值上第11波段大于第10波段，通過表1可以基本確定本文計算得到的地表輻射率是符合ASTER光譜庫數(shù)據(jù)的。

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

表2 Landsat8與MODIS熱紅外波段參數(shù)

為了進(jìn)一步驗證文中計算得到的地表比輻射率的精度，采用平均誤差Bias、均方根誤差Rmse與相關(guān)系數(shù)R來評價本文計算得出的地表比輻射率與MODIS LSE地表比輻射率之間的差異，分為Landsat8 第10波段與MODIS第31波段、Landsat8 第11波段與MODIS第32波段兩組進(jìn)行統(tǒng)計對比分析。結(jié)果如表3所示。從表3結(jié)果中可以得出，第10波段與第11波段同對應(yīng)MODIS波段的平均誤差均小于0.004，均方根誤差均小于0.005，相關(guān)系數(shù)均大于0.8，為顯著相關(guān)。對比結(jié)果表明本文方法計算的地表比輻射率與MODIS LSE產(chǎn)品具有較好的一致性。第10波段光譜范圍為10.6μm～11.19μm,第11波段光譜范圍為11.5μm～12.51μm，而地表比輻射率在10～12um區(qū)間內(nèi)較為平穩(wěn)，在此區(qū)間以外則變化幅度較大，故第11波段數(shù)據(jù)相關(guān)參數(shù)稍大于第10波段。

表3 Landsat8與MODIS熱紅外波段參數(shù)對比

單獨從以上指標(biāo)并不能完全說明本文所反演的地表比輻射率精度，為了驗證本文方法的穩(wěn)定性，利用計算出的Landsat8熱紅外波段地表比輻射率重采樣后與MODIS LSE產(chǎn)品進(jìn)行差值運算，并把差值分為三種類型：0～0.005，精度較高；0.005～0.01，精度一般；>0.01，精度較低。差值運算如圖3。

(a)第10波段與MODIS 31差值圖 (b)第11波段與MODIS 32差值圖

圖中灰色區(qū)域是反演精度相對較高的地區(qū)，淺綠色是反演精度一般的地區(qū)，紅色是反演精度較低的地區(qū)。從圖3觀察，Landsat8兩熱紅外波段反演精度相對較高，第10波段最大偏差為0.013，第11波段反演最大偏差為0.026。為了更加直觀觀察反演精度，我們通過計算像元所占比例來判斷反演精度，結(jié)果如表4。

可以發(fā)現(xiàn)，Landsat8熱紅外波段利用本文方法計算得到的地表比輻射率精度較高區(qū)域在80%以上，一般區(qū)域第10波段與第11精度分表為12%與15%，較差區(qū)域所占面積分別為4%、5%?？傮w來說本文估算得到的地表比輻射率精度較高，可以滿足目前研究要求。由表3，表4可知，本文估算地表比輻射率方法與MODIS LSE產(chǎn)品具有較高的一致性，精度較高，可適用于Landsat8熱紅外波段地表比輻射率的估算。

表4 差值結(jié)果比例

3.2 基于地表溫度的驗證

為了進(jìn)一步驗證文中估算地表比輻射率方法的準(zhǔn)確性，運用本文估算出的地表比輻射率反演得到研究區(qū)的地表溫度，并利用地面氣象站點氣溫數(shù)據(jù)與MODIS地表溫度產(chǎn)品(MOD11A1)對反演得到的地表溫度數(shù)據(jù)進(jìn)行對比驗證。目前利用Landsat8 熱紅外數(shù)據(jù)反演地表溫度的方法有輻射傳輸方程法[19](RTE)、普適性單通道算法(SC)、基于影像算法(IB)、覃志豪單窗算法(MW)與劈窗算法(SW)。由于因為SC、IB,MW與SW均需要運用相對較多的大氣參數(shù)，為了防止其他參數(shù)對本文計算結(jié)果產(chǎn)生影響，故選用參數(shù)較少的輻射傳輸方程法進(jìn)行地表溫度反演。第10波段與第11波段反演結(jié)果經(jīng)重采樣后結(jié)果如圖4。

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

利用地表比輻射率反演得到的地表溫度采用最小值、最大值與平均值對反演結(jié)果進(jìn)行評價，結(jié)果如表5，可以知道Landsat8第10波段與第11波段反演結(jié)果相對于MODIS LST產(chǎn)品最小值分別差0.9℃、0.02℃，最大值差別分別是-2.85℃、-1.35℃，平均值差別分別為-0.52℃、-0.58℃。從統(tǒng)計指標(biāo)上來看，Landsat8 熱紅外波段利用本文計算得到的地表比輻射率反演的地表溫度可信度相對較高。

表5 Landsat8地表溫度與MODIS LST產(chǎn)品對比統(tǒng)計表

驗證文中由地表比輻射率計算出的地表溫度精度不可僅僅從統(tǒng)計指標(biāo)進(jìn)行驗證，需進(jìn)一步運用實測溫度數(shù)據(jù)進(jìn)行驗證，利用合肥市6個地面氣象站點氣溫數(shù)據(jù)對反演得到的地表溫度進(jìn)行進(jìn)一步的驗證。結(jié)果表明，MODIS LST產(chǎn)品、反演得到的Landsat8 第10波段和第11波段與地面氣象站點相比的絕對誤差分別為1.35℃、0.82℃與0.84℃。其中MODIS溫度產(chǎn)品與氣溫數(shù)據(jù)絕對誤差相對較大，而Landsat8熱紅外波段反演出的地表溫度同氣溫數(shù)據(jù)絕對誤差相對較小，這是由于MODIS溫度產(chǎn)品為北京時間10時35分過境而Landsat8在北京時間10時52分過境，且當(dāng)日合肥市無降水，天氣晴朗，可知，MODIS LST產(chǎn)品絕對誤差較大是符合當(dāng)日氣候特點的。Landsat8衛(wèi)星過境日期為2018年4月10日，由地-氣溫差規(guī)律[20]可知，地表溫度需大于氣象站點實測氣溫數(shù)據(jù)，而MODIS LST產(chǎn)品與反演得到的地表溫度數(shù)據(jù)平均值均大于氣象站點氣溫平均值，這一點也可以證明地表溫度適用性較好。

通過表5、表6對比驗證，可以說明利用本文方法計算的地表比輻射率反演得到的地表溫度精度與適用性都較好，也可以間接驗證該方法求出的地表比輻射率結(jié)果較為準(zhǔn)確，本文提出的地表比輻射率計算方法可用于Landsat8熱紅外波段估算地表比輻射率，進(jìn)而可用于輸入?yún)?shù)對Landsat8地表溫度進(jìn)行反演。

表6 氣象站點數(shù)據(jù)、反演得到的地表溫度與MODIS溫度產(chǎn)品數(shù)據(jù)對比

4 結(jié)語

針對Landsat8 熱紅外波段地表比輻射率估算方法較復(fù)雜的問題，本文基于植被指數(shù)混合模型，提出一種適合Landsat8熱紅外波段并結(jié)合ASTER光譜庫數(shù)據(jù)計算地表比輻射率的方法。通過MODIS 地表比LSE產(chǎn)品、MODISLST產(chǎn)品與地面氣象站實測氣溫數(shù)據(jù)對本文計算出的地表比輻射率進(jìn)行對比驗證后得出該方法精度較高，可以滿足目前對Landsat8 熱紅外波段估算地表比輻射率的精度要求，利用本文計算得到的地表比輻射率反演出的地表溫度也可以滿足地表溫度反演的精度要求。在今后對Landsat8進(jìn)行地表比輻射率與地表溫度反演時，可以更加簡單精確得到實驗結(jié)果。但只對合肥市進(jìn)行了驗證與分析，在其他地方還需要進(jìn)行更廣泛的對比驗證，此外研究區(qū)地形起伏相對較小,在此后的研究中還需要進(jìn)一步考慮地形起伏等因素對地表比輻射率的影響。