基于Contourlet域隐马尔可夫树模型的SAR图像滤波方法

邓 磊，李家存，朱佳文，孙 萍

（首都师范大学 资源环境与旅游学院，北京100048）

合成孔径雷达（synthetic aperture radar，SAR）具有全天候、全天时、多波段、多极化等特点，已成为地球表面信息的重要信息来源.斑点噪声是由SAR所固有的特性造成的，它的存在严重影响了SAR图像的应用.因此，对斑点噪声抑制技术的研究一直是一项重要课题［1］.

一般来说，滤波技术可分为成像前多视处理［2］和成像后滤波.前者以降低空间分辨率为代价来获得斑点噪声的抑制效果，后者又主要分为空间域滤波技术和频率域滤波技术.空间域滤波器主要有传统方法和局域统计自适应滤波方法.前者包括诸如均值滤波器和中值滤波器［3］，它们的缺点在于对噪声和边缘信息不加区分.后者以 Lee［4］，Frost［5］，Gamma－MAP［6］等为代表，这类滤波器在对噪声进行较强平滑的同时对边缘尽量予以保留.频率域滤波器方面，小波分析由于具有灵活的多尺度特征，在SAR图像滤波方面的应用越来越突出，软阈值和硬阈值等［7］方法是此类型滤波器的典型代表.Contourlet变换［8］是对小波变换的一种新扩展，具有多分辨率、局部定位、多方向性和各向异性等性质.隐马尔可夫树（hidden markov tree，HMT）模型［9］是隐马尔可夫模型（hidden markov model，HMM）的一种树状表达，将contourlet变换和HMT结合产生了一种新的统计建模数学工具，称为Contourlet域隐马尔可夫树（contourlet－domain hidden markov tree，CHMT）模型［10］.

1 Contourlet域隐马尔可夫树模型

1.1 Contourlet变换

Contourlet变 换［8，11］可 有 效 弥 补 小 波、脊 波［12］和曲波［13］等的不足［14］，不仅可提供任意方向上的信息，而且采样冗余度小（≤4／3），其计算复杂度远低于曲波变换，在图像去噪、纹理检测和压缩编码等方面有很好的应用前景.Contourlet变换将多尺度分析和方向分析分开进行，先用拉普拉斯金字塔［15］结构对图像进行多尺度分解，获得不连续的点，接着由方向滤波器组将分布在同方向上的奇异点合成为一个contourlet系数.

1.2 高斯混合模型

自然图像经Contourlet变换后，其系数分布函数可以用一个具有2个状态的高斯混合模型来很好地表示［11］，即：Contourlet系数（W）由大量的小值系数和少量的大值系数所组成.系数分布的概率密度函数表示如下：

式中：S表示系数可能的状态，可以取l或者是s；P（S）表示状态S的概率分布函数；f（W｜S）是一个高斯函数，表示状态为S时，Contourlet系数的概率密度函数.

图1是真实影像Peppers（辣椒）经过Contourlet变换后，取其中某一子波段，然后利用期望最大化（expectation maximum，EM）算法［16］估计其统计参数，得到两个高斯分布.从图1中可以明显看出，该高斯混合分布（图中实线所示）可以由两个高斯分布（图中虚线所示）来合成，两个高斯分布的均值分别为3和5，方差分别为26和716.从参数数值上来看，两个高斯分布的均值都接近于零，其中一个的方差（716）远大于另一个的方差（26），符合前面所进行的假设.

1.3 隐马尔可夫树模型

图1 经contourlet变换后子波段的概率分布函数Fig.1 The PDF of the Peppers image’s sub－band image decomposed with contourlet

隐马尔可夫模型最重要的特点是采用Contourlet系数对应的隐状态变量而不是系数本身来体现系数之间的马尔可夫相关性.隐马尔可夫树模型是一种利用树结构建立起来的HMM模型，该模型有一种自然的父子依赖关系，一个节点的状态仅仅依赖于它的直接父节点和子节点，状态变量的依赖性通过父节点的状态到其子节点状态之间的转移概率来构建.因此，HMT抓住了多尺度分析中跨尺度关联的最重要的特征.

