CN107392861B - 一种基于高斯比例混合模型的稀疏表示sar图像降斑方法 - Google Patents

Abstract

本发明涉及一种基于高斯比例混合模型的稀疏表示SAR图像降斑方法，属于图像滤波研究领域。本方法首先，建立单个图像块的稀疏表示模型；然后根据相干斑的统计特性与贝叶斯估计原理，将稀疏系数

Description

一种基于高斯比例混合模型的稀疏表示SAR图像降斑方法

技术领域

本发明涉及一种基于高斯比例混合模型的稀疏表示SAR图像降斑方法，属于图像滤波研究领域。

背景技术

合成孔径雷达(Synthetic Aperture Radar，SAR)因其在方位向及距离向的高分辨率成像能力成为微波遥感领域具有代表性的雷达系统，在军用及民用方面都受到广泛关注。雷达在发射电磁信号照射目标的过程中，目标的随机散射信号与发射信号之间的干涉会产生相干斑噪声，严重影响了图像质量，因此对于相干斑噪声的抑制尤为重要。消除斑点噪声最直接的方法是采用多视处理，但该方法会降低方位向分辨率，更为合理的方式是使用滤波技术在有效滤除噪声的同时尽可能多地保留图像的辐射特性及纹理特征。

为了有效降低相干斑对SAR图像质量的影响，国内外学者已经提出了各种SAR图像滤波方法，主要可分为基于统计特性的滤波方法、偏微分扩散类算法、基于结构相似度的非局部均值滤波、基于稀疏度的变换域滤波等。早期较经典的Lee滤波、Kuan滤波、Frost滤波、MAP滤波等均是基于统计特性的滤波方法；偏微分扩散类算法利用偏微分方程各向异性的特点将相干斑抑制问题转化为泛函极值问题，通过变分法和数值计算得到噪声抑制后的SAR图像，目前提出的有自蛇扩散、P-M扩散、SRAD算法、DPAD算法等。

近几年随着非局部均值滤波思想及稀疏表示的普及，与之相关的针对SAR图像滤波的算法也相继出现，其中最具有代表性且滤波性能较优的算法主要是PPB算法、SAR-BM3D算法及FANS算法。然而PPB算法虽然能较好地抑制相干斑，但是会引入大量人工伪影，导致图像细节信息大量丢失，SAR-BM3D及FANS算法在细节信息的保持能力上相对较优，但是降斑效果不如PPB，无法尽可能多的抑制相干斑噪声。因此，如何在尽可能多地保留图像的局部细节信息及纹理特征的同时实现对相干斑的良好抑制成为当前亟待解决的问题。

发明内容

为了在尽可能多地保留原SAR图像细节信息及纹理特征的情况下，滤除图像中的相干斑噪声，提高图像局部平滑区域的等效视数，本发明提出了一种基于高斯比例混合模型的稀疏表示SAR图像降斑方法，该方法能够实现在有效抑制相干斑的同时较好地保留SAR图像的细节信息，达到良好的图像重建质量。

本发明为解决其技术问题采用如下技术方案：

一种基于高斯比例混合模型的稀疏表示SAR图像降斑方法，首先，建立单个图像块的稀疏表示模型；然后根据相干斑的统计特性与贝叶斯估计原理，将稀疏系数α用GSM模型进行表示，得到优化模型；同时对SAR图像进行分类，根据分类结果建立稀疏模型；最后，利用凸优化方法对上述模型求解，得到最优的稀疏表示，进而得到降噪图像。

一种基于高斯比例混合模型的稀疏表示SAR图像降斑方法，具体步骤如下：

步骤1，对单一图像块进行稀疏表示建模分析，给出最基本的凸优化数学模型；

步骤2，对步骤1中提出的数学模型进行分析，对稀疏系数α进行GSM建模，结合相干斑的统计特性与贝叶斯估计将该GSM模型代入到原有的凸优化模型中，由此转换求解域，得到更易于求解的凸优化方程；

步骤3，将SAR图像进行分类，得到相似图像块的集合，将步骤2中对于单一图像块建立的数学模型推广到图像块的集合中，得到图像块集合的稀疏表示模型；

步骤4，对步骤3中提出的稀疏表示模型，使用凸优化方法进行求解，得到最优解

分别对图像块的集合进行同样的求解过程，根据所得到的最优解对图像进行重建，实现噪声的滤除。

所述步骤1中对单一图像块进行稀疏表示建模分析，给出最基本的凸优化数学模型的具体方法为：

步骤11，给出稀疏表示的基本数学模型：

假设图像块的大小为

像素，按照字典序排列形成列向量x∈Rⁿ，针对稀疏域建模，定义一个字典矩阵D∈R^n×k，其中k≥n，表明字典是冗余的，由此根据所提出的模型使用该字典稀疏表示每个图像块x，表示形式如下：

