雷令婷，高金萍，張曉麗，高顯連

(1.北京林業(yè)大學 林學院，北京 100083；2.國家林業(yè)和草原局 調查規(guī)劃設計院，北京 100714)

森林蓄積量反映一個國家森林資源規(guī)模，是中國森林資源調查的重要內容[1-2]，也是衡量森林質量、森林固碳作用[3-4]以及評價森林經營效果的關鍵性指標[5].遙感技術對于高效準確地獲取森林蓄積量、進行森林資源調查和監(jiān)測等具有重要意義.如何利用遙感技術結合少量地面實測數據進行蓄積量估算[6]，減少野外調查工作量，提高調查的精度和時效性，是目前森林資源調查的重點問題[7].

光學遙感數據憑借多時相、多分辨率、光譜信息豐富、獲取方便等優(yōu)勢[8]，已被廣泛應用于森林蓄積量估算[9-11].基于光學遙感數據的蓄積量估算方法主要通過遙感因子結合地面實測數據，利用多元線性回歸[12]、隨機森林[13]、k-近鄰算法（k-Nearest Neighbor,k-NN）[14]等方法建立反演模型.光學影像提供豐富的冠層信息，但因波長較短、對森林空間結構信息不敏感[15]、無法穿透樹冠且光學信號易飽和[16-17]等局限性影響蓄積量估算精度.合成孔徑雷達（Synthetic Aperture Radar,SAR）等獲取微波遙感數據的方法森林穿透性較好，更好地刻畫森林垂直結構特征、森林密度等信息[18-19]，且不受天氣影響[20]等，在森林蓄積量估算方面具有一定優(yōu)勢.近年來，學者們利用SAR 數據開展了較多森林蓄積量估算研究[21].然而，SAR 數據無法提供與植被生理相關的光譜信息，總體的估算精度并不高.此外，由于SAR 數據還受地形因素影響[22]，一定程度上影響了蓄積量估算效果.僅利用光學影像提取的森林水平結構信息或SAR 垂直結構信息進行蓄積量估算均具有局限性.若將光學和SAR 數據結合估算森林蓄積量，可以優(yōu)勢互補，具有提高估算精度的潛力.因此，研究者們開展大量研究挖掘協同森林水平與垂直結構信息進行森林蓄積量估算的方法[23].總之，協同主動、被動遙感數據進行蓄積量估算，精度具有一定程度的提高.然而，目前國內外利用光學遙感影像與SAR 數據結合的研究多集中于生物量估算和樹高估算中，對森林蓄積量估算研究相對較少.

在吉林省臨江市西小山林場內，本文以微波遙感數據的Sentinel-1 和光學數據的Sentinel-2A 為遙感數據源，通過多元逐步回歸和隨機森林方法構建森林蓄積量估算模型，并利用地面實測數據進行精度驗證，確定森林蓄積量估算最佳模型，同時分析各變量對森林蓄積量估算模型的貢獻，最終獲取西小山林場森林蓄積量的空間分布特征，可以實現森林蓄積量高效精準調查.

1 研究區(qū)與數據

1.1 研究區(qū)概況研究區(qū)位于吉林省臨江市西小山林場（127°18′8″～127°32′52″N，41°45′46″～41°36′25″E），林場總面積為218.869 km2，地處長白山腹地，鴨綠江畔.該區(qū)屬于溫帶大陸性季風氣候，年平均氣溫為2～4 ℃，年平均降水量為750～1 000 mm.西小山林場以人工林為主，樹種結構復雜，包括針葉林和闊葉林，其中針葉樹由云杉（Picea asperata Mast）、紅松（Pseudotsuga brevifolia Cheng et L.K.Fu）、臭松（Abies sibirica）、樟子松（Pinus sylvestris）等組成；闊葉樹包括白樺（Betula platyphylla）、大青楊（Populus ussuriensis Kom）和榆樹（Ulmus pumila L）等.森林覆蓋率達66.7%，成熟林面積達133.220 km2，約占林場總面積的61.11%，平均樹高約為13.3 m.

1.2 數據收集與處理

1.2.1 地面調查數據 地面調查數據主要為樣地數據和小班數據，均來源于2016—2018 年國家林業(yè)和草原局開展的東北內蒙古重點國有林區(qū)二類調查.其中，研究區(qū)按照林場面積計算該林場內布設樣地的數量，從而確定抽樣網格密度，每個網格的中心作為樣地的位置.最終確定共布設56 塊樣地（東西為30 m，南北為20 m）（圖1），記錄了環(huán)境因子數據（坡位、坡向和海拔等）、每木檢尺數據（胸徑、樹高等）以及郁閉度、林齡和株數等信息.

