张爱竹，王 伟，郑雄伟，姚延娟，孙根云，辛 蕾，王 宁，胡 光

(1.中国石油大学(华东)海洋与空间信息学院,青岛 266580；2.青岛海洋国家实验室海洋矿产资源评价与探测技术功能实验室,青岛 266237；3.中国自然资源航空物探遥感中心，北京 100083；4.环境保护部卫星环境应用中心,北京 100094；5.国家海洋局北海预报中心，青岛 266061；6.武汉中测晟图遥感技术有限公司，武汉 430223)

0 引言

长时间序列的高空间分辨率遥感数据是高精度地表快速监测的基础，广泛应用于土地利用制图和变化检测、农作物的生长监测、地表温度监控等领域[1]。然而，遥感数据在空间分辨率和时间分辨率上的相互制约，导致单一传感器无法同时获取高空间和高时间分辨率的遥感数据。时空融合技术可较好地解决这一矛盾，在植被动态监测[2]、叶面积指数[3]、城市热岛监测[4]、地表温度生成[5]、土地覆盖分类[6]、病毒传播监测[7]等多个领域得到应用。

在众多已有的遥感时空融合算法中，简单高效的线性模型得到了最为广泛的应用[8-13]。其本质思想是通过构建高低空间分辨率图像中对应像元反射率之间的线性关系，将低空间分辨率图像的时间变化信息融入到已有的高空间分辨率图像中[14]。Gao等[15]最早提出了一种时空自适应反射率融合模型(spatial and temporal adaptive reflectance fusion model,STARFM)，该模型通过融合MODIS图像和Landsat图像生成了具有MODIS时间分辨率和Landsat空间分辨率的地表反射率数据。虽然STARFM模型能够获取预测数据，但是不适用于反射率变化较大的区域和异质性较高的区域[16]。针对前者，Hilker等[17]提出了一种反射率变化映射的时空自适应算法，该算法首先确定不同时相Landsat的空间变化和MODIS的时间变化，然后选择地表覆盖变化与MODIS最接近的Landsat图像作为数据预测的基准数据，从而提高了Landsat图像的预测精度。针对高异质性问题，Zhu等[18]提出了增强型时空自适应反射率融合模型(enhanced spatial and temporal adaptive reflectance fusion model,ESTARFM)，假设地物端元在短时间内只发生线性变换，通过引入转换系数计算出不同端元的反射率变化，提高遥感图像异质性区域的预测精度。虽然这些改进模型取得了一定的应用效果，但是它们的普遍问题是不能适用于土地覆盖类型发生突变区域的预测[19]。

近些年来，为了进一步提高时空数据融合的预测精度，研究人员尝试将短时间内同一研究区的数据变化关系建模为线性变化模型[20]。其本质思想是将基准和预测时刻的两幅低空间分辨率图像的线性方程的系数直接应用于高空间分辨率图像的预测。例如，Cheng等[21]提出了一种基于非局部滤波的时空融合模型，利用两个回归系数更准确地描述地表覆盖变化信息，提高了对复杂变化地物的预测能力；Wang等[22]提出了一种名为Fit-FC的新型时空融合算法，综合利用回归模型、空间滤波和残差补偿提高预测精度。但是，如何对土地覆盖类型发生突变的区域进行数据融合与预测仍然是个难题。

针对时空融合模型在突变区域预测精度不高的问题，本文提出了一种基于分层策略的时空融合模型(hierarchical spatial-temporal fusion model,H-STFM)。通过分别捕捉渐变和突变的土地覆盖类型变化区域，并设计不同的融合方法，有效提高预测图像异质性区域的预测精度。该模型由线性回归模块和加权滤波模块组成，分别用于预测物候变化和突变区域。首先判别目标像素(待预测像素)的变化性质，然后对物候变化像素进行线性回归预测；对突变像素进行加权滤波处理，在超像素邻域中选择变化同质性像素作为相似像素进行加权滤波预测；最后将物候变化和突变区域的预测结果利用优化的时间加权函数融合生成最后预测图像。旨在充分考虑地表覆盖的突出变化，解决时空融合模型在地表突出变化上预测效果不好的问题，使得预测图像更加接近真实地表数据。

