CN116341724B

CN116341724B - 基于全球气候模式驱动碳循环机理模型的碳吸收预估方法

Info

Publication number: CN116341724B
Application number: CN202310217022.2A
Authority: CN
Inventors: 于敏; 贾小龙; 程春香; 薄宇; 阮多; 徐作敏; 于大江; 张作君
Original assignee: Heilongjiang Provincial Ecological Meteorological Center
Current assignee: Heilongjiang Provincial Ecological Meteorological Center
Priority date: 2023-03-08
Filing date: 2023-03-08
Publication date: 2023-09-15
Anticipated expiration: 2043-03-08
Also published as: CN116341724A

Abstract

基于全球气候模式驱动碳循环机理模型的碳吸收预估方法，为了解决不同排放情景下未来陆地生态系统高分辨率的碳吸收数据难以获取的问题。本发明首先收集目的地历史时段数据以及CMIP6模拟数据，分为历史时段模拟数据和未来时段模拟数据，并对CMIP6中的多种模式数据进行尺度统一处理；基于目的地历史时段数据，对CMIP6中多种模式历史时段模拟效果进行综合评估，筛选每个数据种类各自对应的最优的模式并进行等权重集合；基于同一时间阶段的目的地历史时段数据与CMIP6历史时段集合模拟数据，构建模拟数据订正模型，然后基于未来集合模拟值订正后的数据，利用BEPS模型预估未来碳吸收指标。

Description

基于全球气候模式驱动碳循环机理模型的碳吸收预估方法

技术领域

本发明属于碳循环模拟技术领域，涉及一种碳吸收预估方法。

背景技术

为了更好地实现双碳战略目标，维持全球碳平衡，减缓气候变化的多方面影响，在实行生态系统保护修复和节能减排政策的同时，更应针对碳循环分布格局、增汇潜力等方面展开更加全面系统的研究。科学地评估生态系统碳收支时空分布特征和生态系统服务功能，是其潜力评估，以及适应和减缓气候变化的关键。目前也有很多专家对碳吸收估计和预估展开了研究，但是利用CMIP6模型数据驱动BEPS碳循环模型，对于不同陆地生态系统碳循环规律进行研究的案例较少。基于CMIP6全球气候模式不同碳排放路径情景下，开展碳吸收定量预估，可以为生态系统保护和修复、服务双碳目标生态气象服务提供技术支撑。

发明内容

本发明为了解决不同排放情景下，未来陆地生态系统高分辨率的碳吸收数据难以获取的问题。

基于全球气候模式驱动碳循环机理模型的碳吸收预估方法，包括以下步骤：

S1、利用CMIP6模型预测未来驱动数据变化特征：

步骤1.1、收集目的地历史时段数据：

目的地历史时段数据中的数据种类包括降水量、太阳辐射量、最高温度、最低温度、相对湿度、二氧化碳浓度、叶面积指数，以及氮沉降、地表覆盖、30年平均温度、土壤质地、聚集度指数、经纬度数据；

步骤1.2、获取CMIP6模拟数据，CMIP6模拟数据包括多种模式数据，CMIP6模拟数据中分为历史时段模拟数据和未来时段模拟数据；

CMIP6模拟数据中的数据种类包括降水量、太阳辐射量、最高温度、最低温度、相对湿度、二氧化碳浓度、叶面积指数；将降水量、太阳辐射量、最高温度、最低温度、相对湿度记为气候要素数据；

步骤1.3、对CMIP6中的多种模式数据进行尺度统一处理；

步骤1.4、基于目的地历史时段数据，对CMIP6中多种模式历史时段模拟效果进行综合评估，筛选每个数据种类各自对应的最优的3个模式；

步骤1.5、针对每个数据种类的数据，分别将各自对应的3个模式的模拟值进行等权重集合：即(模式1+模式2+模式3)/3，得到每个数据种类对应模拟值；

将CMIP6中历史时段模拟数据和未来时段模拟数据对应的等权重集合后的每个数据种类模拟值分别记为CMIP6历史时段集合模拟数据和CMIP6未来时段集合模拟数据；

步骤1.6、基于同一时间阶段的目的地历史时段数据与CMIP6历史时段集合模拟数据，构建模拟数据订正模型：

R_future(i，j，k)＝R_CMIP6future(i，j，k)-error(i，j)

其中，是CMIP6历史时段集合模拟数据平均值，其中降水量、太阳辐射量、最高温度、最低温度、相对湿度、叶面积指数为逐日平均格点值，CO₂浓度为逐月平均值；/>是目的地历史时段数据的平均值，降水量、太阳辐射量、最高温度、最低温度、相对湿度、叶面积指数为逐日平均格点值，CO₂浓度为逐月平均值；error(i，j)是二者的偏差；R_future(i，j，k)是未来集合模拟值订正后的结果，R_CMIP6future(i，j，k)是CMIP6未来时段集合模拟数据如果对应降水量、太阳辐射量、最高温度、最低温度、相对湿度、叶面积指数，R_future(i，j，k)是各自对应的逐日平均格点值，如果对应CO₂浓度，R_future(i，j，k)是对应的逐月值；

S2、基于未来集合模拟值订正后的降水量、太阳辐射量、最高温度、最低温度、相对湿度、二氧化碳浓度、叶面积指数数据，利用BEPS模型预估未来碳吸收指标，所述碳吸收指标包括总初级生产力GPP，和/或净初级生产力NPP，和/或净生态系统生产力NEP，和/或土壤水分含量SW，和/或蒸散发ET，和/或释氧量，和/或碳储量。

进一步地，S2中基于未来集合模拟值订正后的降水量、太阳辐射量、最高温度、最低温度、相对湿度、二氧化碳浓度、叶面积指数数据，利用BEPS模型预估未来碳吸收指标之前，需要对BEPS模型进行模型预热，即：

输入目的地的历史数据，进行历史碳吸收指标模拟；

所述目的地的历史数据，为时间尺度范围更大的目的地的历史时段数据，即目的地的历史数据的数据类型与目的地的历史时段数据的数据类型完全相同，目的地的历史数据的时间尺度大于目的地的历史时段的时间尺度。

进一步地，利用BEPS模型预估未来碳吸收指标之后，对碳吸收指标进行时空分布特征评估的具体包括以下步骤：

将碳吸收指标记为Θ，Θ为GPP、NPP、NEP、SW、ET、释氧量或碳储量；

