結合冠層密度的森林凈初級生產(chǎn)力遙感估測

李 陶，李明陽，錢春花

(南京林業(yè)大學林學院，江蘇 南京 210037)

森林凈初級生產(chǎn)力(net primary productivity，NPP)是指森林在單位面積、單位時間內所積累的有機物數(shù)量，表現(xiàn)為光合作用固定的有機碳中扣除植物本身呼吸消耗的部分[1]。NPP作為地表碳循環(huán)的組成部分，能夠直接反映自然環(huán)境條件下森林的生產(chǎn)能力，是評價森林生態(tài)系統(tǒng)可持續(xù)發(fā)展的重要指標之一，是生態(tài)系統(tǒng)碳匯和調節(jié)生態(tài)過程判定的主要因子，在全球變化及碳平衡中有著重要的作用[2-3]。森林凈初級生產(chǎn)力的估算研究有助于了解全球碳循環(huán)，為合理利用自然資源、修復受損森林等提供科學依據(jù)。

利用遙感數(shù)據(jù)可以獲得大區(qū)域尺度內的森林覆蓋狀況、森林光譜特征等信息，同時這些信息與森林生產(chǎn)力密切相關，光學傳感器和主動傳感器已經(jīng)被廣泛應用于森林凈初級生產(chǎn)力估測中[4]。因Landsat系列遙感衛(wèi)星像元的空間分辨率為30 m×30 m，像元大小與森林經(jīng)營單位(小班)的最小區(qū)劃面積接近，經(jīng)常被用來進行區(qū)域森林生產(chǎn)力遙感估測。區(qū)域尺度森林生產(chǎn)力估測中結合地面調查數(shù)據(jù)的遙感估測已成為重要的手段之一，但是在過去的研究中發(fā)現(xiàn)僅僅采用 Landsat系列遙感圖像構建模型估測森林NPP，存在光飽和因子問題，為了更加準確地構建森林NPP遙感估測模型，本研究引入一個新的變量——森林冠層密度(canopy density)。森林生產(chǎn)力與年齡、植被生長狀況有關，而采用傳統(tǒng)方法提取的遙感特征變量卻不包含這些信息，森林冠層密度是從空中觀測到的森林植被與投影面積的比例，在數(shù)值上與郁閉度相同，是反映林分密度的指標，是決定林分結構的重要因子之一，可以用來表示光照、水分、溫度等環(huán)境因子通過林冠進入林分的再分布狀況[5]。

本研究以廣東省重點林區(qū)韶關市為研究對象，結合森林冠層密度采用Landsat-8 OLI數(shù)據(jù)和隨機森林、多元線性回歸、人工神經(jīng)網(wǎng)絡、K-最近鄰分類法等4種模型，進行森林凈初級生產(chǎn)力特征變量選取、參數(shù)建模、精度評價、空間制圖，以期能在區(qū)域尺度上通過遙感反演較為準確地預估森林NPP，通過加入森林冠層密度這一變量以提高NPP遙感估測模型的估測精度。

1 材料與方法

1.1 研究區(qū)概況

韶關市位于廣東省北部，北接湖南，東鄰江西，東南面、南面和西面分別與廣東省河源、惠州、廣州及清遠等市接壤(112°53′～114°45′E，23°53′～25°31′N)，全境直線距離東西跨度為186.3 km，南北跨度為173.4 km，全市土地面積1.84萬km2。韶關地形以山地丘陵為主，河谷盆地分布其中，平原、臺地面積約占20%。地勢北高南低，北部石坑崆海拔1 902 m，為廣東第一高峰，南部最低處海拔僅35 m。河流主要為珠江水系北江流域，北江以湞江為干流，主要支流有武江、墨江、錦江、滃江、南水。

韶關屬于中亞熱帶濕潤季風氣候區(qū)，氣候宜人。年平均氣溫21 ℃，冬季各地氣溫自北向南遞增，夏季各地氣溫較接近。雨量充沛，年均降雨1 700 mm，3—8月為雨季，9—2月為旱季，冬季北部有雪。Landsat-8 OLI遙感圖像受云量過多的限制，為保證研究質量，本研究僅選取2017年10月下旬的條帶號122、行編號43、云量0.7%的一景圖像進行研究，研究區(qū)域面積占韶關市總面積的94.5%。

