CN104778719A - 基于二阶统计量扰动分析的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于时变SAR图像二阶统计量扰动分析的变化检测方法，主要解决时变相干斑噪声影响SAR图像变化检测结果的问题。其实现过程是：(1)将已配准且进行辐射校正和几何校正的两幅待检测SAR图像做差，并取模产生差图像；(2)用两幅待检测图像和差图像构造二阶统计量向量；(3)结合扰动分析计算二阶统计量向量与不变区域二阶统计量的标准形式的相关系数，用相关系数构造差异图；(4)用模糊C均值聚类方法对差异图聚类，实现变化检测。本发明有效地减少了由于时变相干斑噪声造成的虚警，可用于对不同时段同一场景变化区域的检测。

Description

基于二阶统计量扰动分析的SAR图像变化检测方法

技术领域

本发明属于雷达技术领域，涉及SAR图像检测方法，可用于对不同时段同一场景的变化区域的检测。

背景技术

雷达成像技术是20世纪50年代发展起来的，在此后的60年里得到了突飞猛进的发展，目前，已经在军事、农林、地质、海洋、灾害、测绘等诸多方面得到广泛的应用。

SAR作为一种主动微波传感器，具有全天候、全天时、分辨率高以及穿透能力强等特点，在突发事件、自然灾害的检测和评估等方面具有独特优势，利用SAR图像进行变化检测具有重要意义。目前，SAR图像变化检测已经在多个方面取得广泛的应用，如土地分析、森林采伐监测、灾情估计、军事侦察、打击效果评估等。

SAR图像变化检测分为监督和非监督两大类。其中，非监督变化检测不需要先验变化信息的支持，可直接从配准好的前后两个时刻SAR图像获取变化信息，这种检测方法不仅降低人为误差的影响，同时也符合实际中先验变化信息获取困难的现实情况，因此非监督变化检测成为国内外研究的重点。SAR图像非监督变化检测的流程可概括为三个步骤，首先对待检测的两幅图像进行图像配准、斑噪滤波、辐射校正、几何校正等预处理，然后比较两幅图像获得差异图，最后对差异图进行分割得到最终的检测结果。

现有文献提出很多差异图的的获取方法。例如将两幅图像分别按像素取对数，然后将取对数后的结果做差后取模值的结果作为差异图；或将两幅图像按像素互相做商，取商较小的值，与一相减并取模值作为差异图中的该像素的值。对图像取对数的目的是将SAR图像的相干斑噪声由乘性噪声转化为加性噪声，并通过对两幅取对数之后的SAR图像做差来消除相干斑噪声。但实际上，由于相干斑噪声有时变性，SAR图像的相干斑噪声由乘性转化为加性之后直接相减并不能地抑制时变的相干斑噪声。将两幅图像按像素做商的目的也是为了消除乘性的相干斑噪声，同样由于相干斑噪声具有时变特性，按像素做商也不能完全抑制相干斑噪声。由于差异图中相干斑噪声的存在会影响后续的差异图分割的性能，造成变化检测的结果中存在大量虚警。

由于两次成像时飞机飞行的方向存在微小差异，不变区域中同一位置的雷达回波的强度会存在轻微的扰动，这种扰动会导致该不变区域在差异图中呈现较大的幅值，影响差异图分割的性能，导致变化检测的结果中存在虚警。

发明内容

本发明针对上述已有技术的不足，提出一种基于二阶统计量扰动分析的SAR图像变化检测方法，以减小相干斑噪声和不变区域的扰动对差异图分割的影响，提高变化检测的性能。

本发明是这样实现的：

一、技术思路

