CN115952702B - 一种基于forcchn模型和遥感数据的森林nep的计算方法 - Google Patents

Abstract

本发明涉及一种基于FORCCHN模型和遥感数据的森林NEP的计算方法，包括：S1、对目标区域进行栅格化处理；S2、根据目标区域和目标时间区间，获取包括植被数据的多种数据，植被数据是通过基于遥感数据和地面高程数据的反演过程得到的；S3、对预先搭建的FORCCHN模型进行参数优化，得到本地化的FORCCHN模型；S4、将本地化的FORCCHN模型运行至生态平衡态；S5、运行FORCCHN模型进行动态模拟，得到在目标时间区间内，目标区域中每个栅格逐日的GPP_d、R_m、R_g、R_s；S6、基于逐日的GPP_d、R_m、R_g、R_s，计算得到目标区域在目标时间区间内的NEP。本发明的方法融合遥感数据与地面高程数据，能够提高计算所得NEP的精确度。

Description

一种基于FORCCHN模型和遥感数据的森林NEP的计算方法

技术领域

本发明涉及生态环境技术领域，尤其涉及一种基于FORCCHN模型和遥感数据的森林NEP的计算方法。

背景技术

森林作为陆地生态系统的主体，储存了陆地生态系统约45％的有机碳，具有碳汇量大、成本低、生态附加值高等特点。净生态系统生产力(Net Ecosystem Productivity，NEP)表示较大尺度生态系统中碳的净贮存量，是生态系统植被净初级生产力和土壤异氧呼吸之差，在全球尺度上直接揭示陆地生态系统与大气系统之间的二氧化碳交换程度，即碳平衡进程。目前，充分发挥森林生态系统固碳作用能够降低温室气体浓度和抑制全球变暖趋势，被认为是当前减缓全球气候变化的一种较好的途径。随着全球气候变化问题的日益突出，对森林净生态系统生产力的评估已经引起广泛关注。

森林净生态系统生产力的传统计算方法主要有森林定位清查法、涡度相关通量观测法、同位素示踪法等方法。然而，由于受到森林区域环境条件、气候条件和植被类型等因素的影响，以及受限于传统的抽样调查、定点观测法调查、观测数量和资金等原因，传统计算方法将很难在区域和全球尺度上直接和全面地测量森林的净生态系统生产力，从而使模型方法成为陆地碳循环研究中不可替代的手段，但模型方法依然需要采集相关计算区域的植被数据才能进行准确的模拟和计算，因而基于遥感技术采集相关的植被数据来计算森林净生态系统生产力的方法逐渐开始引起人们的关注，但遥感技术采集的植被数据通常是基于二维平面的，而实际的地形则存在高低起伏，使实际的植被采光面积与遥感技术采集的数据出现差异，影响模型方法最终计算的森林净生态系统生产力的准确性。

发明内容

(一)要解决的技术问题

鉴于现有技术的上述缺点、不足，本发明提供一种基于FORCCHN模型(FORestecosystem Carbon budget model For CHiNa，中国森林生态系统碳收支模型)和遥感数据的森林NEP的计算方法，其解决了遥感技术采集的植被数据与实际植被数据存在差异，降低模型方法计算得到的NEP的精确度的技术问题。

(二)技术方案

为了达到上述目的，本发明采用的主要技术方案包括：

本发明实施例提供一种基于FORCCHN模型和遥感数据的森林NEP的计算方法，包括：

S1、对目标区域按照预设尺寸进行栅格化处理；

S2、根据目标区域和目标时间区间，获取本地化数据、初始数据和驱动数据；

所述本地化数据为所述目标区域中，每个栅格任意时间区间内的气象数据、土壤数据、植被数据、地理数据和碳通量数据；

所述初始数据为所述目标时间区间的前一年内，所述目标区域中每个栅格的气象数据、土壤数据、植被数据和地理数据；

所述驱动数据为所述目标时间区间内，所述目标区域中每个栅格的气象数据、土壤数据、植被数据和地理数据；

其中，所述植被数据是通过基于遥感数据和地面高程数据的反演过程得到的；

S3、基于所述本地化数据对预先搭建的FORCCHN模型进行参数优化，得到本地化的FORCCHN模型；

S4、将所述初始数据输入本地化的FORCCHN模型运行，得到平衡态的FORCCHN模型；

S5、将所述驱动数据输入平衡态的FORCCHN模型进行动态模拟，得到在所述目标时间区间内，所述目标区域中每个栅格逐日的森林总初级生产力(Gross PrimaryProductivity，GPP)GPP_d、维持呼吸量R_m、生长呼吸量R_g、土壤异养呼吸量R_s；

S6、基于每个栅格所述逐日的GPP_d、R_m、R_g、R_s，计算得到目标区域在目标时间区间内的森林净生态系统生产力NEP。

本发明实施例提出的方法中，基于气象数据、土壤数据、植被数据和地理数据，使用FORCCHN模型计算逐日的森林总初级生产力GPP_d、维持呼吸量R_m、生长呼吸量R_g、土壤异养呼吸量R_s，进而计算得到森林净生态系统生产力NEP。其中，上述植被数据是通过基于遥感数据和地面高程数据的反演过程得到，从而使基于遥感技术采集的数据与地面高程数据融合，使植被数据更贴合实际情况，为模型提供更为准确的输入数据。此外，上述FORCCHN模型是一个基于个体生长过程的斑块模型，其采用逐个模拟树木碳过程的方法，获得碳循环过程中碳通量和碳库量的计算结果，能够进一步提高NEP的精确度。

