CN106203521B - 基于差异图自步学习的sar图像变化检测方法 - Google Patents

Abstract

本发明提出了一种基于差异图自步学习的SAR图像变化检测方法，用于解决现有SAR图像变化检测方法结果正确率低和受噪声影响大的问题，其实现步骤是：1.输入两幅校准后待检测的SAR图像，并获得差异图DI；2.对差异图像进行模糊C均值聚类，得到隶属度矩阵U；3.根据隶属度矩阵得到初始变化检测结果图R₁；4.利用结果图R₁构建候选样本集X'；5.采用均匀选择从X'中生成最终样本集X，并进行自步学习得到分类器；6.使用分类器获得中间的变化检测结果图R₂，7.根据R₂利用局部邻域信息获得最终变化检测结果图R。本发明提高了变化检测结果的正确率，减少了噪声影响，可用于灾难评估，环境监测，地貌变化预测。

Description

基于差异图自步学习的SAR图像变化检测方法

技术领域

本发明属于SAR图像处理技术领域，特别涉及一种SAR图像的变化检测方法，可用于灾难评估，环境监测，地貌变化预测。

背景技术

由于合成孔径雷达SAR具有全天时，全天候，拍摄范围大的成像特点，对于SAR图像处理的应用与需求不断提高，而SAR图像变化检测是SAR图像处理技术中的一个重要分支，在军事和民用方面都有着重要的应用。SAR图像变化检测是从两幅在相同地点不同时间拍摄得到的SAR图像中，根据某些特征的差异程度，识别出发生变化的区域和没有发生变化的区域的过程。

随着国内外SAR图像技术不断发展，研究人员相继提出了许多有效的变化检测方法，应用最广泛的是先产生差异图再对差异图进行分析的思路。产生差异图主要分为差值方法和比值方法，差值方法是早期获得差异图的方法，该方法运算简单，容易实现，但对噪声十分敏感。比值方法可以有效得抑制SAR图像中富含的乘性噪声，是目前广泛应用的一种方法。对差异图进行分析类似于图像分割，主要分为阈值方法和聚类方法。阈值方法是在差异图中建立一个数学模型，得到一个合适的阈值，该方法运算简单，容易实现，但结果的正确率很低。聚类方法是在差异图中得到两个聚类中心，将所有像素分为两类，模糊C均值聚类是应用最广泛的聚类算法，该方法利用模糊理论对图像进行分割，能够保留图像的细节信息，获得良好的结果，但是该方法没有利用邻域信息，受噪声影响很大，鲁棒性低。

随着现代社会对SAR图像变化检测要求的提高，传统的各种算法存在结果正确率低，对噪声的鲁棒性差的缺点，使得传统的各种方法所得到的变化检测结果已不能满足要求，研究人员针对这些缺陷做出了一些改进，例如对传统的模糊C均值聚类(FCM)进行的改进：

Stelios Krinidis等人在IEEE Transactions On Image Process,19(2010)1328-1337上发表了论文“A robust fuzzy local information C-means clusteringalgorithm”，该论文首先随机初始化模糊隶属度矩阵，计算聚类中心，利用邻域的空间和灰度信息计算每个像素点的邻域信息，对每个像素点的模糊隶属度进行修正，再计算聚类中心，并多次迭代，输出最终的聚类结果。该方法通过利用邻域的空间和灰度信息减少了噪声的影响，提高了结果的正确率。然而对于一些受到噪声污染的SAR图像，该方法得到的结果依然受到严重的斑点噪声影响，产生了较低的分类准确率。

发明内容

本发明的目的在于针对上述现有技术存在的不足，提出一种基于差异图自步学习的SAR图像变化检测方法，以减小斑点噪声的影响，提高分类的准确率。

为实现上述目的，本发明的技术方案包括如下：

(1)输入两幅校准后待检测的SAR图像I₁和I₂，并利用对数率操作器获得这两幅图像的差异图DI；

(2)对差异图DI进行模糊C均值聚类，获得该差异图DI的隶属度矩阵U；

(3)根据隶属度矩阵U，对差异图DI中每个像素进行分类，得到初始的变化检测结果图R₁；

