2001—2018年中国总初级生产力时空变化的遥感研究

张心竹,王鹤松,延 昊,艾金龙

1 北京林业大学林学院, 北京 100083 2 北京林业大学生态与自然保护区学院, 北京 100083 3 国家气象中心, 北京 100081 4 益阳职业技术学院,生物与信息工程系, 益阳 413049

植被总初级生产力GPP(Gross Primary Production)是绿色植物通过光合作用吸收大气CO2制造的有机物总量,驱动着大气CO2含量的季节和年际变化,体现了陆地生态系统在自然条件下的生产能力[1],在全球碳循环中扮演着重要角色[2]。陆地生态系统在全球碳循环中作为一个主要的碳汇,能抵消大量的人为碳排放,所以准确量化GPP对生态系统功能评估和碳平衡研究至关重要,这已成为地球系统科学领域的研究热点[3]。

量化GPP的主要方法有地面观测和模型模拟两种。地面观测法多采用涡度相关技术,通过仪器直接获得相关源区的碳通量,经处理后得到GPP[4]。由于通量站点数量少,分布稀疏,所以地面观测数据常用于GPP模型的验证分析[5]。在区域乃至全球尺度上,生态系统过程模型和光能利用率模型是模拟GPP的有效方法[6]。过程模型考虑植物生理生态过程,理论意义强,并且能够模拟未来生产力变化。然而过程模型拥有不足,其结构复杂,包含大量理想化假设,所需参数较多,尤其是在大尺度应用上,模型精度受到一定影响[7]。

由于卫星遥感技术可以稳定、持续地获取植被的大尺度空间动态变化信息[8],多种遥感植被数据被广泛用于植被生产力估算。基于一层大叶模型假设的光能利用率(LUE)理论,建立了大量遥感光能利用率模型,如MODIS GPP标准产品,VPM模型、EC-LUE模型等[9- 11]。但不同模型对植被光能利用率的参数化方案存在较大差异。如MODIS GPP存在一定程度的低估,尤其是在农田地区,与其算法中没有区分C3、C4植物的最大光能利用率参数有关[5,12]。为此,Yan[13]提出了区分C3、C4植物光能利用率的TEC模型。由于植物冠层对散射辐射较直射辐射具有更高的光能利用率,He[14]提出了区分阴生叶和阳生叶的双叶光能利用率模型。Yan[15]进一步提出改进的双叶光能利用率模型(DTEC),考虑了直接辐射和散射辐射及C3和C4植物差异的影响。

日光诱导叶绿素荧光(solar-induced chlorophyll fluorescence,SIF)是绿色植物在自然光照下吸收太阳能进行光合作用时重新释放的红光和近红外光子,是植被瞬时光合作用活动的理想探针,在植物生长监测中有着巨大的应用潜力[16]。大量研究表明SIF与通量站点GPP显著正相关[17-18]。近年来,GOSAT、OCO- 2、GOME- 2等卫星均能反演植被SIF信息,为大尺度植被GPP研究带来了新的思路和方法[19-20]。但SIF卫星数据分辨率较粗,对GPP的大尺度模拟和验证造成了困难。为解决这一问题,Li和Xiao[21]基于OCO- 2 SIF和气象数据使用立方体集成法(Cubist),开发了时间分辨率为8 d,空间分辨率为0.05°的全球SIF数据集(GOSIF)。

应用各种大叶遥感光能利用率模型在中国已开展了大量研究,在GPP年际变化方面存在较大差异[6,22-26],而采用双叶遥感GPP模型和基于日光诱导叶绿素荧光SIF的GPP数据集的中国区域研究还很少。为了更好地准确量化和分析中国区域植被GPP时空动态,本文使用GOSIF-GPP和DTEC-GPP对中国2001至2018年GPP进行模拟,着重研究：(1)中国陆地生态系统总初级生产力的时空格局及年际变化；(2)采用中国通量网的6个通量站点开展GPP模型验证实验；(3)比较两种模型GPP的差异。

1 材料与方法

1.1 GOSIF-GPP产品

