汪風(fēng)華，文 敏，任藍(lán)翔，胡中岳

(1.桃江縣林業(yè)局，湖南 桃江 413400；2.國家林業(yè)和草原局中南調(diào)查規(guī)劃設(shè)計(jì)院，湖南 長沙 410014)

碳儲(chǔ)量是量化森林生態(tài)系統(tǒng)所需的重要變量之一。已有較多研究方法用于估測局部、區(qū)域以及全球范圍的森林碳儲(chǔ)量，比如野外測量、遙感方法和基于過程的生態(tài)系統(tǒng)模型等[1]。遙感技術(shù)的快速發(fā)展，為大面積、長周期的森林碳密度動(dòng)態(tài)研究提供了新的途徑和方法，遙感影像可以表征大面積土地表面特征。但是，由于光學(xué)和雷達(dá)影像的數(shù)據(jù)飽和，光譜和空間分辨率的限制，以及碳儲(chǔ)量和光譜變量之間的復(fù)雜關(guān)系，使碳儲(chǔ)量估測精度較低，尤其是當(dāng)碳儲(chǔ)量值過高或過低時(shí)[2]。

目前利用遙感估測森林碳密度，一般是通過遙感數(shù)據(jù)和森林地面樣地建立數(shù)學(xué)模型實(shí)現(xiàn)不同空間尺度的估算[3]。常用方法中，經(jīng)驗(yàn)?zāi)Ｐ秃唵巍⒈阌谟?jì)算，但需要大量觀測樣點(diǎn)數(shù)據(jù)，受樹種和區(qū)域背景影響較大[4]；過程模型能揭示碳儲(chǔ)量的形成機(jī)理，但對(duì)觀測數(shù)據(jù)要求較高[5]；半經(jīng)驗(yàn)半物理模型簡單直觀，但缺乏植被生理機(jī)制方面的嚴(yán)密依據(jù)，很少用于森林碳儲(chǔ)量的反演研究[6]；非參數(shù)估計(jì)模型忽略隨機(jī)擾動(dòng)項(xiàng)[7]，精度較高。本研究選取2013年桃江縣森林資源一類調(diào)查數(shù)據(jù)和Landsat 8遙感影像，構(gòu)建徑向基函數(shù)(Radical Basis Function，RBF)神經(jīng)網(wǎng)絡(luò)模型，并與多元逐步回歸模型(Linear Stepwise Regression，LSR)和偏最小二乘回歸模型(Partial Least Square Regression，PLSR)進(jìn)行對(duì)比，估測桃江縣森林碳儲(chǔ)量，探討可以改善其過高和過低估計(jì)的建模方法。

1 研究區(qū)概況

桃江縣位于湖南省益陽市，地處湘中偏北，資江下游，地理位置為111°36′—112°19′E，28°13′—28°41′N。該縣國土總面積2068 km2，平均海拔200m，屬亞熱帶季風(fēng)氣候，年均氣溫16.6 ℃，森林覆蓋率62.98%，常見樹種為杉木、馬尾松、青岡等。

2 材料與方法

2.1 遙感影像及信息處理

研究選用2013年7月的Landsat 8 OLI 遙感影像，軌道號(hào)為124/040。用ENVI 5.3軟件對(duì)其進(jìn)行輻射校正和幾何校正，誤差小于半個(gè)像元。根據(jù)桃江縣行政區(qū)劃圖對(duì)影像進(jìn)行裁剪，得到桃江縣Landsat 8遙感影像圖。植被信息可以通過植被指數(shù)來表征。本研究選用Landsat 8影像多光譜波段以及6種植被指數(shù)及其衍生因子，作為碳儲(chǔ)量遙感反演的可選光譜因子[8]，包括歸一化植被指數(shù)(NDVI)、土壤修正植被指數(shù)(SAVII，I=0.1，0.25，0.3，0.5)、大氣抗阻植被指數(shù)(ARVI)、增強(qiáng)型植被指數(shù)(EVI)、差值植被指數(shù)(DVI)、比值植被指數(shù)(RVI)，從而減少坡度和坡向的影響[9]。同時(shí)，用研究區(qū)DEM計(jì)算高程、坡向、坡度等地理因子，并提取3×3窗口下的紋理因子[10-11]。用ArcGIS 10.4軟件獲取變量因子值，用R軟件corr.test()函數(shù)計(jì)算其與碳儲(chǔ)量的相關(guān)系數(shù)。遙感變量作為模型的獨(dú)立變量可能會(huì)彼此顯著相關(guān)，導(dǎo)致信息重復(fù)并干擾模型性能。因此，需要對(duì)變量因子進(jìn)行篩選。

