CN101398487A - 一种降低合成孔径雷达影像斑点噪声的方法 - Google Patents

Abstract

本发明公开了一种降低合成孔径雷达影像(SAR：Synthetic Aperture Radar)斑点噪声的方法，包括进行图像对数变换，将乘性噪声转变为加性噪声；对变换后的图像进行均值校正，使得校正的结果满足零均值高斯白噪声的假设前提；对经过均值校正的图像进行循环平移，从而消除Contourlet变换可能造成的图像可视效果畸变；对一系列得到的平移影像分别应用Contourlet变换，使用Contourlet域滤波器进行滤波处理；对滤波过的各个平移影像进行逆Contourlet变换，然后将所有经过逆变换的影像分别进行逆平移变换，再将这些经过逆平移变换的影像求平均；对平均影像进行指数变换，得到滤波后的影像，本发明可以极大地降低SAR影像斑点噪声。

Description

一种降低合成孔径雷达影像斑点噪声的方法

技术领域

本发明涉及一种降低斑点噪声的方法，特别是涉及一种降低合成孔径雷达影像斑点噪声的方法。

背景技术

合成孔径雷达(SAR)已经成为获取地球表面信息的重要资源，在越来越的部门得到了广泛的应用，如地质探矿、农业估产、水灾监测、城市扩展、生态环境调查等。由于采用主动式相干波成像方式，在从原始信号到图像的重建过程中，存在一个固有的问题，即斑点噪声(Speckle)的干扰。斑点噪声在降低SAR图像空间分辨率的同时，也干扰了对SAR后向散射系数的测量，严重影响了SAR图像的应用。

SAR图像有很多不足之处，最突出的问题是雷达的成像难以理解。雷达图像形成的是一种距离投影，与我们所常见的可见光影像的成像规律不同，因此人们必须经过训练才能准确辨认雷达图像所传达的信息。虽然现在有人开始否定这种说法，认为现在的SAR图像已经非常完美，与普通照片的相像程度已经达到85％，但是多数专家认为SAR刚开始由抽象图像到直观图像的转变。另一个不足之处是，SAR对平台轨道的精度要求高，对雷达设备和技术要求的门槛也很高，导致SAR的研制和应用的费用很高。回顾SAR发展的历史，其寿命周期费用一直比可见光和红外系统高。但现在这种情况已经开始改变，有专家认为目前SAR的寿命周期费用已经下降到与红外及较复杂的可见光系统相当的程度。

综上所述，可以发现，目前世界各国普遍高度重视SAR传感器的研发与应用。21世纪，合成孔径雷达遥感技术将朝着多平台、高分辨率、多波段、多极化和多工作模式等方向迅速发展。

发明内容

本发明目的在于公开一种降低合成孔径雷达影像斑点噪声的方法。

本发明目的通过如下技术方案实现的：

步骤1：进行图像对数变换，将乘性噪声转变为加性噪声；

步骤2：对变换后的图像进行均值校正，使得校正的结果满足零均值高斯白噪声的假设前提；

步骤3：对经过均值校正的图像进行循环平移，从而消除Contourlet变换可能造成的图像可视效果畸变；

步骤4：对一系列得到的平移影像分别应用Contourlet变换，使用Contourlet域滤波器进行滤波处理；

步骤5：对滤波过的各个平移影像进行逆Contourlet变换，然后将所有经过逆变换的影像分别进行逆平移变换，再将这些经过逆平移变换的影像求平均；

步骤6：对平均影像进行指数变换，得到滤波后的影像。

其中步骤1乘性噪声转化为加性噪声的方法：

乘性噪声模型来描述SAR图像，如下式(i)所示：

g＝f·e         (i)

其中，g代表观测值，也就是SAR图像；f代表真实值；e代表噪声。

对数变换(LOG)是一种有效的方法，如式(ii)：

δ＝K+ε       (ii)

其中，δ＝LOG(g)，K＝LOG(f)，ε＝LOG(e)。

进行完对数变换之后要进行Contourlet变换，由于Contourlet变换是一种线性变换，用下式表示：

y＝x+n        (iii)

其中，y＝CT(δ)代表含噪声的Contourlet系数，x＝CT(κ)表示原始信号的Contourlet系数，n＝CT(ε)表示噪声的Contourlet系数。CT(·)代表Contourlet变换；

