CN103065320A - 基于恒虚警阈值的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于恒虚警阈值的SAR图像变化检测方法，主要解决现有SAR图像变化检测虚警与阈值对应关系不明确的问题。其实现过程是：(1)对输入的两幅幅度SAR图像，计算基于块均值的比值图；(2)根据预期的虚警概率，按照比值分布模型确定所对应的阈值；(3)根据该阈值对比值图进行分割，得到变化检测结果图。本发明能够根据预先设定的虚警率计算出对应的阈值，进而得到符合预期虚警率的SAR图像变化检测结果；同时由于考虑了邻域信息，减少了噪声干扰，提高了检测精度，可用于资源和环境监测、农业研究、自然灾害监测和评估。

Description

基于恒虚警阈值的SAR图像变化检测方法

技术领域

本发明属于图像处理技术领域，涉及一种基于恒虚警阈值的SAR图像变化检测方法，可用于对SAR图像中感兴趣的地物变化信息进行检测处理。

背景技术

SAR图像变化检测旨在研究同一场景、不同时段的两幅或多幅SAR图像间的差异来得到感兴趣的地物变化信息。由于SAR***具有全天候、全天时、高分辨率及强大的穿透能力等独特优势，同时具有固定的重访周期，因此它比光学遥感更适合于变化检测，是很好的变化检测信息源。SAR图像变化检测技术被广泛应用于很多方面，例如资源和环境监测、农业研究、自然灾害监测和评估、军事分析等。

SAR图像变化检测技术是遥感信息处理领域的重要应用之一，它是通过分析同一地区不同时期拍摄的遥感图像间光谱特征差异或空间结构特性的差异，来获得人们所需要的地物类型的转变或内部条件和状态的变化。变化检测技术可以检测出不同时间图像灰度值或局部纹理之间的变化，在此基础上获得需要的目标在形状、位置、数量、及其他属性的变化情况。但是，SAR图像斑点噪声的影响依然是其应用的一个主要瓶颈。

经过近十几年的发展，针对不同类型的遥感图像，学者们提出了许多变化检测的技术和方法，可粗略的归类为：1：分类后比较法；2：基于简单运算法；3：基于图像建模法；4：基于空间领域特征法；5：基于特征的变化检测方法等。其中，近年来变化检测算法较为常用的图像差值法、比值法、均值比值法可归类于简单运算法。差值法：其主要过程是将两时相SAR图像对应像素点灰度值相减得到差异图，然后取0~255的阈值进行分割，得到变化和不变化的区域；图像比值法：其主要过程是计算两时相SAR图像对应像素点灰度值的比值差异图，如果像素未发生变化，则比值应接近于1，反之则远小于或大于1；均值比值法：其主要过程是运用两时相SAR图像对应像素点以及其邻域像素点的均值的比值得到差异图，此方法与比值法类似。差值法算法简单，但易受SAR成像质量等客观条件的影响，易产生“伪变化”信息；比值法和均值比值法，均对SAR图像的乘性噪声不敏感，因此主要应用于乘性噪声的SAR图像。

综上所述的方法，虽然在很多方面都得到了应用，但仍然存在着不少的困难和问题：（1）传统的方法中没有一个虚警和阈值之间相互对应的模型，实际应用中无法得到检测结果的虚警概率；

（2）基于点的方法没有考虑邻域信息，受噪声干扰大，检测精度不高。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于恒虚警阈值的SAR图像变化检测方法，以调整虚警概率，满足实际需求；减小噪声干扰，提高检测精度。

实现本发明目的的技术方案是：根据斑点噪声的比值分布，对比值概率公式进行推导，从而得到关于虚警率和阈值的关系式。通过给定的虚警率确定阈值，并通过调整虚警率满足不同的实际需求，具体步骤包括如下：

(1)输入视数为L的两副幅度SAR图像V₁和V₂，分别以幅度图像V₁、V₂中的像素点v_1,k、v_2,k为中心，取M×M大小的邻域块u_1，k、u_2,k，计算每两个像素点对应的块均值比值r_k，从而得到比值图R；

(2)确定所要求的虚警概率α对应的阈值T；

2a)根据如下比值分布概率公式计算每个比值r_k出现的概率p(r_k)：

