應(yīng)用Landsat8 OLI和GBRT對高山松地上生物量的估測1)

張加龍 胥輝 陸馳

(西南林業(yè)大學(xué)，昆明，650224)

森林生物量是陸地生態(tài)系統(tǒng)碳循環(huán)過程中最基本的數(shù)量特征[1]和主要參數(shù)[2]，其精確估算能大大減小碳儲(chǔ)量估測的不確定性[3]，也為森林可持續(xù)經(jīng)營與發(fā)展、全球氣候變化研究提供科學(xué)依據(jù)。日益發(fā)展的遙感技術(shù)具有快速、實(shí)時(shí)、時(shí)空監(jiān)測的優(yōu)勢[4]，使得遙感監(jiān)測手段成為獲取森林地上生物量的主要途徑。

利用遙感數(shù)據(jù)估測生物量存在樣地大小、數(shù)據(jù)類型和建模方法3個(gè)大的問題，通過文獻(xiàn)總結(jié)可以發(fā)現(xiàn)：生物量一般不直接測量，而是采用樣地調(diào)查，基于樹高和胸徑變量的異速生長方程來計(jì)算；應(yīng)用頻率最高的傳感器依次是機(jī)載LiDAR、星載多光譜、機(jī)載或星載雷達(dá)(RADAR)系統(tǒng)；建模方法依次為線性模型、地統(tǒng)計(jì)平滑法、分類與回歸樹算法、近鄰法、支持向量機(jī)方法、隨機(jī)森林、高斯過程和其他數(shù)據(jù)挖掘方法[3,5]。機(jī)載雷達(dá)獲取成本較高，而Landsat覆蓋了全球范圍中等分辨率最長記錄(1972至今)的對地觀測數(shù)據(jù)，美國國家航空航天局(簡稱NASA)于2009年免費(fèi)開放存檔數(shù)據(jù)[6]，使得獲取高時(shí)空的地表森林參數(shù)成為可能，也成為生物量估測重要的遙感數(shù)據(jù)源。

采用遙感數(shù)據(jù)估測生物量，仍然存在遙感因子選取、生物量建模方法[7]、遙感數(shù)據(jù)飽和[8]等眾多不確定性問題。梯度提升回歸樹(GBRT)是一種非參數(shù)建模方法，是Fridedman在1999年提出的一種組合模型，它的基本思想是通過構(gòu)建M個(gè)弱分類器，經(jīng)過多次迭代最終組合而成一個(gè)強(qiáng)分類器[9-10]。GBRT在生物量建模方面的研究很少，類似的采用隨機(jī)梯度提升(SGB)方法進(jìn)行遙感生物量建模，發(fā)現(xiàn)梯度提升方法要優(yōu)于多元自適應(yīng)回歸樣條(MARS)和Cubist非參數(shù)方法[11]。

滇西北地區(qū)為中國最大的世界自然和文化遺產(chǎn)保護(hù)地，擁有豐富的森林資源[12]、較高的生物量和生產(chǎn)力,這些森林資源不僅為當(dāng)?shù)氐纳鐣?huì)經(jīng)濟(jì)發(fā)展提供了直接的物質(zhì)基礎(chǔ),而且在維護(hù)地區(qū)生態(tài)和碳平衡中起著巨大作用[13]，也為遙感估測生物量提供了較好的地域條件。本研究使用Landsat8 OLI數(shù)據(jù)，以香格里拉市優(yōu)勢樹種高山松為研究對象，通過外業(yè)樣地調(diào)查并進(jìn)行生物量計(jì)算的基礎(chǔ)上，提取對應(yīng)的多種遙感因子，篩選重要性因子，對比參數(shù)和非參數(shù)的建模方法，以尋找較為理想的建模方法，為區(qū)域生物量遙感估測提供技術(shù)和理論更新。

1 研究區(qū)概況