变化检测的关键在于找到一个特征域，在这个特征域中，变化区域与非变化区域有良好的可分性。本发明结合扰动分析，以SAR图像二阶统计量的相关系数域作为特征域得到差异图。由于SAR图像的相干斑噪声有时变性，相干斑噪声对每个像素的幅值造成的变化差异很大，但是相干斑噪声对每个像素造成的功率差异是一定的，本发明结合SAR图像中像素的邻域信息来计算SAR图像的二阶统计量，来描述变化区域与非变化区域的功率信息，以有效地抑制SAR图像的相干斑噪声。在两幅SAR图像中，未变化区域本身可能存在扰动，导致在特征域中这些有扰动的区域比较接近变化区域，本发明通过引入扰动分析来降低未变化区域中扰动的影响，并利用SAR图像二阶的统计量来抑制SAR图像的相干斑噪声，结合扰动分析，获取到有更准确的变化区域以及良好的可分性的差异图。用该差异图进行变化检测以有效减少非变化区域的虚警，提高检测率。

二.技术方案

本发明基于二阶统计量扰动分析的SAR图像变化检测方法，包括如下步骤：

(1)对采集场景初次的SAR图像X和该场景出现变化后的SAR图像Y依次进行配准、几何校正、辐射校正的预处理，并对这两幅图像做差取模值，得到差图像Z＝|X-Y|；

(2)用上述三个图像X、Y、Z中相同坐标的元素构造一组输入向量k＝[X_(x,y),Y_(x,y),Z_(x,y)]^T，并将k的形式简化为k＝[k₁,k₂,k₃]^T，其中X_(x,y)是图像X在(x，y)坐标处的像素值，Y_(x,y),Z_(x,y)分别对应图像Y,Z在(x，y)坐标处的像素值，k₁,k₂,k₃分别与X_(x,y),Y_(x,y),Z_(x,y)一一对应；

(3)根据向量k定义二阶统计量向量t：

t＝[<k₁ ²>,<k₂ ²>,<k₃ ²>,<k₁k₂>,<k₁k₃>,<k₂k₃>]^T/||[<k₁ ²>,<k₂ ²>,<k₃ ²>,<k₁k₂>,<k₁k₃>,<k₂k₃>]||₂

其中<·>代表求集合平均，选择像素点周围3×3邻域的元素作为一个集合，||·||₂代表求向量的二范数；

(4)计算二阶统计量向量标准形式t_m到不变区域二阶标准向量t_u的投影矩阵P：

(4a)在不考虑相干斑噪声的前提下，设不变区域的输入向量k＝a[1,1,0]，a是一个与雷达回波强度有关的标量；

(4b)根据二阶统计量向量t的构造方式，计算出理想情况下不变区域二阶统计量向量t的形式为去掉系数记作二阶统计量向量t的标准形式t_m＝[1,1,0,1,0,0]^T。

(4c)将标准形式t_m投影到[1,0,0,0,0,0]^T方向，记作不变区域二阶标准向量t_u＝[1,0,0,0,0,0]^T；

(4d)计算出从t_m到t_u的投影矩阵P：

P＝E-2(w·w^T)

其中，E是一个6阶单位矩阵，w是一个中间变量，^T为转置操作符，

$E=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right],$ w＝(t_m-|t_m|·[1,0,0,0,0,0])/||t_m-|t_m|·[1,0,0,0,0,0]||_2，

|·|代表取模操作，||·||₂代表求向量的二范数；

(5)通过投影矩阵P对全部二阶统计量向量t进行线性投影，得到一组二阶向量t_p＝P·t，引入扰动参数，依次计算对每个坐标下的t_p与不变区域二阶标准向量t_u的相关系数γ：

$\mathrm{\&gamma;}=\frac{1}{\sqrt{1+R\frac{{\left|{t_p}^{*T}{t}_u\right|}^2}{{t_p}^{*T}{t}_p-{\left|{t_p}^{*T}{t}_u\right|}^2}}}$

其中，R是根据实际数据设置的扰动参数，取值为R＝0.5，^*为共轭操作符；

(6)每个坐标下的t_p都能计算出一个γ，将这些γ按t_p对应的坐标排列，构成差异图D；

(7)对上述差异图D做均值滤波，并对做均值滤波之后的差异图依次进行幂次变换和模糊C均值聚类，得到二值图像I，完成变化检测。

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