上式中α∈R^K为稀疏系数，||α||₀为系数α的0范数，

为利用凸优化工具得到的系数α的最优解；

步骤12，对步骤11中的基本数学模型进行等价变换：

使用l₁范数来替代原始的非凸优化问题，上述模型等价为：

其中λ为正则化参数，||α||₁为系数α的1范数，||x-Dα||₂为原图像与重建后图像的误差的2范数。

所述步骤2的具体过程如下：

步骤21，实现稀疏系数的GSM建模：

使用GSM模型对稀疏系数α进行建模，将向量α分解为一个高斯向量β和一个标量乘性因子θ的逐点乘积，即α_i＝θ_iβ_i，其中θ_i表示概率为P(θ_i)的正标量，α_i为稀疏系数的一个元素，β_i为高斯向量的一个元素，假设θ_i是独立同分布的且与β_i无关，则α的GSM先验信息表示为如下：

其中：P(θ_i)为θ_i的概率，P(α_i|θ_i)为θ_i条件下α_i的概率；

步骤22，联合贝叶斯估计与相干斑的统计特性推导新的稀疏表示数学模型：

对于每一个图像块x∈Rⁿ，其稀疏表示数学模型写为如下形式：

x＝y·u＝y+y·(u-1)

＝y+v＝Dα+v

其中x表示观测图像块，y表示不含噪图像块，u表示相干斑，v＝y·(u-1)表示等效加性噪声，α为稀疏系数,D为字典；在振幅形式上，相干斑的统计特性服从Nakagami分布，概率密度函数表示为如下：

其中L表示等效视数，Γ(·)表示Gamma函数；由于y与u相互独立，且u-1的均值为0，因此v的均值为0；

根据贝叶斯准则，对于一个已知观测信号x＝Dα+v，表示为如下MAP估计：

(α,θ)＝argmaxlogP(x|α,θ)P(α,θ)

＝argmaxlogP(x|α)+logP(α|θ)+logP(θ)

其中：P(x|α,θ)表示在α和θ的条件下x的概率，P(α,θ)表示在θ的条件下α的概率，P(x|α)为α条件下x的概率，P(α|θ)为θ条件下α的概率，P(θ)为θ的概率；上式中，先验项P(α|θ)表示为如下形式：

结合相干斑的概率密度函数得到：

P(x|α)＝P(v)＝P(y·(u-1))＝P(y)·P(u-1)

＝C·(u-1)^2L-1exp(-L(u-1)²)

其中，

Γ(L)L的Gamma函数，P(v)为v的概率，P(y)为y的概率；记

选取

结合上述各项式子，推导出如下稀疏表示模型：

其中：||x||₂为x的2范数，使用log(θ_i+ε)代替logθ_i，ε为一个很小的正常数，记

为log(θ+ε)，上述公式简化为：

注意到原始的GSM模型的矩阵形式为α＝Λβ，其中Λ＝diag(θ_i)∈R^K×K为一个对角矩阵，表示被选图像块的方差域，R^K×K为K×K阶的实数域；因此，上述稀疏编码问题从α域转换到β域，形式如下：

上述稀疏表示模型即为单一图像块的GSM稀疏编码模型，其中||x-DΛβ||₂为误差的2范数，||β||₂为β的2范数。

所述步骤3中将SAR图像进行分类，得到相似图像块的集合，并得到图像块集合的数学模型的具体方法为：

步骤31：将SAR图像进行相似图形块的分类：

使用概率估计来代替欧氏距离作为相似度的评价指标：

上式中，ω(y_i,y_j)是在平均过程中作为权重的相似度，表示了在这两个图像块中隐藏的降噪信号的概率是相同的；y_i和y_j分别表示第i个和第j个图像块的观测值；x_i和x_j分别表示第i个和第j个图像块的降噪值，Ω表示在图像块中的所有像素的集合；上述公式是一种最大似然估计且假设x_i(k)与x_j(k)条件独立，最终简化为：

将SAR图像相干斑的概率分布代入简化后的公式，得到最终的权重表达式：

其中，

z_i(k)为第i个图像块的观测值求根号，z_j(k)为第j个图像块的观测值求根号；选取S个权重值大于其他的图像块组成一个相似块的集合Φ_i，由此将一幅SAR图像分成相似图像块的集合；

步骤32，推导一个图像块集合的稀疏表示模型：

对于包含m个相似块的集合，如果考虑GSM模型的同步稀疏编码，根据步骤22中得到的单个图像块的稀疏表示模型得到群稀疏表达式：

其中，X＝[x₁,x₂…,x_m]表示m个相似块的集合，x₁为第一个图像块，x₂为第二个图像块，x_m为第m个图像块，||x_i||₂为x_i的二范数，A＝ΛB表示GSM模型的稀疏系数的群表示，其中B＝[β₁,β₂…,β_m]∈R^K×m表示稀疏系数对应的高斯向量的集合，β₁为第一个高斯向量，β₂为第二个高斯向量，β_m为第m个高斯向量，R^K×m为K×m阶的实数矩阵,||B||_F表示B的F范数，||X-DΛB||_F表示重建误差的F范数。

