CN103971346A - 基于稀疏域噪声分布约束的sar图像抑斑方法 - Google Patents

Abstract

本发明公开了一种基于稀疏域噪声分布约束的SAR图像抑斑方法，主要解决现有抑斑方法不能同时保持图像纹理细节及辐射特性的问题。主要步骤为：1.对原始SAR图像进行对数变换，将变换后的图像划分成重叠的大小相同子图像块，并用K均值聚类的方法对图像块进行聚类；2.用主成分分析的方法求出每一个聚类的稀疏表示子字典，并用非局部均值算法计算每个子图像块对应的无噪子图的稀疏表示系数的估计值；3.用软阈值收缩算法求得最终的无噪稀疏表示系数，并通过反稀疏变换得到去噪后的子图像块；4.拼接子图像块并使用反对数变换得到抑斑后的SAR图像。本发明相干斑抑制能力强，图像的纹理细节及辐射特性保持好，可用于目标识别与特征提取。

Description

基于稀疏域噪声分布约束的SAR图像抑斑方法

技术领域

本发明属于图像处理技术领域，具体地说是一种合成孔径雷达SAR图像抑斑方法，可用于目标识别与特征提取。

背景技术

合成孔径雷达是一种主动式微波遥感器件,具有全天时,全天候成像能力,在军事和民用方面得到了广泛的应用。但是SAR的相干成像机制导致成像结果含有很强的相干斑噪声，使得目标的识别和特征提取变得相当困难，因此有必要展开SAR图像的相干斑抑制工作。

SAR图像抑斑主要有两个目标，其一是有效地抑制图像同质区域的相干斑噪声，其二是尽可能的保留图像中场景的边缘，纹理，点目标等细节信息。

最初的SAR图像抑斑采用空域的处理方法，在上世纪80年代得到了极大的发展，涌现了很多优秀滤波算法，例如，Lee滤波器，Kuan滤波器，Frost滤波器及其增强版本等。这些降班方法都是利用图像的局部图像块信息进行去噪，其优点是图像的辐射特性保持很好，但是大小固定的滤波窗口容易对图像的细节和纹理产生过平滑现象。

CA Deledalle等人在2009年提出的“基于图像块最大似然概率迭代加权的SAR去噪”简记为PPB方法，是迄今为止最优秀的去噪方法之一。该方法能够在抑制相干斑的同时较好地保持图像的纹理和细节，但是具有非常明显的辐射特性损失。

发明内容

本发明的目的在于针对上述现有的技术问题，提出一种基于稀疏域噪声分布约束的SAR图像抑斑方法，以在保持图像纹理和细节的同时，减小辐射特性损失。

为实现上述目的，本发明的技术方案包括如下步骤：

(1)将原始SAR图像进行对数变换，得到变换后的图像m，并将m划分为大小相同的重叠的子图像块M＝{m₁,m₂,…,m_i,…,m_L},i＝1,2,...,L，其中L为划分的子图像块总数；

(2)用K均值聚类方法将子图像块M聚为K类{S₁,S₂,…,S_q,…,S_K}，q＝1,2,…,K；

(3)用主成分分析的方法计算得到每一个聚类的子字典D_q，并将子字典顺序拼接得到稀疏表示字典D；

(4)计算第i个含噪子图m_i的稀疏表示系数α_i，α_i＝D^Tm_i，并用稀疏域非局部均值的方法计算m_i对应的无噪子图稀疏表示系数的估计值得到稀疏域噪声

(5)利用变换后的图像m，并根据稀疏域噪声ε_i的稀疏特性，构建如下目标函数：

${\widehat{\mathrm{\α}}}_i=\arg \underset{{\mathrm{\α}}_i}{\min }\{{\left|\right|R_{i}m-{\mathrm{D\α}}_i\left|\right|}_2^2+{\mathrm{\κ}}_i{\left|\right|{\mathrm{\α}}_i-{\stackrel{\‾}{\mathrm{\α}}}_i\left|\right|}_1\},i=\mathrm{1,2},...,L,$

其中R_i是抽取图像m中第i个图像块的矩阵，κ_i是稀疏约束项的正则参数，||·||₂表示向量的2范数，||·||₁表示向量的1范数，表示降噪稀疏表示系数；

(6)求解目标函数，得到去噪后的SAR图像：

(6a)初始化目标函数的正则参数κ_i＞0，

(6b)通过迭代收缩算法求解目标函数，得到每一个子图像块最终的降噪稀疏表示系数

(6c)对最终的降噪稀疏表示系数进行反稀疏变换，得到降噪后的子图像块 ${\hat{m}}_i=D{\widehat{\mathrm{\α}}}_i;$

