CN105160666B - 基于非平稳分析与条件随机场的sar图像变化检测方法 - Google Patents

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

Info

Publication number
CN105160666B
CN105160666B CN201510526592.5A CN201510526592A CN105160666B CN 105160666 B CN105160666 B CN 105160666B CN 201510526592 A CN201510526592 A CN 201510526592A CN 105160666 B CN105160666 B CN 105160666B
Authority
CN
China
Prior art keywords
image
pixel point
group
sar image
matrix
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201510526592.5A
Other languages
English (en)
Other versions
CN105160666A (zh
Inventor
李明
朱华慧
张鹏
吴艳
贾璐
安琳
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xidian University
Original Assignee
Xidian University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xidian University filed Critical Xidian University
Priority to CN201510526592.5A priority Critical patent/CN105160666B/zh
Publication of CN105160666A publication Critical patent/CN105160666A/zh
Application granted granted Critical
Publication of CN105160666B publication Critical patent/CN105160666B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00—Scenes; Scene-specific elements
    • G06V20/10—Terrestrial scenes
    • G06V20/13—Satellite images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30108—Industrial image inspection

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Astronomy & Astrophysics (AREA)
  • Remote Sensing (AREA)
  • Multimedia (AREA)
  • Theoretical Computer Science (AREA)
  • Image Analysis (AREA)

Abstract

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

Description

基于非平稳分析与条件随机场的SAR图像变化检测方法
技术领域
本发明涉及图像处理技术领域，具体涉及一种基于非平稳分析与条件随机场的SAR图像变化检测方法，可用于对合成孔径雷达(SAR)图像进行变化检测处理。
背景技术
随着合成孔径雷达(synthetic aperture radar，SAR)技术的逐步成熟和SAR图像分辨率的不断提高，SAR图像的使用逐渐为人们所重视。合成孔径雷达是一种高分辨率成像雷达，其在民用和军用中的应用需要SAR图像变化检测技术作为支撑。SAR图像变化检测通过对不同时期的SAR图像进行比较分析，并根据SAR图像之间的差异分析来获取所需要的地物变化信息。SAR图像变化检测技术可以应用于很多方面，例如在土地资源利用检测方面的应用、在重大自然灾害预防与检测的应用、对海洋变化分析的应用；在军事方面的应用等。
SAR图像变化检测算法一般可分为：基于灰度特征的算法，如直接比较法、统计假设法、预测模型法、相干模型法、主分量法和分类比较法；基于空间特征的算法，如统计纹理特征分析法和面向对象分类法。
近期在SAR图像变化检测上的研究比较多的有：基于多尺度分析的变化检测算法，如Kai-Kuang Ma提出一种基于双树-复小波变换(DT-CWT，Dual-Tree Complex WaveletTransform)的多尺度变化检测方法，它利用DT-CWT对对数比值图进行多尺度分解，但没有考虑到SAR图像的纹理信息，阈值的选取也是一个棘手的问题；变化分量分析算法，将多通道的SAR图像看作向量，使用变化向量来描述SAR图像从时间t0到时间t1变化的方向和大小，通过计算SAR图像变化向量的大小，检测SAR图像有没有发生变化，优点是可以利用较多甚至全部的数据来检测变化像素，并提供变化像素的类型信息。近年来新发展起来的是基于核方法的SAR图像变化检测算法，Gustavo Camps-Valls在2008年首先提出了将核方法应用于SAR图像变化检测，该方法首先提取SAR图像的强度信息和纹理信息，然后构造强度纹理比值差值合成核(RDC_kernel)实现SAR图像变化检测，该方法可以有效的实现SAR图像变化检测。但是目前的变化检测方法都没有明确地从空域结构的角度上考虑SAR图像的非平稳性。
发明内容
针对上述现有技术的问题，本发明的目的在于提供一种基于非平稳分析与条件随机场的SAR图像变化检测方法，以提高SAR图像变化检测的检测效率和检测精度。
为了达到上述目的，本发明采用以下技术方案予以实现。
一种基于非平稳分析与条件随机场的SAR图像变化检测方法，其特征在于，包括以下步骤：
步骤1，输入待检测的原始第一时相SAR图像X1和原始第二时相SAR图像X2；
步骤2，对原始第一时相SAR图像X1和原始第二时相SAR图像X2分别进行归一化处理，并对归一化后的第一时相SAR图像X′1和归一化后的第二时相SAR图像X′2作对数比值运算，得到对数比值图像Xb；对数比值图像Xb的大小与原始第一、第二时相SAR图像的大小均相同；
步骤3，利用滑动窗口分别提取归一化后的第一时相SAR图像的纹理特征矩阵和归一化后的第二时相SAR图像的纹理特征矩阵；并将归一化后的第二时相SAR图像的纹理特征矩阵减去归一化后的第一时相SAR图像的纹理特征矩阵，得到对数比值图像的纹理特征矩阵；所述纹理特征矩阵依次由均值特征矩阵、对比度特征矩阵和半方差特征矩阵组合而成；
步骤4，利用三重马尔科夫场模型算法，将对数比值图像Xb进行平稳性分割，得到两种平稳性区域：A平稳区域与B平稳区域；A平稳区域中像素点的个数与B平稳区域中像素点的个数的和等于对数比值图像中像素点的个数；A平稳区域与B平稳区域的平稳性不同，体现为A平稳区域和B平稳区域中像素点的纹理特征不同。所述纹理特征依次包含均值特征、对比度特征和半方差特征；
步骤5，分别选出A平稳区域和B平稳区域的样本点；
步骤6，提取A平稳区域的所有样本点的纹理特征，组成A平稳区域的纹理特征向量，给定A平稳区域的训练标签lA，结合A平稳区域的训练标签lA和A平稳区域的纹理特征向量构成A平稳区域的训练样本，利用该训练样本训练支撑向量机，得到第一分类标签矩阵Ac；
步骤7，提取B平稳区域的所有样本点的纹理特征，组成B平稳区域的纹理特征向量，给定B平稳区域的训练标签lB，结合B平稳区域的训练标签lB和B平稳区域的纹理特征向量构成B平稳区域的训练样本，利用该训练样本训练支撑向量机，得到第二分类标签矩阵Bc；
步骤8，根据第一分类标签矩阵Ac和第二分类标签矩阵Bc，求得对数比值图像Xb的初始分类标签矩阵Oc，并求得对数比值图像Xb中每个像素点的类条件概率，再对对数比值图像Xb中每个像素点的类条件概率取对数，得到对数比值图像Xb中每个像素点的一元势能函数；
步骤9，利用指数加权均值比率算子提取对数比值图像Xb中每个像素点的边界强度，并构建对数比值图像Xb中每个像素点的二元势能函数；
步骤10，根据对数比值图像Xb中每个像素点的一元势能函数和二元势能函数，得到初步算法模型p(X|Y)，其中，X为标记场，Y为观测场；
所述初步算法模型p(X|Y)的表达式为：
其中，X为标记场，Y为观测场，U为非平稳场，1/Z为初步算法模型p(X|Y)的分布函数，Qn为对数比值图像Xb中第n个像素点的邻域系统，xn为标记场X中对应于对数比值图像Xb的第n个像素点的标记值，xt为标记场X中对应于对数比值图像Xb的像素点t的标记值，n∈{1，2，...，N}，N为对数比值图像Xb的像素点总数，A(xn，Fn，Ub，X)为对数比值图像Xb中第n个像素点的一元势能函数，I(xn，xt，X，Y，U)为对数比值图像Xb中第n个像素点的二元势能函数：
QH为水平邻域系统、QV为垂直邻域系统，rn为对数比值图像Xb中第n个像素点的边界强度，rt为对数比值图像Xb中像素点t的边界强度，edge_C为常量，为二元势能函数参数，Y为观测场，U为非平稳场，g1，g2为非平稳场U的两个不同取值，一般分别取0和1，Un为对数比值图像Xb中第n个像素点的平稳性标记值，Ut为对数比值图像Xb中像素点t的平稳性标记值，δ为阶跃函数，其中：
步骤11，利用条件迭代期望算法计算初步算法模型p(X|Y)中的二元势能函数参数，根据该二元势能函数参数求得对数比值图像Xb的最终分类标签矩阵，即得到原始第二时相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图像X1和原始第二时相SAR图像X2。
输入的待检测的原始第一时相SAR图像X1和原始第二时相SAR图像X2如下：
其中，和分别为原始第一时相SAR图像的第n个像素点的坐标和灰度值，和分别为原始第二时相SAR图像的第n个像素点的坐标和灰度值；I为原始第一时相SAR图像的长度，也是原始第二时相SAR图像的长度；J为原始第一时相SAR图像的宽度，也是原始第二时相SAR图像的宽度；n∈{1，2，...，N}，N为原始第一时相SAR图像的像素点总数，也是原始第二时相SAR图像的像素点总数。
步骤2，对原始第一时相SAR图像X1和原始第二时相SAR图像X2分别进行归一化处理，并对归一化后的第一时相SAR图像X′1和归一化后的第二时相SAR图像X′2作对数比值运算，得到对数比值图像Xb；对数比值图像Xb的大小与原始第一、第二时相SAR图像的大小均相同。
步骤2的具体子步骤为：
2a)采用以下公式分别对原始第一时相SAR图像X1和原始第二时相SAR图像X2进行归一化处理，得到归一化后的第一时相SAR图像X′1和归一化后的第二时相SAR图像X′2：
其中，min表示取最小值，max表示取最大值；
2b)采用以下公式对归一化后的第一时相SAR图像X′1和归一化后的第二时相SAR图像X′2作对数比值运算，得到对数比值图像Xb：
其中，log表示取对数。
步骤3，利用滑动窗口分别提取归一化后的第一时相SAR图像的纹理特征矩阵和归一化后的第二时相SAR图像的纹理特征矩阵；并将归一化后的第二时相SAR图像的纹理特征矩阵减去归一化后的第一时相SAR图像的纹理特征矩阵，得到对数比值图像的纹理特征矩阵；所述纹理特征矩阵依次由均值特征矩阵、对比度特征矩阵和半方差特征矩阵组合而成。
步骤3的具体子步骤为：
3a)采用以下公式分别计算归一化后的第一时相SAR图像的第n个像素点的均值特征和归一化后的第二时相SAR图像的第n个像素点的均值特征
其中，和分别为归一化后的第一时相SAR图像的第n个像素点的坐标和灰度值，和为归一化后的第二时相SAR图像的第n个像素点的坐标和灰度值；I为归一化后的第一时相SAR图像的长度，也是归一化后的第二时相SAR图像的长度，J为归一化后的第一时相SAR图像的宽度，也是归一化后的第二时相SAR图像的宽度；n∈{1，2，...，N}，N为归一化后的第一时相SAR图像的像素点总数，也是归一化后的第二时相SAR图像的像素点总数；n0为滑动窗口包含的像素点数；本发明实例中，滑动窗口的尺寸设定为5×5，则n为25；
归一化后的第一时相SAR图像的均值特征矩阵μ1和归一化后的第二时相SAR图像的均值特征矩阵μ2分别为：
3b)采用以下公式分别计算归一化后的第一时相SAR图像的第n个像素点的对比度特征和归一化后的第二时相SAR图像的第n个像素点的对比度特征
归一化后的第一时相SAR图像的对比度特征矩阵cs1和归一化后的第二时相SAR图像的对比度特征矩阵cs2分别为：
3c)采用以下公式分别计算归一化后的第一时相SAR图像的第n个像素点的半方差特征和归一化后的第二时相SAR图像的第n个像素点的半方差特征
其中，为归一化后的第一时相SAR图像在像素点的灰度值，为归一化后的第二时相SAR图像在像素点处的灰度值；归一化后的第一时相SAR图像的像素点与归一化后的第一时相SAR图像的第n个像素点的欧几里得距离为h，归一化后的第二时相SAR图像的像素点与归一化后的第二时相SAR图像的第n个像素点的欧几里得距离为h；
归一化后的第一时相SAR图像的半方差特征矩阵χ1和归一化后的第二时相SAR图像的半方差特征矩阵χ2分别为：
3d)将归一化后的第二时相SAR图像的均值特征矩阵μ2减去归一化后的第一时相SAR图像的均值特征矩阵μ1，得到对数比值图像的均值特征矩阵μb：
其中，和分别为对数比值图像的第n个像素点的坐标和均值特征；3e)将归一化后的第二时相SAR图像的对比度特征矩阵cs2减去归一化后的第一时相SAR图像的对比度特征矩阵cs1，得到对数比值图像的对比度特征矩阵csb：
其中，为对数比值图像的第n个像素点的对比度特征；
3f)将归一化后的第二时相SAR图像的半方差特征矩阵χ2减去归一化后的第一时相SAR图像的半方差特征矩阵χ1，得到对数比值图像的半方差特征矩阵χb：
其中，为对数比值图像的第n个像素点的半方差特征；
3g)在滑动窗口将对数比值图像的均值特征矩阵、对比度特征矩阵和半方差特征矩阵，构成对数比值图像的纹理特征矩阵gs，gs＝[μb，cab，χb]。
步骤4，利用三重马尔科夫场模型算法，将对数比值图像Xb进行平稳性分割，得到两种平稳性区域：A平稳区域与B平稳区域；A平稳区域中像素点的个数与B平稳区域中像素点的个数的和等于对数比值图像中像素点的个数；A平稳区域与B平稳区域的平稳性不同，体现为A平稳区域和B平稳区域中像素点的纹理特征不同。所述纹理特征依次包含均值特征、对比度特征和半方差特征。
需要说明的是，纹理特征是SAR图像的众多特征中具有较高区分度的特征，它反映了SAR图像中目标的表面结构组织排列、空间结构关系和相互邻域关系。在SAR图像的不同平稳区域的成像场景中，成像目标的表面结构组织排列、空间结构关系和相互邻域关系是有区别的，所以我们可根据SAR图像不同平稳区域的纹理特征间的差别，对SAR图像进行平稳性分割。SAR图像具有非平稳性，可表现出同一平稳区域的纹理特征相同、不同平稳区域的纹理特征不同的特点。本发明实例中，A平稳区域与B平稳区域的区别在于：A平稳区域中像素点的纹理特征和B平稳区域中像素点的纹理特征不同。本发明实例中，采用三重马尔科夫场模型算法(Triplet Markov Field，TMF)将对数比值图像进行平稳性分割，具体地，采用TMF引入非平稳场U，用U的不同取值来描述两种不同的平稳性分布，结合马尔科夫场对对数比值图像进行平稳性分割；具体实现时，也可以采用其他方法对对数比值图像进行平稳性分割，本发明实例对此不做限制。
步骤5，分别选出A平稳区域和B平稳区域的样本点。
步骤5的具体子步骤为：
5a)在对数比值图像中选取N0个原始第二时相SAR图像相对于原始第一时相SAR图像在视觉上未变化的像素点，再选取N1个原始第二时相SAR图像相对于原始第一时相SAR图像在视觉上变化的像素点，N0+N1＜N，N为原始第一时相SAR图像的像素点总数，也是原始第二时相SAR图像的像素点总数；
5b)根据所选出的每个在视觉上未变化的像素点的坐标，判断其位于A平稳区域还是B平稳区域，将所有位于A平稳区域的在视觉上未变化的像素点构成第一集合KA0，该第一集合KA0的大小为NA0，将所有位于B平稳区域的在视觉上未变化的像素点构成第二集合KB0，该第二集合KB0的大小为NB0；
根据所选出的每个在视觉上变化的像素点的坐标，判断其位于A平稳区域还是B平稳区域，将所有位于A平稳区域的在视觉上变化的像素点构成第三集合KA1，该第三集合KA1的大小为NA1，将所有位于B平稳区域的在视觉上未变化的像素点构成第四集合KB1，该第四集合KB1的大小为NB1；
将第一集合KA0和第三集合KA1中的像素点作为A平稳区域的样本点；将第二集合KB0和第四集合KB1中的像素点作为B平稳区域的样本点；
其中，N0＝NA0+NB0，N1＝NA1+NB1。
步骤6，提取A平稳区域的所有样本点的纹理特征，组成A平稳区域的纹理特征向量，给定A平稳区域的训练标签lA，结合A平稳区域的训练标签lA和A平稳区域的纹理特征向量构成A平稳区域的训练样本，利用该训练样本训练支撑向量机，得到第一分类标签矩阵Ac。
步骤6的具体子步骤为：
6a)根据对数比值图像的纹理特征矩阵gs，提取第一集合KA0和第三集合KA1的像素点的纹理特征，组成A平稳区域的纹理特征向量gA：
其中，分别为A平稳区域中在视觉上未变化的NA0个像素点中第1个像素点的均值特征、对比度特征和半方差特征；分别为A平稳区域中在视觉上未变化的NA0个像素点中第NA0个像素点的均值特征、对比度特征和半方差特征；分别为A平稳区域中在视觉上变化的NA1个像素点中第1个像素点的均值特征、对比度特征和半方差特征； 分别为A平稳区域中在视觉上变化的NA1个像素点中第NA1个像素点的均值特征、对比度特征和半方差特征；
再结合A平稳区域的训练标签lA构成A平稳区域的训练样本SpA：
SpA＝{gA，lA}
其中，A平稳区域的训练标签lA是一个长度为N的行向量，其前N0个元素为0，其后N1个元素为1：
6b)将支撑向量机的类型设置为C-SVC模型，并将核函数的类型设置为径向基核函数；
需要说明的是，支撑向量机是一种二类分类模型，其基本定义为特征空间上的间隔最大的线性分类器，其学习策略是间隔最大化；本发明实例中，设置支撑向量机参数s＝0，即将支撑向量机的类型设置为C-SVC模型，设置支撑向量机参数t＝2，即将核函数的类型设置为径向基核函数；
6c)将惩罚因子C和径向基核函数中核参数γ进行交叉验证，得到对应于A平稳区域的最优(C，γ)组合；
具体地，选取惩罚因子C∈{2-8，2-7.5，...，27.5，28}，径向基核函数中核参数γ∈{2-8，2-7.5，...，27.5，28}，则共有33×33个(C，γ)组合，对该33×33个(C，γ)组合进行5层交叉验证，具体步骤为：
将A平稳区域的训练样本SpA随机地平均分成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平稳区域的训练样本SpA对支撑向量机进行训练，得到第一分类标签矩阵Ac，其维数与对数比值图像Xb的维数相同。
步骤7，提取B平稳区域的所有样本点的纹理特征，组成B平稳区域的纹理特征向量，给定B平稳区域的训练标签lB，结合B平稳区域的训练标签lB和B平稳区域的纹理特征向量构成B平稳区域的训练样本，利用该训练样本训练支撑向量机，得到第二分类标签矩阵Bc。
步骤7的具体子步骤为：
7a)根据对数比值图像的纹理特征矩阵gs，提取第二集合KB0和第四集合KB1的像素点的纹理特征，组成B平稳区域的纹理特征向量gB：
其中，分别为B平稳区域中在视觉上未变化的NB0个像素点中第1个像素点的均值特征、对比度特征和半方差特征； 分别为B平稳区域中在视觉上未变化的NB0个像素点中第NB0个像素点的均值特征、对比度特征和半方差特征；分别为B平稳区域中在视觉上变化的NB1个像素点中第1个像素点的均值特征、对比度特征和半方差特征；分别为B平稳区域中在视觉上变化的NB1个像素点中第NB1个像素点的均值特征、对比度特征和半方差特征；
再结合B平稳区域的训练标签lB构成B平稳区域的训练样本SpB：
SpB＝{gB，lB}
其中，B平稳区域的训练标签lB是一个长度为N的行向量，其前N0个元素为0，其后N1个元素为1：
7b)将支撑向量机的类型设置为C-SVC模型，并将核函数的类型设置为径向基核函数；
本发明实例中，设置支撑向量机参数s＝0，即将支撑向量机的类型设置为C-SVC模型，设置支撑向量机参数t＝2，即将核函数的类型设置为径向基核函数；
7c)将惩罚因子C和径向基核函数中核参数γ进行交叉验证，得到对应于A平稳区域的最优的(C，γ)组合；
具体地，选取惩罚因子C∈{2-8，2-7.5，...，27.5，28}，径向基核函数中核参数γ∈{2-8，2-7.5，...，27.5，28}，则共有33×33个(C，γ)组合，对该33×33个(C，γ)组合进行5层交叉验证，具体步骤为：
将B平稳区域的训练样本SpB随机地平均分成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平稳区域的训练样本SpB对支撑向量机进行训练，得到第二分类标签矩阵Bc，其维数与对数比值图像Xb的维数相同。
步骤8，根据第一分类标签矩阵Ac和第二分类标签矩阵Bc，求得对数比值图像Xb的初始分类标签矩阵Oc，并求得对数比值图像Xb中每个像素点的类条件概率，再对对数比值图像Xb中每个像素点的类条件概率取对数，得到对数比值图像Xb中每个像素点的一元势能函数。
步骤8的具体子步骤为：
8a)构建与原始第一时相SAR图像维数相同的训练标签矩阵Oc′，训练标签矩阵Oc′的初始值为全0矩阵，将第一分类标签矩阵Ac中与A平稳区域中的每个像素点分别对应的元素的值，按A平稳区域中的每个像素点的坐标位置，分别填入训练标签矩阵Oc′；将第二分类标签矩阵Bc中与B平稳区域中的每个像素点分别对应的元素的值，按B平稳区域中的每个像素点的坐标位置，分别填入训练标签矩阵Oc′；将训练标签矩阵Oc′作为对数比值图像Xb的初始分类标签矩阵Oc，对数比值图像Xb的初始分类标签矩阵Oc中元素的取值为0或1，0表示原始第二时相SAR图像中的对应像素点相对于原始第一时相SAR图像中的对应像素点未发生变化，1表示原始第二时相SAR图像中的对应像素点相对于原始第一时相SAR图像中的对应像素点发生变化；
8b)将对数比值图像Xb的初始分类标签矩阵Oc拟合为S形生长曲线函数，即求得对数比值图像Xb中每个像素点的类条件概率，并对对数比值图像Xb中每个像素点的类条件概率取对数，得到对数比值图像Xb中每个像素点的一元势能函数，其中，对数比值图像Xb中第n个像素点的一元势能函数A(xn，Fn，Un，X)为：
A(xn，Fn，Un，X)＝log p(xn|dn(F，U))
其中，p(xn|dn(F，U))为对数比值图像Xb中第n个像素点的类条件概率，X为标记场，xn为标记场X中对应于对数比值图像Xb的第n个像素点的标记值，Fn为对数比值图像Xb中第n个像素点的纹理特征向量， 分别为对数比值图像Xb中第n个像素点的均值特征、对比度特征和半方差特征，为对数比值图像Xb中第n个像素点的坐标，Un为对数比值图像Xb中第n个像素点的平稳性标记值，F为纹理特征场，U为非平稳场，dn(F，U)为对数比值图像Xb中第n个像素点的纹理特征场与非平稳场的混合场，n∈{1，2，...，N}，N为对数比值图像Xb的像素点总数，log表示取对数。
需要说明的是，S形生长曲线函数即Sigmoid函数，Sigmoid函数是一个在生物学中常见的S型函数，也称为S形生长曲线。
步骤9，利用指数加权均值比率算子提取对数比值图像Xb中每个像素点的边界强度，并构建对数比值图像Xb中每个像素点的二元势能函数。
指数加权均值比率算子(ROEWA)是基于线性最小均方误差准则的指数滤波器，其计算结果为经过指数加权处理后的均值；
步骤9的具体子步骤为：
9a)定义平滑函数f(λ)由因果滤波器f1(λ)和非因果滤波器f2(λ)组成：
其中，f1(λ)＝cdλu(λ)，f2(λ-1)＝cd-(λ-1)u(-(λ-1))，d为设定常数，且0＜d＜1，c＝1-d，u(λ)为阶跃函数，λ为自变量；
9b)给出指数加权均值比率算子的定义如下：
其中，表示对数比值图像Xb中第n个像素点的坐标，n∈{1，2，...，N}，N为对数比值图像Xb的像素点总数，μI1，μI2，μJ1，μJ2为指数加权值，按下式进行计算：
其中，Fn为对数比值图像Xb中第n个像素点的纹理特征向量，分别为对数比值图像Xb中第n个像素点的均值特征、对比度特征和半方差特征，*表示水平方向上的卷积，则表示垂直方向上的卷积；
9c)通过以下公式提取对数比值图像Xb中第n个像素点的边界强度rn：
9c)通过以下公式构建对数比值图像Xb中第n个像素点的二元势能函数I(xn，xt，X，Y，U)：
其中，QH为水平邻域系统、QV为垂直邻域系统，rn为对数比值图像Xb中第n个像素点的边界强度，rt为对数比值图像Xb中像素点t的边界强度，X为标记场，xn为标记场X中对应于对数比值图像Xb的第n个像素点的标记值，xt为标记场X中对应于对数比值图像Xb的像素点t的标记值，edge_C为常量，为二元势能函数参数，Y为观测场，U为非平稳场，g1，g2为非平稳场U的两个不同取值，一般分别取0和1，Un为对数比值图像Xb中第n个像素点的平稳性标记值，U1为对数比值图像Xb中像素点t的平稳性标记值，δ为阶跃函数，其中：
步骤10，根据对数比值图像Xb中每个像素点的一元势能函数和二元势能函数，得到初步算法模型p(X|Y)，其中，X为标记场，Y为观测场。
所述初步算法模型p(X|Y)的表达式为：
其中，X为标记场，Y为观测场，U为非平稳场，1/Z为初步算法模型p(X|Y)的分布函数，Qn为对数比值图像Xb中第n个像素点的邻域系统，xn为标记场X中对应于对数比值图像Xb的第n个像素点的标记值，xt为标记场X中对应于对数比值图像Xb的像素点t的标记值，n∈{1，2，...，N}，N为对数比值图像Xb的像素点总数，A(xn，Fn，Un，X)对数比值图像Xb中第n个像素点的一元势能函数，I(xn，xt，X，Y，U)为对数比值图像Xb中第n个像素点的二元势能函数：
QH为水平邻域系统、QV为垂直邻域系统，rn为对数比值图像Xb中第n个像素点的边界强度，rt为对数比值图像Xb中像素点t的边界强度，edge_C为常量，为二元势能函数参数，Y为观测场，U为非平稳场，g1，g2为非平稳场U的两个不同取值，一般分别取0和1，Un为对数比值图像Xb中第n个像素点的平稳性标记值，Ut为对数比值图像Xb中像素点t的平稳性标记值，δ为阶跃函数，其中：
步骤11，利用条件迭代期望算法计算初步算法模型p(X|Y)中的二元势能函数参数，根据该二元势能函数参数求得对数比值图像Xb的最终分类标签矩阵，即得到原始第二时相SAR图像相对于原始第一时相SAR图像的变化检测结果。
步骤11的具体子步骤为：
11a)在对数比值图像Xb的初始分类标签矩阵Oc中随机选取出一个大小为M×M的子矩阵作为初步算法模型p(X|Y)中的标记场X的初始值X(0)，并根据初步算法模型p(X|Y)求得二元势能函数参数的初始值，再利用所求出的二元势能函数参数的初始值，根据初步算法模型p(X|Y)求得第1代标记场X1作为第1代分类标签矩阵；设置迭代次数t＝1；
11b)在第t代分类标签矩阵中随机取出一个大小为M×M的子矩阵作为初步算法模型p(X|Y)中的第t次迭代的标记场X(t)，并根据初步算法模型p(X|Y)求得第t次迭代的二元势能函数参数，再利用所求出的第t次迭代的二元势能函数参数，根据初步算法模型p(X|Y)求得第t+1代标记场Xt+1作为第t+1代分类标签矩阵；
11c)如果迭代次数t＞9，则将第t+1代分类标签矩阵作为对数比值图像Xb的最终分类标签矩阵，即将第t+1代分类标签矩阵作为原始第二时相SAR图像相对于原始第一时相SAR图像的变化检测结果；反之，令迭代次数t增加1，返回步骤11b)。
本发明的效果可通过以下仿真实验作进一步证实：
1)实验条件：
仿真实验的环境为：MATLAB R20111b，Libsvm-3.17工具箱，Intel(R)Core(TM)i5-3570 CPU 3.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图像的实际情况。
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (9)

