CN111738329B - 一种面向时间序列遥感影像的土地利用分类方法 - Google Patents

Abstract

本发明公开了一种面向时间序列遥感影像的土地利用分类方法，包括以下步骤：对组成时间序列遥感数据的多光谱影像进行主成分分析，获取三个主成分的影像；对每个三波段影像进行预训练，提取出特征图像；特征图像按时间顺序依次输入半监督卷积长短期记忆网络模型进行训练；利用训练好的模型对最后一时相的影像进行预测分类，获得分类结果。本发明综合考虑了时间序列的时间上下文信息、遥感影像的空间和光谱特征信息，利用预训练模型以及半监督分类学习的方式，较好地解决了现有技术存在的问题，使得所述的土地利用分类方法更加适用于训练样本数据较少，且遥感影像数据缺失或者部分缺失的分类场景，获得了较好的土地利用分类结果。

Description

一种面向时间序列遥感影像的土地利用分类方法

技术领域

本发明属于遥感影像分析技术领域，涉及一种面向时间序列遥感影像的土地利用分类方法。

背景技术

随着遥感技术的不断发展，不断增加的多源遥感数据和累积的历史影像，使得多源、多时相对地观测的遥感数据更容易被获取。不断发展的技术和增长的需求使影像分类技术从单一特征转变为多特征的应用。近年来的研究证明，时间上下文特征有利于提高影像分类的准确性，如Kun J&Shunlin L等人利用时间序列数据对森林覆盖的提取分类优于单时相分类(参考文献：Jia K,Liang S,Zhang L,et al.Forest cover classificationusing Landsat ETM+data and time series MODIS NDVI data[J].InternationalJournal of Applied Earth Observation and Geoinformation,2014,33:32-38.)。由于时间序列数据集的高维度特征及不断智能化处理的需求增长，基于深度学习理论的方法已成为较为热门且表现较好的一类分类方法，如M.Rubwurm等人利用深度学习中的长短期记忆网络模型(LSTM)对多时相遥感数据进行土地利用分类，且分类精度优于卷积神经网络(CNN)及非深度学习的支持向量机(SVM)。基于深度学习的分类方法多种多样，但由于这类方法通常需要丰富多样的样本数据用于网络模型的训练，这对于长时间序列遥感数据集的收集增加了难度，尤其是标签数据，在实际中很难保证时间序列中每张影像都有对应的真实地表标签数据。此外，光学遥感影像的采集难以避免云雾遮挡的情况，由于云雾遮挡的不可避免和不确定性，也会导致连续影像的缺失，不仅阻碍了时间序列遥感影像的收集，还会影响分类的效果和精度。

根据目前的研究背景和技术，对于时间序列遥感数据的分类而言，主要存在以下问题有待解决：(1)时间序列遥感数据集的收集较为困难，尤其是在实际中缺乏相应的标签样本，严重阻碍了依赖训练样本的分类模型的发展和分类效果。(2)光学遥感影像的采集难以避免有云雾遮挡的情况，由于云雾遮挡的不可避免和不确定性，也会导致连续数据的缺失，这不仅阻碍了时间序列遥感影像的收集，也影响了分类的效果和精度。

发明内容

本发明的目的在于提供一种面向时间序列遥感影像的土地利用分类方法，本发明是一种端到端的分类方法，综合考虑了时间序列的时间上下文信息、遥感影像的空间和光谱特征信息，利用深度学习、预训练模型以及半监督分类的联合方式，较好地解决了现有技术存在的问题，使得所述的土地利用分类方法更加适用于训练样本的标签数据较少，且遥感影像数据缺失或者部分缺失的分类场景，获得了较好的土地利用分类结果。

本发明的目的是这样实现的，一种面向时间序列遥感影像的土地利用分类方法，包括以下步骤：

步骤1，对组成时间序列遥感数据的多光谱影像进行主成分分析，获取三个主成分的影像；

步骤2，对每个三波段影像进行预训练，提取出特征图像；

步骤3，特征图像按时间顺序依次输入半监督卷积长短期记忆网络模型进行训练；

步骤4，利用训练好的模型对最后一时相的影像进行预测分类，获得分类结果；

在步骤1中，所述的时间序列遥感数据集Γ＝(X₁,X₂,X₃,...,X_N)中的所有多光谱影像X_n(n＝1,2,3,...,N)的尺寸为h×w×c，h、w、c分别表示影像的长、宽、和波段数，经过主成分分析后的尺寸为h×w×3，成为三波段影像。

