CN104537624B - 基于ssim校正聚类稀疏表示的sar图像降斑方法 - Google Patents

Abstract

本发明公开了一种基于SSIM校正聚类的稀疏表示SAR图像降斑方法，主要解决现有技术中不能保留更多的点目标和纹理细节的问题。其实现步骤为：(1)输入图像；(2)估计输入SAR图像中噪声的方差；(3)对输入图像进行分类；(4)对分类后图像进行初始化字典；(5)更新分类后图像的稀疏表示矩阵和字典；(6)利用更新后的稀疏表示矩阵和字典重构图像，得到降斑后的SAR图像。本发明不仅能够提高SAR图像中斑点噪声的抑制效果，而且能很好保留图像中点目标和纹理细节信息，可用于SAR图像分割和目标识别的预处理。

Description

基于SSIM校正聚类稀疏表示的SAR图像降斑方法

技术领域

本发明属于图像处理技术领域，更进一步涉及合成孔径雷达SAR图像的降斑方法，可用于SAR图像分割和目标识别的预处理。

背景技术

合成孔径雷达SAR具有全天候、全天时、高分辨率和穿透力强等特点，因此SAR在军事方面和民用方面都得到了广泛的应用。由于SAR系统是微波相干成像，因此SAR图像在获取的过程中容易受到斑点噪声的影响。斑点噪声的存在大大降低了SAR图像的分辨率，影响了SAR图像后续的处理与解译，如何有效地抑制SAR图像中的斑点噪声变得尤为重要。

近年来，稀疏表示被广泛的应用到图像去噪中。然而，标准的稀疏表示算法假定稀疏矩阵中非零系数的出现是随机的，并没有考虑到图像块之间的相似性。因此，一些学者提出了将图像的几何结构相似性与稀疏表示相结合的算法。其中，董伟生等人提出的基于聚类的稀疏表示CSR算法认为图像块之间的稀疏表示系数的分布不是随机的，而是在某种程度上与图像信号的自相似性有关。

西安电子科技大学在其专利申请“基于稀疏域噪声分布约束的SAR图像抑斑方法”(专利申请号：201410229349.2，公开号：CN103971346A)中提出了一种基于稀疏域噪声分布约束的SAR图像抑斑方法。该方法首先对SAR图像进行对数变换，利用K均值聚类的方法对图像块进行聚类。针对每类图像块，采用主成分分析的方法求出对应的子字典。然后，利用非局部均值算法和软阈值收缩算法得到最终的稀疏表示系数，通过反稀疏变换和反对数变换得到最终的SAR图像。该方法虽然解决了现有降斑方法不能同时保持图像纹理细节与辐射特性的问题，能对图像块进行聚类，但是不能保证所有图像块都得到了最好的划分。

西安电子科技大学在其专利申请“基于SAR图像局部统计特性的KSVD相干斑抑制方法”(专利申请号：201110318457.3，公开号：CN102509263A)中提出了一种基于SAR图像局部统计特性的KSVD相干斑抑制方法。该方法从图像中随机抽取训练样本，利用SAR_KSVD算法对字典进行训练，得到最终的训练字典。然后，利用SAR_OMP算法进行稀疏编码的过程，得到稀疏编码系数。最后，根据冗余稀疏表示图像噪声抑制理论，得到最终的降斑图像。该方法能够很好的解决现有相干斑抑制方法中的边缘和纹理等细节信息被模糊的问题。但是由于该方法在提取训练图像块时并没有将图像块进行详细的划分，而是通过随机抽取的方法进行选取，导致训练得到的字典不能充分的包含图像中的纹理细节信息。

发明内容

本发明针对上述现有技术的不足，提出一种基于SSIM校正聚类稀疏表示的SAR图像降斑方法，以在保留点目标和纹理细节信息的前提下，提高SAR图像中的斑点噪声抑制效果。

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

(1)输入任选的一幅合成孔径雷达SAR图像；

(2)估计输入SAR图像中噪声的方差：

