近年来，由于城市化的快速发展，植被、农田等被建筑物、道路和其它不透水地表取代，影响了近地面水和热的自然交换，使城市地区通常有更高的太阳辐射吸收率、热容量和传导率，加剧了城市的热岛效应。国内外学者针对城市土地利用对地表温度(land surface temperature, LST)的影响进行了大量研究，取得了许多代表性成果^[1-8]。研究普遍认为，建设用地对地表具有较强的增温作用，而绿地和水体则具有显著的降温作用，如李文亮等^[9]利用Landsat影像进行哈尔滨市城市扩展对地表热环境的影响研究认为建设用地的扩张对地表热效应的贡献是最大的。孙俊等^[10]研究认为林地比农田和草地的调节温度效果好，在城市中增加林地的面积，对缓解城市热岛效应的作用最大。赵梓淇等^[11]，梁保平等^[12]研究则认为水体和林地在缓解城市热岛效应方面作用十分显著，两者比较，水体的降温效果优于林地。前人的这些研究无疑为解决城市热岛问题、土地合理布局提供详实的数据和可靠的依据。但这些研究大多从全局的角度出发探讨土地利用变化对城市地表温度的影响，将研究区域内地表温度受各土地利用类型的影响程度设定为常数，忽视土地利用影响的空间非平稳性，然而，城市地表温度的变化规律和过程十分复杂，其自身及其影响因子是与特定的空间位置相关联，具有强烈的空间非平稳性，随着空间位置变化，其主导因素及其影响程度可能发生变化^[13-15]。

地理权重回归(geographically weighted regression, GWR)模型是对普通最小二乘法回归(ordinary least squares regression, OLS)进行空间扩展而建立的局部空间回归分析方法，自1996年FORTHERINGHAM et al^[16]提出后，已被广泛应用于空间非平稳性关系探讨^[17-20]，但目前鲜少将该模型应用到土地利用与地表温度二者关系的研究中。本研究选择福建省福州市为研究对象，在360 m×360 m格网尺度下(栅格总数为6 210)，以地表温度为因变量，各土地利用比例为自变量，构建揭示土地利用与地表温度空间关系的GWR模型，从格网尺度上探讨城市土地利用变化对地表温度影响的空间非平稳性，为模型推广提供理论和技术支持。

1 研究区域与数据 1.1 研究区概况

福州是福建省的省会城市，地处我国东南沿海(东经118°24′~ 120°30′，北纬25°16′~ 26°39′)，属于亚热带海洋性季风气候，境内多山，濒临海洋，地势自东向西倾斜，地形起伏大。随着城镇化水平的不断提高，福州城市范围急剧扩大，从而导致了土地利用发生了显著变化，使得城市热岛效应日益严重，已成为福州可持续发展迫切需要解决的生态问题^[21]，本文研究区选择面积为795 km²的福州建成区及其周边区域。

1.2 数据来源及预处理

遥感数据主要包括Landsat 5 TM影像，获取时间为1996年8月21日和2006年8月17日; Landsat 8影像, 获取时间为2016年7月27日，行列号均为119/42。3期影像的成像时间均为10:30左右，研究区上空无云覆盖，数据质量良好，地物识别度较高。其它数据主要有用于辅助土地利用分类和精度检验的福州市第2次土地调查数据。遥感数据预处理包括辐射定标、大气校正、几何校正及影像裁剪等。

1.3 土地利用信息提取

参照中国土地资源分类系统的土地分类标准^[22]，并结合研究区的土地利用特征和遥感影像的可分辨程度，将研究区土地利用类型划分为耕地(水田和旱地)、林草地(有林地、灌木林地、疏林地、草地和其它林地)、建设用地(城镇用地、农村居民地用地、工矿交通用地)、水体(湖泊、河流、水库、池塘)和裸地(裸土地和沙滩) 5类。并采用人机交互式方式对分类数据进行校正和修改。

1.4 地表温度反演

目前Landsat 8第11波段(TIRS 11)值仍不稳定，故本文利用Landsat 8的第10波段和Landsat TM数据的第6波段进行地表温度反演，具体反演步骤如下^{[21, 23]}。