可选地，在S2中，所述气象数据包括：每个栅格上逐日的平均气温、最高气温、最低气温、降水量、日照时数、平均气压、风速、相对湿度、地表太阳辐射，以及气象站点的经纬度和海拔高度；

所述土壤数据包括：每个栅格上的土壤砂粒含量、土壤粉粒含量、土壤粘粒含量、土壤容重、土壤田间持水量、土壤萎蔫含水量、土壤碳密度、土壤氮密度；

所述植被数据包括：每个栅格上的森林类型、逐日的叶面积指数LAI(LeafAreaIndex)、最大叶面积指数LAI_max；

所述地理数据包括：所述目标区域栅格化后每个格点的经纬度和地面高程数据；

所述碳通量数据包括已公开的碳通量实测数据，和/或公开的模拟数据。

可选地，所述植被数据中，所述逐日的叶面积指数、最大叶面积指数LAI_max为通过基于遥感数据和地面高程数据的反演过程得到的；

所述遥感数据为公开的NDVI(Normalized Difference Vegetation Index，归一化植被指数)数据，所述基于遥感数据和地面高程数据的反演过程包括：

A1、获取目标时间区间内，每一栅格上逐日的NDVI数据；

A2、基于地理数据中的地面高程数据，根据公式(1)计算每一栅格上的地形起伏度指数R，所述公式(1)为：

公式(1)中，H_max表示该栅格中的最大高程值，H_min表示该栅格中的最小高程值，ΔH表示目标区域内的最大高程值与最小高程值的差值，m表示一个与地形起伏度指数有关的经验常数；

A3、基于一栅格上的R和逐日的NDVI数据，根据公式(2)得到每一栅格上逐日的光合有效辐射吸收比例FPAR(Fraction of Photosynthetically Active Radiation)，所述公式(2)为：

公式(2)中，NDVI表示每一栅格上逐日的NDVI数据，NDVI_max表示与森林类型相关的最大NDVI值，NDVI_min表示与森林类型相关的最小NDVI值，NDVI_max和NDVI_min的值是与森林类型相关的一个经验值；

A4、基于逐日的FPAR，比较得到目标时间区间中第一自然年中最大的FPAR值FPAR_max和最小的FPAR值FPAR_min；

A5、基于FPAR_max、FPAR_min和逐日的FPAR，根据公式(3)计算LAI_max、LAI_min、逐日的LAI，所述公式(3)为：

公式(5)中，LAI_v表示与森林类型相关的一个经验值；

将FPAR_max和逐日的FPAR带入公式(3)，得到逐日的LAI；

将FPAR_max和FPAR＝FPAR_min带入公式(3)，得到LAI_max。

可选地，所述气象数据为基于中国气象局公开的实测气象数据，根据自然邻域法和反向距离加权法，以预设尺寸为基本空间单元进行空间插值得到的与栅格一一对应的气象数据。

可选地，所述土壤数据为基于公开的全国1：100万数字化土壤图，以预设尺寸为基本空间单元，经过空间化处理得到的与栅格一一对应的土壤数据。

优选地，在S3中，所述预先搭建的FORCCHN模型包括：

初始化过程，用于根据植被数据确定每一栅格上的树木株数、森林类型、树木初始基径、树木初始树高和树木生理生态参数；

光合作用过程，用于基于气象数据和植被数据模拟树木的光合作用过程，逐日计算每个栅格上树木的总光合，得到每个栅格逐日的森林总初级生产力GPP_d；

呼吸作用过程，用于基于气象数据和植被数据模拟树木的呼吸作用过程，逐日计算得到每个栅格上树木的维持呼吸量R_m和生长呼吸量R_g；

凋落过程，用于计算树木的凋落通量；

光合产物分配过程，用于基于光合作用过程、呼吸作用过程和凋落过程，计算光合产物向多个碳库的分配量；

土壤有机物呼吸和转移过程，用于基于气象数据、土壤数据，以及凋落过程，模拟土壤中逐日的呼吸、转移和氮的矿化过程，得到土壤逐日的呼吸排量，基于土壤逐日的呼吸排量计算得到土壤逐日的土壤异养呼吸量R_s；

年生长过程，用于基于光合产物分配过程和凋落过程，模拟树木的年生长过程，计算树木基径增量和树木树高增量，将树木初始基径与树木基径增量之和作为下一自然年的树木初始基径，将树木初始树高与树木树高增量之和作为下一自然年的树木初始树高。

优选地，所述S3包括：

B1、将本地化数据中的气象数据、土壤数据、植被数据和地理数据输入预先搭建的FORCCHN模型运行，基于FORCCHN模型输出的结果计算得到模拟结果数据；

其中，所述模拟结果数据包括模拟的森林净初级生产力NPP_simulated、模拟的森林净生态系统生产力NEP_simulated、模拟的森林净生态系统碳交换量NEE_simulated中的至少一种；

B2、基于本地化数据中的碳通量数据得到与模拟结果数据对应的实际结果数据；

其中，所述实际结果数据包括实际的森林净初级生产力NPP_true、净生态系统生产力NEP_true、净生态系统碳交换量NEE_true中的至少一种；

B3、将所述模拟结果数据与实际结果数据比较；若两者的差异在预设范围内，则将当前的FORCCHN模型作为本地化的FORCCHN模型；

若两者的差异超过预设范围，则循环执行调参操作，直至两者的差异在预设范围内，得到本地化的FORCCHN模型；

所述调参操作包括：

调整FORCCHN模型中的树木生理生态参数，将本地化数据中的气象数据、土壤数据、植被数据和地理数据输入调整树木生理生态参数后的FORCCHN模型中运行，基于输出的结果计算得到新的模拟结果数据，将所述新的模拟结果数据与实际结果数据比较。

