1. 引言

青藏高原是中国生态学研究的热点地区 [1] [2] ，该区域分布的草地生态系统是中国草地生态系统的主要组成部分 [3] [4] 。青藏高原大面积的高寒山区草地生态系统与中亚 [5] 、蒙古 [6] [7] 干旱型草地生态系统、美国 [8] [9] 中部草地生态系统和欧洲 [10] [11] 高寒山区草甸生态系统有着显著的差异，因此青藏高原高寒山区草地生态系统可以补充大尺度高寒草地生态系统的认知。草地是我国碳循环系统研究的重点生态系统类型 [12] [13] 。近年来全球气温变化显著，青藏高原温度升高相对迅速 [14] ，研究青藏高原地区GPP变化，有利于解决全球变暖的问题，加快我国实现碳中和与碳达峰的目标 [15] 。

生态系统生产力的研究一般采用实测、器测、遥感反演和生态系统过程模型模拟等方法进行。实测受制于人工采样效率，很难进行整个青藏高原生态系统的研究，仅能利用文献进行meta荟萃分析 [16] [17] [18] 。在青藏高原上仅有少数的涡度观测通量塔，而且空间分布不均匀 [19] [20] 。在青藏高原已有大量使用遥感方法研究青藏高原GPP的变化研究工作，如He等 [21] 采用MODIS和GLASS的GPP遥感产品分析青藏高原的GPP变化 [22] [23] 。遥感方法能够很好地捕捉到青藏高原GPP空间分布格局，但不同产品的GPP年均值和时间序列上的年代际变化均不同，很难衡量气候变化背景下青藏高原生态系统演变规律。

模型可以很好的补充遥感等研究方法的不足。Piao等 [24] 基于IPCC 2013年的报告提出采取模型的方法由于数据样本种类多具有不确定性。然而现有的研究在使用模型模拟所得到的结果往往空间分辨率较低、时间序列较短，对于GPP预测与分析的准确性不高 [25] ，难以满足更高精度研究的需求 [26] [27] 。Biome-BGC模型是一种可以反映生态系统内的各种成分循环过程，以及其与外界环境相互作用所产生变化的一种过程模型 [28] 。通过Forest-BGC模型发展而来的Biome-BGC可以较好地拟合青藏地区的碳格局 [29] [30] [31] 。可以模拟不同区域GPP的变化，直接反映GPP在生态系统中的时空特征 [32] 对研究青藏地区碳排放有着重要意义。但是分辨率不够，时间序列不长。

针对以上问题，本文采用Biome-BGC模型，利用高分辨率的气象数据和地理信息数据，在青藏高原进行1979~2018年0.1˚分辨率的模拟。通过研究青藏高原不同地区的GPP的变化，以及各个环境因子对GPP的影响，并对GPP变化趋势进行研究分析。

2. 研究数据与方法

青藏高原，被誉为“世界屋脊”，其多年平均温度在−6.5℃到12.9℃之间，多年平均降水量则介于50到2000毫米 [33] 。在本研究的模拟区域，共涉及25,022个格点，覆盖经度73.05˚到104.2˚，纬度26.25˚到39.95˚。如图1(a)所示该区域的平均海拔高度超过4000米，地形自西向东逐渐降低，最高点可达8000米以上。根据图1(b)数据显示，青藏高原大部分区域较为贫瘠，草地面积约160万平方公里 [34] ，其他区域包括裸地、灌木和常绿阔叶林等生态类型。

在本研究中，我们运用了Biome-BGC模型来模拟青藏高原的碳循环过程。Biome-BGC是一个综合利用气象数据和植被生理生态参数的生态系统模拟模型，能够处理从单个点到整个陆地生态系统的空间尺度。该模型在中国的草原和森林生态系统中已获得广泛验证 [35] [36] [37] 。模型输入包括日尺度气温、降水和向下短波辐射数据。本研究使用了中国科学院青藏高原研究所制作的CMFD数据集(China Meteorological Forcing Dataset, CMFD)作为气象驱动数据 [38] ，时间分辨率为3小时，空间分辨率为0.1˚，覆盖1979至2018年。这套数据已经在中国多个地区和青藏高原进行过了大量评估有较好的结果 [39] [40] 。此外，我们还使用了MODIS的IGBP土地利用分类，GTOPO30S的DEM数据(空间分辨率1 km，双线性插值至0.1˚)以及FAO的土壤数据(空间分辨率0.1˚)，以此来模拟青藏高原在1979~2018年间0.1˚分辨率的GPP时空变化。

