基于光谱特征提取永嘉县茶园分布

潘永地， 黄鹏， 梁金晨

(1.温州市气象局 温州台风监测预报技术重点实验室，浙江 温州 325027； 2.吉林省卫星遥感应用技术重点实验室，吉林 长春 130102；3.瑞安市气象局，浙江 瑞安 325200)

茶园分布对于政府部门开展茶叶生产技术指导、防灾减灾以及掌握茶叶整体生产形势具有极为重要的意义。在传统上通过各地上报可以得出粗略的各行政区域茶园面积，但是分布位置模糊，无地块形状等信息，更新困难，在使用上存在很多不足。随着遥感技术的发展，快速客观提取茶园分布成为可能，国内外许多学者开展了相关研究。Fauziana等[1]运用归一化差值植被指数(NDVI)估计茶冠密度，进而估算茶叶产量。Akar等[2]将NDVI、灰度共生矩阵(GLCM)和带通(Gabor)滤波等纹理提取方法与随机森林(RF)算法相结合，提高了分类器分类的精度。朱泽润[3]以光谱特征和Gabor纹理特征作为底层特征，从中提取场景的视觉词袋特征，并使用支持向量机完成茶园提取。王斌等[4]利用RF对茶园影像特征进行重要性评估、排序和特征选择，设计单季节初始特征集、单季节优选特征集、多季节优选特征集进行茶园提取，表明多季节优选特征集具有最好的性能表现。许光明[5]将面向对象和多源数据融合的多层次规则分类方法应用于茶园遥感提取研究。马超等[6]使用面向对象方法和决策树分类模型对Landsat影像提取初步分类结果，再利用MODIS增强植被指数时序提取茶园。杨艳魁[7]利用GLCM、Gabor滤波器以及局部二值模式(LBP)，分别构建了适合于World view-2影像和GF-2影像的纹理特征，对比5种不同光谱、纹理特征组合方案，结果表明，结合光谱信息和LBP_GABOR纹理的分类方法提取茶园分布的效果最好。吴锦玉[8]基于空间点模式分析增强茶园纹理，并基于纹理模式以面向对象非参数化方法实现茶园识别。黄艳红[9]分析了地形校正和不同分类器对茶园面积提取的影响，比较了Landsat8 OLI、GF-1 WFV和Sentinel-2 MSI等3种影像提取茶园的效果，结果表明，Sentinel-2 MSI数据和红边波段可以有效提高茶园遥感分类的精度，利用最小二乘法将HJ-1 CCD数据、GF-1 WFV数据和Sentinel-2 MSI数据的分类特征相对校正到Landsat8 OLI数据集，并构建多时相的Landsat8 OLI数据分类特征数据集，再使用RF分类器分类提高提取精度。Rajapakse等[10]对通过分光光度计测得的NDVI与茶园的LAI做线性回归，发现其有较强的相关性。赵晓晴等[11]利用Sentinel-2A影像分析了7种典型地物的时序光谱变化与NDVI变化，结果表明，分类精度较高的是SR-SWIR2-NIR_May。陈慧等[12]以GF-2影像为数据源，确定茶园最优分割尺度为170，并利用面向对象分类方法实现茶园提取。蒙良莉等[13]基于Sentinel遥感数据，通过多特征耦合优化模式的分类，实现对红树林信息提取。王亚琴等[14]基于MODIS、Landsat TM、Quick Bird数据，对山西芦芽山2000、2005和2010年植被变化作了定量和定性分析。林娜等[15]基于GF-1号影像提取了南方水稻种植信息。

上述研究主要利用茶园或其他植被的纹理和光谱特征进行提取，较好地克服了传统人工调查统计中耗时长、位置不精确等方面的不足，但是因为品种、季节、栽培管理方式的不同，茶园的光谱特征和纹理特征都可能不同，纹理识别需要大量图片作为训练样本，所以研究一套普适性的茶园提取方法还存在一定困难。本文通过分析浙江省永嘉县茶叶的管理方式，总结永嘉县茶园与其他类似度或出现频率较高植被的光谱特征，利用茶园和其他植被在光谱特征上的差异及管理方式造成光谱时序变化差异进行永嘉县茶园的提取，结果表明，利用管理方式造成的NDVI时序变化特征结合茶园固有光谱特征提取出的茶园在面积、位置、形状上都较接近于实际情况。