优选地，所述S4包括：

将所述初始数据输入本地化的FORCCHN模型运行，直至所述本地化的FORCCHN模型达到生态系统平衡态；其中，所述生态系统平衡态指基于FORCCHN模型运行并计算得到的NEP为零的状态。

优选地，所述目标时间区间为一个或多个连续的自然年，所述S6包括：基于每个栅格所述逐日的GPP_d、R_m、R_g、R_s，根据公式(4)得到每个栅格逐年的年森林净生态系统生产力NEP_grid，对所有栅格相同自然年的NEP_grid求和，得到目标区域在目标时间区间内逐年的NEP；

所述公式(4)为：

公式(4)中，NPP表示每个栅格的年森林净初级生产力(Net PrimaryProductivity，NPP)，NPP根据公式(5)得到，所述公式(5)为：

(三)有益效果

在本发明实施例中，通过基于遥感数据和地面高程数据的反演过程得到植被数据，结合气象数据、土壤数据和地理数据等背景数据，为FORCCHN模型提供更为准确的输入数据；由于遥感数据本身具有较高的分辨率，再与地面高程数据结合，根据地面的地形起伏度指数调整基于遥感数据的计算得到的光合有效辐射吸收比例，从而使得反演得到的叶面积指数、最大叶面积指数更为贴合实际的植被光合有效辐射吸收情况，进而为FORCCHN模型提供了更为准确的输入数据。

本发明实施例使用的FORCCHN模型是一个基于个体生长过程的斑块模型，其采用逐个模拟树木碳过程的方法，获得碳循环过程中碳通量和碳库量的计算结果；其基于气象数据、土壤数据、植被数据和地理数据，计算逐日的森林总初级生产力GPP_d、维持呼吸量R_m、生长呼吸量R_g、土壤异养呼吸量R_s，进而计算得到森林净生态系统生产力NEP，能够进一步提高NEP的精确度。

附图说明

图1为实施例中提供的一种基于FORCCHN模型和遥感数据的森林NEP的计算方法的流程示意图；

图2为实施例中FORCCHN模型的框架结构示意图；

图3为实施例中东北地区NPP实际观测值与若干次调参操作后的FORCCHN模型模拟计算得到的NPP模拟值的对比示意图；

图4为实施例中东北地区NEE观测值与本地化的FORCCHN模型模拟计算得到的NEE模拟值的对比示意图；

图5为实施例中模拟计算得到的东北地区森林NEP空间分布和变化趋势图。

具体实施方式

为了更好的理解上述技术方案，下面将参照附图更详细地描述本发明的示例性实施例。虽然附图中显示了本发明的示例性实施例，然而应当理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更清楚、透彻地理解本发明，并且能够将本发明的范围完整的传达给本领域的技术人员。

实施例一

如图1所示，本实施例提供一种基于FORCCHN模型和遥感数据的森林NEP的计算方法，包括：

S1、对目标区域按照预设尺寸进行栅格化处理。

具体地，所述栅格化处理可基于经纬度进行划分，也可基于直线距离进行划分，所述预设尺寸可根据目标区域的大小进行设置，栅格的预设尺寸越小，则对目标区域的划分越精细，本方法最终计算得到的NEP越精确，但计算量也会随预设尺寸的减小而增加，较佳地，所述预设尺寸可为1km×1km、2km×2km、5km×5km、10km×10km、20km×20km、50km×50km、100km×100km。

S2、根据目标区域和目标时间区间，获取本地化数据、初始数据和驱动数据。

所述本地化数据为所述目标区域中，每个栅格任意时间区间内的气象数据、土壤数据、植被数据、地理数据和碳通量数据；所述初始数据为所述目标时间区间的前一年内，所述目标区域中每个栅格的气象数据、土壤数据、植被数据和地理数据；所述驱动数据为所述目标时间区间内，所述目标区域中每个栅格的气象数据、土壤数据、植被数据和地理数据。

具体地，所述气象数据包括：每个栅格上逐日的平均气温、最高气温、最低气温、降水量、日照时数、平均气压、风速、相对湿度、地表太阳辐射，以及气象站点的经纬度和海拔高度。

所述土壤数据包括：每个栅格上的土壤砂粒含量、土壤粉粒含量、土壤粘粒含量、土壤容重、土壤田间持水量、土壤萎蔫含水量、土壤碳密度、土壤氮密度。

所述植被数据包括：每个栅格上的森林类型、逐日的叶面积指数LAI、最大叶面积指数LAI_max。

所述地理数据包括：所述目标区域栅格化后的地面高程数据和每个格点的经纬度。

所述碳通量数据包括：已公开的碳通量实测数据，和/或公开的模拟数据。在实际应用中，所述已公开的碳通量实测数据通常包括实测区域的NPP数据、NEP数据、NEE(NetEcosystem Exchange，净生态系统碳交换量)数据中的一种或多种；所述公开的模拟数据是指其他研究者公开的、基于其他模型方法模拟得到的NPP数据、NEP数据或NEE数据。

需要说明的，上述气象数据、土壤数据、地理数据、碳通量数据、遥感数据可从公开的渠道下载、购买获得，或通过公开的文献获得。其中，所述遥感数据可为NDVI数据、FPAR数据或其他反映植被特征的数据。

其中，所述植被数据是通过基于遥感数据和地面高程数据的反演过程得到的。

S3、基于所述本地化数据对预先搭建的FORCCHN模型进行参数优化，得到本地化的FORCCHN模型。

