0 引言

合成孔径雷达图像为许多应用提供有用信息，如遥感测绘、地面监测等.但由于成像传感器中大量随机分布的散射体反射的雷达回波相干叠加，从而不可避免地在SAR图像中产生相干斑乘性噪声^[1].由于斑点噪声的存在，视觉上观测图像质量下降，并严重削弱了其自动场景分割的性能.因此，散斑噪声去除是预处理步骤中一个关键的任务^[2-3].SAR图像相干斑噪声抑制算法主要是空域滤波算法、变换域滤波算法两类.其中空域滤波^[4-5]计算速度快，实时性好，但噪声平滑和图像细节保护受窗的选择影响大；变换域滤波是在各种变换域中结合系数的统计特性对图像进行抑斑工作，其中NSST变换结合了contourlet^[6]和curvelet^[6]各自的优点，具有局部特性、多尺度特性和各向异性的特点，并增强了shearlet^[7]变换的方向选择性和平移不变性，增强了细节结构的保持能力.

非局部自相似性最早在文献[8]中应用到图像去噪中，提出了非局部图像去噪算法，它能够较好地保留图像细节.Dabov提出的三维块匹配(block matching 3D, BM3D)算法^[9]同样基于非局部自相似性，是当前公认的最佳去噪方法之一.低秩矩阵逼近^[10-11]是一种矩阵秩最小问题，可等价为线性约束条件下的核范数最小化(nuclear norm minimization，NNM)^[8]问题.文献[12]证明低秩矩阵可以用NNM近似，NNM又可通过奇异值的软阈值法实现.但由于NNM算法在计算过程中将奇异值等同对待，这样在实际的图像处理中造成一定的偏差；文献[13]提出了WNNM的思路, 奇异值越大，权重越大，所占比重越大.该方法能够取得较高的峰值信噪比，且视觉效果很好.

SAR图像具有相似的结构特征和数据特征，其相似块构成的矩阵被认为类似低秩的.由于NSST变换具有良好的多分辨率特性和时频局部化特性，本文在NSST变换域中将WNNM算法和KAD^[14]扩展到SAR图像去噪中，进行NSST变换之后，高频分量含有图像边缘信息和大部分噪声，采用能够有效分离噪声和图像边缘的KAD算法进行去噪.低频分量包含图像的主要能量和少量噪声，采用WNNM算法能够有效去除噪声.本文首先使用小波鲁棒中值估计全局噪声方差^[14]，对输入的SAR图像进行NSST变换，然后对变换后的高频分量和低频分量分别采用KAD算法和WNNM算法进行去噪处理，最后，通过逆变换重构得到去噪后的图像，为了方便，把该方法记为NSST-WNNM.

1 相干噪声模型和NSST、WNNM、KAD理论 1.1 相干噪声模型

相干噪声是一种乘性噪声，即噪声水平随着目标后向散射强度的增强而增强，因此，乘性模型比较适合于对SAR图像的观测信号进行建模，即

$ Y(x, y)=X(x, y)N(x, y), $

(1)

其中：Y(x, y)表示强度格式的观测信号；X(x, y)表示被照射目标的雷达截面积；N(x, y)为衰落过程中所引起的相干斑噪声.并认为X(x, y)与N(x, y)是相互独立的随机过程，N(x, y)是独立同分布噪声图像，并且服从均值为1、方差为1/L的Gamma分布，L为等效视数，为了便于去噪处理，一般对式(1) 两边取对数变换，将乘性噪声转化为高斯噪声：

$ \ln \left( {Y\left( {x, y} \right)} \right) = \ln \left( {X\left( {x, y} \right)} \right) + \ln \left( {N\left( {x, y} \right)} \right). $

(2)

1.2 非下采样shearlet变换理论

非下采样的shearlet变换相对传统具有多方向性、平移不变性和稳定性等特点，NSST在对图像的线、边缘和纹理等特征的描述方面更为准确，具有较强的抗噪能力.当维数n=2时，具有合成膨胀的放射系统定义为:

$ {M_{\pmb{AB}}}\left( \psi \right) = \left\{ {{\psi _{i, j, k}}\left( x \right) = |\det \pmb{A}{\pmb{|}^{j/2}}\psi \left( {{\pmb{B}^l}{\pmb{A}^j}x-k} \right):j, l \in {\bf Z}, k \in {{\bf Z}^2}} \right\}, $

(3)

其中:ψ∈L²(R²); A和B是2×2可逆矩阵；|det B|=1.如果M_AB(ψ)具有式(4) 的紧框架条件，则M_AB(ψ)中的元素称为合成小波，A称为各向异性膨胀矩阵, Aⁱ与尺度变换相关联，B为剪切矩阵, B^j与保持面积不变的几何变换相关联.对∀f∈L²(R²)，有

$ \sum\limits_{j, l, k} {| < f, {\psi _{i, l, k}} > {|^2} = {{\left\| f \right\|}^2}.} $

(4)

NSST分解过程分为多尺度分解和多方向分解两部分.多尺度分解沿用非下采样金字塔分解方式加以实现，在经过k级多尺度分解后，每一幅源图像共可衍生出1幅低频微观子图像和k幅高频微观子图像，这些子图像均与源图像具有相同尺寸.NSST的多方向分解是将标准剪切滤波器从伪极化网络系统映射到笛卡尔坐标系统，完全摒弃了下采样环节，使图像的冗余度得到了很大的提高，图像系数冗余度的提高实现了自身的平移不变性，有效避免了类似于离散小波变换、轮廓波变换结果中的振铃效应^[15].图 1为非下采样shearlet变换的分解示意图.

1.3 加权核范数最小化理论

为了方便算法的描述，将公式(2) 的去噪模型转换为

$ y=x+n. $

(5)

其中:y是观测图像(噪声图像); x是需要恢复的清晰图像; n是噪声方差为σ_n的高斯噪声.对观测图像y进行分块，设y_j为其中的第j块，可以通过一些方法如块匹配^[16]来搜索整个图像的非局部相似块.通过堆叠这些非局部的相似块形成矩阵，设为Y_j，令Y_j=X_j+N_j，其中X_j和N_j分别表示清晰图像和噪声相应的矩阵.由先验知识可知，SAR图像具有冗余性和相似性，X_j认为近似低秩^[13].因此X_j可以利用WNNM低秩矩阵恢复算法进行低秩矩阵恢复^[15], 加权核范数最小化可以描述为问题：

$ {{\pmb{\hat X}}_j} = \mathop {\min }\limits_{{\pmb{X}_j}} \left\| {{\pmb{Y}_j}-{\pmb{X}_j}} \right\|_F^2 + {\left\| {{\pmb{X}_j}} \right\|_{w, *}}, $

(6)

其中: ${\left\| {{\pmb{X}_j}} \right\|_{w, *}} = \sum\limits_i {{{\left| {{w_i}{\sigma _i}\left( {{\pmb{X}_j}} \right)} \right|}_1};w = \left[{{w_1}, \cdots, {w_n}} \right]} $和w_i≥0是分配到σ_i(X_j)的一个非负权重.式(6) 通过奇异值阈值分解获得全局最优解：

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\mathit{\boldsymbol{U \boldsymbol{\varLambda} }}{\mathit{\boldsymbol{V}}^{\rm{T}}} = svd\left( {{{\mathit{\boldsymbol{\hat X}}}_j}} \right),\\ {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_w} = {\mathit{\boldsymbol{S}}_w}\left( \mathit{\boldsymbol{ \boldsymbol{\varLambda} }} \right) = \max \left( {{\sigma _i}\left( {{\mathit{\boldsymbol{X}}_j}} \right) - w,0} \right), \end{array} $

(7)

其中权重w_i的计算为