所述步骤4中使用凸优化方法求解群稀疏编码模型的具体方法为：

步骤41，保证θ的值不变，求B：

当θ的值固定时，群稀疏表达式简化为如下形式：

由X＝DA，x_i＝Dα_i，D为正交矩阵，上式简化为如下：

其中：||A-ΛB||_F为误差的2范数，||α_i||₂为α_i的2范数；

令

则f(B)写为：

f(B)＝B^T(Λ^TΛ+dI)B-2B^TΛ^TA+A^TA

当函数f(B)为最小值时相应的B通过f(B)的一阶导数为0时来求解，即：

根据上式求得B的值为：

B＝(Λ^TΛ+dI)^-1Λ^TA

步骤42，保证B的值不变，求θ

当B的值固定时，群稀疏表达式简化为如下形式：

同步骤41上式简化为：

令

则上式写为：

式中α^j和β^j分别表示A和B的第j行，θ_j表示与α^j对应的θ值，

为第j行α的2范数，||β^j||₂为第j行β的2范数，

令

b_j＝-2α^j(β^j)^T，则

上式分解为一系列标量最小化问题来进行求解：

令g(θ_j)＝a_jθ_j ²+b_jθ_j+clog(θ_j+ε)，最小化问题同样用一阶导数为0的方法进行求解，即

根据步骤41和步骤42所求得的B和θ得到X的估计值：

其中

和

分别为B和Λ的估计值。

步骤43，根据步骤41和42的方法求解一幅SAR图像的数学模型：

对于受相干斑影响的SAR图像，将其划分为N个不同的图像块集合，每个集合里面包含m个相似图像块，则图像的稀疏编码优化问题表示为如下：

上式中X_j表示第j个图像块集合，

表示第j个图像块集合里的第i个图像块，||X_j-DΛ_jB_j||_F为误差的F范数，||x_i ^j||₂为

的2范数，||B_j||_F为B_j的F范数，η为正则化参数，根据经验进行设置。

本发明的有益效果如下：

本发明基于GSM(高斯比例混合)模型和稀疏表示原理的数学模型，在模型建立以及求解问题上都易于理解，对于凸优化问题采用迭代正则化的方式，能够很好地抑制相干斑，有效提高图像同质区的等效视数，且能很好地保留异质区的细节信息及纹理特征。产生该优点的原因是在稀疏编码的过程，联合相干斑的统计特性及贝叶斯估计可以很好地将噪声滤除，GSM先验信息的引入考虑到了稀疏系数的全局与局部相关性，使得细节信息得到保留。由于在编码过程中对图像块进行了分类，简化了数学模型。

具体实施方式

下面对本发明创造做进一步详细说明。

实施例

本实施例的一种基于GSM模型的稀疏表示SAR图像降斑方法，首先使用概率估计的方法将SAR图像分成若干个服从相同概率统计特性的相似块的集合，根据稀疏表示原理对某个集合里的图像块进行建模，得到凸优化数学模型。然后联合贝叶斯估计与相干斑的概率密度函数，对该模型里的稀疏系数进行GSM建模，得到基于GSM的稀疏表示模型。由于每个集合里的相似图像块服从相同的数学统计，因此满足同一稀疏表示。用以上建模方法对若干个图像块的集合分别建立数学模型，模型的求解使用迭代正则化的方式，得到该凸优化模型的最优解之后，即可实现SAR图像的重建。

具体实现步骤如下：

1、对单一图像块进行稀疏表示建模分析，给出最基本的凸优化数学模型

步骤11，给出稀疏表示的基本数学模型：

假设图像块的大小为

像素，按照字典序排列形成列向量x∈Rⁿ，针对稀疏域建模，定义一个字典矩阵D∈R^n×k(k≥n，表明字典是冗余的)，n为图像块像素个数，由此根据所提出的模型可以使用该字典稀疏表示每个图像块x，表示形式如下：

满足Dα≈x

上式中α∈R^K为稀疏系数，||α||₀为系数α的0范数，

为利用凸优化工具得到的系数α的最优解。

步骤12，对步骤11中的基本数学模型进行等价变换：

由于步骤11中的l₀范数求解问题难以实现，因此使用l₁范数来替代原始的非凸优化问题，则上述模型等价为：

其中λ为正则化参数，||α||₁为系数α的1范数，||x-Dα||₂为原图像与重建后图像的误差的2范数。

2、对稀疏系数α进行GSM建模，得到更易于求解的凸优化方程

步骤21，实现稀疏系数的GSM建模：

