基于深度学习的草地生态系统净碳交换模拟

2020-06-29齐建东黄俊尧

农业机械学报 2020年6期

齐建东黄俊尧

(北京林业大学信息学院，北京 100083)

0 引言

草地生态系统覆盖陆地面积约20%，其净初级生产力约占全球陆地生物区净初级生产力的1/3，在地球碳循环过程中具有重要的作用[1]。受人类生产、生活以及环境变化影响，草地生态系统具有较大的不确定性，在生态系统生产力、群落结构和土壤的碳氮储量等方面，不同年份表现出较大的差异。因此，研究草原生态系统碳循环的主要过程以及影响机理是全球碳循环研究的关键问题之一[2]。

迄今，大多数生态系统碳通量模型模拟及有关环境因子响应的研究主要依赖于复杂的生物和物理过程的模型[3-4]，但生态系统潜在的非线性和非平稳特征会降低这些机理模型的精确性和可靠性。文献[5]认为，基于数据驱动的模型比基于过程的生态学模型更具优势。伴随着生态监测数据的不断累积，机器学习理论及技术在生态模型模拟上的有效性已获得广泛认可[6-8]，机器学习在估算不同空间和时间尺度的碳通量方面比基于传统生态学过程的模型表现更好[9]。如使用人工神经网络(Artificial neural network，ANN)对城市绿地生态系统的净碳交换(Net ecosystem exchange，NEE)进行模拟与分析[10]，使用支持向量机(Support vector machine，SVM)估算亚洲陆地生态系统的54个通量塔站点初级生产量(Gross primary productivity，GPP)和NEE[11]，使用随机森林(Random forest，RF)分析光利用效率在模拟生态系统初级生产量GPP中的重要性[12]，使用自适应神经模糊推理系统(Adaptive network-based fuzzy inference system，ANFIS)和广义回归神经网络(General regression neural network，GRNN)预测森林生态系统的日碳通量[13]，使用模型树集合(Model tree ensemble，MTE)进行土地生物圈GPP的模拟[14]等。

目前，已有研究主要使用传统机器学习算法，大多针对森林、农田生态系统，针对草地生态系统的研究较少。深度学习理论及技术作为当前机器学习领域最活跃的研究分支，在层次特征学习与表征上更具优势，但其在生态模型模拟领域的应用几乎为空白[15]。从计算科学角度出发，生态领域中碳通量预测、插补等问题属于时间序列分析问题，而在时间序列预测方面，深度学习领域中长短期记忆网络(Long short-term memory，LSTM)的优秀性能已被证明[16]。其中，注意力机制[17]的引入进一步显著提升了时间序列分析效果。

本文对高寒地区草原生态系统NEE进行模拟，使用随机森林算法对影响其变化的各环境因子重要程度进行计算及分析，并验证深度学习理论及技术在生态模型模拟中的有效性。

1 材料和方法

1.1 研究区概况

研究地点位于内蒙古自治区与华北平原之间(42°27′ N，116°41′ E，海拔1 380 m)，属于半干旱农牧交错带典型地区。根据多伦气象站的长期气候资料(1994—2004年)显示，多伦年平均气温为3.3℃。年平均降水量为399 mm，集中分布在7—8月。年蒸发量为1 600～1 800 mm，平均蒸发量为1 748.0 mm，是降水量的4～5倍。最低、最高的月平均气温为1月的-17.5℃和7月的18.9℃。土壤类型为栗钙土，土壤pH值为7.12，有机质含量(质量比)为12.28 g/kg，主要物种为克氏针茅(Stipakrylovii)[18]。

