CN111598019A - 基于多源遥感数据的作物类型与种植模式识别方法 - Google Patents

Abstract

本发明公开了一种基于多源遥感数据的作物类型与种植模式识别方法，包括如下步骤：1)输入单时相高空间分辨率遥感影像并进行数据预处理；2)输入多时相Sentinel‑2遥感影像并进行数据预处理；3)基于单时相影像进行面向对象多尺度分割，生成对象尺度矢量数据集；4)基于多时相影像进行数据质量优化，计算植被敏感指数，生成像元尺度植被指数时间序列数据集；5)构建对象尺度植被指数时间序列数据集；6)构建分类样本库并进行分层抽样；7)利用时间加权动态时间规整方法，对作物类型与种植模式进行识别。本发明以地块作为最小分类单元提取作物类型和种植模式，减少了因光谱异质性产生的提取误差，实现了作物类型的与种植模式的精准识别。

Description

基于多源遥感数据的作物类型与种植模式识别方法

技术领域

本发明涉及一种遥感农情监测方法，特别是指一种基于多源遥感数据的作物类型与种植模式识别方法。

背景技术

实时准确的作物监测对于辅助作物管理与适宜性评价、作物产量预估、作物灾害预警以及作物种植模式规划有重要意义。在区域尺度上进行精确作物制图可以为政府宏观调控提供新的数据支撑。传统的农作物实地调查监测的方式耗时耗力，满足不了大范围、快速及时的农情监测的需要。而随着遥感卫星技术的发展与智能化处理技术的提高，凭借遥感技术对地表信息获取的周期性、宏观性、时效性和经济性，以遥感技术为基础的农情监测成为未来农业发展与改革的重要方向。

当前，基于遥感技术的作物类型识别受到地面观测诸多因素的影响，如作物信息与影像像元所含植被信息不匹配、不同作物物候期具有相似性、相同作物存在耕作与生长周期的差异性、不同植被指数差异性与数据成本高等，精度和准确度有待提高。传统作物分类以像元为基本单位，利用农作物的光谱信息与植被物候信息，增强不同地物、作物的差异实现作物类型的分类。基于像元的作物识别因数据源的空间分辨率与技术等因素，产生的混合像元和分类结果“椒盐现象”往往限制了分类精度。

随着遥感影像空间分辨率的提高与面向对象技术的产生，面向对象高分辨率耕地信息提取的技术与方法得到广泛的探索与研究。尽管面向对象技术可以有效处理数据中“同质异谱”与“异质同谱”现象，减少椒盐噪声，提高分类精度，但由于影像时间分辨率低，无法详细刻画不同作物的生长物候信息。

如何在受到遥感数据时空分辨率的限制、数据质量、获取手段等因素的影响下，实现高精度、低成本、大范围的农业作物类型和种植模式遥感监测研究，是当前农业遥感研究和应用的前沿性和关键性问题。

发明内容

本发明的目的在于提供一种能够有效进行作物类型与种植模式监测的基于多源遥感数据的作物类型与种植模式识别方法。

为实现上述目的，本发明所提供的基于多源遥感数据的作物类型与种植模式识别方法，包括如下步骤：

1)输入单时相高空间分辨率遥感影像并对其进行预处理；

2)输入L1C级别的多时相时间序列Sentinel-2影像(大气表观反射率产品)并对其进行批量预处理，即利用欧空局Sen2cor插件对该影像数据集进行批量辐射校正至L2A级别(大气底层反射率产品)，并进行批量裁剪处理；

3)对于步骤1)中经过预处理的单时相高空间分辨率遥感影像，采用面向对象技术进行地块信息提取与分割，生成对象尺度矢量数据集；

4)对于步骤2)中经过批量预处理的Sentinel-2影像，先设置云含量阈值，利用云质量文件对影像进行掩膜处理；其次，对云含量大于阈值的区域进行数据重建；最终，批量计算植被敏感指数，生成像元尺度植被指数时间序列数据集；

5)结合步骤3)得到的对象尺度矢量数据集，对步骤4)得到的像元尺度植被指数时间序列数据集进行分区统计，将每个对象范围内的像元尺度植被指数时间序列数据集中的像元均值赋予给该对象范围内的所有像元，生成对象尺度植被指数时间序列数据集；

6)根据步骤5)得到的对象尺度植被指数时间序列数据集，结合野外采样各地物的样本点，构建分类样本库并进行分层抽样；