$p\left(r_k\right)=\frac{2\mathrm{\Γ}\left(\mathrm{NL}\right)}{{\mathrm{\Γ}}^2\left(\mathrm{NL}\right)}\frac{{2r}_k^{2\mathrm{NL}-1}}{{(r_k^2+1)}^{2\mathrm{NL}}},$ r_k∈[0,1]

其中，Γ(·)表示gamma函数，点数N＝M×M，NL表示块内点数N与视数L的乘积；

2b)根据概率p(r_k)，得到比值概率函数P(T)：

$P\left(T\right)={\∫}_0^{T}p\left(r_k\right)d{r}_k$

其中，T是门限阈值；

2c)根据比值概率函数P(T)，由概率统计性质得到虚警概率α：

$\mathrm{\α}=\frac{{\∫}_0^{T}p\left(r_k\right)d{r}_k}{{\∫}_0^{1}p\left(r_k\right)d{r}_k}$

对上式进行化简得到T和α的线性表示：

α＝bT

其中常量b满足关系式：

$\mathrm{bT}=\frac{{\∫}_0^T\frac{r_k^{2\mathrm{NL}-1}d{r}_k}{{(r_k^2+1)}^{2\mathrm{NL}}}}{{\∫}_0^1\frac{r_k^{2\mathrm{NL}-1}d{r}_k}{{(r_k^2+1)}^{2\mathrm{NL}}}}$

且NL已由步骤2a)求得；

2d)再根据给定的虚警概率α，计算所需要的门限阈值T：

$T=\frac{\mathrm{\α}}{b}$

(3)根据所求得的阈值T，遍历整个比值图像R，将R中像素值小于阈值T的点置为255，将像素值大于等于阈值T的点置为0，从而得到变化检测结果图。

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

1.本发明能在空域中进行，实现过程简单，且可以并行实现；

2.本发明通过虚警率确定阈值，进而通过所确定的阈值得到检测结果，能够满足不同实际需求，且检测结果更具实用价值，更加灵活；

3.本发明考虑了邻域信息，减少了噪声干扰，检测结果更加稳定，精度更高。

附图说明

图1是本发明的流程图；

图2是本发明仿真使用的两幅待合成的lena原始图和标准图；

图3是用本发明方法对基于图2合成的SAR图像进行变化检测的结果图；

图4是本发明仿真使用的两幅Ottawa真实SAR影像图；

图5是用本发明方法对基于图4进行变化检测的结果图；

图6是本发明仿真使用的两幅待合成的box原始图和标准图；

图7是用本发明方法对基于图6合成的SAR图像进行变化检测的结果图。

具体实施方式

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

步骤1，输入视数为L的两副幅度SAR图像V₁和V₂，计算每两个像素点对应的块均值比值r_k，得到比值图R。

1a)分别以像素点v_1，k、v_2,k为中心，各取M×M大小的领域块u_1，k、u_2,k；

1b)分别计算上述块u_1，k、u_2,k的均值

称为块均值；其中，

${\stackrel{\‾}{v}}_{1,k}=\mathrm{mean}\left(u_{1,k}\right),$

${\stackrel{\‾}{v}}_{2,k}=\mathrm{mean}\left(u_{2,k}\right)$

式中，mean(·)表示均值函数；

1c)计算上述块均值的比值r_k：

$r_k=\min \{\frac{{\stackrel{\‾}{v}}_{1,k}}{{\stackrel{\‾}{v}}_{2,k}},\frac{{\stackrel{\‾}{v}}_{2,k}}{{\stackrel{\‾}{v}}_{1,k}}\},$ r_k∈[0,1]；

1d)将上述比值r_k作为比值图像R中第k个像素点的值。

步骤2，确定所要求的虚警概率α对应的阈值T；

2a)根据如下比值分布概率公式计算每个比值r_k出现的概率p(r_k)：

$p\left(r_k\right)=\frac{2\mathrm{\Γ}\left(\mathrm{NL}\right)}{{\mathrm{\Γ}}^2\left(\mathrm{NL}\right)}\frac{{2r}_k^{2\mathrm{NL}-1}}{{(r_k^2+1)}^{2\mathrm{NL}}},$ r_k∈[0,1]