(2a)对SAR图像进行方向波变换，得到低频子带系数矩阵和高频子带系数矩阵；

(2b)按照下式，采用基于方向波域的斑点噪声估计方法，估计SAR图像中噪声的方差：

其中，σ表示SAR图像中噪声的标准差，C表示噪声的归一化标准差，Ψ表示控制参数，μ_y表示SAR图像的均值，y表示SAR图像矩阵，表示SAR图像的高频子带系数矩阵的方差，D_y表示SAR图像的高频子带系数矩阵；

(3)图像块分类：

(3a)采用KSVD算法，对SAR图像进行预滤波，得到滤波后的SAR图像y′；

(3b)以滤波后的SAR图像y′中任意像素为中心，抽取大小为8×8的图像块，得到SAR图像像素图像块集合B₁，以SAR图像的低频子带系数矩阵中任意系数为中心，抽取大小为8×8的图像块，得到SAR图像低频子带图像块集合B₂，合并SAR图像像素图像块集合B₁和SAR图像低频子带图像块集合B₂，得到总图像块集合B；

(3c)计算总图像块集合B中每个图像块的方差，根据得到的图像块的方差大小，将总图像块集合B分为光滑图像块集合E和非光滑图像块集合F；

(4)初始化字典：

(4a)采用基于SSIM校对的K均值聚类算法，对非光滑图像块集合F中的所有图像块进行聚类，得到N类相似图像块集合S_k，1≤k≤N，N为K均值聚类所得到的类别总数；

(4b)采用主成分分析方法，求解相似图像块集合S_k对应的64个特征向量，将这64个特征向量按列排列，得到S_k对应的字典D_k，1≤k≤N；

(4c)将字典D_k作为相似图像块集合S_k的初始字典，将大小为64×256的离散余弦变换DCT字典作为光滑图像块集合E的初始字典；

(5)根据基于结构聚类的稀疏表示方法和字典更新方法，对相似图像块集合S_k和光滑图像块集合E进行稀疏表示矩阵更新和字典更新，得到最终的稀疏表示矩阵和最终的字典

(6)根据上述得到的总图像块集合B、最终的稀疏表示矩阵和最终的字典利用下式得到降斑后的SAR图像x：

其中，λ表示拉格朗日因子，I表示与SAR图像大小一样的单位矩阵，R表示总图像块集合B中任意图像块矩阵，R^T表示图像块矩阵R的转置。

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

第一，由于本发明采用了方向波变换来获取不同尺度下的图像，克服了现有技术中不能充分结合不同尺度下图像块的自相似性而导致训练得到的字典包含信息量少的缺点，使得本发明充分的利用了图像的细节信息。

第二，由于本发明采用了基于方向波域的斑点噪声估计方法，克服了现有技术中不能更精确的估计SAR图像中斑点噪声方差而导致降斑结果不理想的问题，使得本发明能够更好的去除图像中的斑点噪声。

第三，由于本发明采用了基于SSIM校正的K均值聚类方法，克服了现有技术中利用K均值聚类方法对图像块进行聚类时而导致的图像块划分不准确的问题，使得本发明能够更好的保留图像中的细节信息。

附图说明

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

图2为本发明与现有技术对haigang合成孔径雷达SAR图像的降斑效果对比图；

图3为本发明与现有技术对field合成孔径雷达SAR图像的降斑效果对比图。

具体实施方式

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

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

步骤1：输入图像。

输入任选的一幅合成孔径雷达SAR图像。本发明实施例中使用的合成孔径雷达SAR图像分别如附图2(a)和附图3(a)所示。其中，附图2(a)是合成孔径雷达SAR图像海港图haigang，大小为256×256，附图3(a)是合成孔径雷达SAR图像田野图field，大小为512×512。

步骤2：估计输入SAR图像中的噪声方差。

(2a)对SAR图像进行方向波变换，得到低频子带系数矩阵和高频子带系数矩阵；

