黄河源永曲河流域地表温度与土地覆盖的时空变化特征分析

丁圆圆， 赵健赟*, 杨静, 赵沁浩, 王祖顺, 赵利江

(1.青海大学地质工程系， 西宁 810016； 2.青海省刚察县气象站， 刚察 812300； 3.青海省基础测绘院, 西宁 810001)

地表温度是地表与大气相互作用过程中最重要的物理学参量之一[1]，在全球和局部气候变化、生态环境、城市热岛效应等研究中具有重要意义[2-3]，其时空分布与生态环境、地形因子、土地覆盖类型等紧密相关，因此，关于地表温度的反演及其影响因素的研究受到了越来越多学者的关注,已经成为全球变化领域的研究热点之一。

目前中外已经开展了较多地表温度领域的研究，大范围地表温度数据一般通过遥感影像进行反演，方法主要有辐射传输方程法[4-7]、单窗算法[8-9]、劈窗算法[10-11]。对地表温度变化及其影响因素方面的研究也已取得了一定的成果，例如，王佳等[12]分别利用普通最小二乘法(ordinary least square, OLS)线性回归、地理加权回归(geographically weighted regression, GWR)模型拟合土地覆盖比例和地表温度的关系，结果表明不同土地覆盖类型的地表温度存在显著差异，且两者之间的定量关系受到空间异质性的影响。毛洋等[13]和李军等[14]认为地表温度的时空分布与地表类型密切相关。另外，有不少学者对地表温度随地形因子的时空变化规律也进行了研究，赵伟等[15]得到了地表温度随海拔高度、坡度的增加而逐渐增加的结论。

李春波等[16]运用统计学方法分析阜平县地表温度与坡向之间的关系，结果表明不同坡向的地表温度分布具有一定的差异,阳坡的地表温度高于阴坡。关于归一化植被指数(normalized difference vegetation index, NDVI)与地表温度的关系也引起了不少学者的关注，但是得到的结论不完全一致，许剑辉等[17]通过GWR模型拟合二者之间的相关关系，得到了NDVI聚集密度对地表温度有着从正到负的影响的结论[17]；Mishra等[18]和Malik等[19]认为地表温度随植被覆盖度的变化而变化。张晓娟等[20]则认为地表温度与NDVI呈负相关。在青藏高原气候暖湿化的背景下，黄河源作为高寒生态系统敏感和脆弱地带，其地表覆盖、温度等环境气候因子的变化引起了广泛的关注，而以上研究主要集中于低海拔地区地表温度的相关研究中，缺少黄河源高寒草甸覆盖区域地表温度时空变化特征与影响因素的研究。本研究以黄河源区永曲河流域为例，同时兼顾该区域植被生长特征、时间尺度、遥感数据源的获取与质量等问题，选用2008年、2020年两期Landat遥感影像反演地表温度，并与土地覆盖、地形因子和NDVI的关系进行分析和探讨，以获得研究区地表温度的时空变化特征及其规律。

1 研究区概况

研究区为青海省东南部黄河源地区的永曲河流域，处于101°20′～101°40′E，34°10′～34°30′N，如图1所示，位于青海省东南部，隶属黄南藏族自治州，是黄河流域重要的水源地，总面积约为745.76 km2。海拔在3 317～4 451 m，总体地势表现为东北高西南低，且东部小部分地区常年有冻土积雪分布。流域气候为高原大陆性气候，属于高原亚寒带半湿润气候区，由于海拔较高，地势复杂和受季风影响，高原大陆性气候特点比较明显。每年5—10月份温暖、多雨，11月—次年4月份寒冷、干燥、多大风天气。常年风向为西北风，最大风速可达到23.7 m/s，年平均风速约为2.2 m/s，年平均降水量在400 mm以上，植被类型以草地为主，且仅在高海拔地区零星分布小灌木丛。

图1 研究区遥感影像图和高程Fig.1 Remote sensing image and altitude in the study area

2 数据与研究方法

2.1 数据来源与预处理

本文中使用的数据主要包括 Landsat、数字高程模型(digital elevation model, DEM)、土地利用类型和气象站数据。其中Landsat影像来源于地理空间数据云，轨道号为132/36，空间分辨率为30 m，TM5和TM8成像时间分别为2008年8月 17日和2020年8月2日;DEM数据来源于地理空间数据云，经纬度为N34°E101°，空间分辨率为30 m；土地覆盖数据来源于中国科学院资源环境科学与数据中心，空间分辨率为30 m。

对收集到的数据进行了去云、辐射定标、大气校正、投影、重分类等处理，获得了研究区2008、2020年的30 m分辨率影像数据和DEM数据。