第一，本发明在获取差异图的时候,利用邻域信息引入了二阶统计量，有效地消除了由于存在相干斑噪声造成的虚警。

第二，本发明在获取差异图的时候引入了扰动参数，降低了由于非变化区域存在扰动造成的虚警。

附图说明

图1为本发明的实现流程图；

图2为Bern地区发生洪水前采集的SAR图像；

图3为Bern地区发生洪水后采集的SAR图像；

图4为参考Bern地区发生洪水前后实际地形标记的标准变化图；

图5为用现有的基于高斯统计模型方法对图2和图3检测得到的变化图；

图6为用现有的基于模糊C均值FCM方法对图2和图3检测得到的变化图；

图7为用本发明方法对图2和图3检测得到的变化图。

具体实施方式

参照图1，本发明的具体实现步骤包括如下：

步骤1，图像预处理。

对已进行配准、几何校正、辐射校正的某场景初次采集的SAR图像X和该场景出现变化后采集的SAR图像Y做差后，取模得到差图像Z。

Z＝|X-Y|，其中，|·|代表取模操作。

本实例使用图2做为图像X，图3做为图像Y。

步骤2，构造输入向量k。

根据初次采集的SAR图像X在(x，y)坐标处的像素值X_(x,y)，场景出现变化后采集的SAR图像Y在(x，y)坐标处的像素值Y_(x,y)，差图像Z在(x，y)坐标处的像素值Z_(x,y)，构造输入向量k：

k＝[X_(x,y),Y_(x,y),Z_(x,y)]^T；

令k₁,k₂,k₃分别与X_(x,y),Y_(x,y),Z_(x,y)一一对应，将k表示为：

k＝[k₁,k₂,k₃]^T。

步骤3，用每一个坐标(x，y)下的向量k，构造出一个二阶统计量向量t：

t＝[<k₁ ²>,<k₂ ²>,<k₃ ²>,<k₁k₂>,<k₁k₃>,<k₂k₃>]^T/||[<k₁ ²>,<k₂ ²>,<k₃ ²>,<k₁k₂>,<k₁k₃>,<k₂k₃>]||₂

其中<·>代表求集合平均，本实例用(x，y)坐标周围的3×3邻域中的输入向量k组成一个集合，||·||₂代表求向量的二范数。

步骤4，计算出不变区域二阶统计量向量的标准形式t_m。

在不考虑相干斑噪声的前提下，不变区域的输入向量k＝a[1,1,0]^T，a是一个与雷达回波强度有关的标量；

将k＝a[1,1,0]^T代入步骤3中t的表达式，计算出二阶统计量向量t的标准形式为由于本发明中只需要使用标准形式的方向信息，因此忽略幅度项得到二阶统计量向量t的标准形式：t_m＝[1,1,0,1,0,0]^T。

步骤5，求得一组二阶向量t_pr。

(5a)定义不变区域二阶标准向量t_u＝[1,0,0,0,0,0]^T，计算二阶统计量向量标准形式t_m＝[1,1,0,1,0,0]^T到不变区域二阶标准向量t_u＝[1,0,0,0,0,0]^T的投影矩阵P:

P＝E-2(w·w^T)

其中，E是一个6阶单位矩阵，w是一个中间变量，^T为转置操作符，

$E=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right],$ w＝(t_m-|t_m|·[1,0,0,0,0,0])/||t_m-|t_m|·[1,0,0,0,0,0]||_2，

|·|代表取模操作；

(5b)通过投影矩阵P对每一个坐标(x，y)下的二阶统计量向量t进行线性投影，得到每一个坐标(x，y)下的二阶向量t_pr：

t_pr＝P·t。

步骤6，引入扰动参数，生成差异图D。

(6a)引入扰动参数R，调节R可以在生成相关系数图时改变对扰动的容忍程度，较大的R表示可以容忍较大的扰动，较小的R表示可以容忍较小的扰动，本实例设置扰动参数R＝0.5。

(6b)计算每一个坐标(x，y)下的二阶向量t_pr与不变区域二阶标准向量t_u的相关系数γ_(x,y)：

