甘肃石羊河流域上游地表温度反演及空间特征分析

丁国民，汪有奎，邸 华，车宗玺，李硕军

(甘肃祁连山国家级自然保护区管理局，甘肃 张掖 734000)

甘肃石羊河流域上游地表温度反演及空间特征分析

丁国民，汪有奎，邸 华，车宗玺，李硕军

(甘肃祁连山国家级自然保护区管理局，甘肃 张掖 734000)

采用大气校正法反演研究区地表温度，用地表温度数据分别与经度、纬度、海拔、坡向、坡度进行相关、回归分析，研究地表温度空间分布特征。以2014年7月的2景Landsat8影像为数据源，以甘肃石羊河流域上游为研究对象，应用大气校正法反演了研究区地表温度。结果表明，研究区地表温度在-5.23～50.77℃之间，60.94%区域的地表温度集中在20～30℃之间；研究区地表温度随海拔的升高而降低，呈线性负相关,相关极显著；地表温度随纬度的增大而升高，呈线性正相关,相关极显著；地表温度随坡向的变化呈周期性变化趋势，呈多项式相关,相关极显著；地表温度随经度的增大而降低，相关不显著；地表温度随坡度的增大而降低，线性相关不显著。

大气校正法；地表温度；反演；空间分布特征；环境因子；甘肃石羊河流域

陆地表面温度是地球表面能量平衡和温室效益的一个重要指标。传统获取地表温度的途径只能依靠人工测量，只能获取小范围的点数据；遥感数据范围大、全天候等特点为获取地表温度数据提供了便利手段[1]。地表温度反演算法主要有3种：大气校正法(也称为辐射传输方程：Radiative Transfer Equation——RTE)、单通道算法和分裂窗算法。辐射传输方程描述了卫星的微波辐射计所观测到的辐射总强度,不仅有来自地表的辐射,而且还有来自大气的向上和向下的路径辐射。这些辐射成分在穿过大气层到达遥感器的过程中,还受到大气层吸收作用的影响而削减。同时,地表和大气的辐射也在这一过程中产生不可忽略的影响。因此,地表温度的演算实际上是一个复杂的求解问题[2]。

Sobrino对大气校正法、单窗算法和普适性单通道算法进行了比较分析，结果表明，当使用实时大气廓线数据时，大气校正法的均方根误差(RMSE)是0.6K，单窗算法和普适性单通道算法的均方根误差都是0.9K[3，4]；大气校正法反演地表温度具有相对较高的精度，故本文基于大气校正法，利用Landsat8 TIRS反演研究区地表温度。

1 研究区概况与研究数据

1.1 研究区概况

石羊河是发源于祁连山区的三大内陆河之一，石羊河流域上游的祁连山区属高寒半干旱湿润区：海拔2 000～5 000 m，年降水量300～600 mm，年蒸发量700～1 200 mm，干旱指数1～4。其上游全部分布在祁连山自然保护区辖区内，在保护区内石羊河上游面积约8 851.5 km2，主要由森林、草原、草甸、荒漠、农田、冰川积雪和湿地等自然生态系统构成，是石羊河上游重要水源涵养区，应用遥感技术进行石羊河上游土地地表温度反演研究，具有重要科研价值和现实意义。

1.2 影像来源与处理

目前用于反演地表温度的遥感数据主要集中在NOAA、MODIS、TM等数据中，这些数据空间分辨率都较低，已经不能满足当前应用的需求[4，5]。本文选取Landsat8遥感数据为数据源，成像时间为2014

-07-17和2014-07-26，轨道号是 131/034，132/034。Landsat 8卫星刚发射运行不久，TIRS B11波段暂时存在定标不稳定性，故本文主要针对 Landsat 8 TIRS B10数据讨论大气校正法算法[6]。

1.3 DEM及地形数据

采用“地理空间数据云”公布的GDEMV2 30M分辨率数字高程数据，通过env 4.8软件地形分析功能，用DEM数据提取研究区的高程、坡向、坡度数据。

1.4 数据预处理

