基于XGBoost-ANN的城市绿地净碳交换模拟与特征响应

齐建东 黄金泽 贾 昕

(1.北京林业大学信息学院， 北京 100083； 2.北京林业大学水土保持学院， 北京 100083)

0 引言

近年来，随着城市化进程的不断推进，城市绿地面积也在不断增加，截止2016年，北京城市绿地总面积约为3×104hm2，人均公园绿地面积为16.1 m2，城市绿地覆盖率达48.4%[1]。人工林植被作为城市绿地的重要组成部分，对气候具有重要的调控作用。定量分析城市绿地净生态系统碳交换(Net ecosystem exchange，NEE)数据不仅能够促进人们对区域碳源、汇功能的理解，还可为研究不同生态系统对于全球气候变化的反馈机制、预测区域气候变化提供参考。

NEE数据具有明显的季节性特征，并且与温度、光合有效辐射等各类气象、环境因子存在复杂的非线性关系[2-3]，因此，模型模拟难度较大，难以保证模拟效果。虽然通过试凑法或者经验选择法选择环境因子往往能获取较高的精度，但过于复杂的生态学意义不明晰，且不利于模型推广。随着机器学习技术的不断进步，以随机森林、XGBoost和人工神经网络(Artificial neural network，ANN)为代表的机器学习算法能够有效地克服上述缺陷，广泛适用于生态学领域。在城市绿地净碳交换模拟中，王宏莹[4]使用BP-ANN对NEE进行插补，评价徐州南部城区碳源、汇情况，MENZER等[5]通过对比3种不同的ANN模型，探究风速与风向对于城市生态系统中NEE的影响。虽然由于ANN强大的非线性映射能力而被广泛使用，但是在输入因子选择等方面则需要研究人员凭经验确定[3]，因此具有不确定性。XGBoost模型是CHEN等[6]于2016年提出的模型，该算法不仅可以通过计算输入因子的相对重要性对结果进行解释，而且通过内置交叉验证等方法有效防止模型过拟合，使计算结果更为科学可靠，目前已经被广泛用于空气质量预报[7]等领域。充分利用XGBoost对结果的可解释性，将其与ANN模型相结合，可以有效弥补ANN在因子选择方面的缺陷。

在城市生态系统中，地表接收的太阳辐射强度和空气温度和湿度、土壤含水率以及风速等环境因素受到公园内植物多样性、人工林管理方式以及建筑布局等人类活动的影响。这些城市特征性使生态系统对环境因子的响应发生了巨大变化，已经引起科学家广泛关注。目前已经有大量学者开始探究城市生态系统中影响NEE的主要因子，认为NEE对环境因子的响应主要与光合有效辐射、温度、风速、风向有关[8-9]，但是对于城市生态系统中NEE对环境因子的响应尚缺乏深入探讨。鉴于此，本文基于北京奥林匹克森林公园2013—2016年生长季白天利用涡度相关法连续观测的NEE数据，采用XGBoost方法分析空气温度(Ta)、土壤温度(Ts)、光合有效辐射(PAR)、风速(WS)、相对湿度(RH)、饱和水汽压差(VPD)、10 cm深度土壤含水率(VWC10)7个环境因子对NEE的影响程度，并对其进行选择和评价，利用ANN神经网络模型解释NEE对主要环境因素的响应结果。

1 材料与方法

1.1 研究区概况

研究区位于北京奥林匹克森林公园(40°01′N，116°23′E)，海拔51 m，主要土壤类型为潮褐土，4年均温变化范围为11.69～13.01℃，1月与7月均温变化范围分别为-5.22～-1.01℃与25.45～27.82℃，极端低温变化范围-19.00～-12.01℃，极端高温变化范围为36.58～39.82℃。该地区气候属于暖温带半湿润大陆性季风气候，四季分明，4年均降水量变化范围458.61～669.12 mm，无霜期208～225 d。降水季节分配不均，主要集中于6、7、8月(2013—2016年北京统计年鉴)。观测地内主要植被为人工营造的乔冠草复层景观林，乔木类代表物种包括油松(Pinustabulaeformis)、侧柏(Platycladusorientalis)、国槐(Sophorajaponica)、白蜡(Fraxinuschinensis)、银杏(Ginkgobiloba)，灌木主要为山桃(Prunusdavidiana)，丛生灌木主要为丁香(Syzygiumaromaticum)、地被植物主要为石竹(Dianthuschinensis)等[10]。