综合考虑高斯混合模型和Contourlet变换后父、子节点尺度间的关系，可以建立由参数集θHMT＝…，M｝组成的CHMT模型，其中，ps1（m）代表根结点S1的状态概率分布函数表示所有父、子节点之间的状态转移概率，即在父节点状态为r时，子节点状态为m的概率；μi，m和分别代表每个Contourlet系数i对应的在不同状态m下的均值和方差；P代表系数的个数；M代表状态数.这样每个Contourlet系数都可以用4个参数进行表示.通过EM算法求得这些参数，就建立了用参数集θHMT表示的CHMT模型.

1.4 基于粗分类的CHMT模型参数训练

利用EM算法求解模型参数的思路虽然较为简单，但在实际计算过程中，只有先解决参数量过大的问题才能求解模型参数.一般使用绑定的方法来解决这个问题，其实质就是将节点进行归类［17］.常用的方法包括树间绑定和树内绑定［9］.树间绑定（CHMT－M）又称多树算法，它认为把子波段用树结构表示后，在同一位置上相应节点具有相同的统计特性；树内绑定（CHMT－S）又称单树算法，它假设在同一树内同一尺度上的所有节点具有相同的统计特性.图2是树间绑定与树内绑定的示意图.

图2 系数间不同的绑定方式Fig.2 The different tying method between the contourlet coefficients

绑定的目的是为了给模型训练提供足够的训练数据，绑定时假设被绑定的点具有相同的属性或统计特性.如果在绑定前能够预先获知各个点的属性，将属于同一类的点进行绑定，则可以极大地提高模型训练的精度.为了获得这种分类信息，本文采用K均值方法对Contourlet变换中粗尺度的方向子波段系数进行粗分类（CC－CHMT），根据分类的结果，决定哪些树需要进行绑定，然后利用多树算法，对整个图像进行训练.由于事先对粗尺度的系数进行了分类，因此在实现绑定后，被绑定的系数之间有极大的“相似性”，利用这种相似性，使得模型训练收敛的速度得到提高.本文选用两状态零均值的Gauss混合模型，对尺寸为256像素×256像素的Barbara（芭芭拉）图像进行3层分解，每层的方向数分别为4，8，16，选用‘9－7’作为拉普拉斯滤波器组，‘PKVA’作为方向滤波器组［8］，收敛阈值为0.01.各种绑定方法的训练迭代次数和CPU耗时列于表1.可以看到，CCCHMT算法只需迭代66次，耗时10s即可达到收敛，其速度远高于其他两种方法.

表1 各种绑定方法的收敛速度Tab.1 The convergence speed of different tying method

2 基于SAR滤波处理统一框架和CHMT模型的滤波算法

本文建立的SAR滤波处理统一框架可表示为

式中：Y代表含噪SAR图像；X代表滤波处理后的图像；lg代表对数变换，目的是将乘性噪声转变为加性噪声；C［·］和C－1［·］分别代表Contourlet变换和逆变换；L表示循环平移的次数.

Zi和Z－i分别代表步长为i的循环平移变换和循环平移逆变换，从而消除Contourlet变换可能造成的图像可视效果畸变.循环平移可以有效地解决不满足平移不变性变换的问题，提高结果的视觉效果［18］.循环平移方法的实现过程如下：假设有一幅图像，让图像在水平方向以一定的步长S重复移动K1次，同样在垂直方向重复移动K2次，那么，可以得到K1·K2个循环平移的图像，这些图像被称为原始图像的循环平移图像；

R（·）代表均值校正处理，从而使校正的结果满足零均值高斯白噪声的假设前提.这是由于对数变换后噪声的均值不为零［19］，而通常滤波算法均建立在零均值高斯白噪声的假设前提下，特别是对于CHMT模型来说，经常假设Contourlet系数是由两个零均值的高斯混合分布构成的，因此，非常有必要对变换后的图像进行均值校正.均值校正应在进行指数变换之前，只要从Contourlet逆变换后的影像中减去经过对数变换的噪声影像的均值即可.