(6d)拼接降噪后的子图像块得到完整的图像，并对该图像进行反对数变换，得到去噪后的SAR图像。

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

第一，抑斑能力强：以等效视数ENL为评价指标，本发明的SAR图像同质区域抑斑能力明显好于传统的经典Frost滤波器。

第二，辐射特性保持度好：以比值图像的均值为评价指标，本发明的SAR图像的辐射特性保持能力优于传统的Frost滤波方法和最近提出的经典迭代版PPB滤波方法。

第三，图像的纹理和细节保持较好：以比值图像的视觉效果展示，用本发明的方法抑斑后，SAR图像的纹理和细节的保持能力明显好于传统的Frost滤波方法，并且优于迭代版PPB滤波方法。

附图说明

图1是本发明的总流程图；

图2是原始的SAR图像；

图3是采用现有Frost方法对图2抑斑后的图像；

图4是采用现有Frost方法对图2抑斑后的图像的比值图；

图5是采用现有迭代版PPB方法对图2抑斑后的图像；

图6是采用现有迭代版PPB方法对图2抑斑后的图像的比值图；

图7是采用本发明方法对图2抑斑后的图像；

图8是采用本发明方法对图2抑斑后的图像的比值图。

具体实施方式

下面结合附图对本发明做进一步的详细描述。

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

步骤1，对原始的SAR图像进行对数变换并划分为重叠的子图像块。

将原始SAR图像先进行对数变换，得到变换后的图像m；再将变换后的图像m划分为L个大小相同的重叠的子图像块M＝{m₁,m₂,…,m_i,…,m_L}，m_i表示i个子图像块，i＝1,2,...,L，各子图像块的大小均为5×5，子图像块的像素总数为T＝25。

步骤2，对子图像块M进行K均值聚类。

(2a)先对子图像块M进行高通滤波处理，得到滤波后的子图像块为其中是第i个子图m_i滤波后的子图像块,i＝1,2,...,L；

(2b)随机选取滤波后的子图像块M^h中的K个子图作为初始聚类的几何中心，其中K为聚类总数，大小为K＝70；

(2c)计算滤波后的子图像块M^h中的每一个子图到所有几何中心的欧氏距离，并将该子图归入欧氏距离最小的聚类中；

(2d)计算每一个聚类里面所有图像块的均值，得到新的几何中心；

(2e)重复(2c)-(2d)，经过12次迭代，将滤波后的子图像块M^h聚为K类{C₁,C₂,…,C_q,…,C_K}，q＝1,2,…,K；

(2f)根据子图像块M滤波前后的一一对应关系，将相应的子图像块M直接划分为K类{S₁,S₂,…,S_q,…,S_K}，其中S_q表示第q个聚类，q＝1,2,…,K。

步骤3，用主成分分析的方法构造稀疏表示字典D。

(3a)计算子图像块M的第q个聚类S_q的协方差矩阵Ψ_q＝S_q*(S_q)^T，q＝1,2,…,K，

其中T表示矩阵转置，K表示聚类总数；

(3b)对协方差矩阵Ψ_q进行奇异值分解，将其分解为正交矩阵D_q，特征值矩阵Λ_q及正交矩阵V_q，即Ψ_q＝D_qΛ_qV_q ^H，其中H表示共轭转置，取正交矩阵D_q为聚类S_q的子字典；

(3c)将子字典D_q顺序拼接，得到稀疏表示字典D，q＝1,2,…,K。

步骤4，计算各子图m_i对应的无噪子图稀疏表示系数的估计值。

(4a)先计算子图m_i与其周围大小为28×28的正方形邻域内所有子图的欧氏距离，并保留欧氏距离最小的P个子图，P＝13，记为m_ij是子图m_i的第j个相似子图，j＝1,2,…,P，其中是包含于符号，i＝1,2,...,L，L为子图像块总数；

(4b)计算子图m_i的第j个相似子图m_ij在字典D下的稀疏表示系数，α_ij＝D^Tm_ij；

(4c)计算子图m_i对应的无噪子图稀疏表示系数的估计值其中ω_ij表示加权系数， ${\mathrm{\ω}}_{\mathrm{ij}}={\left(2^T\exp \right[{\mathrm{\Σ}}_{t=1}^T(x_t+y_t)]/{\mathrm{\Π}}_{t=1}^T[\exp \left({2x}_t\right)+\exp \left({2y}_t\right)\left]\right)}^{1/\left(\mathrm{Th}\right)},$ W是归一化因子，x_t为子图m_i的第t个像素值，y_t为子图m_i的第j个相似子图m_ij的第t个像素值，T为子图的像素总数，h为衰减参数，h的大小设定为图像m噪声方差的12～15倍。

