近20年来西藏地区NPP时空变化及影响因素

凡晚晴吴 华樊风雷

(1.西藏大学 工学院,西藏 拉萨850012;2.西藏大学 高原地表环境遥感监测联合实验室,西藏 拉萨850000;3.华南师范大学 地理科学学院,广东 广州510631)

植被净初级生产力(net primary productivity,NPP)是指植被在单位地表面积单位时间内所积累的有机物质总量,由植物通过光合作用途径所固定的碳总量(gross primary productivity)减去植物自养呼吸(autotrophic respiration)的碳消耗所得到[1-2]。陆表生态系统是碳循环过程中受人类活动影响最大的部分,陆表生态系统中的碳储量对调节全球碳平衡及温室气体的减少至关重要,如何及时有效地获取碳储量成为当前研究领域的热点。NPP不仅反映了自然环境条件下植物的生产能力,也是研究地表碳循环过程,估算地球物质支持能力,反映陆表生态系统环境质量的重要生态指标之一,它的大小能够直观地反映陆表生态系统碳储量的多少[3-4]。西藏地区地形复杂,气候多变,是全国乃至全球气候变化响应的敏感区,其生态系统的动态变化将直接影响我国碳平衡、气候及生态环境,准确地掌握区域NPP动态变化规律及其影响因素,对生态系统的研究和生态保护建设的实施具有重要意义[5]。

近年来,随着遥感技术的发展,已有学者[6]利用MOD17A3数据针对NPP开展了大量相关研究。李登科等[7]对2000—2015年中国陆地植被NPP变化特征进行了分析,认为我国植被NPP在波动增加,并呈现东南高西北低的趋势;国志兴等[8]分析了我国东北地区植被净初级生产力的时空分布格局,并认为降水及土地类型变化是造成东北地区植被NPP变化的主要因素;王强等[9]分析了2004—2014年横断山区植被NPP的空间分布,并讨论了不同植被生态系统NPP的年际变化;李传华等[10]建立NPP人为影响模型分析石羊河流域植被NPP时空分布,认为降雨与人为影响是植被NPP变化的主导因子;王芳等[11]基于MODIS NPP数据对安徽省NPP变化规律及驱动分布进行了分析,并认为江淮地区植被NPP变化的主要影响因素是降水和人类活动,城市化对植被NPP变化影响日益增加。综上所述,研究区集中在全国尺度[12-16]、区域尺度[17-19],包括东北地区[20]、西部地区[21]、西南地区[22-23]、青藏高原[24-28]、东南地区等[29],主要针对NPP时空变化及气候因子之间响应关系进行研究,基于MODIS NPP数据对西藏地区进行NPP时空变化及自然、人为影响因子的研究较少,且尚未有人利用地理探测器针对西藏地区NPP进行影响因子贡献率计算,探讨区域NPP变化的主要驱动因子及影响因子交互作用探测。本文采用MOD17A3 NPP数据分析西藏地区NPP的时空分布格局,同时对NPP与气候因子(气温、降雨量)、土地利用覆被、地形因子(海拔、坡度)、植被类型进行影响因子分析,从而在一定程度上揭示其时空动态变化的原因,为该地区植被建设管理、生态环境保护、自然资源开发和绿色生态建设提供科学依据。

1 研究区与数据来源

1.1 研究区概况

本文的研究范围覆盖整个西藏自治区。西藏自治区位于78°24′—99°06′E和26°52′—36°32′N之间,全自治区总面积122.84万km2。研究区内地形复杂,地貌丰富,境内有山地、高原、盆地和峡谷等多种地形,整体地势呈西北高东南低,海拔落差大的特点,在63～8 848 m之间;西藏地区南北气候差异较大,藏北主要为大陆性气候,藏南主要为热带山地季风湿润气候,年均温为-31.1～24.5℃,年均降雨量为9～3 105 mm;土壤类型以草原土、草甸土、寒漠土为主;植被类型主要有高寒荒漠、高寒草原、高寒草甸、灌丛、高山植被,藏东南地区有大量的针叶林和阔叶林。

1.2 数据来源与处理

本文所使用的数据主要有:MODIS NPP数据、气象数据和植被类型空间分布数据。