多伦通量测量系统于2005年6月成立，利用涡动协方差技术测量CO2体积分数、含水率、显热通量和潜热通量。涡动协方差传感器包括Li-7500型CO2/H2O红外分析仪和CSAT-3型三维声波风速计。CO2/H2O传感器头安装在声波风速计的下风向处，分析仪每6个月校准一次。用距离地面3.5 m的CNR-1型净辐射计测量净辐射。通过LI-190SB型PAR量子传感器测量光合有效辐射。用HMP45C型温度和相对湿度探头测量空气温度和相对湿度。将CS616型TDR(时域反射)传感器插入0～10 cm土层测量土壤含水率，并用3个深度(5、10、20 cm)的热电偶测量土壤温度。用3个HFT-3型土壤热传感器在土壤表面以下5 cm处测定土壤热通量。用在地面以上4 m处的TE525MM型翻斗式雨量计测量降水量。用CR23XI型数据记录器以30 min的间隔记录小气候数据，以10 Hz对风速、声波温度、CO2体积分数和含水率的波动进行采样，并通过CR5000型数据记录器记录[19]。

1.2 长短期记忆网络

LSTM[20]是一种全互连的神经网络，神经元间存在反馈连接，在学习新特征的同时能够记忆历史信息。

每一个LSTM单元都由输入门、输出门、遗忘门和记忆单元组成。LSTM通过这些门控单元来控制历史信息对当前信息的影响，使得模型能够较长时间地保存并传递信息。

t时刻，LSTM的输入门it、遗忘门ft、输出门ot、记忆单元的状态ct、隐藏层信号ht分别为

it=σ(Wiht-1+WiXt+bi)

(1)

ft=σ(Wfht-1+WfXt+bf)

(2)

ot=σ(Woht-1+WoXt+bo)

(3)

ct=ftct-1+ittanh(Wcht-1+WcXt+bc)

(4)

ht=ottanh(ct)

(5)

其中

Xt=[x1x2…xn]

式中Wi、Wf、Wo、Wc——输入门it、遗忘门ft、输出门ot和记忆单元状态ct的权重矩阵

bi、bf、bo、bc——输入门it、遗忘门ft、输出门ot和记忆单元状态ct的偏移量

Xt——t时刻监测指标所构成的矩阵

结构如图1所示。

1.3 编码器-解码器模型

编码器-解码器网络[21]是一个Seq2Seq模型，编码器接受一个序列作为输入，并将输入序列中的信息编码作为中间表示，解码器将其解码输出另外一个序列。

在编码器层，将通量塔主要监测指标矩阵Xt=[x1x2…xn]输入到编码器网络中

ht=Fencoder(ht-1,Xt-1)

(6)

其中Fencoder代表上述e个LSTM单元连接而成的网络，并将编码器层所有LSTM的输出信号汇总进行全连接输出语义向量

c=f(h1,h2,…,ht)

(7)

在解码器层，将编码器层的语义向量以及t-1时刻的碳通量yt输入到解码器网络中

st=Fdecoder(st-1,yt-1,c)

(8)

其中Fdecoder表示上述d个LSTM单元连接而成的网络，将解码器层LSTM的输出信号以及编码器层的语义向量c进行全连接得到t时刻碳通量

yt=F(st,c)

(9)

1.4 注意力机制

上述过程在式(7)中将历史信息映射为单一向量的行为造成了编解码过程中信息细节的大量丢失。注意力机制[17]根据编码器层e个LSTM的隐藏层信号hj(j=1,2,…,e)计算e个语义向量αj，并根据t-1时刻解码器的输入信号st-1，与t-1时刻的NEE值，得到解码器第j个解码器的输出lj。

(10)

(11)

lj=f(st-1,ht)

(12)

注意力机制打破了编码器-解码器结构在编解码时都依赖于内部一个固定长度向量的限制，结构如图2所示。

1.5 随机森林模型

随机森林[22]是一种由多棵决策树组成的组合模型。首先，对原始输入检测指标形成的样本矩阵利用自助重抽样(Bootstrap re-sampling)形成n个样本；其次，对每个样本分别建立较弱的分类回归树(Classification and regression tree，CART)[23]模型；最后，以n棵回归树结果的平均值作为因变量的预测值。回归过程如图3所示。

