CN103810704B - 基于支持向量机和判别随机场的sar图像变化检测方法 - Google Patents

Abstract

本发明属于SAR图像变化检测技术领域，公开了基于支持向量机和判别随机场的SAR图像变化检测方法。该基于支持向量机和判别随机场的SAR图像变化检测方法包括以下步骤：对原始两时相图像进行灰度值归一化处理，并在处理后的图像中提取对应的灰度特征差值和纹理特征差值；组成差值特征向量；利用加权平均比率算子提取差值图像中每个像素点的边界强度；在差值图像中选取训练样本，将训练样本用对应的差值特征向量进行表示，通过训练支持向量机，得到测试样本的初始分类标签、以及测试样本的分类标签的后验概率；得出初始的支持向量机—判别随机场模型；更新支持向量机—判别随机场模型，得出对应的测试样本的最终分类标签以及变化检测结果。

Description

基于支持向量机和判别随机场的SAR图像变化检测方法

技术领域

本发明属于SAR图像变化检测技术领域，特别涉及基于支持向量机和判别随机场的SAR图像变化检测方法。

背景技术

随着合成孔径雷达(synthetic aperture radar,SAR)技术的逐步成熟和SAR图像分辨率的不断提高，SAR图像的使用逐渐为人们所重视。同光学遥感图像相比，SAR图像不受天气、云层等因素的影响，可以全天候、全天时获得遥感数据，是较好的变化检测信息源。

SAR图像变化检测(change detection)通过对不同时期SAR图像的比较分析，根据图像之间的差异分析来获取所需要的地物变化信息。变化检测技术可以应用于很多方面，例如对地震区域的定位和灾害评估；对农作物生长状况的监测；城区土地使用的监测等等，它在环境、农业、水利和军事等国民经济诸多领域都有着非常广泛的应用。

SAR图像变化检测方法一般可分为：基于直接比较法的变化检测方法，如图像差值法、图像比值法；基于间接比较法的变化检测方法，如特征提取后比较法、分类后比较法；基于多元变量分析的变化检测方法，如主成分分析法。近期在SAR图像变化检测上研究比较多的有：基于统计模型的变化检测方法，如L.Bruzzone对SAR图像对数比值差异图进行广义高斯建模(GGD,Generalized Gaussian Distributions)然后在GGD下应用改进的KI(Kittler–Illingworth)门限选择算法进行分类以获得最终变化影像，该方法取得了较好的检测结果但没有考虑空间信息，且模型参数的选择也是一个难点；基于多尺度分析的变化检测方法，如Kai-Kuang Ma提出一种基于双树-复小波变换(DT-CWT，Dual-Tree Complex Wavelet Transform)的多尺度变化检测方法，它利用DT-CWT对对数比值图进行多尺度分解，但没有考虑到图像的纹理信息，阈值的选取也是一个棘手的问题；近年来新发展起来的是基于核方法的SAR图像变化检测算法，Gustavo Camps-Valls在2008年首先提出了将核方法应用于SAR图像变化检测，该方法首先提取图像的强度信息和纹理信息，然后构造强度纹理比值差值合成核(RDC_kernel)实现SAR图像变化检测，该方法可以有效的实现SAR图像变化检测，但是没有考虑空间信息，且对噪声比较敏感。

发明内容

本发明的目的在于提出基于支持向量机和判别随机场的SAR图像变化检测方法。该基于支持向量机和判别随机场的SAR图像变化检测方法该方法能够很好的结合SAR图像的强度特征和纹理特征，能充分考虑图像的空间信息，具有检测速度快、分类精度高的特点。

为实现上述技术目的，本发明采用如下技术方案予以实现。

基于支持向量机和判别随机场的SAR图像变化检测方法包括以下步骤：

S1：利用合成孔径雷达接收原始两时相图像，原始两时相图像包括第1时刻图像和第2时刻图像；然后分别对第1时刻图像和第2时刻图像进行灰度值归一化处理，得到第1时刻归一化图像X₁和第2时刻归一化图像X₂；所述第k时刻归一化图像X_k中第i行第j列的像素点表示为X'_k(i,j)，k取1和2，i取1至I，j取1至J，I为第1时刻归一化图像X₁的长度，J为第1时刻归一化图像X₁的宽度；提取X'_k(i,j)的灰度值g'_k(i,j)和X'_k(i,j)的纹理特征w'_k(i,j)；按照以下公式获得灰度特征差值△g(i,j)以及纹理特征差值△w(i,j)：△g(i,j)＝g'₂(i,j)-g'₁(i,j)，△w(i,j)＝w'₂(i,j)-w'₁(i,j)；然后将△g(i,j)和△w(i,j)组合成X'_k(i,j)的差值特征向量y(i,j)：y(i,j)＝{△g(i,j),△w(i,j)}；

S2：对第1时刻归一化图像X₁和第2时刻归一化图像X₂按照灰度值作差值运算，得到差值图像△X，利用加权平均比率算子提取差值图像中第t个像素点的边界强度r_t，t取1至M，M＝I×J；

S3：在差值图像中选取训练样本，将训练样本用对应的差值特征向量进行表示，通过训练支持向量机，得到测试样本的初始分类标签、以及测试样本的分类标签的后验概率；

S4：根据所述差值图像中每个像素点的边界强度、以及测试样本的后验概率，得出初始的支持向量机—判别随机场模型；

S5：根据所述初始分类标签和初始的支持向量机—判别随机场模型，更新支持向量机—判别随机场模型的相互势能函数，得出对应的测试样本的最终分类标签；根据所述对应的测试样本的最终分类标签，得出SAR图像的变化检测结果。

本发明的特点和进一步改进在于：

在步骤S1中，X'_k(i,j)的灰度值g'_k(i,j)为：

${g^{\′}}_k(i,j)=\frac{g_k(i,j)-\min \left(g_k\right)}{max\left(g_k\right)-\min \left(g_k\right)}$

其中，g_k(i,j)为所述第k时刻图像中第i行第j列的像素点的灰度值，min(g_k)为所述第k时刻图像中所有像素点的灰度值的最小值，max(g_k)为所述第k时刻图像中所有像素点的灰度值的最大值；