其中t_pr(x,y)是在坐标(x,y)处求得的二阶向量t_pr；

将相关系数γ_(x,y)按坐标(x，y)排列，生成相关系数图像γ；

(6c)通过中值滤波器M_f对相关系数图像γ做中值滤波，得到差异图D，其中，中值滤波器M_f的形式为全1方阵：本实例取n＝3。

步骤7，应用模糊C均值聚类的方法对差异图D进行二值分割，实现变化检测。

模糊C均值聚类FCM，是用隶属度确定每个数据点属于某个聚类的程度的一种聚类算法，FCM采用模糊划分，使得每个给定数据点用值在0，1间的隶属度数值来确定其属于各类的程度，令每个数据点属于隶属度数值较大的类，实现聚类。具体步骤如下：

(7a)对差异图D进行幂次变换，拉伸数据的动态范围，其中幂次变换之后的差异图D_p在坐标(x,y)处的像素值为D_(x,y)代表差异图D在坐标(x,y)处的像素值，n是幂次数，本实例取值n＝6，得到变换后的差异图D_p；

(7b)使用Matlab软件中的fcm函数对变换后的差异图D_p中的像素点进行模糊2均值聚类操作，得到每个像素点对变化区域的隶属度数值和对不变区域的隶属度数值，并将这两个区域的隶属度数值进行比较：若像素点对变化区域的隶属度数值大于对不变区域的隶属度数值，则判定该像素点属于变化区域，反之则判定该像素点属于不变区域；

(7c)给属于变化区域的像素点标号为1，给属于不变区域标号为0，使用每个像素点的标号按像素点对应的位置组成一幅二值图像I，该图像I即变化检测的最终结果。

本发明的效果通过以下对实测数据的实验进一步说明：

1、实验场景：

本发明所用的实测数据为Bern地区发生洪水前采集的SAR图像，如图2所示，Bern地区发生洪水后采集的SAR图像，如图3所示，以及参考实际地形给出的标准变化图，如图4所示。

2、实验内容：

实验1，采用现有基于高斯统计模型检测方法对图2和图3进行检测，其变化检测结果如图5所示；

实验2，采用现有基于模糊C均值FCM聚类检测方法对图2和图3进行检测，其变化检测结果如图6所示；

实验3，采用本发明方法对图2和图3进行检测，其变化检测结果如图7所示。

3、实验结果分析

从图5、图6和图7可以看出，与现有的应用传统差异图的基于高斯统计模型检测方法和基于模糊C均值FCM聚类检测方法相比，本发明中利用SAR图像二阶统计量的相关系数构造差异图，有效地抑制了相干斑噪声，并引入扰动参数，能检测到更准确的变化区域。

将采用上述三种方法得到的变化图与参考实际地形给出的标准变化图4相比较，找出检测得到的变化图与标准变化图相同的部分，计算检测率，检测率的计算公式如下：

通过计算得到上述三种方法检测率比较，如表1所示：

表1 本发明方法与其他方法检测率的比较

由表1可以看出，本发明与现有的两种检测方法相比，其检测率具有提高。

Claims (3)

1.一种基于二阶统计量扰动分析的SAR图像变化检测方法，包括如下步骤：

(1)对采集场景初次的SAR图像X和该场景出现变化后的SAR图像Y依次进行配准、几何校正、辐射校正的预处理，并对这两幅图像做差取模值，得到差图像Z＝|X-Y|；

(2)用上述三个图像X、Y、Z中相同坐标的元素构造一组输入向量k＝[X_(x,y),Y_(x,y),Z_(x,y)]^T，并将k的形式简化为k＝[k₁,k₂,k₃]^T，

其中X_(x,y)是图像X在(x，y)坐标处的像素值，Y_(x,y),Z_(x,y)分别对应图像Y,Z在(x，y)坐标处的像素值，k₁,k₂,k₃分别与X_(x,y),Y_(x,y),Z_(x,y)一一对应，^T为转置操作符；