图2 基于注意力机制的编码器-解码器网络 Fig.2 Attention-based encoder-decoder

图3 随机森林回归架构 Fig.3 Random forest regression framework

图4 基于注意力机制的NEE模拟分析模型Fig.4 Attention based net ecosystem exchange model for forecasting and analysis

随机森林中使用二分递归分割技术CART作为基分类器，在每个节点(除叶节点外)将当前样本集分割为两个子集节点。通过每个节点的基尼杂质指数(Gini impurity)来计算监测指标的重要性评分(Variable importance measure，VIM)[24]。

基尼杂质指数是一种基于杂度的属性分裂方法。杂度与可获取的有用信息呈负相关。杂度越小，变量的离散程度越差，得到的信息量越大。基尼杂质指数的计算公式为

(13)

式中c——样本类别数

pi——特征中属性a对应的样本属于类别i的概率

1.6 基于注意力机制的NEE模拟分析模型

图4为基于注意力机制的NEE预测分析模型，包括输入层、模型层和输出层。输入层以内蒙古多伦草原通量塔各项监测指标为自变量，利用随机森林实现因子选择，并对数据归一化，划分训练集、测试集等；输出层实现数据反归一化，计算模型输出等；编码层将通量塔检测指标编码为若干个语义向量和历史碳通量数据共同输入解码层，自动挖掘监测指标和历史碳通量数据的相关特征并进行模拟，之后利用随机森林模型分析其各监测指标影响NEE的重要性。

2 实验结果与分析

实验数据来源于全球通量观测网络FLUXNET 2015数据集[25]中国内蒙古多伦草原通量站点2007年1月1日至2008年12月31日的30 min监测数据，采用边缘分布抽样[26](Marginal distribution sampling，MDS)进行数据插补，共计731 d，35 088条数据，其中使用23 494条数据作为训练集，11 572条数据作为测试集，去除无关指标、标志位指标以及不同计量的重复指标后，通过随机森林模型计算剩余监测指标的重要性得分，过滤掉潜热通量(H)、散热通量(LE)、大气压力(PA)和风速(WS)等重要性得分较低的监测指标后，选取重要性得分最高的6个监测指标，分别为土壤温度(Ts)(27.0)、光量子通量密度(PPFD)(19.9)、降雨量(P)(15.4)、土壤含水率(SWC)(17.6)、饱和水汽压差(VPD)(18.3)以及空气温度(Ta)(10.6)，作为模型的输入。

为消除不同性质、不同量纲的指标对综合评价结果的影响，需对各指标值采用(0，1)标准化方法对数据进行归一化处理，即

(14)

式中xnormal——归一化后数据

x——原始数据

xmin——x最小值xmax——x最大值

采用已在NEE预测方面广为使用的两种机器学习模型人工神经网络(ANN)与支持向量机(SVM)，与3种在时间序列预测方面取得较好效果的深度学习模型递归神经网络(Recurrent neural networks，RNN)、LSTM、编码器-解码器网络(Encoder-decoder)与注意力机制(Attention)进行模型预测效果的对比。

实验运行环境为：操作系统使用Ubuntu 16.04 Server，GPU为GTX 1050ti，处理器为Intel i7-8750H，内存为8.0 GB。在注意力机制模型中，设置输入历史序列窗口的大小T=10，编码器层LSTM隐藏层的大小为128，解码器层LSTM隐藏层的大小为64，初始学习率为0.000 1，共训练100个epochs，每10个epoch学习率衰减为原来的90%，batch size设置为128。

以预测值和真实值的均方误差(Mean square error，MSE)、平均绝对误差(Mean absolute error，MAE)以及决定系数R2为损失函数分别执行30次实验，结果表明，以均方误差为损失函数的效果最好。因此选择均方误差为损失函数。

在ANN中，为避免梯度消失，设置隐藏层大小为6，激活函数为sigmoid函数，训练50 epochs。在SVM中，设置惩罚因子为1.0，选择RBF核函数，设置多项式次数为3。

