孫忠秋,吳發(fā)云,胡 楊,高顯連,高金萍

(1.國(guó)家林業(yè)和草原局調(diào)查規(guī)劃設(shè)計(jì)院,北京 100714；2.寧夏大學(xué),銀川 750021)

基于遙感數(shù)據(jù)進(jìn)行森林參數(shù)提取是林業(yè)領(lǐng)域中進(jìn)行森林資源監(jiān)測(cè)的一項(xiàng)重要手段,通過(guò)遙感數(shù)據(jù)的光譜信息,聯(lián)合地面樣地調(diào)查點(diǎn)的數(shù)據(jù)進(jìn)行建模與反演,便可獲得一定精度范圍內(nèi)的森林資源估測(cè)結(jié)果。目前基于遙感數(shù)據(jù)進(jìn)行森林資源估測(cè)的研究有很多,如毛學(xué)剛等[1]使用Landsat長(zhǎng)時(shí)間序列數(shù)據(jù)估測(cè)了我國(guó)東北部分區(qū)域的樹(shù)高,結(jié)果表明使用長(zhǎng)時(shí)間序列的數(shù)據(jù)比單時(shí)相數(shù)據(jù)估測(cè)效果好,RMSE降低了0.50 m而模型估測(cè)精度提高了6.34%；周小成等[2]基于兩期高分辨率無(wú)人機(jī)影像估算了福建省將樂(lè)縣的杉木馬尾松林蓄積量,通過(guò)布料模擬濾波算法生成了研究區(qū)的數(shù)字表面模型、數(shù)字高程模型和冠層高度模型,對(duì)林分平均高和蓄積量的估測(cè)精度分別為94.74%和91.46%；韓宗濤等[3]基于Landsat-8波段信息、植被指數(shù)以及地形因子、紋理因子和SAR數(shù)據(jù)估測(cè)了大興安嶺根河地區(qū)的森林生物量,通過(guò)對(duì)比KNN-FIFS法和多元逐步回歸方法,發(fā)現(xiàn)KNN-FIFS方法估測(cè)效果更好(R2=0.77,RMSE=22.74 t/hm2)；龐勇等[4]基于機(jī)載高分辨率航空影像估測(cè)了小興安嶺地區(qū)的森林生物量,其基于FOTO算法提取的紋理因子聯(lián)合多元線性逐步回歸方法在估測(cè)森林生物量時(shí)發(fā)現(xiàn)無(wú)飽和現(xiàn)象,模型R2達(dá)到0.81,RMSE為46.78 t/hm2；嚴(yán)恩萍等[5]基于Landsat5和MODIS數(shù)據(jù)估測(cè)了湖南省攸縣森林的碳密度分布,發(fā)現(xiàn)Landsat5在碳密度估測(cè)時(shí)表現(xiàn)更優(yōu),估測(cè)精度達(dá)到82.02%；祖笑鋒等[6]基于MODIS產(chǎn)品建立了大區(qū)域燃燒生物量模型,其研究方法很好地反映了我國(guó)受火災(zāi)影響減少的森林生物量；郭云等[7]基于Landsat5數(shù)據(jù)估測(cè)了黑河流域地區(qū)的森林生物量,使用ASTER GDEM產(chǎn)品為輔助數(shù)據(jù),地形校正能提高生物量的估測(cè)精度,在多元線性逐步回歸算法下,地形校正前后的R2達(dá)到分別為0.31和0.46,RMSE分別為34.41,30.51 t/hm2；戚玉嬌等[8]基于Landsat5數(shù)據(jù)估測(cè)了黑龍江大興安嶺地區(qū)的森林生物量,發(fā)現(xiàn)存在高值低估和低值高估的現(xiàn)象；李亦秋等[9]基于Landsat TM影像和多元線性逐步回歸方法估測(cè)山東省的蓄積量,其建立的模型估測(cè)精度達(dá)到87.35%。在國(guó)際上,許多學(xué)者同樣基于遙感數(shù)據(jù)對(duì)森林參數(shù)開(kāi)展了相關(guān)估測(cè)研究[10-15]。

