基于2D-PWVD的SAR图像相干斑噪声抑制

刘 聪，李言俊，刘保元

(1.中国人民解放军驻西飞公司军事代表室， 西安710089;2.西北工业大学航天学院， 西安710072)

0 引言

由于合成孔径雷达(SAR)图像的成像特点，相干斑噪声是SAR图像的固有特性，为实现对SAR图像的准确解译，研究SAR图像的噪声抑制非常有意义。目前，已提出的SAR图像相干斑噪声抑制方法主要分为时域滤波方法和变换域滤波方法两类。时域滤波方法主要包括:Lee滤波算法、Kuan滤波算法、Frost滤波及其增强算法、Gamma-MAP滤波算法等。变换域滤波算法主要是基于小波变换的滤波算法。

基于时域滤波的所有方法都是利用SAR图像的统计特性，利用一个滑动窗口与原图像做卷积，达到对图像平滑的目的。这类方法实际上是以牺牲图像细节为代价，达到噪声抑制的目的，这就造成了SAR图像斑块抑制效果与图像的边缘细节保留之间的矛盾，通常增大滑动窗口的尺寸，能提高斑点噪声抑制的效果，但同时也造成更多图像结构信息的丢失，为保证二者的兼顾，窗口尺寸的合理选择非常困难。

基于小波变换的相干斑噪声抑制方法是将斑点噪声转换为高频各层次的小波系数，根据高频子域系数特征设计不同的处理策略，达到相干斑噪声抑制的目的。相比于时域滤波方法，基于小波变换的滤波方法能比较好地实现区域平滑和纹理保留的折中。围绕高频子域系数特征的处理策略，出现了多种基于小波变换的噪声抑制方法。Donodo、Johnstone［1］等提出了称为“小波收缩”的信号去噪方法，在此基础上，Devore与 Lucier［2］、B.Sanker［3］等提出了大量优化方法［4］。Romberg［5］基于多种隐马尔科夫模型(HMM)，利用小波系数间的相关依存关系，对小波系数进行了重新估计，使其逼近原始图像的小波系数，达到噪声抑制的目的。

利用伪魏格纳分布对图像进行分解能得到图像的低频和高频能量谱图，一般认为图像的边缘噪声能量主要集中在高频区域，受小波软阈值SAR图像噪声抑制方法的启发，结合魏格纳分布的特性，本文阐述了基于二维伪魏格纳分布的SAR图像噪声抑制方法。

1 二维魏格纳分布理论

Wigner于1932年首先提出了Wigner分布的概念，并应用于量子力学领域。1948年Ville将Wigner分布引入信号分析领域，应用于信号的检测和信号细节分类。因此Wigner分布又称为Wigner-Ville分布(WVD)，Wigner-Ville分布属于双线性时频分析Cohen＇s类的特殊形式，具体表达式为

式中:*代表复数的共轭;t为时间分量;ω为频率分量;W(t，ω)为信号f(t)的魏格纳分布。

Wigner-Ville分布具有一系列好的性质，如奇偶性、虚实性、时间边缘性和频率边缘性等。Jacobson和Wechsler在20世纪80年代将Wigner-Ville分布拓展到二维信号和三维信号的处理中，并阐述了多维Wigner-Ville分布与一维Wigner-Ville分布具有相同的特性［6-8］。二维 Wigner-Ville分布(2D-WVD)定义如下

式中:Wf(x，y，u，v)为信号 f(x，y)的二维维格纳分布;(x，y)为时间分量;(u，v)为频率分量。

离散形式定义为

信号的WVD变换存在交叉项的问题，交叉项是信号的N元成分经WVD变换后形成的N(N-1)/2项幅度，为各信号成分两倍的干扰项，干扰项的存在严重影响对信号的准确认识，伪魏格纳分布(PWVD)通过在时域加窗能有效抑制多元成分信号的交叉项，达到正确认识信号的目的，二维PWVD(2D-PWVD)定义为

由式(5)可以看出，选择大小为(2N1-1)×(2N2-1)大小的时域窗口，针对二维信号的每个像素点进行2D-PWVD，可以得到整个二维信号的PWVD变换图。

2 图像的二维伪魏格纳分解

