碳通量计算方法及装置

技术领域

本发明涉及遥感和大气交叉技术领域，具体而言，涉及一种碳通量计算方法及装置。

背景技术

发展可靠的观测方法和观测技术是开展湍流通量研究的前提条件。传统的通量观测的方法主要有：空气动力学法、波文比能量平衡法和涡度相关法等。

其中，涡度相关法为直接计算方法，其利用快速响应仪器直接测量湍流脉动量得到通量，具有理论最为完备可靠、观测精度相对较高等优势，因而在多数研究中，都是将其作为一种“准”标准方法，但其缺点是观测仪器价格昂贵，操作复杂，不易推广普及。同时其测量的准确性受一些环境条件的限制，导致能量不能闭合，且观测高度通常只有十几米至几百米左右。空气动力学法和波文比能量平衡法为间接计算方法，其利用近地面层的风速、温度、湿度梯度观测资料及净辐射和土壤热通量的观测资料，间接计算湍流通量。由于近地面层的气象要素梯度和净辐射观测资料相对容易获得，而且观测设备成本较低，因此间接计算方法被广泛应用。

现有湍流通量的间接计算方法，存在依赖于大气稳定度以及局限于近地面层的问题。

发明内容

本发明解决的是现有湍流通量的间接计算方法，存在依赖于大气稳定度以及局限于近地面层的问题。

为解决上述问题，本发明提供一种碳通量计算方法，所述方法包括：获取气象站点的气象要素信息，及根据所述气象要素信息计算得到风廓线和二氧化碳浓度廓线，所述气象要素信息包括平均风速、二氧化碳平均浓度；根据所述风廓线及全边界层通量模型，通过多点非迭代算法计算摩擦速度；所述全边界层通量模型由近地层湍流特征尺度及边界层中上部的湍流特征尺度共同确定，所述边界层中上部的湍流特征尺度与所述边界层的厚度相关；根据所述摩擦速度计算得到湍流扩散系数；根据所述二氧化碳浓度廓线计算得到二氧化碳浓度梯度；将不同高度的湍流扩散系数乘以所述二氧化碳浓度梯度得到二氧化碳通量。

可选地，所述全边界层通量模型的表达式如下：

其中，u

可选地，所述根据所述风廓线及全边界层通量模型，通过多点非迭代算法计算摩擦速度，包括：根据预先拟合的莫宁长度与稳定度函数、测量高度对应平均风速的变化关系式，确定多个目标高度对应的平均风速对应的目标莫宁长度；根据所述目标莫宁长度确定稳定度函数，以及将所述稳定度函数输入所述全边界层通量模型的表达式，计算得到摩擦速度。

可选地，若已知z

其中，

若满足z

拟合莫宁长度L随

可选地，所述湍流扩散系数的计算公式如下：

其中，K

可选地，所述二氧化碳通量的计算公式如下：

其中，F

本发明提供一种碳通量计算装置，所述装置包括：廓线计算模块，用于获取气象站点的气象要素信息，及根据所述气象要素信息计算得到风廓线和二氧化碳浓度廓线，所述气象要素信息包括平均风速、二氧化碳平均浓度；摩擦速度计算模块，用于根据所述风廓线及全边界层通量模型，通过多点非迭代算法计算摩擦速度；所述全边界层通量模型由近地层湍流特征尺度及边界层中上部的湍流特征尺度共同确定，所述边界层中上部的湍流特征尺度与所述边界层的厚度相关；湍流扩散系数计算模块，用于根据所述摩擦速度计算得到湍流扩散系数；浓度梯度计算模块，用于根据所述二氧化碳浓度廓线计算得到二氧化碳浓度梯度；二氧化碳通量计算模块，用于将不同高度的湍流扩散系数乘以所述二氧化碳浓度梯度得到二氧化碳通量。

可选地，所述全边界层通量模型的表达式如下：

其中，u

可选地，所述摩擦速度计算模块，具体用于：

根据预先拟合的莫宁长度与稳定度函数、测量高度对应平均风速的变化关系式，确定多个目标高度对应的平均风速对应的目标莫宁长度；

根据所述目标莫宁长度确定稳定度函数，以及将所述稳定度函数输入所述全边界层通量模型的表达式，计算得到摩擦速度。

可选地，若已知z

其中，

若满足z

拟合莫宁长度L随

将

根据所述目标莫宁长度及所述全边界层通量模型的表达式，计算得到摩擦速度u

本发明提供的基于全边界层廓线数据测量的碳通量计算方法及装置，通过构建全边界层通量模型拓展了近地层通量模型的使用范围，采用多点法非迭代算法计算莫宁长度和摩擦速度，减小了计算量，而且可以避免事前对大气稳定度进行判断，以及可以直接计算湍流扩散系数。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。

图1为本发明的一个实施例中一种基于全边界层廓线数据测量的碳通量计算方法的示意性流程图；

图2为本发明的一个实施例中黑河流域大满站(2019年)的摩擦速度的对比示意图；

图3为本发明的一个实施例中黑河流域四道桥站(2019年)的摩擦速度的对比示意图；

