多時相Sentinel-2影像在浙西北茶園信息提取中的應用

李龍偉， 李 楠， 陸燈盛

（1.浙江農林大學 浙江省森林生態系統碳循環與固碳減排重點實驗室，浙江 杭州311300；2.浙江農林大學 環境與資源學院，浙江 杭州311300；3.南京林業大學 生物與環境學院，江蘇 南京 210037；4.福建師范大學 地理科學學院，福建 福州350007）

隨著茶葉價格飄升，茶農積極性高漲，大面積的山地被開發成茶園。茶葉大面積栽培，促進當地經濟發展的同時引發了一系列生態環境問題。茶園面積過度擴張，種植面積明顯超過理想面積，導致森林資源破壞，水土流失嚴重［1］。茶園除草、施肥等經營措施，又導致土壤肥力下降、生物多樣性降低、土地退化等系列問題。明確茶園的空間分布和栽培面積，是政府部門管理決策的關鍵。但茶園分布廣泛，面積增加迅速，實地調查準確的茶園面積及空間信息需耗費大量人力物力。遙感技術可有效地獲取地表信息，既可以大面積、周期性的重復觀測，也可以節省大量的人力、物力和時間［2］。過去的幾十年中，有學者利用遙感技術進行經濟林專題信息提取等相關研究［3-4］。梁守真等［5］以MODIS和Landsat數據為基礎構建分類模型，開展橡膠Hevea brasiliensis林提取研究。XI等［6］基于多時相Landsat影像使用混合像元分解方法提取山核桃Carya cathayensis分布范圍，并檢測其受干旱干擾程度。用傳統Landsat等中等分辨率影像進行經濟林信息提取，受空間分辨率限制，精度不高。隨著遙感影像空間分辨率不斷提高，高分辨率影像越來越多被應用于林業遙感中［7］。WANG等［8］利用GF-1和ZY-3高分辨率影像，使用專家規則方法提取了香榧Torreya grandis‘Merrillii’在浙江的分布信息。梁文海等［9］基于面向對象方法使用GF-2影像提取廣西橫縣桉樹經濟林的分布。隨著影像空間分辨率的提高，地物的空間信息更加豐富，但 “同物異譜”現象也更加嚴重。由于高分辨率影像一般僅含紅、綠、藍、近紅外等4個波段，受光譜分辨率限制，在區分植被類型時有一定的局限性。2015年12月發射的Sentinel-2衛星包含3個紅邊波段，是綠色植物生長狀況的敏感性波段。如IMMITZER等［10］使用Sentinel-2數據對中歐農作物及樹木進行分類，表明紅邊及短波紅外波段可提取植被信息。茶園主要分布在景觀復雜的丘陵山區，且茶園與灌叢等木本植被的光譜特征相似［11］，這些因素為茶園遙感提取帶來一定困難，僅見個別學者應用遙感技術提取茶園信息［11-13］，尚無利用Sentinel-2數據提取茶園分布的研究。鑒于此，本研究以浙江西北部地區為研究區，基于Sentinel-2多光譜影像，結合茶園物候信息和經營管理模式，分析不同時間下茶園與其他地類的光譜特征，以紅邊波段構建歸一化茶園指數，使用決策樹方法對茶園信息進行提取，為Sentinel-2數據應用于植被提取及監測提供參考。

1 研究區概況

研究區位于浙江省西北部（30°39′～30°52′N， 119°30′～119°58′E）（圖 1）， 東西距離為 46.43 km， 南北距離為28.77 km，總面積為1335.79 km2。研究區地形以山地丘陵為主，海拔為0～700 m。該地屬于亞熱帶季風氣候區，光照充足、氣候溫和、雨水充沛、四季分明，適宜茶樹生長。該地年降水量為1512.4 mm，年日照時數為1961.5 h，年平均氣溫為16.9℃。研究區內的茶園類型主要為白茶Koilodepas hainanense樹，屬灌木型，中葉類，主干明顯，為安吉縣所特有［14］。根據安吉縣森林資源規劃設計調查成果報告，2016年安吉縣有白茶12058 hm2，占經濟林面積的60.13%，主要分布在梅溪、溪龍、遞鋪、天子湖、孝源、孝豐、杭垓鎮等鄉鎮（街道）。白茶在1998年不足100 hm2，2007年增加到5958 hm2，2016年增加到12058 hm2，增長迅速。

