时空融合技术在城市精细变化检测中的应用研究

古春霞，李大成，2*，赫华颖，杨 毅

(1.太原理工大学 矿业工程学院，山西 太原 030024；2.高分辨率对地观测系统山西数据与应用中心，山西 晋中030600；3.中国资源卫星应用中心，北京 100094；4.太原理工大学 物理与光电工程学院,山西 太原 030024)

0 引言

遥感图像变化检测作为一项关键技术，在自然资源管理、生态环境监测、基础地理数据库更新、目标动态监测和可持续发展规划等诸多领域中有着非常广泛的应用[1-3]。高空间分辨率遥感影像变化检测可提供丰富的空间细节信息，得到精细的变化检测结果。但由于天气、卫星重访周期及传感器自身物理机能对遥感影像质量的影响，难以获得高时序的高空间分辨率影像，无法为时序性变化检测提供支持。针对高分辨率时序影像缺失问题，可以使用时空融合算法重建缺失遥感数据，为变化检测提供影像支撑。现今时空融合技术已在各个领域得到广泛的应用，包括特征识别、目标监测、城市扩张、森林砍伐监测、作物制图和产量估算以及洪水和湿地制图或气候变化等[4]。

遥感图像时空融合的基本思想是利用特定的算法合成2幅或多幅遥感图像并对其进行相关运算，从而重构出新的预测时期的图像。在融合过程中，通常使用一对或多对高低分辨率影像融合生成预测时刻的高分辨率影像。其中，高空间分辨率影像可以提取空间上的细节信息，低空间分辨率影像可以较为精确地描述地物的时序变化情况，以此获取随时间变化的地物空间特征。遥感影像时空融合可以得到同时具有高空间分辨率和高时间分辨率特征的遥感影像序列，实现对土地利用类型变化的精确记录。基于目前已有时空融合模型，可以大致将其分为3类：基于变换的模型、基于重建的模型和基于学习的模型[5]。

在基于重建的时空融合模型中应用最广的是Gao等[6]2006年提出的时空自适应融合模型(Spatial and Temporal Adaptive Reflectance Fusion Model,STARFM)，该方法通过邻域加权过程合成每天的高分辨率影像，融合时间序列影像的质量取决于高时间分辨率影像集的数目和无云数据对的可用性。Zhu等[7]在STARFM的基础上提出了改进型时空自适应融合模型(Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model,ESTARFM)，进一步提高了预测精度。Emelyanova等[8]发现在空间变异为主的地区ESTARFM比STARFM融合效果更好，但在以时间变异为主的地区ESTARFM融合效果一般。Zhu等[9]基于解混和样条插值提出了灵活的空间数据融合算法 (Flexible Spatiotemporal Data Fusion Model，FSDAF)，该算法能预测异质地表区域的地物平缓变化与类型变化。

在基于变换的时空融合模型中主要是小波变换，即对影像进行小波分解，然后对分解后的各层进行融合处理，以此实现影像的融合[10]。在基于学习的时空融合模型方面，Huang[11]等最早将稀疏表达与字典学习用于时空融合，以此提出了基于稀疏表示和稀疏编码的时空反射率融合模型(SParse-representation-based Spatio Temporal Reflectance Fusion Model，SPSTFM)，通过利用从2种影像类型生成的影像块创建的字典来学习高空间分辨率图像和高时间分辨率影像之间的差异。

随着时空融合技术的不断成熟，不少学者将时空融合技术应用在变化检测方面。其中Wang等[12]提出了一种基于多尺度融合的面向方向的变化检测方法，该实验是在高分辨率机载和SPOT 5遥感影像上进行的，与不同的面向对象和面向像素的检测方法相比，实验结果验证了该方法的有效性和可靠性，但该研究仅是针对较小区域进行研究分析。Lal等[13]通过组合稀疏字典融合和K均值聚类算法对多时相遥感影像进行半监督变化检测，结果得出该方法可以减少一定错误分类，但执行速度较慢。Hilker等[14]使用STAARCH对Landsat TM/ETM和MODIS反射率数据进行时空融合并进行变化检测，该方法主要针对于植被变化检测，结果表明结合STAARCH进行的变化检测能较精细地识别景观中的时空变化。Guo等[15]使用STARFM和ESTARFM构建Landsat NDVI时间序列，提高了Landsat NDVI时间序列的质量，并弥补近实时环境监测任务与数据稀缺性之间的差距，同时提出需考虑不同融合算法的有效性在生长季节和非生长季节存在一定差异。综上所述，运用时空融合技术进行变化检测已成为一种可行且有效的手段。通过时空融合技术得到遥感数据集，很大程度上使变化检测更为准确。

