中巴經(jīng)濟(jì)走廊2013-2018年地表溫度高分辨率反演數(shù)據(jù)集

趙國(guó)輝，張耀南，康建芳

1. 中國(guó)科學(xué)院西北生態(tài)環(huán)境資源研究院，寒旱區(qū)科學(xué)大數(shù)據(jù)中心，蘭州 730000

2. 中國(guó)科學(xué)院大學(xué)，北京 100049

3. 中國(guó)科學(xué)院超級(jí)計(jì)算蘭州分中心，蘭州 730000

數(shù)據(jù)庫(kù)（集）基本信息簡(jiǎn)介

數(shù)據(jù)庫(kù)（集）組成國(guó)家科技基礎(chǔ)條件平臺(tái)“特殊環(huán)境特殊功能觀測(cè)研究臺(tái)站共享服務(wù)平臺(tái)”項(xiàng)目（Y719H71001）；中國(guó)科學(xué)院信息化專項(xiàng)“寒旱區(qū)環(huán)境演變研究‘科技領(lǐng)域云’的建設(shè)與應(yīng)用”項(xiàng)目（XXH13506）。數(shù)據(jù)集按照年份共分為6個(gè)目錄，在每個(gè)年份目錄下按照月份劃分為12個(gè)子目錄，數(shù)據(jù)命名規(guī)則為 CPEC_LST_Xm_PPPRRR _YYYYMMDD.tif，其中 X為分辨率，PPPRRR為 WRS2行列號(hào)，YYYYMMDD 為年月日標(biāo)識(shí)，例如CPEC_LST_30m_152042_20170603.tif。

引 言

地表溫度（Land Surface Temperature，LST）作為陸地表面熱輻射狀態(tài)的度量，綜合了地表與大氣相互作用以及大氣和陸地之間能量交換的全部結(jié)果，是表征地表水熱平衡和氣候變化的關(guān)鍵參數(shù)、地表循環(huán)過程研究的重要組成部分。準(zhǔn)確地反演區(qū)域地表溫度時(shí)空變化不僅對(duì)氣候變化、生態(tài)保護(hù)以及資源評(píng)價(jià)等相關(guān)領(lǐng)域的研究至關(guān)重要，而且對(duì)環(huán)境監(jiān)測(cè)、災(zāi)害防治、工程建設(shè)等也具有重要應(yīng)用價(jià)值[1]。

傳統(tǒng)地表溫度監(jiān)測(cè)技術(shù)很難獲取其時(shí)空分布特征數(shù)據(jù)，在區(qū)域研究和應(yīng)用中受極大的限制，而熱紅外遙感探測(cè)技術(shù)能夠探測(cè)到地表直接發(fā)射能量，成為獲取區(qū)域地溫的重要途徑[2]?；谶b感手段的陸表參數(shù)反演為地表溫度時(shí)空數(shù)據(jù)制備提供了直接的技術(shù)手段。目前常用的地表溫度遙感數(shù)據(jù)多數(shù)來(lái)自MODIS衛(wèi)星，其雖具有較好的區(qū)域性和實(shí)效性，對(duì)于地表變化連續(xù)觀測(cè)具有明顯優(yōu)勢(shì)，但是無(wú)法滿足局部區(qū)域環(huán)境變化監(jiān)測(cè)的空間分辨率。Landsat系列衛(wèi)星數(shù)據(jù)憑借其長(zhǎng)期連續(xù)、全球覆蓋、高分辨率、科學(xué)的數(shù)據(jù)存檔與分發(fā)策略等優(yōu)勢(shì)，成為地表特征和地球系統(tǒng)科學(xué)研究中最有效的遙感數(shù)據(jù)之一[3]。自2013年Landsat 8衛(wèi)星成功發(fā)射以來(lái)，提供免費(fèi)的中分辨率數(shù)據(jù)以及16天的重放周期，成為區(qū)域熱紅外遙感的重要數(shù)據(jù)源，目前已擁有5年多的熱紅外波段存檔數(shù)據(jù)，為地表溫度反演提供了重要的數(shù)據(jù)支撐。

目前針對(duì)不同的應(yīng)用領(lǐng)域以及區(qū)域，國(guó)內(nèi)外學(xué)者發(fā)展了不同的地表溫度反演模型。為了準(zhǔn)確地反演地表溫度的時(shí)空變化，兼顧模型的適用性，構(gòu)建了地表溫度反演模型集，其中包括單通道反演模型、雙通道劈窗反演模型以及數(shù)據(jù)融合模型。依據(jù)遴選標(biāo)準(zhǔn)，自動(dòng)生成制備了中巴經(jīng)濟(jì)走廊區(qū)域2013-2018年Landsat 8地表溫度遙感反演數(shù)據(jù)集（圖1），以期為中巴經(jīng)濟(jì)走廊自然環(huán)境變時(shí)空化監(jiān)測(cè)研究、工程建設(shè)以及社會(huì)服務(wù)提供基礎(chǔ)數(shù)據(jù)資料。