2 數據收集及預處理

2.1 數據收集

本研究使用的遙感數據是從歐洲太空局的Copernicus Open Access Hub網站上獲取Sentinel-2衛星數據，是大氣表觀反射率一級產品（L1C），共4景無云數據，時間分別為2017年10月31日，2018年2月13日，2018年5月4日和2018年7月18日。Sentinel-2衛星搭載多光譜影像儀（multi-spectral instrument,MSI），拍攝影像包含13個波段，其中藍、綠、紅和近紅外波段空間分辨率為10 m，3個紅邊波段、窄波近紅外波段、2個短波紅外波段空間分辨率為20 m，沿海氣溶膠波段、水蒸氣波段和卷云波段的空間分辨率為60 m。本研究選用分辨率為10和20 m共10個多光譜波段。數字高程模型（DEM）下載于美國航空航天局（NASA）官網，空間分辨率為12.5 m。

圖1 研究區位置及驗證樣點分布示意圖Figure 1 Location of study area and validation points

2016-2018年通過野外調查收集了研究區內主要綠色植被類型的地面數據，如茶園、毛竹Phyllostachys edulis林、闊葉林和針葉林等，其中的毛竹林被分為大年和小年分別調查。通過佳能D7000 GPSCAMERA，獲取了帶有坐標信息的地類照片300張。同時，通過野外調查和室內勾繪的方法在奧維互動地圖上標定了517個不規則多邊形地塊，將調研獲取的真實地塊信息轉化為矢量格式，以此作為本研究使用的樣本數據。

2.2 數據預處理

數據預處理主要包括投影系統轉換、大氣校正、幾何校正、地形校正4部分。首先將本研究的空間數據統一轉換為UTM 50N投影坐標系統，大地基準為WGS 1984。使用歐洲太空局的SNAP軟件中的Sen2Cor插件對原始數據進行大氣校正，獲取地表反射率數據［15］。基于雙3次卷積插值法，將校正后的10個波段重采樣到10 m。在Sentinel-2影像和DEM上選取均勻分布的30個控制點，采用3次卷積多項式模型進行幾何校正，均方根誤差控制在0.5個像元之內。

為了減少地形對分類的影響，基于DEM數據，采用分坡度的C校正模型［16］對Sentinel-2影像進行地形校正。C校正模型是基于影像像元值與太陽入射角余弦值構建線性關系，并加入經驗參數對影像進行校正。其計算公式如下：

式（1）～（2）中：ρm為校正后像元值；ρ為校正前像元值；c為校正系數；i為太陽入射角；β為太陽天頂角；λ為太陽方位角；θ為坡度；ω為坡向。

從影像頭文件中獲取太陽天頂角和高度角，基于DEM計算坡度和坡向。將坡度分為4個等級，分別對影像的像元值和太陽入射角余弦值進行擬合，獲取校正系數c。結果顯示：各波段的校正系數均小于1，校正后的影像凹凸立體感明顯消除并趨于扁平（圖2），陰坡和陽坡亮度值差異得到改善，同時對比校正前后影像的均值和標準差（表1），校正后各波段均值增加，標準差減小，地形校正效果較好。

3 研究方法

茶樹屬于灌木，四季常綠，與灌叢等植被的光譜特征相似。但研究區內茶樹有其特殊性，3-4月是茶葉采摘期，4月下旬和5月初進行修剪，受人為經營管理影響，茶園結構單一，林下無植被。這些特征為區分茶樹與其他植被提供了可能。本研究分別對不同地類的光譜特征、季節特征、茶園指數構建和茶園信息提取與精度驗證進行分析。

表1 地形校正前后影像反射率的均值和標準差Table 1 Mean and standard of image reflectivity before and after topographic