(4)根据初始变化检测结果图R₁，选择满足的像素点，构成候选样本集X'，其中，t表示像素点邻域内与中心像素类标一致的像素点数目，s表示像素点的邻域范围，α表示选择阈值，取值为0.7；

(5)从候选样本集X'中均匀选择出数目为0.1×N的样本，得到各类样本数目相同的最终样本集X，其中N为差异图DI中的像素总个数；

(6)根据最终样本集X，使用逻辑斯谛回归进行自步学习，得到分类器的参数w_end；

(7)根据所获取的参数w_end，对差异图DI中每一个像素进行分类，得到中间的变化检测结果图R₂；

(8)根据中间的变化检测结果图R₂，利用局部邻域信息重新对差异图DI中的每一个像素进行分类，得到最终的变化检测结果图R：

8a)根据中间的变化检测结果图R₂，统计每个像素点邻域内类标分别为0和1的像素点数目m₁和m₂；

8b)判断m₁是否大于m₂，若是，则将像素的类标设置为0，否则，将像素的类标设置为1；

8c)将得到的类标作为差异图DI中各个像素点的灰度值，生成最终的变化检测结果图R，并通过imshow函数显示出来。

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

1.本发明由于在生成样本集时，采用均匀选择的策略，保持两个类的样本数目相同，样本的特征采用差异图中的空间邻域信息，与现有技术随机选择样本的算法相比，有效提升了样本的多样性，提高了样本集的质量，减少了噪声对结果的影响，同时提高了分类器的准确率。

2.本发明由于在训练分类器时，采用了自步学习框架，与现有技术使用全部样本学习分类器的方法相比，本发明使训练过程从易到难得进行，减少了随机初始化参数对分类器得影响，有效避免了训练过程陷入局部最优解，提升了训练结果的稳定性，同时进一步提高了分类的准确率。

3.本发明由于在得到最终的变化检测结果图像时，利用了图像的邻域相似性，与现有技术直接使用分类器的分类结果的方法相比，有效得减少了噪声对图像分类结果的影响，提升了对噪声的鲁棒性，同时进一步提高了变化检测结果的准确率。

附图说明

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

图2是用本发明和现有FCM，MRFFCM、FLICM三个SAR图像变化检测方法对Bern数据集进行变化检测的结果图。

具体实施方式

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

参照图1，本发明的实现步骤如下：

步骤1，输入两幅校准后待检测的SAR图像I₁和I₂，并利用对数率操作器获得这两幅图像的差异图DI；

输入的两幅图像是任意相同尺寸且在同一地点拍摄的两幅SAR图像，在本实施例中，采用的是在瑞士伯尔尼地区拍摄的尺寸为301×301的SAR图像数据集；

使用对数率操作器，即将比值图像变换到对数尺度，将乘性噪声转变为加性噪声，以减少相干斑噪声对变化检测结果的影响，较好得保留图像的细节信息，具体实现步骤如下：

1a)根据输入图像计算两幅图像中相同位置像素点灰度值的比值对数差异值dif＝{dif₁,…,dif_i,…,dif_N}，其中N表示每幅图像中像素点的总个数，N＝90601，SAR图像I₁和I₂中第i个位置像素点灰度值的比值对数差异值i＝1,2,…,N，x_1i表示输入图像I₁中第i个像素点的灰度值，x_2i表示输入图像I₂中第i个像素点的灰度值，|·|为取绝对值运算符，ln(·)表示取自然对数运算符；

1b)将得到的相同位置像素点灰度值的比值对数差异值dif_i作为灰度值，生成和输入图像尺寸相同的差异图DI。

步骤2，对差异图DI进行聚类，获得该差异图DI的隶属度矩阵U。

现有技术中对图像的聚类方法有模糊C均值聚类，k均值聚类，FCM-S等，由于模糊C均值聚类方法实现简单，对噪声有较高的鲁棒性，因此本实例采用模糊C均值聚类方法实现对差异图DI的聚类，其实现步骤如下：

2a)随机初始化差异图DI的聚类中心为V＝{V₁,V₂,…,V_c}，其中c＝2表示聚类数目，设初始迭代次数为1；

