本发明属于数字图像处理技术领域，具体涉及一种基于总曲率的SAR图像变分去噪方法。

背景技术：

：

相干斑噪声是合成孔径雷达(Synthetic Aperture Radar，简称SAR)图像的重要特征，严重影响SAR图像的可解译性。相干斑噪声通常作为乘性噪声来建模f＝uη，f为观察到的退化图像，u为原始清晰图像，η为噪声。理想的SAR图像去噪方法是能在去除斑点噪声的同时保留图像的边缘和细节信息。去噪方法为定义一个滤波器窗口估计相干斑的局部噪声方差，利用估计值进行滤波处理。现有技术中，常用滤波算法有均值滤波、中值滤波、局部滤波、Lee滤波、Lee-Sigma滤波、Frost滤波和Gamma-MAP滤波。研究表明，在均匀图像区域，这些滤波方法能够较好的削减噪声，而在非均匀区域，图像过于平衡或模糊，不能很好的保持边缘细节信息；现有技术中，全变分方法将图像去噪构建为能量函数的最小化问题，引入各项异性扩散方程，在平滑噪声的同时保持边缘。变分模型包括数据保真项和规则项，基于全变分TV规则项，AA模型是最早的SAR图像Gamma分布的乘性噪声去除模型，SST模型为Poisson分布的乘性噪声去除模型，DTDS模型为Rayleigh分布的乘性噪声去除模型，SO模型为综合的乘性噪声去除模型。常用的梯度下降法在求解乘性噪声去噪模型时存在计算速度慢的问题，通常在求解过程中引入一些快速算法包括Split Bregman算法，对偶算法等。TV规则项能够较好的保持边缘，但阶梯效应是其主要缺点，通常引入高阶项来克服这一缺点，主要包括Hessian矩阵规则项、拉普拉斯Δu规则项和总曲率规则项。由于计算的复杂性和非线性，目前乘性噪声变分模型还未引入高阶规则项，因此设计一种基于总曲率的SAR图像变分去噪方法，能够将待处理的图像既能平滑乘性噪声又能保图像边缘细节信息。

技术实现要素：

：

本发明的目的在于克服现有方法存在的缺陷，寻求设计提供一种基于总曲率的SAR图像变分去噪方法，该方法涉及的变分能量方程包括数据保真项和总曲率规则项，并且基于交替方向乘子法(Alternating Direction Method of Multipliers，缩写为ADMM)巧妙设计辅助变量，通过L2范数约束，实现能量方程最小化极值问题的求解，求解的图像既能平滑乘性噪声又能保图像边缘细节信息。

为了实现上述目的，本发明涉及的基于总曲率的SAR图像变分去噪方法的具体操作方法按照如下步骤进行：

a.选择一幅待处理的原始SAR图像f并根据该图像f建立基于总曲率的SAR图像变分去噪能量方程，对于输入的原始超声图像f，期望得到的去噪后的图像为u，基于总曲率的能量方程为：

其中，Ω为SAR图像区域，α为权重系数，a、b和c为成型噪声一阶项、平方项和对数项的参数；曲率规则项的使用能够在SAR图像去噪过程中更好的保持边缘细节信息；

b.对步骤a中所述的总曲率的能量方程进行转换，步骤a建立的能量方程的数据项和规则项均为非凸非线性，因此引入u＝ez进行变量替换，能量方程转换如下：

c.步骤b建立的能量方程具有高阶、非凸性，难以进行求解，引入分裂算子w、和q简化总曲率规则项，步骤b的能量方程形式化为带约束的极小值问题：

d.步骤c中所述的约束能够转换为和两个等价约束，设计约束因此，又被转化为和具有约束的变量是松弛的，至此，能量方程转化为可使用增广拉格朗日方法求解的方程：

e.步骤d中所述的约束w＝z、和采用L2惩罚项，由能够推导出因此使用L1惩罚约束这样极小化问题转换为以下子问题的交替优化问题：

其中，β1、β2、β3、β4和β5是正的惩罚参数，λ1、λ2、λ4和是拉格朗日乘子，能够根据相应规则更新；

f.利用变量交替迭代优化求解分别计算步骤e中的变量z,w,q,将步骤e的极小化问题转换为以下6个子问题：

g.分别求解步骤f中的ε1(z)、ε2(w)、ε4(q)、和的欧拉方程；ε2(w)的欧拉方程采用梯度降方法直接求解，和ε4(q)的欧拉方程采用广义软阈值公式求解，的欧拉方程能够直接采用投影方法，而ε1(z)和的欧拉方程为非线性，采用快速傅里叶变换方法进行求解；

h.对步骤g中的z,w,q,进行迭代求解，当相邻两次迭代的能量差小于设定的阈值时停止；

i.采用u＝ez得到的u即为去噪后的SAR图像。

本发明与现有技术相比，利用总曲率规则项进行SAR图像去噪，对于利用总曲率项建立的能量方程为了避免在求解时所产生的复杂运算，同时巧妙设计约束，引入辅助变量进行求解，不但提高了效率，而且减少了计算的复杂度，同时本发明提出的基于总曲率的SAR图像去噪方法具有非常好的实际应用价值，对于提高SAR图像的清晰度，提高图像的解译度起到了非常重要的作用，应用价值极高，市场前景广阔。

附图说明：

图1为本发明涉及的基于总曲率的SAR图像去噪方法流程图。

