CN102346908B - 基于稀疏表示的sar图像相干斑抑制方法 - Google Patents

Abstract

本发明公开了一种基于稀疏表示的SAR图像相干斑抑制的方法，主要解决现有技术对SAR图像去斑的细节信息不够完整，同质区域不够平滑的问题。其实现过程为：对原始SAR图像Y进行对数变换；将对数变换后的图像Y′分割为

的重叠分块；将图像块的自身信息作为控制因子对其进行稀疏表示；接着应用近似KSVD算法进行字典学习，得到自适应字典和更新后的稀疏表示系数；利用自适应字典和更新后的稀疏表示系数得到图像W，对W进行指数变换得到图像R；对由原始SAR图像Y和指数变换后的图像R得到的差值图像V进行非线性各向异性扩散，得到最终的去斑图像。本发明具有去斑结果同质区域光滑，点目标清晰可见以及边缘信息保持完整等优点，可用于SAR图像理解的预处理过程。

Description

基于稀疏表示的SAR图像相干斑抑制方法

技术领域

本发明属于图像处理技术领域，涉及一种SAR图像相干斑抑制的方法，可用于SAR图像的场景分析和图像理解的预处理。

背景技术

合成孔径雷达SAR由于具有全天候、全天时成像和穿透性强等诸多优点，在军用和民用领域有着广泛的应用。在SAR的成像过程中，由于电磁波的相干作用而产生的衰落效应的影响，使图像中产生随机分布的斑点噪声，即相干斑。相干斑的存在严重影响着对图像的理解与自动解译，降低了图像分割、目标识别提取的有效性。因此对SAR图像进行相干斑抑制是进行其他操作的必要前提。

SAR图像抑斑的一类主要方法是在空域进行处理的。空域滤波方法主要是基于SAR图像相干斑噪声的统计特性，如Lee滤波，Gamma MAP滤波等方法。空域的这些方法主要是通过在图像上取一个滑动窗口对图像进行滤波。其滤波能力与滑动窗的大小成正比，随着窗口的增大，其平滑效果变好，但是会导致图像边缘的模糊，使得细节信息丢失严重，在抑斑的同时很难较好的保留边缘和纹理细节。频域滤波主要是利用多分辨分析的方法，常见的有基于小波和众多的后小波变换分析的方法。

近年来，信号的稀疏表示在众多领域中取得了广泛的应用。设定一个过完备字典D∈R^n×L，其中字典D的每一列

称为一个字典原子，信号y∈Rⁿ可以表示为少数字典原子d_j的线性组合。信号的稀疏表示模型可以描述为：

s.t.Dα≈y，其中α是信号y的稀疏表示系数，||α||₀代表向量α的非零个数，并有||α||₀＜＜n。在该过程中，对于字典D的选择是一个关键的问题。对字典D的学习主要可分为两个方面：解析法和基于学习的方法。在解析法中，先设定数据遵从一定的数学模型，在该模型下推导解的形式来有效的表示信号。该方法因为没有明确的字典形式，而是通过一定的算法来描述其理念，所以称之为隐式字典法。这种方法包括小波分析及在小波基础上发展而来的众多多尺度分析，如Curvelets，Contourlets，Bandelets等。基于学习的方法则是根据机器学习的理论，从样本中训练一个字典矩阵。这种方法可以得到一个显式的字典模型，在该字典下，可以得到信号的表示系数。

前面所说的频域处理方法正是隐式字典法的应用，其在SAR图像的相干斑抑制方面已取得了广泛的应用，而基于学习训练的方法在SAR图像抑斑问题上的研究才刚刚起步。S.Foucher应用稀疏冗余表示对SAR相干斑抑制做了初步的尝试，他通过修正目标函数来重建对数变换后的图像，但是其处理后的图像细节丢失较为严重，边缘也出现了一定的模糊。

发明内容

本发明的目的在于克服上述已有技术的不足，提出了一种基于稀疏表示的SAR图像相干斑抑制方法，以有效地保持图像的边缘清晰度，提高SAR图像降斑的质量。

