CN101634709B - 基于多尺度积和主成分分析的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于多尺度积和主成分分析的SAR图像变化检测方法，主要解决现有技术自适应性差、适用范围窄以及变化检测结果受图像误配准影响的问题。具体实现过程是：首先对输入的两时相SAR图像进行对数比运算得到差异图，对其进行小波变换，对每一分解层高频信息都进行多尺度积去噪；然后将每一层去噪后的图像组合起来进行主成分分析变换，并将第一主成份图作为新的差异图；最后采用广义高斯模型的最小错误率阈值对新的差异图分类，得到最后的变化结果图。实验表明本发明能够增强变化信息，抗噪声性强，能够减弱图像误配准的影响，具有很好的适用性，可用于SAR图像的灾情检测。

Description

基于多尺度积和主成分分析的SAR图像变化检测方法

技术领域

本发明属于数字图像处理技术领域，具体的说是一种基于多尺度积和主成分分析的SAR图像变化检测的方法，适用于遥感图像处理与分析。

背景技术

遥感图像的变化检测是指对同一地理位置不同时期的遥感图像进行分析获得其中的变化信息。由于合成孔径雷达(Synthetic Aperture Radar，SAR)具有全天候、全天时获取地表信息的能力，并具有一定的穿透能力，所以SAR图像的解译和应用已成为现代遥感技术研究的一个热点问题。SAR图像变化检测技术对地震区域的定位和灾害评估等民用领域，以及战场区域或重点监视目标实现连续侦察，可用于目标打击效果评估、战场信息动态感知、军事目标和***监测等军事方面得到了广泛应用。

目前，基于非监督的SAR图像变化检测方法首先对两时相SAR图像分别做对数运算，将两时相SAR图像中的斑性噪声转化为加性噪声，再将两时相图做比值较得到两时相图的差异图，然后对该差异图进行有效的分析处理得到图像的“变化”和“非变化”区域。这类方法由于对图像几何配准误差以及斑点噪声比较敏感，因而影响了提取变化信息的质量。

近几年，为了降低SAR图像中斑噪的影响，有关学者提出以下几种方法：

一、Francesca Bovolo等学者在文章“A detail-preserving scale-driven approach tochange detection in multitemporal SAR images”提出了基于小波变换的变化检测方法，该方法对差异图进行小波变换，并对小波变换的每一尺度的低频信息进行人工阈值初始分割，然后将初始分割的结果进行决策得到最后的变化检测结果。该方法优点是抗噪声性强，但该方法缺点是最后的变化检测结果依赖于初始分割的效果，初始分割效果不好，则最后的变化检测结果也会受到影响，而且该方法阈值的选取需要人工干预。

二、Luis M.T.de Carvalho等学者在文章“Wavechange：a Procedure for ChangeDetection Based on Wavelet Product Spaces”提出了小波多尺度积的变化检测方法，该方法对差异图采用小波变换，并将某分解层与其相邻层进行相乘，得到该层的多尺度积来增强变化信息，去除噪声，并利用某分解层的低频信息对变化区域定位，并根据变化区域位置计算局部最大，去除一些伪变化信息，并将得到的局部最大点作为种子点进行区域生长得到很好的变化结果，证明了计算多尺度积能够去除由于几何辐射错配准引起的误差。但是该方法的缺点是最后的变化检测结果受某一分解层低频中对变化信息的定位以及对局部最大值的定位准确性的影响，从而限制了该方法的适用性。

三、张辉等学者提出了基于主成分分析的SAR图像变化检测方法，其中将两时相图像进行主成分分析分解，通过数学推导和理论分析，发现其主分量主要表征了两图的未变化部分，而次分量则可用于表征图像的变化部分，此方法虽说抗斑噪性强，但是结果易受图像误配准的影响，使其检测的变化结果中存在大量的伪变化信息。

发明内容

本发明的目的在于克服上述已有的SAR图像变化检测技术的不足，提出了一种基于多尺度积和主成分分析的SAR图像变化检测方法，以减少两时相SAR图像误配准的影响，并采用自适应阈值分割提高自适应性，并通过考虑不同尺度的信息提高其适用性。

为实现上述目的，本发明的检测方法包括如下步骤：

(1)对输入的两时相SAR图像进行对数比运算，得到一幅对数比差异图X_LR；

(2)对对数比差异图X_LR进行二维平稳小波变换，分解为4层，对每一层的水平、对角、垂直方向的高频信息分别进行多尺度积去噪，并将每一层去噪后的高频系数和该尺度的低频系数进行逆小波变换，得到每个层去噪之后的对数比差异图X_LR ^j(j＝1，2，...4)；