步骤2中所述的预训练过程为：选择在开源的ImageNet大型自然影像数据库上完成训练后的ResNet模型作为预训练模型，提取出预训练模型中浅层的特征图像，尺寸为h'×w'×c'，其中，h'＜h，w'＜w，c'＞c，为了恢复至原始影像大小的长和宽，采用双线性差值算法对提取后的特征图进行上采样，得到尺寸为h×w×c′的特征图像，再将上采样后的特征图像与原始影像的多光谱波段进行融合，以得到尺寸为h×w×(c+c′)的数据；

步骤3中所述的半监督卷积长短期记忆网络模型是基于长短期记忆网络的基础上，将输入到状态和状态到状态之间的全连接传输过程转换为卷积运算，即ConvLSTM单元，在ConvLSTM单元中，随着时间序列数据中每个时刻数据的输入X_n(n＝1,2,3,...,N)，网络记忆单元状态S_n，以及当前隐藏层状态值H_n，由三个门机制控制，即输入门i_n、遗忘门f_n和输出门o_n，并不断被保存和更新，计算公式为：

i_n＝σ(W_i,X*X_n+W_i,H*H_n-1+W_i,S·S_n-1+b_i)

f_n＝σ(W_f,X*X_n+W_f,H*H_n-1+W_f,S·S_n-1+b_f)

o_n＝σ(W_o,X*X_n+W_o,H*H_n-1+W_o,S·S_n+b_o)

H_n＝o_n·tanh(S_n)

其中，S_n-1,H_n-1分别表示上一时刻的记忆单元状态和隐藏层状态，函数σ和tanh分别表示sigmoid和双曲正切函数，W_i，X，W_i，H，W_i，S表示输入门i_n分别与当前的输入X_n、隐藏层H_n和记忆单元状态S_n之间的权重矩阵，W_f，X，W_f，H，W_f，S表示遗忘门f_n分别与当前的输入、隐藏层和记忆单元状态之间的权重矩阵，W_o，X，W_o，H，W_o，S表示输出门o_n分别与当前的输入、隐藏层和记忆单元状态之间的权重矩阵，W_S，X，W_S，H分别表示记忆状态与输入和隐藏层之间的权重矩阵，b_i，b_f，b_o，b_S对应各个门和网络记忆单元的偏置项，

表示当前输入激活的记忆状态，符号‘*’及‘·’分别表示卷积运算和Hadamard矩阵乘积。

在训练过程中，由隐藏层状态计算出当前时间序列数据输入所对应的预测值，并利用已知和未知标签的联合信息计算损失函数，使其收敛以训练并优化网络模型，其损失函数由两部分组成，即利用已知标签信息的标准交叉熵，和联合已知和未知标签信息的无监督损失函数，公式可表示为：

其中，B表示每次输入的批次量，L代表有标签的训练样本，

表示该训练样本所对应的真实标签，z_i表示时间序列输入所对应的预测标签经softmax函数转化后的输出概率，

表示与同一个时间序列数据输入经过数据增强或时序变化后所对应的预测标签概率，C为类别数量，w(t)为一个随时间变化的变量，以权衡两部分损失函数的占比。

在步骤4中，利用训练好的模型计算出待测试影像的预测标签概率，并采用argmax函数取概率最大值的类别作为分类结果。

本发明方法提供了一种面向时间序列遥感影像的土地利用分类方法，所述方法综合考虑了时间序列的时间上下文信息，数据的时间、空间和光谱特征信息，完成了在少量已知标签的条件下有效的训练过程，同时，在一定程度上能抵抗云覆盖对光学遥感数据的影响，以实现精准的土地利用分类结果。

附图说明

图1本发明方法的流程示意图；

图2是本发明实施例中网络模型的整体结构流程图；

图3是卷积长短期记忆网络模型结构图，(a)可视为将网络模型展开的流程图，(b)图为一个神经单元的结构图；

图4(a)、(b)、(c)分别对应Data_1、Data_2、Data_3实验数据集下不同方法的分类结果，(d)图展示了分类结果的图例；

图5(a)、(b)、(c)分别对应实验区域的时间序列影像数据集；

图6是本实施例中三组实验数据集的时间分布情况；

图7是本实施例中改变时间序列长度值对分类的影响示意图；

图8是本实施例中不同标签数据占比进行训练的实验结果(以Kappa系数为例)；

