一种无监督SAR图像变化检测方法

摘要

本发明提出了一种无监督SAR图像变化检测方法，用于解决现有技术中存在的检测准确性和计算效率较低的技术问题。实现步骤为：对配准后的两时相SAR图像滤波；产生初始差异图；对初始差异图进行显著性检测获得显著性差异图；利用模糊C均值聚类算法对显著性差异图进行预分类获得候选的训练样本以及不确定样本；对候选的候选训练样本集做均衡化处理获得训练样本；用PCA滤波器提取训练像样本征和不确定样本特征；用训练样本特征训练支持向量机，用训练完成的支持向量机对不确定样本分类，获得变化检测结果。本发明能够提高算检测准确性同时减少运算时间，可用于评估灾害和预测灾害发展趋势、目标探测以及土地覆盖与利用监测等。

说明书

技术领域

本发明属于图像处理技术领域，涉及一种图像变化检测方法，具体涉及一种改进的PCANet无监督SAR图像变化检测方法，可用于评估灾害和预测灾害发展趋势、目标探测以及土地覆盖与利用监测等。

背景技术

图像变化检测指利用多时相获取的同一地表区域的遥感影像及其它辅助数据来确定和分析地表变化，是更新地理数据的关键技术，在评估灾害和预测灾害发展趋势以及土地覆盖和利用监测方面发挥着重要作用。以合成孔径雷达(SAR)为代表的微波遥感是有源微波成像传感器。与光学遥感等微波遥感相比，SAR成像技术不仅具有可变侧视角，并且可以在不受天气影响的情况下获得详细的地理信息，因此在变化检测技术上有着广泛的应用。根据人工标记样本的需要，可将图像变化检测方法分为有监督式、半监督式和无监督式三种。有监督和半监督式方法对标记样本具有很强的依赖性，在SAR图像处理中，标记样本不易获取。无监督方法不需要标记样本和人工干预，因此在SAR图像变化检测方法中更常用。

无监督方法在变化检测的应用中，需要先将所有样本自动分为不同的类别，再选出一些具有代表性的样本作为训练样本，用于分类器的训练。这一过程中，训练样本的选取会影响变化检测准确度。

由于深度学习方法在图像分类和人脸识别领域具有良好表现，SAR图像变化检测中的分类器可以利用深度学习方法。Gao F.等人于2016年在《Geoscience Remote SensingLetters》上发表一篇题为“Automatic change detection in synthetic aperture radarimages based on PCANet”的文章首次利用PCANet对SAR图像进行变化检测。该方法首先用Gabor滤波器对原始图像提取Gabor特征，用模糊C均值聚类方法根据图像Gabor特征将像素预先分为变化类、非变化类和不确定类，然后使用PCANet模型对属于不确定类的像素进行分类。该方法能够有效对不确定像素分类，获得准确的变化检测结果。然而仍然存在以下不足：(1)因为SAR图像受相干斑噪声影响，该方法从原图像选取的训练样本不可靠，且特征提取时该方法没有考虑训练样本不均衡的问题，导致对各类样本特征提取的不完整，另外该方法采用基于像素的方式，没有考虑像素邻域信息，会造成部分像素错分，从而使变化检测准确性降低；(2)该方法在所有像素中选择训练样本，导致处理时间较长。

发明内容

本发明的目的在于克服上述已有技术的不足，提出了一种无监督SAR图像变化检测方法，用于解决现有技术中存在的检测准确性和效率较低的技术问题。

本发明的技术思路是：用三维块匹配滤波对配准后的两时相SAR图像进行滤波处理；产生初始差异图；对初始差异图进行显著性检测获得显著性差异图；利用模糊C均值聚类算法对显著性差异图进行预分类获得候选的训练样本以及不确定样本；利用一致选取法对候选的候选训练样本集做均衡化处理获得训练样本；用PCA滤波器提取训练像素特征和不确定像素特征；用训练像素特征训练支持向量机，用训练完成的支持向量机对不确定像素分类，获得变化检测结果。

根据上述技术思路，实现本发明目的采取的技术方案包括如下步骤：

(1)对两幅待检测的SAR图像进行滤波：

对拍摄于同一地点不同时刻的两幅SAR图像进行配准，并对配准后的SAR图像进行滤波，得到滤波后的SAR图像im₁和im₂；

(2)计算滤波后的SAR图像im₁和im₂的初始差异图D₁；

(3)对D₁进行显著性检测，并对检测结果进行二值化：

对D₁进行显著性检测，得到D₁的显著图D_s，并对D_s进行二值化，得到二值化显著图D_s'；

(4)获取D₁的显著性差异图D₂：

利用二值化显著图D_s'，提取初始差异图D₁中的显著性部分，得到D₁的显著性差异图D₂：