w（·）代表滤波处理算法.w（·）既可以是下文介绍的CHMT方法，也可以是其他不满足平移不变性变换的算法（如小波软阈值法、硬阈值法等）.这也是为该框架冠以“统一”一词的原因.该框架的流程如图3所示.

图3 SAR滤波处理统一框架示意图Fig.3 The sketch map of the universal SAR despeckle method

在建立SAR滤波处理统一框架之后，利用CHMT模型来设计滤波处理算法w（）.由于Contourlet变换是线性变换，因此，含加性噪声的影像经Contourlet变换后，仍然可以用加性噪声模型来表示，即

式中：y＝C（Y），代表含噪声的Contourlet系数；x＝C（X），表示原始信号的Contourlet系数；n＝C（N），表示噪声的Contourlet系数.去噪问题可以表述为：在已知含噪影像y的情况下对原始影像x进行估计.首先，对y建立CHMT模型，并运用基于粗分类绑定的EM算法对其参数进行估计，获取其模型参数（假设 Contourlet系数服从由m个零均值的混合高斯分布，m＝2），然后，通过将噪声的方差从θy中去除，可以获得去噪的Contourlet系数的参数θx，即

式中：［x］＋表示如果x≥0，则［x］＋＝x，否则［x］＋＝0；（j，k，i）表示在尺度为j、方向为k的子波段中的第i个Contourlet系数；m为Contourlet系数所属的状态.其中，噪声系数的方差（）2可以通过蒙特卡罗方法获得［10］.需要指出的是，通过不同的绑定方法，获得的模型参数）2必然不同.获得去噪影像的CHMT参数θx之后，就可以利用该参数对去噪影像的Contourlet系数进行估计了.假设噪声是均值为零的高斯白噪声，并且与Contourlet系数无关，因此，可以得到

其中，Sj，k，i表示在尺度为j、方向为k的子波段中的第i个Contourlet系数的状态.然后，利用在EM算法时获得的条件概率p（Sj，k，i＝m｜yj，k，i，θx），最终可以获得去噪影像x的估计为［9］

3 实验结果与分析

本文分别选用模拟图像和真实SAR图像，对多种滤波器的效果进行评估.模拟影像是向一幅快鸟（quickbird）影像中加入方差为0.05（该值已经过量纲一化，即在加入噪声时，假设图像的值在0～1之间）的服从伽马（Gamma）分布的斑点（speckle）噪声生成.SAR影像选用Radarsat影像（C波段，HH极化，空间分辨率25m×25m）.图像大小均为256像素×256像素.采用的滤波方法包括基于斑点噪声模型的Lee方法、小波软阈值法（WST）和3种采用不同绑定算法求解的CHMT方法，分别为多树绑定法（CHMT－M）、单树绑定法（CHMT－S）和基于K均值粗分类的绑定法（CC－CHMT）.CHMT模型训练中均使用零均值、两状态的混合高斯分布，拉普拉斯分解选用‘9－7’滤波器，方向滤波器组均选用PKVA，Contourlet分解层数为3层，每层的方向数分别为4，8，16；小波域滤波器中，为了避免使用不同的小波基对结果产生的影响，小波基选用与Contourlet变换相同的‘9－7’，分解层数也为3层.

图4是对模拟影像使用不同滤波器的处理结果.由于CHMT－M和CHMT－S方法的目视效果较为一致，为节省篇幅此处未列出CHMT－S的滤波结果.由图可见，快鸟影像加入噪声后退化得比较严重.从滤波结果来看，Lee方法虽然细节信息保持得比较好，但抑噪的能力明显低于其他方法；WST方法产生了较严重的“振铃”效应，丧失了很多细节信息；CHMT－M方法与CC－CHMT的结果比较接近.从总体上来说，基于Contourlet域的HMM方法（CHMT－M和CC－CHMT）得到的结果可视效果好于Lee方法和WST方法，其结果在抑噪和细节保留上均取得了较好的效果，其可视效果也更接近于快鸟影像.

