CN115512233A - 地块尺度耕地种植属性多时相遥感提取方法 - Google Patents

Abstract

本发明提供了一种地块尺度耕地种植属性多时相遥感提取方法，包括：第一步骤：执行数据收集及预处理；第二步骤：10m多时相NDVI时序数据构建；第三步骤：利用作物NDVI多时相变化特征，执行时序NDVI聚类分析；第四步骤：基于前期土地利用调查数据，结合卫星遥感数据聚类约束，生成地表覆盖类型时空标准样本集，并构造地物NDVI标准时序曲线；第五步骤：利用地物NDVI标准时序曲线，对10m多时相NDVI聚类结果进行类别标定，在前期土地利用图斑的基础上，对耕地地块进行统合约束，提取种植属性信息。

Description

地块尺度耕地种植属性多时相遥感提取方法

技术领域

本发明涉及耕地遥感监测技术领域，具体涉及一种地块尺度耕地种植属性多时相遥感提取方法。

背景技术

耕地是人类生存和发展的基本资源，随着社会经济发展和农业经营方式的变化，人类在耕地上实施的作物种植活动逐渐多样化，导致耕地中存在多种种植属性，如种植粮食、牧草、林木及温室农业等。这些利用方式不但减少了作物种植面积，而且影响了土壤肥力等耕地环境，制约农业规模化发展。第三次全国国土调查中明确要求对耕地现状进行详细调查，要求对耕地图斑内的利用状况进行清理及标注，如“临时种植园木”、“观赏园艺”、“速生林木”和“绿化草地”，归为其他地类；同时，按照实际耕作情况，逐图斑标注耕地种植属性，包括“未耕种”、“林粮间作”、“种植粮食作物”、“种植非粮作物”、“粮与非粮轮作”、“休耕”等。因此，在农业生产管理过程中迫切需要一种快速准确的获取耕地种植属性信息的方法，及时了解各种种植利用及变化情况，更加科学合理的制定农业发展和管理策略。

卫星遥感技术凭借其多尺度、多时相、多谱段的对地观测能力能够在空间和时间上提供准确和客观的地物信息，已经被广泛用于区域尺度的农田制图。Terra/MODIS逐日观测数据能够获取作物生长的动态变化信息，成为大尺度耕地利用监测中最常用的数据。与较低空间分辨率500m(可见光和近红外波段250米)的MODIS数据相比，Landsat系列卫星(15-30m)可以更准确的监测不同田块的种植差异，然而其重访周期较长，加之受云等气象条件影响，难以获取高时空分辨率的遥感数据。欧空局发射的具有5天重访周期10m分辨率的Sentinel-2卫星，为大区域田块尺度的时序遥感监测提供了数据基础。已有众多研究利用Sentinel-2卫星数据开展大区域、高精度的耕地作物监测，研究结果表明该卫星观测数据极大地提升了区域尺度作物监测的精度。Sentinel-2的10m空间分辨率观测数据一方面能够开展1；5万大比例尺土地利用制图、提升变化检测精度，另一方面也为拓展耕地质量的精准监测应用提供了契机，为面向土地管理需求的大区域田块尺度的耕地种植属性遥感监测提供了数据基础。

近年来，为了解决多源卫星观测数据收集和数据处理效率等问题，谷歌公司推出了谷歌地球云引擎(Google Earth Engine,GEE)遥感大数据平台。GEE存储了近40年来几乎所有公开的卫星遥感和地理空间数据集，数据总量达到BP级别；同时GEE提供了较为友好的人机交互操作方式。操作人员无需具备强大的计算机或专业级别的代码能力，就能够利用平台进行多源/多时相卫星观测数据的预处理，下载经过预处理后的数据，而无需下载大量的原始数据。很多基于GEE的区域甚至全球尺度的研究证明，GEE极大的提升了大区域、长时序、高空间分辨率遥感监测的能力，使得大区域，精细化的耕地种植属性监测成为可能。

