为弥补遥感图像受硬件制约导致的“时空矛盾”问题，近年来国内外研究学者针对能够获取到的高空间低时间分辨率图像和低空间高时间分辨率图像，逐步发展了几种核心理论不同的时空融合算法，以生成高空间高时间分辨率图像。模型主要分为3类：基于权重滤波的时空融合模型、基于稀疏表达的时空融合模型以及基于混合像元分解的时空融合模型^[1,2]。其中基于权重滤波的时空融合模型原理简单，目前发展最为成熟，典型算法包括STARFM(Spatial and Temporal Adaptive Reflectance Fusion Model)^[3]、ESTARFM(Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model)^[4,5]、STAARCH(Spatial Temporal Adaptive Algorithm for mapping Reflectance Change)^[6]、STNLFFM(A Spatial and Temporal Nonlocal Filter-Based Data Fusion Method)^[7]等。它们在植被物候分析^[8]、水污染监测^[9]、地表蒸散发研究^[10,11,12]等方面得到广泛应用。

基于权重滤波的时空融合算法模型基本流程：以目标像元为中心像元，在邻域范围内设定一个搜索框，并按照一定的搜索规则在搜索框内逐像元选取与中心像元具有光谱和时相变化一致性的像元作为相似像元，根据相似像元与中心像元在光谱维和空间维的差异计算每个相似像元所能提供辅助信息的权重比，最后根据设定的融合规则计算预测时刻目标像元值，对图像中所有像元逐一计算生成预测时刻高空间分辨率图像。从上可以看出，相似像元的选取直接影响融合算法性能，是其中的一个重要环节。

目前在基于权重滤波的融合算法研究中，相似像元选取方法一般采用阈值法，基于阈值法的相似像元选取是在高空间分辨率图像上的搜索框内，逐像元地计算其与中心像元的差异，并给出一个判定条件作为阈值，当差值在阈值范围内时则归为相似像元。 ESTARFM模型假定地物时相变化遵循线性变化，并引入能够表示低空间分辨率混合像元与高空间分辨率纯净像元之间转换关系的转换系数，使得该模型在区域异质性程度较高的应用中能够获取得到高精度的融合结果，得到广泛应用^[13,14,15]。ESTARFM模型相似像元选择是利用图像标准差和给定的图像分类数作为阈值进行判定，其中人为设定的图像分类数对相似像元选取结果影响较大，存在一定的不确定性，分类数的大小直接影响相似像元选取的质量和数量，进而影响ESTARFM算法融合精度和算法运行效率^[16]。而由Cheng等^[7]提出的STNLFFM时空融合模型，阈值是以中心像元值为基础建立判定函数，该方法考虑到了中心像元自身的影响与差异性，减少了人为设定分类数给算法融合性能带来的不确定性。

因此，本文将STNLFFM中相似像元选定规则与ESTARFM模型相结合，对ESTRAFM模型进行改进，提出ESTARFM_NL算法，在保持精度的同时降低了其实际应用过程中的不确定性，并大幅提升算法运行效率，这对于大区域海量数据时空融合的应用具有重要意义^[17]。

ESTARFM算法中要求至少输入两组或多于两组参考时刻t_k(k=1,2,…l)的高空间分辨率图像和低空间分辨率图像，以及预测时刻t_p的低空间分辨率图像。算法主要包括以下3个步骤。

ESTARFM模型中首先针对每一参考时刻进行相似像元选取，进而对其取交集，作为最终相似像元选取的结果。单一参考时刻的相似像元选取计算方式如下：

其中：F(xi,yj,tk,B)和F(xw/2,yw/2,tk,B)分别表示在参考时刻tk高空间分辨率图像搜索框内相似像元(xi,yj)和中心像元(xw/2,yw/2)在第B个波段的像元值，w表示邻域内搜索窗口大小，ω(B)表示高空间整幅图像在第B个波段的标准差，m表示搜索窗口内分类数目通常全图设同一个值。可见，相似像元选取时基于每个参考时刻内的高空间分辨率图像，以目标像元为中心，在邻域范围内设定一定大小的搜索框，判断框内每个像元与中心像元值的差值，若该值小于阈值，则该像元为相似像元。

