應用Sentinel-2A衛(wèi)星遙感影像估測森林蓄積量

龐曉燕 劉海松 年學東 張志超 李崇貴 董相軍

(內(nèi)蒙古自治區(qū)林業(yè)監(jiān)測規(guī)劃院，呼和浩特，010031) (西安科技大學) (內(nèi)蒙古自治區(qū)克什克騰旗白音敖包機械防火站)

森林蓄積量是衡量一個國家或地區(qū)森林資源的重要指標，也是制定林業(yè)經(jīng)營方案、實現(xiàn)森林資源可持續(xù)發(fā)展的重要依據(jù)。及時掌握森林資源現(xiàn)狀及其發(fā)展變化趨勢是林業(yè)科技工作者關注的熱點。傳統(tǒng)森林蓄積量調(diào)查方法主要有標準木法、采集表法等，需要進行大量的野外調(diào)查，費用昂貴、耗時費力[1]。隨著3S技術(shù)的快速發(fā)展及應用，國內(nèi)外林業(yè)科技工作者對森林蓄積量估測進行了大量的研究試驗，但多應用在縣/林業(yè)局或林場級尺度上，而在小班尺度的應用研究較少[2-5]。國內(nèi)外學者多用Landsat TM系列、Modis、SPOT5、資源三號、高分系列等遙感數(shù)據(jù)進行森林蓄積量估測，利用Sentinel-2A遙感影像進行蓄積量估測只有少量的報道。

Sentinel-2A衛(wèi)星是全球環(huán)境與安全監(jiān)測(GMES)計劃的第2顆衛(wèi)星，采用的多光譜成像儀。衛(wèi)星遙感影像依據(jù)光譜帶的不同，空間分辨率分為10、20和60 m，從可見光、近紅外到短波紅外覆蓋了13個譜段，其中包含的3個紅邊波段是植被長勢監(jiān)測的理想數(shù)據(jù)源。3個紅邊波段在全球環(huán)境變化、農(nóng)作物生長監(jiān)測、森林監(jiān)測、災害應急與救援等方面具有重要價值[6]。何云等[7]利用Sentinel-2A影像結(jié)合隨機森林方法對中南半島穆河流域的東南部區(qū)域?qū)ν恋馗采w進行了分類，分類精度達到87.53%，Kappa系數(shù)達到0.8461；陳瑜云[8]應用三個時期Sentinel-2影像數(shù)據(jù)對毛竹林區(qū)分大小年進行生物量了估測，實驗表明7月數(shù)據(jù)的RF模型精度最高，總體樣本驗證R2為0.46，均方根誤差為10.68 Mg/ha，相對均方根誤差為19.09%；韓濤等[9]針對Sentinel-2A與Landsat-8數(shù)據(jù)在油菜識別方面的差異性進行了研究，結(jié)果表明Sentinel-2A數(shù)據(jù)能夠在種植結(jié)構(gòu)復雜的小尺度區(qū)域提取更高精度的作物分布信息；蘇偉等[10]利用Sentinel-2A影像對玉米冠層葉綠素進行估算，表明引入紅邊波段后玉米冠層葉綠素含量的估算精度的決定系數(shù)為0.796；高璐璐[11]基于Sentinel-2A數(shù)據(jù)反演低山丘陵區(qū)蘋果樹冠層葉綠素含量，研究表明以植被指數(shù)建立的支持向量機反演模型較好，決定系數(shù)為0.729。綜上所述，Sentinel-2A數(shù)據(jù)在不同的領域的應用均取得了良好效果。

本研究以內(nèi)蒙古自治區(qū)某林業(yè)局的數(shù)據(jù)(包括：一類清查樣地數(shù)據(jù)、二類調(diào)查小班數(shù)據(jù)、數(shù)字高程模型(DEM)以及林地數(shù)據(jù))和Sentinel-2A影像為數(shù)據(jù)源，采用k-近鄰(k-NN)法、穩(wěn)健估計法及偏最小二乘法，以樣地大小為單位，分別按照針葉林、闊葉林和針闊混交林等3種林型構(gòu)建的森林蓄積量估測模型，并從林業(yè)局、林場和小班等3個尺度討論應用Sentinel-2A影像進行森林蓄積量估測的可行性。

1 研究區(qū)概況

研究區(qū)位于內(nèi)蒙古自治區(qū)東北部，地理坐標為東經(jīng)120°36′～121°46′，北緯47°5′～47°45′，屬于大興安嶺興安落葉松林區(qū)的南延部分，森林覆蓋率為75%。該區(qū)域主要以興安落葉松(Larixgmelinii(Rupr.) Kuzen)為主，同時生長有白樺(BetulaplatyphyllaSuk.)、山楊(Populusdavidiana)、黑樺(BetuladahuricaPall.)、蒙古櫟(QuercusmongolicaFisch. ex Ledeb)等樹種。森林類型主要以蒙古櫟、黑樺、白樺和興安落葉松組成的針闊混交林。