1.一种基于非平稳分析与条件随机场的SAR图像变化检测方法，其特征在于，包括以下步骤：
步骤1，输入待检测的原始第一时相SAR图像X1和原始第二时相SAR图像X2；
步骤2，对原始第一时相SAR图像X1和原始第二时相SAR图像X2分别进行归一化处理，并对归一化后的第一时相SAR图像X′1和归一化后的第二时相SAR图像X′2作对数比值运算，得到对数比值图像Xb；对数比值图像Xb的大小与原始第一、第二时相SAR图像的大小均相同；
步骤3，利用滑动窗口分别提取归一化后的第一时相SAR图像的纹理特征矩阵和归一化后的第二时相SAR图像的纹理特征矩阵；并将归一化后的第二时相SAR图像的纹理特征矩阵减去归一化后的第一时相SAR图像的纹理特征矩阵，得到对数比值图像的纹理特征矩阵；所述纹理特征矩阵依次由均值特征矩阵、对比度特征矩阵和半方差特征矩阵组合而成；
步骤4，将对数比值图像Xb进行平稳性分割，得到两种平稳性区域：A平稳区域与B平稳区域；A平稳区域中像素点的个数与B平稳区域中像素点的个数的和等于对数比值图像中像素点的个数；
步骤5，分别选出A平稳区域和B平稳区域的样本点；
步骤6，提取A平稳区域的所有样本点的纹理特征，组成A平稳区域的纹理特征向量，给定A平稳区域的训练标签lA，结合A平稳区域的训练标签lA和A平稳区域的纹理特征向量构成A平稳区域的训练样本，利用该训练样本训练支撑向量机，得到第一分类标签矩阵Ac；
步骤7，提取B平稳区域的所有样本点的纹理特征，组成B平稳区域的纹理特征向量，给定B平稳区域的训练标签lB，结合B平稳区域的训练标签lB和B平稳区域的纹理特征向量构成B平稳区域的训练样本，利用该训练样本训练支撑向量机，得到第二分类标签矩阵Bc；
步骤8，根据第一分类标签矩阵Ac和第二分类标签矩阵Bc，求得对数比值图像Xb的初始分类标签矩阵Oc，并求得对数比值图像Xb中每个像素点的类条件概率，再对对数比值图像Xb中每个像素点的类条件概率取对数，得到对数比值图像Xb中每个像素点的一元势能函数；
步骤9，提取对数比值图像Xb中每个像素点的边界强度，并构建对数比值图像Xb中每个像素点的二元势能函数；
步骤10，根据对数比值图像Xb中每个像素点的一元势能函数和二元势能函数，得到初步算法模型p(X|Y)，其中，X为标记场，Y为观测场；
所述初步算法模型p(X|Y)的表达式为：
其中，U为非平稳场，1/Z为初步算法模型p(X|Y)的分布函数，Qn为对数比值图像Xb中第n个像素点的邻域系统，xn为标记场X中对应于对数比值图像Xb的第n个像素点的标记值，xt为标记场X中对应于对数比值图像Xb的像素点t的标记值，n∈{1，2，...，N}，N为对数比值图像Xb的像素点总数，A(xn，Fn，Un，X)为对数比值图像Xb中第n个像素点的一元势能函数，Fn为对数比值图像Xb中第n个像素点的纹理特征向量，I(xn，xt，X，Y，U)为对数比值图像Xb中第n个像素点的二元势能函数：
QH为水平邻域系统、QV为垂直邻域系统，rn为对数比值图像Xb中第n个像素点的边界强度，rt为对数比值图像Xb中像素点t的边界强度，edge_C为常量，为二元势能函数参数，g1，g2为非平稳场U的两个不同取值，分别取0和1，Un为对数比值图像Xb中第n个像素点的平稳性标记值，Ut为对数比值图像Xb中像素点t的平稳性标记值，δ为阶跃函数，其中：
步骤11，计算初步算法模型p(X|Y)中的二元势能函数参数，根据该二元势能函数参数求得对数比值图像Xb的最终分类标签矩阵，即得到原始第二时相SAR图像相对于原始第一时相SAR图像的变化检测结果。
2.如权利要求1所述的基于非平稳分析与条件随机场的SAR图像变化检测方法，其特征在于，所述步骤2的具体子步骤为：
2a)采用以下公式分别对原始第一时相SAR图像X1和原始第二时相SAR图像X2进行归一化处理，得到归一化后的第一时相SAR图像X′1和归一化后的第二时相SAR图像X′2：
其中，min表示取最小值，max表示取最大值；
2b)采用以下公式对归一化后的第一时相SAR图像X′1和归一化后的第二时相SAR图像X′2作对数比值运算，得到对数比值图像Xb：
其中，log表示取对数。
3.如权利要求1所述的基于非平稳分析与条件随机场的SAR图像变化检测方法，其特征在于，所述步骤3的具体子步骤为：
3a)采用以下公式分别计算归一化后的第一时相SAR图像的第n个像素点的均值特征和归一化后的第二时相SAR图像的第n个像素点的均值特征
其中，和分别为归一化后的第一时相SAR图像的第n个像素点的坐标和灰度值，和为归一化后的第二时相SAR图像的第n个像素点的坐标和灰度值；I为归一化后的第一时相SAR图像的长度，也是归一化后的第二时相SAR图像的长度，J为归一化后的第一时相SAR图像的宽度，也是归一化后的第二时相SAR图像的宽度；n0为滑动窗口包含的像素点数；
归一化后的第一时相SAR图像的均值特征矩阵μ1和归一化后的第二时相SAR图像的均值特征矩阵μ2分别为：
3b)采用以下公式分别计算归一化后的第一时相SAR图像的第n个像素点的对比度特征和归一化后的第二时相SAR图像的第n个像素点的对比度特征
归一化后的第一时相SAR图像的对比度特征矩阵cs1和归一化后的第二时相SAR图像的对比度特征矩阵cs2分别为：
3c)采用以下公式分别计算归一化后的第一时相SAR图像的第n个像素点的半方差特征和归一化后的第二时相SAR图像的第n个像素点的半方差特征
归一化后的第一时相SAR图像的半方差特征矩阵χ1和归一化后的第二时相SAR图像的半方差特征矩阵χ2分别为：
3d)将归一化后的第二时相SAR图像的均值特征矩阵μ2减去归一化后的第一时相SAR图像的均值特征矩阵μ1，得到对数比值图像的均值特征矩阵μb：
其中，和分别为对数比值图像的第n个像素点的坐标和均值特征；
3e)将归一化后的第二时相SAR图像的对比度特征矩阵cs2减去归一化后的第一时相SAR图像的对比度特征矩阵cs1，得到对数比值图像的对比度特征矩阵csb：
其中，为对数比值图像的第n个像素点的对比度特征；
3f)将归一化后的第二时相SAR图像的半方差特征矩阵χ2减去归一化后的第一时相SAR图像的半方差特征矩阵χ1，得到对数比值图像的半方差特征矩阵χb：
其中，为对数比值图像的第n个像素点的半方差特征；
3g)在滑动窗口将对数比值图像的均值特征矩阵、对比度特征矩阵和半方差特征矩阵，构成对数比值图像的纹理特征矩阵gs，gs＝[μb，csb，χb]。
4.如权利要求1所述的基于非平稳分析与条件随机场的SAR图像变化检测方法，其特征在于，所述步骤5的具体子步骤为：
5a)在对数比值图像中选取N0个原始第二时相SAR图像相对于原始第一时相SAR图像在视觉上未变化的像素点，再选取N1个原始第二时相SAR图像相对于原始第一时相SAR图像在视觉上变化的像素点，N0+N1＜N；
5b)根据所选出的每个在视觉上未变化的像素点的坐标，判断其位于A平稳区域还是B平稳区域，将所有位于A平稳区域的在视觉上未变化的像素点构成第一集合KA0，该第一集合KA0的大小为NA0，将所有位于B平稳区域的在视觉上未变化的像素点构成第二集合KB0，该第二集合KB0的大小为NB0；
根据所选出的每个在视觉上变化的像素点的坐标，判断其位于A平稳区域还是B平稳区域，将所有位于A平稳区域的在视觉上变化的像素点构成第三集合KA1，该第三集合KA1的大小为NA1，将所有位于B平稳区域的在视觉上变化的像素点构成第四集合KB1，该第四集合KB1的大小为NB1；
将第一集合KA0和第三集合KA1中的像素点作为A平稳区域的样本点；将第二集合KB0和第四集合KB1中的像素点作为B平稳区域的样本点；
其中，N0＝NA0+NB0，N1＝NA1+NB1。
5.如权利要求4所述的基于非平稳分析与条件随机场的SAR图像变化检测方法，其特征在于，所述步骤6的具体子步骤为：
6a)根据对数比值图像的纹理特征矩阵gs，提取第一集合KA0和第三集合KA1的像素点的纹理特征，组成A平稳区域的纹理特征向量gA：
其中，分别为A平稳区域中在视觉上未变化的NA0个像素点中第1个像素点的均值特征、对比度特征和半方差特征；分别为A平稳区域中在视觉上未变化的NA0个像素点中第NA0个像素点的均值特征、对比度特征和半方差特征；分别为A平稳区域中在视觉上变化的NA1个像素点中第1个像素点的均值特征、对比度特征和半方差特征；分别为A平稳区域中在视觉上变化的NA1个像素点中第NA1个像素点的均值特征、对比度特征和半方差特征；
再结合A平稳区域的训练标签lA构成A平稳区域的训练样本SpA：
SpA＝{gA，lA}
其中，A平稳区域的训练标签lA是一个长度为N的行向量，其前N0个元素为0，其后N1个元素为1：
6b)将支撑向量机的类型设置为C-SVC模型，并将核函数的类型设置为径向基核函数；
6c)将惩罚因子C和径向基核函数中核参数γ进行交叉验证，得到对应于A平稳区域的最优(C，γ)组合；
选取惩罚因子C∈{2-8，2-7.5，...，27.5，28}，径向基核函数中核参数γ∈{2-8，2-7.5，...，27.5，28}，则共有33×33个(C，γ)组合，对该33×33个(C，γ)组合进行5层交叉验证，具体步骤为：
将A平稳区域的训练样本SpA随机地平均分成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平稳区域的训练样本SpA对支撑向量机进行训练，得到第一分类标签矩阵Ac，其维数与对数比值图像Xb的维数相同。
6.如权利要求4所述的基于非平稳分析与条件随机场的SAR图像变化检测方法，其特征在于，所述步骤7的具体子步骤为：
7a)根据对数比值图像的纹理特征矩阵gs，提取第二集合KB0和第四集合KB1的像素点的纹理特征，组成B平稳区域的纹理特征向量gB：
其中，分别为B平稳区域中在视觉上未变化的NB0个像素点中第1个像素点的均值特征、对比度特征和半方差特征；分别为B平稳区域中在视觉上未变化的NB0个像素点中第NB0个像素点的均值特征、对比度特征和半方差特征；分别为B平稳区域中在视觉上变化的NB1个像素点中第1个像素点的均值特征、对比度特征和半方差特征；分别为B平稳区域中在视觉上变化的NB1个像素点中第NB1个像素点的均值特征、对比度特征和半方差特征；
再结合B平稳区域的训练标签lB构成B平稳区域的训练样本SpB：
SpB＝{gB，lB}
其中，B平稳区域的训练标签lB是一个长度为N的行向量，其前N0个元素为0，其后N1个元素为1：
7b)将支撑向量机的类型设置为C-SVC模型，并将核函数的类型设置为径向基核函数；
7c)将惩罚因子C和径向基核函数中核参数γ进行交叉验证，得到对应于B平稳区域的最优的(C，γ)组合；
选取惩罚因子C∈{2-8，2-7.5，...，27.5，28}，径向基核函数中核参数γ∈{2-8，2-7.5，...，27.5，28}，则共有33×33个(C，γ)组合，对该33×33个(C，γ)组合进行5层交叉验证，具体步骤为：
将B平稳区域的训练样本SpB随机地平均分成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平稳区域的训练样本SpB对支撑向量机进行训练，得到第二分类标签矩阵Bc，其维数与对数比值图像Xb的维数相同。
7.如权利要求1所述的基于非平稳分析与条件随机场的SAR图像变化检测方法，其特征在于，所述步骤8的具体子步骤为：
8a)构建与原始第一时相SAR图像维数相同的训练标签矩阵Oc′，训练标签矩阵Oc′的初始值为全0矩阵，将第一分类标签矩阵Ac中与A平稳区域中的每个像素点分别对应的元素的值，按A平稳区域中的每个像素点的坐标位置，分别填入训练标签矩阵Oc′；将第二分类标签矩阵Bc中与B平稳区域中的每个像素点分别对应的元素的值，按B平稳区域中的每个像素点的坐标位置，分别填入训练标签矩阵Oc′；将训练标签矩阵Oc′作为对数比值图像Xb的初始分类标签矩阵Oc，对数比值图像Xb的初始分类标签矩阵Oc中元素的取值为0或1，0表示原始第二时相SAR图像中的对应像素点相对于原始第一时相SAR图像中的对应像素点未发生变化，1表示原始第二时相SAR图像中的对应像素点相对于原始第一时相SAR图像中的对应像素点发生变化；
8b)将对数比值图像Xb的初始分类标签矩阵Oc拟合为S形生长曲线函数，即求得对数比值图像Xb中每个像素点的类条件概率，并对对数比值图像Xb中每个像素点的类条件概率取对数，得到对数比值图像Xb中每个像素点的一元势能函数，其中，对数比值图像Xb中第n个像素点的一元势能函数A(xn，Fn，Un，X)为：
A(xn，Fn，Un，X)＝logp(xn|dn(F，U))
其中，p(xn|dn(F，U))为对数比值图像Xb中第n个像素点的类条件概率，Fn为对数比值图像Xb中第n个像素点的纹理特征向量， 分别为对数比值图像Xb中第n个像素点的均值特征、对比度特征和半方差特征，为对数比值图像Xb中第n个像素点的坐标，F为纹理特征场，dn(F，U)为对数比值图像Xb中第n个像素点的纹理特征场与非平稳场的混合场，log表示取对数。
8.如权利要求1所述的基于非平稳分析与条件随机场的SAR图像变化检测方法，其特征在于，所述步骤9的具体子步骤为：
9a)定义平滑函数f(λ)由因果滤波器f1(λ)和非因果滤波器f2(λ)组成：
其中，f1(λ)＝cdλu(λ)，f2(λ-1)＝cd-(λ-1)u(-(λ-1))，d为设定常数，且0＜d＜1，c＝1-d，u(λ)为阶跃函数，λ为自变量；
9b)给出指数加权均值比率算子的定义如下：
其中，表示对数比值图像Xb中第n个像素点的坐标，μI1，μI2，μJ1，μJ2为指数加权值，按下式进行计算：
其中，Fn为对数比值图像Xb中第n个像素点的纹理特征向量， 分别为对数比值图像Xb中第n个像素点的均值特征、对比度特征和半方差特征，*表示水平方向上的卷积，则表示垂直方向上的卷积；
9c)通过以下公式提取对数比值图像Xb中第n个像素点的边界强度rn：
9c)通过以下公式构建对数比值图像Xb中第n个像素点的二元势能函数I(xn，xt，X，Y，U)：
9.如权利要求1所述的基于非平稳分析与条件随机场的SAR图像变化检测方法，其特征在于，所述步骤11的具体子步骤为：
11a)在对数比值图像Xb的初始分类标签矩阵Oc中随机选取出一个大小为M×M的子矩阵作为初步算法模型p(X|Y)中的标记场X的初始值X(0)，并根据初步算法模型p(X|Y)求得二元势能函数参数的初始值，再利用所求出的二元势能函数参数的初始值，根据初步算法模型p(X|Y)求得第1代标记场X1作为第1代分类标签矩阵；设置迭代次数k＝1：
11b)在第k代分类标签矩阵中随机取出一个大小为M×M的子矩阵作为初步算法模型p(X|Y)中的第k次迭代的标记场X(k)，并根据初步算法模型p(X|Y)求得第k次迭代的二元势能函数参数，再利用所求出的第k次迭代的二元势能函数参数，根据初步算法模型p(X|Y)求得第k+1代标记场Xk+1作为第k+1代分类标签矩阵；
11c)如果迭代次数k＞9，则将第k+1代分类标签矩阵作为对数比值图像Xb的最终分类标签矩阵，即将第k+1代分类标签矩阵作为原始第二时相SAR图像相对于原始第一时相SAR图像的变化检测结果；反之，令迭代次数k增加1，返回步骤11b)。
CN201510526592.5A 2015-08-25 2015-08-25 基于非平稳分析与条件随机场的sar图像变化检测方法 Expired - Fee Related CN105160666B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510526592.5A CN105160666B (zh) 2015-08-25 2015-08-25 基于非平稳分析与条件随机场的sar图像变化检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510526592.5A CN105160666B (zh) 2015-08-25 2015-08-25 基于非平稳分析与条件随机场的sar图像变化检测方法