从公式(1)可见，搜索框大小和分类数目共同影响了相似像元选取结果。基于ESTARFM算法的应用研究中，关于搜索框大小设定，利用Landsat和MODIS反射率数据进行时空融合时，通常将搜索框区域面积设定为1 500 m×1 500 m，覆盖了空间分辨率为500 m 的MODIS图像3×3个像元，对应到空间分辨率为30 m的 Landsat图像数据时搜索框大小为51×51个像元^[4,18]；考虑到空间分辨率的倍数关系，Irina等^[19]在使用时将Landsat数据重采样到25 m，此时对应到Landsat图像上为61×61个像元。Ma等^[16]对ESTARFM算法搜索框大小的影响分析表明，对于25 m的Landsat数据，搜索框适宜设置为61×61个像元，小于该值时，融合结果精度随搜索框增大而提高，但大于该值时搜索框增大时精度上升不再明显。同时，分类数m的大小决定阈值大小，进而影响选取的相似像元结果，分类数一般是基于高空间分辨率图像中的地物类别数进行设置，Zhu等^[4]在ESTARFM算法研究中将m设置为4;Fu等^[18] 提出利用Landsat分类产品来确定每个搜索框内的地物类别数;Ma等^[16]针对不同分类数对ESATRFM融合精度分析表明，当分类数为3时融合精度最高，但计算效率较低。

可见，在相似像元选择中，现阶段对于搜索框大小的设定较为统一，一般都设置成空间范围为1 500 m×1 500 m。而对于分类数，现有分析表明其对相似像元选取有很大影响，但缺乏相对统一的标准，上述分析表明分类数为3时，ESTARFM模型算法结果精度最高。因此，本文基于Landsat与MODIS图像数据进行实验与对比分析时，对于ESTARFM模型设置的搜索框为1 500 m×1 500 m，分类数为3。

对于选取的相似像元，要对其进行权重分配，以用于预测时刻目标像元值的计算。ESTARFM模型中权重计算步骤如下：