研究區(qū)香格里拉市(如圖1)隸屬迪慶藏族自治州，位于云南西北部，其地理范圍為北緯26°52′～28°52′，東經(jīng)99°20′～100°19′，總面積11 613 km2。香格里拉地形起伏較大，海拔高差4 042 m(最低海拔1 503 m)，平均海拔3 459 m。其森林覆蓋率較高，達(dá)到75%[14]。主要植被類型為寒溫性針葉林，優(yōu)勢樹種依次為云冷杉(Piceaasperata)、高山松(Pinusdensata)、云南松(Pinusyunnanensis)、高山櫟(Quercussemicarpifolia)等。

圖1 研究區(qū)概況及外業(yè)樣地分布圖

2 研究方法

2.1 外業(yè)樣地調(diào)查及地上生物量計(jì)算

先后于2015年11月和2016年3月共調(diào)查58塊30 m×30 m樣地，在每個(gè)鄉(xiāng)鎮(zhèn)高山松純林區(qū)域隨機(jī)布設(shè)，且任意兩樣地間相隔3 km以上。采用羅盤儀進(jìn)行布設(shè)樣地邊界后，利用靜態(tài)差分GPS(精度為5 cm)進(jìn)行準(zhǔn)確定位，對胸徑大于5 cm的高山松進(jìn)行每木檢尺[15]，記錄每木樹高和胸徑。樣木的調(diào)查、取樣、測定、單木生物量構(gòu)建過程詳見文獻(xiàn)[16]，最后得到的單木生物量模型見式1。

W=0.073×D1.739×H0.880。

(1)

式中：D為胸徑；H為樹高。應(yīng)用樣地每木檢尺胸徑和樹高數(shù)據(jù)，計(jì)算單木生物量，匯總得到樣地地上生物量。

2.2 Landsat數(shù)據(jù)及處理

Landsat8 OLI來源于http://glovis.usgs.gov/，三景數(shù)據(jù)編號(hào)分別為LC81310412015354LGN00、LC81320402015313LGN00、LC81320412015313LGN00，均為2015年11—12月份成像，且云量小于2%。對三景數(shù)據(jù)進(jìn)行輻射定標(biāo)，之后采用FLAASH方法進(jìn)行大氣校正。以研究區(qū)SPOT-5影像數(shù)據(jù)為參考數(shù)據(jù)進(jìn)行幾何校正，選擇了100個(gè)地面控制點(diǎn)和二階多項(xiàng)式來校正Landsat影像，誤差控制在一個(gè)像元內(nèi)。采用坡度匹配模型[17]進(jìn)行地形校正，影像校正(第二次校正)后陰坡和陽坡的反射率平均值差異均較小。

提取原始波段、簡單比值植被指數(shù)、主成分變換(PCA)、紋理因子，具體如下。為了方便引用文獻(xiàn)中的公式計(jì)算，將Landsat 8第1—第7個(gè)波段(Band)按順序分別對應(yīng)命名：海岸帶(C)、藍(lán)(B1)、綠(B2)、紅(B3)、近紅(B4)、短波紅外(B5)、短波紅外2(B7)。

①原始單波段：C、B1-B5、B7。

②簡單比值植被指數(shù)[18-20]：B4/B3、B4/B2、B5/B3、B5/B4、B5/B7、B7/B3、B3/Albedo、B4×B3/B7，其中Albedo為B1-B5、B7共6個(gè)波段的和。

③植被指數(shù)[11,21]：NDVI、ND32、ND54、ND53、ND57、ND452、DVI。

④信息增強(qiáng)[22-23]：主成分變換3個(gè)，圖像變換特征(VIS123)、Albedo、MID57。

⑤紋理信息[24]：均值(ME)、方差(VA)、均一性(HO)、反差(CO)、相異(DI)、熵(EN)、角二階矩(SM)、相關(guān)性(CC)8個(gè)，窗口大小有5×5和9×9兩種(用R5和R9分別代表5×5和9×9窗口，B1、B2、…、B7分別代表波段，紋理因子采用字母縮寫)。

為了更為準(zhǔn)確地將地面樣地值與遙感像元值對應(yīng)，首先利用像元的邊界將樣地劃分開，如圖2所示。然后利用每個(gè)地塊重心位置的遙感因子值乘以所在地塊的面積，再將4個(gè)值相加除以總面積即為此樣地對應(yīng)的遙感因子值。樣地與像元的匹配誤差可能對于估算結(jié)果的影響不大，但此方法可以更為客觀地將樣地的生物量值對應(yīng)到影像像元上。

