一种多时相SAR图像城市区域变迁变化检测方法

摘要

本发明公开了一种多时相SAR图像城市区域变迁变化检测方法，首先采集同一地域不同时期的两幅多时相SAR图像，进行几何校正和辐射校正；其次对校正后的同一地域不同时期的两幅多时相SAR图像进行处理，得到变化检测结果。本发明检测方法，路线简单、切实可行，易于实现，是专门针对城市遥感图像的应用而提出的。对于呈现复杂多结构化特征的城市区域的空时变迁变化检测，本发明方法具有独特的图像表示优势。

说明书

技术领域

本发明属于电子科学技术领域，具体涉及一种多时相SAR图像城市区域变迁变化检测方法。

背景技术

变化检测是通过分析不同时间获得的图像确定地物分布变化情况的一种过程。在民用和军用领域，变化检测都有着重要的研究意义。例如，通过对农田的监视，可以分析作物的生长态势；通过对城市发展用地的调查，可以评估城市化的发展情况。在军用方面，变化检测技术可以用于战场态势的估计，有助于及时地改变战略部署。

随着SAR成像技术(特别是星载SAR成像技术)和计算机技术的发展，获取SAR图像变得越来越便捷。而且，全天时全天候的工作特性也体现了通过SAR获取多时相遥感图像的独特优势，人们可以在短时间内不受天气状况和光照强度的影响下获得同一地理位置不同时段的SAR图像。但从另一方面来说，SAR图像中固有斑点噪声的存在使得变化检测的难度要高于基于多光谱图像的检测方法。变化检测分为有监督变化检测和非监督变化检测两大类。由于SAR图像的非平稳程度较高，非监督变化检测就成为变化检测的主要方式。近年来，随着SAR成像技术的不断成熟和成像质量的日益提高，基于SAR图像的多时相变化检测研究无论在理论研究上还是在具体应用上都得到了快速的发展。

随着图像信号分解的滤波器设计理论和相关技术的发展，近年来在变化域对SAR图像进行处理的研究逐渐增多，特别是基于小波分析的SAR图像处理。但是，人们在应用小波分析对SAR图像进行处理时，发现小波在表示图像这种二维信号时存在根本性的不足，因此利用这种工具提高SAR图像处理的性能空间比较有限。近年来，许多学者不断探索和研究，提出了超越传统小波分析能力的一系列图像表示工具来弥补小波分析的内在缺陷，这些工具统称为多尺度几何分析(Multiscale Geometrical Analysis，MGA)。

MGA是一门新兴的信号处理技术，它是一批推动小波分析发展的信号处理界的学者不断提出和发展的图像多尺度表示形式，它包括许多种方法：Ridgelet、Brushlet、Curvelet、Contourlet、Bandelet、Beamlet、Directionlet、和Shearlet，等等。这些MGA方法的提出和发展正是致力于探索一种崭新的高维函数最优表示方法，弥补小波分析的不足。这些方法除了具有小波分析的优良特性之外，更重要的是在构造的过程中使用了具有各向异性的基函数，因此能够稀疏地逼近各种奇异曲线或曲面。MGA作为最优表示二维图像信号的有效工具，已经受到了众多学科领域学者的关重，已经开始探索把MGA方法运用到图像处理的不同领域，并取得了很好的效果。

上述提到的各种MGA工具的提出在对图像这种二维信号进行表示上均具有小波分析无可比拟的优点-可以同时对图像进行局部化、多尺度、多方向的表示，但它们之间也存在固有的差别，在不同的图像处理应用中所表现出来的特点互有各异。

在对图像的边缘或脊这些奇异性信息的表示上Shearlet方法在性能上更具有独特的优势，它通过具有高度各向异性的基以最优稀疏表示的形式更细致、更准确地刻画图像中的边缘或脊，对不同粗细程度的边缘或脊描述得更有效，对这些奇异性信息的方向敏感性更强。换句话说，Shearlet方法在对图像中的边缘或脊这类奇异性信息在尺度、位置和方向上能提供更准确更有效的表达效果。图像经过Shearlet变换后每级分解先将图像分解成一个逼近系数子带(主要表现图像的低频信息，包含图像的绝大多数能量)和多个方向子带(主要反映图像在不同方向上的高频信息，包含图像各个方向的细节特征)。而且，由于Shearlet方法是一种对图像进行的多尺度变换，因而也特别适用于从含有噪声的图像中完成结构化信息的提取。