步骤5，构造目标函数。

利用对数变换后的图像m，稀疏表示字典D，并根据稀疏域噪声i＝1,2,...,L的稀疏特性，构建如下目标函数：

${\widehat{\mathrm{\α}}}_i=\arg \underset{{\mathrm{\α}}_i}{\min }\{{\left|\right|R_{i}m-{\mathrm{D\α}}_i\left|\right|}_2^2+{\mathrm{\κ}}_i{\left|\right|{\mathrm{\α}}_i-{\stackrel{\‾}{\mathrm{\α}}}_i\left|\right|}_1\},i=\mathrm{1,2},...,L,$

其中R_i是抽取图像m中第i个图像块的矩阵，κ_i是稀疏约束项的正则参数，||·||₂表示向量的2范数，||·||₁表示向量的1范数，表示降噪稀疏表示系数。

步骤6，求解目标函数，得到各子图像块m_i的最终降噪稀疏表示系数。

(6a)初始化目标函数的正则参数κ_i＞0，设总迭代次数为N＝15，迭代次数初始值为n＝1；

(6b)用K均值聚类方法更新子图像块M的聚类，记更新后的聚类为其中表示第n次迭代中的第q个聚类，q＝1,2,…,K，K表示聚类总数；

(6c)用主成分分析方法更新字典D，记第n次更新后的字典为D⁽ⁿ⁾；

(6d)用稀疏表示系数估计方法更新第i个子图m_i对应的无噪子图稀疏表示系数的估计值记第n次更新后的估计值为i＝1,2,...,L，L表示子图像块的总数；

(6e)用软阈值收缩算子计算并更新降噪稀疏表示系数得到第n次更新后的第i个降噪稀疏表示系数为其中‘soft’表示表示经典的软阈值收缩算子，表示第i个子图m_i在字典D⁽ⁿ⁾下的稀疏表示系数，κ_i是目标函数稀疏约束项的正则参数；

(6f)用反稀疏变换更新降噪后的子图像块记第n次更新后的第i个降噪子图像块为其中D⁽ⁿ⁾为第n次更新后的字典；

(6g)用第n次更新后的第i个降噪子图像块替换子图像块M中第i个子图像块m_i；

(6h)判断是否满足总迭代次数N，若不满足，将迭代次数n加1，重复(6b)-(6f)继续迭代；若满足总迭代次数，则用最后一次迭代得到的降噪稀疏表示系数作为最终的降噪稀疏表示系数

步骤7，获得抑斑后的SAR图像。

(7a)对最终的降噪稀疏表示系数进行反稀疏变换，得到降噪后的子图像块i＝1,2,...,L；

(7b)将降噪后的各个子图像块进行拼接，得到完整的图像，并对该图像进行反对数变换，得到去噪后的SAR图像。

本发明的仿真效果可以通过下述仿真实验加以说明

1.仿真条件

运行平台配置：

CPU:Inter(R)Core(TM)i56503.20GHz；

内存：8GB(三星DDR31600MHz)；

操作系统:Windows7旗舰版64位SP1；

仿真软件:MATLAB R(2011b)。

2.仿真内容与结果

仿真1，用现有Frost抑斑方法对图2进行抑斑，并获取抑斑前后的比值图，结果如图3和图4。其中：

图3(a)为对图2(a)所示的梯田图像抑斑后的结果，

图3(b)为对图2(b)所示的港口图像抑斑后结果，

图4(a)为图2(a)与图3(a)各像素点的比值图，

图4(b)为图2(b)与图3(b)各像素点的比值图。

仿真2，用现有PPB抑斑方法对图2进行抑斑，并获取抑斑前后的比值图，结果如图5和图6。其中：

图5(a)为对图2(a)所示的梯田图像抑斑后的结果，

图5(b)为对图2(b)所示的港口图像抑斑后结果，

图6(a)为图2(a)与图5(a)各像素点的比值图，

图6(b)为图2(b)与图5(b)各像素点的比值图。

仿真3，用本发明的抑斑方法对图2进行抑斑，并获取抑斑前后的比值图，结果如图7和图8。其中：

图7(a)为对图2(a)所示的梯田图像抑斑后的结果，

图7(b)为对图2(b)所示的港口图像抑斑后结果，

图8(a)为图2(a)与图7(a)各像素点的比值图，

图8(b)为图2(b)与图7(b)各像素点的比值图。

3.仿真结果分析

