基于非平稳分析与条件随机场的SAR图像变化检测方法

摘要

本发明公开了一种基于非平稳分析与条件随机场的SAR图像变化检测方法，包括以下步骤：(1)对原始第一、第二时相SAR图像分别进行归一化处理，再作对数比值运算，得到对数比值图像，并求取对数比值图像的纹理特征矩阵；(2)将对数比值图像进行平稳性分割，得到A、B平稳区域，构建A、B平稳区域的训练样本，分别利用A、B平稳区域的训练样本训练支撑向量机，得到第一、第二分类标签矩阵；(3)根据第一、第二分类标签矩阵求得对数比值图像的初始分类标签矩阵，并求得对数比值图像中每个像素点的一元势能函数和二元势能函数，进而得到初步算法模型，再计算初步算法模型中的二元势能函数参数，进而求得对数比值图像的最终分类标签矩阵。

说明书

技术领域

本发明涉及图像处理技术领域，具体涉及一种基于非平稳分析与条件随机场的SAR图像变化检测方法，可用于对合成孔径雷达(SAR)图像进行变化检测处理。

背景技术

随着合成孔径雷达(syntheticapertureradar，SAR)技术的逐步成熟和SAR图像分辨率的不断提高，SAR图像的使用逐渐为人们所重视。合成孔径雷达是一种高分辨率成像雷达，其在民用和军用中的应用需要SAR图像变化检测技术作为支撑。SAR图像变化检测通过对不同时期的SAR图像进行比较分析，并根据SAR图像之间的差异分析来获取所需要的地物变化信息。SAR图像变化检测技术可以应用于很多方面，例如在土地资源利用检测方面的应用、在重大自然灾害预防与检测的应用、对海洋变化分析的应用；在军事方面的应用等。

SAR图像变化检测算法一般可分为：基于灰度特征的算法，如直接比较法、统计假设法、预测模型法、相干模型法、主分量法和分类比较法；基于空间特征的算法，如统计纹理特征分析法和面向对象分类法。

近期在SAR图像变化检测上的研究比较多的有：基于多尺度分析的变化检测算法，如Kai-KuangMa提出一种基于双树-复小波变换(DT-CWT，Dual-TreeComplexWaveletTransform)的多尺度变化检测方法，它利用DT-CWT对对数比值图进行多尺度分解，但没有考虑到SAR图像的纹理信息，阈值的选取也是一个棘手的问题；变化分量分析算法，将多通道的SAR图像看作向量，使用变化向量来描述SAR图像从时间t₀到时间t₁变化的方向和大小，通过计算SAR图像变化向量的大小，检测SAR图像有没有发生变化，优点是可以利用较多甚至全部的数据来检测变化像素，并提供变化像素的类型信息。近年来新发展起来的是基于核方法的SAR图像变化检测算法，GustavoCamps-Valls在2008年首先提出了将核方法应用于SAR图像变化检测，该方法首先提取SAR图像的强度信息和纹理信息，然后构造强度纹理比值差值合成核(RDC_kernel)实现SAR图像变化检测，该方法可以有效的实现SAR图像变化检测。但是目前的变化检测方法都没有明确地从空域结构的角度上考虑SAR图像的非平稳性。

发明内容

针对上述现有技术的问题，本发明的目的在于提供一种基于非平稳分析与条件随机场的SAR图像变化检测方法，以提高SAR图像变化检测的检测效率和检测精度。

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

一种基于非平稳分析与条件随机场的SAR图像变化检测方法，其特征在于，包括以下步骤：

步骤1，输入待检测的原始第一时相SAR图像X₁和原始第二时相SAR图像X₂；

步骤2，对原始第一时相SAR图像X₁和原始第二时相SAR图像X₂分别进行归一化处理，并对归一化后的第一时相SAR图像X′₁和归一化后的第二时相SAR图像X′₂作对数比值运算，得到对数比值图像X_b；对数比值图像X_b的大小与原始第一、第二时相SAR图像的大小均相同；

步骤3，利用滑动窗口分别提取归一化后的第一时相SAR图像的纹理特征矩阵和归一化后的第二时相SAR图像的纹理特征矩阵；并将归一化后的第二时相SAR图像的纹理特征矩阵减去归一化后的第一时相SAR图像的纹理特征矩阵，得到对数比值图像的纹理特征矩阵；所述纹理特征矩阵依次由均值特征矩阵、对比度特征矩阵和半方差特征矩阵组合而成；

步骤4，利用三重马尔科夫场模型算法，将对数比值图像X_b进行平稳性分割，得到两种平稳性区域：A平稳区域与B平稳区域；A平稳区域中像素点的个数与B平稳区域中像素点的个数的和等于对数比值图像中像素点的个数；A平稳区域与B平稳区域的平稳性不同，体现为A平稳区域和B平稳区域中像素点的纹理特征不同。所述纹理特征依次包含均值特征、对比度特征和半方差特征；

步骤5，分别选出A平稳区域和B平稳区域的样本点；

步骤6，提取A平稳区域的所有样本点的纹理特征，组成A平稳区域的纹理特征向量，给定A平稳区域的训练标签l_A，结合A平稳区域的训练标签l_A和A平稳区域的纹理特征向量构成A平稳区域的训练样本，利用该训练样本训练支撑向量机，得到第一分类标签矩阵Ac；

步骤7，提取B平稳区域的所有样本点的纹理特征，组成B平稳区域的纹理特征向量，给定B平稳区域的训练标签l_B，结合B平稳区域的训练标签l_B和B平稳区域的纹理特征向量构成B平稳区域的训练样本，利用该训练样本训练支撑向量机，得到第二分类标签矩阵Bc；

步骤8，根据第一分类标签矩阵Ac和第二分类标签矩阵Bc，求得对数比值图像X_b的初始分类标签矩阵Oc，并求得对数比值图像X_b中每个像素点的类条件概率，再对对数比值图像X_b中每个像素点的类条件概率取对数，得到对数比值图像X_b中每个像素点的一元势能函数；