2b)根据聚类中心V，计算当前迭代中差异图DI的每个像素点对每个聚类中心的隶属度u_ij，将所有像素点的隶属度组成隶属度矩阵U，其中隶属度u_ij利用下式得到：

其中，u_ij表示差异图DI中第i个像素对第j个聚类中心的隶属度，u_ij满足约束条件：m＝2表示模糊指数；d_ij和d_il分别表示差异图DI中第i个像素与第j个和第l个聚类中心的欧式距离；

2c)根据隶属度矩阵U和聚类中心V_j，计算当前迭代的目标函数值J，并得到相邻两次迭代中目标函数值J的改变量δ：

2c1)根据隶属度矩阵U和聚类中心V，计算当前迭代的目标函数值J：

其中N表示差异图DI中的像素总个数；

2c2)根据获得的当前迭代的目标函数值J以及上次迭代的目标函数值J₁，利用δ＝|J-J1|得到相邻两次迭代中目标函数值J的改变量δ；

2d)判断改变量δ是否小于收敛阈值ε＝0.00001，若是，则输出当前的隶属度矩阵U，否则，迭代次数加1，并采用梯度下降法，获得下一次迭代的聚类中心V＝{V₁,…,V_j,…,V_c}，返回步骤2b)，其中下一次迭代的第j个类的聚类中心V_j的计算为：

其中，x_i表示差异图DI中第i个像素点的灰度值。

步骤3，根据隶属度矩阵U，对差异图DI中每个像素进行分类，得到初始的变化检测结果图R₁。

3a)根据得到的隶属度矩阵U，找出差异图DI中每个像素所在列中的最大隶属度，并将最大隶属度在隶属度矩阵U中的行标号作为对应像素的类标；

3b)将所获得的差异图DI中每个像素的类标作为灰度值，生成初始的结果图像R₁。

步骤4，构建候选样本集X'。

4a)设置选择样本的条件其中，t表示像素点邻域内与中心像素类标一致的像素点数目，s表示像素点的邻域范围，取值为s＝3，α表示选择阈值，取值为α＝0.7；

4b)对初始变化检测结果图R₁中的每一个像素点，判断其是否满足若是，则将该像素点包含入候选样本集X'，否则，不加入。

按照条件选择出来的像素点同时具备准确性和多样性，提升了样本集的质量，从而增加样本集的丰富性，提高分类器的正确率。

步骤5，从候选样本集X'中均匀选择出数目为0.1×N的样本，得到各类样本数目相同的最终样本集X。

本步骤采用均匀选择和使用差异图作为样本特征，以减少样本集中各个类别样本数目的差异，增加各类样本的多样性，进一步提高分类的正确率，其实现步骤为：

5a)判断候选样本集X'中两个类的像素点数目n₁和n₂是否有差异，若是，则复制候选样本集X'中数目小的一类的像素点加入候选样本集X'，使n₁和n₂一致，更新候选样本集X'，否则，执行步骤5b)；

5b)在更新后的候选样本集X'中随机选择num＝0.1×N＝9060个像素点，作为最终的样本；

5c)根据得到的9060个像素点，利用初始的变化检测结果图R₁和差异图DI，生成最终的样本集X：

5c1)利用9060个像素点和差异图DI，得到第i个像素点的特征f_i，

f_i＝{f_i1,…,f_ik,…,f_ip},f_ik∈Q_i,k＝1,2,…,p，i＝1,2,…,9060

其中，Q_i表示差异图DI中以第i个像素点为中心像素的邻域，尺寸为3，p＝3×3＝9表示邻域中像素点的数目，f_ik表示邻域Q_i中第k个像素点在差异图DI中的灰度值；

5c2)将9060个像素点在初始的变化检测结果图R₁中对应的灰度值设为各个像素点的特征所对应的类标；

5c3)将得到的9060个像素点的特征和类标作为最终样本，输出最终的样本集X。

步骤6，根据最终样本集X，使用逻辑斯谛回归进行自步学习，得到分类器的参数w_end。

现有技术中分类器的模型有感知器，贝叶斯分类器，逻辑斯谛回归等，由于逻辑斯谛回归利用逻辑斯谛函数和极大似然估计，可以产生性能良好的分类器。因此本实例采用逻辑斯谛回归的分类器模型。