2.2 樣地信息獲取與處理

研究以2013年桃江縣一類調(diào)查數(shù)據(jù)為實(shí)測碳儲(chǔ)量的數(shù)據(jù)來源。在桃江縣的88個(gè)樣地中，保留48個(gè)地類為林地的參與模型構(gòu)建。目前對(duì)碳儲(chǔ)量的估算主要根據(jù)生物量與含碳率的乘積，因此獲取各樹種生物量非常關(guān)鍵。竹林、灌木、經(jīng)濟(jì)林和混交林的含碳率均為0.5000，其他主要樹種為：杉木0.5201，馬尾松0.4596，硬闊葉類0.4834，軟闊葉類0.4956，桉樹0.5253[12]。

喬木林生物量采用李?？萚12]2010年開創(chuàng)的模型(表1)，經(jīng)濟(jì)林和灌木林生物量分別為23.7 t·hm-2、19.76 t·hm-2，混交林按比例計(jì)。

2.3 多元逐步回歸模型

多元逐步回歸模型根據(jù)自變量因子是否同時(shí)符合兩個(gè)條件來篩選：一是對(duì)因變量有足夠影響力，二是因子之間相關(guān)性較低[13]。在不斷的引入和刪除自變量的過程中，既保證了多元逐步回歸模型不會(huì)漏選對(duì)因變量影響顯著的自變量因子，又可以避免多重共線性問題。

2.4 偏最小二乘回歸模型

偏最小二乘回歸模型可以從總體中篩選出具有最佳解釋作用的綜合自變量，可以解決變量之間的多重相關(guān)性[14]。首先計(jì)算交叉有效性和累計(jì)解釋量，確定潛在建模因子個(gè)數(shù)，然后用投影重要性指標(biāo)VIPj篩選出因子構(gòu)建模型(式1)。

(1)

表1 喬木樹種生物量計(jì)算方程Tab.1 Biomass regression equation and carbon content of different tree species樹種生物量計(jì)算方程杉木和其他杉類BS=0.073 429(D2H)0.862 62;BP=0.013 775(D2H)0.844 63;BB=0.000 482(D2H)1.233 14; BL=0.019 638(D2H)0.789 69;BT=BS+BP+BB+BL馬尾松和其他松類BT=0.071 556(D2H)0.857 209硬闊葉類軟闊葉類BT=0.049 550 2(D2H)0.952 453BT=0.049 550 2(D2H)0.952 453桉樹柏木BS=0.090 252 6D2.448 15; BB=0.004 916 3D2.817 79; BL=0.012 394D2.268 39; BT=BS+BB+BLBS=0.125 31(D2H)0.733; BB=0.137 403+0.012 887 D2H; BL=0.053 49+0.009 97D2H; BT=BS+BB+BL竹林BT=0.643 9D1.537 3 注:BS—樹干生物量;Bp—樹皮生物量;BB—樹枝生物量;BL—樹葉生物量;BT—地上部分總生物量;D—樣地平均胸徑;H—樣地平均樹高。

式中：n為自變量數(shù)，Ph為相關(guān)自變量主成分，R(Y，Ph)為因變量Y和相關(guān)自變量主成分Ph的相關(guān)系數(shù)，ωhj為自變量在主成分上的權(quán)重。

當(dāng)VIPj值大于1，說明該自變量因子對(duì)因變量有強(qiáng)解釋作用，對(duì)模型有重要貢獻(xiàn)；當(dāng)VIPj值小于1時(shí)，自變量的解釋作用弱，對(duì)模型的貢獻(xiàn)性小，VIPj<0.8的變量可以考慮剔除。

2.5 RBF神經(jīng)網(wǎng)絡(luò)模型

RBF神經(jīng)網(wǎng)絡(luò)模型可以有效避免局部極小問題，一般有三層：輸入層、具有非線性激活功能的隱藏層和線性輸出層[15]。本研究將通過K-means聚類得到的n個(gè)中心作為基函數(shù)的n個(gè)中心，隱藏層以高斯函數(shù)作為基函數(shù)，方差計(jì)算見式(2)，隱藏層到輸出層的權(quán)重可以用最小均方根誤差求得。該模型利用Matlab軟件的newrb()函數(shù)來實(shí)現(xiàn)。

(2)

式中：σi為方差，cmax為聚類中心的最大距離，n為隱含層的節(jié)點(diǎn)數(shù)。

2.6 精度驗(yàn)證