图9本实施例中三组实验数据在不同云覆盖量下的分类精度结果。

具体实施方式

下面结合实施例和附图对本发明作进一步的说明，但不以任何方式对本发明加以限制，基于本发明教导所作的任何变换或替换，均属于本发明的保护范围。

如图1所示，一种面向时间序列遥感影像的土地利用分类方法，包括以下步骤：

步骤1，对组成时间序列遥感数据的多光谱影像进行主成分分析，获取三个主成分的影像；

步骤2，对每个三波段影像进行预训练，提取出特征图像；

步骤3，特征图像按时间顺序依次输入半监督卷积长短期记忆网络模型进行训练；

步骤4，利用训练好的模型对最后一时相的影像进行预测分类，获得分类结果；

在步骤1中，所述的时间序列遥感数据集Γ＝(X₁,X₂,X₃,...,X_N)中的所有多光谱影像X_n(n＝1,2,3,...,N)的尺寸为h×w×c，h、w、c分别表示影像的长、宽、和波段数，经过主成分分析(PCA)后的尺寸为h×w×3，成为三波段影像。

步骤2中所述的预训练过程为：选择在开源的ImageNet大型自然影像数据库上完成训练后的ResNet模型作为预训练模型，提取出预训练模型中浅层的特征图像，尺寸为h'×w'×c'，其中，h'＜h，w'＜w，c'＞c，为了恢复至原始影像大小的长和宽，采用双线性差值算法对提取后的特征图进行上采样，得到尺寸为h×w×c′的特征图像，再将上采样后的特征图像与原始影像的多光谱波段进行融合，以得到尺寸为h×w×(c+c′)的数据；

步骤3中所述的半监督卷积长短期记忆网络模型是基于长短期记忆网络的基础上，将输入到状态和状态到状态之间的全连接传输过程转换为卷积运算，即ConvLSTM单元，在ConvLSTM单元中，随着时间序列数据中每个时刻数据的输入X_n(n＝1,2,3,...,N)，网络记忆单元状态S_n(即长期记忆)，以及当前隐藏层状态值H_n(即短期记忆)，由三个门机制控制，即输入门i_n、遗忘门f_n和输出门o_n，并不断被保存和更新，计算公式为：

i_n＝σ(W_i,X*X_n+W_i,H*H_n-1+W_i,S·S_n-1+b_i)

f_n＝σ(W_f,X*X_n+W_f,H*H_n-1+W_f,S·S_n-1+b_f)

o_n＝σ(W_o,X*X_n+W_o,H*H_n-1+W_o,S·S_n+b_o)

H_n＝o_n·tanh(S_n)

其中，S_n-1,H_n-1分别表示上一时刻的记忆单元状态和隐藏层状态，函数σ和tanh分别表示sigmoid和双曲正切函数，W_i，X，W_i，H，W_i，S表示输入门i_n分别与当前的输入X_n、隐藏层H_n和记忆单元状态S_n之间的权重矩阵，W_f，X，W_f，H，W_f，S表示遗忘门f_n分别与当前的输入、隐藏层和记忆单元状态之间的权重矩阵，W_o，X，W_o，H，W_o，S表示输出门o_n分别与当前的输入、隐藏层和记忆单元状态之间的权重矩阵，W_S，X，W_S，H分别表示记忆状态与输入和隐藏层之间的权重矩阵，b_i，b_f，b_o，b_S对应各个门和网络记忆单元的偏置项，

表示当前输入激活的记忆状态，符号‘*’及‘·’分别表示卷积运算和Hadamard矩阵乘积。

在训练过程中，由隐藏层状态计算出当前时间序列数据输入所对应的预测值，并利用已知和未知标签的联合信息计算损失函数，使其收敛以训练并优化网络模型，其损失函数由两部分组成，即利用已知标签信息的标准交叉熵，和联合已知和未知标签信息的无监督损失函数，公式可表示为：

其中，B表示每次输入的批次量，L代表有标签的训练样本，

表示该训练样本所对应的真实标签，z_i表示时间序列输入所对应的预测标签经softmax转化后的输出概率，

表示与同一时间序列数据经过数据增强或时序变化后所对应的预测标签概率，C为类别数量，w(t)为一个随时间变化的变量，以权衡两部分损失函数的占比。

在步骤4中，利用训练好的模型计算出待测试影像的预测标签概率，并采用argmax函数取概率最大值的类别作为分类结果。