现有技术中对分类器的学习有普通直接学习，强化学习，自步学习等，由于自步学习可以有效得避免梯度下降方法受初始值影响大和容易陷入局部最优解的问题。因此本实例采用自步学习来得到分类器，从而提高了梯度下降方法找到最优解的可能性，增加了结果的稳定性，提升了分类器的正确率，其实现步骤如下：

6a)随机初始化分类器的参数w，初始自步迭代次数设为1，设定初始阈值λ＝0.1；

6b)根据分类器的参数w和阈值λ，选择对数似然函数值小于λ的样本构成当前自步迭代中参与训练的样本集A：

6b1)利用当前自步迭代下的参数w，计算样本集X中每个样本的对数似然函数值l_i：

l_i＝y_i×ln(g_i)+(1-y_i)×ln(1-g_i)

其中，y_i为第i个样本的类标，ln(·)为取自然对数运算符，g_i为第i个样本属于第一个类的概率，计算如下：

其中，e^(·)为指数运算符，(·)^T为取转置运算符；

6b2)判断l_i是否小于阈值λ，若是，则将该样本加入当前自步迭代的样本集A中，否则不加入；

6c)根据当前自步迭代的参数w和样本集A，利用逻辑斯谛回归和梯度下降法，得到当前自步迭代的新参数w^*：

6c1)根据当前自步迭代的参数w和样本集A，设定初始梯度下降迭代次数为1，得到当前自步迭代的总对数似然函数l：

其中h为当前自步迭代的样本集A中样本的数目，∑(·)为求和运算符；

6c2)利用梯度下降法优化总对数似然函数l，获得当前梯度下降迭代的新的权值w_g：

其中，a为梯度下降的步长，

6c3)判断当前梯度下降迭代次数是否大于最大梯度下降迭代次数t＝50，若是，令w^*＝w_g，输出获得的参数w^*，否则，令参数w＝w_g，梯度下降迭代次数加1，返回步骤6c1)；

6d)判断当前自步迭代次数是否大于最大自步迭代次数T＝15，若是，则将当前自步迭代的新参数w^*作为最终的分类器参数w_end输出，否则，令参数w＝w^*，当前自步迭代次数加1，利用λ＝λ×β增大阈值λ，开始下一次迭代，返回步骤6b)，其中，β＝1.1为自步学习的步长。

步骤7，根据最终的分类器参数w_end，对差异图DI中每一个像素进行分类，得到中间的变化检测结果图R₂。

7a)根据差异图DI，获得所有像素点的特征：

f_i＝{f_i1,…,f_ik,…,f_ip},f_ik∈Q_i,k＝1,2,…,p,i＝1,2,…,N

其中f_i表示差异图DI中第i个像素点的特征。

7b)根据所获取的权值w和像素点的特征，得到每一个像素点属于第一个类的概率，其中第i个像素点属于第一个类的概率P_i为：

7c)判断第i个像素点属于第一个类的概率P_i是否大于0.5，若是，则将测试像素点分为第一类，类标设为0，否则，将像素点分为第二类，类标设为1；

7d)将得到的每个像素的类标作为灰度值，生成中间的结果图像R₂。

步骤8，根据中间的变化检测结果图R₂，利用局部邻域信息重新对差异图DI中的每一个像素进行分类，得到最终的变化检测结果图R。

利用图像邻域信息，可以减少噪声对结果的影响，增强对噪声的鲁棒性，进一步提升分类的准确率，实现步骤如下：

8a)根据中间的变化检测结果图R₂，统计每个像素点邻域内类标分别为0和1的像素点数目m₁和m₂；

8b)判断m₁是否大于m₂，若是，则将像素的类标设置为0，否则，将像素的类标设置为1；

8c)将得到的类标作为差异图DI中各个像素点的灰度值，生成最终的变化检测结果图R，并通过imshow函数显示出来。

以下结合附图和仿真实验结果，对本发明的技术效果作进一步描述：

1.仿真条件：

本发明采用Matlab R2013软件在配置为core i3 2.30GHZ，内存4GB，WINDOWS 7系统的计算机上进行。

2、仿真内容：