1.2 数据获取及预处理

实验数据通过涡度通量观测仪测量获取。仪器主要包括三维超声风速仪(CSAT3型，Campbell Scientific Ltd.，美国)、红外气体分析仪(EC155型， Campbell Scientific Ltd.，美国)、净辐射仪(CNR-4型, Kipp & Zonen Inc.，荷兰)、光量子传感器(PARLITE型, Kipp & Zonen Inc.，荷兰)、空气温湿度传感器(HMP45C型，Campbell Scientific Ltd.，美国)、土壤温度传感器(Campbell-109型，Campbell Scientific Ltd.，美国)及土壤含水率传感器(CS616型，Campbell Scientific Ltd.，美国)，分别用于测量风速、CO2/H2O密度脉动、辐射、空气温湿度以及土壤温度等微气象数据。数据采集器(CR3000型，Campbell Scientific Ltd., 美国)以10 Hz频率记录涡度通量观测仪的数据，经过野点剔除、二次坐标轴旋转、频率响应校正和WPL校正等操作后，在线计算30 min通量值[10]。

采用2013—2016年30 min通量塔NEE与微气象数据，由于恶劣天气、硬件设施故障以及人为因素的影响，数据存在异常值和缺失值。针对上述情况，参考中国通量网通量数据标准处理流程与净生态系统交换量标准处理流程[4,11]，将数据进行如下操作：

(1)利用3倍标准差法剔除野点，检查数据范围。

(2)存储通量计算。考虑到奥林匹克森林公园植被类型空间分布均匀，植被较高，具有涡度相关系统高度下冠层内存储通量不为零的特点。NEE定义为

PNEE=Fc+Fs

(1)

式中PNEE——净生态系统碳交换量

Fc——通量观测塔在植被上部观测值

Fs——涡度通量观测仪安装高度下冠层内存储通量

(3)使用分段平均值检验法计算的摩擦风速阈值为0.2 m/s，因此剔除夜间NEE数据中摩擦风速小于0.2 m/s的碳通量数据。

(4)考虑到不同量纲的数据序列会增加数据处理成本与模型拟合时间，对数据进行归一化处理。最后，选取4年中生长季(6—9月)白天[14]，即PAR大于10 μmol/(m2·s)的数据作为研究对象[12]。各年有效数据统计结果如表1所示。

1.3 模拟方法

1.3.1XGBoost模型与因子选择

XGBoost模型属于集成学习模型，通过构建并结合多个回归树模型来完成学习任务。首先通过自举法(Bootstrap)生成N个训练集。其次，对于每个训练集均建立回归树模型并进行训练。最后，计算所有回归树模拟结果，加权后作为输入变量的预测值。具体流程如图1所示。由于XGBoost模型在训练过程中根据每个输入因子的信息增益选择最好的特征进行分裂，因此，通过计算所有回归树中第i个输入因子xi出现的次数，可得到该输入因子在整个XGBoost模型中的重要性得分。计算公式为

表1 各年生长季白天NEE有效数据统计Tab.1 Summary of effective growing season daytime NEE data for each year 条

(2)

式中n——XGBoost中树的编号

Nt——XGBoost中树的数量

图1 XGBoost 回归原理图Fig.1 XGBoost regression principle

1.3.2ANN神经网络及其偏导数