(3.1)计算仿真1，仿真2和仿真3三种不同抑斑方法得到的各个抑斑后的图像的同质区域的等效视数，结果如表1所示。表中等效视数值越大说明抑斑能力越强。

表1三种不同方法得到抑斑后的图像的同质区域的等效视数

	Frost	迭代版PPB	本发明方法
				梯田	239.55	798.00	668.67

港口

43.79

50.10

67.91

从表1中数据可以看出，本发明的方法在抑斑能力上明显优于现有的Frost方法，而且对于具有丰富细节特性的港口图像，本发明的抑斑能力优于现有的迭代版PPB方法。

(3.2)计算仿真1，仿真2和仿真3三种不同抑斑方法得到的各个抑斑后的图像的比值图的均值，结果如表2所示。表中均值越接近1说明抑斑前后图像辐射特性保持越好。

表2三种不同方法得到抑斑后的图像的比值图的均值

	Frost	迭代版PPB	本发明方法
				梯田	0.9974	0.9915	1.0008
港口	0.9851	0.9609	0.9993

从表2中数据可以看出，本发明的SAR图像抑斑方法相比其他两种抑斑方法，比值图像的均值更接近于1，说明在抑斑的过程中原始SAR图像没有明显增加或损失辐射特性，因此本发明的抑斑方法，辐射特性保持能力要优于Frost和PPB抑斑方法。

(3.3)观察不同抑斑方法得到的图像的比值图的视觉效果，发现现有的Frost抑斑方法得到的比值图像含有很强的纹理结构残留，说明在抑斑过程中，原始SAR图像的细节纹理信息被过平滑掉，有信息的损失；此外现有的迭代版PPB抑斑方法得到的比值图像也含有一定程度的纹理结构残留，同样有信息损失；而本发明方法几乎不存在纹理结构残留，因此在边缘和纹理信息保持方面要优于传统的经典Frost方法和最近提出的PPB方法。

Claims (5)

1.一种基于稀疏域噪声分布约束的SAR图像抑斑方法，包括如下步骤：

(1)将原始SAR图像进行对数变换，得到变换后的图像m，并将m划分为大小相同的重叠的子图像块M＝{m₁,m₂,…,m_i,…,m_L},i＝1,2,...,L，其中L为划分的子图像块总数；

(2)用K均值聚类方法将子图像块M聚为K类{S₁,S₂,…,S_q,…,S_K}，q＝1,2,…,K；

(3)用主成分分析的方法计算得到每一个聚类的子字典D_q，并将子字典顺序拼接得到稀疏表示字典D；

(4)计算第i个含噪子图m_i的稀疏表示系数α_i，α_i＝D^Tm_i，并用稀疏域非局部均值方法计算m_i对应的无噪子图稀疏表示系数的估计值，得到稀疏域噪声i＝1,2,...,L；

(5)利用变换后的图像m，并根据稀疏域噪声ε_i的稀疏特性，构建如下目标函数：

${\widehat{\mathrm{\α}}}_i=\arg \underset{{\mathrm{\α}}_i}{\min }\{{\left|\right|R_{i}m-{\mathrm{D\α}}_i\left|\right|}_2^2+{\mathrm{\κ}}_i{\left|\right|{\mathrm{\α}}_i-{\stackrel{\‾}{\mathrm{\α}}}_i\left|\right|}_1\},i=\mathrm{1,2},...,L,$

其中R_i是抽取图像m中第i个图像块的矩阵，κ_i是稀疏约束项的正则参数，||·||₂表示向量的2范数，||·||₁表示向量的1范数，表示降噪稀疏表示系数；

(6)求解目标函数，得到去噪后的SAR图像：

(6a)初始化目标函数的正则参数κ_i＞0，

(6b)通过迭代收缩算法求解目标函数，得到每一个子图像块最终的降噪稀疏表示系数；

(6c)对最终的降噪稀疏表示系数进行反稀疏变换，得到降噪后的子图像块 ${\hat{m}}_i=D{\widehat{\mathrm{\α}}}_i;$

(6d)拼接降噪后的子图像块得到完整的图像，并对该图像进行反对数变换，得到去噪后的SAR图像。

2.根据权利要求1所述的基于稀疏域噪声分布约束的SAR图像抑斑方法，其中第(2)步所述的用K均值聚类方法将子图像块M聚为K类，按如下步骤进行：

(2a)对子图像块M进行高通滤波处理，记滤波后的子图像块为其中是第i个子图m_i滤波后的子图像块,i＝1,2,...,L；

(2b)用K均值聚类的方法对滤波后的子图像块M^h进行聚类：

