基于Sentinel-2影像的厦门市茶园遥感提取

李 艳，张 帆

(河口生态安全与环境健康福建省高校重点实验室 / 厦门大学嘉庚学院环境科学与工程学院，福建 漳州 363105)

茶Camelliasinensis是常绿多年生木本植物，也是一种重要的经济作物，主要分布在我国热带和亚热带山区。福建省是我国茶叶主要产区之一，近年来，随着茶叶经济价值的不断提升和低碳农业发展理念的不断深入，茶园种植面积发生了明显的变化，茶园碳汇也成为农林业生态环境的重要组成部分[1]。因此，对茶园范围的监测具有重要的生态环境保护意义。

由于茶园分布广泛，且多分布于山区和丘陵地带，通过传统的实地调查获取茶树空间分布需花费大量时间和人力，具有较大的局限性。遥感技术具有宏观性、动态性、速度快、成本低等特点，是实现大面积、大尺度茶园种植区提取的重要手段[2—3]。光学遥感影像具有较高的空间分辨率和时间分辨率以及合适的光谱分辨率等优势，在植被遥感分类中被广泛使用[4—5]。马超等[6]基于Landsat8 和时序的MODIS 影像的光谱特征、植被指数和地形特征，采用决策树分类模型对研究区的茶园进行提取，总体分类精度达85.71%。赵晓晴等[7]利用时序Sentinel-2A 影像分析研究区典型地物的时序光谱特征与归一化植被指数(Normalized Difference Vegetation Index,NDVI)的时序特征，获得可用于茶园提取的特征波段，并通过这些特征波段结合NDVI 作为特征参数构建决策树实现研究区茶园的提取，但对于种植分散的茶园，该方法的普适性有待考证。李龙伟等[8]根据多时相Sentinel-2 数据的红边与短波红外波段构建了归一化茶园指数(Normalized Difference Tea Plantations Index,NDTI)，并基于NDTI 通过决策树模型提取了浙江西北部的茶园，提取总精度达93.83%。

由于茶树的光谱反射率值往往与其他一些植物如灌木、果园等重叠，引起提取误差，而使用基于空间信息的纹理特征可有效解决光谱混合问题[9]。杨艳魁等[10]利用高分二号(GF-2)影像数据构建差异归一化差分植被指数(Modified Normalized Difference Vegetation Index,DNDVI)和灰度共生纹理(Gray-level Co-occurrence Matrix,GLCT)，然后结合光谱和纹理特征运用支持向量机(Support Vector Machine,SVM)的方法对福建省安溪县茶园进行分类。结果表明，结合光谱特征和纹理特征提取的茶园精度(89.8%)高于仅用光谱特征提取的茶园精度(79.6%)。柏佳等[11]利用高分1 号(GF-1)的光谱、纹理特征和 Sentinel-2 时序特征，采用随机森林(Random Forests,RF)的算法提取了茶园，分类结果的准确率达到96.91%。

基于前人的研究，本研究以厦门市同安区高海拔地区的莲花镇为研究区，采用多时相的Sentinel-2遥感影像数据，通过光谱特征选择植被指数，并结合纹理特征，研究茶园提取方法，以期为茶园的管理及政策制定提供参考。

1 研究区及数据

1.1 研究区概况

厦门市(24°23'～24°54'N、117°53'～118°26'E)位于福建省东北部，总面积1700.61 km2，地形以滨海平原、台地和丘陵为主，西北部多中低山，由西北向东南倾斜。研究区位于同安区莲花镇，地处厦门西北部，最高海拔近千米，属于山区，全域植被覆盖率较高，茶树为当地主要经济作物。通过实地野外调查并结合同时期Google Earth 影像判读，研究区主要地物类型大致可分为茶园、林地、农业用地、建筑用地、裸地以及水体等，地物类型丰富。

1.2 遥感数据及预处理