图4为本发明的一个实施例中黑河流域阿柔站(2019年)的摩擦速度的对比示意图；

图5为本发明的一个实施例中黑河流域阿柔站(2019年)的每月平均摩擦速度的对比示意图；

图6为本发明的一个实施例中北京325m观测站点(2015年春季)的二氧化碳通量的计算值与测量值的对比示意图；

图7为本发明的一个实施例中北京325m观测站点(2015年夏季)的二氧化碳通量的计算值与测量值的对比示意图；

图8为本发明的一个实施例中北京325m观测站点(2015年秋季)的二氧化碳通量的计算值与测量值的对比示意图；

图9为本发明的一个实施例中北京325m观测站点(2015年冬季)的二氧化碳通量的计算值与测量值的对比示意图；

图10为本发明的一个实施例中北京325m观测站点(2015年全年)的二氧化碳通量的计算值与测量值的对比示意图；

图11为本发明的一个实施例中一种基于全边界层廓线数据测量的碳通量计算装置的结构示意图。

具体实施方式

为使本发明的上述目的、特征和优点能够更为明显易懂，下面结合附图对本发明的具体实施例做详细的说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

图1是本发明的一个实施例中一种基于全边界层廓线数据测量的碳通量计算方法的示意性流程图，该方法包括以下步骤：

S102，获取气象站点的气象要素信息，及根据气象要素信息计算得到风廓线和二氧化碳浓度廓线。

根据实际需求选取气象站点的各个气象要素，例如多个不同垂直高度对应的风速、二氧化碳浓度等，然后计算各气象要素的时均量，例如平均风速、二氧化碳平均浓度等。在得到上述时均量后计算得到风廓线和二氧化碳浓度廓线。反演风廓线和二氧化碳浓度廓线的方法采用现有的反演方法，在本实施例中对此不作限定。

S104，根据上述风廓线及全边界层通量模型，通过多点非迭代算法计算摩擦速度。

其中，全边界层通量模型由近地层湍流特征尺度及边界层中上部的湍流特征尺度共同确定，边界层中上部的湍流特征尺度与边界层的厚度相关。

可选地，全边界层通量模型的表达式如下：

其中，u

一般而言，近地层(100m以下)通量模型的表达式为：

其中l＝z/φ

可选地，首先，根据预先拟合的莫宁长度与稳定度函数、测量高度对应平均风速的变化关系式，确定多个目标高度对应的平均风速对应的目标莫宁长度；

然后，根据上述目标莫宁长度确定稳定度函数，以及将该稳定度函数输入全边界层通量模型的表达式，计算得到摩擦速度。以下以近地层通量模型为例，阐明多点非迭代算法的计算原理。

例如，近地层通量模型为

其中，

若满足z

拟合莫宁长度L随

将

S106，根据上述摩擦速度计算得到湍流扩散系数。

湍流扩散系数的计算公式如下：

其中，K

S108，根据上述二氧化碳浓度廓线计算得到二氧化碳浓度梯度。该浓度梯度表示为

S110，将不同高度的湍流扩散系数乘以二氧化碳浓度梯度得到二氧化碳通量。

可选地，二氧化碳通量的计算公式如下：

其中，F

本实施例提供的基于全边界层廓线数据测量的碳通量计算方法，通过构建全边界层通量模型拓展了近地层通量模型的使用范围，采用多点法非迭代算法计算莫宁长度和摩擦速度，减小了计算量，而且可以避免事前对大气稳定度进行判断，以及可以直接计算湍流扩散系数。

为验证上述基于全边界层廓线数据测量的碳通量计算方法的有效性，还提供了以下验证数据。

其中，边界层通量模型验证主要包含湍流通量和二氧化碳通量的综合验证。湍流通量主要受摩擦速度影响，其决定着湍流扩散系数，进而影响二氧化碳通量的计算结果，因此首先需要对摩擦速度进行验证，通过黑河流域气象站点各高度气象要素数据(主要是水平风速)进行统计平均，通过模型计算摩擦速度并与涡动相关法(Eddy Correlation，EC)直接测量的摩擦速度进行对比验证。

黑河流域地表过程综合观测网目前业务运行站点为11个(3个超级站和8个普通站)，覆盖了黑河流域主要下垫面类型。发布的数据包含地面观测数据集，包括感热通量、潜热通量和碳通量等通量观测数据集；风温湿压、降水量、辐射、土壤温湿度廓线、土壤热通量、光合有效辐射和地表辐射温度等水文气象观测数据集。其中超级站主要包含阿柔、大满和四道桥这三个站点。三大超级站的基本情况如下：

黑河流域上游阿柔超级站位于青海省祁连县阿柔乡草达坂村，下垫面是亚高山山地草甸，经纬度100.4643E，38.0473N，海拔3033m；中游大满超级站位于甘肃省张掖市大满灌区内，下垫面玉米，经纬度100.37223E，38.85551N，海拔1556.06m；下游四道桥超级站位于内蒙古额济纳旗四道桥，下垫面柽柳，经纬度101.1374E，42.0012N，海拔873m。

