CN111161229A - 一种基于几何主动轮廓模型和稀疏自编码的变化检测方法 - Google Patents

Abstract

一种基于几何主动轮廓模型和稀疏自编码的变化检测方法，主要针对监督变化检测需要人为收集样本以及直接对差异图处理获得变化检测图出现的“椒盐”噪声等缺点。实现步骤：(1)图像读入；(2)差异图生成；(3)感兴趣区域提取；(4)候选样本提取；(5)训练和待分类样本生成；(6)稀疏自编码网络学习；(7)变化检测结果生成。本发明采用几何主动轮廓模型自动提取差异图的感兴趣区域作为候选样本集，并利用模糊C均值聚类生成训练样本和待分类样本。稀疏自编码网络通过学习样本的差异表示对待分类样本进行分类，最后通过标签融合得到变化检测图。整个变化检测过程结合无监督和监督的优势，避免了人为收集训练样本的流程，降低了“椒盐”噪声，提高变化检测精度。

Description

一种基于几何主动轮廓模型和稀疏自编码的变化检测方法

技术领域

本发明属于图像处理技术领域，具体涉及遥感图像变化检测技术领域中的一种基于几何主动轮廓模型和稀疏自编码的变化检测方法。本发明可用于检测两幅不同时相多光谱遥感图像中农田、植被、水体、建筑等变化，可广泛应用在土地利用和覆盖、城市规划、农业监测等。

背景技术

随着遥感技术的进步，卫星图像逐渐向高空间分辨率发展，其数据规模也日益庞大。在众多卫星数据中，多光谱卫星图像以其包含丰富的地表信息和良好的可视化特征广泛应用在对地观测领域。其中，通过多光谱图像分析和识别地表变化有助于人们监控和管理自然资源的利用，这对可持续发展起到了至关重要的作用。然而，在多光谱图像中会出现不同地物具有相同光谱反射率或相同地物具有不同光谱反射率的现象，并且图像分辨率越高，此种现象会越明显，进而增加了多光谱变化检测的难度。同时，利用机器学习和图像处理方法降低变化检测中人力和资源的消耗，实现自动变化检测是值得研究的问题。

现有的变化检测技术主要分为监督变化检测和无监督变化检测。监督变化检测精度较高，但是需要人为收集真实的地物变化标签，有时地物变化标签很难甚至不可获得。无监督变化检测的优点是降低了人工和成本，具有实际的应用价值，然而变化检测精度相比监督变化检测较低。

主动轮廓模型的思想是利用连续的演化曲线逼近图像的目标区域来实现图像分割，其可分为参数主动轮廓模型和几何主动轮廓模型。前者将图像的特征信息包含在能量泛函中，通过曲线内力和外力共同作用下收敛到目标区域，但是分割结果受到初始轮廓的影响。为了解决前者的缺点，后者通过水平集函数演化并取零水平集获得目标分割结果，并在图像分割和变化检测领域中得到了应用。

西安电子科技大学在申请的专利“基于稀疏自动编码机的遥感图像变化检测方法”(专利申请号201510512145.4，公开号：CN105205807A)中公开了一种基于稀疏自动编码机的遥感图像变化检测方法。该方法将读入的合成孔径雷达图像利用对数比值构造差异图，采用差异图对稀疏自编码网络进行无监督训练，最后通过模糊C均值对稀疏自编码网络提取的差异图像素特征聚类得到变化检测结果。该方法不足之处在于利用稀疏自编码网络提取差异图特征过程中会学习到噪声数据，进而会对实验结果产生影响。

CelikTurgay和Kai-Kuang Ma在论文“Multitemporal image change detectionusing undecimated discrete wavelet transform and activecontours.”(IEEETransactions on Geoscience and Remote Sensing,2010,49(2):706-716.)中公开了基于无抽取离散小波变换和主动轮廓模型的遥感图像变化检测方法。该方法将光学图像或合成孔径雷达图像的差异图采用无抽取离散小波变化得到差异图的多分辨率表示，通过主动轮廓模型分割获得最终的变化检测结果。该方法的不足之处在于直接采用主动轮廓模型得到变化检测结果受到差异图生成方式的影响。