2.1 NEE预测模拟结果

图5 各模型预测值与实际值对比Fig.5 Comparison of original and predictive values

如图5所示，深度学习模型均能在不同程度上较好地拟合NEE，预测变化趋势与实际值变化趋势较为一致，但预测效果存在一定的差异。就整体预测结果而言，加入注意力机制的预测结果比其他模型的拟合效果更好，变化起伏更小，预测效果更稳定，全局无较大起伏的波动点。

加入注意力机制后预测效果比编码器-解码器网络略有提升，长短期记忆网络只能够大致模拟出碳通量实际值的趋势。如图6所示，各预测模型均不同程度地在04:00—08:00和16:00—20:00时段内呈现相对较大的起伏，这是因为在这两个时段植物发生光合作用，碳通量的变化幅度最大。由于碳通量随土壤温度等的变化具有延时性，在变化幅度较大的时间段中，预测的难度相对较大。

图6 加入注意力机制后各模型预测值与实际值对比 Fig.6 Comparison of original and predictive values after adding Attention mechanism

本文选择均方根误差(Root mean square error，RMSE)、MAE、平均绝对百分误差(Mean absolute percent error，MAPE)、R24项指标评估模型预测的准确性。

如表1所示，深度学习模型的精度指标相比传统机器学习模型ANN、SVM而言提升显著。而注意力机制模型RMSE为0.28 μmol/(m2·s)，MAE是0.18 μmol/(m2·s)，R2为0.93，与普通深度学习模型长短期记忆网络和普通编码器-解码器网络相比更精确。

表1 不同模型预测结果Tab.1 Evaluation indices of different models

2.2 NEE预测稳定性结果

图7表明，注意力机制的预测精度比普通深度学习模型编码器-解码器网络与长短期记忆网络更稳定。以多伦草原的10次预测结果RMSE为例，标准差为0.002 μmol/(m2·s)，相比编码器-解码器网络和长短期记忆网络下降0.005、0.036 μmol/(m2·s)。因此，基于注意力机制的碳通量预测模型能够更稳定地进行预测。

图7 模型稳定性对比 Fig.7 Comparison of stability of models

综上所述，基于注意力机制的模拟模型在长期NEE预测上不仅在预测精度上表现出良好的性能，并且具有较强的稳定性，可以挖掘相对于传统机器学习模型更深尺度的特征信息。

2.3 NEE主要环境因子响应过程

2.3.1NEE日变化及季节变化

本文选取生长季5—9月NEE的数据进行分析。如图8所示，在生长季NEE日进程呈现为单峰型，而在非生长季呈现为双峰型趋势[27]。昼间的NEE为负值(碳汇)，夜间的NEE为正值(碳源)。各月之间存在显著差异，生长季初期5月昼间生态系统碳吸收与夜间生态系统的呼吸都较弱。在生长旺季的7、8月，昼间生态系统碳吸收与夜间的呼吸都达到最强。另外，生态系统碳吸收的最强值均出现在10:00—12:00，随后下降，18:00以后，NEE逐渐变为正值，由碳吸收转变为碳排放。生态系统碳吸收的日最强值为净碳交换量-3.51 μmol/(m2·s)，夜间NEE的最大值为2.05 μmol/(m2·s)。

图8 生长季与非生长季不同时期NEE的日变化曲线Fig.8 Net ecosystem exchange in day-scale of growing and non-growing seasons in 2007

各月NEE的变化曲线如图9所示，1—3月均表现为碳吸收，4—7月表现为碳排放，在7月初时达到一年中的最大值，为1.90 μmol/(m2·s)，7—8月间生态系统的NEE开始下降，逐渐由碳源转换为碳汇，8月初达到碳吸收的最大值，净碳交换量为-1.01 μmol/(m2·s)。9月初出现碳排放的最高值，之后随着温度下降碳排放也逐渐降低。从全年来看，该草原生态系统表现为碳汇，全年的碳吸收量为19.06 g/m2。