7)利用步骤6)中构建的分类样本库，结合野外采样作物的生长日历信息，利用时间加权动态时间规整的方法进行训练和预测，实现作物类型与种植模式的识别。

优选地，所述步骤1)中，预处理包括辐射校正、正射校正、几何校正、影像融合以及影像裁剪。

优选地，所述步骤3)具体包括如下步骤：

3.1)在高空间分辨率遥感影像中选择典型区域进行多层多尺度分割实验，所得矢量结果与野外采集地块数据进行叠加分析，确定面向对象多尺度分割最优尺度范围；叠加分析的主要原理是将调查结果与分割结果进行栅格化，统计调查结果内的像元数量落至分割结果内的个数占总像元个数的比例；

3.2)将步骤3.1)得到的最优尺度范围端点作为尺度评估的起始尺度与终止尺度，利用最优尺度评估工具进行最优尺度参数估计得到最优分割尺度；

3.3)根据步骤3.2)得到的最优分割尺度，基于eCognition Developer平台，利用多尺度分割算法对整个高空间分辨率遥感影像进行面向对象多尺度分割，生成对象尺度矢量数据集。

优选地，所述步骤3.2)中，利用最优尺度评估工具ESP2计算场景内的局部方差来进行最优尺度参数估计，局部方差变化率峰值所对应的分割尺度即为最优分割尺度；

优选地，所述步骤3.3)中，所述多尺度分割算法的参数包括对象异质性指标，所述对象异质性指标包括光谱异质性指标和形状异质性指标，其计算公式如下：

光谱异质性指标h_color：

式中：c为总波段数，W_c为波段权重值，δ_c为波段光谱值的标准差；

形状异质性指标h_shape：

h_shape＝W_smooth*h_smooth+W_compactness*h_compactness

式中：W_smooth为光滑度权重，W_compactness为紧致度权重，W_smooth+W_compactness＝1，h_smooth为光滑度异质性，h_compactness为紧致度异质性，E为分割后对象的周长，B为对象在水平方向上外接矩形的最短周长，N为对象内像元的总数；

对象异质性指标d：

d＝W_color*h_color+W_shape*h_shape

式中，W_color为光谱权重，W_shape为形状权重，W_color+W_shape＝1；h_color为光谱异质性指标，h_shape为形状异质性指标。

一般情况下，根据不同的研究区地形复杂度，地块规整度等特征的不同，各波段权重均设置为相等，光谱权重、形状权重、紧致度权重及光滑度权重的取值范围均为0～1之间，对象异质性指标的取值标准依据多次试验对比结果，取分割效果较为准确的参数组合。

优选地，所述步骤4)具体包括如下步骤：

4.1)设置云含量阈值为20％，利用L2A级别影像中的云质量文件，按照如下公式对影像进行掩膜处理：

式中，F(i,j)为影像像元值，NAN为空值，Cloud(i,j)为云质量文件中的像元值的百分比，i与j为对应像元的行列号；

4.2)基于MATLAB程序，对云含量大于20％，即像元值为NAN值的像元按照如下公式进行重建：

F_k(i,j)＝(F_k-1(i,j)+F_k+1(i,j))/2

式中，F_k(i,j)为第k景影像重建后的像元值，F_k-1(i,j)为第k-1景影像的像元值，F_k+1(i,j)为第k+1景影像的像元值，i与j为对应像元的行列号；当F_k-1(i,j)为NAN时，F_k-1(i,j)取时间轴左侧最邻近的像元值为非NAN时间点对应的影像像元值；当F_k+1(i,j)为NAN时，F_k+1(i,j)取时间轴右侧最邻近的像元值为非NAN时间点对应的影像像元值；

4.3)基于步骤4.2)重建后的影像，利用植被敏感计算公式，计算每个像元对应的植被指数值，获得由每个像元植被指数值构成的指数影像，得到像元尺度植被指数时间序列数据集；所述的植被指数包含NDVI与EVI，其计算公式如下：

NDVI＝(Nir-Red)/(Nir+Red)

EVI＝2.5*(Nir-Red)/(Nir+6*Red-7.5*Blue+1)

式中，Nir为近红波段地表反射率值，Red为红波段地表反射率值，Blue为蓝波段地表反射率值。

优选地，所述步骤5)中，按照如下公式将每一对象范围内的全部植被指数时间序列数据集中的像元的均值赋予到该对象范围内的全部像元，生成对象尺度植被指数时间序列数据集：