针对现有变化检测技术的不足，目前需要一种结合监督和无监督变化检测优势、提高变化检测精度的同时避免人工参与的变化检测方法。

本发明由国家自然自然科学基金项目(61773087)资助。

发明内容

本发明的目的是克服上述现有技术的不足，提出一种基于几何主动轮廓模型和稀疏自编码的变化检测方法。本发明与现有技术相比通过提取差异图感兴趣区域和模糊C均值聚类减少了噪声干扰，避免了人工收集样本的流程，采用稀疏自编码网络提高了变化检测精度。

本发明解决上述问题的技术方案如下：

一种基于主动轮廓模型和稀疏自编码网络的变化检测方法，包括以下步骤：

(1)图像读入：

(1a)读入两幅已配准的不同时期同一地理区域的多光谱遥感图像；

(2)差异图生成：

(2a)利用差异变化向量分析方法生成两时相多光谱图像差异图d，如下式所示：

其中，B是多光谱图像的通道数，n是多光谱图像像素总数，

和

分别代表t₁、t₂不同时期图像第i个通道的第j个像素值。

(2b)将差异图d的像素值归一化到[0,1]之间，如下式生成归一化后的差异图D：

这里，min(·)用来计算差异图d像素中的最小值，max(·)用来计算差异图d像素中的最大值。

(3)感兴趣区域提取：

(3a)根据差异图D的灰度特性，自动生成几何主动轮廓模型的初始轮廓C_area，如下所示：

C_area＝{(x,y)|(x-x₀)²+(y-y₀)²＜r²}

其中，(x,y)是差异图每个像素的坐标位置，r是初始轮廓半径，中心坐标x₀和y₀满足如下条件：

式中，*是卷积算子，h是卷积模板，|·|代表绝对值运算。

(3b)设置迭代次数iter_cv、长度项系数μ和正则化参数λ₁、λ₂，通过最小化几何主动轮廓模型的能量泛函E(C,c₁,c₂)得到感兴趣区域轮廓C：

式中，D(x,y)表示差异图中位置为(x,y)的像素值，L(·)用来约束感兴趣区域轮廓C的长度，c₁和c₂分别表示感兴趣区域Ω_in和背景区域Ω_out的平均灰度。为了得到几何主动轮廓模型的数值解，引入Heaviside函数、Dirac测度和水平集函数φ(x,y)将能量泛函E(C,c₁,c₂)改写为水平集形式E(φ,c₁,c₂)：

其中，▽是梯度算子，H(x)代表Heaviside函数：

δ(x)是Dirac测度：

c₁和c₂分别是感兴趣区域和背景区域平均灰度，计算如下：

(3c)按照(3b)的公式更新水平集演化方程，在满足迭代次数iter，轮廓C包围部分即是差异图感兴趣区域。根据水平集函数φ每个位置(x,y)的数值可以定量确定轮廓C的位置，如下所示：

式中，Ω_in是感兴趣区域，Ω_out是背景区域。

将感兴趣区域与背景区域进行二值化操作得到掩模图像M，如下式所示：

(4)候选样本提取：

(4a)差异图每个像素D(x,y)经过掩模图像对应位置的像素M(x,y)如下式约束得到候选样本集合S：S＝{D(x,y)|M(x,y)＝1}。

(5)训练和待分类样本生成：

(5a)设置最大迭代次数iter_fcm、迭代停止阈值T和聚类类别i_fcm，利用模糊C均值对候选样本集合S＝{s₁,s₂,s₃,...,s_n}聚类，其中，s₁,s₂,s₃,...,s_n是差异图的像素，n是样本总数。通过最小化如下目标函数J(u,c)得到样本的类别标签：

其中，α为模糊权重指数，u_ij是第j个样本属于第i类的隶属度，约束条件为

||·||为取模算子，s_j是集合S的第j个样本，c_i为第i类的聚类中心。最小化目标函数过程中，通过下式更新隶属度u_ij和聚类中心c_i：

式中，c_k是第k类的聚类中心。

(5b)在满足最大迭代次数iter_fcm或迭代停止阈值T后，利用下式得到模糊C均值聚类后变化、未变化和待分类样本标签集合L_s：

根据标签集合