由于Shearlet是一种能够高效表示图像细节信息的MGA方法，在刻画图像的各种结构化特征方面具有独特的优势，因此它特别适合于像城市区域这种具有复杂结构的地物类目标的变化检测。

发明内容

本发明的目的是提供一种多时相SAR图像城市区域变迁变化检测方法，解决了现有检测方法不适合于城市区域这种具有复杂结构的地物类目标的变化检测的问题。

本发明所采用的技术方案是，一种多时相SAR图像城市区域变迁变化检测方法，具体按照以下步骤实施：

步骤1：采集同一地域不同时期的两幅多时相SAR图像，进行几何校正和辐射校正；

步骤2：对步骤1得到的校正后的同一地域不同时期的两幅多时相SAR图像进行处理，得到变化检测结果。

本发明的特点还在于，

其中步骤2对步骤1得到的校正后的同一地域不同时期的两幅多时相SAR图像进行处理，具体按照以下步骤实施：

步骤a：多时相图像的Shearlet变换，

采用Laplacian金字塔分解和带通滤波器滤波，采用3-4级分解，每一级分解为4-8个子带方向，变换后得到图像在不同尺度、不同方向子带内的Shearlet系数图像；

步骤b：子带差别图的确立，

将两幅SAR图像的Shearlet系数在各个方向子带对应地相减，然后取绝对值，具体按照如下表达式来执行：

${\mathrm{DM}}_{l,d}(i,j)=|S_{l,d}^{I_1}(i,j)-S_{l,d}^{I_2}(i,j)|,$

其中，为SAR图像I_k，k＝1，2在尺度l方向子带d中对应位置为(i，j)的Shearlet系数，DM_l，d(i，j)为尺度l方向子带d中对应位置为(i，j)的差别值；

步骤c：通过统计分类获取子带变化图，

采用Bayesian分类来完成，将DM_l，d(i，j)对应的变化类属记为C_l，d(i，j)，按照模式分类中最大后验分类的准则来确定变化的类属情况：变化类“w_changed”或无变化类“w_unchanged”，其依据的表达式为：

$C_{l,d}(i,j)=\underset{w_m\in \{w_{\mathrm{changed}},w_{\mathrm{unchanged}}\}}{\arg \max }P\left(w_m\right|{\mathrm{DM}}_{l,d}(i,j)),$

在上述表达式中，

$P\left(w_m\right|{\mathrm{DM}}_{l,d}(i,j))=\frac{P_{l,d}\left(w_m\right)p\left({\mathrm{DM}}_{l,d}(i,j)\right|w_m)}{\underset{w_m}{\Sigma }{P}_{l,d}\left(w_m\right)p\left({\mathrm{DM}}_{l,d}(i,j)\right|w_m)},w_m\in \{w_{\mathrm{changed}},w_{\mathrm{unchanged}}\},$

其中，P_l，d(w_m)为类属w_m的先验概率且P_l，d(w_changed)+P_l，d(w_unchanged)＝1，p(DM_l，d(i，j)|w_m)为在类属w_m给定的条件下差别值DM_l，d(i，j)对应的似然性值；

确定“w_changed”和“w_unchanged”类的先验概率P_l，d(w_m)以及计算p(DM_l，d(i，j)|w_m)：假设p(DM_l，d(i，j)|w_m)近似服从高斯分布，子带差别图的直方图呈现一种“尖峰值、长拖尾”的非高斯分布形式，将这种分布形式用高斯混合模型来表达：

p(DM_l，d(i，j))＝P_l，d(w_changed)p(DM_l，d(i，j)|w_changed)+P_l，d(w_unchanged)p(DM_l，d(i，j)|w_unchanged)，w_changed类和w_unchanged类所对应的p(DM_l，d(i，j)|w_m)通过各自对应的均值和方差来计算；

P_l，d(w_m)和p(DM_l，d(i，j)|w_m)对应的均值与方差参数的估计：采用期望最大化算法；

接着，按照如下规则确立尺度l方向子带d中对应位置为(i，j)的变化类属值，从而获得多尺度子带变化图：

${\mathrm{CM}}_{l,d}(i,j)=\left(\begin{array}{cc}1,& \mathrm{if}{C}_{l,d}(i,j)=w_{\mathrm{changed}}\\ 0,& \mathrm{if}{C}_{l,d}(i,j)=w_{\mathrm{unchanged}}\end{array}\right),$