(2b)按照下式，采用基于方向波域的斑点噪声估计方法，估计SAR图像中噪声的方差：

其中，σ表示SAR图像中噪声的标准差，C表示噪声的归一化标准差，Ψ表示控制参数，μ_y表示SAR图像的均值，y表示SAR图像矩阵，表示SAR图像的高频子带系数矩阵的方差，D_y表示SAR图像的高频子带系数矩阵。

步骤3：图像块分类。

(3a)采用KSVD算法，对SAR图像进行预滤波，得到滤波后的SAR图像y′，其中，KSVD算法的具体步骤如下：

(3a1)以SAR图像y中任意像素为中心，抽取大小为8×8的图像块，得到的SAR图像像素图像块集合M，初始化字典D为离散余弦变换DCT字典，其大小为64×256；

(3a2)按照下式，更新稀疏表示矩阵，得到更新后的稀疏表示矩阵a′：

其中，a表示稀疏表示矩阵，μ表示控制参数，表示取二范数的平方操作，||·||₀表示取零范数操作；

(3a3)按照下式更新字典，得到更新后的字典D′：

其中，λ表示拉格朗日因子，||·||₁表示取一范数操作；

(3a4)重复执行步骤(3a2)和(3a3)共10次，得到KSVD算法的目标稀疏表示矩阵a″和目标字典D″；

(3a5)按照下式重构得到KSVD算法滤波后的SAR图像y′：

其中，I表示与SAR图像大小一样的单位矩阵，R_M表示SAR图像像素图像块集合M中任意图像块矩阵，表示R_M的转置；

(3b)以滤波后的SAR图像y′中任意像素为中心，抽取大小为8×8的图像块，得到SAR图像像素图像块集合B₁，以SAR图像的低频子带系数矩阵中任意系数为中心，抽取大小为8×8的图像块，得到SAR图像低频子带图像块集合B₂，合并SAR图像像素图像块集合B₁和SAR图像低频子带图像块集合B₂，得到总图像块集合B；

(3c)计算总图像块集合B中每个图像块的方差，并将图像块的方差与设定的阈值进行比较，若图像块的方差大于设定的阈值τ，则将该图像块划分到非光滑图像块集合F中，否则，将该图像块划分到光滑图像块集合E中，τ为最大图像块方差的三分之一。

步骤4：初始化字典。

(4a)采用基于SSIM校对的K均值聚类算法，对非光滑图像块集合F中的所有图像块进行聚类：

(4a1)利用K均值聚类算法，对非光滑图像块集合F中的所有图像块进行N类聚类，得到每个图像块的类别和每个类别的聚类中心；

(4a2)计算每个图像块及其所属类别的聚类中心的结构相似度SSIM值，判断该SSIM值是否小于0.85，若是，则执行(4a3)，否则，执行(4a4)；

(4a3)计算该图像块与其他类别的聚类中心的SSIM值，将该图像块划分到其最大的SSIM值所对应的类别中；

(4a4)将所属类别为k的所有图像块聚集起来，得到相似图像块集合S_k，1≤k≤N，N为K均值聚类所得到的类别总数；

(4b)采用主成分分析方法，求解相似图像块集合S_k对应的64个特征向量，将这64个特征向量按列排列，得到S_k对应的字典D_k；

(4c)将字典D_k作为相似图像块集合S_k的初始字典，将大小为64×256的离散余弦变换DCT字典D_l作为光滑图像块集合E的初始字典。

步骤5：对相似图像块集合S_k和光滑图像块集合E进行稀疏表示矩阵更新和字典更新。

(5a)根据基于结构聚类的稀疏表示方法和字典更新方法，按照下式，更新稀疏表示矩阵和聚类中心所对应的稀疏表示矩阵，得到更新后的稀疏表示矩阵和聚类中心所对应的稀疏表示矩阵

