卫星图像同时具有高空间分辨率和高时间分辨率，是研究地表动态的重要数据，在监测土地变化、灾害检测和估算农作物生长态势等方面具有广泛应用(He等，1999)。然而，由于技术限制等原因，目前没有传感器可以同时具有高精度的空间分辨率和频繁的重访周期，限制了那些同时需要高空间分辨率和高时间分辨率数据的应用(Luo等，2000)。例如，Landsat ETM+(landsat enhanced thematic mapper plus)传感器具有精细的空间分辨率(15~60 m)，但重访频率较长(16 d)；相比之下，中分辨率成像光谱仪MODIS(moderate resolution imaging spectroradiometer)传感器的重访频率仅为1 d，但空间分辨率比较粗略(250~1 000 m)。为了获取同时具有高空间分辨率和高时间分辨率的数据，将在时间和空间分辨率上互补的多个传感器卫星图像进行融合。以两个传感器为例，其中传感器Ⅰ(如MODIS和Sentinel-3)具有高的重访频率但是空间分辨率很低，传感器Ⅱ(如Landsat ETM+和Sentinel-2)具有高的空间分辨率但重访频率较低。时空融合的目的是输出具有传感器Ⅰ的时间频率和传感器Ⅱ的空间分辨率的合成图像(Zhu等，2018)。

近年来，时空融合经历了快速发展，在该领域已经诞生了多种方法来融合不同类型的遥感图像。现有的时空数据融合方法主要分为以下3类：基于加权函数的方法(Chen等，2018；Cheng等，2017；Gao等，2006；Wang和Huang，2017；Wang等，2015；Wang等，2017；Zhao等，2018；Zhu等，2010)、基于解混的方法(Amorós-López等，2013；Gevaert和García-Haro，2015；Zhukov等，1999；Zurita-Milla等，2008)和基于字典的学习方法(Huang和Song，2012；Song和Huang，2013)。此外，还有基于贝叶斯的方法，基于神经网络的方法(Song等，2018)等。

在基于加权函数的方法中，Gao等人(2006)提出一种遥感图像时空自适应反射率融合模型STARFM(spatial and temporal adaptive reflectance fusion model)，假设在相同的时间段内，高空间分辨率低时间分辨率图像(以下称为精细图像)的变化与低空间分辨率高时间分辨率图像(以下称为粗略图像)的变化是一致且可比较的，从而可以将从粗略像素导出的时间变化直接添加到精细图像中的像素上以获得预测。基于STARFM模型，Hilker等人(2009)提出了自适应时空融合变化检测方法STAARCH(spatial temporal adaptive algorithm for mapping reflectance change)，能够精确检测土地覆盖类型的变化，但需要多组精细—粗略空间分辨率图像作为输入。Zhu等人(2010)提出了一种增强型的时空自适应反射率融合模型ESTARFM(enhanced spatial and temporal adaptive reflectance fusion model)，认为精细图像间的时间变化与粗略图像间的时间变化存在线性关系，并用线性回归的方法来拟合这种关系，与原始的STARFM模型相比，该方法在异质区域取得了更好的效果。

基于解混的方法中，Zhukov等人(1999)提出的多传感器多分辨率技术MMT(multisensor multiresolution technique)是最早的用解混来获得高时间分辨率和空间分辨率图像的方法，主要通过4个步骤预测精细分辨率图像：1)对输入精细图像进行分类以定义粗略图像的端元；2)计算每个粗略像素的端元丰度；3)在移动窗口内解混预测时间的粗略像素；4)将解混的反射率分配给精细像素。近年来，许多研究者对MMT进行改进，Zurita-Milla等人(2008)在线性分离过程中引入约束，以确保求解的反射率值为正且在适当的范围内。Wu等人(2012)通过加权线性混合方法计算反射率变化，并将估计的变化加到基准日期的精细分辨率图像上以得到融合结果。Amorós-López等人(2013)修改代价函数以避免不切实际的光谱估计同时减少解混误差。Gevaert和García-Haro (2015)直接对粗略像素的变化进行解混，以估计端元的变化。Zhu等人(2016)将解混技术与权重函数相结合，提出一种灵活的时空数据融合算法FSDAF(flexible spatiotemporal data fusion)，在地物发生剧烈变化时仍能取得良好的预测效果。