将候选样本集S分为变化S_change、未变化S_unchange和待分类S_unknown样本集合。

(6)稀疏自编码网络学习：

(6a)设置稀疏自编码网络结构、每层迭代次数iter_sae、稀疏系数κ、惩罚权重β和正则化参数λ；

(6b)初始化网络权重，将变化样本集合S_change和未变化样本集合S_unchange堆叠作为稀疏自编码网络的输入，通过下式逐层计算最小化损失函数J(W,b)：

其中，N是训练样本数目，s_i和

分别代表输入样本和重构样本，l₂是隐层网络节点个数，

是第i个节点的平均激活值，l_m和l_m+1分别是第m和第m+1层神经元节点个数，

为第m层节点i和第m+1层节点j的连接权值。

(6c)在满足迭代次数iter_sae，通过下式更新网络的权重

和偏置

其中，

和

分别为更新前和更新后第m层节点i和第m+1层节点j的连接权值，

和

分别是更新前和更新后第m层节点i和第m+1层节点j之间的偏置。η是学习率，

和

利用反向传播算法计算偏导数如下：

式中，

为第m层第j个节点的输出值，

是第m+1层第i个节点的残差。根据(6b)中的公式，m＝1,2，则输出层第i个神经元的残差

和隐藏层第i个节点的残差

如下式更新：

通过更新(6b)和(6c)的公式完成稀疏自编码网络的无监督学习。

(6d)将稀疏自编码网络最后一层的输出作为Softmax分类层的输入，通过样本标签

和

利用反向传播算法微调整个网络的权重和偏置，并对待分类S_unknown样本分类得到分类结果

(7)变化检测标签融合：

(7a)将分类结果

与样本标签

按照下式融合得到变化检测标签L：

(8)变化检测结果生成：

变化检测标签为0和1，将其分别映射到0和255，从而得到变化检测图像。

本发明的益处与优点在于：

第一，利用几何主动轮廓模型自动提取差异图感兴趣目标，减少了差异图噪声干扰，降低了大部分未变化区域的数据量。第二，无监督生成变化和未变化训练样本，避免了监督变化检测需要人为收集样本的缺点。第三，结合稀疏自编码无监督特征提取和监督学习的优势，使用少量样本数据学习提高无监督变化检测的精度，具有实际应用价值。

附图说明

图1是基于几何主动轮廓模型和稀疏自编码的变化检测流程图。

图2是两幅不同时期辽宁省大连市荣兴水库附近landsat8遥感图像。(a)为2014年5月26日获取的荣幸水库附近Landsat8图像；(b)为2017年6月3日获取的荣兴水库附近Landsat8图像。

图3是几何主动轮廓模型初始轮廓和差异图感兴趣目标分割结果。(a)为几何主动轮廓模型初始轮廓图；(b)为感兴趣区域轮廓图。

图4是变化检测仿真实验和标准变化检测图。(a)为CVA-Kmeans方法变化检测图；(b)为PCA-Kmeans方法变化检测图；(c)为本发明变化检测图；(d)为变化检测标准图。

具体实施方式

以下结合技术方案和附图详细叙述本发明的具体实施方式。

参照附图1，一种基于几何主动轮廓模型和稀疏自编码网络的变化检测方法，包括以下步骤：

1)图像读入；

读入两幅已配准的不同时期同一地理区域的多光谱遥感图像。

2)差异图生成；

2.1)利用差异变化向量分析方法生成差异图；

对于通道数目B＝7，像素总数n＝284200的两幅图像，通过差异变化向量分析方法如下式生成差异图：

其中，

和

分别代表t₁、t₂不同时期图像第i个通道的第j个像素值。

2.2)差异图归一化；

将差异图d的像素值归一化到[0,1]之间，如下式生成归一化后的差异图D：

这里，min(·)用来计算差异图d像素中的最小值，max(·)用来计算差异图d像素中的最大值。

3)感兴趣区域提取；

3.1)自动生成初始轮廓；

根据差异图D的灰度特性，设置初始轮廓半径r＝10，自动生成几何主动轮廓模型的初始轮廓C_area，如下所示：

C_area＝{(x,y)|(x-x₀)²+(y-y₀)²＜10²}