圖1 研究區(qū)位置圖Fig.1 Location of study area

樣地蓄積量計算是根據每木檢尺結果，計算樣地林木株數、胸徑斷面積，采用立木材積表[24]分樹種計算各林木材積，各樹種材積之和得到樣地材積，樣地材積除以樣地面積得到每公頃蓄積量.

小班數據主要包括小班面積、平均胸徑、平均蓄積量和樹種類型等，主要用于生成林場小班蓄積量分布圖.

1.2.2 影像數據 本研究影像采用歐洲“哥白尼計劃”中的Sentinel 系列衛(wèi)星數據，分別是Sentinel-1 和Sentinel-2.其中Sentinel-1 SAR 是一個基于C波段的成像系統(tǒng)，分辨率最高5 m，幅寬達到400 km[25].Sentinel-2 衛(wèi)星攜載MSI 多光譜成像儀，覆蓋13 個波段[26]，觀測幅寬達到290 km，空間分辨率分別為10、20 m 和60 m，重返周期為5 d[27].通過美國ASF 網站（https://vertex.daac.asf.alaska.edu/）下載2017 年9 月15 日sentinel-1 干涉寬模式（Interferometric Wide swath,IW）S1B GRDH 影像，極化方式為 VV、VH，空間分辨率為 5 m×20 m；在USGS 網站（https://earthexplorer.usgs.gov/）下載2017年9 月23 日Sentinel-2A 的 Level-1C 數據.由于SAR 影像數據存在因地形、斑點噪聲引起的幾何和輻射畸變，為了能夠實現地面實測數據和SAR影像精確匹配，對數據進行軌道校正、輻射定標、多視化、地形校正、分貝化和有效散射面積校正以及重采樣（采樣大小為10 m）等處理，得到后向散射系數.由于Sentinel-2A Level-1C 數據是經過幾何精校正的正射影像，并沒有進行輻射定標和大氣校正，因此利用Sen2Cor 大氣校正處理器（2.5.5 版本）進行大氣校正得到Level-2A 大氣底層反射率數據，同時采用C 模型進行地形校正，并將所有波段影響重采樣到10 m.

2 研究方法

2.1 紋理變量紋理信息可以反映森林冠層信息，避免由于樹木相互遮擋造成的影像空間色調變化，在一定程度上反映森林蓄積量信息[3].本研究利用灰度共生矩陣提取Sentinel-2A 波段的紋理信息，窗口大小設置為3×3.具體見表1.

表1 Sentinel-1 紋理信息列表Tab.1 The list of texture information from Sentinel-1

2.2 光學變量由于光學影像的波段光譜信息及其相關植被指數能夠反映植被生長狀況，可以作為反映森林蓄積量的直接變量.本研究從Sentinel-2A 影像中提取10 個波段的光譜信息和14 個植被指數，具體見表2.

表2 Sentinel-2A 特征變量列表Tab.2 The list of variables derived from Sentinel-2A

2.3 SAR 的投影角校正的后向散射系數合成孔徑雷達（SAR）數據對于森林垂直結構信息較為敏感，在一定程度上可以反映森林蓄積量狀況.然而，由于SAR 傳感器側視成像的特點，后向散射特征受到局部地形起伏的影響.因此，我們對SAR 數據進行投影角校正，并提取VV 極化方式下的變量VH 極化方式下的變量和VH/VV 的變量

2.4 森林蓄積量模型本研究采用多元逐步回歸方法，分別以提取的紋理信息、光學變量和后向散射系數以及不同類型變量組合為自變量，以實測樣地蓄積量為因變量，構建森林蓄積量估算模型.建模過程中，每個自變量隨著對回歸方程貢獻的變化，將自變量引入或剔除出回歸方程，最終在回歸方程中的自變量均為顯著的變量.此外，利用隨機森林方法構建森林蓄積量估算模型.隨機森林的流程主要為：①輸入樣本集；②隨機選擇訓練集；③從M個特征變量中選擇m個特征變量（m＜M）；④基于以上步驟，構建n個決策樹并求取平均值從而得到最終預測結果.該方法可以克服回歸方程可能造成的過擬合問題，其建立多個決策樹，并將其融合從而得到一個較為準確的模型[1].