求解步骤12中的l₁范数最小化问题相当于推导服从独立同分布的拉普拉斯先验信息的α的最大后验估计(MAP)。使用GSM模型对稀疏系数α进行建模，将向量α分解为一个高斯向量β和一个标量乘性因子θ的逐点乘积，即α_i＝θ_iβ_i，其中θ_i表示概率为P(θ_i)的正标量，α_i为稀疏系数的一个元素，β_i为高斯向量的一个元素。假设θ_i是独立同分布的且与β_i无关，则α的GSM先验信息可以表示为如下：

其中α_i表示单个像素的稀疏系数，P(θ_i)为θ_i的概率，P(α_i|θ_i)为θ_i条件下α_i的概率。

步骤22，联合贝叶斯估计与相干斑的统计特性推导新的稀疏表示数学模型：

一般而言，假设一幅SAR图像受到相干斑影响，则反向散射信号会被乘性噪声污染，对于每一个图像块x∈Rⁿ，其稀疏表示数学模型可以写为如下形式：

x＝y·u＝y+y·(u-1)

＝y+v＝Dα+v

其中x表示观测图像块，y表示不含噪图像块，u表示相干斑，v＝y·(u-1)表示等效加性噪声，α为稀疏系数,D为字典。在振幅形式上，相干斑的统计特性服从Nakagami分布，概率密度函数表示为如下：

其中L表示等效视数，Γ(·)表示Gamma函数，u为相干斑。由于y与u相互独立，且u-1的均值为0，因此v的均值为0。

根据贝叶斯准则，对于一个已知观测信号x＝Dα+v，可以表示为如下MAP估计：

(α,θ)＝argmaxlogP(x|α,θ)P(α,θ)

＝argmaxlogP(x|α)+logP(α|θ)+logP(θ)

其中：P(x|α,θ)表示在α和θ的条件下x的概率，P(α,θ)表示在θ的条件下α的概率，P(x|α)为α条件下x的概率，P(α|θ)为θ条件下α的概率，P(θ)为θ的概率。上式中，先验项P(α|θ)可以表示为如下形式：

结合相干斑的概率密度函数可以得到：

P(x|α)＝P(v)＝P(y·(u-1))＝P(y)·P(u-1)

＝C·(u-1)^2L-1exp(-L(u-1)²)

其中，

Γ(L)L的Gamma函数，P(v)为v的概率，P(y)为y的概率。记

选取

结合上述各项式子，可以推导出如下稀疏表示模型：

其中：||x||₂为x的2范数，使用log(θ_i+ε)代替logθ_i，ε为一个很小的正常数，记

为log(θ+ε)，上述公式可以简化为：

注意到原始的GSM模型的矩阵形式为α＝Λβ其中Λ＝diag(θ_i)∈R^K×K为一个对角矩阵，表示一个被选图像块的方差域，R^K×K为K×K阶的实数域。因此，上述稀疏编码问题可以从α域转换到β域，形式如下：

上述稀疏表示模型即为需要求解的单一图像块的基于GSM的稀疏编码模型，其中||x-DΛβ||₂为误差的2范数，||β||₂为β的2范数。

3、将SAR图像进行分类，并推导图像块集合的数学模型

步骤31：将SAR图像进行相似图形块的分类：

通常，对于加性噪声，块之间的相似度通过欧氏距离的均值进行估计，相似度高的块才能被选择进行加权平均。但是对于SAR图像块相似度的衡量并不能使用欧氏距离，在此使用概率估计来代替欧氏距离作为相似度的评价指标。

上式中，ω(y_i,y_j)是在平均过程中作为权重的相似度，表示了在这两个图像块中隐藏的降噪信号的概率是相同的。y_i(k)和y_j(k)分别表示第i个和第j个图像块的观测值。x_i(k)和x_j(k)分别表示第i个和第j个图像块的降噪值。Ω表示在图像块中的所有像素的集合。上述公式是一种最大似然估计且假设x_i(k)与x_j(k)条件独立，最终可以简化为：

将SAR图像相干斑的概率分布代入简化后的公式，可以得到最终的权重表达式：

其中，

z_i(k)为第i个图像块的观测值求根号，z_j(k)为第j个图像块的观测值求根号。选取S个权重值大于其他的图像块组成相似块的集合Φ_i，由此可将一幅SAR图像分成若干个不同类型的相似图像块的集合。

步骤32，推导图像块集合的稀疏表示模型：

对于相似块的集合，里面的每个块所对应的稀疏系数α应该服从相同的先验信息，即它的概率密度函数包含相同的θ。因此，对于包含m个相似块的集合，如果考虑GSM模型的同步稀疏编码，根据步骤22中得到的单个图像块的稀疏表示模型可以得到群稀疏表达式：