3.1 光譜特征分析

綜合考慮茶園的物候特征、人為干擾和遙感數據的可獲取性，本研究選擇了茶園生長減緩期（12月）、休眠期（2月）、修剪期（5月）和生長旺盛期（7月）4個時期的Sentinel-2數據的地表反射率圖像，分別進行地物光譜特征分析。基于野外調查的樣本數據，從Sentinel-2影像上選擇2017-2018年沒有變化的6類典型植被樣本。其中毛竹林選取了小年830個像元，大年960個像元，闊葉林940個像元，針葉林800個像元，茶園680個像元，農田780個像元，每個像元大小為10 m×10 m。從影像上提取每個地類在不同光譜波段上的光譜反射率，繪制不同季節不同地類的光譜曲線圖（圖3），對比分析茶園的物候特征，以及茶園與其他地類在光譜曲線上的差異。

圖2 地形校正前后對比示意圖Figure 2 Images comparison before and after terrain correction

在茶園生長減緩期（12月），茶園和其他植被類型均表現出健康植被的光譜特征，針葉林在各個波段上反射率均低于其他地類，茶園僅在第3個紅邊波段和2個近紅外波段的反射率比其他植被略高，但差異不明顯。在冬末春初的2月，針葉林和闊葉林在各個波段的反射率低于其他地類，而農田還未種植作物，呈現裸地光譜特征，即農田在紅邊和近紅外的反射率明顯降低，而在短波紅外波段上的反射率則明顯高于其他地類，而茶園和竹林在冬季常綠，光譜特征接近，這個時期的茶園與竹林易混淆。在春末夏初的5月，闊葉林、小年竹林開始換葉，它們在紅邊和近紅外的反射率明顯升高，而茶園由于統一修剪，枝葉減少，裸地露出，此時所表現出的光譜特征與冬季的農田類似，即茶園在紅邊和近紅外波段的反射率明顯下降。除此之外，茶園在紅邊3波段上的反射率低于其他所有植被地類，而在短波紅外2波段上的反射率則高于其他所有植被。在7月，各個植被類型都表現出生長旺盛期的植被光譜特征，各個地類之間的光譜反射率差異較小。

總體上，研究區內的植被在不同時期表現出明顯的季節特征和光譜差異，茶園和其他典型植被在12月、2月和7月表現的光譜特征較一致，區分性不大，而5月茶園與其他植被類型光譜曲線特征明顯不一致，因此可確定茶園提取的最佳時間在5月。由于農田存在休耕期，因此在5月茶園最容易混淆的地類為農田和裸地。

3.2 季節指數構建與對比

歸一化植被指數NDVI（normalized difference vegetation index）是反映植被生長狀態的指示因子［17］。本研究首先計算了各個月的歸一化植被指數。計算公式為：

圖3 研究區內Sentinel-2影像真彩色組合圖及主要植被類型光譜曲線示意圖Figure 3 True color combination of Sentinel-2 images and reflectance curves of vegetation

式（3）中：m表示月份，如NDVI_5表示5月的歸一化植被指數；ρnir和ρred分別為Sentinel-2數據的近紅外波段（中心波長842 nm）和紅波段（中心波長665 nm）的反射率。

根據多時相Sentinel-2數據的光譜特征曲線分析結果，在2和5月，茶園和農田在紅邊波段、近紅外以及短波紅外波段與其他地類的明顯差異。計算短波紅外2和紅邊波段3的歸一化指數時可以突出茶園和農田在2個波段上的差異，因此構建歸一化茶園指數NDTI（normalized difference tea garden index）。表達式如下：

