一种两时相遥感图像变化检测方法与流程

本发明涉及一种有效的两时相遥感图像变化检测方法，属于图像处理技术领域。

背景技术：

遥感图像变化检测是指通过比较分析同一地理区域不同时相的遥感图像，获得该地区这段时间的地物变化信息。近几年，遥感图像变换检测技术已广泛应用于众多领域，如环境监测、农业研究，自然灾害评估、森林植被变化监测等。

目前常用的遥感图像变化检测方法主要是基于遥感图像的纹理特征、边缘特征或形状特征等特征来构造差异图像，不同类型的特征侧重描述图像不同方面的细节信息。但是仅仅依靠单一类型特征进行变化检测往往具有局限性，也容易在运算中丢失有用信息。因此，本发明综合利用遥感图像的灰度值信息和纹理特征信息，通过叠加两时相遥感图像的灰度差值图和纹理特征统计量差值图构造差异图像。叠加后的差异图像，综合考虑了遥感图像灰度值信息和纹理特征信息，改善了利用单一特征构造差异图像时信息量不足的缺点，为后续的检测操作提供了更丰富的原始图像信息。

模糊c均值聚类是众多模糊聚类算法中，应用最广泛且效果较好的聚类算法之一。模糊c均值算法不再将待分类的对象严格的划分到某一特定类别中，而是运用概率的理论来描述对象隶属于不同类别的程度。这种运用模糊理论的软划分，更能客观的描述现实世界，从而得到了广泛的应用。

模糊c均值聚类方法的主要原理是通过迭代使目标函数最小化来得到最优解。在每次迭代运算中，需要计算所有样本点的隶属度，并更新聚类中心，收敛速度慢。同时，模糊c均值聚类方法的聚类效果依赖于初始值的选取，当初始聚类中心选取不当时，也会影响算法的收敛速度，结果容易陷入局部极小值，难以获得全局最优解的。考虑到模糊c均值聚类方法的这两点不足，本发明在模糊c均值方法的基础上，提出了一种改进的方法来检测差异图像中的变化区域和未变化区域。

技术实现要素：

发明目的：针对现有技术中存在的问题，本发明提供一种有效的两时相遥感图像变化检测方法，该方法基于特征融合和模糊c均值(fuzzyc-means，fcm)，有效结合了不同特征差异图像的信息，提高了抗噪性能和变化检测的精确度；同时，进一步提高了fcm算法的收敛速度。

技术方案：一种两时相遥感图像变化检测方法，包括如下步骤：

步骤一：给定待检测的两时相遥感图像，设计一种融合灰度特征和纹理特征的方法来构造差异图像。

(1)对待检测的两时相遥感图像i1，i2(其中i1为第一时相图像，i2为第二时相图像，大小均为m×n)分别进行辐射校正、几何校正等预处理操作；

(2)算图像i1，i2的灰度值差值图x0：

x0(i,j)＝|i1(i,j)-i2(i,j)|；(1≤i≤m,1≤j≤n)

其中，i1(i,j)和i2(i,j)为两时相图像i1和i2对应像素点的灰度值；

(3)针对遥感图像i1和i2，提取它们的纹理特征。本发明在综合考虑算法计算量和检测效果等因素的基础上，选取了4种纹理特征统计量，即能量、对比度、相关性和熵进行计算。

(4)分别计算图像i1和i2的4种纹理特征统计量的差值矩阵。

为了使灰度值差值图x0和纹理特征差值图d的数值有相同的取值区间，需要对纹理特征差值图d进行归一化处理，以得到归一化后的纹理特征差值图d'：

其中，dmax为矩阵d中的最大值；

(5)融合灰度值差值图x0和纹理特征差值图d'，得到最终差异图像x：

差异图像x包含了两时相遥感图像i1和i2的灰度值信息和纹理特征信息，它将作为后续改进的模糊c均值方法的输入图像，以得到最终的变化检测结果。

步骤二：提出了一种快速模糊c均值方法，对融合后的差异图像进行二分类(变化类和非变化类)。