图9 2007年不同时期净NEE的月份变化曲线Fig.9 Net ecosystem exchange in month-scale of growing and non-growing seasons in 2007

2.3.2环境因子的季节变化

2007年全年观测的光量子通量密度(PPFD)、空气温度(Ta)、土壤温度(Ts)、饱和水汽压差(VPD)、土壤含水率(SWC)和降雨量(P)的变化如图10、11所示。通常，PPFD受太阳高度角、天气、海拔影响。5—7月太阳高度角较高，PPFD达到峰值，最大达到700 μmol/(m2·s)。6月受到降水影响，PPFD略有降低。8月之后PPFD随着太阳高度角的降低而逐渐下降，在12月末达到最小值。空气温度和土壤温度具有相同的季节变化趋势，最高值出现在7月，最低值出现在1月，土壤温度比空气温度平均高3.14℃，空气温度和土壤温度变化幅度分别为-20.90～23.82℃与-12.48～26.16℃，年平均值分别为3.71、6.85℃。年总降雨量为258.01 mm，其中5—9月生长季的总降雨量为227.64 mm，占全年降雨量的88%左右。土壤含水率的变化主要受到降雨量的影响。在3—9月保持着较高的含水率。VPD同样具有明显的季节趋势，在6月达到最高值，为25.34 hPa，在1月达到最低值，为0.32 hPa。

2.3.3NEE对辐射的响应

图10 PPFD、Ta、Ts、VPD、SWC和P的月份变化曲线Fig.10 Photosynthetic photon flux density, air temperature, soil temperature, vapor pressure deficit, soil water content and precipitationin month-scale of growing and non-growing seasons in 2007

图11 PPFD、Ts、Ta、VPD的日变化曲线Fig.11 Photosynthetic photon flux density, soil temperature, air temperature and vapor pressure deficit in day-scale of growing and non-growing seasons in 2007

图12 NEE对PPFD的响应Fig.12 Response of net ecosystem exchange to photosynthetic photon flux density

如图12所示，4—9月上、下午草原生态系统的NEE与PPFD的关系表明，草原生态系统的NEE与PPFD之间有很好的直角双曲线关系[28]。而非生长季上、下午生态系统的NEE与PPFD之间没有明显的曲线关系。生长季各月上、下午NEE随PPFD的变化速率存在明显差异，各月都是上午的NEE对PPFD的响应强于下午。在植物生长初期的5月，NEE随PPFD的变化较缓慢，上、下午生态系统净光合速率的差异较小。6月的净光合速率明显高于5月，上、下午的净光合速率差值与5月相当，但当PPFD超过1 250 μmol/(m2·s)时，上、下午的差值略有增大。然而，在生长旺盛季节的7、8月，生态系统的净光合速率达到最大。9月上、下午生态系统的净光合速率与6月相似。

2.3.4NEE对温度的响应

NEE与土壤温度的关系如图13所示。当土壤温度在0℃以下时，土壤微生物活性很低，NEE随土壤温度的变化不大；土壤温度在0～15℃时，碳排放整体随土壤温度的升高而增大；土壤温度在15℃以上时，NEE随着土壤温度的升高而急剧增加。如图14所示，非生长季温度系数Q10明显大于生长季，在1月达到最大值，为5.96。

图13 NEE对Ts的响应Fig.13 Response of net ecosystem exchange to soil temperature

图14 温度系数Q10的月度变化Fig.14 Change of temperature coefficient Q10 in month-scale

2.3.5NEE对水分的响应

图15为NEE与VPD的关系。在VPD小于5 hPa时，NEE增加速率较大。VPD为5～20 hPa时，NEE增长速率稍缓。VPD大于20 hPa时，NEE随着VPD的增加而减少。

图15 NEE对VPD的响应Fig.15 Response of net ecosystem exchange to vapor pressure deficit

3 讨论

3.1 辐射对NEE的影响

