遥感GPP模型在中亚干旱区4个典型生态系统的适用性评价

许世贤,井长青,*,高胜寒,邬昌林

1 新疆农业大学新疆草地资源与生态重点实验室,乌鲁木齐 830052 2 中国科学院水利部成都山地灾害与环境研究所,成都 610041

总初级生产力(Gross Primary Productivity,GPP)是评估陆地生态系统碳通量的重要指标,定义为植被通过光合作用固定的有机碳量[1]。准确估算GPP有助于理解大气和陆地生态系统的碳循环过程。涡度通量技术(Eddy Covariance,EC)可以在站点尺度监测地面生态系统的二氧化碳交换量,从而准确估算GPP[2],但涡度通量观测得到的数据所代表的空间尺度较小,受空间异质性及下垫面复杂的影响,难以将估算结果外推到更大尺度[3]。相较于站点观测二氧化碳通量,遥感模型具有长时间序列大区域范围监测的优势[4]。中亚干旱、半干旱区是全球八大干旱区之一,面积广阔[5],其碳通量是全球碳循环的重要组成部分,而区域内生态系统种类较多且易受全球气候变化干扰,对外界环境因素的变化有很快速的响应[6]。比较遥感GPP模型在区域内不同生态系统上的适用性,有助于更准确的评估干旱区生态系统的生产力和理解全球生态系统碳循环。

目前遥感GPP估算模型主要包括光能利用率模型、光合作用机理过程模型、遥感植被指数模型。过程模型的代表有BEPS模型[7]、BESS模型等[8],模型从植被吸收二氧化碳的机理过程出发,利用算法模拟,具有较高的精度,但模型参数众多且结构复杂,精度依赖于模型结构和众多输入参数的可靠性。光能利用率模型基于光能转化率理论[9],理论认为植被实际光合作用能力的影响因素为光合有效辐射量(PAR)、植被光能最大吸收效率及其环境压力对植被光能吸收效率的限制。光能利用率模型考虑了光能吸收机理的同时模型结构较为简单,因此得到了很大程度的发展。CASA(Carnegie Ames Stanford Approach)[10]、MOD17模型[11]、EC-LUE(Eddy Covariance Light Use Efficiency)[12]、VPM(Vegetation Photosynthesis Model)[13]等光能利用率模型被建立并得到了广泛应用[14]。植被指数模型同样根据光能转化率理论建立而来,主要通过遥感植被指数代替光能吸收效率,从而实现模型的简化,通常植被指数模型相较于光能利用率模型输入参数更少,模型结构也更简单。Sims等发现增强型植被指数(EVI)与GPP有很强的相关性[15],之后在此发现的基础上引入对植被生长有重要影响的因子温度,进而建立起温度与绿度模型(TG,Temperature And Greenness Model)[16]；2010年Wu提出了只依赖EVI与PAR的VI模型(vegetation index model)[17]；SANIRv模型考虑到荒漠植被植被覆盖度低的特点,建立了植被近红外反射率指数(NIRv)和PAR与GPP的函数关系[18]。目前,也有利用微波遥感[19]或日光诱导植被荧光(solar-induced fluorescence,SIF)观测技术建立的GPP模型出现[20]。

多种遥感模型在全球或区域尺度上被建立应用,但不同遥感模型在同一区域内得出的结果可能存在较大差异,因此在应用前需要评估模型在区域的适用性。例如Gebremichael等人在墨西哥半干旱区对MODIS GPP产品进行了验证[21]。Chen等人对比三种遥感模型在中国不同气候区的区域差异[22]。Wang等人通过中国华北10个通量站点数据验证MOD17模型的精度[23]。由于中亚地区地理位置和气候条件的特殊性,不能简单的将其他区域的模型评价结果“迁移”至中亚地区,因此需要针对中亚干旱半干旱区评估遥感GPP模型的适用性。

本研究以中亚四个典型生态系统涡度相关观测站点作为研究对象,选取遥感GPP模型中具有代表性的光能利用率模型：MOD17模型、VPM模型；植被指数模型：TG模型、SANIRv模型；微波遥感模型：VODCA2模型,将基于EC技术观测的通量数据作为模型验证数据,比较五种遥感GPP模型在中亚四个典型地表生态系统的适用性,为遥感GPP模型的发展及其在干旱区的应用提供理论依据。

