由卫星NOAA AVHRR数据得到的NDVI (nor-mal difference vegetation index)时间序列数据可用于长时间的土地覆被研究和多尺度的陆地生态系统模拟。NDVI与植被覆盖度、叶面积指数、生物量和生产力等性状关系密切^[1-2], 同时由于AVHRR数据的时间分辨率较高, AVHRR NDVI数据广泛应用于农作物种植面积估算、产量估算、作物长势监测以及森林草地等自然植被监测等大范围植被活动的研究^[3-4], 另外, 植被的生长过程与气温和降水等气候条件也有密切关系^{[1, 5]}。

20世纪90年代以来, 全球环境以前所未有的速度发生变化, 一系列全球性重大环境问题对人类的生存和发展构成严重威胁, NDVI动态变化及与气候因子的相互关系成为地理学和生态学的研究热点。相关学者进行了不同尺度和视角的研究^[4-6], 文献[7-10]探讨不同时间序列全国及局部生态敏感区(如黄河三角洲、京津冀风沙源区、青藏高原和黄土高原等地区) NDVI与气候因子的关系, 文献[11-13]分析了不同植被类型的NDVI对气候因子的响应机制。结果表明, 植被NDVI的变化特征及与气候因子的关系随着气候区、植被类型和研究时期的不同而有所差异^{[4, 6, 14]}。

以上研究主要基于统计回归和相关性分析说明NDVI与气候因子时间序列的相关关系, 并讨论不同区域、不同植被类型中相关性的差异, 即使在相似的覆被条件或植被类型中, 这种相关性仍具有空间非平稳性^[15-16]、非线性和时空异质依赖性。NDVI在中国的变化也有很强的空间异质性^[17]。例如, 在东部沿海地区NDVI有下降趋势, 在西部地区和一些农业区有上升趋势^[3]。传统的基于线性假设的研究方法与NDVI特性冲突, 难以发掘其动态变化与气候因子间的因果关系。近年来, 作为传统回归技术的延伸, 地理加权回归越来越多地应用于探索空间变化关系^[18]。

地理加权回归(geographically weighted regre-ssion, GWR)方法通过计算回归模型的局部参数, 可以很好地表明NDVI与气候因子关系的空间非平稳性和区域特性, 从而大幅度提高模型的拟合优度, 取得更好的拟合效果^{[16, 19]}。但是, 目前此方法多应用在局地区域研究^[20-23]。鉴于此, 本文以AVHRR NDVI和全国756个气象台站29年时间跨度的长序列数据源为基础, 分别运用GWR模型和OLS模型, 在国家尺度上拟合NDVI与气温和降水的关系, 并对结果进行对比, 旨在探索NDVI对气候因子的时空响应机制, 揭示其关系的空间非平稳性和尺度依存性。

1 数据来源与研究方法 1.1 气象数据

