0 引 言
随着现代农业的发展，蔬菜生产呈现出集约化与高复种率的特征，一年多季的种植模式不仅提高了土地利用效率，也在稳定供应和保障“菜篮子”安全中发挥了关键作用。适度提升复种率有助于促进蔬菜产业发展，而过度复种则可能引发土壤肥力衰退，从而制约农业的可持续发展［1］。耕地复种指数（Cropping Index, CI）作为衡量耕地利用程度的重要指标，对评估农业生产潜力与可持续发展水平具有重要意义。
与传统基于农业普查数据的方法相比，遥感技术因其高效性、大规模监测和低成本等优点，被广泛应用在精准农业监测中［2］，成为耕地复种指数提取的普遍手段。国内外学者已在全球［3］及区域尺度（如欧洲［4］、南亚［5］、湄公河流域［6］），以及中国全国范围［7］、全国粮食主产区［8］、黄淮海地区［9］、湖南省［10］、华南水稻种植区［11］等区域开展了不同尺度的耕地复种指数监测研究。相关研究多聚焦于玉米、水稻、小麦等传统作物，少有研究针对蔬菜这类高复种作物。所用数据以MODIS［12, 13］、Landsat［14］、Sentinel系列［15］等的卫星数据为主，主流方法可以概括为峰值点探测法［16］、生长周期判断法［17, 18］、小波变换法［19］、时域混合模型分解法［12, 20］，蔬菜作物因其高复杂的特性，是否可以直接套用粮食作物复种指数提取的流程和方法有待考证，特别是在时序曲线重构方面。
云南省玉溪市通海县依托得天独厚的气候条件，不断完善基础设施配套，建立较为完善的市场体系，规模化、标准化、专业化、科技化水平持续提升，已发展成为云南省重要的蔬菜生产基地和西南地区主要的蔬菜集散地。长期以来，当地形成了以露地蔬菜为主、设施栽培为辅的多茬高复种生产格局，可全年轮作3~4茬不同品种蔬菜，4季外销、全年无休。对该地区开展复种指数监测，不仅有助于准确量化多季种植格局、揭示耕地利用强度及其变化趋势，也为农业生产布局优化、蔬菜稳产保供、生态环境保护和区域农业可持续发展提供科学支撑。
因此，本研究针对通海县这一典型露天蔬菜种植县，使用了一种自动识别覆盖研究区最优的网格和轨道编号的方法，基于所识别出特定的网格和轨道编号的长时序的Sentinel-2卫星影像数据，采用Whittaker Smoothing（WS）对归一化植被指数（Normalized Difference Vegetation Index, NDVI）时间序列曲线进行平滑重构，结合二次差分法提取时间序列曲线上的波峰，构建了通海县耕地复种指数遥感提取的方法，提取了2020—2024年通海县耕地复种指数，并分析了其空间分布和时空演变特征，为通海县耕地资源管理和可持续利用提供理论和技术参考。本研究主体工作于Google Earth Engine（GEE）中完成，该平台支持海量时空数据的快速获取与高效计算，显著降低了大尺度、长时间序列监测的技术门槛，众多学者基于此平台开展了有关农业方面的研究［21, 22］。
1 研究区及数据来源
1.1 研究区概况
本研究选取云南省玉溪市通海县（23°55′11″~24°14′49″N，102°30′25″~102°52′53″E）为研究区（图1），其位于玉溪市东部，杞麓湖畔的湖积平原上，区域总面积721 km2。县域地形整体呈南低北高格局，地貌类型主要包括盆地、中山和河谷，海拔范围1 350~2 443 m。研究区属亚热带半湿润高原季风气候区，多年平均降水量898 mm，平均气温15.70 ℃，年均日照时间2 273.80 h。具有良好的农业生产基础，且农业规划良好，地块形状规则、面积较大，为遥感监测提供了有利条件，县域内耕地分布广泛，全县耕地面积112.60 km2，常住人口28.97万人，按农业人口平均473.33 m2，为云南省典型的露天蔬菜种植县，蔬菜轮作制度显著。2024年蔬菜总产量915 000 t，被列为“国家农业可持续发展试验示范区”（资料来源：通海县人民政府网站，综合发布于2024—2025年）［23］。
1.2 数据来源及预处理
1.2.1 遥感影像数据
本研究采用由欧空局发布的Sentinel-2MSI（SR）2020—2024年影像共362景，该卫星提供从可见光到短波红外的13个光谱波段，已被广泛应用于农业、林业、水文、城市与环境变化等的各类研究中。根据研究区的范围在GEE中自动筛选出研究区内网格编号为47QRG、轨道号为61的影像，利用SCL波段进行云掩膜处理，进一步计算NDVI和APGI（Advanced Plastic Greenhouse Index），生成5 d间隔、10 m空间分辨率的时间序列数据，用于复种指数提取及塑料大棚识别。2020—2024年去云后耕地像元的有效观测次数如图2所示，研究区受云遮挡导致数据缺失严重，但在对其进行耕地掩膜后，发现数据缺失严重区域多为山地，耕地区域数据质量良好，约80%的耕地像元有效观测次数大于140次。另外，本研究对年均每两月的有效观测次数进行统计，如图2c所示，7~10月因云雨的影响数据缺失较为严重，但1~6月、11~12月有效观测次数良好。
1.2.2 实地调查验证数据集
复种指数验证数据主要通过实地走访农户调查获取，本研究于2025年9月开展实地调查，走访农户100余人，询问地块2024年的种植次数，得到338例验证点，其中单季种植耕地0例、双季种植耕地35例、3季种植耕地149例、4季种植耕地146例，5季种植耕地8例，主要分布于通海县杞麓湖畔周围的城区，如图3所示。
1.2.3 土地利用数据集
本研究采用了GLC_FCS10［22］土地覆盖数据集，这是1份全球10 m分辨率、30类精细土地覆盖产品，可在以下网址中免费获取（https://zenodo.org/records/14729665）。对于耕地，GLC_FCS10将其细分为草本覆盖耕地、树木或灌木覆盖耕地、灌溉耕地3种，本研究将这3种耕地都纳入到耕地的范围中，其与ESA WorldCover、ESRI LC及其他10 m或30 m产品相比，耕地、森林及湿地等类别表现更优，具有更高的精度。
2 研究方法
在生成年度复种指数产品的数据处理和工作流程中，首先使用线性插值填补了时间序列数据中空缺的部分，并采用WS对时间序列曲线进行平滑重构，运用二次差分法结合物候信息提取时间序列曲线中峰值的个数，进而确定复种次数，得到年度耕地复种指数产品，最后，通过实地调查数据对结果进行精度验证，并结合混淆矩阵和Kappa系数评估复种指数提取的准确性。
2.1 遥感数据筛选方法
本研究使用了一种自动识别覆盖研究区最优的网格和轨道编号的方法，该方法流程如下：1）根据研究区筛选指定时间内的所有影像，并提取影像的网格和轨道编号；2）计算研究区面积、不同网格和轨道编号影像与研究区相交的面积，以此计算出不同网格和轨道编号影像填充该研究区的面积占比；3）根据贪心算法的逻辑，在每一步迭代中优先选取对剩余区域贡献覆盖面积最大的要素，并不断更新未覆盖区域，直至达到完全覆盖或无增益为止。最后将完全覆盖研究区的网格和轨道编号视为最优。因本研究区较小，仅1次迭代即网格编号为47 QRG、轨道号为61的影像满足完全覆盖研究区，从而选择该网格及轨道编号的影像作为原始数据构建时间序列曲线。
2.2 时间序列曲线重构方法
式中： 为目标函数，由保真度项和粗糙度项组成； 为保真度，表示平滑曲线与原始数据之间的偏差； 为粗糙度，表示平滑曲线的平滑程度； 为原始数据； 为平滑后的数据； 为差分算子， 为平滑参数， 越大曲线越光滑，越接近0曲线越贴近原始数据。
2.3 活跃耕地识别
为了剔除耕地数据中潜在的非作物种植耕地，本研究基于植被-土壤-色素时序耦合度和后向散射系数时序变异性指标（Vegetation-Soil-Pigment Indices and Synthetic-Aperture Radar （SAR） Time-Series images, VSPS）区分并识别了活跃耕地［26］，该方法运用3个指标：第1个是基于Sentinel-1雷达数据后向散射系数VV波段的时间方差（Temporal Variance of VV Backscatter Coefficient, VVV）；第2个指标是基于Sentinel-2数据提出的植被-土壤差异时间变化（Vegetation-Soil Difference in Temporal Variation, VSD），揭示了植被变化与土壤覆盖之间的同步程度；第3个指标是根据Sentinel-2的红边波段提出的植被-色素差异时间变化（Vegetation-Pigment Difference in Temporal Variation, VPD），用于表征植被与色素指数之间的相对变化趋势。另外，通海县虽以露天种植蔬菜为主，但少部分仍采用塑料大棚模式种植，本研究依据APGI指数识别并排除了塑料大棚［27］，由于塑料大棚覆盖影响光谱反射特征，其作物生长情况难以通过遥感反射率产品有效监测。
2.4 复种指数提取算法
本研究采用二次差分法识别峰值提取作物的耕地复种指数，该方法的原理如下：将每个像元的时间序列植被指数看作1个数列，即N景影像的植被指数对应N个数列元素，依次用后1个元素的值减去与其相邻的前1个元素的值，形成包含N-1个元素的新序列。同时，对新序列重新赋值，负数赋值为-1，正数赋值为1，完成一阶差分；对一阶差分得到的结果进行第2次差分，与一阶差分做法类似，用后1个元素减去与其相邻的前1个元素的值，得到包含N-2个元素的新序列，新序列中值为-2，且前后位置的值均为0的点即为极大值点。
尽管时序曲线经过了WS的平滑重构，但所提取出的峰值仍包含一些伪峰（非作物波峰），存在伪峰的原因有两个：其一是非作物生长期的杂草，以及其他草本植物的生长，会出现非作物的波峰；其二是平滑重构前NDVI如果连续较低或较高时，细微波动在平滑后可能被放大，从而生成非真实的峰值。
统计排除伪波峰后耕地像元波峰的数量，依据公式（4） 计算了研究区内耕地复种指数。
式中：M表示研究区内统计出的峰值总数量； 表示研究区内的耕地总像元数； 表示研究区内平均耕地复种指数。
2.5 准确性评估
本研究采用定量评估的方法评估复种指数提取结果的准确性，定量评估基于研究区内走访农户调查的338个样本点，以此构建了误差矩阵，计算了用户精度、生产者精度、总体精度、Kappa系数、平均绝对误差（Mean Absolute Error, MAE），以及均方根误差（Root Mean Square Error, RMSE）。
3 结果与分析
3.1 时间序列曲线重构结果
基于特定网格和轨道编号的时间序列曲线重构结果如图4a、图4d、图4g、图4j所示，不同种植制度的NDVI时间序列重构曲线的作物生长周期显著，能够反映作物的生长状况，波峰特征明显，便于进行下一步的波峰提取。与常规方法（按10 d以最大值或中值复合影像）相比，本研究所采用的方法显著更优，常规方法因影像数量减少，在波峰间隔时间较短时重构曲线的波峰不明显，如图4e、图4f所示，考虑平滑参数λ与之不适配，调整平滑参数为λ=2后，如图4h、图4i、图4k、图4l所示，3季和4季作物的重构曲线依然存在此问题，可以得出本研究基于特定网格和轨道编号的时间序列曲线重构方法对于复种频率高波峰时间间隔短的情况具有更好的表现。
3.2 定量复种指数评估
本研究所使用的方法的定量评估如表1所示。整体准确率为89.94%，95%置信区间为86.70%~93.20%，Kappa系数为0.84，95%置信区间为0.82~0.85，MAE为0.11，RMSE为0.36。在生产者精度方面，双、3、4、5季种植耕地分别为94.29%、93.96%、86.30%、62.50%；用户精度方面，双、3、4、5季种植耕地分别为84.62%、89.17%、93.33%、83.33%，其中4季种植耕地的生产者精度较低，是由于4季作物波峰时间间隔小，其真实波峰的位置会受云遮掩导致波峰时间提前或延后，造成两真实波峰的时间间隔不足60 d，从而在排除伪波峰环节会出现少量误判。单季和5季种植耕地由于其类别比例过低（<3%），涉及样本数量太小，其精度参考价值不大。从验证结果看，利用WS对NDVI时间序列曲线平滑重构并运用二次差分算法进行通海县耕地复种指数提取是可行的且精度良好。
表1 通海县耕地复种指数提取精度验证结果Table 1 Accuracy validation results of cropping index extraction in Tonghai county |
| 种植制度 | 单季 | 双季 | 3季 | 4季 | 5季 | 生产者精度/% |
|---|---|---|---|---|---|---|
| 总体精度（OA） | 89.94%， | Kappa系数 | 0.835 | |||
| 单季 | 0 | 0 | 0 | 0 | 0 | / |
| 双季 | 1 | 33 | 1 | 0 | 0 | 94.29 |
| 3季 | 0 | 3 | 140 | 6 | 0 | 93.96 |
| 4季 | 0 | 3 | 16 | 126 | 1 | 86.30 |
| 5季 | 0 | 0 | 0 | 3 | 5 | 62.50 |
| 用户精度/% | / | 84.62 | 89.17 | 93.33 | 83.33 | / |
|
3.3 通海县2020—2024年耕地复种指数空间分布分析
基于上述耕地复种指数提取方法，对2020—2024年通海县耕地复种指数进行了提取并进行空间分布绘图，结果如图5a~图5e所示，2020—2024年通海县耕地复种指数分别为221.45%、217.80%、275.37%、232.41%、237.50%，其中各年份不同种植制度耕地占比如图5f所示。2020、2021、2023、2024年都以双季种植为主，3季种植次之，2022年以3季种植为主，双季种植次之。从地域分布上看，多季种植（≥3）多分布于沿杞麓湖的城区，核心原因是：1）该区域水资源丰富，自然条件优渥；2）该区域耕地范围集中，能够开展大规模集中作业；3）城区交通便利，相关市场和产业链都集中在城区；多种有利条件共同支撑了该区域蔬菜多季种植，而分布在其他位置的小宗耕地，多种植传统粮食作物，复种指数就会明显低于以蔬菜种植为主的城区。
3.4 通海县2020—2024年耕地复种指数年度变化分析
2020—2021年，复种指数下降3.67%，其中复种指数下降、不变和上升的耕地面积分别占总耕地面积的31.08%、40.23%和28.69%，其中10.45%的耕地由3季转为双季，占主导地位。
2021—2022年，复种指数大幅上升57.57%，下降、不变和上升的耕地面积分别占13.79%、33.02%和53.18%，其中17.60%的耕地由双季转为3季，占主导地位。
2022—2023年，复种指数下降42.96%，下降、不变和上升的耕地面积占比分别为47.96%、35.06%和16.99%，其中16.62%的耕地由3季转为双季，占主导地位。
2023—2024年，复种指数小幅上升5.09%，下降、不变和上升的耕地面积占比分别为27.49%、39.19%和33.32%，其中11.21%的耕地由双季转为3季，占主导地位。
4 讨 论
4.1 与现有耕地复种指数算法相比
时序曲线构建方法：对于遥感时间序列曲线的构建方法，WS因涉及矩阵运算要求参与运算的所有影像空间范围一致。常规方法通过影像复合（如10 d或月尺度的中值/最大值）来构建时序曲线，这削弱了Sentinel-2的5 d重访周期的优势，同时，对于蔬菜作物此类复种频率高、波峰时间间隔短的情况，复合影像会导致重构后的时间序列曲线波峰特征减弱，进而影响复种指数识别的准确性。本研究根据研究区自动识别了覆盖研究区的网格和轨道编号，基于所识别出特定的网格和轨道编号的长时序的Sentinel-2卫星影像数据以5 d时间分辨率平滑重构得到了时序植被指数曲线，更好地保留了作物生长的动态变化特征。与直接按行政区划选择覆盖研究区影像的方法相比，该策略有效避免了1 d内多景不同网格及轨道编号影像同时覆盖同一像元所造成的辐射差异问题。由于不同景影像的拍摄位置和角度存在差异，后续复合会导致同一位置选取到不同观测角度下的波段值，所构建的时序曲线质量受到一定程度的影响。本研究通过保持网格与轨道编号一致，确保影像拍摄位置与角度相同，在保留构建时序曲线所必要的影像的基础上，也减少了数据冗余与计算复杂度。考虑研究区不能被单一网格和轨道完全覆盖的情况，识别网格和轨道编号的方法会根据该影像在研究区内的覆盖面积占比排序，从大到小依次输出网格和轨道编号的组合，直到完全覆盖研究区为止。这也就导致大范围的研究区无法适用该方法，因为范围越大，完全覆盖研究区的网格和轨道编号的数量越多，此操作就变得没有意义，故只适用于小尺度（县级尺度）研究。
4.2 通海县耕地复种指数特性
与本研究结果相比，现有产品显著低估了通海县的复种水平。这是由于全国尺度产品在算法设计上主要面向粮食作物区，缺乏对露天蔬菜多季轮作制度的针对性考虑，导致难以反映该类地区的真实耕作特征。
因此，本研究结合通海县的种植制度特点，利用Sentinel-2时序影像开展了复种指数监测，构建了2020—2024年通海县的复种指数产品。该成果更符合通海县实际耕作状况，可为后续农业精细化管理与可持续发展研究提供科学支撑。
5 结论与展望
本研究根据研究区自动识别了覆盖研究区最优的网格和轨道编号，基于所识别出特定的网格和轨道编号的长时序的Sentinel-2卫星影像数据，结合WS平滑重构、APGI指数与VSPS算法，排除了非活跃耕地和塑料大棚干扰，采用二阶差分法以10 m的空间分辨率提取了2020—2024年通海县的耕地复种指数，系统揭示了其时空演变特征。结果表明，通过对2024年实地调查结果的验证，通海县耕地复种指数提取总体精度达89.94%，Kappa系数为0.84，MAE为0.11，RMSE为0.36。通海县耕地复种指数从2020年221.45%先提升到2022年275.37%后降低到2024年237.50%，其中耕地复种指数的变化由双季与3季种植之间的互相转换为主导。2020、2021、2023、2024年都以双季种植为主，3季种植次之，2022年以3季种植为主，双季种植次之。从地域分布上看，多季种植（≥3）多分布于沿杞麓湖的城区，该区域水资源丰富，为蔬菜等农作物露天种植提供了必要条件。本研究成果可为区域农业种植结构优化与耕地利用强度评估提供可靠的数据支持，并为政府制定与调整农业相关政策提供科学依据。
尽管通海县为典型的露天蔬菜种植区，但蔬菜多分布于杞麓湖畔周边城区，种类繁多，而其他区域仍以粮食作物为主。由于 Sentinel-2 数据不足以区分不同类型蔬菜，且实地调查发现部分农户在蔬菜收获后种植玉米等粮食作物，导致蔬菜与粮食作物难以按地块精确划分，因此本研究未对蔬菜与粮食作物进行区分，所得到的复种指数为通海县整体的复种指数，包括了所有的作物。这一局限在一定程度上限制了对特定作物种植制度的深入分析。未来的研究可在以下几个方向上进一步拓展：1）结合多源遥感数据，尤其是光学与雷达数据的互补融合，以提升多云多雨地区关键作物生长期的时序连续性与监测精度；2）在方法层面引入智能化算法，如基于机器学习与深度学习的峰值识别与曲线重构方法，更有效地区分作物真实生长信号与非作物干扰，从而提升复种指数的提取精度；3）在研究范围上，将本方法推广至与通海县具有相似种植特征的其他露天蔬菜主产区，以验证模型的普适性与稳定性。本研究并未分析不同年耕地复种指数变化的原因，后续将结合气候、土壤与市场等多维驱动因素，深入揭示复种指数的时空变化机理。这些探索将为耕地资源的优化利用与农业可持续发展提供更加可靠的科学支撑。
