基于機器學習算法的九寨溝自然保護區植被生態水儲量定量反演研究

周湘山，楊武年，羅 可，張宇航，文 艷，唐曉鹿

（1.中國電建集團成都勘測設計研究院有限公司，四川 成都 610072；2.成都理工大學，四川 成都 610059）

0 前 言

森林生態系統約占地球陸地表面積的33%，被稱為“地球之肺”，是地球上最大的陸地生態系統，是全球生物圈中重要的一環，也是地球上重要的基因庫、碳貯庫、蓄水庫和能源庫，對維系整個地球的生態平衡有著不可替代的作用［1］。森林生態系統具有一項重要的生態服務功能就是涵養水源［2］，其涵養著陸地90% 以上的淡水資源［3］。林冠層是森林生態系統降水調節過程的第一個環節，也是森林生態系統中水分重新分配的第一個界面層［4］，林冠層對降水的調節作用主要表現為林冠截留和穿透雨兩部分；其次是樹干徑流，是指降雨沿著樹干流至樹木根部的過程，樹干徑流量在一次降雨中所占的比率比較低，一般低于5%；枯落物層作為森林生態系統的第二個水文作用層，主要功能是增加土壤的有機質含量、改善土壤的理化性質，提高土壤孔隙率，防止土壤嚴重板結、增強土壤入滲性能，減少土壤的無效蒸發［5］，使更多的水分能夠儲存在土壤中供植物生長；土壤層是森林生態系統的第三個水文作用層，在整個森林生態系統中水源涵養能力最強，不僅能夠調節降水資源的分配，也是生態系統水分和元素循環的主要儲蓄庫，是評價森林生態系統水源涵養能力的重要指標。研究不同森林植被類型儲水量及其水文功能，是生態修復、重建過程中的重要問題，也是需要迫切解決的基礎理論問題。

目前已有研究成果中，基本采用單一遙感影像，通過各類植被指數、紋理指數等地面參數與野外實測數據建立回歸關系反演或估算植被含水量。部分研究增加了LIDAR、INSAR、高光譜等多源遙感融合反演的成果，但利用機器學習算法建立回歸模型的應用較少。常用的機器學習算法包括：隨機森林、支持向量機、人工神經網絡等。機器學習算法較以往的線性參數模型有更好的擬合精度，對植被含水量數據需要滿足的統計條件較低，能夠有效地解決高維數據（變量個數較多、變量類型多樣）在建立模型過程中所帶來的變量選擇困難。隨著機器學習算法的發展及森林植被生態水相關遙感研究的深入，利用機器學習算法對森林植被進行分類、植被生態水反演模型構建，將進一步提高植被生態水儲量的估算精度和準確性。

1 研究區概況

研究區位于九寨溝漳扎鎮，漳扎鎮位于九寨溝縣南部，原名九寨鎮，地處N33°02′N～33°21′N，108°38′E～104°03′E，位于長江系嘉陵江上游白水江源頭。地勢整體呈南高北低，地貌以高山為主，海拔在1 800～4 800 m之間；氣候溫和，年平均氣溫為12.7℃。春季氣溫較低且溫差較大，平均氣溫位于9～18℃；夏季較涼爽；秋季氣候宜人但晝夜溫差較大；冬季較寒冷，日照充足，年平均日照約1 600 h。保護區內土壤大部分以半濕潤、半干旱暖的溫帶基帶山地淋溶褐土為代表，林下土壤則多為山地棕壤和山地暗棕壤，腐殖質含量較高。保護區內野生動植物資源非常豐富，已發現的野生動物有600余種，其中屬于國家重點保護動物有20多種。此外，保護區還是四川省的第二大林區，森林覆蓋率超過70%，有四川紅杉、白皮杉等多種國家級珍稀保護植物。

2 數據采集及處理

2.1 實測數據分析

2.1.1 野外數據采集

由于九寨溝自然保護區范圍的氣候條件復雜，云層覆蓋較多，考慮不同季節植被的含水量會有較大變化，而不同年份的相同季節植被含水量較之變化不大的因素，盡量保持下載的Sentinel-2B和Landsat8 OLI影像成像時間與所采集樣點實測數據的季節相契合，布點方式考慮了植被類型多樣和交通便利兩種因素。野外工作由A、B兩個小組一起完成，共采集54個樣地數據。樣地采用圓形樣地方式，以一點為圓心，10 m距離為半徑，對圓內所有胸徑大于5 cm的樹木進行單木檢尺，主要記錄樣地編號、樹種名稱、1.3 m高度處胸徑和樹高。其中針葉林類型14個，其主要樹種有油松、云杉、馬尾松、冷杉、云杉等；闊葉林類型16個，其主要樹種有櫟樹、樺樹、楊樹等；混交林類型24個。

2.1.2 生物量計算

本文采用模型法對樣地林木生物量進行估算，在樣地內測得的數據有胸徑和樹高，利用已有的生物量方程計算樣地內單棵樹木干枝葉各部位的生物量。