(2b1)随机选取滤波后的子图像块M^h中的K个子图作为初始聚类的几何中心,其中K为聚类总数；

(2b2)计算滤波后的子图像块M^h中的每一个子图到所有几何中心的欧氏距离，并将该子图归入欧氏距离最小的聚类中；

(2b3)计算每一个聚类里面所有图像块的均值，得到新的几何中心；

(2b4)重复(2b2)-(2b3)，经过12次迭代，将滤波后的子图像块M^h聚为K类{C₁,C₂,…,C_q,…,C_K}，q＝1,2,…,K；

(2c)根据子图像块M滤波前后的一一对应关系，将相应的子图像块M直接划分为K类{S₁,S₂,…,S_q,…,S_K}，其中S_q表示第q个聚类，q＝1,2,…,K。

3.根据权利要求1所述的基于稀疏域噪声分布约束的SAR图像抑斑方法，其中步骤(3)所述的用主成分分析的方法计算得到每一个聚类的子字典D_q，按如下步骤进行：

(3a)计算子图像块M的第q个聚类S_q的协方差矩阵Ψ_q＝S_q*(S_q)^T，q＝1,2,…,K，其中T表示矩阵转置，K表示聚类总数；

(3b)对协方差矩阵Ψ_q进行奇异值分解，将其分解为正交矩阵D_q，特征值矩阵Λ_q及正交矩阵V_q，即其中H表示共轭转置，取正交矩阵D_q为聚类S_q的子字典。

4.根据权利要求1所述的基于稀疏域噪声分布约束的SAR图像抑斑方法，其中步骤(4)所述的用稀疏域非局部均值的方法计算子图m_i对应的无噪子图稀疏表示系数的估计值，按如下步骤进行：

(4a)计算子图m_i与其周围大小为w×w的正方形邻域内所有子图的欧氏距离，并保留欧氏距离最小的P个子图，记为，m_ij是子图m_i的第j个相似子图，j＝1,2,…,P，其中w是正方形邻域的边长，是包含于符号；

(4b)计算子图m_i的第j个相似子图m_ij在字典D下的稀疏表示系数，α_ij＝D^Tm_ij；

(4c)计算子图m_i对应的无噪子图稀疏表示系数的估计值其中ω_ij表示加权系数， ${\mathrm{\ω}}_{\mathrm{ij}}={\left(2^T\exp \right[{\mathrm{\Σ}}_{t=1}^T(x_t+y_t)]/{\mathrm{\Π}}_{t=1}^T[\exp \left({2x}_t\right)+\exp \left({2y}_t\right)\left]\right)}^{1/\left(\mathrm{Th}\right)},$ W是归一化因子，x_t为子图m_i的第t个像素值，y_t为子图m_i的第j个相似子图m_ij的第t个像素值，T为子图的像素总数，h为衰减参数，h的大小设定为图像m噪声方差的12～15倍。

5.根据权利要求1所述的基于稀疏域噪声分布约束的SAR图像抑斑方法，其中步骤(6)中的(6b)所述的通过迭代收缩算法求解目标函数，得到每一个子图像块最终的降噪稀疏表示系数按如下步骤进行：

(6b1)总迭代次数设为N，迭代次数初始值设为n＝1；

(6b2)用K均值聚类方法更新子图像块M的聚类，记更新后的聚类为其中表示第n次迭代中的第q个聚类，q＝1,2,…,K，K表示聚类总数；

(6b3)用主成分分析方法更新字典D，记第n次更新后的字典为D⁽ⁿ⁾；

(6b4)用稀疏域非局部均值方法更新第i个子图mi对应的无噪子图稀疏表示系数的估计值记第n次更新后的估计值为i＝1,2,...,L，L表示子图像块的总数；

(6b5)用软阈值收缩算子计算并更新降噪稀疏表示系数得到第n次更新后的第i个降噪稀疏表示系数为其中‘soft’表示表示经典的软阈值收缩算子，表示第i个子图m_i在字典D⁽ⁿ⁾下的稀疏表示系数，κ_i是目标函数稀疏约束项的正则参数；

(6b6)用反稀疏变换更新降噪后的子图像块记第n次更新后的第i个降噪子图像块为其中D⁽ⁿ⁾为第n次更新后的字典；

(6b7)用第n次更新后的第i个降噪子图像块替换子图像块M中第i个子图像块mi；

(6b8)判断是否满足总迭代次数N，若不满足，将迭代次数n加1，重复(6b2)-(6b6)继续迭代；若满足总迭代次数，则用最后一次迭代得到的降噪稀疏表示系数作为最终的降噪稀疏表示系数