联合显著性和多方法差异影像融合的遥感影像变化检测

2021-09-24王译著陈朋弟李文国余晓娜

自然资源遥感 2021年3期

王译著，黄亮,2，陈朋弟，李文国，余晓娜

(1.昆明理工大学国土资源工程学院，昆明 650093；2.云南省高校高原山区空间信息测绘技术应用工程研究中心，昆明 650093；3.昆明工业职业技术学院，昆明 650302)

0 引言

多时相遥感影像变化检测是宏观、动态地监测与分析地表变化的主要技术手段之一[1],该技术通过对遥感影像多次观测并识别地表变化特征，得到前后地物的变化信息，进而分析出各类地物的性质、位置、和范围变化[1-2]。目前变化检测在城市扩张、农业调查和灾害分析等领域得到广泛应用[3-4]。

随着卫星遥感影像分辨率不断提高，影像中地物更加多样和复杂，这对变化检测方法提出了更高的要求[4]。目前根据是否分类，可以将变化检测方法分为两种：分类后比较法和直接比较法。分类后比较法需对不同时相遥感影像独立分类，多次分类误差累计会导致变化检测精度变低[5-6]，相比之下直接比较法简单易操作而得到广泛应用。差异影像构建和变化区域提取是直接比较法的两个关键步骤，二者直接决定了后续变化检测的结果精度。其中，差异影像构建是直接比较法的基础和关键。目前，国内外研究人员围绕差异影像构建进行大量研究，取得了一定的成果。如文献[7]和文献[8]通过变化矢量分析(change vector analysis，CVA)构建差异影像，然后分别使用机器学习方法和OTSU法对差异影像进行处理得到最终的变化区域。采用CVA构建差异影像只考虑了光谱特征，对“同物异谱”“异物同谱”较为显著的高空间分辨率遥感影像易造成错检和漏检，难以得到很好的检测结果；文献[9]采用CVA得到各个尺度差异影像图，然后引入模糊融合和两种决策级融合方法进行多尺度融合。但是该方法对融合的尺度选取的要求较高；文献[10]建立一个链式模型来表示光谱斜率差异(spectral gradient difference，SGD)模型，通过该SGD模型知识库对遥感影像进行变化检测。采用SGD构建差异影像不适合分析不同生长时期的植被，同时需要大量的训练样本；文献[11]则将CVA和SGD加权融合构造差异影像，然后采用OTSU方法获得最终变化区域。融合CVA和SGD的差异影像构造方法忽略了空间信息，因此“同物异谱”和“异物同谱”问题仍未得到很好的解决。

综上，CVA和SGD的方法仅考虑光谱和纹理特征，忽略了空间和对比度信息，对“同物异谱”和“异物同谱”的检测精度不高，通过CVA和SGD加权融合获取差异影像又容易将CVA和SGD引起的误差累积，对尺度选取的要求较高。根据上述问题，本文引入基于聚类的联合显著性检测技术(co-saliency)[14]，该技术不仅考虑了对比度信息和空间信息，同时综合考虑了影像间的相关性信息。因此本文提出了一种联合显著性和多方法差异影像融合的多时相遥感影像变化检测方法，该方法先分别采用CVA和SGD构建差异影像；然后将差异影像进行基于聚类的联合显著性检测[13-14]并进行显著性差异影像融合得到联合显著性图；最后采用OSTU法阈值分割[16-17]和闭运算处理得到变化影像，以期获得更优异的变化检测结果。

1 研究方法

本文所提出方法的具体流程图如图1所示。

图1 方法流程图Fig.1 Flow chart of the method

主要包括以下步骤：①利用CVA和SGD算法分别对两期遥感影像进行差异影像构造；②利用基于聚类的联合显著性检测方法的对比度信息和空间信息，对CVA和SGD构建的差异影像进行联合显著性检测，分别得到两幅差异影像的共同显著区域，检测结果可以认为是CVA和SGD共同的变化区域；③利用两幅差异影像的相关性，将CVA和SGD的联合显著性图进行融合得到联合显著性差异影像；④通过OTSU和闭运算步骤处理联合显著性差异影像，得到最终的变化区域和未变化区域；⑤对变化检测结果进行定量评价。