图像为二维数据，根据表达式(3)，对M×N图像进行2D-WVD变换，需要建立一个M×N×M×N大小的四维数组，若再对M×N×M×N四维数组进行相关处理，计算量相当大，不利于现实应用，且图像多元成分的交叉项影响也相当严重。根据局部频谱思想，本文采用表达式(5)对原图像进行2D-PWVD变换。在时域上，选取大小为N1×N2(要求N1=N2，且为奇数)的窗口函数在原图像进行滑动，得到原图像每个像素点大小为N1×N2的频谱图像。定义I(i，j)为数字图像空间坐标(i，j)处的像素值，Xij为 I(i，j)的根据式(5)计算得到的2D-PWVD，其中

N1×N2为二维滑动窗口大小。代表像素I(i，j)在(k，j)频率点的PWVD 值，其中k=1，2，…，N1，l=1，2，…，N2。提取图像各空间位置的‖‖，按照图像对应空间位置组成“能量谱图”，一幅图像能得到N1×N2幅与原图像同样大小的“能量谱图”。分解流程如图1所示。分解图像的排列方式如图2所示。

图1 图像数据的2D-PWVD分解流程

图2 分解图像的排列方式

选取时域窗口函数h(k，l)=1(矩形窗)，大小为5×5，通过图1的变换重组流程，将lenna原始图像分解为25个不同的频段图(图3)。从分解重组图可以看出，图3(11)相对于lenna原始图像，图像边缘模糊，可以理解为是原始图像经过2D-PWVD变换重组的低通滤波图像，分解图像(12)～(55)是原始图像边缘的信息体现，可以理解为是原始图像经过2D-PWVD变换重组的高通滤波图像。仔细分析lenna图像的分解图像，可以看出:(12)和(15)、(13)和(14)、(21)和(51)、(31)和(41)、(22)和(55)、(23)和(54)、(24)和(53)、(25)和(52)、(32)和(45)、(33)和(44)、(34)和(43)、(35)和(42)分别是完全相同的。因此，lenna图像经过2D-PWVD变换重组，实际上形成了13个分解图像。推广到所有时域窗口的情况下，分解重组后形成的图像数为(N1×N2+1)/2。

图3 lenna图像的2D-PWVD 5×5分解图

在图3的高频分解图像中，(12)和(13)对原图像的水平向变化剧烈的边缘敏感，(21)和(31)对原图像的垂直向变化剧烈的边缘敏感，(22)和(33)对原图像反对角线和对角线方向变化剧烈边缘敏感，对角线方向敏感更强，(23)和(32)也对原图像反对角线和对角线方向变化剧烈边缘敏感，但反对角线方向敏感更强。通过增加时域窗口的大小为 7×7，9×9，11×11，13×13(由于分解图像众多，在此只表示5×5的情况，如图3所示)，并分析分解重组图像，得出结论(按照图2分解图像排列方式，且不包括11频段图像):(1)第一行表示原图像的水平向高频图像，第一列表示原图像的垂直向高频图像，剩余分解高频图像表示方向如图4所示;(2)高频分解方向数为2×N1-2，其中N1为时域窗口大小;(3)随着时域窗口的增加，多元成分的干扰项越严重。

图4 图像2D-PWVD分解图像方向表示图

3 基于2D伪魏格纳分布的SAR图像噪声抑制

Fukuda等［9］根据小波分解高频系数具有不同的方向性特点，构造了3个如图5所示的滤波窗口，来判断高频细节中图像边缘和噪声，达到抑制噪声的目的。

图5 高频子带定向滤波器

在本文第2节阐述二维伪魏格纳威尔分布分解后具有2*N-1个方向(N表示时域窗口的大小)，在Fukuda提出算法的启示下，考虑高频图谱的方向性，使SAR图像噪声抑制后图像边缘特性得到保持。在本算法中，同样采用图5的滤波窗口，在分解的图像的垂直向高频图像中，利用图5a)滤波窗口，水平高频图像利用图5b)滤波窗口，其他高频段图像利用图5c)滤波窗口。具体算法如下:

(1)选择时域窗口大小为N×N。

(2)将原图像进行二维魏格纳威尔变换，形成N×N幅分解图像