本研究使用的气象数据是从国家气象数据共享网(http://cdc.nmic.cn/home.do)获得的中国地面气候资料日值数据集, 选取中国756个气象台站的日均温和日降水量, 时间跨度为29年(1982年1月-2010年12月)。剔除部分长时间缺测站点后, 对各年的日均温求算术平均, 对日降水量求和, 得到604个站点(图 1) 29年的年均气温和年降水量。

1.2 AVHRR NDVI数据

本研究使用的遥感数据为AVHRR NDVI月值图, 来自NASA戈达德航天中心(Goddard Space Flight Center)全球库存建模和映射研究组(Global Inventory Monitoring and Modelling Studies Groups)的NDVI3g数据集。图像空间分辨率为8 km × 8 km, 时间分辨率为月, 时间跨度为1982年1月-2010年12月。为了减小误差, 以台站周围15 km内的NDVI均值作为台站的NDVI值。

1.3 研究方法

对于空间现象而言, 传统的统计分析方法大多忽视空间数据之间的相互依赖性, 因而其结果说服力不够^[24]。某些情况下, 空间上的依赖性比时间上的更复杂, 因此在传统统计方法基础上引入空间因素成为地理学研究方法的热点。GWR模型考虑了距离因素影响, 是研究地理空间异质性现象的回归模型。其特定位置的回归系数不是利用全局信息获得, 而是在考虑了临近空间要素影响的情况下, 利用临近观测值的子样本信息进行局部回归得到^[6]。GWR模型在普通线性回归模型的基础上, 将空间因素作为回归参数。因此, GWR模型可以视为普通线性回归模型的改进。

本研究将NDVI与气温和降水的关系分别用全局性的普通最小二乘法(ordinary least square, OLS)和地理加权回归(GWR)方法拟合, 并对结果进行对比分析。

1.3.1 经典的线性回归模型

经典的线性回归模型如下:

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

式中, y为因变量, 是自变量x_k(k=1, …, p)的一个线性组合; e_i为符合正态分布的随机误差项, i=1, …, n是观测值的个数。回归系数β_k(k=1, …, p)假定为常数, 一般采用OLS法估计。

1.3.2 地理加权回归模型

GWR模型对一般的线性模型进行拓展, 模型参数是位置i的函数^[5]。GWR模型表示如下:

$ {y_i} = {\beta _{i0}}\left( {{u_i}, {v_i}} \right) + \mathop \sum \limits_{k = 1}^p \;{\beta _{ik}}\left( {{u_i}, {v_i}} \right){x_{ik}} + {\varepsilon _i} $

式中, (u_i, v_i)为第i个样本点的空间位置, β_i0(u_i, v_i)和β_ik(u_i, v_i)x_ik分别是第i个样本点的常数估计值和参数估计值, e_i(i=1, 2, 3, …, n)是第i个样本点的随机误差项。若β_1k=β_2k=…=β_nk, 则地理加权回归模型就蜕变为普通线性回归模型。 由于位置i周围的观测值对i点的参数估计的影响随着距离的增加逐渐减小, 具有距离衰减效应, 因此利用距离加权OLS法估计参数^[25]:

$ \boldsymbol{\widehat \beta} ({u_i}, {v_i}) = {({\boldsymbol{X}^{\text{T}}}\boldsymbol{W}({u_i}, {v_i})X)^{-1}}{\boldsymbol{X}^{\text{T}}}\boldsymbol{W}({u_i}, {v_i})\boldsymbol{Y} $

式中, $\boldsymbol{\hat \beta} $为参数β的估计值, 是自变量观测值构成的矩阵, 其第一列元素为1; Y为由自变量观测值构成的变量; W为空间权值矩阵, 保证距离每一个样本点更近的观测值具有更大的权重。实证分析中经常用到的空间权值矩阵计算方法有高斯距离权值(Gaussian distance)、指数距离权值(exponential distance)和三次方距离权值(tricube distance)等, 本研究采用高斯距离权值方法确定W:

$ \boldsymbol{W}({u_i}, {v_i}) = {{\text{e}}^{ - \frac{1}{2}{{\left( {\frac{{{d_{ij}}}}{b}} \right)}^2}}} $

式中, W(u_i, v_i)为观测点j对第i个样本点的空间权重, b为核函数的基带宽度, d_ij为回归点i到观测点j的欧氏距离。 当因变量y和自变量x(x₁, x₂, …, x_n)没有空间差异时, 则称数据具有空间平稳性, 即参数β_k(u_i, v_i)不会随着空间位置的变化而变化。 Brunsdon等^[5]引入平稳性指数来评估地理加权回归模型中的数据平稳性, 计算公式如下:

$ {\text{SI}} = \frac{{{\beta _{{\text{GWR_iqr}}}}}}{{2 \times {\text{GLM_se}}}} $

式中, SI是平稳性指数, β_{GWR_iqr}是x_n系数的标准误的四分位距, GLM_se是基于全局性的普通最小二乘法回归的标准误。如果SI≤1, 则因变量y和自变量x(x₁, x₂, …, x_n)的关系具有平稳性。 对地理加权回归模型进行模型显著性检验和模型误差项空间自相关的显著性检验, 常用的检验方法包括蒙特卡罗检验^[5]、AIC信息准则检验^[26]和Leung等^[14]的检验方法。AIC检验用于对比不同模型的显著性, 较低的AIC值表明模型的模拟效果更好。表达式如下:

${\text{AIC}} = 2n\ln (\hat \sigma ) + n\ln (2{\text{π }}) + n\left[ {\frac{{n + {\text{tr}}(S)}}{{n - 2 - {\text{tr}}(S)}}} \right]$

其中, n是样本的大小, $\hat \sigma $是误差项估计的标准离差, tr (S)是地理加权回归模型的投影矩阵S的迹, 是带宽的函数。当两个模型对的AIC值之差大于3, 则具有较小AIC值的模型对数据的模拟效果更好。

2 研究结果 2.1 NDVI与气候因子关系的带宽选择

尽管GWR模型可以研究因变量与自变量之间关系的空间分布, 然而核函数的带宽选择会对回归结果的空间表现有重要影响。带宽过小, 估计值接近真实值, 偏差大; 带宽过大, 估计值接近普通最小二乘回归的结果, 偏差大。本文采用高斯距离权值法确定权值矩阵, 结合自变量的平稳性SI指数和模型的AIC值确定最优带宽, 并用局部加权最小二乘法估计参数。

从图 2可以看出, NDVI与气温的平稳性指数SI在空间范围350 km时下降到1.03, 说明NDVI与气温的关系在带宽350 km以上具有平稳性。NDVI与降水的平稳性指数SI在空间范围450 km时下降到1.09, 说明NDVI与降水的关系在带宽450 km以上时具有平稳性。如果在全国范围对NDVI与气温和降水建立空间平稳性的回归关系, 带宽应该在350~450 km之间。根据AIC法则, 本文选取350 km为最优带宽。

2.2 NDVI与气温和降水关系的空间异质性

整体来看, 基于GWR模型的NDVI与气温的回归系数(图 3)比NDVI与降水的回归系数(图 4)大, 表明气温对NDVI的影响大于降水。这一结果与已有的研究结果基本上一致, 如李本纲等^[27]认为, 对中国大部分地区而言, 气温对植被的影响超过降水。NDVI与气温和降水的回归系数及相关系数空间格局分析表明, 存在着从北到南和从东到西的空间变化。整体上表现为西北部高、东南部低。这与前人的研究结论类似。如李晓兵等^[28]对1983-1992年中国主要植被类型NDVI变化与气温和降水的关系的分析显示。从北到南, NDVI与气候因子的相关性逐渐降低; 从东南到西北, NDVI与气候因子的相关性逐渐增加。

中国地域广阔, 自然环境条件具有显著的空间差异, 因而NDVI与气温和降水的关系也呈现显著的空间异质性(图 3和4)。NDVI与气温的回归关系在全国大部分地区均为正值, 表明气温对于NDVI有正向影响, 即温度越高, NDVI值越大。在华中、华北、西南、西北和青藏高原地区, 这种正向的关系表现得尤为显著。然而, 在东北和华北地区, NDVI与气温的回归系数出现负值, 表明气温在这些地区对NDVI产生负向控制作用。

NDVI与降水的回归系数呈现更加复杂的空间格局, 在全国大部分地区均为正值, 但数值较小, 表明降水对NDVI具有弱的正向影响。在东北西部、西北和东南沿海地区有较弱的正向关系, 而在华北、华中、华南和西南地区, 两者的正向关系更弱甚至转变负向关系, 这与其他研究成果基本上一致^{[24, 29]}。

GWR模型的局部拟合优度(调整的R²)的空间分布如图 5所示。可以看出受人类活动干扰较小的中西部和东北部等以森林、草原和荒漠为主要覆被类型的地区, NDVI主要受气温和降水等气候因素的控制, 模型的拟合优度较高。东南地区的植被除受气候等自然因素作用外, 人类活动的扰动更强, NDVI受到更多复杂因素的影响, 因此模型的拟合优度较小。

2.3 OLS与GWR模型拟合效果对比

本研究分别应用OLS模型和GWR模型对NDVI进行拟合, 并对拟合效果进行对比, 结果如表 1所示。基于OLS模型回归的拟合优度为0.33, 调整的拟合优度为0.32, AIC值为5537.55;基于GWR模型回归的拟合优度为0.65, 调整的拟合优度为0.59, AIC值为5142.42。GWR模型拟合优度显著大于OLS模型, 因此GWR模型的拟合效果更好。GWR模型的AIC值小于OLS模型, 且两者之差大于3, 可以判定地理加权回归模型比普通线性回归模型更接近真实状态, 即NDVI与气候因子之间存在明显的空间非平稳性^[5]。

在空间回归分析中, 为了保证回归模型参数的可靠性, 空间样本点之间不应具有显著相关性。通常回归残差是否随机分布可以作为检验样本点空间自相关的标准。GWR模型回归残差的空间分布(图 6)比OLS模型更离散。应用最近邻法对GWR模型回归残差进行自相关检验, 残差满足随机分布。

2.4 NDVI与气候因子关系的区域差异

由于中国国土面积较大, 因此自然环境条件复杂多变。利用地理加权回归的最优带宽、拟合优度和AIC值, 可以分析NDVI与气候因子关系的空间异质性。根据温度、降水、海拔、植被和地形地貌等数据, 将中国分为4个生态地理区, 即东北部湿润半湿润生态地理区、北部干旱半干旱生态地理区、南部湿润生态地理区和青藏高原生态地理区^[30]。对不同生态地理区的样本点分别计算平稳性指数(图 7), 可以得到以下结论。

1) 各个生态地理区NDVI与气温和降水平稳性关系的空间尺度不同。在北部干旱半干旱生态地理区, NDVI与气温之间关系达到平稳性的空间尺度是225 km, 降水则是725 km (图 7(a)); 在东北部湿润半湿润生态地理区, NDVI与气温之间平稳性关系的最小空间尺度是234 km, 而降水则是242 km (图 7(b)); 在南部湿润生态地理区, NDVI与气温之间的关系在130 km以上达到平稳性状态, 降水则为93 km (图 7(c)); 在青藏高原生态地理区, NDVI与气温之间的关系在109 km以上达到稳定, 而降水则为152 km (图 7(d))。整体上, 除南部湿润生态地理区外, 气温与NDVI达到平稳关系的空间尺度小于降水与NDVI之间的关系; 随着带宽尺度减小, 数据平稳性下降, 模型回归结果可信性小; 随着带宽尺度增加, 数据平稳性增加, 模型回归结果可信性大。一般情况下, 自然环境控制因素越单调, NDVI与气候因素达到平稳性的空间尺度就越小, 反之就越大。因而, 各区NDVI与气候因子平稳性关系的空间尺度是区域面积、地形复杂度、人类扰动强度等多种因素综合作用的结果。