式中，H_g(i,j)为对象内第g个像元的植被指数值，q为对象内植被指数像元的个数，i与j分别代表对象范围内植被指数像元的行列号。

优选地，所述步骤6)中，结合野外采样各地物的样本点，构建分类样本库，先将总体样本按其类型进行分层，然后在每层中随机抽取样本，并将抽取的样本按一定比例分为训练样本与验证样本。其中训练样本与验证样本的抽取比例优选为2:1。

优选地，所述步骤7)具体包含以下步骤：

7.1)利用广义相加模型对各地物的植被指数物候曲线(如NDVI曲线、EVI曲线)进行拟合；

7.2)将步骤7.1)中训练样本的拟合曲线应用至时间加权动态时间规整算法中，进行作物类型与种植模式的识别；

7.3)利用步骤6)中的分层抽样后的验证样本，对步骤7.2)的分类结果进行精度验证；

7.4)如精度验证结果如果与期望精度不符，则返回步骤6)重新进行分层抽样，直至与期望精度相符，最终完成作物类型与种植模式的识别。

优选地，所述步骤7.3)中，精度验证的指标选择总体准确率(Overall accuracy)、精确率(Precision)、召回率(Recall)以及F-score四个指标。

与现有技术相比，本发明的有益效果在于：本发明基于单时相高空间分辨率数据与多时相时间序列Sentinel-2数据，突破了单一数据源限制，减少了农田因光谱异质性产生的提取误差，既满足了遥感作物监测所需要的生长规律信息，又实现了农田地块中作物类型的精准识别；同时结合时间序列分析与面向对象分割技术，以具有物候信息的时间序列对象为最小分类单元，利用加权动态时间规整的方法实现作物类型与种植模式的识别，解决了传统分别基于像元尺度与基于对象尺度进行相关监测存在的瓶颈，从时序对象角度降低了数据计算维度，优化了精准农业监测手段，提高了作物类型和种植模式识别精度。

附图说明

图1为本发明实施例中所提供的作物类型与种植模式识别的流程示意图。

图2为本发明实施例中输入数据成像时间示意图。

图3为本发明实施例中研究区A最优尺度范围对比图。

图4为本发明实施例中研究区B最优尺度范围对比图。

图5为本发明实施例中研究区A最优尺度计算(ESP2)示意图。

图6为本发明实施例中研究区B最优尺度计算(ESP2)示意图。

图7为本发明实施例中研究区A主要作物生长日历示意图。

图8为本发明实施例中研究区B主要作物生长日历示意图。

图9为本发明实施例中研究区A主要地类NDVI拟合曲线示意图。

图10为本发明实施例中研究区B主要地类NDVI拟合曲线示意图。

图11为本发明实施例中研究区A作物类型与种植模式分类结果示意图，原图为彩色，由于专利附图限制改成了灰度图片，原始颜色名称在图例中进行了标注。

图12为本发明实施例中研究区B作物类型与种植模式分类结果示意图，原图为彩色，由于专利附图限制改成了灰度图片，原始颜色名称在图例中进行了标注。

具体实施方式

下面结合附图和具体实施例对本发明作进一步的详细说明。

如图1所示，本发明所提供的一种基于多源遥感数据的作物类型与种植模式识别方法，包括如下步骤：

1)输入单时相高空间分辨率地表反射率遥感影像，对该数据进行辐射校正、正射校正、几何校正、影像融合以及影像裁剪工作；其中，单时相高空间分辨率遥感影像为空间分辨率1米及1米以内的四波段影像，并且获取时间要在年周期内。

实施例中，研究区A中所采用的单时相高空间分辨率数据来源于国产GF-2号数据，数据成像时间为2018年9月16日，经过预处理之后，生成数据为空间分辨率1米的四波段数据；研究区B中所采用的单时相高空间分辨率数据来源于国产GF-2号数据，数据成像时间为2017年9月16日，经过预处理之后，生成数据为空间分辨率1米的四波段数据。步骤1中涉及的预处理工作基于ENVI5.5软件完成。

2)输入时间序列L1C级别的Sentinel-2影像，利用欧空局Sen2cor插件对该数据集进行批量辐射校正得到L2A级别影像，并进行批量影像裁剪处理；其中，多时相时间序列影像是指在年周期内获取不同月份的遥感影像，时相范围应考虑研究区主要作物的生长物候时间。数据的起始时相点与终止时相点对应的数据的成像质量应该最佳，即不受到云、雨天气的影响。数据的空间分辨率为10米或10米以内并且由四波段组成。

