潘永地， 黃鵬， 梁金晨

(1.溫州市氣象局 溫州臺風監(jiān)測預報技術重點實驗室，浙江 溫州 325027； 2.吉林省衛(wèi)星遙感應用技術重點實驗室，吉林 長春 130102；3.瑞安市氣象局，浙江 瑞安 325200)

茶園分布對于政府部門開展茶葉生產(chǎn)技術指導、防災減災以及掌握茶葉整體生產(chǎn)形勢具有極為重要的意義。在傳統(tǒng)上通過各地上報可以得出粗略的各行政區(qū)域茶園面積，但是分布位置模糊，無地塊形狀等信息，更新困難，在使用上存在很多不足。隨著遙感技術的發(fā)展，快速客觀提取茶園分布成為可能，國內外許多學者開展了相關研究。Fauziana等[1]運用歸一化差值植被指數(shù)(NDVI)估計茶冠密度，進而估算茶葉產(chǎn)量。Akar等[2]將NDVI、灰度共生矩陣(GLCM)和帶通(Gabor)濾波等紋理提取方法與隨機森林(RF)算法相結合，提高了分類器分類的精度。朱澤潤[3]以光譜特征和Gabor紋理特征作為底層特征，從中提取場景的視覺詞袋特征，并使用支持向量機完成茶園提取。王斌等[4]利用RF對茶園影像特征進行重要性評估、排序和特征選擇，設計單季節(jié)初始特征集、單季節(jié)優(yōu)選特征集、多季節(jié)優(yōu)選特征集進行茶園提取，表明多季節(jié)優(yōu)選特征集具有最好的性能表現(xiàn)。許光明[5]將面向對象和多源數(shù)據(jù)融合的多層次規(guī)則分類方法應用于茶園遙感提取研究。馬超等[6]使用面向對象方法和決策樹分類模型對Landsat影像提取初步分類結果，再利用MODIS增強植被指數(shù)時序提取茶園。楊艷魁[7]利用GLCM、Gabor濾波器以及局部二值模式(LBP)，分別構建了適合于World view-2影像和GF-2影像的紋理特征，對比5種不同光譜、紋理特征組合方案，結果表明，結合光譜信息和LBP_GABOR紋理的分類方法提取茶園分布的效果最好。吳錦玉[8]基于空間點模式分析增強茶園紋理，并基于紋理模式以面向對象非參數(shù)化方法實現(xiàn)茶園識別。黃艷紅[9]分析了地形校正和不同分類器對茶園面積提取的影響，比較了Landsat8 OLI、GF-1 WFV和Sentinel-2 MSI等3種影像提取茶園的效果，結果表明，Sentinel-2 MSI數(shù)據(jù)和紅邊波段可以有效提高茶園遙感分類的精度，利用最小二乘法將HJ-1 CCD數(shù)據(jù)、GF-1 WFV數(shù)據(jù)和Sentinel-2 MSI數(shù)據(jù)的分類特征相對校正到Landsat8 OLI數(shù)據(jù)集，并構建多時相的Landsat8 OLI數(shù)據(jù)分類特征數(shù)據(jù)集，再使用RF分類器分類提高提取精度。Rajapakse等[10]對通過分光光度計測得的NDVI與茶園的LAI做線性回歸，發(fā)現(xiàn)其有較強的相關性。趙曉晴等[11]利用Sentinel-2A影像分析了7種典型地物的時序光譜變化與NDVI變化，結果表明，分類精度較高的是SR-SWIR2-NIR_May。陳慧等[12]以GF-2影像為數(shù)據(jù)源，確定茶園最優(yōu)分割尺度為170，并利用面向對象分類方法實現(xiàn)茶園提取。蒙良莉等[13]基于Sentinel遙感數(shù)據(jù)，通過多特征耦合優(yōu)化模式的分類，實現(xiàn)對紅樹林信息提取。王亞琴等[14]基于MODIS、Landsat TM、Quick Bird數(shù)據(jù)，對山西蘆芽山2000、2005和2010年植被變化作了定量和定性分析。林娜等[15]基于GF-1號影像提取了南方水稻種植信息。