数据源来自欧洲航天局研发的Sentinel 系列卫星之一的Sentinel-2。Sentinel-2 采用天体平台-L(AstroBus-L)，由Sentinel-2A 和Sentinel-2B 两颗卫星联合工作，两颗卫星运行于同一条轨道上，相位相差180°，双星座重访周期可达5 d。Sentinel-2 携带高分辨率多光谱成像仪(MSI)，采用堆扫式成像模式，拥有从可见光到短波红外的13 个光谱波段。空间分辨率高达10 m，并在红边范围具有3 个波段(表1)，为区域植被生态环境特征信息的提取分析提供了全新的解决方案[12]。

表1 Sentinel-2 多光谱影像主要参数Table 1 Main parameters of Sentinel-2 multispectral images

Sentinel-2 L2A 影像经过辐射定标和大气校正，因此采用欧洲航天局的SNAP 软件对原始影像进行重采样与格式转换。

1.3 实测数据

于2022年5月、10月对研究区进行野外考察。利用GPS 手持机和谷歌地球影像，进行野外地面样点布设及环境要素的考察与识别，标定764 个不规则多边形地块，样本分布均匀，其中80%作为训练样本，剩余作为验证样本(图1)。

图1 研究区位置及样本分布Fig.1 Study area location and the distribution of samples

2 研究方法

基于Sentinel-2 多光谱影像的光谱数据和纹理特征，采用随机森林(RF)算法进行分类，进而对茶园种植区进行提取。首先分析多时相的Sentinel-2光谱在不同地物的光谱特征差异，通过敏感光谱建立光谱指数并分析，然后利用灰度共生矩阵(GLCM)构建影像相关纹理特征，最后应用RF 分类方法对图像进行分类，并进行精度检验和分析。

2.1 分类特征提取

2.1.1 光谱特征

不同地物的光谱曲线具有独特的光谱特征，广泛应用于遥感影像分析[13]。通过Sentinel-2 影像对研究区各地物样本统计分析得到光谱曲线特征，为茶园提取提供依据。

综合考虑研究区生态茶园的物候特征以及可获取的高质量遥感数据，选择茶树春季采摘修剪期(2022年4月8日)、生长旺盛期(2022年7月22日)、秋季采摘修剪期(2022年10月20日)和生长减缓期(2022年12月9日)4 个时期的 Sentinel-2 数据中地表反射率图像，提取研究区内6 种典型地物样本(茶园、林地、农业用地、建筑用地、裸地、水体)在不同光谱波段上的光谱反射率，绘制光谱特征曲线图(图2)，并进行光谱特征分析。

图2 研究区典型地物不同时相光谱曲线Fig.2 The spectral curves of features in several times

2.1.2 光谱指数

利用光学遥感的两个波段或多个波段组合建立的植被指数不仅形式简单，而且能够有效地描述植被状态，被广泛应用于植被统计和提取、植被时空变化分析等[14]。

根据光谱在波段上的反映，选择三种光谱指数，即归一化植被指数(NDVI)、归一化差值水体指数(Normalized Difference Water Index,NDWI)和倒红边叶绿素指数(Inverted Red-Edge Chlorophyll Index,IRECI)。

NDVI 是近红外波段和红光波段光谱反射率之差与之和的比率，是反映植被生长状态的重要指示因子[15]，表达式为：

式中，ρNIR和ρred分别为近红外波段和红光波段的光谱反射率。

NDWI 是短波红外波段和近红外波段的归一化比值指数[16]，表达式为：

式中，ρSWIR为短波红外波段光谱反射率。根据典型地物在短波红外波段的表现，本研究采用的是Sentinel-2 的短波红外1 波段，即SWIR1。

Sentinel-2 在红光与近红外波段的红边区域(670～760 nm)内具有3 个波段，IRECI 包含了Sentinel-2 的3 个红边波段[17]，表达式为：

式中，ρ783、ρ665、ρ705和ρ740分别对应Sentinel-2 的Band7、Band4、Band5 和Band6 的光谱反射率。

2.1.3 纹理特征