(1)MODIS NPP数据为美国地质勘探局官网(http:∥www.ntsg.umt.edu/)发布的MOD17A3逐年合成产品数据,空间分辨率为500 m,时间分辨率为1 a。由于研究区面积较大,每年选取4景MODIS影像覆盖整个研究区,时间序列为2000—2020年,一共选取84景影像;利用MRT(MODIS reprojection tool)软件对影像数据进行镶嵌和等重投影预处理,剔除异常值,利用研究区矢量边界裁剪影像,并对其进行重采样,将分辨率设置为1 km,获得该区域2000—2020年NPP数据。

(2)气象数据来源于国家地球系统科学数据共享服务平台(http:∥www.geodata.cn),时间为2000—2020年,数据包括中国1 km分辨率栅格逐月降水量数据集和气温数据集,数据集均经过验证,验证结果可靠,对月数据进行合成得到逐年气象数据,取气温均值和年降雨量数据进行分析[11]。

(3)DEM数据为地理空间数据云(https:∥www.gscloud.cn/)的ASTER GDEM V2数字高程数据,空间分辨率为30 m,同样对高程数据进行重采样等预处理,最后统一为1 km空间分辨率。利用ArcGIS软件提取坡度,并对坡度、海拔进行分级。

(4)2000,2005,2010,2015和2020年土地利用覆被数据,空间分辨率为1 km;中国1∶100万植被类型空间分布数据来自中科院资源环境科学与数据中心(http:∥www.resdc.cn/);依据土地利用类型将研究区划分为耕地、林地、草地、水域、城镇和未利用土地6个一级地类;为了定量分析不同植被类型NPP之间的差异,依据研究区主要植被类型将其划分为荒漠、草原、草甸、灌丛、阔叶林、针叶林、高山植被和栽培植被8种植被类型;根据不同生态区的温度和干湿状况等自然条件,将研究区划分为7个自然生态区。

1.3 研究方法

1.3.1 偏离分析 对研究区NPP数据取年平均值,利用每年的NPP与年均NPP之差反映该年NPP的增减程度,若偏离值大于0,则表示该年NPP大于近20 a均NPP水平;若偏离值小于0,则表示该年NPP小于近20 a均NPP水平[11]。计算公式为:

式中:D为NPP偏离值;NPPi为第i年的NPP值;为多年均NPP值;i为年份。

1.3.2 趋势分析 基于一元线性回归趋势分析法拟合出线性倾向值,对2000—2020年西藏地区NPP时间变化趋势进行估算[9]。计算公式为:

式中:θSlope为线性倾向值;NPPi为第i年的NPP;i为年份;n=21,当θSlope＞0时,表示随时间推移,NPP的变化趋势是上升的,θSlope＜0时,表示随时间推移,NPP的变化趋势是下降的。

1.3.3 变异系数 采用变异系数(coefficient of variation,CV)反映NPP空间变异程度[30]。计算公式为:

式中:Cv为NPP变异系数;NPPi为第i年的NPP;为多年均年NPP值。

1.3.4 相关性分析 基于像元对NPP与降雨、气温的相关系数分别进行计算,分析气候因子与NPP之间响应关系的密切程度,同时利用相关系数T检验法对结果进行显著性检验,确定NPP与气候因子显著相关区域[9]。计算公式为:

式中:Rxy为x,y两变量的相关系数;xi为第i年x变量的值;yi为第i年y变量的值;ˉx为变量x均值;ˉy为变量y均值,n=21。其偏相关系数计算公式为:

式中:Rxy,z为固定自变量z后因变量x与自变量y的偏相关系数。

1.3.5 地理探测器 基于地理探测器分析影响因子与NPP变化之间的非线性关联,地理探测器是分析空间分层异质性的工具,包括因子探测、风险区探测、生态探测和交互作用探测共4个探测器[31]。利用ArcGIS软件对研究区以1 m×1 m为一个单元布设随机采样点,共随机采样1 000个点,并提取采样点的NPP值及影响因子相关值,利用GeoDetector软件对因变量Y和自变量X进行交互分析[32]。计算公式为:

式中:h=1,…,L为X或Y的分区;Nh为层h的单元数,N为全区的单元数为Y在层h和全区的方差;SSW为层内方差之和;SST为全区总方差。q值越大表示因变量的空间分异性越大,自变量对因变量的解释力越强,反之越弱[31]。