另外一方面，为准确获取耕地种植属性信息，需要大量的训练和验证样本。目前以实地调查或目视解译方式难以获取大量的地面样本。研究表明训练样本的数量和质量对最终分类结果的影响最大，且制图精度随样本数量的增加而不断提高。而在区域尺度的制图过程中，需要耗费大量人力和时间制作训练样本集，并且难以对生成的样本进行大面积的更新。为提高样本获取的效率，众多国家基于已有土地利用调查产品进行样本筛选，如美国的农田数据层数据(CDL)，加拿大的全国作物底图数据，欧洲的LPIS数据等。中国第二次土地利用调查数据是地类最全、精度最高的1:1万土地利用制图产品。然而目前还没有利用该土地利用数据结合高时空分辨率遥感数据获取大量地物标准样本的应用案例。

发明内容

本发明所要解决的技术问题是针对现有技术中存在的上述缺陷，提供一种地块尺度耕地种植属性多时相遥感提取方法，能够利用已有土地利用数据结合高时空分辨率遥感数据获取大量标准样本，并在已有土地利用数据的基础上通过时序NDVI统合约束的方法提取耕地地块尺度的种植属性。根据本发明，提供了一种地块尺度耕地种植属性多时相遥感提取方法，包括：

第一步骤：执行数据收集及预处理；

第二步骤：基于收集的数据，执行10m多时相NDVI时序数据构建；

第三步骤：基于构建的10m多时相NDVI时序数据，利用作物NDVI多时相变化特征，执行时序NDVI聚类分析；

第四步骤：基于前期土地利用调查数据，结合聚类约束，生成地表覆盖类型标准样本集，并构造地物NDVI标准时序曲线；

第五步骤：利用地物NDVI标准时序曲线，对10m多时相NDVI聚类结果进行类别标定，在前期土地利用图斑的基础上，对耕地地块进行统合约束，提取种植属性信息。

优选地，第五步骤利用地物NDVI标准时序曲线，应用动态时间归整算法对聚类分析结果进行类别匹配，在前期土地利用数据基础上，针对耕地区域提取耕地种植属性信息。

优选地，在第五步骤中，针对待归类对象T，按照以下公式分别计算其与多个标准类型N的Wk值，最小W的标准类m即为待归类对象：

γ(i,j)＝∥Yi–Sj∥；

式中Yi为第i月待归类对象的NDVI(i＝1,2,…,12)；Sj为第j天标准NDVI(j＝1,2,…,365)；

Di＝min[γ(i,j)],j＝1,2,…,365；

式中Di为规整路径，i＝1,2,…12

Wk＝∑Di i＝1,2,…,12；

待归类对象与标准类型中最小距离的类型即为待归类对象的类型：

T∈m when Wm＝min(Wk)，k＝1,..,m,…,N。

优选地，收集的数据包括10m及以上空间分辨率的多时相卫星遥感数据以及已有土地利用调查数据。

优选地，收集的前期土地用数据精度需满足1；5万及以上的数据精度，例如第二次全国土地调查数据、第三次全国国土调查数据等。

优选地，10m多时相NDVI构建，需要采用但不限于月度最大值、均值、中值等合成方法，但数据的时间分辨率必须要能够覆盖作物的生长季且可以刻画作物的物候变化规律。

优选地，土地利用调查数据不包含耕地种植属性信息。

优选地，地表覆盖类型包括：林地、草地、裸地、水体、建设用地、设施农用地、园地、种植粮食作物1季耕地、种植粮食作物2季耕地、其他作物耕地。

优选地，多时相卫星遥感数据的收集、数据的预处理、时序NDVI的构建、时序NDVI的聚类分析、生成地物NDVI标准时序曲线可以但不限于在云平台中执行，其余处理在本地计算机中执行。