2.2 LST反演算法与原理

研究采用辐射传输法反演地表温度，该方法可用于任何热红外遥感波段，其基础理论均基于普朗克定律的热辐射传输方程。为提高地表温度反演精度，中外学者们开展了大量的研究，例如，为解决地面物体某些性质不能确定的问题，将“黑体”引入辐射传输理论；进而又考虑到灰体在自然界的普遍存在，将地表比辐射率应用于辐射传输模型，所需参数易于获得，计算过程较为简单，反演精度高，结果可靠[21-22]。辐射传输法先估计大气对地表热辐射的影响，然后把这部分大气影响从卫星传感器所观测到的热辐射总量中减去，进而得到地表热辐射强度，再把这一热辐射强度转化为相应的地表温度。

2.2.1 热红外波段辐射亮度的计算

对Landsat热红外波段进行辐射定标处理，得到热红外波段的辐射亮度Lλ[23]为

Lλ=GainDN+Bias

(1)

式(1)中：DN表示传感器记录的像元亮度值；Gain表示增益，W/(m2·μm·sr)；Bias表示偏移，W/(m2·μm·sr)。

2.2.2 地表比辐射率的计算

(1)归一化植被指数的获取。归一化植被指数NDVI是反映植被生长状态、地表覆盖植被状况的一种遥感指标，可以通过近红外和红光波段获得，即

(2)

式(2)中：NIR为遥感多波段图像中的近红外波段；R为红波段；NDVI位于-1～1，负值表示地面覆盖为云、水、雪等，0表示有岩石或裸土等，正值表示有植被覆盖，且随覆盖度增大而增大。

(2)植被覆盖度的计算。将整个影像的地类大致分为水体、植被和建筑，采用混合像元分解法计算植被覆盖度FV，即

(3)

式(3)中：NDVIV和NDVIS分别表示植被和裸地的NDVI值。NDVIV=0.7，NDVIS=0.05。当NDVI>NDVIV时，FV取值为1；NDVI

(3)地表比辐射率的计算。将水体像元的比辐射率赋值为0.995，自然表面和城镇像元的比辐射率[24]为

ε自然表面=0.962 5+0.061 4FV-0.046 1FV2

(4)

ε城镇=0.958 9+0.086FV-0.067 1FV2

(5)

2.2.3 相同温度下黑体辐射亮度的计算

依据Landsat的成像时间以及中心经纬度得到大气剖面参数的向上辐射亮度Lup、大气向下辐射亮度Ldown、大气在热红外波段的透过率τ，计算B(TS)为

B(TS)=[Lλ-Lup-τ(1-ε)Ldown]/(τε)

(6)

式(6)中：B(TS)表示温度为TS的黑体辐射亮度。

2.2.4 地表温度的计算

由普朗克公式的反函数反演获得真实地表温度TS为

TS=k2/ln[k1/B(TS)+1]

(7)

式(7)中：k1、k2为是辐射定标参数，随传感器的不同，而有不同的取值。对于Landsat 5数据,k1=607.76 W/(m2·μm·sr),k2=1 260.56k；k是热力学温度，K。对于Landsat 8数据，k1=774.89 W(m2·μm·sr)，k2=1 321.08k。

2.3 地理加权回归分析

地理加权回归模型通过建立空间范围内每个点处的局部回归方程，来探索研究对象在某一尺度下的空间变化及相关驱动因素，并可用于对未来结果的预测。它是对普通线性回归模型的扩展，其显著的特征是考虑到了空间异质性，将数据的地理位置嵌入到回归参数中，因此具有更高的准确性，其公式为

i=1,2,…，n

(8)

式(8)中：(ui,vi)是第i个样点的坐标；εi是第i个样点的随机误差；β0是第i个样点的截距常量，βk(ui,vi)是第i个样点的第k个回归参数，是关于地理位置的函数，在估算的过程中采用权函数的方法得到。本文中从局部出发，以地表温度(land surface temperature, LST)为因变量，NDVI为解释变量，并采用高斯函数作为空间权重矩阵函数和修正后的Akaike信息量准则(Akaike information criterion, AIC)确定带宽，分析空间异质性分布的变量之间的线性关系。

3 研究结果与分析

3.1 地表温度时空变化特征

利用辐射传输法对永曲河流域2期Landsat影像做地表温度反演，结果如图2所示，并统计反演结果的最大值、最小值和均值，结果如表1所示。从表1可知，2020年地表温度的最大值、最小值和均值都要高于2008年，其中最小值增加了5.147 ℃，变化最大，这表明2020年的地表温度相对于2008年整体呈显著上升趋势。

