卫星遥感地表温度降尺度的光谱归一化指数法

地表温度是地-气系统相互作用过程的重要物理参数，在地表蒸散量估算[1]、土壤水分估计[2]及城市热岛[3]等研究领域中应用广泛。但是，受到成像条件的制约，目前常用的卫星热红外传感器均存在时空分辨率之间的矛盾[4]，或时间分辨率较高而空间分辨率较低 (如MODIS)，或时间分辨率较低而空间分辨率较高 (如Landsat TM)。这种矛盾使得LST数据的使用受到极大限制[5]，为获取满足应用需求的LST产品，地表温度降尺度技术应运而生。

目前LST降尺度方法有很多，基于像元尺度的LST尺度转换方法主要分为3类：数理统计方法、光谱混合模型法和调制分配方法。数理统计回归法的基础假设是“关系尺度不变性”，即LST与回归核的统计关系在各个尺度上保持不变。回归核包括反照率、组分权重和光谱指数等[1, 5-6]。光谱混合模型的基本原理是线性光谱混合模型，通过关联不同尺度上的LST，进而回归求解低尺度空间 (即高分辨) 上的LST[7]。基于调制分配的方法是将高尺度空间的LST按照一定的权重分配给子像元，分配因子包括发射率、同一时刻获取的更高分辨率传感器的LST以及二者的组合等[5]。按照LST降尺度方法中涉及的尺度因子个数，降尺度方法分为单一尺度因子法[1, 6, 8]和多尺度因子法[9-12]。多尺度因子法主要建立LST与多个相关因子间的线性关系，如逐步回归法[9]和主成分分析[10]，或更为复杂的关系，如支持向量机[11]或神经网络[12]。单一尺度因子法因其物理意义明确并且模型简单而得到广泛应用。文献[1]首次提出利用LST与NDVI (normalized difference vegetation index) 的关系来实现LST降尺度的DisTrad (disaggregation procedure for radiometric surface temperature) 算法，通过建立两者间的二次回归模型，成功实现对千米级到百米级LST的降尺度。文献[6]对DisTrad算法进行了改进，通过对LST、植被覆盖度和植被指数进行回归分析的基础上提出了TsHARP (an algorithm for sharpening thermal imagery) 算法，该算法认为LST与植被覆盖度之间的一元线性关系能够达到更好的降尺度效果。文献[8]利用TsHARP法的5种变形讨论了全局、分段回归、局域和土地利用分层的降尺度模型，并认为分层回归降尺度法的精度要高于原始全局方法。但是该方法中涉及的尺度因子均为NDVI，对于下垫面复杂的区域，不同类型地物的LST最佳指示因子并不一定是NDVI指数。

上述研究已经取得突破性的进展，但仍存在一些问题需要解决：①多数研究是先将高分辨率的LST数据聚合至低分辨率而后进行降尺度处理，并再利用原始LST对降尺度结果进行验证，原始LST数据在重采样或升尺度过程中会产生新的误差，会对降尺度结果产生大的干扰，因此其结果缺乏说服力。②NDVI虽然被人们最早用于LST降尺度研究，但是对于下垫面复杂的城市地区，不能较好地反映建筑物和裸地的高温区以及水体造成的低温区。

针对上述问题，本文以TsHARP降尺度方法为基础，直接对原始MODIS LST产品进行降尺度研究，根据地表覆盖类型的不同，选择与LST相关性更好的光谱指数建立二者之间的关系实现降尺度来弥补问题2的不足，由于Landsat 8热红外波段的分辨率明显高于MODIS数据，并且文献[13—14]提出的单通道算法进行Landsat 8热红外传感器TIRS (thermal infrared sensor) 影像地表温度反演具有很好的精度和敏感性，TIRS和MODIS成像时间也非常接近 (本文选择的数据二者相差小于35 min)，因此可以把TIRS反演的地表温度作为当日地表温度的参考值，对不同方法的降尺度结果进行验证。

研究区包含北京北以及河北部分地区，如图 1所示，其地理位置为116°6′15″E—117°45′12″E，40°14′52″N—40°59′19″N，属于暖温带大陆性气候；区内地势平坦，只有西北部和东部有部分山脉，东南部为平原，属于华北平原的西北边缘区。中下部为北京郊区，包括怀柔城区等下垫面为不透水层覆盖地区以及水系比较发达的密云区；周边主要以林地和耕地覆盖，整个研究区地表覆盖类型复杂，空间异质性高，对进行地表温度降尺度研究具有代表性。

