植被净初级生产力(Net Primary Productivity, NPP)是陆地生态系统碳循环的重要组成, 也是揭示碳源/汇的关键环节, NPP估算不仅可以反映植被的物质生产能力, 也是揭示大气二氧化碳收支和气候环境变化的重要指标^[1-3]。

NPP估算模型较多, CASA模型其原型由Monteith提出, 之后Potter和Field等^[4]在Monteith基础上改进, 因结构简单, 参数可通过遥感获得, 广泛应用于全球和区域尺度上的NPP估算^[5-7]。CASA模型中需要输入气象数据, 气象数据一般来源于气象站点, 而气象站点通常选择在空旷平坦的地带, 因此, 通过站点数据直接空间插值的气象数据较难反应复杂地形地区的气候状况, 在估算NPP时会产生较大的误差。基于CASA模型估算NPP与气候因子相关的参数有太阳辐射、温度胁迫系数和水分胁迫系数。太阳辐射校正方法比较成熟^[8]；温度胁迫因子是通过气温得到的, 大部分研究使用如Kriging、IDW和Spline等方法对气温进行空间插值, 效果较差^[9-10], 有研究表明, 在甘肃省海拔、经度和纬度对气温影响也很大^[11]；水分胁迫系数在CASA模型中是实际蒸散量与潜在蒸散量的比值, 这两个变量是根据经验公式计算, 没有考虑到地表植被差异, 即相同的比值对不同植被类型的水分胁迫是一样的, 显然, 这种计算存在较大的不确定性。并且, 原CASA模型中全球植被的最大光能利用率为0.389 gC/MJ^[12], 对于干旱区并不适用。目前, 有很多学者都对CASA模型进行了校正^[13-15], 但对模型校正前后NPP估算的影响却鲜有研究。

河西走廊位于中国西北地区, 属于干旱半干旱区, 是国家丝绸之路经济带的重要通道, 祁连山国家森林公园也位于此, 区内景观多样, 在自然环境和人类活动双重作用下, 该区生态系统非常脆弱, 植被退化非常严重, 荒漠化发展速度较快^[16]。本文利用校正的CASA模型, 以气象数据、遥感数据、DEM数据以及土地利用数据等为基础, 针对CASA模型的以上缺点, 利用海拔、经度和纬度对气温进行校正, 并使用地表水分指数(Land Surface Water Index, LSWI)反演水分胁迫系数, 最大光能利用率采用冯益明等^[17]针对干旱区得到的不同植被类型上的相应值, 估算了2015年河西走廊NPP。旨在提高基于CASA模型的NPP模拟精度, 并着重分析校正前后的CASA模型对NPP估算的差异及原因, 同时也为该区植物生产力评价, 生态资产核算, 以及生态补偿机制建立提供科学参考。

1 研究区概况及数据 1.1 研究区概况

河西走廊地处甘肃省境内, 位于黄河以西, 南部被祁连山和阿尔金山、北部被合黎山、龙首山以及马鬃山夹持, 是古丝绸之路的交通要道, 地理位置为37°17′—42°48′ N, 92°12′—103°48′ E, 行政区划包括武威市、张掖市、金昌市、酒泉市和嘉峪关市(图 1)。河西走廊地势南高北低, 南部山地海拔大多在3000—3500 m以上；中部走廊平原海拔为1000—2600 m左右；北部山地断续分布。河西走廊属于温带大陆性气候, 光照强烈, 年降水量大部分地区在35—200 mm, 大部分地区的年蒸发量在1500 mm以上, 年均温5.8—9.3℃。该区域内植被类型主要有林地、耕地、草地和荒漠。河西走廊由三大流域组成, 从东向西依次为石羊河流域、黑河流域和疏勒河流域, 三大河流都发源于祁连山, 区域内绿洲农业发达。

1.2 数据来源与预处理 1.2.1 遥感数据