本申请提出了一种基于10m时序遥感数据的区域耕地种植属性信息遥感提取技术方法，该方法利用作物NDVI多时相变化特征，基于已有土地利用调查数据，结合卫星遥感数据聚类约束，选取典型类别标准样本，从而快速获取种植属性类型的地物样本集，并通过耕地地块尺度(图斑)的统合约束，快速提取耕地种植属性信息，总体精度可以达88％以上，解决了大区域耕地种植属性大量样本获取及耕地种植属性快速提取的问题，为我国耕地种植属性监测提供了一种全新的技术方法。

附图说明

结合附图，并通过参考下面的详细描述，将会更容易地对本发明有更完整的理解并且更容易地理解其伴随的优点和特征，其中：

图1示意性地示出了根据本发明优选实施例的地块尺度耕地种植属性多时相遥感提取方法的流程图。

图2示意性地示出了根据本发明优选实施例的地块尺度耕地种植属性多时相遥感提取方法的研究区位置及实地调查样本点示例的分布示意。

图3示意性地示出了2019年京津冀Sentinel-2A/B去云后有效观测次数。

图4示意性地示出了京津冀地物标准样本空间分布。

图5示意性地示出了京津冀北部区域标准样本2019年时序NDVI中值。

图6示意性地示出了京津冀南部区域标准样本2019年时序NDVI中值。

图7示意性地示出了根据本发明优选实施例的地块尺度耕地种植属性多时相遥感提取方法的耕地种植属性提取结果空间分布；其中a1和a2表示区域a的图像，b1和b2表示区域b的图像，c1和c2表示区域c的图像，d1和d2表示区域d的图像。

图8示意性地示出了根据本发明优选实施例的地块尺度耕地种植属性多时相遥感提取方法的耕地种植属性提取结果精度矩阵。

需要说明的是，附图用于说明本发明，而非限制本发明。注意，表示结构的附图可能并非按比例绘制。并且，附图中，相同或者类似的元件标有相同或者类似的标号。

具体实施方式

为了使本发明的内容更加清楚和易懂，下面结合具体实施例和附图对本发明的内容进行详细描述。

高精度的耕地种植属性制图是耕地数量及质量监测的基础。针对大区域耕地种植属性提取方法缺失的问题，本申请提出一种基于聚类及标准样本匹配的快速提取耕地种植属性方法(MSC-Basing on Matching Samples with Clustering Analysis)。首先，利用已有土地利用调查数据，结合卫星遥感数据聚类约束，选取典型类别标准样本，从而快速获取到包括种植粮食作物(1季)、种植粮食作物(2季)、草地、水域、裸地、大棚、林地、园地等类别的地物样本集，解决了大区域范围大量样本获取的问题。其次，利用样本NDVI(归一化植被指数)标准时序曲线，应用动态时间归整算法(Dynamic Time Warping，DTW)对聚类分析结果进行类别匹配，在耕地区域内提取耕地种植属性信息。测试结果表明，以前期土地利用数据为基础，利用时序遥感数据可以快速获取大量地物样本，快速提取耕地种植属性信息。本发明的方法可以实现精确、快速的大区域耕地种植属性信息提取，总体精度达0.91，满足区域耕地种植属性监测需求。

下面将结合附图来描述本发明的具体优选实施例。

图1示意性地示出了根据本发明优选实施例的地块尺度耕地种植属性多时相遥感提取方法的流程图。

如图1所示，根据本发明优选实施例的地块尺度耕地种植属性多时相遥感提取方法包括：

第一步骤S1：执行数据收集及预处理；

例如，收集的数据包括多时相卫星遥感数据以及已有土地利用调查数据。

在此，数据预处理的方法可以采用任何适当的方式进行。

第二步骤：基于收集的数据，执行10m多时相NDVI时序数据构建；

在此，多时相NDVI构建，需要采用但不限于月度最大值、均值、中值等合成方法，但数据的时间分辨率必须要能够覆盖作物的生长季且可以刻画作物的物候变化规律。

