基于CASA模型长江流域植被NPP时空演变及与地形因子的关系

李俊豪,梁娟珠

(福州大学 数字中国研究院(福建),福建 福州 350002)

全球气候变暖、温室气体排放和极地冰川融化等环境问题已变成全人类共同面临的不同寻常的挑战[1]。我国提出力争2030年前实现碳达峰,2060年前努力实现碳中和的目标,由相对减排到绝对减排,进而实现零排放[2]。植被净初级生产力(net primary productivity, NPP)是指在单位面积上植物所积攒的有机物的数量,是由光合作用所固定的有机物总量扣除自身呼吸作用消耗量后的剩余部分,同时,植被NPP是定量描述植被生产能力特征的指标,是衡量区域碳汇能力高低和调节生态过程的重要指标[3],也是陆地生态系统中物质与能量转换研究的基石。总而言之,植被NPP就是判定植被生长状况和固定碳及积累能量的能力。

植被NPP的反演方法主要包括站点实测法和模型法[4-5]。实测法主要在小区域上进行,选取站点的生物量推算NPP,虽然精度高但是费时费力,难以在大范围区域开展实时动态监测。因此,在全球或区域尺度,模型法成为一种重要的主流的研究方法,而站点实测数据更多用于NPP反演结果的精度验证。目前,反演植被NPP发展了一系列的模型反演方法,主要包括气候模型、过程模型和光能利用率模型[6]。

气候模型主要用于植被NPP早期研究,该类模型相对简单且所需要的气象数据易获取,反映的是理想情况下的植被NPP,即一种区域潜在的植被NPP[7]。而过程模型[8]综合考虑了植物光合作用、呼吸作用、有机物的分解以及合成作用等多种因素,但是模型本身比较复杂,且输入的数据较难获得,区域尺度转换难以实现。CASA(Carnegie-Ames-Stanford approach)模型是一个基于光能利用率的具有代表性的模型,在全球以及区域生产力的估算中有着较大的影响力[9]。国内外学者利用CASA模型对草地、植被等不同研究对象在全国、省、市群等不同空间尺度开展了大量研究。LIU等[10]通过CASA模型反演,结果表明,1982—2008年全球草地NPP总量为23.98～25.73 PgC,年均增速为0.03 PgC/a。PIAO等[11]利用CASA模型反演了我国1982—1999年植被NPP,结果表明,植被NPP整体上呈增加趋势,增加速率为0.015 PgC/a。周伟等[12]利用CASA模型模拟我国1982—2010年草地NPP总量为988.3 TgC,其中超过60%的草地表现为增加的趋势。李军玲[13]、戴靓[14]等采用CASA模型并借助其他辅助数据分别对河南省和长江三角洲地区的NPP进行分析。贺倩等[15]基于改进的CASA模型反演得到2010—2015年三江源地区植被NPP均值降低51.99 gC/(m2·a)。YANG等[16]基于CASA模型模拟2001—2016年安徽省植被NPP,均值约为228.27 gC/(m2·a),不同土地利用类型的植被NPP均值从高到低为耕地、草地、林地、农用地、水体、湿地、低植被覆盖区。

在全球植被NPP模拟和碳循环研究中,CASA模型在我国和区域植被生产力的模拟中发挥了很大作用[17];但是由于需要的参数多,精确赋值困难,该模型仍缺乏广泛和细致的应用,尤其是在流域尺度上,比如我国长江流域地区。因此,本研究选用2001—2018年的MODIS数据和中国气象数据驱动CASA模型反演长江流域植被NPP,探究长江流域不同时空尺度的植被NPP分布现状和变化趋势,以期为流域植被恢复及流域生态环境保护提供科学参考。

1 研究区概况、数据来源

1.1 研究区概况