1 H-STFM模型

H-STFM模型主要包括4个步骤：变化像素的判断、加权滤波预测、线性回归预测、时间加权融合，具体流程如图1所示。为了便于讨论，本文将低空间分辨率且重访周期短的遥感图像称为“低分辨率图像”，将高空间分辨率但重访周期长的遥感图像称为“高分辨率图像”。

图1 算法流程图Fig.1 Algorithm flowchart

1.1 判断变化像素

H-STFM模型对变化像素进行判断，需要两个基准时刻(图1所示t1和t3时刻)的两幅高、低分辨率图像对和预测时刻(图1所示t2时刻)的低分辨率图像。首先，对上述5幅遥感图像进行数据预处理，将3幅低分辨率图像进行重采样和剪切处理，使其与高分辨率图像具有相同的空间分辨率和分布范围。然后，计算t1和t3时刻低分辨率图像分别与t2时刻低分辨率图像的反射率差值并进行阈值判断，将大于阈值的数据选出作为突变数据，其他作为物候变化数据。以t1与t2时刻低分辨率数据为例，判断方法如式(1)所示：

(1)

式中：M为已知的低分辨率图像；(x,y)为目标像素的坐标；B为图像的波段；σ(M2-M1)为两期低分辨率图像反射率差值的标准差；N为土地覆盖类型的数量。即阈值由反射率差值的标准差与土地覆盖类型数量决定。基于以上方法，可以分别得到t1与t2时刻以及t3与t2时刻低分辨率图像的突变数据和物候变化数据。

1.2 基于线性回归的预测

本文采取两层处理的策略分别对物候变化数据和突变数据进行预测，即对整体物候变化和局部突变分别进行预测。对于物候变化数据，本文采用全局线性回归模型，即对两幅低分辨率图像使用全局窗口来对每个波段中像素进行线性回归拟合，如式(2)所示：

M(x,y,t2,B)=a×M(x,y,t1,B)+b+R，

(2)

式中：a和b为线性回归的两个系数，可用最小二乘估计法得到[10]；R为系统残差。

线性回归模型的理论基础是在同一时刻的遥感图像中，低分辨率图像和高分辨率图像在相同时间范围内的变化是一致的。因此，可以直接将基于t1时刻低分辨率图像构建的回归模型，应用于t2时刻高分辨率图像的预测中。进而得到物候变化区域t2时刻的高分辨率图像，如式(3)所示：

LR(x,y,t2-1,B)=a×L(x,y,t1,B)+b+R，

(3)

式中：LR(x,y,t2-1,B)为基于回归模型基于t1时刻预测的t2时刻高分辨率图像；L(x,y,t1,B)为t1时刻的高分辨率图像。同理，可以基于t3时刻图像预测得到t2时刻的高分辨率图像LR(x,y,t2-3,B)。

1.3 基于加权滤波的预测

对于突变数据采用加权滤波的方法进行预测。经典时空融合算法仅基于光谱信息在矩形窗口中选择相似性像素，不能充分利用高分辨率图像中丰富的空间特征[23]，不适用于突变数据选择相似性像素。为充分利用图像的空-谱信息，本文先对t1时刻的高分辨率图像利用分形网络分割算法进行超像素分割[24]，然后将超像素作为邻域窗口进一步约束相似像素的选择范围。此外，由于目标像素是突变数据，即在该时间段内地表覆盖发生了较大的变化。为了更准确地选择相似像素，本文提出基于混合像元分解的相似像素选择方法：先计算t1与t2时刻低分辨率图像目标像素的光谱反射率差值D，根据高分辨率图像的标准差和土地覆盖类别设置阈值Tp，得到阈值范围为[D-Tp,D+Tp]，如式(4)所示：

Tp=σ(L1)/2×N，

(4)