D₂＝DotM(D_s'，D₁)

其中，DotM(D_s'，D₁)表示D_s'和D₁同一位置像素的灰度值相乘；

(5)对D₁的显著性差异图D₂中的像素进行预分类：

对显著性差异图D₂中的像素进行聚类，得到正样本集、负样本集和不确定样本集，并将正样本集和负样本集合并为候选训练样本集；

(6)对候选训练样本集进行均衡化：

对候选训练样本集中的正样本集和负样本集进行均衡化处理，得到均衡化的训练样本集，并从其中随机选取30％的样本作为训练样本集；

(7)提取训练样本特征和不确定样本特征：

(7a)对im₁进行逐像素有重叠划分，得到由大小为s×s的多个图像块组成的集合I¹，同时对im₂进行逐像素有重叠划分，得到大小为s×s的多个图像块组成的集合I²：

其中，s为大于等于1的奇数，pos表示图像块序号，pos＝1,2,...,T，T表示SAR图像im₁或im₂中图像块的个数；

(7b)将训练样本集和不确定样本集合并为样本集，并从I¹中选择与样本集中样本位置x_a相对应的图像块同时从I²中选择与样本集中样本位置x_a相对应的图像块a表示样本序号，且a＝1,2,...,N_E，N_E为样本集中样本的总数；

(7c)对进行逐像素有重叠划分，得到由大小为m×m的多个邻域块P_i¹，同时对进行逐像素有重叠划分，得到大小为m×m的多个邻域块P_i²，并将邻域块P_i¹和P_i²合并成的邻域块P_i组合到集合P_temp中：

其中，m为大于等于1小于等于s的奇数，i表示邻域块序号，且i＝1,2,...,n，n表示邻域块个数；

(7d)对P_temp中的邻域块P_i进行去均值后再进行向量化，得到包括由组成的向量矩阵P：

(7e)计算PP^T的特征向量，并将PP^T前N₁个特征向量作为第一层PCA滤波器W_l¹：

W_l¹＝mat(q_l(PP^T))∈R^2m×m>1

其中，N₁为正整数，q_l(PP^T)表示PP^T的第l个特征向量，mat(v)表示把向量映射成为一个矩阵W∈R^2m×m；

(7f)通过第一层PCA滤波器W_l¹计算样本P_i的特征信息P_i^l：

P_i^l＝P_i*W_l¹

其中，*表示三维卷积操作；

(7g)对第一层滤波器W_l¹输出的特征信息P_i^l进行去均值后再向量化，得到包括由组成的向量矩阵Q：

(7h)计算QQ^T的特征向量，并将QQ^T前N₂个特征向量作为第二层PCA滤波器

其中，q_k(QQ^T)表示QQ^T的第k个特征向量，k＝1,2,...,N₂，N₂为正整数，mat(v)表示把向量映射成为一个矩阵W∈R^2m×m；

(7i)通过第二层PCA滤波器计算特征信息P_i^l的二次特征信息

(7j)采用赫维赛德阶跃函数对特征信息矩阵进行二值化，并将得到的二值化信息矩阵转换为每个位置的数值为的矩阵O后，计算P_i的特征f_i，并将作为训练样本特征，作为不确定样本特征，其中：

其中，表示对的直方图统计，N_T表示训练样本集中样本的个数；

(8)获取SAR图像的变化检测结果：

将训练样本特征作为支持向量机的输入，得到训练好的支持向量机，并将不确定样本特征作为训练好的支持向量机的输入，得到SAR图像的变化检测结果图。

本发明与现有技术相比，具有以下优点：

(1)本发明在获取训练样本过程中，采用以邻域块的灰度值为聚类特征，避免了现有技术以原图像的Gabor特征为聚类特征造成训练样本不可靠的缺陷；同时在特征提取过程中，一方面采用一致选取法对候选训练样本集进行均衡化，能够有效消除样本不均衡问题，避免了现有技术直接对样本提取特征造成的提取的各类样本特征提取不完整的缺陷，另一方面采用基于邻域块的方式，避免了现有技术采用基于像素的方式造成的部分像素错分的缺陷，实验结果表明本发明能够有效提高变化检测准确性。

(2)本发明用Context-aware显著性检测方法提取显著性区域，在显著性区域中选择训练样本和不确定样本，避免了现有技术从整幅图像中选取训练样本和不确定样本造成样本数量大的缺陷，实验结果表面本发明能够有效降低运算时间，提高图像变化检测的效率。

附图说明

图1为本发明的实现流程图；

图2为本发明和现有技术应用于两时相Bern数据Real SAR图像的变化检测结果对比图；

图3为本发明和现有技术应用于两时相黄河入海口数据Real SAR图像的变化检测结果对比图；

图4为本发明和现有技术应用于两时相Ottawa数据Real SAR图像的变化检测结果对比图。

具体实施方式

下面结合附图和具体实施例，对本发明作进一步详细描述：

参照图1.一种无监督SAR图像变化检测方法，包括如下步骤：

步骤1)对两幅待检测SAR图像进行滤波：

图像im^o₁为时刻1拍摄的SAR图像，图像im^o₂为时刻2拍摄的SAR图像，im^o₁和im^o₂拍摄与同一地点，对im^o₁和im^o₂进行配准，采用三位块匹配滤波方法，对配准后的SAR图像用进行滤波，去除图像中的相干斑噪声并保持图像边界信息，得到滤波后的图像im₁和im₂：

im₁＝{im₁(i,j)|1≤i≤H,1≤j≤W}，

im₂＝{im₂(i,j)|1≤i≤H,1≤j≤W}，

H表示已获得的SAR图像的高度，W表示已获得的SAR图像的宽度，(i,j)表示图像中像素的位置。

步骤2)获取滤波后的SAR图像im₁和im₂的初始差异图D₁：