第三步骤：基于构建的10m多时相NDVI时序数据，执行时序NDVI聚类分析；

在此，聚类的类别数目一般设置为研究区地表覆盖类型数目的2-3倍。

第四步骤：基于前期土地利用调查数据，结合聚类约束，生成地表覆盖类型标准样本集，并构造地物NDVI标准时序曲线；

第五步骤：利用地物NDVI标准时序曲线，对10m多时相NDVI聚类结果进行类别标定，在前期土地利用图斑的基础上，对耕地地块进行统合约束，提取种植属性信息。

优选地，多时相卫星遥感数据的收集、数据的预处理、时序NDVI的构建、时序NDVI的聚类分析、生成地物NDVI标准时序曲线可以但不限于在云平台中执行，其余处理在本地计算机中执行

下面，以北京、天津、河北省(简称京津冀)为研究区具体描述本发明的实施方式。

利用GEE云平台，挖掘Sentinel-2A/B卫星多时相观测数优势，利用多时相NDVI数据结合已有土地利用图的约束，获取大量训练和验证本，提出利用10m多时相卫星观测实现区域尺度耕地种植属性快速监测方法，为区域的耕地精细管理提供决策依据。

京津冀地区(113°E-120°E，36°N-43°N)是中国北方经济总量最大，人口最密集的地区，总面积217.156km²。地势呈现西北高、东南低地的特点。区域南北气温、降水、风速等气象条件差异较大，北部以山区为主，呈中温带气候特点，南部以平原为主，呈暖温带气候特点。年雨量在400～800mm之间，降雨主要集中在夏季。该地区地利用变化剧烈、人地交互密集，存在多种耕地利用形式，包括：种植粮食作物，池塘养殖，温室农业等，是耕地种植属性研究的典型区域。

区域的热量资源决定着作物全年的生育时长，常用活动积温来表示。活动积温是指当某一地区的日平均气温达到10℃以上，且持续大于等于10℃这段时间的日平均气温之和。京津冀作物种植制度随3400℃年活动积温线而变化，北部为1年1季，主要种植玉米，莜麦，粟、高粱等粮食作物，一般为4月播种，10月收割。南部为1年2季，主要以冬小麦-夏玉米1年2季的轮作方式为主，每年10月前后播种冬小麦，次年3月到4月开始返青，植被覆盖度迅速增加。6月成熟收割后，继续种植玉米，于10月前后收割后部分地块休耕而另一部分将继续种植冬小麦。同时兼有部分一年一季的水稻、春玉米、藕等作物。均是在五月前后播种，9月前后达到生长高峰，10月收获。

数据收集及预处理

(1)多时相卫星遥感数据

基于Google Earth Engine(GEE)遥感大数据平台，收集了2019年覆盖京津冀范围的10m空间分辨率Sentinel-2A/B数据。数据全部为经过辐射校正的大气顶层反射率(L1C级)数据(TOA)，共计5856景。Sentinel-2A/B数据来自于欧洲航天局，是高分辨率，多光谱成像卫星，重返周期为5天，包括12个波段：其中4个10m可见光和NIR波段，6个为20m的红边和SWIR波段，2个60m的大气波段。在此使用了L1C级产品中10m可见光和NIR波段数据的QA60波段。

为了验证本发明提出的耕地种植属性遥感提取方法，研究收集了2019年覆盖研究区的200景高分1号及高分2号数据。高分1号数据拥有60km和800km两种幅宽，拥有2m分辨率全色/8m分辨率多光谱和16m分辨率多光谱的空间分辨率，重返周期41天(不侧摆时)。高分2号数据，幅宽45km，拥有多光谱4m和全色波段1m的空间分辨率，重返周期69天(不侧摆时)。本发明利用HSI变换融合方法，生成了2m分辨率高分1号多光谱数据和1m分辨率高分2号多光谱数据。详细数据情况见下表。

(2)土地利用及地面调查数据