式中:X(i，j)、Y(i，j)分别为计算后各分解图像的各点像素值和计算前各分解图像各像素值。

(3)针对各高频分解图像，计算各分解图像的均值μ，根据如下原则

计算阈值T，将各分解图像分成噪声区域和可能细节区域。区分原则为:X(i，j)＜T为噪声区域，X(i，j)≥T认为为可能细节区域。式中，X(i，j)为各分解高频图像的各点像素值。

(4)在噪声区域

Xchu(i，j)=k1× X(i，j)

式中:Xchu(i，j)为处理后各点像素值;k1为噪声抑制系数。

(5)在可能细节区域，当除中心点之外坐标中有一个点的值满足X(i，j)≥T(图5中深色的点)，认为当前像素值为目标边缘细节，则 Xchu(i，j)=X(i，j)，否则认为是噪声，Xchu(i，j)=k1×X(i，j)。

(6)将处理后的各分解图像相加，形成SAR噪声抑制后图像。

4 SAR图像相干斑抑制性能的评估

SAR图像经过相干斑抑制算法处理后图像质量是SAR图像处理关心的重点，因此必须用一定的标准来衡量SAR图像的质量。SAR图像的应用十分广泛，不同的用途对SAR图像质量的要求也不完全一样。一般来说，都要求图像具有较好的区分临近目标的能力和检测目标的能力，要求图像具有较丰富的细节信息［10］。

目前，已提出的SAR图像斑点噪声抑制质量定量评估指标主要包括等效视数和边缘保持指数，下面分别介绍各个指标的定义。

1)等效视数［11］

等效视数是衡量一幅图像相干斑噪声相对强度的一个有效指标，等效视数越大，表明图像上的相干斑越弱，可解译性越好。定义为

2)边缘保持指数［12］

边缘保持指数定义为

式中:N为图像的像素点数;ps，psn分别为滤波后图像的像素值和该像素点水平或垂直方向的领域值;po，pon分别为原始图像的像素值和该像素点水平或垂直方向的邻域值。因此，边缘保持指数可分为水平边缘保持指数和垂直方向边缘保持指数。

5 噪声抑制试验与结果分析

为验证本文提出SAR图像噪声抑制算法的有效性，利用圣地亚哥国家实验室提供的MiniSAR真实合成孔径雷达数据进行试验。

MiniSAR数据(目前公布19幅图像)，在2005年由圣地亚哥国家实验室获取，采用Ka和Ku波段聚束式成像方式，最高分辨率为0.1 m×0.1 m，图像未进行多视处理，图像为复杂场景。

针对均值滤波算法［11］、LEE 滤波算法［13］、软阈值小波滤波算法［1］和本文提出的噪声抑制算法进行对比试验，性能判别指标选用等效视数和边缘保持系数。在试验中，均值滤波算法采用3×3的矩形窗口，LEE滤波算法窗口大小7×7(MiniSAR数据)，软阈值小波滤波算法中的α和β选择为0.5，本文算法的k1=0.5，h(k，l)=1，大小为3×3的矩形窗口。由于篇幅有限，列出第9幅图像(如图6所示)的实验结果(见图7和表1)。

图6 原始图像

图7为图6图像采用不同方法噪声抑制后的效果图。表1是采用不同方法的噪声抑制性能结果。

图7 MiniSAR数据的第9幅图像噪声抑制后的效果图

表1 MiniSAR数据噪声抑制方法性能参数表

从图7可以看出，采用本文算法和软阈值小波滤波算法抑制图像从视觉效果上基本一致。

从表1可以看出，采用本文算法得到的等效视数相对于原图像和软阈值小波滤波算法有明显提高，但略低于采用均值滤波和LEE滤波得到的对应参数;本文算法计算的边缘保持系数明显高于均值滤波和LEE滤波算法，与小波滤波算法相当。

为验证本文算法的鲁棒性，针对等效视数和边缘保持系数指标，表2统计了本文算法相对于均值滤波算法、LEE滤波算法和软阈值小波滤波算法对19幅MiniSAR数据实验结果。

