随着卫星遥感技术的发展，多传感器遥感影像为地质灾害分析及预防、水体环境变化监测和评估等方面提供重要的数据支撑^［1］。但由于技术条件限制，兼具高空间分辨率、高时间分辨率等特征的遥感数据通常不易获取，这极大地影响了其在地表信息动态监测研究中的应用^［2］。此外，云覆盖等不利天气条件会遮挡传感器接收下垫面电磁波信号，造成信息缺失问题，降低了遥感影像的可用性。针对这些问题，在不改变现有观测条件的情况下，结合多传感器影像优势特征获得兼具高空间分辨率和高时间分辨率特征的时空融合方法应运而生。

目前时空融合方法主要分为4类：基于像元解混的融合、基于滤波的融合结果、组合型的融合方法和基于学习模型的融合方法^［3］。基于像元解混的方法最早由Zhukov等^［4］提出，该方法引入窗口技术，基于目标像元和窗口内同类地物反射率空间相似性假设进行时空融合，但未能彻底解决同类地物光谱空间变化的问题。此后，众多学者针对反射率异常值^［5］、地物类型变化^［6］、窗口内相邻像素空间可变性^［7］等开展了更深入的研究。总体来说，这类方法因建立在线性光谱混合模型的基础上，而遥感影像因分辨率限制常常存在大量非线性混合像元，在像元分解过程中得到的反射率值误差明显因而不适用于异质性较强的区域。

基于滤波的融合方法中，应用最广泛的是Gao等^［8］提出的时空自适应反射模型（Spatial and Temporal Adaptive Reflectance Fusion Model， STARFM）和Zhu等^［9］提出的增强型时空自适应融合模型（Enhanced Spatial and Temporal Adaptive Reflectance Fusion Mode，ESTARFM），这两种方法的基本思想是基于加权函数构建高—低空间分辨率相关关系，STARFM方法主要建立在低空间分辨率影像中一个像素内只有一种土地覆盖类型这一假设上，因此在异质景观上并不适用从而无法实现准确预测；ESTARFM通过加入转换系数改进相似像元权重计算方式，并在基础时间的两组高分辨率影像上同时搜索移动窗口内的相似像元，有效地提高了模型在异质性区域的融合精度。其他学者在景观异质性^［10-11］和权重函数^［12-13］等方面提出了一系列改进算法。此外，Wang等^［14］综合回归模型拟合、空间滤波和残差补偿的三步法生成哨兵2号卫星时间序列影像，该方法在环境动态变化显著的监测上具有很大的应用潜力。这类方法目前具有较多的研究，在时空融合结果上取得较大成果但在不同时相影像的地物类型变化预测效果上仍有一定的局限性。

组合型的融合方法包含不同融合方法的优势，从高/低分辨率数据处理方式和空间信息变化等多方面因素来综合考虑。如Zhu等^［15］提出灵活时空融合模型（Flexible Spatiotemporal Data Fusion，FSDAF），在地表覆盖异质性明显的区域适应性更高，能够较好的预测土地覆盖类型变化。进一步地，Zhao等^［16］提出鲁棒自适应时空融合模型（Robust Adaptive Spatial and Temporal Fusion Model， RASTFM），通过组合像素级和特征级融合提高模型的准确性。

基于学习模型的融合方法，一般是在高分辨率影像和低分辨率影像对之间进行学习。代表性研究有基于纯字典学习^［17］、稀疏表达与字典学习结合^［18］、压缩感知学习^［19］的时空融合方法。总体上，这类方法虽然在土地覆盖类型变化区域具有良好的处理能力，但是存在深层特征挖掘不足、计算复杂、效率较低等问题。