实施例中，研究区A中所采用的时间序列Sentinel-2数据来源于欧洲航天局(https://scihub.copernicus.eu/)，考虑到该研究区11月初为油菜初生长阶段，因此成像时间为2017年11月09日至2018年11月09日，共19景数据，经过预处理之后，生成数据集为空间分辨率10米的多波段数据；研究区B中所采用的时间序列Sentinel-2数据来源于欧空局开源网站，考虑到该研究区11月初为油菜的初生长阶段，因此成像时间为2017年9月30日至2018年9月30日，共21景数据，经过预处理之后，改成数据集为空间分辨率为10米的多波段数据；步骤2中涉及的预处理工作是基于SNAP和ENVI5.5软件完成。有关步骤1)与步骤2)中涉及输入数据的成像时间见图2。

3)基于eCognition Developer8.9平台，利用高空间分辨率影像数据，结合野外采集地块信息与尺度评估工具进行面向对象多尺度分割，生成对象尺度矢量数据集，包括如下步骤：

3.1)选择比较有代表性的区域，利用野外采集地块信息，与该区域的多层多尺度分割后的矢量结果进行叠加分析，确定面向对象多尺度分割最优尺度范围。叠加分析的主要原理是将调查结果与分割结果进行栅格化，统计调查结果内的像元数量落至分割结果内的个数占总像元个数的比例。

3.2)取步骤3.1)中的最优尺度范围端点作为尺度评估的起始尺度与终止尺度，利用最优尺度评估工具ESP2计算场景内的局部方差来进行最优尺度参数估计，局部方差变化率峰值所对应的分割尺度即为最优尺度。

3.3)根据步骤3.2)得到的最优分割尺度，基于eCognition Developer 8.9平台，利用多尺度分割算法对整个高空间分辨率遥感影像进行面向对象多尺度分割，生成对象尺度矢量数据集。除分割尺度外，影响多尺度分割算法的因子还有对象异质性因子，后者用对象异质性指标来衡量，包括光谱异质性指标和形状异质性指标，而形状异质性指标又包括光滑度和紧致度，各指标的具体计算公式如下：

光谱异质性指标h_color：

式中：c为总波段数，W_c为波段权重值，δ_c为波段光谱值的标准差。

形状异质性指标h_shape：

h_shape＝W_smooth*h_smooth+W_compactness*h_compactness

式中：W_smooth为光滑度权重，W_compactness为紧致度权重，h_smooth为光滑度异质性，h_compactness为紧致度异质性，E为分割后对象的周长，B为对象在水平方向上外接矩形的最短周长，N为对象内像元的总数。其中：W_smooth+W_compactness＝1。

影像对象总异质性指标d：

d＝W_color*h_color+W_shape*h_shape

式中，W_color为光谱权重，W_shape为形状权重，h_color为光谱异质性，h_shape为形状异质性。其中W_color+W_shape＝1。

一般情况下，根据不同的研究区地形复杂度，地块规整度等特征的不同，各波段权重均设置为相等，光谱权重、形状权重、紧致度权重及光滑度权重的取值范围均为0～1之间，对象异质性指标的取值标准依据多次试验对比结果，取分割效果较为准确的参数组合。

实施例中，经过野外实地调查，利用地块大小、形状、范围信息确定研究区A与研究区B的地块最优尺度范围，参考标准依据调查结果内的像元数量落至分割结果内的个数占总像元个数的比例大于90％，则该分割尺度隶属于最优尺度范围。图3代表研究区A的最优尺度范围，图3(c)为欠分割状态，图3(h)为过分割状态，因此，研究区A的最优尺度范围为90～120。图4代表研究区B的最优尺度范围，图4(c)为欠分割状态，图4(h)为过分割状态，因此，研究区B的最优尺度范围为70～100。其次，利用最优尺度评估工具ESP2进行最优分割尺度计算，选择局部方差峰值所对应的分割尺度为最优分割尺度，两个实例的ESP2计算结果示意图见图5与图6，其中，图5的局部方差变化率为分割尺度100时最高，而图6的局部方差变化率为分割尺度80时最高，因此研究区A的最优分割度为100，研究区B的最优分割尺度为80。利用多尺度分割算法进行面向对象影像分割，除分割尺度外，研究区A与研究区B的光谱权重、形状权重、紧致度权重以及光滑度权重的设定依据以往研究经验与多次试验，均设置为0.5，各波段的波段权重均设置为1，最终生成对象尺度矢量数据集。