2 结果及分析

2.1 NPP的空间分布特征

采用2000—2020年的NPP数据对研究区近20 a的平均NPP进行计算(图1)。西藏地区NPP分布具有显著的空间异质性。年平均值NPP为0～1 930.48 g/(m2·a),总体均值为137.36 g/(m2·a),多数区域年平均植被NPP处于0～900 g/(m2·a)之间,占研究区总面积的96.79%。

图1 西藏地区2000—2020年平均NPP空间分布

西藏生态地理分区自东南到西北分别为中亚热带湿润的东喜马拉雅南翼地带(VA6)和高原温带湿润的川西藏东高山深谷地带(HIIAB1),高原温带半干旱的藏南山地地带(HIIC2),高寒半湿润的果洛那曲丘状高原地带(HIIB1),高寒半干旱的羌塘高原湖盆地带(HIC2),高原温带干旱的阿里山地地带(HIID3)和高寒干旱的昆仑高山高原地带(HID1)。植被类型自东南到西北的分布分别为常绿阔叶林、针叶林、山地灌丛及草甸、高寒灌丛、高寒草甸、高寒草原和高寒荒漠。故不同生态地理分区的植被NPP空间分布差异明显,且自东南向西北呈逐渐降低的趋势。

将NPP划分为6个等级:①0～100 g/(m2·a)主要分布在昆仑高山高原地带、阿里山地地带、羌塘高原湖盆地带的和藏南山地地带的西北部,植被类型主要以高寒荒漠、山地荒漠和高寒草原为主,约占研究区的69.11%;②100～200 g/(m2·a)主要分布在羌塘高原湖盆地带的东南部及藏南山地地带的中部,植被类型主要以高寒草原和山地草原为主,约占研究区的15.01%;③200～300 g/(m2·a)主要分布在果洛那曲丘状高原地带及藏南山地地带的东北部,川西藏东高山深谷地带也有零星分布,植被类型多为高寒草甸,约占研究区的7.11%;④300～600 g/(m2·a)主要分布在川西藏东高山深谷地带的东部和西南部,植被类型以高山灌丛为主,约占研究区的3.95%;⑤600～900 g/(m2·a)主要分布在东喜马拉雅南翼地带,藏南山地地带和川西藏东高山深谷地带的南部也有零星分布,植被类型多为针叶林和灌丛,约占研究区的1.62%;⑥＞900 g/(m2·a)主要分布在东喜马拉雅南翼地带,植被类型主要为常绿阔叶林,约占研究区的3.21%。对研究区按生态地理分区进行划分,2000—2020年东喜马拉雅南翼地带的植被NPP平均值最大,为911.97 g/(m2·a);昆仑高山高原地点的最小,为29.20 g/(m2·a),按均值大小顺序排列为:东喜马拉雅南翼地带〔911.97 g/(m2·a)〕＞川西藏东高山深谷地带〔228.09 g/(m2·a)〕＞果洛那曲丘状高原地带〔140.26 g/(m2·a)〕＞藏南山地地带〔109.02 g/(m2·a)﹞＞羌塘高原湖盆地带〔53.09 g/(m2·a)〕＞阿里山地地带〔33.20 g/(m2·a)〕＞昆仑高山高原地带〔29.20 g/(m2·a)〕。

2.2 NPP的时间变化特征

2000—2020年西藏地区年均NPP整体上呈波动上升趋势(图2)。近20 a间年均NPP主要在128～149 g/(m2·a)之间波动,总体平均值为137.36 g/(m2·a)。2006年NPP为近20 a间的最大值,达到148.06 g/(m2·a),比年平均值高7.79%;2014年NPP为最小值,仅有128.04 g/(m2·a),比年平均值低6.78%。由图2可知,近20 a来西藏地区NPP变化大体可分为3个阶段,2000—2006年上升,2006—2014年围绕年平均值波动下降,2014—2020年再次上升,且2020年NPP高于2000年。2006—2014年,由于年均温的影响NPP产生波动,其中2006年和2009年NPP增加可能与年均温较高有关。

图2 西藏地区2000—2020年NPP年际变化