现有的时空融合技术应用于变化检测领域主要是对某一特定变化区域进行试验或对区域植被变化进行检测，且融合生成数据基本为中分辨率数据，缺乏高分辨率数据对精细变化检测的应用。针对城市精细变化检测，采用时空融合技术重建高时空分辨率数据，并利用迭代加权多元变化检测算法(Iteration Reweight Multivariate Alteration Detection,IR-MAD)[16]对原始影像和重建影像进行变化检测，分析二者变化检测精度差异，探究时空融合模型在城市精细变化检测中的可用性，验证其在变化检测中的应用精度。

1 研究区和数据源

本文选取山西省太原市中部地区为研究区，地理覆盖范围为37°44′N～38°00′N，112°19′E～112°41′E，如图1所示。该研究区属于温带半干旱季风性气候，研究区域内主要土地覆盖类型包括人工建筑、草地、林地、水体和耕地。本文研究区以城区为主，研究区西部和东北部有部分林地，研究区东南部地物类型多样，农田、林地与居民地混合分布。

图1 研究区地理位置Fig.1 Geographical location of study area

本文选取研究区5个时刻的GF-1 WFV和GF-2多光谱无云遥感影像为研究数据，具体如表1所示。其中，因在本文研究区中2018年没有时相完全重合的两景无云影像，故2018年所选取的GF-1WFV与GF-2多光谱影像的时间相差2天。

表1 遥感数据信息Tab.1 Remote sensing data information

2 实验方法

2.1 基于单数据对的时空融合模型

2.1.1 STARFM模型

STARFM是Gao等在2006年提出的时空融合模型，该模型是输入基准时相的高低分辨率数据和预测时相的低分辨率数据，确定移动窗口大小，搜索相似像元，继而计算邻近相元的距离权重和转换系数，最后得出预测时刻的高分辨率影像L(xw/2,yw/2,t0):

L(xi,yi,tk)-M(xi,yi,tk))，

(1)

式中，(xw/2,yw/2)是移动搜索框的中心像元；M(xi,yi,t0)是低空间分辨率数据在t0时刻的反射率值；L(xi,yi,tk)和M(xi,yi,tk)分别是高空间分辨率和低空间分辨率数据在tk时刻的反射率值;Wijk是权重函数：

(2)

式中，Wijk是权重系数，需要由光谱差异Sijk，时相差异Tijk，中心像元(xw/2,yw/2)和候选像元(xi,yi)在tk时的空间距离dijk来计算:

Sijk=|L(xi,yi,tk)-M(xi,yi,tk)|，

(3)

Tijk=|M(xi,yi,tk)-M(xi,yi,t0)|，

(4)

(5)

2.1.2 FSDAF模型

(6)

(7)

R(xi,yi)=ΔC(xi,yi)-

(8)

(9)

CW(xij,yij)=Eho(xij,yij)×HI(xij,yij)+

R(xi,yi)×[1-HI(xij,yij)]，

(10)

(11)

(12)

2.2 迭代加权多元变化检测算法

Morton Canty和Allan Nielsen在MAD算法[17]基础上，结合EM算法[18]，提出了IR-MAD算法。主要计算方法如下：

将多光谱图像在不同波段中的观测值表示为随机向量F和G，生成随机变量U，V计算MAD变量为：

(13)

式中，ai，bi是对应的特征向量。

再由随机变量Z代表标准化MAD变量的平方和：

(14)

而随机变量Z是具有K自由度的卡方分布，因此得出：

Pr(change)=Px2;K(z)，

(15)

式中，x2代表卡方分布，z值越小，变化概率越高。不变概率可以用作观测值的权重，并且整个过程将反复进行，直到满足停止标准为止。

3 实验结果与分析

先是通过STARFM和FSDAF两种时空融合模型融合得到其余4个时相的遥感影像，对融合后影像与原始影像进行定性与定量分析，分析融合影像质量；然后使用IR-MAD算法对重建后的影像进行变化检测，并与真实变化检测进行对比，分析其实用性及准确性。

3.1 时空融合结果与分析

在融合重建过程中，以2017年5月5日的高低分辨率数据作为STARFM模型与FSDAF模型的基础时相数据，对其他4个时相的影像数据进行了重建，重建影像结果如图2和图3所示。在时空融合结果质量定性评价中，通过目视对比，可以看出融合后的影像色彩失真较小，纹理特征明显。