受到卷积神经网络（Convolutional Neural Network， CNN）能够自动提取大量遥感影像结构和纹理特征信息优势的启发，近年来众多学者将其引入时空融合方法中^［20-24］。如Song等^［20］提出基于深度卷积神经网络的遥感时空融合方法（Spatiotemporal Satellite Image Fusion Using Deep Convolutional Neural Networks， STFDCNN）、Tan等^［21］提出的深度卷积时空融合模型（The Deep Convolutional Spatiotemporal， DCSTFN）、Liu等^［22］提出双流卷积神经网络（A Two-stream CNN for Spatiotemporal Image Fusion， StfNet）、姚凯旋等^［23］提出密集连接块改进了CNN网络、Zhang等^［24］提出了基于生成对抗网络的时空融合模型。这类方法对时空融合中尺度转换和时相变化等复杂非线性变换具有较好的建模能力，因此越来越受到欢迎。

相比已有融合方法，尽管基于深度学习的方法显著提高了复杂异质性景观信息预测能力，但均假定输入数据无云覆盖或进行云掩膜处理。根据时空融合基本原理，基础时间高/低分辨率数据对可以择优晴空条件的高质量数据，但预测时间低分辨率数据往往具有不可替代性。由于无法保证预测时间低空间分辨率数据不受云覆盖影响，现有方法往往不能融合得到所需时间的高分辨率数据^［25-26］。因此，对于预测时间存在云覆盖的低分辨率数据来说，研究面向云覆盖的遥感影像时空融合方法具有重要意义。目前已有大量文献对云覆盖像元缺失值重建进行了研究^［27-28］，特别是深度学习框架下遥感信息重建取得了很大进展^［29-32］。因此本文将遥感时空融合与云下信息重建结合起来处理，提出一种面向云覆盖的遥感影像时空融合深度学习模型。首先，基于基础时间MODIS数据，设计重建子网络恢复预测时刻存在缺失像素的MODIS影像；在此基础上，构建时空融合子网络将重建后的MODIS影像和前后相邻时刻的高、低分辨率数据对进行融合，最终得到所需时间的高分辨率融合结果，进而用于水体信息提取等应用研究。

图1给出了面向云覆盖的遥感影像时空融合过程的示意图，已知连续时刻的低空间分辨率数据序列（如t₁、t₂、t₃的MODIS数据称之为低分辨影像，t₂时刻MODIS影像存在缺失像素值），与之对应时刻下的高空间分辨率数据序列（如t₁、t₃的Landsat数据，称之为高分辨影像）。利用t₁、t₃时刻无云覆盖的MODIS影像及其对应Landsat影像与t₂时刻含云覆盖的MODIS影像可融合出t₂时刻具有类似Landsat分辨率的影像，即时刻t₂的高空间分辨率影像。此融合过程在时间域和空间域同时进行，达到充分利用高—低分辨率影像序列时空变化信息的目的。

实验采用的深度学习模型包含低分辨率影像缺失信息重建子网络和时空融合子网络。在低分辨率影像缺失信息重建子网络中，利用局部残差连接的特性，连接对应尺度的浅层结构信息和高层网络中的语义信息，提高网络训练中信息提取的丰富性；在时空融合子网络中，基于解码和编码结构，利用多时相高—低空间分辨率影像数据组合以达到充分利用时空特征信息的目的，从而提高融合模型的性能。两部分结构均采用上采样和下采样处理进行特征提取和细节特征恢复。该网络总体结构图如图2所示。

重建子网络主要是针对低空间分辨率影像的信息处理，共包括16个卷积层，卷积核从64增加到512个。如图2（a）所示，首先将t₁、t₃时刻的高分辨率影像分别与中间时刻t₂局部存在缺失像素值的低分辨率影像（M_t）合成并相加，再将合成与相加的结果对应组合，使用卷积核大小为3×3的卷积层进行特征提取，卷积运算的数学公式：

其中：Y输出特征图；w、b表示卷积核参数；x是输入图像。

