从气候学意义上看，影响温度空间分布的因素主要包括宏观地理条件(经、纬度以及宏观的气候背景条件等)、局地地形条件(海拔、坡度、坡向、遮蔽度等)以及下垫面性质(土壤、植被状况)等[1]。目前在下垫面性质影响地表温度特征的研究较多。地表热特征主要由土壤湿度和植被覆盖确定[2, 3]。因此，许多学者应用植被指数(NDVI)来定量分析地表温度(LST)[4, 5]，并大都发现两者之间存在着明显负相关关系[6, 7, 8, 9]。但不少研究同时发现NDVI在表征地表温度上有其局限性，地表温度与NDVI特征空间有时呈三角形，难以用简单的函数关系加以描述[10, 11, 12]。为了能够更好的量化与表征地表植被特征，有些学者开始探寻新的替代性指标。Wilson[13]在森林管理与采伐的研究中发现应用归一化水汽指数(NDMI)比NDVI指数更容易鉴别出不同的采伐或森林干扰类型。钱乐祥[3]在城市热岛的研究中也发现NDMI比NDVI更适合表达地表热特征。然而，以山区作为研究区比较NDVI和NDMI对地表温度影响的研究目前并不多见。

除了下垫面性质之外，地形也是地表温度的重要影响因素[14, 15]。Bailey[16]研究发现海拔每升高1000 m，温度平均下降6.4 ℃，但这一规律具有地域特征。Wilson[17]指出，除了海拔，地表温度可以表示成纬度、坡度、坡向、地形阴影和时间的函数。因为太阳入射能量对维持地表能量平衡具有重要影响，特别是在地形复杂的地区，影响更加显著[18]。Martin[15]将山区植被冠层温度表示成海拔与太阳入射角的函数，并获得了较理想的结果，为研究地形对地表温度的影响提供了依据。然而，目前国内对于山区地表温度反演及影响因素特征分析的研究还不多见。

本文将以泰山作为研究对象，利用泰山的TM卫星遥感影像数据、ASTGTM2高程数据、MODIS影像，在综合考虑了泰山地区的下垫面特征、局地地形条件、太阳入射能量的基础上，首先分别探讨了每类因子与地表温度的统计学特征，然后利用逐步回归法综合分析了地表温度和归一化指数、地形因子之间的关系，以期为了解山区地表温度空间分异特征及形成机制提供科学的理解与借鉴。

泰山是中国五岳之首，主峰玉皇顶海拔1545 m，位于泰安市城北，地理位置117°06′15″E，36°15′17″N。泰山地势差异显著，地形起伏大，总体地势呈现北高南低、西高东低的特征，在不到10 km的水平距离内，玉皇顶与其山前平原相对高差达1300 m以上。泰山的山麓区域气候属于暖温带半湿润大陆季风气候类型，四季分明，夏季炎热多雨，冬季寒冷干燥。气温日、年较差比较平缓，具有高山气候的特点。气温垂直变化明显，山下1月均温-3 ℃，山顶-9 ℃，山下7月均温26 ℃，山顶18 ℃。泰山森林类型主要包括针叶林、落叶阔叶林、针阔混交林等，植被覆盖率90%以上。主要植被种类有麻栎(Quercus acutissima)、栓皮栎(Quercus variabilis)、 刺槐(Robinia pseudoacacia L)、枫杨(Pterocarya stenoptera)、 侧柏(Platycladus orientalis)、油松(pinus tabuliformis)、 落叶松(Larix gmelinii)、苔草(Carex spp)等。本文选取的研究区位于济南市的历城区、长清区与泰安市泰山区、岱岳区4个市辖区交界处，研究区边界南面基本上与200 m等高线重合，北面沿山谷划定，总面积147.85 km2。

本文的主要基础数据包括：1)2006年6月19日的Landsat 5 TM影像，条带号122—35，其中TM6热红外波段分辨率120 m，从USGS网站获取时，已被重采样至60 m，其余波段30 m，云量0%；2)MODIS影像，MOD07 L2为MODIS每天的大气剖面数据，空间分辨率为5 km，影像拍摄时间比TM影像晚15 min，大气状况短时间内的影响可以忽略，因此利用MOD07估计单窗算法的大气含水量；3)ASTGTM2 DEM数据，于2011年10月发布的，全球水平空间分辨率为30 m，海拔精度17 m。