为验证Biome-BGC模型模拟结果的可靠性，本研究使用了两种方法：一是直接利用ChinaFLUX通量塔数据，二是参考相关文献数据进行验证。ChinaFLUX是一个监测碳、氮等通量的网络 [41] 。本研究选取了青藏高原上的三个通量塔(海北站、海北南、海北) 2003~2008年的数据，以此来验证Biome-BGC模型的拟合程度。此外，还通过对比其他文献中的数据与本研究的计算结果，进行了拟合分析。评估模型时，我们采用决定系数R²和均方根误差RMSE作为统计指标，其中R²用于衡量数据与观测数据和文献数据的趋势线拟合程度，而RMSE则表示数据间的平均差值。除点尺度评估外，如表1所示本研究还利用遥感数据可以较好的捕捉到GPP空间分布格局的特性 [42] ，在更大范围区域评估Biome-BGC模型模拟结果，进一步证实了其可靠性。

在本研究中，还研究了青藏高原上不同生态系统之间的差异和气温、降雨、向下短波辐射强度与GPP之间的关系，分析青藏高原不同生态对环境变化的响应，从而更加精确的观察GPP的变化趋势。

3. 结果

3.1. 模拟结果验证

在表2中，模型预测的GPP值与实际观察到的GPP值之间存在一定的差异。这种差异可能源于模型的不完善，或者是观察数据的误差。在一些站点，如Ali站点，模型预测的GPP值为1.39 gC·m⁻²·yr⁻¹，而实际观察到的GPP值为15 gC·m⁻²·yr⁻¹，差值为−13.61 gC·m⁻²·yr⁻¹。这说明在这个地区，模型严重低估了生态系统的生产力。这可能是因为该地区存在特殊的地形或气候条件，导致模型无法准确预测其生产力。在阿里站，阿柔站，当雄站，大野口站，海北(南)站，海北(西)站，玛曲站，大沙龙站，改则站，温泉站站点，模型预测的GPP值高于实际观察到的值。而在纳木错站点，模型预测的GPP值低于实际观察到的值。

如图2(a)~(c)所示对海北灌丛、海北湿地、当雄草地三个站点的模拟数据结果进行验证。当雄草地、海北灌丛和海北湿地的R²分别是0.71、0.65、0.55；三者数据均p < 0.01可信度较高；海北灌丛、海北湿地、当雄草地三个站点均方根误差RMSE分别为0.90 gC∙d⁻¹、1.83 gC∙d⁻¹、2.41 gC∙d⁻¹。误差较小模拟结果较好，模型数据可靠。

本研究选取多个学者实地测量的多个地区站点(见表2)的数据与Biome-BGC模拟结果进行分析。如图2(d)所示，通过拟合结果得出R²为0.88由此可知每对数据间的线性程度较强，RMSE = 240.98 gC∙d⁻¹。且本文通过BGC模拟青藏地区多年GPP的变化与大多学者研究过结果阶级且p < 0.001表明置信度较好。因此本研究采取的研究方法较为准确、可信度高。

3.2. 青藏地区40年多年平均GPP空间格局及空间趋势变化

在图3(a)中，青藏高原昆仑山脉、柴达木盆地和南部某些区域的GPP主要在0到250 gC∙m⁻²∙yr⁻¹范围内。而青海湖、金沙江和横断山脉等东部地区的GPP较高，大多在1000到1250 gC∙m⁻²∙yr⁻¹之间，尤其是青海湖、横断山脉和四川省中部地区，GPP最高可达1250到1500 gC∙m⁻²∙yr⁻¹。在图3(b)中，1979至2018年间，青藏高原大部分地区GPP呈现显著增长趋势，中部和东北部增长为0到6 gC·m⁻²·yr⁻²，南部某些地区增长高达0到30 gC·m⁻²·yr⁻²。而西北部和东南部则呈现下降趋势，约为−6到0 gC·m⁻²·yr⁻²，与Deng等的研究结果相似 [49] ，证明模拟结果的可靠性。

3.3. 基于模型与遥感数据下的青藏地区GPP变化趋势

在图4中，过去40年中，GPP整体呈上升趋势。2000年前GPP较低，但之后显著提升。1979至2018年间，BGC、GLASS和微波数据的GPP最大值均约为350 gC·m⁻²·yr⁻¹，尽管出现年份不同，但变化趋势一致。EC-LUE的GPP观测数据较低，最大值约300 gC·m⁻²·yr⁻¹，但变化趋势相似。1993年微波数据的GPP异常低，被排除。EC-LUE数据从1982年的211.92 gC·m⁻²·yr⁻¹增至2018年的292.08 gC·m⁻²·yr⁻¹，约每年增长0.28%。GLASS数据从1982年的238.41 gC·m⁻²·yr⁻¹增至335.25 gC·m⁻²·yr⁻¹，约每年增长0.33%。VOD数据呈现复杂趋势，但总体上升。Biome-BGC数据从1979年的289.86 gC·m⁻²·yr⁻¹增至356.47 gC·m⁻²·yr⁻¹，年增长率约0.23%。EC-LUE数据在1998年后增长显著，GLASS数据增长较平稳。

