吉林省落叶松林净初级生产力时空特征及其对气候变化的响应

李 芸,王轶夫,*,孙玉军,雷渊才,邵卫才,李 杰

1 北京林业大学森林资源和环境管理国家林业和草原局重点实验室, 北京 100083 2 中国林业科学研究院资源信息研究所,北京 100091 3 北京市通州区园林绿化局,北京 101121 4 北京市十三陵林场,北京 102200

净初级生产力(Net Primary Productivity,NPP)是绿色植物呼吸消耗后剩下的单位面积单位时间内所固定的能量或所生产的有机物质,作为表征陆地生态过程的关键变量,能反映植物在自然状态下的生产能力,是全球碳循环过程中的重要组成部分[1—4]。由于地形地貌和生物参数的复杂性,准确估算大规模植被NPP一直较为困难,往往会采用生物量转换成NPP作为实测值。目前,基于林木测树因子构建的异速生长模型求算森林生物量和生产力被广泛应用[5]。在建立模型时,通常将林木各组分(树叶、树枝、树干、树根)生物量与实际测树因子(如胸径D)分别拟合各组分生物量模型,并采用似乎不相关回归建立总量及各组分可加性生物量模型,以解决各组分生物量模型间不相容性或不可加性[5]。随着遥感技术的发展,高分辨率遥感卫星影像估算NPP的方法被广泛应用,Thakur等[6]利用印度遥感卫星(IRS-ID)LISS III传感器数据,选取归一化植被指数等为自变量,建立了NPP光谱模型,用于估算热带森林NPP和生物量；Ji等[7]利用5个主要森林生态系统的NPP观测数据、遥感数据建立基于统计学的多元回归模型对中国森林NPP进行估算,并分析了2000—2018年中国森林净初级生产力变化及其驱动因素。遥感数据往往只能估算较短时间的NPP,相对而言过程模型的优势在于能够从机理上模拟植被生命过程及生物地球化学循环过程,已被国内外广泛应用于NPP的模拟[8]。崔博超等[9]利用CASA模型并结合MODIS数据对2006—2016年塔里木河流域草地净初级生产力时空分异特征进行了研究；孙静等[10]综合考虑InTEC模型中干扰和非干扰因子,并嵌入林火数据,重新模拟了塔河地区不同地位指数下的森林NPP及其与林龄的关系。Jiang等[11]将土壤热模型(STM)整合到LPJ-DGVM模型中改进了地表3 m内的土壤温度变化后对北极圈内碳储量和生产力的模拟；孙国栋[12]利用LPJ-DGVM模型对1981—1998 年中国区域潜在植被分布和碳通量进行模拟；Sallaba等[1]利用LPJ-GUESS模拟了欧洲潜在自然植被当前和未来气候变化和CO2情景下的NPP。LPJ-DGVM全球动态植被模型(Lund-Potsdam-Jena Dynamic Global Vegetation Model,下文简称“LPJ模型”)作为目前被广泛应用的过程模型,已经成为模拟大尺度(区域至全球)的植被地理分布、生产力和碳平衡以及预测气候变化对陆地生态系统潜在影响的有效工具[13—16]。

吉林省位于全球变暖较为显著的中高纬度带,是中国重要的林业省份之一,在生态环境建设、碳汇管理乃至全球碳循环研究中占有重要地位[17—20]。落叶松(Larix)林是吉林省的代表性森林类型之一,也是我国重要的用材林。《中国森林资源报告(2014—2018)》显示,吉林省落叶松林面积和蓄积分别占全省森林总面积和总蓄积的8.20%和5.60%。近年来,我国对落叶松林的生长、结构和生物量等开展了大量的调查和研究工作[21—23]。然而,关于吉林省落叶松林生产力动态变化规律,目前尚缺乏全面系统的研究,本文旨在利用LPJ模型模拟吉林省落叶松林的生产力现状,分析其时空格局,为吉林省落叶松林生产力和植被监测研究提供参考依据。

1 数据与方法

1.1 研究区概况