(a)2015-04-23

(a)2015-04-23

时空融合结果质量定性评价还可以通过影像波段散点图(如图4和图5)分析，从图4和图5可以看出，融合后的影像与真实影像的离散度较好，二者均说明重建影像的可靠性。由于重建影像波段较多，文中只对其中一个波段的散点图进行展示。

在时空融合结果质量定量评价中，通过绝对平均偏差AAD、均方根误差RMSE、结构相似性SSIM[19]、峰值信噪比PSNR以及相对整体位数综合误差ERGAS[20]分析重建影像质量。其中AAD、RMSE、ERGAS指标值越小，说明融合结果的光谱误差越小；PSNR值越大，则表明融合结果失真程度越小；SSIM值越大，表明融合结果与真实结果匹配程度越高[21]。融合影像评价指标如表2所示。

(a)2015-04-23

(a)2015-04-23

表2 融合影像评价指标Tab.2 Fusion image evaluation index

分析表2可知，2种时空融合模型重建结果AAD均在0.027左右，RMSE均在0.029左右，ERGAS均在37左右，可以得出融合后影像误差值较低；PSNR均在25以上，SSIM均在0.82以上，可以得出融合后影像失真程度较小，与真实影像匹配程度较高；5种指标值均表明融合后影像质量较高，可用于进一步研究；分析不同融合模型的融合精度，从各个指标中可以看出，FSDAF模型融合精度略高于STARFM。在本文中使用的2种时空融合算法，融合得到的影像与真实影像匹配度较高，经过精度评定，证明了2种时空融合模型的实用性。

3.2 变化检测结果与分析

本文使用IR-MAD算法分别对2015年和2018年、2016年和2019年2组时间段融合后的影像和原始影像进行变化检测分析，对比融合后的变化检测图和原始影像变化检测图之间的差异，通过对变化检测精度进行量化分析，以此来评价变化检测算法性能。在研究区域内均匀随机选取1 000个验证点验证变化检测的准确性，构建变化混淆矩阵计算5种量化指标分析变化检测精度，分别为变化检测成功率、错检率、漏检率、变化检测总精度和Kappa系数。其中，变化检测成功率是检测正确变化像元数与检测到的变化像元个数的比值；错检率为实际未变化却检测为变化的像元数与检测到的变化像元个数的比值；漏检率为实际变化却检测为未变化的像元数与实际变化像元数的比值；变化检测总精度为检测正确像元数与总像元数的比值；Kappa系数为离散程度。不同时相变化检测图如图6和图7所示，图中红色为变化区域，黑色为未变化区域。对应的变化检测精度评定结果如表3和表4所示。

(a)原始影像变化检测结果

(a)原始影像变化检测结果

对比原始影像变化检测结果与重建影像变化检测结果，二者相差较小，证明了时空融合方法在区域变化检测中的可用性。从表3和表4可以看出，整体变化检测总精度均在87%以上，Kappa系数均在0.75以上。对比时空融合重建影像的变化检测结果精度与原始影像变化检测结果精度，二者整体差值较小，总精度最大相差2%，Kappa系数最大相差0.05，证明了使用时空融合重建遥感影像进行变化检测的可靠性。不论研究区变化大小，从各个变化检测结果指标均可看出，检测总精度相差较小，说明了使用时空融合重建遥感影像对不同变化程度的区域进行变化检测的可靠性。对比不同融合模型所得影像的变化检测结果指标，可以看出FSDAF重建影像变化检测结果指标略大于STARFM，与融合影像精度一致，证明变化检测精度受时空融合模型的影响，其随时空融合精度的增大而增大。

表3 2015-2018年不同影像变化检测结果指标Tab.3 2015-2018 different image change detection result indicators

表4 2016-2019年不同影像变化检测结果指标Tab.4 2016-2019 different image change detection result indicators

4 结束语

对城市实施精细变化检测有利于城市规划和土地资源管理。本文提出使用时空融合技术重建得到的高分辨率影像进行城市精细变化检测，该方法弥补了高质量高时空分辨率影像的缺失问题。研究表明，时空融合技术为城市精细变化检测提供了有力的数据支撑，对不同变化程度区域的变化检测都有较高的精度，分析重建影像的变化检测结果精度，得到整体变化检测总精度均在87%以上，Kappa系数均在0.75以上，且变化检测精度受时空融合模型精度的影响，其随时空融合精度的增大而增大，证实了采用时空融合重建影像进行变化检测的可用性。