所述本地化数据用于输入预先搭建的FORCCHN模型运行，得到模拟输出结果，通过比较模拟输出结果和实测结果评价FORCCHN模型，根据评价结果调整FORCCHN模型的参数，直至评价结果达到预设的要求，即，本步骤实际上是一个对FORCCHN模型进行参数优化的过程，用于得到一个具有适配参数的本地化的FORCCHN模型。

需要说明的是，本实施例所使用的预先搭建的FORCCHN模型为已经公开的模型，具体地，如图2所示，所述预先搭建的FORCCHN模型包括：

初始化过程，用于根据植被数据确定每一栅格上的树木株数、森林类型、树木初始基径、树木初始树高和树木生理生态参数；其中，所述森林类型可根据公开的中国植被类型空间分布数据得到，所述树木株数、树木初始基径、树木初始树高则是基于植被数据中的最大叶面积指数LAI_max和经验方程法推算得到的。具体地，所述经验方程可为异速生长方程法，异速生长方程法(allometric equation)是指选择具有代表性的各个树种，测定如胸径、基面积、边材面积、树高、树冠体积和冠层郁闭等因子，通过建立它们与LAI之间的经验关系来估算LAI。本实施例采用随机数法或自然稀疏-3/2规律，针对栅格中的树木一株一株地产生其胸径，基于异速生长方程法，根据森林类型和胸径计算一株一株树木的叶面积指数，直至该栅格中所有树木的叶面积指数之和等于植被数据中的最大叶面积指数LAI_max，则当前树木的株数即为该栅格中的树木株树；此外，树木初始基径、树木初始树高也可基于上述异速生长方程法求得。

光合作用过程，用于基于气象数据和植被数据模拟树木的光合作用过程，逐日计算每个栅格上树木的总光合，得到每个栅格逐日的森林总初级生产力GPP_d。

呼吸作用过程，用于基于气象数据和植被数据模拟树木的呼吸作用过程，逐日计算得到每个栅格上树木的维持呼吸量R_m和生长呼吸量R_g。

凋落过程，用于计算树木的凋落通量。

光合产物分配过程，用于基于光合作用过程、呼吸作用过程和凋落过程，计算光合产物向多个碳库的分配量。

土壤有机物呼吸和转移过程，用于基于气象数据、土壤数据，以及凋落过程，模拟土壤中逐日的呼吸、转移和氮的矿化过程，得到土壤逐日的呼吸排量，基于土壤逐日的呼吸排量计算得到土壤逐日的土壤异养呼吸量R_s。

年生长过程，用于基于光合产物分配过程和凋落过程，模拟树木的年生长过程，计算树木基径增量和树木树高增量，将树木初始基径与树木基径增量之和作为下一自然年的树木初始基径，将树木初始树高与树木树高增量之和作为下一自然年的树木初始树高。

上述FORCCHN模型通过两种时间步长进行计算，对于每株树木而言，每天计算一次总光合、维持呼吸量、生长呼吸量、以及部分凋落物通量，光合产物分配和土壤有机物呼吸和转移过程也都是以日为步长发生的。FORCCHN模型假定每株树木每日内的维持呼吸量、部分叶、细根凋落、土壤有机质分解和氮的矿化速率恒定不变，从而避免在较短的步长上频繁计算相关通量。FORCCHN模型所计算的树木个体碳氮的年净增加量，需要加上前一年储存的碳量，在每年最后一天分配到不同的树木成分内。当树木没有足够的碳用于光合作用产生叶片时，认为树木已经死亡并停止生长，从而树木的碳储存下降。死亡的树木所包含的碳氮含量，以及活的木质部分，花，果，叶和根的凋落物通量在每年最后一天加入到土壤枯枝落叶库中。FORCCHN模型以个体为基础的斑块模型模拟所有时间内每株树木的出生、生长和死亡的全过程，每株树木的生长是通过对光照、水分及土壤的养分竞争获得，模拟的森林动态演替过程及其碳收支变化相对合理，能够较为准确地计算GPP_d、R_m、R_g、R_s，从而提高最终计算的NEP的精确度。

由于FORCCHN模型是一个逐日模拟森林生态系统碳收支过程的模型，因而基于该模型，可根据用户不同的需求计算不同时间步长的数值，具体地，可根据上述逐日的GPP_d、R_m、R_g、R_s数据，计算逐日、逐月、逐年的GPP_d、R_m、R_g、R_s、NPP、NEP、NEE等数据。

S4、将所述初始数据输入本地化的FORCCHN模型运行，得到平衡态的FORCCHN模型。

本步骤使用初始数据使FORCCHN模型运行至平衡态。区别于常规的基于平衡态假设设定的生态系统状态变量初始值，本实施例中的所述初始数据为目标时间区间的前一年内，所述目标区域中每个栅格的气象数据、土壤数据、植被数据和地理数据，将上述初始数据输入FORCCHN模型中运行，直至FORCCHN模型的生态系统达到平衡态，即，NEP＝0的状态，因而能够使平衡态的FORCCHN模型在后继基于驱动数据进行模拟时更为贴合生态系统的实际情况，进一步提升FORCCHN模型模拟得到的结果的准确性。

S5、将所述驱动数据输入平衡态的FORCCHN模型运行，得到在所述目标时间区间内，所述目标区域中每个栅格逐日的森林总初级生产力GPP_d、维持呼吸量R_m、生长呼吸量R_g、土壤异养呼吸量R_s。

S6、基于每个栅格逐日的GPP_d、R_m、R_g、R_s，计算得到目标区域在目标时间区间内的森林净生态系统生产力NEP。