1 材料与方法

1.1 研究区概况

中国新疆阜康荒漠生态国家野外科学观测研究站(简称CN-Fuk)位于古尔班通古特沙漠南缘绿洲(44°17′N,87°56′E)[24],属于典型荒漠生态系统,通量观测站点属于中国通量观测联盟(ChinaFLUX),地表优势种群为荒漠灌木怪柳科植物以及旱生草本植物,植被盖度在20%—30%。中国乌兰乌苏绿洲农田生态与农业气象试验站(简称CN-Wul),位于新疆天山北麓中段、准噶尔盆地南缘(44°17′N,85°51′E)。属于典型农田生态系统,农作物主要为棉花、玉米和冬小麦等。土壤类型为沙壤土,肥力中等[25]。哈萨克斯坦咸海站(简称KZ-Ara)位于咸海东北部(44°58′N,61°05′E)[26],为典型荒漠生态系统,地表优势物种为为灌木、柽柳等,并伴生多种旱生草本植物、苦豆子、芦苇等。植被覆盖度约为30%。哈萨克斯坦巴尔喀什湖站(简称KZ-Bal)位于巴尔喀什湖与卡普恰盖水库之间(44°34′N,76°39′E)[26],属于典型草地生态系统,地表覆盖类型为草地,并伴生有少部分的农田与灌木,植被覆盖度约为50%(图1)。本研究所采用的四个通量站点均位于中纬度欧亚大陆腹地,均属于温带大陆性气候,夏季炎热,冬季寒冷,年降雨量少且主要集中在夏季。四个站点地形起伏均较小。咸海站和阜康站气候干旱,年均温在5.7—8.3 ℃之间,降雨稀少,年降雨量在140—150 mm之间[27],属于典型的干旱区荒漠生态系统；巴尔喀什湖站气候条件与咸海站、阜康站类似,但因其毗邻河流和绿洲农田,地表水分条件较荒漠站更好；乌兰乌苏站年均温在7℃左右,年降雨量在210 mm左右,农作物生长所需水分主要依赖人工灌溉。CN-Fuk、KZ-Ara、KZ-Bal三个站点数据来源于中国科学院新疆生态与地理研究所,CN-Wul数据来源于乌兰乌苏绿洲农田生态与农业气象试验站。

图1 通量站点分布及土地利用类型图Fig.1 Flux station distribution and land use map

1.2 数据来源与预处理

1.2.1 遥感数据