在基于字典的学习方法中，Huang和Song(2012)最先将自然图像超分辨率的字典学习技术引入时空数据融合，提出了基于稀疏表达的时空反射率融合模型SPSTFM(sparse-representation-based spatio temporal reflectance fusion model)，使用不同观测时间得到的图像间的差分图像块来训练字典，通过稀疏表达建立精细分辨率图像和粗略分辨率图像之间的对应关系，并使用稀疏编码技术重建精细分辨率图像。在SPSTFM之后，Song和Huang(2013)开发了另一种基于单对图像学习的融合方法，仅使用一对精细—粗略图像来训练字典，在输入的精细和粗略分辨率图像对上训练耦合字典，然后通过稀疏编码技术在目标日期重建精细分辨率图像。由于MODIS图像和Landsat图像之间巨大的分辨率差异，该方法使用了双层架构，即首先预测精细分辨率和粗略分辨率之间的中等分辨率图像，然后基于中等分辨率图像预测精细分辨率图像。

在目前的时空融合领域，基于权重函数的融合框架是最主流的方法，该类方法沿用了STARFM算法中的框架。但是在只有一对精细—粗略图像作为先验的条件下，只能通过不同观测时间粗略图像间的差分图像来表达图像的时间变化。然而由于粗略图像的低分辨率，这种方法存在严重的模块效应，并不能准确表达地物的时间变化。因此在地物发生剧烈变化时，难以取得令人满意的结果。针对这种情况，本文提出一种新的遥感图像时空融合方法，主要有以下特点：1)使用线性关系表示遥感图像之间的时间模型，根据时间变化的特性，从局部和全局两个方面对模型进行约束，使得求解的模型能够更准确地表示图像间的时间变化。2)提出一种多时相的相似像素搜寻策略，在搜寻相似像素时更加灵活，能够克服以往的时空融合算法在搜寻相似像素时存在的模块效应，从而提高预测结果的准确性。

在时空数据融合中，不同时间观测到的图像之间的关系称为时间模型。本文假设时间模型与传感器无关，且可以用线性模型来表示，即

$ C_{2}(x, y, B)=a(x, y, B) \times C_{1}(x, y, B)+b(x, y, B) $

(1)

式中，$C_{1}(x, y, B)$和$C_{2}(x, y, B)$分别表示$t_{1}$和$t_{2}$时刻坐标$(x, y)$的粗略像素点，$B$是指定的波段，$a$和$b$是对应的线性模型系数，$a$称为增益系数，$b$称为偏差。

忽略大气校正和成像角度的差异，$t_{1}$与$t_{2}$时刻精细图像之间同样满足上述关系，即

$ F_{2}(x, y, B)=a(x, y, B) \times F_{1}(x, y, B)+b(x, y, B) $

(2)

式中，$F$表示精细图像像素。

就本文的时空融合问题而言，已知观测的粗略图像$\boldsymbol{C}_{1}$、$\boldsymbol{C}_{2}$和精细图像$\boldsymbol{F}_{1}$，目标是求解精细图像$\boldsymbol{F}_{2}$。依据观测数据$\boldsymbol{C}_{1}$和$\boldsymbol{C}_{2}$，估计式(1)中的线性模型系数$a$和$b$，然后通过式(2)求解出精细图像$\boldsymbol{F}_{2}$，将时空融合转化为估计时间模型参数。

依据1.1节对时间模型的分析，将时空融合问题转化为估计模型参数，分析模型参数$a$和$b$的特性。