本发明实施例提出的方法中，基于气象数据、土壤数据、植被数据和地理数据，使用FORCCHN模型计算逐日的森林总初级生产力GPP_d、维持呼吸量R_m、生长呼吸量R_g、土壤异养呼吸量R_s，进而计算得到森林净生态系统生产力NEP。其中，上述植被数据是通过基于遥感数据和地面高程数据的反演过程得到，从而使基于遥感技术采集的数据与地面高程数据融合，使植被数据更贴合实际情况，为模型提供更为准确的输入数据。此外，上述FORCCHN模型是一个基于个体生长过程的斑块模型，其采用逐个模拟树木碳过程的方法，获得碳循环过程中碳通量和碳库量的计算结果，能够进一步提高NEP的精确度。

实施例二

为了更好地理解实施例一，本实施例以中国东北地区为目标区域，以1981～2002年为目标时间区间，结合具体的步骤进行详细说明。

需要说明的是，本实施例中的中国东北地区包括辽宁、吉林、黑龙江三省和内蒙古自治区东部的呼伦贝尔市、兴安盟、通辽市、赤峰市。北起黑龙江主航道中心线(53°34′N)，南抵辽东半岛南端(38°43′N)，跨纬度14°51，西起内蒙古新巴尔虎右旗以西与蒙古国交界处(115°37′E)，东至黑龙江省抚远以东、乌苏里江汇入黑龙江处的耶字碑东角(135°5′E)，跨经度19°28′。

本实施例提供一种基于FORCCHN模型和遥感数据的森林NEP的计算方法，包括下述步骤：

S1、对目标区域按照10km×10km进行栅格化处理。

S2、根据目标区域和目标时间区间，获取本地化数据、初始数据和驱动数据。

所述本地化数据为目标区域中，每个栅格任意时间区间内的气象数据、土壤数据、植被数据、地理数据和碳通量数据；本实施例中，采用1989～1993年中国林业部森林普查所获得的690个点的NPP数据(罗天祥，1996)，从中选出属于东北地区的站点的NPP数据作为本地化数据中的碳通量数据，同时将1989～1993年的气象数据作为本地化数据中的气象数据，将1989～1993年的植被数据作为本地化数据中的植被数据。需要说明的是，受目标时间区间实测碳通量数据稀少的限制，本实施中采用NPP数据作为对FORCCHN模型进行参数优化的验证数据，也即S3中的实际结果数据，若目标时间区间的实测碳通量数据较为丰富时，也可采用NEP数据、NEE数据作为FORCCHN模型进行参数优化的验证数据。

所述初始数据为所述目标时间区间的前一年内，所述目标区域中每个栅格的气象数据、土壤数据、植被数据和地理数据；也即1980年的气象数据、土壤数据、植被数据和地理数据。

所述驱动数据为所述目标时间区间内，所述目标区域中每个栅格的气象数据、土壤数据、植被数据和地理数据；

其中，所述植被数据是通过基于遥感数据和地面高程数据的反演过程得到的。

在上述获取本地化数据、初始数据和驱动数据的过程中，每种数据的具体内容和获取方法如下：

所述气象数据包括：每个栅格上逐日的平均气温、最高气温、最低气温、降水量、日照时数、平均气压、风速、相对湿度、地表太阳辐射，以及气象站点的经纬度和海拔高度。具体地，所述气象数据为基于中国气象局公开的实测气象数据，根据自然邻域法和反向距离加权法，以10km×10km为基本空间单元进行空间插值得到的与栅格一一对应的气象数据。

所述土壤数据包括：每个栅格上的土壤砂粒含量、土壤粉粒含量、土壤粘粒含量、土壤容重、土壤田间持水量、土壤萎蔫含水量、土壤碳密度、土壤氮密度；所述土壤数据为基于公开的全国1：100万数字化土壤图，以预设尺寸为基本空间单元，经过空间化处理得到的与栅格一一对应的土壤数据。

所述地理数据包括：所述目标区域栅格化后每个格点的经纬度和地面高程数据，其中，所述地面高程数据为基于栅格的1km×1km数据，即，一个栅格内包括100个地面高程值；所述地理数据可基于现有的GIS工具生成。

所述碳通量数据为1989～1993年中国林业部森林普查所获得的690个点的NPP数据中选出的属于东北地区的站点的NPP数据。

所述植被数据包括：每个栅格上的森林类型、逐日的叶面积指数、以及最大叶面积指数LAI_max；其中，所述森林类型可根据公开的中国植被类型空间分布数据得到，所述逐日的叶面积指数、最大叶面积指数LAI_max为通过基于遥感数据和地面高程数据的反演过程得到的；

所述遥感数据为公开的NDVI数据，所述基于遥感数据和地面高程数据的反演过程包括下述子步骤：

A1、获取目标时间区间内，每一栅格上的森林类型和逐日的NDVI数据。

A2、基于地理数据中的地面高程数据，根据公式(1)计算每一栅格上的地形起伏度指数R，所述公式(1)为：

公式(1)中，H_max表示该栅格中的最大高程值，H_min表示该栅格中的最小高程值，ΔH表示目标区域内的最大高程值与最小高程值的差值，m表示一个与地形起伏度指数有关的经验常数。m用于平衡R值的大小，使R的取值在0.8～1.2之间浮动，以免R值过大或过小而影响FPAR的精确度。