2) 各区GWR模拟性能有较大差异。根据各个生态地理区NDVI与气温和降水平稳性关系空间尺度的分析结果, 进行分区地理加权回归计算, 得到结果如表 2所示。可以看出, 不同生态地理区的地理加权回归结果有显著差异。GWR模型在青藏高原生态区拟合优度最高, 达到0.98; AIC值较小, 为577.58;在东北部湿润半湿润生态区拟合优度最低, 为0.39; AIC值最大, 为1194.77。

3 结论与讨论 3.1 结论

本文应用地理加权回归模型(GWR), 以1982-2010年中国多年平均NDVI为因变量, 年降水量和年均气温为自变量, 对NDVI与气候因子关系的空间非平稳性和区域特征进行分析, 结论如下。

1)  NDVI与年均温和年降水量之间的关系分别在350 km和450 km尺度上趋于平稳。根据AIC法则和平稳性指数的变化趋势, 确定350 km为最优带宽, 在这一尺度上建立的NDVI与气候因子的回归关系较为可靠。

2) 与一般全局性回归模型比较, 不论是模拟精度还是空间结构, GWR模拟结果都显著优于一般全局性回归模型。GWR模型的R²为65, 调整的R²为0.59, 而OLS方法的R²为0.33, 调整的R²为0.32; GWR模型的AIC值为5142.42, 而OLS模型的AIC值为5537.55。