在具体实施例中，针对时间序列遥感影像分类问题，在模型的训练过程中，输入为按照时间顺序排列的遥感影像Γ＝(X₁,X₂,X₃,...,X_N)组成的时间序列数据，和与之所对应的标签数据集(D₁,D₂,D₃,...,D_N)。其中，时间序列数据的时间总长度为N，X_n表示以时间排列顺序的第n时刻下的影像，D_n表示在n时刻影像所对应的标签数据。时间序列数据是一个h×w×N×c的四维数据，h，w，c分别表示影像的长、宽、和深度(即遥感影像的波段数)。与时间序列影像所对应的标签数据集(D₁,D₂,D₃,...,D_N)是有混合了有类别标签和无标签的序列数据。有类别标签的数据可以根据国际分类标准和实际地区情况确定并标注，或利用已有的土地利用数据产品，使影像中的每一个像素都对应一个真实类别值，而无标签的数据用0值填充，因此，每个标签数据都是一个h×w的二维数据。为了验证本方法能面向时间序列遥感影像Γ＝(X₁,X₂,X₃,...,X_N)，利用少量标签样本，对最后时刻的影像X_N实现像素级别的准确分类，将对输入的时间序列数据将划分为两个部分，即除去最后时刻影像的其余影像作为训练集，而最后时相的数据作为测试数据。本网络模型的训练过程可拆分为两个阶段，如图2所示，具体原理将详细阐述。

(1)预训练模型

在将深度学习技术应用到遥感影像处理时，由于地表覆盖类型多样且变化复杂，直接应用于复杂地表覆盖的结果不够稳定。因此，为了加强数据输入的特征表达，使用了在其他图像处理任务下已经完成训练的神经网络模型作为特征提取器，以实现预训练过程。此处选用了在ImageNet大型自然影像数据库上已完成训练过程的ResNet模型作为预训练特征模型。预训练特征的过程可以用以下公式表示。首先，需要对时间序列数据中的每张原始影像X_n(n＝1,2,3,...,N)的多光谱维度作降维处理，利用主成分分析从原始多光谱影像提取三个主成分以作为影像的三个波段。接着，送入选定的已经完成训练的网络模型中，并提取出模型中间某一层的特征图。本发明方法选择了ResNet神经网络模型第一层卷积后的特征图，即利用7×7的卷积核由主成分分析后的3个波段数增加至64，并且，此时输出的特征图的长和宽也都发生了变化。为了恢复至与原始影像输入相同的长和宽，将提取的特征图利用双线性插值法进行上采样过程F_up，得到上采样后的特征图。最后，将上采样后的特征图与原始多光谱影像X_n进行通道融合F_cat，以得到预训模型后的特征值

由于选定的ResNet网络在大型自然影像数据库的训练下具有了一定的空间表达能力，经过预训练处理后的特征值不仅保持了原始影像的空间结构，即相同尺度的长和宽，其深度信息在原有基础上还得到了扩展和增强，即深度为(c+64)。如此一来，能有效地利用遥感影像的光谱、空间和纹理信息，有利于遥感影像的分类。

其中，ResNet¹表示从ResNet网络第一层卷积层中提取的特征值。

(2)基于半监督的卷积长短期记忆网络模型

大量的研究表明，长短期记忆网络(LSTM)模型是一种特殊的循环神经网络结构，具有三个门机制，即输入门、遗忘门和输出门，可以有效地解决时序化问题，已经广泛地应用于自然语言处理、机器翻译、视频分析等领域。而针对时间序列遥感影像的分类模型而言，引用目前已提出的具有卷积结构的长短期记忆网络(ConvLSTM)可以充分利用数据的多特征时空信息，不仅有效地提取时间上下文信息，还能有效地表达光学遥感影像的空间、纹理和光谱特征信息。而ConvLSTM是在LSTM的基础上，将LSTM结构中输入与状态、状态与状态之间的全链接结构替换成可处理三维影像的卷积结构，使得ConvLSTM既能有效处理时间特征，也具有一定的空间表达能力。

为了计算机的运算效率和模型的训练效果，需要对总时间长度为N的时间序列数据集进行数据的划分：以时间长度为单位，按时间顺序依次往后移动，从而划分为若干个子时间序列数据集：