收集了2017年京津冀土地利用调查数据，来自于第二次全国土地利用调查，包括耕地，园地，林地，草地，水域等8个1级类和38个2级类，制图比例尺达1:10000。国家每年要开展一次土地利用变更调查工作，投入大量的人力和物力，通过高分辨率遥感影像(高分系列(GF)、Worldview、Quickbird及航空影像等)的目视判读和地面核查，在前一年土地利用数据基础上，更新土地利用数据。在当前的具体示例中，数据来自中国国土勘测规划院，土地利用数据处理和分析结果(文中所示结果图)均为无坐标信息的图片以及文字，以JPEG及txt格式保存，所有考出成果均经过严格审查。

通过地面调查进一步验证耕地种植属性提取结果，于2018年和2019年对研究区进行了地面调查和验证。在调查期间，记录了各种耕地种植属性并拍摄现场照片，包括种植粮食作物(1季)、种植粮食作物(2季)、经济作物，草地，温室农业等，共246个典型样点。

具体处理操作

耕地内不同植被的生长、发育、收割或凋谢等生长过程均能由多时相NDVI数据刻画揭示。本发明提出直接利用多时相NDVI时序变化信息提取耕地种植属性的方法，提高信息提取效率。利用地物NDVI标准时序曲线，应用动态时间归整算法(Dynamic TimeWarping，DTW)对聚类分析结果进行类别匹配，针对耕地区域提取耕地种植属性信息，提出的该方法称为MSC(basing on Matching Samples with Clustering analysis)方法。主要流程包括时序NDVI构建与聚类分析、地物标准样本的建立、种植属性提取、方法的对比与验证，方法流程如图1所示。

时序NDVI构建与聚类分析

统计分析研究区Sentinel-2A/B多时相云掩膜后数据产品覆盖频率后发现，虽然Sentinel-2A/B观测周期是5天，但由于近邻轨道观测的叠加，一些区域可以2天重复一次观测，总体上至少可以每月获取到无云覆盖的有效数据(图3)。由此，利用GEE平台，收集研究区2019年Sentinel-2A/B数据，计算各景影像的归一化植被指数(NDVI＝(NIR-Red)/(NIR+Red))，以月为单位进行NDVI最大值合成处理，然后从平台下载处理后的多时相NDVI数据。

在多时相NDVI变化中各地表覆盖类型具有不同的变化特征，利用该特征可以识别不同的地表覆盖。利用K-mean聚类方法对2019年12个月的NDVI时序数据进行聚类分析。K-mean聚类的分类数目一般设置为研究区地表覆盖类型数目的2-3倍，处理中设置初始聚类20类，迭代次数40次，变换阈值为0.05％。为解决区域范围太大导致聚类运算困难问题，本研究按照京津冀行政区划(13个区域)分区进行了聚类处理。

标准样本与NDVI时序数据的生成

准确获取耕地种植属性信息，需要大量的训练和验证样本。本发明利用京津冀地区2017年土地利用图及2019年时序NDVI数据，快速提取种植属性样本，以提高耕地种植属性信息提取效率。处理过程中的关键技术方法包括时序NDVI构建与聚类分析、图斑空间叠加提取纯类别等。土地利用调查数据是经过从卫星影像目视判读、地面调查等多层验证的高精度数据。

这里所利用的土地利用调查数据并没有包含耕地种植属性信息如粮食作物种植季数、温室等地表实时覆盖信息；且由于是前一年数据，当年其土地利用也会发生变化。为此，本发明叠加当年多时相NDVI聚类分析结果，提取了地表覆盖类型标准样本，地表覆盖类型包括10类：林地(包含森林和灌木林)、草地、裸地、水体、建设用地、设施农用地(蔬菜大棚及畜舍等)、园地、种植粮食作物(1季)耕地、种植粮食作物(2季)耕地、其他作物耕地(经济作物、蔬菜等)。进一步利用GEE，针对各类样本(polygon)，以Sentinel-2A/B数据像元为单位提取像元在全年中不同观测时间获取的无云NDVI值，并进行时序平滑；例如以1天为单位计算样本中所有平滑后像元NDVI的中值，得到各类标准NDVI时序值。