1 材料与方法

1.1 试验流程

我们通过实地考查确定出茶园及易干扰茶园提取的典型植被代表地块作为分析点，测定经、纬度。下载历史遥感数据，剔除分析点受薄云、碎云、卷云和雾等影响的资料，读取各分析点上各波段数值。计算多种指数，形成单波段值及各指数的时间变化序列，确定出茶园和其他点差异较大的光谱特征，根据这些特征制定茶园提取方案。对提取出的茶园分布总面积与行政统计面积进行比较，选择重点区域分析提取出的茶园位置和形状与经人工判读结合实地调查修正的同区域实际茶园分布进行比较，分别考察提取茶园分布的误差，证明提取方案的可行性。

1.2 分析的指数

归一化差值植被指数NDVI：

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

式中，Nir、Red分别为近红外波段、红波段的反射率，在本研究中对应波段编号B7、B3。

红绿比率指数RGRI：

RGRI=Green/Red=B3/B2。

式中，Red、Green分别表示红光、绿光波段反射率，在本研究中的波段编号依次为B3和B2。

近红外绿比率指数NGRI：

NGRI=Nir/Green=B7/B2。

该指数系本文参照红绿比指数定义。

1.3 分析点地理位置

本研究选取茶园5处、灌木丛1处、毛竹林1处、松树林1处、常绿阔叶林1处，其中常绿阔叶林主要是樟树。各分析点具体经、纬度见表1。

表1 各分析点的地理位置

1.4 遥感影像数据

遥感影像数据采用Sentinel-2 L2A(L2A级产品是经过大气校正的BOA反射率数据)2019、2020年影像数据。Sentinel-2为高分辨率多光谱成像卫星，是欧洲空间局(European space agency，ESA)的哥白尼计划系列卫星的重要组成部分，包括Sentinel-2A和Sentinel-2B卫星。单星重访周期为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遥感影像产品日期为：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波段的反射率变化，其中茶园有多个分析点数据，本研究采用平均值。Red在4—6月存在一个最大值，在8—9月存在一个最小值。经调查了解，永嘉茶叶普遍在4月份开展修剪，晚一点的可能在5月份修剪，修剪基本采用重剪(将茶叶的冠层全部砍掉)，修剪后的枝条均匀铺在株间。刚修剪后铺在株间的枝条与修剪前在Red波段反射率上无大的差别，之后逐渐干枯，Red反射率逐渐达到最大值。经过一段时间，修剪后的茶树植株逐渐产生嫩梢并逐渐生长扩大叶面积，在7月份基本恢复至修剪前水平，9月份叶片生长和绿度基本达到最大，Red波段反射率逐渐达到最低，之后随着气温的逐渐降低，叶片绿度略有降低，Red波段反射率略有上升。其他植被分析点由于没有茶园的人为修剪，所以未表现出4—6月的明显峰值，在9月基本上也是其他常绿植物最绿的时期，所以也与茶园具有基本一致的谷值出现。茶园Nir在变化趋势上与Red相反，在4—6月有一个谷值，在8—9月有一个峰值。4—6月以外的时间内各植被在Nir反射率变化上未表现出明显的差异，Nir在8月份后稳定地表现为茶园最大。此外，通过对Re2、Swir1波段的分析，Re2与Nir变化趋势基本一致，Swir1与Red变化趋势基本一致，空间分辨率上Red、Nir较Swir1、Re2高一些，故选取Red、Nir做差异性研究。Blue、Green、Gb等其他波段反射率均没有明显体现出茶叶与其他植物的差异。

图1 茶园及其他植被分析点Red波段反射率时序变化

图2 茶园及其他植被分析点Nir波段反射率时序变化

