茉莉花(Jasminum sambac)最早源于古罗马帝国，在汉代通过海上丝绸之路传入古波斯和天竺，从印度随佛教传入中国已有2 200多年的栽培历史^[1]。茉莉花得益于独特的香气成为了重要的香精原料，同时也是花茶窨制的重要原材料。经过长达几千年的栽培和生产发展，茉莉花已成为我国南方山地丘陵地区农业经济发展中不可或缺的经济来源^[2]。

如何有效提取茉莉花的种植范围，巩固提升当前的乡村振兴成果，为地区产业发展规划、精准政策布局提供参考，成为当下首要解决的难题。随着近年来遥感技术的不断发展，具有覆盖面积广、时效性强、获取成本低及数据来源广等优势的卫星图像数据已被广泛应用于“三农”领域^[3-5]。如何充分利用不同遥感影像数据集之间的优势，实现对茉莉花种植区的遥感识别提取成为当前研究的主要方向^[6]。徐伟燕等^[7]基于资源三号卫星数据，使用归一化植被指数(NDVI)时像差异及光谱特征，结合决策树及神经网络分类方法对平原地区茶树种植区进行提取；马超等^[8]运用MODIS数据及Landsat数据，在考虑光谱、物候、植被指数特征下，采用决策树分类方法对南方丘陵山区茶园进行提取；Akar等^[9]考虑到茶园和榛子园与周围植被的光谱反射率值相似的特性，使用NDVI结合灰度共生矩阵(GLCM)和Gabor滤波器的方法对土耳其黑海东部的茶园和榛子园进行范围提取；王斌等^[10]基于Landsat 8 OLI影像获取不同季节影像数据的光谱、纹理、物候三大特征参数，分析优选不同季节影像特征对浙江省安吉县茶园进行提取。当前茶园提取研究已由单一遥感数据源扩展到多源遥感数据源相融合的发展趋势，进一步发挥了不同遥感数据图像信息的优势，但是如何充分发挥当前遥感图像信息对茉莉花种植区进行时空信息提取的优势仍是当前急需解决的问题^{[11, 12]}。

遥感时空数据融合技术是对“高空间-低时间”与“高时间-低空间”分辨率的遥感卫星影像进行融合运算及处理，综合两者高空间、高时间的优势，解决遥感卫星数据“时-空矛盾”问题的研究技术^[13]。针对Landsat数据与MODIS数据相互融合的问题，Gao等^[14]在2006年提出时空自适应反射率融合模型(STARFM)，这也是当前影响最大、运用最为广泛的时空数据融合模型；随后Zhu等^[15]于2010年在原有STARFM的基础上提出增强型时空自适应反射率融合模型(ESTARFM)，有效提升了时空数据融合的精度，并且在农作物监测、土地利用/覆盖、植被分析等领域中得到广泛应用。Chen等^[16]基于ESTARFM得到高空间分辨率的涡度通量数据，提高了区域景观尺度下的常绿针叶林生物量估算精度；柳文杰等^[17]通过ESTARFM融合环境星数据与MODIS数据结合物候特征参数对洞庭湖水稻种植区域进行提取；赵亚杰等^[18]基于ESTARFM融合Landsat 8数据与MODIS数据，构建NDVI时序数据及水稻物候特征参数，对江汉平原地区水稻种植区域进行提取研究。遥感时空数据融合技术能够有效提高图像数据质量，提升图像信息提取精度^[19]。

为此，本研究以广西横州市为研究区，基于Landsat和MODIS遥感数据对横州市茉莉花种植区进行提取，充分挖掘其图像时序物候信息、几何空间、纹理特征，为解决遥感数据时空信息的缺失、无法有效监测农作物生长的关键物候信息以及遥感图像地表空间异质性等问题提供切实可行的技术方法，同时也为地区特色农业作物种植区提取提供参考方案，为促进地区乡村产业振兴发展、制定农业产业规划提供科学支持和技术支撑。

1 研究区与数据 1.1 研究区概况