步骤9，利用指数加权均值比率算子提取对数比值图像X_b中每个像素点的边界强度，并构建对数比值图像X_b中每个像素点的二元势能函数；

步骤10，根据对数比值图像X_b中每个像素点的一元势能函数和二元势能函数，得到初步算法模型p(X|Y)，其中，X为标记场，Y为观测场；

所述初步算法模型p(X|Y)的表达式为：

$>p\left(X\right|Y)=\frac{1}{Z}\exp \left(\underset{n=1}{\overset{N}{\Sigma }}A\right(x_n,F_n,U_n,X)+\underset{n=1}{\overset{N}{\Sigma }}\underset{t\in Q_n}{\Sigma }I\left(x_n,x_t,X,Y,U\right))$>

其中，X为标记场，Y为观测场，U为非平稳场，1/Z为初步算法模型p(X|Y)的分布函数，Q_n为对数比值图像X_b中第n个像素点的邻域系统，x_n为标记场X中对应于对数比值图像X_b的第n个像素点的标记值，x_t为标记场X中对应于对数比值图像X_b的像素点t的标记值，n∈{1，2，...，N}，N为对数比值图像X_b的像素点总数，A(x_n，F_n，U_b，X)为对数比值图像X_b中第n个像素点的一元势能函数，I(x_n，x_t，X，Y，U)为对数比值图像X_b中第n个像素点的二元势能函数：

$>\left(\begin{array}{c}I\left(x_n,x_t,X,Y,U\right)=\underset{\left(n,t\right)\in Q_H}{\Sigma }{\alpha }_H^1\left(1-2\delta \left(x_n,x_t\right)\right)-\left({\alpha }_{aH}^2\delta \left(U_n,U_t,g_1\right)+{\alpha }_{bH}^2\delta \left(U_n,U_t,g_2\right)\right)\\ \times \left(1-\delta \left(x_n,x_t\right)\right)\times \exp \left(-\left(r_n-r_t\right)/edge\_C^2\right)\\ +\underset{\left(n,t\right)\in Q_V}{\Sigma }{\alpha }_V^1\left(1-2\delta \left(x_n,x_t\right)\right)-\left({\alpha }_{aV}^2\delta \left(U_n,U_t,g_1\right)+{\alpha }_{bV}^2\delta \left(U_n,U_t,g_2\right)\right)\\ \times \left(1-\delta \left(x_n,x_t\right)\right)\times \exp \left(-\left(r_n-r_t\right)/edge\_C^2\right)\end{array}\right)$>

Q_H为水平邻域系统、Q_V为垂直邻域系统，r_n为对数比值图像X_b中第n个像素点的边界强度，r_t为对数比值图像X_b中像素点t的边界强度，edge_C为常量，为二元势能函数参数，Y为观测场，U为非平稳场，g₁，g₂为非平稳场U的两个不同取值，一般分别取0和1，U_n为对数比值图像X_b中第n个像素点的平稳性标记值，U_t为对数比值图像X_b中像素点t的平稳性标记值，δ为阶跃函数，其中：

步骤11，利用条件迭代期望算法计算初步算法模型p(X|Y)中的二元势能函数参数，根据该二元势能函数参数求得对数比值图像X_b的最终分类标签矩阵，即得到原始第二时相SAR图像相对于原始第一时相SAR图像的变化检测结果。

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

第一，本发明对SAR图像进行变化检测时无需将SAR图像数据进行降维处理，所以本发明方法在收敛度、训练速度等方面有较高的性能。

第二，本发明建立初步算法模型时引入了SAR图像的非平稳信息，克服了条件随机场不能明确地从空域结构的角度上考虑SAR图像的非平稳性的缺点，因而本发明对SAR图像的模型拟合性得到了改善。

第三，本发明将通过三重马尔科夫场模型算法引入SAR图像的非平稳信息很好地与SAR图像的纹理特征进行结合，能够充分考虑SAR图像的空间信息，所以本发明在SAR图像的变化检测中具有很好的抗噪性能和检测性能。

附图说明

下面结合附图说明和具体实施方式对本发明作进一步详细说明。

图1为本发明的流程图。

图2是本发明方法与基于条件随机场的SAR图像变化检测方法对印度尼西亚稻田SAR图像的变化检测结果对比图，其中：

图2a是实测的印度尼西亚稻田受洪水灾害前的第一时相SAR图像；

图2b是实测的印度尼西亚稻田受洪水灾害后的第二时相SAR图像；

图2c是印度尼西亚稻田SAR图像的变化检测结果参考图；

图2d是基于条件随机场的SAR图像变化检测方法对印度尼西亚稻田SAR图像的变化检测结果图；

图2e是本发明方法对印度尼西亚稻田SAR图像的变化检测结果图。

图3是本发明方法与基于条件随机场的SAR图像变化检测方法对澳大利亚哥尼达机场SAR图像的变化检测结果对比图，其中：

图3a是实测的澳大利亚哥尼达机场受洪水灾害前的第一时相SAR图像；

图3b是实测的澳大利亚哥尼达机场受洪水灾害后的第二时相SAR图像；

图3c是澳大利亚哥尼达机场SAR图像的变化检测结果参考图；

图3d是基于条件随机场的SAR图像变化检测方法对澳大利亚哥尼达机场SAR图像的变化检测结果图；

图3e是本发明方法对澳大利亚哥尼达机场SAR图像的变化检测结果图。

图4是本发明方法与基于条件随机场的SAR图像变化检测方法对日本城市SAR图像的变化检测结果对比图，其中：

图4a是实测的日本城市的第一时相SAR图像；

图4b是实测的日本城市的第二时相SAR图像；

图4c是日本城市SAR图像的变化检测结果参考图；

图4d是基于条件随机场的SAR图像变化检测方法对日本城市SAR图像的变化检测结果图；

图4e是本发明方法对日本城市SAR图像的变化检测结果图。

图5是本发明方法与基于条件随机场的SAR图像变化检测方法对加拿大渥太华地区SAR图像的变化检测结果对比图，其中：

图5a是实测的加拿大渥太华地区的第一时相SAR图像；