根据Red、Nir的时间变化特点和NDVI定义，使用NDVI必然较单独使用Red、Nir单波段值更能强化茶园与其他植被之间的差异。但是，该特征包含一些干扰信息，主要是一些山区在4—6月可能会出现森林火灾或开垦，使NDVI达到最小值，在9—10月后逐渐恢复植被，NDVI又达到最大值。Nir波段在8月之后具有茶园分析点一直高于其他植被分析点的特征，故可采用NIR波段来进一步提高茶园分布提取的精度，同时为了减少波动，与Green构成NGRI，用以区分茶园和其他人为因素导致在4—6月NDVI降低的植被。

2.2 NDVI和NGRI的时间变化特征

我们分别对5个茶园分析点的各指数值绘制时间变化曲线，5个点基本一致，故以5个点的平均值进行探讨。图3、4是NDVI、NGRI在茶园、常绿阔叶林、毛竹林、松树林、灌木丛等分析点的时间变化。

图3 茶园及其他植被NDVI时序变化

由图3、4可见，NDVI在茶园具有与Nir类似的变化特征，但变化幅度更大；NGRI表现为9月份后常绿阔叶林、茶园较其他植被大的特点。根据NDVI和NGRI的变化特征制定出茶园提取方案：利用2020年9月NDVI剔除城市用地、水体等区域(本文采用NDVI阈值为0.2)；计算2020年1—3月最大NDVI和4—6月最小NDVI，利用两者的差值(本文采用NDVI差阈值为0.13)进一步提取出茶园的基础区域；利用9—10月最大NGRI进一步去除近红外绿比较小区域(本文采用6.0)。

2.3 山区茶园分布的提取结果及误差分析

我们根据2.2节制定方案提取出永嘉县茶园的分布，并统计出面积，提取结果见图5。

图4 茶园及其他植被NGRI时序变化

图5 永嘉茶园提取结果

误差分析从两方面考察：提取的茶园面积与统计年鉴上的茶业面积相比较；通过局部的提取茶园与同区域通过高分影像人工判读得出的茶园在位置、形状上的吻合情况来分析。根据温州市统计年鉴，永嘉县2020年茶园面积为3 315 hm2，本文提取茶园平面(正射投影)面积为2 997.38 hm2，再利用ALOS卫星生产的12.5 m分辩率DEM折算到坡度上得到斜面面积为3 210.23 hm2，面积拟合率达96.8%，误差为3.2%。为了更好地判断提取茶园与实际茶园在位置、形状上的吻合情况，我们对永嘉县东南部茶叶分布最多的三江、乌牛街道进行提取区域和高分影像人工判读出的区域进行对比(图6、7)。从图6、7可知，根据2.2中方案提取出的茶园分布图与高分影像人工判读出的茶园分布图在位置分布上非常一致，在形状上也比较接近。

图6 自动提取(左)与人工判读(右)空间分布对比

3 小结与讨论

永嘉茶叶以乌牛早为主，近两年白茶有所增加，在4—5月进行重剪，修剪后的枝条铺于原植株行间，该管理方式导致茶园NDVI在4—6月期间逐渐降低达到最小值后又逐渐升高，在9月达到最大值，这一特征显著区别于其他植被。茶园和常绿阔叶林的NGRI在9—10月较其他植被偏高，可以用于剔除一些如森林火灾后恢复植被的情况，从而提高茶园提取的精度。根据茶叶修剪造成的茶园NDVI时序变化显著特征进行NDVI变化幅度阈值茶园提取，并以NGRI阈值剔除一些具有类似的NDVI变化特征的非茶园区域，最终可以得到面积精度高于90%，位置、形状上也与实际情况较为吻合的茶园分布图。

本方法具有影像资源获取方便、提取快速的优势，较人工判读大大提高效率，能够完成人工判读较难完成的大范围判读。本方法主要根据光谱特征提取茶园，仍会受到山体、云的阴影等因素影响，开展各种阴影条件下的纠正值得进一步研究。

图7 自动提取(黄线)与人工判读(绿线)形状对比