按照每k’年作为一个年代纪宽度，将预估时间段范围，以及预估时间段范围内第一年之前的一段时间范围，进行年代纪划分，共计得到k个年代纪；

预估时间段范围内第一年之前的一段时间范围应小于一个年代纪宽度，即将预估时间段范围内第一年所在年代纪记为初始年代纪；

在一个年代纪宽度内，将Θ栅格数据中的某一像元点的位置(i，j)对应的值求平均值，记为Θ_k(i，j)；初始化Θ_k(i，j)对应的年代纪特征Z(i，j)，并计算每个年代纪中Θ_k(i，j)对应的年代纪特征Z(i，j)：

最后得到所有年代纪对应的年代纪特征Z(i，j)，Z(i，j)即为空间分布特征；进而根据空间分布特征分析碳吸收指标的年代际变化特征。

进一步地，利用BEPS模型预估未来碳吸收指标之后，还要对模型模拟结果的准确度进行评估，具体包括以下步骤：

首先采用一元线性回归法，进行预估碳吸收指标变化率计算和显著性分析；然后利用Mann-Kendall趋势检验法进行突变分析实现显著性的检验，评估模型模拟结果的准确度。

进一步地，步骤1.4中采用泰勒图对CMIP6中多种模式历史时段模拟效果进行综合评估。

进一步地，采用泰勒图对CMIP6中多种模式历史时段模拟效果进行综合评估的过程包括以下步骤：

CMIP6中的数据存在多个模式，设气候要素数据的模式为N1个；针对气候要素数，将降水量、太阳辐射量、最高温度、最低温度、相对湿度分别进行如下处理：

将目的地历史时段数据中的气候要素数的COR和/或RMSE和/或STD记为基准数据BUOY；利用泰勒图，将BUOY以及CMIP6中同种气候要素N1个模式数据的COR和/或RMSE和/或STD绘制在一张极坐标图上进行比较，选择与BUOY最接近的3个模式；

其中，COR为空间相关系数，RMSE为均方根误差，STD为标准差；

针对二氧化碳浓度和叶面积指数采用相同的方式得到各自对应的3个模式。

进一步地，针对目的地的历史时段数据中的叶面积指数采用源于AVHRR和MODIS的叶面积指数；源于AVHRR和MODIS的叶面积指数需要经过异源订正处理，方法如下：

R_NEWLAIMODIS(i，j，k)＝R_LAIMODIS(i，j，k)+error(i，j)

其中是第一历史时段的AVHRR中的逐格点叶面积指数数据的平均值，/>是第一历史时段的MODIS中逐格点叶面积指数数据的平均值，error(i，j)是二者的偏差；R_NEWLAIMODIS(i，j，k)是订正后的MODIS叶面积指数的结果；R_LAIMODIS(i，j，k)为订正前的MODIS叶面积指数的结果；

所述第一历史时段为历史时段期间内部分历史时段。

进一步地，步骤1.3对CMIP6中的多种模式数据进行尺度统一处理时，气候要素数据的尺度：空间分辨率为1km×1km，时间分辨率为日；叶面积指数的尺度：空间分辨率为1km×1km，时间分辨率为8日；二氧化碳浓度的尺度：时间分辨率为月。

进一步地，步骤1.3对CMIP6中的多种模式数据进行尺度统一处理的过程中，采取最近邻插值法对空间尺度进行统一，即最近邻插值法统一空间分辨率。

本发明利用CMIP6模型数据驱动BEPS碳循环模型，实现了未来碳排放情景下(SSP126、SSP245、SSP585)GPP、NPP、NEP等的时空变化特征的预估分析。具体具有以下特点：

1.利用CMIP6模型预测未来驱动数据变化

目前，针对一定区域利用气候模型，对未来气候变化进行预测，实现成体系的筛选、集合、订正、评估的方法较少。具体的选取方法和订正手段有待于进一步探讨，尤其针对CMIP6气候模型数据进行区域化应用的范例较少。本发明改进了预估数据的筛选和订正方法，增加了定量化的指标来评估模拟数据的精度，弥补了以往气候模型模拟偏差较大的问题。此外，本发明还开创性的选取CMIP6模型中的二氧化碳浓度数据和叶面积指数数据进行预测分析，能够有效地提高模拟碳吸收相关指标的准确性。

2.利用BEPS模型模拟历史碳吸收变化，能够有效提高对未来碳吸收指标预估的准确度。BEPS模型模拟需假设模拟起始时间处于碳平衡状态，因此碳吸收预估也需要从早期历史时间点(人为影响有限，可以假设该时间点处于碳平衡状态)开始。

3.利用BEPS模型未来预估碳吸收变化特征。

由于目前还没有基于碳循环机理模型，结合未来气候数据，预估区域演变规律分析的先例。本发明将气候模型与碳循环模型相结合，选取适用于输入未来气候情景数据的模型，发明了一种预估碳吸收的方法，改进了碳吸收模拟的准确性，从无到有弥补了以往预估结果缺失，难以获取的问题。

附图说明

图1为降水的泰勒图。

图2为太阳辐射的泰勒图。

图3为平均温度的泰勒图。

图4为最高温度的泰勒图。

图5为最低温度的泰勒图。

图6为相对湿度的泰勒图。

图7为叶面积指数产品ACCESS-ESM1-5的效果。

图8为叶面积指数产品CMCC-ESM2的效果。

图9为叶面积指数产品MPI-ESM1-2-HR的效果。

图10为叶面积指数产品MPI-ESM1-2-LR的效果。

图11为本发明叶面积指数的效果。

图12为ACCESS-ESM1-5、CMCC-ESM2、MPI-ESM1-2-HR、MPI-ESM1-2-LR和本发明的叶面积指数曲线图。

图13为本发明CO₂模拟结果。

图14为黑龙江历史段订正前太阳辐射量的效果图。

图15为黑龙江历史段订正后太阳辐射量的效果图。

图16为黑龙江预估段太阳辐射量的效果图。

图17为黑龙江2015年太阳辐射量对比曲线图。

图18为黑龙江2001-2014年评估段及2015-2060年预估段的逐年太阳辐射量对比曲线。

图19为黑龙江2015年对应的订正前的降水量的效果图。

图20为黑龙江2015年对应的订正后的降水量的效果图。

