植被总初级生产力估算方法、装置、电子设备及存储介质与流程

1.本技术涉及遥感图像处理技术领域，尤其是涉及植被总初级生产力估算方法、装置、电子设备及存储介质。


背景技术：

2.植被总初级生产力（gross primary productivity，gpp）是估算地球支撑能力和评价陆地生态系统可持续发展的重要因子，在全球变化以及碳循环中扮演着重要角色，对植被gpp 的准确评估显得尤为重要。gpp是指生态系统中绿色植物单位时间通过光合作用途径固定有机碳的总量，决定了进入陆地生态系统的初始能量与物质总量。国际上通用的估算gpp方法主要有通量站连续观测、陆地生态过程模型估测等方法。结合遥感数据的gpp估测模型实现了空间连续、不破坏植被的估测植被gpp。遥感估测gpp模型主要分为3类：经验型植被指数模型、植被生态过程模型和机器学习模型。
3.其中遥感光能利用率模型利用植被吸收的光合有效辐射（absorbed photosynthetic active radiation，apar）和光能利用率(light use efficiency, lue)估算陆地净初级生产力，并考虑养分、水分和温度等环境胁迫因素对植被光合作用的影响，是应用较多的结合植被生态过程生产gpp的模型。光能利用率模型输入数据主要来自于遥感数据，降低或避免了对实地观测数据的依赖，被广泛的应用于区域和全球尺度的 gpp 估算研究。当前在光能利用率模型基础上，对模型中的一些参数进行全球参数化，结合长时间序列数据，已经能够生产全球gpp遥感产品。但是，现有gpp估算模型和算法中存在一些不足，现有产品中全球不同区域产品精度不一，难以满足实际应用需求。


技术实现要素：