图5b是实测的加拿大渥太华地区的第二时相SAR图像；

图5c是加拿大渥太华地区SAR图像的变化检测结果参考图；

图5d是基于条件随机场的变化检测方法对加拿大渥太华地区SAR图像的变化检测结果图；

图5e是本发明方法对加拿大渥太华地区SAR图像的变化检测结果图。

具体实施方式

参照图1，本发明的一种基于非平稳分析与条件随机场的SAR图像变化检测方法，包括以下具体步骤：

步骤1，输入待检测的原始第一时相SAR图像X₁和原始第二时相SAR图像X₂。

输入的待检测的原始第一时相SAR图像X₁和原始第二时相SAR图像X₂如下：

$>X_1=\left\{x(i_1^n,j_1^n)\right|i_1^n\in \{1,2,...,I\},j_1^n\in \{1,2,...,J\}\}$>

$>X_2=\left\{x(i_2^n,j_2^n)\right|i_2^n\in \{1,2,...,I\},j_2^n\in \{1,2,...,J\}\}$>

其中，和分别为原始第一时相SAR图像的第n个像素点的坐标和灰度值，和分别为原始第二时相SAR图像的第n个像素点的坐标和灰度值；I为原始第一时相SAR图像的长度，也是原始第二时相SAR图像的长度；J为原始第一时相SAR图像的宽度，也是原始第二时相SAR图像的宽度；n∈{1，2，...，N}，N为原始第一时相SAR图像的像素点总数，也是原始第二时相SAR图像的像素点总数。

步骤2，对原始第一时相SAR图像X₁和原始第二时相SAR图像X₂分别进行归一化处理，并对归一化后的第一时相SAR图像X′₁和归一化后的第二时相SAR图像X′₂作对数比值运算，得到对数比值图像X_b；对数比值图像X_b的大小与原始第一、第二时相SAR图像的大小均相同。

步骤2的具体子步骤为：

2a)采用以下公式分别对原始第一时相SAR图像X₁和原始第二时相SAR图像X₂进行归一化处理，得到归一化后的第一时相SAR图像X′₁和归一化后的第二时相SAR图像X′₂：

$>X_1^{\prime }=\frac{X_1-\min \left(\min \left(X_1\right)\right)}{\max \left(\max \left(X_1\right)\right)-\min \left(\min \left(X_1\right)\right)}$>

$>X_2^{\prime }=\frac{X_2-min\left(min\right(X_2\left)\right)}{max\left(max\right(X_2\left)\right)-min\left(\min \right(X_2\left)\right)}$>

其中，min表示取最小值，max表示取最大值；

2b)采用以下公式对归一化后的第一时相SAR图像X′₁和归一化后的第二时相SAR图像X′₂作对数比值运算，得到对数比值图像X_b：

$>X_b=log\left(\frac{X_1^{\prime }}{X_2^{\prime }}\right)$>

其中，log表示取对数。

步骤3，利用滑动窗口分别提取归一化后的第一时相SAR图像的纹理特征矩阵和归一化后的第二时相SAR图像的纹理特征矩阵；并将归一化后的第二时相SAR图像的纹理特征矩阵减去归一化后的第一时相SAR图像的纹理特征矩阵，得到对数比值图像的纹理特征矩阵；所述纹理特征矩阵依次由均值特征矩阵、对比度特征矩阵和半方差特征矩阵组合而成。

步骤3的具体子步骤为：

3a)采用以下公式分别计算归一化后的第一时相SAR图像的第n个像素点的均值特征和归一化后的第二时相SAR图像的第n个像素点的均值特征

$>\mu \left(i_1^n,j_1^n\right)=\frac{\underset{i_1^n,j_1^n}{\Sigma }\hat{x}\left(i_1^n,j_1^n\right)}{n_0},\mu \left(i_2^n,j_2^n\right)=\frac{\underset{i_2^n,j_2^n}{\Sigma }\hat{x}\left(i_2^n,j_2^n\right)}{n_0}$>

其中，和分别为归一化后的第一时相SAR图像的第n个像素点的坐标和灰度值，和为归一化后的第二时相SAR图像的第n个像素点的坐标和灰度值；$>i_1^n\in \{1,2,...,I\},j_1^n\in \{1,2,...,J\},i_2^n\in \{1,2,...,I\},j_2^n\in \{1,2,...,J\},$>I为归一化后的第一时相SAR图像的长度，也是归一化后的第二时相SAR图像的长度，J为归一化后的第一时相SAR图像的宽度，也是归一化后的第二时相SAR图像的宽度；n∈{1，2，...，N}，N为归一化后的第一时相SAR图像的像素点总数，也是归一化后的第二时相SAR图像的像素点总数；n₀为滑动窗口包含的像素点数；本发明实例中，滑动窗口的尺寸设定为5×5，则n为25；

归一化后的第一时相SAR图像的均值特征矩阵μ₁和归一化后的第二时相SAR图像的均值特征矩阵μ₂分别为：

$>{\mu }_1=\left\{\mu (i_1^n,j_1^n)\right|i_1^n\in \{1,2,...,I\},j_1^n\in \{1,2,...,J\}\}$>

$>{\mu }_2=\left\{\mu (i_2^n,j_2^n)\right|i_2^n\in \{1,2,...,I\},j_2^n\in \{1,2,...,J\}\};$>

3b)采用以下公式分别计算归一化后的第一时相SAR图像的第n个像素点的对比度特征和归一化后的第二时相SAR图像的第n个像素点的对比度特征

$>cs\left(i_1^n,j_1^n\right)=\underset{i_1^n}{\Sigma }\underset{j_1^n}{\Sigma }{\left(i_1^n-j_1^n\right)}^2\frac{\hat{x}\left(i_1^n,j_1^n\right)}{\underset{i_1^n,j_1^n}{\Sigma }\hat{x}\left(i_1^n,j_1^n\right)},cs\left(i_2^n,j_2^n\right)=\underset{i_2^n}{\Sigma }\underset{j_2^n}{\Sigma }{\left(i_2^n-j_2^n\right)}^2\frac{\hat{x}\left(i_2^n,j_2^n\right)}{\underset{i_2^n,j_2^n}{\Sigma }\hat{x}\left(i_2^n,j_2^n\right)};$>

