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

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

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

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

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

1 研究區(qū)概況

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

2 研究方法

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

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

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

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

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

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

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

按照一定的比例，隨機(jī)抽取m個(gè)樣地參與建立估測模型，用其余n-m個(gè)樣地來檢驗(yàn)?zāi)Ｐ偷念A(yù)報(bào)精度，假設(shè)單個(gè)樣地蓄積量v與其對應(yīng)的q個(gè)自變量之間存在如下線性關(guān)系：

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

故m個(gè)樣地蓄積v與q個(gè)自變量之間的關(guān)系可以表示為：

V=Xβ+eV,

E(eV)=0,

Cov(eV)=σ2I。

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

則誤差方程表示為：

eV=Xβ-V。

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

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

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

Xiqβq-vi|。

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

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

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

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

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

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

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

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

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

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

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

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

3 模型評價(jià)及精度檢驗(yàn)

3.1 不同建模方法的綜合評價(jià)

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

表2 3種估測方法建模綜合評價(jià)

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

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

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

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

圖4 各模型預(yù)報(bào)偏差相對均方根誤差

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

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

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

表3 各模型精度

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

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

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

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

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

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

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

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

式中：Py表示陰坡。

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

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

圖8 森林蓄積量估測等級(jí)分布圖

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

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

4 結(jié)論

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