基于Sentinel-2時序多特征的植被分類

郭文婷， 張曉麗

（1.北京林業大學 省部共建森林培育與保護教育部重點實驗室，北京100083；2.北京林業大學 精準林業北京市重點實驗室，北京100083）

森林資源調查可以統計不同林地的生物學特性、生態特征和經營價值，為森林的合理經營規劃提供有效的數據支持［1-2］。由于中國各地區之間自然條件差異較大，森林類型多種多樣，林分結構較為復雜，森林資源調查工作難度很大［3］。植被分類是研究森林資源狀況和動態變化規律的基礎，利用遙感手段可以更加快速、準確地識別植被類型［4］。因此傳統的依靠人力的地面調查已經逐漸轉變為基于遙感手段的植被識別和植被分類［5-6］。以往許多研究學者基于植被的光譜特征進行了植被分類［7-8］。由于許多地物光譜特征相似度高及環境因素的影響，許多地物難以分離。針對這一問題，有學者采用神經網絡［9-10］、支持向量機［11-12］、隨機森林［13-14］、決策樹等［15-17］分類方法對植被進行分類并取得較好的分類結果。但是植被分布雜亂無序，且存在 “同物異譜”和 “同譜異物”的現象［18］，只利用光譜特征進行分類很難取得高精度的分類結果。有的學者采用結合紋理特征進行分類的方法并得到了更高的分類精度，證明了紋理特征有助于提高分類精度［19］。近年來，為了進一步提高植被的分類精度，有學者對植被的物候信息進行研究，發現有些植被的物候特征差異很大，利用多時相的遙感數據可以提高物候差異明顯的植被區分能力［20］。不同植被的物候特征、紋理特征、光譜特征各不相同，因此其歸一化植被指數（normalized difference vegetation index，NDVI）時間序列、紋理和光譜響應存在一定的差異。利用單一特征進行植被分類，精度往往較低。由于過去遙感數據源的限制，數據很難同時滿足高時間分辨率、高空間分辨率和豐富的光譜信息這3個特點，例如MODIS數據具有較高的時間分辨率和豐富的光譜信息，但空間分辨率較低，只適合大尺度范圍的植被分類［21］；Landsat系列數據的空間分辨率也較低，且其光譜信息中缺少植被十分敏感的紅邊波段［22］；而Sentinel-2數據則具有較高的時間分辨率、空間分辨率和豐富的光譜信息［23-24］，尤其是其多個紅邊波段對植被區分具有很大幫助。由于不同植被的物候期存在一定的差異，因此其NDVI時間序列變化特征也存在一定差異，利用多期遙感影像構建NDVI時間序列有助于將物候期差異明顯的植被區分。本研究基于多時相的Sentinel-2遙感影像，首先選取植被生長旺盛時期的影像，計算NDVI并設定合適的閾值將研究區內的植被提取出來，剩余部分歸為非植被。然后選取NDVI時間序列、最佳時相的Sentinel-2數據中10個波段的光譜反射率特征和主成分分析前3個分量的紋理特征作為分類特征，將研究區內的植被類型分為耕地、草地、常綠針葉林、落葉針葉林和落葉闊葉林五大類。

1 研究區概況

旺業甸實驗林場（圖1）位于內蒙古赤峰市喀喇沁旗西南部，地屬燕山山脈北麓，屬茅荊壩次生林區的一部分，為淺山丘陵地貌，地勢西南高東北低，一般坡度為10°～30°，平均海拔為1300 m，最高峰是翠云峰， 海拔為 1890.9 m。41°21′～41°39′N， 118°09′～118°30′E，東西長 41 km，南北寬 20 km。林場總面積為253.07 km2，有林地面積為219.6 km2，林木總蓄積量為157萬m3，森林覆被率為93.24%。在有林地面積中：天然次生林為115.33 km2，人工林為104.27 km2。其中人工落葉松Larix gmelinii面積為 45.33 km2， 油松Pinustabuliformis為 56.47 km2， 樟子松Pinus sylvestrisvar.mongholica為 2.15 km2，云杉Picea asperata為0.33 km2。旺業甸林區的生物資源復雜多樣，具有很強的代表性、典型性、特有性及脆弱性，在全國的多樣性保護中占有重要地位。林場有林地有落葉松，云杉，樟子松，紅松Pinuskoraiensis，油松，黑樺Betula dahurica，白樺Betula platyphylla，山楊Populus davidiana，核桃楸Juglansmandshurica等12個優勢樹種。