公式(1)求得的R用于后继的基于NDVI数据的LAI反演过程，用于平衡由于地形起伏对NDVI数据准确度的影响：若R值越大，说明地形起伏越大，则该栅格内植被接收阳光辐射的面积越大，则FPAR越大；反之，若R值越小，说明地形起伏越小，则该栅格内植被接收阳光辐射的面积越接近栅格面积，FPAR相对减小。

A3、基于一栅格上的R和逐日的NDVI数据，根据公式(2)得到每一栅格上逐日的光合有效辐射吸收比例FPAR，所述公式(2)为。

公式(2)中，NDVI表示每一栅格上逐日的NDVI数据，NDVI_max表示与森林类型相关的最大NDVI值，NDVI_min表示与森林类型相关的最小NDVI值，NDVI_max和NDVI_min的值是与森林类型相关的经验值。

具体地，针对东北地区，主要森林类型为常绿针叶林、落叶针叶林、落叶阔叶林和针阔混交林，NDVI_max和NDVI_min的值与森林类型的对应关系如表1。

表1东北地区森林类型与NDVI_max、NDVI_min、LAI_v对应关系表

A4、基于逐日的FPAR，比较得到目标时间区间中1981年中最大的FPAR值FPAR_max和最小的FPAR值FPAR_min。

A5、基于FPAR_max、FPAR_min和逐日的FPAR，根据公式(3)计算LAI_max、LAI_min、逐日的LAI，所述公式(3)为：

公式(3)中，LAI_v表示与森林类型相关的一个经验值；

将FPAR_max和逐日的FPAR带入公式(3)，得到逐日的LAI；

将FPAR_max和FPAR＝FPAR_max带入公式(3)，得到LAI_min；

将FPAR_max和FPAR＝FPAR_min带入公式(3)，得到LAI_max。

S3、基于所述本地化数据对预先搭建的FORCCHN模型进行参数优化，得到本地化的FORCCHN模型。

具体地，S3包括下述子步骤：

B1、将本地化数据中的气象数据、土壤数据、植被数据和地理数据输入预先搭建的FORCCHN模型运行，基于FORCCHN模型输出的结果计算得到模拟结果数据；

其中，所述模拟结果数据包括模拟的森林净初级生产力NPP_simulated、模拟的森林净生态系统生产力NEP_simulated、模拟的森林净生态系统碳交换量NEE_simulated中的至少一种；具体的，本实施例中采用模拟的森林净初级生产力NPP_simulated作为模拟结果数据。

B2、基于本地化数据中的碳通量数据得到与模拟结果数据对应的实际结果数据；

其中，所述实际结果数据包括实际的森林净初级生产力NPP_true、净生态系统生产力NEP_true、净生态系统碳交换量NEE_true中的至少一种；具体地，本实施中使用S2中获取的1989～1993年中国林业部森林普查所获得的东北地区的站点的NPP数据NPP_true作为实际结果数据。

B3、将所述模拟结果数据与实际结果数据比较；若两者的差异在预设范围内，则将当前的FORCCHN模型作为本地化的FORCCHN模型；

若两者的差异超过预设范围，则循环执行调参操作，直至两者的差异在预设范围内，得到本地化的FORCCHN模型。

所述调参操作包括：

调整FORCCHN模型中的树木生理生态参数，将本地化数据中的气象数据、土壤数据、植被数据和地理数据输入调整树木生理生态参数后的FORCCHN模型中运行，基于输出的结果计算得到新的模拟结果数据，将所述新的模拟结果数据与实际结果数据比较。

具体地，在B3中，所述模拟结果数据与实际结果数据比较的指标以及两者差异的预设范围可根据实际需求确定，例如，所述比较的指标可为相关系数Corr、均方根误差RMSE、以及偏差Bias；又例如，可基于NPP_simulated与NPP_true的比值进行比较，如图3所示，坐标系中的每一个点的横坐标为NPP实际观测值NPP_true，纵坐标为FORCCHN模型的模拟值NPP_simulated，在NPP_true＝NPP_simulated的理想情况下，则所有点均应落在y＝x线上，即所有点的回归直线方程的斜率为1。因而在实际的参数优化过程中，可基于NPP_simulated与NPP_true计算回归直线方程，回归直线方程的斜率越接近1，表示FORCCHN模型的参数优化效果越好，因此，可将两者差异的预设范围设置为上述回归直线方程的斜率在0.8～1.2之间。在图3中，左图为基于刚搭建好的原始FORCCHN模型模拟得到的于NPP_simulated与NPP_true的比较结果，右图为本实施例中经过若干次调参操作后的NPP_simulated与NPP_true的比较结果，对比左图和右图可以看出，基于经过调参操作后的FORCCHN模型得到的回归直线方程的斜率明显向1接近，其回归直线方程也更接近于y＝x。

此外，为了进一步验证上述本地化的FORCCHN模型的性能，本实施例还基于本地化数据中1993年长白山站点逐月的实测NPP数据，以及该站点的土壤异养呼吸量R_s/month，根据公式(6)推算得到1993的实际NEE。

公式(6)为：

NEE＝-(NPP-R_s/month) (6)

基于上述本地化的FORCCHN模型模拟运算得到的模拟的NEE数据与上述计算得到实际的NEE数据如图4所示。需要说明的是，NEE和NEP都是常用的用于反映生态系统碳通量的指标，但其量纲不同，实际应用时两者间的转换方法是本领域技术人员所公知的。NEE是指陆地和大气之间的CO₂通量，通常与观测到的生态系统CO₂通量一致，也就是说，NEE与NEP的变化趋势是相关的，能够用于评价FORCCHN模型的性能。从图4可以看出，除个别月份外，本地化的FORCCHN模型模拟的结果基本上可以反映森林NEE的月变化特征。从上述分析可以看出，本地化的FORCCHN模型能够很好地模拟东北地区森林生态系统NEP的时空变化动态过程。