其中步骤2：对变换后的图像进行均值校正：在进行指数变换之前，从Contourlet逆变换后的影像中减去经过对数变换的噪声影像的均值即得；

其中步骤3：循环平移方法的实现过程如下，让图像在水平方向以一定的步长S重复移动K₁次，同样在垂直方向重复移动K₂次，得到K₁×K₂个循环平移的图像。

其中步骤4滤波算法：阈值(threshold)方法(硬阈值或软阈值)或基于贝叶斯估计的系数收缩(shrinkage)方法或基于隐马尔可夫树模型(HMT)的滤波算法或基于上下文的隐马尔可夫模型的滤波算法。

其中基于隐马尔可夫树模型(HMT)的滤波算法为：

首先，对y建立HMT模型，并运用EM算法对其参数进行估计，获取其模型参数 ${\mathrm{\θ}}_y=\{p_{S_1}\left(m\right),{\mathrm{\ϵ}}_{i,m}^{m,r},{\mathrm{\σ}}_{i,m}^2\}$ (假设Contourlet系数服从由M个0均值的混合高斯分布)，然后，通过将噪声的方差从θ_y中去除，可以获得去噪的Contourlet系数的参数θ_x，即：

${\left({\mathrm{\σ}}_{(j,k,i),m}^{\left(x\right)}\right)}^2={[{\left({\mathrm{\σ}}_{(j,k,i),m}^{\left(y\right)}\right)}^2-{\left({\mathrm{\σ}}_{(j,k,i)}^{\left(n\right)}\right)}^2]}_{+}---\left(\mathrm{IV}\right)$

其中，[x]₊表示如果x大于等于0，则[x]₊＝x，否则[x]₊＝0。(j，k，i)表示在尺度为j方向为k的子波段中的第i个Contourlet系数，m为Contourlet系数所属的状态。其中，噪声的系数的方差

可以通过对最精细尺度的Contourlet系数进行中值估计获得，或者通过蒙特卡罗(Monte-Carlo)方法获得。

当获得了去噪影像的HMT参数θ_x之后，利用该参数对去噪影像的Contourlet系数进行估计得到：

$E\left[x_{j,k,i}\right|y_{j,k,i,}{\mathrm{\θ}}_x]=\underset{m}{\mathrm{\Σ}}p(S_{j,k,i}=m|y_{j,k,i},{\mathrm{\θ}}_x)\×\frac{{\left({\mathrm{\σ}}_{(j,k,i),m}^{\left(x\right)}\right)}^2}{{\left({\mathrm{\σ}}_{(j,k,i),m}^{\left(x\right)}\right)}^2+{\left({\mathrm{\σ}}_{(j,k,i)}^{\left(n\right)}\right)}^2}\×y_{j,k,i}---\left(V\right)$

本发明在图像(SAR)乘性斑点噪声模型的基础上，通过分别应用对数变换将乘性噪声转化为加性噪声，应用循环平移方法消除图像可视效果畸变，使用均值校正解决对数变换后的均值偏移，为使用各种不满足平移不变性变换的算法(如wavelet，contourlet)进行SAR图像滤波处理建立了统一的滤波处理框架；在此框架之下，应用几种常用的SAR滤波算法，对模拟影像和真实SAR影像进行了滤波处理，通过比较各种方法的可视效果和统计指标，表明在SAR滤波处理统一框架下，基于contourlet域隐马尔可夫树模型的去噪方法在滤波结果在可视效果和统计指标方面均优于常用的SAR滤波算法。

下述实验例和实施例用于进一步说明但不限于本发明；

实验例1 模拟图像降低影像斑点噪声实验

选用模拟图像对各种滤波器的效果进行评估。模拟影像是向一幅QuickBird影像中加入方差为0.05(该值已经过归一化，即在加入噪声时，假设图像的值在0～1之间)的服从Gamma分布的speckle噪声生成。图像大小为256像素×256像素。滤波方法包括基于斑点噪声模型的Lee滤波算法(Lee J S.Speckle Suppression and Analysis for Synthetic ApertureRadar Images[J].Opt.Eng.1986，25(5)：636-643.)(滤波器窗口为5×5)、小波软阈值法(WST)和CHMT方法---基于K均值粗分类的绑定法(CC-CHMT，本发明实施例1中所述的滤波算法)。CHMT模型训练中均使用零均值、2状态的混合高斯分布，LP分解选用’9-7’，方向滤波器组均选用pkva(Do M N，Vetterli M.The Contourlet Transform：An EfficientDirectional Multiresolution Image Representation[J].IEEETransactions on Image Processing.2005，14(12)：2091-2106.)，Contourlet分解层数为3，每层的方向数分别为4，8，16；小波域滤波器中，为了避免使用不同的小波基对结果产生的影响，小波基选用与Contourlet变换相同的’9-7’，分解层数也为3层。

