基于Landsat 的高山松地上生物量動態變化估測模型構建

廖 易 張加龍 鮑 瑞 許冬凡

（西南林業大學林學院，云南 昆明 650233）

森林生物量是森林生態系統最基本的數量特征[1]，它是森林生態系統長期代謝過程的積累，在一定程度上體現了森林的固碳能力。區域森林生物量的精確估測對該地區的森林資源經營管理有重要指示作用，對區域森林生態系統的研究和碳匯評價有重要意義[2]。森林地上生物量（AGB）是森林生物量的重要組成部分。當前，遙感技術已被廣泛運用于森林AGB 的估測，成為森林AGB估測的主要方法[3]。

林分的變化會引起遙感圖像光譜特征的變化，而光譜特征的變化對林分生物量值的改變有一定的指示作用。與光學影像相比，雖然激光雷達在捕捉森林結構方面有較大優勢，但其采集成本高，計算要求高，難以獲取大面積的數據[4]。此外，長時間序列的激光雷達數據有限，難以全面獲取森林生長變化[5]。目前大面積、長時間的森林生物量監測系統主要采用Landsat 時間序列等多光譜衛星和實地調查數據[6]。Landsat 存檔數據和預處理產品的開放，促進了其在森林監測活動中的研究[7]。當前，森林生物量估測方法應用較多的有線性回歸[8]、隨機森林[9]、支持向量機[10]等。隨著基于單期影像估測森林生物量方法的成熟，出現了大量的估測生物量動態的研究，其中大多采用（定期）變化量建模[11]。

變化量建模相對于靜態數據建模有一定的優勢。陸馳[12]通過研究證明定期變化量模型提升了預測精度。張加龍等[13-14]對比靜態和動態（即定期變化量）的森林生物量建模，也得出動態模型比靜態模型效果更好。國外學者Gómez 等[15]將變化過程分為狀態變量（靜態）和過程變量（動態）構建模型，結果表明過程變量比狀態變量對森林AGB 具有更強的預測能力。Puliti 等[16]使用雙時相（AGB 變化量）和單時數據估測AGB 變化量，結果表明雙時相方法可以得到更精確的AGB 變化值。Boisvenue 等[17]以樹木的年齡和截距擬合變化量模型，并結合遙感信息的變化估測生物量，結果表明該模型能夠更好的估計生物量及其變化。這些研究雖然都構建了動態變化模型，但其中大多只采用定期變化量建模，Gómez等雖然采用變化率進行研究，但并未比較定期變化量和變化率的效果，且這些研究的建模效果仍然有一定的提升空間。

故本研究擬基于1987—2017 年的森林資源連續清查樣地和Landsat TM/OLI 數據，采用多元線性回歸和隨機森林建模方法，在年平均變化量和定期變化量建模研究的基礎上加入變化率進行香格里拉市高山松AGB 模型構建，以期運用不同類型的遙感因子變化量，來提高高山松森林AGB 動態變化遙感估測的精度，為森林AGB 的動態估測研究提供技術支持。

1 研究區概況

研究區香格里拉市位于云南省西北部，地處北緯26°52′～28°52′，東經99°20′～100°19′，總面積為11 613 km2。該區域年平均氣溫為6.3 ℃，年平均降水量為651.1 mm，日照充足、太陽輻射強，屬寒溫帶季風氣候。境內水系發達，水資源豐富，全境河流均屬金沙江水系。區域內地形起伏較大，最高點海拔5 545 m，最低點海拔1 503 m，平均海拔3 459 m。香格里拉地處“三江并流”的中心地帶，有豐富的自然和生物資源，且風光獨特，是西南地區重要自然資源儲存地。境內有包括褐土、紅壤和高山草甸土等在內的11 個森林土壤類型，森林覆蓋率較高，達到75%[18]，高山松（Pinus densata）、云南松（Pinus yunnanensis）與云冷杉（Picea asperata）為優勢樹種[19]。

2 材料與方法

2.1 數據來源與處理

2.1.1 樣地數據

研究所采用的樣地數據為1987—2017 年共7 期的國家森林資源連續清查固定樣地數據，各年份樣地數量分別為19、22、23、16、16、17、23，共計136 塊，其中包含部分復測樣地。各年份樣地胸徑、樹高統計結果見表1。樣地大小為28.28 m × 28.28 m，面積為0.08 hm2，測樹因子包含胸徑、樹高和株數等。

