植被类型对树轮重建贺兰山NDVI的影响

李明明，李刚

(1.中国科学院西北生态环境资源研究院，沙漠与沙漠化重点实验室，甘肃 兰州 730000；2，中国科学院大学，北京 100049；3.甘肃祁连山国家级自然保护管理局，甘肃 张掖 734000)

树木作为组成森林中独立的生命有机体，对气候变化极为敏感且为减缓气候变化带来的影响发挥了重要作用，在其生长过程中树木年轮记录了大量的气候、环境等方面的信息[1-2].然而，树木的生长又分为初级生长和次级生长[3-4]，初级生长包括发芽、枝叶的生长及根系的衍生，次级生长包括树干、枝及根系的生长.一方面树木冠层的叶片通过光合作用合成有机质并通过树木的韧皮部输送到树干以及树根的各部位，另一方面，树干将根系吸收的水和无机盐等营养物质通过边材运送到树冠，并与根系一起支撑整棵树木[5].因此，研究植被冠层与主干间的生长关系，有利于我们明确气候变化背景下树木生长的动态变化，更好地了解树木的生理机制以及为应对气候变化提供参考依据.

随着遥感技术的进步与发展，利用遥感卫星监测探究植被的地表覆盖度、生物量、净初级生产力变化的研究日益增多[6-8]，也为研究气候变化对植被动态的影响提供重要参考.归一化植被指数(normalized difference vegetation index，NDVI)作为使用较普遍的遥感数据之一，它是植物生长状态及植被空间分布密度的最佳指示因子，能反映冠层繁茂程度和植被覆盖度的动态变化[9].但NDVI仅代表树木冠层的变化情况，而生长过程中树木冠层与主干之间存在差异，且由于遥感卫星数据记录较短，限制了我们对过去植被变化历史的认识与理解.利用树轮资料对NDVI进行重建则可有效解决这一问题，因为树轮资料具有定年准确、分辨率高、样本量大等优点[10]，被广泛用于重建过去不同时期的温度、降水、干湿指数、径流量等环境信息[11-15]，同样可以利用树轮资料重建NDVI，不仅可以弥补NDVI时间序列长度上的不足，且对揭示年代至百年尺度的植被生长变化规律、预估未来气候情景下植被生长变化具有重要意义.

国内外学者已经开展了NDVI与树轮资料的相关研究[16-22]，汪青春等[16]利用树轮资料重建了青海布哈河流域的NDVI，王文志等[17]利用树轮资料重建了祁连山地区的NDVI.此外，Kaufmann等[18]发现在北半球中高纬度地区两者的相关性受不同纬度和树种影响；D’Arrigo等[19]研究发现美国北部地区树轮最大晚材密度与NDVI呈显著相关，然而在北美高纬度地区、欧洲地区和加拿大的部分地区并没有发现显著相关性[20-22]，表明树轮资料和NDVI之间的响应关系可能因研究区域和树种或植被类型不同有所差异.以上研究在不同空间尺度上证实了NDVI与树轮资料存在相关性，但结果存在差异.NDVI重建过程有诸多不确定性，可能是因为NDVI包含区域内所有植被的生长信息，不同植被类型的生理过程迥异，而现有的研究中几乎均未考虑NDVI中植被类型对重建结果的影响[16-22]，因此直接将NDVI结果与树轮资料结合进行重建，可能导致重建的NDVI序列出现偏差.如何提高重建结果的可靠性，使重建的NDVI更加准确反映植被变化是一个值得探讨的问题.本研究选取贺兰山为研究区域，尝试将NDVI中林地与草地信息分离，并单独使用林地NDVI与树木径向生长的关系对该地区NDVI进行重建，以期更准确地了解贺兰山地区林地NDVI的变化历史，为预测未来森林动态提供理论.

1 材料与方法

1.1 试验区概况

