基于哨兵2/3衛星高分辨率地表溫度的降尺度方法

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

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

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

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

1 研究區概況及研究方法

1.1 研究區概況

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

1.2 研究方法

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

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

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

(1)

其中，TS1 000為1 000 m尺度上根據NDVI與LST的關系，計算得到的LST真實值；a，b為回歸系數；INDV1 000為1 000 m尺度上NDVI值；ε為回歸殘差。式中f(INDV)=a+bINDV，f(INDV)是NDVI和LST之間的函數關系，是NDVI的轉換函數，也適用于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轉換函數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衛星數據的地表溫度空間降尺度方法技術路線

1.3 地表光譜數據及預處理

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

表1 研究中所用衛星數據信息列表

(5)

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

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

TS=-27.584INDV+325.070。

(6)

式中，TS為地表溫度值。

1.4 驗證方法

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

圖2 研究區地面氣象站點分布

2 結果分析

2.1 降尺度結果分析

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

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

2.2 地面點驗證

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

3 結論

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

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

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

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