袁 媛,郭毅春,杜 婧,屈 妍

(1.商洛市氣象局，陜西商洛 726000；2.陜西省氣象局秦嶺和黃土高原生態(tài)環(huán)境氣象重點實驗室，西安 710016)

“哨兵(Sentinel)”系列衛(wèi)星是歐洲“哥白尼”對地觀測計劃部分的專用衛(wèi)星系列。其中，哨兵2衛(wèi)星承擔著多光譜高分辨率成像任務，用于陸地監(jiān)測，可提供植被、土壤和水覆蓋、內(nèi)陸水路及海岸區(qū)域等圖像，還可用于緊急救援服務。哨兵3衛(wèi)星攜帶多種有效載荷，用于高精度測量海面地形、海面和地表溫度、海洋水色和土壤特性，還將支持海洋預報系統(tǒng)及環(huán)境與氣候監(jiān)測。哨兵3-A衛(wèi)星上搭載的海陸表面溫度輻射計(SLSTR)擁有較高時間分辨率和2個熱紅外通道，提供垂直觀測(0°)和前向觀測(55°)，可獲取高精度的海面和地表溫度[1]。

地表溫度(land surface temperature, LST)就是地面的溫度。LST作為一種重要的氣候參數(shù)，在城市熱島效應、火情監(jiān)測等領域被廣泛應用[2-4]。但是，熱遙感系統(tǒng)存在時空分辨率之間的矛盾，不能同時滿足時間和空間的高分辨率需求[5]，比如靜止氣象衛(wèi)星(FY-4)可連續(xù)觀測地表熱特征，但是空間分辨率粗糙，不能很好地反應地物特征[6]，MODIS Terra/Aqua 和哨兵3每天都可以采集數(shù)據(jù)，同樣空間分辨率較低[7]。因此，對低空間分辨率的LST進行降尺度，提高其空間分辨率很有必要[8]。目前，國內(nèi)外大多數(shù)學者對于地表溫度降尺度的研究是基于MODIS地表溫度產(chǎn)品和Landsat 影像[9-10]，鮮少有學者使用哨兵2/3影像數(shù)據(jù)。為填補國內(nèi)空白，本研究系統(tǒng)介紹了運用哨兵2影像對哨兵3地表溫度進行降尺度的方法，為哨兵數(shù)據(jù)用戶的相關研究提供參考。

1 研究區(qū)概況及研究方法

1.1 研究區(qū)概況

研究區(qū)域為商洛市周邊地區(qū)，隨機選取2019年一張成像質量高、天空晴朗無云的影像，將影像中的區(qū)域作為研究區(qū)，該區(qū)橫跨河南省三門峽市、南陽市少部分地域，陜西省渭南市南部及商洛市中東部大部分地區(qū)。該區(qū)域橫跨亞熱帶和暖溫帶兩個過渡性季風氣候帶，四季分明，冬無嚴寒，夏無酷暑。區(qū)域內(nèi)水資源豐富，森林植被覆蓋率分布不均，地形、地貌差異較大，因此，對該區(qū)域進行LST降尺度方法研究有較好的典型性。

1.2 研究方法

降尺度方法可分為熱銳化和溫度分解兩種方法[11-12]。本研究使用的統(tǒng)計降尺度方法是熱銳化方法的一種，它是通過使用多源衛(wèi)星探測器，將低空間分辨率采集的數(shù)據(jù)應用在高空間分辨率的探測器上，進而獲得高時空分辨率影像的過程。本研究是利用歸一化植被指數(shù)(Normalized Differential Vegetation Index, NDVI)和LST之間存在的相關關系，結合哨兵2/3的觀測數(shù)據(jù)實現(xiàn)的。

首先提取哨兵3的1 000 m尺度上LST與NDVI，建立兩者之間的關系，得到簡單的線性回歸模型，計算公式如下

TS1 000=a+bINDV1 000+ε=f(INDV1 000)+ε。

(1)

其中，TS1 000為1 000 m尺度上根據(jù)NDVI與LST的關系，計算得到的LST真實值；a，b為回歸系數(shù)；INDV1 000為1 000 m尺度上NDVI值；ε為回歸殘差。式中f(INDV)=a+bINDV，f(INDV)是NDVI和LST之間的函數(shù)關系，是NDVI的轉換函數(shù)，也適用于10 m尺度上NDVI與LST之間的關系轉換，它的這一特性是能夠將哨兵3的1 000 m空間分辨率上的LST降尺度到哨兵2的10 m空間分辨率的關鍵所在。另外，應用相關關系的回歸模型對低空間分辨率NDVI影像LST值也可進行模擬，T′S1 000表示1 000 m尺度上LST的估算值，如公式(2)所示

T′S1 000=f(INDV1 000)。

(2)

但在這過程中，由于受到不同地形、地貌等因素的影響，運用NDVI很難精準地計算出LST實際值，由公式(1)和(2)可得到在1 000 m空間分辨率上每個像元的殘差值ΔT′S1 000

ΔT′S1 000=TS1 000-T′S1 000=ε。

(3)

將殘差值ΔT′S1 000和在哨兵2的10 m空間分辨率上提取的NDVI值，代入由1 000 m空間分辨率建立的NDVI轉換函數(shù)f(INDV)，得到10 m尺度上的LST值。公式如下所示