图1-5是对模拟影像使用不同滤波器的处理结果。图1为原始QuickBird遥感影像，图2为加入噪声后的QuickBird影像，图3为使用Lee滤波器进行处理得到的结果，图4为使用小波软阈值法(WST)处理的结果，图5是使用基于K均值粗分类绑定方法(CC-CHMT)得到的结果。原始QuickBird影像加入噪声后退化的比较严重。从滤波结果来看，Lee方法虽然细节信息保持的比较好，但噪声压制的能力明显低于其他方法；WST方法产生了较严重的“振铃”效应，丧失了很多细节信息。

从总体上来说，基于Contourlet域的HMM方法得到的结果可视效果好于Lee方法和WST方法，其结果在斑点滤除和细节的保留上均取得了较好的效果，其可视效果也更接近于原始的Quickbird影像。

为了进一步研究各种滤波方法的滤波能力，在图像上任取一行，通过其剖面(profile)对滤波后影像与原始影像进行比较。Profile的位置如图1中的直线所示。这条线经过了很多边缘，便于我们比较不同的滤波器在处理边缘时的能力。图6显示了原始影像、Lee方法、WST方法、CC-CHMT方法的profile。可以看到，在Lee方法的曲线上产生了很多小的“尖峰”，与原始图像的曲线不能很好的吻合，表现在图像上就是噪声压制的效果不强；CST方法是与原始图像重叠的最好的，但从图4中可以发现，该方法得到的图像有较强的“振铃”效应。

实验例2 SAR影像利用不同滤波器实验

SAR影像选用Radarsat影像(C波段，HH极化，空间分辨率25m×25m)。图像大小为256像素×256像素。图7(a)为Radarsat影像；图7(b)为使用Lee滤波器进行处理得到的结果，图7(c)为使用小波软阈值法(WST)处理的结果，图7(d)是使用CHMT得到的结果。总得来说，各种方法的表现与使用模拟影像的结果较为近似。Lee方法在滤除噪声的同时，较大地模糊了影像(比如左上角的高反射区域经过滤波后，几乎合并成一片)，纹理细节信息丢失较严重；WST方法的视觉效果较差，有明显的“振铃”现象；CHMT方法的在有效滤除噪声的同时，都具备较强的细节保留能力，滤波后的图像比较自然、清晰。

我们使用等效视数(Equivalent Number of Looks，ENL)和信噪比(Signal Noise Ration，SNR)这2个参数来定量评价滤波器的性能。ENL可以评价一个滤波器的斑点抑制能力，ENL越大，表明滤波器对噪声的压制效果越强，在匀质区域，ENL定义为该区域均值的平方除以该区域的方差。SNR表示含噪声图像或者滤波后图像与原始图像相比，质量变化的情况。信噪比越高，表示滤波的效果越好。

表1 滤波器性能比较

 

	Quickbird	Radarsat
			ImageNoisy_imageLeeWSTCHMT	ENL        SNR16.0650    N/A45.8598    9.2614241.4610   13.6235180.1253   11.5247	ENL       SNR9.5804    N/A18.5187   29.758413.3264   32.036758.2718   26.6304

表1列出了对模拟影像和SAR影像利用不同滤波器进行滤波后的参数指标。从模拟影像的统计结果可以看到，Lee方法不论是ENL还是SNR都远低于其他几种方法；WST方法有较高的ENL和SNR，但是通过观察发现，这两种方法都在很大程度上模糊了影像，也就是说，由于该方法具有很强的“平滑”能力，才导致了较高的ENL和SNR值；从指标值来看，CHMT的方法能取得较好的效果。

下述实施例均能够实现上述实验例所述的效果

图1为原始QuickBird遥感影像

图2为加入噪声后的QuickBird影像

图3为使用Lee滤波器进行处理得到的结果