NDVI与NDMI的物理依据是地物反射率的差异，而TM影像存储的为DN值，所以需要对TM影像进行辐射校正，将DN值转换成真实的地表反射率。MODIS产品原始的投影为正弦投影，本文利用ENVI地理校正工具将投影转换为经纬度投影。海拔数据直接从DEM中提取，坡度和坡向数据由ArcGIS软件生成。由于直射像元接受的辐射量正比于太阳入射角i的余弦值[19]，所以用它的余弦值表示其大小，其计算公式如下[20]：

式中,θp、Φ0为坡度、坡向，θz、Φa分别为TM影像拍摄时的太阳天顶角与方位角。θz、Φa可通过查找TM影像的头文件得到，分别为23.58°，116.03°。

所用TM为L1T数据，已经经过了地形校正，并且TM影像与ASTGTM2投影相同，因此，无需再对两者进行匹配，可以直接使用。

地表温度反映的是传感器在卫星高度观测到的热辐射强度相对应的地面温度[21]。本研究应用TM6的热红外波段并主要根据覃志豪等[22]推导出的单窗算法来反演泰山的地表温度。单窗算法依据地表热辐射传导方程，利用地表辐射率、大气透射率和大气平均温度3个参数推导出的一个简单易行的演算方法，来进行地表温度的演算。大气透过率采用MOD07的大气廓线数据进行估计，其他的演算公式以及参数取值详见文献22，此处不再累述。当基本参数估计没有误差时，演算绝对精度 < 0.4 ℃，当参数有适度误差时，演算精度仍达 < 1.1 ℃[22]。当利用MODIS反演的大气参数应用于TM地表温度反演时，单窗算法的反演精度为0.76 ℃[23]。与1 km的MODIS地温产品相比，通过单窗算法反演的地表温度具有高空间分辨率的优势。将NDVI、遥感影像分类图、TM6波段以及各个参数代入演算公式中，借助ArcGIS空间建模工具，获得研究区的地表温度(图 2)。

由于所选研究区属于区域性山脉，因此，只考虑局地地形条件。Martin[15]在研究地形对植被冠层温度影响时，选取了海拔与太阳入射角两个地形因素。Martin研究认为太阳入射角是由坡度、坡向以及太阳高度角、方位角计算得出，因而太阳入射角就包含了这些因素所携带的信息。但是这种做法的科学性值得商榷，因为坡向、坡度是决定植被空间分布的主要原因[24]，而植被可以引起地表温度空间分异，即坡度、坡向可以通过影响地表性质，间接影响到地表温度。而太阳入射角只是反映了栅格单元间所接收太阳入射能量的差异，不能反映地表性质的变化。综合考虑，本文中选定了4个地形因子：海拔、太阳入射角、坡度、坡向。

地表热特征与植被紧密相关，而NDVI是表征地表植被覆盖特征的一种常用的植被指数。因此，NDVI在一定程度上可以反映地表温度的特征。其值域范围[-1,1]，正值对应着植被冠层的密度和绿色，值越大代表植被的覆盖度越高；0代表该区域基本没有植被生长；负值代表非植被覆盖的区域。其计算公式为[25]：

式中，NIR、dsat5Rdsat5表示TM4近红外波段、TM3红光波段。

NDMI是利用近红外与短波红外之间的差异来表达的，与NDVI相比，NDMI与植被冠层水汽含量高度相关[26]，表征地表水汽特征。值域范围[0, 1]，高值对应着水汽含量高的植被冠层、水体，低值对应着岩石、建筑物、构筑物等表面水汽含量低的地表物质[3]。其计算公式[27]：

式中,MIR表示TM5短波红外波段。