{(X₁,X₂,...,X_t),(X₂,X₃,...,X_(t+1)),...,(X_(n-t+1),X_(n-t+2),...,X_n),...,(X_(N-t+1),X_(N-t+2),...,X_N)}

如图2所示，由于“多对一”的分类网络模型，一个时间序列输入只有一个输出，即每一个子序列数据将对应一个输出标签数据(子序列数据中最后时刻影像的预测标签值)。相应地，上述划分的测试集和训练集用子序列数据表示更清晰，即包含最终时刻影像的子序列数据(X_(N-t+1),X_(N-t+2),...,X_N)为测试集，其余以t为子序列数据为训练集的时间单位长度。

以子序列数据(X_(n-t+1),X_(n-t+2),...,X_n)为例，ConvLSTM网络模型中的主要运算原理如下公式所示，其中，i_n,f_n,o_n分别表示当前n时刻的输入门、遗忘门和输出门结构，随着子序列数据中每个数据的依次输入，网络记忆单元状态S_n(长期记忆)和当前隐藏层状态值H_n(短期记忆)由三个门机制控制，并不断被保存和更新。不同的门与输入、隐藏层之间的权重矩阵是共享的，使得在输入下一时刻数据后的记忆状态得到了更新，也保留了上一时刻运算的记忆状态。H_n-1，S_n-1分别表示上一时刻的隐藏层状态和记忆状态，函数σ和tanhO分别表示sigmoid和双曲正切函数，W_i，X，W_i，H，W_i，S表示输入门i_n分别与当前的输入X_n、隐藏层H_n和记忆单元状态S_n之间的权重矩阵，W_f，X，W_f，H，W_f，S表示遗忘门f_n分别与当前的输入、隐藏层和记忆单元状态之间的权重矩阵，W_o，X，W_o，H，W_o，S表示输出门o_n分别与当前的输入、隐藏层和记忆单元状态之间的权重矩阵，W_S，X，W_S，H分别表示记忆状态与输入和隐藏层之间的权重矩阵，b_i，b_f，b_o，b_S对应各个门和网络记忆单元的偏置项，

表示当前输入激活的记忆状态，符号‘*’及‘·’分别表示卷积运算和Hadamard矩阵乘积。

i_n＝σ(W_i,X*X_n+W_i,H*H_n-1+W_i,S·S_n-1+b_i)

f_n＝σ(W_f,X*X_n+W_f,H*H_n-1+W_f,S·S_n-1+b_f)

o_n＝σ(W_o,X*X_n+W_o,H*H_n-1+W_o,S·S_n+b_o)

H_n＝o_n·tanh(S_n)

由于同一个神经单元共享权重，ConvLSTM网络可展开为图3结构所示，随着时间序列中影像的连续输入，长期记忆S_n将不断更新和保存，并由长期记忆单元状态计算出子时间序列中最后一时刻影像的预测标签。网络模型的训练过程实际上是预测标签与真实标签之间定义的损失函数loss，通过反向传播优化模型参数，是模型在训练数据上的损失函数尽可能收敛，以得到可以实现准确分类的网络模型。为了便利计算机对多分类任务的运算，在训练过程中，标签数据都采取了one-hot编码，将不同类别变量进行二进制向量的转变。而在测试过程为例得到最终的预测分类结果，需要通过softmax函数对测试预测值进行解码，得到最终像素级别的分类预测结果

且每一个像素对应一个类别值。

在实际应用中，用于模型训练的长时间序列数据中，获取每张影像所对应的标签数据成本和工作量较大。为了降低标签制作的工作量和成本，基于半监督学习的启发，设计了与原始ConvLSTM网络模型不同的训练函数。原始模型属于监督分类，是利用模型前向传播计算的预测标签与相应的真实标签之间的损失函数，在不断迭代计算过程中优化网络模型参数，使损失函数收敛。而本发明定义的损失函数loss是由监督与无监督两部分组成，如下公式所示。属于监督分类的损失函数是由影像X_i作为最后时相的子时间序列数据的输出z_i，与i时刻影像对应的真实土地利用标签

之间的标准交叉熵。而另一部分的损失函数是由最后时刻影像为无标签的子时间序列数据计算，将同一子时间序列数据输入网络两次，由于第二次输入有随机的时序变化，使得两次隐藏层的输出z_i和

不同，并作平方差。通过设置一个随时间变化的变量w(t)，以权衡两部分损失函数的占比，并获取较为准确的训练模型。w(t)的初始值设为0，即不计算无监督部分的损失值，而随着新数据的输入和不断的迭代计算，w(t)也随之增大，直至损失函数收敛，模型各参数值也趋于稳定，最终完成训练过程得到模型参数。