为了进一步研究各种滤波方法的滤波能力，在图像上任取一行，通过其剖面对滤波后影像与原始影像进行比较.剖面的位置如图4a中的直线所示.这条线经过了很多边缘，便于比较不同的滤波器在处理边缘时的能力.

图5显示了原始影像、Lee方法、WST方法和CC－CHMT方法的剖面曲线（因为CHMT－M方法与CC－CHMT方法的曲线非常接近，为了让结果更加清晰，没有在图中标出CHMT－M）.可以看到，在Lee方法的曲线上产生了很多小的“尖峰”，与原始图像的曲线不能很好地吻合，表现在图像上就是噪声压制的效果不强；WST方法是与原始图像重叠的最好的，但从图4d中可以发现，该方法得到的图像有较强的“振铃”效应；CC－CHMT方法在大多数地方能够与原始图像的曲线较好地重叠.

图5 不同滤波方法的剖面图Fig.5 The profile of the different filters

图6 a为Radarsat影像滤波结果.由图可见，Lee方法在滤除噪声的同时，较大地模糊了影像（比如左上角以及中央的高反射区域经过滤波后，几乎合并成一片），纹理细节信息丢失较严重，这可能与Lee滤波器判定这片区域为匀质区域而进行了较强平滑有关；WST方法的视觉效果较差，有明显的“振铃”现象，这是由于小波变换采用了下采样技术而丧失了平移不变性，因此在这幅边缘信息丰富的图像上产生了较严重的伪吉布斯效应；由于CHMT模型可以有效地“捕获”影像的方向性信息，再通过利用SAR统一滤波框架弥补了Contourlet变换下采样引起的非平移不变性对视觉效果的影响，3种CHMT方法在有效滤除噪声的同时，都具备较强的细节保留能力，滤波结果比Lee方法和WST方法自然、清晰.

使用等效视数和信噪比这两个参数来定量评价滤波器的性能.等效视数可以评价一个滤波器的斑点抑制能力，该值越大，表明滤波器对噪声的压制效果越强，在匀质区域，等效视数定义为该区域均值的平方除以该区域的方差.信噪比表示含噪声图像或者滤波后图像与原始图像相比，质量变化的情况.信噪比越高，表示滤波的效果越好.

表2列出了对噪声影像利用不同滤波器进行滤波后的参数指标.从模拟影像的统计结果可以看到，Lee方法不论是等效视数还是信噪比都远低于其他几种方法；WST方法有较高的等效视数和信噪比，但是通过观察发现，这两种方法都在很大程度上模糊了影像，也就是说，由于该方法具有很强的“平滑”能力，才导致了较高的等效视数和信噪比值；从指标值来看，几种CHMT方法都能取得较好效果，且没有很大的差别，CC－CHMT方法在等效视数上略高于其他两种CHMT方法，其信噪比与另两种CHMT方法基本相当.对Radarsat影像的统计结果也有相似的结论，不再赘述.

表2 滤波器性能比较Tab.2 The statistical properties of the different filters

4 结论

在SAR图像乘性斑点噪声模型的基础上，通过分别应用对数变换将乘性噪声转化为加性噪声，应用循环平移方法消除图像可视效果畸变，使用均值校正解决对数变换后的均值偏移，建立了统一的SAR图像滤波处理框架；在此框架之下，建立了基于两状态零均值的Contourlet域隐马尔可夫树模型滤波方法；针对模型参数求解速度慢的问题，提出了基于粗分类的系数绑定训练算法，实验证明该方法的收敛速度高于多树算法和单树算法；最后应用几种常用的SAR滤波算法，对模拟影像和真实SAR影像进行了滤波处理，通过比较各种方法的可视效果和统计指标，表明在SAR滤波处理统一框架下，应用基于粗分类绑定训练算法的Contourlet域隐马尔可夫树模型的滤波方法，可以极大地提高CHMT模型训练算法的速度，其滤波结果在可视效果和统计指标方面均优于常用的几种SAR滤波算法.

