基于深度玻尔兹曼机的遥感影像变化检测方法

于 杰，李大成*，和栋材，杨 毅

(1.太原理工大学 矿业工程学院，山西 太原 030024；2.太原理工大学 光电工程学院，山西 太原 030024)

0 引言

遥感影像变化检测是基于遥感技术在不同时间对同一区域成像而判断地表类型是否发生变化的手段,其中双时相变化检测是利用地物变化前后的2期影像进行判别[1]。在国土调查、卫片执法、土地利用、农业生产、城市规划、环境监测和应急管理等领域,遥感影像变化检测具有重要作用[2-3]。随着遥感技术的发展,各种多光谱遥感影像中蕴含大量且丰富的地物信息,使得用于变化检测的特征提取手段越来越多元化[4]。传统基于像元和面向对象的方法很难对各种复杂的地表变化进行充分描述和表达[5],同时基于深度学习的算法在不断地探索和发展[3,5],因此如何提升对遥感影像的特征提取对于准确检测地物变化具有现实意义。

变化检测过程通常包含数据预处理、变化检测方法选择和精度评价3个阶段,而变化检测方法的选择是目前多数研究人员关注的重点,其对检测结果的质量影响较大[6]。基于多光谱遥感影像的变化检测方法研究[7],主要包括代数法、分类法、变换法、深度学习法以及多种方法组合的混合方法等。代数法对于变化差异较大的区域较为敏感,但不易识别地物类型的变化[8];分类法能够识别变化类型和较为深层的变化信息,但结果又依赖于分类精度[9];变换法可以很好地突出变化信息,但降维过程影响变化区域的定位和变化类型的确定[10-11];深度学习方法能够实现自动、多层次、多维度地提取影像更深层次的语义特征,且具有很好的鲁棒性和适应性,能够消除检测结果对变化影像的依赖影响,使变化检测的精度大幅提升[3,5]。

应用深度学习进行变化检测的网络结构具有如下特征:自编码器(Auto Encoder,AE)通常作为特征提取器在无监督的情况下实现变化检测且性能良好[12];深度波尔兹曼机(Deep Boltzmann Machine,DBM)作为一种类似深度信念网络(Deep Belief Network,DBN)但无向的模型,可以通过贪婪学习,从隐藏层单元中提取深层特征[13];卷积神经网络(Convo-lutional Neural Network,CNN)通过构建深度神经网络提取遥感影像的高级特征,实现结合监督学习的方法完成变化检测,许多改进的CNN作为分类器或特征提取器也常用于变化检测,如VGGNet[14]、UNet[15]和ResNet[16]等。虽然大多数深度学习方法能够在无监督或弱监督的情况下获得影像的高阶特征,但变化检测精度还主要依赖于先验知识和标记训练[2]。

基于现有研究,考虑到深度学习在影像特征提取领域的应用优势,提出了一种基于DBM的双时相遥感影像自监督训练变化检测方法。该方法融合了前后时相遥感影像的光谱和纹理变化提取变化区域,将变化标记作为监督样本输入DBM模型进行训练,最终完成变化检测。检测结果对变化地物的边缘提取敏感,检测过程实现了流程化应用。

1 测试数据与预处理

变化检测数据集CD_Data_GZ是一个用于变化检测研究的标准数据集,用于测试和评估变化检测算法的性能[17]。该数据集是2006—2019年期间覆盖中国广州市郊区部分区域的谷歌遥感影像,其中包含红、绿和蓝3个波段,空间分辨率为0.55 m,影像大小为1 006 pixel×1168 pixel～4 936 pixel×5 224 pixel;同时还提供了由人工标注的前后时相影像对应的变化标记,包含建筑物、道路、农田、水体、林地和裸地等多种变化类型,但其更侧重于建筑物的变化标记。其中部分典型示例如图1所示。

图1 CD_Data_GZ变化检测数据集示例Fig.1 Samples of “CD_Data_GZ” change detection datasets

一般在对遥感影像进行变化检测之前,从多光谱原始数据开始需要进行辐射定标、大气校正、正射校正、影像融合、配准纠正、影像裁剪、相对辐射校正和影像向量化等预处理步骤。而使用标准数据集作为测试数据,仅需要完成相对辐射校正和影像向量化2个预处理流程。首先,将每组影像对以T1时刻(前时相)影像为基准,使用直方图匹配方法对T2时刻(后时相)影像完成相对辐射校正,消除前后时相遥感影像因传感器差异、成像角度和大气条件等因素造成的辐射差异;然后使用3×3移动窗口将遥感影像的主成分分析(Principal Component Analysis,PCA)的第一分量(PC1)转换为向量组,解决遥感影像不能直接输入DBM模型参与训练的局限性。影像向量化方法示例[18]如图2所示。