其中，u_k表示类别为k的聚类中心向量，k＝1,2,...,N，λ₁和λ₂表示两个不同数值的常数，N表示K均值聚类的类别总数，C_k表示类别为k的稀疏表示向量矩阵，i表示稀疏表示向量的列数，a_i表示列数为i的稀疏表示向量，β_k表示类别为k的聚类中心对应的稀疏表示向量；，

(5b)利用步骤(3a)中KSVD算法的字典更新方法，更新字典D_k和D_l；

(5c)重复执行步骤(5a)和(5b)共3次，得到最终的稀疏表示矩阵和最终的字典

步骤6：根据上述得到的总图像块集合B、最终的稀疏表示矩阵和最终的字典利用下式得到降斑后的SAR图像x：

其中，R表示总图像块集合B中任意图像块矩阵，R^T表示图像块矩阵R的转置。

下面结合附图2和附图3的仿真图对本发明的效果做进一步说明。

1.仿真实验条件：

本发明的硬件测试平台是：处理器为Inter Core 2 Duo CPU E8200，主频为2.67GHz，内存2GB，软件平台为：Windows 7旗舰版32位操作系统和Matlab R2012b。本发明的输入图像分别为合成孔径雷达SAR图像海港图haigang和合成孔径雷达SAR图像田野图field，大小分别为256×256和512×512，格式都为PNG。

仿真使用的方法是：本发明和如下的两种现有方法：

第一种现有方法：S.Foucher等人在文献“Multiscale MAP filtering of SARimages[J].IEEE Trans.on Image Processing,2001,10(1):49-60”中提到的Gamma-MAP滤波方法。

第二种现有方法：Deledalle等人在文献“Iterative weighted maximumlikelihood denoising with probabilistic patch-based weights[J].IEEETransactions on Image Processing,2009,18(12):2661-2672.”中提出的合成孔径雷达SAR图像去斑方法，简称PPB方法。

2.仿真内容与结果：

仿真1，用本发明和现有两种方法对合成孔径雷达SAR海港图像haigang进行去斑处理。结果如图2，其中，

图2(a)为大小为256×256的原始合成孔径雷达SAR海港图像haigang；

图2(b)为用第一种现有方法对图2(a)进行去斑的结果图；

图2(c)为用第二种现有方法对图2(a)进行去斑的结果图；

图2(d)为用本发明对图2(a)进行去斑的结果图。

仿真2，用本发明和现有两种方法对合成孔径雷达SAR田野图像field进行去斑处理。结果如图3，其中：

图3(a)为大小为512×512的原始合成孔径雷达SAR田野图像field；

图3(b)为用第一种现有方法对图3(a)进行去斑的结果图；

图3(c)为用第二种现有方法对图3(a)进行去斑的结果图；

图3(d)为用本发明对图3(a)进行去斑的结果图。

3.仿真结果分析：

从图2(b)和图3(b)中可以看到，采用第一种现有方法得到的结果图中仍然存在大量的斑点噪声，且图像中纹理细节区域比较模糊。

从图2(c)和图3(c)中可以看出，采用第二种现有方法得到的结果图中斑点噪声得到了很好的抑制，均匀区域非常光滑，但是不能很好的保留图像中点目标和纹理细节。

从图2(d)和图3(d)中可以看出，本发明的方法不仅能够提高SAR图像中斑点噪声的抑制效果，且能够很好的保留图像中的点目标和纹理细节信息。

采用本发明和现有两种方法分别对合成孔径雷达SAR海港图像haigang和合成孔径雷达SAR田野图像field进行降斑的性能指标如表1所示：

表1对海港图像haigang和田野图像field去斑结果的评价指标

表1中均值mean用来衡量图像灰度值的保持能力，降斑后图像的均值越接近原始图像的均值越好；标准差std用来衡量降斑方法的平滑能力，标准差越小，表示其平滑能力越强；等效视数ENL是针对合成孔径雷达SAR图像设计的度量标准，等效视数越大，降斑效果越理想。