式中：Tp为判别异常像元的阈值;σ(L1)为t1时刻高分辨率图像反射率的标准差。

然后计算t2时刻高分辨率图像目标像素与超像素内其他像素的反射率差值，如果某一像素对应的差值在阈值范围内，则被选定为备选相似像素；接着利用基于纯净像元指数的端元提取方法[25]与完全约束最小二乘法[26]完成基准时刻高分辨率图像的混合像元分解，得到其端元丰度图；最后，根据实验经验设置端元的变化阈值，将备选相似像素与目标像素反射率差值大于阈值的像素作为最终的相似像素，相似像素个数标记为N(x,v)。该方法基于地表覆盖变化信息选择相似像素，即选择该像素经过土地覆盖变化后的光谱相似像素，而非原始像素的相似像素，与传统方法差异较大。为了避免混淆，本文将相似像素定义为变化同质性像素。

图2以两个时刻Landsat与MODIS图像中某一超像素为例，展示了本文变化同质性像素的选择方法。如图2(a)所示，t1到t2时刻的MODIS图像目标像素的反射率变化为5(反射率从5变化到10)。阈值Tp根据式(4)计算得到，数值为0.5，因此本文选择变化同质性像素的阈值范围为[4.5,5.5]。即在t1时刻Landsat的超像素邻域中，选择与该目标像素反射率之差在阈值范围的像素作为变化同质性像素的备选像素，如图2(a)中标记为红色的反射率为9.6,9.7,10.4和10.5的像素。图2(b)为t1时刻的Landsat图像混合像元分解之后得到的3幅端元丰度图。由于本小节处理的数据为突变数据，假设目标像素已发生较大地物变化，所以应该选择在t1时刻与目标像素分度值差异较大的像素。本文根据实验经验设置端元变化阈值为0.5，选择与目标像素差值大于端元变化阈值的备选像素，最后将3幅丰度图的备选像素取交集作为最终的变化同质性像素，如图2(b)所示。

(a)基于阈值的变化同质像素初步选择 (b)基于混合像元解混的变化同质像素约束图2 基于超像素的变化同质性像素筛选Fig.2 Superpixel-based change pixel filtering

获取变化同质性像素后，需要据其构建加权滤波的权重函数。为充分利用图像的空间与光谱信息，本文综合利用目标像素与变化同质性像素的光谱差异和距离差异构建权重函数。对于坐标为(x,y)的目标像素，其变化同质性像素的坐标可以定义为(xi,yi)，其中i∈[1,2，…，Nx,y]，j∈[1,2，…，Nx,y]。目标像素(x,y)与变化同质性像素(xi,yi)之间的光谱差异Sij表示为：

Sij=|[L(x,y,t1,B)-L(xi,yi,t1,B)]-[M(x,y,t2,B)-M(x,y,t1,B)]|。

(5)

对应的距离权重Dij由变化同质性像素距离目标像素的空间距离定义，公式为：

(6)

式中A为比例系数，调整距离权重与综合权重的限制常数。

变化同质性像素与目标像素光谱差异越小、空间距离越接近，表示与目标像素的相似程度越高，因此应当被赋予更高的权值。基于此，本文将光谱差异和距离差异进行归一化，得到每个变化同质性像素对目标像素的贡献权重Wij，公式为：

(7)

求得归一化权重后，对于突变数据，利用基准时刻t1的低分辨率图像计算预测时刻t2的高分辨率图像，公式为：

(8)

式中LF(x,y,t2-1,B)为基于滤波模型利用t1时刻预测的t2时刻高分辨率图像。同样，可以利用上述方法得到基于t3时刻的基准图像预测得到t2时刻的高分辨率图像LF(x,y,t2-3,B)。

基于1.2和1.3小节的方法，可以分别针对物候变化区域和突变区域得到基于t1时刻图像预测的t2时刻图像，如式(3)与式(8)所示。两者综合得到整体的t2时刻预测图像，公式为：

L(x,y,t2-1,B)=LR(x,y,t2-1,B)+LF(x,y,t2-1,B)。