归一化后的第一时相SAR图像的对比度特征矩阵cs₁和归一化后的第二时相SAR图像的对比度特征矩阵cs₂分别为：

$>{\mathrm{cs}}_1=\left\{cs(i_1^n,j_1^n)\right|i_1^n\in \{1,2,...,I\},j_1^n\in \{1,2,...,J\}\}$>

$>{\mathrm{cs}}_2=\left\{cs(i_2^n,j_2^n)\right|i_2^n\in \{1,2,...,I\},j_2^n\in \{1,2,...,J\}\};$>

3c)采用以下公式分别计算归一化后的第一时相SAR图像的第n个像素点的半方差特征和归一化后的第二时相SAR图像的第n个像素点的半方差特征

$>\chi \left(i_1^n,j_1^n\right)=\frac{1}{2N}\underset{i_1^n,j_1^n}{\Sigma }|\hat{x}\left(i_1^n,j_1^n\right)-{\hat{x}}^{\prime }\left(i_1^n,j_1^n\right)|,\chi \left(i_2^n,j_2^n\right)=\frac{1}{2N}\underset{i_2^n,j_2^n}{\Sigma }|\hat{x}\left(i_2^n,j_2^n\right)-{\hat{x}}^{\prime }\left(i_2^n,j_2^n\right)|$>

其中，为归一化后的第一时相SAR图像在像素点的灰度值，为归一化后的第二时相SAR图像在像素点处的灰度值；归一化后的第一时相SAR图像的像素点与归一化后的第一时相SAR图像的第n个像素点的欧几里得距离为h，归一化后的第二时相SAR图像的像素点与归一化后的第二时相SAR图像的第n个像素点的欧几里得距离为h；

归一化后的第一时相SAR图像的半方差特征矩阵χ₁和归一化后的第二时相SAR图像的半方差特征矩阵χ₂分别为：

$>{\chi }_1=\left\{\chi (i_1^n,j_1^n)\right|i_1^n\in \{1,2,...,I\},j_1^n\in \{1,2,...,J\}\}$>

$>{\chi }_2=\left\{\chi (i_2^n,j_2^n)\right|i_2^n\in \{1,2,...,I\},j_2^n\in \{1,2,...,J\}\};$>

3d)将归一化后的第二时相SAR图像的均值特征矩阵μ₂减去归一化后的第一时相SAR图像的均值特征矩阵μ₁，得到对数比值图像的均值特征矩阵μ_b：

$>{\mu }_b={\mu }_2-{\mu }_1=\{{\mu }_b(i_b^n,j_b^n)=\mu (i_2^n,j_2^n)-\mu (i_1^n,j_1^n)|i_b^n\in \{1,2,...,I\},j_b^n\in \{1,2,...,J\}\}$>

其中，和分别为对数比值图像的第n个像素点的坐标和均值特征；3e)将归一化后的第二时相SAR图像的对比度特征矩阵cs₂减去归一化后的第一时相SAR图像的对比度特征矩阵cs₁，得到对数比值图像的对比度特征矩阵cs_b：

$>{\mathrm{cs}}_b={\mathrm{cs}}_2-{\mathrm{cs}}_1=\{{\mathrm{cs}}_b(i_b^n,j_b^n)=cs(i_2^n,j_2^n)-cs(i_1^n,j_1^n)|i_b^n\in \{1,2,...,I\},j_b^n\in \{1,2,...,J\}\};$>

其中，为对数比值图像的第n个像素点的对比度特征；

3f)将归一化后的第二时相SAR图像的半方差特征矩阵χ₂减去归一化后的第一时相SAR图像的半方差特征矩阵χ₁，得到对数比值图像的半方差特征矩阵χ_b：

$>{\chi }_b={\chi }_2-{\chi }_1=\{{\chi }_b(i_b^n,j_b^n)=\chi (i_2^n,j_2^n)-\chi (i_1^n,j_1^n)|i_b^n\in \{1,2,...,I\},j_b^n\in \{1,2,...,J\}\};$>

其中，为对数比值图像的第n个像素点的半方差特征；

3g)在滑动窗口将对数比值图像的均值特征矩阵、对比度特征矩阵和半方差特征矩阵，构成对数比值图像的纹理特征矩阵g_s，g_s＝[μ_b，ca_b，χ_b]。

步骤4，利用三重马尔科夫场模型算法，将对数比值图像X_b进行平稳性分割，得到两种平稳性区域：A平稳区域与B平稳区域；A平稳区域中像素点的个数与B平稳区域中像素点的个数的和等于对数比值图像中像素点的个数；A平稳区域与B平稳区域的平稳性不同，体现为A平稳区域和B平稳区域中像素点的纹理特征不同。所述纹理特征依次包含均值特征、对比度特征和半方差特征。

需要说明的是，纹理特征是SAR图像的众多特征中具有较高区分度的特征，它反映了SAR图像中目标的表面结构组织排列、空间结构关系和相互邻域关系。在SAR图像的不同平稳区域的成像场景中，成像目标的表面结构组织排列、空间结构关系和相互邻域关系是有区别的，所以我们可根据SAR图像不同平稳区域的纹理特征间的差别，对SAR图像进行平稳性分割。SAR图像具有非平稳性，可表现出同一平稳区域的纹理特征相同、不同平稳区域的纹理特征不同的特点。本发明实例中，A平稳区域与B平稳区域的区别在于：A平稳区域中像素点的纹理特征和B平稳区域中像素点的纹理特征不同。本发明实例中，采用三重马尔科夫场模型算法(TripletMarkovField，TMF)将对数比值图像进行平稳性分割，具体地，采用TMF引入非平稳场U，用U的不同取值来描述两种不同的平稳性分布，结合马尔科夫场对对数比值图像进行平稳性分割；具体实现时，也可以采用其他方法对对数比值图像进行平稳性分割，本发明实例对此不做限制。