实现本发明目的的技术方案是：首先对原始SAR图像进行重叠分块，在对分割后的图像块的表示过程中，引入图像块的自身信息作为控制因子，针对不同的图像块内容产生稀疏度不同的表示。其次利用近似KSVD算法更新字典原子和稀疏表示系数，利用更新后的字典和稀疏表示系数得到图像块的稀疏表示。接下来对具有多个估计值的像素点进行平均，并对其做指数变换得到处理后的图像。在此基础上，利用非线性各向异性扩散方法平滑处理后的图像与原图像之间的差值图像得到点目标图像，把点目标图像加到处理后的图像上得到最终的去斑结果。具体步骤包括如下：

(1)对要进行相干斑抑制的SAR图像Y进行对数变换，得到对数变换后的图像Y′；

(2)将对数变换后的图像Y′分割为大小为

的重叠块；

(3)对分割后的每个图像块进行稀疏编码，在编码过程中引入该图像块的自身信息作为控制因子控制编码：

(3a)将待编码的每个图像块拉成列向量y，所有的列向量构成一个矩阵X；

(3b)选择初始字典D为离散余弦变换DCT字典；

(3c)根据下列公式应用正交匹配追踪算法OMP，实现对y的稀疏编码：

$\widehat{\mathrm{\α}}=\arg \underset{\mathrm{\α}}{\min }{\left|\right|\mathrm{D\α}-y\left|\right|}_2^2+f\left(y\right){\left|\right|\mathrm{\α}\left|\right|}_0,$

其中α是列向量y在字典D下的稀疏表示向量，

为迭代后的最终表示向量，||α||₀代表向量α的非零个数，

为设定的惩罚项，用来控制编码过程，

是列向量y的归一化L₂范数，σ为噪声标准差，其计算公式为：

σ＝median(|W_Y|-median(W_Y))/0.6745，

其中W_Y表示原始SAR图像Y经过小波分解之后的高频子带系数；

(3d)对矩阵X中所有的列向量进行步骤(3c)的编码操作，得到稀疏表示系数矩阵Λ；

(4)应用近似KSVD算法进行字典更新，得到更新后的字典D′和稀疏表示系数矩阵Λ′；

(5)利用更新后的字典D′和稀疏表示系数矩阵Λ′，得到处理后的图像W；

(6)对处理后的图像W进行指数变换，得到指数变换后的图像R；

(7)将原始SAR图像Y与指数变换后的图像R相减得到差值图像：V＝Y-R，对差值图像V进行非线性各向异性扩散，进行k次迭代得到点目标图像v^k；

(8)将点目标图像v^k加到指数变换后的图像R上，得到相干斑抑制后的图像I＝R+v^k。

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

(1)本发明由于在稀疏表示过程中引入了图像块的自身信息作为控制因子，这样对不同的图像块内容产生自适应的控制因子，在控制因子的作用下，当表示边缘细节信息时，控制因子较小，能使细节信息得到较好的保留；而在同质区域，较大的控制因子能使噪声得到较好的抑制。

(2)本发明由于对图像进行了重叠分块，利用重叠块的冗余信息，能够对SAR图像纹理区域的方向性细节信息进行提取，从而保证了图像纹理信息的完整性。

(3)本发明由于采用了对差值图像进行非线性各向异性扩散的方法，能有效地提取差值图像中的点目标，将提取的点目标图像加到处理后的图像上，有效地保留了图像中的点目标。

仿真结果表明，本发明对SAR图像的相干斑抑制能得到比现有技术更好的效果。

附图说明

图1是本发明的流程图；

图2是用本发明与已有方法对一幅2视幅度SAR图像Field的相干斑抑制效果对比图；

图3是用本发明与已有方法对一幅4视强度SAR图像Horsetrack的相干斑抑制效果对比图；

图4是用本发明与已有方法对一幅6视强度SAR图像Bedford的相干斑抑制效果对比图。

具体实施方式

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

步骤1.输入要进行相干斑抑制的SAR图像Y，并对其作对数运算，得到对数变换后的图像Y′。

步骤2.将对数变换后的图像Y′分割为大小为

的重叠块。

选取n的大小为64，即把图像分割为大小为8×8的重叠块，以图像中的每个像素点为起始点分别向右和向下提取大小为8×8的图像块，并以起始像素点的坐标来标记该块在原始图像中的位置。实验中采用的SAR图像大小为256×256，则分割后的图像块数目为(256-7)×(256-7)＝249²＝62001个。