本研究所采用的MODIS遥感数据有2005—2006、2009—2010、2012年的h23v04和h24v04两景的MOD09A1 8d 500m地表反射率合成产品、MOD11A2 8d 1km地表温度合成产品。以上数据均来源于美国国家航空航天局的分布式档案中心(https://ladsweb.modaps.eosdis.nasa.gov/)。首先采用MRT工具将所需数据从HDF数据集中提出并进行批量合成、重投影、重采样等处理,再根据数据质量控制文件,选取可靠像元,并插补缺失及不可靠像元,最后对计算得到的植被指数进行平滑处理,以消除云、气溶胶等因素的影响。

本研究遥感GPP模型的输入参数有：增强型植被指数EVI、归一化植被指数NDVI和陆地水分指数LSWI使用MOD09A1地表反射率产品计算获得,各植被指数的计算公式为：

(1)

(2)

(3)

式中,RED为红波段(620—670 mm)、NIR表示近红外波段(841—875 mm)、BLUE为蓝波段(459—479 mm)、SWIR为短波红外(841—875 mm)。

1.2.2 气象数据

本研究所采用的气温数据来自欧洲中期天气预报中心(ECMWF)的ERA5-land气象再分析数据,使用MATLAB软件提取站点坐标周边3 km×3 km地表2 m的气温数据。光合有效辐射(PAR)数据来自全球陆表特征参量产品数据集(http://www.glass.umd.edu/)[28]。使用MATLAB软件将获取到的HDF格式数据输出为tif格式,并提取站点周边3 km×3 km的数据。

1.2.3 GPP产品数据集

MOD17模型属于光能利用率模型,该模型结构可以表述为：

GPPmodis=0.45×εmax×m(Tmin)×m(VPD)×fPAR×SWrad

(4)

式中,εmax表示为最大光能利用率(LUE),m(VPD)与m(Tmin)分别为饱和水气压差和低温条件降低εmax的标量,fPAR为PAR的分量,SWrad为太阳短波辐射分量。MOD17全球GPP数据可由MOD17产品获得,来源为美国国家航空航天局的分布式档案中心(https://ladsweb.modaps.eosdis.nasa.gov/)。

VODCA2模型是基于被动微波遥感估算GPP的方法。模型的主要参数为微波遥感植被光学深度(VOD),植被光学深度主要代表了植被含水量,模型通过反演植被含水量模拟GPP,模型原理详见[29],模型数据可由期刊Earth System Science Data获得[19]。模型主要结构可表述为：

GPPVOD=sVOD+s(Δ(VOD))+s(mdn(VOD))

(5)

1.2.4 GPP验证数据

CN-Fuk站的GPP验证数据来自于刘冉发表在期刊上的2005—2006年在阜康站采用涡度相关系统测得的5—10月的月值数据[30]。KZ-Ara、KZ-Bal、CN-Wul三个站的GPP验证数据来源于涡度相关系统观测(表1)。3个站的涡度相关系统主要由开路式CO2/H2O分析仪(LI-7500, LI-COR, Lincoln, USA)、三维超声风速仪(CSAT3, Campbell Scientific, Logan, USA)和数据采集器(CR5000, Campbell Scientific, Logan, USA)组成。CO2/H2O分析仪与三维超声风速仪观测高度为2 m,采样频率为10 Hz, 分别用来测定CO2、水汽脉动和风速、温度脉动,并通过数据采集器CR5000 (Campbell Scientific, Logan, USA)记录和存储10 Hz的原始通量数据。EddyPro软件可以将生态系统净交换碳通量(NEE)数据按NEE=Re-GPP的关系拆分,夜间NEE等于生态系统呼吸(Re),利用夜间数据建立Re与温度的关系从而推导日间的Re,进而计算得到日间的GPP[31],为与遥感数据的时间步长进行匹配,计算每8天的累积GPP。

表1 涡度通量数据基本信息

1.3 研究方法

1.3.1 VPM模型

VPM模型为光能利用率模型,该模型将温度胁迫、水分胁迫、叶片成熟度作为限制光能利用效率的因子,分别使用Tscalar、Wscalar、Pscalar作为量化指标。该模型结构表述为：

GPPVPM=ε×FPARchl×PAR

(6)

ε=εmax×Tscalar×Wscalar×Pscalar

(7)

(8)

(9)

Pscalar=1全叶时期

(10)

(11)

FPARchl=a×EVI

(12)

式中,PAR为光合有效辐射,单位为MJ,Tmin、Tmax、Topt分别为植被光合作用活性最低温、最高温、最适温度,a的值为1。

1.3.2 TG模型

TG模型为植被指数模型,遥感驱动数据为增强型植被指数EVI和地表温度LST,公式为：

GPPTG=mTG×EVIscaled×LSTscaled

(13)

EVIscaled=EVI-0.1

(14)

(15)

式中,m为缩放因子,表示为GPP与温度-绿度之间的斜率,需要通过遥感数据与站点观测数据标定。

1.3.3 SANIRv模型

SANIRv模型基于归一化植被指数(NDVI)和近红外波段(NIR)与光合有效辐射(PAR)的GPP估算模型,该模型结构可表述为：

GPPSANIRv=c×PAR×SANIRv

(16)

(17)

NIRv=NDVI×NIR

(18)

式中,c是GPP与PAR和SANIRv的斜率,通过遥感数据与站点实测数据标定。NIRvmax为单个像元在时间序列上的最大值。NIRvsoil是对土壤背景的NIRv值的初始估计,假设NIRv值在0—0.1区间内,将此区间分为十个步长为0.01的区间,统计每个像元值在时间序列上落在每个区间内的次数,并计算落在频次最高及其相邻的两个区间的像元平均值,将像元值作为NIRvsoil(为去除生长季对土壤背景值的影响,只统计小于NIRv在时间序列上的平均值的数)[32]。

2 结果与分析

2.1 模型关键参数的标定

不同遥感GPP模型中存在着不同的敏感参数,而敏感参数的准确选取很大程度上影响着模型的模拟精度。敏感参数通常受植被类型或气候条件的影响,因此,在应用遥感GPP模型的模拟估算时,需要对敏感参数进行标定[33]。TG模型中的敏感参数m和SANIRv模型中的参数c为模型其他参数与GPP之间的斜率,参数值通过站点数据与模型输入数据进行标定。VPM模型中的敏感参数εmax是最大光能利用率,三个模型在不同站点上的敏感参数数值如表2所示：

表2 遥感模型关键参数

VPM模型中,最大光能利用率εmax的选取所需植被类型只有C3或C4植物之分,经前人根据FLUXNET各通量站数据对εmax的标定[34],当植被类型为C3和C4时,εmax分别为1.911 gC MJ-1和2.86 gC MJ-1。本研究的研究区选取了中亚干旱区内四个典型生态系统,而中亚地处温带,生长的植被大部分为C3植物,四个典型生态系统的优势种均为C3植物,因此四个站点的最大光能利用率确定为1.911 gC MJ-1。

2.2 模型在四个站点上的适用性

2.2.1 模型在荒漠生态系统的精度

图2为KZ-Ara荒漠站五种模型模拟值与观测值的拟合结果,通过图2和表3可以看出五种模型在咸海地区荒漠灌丛生态系统中对GPP的模拟精度均较差。五种模型的决定系数排序为MOD17>VODCA2>TG>VPM>SANIRv,精度最高的模型MOD17模型的R2仅为0.52,五种模型的RMSE均在1 gC m-2d-1左右。CN-Fuk站与KZ-Ara站的下垫面条件相似,同样为荒漠灌丛生态系统。在此区域内五种模型的R2排序为VPM>SANIRv>MOD17>VODCA2>TG,其中VPM模型的R2达到了0.80,五种模型的RMSE均在2 gC m-2d-1以下,其中VPM模型的RMSE为0.55 gC m-2d-1。光能利用率模型在模拟荒漠灌丛GPP的精度要优于其他种类的模型,同时模型对短时间(8天)的GPP变化不敏感,更适宜模拟长时间(月)的GPP变化。

图2 KZ-Ara站和CN-Fuk站实测GPP与模型模拟GPP散点图Fig.2 Scatterplot of measured GPP vs model simulated GPP at KZ-Ara and CN-Fuk stationsGPP：总初级生产力 Gross primary productivity;TG：温度与绿度模型；VPM：植被光合作用模型

2.2.2 模型在农田生态系统的精度

图3为CN-Wul农田站五种模型模拟值与观测值的拟合结果,通过图3和表3可以看出,五种模型对农田生态系统GPP的模拟精度均较好。R2排序为SANIRv>TG>VPM>MOD17>VODCA2,R2均在0.5以上,其中TG模型、VPM模型、SANIRv模型的R2达到了0.9以上；植被指数模型的RMSE均在1 gC m-2d-1左右,光能利用率模型和微波模型的RMSE均在2 gC m-2d-1左右,植被指数模型较其他种类模型精度更高。

2.2.3 模型在草地生态系统的精度

图4为KZ-Bal草地站五种模型模拟值与观测值的拟合结果,通过图4和表3可以看出,五种模型对草地生态系统GPP的R2排序为VPM>SANIRv>MOD17>TG>VODCA2,R2较高的两个模型SANIRv和VPM的R2均在0.8以上,而VODCA2模型的R2为0.03,植被指数模型的RMSE均在2 gC m-2d-1左右,而光能利用率模型和微波模型的RMSE均在5 gC m-2d-1以上。尽管VPM模型的R2达到了0.88,但是VPM模型对GPP出现了明显的低估。因此,在KZ-Bal站点区域对GPP模拟精度最高的模型为植被指数模型中的SANIRv模型。

图3 CN-Wul站实测GPP与模型模拟GPP散点图 Fig.3 Scatterplot of measured GPP vs model simulated GPP at CN-Wul stations

图4 KZ-Bal站实测GPP与模型模拟GPP散点图 Fig.4 Scatterplot of measured GPP vs model simulated GPP at KZ-Bal stations

表3 涡度观测GPP与模型模拟GPP的比较

2.3 模拟GPP与地面观测GPP的季节动态

2.3.1 荒漠生态系统的季节动态

图5和图6分别为KZ-Ara站和CN-Fuk站模型模拟值和观测值的季节动态,两个站点的GPP的峰值均在4 gC m-2d-1左右。两个站点上MOD17、SANIRv、TG、VPM模型均与EC GPP保持了较好的一致性,而VODCA2模型在生长季的开始前和结束后出现了峰值,对GPP的动态模拟能力最弱。

KZ-Ara站在生长季的观测期内,EC GPP的累计量为301.22 gC m-2,MOD17、VPM、TG、SANIRv、VODCA2模型的模拟值分别为268.87 gC m-2、332.15 gC m-2、284.71 gC m-2、285.52 gC m-2和240.03 gC m-2。SANIRv模型的相对误差最低。CN-Fuk站在两年的生长季观测期内,EC GPP的累积量为845.98 gC m-2,MOD17、VPM、TG、SANIRv、VODCA2模型的模拟值分别为516.49 gC m-2、815.11 gC m-2、542.57 gC m-2、567.32 gC m-2和864.34 gC m-2。VODCA2模型的相对误差最低(表4)。

图5 KZ-Ara站涡度通量GPP与模型模拟GPP季节动态Fig.5 Seasonal dynamics of EC GPP and model simulated GPP at KZ-Ara stationEC GPP表示涡度观测GPP,其他图例分别为各个模型模拟的GPP

图6 CN-Fuk站涡度通量GPP与模型模拟GPP季节动态Fig.6 Seasonal dynamics of EC GPP and model simulated GPP at CN-Fuk station

2.3.2 农田生态系统的季节动态

图7为CN-Wul站模型模拟值和观测值的季节动态,农田GPP的峰值在12 gC m-2d-1左右。SANIRv、TG模型与EC GPP有较好的一致性,MOD17、VPM模型出现了较大程度的低估,VODCA2模型在生长季的开始前就出现了峰值,与EC GPP的一致性最差。在两年的观测期内,EC GPP的累积量为1750.52 gC m-2,MOD17、VPM、TG、SANIRv、VODCA2模型的模拟值分别为1243.97 gC m-2、1244.03 gC m-2、1731.31 gC m-2、1796.29 gC m-2和1579.84 gC m-2,TG模型对CN-Wul农田站的模拟相对误差最低(表4)。

图7 CN-Wul站涡度通量GPP与模型模拟GPP季节动态Fig.7 Seasonal dynamics of EC GPP and model simulated GPP at CN-Wul station

2.3.3 草地生态系统的季节动态

图8为KZ-Bal站模型模拟值和观测值在时间序列上的特征,草地生态系统GPP的峰值在11 gC m-2d-1左右。TG、SANIRv模型与EC GPP有较好的一致性,MOD17、VODCA2、VPM模型则对GPP有较大程度的低估。观测期内,EC GPP的累积量为927.55 gC m-2,MOD17、VPM、TG、SANIRv、VODCA2模型的模拟值分别为206.87 gC m-2、350.32 gC m-2、909.18 gC m-2、948.19 gC m-2和211.10 gC m-2,MOD17、VPM、VODCA2模型的模拟结果对GPP的低估程度较大,TG模型对KZ-Bal草地站的模拟相对误差最低(表4)。

图8 KZ-Bal站涡度通量GPP与模型模拟GPP季节动态Fig.8 Seasonal dynamics of EC GPP and model simulated GPP at KZ-Bal station

表4 涡度观测GPP总量与模型模拟GPP总量的相对误差RE/%

3 讨论

本研究通过比较涡度协相关法观测的地面数据与五种遥感GPP模型模拟的GPP数据,分析遥感模型在4个典型生态系统上的适用性。将五种GPP模型的关键输入参数与GPP进行对比分析,可以发现在不同的站点上,GPP与PAR的决定系数均较高且大致相等,说明四个站点上GPP受PAR影响的程度较为一致,可以看出GPP并未受光照条件限制。在KZ-Ara和CN-Fuk两个荒漠植被站点上GPP与LST的决定系数较低,与气温的决定系数较高,这可能是因为荒漠植被地表覆盖度较低,地表温度较气温更高,不能很好的代表植被实际环境温度[35]；而CN-Wul农田站和KZ-Bal草地站植被覆盖较高,地表温度对植被实际环境温度的代表性较好,因此在这两个站点上GPP与LST的决定系数较高。除CN-Wul农田站外,其余三个站点植被指数NDVI、EVI与GPP的决定系数偏低,NDVI、EVI可以反映植被生长状况,但由于中亚地处干旱区,自然生态系统的植被覆盖度往往较低,受土壤背景等因素影响,植被指数会偏低,这可能会对GPP的模拟带来误差。KZ-Ara荒漠站LSWI与GPP的决定系数达到0.9以上(表5)。这说明荒漠植被GPP的主要限制因素为水分胁迫,而CN-Wul农田站LSWI与GPP的相关系数非常低,这可能是因为农田作物生长所需要的水分来自于浇灌,而遥感数据为瞬时值,难以捕捉八天内植被、土壤中水分的动态变化,这可能会对GPP的模拟估计造成误差[36]。

表5 涡度通量GPP与模型输入因子的决定系数R2

本研究比较的两种光能利用率模型MOD17、VPM模型对农田、草地的GPP有较为明显的低估,这可能是由于MOD17定量水分胁迫的指标为VPD,主要体现了大气水分亏缺量,而农田与草地站地处干旱区,空气较为干燥,植被生长所需的水主要来源于浇灌或地下水,这可能导致MOD17模型过分估计了水分胁迫的程度,最终导致模型对GPP的低估[37]；而VPM模型定量水分胁迫的指标为LSWI,基于前文所述,遥感数据可能对植被、土壤中水分的动态估计存在一定偏差,这可能导致了VPM模型对两个站点GPP的低估。以往的研究结果表明：VPM模型对各土地覆盖类型GPP的整体模拟精度较高,但对农田的GPP也存在一定程度的低估(～15%)[38],MOD17模型模拟玉米、大豆的GPP时也出现了低估[39]。王旭峰等人在黑河上游阿柔冻融观测站利用VPM模型高寒草甸GPP,结果表明VPM模型可以很好的模拟高寒草甸的GPP[40]。杨明兴等人验证了MOD17模型在亚热带针叶林和常绿阔叶林的GPP模型精度,结果表明MOD17模型对常绿阔叶林的模拟效果较好,而对针叶林的GPP有较大程度的低估,全年低估量达到了41.27%[41]。两种植被指数模型TG、SANIRv模型因其输入参数较少,只通过植被指数表征植被生长状况,而“固碳”能力通过站点数据与遥感数据进行标定,减少了输入数据不准确对结果的干扰。前人研究发现KZ-Bal草地生态系统的固碳能力较KZ-Ara荒漠生态系统的能力更高[27],这与本文对TG、SANIRv模型参数m和c的标定结果呈现的特征一致。微波遥感模型VODCA2模型在四个站点的模拟上精度均较低,这可能是由于数据集空间分辨率较低,导致观测数据与模型模拟数据之间存在尺度不匹配的问题。由于中亚干旱区内分布的通量站点较少,本研究每种植被类型只有一个通量站的观测数据,可能会使结果产生一定的偏差,有待于进一步研究。

4 结论

本研究针对中亚干旱区四个典型生态系统,将五种遥感GPP模型的模拟结果与地面观测结果进行对比分析,初步揭示了五种遥感GPP模型在中亚地区对GPP的区域化模拟的应用潜力,主要结论如下：

(1)综合三个评价指标,遥感GPP模型中对咸海和阜康荒漠灌丛生态系统的GPP模拟效果最好的模型是光能利用率模型；对中国乌兰乌苏农田生态系统模拟效果最好的模型是植被指数模型中的TG模型；对哈萨克斯坦巴尔喀什湖草地生态系统模拟效果最好的模型是植被指数模型中的SANIRv模型。

(2)对比中亚干旱区三种生态系统GPP的峰值发现,农田(12 gC m-2d-1)>草地(11 gC m-2d-1)>荒漠植被(4 gC m-2d-1)。四个涡度通量观测站点地表植被GPP与LSWI的决定系数在荒漠植被上较高而在农田、草地植被上较低,说明水分亏缺是干旱区植被GPP的主要限制因素。