基于Sentinel-2光譜與地形特征的山區森林分類
——以武夷山國家公園為例

張春瑩，江 洪，，，林敬蘭，岳 輝

（1.福州大學 空間數據挖掘和信息共享教育部重點實驗室，衛星空間信息技術綜合應用國家地方聯合工程研究中心，數字中國研究院（福建），福建 福州 350108；2.福建省水土保持試驗站，福建 福州 350002；3.長汀縣水土保持中心，福建 長汀 366300）

森林是全球陸地生態系統的重要組成部分［1］.傳統的森林生態系統調查依靠實地調查來收集相關信息，但是需要大量的人力、物力和財力，并且容易受到區域的限制［1-2］.遙感技術在林業資源管理中具有監測范圍廣、成像速度快、重訪周期短、數據成本低等優勢［3-4］.傳統基于像元的分類方法僅利用了地物的光譜信息，因此難以區分光譜信息相似的地物，而面向對象分類方法的基本處理單元是對象，可在此基礎上充分利用不同對象的光譜、空間和紋理等相關特征［5］.光學圖像中豐富的光譜變量是遙感圖像分類中最重要的變量，研究表明，近紅外和紅邊波段是植被較為敏感的波段［6］.例如，李丹等［7］通過采用WorldView-2 影像來識別典型喬木樹種，研究發現紅邊波段明顯改善了識別結果；Tigges 等［8］研究發現RapidEye 影像的紅邊波段能顯著提升植被的分類結果的準確性.

歐洲空間局（European Space Agency，ESA）分別于2015 年6 月和2017 年3 月成功發射了Sentinel-2 系列衛星，該衛星為森林植被的識別提供了新的數據源，其搭載的多光譜儀（MSI）已被應用于樹種分類實驗中［9-10］.例如，于婉婉等［11］研究發現不同生長期的優勢樹種在Sentinel-2數據的紅邊波段表現有不同差異；黃雙燕等［12］利用Sentinel-2影像提取干旱區農作物的分類信息，研究發現紅邊特征對于識別不同作物間的物候差異表現更為敏感.此外，在地形復雜的山區，山區地物通常形成地帶性分布，輔助以地形因子進行地類信息提取可顯著提升分類精度.例如，陳元鵬等［13］利用GF-1影像在山地工礦復墾區進行地類信息提取時加入地形因子得到了較高精度提取結果；Dubeau等［14］利用Landsat影像光學、雷達和地形變量的組合方式提取達布斯濕地信息，分類效果較好.考慮到植被生長茂密的區域多位于山地丘陵區，在進行地物提取時復雜的地形條件會對結果造成較大干擾，導致在開展森林植被信息提取時精度難以滿足林業資源管理的需求.因此，采用Sentinel-2 影像的紅邊光譜特征輔助地形因子開展山區森林分類研究很有必要.鑒于此，筆者以武夷山國家公園為研究區，以較高分辨率的Sentinel-2 影像和ALOS PALSAR 數據為數據源，采用多尺度分割方法對多光譜影像進行分割以構建影像對象，提取27個光譜特征變量并計算其最優特征空間，將地形信息與最優特征空間相結合組成3種不同的特征變量子集，在此基礎上對武夷山國家公園進行森林分類，以期為地形復雜的山區高精度森林分類提供有效的技術手段.

1 研究區概況與數據源

1.1 研究區概況武夷山國家公園（117 °24 ＇13 ″～117 °59 ＇19 ″ E，27 °31 ＇20 ″～27 °55 ＇49 ″ N）（如圖1所示），總規劃面積1 001.41 km2，屬于中亞熱帶季風氣候區，地貌以山地和丘陵為主，海拔高度差異懸殊，最高處海拔為2 160.8 m，最低處海拔為176.1 m，植被類型多樣，森林覆蓋率為93.74%.

圖1 研究區示意圖

1.2 數據源與預處理

1.2.1 數據源本文所用數據包括Sentinel-2 數據、ALOS PALSAR 高程數據、植被分布圖，Sentinel-2 數據下載于歐洲航天局（ESA）官方網站（http：//scihub.copernicus.eu/），數據獲取時間為2019 年9 月24 日，ALOS PALSAR高程數據下載于美國阿拉斯加衛星設備分布式活動檔案中心（https：//search.asf.Alaska.edu/#/）.植被分布圖來自國家林業和草原局昆明勘察設計院.所采用的Sentinel-2相關波段參數如表1所示.

表1 Sentinel-2波段介紹

