CN115907178A - 一种净生态***co2交换量的预测方法 - Google Patents

Abstract

本发明公开了一种净生态***CO₂交换量的预测方法，针对当前NEE的环境影响因子众多且多为高度非线性关系而无法精确预测的问题，本发明将多种环境驱动因子考虑在内，采用最大互信息系数MIC筛选影响NEE的主控因子，构建以主控因子为自变量以NEE为目标变量的数据库，采用支持向量回归机理论和人工蜂群优化算法，建立了NEE最优化预测模型，有效地提高了NEE的预测精度。

Description

一种净生态***CO2交换量的预测方法

技术领域

本发明涉及环境监测领域，尤其涉及一种净生态***CO₂交换量的预测方法。

背景技术

CO₂净生态***交换量NEE收到广泛关注，NEE的精确预测极为重要，但是NEE的影响因素众多且为高度非线性关系，目前已有的方法无法对其进行精确的预测本发明将多种环境驱动因子考虑在内，本发明将多种环境驱动因子考虑在内，采用最大互信息系数MIC筛选影响NEE的主控因子，构建以主控因子为自变量以NEE为目标变量的数据库，采用支持向量回归机理论和人工蜂群优化算法，建立了NEE最优化预测模型，有效地提高了NEE的预测精度。

发明内容

为了针对现有NEE环境影响因子众多且为高度非线性关系而无法精确预测的问题，本发明提供一种净生态***CO₂交换量的预测方法，方法包括以下步骤：

S1、获取涡度相关***连续性观测数据；所述涡度相关***连续性观测数据包括：净生态***CO₂交换量和环境驱动变量数据；

S2、对涡度相关***连续性观测数据进行数据预处理，得到有效的净生态***CO₂交换量和环境驱动变量数据；

S3、利用最大互信息系数法MIC析有效的净生态***CO₂交换量与环境驱动变量数据的相关性，从环境驱动变量数据中筛选与净生态***CO₂交换量相关性较大的前N个主控因子；以N个主控因子为输入变量，对应的净生态***CO₂交换量为输出变量，构建有效数据库；

S4、对有效数据库中的每一项数据作归一化处理，得到归一化后的数据；将将归一化后的数据分为训练集和测试集；

S5、建立净生态***CO₂交换量支持向量机预测模型，利用训练集训练所述支持向量机模型，得到训练完成的模型；

S6、利用人工蜂群优化算法优化所述训练完成的模型，得到最优预测模型；

S7、利用测试集测试最优预测模型，并对最优预测模型进行指标评价，得到最终模型；

S8、利用最终模型进行净生态***CO₂交换量预测。

本发明提供的有益效果是：具体考虑了NEE与其各主控环境驱动变量之间的关系，采用筛选的主控因子作为NEE模型的输入变量，有效改善了模型的预测性能。

附图说明

图1是本发明方法的流程图；

图2是训练数据库预测结果；

图3是验证数据库预测结果。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地描述。

请参考图1，图1是本发明方法流程图。本发明提供了一种净生态***CO₂交换量的预测方法，包括以下步骤：

S1、获取涡度相关***连续性观测数据；所述涡度相关***连续性观测数据包括：净生态***CO₂交换量和环境驱动变量数据；

具体的说，本发明通过神农架大九湖亚高山泥炭湿地通量观测塔的开路涡度相关***获取净生态***CO₂交换量(NEE)数据及环境驱动变量数据(气象因子数据)，气象因子主要包括空气温度(Ta)、相对湿度(RH)、光量子通量密度(PPFD)、净辐射(RN)、向下短波辐射(SWIN)、三层土壤温度(Ts1、Ts2、Ts3)、三层土壤含水量(SWC₁、SWC₂、SWC₃)、饱和水汽压差(VPD)等12个变量；

S2、对涡度相关***连续性观测数据进行数据预处理，得到有效的净生态***CO₂交换量和环境驱动变量数据；

步骤S2中数据预处理具体为：根据预设阈值区间，对处于阈值区间范围外数据异常值进行剔除。作为一种实施例而言：本发明30min的NEE由涡动相关通量计算软件EddyPro6.1对10Hz的湍流原始资料进行处理得到，其主要步骤包括：野点剔除、时间延迟校正、坐标旋转、超声虚温修正、空气密度效应修正和频率响应修正。对计算得到的30min通量资料，进一步根据质量等级标志剔除低质量的数据，并剔除由于红外气体分析仪镜面污染所造成的不可信数据。进一步对NEE资料进行夜间摩擦风速修正，即剔除夜间摩擦风速小0.2m/s的NEE数据。所形成的无缺失的数据将作为建模的有效数据集；