根據野外采集的胸徑D計算得出單棵林木的枝、干、葉生物量（干重），分別為BM枝、BM葉、BM干。由野外采集的林木的枝、干、葉少量樣本，可以得到分別對應器官的樣本鮮重FM枝、FM葉、FM干，帶回室內后將樣品放置烘箱中，將葉類樣本在70℃條件下烘干，枝干部分樣本在80℃條件下烘干，反復稱量至恒重，分別得到對應林木器官的干重DM枝、DM葉、DM干。

據生物量（干重）和含水率的關系，可以計算出林木的枝、干、葉植被含水量，并將森林植物樣地的含水量統一換算到標準單位噸/公頃（Mg/hm2）。

2.1.3 研究區植被含水量

按照54個樣地的實測數據，根據各類樹種生物量模型計算出各類樹種的干、枝、葉等部位的生物量，通過含水率和含水量計算出各類樣地所有樹種干、枝、葉和總的含水量。

2.2 遙感影像分類

本文主要提取研究區植被林種分布及面積，選擇計算機“解譯為主，目視解譯為輔”的信息提取方法。將研究區地物類型分為針葉林、闊葉林、針闊混交林、草地、灌木、水域、裸地和雪地，如圖1所示。對分類結果進行混淆矩陣精度驗證，其中總分類精度為81.73%，Kappa系數為0.795 6，滿足分類要求。

圖1 研究區土地利用類型分類

2.3 特征波段選取

本文結合植被光譜特性和Sentinel-2B、Landsat8OLI衛星傳感器參數參數，選取了Sentinel-2B影像 中 的 Band1、Band2、Band3、Band4、Band5、Band6、Band7、Band8、Band8b、Band9、Band11、Band12和Landsat8 OLI衛星影像中的Band1、Band2、Band3、Band4、Band5、Band6、Band7、Band9等原始波段地表反射率作為可選的特征因子。

2.4 植被指數提取

本研究主要選取了基于代表綠度的垂直植被指數（Perpendicular Vegetation Index）、三角植被指數（Triangle Vegetation Index）、差值 環 境植被指數（Difference Vegetation Index）、歸一化綠藍差異指數（Normalized Green-Blue Difference Index）、綠藍比值指數（Green-Blue Ratio Index）、可見光抗大氣指數（Atmospherically Resistant Vegetation Index）、綠紅比值指數（Green-Red Ratio Index）、有效葉面積指數（Specific Leaf Area Vegetation Index）、歸一化水分指數（Normalized Different Moisture Index）、歸一化紅外指數（Normalized Difference Infrared Index）等10種植被指數作為遙感特征因子。

2.5 紋理指數提取

Haralick（1979）定義了8種常用的紋理特征，即：均值、方差、協同性、對比度、相異性、熵、二階矩、相關性［6］。

3 模型算法

3.1 XGBoost

XGBoost算法是由陳天奇等人［7］提出來的，與一般的梯度boost方法相比，XGBoost對目標函數進行了二階泰勒展開，并在訓練過程中利用二階導數加快了模型的收斂速度。其核心思想就是不斷地添加樹，不斷地進行特征分裂來生長一棵樹，每次添加一個樹，其實是學習一個新函數f（x），去擬合上次預測的殘差。最后只需要將每棵樹對應的分數加起來就是該樣本的預測值。

3.2 MARS

多元自適應回歸樣條（Multivariate Adaptive Regression Splines，MARS）是一種數據分析方法。該方法以樣條函數的張量積作為基函數，分為前向過程、后向剪枝過程與模型選取三個步驟。其優勢在于能夠處理數據量大、維度高的數據，而且計算快捷、模型精確。

3.3 隨機森林

隨機森林算法是通過集成學習的思想將多棵決策樹集成的一種算法，它的基礎組成單元是決策樹，而隨機森林可以看作是若干棵決策樹的集成。基本組成單元采用CART算法［8］，其本質屬于機器學習的一大分支——集成學習（Ensemble Learning）方法，其工作原理是生成多個分類器或者模型，各自獨立地學習和做出預測。

4 模型精度評價

本研究選取交叉驗證（Cross Validation，CV）的方式來檢查模型的精度，采用決定系數（R-squared，R2）和均方根誤差（Root Mean Square Error，RMSE）進行評價。

4.1 Sentinel-2B

采用三種機器學習方法建立生態水反演模型，并通過十倍交叉驗證法進行精度驗證。其中XGBoost模型決定系數R2為0.45，均方根誤差RMSE為54；MARS模型決定系數R2為0.37，均方根誤差RMSE為57；隨機森林模型決定系數R2為0.40，均方根誤差RMSE為56。

4.2 Landsat8 OLI