韶關是全國重點林區(qū)，為廣東省用材林、水源林和重點毛竹林基地，被譽為華南生物基因庫和珠江三角洲的生態(tài)屏障。全市林業(yè)用地面積144.8萬hm2，森林覆蓋率為73.8%，林木綠化率為74.6%，森林蓄積量8 917萬m3。省級以上生態(tài)公益林面積占林業(yè)用地面積比達45.6%，其中，國家級公益林占46.3%，省級公益林占53.7%。韶關市森林闊葉混交林和闊葉純林居多，闊葉純林主要有櫟類(Quercusspp.)、桉樹(Eucalyptusspp.)、樟樹(Cinnamomumcamphora)等樹種；其次是以馬尾松(Pinusmassoniana)、杉木(Cunninghamialanceolata)等為代表的針葉純林、針闊混交林以及以毛竹(Phyllostachysedulis)等為代表的竹林等。

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

本研究主要數(shù)據(jù)包括：①自地理空間數(shù)據(jù)云網(wǎng)站(http://www. gscloud. cn/)下載的廣東省韶關市2017年10月Landsat-8 OLI遙感影像，條帶號122，行編號43，云量0.7%，多光譜波段空間分辨率為30 m×30 m；②地理空間數(shù)據(jù)云網(wǎng)站下載的韶關地區(qū)數(shù)字高程模型(DEM)，空間分辨率為30 m×30 m；③韶關市2017年357塊森林資源連續(xù)清查固定樣地數(shù)據(jù)，樣地間距為6 km×8 km，每塊樣地面積為0.06 hm2，樣地屬性包括地形、地貌、土壤名稱、優(yōu)勢樹種、平均胸徑、平均樹高、齡組、郁閉度等60多個調查因子。

1.2.1 影像數(shù)據(jù)預處理

1)Landsat-8 OLI影像預處理。使用ENVI 5.3軟件對Landsat-8 OLI影像進行預處理。為消除大氣和光照等因素對地物反射的影響，首先對影像進行輻射定標，然后使用ENVI軟件下的FLAASH模塊進行快速大氣校正。研究區(qū)內地形以山地丘陵為主，南北海拔差異較大，遙感影像受傳感器方位與太陽高度、方位等影響，造成因陰坡、陽坡接受的照度不同而引起亮度值出現(xiàn)差異，因此采用基于余弦校正模型的C校正算法對遙感圖像進行地形校正，消除由于地形起伏而引起的影像輻射亮度值的變化，使影像更好地反映地物的光譜特征。

2)特征變量的提取。圖像預處理結束后，根據(jù)肖興威[6]的研究，森林生產(chǎn)力與林分、土壤、地形等狀況密切相關，所以使用ENVI軟件提取歸一化植被指數(shù)(NDVI)、比值植被指數(shù)(RVI)、差值植被指數(shù)(DVI)、增強型植被指數(shù)(EVI)、綠度植被指數(shù)(GVI)、垂直植被指數(shù)(PVI)等6種植被指數(shù)以及葉面積指數(shù)(LAI)。對遙感圖像進行纓帽變換，提取前3個分別反映土壤巖石、植被及土壤和植被中水分信息的分量?；贒EM提取地形因子——坡度作為變量之一。依據(jù)Li等[7]的研究結果，本研究提取紋理特征時窗口大小選擇3×3，灰度量化等級為64，計算對比度、相異性、平均值、均一性、角二階矩、熵、偏度、熵相關性等。預測變量類型、名稱、數(shù)量見表1。

表1 NPP遙感估測模型備選自變量Table 1 Independent variables set for model of remote sensing estimation of NPP

1.2.2 樣地森林凈初級生產(chǎn)力計算