在進行樣地調查時，由于有4 個樣地中蓄積量為零，即在這4 個樣地中沒有樹木予以剔除，以排除蓄積量為零對建模產生影響，從而導致蓄積量估算模型不準確.

為了驗證Sentinel-1 以及Sentinel-2A 影像在森林蓄積量估算方面的潛力并分析模型的估算精度，利用留一交叉驗證方法對5 個模型結果分別進行精度評價，并通過決定系數（coefficient of determination，R2）、留一交叉驗證均方根誤差(Root Mean Square Error of leave-one-out Cross-Validation,RMSECV)和相對留一交叉驗證均方根誤差(relative Root Mean Square Error of leave-one-out Cross-Validation，rRMSECV)驗證模型的預測能力.

3 結果與分析

3.1 特征變量篩選將基于Sentinel-1 和Sentinel-2A 所提取的因子分成4 組建立森林蓄積量估算模型：①僅有紋理信息；②僅有光學特征變量；③僅有后向散射系數；④3 組變量組合；對每組均采用逐步回歸方法篩選出森林蓄積量估算的最佳變量（表3），進而完成對森林蓄積量的估算.

4 組模型特征變量篩選結果如表3 所示，R2介于0.115～0.513，不同類型變量建立的森林蓄積量估算效果差異較大.其中，在紋理變量中篩選出CorB4、VarB5和VarB11與森林蓄積量相關性較高；在光學特征變量中篩選出B2、B11、PSSRa和REIP等4 個變量相關性較高；在后向散射系數中，等3 個變量相關性均較好.利用3 類變量的模型（d）具有較好的表現，R2為0.513，僅使用單獨紋理信息、光學特征變量以及后向散射系數等分別建立的模型，R2較低.

表3 特征變量篩選結果Tab.3 The results of feature variable screening

3.2 森林蓄積量建模本文將篩選的變量分別利用多元逐步回歸和隨機森林方法構建模型，如表4所示.結合3 種類型變量建立森林蓄積量模型（d）表現出明顯的優(yōu)勢，其中模型（d）的RMSECV為49.70 m3·hm-2，rRMSECV為0.26.而僅利用單一類型的變量進行森林蓄積量估算時，表現均低于模型（d），RMSECV相對較高為53.91 m3·hm-2，其中基于紋理信息模型表現相對較好.基于后向散射系數模型表現最差，RMSECV達到57.78 m3·hm-2，無法進行森林蓄積量估算.利用隨機森林的模型5 表現相對于線性模型較差，RMSECV達67.45 m3·hm-2.

表4 模型精度評價表Tab.4 The estimation accuracy of different models

利用留一交叉驗證對5 個模型估算性能進行了驗證，建立預測蓄積量與實測蓄積量的散點圖（圖2）.從圖2 可以看出，模型（d）表現相對較好，散點均勻分布在兩側，表明這個模型具有一定的適用性.而單一類型變量構建模型均相對較差，RMSECV均較高.此外，利用SAR 數據提取變量建立的模型（c）（圖2（c））由于變量信息較少，無法較好估算森林蓄積量，RMSECV相對較高.綜上說明，模型（d）相較于單一變量模型和隨機森林模型更適合研究森林蓄積量估算與區(qū)域制圖.

圖2 不同模型預測與實測蓄積量散點圖Fig.2 Scatter plot of predicted and measured stock volume from different models.

3.3 森林蓄積量空間分布比較分析5 種森林蓄積量模型可知，結合紋理變量、光學特征變量和后向散射系數構建的模型（d）在研究區(qū)森林蓄積量估算精度最高.因此，本研究以小班為單位，利用模型（d）對林場森林蓄積量進行反演，結果如圖3（a）所示，并以小班蓄積量實測圖進行對比（圖3（b））.從圖3 可以看出，西小山林場蓄積量呈現中南部較低，東部和中部偏西蓄積量較高的空間格局.小班實測平均蓄積量為147.14 m3·hm-2，預測平均蓄積量為188.31 m3·hm-2，其中，蓄積量在100～200 m3·hm-2分布最多，約占總面積的41.96%，蓄積量0～100 m3·hm-2和300～400 m3·hm-2均較少，分別占總面積的10.44%和7.55%，在預測過程中出現蓄積量值高估的現象，高估的原因可能是建模樣本數量不夠多.