步骤5，分别选出A平稳区域和B平稳区域的样本点。

步骤5的具体子步骤为：

5a)在对数比值图像中选取N₀个原始第二时相SAR图像相对于原始第一时相SAR图像在视觉上未变化的像素点，再选取N₁个原始第二时相SAR图像相对于原始第一时相SAR图像在视觉上变化的像素点，N₀+N₁＜N，N为原始第一时相SAR图像的像素点总数，也是原始第二时相SAR图像的像素点总数；

5b)根据所选出的每个在视觉上未变化的像素点的坐标，判断其位于A平稳区域还是B平稳区域，将所有位于A平稳区域的在视觉上未变化的像素点构成第一集合K_A0，该第一集合K_A0的大小为N_A0，将所有位于B平稳区域的在视觉上未变化的像素点构成第二集合K_B0，该第二集合K_B0的大小为N_B0；

根据所选出的每个在视觉上变化的像素点的坐标，判断其位于A平稳区域还是B平稳区域，将所有位于A平稳区域的在视觉上变化的像素点构成第三集合K_A1，该第三集合K_A1的大小为N_A1，将所有位于B平稳区域的在视觉上未变化的像素点构成第四集合K_B1，该第四集合K_B1的大小为N_B1；

将第一集合K_A0和第三集合K_A1中的像素点作为A平稳区域的样本点；将第二集合K_B0和第四集合K_B1中的像素点作为B平稳区域的样本点；

其中，N₀＝N_A0+N_B0，N₁＝N_A1+N_B1。

步骤6，提取A平稳区域的所有样本点的纹理特征，组成A平稳区域的纹理特征向量，给定A平稳区域的训练标签l_A，结合A平稳区域的训练标签l_A和A平稳区域的纹理特征向量构成A平稳区域的训练样本，利用该训练样本训练支撑向量机，得到第一分类标签矩阵Ac。

步骤6的具体子步骤为：

6a)根据对数比值图像的纹理特征矩阵g_s，提取第一集合K_A0和第三集合K_A1的像素点的纹理特征，组成A平稳区域的纹理特征向量g_A：

其中，分别为A平稳区域中在视觉上未变化的N_A0个像素点中第1个像素点的均值特征、对比度特征和半方差特征；分别为A平稳区域中在视觉上未变化的N_A0个像素点中第N_A0个像素点的均值特征、对比度特征和半方差特征；分别为A平稳区域中在视觉上变化的N_A1个像素点中第1个像素点的均值特征、对比度特征和半方差特征；分别为A平稳区域中在视觉上变化的N_A1个像素点中第NA₁个像素点的均值特征、对比度特征和半方差特征；

再结合A平稳区域的训练标签l_A构成A平稳区域的训练样本S_pA：

S_pA＝{g_A，l_A}

其中，A平稳区域的训练标签l_A是一个长度为N的行向量，其前N₀个元素为0，其后N₁个元素为1：

6b)将支撑向量机的类型设置为C-SVC模型，并将核函数的类型设置为径向基核函数；

需要说明的是，支撑向量机是一种二类分类模型，其基本定义为特征空间上的间隔最大的线性分类器，其学习策略是间隔最大化；本发明实例中，设置支撑向量机参数s＝0，即将支撑向量机的类型设置为C-SVC模型，设置支撑向量机参数t＝2，即将核函数的类型设置为径向基核函数；

6c)将惩罚因子C和径向基核函数中核参数γ进行交叉验证，得到对应于A平稳区域的最优(C，γ)组合；

具体地，选取惩罚因子C∈{2^-8，2^-7.5，...，2^7.5，2⁸}，径向基核函数中核参数γ∈{2^-8，2^-7.5，...，2^7.5，2⁸}，则共有33×33个(C，γ)组合，对该33×33个(C，γ)组合进行5层交叉验证，具体步骤为：

将A平稳区域的训练样本Sp_A随机地平均分成5组数据：组A1、组A2、组A3、组A4、组A5，将5组数据分别作为一次验证集，并将5组数据分别对应的其余4组数据进行顺序组合作为训练集；首先，将组A1作为验证集，将组A2、组A3、组A4、组A5进行顺序组合作为训练集，将该训练集和每一个(C，γ)组合进行分类，得到33×33个分类结果，将该33×33个分类结果分别与组A1进行比较，求得对应于组A1的每一个(C，γ)组合的分类准确率，选取其中分类准确率最高的(C，γ)组合作为对应于组A1的最优(C，γ)组合；然后，分别将组A2、组A3、组A4、组A5作为验证集，选出对应于组A2、组A3、组A4、组A5的最优(C，γ)组合；最后，比较对应于组A1、组A2、组A3、组A4、组A5的最优(C，γ)组合，将其中分类准确率最高的(C，γ)组合作为对应于A平稳区域的最优(C，γ)组合；在本发明实例中，若对应于组A1、组A2、组A3、组A4、组A5的最优(C，γ)组合中，分类准确率最高的(C，γ)组合有多个，则选取其中惩罚因子C为最小的(C，γ)组合作为对应于A平稳区域的最优(C，γ)组合；

6d)利用对应于A平稳区域的最优(C，γ)组合和A平稳区域的训练样本Sp_A对支撑向量机进行训练，得到第一分类标签矩阵Ac，其维数与对数比值图像X_b的维数相同。

步骤7，提取B平稳区域的所有样本点的纹理特征，组成B平稳区域的纹理特征向量，给定B平稳区域的训练标签l_B，结合B平稳区域的训练标签l_B和B平稳区域的纹理特征向量构成B平稳区域的训练样本，利用该训练样本训练支撑向量机，得到第二分类标签矩阵Bc。

步骤7的具体子步骤为：

7a)根据对数比值图像的纹理特征矩阵g_s，提取第二集合K_B0和第四集合K_B1的像素点的纹理特征，组成B平稳区域的纹理特征向量g_B：