式（4）中：m表示月份，如NDTI_5表示5月的歸一化茶園指數，ρswir2和ρred-edge3分別為Sentinel-2數據的短波紅外波段（中心波長2190 nm）和紅邊波段（中心波長783 nm）的反射率。根據式（3）和式（4），分別計算不同月份的歸一化植被指數和歸一化茶園指數，統計分析研究區內主要地類的歸一化植被指數和歸一化茶園指數值，對比分析茶園與其他地類的可區分性。如圖4所示：不透水地表和水體的歸一化植被指數和歸一化茶園指數表現出完全相反的指數值，2個指數都可以將不透水地表和水體與其他植被區分，其中NDVI_2對這2個地類的區分性更好，而且合適的取值還可以將與茶園最易混淆的農田信息剔除。茶園與其他典型的植被在10，2和7月的歸一化茶園指數值比較接近，區分度很低，這是由于它們光譜特征接近，而在5月，茶園歸一化茶園指數值明顯大于其他植被，因此NDTI_5是區分茶園與其他植被的最優指數。

3.3 茶園信息提取與精度驗證

本研究基于光譜特征分析和建立的歸一化茶園指數，計算了4個不同季節的歸一化植被指數和歸一化茶園指數，對比了研究區內主要地類在不同季節指數上的特征，最大程度地提高了茶園與其他地類的可分離度。最終基于歸一化植被指數NDVI_2和歸一化茶園指數NDTI_5構建決策樹，進行茶園信息的提取。

本研究精度驗證的方法是混淆矩陣法。分類數據是基于Sentinel-2影像的決策樹分類結果，參考數據是空間分辨率為0.6 m的2018年谷歌地球數據和后續實地調查的真實地塊數據。分類結果為2類：一類是茶園；另一類是其他地類，主要有農田、不透水地表、植被和裸土等。為保證精度驗證點的客觀性，從分類結果中分層隨機選取600個驗證點，即茶園和其他地類各300個驗證點，結合谷歌地球影像數據和實地調查數據進行精度評價，利用混淆矩陣和Kappa系數方法，對茶園提取結果進行精度評價。

圖4 研究區內主要地表類型的歸一化植被指數（NDVI）與歸一化茶園指數（NDTI）Figure 4 Normalized difference vegetation index（NDVI） and normalized tea garden index（NDTI） of typical land cover

4 結果與分析

4.1 各地類在植被指數上的特征

從圖5可以看出：各個地類在不同季節的指數上反映的特征不同。在2月歸一化植被指數圖（圖5A）上，水體呈現黑色，裸土、農田和不透水地表呈灰黑色，而茶園和植被呈亮白色。在2月歸一化茶園指數圖（圖5B）上，茶園、部分植被和水體的值最低，呈黑色，裸土、農田、部分植被呈淺灰色，部分不透水地表呈白色。在5月歸一化植被指數圖（圖5C）上，水體值最低，呈黑色，農田、茶園、裸地和不透水地表呈灰色，而植被呈亮白色。在5月歸一化茶園指數圖（圖5D）中，水體和部分不透水地表的值最高，呈亮白色，茶園、部分不透水地表、裸地的值較低，呈灰色，而農田和植被類型呈黑色。

圖5 各地類在指數上的特征示意圖Figure 5 Characteristic graphs of main land cover in index

4.2 茶園信息提取的指數選擇

結合研究區植被類型和其他地類在2個歸一化指數上的表現，茶園信息提取的關鍵在與找到它與其他植被類型的區別，因此本研究選取了冬季的歸一化植被指數和5月的歸一化茶園指數2個最優變量進行茶園信息的提取。首先，選用了冬季的歸一化植被指數NDVI_2區分植被和非植被地類，在決策樹中選取閾值0.3，將植被和非植被區分開來，即當NDVI_2＞0.3的像元為植被地類，而NDVI_2＜0.3的像元則為其他地類，如裸土、不透水地表、水體等。結合圖4和圖5可知：5月歸一化植被指數和茶園指數區分茶園與其他植被的效果較好，且5月歸一化茶園指數優于5月歸一化植被指數。因此，本研究選擇基于歸一化茶園指數對茶園與其他植被類型進行區分，選取閾值0.4，當NDTI_5＞0.4的像元皆為茶園信息，最終將決策樹分類結果匯總為茶園和非茶園2個大類，并獲取研究區的茶園空間分布數據。

4.3 茶園空間分布