图2 影像向量化方法示例Fig.2 Example of image vectorization method

图2表明,对于同一区域2期遥感影像,使用移动窗口逐像元读取窗口内各像元值并转换为向量,再合并为向量组。若前后时相影像大小为(rows,cols),移动窗口半径为R,中心像元坐标为(i,j),则窗口内各像元坐标取值为(i±R,j±R);影像向量化后得到向量组维度为(M,N),其中M=(rows-2×R)×(cols-2×R),N=2×(2×R+1)2。

本文变化检测方法是根据前后2期遥感影像判别地表是否发生变化,故得到的变化检测结果是二值图。对于变化影像的向量化,得到的向量组维度为(M,2),其中每个行向量第一列数值为0且第二列数值为1时,表示该像元未发生变化;而每个行向量第一列数值为1且第二列数值为0时,表示该像元发生变化。变化影像向量化有且仅有该2类情形。

2 研究方法

基于DBM的遥感影像变化检测方法的整体思路为:相对辐射校正后,将前后时相遥感影像进行PCA,然后使用2期影像的第一主分量分别计算余弦相似角和等价模式局部二值模式(Local Binary Pattern,LBP),将余弦相似角乘以等价LBP的差值,接着利用大津阈值算法分割上述乘积结果,将得到的变化标记作为监督样本输入DBM,最后基于DBM实现对2期遥感影像的自监督训练变化检测,获得最终的变化标记二值图。具体变化检测方案如图3所示。

图3 基于DBM的变化检测方案Fig.3 Change detection scheme based on DBM

2.1 融合光谱和纹理变化的样本生成方法

PCA是一种数据集中和降维的方法,其核心是构造新的坐标系,求解协方差矩阵的特征值和特征向量。在遥感影像中,由于各波段之间相关性较高,通过PCA得到第一分量是将各波段大部分的有用信息集中到第一个主分量中,提取影像光谱信息,从而达到数据增强和降维的作用[10]。

相似性度量一般用于分析2个向量的相似程度,余弦相似角即为一种常见的相似度测量方法[19],是通过计算2个向量夹角的余弦值来衡量其相似程度。得到的余弦值越大,说明2个向量之间的夹角越小,表示向量越相似。在遥感影像中,余弦相似角被应用于图像匹配和变化检测等方面。其具体计算公式为:

(1)

式中:Ai和Bi分别表示向量A和B的分量。

本文使用3×3大小的移动窗口将前后时相影像的第一主成分分量分别向量化,在该过程中,提取移动窗口内对应的2个行向量,计算余弦相似角并赋值给该移动窗口的中心像元,得到余弦相似角灰度图,即利用2期影像的光谱相似性表征光谱变化,其绝对值越大,像元相似性越强,变化性越弱。

LBP是一种纹理特征提取方法,可以有效地对影像的纹理特征进行描述和表达,其基本原理是基于像元的灰度值差异,在移动窗口内将邻域像元按序逐个与中心像元值进行比较,并将比较结果编码为二进制数,最终得到一个二进制模式(特征值),转换为十进制后赋值给LBP对应窗口的中心像元。移动窗口设置为3×3,LBP特征值计算方法示意如图4所示。

图4 LBP特征值计算方法示意Fig.4 Schematic diagram of LBP eigenvalue calculation method

图4中,移动窗口内包含1个中心像元和8个邻域像元,共计可产生28=256种LBP。随着移动窗口的半径扩大,LBP的数量也会急剧增加,使得其直方图变得稀疏分散,这对于纹理特征提取和分析是不利的。故本文采用Ojala等[20]改进的等价模式LBP分别对2期影像进行提取,得到降维的LBP,其具有灰度不变性和旋转不变性的特征,在纹理提取和地物分类等应用领域具有较高的鲁棒性和可靠性。等价LBP认为LBP的分布是不均匀的,而且仅有部分二进制模式有较高的频率。在此基础上定义二进制模式中仅有一次0～1或1～0的跳变时,由于二进制模式对应邻域的闭合性,会出现2次跳变,则将具有相同类型的模式定义为等价模式类中的一种,其区别在于二进制模式1出现的位置不同。当未发生二进制跳变时,出现(00000000)2和(11111111)2这2种等价模式类,除此之外,超过2次跳变的其他所有模式均归为混合模式类,例如(10100101)2发生6次二进制跳变。表1描述了8邻域下的等价LBP模式。