吉林省(121.63°—131.32°E,40.87°—46.30°N,图1)位于我国东北地区中部,地处北温带,气候多样、地貌复杂,月均气温在-15—22 ℃之间,年降水量为400—800 mm,自东向西呈明显的由湿润到半湿润再到半干旱气候的变化趋势。境内主要山脉为大黑山、张广才岭、吉林哈达岭等,以中部大黑山为界将全省分为中西部平原和东部山地,具有明显的东南高、西北低的地势。森林集中分布于中、东部山区和低山丘陵区,乔木树种以落叶松(Larix)、红松(Pinuskoraiensis)、云杉(Piceaasperata)、樟子松(Pinussylvestrisvar.mongolica)等针叶树种和蒙古栎(Quercusmongolica)、白桦(Betulaplatyphylla)、杨树(Popular)等阔叶树种为主,灌木主要有粉枝柳(Salixrorida)、绣线菊(Spiraeasalicifolia)、忍冬(Lonicerajaponica)、榛子(Corylusheterophylla)等,草本以禾本科(Gramineae)、菊科(Compositae)为主。

图1 研究区气象站点及筛选落叶松样地分布Fig.1 Distribution of weather stations and Larix sampling plots in study area

1.2 气候数据

气象数据来源于中国气象数据网(http://data.cma.cn),包括2000—2019年吉林省及其周边100 km范围内的共41个气象站点逐月数据,采用Kriging插值法将离散的站点数据插值为连续的栅格数据。

图2 研究区2000—2019年平均气温和年降水量的变化趋势Fig.2 Trends of mean temperature and annual precipitation during 2000—2019 in study area

1.3 NPP遥感观测数据

NPP遥感数据来源于2000—2019年的MOD17A3H数据产品(https://e4ftl01.cr.usgs.gov/),空间分辨率为500 m,单位为kg C/m2(下文统一换算为g C/m2)。

1.4 固定样地实测数据

落叶松林固定样地实测数据来自于吉林省第八次和第九次国家森林资源连续清查数据(下文简称“连清”数据),调查时间为2009年和2014年,分布于幼龄林(25个)、中龄林(37个)、近熟林(14个)和成熟林(20个)共96个样地。根据Dong等[24]和Wang[25]研究中的中国东北地区主要树种生物量方程计算样地内样木生物量,汇总并除以样地面积得到样地单位面积生物量。按照样地前后两时期的单位面积生物量之差乘以含碳率(统一采用0.5[26—27])再除以间隔期计算NPP年平均值。

1.5 LPJ-DGVM全球动态植被模型

LPJ模型是模拟陆地生态系统的过程模型,从植被动力学原理出发,以植物光合生化反应、冠层能量平衡、生物量异速生长分配以及土壤水平衡等为基础,在同一个模型框架中逐栅格逐年地模拟生态系统植物的光合作用、呼吸作用、蒸发散等快反应过程与资源竞争、组织周转、异速分配、凋落物分解等慢反应过程,并且在模拟种群的建立与个体死亡过程中考虑了生态系统的自然干扰因素和自然死亡影响[4,15—16,28]。这些过程主要受模拟区域的环境条件、物种组成以及物种属性的影响。根据物候、生理、形态特征不同,LPJ模型可以模拟8 种木本植物功能型(Plant Function Type,PFT)和2种草本植被功能类型。为了实现模拟单元中植被个体到种群的尺度转换,模型以PFT平均个体为计算单元,每个模拟单元由一个或多个PFT组成,木本植物碳储量按尺度转换规则分配到4个植被组织库(叶、边材、心材和细根)以及地面凋落物库、地下凋落物库和两个土壤碳库中。植被光合作用是在冠层尺度上进行计算的,并按照植被光合作用积累的有机物质(即总初级生产力)扣除自身呼吸消耗量计算NPP(Sitch等[15]),计算公式如下:

NPP=GPP-Rm

(1)

式中,GPP为总初级生产力,Rm为植被自身呼吸消耗。

模型由气候数据(包括月均温、月降水量、月湿润天数、月日照百分率)驱动,同时需要输入模拟区域的土壤属性数据、大气CO2浓度数据。气象数据的获取与处理见1.2,土壤属性数据来自国际粮农组织(Food and Agriculture Organization, FAO)的土壤数据集。由于模拟时间区间较短,研究区CO2浓度是利用夏威夷Mauna Loa站1951年至今观测的全球平均大气CO2浓度作为该研究区历史平均浓度。植被分布最敏感的气候因素为最冷月月均温,只有实际月均温在生长温度范围内并且达到生长所需的积温下限后该物种才会出现,根据《植被图》分类系统的二级单位(植被型),原LPJ模型将其定义为北方落叶针叶林与实际分布不符,本研究将落叶松林划分为温带和寒温带落叶针叶林,并对其生物气候限制参数进行修正(见附表1)。

模型从裸地(无植被覆盖)开始模拟,经过1000年的“spin-up”过程使植被覆盖和土壤碳库达到平衡态,本研究使用2000—2019年20年的气候数据驱动这一过程。考虑到研究区的落叶松林均为天然次生林和人工林,在1001年时对模型进行了“归零”处理,即对4个植被组织库(叶、边材、心材和细根)的生物量归零,凋落物碳库和土壤碳库均照常模拟,最后分龄组获取落叶松林2000—2019年NPP进行分析。

1.6 精度评价指标

采用相关系数(R)、平均相对误差(MRE)和平均绝对误差(MAE)对LPJ模型模拟值与实测数据、遥感数据分别在NPP年均值和样地尺度上的年均NPP进行精度验证,其中高优指标R越大越好,低优指标MRE、MAE越小越好,计算公式如下:

(2)

(3)

(4)

1.7 变异系数

变异系数(Coefficient of Variation,CV)能够反映地理要素随时间的波动程度[29],可以用来评估2001—2019年吉林省落叶松林NPP在时间序列上的稳定性(式(5—6)):

(5)

(6)

1.8 变化趋势分析

本文采用相关性分析法探究温度、降水与吉林省落叶松林NPP的相关性；采用线性回归分析,研究吉林省落叶松林NPP和气候因子的年际变化趋势[30—31],回归方程斜率(β)代表研究区内落叶松林NPP随时间的变化趋势(式(7),参数含义同上)。当β>0时,表示NPP随时间变化呈上升趋势,反之,NPP则呈下降趋势,其值的大小反映了NPP上升或下降的速率。并采用P值检验其相关性是否显著,P>0.05为不显著,P<0.05为显著,P<0.01为极显著。

(7)

(8)

其累计离差X(i,t)为:

(9)

极差R(t)定义为:

(10)

标准差为S(t)定义为:

(11)

R(t)、S(t)和t满足一般关系式(式(12)):

R(t)/S(t)=c×tH

(12)

式中,c为常数,R(t)/S(t)为重标极差,H为Hurst指数,能够反映NPP时间序列均值的累计离差随时间变化,其值介于0—1之间,可采用最小二乘法拟合求得。若H>0.5时,表明该时间序列具有持续性,即时间序列过去与未来有正的长期依赖性,且H越接近1持续性越强；若H=0.5,表明该序列为随机序列,即过去与未来没有任何相关性；若H<0.5,说明该时间序列具有反向持续性,即过去减小的趋势可能会使未来出现增长趋势,且H越接近0反持续性越强。

2 结果与分析

2.1 模型精度验证

将LPJ模型模拟的吉林省落叶松林NPP与NPP遥感数据和样地实测数据进行比较分析,结果表明,LPJ模型模拟的2000—2019年吉林省落叶松林NPP为450—788 g C m-2a-1,年均值为592 g C m-2a-1,年均增长率为2.81%,略大于NPP遥感数据的年均增长率(1.86%),均值为626 g C m-2a-1,但无显著差异(P>0.05)。由表1可知,2000—2019年LPJ模拟值与NPP遥感数据的MRE(10%)和MAE(50 g C m-2a-1)均比与样地实测值的MRE(41%)和MAE(168 g C m-2a-1)小,表明LPJ模型模拟值在年均值上与NPP遥感数据较接近。

表1 精度评价表

对LPJ模型模拟值与实测值、NPP遥感数据在样地尺度上的年均NPP进行精度验证(见表1),2000—2019年各样地LPJ模型模拟年均NPP值与遥感数据有极显著相关性(R=0.662,P<0.01),2009—2014年LPJ模型模拟年均NPP值与实测NPP亦极显著相关(R=0.748,P<0.01)；2000—2019年LPJ模拟与遥感数据的MRE和MAE较小,两者结果相近,而2009—2014年LPJ模拟值与实测值的MRE和MAE较大,表明两者估算结果差异较大。这是因为实测值中未包含灌草的生产力且LPJ模拟中未考虑人类活动如新植和采伐等,另一方面该实测值为2009—2014年各样地的平均值,并非当年的实测值。

2.2 吉林省落叶松林NPP空间格局分布

2019年吉林省落叶松林NPP空间分布存在较大差异(图3),吉林省北部和南部即吉林市和白山市落叶松林NPP较高,高于800 g C m-2a-1,而东北部和东南部NPP较低,普遍小于600 g C m-2a-1,一方面因为该区域海拔较高,水热条件较差,另一方面可能与该区域落叶松林多为幼龄林有关。LPJ模拟吉林省落叶松林NPP年均值变化范围为-13—47 g C m-2a-1(图3),存在明显的空间异质性。2000—2019年NPP增加的区域占落叶松林面积的63.5%,北部区域近20年NPP增长较快,高于40 g C m-2a-1,而抚松县、安图县等东南部少数样地NPP略有下降。

图3 2019年吉林省落叶松林NPP空间分布及2000—2019落叶松林NPP空间变化Fig.3 Spatial distribution of Larix forests NPP in 2019 and spatial variation from 2000 to 2019NPP: 净初级生产力Net primary productivity

2.3 吉林省落叶松林NPP年际变化

图4 LPJ模型模拟的2000—2019年吉林省落叶松林NPP与NPP遥感数据年际变化Fig.4 Interannual variation of Larix forests NPP simulated by LPJ model and NPP remote sensing data from 2000 to 2019

分析2000—2019年吉林省落叶松林年均NPP(如图4),与NPP遥感数据变化趋势一致,LPJ模型模拟值随时间变化呈波动上升趋势,线性趋势回归方程的R2分别为0.784和0.386。NPP随时间变化整体上先减小后增加,最小值出现在2002年,为450 g C m-2a-1,可能是因为此前人类陡坡开荒以及城镇建设等无序活动导致植被退化严重,此外,落叶松正处于幼龄林阶段也是可能原因之一[33]。

为了分析吉林省落叶松林NPP稳定性,计算2000—2019年各样地落叶松林NPP的变异系数,系统分析各样地落叶松林NPP时间波动特征。总体来看,NPP变异系数为0.07—2.33,均值为0.48,吉林省落叶松林NPP整体上趋于稳定,生态系统受自然灾害和人为干扰小,处于正向演替状态。从局部看,位于敦化市和汪清县范围内的落叶松林NPP变异系数较大(图5),从图2可知,该范围内落叶松林多为幼龄林,平均CV(1.51)明显大于其他各龄组平均CV(约0.1),故其前后时期年均NPP差异较大。

图5 2000—2019年吉林省落叶松林NPP变异系数Fig.5 Coefficient of Variation of Larix forests NPP in Jilin Province from 2000 to 2019

图6 2000—2019年吉林省落叶松林不同龄组NPP年际变化Fig.6 Interannual variation of different age groups of Larix forests NPP in Jilin Province from 2000 to 2019

为了解吉林省落叶松林NPP在2000—2019年的生长趋势及其强度,本研究通过β值来分析NPP变化,并进行F检验。结果表明,吉林省落叶松林NPP近20年总体上处于显著增长趋势(β=14.55,R2=0.784,P<0.05,图4),年均增长量为12.81 g C m-2a-1。这种增长趋势在不同龄组森林中不尽相同(如图6),幼龄林(β=48.43,R2=0.837,P<0.05)和中龄林NPP(β= 7.18,R2=0.379,P<0.05)增长趋势显著,也佐证了幼龄林生长向好而导致其平均CV较大；此外,幼龄林NPP在近20年来小于其他龄组,这也反映出吉林省东北部和东南部NPP较低的格局。而近熟林(β=-1.95,R2=0.046,P>0.05)和成熟林(β=-0.55,R2=0.004,P>0.05)NPP均稍有下降,结合龄组分布,可以解释空间格局中安图县等少数东南地区落叶松林的NPP有下降趋势这一现象。综上,幼、中龄林是吉林省落叶松林NPP增长的主要贡献者。而根据《中国森林资源报告(2014—2018)》,吉林省落叶松林幼、中龄林面积占全省落叶松面积的61.38%,因此,加强对幼、中龄林的经营管理,进一步提高森林质量至关重要。

吉林省落叶松林NPP的Hurst指数介于0.441—0.849之间,均值为0.612,将Hurst指数结果划分为强持续性(0.75幼龄林(0.618)>近熟林(0.462)>成熟林(0.441),可知中幼龄林NPP未来将持续上升,而成熟林和近熟林NPP未来变化趋势呈弱反持续性,由此可以判断未来吉林省落叶松林NPP呈增加趋势。

2.4 吉林省落叶松林NPP对气候变化的响应

分析2000—2019年吉林省落叶松林NPP与温度、降水的年际变化,年均温和年降水均随时间变化稍有上升趋势(图7)。2000—2019降水量波动较大,在2010年出现最大值；而温度是在2010年之前有较大波动,随后逐渐上升。

图7 2000—2019年吉林省落叶松林NPP与降水、温度的年际变化Fig.7 Interannual variation of Larix forests NPP and precipitation and temperature from 2000 to 2019

分析2000—2019年吉林省落叶松林NPP对当年气候因子(温度、降水)以及生长季(5—10月)气候之间的关系表明(如图8),年均NPP与年总降水之间无显著相关性(R=0.063,P>0.05),与年均温之间为显著正相关(R=0.446,P<0.05)。NPP与生长季降水量亦无显著相关性,但与生长季均温呈极显著正相关(R=0.673,P<0.01)。总体上,在2000—2019年期间,温度比降水更能对模型模拟的吉林省落叶松林NPP的年际变化产生影响。

图8 吉林省落叶松林NPP与降水、温度的相关性Fig.8 Correlation between NPP and precipitation and temperature of Larix forests in Jilin Province

3 讨论

3.1 吉林省落叶松林NPP模拟

LPJ模型模拟的吉林省植被分布与实际分布基本一致,落叶松林主要分布于吉林省东部山区和中部低山丘陵区,接近NPP实测值,模型模拟效果良好。对比其他方法对落叶松林NPP的模拟(见表2),与NPP遥感数据接近且无显著差异,但略低于何丽鸿[22]采用BIOME-BGC模型对东北地区长白落叶松林1998—2010年净初级生产力的模拟值,同时其研究也表明,BIOME-BGC模型模拟值较其样地实测NPP偏高,且模型的研究范围、参数设置以及数据源不同也会对模拟结果产生差异。此外,对比赵东升等[3]利用LPJ模型对中国自然植被净初级生产力研究中长白山等森林地区的NPP为600—700 g C m-2a-1,与本研究结果基本一致。

表2 2000—2019年LPJ模型模拟落叶松林NPP与其他研究对比/(g C m-2 a-1)

3.2 吉林省落叶松林NPP对气候变化的响应

本研究发现吉林省落叶松林NPP与年均温、生长季均温呈正相关,与降水则不显著,与何丽鸿[22]利用BIOME-BGC模型对东北地区长白落叶松林NPP对气候变化的响应一致,但李洁等[34]利用CEVSA模型模拟东北地区平均NPP与近50年的气候研究表明其与降水呈显著正相关。对比张凤英等[35]对长江流域森林NPP对气候变化的响应发现,该区域内气温比降水更能对NPP产生影响,表明不同的地形地貌以及植被类型下的NPP对气候的响应不同。本研究以年为步长输出的NPP模拟值在与气候因子进行相关性分析时,将NPP年均值与温度、降水年平均值以及生长季的均温和降水进行分析,在相关性上并不十分显著,研究表明[22],当月NPP与当月的气候存在显著相关性,可以考虑在以后的研究中将模拟的后几十年NPP以月值输出,以便更好的与气候因子进行分析。Ji等[7]对中国森林NPP与气候因子响应研究表明,NPP与降水的相关性较大,其次是温度；而孔蕊等[36]和Liu等[37]分别对长江流域和广西森林NPP的研究中发现温度比降水更能对NPP产生重要影响,因此在研究NPP与气候的相关性时区分研究区是必要的。此外,模型中NPP的输出不光与温度、降水两个重要气候因子相关,还与最冷月气温、最暖月气温、年均相对湿度、日照时数、日照百分率等有关,研究表明[38],利用主成分分析法分析草地生态系统碳储量的气候敏感性时,其敏感性顺序为年降水量>日照百分率>日照时数>3—5月的降水量>最暖月气温>最冷月气温>年均温>6—8月降水量,可见,对于草地生态系统来说,日照百分率比降水和温度对NPP的贡献率更大,但对于不同的森林生态系统是否也是这个贡献率排序,有待研究。

3.3 模拟的不确定性分析

由于数据的限制,本研究的结果可能存在一定的不确定性。首先,用于检验的实测数据在筛选的过程中,可能存在一定的误差；其次,LPJ模型中部分生理生态参数参考了前人研究中的设置,在一定程度上对模拟结果产生影响；第三,模型未考虑人类活动对NPP模拟的影响,从而忽略了天然林保护工程等生态建设对NPP增长的积极作用,而考虑了人为因素的NPP遥感数据在近20年的均值大于LPJ模型模拟均值,也说明了这一点。第四,LPJ模型模拟结果的分辨率为0.1°×0.1°,而NPP遥感数据的分辨率为500 m×500 m,样地的大小为0.06 hm2,不同尺度数据之间的转换也给模拟结果带来一定差异。

LPJ模型作为一个含有众多生理生态参数的生态过程模型,在针对更小尺度植被净初级生产力进行估算时,需考虑研究区的地形地貌、气候属性以及研究物种本身生理生态属性因地因种的改进模型,例如,根据区域内不同优势树种针对性的改进其平均个体属性,如平均单株树木叶面积与边材横截面积的关系,胸径与树高、冠幅的关系,以及四个活组织碳库(叶、细根、边材、心材)的生物量分配,以使模型有更好的局地适用性和更高的精确度,也是今后研究方向之一。

4 结论

本文基于2000—2019年吉林省及其100 km缓冲区内41个气象站点资料采用LPJ模型模拟了2000—2019年吉林省落叶松林近20年的净初级生产力时空变化及其对气候因子的响应,结论如下:

(1)LPJ模型模拟吉林省落叶松林2000—2019年NPP与样地实测值极显著相关(R=0.748,P<0.01),可以用于模拟吉林省落叶松林的NPP。2000—2019年年均NPP为592 g C m-2a-1,年均增长率为2.81%,随时间推移呈现波动增长的趋势(β=14.55,R2=0.784)。最大值出现在2018年,为788 g C m-2a-1,最小值为2002年的450 g C m-2a-1,与遥感数据无显著差异。

(2)除幼龄林外,2000—2019年吉林省落叶松林NPP在时间上趋于稳定,NPP变异系数为0.07—2.33,均值为0.48,整体波动较小。Hurst指数介于0.441—0.849之间,均值为0.612,未来吉林省落叶松林NPP呈增加趋势。

(3)吉林省落叶松林NPP存在明显的空间异质性,吉林省北部和南部区域NPP较高(大于800 g C m-2a-1),也是近20年NPP增长较快的区域,东北部和东南部边界NPP则较低,普遍小于600 g C m-2a-1。

(4)2000—2019年吉林省落叶松林年均NPP与年总降水和生长季降水量之间均不显著(P>0.05),与年均温之间为显著正相关(P<0.05),与生长季均温之间为极显著正相关(P<0.01)。该阶段内年均温比年降水量更能对吉林省落叶松林NPP的年际变化产生影响。