表1 樣地測樹因子統計情況Table 1 Statistical table of tree measuring factors in sample plots

2.1.2 遙感數據

研究獲取的Landsat 影像數據共7 期21 景，從美國地質勘探局（http://glovis.usgs.gov/）網站獲取。所有影像的獲取年份與該地區連清高山松固定樣地調查數據相對應，影像經過輻射定標和FLAASH 大氣校正預處理。因研究區域內地形起伏較大，使用坡度匹配法進行地形校正。參照已有的SPOT-5 影像（北京1954 坐標系）對選取的Landsat 影像進行幾何校正，誤差控制在1 個像元內。在ENVI 5.1 軟件中對各部分影像進行拼接，然后使用香格里拉市行政邊界數據進行裁剪，分別得到各年份的香格里拉市影像。

2.2 遙感因子提取

提取5 種類型573 種遙感因子用于篩選自變量，詳見表2。因子提取在ENVI 5.1 中進行，提取因子時以樣地中心點為基礎，落入Landsat 影像的對應位置并提取因子值。

表2 遙感因子Table 2 Remote sensing factors

2.3 高山松樣地地上生物量計算

根據單木生物量模型[20]來計算高山松單木生物量（W），該模型為項目組研究成果，是利用實測的116 株高山松，選用胸徑（D）、樹高（H）作為模型的自變量建立，模型見公式（1）。使用連清數據中每塊高山松樣地的平均胸徑、平均樹高代入該單木模型計算平均單木生物量，然后乘以樣地高山松株數得到樣地生物量的近似值。

2.4 生物量和遙感因子變化量的計算

本次研究中，提出將定期變化量和變化率分別按照5 a 和10 a 進行研究。年平均和定期變化量計算時，剔除未進行重復觀測和變化量為負值的數據，變化率數據采用7 個年份的共同樣地計算。計算變化量后采用皮爾遜相關性法進行遙感因子篩選，每種變化量均選擇10 個相關性最強且顯著的遙感因子用于建模。計算得到的生物量變化量統計值見表3。

表3 地上生物量變化量統計表Table 3 Statistical table of variations of AGB

2.4.1 定期變化值

定期變化量是指從1987 年、1992 年、1997 年、2002 年、2007 年、2012 年、2017 年 的 樣 地 數 據中選取2 個年份，首先篩選出這2 個年份的重復觀測樣地，用較大年份的AGB 和遙感因子值減去較小年份的對應值即得到變化量值。本研究計算了5 a 和10 a 的定期變化量，計算公式如下：

式中： Δb為定期變化量，b為高山松AGB 或對應的遙感因子值，n和m為用來計算變化量的2 個不同年份。

2.4.2 年平均變化值

年平均變化量包括5、10、15、20、25、30 a變化的年平均值，將定期變化量值比上對應的變化年份即可得到年平均變化量。公式如下：

式 中： Δg為年平均變化量，Δb為定期變化量，n和m為用來計算變化量的2 個不同年份。使用該公式計算AGB 及對應的遙感因子變化值。

2.4.3 變化率值

當 Δx→0時，Δy/Δx有極限，則稱函數y=f(x)在點x0處可導，這個極限叫做f(x)在x0處的導數，這個導數即為該點處的瞬時變化率[21]，簡稱變化率，通常記為f′(x0)，表達式為：

將高山松AGB 取重復觀測樣地數據按年份進行排列，以年份和生物量分別作為橫縱坐標進行擬合并生成曲線，對曲線每個年份處求導得到此處導數作為變化率，即表示生物量在該年份處5 a 或10 a 的瞬時變化速率。本研究中得到變化率建模數據的過程中需進行2 次相關性篩選，第1 次篩選出20 個與原始樣地AGB 相關性最高的遙感因子，計算其與對應AGB 的變化率后，進行第2 次相關性篩選，得到10 個相關性最高的遙感因子變化率用于建模。

2.5 建模方法

本研究隨機抽取80%的數據采用多元線性回歸和隨機森林方法建模，剩余20%數據用于檢驗。本研究中同一種變化類型中多元線性回歸（MLR）與隨機森林（RF）模型輸入的建模因子相同。

2.5.1 多元線性回歸

當多個自變量與因變量之間是線性關系時，可以進行多元線性回歸分析。多元線性回歸模型的基本公式為：