上述研究主要利用茶園或其他植被的紋理和光譜特征進行提取，較好地克服了傳統(tǒng)人工調查統(tǒng)計中耗時長、位置不精確等方面的不足，但是因為品種、季節(jié)、栽培管理方式的不同，茶園的光譜特征和紋理特征都可能不同，紋理識別需要大量圖片作為訓練樣本，所以研究一套普適性的茶園提取方法還存在一定困難。本文通過分析浙江省永嘉縣茶葉的管理方式，總結永嘉縣茶園與其他類似度或出現(xiàn)頻率較高植被的光譜特征，利用茶園和其他植被在光譜特征上的差異及管理方式造成光譜時序變化差異進行永嘉縣茶園的提取，結果表明，利用管理方式造成的NDVI時序變化特征結合茶園固有光譜特征提取出的茶園在面積、位置、形狀上都較接近于實際情況。

1 材料與方法

1.1 試驗流程

我們通過實地考查確定出茶園及易干擾茶園提取的典型植被代表地塊作為分析點，測定經(jīng)、緯度。下載歷史遙感數(shù)據(jù)，剔除分析點受薄云、碎云、卷云和霧等影響的資料，讀取各分析點上各波段數(shù)值。計算多種指數(shù)，形成單波段值及各指數(shù)的時間變化序列，確定出茶園和其他點差異較大的光譜特征，根據(jù)這些特征制定茶園提取方案。對提取出的茶園分布總面積與行政統(tǒng)計面積進行比較，選擇重點區(qū)域分析提取出的茶園位置和形狀與經(jīng)人工判讀結合實地調查修正的同區(qū)域實際茶園分布進行比較，分別考察提取茶園分布的誤差，證明提取方案的可行性。

1.2 分析的指數(shù)

歸一化差值植被指數(shù)NDVI：

NDVI=(Nir-Red)/(Nir+Red)=(B7-B3)/(B7+B3)。

式中，Nir、Red分別為近紅外波段、紅波段的反射率，在本研究中對應波段編號B7、B3。

紅綠比率指數(shù)RGRI：

RGRI=Green/Red=B3/B2。

式中，Red、Green分別表示紅光、綠光波段反射率，在本研究中的波段編號依次為B3和B2。

近紅外綠比率指數(shù)NGRI：

NGRI=Nir/Green=B7/B2。

該指數(shù)系本文參照紅綠比指數(shù)定義。

1.3 分析點地理位置

本研究選取茶園5處、灌木叢1處、毛竹林1處、松樹林1處、常綠闊葉林1處，其中常綠闊葉林主要是樟樹。各分析點具體經(jīng)、緯度見表1。

表1 各分析點的地理位置

1.4 遙感影像數(shù)據(jù)

遙感影像數(shù)據(jù)采用Sentinel-2 L2A(L2A級產(chǎn)品是經(jīng)過大氣校正的BOA反射率數(shù)據(jù))2019、2020年影像數(shù)據(jù)。Sentinel-2為高分辨率多光譜成像衛(wèi)星，是歐洲空間局(European space agency，ESA)的哥白尼計劃系列衛(wèi)星的重要組成部分，包括Sentinel-2A和Sentinel-2B衛(wèi)星。單星重訪周期為10 d，雙星重訪周期為5 d。在兩星覆蓋條帶的重疊部分，重訪周期可達2～3 d。其主要有效載荷為MSI(多光譜成像儀)，共有13個波段，光譜范圍在400～2 400 nm，空間分辨率可見光10 m，近紅外20 m，短波紅外60 m，成像幅寬290 km。本研究中下載的波段及命名如表2所示。