采用三種機器學習方法建立生態水反演模型，并通過十倍交叉驗證法進行精度驗證。其中XGBoost模型決定系數R2為0.38，均方根誤差RMSE為57；MARS模型決定系數R2為0.35，均方根誤差RMSE為59；隨機森林模型決定系數R2為0.32，均方根誤差RMSE為60。

4.3 聯合反演

采用三種機器學習方法建立生態水反演模型，并通過十倍交叉驗證法進行精度驗證。其中XGBoost模型決定系數R2為0.57，均方根誤差RMSE為48；MARS模型決定系數R2為0.46，均方根誤差RMSE為53；隨機森林決定系數R2為0.43，均方根誤差RMSE為55。

4.4 模型評價

根據分析結果得出，通過Sentinel-2B影像和Landsat8 OLI影像分別提取的54個樣本點的各類特征波段、植被指數和紋理指數指標進行重要性分析，選取相對重要性和節點純度高的指標帶入XGBoost模型、MARS和RF三種模型進行相關性計算，發現提取指標聯合反演的相關性最高，Sentinel-2B影像提取的指標單獨反演的相關性其次，Landsat8 OLI影像提取的指標單獨反演的相關性最低。同類影像提取指標中XGBoost反演的相關性最高，MARS反演的相關性其次，RF反演的相關性最低。

5 植被生態水定量反演

利用森林生物量方程計算了54個樣地中的不同樹種的地上生物量，通過不同樹種地上生物量和含水率數據計算得到54個樣地植被生態水含量。對植被生態水含量反演方法進行研究，選擇特征波段、植被指數和紋理指數進行植被生態水含量反演，利用所建立的XGBoost模型植被生態水反演模型進行回歸分析，對研究區2019年12月份的植被生態水含量進行反演，反演結果如圖2所示。

圖2 研究區生態水儲量分布

根據本研究區的土地利用類型分類結果，研究區植被生態水儲量主要分布在針葉林、闊葉林、混交林和灌木林等四種森林類型，其中針葉林生態水儲量為351.94×104Mg，平均值210.33 Mg/hm2；闊葉林生態水儲量132.52×104Mg，平均值127.88 Mg/hm2；針闊混交林生態水儲量14.48×104Mg，平均值為132.52 Mg/hm2；灌木林生態水儲量為25.02×104Mg，平均值為150.54 Mg/hm2。研究區植被生態水儲量為523.97×104Mg，平均值為175.53 Mg/hm2。由于本研究基于的是光學遙感影像的定量反演，不具備LIDAR和INSAR等影像數據的穿透性，無法獲取枯落物和土壤層的生態水含量。同時基于生物量公式及含水率進行計算，因此研究區植被生態水儲量主要代表各森林類型的枝、干、葉等器官的地上植被生態水含量。

6 結論及展望

6.1 結 論

（1）根據特征選擇結果，在同一Sentinel-2B變量三種模型下的重要性分析中，紋理指數指標占比最多，其次為特征波段，植被指數占比最低。在同一Landsat8 OLI變量三種模型下的重要性分析中，紋理指數指標占比最多，其次為植被指數。在同一聯合變量三種模型下的重要性分析中，紋理指數指標占比最多，其次為特征波段，植被指數占比最低。可以看出三種模型下紋理指數指標的重要性占比均最多，是九寨溝自然保護區及川西高原生態水定量反演的重要參數指標。

（2）根據模型反演結果，同類影像提取指標中XGBoost模型反演的相關性最高，指標構建中聯合反演的相關性最高。在XGBoost模型中，結合Landsat8 OLI和Sentinel-2B影像作為預測變量提供了最好的生態水儲量估計。通過不同傳感器數據集和建模算法的協同作用，為基于遙感的九寨溝自然保護區及川西高原生態水定量反演提供了一種新的方法。

（3）根據植被生態水含量定量反演結果，按照不同植被類型生態水儲量總量占比，從大到小依次為針葉林、闊葉林、灌木林和針闊混交林。其中針葉林分布面積最大；生態水儲量平均值最大；生態水儲量總量最大，占比高達67.17%。

6.2 展 望

（1）對于植被生態水的時序監測，應該以不同時間（季節/年份）、不同植被類型、固定樣地的野外數據作分析研究，既簡化了野外工作，又方便建立研究區植被生態水數據庫。

（2）分析地形因子對植被生態水分布的影響。經野外采集數據時發現不同海拔、坡度條件下，植被生長狀況和植被類型都不完全一致，說明海拔、坡度可能是影響植被生態水含量的因素，具體相關性有待進一步研究。

（3）本研究中植被生態水反演采用了光學影像提取地面反演參量建立模型，未采用高光譜影像、SAR影像及其結合的研究。由于高光學影像對植被表層信息捕獲能力更強，SAR影像的后向散射信息還包含了枝干之間更為復雜的回波信號，可以作為植被垂直空間信息的補充。后續可以采用高光譜影像、SAR影像開展相關性反演工作，進一步提升植被生態水反演相關性及精度。