步骤3.对分割后的每个重叠块进行稀疏编码，在编码过程中引入该图像块的自身像素信息作为控制因子控制编码。

(3a)将待编码的每个图像块拉成列向量y，所有的列向量构成一个矩阵X；

(3b)选择初始字典D为离散余弦变换DCT字典；

(3c)根据下列公式应用正交匹配追踪算法OMP，实现对y的稀疏编码：

$\widehat{\mathrm{\α}}=\arg \min {\left|\right|\mathrm{D\α}-y\left|\right|}_2^2+f\left(y\right){\left|\right|\mathrm{\α}\left|\right|}_0,$

其中α是列向量y在字典D下的稀疏表示向量，

为迭代后的稀疏表示向量，||α||₀代表向量α的非零个数。

为设定的惩罚项，用来控制编码过程。

是列向量y的归一化L₂范数，σ为噪声标准差，其计算公式为：

σ＝median(|W_Y|-median(W_Y))/0.6745，

其中W_Y表示原始SAR图像经过小波分解之后的高频子带系数；

(3d)对所有的列向量进行(3c)操作，得到稀疏表示的系数矩阵Λ，并有如下关系式成立：X≈DΛ。

步骤4.应用近似KSVD算法进行字典更新，得到更新后的字典D′和稀疏表示系数矩阵Λ′。

(4a)将字典D的每一列称为一个字典原子d_i，i＝1，2，…，L，L为字典原子数目，更新第i个字典原子d_i时，从矩阵X中抽取出应用字典原子d_i进行稀疏表示的信号子矩阵X_I及其对应的稀疏表示系数子矩阵Λ_I，此时信号子矩阵X_I与字典D和稀疏表示系数子矩阵Λ_I之间的关系为：X_I≈DΛ_I；

(4b)将稀疏表示系数子矩阵Λ_I中的第i行元素置为零得到

计算此时信号子矩阵X_I与

之间的残差E_I：

$E_I=X_I-D{\stackrel{\‾}{\mathrm{\Λ}}}_I;$

(4c)通过将残差E_I最小化，即得到更新后的字典原子d：

d＝E_Iα/||E_Iα||₂，

其中，α^T是α的转置向量，得到更新后的字典原子d之后，通过下式得到对应的更新后的稀疏表示系数：

${\mathrm{\α}}^{\′}=E_I^{T}d,$

其中

是E_I的转置矩阵；

(4d)对所有字典原子重复进行步骤(4a)-(4c)，得到更新后的字典D′和稀疏表示系数矩阵Λ′。

步骤5.利用更新后的字典D′和稀疏表示系数矩阵Λ′，得到处理后的图像W：

(5a)应用更新后的字典D′和稀疏表示系数矩阵Λ′，得到估计的矩阵X′，即X′＝D′Λ′；

(5b)将所述X′中的每一列x反变换为大小为

的图像块x′，放回到原图像初始分割的位置处，即x′＝R_ijx，其中R_ij代表抽取矩阵，将抽取的图像块再放置到原图像中初始分割时的位置处；

(5c)对每个像素点处的多个估计值进行平均，得到处理后的图像W。

步骤6.对处理后的图像W进行指数变换，即对图像W中的每个像素值进行指数运算，得到指数变换后的图像R。

步骤7.将原始SAR图像Y与指数变换后的图像R相减得到差值图像：V＝Y-R，对差值图像V进行非线性各向异性扩散，进行k次迭代得到点目标图像v^k。

首先令点目标图像为差值图像，即：v⁰＝V，然后按照如下的公式对点目标图像进行迭代：

$v_{i,j}^k=v_{i,j}^{k-1}+t\·\underset{\stackrel{r,s=-1}{(r,s)\≠0}}{\overset{1}{\mathrm{\Σ}}}\frac{g\left({\sqrt{2}}^{1-\left|r\right|-\left|s\right|}\right|v_{i+r,j+s}^{k-1}-v_{i,j}^{k-1}\left|\right)(v_{i+r,j+s}^{k-1}-v_{i,j}^{k-1})}{r^2+s^2},$

其中，

是点目标图像v^k中的任一像素值，t为平滑尺度参数，k是迭代次数，

