基于Landsat-8 OLI數據的馬尾松林蓄積量飽和點確定及估測

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

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

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

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

1 研究區(qū)概況

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

2 研究方法

2.1 數據源

2.1.1樣地調查數據

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

V=0.95682492×D1.8551497×H0.000062341803

(1)

式中:V為材積；D為胸徑；H為樹高。

表1 樣地數據統(tǒng)計信息

2.1.2Landsat-8 OLI數據

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

依據有關研究[17],植被指數選擇差異植被指數(Difference Vegetation Index,DVI)、調整差異植被指數(Renormalized Difference Vegetation Index,RDVI)、增強植被指數(Enhanced Vegetation Index,

表2 Landsat-8 OLI波段信息

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

很多時候，換一種說話方式，彼此的心情會截然不同。夫妻相伴多年，更要好好說話，多表達關心，少一些指責，只有這樣，那個伴兒才會長久地陪你走下去。

DVI=NIR-Red

(2)

(3)

(4)

(5)

(6)

LAI=3.618EVI-0.118

(7)

(8)

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

2.2 飽和值確定

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

(9)

式中:c0為塊金常數；c為拱高；(c0+c)為基臺值；a為變程。當c0=0,c=1時,稱為標準的球狀模型。

y=a0+a1x+a2x3

(10)

通過對式(10)求解即可得到各參數值。

該簡化模型是一種曲線擬合模型,本文嘗試提出使用二項式模型進行飽和點求解,并與前人提出的方法進行對比,二項式模型的具體形式為:

y=a0+a1x+a2x2

(11)

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

2.3 蓄積量估測

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

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

(12)

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

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

飽和點估測模型是本研究基于飽和點估測方法所提出的蓄積量估測方法,即基于二項式模型對蓄積量進行估測。首先,對建模變量進行相關性分析,得出與蓄積量相關性最大的因子,基于該最相關因子對蓄積量進行建模預測,其模型可表示為:

y=b0+b1x+b2x2

(13)

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

2.4 模型評價指標

基于模型評價指標[22-24],本文選擇決定系數R2,平均絕對誤差MAE和均方根誤差RMSE作為模型評價指標,具體計算公式為:

(14)

(15)

(16)

3 結果與分析

3.1 波段飽和點

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

表3 波段變量相關性分析

表4 模型估測參數

圖1 基于B3的蓄積量飽和點確定

3.2 植被指數飽和點

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

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

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

表5 植被指數變量相關性分析

表6 模型估測參數

圖2 基于GNDVI的蓄積量飽和點確定

表7 多元線性逐步回歸建模結果

著性水平設置為0.05,使用4個自由度的全球檢驗評估線性模型假設,假設檢驗全部通過。對所建模型進行精度驗證,驗證結果如圖3所示,模型預測蓄積量值與樣地測量值的R2為0.29,MAE為63.35 m3/hm2,RMSE為69.53 m3/hm2。

圖3 多元線性逐步回歸模型驗證結果

3.4 基于二項式(飽和點估測)模型的蓄積量估測

對隨機分好組的建模數據(n=51)使用二項式模型進行蓄積量估測建模,建模變量僅選擇B3,建模結果如表8所示。對所建模型進行精度驗證,驗證結果如圖4所示,該模型的預測蓄積量值與樣地測量值的R2為0.49,MAE為53.76 m3/hm2,RMSE為58.71 m3/hm2。

表8 二項式模型蓄積量估測建模結果

圖4 二項式模型驗證結果

4 結論與討論

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

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

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