Li和Xiao[27]基于GOSIF数据和通量站观测GPP建立起八种不同植被类型的由SIF反演GPP的统计模型,进而由GOSIF数据反演得到2001—2018年全球GPP数据集。本研究使用的是2001—2018年中国区域的GOSIF GPP月尺度产品(http://globalecology.unh.edu)。

1.2 DTEC-GPP模型

DTEC模型是气象和遥感数据驱动的双叶光能利用率模型,区分了阴生叶、阳生叶对GPP的贡献和直射辐射、散射辐射影响[14],并采用了一个改进的水分胁迫系数[13],模型主要结构如下：

GPP=(εmsu×APARsun+εmsh×APARshd)×Wε×Tε

(1)

式中,εmsu和εmsh分别是阳生叶和阴生叶的最大光能利用率,APARsun和APARshd分别是被阳生叶和阴生叶吸收的光合有效辐射[14],Wε和Tε分别是水分和温度胁迫因子。对于C3植物：εmsh=3.78×Df1.8,εmsu=1.67 gC/MJ。对于C4植物：εmsh=5.78×Df1.8,εmsu=2.56 gC/MJ。

(2)

(3)

(4)

LAIshd=LAI-LAIsun

(5)

PARdif=PAR×Df

(6)

式中,PARdir和PARdif分别是直射光合有效辐射和散射光合有效辐射,LAIsun和LAIshd分别是阳生叶和阴生叶的叶面积指数,PAR是光合有效辐射,Df是散射比[15],其余参数的含义及计算方法见文献[15]。

水分胁迫因子Wε的定义来自于TEC模型[13]：

Wε=E/EPT

(7)

式中,EPT是潜在蒸发散[28],E是实际蒸散[29]。

1.3 气象数据

驱动DTEC-GPP模型的气象因子包括空气温度、水汽压、空气相对湿度、降水量、风速和日照时数。本研究使用了来自于国家气象信息中心(http://data.cma.cn)的中国区域2000个高密度台站2001—2018年的月气象数据,采用反距离加权(IDW)插值方法将2000个台站的气象数据插值到0.05 ×0.05°的空间分辨率。并根据联合国粮食及农业组织(FAO)的方法[30]由日照时数等气象数据计算出月地表总辐射和净辐射。

1.4 遥感数据

驱动DTEC-GPP模型的LAI遥感数据来自于MOD15A2H LAI/FPAR产品(https://lpdaac.usgs.gov/products/mod15a2hv006/),这是Terra卫星上搭载的中分辨率成像光谱仪(MODIS)数据反演的植被参数产品,时间分辨率是8d,空间分辨率是500 m。本文将通过质量控制后的LAI数据重采样到月尺度5 km空间分辨率,与月尺度的气象数据结合驱动DTEC GPP模型。

1.5 涡度通量数据

涡度协方差(EC)方法可以测量陆地生态系统和大气之间的CO2、水和能量通量,因此EC观测结果常被用于评估各种GPP模型[4]。本文用来验证GPP的通量数据来自于中国通量网(http://chinaflux.org)提供的六个通量站的GPP数据,包括森林、草地、农田三种生态系统类型,具体站点信息见表1[31]。

表1 中国通量网六个站点的简要信息

Table 1 The brief description of six ux sites of Chinaflux

表1 中国通量网六个站点的简要信息

站点名称Site name缩写Abbreviation数据年份Year经度/(E)Latitude纬度/(N)Longitude海拔/mAltitude植被类型Biome type优势种Dominant species 长白山站CBS2003—2008128°05'45″42°24'9″738温带针阔混交林红松、椴树千烟州站QYZ2003—2008115°03'29.2″26°44'29.1″100亚热带常绿针叶林马尾松、湿地松海北站HB2003—2008101°19'52″37°39'55″3358高寒金露梅灌丛草甸金露梅、矮蓑草当雄站DX2004—200891°03'58.90″30°29'50.21″4333草原化高寒草甸高山蓑草内蒙古站NMG2004—2008116°24'14.4″43°19'31.8″1200温带羊草草原羊草、冰草禹城站YC2003—2008116°34'12.72″36°49'44.4″28农田冬小麦、夏玉米

1.6 研究方法

对2001—2018年的DTEC和GOSIF GPP年数据,进行时空变化分析。采用一元线性回归法在全国尺度和象元尺度分析18年间GPP值的趋势倾向率,具体采用线性拟合方程(y=k×x+b)的斜率k表征18年间GPP的变化趋势和幅度。k>0表示GPP为增长趋势；k<0表示GPP为减小趋势。

采用六个通量站观测的GPP数据对DTEC和GOSIF估计的GPP进行了月尺度验证,使用决定系数(R2)、均方根误差(RMSE)和偏差(Bias)进行精度评估。

2 结果

2.1 GPP模拟验证

图1显示DTEC模型和GOSIF模型都能较好模拟通量站点GPP的季节动态。在长白山站点,DTEC和GOSIF模型表现十分优异,有与站点GPP几乎相同的季节动态。在千烟州站,二者在冬季均有一定程度的低估,这可能是由于千烟州站地面异质性大,空间代表性不佳。在当雄和海北草地站点上,两种GPP模型在夏季生长季都存在高估,尤其是GOSIF GPP,Xiao和 Li[27]也指出GOSIF GPP在生产力偏低的地区存在一定程度的高估。禹城农田站由于种植冬小麦和夏玉米,一年有两个生长高峰,DTEC和GOSIF模型均能捕捉到两个生长高峰但同时也均存在低估。值得注意的是,遥感象元数据与地面通量塔观测数据存在尺度差异,而且地面的高异质性也会使通量站数据的代表性较差,这都会影响验证效果。

图1 月尺度上模拟的GOSIF GPP、DTEC GPP和通量观测GPP的季节动态Fig.1 Seasonal variations of monthly estimated GOSIF GPP, DTEC GPP and flux tower-observed GPP at six flux tower sitesGPP: 总初级生产车Gross primary productivity

对两种GPP模拟结果的验证统计见表2,GPP模拟精度随着通量站点和生态系统类型变化。在森林站点的验证结果最好,R2为0.88—0.95。在草地站点中,海北站、当雄站的精度较高,R2为0.83—0.91,而内蒙古站R2较低,小于0.6。在农田站点,GOSIF的R2明显高于DTEC,但两者的RMSE和 Bias均相对较大,存在明显低估。GOSIF GPP在3个草地站均存在高估,在森林站和农田存在低估,而DTEC GPP仅在当雄站高估,在其余5个站点均为低估。此外,在各个站点GOSIF GPP精度略高于DTEC GPP。

所有站点的总体验证结果(图2)表明,GOSIF GPP的精度(R2=0.78; RMSE=54.67 gC m-2月-1; Bias=-5.15 gC m-2月-1)略好于DTEC模型(R2=0.74; RMSE=66.07 gC m-2月-1;bias=-30.35 gC m-2月-1),DTEC低估程度更大。王克清[32]在与本文相同的六个通量站点处对两种GPP模型进行验证,R2分别在0.72—0.89与0.67—0.94之间。Jia等人[33]在8个草地站点对8个模型进行验证,R2在0.64—0.89之间。Li[24]在中国32个站点对EC-LUE模型进行验证,C3站点和C4站点的R2分别为0.79和0.62。与上述研究相比,本文的两种模型精度与之相近。

2.2 中国主要生态系统类型的GPP年际变化特征

2001—2018年主要生态系统类型GPP年值变化曲线见图3。GOSIF GPP和DTEC GPP模拟的不同生态系统年均值大小顺序相同,常绿阔叶林的GPP年均值最高,分别为2140和2014 gC m-2a-1,其次是常绿针叶林和混交林,分别是1713、1717 gC m-2a-1和1464、1452 gC m-2a-1,然后是落叶阔叶林、农田、灌丛、落叶针叶林、草地。草地生态系统平均生产力最低,两种模型模拟结果均不到300 gC m-2a-1。此外,GOSIF GPP除在常绿针叶林生态系统稍低于DTEC GPP,在其余生态系统类型上均高于DTEC GPP。

从变化趋势上看,两个模型模拟的主要生态系统GPP年值从2001至2018年都呈显著上升趋势(P<0.01),仅GOSIF GPP在落叶针叶林生态系统中的上升趋势不显著。两种模型的波动幅度有差异,GOSIF GPP波动较小,呈平稳上升,而DTEC GPP波动幅度较大。

表2 模拟的月尺度GOSIF GPP、DTEC GPP在6个通量站点的验证统计

图2 GOSIF GPP, DTEC GPP与通量站点GPP月数据的比较Fig.2 Scatter plot of GOSIF GPP and DTEC GPP vs. observed GPP at six flux tower sites on a monthly scale

图3 2001—2018年中国不同植被类型 GPP 年际变化Fig.3 Interannual variations of GPP of different vegetation types from 2001 to 2018

2.3 中国陆地生态系统GPP时空变化特征

2.3.1中国陆地生态系统GPP多年空间分布格局

从2001—2018年GPP多年均值空间分布(图4)来看,GOSIF GPP和DTEC GPP的空间分布格局基本相同,都展示了巨大的空间差异性,GPP高值主要分布在东南地区,低值主要分布在西北地区,总体上呈现由东南向西北、从沿海向内陆递减的分布趋势,这与已有研究结果[6,23-24,26]基本一致。

其中华南地区、华东地区、云南中南部及青藏高原南部GPP多年均值高于2000 gC m-2a-1,西南地区、长江中下游地区GPP多年均值为1500 gC m-2a-1左右。南方地区生产力高主要是由于气候温暖湿润、光照充足和植被覆盖度高。此外,东北森林地区也具有较高的GPP年均值,约为1000 gC m-2a-1。而西北的沙漠和荒漠地区、青藏高原中北部以及内蒙古中西部的草原区,年均GPP不足500 gC m-2a-1,这些地区由于水分胁迫和低温制约,植物生长季偏短,导致了生产力低下。内蒙东部及青藏高原东南部为过渡地带,GPP多在500 gC m-2a-1左右,与我国半干旱半湿润分界线基本一致,说明降水是影响我国植被生产力分布的一个重要因素。华北、黄淮大部分地区属于农业区,GPP多年均值在800—1500 gC m-2a-1之间。此外,新疆天山山脉附近GPP能达到500 gC m-2a-1左右,这是由于来自大西洋的盛行西风和来自北冰洋的气流,碰到山坡抬升,产生较强降水,促使天山西部及西北部具有较高的植被覆盖度[34],同时天山的高山融雪也会促使山前地带绿洲的形成。

图4 2001—2018年中国年均GPP空间分布格局Fig.4 Spatial pattern of annual GPP in China during 2001 to 2018

图5 2001—2018年中国植被 GPP 年际变化Fig.5 Interannual variations of GPP in China from 2001 to 2018

2.3.22001—2018年中国陆地生态系统GPP总量年际变化

2001—2018年GOSIF GPP和DTEC GPP的年际变化如图5,总体均呈显著增加趋势,但年际间有波动,其中DTEC GPP变化波动更剧烈。使用GOSIF数据得到的GPP年总量变化范围为6.35—8.08 PgC/a,均值为7.23 PgC/a,年增幅约为0.094 PgC/a(P<0.01)。 DTEC模型得到的GPP年总量变化范围为6.20—7.79 PgC/a,均值为6.93 PgC/a,年增幅约为0.073 PgC/a (P<0.01)。无论是GPP总量还是增长幅度,GOSIF GPP均大于DTEC GPP。DTEC GPP的年值曲线在2003、2005、2009、2014年分别出现过较大的波动,这与刘刚等[35]模拟的2001—2014年中国植被净第一性生产力的结果类似,在2003、2005、2009、2014年都出现GPP低值,主要是这些年份遭受到不同程度的气象灾害。本研究期间,2001—2018年,中国GPP显著上升,这是由适宜的气候变化、CO2的施肥效应、氮沉降、人类活动(如植树造林、农业灌溉)等多方面因素共同作用的结果[1,36-38]。

2.3.3中国陆地生态系统GPP多年时空变化特征

2001—2018年我国GPP时空变化见图6,近18年间GOSIF GPP和DTEC GPP年变化趋势基本相同,上升趋势和下降趋势并存,但呈上升趋势的面积更大。变化趋势率k值总体上呈东南地区大于西北地区、沿海地区大于内陆地区的分布。年GPP呈上升趋势最高的地区包括华南、华中、华东地区、西南地区东南部以及海南岛,趋势率k能够达到20—30 gC m-2a-1,其次为陕西、山西、河北、东北大部地区,上升幅度约15 gC m-2a-1左右。此外,GOSIF GPP数据结果显示新疆天山山脉附近区域GPP也有明显的上升趋势,个别地区上升幅度甚至能达到20 gC m-2a-1。上述地区的GPP增长趋势与已有的研究基本一致[22,35,39]。西北地区、西藏、内蒙古中东部变化幅度较小,基本在-5—5 gC m-2a-1之间。GPP下降趋势比较严重的地区则分布在青藏高原南部、云南中北部和上海、天津等城市区域,下降幅度约为-10 gC m-2a-1。此外,DTEC GPP显示河南南部地区有较明显的下降趋势,大约在-20 — -5 gC m-2a-1之间。青藏高原南部、云南中北部和河南南部地区主要是由于降雨偏少造成干旱事件引发GPP下降[25,40]。而上海、天津等城市GPP的下降趋势是城市建设用地扩张造成植被生产力下降。总体上,我国长江以南大部分地区GOSIF GPP的年变化趋势率k高于DTEC GPP。

图6 2001—2018年中国年GPP变化趋势Fig.6 Trend of annual GPP in China during 2001 to 2018

3 讨论

3.1 中国GPP年际变化

在已有的中国GPP研究中,光能利用率模型、过程模型和以机器学习法为主的统计模型均有广泛应用。由于时间段选取、模型选择、参数确定、数据集来源不同以及各种误差,GPP模拟结果也有较大差异(表3)。其中多年平均GPP最大值为使用支持向量回归方法[41]得到的2000—2015年间的7.81 PgC/a,最小值为使用EC-LUE模型[42]得到的的2000—2009年间的5.38 PgC/a。本研究使用GOSIF和DTEC双叶模型得到的GPP年总量分别为7.23PgC/a和 6.93 PgC/a,均处在已有结果的范围之内。采用多种类型的GPP模型从多角度开展大尺度GPP研究,可以减小陆地生态系统碳循环研究的不确定性。

从年际变化趋势上看,本研究发现中国植被生产力从2001至2018年呈上升趋势,与已有基于卫星观测和生态系统模型的研究结论基本一致[46]。但本研究中GOSIF和DTEC GPP上升趋势分别为0.094 PgC/a和0.073 PgC/a,GOSIF GPP上升幅度更大。而已有研究得到的GPP上升趋势为0.02 PgC/a至0.057 PgC/a(表3),低估了中国GPP增长的趋势。有研究表明,目前对GPP年际上升趋势的估算存在明显低估的现象[6,43],与本研究结论一致。本文选取的GOSIF算法和DTEC模型是目前最新的研究方法,在机理上都具有深刻意义,以中国地区为研究对象也是首次。

表3 不同模型估计的中国年均GPP

3.2 GOSIF算法与DTEC模型对比

GOSIF算法与DTEC模型是两类不同的模型,在不同方面拥有优势。从机理上看,日光诱导叶绿素荧光SIF包含大量光合信息,能准确地反映植物光合作用,与植物生理过程高度耦合[16],在叶片、冠层、植株、生态系统等不同尺度上均与GPP显著相关[17-18]。GOSIF是目前时间和空间分辨率最高的SIF全球产品,GOSIF GPP根据植物叶绿素进行光合作用时释放的SIF辐射直接估计GPP,已经包含各种环境胁迫和人为因素的影响[21,27,47],能够更加客观的反应植物的真实生长状态,捕捉更多空间分布上的细节。但是GOSIF GPP在生产力较低地区出现高估,精度受统计模型和训练数据的影响。DTEC GPP根据温度和水分胁迫下植物对光合有效辐射的利用率进行GPP估计,考虑了散射辐射的作用,采用改进后的水分胁迫参数,区分了C3、C4植物参数化,在冠层复杂的森林和干旱少雨的地区有良好表现[15],较传统大叶LUE模型有较大改进。但由于模型输入参数较多,并且没有考虑人为灌溉和西部地区高山融雪的影响,可能会引入误差。例如DTEC GPP主要是通过降雨量驱动土壤水分平衡方程估计水分胁迫,这导致DTEC GPP在人为经营程度高的地区,如河南和河北的农田、新疆天山附近的绿洲和一些集约化经营管理的林地会低估GPP。而GOSIF GPP包含了自然降雨和人工灌溉的影响,能较客观的反映出自然和人为影响下的植被生产力。因此,使用这两类GPP模型有助于理解自然和人为因素对植被生产力的不同影响。

4 结论

DTEC和GOSIF模型都揭示出2001至2018年中国陆地植被GPP呈显著增加趋势,平均每年分别增加0.073Pg C和0.094Pg C,在植被稠密的东南地区GPP增幅高于植被稀疏的西北地区。中国DTEC和GOSIF GPP多年平均值分别为6.93 Pg C和7.23 Pg C,GPP空间格局呈现出明显异质性,从南到北、从东到西GPP年均值逐渐减小。在6个通量站的验证表明DTEC和GOSIF模型表现良好,具有可用性。

比较而言,GOSIF GPP算法能客观地反映植被生产力状况,包括人为影响严重的灌溉农田和绿洲地区；而DTEC模型则更适合自然条件下植被生产力的模拟,研究光、温、水等气象因子对植被生产力的影响。未来仍需要改进GPP估计模型和方法,采用多种GPP模型研究方法,探索气候和人为因素对植被生产力的影响以及植被对气候变化的反馈。