目前多數(shù)研究基于多元線性逐步回歸模型方法使用遙感數(shù)據(jù)估測(cè)森林參數(shù),而針對(duì)蓄積量飽和點(diǎn)的估測(cè)的研究則很少。本文選擇湖南省西北部馬尾松林為研究對(duì)象,對(duì)Landsat-8 OLI數(shù)據(jù)估測(cè)馬尾松林蓄積量的飽和點(diǎn)進(jìn)行定量化估測(cè),從根本上了解Landsat-8 OLI數(shù)據(jù)估測(cè)馬尾松林蓄積量的能力。此外,基于飽和點(diǎn)估測(cè)方法,發(fā)展二項(xiàng)式估測(cè)模型對(duì)馬尾松林蓄積量進(jìn)行建模估測(cè)并與多元逐步回歸模型進(jìn)行對(duì)比,以期評(píng)價(jià)二項(xiàng)式模型和多元線性回歸模型估測(cè)森林蓄積量的能力,為森林資源的快速獲取提供新的估測(cè)方法。

1 研究區(qū)概況

研究區(qū)位于湖南省,該區(qū)森林資源豐富,尤其是馬尾松林 (Pinusmassoniana),在該地區(qū)的森林資源中占有非常大的比重[16]。在第八次全國(guó)森林資源清查中,湖南省馬尾松林的蓄積量占全省森林總蓄積量的14.02%。

2 研究方法

2.1 數(shù)據(jù)源

2.1.1樣地調(diào)查數(shù)據(jù)

地面樣地調(diào)查在2017年進(jìn)行,為了保證調(diào)查的樣地具有代表性,首先對(duì)樣地的郁閉度進(jìn)行等級(jí)劃分:0.20～0.39為第一郁閉度等級(jí)；0.40～0.69為第二郁閉度等級(jí)；0.70～1.00為第三郁閉度等級(jí)。各郁閉度等級(jí)分別調(diào)查了30,24,18塊樣地,最終共調(diào)查馬尾松林地面樣地72塊。其中,幼齡林樣地2塊,占比2.78%；中齡林樣地26塊,占比36.11%；近熟林樣地25塊,占比34.72%；成熟林樣地18塊,占比25%；過(guò)熟林樣地1塊,占比1.39%。對(duì)所有的樣地內(nèi)胸徑大于5 cm的樣木進(jìn)行每木檢查,樹(shù)高測(cè)量使用VL5超聲波測(cè)高器,胸徑測(cè)量使用胸徑尺,獲得樣地內(nèi)每株樣木的樹(shù)高、胸徑等單木因子?；谛罘e量估測(cè)方程計(jì)算每株樣木的材積,然后求和得到樣地尺度的蓄積量,通過(guò)尺度擴(kuò)展,得到調(diào)查樣地的每公頃蓄積量值。使用GNSS RTK技術(shù)對(duì)樣地和樣木進(jìn)行精確定位,以保證調(diào)查數(shù)據(jù)可以和遙感像元進(jìn)行精確的匹配。最后,基于樣地中心和樣木坐標(biāo)把圓形樣地重采至25 m×25 m的樣地?cái)?shù)據(jù)。使用所有的數(shù)據(jù)進(jìn)行飽和點(diǎn)計(jì)算,在蓄積量估測(cè)階段,把樣地?cái)?shù)據(jù)進(jìn)行隨機(jī)分組,70%的樣地?cái)?shù)據(jù)用于建模,30%的樣地?cái)?shù)據(jù)用于模型驗(yàn)證,樣地?cái)?shù)據(jù)具體描述如表1所示。樣木材積計(jì)算公式如下:

V=0.95682492×D1.8551497×H0.000062341803

(1)

式中:V為材積；D為胸徑；H為樹(shù)高。

表1 樣地?cái)?shù)據(jù)統(tǒng)計(jì)信息

2.1.2Landsat-8 OLI數(shù)據(jù)

Landsat-8 OLI數(shù)據(jù)基于GEE平臺(tái)獲取,通過(guò)GEE平臺(tái)直接對(duì)處理好的地表反射率產(chǎn)品進(jìn)行提取,該產(chǎn)品為L(zhǎng)ANDSAT/LC08/C01/T1_SR,為了匹配地面樣地調(diào)查數(shù)據(jù)日期,同時(shí)兼顧對(duì)研究區(qū)進(jìn)行覆蓋,產(chǎn)品日期選擇20170501—20171031和20180501—20181031?；谠摦a(chǎn)品的pixel_qa質(zhì)量屬性波段,對(duì)該產(chǎn)品進(jìn)行去云處理,去云處理使用GEE官方公布的方法,直接掩膜掉pixel_qa波段像元屬性為3和5的像元。去云后,使用中值函數(shù)對(duì)覆蓋的影像進(jìn)行中值處理,之后進(jìn)行重采樣(25 m×25 m)和坐標(biāo)系統(tǒng)匹配。最后,使用地面樣地矢量文件對(duì)處理好的數(shù)據(jù)進(jìn)行波段信息和植被指數(shù)提取。其中,波段信息如表2所示。