森林凈初級生產(chǎn)力由群落生長量和年凋落量組成。根據(jù)2017年森林資源連續(xù)清查數(shù)據(jù)中樣地林木蓄積量數(shù)據(jù)及余超等[8]提供的中國不同森林類型生物量與蓄積量、生物量與群落生長量和年凋落量之間的函數(shù)關系，以蓄積量為基礎計算不同森林類型的生物量，依據(jù)生物量分別計算群落生長量與凋落量，求得每塊樣地的森林凈初級生產(chǎn)力(NPP)[9]。

1.2.3 森林冠層密度的計算

引用森林冠層密度制圖器(FCD)進行森林冠層密度的提取。根據(jù)Joshi等[10]的研究可知，森林冠層密度的提取有人工神經(jīng)網(wǎng)絡(ANN)、多元線性回歸技術(MLR)、森林冠層密度制圖器(FCD)、最大似然分類(MLC)4種方法。在4種方法中，F(xiàn)CD在3個東南亞國家獲得的森林冠層密度準確率平均值為92%[11]。FCD模型以森林生長現(xiàn)象為基礎，能夠檢測森林冠層密度時間上的變化。FCD模型基于Landsat-8提取植被指數(shù)(IV)、裸土指數(shù)(IB)、陰影指數(shù)(IS)3個指數(shù)來計算森林冠層密度，繪制森林冠層密度圖。與NDVI相比，IV對植被數(shù)量的反應更為敏感，IB隨地表裸露程度的增加而增大，IS隨森林密度的增加而增大。3個指數(shù)及FCD(式中記為FCD)計算公式如下：

IV=[(B5+1)(65 536-B4)(B5-B4)]1/3;

(1)

(2)

IS=[(65 536-B2)(65 536-B3)(65 536-B4)]1/3;

(3)

FCD=(VD×ISS+1)1/2-1。

(4)

式中：B2—B6分別表示藍、綠、紅、近紅外和短波紅外波段的亮度值;當(B5-B4) 小于0時，IV=0；VD表示植被密度值，由IV和IB利用主成分分析合成得到；ISS表示尺度陰影指數(shù)，通過IS的線性變換得到[12]。

1.3 遙感估測模型

隨機森林(random forest，RF)是一種組合分類器，它利用一種重抽樣(bootstrap)方法從原始樣本中抽取多個樣本，對每個重抽樣的樣本進行決策樹建模，然后將決策樹組合，通過投票得到最終的預測結果[13]。大量理論與實際研究證明隨機森林預測準確率較高，對于異常值和噪音有很好的容忍性，由于隨機森林生成決策樹的樣本是隨機選擇的，能夠減少過擬合現(xiàn)象[14]。隨機森林模型通過R語言randomForest數(shù)據(jù)包中randomForest函數(shù)執(zhí)行，在使用時需要調整建立的決策樹數(shù)量(ntree)和決策樹分裂時抽取的變量個數(shù)(mtry)兩個參數(shù)。mtry值為構建隨機森林回歸模型時自變量個數(shù)的1/3，當自變量個數(shù)小于3時，mtry值默認為1[15]。

多元線性回歸(multiple linear regression，MLR)有兩個或兩個以上自變量。多元線性回歸模型以植被指數(shù)、紋理特征、原始波段、森林冠層密度等66個因子作為自變量，以2017年韶關樣地NPP作為因變量，建立遙感數(shù)據(jù)與樣地實測NPP之間的相關關系，估算森林NPP。自變量的篩選選擇逐步回歸，其基本思想是將變量逐個引入模型，每引入一個新變量進行一次F檢驗，入選的模型變量再逐一進行t檢驗，當原來引入的變量由于后面變量的引入變得不再顯著時，則將其刪除，直至既沒有顯著的變量入選回歸方程、也沒有不顯著的變量從回歸方程中剔除為止，此時回歸模型中變量都與因變量呈顯著相關。逐步回歸通過SPSS軟件執(zhí)行，步進標準使用F概率，進入與刪除分別設置為0.05與0.1。