2 研究數據

2.1 遙感數據與預處理

使用的Sentinel-2數據是從歐洲航天局（European Space Agency，ESA）的數據共享網站（https://scihub.copernicus.eu/dhus/#/home）免費下載。Sentinel-2A是歐洲航天局（ESA）于2015年6月23日發射的環境監測衛星，重訪周期為10 d，其雙星Sentinel-2B于2016年發射后，將重訪周期縮短到5 d。Sentinel-2A數據的具體波段信息見表1。由于植被自身的物候特征，植被的光譜反射率在時間序列上有一定的變化規律，因此，在2017年內具有可用影像的月份均選取一景云量較少、質量較高的Sentinel-2影像，共9景，成像時間分別 為 2017-01-15， 2017-02-14， 2017-04-15， 2017-05-25， 2017-06-14，2017-09-22， 2017-10-27， 2017-11-21，2017-12-31。 利用 ESA提供的SNAP軟件對數據進行輻射定標、大氣校正、幾何校正和圖像裁剪。將像元較大的波段重采樣為10 m×10 m的像元，從而使各個波段的像元大小一致。

圖1 研究區位置示意圖Figure 1 Location of the study area

表1 Sentinel-2波段信息表Table 1 Band information of Sentinel-2

2.2 外業數據

本研究中的外業數據來源于2017年9月研究區的實地調查，共選取300個適合本次研究的樣點數據。實地調查遵循樣本具有代表性、樣本均勻分布的原則，在研究區內選取具有代表性的植被類型進行采樣，且保證采集的樣本點均勻分布在整個研究區內。樣本點的數據屬性包括樹種類別和用手持GPS定位的經緯度坐標。由于采集得到的每個類別樣本點數量不同，因此分別選取每個類別的訓練樣本，使其能夠均勻分布于研究區內。共選取訓練樣本數據122個，其中非植被為18個，耕地20個，常綠針葉林22個，落葉闊葉林24個，落葉針葉林22個，草地16個。剩余178個樣點數據作為驗證樣本，其中非植被為19個，耕地13個，常綠針葉林46個，落葉闊葉林41個，落葉針葉林43個，草地16個。

3 研究方法

采用以下3種方法對研究區進行植被分類：①最大似然法，以外業采集數據作為樣本，基于Sentinel-2數據的 10個波段（B02， B03， B04， B05， B06， B07， B08， B08a， B11， B12）進行分類， 且各個類別的樣本數據光譜特征均服從正態分布。②根據不同地物的歸一化植被指數（NDVI）時序變化特征和光譜特征（NDVI時序＋光譜特征），選取最佳分類閾值，建立基于決策樹的分層分類模型，將研究區分為非植被、耕地、草地、常綠針葉林、落葉針葉林、落葉闊葉林。③NDVI時序特征＋光譜特征＋紋理特征（Sentinel-2時序多特征），采用支持向量機分類器，經過多次實驗核函數選取高斯函數，核函數系數γ為0.2，懲罰因子C為7，對在植被生長旺盛期，利用NDVI閾值提取出的植被進行植被類型分類，將其分為耕地、草地、常綠針葉林、落葉針葉林、落葉闊葉林。

NDVI對植被反應敏感，能夠反映出植被與其他地物在可見光、近紅外波段的明顯差異。因此，首先選取植被生長茂盛時期的影像，計算其NDVI，設置合適的閾值將植被提取出來。NDVI計算公式如下： INDV=（B08-B04）/（B08＋B04）。 其中： INDV為歸一化植被指數， B08為近紅外波段的反射率，B04為紅光波段的反射率。

3.1 構建NDVI時間序列

不同時間的NDVI能夠反映植被的物候信息，且兩者有一定的規律性，是區分植被類型的重要特征。例如，常綠針葉林常年保持綠色，NDVI在1 a中的變化較小；農作物在6-9月為生長期，此時耕地的NDVI也呈增長趨勢，到9月NDVI達到最大值；落葉類植被一般在5月開始長葉10月開始落葉，NDVI也隨之在5月開始上升10月開始下降。因此，分別計算9個時相影像的NDVI，構建6種地物的NDVI時序曲線，描述6種地物的動態變化過程。

3.2 最佳時相光譜反射率特征分析

