浑善达克沙地位于内蒙古中部锡林郭勒草原南侧(41°46'—44°24'N,112°22'—117°57'E),地势西南高,东北低,平均海拔1 300 m^[11]。分布在锡林郭勒盟的苏尼特左旗、苏尼特右旗、阿巴嘎旗、锡林浩特市、镶黄旗、正镶白旗、正蓝旗、多伦县、二连浩特市少部分和赤峰市的克什克腾旗10个旗县(市)(图1)。气候温和,属中温带大陆性气候,年平均气温为1.5℃,1月份平均气温-18.3℃,7月份平均气温18.7℃,极端最高温度35.9℃,极端最低气温-36.6℃。全年降雨量为365.1 mm,而且主要集中在7月、8月、9月份,约占全年降雨量的80%～90%。植被类型以灌木、半灌木草原、草原化荒漠和草原地带的沙地植被为主体,其生态系统对外界干扰非常敏感、脆弱^[12]。

图1 浑善达克沙地地理位置

气象数据来自中国气象科学数据共享服务网(http:∥cdc.cma.gov.cn)提供的2000—2016年沙地境内10个旗县气象站点的月平均温度数据。太阳辐射数据来自于内蒙古气象局提供的覆盖全内蒙古的8个太阳总辐射站2000—2016年的观测数据。对气象数据采用克里金插值法进行空间插值,获取与NDVI数据投影、像元大小相同的250 m分辨率的栅格数据。太阳辐射数据通过反距离插值法对全内蒙古的8个太阳辐射站的月值数据插值为1 000 m的栅格数据,然后进行数据重采样、掩膜获取与NDVI数据分辨率相同的栅格数据。

遥感数据来自NASA对地观测系统数据(http:∥modis.gsfc.nasa.gov/)共享平台。本文使用涵盖浑善达克沙地16 d合成的MOD13Q1数据产品以及8 d合成的MOD09A1的数据产品。MOD13Q1产品的空间分辨率为250 m,因内蒙古大部分地区的植物在冬天不再生长或者被冰雪覆盖等原因,选取生长季(4—10月)分析NPP的变化趋势,运用最大值合成法(MVC),生成月NDVI数据,降低大气等外界因素的干扰^[12]。选用分辨率为500 m的MOD09A1产品反演地表湿润指数(LSWI),LSWI是近红外波段(ρ_nir)和短波红外波段(ρ_swir)的归一化指数^[13],并根据地表湿润指数来估算水分胁迫因子。将获取的月NDVI数据和插值得到的月气象数据、月太阳辐射数据通过CASA模型求出月NPP数据,并将每年4—10月的7个月数据进行累加获取2000—2016年NPP的时空数据集。

以2016年8月为例,依据CASA模型建模思路,通过对SOL,FPAR,T_ε,W_ε众因子的估算,得出2016年8月的浑善达克沙地NPP取值范围为0.057～184.098 g/m²,平均值为73.434 g/m²(图2)。

图2 浑善达克沙地2016年8月NPP估算结果

由于本文缺乏实测数据,故采用间接验证的方式对估算结果进行验证。将本文CASA模型估算的NPP值与其他人估算的NPP的均值或部分数据进行对比验证^{[12,15,18-22]},来检测估算数据的准确性与可行性(表1)。

表1 浑善达克沙地NPP模拟值与其他模拟结果比较

参考文献中无月值实测数据。

本文估算的2000—2016年浑善达克沙地植被NPP多年均值为282.42 g/(m²·a),与表1对比分析,估算结果介于内蒙古地区植被NPP的模拟结果之间。且7月、8月份NPP均值为78.31,67.31 g/m²,与李刚等^[18]、穆少杰等^[20]、元志辉^[12]7月、8月计算的实测数据相符合。此外,由于生态治理工程的推进,2010年后研究区内植被NPP值不断增大,本研究结果高于元志辉^[12]估算的浑善达克沙地NPP值,该估算结果在一个合理的区间内,因此可以用于后续的结果与分析。

利用CASA模型估算的浑善达克沙地17 a生长季NPP的变化趋势,其结果表明沙地生长季NPP呈波动上升趋势,年均增长率为2.136 2 g/(m²·a),从2000年的276.82 g/(m²·a),增长到2016年的292.7 g/(m²·a)。在2000—2010年内,NPP值波动变化明显,2003年、2004年、2006年出现高值,2002年、2005年、2007年出现低值。而在2010年以后,每年NPP值不断增大,在2013年达到峰值,NPP值为339.96 g/(m²·a)(图3A)。

图3 2000-2016年生长季NPP年际变化和偏差分析