其中，(x,y)是差异图每个像素的坐标位置，中心坐标x₀和y₀满足如下条件：

式中，*是卷积算子，h是卷积模板，|·|代表绝对值运算。

3.2)感兴趣区域轮廓演化；

设置迭代次数iter_cv＝50、长度项系数μ＝0.5和正则化参数λ₁＝1，λ₂＝1，通过最小化几何主动轮廓模型的能量泛函E(C,c₁,c₂)可得到感兴趣区域轮廓C：

式中，D(x,y)表示差异图中位置为(x,y)的像素值，L(·)用来约束感兴趣区域轮廓C的长度，c₁和c₂分别表示感兴趣区域Ω_in和背景区域Ω_out的平均灰度。为了得到几何主动轮廓模型的数值解，引入Heaviside函数、Dirac测度和水平集函数φ(x,y)将能量泛函E(C,c₁,c₂)改写为水平集形式E(φ,c₁,c₂)：

其中，▽是梯度算子，H(x)代表Heaviside函数：

δ(x)是Dirac测度：

c₁和c₂分别是感兴趣区域和背景区域平均灰度，计算如下：

3.3)掩模图像生成；

按3.2)的公式更新水平集演化方程，在满足迭代次数iter＝50，轮廓C包围部分即是差异图感兴趣区域。根据水平集函数φ每个位置(x,y)的数值可以定量确定轮廓C的位置，如下所示：

式中，Ω_in是感兴趣区域，Ω_out是背景区域。

将感兴趣区域与背景区域进行二值化操作得到掩模图像M，如下式所示：

4)候选样本提取；

差异图每个像素D(x,y)经过掩模图像对应位置的像素M(x,y)如下式约束得到候选样本集合S：

S＝{D(x,y)|M(x,y)＝1}

5)训练和待分类样本生成；

5.1)候选样本聚类；

设置最大迭代次数iter_fcm＝100、迭代停止阈值T＝1e-5，聚类类别i_fcm＝3和模糊权重指数α＝2，利用模糊C均值对候选样本集合S＝{s₁,s₂,s₃,...,s₂₈₄₂₀₀}聚类，其中，s₁,s₂,s₃,...,s₂₈₄₂₀₀是差异图的像素。通过最小化如下目标函数J(u,c)得到样本的类别标签：

其中，u_ij是第j个样本属于第i类的隶属度，约束条件为

||·||为取模算子，s_j是集合S的第j个样本，c_i为第i类的聚类中心。最小化目标函数过程中，通过下式更新隶属度u_ij和聚类中心c_i：

式中，c_k是第k类的聚类中心。

5.2)获得训练和待分类样本；

在满足最大迭代次数iter_fcm＝100或迭代停止阈值T＝1e-5后，利用下式得到模糊C均值聚类后变化、未变化和待分类样本标签集合L_s：

根据标签集合

将候选样本集S分为变化S_change、未变化S_unchange和待分类S_unknown样本集合。

6)稀疏自编码网络学习；

6.1)网络参数设置；

设置稀疏自编码网络结构1-80-50-2、样本数N＝24580、每层迭代次数iter_sae＝18、稀疏系数κ＝0.1、惩罚权重β＝3和正则化参数λ＝3e-3。

6.2)最小化损失函数；

初始化网络权重，将变化样本集合S_change和未变化样本集合S_unchange堆叠作为稀疏自编码网络的输入，通过下式逐层计算最小化损失函数J(W,b)：

其中，s_i和

分别代表输入样本和重构样本，l₂是隐层网络节点个数，

是第i个节点的平均激活值，l_m和l_m+1分别是第m和第m+1层神经元节点个数，

为第m层节点i和第m+1层节点j的连接权值。

6.3)网络更新权重和偏置；

在满足迭代次数iter_sae＝18，学习率η＝0.01，通过下式更新网络的权重

和偏置

其中，

和

分别为更新前和更新后第m层节点i和第m+1层节点j的连接权值，

和

分别是更新前和更新后第m层节点i和第m+1层节点j之间的偏置。

和

利用反向传播算法计算偏导数如下：

式中，

为第m层第j个节点的输出值，