神经网络主要是通过自学习寻找目标值与输入变量之间的映射关系。一个神经网络主要由输入层、输出层以及中间的隐含层构成(图2)。输入层负责接收输入数据，输出层负责输出整个神经网络的计算结果，隐含层则负责描述问题的层次关系。神经网络上每个节点称之为神经元，各层之间的神经元通过一定的权重相互连接。其基本原理为：当网络输出层的计算结果与期望结果偏差过大时，通过优化算法对每个神经元的权重进行更新。通常ANN神经网络使用反向传播算法作为优化算法，其本质是计算目标函数的梯度来寻求最优值。在数学上，梯度值表示目标函数在该点变化率最大的方向。在生态学中，则可以表示为NEE对于环境因子的响应速率[13]。

图2 ANN分析原理图Fig.2 Diagram of ANN analysis principle

本文采取自动调整学习率以及早停策略防止ANN过拟合或者欠拟合。

1.3.3基于XGBoost与ANN神经网络的NEE分析模型

基于上述算法，为了能够更为准确地分析影响城市生态系统NEE的主要影响因子与其响应关系，本文将两种算法结合进行分析。其分析流程如图3所示，首先对采集到的NEE数据进行质量控制和归一化预处理；其次利用XGBoost模型筛选出NEE主要影响因子，并将其作为ANN模型的输入因子；通过ANN模型对各输入因子的偏导数，探究NEE对各环境因子的响应关系，实现城市生态系统NEE对主要环境因子的响应分析。

图3 基于XGBoost和ANN的NEE分析技术路线Fig.3 CO2 flux analysis flow chart based on XGBoost and ANN

1.4 性能评估

为评估ANN模型的拟合效果，本研究中使用了目前最常用的评估标准[14]：决定系数R2、平均绝对误差(Mean absolute error, MAE)、均方根误差(Root mean square error, RMSE)和一致性系数(Index of agreement, IA)。

研究所用程序设计语言为Python 3.6(64-bit)，集成开发环境为Anaconda 3。程序设计中，XGBoost模型基于xgboost包实现，ANN模型则由Keras 2.2.0和TensorFlow 1.6.0完成编写。为保证结果可靠性，每组实验在相同条件下重复100次，实验结果取平均值。

2 结果与分析

2.1 模型参数优化

2.1.1XGBoost模型

在XGBoost模型中，主要参数有3个：每棵树的最大深度、树的数量、最小叶子节点权重之和。本文通过网格寻优策略测试了3种参数的240种组合，最终确定当树的最大深度为8、树的数量为1 500、最小叶子节点权重之和为200时，XGBoost模型的目标损失函数值最小，达到最优。

2.1.2ANN神经网络

本文设定初始学习率为1.5，最小学习率为0.001；迭代次数经早停策略优化后平均次数为473次；批大小为16，所用隐含层神经元数量为4，输入层与隐含层共用bias单元，如图2所示。另外，激活函数采用sigmoid函数。

2.2 环境因子变化特征

4年内生长季白天时间段内环境变化的月平均日变化基本特征如图4所示。本研究数据时间段，北京奥林匹克森林公园的PAR、Ta、Ts、VPD在12:00—14:00达到极值，PAR极大值出现在2013年7月16日，为1 533 μmol/(m2·s)，月平均日变化极值出现在2014年6月，为1 045 μmol/(m2·s)；气温的月平均日变化极小值为24.12℃，极大值为33.35℃；由于VPD主要受到水热条件以及植物蒸腾作用的影响，生长季期间饱和水汽压差明显波动，且变化剧烈。2014年的VPD显著高于其他年份；土壤含水率反映了区域降水情况[15]，2014年9月的降水使得VWC10指标与其他3年呈明显差异，VWC10与WS 4年内月平均日变化范围分别为16.54%～35.36%与0.93～1.56 m/s。

2.3 环境因子重要性得分

图5为计算得到的环境因子对NEE影响的重要性得分。可以看出，环境因子对NEE影响的重要性得分由大到小表现依次为PAR、VPD、Ta、RH、Ts、WS、VWC10。由此可知，PAR、VPD、Ta是影响奥林匹克森林公园植物生长季NEE变化的重要因素。