为了使训练过程更快更稳定，使用批量归一化层（Batch Normal， BN）将输入进行归一化，删除冗余参数加速网络训练，BN通常作用于非线性层之前。由于校正线性单元（Rectified Linear Unit， RELU）能够增强模型非线性表达能力并且避免训练过程中出现梯度消失现象，因此选用RELU作为激活函数。此后通过建立大范围内不同像素之间的关系恢复缺失区域的信息，网络中使用平均池化层进一步获取特征信息及背景信息。模型采用局部残差连接充分利用浅层影像的细节特征信息，更好地挖掘低分辨率和高分辨率影像对之间的深层关系，丰富的信息利用提高了模型训练的效果。在上采样中引入转置卷积层，进行端对端地学习恢复高分辨率影像；每个阶段对前一次输出的影像进行上采样，得到更高的分辨率，同时与下采样过程中相同尺寸的特征融合，该过程能有效保留辅助影像的语义信息，最后通过卷积层进行精调得到更为精细化的MODIS影像，以达到预测时刻缺失信息重建的目的。

时空融合子网络的输入影像分别是t₁时刻的精细—粗分辨率影像对（L₁， M₁），t₃时刻的精细—粗分辨率影像对（L₃， M₃）以及重建后的t₂时刻粗分辨率影像（称之t₂过渡影像，H_t）。如图2（b）图所示，为了充分利用Landsat影像的高空间分辨率特性，将H_t与L₁-M₁、L₃-M₃影像对进行线性组合，多次利用卷积层进行特征提取，该部分上采样处理使得卷积核数目从64依次增加到128、256、512，同时影像尺寸对应缩小2倍，最后使用反卷积层逐步恢复影像尺寸，此时卷积核数目从512递减到256、128、64，影像尺寸对应扩大2倍。为利用浅层特征信息，使用3个局部残差连接将收缩和扩张层对应尺寸的特征影像进行融合，得到最终的t₂时刻高空间分辨率影像（L_t）。该子网络包含25个卷积层，每个卷积核大小为3×3。此过程中将连续时刻辅助影像的空间信息、时间信息在网络中始终同步提取，保证了重建与融合过程的信息补充。

为了测试提出方法的应用潜力，以融合结果的水体提取为应用对象，根据研究区域地处平原、无山体阴影干扰、沉陷区周围建筑物较少等特点，采用改进的归一化水体指数模型^［33］（Modified Normalized Difference Water Index， MNDWI）进行水体提取。该方法不仅能够较好地抑制植被的影响突出水体特征，而且在区分建筑物特征和消除阴影方面得到更好的效果。MNDWI的计算公式为：

其中：Green、MIR分别为绿光、中红外波段的反射率值，即Landsat-7 ETM+影像的B2、B5波段，Landsat-8 OLI影像的B3、B6波段。通过与Landsat影像进行目视对比、直方图分析来确定合适的阈值，完成研究区的水体信息提取。

研究区地跨安徽省阜阳东部颍上县及淮南西部凤台县，区域范围界限在116°08′01″—116°44′51″ E，32°41′10″—32°56′05″ N之间。该区域地势平坦，水系分布较多，近些年煤炭资源开采频繁导致该区域分布了较多因矿区地表沉陷而形成的大面积沉陷水体。数据集由2014年12月至2015年10月之间的10对无云Landsat影像和MODIS影像组成，如表1所示。高分辨率影像由陆地卫星影像中的Landsat-7 ETM+/Landsat-8 OLI传感器获得，并使用ENVI补丁进行条带修复和ENVI-FLAASH模块进行大气校正，采用最邻近法重采样到25 m空间分辨率。低分辨率影像由MODIS地表反射率日产品（MOD09GA）产品获得，所有MODIS数据均使用MODIS重投影工具（MODIS Reprojection Tools，MRT）重新投影到WGS84坐标系，并采用最邻近法重采样到25 m空间分辨率，使其与Landsat影像保持一致。实验影像大小为2 400×1 200，选取常用于区分水体边界的4个波段（绿波段、红波段、近红外波段、中红外波段）进行实验。详细的光谱范围如表2所示。