为扩散方程，ρ是扩散方程中的一个常数，

代表第k次迭代时点目标图像v^k中位于第i行，第j列的像素值，(r，s)代表以像素点v(i，j)为中心选取的一个方形窗的坐标，r＝-1，…，1，s＝-1，…，1且r和s不能同时为零，在上式中，选择t＝1/8，ρ＝5。

步骤8.将点目标图像v^k加到指数变换后的图像R上，得到相干斑抑制后的图像I＝R+v^k。

本发明的效果可以通过如下仿真实验进一步说明。

1.仿真条件

硬件平台为：Intel Core2 Duo CPU E65502.33GHZ、2GB RAM

软件平台为：Matlab 7.1

2.仿真内容

应用本发明和现有的增强Lee滤波，Gamma MAP滤波方法分别对三幅SAR图像进行去斑实验，从同质区域的内部一致性，边缘信息和纹理信息的完整性以及点目标的清晰性来评价这些方法的性能。

3.仿真实验结果

A)2视幅度SAR图像Field的实验结果

应用本发明以及增强Lee滤波和Gamma MAP滤波方法对图2(a)进行去斑仿真，其效果比较图如图2(b)-(d)所示，其中：图2(b)为增强Lee滤波方法对图2(a)进行去斑的结果；图2(c)为Gamma MAP滤波方法对图2(a)进行去斑的结果；图2(d)为本发明对图2(a)进行去斑的结果。由图2(d)可见，本发明中，斑点噪声得到了有效的抑制，边缘细节信息得到了较好的保留，抑斑后的图像整体效果清晰。相比而言，增强Lee滤波方法处理后的图像斑点噪声抑制的不够彻底，Gamma MAP滤波方法在一定程度上抑制了斑点噪声，但边缘比较模糊，点目标也出现了严重的模糊。

B)4视强度图像Horsetrack的实验结果

应用本发明以及增强Lee滤波和Gamma MAP滤波方法对图3(a)进行去斑仿真，其效果比较图如图3(b)-(d)所示，其中：图3(b)为增强Lee滤波方法对图3(a)进行去斑的结果；图3(c)为Gamma MAP滤波方法对图3(a)进行去斑的结果；图3(d)为本发明对图3(a)进行去斑的结果。由图3(d)可见，本发明处理后的图像边缘得到了很好的保留，同质区域的斑点噪声得到了理想的抑制。相比而言，增强Lee滤波方法处理后的图像斑点噪声依然较为明显，Gamma MAP滤波方法在一定程度上抑制了斑点噪声，但图像整体效果较为模糊。

C)6视强度图像Bedford的实验结果

应用本发明以及增强Lee滤波和Gamma MAP滤波方法对图4(a)进行去斑仿真，其效果比较图如图4(b)-(d)所示，其中：图4(b)为增强Lee滤波方法对图4(a)进行去斑的结果；图4(c)为Gamma MAP滤波方法对图4(a)进行去斑的结果；图4(d)为本发明对图4(a)进行去斑的结果。由图4(d)可见，本发明处理后的图像整体效果清晰，同质区域的斑点噪声得到了很好的抑制。相比而言，增强Lee滤波方法并没有较好的滤除图像的斑点噪声，GammaMAP滤波方法处理后的图像整体效果模糊，细节信息丢失较为严重。

Claims (3)

1.一种基于稀疏表示的SAR图像相干斑抑制方法，包括以下步骤：

（1）对要进行相干斑抑制的SAR图像Y进行对数变换，得到对数变换后的图像Y′；

（2）将对数变换后的图像Y′分割为大小为的重叠图像块,n的大小为64；

（3）对分割后的每个图像块进行稀疏编码：

（3a）将待编码的每个图像块拉成列向量y，所有的列向量构成一个矩阵X；

（3b）选择初始字典D为离散余弦变换DCT字典；

（3c）根据下列公式应用正交匹配追踪算法OMP，实现对y的稀疏编码：

$\widehat{\mathrm{\α}}=\arg \underset{\mathrm{\α}}{\min }{\left|\right|\mathrm{D\α}-y\left|\right|}_2^2+f\left(y\right){\left|\right|\mathrm{\α}\left|\right|}_0,$

其中α是列向量y在字典D下的稀疏表示向量，