同一场景不同时刻拍摄的图像发生变化可能有以下两种原因：1)由于拍摄时的大气状态、拍摄角度等条件的差异导致的改变，这种改变可以认为是全局一致的，且变化相对比较平缓；2)由于地物目标发生变化导致的改变，这种改变是局部性的，且改变相对比较剧烈。

基于以上两个原因，模型参数$a$有以下3个特点：1)就全局而言，$a$是关于空间的变量，不同位置的$a$可能会有较大的差异；2)就局部块而言，$a$满足局部一致性；3)在地物没有发生变化的区域，$a$在局部是一致的；在地物发生变化的区域，$a$也是局部一致的；在变化区域与非变化区域邻接的边缘，$a$会有较剧烈的变化。

对于模型参数$b$，其含义是模型残差项，用来保证时间模型的准确性是由于准确捕获了时间变化规律，而不是通过巨大的偏差实现的。因此，$b$应该在0周围变化。

图 1展示了本文方法的流程图。大体来说，本文方法主要分为3个部分：搜寻相似像素，参数估计，精细图像重建。

1.2节中的分析是基于局部一致性前提下进行的，然而在异质区域进行参数估计时，并不能完全满足局部一致性，尤其随着局部块增大，违背局部一致性的可能性加大。为了准确地进行参数估计，需要在局部区域内搜寻相似像素，确保用来做参数估计的数据满足局部一致性假设。另外，考虑到$t_{1}$和$t_{2}$时刻地物目标可能会发生变化，两个时刻相似像素空间分布也可能会发生变化。

基于以上分析，提出一种多时相的相似像素搜寻策略，只有在基础时间和预测时间都与中心像素有着相似光谱信息的像素，才会被认为是中心像素的相似像素。由于$t_{2}$时刻可用的光谱信息都在粗略图像中，如果直接在粗略图像中搜寻相似像素，可能会产生严重的模块效应，因此需要对粗略图像做上采样处理。考虑到粗略图像与精细图像之间巨大的分辨率差异，如果直接将粗略图像上采样到精细分辨率，可能会引入巨大的采样误差，因此将粗略图像$\boldsymbol{C}_{1}$与$\boldsymbol{C}_{2}$上采样到介于精细分辨率与粗略分辨率之间的中等分辨率，得到中等分辨率图像$\boldsymbol{M}_{1}$与$\boldsymbol{M}_{2}$。然而上采样得到的图像缺少高频信息，无法提供充足的光谱信息，因此引入了高通调制操作，具体为

$ \boldsymbol{F M}_{2}=\boldsymbol{M}_{2}+\frac{\boldsymbol{M}_{2}}{\boldsymbol{M}_{1}}\left(\boldsymbol{F}_{1}-\boldsymbol{M}_{1}\right) $

(3)

式中，$\boldsymbol{FM}_{2}$是将$\boldsymbol{F}_{1}$中的高频信息使用$\boldsymbol{M}_{1}$与$\boldsymbol{M}_{2}$的转换系数调制后注入到$\boldsymbol{M}_{2}$后生成的图像。通过式(3)的高通调制操作，$\boldsymbol{FM}_{2}$可以提供更丰富的光谱信息。使用ISODATA(iterative self-organizing data analysis techniques algorithm)算法分别对$\boldsymbol{F}_{1}$和$\boldsymbol{FM}_{2}$做分类处理。值得注意的是，本文认为中心像素的相邻相似像素应该在多个时相与中心像素拥有相似的光谱特性。因此基于以下3个条件来寻找相似像素：1)该像素处于中心像素邻域内；2)在$\boldsymbol{F}_{1}$分类结果中与中心像素属于同一类别；3)在$\boldsymbol{FM}_{2}$分类结果中与中心像素属于同一类别。