S4、将所述初始数据输入本地化的FORCCHN模型运行，得到平衡态的FORCCHN模型。

具体地，将所述初始数据输入本地化的FORCCHN模型运行，直至所述本地化的FORCCHN模型达到生态系统平衡态；其中，所述生态系统平衡态指基于FORCCHN模型运行并计算得到的NEP为零时的状态。

S5、将所述驱动数据输入平衡态的FORCCHN模型运行，得到在所述目标时间区间内，所述目标区域中每个栅格逐日的森林总初级生产力GPP_d、维持呼吸量R_m、生长呼吸量R_g、土壤异养呼吸量R_s。

S6、基于每个栅格所述逐日的GPP_d、R_m、R_g、R_s，计算得到目标区域在目标时间区间内的森林净生态系统生产力NEP。

具体地，基于每个栅格所述逐日的GPP_d、R_m、R_g、R_s，根据公式(4)得到每个栅格逐年的年森林净生态系统生产力NEP_grid，对所有栅格相同自然年的NEP_grid求和，得到目标区域在目标时间区间内逐年的NEP。

所述公式(4)为：

NEP_grid＝NPP-R_s (4)

公式(4)中，NPP表示每个栅格的年森林净初级生产力，NPP根据公式(5)得到，所述公式(5)为：

结果表明，东北地区森林生态系统1981～2002森林NEP变化于0.058～0.114PgC·a^-1之间，平均约为0.087Pg C·a^-1；在目标时间区间内东北大部分地区是碳汇区(NEP﹥0)，这也与公开文献中记载的东北地区森林生态系统的NEP变化趋势一致。具体地，基于本实施例模拟计算得到的东北地区森林NEP空间分布和变化趋势如图5所示。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例，或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本发明是参照本发明实施例的方法、设备(系统)和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。

应当注意的是，在权利要求中，不应将位于括号之间的任何附图标记理解成对权利要求的限制。词语“包含”不排除存在未列在权利要求中的部件或步骤。位于部件之前的词语“一”或“一个”不排除存在多个这样的部件。本发明可以借助于包括有若干不同部件的硬件以及借助于适当编程的计算机来实现。在列举了若干装置的权利要求中，这些装置中的若干个可以是通过同一个硬件来具体体现。词语第一、第二、第三等的使用，仅是为了表述方便，而不表示任何顺序。可将这些词语理解为部件名称的一部分。

此外，需要说明的是，在本说明书的描述中，术语“一个实施例”、“一些实施例”、“实施例”、“示例”、“具体示例”或“一些示例”等的描述，是指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。

尽管已描述了本发明的优选实施例，但本领域的技术人员在得知了基本创造性概念后，则可对这些实施例作出另外的变更和修改。所以，权利要求应该解释为包括优选实施例以及落入本发明范围的所有变更和修改。

显然，本领域的技术人员可以对本发明进行各种修改和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也应该包含这些修改和变型在内。

Claims (8)

1.一种基于FORCCHN模型和遥感数据的森林NEP的计算方法，其特征在于，包括：

S1、对目标区域按照预设尺寸进行栅格化处理；

S2、根据目标区域和目标时间区间，获取本地化数据、初始数据和驱动数据；

所述本地化数据为所述目标区域中，每个栅格任意时间区间内的气象数据、土壤数据、植被数据、地理数据和碳通量数据；

所述初始数据为所述目标时间区间的前一年内，所述目标区域中每个栅格的气象数据、土壤数据、植被数据和地理数据；

所述驱动数据为所述目标时间区间内，所述目标区域中每个栅格的气象数据、土壤数据、植被数据和地理数据；

其中，所述植被数据是通过基于遥感数据和地面高程数据的反演过程得到的；

S3、基于所述本地化数据对预先搭建的FORCCHN模型进行参数优化，得到本地化的FORCCHN模型；

S4、将所述初始数据输入本地化的FORCCHN模型运行，得到平衡态的FORCCHN模型；

S5、将所述驱动数据输入平衡态的FORCCHN模型进行动态模拟，得到在所述目标时间区间内，所述目标区域中每个栅格逐日的森林总初级生产力GPP_d、维持呼吸量R_m、生长呼吸量R_g、土壤异养呼吸量R_s；

S6、基于每个栅格逐日的GPP_d、R_m、R_g、R_s，计算得到目标区域在目标时间区间内的森林净生态系统生产力NEP；

在S3中，所述预先搭建的FORCCHN模型包括：

初始化过程，用于根据植被数据确定每一栅格上的树木株数、森林类型、树木初始基径、树木初始树高和树木生理生态参数；

光合作用过程，用于基于气象数据和植被数据模拟树木的光合作用过程，逐日计算每个栅格上树木的总光合，得到每个栅格逐日的森林总初级生产力GPP_d；

呼吸作用过程，用于基于气象数据和植被数据模拟树木的呼吸作用过程，逐日计算得到每个栅格上树木的维持呼吸量R_m和生长呼吸量R_g；

凋落过程，用于计算树木的凋落通量；

光合产物分配过程，用于基于光合作用过程、呼吸作用过程和凋落过程，计算光合产物向多个碳库的分配量；

土壤有机物呼吸和转移过程，用于基于气象数据、土壤数据，以及凋落过程，模拟土壤中逐日的呼吸、转移和氮的矿化过程，得到土壤逐日的呼吸排量，基于土壤逐日的呼吸排量计算得到土壤逐日的土壤异养呼吸量R_s；