横州市位于广西壮族自治区东南部，南宁市东部，享有“中国茉莉花之乡”“世界茉莉花都”的美誉，地处东经108°48′-109°37′，北纬22°08′-23°30′，东临广西贵港市，南接广西钦州市灵山县，北部与南宁市宾阳县接壤，辖区土地总面积3 464 km²(图 1)。土地类型主要以丘陵、平原为主，地势四周高、中部低，形似一个广阔的盆地。受地形及地质构造的影响，土壤多以沉积、冲积形成，尤其在中部广阔平缓地带，土层厚度大，土壤熟化程度高，且多年平均气温稳定在22 ℃，年降水量稳定在1 600-1 700 mm，气候舒适宜人，日照充足，良好的气候条件孕育着丰富的农业资源，是南宁乃至广西的农产品种植基地。横州市种植的农作物主要有甘蔗、水稻、玉米、茉莉花等，其中茉莉花是横州市最负盛名的经济作物，茉莉花种植产量占全国总产量的80%，占世界总产量的60%，是中国最大的茉莉花生产和茉莉花茶加工基地^[2]。

1.2 研究数据与预处理

研究区使用的遥感数据为Landsat 8 OLI数据和中分辨率成像光谱仪数据产品MOD13Q1(16d_NDVI)。其中，覆盖研究区Landsat 8数据的行列号为125/44，研究选取了时间跨度在2019年12月1日-2021年4月1日内图像质量最好、少云的7景影像，空间分辨率为30 m，包含红、绿、蓝及近红外等波段；MOD13Q1数据的行列号为h28v06，遥感数据时间跨度为2019年12月1日-2021年4月1日，空间分辨率为250 m，数据经过16 d的最大值合成，共计23处影像，MODIS数据主要利用时间采集上的优势与Landsat数据进行融合，构建年度内16 d/景30 m×30 m的NDVI数据，用于提取主要作物物候信息。本研究所获取的数据来源于美国地质调查局(USGS)、美国航空航天局(NASA)，遥感图像数据成像时间及配对信息如表 1所示。

对获取得到的Landsat 8多光谱数据在ENVI5.3软件平台上分别进行辐射校正(Radiometric calibration)和大气校正(FLAASH)，去除大气、土壤、传感器角度等对影像的影响，然后裁剪研究区范围，并进行波段计算得到NDVI数据。MOD13Q1数据处理则使用ENVI5.3平台拓展工具MTCK插件进行格式转换和重投影，投影坐标系与Landsat 8数据一致，随后以Landsat 8数据为基准对MOD13Q1数据进行几何精校正。最后，提取出MOD13Q1数据中的16d_NDVI波段，并重采样使其空间分辨率与Landsat 8数据一致(30 m)，以便在后续的ESTARFM中进行时空数据融合，重建年度NDVI数据。

2 研究方法

本研究主要分为以下步骤：首先采用ESTARFM重建具有完整年度时间序列的ESTARFM-NDVI数据；其次对重建得到的ESTARFM-NDVI数据进行精度评价，再引入Savitzky-Golay (S-G)滤波对时序NDVI数据进行平滑重构，消除数据噪声的影响，进而获取物候特征信息，基于GLCM获取研究区遥感影像的8个纹理特征参数；最后结合不同数据集构建分类方案，采用随机森林(Random Forest，RF)分类方法对横州市茉莉花种植区进行提取，并以2020年横州市野外实地采样数据结合天地图高分辨率影像数据作为精度检验样本。本研究的技术路线如图 2所示。

2.1 时空数据融合模型

本研究使用ESTARFM对Landsat和MODIS数据进行融合，实现年度内长时间序列NDVI数据重建。

ESTARFM的主要思路与STARFM相似，该算法是在Gao等^[14]提出的STARFM的基础上进行改进，考虑到地表地物类型的复杂性，且低分辨率图像数据像元覆盖面积大，往往存在大量的混合像元，因此需要考虑邻近像元之间的相似反射率信息。为此，本研究使用Gao等^[14]提出的移动窗口法对相邻像元的反射率信息进行预测，在充分考虑地表地物复杂性的情况下实现对高分辨率影像反射率的计算，计算公式如式(1)所示。