［1］徐新，廖明生，卜方玲.一种基于相对标准差的SAR图像Speckle滤波方法［J］.遥感学报，2000，4（3）：214.XU Xin，LIAO Mingsheng，PU Fangling.An adaptive method for speckle filtering of SAR image based on local relative standard deviation［J］.Journal of Remote Sensing.2000，4（3）：214.

［2］Porcello L J，Massey N G，Innes R B，et al.Speckle reduction in synthetic radars［J］.Journal of the Optical Society of America，1976，66：1305.

［3］Rees W G，Satchell M J.The effect of median filtering on synthetic aperture radar images［J］.International Journal of Remote Sensing，1997，18（13）：2887.

［4］Lee J S.Speckle suppression and analysis for synthetic aperture radar images［J］.Optical Engineering，1986，25（5）：636.

［5］Frost V S.An adaptive filter for smoothing noisy radar images［J］.IEEE Proceedings，1981，69（1）：133.

［6］Lopes A，Nezry E，Touzi R，et al.Maximum a posteriori speckle filtering and first order texture models in SAR images［C］／／Proceedings of IGARSS＇90.［S.l.］：Institute of Electrical and Electronics Engineers，1990：2409－2412.

［7］Donoho D L.Denoising by soft－thresholding［J］.IEEE Transactions on Information Theory，1995，41（3）：613.

［8］Po D D Y，Do M N.Directional multiscale modeling of images using the contourlet transform［J］.IEEE Transactions on Image Processing，2005，15（16）：2691.

［9］Crouse M S，Nowak R D，Baraniuk R G.Wavelet－based statistical signal processing using hidden markov models［J］.IEEE Transactions on Signal Processing，1998，46（4）：886.

［10］Li J，Gray R M， Olshen R A.Multiresolution image classification by hierarchical modeling with two－dimensional hidden Markov models［J］.IEEE Transactions on Information Theory，2000，46（5）：1826.

［11］Po D D Y，Do M N.Directional multiscale modeling of images using the contourlet transform［J］.IEEE Transactions on Image Processing，2006，15（6）：1610.

［12］Candes E J.Ridgelet：theory and applications［D］.Stanford：Stanford University，1998.

［13］Choi M，Rae Y K，Nam M R.Fusion of multispectral and panchromatic satellite images using the curvelet transform［J］.IEEE Geoscience and Remote Sensing Letters，2005，2（2）：136.

［14］易文娟，郁梅，蒋刚毅.Contourlet：一种有效的方向多尺度变换分析方法［J］.计算机应用研究，2006（9）：18.YI Wenjuan，YU Mei，JIANG Yigang.Contourlet：efficient directional and multiresolution analytic tool［J］.Application Research of Computer，2006（9）：18.

［15］Burt P J，Adelson E H.The laplacian pyramid as a compact image［J］.IEEE Transactions on Communications，1983，31（4）：532.

［16］Tomasi C.Estimating gaussian mixture densities with EM－a Tutorial［R／OL］.Durham：Duke University Department of Computer Science，2011［2011－02－10］.http：／／www.cs.duke.edu／courses／spring04／cps196.1／handouts／EM／tomasiEM.pdf.

［17］Rabiner L R.A tutorial on hidden Markov models and selected applications in speech recognition［J］.Proceedings of the IEEE，1989，77（2）：257.

［18］Guoliang F，Xianggen X.Image denoising using a local contextual hidden Markov model in the wavelet domain［J］.IEEE Signal Processing Letters，2001，8（5）：125.

［19］Xie H，Pierce L E，Ulaby F T.SAR speckle reduction using wavelet denoising and markov random field modeling［J］.IEEE Transactions on Geoscience and Remote Sensing，2002，40（10）：2196.