首先将2景影像分别进行辐射定标、大气校正，并进行几何精校正。校正后利用envi 4.8软件自定义坐标系，将2景影像、地形数据、石羊河上游矢量边界矢量数据的投影和坐标系统一为transverse mercator，D-Beijing-1954坐标系，应用石羊河流域上游矢量数据对DEM及地形数据进行裁剪。应用波段运算得到归一化植被指数(NDVI)。因2景影像成像的时间不同，大气上行辐射、下行辐射以及大气透过率也不同，应分别进行地表温度的反演，最后将反演的地表温度影像进行镶嵌、裁剪，得到研究区地表温度反演影像数据。

2 研究方法

2.1 LST模型构建与参数计算

在地—气的辐射传输中，卫星接收到的热红外辐射能量Lλ包含3部分内容：地面真实辐射经大气衰减之后被卫星传感器接收到的热辐射能量、大气的上行辐射亮度L↑、大气下行辐射亮度L↓(大气向地面热辐射)。卫星传感器接收到的热红外辐射亮度值Lλ的表达式(辐射传输方程)可写为：

Lλ= [εB(TS) + (1-ε)L↓]τ + L↑

(1)

式中：ε为地表辐射率；TS为地面真实温度，单位为k；τ为大气在热红外波段的透过率。L↓大气下行辐射亮度和L↑大气上行辐射亮度，其单位均为W/(m2·sr·μm)。B(TS)为普朗克定律推算得到的黑体在TS的热辐射亮度。

则温度为T的黑体在热红外波段的辐射亮度B(TS)为：

B(TS) = [Lλ- L↑-τ(1-ε)L↓]/τε

(2)