图2为本发明涉及的在图像SAR-1得到的结果与AA模型的比较，其中图2(a)为原始SAR-1图像，图2(b)为基于本发明α＝0.5得到的去噪结果图，图2(c)为基于本发明α＝1得到的去噪结果，图2(d)为基于AA模型α＝0.5得到的去噪结果图。

图3为本发明涉及的在图像SAR-2得到的结果与SST模型的比较，其中图3(a)为原始SAR-2图像；图3(b)为基于本发明α＝0.5得到的去噪结果图，图3(c)为基于本发明α＝1得到的去噪结果图，图3(d)为基于SST模型α＝0.5得到的去噪结果图。

图4为本发明涉及的在图像SAR-3得到的结果与DTDS模型的比较图，其中图4(a)为原始SAR-3图像，图4(b)为基于本发明α＝0.5得到的去噪结果图，图4(c)为基于本发明α＝1得到的去噪结果图，图4(d)为基于DTDS模型α＝0.5得到的去噪结果图。

图5为基于本发明开发的图像去噪应用程序，程序运行包括图像灰度值动态模拟和图像|能量值动态模拟两种方式，其中图5(a)为主界面图，图5(b)为图像灰度值动态模拟图，图5(c)图像|能量动态模拟图，图5(d)结果输出图。

图6为本发明涉及的图像SAR-1的灰度值时空变化结果图，其中图6(a)为原始SAR-1图像，图6(b)为基于本发明α＝0.5得到的去噪结果的灰度值三维图，图6(c)为基于本发明α＝1得到的去噪结果灰度值三维图，图6(d)为基于AA模型α＝0.5得到的去噪结果的灰度值三维图。

具体实施方式：

下面结合附图和具体实施方式对本发明做进一步说明：

实施例1：

本实施对SAR图像变分去噪时，具体操作方法按照如下步骤进行：

a.选择一幅待处理的原始超声图像f并根据该图像f建立基于总曲率的SAR图像变分去噪能量方程，对于输入的原始超声图像f，期望得到的去噪后的图像为u，基于总曲率的能量方程为：

其中，Ω为SAR图像区域，α为权重系数，a、b和c为成型噪声一阶项、平方项和对数项的参数；

b.对步骤a中所述的总曲率的能量方程进行转换，步骤a建立的能量方程的数据项和规则项均为非凸非线性，因此引入u＝ez进行变量替换，能量方程转换如下：

c.将步骤b建立的能量方程进行求解，引入分裂算子w、和q简化总曲率规则项，步骤b的能量方程形式化为带约束的极小值问题：

d.步骤c中所述的约束能够转换为和两个等价约束，设计约束因此，又被转化为和具有约束的变量是松弛的，至此，能量方程转化为可使用增广拉格朗日方法求解的方程：

e.步骤d中所述的约束w＝z、和采用L2惩罚项，由能够推导出因此使用L1惩罚约束这样极小化问题转换为以下子问题的交替优化问题：

其中，β1、β2、β3、β4和β5是正的惩罚参数，λ1、λ2、λ4和是拉格朗日乘子，能够根据相应规则更新；

f.利用变量交替迭代优化求解分别计算步骤e中的变量z,w,q,将步骤e的极小化问题转换为以下6个子问题：

g.分别求解步骤f中的ε1(z)、ε2(w)、ε4(q)、和的欧拉方程；ε2(w)的欧拉方程采用梯度降方法直接求解，和ε4(q)的欧拉方程采用广义软阈值公式求解，的欧拉方程能够直接采用投影方法，而ε1(z)和的欧拉方程为非线性，采用快速傅里叶变换方法进行求解；

h.对步骤g中的z,w,q,进行迭代求解，当相邻两次迭代的能量差小于设定的阈值时停止，具体迭代步骤如下：

(1)初始化参数z＝logf，w＝z，(β1,β2,β3,β4,β5,Δt,iteration)＞0，根据噪声分布函数确定参数a，b和c；

(2)固定wk，和求解ε1(z)的欧拉方程，采用快速傅里叶变换(Fast Fourier Transform,缩写为FFT)方法求z；欧拉方程为

上述方程能够形式化为如下方程：

步骤(2)中的欧拉方程的离散形式为：

(β5-β2(S1++S1--2I+S2++S2--2I))z(i,j)＝g(i,j),

将离散后的欧拉方程采用离散傅里叶变换(Discrete Fourier Transform，缩写为DFT)变换，得到如下方程：

由于dz＝(β1-2β3(coszi+coszj-2))＞0，采用反傅里叶变换求得z，其中

(3)固定zk+1和λ1k，求解ε2(w)的欧拉方程，先采用梯度降方法求w，其中欧拉方程为α(afe-w+bf2e-2w-c)+β1(zk+1-w)+λ1k＝0；在采用梯度降求解，具体方式如下：

(4)固定zk+1、λ2k和求解的欧拉方程，采用广义软阈值公式求其中欧拉方程的计算方式如下：

广义软阈值公式求解

(5)固定和λ4k.，求解ε4(q)的欧拉方程，采用广义软阈值公式求q；欧拉方程为广义软阈值公式求解

(6)固定qk+1,和λ4k，求解的欧拉方程，采用FFT求欧拉方程为

上述欧拉方程依据移位算子表达如下：

其中

再将上述移位算子表达后的方程采用FFT变换，得到如下方程：

系数为：

确保β4β5＞0,行列式D＝β5-2β4β5(coszi+coszj-2)即大于0，采用反傅里叶变换求得

(7)固定λ2k和求解的欧拉方程，采用投影法求

欧拉方程为

投影法求解

i.采用u＝ez得到的u即为去噪后的SAR图像。