全球变化与陆地生态系统的关系是全球变化科学研究中的核心内容之一^[1]。植被作为陆地生态系统的重要组成部分, 不仅在大气圈、生物圈和全球碳循环中扮演着重要的角色, 也是连接土壤、水分和能量等自然要素的纽带, 同时对全球变化具有重要指示作用^[2—3]。探究植被绿度时空变化规律及其对气候变化的响应, 可为全球变化研究提供科学依据。其中, 叶面积指数(Leaf area index, LAI)对植被绿度具有良好的表征作用, 它指单位土地面积上植物叶片面积占土地面积的倍数, 能够反映出水平覆盖状况和垂直结构, 且LAI与植被光合、蒸腾和碳氮循环等生态过程关系密切^[4—5]。因此LAI评估对于全球变化背景下物质循环与能量流动研究具有重要意义。

影响植被绿度的气候要素主要是气温和降水^[6—9], 植被绿度对气候变化的响应具有一定的时滞效应与非线性特征^[10—15], 且影响过程与机理研究已相对成熟。越来越多的研究表明地表温度和土壤湿度在植被生长及格局变化中起着关键作用：一方面, 地表温度的升降影响土壤和大气水分交换；另一方面, 土壤湿度变化直接影响土壤对植被生长用水的补给^[16]。相较于气温和降水, 地表温度和土壤湿度作用于植物根部, 对植物种子萌发与植株生长的作用更加直接^[17]。植被绿度受到地表水热综合作用, 这种作用是正向的, 还是负向的？是即时的, 还是滞后的？是线性的, 还是非线性的？这些问题尚不十分清楚。

秦岭-大巴山是中国北亚热带向暖温带过渡的重要地段, 兼具两种气候带自然地理特征, 地表水热差异悬殊, 加之地形起伏较大, 随着地形的变化, 气候也会发生明显变化^[18], 加剧了地表水热条件的不确定性和复杂性。且以往对于地表水热要素的研究主要集中在样方尺度和点尺度, 在区域尺度植被对地表水热的响应过程有待进一步研究^[17]。因此, 在全球变化大背景下, 以地形地貌特征为基础, 探讨地表水热如何影响山地植被变化过程及未来空间格局, 成为全球变化和山地地理共同关注的问题^[7]。

本文采用MODIS-LAI数据集, 结合DEM、坡向、植被覆盖、地表温度数据, 以及基于温度植被干旱指数法反演的土壤湿度数据, 运用均值、趋势分析和Hurst指数分析了秦巴山地植被绿度空间格局与变化趋势, 并探讨了植被绿度与地表水热在不同地形上的空间关系, 借助相关性分析、主导因素分析、时滞性及广义可加模型(Generalized additive models, GAM), 进一步探究了LAI对地表水热的时滞及非线性响应过程等。本研究利于加深对秦巴山地气候和植被变化机理的认识, 对了解南北过渡带植被对全球变化的响应以及防范山地生态风险具有一定价值。

1 数据处理与研究方法 1.1 研究区概况

秦巴山地地跨甘肃、四川、陕西、重庆、河南、湖北六省市, 东西长约1000 km, 南北宽约300 km, 总面积约为30万km², 是北亚热带和暖温带的分界线(图 1)。地理位置的特殊性使得研究区植被区系复杂, 秦岭南坡主要为混有常绿树种的落叶阔叶林, 北坡则主要为暖温带落叶阔叶林, 具有明显过渡特点。研究区年平均气温为12—16℃, 年平均降水量为709—1500 mm, 由于受大陆性气候和季风气候的双重影响, 该地区气候对全球变化较为敏感。秦巴山地也是生物物种东西向和南北向迁移的通道, 物种分布具有明显的多样性、复杂性和过渡性^[19], 该地区在气候、植被、土壤等诸多自然地理要素方面都表现出垂直地带性、经度地带性、纬度地带性、坡向分异等特点。

1.2 数据来源及处理