(9)

同理，可以基于t3时刻的高低分辨率图像，预测得到t2时刻的高分辨图像L(x,y,t2-3,B)。

1.4 时间加权

获得两幅t2时刻的预测图像后，对其进行时间加权融合得到最终的预测结果。本文选用余弦相似度作为权重的赋值标准：将2期低分辨率图像看作向量，计算2个向量的夹角余弦值来评估它们的相似度，即余弦值越大相似性越高，设置的权重越大。余弦相似度不仅可以体现光谱值差异，还可以体现光谱曲线形状上的差异，计算的相似度量更为可靠。具体的时间权重函数Thk表达式为：

(10)

其中，余弦相似度表达式为：

(11)

式中：k为不同时刻的基准图像，值为1或3；cos为夹角余弦值；θ1为t1和t2时刻低分辨率像素向量夹角；θ3为t3和t2时刻低分辨率像素向量夹角。基于以上权重，可综合求得t2时刻的高分辨率图像L(x,y,t2,B)，公式为：

L(x,y,t2,B)=L(x,y,t2-1,B)×T1+L(x,y,t2-3,B)×T3。

(12)

2 实验

2.1 实验数据

本文采用两个Landsat7 ETM+和MODIS数据集对算法性能进行测试，每个数据集包含3幅Landsat图像和3幅MODIS图像。经过投影校正、大气校正、重采样和裁剪等预处理，Landsat图像分辨率为30 m，MODIS图像分辨率为500 m。本文采用绿、红、近红外3个波段进行预测实验，即Landsat的4，3，2波段和MODIS的3，1，2波段。