表2 Sentinel-2波段與下載影像命名

獲取的Sentinel-2遙感影像產(chǎn)品日期為：2019年的1月31日、3月12日、3月15日、3月30日、4月1日、4月24日、5月1日、5月24日、6月5日、6月15日、6月30日、7月25日、8月2日、8月12日、8月27日、9月8日、9月11日、9月16日、9月18日、9月23日、9月26日、10月3日、10月11日、10月16日、10月18日、10月21日、10月23日、10月31日、11月5日、11月10日、11月15日、12月2日、12月7日、12月10日、12月12日、12月17日、12月27日；2020年的1月31日、3月12日、3月15日、3月30日、4月1日、4月24日、5月1日、5月24日、6月5日、6月15日、6月30日、7月25日、8月2日、8月12日、8月14日、8月17日、8月27日、9月8日、9月11日、9月16日、9月18日、9月23日、9月26日、10月3日、10月11日、10月16日、10月18日、10月21日、10月23日、10月31日、11月5日、11月10日、11月15日、12月2日、12月7日、12月10日、12月12日、12月17日、12月27日。

用于局部對照的是分辨率為0.75 m的吉林一號寬幅01a星2021年1月9日高分影像。

2 結果與分析

2.1 各波段反射率變化特征

圖1、圖2分別是各分析點在Red、Nir波段的反射率變化，其中茶園有多個分析點數(shù)據(jù)，本研究采用平均值。Red在4—6月存在一個最大值，在8—9月存在一個最小值。經(jīng)調查了解，永嘉茶葉普遍在4月份開展修剪，晚一點的可能在5月份修剪，修剪基本采用重剪(將茶葉的冠層全部砍掉)，修剪后的枝條均勻鋪在株間。剛修剪后鋪在株間的枝條與修剪前在Red波段反射率上無大的差別，之后逐漸干枯，Red反射率逐漸達到最大值。經(jīng)過一段時間，修剪后的茶樹植株逐漸產(chǎn)生嫩梢并逐漸生長擴大葉面積，在7月份基本恢復至修剪前水平，9月份葉片生長和綠度基本達到最大，Red波段反射率逐漸達到最低，之后隨著氣溫的逐漸降低，葉片綠度略有降低，Red波段反射率略有上升。其他植被分析點由于沒有茶園的人為修剪，所以未表現(xiàn)出4—6月的明顯峰值，在9月基本上也是其他常綠植物最綠的時期，所以也與茶園具有基本一致的谷值出現(xiàn)。茶園Nir在變化趨勢上與Red相反，在4—6月有一個谷值，在8—9月有一個峰值。4—6月以外的時間內各植被在Nir反射率變化上未表現(xiàn)出明顯的差異，Nir在8月份后穩(wěn)定地表現(xiàn)為茶園最大。此外，通過對Re2、Swir1波段的分析，Re2與Nir變化趨勢基本一致，Swir1與Red變化趨勢基本一致，空間分辨率上Red、Nir較Swir1、Re2高一些，故選取Red、Nir做差異性研究。Blue、Green、Gb等其他波段反射率均沒有明顯體現(xiàn)出茶葉與其他植物的差異。

圖1 茶園及其他植被分析點Red波段反射率時序變化

圖2 茶園及其他植被分析點Nir波段反射率時序變化

根據(jù)Red、Nir的時間變化特點和NDVI定義，使用NDVI必然較單獨使用Red、Nir單波段值更能強化茶園與其他植被之間的差異。但是，該特征包含一些干擾信息，主要是一些山區(qū)在4—6月可能會出現(xiàn)森林火災或開墾，使NDVI達到最小值，在9—10月后逐漸恢復植被，NDVI又達到最大值。Nir波段在8月之后具有茶園分析點一直高于其他植被分析點的特征，故可采用NIR波段來進一步提高茶園分布提取的精度，同時為了減少波動，與Green構成NGRI，用以區(qū)分茶園和其他人為因素導致在4—6月NDVI降低的植被。