1 數(shù)據(jù)采集和處理方法

1.1 基礎(chǔ)數(shù)據(jù)

Landsat 8衛(wèi)星攜帶的2個(gè)主要傳感器陸地成像儀（OLI）和熱紅外傳感器（TIRS），在掃描方式、波段設(shè)置、輻射分辨性能等方面較之前的 Landsat系列衛(wèi)星均有很大的改進(jìn)，尤其是熱紅外遙感器TIRS，可探測(cè)10.6-11.2 μm和11.5-12.5 μm兩個(gè)熱紅外波段的輻射，瞬時(shí)視場(chǎng)角15°，地面幅寬185 km，可以更好地對(duì)熱紅外數(shù)據(jù)進(jìn)行大氣校正，進(jìn)而能夠?qū)Φ乇頊囟冗M(jìn)行準(zhǔn)確地反演。為此，本數(shù)據(jù)集基于Landsat 8的兩個(gè)熱紅外通道（第10與11波段），反演制備地表溫度產(chǎn)品?；A(chǔ)數(shù)據(jù)來(lái)源于美國(guó)USGS EarthExplorer（https://earthexplorer.usgs.gov/）官方網(wǎng)站，獲取了自2013年到目前所有的Landsat 8 OLI/TIRS C1 Level-1數(shù)據(jù)產(chǎn)品（表1），具體行列號(hào)見圖1。

圖1 中巴經(jīng)濟(jì)走廊區(qū)域2013-2018年地表溫度遙感反演數(shù)據(jù)集覆蓋范圍

表1 Landsat 8 OLI/TIRS C1 Level-1數(shù)據(jù)產(chǎn)品描述

1.2 處理方法

本數(shù)據(jù)集的數(shù)據(jù)處理流程如圖2，具體步驟如下：

圖2 中巴經(jīng)濟(jì)走廊地表溫度數(shù)據(jù)集制備流程

步驟 1首先對(duì)Landsat 8的基礎(chǔ)數(shù)據(jù)按照數(shù)據(jù)質(zhì)量進(jìn)行篩選，把云量大于20%的數(shù)據(jù)剔除，然后分別對(duì)OLI和TIRS傳感器數(shù)據(jù)進(jìn)行處理。

步驟 2按照Landsat 8數(shù)據(jù)手冊(cè)[4]，分別計(jì)算TIRS第10、11波段的輻射亮度值Lλ，并基于輻射亮度值計(jì)算星上亮度溫度BT。依據(jù)裂窗協(xié)方差比算法[5,10]反演大氣水汽含量ω，以修正大氣和地表比輻射率的影響，算法如下：

式中τi，τj分別為第10波段和第11波段的大氣透過率，εi與εj為其對(duì)應(yīng)的比輻射率，為其移動(dòng)窗口內(nèi)第10波段、第11波段的亮溫均值，k為窗口內(nèi)像元的索引，N為窗口大小，c0、c1與c2為基于大氣輻射傳輸模型 MODTRAN和大氣廓線數(shù)據(jù)庫(kù) TIGR（Thermodynamic Initial Guess Retrieval database）模擬產(chǎn)生的回歸系數(shù)。

步驟 3地表發(fā)射率是衡量表面散熱效率的指標(biāo)，不同的下墊面有不同的發(fā)射率值（0-1）。本數(shù)據(jù)集采用植被覆蓋度加權(quán)法[6-7]，基于Landsat 8 OLI第4、5波段獲取地表發(fā)射率ε：

式中：植被組分發(fā)射率εv與地表背景發(fā)射率εs數(shù)據(jù)來(lái)源于光譜數(shù)據(jù)庫(kù)；表示像元內(nèi)組分間多次散射而形成的腔體效應(yīng)參數(shù)。FVC為植被覆蓋度，NDVIs和NDVIv分別為裸土和濃密植被的歸一化植被指數(shù)NDVI值，ρ4與ρ5為經(jīng)過大氣輻射校正模型COST（Cosine of solar zenith angle）修正過的Landsat 8紅波段Band 4與近紅外波段Band 5的地表反射率。

步驟 4地表熱輻射傳輸方程是遙感反演地表溫度的基礎(chǔ)，依據(jù)地表熱輻射理論，建立地表溫度LST的反演計(jì)算模型集，包括輻射傳輸方程單通道反演模型（SC1）[8]、單通道發(fā)射率校正模型（SC2）[9]、雙通道劈窗算法模型（SW）[10]以及數(shù)據(jù)融合算法（DF），其具體推導(dǎo)過程及參數(shù)不再一一贅述，用戶可具體參考每個(gè)算法對(duì)應(yīng)的文獻(xiàn)資料，主要公式如下：