其中，B表示每次输入的批次量，L代表有标签的训练样本，C表示类别数量，z_i表示时间序列输入所对应的预测标签经softmax函数转化后的输出概率，

表示与同一个时间序列数据输入经过数据增强或时序变化后所对应的预测标签概率。

为了验证本分类方法的可行性，选取了中国黑龙江省佳木斯市中的三块区域作为实验研究区域(尺寸均为256*256像素)，如图5所示，该区域主要以城市扩张和植被物候变化为主体的土地利用类型，分别表示为Data_1、Data_2和Data_3。使用了历史Landsat 8/OLI影像，分辨率为30m，收集了2015～2016年间的无云及部分云遮挡的影像数据(即去掉了全云或多云遮挡的影像)，并进行辐射校正、大气校正、图像配准及裁剪等影像预处理工作。最终，每张影像有7个波段数据(即海岸波段、蓝B、绿G、红R、近红外NIR、两个短波红外波段SWIR)，且三组时间序列数据的时间跨度为两年，分别共有26、38和36张影像，影像的时间分布情况如图6所示，平均每组数据每月均有1～2景影像。对收集到的影像进行辐射校正、大气校正、图像配准及裁剪等影像预处理工作，并根据标准土地利用分类规范和实际地表覆盖情况，确定本实验任务中可分为6类标签(如图4(d))，即耕地、林地、建设用地、水域、云和阴影，从而制作出相应的标签数据集。每组数据集将划分为两个部分，一个是用于模型训练的样本，即除去最后一个时刻的影像不参与训练，其他剩余的所有时序数据都将用于网络模型的训练；另一个是利用包含最后一个时刻的时间序列输入用于测试。

同时，为了验证本发明方法在时间序列影像上的分类效果，将本发明方法与其他一些表现较优的分类方法进行了对比实验，即深度学习的LSTM、ConvLSTM网络及非深度学习的SVM分类器，实验结果如图4所示。并且，使用了三种分类精度评价指标，即总体精度(OA)、Kappa系数(K)、及宏平均F1值(Macro-F1)。如表1所示，将不同分类模型应用于三组不同的长时间序列遥感影像数据集中，结果表明，SVM分类器的分类精度最低，且运算耗时长，而在深度学习方法中，由于卷积结构能有效地处理影像的空间、纹理和光谱特征，使得ConvLSTM分类效果优于LSTM方法(三种分类精度指标均有提升，且OA值可平均达到93.50％以上)。而基于半监督学习的优势，本方法在全标签样本集的分类精度上提升不明显，但在复杂情况下具有突出优势(如缺少已知标签样本，或存在云雾遮挡等情况)，并大量减少了标签标记工作，降低了对长时间序列影像数据集的要求。

表1

为了验证并分析本发明方法在复杂情况下的有效性，还设计了三组不同的实验进行分析和探讨：(1)改变输入数据的时间序列长度，从而探讨时间上下文信息对分类的重要性；(2)通过改变训练标签的数量，从而探究本发明方法在少量标签样本情况下的分类有效性；(3)通过改变训练影像样本中有云遮挡的影像数据占比，从而探讨本发明方法是否能处理含有云雾遮挡所影响的时间序列影像分类问题。三组实验的具体实施过程及实验结果与分析将详细阐述：

实验一：时间序列长度的探讨

时间上下文特征有利于提高遥感影像的分类精度，对于遥感领域等方面的应用非常重要。而针对本方法处理时间序列遥感影像的分类而言，设置时间序列长度为多少才有利于影像分类的精度和效率是一大关键性问题。因此，通过仅改变模型中时间序列长度(参数t)的设置，即图2中经过预训练特征模型后的时间序列数据，由于改变了时间序列长度t从而划分出不同时序长度的子时间序列数据，再流入步骤二的神经网络模型中进行训练。时间序列长度设置为t＝{2,6,10,14,18,20,22}，其他设置保持不变进行对比实验，实验结果如图7所示。实验结果表明，随着时间序列长度的增加，分类精度得到了提升，尤其是在t＝20时表现较优。然而，随着t的值增加至22时，分类效果反而变差，这是因为时间序列总长度是有限的，在时间序列长度增加的同时，时序数据的输入量也随之减小，不利于完全有效地训练模型。同时，时间序列越长，深度神经网络的记忆可能会在不断更新的过程中丢失，就像人类对更长时间以前的记忆更模糊一样。因此，通过本实验可以看出合适的时间序列长度对于时间序列影像的分类而言有着积极的作用，较长的时序长度可以提高分类精度。而在本分类任务中，本发明方法可以充分的利用时序遥感影像的时间、空间和光谱信息，以完成较精准的分类结果。