圖2 樣地及Landsat像元匹配示意圖

2.3 因子篩選方法

分析因子的重要性，采用散點(diǎn)圖矩陣和皮爾遜相關(guān)系數(shù)來篩選建模因子。利用散點(diǎn)圖矩陣來繪制各自變量間的散點(diǎn)圖，可以快速發(fā)現(xiàn)生物量和遙感因子多個(gè)變量間的主要相關(guān)性。通過R語言下的“colSelection”函數(shù)選擇變量因子，“featurePlot”來繪制生物量和變量之間的散點(diǎn)圖矩陣。同時(shí)，采用皮爾遜相關(guān)系數(shù)度量生物量和遙感因子之間的相關(guān)性。

2.4 建模方法

為了更好地闡明梯度提升回歸樹在森林地上生物量估測表現(xiàn)出的優(yōu)勢，與多元線性回歸、偏最小二乘法、隨機(jī)森林建模方法和結(jié)果進(jìn)行了對比分析。

2.4.1 多元線性回歸(MLR)

本研究中，變量選入的顯著性水平設(shè)定為p≤0.05，變量剔除的水平設(shè)置為p≥0.1。同時(shí)，為了客服變量之間的共線性問題，使用方差膨脹因子(VIF)[25]來評(píng)價(jià)；當(dāng)VIF值大于10時(shí)，對應(yīng)的自變量將舍去。最后，有8個(gè)自變量選入到MLR模型中。

2.4.2 偏最小二乘法(PLSR)

利用SIMCA-P分析軟件和SPSS Statistics Python Integration Plug-in插件，構(gòu)建高山松地上生物量和遙感因子的偏最小二乘模型。

2.4.3 隨機(jī)森林(RF)

本研究使用R平臺(tái)下的語言包RF進(jìn)行高山松地上生物量和遙感因子的建模分析。該語言包提供了因子重要性篩選，回歸樹與回歸誤差的判定，生物量預(yù)測及交叉驗(yàn)證[26]。

2.4.4 梯度提升回歸樹

2.4.4.1模型基本思想

梯度提升決策樹(GBRT)[9-10]通過構(gòu)建M個(gè)弱分類器，經(jīng)過多次迭代最終組合而成一個(gè)強(qiáng)分類器。它在原有的模型基礎(chǔ)上進(jìn)行迭代改進(jìn)，使得下一次生成的模型比上次的模型有更小的誤差，在殘差減小的梯度方向上建立新的組合模型。決策樹的類型有回歸和分類，GBRT可以累加所有決策樹的輸出為最終輸出，而分類樹的特性決定無法累加，因此在GBDT中，分類器基本都采用回歸樹。

(2)

(3)

對于m=1、2、…、M，將上式進(jìn)行轉(zhuǎn)換：

(4)

fm(x)=-ρmgm(x)。

(5)

其中，gm(xi)、Fm-1(x)和乘數(shù)ρm用下式表示

(6)

(7)

(8)

Fm(x)是之前取得Fm-1(x)模型的梯度方向下降的損失函數(shù)。使用最小二乘法得

(9)

進(jìn)而得到

(10)

通過迭代M次，得到P的參數(shù)集合；通過合并，最后得到模型

Fm(x)=Fm-1(x)+ρmhm(x,αm)。

(11)

2.4.4.2規(guī)則化

在進(jìn)行預(yù)測時(shí)，擬合值過于接近訓(xùn)練值時(shí)會(huì)適得其反，規(guī)則化方法通過限制擬合過程來防止“過度擬合”，通過增加一個(gè)參數(shù)v將式11替換為

Fm(x)=Fm-1(x)+vρmhm(x)，0

(12)

式中：ρm表示每次迭代的步長；v稱為學(xué)習(xí)速率，用以設(shè)置每次迭代學(xué)習(xí)的速度；hm(x)表示決策樹，又稱為弱分類器，其大小由決策樹深度S表示，它用以修正當(dāng)前模型Fm-1(x)(由前m-1次迭代產(chǎn)生)的誤差來產(chǎn)生更好的模型Fm(x)。

