多时相SAR影像变化检测方法

多时相sar影像变化检测方法
技术领域
1.本公开涉及雷达影像处理技术领域，尤其涉及一种多时相sar影像变化检测方法。


背景技术：

2.合成孔径雷达(synthetic aperture radar，sar)影像变化检测是利用覆盖同一地区不同时期的sar影像进行比较分析，来获取地表覆盖变化信息。变化检测的基本流程主要包括三大技术环节，依次为：预处理
→
差异图生成
→
差异图分析(分割)。差异图生成和差异图分析是变化检测过程的重中之重。然而，现有技术中，基本围绕单极化或全极化雷达影像数据开展变化检测，针对双极化雷达影像数据的检测极少。


技术实现要素：

3.有鉴于此，本公开一方面提供一种多时相sar影像变化检测方法，包括：构建第一时相双极化sar影像的第一协方差矩阵与第二时相双极化sar影像的第二协方差矩阵之间的相似矩阵；根据所述相似矩阵计算所述第一时相双极化sar影像与所述第二时相双极化sar影像之间的差异图；对所述差异图进行初始分割，得到所述差异图中各个像素的初始标记，其中，所述初始标记用于区分所述差异图中的变化类及未变化类；根据所述初始标记，基于广义高斯分布和各个像素在八邻域方向上与最邻近像素的上下文关系对所述差异图中的变化类及未变化类的初始标记进行优化，以得到变化二值图。
4.根据本公开的实施例，其中，所述构建第一时相双极化sar影像的协方差矩阵与第二时相双极化sar影像的协方差矩阵之间的相似矩阵包括：根据所述第一时相双极化sar影像对应的共极化通道的后向散射信息和交叉极化通道的后向散射信息构建所述第一协方差矩阵；根据所述第二时相双极化sar影像对应的共极化通道的后向散射信息和交叉极化通道的后向散射信息构建所述第二协方差矩阵；根据所述第一协方差矩阵与所述第二协方差矩阵构建所述相似矩阵。
5.根据本公开的实施例，其中，所述相似矩阵r为：
[0006][0007][0008]
其中，s
xx
为共极化通道的后向散射信息，s
xy
为交叉极化通道的后向散射信息，*表示复数共轭计算；《
·
》表示窗口取平均计算，|
·
|表示复数取模计算，表示矩阵的共轭转置，c1为所述第一协方差矩阵，c2为所述第二协方差矩阵。
[0009]
根据本公开的实施例，其中，所述根据所述相似矩阵计算所述第一时相双极化sar影像与所述第二时相双极化sar影像之间的差异图包括：求解所述相似矩阵的最大特征值
和最小特征值；根据所述最大特征值和所述最小特征值计算所述第一时相双极化sar影像与所述第二时相双极化sar影像之间的相似性测度；根据所述相似性测度得到所述差异图。
[0010]
根据本公开的实施例，其中，所述相似测度sim为：
[0011][0012]
其中，λ
max
、λ
min
分别为所述相似矩阵的最大特征值和最小特征值。
[0013]
根据本公开的实施例，其中，所述对所述差异图进行初始分割，得到所述差异图中各个像素的初始标记包括：基于广义高斯分布模型对所述差异图中的变化类和未变化类进行统计建模，得到所述变化类和未变化类的统计特性；基于所述统计特性，采用kittler-illingworth阈值选取准则确定阈值对所述差异图进行分割，得到各个像素的初始标记。
[0014]
根据本公开的实施例，其中，根据
[0015][0016][0017][0018]
构建所述广义高斯分布模型，其中，x
l
为所述差异图的灰度影像的灰阶，x
l
＝0，1，2，...，l-1，l为所述灰度影像的最大像素值，mi为所述差异图中变化类或未变化类像素的均值，βi为广义高斯分布模型的形状参数，σi为差异图中变化类或未变化类的标准差，γ(
·
)为伽玛函数，ωi为变化类或为变化类，i表示变化类或未变化类的编号，c表示变化类，u表示未变化类；其中，所述kittler-illingworth阈值选取准则以满足代价函数最小为准则，在灰阶范围[0，l-1]确定阈值，根据确定阈值获取所述mi、σi和βi。
[0019]
根据本公开的实施例，其中，根据所述初始标记，基于所述广义高斯分布和各个像素在八邻域方向上与最邻近像素的上下文关系对所述差异图中的变化类及未变化类的初始标记进行优化包括：根据所述初始标记，基于广义高斯分布及各个像素在八邻域方向上与最邻近像素的上下文关系构建所述差异图的马尔科夫随机场模型能量函数；求解最小化的马尔科夫随机场模型能量函数，得到所述差异图中各个像素的最优标记；根据所述最优标记输出变化二值图。
[0020]
根据本公开的实施例，其中，所述马尔科夫随机场模型能量函数为：
[0021][0022]dp
(f
p
)＝ln(g(x
l
|f
p
))，
[0023][0024]
其中，p表示所述差异图中所有像素的集合，p表示像素编号，q表示第p个像素在八邻域方向上最邻近像素的编号，n为第p个像素与最邻近像素组成的集合，f
p
表示第p个像素
的初始标记，fq表示第q个像素的初始标记，φ为平滑参数，λ为平衡因子。
[0025]
根据本公开的实施例，其中，所述求解最小化的马尔科夫随机场模型能量函数包括：随机从所述差异图中选择一像素，在该像素的类别标记为变化类的情况下，依据初始标记获取一未变化类；基于所述马尔科夫随机场模型能量函数计算变化类的初始标记与未变化类的初始标记对应的能量之间的差值；若所述差值满足预设条件，则将当前像素的类别标记由变化类更新为未变化类；以此进行迭代，遍历所述差异图中的其它像素，得到所述差异图对应的能量之差的总和，直至所述能量之差的总和小于预设阈值或迭代次数达到预设次数，获得所述差异图中各像素的最优标记。
附图说明
[0026]
通过以下参照附图对本公开实施例的描述，本公开的上述以及其他目的、特征和优点将更为清楚，在附图中：
[0027]
图1示意性示出了本公开实施例提供的多时相sar影像变化检测方法的流程图。
具体实施方式
[0028]
为使本公开的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本公开进一步详细说明。显然，所描述的实施例是本公开一部分实施例，而不是全部的实施例。基于本公开中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本公开保护的范围。
[0029]
在此使用的术语仅仅是为了描述具体实施例，而并非意在限制本公开。在此使用的术语“包括”、“包含”等表明了所述特征、步骤、操作和/或部件的存在，但是并不排除存在或添加一个或多个其他特征、步骤、操作或部件。
[0030]
在本公开中，除非另有明确的规定和限定，术语“安装”、“相连”“连接”、“固定”等术语应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或成一体；可以是机械连接，也可以是电连接或可以互相通讯；可以是直接连接，也可以通过中间媒介间接相连，可以是两个元件内部的连通或两个元件的相互作用关系。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本公开中的具体含义。
[0031]
在本公开的描述中，需要理解的是，术语“纵向”、“长度”、“周向”、“前”、“后”、“左”、“右”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本公开和简化描述，而不是指示或暗示所指的子系统或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本公开的限制。
[0032]
贯穿附图，相同的元素由相同或相近的附图标记来表示。可能导致本公开的理解造成混淆时，将省略常规结构或构造。并且图中各部件的形状、尺寸、位置关系不反映真实大小、比例和实际位置关系。另外，在权利要求中，不应将位于括号之间的任何参考符号构造成对权利要求的限制。
[0033]
类似地，为了精简本公开并帮助理解各个公开方面中的一个或多个，在上面对本公开示例性实施例的描述中，本公开的各个特征有时被一起分到单个实施例、图或者对其描述中。参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本公开的至少一
个实施例或示例中。本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或者多个实施例或示例中以合适的方式结合。
[0034]
此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。因此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。在本公开的描述中，“多个”的含义是至少两个，例如两个、三个等，除非另有明确具体的限定。
[0035]
本公开的目的在于针对双极化sar数据特性，提出一种基于相似矩阵特征值和马尔科夫随机场模型的多时相双极化sar自动变化检测方法。该方法首先构造前后时相双极化sar数据协方差矩阵之间的相似矩阵，并通过计算相似矩阵的最大、最小特征值来构造一种相似性测度，用来计算前后时相影像之间的差异图。然后，联合基于广义高斯分布的kittler-illingworth自动阈值法和基于马尔科夫随机场模型的影像分割方法，实现差异图中变化类和未变化类的自动精确分割。利用广义高斯分布能够较好地拟合差异图中变化类与未变化类的统计特性，并结合kittler-illingworth阈值选取准则自动确定阈值，实现差异图中变化区域与未变化区域较好的初始标记。但sar采用雷达回波信号相干成像，获得影像中往往存在着较为严重的“颗粒状”斑点噪声，易引起误检。因此，综合利用广义高斯分布和八邻域空间结构来构建差异图的马尔科夫随机场模型能量函数，利用八邻域空间结构能够很好描述差异图中各像素上下文信息的优势，以便有效抑制斑点噪声产生的误检。同时，采用基于metropolis准则的模拟退火全局最优优化算法，不断迭代获得差异图各像素的最优标记，实现变化类和未变化类的精确分割。
[0036]
图1示意性示出了本公开实施例提供的多时相sar影像变化检测方法的流程图。
[0037]
如图1所示，该方法例如可包括操作s101～操作s104。
[0038]
在操作s101，构建第一时相双极化sar影像的第一协方差矩阵与第二时相双极化sar影像的第二协方差矩阵之间的相似矩阵。
[0039]
在本公开实施例中，多时相通常指反映一组遥感影像在时间系列上具有的特征，也即在不同时间获取的同一地域的一组影像。第一时相可以在第二时相之前。
[0040]
在本公开实施例中，可以根据第一时相双极化sar影像对应的共极化通道的后向散射信息和交叉极化通道的后向散射信息构建所述第一协方差矩阵。根据第二时相双极化sar影像对应的共极化通道的后向散射信息和交叉极化通道的后向散射信息构建第二协方差矩阵。其中，后向散射信息指的是雷达向目标(如地面、海洋)发出探测波后对目标反射回来的波进行收集而获取雷达影像中包含的目标的相关信息。
[0041]
示例性的，协方差矩阵c可以定义为：
[0042][0043]
其中，s
xx
为双极化sar中共极化通道的后向散射信息(一般是复数数据)，s
xy
为双极化sar中交叉极化通道的后向散射信息，下标x表示水平极化或垂直极化，下标y表示垂直极化或水平极化，*表示复数共轭计算；《
·
》表示窗口取平均计算，|
·
|表示复数取模计算。
[0044]
再根据第一协方差矩阵与第二协方差矩阵构建相似矩阵。相似矩阵r可以定义为：
[0045][0046]
其中，表示矩阵的共轭转置，c1为第一协方差矩阵，c2为第二协方差矩阵。
[0047]
在操作s102，根据相似矩阵计算第一时相双极化sar影像与第二时相双极化sar影像之间的差异图。
[0048]
在本公开实施例中，可以通过求解相似矩阵的最大特征值和最小特征值。再根据最大特征值和最小特征值计算第一时相双极化sar影像与第二时相双极化sar影像之间的相似性测度。最后根据相似性测度得到差异图。
[0049]
示例性的，相似测度sim可以为：
[0050][0051]
其中，λ
max
、λ
min
分别为所述相似矩阵的最大特征值和最小特征值。
[0052]
在操作s103，对差异图进行初始分割，得到差异图中各个像素的初始标记，其中，初始标记用于区分差异图中的变化类及未变化类。
[0053]
在本公开实施例中，可以利用基于广义高斯分布的kittler-illingworth自动阈值法对差异图进行初始分割。具体可以为：首先，基于广义高斯分布模型对差异图中的变化类ωc和未变化类ωu进行统计建模，得到变化类ωc和未变化类ωu的统计特性。再基于统计特性，采用kittler-illingworth阈值选取准则确定阈值对差异图进行分割，得到各个像素p的初始标记p的初始标记由于广义高斯分布相对于高斯分布仅多出一个形状参数β。通过形状参数β的调整，广义高斯分布能够近似多种分布，如当β＝2时，近似为高斯分布；当β＝1时近似为拉普拉斯分布。因此，利用广义高斯分布能够较好地拟合差异图中变化类与未变化类的统计特性，并结合kittler-illingworth阈值选取准则自动确定阈值，可以实现差异图中变化区域与未变化区域较好的初始标记。
[0054]
广义高斯分布的模型可以定义为：
[0055][0056][0057][0058]
其中，x
l
为差异图的灰度影像的灰阶，x
l
＝0，1，2，...，l-1，l为灰度影像的最大像素值，mi为差异图中变化类或未变化类像素的均值，βi为广义高斯分布模型的形状参数，βi≥0，σi为差异图中变化类或未变化类的标准差，γ(
·
)为伽玛函数，ωi为变化类或为变化类，i表示变化类或未变化类的编号，c表示变化类，u表示未变化类。
[0059]
kittler-illingworth阈值选取准则以满足代价函数j(t)最小为准则，在灰阶范围[0，l-1]确定阈值t，根据确定阈值获取变化类或未变化类对应mi、σi和βi。具体可以为：
[0060][0061][0062][0063]
hc(t)＝1-hu(t)
[0064][0065][0066][0067][0068]
其中，h(
·
)表示差异图的灰阶直方图。
[0069]
在操作s104，根据初始标记，基于广义高斯分布和各个像素在八邻域方向上与最邻近像素的上下文关系对差异图中的变化类及未变化类的初始标记进行优化，以得到变化二值图。
[0070]
在本公开实施例中，可以根据初始标记，基于广义高斯分布及各个像素在八邻域方向上与最邻近像素的上下文关系构建差异图的马尔科夫随机场模型能量函数。求解最小化的马尔科夫随机场模型能量函数，得到差异图中各个像素的最优标记。根据最优标记输出变化二值图。
[0071]
具体的，马尔科夫随机场模型e定义为数据能量ed和平滑能量es的线性组合，数据能量ed利用广义高斯分布来刻画差异图中变化类与未变化类的统计特性，平滑能量es采用八邻域空间结构来描述差异图中各像素p在八邻域方向上与最邻近像素q的上下文依赖关系，即：
[0072][0073]dp
(f
p
)＝ln(g(x
l
|f
p
))，
[0074][0075]
其中，p表示差异图中所有像素的集合，p表示像素编号，q表示第p个像素在八邻域方向上最邻近像素的编号，n为第p个像素与最邻近像素组成的集合，f
p
表示第p个像素的初始标记，fq表示第q个像素的初始标记，φ为平滑参数，一般默认为0.9，λ为数据能量ed和平滑能量es的平衡因子，一般默认为1。
[0076]
可以采用metropolis准则的模拟退火优化算法迭代最小化马尔科夫随机场模型能量函数，获得差异图各像素的最优标记f
′
p
(f
′
p
∈{ωc，ωu})。具体过程可以为：
[0077]
首先，预设模拟退火控制参数，包括初始退火温度t0、迭代终止阈值t和迭代总次数k。
[0078]
接下来，随机从差异图中选择一像素，根据差异图初始标记接下来，随机从差异图中选择一像素，根据差异图初始标记判断当前像素的类别标记，在该像素的类别标记为变化类ωc的情况下，从初始标记组合中获取一未变化类ωu，基于马尔科夫随机场模型能量函数，利用广义高斯分布模型参数mi、σi、βi，以及平滑参数φ，计算前后这两种类别标记对应能量之间的差值δe。
[0079]
接下来，采用metropolis准则判断差值是否满足预设条件，若exp(δe/tk)≥α，α∈(0，1)之间的随机数，则将当前像素的类别标记由变化类ωc更新为未变化类ωu，并计算更新前后能量差值绝对值|δe|，
[0080]
以此进行迭代，采用上述方法遍历差异图中的其它像素，更新各像素类别标记，同时计算整幅差异图满足metropolis准则的总能量∑|δe|。
[0081]
最后，判断迭代是否达到终止条件，包括：总能量∑|δe|是否小于迭代终止阈值t或迭代次数是否达到预设次数k。若当前迭代k不满足终止条件，则降温处理t
k 1
＝λtk，k∈[0，k-1]，λ默认为0.98。重复前述过程，直到满足迭代终止条件，获得差异图各像素的最优标记，实现变化类和未变化类的精确分割，输出变化二值图。
[0082]
综上所述，本公开实施例提供的多时相sar影像变化检测方法，通过构造前后时相双极化sar数据协方差矩阵之间的相似矩阵，并利用相似矩阵的最大、最小特征值来构造一种相似性测度，用来计算前后时相影像之间的差异图，能够有效突出显示前后时相影像间的变化区域；利用“基于广义高斯分布的kittler-illingworth自动阈值初始分割
→
基于马尔科夫随机场模型的精确分割”的层级影像分割策略，利用广义高斯分布能够较好地拟合差异图中变化类与未变化类的统计特性，利用八邻域空间结构能够很好描述差异图中各像素的上下文信息，从而有效抑制斑点噪声产生的误检，实现差异图中变化类和未变化类的自动精确分割，在sar变化检测中起到重要的作用，对于推进双极化sar变化检测实用化具有重要意义。对于推进双极化sar变化检测实用化具有重要意义。
[0083]
以上所述的具体实施例，对本公开的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本公开的具体实施例而已，并不用于限制本公开，凡在本公开的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。