(3)根据向量k定义二阶统计量向量t：

t＝[<k₁ ²>,<k₂ ²>,<k₃ ²>,<k₁k₂>,<k₁k₃>,<k₂k₃>]^T/||[<k₁ ²>,<k₂ ²>,<k₃ ²>,<k₁k₂>,<k₁k₃>,<k₂k₃>]||₂

其中<·>代表求集合平均，||·||₂代表求向量的二范数；

(4)计算二阶统计量向量标准形式t_m到不变区域二阶标准向量t_u的投影矩阵P：

(4a)在不考虑相干斑噪声的前提下，设不变区域的输入向量k＝a[1,1,0]^T，a是一个与雷达回波强度有关的标量；

(4b)根据二阶统计量向量t的构造方式，计算出理想情况下不变区域二阶统计量向量t的形式为去掉系数记作二阶统计量向量t的标准形式t_m＝[1,1,0,1,0,0]^T。

(4c)将标准形式t_m投影到[1,0,0,0,0,0]^T方向，将[1,0,0,0,0,0]^T记作不变区域二阶标准向量t_u；

(4d)计算出从t_m到t_u的投影矩阵P：

P＝E-2(w·w^T)

其中，E是一个6阶单位矩阵，w是一个中间变量，^T为转置操作符，

$E=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right],$ w＝(t_m-|t_m|·[1,0,0,0,0,0])/||t_m-|t_m|·[1,0,0,0,0,0]||₂，

|·|代表取模操作，||·||₂代表求向量的二范数；

(5)通过投影矩阵P对全部二阶统计量向量t进行线性投影，得到一组二阶向量t_pr＝P·t，结合扰动分析，依次计算对每个坐标下的t_pr与不变区域二阶标准向量t_u的相关系数γ：

$\mathrm{\&gamma;}=\frac{1}{\sqrt{1+R\frac{{\left|{t_{\mathrm{pr}}}^{*T}{t}_u\right|}^2}{{t_{\mathrm{pr}}}^{*T}{t}_{\mathrm{pr}}-{\left|{t_{\mathrm{pr}}}^{*T}{t}_u\right|}^2}}}$

其中，R是根据实际数据设置的扰动分析参数，取值为R＝0.5，^*为共轭操作符，

(6)每个坐标下的t_pr都能计算出一个γ，将这些γ按t_pr对应的坐标排列，构成差异图D；

(7)对上述差异图D做均值滤波，并对做均值滤波之后的差异图依次进行幂次变换和模糊C均值聚类，得到二值图像I，完成变化检测。

2.根据权利要求1所述的方法，其中(4b)所述的根据二阶统计量向量t的构造方式，计算出理想情况下不变区域二阶统计量向量t的形式为是将步骤(4a)中所述的不变区域的输入向量k＝a[1,1,0]代入到步骤(3)中的二阶统计量向量公式：t＝[<k₁ ²>，<k₂ ²>,<k₃ ²>,<k₁k₂>,<k₁k₃>,<k₂k₃>]^T/||[<k₁ ²>,<k₂ ²>,<k₃ ²>,<k₁k₂>,<k₁k₃>,<k₂k₃>]||₂中，计算出二阶统计量向量

3.根据权利要求1所述的方法，其中所述(7)对均值滤波之后的差异图依次进行幂次变换和模糊C均值聚类，得到二值图像I，按如下步骤进行：

(7a)对差异图D进行幂次变换，得到变换后的差异图D_p，其中幂次变换之后的差异图D_p在坐标(x,y)处的像素值为：D_(x,y)代表差异图D在坐标(x,y)处的像素值，n是幂次数，取值为n＝6；

(7b)使用Matlab软件中的fcm函数对上述幂次变换之后的差异图D_p进行聚类操作，得到每个像素点对变化区域和不变区域两类区域的隶属度数值；

(7c)令每个像素点属于隶属度数值较大的那一类区域，给定不变区域的标号为0和变化区域的标号为1，使用每个像素点的标号组成二值图像I。