4)根据步骤2)中经过批量预处理的时间序列Sentinel-2影像，先设置阈值，利用云质量文件对影像进行掩膜处理；其次，对云含量大于阈值的区域进行数据重建；最终，批量计算植被敏感指数，生成像元尺度植被指数时间序列数据集。具体包括如下步骤：

4.1)设置云含量阈值为20％，利用L2A级别影像中的云质量文件，按照如下公式对时间序列Sentinel-2影像进行掩膜处理：

式中，F(i,j)为Sentinel-2影像像元值，NAN为空值，Cloud(i,j)为云质量文件中的像元值的百分比，i与j为对应像元的行列号。

4.2)基于MATLAB程序，对云含量大于20％，即像元值为NAN值的像元进行重建，重建的标准按照如下公式进行：

F_k(i,j)＝(F_k-1(i,j)+F_k+1(i,j))/2

式中，F_k(i,j)为第k景Sentinel-2影像重建后的像元值，F_k-1(i,j)为第k-1景Sentinel-2影像的像元值，F_k+1(i,j)为第k+1景Sentinel-2影像的像元值，i与j为对应像元的行列号。当F_k-1(i,j)为NAN时，F_k-1(i,j)取时间轴左侧最邻近的像元值为非NAN时间点对应的影像像元值；当F_k+1(i,j)为NAN时，F_k+1(i,j)取时间轴右侧最邻近的像元值为非NAN时间点对应的影像像元值。

4.3)基于步骤4.2)重建后的Sentinel-2影像，利用植被敏感计算公式，计算每个像元对应的植被指数值，获得由每个像元植被指数值组成的指数影像，所述的植被指数包含NDVI与EVI，公式如下：

NDVI＝(Nir-Red)/(Nir+Red)

EVI＝2.5*(Nir-Red)/(Nir+6*Red-7.5*Blue+1)式中，Nir为近红波段地表反射率值，Red为红波段地表反射率值，Blue为蓝波段地表反射率值。

实施例中，根据研究区A的19景Sentinel-2数据与研究区B中的21景Sentinel-2数据，以L2A文件中云质量文件为参考，20％的云含量为阈值，进行掩膜处理；利用MATLAB平台，对像元值为NAN值的像元进行重建的标准利用相邻时相此像元对应位置像元的均值进行重建；植被敏感指数的种类有很多，本实例仅以归一化植被指数NDVI为例，计算每个像元对应的NDVI值，获得由每个像元NDVI数值组成的指数影像，生成像元尺度植被指数时间序列数据集。

5)将步骤3)得到的对象尺度矢量数据集确定的对象范围，应用到步骤4)得到的像元尺度植被指数时间序列数据集，基于ArcGIS10.7软件的空间分析模块，进行分区统计，将某对象范围内的植被指数像元均值赋予该对象范围内的所有像元，生成对象尺度植被指数时间序列数据集，该步骤按照如下公式将每一对象范围内的全部植被指数时间序列数据集中的像元的均值赋予到该对象范围内的全部像元，生成对象尺度植被指数时间序列数据集：

式中，H_g(i,j)为对象内第g个像元的植被指数值，q为对象内植被指数像元的个数，i与j分别代表对象范围内植被指数像元的行列号。

6)结合野外采样各地物的样本点，构建分类样本库；综合考虑不同地物的物候曲线特征差异性较大，先将总体样本按其类型进行分层，然后在每层中随机抽取样本，并将抽取的样本按一定比例分为训练样本与验证样本。其中训练样本与验证样本的抽取比例优选为2:1。

实例中，研究区A的分类样本库的容量为300，共六类；研究区B的分类样本库的容量为300，共六类。对研究区A与研究区B的样本按照训练样本与验证样本2:1的比例进行分层(均为6层)抽取与训练，其中研究区A训练样本和验证样本的数量分别为200和100，研究区B训练样本和验证样本的数量分别为200与100。

7)结合野外采样主要作物的生长日历信息，利用时间加权动态时间规整的方法进行训练和预测，实现作物类型与种植模式的识别。主要包含以下步骤：

7.1)基于R语言平台，利用广义相加模型对各地物的植被指数物候曲线进行拟合，去除物候曲线不规则的样本，构建准确的分类样本库。广义相加模型不仅可以较好消除噪声、保持曲线形状与宽度，而且对影像时相特征较为不敏感，适合用于影像时相缺失等场景。