图 2展示了寻找相似像素的流程。图 2(a)(b)分别是在$\boldsymbol{F}_{1}$和$\boldsymbol{FM}_{2}$分类图中寻找光谱相似像素的结果，其中黑色像素表示只在某一个时刻被认定为光谱相似的像素，橙色像素表示在多个时刻均被认定为光谱相似的像素，即最终要寻找的像素，图 2(c)展示了最终的相似像素。

在完成相似像素的搜寻后，在局部区域内使用不同时间对应的相似像素进行参数估计。

基于之前对于时间变化的分析，对参数$a$的约束为

$ \min \sum\limits_{i}\left[a \times C_{1}\left(x_{i}, y_{i}\right)+b-\bar{r}\right]^{2} $

(4)

式中，$C_{1}\left(x_{i}, y_{i}\right)$和$C_{2}\left(x_{i}, y_{i}\right)$)分别表示$t_{1}$和$t_{2}$时刻的第$i$个像素，($x_{i}$, $y_{i}$)表示第$i$个像素的坐标，$\bar{r}=\frac{1}{N} \sum\limits_{i} C_{2}\left(x_{i}, y_{i}\right)$表示$t_{2}$时刻第$i$个像素对应的局部均值，$N$表示局部块像素的数量。在满足图像局部一致性假设的情况下，当局部地物没有发生改变时，$t_{1}$与$t_{2}$时刻对应的局部均值也相似，因此为了满足式(4)最小化，$a$的值会在1周围变化；当局部地物发生变化时，$t_{1}$与$t_{2}$时刻对应的局部均值会发生较大变化，为了满足式(4)最小化，$a$值会有较大幅度的变化。

根据之前的分析，希望$b$的影响尽可能小。因此，对$b$的约束为

$ \min |b|^{2} $

(5)

综上所述，通过优化代价函数

$ \begin{array}{c} \min \frac{1}{2} \sum\limits_{i}\left[a \times C_{1}\left(x_{i}, y_{i}\right)+b-C_{2}\left(x_{i}, y_{i}\right)\right]^{2}+ \\ \frac{\gamma}{2 N} \sum\limits_{i}\left[a \times C_{1}\left(x_{i}, y_{i}\right)+b-\bar{r}\right]^{2}+\frac{\beta}{2}|b|^{2} \end{array} $

(6)

来估计时间模型参数。式(6)中，$γ$和$β$代表约束项权重系数。值得注意的是，以上是从局部角度出发进行参数估计，当参与估计的像素存在噪声时，求解的参数很可能会引入误差，不能很准确地表示图像的时间变化。根据之前的分析，参数$a$在同一性质的区域内部(变化区域或者非变化区域)是局部一致的，而在两种区域邻接的边缘则可能会发生突变，因此引入L0范数对参数$a$进行全局优化(Xu等，2011)，具体为

$ \begin{array}{c} \min \sum\limits_{i}\left(a^{*}\left(x_{i}, y_{i}\right)-a\left(x_{i}, y_{i}\right)\right)^{2}+ \\ \lambda \cdot c\left(\partial_{x} a^{*}, \partial_{y} a^{*}\right)+\delta \cdot \sum\limits_{i}\left(\left(\partial_{x} a^{*}\left(x_{i}, y_{i}\right)-\right.\right. \\ \left.\left.h\left(x_{i}, y_{i}\right)\right)^{2}+\left(\partial_{y} a^{*}\left(x_{i}, y_{i}\right)-v\left(x_{i}, y_{i}\right)\right)^{2}\right) \end{array} $

(7)

式中，$c\left(\partial_{x} a^{*}, \partial_{y} a^{*}\right)=\#\left\{i|| \partial_{x} a^{*}\left(x_{i}, y_{i}\right) |+\right.$$\left.\left|\partial_{y} a^{*}\left(x_{i}, y_{i}\right)\right| \neq 0\right\}, \left(x_{i}, y_{i}\right)$表示参与计算像素的坐标，由于是进行全局优化，因此所有像素均会参与计算。$a$表示之前估计得到的增益系数，$a^*$表示优化后的增益系数，$\partial_{x}$与$\partial_{y}$分别表示水平和竖直方向的梯度算子，$h$与$v$是引入的辅助变量，分别对应$\partial_{x}$与$\partial_{y}$。式(7)的求解过程参见Xu等人(2011)的方法。