Publications (2)

Publication Number Publication Date
CN105160666A CN105160666A (zh) 2015-12-16
CN105160666B true CN105160666B (zh) 2018-03-06

Family

ID=54801508

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510526592.5A Expired - Fee Related CN105160666B (zh) 2015-08-25 2015-08-25 基于非平稳分析与条件随机场的sar图像变化检测方法

Country Status (1)

Country Link
CN (1) CN105160666B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109191503B (zh) * 2018-08-23 2021-08-27 河海大学 基于条件随机场的遥感影像变化检测方法及系统
CN109242832B (zh) * 2018-08-23 2021-08-27 河海大学 一种多时相多光谱遥感影像变化检测方法及系统
CN109376787B (zh) * 2018-10-31 2021-02-26 聚时科技(上海)有限公司 流形学习网络及基于其的计算机视觉图像集分类方法
CN113326724B (zh) * 2021-02-07 2024-02-02 海南长光卫星信息技术有限公司 一种遥感影像变化检测方法、装置、设备及可读存储介质

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101634705A (zh) * 2009-08-19 2010-01-27 西安电子科技大学 基于方向信息测度的sar图像目标变化检测方法
EP2431764A1 (en) * 2010-09-17 2012-03-21 BAE Systems PLC Processing SAR imagery
CN102779346A (zh) * 2012-07-05 2012-11-14 西安电子科技大学 基于改进c-v模型的sar图像变化检测方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101634705A (zh) * 2009-08-19 2010-01-27 西安电子科技大学 基于方向信息测度的sar图像目标变化检测方法
EP2431764A1 (en) * 2010-09-17 2012-03-21 BAE Systems PLC Processing SAR imagery
CN102779346A (zh) * 2012-07-05 2012-11-14 西安电子科技大学 基于改进c-v模型的sar图像变化检测方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于图像分割的SAR图像变化检测算法及实现;黄勇 等;《信号处理》;20050430;第21卷(第2期);第149-152页 *
混合的SAR图像变化检测算法;李玲玲等;《计算机工程与设计》;20150531;第36卷(第5期);第1256-1260页 *