实施例中，结合研究区A与研究区B的主要作物类型物候的生长日历信息，利用广义相加模型的拟合方法，对研究区A的19景植被指数数据与研究区B的21景植被指数数据进行NDVI曲线拟合，构建准确的分类样本库，本步骤是基于R语言平台中dtwSat开源包完成。

实施例中，研究区A与研究区B的主要作物生长日历见图7与图8。其中，研究区A的主要作物包括单季中稻、油菜(一般与单季中稻轮作)、莲藕等，其中，单季中稻的时间为当年5月中下旬至9月下旬，生长周期约为130～140天；油菜的时间为前一年10月上中旬至当年5月中旬，生长周期约为250～270天；莲藕的时间为当年4月上旬至11月中下旬，持续时间较长；研究区B的主要作物包括：双季早稻、单季中稻、双季晚稻、油菜(一般与单季中稻轮作)等，其中，双季早稻时间为当年4月上中旬至7月中下旬，生长期约为90～110天左右；单季中稻时间为当年6月上旬至10月中上旬，生长期约为120～130天；双季晚稻的时间为当年7月下旬至10月中下旬，生长期约90～110天；油菜的时间为前一年10月中下旬至当年5月中下旬，生长期约为240天～250天。

实施例中，研究区A与研究区B各地类曲线结果见图9与图10。从拟合的曲线中可以看出，虽然研究区A和B的地物种类有差别，但是相同作物的物候曲线也十分类似，这是由于两个研究区均处于亚热带季风区，水、光、热等条件大致相同。其中，研究区A的主要代表作物是单季中稻(也称单季稻)与单季中稻-油菜(也称单季稻-油菜)，有轮作的单季中稻的物候曲线要比无轮作的物候曲线的时间点稍微后移，但总体相差不大，基本是五月底六月初开始生长，八月的NDVI达到峰值，九月底左右已完成收割。轮作的油菜属于过冬型作物，第一个峰值出现在一月中旬，待油菜生长到开花阶段，NDVI值就会呈现下降的趋势，三月中上旬达到花期的低谷，过了花期，油菜的NDVI值又会出现新的峰值并持续到五月中旬完成收割。此外，莲藕作为研究区A的特殊作物，每年的生长期较长，而种植莲藕的池塘水面常年覆盖浮游植物，因此全年的NDVI均值均在0.4及以上；研究区B的单季中稻与单季中稻—油菜轮作的物候时间与研究区A的类似，而研究区B独有的双季稻则呈现典型的双峰，双季稻中的早稻从四月中上旬一直持续至七月中旬，并在六月初左右出现第一个峰值，伴随着早稻的收割与晚稻的播种，晚稻的NDVI曲线会在九月初左右呈现第二个峰值，直至晚稻被收割。研究区A与研究区B的灌林地的NDVI曲线普遍较高，而不透水表层的NDVI曲线值很低但总体位于零值以上，水体的NDVI曲线则介于零值上下波动，波动范围不大。

7.2)基于R语言平台，将训练样本的拟合曲线应用至时间加权动态时间规整算法中，并实现作物类型与种植模式的识别。

实施例中，基于R语言平台的dtwSat开源包，需输入数据并调整相关参数即可。其中，时间加权动态时间规整(TWDTW)可通过logistic TWDTW和linear TWDTW两个模式来计算时间成本，由于logistic TWDTW在处理小时间规整方面优于linear TWDTW从而获取更高的分类精度，因此我们选用logistic TWDTW模式进行分类，选用参数α＝0.1与β＝50，其中，α代表梯度参数，β代表中点。α作为引入的时间权重参数，而β作为最大时间延迟，即曲线由于匹配而规整的最大时间为50天，这很大程度上避免了由于物候期微弱差异而导致的相同作物的曲线异质性。研究区A与研究区B的分类结果见图11与图12。

7.3)利用步骤6)中的验证样本，对步骤7.2)的分类结果进行精度验证，精度验证的指标选择总体准确率(Overall accuracy)、精确率(Precision)、召回率(Recall)以及F-score四个指标，其中计算公式如下：

accuracy＝(TP+TN)/(TP+FN+FP+TN)

precision＝TP/(TP+FP)

recall＝TP/(TP+FN)

F-score＝2*(precision*recall)/(precision+recall)