为了便于比较，裁剪具有代表性的子区域（影像大小为700×220个像素）开展定性和定量比较；该区域地表信息中水体斑块包括采煤沉陷区块状水体和细碎水体，以及狭长河流自然水体和其他非水体地物信息。图3给出了2015年1月19日由“中红外、近红外、红色”三波段组合的假彩色影像，三对连续时刻的Landsat［图3（a）—图3（c）］及对应的MODIS［图3（d）—图3（f）］实验数据、不同缺失率的掩膜类型数据［图3（g）—图3（i）］。

训练数据选用与待预测时刻相邻的前后各3个时刻的Landsat和MODIS影像组成高—低分辨率数据对。数据分组的基本规则如表3所示，分为两组实验进行，第一组实验数据，以2015年1月19日作为待预测时刻，第二组实验数据，以2015年5月19日作为待预测时刻。其中两组实验的训练数据选择均遵循如下规定：基础时间数据尽可能接近预测时间数据，同时为保证训练阶段预测时间前后时刻的数据分布均衡，基础时间t₁的数据偏重于选取预测时间t₂之前的影像，基础时间t₃的数据偏重于选取预测时间t₂之后的影像，预测时间t₂数据则选取的是该时刻前后两个时刻的影像。具体影像对分组情况以第一组实验数据选择为例进行说明：选用2014年12月2日、12月18日、2015年1月11日、2月12日Landsat和MODIS影像作为t₁时刻的高—低分辨率数据对，2015年1月11日、2月12日、4月25日、5月19日的Landsat和MODIS影像作为t₃时刻的高低分辨率数据对。2014年12月18日、2015年1月11日、2月12日、4月25日Landsat和未受云污染的MODIS影像分别作为训练阶段中待预测时刻t₂的标签数据及低分辨率数据，为了能够充分训练网络，将大小为2 400×1 200像元的Landsat影像和MODIS影像对应裁剪成160×160像元的影像块，重复裁剪滑动步长为80像元。

为了模拟不同程度像素缺失情况，制作了具有不同像素缺失率的掩膜，训练过程中，对待预测时刻的MODIS影像的160×160像元影像块进行15%和20%两种缺失率的掩膜处理，得到每一时刻近3 000对高低分辨率影像块用于输入网络进行训练。

在已有的研究中^{［21，24］}，多使用均方误差（Mean Square Error，MSE）作为损失函数，本文考虑到云覆盖区域的异常值影响，损失函数使用具有对异常值敏感度高等特点的回归带参损失（huber loss）^［27］，将重建网络输出的t₂低分辨率过渡影像（H_t）与t₂真实观测MODIS影像（M₂）进行损失计算，记为L¯1。融合得到t₂的高分辨率影像（L_t）与实际观测Landsat影像（L₂）的损失计算，记为L¯2。将两个损失函数值相加之和作为整个网络的损失函数，记为L¯。F（·）表示面向云覆盖的遥感影像深度学习时空融合网络，W₁、W₂分别表示训练过程中重建子网络和时空融合子网络的参数，重建的t₂低分辨率过渡影像和最终融合的高分辨率影像表示如式（3）、（4）所示：

使用以下损失函数公式：

其中：δ是超参数，当真实影像与预测结果之间的差值小于δ时，采用均方误差公式，当真实影像与预测结果之间差值大于δ时，采用线性误差公式。式（5）中θ₁表示重建子网络权重和偏置项；式（6）中θ₂表示时空融合子网络权重和偏置项，根据经验将δ设为1.0。式（7）中L¯1与L¯2之和作为整体网络的损失函数L¯。选用自适应矩估计方法（Adaptive Moment， Adam）对两次损失进行优化学习并更新网络权重直到loss值无明显变化，即达到收敛状态。Adam中将一阶矩估计的指数衰减设置为0.9。二阶矩估计的指数衰减设置为0.999。学习率最初设置为0.000 1，每十轮迭代乘以衰减因子0.1，逐步减小学习率。