1.2.2 樣本數據通過對武夷山國家公園植被分布圖的讀取，參考國家林業和草原局出臺的《國家森林資源連續清查技術規定 2014》，并結合研究需求，得到本文實驗區所涵蓋的地物類別，包括闊葉林、針葉林、竹林、農業用地、建設用地、水域.通過Google Earth 高分辨率影像（http：//earth.google.com/）選取樣本點，遵循樣本點的選取盡可能均勻分布的原則，樣本點的空間分布如圖2 所示，共選取750 個樣本點，158個針葉林樣本，150個闊葉林樣本，183個竹林樣本，62個農業用地樣本，91個建設用地樣本和61個水域樣本，其中訓練樣本533個，驗證樣本217個.

圖2 樣本點空間分布

1.2.3 數據預處理筆者使用的是Sentinel-2 影像中除沿海大氣氣溶膠波段和卷積云波段以外的10 個波段，首先采用歐空局提供的Sen2Cor插件將下載的數據級別由L1C級轉化為L2A級，即將10個波段的大氣表觀反射率處理成大氣底層反射率，然后在SNAP軟件中全部重采樣為10 m，最終得到實驗所需的包含10個波段的L2A級Sentinel-2影像數據.

2 研究方法

采用面向對象的方法進行分類，影像分割是面向對象分類方法中的重要步驟，分割是從子像素開始，從上到下合并區域，直到滿足異質性標準［15］，特征提取是利用影像的光譜和空間特征等來提取分割對象所需要的類別信息.基于eCognition9.0遙感處理平臺對預處理后的影像進行面向對象分類實驗.為了檢驗Sentinel-2 光譜特征和地形因子對森林分類的能力和貢獻，構建了3 個特征組合方案：方案1：參考光譜指數；方案2：參考光譜指數+紅邊指數；方案3：參考光譜指數+紅邊指數+地形因子.3 個特征組合方案作為隨機森林分類器的輸入特征，比較不同特征組合方案在分類精度提升上的表現，技術流程如圖3所示.

圖3 技術流程

2.1 影像分割與最優分割尺度評價影像分割是面向對象分類的一個重要部分，分割尺度是決定分割效果的關鍵.選用基于可見光和近紅外波段的面向對象多尺度分割方法進行影像分割.為避免由于分割尺度不合適所出現的欠分割或過分割的弊端，采用ESP2（Estimation Scale Parameter）工具輔助目視解譯來確定影像的最優分割尺度.ESP2 來源于Dragut 等［16］提出的ESP 算法，該工具計算影像對象異質性局部方差的變化率，生成局部方差變化率曲線，當曲線出現峰值時，表示該點變化率最大，則該點相對應的分割尺度為最優分割尺度.

影像多尺度分割在可見光和近紅外4 個波段上進行，經過多次試驗，最終形狀因子確定為0.2，緊致度因子確定為0.5，在此基礎上，使用ESP2工具計算最優分割尺度，將分割尺度的步長設為1，起始分割尺度設為30，得到ESP 分割尺度評價結果（如圖4 所示），變化率曲線有多個峰值，表示影像的最優分割尺度參數有多個，選擇分割尺度40、50、89、111進行分割尺度效果的目視解譯對比.

圖4 ESP分割尺度評價

分割尺度為40、50、89、111 下的分割效果如圖5 所示，圖5a 和e 表示分割尺度為40 時，闊葉林、針葉林、農業用地等分布區域面積大的地物類型被分割得較為細碎，產生了過分割現象，不滿足實際地物的分割需求；如果分割尺度過大，小面積地物會與大面積地物混雜在一起，導致分割地物的混淆，如圖5c～d 和g～h，當分割尺度為89和111時，同一對象內包含了不同的分割地物，6種地物類型都表現出了不同程度的欠分割現象.而在最佳分割尺度下，既合并了同類地物又與其他地物明顯區分開來，地物分割都達到了較好的分割效果，如圖5b 和f，當分割尺度為50 時，各地物類型分割效果較好，分割對象符合地物實際空間分布邊界，可在此基礎上進行后續實驗.因此，將最優分割尺度選定為50.

圖5 不同分割尺度下影像分割效果

2.2 Sentinel-2光譜特征提取根據參考的相關文獻［9］，［17-26］，構建27個光譜特征，包含Sentinel-2影像的10個原始波段和17個光譜指數（如表2所示），其中包括11個紅邊指數.

表2 基于Sentinel-2影像的光譜指數

2.3 地形因子提取地形因子中包含豐富的地形表面形態特征信息，將各種地形因子結合在一起，能夠較好地刻畫地表起伏的變化.因此，根據本文需求，將高程、坡度和坡向作為地形因子變量，ALOS PALSAR 數據即DEM（Digital Elevation Model）數字高程模型含有該區域海拔高度信息；坡度、坡向特征可使用ArcMap軟件從DEM中計算獲得.