其中，X＝[x₁,x₂…,x_m]表示m个相似块的集合，x₁为第一个图像块，x₂为第二个图像块，x_m为第m个图像块，||x_i||₂为x_i的2范数，A＝ΛB表示GSM模型的稀疏系数的群表示，其中B＝[β₁,β₂…,β_m]∈R^K×m表示稀疏系数对应的高斯向量的集合，β₁为第一个高斯向量，β₂为第二个高斯向量，β_m为第m个高斯向量，R^K×m为K×m阶的实数矩阵,||B||_F表示B的F范数，||X-DΛB||_F表示重建误差的F范数。

4、求解SAR图像的稀疏表示优化模型

步骤41，保证θ的值不变，求B：

当θ的值固定时，群稀疏表达式可以简化为如下形式：

由X＝DA，x_i＝Dα_i，D为正交矩阵，上式可以简化为如下：

其中：||A-ΛB||_F为误差的2范数，||α_i||₂为α_i的2范数。

令

则f(B)可以写为：

f(B)＝B^T(Λ^TΛ+dI)B-2B^TΛ^TA+A^TA

其中：Λ为θ的集合，Λ^T为Λ的转置，A^T为A的转置，B^T为B的转置。

当函数f(B)为最小值时相应的B可以通过f(B)的一阶导数为0时来求解，即：

根据上式求得B的值为：

B＝(Λ^TΛ+dI)^-1Λ^TA

步骤42，保证B的值不变，求θ

当B的值固定时，群稀疏表达式可以简化为如下形式：

其中：ε为一个很小的常数。

同步骤41上式可以简化为：

令

则上式可以写为：

式中α^j和β^j分别表示A和B的第j行，θ_j表示与α^j对应的θ值，||α^j||₂为第j行α的2范数，||β^j||₂为第j行β的2范数，

令

b_j＝-2α^j(β^j)^T，则

上式可以分解为一系列标量最小化问题来进行求解：

令g(θ_j)＝a_jθ_j ²+b_jθ_j+clog(θ_j+ε)，最小化问题同样可以用一阶导数为0的方法进行求解，即

根据步骤41和步骤42所求得的B和θ可以得到X的估计值：

其中

和

分别为B和Λ的估计值。

步骤43，根据步骤41和42的方法求解一幅SAR图像的数学模型：

对于受相干斑影响的SAR图像，可以将其划分为N个不同的图像块集合，每个集合里面包含m个相似图像块，则一幅图像的稀疏编码优化问题可以表示为如下：

上式中X_j表示第j个图像块集合，

表示第j个图像块集合里的第i个图像块，||X_j-DΛ_jB_j||_F为误差的F范数，||x_i ^j||₂为

的2范数，||B_j||_F为B_j的F范数，η为正则化参数，根据经验进行设置。

根据上述模型，可以将图像的稀疏编码优化问题分解为N个图像块集合的优化问题，分别对这N个优化方程进行求解，最终整合起来得到整幅图像的稀疏编码最优解。具体实现流程如下：

①初始化：

1.设置图像的初始估计：

2.设置正则化参数η；

②外循环：迭代k＝1,2,…,k_max次

1.获取N个图像块集合{X_j}，计算每个X_j对应的字典基D_j并初始化θ_j、B_j；

2.内循环：迭代J＝1,2,…,J_max次

(I)对于固定的B_j，更新θ_j；

(II)对于固定的θ_j，更新B_j；

(III)使用B_j和θ_j对X_j进行重建；

(IV)根据更新后的X_j计算下一次迭代所需的θ_j和B_j；

循环结束

3.如果mod(k,k₀)＝0，为每个X_j更新字典基D_j；

循环结束

③输出

5、仿真结果

由于在实际工程应用中无法获取到不含噪的SAR图像，此处选择用合成场景的SAR图像以及真实场景的SAR图像分别对本发明所提出的方法进行实验验证，并将该方法与已经提出的较先进的SAR图像降斑算法进行比较，实验结果在如下部分给出。

合成场景SAR图像实验结果如下所示：

在合成场景SAR图像实验中，选取不同类型的3幅高分辨率光学图像，设置等效视数L为1、4、8以及16，分别对所选图像进行相干斑噪声的乘性添加，得到合成场景SAR图像。使用PSNR(Peak Signal to Noise Ratio，峰值信噪比)、SSIM(Structural Similarity，结构相似度)以及EPI(Edge Preserve Index，边缘保持指数)这三个指标对降噪算法进行性能评估，实验结果如下所示。表1为不同等效视数下各种算法进行降斑后所得的峰值信噪比，表2为相应的SSIM值。