根據地物的物候特征可知：在5-9月均為5種地物的最佳時相，因此選取這段時間內質量最好的2017-06-14影像繪制光譜反射率曲線進行單波段光譜分析。對于NDVI時間序列曲線相似、難以區分的地物，利用不同地物的光譜反射率差異可以增加這種相似地物的區分度，但是仍可能存在混分現象。例如落葉針葉林和草地，兩者的NDVI時序變化規律相似，光譜反射率的變化規律也相似。因此，對紋理特征進行分析后，在最終分類時應用紋理特征以提高分類精度。

3.3 紋理特征分析

紋理特征是一種全局特征，它描述了物體表面的特性。與基于像素點的顏色特征不同，紋理特征需要在包含多個像素點的像素區域中進行統計計算。而且紋理特征常具有旋轉不變性，對于噪聲也具有較強的抵抗能力。本研究共選取了Sentinel-2數據中的10個波段進行分類研究，若提取每個波段的紋理特征并參與最終分類，將會產生大量數據，計算量太大。因此，在提取紋理特征前首先對數據進行主成分分析處理，然后選取前3個主成分分量用來進行后續的特征提取。1973年HARALICK提出的灰度共生矩陣是遙感圖像中常采用的紋理特征提取方法［25-26］。灰度共生矩陣是通過研究圖像上像元灰度的空間相關特性來描述紋理。通過研究分析選取5×5大小的移動窗口，利用灰度共生矩陣分別計算前3個分量的8種紋理特征，即均值、方差、同質性、對比度、差異性、熵值、二階矩、相關性。

綜上所述，共選取 43個特征用于分類，包括9個時相（2017-01-15，2017-02-14，2017-04-15，2017-05-25，2017-06-14，2017-09-22， 2017-10-27， 2017-11-21， 2017-12-31）的 NDVI特征，10 個波段（B02，B03，B04，B05，B06，B07，B08，B08a，B11，B12）的光譜反射率特征和24個紋理特征。

4 結果與分析

4.1 特征提取結果

由圖2分析可得：①常綠針葉林在1 a中都是綠色植被，NDVI變化相對較小。在1，2，3，11，12月其他植物多處于非生長期，常綠針葉林的NDVI明顯高于其他地物，容易被區分；②研究區內的耕地在5-9月屬于生長期，NDVI處于上升趨勢，但同一時期的其他植物也處于生長旺盛期，5種地物的NDVI均較高，相比之下耕地的NDVI值反而是最低的，容易被區分；③落葉闊葉林在5-9月屬于1 a中生長最旺盛的時期，其NDVI值達0.80以上，高于同一時期的其他地物，可以進行區分。

通過分析圖3發現：①耕地在B11波段的反射率上升相比于其他地物的反射率下降具有顯著差異，利用B11波段很容易將耕地與其他地物區分出來；②草地在B06，B07，B08，B08a波段的反射率與其他地物均存在差異，可以利用這幾個波段將草地與其他地物區分；③落葉針葉林在B05波段的反射率低于落葉闊葉林，而在B06波段的反射率卻高于落葉闊葉林，可以利用這一特征將兩者進行區分。

4.2 分層分類閾值設置

通過分析典型地物的NDVI時間序列曲線與典型地物的光譜反射率曲線，得出結合NDVI時序與光譜特征分類方法的最佳分類閾值。在2017-09-22的影像中，INDV＜0.37為非植被；在2017-12-31的影像中，INDV＞0.39為常綠針葉林；在2017-06-14的影像中，B11波段的反射率B11＞25%為耕地（圖3）；在2017-06-14的影像中，B08波段的反射率B08＞21%為草地；在2017-05-25的影像中，INDV＞0.8為落葉闊葉林；在2017-05-25影像中，0.7＜INDV＜0.8為落葉針葉林。

圖2 典型地物的NDVI時間序列曲線Figure 2 NDVI temporal profiles of typical land cover types

圖3 典型地物的光譜反射率Figure 3 Spectral reflectance curve of typical land cover types

4.3 分類結果