如果两幅遥感影像地物未变化，可认为影像的显著性相同，该技术充分考虑了对比度信息、空间信息和相关性信息。因此，通过将两幅差异影像进行联合显著性检测，检测结果可以认为是CVA和SGD共同的变化区域，结果更具有鲁棒性，在复杂场景中也能将显著性变化目标检测出来[12-15]。

1.1 结合CVA和SGD的联合显著性差异图构建

设两个时相T1和T2的遥感影像分别为G1和G2，则某一像素Xij的灰度值G1=[Xij1(T1)，Xij2(T1)，…，Xijn(T1)]T和G2=[Xij1(T2)，Xij2(T2)，…，Xijn(T2)]T，其中，n是计算过程中相互独立的波段数，X是所对应像元的灰度值。分别通过CVA和SGD构建差异影像DCVA和DSGD，然后通过co-saliency技术获取DCVA和DSGD的联合显著图SCVA和SSGD，并通过将二者融合得到最终的差异图DFusion，具体步骤如下：

1)构建DCVA。计算G1和G2影像的变化向量ΔG，公式为：

(1)

然后计算G1和G2两幅影像的变化强度‖G‖，公式为：

(2)

式中‖G‖为影像的变化强度，数值越大，表明该像素发生变化的可能性大；e为波段。本文采用4个波段的遥感影像数据，即在式(2)中，当n=4时，DCVA=‖G‖。

2)构建DSGD。SGD方法主要是通过计算光谱斜率空间，比较两期遥感影像光谱斜率的强度差异来判断地物是否发生了变化。假设遥感数据有n个波段，波段e与e+1(e≤n-1)之间的斜率g是：

(3)

式中：Re+1和Re分别为波段e+1和e的反射率值；ΔR(e,e+1)为Re+1和Re之间的差值；λe+1和λe为波段e+1和e的波长;Δλ为二者的差。

光谱斜率描述了两个相邻波段之间的趋势。如果g(e,e+1)>0，这表明光谱值从波段e增加到e+1;相比之下，如果g(e,e+1)<0，则表明光谱值从波段e减小到e+1;如果g(e,e+1)=0，反射率保持不变，而更大的Δg(e,e+1)则表示两个相邻波段之间的变化更大。

光谱曲线的形状用光谱梯度向量(spectral gradient vector，SGV)描述，当采用4个波段遥感数据时，SGV将光谱信息从传统的光谱空间投影到梯度空间，公式为：

A=[g(1,2),g(2,3),g(3,4)]T。

(4)

假设A1和A2分别是在时间T1和T2的一个像素的SGV，那么它们之间的差值ΔA可以计算为：

ΔA=A2-A1=[g2(1,2),g2(2,3),g2(3,4)]T-[g1(1,2),g1(2,3),g1(3,4)]T。

(5)

在式(5)中，DSGD等于ΔA的绝对值。作为基于SGD的变化幅度，较大的ΔA表示变化的可能性较高。在计算光谱梯度的变化幅度图像后，建立特定的阈值来检测变化和未变化区域。

3)基于co-saliency技术构建SCVA和SSGD。基于聚类的联合显著性检测方法能构建多特征的相关性，将遥感影像的对比度信息和空间特性检测出来，如果一个地物在多张差异图像上表现出很高的相似性，那么就可以被视为一个共同显著性物体。具体步骤如下：

①分别输入差异影像DCVA和DSGD。

②采用K均值聚类方法分别将DCVA和DSGD分为K个聚类，其中k代表某一聚类。

③计算每个聚类的对比度线索(ωc(k))和空间线索(ωs(k))，公式分别为：

(6)

(7)