通过计算NPP的变化趋势斜率和变异系数来反映西藏地区2000—2020年NPP的年际变化和空间分布特征。由图3可以看出,研究区NPP的平均变化趋势斜率为0.04,表明近20 a西藏地区NPP整体呈上升趋势,且平均以0.04 g/(m2·a)的年变化速率增长。NPP的线性倾向斜率范围为-65.28～42.36,其中83.09%的研究区介于-1～1之间。东喜马拉雅南翼地带NPP明显下降,下降趋势大于-5;川西藏东高山深谷地带的南部NPP呈下降趋势,在-5～-2之间;而在果洛那曲丘状高原地带、川西藏东高山深谷地带的北部及藏南山地地带的东北部NPP呈上升趋势,在1～2之间。研究区NPP的变异系数范围为0.025～2.395,大部分分布在0.1～0.2之间,约占研究区的67.32%。东喜马拉雅南翼地带的东部及川西藏东高山深谷地带的西部NPP波动较大,变异系数多数为0.25;而藏南山地地带、果洛那曲丘状高原地带、羌塘高原湖盆地带、阿里山地地带及昆仑高山高原地带的大部分地区的NPP波动较小,在0～0.2之间。

图3 西藏地区2000—2020年NPP变化趋势斜率及变异系数

2.3 不同植被类型的NPP变化

本文将植被类型划分为8种对NPP进行分析(图4),不同植被类型20 a的平均NPP差别较大,其中阔叶林的NPP最高,为1 107.09 g/(m2·a);其次是栽培植被,为762.04 g/(m2·a);荒漠的NPP最低,为29.76 g/(m2·a),按NPP均值大小顺序排列为:阔叶林〔1 107.10 g/(m2·a)〕＞栽培植被〔762.04 g/(m2·a)〕＞针叶林〔452.24 g/(m2·a)〕＞灌丛〔203.71 g/(m2·a)〕＞草甸〔113.91 g/(m2·a)〕＞高山植被〔60.92 g/(m2·a)〕＞草原〔49.26 g/(m2·a)〕＞荒漠〔29.76 g/(m2·a)〕。对研究区内不同植被类型的变化趋势斜率及变异系数进行统计(表1),大部分植被类型NPP的趋势斜率处于增加趋势,其中阔叶林、针叶林和栽培植被的趋势斜率均小于0,分别为-6.606,-1.589,-3.790;所有植被类型的NPP变异系数均小于0.2,表示NPP空间变化较为稳定。

表1 不同植被类型NPP变化趋势斜率及变异系数

图4 不同植被类型的年均NPP

2.4 不同地形条件的NPP变化

地形因素是造成高原地区植被净初级生产力空间异质性的主要影响因素之一,它通过改变区域的水热条件、土壤湿度和太阳辐射等,进而影响NPP空间格局。研究区整体上表现出东南低西北高、地形起伏大的特征,海拔最高处为8 848 m,坡度最大值为88.47°。本文基于数字高程模型提取西藏地区海拔和坡度两个地形因子,分析不同地形条件下西藏地区NPP的空间分布特征。将海拔分为10个等级(表2),结果显示,研究区内有90%的区域海拔在4 000～6 000 m之间。西藏地区的NPP垂直分布差异明显,随海拔的增加NPP表现出先增加后减少的趋势;500～1 000 m海拔范围的NPP均值最高,为1 397.41 g/(m2·a),主要位于东喜马拉雅南翼地带,植被类型多为常绿阔叶林。海拔4 000 m以上的高原地区NPP均值较低,其中海拔6 000 m以上NPP均值最低,仅为39.62 g/(m2·a)。利用ArcGIS软件对数字高程模型进行坡度分析,将研究区的坡度分为9个等级统计不同坡度等级的NPP变化。根据统计结果显示(表3),研究区内大部分区域坡度小于20°,占研究区总面积的67.13%。西藏地区NPP随坡度增加呈增长趋势,且坡度大于50°之后NPP增长变缓。

表2 不同海拔分级的年均NPP

表3 不同坡度分级的年均NPP

2.5 NPP与气象因子相关性分析