表1 8邻域等价LBP模式Tab.1 Eight neighborhood uniform LBP model

由表1可知,8邻域等价LBP模式共对应58种LBP模式,另加1种混合模式类,故8邻域的等价LBP模式实现了256～59的集中降维。

在2期遥感影像的PC1基础上使用3×3大小的移动窗口分别计算其等价模式LBP,然后使用波段差值计算2期等价LBP灰度图的差异,得到纹理变化灰度图,该变化灰度图描述了前后时相影像的纹理变化。对描述2期影像光谱变化的余弦相似角灰度图和描述纹理变化的等价LBP变化灰度图进行波段乘积运算,获得融合光谱和纹理变化的灰度图;应用大津阈值算法[21]对该变化灰度图进行自动阈值分割,向量化后将其作为DBM模型训练的监督样本。

本文使用变化检测数据集中共包含20组季候性差异的遥感影像。在对每组数据进行监督样本生成时,仅考虑该组数据前后时相影像的光谱和纹理变化,故每个影像组均对应一个独立的监督样本作为训练输入。值得注意的是,本文基于DBM的监督训练是一种自监督策略,即每个监督样本仅应用于各自前后时相遥感影像的监督训练进行变化检测,而各组影像对互相之间独立且无关联。

2.2 基于DBM的变化检测方法

DBM是由多个受限玻尔兹曼机(Restricted Boltzmann Machine,RBM)堆叠构成的深度神经网络,各层之间无向连接。通过逐层贪婪训练的策略,直到整个模型收敛,从2期遥感影像中学习到深层特征;再应用融合光谱和纹理变化得到的变化标记作为监督样本输入模型对网络参数进行微调,进一步提高网络的检测精度[22]。在DBM中,不仅包含自底而上的传播训练,而且包含自顶向下的参数反馈,其具有较高的数据泛化能力和良好的鲁棒性,能够很好地传播模糊输入的不确定性。

RBM包含可见层v和隐藏层h,可见层是输入数据的特征表征,本文中为2期遥感影像第一主成分向量化得到的向量组,即影像的灰度特征;隐藏层是描述可见层与隐藏层对应向量之间的依赖关系,2层节点之间存在权重和偏置项,是进行特征提取的关键学习参数。图5展示了一个由3层RBM构成的DBM网络结构,各层内部神经元之间互相无连接,层与层之间无向连接。

图5 由3层RBM构成的DBM网络结构Fig.5 DBM network structure composed of three-layer RBM

在RBM模型中,给定一组状态(v,h),定义其能量函数为:

(2)

式中:vi、ai分别为可见层第i个向量的状态和偏置,hj、bj分别为隐藏层第j个向量的状态和偏置,wj,i为隐藏层第j个向量与可见层第i个向量的连接权重,θ=(W,a,b)表示RBM中的学习参数。由式(2)可得,可见层向量v的边缘分布如式(3)所示:

(3)

当输入可见层时,各隐藏层向量的激活状态彼此独立,反之给定隐藏层状态时,可见层各向量的激活状态也独立。故根据式(2)和式(3),隐藏层的条件概率(即隐藏层向量的激活概率)可表示为:

(4)

式中:sigmoid函数为激活函数,wk表示权重W的第k组向量。

最大化似然函数是训练RBM的目标,即为了得到最优的学习参数θ,需要将边缘分布函数对每个参数求偏导数。在本文RBM的训练过程中,对最优参数的逼近采用对比散度(Contrastive Divergence,CD)方法,即进行k步Gibbs采样完成权重的更新,使RBM模型收敛。而在Hinton[23]的验证中,一般取k=1即可得到良好的训练效果。为计算得到最优的参数θ,采用梯度上升的方法通过迭代更新各参数进而使模型收敛,各参数迭代更新方法如下式:

(5)

(6)

Δbi+[P(hi=1|v(0))-P(hi=1|v(k))]。

(7)