实验二：训练标签数据量的探讨

由于采用了半监督分类的训练方式，使得本发明方法可以应用于少量标签数据的情况下进行模型的训练。为了探究本发明方法在少量标签训练样本的条件下的分类效果，通过改变真实标签样本在模型训练中的数量，并使用了不同分类方法进行对比实验。

以所有的时间序列影像所对应的标签为基础，去除最后一个时刻影像的剩余数据作为训练数据集，再按时间先后顺序去除一定数量(0、30％、50％、80％、90％)的真实标签数据，而去除的部分用无标签(值为0)代替，最终利用100％的有标签数据、和保留剩余的70％、50％、20％、10％真实标签数据分别在三组研究区域中进行对比实验。由于标签数量为整数值，因此，在保留标签时统一选择了向下取整的方式用于此次实验，而用于此次实验的实际有标签数据的数量如表2统计。实验结果以分类评价指标Kappa系数为例，如图8所示，本方法相比于其他分类器而言，在不同的有标签数量下都能表现出最佳的分类效果，尤其是有标签数据占比仅10％，即有标签数据仅有2～3张影像可用于训练的情况下，依旧有较高的分类精度(OA均为90％以上，Kappa可高达0.9)。实验一已探讨了本分类任务中最佳的时间序列长度t设置为20，由于“多对一”形式的网络模型，输入为20张影像组成的时间序列数据，输出为最后时刻的一张预测标签，而在训练过程中，还将输入最后时刻对应的标签(可能是有标签数据，也有可能是无标签数据)。因此，将总长度为n的影像按时间先后顺序排列后依次以t＝20为单位输入时序数据时，在t＝20这一时刻以前的标签数据并没有得到实际应用。在三组50％有标签数据中分别去除了前13、19和18张的影像，均小于t值的设置，正因为此，有标签数据占总训练标签数据的比例不少于50％的情况下不同的深度学习分类方法也能保持较为稳定不变的分类精度。而随着训练标签数量的减少，相应的分类精度也必然受到一定程度的影响，但本发明方法相比于LSTM和ConvLSTM神经网络而言精度下降幅度较小，且仍旧能保持较优的分类效果。

表2

实验三：探讨云雾遮挡对时间序列影像分类的影响

越来越多的卫星以几天或几小时为间隔可以定期监测地球表面的动态时空过程。然而，因为地球表面通常被云层完全或部分覆盖，卫星影像不可避免地会存在缺失。这限制了大多数遥感方法的广泛研究和应用，并使一些针对无云影像而设计的方法面临着无处不在的挑战。在此次实验中，在原有数据集的基础上，扩充了可收集到的全云或多云覆盖的卫星影像，并统计每张影像中云覆盖量所占的比例后进行归类，以创建此次实验的子数据集，如表3所示，同时，本次实验中训练使用的标签数据都为无云覆盖影像(即云覆盖量小于1％的影像)所对应的标签，仅包含4个类别(耕地、林地、建设用地和水域)。图9展示了本实验结果，可见在所有影像(包含不同程度云覆盖量的影像)和无云遮挡影像的分类精度基本一致，但利用无云影像(云覆盖量<1％)作为训练数据的分类效果最佳，这说明对于时间序列影像分类而言，即使充分利用时间上下文信息，云雾遮挡仍旧会影响分类的精度。另外特别的发现是Data_3数据集在低于25％和低于50％云覆盖量的情况下分类效果较差。这是因为三个区域的数据虽然在同一云覆盖量下的影像数量相近，但Data_3中由于云雾遮挡导致连续数据缺失的情况较为严重，长达5个月没有低云或无云遮挡的影像数据，而其他两组实验区域最多仅连续1～2个月没有低于或无云遮挡的影像数据。由此说明，连续数据的缺失会很大程度的影响了时间序列影像的分类效果。因此，本方法在连续数据缺失但无云遮挡、和有云遮挡但无连续数据缺失的情况下都可以有较为不错的分类效果(OA均在90％以上)，如此一来可以降低对时间序列影像集的要求，使本方法可以更广泛地应用于某季度常有云雾覆盖的区域的遥感研究，如亚热带区域的遥感研究。

