基于哨兵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个，建议在乡镇区域站气象观测要素中可添加地表温度，或在今后的研究中加入人工实地观测，增加地面验证点密度，进一步提高验证的精度。