种植属性信息提取与验证

利用地物NDVI标准时序曲线，应用动态时间归整算法(Dynamic Time Warping，DTW)对聚类分析结果进行类别匹配，针对耕地区域提取耕地种植属性信息。DTW算法可以计算多类排列组合的2条时间序列曲线之间的距离，以最短距离组成的曲线称为归整路径，规整路径的和称为DTW距离(W)。针对待归类对象T，按照以下公式分别计算其与多个标准类型(N＝10)的Wk值，最小W的标准类m即为待归类对象。

γ(i,j)＝∥Yi–Sj∥…(1)

式中Yi为第i月待归类对象的NDVI(i＝1,2,…,12)；Sj为第j天标准NDVI(j＝1,2,…,365)

Di＝min[γ(i,j)],j＝1,2,…,365…(2)

式中Di为规整路径，i＝1,2,…12

Wk＝∑Di i＝1,2,…,12…(3)

式中Di规整路径的和，k＝1,..,N

待归类对象与标准类型中最小距离的类型(公式(4))即为待归类对象的类型。

T∈m when Wm＝min(Wk)，k＝1,..,m,…,N…(4)

应用精度混淆矩阵计算出各类指标对比验证了两种方法的精度，包括：总体精度、制图精度、用户精度及Kappa系数等。

地物标准样本集时空分布

应用上述土地利用调查数据和NDVI聚类分析的组合方法，共计筛选提取出47120个标准样本(见图4)，其中林地10217个、草地4434个、裸地394个、建设用地4098个、水体1230个、园地2381个、设施农用地4617个、种植粮食作物(1季)耕地8616个、种植粮食作物(2季)耕地7756个、其他作物耕地3377个。

根据样本的空间位置，利用GEE平台统计每类样本2019年时序NDVI中值，获取地表类型NDVI时序样本数据(图5)。京津冀研究区范围面积大，南北气温的差异导致作物种植和收割时间等植被物候不同，由此显示出不同的NDVI时序变化特征。为此在提取NDVI时序样本时，参考10℃以上活动积温气象数据，以作物种植的积温需求节点3400℃活动积温为界限，将京津冀区域分割南北两个区，针对10类在两个区分别建立了2019年NDVI时序数据标准样本(图5)。图5显示出各地表覆盖类型的NDVI时序变化特征，各类型间的NDVI时序变化体现了明显差异，其中耕地中作物的NDVI时序变化与其他地物类型的差异最为显著，由此可以利用该NDVI时序变化作为地表覆盖分类的参照标准。

耕地种植属性提取结果

京津冀区域耕地面积共75407km²，2019年京津冀地区耕地实际种植粮食作物面积49840km²，其中种植粮食作物(1季)和种植粮食作物(2季)种植面积占比分别为36％和29％，其他作物种植和温室大棚等设施农用地分别占13％和4％，而非作物种植的利用占据了18％，其中园林木、休耕地、建设开发用地和坑塘水体分别占7％、7％、2.7％、1.3％。在非作物种植的利用中，休耕地以张家口面积最大，占京津冀全区休耕地的47％，这与该区域半干旱脆弱的地理环境有关，为保持土壤肥力，当地农民多以休耕轮作方式利用耕地。

从空间上，粮食作物(2季)主要分布在保定、衡水、邯郸、邢台、沧州等南部区域，其余区域以粮食作物(1季)为主，符合3400℃积温线的地理格局。园地占用耕地的情况主要发生在北京及石家庄的东南地区，存在耕地“非粮化”的问题。此外，为保护土壤耕作层，维持土壤肥力，近年来耕地的休耕轮作比例逐渐增加，特别是在研究区北部张家口地区最为明显，休耕地主要分布在这一区域，这些区域应是耕地“非农化”及“非粮化”监测的重点。