依據(jù)有關(guān)研究[17],植被指數(shù)選擇差異植被指數(shù)(Difference Vegetation Index,DVI)、調(diào)整差異植被指數(shù)(Renormalized Difference Vegetation Index,RDVI)、增強(qiáng)植被指數(shù)(Enhanced Vegetation Index,

表2 Landsat-8 OLI波段信息

EVI)、歸一化植被指數(shù)(Normalized Difference Vegetation Index,NDVI)、綠度歸一化植被指數(shù)(Green Normalized Difference Vegetation Index,GNDVI)、葉面積指數(shù)(Leaf Area Index,LAI)、簡(jiǎn)單比值植被指數(shù)(Simple Ratio Vegetation Index,RVI),各植被指數(shù)的具體計(jì)算公式如下:

很多時(shí)候，換一種說(shuō)話方式，彼此的心情會(huì)截然不同。夫妻相伴多年，更要好好說(shuō)話，多表達(dá)關(guān)心，少一些指責(zé)，只有這樣，那個(gè)伴兒才會(huì)長(zhǎng)久地陪你走下去。

DVI=NIR-Red

(2)

(3)

(4)

(5)

(6)

LAI=3.618EVI-0.118

(7)

(8)

式中:NIR為近紅外區(qū)反射率；Red為紅光區(qū)反射率；Blue為藍(lán)光區(qū)反射率；Green為綠光區(qū)反射率。

2.2 飽和值確定

光譜信息只能反映一定數(shù)據(jù)變化范圍內(nèi)的森林參數(shù)(葉面積指數(shù)、蓄積量、生物量、碳儲(chǔ)量等),當(dāng)某一森林參數(shù)超出某一范圍時(shí),光譜信息并不能對(duì)其進(jìn)行很好的表達(dá),通常把這個(gè)臨界點(diǎn)定義為飽和點(diǎn),飽和點(diǎn)是定量化確定遙感數(shù)據(jù)估測(cè)森林參數(shù)能力的重要指標(biāo)。有學(xué)者提出基于地統(tǒng)計(jì)學(xué)模型中半方差函數(shù)理論來(lái)求解森林生物量飽和點(diǎn),并用此方法計(jì)算了Landsat-8數(shù)據(jù)估測(cè)浙江省森林生物量的飽和點(diǎn)[18]。該模型具體公式為:

(9)

式中:c0為塊金常數(shù)；c為拱高；(c0+c)為基臺(tái)值；a為變程。當(dāng)c0=0,c=1時(shí),稱為標(biāo)準(zhǔn)的球狀模型。

y=a0+a1x+a2x3

(10)

通過(guò)對(duì)式(10)求解即可得到各參數(shù)值。

該簡(jiǎn)化模型是一種曲線擬合模型,本文嘗試提出使用二項(xiàng)式模型進(jìn)行飽和點(diǎn)求解,并與前人提出的方法進(jìn)行對(duì)比,二項(xiàng)式模型的具體形式為:

y=a0+a1x+a2x2

(11)

式中:y表示波段光譜反射率；x表示蓄積量；a0,a1,a2為模型參數(shù)。

2.3 蓄積量估測(cè)

確定Landsat-8 OLI數(shù)據(jù)估測(cè)湖南省馬尾松林的光譜飽和點(diǎn)后,使用該光譜數(shù)據(jù)對(duì)馬尾松林進(jìn)行建模預(yù)測(cè),建模方法選擇多元逐步回歸方法和基于飽和點(diǎn)確定模型的估測(cè)方法,對(duì)比2種方法估測(cè)馬尾松林蓄積量的能力。