相似像素具有近似的光谱特性和时间变化，引入相似像素的光谱信息，可以让预测模型具有更好的鲁棒性。因此引入权重函数修改式(2)，重建预测日期的精细图像，具体为

$ \begin{array}{c} F_{2}\left(x_{w / 2}, y_{w / 2}, B\right)= \\ \sum\limits_{i} W_{i} \times\left(a\left(x_{i}, y_{i}, B\right) * F_{1}\left(x_{i}, y_{i}, B\right)+b\left(x_{i}, y_{i}, B\right)\right) \end{array} $

(8)

式中，$w$表示滑动窗口，$\left(x_{w / 2}, y_{w / 2}\right)$表示滑动窗口的中心像素，$W_{i}$表示第$i$个相似像素的权重。为了不引入可能会造成误差的信息，只使用2.1节寻找到的相邻相似像素来计算式(8)。$W_{i}$表明每个相邻相似像素对中心像素反射率估计的贡献程度，考虑到应在多方面做出约束以确保算法的鲁棒性，$W$由以下因素确定：

1) 相似像素与中心像素的位置距离，可以使用像素的几何坐标来计算它们之间的几何距离，具体为

$ {dis}\left(x_{i}, y_{i}\right)=\sqrt{\left(x_{w / 2}-x_{i}\right)^{2}+\left(y_{w / 2}-y_{i}\right)^{2}} $

(9)

为了避免几何距离差异过大降低算法的鲁棒性，将几何距离转化为相对几何距离，具体为

$ D\left(x_{i}, y_{i}\right)=1+{dis}\left(x_{i}, y_{i}\right) / R $

(10)

式中，$R$表示滑动窗口的半径。式(10)可以将像素间的几何距离约束在$[1, 1+\sqrt{2}]$的范围内，避免了过大的几何距离差异可能带来的不稳定。

2)$t_{1}$时刻相似像素与中心像素的光谱差异，可以表示为

$ S_{t_{1}}\left(x_{i}, y_{i}, B\right)=\left|F_{1}\left(x_{w / 2}, y_{w / 2}, B\right)-F_{1}\left(x_{i}, y_{i}, B\right)\right| $

(11)

由此可以判断第$i$个相似像素与中心像素在$t_{1}$时刻的光谱相似性，差异更小的像素有更高的相似性。

3)$t_{2}$时刻相似像素与中心像素的光谱差异，由于没有$t_{2}$时刻的精细图像，因此只能从$FM_{2}$中计算$t_{2}$时刻的光谱差异，具体为

$ \begin{array}{c} S_{t_{2}}\left(x_{i}, y_{i}, B\right)= \\ \left|F M_{2}\left(x_{w / 2}, y_{w / 2}, B\right)-F M_{2}\left(x_{i}, y_{i}, B\right)\right| \end{array} $

(12)

由此可以判断第$i$个相似像素与中心像素在$t_{2}$时刻的光谱相似性，差异更小的像素有更高的相似性。

通过以上3个因素，可以在空间几何距离以及两个时间点的光谱相似性上做出有效约束，使得相似性更高的像素得到更高的权重。综合3个因子，可以计算相似像素的权重，即

$ \begin{array}{c} Q\left(x_{i}, y_{i}, B\right)=\\ \exp \left(-\frac{D\left(x_{i}, y_{i}\right)+S_{t_{1}}\left(x_{i}, y_{i}, B\right)+S_{t_{2}}\left(x_{i}, y_{i}, B\right)}{h^{2}}\right) \end{array} $

(13)