偏差分析法直观展示每年生长季NPP平均值与17年间NPP多年均值的偏离情况^[22],结果表明:2000—2016年,生长季NPP的偏差值呈上升趋势。2010年以前,除2003年、2004年、2006年以外,其余年份的生长季NPP值均低于多年NPP均值。2010年后,每年的生长季NPP值均大于多年均值,2014年NPP值低于2013年,但总体处于上升状态(图3B),说明浑善达克沙地植被状况得到改善和好转,与杭玉玲等^[23]关于锡林郭勒草原植被覆盖状况转好的研究结果趋于一致。

从季节尺度分析,浑善达克沙地2000—2016年春季NPP值呈下降趋势,夏季和秋季均呈上升趋势(图4)。春季气温升高,植物开始返青,NPP值逐渐增大,浑善达克沙地春季NPP值波动较大,2003年、2009年、2013年NPP值达到高值期,其中2003年达到峰值,值为57.5 g/(m²·a),研究时间内多数年份NPP值低于平均水平,NPP值以0.158 8 g/(m²·a)的速率缓慢下降(图4A),比较符合近年来内蒙古植被返青期推迟的趋势^[24]。夏季,气温、降水、太阳辐射等气候因子都为植物的生长提供良好的环境,NPP值达到最高,以2.364 6 g/(m²·a)的速率增长,2010年以后NPP上升趋势呈指数型增长,到2013年达到峰值,值为248.83 g/(m²·a),2014年以后又呈现下降趋势,但总体水平还是呈增长状态,这也与全年NPP变化趋势相一致(图4B)。秋季,气温降低降水减少,植被开始枯黄凋零,NPP值减小。2000—2016年,秋季NPP值呈现波动上升趋势,每年以0.150 3 g/m²的速率缓慢增长,秋季NPP值在2000年、2004年、2012年、2013年出现明显的高值期(图4C)。

从月变化来看,浑善达克沙地17 a内月平均NPP值变化明显(图4D)。由于不同月份的水热组合不同,各月的NPP值也各不相同,4—7月NPP值呈上升趋势,到7月达到峰值,值为78.31 g/m²,8—10月份,NPP值呈下降趋势,10月份,NPP值降到生长季最低值,夏季6—8月NPP值整体处于高值。其中,7月、8月份积累的NPP值占总量的51.57%,4月和10月的NPP值仅占总量的4.09%,3.05%。

图4 NPP季节与月变化特征

将2000—2016年的NPP值取平均值来反映浑善达克沙地17 a来NPP的整体空间分布特征,浑善达克沙地NPP值整体呈由西向东逐渐递增的趋势(附图2),东部的克什克腾旗、多伦县和正蓝旗处于高值区,且大部分地区NPP值高于350 g/(m²·a)。西部的二连浩特市、苏尼特右旗、苏尼特左旗和镶黄旗处于低值区,占地面积为研究区的46.4%,NPP值均低于250 g/(m²·a)。从不同等级NPP的空间分布来看,NPP值为0～50 g/(m²·a)的面积占研究区的0.28%,主要分布在苏尼特右旗、苏尼特左旗、阿巴嘎旗、正镶白旗和正蓝旗的湖泊水域部分。NPP值大于650 g/(m²·a)的面积占沙地面积的0.25%,主要分布在克什克腾旗的东部以及东南部。NPP值介于250～350 g/(m²·a)的区域占总面积的26.06%,属于研究区的均值区,分布在阿巴嘎旗的中西部地区、正镶白旗的南部区域、正蓝旗的西部地区和锡林浩特市的南部区域。NPP值介于350～450 g/(m²·a)的面积占沙地总面积的17.65%,集中分布在正蓝旗、锡林浩特市北部、克什克腾旗和多伦县。NPP值为450～650 g/(m²·a)的高值区占总面积的9.35%,集中分布在沙地的克什克腾旗和多伦县。

为进一步分析浑善达克沙地植被NPP的空间变化规律,本文采用一元线性回归法,基于像元尺度分析2000—2016年生长季NPP的空间变化趋势(附图3A),结果表明:沙地内78.2%的面积呈现增长趋势,21.8%的面积呈现降低趋势。NPP上升最明显的区域主要分布在多伦县、正蓝旗东南部、阿巴嘎旗北部以及锡林浩特市,变化速率超过15 g/(m²·a),研究区大多数区域NPP值的增长趋势为0～5 g/(m²·a)。苏尼特右旗东部、苏尼特左旗中部、锡林浩特市以及克什克腾旗中部区域NPP值呈现下降趋势,其中克什克腾旗中部地区植被NPP下降十分明显,值介于-10～-5 g/(m²·a)。