將48個(gè)樣地隨機(jī)分成兩部分：2/3(32個(gè))作為訓(xùn)練數(shù)據(jù)，1/3(16個(gè))作為驗(yàn)證數(shù)據(jù)，對(duì)3種模型得到的森林碳儲(chǔ)量值開展精度驗(yàn)證。模型效果的評(píng)估選取決定系數(shù)R2和相對(duì)均方根誤差RRMSE，計(jì)算公式分別見式(3)和式(4)。

(3)

(4)

3 結(jié)果與分析

3.1 相關(guān)性分析

研究用ArcGIS 10.4軟件提取Landsat 8的194個(gè)光譜、地理和紋理因子，作為建?？蛇x的自變量因子，用R軟件計(jì)算相關(guān)性(表2)。

表2 森林碳儲(chǔ)量與自變量的相關(guān)系數(shù)分析Tab.2 Analysis of correlation coefficients between forest carbon stock and independent variables因子相關(guān)系數(shù)因子相關(guān)系數(shù)因子相關(guān)系數(shù)因子相關(guān)系數(shù)SR34-0.632**SR460.511**SR260.415**DVI130.358**SAVI0250.617**NDVI0.509**EVI0.408**SR74-0.353**SAVI0350.523**SAVI050.501**DVI360.377**SR160.313**SAVI010.520**SR24-0.432**DVI120.365**DVI230.308**ARVI0.515**SR14-0.419**SR32-0.361**SR350.302** 注:**表示相關(guān)性檢驗(yàn)顯著。

由表2可知，有20個(gè)自變量因子在0.01水平與碳儲(chǔ)量顯著相關(guān)，其中SR34與碳儲(chǔ)量的相關(guān)性最高，相關(guān)系數(shù)達(dá)到-0.632。光譜變量與碳儲(chǔ)量的相關(guān)系數(shù)比碳儲(chǔ)量與紋理因子和地理因子的相關(guān)系數(shù)高。

3.2 多元逐步回歸模型

一般來說，擬合模型自變量越多，估測值的誤差越小，模型性能越好。然而自變量增多會(huì)加大工作量，其中一些不顯著因子會(huì)影響模型效果。因此，選擇合適的自變量數(shù)目非常重要。研究用SPSS 22.0軟件對(duì)樣本數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化。由表3可知，經(jīng)過4次擬合，4個(gè)自變量均進(jìn)入了回歸方程。隨著自變量的增加，方程的決定系數(shù)R2隨之增大，估計(jì)標(biāo)準(zhǔn)誤差隨之減小。當(dāng)變量個(gè)數(shù)為4個(gè)時(shí)，R2達(dá)到0.608，校正R2達(dá)到0.607，模型效果較好。

因此，森林碳密度的最優(yōu)回歸方程包括SR34、SR46、NDVI、EVI共4個(gè)自變量，具體見式(5)。

y=-3.488+0.412x1+1.577x2+0.009x3-

1.707x4

(5)

式中：x1為SR34；x2為SR46；x3為NDVI；x4為EVI；y為碳儲(chǔ)量。

表3 基于Landsat 8的多元逐步回歸分析模型統(tǒng)計(jì)Tab.3 Statistics of LSR model based on Landsat 8 image因變量模型自變量RR2校正R2估計(jì)標(biāo)準(zhǔn)誤差F顯著水平1SR340.4970.2470.2397.02821.5320.000 2SR34, SR460.6730.4530.4416.02418.5230.000 森林碳儲(chǔ)量3SR34, SR46,NDVI0.7630.5820.5719.41612.1650.0004SR34, SR46,NDVI,EVI0.780 0.6080.6075.1479.1370.001

3.3 偏最小二乘回歸模型

研究篩選20個(gè)顯著相關(guān)的因子參與回歸分析，設(shè)置潛在因子最大數(shù)為20，逐次分析隨潛在因子量的增大自變量和因變量的累積解釋量發(fā)生的變化，選取最優(yōu)偏最小二乘模型。累積解釋量隨潛在因子數(shù)量變化過程如圖1所示。

由圖1可知，使用20個(gè)自變量參與建模，當(dāng)潛在因子數(shù)量不同時(shí)，累積解釋量也會(huì)發(fā)生變化。大于8個(gè)時(shí)，自變量累積解釋量增長幅度趨于平緩，因變量累積解釋量則出現(xiàn)下降趨勢。說明PLSR分析的最佳潛在因子數(shù)量為8個(gè)。計(jì)算20個(gè)自變量因子的投影重要性指標(biāo)，結(jié)果見圖2。