$>g_B=[{\mu }_b({\tilde{i}}_b^1,{\tilde{j}}_b^1),...,{\mu }_b({\tilde{i}}_b^{N_{B0}},{\tilde{j}}_b^{N_{B0}}),{\mu }_b({\hat{i}}_b^1,{\hat{j}}_b^1),...,{\mu }_b({\hat{i}}_b^{N_{B1}},{\hat{j}}_b^{N_{B1}}),$>

$>{\mathrm{cs}}_b({\tilde{i}}_b^1,{\tilde{j}}_b^1),...,{\mathrm{cs}}_b({\tilde{i}}_b^{N_{B0}},{\tilde{j}}_b^{N_{B0}}),{\mathrm{cs}}_b({\hat{i}}_b^1,{\hat{j}}_b^1),...,{\mathrm{cs}}_b({\hat{i}}_b^{N_{B1}},{\hat{j}}_b^{N_{B1}}),$>

$>{\chi }_b({\tilde{i}}_b^1,{\tilde{j}}_b^1),...,{\chi }_b({\tilde{i}}_b^{N_{B0}},{\tilde{j}}_b^{N_{B0}}),{\chi }_b({\hat{i}}_b^1,{\hat{j}}_b^1),...,{\chi }_b({\hat{i}}_b^{N_{B1}},{\hat{j}}_b^{N_{B1}})]$>

其中，分别为B平稳区域中在视觉上未变化的N_B0个像素点中第1个像素点的均值特征、对比度特征和半方差特征；分别为B平稳区域中在视觉上未变化的N_B0个像素点中第N_B0个像素点的均值特征、对比度特征和半方差特征；分别为B平稳区域中在视觉上变化的N_B1个像素点中第1个像素点的均值特征、对比度特征和半方差特征；分别为B平稳区域中在视觉上变化的N_B1个像素点中第N_B1个像素点的均值特征、对比度特征和半方差特征；

再结合B平稳区域的训练标签l_B构成B平稳区域的训练样本Sp_B：

Sp_B＝{g_B，l_B}

其中，B平稳区域的训练标签l_B是一个长度为N的行向量，其前N₀个元素为0，其后N₁个元素为1：

7b)将支撑向量机的类型设置为C-SVC模型，并将核函数的类型设置为径向基核函数；

本发明实例中，设置支撑向量机参数s＝0，即将支撑向量机的类型设置为C-SVC模型，设置支撑向量机参数t＝2，即将核函数的类型设置为径向基核函数；

7c)将惩罚因子C和径向基核函数中核参数γ进行交叉验证，得到对应于A平稳区域的最优的(C，γ)组合；

具体地，选取惩罚因子C∈{2^-8，2^-7.5，...，2^7.5，2⁸}，径向基核函数中核参数γ∈{2^-8，2^-7.5，...，2^7.5，2⁸}，则共有33×33个(C，γ)组合，对该33×33个(C，γ)组合进行5层交叉验证，具体步骤为：

将B平稳区域的训练样本Sp_B随机地平均分成5组数据：组B1、组B2、组B3、组B4、组B5，将5组数据分别作为一次验证集，并将5组数据分别对应的其余4组数据进行顺序组合作为训练集；首先，将组B1作为验证集，将组B2、组B3、组B4、组B5进行顺序组合作为训练集，将该训练集和每一个(C，γ)组合进行分类，得到33×33个分类结果，将该33×33个分类结果分别与组B1进行比较，求得对应于组B1的每一个(C，γ)组合的分类准确率，选取其中分类准确率最高的(C，γ)组合作为对应于组B1的最优(C，γ)组合；然后，分别将组B2、组B3、组B4、组B5作为验证集，选出对应于组B2、组B3、组B4、组B5的最优(C，γ)组合；最后，比较对应于组B1、组B2、组B3、组B4、组B5的最优(C，γ)组合，将其中分类准确率最高的(C，γ)组合作为对应于B平稳区域的最优(C，γ)组合；在本发明实例中，若对应于组B1、组B2、组B3、组B4、组B5的最优(C，γ)组合中，分类准确率最高的(C，γ)组合有多个，则选取其中惩罚因子C为最小的(C，γ)组合作为对应于B平稳区域的最优(C，γ)组合；

7d)利用对应于B平稳区域的最优(C，γ)组合和B平稳区域的训练样本Sp_B对支撑向量机进行训练，得到第二分类标签矩阵Bc，其维数与对数比值图像X_b的维数相同。

步骤8，根据第一分类标签矩阵Ac和第二分类标签矩阵Bc，求得对数比值图像X_b的初始分类标签矩阵Oc，并求得对数比值图像X_b中每个像素点的类条件概率，再对对数比值图像X_b中每个像素点的类条件概率取对数，得到对数比值图像X_b中每个像素点的一元势能函数。

步骤8的具体子步骤为：

8a)构建与原始第一时相SAR图像维数相同的训练标签矩阵Oc′，训练标签矩阵Oc′的初始值为全0矩阵，将第一分类标签矩阵Ac中与A平稳区域中的每个像素点分别对应的元素的值，按A平稳区域中的每个像素点的坐标位置，分别填入训练标签矩阵Oc′；将第二分类标签矩阵Bc中与B平稳区域中的每个像素点分别对应的元素的值，按B平稳区域中的每个像素点的坐标位置，分别填入训练标签矩阵Oc′；将训练标签矩阵Oc′作为对数比值图像X_b的初始分类标签矩阵Oc，对数比值图像X_b的初始分类标签矩阵Oc中元素的取值为0或1，0表示原始第二时相SAR图像中的对应像素点相对于原始第一时相SAR图像中的对应像素点未发生变化，1表示原始第二时相SAR图像中的对应像素点相对于原始第一时相SAR图像中的对应像素点发生变化；

8b)将对数比值图像X_b的初始分类标签矩阵Oc拟合为S形生长曲线函数，即求得对数比值图像X_b中每个像素点的类条件概率，并对对数比值图像X_b中每个像素点的类条件概率取对数，得到对数比值图像X_b中每个像素点的一元势能函数，其中，对数比值图像X_b中第n个像素点的一元势能函数A(x_n，F_n，U_n，X)为：