表1中不同图像在不同的等效视数下经过降斑算法对图像进行重建之后所得到的最优PSNR用粗体字标出，PSNR越大证明降斑效果越好。通过比较表1中各算法降斑后的PSNR可以看出，Lee滤波算法作为经典降斑算法在性能上已经达到较好效果，滤波后PSNR相对于原始噪声图像提高了接近10个分贝，然而当前已经提出的较复杂的算法得到了更优的降斑效果。当L＝1时，图像受相干斑污染最严重，SAR-BM3D(三维转换域的非局部SAR图像降斑算法)算法所得到的PSNR在Lee滤波效果的基础上又增加了6～7个分贝。对于低视数的SAR图像，迭代PPB(基于块概率权重的降噪算法)算法比非迭代PPB(基于块概率权重的降噪算法)算法的降斑效果稍好，对于高视数的SAR图像，非迭代PPB的降斑效果更优。SAR-BM3D算法与FANS(快速非局部SAR降斑算法)算法的降斑性能几乎一致，均能得到较高的PSNR，在所比较的几种算法中性能最优。

从表1中可以看出，随着L的增大，本发明所提出的方法降噪能力越强，对于低PSNR的图像，由于图像被相干斑污染较严重，在稀疏编码过程中无法对稀疏系数进行较好的估计，因此降斑性能得不到很好的提高。当L较大时(L＞4)，本发明所提方法的PSNR可以达到PPB算法的效果，甚至优于PPB算法，但是相比于SAR-BM3D和FANS算法，依旧有轻微的差距。

表1各算法所得PSNR(dB)值

表2中列举了各算法对不同图像进行降斑之后计算所得的SSIM值，最优SSIM用粗体字标出。SSIM值越接近1，则表明降斑之后的图像越接近原图像。通过比较可得SAR-BM3D(三维转换域的非局部SAR图像降斑算法)与FANS(快速非局部SAR降斑算法)算法降斑之后的SSIM值相比其他算法更大，与PSNR所得的结果基本一致。对于不同视数的受污染图像，本发明所提出的方法所得的SSIM值优于PPB(基于块概率权重的降噪算法)算法(迭代以及非迭代)，这是由于在稀疏编码过程中尽可能多地保留了图像的细节信息，相比于SAR-BM3D以及FANS算法，SSIM值低0.02～0.03。

表2各算法所得SSIM值

为了进一步分析这几种方法的性能，表3给出了等效视数为1时，降斑后的图像相对于原图的边缘保持能力，EPI值越接近于1，则表明算法的边缘保持能力越强。由于Lee滤波算法相比其余算法的参考意义不大，因此表中只展示了当前降斑性能较优的几种算法的结果。通过比较表3所示的EPI值，可以看出当L＝1时，SAR-BM3D算法的边缘保持能力最强，其次是FANS算法，本发明所提出的方法的边缘保持能力略低于FANS，PPB算法最差，这与上述的SSIM值所得出的结果基本保持一致，且与降斑后的图像所呈现出来的视觉效果一致。

表3L＝1时各算法所得EPI值

	PPB non-it.	PPB 30it.	SAR-BM3D	FANS	本发明方法
						Lena	0.5019	0.5693	0.7184	0.6679	0.6172
Boat	0.4777	0.6260	0.7744	0.7300	0.6654
						House	0.4257	0.4499	0.7263	0.6308	0.7176

通过对上述合成场景SAR图像实验结果的分析可知，本发明所提出的方法在降斑性能上优于PPB算法，在降斑能力以及对图像的边缘特性与纹理特征的保持上均有明显优势，与目前性能较优的SAR-BM3D算法及FANS算法几乎没有差异。

真实场景SAR图像实验结果如下所示：

选取5幅真实场景SAR图像，这5幅图像包含了城市区域、农田、树木、河流这几类不同场景。分别使用已提出的PPB(迭代与非迭代)算法、SAR-BM3D算法、FANS算法以及本发明方法进行降斑处理，用ENL(Equivalent Number of Looks，等效视数)对降斑后的图像进行评估，ENL值越大表明相干斑抑制的效果越明显。

表4为降斑前后所选取的5幅图像的ENL值，表中最优的ENL用粗体字标注，可以看出使用本发明方法所得到的ENL值优于SAR-BM3D和FANS算法，PPB非迭代算法得到的降斑后的ENL值最大，对于相干斑的抑制效果最显著。根据重建后的图像效果可以看出PPB算法的降斑能力最强，能够很明显区分出图像的不同区域，但是对于异质区的降斑会出现过平滑现象；FANS算法的降斑能力优于SAR-BM3D算法，对细节的保持能力优于PPB算法；本发明方法在相干斑的抑制上优于FANS及SAR-BM3D，相比于PPB算法，在滤波过程中没有引入人工伪影，很好地保留了图像的纹理特征，且该方法原理简单，易于理解。综上所述，本发明提出的方法在对SAR图像的降斑处理上有很大的应用价值。

表4真实场景SAR图像降斑后的ENL值

	图1	图2	图3	图4	图5
						Noisy	8.5006	1.6435	21.0079	14.0650	6.0312
PPB non-it.	121.1494	10.7461	396.4247	325.1579	208.1072
						PPB 30it.	27.0814	3.3660	395.7035	317.5623	126.6741