总体说，本发明步骤S2剔除异常值时，共进行三个方面的处理，(1)剔除夜间临界摩擦风速取0.2m/s的数据；(2)剔除降水量大于0的数据；(3)剔除-22.72μmol·s^-1·m^-2≤NEE≤22.72μmol·s^-1·m^-2以外的数据；(4)剔除NEE与平均值相差超过4倍标准差的数据。去除后获得的有效数据率为36.32％。高于国际通量网(FLUXNET)35％的平均值。

S3、利用最大互信息系数法分析有效的净生态***CO₂交换量与环境驱动变量数据的相关性，从环境驱动变量数据中筛选与净生态***CO₂交换量相关性较大的前N个主控因子；以N个主控因子为输入量，对应的净生态***CO₂交换量为输出量，构建有效数据库；

6、步骤S3具体为：

S31、将净生态***CO₂交换量X和任一环境驱动变量数据Y按照数据对的形式D＝{(X₁,Y₁),(X₂,Y₂),···,(X_n,Y_n)构建一个散点图，然后在水平和竖直方向上分别划分a和b个网格，一组a和b的值对应一种划分尺度，而每一种划分尺度对应多种划分方案，每一种划分方案下变量X和Y的互信息值为：

S32、通过步骤S31计算所有的划分方案下的互信息值，取最大值即为该划分尺度的最大互信息系数：

I^*(X,Y)＝max(I(X,Y))；

S33：对步骤S32所述的公式采用log₂ min(a,b)进行归一化处理，即可获得X和Y在该划分尺度上的归一化信息值：

S33、不同的a和b对应不同的M(D)，计算所有划分尺度的M(D)构成信息值特征矩阵M(D)_a,b，则X和Y的最大互信息MIC的值即可通过计算信息值的特征矩阵的最大值获得：

其中，B*(n)＝n^0.6是网格数的上限，n是数据量的大小；

S34、计算所有环境驱动变量数据与净生态***CO₂交换量之间的MIC值，优先选择与净生态***CO₂交换量之间MIC值较大的环境驱动变量数据；MIC的物理意义等同于传统回归分析方法中的决定系数(R²)，MIC(X,Y)的取值范围为[0,1]，当MIC(X,Y)值为1时，表示变量Y与变量X完全相关，当MIC(X,Y)值为0时，表示变量Y与变量X完全独立。

S35、计算所有环境驱动变量数据之间的MIC值，最终确定N个NEE的主控因子[x_i|i＝1,…,N]。

S4、对有效数据库中的每一项数据作归一化处理，得到归一化后的数据；将将归一化后的数据分为训练集和测试集；

本发明中以归一化后的数据集的75％作为训练数据库，剩余的25％作为验证数据库；

需要说明的是，步骤S4中在进行归一化处理时采用以下公式：

x'＝(x-x_min)/(x_max-x_min)

其中：x’为归一化后样本值，x为样本真实值，x_max为样本最大值，x_min为样本最小值；所述样本真实值包括从数据库中的带入数据。

S5、建立净生态***CO₂交换量支持向量机预测模型，利用训练集训练所述支持向量机模型，得到训练完成的模型；

需要说明的是，结合支持向量机理论和步骤S4得到的训练数据库[(x_i,j,y_j)|i＝1,…,N；j＝1,…,l₁](x_i是得到的影响NEE的主控因子，y_i是对应的NEE值，l₁是训练数据库的样本总量)，可以建立NEE的预测模型，其表达式如下：

其中f(x)是预测值，

是非线性映射函数，w和b是待调整系数，可通过求解如下被约束的二次最优化问题得到待调整系数(w和b)：

约束条件为：

式中，c是惩罚因子，p为描述模型解的稀疏性的敏感损失系数，ξ和ξ*是松弛变量。引入拉格朗日乘子(α，α*)可以将所述最优化问题转化为双偶问题：

约束条件为：

式中，K(x_i,x_j)为核函数。

支持向量机核函数选择最适用于非线性预测的高斯径向基(RBF)核函数：

g为核函数参数。

待调整参数w可以由拉格朗日乘子和核函数计算得到，其表达式为：

最终，NEE的SVR预测模型为：

S6、利用人工蜂群优化算法优化所述训练完成的模型，得到最优预测模型；

上述三个参数(p，c，g)对SVR模型的预测精度具有重要的影响，需要对三个参数进行同时寻优得到最优参数组合，通过所述的最优参数组合建立最优化的NEE预测模型，本发明采用人工蜂群优化算法进行参数寻优，具体步骤如下：

S61、采用十折交叉验证方法，构建人工蜂群优化算法的优化目标函数(f_u)，计算公式如下：

u是三个参数(p，c，g)的取值，f_u是训练数据库的交叉验证均方误差(CVMSE)，M是训练数据库的数量，S是交叉验证将训练数据库划分的分数，G_S是用于验证的第S份数据，y_i是通量塔测得的实际值，f(x_i)|u是当(p，c，g)等于u时通过SVR模型计算得到的NEE的预测值。

S62：对算法进行初始化设置，设置最大迭代次数maxCycle、蜜蜂总数NP、食物源数量SN(NP/2)、维数D、个体最大更新次数limit以及待优化参数c，g，p的取值范围，采用人工蜂群优化算法即可获得参数组合(p，c，g)的最优解；

S63：人工蜂群优化算法共包含三个阶段，即雇佣蜂阶段、观察蜂阶段和侦查蜂阶段。在雇佣蜂阶段，根据公式v_ij＝u_ij+θ_ij(u_ij-u_kj)生成NP个备选解，并通过步骤S61得到的f_u对备选解进行评估，然后通过轮盘赌的方式更新雇佣蜂的解；在观察蜂阶段，从根据

公式选择的解中生成NP个备选解，并通过步骤S61得到的f_u对备选解进行评估，然后通过轮盘赌的方式更新观察蜂的解；在侦查蜂阶段，如果存在某个解在迭代limit次后仍没有被更新，则采用公式u_ij＝u_imin+rand(0,1)(u_jmax-u_jmin)对其进行更新。

S64：记录每个循环中最小的CVMSE对应的参数组合(p，c，g)；

S65：重复步骤S63和步骤S64，直到达到终止条件，输出所有循环中最小的fu对应的参数组合(p，c，g)，所述参数组合即为最优参数组合.

S7、利用测试集测试最优预测模型，并对最优预测模型进行指标评价，得到最终模型；

S8、利用最终模型进行净生态***CO₂交换量预测。

需要说明的是，在最后进行预测时，还需要将归一化后的预测结果进行反归一化，反归一化处理时采用以下公式：x＝x'(x_max-x_min)+x_min。

本发明所述的一种净生态***CO₂交换量的预测方法中，通过收集通量观测塔的开路涡度相关***监测到环境驱动变量和NEE数据，采用MIC–SVR预测方法，有效的筛选预测模型的特征变量，构建环境驱动变量与NEE之间SVR模型，提高对NEE预测精度。

最后，作为一种实施例，本发明具体实验过程如下：

L1、对建立在神农架大九湖3号湖附近的通量观测塔的数据进行收集，数据时间从2018年1月–2018年12月。数据主要包括CO₂净生态***交换量NEE原始数据及12个环境变量(Ta、RH、PPFD、RN、SWIN、VPD、Ts₁、Ts₂、Ts₃、SWC₁、SWC₂、SWC₃)的原始数据。

L2、对采集到的原始数据进行质量控制，采用MATLAB按照前述方法对原始数据进行剔除，最终获得2018年的有效数据集。

表1 2018年神农架大九湖亚高山泥碳湿地半小时尺度NEE通量缺失情况

L3、基于上述得到的有效数据集，通过最大互信息系数MIC方法分析环境变量与NEE的相关性，根据计算得到的MIC值并结合权利要求3所述的原则选择NEE的主控因子，最终确定六个主控因子：PPFD、Ta、Ts₃、VPD、RH、SWC₃。

表2最大互信息系数(MIC)计算结果

L4、以L3确定的主控因子为输入变量，以对应的NEE为输出变量，构建用于SVR建模的有效数据库[(PPFD_j,Ta_j,Ts_3j,VPD_j,RH_j,SWC_3j,y_j)|j＝1,…,l]，l是数据库的样本数量，共6364个样本；

L5、对所述数据库中的每一项数据进行归一化处理，以有效数据集的75％作为训练数据库(4773个样本)，剩余的25％作为验证数据库(1591个样本)；

L6、采用支持向量机机器学习算法，将步骤L5归一化处理后的数据作为输入变量，带入到支持向量回归公式，建立支持向量回归机模型；所述向量回归公式为：

L7、采用人工蜂群优化算法，对步骤L6建立的支持向量回归模型进行优化，建立最优化的NEE预测模型。

L8、请参考图2，图2是训练数据库预测结果；将L5的训练数据库带入L7建立的NEE预测模型，决定性系数R²为0.91，预测均方根误差Rmse为0.011。

L9、请参考图3，图3是验证数据库预测结果；将L5的训练数据库带入L7建立的NEE模型进行预测并反归一化，即可获得NEE的预测值，决定性系数R²为0.85，预测均方根误差Rmse为0.025。

综合来看，本发明的有益效果是：具体考虑了NEE与其各主控环境驱动变量之间的关系，采用筛选的主控因子作为NEE模型的输入变量，有效改善了模型的预测性能。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种净生态***CO₂交换量的预测方法，其特征在于：包括以下步骤：

S1、获取涡度相关***连续性观测数据；所述涡度相关***连续性观测数据包括：净生态***CO₂交换量和环境驱动变量数据；

S2、对涡度相关***连续性观测数据进行数据预处理，得到有效的净生态***CO₂交换量和环境驱动变量数据；

S3、利用最大互信息系数法MIC析有效的净生态***CO₂交换量与环境驱动变量数据的相关性，从环境驱动变量数据中筛选与净生态***CO₂交换量相关性较大的前N个主控因子；以N个主控因子为输入变量，对应的净生态***CO₂交换量为输出变量，构建有效数据库；

S4、对有效数据库中的每一项数据作归一化处理，得到归一化后的数据；将将归一化后的数据分为训练集和测试集；

S5、建立净生态***CO₂交换量支持向量机预测模型，利用训练集训练所述支持向量机模型，得到训练完成的模型；

S6、利用人工蜂群优化算法优化所述训练完成的模型，得到最优预测模型；

S7、利用测试集测试最优预测模型，并对最优预测模型进行指标评价，得到最终模型；

S8、利用最终模型进行净生态***CO₂交换量预测。

2.如权利要求1所述的一种净生态***CO₂交换量的预测方法，其特征在于：步骤S2中数据预处理具体为：根据预设阈值区间，对处于阈值区间范围外数据异常值进行剔除。

3.如权利要求1所述的一种净生态***CO₂交换量的预测方法，其特征在于：步骤S3具体为：

S31、将净生态***CO₂交换量X和任一环境驱动变量数据Y按照数据对的形式D＝{(X₁,Y₁),(X₂,Y₂),···,(X_n,Y_n)}构建一个散点图，然后在水平和竖直方向上分别划分a和b个网格，一组a和b的值对应一种划分尺度，而每一种划分尺度对应多种划分方案，每一种划分方案下变量X和Y的互信息值为：

S32、通过步骤S31计算所有的划分方案下的互信息值，取最大值即为该划分尺度的最大互信息系数：

I^*(X,Y)＝max(I(X,Y))；

S33：对步骤S32所述的公式采用log₂ min(a,b)进行归一化处理，即可获得X和Y在该划分尺度上的归一化信息值：

S33、不同的a和b对应不同的M(D)，计算所有划分尺度的M(D)构成信息值特征矩阵M(D)_a,b，则X和Y的最大互信息MIC的值即可通过计算信息值的特征矩阵的最大值获得：

其中，B*(n)＝n^0.6是网格数的上限，n是数据量的大小；

S34、计算所有环境驱动变量数据与净生态***CO₂交换量之间的MIC值，优先选择与净生态***CO₂交换量之间MIC值较大的环境驱动变量数据；

S35、计算所有环境驱动变量数据之间的MIC值，遵循选择的环境驱动变量数据之间相互独立的原则，最终确定N个NEE的主控因子[x_i|i＝1,…,N]。

4.如权利要求1所述的一种净生态***CO₂交换量的预测方法，其特征在于：步骤S5中净生态***CO₂交换量支持向量机预测模型的构建过程如下：

S51、构建初始净生态***CO₂交换量支持向量机预测模型表达式：

其中f(x)是预测值，

是非线性映射函数，w和b是待调整系数；

S52、通过求解如下被约束的二次最优化问题得到待调整系数w和b；

约束条件为：

式中，c是惩罚因子，p为不敏感损失函数，ξ和ξ*是松弛变量；

S53、引入拉格朗日乘子(α，α*)将所述二次最优化问题转化为双偶问题：

约束条件为：

式中，K(x_i,x_j)为核函数；支持向量机核函数选择最适用于非线性预测的高斯径向基(RBF)核函数：

g为核函数参数；

S54、待调整参数w由拉格朗日乘子和核函数计算得到，其表达式为：

最终，预测模型为：

5.如权利要求4所述的一种净生态***CO₂交换量的预测方法，其特征在于：步骤S6中得到最优预测模型的具体过程为：

S61、采用十折交叉验证方法，构建人工蜂群优化算法的优化目标函数(f_u)，计算公式如下：

u是三个参数(p，c，g)的取值，f_u是训练数据库的交叉验证均方误差(CVMSE)，M是训练数据库的数量，S是交叉验证将训练数据库划分的分数，G_S是用于验证的第S份数据，y_i是通量塔测得的实际值，f(x_i)|u是当(p，c，g)等于u时通过预测模型计算得到的净生态***CO₂交换量的预测值；

S62、对算法进行初始化设置，设置最大迭代次数maxCycle、蜜蜂总数NP、食物源数量SN(NP/2)、维数D、个体最大更新次数limit以及待优化参数p，c，g的取值范围，采用人工蜂群优化算法即可获得参数组合(p，c，g)的最优解。

Images