由圖2可知，有6個(gè)自變量因子的VIP值大于1，說明對(duì)模型的貢獻(xiàn)性強(qiáng)；2個(gè)因子VIP值小于1而被剔除。排在前8位的自變量因子依次為SR34、DVI12、NDVI、SR32、SAVI025、ARVI、SAVI01、SR46，偏最小二乘回歸模型見式(6)。

y=5.162+0.097x1+0.677x2+0.007x3+

0.635x4-0.501x5-1.213x6-0.412x7-

0.144x8

(6)

式中：x1為SR34；x2為SR46；x3為NDVI；x4為SR32;x5為SAVI025;x6為 ARVI;x7為SAVI01;x8為DVI12；y為碳儲(chǔ)量。

3.4 RBF神經(jīng)網(wǎng)絡(luò)

將2/3的樣本(32個(gè))作為訓(xùn)練數(shù)據(jù)，篩選出顯著相關(guān)的變量20個(gè)，組成32×20的矩陣作為輸入向量，碳儲(chǔ)量作為輸出向量。調(diào)用Matlab軟件中的newrb函數(shù)創(chuàng)建神經(jīng)網(wǎng)絡(luò)，均方根誤差為0.001，最大神經(jīng)元個(gè)數(shù)為200，創(chuàng)建函數(shù)分布密度為0.1、0.2、0.3、0.4、0.5的5個(gè)神經(jīng)網(wǎng)絡(luò)模型，當(dāng)分布密度為0.1時(shí)，擬合效果最好。此時(shí)模型擬合的決定系數(shù)達(dá)到0.633，相對(duì)均方根誤差為15.250 t·hm-2，效果較好。

3.5 模型比較

以剩下的1/3樣本(16個(gè))對(duì)3種方法得到的森林碳儲(chǔ)量值進(jìn)行精度驗(yàn)證，分別計(jì)算R2和RRMSE，結(jié)果如表4所示。

由表4可知，RBF神經(jīng)網(wǎng)絡(luò)模型預(yù)測精度最高，R2達(dá)到0.645，RRMSE為15.582 t·hm-2；其次偏最小二乘回歸模型(PLSR)；多元逐步回歸模型(LSR)精度最低。在本研究區(qū)，神經(jīng)網(wǎng)絡(luò)模型在估測碳儲(chǔ)量時(shí)表現(xiàn)較好。

表4 3種建模方法的結(jié)果對(duì)比Tab.4 Comparisons of three models遙感影像模型R2RRMSE/(t·hm-2)LSR0.43118.105landsat 8PLSR0.51117.135RBF0.64515.582

3.6 桃江縣森林碳儲(chǔ)量空間分布

根據(jù)RBF神經(jīng)網(wǎng)絡(luò)模型繪制了桃江縣碳儲(chǔ)量空間分布圖(圖3)。由圖3可知，桃江縣海拔較高的地方碳儲(chǔ)量較大，分布在30～60 t·hm-2范圍內(nèi)；城區(qū)碳儲(chǔ)量較低，分布在0～30 t·hm-2范圍內(nèi)。

4 結(jié)論與討論

結(jié)合Landsat 8影像和2013年桃江縣一類調(diào)查數(shù)據(jù)，建立多元逐步回歸、偏最小二乘回歸和RBF神經(jīng)網(wǎng)絡(luò)模型，反演桃江縣碳儲(chǔ)量。

(1)利用遙感影像原始波段計(jì)算的植被指數(shù)與森林碳儲(chǔ)量相關(guān)性較高，可以準(zhǔn)確高效地預(yù)測森林碳儲(chǔ)量。

(2)RBF神經(jīng)網(wǎng)絡(luò)模型效果最好，決定系數(shù)和相對(duì)均方根誤差分別為0.645和15.582t·hm-2；其次是偏最小二乘回歸模型，分別是0.511和17.135 t·hm-2。多元逐步回歸模型估算精度最低，分別為0.431和18.105 t·hm-2。

(3)高海拔地區(qū)碳儲(chǔ)量大，城市地區(qū)碳儲(chǔ)量較小，符合桃江縣植被覆蓋實(shí)際情況。

與現(xiàn)有研究[16-18]相比，本研究引入RBF神經(jīng)網(wǎng)絡(luò)優(yōu)于傳統(tǒng)經(jīng)驗(yàn)?zāi)Ｐ?，擬合精度高，估計(jì)誤差小，無局部極小問題存在，且學(xué)習(xí)過程收斂速度快。但仍存在一些不足，比如參數(shù)選定較為困難，不能很好的解釋推理過程和推理依據(jù)。下一步研究重點(diǎn)將集中在優(yōu)化模型參數(shù)選取。