首先，通过公式T=K₂/ln(K₁/L+1)计算亮温。其中，T为地表亮温(K)；L为热红外波段辐射校正后的光谱辐射值(W · m^-2 · sr^-1 · μm^-1)；K₁和K₂为热红外波段的定标常数，TM6的K₁和K₂分别为607.76 W · m^-2 · sr^-1 · μm^-1和1 260.56 K，TIRS10的K₁和K₂分别为774.89 W · m^-2 · sr^-1 · μm^-1和1 321.08 K。

其次，采用覃志豪等^[24]提出的经验公式法进行地表比辐射率ε的估算，将地类分为水面、城镇和自然表面3类，计算各地类的植被覆盖度，进而计算各地类的地表比辐射率。

最后，根据公式$ {T_{{\rm{LS}}}} = \frac{T}{{1 + \lambda T/\rho {\rm{ln}}\varepsilon }}$将求出的地表亮温转换为地表温度。其中，T_LS为地表温度(℃)；λ为波段的中心波长(μm)，TM 6和TIRS 10的中心波长分别取值11.5和10.9 μm；ρ=hc/σ= 1.438×10^-2 mK，其中，σ=1.38×10^-23 J · K^-1 (斯特潘-玻尔兹曼常数)，h=6.626×10^-34 Js (普朗克常数)，c=2.998×10⁸ m · s^-1 (光速)；ε为地物比辐射率。

2 研究方法

GWR模型是在OLS模型[公式(1)]的基础上，假定其回归系数为变量空间位置的函数而建立的一种空间变参数模型[公式(2)]^[16]。

$ {y_i} = {\beta _0} + \sum\limits_{i = 1}^k {{\beta _k}{x_{ik}} + {\varepsilon _i}} $

(1)

$ {y_i} = {\beta _0}\left( {{\mu _i}, {\nu _i}} \right) + \sum\limits_{i = 1}^k {{\beta _k}\left( {{\mu _i}, {\nu _i}} \right){x_{ik}} + {\varepsilon _i}} $

(2)

式中：y_i为样本点i的因变量，β₀为回归的常数项，β_k为变量k的回归参数，x_ik为样本点i的第k个变量观测值，ε_i为误差项，(μ_i, ν_i)为样本点i的地理坐标，β₀ (μ_i, ν_i)为样本点i的回归常数项，β_k (μ_i, ν_i)为x_k在样本点处i的回归参数。

可采用加权最小二乘法实现GWR模型中β_k (μ_i, ν_i)的估算，矩阵形式表示如公式(3)：

$ {{\hat \beta }_k}\left( {{\mu _i}, {\nu _i}} \right) = {[{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{W}}\left( {{\mu _i}, {\nu _i}} \right)\mathit{\boldsymbol{X}}]^{ - 1}}{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{W}}\left( {{\mu _i}, {\nu _i}} \right)\mathit{\boldsymbol{Y}} $

(3)

式中：X为自变量矩阵，X^T为矩阵X的转置，Y为因变量矩阵，W (μ_i, ν_i)为空间权重矩阵，由回归点与其周围样本点之间地理距离的衰减函数，假设距离回归点i越近的样本点对局部回归参数的影响越大。空间权重矩阵是估算β_k (μ_i, ν_i)的关键，本文采用高斯函数来实现[公式(4)]。

$ {W_{ij}} = {\rm{ex}}{{\rm{p}}^{ - \frac{1}{2}\left( {\frac{{{d_{ij}}}}{h}} \right)2}} $

式中：W_ij为利用空间样点j估计待测点i时的权重，d_ij为样点j与待估测点i间的欧氏距离，h为带宽，利用最小赤池信息准则(Akaike information criterion，AIC)进行确定。

由于ArcGIS 10.2 GWR模型建立的门槛限制，如采用热红外波段的空间分辨率进行模型拟合时，格网内土地类型单一，土地利用比例出现较多的0值和1值，导致模型无法建立，故本文以120 m为间隔重采样，采用局部方差法^[25]，即平均局部方差最大时所对应的空间分辨率为研究的最佳格网尺度，最终选择的研究格网尺度为360 m×360 m。在ArcGIS 10.2软件支持下，统计360 m×360 m格网的平均地表温度和各土地利用比例，以地表温度为因变量，土地利用比例为自变量，分别建立OLS模型和GWR模型，定量揭示土地利用与地表温度的关系。并采用AIC、决定系数(R²)、残差及残差的空间自相关性等指标，将GWR模型拟合结果与OLS结果进行比较分析，探讨GWR模型在土地利用与地表温度空间关系的非平稳性研究中的有效性。