式中，$h$是衰减因子，对式(13)的计算结果进行归一化，得到最终的相似像素权重，即

$ W\left(x_{i}, y_{i}, B\right)=\frac{Q\left(x_{i}, y_{i}, B\right)}{\sum\limits_{i=1} Q\left(x_{i}, y_{i}, B\right)} $

(14)

本文提到的参数，可以根据处理图像的特性做出调整。在实验中，为了将权重函数约束在一个合理区间内，设置衰减因子$h=\sqrt{2}$，且在式(14)中进行归一化处理，可以适用于大多数情况。线性参数拟合正则项权重系数$γ$过大会导致求解的参数整体过于平滑，过小则无法约束参数的局部特性。参数$β$过大可能会导致过拟合，过小则导致求解的参数存在巨大的误差。本文设置$γ=1$、$\beta=\sqrt{2}$。全局优化参数$λ$如果太大，会导致全局优化后的$a$全部游荡在某一个值周围，无法准确表示图像的时间变化，如果设置太小，则优化效果不明显，也无法准确表示图像的时间变化，在前两组实验和后两组实验中分别设置为0.02和0.01。$δ$控制优化迭代的次数，初始值和最大值分别设置为2$λ$和10⁵。滑动窗口大小可以根据所处理图像的特点来设置，在像素光谱差异比较小的同质区域，可以适当缩小滑动窗口以加快计算速度。反之，在像素光谱差异较大的异质区域，则需要扩大滑动窗口寻找足够多的相似像素。在本文实验中，为了找到足够的相似像素以保证预测结果，设置滑动窗口大小为51×51像素。

为了避免引入大气校正和图像配准带来的误差，本文中的粗略图像是由精细图像聚合而成的, 该方法在时空融合研究中得到广泛应用(Gevaert和García-Haro，2015；Zhu等，2016；Wang和Atkinson，2018)。为方便计算，实验中的所有图像均映射为精细图像的大小，即这些图像都具有相同的空间分辨率和覆盖范围。

为了展示本文算法的特性，选取4个场景在不同日期的精细—粗略图像对进行实验。

前两组实验使用Landsat ETM+卫星图像和以此聚合而成的类MODIS图像，精细图像空间分辨率为30 m，粗略图像空间分辨率为480 m，聚合尺度系数为16。场景1是2001年5月24日和2001年7月11日在54°N、104°W地区获取的$t_{1}$和$t_{2}$时刻图像，在研究时段内该区域主要发生了地表事物的物候变化，如图 3所示。场景2是2004年11月26日和2004年12月28日在澳大利亚新南威尔士州获取的$t_{1}$和$t_{2}$时刻图像，在这段时间内该区域发生了剧烈的土地覆盖类型变化，如图 4所示。

后两组实验使用Sentinel-2卫星图像和以此聚合而成的类Sentinel-3图像。与前两组实验所用卫星图像相比，后两组实验所用图像的空间分辨率差异更大，Sentinel-2图像空间分辨率为10 m，类Sentinel-3图像空间分辨率为300 m，聚合尺度系数为30。场景3为2019年4月17日和2019年6月11日在安徽庐江获取的$t_{1}$和$t_{2}$时刻的图像，该地区在研究时段内发生了一定程度的物候和土地覆盖类型变化，如图 5所示。场景4为2019年4月8日和2019年6月27日在湖南省洞庭湖周边获取的$t_{1}$和$t_{2}$时刻图像，该地区在研究时段内发生了大幅的物候变化，且在局部地区发生了剧烈的土地覆盖类型变化，对当前的时空融合算法是一个巨大的挑战，如图 6所示。