2 研究方法

2.1 數(shù)據(jù)獲取與處理

數(shù)據(jù)來源于研究區(qū)2018年森林資源連續(xù)清查固定樣地數(shù)據(jù)和依據(jù)2016年二類調(diào)查小班數(shù)據(jù)系統(tǒng)布設的解譯樣地。2018年一類清查樣地數(shù)據(jù)為0.06 hm2的長方形(10 m×60 m)樣地，樣地長邊為南北方向，樣地坐標采集點位于樣地中心，樣地總數(shù)為55個，有林地43個，無法滿足依據(jù)林地類型進行建模估測的需求，故參照一類樣地的大小和形狀在試驗區(qū)內(nèi)按照4 km×4 km的間距布設解譯樣地，將解譯樣地所落小班的公頃蓄積乘以0.06推算各個樣地的蓄積。綜合兩類樣地數(shù)據(jù)共得到175個有林地樣地，其中針葉林樣地63個，闊葉林樣地80個，針闊混交林樣地32個。

遙感影像采用Sentinel-2A衛(wèi)星遙感影像，成像時間為2016年8月28日。Sentinel-2A的Level-1C數(shù)據(jù)經(jīng)過了輻射校正、幾何校正、正射校正和亞像素精度全球參考系統(tǒng)的空間配準[12]。先利用Sen2Cor對影像數(shù)據(jù)進行大氣校正，在利用SNAP軟件進行重采樣處理，將所有波段全部重采樣至分辨率為10 m，最后利用ENVI5.3SP1的Layer Stacking將影像數(shù)據(jù)波段2～波段8A進行組合，獲得具有8個波段的遙感影像[11-12]。將樣地數(shù)據(jù)疊加到遙感影像上，即可提取樣地范圍對應的遙感信息(見圖1)。

數(shù)字高程模型采用在地理空間數(shù)據(jù)云下載的30 m分辨率的DEM數(shù)據(jù)，用于獲取研究區(qū)域的坡度和坡向信息，并作為可能影響蓄積量估測的地形因子參與建模。利用ArcGIS10.1將研究區(qū)2017年林地數(shù)據(jù)生成林型分布圖，包括針葉林、闊葉林和針闊混交林等3種類型(見圖2)。

2.2 蓄積量估測模型構(gòu)建方法

設在監(jiān)測區(qū)域共計布設了n個樣地，每個樣地設置了可能影響蓄積量估測的p個自變量X1,X2，…，Xp，每個樣地就是p維歐式空間中的一點。將n個樣地展繪在校正好的遙感影像上，提取樣地對應的遙感影像各波段的灰度信息，并構(gòu)造比值波段。考慮地形信息，得到可能影響蓄積量估測p個自變量因子的值。經(jīng)過最優(yōu)變量選擇和建模樣地抽樣后，依據(jù)所抽樣地分林地類型建立以樣地大小為單位的蓄積量估測模型。將遙感影像轉(zhuǎn)換為以樣地大小為單位的影像，結(jié)合DEM數(shù)據(jù)獲取各樣地大小單元對應的自變量值，代入估測模型計算每個樣地大小單元對應的蓄積量[4-5，13]。在針對不同林地類型構(gòu)建蓄積量估測模型時，本文選用穩(wěn)健估計法、偏最小二乘估計法和k-近鄰法等3種方法建立模型。

穩(wěn)健估計法：穩(wěn)健估計對異常樣地具有較好的抗干擾能力，能夠有效消除異常值的影響，確保模型具有良好的穩(wěn)定性。

按照一定的比例，隨機抽取m個樣地參與建立估測模型，用其余n-m個樣地來檢驗模型的預報精度，假設單個樣地蓄積量v與其對應的q個自變量之間存在如下線性關系：

v=β0+β1X1+β2X2+…+βqXq。

故m個樣地蓄積v與q個自變量之間的關系可以表示為：

V=Xβ+eV,

E(eV)=0,

Cov(eV)=σ2I。

式中：V為因變量矩陣，V=(v1,v2，…，vm)′；β是待定參數(shù)向量，β=(β0，β1，…，βq)′；σ2為樣地蓄積量估測方差；eV是m個樣地蓄積量的觀測誤差向量，eV=(eV1,eV2,…,eVn)′；I為m階單位陣;X是m個樣地對應q個變量因子組成的自變量矩陣

則誤差方程表示為：

eV=Xβ-V。