人工神經(jīng)網(wǎng)絡(artificial neural network，ANN)以數(shù)學模型模擬神經(jīng)元活動，是基于模仿大腦神經(jīng)網(wǎng)絡結構和功能而建立的一種信息處理系統(tǒng)。目前，人工神經(jīng)網(wǎng)絡模型中多層感知器類型的反向傳播網(wǎng)絡(BP網(wǎng)絡)已成為應用最為廣泛的模型之一，本研究使用3層結構的BP神經(jīng)網(wǎng)絡模型，包括輸入層、隱藏層和輸出層[16]。BP神經(jīng)網(wǎng)絡模型用R語言neuralnet函數(shù)包進行構建，首先利用系統(tǒng)默認節(jié)點數(shù)和隱含層數(shù)建立模型，根據(jù)結果誤差情況，進一步增加隱含層數(shù)目，以此提高模型精度。

K-最近鄰分類法(K-nearest neighbor，K-NN)是基于類比學習，通過給定的檢驗元組與和它相似的訓練元組進行比較來學習[17]。根據(jù)給定的測試樣本，計算測試樣本與訓練集中每個對象的距離，圈定距離最近的K個訓練對象作為最近鄰，這K個最近鄰中選出此次數(shù)最多的類標號作為測試樣本的類標號值。K-最近鄰分類模型用R語言class包中的knn函數(shù)實現(xiàn)，依據(jù)陳爾學等[18]的研究，本研究設置K=10。

1.4 模型特征變量選擇

本研究基礎變量的優(yōu)化選擇統(tǒng)計檢驗方法，即相關性分析，是指對兩個或多個具備相關性的變量元素進行分析，從而衡量兩個變量因素的相關密切程度。本研究共有66個備選變量，使用SPSS進行Pearson相關性分析，選擇相關性顯著的特征變量加入模型構建。經(jīng)Pearson相關性分析，共有10個特征變量與NPP具有顯著相關性(表2)。

表2 特征因子與NPP相關性Table 2 The correlation between characteristic factors and NPP

1.5 模型精度驗證與評價

本研究使用十折交叉驗證進行模型精度驗證。該方法將數(shù)據(jù)集分成10份，輪流將其中9份作為訓練數(shù)據(jù)，1份作為測試數(shù)據(jù)進行試驗，10次結果均值將作為對算法精度的估計，本研究進行10次十折交叉驗證求取平均值，即為該模型最終精度。

遙感估測模型評價指標有很多，如決定系數(shù)(R2)[19]、均方根誤差(RMSE)、平均絕對誤差(MAE]、總體相對誤差(RS)、平均相對誤差(E1)、總體相對誤差絕對值(E2)、預估精度(P)等。采用R2、RMSE和MAE對模型精度進行評價。R2反映變量的全部變異能通過回歸關系被自變量解釋的比例，其區(qū)間通常為(0，1)，R2的值越大，表示預測值與實際值之間相關性越強。RMAE是預測值與真實值的誤差平方根的均值，預測值與實際值之間的RMSE越小，表示模型預測的效果越好。MAE是絕對誤差的平均值，可以更好地反映預測誤差的實際情況，當預測值與真實值完全吻合時MAE等于0，誤差越大，該值越大。

本研究模型擬合主要是通過計算模型的決定系數(shù)(R2)、均方根誤差[RMSE，式中記為σ(RMSE)]、平均絕對誤差[MAE，式中記為σ(MAE)]進行精度評價。

(5)

(6)

(7)

2 結果與分析

2.1 不同森林類型NPP比較分析

不同森林類型的NPP對比見圖1。從圖1可看出：2017年研究區(qū)內風景林的NPP最高[16.95 t/(hm2·a)]，原因為風景林主要位于研究區(qū)人口密集、經(jīng)濟發(fā)達的城鎮(zhèn)及周邊地區(qū)，經(jīng)營者投資力度大，經(jīng)營強度高；自然保護林的NPP為16.91 t/(hm2·a)，排名第2，與遠離城鎮(zhèn)、保護力度大、人為干擾少存在密切聯(lián)系；而水土保持林的NPP最低[10.64 t/(hm2·a)]，原因在于水土保持林大多地處海拔較高、坡度較陡、立地條件較差的地段；短輪伐期用材林的NPP[14.49 t/(hm2·a)]高于一般用材林[13.69 t/(hm2·a)]，原因在于短輪伐期的樹種以速生豐產(chǎn)的桉樹為主，而一般用材林則以輪伐期較長、生長較慢的杉木、馬尾松為主。