预测阶段，直接将预测时刻的前后时刻t₁、t₃的高—低分辨率数据对（Landsat和MODIS数据对），预测时刻的模拟存在缺失信息的低分辨率数据（局部区域经过掩膜处理，模拟缺失信息情况）输入到已训练好的网络中，如在第一组实验中，将2015年1月19日作为预测时刻，t₁、t₃分别是2015年1月11日、2月12日的Landsat-7 ETM+和MODIS数据对，第二组实验中将2015年5月19日作为预测时刻，t₁、t₃分别是2015年4月25日、7月30日的Landsat-8 OLI和MODIS数据对，融合得到预测t₂时刻的高分辨率影像。一般来说，影像质量评估方法包括基于人类感知的主观方法和客观计算方法，主观评价上主要通过视觉效果上的纹理细节、颜色等方面评价。客观上，将预测日期的融合结果与真实观测的Landsat影像进行比较，选用以下两个评价指标进行方法性能的定量评价。第一个指标是均方根误差（Root Mean Square Error，RMSE），用来衡量融合结果与真实观测影像之间整体差异，RMSE越小，说明影像质量越好。第二个统计指标是结构相似性（Structural Similarity， SSIM），该指标重在比较融合结果和真实观测影像之间的整体结构的相似性来评估融合结果的空间信息恢复效果。SSIM有效范围在［-1，1］之间，SSIM值越大，说明影像质量越好。

为了有效验证本文方法，对原本无云的预测时间MODIS数据模拟不同面积大小的云覆盖，得到不同缺失率的预测时刻MODIS数据。此外，采用传统经典时空融合方法ESTARFM作为对照，并且在ESTARFM运行前，对不同缺失率的预测时刻MODIS做样条插值（Spline Spatial Interpolation Method，SSIM）处理来重建缺失信息。针对2015年1月19日和5月19日两个预测时间融合结果，表4给出了两种方法在不同缺失率情况下典型区域的定量对比，即两种方法分别在绿波段、红波段、近红外、中红外波段上的综合定量评价结果。表中黑色加粗数字表示最佳结果。由表4可见，本方法融合结果在4个波段上的RMSE值大多数低于ESTARFM预测结果，而SSIM值大多数高于ESTARFM预测结果。相比ESTARFM，不同像素缺失率条件下本方法取得较小的RMSE平均值和较高的SSIM平均值，且变化幅度较小。这表明本方法在纹理细节上重建精度更高，空间信息恢复效果更好，本文所提出的融合方法在地表反射率方面提取结构细节能力明显优于ESTARFM方法，模型稳定性更强。

为了比较本文方法和ESTARFM方法的视觉效果，图4展示了2015年1月19日和5月19日真实Landsat局部影像及两种方法的对应结果。第一、三行显示真实Landsat影像、两种方法融合结果，二、四行对应展示了红色、黄色矩形框中影像局部放大的细节图，如图4（d）—图4（f）和图4（j）—图4（l）。由图4（a）—图4（c）和图4（g）—图4（i）可以直观地看出，两种方法的融合结果均能够预测出接近于真实Landsat影像的地物信息，但仍存在一定程度上的光谱失真；具体地，如图4（f）和图4（l），ESTARFM在部分水体区域，影像纹理信息缺失严重，水体边界模糊，存在明显的噪声斑点、地物信息颜色过亮而失真；而本文方法预测影像的纹理信息更接近真实影像，水体边界信息丰富完整、整体较为平滑，狭长细小水体轮廓清晰，例如图4（d）红色矩形放大细节图和图4（j）黄色矩形放大细节图。因此本文融合方法重建及恢复细节信息能力更强。