3.4. 不同土地利用下GPP的变化趋势

在图5中，各类土地的空间平均初级生产力(GPP)均为正值，其中森林GPP值最高(974.45 gC∙m⁻²∙yr⁻¹)，其次是草原(543.29 gC·m⁻²·yr⁻¹)，沙漠最低(28.07 gC·m⁻²·yr⁻¹)。所有土地类型的GPP均呈现上升趋势，特别是森林的增长最快(2.99 gC·m⁻²·yr⁻²)，而沙漠增长最慢(0.61 gC·m⁻²·yr⁻²)。草地和草甸GPP增长趋势高于灌木林，可能因CMFD数据集中降雨量和温度增加。四十年间，森林和草地GPP年际变化大，显示对气象因素敏感；灌木林地和荒漠在1995~2005年间也显著增长。这些结果揭示了气候因素与植被动态之间的复杂相互作用。

3.5. 生态因子对GPP的影响

如图6所示，不同生态系统类型的初级生产力对环境变化反应存在差异性。研究发现，当温度升高时，开阔灌木丛和草地的GPP值增加，在10℃左右时可达到1200 gC·m⁻²·yr⁻¹的最高值，而裸地的GPP变化趋势较为复杂；当向下短波辐射超过200 W/m²时，开阔灌木丛和草地的GPP出现下降，裸地GPP变化不明显；随着降雨量增加，各类生态系统的GPP均出现不同程度的增长，其中草地的GPP最大可达1000 gC·m⁻²·yr⁻¹。这说明不同系统对环境变化的敏感性存在差异。另外研究还发现，温度和辐射超过一定阈值，或降雨过多时，生态系统GPP与环境因子的关系变得复杂。总体而言，当环境变化较小时，GPP与温度、辐射量和降雨量保持较好的线性关系，这可用于特征化不同生态系统的韧性变化规律，评估气候变化的影响效应。

4. 讨论和结论

近年来，对青藏高原总初级生产力(GPP)变化的研究引起了广泛关注。该领域的研究内容多样，从具体地区(如青海湖 [50] 、Arou地区 [51] )到更广泛的空间(如全中国 [52] 、全球GPP分析 [53] )。然而，现有研究并未达成广泛一致性，对青藏地区GPP变化的理解存在差异。本研究旨在填补这一空白，通过运用Biome-BGC模型，分析了青藏高原GPP的变化及其对环境因子的响应。结果表明，模型模拟的GPP通常高于观测值，但总体上具有参考价值。同时，我们发现不同类型植被对环境因子的反应存在规律性变化，但GPP的变化受多种环境因子的影响，这些因子之间的相互作用使得GPP变化难以预测。模型存在一些问题，需要对参数进行更多测试和改进。此外，数据的空间分辨率也对研究结果有所影响。因此，对青藏地区GPP变化的研究需要更全面、精确和多元化的方法，包括改进模型、考虑更多环境因子，以及探索新的分析和计算方法，以更好地理解和预测青藏地区的环境变化。

本研究通过使用Biome-BGC模型在青藏高原1979~2018年对生态系统GPP模拟，得出如下结论：

1) 本研究使用的方法和数据具有很高的可靠性，模拟的GPP与观测数据拟合度非常好，三种区域数据的决定系数R²分别为0.71、0.65、0.55，均方根误差较小，且与模拟数据的线性拟合度极高(R²= 0.93)。

2) 在青藏高原，大部分区域的GPP值在0~250 gC·m⁻²·yr⁻¹之间，而青海湖和横断山脉等地区的GPP最高，达到1250~1500 gC·m⁻²·yr⁻¹。40年间，青藏高原GPP呈增长趋势，尤其在北部地区。

3) 使用BGC模型模拟的1979~2018年间的GPP最大值约为350 gC·m⁻²·yr⁻¹，与GLASS及微波数据相似，而ECLUE数据观测的GPP值略低，最大值约为300 gC·m⁻²·yr⁻¹。

4) 各个植被类型的GPP均呈增长趋势。森林的GPP值最大，为889.5 gC·m⁻²·yr⁻¹，增长趋势为2.99 gC·m⁻²·yr⁻¹；灌木的GPP值较大，为657.6 gC·m⁻²·yr⁻¹；草地GPP增长较快，为2.06 gC·m⁻²·yr⁻²；沙漠的GPP值为39.15 gC·m⁻²·yr⁻¹，增长趋势为0.61 gC·m⁻²·yr⁻²。

基金项目

国家自然科学基金(批准号：32260294)；海南大学科研基金项目(批准号：KYQD(ZR)-22083)。

NOTES

^*通讯作者。

参考文献