圖1 不同森林類型的NPP比較Fig.1 Comparison of NPP of different forest category

不同林型的NPP比較見圖2。從圖2可看出，2017年研究區(qū)內竹林的NPP最高[26.35 t/(hm2·a)]，這與竹類生長迅速的特性相吻合；闊葉林的NPP[12.39 t/(hm2·a)]高于針葉林的[11.41 t/(hm2·a)]；闊葉混交林的NPP[17.71 t/(hm2·a)]高于針葉混交林的[17.50 t/(hm2·a)]。針闊混交林的NPP[12.56 t/(hm2·a)]低于針葉混交林的原因是針闊混交林以天然次生林為主，其經(jīng)濟效益差，經(jīng)營強度低，林木生長緩慢。而針葉混交林以人工營造的水源涵養(yǎng)林和一般用材林為主，經(jīng)營強度較高，受人為不良干擾較少，林木生長狀況較好，生長速度較快。

圖2 不同林型的NPP比較Fig.2 Comparison of NPP of different forest type

2.2 森林冠層密度的精度驗證

在357塊樣地中，篩選出林地中郁閉度不為0的233塊樣地，對FCD模型估測的森林冠層密度進行驗證，郁閉度與森林冠層密度的相關系數(shù)為0.87，RMSE為0.12，由此可見FCD模型估測精度較高，可用此模型對森林冠層密度進行估測。

依據(jù)FCD模型計算得到研究區(qū)森林冠層密度圖(圖3)。從圖3可以看出，F(xiàn)CD的空間分布與海拔、人口密度等地形條件及經(jīng)濟發(fā)展水平存在密切的關系：研究區(qū)中部、西部及東北部，地形以河谷盆地為主，坡度平緩，城鎮(zhèn)密集，經(jīng)濟較為發(fā)達，人口密度大，森林覆蓋率低，森林冠層密度較低；研究區(qū)北部、南部及西南部，海拔較高，坡度較陡，人口密度較低，森林覆蓋率高，人為干擾少，林木生長狀況良好，森林冠層密度較高。

圖3 研究區(qū)森林冠層密度圖Fig.3 Forest canopy density map of the study area

2.3 模型精度評價

模型精度評價分為是否加入森林冠層密度兩種情形，將表2中的特征變量，分別代入隨機森林、人工神經(jīng)網(wǎng)絡、多元線性回歸、K-最近鄰分類法4種模型中，建立NPP遙感估測模型，使用十折交叉驗證法進行驗證。具體遙感預測精度評價結果見表3。

表3顯示，在4種遙感估測模型中，無論是否加入森林冠層密度，隨機森林算法預測精度最高，其次是人工神經(jīng)網(wǎng)絡、多元線性回歸，K-最近鄰分類法預測精度最低。隨機森林是進行兩次取樣，首先是有放回的隨機抽樣得到訓練樣本的采樣集，然后從所有特征中隨機選擇特征，并選擇最佳分割特征作為節(jié)點構建分類回歸樹，這使得最終建立的模型有較強的泛化性。研究區(qū)海拔差異大，地形復雜，因變量與大多數(shù)特征變量線性相關程度較低，因此在加入森林冠層密度之前多元線性回歸在4種估測模型中估測精度最低。在多元線性回歸模型中，加入冠層密度之后，由于冠層密度與NPP相關系數(shù)較高，并且冠層密度與NPP直接相關的年齡、胸徑、樹高相關系數(shù)較高，較大幅度提升了多元線性回歸模型的預測精度。人工神經(jīng)網(wǎng)絡相對于傳統(tǒng)的機器學習算法來說，需要更多的數(shù)據(jù)，本研究僅有樣地數(shù)據(jù)357塊，因此預測精度較低。K-最近鄰分類法預測精度最差的原因是分類預測時一般選擇多數(shù)表決法,當樣本容量較小時較容易產(chǎn)生誤分。從表3還可以看出，森林冠層密度加入4種遙感估測模型后預測精度明顯提升，這是由于森林冠層密度與年齡、林分密度、林分結構等有關，而以上因素均會影響森林NPP的大小。