A(x_n，F_n，U_n，X)＝logp(x_n|d_n(F，U))

其中，p(x_n|d_n(F，U))为对数比值图像X_b中第n个像素点的类条件概率，X为标记场，x_n为标记场X中对应于对数比值图像X_b的第n个像素点的标记值，F_n为对数比值图像X_b中第n个像素点的纹理特征向量，分别为对数比值图像X_b中第n个像素点的均值特征、对比度特征和半方差特征，为对数比值图像X_b中第n个像素点的坐标，U_n为对数比值图像X_b中第n个像素点的平稳性标记值，F为纹理特征场，U为非平稳场，d_n(F，U)为对数比值图像X_b中第n个像素点的纹理特征场与非平稳场的混合场，n∈{1，2，...，N}，N为对数比值图像X_b的像素点总数，log表示取对数。

需要说明的是，S形生长曲线函数即Sigmoid函数，Sigmoid函数是一个在生物学中常见的S型函数，也称为S形生长曲线。

步骤9，利用指数加权均值比率算子提取对数比值图像X_b中每个像素点的边界强度，并构建对数比值图像X_b中每个像素点的二元势能函数。

指数加权均值比率算子(ROEWA)是基于线性最小均方误差准则的指数滤波器，其计算结果为经过指数加权处理后的均值；

步骤9的具体子步骤为：

9a)定义平滑函数f(λ)由因果滤波器f₁(λ)和非因果滤波器f₂(λ)组成：

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

其中，f₁(λ)＝cd^λu(λ)，f₂(λ-1)＝cd^-(λ-1)u(-(λ-1))，d为设定常数，且0＜d＜1，c＝1-d，u(λ)为阶跃函数，λ为自变量；

9b)给出指数加权均值比率算子的定义如下：

$>\left(\begin{array}{c}{r}_{\mathrm{Im}ax}\left(i_b^n,j_b^n\right)=\max \left\{\frac{{\mu }_{I1}\left(i_b^n-1,j_b^n\right)}{{\mu }_{I2}\left(i_b^n+1,j_b^n\right)},\frac{{\mu }_{I2}\left(i_b^n+1,j_b^n\right)}{{\mu }_{I1}\left(i_b^n-1,j_b^n\right)}\right\}\\ r_{Jmax}\left(i_b^n,j_b^n\right)=\max \left\{\frac{{\mu }_{J1}\left(i_b^n,j_b^n-1\right)}{{\mu }_{J2}\left(i_b^n,j_b^n+1\right)},\frac{{\mu }_{J2}\left(i_b^n,j_b^n+1\right)}{{\mu }_{J1}\left(i_b^n,j_b^n-1\right)}\right\}\end{array}\right)$>

其中，表示对数比值图像X_b中第n个像素点的坐标，n∈{1，2，...，N}，N为对数比值图像X_b的像素点总数，μ_I1，μ_I2，μ_J1，μ_J2为指数加权值，按下式进行计算：

$>\left(\begin{array}{c}{\mu }_{I1}(i_b^n,j_b^n)=f_1\left(i_b^n\right)*(f\left(j_b^n\right)\otimes F_n)\\ {\mu }_{I2}(i_b^n,j_b^n)=f_2\left(i_b^n\right)*(f\left(j_b^n\right)\otimes F_n)\\ {\mu }_{J1}(i_b^n,j_b^n)=f_1\left(j_b^n\right)*(f\left(i_b^n\right)*F_n)\\ {\mu }_{J2}(i_b^n,j_b^n)=f_2\left(j_b^n\right)*(f\left(i_b^n\right)*F_n)\end{array}\right)$>

其中，F_n为对数比值图像X_b中第n个像素点的纹理特征向量，$>F_n=[{\mu }_b(i_b^n,j_b^n),{\mathrm{cs}}_b(i_b^n,j_b^n),{\chi }_b(i_b^n,j_b^n)],{\mu }_b(i_b^n,j_b^n),{\mathrm{cs}}_b(i_b^n,j_b^n),{\chi }_b(i_b^n,j_b^n)$>分别为对数比值图像X_b中第n个像素点的均值特征、对比度特征和半方差特征，*表示水平方向上的卷积，则表示垂直方向上的卷积；

9c)通过以下公式提取对数比值图像X_b中第n个像素点的边界强度r_n：

$>r_n=\sqrt{r_{\mathrm{Im}ax}^2(i_b^n,j_b^n)+r_{Jmax}^2(i_b^n,j_b^n)};$>

9c)通过以下公式构建对数比值图像X_b中第n个像素点的二元势能函数I(x_n，x_t，X，Y，U)：

$>\left(\begin{array}{c}I\left(x_n,x_t,X,Y,U\right)=\underset{\left(n,t\right)\in Q_H}{\Sigma }{\alpha }_H^1\left(1-2\delta \left(x_n,x_t\right)\right)-\left({\alpha }_{aH}^2\delta \left(U_n,U_t,g_1\right)+{\alpha }_{bH}^2\delta \left(U_n,U_t,g_2\right)\right)\\ \times \left(1-\delta \left(x_n,x_t\right)\right)\times \exp \left(-\left(r_n-r_t\right)/edge\_C^2\right)\\ +\underset{\left(n,t\right)\in Q_V}{\Sigma }{\alpha }_V^1\left(1-2\delta \left(x_n,x_t\right)\right)-\left({\alpha }_{aV}^2\delta \left(U_n,U_t,g_1\right)+{\alpha }_{bV}^2\delta \left(U_n,U_t,g_2\right)\right)\\ \times \left(1-\delta \left(x_n,x_t\right)\right)\times \exp \left(-\left(r_n-r_t\right)/edge\_C^2\right)\end{array}\right)$>

其中，Q_H为水平邻域系统、Q_V为垂直邻域系统，r_n为对数比值图像X_b中第n个像素点的边界强度，r_t为对数比值图像X_b中像素点t的边界强度，X为标记场，x_n为标记场X中对应于对数比值图像X_b的第n个像素点的标记值，x_t为标记场X中对应于对数比值图像X_b的像素点t的标记值，edge_C为常量，为二元势能函数参数，Y为观测场，U为非平稳场，g₁，g₂为非平稳场U的两个不同取值，一般分别取0和1，U_n为对数比值图像X_b中第n个像素点的平稳性标记值，U₁为对数比值图像X_b中像素点t的平稳性标记值，δ为阶跃函数，其中：