SAR-BM3D	41.3188	3.3091	210.5339	104.1148	23.2862
						FANS	44.5532	2.8572	304.5839	165.0898	97.4936
本发明方法	80.1609	10.9697	241.3223	120.0762	161.2039

Claims (6)

1.一种基于高斯比例混合模型的稀疏表示SAR图像降斑方法，其特征在于：首先，建立单个图像块的稀疏表示模型；然后根据相干斑的统计特性与贝叶斯估计原理，将稀疏系数α用GSM模型进行表示，得到优化模型；同时对SAR图像进行分类，根据分类结果建立稀疏模型；最后，利用凸优化方法对优化模型求解，得到最优的稀疏表示，进而得到降噪图像。

2.根据权利要求1所述的一种基于高斯比例混合模型的稀疏表示SAR图像降斑方法，其特征在于，具体步骤如下：

步骤1，对单一图像块进行稀疏表示建模分析，给出最基本的凸优化数学模型；

步骤2，对步骤1中提出的数学模型进行分析，对稀疏系数α进行GSM建模，结合相干斑的统计特性与贝叶斯估计将该GSM模型代入到原有的凸优化模型中，由此转换求解域，得到更易于求解的凸优化方程；

步骤3，将SAR图像进行分类，得到相似图像块的集合，将步骤2中对于单一图像块建立的数学模型推广到图像块的集合中，得到图像块集合的稀疏表示模型；

步骤4，对步骤3中提出的稀疏表示模型，使用凸优化方法进行求解，得到最优解

分别对图像块的集合进行同样的求解过程，根据所得到的最优解对图像进行重建，实现噪声的滤除。

3.根据权利要求2所述的一种基于高斯比例混合模型的稀疏表示SAR图像降斑方法，其特征在于，所述步骤1中对单一图像块进行稀疏表示建模分析，给出最基本的凸优化数学模型的具体方法为：

步骤11，给出稀疏表示的基本数学模型：

假设图像块的大小为

像素，按照字典序排列形成列向量x∈Rⁿ，针对稀疏域建模，定义一个字典矩阵D∈R^n×k，其中k≥n，表明字典是冗余的，由此根据所提出的模型使用该字典稀疏表示每个图像块x，表示形式如下：

subject to Dα≈x

上式中α∈R^K为稀疏系数，||α||₀为系数α的0范数，

为利用凸优化工具得到的系数α的最优解；

步骤12，对步骤11中的基本数学模型进行等价变换：

使用l₁范数来替代原始的非凸优化问题，上述模型等价为：

其中λ为正则化参数，||α||₁为系数α的1范数，||x-Dα||₂为原图像与重建后图像的误差的2范数。

4.根据权利要求2所述的一种基于高斯比例混合模型的稀疏表示SAR图像降斑方法，其特征在于，所述步骤2的具体过程如下：

步骤21，实现稀疏系数的GSM建模：

使用GSM模型对稀疏系数α进行建模，将向量α分解为一个高斯向量β和一个标量乘性因子θ的逐点乘积，即α_i＝θ_iβ_i，其中θ_i表示概率为P(θ_i)的正标量，α_i为稀疏系数的一个元素，β_i为高斯向量的一个元素，假设θ_i是独立同分布的且与β_i无关，则α的GSM先验信息表示为如下：

其中：P(θ_i)为θ_i的概率，P(α_i|θ_i)为θ_i条件下α_i的概率；

步骤22，联合贝叶斯估计与相干斑的统计特性推导新的稀疏表示数学模型：

对于每一个图像块x∈Rⁿ，其稀疏表示数学模型写为如下形式：

x＝y·u＝y+y·(u-1)

＝y+v＝Dα+v

其中x表示观测图像块，y表示不含噪图像块，u表示相干斑，v＝y·(u-1)表示等效加性噪声，α为稀疏系数,D为字典；在振幅形式上，相干斑的统计特性服从Nakagami分布，概率密度函数表示为如下：

其中L表示等效视数，Γ(·)表示Gamma函数；由于y与u相互独立，且u-1的均值为0，因此v的均值为0；

根据贝叶斯准则，对于一个已知观测信号x＝Dα+v，表示为如下MAP估计：

(α,θ)＝arg max log P(x|α,θ)P(α,θ)

＝arg max log P(x|α)+log P(α|θ)+log P(θ)

其中：P(x|α,θ)表示在α和θ的条件下x的概率，P(α,θ)表示在θ的条件下α的概率，P(x|α)为α条件下x的概率，P(α|θ)为θ条件下α的概率，P(θ)为θ的概率；上式中，先验项P(α|θ)表示为如下形式：

结合相干斑的概率密度函数得到：

P(x|α)＝P(v)＝P(y·(u-1))＝P(y)·P(u-1)

＝C·(u-1)^2L-1exp(-L(u-1)²)

其中，