2.4 NEE模拟结果以及对主要环境因子的响应过程

表2为不同组合的输入指标在测试数据集上的评价结果。计算结果表明，随着环境因子输入数量的增加，R2总体逐渐增加，当输入因子为PAR、VPD、Ta、RH、Ts、WS、VWC10时，训练集R2为0.712，MAE为3.129 μmol/(m2·s)，RMSE为4.349 μmol/(m2·s)，一致性指数为0.911，测试集决定系数R2为0.748，RMSE与MAE分别为4.253、2.971 μmol/(m2·s)，IA为0.920，相较于其他组合为最优结果。训练集模拟值与观测值的拟合结果如图6所示，测试集的拟合结果如图7所示。

观察图5与表2结果可知，当输入因子依次加入VPD以及Ta时，模型测试集R2增加显著(R2从0.587提升到0.741)；若继续增加环境因子数量，R2提升效果并不明显。因此，本文选取PAR、VPD、Ta 3个环境因子进行分析讨论。

2.4.1PAR对NEP的影响

NEP(Net ecosystem productivity)为生态系统净生产力。由图5可知，在生长季，PAR是影响白天净生态系统交换量的决定性因素。图8显示了NEP随PAR的变化情况以及对PAR的偏导数，即表观量子效率。从图中可知，∂PNEP/∂PPAR>0(PNEP为生态系统净生产力，PPAR为光合有效辐射)，光合作用随着PAR的增大而逐渐增强，生态系统对于CO2的吸收量逐渐增大，碳汇能力增强。当PAR大于1 200 μmol/(m2·s)时，表观量子效率逐渐趋于0并保持稳定，说明此时光合有效辐射已经不是部分植被进行光合作用的主要因素。图8中，最大表观量子效率为0.087。

表2 不同输入因子模拟结果评估Tab.2 Evaluation indices of different combinations

注：Y为输入因子中包括该因子，N为未包括该因子。

图6 训练集观测值与模拟值的关系Fig.6 Relationship between observed and simulation results in train dataset

图7 测试集观测值与模拟值的关系Fig.7 Relationship between observed and simulation results in test dataset

2.4.2饱和水汽压差对NEP的影响

VPD可由RH、Ta通过公式估算得出，该指标能够体现出温度与相对湿度的特点[16]，因此，其重要性得分要高于RH、Ta。当VPD过高时，会对植物叶片气孔闭合产生压力，影响植物的光合和蒸腾作用。由图9可知，当VPD较小时，植被生长环境适宜的情况下，∂PNEP/∂PVPD>0(PVPD为饱和水汽压差)，NEP与VPD正相关，促进生态系统碳汇，但促进能力逐渐减弱。随着VPD的增大，研究区域内大部分植物光合作用受到抑制，但仍有一部分时刻的偏导数大于0，对NEP的变化起促进作用，但总体趋近于0。这说明即使在较高VPD时，由于研究区域内植物种类多样，其中包括侧柏等耐高温、抗旱性强植物，生态系统总体依然呈碳汇现象，但是随着VPD的继续增加，对植物叶片气孔闭合产生压力，影响植物的光合和蒸腾作用[8,22]。

图8 NEP对PAR的偏导数Fig.8 Numerical partial derivatives of daytime NEP response to PAR for each half-hourly data

图9 NEP对VPD的偏导数Fig.9 Numerical partial derivatives of daytime NEP response to VPD for each half-hourly data

2.4.3空气温度对NEP的影响

空气温度主要通过影响生态系统呼吸和光合作用来影响生态系统CO2的交换。在本文中，当温度大于15℃时，大部分时刻空气温度与NEP呈正相关，小部分呈负相关。在合适的温度范围内，随着温度的增加，逐渐转变为正相关(图10，PTa为空气温度)。该现象说明温度较低时，生态系统光合作用弱于呼吸作用，随着温度的逐渐升高，光合作用速率逐渐加快，大于呼吸作用速率。这与陈文婧等[10]关于奥林匹克森林公园温度与生态系统碳交换影响的研究结果类似。