在给定可见层输入数据v0时,计算隐藏层神经元被激活的概率,以此再反向确定可见层各神经元被激活的概率,从而实现可见层单元的重构,而在DBM中的上层RBM隐藏层(输出层)即为下层RBM的输入层。DBM通过初始化权重设置,使用贪婪学习逐层在RBM中训练优化模型的学习参数,故中间层均受其上下2层影响,且从底层到高层尽可能多地保留特征信息;再将上节生成的监督样本变化标记向量化后作为监督样本输入模型。根据标记对输出层状态进行分类,反向对每层参数进行微调,方法如式(8)所示,使得网络达到更优的状态。在DBM中可见层和隐藏层的状态均为二值状态,故最终的检测结果是变化标记的二值图。

Δwi,j=(vhT-v′h′T),

(8)

经过实验,本文方法对于DBM参数设置如下:采用18-90-36-2的3层网络结构,即通过3×3移动窗口向量化测试数据影像组PC1和变化标记;迭代次数设置为30;因测试数据对应变化标记分布随机且不均匀,选择批处理大小为2 000;初始学习率为0.01,动量学习率为0.5。在此基础上,本文方法实现了变化检测过程的自动化过程,为多光谱遥感影像变化检测流程化应用提供了一种解决方案。

2.3 变化检测精度评价

为验证基于本文方法的变化检测结果精度,将检测结果与数据集标准参考变化影像进行像元统计,并采用总体准确率(Overall Accuracy,OA)、误检率(False Positive Rate,FPR)、漏检率(False Negative Rate,FNR)和召回率(Recall,RE)等4个指标对检测结果进行精度评价,计算公式中各参量见表2所示的变化检测误差矩阵[24]。

表2 变化检测误差矩阵Tab.2 Error matrix of change detection

TP、TN、FP和FN分别表示正确检测变化、正确检测未变化、错误检测变化和错误检测未变化的像元数目。通常,指标OA和RE的值越高表示变化检测精度越高,而FPR和FNR的值越低越好。

3 实验结果与分析

为验证本文方法的可行性与有效性,分别设置变化向量分析方法(Change Vector Analysis,CVA)[8]、迭代加权多元变化检测方法(Iterative Reweighted Multivariate Alteration Detection,IRMAD)[25]、ISODATA(Iterative Self-Organizing Data Analysis Techniques Algo-rithm)分类后变化检测方法[26]和ResNet变化检测方法[16]等4组对照实验。通过对测试数据集所有影像对的变化检测结果与参考变化影像按照变化检测误差矩阵进行像元统计,并计算各方法检测结果对应的各项精度评价指标,结果如表3所示。另外,增加输入DBM的监督样本的精度指标,分析监督样本的可靠性。

实验中传统方法的变化二值图是在变化灰度图的基础上应用大津阈值算法进行自动阈值分割,并使用9×9大小的移动窗口进行平滑和聚合得到。IRMAD方法迭代次数和阈值分别为30、0.001;ISODATA方法类别数目设置为10。而基于深度学习的ResNet方法首先将数据集裁剪成256 pixel×256 pixel大小的样本,重叠度50%;然后按照9∶1的比例划分训练集和验证集,选择ResNet-50作为骨干模型,批处理大小为8,初始学习率为0.01,迭代次数为200。

表3 对照实验变化检测结果精度Tab.3 Accuracy of change detection for control experiments

实验选取8组具有代表性的数据样本,见图6(a)～图6(h),每一行表示一组样本,列分别展示其前后时相影像及相应的变化影像、对照实验各方法的检测结果、本文方法检测结果和输入DBM的监督样本标记;灰度图中白色表示检测变化,黑色表示检测未变化。

监督样本的总体准确率为89.75%,误检率为4.06%,漏检率和召回率分别为72.07%和27.93%。将其输入DBM进行训练后得到最终的变化检测结果,漏检率和召回率精度指标有很大提升,提升幅度均超过了33%。但其误检率由4.06%增加到6.42%,说明输入DBM的监督样本的精度直接影响最终的变化检测结果;总体准确率指标提升0.94%。

由表3可以看出,相较传统方法,基于深度学习的ResNet变化检测方法和本文方法的各项精度指标都表现很好,其中ResNet方法的总体检测准确率取得了最优值,为93.51%;而本文方法在漏检率和召回率指标具有优势,总体的准确率也超过了90%,具有较高的检测能力。虽然基于变换的IRMAD方法的总体检测准确率和误检率指标表现出良好的性能,但其遗漏检测和召回率是最差的。而基于代数的CVA方法和基于分类的ISODATA方法的各项评价指标结果基本一致,表现出较高的错误和遗漏检测。本文方法考虑融合遥感影像的光谱和纹理变化,改善了传统方法仅关注某方面特征表达能力的缺点,并结合DBM通过监督训练的方式进一步提取变化特征,得到的变化检测结果漏检率最低且召回率最高,验证了该方法的可行性和有效性。