图21为黑龙江2015年对应的预估段的降水量的效果图。

图22为黑龙江2015年降水量对比曲线图。

图23为黑龙江2001-2014年评估段及2015-2060年预估段的逐年降水量对比曲线。

图24为2020年NPP效果图。

图25为2020年NEP效果图。

图26为1961-2019年NPP曲线图。

图27为1961-2019年NEP曲线图。

图28为GPP的估计曲线图。

图29为NPP的估计曲线图。

图30为NEP的估计曲线图。

具体实施方式

本实施方式为一种基于全球气候模式驱动碳循环机理模型的碳吸收预估方法，经过研究发现降水量、太阳辐射量、最高温度、最低温度、相对湿度、二氧化碳浓度和叶面积指数对BEPS模型影响强于其他因素，因此本发明利用这些因素进行未来数据预估。同时想要得到未来的GPP(总初级生产力)、NPP(净初级生产力)、NEP(净生态系统生产力)，需要向碳汇模型里面输入未来的降水量、太阳辐射量、最高温度、最低温度、相对湿度、二氧化碳浓度、叶面积指数。未来的降水量、太阳辐射量、最高温度、最低温度、相对湿度、二氧化碳浓度、叶面积指数通过CMIP6来获得，CMIP6的数据可以从网站上直接下载，但是网站上下载的模拟数据与预估目的地的数据是有差异的，因此本发明需要结合目的地的历史时段数据对CMIP6模拟数据进行调整和修正，所以本发明首先需要结合目的地的历史时段数据对在网上下载的CMIP6模拟数据进行调整，即以目的地的历史时段数据为准，将下载的CMIP6模拟数据处理为近似于目的地的历史时段数据。这样就可以得到一套比较准确的、适合BEPS模型的数据集；最后再利用BEPS模型实现碳吸收预估。

本实施方式所述的一种基于全球气候模式驱动碳循环机理模型的碳吸收预估方法，包括以下步骤：

步骤一、利用CMIP6模型预测未来驱动数据变化特征：

步骤1.1、收集历史时段驱动碳循环机理模型所需驱动数据：

碳循环机理模型需要的驱动数据主要有降水量、太阳辐射量、最高温度、最低温度、相对湿度、二氧化碳浓度、叶面积指数，以及氮沉降、地表覆盖、30年平均温度、土壤质地、聚集度指数、经纬度数据。

碳循环机理模型采用BEPS中的BEPS_NJU模型，BEPS_NJU模型的输入数据包括时空变化的数据、空间变化的数据和时间变化的数据三大类。

(1)时空变化的输入数据

①气候数据

逐日气候数据：包括逐日最高气温、最低气温、平均相对湿度、日降水量和日总辐射5个变量。数据格式为短整型二进制文件(每个格点的数据占用2个字节)存储。

变量单位：日总辐射为0.1MJ/m²；相对湿度为10％(扩大10倍)；日降水量为0.1mm；最高气温为0.1°k；最低气温为0.1°k。

②叶面积指数数据

生态系统的碳蓄积是一个长期过程，为了更好地模拟现在的碳通量，需要进行陆地生态系统碳循环过程的长期模拟。驱动BEPS_NJU模型长时间序列的叶面积指数(LAI)作为输入，包括利用AVHRR数据生成的1981-1999年的LAI数据、2000年之后利用MODIS数据反演生成的LAI数据，以及利用近年的FY数据反演的LAI数据。这些数据都组织为短整型(每个格点的数据占用2个字节)的二进制文件。

需要特别注意的是：由于不同LAI产品存在差异，所以在连续使用不同的LAI数据产品时，应该进行校正处理，实现不同的LAI数据产品的“无缝”对接，以防止模拟的碳通量存在虚假“跃变”。

③氮沉降数据

BEPS模型需年氮沉降数据作为输入。数据为浮点型二进制文件(每个格点的数据占用4个字节)。数据单位：g N m-2yr-1。

(2)空间变化的输入数据

①地表覆盖数据

BEPS_NJU模型使用的地表覆盖数据采用国际地圈生物圈计划(IGBP)定义的产品。数据为字符型二进制文件。

BEPS模型将植被类型分为：水(0)、常绿针叶林(1)、常绿阔叶林(2)、落叶针叶林(3)、落叶阔叶林(4)、混交林(5)、郁闭灌丛(6)、开放灌丛(7)、多树的草原(8)、稀树草原(9)、草原(10)、永久湿地(11)、作物(12)、城市与建成区(13)、作物和自然植被的镶嵌体(14)、雪冰(15)、裸地或低植被覆盖地(16)、未分类区(254)、填充值(255)。

②土壤质地数据

BEPS_NJU模型根据土壤质地数据(砂粒、粉粒和粘粒比例)生成饱和持水量、田间持水量、凋萎系数和土壤水分饱和传导率等水分参数，用于模拟土壤水分动态模拟。假设第2和3层的土壤质地相同，第3层土壤的水文参数用第2层的土壤质地数据估算。数据为字符型二进制文件。

③聚集度指数数据

BEPS_NJU模型根据聚集度指数和冠层LAI估算阳叶和阴叶的LAI，模拟冠层辐射传输过程。BEPS_NJU模型使用由MODIS BRDF数据产品反演生成的聚集度数数据。

注意：叶面积指数反演和运行BEPS_NJU模型应该使用相同的聚集度数据。

④30年平均温度数据

模型利用30年平均温度表示每个格点的热量条件，修正自养呼吸对温度的敏感系数。数据存储为浮点型二进制文件。

⑤经纬度数据

每个格点的经度和纬度数据，数据存储为浮点型二进制文件。

⑥时间变化的输入数据

现有版本的BEPS_NJU模型仅考虑了CO₂浓度随时间变化，未考虑其空间差异，也就是全球所有地区每一年使用一个CO₂浓度数据。数据存储为文本文件。

本实施方式中构建BEPS_NJU模型的输入数据主要分为三大类：

a、时空变化的气候(逐日)、植被(逐八日)、氮沉降数据(逐年)；

b、空间变化的地表覆盖、土壤质地、聚集度指数、经纬度、30年气候态数据；

c、时间变化的二氧化碳浓度数据(逐月)。