(3)将每层去噪之后的对数比差异图X_LR ^j分别列向量化并将四层列向量组成矩阵，对该矩阵进行主成分分析变换，将得到的第一主成分图作为新的差异图X’_LR，本发明称为第一主成分差异图；

(4)对第一主成分差异图X’_LR利用广义高斯模型最小错误准则计算其阈值并进行分类，大于阈值的像素点属于变化类，否则属于非变化类，得到变化结果图。

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

(1)本发明由于考虑了小波变换分解的每个分解层的高频、低频信息，使得变化检测结果不受某一分解层的影响，因此具有很好的适用性。

(2)本发明由于对每个分解层的高频信息采用多尺度积去噪，很好地抑制了斑点噪声。

(3)本发明由于采用主成分分析融合了差异图的每个分解层去噪之后的图的有用信息，去除了分解层图间的冗余信息，减少了图像误配准的影响，增强了变化信息并达到了降维的效果。

(4)本发明由于采用广义高斯模型能够更准确地拟合差异图直方图，并根据贝叶斯最小错误准则求得差异图阈值，该方法具有自适应性，并且根据该方法求得的阈值对差异图分类能够得到较精确的变化检测结果。

附图说明

图1是本发明的实现流程示意图；

图2是本发明输入的前后两个时相的SAR图像；

图3是本发明的对数比差异图；

图4是本发明主成分分析变换得到的第一主成分差异图；

图5是本发明实验中的变化检测结果参考图；

图6是本发明与现有的三种方法得到的变化检测结果图。

具体实施方式

参照图1，本发明的实施如下：

步骤1，对输入的同一地理位置的两个不同时相的大小为m*n的SAR图像进行对数比运算，得到一幅差异图X_LR，如图3所示。

步骤2，对对数比差异图X_LR进行二维平稳小波变换，分解为4层，对每一层的水平、对角、垂直方向的高频信息分别进行多尺度积去噪，并将每一层去噪后的高频系数和该尺度的低频系数进行逆小波变换，得到每个分解层去噪之后的对数比差异图，具体步骤如下：

(2a)对对数比差异图X_LR进行二维平稳小波(2D-SWT)变换。本发明设置小波变换的分解层数为四层，其中每一分解层图中包括高频系数图和低频系数图，高频系数图包括水平、垂直、对角方向的高频系数图；

(2b)对每一分解层的水平、垂直、对角方向的高频系数进行多尺度积去噪，令W(m，n)为小波变换系数，对每个分解层的高频系数计算多尺度积，并根据空间选择性的噪声滤波方法对每个分解层的高频系数去噪，其具体步骤如下：

2b1)计算当前分解层的多尺度积，并将其定义为相邻层间小波系数W_j ^d(m，n)的相关性C_j ^d(m，n)，即

${C^d}_j(m,n)=\underset{j}{\overset{j+1}{\mathrm{\Π}}}{W^d}_j(m,n),j=\mathrm{1,2},...4---\left(1\right)$

其中j为小波分解当前层，d为水平、垂直、对角方向；

2b2)计算小波相关性系数C_j ^d(m，n)和小波系数W_j ^d(m，n)的能量，即

${\mathrm{PW}}_j^d\left(m\right)=\underset{m,n}{\mathrm{\Σ}}{\left(W_j^d\right)}^2(m,n)---\left(2\right)$

${\mathrm{PC}}_j^d=\underset{m,n}{\mathrm{\Σ}}{\left(C_j^d\right)}^2(m,n);---\left(3\right)$

2b3)对小波系数相关性系数进行归一化：

${\mathrm{NC}}_j^d(m,n)=C_j^d(m,n)*\left(\sqrt{{\mathrm{PW}}_j^d/{\mathrm{PC}}_j^d}\right);---\left(4\right)$

2b4)比较归一化后的相关性系数和小波能量系数，如果 $\left|{\mathrm{NC}}_j^d(m,n)\right|>\left|W_j^d(m,n)\right|,$ 则认为该像素是有用信息，将该像素对应的小波系数保留，否则认为是噪声，并将该像素的小波系数置零，即