是第m+1层第i个节点的残差。根据6.2)公式，m＝1,2，则输出层第i个神经元的残差

和隐藏层第i个节点的残差

如下式更新：

通过更新6.2)和6.3)的公式完成稀疏自编码网络的无监督学习。

6.4)稀疏自编码网络分类；

将稀疏自编码网络最后一层的输出作为Softmax分类层的输入，通过样本标签

和

利用反向传播算法微调整个网络的权重和偏置，并对待分类S_unknown样本分类得到分类结果

7)变化检测标签融合；

将分类结果

与样本标签

按照下式融合得到变化检测标签L：

8)变化检测结果生成；

变化检测标签为0和1，将其分别映射到0和255，从而得到变化检测图像。下面结合附图2、附图3和附图4对本发明的效果进一步说明。

仿真实验：

本发明仿真实验的硬件环境为：Intel Core i5处理器，2.5GHz主频，8G内存；软件环境为：Windows7操作系统，Matlab R2016b。

本发明处理的是辽宁省大连市荣兴水库附件Lansat 8图像，如图2中(a)和(b)所示，图像获取时间分别为2014年5月26和2017年6月3日，分辨率为15米，大小均为490×580。

通过以下指标定量评估本发明和对比实验的变化检测结果：(1)FP：真实变化被识别为未变化的像素个数；(2)FN：真实未变化被识别为变化的像素个数；(3)MA：FN占真实变化像素数目C_t的比例；(4)FA：FP占实验中变化像素数目C_e的比例；(5)OA：实验中正确识别出的变化和非变化像素数目T_e占总像素数目n的比例；(6)Kappa系数：实验变化检测图与标准图一致性程度，数值越接近于1一致性越好，如下式计算：

其中，U_t和U_e分别代表真实和实验中未变化像素的个数。

图3是本发明感兴趣区域提取实验，(a)是几何主动轮廓模型初始轮廓，(b)是感兴趣区域轮廓。

图4为变化检测仿真实验，将无监督变化检测CVA-Kmeans和PCA-Kmeans与本发明作对比，(a)是CVA-Kmeans，(b)是PCA-Kmeans，(c)是本发明仿真结果，(d)是标准图。

仿真结果分析：

从图中可以看出，(a)和(b)结果能识别出大部分变化区域，但是均出现了“椒盐”现象，且(b)相对于(a)具有较少的错误像素，其“椒盐”问题主要是无监督聚类和差异图中的噪声所引起的。本发明相比(a)和(b)显著降低了“椒盐”噪声，同时与标准图具有良好的一致性，但是一些区域会有漏检测存在。仿真实验定量评估如表1所示。

表1变化检测精度表

从表中可以看出，本发明Kappa系数和OA较高且具有较低的FA和FP，说明其与标准图的一致性较高。然而，本发明的FN相对较多。CVA-Kmeans和PCA-Kmeans的FP达到了6000以上，具有较高的FA，但是漏检率较低。整体上看，本发明降低了“椒盐”现象，错误变化像素明显减少，在不需要人为干预样本的同时提高了变化检测精度。

以上所述实施例仅表达本发明的实施方式，但并不能因此而理解为对本发明专利的范围的限制。应当指出，对于本领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (1)

1.一种基于几何主动轮廓模型和稀疏自编码的变化检测方法，其特征在于，包括以下步骤：

(1)图像读入：读入两幅已配准的不同时期同一地理区域的多光谱遥感图像；

(2)差异图生成：

(2a)利用差异变化向量分析法生成两时相多光谱图像差异图d，如下所示：

其中，B是多光谱图像的通道数，n是多光谱图像像素总数，

和

分别代表t₁、t₂不同时期图像第i个通道的第j个像素值；

(2b)将差异图d的像素值归一化到[0,1]之间，生成归一化后的差异图D，如下所示：

(3)感兴趣区域提取：

(3a)根据差异图D的灰度特性，自动生成几何主动轮廓模型的初始轮廓C_area，如下所示：

C_area＝{(x,y)|(x-x₀)²+(y-y₀)²＜r²}

其中，(x,y)是差异图每个像素的坐标位置，r是初始轮廓半径，x₀、y₀为中心坐标；