BEPS模型历史驱动数据要求【降水量、太阳辐射量、最高温度、最低温度、相对湿度】(1901-2020分辨率1km、日)、二氧化碳浓度(1871-2020年、全省平均、月)、叶面积指数(1981-2020年、1km、2000年前逐16日，2000年后逐8日)。

针对目的地的历史时段数据中的降水量、太阳辐射量、最高温度、最低温度、相对湿度，本实施方式中的数据来源于东安格利亚大学(UEA)气候研究小组(CRU)和日本气象厅(JMA)的日本再分析数据(JRA)。

基于东安格利亚大学(UEA)气候研究小组(CRU)和日本气象厅(JMA)的日本再分析数据(JRA)，简称CRUJRAV 2.1。CRU JRAV2.1数据集是由东安格利亚大学(UEA)气候研究小组(CRU)制作的10个气象变量(2米温度、2米最高和最低温度、总降水量、比湿度、向下的太阳辐射通量、向下的长波辐射通量、压力以及风速的纬向和经向分量)的6小时陆地表面网格时间序列，旨在用于驱动模型。变量是在0.5°x 0.5°网格上提供的，该网格接近全球，但不包括南极洲。来源网址https://catalogue.ceda.ac.uk/uuid/10d2c73e5a7d46f4ada08 b0a26302ef7。

针对目的地的历史时段数据中的叶面积指数，本实施方式中的数据来源于1981-1999年之间来源于AVHRR，以及2000-2020年来源于MODIS。这是因为MODIS的叶面积指数比较准，但是2000年前没有数，除了MODIS之外，AVHRR也得到广泛应用，为了保证数据空间的连续性和准确性，因此发明需要在消除异源数据差异性的基础上使用两套数据。

由于2000年前后数据不同源，因此有必要在对于输入叶面积指数数据进行处理，处理方法如下：

R_NEWLAIMODIS(i，j，k)＝R_LAIMODIS(i，j，k)+error(i，j)

其中，是2001-2020年AVHRR中的逐格点叶面积指数数据的平均值，是2001-2020年MODIS中逐格点叶面积指数数据的平均值，error(i，j)是二者的偏差；R_NEWLAIMODIS(i，j，k)是订正后的MODIS叶面积指数的结果；R_LAIMODIS(i，j，k)为订正前的MODIS叶面积指数的结果。

AVHRR数据已被广泛用于推导陆面植被参数，并应用在各种陆面模式中，AVHRR的植被参数(即LAI、LCC和FVC)都是从NDVI数据导出的，其具有有限的可用信息和无用的冗余变量。AVHRR的LAI数据来自美国伊利诺州立大学水文研究中心，该数据时间分辨率是16天，空间分辨率是1km，覆盖全球的AVHRR数据。

MODIS：1999年12月18日，美国国家航空航天局(NASA)成功地发射了地球观测系统(E0S)的第一颗先进的极地轨道环境遥感卫星Terra，2002年4月18日NASA又发射了第二颗Aqua(PM 1)。Terra是过境时间为地方时10:30左右，Aqua过境时间为地方时14:30左右。中分辨率成像光谱仪——MODIS(Moderate-resolution imaging Spectroradiometer)是Terra和Aqua上搭载的最主要和最有特色的传感器。Terra/MODIS和Aqua/MODIS的成功发射和投入运行把遥感技术应用在全球变化研究中推进到一个新的时代。它的主要目标是实现从单系列极轨空间平台上对太阳辐射、大气、海洋和陆地进行综合观测，获取有关海洋、陆地、冰雪圈和太阳动力系统等观测信息，用于土地利用和土地覆盖研究、气候季节和年际变化研究、自然灾害监测的分析研究、长期气候变率的变化研究以及大气臭氧变化研究，实现对大气和地球环境进行长期观测和进一步的理解。为了充分了解地球系统的变化，保证数据源的连续性，美国今后的业务卫星观测系统将和EOS观测系统接轨，以提供系统的、连续的地球观测信息。MODIS数据有3个显著特点：(1)空间分辨率和光谱分辨率大幅提高。从0.4～14.0μm，MODIS有36个波段，其中2个波段(可见光0.62～0.67μm和近红外0.841～0.876μm)的空间分辨率为250m；5个可见光、远红外波段空间分辨率为500m；其余29个波段空间分辨率为1km。(2)回访周期短。时间分辨率较高。Terra和Aqua都是太阳同步极轨卫星，可以得到每天最少2次白天和2次黑夜更新数据。(3)MODIS数据的全球免费接收政策，这样的数据接收和使用政策为科学研究提供了廉价并且实用的数据资源。本发明使用的MODIS产品是MOD15A2(空间分辨率为1km的逐8日产品数据)。

步骤1.2、获取CMIP6历史时段模拟数据和未来时段模拟数据：

在https://esgf-node.llnl.gov/search/cmip6/上下载CMIP6数据，CMIP6数据是模拟数据，其包括历史时段模拟数据和未来时段模拟数据；然后进行筛选、集合、订正；

CMIP6历史时段模拟数据和未来时段模拟数据是后续计算必须要获取到的数据。CMIP6模型内模式模拟的变量个数、时空分辨率都不同。

气候要素：为了降低降水量、太阳辐射量、最高温度、最低温度、相对湿度模拟结果的差异性，我们选择满足2个要求的模式下载：一是模式对降水量、太阳辐射量、最高温度、最低温度、相对湿度全部变量的逐日数据都有模拟且能下载下来；二是模式对历史时段模拟数据、未来(ssp126、ssp245、ssp585)情景的模拟数据都全能下载下来。

表1

叶面积指数：本实施方式选择了四种CMIP6模型叶面积指数产品(ACCESS-ESM1-5、CMCC-ESM2、MPI-ESM1-2-HR、MPI-ESM1-2-LR)参与对比计算，这是因为在当时研究阶段，具有历史时段模拟数据、未来(ssp126、ssp245、ssp585)情景的逐日模拟数据的只有这四个模式。

CO₂：选择了五种CMIP6模型CO₂产品(BCC-CSM2-MR、CESM2-WACCM、MRI-ESM2-0、NorESM2-LM、NorESM2-MM)参与对比计算，这是因为在当时的研究阶段具有历史时段模拟数据、未来(ssp126、ssp245、ssp585)情景的逐月模拟数据的只有这五个模式。