本文选取了两类影响地表温度的因子：地形因子与地表特征指数NDVI、NDMI。首先分别分析这两类因子与地表温度之间的关系，每个因子与地表温度之间的相关程度采用Pearson相关系数来描述。由于很多研究表明水体具有明显的降温效应，往往在地表形成绝对“冷岛”[28, 29]，因此在研究地表温度与地形的关系时，对研究区水体进行了掩模处理，由于研究区像元数太多，为较清晰地反映地形因子与地表温度间的关系，所以采用300 m×300 m渔网采样。NDVI与NDMI都是反应地表特征的指数，由于研究区地表类型分布不均匀，水体较少，使用渔网采样不能较好地反映整体的特征。因此，计算逐像元条件下的LST-NDVI与LST-NDMI的Pearson相关系数，进行对比分析。然后，依据地表温度与每个因子相关系数的大小选取重要因子，与地表温度间进行单独分析。最后综合所有可选因子采用逐步回归分析法构建最优的地表温度回归方程，回归方程中由于变量较多，使用简单相关关系(Pearson相关系数)已不能真实地反映单个因子与地表温度的关系，采用偏相关系数来判断变量的相对重要性[30]。分析统计过程在SPSS软件中实现。

为直观地展现泰山地表温度的分布特征，将地表温度反演图与ASTGTM2高程图叠加生成LST随地形变化的三维图 (图 3)。

从图 3可以看出，由边缘向中心，颜色逐渐由红色向绿色过渡，即地表温度随着海拔呈逐渐降低的趋势。为清晰准确地表示地表温度与地形之间的关系，将基于回归分析，进一步揭示它们之间的统计学规律。

将选定的4个地形因子与地表温度间进行相关分析，分析结果见表 1。

从表 1可以看出地表温度与坡度之间的相关性强于与太阳入射角(cosi)之间的相关性，这也说明提前剔除掉坡度、坡向因子是不合理的。地表温度与坡向之间的皮尔森相关系数最差，与海拔之间相关性最好。对除水体之外整个研究区用渔网采样后，构建了地表温度与海拔之间的散点图(图 4)。

从图 4散点图可以看出，随着海拔的升高，地表温度呈明显的降低趋势。为定量分析地表温度与海拔之间的关系，又对两者进行了线性与自然对数拟合。如图 4所示，在海拔低于900 m左右，线性关系能很好地反映两者之间的关系，然而随着海拔继续升高，地表温度随海拔的变化率在逐渐的减小，自然对数的拟合分析整体效果较好。由此，可以得出，地表温度随着海拔的升高，并非是简单地线性降低，而是呈自然对数形式递减。

为进一步分析地表温度与地形因子之间的关系，应用逐步回归分析法建立地表温度与地形因子的回归方程，其中海拔取自然对数值，结果如表 2所示。

由表 2中可以看出，常数及其各变量都通过了T检验，方差膨胀因子VIF都较小，因子间没有多重共线性。由偏相关系数一列，可以得出：地表温度与海拔、坡度间呈负相关关系，与太阳入射角(cosi)、坡度呈正相关关系；影响地表温度的主要地形因素是海拔，坡度、坡向以及太阳入射能量对地表温度影响均非常小。

将表 2中的“非标准回归系数”列数据代入多元回归模型得到地形因素影响下的最优的逐步回归模型：

对研究区内的所有像元进行地表温度与NDVI、NDMI的相关分析，Pearson相关系数分别为-0.742、-0.830，说明地表温度与两者之间均存在着高度负相关关系。为直观描述地表温度随NDVI、NDMI的变化规律，进一步统计了LST-NDVI与LST-NDMI的散点图(图 5)。

从图 5可以看出这两个指数与地表温度的散点图中都存在“离群值点”，不同的是NDVI位于图 5的低值处，而NDMI位于图 5的高值处。将地表温度图与遥感影像图叠加进行对比，可以发现LST-NDMI之间的“离群点”都是水体纯净像元(图 5黑色椭圆实线内)；而LST-NDVI的“离群点”不但包括水体纯净像元(图 5黑色椭圆实线内)，还有水体与自然地表混合像元(图 5黑色椭圆虚线内)。水体与自然地表的混合像元如图 6红色框所示，1个热红外像元包含4个可见光像元。图 6中红色框内，NDVI为-0.067，LST为27.81 ℃，处于图 5黑色虚线内。另外，由图 5还可以看出水体以及水体与地表混合像元的LST与NDVI之间存在着正相关关系，这与梁保平[21]结论一致。