本研究使用的LAI数据为MOD15A2H产品, 其空间分辨率为500 m, 时间分辨率为8 d(https://ladsweb.modaps.eosdis.nasa.gov/search/), 研究时段跨度为2000—2020年, 研究区涉及h26v05和h27v05两个影像分带。研究表明该数据集在全球主要的六套LAI产品的交叉验证中显示不确定性最低^[20], 且其时空分辨率较其它植被指数产品更具优势。数据预处理上, 首先使用MRT(MODIS Reprojection Tool)工具进行拼接和提取单波段处理, 然后进行最大值合成(Maximum Value Composition, MVC)以最大程度消除云遮挡^[21], 随后在ArcGIS中进行裁剪、单位换算、重采样、投影、格式转换等处理, 最终生成空间分辨率为1000 m、投影为WGS1984 Albers的逐月LAI数据。秦巴山地植被在冬季大多处于休眠状态, 并且有时被积雪覆盖, 这一时段LAI值对植被动态指示意义不大, 参照以往有关秦岭地区植被的研究^[22], 本文选取每年4—10月份作为植被生长季, 其均值作为年LAI数据。

LST数据为2000—2020年MOD11A2产品, 空间分辨率为1000 m, 时间分辨率为8 d(https://ladsweb.modaps.eosdis.nasa.gov/search/)。在ArcGIS中将其处理成与LAI数据具有相同分辨率与投影的2000—2020年逐月LST数据。归一化植被指数(Normalized difference vegetation index, NDVI)数据为2000—2020年MOD13A2产品, 空间分辨率为1000 m, 时间分辨率为16 d(https://ladsweb.modaps.eosdis.nasa.gov/search/)。NDVI数据结合LST数据用来反演温度植被干旱指数TVDI, 作为土壤湿度的衡量指标, 具体反演方法见1.3.1。为与LAI数据相对应, 计算4—10月份LST及TVDI数据均值, 分别作为LST和TVDI的年值数据。

使用2000—2018年气象站点地表温度日值数据集作为MODIS地表温度的验证数据, 站点数据来源于国家气象科学数据中心(http://data.cma.cn/)。将秦巴山地75个气象站点的多点月均数据与MODIS月均地表温度数据绘制散点图(图 2), 使用曲线拟合优度R²和均方根误差(Root Mean Square Error, RMSE)表征两组数据的差异程度, R²越高, RMSE越低, 表示两组数据差异越小, 即MODIS地表温度数据的精度越高。结果表明, R²为0.9442, RMSE为2.56, 说明MODIS数据能很好地反映秦巴山地地温状况, 可用于后续分析。RMSE计算方法见相关文献^[23]。

植被类型数据为2001和2019年MCD12Q1产品, 空间分辨率为500 m, 时间分辨率为1 a(https://ladsweb.modaps.eosdis.nasa.gov/search/)。该数据基于国际地圈生物圈计划(IGBP)17类全球植被分类方案, 包括11个自然植被类型、3个土地开发和镶嵌的地类和3个非草木土地类型定义类^[24]。数据经过严格的质量控制, 对地表覆被变化反应敏感且空间连续性较好, 适用于较大空间尺度的植被类型变化监测^[25]。在ArcGIS中筛选出两期数据中未变化部分作为本文植被类型数据, 此方法能够有效消除单期数据的偶然性和不确定性。相较于自然植被, 作物对地表水热的响应受人为干扰较大, 因此在后续分析时剔除作物及作物和自然植被体的镶嵌体, 保留落叶阔叶林、常绿和落叶混交林(以下简称混交林)、多树草原、稀树草原、草原5种植被类型(图 3)。植被类型数据用来分析地表水热对不同植被类型LAI的影响程度。

参考赵芳等基于秦岭、大巴山主山脊线和主要河流(嘉陵江、汉江)位置的划分方法, 将秦巴山地按坡向分为秦岭南坡、秦岭北坡、大巴山南坡、大巴山北坡、西秦岭^[26], 坡向数据用于分区统计秦巴山地不同坡向LAI与地表水热的分布特征。DEM数据来源于地理空间数据云平台STRM DEM 90 m分辨率高程数据(http://www.gsclod.cn)。为反映秦巴山地植被绿度与地表水热的垂直分布规律, 本文将DEM以500 m为间隔划分海拔梯度：< 500 m、500—1000 m、1000—1500 m、1500—2000 m、2000—2500 m、2500—3000 m、3000—3500 m、3500—4000 m及>4000 m, 由于4000 m以上面积较小, 因此将>4000 m区域分为同一海拔梯度。

1.3 研究方法 1.3.1 基于NDVI-LST特征空间的TVDI模型构建

Sandholt^[27]等基于植被指数和地表温度之间的二维特征空间模型提出了温度植被干旱指数(Temperature vegetation drought index, TVDI), 它是一种基于光学与热红外遥感通道数据进行植被覆盖区域表层土壤水分反演的方法^[28]。其公式为：

式中, Ts为地表温度；Ts_max为NDVI对应的最高地表温度, 即干边；Ts_min为NDVI对应的最低地表温度, 即湿边；a₁, b₁, a₂, b₂分别为干、湿边线性拟合方程的系数。

本研究筛选相同NDVI下研究区内所有像元对应的最高和最低地表温度, 得到252个(每年12个, 共21年)时间段的LST-NDVI特征空间。TVDI取值范围为[0, 1], 湿边TVDI值为0, 干边TVDI值为1。土壤湿度与TVDI值呈负相关关系, TVDI越小土壤湿度越大, TVDI越大土壤湿度越小。以第230期为例(图 4), 干边拟合R²=0.5937, 湿边拟合R²=0.6530, 模型拟合效果良好, 可用于后续分析。

1.3.2 秦巴山地LAI变化趋势分析

Slope趋势分析可反映某一时段植被变化趋势, S>0代表植被呈改善趋势, S < 0代表植被呈退化趋势^[29]。Hurst指数是定量描述时间序列长程依赖性的主要方法之一, 可以反映现在与未来变化的一致性^[30], 0.5 < H < 1, 代表时间序列是一个持续性序列, 即未来变化趋势与过去趋势保持一致；如果H=0.5, 代表时间序列为随机序列, 即未来变化趋势与过去无关；如果0 < H < 0.5, 代表时间序列是反持续性序列, 即未来变化趋势与过去趋势相反。两者结合可反映秦巴山地LAI现在及未来变化趋势, 将秦巴山地植被趋势分为5类：持续改善区(S>0, 0.5 < H < 1)、由增变减区(S>0, 0 < H < 0.5)、由减变增区(S < 0, 0 < H < 0.5)、持续退化区(S < 0, 0.5 < H < 1)和随机变化区(H=0.5)。Slope和Hurst的计算方法见相关参考文献^[29—30]。

1.3.3 秦巴山地LAI影响因素分析

(1) 相关性分析。通过像元尺度的皮尔森相关性分析, 探究地表水热对LAI的影响。参考相关文献^[10], 设置当月(即生长季同期)、提前1个月(3—9月份)、提前2个月(2—8月份)3种地表水热情景, 分别与生长季(4—10月份)LAI计算相关系数, 探讨LAI对地表水热是否存在时间滞后响应。

(2) 主导因子分析。以LAI为因变量, LST、TVDI为自变量建立像元尺度的多元线性回归方程, 然后对回归系数进行标准化, 标准化回归系数绝对值最大的因素即为该像元影响LAI的主导因素, 计算方法见相关参考文献^[17]。

(3) 广义可加模型(Generalized additive models, GAM)。传统的归因分析假设自变量与因变量为线性关系, 然而实际中因变量往往不随自变量而单调变化, 而是呈现非线性相关。GAM模型是由数据驱动而非统计分布模型驱动的非参数回归模型, 可对部分解释变量进行线性拟合, 对其他因子进行光滑函数拟合。模型不需要预先设定参数模型, 模型通过解释变量的平滑函数建立, 能够自动选择合适的多项式, 适用于响应变量与解释变量之间的关系是非线性或非单调的数据分析^[31—32]。本研究基于Rstudio中的mgcv包, 探讨响应变量LAI与解释变量LST和TVDI的非线性关系。

2 结果与分析 2.1 秦巴山地植被绿度时空分布及其与地表水热的关系 2.1.1 植被绿度空间分布及其与地表水热的关系

LAI、LST以及TVDI年均值在不同海拔、坡向上具有空间差异(图 5)。秦巴山地LAI值范围在0—7, 其中LAI值在1—5的区域面积占比94.67%。LST范围在2.4—32.9℃, TVDI范围在0.13—0.95。从不同海拔分析, 随海拔升高LAI均值先增后减, < 1500 m区域LAI随海拔升高而增加, ＞1500 m区域LAI随海拔升高而减少, LST及TVDI均随海拔升高呈减少趋势。< 500 m地区LAI值较低, LST和TVDI值较高；500—3000 m地区LAI值相对较高, 其中1000—1500 m地区LAI值达到最大, 该区域LST和TVDI值适中；>3000 m地区LAI、LST和TVDI值均较低。从不同坡向分析, LAI均值大小为：大巴山北坡>大巴山南坡>秦岭南坡>秦岭北坡>西秦岭。大巴山为LAI高值区, 该区域的LST值相对较高, TVDI值相对较低；西秦岭为LAI低值区, 对应LST和TVDI的低值区。因此, LAI的分布格局并不是单一地表水热要素作用的结果, 而是水热组合的综合调控的产物。LAI高值往往出现在地表水热条件较好(LST较高/适中, TVDI适中/较低)的地区。

2.1.2 植被绿度变化趋势及其与地表水热的关系

秦巴山地LAI由增变减区占比最大(72.99%), 其次为持续改善区(21.14%), 持续退化、由减变增和随机变化区累计占比不到6%(图 6)。表示目前秦巴山地大部分区域LAI处于增加趋势, 但未来变化的反向特征强于同向特征, 以减少为主。不同海拔带LAI变化趋势存在明显差异(图 6)：< 2500 m及3500—4000 m地区LAI以由增变减为主, 持续改善次之；2500—3500 m及>4000 m地区以由减变增为主。从不同坡向分析(图 6), 各坡向LAI均以由增变减最多, 持续改善次之, 两者加和超过90%, 说明LAI变化趋势的坡向分异不明显。

为阐释秦巴山地LAI趋势与地表水热趋势的关联性, 统计了LAI趋势面积占比最大区域(由增变减区)的同期地表水热变化情况(图 6)。结果表明, LST以由增变减(56.04%)和持续增加(28.18%)为主, TVDI以由减变增(59.25%)和持续减少(35.23%)为主。表明当前大部分地区的地表水热呈现增加趋势, 但未来变化的反向特征强于同向特征(以减少为主), 这与LAI趋势具有较好的一致性。但小部分区域的地表水热呈现持续增加趋势, 即增加趋势在未来具有可持续性, 这与LAI未来的趋势相反, 说明该区域未来地表水热不断增加, 但LAI却有所减少。

2.2 秦巴山地LAI与地表水热的相关性及其主导因子

在像元尺度上计算了LAI与地表水热的相关性系数, 并进行显著性检验。结果表明LAI与同期LST(当月LST)以正相关为主(图 7), 正相关通过0.05显著性检验的区域占研究区面积的88.38%, 正相关均值为0.52, 两者呈负相关面积较少, 零星分布在整个研究区。LAI与提前1个月、提前2个月的LST主要呈正相关(图 7), 正相关均值均为0.28, 通过0.05显著性检验的面积占比分别为47.52%、35.20%。说明LST对秦巴山地LAI影响的时滞性不明显。

LAI与同期TVDI(当月TVDI)以负相关为主(图 7), 通过0.05显著性的负相关像元占比45.45%, 负相关均值为-0.33。LAI与提前1个月TVDI的负相关像元占比44.03%(P < 0.05), 负相关均值为-0.28, 低于与当月TVDI的相关性(图 7)。LAI与提前2个月TVDI相关系数则以正值为主(10.87%, P < 0.05), 但面积较小, 正相关均值为0.19, 负相关占比2.73%(图 7)。因此, TVDI与LAI以负相关为主, 即土壤湿度越高, 植被生长状况越好, LAI与同期TVDI相关性最好, 强于与提前1个月与提前2个月的相关性, 表明TVDI对LST影响的时滞效应不明显。

地表水热对秦巴山地绝大部分地区LAI均有较大影响, 相比之下, TVDI对LAI影响范围更加广泛, 主导的区域占比72.22%, LST主导区域仅占25.49%(图 8)。地表水热对不同植被类型的主导面积有一定差异(图 8)。TVDI在落叶阔叶林、混交林、多树草原、稀树草原和草原中主导面积的比例分别为78.87%、57.34%、87.09%、74.90%、49.15%。LST在落叶阔叶林、混交林、多树草原、稀树草原和草原中主导面积比例分别为21.13%、42.66%、12.91%、25.10%、50.85%。表明TVDI对落叶阔叶林、混交林、多树草原和稀树草原的LAI占主导地位, TVDI和LST对草原LAI的主导面积接近。

2.3 LAI对地表水热的非线性响应

通过对地表水热与LAI间建立GAM模型, 获得解释变量的平滑回归函数, 得到影响因素对LAI影响效应图(图 9)。GAM模型检验结果如表 1所示, 各植被类型模型调整判定系数均在0.7以上, 方差解释率均超过70%, 说明建模结果良好, 模型具有很好的解释效果。自由度接近1时, 表示是线性关系, 自由度越大时, 则表示曲线关系越强。另外, 为消除自变量共线性的影响, 计算了LST与TVDI的相关系数, 结果表明相关系数绝对值在各植被类型中均小于0.4, 自变量之间不存在共线性, 建模结果可靠。

LST和TVDI对各植被类型LAI变化影响显著(P < 0.001)。(表 1、图 9), 结果表明, 在落叶阔叶林地区, LAI与LST呈线性正向关系(自由度=1), LAI随LST增加而单调递增。LAI与TVDI呈非线性负向关系(自由度=2.670)：TVDI < 0.45时, 随TVDI增加LAI变化平缓, TVDI>0.45时, LAI对TVDI响应的敏感性增强, 随TVDI增加迅速降低。说明落叶阔叶林LAI对TVDI的响应存在一定阈值, 使得LAI的响应规律在阈值两侧呈现明显差异。在混交林地区, LAI与LST大致呈线性正相关(自由度=1.845), 与TVDI呈非线性负相关(自由度=2.750), TVDI < 0.35时, 随TVDI增加LAI变化平缓, 有轻微增加趋势, TVDI>0.35时, 随TVDI增加LAI迅速降低。在多树草原地区, LAI对LST和TVDI具有非线性的响应特征, 当LST和TVDI达到一定阈值, LAI变化不再明显, LST的阈值在25℃左右, TVDI的阈值在0.6左右。在稀树草原地区, LAI对LST和TVDI具有非线性的响应特征, 当LST>27℃, TVDI < 0.7时, LAI变化较为平缓。在草原地区, LAI对LST和TVDI为非线性响应, 随LST增加, 曲线斜率增大, 说明LST对LAI的影响增强, 而TVDI与LAI为波动负相关(自由度=5.731), TVDI在0.4左右出现拐点。

综上所述, 除落叶阔叶林和混交林LAI随LST增加而线性递增外, 各植被类型LAI对地表水热均为非线性响应。LAI随LST增加主要呈非线性增加趋势, 随TVDI增加呈非线性降低趋势, 存在一定阈值界限, 使得阈值两侧LAI的响应规律呈现明显差异, 且不同植被类型的阈值不同。

3 讨论

秦巴山地LAI具有海拔和坡向分异, 一方面是因为地形不同直接导致植被类型的差异^[33], 另一方面, 地形分异导致了地表水热条件的再分配, 从而对植被绿度起到间接的调控作用。LAI高值分布于地表水热较好的1000—1500 m中等海拔段以及大巴山地区。其中1000—1500 m海拔段植被类型以混交林、落叶阔叶林和多树草原为主, 且地表水热状况较好, 无极端地表水热条件制约, 植被覆盖度较高^[34], 该区域为LAI高值中心。另外, 大巴山作为南北过渡带的南部, 植被条件更接近亚热带地区, 优于位于北部的秦岭^[35], 且大巴山的地表温度和土壤湿度高于秦岭地区, 因此该地区LAI值也相对较高。LAI低值区分布在 < 500 m海拔段、>3000 m海拔段以及西秦岭地区, 造成这种分布的原因有所不同：< 500 m的低海拔区植被受人类活动影响较大, 城镇化水平较高, 建筑及道路等基础设施建设占用了大量生态用地^[34]；此外该海拔地区TVDI值较高, 即土壤湿度较低, 受水分胁迫, 该地区植被生长状况较差。而西秦岭和>3000 m高海拔地区以草原为主, 植被条件较差, 虽然土壤湿度条件较好, 但地表温度较低, 极大程度上限制了植被生长, 这与陈超男等基于NDVI的结果较为一致^[34]。因此, 地表温度和土壤湿度均影响植被生长, 但两者交互叠加对驱动植被动态响应具有更重要意义^[36]。

秦巴山地LAI以由增变减为主, 未来变化的反向特征强于同向特征, 这与刘宪峰等的研究结果一致^[29]。在LAI的由增变减区, 大部分区域地表水热也呈由增变减趋势, 两者的变化趋势具有较好的一致性, 小部分区域地表水热呈现持续增加趋势, 这表明一定程度的增温和增湿有利于植被活动增强, 但超过一定幅度则可能会对植被产生不利影响^[14]。

植被对气温和降水往往存在时滞效应^{[10, 37—38]}, 而本文结果表明LAI对地表水热的时滞效应不明显, 与同期地表水热的相关性更强。表明相比气温和降水, 地表水热对植被生长的作用更具直接性。秦巴山地LAI受土壤湿度的主导范围大于地表温度, 在蒙古高原地区也有同样的结论^[39], 而在青藏高原等高寒地区, 地温则为植被生长的主导因素^[17]。造成这种差异的原因是高海拔地区的地表温度较低, 对植被生长的限制作用尤为突出, 在本研究中西部高海拔地区地温主导面积明显上升, 也印证了这一点。

研究发现植被动态对气候具有非线性响应且存在响应阈值, 在阈值两侧响应规律不同^[14], 该结论已在天津地区^[40]、长江流域上游^[41]和海河上游^[36]等地的研究中得到证实。本文结果表明, 秦巴山地各植被类型LAI对地表水热同样以非线性响应为主。首先, 气温和地温适当升高可以增加净光合速率、蒸腾速率和气孔导度^[42], 但增温过高可能会加速土壤水分蒸发, 形成干旱趋势, 植被还可能通过减小叶面积和降低光饱和点等方式防止自身失水, 对植被生长的促进作用消失, 甚至限制植被生长^[43]。其次, 土壤中许多养分和矿物质要在溶于水的条件下才能被植物吸收, 因此适宜的土壤水分对植被生长起促进作用^[44], 然而, 当水分的增加超过植被所需时, 会影响根系呼吸, 对植被的生长正向作用减弱^[45]。

本文主要探讨了秦巴山地LAI对地表水热的响应。但随着人类活动的加剧, 剥离出人类活动和气候因素对植被的影响成为未来需要解决的问题。此外, 植被变化对地表水热的反馈作用也是今后研究的主要内容。

4 结论

本文首先分析了秦巴山地LAI均值和趋势在海拔、坡向等尺度的空间分布格局, 探讨了LAI与地表水热在空间格局上的关联性, 然后分析了LAI与不同时滞情景下地表水热的相关性以及影响LAI变化的主导因素, 最后进一步探究了不同植被类型LAI对地表水热的非线性响应情况。主要得到以下几点结论：

(1) 秦巴山地LAI集中在1—5, 在海拔分布上, LAI均值随海拔先增后减, 最大值出现在1000—1500 m地区, 在坡向分布上, 大巴山南北坡>秦岭南北坡>西秦岭。LAI高值区与地表水热条件较好的区域具有较好的空间一致性, LAI低值区主要分布在LST较低或TVDI较高的地区。

(2) 秦巴山地LAI的变化趋势以由增变减为主。除2500—3500 m及4000 m以上地区以由减变增为主外, 其他海拔均以由增变减为主, 持续改善次之, 坡向对LAI趋势的影响差异不明显。在LAI由增变减区, 地表水热均以由增变减为主, 持续增加次之。

(3) 秦巴山地LAI与LST为显著正相关, 与TVDI为显著负相关。LAI对地表水热以同期响应最为显著, 时滞效应不明显。在主导因素上, 相较于LST, TVDI对LAI影响范围更加广。在混交林、落叶阔叶林、多树草原和稀树草原地区, TVDI主导面积广, 在草原区, LST与TVDI主导面积接近。

(3) 非线性分析结果显示, 除落叶阔叶林和混交林LAI与LST呈线性正向关系外, 各植被类型LAI与地表水热均为非线性关系。随LST增加, LAI主要呈现非线性增加的特征, 随TVDI增加, LAI则呈非线性降低趋势。存在一定阈值, 使得阈值两侧LAI对地表水热的响应存在差别, 且响应阈值因植被类型不同而存在差别。