圖3 西小山林場蓄積量反演圖和小班蓄積量實測圖Fig.3 The inversion map of stock volume and the measure map of stock volume in subcompartment of Xixiaoshan Forest Farm

4 討論與結論

4.1 討論本研究將Sentinel-1 和Sentinel-2A 提取的紋理信息、光學特征變量和后向散射系數作為特征變量，通過逐步回歸方法最終篩選出VarB5、VarB11、B11、PSSRa、REIP、作為建模變量，利用多元逐步回歸和隨機森林方法分別對森林蓄積量進行估算.在紋理變量中方差（VarB5和VarB11）對森林蓄積量的貢獻較大，主要反映給定窗口內光譜異質程度和像元間的相似程度，是影像灰度統(tǒng)計和空間特征信息的有效反映，可以有效地減少影像中陰影對森林參數提取的影響，同時可以提高森林與其他地物的差別[28]，從而使得紋理變量與森林蓄積量的相關性提高；短波紅外的B11 波段為蓄積量估算模型的第一貢獻者，主要原因可能為短波紅外對生物量較為敏感，而生物量與蓄積量存在一定關系[29]；PSSRa 和REIP 對于蓄積量估算的貢獻較大，主要原因可能是Sentinel-2A 數據中包含3 個與植被各類理化參數相關的紅邊波段[30]，而PSSRa和REIP 兩個參數是由紅邊波段構成[31]，可以反映森林蓄積量的變化；SAR 數據后向散射系數主要由于SAR 數據對于樹高等因子較為敏感[32-33]，而樹高與森林蓄積量關系密切[34-35]，因此也作為森林蓄積量建模的變量之一.

篩選特征變量時，利用紋理信息、光學特征變量、投影角[36]的后向散射系數單一特征變量建立模型R2較低，而基于3 類變量建立的模型（d）精度較高，這是由于森林蓄積量對森林水平結構和垂直結構信息均具有較強的依賴性，僅使用單獨結構的變量不足以充分描繪.此外，不同森林結構的特征變量為森林蓄積量估算提供不同權重信息，因此需要進一步明確各個結構特征變量對森林蓄積量的估算能力的影響，從而提高森林蓄積量估算精度.

對比5 個模型發(fā)現，模型（a）篩選的紋理變量能反映影像中森林紋理信息的分布和變化情況；模型（b）篩選出的變量能夠較為全面地代表植被信息，起到圖像增強的作用；模型（c）變量僅反映森林垂直信息.總之，以上模型篩選出的變量較為單一.利用不同類型的變量建立的模型（d），在一定程度上可以較為完整地反映森林蓄積量信息，同時較好地抑制了僅利用光學遙感和微波遙感所造成的蓄積量估算不準確的問題.基于光學和微波遙感數據利用隨機森林方法進行森林蓄積量估算效果相比于多元線性回歸方法較差.盡管目前利用機器學習等方法構建的非線性模型應用較多，但是機器學習給人的感覺就像一個黑盒子，我們無法控制模型內部的運行，只能在不同的參數和隨機種子之間進行嘗試.同時對因變量的解釋性較低，即無法直接解釋測量變量與蓄積量估算模型的關系以及各變量對模型的貢獻[37].

值得注意的是，在本研究中存在蓄積量高估的情況，主要原因可能是樣地數據相對較少，且56 個樣地實測蓄積量值多集中于100 m3·hm-2以上，在森林蓄積量建模時未考慮蓄積量較小值，影響了森林蓄積量精度.目前，廣義神經網絡和深度學習在森林蓄積量估算中的應用取得較好效果，在今后的研究中將嘗試利用此類方法并與線性模型對比，以實現快速準確的森林蓄積量動態(tài)監(jiān)測和定量評價.

4.2 結論

（1）通過Sentine-1 和Sentinel-2A 為數據源提取的80 個紋理變量、24 個光學特征變量和3 個后向散射系數3 種類型變量，CorB4、VarB5、VarB11、B11、PSSRa、REIP、等7 個變量對蓄積量估算模型貢獻較大.

（2）5 組森林蓄積量估算模型中，結合不同類型變量構建的模型（d）效果最好，RMSECV為49.70 m3·hm-2，森林蓄積量估算精度達74%.

（3）西小山林場森林蓄積量整體呈現中南部較低，東部和中部偏西蓄積量較高的空間格局，預測平均蓄積量達到188.31 m3·hm-2.