式中，TP为实际是正样本预测为正样本的样本数，FP为实际是负样本预测为正样本的样本数，TN为实际是负样本预测为负样本的样本数，FN为实际是正样本预测为负样本的样本数。

7.4)精度评估结果如果与期望精度不符，则返回步骤6)重新进行分层抽样，再重复步骤7即可，最终完成作物类型与种植模式精准农业制图与时空特征分析。

实施例中，在MATLAB程序中，对研究区A与研究区B的分类结果进行验证，研究区A中的总体准确率达到83.42％，研究区B中的总体准确率达80.32％，相关地类的召回率、精准率以及F-score指标见下表1、表2：

表1研究区A精度验证结果

表2研究区B精度验证结果

研究区A与研究区B的总体分类精度较好，准确率均大于80％，因此，不需要重复进行样本训练与分类；受地形地貌、光照、温度、劳动力成本及相关农业政策的影响，两个研究区的作物类型、种植模式及其分布存在差异。其中，研究区A以单季稻、油菜、莲藕等作物为主，一年一熟制的单季稻和一年两熟制的单季稻-油菜为主要种植模式，作物主要分布在河流、道路周边以及平坦区域；研究区B以单季稻、双季稻、油菜等作物为主，一年一熟制的单季稻和一年两熟制的单季稻-油菜、双季早晚稻为主要种植模式，其中单季稻-油菜及双季早晚稻主要分布在河流两岸、道路两边及延伸缓冲区，单季稻较为分散地分布在地块破碎程度高的低山丘陵地带。

Claims (10)

1.一种基于多源遥感数据的作物类型与种植模式识别方法，其特征在于：包括如下步骤：

1)输入单时相高空间分辨率遥感影像并对其进行预处理；

2)输入L1C级别的多时相时间序列Sentinel-2影像并对其进行批量预处理，即利用欧空局Sen2cor插件对该影像数据集进行批量辐射校正至L2A级别，并进行批量裁剪处理；

3)对于步骤1)中经过预处理的单时相高空间分辨率遥感影像，采用面向对象技术进行地块信息提取与分割，生成对象尺度矢量数据集；

4)对于步骤2)中经过批量预处理的Sentinel-2影像，先设置云含量阈值，利用云质量文件对影像进行掩膜处理；其次，对云含量大于阈值的区域进行数据重建；最终，批量计算植被敏感指数，生成像元尺度植被指数时间序列数据集；

5)结合步骤3)得到的对象尺度矢量数据集，对步骤4)得到的像元尺度植被指数时间序列数据集进行分区统计，将每个对象范围内的像元尺度植被指数时间序列数据集中的像元均值赋予给该对象范围内的所有像元，生成对象尺度植被指数时间序列数据集；

6)根据步骤5)得到的对象尺度植被指数时间序列数据集，结合野外采样各地物的样本点，构建分类样本库并进行分层抽样；

7)利用步骤6)中构建的分类样本库，结合野外采样作物的生长日历信息，利用时间加权动态时间规整的方法进行训练和预测，实现作物类型与种植模式的识别。

2.根据权利要求1所述的基于多源遥感数据的作物类型与种植模式识别方法，其特征在于：所述步骤1)中，预处理包括辐射校正、正射校正、几何校正、影像融合以及影像裁剪。

3.根据权利要求1所述的基于多源遥感数据的作物类型与种植模式识别方法，其特征在于：所述步骤3)具体包括如下步骤：

3.1)在高空间分辨率遥感影像中选择典型区域进行多层多尺度分割实验，所得矢量结果与野外采集地块数据进行叠加分析，确定面向对象多尺度分割最优尺度范围；

3.2)将步骤3.1)得到的最优尺度范围端点作为尺度评估的起始尺度与终止尺度，利用最优尺度评估工具进行最优尺度参数估计得到最优分割尺度；

3.3)根据步骤3.2)得到的最优分割尺度，基于eCognition Developer平台，利用多尺度分割算法对整个高空间分辨率遥感影像进行面向对象多尺度分割，生成对象尺度矢量数据集。

4.根据权利要求3所述的基于多源遥感数据的作物类型与种植模式识别方法，其特征在于：所述步骤3.2)中，利用最优尺度评估工具ESP2计算场景内的局部方差来进行最优尺度参数估计，局部方差变化率峰值所对应的分割尺度即为最优分割尺度。

5.根据权利要求3所述的基于多源遥感数据的作物类型与种植模式识别方法，其特征在于：