2.4.4.3模型中的基本參數(shù)及實(shí)現(xiàn)

在模型學(xué)習(xí)前，需要設(shè)置迭代次數(shù)(M)、決策樹深度(S)和學(xué)習(xí)速率(V)這些參數(shù)。M應(yīng)使模型誤差在收斂的情況下，越小越好；較小的v值(通常取v<0.1)可以減少模型的測試誤差；而S為了避免過擬合情況的出現(xiàn)通常不應(yīng)過大。本研究在Anaconda3平臺(tái)下編程實(shí)現(xiàn)。

①讀取數(shù)據(jù)。

②將自變量數(shù)據(jù)標(biāo)準(zhǔn)化。

③模型參數(shù)選擇。需要確定以下模型參數(shù)：弱分類器(決策樹)數(shù)量，學(xué)習(xí)速率，決策樹深度。模型重要參數(shù)的確定往往需要遵循一定的經(jīng)驗(yàn)，比如M通常在預(yù)測值收斂的情況下越小越好，v通常小于0.1，S通常不應(yīng)大于15。依據(jù)這些經(jīng)驗(yàn)，設(shè)置這些參數(shù)的范圍分別為‘決策樹數(shù)量’:[300、500、600、800、1 000、2 000]、‘學(xué)習(xí)速率’:[0.005、0.010、0.050、0.001]、‘決策樹深度’:[2、3、4、…、12]，然后使用格網(wǎng)搜索的方式進(jìn)行最優(yōu)模型參數(shù)的確定。格網(wǎng)搜索窮盡比較每種參數(shù)值組合下模型學(xué)習(xí)的效果，以得到最優(yōu)參數(shù)的模型。

④特征變量選擇。利用上一步得到的最優(yōu)模型參數(shù)進(jìn)行變量重要性計(jì)算。根據(jù)重要性排序圖，取重要性大的作為建模因子。

⑤用選出的特征因子建模。把選出的特征因子及生物量數(shù)據(jù)整理成新的數(shù)據(jù)集，并重復(fù)進(jìn)行①、②、③步驟，得到新的模型參數(shù)。

⑥建立GBRT生物量估測模型。利用第5步的模型參數(shù)和15個(gè)變量因子建立估測模型，得到擬合結(jié)果。

⑦預(yù)測。利用上步建立的GBRT模型，讀取檢驗(yàn)數(shù)據(jù)集，輸出預(yù)測結(jié)果。

2.5 精度評(píng)價(jià)方法

①相對均方根誤差(rRMSE)。它是預(yù)測值與真實(shí)值的相對偏離程度的度量。其定義為

(13)

②決定系數(shù)。實(shí)測值與預(yù)測值間的決定系數(shù)(R2)，也稱為擬合優(yōu)度：

(14)

③預(yù)估精度(P)。它能夠反映出模型的平均預(yù)估能力，定義為

(15)

3 結(jié)果與分析

3.1 高山松樣地地上生物量

表1 高山松樣地地上生物量統(tǒng)計(jì)

3.2 建模因子比較

MLR模型建模因子中有2個(gè)紋理因子，1個(gè)植被指數(shù)因子，紋理因子占了2/3；PLSR模型有7個(gè)紋理因子，占7/15；RF模型中紋理因子有11個(gè)，占11/15；GBRT模型有12個(gè)紋理因子，占4/5。紋理因子相對于其他因子來說是比較重要的因子，即紋理信息對生物量有重要影響。R5B4EN在4個(gè)模型的建模因子中均有出現(xiàn)，說明R5B4EN因子對生物量具有最為重要的影響。從對生物量有重要影響的因子中觀察發(fā)現(xiàn)，由第四波段衍生出的因子有4個(gè)，說明影像第四波段(近紅外波段)的信息對生物量具有最為重要的影響。在6個(gè)重要影響力的因子中，有2個(gè)因子為紋理信息中的相關(guān)性(CC)，說明相關(guān)性因子對生物量有重要影響(見表2)。

3.3 GBRT建模