3 结果与分析 3.1 福州市的土地利用变化特征

分别在1996、2006、2016年的土地利用分类图上，每个地类各随机选取50个点(共250个)，以相应年份遥感影像的目视解译结果为参照数据，运用混淆矩阵进行精度检验。1996、2006、2016年土地利用分类图的总体分类精度分别为86.39%、91.44%和92.11%，Kappa系数分别为0.80、0.89和0.90，满足分类精度要求。对土地利用图进行统计，福州市土地利用及其变化特征如表 1所示。

由表 1可知，福州市土地利用以林草地、建设用地、耕地和水体为主，但随着时间推移，建设用地在研究期间呈现快速增长趋势，占全区比重由1996年的12.59%增加到2016年的36.32%，成为面积仅次于林草地的土地利用类型；而林草地和耕地面积则快速缩减，其中，耕地面积降幅最大，占全区比重由1996年的17.33%减少到2016年的5.96%；林草地面积占全区比重也从1996年的55.13%减少至2016年的46.10%。而水体的面积在研究期间略有下降，变化相对较小。

3.2 福州市不同土地利用类型的地表温度差异

图 1为福州市不同土地利用类型的地表温度统计特征。其中，2006年反演的地表温度总体偏低主要是由于影像获取时间是2006年8月17日，受8月中上旬(8月5日—8月11日)超强台风桑美的持续影响，相较于其它年份相同月份气温偏低，相应地地表温度总体偏低。1996、2006、2016年各土地利用类型的地表温度差异明显，并表现出较为一致的规律，平均地表温度最高的为裸地，其次为建设用地，接着依次为耕地、林草地和水体。其中，裸地和建设用地的平均温度较为接近，明显高于全区平均温度，耕地的平均温度略高于全区均温，而水体和林草地的平均温度则明显低于全区平均温度。

3.3 福州市土地利用对地表温度的影响 3.3.1 土地利用对地表温度影响的普通最小二乘法回归分析

土地利用比例与地表温度的普通最小二乘法回归拟合结果如表 2所示，所有方程均通过P=0.00水平的显著性检验，说明各土地利用类型对地表温度均有显著影响。其中，增加林草地和水体的比例能对研究区的地表温度起到很好的降温作用(方程拟合系数为负数)，随时间推移，林草地和水体的回归系数均呈现下降趋势，说明该两种地类的降温作用逐渐增强，尤其是水体，到了2006和2016年，其降温作用已超过了林草地，每增加10%的水体，温度分别降低0.33和0.49 ℃，超过了林草地的0.24和0.44 ℃；而增加建设用地、裸地和耕地的比例则对研究区地表温度有明显的增温作用(方程拟合系数为正数)，其中，以裸地的增温作用最为突出，每增加10%的裸地面积，地表温度可增加0.61~0.96 ℃ (1996、2006、2016年拟合回归系数分别为0.091、0.061和0.096)。建设用地的增温作用仅次于裸地，每增加10%的建设用地，地表温度可增加0.52~0.76 ℃。而耕地相较于建设用地与裸地，其增温作用较小，每增加10%的建设用地，地表温度可增加0.06~0.24 ℃，而且这种增温作用随着时间的推移呈现逐渐递减趋势，3a的回归系数分别为0.024、0.015和0.006。

从解释能力来看，建设用地对研究区地表温度的解释能力最强，1996、2006和2016年的R²均为同期所有地类中的最大值，其次为林草地和水体，而耕地对研究区的地表温度的解释能力则最小。随着时间推移，建设用地、林草地和水体的R²呈现持续增大趋势，尤以建设用地最为显著，其R²从1996年的0.513 5提高至2006年的0.725 9至2016年的0.832 1，而耕地和裸地的R²则呈波动式下降趋势。这说明了在研究期间，建设用地、林草地和水体对研究区地表温度的影响越来越大，越来越成为影响研究区地表温度的关键因素。故本文重点探讨建设用地、林草地和水体3种地类与地表温度关系的空间非平稳性。

3.3.2 土地利用对地表温度影响的空间非平稳性

建设用地、林草地和水体与地表温度空间关系的GWR拟合结果见表 3和图 2~图 4。

