基于稀疏域噪声分布约束的sar图像抑斑方法
【专利摘要】本发明公开了一种基于稀疏域噪声分布约束的SAR图像抑斑方法，主要解决现有抑斑方法不能同时保持图像纹理细节及辐射特性的问题。主要步骤为：1.对原始SAR图像进行对数变换，将变换后的图像划分成重叠的大小相同子图像块，并用K均值聚类的方法对图像块进行聚类；2.用主成分分析的方法求出每一个聚类的稀疏表示子字典，并用非局部均值算法计算每个子图像块对应的无噪子图的稀疏表示系数的估计值；3.用软阈值收缩算法求得最终的无噪稀疏表示系数，并通过反稀疏变换得到去噪后的子图像块；4.拼接子图像块并使用反对数变换得到抑斑后的SAR图像。本发明相干斑抑制能力强，图像的纹理细节及辐射特性保持好，可用于目标识别与特征提取。
【专利说明】基于稀疏域噪声分布约束的SAR图像抑斑方法
【技术领域】
[0001]本发明属于图像处理【技术领域】，具体地说是一种合成孔径雷达SAR图像抑斑方法，可用于目标识别与特征提取。
【背景技术】
[0002]合成孔径雷达是一种主动式微波遥感器件，具有全天时，全天候成像能力，在军事和民用方面得到了广泛的应用。但是SAR的相干成像机制导致成像结果含有很强的相干斑噪声，使得目标的识别和特征提取变得相当困难，因此有必要展开SAR图像的相干斑抑制工作。
[0003]SAR图像抑斑主要有两个目标，其一是有效地抑制图像同质区域的相干斑噪声，其二是尽可能的保留图像中场景的边缘，纹理，点目标等细节信息。
[0004]最初的SAR图像抑斑采用空域的处理方法，在上世纪80年代得到了极大的发展，涌现了很多优秀滤波算 法，例如，Lee滤波器，Kuan滤波器，Frost滤波器及其增强版本等。这些降班方法都是利用图像的局部 图像块信息进行去噪，其优点是图像的辐射特性保持很好，但是大小固定的滤波窗口容易对图像的细节和纹理产生过平滑现象。
[0005]CA Deledalle等人在2009年提出的“基于图像块最大似然概率迭代加权的SAR去噪”简记为PPB方法，是迄今为止最优秀的去噪方法之一。该方法能够在抑制相干斑的同时较好地保持图像的纹理和细节，但是具有非常明显的辐射特性损失。
【发明内容】
[0006]本发明的目的在于针对上述现有的技术问题，提出一种基于稀疏域噪声分布约束的SAR图像抑斑方法，以在保持图像纹理和细节的同时，减小辐射特性损失。
[0007]为实现上述目的，本发明的技术方案包括如下步骤:
[0008](I)将原始SAR图像进行对数变换，得到变换后的图像m，并将m划分为大小相同的重叠的子图像块M = Im1, m2，- ,Iiii, - , mL}, i = I, 2,..., L，其中L为划分的子图像块总数；
[0009](2)用K均值聚类方法将子图像块M聚为K类(S1, S2，-,Sq,…，SK}，q= 1，2，…，K ;
[0010](3)用主成分分析的方法计算得到每一个聚类的子字典D,，并将子字典顺序拼接得到稀疏表不字典D ；
[0011](4)计算第i个含噪子图HIi的稀疏表示系数QyCIi = Dt1v并用稀疏域非局部均值的方法计算Hli对应的无噪子图稀疏表示系数的估计值焉4导到稀疏域噪声《( =-- ；
[0012](5)利用变换后的图像m，并根据稀疏域噪声ε ,的稀疏特性，构建如下目标函数:
[0013]at = argmin||/?,-? - Dai|~ + Ki |a(.-Si| J, 1 = 1,2,.,,?L ,
[0014]其中Ri是抽取图像m中第i个图像块的矩阵，K i是稀疏约束项的正则参数，11.112表示向量的2范数，Il.M1表示向量的I范数，爲表示降噪稀疏表示系数；
[0015](6)求解目标函数，得到去噪后的SAR图像:
[0016](6a)初始化目标函数的正则参数KiSO,
[0017](6b)通过迭代收缩算法求解目标函数，得到每一个子图像块最终的降噪稀疏表示
系数忒；
[0018](6c)对最终的降噪稀疏表示系数戎进行反稀疏变换，得到降噪后的子图像块
m, = Dai?
[0019](6d)拼接降噪后的子图像块得到完整的图像，并对该图像进行反对数变换，得到去噪后的SAR图像。
[0020]本发明与现有技术相比具有如下优点:
[0021]第一，抑斑能力强:以等效视数ENL为评价指标，本发明的SAR图像同质区域抑斑能力明显好于传统的经典Frost滤波器。
[0022]第二，辐射特性保持度好:以比值图像的均值为评价指标，本发明的SAR图像的辐射特性保持能力优于传统的Frost滤波方法和最近提出的经典迭代版PPB滤波方法。
[0023]第三，图像的纹理和细节保持较好:以比值图像的视觉效果展示，用本发明的方法抑斑后，SAR图像的纹理和细节的保持能力明显好于传统的Frost滤波方法，并且优于迭代版PPB滤波方法。
【专利附图】
【附图说明】
[0024]图1是本发明的总流程图；
[0025]图2是原始的SAR图像；
[0026]图3是采用现有Frost方法对图2抑斑后的图像；
[0027]图4是采用现有Frost方法对图2抑斑后的图像的比值图；
[0028]图5是采用现有迭代版PPB方法对图2抑斑后的图像；
[0029]图6是采用现有迭代版PPB方法对图2抑斑后的图像的比值图；
[0030]图7是采用本发明方法对图2抑斑后的图像；
[0031]图8是采用本发明方法对图2抑斑后的图像的比值图。
【具体实施方式】
[0032]下面结合附图对本发明做进一步的详细描述。
[0033]参照图1，本发明的具体实施步骤如下:
[0034]步骤1，对原始的SAR图像进行对数变换并划分为重叠的子图像块。
[0035]将原始SAR图像先进行对数变换，得到变换后的图像m ;再将变换后的图像m划分为L个大小相同的重叠的子图像块M = Im1, m2，…，Iiii,…，mj , Hii表示i个子图像块，i =1，2，...，L，各子图像块的大小均为5 X 5，子图像块的像素总数为T = 25。
[0036]步骤2，对子图像块M进行K均值聚类。
[0037](2a)先对子图像 块M进行高通滤波处理，得到滤波后的子图像块为Λ#,( ={μ；',λ4,...--*?..%<},其中》f是第i个子图Hii滤波后的子图像块，i = 1，2，...，L ;
[0038](2b)随机选取滤波后的子图像块Mh中的K个子图作为初始聚类的几何中心，其中K为聚类总数，大小为K = 70;
[0039](2c)计算滤波后的子图像块Mh中的每一个子图到所有几何中心的欧氏距离，并将该子图归入欧氏距离最小的聚类中；
[0040](2d)计算每一个聚类里面所有图像块的均值，得到新的几何中心；
[0041](2e)重复(2c)_(2d)，经过12次迭代，将滤波后的子图像块Mh聚为K类IC1, C2，…，Cq，…，CK}，q= 1，2，…，K;
[0042](2f)根据子图像块M滤波前后的一一对应关系，将相应的子图像块M直接划分为K类{Si，S2,…，S,，…，SK}，其中Sq表示第q个聚类，q = 1,2,…，K。 [0043]步骤3，用主成分分析的方法构造稀疏表示字典D。
[0044](3a)计算子图像块M的第q个聚类Sq的协方差矩阵' = SJ(Sq)U= 1，2，…，K，
[0045]其中T表示矩阵转置，K表示聚类总数；
[0046](3b)对协方差矩阵Ψ,进行奇异值分解，将其分解为正交矩阵Dq，特征值矩阵Λ q及正交矩阵\，即Ψ, = Dq Λ JqH，其中H表示共轭转置，取正交矩阵Dq为聚类Sq的子字典；
[0047](3c)将子字典Dq顺序拼接，得到稀疏表示字典D，q = 1，2，…，K。
[0048]步骤4，计算各子图Hii对应的无噪子图稀疏表示系数的估计值。
[0049](4a)先计算子图Hii与其周围大小为28X28的正方形邻域内所有子图的欧氏距
离，并保留欧氏距离最小的P个子图，P = 13,记为Ili = {ma,mi2,...,ZKijmiPj c M , Inij是
子图Hii的第j个相似子图，j = 1，2，".，Ρ，其中g是包含于符号，i = 1，2，...，L，L为子图像块总数；
[0050](4b)计算子图Hii的第j个相似子图Hiij在字典D下的稀疏表示系数，a = D1Hiij ；
[0051](4c)计算子图Hii对应的无噪子图稀疏表示系数的估计值5；,其中
W / I
__/ r*| t~~wf r\ l/{H)
Qij表示加权系数，% =(2 e零[Ls|(xf+,f)」/[[_[exp(2xJ+acp(2yi)jj , W是归一化因
P
子，，Xt为子图Hli的第t个像素值，yt为子图Hli的第j个相似子图Hlij的第t个像
J I
素值，T为子图的像素总数，h为衰减参数，h的大小设定为图像m噪声方差的12~15倍。
[0052]步骤5，构造目标函数。
[0053]利用对数变换后的图像m，稀疏表示字典D，并根据稀疏域噪声i =1，2，...，L的稀疏特性，构建如下目标函数:
[0054]4 = arS {||尺谓一Dai ||; + ||?i —戌 |} - 1...,L ,
[0055]其中Ri是抽取图像m中第i个图像块的矩阵，K i是稀疏约束项的正则参数，
I.112表示向量的2范数，11.111表示向量的I范数，眞表示降噪稀疏表示系数。
[0056]步骤6，求解目标函数，得到各子图像块HIi的最终降噪稀疏表示系数。[0057](6a)初始化目标函数的正则参数KiX),设总迭代次数为N = 15,迭代次数初始值为η = I ;
[0058](6b)用K均值聚类方法更新子图像块M的聚类，记更新后的聚类为
…Spj,其中Sf表示第n次迭代中的第q个聚类，q= 1，2，…，K，K表示
聚类总数；
[0059](6c)用主成分分析方法更新字典D，记第η次更新后的字典为D(n)；
[0060](6d)用稀疏表示系数估计方法更新第i个子图Hii对应的无噪子图稀疏表示系数
的估计值5；,记第η次更新后的估计值为i = 1，2，...，L，L表示子图像块的总数；
[0061](6e)用软阈值收缩算子计算并更新降噪稀疏表示系数离》得到第η次更新后的第i个降噪稀疏表示系数为Sf=SOfthw -却,)4 ,其中‘soft’表示表示经典的软阈值收缩
算子，#表示第i个子图Hii在字典D(n)下的稀疏表示系数，K i是目标函数稀
疏约束项的正则参数；
[0062](6f)用反稀疏变换更新降噪后的子图像块龟，记第η次更新后的第i个降噪子图像块为却?;η>,其中D(n)为第η次更新后的字典；
[0063](6g)用第η次更新后的第i个降噪子图像块<>替换子图像块M中第i个子图像块叫；
[0064](6h)判断是否满足总迭代次数N，若不满足，将迭代次数η加1，重复^b)-^f)继续迭代；若满足总迭代次数，则用最后一次迭代得到的降噪稀疏表示系数(if1作为最终的
降噪稀疏表示系数:為。
[0065]步骤7，获得抑斑后的SAR图像。
[0066](7a)对最终的降噪稀疏表示系数《进行反稀疏变换，得到降噪后的子图像块Ai = Dai，i = 1，2，...，L ;
[0067](7b)将降噪后的各个子图像块进行拼接，得到完整的图像，并对该图像进行反对数变换，得到去噪后的SAR图像。
[0068]本发明的仿真效果可以通过下述仿真实验加以说明
[0069]1.仿真条件
[0070]运行平台配置:
[0071]CPU:1nter (R)Core(TM)i565003.20GHz ；
[0072]内存:8GB(三星 DDR31600MHz)；
[0073]操作系统:Windows7旗舰版64位SPl ;
[0074]仿真软件: MATLABR (20 Ilb)。
[0075]2.仿真内容与结果
[0076]仿真1，用现有Frost抑斑方法对图2进行抑斑，并获取抑斑前后的比值图，结果如图3和图4。其中:
[0077]图3(a)为对图2(a)所示的梯田图像抑斑后的结果，
[0078]图3(b)为对图2(b)所示的港口图像抑斑后结果，
[0079]图4(a)为图2(a)与图3(a)各像素点的比值图，
[0080]图4(b)为图2(b)与图3(b)各像素点的比值图。
[0081]仿真2，用现有PPB抑斑方法对图2进行抑斑，并获取抑斑前后的比值图，结果如图5和图6。其中:
[0082]图5(a)为对图2(a)所示的梯田图像抑斑后的结果，
[0083]图5 (b)为对图2 (b)所示的港口图像抑斑后结果，
[0084]图6 (a)为图2 (a)与图5 (a)各像素点的比值图，
[0085]图6(b)为图2(b)与图5(b)各像素点的比值图。
[0086]仿真3，用本发明的抑斑方法对图2进行抑斑，并获取抑斑前后的比值图，结果如图7和图8。其中:
[0087]图7(a)为对图2(a)所示的梯田图像抑斑后的结果，
[0088]图7 (b)为对图2 (b)所示的港口图像抑斑后结果，
[0089]图8(a)为图2(a)与图7(a)各像素点的比值图，
[0090]图8(b)为图2(b)与图7(b)各像素点的比值图。
[0091]3.仿真结果分析
[0092](3.1)计算仿真1，仿真2和仿真3三种不同抑斑方法得到的各个抑斑后的图像的同质区域的等效视数，结果如表1所示。表中等效视数值越大说明抑斑能力越强。
[0093]表1三种不同方法得到抑斑后的图像的同质区域的等效视数
[0094]
【权利要求】
1.一种基于稀疏域噪声分布约束的SAR图像抑斑方法，包括如下步骤: (1)将原始SAR图像进行对数变换，得到变换后的图像m，并将m划分为大小相同的重叠的子图像块M = Im1, m2,...,Iiii,...，mL}，i = 1,2,..., L,其中L为划分的子图像块总数；(2)用K均值聚类方法将子图像块M聚为K类以，S2,…，S,，…，SK}，q= 1，2，…，K ; (3)用主成分分析的方法计算得到每一个聚类的子字典D,，并将子字典顺序拼接得到稀疏表不字典D ； (4)计算第i个含噪子图HIi的稀疏表示系数ai，Qi = Dy，并用稀疏域非局部均值方法计算Ini对应的无噪子图稀疏表示系数的估计值Sj，得到稀疏域噪声€.=ai ，i =12I.(5)利用变换后的图像m，并根据稀疏域噪声εi的稀疏特性，构建如下目标函数:
a, = argmm— Dai|丨;+ Kj |α;, i = 1,2,,..,1,, 其中Ri是抽取图像m中第i个图像块的矩阵，Ki是稀疏约束项的正则参数，Il.||2表示向量的2范数，Μ.M1 表示向量的I范数，必表示降噪稀疏表示系数； (6)求解目标函数，得到去噪后的SAR图像: (6a)初始化目标函数的正则参数K i > O, (6b)通过迭代收缩算法求解目标函数，得到每一个子图像块最终的降噪稀疏表示系数?..(6c)对最终的降噪稀疏表示系数忒进行反稀疏变换，得到降噪后的子图像块Mi = Dai ； (6d)拼接降噪后的子图像块得到完整的图像，并对该图像进行反对数变换，得到去噪后的SAR图像。
2.根据权利要求1所述的基于稀疏域噪声分布约束的SAR图像抑斑方法，其中第(2)步所述的用K均值聚类方法将子图像块M聚为K类，按如下步骤进行: (2a)对子图像块M进行高通滤波处理，记滤波后的子图像块为 = …,《f,…其中BIf是第i个子图Hii滤波后的子图像块，i = 1，2，...，L ; (2b)用K均值聚类的方法对滤波后的子图像块Mh进行聚类: (2bl)随机选取滤波后的子图像块Mh中的K个子图作为初始聚类的几何中心，其中K为聚类总数； (2b2)计算滤波后的子图像块Mh中的每一个子图到所有几何中心的欧氏距离，并将该子图归入欧氏距离最小的聚类中； (2b3)计算每一个聚类里面所有图像块的均值，得到新的几何中心； (2b4)重复(2b2)-(2b3)，经过12次迭代，将滤波后的子图像块Mh聚为K类IC1, C2，…，Cq，…，CK}，q= 1，2，…，K; (2c)根据子图像块M滤波前后的一一对应关系，将相应的子图像块M直接划分为K类(S1, S2,…，Sq,…，SK},其中Sq表不第q个聚类，q = I, 2,…，K。
3.根据权利要求1所述的基于稀疏域噪声分布约束的SAR图像抑斑方法，其中步骤(3)所述的用主成分分析的方法计算得到每一个聚类的子字典D,，按如下步骤进行: (3a)计算子图像块M的第q个聚类Sq的协方差矩阵Ψ, = Sq* (Sq)T, q = 1，2，…，K，其中T表示矩阵转置，K表示聚类总数；(3b)对协方差矩阵Ψ,进行奇异值分解，将其分解为正交矩阵Dq，特征值矩阵Aq及正交矩阵\，即
4.根据权利要求1所述的基于稀疏域噪声分布约束的SAR图像抑斑方法，其中步骤(4)所述的用稀疏域非局部均值的方法计算子图!1^对应的无噪子图稀疏表示系数的估计值g，按如下步骤进行:(4a)计算子图Hii与其周围大小为wXw的正方形邻域内所有子图的欧氏距离，并保留欧氏距离最小的P个子图，记为
5.根据权利要求1所述的基于稀疏域噪声分布约束的SAR图像抑斑方法，其中步骤(6)中的^b)所述的通过迭代收缩算法求解目标函数，得到每一个子图像块最终的降噪稀疏表示系数名,按如下步骤进行: (6bl)总迭代次数设为N，迭代次数初始值设为η = I ;(6b2)用K均值聚类方法更新子图像块M的聚类，记更新后的聚类为.…，Sf.SM ,其中表示第η次迭代中的第q个聚类
【文档编号】G06T5/00GK103971346SQ201410229349
【公开日】2014年8月6日 申请日期:2014年5月28日 优先权日:2014年5月28日
【发明者】王成龙, 赵光辉, 周丹凤, 石光明, 沈方芳 申请人:西安电子科技大学