多元逐步回歸是多元線性回歸的一種,其通過(guò)逐個(gè)對(duì)建模變量進(jìn)行引入,不斷剔除對(duì)建模結(jié)果無(wú)意義的變量,最終達(dá)到優(yōu)選變量進(jìn)行建模預(yù)測(cè)的目的[19-21]。其原理是:每一步只引入或剔除一個(gè)自變量,自變量是否引入或剔除取決于其偏回歸平方和的F檢驗(yàn)或校正決定系數(shù)。如方程中已引入了(m-1)個(gè)自變量,在此基礎(chǔ)上考慮再引入變量Xi。記引入Xi后方程(即含m個(gè)自變量)的回歸平方和為SS回歸,殘差為SS殘差；之前含(m-1)個(gè)自變量(不包含Xi)方程的回歸平方和為SS回歸(-i),則Xi的偏回歸平方和為U=SS回歸-SS回歸(-i),檢驗(yàn)統(tǒng)計(jì)量為(12):

(12)

式中：Fi服從Fα(1,n-m-1)分布。

如果Fi>Fα(1,n-m-1),則Xi選入模型；否則,不入選。從方程中剔除無(wú)統(tǒng)計(jì)學(xué)作用的自變量,過(guò)程則相反,但檢驗(yàn)一樣。

飽和點(diǎn)估測(cè)模型是本研究基于飽和點(diǎn)估測(cè)方法所提出的蓄積量估測(cè)方法,即基于二項(xiàng)式模型對(duì)蓄積量進(jìn)行估測(cè)。首先,對(duì)建模變量進(jìn)行相關(guān)性分析,得出與蓄積量相關(guān)性最大的因子,基于該最相關(guān)因子對(duì)蓄積量進(jìn)行建模預(yù)測(cè),其模型可表示為:

y=b0+b1x+b2x2

(13)

式中:y表示蓄積量；x表示波段光譜反射率；a0,a1,a2為模型參數(shù)。

2.4 模型評(píng)價(jià)指標(biāo)

基于模型評(píng)價(jià)指標(biāo)[22-24],本文選擇決定系數(shù)R2,平均絕對(duì)誤差MAE和均方根誤差RMSE作為模型評(píng)價(jià)指標(biāo),具體計(jì)算公式為:

(14)

(15)

(16)

3 結(jié)果與分析

3.1 波段飽和點(diǎn)

基于SPSS 25.0對(duì)Landsat-8提取的波段變量進(jìn)行皮爾遜(Person)相關(guān)性分析,并進(jìn)行雙尾相關(guān)性顯著檢驗(yàn),結(jié)果如表3所示。所有變量均在0.01水平上顯著相關(guān),且B3波段與森林蓄積量(FSV)的相關(guān)性最好,為-0.625?；贐3變量求解馬尾松林蓄積量的估測(cè)飽和點(diǎn),基于球狀模型和二項(xiàng)式的模型估測(cè)參數(shù)如表4所示,得到基于球狀模型估測(cè)的最大蓄積量飽和點(diǎn)為217.05 m3/hm2(圖1紅線),基于二項(xiàng)式模型估測(cè)的最大蓄積量飽和點(diǎn)為206.86 m3/hm2(圖1綠線)。

表3 波段變量相關(guān)性分析

表4 模型估測(cè)參數(shù)

圖1 基于B3的蓄積量飽和點(diǎn)確定

3.2 植被指數(shù)飽和點(diǎn)

基于SPSS 25.0對(duì)Landsat-8提取的波段變量進(jìn)行皮爾遜(Person)相關(guān)性分析,并進(jìn)行雙尾相關(guān)性顯著檢驗(yàn),結(jié)果如表5所示。除了RDVI,EVI和LAI變量與FSV不顯著相關(guān)外,其余所有變量均在0.01和0.05水平上顯著相關(guān),且GNDVI波段與FSV的相關(guān)性最好,為0.451?；贕NDVI變量求解馬尾松林蓄積量的估測(cè)飽和點(diǎn),方案1和方案2的模型估測(cè)參數(shù)如表6所示,得到方案1估測(cè)最大蓄積量飽和點(diǎn)為196.95 m3/hm2(圖2紅線),方案2估測(cè)的最大蓄積量飽和點(diǎn)為183.06 m3/hm2(圖2綠線)。

3.3 基于多元線性逐步回歸的蓄積量估測(cè)

對(duì)隨機(jī)分好組的建模數(shù)據(jù)(n=51)進(jìn)行多元線性逐步回歸建模,波段變量和植被指數(shù)變量均逐步進(jìn)入回歸模型,建模結(jié)果如表7所示。結(jié)果表明,與波段信息相比,植被指數(shù)在估測(cè)湖南省馬尾松林蓄積量的表現(xiàn)較差,沒(méi)有變量進(jìn)入到預(yù)測(cè)模型中。顯