Haralick 等[18]基于统计的方法提出灰度共生矩阵(GLCM)，可以反映图像灰度关于方向、相邻间隔、变化幅度以及快慢的综合信息。由于Sentinel-2 包含波段较多，若对每个波段都提取纹理特征图像，容易产生数据信息的冗余。因此，本研究对反射率数据采用主成分分析(Principal Components Analysis,PCA)，选取前3 个主成分。

纹理提取时窗口尺寸大小会影响GLCM性能和纹理特征完整性，若窗口较大，会降低区分不同地物的灵敏度，而窗口较小会受像元亮度不均匀性的影响。7×7 像素和9×9 像素的窗口尺寸对于茶园提取是较合适的[19]。通过实验采用7×7 像元窗口，计算纹理特征值中常用的6 个纹理：方差(Variance)、同质性(Homogeneity)、对比度(Contrast)、相关性(Correlation)、角二阶矩(Angular Second Moment,ASM)和熵(Entropy)。

由于GLCM 还具有方向性，本研究获取的纹理值是0°、45°、90°和135°四个方向的平均值。

综上所述，研究中将每个时期的3 种光谱指数和主成分变换第一、第二、第三分量的6 种纹理特征作为分类提取特征。

2.2 茶园提取方法

结合实地调研及更高分辨率的 Google Earth图像，选取茶园(183 个地块，共3796 个像元)与非茶园(423 个地块，共28838 个像元)训练样本采用随机森林(RF)算法进行分类，分类工作基于SNAP 实现。

3 结果与分析

3.1 典型地物光谱特征

将研究区6 种主要典型地物样本分为植被样本(茶园、林地、农业用地)和非植被样本(建筑用地、裸地、水体)。从图2 可以看出，植被样本和非植被样本的光谱曲线在不同时相下，同一地物的变化趋势相似，这与研究区气候特征有关，研究区植被多是常青植被。但在同一时相下，植被样本和非植被样本的光谱曲线有较大差异。

非植被样本在同一时相下，不同样本的光谱曲线趋势较为一致，且不同样本在相同波段的光谱值相差较大，如图2(b)在近红外(NIR)波段，裸地的光谱值达到5300，建筑用地的光谱值接近4000，而水体的光谱值在1100 左右。

植被样本在同一时相下，不同样本的光谱曲线趋势较为一致，呈现先上升后下降的趋势，在蓝(B)、绿(G)、红(R)波段的变化比较平缓，且光谱值比较接近；从植被红边1(RE1)波段开始，植被样本的光谱值呈明显上升趋势，茶园在近红外(NIR)波段处达最大，随后光谱值下降。林地和农业用地在不同时相下，最大光谱值的波段并不一致。在春季采摘修剪期(4月)，茶园与林地在多个波段的光谱值比较接近；在生长旺盛期(7月)和秋季采摘修剪期(10月)，除短波红外2(SWIR2)波段，茶园在其他波段的光谱值均高于其他植被，且有一定的差值；生长减缓期(12月)，茶园从植被红边3(RE3)到水汽(WV)波段的光谱值高于林地和农业用地。总体来看，茶园与其他植被具有一定光谱可分性的波段主要有 RE2、RE3、NIR 和SWIR1。

3.2 典型地物光谱指数

研究区典型地物在多时相NDVI、NDWI 和IRECI 的特征图如图3～图5所示。

图3 研究区典型地物多时相的NDVI 特征Fig.3 Characteristics of typical land cover in multi-temporal NDVI

在NDVI 特征图中，非植被类地物(水体、建筑用地)在各个时期均呈现暗黑色，植被类地物(茶园、林地、农业用地)各个时期均呈现亮白色或灰色(图3)。建筑用地附近的农业用地在7月和12月特征值较低，呈灰色。大部分林地由于较密集而呈高亮色，部分林地在坡度较大的区域也呈暗灰色。而茶园在各个时期的特征值变化不大，多呈现高亮色。与林地易混淆，无法实现很好的区分。

在NDWI 特征图中，非植被类与NDVI 特征图类似，呈暗黑色，但植被类也有部分区域呈现暗灰色(图4)。这是由于植被在短波红外波段的反射率与其叶片含水量相关，而研究区水资源较为匮乏，林地、茶园和部分农田无灌溉条件。茶树与高大林木相比较，根系较浅，在雨水较为匮乏的4月和12月，茶树冠层含水量低于高大林木，茶园所在区域灰度值较低。