圖6可知：研究區內茶園主要分布在研究區的中下部，總體呈西南—東北方向，片狀分布。結合研究區的地形圖可以發現：茶園主要分布在山地與平原交界處的緩坡丘陵地帶，坡度主要集中在5°～25°，這與茶園的生長環境要求相符。利用空間自相關分析發現：茶園的空間范圍與人類居住范圍空間相關性較強，這是由于大多數茶園需要人工養護，適當的距離可以節約成本。研究區內茶園總面積99.9 km2，約占研究區總面積的7.7%，按行政區劃統計，和平鎮的茶園分布最廣，達31.32 km2，其次為遞鋪鎮和梅溪鎮，分別為18.48和16.34 km2。按照茶園分布面積占行政區劃面積的比例來看，溪龍鄉最高，茶園比例接近33.0%，和平鎮次之，接近25.0%，梅溪鎮和皈山鄉的茶園比例在10.0%左右，遞鋪鎮和天子湖鎮則僅為6.0%。

分別隨機選取了研究區內茶園和其他地類各300個驗證點，在谷歌地球中確認其位置的正確性，對茶園提取結果進行精度驗證。經統計，在茶園的300個驗證點中，276個為茶園，24個為其他地類；而在其他地類的300個驗證點中，287個為其他地類，13個為茶園地類。基于NDTI_5提取的茶園生產者精度為95.50%，用戶精度為92.00%；其他地類的生產者精度為92.28%，用戶精度為95.67%。總體分類精度達93.83%，Kappa系數為0.917。

圖6 基于歸一化茶園指數的茶園空間分布示意圖Figure 6 Spatial distribution of tea garden in study area by normalized tea garden index

5 結論與討論

本研究以Sentinel-2數據為基礎，分析探討了基于物候特征構建指數提取茶園的方法。由于茶園常綠，容易和其他綠色植被混淆，因此本研究根據茶園生長及經營管理特征，分析了茶園和主要綠色植被在不同季節的光譜特征與差異，確定了提取茶園的最佳時間在5月，分別計算了歸一化植被指數和構建了歸一化茶園指數，對比分析了這2個指數在不同季節上茶園與其他植被的可分離性，確定了茶園指數提取的2個最優指數組合為冬季的歸一化植被指數和5月的歸一化茶園指數，利用決策樹分類實現了研究區內茶園信息的準確提取。結果表明：基于Sentinle-2紅邊波段構建的歸一化茶園指數能夠將茶園與易混淆的綠色植被區分開，分類總精度達93.83%，Kappa系數為0.917。

有研究表明：基于高分辨率數據使用面向對象、最大似然、支持向量機、隨機森林4種方法提取茶園，其中面向對象分類方法精度最高，為88.07%［11］。高分辨率影像分割的最佳尺度難確定，獲取對象高維特征后再優化，分類過程較繁瑣。研究表明：結合非監督分類和卷積神經網絡方法提取茶園，精度在88%以上［13］。但卷積神經網絡參數較復雜，需要編寫算法調整參數，解釋性較差。本研究根據植被物候特征，基于Sentinel-2紅邊波段構建茶園指數來提取茶園，易于理解和實現，為其他地區茶園提供借鑒，具有較高的推廣價值。

據安吉縣二類調查統計報告顯示：茶園在2006-2016年面積增加了1倍，大量森林被開發成茶園。由于茶園內植被覆蓋低，水土流失風險高，森林轉化成茶園，不僅增大了水土流失的風險，還間接造成了森林碳儲量的減少。由Sentinel-2紅邊波段構建的茶園指數能夠很好地區分出茶園與其他地類，但是Sentinel-2缺乏歷史數據，研究近幾十年的茶園動態變化，Landsat數據是主要選擇，因此如何用Landsat數據更準確地提取茶園、監測茶園變化是下一步研究重點。本研究表明：5月的歸一化植被指數也可以區分茶園與其他地類，但茶園和休耕農田存在一定程度的混淆，由于農田主要分布在平原地區，坡度較緩，而茶園主要分布在坡度小于25°以下的地區，后續可以加入坡度數據來剔除農田，進一步提高茶園的提取精度。