步骤1.3、对CMIP6模型中的多种模式数据进行尺度统一处理，得到空间分辨率(1km×1km)，时间分辨率(日)的气候模拟数据，空间分辨率(1km×1km)，时间分辨率(8日)的叶面积指数模拟数据，时间分辨率(月)的CO₂模拟数据：

网上下载下来的CMIP6数据时空分辨率有逐小时、日、月、年时间分辨率，空间上1000km、250km、50km、25km、1°，目前具有的那套历史时段模拟数据时空分辨率也都不一致，但是BEPS模型对于输入数据有要求，就需要把这些数据都统一成一样的，满足BEPS模型要求的格式，这样做有两个目的：一是不管是历史还是未来数据，它都是要输入到模型里，得满足模型的分辨率要求；二是对这些数据进行时空尺度对比的时候，分辨率一致才有可比性，才能在后期处理的时候，直接对他们进行运算和处理。

可以采取最近邻插值法对降水量、太阳辐射量、最高温度、最低温度、相对湿度、叶面积指数模式模拟结果进行空间插值处理；优点是计算量很小，算法也简单，因此运算速度较快，效率较高。

最邻近插值法是插值点的四个邻近点最近的点的数值作为插值点的数值。设插值点g(i，j)与周围四个距离最近的邻点f_k(i，j)，的距离为d_k，k＝1，2，3，4，当d_k＝min{d₁，d₂，d₃，d₄}时，g(i，j)＝f_k(i，j)。最邻近插值法假设之一是最接近于网格点p(x，y)的点的属性值就是该点的属性，要处理的节点值可以使用网格节点的最邻近值。同样对一些区域，比如填充无取值的数据也很有用。

CO₂的模拟结果不用去做插值，因为BEPS模型需要的就是一个月一个数，所以这一步就读取CO₂模拟格点数据，取在黑龙江省范围内的格点做平均计算即可。

步骤1.4、对CMIP6中多种模式历史时段模拟效果进行综合评估：

CMIP6中的数据存在多个模式，本实施方式的气候找到9个模式、叶面积指数找到4个模式、CO₂找到5个模式，如果只选一个模式往后计算，单一模式偏差大，如果全部拉进来平均，那就都给平滑了，这个因素平均9个模式的值，那个因素平均4个模式的值也不合适，而且这个数据如果一个模拟效果低于其他大部分模式的模式被放进来，一下子就影响模拟的准确性了。

应用泰勒图分析CMIP6中多模式模拟区域降水量、太阳辐射量、最高温度、最低温度、相对湿度、二氧化碳浓度、叶面积指数的适用性。泰勒图可以将模拟值与日本再分析数据的COR、RMSE、STD绘制在一张极坐标图上进行比较，评价各个模式的模拟效果。优点是可以直观的各个统计量的差异，一张图就集合了3个统计量，考虑的会比较全面，在对于多个模式进行对比时，直接看时空分布图，让人眼花缭乱看不出区别，选用泰勒图更加直观清晰。

三种统计量的定义公式为：

其中，COR代表空间相关系数，通过方位角表示；RMSE为均方根误差，通过与参考点的径向距离表示，均方根误差越小模拟效果越好；STD表示标准差，和/>为参考值和模拟值平均值，x_i，y_i为参考点和模拟点样本值，n为样本点数量。

本实施方式选取黑龙江省作为研究区域，使用线性插值方法，将所有数据插值到1km水平分辨率上，采用泰勒图评估模式的模拟能力，对于模拟较好的模型进行集合，并且进一步根据实测数据进行订正，用于对黑龙江省气候条件进行模拟。

图1-图6分别为降水、太阳辐射、平均温度、最高温度、最低温度、相对湿度的泰勒图；

BUOY指的是基准数据(即日本再分析数据的气候数据，AVHRR与MODIS的叶面积指数数据，CO₂数据)，基于表1中A-I对应的模式，进行综合评估判断：

降水：图1中各模型模拟结果的相关系数接近不纳入评价，A、C、D、H的标准差与BUOY最接近，对比均方根误差A、C、D模拟效果最好；

太阳辐射：图2中A、C、D的相关系数最高，均方根误差最小，C、D、E标准差与BUOY最接近，A、C、D模拟效果最好；

平均温度：图3中均方根误差、相关系数相接近，A、B、G标准差与BUOY最接近，A、B、G模拟效果最好；

最高温度：图4中相关系数接近，A、B均方根误差最小，F的标准差与BUOY最接近，ABI的模拟效果最好

最低温度：图5中相关系数接近，A、B均方根误差最小，G的标准差与BUOY最接近，ABG的模拟效果最好

相对湿度：图6中A、B、D、E、G均方根误差最小，C、D标准差与BUOY最接近，D、A、F、E相关系数最高，A、D、E的模拟效果最好

A出现6次，B出现3次，C出现2次，D出现3次，E出现1次，G出现2次，I出现1次，按照出现的次数，模拟效果最好的模型是A、B、D、C、G、E、I，前三名是A、B、D。

表2模拟效果较好的模式名称

叶面积指数和CO₂可用模式较少可以针对时空分布特征，进行筛选；

图7-图11为不同叶面积指数产品(ACCESS-ESM1-5、CMCC-ESM2、MPI-ESM1-2-HR、MPI-ESM1-2-LR)以及本发明的面积指数效果。图12为不同叶面积指数产品(ACCESS-ESM1-5、CMCC-ESM2、MPI-ESM1-2-HR、MPI-ESM1-2-LR)和本发明的叶面积指数曲线图。

图13为CO₂模拟结果(1961年1月-2015年12月)。

步骤1.5、选取CMIP6模型的最优模式进行集合：

基于前述步骤对于CMIP6模型内各模式的模拟效果有了一个基本的认识。接下来我们就要制定一个选哪个模式进行集合，怎么集合的规则，由此得到一个集合后的历史时段模拟数据、未来时段模拟数据，为下一步的订正做准备。

根据泰勒图列出模拟降水量、太阳辐射量、最高温度、最低温度、相对湿度、二氧化碳浓度、叶面积指数较优的3个模式，筛选出3个出现次数较高的模式，再分别对历史时段3个模式的模拟值。

降水量、太阳辐射量、最高温度、最低温度、相对湿度选择ACCESS-CM2、ACCESS-ESM1-5、FGOALS-g3做集合。

叶面积指数选择CMCC-ESM2、MPI-ESM1-2-HR、MPI-ESM1-2-LR做集合。