(3b)设置迭代次数iter_cv、长度项系数μ和正则化参数λ₁、λ₂，通过最小化几何主动轮廓模型的能量泛函E(C,c₁,c₂)得到感兴趣区域轮廓C：

式中，D(x,y)表示差异图中位置为(x,y)的像素值，L(·)用来约束感兴趣区域轮廓C的长度，c₁和c₂分别表示感兴趣区域Ω_in和背景区域Ω_out的平均灰度；

为了得到几何主动轮廓模型的数值解，引入Heaviside函数、Dirac测度和水平集函数φ(x,y)将能量泛函E(C,c₁,c₂)改写为水平集形式E(φ,c₁,c₂)：

其中，

是梯度算子，H(x)代表Heaviside函数：

δ(x)是Dirac测度：

c₁和c₂分别是感兴趣区域和背景区域平均灰度；

(3c)按照步骤(3b)更新水平集演化方程，在满足迭代次数iter，轮廓C包围部分即是差异图感兴趣区域；根据水平集函数φ每个位置(x,y)的数值定量确定轮廓C的位置，如下所示：

式中，Ω_in是感兴趣区域，Ω_out是背景区域；

将感兴趣区域与背景区域进行二值化操作得到掩模图像M，如下式所示：

(4)候选样本提取：

(4a)差异图每个像素D(x,y)经过掩模图像对应位置的像素M(x,y)如下式约束得到候选样本集合S：S＝{D(x,y)|M(x,y)＝1}；

(5)训练和待分类样本生成：

(5a)设置最大迭代次数iter_fcm、迭代停止阈值T和聚类类别i_fcm，利用模糊C均值对候选样本集合S＝{s₁,s₂,s₃,...,s_n}聚类，其中，s₁,s₂,s₃,...,s_n是差异图的像素，n是样本总数；通过最小化如下所示的目标函数J(u,c)得到样本的类别标签：

其中，α为模糊权重指数，u_ij是第j个样本属于第i类的隶属度，约束条件为

u_ij∈[0,1]；s_j是集合S的第j个样本，c_i为第i类的聚类中心；

最小化目标函数过程中，通过下式更新隶属度u_ij和聚类中心c_i：

式中，c_k是第k类的聚类中心；

(5b)在满足最大迭代次数iter_fcm或迭代停止阈值T后，利用下式得到模糊C均值聚类后变化、未变化和待分类样本标签集合L_s：

根据标签集合

将候选样本集S分为变化S_change、未变化S_unchange和待分类S_unknown样本集合；

(6)稀疏自编码网络学习：

(6a)设置稀疏自编码网络结构、每层迭代次数iter_sae、稀疏系数κ、惩罚权重β和正则化参数λ；

(6b)初始化网络权重，将变化样本集合S_change和未变化样本集合S_unchange堆叠作为稀疏自编码网络的输入，通过下式逐层计算最小化损失函数J(W,b)：

其中，N是训练样本数目，s_i和

分别代表输入样本和重构样本，l₂是隐层网络节点个数，

是第i个节点的平均激活值，l_m和l_m+1分别是第m和第m+1层神经元节点个数，

为第m层节点i和第m+1层节点j的连接权值；

(6c)在满足迭代次数iter_sae，通过下式更新网络的权重

和偏置

其中，

和

分别为更新前和更新后第m层节点i和第m+1层节点j的连接权值，

和

分别是更新前和更新后第m层节点i和第m+1层节点j之间的偏置；η是学习率；

根据步骤(6b)中的公式，m＝1,2，则输出层第i个神经元的残差

和隐藏层第i个节点的残差

如下式更新：

通过更新步骤(6b)和步骤(6c)的公式完成稀疏自编码网络的无监督学习；

(6d)将稀疏自编码网络最后一层的输出作为Softmax分类层的输入，通过样本标签

和

利用反向传播算法微调整个网络的权重和偏置，并对待分类S_unknown样本分类得到分类结果

(7)变化检测标签融合：

(7a)将分类结果

与样本标签

按照下式融合得到变化检测标签L：

(8)变化检测结果生成：

变化检测标签为0和1，将其分别映射到0和255，最终得到变化检测图像。

Images