﻿ 基于区域分割的非局部全变差SAR相干斑滤波
 舰船科学技术  2016, Vol. 38 Issue (5): 105-110 PDF

Non-local total variation speckle filtering based on region segmentation
FAN Hui-ling, QU Chang-wen, LI Jian-wei
Department of Electronic and Information Engineering, Naval Aeronautical and Astronautical University, Yantai 264001, China
Abstract: Due to the limit of imaging system, synthetic aperture radar images are often corrupted with non-Gaussian multiplicative speckle noise. For the purpose of noise suppression, an algorithm called NLM-TV is proposed, which integrates non-local means filter and total variation regularization. It involves three steps. First, convert the multiplicative noise into signal-dependent additive noise, dividing the image into three categories based on noise level, the edge, the strong scattering region and the weak scattering region. Non-local means filter was then applied. To maintain the edge structure effectively, the smooth parameter must be small. In the strong scattering region, TV regularization was used because of insufficient reduction of speckles. In this paper, several simulations were conducted in RADARSAT-2 and TerraSAR-X images, the results showed that compared to a variety of filtering algorithms, NLM-TV algorithm could significantly increase the equivalent number of looks both in the weak scattering area and strong scattering regions. At the same time, edge keeping index could be increased by more than 10%.
Key words: SAR images     speckle     non-local means filter     total variation regularization
0 引 言

1 加性噪声模型与区域划分 1.1 依赖于信号的加性噪声