1) SC1模型

式中：ψ1、ψ2與ψ3是大氣函數(shù)，利用大氣水汽含量ω近似獲取，Lλ為星上輻射亮度值，bi,0、bi,1與bi,2為回歸參數(shù)，K為常數(shù)，其他符號(hào)同上。

2) SC2模型

3) SW模型

式中：BT11為L(zhǎng)andsat 8第11波段亮度溫度，bn（n= 0,1,...7）為模型回歸系數(shù)，ε和?ε分別為TIRS兩個(gè)通道的發(fā)射率均值與差值。

4) DF模型

依據(jù)空間相似性，每個(gè)MOD11A1數(shù)據(jù)產(chǎn)品像元對(duì)應(yīng)不同地表類型Landsat 8像元的地表溫度表達(dá)式為：

式中：PM為MOD11A1像元值，N為一個(gè)MOD11A1像元對(duì)應(yīng)Landsat像元的個(gè)數(shù)，m為N個(gè)像元中所對(duì)應(yīng)的高分辨率地表覆蓋類型總數(shù)，PL,j為m個(gè)地物類別中j類別的像元值，fj標(biāo)識(shí)不同地表類型之間高低分辨率之間的對(duì)應(yīng)關(guān)系，ai,j(i= 0,1,...4) 為MODIS與Landsat混合像元法求解出的回歸參數(shù)。

步驟 5遴選地表溫度產(chǎn)品。以上4種反演模型產(chǎn)品以及第10與第11波段反演得到的亮溫均值產(chǎn)品共計(jì)5種地表溫度產(chǎn)品輸出。為了與MODIS MOD11A1地溫?cái)?shù)據(jù)產(chǎn)品空間網(wǎng)格單元一一對(duì)應(yīng)，每種模型反演產(chǎn)品分別進(jìn)行坐標(biāo)轉(zhuǎn)換與重采樣處理，然后逐像元進(jìn)行統(tǒng)計(jì)比較，依據(jù)評(píng)價(jià)指標(biāo)，選擇最優(yōu)反演結(jié)果納入數(shù)據(jù)集。

2 數(shù)據(jù)樣本描述

地表溫度產(chǎn)品數(shù)據(jù)集命名格式為CPEC_LST_Xm_PPPRRR_YYYYMMDD.tif，其中X表示地表溫度產(chǎn)品的分辨率，Landsat反演的地溫產(chǎn)品為30 m，YYYY表示年份，MM表示月份，DD表示天，PPP是 Landsat WRS2全球參考系統(tǒng)的條帶號(hào)，RRR 為其行號(hào)，例CPEC_LST_30m_152042_20170603.tif，中巴經(jīng)濟(jì)走廊區(qū)域中的數(shù)據(jù)如圖3所示。目前數(shù)據(jù)起止時(shí)間為2013-2018年，并實(shí)時(shí)更新，目前數(shù)據(jù)共17 018景，數(shù)據(jù)容量3 822 GB，數(shù)據(jù)文件按照時(shí)間進(jìn)行分類，數(shù)據(jù)格式為Geotiff文件，坐標(biāo)系統(tǒng)為UTM坐標(biāo)系統(tǒng)，數(shù)據(jù)類型為Single，溫度單位為℃，可使用通過空間分析軟件直接讀取顯示、融合以及分析（圖4）。

圖3 中巴經(jīng)濟(jì)走廊地表溫度產(chǎn)品單幅數(shù)據(jù)示例（CPEC_LST_30m_149038_20170918.tif）

3 數(shù)據(jù)質(zhì)量控制和評(píng)估

3.1 質(zhì)量控制

數(shù)據(jù)集生產(chǎn)過程中對(duì)模型反演結(jié)果進(jìn)行了逐景比較，每種模型反演的數(shù)據(jù)產(chǎn)品進(jìn)行空間融合后與MODIS MOD11A1 V6地表溫度產(chǎn)品逐像元比較，按照評(píng)價(jià)指標(biāo)從中選擇最優(yōu)反演產(chǎn)品。本數(shù)據(jù)集的數(shù)據(jù)遴選指標(biāo)包括納什效率系數(shù)（NSE）、一致性指數(shù)（d）、Kling-Gupta系數(shù)（KGE）[11]、相關(guān)系數(shù)（R）與百分比偏差（Pbias）5種指標(biāo)對(duì)數(shù)據(jù)產(chǎn)品進(jìn)行質(zhì)量控制。