式中的大气在热红外波段的透过率 τ、大气上行辐射亮度L↑和大气下行辐射亮度L↓单位均为W/(m2·sr· μm)，3个参数可以通过NASA官网(http：//atmcorr．gsfc．nasa．gov/)输入影像的成像时间、中心经纬度等相关信息生成，黑体的辐射亮度B(TS)单位为W/(m2·sr· μm) 。

由于NASA官网暂时只能获取B10波段的参数，而不能获取B11波段的参数，因此统一使用B10波段反演地表温度。估算出地表真实温度相同的黑体的辐射亮度B(TS)后，根据普朗克定律反函数，得出地面真实温度，公式为：

TS=K2/ln(K1/B(TS)+ 1)

(3)

式中：TS为传感器处的地表温度(K)，B(TS) 为黑体在热红外波段的辐射亮度，K1和K2为热红外波段的定标常数，对于TIRS 10波段，K1=774.89W/(m2·sr·μm)，K2=1321.08W/(m2·sr·μm)。K1和K2可在下载影像头文件中查得。

2.2 地表比辐射率及参数估算2.2.1 NDVI及植被覆盖度估算

NDVI指数采用envi软件直接计算法获得。植被覆盖度(Pv)以NDVI数据为基础，采用下式计算。

Pv=[(NDVI- NDVISoil)/(NDVIVeg-NDVISoil)]

(4)

式中：NDVI为归一化植被指数，NDVISoil为完全是裸土或无植被覆盖区域的NDVI值，NDVIVeg则代表完全被植被所覆盖的像元的NDVI值，即纯植被像元的NDVI值，取经验值NDVIVeg= 0.70和NDVISoil= 0.05，即当某个像元的NDVI大于0.70时，Pv取值为1；当NDVI小于0.05，Pv取值为0。

2.2.2 地表比辐射率估算

地表比辐射率是地表温度反演中不可或缺的重要参数，目前对地表比辐射率的估计方法主要是基于NDVI的地表比辐射率法。本文采用覃志豪等提出的地表比辐射率计算方法，先将地表分成水体、自然表面和自然、人工混合表面，分别针对3种地表类型计算地表比辐射率[9]：

ε=0.995 (NDVI≤NDVIv) (水体、冰雪覆盖区)

ε= 0.9589 + 0.086Pv-0.0671Pv2(NDVIv

ε=0.9625 + 0.0614Pv-0.0461Pv2(NDVI>NDVIs)(自然表面)

式中：Pv是植被覆盖度；NDVIs 为完全被裸土或无植被覆盖区域的NDVI 值，NDVIv 则代表完全被植被所覆盖的像元的NDVI 值，即纯植被像元的NDVI值。

2.3 黑体辐射亮度计算

大气校正法进行地表温度的反演，其中涉及到大气上行辐射、下行辐射以及大气透过率数据，这些数据可以在NASA官网(http://atmcorr.gsfc.nasa.gov/)中输入成像的时间及其中心经度，则会得到所需要的参数。将相应的参数带入公式2，即可得到同温度下的黑体辐射亮度图像。

2．4 地表温度反演

将Landsat 8同温度下的TIRS 10黑体辐射亮度图像代入公式(3)，即可获得研究区地表温度。应用envi软件对地表温度进行分级,得到地表温度分级图(图1)。

图1 研究区地表温度密度分级Fig.1 Classification of land surface temperature and density in the research area

2.5 地表温度空间分布特征分析方法

应用ArcGIS10.2软件，对镶嵌、裁剪后的地表温度反演图进行鱼网取样，取样的间隔为0.8 km，获得鱼网点经度、纬度、地表温度数据。用同样的方法应用ArcGIS10.2对海拔、坡度和坡向影像图进行鱼网采样，取样精度小于1个像素，获得相同经度、纬度下的海拔、坡向、坡度数据，用于地表温度空间分布特征分析。用地表温度数据分别与经度、纬度、海拔、坡向、坡度进行相关、回归分析，得到地表温度与经度、纬度、海拔、坡向、坡度的回归模型。

3 结果与分析

3.1 地表温度分级特征

应用envi软件对获取的研究区地表温度进行密度分割(图1)，研究区各温度区间的统计特征(表1)，研究区地表温度在-5.23～50.77℃之间，统计反演结果得出0.84%区域的地表温度在10℃以下，1.98%区域的地表温度在10～20℃之间，60.94%区域的地表温度集中在20～30℃之间，29.95%区域的地表温度在30～40℃之间，6.28%区域的地表温度在40℃以上。

3.2 地表温度与环境因子的相关模型

采用133个有效取样点的地表温度分别与相同点经度、纬度、海拔、坡向、坡度作回归分析，在线性、指数、对数、多项式等模型中选择相关系数最大模型作为拟合模型。结果表明，地表温度与相同位置的海拔、纬度、坡向相关关系极显著，与相同位置的经度、坡度相关关系不显著。并对相关极显著的进行模型拟合，拟合结果见表2，图2 ～图4。

表1 研究区地表温度密度分割区间统计Tab.1 Interval partition of land surface temperature and density in the research area

表2 地表温度与生态因子的相关模型Tab.2 Model of land surface temperature and ecological factor

图2 研究区地表温度—海拔散点图Fig.2 Scatter diagram of land surface temperature- altitude in the research area

图3 研究区地表温度—纬度散点图Fig.3 Scatter diagram of land surface temperature- latitude in the research area

图4 研究区地表温度—坡向散点图Fig4 Scatterdiagramoflandsurfacetemperature-aspectintheresearcharea

地表温度随经度的增加而降低，相关系数r=0.151 0

由图2可知，研究区地表温度随海拔的升高而降低，相关系数r=0.3191>r0.01=0.222,相关极显著，表明线性相关模型成立，地表温度的高低与地表气温相关，地表气温低则地表温度也低，且与实际相吻合。

由图3可知，研究区地表温度随纬度的增大而升高，相关系数r=0.413 9>r0.01=0.222,相关极显著，表明线性相关模型成立。研究区随着纬度增加、海拔降低，地表温度逐步升高，这是海拔、纬度变化双重作用所致。祁连山纬度对气温的影响比较特殊，一般纬度增加，气温会降低，其中由于海拔高度变化的影响，祁连山区出现了纬度增加，气温升高的“逆反”现象。

由图4可知，研究区地表温度随坡向的变化呈现周期性变化趋势，相关系数r=0.319 7>r0.01=0.222,相关极显著，表明多项式相关模型成立。研究区随着坡向的变化，山区阳坡光照时数长,强度大;阴坡光照时数短,强度小[10]。因而阳坡的地表温度、气温变幅及水分蒸发均超过阴坡。阳坡干旱程度超过阴坡，地表植被覆盖比阴坡差；东、西坡接受的光照时长、强度、土壤水分、植被覆盖均介于南、北坡之间，造成地表温度随坡向出现周期性的变化。

4 结语

利用研究区域Landsat 8数据(使用B10波段)结合辐射传导方程法得到石羊河流域上游地表温度，反演结果表明，研究区60.94%的区域处于中温度区域。高温区、低温区所占的比例分别为6.28%和0.84%，所占的面积均较低。

用地表温度数据分别与经度、纬度、海拔、坡向、坡度进行相关、回归分析。结果表明，地表温度随海拔的升高而降低，呈线性负相关,相关极显著；地表温度随纬度的增大而升高，呈线性正相关,相关极显著；地表温度随坡向的变化呈周期性变化趋势，呈多项式相关,相关极显著；地表温度随经度的增大而降低，相关不显著；地表温度随坡度的增大而降低，线性相关不显著。

由于没有同步实测数据，无法对反演的地表温度作出定量的评价，但研究得出地表温度与地形因子的经度、纬度、海拔、坡向、坡度之间的相关关系和回归模型及地表温度空间分布特征，仍有一定的参考价值。

[1] 蒋大林，匡鸿海，曹小峰，等.基于Landsat 8的地表温度反演算法研究—以滇池流域为例[J].遥感技术与应用，2015,30(3)：448-454.

[2] 宋金红,陈圣波,包书新,等.地表温度反演的数值方法[J].吉林大学学报(地球科学版)，2007(增刊)：198-201.

[3] 历华，曾永年，负培东,等.利用多源遥感数据反演城市地表温度[J].遥感学报，2007，11(6)：891-898.

[4] 樊辉．基于Landsat TM的城市热岛效应与地表特征参数稳健关系模型[J]．国土资源遥感，2008，77(3)：52-57.

[5] 邹蒲，李婷，邹娟.基于Landsat8卫星数据的城市热岛效应分析—以湘潭市为例[J]．环球人文地理(评论版)，2015(3)：32-34.

[6] 胡德勇，乔琨，王兴玲,等.单窗算法结合 Landsat 8 热红外数据反演地表温度[J].遥感学报，2015，19(6)：964-976.

[7] 吴志刚，江滔，樊艳磊,等.基于Landsat 8数据的地表温度反演及分析研究—以武汉市为例[J].工程地球物理学报，2016，13(1)：135-142.

[8] 游绚，晏路明.基于ETM+ 影像辐射传导方程算法的地表温度反演[J].科技情报开发与经济，2009，19(27)：134-136．

[9] 毛文婷,王旭红,祝明英,等.城市地表温度反演及其与下垫面定量关系分析—以西安市为例[J].山东农业大学学报(自然科学版),2015,46(5):708-714.

[10] 徐化成.景观生态学[M].北京：中国林业出版社,1995：1-83.

Inversion and Spatial Characteristics of Land Surface Temperatures inthe Upper Shiyang River Basin of Gansu

DING Guomin, WANG Youkui, DI Hua, CHE Zongxi LI Shuojun

(Gansu Qilian Mountains National Natural Reserve, Zhangye, Gansu 734000,China)

In order to study the spatial distribution characteristics of land surface temperatures, we usually study the land surface temperature by method of atmospheric correction, and make correlation and regression analysis based on the land surface temperature data with longitude, latitude, altitude, aspect, and gradient respectively. Taking 2 Landsat 8 images of July 2014 as the data resource, the upper Shiyang river basin of Gansu as the research object, surface temperatures were displayed by the method of atmospheric correction. The results showed that land surface temperatures in the research area ranged from -5.23 to 50.77, among which land surface temperatures of 60.94% areas were between 20 and 30 degrees; the land surface temperature decreased with the increase of altitude, presenting the remarkable linear negative correlation, heightened with the increase of latitude, presenting the remarkable linear positive correlation, showed periodic variation trend with the change of aspect, markedly presenting the polynomial correlation, decreased with the increase of longitude, inapparently presenting the correlation; and decreased with the increase of gradient, indistinctively presenting the linear correlation.

atmospheric correction; land surface temperatures; inversion; spatial distribution characteristics; environmental factors ;Shiyang river basin of Gansu

2017-03-30.

2015甘肃省林业科技项目(项目编号：2015kj043).

丁国民(1966-)，男，甘肃人，高级工程师.从事保护区生态保护及科研工作.

10.3969/j.issn.1671-3168.2017.03.005

S716.2；P423

A

1671-3168(2017)03-0017-05