$W_j^d(m,n)=\left\{\begin{array}{cc}{W}_j^d(m,n),& \mathrm{if}\left|{\mathrm{NC}}_j^d(m,n)\right|>\left|W_j^d(m,n)\right|\\ 0,& \mathrm{else}\end{array}\right.;---\left(5\right)$

(2c)将处理之后的每个分解层的小波高频系数和该层的低频系数进行二维平稳小波逆变换，得到差异图每个分解层去噪之后的图像，记为X_LR ^j。

步骤3，将每层去噪之后的对数比差异图□分别列向量化并将其四层列向量组成矩阵A，对该矩阵进行主成分分析变换，并将得到的第一主成分图Z₁作为新的差异图，本发明称为第一主成分差异图，其具体步骤如下：

(3a)将每个分解层去噪之后的图像X_LR ^j列向量化，表示为 $A_j(N,1)=X_{\mathrm{LR}}^j(m,n),$ 其中j＝1，2，...，4，N＝m*n，并组成矩阵A＝(A₁，A₂，A₃，A₄)；

(3b)计算矩阵A的相关系数矩阵R

$R=\left[\begin{array}{cccc}{r}_{11}& r_{12}& r_{13}& r_{14}\\ r_{21}& r_{22}& r_{23}& r_{24}\\ r_{31}& r_{32}& r_{33}& r_{34}\\ r_{41}& r_{42}& r_{43}& r_{44}\end{array}\right]---\left(6\right)$

其中矩阵元素r_ij

$r_{\mathrm{ij}}=\frac{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}(A_k-\stackrel{\‾}{A_i})(A_k-\stackrel{\‾}{A_j})}{\sqrt{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{(A_k-\stackrel{\‾}{A_i})}^2\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{(A_k-\stackrel{\‾}{A_j})}^2}}---\left(7\right)$

其中A_i为矩阵A的第i各向量的均值，A_k为矩阵A中第k个向量，n为矩阵A中总的向量数，n＝4；

(3c)根据相关系数矩阵R得到特征方程为：

|λI-R|＝0；(8)

(3d)采用雅可比法求特征方程的特征值λ_i(i＝1.....4)并按从大到小排列，计算特征值λ_i(i＝1.....4)对应的特征向量u_i(i＝1.....4)；

(3e)根据特征值λ_i(i＝1.....4)计算累计贡献率s，即

$s=\frac{{\mathrm{\Σ}}_{i=1}^M{\mathrm{\λ}}_i}{{\mathrm{\Σ}}_{j=1}^p{\mathrm{\λ}}_j},(M=1....4)---\left(9\right)$

其中M为待确定的主成分个数，p为总的特征值的个数，且M≤p，i＝1，2，...M，j＝1，2，...，p；

(3f)确定主成分个数M，一般取使累计贡献率s达到85％-95％的特征值λ₁，λ₂.....λ_M所对应的第一、第二、……、第M个主成分；

(3g)列向量化后的图像序列通过主成分分析逆变换，得到的主成分图像为：

Z_i＝Au_i(i＝1，2，...M)(10)

并将第一主成分图Z₁作为新的差异图X’_LR第一主成分差异图如图4所示。

步骤4，对第一主成分差异图X’_LR采用广义高斯模型拟合其直方图，估计X’_LR中变化类和非变化类的概率密度函数，并根据最小错误率准则计算新差异图X’_LR的阈值，然后根据得到的阈值对新的差异图X_LR’分类，图X_LR’中大于阈值的像素点属于变化类，否则属于非变化类，得到变化检测结果图，图6(a)所示。

本发明的效果可以通过以下内容进行说明：

1.实验数据

为了验证本发明所提出SAR变化检测方法的性能，对一组仿真SAR图像数据集和两组真实SAR图像数据集总共三组数据集进行实验，每一组数据集都带有检测结果参考图。第一组的仿真实验数据集是图像大小为256*256的两时相SAR图像，其中第一时相SAR图像仿真过程：首先按照吉布斯随机场分布形式获得具有两类场景的参考图像，然后按照SAR图像的成像机理获得带有相干斑的SAR图像，将仿真的变化区域嵌入到参考图像中，第二时相SAR图像按照第一时相SAR图像类似的形成过程获得。第二组真实实验数据集是图像的大小为704×704的真实的两时相SAR图像，第一时相SAR图像如图2(a)所示，第二时相SAR图像如图2(b)所示，这两幅图是通过ERS-2获得的意大利帕维亚地区两幅SAR图像，第一时相SAR图像是在2000年10月20日获得，该图是在洪水灾害刚刚发生后获得的，图像中阴暗部分为首洪水影响的区域。第二时相SAR图像是在2000年10月28日获得，此时洪水几乎完全消失。第三组真实实验数据集是图像大小为288×288的通过ERS-2获得的Bern城市的两时相SAR图像，第一时相SAR图像是1999年4月获得的，第二组实验数据集是1999年5月获得的，在两时相之间，该地区发生了水灾，第二时相图像中阴暗区域即为洪水影响的区域。

2.对比实验

为了说明本发明采用第一主成分作为差异图对变化区和非变化区的增强效果，将采用第一主成分的差异图与对数比差异图进行了对比。