2.4 特征優化與分類隨機森林（Random Forest，RF）算法是Leo 等［27］于2001 年提出的由多個決策樹分類器集成的一種機器學習算法.采用RF 模型用來作為分類與計算和選擇特征重要性的工具.該算法采用自助法重采樣技術，將原始樣本中有放回地隨機抽取2/3的數據作為訓練集，訓練集中的每個樣本都由多個決策樹分類，產生多個分類結果，投票決定每個樣本的分類類別.未被抽取的樣本稱為袋外數據（Out of Bag，OOB），OOB誤差越低，RF模型準確性越高.

RF 算法在完成特定分類任務的同時，還可評估模型中每個特征變量的重要性.由于Sentinel-2 影像的原始光譜波段豐富，在后續研究中眾多學者研究出了相關光譜指數，面向對象的分類方法可輸入眾多特征參與分類，但較多特征會導致信息的冗余，降低分類的速度和精度.因此，為充分利用影像對象豐富的光譜與空間特征，采用RF算法中的基尼系數（Gini）評估各特征的重要性，確定最優特征空間.Gini系數數值越大，則表明該特征的重要性也越高.最終根據重要性排序以及OOB 誤差進行特征優選，確定最終的特征變量子集.

2.5 精度評價本文中3 種特征組合方案的森林分類精度采用混淆矩陣進行評價，評價指標有總體精度（Overall Accuracy，OA）、Kappa 系數、生產者精度（Producer’s Accuracy，PA）、用戶精度（User’s Accuracy，UA）以及調和平均值（F1）.其中，F1的計算公式［28］為

3 結果與分析

3.1 特征重要性排序及評價采用RF Gini 系數計算得到Sentinel-2 影像27 個光譜特征的重要性排序（如圖6 所示）.其中重要性排序前五的特征是LSWI、B2（藍波段）、B4（紅波段）、B12（短波紅外2 波段）和Ndre1.重要性排名最高的特征為LSWI，為11.79%，該特征由近紅外波段B8 與短波紅外1 波段B11 計算得到.其次為藍波段B2和紅波段B4，重要性分別為10.35%和7.50%，其中，藍波段B2的重要性得分遠高于其他波段，可能受到其他相關特征的影響［29］.B12、Ndre1和GNDVI的重要性較高，分別為4.48%、4.37%和4.30%，由于研究區植被覆蓋度高，植被在此處具有較強的吸收性，所以紅光波段達到貢獻度高.貢獻度最低的為NDVIre3，重要性得分僅有0.81%.

圖6 特征變量重要性排序

通過改變特征數量得到27 個光譜特征數量與分類精度的關系如圖7a 所示，前期精度值隨特征數的增多快速增加，此后持續保持上下浮動的趨勢，當特征數為17 時，精度達到最大值0.911 0.因此，在森林分類中選擇重要性排序較高的前17個光譜特征構建最優特征空間.對這17個特征中所包含的Sentinel-2原始波段數量進行統計（如圖7（b）所示），數量最多的是B4（紅波段），8個；其次是B8（近紅外波段），5個；B5（紅邊1 波段），5 個.此外，B6（紅邊2 波段）和B7（紅邊3 波段）也較多，進一步表明紅邊波段在森林分類中具有較高的貢獻度.

圖7 特征數量與分類精度的關系和重要特征波段統計

3.2 分類結果與精度評價采用面向對象多尺度分割方法結合不同特征組合方案進行森林分類，如圖8所示.在方案1中僅利用參考光譜特征分類（圖8b）時，“椒鹽現象”十分明顯；在方案2中加入紅邊特征參與分類（圖8c）后效果較好，各地類的圖斑破碎度有所降低；在方案3 中將光譜特征與地形因子結合進行分類（圖8d）后效果最好，各地類“椒鹽現象”得到改善.此外，將分類結果與Sentinel-2影像對比可以看出，闊葉林、針葉林與竹林交錯分布，闊葉林分布面積最大，最連續；竹林與闊葉林鑲嵌分布，且與農業用地交錯分布；針葉林分布較為集中.

圖8 分類結果

不同特征組合方案的分類結果精度對比如圖9 所示.方案1、2 和3 的總體精度（OA）與Kappa 系數分別為88.13%、0.856 4，89.50%、0.871 0 和90.87%、0.887 8.對比分析3 種特征組合的森林分類精度，僅使用參考光譜特征分類時，闊葉林、竹林和農業用地的生產者精度（PA）較低，三者的F1分別為0.81、0.88和0.67.加入紅邊特征進行分類后，闊葉林和竹林的生產者精度（PA）分別提升至0.88和0.90，其F1都提升了0.03，闊葉林和農業用地的用戶精度（UA）分別從方案1 的0.76 和0.61 提升至0.79 和0.70.在光譜特征基礎上加入地形特征參與森林分類時，竹林和農業用地的生產者精度（PA）分別提升至0.92和0.77，相對于方案1來說分別提升了0.08和0.05，相對于方案2來說提升了0.02和0.15，F1分別提高至0.92和0.79，比只使用參考光譜指數提升了0.03 和0.12，闊葉林和農業用地的用戶精度（UA）提升至0.82 和0.81，相對于方案1 提升了0.06 和0.20.由此可見，方案3 的總體分類精度（OA）和Kappa 系數最高，表明將Sentinel-2的光譜特征與研究區地形因子相結合作為分類器的輸入變量參與分類時，能在一定程度上提高森林分類精度，得到最佳的分類效果.