式中：‖μk-μi‖2用于计算特征空间；ni为聚类Ci中的像素数量；N为输入图像中的像素数；μk和μi为Ck和Ci类中心；oj为输入图像Ij的中心；高斯核函数Υ用于计算像素zij和图像中心oj的欧几里得距离；方差σ2为输入图像的归一化半径；归一化参数nk为Ck的像素数；pij为输入图像Ij中的像素i；Nj为图像Ij的第j个图像像点；M为输入图像的数量；b(pij)为聚类指数。

④通过式(8)融合对比度线索和空间线索，ωi(k)表示显著性线索，公式为：

p(Ck)=Πωi(k)。

(8)

⑤通过式(9)获得最终的显著性效果图，x为输入图像中的像素，公式为：

(9)

4)融合SCVA和SSGD构建最终的差异影像DFusion。将SCVA和SSGD融合构造差异影像DFusion。使用该方法对CVA和SGD的差异影像联合进行显著性变化检测，不仅包含了对比度信息和空间特性，而且检测结果可以认为是共同显著的。

1.2 基于OTSU的差异影像分析

OTSU算法是按灰度特性将图像分成背景和前景两部分，方差作为灰度分布均匀性的度量，其值越大说明构成图像两部分差别越大，即类间方差越大。任何混淆错分现象的发生均会导致错分概率变小，可将图像中两部分进行最理想分开。

将图像的灰度范围平均分为L级，Ni表示第i灰度级的像素点数，像素点总数及其概率密度分布分别为：

(10)

(11)

将图像分割成两类，目标类为C0={0,1,…,T}，背景类为C1={T+1,T+2,…,L-1}，σ2为类间方差，则最佳阈值T′满足类间方差最大，公式为：

σ2(T′)=max[σ2(T)]，0≤T≤L。

(12)

2 实验及分析

2.1 实验数据与精度评价

实验计算机的配置为英特尔酷睿i7-8700K和英伟达GTX 1050Ti 4G，实验过程由MATLAB2018b实现。本文选择了昆明市呈贡区部分地区的GF-1遥感影像，遥感影像包含蓝、绿、红和近红4个波段，获取时间为2013年12月19日和2015年12月10日，将全色与多光谱数据融合后裁剪成3组数据，空间分辨率为2 m，大小分别为：243像素×396像素，301像素×526像素和296像素×405像素。

为了验证本文方法的有效性和可靠性，设计了5组对比实验，分别为CVA与OTSU结合的方法(CVA-OTSU)、SGD与OTSU结合的方法(SGD-OTSU)、卷积小波神经网络(convolutional-wavelet neural network，CWNN)[21]方法、模糊局部C均值聚类(fuzzy local information C-means，FLICM)[22]方法和主成分分析(principal component analysis，PCA)与K均值(K-means)结合方法(PCA-Kmeans)[23]。

在精度评价方面，本文采用错检率(false alarm，FA)[18]、漏检率(missing alarm，MA)[18]、总体精度(overall accuracy，OA)[18]、Kappa系数[19]、和F-measure系数[20]来定量分析试验结果，同时采用直观比较法定性评价。

FA表示实际为未变化的样本被检测成为发生变化的地物的样本比例，公式为：

(13)

式中：FP为实际正确、变化检测结果错误的样本；TN为实际错误、变化检测结果正确的样本。

MA表示实际为发生变化的样本被检测为未变化的地物的样本比例，公式为：

(14)

式中：FN为实际错误、变化检测结果错误的样本；TP为实际正确、变化检测结果正确的样本。

OA是所有的分类正确的样本数除以总得样本数，公式为：

(15)

Kappa系数用于一致性检验，其计算公式基于混淆矩阵，结果用于衡量分类精度，代表着分类与完全随机的分类产生错误减少的比例，公式为：

(16)

(17)

式中PRE为理论检测结果的一致率。Kappa值越高，代表检测结果精度越高。

F-measure系数为算数平均数除以几何平均数，其数值越大代表变化检测的结果越好，公式为：

(18)

2.2 第一组试验结果