本文2015年的遥感数据来自美国国家航空航天局(https://search.earthdata.nasa.gov/search), 产品是MOD13Q1(从中选取归一化植被指数——NDVI(Normalized Difference Vegetation Index), 为16日合成数据, 空间分辨率250 m)、MOD09A1(8天合成的地表反射率数据, 选取近红外波段(841—875 nm)和短波红外波段(1628—1652 nm), 空间分辨率为500 m)和MOD17A3(选取NPP, 分辨率为1 km, 时间分辨率为1年)。利用MRT批处理工具对遥感数据进行了拼接、格式转换和重投影, 并在ArcGIS中完成掩膜, 使用最大合成法(Maximum Value Composite, MVC)将MOD13Q1和MOD09A1产品分别合成月数据。

1.2.2 气象数据

2015年的气象数据来源于中国气象局气象数据中心(http://data.cma.cn/site/index.html), 剔除数据缺失和有异常值的站点, 最终选取33个气象站点的月数据集, 包括：日照百分率、月平均气温和月降水量数据。

1.2.3 土地利用数据

100 m分辨率的土地利用数据来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)。本文参照冯益明等^[17]针对甘肃省提出的土地利用/覆盖类型的划分, 将研究区的土地利用/覆盖类型划分为9种类型：耕地、林地、草地、建筑用地、水域、沼泽地、裸沙地、戈壁和其他用地, 并将其分辨率重采样为250 m。

1.2.4 其他数据

河西走廊行政区划图来源于国家基础地理信息中心(http://ngcc.sbsm.gov.cn)。DEM数据获取自地理空间数据云(http://www.gscloud.cn), 在DEM基础上提取各栅格单元的海拔、坡度和坡向信息。

2 研究方法

CASA模型估算植被NPP的基本思想是利用植被获取太阳辐射, 加上植被自身利用的情况, 从而估算出植被净生长状况^[18]。模型中所估算的NPP可以由植被吸收的光合有效辐射(APAR)和实际光能利用率(ε)两个因子来表示^{[1, 4]}, 公式如下：

(1)

式中, x代表单个像元, t表示月份, APAR(x, t)则表示像元x在t月吸收的光合有效辐射(gC/m²), ε(x, t)表示单个像元x在t月的实际光能利用率(gC/MJ)。

2.1 APAR(x, t)的估算

APAR(x, t)是通过植被所能吸收的太阳有效辐射与植被对入射光合有效辐射的吸收比例来确定, 公式如下：

(2)

式中, SOL(x, t)表示t月在单个像元x处的太阳总辐射(gC/m²), FPAR(x, t)表示植被层对入射光合有效辐射的吸收比例, 常数0.5表示植被所能利用的太阳有效辐射(波长为0.4—0.7 μm)占太阳总辐射的比例。

由于研究区辐射站点较少, 直接利用站点太阳辐射数据插值误差较大, 故本文SOL(x, t)的估算采用和清华等^[8]将天文辐射作为起始值计算西部地区太阳总辐射的公式,

(3)

式中, Q为太阳总辐射；a、b的取值分别为0.185、0.595；s表示日照百分率；Q_A表示天文辐射, 鉴于篇幅有限, 这里不再赘述, 详细计算请参考文献[8]。

植被对太阳有效辐射的吸收比例取决于植被类型和植被覆盖状况^[19]。研究证明, 由遥感数据得到的归一化植被指数(NDVI)能很好地反映植被覆盖状况^[4]。本文利用NASA-MOD15算法中设计的NDVI-FPAR查找表(Mynenietal., 1999)计算每月截取的光合有效辐射比例(FPAR), 选取的是MODIS-NDVI数据, 其数据质量和分辨率都高于原模型中运用的NOAA/AVHRR NDVI数据^[20-21], 计算式如下：

(4)

2.2 光能利用率的估算

光能利用率是在一定时期单位面积上生产的干物质中所包含的化学潜能与同一时间投射到该面积上的光合有效辐射能之比。环境因子如气温、土壤水分状况以及大气水汽压差等会通过影响植被的光合能力而调节植被的NPP, 计算公式如下：

(5)

式中, T_ε1(x, t)和T_ε2(x, t)表示低温和高温对光能利用率的胁迫作用；W_ε(x, t)为水分胁迫系数, 反映水分条件对光能利用率的影响；ε_max是理想条件下的最大光能利用率(gC/MJ), 取值因植被类型而不同, 取值参考文献[17]。

2.2.1 温度胁迫系数估算

T_ε1(x, t)的估算：其反映在低温和高温时植被内在的生化作用对光合的限制而降低净初级生产力。

(6)

式中, T_opt(x)表示某一区域一年内NDVI值达到最高时的当月平均气温(℃), 即植被生长的最适温度。

T_ε2(x, t)的估算：表示环境温度从最适温度T_opt(x)向高温或低温变化时植被光能利用率逐渐变小的趋势, 这是因为低温和高温时高的呼吸消耗必将会降低光能利用率, 生长在偏离最适温度的条件下, 其光能利用率也一定会降低, 计算公式如下：

(7)

式中, T(x, t)表示月均温, T_opt(x)同上。当月均温T(x, t)比最适温度T_opt(x)高10℃或低13℃时, 该月的T_ε2(x, t)值等于月平均温度T(x, t)为最适温度T_opt(x)时T_ε2(x, t)值的一半。

郭婧等提出AMMRR插值方法(“多元回归分析+残差法”)能充分体现复杂多变的地形特点^[22]。因此, 本文通过经度, 纬度和海拔3个变量, 利用AMMRR方法对气温进行插值, 该插值方法不仅考虑了海拔对气象要素的影响, 而且还利用站点数据与插值数据之间的差值进行了修正, 提高了插值精度^[23-24]。

2.2.2 水分胁迫系数W_ε(x, t)的计算

水分胁迫系数W_ε(x, t)反映了植被所能利用的有效水分条件对光能利用率的影响, 随着环境中有效水分的增加, W_ε(x, t)逐渐增大, CASA模型中它的取值范围为0.5(在极端干旱条件下)到1(非常湿润条件下)。

由于W_ε(x, t)与植被水分含量直接相关, 地表水分指数(LSWI), 是描述植被叶片水分含量的一项指标^[25-27]。短波红外(SWIR)对植被水分含量敏感, 近红外和短波红外波段已经被用来获取对水分敏感的植被指数(LSWI), 计算公式如下：

(8)

式中, LSWI计算公式见(9), LSWI_max表示单个像元生长期内LSWI的最大值, 使用MVC方法估算。

(9)

这里, ρ_nir和ρ_swir分别表示MOD09A1数据的近红外和短波红外波段。

Xiao等^[25]在VPM模型中首次运用该方法, 目前已有许多学者将LSWI引入CASA模型来计算W_ε(x, t)^[28-29]。值的注意的是, 上述公式计算出来的W_ε(x, t)的取值范围为0(非常湿润)到1(极端干旱), 与CASA模型中W_ε(x, t)的取值正好相反, 为了适应模型算法, Bao等^[15]提出了修改算法, 并证明该方法适用于干旱半干旱区。因此, 本文采用地表反射率产品MOD09A1计算W_ε(x, t), 这不仅增强了数据的准确性, 也提高了数据的分辨率, 同时遥感数据也包含了地形信息^[30], 公式如下：

(10)

3 结果与讨论 3.1 精度评价

对NPP的验证可采用实测值验证法和相对验证法两种方法, 考虑到研究区范围较大且植被类型复杂, NPP实测值短时间内较难获取, 故本文采用与NPP产品比较和与前人的研究成果比较的方法。

MOD17A3产品已经被验证适用于区域和全球尺度的NPP研究^[31-32], 也有很多学者将产品与估算值进行比较来判断估算值的可靠性^[33-36]。由于2015年的MOD17A3产品未更新, 只能使用2014年产品, 本文采用同样的方法估算出2014年NPP, 两者随机抽样比较, 结果见图 2。可以看出, 校正后的NPP与MOD17A3产品的相关性较强, R²为0.803(P < 0.01), 明显高于校正前的R²(0.681)。

本文也与前人的相关研究进行了对比, 河西走廊NPP的估算鲜有研究, 因而采用干旱区其他区域的NPP均值进行比较。高志强等^[37]利用CASA模型估算的2000年西北地区年均NPP为102.52 gC m^-2 a^-1; 裘骏一^[38]利用CASA模型估算中卫市沙坡头自然保护区植被的年均NPP为60—135 gC m^-2 a^-1；焦伟等^[39]利用CASA模型估算了西北干旱区2000—2014年的植被NPP, 年均值为191.63 gC m^-2 a^-1；本文基于校正的河西走廊年均NPP为151.51 gC m^-2 a^-1, 与前人的研究结果相近。各植被类型的比较见表 1, 分析发现, 草地、耕地和荒漠的模拟结果与前人研究的结果相对较好。林地、水域和沼泽与朱文泉和焦伟的结果相差较大。

造成差异的原因较多：一是数据源的差异；二是研究区空间尺度和空间分辨率的差异, 可能存在混合像元；三是研究时段不同；四还存在算法的差异等。总的来说, 校正后的NPP与MOD17A3产品相关性更高, 与前人的计算结果更加接近, 由此判断, 校正后的NPP估算结果更加准确。

3.2 河西走廊NPP的时空分布特征

2015年河西走廊NPP空间分布见图 3, NPP自东南向西北递减, 介于0—1699.31 gC m^-2 a^-1之间, 年均值为151.51 gC m^-2 a^-1, 总量为34.29 TgC/a, 高值区主要分布在祁连山和绿洲, 低值区分布在北部荒漠, 主要植被类型NPP见表 1。

从流域来看, 石羊河流域NPP总量为13.30 TgC/a, 均值为307.19 gC m^-2 a^-1, 黑河流域总量为13.53 TgC/a、均值为243.59 gC m^-2 a^-1, 疏勒河流域总量为7.45 TgC/a, 均值为58.44 gC m^-2 a^-1。石羊河流域中上游以耕地和草地为主, NPP值高, 下游荒漠和戈壁面积大, NPP值低。黑河流域上游祁连山的东南部及南部NPP较高, 植被类型主要为森林和草地, 中游绿洲、荒漠和戈壁相间分布。疏勒河流域上游祁连山区的植被类型主要为高山草甸, NPP较高, 下游荒漠、戈壁广布, NPP较低, 年累积NPP高值区主要集中在绿洲, 范围较小。

3.3 校正结果的比较 3.3.1 校正前后平均气温比较

图 4为2015年7月校正前后的月平均气温, 由图可以看出校正后的月平均气温与实际月平均气温更吻合, 更能反映河西走廊复杂的地形特点。校正前的7月平均气温都在10℃以上, 而河西走廊南部山地的海拔大多在3000 m以上, 海拔最高达到了5641 m, 因此, 虽然7月是全年气温最高的月份, 但在海拔较高的地区气温存在低于0℃是合理的。总体上, 河西走廊月平均气温由于纬度地带性差异自南向北递减, 同时受海拔的影响也非常显著, 南部祁连山区为气温的低值区, 安敦盆地为全区域气温最高区, 与封志明等^[43]的研究结果一致。

3.3.2 校正前后水分胁迫系数比较

图 5是2015年1月和7月河西走廊校正前后水分胁迫系数空间分布。河西走廊除祁连山区外属于干旱区, 植被基本都处于水分胁迫状态。校正前1月和7月的水分胁迫系数在0.50—0.58之间, 差异很小, 全区域均处于高度水分胁迫状态, 与祁连山区植物生长的水分胁迫程度并不相符, 王海军等^[44]的研究结果表明祁连山区的降水具有明显的季节性, 冬季降水均在13 mm以下, 而在夏季降水量最高可达247 mm, 因而在夏季祁连山区的水分胁迫程度明显降低。在干旱半干旱地区植被生产力的限制因子随海拔会发生转变, 低海拔地区植被主要受到低降水导致的干旱胁迫, 同时较高的温度导致蒸散增加, 进一步加剧了植被的水分胁迫^[45]。校正后1月和7月的水分胁迫系数分布在0.5—1之间, 更能反映河西走廊制约植被生长的水分胁迫情况, 区别出地形变化对水分胁迫的影响差异, 同时还能较好的识别绿洲等人类活动区。1月份虽然降水较少, 但植被停止生长并不需要太多水分, 7月份降水较多, 但植被在生长季需水也最多, 因此, 水分胁迫在生长季(7月)比1月份更高(除祁连山外)是合理的。值的一提的是, 水分胁迫系数校正前后的计算原理差异很大, 前者是实际蒸散量与潜在蒸散量的比值, 只与降水、太阳净辐射和温度有关, 与植被无关, 后者是根据地表反射率反演水分指数而得, 与植被有关, 即同样的水分状况, 针对不同的植被和不同的生长阶段是不一样的, 因此, 校正更能反映植被的水分胁迫程度, 这种情况在中下游绿洲区表现尤为明显。位于图西北部有部分高值, 这是由于MOD09A1影像中的云量所致。

3.3.3 校正前后NPP比较

2015年河西走廊校正前后的NPP空间分布基本一致, 也存在较为明显的差异, 见图 6。可以看出, 地形平坦的下游荒漠区吻合较好, 地形起伏较大的南部高山区, 绿洲区存在明显差异, 校正前前者低估, 后两者高估。整体来看, 差值介于-526.891—296.5 gC m^-2 a^-1之间, 在山区和绿洲相差较大, 平原地区差异较小。

校正对南部高山区的温度胁迫和水分胁迫均有影响, 由图 4、5可以看出, 校正前气温高估和水分低估, 气温高估导致蒸发量增加, 增强该地区的水分胁迫, 这个效应与校正前水分低估叠加, 加剧水分胁迫导致该区域NPP低估。绿洲区地势平坦, 校正对气温的影响较小, 基于地表反射率的水分胁迫校正能准确反应植被的需水状态, 主要植被类型是耕地, 在生长季水分需求大, 校正后水分胁迫程度更高, 因此校正前NPP存在高估。此现象可进一步拓展到人类活动区域, 即改变了原始植被类型的区域, 校正对水分胁迫系数的影响非常大, 进而对NPP估算产生影响。除绿洲区外的其他荒漠区地形起伏较小, 植被分布稀少, 吻合较好。

Sun等^[46]利用地形校正的BEPS(Boreal Ecosystem Productivity Simulator)模型模拟了武陵山区NPP, 得出校正后NPP低于校正前的地区主要分布在山区, 与本研究的高山地区校正前NPP低估相反, 其原因是武夷山区植被的主要胁迫因子是温度, 校正后温度降低使得温度胁迫加强, 校正后NPP减少。平地区校正前后NPP偏差较小, 在10 gC m^-2 a^-1以内, 高山地区NPP偏差可能会达到500 gC m^-2 a^-1, 与本文的结论一致。Chen等^[47]对秦岭西南部山区NPP研究发现地形因子对NPP的影响很大, 尤其是在山区, 也与本文的结果一致。

3.4 讨论 3.4.1 校正对月NPP的影响

从2015年河西走廊校正前后各月NPP变化曲线(图 7)可以看出, 各月NPP差异并不相同, 校正后的NPP除7月份外均高于校正前, 5月份NPP校正的影响相差最大(低估14.03 gC/m²)。春季(3、4和5月)是植被开始生长期, 对水热的变化最敏感, 地形对NPP的影响较大, 校正前NPP低估了21.27 gC/m²。夏季(6、7和8月)生物量大量累积, 也是水热需求量最大的时期, 校正前低估了6.75 gC/m²。秋季(9、10和11月)气温下降, 植被基本停止生长, 此时植被对水热的需要是维持性的, 校正对NPP影响较小；冬季(12、1和2月)为非生长季, 校正对NPP影响可忽略不计。

为了展示校正对月NPP的影响细节(图 8), 本文在肃南裕固族自治县西部截取高低估和吻合区域相间分布的一部分地区(图 1), 海拔介于1737—5332 m, 地形起伏较大, 降水为140—216 mm。校正在非生长季的1、2、3、11和12月对NPP的影响较小, 大多处于-10—10 gC/m², 而在生长季的4、5、6、7、8、9和10月, 校正对NPP的影响较大, 介于-290—161 gC/m²之间。校正前6、7、8月份出现较多高估区, 该区域海拔较高, 气温较低降水比较充沛, 主要胁迫因子为温度, 校正后温度降低, 因此校正前是高估的。

3.4.2 校正对地形因子与NPP关系的影响

地形能够导致局部的水热再分配, 尤其是干旱地区, 水热组合差异影响植物生长甚至造成植被种类不同^[48]。本文由DEM提取出海拔、坡度、坡向信息, 利用等差分级法(海拔级差为30 m^[49]、坡度级差为3°^[50]、坡向分为9级^[51]), 校正对NPP与地形关系的影响见图 9。

校正对NPP在不同海拔范围内的影响各异(图 9)。在2200 m以下荒漠所占面积大, 地势平坦, 校正对NPP的影响较小, 与上文平原区变化不大结论相符。2200—3500 m范围内, 部分绿洲区处于该高程带, 校正前NPP高于校正后, 其原因是基于LSWI的校正更能反映其水分胁迫状况。3500—5000 m海拔较高, 气温下降, 植被蒸腾作用减弱, 校正后水分更能满足植被生长需求, 校正前低估了NPP。海拔5000 m以上的面积非常小, 而且都是冰川雪地, NPP基本为0, 校正的影响可忽略不计。

校正对坡度与NPP关系的影响也较大, 主要体现在土壤含水量上, 坡度与土壤含水量呈负相关, 即随坡度的增大而降低^[52]。小于15°属于缓坡, 坡度对土壤水分持有量的影响不显著^[53], 所以校正对该区间影响不大(图 9)；大于15°区间, 校正后NPP减少, 其原因是在重力作用下土壤含水量流失的坡度效应, 并且, 坡度越大两者的差异越大, 因此, 随坡度增加校正对NPP的影响也越大。

坡向导致地表接受的太阳辐射不同, 使得不同坡向水热有显著差异^[54], 从而影响植物生长。各个坡向的NPP变化见图 9, 东坡NPP最高, 南坡NPP最小。校正后北坡, 东北坡、西北坡和西坡高于校正前, 平地、东坡、东南坡、南坡、西南坡低于校正前。北坡差异最大, 校正前低估了7.12 gC m^-2 a^-1, 其次是东南坡, 校正前高估了3.38 gC m^-2 a^-1, 其他坡向的差值都低于2 gC m^-2 a^-1。坡向校正前不能够识别水热在坡向上的变化, 校正后阳坡太阳辐射强烈气温高蒸发强烈, 水分减少, 导致NPP减少, 阴坡气温低, 蒸发较弱, NPP增加, 与Sun等^[46]在武陵山区的研究结果一致。

4 结论

(1) 温度胁迫系数基于地理因子回归校正, 水分胁迫系数基于遥感数据校正能有效提高CASA模型的估算精度。河西走廊NPP在原CASA模型中存在高估, 校正后NPP总量为34.29 TgC/a, 原CASA模型高估了0.23 TgC/a。

(2) 校正对高海拔和地形起伏较大区域, 以及人类活动地区的NPP估算影响较大。在河西走廊, 前者主要是温度胁迫的影响造成, 后者主要受水分胁迫的影响。

(3) 绿洲区校正结果表明, 在人类活动区域, 基于遥感数据的水分胁迫计算比原模型中基于蒸散量的计算更加可靠, 后者存在高估。

(4) 校正对月份, 海拔、坡度以及坡向与NPP关系的影响也较大。校正对生长季的影响大于非生长季；海拔上, 2200—3500 m范围内, 校正前高估NPP, 在3500—5000 m, 校正前低估了NPP；坡度上, 校正前高估了NPP, 坡度越大两者的差异越大；坡向上, 校正前高估了阳坡NPP, 低估了阴坡NPP。

本研究的不足之处, 一是地表反射率产品存在有云现象, 导致部分区域估算不准确; 二是太阳辐射是通过克吕格插值得到, 精度有限; 三是缺乏NPP实地采样验证, 今后都有待深入研究。