通过对皮尔逊相关系数进行显著性检验来分析浑善达克沙地植被NPP的变化趋势,r_xt>0表示NPP呈上升趋势,r_xt<0表示NPP呈下降趋势,若r_xt通过0.05的显著性水平检验(p<0.05),则表示NPP呈显著增加或显著减少,若r_xt通过0.01的显著性水平检验(p<0.01),则表示NPP呈增长极显著或降低极显著。未通过检验的为不显著^[1],基于此6个等级进一步分析NPP的变化情况(附图3B),研究区内12.17%的区域增长显著或增长极显著,主要分布在多伦县、正镶白旗北部、正蓝旗东部以及克什克腾旗东北区域,表明这些区域植被生长环境得到改善。研究区内2.01%的区域降低显著或极显著,主要零星分布在克什克腾旗的中部以及苏尼特左旗的东部区域。浑善达克沙地56.07%的面积属于增长不显著,19.79%的面积属于下降不显著,主要分布在沙地的中西部地区。

本文研究结果中浑善达克沙地植被整体呈现上升趋势,该研究结果与国内许多学者相一致,高晓霞等^[25]基于logistic回归模型分析浑善达克沙地的动态变化趋势,整体呈现上升趋势,孙建芸^[26]基于遥感数据分析浑善达克沙地的荒漠化情况,表明近年来浑善达克沙地植被覆盖度明显提高,荒漠化状况得到改善,元志辉^[12]分析2000—2013年浑善达克沙地植被净初级生产力也表明浑善达克沙地NPP呈上升趋势。从NPP的年内变化来看,夏季NPP值占全年比例最大,且7月、8月份积累的NPP值最高,表明NPP的变化主要与植被的生长旺季有关。研究区植被东部NPP值大于中西部,原因可能是沙地呈东西走向,具有一定的地带性,东部降水多,植被类型丰富,另外最初2000年实施生态工程后,多伦县、正蓝旗和锡林浩特市的白音锡勒牧场是最早的防沙治沙试验区^[11],因而植被NPP值较高,中西部地区降水稀少,且处于沙地腹部,生态工程治理难度大,NPP值较低。

此外,在CASA模型选取的参数因子中,光能利用率ε对NPP的估算结果影响较大^[13]。然而准确建立水分胁迫因子对ε的影响函数是艰难的,传统的W_ε计算都是基于蒸散量和土壤系数的函数,其局限性在于空间降水和土壤质地数据的空间异质性表达^[16]。本文考虑到利用MOD09A1数据提取近红外波段和短红外波段计算地表湿润指数可以减少较为复杂且涉及到较多气象数据资料的过程^[27]。此外,基于MOD09A1计算的地表湿润指数能够很好地表现空间异质性,比较准确地反映地表的湿润状况,有效减少了蒸散量和土壤系数估算的人为误差^[13]。

最大光能利用率ε_max的取值对最终的NPP估算结果也有较大影响^[17],目前,由于研究区的植被类型不同和气候环境有所差异,以往研究在区域尺度和全球尺度上求取的ε_max值各不相同,主要波动在0.389～1.044,如王莺等^[16]将甘南地区草地的ε_max取值为0.608 g C/MJ; 卫亚星等^[28]计算青海省ε_max为0.403～0.908 g C/MJ; 朱文泉等^[17]根据误差最小的原则,利用中国690个观测点NPP实测数据,模拟各植被类型的ε_max,该研究结果得到许多学者的应用^[8,19]。本文对ε_max的取值也参照这一研究成果,取值为0.542 g C/MJ。但是由于ε_max的取值随着植被类型的不同而变化^[15],沙地实施生态治理工程后,植被类型和群落组成均发生了重大的变化^[11],在后续的研究中有必要实地考察采样,获取最新的沙地植被类型分布,结合地面实测NPP数据,计算研究区主要植被建群种类型的最大光能利用率,提高NPP模拟的精度。

另外,本文方法中使用的气象数据来自气象站的离散点数据,该方法的优点在于其简单且易于实施,缺点是随着距观测点的距离增加精度降低。而且插值方法本身可能会影响结果。因此,在后续的研究中可以考虑结合遥感数据反演光合有效辐射值^[29],提高研究结果的精确性。

(1)浑善达克沙地2000—2016年NPP呈现波动上升趋势,NPP多年均值为282.42 g/(m²·a),年均增长量为2.136 2 g/(m²·a),从2000年的276.82 g/(m²·a)增长到2016年的292.7 g/(m²·a),表明近17 a浑善达克沙地植被固碳能力增加,植被长势趋于良好。

(2)NPP季节变化明显,春季NPP值呈下降趋势,夏季和秋季呈上升趋势,夏季NPP变化趋势与全年NPP变化趋势基本一致。从月变化来看,NPP值在7月达到最高值,7月、8月份积累的NPP值占总量的51.57%。

(3)空间上,NPP整体呈由西向东递增趋势,沙地内78.2%的面积呈增长趋势,21.8%的面积呈现降低趋势,2.17%的区域增长显著或增长极显著,2.01%的区域降低显著或极显著,表明国家在对浑善达克沙地进行生态治理工程后取得了一定的效果,但是浑善达克沙地局部地区植被恶化状况仍然存在,下一步有必要加强对沙地中西部地区生态工程建设的力度。