所述步骤3.3)中，所述多尺度分割算法的参数包括对象异质性指标，所述对象异质性指标包括光谱异质性指标和形状异质性指标，其计算公式如下：

光谱异质性指标h_color：

式中：c为总波段数，W_c为波段权重值，δ_c为波段光谱值的标准差；

形状异质性指标h_shape：

h_shape＝W_smooth*h_smooth+W_compactness*h_compactness

式中：W_smooth为光滑度权重，W_compactness为紧致度权重，W_smooth+W_compactness＝1，h_smooth为光滑度异质性，h_compactness为紧致度异质性，E为分割后对象的周长，B为对象在水平方向上外接矩形的最短周长，N为对象内像元的总数；

对象异质性指标d：

d＝W_color*h_color+W_shape*h_shape

式中，W_color为光谱权重，W_shape为形状权重，W_color+W_shape＝1；h_color为光谱异质性指标，h_shape为形状异质性指标。

6.根据权利要求1～5中任一项所述的基于多源遥感数据的作物类型与种植模式识别方法，其特征在于：

所述步骤4)具体包括如下步骤：

4.1)设置云含量阈值为20％，利用L2A级别影像中的云质量文件，按照如下公式对影像进行掩膜处理：

式中，F(i,j)为影像像元值，NAN为空值，Cloud(i,j)为云质量文件中的像元值的百分比，i与j为对应像元的行列号；

4.2)基于MATLAB程序，对云含量大于20％，即像元值为NAN值的像元按照如下公式进行重建：

F_k(i,j)＝(F_k-1(i,j)+F_k+1(i,j))/2

式中，F_k(i,j)为第k景影像重建后的像元值，F_k-1(i,j)为第k-1景影像的像元值，F_k+1(i,j)为第k+1景影像的像元值，i与j为对应像元的行列号；当F_k-1(i,j)为NAN时，F_k-1(i,j)取时间轴左侧最邻近的像元值为非NAN时间点对应的影像像元值；当F_k+1(i,j)为NAN时，F_k+1(i,j)取时间轴右侧最邻近的像元值为非NAN时间点对应的影像像元值；

4.3)基于步骤4.2)重建后的影像，利用植被敏感计算公式，计算每个像元对应的植被指数值，获得由每个像元植被指数值构成的指数影像，得到像元尺度植被指数时间序列数据集；所述的植被指数包含NDVI与EVI，其计算公式如下：

NDVI＝(Nir-Red)/(Nir+Red)

EVI＝2.5*(Nir-Red)/(Nir+6*Red-7.5*Blue+1)

式中，Nir为近红波段地表反射率值，Red为红波段地表反射率值，Blue为蓝波段地表反射率值。

7.根据权利要求1～5中任一项所述的基于多源遥感数据的作物类型与种植模式识别方法，其特征在于：

所述步骤5)中，按照如下公式将每一对象范围内的全部植被指数时间序列数据集中的像元的均值赋予到该对象范围内的全部像元，生成对象尺度植被指数时间序列数据集：

式中，H_g(i,j)为对象内第g个像元的植被指数值，q为对象内植被指数像元的个数，i与j分别代表对象范围内植被指数像元的行列号。

8.根据权利要求1～5中任一项所述的基于多源遥感数据的作物类型与种植模式识别方法，其特征在于：所述步骤6)中，结合野外采样各地物的样本点，构建分类样本库，先将总体样本按其类型进行分层，然后在每层中随机抽取样本，并将抽取的样本按一定比例分为训练样本与验证样本。

9.根据权利要求1～5中任一项所述的基于多源遥感数据的作物类型与种植模式识别方法，其特征在于：所述步骤7)具体包含以下步骤：

7.1)利用广义相加模型对各地物的植被指数物候曲线进行拟合；

7.2)将步骤7.1)中训练样本的拟合曲线应用至时间加权动态时间规整算法中，进行作物类型与种植模式的识别；

7.3)利用步骤6)中的分层抽样后的验证样本，对步骤7.2)的分类结果进行精度验证；

7.4)如精度验证结果如果与期望精度不符，则返回步骤6)重新进行分层抽样，直至与期望精度相符，最终完成作物类型与种植模式的识别。

10.根据权利要求9所述的基于多源遥感数据的作物类型与种植模式识别方法，其特征在于：所述步骤7.3)中，精度验证的指标选择总体准确率、精确率、召回率以及F-score四个指标。

Images