表5 植被指數(shù)變量相關(guān)性分析

表6 模型估測(cè)參數(shù)

圖2 基于GNDVI的蓄積量飽和點(diǎn)確定

表7 多元線性逐步回歸建模結(jié)果

著性水平設(shè)置為0.05,使用4個(gè)自由度的全球檢驗(yàn)評(píng)估線性模型假設(shè),假設(shè)檢驗(yàn)全部通過(guò)。對(duì)所建模型進(jìn)行精度驗(yàn)證,驗(yàn)證結(jié)果如圖3所示,模型預(yù)測(cè)蓄積量值與樣地測(cè)量值的R2為0.29,MAE為63.35 m3/hm2,RMSE為69.53 m3/hm2。

圖3 多元線性逐步回歸模型驗(yàn)證結(jié)果

3.4 基于二項(xiàng)式(飽和點(diǎn)估測(cè))模型的蓄積量估測(cè)

對(duì)隨機(jī)分好組的建模數(shù)據(jù)(n=51)使用二項(xiàng)式模型進(jìn)行蓄積量估測(cè)建模,建模變量?jī)H選擇B3,建模結(jié)果如表8所示。對(duì)所建模型進(jìn)行精度驗(yàn)證,驗(yàn)證結(jié)果如圖4所示,該模型的預(yù)測(cè)蓄積量值與樣地測(cè)量值的R2為0.49,MAE為53.76 m3/hm2,RMSE為58.71 m3/hm2。

表8 二項(xiàng)式模型蓄積量估測(cè)建模結(jié)果

圖4 二項(xiàng)式模型驗(yàn)證結(jié)果

4 結(jié)論與討論

1) 基于本文提出的方法,在對(duì)某一區(qū)域的森林蓄積量進(jìn)行估測(cè)時(shí),就可以根據(jù)其蓄積量的大致分布范圍選擇適合對(duì)其進(jìn)行估測(cè)的遙感數(shù)據(jù)源。

2) 基于前人研究中提出的基于球狀模型的飽和點(diǎn)估測(cè)方法,本研究提出的飽和點(diǎn)估測(cè)方法更好,原因主要有以下兩點(diǎn):一是球狀模型是曲線模型,而二項(xiàng)式模型也是曲線模型,在一定的數(shù)據(jù)范圍內(nèi)兩個(gè)數(shù)據(jù)擬合曲線基本重合,具有較高的相似性,從而在一定程度上具有可比性；二是相較于球狀模型,二項(xiàng)式模型對(duì)數(shù)據(jù)的擬合結(jié)果更好,在波段信息飽和點(diǎn)估測(cè)階段,球狀模型和二項(xiàng)式模型的擬合R2分別為0.53和0.56；在植被指數(shù)估測(cè)階段,球狀模型和二項(xiàng)式模型的擬合R2分別為0.39和0.43。在蓄積量估測(cè)建模階段,盡管有許多研究基于多元線性回歸模型進(jìn)行森林參數(shù)建模預(yù)測(cè),本文研究表明基于多元線性回歸模型對(duì)遙感蓄積量估測(cè)研究并不合適(R2=0.29,MAE=63.35 m3/hm2,RMSE=69.53 m3/hm2),即遙感變量與蓄積量并不是一種簡(jiǎn)單的線性關(guān)系,基于二項(xiàng)式模型的單變量估測(cè)方法的估測(cè)精度更高(R2=0.49,MAE=53.76 m3/hm2,RMSE=58.71 m3/hm2)。

3) 由于湖南省全年受云量的影響,對(duì)Landsat-8遙感數(shù)據(jù)的選擇提出了很高的要求,本文使用的Landsat-8數(shù)據(jù)是基于2017年和2018年2年的數(shù)據(jù)合成產(chǎn)品,這可能與同一數(shù)據(jù)源單時(shí)相的遙感產(chǎn)品的估測(cè)結(jié)果存在一定的差異。對(duì)特定地區(qū)的特定森林類型應(yīng)進(jìn)行特定的分析,研究區(qū)不同、森林類型的差異對(duì)估測(cè)結(jié)果都會(huì)產(chǎn)生一定的影響。