本发明与现有的三个变化检测方法进行对比，对比实验一是Francesca Bovolo等学者在文章“A detail-preserving scale-driven approach to change detection in multitemporalSAR images”提出的基于小波变换的变化检测方法，其中只考虑小波变换的每个分解层的低频信息，并采用主成分分析将每一层的低频信息融合起来，基于广义高斯模型最小错误阈值对主成分分析分解后的第一主成分图进行分割得到最后的变化检测结果图，该方法的变化结果图见图6(b)。对比实验二是Luis M.T.de Carvalho等学者在文章“Wavechange：a Procedure for Change Detection Based on Wavelet Product Spaces”提出了小波多尺度积的变化检测方法，其中选择对数比差异图小波变换的第三分解层进行尺度积去噪，然后用并基于广义高斯模型最小错误阈值对差异图去噪后的图进行分割，该方法的变化检测结果见图6(c)。对比实验三是张辉等学者在文章“基于主成分分析的SAR图像变化检测方法”提出的基于主成分分析的SAR图像变化检测方法，其中将通过主成分分析变换后的次主成分图作为差异图，并基于广义高斯模型最小错误阈值对差异图图进行分割，该方法的变化检测结果见图6(d)。两时相SAR图像的变化检测结果参考图如图5所示。

通过与两时相SAR图像的变化检测结果参考图进行比较，定量评价了本发明与两个对比实验的变化检测方法的性能，评价指标包括虚警数，漏检数和总错误数，其中虚警数是变化检测结果图如图6(a)、图6(b)、图6(c)和图6(d)中为变化类而参考图中为非变化类的个数，漏检数是参考图中为变化类而变化检测结果图如图6(a)、图6(b)、图6(c)和图6(d)中为非变化类的个数，总错误数为虚警数和漏检数之和。表格中给出了本发明与两个对比实验的变化检测方法的性能评定结果。

3.实验结果与分析

首先是采用第一主成分差异图与不采用主成分分析变换的差异图比较，第一主成分差异图如图4所示，两时相SAR图像的对数比差异图如图3所示，可以看出图4比图3减弱了伪变化信息，同时也增强了变化信息。

本发明与上文提及的三个方法对三组实验数据集的变化检测结果进行了比较，结果分别如图6(a)、图6(b)、图6(c)和图6(d)所示。与两时相SAR图像的变化检测结果参考图相比如图5所示，对比实验三的变化检测结果图，如图6(d)所示，其中存在很多由于配准误差导致的伪变化信息以及由于斑性噪声引起的杂点，对比实验一变化检测结果图，如图6(b)所示，其中虽然伪变化信息较少，但是仍存在一些杂点。对比实验二的变化检测结果图，如图6(c)所示，伪变化信息很少，但是变化部分细节信息保持不好。而本发明的变化检测结果图，如图6(a)所示，其中伪变化信息和杂点都得到大大减少，细节信息保持较好，而且可以看出本发明的变化检测结果图最接近于参考图如图5所示，且本发明的变化检测结果图中变化部分的边缘保持较好，证明了本发明抗噪声性强，同时能够降低图像误配准的影响。从表1，表2，表3中可以看出，与其他两个对比实验相比，本发明的得到变化检测结果中总错误数都大大减少，远远好于其他两种方法的变化检测结果，证实了本发明算法稳定不受噪声和图像误配准的影响，得到更加准确的变化检测结果，而且在实际SAR图像中也能够得到准确的变化检测结果，证明了本发明具有很好的适用性。

表1仿真SAR图像变化检测结果性能评价

表2第一组真实SAR图像变化检测结果性能评价

表3第二组真实SAR图像变化检测结果性能评价

Claims (2)

1.一种基于多尺度积和主成分分析的SAR图像变化检测方法，包括如下步骤：

1)对输入的两时相SAR图像进行对数比运算，得到一幅对数比值差异图X_LR；

2)对差异图X_LR进行二维平稳小波变换，分解为4层，对每一层的水平、对角、垂直方向高频信息分别进行多尺度积去噪，并将每一层去噪后的高频系数和该尺度的低频系数进行逆小波变换，得到对数比值差异图X_LR的每个层去噪之后的图

j＝1，2，...4；

3)将每层去噪之后的对数比差异图

分别列向量化并将四层列向量组成矩阵，对该矩阵进行主成分分析变换，将得到的第一主成分图作为新的差异图

，称为第一主成分差异图；

4)对第一主成分差异图

利用广义高斯模型最小错误准则计算其阈值并进行分类，大于阈值的像素点属于变化类，否则属于非变化类，得到变化结果图。

2.根据权利要求1所述的SAR图像变化检测方法，其中步骤2)所述的对每一层的水平、对角、垂直方向高频信息分别进行多尺度积去噪，采用空间选择性的噪声滤波方法去噪。

Images