式中：β0為 常數項，β1～βn表示回歸系數，?為隨機誤差，xn表 示自變量，在本研究中代表遙感因子，y表示因變量，在本研究中代表AGB。使用方差膨脹因子（VIF）[22]來克服變量之間的共線性問題。

2.5.2 隨機森林（RF）

隨機森林[23]是一個集成學習模型，其中含有多個決策樹，常被應用于多領域的數據分類和非參數回歸[24]。本研究中隨機森林模型算法的實現基于Python 語言“Sklearn”包中提供的“Random Forest Regressor”算法。建模過程中需對模型參數進行調參，以得到較好的模型結果。主要進行調整的參數為：模型的最大迭代次數n_estimators，決策樹的最大深度max_depth，葉節點的最小樣本數min_samples_leaf，內部節點再劃分所需最小樣本數min_samples_split。調參時根據樣本量因素，min_samples_leaf 和min_samples_split設為默認值。

2.6 模型精度評價方法

本研究采用的精度評價指標為決定系數（R2）、均方根誤差（RMSE）及相對均方根誤差（rRMSE），計算公式如下：

式中：yi為真實值，為模型回歸值，為真實值均值，n為樣本個數。

3 結果與分析

3.1 建模因子分析

本研究通過皮爾遜相關性分析篩選用于建模的遙感因子，由于因子較多，每種變化量選出10 個因子用于建模。5 種變化量最終選出的遙感因子見表4。由表4 可知，采用RF 方法進行5 種類型的變化量建模，選入的50 個（每種變化量選入10 個）建模因子中48 個為紋理因子，其中SK 有16 個，出現頻率高，HO、EN、DI 各有6個，ME、SM 各有4 個、VA 有3 個、CO 有2 個、CC 只有1 個，另有植被指數因子ND53 有2 個，表明紋理因子在生物量建模中具有較強的敏感性。

表4 5 種變化量的建模因子Table 4 Modeling factors of 5 variables

3.1.1 多元線性回歸模型

采用篩選出的遙感因子變化量進行MLR 建模，得到各類型的MLR 模型結果見表5。3 類5 種變化量的MLR 模型均通過方差檢驗。由表5可知，最終選入模型的因子均為紋理因子，這表明紋理因子對MLR 的變化率建模較為敏感。

表5 5 種變化量的MLR 模型Table 5 MLR model with 5 type variations

3.1.2 隨機森林模型

調用Python 語言“Sklearn”包中提供的“Random Forest Regressor”算法建模，得到RF的最佳模型參數見表6。隨機森林會將輸入的因子根據重要性進行排序，并輸出各因子的建模貢獻度值。本研究中5 種變化量類型的遙感因子貢獻度（%）見圖1。圖中5 種變化量建模中貢獻度最高的因子分別為R17B5ME、R11B7SK、ND53、R19B7HO 及R11B5SK，其中R19B7HO 在所有因子中貢獻值最高，達50.26%。

表6 RF 模型最佳參數Table 6 Optimal parameters of RF model

圖1 5 種變化量類型的RF 建模因子貢獻度Fig. 1 Contribution degree of RF modeling factors for 5 variation types

3.2 模型精度評價

3 類變化量數據分別采用MLR 與RF 方法構建AGB 變化估測模型，得出結果的對比見表7，表中的模型編號及其變化量類型與表5、表6 相對應。

表7 3 類5 種變化量的MLR 與RF 建模精度評價Table 7 Accuracy evaluation of MLR and RF modeling based on 3 types and 5 kinds

由表7 對比所有模型的擬合效果，發現模型3、2、1、5、4、8、7、6、10、9 的R2依次增大，其中MLR 模型4 的R2最高（0.465），RF 模型9 的R2最高（0.956）。采用RF 建模，模型9 的擬合RMSE 最小，為0.664；采用MLR 建模，模型3 的RMSE 最小，為2.318。MLR 模型擬合的rRMSE 值最小為48.742%（模型4）；RF 模型擬合的rRMSE 值最小為35.883%（模型9）。以上結果表明模型9 的擬合效果最好，即RF 模型的5 a 變化率模型擬合效果最好。MLR 和RF 模型中，各變化量的建模擬合效果從劣到較優依次為年平均變化量、10 a 定期變化量、5 a 定期變化量、10 a 和5 a 變化率。表7 中，模型3 和模型8 的RMSE 在同種模型（MLR 或RF）中較小，這是變化的時間尺度造成的，模型3 和8 均為年平均變化量建模，表示地上生物量一年的平均變化量，而其他4 種變化量表示的是地上生物量至少5 年的變化，因此其變化數值相對較小，最終導致模型的RMSE 也較小。