穩(wěn)健估計用誤差的絕對最小值作為穩(wěn)健估計的極值函數(shù)，以便減弱異常樣地的影響，此時需滿足下式要求。

式中：ρ(ei)是極值函數(shù)。

建立穩(wěn)健估計蓄積量估測模型的重點是選用增長速度緩慢的有界函數(shù)作為極值函數(shù)，從而有效控制異常樣地對蓄積量估測模型中待定參數(shù)的影響。以ρ(e)=|e|為例，采用1-范數(shù)作為極值函數(shù)。本文采用如下迭代函數(shù)法進行求解。

Xiqβq-vi|。

利用迭代法求解穩(wěn)健估計模型中待定參數(shù)最優(yōu)解，其步驟如下：

③構(gòu)建法方程(X′W)(1)Xβ=X′W(1)V；

偏最小二乘法：偏最小二乘估計集多元逐步回歸分析，主成分分析和典型相關性分析的基本功能于一體，不僅能夠從自變量集中選取具有最佳解釋能力的新綜合變量，克服多重相關性造成的不利影響，還能夠有效區(qū)分系統(tǒng)中的信息和噪聲，提高型模的可靠性[15-16]。

偏最小二乘估計建模步驟如下：(1)設參與建模的m個被抽樣地蓄積量組成的向量為V，每個樣地p個自變量組成的觀測矩陣為Xm×p；(2)提取觀測矩陣X的一個主成分t1(t1是x1,x2,…,xp的線性組合)，要求t1不僅需要盡可能多地攜帶X中的變異信息，還要與V的相關性最大；(3)提取第1個主成分后，求出V和X對t1的回歸方程。若此時回歸方程達到精度要求，則算法停止，否則，需利用X和V被t1解釋后的殘余信息進行第2輪主成分提取，如此反復，直至滿足精度要求為止。若最終對X提取了q個主成分t1,t2,…,tq，偏最小二乘回歸將變?yōu)閂對q個主成分的回歸，據(jù)此再求出V對原變量X的回歸方程[14，17]。

k-近鄰法(k-NN)：k-NN法是一種典型的非參數(shù)化方法，核心思想是根據(jù)多維特征空間中待估測點與已知點之間的空間相似性進行單變量或多變量估測，相比于通過參數(shù)化方法的回歸分析建模，k-NN方法能夠更多的考慮蓄積量與自變量之間的非線性依賴關系，不受自變量間存在復共線性的影響[17-18]。估測公式如下：

在光譜空間中找出與待估測單元最鄰近的k個樣地，其中p1,p2,…,pk，其中Dp1,q

2.3 自變量設置與優(yōu)選

可能影響蓄積量估測的自變量主要包括影像波段灰度信息、比值波段和地形信息三種類型。影像波段灰度信息包括Sentinel-2A影像B2～B8、B8A共8個波段光譜的反射率信息。地形信息包括陰坡、陽坡和坡度。比值波段分為傳統(tǒng)光譜特征指數(shù)(歸一化植被指數(shù)、增強型植被指數(shù)、差值植被指數(shù)、改進紅邊歸一化植被指數(shù)、簡單比值植被指數(shù))和紅邊光譜特征指數(shù)(紅邊差值植被指數(shù)、紅邊歸一化植被指數(shù)、修正型葉綠素吸收反射率指數(shù)、特征色素簡單比值指數(shù)、差值植被指數(shù)45)[8-12，20-21]。

Sentinel-2A數(shù)據(jù)包含3個紅邊波段，參考現(xiàn)有資料，中心波段在705 nm處的反射率與葉綠素含量間的相關性優(yōu)于740 nm和783 nm處的反射率，因此本文紅邊光譜特征指數(shù)除PSSRa外，均選取波長為705 nm的波段作為紅邊波段[11-12]。比值波段構(gòu)造公式如表1所示。

表1 植被指數(shù)的計算公式

參與建模的變量并不是越多越好，有些變量可能對蓄積量估測影響較小甚至沒有影響，有些變量之間可能存在復共線性。過多的變量會造成數(shù)據(jù)冗余，在增大計算量的同時，也會影響模型預報精度。本文采用平均殘差平方和準則和窮舉算法對以上提及的21個變量進行建模變量優(yōu)選[15，17，22-23]。平均殘差平方和表達式為：

式中：q為建模變量個數(shù);V為因變量向量;X為觀測陣。

3 模型評價及精度檢驗

3.1 不同建模方法的綜合評價

按照建模樣本與驗證樣本7∶3的比例對樣地進行隨機抽樣，然后依據(jù)建模樣地進行最優(yōu)變量選擇，最后采用k-NN法、穩(wěn)健估計以及偏最小二乘估計等3種方法分林型進行建模。建模過程重復進行15次，以此來減少隨機抽樣給建模效果帶來的偶然性影響(見表2和圖3、4)。