长江流域是世界三大流域之一,总面积约为180万平方公里,约占中国陆地总面积的1/5,流域内植被类型丰富[17],主要为常绿阔叶林、常绿针叶林、落叶针叶林、混交林和灌丛等。流域地形地貌类型复杂多样,地势西高东低,横跨我国地势的三大阶梯,主要流经青藏高原、云贵高原、四川盆地、长江中下游平原。流域内气候类型多样,其中大多数地区处在亚热带季风气候区。流域内气温空间呈东高西低、南高北低的分布趋势,多年平均气温为13 ℃,平均年降水量为1 100 mm,雨量充沛但空间分布不均,呈东南向西北递减的趋势。如图1所示,2001—2018年归一化植被指数(normalized differential vegetation index,NDVI)多年均值范围为0～0.95。

图1 长江流域气象站点分布及二级流域示意图Fig.1 Schematic diagram of meteorological station distribution and secondary river basin in the Yangtze River Basin

1.2 数据来源及预处理

遥感数据:NDVI数据通过NASA官网发布的MOD13A1数据集(https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/products/),空间分辨率500 m,其行列号分别为h25v05、h26v05、h26v06、h27v05、h27v06、h28v05、h28v06,共7景影像,并对其进行格式转换、投影变换、镶嵌、拼接,通过最大值合成法(maximum value composite, MVC)进一步抵消环境因素的干扰,从而获得逐月NDVI数据。2005年、2010年、2015年、2018年植被覆盖数据MCD12Q1,空间分辨率为500 m,提取国际地圈-生物圈计划(International Geosphere-Biosphere Programme,IGBP)全球植被分类图。 MOD17A3HGF数据是NPP的数据产品,来源于NASA官网。

