本发明属于卫星遥感图像处理与应用的技术领域。
背景技术：
农作物分类信息的准确获取在预估产量、加强作物生产管理、农业保险与灾害评估等诸多方面中具有十分重要的科学与实践意义，近年来是遥感领域的研究热点之一。在诸多基于遥感影像的农作物分类算法中，深度学习方法被认为是一项突破性技术并以高效学习的特点在遥感农业识别领域掀起一阵热潮。近年来，国内外的研究学者提出了许多基于深度学习高分辨率遥感影像的作物分类方法。代表研究如下：
技术实现要素：
本发明采用的技术方案具体步骤如下：
步骤一、影像预处理：获取到的多时相高分辨率遥感影像是进行过几何精校正的正射影像，对该影像进行辐射定标与大气校正进而获取预处理之后的遥感影像。
(a)辐射定标的方法：辐射定标就是将图像的数字量化值(dn)转化为辐射亮度值、反射率或者表面温度等物理量的处理过程。辐射定标参数一般放在元数据文件中，利用完整的遥感图像处理平台(envi)的辐射定标工具(radiometriccalibration)能自动从元数据文件中读取参数，从而完成辐射定标。
(b)大气校正的方法：利用完整的遥感图像处理平台(envi)的大气校正工具(flassh)来消除由于大气影响所造成的辐射误差，反演地物真实的表面反射率。
步骤二、信息源：为了保证训练样本提取的准确性，基于原始多时相遥感影像提取光谱指数或特殊的光谱信息作为信息源，突出典型作物表现出的光谱特性。根据农作物不同的光谱特性，选择的信息源如下。
(a)可变大气阻力指数variableatmosphericallyresistantindex(varigreen)：光谱指数可以体现不同作物在可见光以及近红外波段的光谱反射特性，采用不同的光谱指数，可以获得表达作物不同生长发育阶段的信息。每个光谱指数都表达了一个特定的作物状态特征。因此，本发明在多时相遥感影像上对比了15种常用的光谱指数，发现7月份的varigreen可以对水稻部分进行较好的表达。考虑到水稻的下面是水，水在蓝绿波段的反射率较高并且作物在可见光通道具有不同的反射率。因此，可以增强对水稻的识别，varigreen的公式如式(1)所示：
(b)比值建筑用地指数ratiobuildingindex(rbi)：考虑到建筑物在遥感影像上的反射率较高，与作物在遥感影像上的反射率相差较多。因此，可以利用蓝波段与近红波段的比值来突出建筑物，增强对建筑物的识别，rbi的公式如式(2)所示：
(c)短波红外short-waveinfrared(swir)：通过大量的实验我们发现短波红外波段可以对不易区分的玉米与大豆两种作物进行有效的表达与区分。因此，本发明提取了多时相遥感影像的b11、b12，并通过pca降维，提取了多时相遥感影像短波红外波段的第一主成分。在该信息源上，城镇的反射率较高，玉米的反射率较低，大豆的反射率介于玉米与城镇之间，可以对玉米与大豆进行有效的区分。
(a)从原始多时相遥感影像中提取信息源，包括7月份的vraigreen，swirpca1与rbi。利用otsu分割算法对获取的信息源进行分割处理。otsu是一种确定图像分割阈值的算法，使用了聚类的思想，该方法又称最大类间方差法，有时也称为大津算法。该算法把图像的灰度数按照灰度级分为两个部分，使两个部分的灰度值差异最大，每个部分之间的灰度差异最小，通过方差的计算寻找合适的灰度级别进行划分。
(b)根据分割结果在原始信息源上采用直方图统计的方法并根据设置的置信区间自动确定各类作物提取的阈值范围，进而获取粗训练样本。
(c)人工选择的样本通常偏向于均质区域并且忽略异常值。在本发明中，为每类自动挑选合适的面积阈值用于候选样本，删除面积小于设定阈值的对象，用于去除小的异构点。
(d)去除各类样本的重叠区域。考虑到各类提取样本重叠部分的不确定性，将多于一类标签的训练样本删除，进而获取初步提纯的样本数据。
(e)为了进一步提纯样本数据集，保证样本挑选的准确性同时出于对各类样本边缘混合点的考虑引入滑窗的设计，自动化的挑选出合适的滑窗尺寸。在挑选样本时我们首先保证滑窗内数据标签一致，不包括外围一圈的数据，本发明选择内部数据用于构建样本数据集。相当于在样本周围扩充了一圈具有相同标签的数据，用来保证边缘混合点不会对纯净样本的提取造成影响。
自动化影像块训练样本提取规则：首先设计不同的滑窗尺寸与影像块大小，通过移动滑窗获取满足条件的训练样本，滑窗的移动步长stride＝1，影像块大小与影像块样本滑窗尺寸之间遵循的关系如下：(s×b+p)×(s×b+p)＝n×n。p＝2代表在影像块样本周围扩充了一圈具有相同标签的数据，b×b代表影像块大小，s＝1,2,3,…代表滑窗内部的样本数据被划分为s×s个大小为b×b的影像块，n×n代表挑选影像块样本的滑窗尺寸。
自动化像素点训练样本提取规则：挑选像素点样本的滑窗尺寸为m×m，保证滑窗内部的数据具有相同的标签，选择(m-1)×(m-1)大小的影像块制作像素点训练样本，为了保证像素点样本选择的不重复性，设计规则如下。用(i,j)(i＝1,2,…,(m-m+1)；j＝1,2,…,(n-m+1))表示输入影像某位置的横纵坐标，输入影像高度m，宽度n。
(a)如果i＝1，j＝1，提取(i+1:i+(m-1),j+1:j+(m-1))范围内的像素点样本
(b)如果i＝1，j≠1，提取(i+1:i+(m-1),j+(m-1))范围内的像素点样本
(c)如果i≠1，j＝1，提取(i+(m-1),j+1:j+(m-1))范围内的像素点样本
(d)如果i≠1，j≠1，提取(i+(m-1),j+(m-1))位置的像素点样本
上述滑窗规则的设计实现了全自动化影像块训练样本与像素点训练样本的提取过程，为实现最优的作物制图结果提供了可靠的样本数据。
步骤四、基于影像块分类的深度学习网络：采用全自动化方法获取的训练样本与geo-3dcnn网络对多时相的高分辨率遥感影像基于影像块分类，得到作物制图结果。
geo-3dcnn网络结构中包括输入层、卷积层、池化层、全连接层与输出层。我们采用广泛被应用的vgg网络为模板去训练一个深度卷积神经网络，其中的所有二维卷积核将被三维卷积核替代。网络的层数与每层的内核数常取决于输入的容量和问题的复杂度。
输出特征图大小的计算方式有两种，一种是padding＝'same'，另一种是padding＝'valid'。若采用'same'方式，若采用'valid'方式，其中，win是输入特征图边的大小，wout是输出特征图边的大小，padding为特征图填充的圈数，strides代表步长，kernel_size代表内核大小。
所述的geo-3dcnn网络结构如下：
(a)geo-3dcnn网络的第0层为输入层，输入图像大小为10@3×3×3，10代表输入遥感影像数据包括10个通道，3代表输入影像的时相数，3×3代表输入网络分类的影像块大小。
(b)输入层之后为geo-3dcnn网络的第一部分，geo-3dcnn网络的第一部分由卷积层conv3d、激活层activation和池化层maxpool3d组成，卷积层设定的参数如下：
滤波器的个数filters＝32，
内核大小kernel_size＝[2，2，2]，
输出特征图大小的计算方式为padding＝'same'，
填充padding＝1，
步长strides＝[1，1，1]，
输入影像经过卷积层后，时相方向：高度与宽度方向：即该卷积层输出特征响应图的尺寸为32@3×3×3，
接下来为激活层activation，使用的激活函数是relu，
激活层activation之后为池化层maxpool3d，池化层设定参数如下：
步长strides＝[1，2，2]，
池化方式为最大池化，时相方向：高度与宽度方向：即池化层输出特征响应图的尺寸为32@3×2×2。
(c)geo-3dcnn网络的第二部分由卷积层conv3d、激活层activation和池化层maxpool3d组成，卷积层设定的参数如下：
滤波器的个数filters＝64，
输入影像经过卷积层后，时相方向：高度与宽度方向：即该卷积层输出特征响应图的尺寸为64@3×2×2，
池化方式为最大池化，时相方向：高度与宽度方向：即池化层输出特征响应图的尺寸为64@3×1×1。
(d)geo-3dcnn网络的第三部分由卷积层conv3d、激活层activation和池化层maxpool3d组成，卷积层设定的参数如下：
滤波器的个数filters＝128，
输入影像经过卷积层后，时相方向：高度与宽度方向：即该卷积层输出特征响应图的尺寸为128@3×1×1，
池化方式为最大池化，时相方向：高度与宽度方向：即池化层输出特征响应图的尺寸为128@3×1×1。
(e)geo-3dcnn网络的第四部分是flatten层，flatten层的大小为384×1，提取每个影像块训练样本的中心经纬度坐标，并与flatten层进行拼接，融合空间地理信息，经过融合操作的flatten层输出为386×1。
(f)geo-3dcnn网络的第五部分由两个连续的全连接层dense组成，最后一个全连接层(输出层)输出4×1特征向量(对应于水稻，城镇，玉米和大豆这四种类别的概率)。其中，第一个全连接层dense的输入个数为386×1，权重系数w为386×128，偏置b为128，输出个数为128×1；第二个全连接层dense的输入个数为128×1，权重系数w为128×4，偏置b为4，输出个数为4×1。
步骤五、基于像素点分类的深度学习网络：采用全自动化方法获取的训练样本与geo-conv1d网络对多时相的高分辨率遥感影像基于像素点分类，得到作物制图结果。
geo-conv1d网络是cnn的一种特殊形式，它通过卷积层、池化层、全连接层、dropout层与输出层来实现。
所述的geo-conv1d网络结构如下：
(a)geo-conv1d网络的第0层为输入层，输入图像大小为30×1，30代表了多时相遥感影像输入的总波段数。
(b)输入层之后为geo-conv1d网络的第一部分，geo-conv1d网络的第一部分由卷积层conv1d和激活层activation组成，卷积层设定的参数如下：
内核大小kernel_size＝[3]，
输出特征图大小的计算方式为padding＝'valid'，
步长strides＝[1]，
输入影像经过卷积层后，即该卷积层输出特征响应图的尺寸为28×64，
接下来为激活层activation，使用的激活函数是relu。
(c)geo-conv1d网络的第二部分是inception模块与dropout层，inception模块由卷积层conv1d、激活层activation、池化层maxpool1d与concatenate层组成，卷积层设定的参数如下：
内核大小kernel_size＝[3]和[5]，
输入影像经过卷积层后，和即卷积层输出特征响应图的尺寸为26×128和24×128，
池化层设定参数如下：
内核大小kernel_size＝[2]，
步长strides＝[2]，
池化方式为最大池化，即池化层输出特征响应图的尺寸为14×128，
接下来是concatenate层，将两个卷积层conv1d与一个池化层maxpool1d进行集成，concatenate层输出特征响应图的尺寸为64×128，
然后连接dropout层，dropout是一种随机丢弃一些神经元的正则化技术，丢弃神经元的比例设置为40％。
(d)geo-conv1d网络的第三部分由卷积层conv1d、激活层activation和dropout层组成，卷积层设定的参数如下：
输入影像经过卷积层后，即该卷积层输出特征响应图的尺寸为62×128，
(e)geo-conv1d网络的第四部分由卷积层conv1d、激活层activation和dropout层组成，卷积层设定的参数如下：
滤波器的个数filters＝256，
输入影像经过卷积层后，即该卷积层输出特征响应图的尺寸为60×256，
(f)geo-conv1d网络的第五部分是flatten层，flatten层的大小为15360×1，提取每个影像块训练样本的中心经纬度坐标，并与flatten层进行拼接，融合空间地理信息，经过融合操作的flatten层输出为15362×1。
(g)geo-conv1d网络的第六部分由两个全连接层dense与dropout层组成，最后一个全连接层(输出层)输出4×1特征向量(对应于水稻，城镇，玉米和大豆这四种类别的概率)。其中，第一个全连接层dense的输入个数为15362×1，权重系数w为15362×512，偏置b为512，输出个数为512×1，
然后连接dropout层，dropout是一种随机丢弃一些神经元的正则化技术，丢弃神经元的比例设置为50％，
第二个全连接层dense的输入个数为512×1，权重系数w为512×4，偏置b为4，输出个数为4×1。
步骤六、分类融合方法：基于像素点的作物分类方法可以实现较高分辨率的作物制图，但是伴随着椒盐噪声现象的产生，而基于影像块的分类方法可以减少椒盐噪声现象，但是作物制图的分辨率会降低。为了充分利用两种分类方法的优势，提出一种分类融合方法。
(a)根据本发明提出的全自动化训练样本选择方法准备好样本数据集，包括像素点纯样本与影像块纯样本。
(b)使用当前的影像块样本与像素点样本分别训练geo-3dcnn网络与geo-conv1d网络。
(c)利用训练好的geo-3dcnn网络标记所有待分类的影像块。
(d)根据阈值条件(3)为每类挑选出分类不准确的影像块样本。满足阈值条件的影像块样本被保留，不满足阈值条件的影像块输入geo-conv1d网络进行基于像素点的分类，实现两种分类结果的融合。
b/a≤t1且a≥t2(3)
在式(3)中，a和b代表当前模型给出任意训练样本不同类别的最大得分与第二大得分。如果b/a的比值大于给定的阈值t1，则表明分类器不能对两种类别进行可靠的区分。如果a小于给定阈值t2，则表示对类别的识别不准确。
步骤七、精度评估：在测试数据集上评估每个深度学习网络的农作物分类准确性。本发明应用了基于像素评估的六个评价参数，分别是总体分类精度(oa)，kappa系数，准确率(precision)，召回率(recall)，f1和iou来评估作物分类的准确性。oa指被正确分类的类别像元数与总的类别像元数的比值，oa可以较好的表征分类精度，其计算公式如式(4)所示。
其中，pi,j代表属于第i类而被分配为第j类的像元总数，n代表类别的数目。
kappa系数用于一致性检验，同样用来衡量分类精度，kappa系数的计算是基于混淆矩阵的，其计算公式如式(5)所示。
其中，n代表样本总数，a1,a2,...,an是每种类型的实际样本数，b1,b2,...,bn是每种类型预测的样本数。
通过与地面真实参考数据(gt)进行对比，真阳性(tp)，假阳性(fp)和假阴性(fn)分别代表正确提取的类，错误提取的类和漏分类的数量。precision用于表示正确分类的类别数目在该类别提取结果中所占的比例，recall用于表示被正确提取的类别数占该类别真实目标像素的比例。计算公式分别如式(6)与式(7)所示。
f1得分与iou是综合性的评价指标，f1得分是precision与recall的调和平均数，iou描述了作物分类结果和gt之间的重叠率。计算公式分别如式(8)与式(9)所示。
本发明的有益效果：
附图说明
图1是本发明使用的研究区。
图2是本发明基于时空深度学习融合技术的农作物全自动化分类方法的流程图。
图3是本发明实施例1的信息源otsu分割图。
图4是本发明实施例1的样本提取结果。
图5是本发明实施例1的geo-3dcnn网络结构示意图。
图6是本发明实施例1的geo-conv1d网络结构示意图。
图7是本发明实施例1的深度学习网络融合前后的作物识别结果图。
具体实施方式
实施例1：
表1
图2展示了本发明方法的基本步骤。
步骤一，多时相遥感影像预处理。
步骤二，从原始的多时相sentinel-2影像中提取有用的信息源用于自动化训练样本的选择。
步骤三，全自动化训练样本选择，包括otsu分割、直方图统计、设计面积阈值、去除样本的重叠区域以及对滑窗与影像块大小的设计与选择。
步骤四，设计融合地理信息的geo-3dcnn网络，并利用自动化选择的影像块训练样本训练网络。
步骤五，设计融合地理信息的geo-conv1d网络，并利用自动化选择的像素点训练样本训练网络。
步骤六，为了集成基于像素点分类方法与基于影像块分类方法的优势，采用activelearning策略实现两种分类结果的融合。
步骤七，对深度学习网络分类融合前后的作物识别结果进行精度分析。
本发明方法的具体步骤如下：
步骤一：影像预处理
获取到的多时相高分辨率sentinel-2影像是经过几何精校正的正射影像，对该影像进行辐射定标与大气校正进而获取预处理之后的遥感影像。
(a)辐射定标：利用完整的遥感图像处理平台(envi)的辐射定标工具(radiometriccalibration)完成数字量化值(dn)到辐射亮度值的转化。
(b)大气校正：利用完整的遥感图像处理平台(envi)的大气校正工具(flassh)来消除由于大气影响所造成的辐射误差。
步骤二、信息源
为了保证训练样本提取的准确性，基于原始多时相遥感影像提取光谱指数或特殊的光谱信息作为信息源，突出典型作物表现出的光谱特性。根据农作物不同的光谱特性，本发明发现7月份的varigreen可以突出水稻的特性，rbi可以突出建筑物的特性，而多时相遥感影像短波红外波段的第一主成分(swirpca1)可以突出玉米与大豆的特性。因此，本发明选择了上述三种信息源用于后续训练样本的自动化提取。
步骤三、全自动化训练样本选择
(a)从原始多时相遥感影像中提取信息源，包括7月份的vraigreen，swirpca1与rbi。利用otsu分割算法对获取的信息源进行分割处理，信息源otsu分割结果图如图3所示，(a)是7vraigreen的分割图，(b)是swirpca1的分割图，(c)是rbi的分割图。
(b)根据分割结果在原始信息源上采用直方图统计的方法并根据设置的置信区间自动确定各类作物提取的阈值范围(tr,tu,tc,ts)，进而获取粗训练样本。
(c)人工选择的样本通常偏向于均质区域并且忽略异常值。在本发明中，为每类自动挑选合适的面积阈值(ar、au、ac、as)用于候选样本，删除面积小于设定阈值的对象，用于去除小的异构点。
影像块训练样本最优滑窗选择结果：滑窗移动步长stride＝1，影像块大小与影像块样本滑窗尺寸之间遵循的关系如下：(s×b+p)×(s×b+p)＝n×n(s＝1,2,…,6；l＝2,3,...,6；p＝2)。本发明通过实验发现基于影像块分类的最优滑窗尺寸nopt×nopt为17×17，最优影像块大小bopt×bopt为3×3。
像素点训练样本最优滑窗选择结果：基于像素点的训练样本选择，根据实际情况设计滑窗大小m×m(m＝6,7,…,14)，用来挑选基于像素点分类的最优滑窗尺寸。本发明通过实验发现基于像素点分类的最优滑窗尺寸mopt×mopt为10×10。
图4是多种信息源经过otsu分割与直方图统计获取到的粗样本如(a)所示，经过面积阈值的限定与去除重叠区域获得的初步提纯样本如(b)所示，最后通过最优滑窗与最优影像块尺寸的设计挑选出的影像块纯净样本如(c)所示，像素点纯样本如(d)所示，我们可以看到通过一种自动化的方式正确的标记了各类的纯样本。
步骤四：geo-3dcnn网络
geo-3dcnn网络结构中包括输入层、3个卷积层、3个池化层、2个全连接层(最后1个全连接层对应输出层)。我们采用广泛被应用的vgg网络为模板去训练一个深度卷积神经网络，其中的所有二维卷积核将被三维卷积核替代。研究了一种基于影像块的作物分类方法，图5是geo-3dcnn网络结构示意图。
geo-3dcnn网络结构如下：
步骤五：geo-conv1d网络
geo-conv1d网络包括输入层、3个卷积层、1个inception模块(包括2个卷积层、1个池化层与1个concatenate层)、dropout层、2个全连接层(最后1个全连接层对应输出层)。图6是geo-conv1d网络结构示意图。
geo-conv1d网络结构如下：
步骤六、分类融合方法
利用activelearning策略对geo-3dcnn网络与geo-conv1d网络获得的作物识别结果进行融合，具体实现流程如下：
在式(3)中，a和b代表当前模型给出任意训练样本不同类别的最大得分与第二大得分。如果b/a的比值大于给定的阈值t1(例如：0.1)，则表明分类器不能对两种类别进行可靠的区分。如果a小于给定阈值t2(例如：0.9)，则表示对类别的识别不准确。
步骤七、精度评估
采用总体分类精度(oa)，kappa系数，准确率(precision)，召回率(recall)，f1和iou这六个评估参数对深度学习网络融合前后的作物识别结果进行精度分析。