CO₂选择BCC-CSM2-MR、CESM2-WACCM、MRI-ESM2-0做集合。

针对CMIP6数据(历史时段模拟数据和未来时段模拟数据)，将3个模式的模拟值进行等权重集合(模式1+模式2+模式3)/3。多模式集合能较好地反映要素的时空分布特征，弥补了以往单一模式在评估时偏差较大的问题。

由于模型存在不确定性，在模拟的数值上与真实值存在偏差，这部分偏差不仅会导致输入数据不准、不适用于BEPS内部参数，使得最后输出的变量差别很大，因此要消除这部分偏差。本发明采用历史时段数据和历史时段集合模拟值(由CMIP6数据中的历史时段模拟数据计算得到)算一个差，把这个差叠加到未来集合模拟值，这样就得到了未来集合数据，这一数据将用于下一部分输入BEPS模型。

具体的计算公式是：

R_future(i，j，k)＝R_CMIP6future(i，j，k)-error(i，j)

其中，是1986-2015年CMIP6模拟逐日平均格点值(降水量、太阳辐射量、最高温度、最低温度、相对湿度、叶面积指数)或者逐月平均值(CO₂)，/>是1986-2015年历史时段数据的逐日平均格点值(降水量、太阳辐射量、最高温度、最低温度、相对湿度、叶面积指数)或者逐月平均值(CO₂)，error是二者的偏差；R_future(i，j，k)是未来集合模拟值最终订正后的结果，R_CMIP6future(i，j，k)是CMIP6中的未来模拟数据(2016-2060年)对应的逐日逐格点数据或逐月数据。

CMIP6模型预测过程中，由于排放场景的不确定以及模式本身的不确定性，导致预测结果存在较大的不确定性，在外强迫不变的情况下，耦合模式存在虚假的长期变化，为了减少这部分误差，需要对多模式集合平均数据进行漂移订正。具体的操作是将多模式集合平均的历史时期的模拟数据减去日本再分析数据的多年平均值，计算结果称为漂移误差，将未来各模式模拟的数据分别减去漂移误差，进行漂移订正。

订正处理后的未来气候数据模拟能力增强，以降水、辐射为例进行证明：

图14-图15分别为黑龙江2015年历史段对应的订正前、订正后的太阳辐射量的效果图。图16为黑龙江预估段太阳辐射量的效果图。

图17为黑龙江2015年太阳辐射量对比曲线图；图18为黑龙江2001-2014年评估段(仿真时间节点前的数据处理结果称为是评估)及2015-2060年预估段(仿真时间节点后的数据处理结果称为是预估)的逐年太阳辐射量对比曲线。

图19-图20分别为黑龙江2015年对应的订正前、订正后的降水量的效果图。图21为黑龙江2015年对应的预估段的降水量的效果图。

图22为黑龙江2015年降水量对比曲线图；图23为黑龙江2001-2014年评估段(仿真时间节点前的数据处理结果称为是评估)及2015-2060年预估段(仿真时间节点后的数据处理结果称为是预估)的逐年降水量对比曲线。

二、对BEPS模型输入目的地的历史数据，进行模型预热，并对历史GPP、NPP、NEP、SW、ET、释氧量、碳储量模拟，对历史碳吸收指标进行评估。所述目的地的历史数据，本质上就是目的地的历史时段数据，数据类型完全相同，只不过相对于步骤一的历史时段数据是一个时间范围尺度更广的历史时段，即步骤一所述历史时段数据是历史数据中靠近当前时间点的部分历史时段所对应的数据，例如历史时段数据为从2001开始的数据，而历史数据是从1901年开始的数据。

实际上本发明可以直接利用步骤一得到的R_future(i，j，k)进行碳吸收指标(包括部分表示碳循环的指标)预估，但是为了消除距离当前时间点一定时间范围内气候等因素对平衡态误差的影响，本发明选择以历史数据对碳库的平衡进行处理，能够有效提高对未来碳吸收指标预估的准确度。

本实施方式中，根据黑龙江省区域特征，进行碳循环机理模型参数的调整与架构优化，提高模型运行效率；参数调整与架构优化：修改代码及配置文件(预热和正式计算的起止年份，计算的起止行号和列号以及模型的运行场景)，将计算区域从全国改到黑龙江省提高了效率，根据黑龙江省的输入数据进行了调参，开发二进制结果转为可视化产品的功能，实现GPP、NPP、NEP、释氧量、ET、碳储量和土壤水分含量产品的可视化，包含栅格数据，统计表格，专题图。实现原始数据自动化加工成BEPS模型所需要的二进制格式。

图24-图25分别为2020年NPP、NEP效果图。图26-图27分别为1961-2019年NPP、NEP曲线图。

本实施方式中，对针对黑龙江省的碳循环机理模型，进行了改进。黑龙江省碳循环机理模型继承了碳循环机理模型模拟GPP的主要算法，改进增加了叶氮变化对最大羧化速率的动态修正(V_m＝V_cmax,25·f(T).f(N))。利用Ball_Berry模型通过对气孔导度的二次计算实现了生态系统碳水循环的耦合。通过自养呼吸Ra和异养呼吸Rh的计算实现了NPP和NEP的模拟。下面结合碳循环机理模型及改进进行详细说明：

(A)气孔导度参与下的GPP模拟

模型中以叶片尺度的Farquhar光化学模型为基础进行光合作用模拟，叶片瞬时光合速率为：

A＝min(W_c,W_j)-R_d

R_d＝0.0015V_m

其中，A为叶片光合速率；Wc和Wj分别为受羧化酶活性限制和光子传导限制的光合速率；Rd为叶片白天的暗呼吸；Vm为最大羧化速率，与植被类型、温度和叶片氮含量有关；Ci为叶肉细胞CO₂浓度；Г为暗呼吸时CO₂补偿点；K为酶促反应速度常数；J为电子传递速度。

V_m＝V_cmax,25·f(T).f(N)

其中，Vcmax，25为25℃时的最大羧化速率，f(T)和f(N)分别为温度和叶片氮含量的订正函数。

由于Farquhar模型为瞬时叶片尺度的光化学模型，不能直接应用于模拟每日的冠层光合作用总量(GPP)，黑龙江省碳循环机理模型基于光合速率和气孔导度在一天中呈现正弦变化的假设，对Farquhar模型进行时间尺度转换，然后计算每日的GPP。

