一种结合DT‑CWT和MRF的遥感图像变化检测方法与装置与流程

本发明涉及图像处理技术，特别涉及一种结合dt-cwt和mrf的遥感图像变化检测方法与装置。

背景技术：

随着遥感数据获取手段不断提高和更新周期不断缩短，变化检测技术作为遥感图像处理与分析中的一项重要应用，已广泛应用于土地利用和土地覆盖变化监测，森林或植被变化检测，灾害救援与治理、城镇和道路变化、地理数据库更新等众多领域。由于无监督变化检测方法除了利用两时相的原始图像外，不需要任何额外的先验知识就能得到最终变化检测结果，方法较为高效，故成为研究的热点。通常无监督变化检测是基于像素差异的变化检测，可分为阈值法、变换法、图像建模法等。阈值法是通过选择合适的阈值来区分差异图像中的变化类与非变化类像素，从而得到最终的变化区域，具有简单直观的优点，但是往往只进行单一尺度的分割，难免会引起过度分割或不完全分割，同时检测结果严重依赖于选取的阈值。变换法主要是多尺度变换法方法，如小波变换，这类方法可以进行多尺度分析，能较好地克服了传感器噪声、配准误差等因素的影响，但没有考虑像素间的相关性，检测边缘粗糙，误检像素较多。图像建模法是通过建立模型来模拟空间上下文相关的分布信息，如马尔科夫随机场(markovrandomfield，mrf)充分考虑了邻域内像素之间的相互影响，较好地克服了孤立像素问题，但突变和细节信息表达较弱，伪变化检测较多。为了解决上述方法存在的问题，已经提出了将小波变换与mrf相结合的遥感图像变化检测方法，这些方法利用了小波变换的时频分析能力并且考虑了像素的空间相关性，较好地克服了单一像素孤立性、噪声和配准误差等因素引起的伪变化，但是由于在融合过程中直接舍弃了高频信息，只取低频信息进行试验，难免造成信息的损失，加之小波变换存在着方向选择性和平移敏感性，使得变化检测细节仍表达不足。

技术实现要素：

本发明旨在提供一种结合双树复小波变换(dual-treecomplexwavelettransform，dt-cwt)和马尔科夫随机场(mrf)的遥感图像变化检测方法与装置，其不仅能有效平衡高频信息的保留和噪声的去除，而且边缘检测更加平滑，变化检测结果具有较好的区域一致性，提高遥感图像变化检测的精度，可充分地消除由于现有技术的限制和缺陷导致的一个或多个问题。

本发明另外的优点、目的和特性，一部分将在下面的说明书中得到阐明，而另一部分对于本领域的普通技术人员通过对下面的说明的考察将是明显的或从本发明的实施中学到。通过在文字的说明书和权利要求书及附图中特别地指出的结构可实现和获得本发明目的和优点。

本发明提供了一种结合dt-cwt和mrf的遥感图像变化检测方法，其特征在于，所述方法具体包括以下步骤：

步骤1，输入大小相同的经过精确配准的两期同区域的遥感影像对应的两幅灰度图像x1、x2；

步骤2，求所述两幅灰度图像x1、x2的差值图像；

步骤3，对所述差值图像进行3层dt-cwt分解；

步骤4，利用基于fcm的mrf分割算法，对经dt-cwt分解后得到的高频子带进行变化特征提取；

步骤5，对通过步骤3的dt-cwt分解得到的低频分量和通过步骤4的变化特征提取得到的高频分量进行重构；

步骤6，利用所述基于fcm的mrf分割算法，对通过步骤5重构后的结果进行变化检测；

步骤7，基于相似度融合规则，对通过步骤6得到的检测结果进行融合；

步骤8，输出融合后的变化检测结果。

优选的，步骤2具体包括：

根据公式d(i,j)＝|x2(i,j)-x1(i,j)|求所述差值图像d，其中，d(i,j)为差值图像d在位置(i，j)处的像素值，x1(i,j)和x2(i,j)分别是灰度图像x1和x2的在位置(i，j)处的像素值，i，j分别为1≤i≤p，1≤j≤q的整数，p为所述输入图像的行数，q为所述输入图像的列数。