圖9 不同特征組合的分類精度對比

不同特征組合分類混淆矩陣如圖10所示.使用方案1進行分類的效果較差，如圖10a所示，特別是農業用地與竹林和闊葉林之間存在一定程度上的地物混淆，導致分類結果整體來看效果欠佳，這是因為竹林、闊葉林和農業用地在Sentinel-2 影像上視覺差較小，因此僅利用參考光譜指數很難將三者區分開.如圖10b所示，在方案2中加入紅邊特征后，竹林精度有了明顯提高，說明竹林對于紅邊數據敏感度更高.方案3中加入地形因子后，如圖10c所示，闊葉林與農業用地2類地物相比方案1和方案2來說可更為有效的區分開，說明地物的分布與地形有一定程度上的相關性.

圖10 不同特征組合的分類混淆矩陣

從3 種分類特征組合結果可見，方案3 使用Sentinel-2 的光譜特征輔助地形信息的總體分類精度（OA）、Kappa 系數和F1最高，其次為方案2，方案1 的分類精度最低，效果最不好.闊葉林、針葉林、竹林和農業用地在方案3中取得較好的效果.由此可見，由于植被對紅邊波段的敏感度更高，紅邊波段對于森林分類也顯得更為重要，此外，加入地形因子可提高地形復雜區域地類信息提取的精度，這與陳元鵬等［13］、Dubeau等［14］的研究結論相同.

3.3 分類結果與地形因子分異統計分析為分析各地類在高程、坡向和坡度上的分布特點，通過自然間斷法對高程、坡向、坡度3種地形因子重分類［30］，利用ArcMap分區統計工具統計不同地形因子等級內森林類型面積所占比例，不同森林類型在高程、坡向、坡度上的統計結果如圖11所示.各森林類型面積比例隨著高程的不同而出現明顯差異如圖11a 所示，其中，針葉林的面積比例隨著高程的升高波動起伏較大，在1 427～1 661 m 和1 662～2 162 m 范圍內占比分別為63%和57%；農業用地、水域和建設用地在178～446 m范圍內占比為46%.各森林類型面積比例在不同的坡向上也有明顯的地形響應如圖11b所示，其中，農業用地、水域和建設用地在平面上面積占比為79%，闊葉林在東坡、東南坡和南坡的分布面積比例遠高于其他森林類型.各個森林類型的面積比例隨著坡度的升高而出現明顯且較為規律的變化如圖11c 所示，闊葉林、竹林和針葉林的面積比例隨著坡度的增大逐漸升高，農業用地、建設用地和水域隨著坡向的增加面積占比逐漸減小.由此可見，不同森林類型的分布與地形因子之間存在一定的規律性，所以在進行山區森林分類時加入地形因子是提高分類精度的一種有效手段.

圖11 分類結果與地形因子統計

4 結束語

利用多尺度分割的方法對Sentinel-2 影像進行分割以構建影像對象，在此基礎上提取影像的光譜特征，將其與研究區地形因子相結合，構建3種特征組合方案對武夷山國家公園進行森林分類.結果表明：

1）為避免特征數量過多導致的信息冗余，采用RF Gini系數對特征空間進行重要性評估，在Sentinel-2影像上重要性排名前17 的最優特征空間中，紅波段、近紅外波段和新增的紅邊波段在森林分類中的重要性較高；

2）不同特征組合方案參與分類時，光譜特征與地形因子組合變量的分類方案精度最高，總體分類精度為90.87%，Kappa系數為0.887 8.因此，采用Sentinel-2紅邊光譜特征結合地形因子的多特征組合方法是提高山區森林分類精度的一種有效手段.

隨著遙感影像空間、光譜及時間分辨率的提高，充分利用Sentinel-2 影像紅邊波段光譜信息并將其與地形信息相結合的方式是實現森林高精度分類的一種有效手段.但筆者僅討論了單一時相下將Sentinel-2影像的光譜特征和地形特征相結合的森林分類方法，下一步研究如何融合多時相遙感影像，將其應用于山區復雜地形條件下的森林制圖.