在步骤S1中，以X'_k(i,j)为中心像素点，建立对应的正方形像素窗口，所述正方形像素窗口的边长为η个像素点，η为大于1的奇数；则X'_k(i,j)的纹理特征w'_k(i,j)包括：对应的正方形像素窗口内像素点灰度值的均值μ'_k(i,j)、对应的正方形像素窗口内像素点灰度值的方差σ²'_k(i,j)、对应的正方形像素窗口内像素点灰度值的峰态ku'_k(i,j)、对应的正方形像素窗口内像素点灰度值的三阶矩sk'_k(i,j)、对应的正方形像素窗口内像素点灰度值的能量en'_k(i,j)、以及对应的正方形像素窗口内像素点灰度值的熵ent'_k(i,j)。

在步骤S2中，所述差值图像的边界强度包括差值图像中每个像素点的边界强度，定义平滑函数f(ρ)、因果滤波器f₁(ρ)和非因果滤波器f₂(ρ)：

$f\left(\mathrm{\ρ}\right)=\frac{1}{1+d}{f}_1\left(\mathrm{\ρ}\right)+\frac{d}{1+d}{f}_2(\mathrm{\ρ}-1)$

其中，f₁(ρ)＝cd^ρu(ρ)，f₂(ρ)＝cd^-ρu(-ρ)，d为设定常数且0<d<1，u(·)表示Heaviside函数，ρ为自变量；

然后，根据以下公式计算差值图像中第i行第j列的像素点△X(i,j)的边界强度|r_max(i,j)|为：

$\left|r_{max}(i,j)\right|=\sqrt{r_{\mathrm{Im}ax}^2(i,j)+r_{Jmax}^2(i,j)}$

$r_{Imax}(i,j)=max\{\frac{{\mathrm{\&mu;}}_{I1}(i-1,j)}{{\mathrm{\&mu;}}_{I2}(i+1,j)},\frac{{\mathrm{\&mu;}}_{I2}(i+1,j)}{{\mathrm{\&mu;}}_{I1}(i-1,j)}\}$

$r_{Jmax}(i,j)=max\{\frac{{\mathrm{\&mu;}}_{J1}(i,j-1)}{{\mathrm{\&mu;}}_{J2}(i,j+1)},\frac{{\mathrm{\&mu;}}_{J2}(i,j+1)}{{\mathrm{\&mu;}}_{J1}(i,j-1)}\}$

${\mathrm{\&mu;}}_{I1}(i,j)=f_1\left(i\right)*\left(f\right(j)\&CircleTimes;y(i,j\left)\right)$

${\mathrm{\&mu;}}_{I2}(i,j)=f_2\left(i\right)*\left(f\right(j)\&CircleTimes;y(i,j\left)\right)$

μ_J1(i,j)＝f₁(j)*(f(i)*y(i,j))

μ_J2(i,j)＝f₂(j)*(f(i)*y(i,j))

其中，*表示水平方向上的卷积，表示垂直方向上的卷积。

在步骤S3中，首先根据原始两时相图像中N组像素点，选取对应的N个有标签的训练样本；每组像素点包括：在第1时刻图像和第2时刻图像中处于相同位置的两个像素点；每个训练样本的标签的设置过程如下：通过对第1时刻图像和第2时刻图像进行观察对比，将所述N组像素点分为变化类像素点组和非变化类像素点组，根据N组像素点的分类情况来设置对应N个训练样本的标签；

将所述N个有标签的训练样本表示为其中，第s个有标签的训练样本表示为(x_s,l_s)，其中，s取1至N；x_s＝y(s_i,s_j)，s_i为第s个训练样本对应的像素点的横坐标，s_j为第s个训练样本对应的像素点的纵坐标；l_s表示第s个训练样本的分类标签，当第s个训练样本对应的一组像素点为变化类像素点组时，l_s＝1；当第s个训练样本对应的一组像素点为非变化类像素点组时，l_s＝0；

在支持向量机中建立如下C-SVC模型：

$\begin{array}{cc}\underset{\mathrm{\&alpha;}}{\min }& \frac{1}{2}{\mathrm{\&alpha;}}^{T}Q\mathrm{\&alpha;}-{\mathrm{\&Theta;}}^T\mathrm{\&alpha;}\end{array}$

s.t.l^Tα＝0

0≤α_s≤C,s＝1,...,N

其中，α＝[α₁,...,α_N]^T，α_s为待求的第s个训练样本对应的权重，Q为N×N维半正定矩阵，且Q中第p行第q列的元素Q_pq＝l_pl_qK(x_p,x_q)，p取1至N，q取1至N；K(x_p,x_q)为核函数，l＝[l₁,...,l_N]^T，Θ为N维列向量，Θ中的元素均为1；C和γ的取值交叉验证确定；

在步骤S3中，将差值图像中的每个像素点作为对应的一个测试样本；第t个像素点表示为(x'_t,l'_t)，t取1至M，M为差值图像中像素点的个数；x'_t＝y(t_i,t_j)，t_i为第t个像素点的横坐标，t_j为第t个像素点的纵坐标；l'_t表示待求的第t个像素点的分类标签；

根据所述C-SVC模型得出决策值

${\mathrm{dec}}_{{x^{\&prime;}}_t}=\underset{s=1}{\overset{N}{\mathrm{\&Sigma;}}}{l}_s{\mathrm{\&alpha;}}_{s}K(x_s,{x^{\&prime;}}_t)$

其中，在得出之后，得出第t个像素点的初始分类标签：当时，第t个像素点的初始分类标签为0；反之，第t个像素点的初始分类标签为1；