图4 研究区典型地物多时相的NDWI 特征Fig.4 Characteristics of typical land cover in multi-temporal NDWI

IRECI 对叶绿素含量变化比较敏感，植被区域也呈现亮白色和灰色(图5)。在10月和12月，随着植被光合作用降低，IRECI 指数特征值下降，与NDVI 和NDWI 特征图相似，较难实现茶园和林木的区分。

图5 研究区典型地物多时相的IRECI 特征Fig.5 Characteristics of typical land cover in multi-temporal IRECI

3.3 茶园分布

图6 显示多时相下茶园的分类。总体而言，4个时期提取的茶园空间分布比较一致。研究区茶园主要分布在北部，中部地区有零散的茶园。经实地考察，研究区北部与福建省泉州市的安溪县交界，茶园种植范围广泛。受影像分辨率以及植被光谱特征相似性的限制，在研究区东南部有部分耕地被错分为茶园。

图6 研究区多时相茶园分布图Fig.6 Extraction result of tea plantations distribution

3.4 精度评价

为了验证分类精度，利用目视判读结果随机生成921 个茶园和 8321 个非茶园验证像元，通过建立的混淆矩阵，计算总体精度和总体Kappa 系数(表2)。研究区内6 种地物在4～12月4 个时期的分类总体精度依次为 94.28%、91.05%、95.53%和94.41%；Kappa 系数依次为0.8889、0.8254、0.9088和0.8881。茶园在4月和10月的用户精度较好，分别为82.92%和82.25%，这两个时期也是研究区采茶期；茶园在7月的用户精度较低，仅51.57%。总体来说，茶园提取的准确率较高。

表2 茶园提取精度Table 2 Accuracy assessment for extraction of tea plantations

4 结论与讨论

4.1 讨论

研究区植被类型具有相似的光谱特征，因此本研究基于Sentinel-2 丰富的光谱信息，结合茶园纹理特征，进行多时相的地物分类和茶园提取分析。结果表明，将光谱特征、植被指数和纹理特征结合可有效提取茶园，同时提取精度受季节因素影响。赵晓晴等[7]利用Sentinel-2 多时相影像提取茶园时，发现5月和12月是杭州市西湖区茶园提取特征的最佳组合。Sentinel-2 具有多个适用于植被监测的红边波段，李龙伟等[8]采用Sentinel-2 的红边波段构建植被指数(NDTI)作为茶园提取的分类特征。本研究也利用红边波段构造的植被指数(IRECI)作为茶园的分类特征之一，由于研究区多云雾和阴雨时期较多，本研究未使用Sentinel-2 时间序列信息进行茶园提取，但柏佳等[11]基于 Sentinel-2 与MODIS 数据融合获得完整的 NDVI 时序数据，采用随机森林算法提取茶园准确率高达96.91%。

尽管本研究结果可以有效提取茶园范围，但存在一些不确定性。如林地和农业用地与茶园具有相似的光谱特征，“异物同谱”现象增加了高精度提取茶园的难度；Sentinel-2 的空间分辨率为10 m，对于小范围的茶园，无法准确辨识。

4.2 结论

本文对厦门市同安区莲花镇茶园进行提取，基于影像的光谱特征、植被指数和纹理特征，构造随机森林分类模型，研究发现：

(1)研究区植被的光谱曲线呈现先上升后下降的趋势，而非植被的光谱曲线趋势较平稳。根据光谱特征可知，茶园与其他植被具有一定光谱可分性的波段主要有RE2、RE3、NIR 和SWIR1。

(2)在植被类型容易混淆情况下，基于随机森林分类方法提取茶园，能够较好地满足分类要求。

受多种因素影响，“异物同谱”现象仍然存在，茶园与研究区其他地物有混淆的情况，如少量林地和农地。下一步将结合高光谱数据对研究区分类精度开展更加深入的分析和探讨。