图4为使用小波软阈值法(WST)处理的结果

图5是使用基于K均值粗分类绑定方法(CC-CHMT)得到的结果。

图6不同滤波方法的剖面图

图7(a)为Radarsat影像

图7(b)为使用Lee滤波器进行处理得到的结果

图7(c)为使用小波软阈值法(WST)处理的结果

图7(d)是使用CHMT得到的结果

图8是本发明实施例1的流程图

Y代表含噪SAR图像，X代表滤波处理后的图像；LOG代表对数变换，R代表均值校正处理，Z_i和Z_-i分别代表进行步长为i的循环平移变换和循环平移逆变换，CT和CT^-1分别代表Contourlet变换和逆变换，w代表滤波处理算法，EXP代表指数变换。

具体实施方式

实施例1：

步骤1乘性噪声转化为加性噪声的方法：

乘性噪声模型来描述SAR图像，如下式(i)所示：

g＝f·e    (i)

其中，g代表观测值，也就是SAR图像；f代表真实值；e代表噪声。

对数变换(LOG)是一种有效的方法，如式(ii)：

δ＝K+ε      (ii)

其中，δ＝LOG(g)，K＝LOG(f)，ε＝LOG(e)。

进行完对数变换之后要进行Contourlet变换，由于Contourlet变换是一种线性变换，用下式表示：

y＝x+n        (iii)

其中，y＝CT(δ)代表含噪声的Contourlet系数，x＝CT(κ)表示原始信号的Contourlet系数，n＝CT(ε)表示噪声的Contourlet系数。CT(·)代表Contourlet变换；

步骤2：对变换后的图像进行均值校正：在进行指数变换之前，从Contourlet逆变换后的影像中减去经过对数变换的噪声影像的均值即得；

步骤3：循环平移方法的实现过程如下，让图像在水平方向以一定的步长S重复移动K₁次，同样在垂直方向重复移动K₂次，得到K₁×K₂个循环平移的图像。

步骤4滤波算法：其中基于隐马尔可夫树模型(HMT)的滤波算法为：

首先，对y建立HMT模型，并运用EM算法对其参数进行估计，获取其模型参数 ${\mathrm{\θ}}_y=\{p_{S_1}\left(m\right),{\mathrm{\ϵ}}_{i,m}^{m,r},{\mathrm{\σ}}_{i,m}^2\}$ (假设Contourlet系数服从由M个0均值的混合高斯分布)，然后，通过将噪声的方差从θ_y中去除，可以获得去噪的Contourlet系数的参数θ_x，即：

${\left({\mathrm{\σ}}_{(j,k,i),m}^{\left(x\right)}\right)}^2={[{\left({\mathrm{\σ}}_{(j,k,i),m}^{\left(y\right)}\right)}^2-{\left({\mathrm{\σ}}_{(j,k,i)}^{\left(n\right)}\right)}^2]}_{+}---\left(\mathrm{IV}\right)$

其中，[x]₊表示如果x大于等于0，则[x]₊＝x，否则[x]₊＝0。(j，k，i)表示在尺度为j方向为k的子波段中的第i个Contourlet系数，m为Contourlet系数所属的状态。其中，噪声的系数的方差

可以通过对最精细尺度的Contourlet系数进行中值估计获得，或者通过蒙特卡罗(Monte-Carlo)方法获得。

当获得了去噪影像的HMT参数θ_x之后，利用该参数对去噪影像的Contourlet系数进行估计得到：

$E\left[x_{j,k,i}\right|y_{j,k,i,}{\mathrm{\θ}}_x]=\underset{m}{\mathrm{\Σ}}p(S_{j,k,i}=m|y_{j,k,i},{\mathrm{\θ}}_x)\×\frac{{\left({\mathrm{\σ}}_{(j,k,i),m}^{\left(x\right)}\right)}^2}{{\left({\mathrm{\σ}}_{(j,k,i),m}^{\left(x\right)}\right)}^2+{\left({\mathrm{\σ}}_{(j,k,i)}^{\left(n\right)}\right)}^2}\×y_{j,k,i}---\left(V\right)$

步骤5：对滤波过的各个平移影像进行逆Contourlet变换，然后将所有经过逆变换的影像分别进行逆平移变换，再将这些经过逆平移变换的影像求平均；

步骤6：对平均影像进行指数变换，得到滤波后的影像。