模型預測檢驗結果與擬合效果對應，即模型3、2、1、5、4、8、7、6、10、9 的檢驗效果依次變優。RF 模型8、7、6、10、9 的RMSE 分別為2.082、20.547、17.244、4.422、2.285，MLR 模型3、2、1、5、4 的RMSE 分 別為2.342、21.491、21.277、5.301、4.865。預測效果的rRMSE 最小為42.718%（模型9）。表明5 a 變化率RF 建模的預測效果最好。

經對比，同一變化量類型的RF 建模效果優于MLR 建模效果，其中同種變化量類型MLR 與RF 模型的R2的最大差值為0.695，最小差值為0.491，且5 種變化量中5 a 變化率的RF 和MLR建模效果均為最優。綜上，認為其符合一定規律，即3 類5 種變化量中5 a 變化率建模效果最好，年平均變化量最差，具體優劣排序為年平均變化量＜10 a 定期變化量＜5 a 定期變化量＜10 a 變化率＜5 a 變化率。

3.3 地上生物量反演

根據模型精度評價結果，采用最優的5 a 變化率RF 模型反演，先反演出各年份的變化率，再通過變化率公式反推得出AGB，結合已有的高山松分布面積范圍數據[14]得出香格里拉市各年份變化率和AGB 空間分布圖（圖2、圖3），高山松AGB 統計見表8。由表8 可知，1987—1997 年香格里拉市高山松AGB 總量逐漸增長，最高為929.02 萬t，1997—2017 年AGB 總量逐漸減少，最低為872.25 萬t。

圖2 各年份變化率分布圖Fig. 2 Distribution of change rate in each year

圖3 各年份高山松AGB 分布圖Fig. 3 AGB distribution of P. densata in each year

表8 各年份高山松AGB 反演結果Table 8 Result of the estimating AGB of P. densata in each year

4 結論與討論

本研究以1987—2017 年每5 年1 期的香格里拉市高山松固定樣地為基礎，提取對應樣地的遙感因子，經過處理得到3 類5 種變化量數據，分別為定期變化量（5、10 a 定期變化量）、年平均變化量、變化率（5、10 a 變化率），分別以這5 種變化量進行多元線性回歸（MLR）和隨機森林（RF）建模并比較各模型的精度。結果表明：

1）在同種變化量類型中，RF 模型的效果均優于對應的MLR 模型，其中RF 模型和MLR 模型的R2最大差值達0.705，RMSE 最多降低了11.236，可見非參數模型在遙感森林AGB 的動態變化估測中相較于參數模型具有一定的優勢。

2）所有類型的變化量中，無論采用MLR 方法還是RF 方法建模，變化率（5、10 a 變化率）的模型效果均優于年平均變化量和定期變化量，其中5 a 變化率的RF 建模效果最好，其MLR 和RF 模 型 的R2分 別 為0.454 和0.956，RMSE 分 別為2.543 和0.664。采用變化率建模能夠提高模型的精度。

3）紋理因子在本研究的5 類遙感因子中為出現頻次最多的因子。本研究共篩選出相關性高的50 個建模因子，其中48 個為紋理因子，因此其具有最強的敏感性。在各種9 種紋理因子中，偏度（SK）出現次數最多，最具敏感性，相關性（CC）出現次數最少，最不具敏感性。

本研究在前人的基礎上進一步提出采用年平均變化量與瞬時變化率建模，對比了3 種變化量（變化率、年平均變化量、定期變化量）模型的建模效果，得出結果表明瞬時變化率模型相較于年平均變化量和定期變化量模型在估測精度上有一定的提升，可以為地上生物量的精確估測提供一定的參考。

本研究樣地數據時間跨度大，調查技術手段的改變也可能造成調查結果的不一致性，最終可能會對模型的評價結果產生一定的偏差。變化量模型精度較高，有利于區域AGB 的動態變化研究，在今后的研究中可以進一步挖掘更多類型的變化量用以提升模型精度。研究區香格里拉市地處高原地區且地形起伏明顯，海拔、坡度、坡向等地形因素可能對AGB 有較大影響，但本研究暫未對該方面進行探究，今后的研究可以從此方面著手以探究地形對AGB 變化量模型的影響。