$ {w_i} = c\sqrt q /\left( {{\sigma _i}\left( {{\pmb{X}_j}} \right) + \varepsilon } \right), $

(8)

其中:c>0是一个常数; q是Y_j中相似块的数量; ε=10^-16是避免被零整除的数.

假设噪声是均匀分布在U和V两个子空间中，则σ_i(X_j)可以估计为

$ {{\hat \sigma }_i}\left( {{\pmb{X}_j}} \right) = \sqrt {\max \left( {\sigma _i^2\left( {{\pmb{Y}_j}} \right)-n\sigma _n^2, 0} \right)}, $

(9)

其中：σ_i(Y_j)是Y_j的第i个奇异值，σ_n²为噪声方差.

1.4 核各向异性扩散去噪算法理论^[14-15]

KAD算法的偏微分方程为：

$ \left\{ \begin{array}{l} \frac{{\partial I}}{{\partial t}} = {\rm{div}}\left[{c\left( {\left\| {\nabla \left[{\phi \left( I \right) \cdot \nabla I} \right]} \right\|} \right)} \right], \\ I\left( {t = 0} \right) = {I_0}, \end{array} \right. $

(10)

其中：div为散度算子；∇为梯度算子；t为时间；‖∇ϕ(I)‖为映射到高维特征空间后的梯度模.

通过核函数k(x, y)映射后，高维特征空间中某一特定方向上的灰度差分为

$ \left\| {\phi \left( {{I_p}} \right)- \phi \left( {{I_q}} \right)} \right\| = {\left[{k\left( {{I_p}, {I_p}} \right) + k\left( {{I_q}, {I_q}} \right)-2k\left( {{I_p}, {I_q}} \right)} \right]^{0.5}}, $

(11)

则高维特征空间的梯度模为

$ {\left\| {\nabla \phi \left( I \right)} \right\|_p} = {\left\{ {\frac{1}{{\left| {{\xi _p}} \right|}}\sum\limits_{q \in {\xi _p}}^{} {\left[{k\left( {{I_p}, {I_p}} \right) + k\left( {{I_q}, {I_q}} \right)-2k\left( {{I_p}, {I_q}} \right)} \right]} } \right\}^{0.5}}, $

(12)

式中，ξ_p为像素p的8邻域集，|ξ_p|为集的势，本文选用的扩散方程为

$ c = \left( {\left\| {\nabla \phi \left( I \right)} \right\|} \right) = \exp \left[{-{{\left( {\left\| {\nabla \phi \left( I \right)} \right\|/k} \right)}^2}} \right], $

(13)

式中，k为扩散门限，采用绝对偏差中值(median absolute deviation, MAD)，即

$ k = \frac{1}{{0.674\;5}}{f_{{\rm{med}}}}\left\{ {\left\| {\nabla \phi \left[I \right] - {f_{{\rm{med}}}}\left( {\left\| {\nabla \phi \left[I \right]} \right\|} \right)} \right\|} \right\}, $

(14)

式中，f_med表示取中值函数.

2 基于NSST变换域的WNNM和KAD的SAR图像去噪算法

在首次对图像进行NSST变换之前，为了确定初次分块操作的参数和窗口大小，首先采用小波鲁棒中值估计方法^[16]对SAR图像进行噪声预估计，

$ \sigma _n^2 = \frac{{Median(\left| \pmb{Y} \right|)}}{{0.674\;5}}, \pmb{Y} \in subband\left( {H{H_1}} \right), $

(15)

通过式(15) 得到SAR图像的预估计噪声，之后根据噪声方差的大小用经验值确定参数.采用NSST变换获得一系列多尺度多方向下的子带图像；对低频子带采用WNNM算法进行去噪处理，得到去噪后的系数，对6个高频子带采用KAD算法；最后进行NSST反变换获得最终结果.

本文提出的基于NSST变换域WNNM算法的SAR图像去噪过程如图 2所示.去噪步骤如下：

步骤1  对SAR图像进行对数变换，得到图像y；

步骤2  对y采用非下采样拉普拉斯金字塔算法进行多尺度分解，得到低频分量和各个尺度下的高频分量，在本文中，我们使用单层NSST变换.

步骤3  采用带方向和尺度变化的窗函数对高频分量进行方向剖分，得到多个频率分量，之后进行坐标映射，获得高频分量的shearlet系数.

步骤4  将WNNM应用到低频分量去噪，得到去噪后的低频分量，能够尽可能地保留图像细节特征.

步骤5  由于高频分量中含有图像边缘信息和大部分噪声，采用KAD算法分别对6个高频系数进行去噪，得到去噪后的高频分量.

步骤6  用NSST逆变换逐层重构得到去噪图像.

3 实验结果与分析

为了测试算法的真实效果，我们将算法运用到SAR图像上进行相干斑去噪，图 3和图 4中使用的方法(a)~(d)，依次是KAD去噪、NSST变换去噪、WNNM去噪算法和本文算法.

从实验结果可以看出：图 3(a)和图 4(a)中单独采用核各向异性算法噪声抑制不彻底，存在残留噪声的现象，且边缘比较模糊；图 3(b)和图 4(b)中的NSST变换滤波因为具有平移不变性且方向选择性较好，因此在图中残留的噪声较少，去噪后的图像也较为清晰，但由于采用的阈值函数具有不连续性，仍存在伪吉布斯效应；图 3(c)和图 4(c)的WNNM算法噪声残留较少，去噪后图像具有良好的视觉效果，但是背景过于平滑；本文方法采用的NSST具有平移不变性，弥补了剪切波的不足，KAD算法对于分离噪声和边缘具有很好的效果，结合WNNM算法分别在高频子带和低频子带中使用，使去噪效果更优.

为了更好地展示本文所提到的算法的优越性，本文计算了各种去噪算法的几种常用的去噪性能参数，其中包括峰值信噪比(peak signal to noise ratio，PSNR)^[17]、等效视数(equivalent numbers of looks, ENL)^[18]和边缘保持指数(edge preservation index, EPI)^[18]，其中PSNR越大表明算法的去噪能力越强，ENL越大表明算法的去噪后的视觉效果越好，而EPI越大表明算法的边界保持能力越强.表 1给出了图 3和图 4中的3种去噪算法对原始SAR图像进行去噪后的PSNR、ENL和EPI.从表 1所示的各种去噪方法的性能参数来看基于NSST变换域的WNNM算法的SAR图像去噪算法是一种较好的去噪算法.而在ENL和EPI方面，本文算法也有了明显的提高，这表明本文算法去噪的视觉效果更好，边界保持能力更强.NSST-WNNM算法相较于WNNM算法，PSNR提高了0.3~0.4 dB，相较于NSST变换算法，PSNR提高了1.2~1.5 dB，相较于KAD算法，PSNR提高了3 dB, NSST-WNNM算法相较于其他3种算法，ENL和EPI有了些许提高.综合上面的分析，可知本文的算法对SAR图像的去噪不仅去噪能力很强，且去噪后的视觉效果也很好，能更好保留SAR图像的纹理信息.

4 结论

本文提出一种基于非下采样shearlet变换域WNNM和KAD算法的SAR图像去噪，该算法利用NSST变换的多分辨率特性、时频局部化特性、KAD算法在边缘和噪声中较强的分辨性能以及WNNM算法处理图像过程中的全局特性.首先分析了相干斑噪声的模型，对SAR图像进行对数变换，将相干噪声转化为高斯噪声，其次对SAR图像进行NSST变换波分解，对分解后的低频分量用WNNM算法进行去噪处理，高频分量采用KAD算法进行处理，最后用处理后的结果进行NSST重构得到去噪图像，生成相干斑抑制后的SAR图像.实验结果表明本文算法能有效地抑制斑点噪声、同时又能较好保持图像的纹理细节和边缘信息，能够满足后续图像处理任务的要求.