为验证两种方法融合结果在典型区域的水体提取应用效果，以2015年1月19日和5月19日预测结果为例，采用MNDWI分别对缺失率33%、59%的本方法预测结果、ESTARFM预测结果和真实Landsat影像进行水体提取，图5展示了MNDWI使用绿波段和中红外波段进行水体信息提取的分类图，其中白色为水体、黑色为非水体。从图5可以直观看出，两种方法的融合结果在水体提取效果上均能够提取出无干扰信息存在的水体，但是ESTARFM方法因部分存在水体的区域融合效果较差，水体提取难度增大从而导致未能完整的识别水体信息，因此提取结果中存在较多漏提和明显的细碎误提水体信息现象。而本文方法预测结果的水体提取效果更加接近真实影像的水体提取效果，水体斑块清晰、能够较为准确地捕捉到水体边界，漏提和误提现象不明显，如图5中白色矩形框中所示。综合以上水体提取结果，可以认为本文方法的融合结果在该区域水体分类的应用准确度更高、提取的水体信息更加接近真实Landsat影像水体分类结果。表5为两种融合方法融合影像的水体分类精度评价结果。从表5中可以得出，在不同缺失率情况下两种融合方法水体总体精度均达到92%以上，这说明使用MNDWI在两种融合影像上均能够很好地提取水体信息的主要轮廓，但是由于细碎水体的提取存在一定差异，出现了漏分和错分。总体而言本文融合方法水体边界范围提取结果更加完整且接近真实影像水体分类图。另外表5中两种方法的Kappa系数表明，本文融合方法在水体提取监测方面的精度优于ESTRAFM方法。

通过对融合结果及水体提取精度分析发现，本文方法能够获得优于ESTARFM方法的融合结果。这是由于在预测时刻MODIS影像存在缺失信息的情况下，本文方法结合了重建子网络和时空融合网络块，在时间域和空间域同时进行优化影像特征信息，达到了充分利用Landsat影像的空间信息和MODIS影像的时间变化信息的效果，这是传统融合方法所不能实现的。此外，本文方法采用与待预测时刻相邻的前后各3个时刻的Landsat和MODIS影像组成的高-低分辨率数据对，训练数据充分，进而能够更好地训练网络。但是作为对照方法的ESTARFM，一方面，尽管使用样条插值处理重建出缺失信息，其重建的准确性仍然较低，影响了后续融合效果。另一方面，该传统融合方法仅仅将相邻时刻的高—低空间分辨率影像对和中间预测时刻的低空间分辨率影像作为输入数据，未能充分获取Landsat影像的空间信息，故而无法有效地实现云覆盖情况下的遥感影像时空融合。在对应融合影像的水体提取方面，本文融合影像纹理细节信息丰富，因而其水体提取效果更加接近真实影像的水体提取效果。但实验仍然存在一定不足：①实验是分波段进行训练及预测，这大大增加了方法的计算成本，且训练数据采用多对高-低分辨率影像组合，在数据量较少的情况下，限制了其应用；而ESTARFM的输入数据量较少，整体效率高于本文方法，因此在今后的研究中，还需要提高本文方法的运行效率。②本文方法融合结果在颜色上稍偏离真实影像，出现一定的光谱失真现象，初步分析是优于训练数据选取的时间跨度过大，可能包含了多个季节不同物候期影像，造成一定的光谱失真。

因此，后续研究还需要进一步改进，尝试使用一对数据结合多个波段的信息进行融合，吸收不同波段的特征信息；同时进一步优化时空融合模型的结构，以此来拓展遥感影像数据信息的可用性。

研究提出了面向云覆盖的遥感影像时空融合深度学习模型，并应用于安徽淮南采煤沉陷区的水体面积提取研究。从预测时刻融合影像的定量对比结果、视觉效果以及水体提取效果中可以得出以下结论：

（1）在预测时刻MODIS影像存在缺失信息时，该方法融合影像更加接近于原始影像，融合结果各波段的RMSE和SSIM均取得较好地定量评价精度，融合影像的细节重建效果明显优于传统经典方法。

（2）由典型区域的水体提取应用效果分析可知，该方法在水体信息提取上取得更好的效果，总体精度为95%，Kappa系数0.83，进一步说明本文方法具有更高的融合精度，降低了时空融合对数据的限制要求，并具有更有效的应用性。