每日的平均光合速率计算为：

其中，A是W_c或W_j。

对于W_c而言：

a＝(K+C_a)² (9)

b＝2(2Γ+K-C_a)·V_m+2(C_a+K)·R_d (10)

c＝(V_m-R_d)² (11)

对于W_j而言：

a＝(2.3Γ+C_a)² (12)

b＝0.4(4.3Γ-C_a)·J+2(C_a+2.3Γ)·R_d (13)

c＝(0.2J-R_d)² (14)

其中，C_a为大气中CO₂的浓度；g_min是最小气孔导度；g_n是中午时的气孔导度。

黑龙江省碳循环机理模型分两步进行气孔导度的模拟计算，以实现碳水循环的耦合。首先利用Jarvis模型计算每天气孔导度的初始值，进而计算光合速率。然后将计算的光合速率代入Ball-Berry模型得到土壤水分修正后的最终气孔导度。

其中，a为与植被类型有关的参数，表示气孔导度与光合速率的耦合程度；f(w)为土壤水分修正因子；VPD为叶片表面饱和水气压差；D为气孔导度对VPD变化的敏感系数。

(B)、自养呼吸计算与NPP的模拟

NPP计算为GPP与自养呼吸Ra之差。其中，Ra计算为生长呼吸Rg和维持性呼吸Rm之和，Rm表示为碳库大小、参考温度下的呼吸速率和温度的函数，即：

NPP＝GPP-R_g-R_m

其中，Rg是生长呼吸，取为0.25GPP；C_α是4个植被碳库的大小，每天更新；α＝1，2，3，4；1代表叶、2代表茎、3代表粗根、4代表细根；r_α是参考呼吸速率，为模型参数，在子程序readb()中定义；Ts表示温度；Tb表示参考温度；Q表示呼吸对温度的响应函数。

(C)、异养呼吸计算与NEP的模拟

NEP计算为NPP与自养呼吸Rh之差，Rh的计算采用CENTURY模型的方法进行：

NEP＝NPP-R_h

其中，C_β是9个土壤碳库的大小；β＝1，2，3，…，9，分别为①地表凋落物结构性部分碳库，②地表凋落物易分解部分碳库，③地表微生物碳库，④茎和粗根凋落物碳库，⑤土壤凋落物结构性部分碳库，⑥土壤凋落物易分解部分碳库，⑦土壤微生物碳库，⑧土壤慢分解碳库，⑨土壤惰性碳库；k_β为各个碳库的参考呼吸速率；Ts表示温度；f(θ)土壤水分对凋落物和土壤有机碳分解影响的调节因子；Te表示土壤质地。

针对上述植物碳库和土壤碳库进行动态更新：

基于平衡态假设驱动碳循环机理模型进行土壤碳库的初始化，确定碳库初始值。通过计算的碳通量和异养呼吸环境调节因子，进行碳库更新。

每年的碳库更新为：

C_j(i)＝C_j(i-1)+ΔC_j(i)

对于植被碳库：

其中，ΔCj(i)为第j个植被碳库在第i年的变化量；NPP(i)为第i年森林的NPP；fj表示NPP分配到第j个植被碳库的比例系数；kj表示第j个生物量碳库的周转率；Dj(i)为森林扰动导致的第j个植被碳库的直接碳排放量；ΔCj(i)为第j个植被碳库在第i年的变化量。

对于土壤碳库：

其中，ΔC_k(i)为第k个植被碳库在第i年的变化量；ζ_m,j为从第m个碳库分解后向第j个碳库的转移速率；ζ_j为第j个土壤碳库的分解系数；D_j(i)为森林扰动导致的第k个土壤碳库的直接碳排放量。

在模拟土壤碳循环过程的同时，模拟土壤氮循环过程。氮随碳在不同的碳库中转移，在这一过程中部分氮矿化为有机氮进入周转碳库，供植被吸收利用。

三、利用BEPS模型未来预估碳吸收变化特征

步骤一获取了未来预估的降水量、太阳辐射量、最高温度、最低温度、相对湿度、二氧化碳浓度、叶面积指数数据，并对其进行了漂移订正，就是为了这一步我们把这些数据输入BEPS模型，得到一套新的未来的GPP(总初级生产力)、NPP(净初级生产力)、NEP(净生态系统生产力)、SW(土壤水分含量)、ET(蒸散发)、释氧量、碳储量。

步骤3.1、GPP、NPP、NEP、SW、ET、释氧量、碳储量预估。

将漂移订正后的未来预估的降水量、太阳辐射量、最高温度、最低温度、相对湿度、二氧化碳浓度、叶面积指数数据输入BEPS模型，实际上是基于历史数据模型输入只修改这些数据的值，其余的氮沉降、地表覆盖、30年平均温度、土壤质地、聚集度指数、经纬度数据使用的是预测时段前一年(或者说历史数据中最后一年)的数据。根据黑龙江省区域特征，输入预测的气象、二氧化碳浓度与叶面积指数驱动数据，其他对模型影响较小的驱动数据，复用2020年结果，基于上述驱动数据模拟获取黑龙江省精细化GPP、NPP、NEP、SW、ET、释氧量、碳储量预估数据。

图28为GPP的估计曲线图，图29为NPP的估计曲线图，图30为NEP的估计曲线图。

步骤3.2、黑龙江省GPP、NPP、NEP、SW、ET、释氧量、碳储量预估数据时空分布特征评估：

前面得到了一套新的未来的GPP、NPP、NEP、SW、ET、释氧量、碳储量，但是他们的变化规律，突变特征尚没有办法确定，所以势必要采取一些措施，我们来评估一下这些数据的时空变化特征。

本发明采用一元线性回归法，可直接进行预估碳吸收指标变化率计算和显著性分析；并利用Mann-Kendall(简称MK)趋势检验法进行突变分析实现显著性的检验。对比碳循环模型预估数据与评估数据，分析未来黑龙江省碳收支的变化特征，评估模型模拟结果的准确度。

本发明采用一元线性回归法，可直接进行碳吸收指标数据变化率计算和显著性分析。

通过最小二乘法计算碳吸收指标数据与时间的一元线性回归系数的方法：

y＝ax+b