计算出第t个像素点的分类标签的后验概率p(l'_t|y(t_i,t_j))：

$p\left({l^{\&prime;}}_t\right|y\left(t_i,t_j\right))=\frac{1}{1+e^{{\mathrm{Adec}}_{{x^{\&prime;}}_t}+B}}$

其中，A和B通过以下公式确定：

$\underset{z=(A,B)}{min}-\underset{s=1}{\overset{N}{\mathrm{\&Sigma;}}}\left(t_{s}ln\right(p_s)+(1-t_s\left)ln\right(1-p_s\left)\right)$

$t_s=\left\{\begin{array}{cc}\frac{N_{+}+1}{N_{+}+2}& l_s=0\\ \frac{1}{N_{-}+2}& l_s=1\end{array}\right.,p_s=\frac{1}{1+e^{{\mathrm{Adec}}_{{x^{\&prime;}}_t}+B}}$

其中，N₊为：分类标签为1的训练样本的个数，N_-为：分类标签为0的训练样本的个数。

在步骤S4中，根据所述差值图像中每个像素点的边界强度，构建支持向量机—判别随机场模型中的相互势能函数I(l'_a,l'_b,r)：

$\begin{array}{c}I({l^{\&prime;}}_a,{l^{\&prime;}}_b,r)=\underset{(a,b)\&Element;N_H}{\mathrm{\&Sigma;}}\left(2\mathrm{\&delta;}\right({l^{\&prime;}}_a,{l^{\&prime;}}_b)-1)\&times;{\mathrm{\&alpha;}}_H\&times;\exp (-(r_a-r_b)/e\mathrm{deg}e\_C^2)\\ +\underset{(a,b)\&Element;N_V}{\mathrm{\&Sigma;}}\left(2\mathrm{\&delta;}\right({l^{\&prime;}}_a,{l^{\&prime;}}_b)-1)\&times;{\mathrm{\&alpha;}}_V\&times;\exp (-(r_a-r_b)/e\mathrm{deg}e\_C^2)\end{array}$

其中，(a,b)∈N_H的含义为：第a个像素点和第b个像素点水平相邻；(a,b)∈N_V的含义为：第a个像素点和第b个像素点垂直相邻，；edege_C为设定常数，α_H和α_V为I(l'_a,l'_b,r)的两个参数；当l'_a＝l'_b时，δ(l'_a,l'_b)＝1，反之，δ(l'_a,l'_b)＝0；

构建支持向量机—判别随机场模型中的联合势能函数A(l'_a,y(a_i,a_j))，A(l'_a,y(a_i,a_j))＝p(l'_a|y(a_i,a_j))，其中，p(l'_a|y(a_i,a_j))为第a个像素点的分类标签的后验概率；

构建初始的支持向量机—判别随机场模型p(l'|y,r)：

$p\left(l^{\&prime;}\right|y,r)=\frac{1}{Z}\exp \left(\underset{a\&Element;S}{\mathrm{\&Sigma;}}A\right({l^{\&prime;}}_a,y\left(a_i,a_j\right))+\underset{a\&Element;S}{\mathrm{\&Sigma;}}I({l^{\&prime;}}_a,{l^{\&prime;}}_b,r\left)\right)$

其中，Z为设定常数，S表示差值图像中所有像素点的集合；

将α_H和α_V用参数θ进行表示，即θ＝{α_H,α_V}，然后利用最小二乘法估计出θ的初始值θ₀。

步骤S5具体包括以下步骤：

S51：利用差值图像的每个像素点的初始分类标签构成原始标记场，设定k＝1；

S52：将第a个像素点的分类标签l'_a设为0，将θ的当前取值以及l'_a代入初始的支持向量机—判别随机场模型p(l'|y,r)中，计算出将第a个像素点的分类标签l'_a设为1，将θ的当前取值以及l'_a代入初始的支持向量机—判别随机场模型p(l'|y,r)中，计算出

将对应的分类标签作为第a个像素点更新后的分类标签；

S53：采用ICE迭代算法对参数θ和标记场进行更新，所述标记场指差值图像中每个像素点的分类标签；

S54：令k＝k+1，判断k是否小于K，K为设定值且K为大于1的自然数，如果k小于K，返回执行步骤S52；如果k＝K，则将当前的标记场作为最终的标记场，然后按照最终的标记场得出SAR图像的变化检测结果。

本发明的有益效果为：本发明在利用上下文信息的能力上、降噪上和检测精度上具有很大的优势。分类时无需数据降维，在检测速度方面有较高的性能，有效地减少了误分的出现，极大的提高了变化检测的检测精度。

附图说明

图1为本发明的基于支持向量机和判别随机场的SAR图像变化检测方法的流程示意图；

图2为本发明中最终输出分类标签的流程示意图；

图3为ICE迭代算法的流程示意图；

图4a为实测稻田受洪水灾害的第一时刻ERS-1SAR图像；

图4b为实测稻田受洪水灾害的第二时刻ERS-1SAR图像；

图4c为实测稻田受洪水灾害的变化检测参考图；

图4d为实测稻田受洪水灾害的RDC_Kernel变化检测结果示意图；

图4e为实测稻田受洪水灾害的本发明的变化检测结果示意图；

图5a为实测机场受洪水灾害的第一时刻JERS SAR channel1图像；

图5b为实测机场受洪水灾害的第二时刻JERS SAR channel1图像；

图5c为实测机场受洪水灾害的变化检测参考图；

图5d为实测机场受洪水灾害的RDC_Kernel变化检测结果示意图；

图5e为实测机场受洪水灾害的本发明的变化检测结果示意图；

图6a为实测城市受洪水灾害的第一时刻ERS-2SAR图像；

图6b为实测城市受洪水灾害的第二时刻ERS-2SAR图像；

图6c为实测城市受洪水灾害的变化检测参考图；

图6d为实测城市受洪水灾害的RDC_Kernel变化检测结果示意图；

图6e为实测城市受洪水灾害的本发明的变化检测结果示意图；

图7a为实测农田的第一时刻SAR图像；

图7b为实测农田的第二时刻SAR图像；

图7c为实测农田的变化检测参考图；

图7d为实测农田的RDC_Kernel变化检测结果示意图；

图7e实测农田的本发明的变化检测结果示意图。

具体实施方式

下面结合附图对本发明作进一步说明：

参照图1，为本发明的基于支持向量机和判别随机场的SAR图像变化检测方法的流程示意图。该基于支持向量机和判别随机场的SAR图像变化检测方法包括以下步骤：

S1：利用合成孔径雷达接收原始两时相图像，原始两时相图像包括第1时刻图像和第2时刻图像，第1时刻图像和第2时刻图像为同一场景大小相同、时段不同的两幅SAR图像；然后分别对第1时刻图像和第2时刻图像进行灰度值归一化处理，得到第1时刻归一化图像X₁和第2时刻归一化图像X₂；所述第k时刻归一化图像X_k中第i行第j列的像素点表示为X'_k(i,j)，k取1和2，i取1至I，j取1至J，I为第1时刻归一化图像X₁的长度，J为第1时刻归一化图像X₁的宽度；提取X'_k(i,j)的灰度值g'_k(i,j)和X'_k(i,j)的纹理特征w'_k(i,j)。

X'_k(i,j)的灰度值g'_k(i,j)为：

${g^{\&prime;}}_k(i,j)=\frac{g_k(i,j)-\min \left(g_k\right)}{max\left(g_k\right)-\min \left(g_k\right)}$

其中，g_k(i,j)为所述第k时刻图像中第i行第j列的像素点的灰度值，min(g_k)为所述第k时刻图像中所有像素点的灰度值的最小值，max(g_k)为所述第k时刻图像中所有像素点的灰度值的最大值。

以X'_k(i,j)为中心像素点，建立对应的正方形像素窗口，上述正方形像素窗口的边长为η个像素点，η为大于1的奇数。如果X'_k(i,j)位于对应的归一化图像的边缘处，则需要对归一化图像进行向外扩展。此时，X'_k(i,j)的纹理特征w'_k(i,j)包括：对应的正方形像素窗口内像素点灰度值的均值μ'_k(i,j)、对应的正方形像素窗口内像素点灰度值的方差σ²'_k(i,j)、对应的正方形像素窗口内像素点灰度值的峰态ku'_k(i,j)、对应的正方形像素窗口内像素点灰度值的三阶矩sk'_k(i,j)、对应的正方形像素窗口内像素点灰度值的能量en'_k(i,j)、以及对应的正方形像素窗口内像素点灰度值的熵ent'_k(i,j)。X'_k(i,j)的纹理特征w'_k(i,j)按如下公式进行计算：

${{\mathrm{\&mu;}}^{\&prime;}}_k(i,j)=\frac{{\mathrm{\&Sigma;}}_{i,j}{g^{\&prime;}}_k(i,j)}{{\mathrm{\&eta;}}^2}$

${{\mathrm{\&sigma;}}^{2^{\&prime;}}}_k(i,j)=\frac{{\mathrm{\&Sigma;}}_{i,j}\left({g^{\&prime;}}_k\right(i,j)-{{\mathrm{\&mu;}}^{\&prime;}}_k(i,j\left)\right)}{{\mathrm{\&eta;}}^2-1}$

${{\mathrm{ku}}^{\&prime;}}_k(i,j)=\frac{{\mathrm{\&Sigma;}}_{i,j}{\left({g^{\&prime;}}_k\right(i,j)-{{\mathrm{\&mu;}}^{\&prime;}}_k(i,j\left)\right)}^4}{({\mathrm{\&eta;}}^2-1){{\mathrm{\&sigma;}}^{4^{\&prime;}}}_k(i,j)}$

${{\mathrm{sk}}^{\&prime;}}_k(i,j)=\frac{{\mathrm{\&Sigma;}}_{i,j}{\left({g^{\&prime;}}_k\right(i,j)-{{\mathrm{\&mu;}}^{\&prime;}}_k(i,j\left)\right)}^3}{({\mathrm{\&eta;}}^2-1){{\mathrm{\&sigma;}}^{3^{\&prime;}}}_k(i,j)}$

${{\mathrm{en}}^{\&prime;}}_k(i,j)=\underset{i,j}{\mathrm{\&Sigma;}}{\left({{\mathrm{\&mu;}}^{\&prime;}}_k\right(i,j\left)\right)}^2$

${{\mathrm{ent}}^{\&prime;}}_k(i,j)=-\underset{i,j}{\mathrm{\&Sigma;}}{p}_{ij}ln\left(p_{ij}\right),p_{ij}=\frac{{{\mathrm{\&mu;}}^{\&prime;}}_k(i,j)}{\underset{i,j}{\mathrm{\&Sigma;}}{{\mathrm{\&mu;}}^{\&prime;}}_k(i,j)}$

然后按照以下公式获得灰度特征差值△g(i,j)以及纹理特征差值△w(i,j)：△g(i,j)＝g'₂(i,j)-g'₁(i,j)，△w(i,j)＝w'₂(i,j)-w'₁(i,j)；然后将△g(i,j)和△w(i,j)组合成X'_k(i,j)的差值特征向量y(i,j)：y(i,j)＝{△g(i,j),△w(i,j)}。

S2：对第1时刻归一化图像X₁和第2时刻归一化图像X₂按照灰度值作差值运算，得到差值图像△X，利用加权平均比率(ROEWA)算子提取差值图像中第t个像素点的边界强度r_t，t取1至M，M＝I×J。ROEWA算子是基于线性最小均方误差准则的指数滤波器，其计算结果为经过指数加权处理后的均值。具体说明如下：

上述差值图像的边界强度包括差值图像中每个像素点的边界强度，首先定义平滑函数f(ρ)、因果滤波器f₁(ρ)和非因果滤波器f₂(ρ)：

$f\left(\mathrm{\&rho;}\right)=\frac{1}{1+d}{f}_1\left(\mathrm{\&rho;}\right)+\frac{d}{1+d}{f}_2(\mathrm{\&rho;}-1)$

其中，f₁(ρ)＝cd^ρu(ρ)，f₂(ρ)＝cd^-ρu(-ρ)，d为设定常数且0<d<1，u(·)表示Heaviside函数，ρ为自变量；

然后，将ROEWA算子定义为：

$r_{Imax}(i,j)=max\{\frac{{\mathrm{\&mu;}}_{I1}(i-1,j)}{{\mathrm{\&mu;}}_{I2}(i+1,j)},\frac{{\mathrm{\&mu;}}_{I2}(i+1,j)}{{\mathrm{\&mu;}}_{I1}(i-1,j)}\}$

$r_{Jmax}(i,j)=max\{\frac{{\mathrm{\&mu;}}_{J1}(i,j-1)}{{\mathrm{\&mu;}}_{J2}(i,j+1)},\frac{{\mathrm{\&mu;}}_{J2}(i,j+1)}{{\mathrm{\&mu;}}_{J1}(i,j-1)}\}$

其中，μ_I1，μ_I1，μ_I1，μ_I1为指数加权值，可按下式进行计算：

${\mathrm{\&mu;}}_{I1}(i,j)=f_1\left(i\right)*\left(f\right(j)\&CircleTimes;y(i,j\left)\right)$

${\mathrm{\&mu;}}_{I2}(i,j)=f_2\left(i\right)*\left(f\right(j)\&CircleTimes;y(i,j\left)\right)$

μ_J1(i,j)＝f₁(j)*(f(i)*y(i,j))

μ_J2(i,j)＝f₂(j)*(f(i)*y(i,j))

其中，*表示水平方向上的卷积，表示垂直方向上的卷积。

综上所述，我们可以得到ROEWA算子所定义的边界强度|r_max(i,j)|为：

$\left|r_{max}(i,j)\right|=\sqrt{r_{\mathrm{Im}ax}^2(i,j)+r_{Jmax}^2(i,j)}.$

S3：在差值图像中选取训练样本，将训练样本用对应的差值特征向量进行表示，通过训练支持向量机，得到测试样本的初始分类标签、以及测试样本的分类标签的后验概率；具体说明如下：

首先根据原始两时相图像中N组像素点，选取对应的N个有标签的训练样本；每组像素点包括：在第1时刻图像和第2时刻图像中处于相同位置的两个像素点；每个训练样本的标签的设置过程如下：通过对第1时刻图像和第2时刻图像进行观察对比，将所述N组像素点分为变化类像素点组和非变化类像素点组，根据N组像素点的分类情况来设置对应N个训练样本的标签。

将上述N个有标签的训练样本表示为其中，第s个有标签的训练样本表示为(x_s,l_s)，其中，s取1至N；x_s＝y(s_i,s_j)，s_i为第s个训练样本对应的像素点的横坐标，s_j为第s个训练样本对应的像素点的纵坐标；l_s表示第s个训练样本的分类标签，当第s个训练样本对应的一组像素点为变化类像素点组时，l_s＝1；当第s个训练样本对应的一组像素点为非变化类像素点组时，l_s＝0；

在支持向量机(SVM)中建立如下C-SVC模型：

$\begin{array}{cc}\underset{\mathrm{\&alpha;}}{min}& \frac{1}{2}{\mathrm{\&alpha;}}^{T}Q\mathrm{\&alpha;}-{\mathrm{\&Theta;}}^T\mathrm{\&alpha;}\end{array}$

s.t.l^Tα＝0

0≤α_s≤C,s＝1,...,N

其中，α＝[α₁,...,α_N]^T，α_s为待求的第s个训练样本对应的权重，Q为N×N维半正定矩阵，且Q中第p行第q列的元素Q_pq＝l_pl_qK(x_p,x_q)，p取1至N，q取1至N；K(x_p,x_q)为RBF(radial basis function)核函数，l＝[l₁,...,l_N]^T，Θ为N×1维列向量，Θ中的元素均为1；C和γ的取值由交叉验证确定。

下面举例说明C和γ的取值的确定过程：设定C的取值范围和γ的取值范围：C∈(2^-8，2^-7.5,...,2^7.5，2⁸)，γ∈(2^-8，2^-7.5,...,2^7.5，2⁸)，即：C取2^ψ，γ取ψ从-8开始以0.5的间隔进行取值，共有17个值。从-8开始以0.5的间隔进行取值，共有17个值。然后采用5层交叉验证，即将训练样本分成5组，将每组训练样本分别做一次验证集，其余的4组训练样本作为训练集，这样会得到5个模型，比较这5个模型的最终的验证集的分类准确率。通过比较交叉验证精度，选择最优的(C,γ)组合，对于交叉验证精度相同的(C,γ)组合，选择C值最小的(C,γ)组合；

将差值图像中的每个像素点作为对应的一个测试样本；第t个像素点表示为(x'_t,l'_t)，t取1至M，M为差值图像中像素点的个数，M＝I×J；x'_t＝y(t_i,t_j)，t_i为第t个像素点的横坐标，t_j为第t个像素点的纵坐标；l'_t表示待求的第t个像素点的分类标签。

根据上述C-SVC模型得出决策值

${\mathrm{dec}}_{{x^{\&prime;}}_t}=\underset{s=1}{\overset{N}{\mathrm{\&Sigma;}}}{l}_s{\mathrm{\&alpha;}}_{s}K(x_s,{x^{\&prime;}}_t)$

其中，在得出之后，得出第t个像素点的初始分类标签：当时，第t个像素点的初始分类标签为0；反之，第t个像素点的初始分类标签为1。

然后通过拟合sigmoid函数计算出第t个像素点的分类标签的后验概率p(l'_t|y(t_i,t_j))：

$p\left({l^{\&prime;}}_t\right|y\left(t_i,t_j\right))=\frac{1}{1+e^{{\mathrm{Adec}}_{{x^{\&prime;}}_t}+B}}$

其中，A和B通过以下公式确定：

$\underset{z=(A,B)}{min}-\underset{s=1}{\overset{N}{\mathrm{\&Sigma;}}}\left(t_{s}ln\right(p_s)+(1-t_s\left)ln\right(1-p_s\left)\right)$

$t_s=\left\{\begin{array}{cc}\frac{N_{+}+1}{N_{+}+2}& l_s=0\\ \frac{1}{N_{-}+2}& l_s=1\end{array}\right.,p_s=\frac{1}{1+e^{{\mathrm{Adec}}_{{x^{\&prime;}}_t}+B}}$

其中，N₊为：分类标签为1的训练样本的个数，N_-为：分类标签为0的训练样本的个数。

S4：根据差值图像中每个像素点的边界强度、以及测试样本的后验概率，得出初始的支持向量机—判别随机场模型。具体说明如下：

在给定观测场的条件下，标记场(指每个像素点对应的分类标签)L＝(l'_t)_t∈S满足马尔科夫性，即后验概率分布满足如下性质：

$\begin{array}{cc}\&ForAll;x,& p\left({l^{\&prime;}}_t\right|y\left(t_i,t_j\right),{l^{\&prime;}}_{S-\left\{t\right\}})=p\left({l^{\&prime;}}_t\right|y\left(t_i,t_j\right),{l^{\&prime;}}_{N_t})\end{array}$

$\begin{array}{cc}\&ForAll;x,& p\left(l^{\&prime;}\right|y)>0\end{array}$

其中S-{t}为差值图像中除第t个像素点之外的像素点集，l'_S-{t}为位于像素点集S-{t}上的标记场，N_t为第t个像素点的邻域系统(与第t像素点相邻的像素点的集合)，为位于第t个像素点的邻域系统内的标记场。

依据Hammersley-Clifford理论，在只考虑双基团势能的情况下，DRF模型的后验概率分布可以表示为：

$p\left(l^{\&prime;}\right|y,r)=\frac{1}{Z}\exp \left(\underset{a\&Element;S}{\mathrm{\&Sigma;}}A\right({l^{\&prime;}}_a,y\left(a_i,a_j\right))+\underset{a\&Element;S}{\mathrm{\&Sigma;}}I({l^{\&prime;}}_a,{l^{\&prime;}}_b,r\left)\right)$

其中，Z为常数，A为支持向量机—判别随机场模型中的联合势能函数，I为支持向量机—判别随机场模型中的相互势能函数，S表示差值图像中所有像素点的集合，a∈S的含义为：第a个像素点位于S中；l'_a表示第a个像素点的分类标签，l'_b表示第b个像素点的分类标签，第b个像素点位于第a个像素点的邻域系统中。

根据上述差值图像中每个像素点的边界强度，构建支持向量机—判别随机场(SVM-DRF)模型中的相互势能函数I(l'_a,l'_b,r)：

$\begin{array}{c}I({l^{\&prime;}}_a,{l^{\&prime;}}_b,r)=\underset{(a,b)\&Element;N_H}{\mathrm{\&Sigma;}}\left(2\mathrm{\&delta;}\right({l^{\&prime;}}_a,{l^{\&prime;}}_b)-1)\&times;{\mathrm{\&alpha;}}_H\&times;\exp (-(r_a-r_b)/e\mathrm{deg}e\_C^2)\\ +\underset{(a,b)\&Element;N_V}{\mathrm{\&Sigma;}}\left(2\mathrm{\&delta;}\right({l^{\&prime;}}_a,{l^{\&prime;}}_b)-1)\&times;{\mathrm{\&alpha;}}_V\&times;\exp (-(r_a-r_b)/e\mathrm{deg}e\_C^2)\end{array}$

其中，(a,b)∈N_H的含义为：第a个像素点和第b个像素点水平相邻，N_H表示水平邻域系统；l'_a表示待求的第a个像素点对应的分类标签，l'_b表示待求的第b个像素点对应的分类标签；(a,b)∈N_V的含义为：第a个像素点和第b个像素点垂直相邻，N_V表示水平邻域系统；r_a表示差值图像中第a个像素点的边界强度，r_b表示差值图像中第b个像素点的边界强度；edege_C为设定常数，α_H和α_V为I(l'_a,l'_b,r)的两个参数；当l'_a＝l'_b时，δ(l'_a,l'_b)＝1，反之，δ(l'_a,l'_b)＝0。

构建支持向量机—判别随机场(SVM-DRF)模型中的联合势能函数A(l'_a,y(a_i,a_j))，A(l'_a,y(a_i,a_j))＝p(l'_a|y(a_i,a_j))，其中，p(l'_a|y(a_i,a_j))为第a个像素点的分类标签的后验概率。

从而构建初始的支持向量机—判别随机场模型p(l'|y,r)：

$p\left(l^{\&prime;}\right|y,r)=\frac{1}{Z}\exp \left(\underset{a\&Element;S}{\mathrm{\&Sigma;}}A\right({l^{\&prime;}}_a,y\left(a_i,a_j\right))+\underset{a\&Element;S}{\mathrm{\&Sigma;}}I({l^{\&prime;}}_a,{l^{\&prime;}}_b,r\left)\right)$

将α_H和α_V用参数θ进行表示，即θ＝{α_H,α_V}，然后利用最小二乘法估计出θ的初始值θ₀，利用最小二乘法估计出θ的初始值θ₀包括以下步骤：

首先利用直方图工具估计在每个像素点的邻域系统中，都有一组对应的邻域系统标记场；假设差值图像中存在K₁个完全不同的邻域系统标记场。如果在差值图像中第a个像素点对应的邻域系统标记场出现了次，则按照以下公式计算 其中，N_a表示第a个像素点的邻域系统。

根据如下公式即可求得参数θ：

${\left(\mathrm{\&phi;}\right({l^{\&prime;}}_k,{l^{\&prime;}}_{N_k})-\mathrm{\&phi;}({l^{\&prime;}}_h,{l^{\&prime;}}_{N_h}\left)\right)}^T\mathrm{\&theta;}=\ln \left(\frac{P({l^{\&prime;}}_h,{l^{\&prime;}}_{N_h})}{P({l^{\&prime;}}_k,{l^{\&prime;}}_{N_k})}\right)$

$\mathrm{\&phi;}({l^{\&prime;}}_a,{l^{\&prime;}}_{N_a})={\&lsqb;\left(Q\right({l^{\&prime;}}_a,{l^{\&prime;}}_{u_1})+Q({l^{\&prime;}}_a,{l^{\&prime;}}_{u_3}\left)\right),\left(Q\right({l^{\&prime;}}_a,{l^{\&prime;}}_{u_2})+Q({l^{\&prime;}}_a,{l^{\&prime;}}_{u_4}\left)\right)\&rsqb;}^T$

其中，N_k表示差值图像中第k个像素点的邻域系统，N_h表示差值图像中第h个像素点的邻域系统，第k个像素点和第h个像素点表示差值图像中任意两个不同位置的像素点。u₁和u₃分别为与差值图像中第a个像素点水平相邻的两个像素点，u₂和u₄分别为与差值图像中第a个像素点垂直相邻的两个像素点，d取1至4。

S5：根据所述初始分类标签和初始的支持向量机—判别随机场模型，更新支持向量机—判别随机场模型的相互势能函数，得出对应的测试样本的最终分类标签；根据上述对应的测试样本的最终分类标签，得出SAR图像的变化检测结果。参照图2，为本发明中最终输出分类标签的流程示意图。步骤S5具体包括以下步骤：

S51：利用差值图像的每个像素点的初始分类标签构成原始标记场，设定k＝1。

S52：将第a个像素点的分类标签l'_a设为0，将θ的当前取值以及l'_a代入初始的支持向量机—判别随机场模型p(l'|y,r)中，计算出 表示l'_a为0时的p(l'|y,r)；将第a个像素点的分类标签l'_a设为1，将θ的当前取值以及l'_a代入初始的支持向量机—判别随机场模型p(l'|y,r)中，计算出表示l'_a为1时的p(l'|y,r)；将对应的分类标签作为第a个像素点更新后的分类标签。

S53：采用ICE迭代算法对参数θ和标记场进行更新，上述标记场指差值图像中每个像素点的分类标签。参照图3，为ICE迭代算法的流程示意图。采用ICE迭代算法对参数θ和标记场进行更新包括以下步骤：

S531：根据当前的标记场和θ的当前取值，根据Gibbs采样定理，求出新的标记场，即分类标签。具体地说，将θ的当前取值以及l'_a代入初始的支持向量机—判别随机场模型p(l'|y,r)中，计算出表示l'_a为0时的p(l'|y,r)；将第a个像素点的分类标签l'_a设为1，将θ的当前取值以及l'_a代入初始的支持向量机—判别随机场模型p(l'|y,r)中，计算出表示l'_a为1时的p(l'|y,r)；将对应的分类标签作为第a个像素点更新后的分类标签。设定ICE迭代参数τ＝1。

S532：根据第a个像素点当前的分类标签，利用最小二乘法对θ进行重新估计，得出θ的重新估计值θ^(τ)；得出θ的重新估计值θ^(τ)的过程与利用最小二乘法估计出θ的初始值θ₀的过程类似，在此不再重复。

S533：令τ＝τ+1，判断τ是否小于T，T为设定值且T为大于1的自然数，如果τ小于T，返回执行步骤S531；τ＝τ+1的含义为：利用τ+1对τ赋值。如果τ＝T，则得出

S54：令k＝k+1，判断k是否小于K，K为设定值且K为大于1的自然数，如果k小于K，返回执行步骤S52；k＝k+1的含义为：利用k+1对k赋值。如果k＝K，则将当前的标记场作为最终的标记场，然后按照最终的标记场得出SAR图像的变化检测结果。

本发明效果可以通过以下实验进一步证实：

实验内容：

分别利用比值差值合成核(RDC_Kernel)变化检测方法和本发明对SAR图像进行变化检测，以验证本发明有很强的利用上下文信息的能力和多特征融合能力，并验证本发明能够有效地提高变化检测的准确率。

为了验证本发明相对RDC_Kernel变化检测方法在SAR图像变化检测的优势，本发明选取变化检测精度和Kappa系数作为性能指标参数，评价检测结果。该实验结果如图4、图5、图6、图7所示，参照图4a，为实测稻田受洪水灾害的第一时刻ERS-1SAR图像；参照图4b，为实测稻田受洪水灾害的第二时刻ERS-1SAR图像；参照图4c，为实测稻田受洪水灾害的变化检测参考图；参照图4d，为实测稻田受洪水灾害的RDC_Kernel变化检测结果示意图；参照图4e，为实测稻田受洪水灾害的本发明的变化检测结果示意图。参照图5a，为实测机场受洪水灾害的第一时刻JERS SAR channel1图像；参照图5b，为实测机场受洪水灾害的第二时刻JERS SAR channel1图像；参照图5c，为实测机场受洪水灾害的变化检测参考图；参照图5d，为实测机场受洪水灾害的RDC_Kernel变化检测结果示意图；参照图5e，为实测机场受洪水灾害的本发明的变化检测结果示意图。参照图6a，为实测城市受洪水灾害的第一时刻ERS-2SAR图像；参照图6b，为实测城市受洪水灾害的第二时刻ERS-2SAR图像；参照图6c，为实测城市受洪水灾害的变化检测参考图；参照图6d，为实测城市受洪水灾害的RDC_Kernel变化检测结果示意图；参照图6e，为实测城市受洪水灾害的本发明的变化检测结果示意图。参照图7a，为实测农田的第一时刻SAR图像；参照图7b，为实测农田的第二时刻SAR图像；参照图7c，为实测农田的变化检测参考图；参照图7d，为实测农田的RDC_Kernel变化检测结果示意图；参照图7e，实测农田的本发明的变化检测结果示意图。从图4至图7可以看出，与RDC_Kernel变化检测相比，本发明的基于支持向量机和判别随机场的SAR图像变化检测方法有效地提高了变化检测精度。

表1两种方法对真实SAR图像变化检测结果精度比较

实验结果分析：

从表1的实验结果表明本发明的基于支持向量机和判别随机场的SAR图像变化检测方法，相比于比值差值核的检测方法，在利用上下文信息的能力上、抗噪性能和检测精度上更具优势。因为判别随机场模型，联合势能函数通过SVM分类器构造分类时考虑了图像了纹理信息和强度特征，无需数据降维，在训练速度等方面有较高的性能；相互势能函数通过相邻像素点的梯度信息构造，可充分调节标记信息的相互作用强度，减少了误分，从而有效地提高了分类精度。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (4)

1.基于支持向量机和判别随机场的SAR图像变化检测方法，其特征在于，包括以下步骤：

S1：利用合成孔径雷达接收原始两时相图像，原始两时相图像包括第1时刻图像和第2时刻图像；然后分别对第1时刻图像和第2时刻图像进行灰度值归一化处理，得到第1时刻归一化图像X₁和第2时刻归一化图像X₂；第k时刻归一化图像X_k中第i行第j列的像素点表示为X'_k(i,j)，k取1和2，i取1至I，j取1至J，I为第1时刻归一化图像X₁的长度，J为第1时刻归一化图像X₁的宽度；提取X'_k(i,j)的灰度值g'_k(i,j)和X'_k(i,j)的纹理特征w'_k(i,j)；按照以下公式获得灰度特征差值Δg(i,j)以及纹理特征差值Δw(i,j)：Δg(i,j)＝g'₂(i,j)-g'₁(i,j)，Δw(i,j)＝w'₂(i,j)-w'₁(i,j)；然后将Δg(i,j)和Δw(i,j)组合成X'_k(i,j)的差值特征向量y(i,j)：y(i,j)＝{Δg(i,j),Δw(i,j)}；

S2：对第1时刻归一化图像X₁和第2时刻归一化图像X₂按照灰度值作差值运算，得到差值图像ΔX，利用加权平均比率算子提取差值图像中第t个像素点的边界强度r_t，t取1至M，M＝I×J；

S3：在差值图像中选取训练样本，将训练样本用对应的差值特征向量进行表示，通过训练支持向量机，得到测试样本的初始分类标签、以及测试样本的分类标签的后验概率；

首先根据原始两时相图像中N组像素点，选取对应的N个有标签的训练样本；每组像素点包括：在第1时刻图像和第2时刻图像中处于相同位置的两个像素点；每个训练样本的标签的设置过程如下：通过对第1时刻图像和第2时刻图像进行观察对比，将所述N组像素点分为变化类像素点组和非变化类像素点组，根据N组像素点的分类情况来设置对应N个训练样本的标签；

将所述N个有标签的训练样本表示为其中，第s个有标签的训练样本表示为(x_s,l_s)，其中，s取1至N；x_s＝y(s_i,s_j)，s_i为第s个训练样本 对应的像素点的横坐标，s_j为第s个训练样本对应的像素点的纵坐标；ls表示第s个训练样本的分类标签，当第s个训练样本对应的一组像素点为变化类像素点组时，l_s＝1；当第s个训练样本对应的一组像素点为非变化类像素点组时，l_s＝0；

在支持向量机中建立如下C-SVC模型：

s.t.l^Tα＝0

0≤α_s≤C,s＝1,...,N

其中，α＝[α₁,...,α_N]^T，α_s为待求的第s个训练样本对应的权重，Q为N×N维半正定矩阵，且Q中第p行第q列的元素Q_pq＝l_pl_qK(x_p,x_q)，p取1至N，q取1至N；K(x_p,x_q)为核函数，l＝[l₁,...,l_N]T，Θ为N维列向量，Θ中的元素均为1；C和γ的取值交叉验证确定；

在步骤S3中，将差值图像中的每个像素点作为对应的一个测试样本；第t个像素点表示为(x'_t,l'_t)，t取1至M，M为差值图像中像素点的个数；x'_t＝y(t_i,t_j)，t_i为第t个像素点的横坐标，t_j为第t个像素点的纵坐标；l'_t表示待求的第t个像素点的分类标签；

根据所述C-SVC模型得出决策值

其中，在得出之后，得出第t个像素点的初始分类标签：当时，第t个像素点的初始分类标签为0；反之，第t个像素点的初始分类标签为1；

计算出第t个像素点的分类标签的后验概率p(l'_t|y(t_i,t_j))：

其中，A和B通过以下公式确定：

其中，N₊为：分类标签为1的训练样本的个数，N_-为：分类标签为0的训练样本的个数；

S4：根据所述差值图像中每个像素点的边界强度、以及测试样本的后验概率，得出初始的支持向量机—判别随机场模型；

根据所述差值图像中每个像素点的边界强度，构建支持向量机—判别随机场模型中的相互势能函数I(l'_a,l'_b,r)：

其中，(a,b)∈N_H的含义为：第a个像素点和第b个像素点水平相邻；(a,b)∈N_V的含义为：第a个像素点和第b个像素点垂直相邻；edege_C为设定常数，α_H和α_V为I(l'_a,l'_b,r)的两个参数，r_a表示差值图像中第a个像素点的边界强度，r_b表示差值图像中第b个像素点的边界强度；当l'_a＝l'_b时，δ(l'_a,l'_b)＝1，反之，δ(l'_a,l'_b)＝0；

构建支持向量机—判别随机场模型中的联合势能函数A(l'_a,y(a_i,a_j))，A(l'_a,y(a_i,a_j))＝p(l'_a|y(a_i,a_j))，其中，p(l'_a|y(a_i,a_j))为第a个像素点的分类标签的后验概率；

构建初始的支持向量机—判别随机场模型p(l'|y,r)：

其中，Z为设定常数，S表示差值图像中所有像素点的集合；

将α_H和α_V用参数θ进行表示，即θ＝{α_H,α_V}，然后利用最小二乘法估计出θ的初始值θ₀；S5：根据所述初始分类标签和初始的支持向量机—判别随机场模型，更新支持向量机—判别随机场模型的相互势能函数，得出对应的测试样本的最终分类标签；根据所述对应的测试样本的最终分类标签，得出SAR图像的变化检测结果。

2.如权利要求1所述的基于支持向量机和判别随机场的SAR图像变化检测方法，其特征在于，在步骤S1中，X'_k(i,j)的灰度值g'_k(i,j)为：

其中，g_k(i,j)为所述第k时刻图像中第i行第j列的像素点的灰度值，min(g_k)为所述第k时刻图像中所有像素点的灰度值的最小值，max(g_k)为所述第k时刻图像中所有像素点的灰度值的最大值；

在步骤S1中，以X'_k(i,j)为中心像素点，建立对应的正方形像素窗口，所述正方形像素窗口的边长为η个像素点，η为大于1的奇数；则X'_k(i,j)的纹理特征w'_k(i,j)包括：对应的正方形像素窗口内像素点灰度值的均值μ'_k(i,j)、对应的正方形像素窗口内像素点灰度值的方差σ²'_k(i,j)、对应的正方形像素窗口内像素点灰度值的峰态ku'_k(i,j)、对应的正方形像素窗口内像素点灰度值的三阶矩sk'_k(i,j)、对应的正方形像素窗口内像素点灰度值的能量en'_k(i,j)、以及对应的正方形像素窗口内像素点灰度值的熵ent'_k(i,j)。

3.如权利要求1所述的基于支持向量机和判别随机场的SAR图像变化检测方法，其特征在于，在步骤S2中，所述差值图像的边界强度包括差值图像中每个像素点的边界强度，定义平滑函数f(ρ)、因果滤波器f₁(ρ)和非因果 滤波器f₂(ρ)：

其中，f₁(ρ)＝cd^ρu(ρ)，f₂(ρ)＝cd^{- ρ}u(-ρ)，d为设定常数且0<d<1，u(·)表示Heaviside函数，ρ为自变量；

然后，根据以下公式计算差值图像中第i行第j列的像素点ΔX(i,j)的边界强度|r_max(i,j)|为：

μ_J1(i,j)＝f₁(j)*(f(i)*y(i,j))

μ_J2(i,j)＝f₂(j)*(f(i)*y(i,j))

其中，*表示水平方向上的卷积，表示垂直方向上的卷积。

4.如权利要求1所述的基于支持向量机和判别随机场的SAR图像变化检测方法，其特征在于，步骤S5具体包括以下步骤：

S51：利用差值图像的每个像素点的初始分类标签构成原始标记场，设定k＝1；

S52：将第a个像素点的分类标签l'_a设为0，将θ的当前取值以及l'_a代入初始的支持向量机—判别随机场模型p(l'|y,r)中，计算出将第a个像素点的分类标签l'_a设为1，将θ的当前取值以及l'_a代入初始的支持向量机—判别随机场模型p(l'y,r)中，计算出

将

对应的分类标签作为第a个像素点更新后的分类标签；

S53：采用ICE迭代算法对参数θ和标记场进行更新，所述标记场指差值图像中每个像素点的分类标签；

S54：令k＝k+1，判断k是否小于K，K为设定值且K为大于1的自然数，如果k小于K，返回执行步骤S52；如果k＝K，则将当前的标记场作为最终的标记场，然后按照最终的标记场得出SAR图像的变化检测结果。