气象数据:降水和气温数据来源于国家气象科学数据中心(http://data.cma.cn/),数据包括长江流域2001—2018年175个气象站点的地面观测数据(平均气温、降水量等),通过克里金插值法获取月降水和月平均气温数据,并使其空间分辨率与NDVI数据一致。太阳辐射使用全球高分辨率地表太阳辐射数据集,其来源于国家青藏高原科学数据中心(http://data.tpdc.ac.cn),时间单位和空间分辨率分别为月和10 000 m。

地形数据:高程、坡度、坡向数据分别来源于地理空间数据云的SRTMDEM 90M、SRTMSLOPE 90M、SRTMASPECT 90M的数据集。

2 研究方法

2.1 NPP的反演

CASA模型充分考虑了环境条件及植物本身特性,将气象及环境调控因子作为输入参数,有效地反演陆地植被NPP,且模型参数可通过遥感或空间分析方法获得。因此,本文采用CASA模型反演长江流域植被NPP,计算公式[18]如下:

Nx,t=Ax,t×εx,t

(1)

式中:Nx,t为像元x处植被在t月的净初级生产力;Ax,t为像元x处植被在t月内所吸收的光合有效辐射;εx,t为像元x处植被在t月内的光能转化率。

Ax,t=Mx,t×Fx,t×0.5

(2)

式中:Mx,t为像元x处植被在t月的太阳总辐射量,MJ/m2;Fx,t为像元x处植被在t月的光合有效辐射吸收比;常数0.5为植被所能利用的太阳有效辐射(波长范围为0.4～0.7 μm)占太阳总辐射的大致比例。

光能转化率的计算公式如下:

εx,t=T1,x,t×T2,x,t×Wx,t×εmax

(3)

式中:T1,x,t和T2,x,t为像元x处植被在t月内的温度因子;Wx,t为像元x处植被在t月内的降水因子;εmax为最大光能转化率。

由于光合有效辐射吸收比(F)与归一化植被指数(INDV)和比值植被指数(IRV)分别存在线性关系,本文根据朱文泉[19-20]的研究,计算出适合长江流域的每个植被类型所对应的INDV和IRV的最大值,赋予最小值,使用2个参数分别进行计算后,取其平均值作为F的值。结合文献[21],设置最大光能转化率(εmax),并从数据获取层提升模型运用的区域适用度,从而实现对静态参数的本地化(表1)。

表1 长江流域植被类型及其参数设置Tab.1 Vegetation types and parameter setting in the Yangtze River Basin

2.2 植被NPP趋势分析

Theil-Sen Median趋势分析法可用于逐像元分析植被NPP在时间序列上的变化趋势及空间差异性,该方法无需数据服从特定分布且不受异常值的干扰,对离散数据具有较强的规避能力[22]。其计算公式为

(4)

式中:C为植被NPP变化速率,C>0说明植被NPP呈上升趋势,C<0 则植被NPP为下降趋势;Ni为i年植被NPP值;Nj为j年植被NPP值。

由于Theil-Sen Median趋势分析法无法进行显著性检验,与Mann-Kendall非参数检验方法结合,能较好地克服噪声对分析结果的影响[23]。因此,植被NPP变化趋势的显著性通过Mann-Kendall统计检验法判断,相关计算公式如下:

检验值:

(5)

符号函数:

(6)

(7)

标准化统计量:

(8)

式中:sgn是符号函数。趋势检验的方法是零假设,H0:R=0;当|Z|>Z1-α/2时,拒绝零假设。其中,Z1-α/2为标准正态方差;α为显著性检验水平。当|Z|分别大于等于1.28、1.96时,表示分别通过了置信度为90%、95%的显著性检验。

2.3 地形因子的分析

地形因子不仅是可以描述地貌形态的参数,也是植被NPP垂直分布及多样性分布的关键因素。根据长江流域复杂的地形地貌特征,选取的地形因子包括海拔、坡度和坡向。借鉴《水土保持综合治理规划通则》[24],结合长江流域地形情况,使用Arcgis 10.2软件对海拔、坡度和坡向重分类并进行等级划分。地形因子分类定级见表2。将研究区各级地形因子与NPP进行叠加分析,统计各级地形因子对应的NPP均值,分析NPP随地形因子等级变化的特征。本文定量分析NPP与高程、坡度和坡向之间的关系,在各地形因子等级选取若干采样点,通过斯皮尔曼等级相关系数研究在各地形因子等级中变量之间的相关关系。

表2 地形因子分类定级Tab.2 Terrain factor classification and grading

3 结果与分析

3.1 植被NPP反演结果验证

3.1.1反演结果与遥感数据产品对比

由于长江流域植被NPP实测数据较难获取,基于2001—2018年长江流域的MODIS NPP数据产品(MOD17A3HGF),经过处理得到长江流域2001—2018年多年NPP的平均值,在研究区内随机选取675个样点,将CASA模型反演的数据与MOD17A3HGF产品进行拟合(图2),以验证本文研究结果。结果表明,两者之间的拟合程度好,本文反演值要高于MOD17A3HGF产品。其原因是:MOD17A3HGF产品低估了灌丛、草原草甸的生产力,与其他植被类型的拟合程度较好[25]。MOD17A3HGF产品在长江流域的城镇及周边地区,如长三角地区、鄱阳湖水系北部、洞庭湖水系北部、宜昌至湖口中部、汉江东部等出现缺值,突出本文CASA模型用于长江流域研究的优势,并且二者相关性R2为0.717(P<0.01),具有强相关性,说明该模型可以有效地反演长江流域植被NPP,可为之后的研究提供基础数据。

图2 反演值与MOD17A3HGF产品的比较Fig.2 Comparison of the inversion values with the MOD17A3HGF product

3.1.2对比验证

为了进一步验证长江流域植被NPP数据的可靠性,将本文的反演值与其他模型的数据进行比较。由表3可知:本文不同植被NPP的反演值在实测值范围内,并且与前人的反演结果一致,均表现为乔木类NPP较高、稀疏植被NPP较低,说明参数本地化后的CASA模型反演的NPP数据基本可靠,可以很好地反映长江流域植被NPP的实际状态。

表3 植被NPP反演值与实测值的比较 单位:gC/(m2·a)

3.2 长江流域植被NPP时间变化特征

3.2.1植被NPP年际变化特征

由CASA模型反演的NPP数据绘制2001—2018年长江流域植被NPP年际变化图,如图3所示。从图3可以看出:长江流域在这18年间表现出增长的趋势,但是增幅较小,NPP年均值范围为511.92～638.25 gC/(m2·a),NPP多年均值为572.72 gC/(m2·a);峰值638.25 gC/(m2·a)出现在2007年,比NPP多年均值高了11.44%;谷值511.92 gC/(m2·a)出现在2001年,比NPP多年均值低了10.62%。由于受到气候干旱[29]的影响, 2001年NPP均值明显小于其他年份。总体而言,2001—2007年NPP年均值变化明显,呈波动式上升,2008—2018年变化相对趋于平稳,表明长江流域生态环境逐渐恢复,区域固碳能力有所增强。

图3 2001—2018年长江流域植被NPP年际变化图Fig.3 Interannual change of NPP in the Yangtze River Basin from 2001 to 2018

3.2.2植被NPP年内变化特征

对2001—2018年长江流域植被NPP取月平均值合成,并计算月NPP总量占年NPP总量的比例(贡献度),如图4所示。由图4可知:长江流域NPP年内变化层次感比较明显,表现出倒“U”型变化特征和季节性。夏半年(4—9月)NPP月均值处在较高位,均大于40 gC/(m2·a),且夏半年NPP总量占全年的77.61%;冬半年(10月—次年3月)NPP月均值整体较低,其中11月—次年2月NPP月均值小于20 gC/(m2·a)。夏季6—8月NPP月均值都超过80 gC/(m2·a),NPP总量占全年NPP的44.66%;春季3—5月、秋季9—11月NPP月均值较低,分别占NPP总量的27.22%、20.35%。冬季对全年NPP的贡献度仅为7.77%,不足夏季的50%;7、8月NPP月均值最高并且贡献度最大,均值分别为90.85 gC/(m2·a)和83.92 gC/(m2·a),贡献度分别是15.86 %和14.65%;12、1月NPP月均值最低,均值分别为13.81、13.43 gC/(m2·月),贡献度分别为2.35 %和2.41 %。植被NPP表现出的年内变化特征,原因在于:长江流域的植被覆盖度高的地区位于湿润地区,属于亚热带季风气候,夏季高温多雨,冬季寒冷干燥,从而引起植被在夏半年生长迅速,在冬半年生长缓慢。

图4 长江流域植被NPP月均值变化图Fig.4 Monthly mean change chart of NPP of vegetation in the Yangtze River Basin

3.3 长江流域植被NPP空间分布及变化特征

3.3.1植被NPP空间分布格局

对2001—2018年长江流域植被NPP取平均值,绘制空间分布图。如图5所示:长江流域多年植被NPP均值在0～1 700.13 gC/(m2·a)范围内,流域内植被NPP空间分异明显,整体分布格局上呈阶梯状分布,由第一阶梯向第三阶梯递增,自西北向东南递增,植被NPP高值聚集在青藏高原山麓地区和鄱阳湖水系地区,而青藏高原部分地区和长三角地区植被NPP较低。

图5 长江流域植被NPP的空间分布图Fig.5 Spatial distribution diagram of the vegetation NPP in the Yangtze River Basin

为了更好地了解长江流域空间分布格局,分区统计各二级流域2001—2018年植被NPP。由图6可以看出:12个二级流域中,乌江、嘉陵江、汉江、洞庭湖水系、湖口干流、鄱阳湖水系共6个水系植被NPP均值超过600 gC/(m2·a),其中鄱阳湖水系植被NPP均值679.99 gC/(m2·a)最大,得益于鄱阳湖所在的江西省高度的植被覆盖度;而金沙江石鼓以上水系的NPP年均值仅为274.00 gC/(m2·a),远远低于其他流域,其他5个二级流域(金沙江石鼓以下、岷沱江、宜宾至宜昌、宜昌至湖口、太湖水系)植被NPP均值为500～600 gC/(m2·a)。

图6 长江流域二级流域的植被NPP均值统计图Fig.6 Statistical graph of NPP mean value in the secondary level basin of the Yangtze River Basin

3.3.2植被NPP空间变化特征

本文采用Theil-Sen Median趋势分析法和Mann-Kendall统计检验法将长江流域植被变化趋势划分为显著减少(C<0和|Z|≥1.96)、轻微减少(C<0和1.28≤|Z|<1.96)、基本不变(|Z|<1.28)、轻微增加(C>0和1.28≤|Z|<1.96)和显著增加(C>0和|Z|≥1.96)总共5个分类,绘制2001—2018年的变化速率以及变化趋势图。从图7可知:长江流域植被NPP平均变化速率为0.066 gC/(m2·a),有17.76%的面积表现为增加趋势,主要分布在岷沱江南部、嘉陵江西部、乌江西部和洞庭湖东部地区,其中岷沱江南部、嘉陵江西部、乌江西部增幅超过5 gC/(m2·a),说明这些区域植被生长环境较好,为植被生长提供丰富的养料。长江流域植被NPP减少趋势的像元占全域的16.16%,显著减少的像元占比仅5.48%,主要分布于汉江西部与长江流域西南部,这与当地近年来气温下降[30]有关,气温较低会导致植被生长放缓。长江流域植被NPP基本不变的像元占比66.08%,主要分布在长江中部和西北部。由此看出,长江流域整体植被生长状况和固碳能力在改善,但是局部地区有恶化的趋势。

图7 2001—2018年长江流域植被NPP变化趋势图Fig.7 Trend chart of NPP change of vegetation in the Yangtze River Basin from 2001 to 2018

3.4 长江流域植被NPP变化与地形因子的相关性分析

3.4.1植被NPP与高程的关系

为了解植被NPP与高程之间的关系,统计分析2001—2018年各高程带植被NPP年均值变化特征,并在各高程带中选取若干个采样点分析植被NPP与高程的相关性,如图8所示。由图8(a)可以看出:2007年,>500～1 000 m和>1 000～2 000 m高程带的植被NPP均超过>2 000～3 000 m高程带的植被NPP,原因是该年降水量激增,在政府有效的调度下[31],干旱缺水的地区有了水分补给,促进了这些地区植被的生长。其他时间,长江流域各高程带植被NPP变化趋势与长江流域植被NPP年际变化趋势整体上一致,都呈现波动性。高程在3 000 m以下时,植被NPP年均值随着各级高程带的增加而增加;高程在3 000m以上时,随着高程的增加植被NPP年均值减少。因此,2007年植被生长最佳的是>1 000～2 000 m高程带,其他时间是>2 000～3 000 m高程带植被NPP年均值最高。由图8(b)可以看出:<200、200～500、>500～1 000、>1 000～2 000 m的高程带中,植被NPP与高程呈正相关,其中>1 000～2 000 m高程带呈显著正相关;>2 000～3 000 m高程带呈显著不相关;高程>3 000～4 000 m和>4 000 m的区域呈显著负相关;整个研究区植被NPP与高程的相关系数是-0.487(P<0.01)。说明高程在>1 000～2 000 m地区,植被NPP随高程的增加显著增加;高程在>2 000～3 000 m地区植被NPP变化比较稳定,而高程在>3 000 m地区植被NPP随高程的增加急剧减少。综上,高程对NPP的解释较好,其中>2 000～3 000 m高程带上植被生长状况最佳。

3.4.2植被NPP与坡度的关系

2001—2018年长江流域地区植被NPP与坡度的关系如图9所示。从图9(a)可以看出:2001—2018年长江流域植被NPP随坡度的增加而增加,且各坡度范围植被NPP变化趋势与研究区NPP年际变化趋势保持一致。图9(b)中各坡度分级范围与植被NPP均呈正相关性。<6°坡度的区域易受到人类活动的影响,因此该地区的坡度对NPP的解释较差。在6°～15°和>15°～25°坡度范围,NPP与坡度呈显著正相关性,即植被NPP随坡度的增加而增加,这些地区是开展退耕还林、小流域治理、生态移民搬迁等一系列生态工程的重点区域,并且人类活动较少,因而该地区植被生长旺盛,NPP增长。坡度>25°地区多为陡峭的山崖,水分难以保持,气候条件很难控制,植被面临光照不足、温度变化等挑战;相较坡度<6°地区,该地区的坡度对NPP的解释略好。因此,<6°坡度地区的NPP是各坡度区间最低的,并且该地区坡度对NPP的解释较差,在6°～15°和>15°～25°坡度对植被NPP的解释较好。

图9 2001—2018年长江流域植被NPP与坡度的关系Fig.9 Relationship between NPP and slope in the Yangtze River Basin from 2001 to 2018

3.4.3植被NPP与坡向的关系

2001—2018年长江流域地区植被NPP与坡向的关系如图10所示。由图10(a)可知:平坡、阴坡、半阴坡、阳坡、半阳坡植被NPP变化趋势与NPP年际变化趋势一致,但是各坡向NPP年均值相差不大。从图10(b)可以看出:只有半阴坡与NPP多年均值呈现显著正相关,其他坡向与NPP相关性系数较低,均小于0.2。因此,除半阴坡外,其他坡向对NPP影响不大。

图10 2001—2018年长江流域植被NPP与坡向的关系Fig.10 Relationship between NPP and slope direction in the Yangtze River Basin from 2001 to 2018

4 结论与讨论

本文基于CASA模型,针对NDVI和最大光能利用率εmax进行本地化,反演出以月为时间尺度的长江流域植被NPP,并分析2001—2018年研究区植被NPP时空变化规律及与地形因子的关系,得出以下结论:

1)通过与MODIS NPP数据产品以及前人的相关研究进行对比,CASA模型可以有效地反演长江流域植被NPP,反映长江流域植被生长的实际状态。长江流域空间格局显示东高西低,2001—2007年NPP均值呈波动式上升,2008—2018年变化相对趋于平稳;长江流域植被NPP具有年内季节性的表现,夏季NPP对全年NPP总量的贡献接近一半。