采用对数比算法，计算滤波后的SAR图像im₁和im₂的初始差异图D₁，计算公式为：

步骤3)对D₁进行显著性检测，并对检测结果进行二值化：

利用Context-aware显著性检测方法对D₁进行显著性检测，得到D₁的显著图D_s，并对D_s进行二值化，得到二值化显著图D_s'，计算公式为：

其中，p(x,y)表示显著图D_s所在位置为(x,y)的像素的灰度值，τ表示通过大津算法计算获取的阈值，1表示显著性像素，0表示非显著性像素。

步骤4)提取D₁中的显著性的部分：

利用二值化显著图D_s'，对初始差异图D₁中的显著性部分进行提取，得到显著性差异图D₂：

D₂＝DotM(D_s'，D₁)

其中，DotM(D_s'，D₁)表示D_s'和D₁同一位置像素的灰度值相乘。

步骤5)对显著性差异图D₂中的像素进行预分类：

步骤5a)以显著性差异图D₂中的像素x_i为中心，将D₂划分成大小为m×m的邻域块m的取值为正整数，本发明实施例中m设置为3，可以有效抑制噪声并保持图像边界信息；

步骤5b)采用模糊C均值聚类算法，以中所有像素的灰度值为聚类特征，对D₂中的像素进行聚类，得到n类像素，n为大于等于3的整数，为了保证样本可靠性，本发明实施例中n设置为5，则得到5类像素，将这五类根据聚类中心模值的大小按降序排列得到c₁,c₂,...,c₅；

步骤5c)将c₁中的每个像素标记为正样本，得到正样本集c_P，同时将c₅中的每个像素标记为负样本，得到负样本集c_P，将c₂,c₃和c₄中的像素标记为不确定样本得到不确定样本集c_U。

步骤5d)将c_P和c_N合并为候选训练样本集。

步骤6)对候选训练样本集进行均衡化：

步骤6a)假设正样本集c_P中样本个数为n_P，负样本集c_N中样本个数为n_N，且n_P＜n_N，复制比例为N：

其中，[·]表示取整函数；

步骤6b)将c_P中的样本复制N份，组成集合c_P'，对c_P'和c_N中样本进行合并，并从合并结果中随机选取30％的样本作为训练样本。

步骤7)用PCA滤波器提取样本特征：

步骤7a)对im₁进行逐像素有重叠划分，得到由大小为s×s的多个图像块组成的集合I¹，同时对im₂进行逐像素有重叠划分，得到大小为s×s的多个图像块组成的集合I²：

其中，s为大于等于1的奇数，pos表示图像块序号，pos＝1,2,...,T，T表示SAR图像im₁或im₂中图像块的个数，与图像中像素个数相等，T＝W·H,本实发明施例中s的取值为13；

步骤7b)将训练样本集和不确定样本集合并为样本集，并从I¹中选择与样本集中样本位置x_a相对应的图像块同时从I²中选择与样本集中样本位置x_a相对应的图像块a表示样本序号，且a＝1,2,...,N_E，N_E为样本集中样本的总数，其中正样本个数为N_T，不确定样本中样本个数为N_U,则有N_E＝N_T+N_U；

步骤7c)对进行逐像素有重叠划分，得到由大小为m×m的多个邻域块P_i¹，同时对进行逐像素有重叠划分，得到大小为m×m的多个邻域块P_i²，并将邻域块P_i¹和P_i²合并成的邻域块P_i组合到集合P_temp中：

其中，m为大于等于1小于等于s的奇数，i表示邻域块序号，且i＝1,2,...,n，n表示邻域块个数，且n＝s·s·N_E，本发明实施例中m的取值为3；