Also Published As

Publication number Publication date
CN105160666A (zh) 2015-12-16

Similar Documents

Publication Publication Date Title
Yuan et al. Factorization-based texture segmentation
Kim et al. Color–texture segmentation using unsupervised graph cuts
Mur et al. Determination of the optimal number of clusters using a spectral clustering optimization
Ghamisi et al. Multilevel image segmentation based on fractional-order Darwinian particle swarm optimization
Huang et al. Classification and extraction of spatial features in urban areas using high-resolution multispectral imagery
CN107633226B (zh) 一种人体动作跟踪特征处理方法
Moallem et al. Optimal threshold computing in automatic image thresholding using adaptive particle swarm optimization
CN109766858A (zh) 结合双边滤波的三维卷积神经网络高光谱影像分类方法
CN109919241B (zh) 基于概率模型和深度学习的高光谱未知类别目标检测方法
CN110458192B (zh) 基于视觉显著性的高光谱遥感图像分类方法及系统
CN108229551B (zh) 一种基于紧凑字典稀疏表示的高光谱遥感图像分类方法
CN105261004A (zh) 基于均值漂移和邻域信息的模糊c均值图像分割方法
CN101763507A (zh) 人脸识别方法及人脸识别系统
CN103985112B (zh) 基于改进多目标粒子群优化聚类的图像分割方法
CN105160666B (zh) 基于非平稳分析与条件随机场的sar图像变化检测方法
Wang et al. Energy based competitive learning
CN110163294B (zh) 基于降维操作和卷积网络的遥感图像变化区域检测方法
US20150356350A1 (en) unsupervised non-parametric multi-component image segmentation method
CN105095913A (zh) 基于近邻正则联合稀疏表示的遥感图像分类方法及系统
CN111460966B (zh) 基于度量学习和近邻增强的高光谱遥感图像分类方法
CN109034213B (zh) 基于相关熵原则的高光谱图像分类方法和系统
CN112084837A (zh) 一种基于深度网络的遥感影像变化检测方法及系统
CN115690086A (zh) 一种基于对象的高分辨率遥感影像变化检测方法及系统
CN112990314A (zh) 基于改进孤立森林算法的高光谱图像异常检测方法及装置
CN105787505A (zh) 一种结合稀疏编码和空间约束的红外图像聚类分割方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20180306

Termination date: 20180825

CF01 Termination of patent right due to non-payment of annual fee