用本发明和现有FCM，MRFFCM、FLICM三个SAR图像变化检测方法对Bern数据集进行变化检测，其结果如图2所示，其中：

图2(a)为输入待处理的第一幅SAR图像I₁；

图2(b)为输入待处理的第二幅SAR图像I₂；

图2(c)为变化检测的真实结果图像；

图2(d)为FCM方法的变化检测结果图像；

图2(e)为MRFFCM方法的变化检测结果图像；

图2(f)为FLICM方法的变化检测结果图像；

图2(g)为本发明方法的变化检测结果图像。

由图2可以看出：与其它对比方法相比，本发明方法的变化检测结果图像与真实结果图像更加接近，错分点少，边缘清晰，区域一致性高，分类正确率高，对噪声的鲁棒性较好，证明了本发明方法可以得到很好的变化检测结果。

将本发明和其他三种方法的结果图即图2(d),(e),(f)和(g)与真实的变化检测结果图即图2(c)进行对比，统计每种方法的结果图的FP，FN，TN和TP值，其中：

FP表示在图2(c)中没有变化类的像素被所用方法错误分类的像素点个数；

FN表示在图2(c)中变化类的像素被所用方法错误分类的像素点个数；

TP表示在图2(c)中变化类的像素被所用方法正确分类的像素点个数；

TN表示在图2(c)中没有变化类的像素被所用方法正确分类的像素点个数。

用表示所用方法的分类错误率，表示所用方法的分类正确率，计算综合评价指标为KC，可以更好得评价变化检测结果的优劣，计算如下：

其中Nc表示图2(c)中属于变化类的像素点数目，Nu表示图2(c)中属于没有变化类的像素点数目。

将本发明和FCM，MRFFCM、FLICM方法在Bern数据集上的变化检测结果的FP，FN，OE，PCC和KC值进行统计，如表1所示。其中FP，FN和OE值越低表示结果越好，PCC和KC值越高表示结果越好。

表1不同方法对Bern数据集的变化检测结果

方法	FP	FN	OE％	PCC	KC
						FCM	123	350	0.52	0.9948	0.7703
MRFFCM	63	281	0.38	0.9962	0.8337
						FLICM	36	293	0.36	0.9964	0.8379
本发明	131	154	0.31	0.9969	0.8738

从表1中可以看出，对于同一组SAR图像数据集，FCM算法的分类结果最差，其它对比方法的结果相比FCM算法的结果好很多，但是都有较高的FN。而本发明方法有最低的FN值，最低的错误率OE和最高的正确率PCC，以及最高的Kappa系数，比其他三个方法的Kappa系数高3个百分点以上，表明采用本发明方法对SAR图像进行变化检测可以取得理想的结果。

Claims (2)

1.一种基于差异图自步学习的SAR图像变化检测方法，包括：

(1)输入两幅校准后待检测的SAR图像I₁和I₂，并利用对数率操作器获得这两幅图像的差异图DI；

(2)对差异图DI进行模糊C均值聚类，获得该差异图DI的隶属度矩阵U；

(3)根据模糊隶属度矩阵U，对差异图DI中每个像素进行分类，得到初始的变化检测结果图R₁；

(4)根据初始变化检测结果图R₁，选择满足的像素点，构成候选样本集X'，其中，t表示像素点邻域内与中心像素类标一致的像素点数目，s表示像素点的邻域范围，α表示选择阈值，取值为0.7；

(5)从候选样本集X'中均匀选择出数目为0.1×N的像素点，并利用初始的变化检测结果图R₁和差异图DI，生成各类样本数目相同的最终样本集X，其中N为差异图DI中的像素总个数：

5a)判断候选样本集X'中两个类的像素点数目n₁和n₂是否有差异，若是，则复制候选样本集X'中数目小的一类的像素点，使n₁和n₂一致，更新候选样本集X'，否则，执行步骤5b)；

5b)在更新后的候选样本集X'中随机选择num＝0.1×N个像素点，作为最终的样本，其中N＝90601；

5c)根据得到的num个像素点，利用初始的变化检测结果图R₁和差异图DI，生成最终的样本集X：

5c1)利用num个像素点和差异图DI，得到第i个像素点的特征f_i，