图7显示了不同市域各类种植属性面积。除北京外，各市域的耕地均以种植1季或2季的种植方式为主，面积均为最大，而北京则是园地面积最大。在各市域中，邢台的2季作物种植面积最多，超过3000km²；张家口的1季作物种植面积最多，超过5000km²；保定的其他作物种植面积最多，约为1600km²。

精度验证

为了更好的验证本发明所提耕地种植属性提取方法，本发明分别使用246个实地调查样点集和1200个实地样点结合遥感目视判读样点集对耕地种植属性结果进行验证。精度矩阵显示实地调查样点集的验证精度更高，总体精度为91％，kappa系数为0.87。实地样点结合遥感目视判读样点集的验证精度略低(图7)。验证结果表明，本发明提出的MSC方法可以对粮食作物(1季)、粮食作物(2季)、水域、园地、建设用地进行很好的划分，用户精度均大于90％。但是，对于休耕地、设施农用地和其他作物的分类结果精度略低。在用户精度方面，园地和种植粮食作物(2季)的精度最高。而生产者精度方面，水体和建设用地的精度最高。

高精度的耕地种植属性制图是耕地数量及质量监测的基础。卫星遥感技术具备多尺度、多时相、多谱段的对地观测能力，能够在空间和时间上提供准确和客观的地物信息，为区域尺度的高精度耕地种植属性信息提取提供有效途径。本发明采用时序遥感数据，集成非监督聚类、空间叠置分析、时间动态规整等方法提取耕地种植属性信息。在已有土地利用数据的基础上，利用地物的时序NDVI变化特征，快速提取耕地种植属性信息，为实际耕地种植属性监测提供了全新的提取策略。通过地面样点和地面结合遥感判断的2种验证方式，总体精度均大于88％，提取精度较高。

Sentinel-2A/B卫星影像具备较高的时空分辨率，对比以往应用较多的LandsatTM/OLI、MODIS时间序列数据，Sentinel-2A/B数据具有回归周期短，数据时空分辨率高等优势，在大区尺度的农作物精细识别、种植属性提取等方面具备巨大潜力。相较于单一时相的地表覆盖分类，月度融合的NDVI可以更好的刻画作物生长特征，可以大大提高作物的分类精度。通过结合现有土地利用数据，可以解决以往非监督分类的类别划分及类别标注的问题，使分类结果更加精确。

本发明所提方法仍有小范围的“错分”和“漏分”现象，裸地、草地、建筑存在部分错分、误分现象。通过对分类后的结果进行人工复检、修改，可以解决道路、草地和田间附属林带误分成农作物的问题。其次，后续的研究还可以对Sentinel-2A/B影像中包含的空间信息进行深度挖掘，例如地物的纹理特征、形状特征等，来提升作物的分类精度。

总之，本申请提出了一种基于10m时序遥感数据的区域耕地种植属性信息遥感提取技术方法，该方法利用作物NDVI多时相变化特征，基于已有土地利用调查数据，结合卫星遥感数据聚类约束，选取典型类别标准样本，从而快速获取种植属性类型的地物样本集，并通过耕地地块尺度(图斑)的统合约束，快速提取耕地种植属性信息，总体精度可以达88％以上，解决了大区域耕地种植属性大量样本获取及耕地种植属性快速提取的问题，为我国耕地种植属性监测提供了一种全新的技术方法。

需要说明的是，除非特别指出，否则说明书中的术语“第一”、“第二”、“第三”等描述仅仅用于区分说明书中的各个组件、元素、步骤等，而不是用于表示各个组件、元素、步骤之间的逻辑关系或者顺序关系等。