L 视 SAR 图像 I 可表示为：

 $I=R\odot \odot .$ (1)

 ${{I}_{i,j}}={{R}_{i,j}}+{{R}_{i,j}}({{S}_{i,j}}-1).$ (2)

 \left\{ \begin{align} &E[R]E[N]=0\text{,} \\ &E[RN]=E[{{R}^{2}}]E[(N-1)]=0. \\ \end{align} \right. (3)

 \left\{ \begin{align} &E[I]=E[R]\text{,} \\ &D[I]=D[R]+D[N]=D[R]+E[{{R}^{2}}]D[S]. \\ \end{align} \right. (4)

 $E[{{I}^{2}}]=E[{{R}^{2}}](1+D[S]).$ (5)

1.2 区域划分

 $\overline{I_{i,j}^{2}}=\frac{1}{\left| {{\Omega }_{ij}} \right|}\sum\limits_{(p,q)\in {{\Omega }_{ij}}}^{{}}{{}}I_{p,q}^{2}.$ (6)

1）为削弱相干斑对边缘提取的影响，将图像 I 作均值滤波，窗口大小视情况而定，一般设为 $5 \times 5$$7 \times 7，由 “canny” 算子提取 SAR 图像的纹理与边缘。记像元索引集合为 {\Lambda _1} 2）h 算子将 I 映射为 G，窗口大小与均值滤波一致。从 G 中剔除属于 {\Lambda _1} 索引的像元，得到 {{G'}} = {{G}}/{\varOmega _1}。对 {{G'}} 进行最佳直方图分割，得到强散射区像元集合 {\Lambda _2} 和弱散射区像元集合 {\Lambda _3} 2 NLM-TV 降斑算法 2.1 非局部均值滤波 NLMF 最先由 Buades 等在文献[7]中提出，其突出贡献打破了传统的局部邻域框架。记 {\mathbf{\tilde I}} 为滤波后的图像，有  {{\tilde{I}}_{i,j}}=\sum\limits_{(p,q)\in \Theta }{W_{i,j}^{p,q}{{I}_{p,q}}}, (7) 式中：\Theta$${I_{i,j}}$ 的相似元素集合；$W_{i,j}^{p,q}$ 为依赖于 $(i,j)$$(p,q) 的权重。  W_{i,j}^{p,q}=\frac{1}{Z}{{e}^{-\frac{\left\| {{\delta }_{i,j}}-{{\delta }_{p,q}} \right\|_{\text{F}}^{2}}{h}}}. (8) 式中：{{{\delta }}_{{{i,j}}}},{{{\delta }}_{{{p,q}}}} 分别为 (i,j),(p,q) 为中心的邻域块；{\left\| {} \right\|_{\text{F}}} 为矩阵的 Frobenius 范数；Z 为归一化因子；h 为平滑参数。 NLMF 通过设置搜索窗口与相似窗口，寻找搜索窗口内的相似块，根据相似程度进行加权处理，其核心思想是利用冗余信息达到降噪目的。平滑参数 h 控制权重的衰减程度，h 越大获得的图像越平滑，h 越小图像结构信息保持越好。 2.2 离散 TV 模型 TV 模型作为一种正则化方法广泛用于图像处理之中，在一定程度上可以起到保持边缘、图像平滑的作用。式（9）为离散 TV 降噪模型  {{\tilde{I}}_{*}}=\arg \min \left\| \tilde{I}-I \right\|_{\text{F}}^{2}+\lambda \text{TV}(\tilde{I}) (9) 式中：\lambda 为大于 0 的正则化参数；对任意的 {{\tilde I}} \in {R^{m \times n}}$${\text{TV(}} \cdot {\text{)}}$ 为离散全变差算子，

 $\text{TV}(\tilde{I})=\sum\limits_{i=1}^{m}{\sum\limits_{j=1}^{n}{(\left| (\nabla \tilde{I})_{i,j}^{x} \right|+\left| (\nabla \tilde{I})_{i,j}^{y} \right|)}},$ (10)

 \begin{align} &\left\{ \begin{array}{*{35}{l}} (\nabla \tilde{I})_{i,j}^{x}={{{\tilde{I}}}_{i,j}}-{{{\tilde{I}}}_{i,j+1}} &j(\nabla {{\tilde I}})_{i,j}^x$${{\tilde I}}$$(i,j)$处的水平方向差分，$(\nabla {{\tilde I}})_{i,j}^y$${{\tilde I}}$$(i,j)$处的垂直方向差分。 正则化参数的选择在 TV 模型起着均衡截断误差$\left\| {{{\tilde I - I}}} \right\|_F^2$和平滑项${\text{TV}}({{\tilde I}})$的重要作用。$\lambda \to 0$时，截断误差$\left\| {{{\tilde I - I}}} \right\|_F^2$成为控制式（9）的主要因素；$\lambda \to + \infty $时，式（9）的解${{{\tilde I}}_{{*}}}$为一幅充分平滑的图像。考虑到本文只对强散射区${\Lambda _2}$进行平滑，式（9）改为 ${{\left. {{{\tilde{I}}}_{*}} \right|}_{{{\Lambda }_{2}}}}=\arg \min {{\left. \left\| \tilde{I}-I \right\|_{\text{F}}^{2} \right|}_{{{\Lambda }_{2}}}}+\lambda {{\left. \text{TV}(\tilde{I}) \right|}_{{{\Lambda }_{2}}}}.$(12) 式中${\left. \cdot \right|_{{\Lambda _2}}}$为在像元集合${\Lambda _2}$上进行操作。 2.3 NLM-TV 算法流程 基于以上分析，将 SAR 图像先进行 NLMF 处理后，对降噪不完全的强散射区再进行 TV 平滑，我们称之为 NLM-TV 降斑算法。具体步骤如下： 1）将 SAR 图像噪声视为加性，按照 1.2 节相关步骤完成区域划分； 2）根据边缘与弱散射区的像素，计算边缘与弱散射区的 NLMF 平滑参数${h_1},{h_2}$3）对图像进行 NLMF 处理，边缘区域使用参数${h_1}$，其余区域使用参数${h_2}$4）根据处理后图像与强散射区索引，计算 TV 平滑参数$\lambda $，对强散射区进行 TV 平滑处理； 5）融合边缘弱散射区与强散射区的处理结果，得到最终降斑图像。 具体流程如图 1 所示。  图 1 NLM-TV 降斑流程图 Fig. 1 NLM-TV speckle suppression flow chart 3 相关参数选取与迭代求解算法 3.1 平滑参数 h 的选取 平滑参数 h 决定了相似度函数的衰减速率，h 越大，相似度函数的衰减速率越慢，相似图像块的权重较大，获得的图像更加平滑；相反，图像纹理越清晰。h 的大小应与图像的噪声水平相关，杨学志等[12]通过研究得出，当 h 介于噪声方差的 0.4 ~ 0.6 倍之间时，降噪效果较好。在本文中，边缘处取比例系数为 0.4,其余区域取比例系数为 0.6。 3.2 正则化参数 λ 的选取 目前主要的正则化参数选择方法有相容性准则，无偏预风险估计、L曲线及广义交叉校验准则等，具体参见文献[14]。本文采用广义交叉校验准则求取正则化参数 λ，对于文中 TV 模型，λ 的选取可由式（13）确定： ${{\lambda }^{*}}=\min \ \frac{(\lambda +1){{\left. {{\left\| {{G}^{*}}-I \right\|}_{\text{F}}} \right|}_{{{\Lambda }_{2}}}}}{\lambda \left| {{\Lambda }_{2}} \right|}.$(13) 3.3 TV 模型的迭代求解算法 在求解大规模离散正则化问题时，须迭代的收敛性问题。文献中，Amir Beck [18]提出的一种快速投影梯度（Fast projected gradient，FPG）算法解决了 TV 的正则化降噪问题，并且证明该方法具有二阶收敛速率，可以很大程度上节省计算资源。$\varPhi $为矩阵对$({{B,C}})$的集合，其中${{B}} \in {R^{(m - 1) \times n}}{{C}} \in {R^{m \times (n - 1)}}，满足：  \left\{ \begin{align} &B_{i,j}^{2}+C_{i,j}^{2}1\quad i=1,\cdots ,m-1,j=1,\cdots ,n-1, \\ &\left| {{B}_{i,n}} \right|1\quad i=1,\cdots ,m-1, \\ &\left| {{C}_{m,j}} \right|1\quad j=1,\cdots ,n-1; \\ \end{align} \right. (14) 线性算子： \kappa :{R^{(m - 1) \times n}} \times {R^{m \times (n - 1)}} \to {R^{m \times n}}，满足：  \begin{align} &{{\left[ \kappa (\mathbf{B},\mathbf{C}) \right]}_{i,j}}={{B}_{i,j}}+{{C}_{i,j}}-{{B}_{i-1,j}}-{{C}_{i,j-1}}, \\ &i=1,\cdots ,m,j=1,\cdots ,n. \\ \end{align} (15) 式中：\forall i \,=\,1,\cdots ,m;j \,=\,1,\cdots ,n{B_{0,j}} \,=\,{B_{m,j}} = {C_{i,0}}= $${C_{i,n}} = 0 线性算子 {\kappa ^\text{T}}$$\kappa $的伴随算子，${\kappa ^\text{T}}:{R^{m \times n}} \to {\kappa ^\text{T}}:{R^{(m - 1) \times n}} \times {R^{m \times (n - 1)}}，满足：  \left\{ \begin{align} &{{\kappa }^{\text{T}}}(I)=(B,C) \\ &{{B}_{i,j}}={{I}_{i,j}}-{{I}_{i+1,j}}\quad i=1,\cdots ,m-1;j=1,\cdots ,n, \\ &{{C}_{i,j}}={{I}_{i,j}}-{{I}_{i,j+1}}\quad i=1,\cdots ,m;j=1,\cdots ,n-1. \\ \end{align} \right. (16) {{ P}_\varPsi }( \cdot ) 为在集合 \varPsi  上的投影算子，满足：  {{P}_{\Psi }}(I)=\left\{ \begin{align} &l\quad \quad {{I}_{i,j}}u. \\ \end{align} \right. (17) 其中 u,l 为集合 \varPsi  上的最大元素与最小元素，算法中 \varPsi  应设为 SAR 图像的动态范围。 表 1 给出了 FPG 算法的实现，具体参见文献[18] 表 1 FPG 算法 Tab.1 FPG algorithm 4 数值仿真与算法对比 以 RADARSAT-2，TerraSAR-X 实测数据对所提算法进行仿真验证，并与 SRAD，NLMF 及非下采样 contourlet 变换（Nonsubsampled Contourlet Transform，NSCT）降噪算法进行比较，通过等效视数（Equivalent number of looks，ENL）、峰值信噪比（Peak signal to noise ratio，PSNR）与边缘保持指数（Edge keeping index，EKI）这 3 个指标对不同降斑效果进行评价。 4.1 算法仿真 图 2（a）为加拿大 RADARSAT-2 于 2008 年 4 月 19 日获取的加拿大温哥华市的部分图像实验数据，HV 极化，256 × 256 像素；图 2（b）为德国 TerraSAR-X 于 2007 年 8 月 10 日获取的德国诺林根城市的部分图像，HH 极化，7 × 7 像素。  图 2 SAR 图像实验数据 Fig. 2 SAR image experimental data 提取图像的纹理与边缘，并将其划分为弱散射区、强散射区及纹理 3 个区域，分别对应图 3 中的黑、灰、白 3 种灰度，均值滤波窗口 7 × 7。  图 3 SAR 图像区域划分结果 Fig. 3 SAR image segmentation result 图像经标准化处理后，对两景图像使用 NLM-TV 降噪，NLMF 滤波时取搜索窗大小 11 × 11，相似窗大小 5 × 5。RADARSAT-2 图像，弱散射区平滑参数 0.085，边缘平滑参数 0.057，正则化参数为 0.104；TerraSAR-X 图像，弱散射区平滑参数 0.088，边缘平滑参数 0.059，正则化参数 0.093。SRAD 滤波处理，迭代步长 0.05，RADARSAT-2 图像迭代 140 次，TerraSAR-X 图像迭代 180 次。 经对数变换后，NLMF 滤波的处理，取搜索窗大小 11 × 11，相似窗大小 5 × 5。RADARSAT-2 图像平滑参数 0.077，TerraSAR-X 图像平滑参数 0.080。基于 NSCT 降噪]处理结果，分解层数为 3，各层子带数均为 8，估计 RADARSAT-2 图像噪声方差为 0.034，TerraSAR-X 图像噪声方差为 0.039。 4.2 算法对比与性能分析 为了更客观地评价降斑效果，引入 ENL，PSNR 和 EKI 三个指标，对 SAR 图像降斑效果量化分析。图 2（a）图 2（b）中标出了 4 个区域（①，②，③，④）。①，② 区域为匀质区域，用于计算 ENL 与 PSNR；③，④ 区域包含了图像的边缘，用于计算边缘保持指数。 ENL 用来衡量图像斑点噪声的相对强度，其定义式为：  ENL=\frac{{{\mu }^{2}}}{{{\sigma }^{2}}}. (18) 式中 \mu\sigma $分别为匀质区域内统计均值与标准差。PSNR 衡量背景杂波中目标的凸显程度，定义为： $PSNR=10{{\log }_{10}}\left( \frac{\underset{(i,j)\in \Omega }{\mathop{\max }}\,\left( \tilde{I}_{i,j}^{2} \right)}{\frac{1}{\left| \Omega \right|}\sum\limits_{(i,j)\in \Omega }{{{({{I}_{i,j}}-{{{\tilde{I}}}_{i,j}})}^{2}}}} \right).$(19) 式中：$\mathop {\max }\limits_{(i,j) \in \varOmega } ( \cdot )$为区域$\varOmega $内的最大值；I 为原始图像；${{\tilde I}}$为处理后图像。 表 2 列出了不同算法降斑处理后，两景图像 ①，② 区的 ENLPSNR 表 2 不同处理算法 ①，② 区的 ENL 与 PSNR Tab.2 The ENL and PSNR of different methods 分析表 2 数据，SRADNLMFNSCTNLM-TV 均能够对滤除相干斑，平滑图像起到一定的作用。在弱散射区，基于 NSCT 的滤波算法能够获得最高的 ENL，但是其代价是牺牲 PSNR。NLMF 与 NSCT 滤波在弱散射区对 ENL 提升倍数达到 13 倍以上，但在强散射区对 ENL 提升仅为 2 倍左右，原因是二者兼以加性噪声模型为基础，第一节的分析表明，当乘性噪声被当成加性噪声处理时，因噪声依赖于信号，导致了不同区域降斑效果不尽相同。SRAD 滤波通过揭示 Lee 滤波与偏微分方程滤波的内在联系，是一种以乘性噪声为基础的滤波算法，因此对强散射区，弱散射区的相干斑抑制能力较为均衡，在弱散射区与强散射区都能对 ENL 进行较大提升。NLM-TV 是 NLMF 与 TV 的融合，相对于 NLMF 具有明显的优势，具体表现为在弱散射区能够获得更高的 ENLPSNR，在强散射区 ENL 提升较高，但是 PSNR 下降不明显，这表明该算法能够很大程度的滤除相干斑噪声，并保证目标不在降噪过程中被平滑处理而淹没。从表 2ENLPSNR 数据来看 SRAD 滤波与本文算法是 4 个算法中较为优秀的 2 个。 EKI 用来衡量经降斑处理后的边缘和纹理保持程度，EKI 定义为： $EKI=\frac{\sum\limits_{i=1}^{N}{T(I,{{w}_{i}})}}{\sum\limits_{i=1}^{N}{T(\tilde{I},{{w}_{i}})}}.$(20) 式中：$T( \cdot ,{w_i})$为图像窗口${w_i}\$ 内梯度的最大值；N 为窗口数目；EKI 越接近 1，证明边缘保持效果越好。

5 结 语

 [1] 焦李成, 张向荣, 侯彪, 等. 智能SAR图像处理与解译[M]. 北京: 科学出版社, 2008 : 43 -53. [2] GOODMAN J W. Some fundamental properties of speckle[J]. Journal of the Optical Society of America , 1976, 66 (11) :1145–1150. DOI:10.1364/JOSA.66.001145 [3] LEE J S. A simple speckle smoothing algorithm for synthetic aperture radar images[J]. IEEE Transactions on Systems, Man, and Cybernetics , 1983, SMC-13 (1) :85–89. DOI:10.1109/TSMC.1983.6313036 [4] KUAN D, SAWCHUK A, STRAND T, et al. Adaptive restoration of images with speckle[J]. IEEE Transactions on Acoustics, Speech, and Signal Processing , 1987, 35 (3) :373–383. DOI:10.1109/TASSP.1987.1165131 [5] YU Y J, ACTON S T. Speckle reducing anisotropic diffusion[J]. IEEE Transactions on Image Processing , 2002, 11 (11) :1260–1270. DOI:10.1109/TIP.2002.804276 [6] PERONA P, MALIK J. Scale-space and edge detection using anisotropic diffusion[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence , 1990, 12 (7) :629–639. DOI:10.1109/34.56205 [7] BUADES A, COLL B, MOREL J M. A non-local algorithm for image denoising[C]//IEEE Computer Society Conference on Computer Vision and Pattern Recognition. San Diego, USA:IEEE, 2005, 2:60-65. [8] ZHONG H, LI Y W, JIAO L C. SAR image despeckling using Bayesian nonlocal means filter with sigma Preselection[J]. IEEE Geoscience and Remote Sensing Letters , 2011, 8 (4) :809–813. DOI:10.1109/LGRS.2011.2112331 [9] PARRILLI S, PODERICO M, ANGELINO C V, et al. A non-local SAR image denoising algorithm based on LLMMSE wavelet shrinkage[J]. IEEE Transactions on Geoscience and Remote Sensing , 2012, 50 (2) :606–616. DOI:10.1109/TGRS.2011.2161586 [10] JOJY C, NAIR M S, SUBRAHMANYAM G R K S, et al. Discontinuity adaptive non-local means with importance sampling unscented Kalman filter for despeckling SAR images[J]. IEEE Journal of Selected Topics in Applied Earth Obersvations and Remote Sensing , 2013, 6 (4) :1964–1970. DOI:10.1109/JSTARS.2012.2231055 [11] 郑永恒, 程建, 曹宗杰. 改进非局部均值滤波的SAR图像降噪[J]. 中国图象图形学报 , 2012, 17 (7) :886–891. [12] 杨学志, 沈晶, 范良欢. 基于非局部均值滤波的结构保持相干斑噪声抑制方法[J]. 中国图象图形学报 , 2009, 14 (12) :2443–2450. [13] DA CUNHA A L, ZHOU J, DO M N. The nonsubsampled contourlet transform:Theory, design, and applications[J]. IEEE Transactions on Image Processing , 2006, 15 (10) :3089–3101. DOI:10.1109/TIP.2006.877507 [14] 常霞, 焦立成, 刘芳. 基于斑点方差估计的非下采样Contour-let域SAR图像去噪[J]. 电子学报 , 2010, 38 (6) :1328–1333. [15] 贾建, 陈莉. 基于双变量模型和非下采样Contourlet变换的SAR图像相干斑抑制[J]. 电子与信息学报 , 2011, 33 (5) :1088–1094. [16] WEN Y W, Chan R H. Parameter selection for total-variation-based image restoration using discrepancy principle[J]. IEEE Transactions on Image Processing , 2012, 21 (4) :1770–1781. DOI:10.1109/TIP.2011.2181401 [17] SHEN J H, CHAN T. Variational restoration of nonflat image features:Models and algorithms[J]. SIAM Journal on Applied Mathematics , 2001, 61 (4) :1338–1361. DOI:10.1137/S003613999935799X [18] BECK A, TEBOULLE M. Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[J]. IEEE Transactions on Image Processing , 2009, 18 (11) :2419–2434. DOI:10.1109/TIP.2009.2028250