(1)初始化模糊c均值的参数，模糊加权指数m＝2，最大迭代次数t＝30，停止迭代的最小误差ε＝0.0001，初始迭代次数t＝1，聚类总数c＝2；

(2)以差异图像x作为改进模糊c均值方法的输入图像，像素点的总个数为n，c1和c2分别表示图像x中的变化类和未变化类。v1和v2分别表示两类的聚类中心，初始化聚类中心：

v1⁰＝xmax

v2⁰＝xmin

其中，xmax和xmin分别为差异图像x灰度值的最大值和最小值。选择距离比最小的c1个点，这些像素点到变化类中心的距离明显小于到未变化类中心的距离，用于表示属于变化类的点的个数；选择距离比最大的c2个点，这些像素点到变化类中心的距离明显大于到未变化类中心的距离，用于表示属于未变化类的点的个数。

(3)计算所有像素点的隶属度矩阵uik(i＝1,...,c)：

其中，xk表示第k个像素点，vi是第i类的聚类中心，vj是第j类的聚类中心，聚类总数c＝2,n是图像像素的总个数。

(4)分别计算差异图像x中所有像素点到两类c1和c2的聚类中心v1^t-1和v2^t-1的欧式距离m1和m2，并计算所有像素点到两类中心的距离比k：

m1＝{m11,m12,…m1n}，m2＝{m21,m22,…m2n}

(5)对距离比k中所有的值进行排序，并选择距离比最小的c1个点，这些像素点到变化类中心的距离明显小于到未变化类中心的距离，则这些像素点有最大的概率属于变化类，记作xk∈c1，修改这些点的隶属度u1k＝1,u2k＝0；选择距离比最大的c2个点，这些像素点到变化类中心的距离明显大于到未变化类中心的距离，则这些像素点有最大的概率属于未变化类，记作xk∈c2，修改这些点的隶属度u1k＝0,u2k＝1；

(6)计算变化类中心v1^t和未变化类中心v2^t：

并更新第(5)步中选出的c1个和c2个点的灰度值，分别为各自聚类中心值；

(7)若t＞t(t表示预设的迭代次数的阈值，t表示已经进行的迭代次数)，或者(其中，停止迭代的最小误差ε＝0.0001)且则停止迭代，得到最终隶属度矩阵uik(i＝1,2)；否则t＝t+1，增加c1和c2的值，并跳转到第(3)步；

(8)根据最终的隶属度矩阵uik(i＝1,2)，判断所有像素点的类别属性：当u1k＞u2k(k＝1,2,…,n)，则像素点k属于变化类，否则像素点k属于未变化类；

(9)输出最终变化检测结果。

本发明采用上述技术方案，具有以下有益效果：

(1)本发明的方法采用融合两时相图像灰度特征和纹理特征的方法构造差异图像，解决了基于单一类型特征构造差异图像时信息量不足的问题

(2)本方法在模糊c均值算法的迭代过程中，通过修改与类中心距离最近的像素点的隶属度来加快算法收敛速度。另外，采用差异图像中灰度值最高和最低的点作为变化类和非变化类的初始聚类中心，进一步提高了模糊c均值算法的收敛速度。

附图说明

图1为本发明实施例的方法流程图。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图1所示，两时相遥感图像变化检测方法，包括如下步骤：

(1)对待检测的两时相遥感图像i1，i2(其中i1为第一时相图像，i2为第二时相图像，大小均为m×n)分别进行辐射校正、几何校正等预处理操作；

(2)算图像i1，i2的灰度值差值图x0：

x0(i,j)＝|i1(i,j)-i2(i,j)|；(1≤i≤m,1≤j≤n)

其中，i1(i,j)和i2(i,j)为两时相图像i1和i2对应像素点的灰度值；

(3)针对遥感图像i1和i2，提取它们的纹理特征。本发明在综合考虑算法计算量和检测效果等因素的基础上，选取了4种纹理特征统计量，即能量、对比度、相关性和熵进行计算：

能量：

其中，p(i,j)为归一化后的灰度共生矩阵，(i,j)表示矩阵像素坐标。能量体现了图像中纹理的均匀性和紧密性，当灰度共生矩阵中的值相差越大时，能量值就越大。