图2 2008年、2020年地表温度反演结果Fig.2 Inversion results of LST in 2008 and 2020

表1 2008年、2020年地表温度统计结果

3.2 地表温度影响因素的分析

地表温度受多方面因素的影响，为了分析永曲河流域地表温度变化的原因，利用DEM和土地利用类型数据，从高程、坡度、坡向、土地覆盖类型4个方面定量分析地表温和NDVI的时空分布和变化特征。

3.2.1 地表温度的变化特征

根据研究区的地形特征，将各地形因子按以下方式进行等级划分。

(1)将海拔高度每隔400 m分为一个等级，研究区的海拔高度划分为3 317～3 736 m、3 736～3 959 m、3 959～4 451 m共3个等级，如图3所示。

图3 高程分级图Fig.3 Image of altitude classification

(2)将坡度划分为平坡 (0°～5°)、缓坡 (5°～15°)、斜坡(15°～25°)、陡坡(25°～35°)和急坡(35°～90°)共5个等级，如图4所示。

图4 坡度分级图Fig.4 Image of slope classification

(3)坡向以正北方向为0°，将0°～360°划分为平缓坡(0°)、阴坡(0°～45°和315°～360°)、半阴坡(45°～90°和270°～315°)、阳坡(135°～225°)、半阳坡(90°～135°和225°～270°)共5个等级，如图5所示。

图5 坡向分级图Fig.5 Image of aspect classification

对不同高程、坡度、坡向等级下的地表温度和NDVI进行统计，如表2、表3所示,统计表明：

表2 2008年地表温度和归一化植被指数分类统计Table 2 Classification of LST and NDVI in 2008

表3 2020年地表温度和归一化植被指数分类统计Table 3 Classification of LST and NDVI in 2020

(1)海拔与地表温度、NDVI呈现负相关的趋势，即海拔越高的地方，地表温度和NDVI越低；当海拔低于3 959 m时，随着海拔升高，地表温度和NDVI下降幅度较小，海拔大于3 959 m时，地表温度和NDVI下降较快，说明二者受海拔影响更显著。

(2)坡度也会对地表温度和NDVI造成一定的影响：地表温度和NDVI随坡度增加而降低，平坡处的地表温度和NDVI是最高的，急坡处的地表温度和NDVI最低，平坡处温和的环境有利于植物的生长，植物长势良好，相应地，该处的NDVI也较高；

(3)按照阳坡、平缓坡、半阳坡、半阴坡、阴坡的顺序，地表温度和NDVI逐渐降低；2008年阴坡和阳坡之间的地表温度相差1.169 ℃，2020年阴坡和阳坡之间的地表温度相差0.626 ℃，其他坡向之间的差异较小。

3.2.2 不同土地覆盖类型的地表温度变化特征

地表温度与土地覆盖类型紧密相关，为分析两者之间的关系，将研究区的土地覆盖类型与地表温度做分区统计,如表4所示。由于研究区特殊的地理位置和环境，土地覆盖主要为草地和小灌木丛两大类，研究区内低海拔地区以草地覆盖为主，周边高海拔地区仅有小灌木丛零星分布。统计结果表明：低温区的分布与小灌木丛的分布在空间上基本保持一致；小灌木丛主要分布在3 959～4 451 m的高海拔地区，草地主要分布在3 736～3 959 m的低海拔地区，高寒环境下，小灌木丛的生长和覆盖度受到抑制，从而造成了小灌木丛的NDVI要比草地的NDVI低的现象，这说明NDVI受海拔影响很大，且二者呈现出负相关的关系。在2008、2020年，小灌木丛处的地表温度都较低，温度最大相差2.043 ℃，而草地处的地表温度变化则很小，但总体上温度呈现上升的趋势。

表4 不同土地覆盖类型的地表温度和归一化植被指数分类统计Table 4 Classification of LST and NDVI for different land cover types

将反演的地表温度分为低温区(<21 ℃)、常温区(21～27 ℃)、高温区(>27 ℃)，并分析2008年、2020年的不同土地覆盖类型占各温度区间的情况，如表5所示，对比2个年份不同地类的温度区间面积，按照低温区、常温区、高温区的顺序，2008年各温度区间所占研究区的面积比例分别为9.147%、77.163%、13.69%，2020年各温度区间所占研究区的面积比例分别为3.583%、90.374%、6.043%，从这两年不同温度区间所占的比例来看，2020年常温区面积比例增加了13.211%，低温区间和高温区间下降的幅度基本相当，进一步表明研究区地表温度呈现上升的趋势，这与2008年、2020年的研究区地表温度反演结果保持一致。

表5 不同土地覆盖类型在各温度区间所占的比例Table 5 Proportion of different land cover types in each temperature range