可以理解的是，虽然本发明已以较佳实施例披露如上，然而上述实施例并非用以限定本发明。对于任何熟悉本领域的技术人员而言，在不脱离本发明技术方案范围情况下，都可利用上述揭示的技术内容对本发明技术方案作出许多可能的变动和修饰，或修改为等同变化的等效实施例。因此，凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所做的任何简单修改、等同变化及修饰，均仍属于本发明技术方案保护的范围内。

Claims (10)

1.一种地块尺度耕地种植属性多时相遥感提取方法，其特征在于包括：

第一步骤：执行数据收集及预处理；

第二步骤：基于收集的数据，执行10m多时相NDVI时序数据构建；

第三步骤：基于构建的10m多时相NDVI时序数据，利用作物NDVI多时相变化特征，执行10m多时相NDVI聚类分析；

第四步骤：基于前期土地利用调查数据，结合卫星遥感数据聚类约束，生成地表覆盖类型时空标准样本集，并构造地物NDVI标准时序曲线；

第五步骤：利用地物NDVI标准时序曲线，对10m多时相NDVI聚类结果进行类别标定，在前期土地利用图斑的基础上，对耕地地块进行统合约束，提取耕地种植属性信息。

2.根据权利要求1所述的地块尺度耕地种植属性多时相遥感提取方法，其特征在于，第五步骤利用地物NDVI标准时序曲线，应用动态时间归整算法对聚类分析结果进行类别匹配，在前期土地利用数据基础上，针对耕地区域提取耕地种植属性信息。

3.根据权利要求1或2所述的地块尺度耕地种植属性多时相遥感提取方法，其特征在于，在第五步骤中，针对待归类对象T，按照以下公式分别计算其与多个标准类型N的Wk值，最小W的标准类m即为待归类对象：

γ(i,j)＝∥Yi–Sj∥；

式中Yi为第i月待归类对象的NDVI(i＝1,2,…,12)；Sj为第j天标准NDVI(j＝1,2,…,365)；

Di＝min[γ(i,j)],j＝1,2,…,365；

式中Di为规整路径，i＝1,2,…12

Wk＝∑Di i＝1,2,…,12；

待归类对象与标准类型中最小距离的类型即为待归类对象的类型：

T∈m when Wm＝min(Wk)，k＝1,..,m,…,N。

4.根据权利要求1或2所述的地块尺度耕地种植属性多时相遥感提取方法，其特征在于，收集的数据包括10m及以上空间分辨率多时相卫星遥感数据以及已有土地利用调查数据。

5.根据权利要求1或2所述的地块尺度耕地种植属性多时相遥感提取方法，其特征在于，收集的前期土地用数据精度需满足1；5万及以上的数据精度，例如第二次全国土地调查、第三次全国国土调查数据等。

6.根据权利要求1或2所述的地块尺度耕地种植属性多时相遥感提取方法，其特征在于，针对多时相NDVI构建，需要采用但不限于月度最大值、均值、中值等合成方法，但数据的时间分辨率必须要能够覆盖作物的生长季且可以刻画作物生长的物候变化规律。

7.根据权利要求1或2所述的地块尺度耕地种植属性多时相遥感提取方法，其特征在于，利用已有土地利用调查数据结合10m时序NDVI聚类结果，快速提取大量地物样本。

8.根据权利要求4所述的地块尺度耕地种植属性多时相遥感提取方法，其特征在于，土地利用调查数据不包含耕地种植属性信息。

9.根据权利要求1或2所述的地块尺度耕地种植属性多时相遥感提取方法，其特征在于，地表覆盖类型包括：林地、草地、裸地、水体、建设用地、设施农用地、园地、种植粮食作物1季耕地、种植粮食作物2季耕地、其他作物耕地。

10.根据权利要求1或2所述的地块尺度耕地种植属性多时相遥感提取方法，其特征在于，多时相卫星遥感数据的收集、数据的预处理、时序NDVI的构建、时序NDVI的聚类分析、生成地物NDVI标准时序曲线可以但不限于在云平台中执行，其余处理在本地计算机中执行。

Images