2)长江流域在2001—2018年18年中植被NPP增加的面积大于减少的面积,增加的区域主要分布在岷沱江水系南部、嘉陵江水系西部、乌江水系西部和洞庭湖水系东部地区,因此长江流域植被生长状况整体来看有所改善,然而局部地区植被的生长存在恶化的现象。

3)关于地形因素的分析发现,相较于坡度、坡向,高程对长江流域植被NPP的解释较强。长江流域植被NPP与不同高程带的响应不同,在<3 000 m的地区植被NPP随着高程的增加而增加,在高程>3 000 m的地区植被NPP显著减少,其中最适合植被生长的高程带是>2 000～3 000 m。在6°～15°和>15°～25°坡度以及半阴坡对植被NPP的解释较好。

本文反演的NPP数据集时间序列长,能够较好地反应植被在不同时空尺度下的变化。地形因子反应是植被在自然环境条件下的综合表现,高程、坡度和坡向等地形因子的改变,影响降水量、气温高低、土壤肥力等资源在时空中的分配,进而影响植被生产能力。分析地形因子与植被NPP之间的关系能够有效地揭露出不同生长条件下植被的分布。除了地形因子的影响,气温、太阳辐射、降水等因子均会直接影响植被NPP的时空变化,各因子之间相互限制和相互关联。如何深入分析植被NPP的各因素之间的综合影响,是研究植物生态环境和评估区域碳收支的基石,今后将进一步研究各自然因子对植被NPP的综合作用。