图10 NEP对Ta的偏导数Fig.10 Numerical partial derivatives of daytime NEP response to Ta for each half-hourly data

3 讨论

3.1 样本大小及模型参数对模拟精度的影响

与随机森林类似，XGBoost模型的样本容量较小时，难以保证极高的模拟精度与泛化能力，因此机器学习模型需要有足够的样本数量保证结果的合理性。本文模型在训练、测试与独立验证3个阶段样本数量分别为5 544、1 109、2 376，样本容量较大。

模型模拟的结果除了与样本数量有关外，还受到模型参数的影响。在XGBoost模型中，树的深度、数量设置过大会导致模型运行时间过长，效率低下；若设置过小则无法有效地表达NEE的特征重要性，合理地设置最小权重之和可以避免模型学习到局部的特殊样本[6]。本文采用的网格寻优策略，是目前较常用的寻优算法，与机器学习模型相结合，广泛应用于煤炭灰熔点预测[17]，工矿复垦区土地利用分类等各个领域[18]。通过该算法优化XGBoost模型参数，可以较好地避免由于参数设置不恰当所造成的误差。

在神经网络模型中，除网络结构外，学习率、迭代次数对模拟结果有巨大影响。学习率过大，导致模型无法收敛，学习率过小，会导致模型收敛速度过慢等。本文采取学习率衰减的方法，并设定最小学习率为0.001，训练过程中每隔50代学习率降低为原来的1/2，直至达到最小学习率，该方法不仅加快了网络收敛速率，也较好地避免了振荡现象，已经被广泛应用于许多改进的BP算法中[19]；迭代次数对于结果的影响体现为：过多的迭代次数导致模型产生过拟合，过少的迭代次数则会使模型在未收敛时训练结束，导致欠拟合。本文通过观察测试集损失函数，采用早停策略[20]，保证在该参数配置下获取最优的拟合结果。

3.2 数据质量对模拟结果的影响

在数据收集过程中，由于大气降水、大气湍流、传感器自身误差等原因，造成数据出现缺失，极端噪声等问题，此类问题会对模型的精度造成极大的影响。在本文中，主要评估了包括光合有效辐射、饱和水汽压差在内的7个环境因素，未考虑生物因素，如气孔导度、水分利用效率等，因此无法更加有效模拟NEE的变化情况。另外，由于人类活动对城市生态系统的影响要远大于对其他生态系统的影响，因此，在未来的NEE的研究中，应当充分考虑植被生理因素、生长环境的气象因素以及城市交通[21]等因素的特点，选用更加精确的模型，以期提高NEE模拟精度。总体来说，XGBoost模型能准确地选取影响生态环境NEE的主要因素，ANN模型能够较好拟合NEE的大小以及变化趋势。

3.3 共线性变量对模型的影响

由图5可知，Ta、RH、Ts 3个环境因子重要性得分接近。为了进一步探究3个因子的主导性，本文评估了3种组合的模拟效果(表2，CIV.3-1，CIV.3-2，CIV.3-3)，模拟优度从大到小依次为CIV.3-1(PAR、VPD、Ta)、CIV.3-3(PAR、VPD、Ts)、CIV.3-2(PAR、VPD、RH)。说明在输入因子包括PAR和VPD的情况下，Ta最重要，其次为Ts，RH重要性最小。与图5结果存在差异，出现这一问题的原因是，XGBoost模型是决策树模型的组合，在模拟过程中输入因子共线性不会影响其预测能力，但是对数据的解释性影响巨大[22]。即当Ta特征被使用时，XGBoost模型对于Ts特征的权重将会减少，此时Ts相对于NEP来说，增加的有效信息少于RH，因此，在图5显示的重要性得分中，由大到小依次为Ta、RH、Ts。另外，CIV.4-1(PAR、VPD、Ta、RH)与CIV.4-2(PAR、VPD、Ta、Ts)组合的评估结果中，CIV.4-1的精度高于CIV.4-2，这一结果也佐证了上述结论。图11(图中PTs为土壤温度的值)展示了当输入因子组合为CIV.3-3时，NEP对Ts的偏导数，观察可知，其变化趋势与空气温度对于NEP的影响相同，但数值更小，对NEP的响应能力弱于Ta。综上所述，就单因子而言，重要性由大到小依次为Ta、Ts、RH；但综合考虑Ta后，RH对NEP影响大于Ts。