光合辐射是生态系统进行光合作用的能源。总体上，生长季NEE随着PPFD的增加而增加，呈现直线双曲线趋势[29]。5月生长初期，由于植物尚未发育完全，NEE随PPFD的变化较缓慢。在6—7月，上午NEE随着PPFD的增加而增加，在午后，生态系统NEE随着PPFD的增加呈升高趋势，这是温度与辐射共同作用的结果，并且生态系统逐渐从碳排放转化为碳吸收。高原上的强辐射会使植物的光合器官受损[30]，为避免这种状况发生，植物的对策是被动关闭气孔，导致光合速率下降，造成植物的“午休”现象[31]。在7月，PPFD处于一年的峰值，可以看到下午的NEE在午后同比上午有所减少。在8—9月，在相同的PPFD条件下，生态系统光合速率随着温度的升高而升高，碳吸收受到PPFD的作用最强，但与此同时，呼吸作用也随着温度的升高而增强，因此下午的碳吸收相比于上午略减弱。10月后非生长季PPFD开始逐渐回落，并失去双曲线趋势。

3.2 温度对NEE的影响

温度同时影响生态系统的光合作用和呼吸作用。多伦草原海拔较高，生态系统的温度相对较低，在生长旺季的7月日平均温度为20.00℃，而在生长淡季日最低温度可达到-20.91℃。低温环境导致植物分解速率变慢，土壤中有机碳含量通常较高[32]。

除此之外，多伦草原生态系统的昼夜温差大，有利于生态系统的碳积累[33-34]。在生长旺季7—8月，昼夜温差可达到约23℃。其中6月植被未发育完全，相对受PPFD的影响较小，上下午温度差距显著，下午温度可达到25℃左右，NEE相比增加。

午后温度的升高在增加光合作用的同时也增加呼吸作用[35]。其中由图10可知，在9月，此时温度与8月基本持平，但PPFD有下降的趋势，因此相比光合作用呼吸作用的效果体现得更强，因此碳吸收的效果相较8月而言明显减弱。

3.3 水分对NEE的影响

水分是影响碳收支的重要因素之一。本研究的年总降雨量为258.01 mm，且主要集中于生长季5—9月。由图10可知，该期间土壤含水率(SWC)一直维持稳定水平，因此，在生长季水分不是影响生态系统NEE的主要因素。VPD是空气中水汽含量的描述，VPD较少时，NEE随着VPD的增加而增加，而当VPD达到一定程度时，NEE的增加速率放缓。VPD继续增加，致使植物气孔关闭阻塞，影响植物呼吸，NEE随着VPD的增加而迅速减少[36]。

多伦草原辐射较强，而温度和VPD相对较低。一年中生长季5—9月的辐射、温度、VPD处于高值。根据图11可知，午后的温度和VPD均处于较高水平，有利于生态系统的光合作用，但同时温度的升高也加剧系统的呼吸作用。

3.4 NEE环境因子重要性的季节性变化

各环境因子共同影响生态系统的碳收支相互制约。图16为各环境因子之间对NEE影响重要性的月份变化曲线。可以看到，由非生长季向生长季过渡的3—4月间，PPFD(33.5)与VPD(30.0)对NEE的变化起主导作用。进入生长季后的5—6月间，对比图10可以发现此时P(5.4)相对较少，因此SWC(50.5)是NEE的变化的主要影响因素。进入7月，P(3.8)依然较少，因此PPFD(26.8)与SWC(60.1)协同作用NEE的变化。8月PPFD(2.8)与SWC(6.9)由于7月底的降雨而变得相对充足，因此VPD(41.5)与P(42.7)成为影响NEE的主要因素。9月后PPFD与P均急剧减少并维持稳定，并且根据图14可知，非生长季温度系数Q10较生长季高，因此非生长季NEE主要与生态系统的呼吸直接关联，在非生长季1—3月Ts(44.6)与10—12月Ts(54.2)通过影响植物呼吸成为影响NEE的决定性因子。