数据集1实验区位于加拿大北部区域(http://ledaps.nascom.nasa.gov/ledaps/tools/StarFM.htm)，如图3所示。研究区内土地覆盖类型以森林(云杉、松树、白杨)为主，辅之以沼泽和稀疏植被(土壤、岩石、烧伤痕)，植被生长季节短，物候变化迅速。

(a)t1时刻MODIS图像 (b)t2时刻MODIS图像 (c)t3时刻MODIS图像

数据集2位于中国内蒙古区域(https://drive.google.com/open?id=1yzw-4TaY6GcLPIRNF BpchETrFKno30he)。研究区内地物目标丰富，包含农田、牧场、山区、城镇、城市等，如图4所示。可以看出，数据集2相比于数据集1具有更高的光谱异质性和空间异质性，预测难度更高。

(a)t1时刻MODIS图像 (b)t2时刻MODIS图像 (c)t3时刻MODIS图像

2.2 融合评价指标

本文采用的时空融合评价指标包括结构相似性指数(structural similarity,SSIM)、平均绝对差(average absolute deviation,AAD)、方差误差(variance of error,VOE)和相对无量纲全局误差(erreur relative globale adimensionnelle de synthèse,ERGAS)。SSIM表示预测图像与真实图像之间的相似度;AAD代表预测图像与真实图像之间的偏差;VOE表示预测图像的误差波动;ERGAS评估光谱范围内所有波段的光谱质量。SSIM值越大表示算法性能越好，最佳数值为1；而AAD,VOE,ERGAS值越低表示图像融合质量越好。对比算法包括STARFM和ESTARFM两种。

2.3 实验结果和讨论

图5展示了数据集1的预测图像。其中，图5(a)为真实Landsat图像，图5(b)—(d)为3种时空融合算法生成的预测图像。可以看出，3个模型的预测结果整体相似，H-STFM模型获取的预测图像与真实图像整体色调上更接近。而从局部放大图可以看出，虽然STARFM和ESTARFM模型也能够捕捉到图像的整体变化，但是对很多空间细节信息的预测不够准确，且有较严重的颜色失真和色斑，H-STFM模型获取的预测图像更好地保留了地物的颜色和空间细节信息。

(a)真实图像 (b)STARFM (c)ESTARFM (d)H-STFM

图6是以红光波段为例展示了数据集1 的图像反射率预测值与真实值的散点图，其他数据结果展示如表1，其中加粗体为各指标最佳结果。如图6所示，相比于两种对比算法，H-STFM模型得到的散点图中所有的数据都接近1∶1线，即其能够较好地捕捉到物候引起的反射率变化，拟合的效果更好。

(a)STARFM (b)ESTARFM (c)H-STFM

表1 数据集1的定量评估结果Tab.1 Quantitative evaluation results of data set 1

从表1中可以看出，在绿波段和红波段，H-STFM模型得到的预测图像的偏差AAD和误差波动VOE更小，预测效果更好。虽然在近红外波段上H-STFM模型获得的AAD略高于ESTARFM模型，但是H-STFM模型在各个波段上的SSIM和VOE都是最优的，表明H-STFM模型得到的预测图像与真实图像之间整体的相似度最高并且误差波动最小。H-STFM模型得到的ERGAS指标为0.455 5，明显优于两种对比算法，表明H-STFM模型在相关光谱范围内所有波段的光谱质量最高。

图7是数据集2的预测结果图。其中，图7(a)为真实Landsat图像，图7(b)—(d)为3种时空融合算法生成的预测图像。由图7可以看出，ESTARFM和H-STFM模型可以得到较STARFM模型质量更高的预测效果；而在黄色圆圈区域，H-STFM模型的预测效果最好，与真实图像最为接近。在红框标示的局部区域可以看出，H-STFM模型的预测图像更好地保留了真实图像的整体色调，且更好地预测了图像空间信息，而STARFM的预测图像出现了较为明显的斑块效应，ESTARFM模型在色彩饱和度方面的预测效果较差。结果表明，H-STFM模型在异质性区域同样能够获取较好的预测结果，散点图和定量评价结果也从定量的角度验证了这一点。

(a)真实图像 (b)STARFM (c)ESTARFM (d)H-STFM

图8同样以红光波段为例，展示了数据集2的Landsat图像反射率预测值与真实值的散点图，其他数据见表2。与图6类似，H-STFM模型得到的散点图中所有的数据较为聚集且接近1∶1线，数据更为集中，即该模型能够较好地捕捉到突变引起的反射率变化，散点数据拟合的效果最好。

(a)STARFM (b)ESTARFM (c)H-STFM

表2 数据集2的定量评估结果Tab.2 Quantitative evaluation results of data set 2

如表2所示，在红波段，H-STFM模型得到的预测图像的AAD和SSIM更小，表示预测效果最好。H-STFM模型在绿波段和近红外波段上AAD与SSIM表现不稳定，但是在各个波段上的VOE表现稳定且都为最优值，表明通过H-STFM模型获取的预测图像与真实图像之间的相似度最高。与数据集1实验结果类似，H-STFM的ERGAS指标也优于两种对比算法，预测图像的光谱质量最高。

3 结论

针对地表覆盖的正常物候变化和异常突出变化，设计了一种分层预测的策略，提出H-STFM模型。根据在时间段内地表类型的变化情况，将待预测像素分为物候变化像素和突变像素，对两种类型的像素分别进行线性回归和加权滤波预测。同时，对相似像素的选择方法进行改进，本文利用超像素邻域窗口和混合像元分解丰度图进行约束，筛选出更准确的相似像素。最后将时间加权函数进行优化改进，以提高图像融合精度。在两个实验区的实验结果表明，H-STFM模型可以得到更高的预测精度，更有利于预测土地覆盖异常变化。主要研究结论如下：

1)利用分层思想，根据像素变化有针对性构建了不同变化数据的预测方法，预测精度明显提高。

2)引入超像素作为备选邻域筛选变化同质性像素，并利用丰度图约束像素选择，选择相似像素的准确度明显提高。

3)选用余弦相似度作为时间加权函数的度量准则，有效提高了预测精度。

本文算法还存在一定的局限性，如将低分辨率图像建立的回归模型直接应用到了低分辨率图像的求解，存在一定的系统误差。后期研究将进一步考虑不同空间分辨率数据成像系统差异，提高图像预测的精度。