$ \begin{array}{r} h\left(x_{\frac{w}{2}}, y_{\frac{w}{2}}, T_p, P\right)=h\left(x_{\frac{w}{2}}, y_{\frac{w}{2}}, T_0, P\right)+ \\ \sum\limits_{i=1}^n W_i \times V_i \times\left(l\left(x_i, y_i, T_p, P\right)-l\left(x_i, y_i, T_0, P\right)\right), \end{array} $

(1)

其中，h为高分辨率影像数据；x, y为影像像元在图像中的位置；l为低分辨率影像数据；T₀、T_p表示像元反射率变化的时刻；P表示相似像元；w表示移动窗口的大小；n表示相似像元的个数；W_i为第i相似像元的权重；V_i表示为第i相似像元的转换系数^[15]。

2.2 年度NDVI数据重构

卫星传感器在接收遥感数据的过程中，受传感器自身性能、角度及大气云层、水汽等因素的影响，导致NDVI时间序列数据易受到噪声的影响，从而干扰识别植被物候信息。因此，在利用NDVI时间序列数据获取植被物候信息前，需要对其进行拟合去噪，重构NDVI时间序列数据^[20]。本研究采用TIMESAT软件中的S-G滤波算法对ESTARFM-NDVI时间序列数据进行拟合重构。S-G滤波算法在对长时序数据进行滤波平滑处理的同时，能够有效保留数据变化的信息，适用于植被物候观测变化研究^[21]。其表达式为

$ N^*=\frac{\sum\limits_{i=-n}^{i=n} Z_b N_{a+b}}{M}, $

(2)

其中，N^*为NDVI拟合值；i为参与平滑去噪的影像数量；N_a+b为NDVI像元原始值a到第b个的累加值；Z_b为第b个NDVI值进行滤波的系数；n为1/2个滤波器窗口的宽度；M为滤波器长度，等于滑动数组的宽度的2m+1^[22]。

植被物候通常会经历从生长开始到结束的生长规律，而平滑后的Landsat-NDVI时间序列曲线能够有效反映研究区植被物候的动态变化特征。因此，本文基于重构后的Landsat-NDVI数据，采用TIMESAT软件提取植被物候参数特征^[23]，选取生长季开始期(SOS，a)、生长季结束期(EOS, b)、生长季持续时期(LOS, c)、生长季拟合最大值(BOS, d)、生长季振幅(AOS, e)、生长季最大值积分(MOS, f)共6个物候参数特征，具体植被物候参数阶段及特征定义如图 3所示。

2.3 遥感图像纹理特征获取

纹理是复杂的视觉实体或子像元的组合，具有亮度、颜色、陡度、大小等特征。卫星影像的混合像元辨识主要通过纹理特征对其进行区分，考虑到遥感图像“异物同谱、同物异谱”的现象，利用GLCM的纹理提取方法，筛选分化程度较高的波段和图像纹理方法作为图像纹理特性。因此，本研究将遥感图像纹理特征作为分类信息特征变量，选取Landsat 8红光、绿光、蓝光及近红外波段为数据源，采用基于GLCM获取研究区8个纹理特征，包括均值(Average)、方差(Variance)、同质性(Homogeneity)、对比度(Contrast)、非相似性(Dissimilarity)、熵(Entropy)、角二阶矩(Angular second moment)和相关性(Correlation)^[24]，计算公式(3)-(10)如表 2所示。

2.4 分类方案与精度分析

随机森林是基于集成学习的思想将多棵决策树集成的一种算法。通过在原始训练样本中随机有放回地抽取多个训练样本，以不同训练样本集来构造决策树，并对多个决策树的结果进行整合投票，最终解决数据的分类及回归问题。随机森林分类通过随机有放回地构建训练样本，预测特征变量集，进而构建决策树的算法思路，具有稳定性高、学习过程迅速、抗噪能力强等优势，适合于复杂环境下的植被信息分类提取^{[25, 26]}。为此，本研究采用随机森林分类算法对横州市茉莉花种植区进行提取，并设计了3种分类方案：①使用时序NDVI数据、物候特征组合分类；②使用时序NDVI数据、纹理特征组合分类；③使用时序NDVI数据、物候特征、纹理特征组合分类。