贺兰山(N 38°27′～39°30′、E 105°20′～106°41′)位于宁夏回族自治区与内蒙古自治区交接地带(图1)，地处东亚夏季风的西北缘，是我国干旱与半干旱、草原与荒漠的过渡地带，生态环境敏感、脆弱[23].贺兰山地区为典型的温带大陆性气候，优势植被主要为针叶林，分布在海拔2 000 m以上，2 000 m以下主要为草地[24].根据气象站观测资料显示(图3)，贺兰山全年干旱少雨，多年平均气温为9.2 ℃，最冷月(1月)平均气温为-8 ℃；最热月(7月)平均气温为23.8 ℃；年降水量在250～438 mm，主要集中在5～9月，约占全年降水量的80%左右，年蒸发量在2 000 mm以上[25].由于年降水量与蒸发量相差较大，空气干燥；日照强烈，年平均日照达3 000 h以上，无霜期短.

图1 贺兰山土地覆被及树轮采样点Figure 1 Land cover and tree-ring sampling site in Helan Mountains

1.2 样本采集与年表建立

本研究油松样本采集自贺兰山苏峪口国家森林公园，采样点选取受人为干扰较小且远离旅游区的原生林，海拔在2 100～2 300 m.共选取98棵油松样树，分别用生长锥从垂直和平行山坡的方向上在树高1.3 m处采集2根树芯.对采集的样芯按照国际通用的试验处理步骤进行样品前处理[27]：将样芯晾干、固定、打磨至显微镜下树木年轮轮廓清晰可见，并剔除断裂及破损的部分样芯，共筛选得到183根样芯用于本研究.使用LINTAB-6系统对树木年轮宽度进行测量，精度为0.01 mm.之后利用TSAPwin程序完成交叉定年[28]，通过COFECHA程序对交叉定年结果进行检验.对完成定年后的树木年轮宽度数据使用ARSTAN软件来去除与树木年龄有关的生长趋势[29]，并基于负指数函数对数据进行拟合，得到树轮宽度标准年表(图2).选取1920～2018年(EPS>0.85)的贺兰山油松树轮宽度年表用于后续分析.

图2 贺兰山地区树轮宽度标准年表及树芯样本量Figure 2 Tree-ring standard chronology for Helan Mountains and sample size expressed as number of cores

1.3 气象数据

