基于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的增加，會對植物葉片氣孔閉合產生壓力，影響其光合和蒸騰作用；對溫度的偏導數說明，隨著溫度的增加，光合速率逐漸大于呼吸速率。