为迭代后的稀疏表示向量，||α||₀代表向量α的非零个数，

为设定的惩罚项，用来控制编码过程，

是列向量y的归一化L₂范数，σ为噪声标准差，其计算公式为：

σ＝median(|W_Y|-median(W_Y))/0.6745，

其中W_Y表示原始SAR图像Y经过小波分解之后的高频子带系数；

（3d）对矩阵X中所有的列向量进行步骤（3c）的编码操作，得到稀疏表示系数矩阵Λ；

（4）应用近似KSVD算法进行字典更新，得到更新后的字典D′和稀疏表示系数矩阵Λ′；

（5）利用更新后的字典D′和稀疏表示系数矩阵Λ′，得到处理后的图像W：

（5a）应用更新后的字典D′和稀疏表示系数矩阵Λ′，得到估计的矩阵X′，即X′＝D′Λ′；

（5b）将上述X′中的每一列x反变换为大小为

的图像块x′，放回到原图像初始分割的位置处，即x′＝R_ijx，其中R_ij代表抽取矩阵，将抽取的图像块放置到原图像中初始分割时的位置处；

（5c）对每个像素点处的多个估计值进行平均，得到处理后的图像W；

（6）对处理后的图像W进行指数变换，得到指数变换后的图像R；

（7）将原始SAR图像Y与指数变换后的图像R相减得到差值图像：V＝Y-R，对差值图像V进行非线性各向异性扩散，进行k次迭代得到点目标图像v^k；

（8）将点目标图像v^k加到指数变换后的图像R上，得到相干斑抑制后的图像I＝R+v^k。

2.根据权利要求1所述的基于稀疏表示的SAR图像相干斑抑制方法，其中步骤(4)所述的应用近似KSVD进行字典更新，得到更新后的字典D′和稀疏表示系数矩阵Λ′，按如下步骤进行：

（4a）将字典D的每一列称为一个字典原子d_i，i＝1,2,…,L，L为字典的原子数目，更新第i个字典原子d_i时，从矩阵X中抽取出应用字典原子d_i进行稀疏表示的信号子矩阵X_I及其对应的稀疏表示系数子矩阵Λ_I，此时信号子矩阵X_I与字典D和稀疏表示系数子矩阵Λ_I之间的关系为：X_I≈DΛ_I；

（4b）将稀疏表示系数子矩阵Λ^I中的第i行元素置为零得到

计算此时信号子矩阵X_I与

之间的残差E_I：

$E_I=X_I-D{\stackrel{\‾}{\mathrm{\Λ}}}_I;$

（4c）通过将残差E_I最小化，即得到更新后的字典原子d：

d＝E_Iα/||E_Iα||₂，

其中，α^T是α的转置向量，得到更新后的字典原子d之后，通过下式得到对应的更新后的稀疏表示系数：

${\mathrm{\α}}^{\′}=E_I^{T}d,$

其中是E_I的转置矩阵；

（4d）对所有字典原子重复进行步骤（4a）-（4c），得到更新后的字典D′和稀疏表示系数矩阵Λ′。

3.根据权利要求1所述的基于稀疏表示的SAR图像相干斑抑制方法，其中步骤(7)所述的对差值图像V进行非线性各向异性扩散，进行k次迭代得到点目标图像v^k，是按照如下公式进行迭代：

$v_{i,j}^k=v_{i,j}^{k-1}+t\·\underset{(r,s)\≠0}{\overset{1}{\underset{r,s=-1}{\mathrm{\Σ}}}}\frac{g\left({\sqrt{2}}^{1-\left|r\right|-\left|s\right|}\right|v_{i+r,j+s}^{k-1}-v_{i,j}^{k-1}\left|\right)(v_{i+r,j+s}^{k-1}-v_{i,j}^{k-1})}{r^2+s^2},$

其中t为平滑尺度参数，k是迭代次数，

为扩散方程，ρ是扩散方程中的一个常数，v⁰＝V，

代表第k次迭代过程中位于第i行，第j列的像素值，(r,s)代表以像素点v(i,j)为中心选取的一个方形窗的坐标，r＝-1,…,1，s＝-1,…,1且r和s不能同时为零，在上式中，选择t＝1/8，ρ＝5。

Images