从表 3和图 2~图 4可知，GWR拟合结果总体上与OLS模型一致，建设用地与地表温度的拟合系数平均值大于0，且以大于0的区域为主，其对研究区以增温作用为主。而林草地和水体与地表温度的拟合系数平均值均小于0，且以小于0区域为主，该两种地类对研究区以降温作用为主，但该3种地类对地表温度的影响均存在显著的空间非平稳性。

建设用地比例的增加在北部林区的增温作用远大于中部建成区及耕地分布较为集中区域(表 3和图 2)，这种差异在1996年表现得尤为显著，建设用地比例每增加10%，在北部区域的地表温度将上升1.1 ℃以上，而在中部区域的地表温度上升则小于0.7 ℃。但随着时间的推移，回归系数平均值(1996、2006、2016年分别为0.098、0.071、0.068)和标准差(分别为0.171、0.060和0.028)均呈现持续下降的趋势，说明单位建设用地比例的增加对研究区地表温度的影响及其空间非平稳性均有所下降。

增加单位面积林草地的降温效果以建成区及其周边林区效果最好(表 3和图 3)，该区域回归系数小于-0.04，即每增加10%的林草地，可使地表温度下降0.4 ℃以上。以这些区域为中心向外，回归系数逐渐递增，降温作用逐渐减弱，到南部的闽江流域，回归系数甚至大于0，对水域则起增温作用；水体降温作用以福州城市中心区域效果最好(表 3和图 4)，该区域回归系数小于-0.06，每增加10%的水体表面可以使地表温度下降0.6 ℃以上，由市中心区域向外，水体的降温作用逐渐降低，并逐渐过渡到林区的增温作用。研究期间，林草地回归系数标准差呈波动式略有下降，其对研究区地表温度的影响空间非平稳性略有下降，而水体对地表温度的影响却随着时间推移其空间分异越来越明显，水体显著降温区(系数 < -0.06)的范围越来越大，到2016年已包含建成区及闽江流域的大范围区域，成为比重最大的区域，而增温区域(系数>0)集中分布于北部山区及南部的部分山区，比重仅次于显著降温区。

3.4 OLS和GWR模型对比分析

OLS模型和GWR模型均能有效地模拟不同土地利用与地表温度的关系，GWR模型是对土地利用与地表温度间关系的局部拟合，可以体现不同土地利用类型对地表温度作用的空间非平稳性。而OLS模型是对土地利用与地表温度间关系的全局拟合，揭示的仅是研究区土地利用与地表温度间关系的平均状态，在一定程度上存在高估或低估各类土地利用对地表温度的影响程度。在中部建成区，OLS模型高估了建设用地的增温作用(GWR模型系数以大于OLS模型系数区域为主)，却低估了林草地和水体的降温作用(GWR模型系数以小于OLS模型系数区域为主)；在北部林区，OLS模型低估了建设用地的增温作用和林草地的降温作用，但高估了水体对地表的降温作用；在水域，OLS模型低估了建设用地的增温作用和水体的降温作用，但高估了林草地的降温作用；而在部分林区和水域，OLS模型的拟合结果甚至与GWR模型揭示的在水域增加林草地或在林区增加水体会使地表温度升高的结果截然相反。

为了更好地评估GWR对土地利用和地表温度空间关系的预测能力，利用R²、AIC值、残差和残差的空间自相关性等指标对GWR和OLS的结果进行比较分析(表 4)。

从表 4可知，各地类GWR模型的AIC值远小于OLS模型(不同模型中的AIC值若能减少3个单位，认为模型有显著改善^{[20, 26]})，拟合残差较OLS模型也有显著下降，而调整R²则较OLS模型有显著提高，说明相较于OLS模型，各地类与地表温度的GWR模型均有显著改善，能更好地揭示建设用地、林草地以及水体与地表温度之间的关系。其中，水体的GWR模型较OLS模型提升效果尤为显著：GWR模型的解释度较OLS模型提高了65.3%以上，而1996、2006、2016年GWR模型的拟合残差较OLS模型则分别下降了1.158、0.916、1.696，为3种地类中的最大值。从模型拟合残差的Morans′I指数看，各地类的GWR模型较OLS模型均有较为明显的下降，但值仍较大，说明GWR模型能够消除OLS模型中出现的部分空间自相性，但存在一定的局限性，在今后的研究中，还需对模型进一步完善。