选取STARFM(spatial and temporal adaptive reflectance fusion model)算法(Gao等，2006)和FSDAF(flexible spatiotemporal data fusion)算法(Zhu等，2016)与本文算法进行对比实验，比较相关系数CC(correlation coefficient)、峰值信噪比PSNR(peak signal-to-noise ratio)、均方根误差RMSE(root mean squared error)、结构相似性SSIM(structural similarity)、相对全局误差ERGAS(erreur relative globale adimensionnelle de synthèse)和光谱角映射SAM(spectral angle mapper)等6个评价指标。从实验结果来看，本文在物候和土地覆盖类型变化时都能取得良好效果。

场景1主要发生了物候变化，没有剧烈的土地覆盖类型变化，STARFM、FSDAF和本文算法均能取得良好的预测效果，图 7和表 1分别展示了各种算法的视觉比较和量化指标，图 7第2行是第1行红色框内的内容。可以看出，本文算法取得了更好的视觉效果，与STARFM和FSDAF算法相比，本文算法具有更好的评价指标。

场景2的地物在两个日期之间发生了剧烈的土地覆盖类型变化，图 8展示了3种算法在土地覆盖类型发生剧烈变化时的预测效果，第2行是第1行红色框内的内容。从图 8可以看出，STARFM与FSDAF算法在土地覆盖类型发生变化时都没有取得良好的预测效果。由于STARFM算法将从粗略像素导出的时间变化直接添加到精细图像中，因此在预测结果中出现了明显的模块效应。FSDAF算法虽然引入了$t_{2}$时刻粗略图像中的空间信息，但在估计预测结果时只使用了$t_{1}$时刻的光谱分类结果，所以在发生剧烈变化的区域，仍然有比较大的预测误差。本文方法结合了$t_{1}$和$t_{2}$时刻的光谱信息，因此能更准确地预测土地覆盖类型发生的变化。表 2展示了各算法在场景2预测结果的量化指标，可以看出，本文算法在各个量化指标上都取得了更好的效果。

场景3和场景4由于粗略图像和精细图像之间的巨大空间分辨率差异，且两个场景都存在一定程度的土地覆盖类型变化，因此与前两组实验相比，更难得到良好的预测结果。图 9和图 10分别展示了场景3和场景4的视觉比较，第2行是第1行中黄色方框的内容。从结果来看，3种算法均与观测到的真实图像存在较大差异，STARFM算法出现了明显的模块效应，FSDAF算法在预测土地覆盖类型变化时出现了明显误差，相比之下，本文方法更接近真实图像。表 3和表 4是各算法在场景3和场景4预测结果的量化指标。可以看出，与前两组实验相比，后两组实验各个算法在各项指标上均出现了大幅下滑，这是由于更加巨大的空间分辨率差异造成的。不过，与STARFM和FSDAF算法相比，本文方法在大部分指标上均有不同程度的提升。

本文提出一种基于线性模型的遥感图像时空融合算法。通过使用线性模型来表示不同观测时间之间图像的变化，并根据时间变化的特性，从局部和全局两方面对时间模型进行约束，求解出的模型能够准确表示图像之间的时间变化。此外，本文算法使用一种多时相的相似像素搜寻策略，能够更灵活地搜寻相似像素，消除以往算法在搜寻相似像素时存在的模块效应，充分利用相邻相似像素中的光谱信息，提升了预测结果的准确度。从实验结果来看，与传统时空融合算法相比，在视觉比较方面，本文算法能够更准确预测地物发生的变化，预测结果更接近真实图像；在定量评价方面，本文方法在各波段的相关系数、峰值信噪比、均方根误差、结构相似性等指标上均有不同程度的提高。

然而，本文算法仍然存在局限性。为了引入多时相图像的光谱信息，需要对多个时相的图像分别做分类处理，虽然提高了融合的准确性和鲁棒性，但增加了计算的复杂性，降低了算法效率。此外，对$t_{2}$时刻中等分辨率图像$\boldsymbol{FM}_{2}$进行分类处理仍然可能会因为分类不够精确而引入误差。因此，将来的研究工作应尝试着从以上方面继续改进，以期得到更好的效果。