表3 加入森林冠層密度前后模型精度對比分析Table 3 Comparative analysis of model accuracy before and after adding forest canopy density

2.4 NPP空間制圖

通過以上建模方式的對比，最終在4種遙感估測模型中選擇最優(yōu)的結合森林冠層密度的隨機森林模型進行研究區(qū)內NPP制圖(圖4)。

從圖4可以看出，NPP較高的區(qū)域主要分布在海拔較高的北部、南部、西南部，NPP較低的林分主要分布在中部、西部、東北部海拔較低的丘陵、河谷盆地。海拔較高的山地人口密度小，人為干擾活動少，森林植被能夠較好地進行自然生長；海拔低的丘陵河谷地帶人口密度大，森林覆蓋率低，人為干擾活動多，林分質量差。由此可以看出，NPP的空間分布趨勢與研究區(qū)的地貌特征、社會經(jīng)濟狀況吻合度較高。

圖4 研究區(qū)森林凈生產(chǎn)力分級圖Fig.4 Net primary productivity grading map of study area

計算結果顯示，研究區(qū)內森林NPP平均值為10.689 t/(hm2·a)，標準差為7.389 t/(hm2·a)。為進一步了解研究區(qū)內各地區(qū)NPP情況，用平均值±標準差的方法，將NPP劃分為低[<3 t/(hm2·a)]、中[≥3～11 t/(hm2·a)]、較高[≥11～18 t/(hm2·a)]和高[≥18 t/(hm2·a)]4個等級。結果表明，研究區(qū)內森林NPP以低級(面積占比39.22%)、較高級(面積占比24.09%)為主，其次是高級(面積占比19.61%)，最低的是中級(面積占比17.08%)。

研究區(qū)內森林NPP平均值為10.689 t/(hm2·a)，其他研究中廣州市NPP平均值為6.95 t/(hm2·a)[20]，這與姜春[21]的研究中珠江三角區(qū)森林NPP年均值要明顯低于粵北地區(qū)的研究結果相吻合。同期中國森林NPP平均值為 9.50 t/(hm2·a)[22]，廣西的NPP為6.62 t/(hm2·a)[19]，福建的NPP為5.87 t/(hm2·a)[23]，江西的NPP為5.23 t/(hm2·a)[24]。研究區(qū)森林NPP稍高于全國森林NPP，這與韶關為全國重點林區(qū)有關；明顯高于鄰近省份，原因是研究區(qū)內森林以中幼林為主，生長速度較快，同時韶關為全國重點林區(qū)，林地質量較好，森林經(jīng)營強度較高，對于森林保護較為重視。

預測變量的選取包括原始波段、植被指數(shù)、紋理特征、地形特征、森林冠層密度等多方面特征因子，但是4種遙感估測模型均出現(xiàn)預測精度小于擬合精度的現(xiàn)象，說明反演過程中依然存在過擬合現(xiàn)象。這可能與地形、林種、林種結構等有關，如何將林種、樹種分類特征引入特征變量，從而提高遙感估測模型精度，需要進一步研究。

3 結 論

1)在區(qū)域森林NPP遙感估測中，加入冠層密度這一特征變量，可以較大幅度提高模型估測精度，彌補了特征變量中因缺失年齡、植被生長狀況等變量而帶來的估測誤差較大的缺陷。

2)在基于Landsat 8 OLI提取的眾多變量中，紅光波段(B4)、紋理特征、歸一化植被指數(shù)(NDVI)、比值植被指數(shù)(RVI)、葉面積指數(shù)(LAI)、地形因子、森林冠層密度等因子，在區(qū)域森林NPP遙感估測時應予以優(yōu)先考慮。

3)在4種估測模型中，隨機森林模型精度最高，其次是多元線性回歸模型、人工神經(jīng)網(wǎng)絡模型，K-最近鄰分類模型預測精度最差。因此，在地形較為復雜的南方集體林可以選擇隨機森林作為NPP遙感估測模型。