优选的，步骤3具体包括：采用两组并行的基于正交变换的q-shift(quartersampleshift)实数滤波器组来对步骤2中得到的差值图像d进行3层dt-cwt分解。第一层分解采用相差一个单位延迟的奇数长度的滤波器组，分解结果的一树作为第一层分解复小波系数的实部，另一树作为第一层分解复小波系数的虚部；第二、三层分解均采用迟近似为1/4采样间隔的偶数长度的滤波器组，分解结果一树作为第二、三层分解复小波系数的实部，另一树作为第二、三层分解复小波系数的虚部。

优选的，步骤4具体包括以下子步骤：

步骤4.1，对所述经过3层dt-cwt分解得到的高频子带hk,l，(k＝1,2,3；l＝±15°,±45°,±75°)中的第一层高频分量h1,l(l＝±15°,±45°,±75°)置零，即得到第一层高频变化特征提取结果，记为

步骤4.2，利用基于bayes阈值的fcm算法对所述经过3层dt-cwt分解得到的高频子带hk,l(k＝1,2,3；l＝±15°,±45°,±75°)中的第二、三层高频分量hk,l(k＝2,3；l＝±15°,±45°,±75°)进行初始聚类分割，得到初始分割结果

步骤4.3，基于mrf与grf的等价性，利用条件迭代模型icm对通过步骤4.2获得的初始分割结果进行迭代更新，得到最终的最优分割，记为

步骤4.4，对通过步骤4.3得到的最终分割结果进行二进制掩膜，得到掩膜图bk,l(k＝2,3；l＝±15°,±45°,±75°)；

步骤4.5，将所述掩膜图bk,l(k＝2,3；l＝±15°,±45°,±75°)矩阵像素值为1的记为变化类c，值为0的记为未变化类u；

步骤4.6，将c类像素对应位置的二、三层高频子带hk,l(k＝2,3；l＝±15°,±45°,±75°)的像素值保持不变，u类像素对应位置的二、三层高频子带的像素值置零，得到提取变化特征后的二、三层高频子带

优选的，步骤4.2具体包括以下子步骤：

步骤4.2.1，根据最小错误概率准则，求自适应bayes阈值level；