本文使用的数据主要有MODIS传感器、Landsat 8上搭载的多光谱陆地成像仪OLI (operational land imager) 和TIRS 3种传感器的数据，数据列表见表 1。研究采用对应区域 (Path：123；Row：32) 的Landsat 8影像，成像时间为2014-05-15、2014-8-19、2014-10-06、2014-12-25，成像效果好，平均云量均低于6%，研究区内无云量，利于进行本文的研究。MODIS数据来自于美国NASA EOS Data Gateway提供的MODIS产品，包括4个时相覆盖研究区 (h26v04) 的500 m分辨率反射率产品 (MOD09GA)、1 km分辨率的每日地表温度产品 (MOD11A1) 和大气水汽产品 (MOD05_L2)，数据均可在NASA官网 (https://eosweb.larc.nasa.gov/) 免费下载。

获取的Landsat 8 OLI_TIRS已经经过几何校正并转至WGS-84/UTM投影下，故只需对其进行辐射定标、大气校正，将OLI多光谱波段DN值转为地表反射率值，使用文献[13]提出的单通道算法来反演地表真实温度，该方法所需参数较少，并且精度和敏感性较高[14-15]。MODIS产品原始的投影为正弦投影，本文使用NASA研发的MODIS Swath Reprojection Tool (MRT Swath) 和MODIS Reprojection Tool (MRT) 分别对MOD05_L2和另外两种MODIS产品将投影转化为WGS-84/UTM，投影带号为50，最后对上述数据进行裁剪，获得研究区数据。

LST降尺度转换总是面临更高分辨率上信息不足的问题，因此常通过引入额外信息构造趋势面才能实现降尺度转换[16]。TsHARP降尺度方法就是通过引入NDVI或植被覆盖度来构造趋势面，其基本假设是：LST和植被指数的关系在各个尺度上基本一致，土壤湿度分布的空间变化是最主要影响因素[17]。以MODIS数据为例，TsHARP方法的原理如下：建立1000 m尺度上LST与趋势面因子间的关系，如式 (1) 所示

式中，T1000表示1000 m尺度上采用趋势面估算的LST值；SI1000表示1000 m尺度上光谱指数，即趋势面因子；f1000(·) 为趋势面转换函数，是趋势面因子与LST之间的函数关系。f1000(·) 仍适用于500 m尺度上趋势面因子与LST之间的转换，它的不变性是将1000 m尺度上的LST降至500 m尺度的关键，也是TsHARP算法的基础假设。但是由于受到土壤含水量等地理要素的影响，趋势面因子总是很难完全反映LST的分布趋势，表现为1000 m尺度上的转换残差ΔT1000

式中，T1000、T1000分别表示1000 m尺度上LST的真实值和估算值。500 m尺度上的LST值，应该由1000 m尺度上的转换残差和500 m尺度趋势面因子SI500带入1000 m尺度上建立的趋势面转换函数f1000(·) 中计算的LST估算值两部分构成，计算公式如下

式中，T500表示500 m的LST降尺度结果。由500 m尺度上光谱指数的转换值，加上由1000 m尺度重采样至500 m尺度的残差得出。

TsHARP算法建立了多种趋势面转换函数，本文选择其中一种模型进行对比研究 (式 (4))。为解决问题2，文中利用决策树模型将地物分为植被、建筑用地、水体和裸地，在保证最小残差和最大决定系数的前提下，分别选择NDVI、NDBI (normalized difference build-up index)、MNDWI (modified normalized difference water index) 和EBSI (enhanced bare soil index) 来建立不同地物光谱指数与LST之间最佳回归关系式 (模型2) 从而实现降尺度，4种指数的计算方法分别参考文献[18-21]，技术路线如图 2所示，两种模型如下：

模型1

模型2

植被

建筑用地

水体

裸地

对1000 m LST到500 m LST的降尺度效果可以从定性和定量两个角度进行。定性分析通过趋势面因子对地物指示性特征分析、对计算结果进行目视对比分析。定量分析包括以下指标。

决定系数 (R2)

相对平均偏差 (MBE′)

相对均方根误差 (RMSE′)

式中, R表示相关系数；mi为500 m尺度上的LST真实参考值；ei表示500 m尺度上的LST估计值；n为样本总数；MBE表示平均偏差；RMSE表示均方根误差。

降尺度结果定性评价如图 3所示，对比图 3中两种模型 (图 3(c)、(d)) 的降尺度结果与原始1 km MODIS LST (图 3(a)) 可以发现，两种降尺度模型均可以捕捉到大部分原始1 km MODIS LST中的地表温度信息，即两幅500 m LST降尺度影像较好地保持了原始地表温度影像热特征的空间分布格局，且影像中高温区和低温区均与图 3(a)中吻合得较好。此外，两种模型的降尺度结果色调更加丰富 (尤其是模型2)，更加精细地描绘了地表热特征的空间差异性，较好地消除了图 3(a)中的马赛克效应。

升尺度后的Landsat 8 TIRS地表温度 (500 m) 如图 3(b)所示，其与1 km MODIS LST的热特征分布基本相同，并且信息更加丰富，影像热特征分布空间变异性更加明显。对比图 3(b)发现，模型2(图 3(d)) 的降尺度结果要比模型1(图 3(c)) 信息更加丰富，同种类型地物温度分布的空间差异性更为明显，温度过渡之间没有明显的过渡圈。为了更加充分验证降尺度的效果，下面对两种模型的降尺度结果进行更进一步的定量统计分析。

为了对降尺度结果进行定量评价，除了计算上述评价指标 (随机选择5000个点)，还对升尺度Landsat 8 TIRS 500 m LST和两种模型降尺度的500 m MODIS LST进行了统计分析，统计量包括均值 (Mean)、标准差 (SD)，此外对Landsat 8 TIRS LST和两种降尺度结果还分别进行了散点图分析，以求更加直观地显示降尺度结果的优良，如图 4所示。统计结果见表 2，由表可知，两种模型降尺度的500 m MODIS LST均值均接近于Landsat 8 TIRS升尺度的LST均值 (25.865℃)，说明两种模型降尺度的结果与实际验证温度之间的总体偏差较小，其中模型1的平均值偏差较大 (1.160℃)。SD反映了图像像素值的离散程度，对于热红外影像而言，SD反映了地表热特征的空间变异性[17]。两种模型的SD计算结果表明，模型2的SD最大，并且最接近于TIRS 500 m LST的SD，说明与模型1相比，模型2能够更好地体现不同下垫面地物的温度差异以及相同下垫面覆盖地物的细节信息，这一点从图 3中也可以看出。两种模型的R2差异并不明显，其中模型2的R2最大。RMSE和MBE能够反映降尺度影像和验证影像之间的一致性，可以用来指示两种降尺度模型与TIRS LST热特征分布的差异性。模型1的RMSE要大于模型2，说明模型2很好地再现了TIRS 500 m LST地表热特征的空间分布格局，能够较好地保持地表真实温度的热特征，而模型1较差，从MBE可以看出，两种模型的降尺度结果与TIRS LST的平均偏差均较小 (绝对值小于1℃)，降尺度结果非常接近真实地表温度，并且两种降尺度结果的MBE均小于0，印证了两种模型降尺度的500 m MODIS LST均值均低于TIRS 500 m LST的均值，这一点由散点图 (图 4) 也可以得到证实，综合对上述统计量的讨论，两种算法均能够达到较好的降尺度效果，不管是在降尺度结果的空间变异性还是精度方面，模型2均要优于模型1。

本文统计了不同地物覆盖类型的降尺度结果 (表 3)，由表 3可知，原始全局降尺度模型 (模型1) 在水体、裸地和建筑用地这些低植被覆盖区表现出较差的降尺度结果，尤其对于裸地和建筑用地更为明显 (|MBE|＞3℃)。两种降尺度模型均表现出对植被、裸地和建筑用地为低估 (MBE＜0)，而水体为高估 (MBE＞0)，图 4也印证了这一点，对整个研究区而言水体表现为低温，低温区位于1:1直线上方，高温区则位于下方，这是由于原始全局降尺度模型在构建趋势面转换函数时剔除了水体，光谱指数二次项的存在会对水体和建筑用地的降尺度结果产生大的干扰[6]。对比模型2和模型1，在植被地区两种模型的差异并不明显，因为在植被地区模型2与模型1使用了相同形式的趋势面转换函数。而对于其他地物类型，新参数的引进和替换均在一定程度上提高了降尺度结果的精度。对于水体地区，模型2的MBE要比模型1提升1.041℃，这是因为原始全局降尺度模型没有考虑水体，而文中改进的模型2，单独建立水体LST与MNDWI之间的关系，较好地避免了其他地物对水体区域的影响。裸地的MBE由原先的-3.203℃提升到了-2.380℃，建筑用地则由-3.463℃提升到了-2.014℃，提升了41.842%，这是由于模型2不仅避开了不同类型地物构建趋势面转换函数时产生的影响，而且对于建筑用地，NDBI与城市中心区的LST拟合效果较好，利用其与LST之间的线性关系对于城市地表温度的降尺度效果要优于单独使用NDVI[16]。两种模型降尺度的500 m MODIS LST与TIRS 500 m LST差异的空间分布如图 5所示，为详细展示不同模型的降尺度差异，文中从研究区选择出一块大小适中的区域，对比图 5(a)、图 5(b)，结合图 5(c)可以发现，模型2对水体、建筑用地的降尺度结果有着明显改善，模型1对水体的降尺度结果明显存在高估，而建筑用地表现为明显低估，空间差异分布情况与上述定量统计具有一致的结论。

不同时相遥感影像反演的LST和4种光谱指数均存在差异，为避免这种差异对降尺度结果评价的影响，本文还分别应用两种模型对不同季节的MODIS LST进行了降尺度，同样以对应时期反演的Landsat 8 TIRS LST作为当期地表温度的参考值对降尺度结果进行验证，验证结果见表 4。由表 4可知，4个季节两种模型的降尺度结果存在明显差异。总体而言，模型2的降尺度结果要优于模型1，但是不同季节提升效果不同。两种模型均是夏、秋季的降尺度结果优于春、冬季，这是由于NDVI受季节影响比其他3个指数要更加严重，并且NDVI温度降尺度方法用于春、冬季要弱于其他方法[5]。此外，分别从RMSE和MBE的差异幅度来看，模型2相比模型1的降尺度结果在春、冬季的提升效果均要优于夏、秋季，这是因为NDVI温度降尺度方法较适合于夏、秋季，因此分土地利用类型使用不同光谱指数和全局尺度上使用NDVI提升效果并不明显，相反，春、冬季处于弱势地位的NDVI被替换后，降尺度效果有着明显提升，此处得出的结论与文献[5]一致。

针对现有研究中地表温度降尺度存在的问题，本文以TsHARP算法为基础，根据地表覆盖类型的不同，分别选择与LST相关性更好的光谱指数建立二者之间的关系，提出了新的转换模型，并用两种转换模型直接将原始1 km MODIS LST产品降尺度到500 m，将Landsat 8 TIRS反演的LST作为当日地表温度的参考值，从定性和定量两个角度评价了两种降尺度模型的降尺度效果及其精度。结果发现：

(1) 对两种模型降尺度结果和原始1 km MODIS LST的对比发现，降尺度影像较好地保持了原始地表温度影像热特征的空间分布格局。两种模型的降尺度结果色调更加丰富，较好地消除了原始1 km影像中的马赛克效应。同样，降尺度结果可以捕捉到大部分TIRS LST影像中的热特征，模型2的降尺度结果要比模型1(TsHARP法) 信息更加丰富。

(2) 定量验证表明，全局尺度上，两种模型降尺度的500 m MODIS LST均值均接近于TIRS升尺度的LST均值，平均偏差的绝对值都小于1℃，两种模型均能够达到较好的降尺度效果，不管是在降尺度结果的空间变异性还是精度方面，模型2(RMSE：1.635℃) 均要优于模型1(RMSE：2.736℃)。

(3) 不同地物降尺度结果的统计表明，原始全局降尺度模型 (模型1) 在水体、裸地和建筑用地这些低植被覆盖区表现出较差的降尺度结果，尤其对于裸地和建筑用地更为明显 (|MBE|＞3℃)。除了植被地区，模型2均在一定程度上提高了其他地物类型的降尺度精度，究其原因，模型2不仅较好地避开了不同类型地物构建趋势面转换函数时产生的影响，而且针对不同的下垫面覆盖类型选择了比NDVI拟合效果更好的光谱指数。

(4) 不同季节的降尺度结果表明，两种模型都是夏、秋季的降尺度结果优于春、冬季，模型2的降尺度结果四季均好于模型1，其中春、冬季的降尺度提升效果要优于夏、秋季。

需要补充说明的是，受限于地面实测气象数据，文中采用相近时间过境的Landsat 8 TIRS LST作为地表实测温度的参考值对4种转换关系的降尺度结果进行定量探讨，虽然这种验证方法应用广泛[15, 17]，但是TIRS LST在升尺度过程中会存在误差，因此，今后还需进行大量试验完善对该方法的验证，此外，针对不同区域、不同分辨率热红外遥感影像降尺度的应用也是今后研究的重点。