从表1中可以看出，本发明去斑结果的均值比较接近于输入图像，表明本发明对合成孔径雷达SAR图像的辐射特性保持较好。本发明去斑结果的标准差最小，表明本发明的平滑能力最好。本发明去斑结果中匀质区域的等效视数ENL最高，表明本发明去斑效果最好。

Claims (3)

1.一种基于SSIM校正聚类稀疏表示的SAR图像降斑方法，包括如下步骤：

(1)输入任选的一幅合成孔径雷达SAR图像；

(2)估计输入SAR图像中噪声的方差：

(2a)对SAR图像进行方向波变换，得到低频子带系数矩阵和高频子带系数矩阵；

(2b)按照下式，采用基于方向波域的斑点噪声估计方法，估计SAR图像中噪声的方差：

${\mathrm{\σ}}^2=\frac{C^2({\mathrm{\Ψ\μ}}_y^2+{\mathrm{\σ}}_{D_y}^2)}{(1+C^2)},$

其中，σ表示SAR图像中噪声的标准差，C表示噪声的归一化标准差，Ψ表示控制参数，μ_y表示SAR图像的均值，y表示SAR图像矩阵，表示SAR图像的高频子带系数矩阵的方差，D_y表示SAR图像的高频子带系数矩阵；

(3)图像块分类：

(3a)采用KSVD算法，对SAR图像进行预滤波，得到滤波后的SAR图像y′；

(3b)以滤波后的SAR图像y′中任意像素为中心，抽取大小为8×8的图像块，得到滤波后的SAR图像y′像素图像块集合B₁，以SAR图像的低频子带系数矩阵中任意系数为中心，抽取大小为8×8的图像块，得到SAR图像低频子带图像块集合B₂，合并SAR图像像素图像块集合B₁和SAR图像低频子带图像块集合B₂，得到总图像块集合B；

(3c)计算总图像块集合B中每个图像块的方差，根据得到的图像块的方差大小，将总图像块集合B分为光滑图像块集合E和非光滑图像块集合F；

(4)初始化字典：

(4a)采用基于SSIM校对的K均值聚类算法，对非光滑图像块集合F中的所有图像块进行聚类，得到N类相似图像块集合S_k，1≤k≤N，N为K均值聚类所得到的类别总数；

所述的基于SSIM校对的K均值聚类算法，对非光滑图像块集合F中的所有图像块进行聚类，按如下步骤进行：

(4a1)利用K均值聚类算法，对非光滑图像块集合F中的所有图像块进行N类聚类，得到每个图像块的类别和每个类别的聚类中心；

(4a2)计算每个图像块与其所属类别的聚类中心的结构相似度SSIM值，判断该SSIM值是否小于0.85，若是，则执行步骤(4a3)，否则，执行步骤(4a4)；

(4a3)计算该图像块与其他类别的聚类中心的SSIM值，将该图像块划分到其最大的SSIM值所对应的类别中；

(4a4)将所属类别为k的所有图像块聚集起来，得到相似图像块集合S_k，1≤k≤N,N为K均值聚类所得到的类别总数；

(4b)采用主成分分析方法，求解相似图像块集合S_k对应的64个特征向量，将这64个特征向量按列排列，得到S_k对应的字典D_k，1≤k≤N；

(4c)将字典D_k作为相似图像块集合S_k的初始字典，将大小为64×256的离散余弦变换DCT字典D_l作为光滑图像块集合E的初始字典；

(5)根据基于结构聚类的稀疏表示方法和字典更新方法，对相似图像块集合S_k和光滑图像块集合E进行稀疏表示矩阵更新和字典更新，得到最终的稀疏表示矩阵和最终的字典

所述的稀疏表示矩阵更新和字典更新按如下步骤进行：

(5a)按照下式，更新稀疏表示矩阵和聚类中心所对应的稀疏表示矩阵，得到更新后的稀疏表示矩阵和聚类中心所对应的稀疏表示矩阵