对比度：

对比度描述了图像中像素值的亮度反差情况，视觉上表现为图像的清晰度以及图像色彩的鲜艳程度。

相关性：

其中，μx，μy，分别表示和的均值和标准差。相关性描述了矩阵中各行各列的相似程度，体现了图像纹理的一致性。

熵：

熵描述了一幅图像的平均信息量，体现了图像中纹理的复杂程度。

以第一时相图像i1为例，将m×n大小的窗口在图像上滑动，计算每个窗口内的4种纹理特征统计量，作为窗口内中心像素点的纹理特征值。每次移动的距离为1个像素点，依次进行0度扫描，45度扫描，90度扫描和135度扫描。当窗口在整幅图像上移动结束，得到4个扫描方向的纹理数据并叠加，最终得到图像i1的能量矩阵d11，对比度矩阵d12，相关性矩阵d13和熵矩阵d14，大小均为m×n。同理得到图像i2的能量矩阵d21，对比度矩阵d22，相关性矩阵d23和熵矩阵d24；

(4)分别计算图像i1和i2的4种纹理特征统计量的差值矩阵。以能量矩阵为例，计算图像i1，i2的能量差值矩阵d1：

d1(i,j)＝|d11(i,j)-d21(i,j)|；(1≤i≤m,1≤j≤n)

同理得到对比度差值矩阵d2，相关性差值矩阵d3和熵差值矩阵d4，叠加4种特征统计量差值矩阵，得到纹理特征差值矩阵d。为了使灰度值差值图x0和纹理特征差值图d的数值有相同的取值区间，需要对纹理特征差值图d进行归一化处理，以得到归一化后的纹理特征差值图d'：

其中，dmax为矩阵d中的最大值；

(5)融合灰度值差值图x0和纹理特征差值图d'，得到最终差异图像x：

差异图像x包含了两时相遥感图像i1和i2的灰度值信息和纹理特征信息，它将作为后续改进的模糊c均值方法的输入图像，以得到最终的变化检测结果。

(6)初始化模糊c均值的参数，模糊加权指数m＝2，最大迭代次数t＝30，停止迭代的最小误差ε＝0.0001，初始迭代次数t＝1，聚类总数c＝2，c1和c2的初始值根据输入图像的不同来选择；

(7)以差异图像x作为改进模糊c均值方法的输入图像，像素点的总个数为n，c1和c2分别表示图像x中的变化类和未变化类。v1和v2分别表示两类的聚类中心，初始化聚类中心：

v1⁰＝xmax

v2⁰＝xmin

其中，xmax和xmin分别为差异图像x灰度值的最大值和最小值；

(8)计算所有像素点的隶属度矩阵uik(i＝1,2)：

(9)分别计算差异图像x中所有像素点到两类c1和c2的聚类中心v1^t-1和v2^t-1的欧式距离m1和m2，并计算所有像素点到两类中心的距离比k：

m1＝{m11,m12,…m1n}，m2＝{m21,m22,…m2n}

(10)对距离比k中所有的值进行排序，并选择距离比最小的c1个点，这些像素点到变化类中心的距离明显小于到未变化类中心的距离，则这些像素点有最大的概率属于变化类，记作xk∈c1，修改这些点的隶属度u1k＝1,u2k＝0；选择距离比最大的c2个点，这些像素点到变化类中心的距离明显大于到未变化类中心的距离，则这些像素点有最大的概率属于未变化类，记作xk∈c2，修改这些点的隶属度u1k＝0,u2k＝1；

(11)计算变化类中心v1^t和未变化类中心v2^t：

并更新步骤(10)中选出的c1个和c2个点的灰度值，分别为各自聚类中心值；

(12)若t＞t，或者且则停止迭代，得到最终隶属度矩阵uik(i＝1,2)；否则t＝t+1，增加c1和c2的值，并跳转到第(8)步；

(13)根据最终的隶属度矩阵uik(i＝1,2)，判断所有像素点的类别属性：当u1k＞u2k(k＝1,2,…,n)，则像素点k属于变化类，否则像素点k属于未变化类；

(14)输出最终变化检测结果。