基于Landsat 的高山松地上生物量動態(tài)變化估測模型構(gòu)建

廖 易 張加龍 鮑 瑞 許冬凡

（西南林業(yè)大學(xué)林學(xué)院，云南 昆明 650233）

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

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

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

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

1 研究區(qū)概況

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

2 材料與方法

2.1 數(shù)據(jù)來源與處理

2.1.1 樣地數(shù)據(jù)

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

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

2.1.2 遙感數(shù)據(jù)

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

2.2 遙感因子提取

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

表2 遙感因子Table 2 Remote sensing factors

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

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

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

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

表3 地上生物量變化量統(tǒng)計表Table 3 Statistical table of variations of AGB

2.4.1 定期變化值

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

式中： Δb為定期變化量，b為高山松AGB 或?qū)?yīng)的遙感因子值，n和m為用來計算變化量的2 個不同年份。

2.4.2 年平均變化值

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

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

2.4.3 變化率值

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

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

2.5 建模方法

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

2.5.1 多元線性回歸

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

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

2.5.2 隨機(jī)森林（RF）

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

2.6 模型精度評價方法

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

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

3 結(jié)果與分析

3.1 建模因子分析

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

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

3.1.1 多元線性回歸模型

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

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

3.1.2 隨機(jī)森林模型

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

表6 RF 模型最佳參數(shù)Table 6 Optimal parameters of RF model

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

3.2 模型精度評價

3 類變化量數(shù)據(jù)分別采用MLR 與RF 方法構(gòu)建AGB 變化估測模型，得出結(jié)果的對比見表7，表中的模型編號及其變化量類型與表5、表6 相對應(yīng)。

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

由表7 對比所有模型的擬合效果，發(fā)現(xiàn)模型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）。以上結(jié)果表明模型9 的擬合效果最好，即RF 模型的5 a 變化率模型擬合效果最好。MLR 和RF 模型中，各變化量的建模擬合效果從劣到較優(yōu)依次為年平均變化量、10 a 定期變化量、5 a 定期變化量、10 a 和5 a 變化率。表7 中，模型3 和模型8 的RMSE 在同種模型（MLR 或RF）中較小，這是變化的時間尺度造成的，模型3 和8 均為年平均變化量建模，表示地上生物量一年的平均變化量，而其他4 種變化量表示的是地上生物量至少5 年的變化，因此其變化數(shù)值相對較小，最終導(dǎo)致模型的RMSE 也較小。

模型預(yù)測檢驗結(jié)果與擬合效果對應(yīng)，即模型3、2、1、5、4、8、7、6、10、9 的檢驗效果依次變優(yōu)。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。預(yù)測效果的rRMSE 最小為42.718%（模型9）。表明5 a 變化率RF 建模的預(yù)測效果最好。

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

3.3 地上生物量反演

根據(jù)模型精度評價結(jié)果，采用最優(yōu)的5 a 變化率RF 模型反演，先反演出各年份的變化率，再通過變化率公式反推得出AGB，結(jié)合已有的高山松分布面積范圍數(shù)據(jù)[14]得出香格里拉市各年份變化率和AGB 空間分布圖（圖2、圖3），高山松AGB 統(tǒng)計見表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 反演結(jié)果Table 8 Result of the estimating AGB of P. densata in each year

4 結(jié)論與討論

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

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

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

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

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

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