年生长过程，用于基于光合产物分配过程和凋落过程，模拟树木的年生长过程，计算树木基径增量和树木树高增量，将树木初始基径与树木基径增量之和作为下一自然年的树木初始基径，将树木初始树高与树木树高增量之和作为下一自然年的树木初始树高。

2.根据权利要求1所述的方法，其特征在于，在S2中，

所述气象数据包括：每个栅格上逐日的平均气温、最高气温、最低气温、降水量、日照时数、平均气压、风速、相对湿度、地表太阳辐射，以及气象站点的经纬度和海拔高度；

所述土壤数据包括：每个栅格上的土壤砂粒含量、土壤粉粒含量、土壤粘粒含量、土壤容重、土壤田间持水量、土壤萎蔫含水量、土壤碳密度、土壤氮密度；

所述植被数据包括：每个栅格上的森林类型、逐日的叶面积指数LAI、最大叶面积指数LAI_max；

所述地理数据包括：所述目标区域栅格化后每个格点的经纬度和地面高程数据；

所述碳通量数据包括已公开的碳通量实测数据，和/或公开的模拟数据。

3.根据权利要求2所述的方法，其特征在于，所述植被数据中，所述逐日的叶面积指数、最大叶面积指数LAI_max为通过基于遥感数据和地面高程数据的反演过程得到的；

所述遥感数据为公开的NDVI数据，所述基于遥感数据和地面高程数据的反演过程包括：

A1、获取目标时间区间内，每一栅格上逐日的NDVI数据；

A2、基于地理数据中的地面高程数据，根据公式(1)计算每一栅格上的地形起伏度指数R，所述公式(1)为：

公式(1)中，H_max表示该栅格中的最大高程值，H_min表示该栅格中的最小高程值，ΔH表示目标区域内的最大高程值与最小高程值的差值，m表示一个与地形起伏度指数有关的经验常数；

A3、基于一栅格上的R和逐日的NDVI数据，根据公式(2)得到每一栅格上逐日的光合有效辐射吸收比例FPAR，所述公式(2)为：

公式(2)中，NDVI表示每一栅格上逐日的NDVI数据，NDVI_max表示与森林类型相关的最大NDVI值，NDVI_min表示与森林类型相关的最小NDVI值，NDVI_max和NDVI_min的值是与森林类型相关的一个经验值；

A4、基于逐日的FPAR，比较得到目标时间区间中第一自然年中最大的FPAR值FPAR_max和最小的FPAR值FPAR_min；

A5、基于FPAR_min和逐日的FPAR，根据公式(3)计算LAI_max、逐日的LAI，所述公式(3)为：

公式(3)中，LAI_v表示与森林类型相关的一个经验值；

将FPAR_max和逐日的FPAR带入公式(3)，得到逐日的LAI；

将FPAR_max和FPAR＝FPAR_min带入公式(3)，得到LAI_max。

4.根据权利要求2所述的方法，其特征在于，所述气象数据为基于中国气象局公开的实测气象数据，根据自然邻域法和反向距离加权法，以预设尺寸为基本空间单元进行空间插值得到的与栅格一一对应的气象数据。

5.根据权利要求2所述的方法，其特征在于，所述土壤数据为基于公开的全国1：100万数字化土壤图，以预设尺寸为基本空间单元，经过空间化处理得到的与栅格一一对应的土壤数据。

6.根据权利要求1所述的方法，其特征在于，所述S3包括：

B1、将本地化数据中的气象数据、土壤数据、植被数据和地理数据输入预先搭建的FORCCHN模型运行，基于FORCCHN模型输出的结果计算得到模拟结果数据；

其中，所述模拟结果数据包括模拟的森林净初级生产力NPP_simulated、模拟的森林净生态系统生产力NEP_simulated、模拟的森林净生态系统碳交换量NEE_simulated中的至少一种；

B2、基于本地化数据中的碳通量数据得到与模拟结果数据对应的实际结果数据；

其中，所述实际结果数据包括实际的森林净初级生产力NPP_true、净生态系统生产力NEP_true、净生态系统碳交换量NEE_true中的至少一种；

B3、将所述模拟结果数据与实际结果数据比较；若两者的差异在预设范围内，则将当前的FORCCHN模型作为本地化的FORCCHN模型；

若两者的差异超过预设范围，则循环执行调参操作，直至两者的差异在预设范围内，得到本地化的FORCCHN模型；

所述调参操作包括：

调整FORCCHN模型中的树木生理生态参数，将本地化数据中的气象数据、土壤数据、植被数据和地理数据输入调整树木生理生态参数后的FORCCHN模型中运行，基于输出的结果计算得到新的模拟结果数据，将所述新的模拟结果数据与实际结果数据比较。

7.根据权利要求1所述的方法，其特征在于，所述S4包括：

将所述初始数据输入本地化的FORCCHN模型运行，直至所述本地化的FORCCHN模型达到生态系统平衡态；其中，所述生态系统平衡态指基于FORCCHN模型运行并计算得到的NEP为零的状态。

8.根据权利要求1所述的方法，其特征在于，所述目标时间区间为一个或多个连续的自然年，

所述S6包括：基于每个栅格逐日的GPP_d、R_m、R_g、R_s，根据公式(4)得到每个栅格逐年的年森林净生态系统生产力NEP_grid，对所有栅格相同自然年的NEP_grid求和，得到目标区域在目标时间区间内逐年的NEP；

所述公式(4)为：

公式(4)中，NPP表示每个栅格的年森林净初级生产力，NPP根据公式(5)得到，所述公式(5)为：

Images