3) 中国区域广阔, 自然环境空间异质性很高, 故NDVI与气候因子的关系存在显著的区域差异, 表征为:各区NDVI与主导气候因子发生作用的特征尺度大小不同; 模拟各区NDVI与主导气候因子地理加权回归模型的性能不尽相同。

3.2 讨论

本文研究结果表明, NDVI与气候因子之间的关系存在显著的时空差异, 并且具有多尺度特征。这一特性已经逐步被研究者所认识, 如刘绿柳等^[31]通过对黄河流域NDVI与降水、温度的年际变化趋势及相互关系的时空变化规律分析, 发现该流域NDVI与降水、温度相关显著的植被类型以草地、灌木为主, 但相关区域的空间位置随时间变化。郭铌等^[11]的研究也证实, 西北地区NDVI与气温的相关系数大于降水, 天山、阿尔泰山和秦岭的NDVI与气温相关系数最高, 而青海东北部NDVI与降水的相关系数最高。

本研究也存在一些不足。首先, GWR模型虽然对于刻画NDVI与气候因子的空间非平稳性和尺度依存特性具有独特的优势, 但自然界一些环境因子对NDVI的作用范围是全局性的, 应用GWR时会出现一些虚假的回归空间结构。将来的工作可以寻求混合型地理加权回归模型来解决这一问题。其次, 本文仅分析了NDVI与气候因子的回归关系, 实际上, 除气候因子外, 诸多环境因素都影响到NDVI, 如地形和土壤条件等。下一步工作将把这些因子融合到地理加权回归模型中, 以期得到更为客观的分析结果。