气象数据下载自中国气象科学数据共享服务网(http：//cdc.nmic.cn /)，所选台站为银川气象站，其距离采样点50 km，可用时段1951～2018年.选取该站的平均气温、平均最高气温、平均最低气温、降水量、相对湿度等要素用于本研究.

帕尔默干旱指数(palmer drought severity index，PDSI)用于分析温度和降水的综合影响，数据下载自荷兰皇家气象研究所(KNMI)(http://climexp.knmi.nl)，格点位置：N 38°42′，E 105°45′，时段1934～2018年.

1.4 土地覆被数据

本文选用清华大学研发的全球土地覆盖产品[26](FROM-GLC10)，下载网址：http://data.ess.tsinghua.edu.cn/fromglc2017v1.html.该数据通过对覆盖全球的陆地卫星专题制图仪和增强型专题制图仪传感器(Landsat TM/ETM+)影像进行分类，获得10 m分辨率的FROM-GLC10，该数据共分为10个不同土地覆盖类型，包括：1 耕地，2 林地，3 草地，4 灌木丛，5 湿地，6 水体，7 苔原，8 建设用地，9 裸地，10 冰川与永久积雪.在ARCGIS10.2中将FROM-GLC 2017数据与NDVI影像叠加，将不同植被类型进行分类，裁剪出贺兰山地区林地和草地的NDVI影像数据.

1.5 NDVI数据

NDVI数据下载自中国科学院资源科学数据中心(http://www.resdc.cn)，该数据是基于SPOT-VEGETATION NDVI卫星遥感数据合成产品，空间分辨率为1 km，时段为1998～2018年.该数据通过最大值合成法(MVC)得到每月1期的月值数据，进一步消除了云和卫星偏角等噪音影响[30].该数据有效反映了在空间和时间尺度上的植被覆盖分布和变化状况，对植被变化状况监测、植被资源合理利用和其它生态环境相关领域的研究有十分重要的参考意义[31].将分类提取后林地、草地的NDVI影像数据转化得到NDVI时间序列[30].

2 结果与分析

2.1 贺兰山NDVI的变化特征

1998～2018年贺兰山5～9月林地NDVI的平均值为0.277，高于草地(0.119)以及未分类情况(0.17)，且林地NDVI变化幅度略高于草地NDVI(图3-A).对比21年NDVI月均值发现(图3-B)，1～3月NDVI维持在全年较低水平，4月下旬植被缓慢生长NDVI开始上升，5月植被生长旺盛NDVI上升速度较快，5～8月持续增长，Gao等[32]2017年在贺兰山地区微树芯监测油松形成层细胞分生速率的结果显示(图3-B)，油松在6～7月形成层细胞分生速率最快.8月温度开始下降(图3-C)，但8月的降水量全年最高(41.22 mm)，此时水热条件下植被仍可生长，NDVI在8月达到峰值，但不同植被类型的NDVI峰值差异较大，林地为0.30，草地为0.13.随着9月降水开始下降，温度进一步降低，植被冠层光合作用强度降低，生长速率开始减缓趋于停滞，植被逐渐停止生长(图3-B).

图3 贺兰山地区不同植被NDVI年际(A)、年内(B)变化趋势，形成层细胞分生速率及银川气象站温度和降水变化(C)Figure 3 Variations of NDVI of different vegetation types on annual (A) and intra-annual (B) scales and cambium cell growth rate in Helan Mountains；distribution of monthly mean temperature and monthly total precipitation of Yinchuan meteorological station (C)

2.2 不同植被类型NDVI与树轮宽度指数的关系

1998～2018年不同植被类型NDVI与树轮宽度指数的相关性分析结果显示(图4)，上一年7月至当年3月，NDVI与树轮宽度指数呈负相关且不显著；在生长季，未分类NDVI与树轮宽度指数仅在6月和7月呈显著正相关，相关系数分别为0.484和0.528(n=21，P<0.05)；草地NDVI与树轮年表在6月和7月也呈显著正相关(r=0.561，r=0.603，n=21，P<0.01)，且草地的相关系数高于未分类.草地NDVI与树轮宽度指数在6月和7月呈显著正相关，相关系数分别为0.561和0.603(n= 21，P<0.01)，林地NDVI与树轮宽度指数在5～7月均达到0.01显著性水平，7月相关性最高，相关系数为0.685(n=21，P<0.01).林地NDVI与树轮宽度指数的相关系数在6～7月高于草地和未分类的情况.不同月份组合的相关性分析结果也表明，林地NDVI(6～8月)与树轮宽度的相关系数最大(r=0.645，n=21，P<0.01)，5～9月的月份组合两者同样具有较高的统计相关性(r=0.634，n=21，P<0.01).NDVI与树轮宽度指数的显著相关性说明两者对限制因子的响应结果较一致，树轮宽度指数可以反映生长季特别是夏季(6～8月)林地NDVI的变化.

图4 1998～2018年贺兰山地区不同植被类型NDVI与树轮年表的相关系数Figure 4 Correlation coefficients between NDVI of different vegetation types and tree-ring chronology in Helan Mountains during 1998～2018

2.3 树轮宽度与不同植被类型NDVI对气候因子的响应

为进一步探究树木径向生长对气候因子的响应，本研究分析了树轮年表、林地NDVI、草地NDVI与温度、降水以及PDSI的关系(图5-A).结果表明，在1951～2018年树木径向生长与温度呈负相关，与降水呈正相关，显著相关的月份较少，两者在6月相关性最强；与PDSI的相关性高于对温度及降水的响应，几乎在整个水文年均呈现显著正相关，与5月和6月份相关系数最高，分别为0.56和0.566(n= 68，P<0.05).因为PDSI是降水、温度及土壤水分等多种因素的综合指标[33]，是水分供需累积效应的近似度量值，可以更好地指示土壤水分变化，表征土壤对于树木生长的可供水量.

1998～2018年6～8月林地NDVI和温度、降水的相关性结果显示(图5-B)，林地NDVI与降水主要是正相关关系，而与温度主要为负相关关系，6月NDVI与降水的相关系数为0.452(n=21，P<0.05).林地NDVI同样表现出对降水响应的滞后效应.草地NDVI对温度、降水的响应结果(图5-C)，逐月相关性不显著，且相关系数低于林地NDVI与气候因子的相关性.因此，林地NDVI与温度、降水及PDSI的相关性高于草地，表明林地对气候变化可能更敏感.

图5 树轮年表(A)与气候因子、林地NDVI(B)和草地NDVI(C)与气候因子的相关性分析Figure 5 Correlation coefficients between tree-ring chronology (A)，forestland NDVI (B) and grassland NDVI (C) and climate variables

2.4 贺兰山夏季林地NDVI的重建

树轮年表与贺兰山地区夏季林地NDVI的相关系数为0.645(n=21，P<0.01).基于树轮年表重建贺兰山地区1920～2018年的林地NDVI，利用最小二乘法建立线性回归方程：

NDVI=0.662×STD+0.221

(1)

公式(1)中NDVI表示林地NDVI在6～8月的重建值，STD表示贺兰山树轮标准化年表.根据拟合方程，本文重建了贺兰山1920～2018年林地NDVI的变化(图6-A).为检验重建的可靠性，利用逐一剔除法对拟合方程进行稳定性检验，误差缩减值(RE)为0.357(P<0.01，n=99)，说明重建方程是稳定可靠的.

根据1920～2018年贺兰山林地夏季NDVI重建结果(图6-A)，林地生长状况最好的5 a分别是：1950，1956，1963，1970和1979；最差的5 a分别是：1928，1947，1966，1973和1982年.对NDVI进行10 a低通滤波分析结果显示，1950～1979年NDVI高于平均值(0.284)，森林生长状况较好；1922～1933年低于平均值，森林生长状况较差.将重建的林地NDVI结果(图6-A)与银川气象站点6～8月的温度、降水等器测数据进行对比(图6-B～C)，发现1951～1981年该地区的温度低于6～8月的多年平均温度22.5 ℃，降水较平均值37.6 mm/月略高，而NDVI处于整个序列平均值之上，树木生长状况较好.根据重建结果，贺兰山地区植被生长较差的时段是1922～1933年，通过旱涝指数发现(图6-E)，造成植被生长较差的主要原因是持续的干旱.在1920～1930年左右，该地区发生过极端干旱事件，主要是降水的减少导致植被生长状况较差，该阶段NDVI处于低值水平.此外，1982年左右贺兰山地区再次出现极端干旱事件，NDVI出现极端低值年，但NDVI对该极端干旱事件存在一定的滞后性.为进一步探究重建的林地NDVI与气象因子的关系(表1)，对原始数据进行一阶差分去除趋势噪音的影响，发现重建的林地NDVI和降水、PDSI的相关性提高.重建NDVI序列与降水的相关系数为0.586(n=69，P<0.01)；与PDSI的相关系数为0.508(n=84，P<'0.01),可见重建结果是可靠的.

图6 1920～2018年贺兰山夏季林地NDVI变化历史(A)与温度(b)、降水(C)、PDSI(D)、旱涝指数(E)对比Figure 6 Comparison between reconstructed summer NDVI (A) and temperature (B)，precipitation (C)，PDSI (D)，drought and flood index (E) in the Helan Mountains

表1 1920～2018年NDVI与温度、降水、PDSI、旱涝指数的Pearson相关系数

3 讨论

通过NDVI变化特征可以了解贺兰山地区植被的生长动态历史.贺兰山地区不同时段NDVI的增长速率相异，但不同植被类型的年内增长趋势大体一致，与温度、降水基本维持同步变化.因为树木开始生长需要特定的条件，比如当积温和降水到达一定程度后树木形成层才开始活动[34-35].贺兰山地区树木径向生长的开始时间大致为4月下旬或者5月初，局部受地形的影响可能在5月之前就已经开始生长，初期生长较为缓慢，这也是5月之前NDVI值较低的主要原因.将NDVI中的林地类型提取后发现其NDVI值高于草地和未分类，可能是由于贺兰山地区建群种为常绿针叶树种，一年四季并未完全枯黄落叶，而草地多为一年生草本植物，其经历发芽、生长、枯黄等过程，导致其NDVI值不及林地.尤其是贺兰山地区林地冠层的变化，即使该地区针叶林的冠层常绿，但随着生长季前期温度升高和降水增多，进入生长季树木冠层绿度逐渐发生变化，这一变化遥感卫星可以较好的捕获.另外，2017年在贺兰山利用微树芯对形成层细胞分生速率监测结果显示[32]，细胞分生速率最快的时期在6～7月，而NDVI在8月达到最大值，这可能是因为在树木生长初期，仅有一部分物质和能量用于树木径向生长，还有一部分用于冠层中新生枝叶的生长[36]，冠层密度和绿度逐渐增加，NDVI也逐渐上升.

此外，NDVI是植物叶面由红光和近红外2个波段合成计算得到的植被指数[7]，主要反映植物叶片的绿度，与光合作用密切相关.树轮宽度反映树木径向生长状况，主要取决于光合作用与呼吸作用所产生的部分碳水化合物[34]，另一部分碳水化合物用于根、茎、叶生长发育以及呼吸作用等.森林冠层生长与径向生长之间存在联系，例如当干旱发生时，树木径向生长与冠层生长减退，表现出较好的一致性.贺兰山地区植被生长最差的时段是1922～1933年，造成植被生长下降的主要原因是降水的减少和持续的干旱，该极端干旱事件已被较多研究结果所证实[37-39].此外，1982年左右贺兰山地区出现了极端干旱事件，但NDVI对干旱的响应存在滞后性，体现了树木对极端干旱的“遗留效应”[40].陈峰等[33]在贺兰山地区发现油松生长主要受水分的限制，我们通过对比分析多气候序列指标，同样发现该地区植被生长受干旱条件的控制，降水是影响树木径向生长的主要原因.

贺兰山地区植被类型主要为林地和草地，生长季在5～9月，分类后林地NDVI值高于草地NDVI，且林地NDVI和树轮年表的相关系数高于草地NDVI与树轮年表的相关系数.两者NDVI的差异可能与植被类型、叶片形状、叶绿素和水分的含量等因素有关[41-42]，叶绿素对蓝色(470 nm)和红色(670 nm)波段最为敏感，叶绿素作为主要的吸光物质，直接影响到植物光合作用的光能利用效率，含量越低，蓝、红波段吸收减弱，可见光波段反射率升高近红外反射率减弱，反之叶绿素含量越多，蓝、红波段吸收增强，可见光波段反射率降低，近红外反射率增强[35]，且近乎所有的近红外辐射由于被散射掉而很少被吸收.杨红飞等[43]发现不同植被类型反射光谱曲线形态和特征不同，造成林地和草地冠层光谱反射差异较大.贺兰山地区树木在夏季冠层茂密，地面生物量大，相对叶绿素含量较高；该地区草地稀疏，地面生物量低，叶绿素含量和水分含量较低.李明明在贺兰山地区通过提取NDVI中不同植被冠层物候与油松形成层物候进行分析后发现[44]，NDVI中林地、草地的物候参数存在差异，林地的植被结束生长时间受8月最低温影响显著，草地则不显著.这也说明了不同植被类型对温度的不同敏感度，也说明不同植被NDVI所包含的植被信息可能存在差异.因此不同植被的生理过程存在差异，例如林地和草地的根系深浅不同，且对环境因子的变化敏感度也存在差异[45].将林地NDVI单独提取后，树轮年表与林地的NDVI相关性有所上升，揭示了考虑不同植被类型对重建NDVI的必要性.

4 结论

根据贺兰山1998～2018年间林地NDVI数据与树轮年表的分析结果，验证了植被类型对NDVI重建产生影响.对比分析草地、林地以及未分类的NDVI与树轮年表的相关性结果表明，分类后的林地NDVI能更好的表示树木生长信号.根据林地NDVI与树轮年表重建1920～2018年6～8月林地NDVI的变化序列，显示在过去的98 a树木生长状况波动较大，受20世纪20至30年代持续干旱事件的影响，该地区林地生长状况较差；20世纪40至80年代降水增加树木生长状况较好；20世纪80年代至今生长基本保持稳定状态.这种树木生长状况的变化受局地的水分状况影响，极端干旱事件的“遗留效应”影响树木径向生长.本研究证实了不同植被类型会影响NDVI中树木生长信号的反映，并为基于树木年轮重建NDVI的研究提供参考依据.