本发明涉及一种遥感农情监测方法，特别是指一种基于图像锐化的时序数据耕地提取方法。
背景技术：
以土地资源为生产对象的农业作为国家的第一产业，支撑着国民经济的建设与发展，也是联合国可持续发展的重要组成部分。实时准确地监测耕地也是21世纪粮食安全的关键所在。耕地资源要素空间分布是耕地管理、精准农业以及土地确权的基本要求。然而，传统的耕地资源监测主要是以实地调查监测的方式为主，这种方式耗时耗力，满足不了大范围、快速准确监测的要求。随着遥感技术的发展和各种智能化技术的提高，越来越多的生产部门将遥感作为主要是监测方式进行耕地资源监测，其周期性、时效性以及经济性的特点受到很大的利用。因此，以遥感技术为基础的耕地资源监测是今后国情监测、土地权属分配等领域的重要支撑技术。
当前，基于遥感技术的耕地资源监测受到诸多因素的影响：1)复杂自然环境中，耕地要素在区域之间的差异明显，特别是不同区域的地形地貌的不同，导致耕地资源要素在遥感影像上呈现完全不一致的特点。例如：东北地区主要以平原为主，地块规整；而南方低山丘陵区地势连绵起伏，小农特征明显，耕地破碎化程度严重，边界模糊；2)耕地具有典型的季节性特征，相同区域内由于作物类型与种植模式的差异，导致区域内耕地类型异质性较大，复杂的物候特征增加了其监测的难度；3)高时空遥感数据可用性不足。尽管越来越多的数据源能够为农业监测带来契机，但是遥感数据却依然无法在时间和空间分辨率上达到一种平衡
多源遥感影像时空融合虽然可以在一定程度上解决遥感数据在时间分辨率和空间分辨率上难以统一的问题，但是多源数据间跨尺度和跨模态匹配问题仍是一大挑战，并且时空融合结果大多是构建中分辨率数据集，而高空间分辨率时序数据集仍然很匮乏。随着各种智能化技术的提高，以机器学习为主的人工智能方法逐渐适应了遥感领域，此类方法模型具有良好的鲁棒性，并且能产生相对较好的监测结果，但是依然受到农业监测最主要的一个问题：真实样本数据缺乏，特别是季节性描述样本的不足。因此，针对耕地监测以及当今农业监测领域所面临的耕地监测难、监测粒度不精确、数据的高时空特性不足等问题，通过融合多源不同尺度多源遥感数据，构建高空间率时序数据集，基于机器学习的方法实现高精度、低成本、大范围的耕地资源要素遥感监测研究，是当前农业遥感研究和应用重要一项内容。
技术实现要素：
本发明的目的在于提供一种能够有效进行耕地资源监测的基于图像锐化的时序数据耕地要素识别方法。
为实现上述目的，本发明所提供的一种基于图像锐化的时序数据耕地提取方法，包括如下步骤：
1)输入单时相高空间分辨率遥感影像全色波段并进行数据预处理；预处理步骤包括正射校正、辐射校正、几何配准、影像裁剪等。在几何配准时，需要手工选择若干控制点与辅助数据(例如谷歌数据等)进行几何配准。
2)输入l2a级别的多时相时间序列sentinel-2影像并对其进行批量预处理，即通过欧空局snap平台进行格式转换输出与批量裁剪处理；
3)对于步骤1)和步骤2)的结果，基于单时相高空间分辨率全色波段和多时相sentinel-2多光谱波段进行gram-schmidt图像锐化操作，生成高空间分辨率时序数据；
4)对于步骤3)中的高空间分辨率时序数据，计算归一化差异植被指数和土壤调节植被指数，生成高空间分辨率植被指数时间序列数据集；并通过对比锐化前与锐化后的时序数据集进行植被指数相关性分析。
5)根据步骤4)得到的高空间分辨率植被指数时间序列数据集，结合野外采样各地物的样本点，构建分类样本库并进行分层抽样；
6)利用步骤5)中构建的分类样本库，先利用随机森林算法进行模型训练，通过测试集的表现确定最优模型参数；将训练好的模型应用至全区域，实现全区域作物类型的识别，并对比原始时间序列数据提取结果与锐化后的时间序列提取结果；
7)利用步骤6)得到的作物类型的结果，将提取的所有类型的作物进行整合，形成耕地要素数据，将所有的非作物类型进行整合，形成非耕地数据，并进行精度验证；最终生成耕地资源要素空间分布数据。
优选地，所述步骤1)中，预处理包括辐射校正、正射校正、几何校正、影像融合以及影像裁剪。
优选地，所述步骤1)中，高空间分辨率全色波段影像预处理包括辐射校正、正射校正、几何配准以及影像裁剪。在进行几何配准时，需要依托于高精度谷歌数据、天地图等开源数据，手动选择若干控制点进行几何纠正。
优选地，所述步骤2)中，输入时间序列l2a级别的sentinel-2影像，并进行批量影像裁剪处理；其中，输入的sentinel-2的影像必须经过严格的云量控制，例如：通过整景云量小于5％或者研究区云量小于2％来进行数据的筛选；所利用到的时序sentinel-2是由蓝、绿、红以及近红外四个波段组成，空间分辨率为10米。
优选地，所述步骤3)中具体包括以下步骤：
3.1)使用多波段低分辨率影像对单波段高分辨率影像进行模拟。模拟的方法根据光谱响应函数wi进行模拟，即模拟的全色波段影像灰度值
3.2)利用模拟的高分辨率波段影像作为gram-schmidt变换第一分量来对模拟的高分辨波段影像和低分辨率波段影像进行gram-schmidt变换。在进行gram-schmidt变换时，第t个变换分量由前t-1个变换分量构造，即：
其中，gst是gs变换后产生的第t个正交分量，bt是原始低空间分辨率多光谱影像的第t个波段影像，ut是第t个原始低空间分辨率多光谱波段影像灰度值的均值，i,j分别表示低空间分辨率影像的行数和列数，
其中，c是影像的列数，r是影像行数，协方差
其中，标准差σ的计算公式如下：
3.3)根据gs第一分量，即模拟波段的均值u和标准差σ对全色波段进行修改，修改公式如下：
p′＝(p×k1) k2
k2＝uintensity-(k1×upan)
其中：p为原始全色波段影像的灰度值，eintensity和epan分别为gs第一分量和全色波段影像灰度值的均值。
3.4)最后将修改后的全色波段作为第一分量，进行gs逆变换，输出n 1个波段，去除第一个波段，就是融合后的高分辨率多光谱影像结果，逆变换的计算公式如下：
优选地，所述步骤4)具体包括如下步骤：
4.1)基于高空间时序数据进行归一化差异植被指数和土壤调节植被指数计算，得到植被指数高空间分辨率时序数据集；所述的两种植被指数的计算公式如下：
归一化差异植被指数ndvi：
ndvi＝(nir-red)/(nir red)
土壤调节植被指数savi：
savi＝(nir-red)*(1 l)/(nir red l)
式中，nir为近红外波段地表反射率值，red为红波段地表反射率值，l是随着植被密度变化的参数，取值范围是从0～1，当植被覆盖度很高时为0，很低时为1，通常取0.5时savi消除土壤反射率的效果最好。
4.2)将原始的时间序列数据进行归一化差异植被指数和土壤调节植被指数的计算，得到植被指数原始时序数据集；
4.3)随机选取区域内1～5万个点，对锐化前与锐化后不同植被指数进行相关性分析；随机选点是通过arcgis的随机产生点的工具进行，并通过多值提取至点，进行不同时相的点之间的相关性大小；
优选地，所述步骤5)，主要是通过野外采集的各土地利用类型(主要是作物类型)的样本点，建立该区域作物类型样本库，按照一定的比例进行分层抽样，获取训练数据集与验证数据集。
优选地，所述步骤6)具体包括如下步骤：
6.1)将训练数据集进行划分，抽取10％作为训练模型的测试数据，根据测试数据集的准确率，选择合适的训练参数，例如：随机树的数量，叶子节点的数量以及最大叶深度等；
6.2)将训练好的模型应用至整个研究区，完成作物类型分类，包括锐化后的高空间时序数据集和锐化前的时序数据集等。
优选地，所述步骤7)具体包括如下步骤：
7.1)将不同的作物类型进行整合，形成耕地数据；将非作物类型整合，形成非耕地数据。
7.2)将验证集中的不同作物类型的样本点合并成一类，即耕地样本点，并进行精度评价。
7.3)如精度验证结果如果与期望精度不符，则返回步骤5)重新进行分层抽样并重复步骤6)与步骤7)，直至与期望精度相符，最终完成耕地类型识别。
与现有的技术相比，本发明的有益效果在于：本发明基于单时相高空间分辨率全色波段影像和时间序列sentinel-2数据，突破了单一数据源在空间分辨率或时间分辨率上的制约，大大提高了耕地监测的颗粒度，经过图像锐化产生的高时空时序数据集既满足了遥感耕地资源要素监测所需要的季节性信息，又实现了亚米级别的要素类型提取；同时结合鲁棒性好，泛化能力强的随机森林算法，解决了耕地要素监测所面临的高时空数据不足、提取精度较低以及监测颗粒度较粗等瓶颈问题。图像锐化时序数据集的构建从高时空数据融合的角度解决了遥感影像在空间和时间分辨率不平衡的问题，在保证精度的情况下，实现作物类型与耕地要素的精细化识别。
附图说明
图1为本发明所提供的一种基于图像锐化的时序数据耕地提取方法的流程示意图。
具体实施方式
下面结合附图对本发明作进一步的详细说明。
如图1所示，本发明所提供的一种基于图像锐化的时序数据耕地提取方法，包括如下步骤：
1)输入单时相高空间分辨率全色波段地表反射率遥感影像，对该数据进行辐射校正、正射校正、几何配准以及影像裁剪工作；其中，单时相高空间分辨率全色波段遥感影像为空间分辨率1米及1米以内的四波段影像，并且获取时间要在年周期内。辐射校正通过不同传感器的定标系数来完成；正射校正利用全球30米数字高程模型数据辅助完成；几何配准需要基于谷歌地球、天地图等高分辨率无偏移辅助数据，手动选择地面控制点完成。
2)输入时间序列l2a级别的sentinel-2影像，并进行批量影像裁剪处理；其中，多时相时间序列影像是指在年周期内获取不同月份的遥感影像，时相范围应考虑研究区不同作物类型的生长物候时间。其中，输入的sentinel-2的影像必须经过严格的云量控制，例如：通过整景云量小于5％或者研究区云量小于2％来进行数据的筛选；所利用到的时序sentinel-2是由蓝、绿、红以及近红外四个波段组成，空间分辨率为10米；
3)基于单时相高空间分辨率全色波段和多时相sentinel-2多光谱波段进行gram-schmidt图像锐化操作，生成高空间分辨率时序数据，gram-schmidt图像锐化改进主成分变换算法中信息过分集中的问题，不受波段限制，较好的保持空间纹理信息，尤其能高保真保持光谱特征。专为最新高空间分辨率影像设计，能较好保持影像的纹理和光谱信息。gram-schmidt图像锐化包括如下步骤：
3.1)使用多波段低分辨率影像对单波段高分辨率影像进行模拟。模拟的方法根据光谱响应函数wi进行模拟，即模拟的全色波段影像灰度值
3.2)利用模拟的高分辨率波段影像作为gram-schmidt变换第一分量来对模拟的高分辨波段影像和低分辨率波段影像进行gram-schmidt变换。在进行gram-schmidt变换时，第t个变换分量由前t-1个变换分量构造，即：
其中，gst是gs变换后产生的第t个正交分量，bt是原始低空间分辨率多光谱影像的第t个波段影像，ut是第t个原始低空间分辨率多光谱波段影像灰度值的均值，i,j分别表示低空间分辨率影像的行数和列数，
其中，c是影像的列数，r是影像行数，协方差
其中，标准差σ的计算公式如下：
3.3)根据gs第一分量，即模拟波段的均值u和标准差σ对全色波段进行修改，修改公式如下：
p′＝(p×k1) k2
k2＝uintensity-(k1×upan)
其中：p为原始全色波段影像的灰度值，eintensity和epan分别为gs第一分量和全色波段影像灰度值的均值。
3.4)最后将修改后的全色波段作为第一分量，进行gs逆变换，输出n 1个波段，去除第一个波段，就是融合后的高分辨率多光谱影像结果，逆变换的计算公式如下：
4)基于高空间分辨率时序数据，计算归一化差异植被指数和土壤调节植被指数，生成植被指数时间序列数据集。其中，归一化差异植被指数ndvi由于计算简单，较好的消除太阳高度角、地形、天气等辐射变化的影响，已经成为植被生长状态及植被覆盖度最佳指示因子；土壤调节植被指数作为ndvi的后续，是基于ndvi和大量观测数据提出的，主要用于减少土壤背景的影响。具体包括如下步骤：
4.1)基于高空间时序数据进行归一化差异植被指数和土壤调节植被指数计算，得到植被指数高空间分辨率时序数据集；所述的两种植被指数的计算公式如下：
归一化差异植被指数ndvi：
ndvi＝(nir-red)/(nir red)
土壤调节植被指数savi：
savi＝(nir-red)*(1 l)/(nir red l)
式中，nir为近红外波段地表反射率值，red为红波段地表反射率值，l是随着植被密度变化的参数，取值范围是从0～1，当植被覆盖度很高时为0，很低时为1，通常取0.5时savi消除土壤反射率的效果最好。
4.2)将原始的时间序列数据进行归一化差异植被指数和土壤调节植被指数的计算，得到植被指数原始时序数据集；所述的两种植被指数的计算公式如下：
归一化差异植被指数ndvi：
ndvi＝(nir-red)/(nir red)
土壤调节植被指数savi：
savi＝(nir-red)*(1 l)/(nir red l)
式中，nir为近红外波段地表反射率值，red为红波段地表反射率值，l是随着植被密度变化的参数，取值范围是从0～1，当植被覆盖度很高时为0，很低时为1，通常取0.5时savi消除土壤反射率的效果最好。
4.3)基于步骤4.1和4.2的结果，随机选取区域内1～5万个点，对锐化前与锐化后不同植被指数进行相关性分析；随机选点是通过arcgis的随机产生点的工具进行，并通过多值提取至点，进行不同时相的点之间的相关性大小；相关性大小通过相关性系数r2表示，相关系数越大，证明融合前后影像灰度值差异就越小，融合精度越准确，反之亦然。
5)结合野外采样各地物的样本点，构建分类样本库；综合考虑不同地物的物候曲线特征差异性较大，先将总体样本按其类型进行分层，然后在每层中随机抽取样本，并将抽取的样本按一定比例分为训练样本与验证样本。其中训练样本集与验证样本集的抽取比例优选为2:1。
6)利用随机森林的分类方法对作物类型进行识别，并对比原始时间序列数据提取结果与锐化后的时间序列提取结果。随机森林就是通过集成学习的思想将多棵树集成的一种算法，它的基本单元是决策树，而它的本质属于机器学习的一大分支—集成学习，随机森林最主要的优势在于不需要特征降维，适合容量较大的数据集并且具有良好的鲁棒性。主要包括以下步骤：
6.1)将训练数据集进行划分，抽取10％作为训练模型的测试数据，根据测试数据集的准确率，选择合适的训练参数，例如：随机树的数量，叶子节点的数量以及最大叶深度等；
6.2)将训练好的模型应用至整个研究区，完成作物类型分类，包括锐化后的高空间时序数据集和锐化前的时序数据集等；
6.3)对比原始时间序列数据提取结果与锐化后的时间序列提取结果。
7)基于步骤6)得到的不同作物类型的结果，将提取的所有类型的作物进行整合，形成耕地要素数据，将所有的非作物类型进行整合，形成非耕地数据并完成精度评价，主要包括步骤：
7.1)将不同的作物类型进行整合，形成耕地数据；将非作物类型整合，形成非耕地数据；
7.2)将验证集中的不同作物类型的样本点合并成一类，即耕地样本点，并进行精度评价；
7.3)如精度验证结果如果与期望精度不符，则返回步骤5)重新进行分层抽样并重复步骤6)与步骤7)，直至与期望精度相符，最终完成耕地类型识别
7.1)将不同的作物类型进行整合，形成耕地数据；将非作物类型整合，形成非耕地数据。
7.2)将验证集中的不同作物类型的样本点合并成一类，即耕地样本点，并进行精度评价。精度验证的指标选择总体准确率(overallaccuracy)、精确率(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.3)如精度验证结果如果与期望精度不符，则返回步骤5)重新进行分层抽样并重复步骤6与步骤7，直至与期望精度相符，最终完成耕地类型识别。
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。
1.一种基于图像锐化的时序数据耕地提取方法，其特征在于，包括如下步骤：
步骤1)，输入单时相高空间分辨率遥感影像全色波段并进行数据预处理；
步骤2)，输入l2a级别的多时相时间序列sentinel-2影像并对其进行批量预处理；
步骤3)，对于步骤1)和步骤2)的结果，基于单时相高空间分辨率全色波段和多时相sentinel-2多光谱波段进行gram-schmidt图像锐化操作，生成高空间分辨率时序数据；
步骤4)，对于步骤3)中的高空间分辨率时序数据，计算归一化差异植被指数和土壤调节植被指数，生成高空间分辨率植被指数时间序列数据集；并通过对比锐化前与锐化后的时序数据进行植被指数相关性分析；
步骤5)，根据步骤4)得到的高空间分辨率植被指数时间序列数据集，结合野外采样各地物的样本点，构建分类样本库并进行分层抽样；
步骤6)，利用步骤5)中构建的分类样本库，先利用随机森林算法进行模型训练，通过测试集的表现确定最优模型参数；将训练好的模型应用至全区域，实现全区域作物类型的识别，并对比原始时间序列数据提取结果与锐化后的时间序列提取结果；
步骤7)，利用步骤6)得到的作物类型的结果，将提取的所有类型的作物进行整合，形成耕地要素数据，将所有的非作物类型进行整合，形成非耕地数据并完成精度评价；最终生成耕地资源要素空间分布数据。
2.根据权利要求1所述的一种基于图像锐化的时序数据耕地提取方法，其特征在于：所述步骤1)中，高空间分辨率全色波段影像预处理包括辐射校正、正射校正、几何配准以及影像裁剪，在进行几何配准时，依托于高精度谷歌数据、天地图开源数据，手动选择若干控制点进行几何纠正；步骤2中所述预处理是通过欧空局snap平台进行格式转换输出与批量裁剪处理。
3.根据权利要求1所述的一种基于图像锐化的时序数据耕地提取方法，其特征在于：所述步骤2)中，输入时间序列l2a级别的sentinel-2影像，并进行批量影像裁剪处理；其中，输入的sentinel-2的影像经过严格的云量控制，即通过整景云量小于n1％或者研究区云量小于n2％来进行数据的筛选；所利用到的时序sentinel-2是由蓝、绿、红以及近红外四个波段组成。
4.根据权利要求1所述的一种基于图像锐化的时序数据耕地提取方法，其特征在于：所述步骤3)具体包括如下步骤：
3.1)使用多波段低分辨率影像对单波段高分辨率影像进行模拟，模拟的方法根据光谱响应函数wi进行模拟，即模拟的全色波段影像灰度值
3.2)利用模拟的高分辨率波段影像作为gram-schmidt变换第一分量来对模拟的高分辨波段影像和低分辨率波段影像进行gram-schmidt变换，在进行gram-schmidt变换时，第t个变换分量由前t-1个变换分量构造，即：
其中，gst是gs变换后产生的第t个正交分量，bt是原始地空间分辨率多光谱影像的第t个波段影像，ut是第t个原始低空间分辨率多光谱波段影像灰度值的均值，i,j分别表示低空间分辨率影像的行数和列数，
其中，c是影像的列数，r是影像行数，协方差
其中，标准差σ的计算公式如下：
3.3)根据gs第一分量，即模拟波段的均值u和标准差σ对全色波段进行修改，修改公式如下：
p′＝(p×k1) k2
k2＝uintensity-(k1×upan)
其中：p为原始全色波段影像的灰度值，eintensity和epan分别为gs第一分量和全色波段影像灰度值的均值；
3.4)最后将修改后的全色波段作为第一分量，进行gs逆变换，输出n 1个波段，去除第一个波段，就是融合后的高分辨率多光谱影像结果，逆变换的计算公式如下：
5.根据权利要求1所述的一种基于图像锐化的时序数据耕地提取方法，其特征在于：所述步骤4)具体包括如下步骤：
4.1)基于高空间时序数据进行归一化差异植被指数和土壤调节植被指数计算，得到植被指数高空间分辨率时序数据集；两种植被指数的计算公式如下：
归一化差异植被指数ndvi：
ndvi＝(nir-red)/(nir red)
土壤调节植被指数savi：
savi＝(nir-red)*(1 l)/(nir red l)
式中，nir为近红外波段地表反射率值，red为红波段地表反射率值，l是随着植被密度变化的参数，取值范围是从0～1；
4.2)将原始的时间序列数据进行归一化差异植被指数和土壤调节植被指数的计算，得到植被指数原始时序数据集；
4.3)随机选取区域内若干个点，采用相关性系数r2对锐化前与锐化后不同植被指数进行相关性分析；随机选点是通过arcgis的随机产生点的工具进行，并通过多值提取至点，进行不同时相的点之间的相关性大小。
6.根据权利要求1所述的一种基于图像锐化的时序数据耕地提取方法，其特征在于：所述步骤5)是通过野外采集的各土地利用类型的样本点，建立该区域作物类型样本库，按照一定的比例进行分层抽样，获取训练数据集与验证数据集。
7.根据权利要求1所述的一种基于图像锐化的时序数据耕地提取方法，其特征在于：所述步骤6是通过随机森林的算法，实现不同作物类型的识别，并对比原始时间序列数据提取结果与锐化后的时间序列提取结果，包含以下步骤：
6.1)将训练数据集进行划分，抽取n3％作为训练模型的测试数据，根据测试数据集的准确率，选择合适的训练参数，包括随机树的数量，叶子节点的数量以及最大叶深度；
6.2)将训练好的模型应用至整个研究区，完成作物类型分类，包括锐化后的高空间时序数据集和锐化前的时序数据集；
6.3)对比原始时间序列数据提取结果与锐化后的时间序列提取结果。
8.根据权利要求1所述的一种基于图像锐化的时序数据耕地提取方法，其特征在于：所述步骤7具体包含以下步骤：
7.1)将不同的作物类型进行整合，形成耕地数据；将非作物类型整合，形成非耕地数据；
7.2)将验证集中的不同作物类型的样本点合并成一类，即耕地样本点，并进行精度评价；
7.3)若精度验证结果与期望精度不符，则返回步骤5)重新进行分层抽样并重复步骤6)与步骤7)，直至与期望精度相符，最终完成耕地类型识别。
9.根据权利要求8所述的所述的一种基于图像锐化的时序数据耕地提取方法，其特征在于：所述步骤7.3)中，精度验证的指标选择总体准确率、精确率、召回率以及f-score四个指标。
技术总结