采用GBRT進(jìn)行建模，求得最佳模型參數(shù)為n_estimators=500,subsample=0.5,learning_rate=0.05,max_depth=6,min_samples_leaf=2。模型偏差圖如圖3a。對于訓(xùn)練數(shù)據(jù)集，當(dāng)?shù)螖?shù)為1～200次時(shí)，偏差迅速減??；當(dāng)大于200時(shí)，偏差降低減緩；當(dāng)達(dá)到600時(shí)，基本趨于0。對于測試數(shù)據(jù)集偏差，當(dāng)?shù)螖?shù)在0～100時(shí)，值迅速減小，當(dāng)大于150時(shí)，偏差降低減緩；當(dāng)達(dá)到200次時(shí)，基本趨于穩(wěn)定。

注:均值(ME)、方差(VA)、均一性(HO)、反差(CO)、相異(DI)、熵(EN)、角二階矩(SM)、相關(guān)性(CC)。R5和R9分別代表5×5和9×9窗口，B1、B2、…、B7分別代表波段。

采用上述重要性因子再次進(jìn)行特征變量選擇，求得最佳模型參數(shù)為n_estimators=600,subsample=0.5,learning_rate=0.05,max_depth=7,min_samples_leaf=4。模型偏差圖如圖3b，對比參數(shù)參與的模型偏差圖(圖5a)，迭代次數(shù)減少到500次偏差趨于穩(wěn)定，且測試數(shù)據(jù)集的偏差減小到200以下，第二次選擇的特征變量建模降低了模型誤差。

圖3 采用GBRT建模偏差與迭代次數(shù)之間的關(guān)系

利用上述最佳模型參數(shù)和重要性變量進(jìn)行建模，得到擬合結(jié)果。再利用上述建模結(jié)果，采用剩余的20%數(shù)據(jù)集進(jìn)行預(yù)測，得到預(yù)測結(jié)果。

3.4 精度檢驗(yàn)

為了便于對4種回歸模型的比較，評(píng)價(jià)結(jié)果在表3中列出。

通過比較可發(fā)現(xiàn)，梯度提升決策樹的R2最大；相對均方根誤差最小，大約為MLR和PLSR的1/4，RF的1/2，因此GBRT在這4個(gè)模型中擬合效果最好。

在模型預(yù)測效果方面，GBRT具有最小的相對均方根誤差(32.50%)，說明該模型的相對偏差在4個(gè)模型中最??；GBRT具有最高的預(yù)估精度(P=73.88%)，表明該模型的預(yù)測能力最好。

表3 4種建模方法的精度評(píng)價(jià)結(jié)果

4 結(jié)論與討論

依據(jù)遙感對生物量建模的不確定因素還有很多[29-30]。樣地的設(shè)置數(shù)量滿足一定的統(tǒng)計(jì)抽樣要求[31]，其空間分布應(yīng)盡量隨機(jī)分布；樣地的類型應(yīng)盡量覆蓋各種立地類型，達(dá)到此要求在研究區(qū)實(shí)際較困難。

非參數(shù)建模在本研究中優(yōu)于參數(shù)建模，但是非參數(shù)方法無法給出具體的模型形式。梯度提升回歸樹在本次森林地上生物量建模方面表現(xiàn)出較好模擬精度，但參數(shù)選擇運(yùn)行時(shí)間較長，迭代次數(shù)、模型評(píng)價(jià)及適用性方面還需要進(jìn)一步研究。紋理因子相關(guān)性最強(qiáng)，今后可以繼續(xù)研究各因子對于生物量估測的貢獻(xiàn)及精度影響，找出更為有效的反映樹高和胸徑的指標(biāo)因子，為生物量建模提供更多支持。Landsat8 OLI近紅外波段的紋理因子方差和相關(guān)性的重要性較高，回歸分析時(shí)，當(dāng)建模因子過多或過少模型精度都較低，選擇4～6個(gè)因子時(shí)模型擬合精度較高。

梯度提升回歸樹在生物量估測方面的研究還較少，它在本研究中的相對均方根誤差、R2、P值都優(yōu)于隨機(jī)森林方法。在非參數(shù)方法中，梯度提升回歸樹可作為高山松遙感生物量建模的另一新的可靠途徑。