Γ(L)L的Gamma函数，P(v)为v的概率，P(y)为y的概率；记

选取

结合上述各项式子，推导出如下稀疏表示模型：

其中：||x||₂为x的2范数，使用log(θ_i+ε)代替logθ_i，ε为一个很小的正常数，记

为log(θ+ε)，上述公式简化为：

注意到原始的GSM模型的矩阵形式为α＝Λβ，其中Λ＝diag(θ_i)∈R^K×K为一个对角矩阵，表示被选图像块的方差域，R^K×K为K×K阶的实数域；因此，上述稀疏编码问题从α域转换到β域，形式如下：

上述稀疏表示模型即为单一图像块的GSM稀疏编码模型，其中||x-DΛβ||₂为误差的2范数，||β||₂为β的2范数。

5.根据权利要求2所述的一种基于高斯比例混合模型的稀疏表示SAR图像降斑方法，其特征在于，所述步骤3中将SAR图像进行分类，得到相似图像块的集合，并得到图像块集合的数学模型的具体方法为：

步骤31：将SAR图像进行相似图形块的分类：

使用概率估计来代替欧氏距离作为相似度的评价指标：

上式中，ω(y_i,y_j)是在平均过程中作为权重的相似度，表示了在这两个图像块中隐藏的降噪信号的概率是相同的；y_i和y_j分别表示第i个和第j个图像块的观测值；x_i和x_j分别表示第i个和第j个图像块的降噪值，Ω表示在图像块中的所有像素的集合；上述公式是一种最大似然估计且假设x_i(k)与x_j(k)条件独立，最终简化为：

将SAR图像相干斑的概率分布代入简化后的公式，得到最终的权重表达式：

其中，

z_i(k)为第i个图像块的观测值求根号，z_j(k)为第j个图像块的观测值求根号；选取S个权重值大于其他的图像块组成一个相似块的集合Φ_i，由此将一幅SAR图像分成相似图像块的集合；

步骤32，推导一个图像块集合的稀疏表示模型：

对于包含m个相似块的集合，如果考虑GSM模型的同步稀疏编码，根据步骤22中得到的单个图像块的稀疏表示模型得到群稀疏表达式：

其中，X＝[x₁,x₂…,x_m]表示m个相似块的集合，x₁为第一个图像块，x₂为第二个图像块，x_m为第m个图像块，||x_i||₂为x_i的二范数，A＝ΛB表示GSM模型的稀疏系数的群表示，其中B＝[β₁,β₂…,β_m]∈R^K×m表示稀疏系数对应的高斯向量的集合，β₁为第一个高斯向量，β₂为第二个高斯向量，β_m为第m个高斯向量，R^K×m为K×m阶的实数矩阵,||B||_F表示B的F范数，||X-DΛB||_F表示重建误差的F范数。

6.根据权利要求2所述的一种基于高斯比例混合模型的稀疏表示SAR图像降斑方法，其特征在于，所述步骤4中使用凸优化方法求解群稀疏编码模型的具体方法为：

步骤41，保证θ的值不变，求B：

当θ的值固定时，群稀疏表达式简化为如下形式：

由X＝DA，x_i＝Dα_i，D为正交矩阵，上式简化为如下：

其中：||A-ΛB||_F为误差的2范数，||α_i||₂为α_i的2范数；

令

则f(B)写为：

f(B)＝B^T(Λ^TΛ+dI)B-2B^TΛ^TA+A^TA

当函数f(B)为最小值时相应的B通过f(B)的一阶导数为0时来求解，即：

根据上式求得B的值为：

B＝(Λ^TΛ+dI)^-1Λ^TA

步骤42，保证B的值不变，求θ

当B的值固定时，群稀疏表达式简化为如下形式：

同步骤41上式简化为：

令

则上式写为：

式中α^j和β^j分别表示A和B的第j行，θ_j表示与α^j对应的θ值，||α^j||₂为第j行α的2范数，||β^j||₂为第j行β的2范数，

令

b_j＝-2α^j(β^j)^T，则

上式分解为一系列标量最小化问题来进行求解：

令g(θ_j)＝a_jθ_j ²+b_jθ_j+c log(θ_j+ε)，最小化问题同样用一阶导数为0的方法进行求解，即

根据步骤41和步骤42所求得的B和θ得到X的估计值：

其中

和

分别为B和Λ的估计值；

步骤43，根据步骤41和42的方法求解一幅SAR图像的数学模型：

对于受相干斑影响的SAR图像，将其划分为N个不同的图像块集合，每个集合里面包含m个相似图像块，则图像的稀疏编码优化问题表示为如下：

上式中X_j表示第j个图像块集合，

表示第j个图像块集合里的第i个图像块，||X_j-DΛ_jB_j||_F为误差的F范数，||x_i ^j||₂为

的2范数，||B_j||_F为B_j的F范数，η为正则化参数，根据经验进行设置。