$(\tilde{a},\widetilde{\mathrm{\β}})=\underset{a,{\mathrm{\μ}}_k}{\arg min}\frac{1}{2}\left|\right|y^{\′}-Da|{|}_2^2+{\mathrm{\λ}}_1|\left|a\right|{|}_1+{\mathrm{\λ}}_2\underset{k=1}{\overset{N}{\Σ}}\underset{i\∈C_k}{\Σ}\left|\right|a_i-{\mathrm{\β}}_k|{|}_1$

其中，a表示稀疏表示矩阵，μ_k表示类别为k的聚类中心向量，k＝1,2,...,N，λ₁和λ₂表示两个不同数值的常数，N表示K均值聚类的类别总数，C_k表示类别为k的稀疏表示向量矩阵，i表示稀疏表示向量的列数，a_i表示列数为i的稀疏表示向量，β_k表示类别为k的聚类中心对应的稀疏表示向量，表示取二范数的平方操作，||·||₁表示取一范数操作；

(5b)利用步骤(3a)中KSVD算法的字典更新方法，更新字典D_k和D_l；

(5c)重复执行步骤(5a)和步骤(5b)共3次，得到最终的稀疏表示矩阵和最终的字典

(6)根据上述得到的总图像块集合B、最终的稀疏表示矩阵和最终的字典利用下式得到降斑后的SAR图像x：

$x={(\mathrm{\λ}I+{\mathrm{\ΣR}}^{T}R)}^{-1}({\mathrm{\λy}}^{\′}+{\mathrm{\ΣR}}^T\hat{D}\hat{a}),$

其中，λ表示拉格朗日因子，I表示与SAR图像大小一样的单位矩阵，R表示总图像块集合B中任意图像块矩阵，R^T表示图像块矩阵R的转置。

2.根据权利要求1所述基于SSIM校正聚类稀疏表示的SAR图像降斑方法，其特征在于，步骤(3a)所述的采用KSVD算法，对SAR图像进行预滤波，按如下步骤进行：

(3a1)以SAR图像y中任意像素为中心，抽取大小为8×8的图像块，得到的SAR图像像素图像块集合M，初始化字典D为大小是64×256的离散余弦变换DCT字典；

(3a2)按照下式，更新稀疏表示矩阵：

$a^{\′}=\underset{a}{\arg min}\left|\right|Da-y|{|}_2^2+\mathrm{\μ}|\left|a\right|{|}_0,$

其中，a′表示更新后的稀疏表示矩阵，a表示稀疏表示矩阵，μ表示控制参数，表示取二范数的平方操作，||·||₀表示取零范数操作；

(3a3)按照下式，更新字典，得到更新后的字典D′：

$D^{\′}=\underset{D}{\arg \min }\mathrm{\Σ}\underset{a}{min}\left\{\right||Da-y|{|}^2+\mathrm{\λ}\left|\right|a\left|{|}_1\right\},$

其中，||·||₁表示取一范数操作；

(3a4)重复执行步骤(3a2)和(3a3)共10次，得到KSVD算法的目标稀疏表示矩阵a″和目标字典D″；

(3a5)按照下式重构得到KSVD算法滤波后的SAR图像y′：

$y^{\′}={(\mathrm{\λ}I+{\mathrm{\ΣR}}_M^T{R}_M)}^{-1}(\mathrm{\λ}y+{\mathrm{\ΣR}}_M^T{D}^{\′\′}{a}^{\′\′}),$

其中，R_M表示SAR图像像素图像块集合M中任意图像块矩阵，表示图像块矩阵R_M的转置。

3.根据权利要求1所述基于SSIM校正聚类稀疏表示的SAR图像降斑方法，其特征在于，步骤(3c)所述的根据得到的图像块的方差大小，将总图像块集合B分为光滑图像块集合E和非光滑图像块集合F，是将每个图像块的方差与设定的阈值进行比较，若该图像块的方差大于设定的阈值τ，则将该图像块划分到非光滑图像块集合F中，否则，将该图像块划分到光滑图像块集合E中，其中τ的取值为最大图像块方差的三分之一。