通過分析比較3種方法的分類結果（圖4和表2）發現：①最大似然法。非植被與耕地在分類結果中的面積明顯偏少，存在大量的漏分；草地在分類結果圖中沒有明顯的表現出來，也存在嚴重的漏分；常綠針葉林、落葉針葉林和落葉闊葉林在分類結果完全是分散分布且交錯分布的，這與實地調查情況不符，存在明顯的混分。②NDVI時序+光譜特征方法。非植被分類情況與基于多特征的分類結果一致，且明顯可以看出包含了一些道路，與實際調查情況相符；耕地和草地的分類結果與其他2種方法相比面積較多，尤其是草地的面積過多，與常綠針葉林、落葉針葉林和落葉闊葉林產生混分。③NDVI時序+光譜特征+紋理特征方法。非植被與耕地的面積和分布情況均符合實際情況；常綠針葉林主要分布在耕地周邊，集中在研究區的中部；落葉闊葉林主要分布在研究區的西部和西南位置，少數分散分布于整個研究區，與落葉針葉林交錯分布；草地在研究內呈現為一些散碎的斑塊，且面積相對較小。整體的分類結果與實地調查情況有較高的一致性。

4.4 精度評價

圖4 3種方法的分類結果示意圖Figure 4 Classification results of three methods

利用實地調查采集的178個樣點數據作為驗證樣本，其中非植被19個，耕地13個，常綠針葉林46個，落葉闊葉林41個，落葉針葉林43個，草地16個。對比3種方法的分類結果與實地調查結果，并計算出每種方法的類別分類精度及總體分類精度。由表3可見：基于NDVI時序+光譜特征+紋理特征的多特征分類總體精度達87.64%，分別比最大似然法和NDVI時序＋光譜特征的分類總體精度提高了15.73%和14.61%。而且基于NDVI時序+光譜特征+紋理特征的多特征分類中單個類別的分類精度也均較高，其中常綠針葉林的分類精度最高達到了95.65%，這是由于利用了NDVI時序特征，可以通過分析不同植物在1 a中的NDVI變化趨勢，選擇出最佳的分類時相將常綠針葉林與其他地物區分。耕地的分類精度也比較高，其NDVI在1 a中的變化具有一定的規律性，且耕地在B11波段的反射率與其他地物有明顯差異，利用這些特征可與其他地物很好地區分。NDVI時序＋光譜特征的分類精度與NDVI時序+光譜特征+紋理特征的分類精度相比，非植被這一類別的分類精度相同，但其他類別的分類精度則較低。這是由于NDVI時序變化和光譜特征代表的僅僅是一個平均值，并不能完全代表地物的特征變化，因此，根據平均值來選取閾值分層分類的精度較低。

表2 3種分類方法的地類面積Table 2 Land area of three methods

表3 3種分類方法的精度評價Table 3 Accuracy evaluation of three methods

5 結論與討論

本研究基于多時相的Sentinel-2遙感影像，選取NDVI時間序列、最佳時相的Sentinel-2數據中10個波段的光譜反射率特征和主成分分析前3個分量的紋理特征作為分類特征，利用支持向量機分類器，將研究區內的植被類型分為耕地、草地、常綠針葉林、落葉針葉林和落葉闊葉林。得出結論如下：①基于多特征的分類總體精度達87.64%，Kappa系數為0.85，分別比最大似然法的分類結果提高了15.73%和0.20。因此，結合多種特征的分類方法有助于提高分類精度。②通過構建NDVI時間序列，分析地物的NDVI變化趨勢，選擇常綠針葉林提取的最佳時相進行提取，可以大大提高常綠針葉林的分類精度。NDVI的時序特征對于植被的區分具有很大幫助。③耕地在Sentinel-2影像的B11波段上的反射率與其他地物有顯著差異，而且耕地在NDVI時間序列上也與其他地物有明顯差異，利用這些特征極大提高了耕地的分類精度。④采用分層分類的思想，首先將研究區內的植被提取出來，然后再對植被進行類型識別，可以排除非植被因素的干擾，有效提高植被類型的分類精度。

分類結果中常綠針葉林與耕地均得到了很高的分類精度，這是由于本研究考慮了植被的物候信息，構建了不同植被的NDVI時間序列變化曲線，并充分利用植被的物候特征進行分類。分類結果中也顯示草地的分類精度是所有類別中精度最低的。一方面原因是，在研究區內草地分布比較分散而且比較破碎，在實地調查中采集的草地樣本點相比其他類別樣本數量較少，分類時訓練樣本的數量相比其他類別也比較少，從而導致了分類精度較低。另一方面，草地的NDVI時間序列變化規律與落葉針葉林、落葉闊葉林相似，其光譜反射率的變化規律與落葉針葉林也相似，從而導致草地極易與落葉針葉林、落葉闊葉林混分，使其分類精度較低。本研究只選取了NDVI時序特征、光譜特征和紋理特征3種特征作為分類特征，而地形特征中也包含了豐富的信息，可以充分利用地形特征進行深入研究。