图2为第一组影像变化检测结果。图2(a)和图2(b)分别是2013年和2015年获取的昆明市呈贡区部分地区GF-1号卫星遥感影像；图2(c)是目视解译的变化参考影像；图2(d)—(h)分别是影像通过CVA-OTSU，SGD-OTSU，CWNN，FLICM和PCA-Kmeans方法得到的变化检测结果；图2(i)是本文所提出的变化检测方法。

(a)T1时相影像 (b)T2时相影像 (c)参考影像

根据表1，本文方法具有较高的OA系数、Kappa系数和F-measure系数，其中，OA系数达到94.02%，Kappa系数和F-measure系数均达到60%以上，FA系数仅为4.09%，精度明显高于其他传统的变化检测方法。对比第一组数据结果(图2)的视觉效果可得出，本文方法优于其他传统检测方法，更接近于参考基准影像。与参考影像相比，本文方法对植被有一定的错检，主要原因为T1和T2时相原始影像的植被处在不同生长时期，光谱特征不同，所以将植被错检测为变化区域。

表1 第一组影像精度评价结果Tab.1 Accuracy evaluation of the first group image (%)

2.3 第二组试验结果

图3为第三组影像变化检测结果，图3(a)和图3(b)分别是2013年和2015年获取的昆明市呈贡区部分地区GF-1号卫星遥感影像；图3(c)是目视解译的变化检测参考影像；图3(d)—(h)分别是影像通过CVA-OTSU，SGD-OTSU，CWNN，FLICM和PCA-Kmeans方法得到的变化检测结果；图3(i)是本文所提出的变化检测方法。根据表2，本文方法具有较高的OA系数、Kappa系数和F-measure系数，其中，OA系数达到99.01%，Kappa系数和F-measure系数均达到70%以上，FA系数仅为0.59%。明显优于其他传统的变化检测方法。对比第二组数据结果(图3)的视觉效果可看出，本文方法更接近于参考影像，但是与参考影像相比本文方法未能将中间斜体建筑物检测出来，所以漏检率达到22%。主要原因是原始影像中的建筑物光谱不一致，T1原始影像为灰色，T2原始影像为深蓝色，均相融于背景颜色，导致该建筑物漏检。

(a)T1时相影像 (b)T2时相影像 (c)参考影像

表2 第二组影像精度评价结果Tab.2 Accuracy evaluation of the second group image (%)

2.4 第三组试验结果

图4为第三组影像变化检测结果，图4(a)和图4(b)分别是2013年和2015年获取的昆明市呈贡区部分地区GF-1号卫星遥感影像；图4(c)是目视解译的变化检测参考影像；图4(d)—图4(h)分别是影像通过CVA-OTSU，SGD-OTSU，CWNN，FLICM和PCA-Kmeans方法得到的变化检测结果；图4(i)是本文所提出的变化检测方法结果。

(a)T1时相影像 (b)T2时相影像 (c)参考影像

根据表3，本文方法具有较高的OA系数、Kappa系数和F-measure系数，其中，OA系数达到97.80%，Kappa系数和F-measure系数在60%左右，FA系数仅为1.35%。对比第三组数据结果(图4)的视觉效果可看出，本文方法更接近于参考影像，错检区域远少于其他传统方法。但是与参考影像相比，本文方法将一些噪点错检为变化区域，噪点可以后续形态学去除。

表3 第三组影像精度评价结果Tab.3 Accuracy evaluation of the third group image (%)

2.5 实验分析

传统方法中，CVA-OTSU和SGD-OTSU方法仅仅考虑纹理和光谱特征，没有考虑空间和对比度信息，所以有大量错检和漏检；根据文献[6]和[21]，CWNN对合成孔径雷达影像有较好的检测结果，但对光学图像的效果还有待提升；FLICM和PCA-Kmeans对于光谱信息更为复杂的高分辨率遥感影像还需要进一步研究提升检测精度。而本文方法不仅考虑光谱和纹理特征，而且结合对比度线索、空间线索和相关性线索较好地提升了变化检测精度，尤其是本文方法的错检率不足5%，证明了本文方法的检测精度高和可行性，是一种效果优异的变化检测方法。