其中，CM_l，d(i，j)为尺度l方向子带d对应的变化图CM_l，d中位置(i，j)的值，C_l，d(i，j)为尺度l方向子带d位置(i，j)的类属；

步骤d：尺度内不同子带变化图的融合，

对于各个尺度内子带间变化图的融合，直接将尺度l∈{1，ΛL}各个方向对应位置的变化类属值进行“OR”操作，表达式如下：

${\mathrm{CM}}_l(i,j)=\underset{\forall d}{Y}{\mathrm{CM}}_{l,d}(i,j),$

步骤e：尺度间变化图的融合，

对于尺度间变化图的融合，首先将不同尺度获得变化图的尺寸进行扩大，使得这些尺度的变化图的尺寸与原始图像的尺寸一致，采用非线性插值，然后，采用如下的“AND”操作来完成尺度间的变化图融合，从而获得最终的变化检测结果：

$\mathrm{CM}(i,j)=\underset{l\in \{1,...,L\}}{I}{\mathrm{CM}}_l(i,j),$

其中，CM(i，j)为1表示变化图CM中位置(i，j)处发生了变化，为0则没有发生变化。

本发明的有益效果是，该方法路线简单、切实可行，易于实现，是专门针对城市遥感图像的应用而提出的。对于呈现复杂多结构化特征的城市区域的空时变迁变化检测，Shearlet方法具有独特的图像表示优势。

附图说明

图1是本发明多时相SAR图像城市区域变迁变化检测方法的流程图；

图2是某一地区1996年获取的ERS-2SAR图像；

图3是图2中某一地区2002年获取的ERS-2SAR图像；

图4是图2中某一地区区域变迁变化的检测结果。

具体实施方式

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

本发明多时相SAR图像城市区域变迁变化检测方法，采用“先比较-后分类”的处理思路来完成多时相SAR图像的变化检测任务，通过使用多尺度Shearlet变换分别对不同时段获取的SAR图像进行多尺度分解，然后通过取差操作，获取对应子带的差别图，并按照Bayesian分类的思想对各个子带的差别图进行统计分类，获取多尺度变化检测图，接着通过尺度内融合的思想对各个尺度内所有子带的变化图进行融合，最后实施尺度间的变化图融合，从而实现整个变化检测任务。如图1所示，具体按照以下步骤实施：

步骤1：采集同一地域不同时期的两幅多时相SAR图像，通过专业软件进行几何校正和辐射校正；

步骤2：对步骤1得到的校正后的同一地域不同时期的两幅多时相SAR图像进行处理，得到变化检测图。具体按照以下步骤实施：

步骤a：多时相图像的Shearlet变换

两幅多时相图像的Shearlet变换所采用的分解方式是一致的，包括多尺度分解和方向性分解所采用的滤波器组、分解的级数、分解的子带方向数均是相同的。一般情况下，需要采用3-4级分解，每一级分解为4-8个子带方向。Shearlet变换主要需要Laplacian金字塔分解和带通滤波器滤波两个关键步骤来完成，变换后会得到图像在不同尺度、不同方向子带内的Shearlet系数图像。

步骤b：子带差别图的确立

直接将两幅SAR图像的Shearlet系数在各个方向子带对应地相减，然后取绝对值。具体地说，可以按照如下表达式来执行：

${\mathrm{DM}}_{l,d}(i,j)=|S_{l,d}^{I_1}(i,j)-S_{l,d}^{I_2}(i,j)|,---\left(1\right)$

其中，表示SAR图像I_k，k＝1，2在尺度l方向子带d中对应位置为(i，j)的Shearlet系数，DM_l，d(i，j)则为尺度l方向子带d中对应位置为(i，j)的差别值。

步骤c：通过统计分类获取子带变化图

对于统计分类，采用Bayesian分类的思想来完成，将DM_l，d(i，j)对应的变化类属记为C_l，d(i，j)，则可以按照模式分类中最大后验分类的准则来确定变化的类属情况：变化类“w_changed”或无变化类“w_unchanged”，其依据的表达式为：

$C_{l,d}(i,j)=\underset{w_m\in \{w_{\mathrm{changed}},w_{\mathrm{unchanged}}\}}{\arg \max }P\left(w_m\right|{\mathrm{DM}}_{l,d}(i,j))---\left(2\right)$

在上述表达式中，