图2示出了黑河流域大满站(2019年)的摩擦速度的对比示意图，对比了本实施例提供的廓线法与现有的涡动相关法的对比结果，其中横坐标为EC法计算得到的平均摩擦速度u

总体上看黑河流域模型大满站点、四道桥站点、阿柔站点的计算值与实测值吻合较好。阿柔站点每月平均摩擦速度对比廓线法与涡动相关法，发现两者在夏季吻合较好，冬春季有所差别。

二氧化碳通量模型验证采用中科院大气物理所325m气象塔的实测数据。中国科学院大气物理研究所气象塔建成1979年8月，气象塔位于(39°58′N，116°22′E)，海拔高度49m，塔高325m，塔上设有15层观测平台，并装有风速仪、风向仪和温度仪。为减小塔体阻挡气流对测风仪造成的影响，分别在西北和东南两个盛行风向各装一台测风仪。其主要功能为首都空气污染以及大气边界层、大气湍流扩散等提供研究服务。目前气象塔半径2km内地势平坦，地表覆盖以建筑为主，塔的北面和南面为密集居民区，覆盖率为0.65；其次为植被和道路。距离气象塔西部约0.3km处为小型公园，由树木、草坪和河流组成。塔的东面为京藏高速，塔的北面为北辰路，每天有大量车辆行驶，覆盖率分别为0.21、0.14。

图6示出了北京325m观测站点(2015年春季)的二氧化碳通量的计算值与测量值的对比示意图，横坐标为高度，纵坐标为二氧化碳通量；图7示出了北京325m观测站点(2015年夏季)的二氧化碳通量的计算值与测量值的对比示意图；图8示出了北京325m观测站点(2015年秋季)的二氧化碳通量的计算值与测量值的对比示意图；图9示出了北京325m观测站点(2015年冬季)的二氧化碳通量的计算值与测量值的对比示意图；图10示出了北京325m观测站点(2015年全年)的二氧化碳通量的计算值与测量值的对比示意图。

总体上看春季、秋季以及全年的统计平均与模型吻合较好，尤其是二氧化碳通量的极大值点对应的横纵坐标吻合较好。夏季和冬季稍微有所偏差。

图11是本发明的一个实施例中一种基于全边界层廓线数据测量的碳通量计算装置的结构示意图，该装置包括：

廓线计算模块110，用于获取气象站点的气象要素信息，及根据气象要素信息计算得到风廓线和二氧化碳浓度廓线，气象要素信息包括平均风速、二氧化碳平均浓度；

摩擦速度计算模块112，用于根据风廓线及全边界层通量模型，通过多点非迭代算法计算摩擦速度；全边界层通量模型由近地层湍流特征尺度及边界层中上部的湍流特征尺度共同确定，边界层中上部的湍流特征尺度与边界层的厚度相关；

湍流扩散系数计算模块114，用于根据摩擦速度计算得到湍流扩散系数；

浓度梯度计算模块116，用于根据二氧化碳浓度廓线计算得到二氧化碳浓度梯度；

二氧化碳通量计算模块118，用于将不同高度的湍流扩散系数乘以二氧化碳浓度梯度得到二氧化碳通量。

本实施例提供的基于全边界层廓线数据测量的碳通量计算装置，通过构建全边界层通量模型拓展了近地层通量模型的使用范围，采用多点法非迭代算法计算莫宁长度和摩擦速度，减小了计算量，而且可以避免事前对大气稳定度进行判断，以及可以直接计算湍流扩散系数。

可选地，作为一个实施例，全边界层通量模型的表达式如下：

其中，u

可选地，作为一个实施例，所述摩擦速度计算模块112，具体用于：根据预先拟合的莫宁长度与稳定度函数、测量高度对应平均风速的变化关系式，确定多个目标高度对应的平均风速对应的目标莫宁长度；根据目标莫宁长度确定稳定度函数，以及将稳定度函数输入全边界层通量模型的表达式，计算得到摩擦速度。

可选地，作为一个实施例，若已知z

其中，

若满足z

拟合莫宁长度L随

将

根据目标莫宁长度及全边界层通量模型的表达式，计算得到摩擦速度u

可选地，作为一个实施例，所述湍流扩散系数的计算公式如下：

其中，K

可选地，作为一个实施例，所述二氧化碳通量的计算公式如下：

其中，F

上述实施例提供的基于全边界层廓线数据测量的碳通量计算装置能够实现上述基于全边界层廓线数据测量的碳通量计算方法的实施例中的各个过程，为避免重复，这里不再赘述。

本发明实施例还提供一种计算机可读存储介质，计算机可读存储介质上存储有计算机程序，该计算机程序被处理器执行时实现上述基于全边界层廓线数据测量的碳通量计算方法实施例的各个过程，且能达到相同的技术效果，为避免重复，这里不再赘述。其中，所述的计算机可读存储介质，如只读存储器(Read-Only Memory，简称ROM)、随机存取存储器(Random Access Memory，简称RAM)、磁碟或者光盘等。

当然，本领域技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令控制装置来完成，所述的程序可存储于一计算机可读取的存储介质中，所述程序在执行时可包括如上述各方法实施例的流程，其中所述的存储介质可为存储器、磁盘、光盘等。

在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。