2.5.1 气候因子的变化趋势 本文利用降雨、气温数据的年平均值分析气候因子与NPP变化之间的响应关系,分别对近20 a来的气象数据进行一元线性回归分析,以斜率作为变化趋势。由图5可知,年降雨量趋势斜率大于0,为0.355 mm/a,年降雨量均值范围为320.21～419.82 mm,多年均值为379.27 mm,其中2020年均降雨量为最大值,2009年最小;气温趋势斜率大于0,为0.044℃/a,年气温均值范围为-2.91～-1.74℃,多年均值为-2.34℃,年均气温在2009年为最大值,2000年最小。总体上表示2000—2020年西藏地区降雨和气温呈波动上升趋势,研究区气候呈现暖湿化趋势。

图5 西藏地区2000—2020年降雨及气温变化趋势

2.5.2 NPP与气候因子相关性结果 通过计算2000—2020年NPP与降雨、气温之间的偏相关系数,分析气候因子对NPP的影响。由图6可知,NPP与气候因子之间的响应关系在空间分布上差异明显。其中NPP与年降雨偏相关系数分布范围为-0.897～0.836,有45.16%的研究区NPP与降雨呈正相关,主要分布在昆仑高山高原地带的西部、阿里山地地带的南部、羌塘高原湖盆地带的中部及东喜马拉雅南翼地带的大部分地区;54.84%的研究区呈负相关,主要分布在昆仑高山高原地带的东部、果洛那曲丘状高原地带及川西藏东高山深谷地带。对偏相关系数进行显著性检验,有2.33%的区域通过p＜0.05的检验,主要分布在昆仑高山高原地带的中部、阿里山地地带的南部及羌塘高原湖盆地带的中部地区。NPP与年气温的偏相关系数分布范围为-0.852～0.904,研究区有75.80%的区域NPP与气温呈正相关,主要分布在昆仑高山高原地带、阿里山地地带、果洛那曲丘状高原地带、羌塘高原湖盆地带和川西藏东高山深谷地带的大部分地区;24.20%的研究区呈负相关,主要分布在羌塘高原湖盆地带的南部、藏南山地地带的西部及东喜马拉雅南翼地带的部分地区。研究区有5.53%的区域达到了极显著水平(p＜0.01),主要分布在阿里山地地带及果洛那曲丘状高原地带,东喜马拉雅南翼地带和川西藏东高山深谷地带也有零星分布。总体而言,NPP与降雨多呈负相关,与气温多呈正相关。

逐像元对NPP与气候因子的复相关关系进行分析(图6)。NPP与气候因子的复相关系数介于0～0.902之间;其中,昆仑高山高原地带的中部、果洛那曲丘状高原地带及川西藏东高山深谷地带的西北部NPP与气候因子复相关性较强,羌塘高原湖盆地带、阿里山地地带的北部及藏南山地地带的大部分地区复相关性较弱;西藏地区NPP与气候因子之间的复相关性空间分布差异较大,可能是由于植被类型及海拔高度的影响。

图6 NPP与降雨—气温的偏相关系数及复相关系数

2.6 影响因子探测

基于随机采样点的像元值,对2000,2005,2010,2015和2020年的NPP与年降雨量、年气温、海拔、坡度、土地利用类型、植被类型共6个影响因子进行空间探测,分析影响研究区NPP时空变化的主要驱动因子。从表4可以看出,不同影响因子对NPP空间分异性的解释力差异明显。其中降雨量对NPP的解释力最强,其次是海拔,土地利用类型的变化对NPP的空间分异影响最弱,各影响因子解释力大小依次为:降雨量(0.896 8)＞海拔(0.885 4)＞气温(0.808 6)＞坡度(0.743 5)＞植被类型(0.673 4)＞土地利用类型(0.439 5)。从时间上分析(图7),2000—2020年不同影响因子对植被NPP空间分异性的解释力具有不同的趋势变化。植被类型解释力总体呈增加趋势,年气温、坡度和土地利用类型解释力总体呈下降趋势,年降雨量和海拔解释力较稳定,波动幅度不大,且年降雨量为第一影响因子,海拔为第二影响因子。2000—2005年解释力总体呈增加趋势,年气温和坡度解释力呈下降趋势;2005—2010年解释力总体呈下降趋势,坡度和海拔解释力呈上升趋势,第一影响因子变为海拔;2010—2015年除年气温解释力呈下降趋势外,其余各影响因子解释力均呈上升趋势,且坡度变为第三影响因子;2015—2020年坡度和土地利用类型解释力呈显著下降趋势,其余各影响因子解释力呈上升趋势,年降雨量、海拔和年气温分别为第一、第二、第三影响因子。总体而言,研究区自然因子对NPP空间分异性的解释力强于人为因子。