上述整个算法的流程可以用示意图8表示。

Claims (6)

1、一种降低SAR影像斑点噪声的方法，其特征在于包括如下步骤：

步骤1：进行图像对数变换，将乘性噪声转变为加性噪声；

步骤2：对变换后的图像进行均值校正，使得校正的结果满足零均值高斯白噪声的假设前提；

步骤3：对经过均值校正的图像进行循环平移，从而消除Contourlet变换可能造成的图像可视效果畸变；

步骤4：对一系列得到的平移影像分别应用Contourlet变换，使用Contourlet域滤波器进行滤波处理；

步骤5：对滤波过的各个平移影像进行逆Contourlet变换，然后将所有经过逆变换的影像分别进行逆平移变换，再将这些经过逆平移变换的影像求平均；

步骤6：对平均影像进行指数变换，得到滤波后的影像。

2、如权利要求1所述的方法，其中步骤1乘性噪声转化为加性噪声的方法为：

乘性噪声模型来描述图像，如下式(i)所示：

g＝f·e       (i)

其中，g代表观测值，也就是图像；f代表真实值；e代表噪声；

对数变换LOG，如式(ii)：

δ＝κ+ε        (ii)

其中，δ＝LOG(g)，κ＝LOG(f)，ε＝LOG(e)；

进行Contourlet变换，用(iii)式表示：

y＝x+n       (iii)

其中，y＝CT(δ)代表含噪声的Contourlet系数，x＝CT(κ)表示原始信号的Contourlet系数，n＝CT(ε)表示噪声的Contourlet系数；CT(·)代表Contourlet变换。

3、如权利要求1所述的方法，其中步骤2中对变换后的图像进行均值校正，在进行指数变换之前，从Contourlet逆变换后的影像中减去经过对数变换的噪声影像的均值即得。

4、如权利要求1所述的方法，其中步骤3循环平移方法的步骤如下：使图像在水平方向以一定的步长S重复移动K₁次，同样在垂直方向重复移动K₂次，得到K₁×K₂个循环平移的图像。

5、如权利要求1所述的方法，其中步骤4滤波算法为阈值方法或基于贝叶斯估计的系数收缩方法或基于隐马尔可夫树模型的滤波算法或基于上下文的隐马尔可夫模型的滤波算法。

6、如权利要求5所述的方法，其中基于隐马尔可夫树模型的滤波算法为：

对y建立HMT模型，并运用EM算法对其参数进行估计，获取其模型参数 ${\mathrm{\θ}}_y=\{p_{s_1}\left(m\right),{\mathrm{\ϵ}}_{i,m}^{m,r},{\mathrm{\σ}}_{i,m}^2\};$ 通过将噪声的方差从θ_y中去除，获得去噪的Contourlet系数的参数θ_x，即：

${\left({\mathrm{\σ}}_{(j,k,i),m}^{\left(x\right)}\right)}^2={[{\left({\mathrm{\σ}}_{(j,k,i),m}^{\left(y\right)}\right)}^2-{\left({\mathrm{\σ}}_{(j,k,i)}^{\left(n\right)}\right)}^2]}_{+}---\left(\mathrm{IV}\right)$

其中，[x]₊表示如果x大于等于0，则[x]₊＝x，否则[x]₊＝0；(j，k，i)表示在尺度为j方向为k的子波段中的第i个Contourlet系数，m为Contourlet系数所属的状态；

其中，噪声的系数的方差

通过对最精细尺度的Contourlet系数进行中值估计获得，或者通过蒙特卡罗(Monte-Carlo)方法获得；

获得去噪影像的HMT参数θ_x后，利用该参数对去噪影像的Contourlet系数进行估计得到：

$E\left[x_{j,k,i}\right|y_{j,k,i},{\mathrm{\θ}}_x]=\underset{m}{\mathrm{\Σ}}p(S_{j,k,i}=m|y_{j,k,i},{\mathrm{\θ}}_x)\×\frac{{\left({\mathrm{\σ}}_{(j,k,i),m}^{\left(x\right)}\right)}^2}{{\left({\mathrm{\σ}}_{(j,k,i),m}^{\left(x\right)}\right)}^2+{\left({\mathrm{\σ}}_{(j,k,i)}^{\left(n\right)}\right)}^2}\×y_{j,k,i}---\left(V\right)$

Images