f_i＝{f_i1,…,f_ik,…,f_ip},f_ik∈Q_i,k＝1,2,…,p,i＝1,2,…,num

其中，Q_i表示差异图DI中以第i个像素点为中心像素的邻域，尺寸为3，p＝3×3＝9表示邻域中像素点的数目，f_ik表示邻域Q_i中第k个像素点在差异图DI中的灰度值；

5c2)将num个像素点在初始变化检测结果图R₁中对应的灰度值设为各个像素点的特征所对应的类标；

5c3)将得到的num个像素点的特征和类标作为最终样本，输出最终的样本集X；

(6)根据最终样本集X，使用逻辑斯谛回归进行自步学习，得到分类器的参数w_end：

6a)随机初始化分类器的参数w，设初始自步迭代次数为1，设定初始阈值λ＝0.1；

6b)根据分类器的参数w和阈值λ，选择对数似然函数值小于λ的样本构成当前自步迭代中参与训练的样本集A:

6b1)利用当前自步迭代下的参数w，计算样本集X中每个样本的对数似然函数值l_i：

l_i＝y_i×ln(g_i)+(1-y_i)×ln(1-g_i)，

其中，y_i为第i个样本的类标，ln(·)为取自然对数运算符，g_i为第i个样本属于第一个类的概率，计算如下：

其中，e^(·)为指数运算符，(·)^T为取转置运算符；

6b2)判断l_i是否小于阈值λ，若是，则将该样本加入当前自步迭代的样本集A中，否则，不加入；

6c)根据当前自步迭代的参数w和样本集A，得到当前自步迭代的新参数w^*：

6c1)根据当前自步迭代的参数w和样本集A，设定初始梯度下降迭代次数为1，得到当前自步迭代的总对数似然函数l：

其中h为当前自步迭代的样本集A中样本的数目，∑(·)为求和运算符；

6c2)利用梯度下降法优化总对数似然函数l，获得当前梯度下降迭代的新的权值w_g：

其中，a为梯度下降的步长，

6c3)判断当前梯度下降迭代次数是否大于最大梯度下降迭代次数t＝50，若是，则令当前自步迭代的新参数w^*＝w_g，否则，令参数w＝w_g，梯度下降迭代次数加1，返回步骤6c1)；

6d)判断当前自步迭代次数是否大于最大自步迭代次数T＝15，若是，则将当前自步迭代的新参数w^*作为最终的分类器参数w_end输出，否则，令参数w＝w^*，当前自步迭代次数加1，利用λ＝λ×β增大阈值λ，开始下一次迭代，返回步骤6b)，其中，β为自步学习的步长；

(7)根据所获取的参数w_end，对差异图DI中每一个像素进行分类，得到中间的变化检测结果图R₂；

(8)根据中间的变化检测结果图R₂，利用局部邻域信息重新对差异图DI中的每一个像素进行分类，得到最终的变化检测结果图R：

8a)根据中间的变化检测结果图R₂，统计每个像素点邻域内类标分别为0和1的像素点数目m₁和m₂；

8b)判断m₁是否大于m₂，若是，则将像素的类标设置为0，否则，将像素的类标设置为1；

8c)将得到的类标作为差异图DI中各个像素点的灰度值，生成最终的变化检测结果图R，并通过imshow函数显示出来。

2.根据权利要求1所述基于差异图自步学习的SAR图像变化检测方法，其中步骤(2)对差异图DI进行模糊C均值聚类，按如下步骤进行：

2a)随机初始化差异图DI的聚类中心为V，设初始迭代次数为1；

2b)根据聚类中心V，计算当前迭代中差异图DI的每个像素点对每个聚类中心的隶属度u_ij，将所有像素点的隶属度组成隶属度矩阵U；

2c)根据隶属度矩阵U和聚类中心V，计算当前迭代的目标函数值J，并得到相邻两次迭代中目标函数值J的改变量δ；

2d)判断改变量δ是否小于收敛阈值ε＝0.001，若是，则输出当前的隶属度矩阵U，否则，迭代次数加1，并采用梯度下降法，获得下一次迭代的聚类中心V，返回步骤2b)。