图7 2000—2020年西藏地区影响力q值变化趋势

通过风险探测,结果均通过显著性检验(p＜0.05),对NPP空间分布具有明显影响的因子主要有气温、植被类型和坡度。通过生态探测分析影响因素对NPP空间分布的影响差异是否有显著,结果表明气温和土地利用类型的影响力要强于其他因素。经交互探测结果表明(表5),2000—2020年6个影响因子交互作用后的影响作用能力更强,对NPP空间分异性的解释能力更高,且呈非线性增强。结合2000,2005,2010,2015和2020年5 a的双因子交互作用结果综合来看,影响力最强的交互因子均为该年份NPP第一影响力因子与其余影响因子的组合,即降雨量因子与其余因子交互作用之后的q值最高,其中2000年降雨量与气温双因子交互作用q值最高,5个代表年份的最强交互因子分别为降雨量∩气温和降雨量∩海拔交替主导,解释力分别为0.973 9,0.968 1,0.972 1,0.961 2和0.965 3,解释力均在95%以上。这说明西藏地区NPP的空间分异是多因素综合影响作用的结果。

表5 西藏地区2000—2020年NPP交互作用q值

3 讨论与结论

3.1 结论

(1)从时间上分析,2000—2020年西藏地区NPP年平均值为137.36 g/(m2·a),波动范围为128.04～148.06 g/(m2·a),2014年NPP均值最小,2006年为NPP最大值;研究区内NPP以0.04 g/(m2·a)的年变化速率增长,整体上呈波动上升趋势。

(2)从空间上分析,2000—2020年西藏地区年均NPP空间分布表现出较强的空间异质性,整体上呈东南高西北低的趋势。不同植被类型NPP均值差异明显,将NPP均值按大小顺序排列为:阔叶林＞栽培植被＞针叶林＞灌丛＞草甸＞高山植被＞草原＞荒漠。不同的地形条件,NPP空间分布差异明显,随海拔增加NPP呈先增加后减少趋势,随坡度增加NPP呈上升趋势。

(3)NPP变化的主要影响因素包括气候因子、人类活动和环境变化,其中气候变化是高原植被NPP变化的主要影响因素。利用NPP与气象因子进行相关性分析发现,西藏地区NPP与气象因子相关性显著,总体上NPP与降雨多呈负相关,与气温多呈正相关。

(4)西藏地区6个影响因子对NPP空间分异性的解释力排序依次为:降雨量＞海拔＞气温＞坡度＞植被类型＞土地利用类型,降雨量是西藏地区NPP影响力最强因子,自然因子对NPP的影响力强于人为因子。此外,双因子交互作用后的影响作用能力更强,对NPP空间分异性的解释能力更高,且呈非线性增强。

本研究对西藏地区NPP的时空格局特征及其影响因子进行了定量分析,直观地反映西藏地区NPP的变化规律,为区域生态环境评价、陆地生态系统固碳能力评估提供科学参考,为区域自然资源开发及绿色生态建设提供依据。

3.2 讨论

本文尚存在一些不足之处,未来需要进一步地改进,主要有以下几方面。

(1)本文所使用的MOD17A3数据及气象数据的空间分辨率过低,研究所覆盖的范围较广,在估值上容易出现较大的误差,且MOD17A3数据在空间分布上地类特征不明显,未来可使用精度更高的数据基于模型对区域陆表生态系统NPP进行估算,争取优化结果;

(2)通过气候因子与NPP变化曲线可以发现,气候因子年际变化趋势与NPP并不吻合,这就表明气候因子变化对高原植被生长的影响可能存在滞后或累加效应,且除降雨和气温外还有其他影响因素造成高原植被NPP的变化。经地理探测器结果验证了这一想法,但本文在进行影响因子的确定时,人为因子的选择过少,可能会对结果造成一定的偏差,且未对土地利用覆被的变化、人类活动及相关政府政策的实施与NPP变化之间的响应关系进行展开。因此,综合影响因素与NPP变化的响应机制将是下一步研究的重点。