其中，Γ(·)表示gamma函数，点数N＝M×M，NL＝N×L表示块内点数N与视数L的乘积；

2b)根据概率p(r_k)，得到比值概率函数P(T)：

$P\left(T\right)={\∫}_0^{T}p\left(r_k\right)d{r}_k,$

其中，T是门限阈值；

2c)根据比值概率函数P(T)，由概率统计性质得到虚警概率α：

$\mathrm{\α}=\frac{{\∫}_0^{T}p\left(r_k\right)d{r}_k}{{\∫}_0^{1}p\left(r_k\right)d{r}_k},$

对上式进行简化得到T与α的线性表示：

α＝bT，

其中常量b满足关系式：

$\mathrm{bT}=\frac{{\∫}_0^T\frac{r_k^{2\mathrm{NL}-1}d{r}_k}{{(r_k^2+1)}^{2\mathrm{NL}}}}{{\∫}_0^1\frac{r_k^{2\mathrm{NL}-1}d{r}_k}{{(r_k^2+1)}^{2\mathrm{NL}}}},$

且NL已由步骤2a)求得；

2d)根据给定的虚警概率α，计算所需要的门限阈值T：

$T=\frac{\mathrm{\α}}{b},$

步骤3，根据所求得的阈值T，构造变化检测结果图：

利用该阈值T遍历整个比值图像R，将R中比值小于阈值T的点设为255，比值大于等于阈值T的点设为0，进而得到变化检测结果图。

本发明效果可以通过以下实验进一步证实：

一.实验条件和内容

实验条件：实验所使用的仿真图像如下：

如图2所述的测试仿真图像lena，其中图2(a)和2(b)均为待合成的lena图，图2(c)为检测标准图；

基于图2(a)和2(b)合成视数L=1的SAR图像，如图3(a)所述；

如图4所述的真实图像Ottawa，其中图4(a)表示1997年5月该地区地貌信息，图4(b)表示1997年8月该地区地貌信息；

如图6所述的测试仿真图像box，其中图6(a)和6(b)均为待合成的box图，图6(c)为检测标准图；

基于图6(a)和6(b)合成视数L=1的SAR图像，如图7(a)所述。

实验内容：在上述实验条件下，使用本发明方法对合成的SAR图像3(a)和7(a)、真实SAR图像4进行检测，所有实验均在同质区域的假设下进行。

二．实验结果

实验一：使用本发明方法对图3(a)进行处理，其中块大小分别取3×3、5×5，采用虚警率为0.001对应的阈值和虚警率为0.0015对应的阈值。实验检测结果如图3(b)、3(c)所示，其中图3(b)和图3(c)从左到右分别表示块大小为3×3、5×5的检测结果。

由图3(b)和图3(c)可见，对于所给定的虚警率，本发明方法能够有效的检测合成SAR图像的变化部分；并且，在虚警率恒定的情况下，所采用块越小，检测的杂点块越小，在块变大的情况下，杂点虽变少，但杂点块变大，不利于有效检测。

实验二：使用本发明方法对图4进行处理，其中块大小分别取3×3、5×5，采用虚警率为0.001对应的阈值和虚警率为0.0015对应的阈值。实验检测结果如图5所示，其中图5(a)和5(b)从左到右分别表示块大小为3×3、5×5的检测结果。

由图5可见，对于所给定的虚警率，本发明方法能够有效的检测真实SAR图像的变化部分。

实验三：使用本发明方法和现有的点的方法对图7(a)进行处理，本发明方法所采用的虚警率为0.001，块大小为3×3。实验检测结果如图7(b)所示，其中图7(b)从左到右分别为点的方法和本发明方法的检测结果。

由7(b)的结果图可以看出，基于点的方法不能有效的检测出变化区域，而本发明方法能较好的检测出变化的区域。

实验四：使用本发明方法确定不同点数N、视数L、虚警率所对应的阈值,并对图3(a)进行变化检测，检测结果见表1；使用本发明方法得到的基于点和基于块的方法进行变化检测，检测结果对比见表2。

表1  不同点数N、视数L的实验与理论虚警率对比

表2  本发明方法和基于点的方法的检测结果对比

表1和表2的结果均为在同质区域的假设下得到。