步骤4.2.2，根据v1＝(min(hk,l(i,j))+level)/2、v2＝(level+max(hk,l(i,j))/2，分别对所述经过3层dt-cwt分解得到的第二、三层高频分量hk,l，(k＝2,3；l＝±15°,±45°,±75°)求初始聚类中心v1、v2，记vm＝[v1,v1],其中m为迭代次数，level为自适应bayes阈值，max(hk,l(i,j))和min(hk,l(i,j))分别表示每个子带hk,l，(k＝2,3；l＝±15°,±45°,±75°)中的像素最大和最小值，hk,l(i,j)表示位置(i,j)处的像素值，1≤i≤p，1≤j≤q的整数，p为hk,l(i,j)对应的行数，q为hk,l(i,j)对应的列数；

步骤4.2.3，根据式

对每个子带hk,l(k＝2,3；l＝±15°,±45°,±75°)计算隶属度矩阵u，其中，c为聚类中心数，n为权值，||hj-cj||为样本hj到第i类聚类中心的距离度量，||hj-ck||为样本hj到第k类聚类中心的距离度量；

步骤4.2.4，根据步骤4.2.3中得到的隶属度矩阵u，求初始标记场，并更新步骤4.2.2的聚类中心vm，m为迭代次数；

步骤4.2.5，重复步骤4.2.3、4.2.4不断对所述标记场进行迭代更新，直到达到迭代终止条件(优选的，迭代终止条件为15次)；

步骤4.2.6，输出迭代最终的标记场，也即是初始分割结果

优选的，步骤4.3具体包括以下子步骤：

步骤4.3.1，根据步骤4.2的fcm初始分割结果计算特征场参数：均值μm和方差

步骤4.3.2，利用步骤4.3.1计算的特征场参数μm和方差来计算当前分割的每个像素的特征场能量，并计算当前分割的每个像素的标记场能量；

步骤4.3.3，根据特征场和标记场总能量之和最小值所在变化或未变化的类别来更新当前分割结果，并记为m为迭代次数；

步骤4.3.4，重复步骤4.3.1-4.3.3，不断对进行迭代更新，直到达到循环终止条件(优选的，循环终止条件设置为30次)，获得最终的mrf-map最优分割结果m为迭代次数。

步骤4.3.5，输出所述最终的mrf-map最优分割结果

优选的，在步骤5中，对步骤3中dt-cwt分解得到的低频分量{l1,l2,l3}和步骤4中变化特征提取得到的高频分量进行重构，其中，步骤5具体包括：采用与步骤3中相同的两组并行的基于正交变换的q-shift(quartersampleshift)实数滤波器组对步骤3中dt-cwt分解得到的低频{l1,l2,l3}和步骤4经过变化特征提取的高频分量进行重构(即逆变换)。通过对l1和进行逆变换得重构后的第一层结果i1，通过对l2和进行逆变换得重构后的第一层结果i2，通过对l3和进行逆变换得重构后的第一层结果i3。

优选的，步骤6具体包括以下子步骤：

步骤6.1，根据v1＝(min(ia(i,j))+level)/2、v2＝(level+max(ia(i,j))/2，对重构后的结果i1，i2，i3分别求初始聚类中心v1、v2，记vm＝[v1,v1]，其中m表示迭代次数，level为自适应bayes阈值，a＝1,2,3，max(ia(i,j))和min(ia(i,j))分别表示i1，i2或i3的像素最大和最小值，1≤i≤p，1≤j≤q的整数，p为ia(i,j)对应的行数，q为ia(i,j)对应的列数；

步骤6.2，根据式

对i1，i2和i3计算隶属度矩阵u，其中c为聚类中心数，n为权值，||hj-cj||为样本hj到第i类聚类中心的距离度量，||hj-ck||为样本hj到第k类聚类中心的距离度量；

步骤6.3，根据步骤6.2中得到的隶属度矩阵u，分别求初始标记场并更新步骤6.2的聚类中心vm，m表示迭代次数；

步骤6.4，判断迭代终止条件(优选的，迭代终止条件设置为15次)，重复步骤6.2、6.3不断对标记场进行迭代更新，得到最终迭代的标记场，也即是初始分割结果w1^m、m表示迭代次数；

步骤6.5，对初始分割结果w1^m、由和分别求均值μm和方差即为特征场的参数，其中a＝1,2,3，p为输入图像的行数，q为输入图像的列数，i、j分别为1≤i≤p，1≤j≤q的整数，m表示迭代次数；

步骤6.6，采用条件迭代模型icm算法，由能量函数(只考虑二阶邻域系统，其中c为势团，c为二阶势团集合，vc为势函数，(w1,w2)表示二阶邻域系统，m为迭代次数)计算当前分割每个像素的标记场能量um(w)，由能量函数计算当前分割每个像素的特征场能量um(w,f)；

步骤6.7，根据能量最小原则em＝min(um(w)+um(w,f))，即根据迭代终止后特征场和标记场总能量之和最小值所在变化或未变化类来更新当前分割结果，更新后的分割结果记为m表示迭代次数；

步骤6.8，判断循环终止条件(优选的，迭代终止条件设置为30次)，重复步骤6.5、6.6和6.7不断对当前分割进行迭代更新，达到最终的mrf-map最优分割结果也即是三层变化检测结果，并记为z1，z2，z3，其中m表示迭代次数。

优选的，步骤7具体包括以下子步骤：

步骤7.1，统计比较通过步骤6获得的变化检测结果z2、z3与z1中对应位置像素值是否相同，①若z1(i,j)≠z2(i,j)且z1(i,j)≠z3(i,j)，则记u＝0；②若z1(i,j)＝z2(i,j)≠z3(i,j)或z1(i,j)＝z3(i,j)≠z2(i,j)，则记u＝0.5；③若z1(i,j)＝z2(i,j)＝z3(i,j)，则记u＝1，其中z1(i,j)、z2(i,j)、z3(i,j)表示z1、z2、z3对应(i,j)位置的像素值，u表示相似度；

步骤7.2，判断相似度u是否满足u≥0.5，若满足，保持z1(i,j)值不变，若不满足，取z2像素值z2(i,j)对z1像素值z1(i,j)进行更新，更新后的结果即为融合后的变化检测结果。

本发明还提供了一种结合dt-cwt和mrf的遥感图像变化检测的装置，所述装置包括计算机存储介质和处理器，其中，所述计算机存储介质中存有计算机程序，所述处理器被配置为执行所述计算机程序，以执行根据上述任一实施例所述的方法。

本发明有助于图像信息的多方向表达和多尺度分析，充分利用了像素间的相关性，不仅能有效平衡高频信息的保留和噪声的去除，而且边缘检测更加平滑，变化检测结果具有较好的区域一致性，大大降低了误检率和去除配准误差的影响，极大提高了遥感图像变化检测的精度。

附图说明

图1为根据本发明实施例的、结合dt-cwt和mrf的遥感图像变化检测方法的流程图。

图2为根据本发明实施例的、3层dt-cwt分解树状图。

图3为根据本发明实施例的、高频变化特征提取流程图。

图4为根据本发明实施例的、基于自适应bayes阈值的fcm分割流程图。

图5为根据本发明实施例的、icm分割流程图。

图6为根据本发明实施例的、基于相似度的变化检测结果融合流程图。

具体实施方式

下面对本发明进行更全面的描述，其中说明本发明的示例性实施例。

如图1所示，本发明提出的结合双树复小波变换(dual-treecomplexwavelettransform，dt-cwt)和马尔科夫随机场(mrf)的遥感图像变化检测方法包括以下步骤：

步骤1，输入大小相同的经过精确配准的两期同区域的遥感影像对应的两幅灰度图像x1、x2；

其中，精确配准可由软件enviclassic下registration模块并通过选择同名控制点来完成，配准精度小于一个像素，本发明进行的变化检测是基于像素的两期的遥感图像变化检测，故在进行变化检测之前要对图像进行精确的配准，避免因配准引起的误检率上升。根据本发明的优选实施例，所述两幅灰度图像x1、x2的大小均为p×q，其中，p为所述输入图像(即，所述灰度图像)的行数，q为所述输入图像(即，所述灰度图像)的列数。

步骤2，求所述两幅灰度图像x1、x2的差值图像d；

具体的，采用差值法，根据公式d(i,j)＝|x2(i,j)-x1(i,j)|求所述差值图像d，其中，d(i,j)为差值图像d在位置(i，j)处的像素值，x1(i,j)和x2(i,j)分别是灰度图像x1和x2的在位置(i，j)处的像素值，i，j分别为1≤i≤p，1≤j≤q的整数，p为所述输入图像的行数，q为所述输入图像的列数。

发明人利用差值法、比值法、对数法求取差值图像进行了对比实验，发现差值法求取差异图像效果较好，因此，本发明优选使用了差值法，通过求取差值图像来为后续变化检测提供数据准备。

步骤3，对所述差值图像进行dt-cwt分解；

步骤3具体包括：采用两组并行的基于正交变换的q-shift(quartersampleshift)实数滤波器组来对步骤2中得到的差值图像d进行3层dt-cwt分解。第一层分解采用相差一个单位延迟的奇数长度的滤波器组，分解结果的一树作为第一层分解复小波系数的实部，另一树作为第一层分解复小波系数的虚部；第二、三层分解均采用迟近似为1/4采样间隔的偶数长度的滤波器组，分解结果一树作为第二、三层分解复小波系数的实部，另一树作为第二、三层分解复小波系数的虚部。

图2示出了3层dt-cwt分解树状图，x为原始图像，x(n)(n＝1，2，3)为第n层分解的高、低频子带，其中，阴影部分为低频子带，空白部分为六个方向的高频子带(±15°,±45°,±75°)。

dt-cwt的分解原理具体可描述为：第一层分解两树(a树和b树)均使用相差一个单位延迟的奇数长度的滤波器，在高层分解时使用群延迟近似为1/4采样(优选的采用群延迟近为1/4采样)的偶数滤波器组，并且设定a树1/4时延，b树为3/4时延，以便保证两树之间所要求的1/2采样周期的时延，这样做的好处是：b树的采样点刚好位于a树中因二抽取隔点采样而舍弃的采样值的点位(这样就等价于没进行二抽取，而且非常易于实现)。对于上述3层dt-cwt的分解，具体可以包括：分解的第一层两树均使用相差一个单位延迟的奇数长度的双正交滤波器，在第二、三层分解时使用群延迟近似为1/4采样(优选的采用群延迟近为1/4采样)的偶数滤波器组，并且设定一树为1/4时延，一树为3/4时延，以便保证两树之间所要求的1/2采样周期的时延。

另外，dt-cwt是有一定冗余的，对于一维信号，产生了2:1的冗余度，对于n维信号，将产生2n:1的冗余度。为了克服奇/偶(odd/even)滤波器因严格线性相位引起的采样结构差异的缺点，本发明采用基于正交变换的q-shift(quartersampleshift)滤波器组对图像进行二维dt-cwt分解，每层分解可以得到一个低频子带和6个方向(±15°,±45°,±75°)的高频子带，所述低频子带记为{l1,l2,l3}，所述6个方向的高频子带记为hk,l，(k＝1,2,3；l＝±15°,±45°,±75°)，其中，k表示分解层数，l表示子带的方向。由于dt-cwt不仅保持了传统离散小波变换dwt(discretewavelettransform)良好的时频分析能力，而且还具有近似的平移不变性、良好的方向选择性、有限的数据冗余等特点，解决了复小波不能完全重构的问题，故能够较好的表达图像的细节信息。

步骤4，利用基于fcm的mrf分割算法，对经dt-cwt分解后得到的高频子带进行变化特征提取；

具体来说，步骤4是对步骤3中经过3层dt-cwt分解得到的高频子带hk,l，(k＝1,2,3；l＝±15°,±45°,±75°)进行变化特征提取，具体实现方法见图3详解。经过多尺度分解后的遥感图像的低频分量包含图像的主要信息，高频分量包含主要的突变信息(包括细节信息和噪声)，由于噪声绝大部分存在于分解的第一层高频中，且第二层和第三层高频信息同样包含了图像的大部分高频信息，本发明通过舍去第一层高频信息并利用mrf模型分割算法预先提取出其他层高频域的变化特征，来达到在去除噪声的同时较好地保留高频细节信息的目的；其中，如图3所示，步骤4具体包括以下步骤：

步骤4.1，对所述经过3层dt-cwt分解得到的高频子带hk,l，(k＝1,2,3；l＝±15°,±45°,±75°)中的第一层高频分量h1,l(l＝±15°,±45°,±75°)置零，作为第一层高频特征提取结果其中，k表示分解层数，l表示子带的方向；

步骤4.2，利用基于bayes阈值的fcm算法(具体实现方法见图4详解)对所述经过3层dt-cwt分解得到的高频子带hk,l，(k＝1,2,3；l＝±15°,±45°,±75°)中的第二、三层高频分量hk,l(k＝2,3；l＝±15°,±45°,±75°)进行初始聚类分割，得到初始分割结果其中m表示迭代次数；

如图4所示，步骤4.2具体包括以下子步骤：

步骤4.2.1，根据最小错误概率准则，求自适应bayes阈值level；

步骤4.2.2，根据v1＝(min(hk,l(i,j))+level)/2、v2＝(level+max(hk,l(i,j))/2，对所述经过3层dt-cwt分解得到的第二、三层高频分量hk,l，(k＝2,3；l＝±15°,±45°,±75°)分别求初始聚类中心v1、v2，记vm＝[v1,v1],其中m为迭代次数，level为自适应bayes阈值，max(hk,l(i,j))和min(hk,l(i,j))分别表示每个子带hk,l，(k＝2,3；l＝±15°,±45°,±75°)中的像素最大和最小值；

步骤4.2.3，根据式

对每个子带hk,l(k＝2,3；l＝±15°,±45°,±75°)计算隶属度矩阵u，其中，c为聚类中心数，n为权值，||hj-cj||为样本hj到第i类聚类中心的距离度量，||hj-ck||为样本hj到第k类聚类中心的距离度量，m表示迭代次数；

步骤4.2.4，根据步骤4.2.3中得到的隶属度矩阵u，求初始标记场并更新步骤4.2.2的聚类中心vm,其中m为迭代次数；

步骤4.2.5，重复步骤4.2.3、4.2.4不断对标记场进行迭代更新，直到达到迭代终止条件(优选的，迭代终止条件为15次)；

步骤4.2.6，输出最终迭代的标记场，也即为初始分割结果

步骤4.3，基于mrf与grf的等价性，利用条件迭代模型(icm)算法对4.2.6中的初始分割结果进行迭代更新(优选的，icm最大迭代次数设置为30次)，得到最终的最优分割，记为

如图5所示，步骤4.3具体包括以下子步骤：

步骤4.3.1，根据步骤4.2的fcm分割结果计算特征场参数：均值μm和方差

具体为，对基于bayes阈值的fcm分割结果由和分别求特征场参数μm和其中，m为迭代次数，μm为每个子带变化与未变化类的均值，为每个子带变化与未变化类的方差，p为所述输入图像的行数，q为所述输入图像的列数，i、j分别为1≤i≤p，1≤j≤q的整数；

步骤4.3.2，利用步骤4.3.1计算的特征场参数μm和来计算当前分割的每个像素的特征场能量，并计算当前分割的每个像素的标记场能量；

具体为，对基于bayes阈值的fcm分割结果采用条件迭代模型icm算法，由能量函数(只考虑二阶邻域系统，其中，c为势团，c为二阶势团集合，vc为势函数，(w1,w2)表示二阶邻域系统)计算当前分割的每个像素的标记场能量um(w)，由能量函数(其中，s为所有像素集合，f为特征场，m为迭代次数)计算当前分割的每个像素的特征场能量um(w,f)；

步骤4.3.3，根据能量最小原则来更新当前分割结果，更新后的结果记为其中m为迭代次数；

具体的，根据能量最小原则em＝min(um(w)+um(w,f))，即,根据特征场和标记场总能量之和最小值所在变化或未变化类来更新当前分割结果，并记为m为迭代次数；

步骤4.3.4，重复步骤4.3.1-4.3.3，不断对当前分割进行迭代更新，直到达到循环终止条件(迭代终止条件设置为30次)，获得最终的mrf-map最优分割结果m为迭代次数。

步骤4.3.5，输出所述最终的mrf-map最优分割结果m为迭代次数。

步骤4.4，对通过步骤4.3得到的最终分割结果进行二进制掩膜，得到掩膜图bk,l(k＝2,3；l＝±15°,±45°,±75°)，其中m为迭代次数；

步骤4.5，将所述掩膜图bk,l(k＝2,3；l＝±15°,±45°,±75°)矩阵像素值为1的记为变化类c，值为0的记为未变化类u；

步骤4.6，将c类像素对应位置的二、三层高频子带hk,l(k＝2,3；l＝±15°,±45°,±75°)的像素值保持不变，u类像素对应位置的二、三层高频子带的像素值置零，得到提取特征后的二、三层高频分量

步骤5，对通过步骤3的dt-cwt分解得到的低频分量{l1,l2,l3}和通过步骤4的变化特征提取得到的高频分量进行重构；

具体的，采用与步骤3中相同的两组并行的基于正交变换的q-shift(quartersampleshift)实数滤波器组，对步骤3中dt-cwt分解得到的低频{l1,l2,l3}和图3步骤1，图4步骤6经过变化特征提取的高频分量进行重构(即逆变换)得到重构后的三层结果i1，i2，i3，这里运用重构的思想来使得各层分量具有和原图同样地大小，避免了因尺度间插值引起的误差，为后续进行不同尺度间各分量的融合做准备。

步骤6，利用所述基于fcm的mrf分割算法，对通过步骤5重构后的三层结果i1，i2，i3进行变化检测；

本发明选择基于dt-cwt与mrf模型相结合的变化检测方法(同步骤4中的基于fcm初始分割的icm分割方法，即对重构后的三层结果，先进行基于自适应bayes阈值的fcm初始分割，再进行基于fcm的icm最终分割)，利用dt-cwt的多方向表达性、各向异性和多尺度特性，增强图像信息的表达与分析，提高了变化检测的精度，利用mrf能对像素间的相关性较好地描述，来降低检测的误检率和消除配准误差的影响，提高鲁棒性。

其中，步骤6具体包括：

步骤6.1，根据v1＝(min(ia(i,j))+level)/2、v2＝(level+max(ia(i,j))/2，分别对重构后的结果i1，i2，i3求初始聚类中心v1、v2，记vm＝[v1,v1],其中m为迭代次数，level为自适应bayes阈值，a＝1,2,3，max(ia(i,j))和min(ia(i,j))分别表示i1，i2或i3的像素最大和最小值，i、j分别为1≤i≤p，1≤j≤q的整数，p为ia(i,j)的行数，q为ia(i,j)的列数；

步骤6.2，根据式

对i1，i2和i3分别计算隶属度矩阵u，其中c为聚类中心数，n为权值，||hj-cj||为样本hj到第i类聚类中心的距离度量，||hj-ck||为样本hj到第k类聚类中心的距离度量,m表示迭代次数；

步骤6.3，根据步骤6.2中得到的隶属度矩阵u，分别求初始标记场并更新步骤6.2的聚类中心vm，m表示迭代次数；

步骤6.4，判断迭代终止条件(迭代终止条件设置为15次)，重复步骤6.2、6.3不断对标记场进行迭代更新，得到最终迭代的标记场，也即是初始分割结果w1^m、m表示迭代次数；

步骤6.5，对初始分割结果w1^m、由和分别求均值μm和方差即为特征场参数，其中a＝1,2,3；p为输入图像的行数，q为输入图像的列数，i、j分别为1≤i≤p，1≤j≤q的整数,m表示迭代次数；

步骤6.6，采用条件迭代模型icm算法，由能量函数(只考虑二阶邻域系统，其中c为势团，c为二阶势团集合，vc为势函数，(w1,w2)表示二阶邻域系统)计算当前分割每个像素的标记场能量um(w)，由能量函数(其中s为所有像素集合，s属于s，f为特征场，m为迭代次数)计算当前分割每个像素的特征场能量um(w,f)；

步骤6.7，根据能量最小原则em＝min(um(w)+um(w,f))，即根据迭代终止后特征场和标记场总能量之和最小值所在变化或未变化类来更新当前分割结果，更新后的分割结果记为m表示迭代次数；

步骤6.8，判断循环终止条件(迭代终止条件设置为30次)，重复步骤6.5、6.6和6.7不断对当前分割进行迭代更新，达到最终的mrf-map最优分割结果也即是三层变化检测结果，并记为z1，z2，z3，其中m表示迭代次数。

步骤7，基于相似度融合规则，对通过步骤6得到的检测结果进行融合；

具体的，基于相似度融合规则，对步骤6中的三层变化检测结果z1，z2，z3进行融合，来进一步提高变化检测精度，经试验显示，融合后的结果边缘检测更加平滑细腻，检测杂点较少，具有较高的变化检测精度和较强的鲁棒性；其中，步骤7具体包括：

图6示出了基于相似度的变化检测结果融合过程，包括：

步骤7.1，统计比较通过步骤6获得的变化检测结果z2、z3与z1中对应位置像素值是否相同，①若z1(i,j)≠z2(i,j)且z1(i,j)≠z3(i,j)，则记u＝0；②若z1(i,j)＝z2(i,j)≠z3(i,j)或z1(i,j)＝z3(i,j)≠z2(i,j)，则记u＝0.5；③若z1(i,j)＝z2(i,j)＝z3(i,j)，则记u＝1，其中u表示相似度；

步骤7.2，判断相似度u是否满足u≥0.5，若满足，保持z1(i,j)值不变，若不满足，取z2像素值z2(i,j)对z1像素值z1(i,j)进行更新，更新后的结果即为融合后的变化检测结果。

步骤8，输出融合后的变化检测结果，即为最终的变化检测结果。

本发明还提出了一种结合dt-cwt和mrf的遥感图像变化检测的装置，所述装置包括计算机存储介质和处理器，其中，所述计算机存储介质中存有计算机程序，所述处理器被配置为执行所述计算机程序，以执行根据上述任一实施例所述的方法。

以上内容仅为本发明的较佳实施例，对于本领域的普通技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，本说明书内容不应理解为对本发明的限制。