表2 3種估測方法建模綜合評價

由表2、圖3～4可知，k-NN法建模穩(wěn)定性較好，隨機抽樣不會對建模效果帶來較大影響，而穩(wěn)健估計法和偏最小二乘法建模效果波動性較大，隨機抽樣帶來的偶然性影響較大。

圖3 各模型殘差相對均方根誤差

3.2 不同建模方法的最優(yōu)模型

按照模型殘差相對均方根誤差與預報偏差相對均方根誤差相近且兩者都在一倍標準差范圍內(nèi)的原則，對15次抽樣建模結(jié)果進行篩選，分別挑選出3種方法各自最具代表性的蓄積量估測模型。

圖4 各模型預報偏差相對均方根誤差

k-NN法選取第12次所建模型，最優(yōu)變量為：綠光波段(B3)、近紅外波段(B8)、窄邊近紅外波段(B8A)、紅邊差值植被指數(shù)和增強型植被指數(shù)；

穩(wěn)健估計法選取第13次所建模型，最優(yōu)變量為：植被紅邊1波段(B5)、近紅外波段(B8)、歸一化植被指數(shù)、簡單比值植被指數(shù)、坡度和陰坡；

偏最小二乘選取第3次所建模型，最優(yōu)變量為：綠光波段(B3)、紅光波段(B4)、植被紅邊1波段(B5)、近紅外波段(B8)、改進紅邊歸一化植被指數(shù)、歸一化差分指數(shù)、紅邊歸一化植被指數(shù)、簡單比值植被指數(shù)和陰坡，各模型精度如表3所示。

表3 各模型精度

如圖5～7所示，對k取不同值時的均方根誤差計算結(jié)果進行分析，針葉林、闊葉林、針闊混交林蓄積量估測模型的最佳k值分別為29、28和17。

圖5 針葉林蓄積量估測模型k取不同值時的均方根誤差

采用穩(wěn)健估計法所得估測方程如下：

0.087 289×P+0.062 703×ISR-0.0.377 512×Py+2.388 314；

式中：P表示坡度；Py表示陰坡。

圖6 闊葉林蓄積量估測模型k取不同值時的均方根誤差

圖7 針闊混交林蓄積量估測模型k取不同值時的均方根誤差

采用偏最小二乘估計所建估測方程如下：

式中：Py表示陰坡。

3.3 不同模型對蓄積量的估測結(jié)果與不同尺度的精度

利用Sentinel-2A遙感影像和DEM數(shù)據(jù)獲取各樣地單元對應自變量的取值，并代入模型估測各單元對應的蓄積量估測值，獲得研究區(qū)域森林蓄積量分布(見圖8)，再利用積分法統(tǒng)計林業(yè)局、林場和小班等3個尺度的蓄積量估測值。隨機抽取研究區(qū)域內(nèi)5%的有林地小班(699個)用于小班尺度的估測精度分析，并以二類調(diào)查小班數(shù)據(jù)為標準進行檢驗，結(jié)果如表4所示。

圖8 森林蓄積量估測等級分布圖

表4 不同尺度下估測精度比較

由表4可知，k-NN方法在3個尺度中的估測效果均明顯優(yōu)于其他兩種方法，對于3種尺度的森林蓄積量估測均能滿足精度要求，而穩(wěn)健估計法和偏最小二乘法在小班尺度的估測精度較差。

4 結(jié)論

本文運用Sentinel-2A衛(wèi)星影像提取樣地范圍的波段灰度信息并構(gòu)造比值波段，結(jié)合樣地數(shù)據(jù)和地形信息，建立以樣地為單位的森林蓄積量估測模型，對研究區(qū)的森林蓄積量進行反演，所得結(jié)果如下：(1)非參數(shù)化方法k-NN法在林業(yè)局、林場和小班等3個尺度中的估測精度分別達到了97.0%、93.2%和83.6%，均表現(xiàn)出良好的估測效果，明顯優(yōu)于其他兩種參數(shù)化方法。(2)穩(wěn)健估計和偏最小二乘估計的估測精度相差不大，均不能滿足小班尺度蓄積量調(diào)查的精度要求，如何提升這兩種方法在小班級尺度上的蓄積量估測精度有待進一步研究。(3)多次抽樣建模發(fā)現(xiàn)，k-NN法具有較好的建模穩(wěn)定性，而穩(wěn)健估計法和偏最小二乘法的建模效果波動性較大；若采用參數(shù)化方法進行建模，實際生產(chǎn)應進行多次抽樣以提高估測模型的穩(wěn)定性。