$P\left(w_m\right|{\mathrm{DM}}_{l,d}(i,j))=\frac{P_{l,d}\left(w_m\right)p\left({\mathrm{DM}}_{l,d}(i,j)\right|w_m)}{\underset{w_m}{\Sigma }{P}_{l,d}\left(w_m\right)p\left({\mathrm{DM}}_{l,d}(i,j)\right|w_m)},w_m\in \{w_{\mathrm{changed}},w_{\mathrm{unchanged}}\}.---\left(3\right)$

其中，P_l，d(w_m)为类属w_m的先验概率且P_l，d(w_changed)+P_l，d(w_unchanged)＝1，而p(DM_l，d(i，j)|w_m)为在类属w_m给定的条件下差别值DM_l，d(i，j)对应的似然性值。因此，问题的关键就在于如何确定“w_changed”和“w_unchanged”类的先验概率P_l，d(w_m)以及计算p(DM_l，d(i，j)|w_m)。

为解决上述问题，采用的策略是假设p(DM_l，d(i，j)|w_m)近似服从高斯分布，这与Shearlet域中子带差别图的先验分布是保持一致的：前期的研究发现，子带差别图的直方图基本上呈现一种“尖峰值、长拖尾”的非高斯分布形式。这样，就可以将这种分布形式用高斯混合模型来表达：

p(DM_l，d(i，j))＝P_l，d(w_changed)p(DM_l，d(i，j)|w_changed)+P_l，d(w_unchanged)p(DM_l，d(i，j)|w_unchanged)    (4)

这样，w_changed类和w_unchanged类所对应的p(DM_l，d(i，j)|w_m)可以通过各自对应的均值和方差来计算。于是，我们要解决的问题就转化为P_l，d(w_m)和p(DM_l，d(i，j)|w_m)对应的均值与方差参数的估计。对于这类参数估计问题，我们就可以采用统计分析中经典的期望最大化(EM)算法直接来解决。这种方法在基于空域方法的多时相遥感图像变化检测过程已经被采用，被证实是可行有效的。

接着，我们就可以按照如下规则确立尺度l方向子带d中对应位置为(i，j)的变化类属值，从而可以获得多尺度子带变化图。

${\mathrm{CM}}_{l,d}(i,j)=\left(\begin{array}{cc}1,& \mathrm{if}{C}_{l,d}(i,j)=w_{\mathrm{changed}}\\ 0,& \mathrm{if}{C}_{l,d}(i,j)=w_{\mathrm{unchanged}}\end{array}\right),---\left(5\right)$

其中，CM_l，d(i，j)为尺度l方向子带d对应的变化图CM_l，d中位置(i，j)的值，C_l，d(i，j)为尺度l方向子带d位置(i，j)的类属。

步骤d：尺度内不同子带变化图的融合

对于各个尺度内子带间变化图的融合，我们的思路是直接将尺度l∈{1，ΛL}各个方向对应位置的变化类属值进行“OR”操作，表达式如下：

${\mathrm{CM}}_l(i,j)=\underset{\forall d}{Y}{\mathrm{CM}}_{l,d}(i,j)---\left(6\right)$

步骤e：尺度间变化图的融合

对于尺度间变化图的融合，由于Shearlet是一种下采样形式的分解，因此首先将不同尺度获得变化图的尺寸进行扩大，使得这些尺度的变化图的尺寸与原始图像的尺寸一致。对于这个问题，采用最近邻插值或其它非线性插值来解决。当然，也可以进一步考虑设计一种更灵活有效的自适应插值方法。然后，我们将采用如下的“AND”操作来完成尺度间的变化图融合，从而获得最终的变化检测结果：

$\mathrm{CM}(i,j)=\underset{l\in \{1,...,L\}}{I}{\mathrm{CM}}_l(i,j)---\left(7\right)$

CM(i，j)为1表示变化图CM中位置(i，j)处发生了变化，为0则没有发生变化。

本发明多时相SAR图像城市区域变迁变化检测方法，通过使用3级离散Shearlet变换，每一级分解分解的方向数目为4，4，8(从粗到细尺度)，然后按照上述步骤依次完成子带差别图、子带变化图的实现和尺度内以及尺度间变化图的融合过程。图2给出了某一地区1996年获取的ERS-2SAR图像，图3则是在该地区2002年获取的ERS-2SAR图像，图4是通过本发明提出的方法实现的该地区区域变迁变化的检测结果。从该图可以看出，该方法对于发生变化的区域(不论大小和形状)均具有较好的检测精度，具有低的漏警率和合理的虚警率。