2.2 NDVI和NGRI的時間變化特征

我們分別對5個茶園分析點的各指數(shù)值繪制時間變化曲線，5個點基本一致，故以5個點的平均值進行探討。圖3、4是NDVI、NGRI在茶園、常綠闊葉林、毛竹林、松樹林、灌木叢等分析點的時間變化。

圖3 茶園及其他植被NDVI時序變化

由圖3、4可見，NDVI在茶園具有與Nir類似的變化特征，但變化幅度更大；NGRI表現(xiàn)為9月份后常綠闊葉林、茶園較其他植被大的特點。根據(jù)NDVI和NGRI的變化特征制定出茶園提取方案：利用2020年9月NDVI剔除城市用地、水體等區(qū)域(本文采用NDVI閾值為0.2)；計算2020年1—3月最大NDVI和4—6月最小NDVI，利用兩者的差值(本文采用NDVI差閾值為0.13)進一步提取出茶園的基礎區(qū)域；利用9—10月最大NGRI進一步去除近紅外綠比較小區(qū)域(本文采用6.0)。

2.3 山區(qū)茶園分布的提取結果及誤差分析

我們根據(jù)2.2節(jié)制定方案提取出永嘉縣茶園的分布，并統(tǒng)計出面積，提取結果見圖5。

圖4 茶園及其他植被NGRI時序變化

圖5 永嘉茶園提取結果

誤差分析從兩方面考察：提取的茶園面積與統(tǒng)計年鑒上的茶業(yè)面積相比較；通過局部的提取茶園與同區(qū)域通過高分影像人工判讀得出的茶園在位置、形狀上的吻合情況來分析。根據(jù)溫州市統(tǒng)計年鑒，永嘉縣2020年茶園面積為3 315 hm2，本文提取茶園平面(正射投影)面積為2 997.38 hm2，再利用ALOS衛(wèi)星生產(chǎn)的12.5 m分辯率DEM折算到坡度上得到斜面面積為3 210.23 hm2，面積擬合率達96.8%，誤差為3.2%。為了更好地判斷提取茶園與實際茶園在位置、形狀上的吻合情況，我們對永嘉縣東南部茶葉分布最多的三江、烏牛街道進行提取區(qū)域和高分影像人工判讀出的區(qū)域進行對比(圖6、7)。從圖6、7可知，根據(jù)2.2中方案提取出的茶園分布圖與高分影像人工判讀出的茶園分布圖在位置分布上非常一致，在形狀上也比較接近。

圖6 自動提取(左)與人工判讀(右)空間分布對比

3 小結與討論

永嘉茶葉以烏牛早為主，近兩年白茶有所增加，在4—5月進行重剪，修剪后的枝條鋪于原植株行間，該管理方式導致茶園NDVI在4—6月期間逐漸降低達到最小值后又逐漸升高，在9月達到最大值，這一特征顯著區(qū)別于其他植被。茶園和常綠闊葉林的NGRI在9—10月較其他植被偏高，可以用于剔除一些如森林火災后恢復植被的情況，從而提高茶園提取的精度。根據(jù)茶葉修剪造成的茶園NDVI時序變化顯著特征進行NDVI變化幅度閾值茶園提取，并以NGRI閾值剔除一些具有類似的NDVI變化特征的非茶園區(qū)域，最終可以得到面積精度高于90%，位置、形狀上也與實際情況較為吻合的茶園分布圖。

本方法具有影像資源獲取方便、提取快速的優(yōu)勢，較人工判讀大大提高效率，能夠完成人工判讀較難完成的大范圍判讀。本方法主要根據(jù)光譜特征提取茶園，仍會受到山體、云的陰影等因素影響，開展各種陰影條件下的糾正值得進一步研究。

圖7 自動提取(黃線)與人工判讀(綠線)形狀對比