图11 NEP对Ts的偏导数Fig.11 Numerical partial derivatives of daytime NEP response to Ts for each half-hourly data

3.4 环境因子对NEE的影响

在本研究中，当PAR增大时，植物光合作用增强，生态系统中植被固碳能力增强，最大表观量子效率为0.087，高于其他生态系统中的计算值，出现这一现象是因为城市中气溶胶浓度较高使得散射辐射的比例增大，而碳收支对散射辐射敏感度较高所导致的[20]。当PAR大于1 200 μmol/(m2·s)后，其不再是影响光合作用的主要因素，如图8所示。PAR除了在城市生态系统中是NEE的主要影响因素外，在其他生态系统也起主要调控作用，WEN等[23]通过遗传神经网络筛选出PAR是湖南会同杉木人工林生态系统的决定性因子，在半干旱荒漠生态系统中，PAR对NEE的影响最为显著[24]，唐祥等[25]对北京八达岭NEE的研究表明，PAR是影响该生态系统生长季NEE变化的主导因子。除PAR外，VPD、Ta是影响生长季碳吸收的重要控制因素，该结论与陈文婧[8]针对城市绿地神态系统碳水通量研究结论一致，MOFFAT[20]对影响Hainich森林的环境因素重要性进行了等级划分，其研究结果也证明，VPD与Ta是重要的非辐射性环境因子，KUNWOR等[12]通过改进Michaelis-Menten关系式，在PAR与Ta的基础上，引入新的环境变量VPD，不仅提高了数据模拟的精度，还在一定程度上保留了模拟数据与观测数据的方差结构。另外，温度主要通过影响生态系统光合作用以及土壤微生物活性等方式，来影响生态系统NEE[26-27]。这也验证了PAR、VPD、Ta对于生态系统NEE模拟的重要性。综上，准确解释环境因子对于NEE的影响，有利于更好地认识植被对于区域气候变化的响应。

4 结论

(1)XGBoost-ANN模型能够较好地捕捉北京奥林匹克森林公园生长季NEE数据的变化特点，在训练集上R2为0.712，MAE为3.129 μmol/(m2·s)，RMSE为4.349 μmol/(m2·s)，IA为0.911；测试集R2为0.748，RMSE与MAE分别为4.253、2.971 μmol/(m2·s)，IA为0.920，说明其在模拟和分析森林NEE方面具有较好的适用性。

(2)模型测试结果表明，在训练过程中，通过机器学习自动化调参以及网格寻优等操作，可以避免由于参数设置不合理带来的误差。

(3)环境因子重要性得分表明，在考虑各因子间相互作用的情况下，影响生长季白天NEE的主要环境特征重要性由大到小依次为PAR、VPD、Ta、RH、Ts、WS、VWC10。该地区NEE变化主要受PAR、VPD、Ta 3个主要因素调控。

(4)各主要环境因子偏导数表明，北京奥林匹克森林公园生长季白天的光合作用表观量子效率最大值为0.087，并且当PAR大于1 200 μmol/(m2·s)时，PAR不再是影响NEP的主要因素；对VPD的偏导数说明，随着VPD的增加，会对植物叶片气孔闭合产生压力，影响其光合和蒸腾作用；对温度的偏导数说明，随着温度的增加，光合速率逐渐大于呼吸速率。