3.2 數(shù)據(jù)評(píng)估

數(shù)據(jù)集質(zhì)量評(píng)估采用整體評(píng)估和局部評(píng)估相結(jié)合的方法。整體評(píng)估采用隨機(jī)抽樣的方法，從全部數(shù)據(jù)產(chǎn)品中隨機(jī)選取6組數(shù)據(jù)進(jìn)行對(duì)比實(shí)驗(yàn)，數(shù)據(jù)樣本分布如圖5所示，每組數(shù)據(jù)產(chǎn)品的對(duì)比評(píng)價(jià)結(jié)果見表2。

圖4 中巴經(jīng)濟(jì)走廊全區(qū)域地表溫度產(chǎn)品數(shù)據(jù)拼接融合示例（CPEC_LST_1km_201807.tif）

圖5 中巴經(jīng)濟(jì)走廊地表溫度數(shù)據(jù)集隨機(jī)抽樣數(shù)據(jù)空間分布圖（紅色表示）

表2 中巴經(jīng)濟(jì)走廊地表溫度數(shù)據(jù)集抽樣數(shù)據(jù)整體統(tǒng)計(jì)特征

其次，為確保數(shù)據(jù)統(tǒng)計(jì)的隨機(jī)性，開展了局部抽樣評(píng)估，即再次從每組數(shù)據(jù)產(chǎn)品中隨機(jī)選取子區(qū)域，進(jìn)行對(duì)比分析，驗(yàn)證其空間相關(guān)性（圖6），結(jié)果表明R2均在0.66以上，與整體趨勢(shì)吻合。

圖6 中巴經(jīng)濟(jì)走廊地表溫度案例數(shù)據(jù)子區(qū)域隨機(jī)采樣統(tǒng)計(jì)分析

4 數(shù)據(jù)價(jià)值

本數(shù)據(jù)集在整體上與MODIS地表溫度產(chǎn)品保持較高的空間相關(guān)性，由于Landsat 8較高的空間分辨率，可以反演出更加清晰的地表紋理特征，例如水體、居住區(qū)、植被等地物，表達(dá)了更加詳盡的地表溫度細(xì)節(jié)與特征（圖7），可直接用于表征中巴經(jīng)濟(jì)走廊區(qū)域陸表溫度時(shí)空變化，開展探索其時(shí)空分異特征與規(guī)律研究，為研究?jī)鋈跒?zāi)害、氣候變化、災(zāi)害預(yù)測(cè)、熱島效應(yīng)、工程安全等科學(xué)問題提供基礎(chǔ)數(shù)據(jù)資料，服務(wù)于中巴兩國(guó)的線路建設(shè)、生態(tài)保護(hù)、資源調(diào)查、工程建設(shè)以及防災(zāi)減災(zāi)等社會(huì)經(jīng)濟(jì)發(fā)展需求。

圖7 數(shù)據(jù)產(chǎn)品局部對(duì)照?qǐng)D（a為區(qū)域遙感影像，b為MODIS地表溫度產(chǎn)品，c為本數(shù)據(jù)集產(chǎn)品）

5 數(shù)據(jù)使用方法和建議

本數(shù)據(jù)集可通過國(guó)家特殊環(huán)境、特殊功能觀測(cè)研究臺(tái)站共享服務(wù)平臺(tái)（http://www.crensed.ac.cn）申請(qǐng)下載。由于數(shù)據(jù)量比較大也可聯(lián)系數(shù)據(jù)作者，訪問中國(guó)科學(xué)院超級(jí)計(jì)算蘭州分中心的計(jì)算與數(shù)據(jù)平臺(tái)，進(jìn)行在線數(shù)據(jù)下載、共享、分析、處理及應(yīng)用。數(shù)據(jù)文件均為GeoTIFF格式，可通過GIS相關(guān)的軟件如ENVI、GRASS、ArcGIS等，直接進(jìn)行查看與應(yīng)用。

由于地表溫度的遙感反演受地表發(fā)射率、大氣透過率以及下墊面狀況的影響，目前針對(duì)地表溫度的遙感反演產(chǎn)品適宜于基于相對(duì)溫度的科學(xué)應(yīng)用，地表溫度的真實(shí)值需要依據(jù)每個(gè)像元的發(fā)射率、地物分類以及大氣條狀況對(duì)本數(shù)據(jù)產(chǎn)品進(jìn)行修正。

致 謝

感謝國(guó)家特殊環(huán)境、特殊功能觀測(cè)研究臺(tái)站共享服務(wù)平臺(tái)給與的項(xiàng)目支持。感謝 USGS EarthExplorer網(wǎng)站提供免費(fèi)的Landsat 8與MODIS數(shù)據(jù)產(chǎn)品。