(1）计算基于参考时刻高空间分辨率图像的空间距离D_ijk：

上式中D_ijk表示相似像元与中心像元的空间距离，A为常数用于控制空间距离的大小，通常设定值为搜索框宽度^[4]。

(2）计算不同空间分辨率图像上(xi,yj)处相似像元在所有参考时刻和波段间的相关系数R_ij：

其中:E表示均值，G表示方差，F_ij和C_ij分别表示相似像元(xi,yj)在高空间分辨率图像和低空间分辨率图像上各个波段和所有参考时刻的像元值的集合。

（3）相似像元权重计算

综合考虑上述两个因素，并进行归一化处理，形成相似像元的权重Wijk：

其中：Pijk=(1-Rij)×Dijk。可见，相似像元权重计算是考虑了相似像元与中心像元的空间差异性，以及在参考时相上不同空间分辨率图像之间相似像元值的相关性。

首先，针对某一参考时刻t_k的图像数据计算预测时刻t_p的中心像元值，公式为：

其中：C(xi,yj,tk,B)和C(xi,yj,tp,B)分别表示参考时刻t_k和预测时刻t_p低空间分辨率图像搜索框内相似像元(xi,yj)在第B个波段的像元值，V表示转换系数，利用预测时刻前后两组不同空间分辨率图像的相似像元值做线性回归分析，计算得到的斜率即为转换系数，表征了低空间分辨率图像混合像元与高空间分辨率图像纯像元间的转换关系^[4]。

进而，综合多个参考时刻预测结果，得到最终的目标像元值：

其中:

STNLFFM是一种基于非局部滤波的时空融合算法，该算法选取相似像元时阈值计算公式为：

其中：d为一个自由参数，不同传感器数据存在微小差异^[7]，Landsat和MODIS时空融合时算法设定为0.01。

从公式(9)可以看出，阈值判别函数是基于参考时刻t_k高空间分辨率图像中心像元值来计算，在相似像元选取过程中，针对每个中心像元，阈值大小的设定与参考时刻高空间分辨率图像上的中心像元值大小相关，当中心像元值较大时，阈值比较大，而当中心像元值较小时，则对阈值的要求更为严格，体现了选取的相似像元与中心像元值的相对变化性。同时，该计算过程消除了人为设定分类数大小而引起的融合结果不确定性,可见，基于该公式进行阈值选取更为合理。因此，本文针对目前得到广泛应用的ESTARFM模型，将2.1.1节相似像元选取中的阈值计算采用公式(9)进行改进，以增强其在实际应用中的适用性，提出了ESTARFM_NL算法。

选取的实验数据集是遥感图像时空融合研究中的常用的LANDSAT和MODIS数据集^[7,16,19]，该数据集对应新威尔士南部的科利姆巴莱Coleambally（150°48’47”E，34°4’4”S）的灌溉区域（CIA），主要地物覆盖类型有农田、河流、林地和土壤等。实验中用到该数据集的5组图像数据,分别为2001年10月8日(T1)、10月17日(T2)、11月2日(T3)、11月9日(T4)以及11月25日(T5)，图1为各组数据的标准假彩色显示结果。

其中，MODIS数据来源自MOD09GA反射率产品，Landsat数据是利用MODTRAN4进行大气校正得到的反射率数据，数据集涵盖了其对应的7个波段，波段信息详见表1。数据集已完成几何配准、重采样等预处理工作，空间分辨率均已被重采样为25 m，图像覆盖范围43 km×51 km，图像大小为1 720×2 040。

为了对比分析不同数据量下ESTARFM和ESTARFM_NL算法的运行效率和融合精度，本文对实验图像数据集进行裁剪，该区域图像大小为800×800。图1中黑框为裁剪区域。根据2.1.1节的分析，实验中ESTARFM模型的搜索框设置为61×61，对应了1 500 m×1 500 m的搜索区域，分类数设置为3，此时ESTARFM模型的精度最高。ESTARFM模型采用的是作者提供的IDL代码(https://xiaolinzhu.weebly.com/open-source-code.html)，其中ESTARFM_NL由该代码改编而成。

实验设计了以下两组不同地物物候变化规律下的时空融合实验，每一组实验数据均包括整景图像和裁剪后的小区域图像。

实验一：T2和T4作为参考时刻图像组对T3时刻的高空间分辨率图像进行预测，时间跨度相对较小为23 d，该时段内图像时相变化较小。

实验二：T1和T5作为参考时刻图像组对T2时刻的高空间分辨率图像进行预测，研究时间跨度相对较大为48 d，参考时刻与预测时刻之间地表覆盖类型发生了较为明显的植被生长，图像时相变化较大。

利用相关系数（R）和相对误差RE对融合结果精度进行定量评价，其中相关系数R具体计算方式如下：

其中：σP和σO分别表示融合图像和预测时刻Landsat真实观测图像（后文统称观测图像）的方差，σPO表示融合图像和观测图像间的协方差。R值越接近于1精度越高，越接近于0则精度越低。

相对误差RE具体计算方式：

其中：F_P表示预测时刻的融合结果，F_O表示观测图像。在精度评价中展示了观测图像与融合结果间的相对误差图像在全谱段的概率分布曲线，并计算其平均相对误差RE_m。

采用相对差异指标定量化对比两种算法的融合结果精度和运行时间，用于分析两种算法的融合性能。

图2(a)和(d)分别展示了实验一和实验二在预测时刻的Landsat图像，图2(b)、(c)和(e)、(f)分别表示实验一与实验二利用ESTARFM、ESTARFM_NL得到的融合结果图像，两实验组中第一行展示整景图像，第二行展示裁剪后图像，图像均以标准假彩色展示，整景图像融合实验中由于参考时刻图像数据有部分缺失，因此融合结果图像有部分缺失可参考图1中T1、T4、T5的高空间分辨率图像。对比融合结果与观测图像，实验一中两种算法的融合结果均能够充分反映预测时刻的植被物候变化信息，在实验二中两种算法在预测时刻的融合结果与观测图像相比均有光谱畸变发生；对比两种算法的融合结果图像，实验一中观察图2(b)和(c)发现，融合结果在光谱和纹理上一致性较高，实验二中即使两种算法融合结果发生了光谱畸变，但通过对比图2(e)和(f)发现，两者同样表现出了较高程度的光谱和纹理一致性。

在进行精度评价时，对于ESTARFM模型分类数设定为3，它对应了该模型精度最高的情况^[16]。

图3为两组实验的融合结果与观测图像相对误差直方图，其中RE越接近0表征融合结果精度越高。从图中可知：在整景图像融合实验中，两种算法得到的融合结果与观测图像之间的相对误差总体分布情况趋于一致。而对于裁剪后图像融合实验得到的相对误差分布则有所差异，在接近于0的部分，ESTARFM_NL的融合结果更多的分布在0附近，可见对于图像中大多数点而言ESTARFM_NL模型的误差更小。

表 2展示了实验融合结果与观测图像在全谱段范围内的平均相对误差RE_m，以及算法之间的精度相对差异。根据表 2可知，在各组实验两算法的融合结果精度对比中ESTARFM相比于ESTARFM_NL的RE_m较小，ESTARFM融合结果精度较高；对比两算法之间融合精度的相对差异，相对差异最大仅为5.22%，最小为3.81%。

表 3为两种融合算法在实验一与实验二融合结果与观测图像之间相关系数，同时还展示了两种算法融合结果精度间的相对差异。通过对比表 3中各组实验的相关系数可知： ESTARFM_NL融合结果精度始终略低于ESTARFM。从两种融合精度之间的相对差异来看：在裁剪后的图像实验中，实验二绿色波段两种融合算法精度之间的相对差异最高仅为3.45%，在实验一的近红外波段两种算法精度相对差异最低为1.37%；在整景图像中的两组实验，在实验二的绿色波段融合结果精度差异最大仅为2.49%，在实验一的绿色波段精度差异最小为0.69%。

上述分析表明ESTARFM_NL与ESTARFM分类数为3时融合精度相当，由2.1.1节的分析可知，目前分类数的设定存在多种情况，当分类数设置为其他值时，ESTARFM_NL可能精度高于ESTARFM模型，ESTARFM_NL模型改进了ESTARFM中相似像元选取带来的不确定性影响，避免了ESTARFM模型融合结果精度受分类数影响导致的不稳定性。

ESTARFM常用于大区域或长时序的遥感图像时空融合应用，涉及到海量数据处理时，时空融合的运行效率非常重要^[17]。

表4为两种算法的运行时间情况，运行环境为CPU:3.60 GHz，RAM:8 GB。在实验过程中输入了两种数据量的图像：整景图像大小为1 720×2 040（40.16 M），裁剪后的图像为800×800（14.65 M）。

从表4可以看出对于整景图像而言，ESTARFM需要2 h以上，而ESTARFM_NL模型则1个多小时就可以完成融合，而对于裁剪后的图像，其数据量大约减少了5.5倍，ESTARFM运行时间需要20 min左右，而ESTARFM_NL则可以在10 min左右完成。将两种模型方法运行时间对比计算可以发现，ESTARFM_NL的运行效率均提高40%以上。而且通过不同数据量的时间对比发现，当数据量增加时，运行时间不是一种线性变化，而是会表现为更慢，可见对于海量的数据处理应用，ESTARFM的运行效率则会更慢，而ESTARFM_NL在保持精度的前提下，可以实现运行效率的大幅上升，对于海量数据应用具有重要价值。同时，ESTARFM_NL避免了阈值设定的主观性，降低了相似像元选取结果的不确定性，使得其融合结果相对稳定，不受制于分类数的变化。

综合上述实验分析，ESTARFM与ESTAR-FM_NL得到的融合结果，从图像直观表达效果来看，两种算法融合结果之间的在色彩和纹理具有较高的一致性；利用相对误差直方图分析两种算法融合精度时发现，在整景图像的实验中，两种算法融合结果相对误差直方图分布总体分布情况趋近于一致，在裁剪后图像中则ESTARFM_NL表现了较高的融合精度；在相关系数R和平均相对误差RE_m的精度评价中表明，ESTARFM与ESTARFM_NL得到的融合精度相对差异最大分别不超过3.45%和5.22%。以上结论可得出，ESTARFM_NL算法融合结果与ESTARFM能够保持较高程度的一致性；从算法运行时间角度考虑， ESTARFM_NL在融合过程中相比于ESTARFM运行效率能够提升40%以上。

相似像元选取是基于权重滤波算法的重要环节之一，ESTARFM在相似像元选取过程中需要利用图像标准差和给定的分类数进行阈值计算，分类数设定中存在不确定性，直接影响了相似像元选取结果，进而影响算法融合结果精度和运行时间。而在STNLFFM中利用中心像元值计算阈值的方法，考虑了中心像元与相似像元之间的相对变化性，且该方法减弱了人为因素对相似像元选取过程的影响。因此对ESTARFM相似像元选取过程中的阈值判定函数进行改进，提出了ESTARFM_NL算法，消除了相似像元选取过程中的不确定性，提高算法稳定性。

为评价ESTARFM和ESTARFM_NL两种算法的性能，进行了不同时相变化条件的融合实验。对融合结果精度和算法运行效率进行对比分析，对预测时刻真实图像和融合结果进行目视解译，发现两种算法融合结果相似程度较高，且均能反映地物的光谱变化情况。同时利用相对误差直方图、平均相对误差、相关系数对融合结果进行评价发现，ESTARFM融合结果精度较高，同时对比两算法之间融合精度的相对差异，相对差异最大仅为5.22%，最小为0.69%。在运行效率的对比中，ESTARFM_NL算法相比于ESTARFM提升40%以上。

因此，ESTARFM_NL算法融合结果精度能够保持与ESTARFM相当的情况下，大幅缩减算法运行时间，这在处理大区域海量数据或长时间序列的遥感图像数据中有良好的应用前景。