4.有鉴于此，本技术提供了植被总初级生产力估算方法、装置、电子设备及存储介质，能够解决植被总初级生产力估算精度不高的技术问题。
5.第一方面，本技术实施例提供了一种植被总初级生产力估算方法，包括：获取监测时间段的目标区域的遥感影像和气象观测数据；获取遥感影像的各像元的空间位置信息，对气象观测数据进行栅格化处理，得到遥感影像的各像元的气象观测数据；基于遥感影像的各像元的遥感数据和气象观测数据，计算遥感影像每个像元的特征量，所述特征量包括：光合有效辐射比例、水分限制因子和温度限制因子；将遥感影像每个像元的特征量输入预先训练完成的植被总初级生产力估算模型，输出每个像元的植被总初级生产力估计值。
6.进一步的，所述光合有效辐射比例的计算步骤包括：获取遥感影像各像元的增强型植被指数evi，计算第一光合有效辐射比例fpar
evi
：其中，β为参数；
获取遥感影像各像元的叶面积指数lai，计算第二光合有效辐射比例fpar
lai
:其中，k为辐射消光系数。
7.进一步的，所述水分限制因子的计算步骤包括：基于遥感影像各像元的swir 波段信息和 nir 波段信息，获取各像元的陆地表面水分指数lswi，计算第一水分限制因子w
lswi
：其中，lswi
max
是各像元内植被生长季最大 lswi值；获取遥感影像各像元的饱和水汽压差
푉푃퐷
，计算第二水分限制因子w
vpd
：其中，vpd
max
和vpd
min
为饱和水汽压差的最大值和最小值，根据植被类型确定；当
푉푃퐷 》 vpd
max
时，w
vpd
为 0，当
푉푃퐷 《vpd
min
时，w
vpd
为 1。
8.进一步的，所述温度限制因子的计算步骤包括：计算第一温度限制因子t
tem
：当t
min 《 t 《 t
max
时：否则，t
tem
=0；式中，t为监测时间段的平均温度，t
opt
为植被光合作用速率最高时对应的温度，采用每个像元增强型植被指数evi最大时对应的温度；t
min
和t
max
分别是植被进行光合作用时最低和最高温度，t
min
=0℃，t
max
=t
opt
( t
opt
－t
min
)2；采用casa算法计算第二温度限制因子ts：：：其中，t
s 1
为降低植被生产力的因子，t
s 2
表示植物光能利用率逐渐变小的趋势因子。
9.进一步的，所述植被总初级生产力估算模型采用bp神经网络，所述植被总初级生产力估算模型的训练过程包括：获取训练集；其中，训练集中包括多个历史遥感影像样本和对应的气象观测数据样本和通量塔数据；获取各历史遥感影像样本的各像元的空间位置信息，对气象观测数据进行栅格化处理，得到各历史遥感影像样本的各像元的气象观测数据；对通量塔数据进行栅格化处理，
得到各历史遥感影像样本的各像元的通量塔数据；基于各历史遥感影像样本的各像元的遥感数据和气象观测数据，计算每个特征量样本，包括：第一光合有效辐射比例、第二光合有效辐射比例、第一水分限制因子、第二水分限制因子、第一温度限制因子和第二温度限制因子；将每个特征量样本输入bp神经网络，得到每个特征量样本对应的植被总初级生产力预测结果；基于每个特征量样本对应的植被总初级生产力预测结果和对应的通量塔数据，确定损失函数值；基于损失函数值，更新bp神经网络的权重参数。
10.第二方面，本技术实施例提供了一种植被总初级生产力估算装置，包括：获取单元，用于获取监测时间段的目标区域的遥感影像和气象观测数据；预处理单元，用于获取遥感影像的各像元的空间位置信息，对气象观测数据进行栅格化处理，得到遥感影像的各像元的气象观测数据；特征量计算单元，用于基于遥感影像的各像元的遥感数据和气象观测数据，计算遥感影像每个像元的特征量，所述特征量包括：光合有效辐射比例、水分限制因子和温度限制因子；植被总初级生产力估算单元，用于将遥感影像每个像元的特征量输入预先训练完成的植被总初级生产力估算模型，输出每个像元的植被总初级生产力估计值。
11.第三方面，本技术实施例提供了一种电子设备，包括：存储器、处理器和存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现本技术实施例的植被总初级生产力估算方法。
12.第四方面，本技术实施例提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机指令，所述计算机指令被处理器执行时实现本技术实施例的植被总初级生产力估算方法。
13.本技术对于不同植被在不同地区均有较好的适应性，能充分反应植被类型和区域的差异，具有更强的空间异质性和更优的拟合效果，提高了植被总初级生产力估算精度。
附图说明
14.为了更清楚地说明本技术具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本技术的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
15.图1为本技术实施例提供的植被总初级生产力估算方法的流程图；图2为本技术实施例提供的植被总初级生产力估算模型的训练过程示意图；图3为本技术实施例提供的植被总初级生产力估算装置的功能结构图；图4为本技术实施例提供的电子设备的硬件结构示意图。
具体实施方式
16.为使本技术实施例的目的、技术方案和优点更加清楚，下面将结合本技术实施例
中的附图，对本技术实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本技术一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本技术实施例的组件可以以各种不同的配置来布置和设计。
17.因此，以下对在附图中提供的本技术的实施例的详细描述并非旨在限制要求保护的本技术的范围，而是仅仅表示本技术的选定实施例。基于本技术中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本技术保护的范围。
18.首先对本技术实施例涉及的技术用语进行简单介绍。
19.植被总初级生产力gpp（gross primary productivity）：表示植被在单位时间内通过光合作用固定的光合产物量或有机碳总量，又称总第一性生产力或总生态系统生产力。
20.遥感光能利用率模型：假设环境因素简化情况下，植被生产力只与温度，水分，光照以及该地区本身光能利用率有关来估算植被生产力的模型。
21.光合有效辐射吸收比例 fpar（the fraction of absorbed photosynthetically active radiation）: 一般定义为植被对波长在400 至700nm间太阳辐射能量的吸收比率，是表征植被光合作用水平和生长状态的重要参量，也是联合国全球气候观测系统认定的反映全球气候变化的关键参量之一。
22.vpm（vegetation photosynthesis model）模型: vpm模型是利用温度等遥感数据以及涡度观测碳通量数据，并且考虑了植被叶绿素吸收的光合有效辐射来估计总初级生产力的一种光能利用率模型。
23.casa（carnegie-ames-stanford model）模型: casa模型是一个充分考虑环境条件和植被本身特征的光能利用率模型，已被广泛应用于陆地生态系统 npp（net primary productivity）的估算。
24.cflux（carbon flux）碳循环模型：利用通量站和遥感卫星来估计总初级生产力的一种光能利用率模型。
25.mod17(moderate resolution imaging spectroradiometer-gpp)模型：主要利用土地利用和环境因素包括温度、水汽压、光照三者之间的关系来估计总初级生产力的一种光能利用率模型。
26.在介绍了本技术涉及的技术用语后，接下来，对本技术实施例的设计思想进行简单介绍。
27.现有gpp估算模型和算法中普遍存在精度不高的问题，而且产品中全球不同区域的精度不一，难以满足实际应用需求。
28.为解决上述技术问题，本技术利用历史遥感数据和历史气象站数据，计算得到历史特征量：光合有效辐射吸收比例（fpar）、水分对光合作用的影响（水分限制因子）和温度对光合作用的影响(温度限制因子)；利用历史特征参量和历史通量塔数据对植被总初级生产力估算模型进行训练；利用训练完成的植被总初级生产力估算模型对当前特征量进行处理，得到植被总初级生产力估算值。本技术的植被总初级生产力估算模型依靠机器学习的支撑，相较于传统光能利用率模型，对于不同植被在不同地区均有较好的适应性，更能反应植被类型和区域的差异，具有更强的空间异质性和更优的拟合效果，由此提高了植被总初
级生产力估算精度。同时，使利用机器学习算法估算gpp的过程更为清晰，模型相较于传统“黑箱”模式的机器学习更具有物理过程意义。
29.本技术的gpp估算方法相较于传统“黑箱”模式的机器学习，引入具有物理意义的计算公式，模型训练更具有针对性，可提高训练效率和精确度。
30.在介绍了本技术实施例的应用场景和设计思想之后，下面对本技术实施例提供的技术方案进行说明。
31.如图1所示，本技术提供了一种植被总初级生产力估算方法，包括：步骤101：获取监测时间段的目标区域的遥感影像和气象观测数据；其中，监测时间段根据输入数据的时间分辨率确定，一般为天，月或年。
32.步骤102：获取遥感影像的各像元的空间位置信息，对气象观测数据进行栅格化处理，得到遥感影像的各像元的气象观测数据；其中，由于气象观测数据对应的空间位置和遥感影像的各像元的空间位置不一致，因此需要对气象观测数据进行插值处理，得到遥感影像的各像元的空间位置处气象观测数据。
33.步骤103：基于遥感影像的各像元的遥感数据和气象观测数据，计算遥感影像每个像元的特征量，所述特征量包括：光合有效辐射比例、水分限制因子和温度限制因子；首先计算光合有效辐射比例：在vpm模型中，植被光合作用所吸收的光合有效辐射比例 fpar可近似地表示为增强型植被指数 evi 的线性函数，获取遥感影像各像元的增强型植被指数evi，计算第一光合有效辐射比例fpar
evi
：其中，一般情况下，参数β的取值为1；在cflux模型和mod17模型中，也可以通过叶面积指数（leaf area index，lai）来计算得到光合有效辐射比例。
34.获取遥感影像各像元的叶面积指数lai，计算第二光合有效辐射比例fpar
lai
:其中，k为辐射消光系数。
35.然后计算水分限制因子：由于swir 波段信息对地表和植被水分比较敏感，结合 nir 波段得到 陆地表面水分指数lswi，陆地表面水分指数lswi对叶片水分含量比较敏感，因此第一水分限制因子w
lswi
的计算公式为：其中lswi
max
是单个像元内植被生长季最大 lswi 值。
36.在mod17算法中，水分限制因子还可以通过饱和水汽压差（vpd, vaport pressure deficit）得到，通过下述公式计算第二水分限制因子w
vpd
：
其中，vpd
max
和vpd
min
为饱和水汽压差的最大值和最小值，与植被类型相关，可通过 mod17模型 查找表中获取到。当
푉푃퐷 》 vpd
max
时，w
vpd
为 0，当
푉푃퐷 《vpd
min
时，w
vpd
为 1。
37.最后计算温度限制因子：采用陆地生态系统模型（terrestrial ecosystem model,tem）中温度限制因子的算法，计算第一温度限制因子t
tem
：当t
min 《 t 《 t
max
时：当t ≤t
min
或t ≥t
max
时：t
tem = 0式中，t为监测时间段的平均温度（℃）；t
opt
为植被光合作用速率最高时对应的温度，采用生长季每个像元 evi 最大时对应的温度； t
min
和t
max
分别是植被进行光合作用时最低和最高温度，当温度低于最低温度或高于最高温度时，t
tem 设为0。t
min
设为 0℃，t
max
计算公式为t
opt
( t
opt
－t
min
)2。
38.casa模型中温度限制因子的算法为：casa模型中温度限制因子的算法为：casa模型中温度限制因子的算法为：t
s 1
反映在低温和高温时植物内的生化作用对光合作用的限制而降低植被生产力，当均温等于或低于-10℃时，t
s 1
取值为0，此时不发生光合作用；t
s 2
表示环境从最适温度向高温和低温变化时植物光能利用率逐渐变小的趋势，例如，当某个月平均温度t比适宜温度t
opt
高10℃或低于13℃，将月平均温度t取为最适宜温度t
opt
，带入上述公式计算出t
s 2
值，将该值的一半作为月平均温度t对应的t
s 2
的值。
39.步骤104：将遥感影像每个像元的特征量输入预先训练完成的植被总初级生产力估算模型，输出每个像元的植被总初级生产力估计值；在本实施例中，植被总初级生产力估算模型采用bp神经网络，所述植被总初级生产力估算模型的训练过程包括：步骤4a：获取样本集；其中，样本集中包括多个历史遥感影像样本和对应的气象观测数据样本和通量塔数据；将样本数据集的60%作为训练集，20%作为验证集，20%作为测试集。
40.步骤4b：获取各历史遥感影像样本的各像元的空间位置信息，对气象观测数据进行栅格化处理，得到各历史遥感影像样本的各像元的气象观测数据；对通量塔数据进行栅格化处理，得到各历史遥感影像样本的各像元的通量塔数据；步骤4c：基于各历史遥感影像样本的各像元的遥感数据和气象观测数据，计算每
个特征量样本，包括：第一光合有效辐射比例、第二光合有效辐射比例、第一水分限制因子、第二水分限制因子、第一温度限制因子和第二温度限制因子；步骤4d：将每个特征量样本输入bp神经网络，得到每个特征量样本对应的植被总初级生产力预测结果；如图2所示，首先设置bp神经网络不同层权重初始值和偏置初始值（i表示层数）；从训练集中选择n个样本；正向传播：对于不同层输出，根据以下公式：从训练集中选择n个样本；正向传播：对于不同层输出，根据以下公式：从训练集中选择n个样本；正向传播：对于不同层输出，根据以下公式：从训练集中选择n个样本；正向传播：对于不同层输出，根据以下公式：从训练集中选择n个样本；正向传播：对于不同层输出，根据以下公式：从训练集中选择n个样本；正向传播：对于不同层输出，根据以下公式：从训练集中选择n个样本；正向传播：对于不同层输出，根据以下公式：从训练集中选择n个样本；正向传播：对于不同层输出，根据以下公式：从训练集中选择n个样本；正向传播：对于不同层输出，根据以下公式：
……
其中表示的是向量，代表的是权重；b为神经元的偏置；激活函数f使用relu型函数；用表示相应神经元的输出；σ为gpp预测结果。
41.步骤4e：基于每个特征量样本对应的植被总初级生产力预测结果和通量塔数据，确定损失函数值；其中，损失函数l使用均方误差（mean-squared loss，mse）来表示：其中，n是样本量，表示gpp真实值，即通量塔数据，y
pred
是gpp的预测值。训练网络的目的就是使mse最小化。
42.由前述流程可知，当不同层权重和偏置b发生改变时，损失函数l也随之变化，其变化大小可以使用和来表示。
43.由于当n = 1 时，则，对于，为了简便演示计算过程，这里可变化层取为0，此时有：，为了简便演示计算过程，这里可变化层取为0，此时有：，为了简便演示计算过程，这里可变化层取为0，此时有：而只会影响到和，所以，所以对于其他权重和偏置b偏导的计算方式与上述过程相似，不再赘述。
44.步骤4f：基于损失函数值，更新bp神经网络的权重参数；为了实现损失函数的最小化，使用随机梯度下降的方法来优化网络的权重和偏置。对于权重和偏置置。对于权重和偏置其中η为一个常数，用于调整训练的速度。当为正数时，和会减小，l就会下降。当对网络中的每个权重和截距项都这样进行优化，损失就会不断下降，性能会不断上升。
45.l随着不断更新权重和偏置会不断减小，当l处在较低值，且后续训练l趋于平缓，
disks，raid）子系统、磁带驱动器以及数据备份存储子系统等。
54.需要说明的是，图4所示的电子设备300仅仅是一个示例，不应对本技术实施例的功能和使用范围带来任何限制。
55.本技术实施例还提供了一种计算机可读存储介质，该计算机可读存储介质存储有计算机指令，该计算机指令被处理器执行时实现本技术实施例提供的植被总初级生产力估算方法。
56.应当注意，尽管在上文详细描述中提及了装置的若干单元或子单元，但是这种划分仅仅是示例性的并非强制性的。实际上，根据本技术的实施方式，上文描述的两个或更多单元的特征和功能可以在一个单元中具体化。反之，上文描述的一个单元的特征和功能可以进一步划分为由多个单元来具体化。
57.此外，尽管在附图中以特定顺序描述了本技术方法的操作，但是，这并非要求或者暗示必须按照该特定顺序来执行这些操作，或是必须执行全部所示的操作才能实现期望的结果。附加地或备选地，可以省略某些步骤，将多个步骤合并为一个步骤执行，和/或将一个步骤分解为多个步骤执行。
58.尽管已描述了本技术的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本技术范围的所有变更和修改。
59.最后应说明的是：以上各实施例仅用以说明本技术的技术方案，而非对其限制；尽管参照前述各实施例对本技术进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本技术各实施例技术方案的范围。