4 讨论与结论

本研究在360 m×360 m格网尺度下，以地表温度为因变量，各土地利用比例为自变量，运用GWR模型揭示土地利用与地表温度的空间关系，研究结果表明：福州市不同土地利用的地表温度差异明显，其中，裸地和建设用地的平均温度明显高于全区平均温度，耕地的平均温度略高于全区均温，而水体和林草地的平均温度则明显低于全区平均温度，这与戚鹏程等^[5]，青松等^[6]，王敏等^[8]，孙俊等^[10]研究结论基本一致。这主要是由于研究区不同土地利用类型的下垫面性质不同造成：裸地主要由沙滩、裸土、裸石以及待建的开发用地构成，其表明粗糙吸热快且热容量小，同时植被覆盖度低造成植被蒸腾作用弱，一旦吸收太阳辐射后地表温度上升较快，成为地表平均温度最高的土地利用类型；而水体则由于其比热容性质较大，在吸收同样的太阳辐射之后，温度升高较慢，成为地表平均温度最低的土地利用类型；建设用地的下垫面几乎为不透水面，且人为释放热量、绿色植被较少，其地表温度较高，平均温度仅略低于裸地；林草地的植被覆盖率较高，由于植被生理活动的作用，对太阳辐射有很强的吸收作用，植被的蒸腾作用较强，有利于地表和大气的热量交换，平均地表温度明显低于研究区；而耕地一方面由于表面较为粗糙，对太阳辐射有着较强的吸收作用，另一方面由于夏季植被覆盖较高，植被的蒸腾作用较强，使得其平均地表温度3 a分别为28.83、25.12和30.64 ℃，略高于同期研究区的平均地表温度。

从全局看，增加林草地和水体比例可明显降低研究区地表温度，而增加建设用地、裸地和耕地比例则会使研究区地表温度升高。其中，建设用地、林草地和水体是影响研究区地表温度变化的主要土地利用类型。该3种地类对地表温度的影响存在显著的空间非平稳性：增加建设用地比例使得北部林区地表温度的增加幅度远大于中部建成区及耕地分布集中区域；林草地的降温能力以建成区及其周边林区最强，以这些区域为中心向外降温作用逐渐减弱，到南部水域，林草地比例增加将使地表温度升高；水体的降温能力以建设用地密集分布的市中心最强，向外逐渐降低，过渡到林区的增温作用。

以上研究结果说明了地表温度与土地利用变化作为典型的生态过程，具有强烈的空间异质性，二者关系的变化规律和过程十分复杂，是与特定的空间位置相关联。如仅从全局的角度出发，忽视土地利用/覆盖作用的空间差异来探讨土地利用/覆盖对城市地表温度的影响，即将研究区域内城市热岛受各土地利用/覆盖类型的影响程度设定为常数，这与实际情况不相吻合，研究结果的通用性势必将受到影响。如孙俊等^[10]利用环境一号卫星的红外波段反演地表温度，揭示地表温度与地表类型的关系，研究认为林地比农田和草地的调节温度效果好，在城市中增加林地的面积，对缓解城市热岛效应的作用最大。而赵梓淇等^[11]，梁保平和李晓宁^[12]在城市LUCC时空格局对地表温度的影响研究中发现，水体和林地在缓解城市热岛效应方面作用十分显著，两者比较，水体的降温效果要优于林地。这些研究由于忽视土地利用/覆盖作用的空间异质性，揭示的是不同土地利用对地表温度作用的平均状态，有可能高估或低估某一地类的增温或降温作用，使得不同研究区域研究得出的结论可能是不同的。

GWR模型的应用前提是空间非平稳性，是对传统回归模型的扩展，将数据的地理位置嵌入到回归参数中，其结果是局部性的并非全局性的，从而能够探测到空间数据的空间非平稳性。本文运用GWR模型揭示土地利用与地表温度二者的空间关系，允许在不同的位置有不同空间关系的存在，充分体现了二者关系的空间非平稳性，使得GWR模型的调整R²较OLS模型提高显著，而AIC值和拟合残差则下降显著，GWR模型较OLS模型能更好地揭示建设用地、林草地以及水体与地表温度之间的空间关系，可为土地利用和地表温度的关系研究提供参考。