其中，x表示时间变量，y表示样本量为n的因变量(得到的GPP、NPP、NEP、SW、ET、释氧量、碳储量)，a表示回归系数，当a>0表示因变量随时间呈增加趋势，当a<0表示因变量随时间呈减少趋势，a越大表示变化速率越快。也可以使用a的10倍来表示因变量随时间的变化趋势。

利用Mann-Kendall(简称MK)趋势检验法进行突变分析。

在时间序列为随机的假设下，

/>

若|UF_k|＞α，则表明序列有明显的趋势变化。

式中，秩序列s_k是第i时刻数值大于j时刻数值个数的累计数，a_ij表示序列内两个数的大小，z为研究序列(得到的GPP、NPP、NEP、SW、ET、释氧量、碳储量)，UF1＝0，E(s_k)，var(s_k)是累计数s_k的均值和方差。

再根据时间序列x的逆序，重复上述过程，求出UB_k。若UF_k和UB_k出现交点，且交点处介于临界值之间，那么交点所对应时刻即突变开始时间。α＝0.05时，临界值为±1.96。

同时，本发明还提出一种局部时间区域范围内趋势分析的方法，对BEPS输出变量进行短时间内的精细分析，由于全球气候处于持续变暖的影响，为了更好地对短时间内(减小甚至消除长时间范围尺度下的部分因素的影响)的因素变化影响的时空分布特征进行评估，因此以年代纪为为单位进行评估，其中初始年代纪为预估第一年度所在的按照划分年代纪宽度划分的第一年代纪。以GPP序列为例：

以每年12月31日GPP作为年值，计算每10年的空间平均值，得到黑龙江省空间分辨率1km的分布栅格数据，2021-2030年、2031-2040年、……、2051-2060年，得到4组栅格数据，用GPP₁、GPP₂、GPP₃、GPP₄来表示，针对逐栅格进行计算。GPP_k(i，j)表示各栅格数值，(i，j)是表示栅格数据中的某一像元点的位置(经纬度坐标)。

其中，Z(i，j)为表示GPP年代纪特征，针对初始年代纪(本实施例中为2021-2030年)初始化为0。

最后得到Z指标，可以根据Z指标的空间分布特征，直观分析黑龙江省不同地区GPP乃至其他输出变量的年代际变化特征，也可以基于不同土地覆盖类型的Z值评估不同生态系统的GPP年代际变化特征。

本发明的上述算例仅为详细地说明本发明的计算模型和计算流程，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。

Claims

1.基于全球气候模式驱动碳循环机理模型的碳吸收预估方法，其特征在于，包括以下步骤：

S1、利用CMIP6模型预测未来驱动数据变化特征：

步骤1.1、收集目的地历史时段驱动碳循环机理模型所需驱动数据，简记为目的地历史时段数据：

步骤1.3、对CMIP6中的多种模式数据进行尺度统一处理；对CMIP6中的多种模式数据进行尺度统一处理时，气候要素数据的尺度：空间分辨率为1km×1km，时间分辨率为日；叶面积指数的尺度：空间分辨率为1km×1km，时间分辨率为8日；二氧化碳浓度的尺度：时间分辨率为月，采取最近邻插值法对空间尺度进行统一，即最近邻插值法统一空间分辨率；

R_future(i，j，k)＝R_CMIP6future(i，j，k)-error(i，j)

其中，是CMIP6历史时段集合模拟数据平均值，其中降水量、太阳辐射量、最高温度、最低温度、相对湿度、叶面积指数为逐日平均格点值，CO₂浓度为逐月平均值；是目的地历史时段数据的平均值，降水量、太阳辐射量、最高温度、最低温度、相对湿度、叶面积指数为逐日平均格点值，CO₂浓度为逐月平均值；error(i，j)是二者的偏差；R_future(i，j，k)是未来集合模拟值订正后的结果，R_CMIP6future(i，j，k)是CMIP6未来时段集合模拟数据，如果对应降水量、太阳辐射量、最高温度、最低温度、相对湿度、叶面积指数，R_future(i，j，k)是各自对应的逐日平均格点值，如果对应CO₂浓度，R_future(i，j，k)是对应的逐月值；

2.根据权利要求1所述基于全球气候模式驱动碳循环机理模型的碳吸收预估方法，其特征在于，S2中基于未来集合模拟值订正后的降水量、太阳辐射量、最高温度、最低温度、相对湿度、二氧化碳浓度、叶面积指数数据，利用BEPS模型预估未来碳吸收指标之前，需要对BEPS模型进行模型预热，即：

输入目的地的历史数据，进行历史碳吸收指标模拟；

3.根据权利要求1或2所述基于全球气候模式驱动碳循环机理模型的碳吸收预估方法，其特征在于，利用BEPS模型预估未来碳吸收指标之后，对碳吸收指标进行时空分布特征评估，具体包括以下步骤：

4.根据权利要求3所述基于全球气候模式驱动碳循环机理模型的碳吸收预估方法，其特征在于，利用BEPS模型预估未来碳吸收指标之后，还要对模型模拟结果的准确度进行评估，具体包括以下步骤：

5.根据权利要求3所述基于全球气候模式驱动碳循环机理模型的碳吸收预估方法，其特征在于，步骤1.4中采用泰勒图对CMIP6中多种模式历史时段模拟效果进行综合评估。

6.根据权利要求5所述基于全球气候模式驱动碳循环机理模型的碳吸收预估方法，其特征在于，采用泰勒图对CMIP6中多种模式历史时段模拟效果进行综合评估的过程包括以下步骤：

其中，COR为空间相关系数，RMSE为均方根误差，STD为标准差；

7.根据权利要求6所述基于全球气候模式驱动碳循环机理模型的碳吸收预估方法，其特征在于，针对目的地的历史时段数据中的叶面积指数采用源于AVHRR和MODIS的叶面积指数；源于AVHRR和MODIS的叶面积指数需要经过异源订正处理，方法如下：

R_NEWLAIMODIS(i，j，k)＝R_LAIMODIS(i，j，k)+error(i，j)

其中，是第一历史时段的AVHRR中的逐格点叶面积指数数据的平均值，是第一历史时段的MODIS中逐格点叶面积指数数据的平均值，error(i，j)是二者的偏差；R_NEWLAIMODIS(i，j，k)是订正后的MODIS叶面积指数的结果；R_LAIMODIS(i，j，k)为订正前的MODIS叶面积指数的结果；

所述第一历史时段为历史时段期间内部分历史时段。