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

廖 易 張加龍 鮑 瑞 許冬凡

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

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

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

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

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

1 研究區(qū)概況

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

2 材料與方法

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

2.1.1 樣地?cái)?shù)據(jù)

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

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

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

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

2.2 遙感因子提取

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

表2 遙感因子Table 2 Remote sensing factors

2.3 高山松樣地地上生物量計(jì)算

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

2.4 生物量和遙感因子變化量的計(jì)算

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

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

2.4.1 定期變化值

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

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

2.4.2 年平均變化值

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

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

2.4.3 變化率值

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

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

2.5 建模方法

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

2.5.1 多元線性回歸

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

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

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

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

2.6 模型精度評(píng)價(jià)方法

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

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

3 結(jié)果與分析

3.1 建模因子分析

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

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

3.1.1 多元線性回歸模型

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

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

3.1.2 隨機(jī)森林模型

調(diào)用Python 語(yǔ)言“Sklearn”包中提供的“Random Forest Regressor”算法建模，得到RF的最佳模型參數(shù)見表6。隨機(jī)森林會(huì)將輸入的因子根據(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 模型精度評(píng)價(jià)

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

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

由表7 對(duì)比所有模型的擬合效果，發(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）中較小，這是變化的時(shí)間尺度造成的，模型3 和8 均為年平均變化量建模，表示地上生物量一年的平均變化量，而其他4 種變化量表示的是地上生物量至少5 年的變化，因此其變化數(shù)值相對(duì)較小，最終導(dǎo)致模型的RMSE 也較小。

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

經(jīng)對(duì)比，同一變化量類型的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)劣排序?yàn)槟昶骄兓浚?0 a 定期變化量＜5 a 定期變化量＜10 a 變化率＜5 a 變化率。

3.3 地上生物量反演

根據(jù)模型精度評(píng)價(jià)結(jié)果，采用最優(yōu)的5 a 變化率RF 模型反演，先反演出各年份的變化率，再通過變化率公式反推得出AGB，結(jié)合已有的高山松分布面積范圍數(shù)據(jù)[14]得出香格里拉市各年份變化率和AGB 空間分布圖（圖2、圖3），高山松AGB 統(tǒng)計(jì)見表8。由表8 可知，1987—1997 年香格里拉市高山松AGB 總量逐漸增長(zhǎng)，最高為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ǔ)，提取對(duì)應(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)于對(duì)應(yīng)的MLR 模型，其中RF 模型和MLR 模型的R2最大差值達(dá)0.705，RMSE 最多降低了11.236，可見非參數(shù)模型在遙感森林AGB 的動(dòng)態(tài)變化估測(cè)中相較于參數(shù)模型具有一定的優(yōu)勢(shì)。

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 個(gè)建模因子，其中48 個(gè)為紋理因子，因此其具有最強(qiáng)的敏感性。在各種9 種紋理因子中，偏度（SK）出現(xiàn)次數(shù)最多，最具敏感性，相關(guān)性（CC）出現(xiàn)次數(shù)最少，最不具敏感性。

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

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