T′S10=f(INDV10)+ΔT′S1 000。

(4)

式中，T′S10為10 m尺度上的LST值，它是由10 m空間分辨率的NDVI與1 000 m空間分辨率降至10 m空間分辨率的精確計算殘差相加得到的。具體降尺度方法技術路線如圖1所示。

圖1 基于哨兵2/3衛(wèi)星數(shù)據(jù)的地表溫度空間降尺度方法技術路線

1.3 地表光譜數(shù)據(jù)及預處理

本研究使用的衛(wèi)星數(shù)據(jù)是哨兵2 MSI傳感器和哨兵3 SLSTR傳感器的數(shù)據(jù)(表1) ，通過 ESA 官網(wǎng)(https://scihub.copernicus.eu/dhus/#/home)下載數(shù)據(jù)。哨兵2下載到的數(shù)據(jù)L1C，是大氣表觀反射率產(chǎn)品，已經(jīng)經(jīng)過輻射定標、幾何重采樣、大氣表觀反射率轉換和地理配準等處理；使用歐洲空間局自行研發(fā)的Sen2cor大氣校正插件對數(shù)據(jù)進行批量大氣校正，獲得S2A級數(shù)據(jù)；再使用SNAP7.0遙感影像分析處理軟件進行NDVI計算，如公式(5)所示。

表1 研究中所用衛(wèi)星數(shù)據(jù)信息列表

(5)

式中，INDV為歸一化植被指數(shù)值，ρNIR和ρRED分別表示近紅外波段(NIR)和紅光波段(RED)的反射率，ρB8和ρB4分別表示在S2數(shù)據(jù)中近紅外波段(B8)和紅光波段(B4)的反射率。

哨兵3 SLSTR傳感器數(shù)據(jù)同樣通過ESA官網(wǎng)下載，選取與哨兵2影像時間相近的一幅影像提取研究數(shù)據(jù)，數(shù)據(jù)類別為SL_2_LST。在SNAP7.0軟件中創(chuàng)建數(shù)據(jù)子集，包括波段、空間、連接點網(wǎng)格子集，并對數(shù)據(jù)進行重投影(UTM/WGS 84)，在Excel2010軟件中對提取到的NDVI和LST數(shù)值進行擬合，得到的LST和NDVI的擬合關系式如下

TS=-27.584INDV+325.070。

(6)

式中，TS為地表溫度值。

1.4 驗證方法

為驗證生成的高空間分辨率LST的精度，選取商洛市5個地面國家氣象自動觀測站實測0 cm地溫數(shù)據(jù)進行驗證，將數(shù)據(jù)單位由攝氏溫度℃轉換為絕對溫度K，站點分布如圖2所示。

圖2 研究區(qū)地面氣象站點分布

2 結果分析

2.1 降尺度結果分析

利用ArcMap10.2軟件輸出原始1 000 m空間分辨率哨兵3 LST圖，與降尺度到10 m空間分辨率的LST圖對比，可以發(fā)現(xiàn)，原始1 000 m LST與降尺度結果中地表溫度信息特征基本一致，高溫區(qū)和低溫區(qū)均高度吻合，說明降尺度結果較好地保留了原始LST影像熱特征的分布情況。另外，10 m LST降尺度影像地貌紋理、地物特征更加清晰，色調更加豐富，溫度過渡較平滑，可以清楚地看到地表溫度的空間異質性，且原始1 000 m LST影像的“馬賽克”現(xiàn)象消失，降尺度前后的低溫區(qū)和高溫區(qū)分布也高度吻合(如圖3)。

圖3 哨兵3原始LST影像(a)及降尺度所得10 m LST影像(b)

2.2 地面點驗證

選取本研究區(qū)商洛市境內(nèi)的5個國家氣象觀測站與遙感影像相同時間的0 cm地溫分鐘數(shù)據(jù)，根據(jù)氣象觀測站的精準地理位置信息，在降尺度結果中提取對應的LST，發(fā)現(xiàn)誤差結果平均值為2.6 K，誤差很小，說明降尺度結果精度較高。

3 結論

本研究利用哨兵2詳細的地物空間信息和哨兵3影像的高時間分辨率地表變化信息，以NDVI和LST的相關關系為基礎，運用統(tǒng)計降尺度方法將LST空間尺度從1 000 m降至10 m，生成高時間分辨率10 m空間分辨率的LST。通過國家地面氣象站0 cm地溫分鐘觀測數(shù)據(jù)進行驗證，得出以下主要結論。

(1)本研究方法適用于地表空間異質性區(qū)域，所生成的高空間分辨率的地表溫度產(chǎn)品地物特征明顯，地貌紋理清晰，能夠清楚地描繪地表熱特征空間分布情況。

(2)利用國家地面氣象站0 cm地溫分鐘觀測數(shù)據(jù)對降尺度結果進行驗證，誤差平均值很小，說明降尺度結果精度較高，具有廣泛適用性。

(3)驗證所用到的地面氣象站點只有5個，建議在鄉(xiāng)鎮(zhèn)區(qū)域站氣象觀測要素中可添加地表溫度，或在今后的研究中加入人工實地觀測，增加地面驗證點密度，進一步提高驗證的精度。