表3

Claims (2)

1.一种面向时间序列遥感影像的土地利用分类方法，其特征在于，包括以下步骤：

步骤1，对组成时间序列遥感数据的多光谱影像进行主成分分析，获取三个主成分的影像；

步骤2，对每个三波段影像进行预训练，提取出特征图像；

步骤3，特征图像按时间顺序依次输入半监督卷积长短期记忆网络模型进行训练；

步骤4，利用训练好的模型对最后一时相的影像进行预测分类，获得分类结果；

在步骤1中，所述的时间序列遥感数据集Γ＝(X₁，X₂，X₃，...，X_N)中的所有多光谱影像X_n(n＝1，2，3，...，N)的尺寸为h×w×c，h、w、c分别表示影像的长、宽、和波段数，经过主成分分析后的尺寸为h×w×3，成为三波段影像；

步骤2中所述的预训练过程为：选择在开源的ImageNet大型自然影像数据库上完成训练后的ResNet模型作为预训练模型，提取出预训练模型中浅层的特征图像，尺寸为h′×w′×c′，其中，h′＜h，w′＜w，c′＞c，为了恢复至原始影像大小的长和宽，采用双线性差值算法对提取后的特征图进行上采样，得到尺寸为h×w×c′的特征图像，再将上采样后的特征图像与原始影像的多光谱波段进行融合，以得到尺寸为h×w×(c+c′)的数据；

步骤3中所述的半监督卷积长短期记忆网络模型是基于长短期记忆网络的基础上，将输入到状态和状态到状态之间的全连接传输过程转换为卷积运算，即ConvLSTM单元，在ConvLSTM单元中，随着时间序列数据中每个时刻数据的输入X_n(n＝1，2，3，...，N)，网络记忆单元状态S_n，以及当前隐藏层状态值H_n，由三个门机制控制，即输入门i_n、遗忘门f_n和输出门o_n，并不断被保存和更新，计算公式为：

i_n＝σ(W_i，X*X_n+W_i，H*H_n-1+W_i，S·S_n-1+b_i)

f_n＝σ(W_f，X*X_n+W_f，H*H_n-1+W_f，S·S_n-1+b_f)

o_n＝σ(W_o，X*X_n+W_o，H*H_n-1+W_o，S·S_n+b_o)

H_n＝o_n·tanh(S_n)

其中，S_n-1，H_n-1分别表示上一时刻的记忆单元状态和隐藏层状态，函数σ和tanh分别表示sigmoid和双曲正切函数，W_i，X，W_i，H，W_i，s表示输入门i_n分别与当前的输入X_n、隐藏层H_n和记忆单元状态S_n之间的权重矩阵，W_f，X，W_f，H，W_f，s表示遗忘门f_n分别与当前的输入、隐藏层和记忆单元状态之间的权重矩阵，W_o，X，W_o，H，W_o，S表示输出门o_n分别与当前的输入、隐藏层和记忆单元状态之间的权重矩阵，W_S，X，W_S，H分别表示记忆状态与输入和隐藏层之间的权重矩阵，b_i，b_f，b_o，b_S对应各个门和网络记忆单元的偏置项，

表示当前输入激活的记忆状态，符号‘*’及‘·’分别表示卷积运算和Hadamard矩阵乘积；

在训练过程中，由隐藏层状态计算出当前时间序列数据输入所对应的预测值，并利用已知和未知标签的联合信息计算损失函数，使其收敛以训练并优化网络模型，其损失函数由两部分组成，即利用已知标签信息的标准交叉熵，和联合已知和未知标签信息的无监督损失函数，公式可表示为：

其中，B表示每次输入的批次量，L代表有标签的训练样本，

表示该训练样本所对应的真实标签，z_i表示时间序列输入所对应的预测标签经softmax函数转化后的输出概率，

表示与同一个时间序列数据输入经过数据增强或时序变化后所对应的预测标签概率，C为类别数量，w(t)为一个随时间变化的变量，以权衡两部分损失函数的占比。

2.根据权利要求1所述的土地利用分类方法，其特征在于，在步骤4中，利用训练好的模型计算出待测试影像的预测标签概率，并采用argmax函数取概率最大值的类别作为分类结果。

Images