3.3 地表温度与NDVI的地理加权回归分析

3.3.1 NDVI空间格局及变化特征

NDVI是反映地表植被覆盖状况的一种遥感指标，包含了研究区的植被信息、水体信息、建筑信息等。本文中分析了2008—2020年的研究区NDVI的分布状况，如图6所示。由图6中可知， 2个年份的NDVI均具有明显的空间分布特征，即大部分地区的NDVI值都比较高，只有东部地区的NDVI较低；从2个年份的NDVI统计值来看，NDVI呈现出上升的趋势，即2020年研究区植被覆盖度相对于2008年有所增加，且在研究区的东北部和西南部增加显著。

图6 2008年、2020年归一化植被指数空间分布Fig.6 The space distribution of NDVI in 2008 and 2020

3.3.2 NDVI对气温的影响特征

对2008年和2020年的LST与NDVI构建GWR模型，并对回归系数做地图可视化,如图7所示。从图7中可以看出：NDVI回归系数为正值，这表示NDVI对地表温度为正影响；2008年NDVI回归系数在0.37～9.77变化，2020年的NDVI回归系数主要在7.28～8.62，这说明2008年NDVI对地表温度的影响程度的变化较大，2020年NDVI对LST的影响程度比较集中，变化较小；从回归系数的空间变化来看，整体上回归系数从西到东逐渐增加，即NDVI对地表温度的影响逐渐增强。

图7 2008年、2020年归一化植被指数的回归系数Fig.7 Regression coefficients of NDVI in 2008 and 2020

拟合优度(LocalR2)表示自变量对因变量的解释程度，研究获得的LocalR2的空间分布如图8所示,较大值主要分布在东部地区，较小值分布在西部地区，且2个年份的LocalR2值均小于0.5，这表示NDVI对LST的解释程度较低，受空间异质性影响较强烈，也即是除了NDVI影响地表温度外，还有其他因素对地表温度有影响。

图8 2008年、2020年地理加权回归模型局部拟合度Fig.8 Degree of local fit of GWR models in 2008 and 2020

3.3.3 GWR残差分析

以1 km×1 km尺度为例，以LST为因变量，NDVI为解释变量，对2008年、2020年的回归方程残差进行分析，结果如图9所示，可以看出，研究区东南部地区的标准化残差超过2.5倍，拟合效果不太理想。考虑到残差的空间自相关性会对分析结果产生一定的影响，为进一步分析和验证结果，使用全局Moran’sI(莫兰指数)检验残差的空间自相关性，经计算，2008年和2020年的Moran’sI分别为0.536和0.492，这说明LST与NDVI的关系受空间异质性影响十分强烈。

图9 2008年、2020年地理加权回归模型残差分布Fig.9 Distribution of GWR residuals in 2008 and 2020

另外，目前关于地表温度和NDVI的研究主要通过相关系数和普通回归模型(如最小二乘法)，但在研究区不同的空间位置， NDVI的变化对地表温度的影响程度可能并不一致。本文研究使用Moran’sI和地理加权回归分别从全局和局部的角度来探讨地表温度与NDVI的相关关系。由于地理加权回归引入了空间权重，能有效地定量研究空间非平稳性和较好地解决研究对象的空间异质性问题，避免了传统线性回归模型不能反映回归参数的真实空间特征的局限性。

4 结论

利用辐射传输法，反演了永曲河流域2008年、2020年的地表温度，并分析了地表温度的时空变化特征及其与地形因子、土地覆盖类型和NDVI的相关关系，得到了以下结论。

(1)2020年研究区地表温度相对于2008年总体上存在较为显著的上升趋势。

(2)地表温度和NDVI受地形因子影响表现出不同的空间分布特征：平均温度和NDVI随着海拔和坡度的增加而逐渐降低，且海拔越高负相关性越显著；不同坡向处的地表温度和NDVI也存在一定的差异，且二者在阳坡和阴坡处的差异更显著。

(3)研究区2008年、2020年NDVI与地表温度呈正相关关系。地理加权回归残差具有较强的空间自相关性，NDVI与地表温度的关系受空间异质性影响较强烈，而地理加权回归法充分考虑了对象的权重问题，能有效获得对象的非平稳关系和异质性特征，具有一定的优势。

研究结果对黄河源的生物多样性、植被保护、水资源安全问题等有重要意义，也能为改善该区脆弱的生态环境提供一定的参考。本文中采用辐射传输法反演地表温度，仍缺乏地面验证，且GWR模型中只分析了1 km×1 km尺度上NDVI与地表温度的关系，对于二者在不同空间尺度上的定性定量分析有待进一步研究。