图6中数据样本组图6(a)～图6(h)的变化影像主要为新增建筑物变化标记,而部分数据存在未关注其他类别的变化:如图6(a)和图6(h)中存在林草地—裸土的地表变化、图6(b)和图6(c)中存在新增道路和不透水面的地表变化,另外图6(h)中存在新增建筑物遗漏标记的情形。实验结果中部分方法如CVA、ISODATA和监督样本生成方法对于该部分变化同样敏感,能够检测出该类变化,然而用于精度检验的变化影像存在标记缺失问题,使得该部分实际变化被归为错误检测,进而影响其精度评价指标。

从监督样本和本文方法的实验结果示例可以看出,本文方法检测结果依赖于监督样本,同时考虑光谱与纹理变化的监督样本生成方法对于地表变化的边缘敏感,尤其对于建筑轮廓、新增道路和不透水面等具有良好的边缘表达,如图6(d)、图6(e)、图6(f)和图6(g)。CVA方法和ISODATA方法的检测结果具有一定的相似性,对于建筑变化的边缘描述不够具体,密集的变化图斑聚合较多且存在中心孔洞,如图6(d)和图6(e)。而IRMAD方法具有最差的检测效果,存在大量的遗漏检测,这与表3所示的精度指标是一致的,该方法在变化检测中对于地表变化与非变化的表征不具稳定性,主要受前后时相地物光谱差异大小的影响。ResNet方法则是基于训练样本的变化检测,因数据集更多地关注建筑变化,故其检测结果同样在建筑变化检测上具有优势,尤其在图6(h)中将遗漏标记的新增建筑识别出来,但其对于变化边缘的表达同样欠缺,密集变化图斑存在大面积聚合。

通过结合人工判别和定量评价指标对变化检测结果进行分析讨论,本文方法的检测结果在检全和提取变化地物边缘方面具有优势,并且能够提取各类复杂的地表变化,表明本文方法对于地物变化检测具有适用性,同时最终的检测结果依赖于监督样本的准确性。

4 结论

为通过2期多光谱遥感影像判断同一区域的地表类型是否发生变化,考虑充分挖掘遥感影像的深层特征而提出了一种基于DBM的自监督训练变化检测方法,应用标准数据集的测试实验结果表明,本文方法可以应用于双时相遥感影像的变化检测过程,并且取得较高的检测精度。最终可以得出如下结论:

① 本文方法变化检测结果的总体准确率达90%以上,召回率超过60%,表明该方法应用于变化检测过程的可行性和有效性,且整体优于传统的变化检测方法。相比其他基于大量训练的深度学习方法,本文方法在提取各种复杂地表变化时表现更为全面,而非只关注训练样本中的类别变化。同时该方法能够根据前后时相遥感影像自动获得监督样本并输入DBM进行训练,优化网络参数以提高变化特征提取的准确性,为双时相多光谱遥感影像变化检测的流程化应用提供了一种新的方案。

② 综合考虑影像光谱和纹理变化、结合DBM的变化检测方法对于变化地物的边缘敏感,能够获得较为准确清晰的变化地物边界,且变化图斑完整。而本文使用的测试数据集更注重建筑的变化,对于道路、裸地和不透水面变化等关注较少,而该类变化的检测结果在验证中因参考未标记被错分,一定程度上影响了整体精度。由于地物实际变化复杂多样,在深度学习中基于样本监督训练的方法对于变化检测数据集前后时相地表变化的准确标注极其重要,故变化检测标准数据集的变化影像标注和监督样本的确定需要更加准确的先验知识支撑。

③ 本文方法是一种基于DBM的自监督训练变化检测方法。若某一区域存在一组前后时相的2期遥感影像,即可依据其光谱和纹理变化选定该组数据的变化区域作为监督样本输入DBM进行训练,得到最终的变化检测结果。该方法具有广泛的适用性和灵活性,而且摆脱了深度学习中基于大量训练数据进行变化检测方法的数据约束。但监督样本的精度仍然是影响最终变化检测结果的关键,因此如何提升监督样本的检测精度还需要进一步研究。