由表1可以看出，在同视数同点数条件下，两虚警率之差相比基本符合，这与理论分析符合。虽然实验值和理论值有一定的误差，但在允许的波动范围内。这是由于计算基于块均值的比值时，块内存在与中心像素点非同质的像素点，导致门限的计算存在一定的误差。这个误差随块的尺寸增大而变大。本发明使用较小的块尺寸，例如3×3。

由表2可以看出，与基于点的方法相比，本发明方法考虑了邻域信息，减少了噪声干扰，检测结果更加稳定，精度更高。

在以上实验条件下采用本发明所提出的方法进行SAR图像的变化检测。从算法实现上可以看出，与现有的阈值变换检测方法相比，使用本发明方法更快捷简便，并能满足不同实际需求，且检测结果更具实用价值，更加灵活。

Claims (2)

1.一种基于恒虚警阈值的SAR图像变化检测方法，包括如下步骤：

(1)输入视数为L的两副幅度SAR图像V₁和V₂，分别以幅度图像V₁、V₂中的像素点v_1,k、v_2,k为中心，取M×M大小的邻域块u_1,k、u_2，k，计算每两个像素点对应的块均值比值r_k,得到比值图R；

(2)确定所要求的虚警概率α对应的阈值T；

2a)根据如下比值分布概率公式计算每个比值r_k出现的概率p(r_k)：

$p\left(r_k\right)=\frac{2\mathrm{\Γ}\left(\mathrm{NL}\right)}{{\mathrm{\Γ}}^2\left(\mathrm{NL}\right)}\frac{{2r}_k^{2\mathrm{NL}-1}}{{(r_k^2+1)}^{2\mathrm{NL}}},$ r_k∈[0,1]

其中，Γ(·)表示gamma函数，点数N＝M×M，NL表示块内点数N与视数L的乘积；

2b)根据概率p(r_k)，得到比值概率函数P(T)：

$P\left(T\right)={\∫}_0^{T}p\left(r_k\right)d{r}_k$

其中，T是门限阈值；

2c)根据比值概率函数P(T)，由概率统计性质得到虚警概率α：

$\mathrm{\α}=\frac{{\∫}_0^{T}p\left(r_k\right)d{r}_k}{{\∫}_0^{1}p\left(r_k\right)d{r}_k}$

对上式进行化简得到T和α的线性表示：

α＝bT

其中常量b满足关系式：

$\mathrm{bT}=\frac{{\∫}_0^T\frac{r_k^{2\mathrm{NL}-1}d{r}_k}{{(r_k^2+1)}^{2\mathrm{NL}}}}{{\∫}_0^1\frac{r_k^{2\mathrm{NL}-1}d{r}_k}{{(r_k^2+1)}^{2\mathrm{NL}}}}$

且NL已由步骤2a)求得；

2d)再根据给定的虚警概率α，计算所需要的门限阈值T：

$T=\frac{\mathrm{\α}}{b}$

(3)根据所求得的阈值T，遍历整个比值图像R，将R中像素值小于阈值T的点置为255，将像素值大于等于阈值T的点置为0，从而得到变化检测结果图。

2.根据权利要求1所述的基于恒虚警阈值的SAR图像变化检测方法，其特征在于步骤(1)所述的计算每两个像素点对应的块均值比值r_k,按照如下步骤进行：

2a)分别以像素点v_1，k、v_2,k为中心，各取M×M大小的邻域块u_1,k、u_2,k；

2b)分别计算上述邻域块u_1，k、u_2,k的均值

称为块均值；其中，

${\stackrel{\‾}{v}}_{1,k}=\mathrm{mean}\left(u_{1,k}\right),$ ${\stackrel{\‾}{v}}_{2,k}=\mathrm{mean}\left(u_{2,k}\right),$

其中，mean(·)表示均值函数；

2c)计算上述块均值

的比值r_k：

$r_k=\min \{\frac{{\stackrel{\‾}{v}}_{1,k}}{{\stackrel{\‾}{v}}_{2,k}},\frac{{\stackrel{\‾}{v}}_{2,k}}{{\stackrel{\‾}{v}}_{1,k}}\},$ r_k∈[0,1]

2d)将上述比值r_k作为比值图像R中第k个像素点的值。

Images