步骤7d)对P_temp中的邻域块P_i进行去均值后再进行向量化，得到包括由组成的向量矩阵P：

步骤7e)计算PP^T的特征向量，并将PP^T前N₁个特征向量作为第一层PCA滤波器W_l¹：

W_l¹＝mat(q_l(PP^T))∈R^2m×m>1

其中，N₁为正整数，q_l(PP^T)表示PP^T的第l个特征向量，mat(v)表示把向量映射成为一个矩阵W∈R^2m×m，本实施例中N₁的取值为8，表示上、下、左、右、左上、左下、右上、右下8个方向；

步骤7f)通过第一层PCA滤波器W_l¹计算样本P_i的特征信息P_i^l：

P_i^l＝P_i*W_l¹

其中，*表示三维卷积操作；

步骤7g)对第一层滤波器W_l¹输出的特征信息P_i^l进行去均值后再向量化，得到包括由组成的向量矩阵Q：

步骤7h)计算QQ^T的特征向量，并将QQ^T前N₂个特征向量作为第二层PCA滤波器

其中，q_k(QQ^T)表示QQ^T的第k个特征向量，k＝1,2,...,N₂，N₂为正整数，mat(v)表示把向量映射成为一个矩阵W∈R^2m×m，本实施例中N₂的取值为8，表示上、下、左、右、左上、左下、右上、右下8个方向；

步骤7i)通过第二层PCA滤波器计算特征信息P_i^l的二次特征信息

步骤7j)采用赫维赛德阶跃函数对特征信息矩阵进行二值化，计算公式为：

将得到的二值化信息矩阵转换为每个位置的数值为的矩阵O后，计算P_i的特征f_i，并将作为训练样本特征，作为不确定样本特征，不确定样本特征个数为s·s·N_U，其中：

其中，表示对的直方图统计；

步骤8)训练支持向量机并对不确定样本分类：

将训练样本特征作为支持向量机的输入，得到训练好的支持向量机，并将不确定样本特征作为训练好的支持向量机的输入，得到SAR图像的变化检测结果。

以下结合仿真实验，对本发明的技术效果作进一步详细的分析。

1、实验条件

以下仿真是在Core i5-6500CPU，主频3.2GHz，内存4GB的硬件环境和MATLAB2015b的软件环境下进行的。

2、实验内容

实验一：使用本发明与现有的利用PCANet对SAR图像进行变化检测方法对两时相Bern数据Real SAR图像进行变化检测，其结果如图2所示。

实验二：使用本发明与现有的利用PCANet对SAR图像进行变化检测方法对两时相黄河入海口数据Real SAR图像进行变化检测，其结果如图3所示。

实验三：使用本发明与现有的利用PCANet对SAR图像进行变化检测方法对两时相Ottawa数据Real SAR图像进行变化检测，其结果如图4所示。

3、实验结果和分析

参照图2，其中图2(a)是第一时刻Bern数据集SAR图像，图2(b)是第二时刻Bern数据SAR图像，图2(c)是变化检测结果参考图，图2(d)是现有技术的变化检测结果，图2(e)是本发明的变化检测结果。

参照图3，其中图3(a)是第一时刻黄河入海口数据SAR图像，图3(b)是第二时刻黄河入海口数据SAR图像，图3(c)是变化检测结果参考图，图3(d)是现有技术的变化检测结果，图3(e)是本发明的变化检测结果。

参照图4，其中图4(a)是第一时刻Ottawa数据SAR图像，图4(b)是第二时刻Ottawa数据SAR图像，图4(c)是变化检测结果参考图，图4(d)是现有技术的变化检测结果，图4(e)是本发明的变化检测结果。

为验证变化检测结果质量，本发明选取总体精度、Kappa系数和执行时间作为性能指标参数，评价本发明的检测准确性和计算效率。

该实验结果如表1和表2所示，表1给出现有技术和本发明检测准确性的比较，表2给出现有技术和本发明计算效率的比较。

表1现有技术与本发明检测精度比较

表2现有技术与本发明效率比较

由图2、图3、和图4可以看到，本发明得到的变化检测结果与参考图给出结果更加接近。由表1可以看到，本发明在检测精度上优于现有技术。这是由于在获取训练样本过程中，本发明采用以邻域块的灰度值为聚类特征，使训练样本更可靠，用一致选取法消除了样本不均衡，使提取的正负训练样本特征更加完整，利用基于块的处理方式充分考虑像素邻域信息并有效抑制了相干斑噪声，使提取的特征能够更好的反映变化样本和不变化样本的区别，因此本发明的变化检测准确度明显高于现有技术。

由表2可以看到，本发明计算效率更高。这是由于本发明从显著性区域中选取训练样本和不确定样本，有效减少了样本数量，因此本发明计算效率高于现有技术。