对水体进行掩模及排除掉水体与自然地表混合像元后，分别建立LST-NDMI和LST-NDVI的回归方程，结果如图 7所示。

由图 7可以看出，地表温度与NDVI、NDMI之间都具有良好的负相关函数关系。地表温度与NDVI之间关系是非线性的，随着NDVI的增大，地表温度降低的变化率在逐渐减小，两者之间的最佳拟合为二项式拟合(图 7)，这也表明了地表温度在高植被覆盖区(NDVI > 0.5左右)的“饱和”现象[31, 32]。图 7显示的地表温度与NDMI的散点图及回归方程可以看出地表温度与NDMI之间呈稳定而简单的线性关系，且复相关系数(R2 = 0.691)较地表温度与NDVI(R2 = 0.650)之间的高。

综上所述，在山区，当没有水体时，NDVI、NDMI与地表温度均存在着很强的负相关函数关系；当有水体时，NDMI在使用时，需要对水体进行掩模，而NDVI还需要对水体做一定距离的缓冲区，将混合像元也排除掉。NDMI适用范围比NDVI广，且与地表温度的拟合效果优于NDVI。结合钱乐祥[3]平原城市中的研究结果，可以得出在夏季，表征地表水汽特征的NDMI比表征地表植被覆盖状况的NDVI更适合于表达地表热特征。NDVI与NDMI虽然表征的内容上有所不同，但两者之间相关系数为0.783，具有高度的相关性，为避免变量间的共线性问题，因此，选取NDMI作为描述地表性质的参量，参与下面的回归分析。

基于以上分析，对研究内水体进行掩模处理后，选取了海拔、坡度、坡向以及太阳入射角和NDMI 5个因子，其中海拔取自然对数，对山区地表温度进行综合分析。统计量如表 3所示。

由偏相关系数一列可以看出，影响泰山山区地表温度的最大因素是地表自身的水汽特征，其次是海拔，且都与地表温度呈显著负相关关系。地表温度与其他3个地形因子间均是弱相关关系。

将表 3中的“非标准回归系数”列数据代入多元回归模型得到泰山山区的最优的逐步回归模型：

本文采用TM6单窗算法，获得泰山山区的地表温度。并应用统计学原理分析地表温度与地形因子和归一化指数(NDVI、NDMI)之间的关系，得到出以下结论：

(1)地表温度随着海拔的升高以显著自然对数形式递减。在地形因子对地表温度影响分析中发现，海拔是地表温度的最大影响因素，而坡度、坡向以及太阳入射能量等局地地形条件的影响均很小。

(2)反映下垫面特性的NDVI与NDMI对地表温度的影响分析发现地表温度与这两个指数都具有良好的负相关关系。但研究也发现NDVI在表征地表温度时，不但要对水体进行掩膜，还要对水体进行一定范围的缓冲区处理，另外地表温度在高植被覆盖区(NDVI > 0.5左右)会出现“饱和”现象。相比于NDVI与地表温度的二次项负相关，NDMI与地表温度的线性负相关更能有效表达地表温度。

基于以上研究可以发现影响夏季泰山山区地表温度的主要因素是下垫面特征，其次是海拔，而坡向、太阳入射角和坡度对地表温度也会产生不同程度的影响。本研究尝试综合考虑地形因素及下垫面特征分析地表温度空间分布差异性的影响因素，研究结果为科学的理解山区地表温度空间格局特征及机制提供一定的参考。整个研究区内的植被覆盖度较高，所得结论在裸地较多的山区是否适用还需单独分析。本文研究区内海拔差异较大，因此采用NDVI估算地表辐射率，会对地表温度的反演结果产生一定影响；除此之外，地表温度还会受人类活动、局地微气象等因素的影响，本研究都没有考虑，这有待于将来进一步研究与探讨。

致谢： 占文凤副教授对本文写作给予帮助，特此致谢。