从表2可以看出，针对MiniSAR数据，本文算法在等效视数指标上不如均值滤波算法和LEE滤波算法优越，在边缘保持系数指标上明显高于其他3种算法，相对于软阈值小波滤波算法，无论从等效视数指标还是边缘保持系数指标，都具有一定的优越性。

从各种算法的原理上分析，本文算法将原图像按照不同方向进行分解，并对目标细节和噪声中可能目标细节进行分别处理，在保留图像中尽可能多细节的同时，对噪声进行了有效抑制。可以认为，本文算法在保持较好的边缘保持系数的条件下，提高了SAR图像的等效视数，达到了抑制SAR图像噪声的目的。

从算法效率上来分析，由于本文算法需要将原图像分解为N幅(与时域窗口得大小相关)与原图像同尺寸的子图像，并需要考虑目标区域、强噪声中可能的目标细节，因此算法效率上比均值滤波算法、LEE滤波算法和软阈值小波滤波算法要低。

6 结束语

本文根据2D-PWVD变换的频谱特性，提出了一种新的SAR图像噪声抑制算法，利用MiniSAR图像数据进行验证试验，并与均值滤波、LEE滤波、软阈值小波滤波进行对比分析，说明本文算法在抑制SAR图像噪声的同时，能更多保留原图像的细节信息，证明本文算法的有效性和正确性。当然，本算法的执行效率还有待提高，以利于实际工程的具体应用。

［1］ Donoho D L，Johnstone I M.Nonlinear wavelet methods for recorvery of signal densities and spectra from indirect and noisy dat［C］//Proceedings of Symosia in Applied Mathematics.Providence Rhode Island:American Mathematical Socity，1993:173-205.

［2］ Devore R A，Lucier B J.Fast wavelet techniques for nearoptimal image processing［C］//Military Communications Conference on Fusing Command，Control and Intelligence.San Diego，CA:IEEE Press，1992:1129-1135.

［3］ Sanker B，Guler E G，Kahya Y P.Multiresolution biological transient extraction applied to respiratory crackles［J］.Computers in Biology and Medicine，1996，26(1):25-39.

［4］ 朱家兵，陶 亮，江有名，等.一种新的高分辨率SAR图像相干斑噪声抑制算法［J］.现代雷达，2005，27(11):54-57，74.Zhu Jiabing，Tao Liang，Jiang Youming，et al.A new speckle reduction algorithm for high resolution SAR images［J］.Modern Radar，2005 ，27(11):54-57，74.

［5］ Romberg J K，Baraniuk R G.Bayesian tree-structured image mode-ling using wavelet-domain hidden markov model［J］.IEEE Transactions on Image Processing，2001，10(7):1056-1068.

［6］ Jacobson L，Wechsler H.A paradigm for invariant object recognition of brightness，optical flow and binocular disparity images［J］.Pattern Recognition Letters，1982(1):61-68.

［7］ Jacobson L，Wechsler H.A theory for invariant object recognition in the frontoparallel plane［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1984，6(3):325-331.

［8］ Jacobson L，Wechsler H.Derivation of optical flow using a spatiotemporal-frequency approach［J］.Computer Vision，Graphics，and Image Processing，1987，38(1):29-65.

［9］ Fukuda S，Hirosawa H.Supression of speckle in synthetic aperture radar images using wavelets［J］.International Journal of Remote Sensing，1998，19(3):507-510.

［10］ 匡纲要，高 贵，蒋咏梅，等.合成孔径雷达目标检测理论、算法及应用［M］.长沙:国防科技大学出版社，2007.Kuang Gangyao，Gao Gui，Jiang Yongmei，et al.Synthetic aperture radar target detection theory algorithm applications［J］.Changsha:National University of Defense Technology press，2007.

［11］ Chris Olive Shaun Quegan.Understanding synthetic apearture radar images［M］.Raleigh，NC:SciTech Publishing，Inc，2004.

［12］ Sheng Yongwei，Xia Zongguo.A comprehensive evaluation of filters for radar speckle suppression［C］//International Geoscience and Remote Sensing Symposium.Lincoln，NE:IEEE Press，1996:1559-1561.

［13］ Lee J S.Digital image enhancement and noise filtering by using local statistics［J］.IEEE Transactions on Pattern and Machine Intelligence，1980，2(2):165-168.