结合研究区野外监测实际状况，将研究区主要的土地利用类型分为林地、水体、建设用地、甘蔗地、水稻地、茉莉花地、玉米地等7类。以实地调查数据为基准，实地采集的样本均为大范围连片区域，平均采样面积大于30 m²。结合天地图、横州市土地利用现状图选取地类样本(表 3)，其中70%的样本为分类训练样本，30%的样本为分类精度检验样本，生成混淆矩阵，检验分类后的数据的总体分类精度、生产者精度、用户精度及Kappa系数，并进行定量分析评价。

3 结果与分析 3.1 融合时序NDVI数据重构

图 4为研究区典型地物重构NDVI后的曲线与原始NDVI曲线。从图 4可以明显看出，经过S-G滤波拟合重构后的时序NDVI数据能够很好地消除外界因素引起的突变，能够更好地区分不同地物类型的物候特征。茉莉花经过2-3月的剪枝后，4月份起逐渐抽枝、出芽、长叶，6月份茉莉花生长最为茂盛，随后开始进入花期，花期可持续到10月中旬。其中在2-3月剪枝期、4月抽枝出芽期与其他地物时序NDVI曲线差异显著。因此，选择影像期数为4-10 (对应第64-160天)内的数据作为茉莉花种植区提取的最佳时序NDVI数据。

通过时空数据融合模型得到横州市23期时序NDVI数据集，为验证影像数据融合效果，本研究对2020年4月22日MOD13Q1数据融合生成的ESTARFM-NDVI影像与2020年4月27日卫星获取的Landsat-NDVI影像进行对比分析，并在此基础上对真实NDVI与预测NDVI数据两处局部细节进行对比。从目视判别分析上来看，两者整体在空间上差异性不大，空间上具有高度相似性(图 5)。对研究区整体与局部细节进行相关系数R检验(图 5)可知，整体和局部区域的真实NDVI与融合NDVI的相关系数能达到0.83以上，表明真实NDVI与融合NDVI具有较高的相关性。因此，基于时空数据融合得到的ESTARFM-NDVI数据能够较好地预测真实植被的NDVI信息，可以满足研究区茉莉花种植区的信息提取。

3.2 茉莉花物候参数及纹理特征

通过对比研究区(剔除林地、水体、玉米地、建设用地)内地物类型的物候参数(图 6、表 4)可知，生长季开始期，茉莉花SOS参数值略高于甘蔗；生长季结束期，茉莉花EOS参数值仅高于水稻；生长季持续时期，茉莉花LOS参数值仅高于水稻；在基准值中，茉莉花BOS参数值最高；在生长季最大值积分中，茉莉花MOS参数值是所有地物类型中最高的；在生长季振幅中，茉莉花AOS参数值最低。因此，研究所选用的6类物候参数中，茉莉花物候参数与其他地物类型存在明显的差异，表明物候参数能够有效区分出茉莉花种植区与其他地物类型。

3.3 茉莉花种植区提取结果及精度分析

3种随机森林分类方案的提取结果如图 7所示。方案1的空间分布特征基本与研究区的分布概况吻合，但是在细节上表现出明显的漏分、少分现象，分布不均匀，部分玉米、甘蔗等地物被错分为茉莉花，存在较多孔隙，空间不连续、完整[图 7(a)]；方案2在空间上表现出明显的连续性，表明纹理特征能更好地表征图像在概率统计下的信息值，在分类程度上，混淆程度大，部分林地被分类为茉莉花种植区，错误程度明显高于方案1[图 7(b)]；方案3充分考虑了物候和纹理的优势信息，在空间分布上与实际情况较为吻合，局部细节上能够较好地将易混淆的地物分类出来，错分及漏分现象明显减少[图 7(c)]。结合实地调查、资料收集及天地图高清影像分析发现，茉莉花主要种植在横州市中东部的平原地带，其中校椅镇、云表镇等地种植范围最为广阔，西部乡镇则少有种植，分布特征不明显。综上可知，3种分类方案对识别提取横州市茉莉花种植区的分布特征各具优势，基本能够表述横州市茉莉花种植区的空间分布特征。