步骤10，根据对数比值图像X_b中每个像素点的一元势能函数和二元势能函数，得到初步算法模型p(X|Y)，其中，X为标记场，Y为观测场。

所述初步算法模型p(X|Y)的表达式为：

$>p\left(X\right|Y)=\frac{1}{Z}\exp \left(\underset{n=1}{\overset{N}{\Sigma }}A\right(x_n,F_n,U_n,X)+\underset{n=1}{\overset{N}{\Sigma }}\underset{t\in Q_n}{\Sigma }I\left(x_n,x_t,X,Y,U\right))$>

其中，X为标记场，Y为观测场，U为非平稳场，1/Z为初步算法模型p(X|Y)的分布函数，Q_n为对数比值图像X_b中第n个像素点的邻域系统，x_n为标记场X中对应于对数比值图像X_b的第n个像素点的标记值，x_t为标记场X中对应于对数比值图像X_b的像素点t的标记值，n∈{1，2，...，N}，N为对数比值图像X_b的像素点总数，A(x_n，F_n，U_n，X)对数比值图像X_b中第n个像素点的一元势能函数，I(x_n，x_t，X，Y，U)为对数比值图像X_b中第n个像素点的二元势能函数：

$>\left(\begin{array}{c}I\left(x_n,x_t,X,Y,U\right)=\underset{\left(n,t\right)\in Q_H}{\Sigma }{\alpha }_H^1\left(1-2\delta \left(x_n,x_t\right)\right)-\left({\alpha }_{aH}^2\delta \left(U_n,U_t,g_1\right)+{\alpha }_{bH}^2\delta \left(U_n,U_t,g_2\right)\right)\\ \times \left(1-\delta \left(x_n,x_t\right)\right)\times \exp \left(-\left(r_n-r_t\right)/edge\_C^2\right)\\ +\underset{\left(n,t\right)\in Q_V}{\Sigma }{\alpha }_V^1\left(1-2\delta \left(x_n,x_t\right)\right)-\left({\alpha }_{aV}^2\delta \left(U_n,U_t,g_1\right)+{\alpha }_{bV}^2\delta \left(U_n,U_t,g_2\right)\right)\\ \times \left(1-\delta \left(x_n,x_t\right)\right)\times \exp \left(-\left(r_n-r_t\right)/edge\_C^2\right)\end{array}\right)$>

Q_H为水平邻域系统、Q_V为垂直邻域系统，r_n为对数比值图像X_b中第n个像素点的边界强度，r_t为对数比值图像X_b中像素点t的边界强度，edge_C为常量，为二元势能函数参数，Y为观测场，U为非平稳场，g₁，g₂为非平稳场U的两个不同取值，一般分别取0和1，U_n为对数比值图像X_b中第n个像素点的平稳性标记值，U_t为对数比值图像X_b中像素点t的平稳性标记值，δ为阶跃函数，其中：

步骤11，利用条件迭代期望算法计算初步算法模型p(X|Y)中的二元势能函数参数，根据该二元势能函数参数求得对数比值图像X_b的最终分类标签矩阵，即得到原始第二时相SAR图像相对于原始第一时相SAR图像的变化检测结果。

步骤11的具体子步骤为：

11a)在对数比值图像X_b的初始分类标签矩阵Oc中随机选取出一个大小为M×M的子矩阵作为初步算法模型p(X|Y)中的标记场X的初始值X(0)，并根据初步算法模型p(X|Y)求得二元势能函数参数的初始值，再利用所求出的二元势能函数参数的初始值，根据初步算法模型p(X|Y)求得第1代标记场X₁作为第1代分类标签矩阵；设置迭代次数t＝1；

11b)在第t代分类标签矩阵中随机取出一个大小为M×M的子矩阵作为初步算法模型p(X|Y)中的第t次迭代的标记场X(t)，并根据初步算法模型p(X|Y)求得第t次迭代的二元势能函数参数，再利用所求出的第t次迭代的二元势能函数参数，根据初步算法模型p(X|Y)求得第t+1代标记场X_t+1作为第t+1代分类标签矩阵；

11c)如果迭代次数t＞9，则将第t+1代分类标签矩阵作为对数比值图像X_b的最终分类标签矩阵，即将第t+1代分类标签矩阵作为原始第二时相SAR图像相对于原始第一时相SAR图像的变化检测结果；反之，令迭代次数t增加1，返回步骤11b)。

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

1)实验条件：

仿真实验的环境为：MATLABR20111b，Libsvm-3.17工具箱，Intel(R)Core(TM)i5-3570CPU3.4GHz，Window7Professional。

2)实验内容：

为了验证本发明方法在SAR图像变化检测中的优势，针对四种SAR图像，比较本发明方法与基于条件随机场的SAR图像变化检测方法的SAR图像变化检测结果，如图2-图5所示，所述四种SAR图像分别为：印度尼西亚稻田SAR图像、澳大利亚哥尼达机场SAR图像、日本城市SAR图像和加拿大渥太华地区SAR图像；

3)实验结果分析：

以SAR图像变化检测精度和卡帕(Kappa)系数作为性能指标，比较本发明方法与基于条件随机场的SAR图像变化检测方法针对四种SAR图像的变化检测结果，如表1所示，表1中，将基于条件随机场的SAR图像变化检测方法简称为基于条件随机场方法：

表1

从表1可以看出，本发明方法相比于基于条件随机场的SAR图像变化检测方法，能够获得更高的SAR图像变化检测精度和卡帕(Kappa)系数，这表明本发明方法在SAR图像变化检测中具有更好的抗噪性能，因为本发明方法考虑了SAR图像的纹理特征和非平稳性，并对SAR图像建立了更为精确的算法模型，更符合SAR图像的实际情况。

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