多时相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°以下的地区，后续可以加入坡度数据来剔除农田，进一步提高茶园的提取精度。