为对比探究不同分类方案下的分类精度，本研究通过混淆矩阵对3种分类方案进行评价，分类精度评价结果如表 5所示。由表 5可知，方案3的分类精度明显优于方案1和方案2，其总体分类精度达到88.36%，Kappa系数为0.863 3，表明该方案精度能够达到分离茉莉花种植区的效果。

综合而言，物候参数能够修正NDVI时序数据在茉莉花分类过程中引起的错分、漏分现象；纹理特征能够对无物候地物进行有效识别区分，并消除提取过程中的孔隙及孔洞现象，达到提高分类精度的目的。为此，将本研究得到的结果与2020年横州市茉莉花种植的统计年鉴面积进行对比分析(表 6)，结果显示，本研究采用的方案对横州市茉莉花种植区提取的结果与实际分布结果基本相符，总面积统计相差约1.02万亩(1亩≈666.67 m²)。

4 讨论

本研究提出一种结合3种不同特征数据集的茉莉花种植区遥感识别方法，针对中低纬度地区遥感影像数据受云、雨等因素影响的问题，引入时空数据融合ESTARFM算法重建年度NDVI时序数据，构建研究区物候特征和纹理特征数据集，并结合拟合重构后的NDVI时序数据集，对比探究不同数据集对茉莉花种植区遥感识别提取的精度，充分发挥不同遥感数据集在茉莉花种植区提取的优势，拓展遥感数据在农业领域的应用研究价值。本研究表明物候参数与纹理特征可以有效提高复杂环境下的茉莉花种植区提取精度，其中物候参数能够将茉莉花种植区与不同作物种植区加以识别区分，更好地减少数据错分、漏分等问题；纹理特征则依据遥感图像上反映同种类型地物的同质化特征信息，在参与分类过程中能够使茉莉花种植区提取结果更加平滑、连片且无孔隙，空间表达上凸显了连续性。

虽然本研究通过时序NDVI数据、物候参数和纹理特征开展茉莉花种植区提取可以有效获取其空间分布特征，但是仍存在局限性：一是时空数据重建所使用的30 m分辨率数据较少，后期研究可以考虑异源遥感数据(如HJ-1A/1B)进行插补或者考虑更高分辨率的图像(如Sentinel-1/2、GF系列)进行数据融合重建；二是在分类方法上，本研究仅采用一种机器学习分类方法对茉莉花种植区进行提取，单一的机器学习分类方法可能会使结果存在一定的局限性。综上所述，本研究采取的重建时序NDVI数据结合物候参数及纹理特征的茉莉花种植区提取方法，取得了良好的提取效果，同时也对未来的研究方向提供了思路，今后的研究可以考虑对多类型特征参数进行重要性优先分析，引入不同分类方法模型进行对比分析，进一步提高茉莉花及其他作物的分类提取精度，以满足土地利用及生态产业调控等要求。

5 结论

本研究以Landsat和MODIS数据为基础，利用ESTARFM重建年度NDVI时序数据，采用S-G滤波进行平滑去噪并提取研究区物候参数，采用随机森林分类方法结合时序NDVI、物候参数以及纹理特征对横州市茉莉花种植区进行提取，得到以下结论。

① ESTARFM重建得到的ESATRFM-NDVI与真实影像获取的NDVI拟合性较高。通过对研究区真实NDVI与融合NDVI拟合对比可知，相关系数达到0.83以上，进一步验证了采用ESTARFM能够较好地重建南方多云雨山地丘陵区的缺失影像数据。

② 重建生成的NDVI数据集难免会受到外界各类噪声的干扰，本研究利用S-G滤波对出现的异常值进行平滑重构，消除噪声的同时还能更好地展现不同地物类型的时序NDVI变化，进一步提高物候特征参数提取的准确性。

③ 本研究采用的方案在茉莉种植区提取中引入物候参数和纹理特征，与仅使用时间序列NDVI数据集的分类方案相比，有效提高了提取精度，提取精度达到85.09%，总体分类精度达到88.36%，Kappa系数为0.863 3，表明物候参数与纹理特征可以有效提高南方山地丘陵地区茉莉花种植区的分类效果。