结合空时上下文信息的视频SAR图像相干斑滤波

刘燊文, 崔兴超, 陈思伟

(国防科技大学电子科学学院电子信息系统复杂电磁环境效应国家重点实验室, 湖南 长沙 410073)

0 引 言

合成孔径雷达(synthetic aperture radar, SAR)能够全天时全天候工作,已成为对地观测领域的主流传感器之一。传统SAR成像系统往往只能得到场景的单幅图像。近年来,视频SAR技术得到发展。视频SAR系统能够获得同一场景的SAR图像序列,实现对场景的持续观测,并将传统SAR的方位-距离二维空域散射信息拓展至方位-距离-时间三维空时散射信息[1-2]。

由于SAR系统是一个相干成像系统,SAR图像上存在相干斑现象,会影响目标检测、分类和识别等后续应用[3-4]。视频SAR图像和传统SAR图像一样也受相干斑影响,因此相干斑滤波通常也是视频SAR必要的预处理过程。相干斑滤波主要包含两个步骤:一是相似候选样本像素的选取,二是无偏估计器的构建[5]。其中,常用的无偏估计器包括最小均方误差估计器、最大后验概率估计器等,研究较为充分。而相似候选样本像素的选取是决定相干斑滤波性能的关键,是该领域的研究重点。

当前通过对SAR图像的杂波散射特性和相干斑机理进行深入的分析,研究人员提出了一系列SAR图像滤波算法[6-10]。这些滤波算法主要可分为空域滤波、变换域滤波两大类。空域滤波通过计算邻域像素与目标像素相似性来选取同质像素,进而通过估计目标像素实现滤波。代表性方法有Lee滤波器、Frost滤波器、Kuan滤波器、Sigma滤波器等[11]。此外,基于区域生长的局部滤波算法(intensity driven adaptive neighborhood, IDAN)[12]、加权估计的非局部均值(non-local means)滤波方法[13-15]、基于概率性块(probabilistic patch-based, PPB)的迭代滤波算法[16]、基于偏微分方程的滤波方法[17]等则进一步改进了像素间的相似性距离计算方式,不同程度地提升了滤波性能。变换域滤波方法则是将图像域通过特定变换转换到其他域再进行滤波,代表性方法有小波变换[18]、Contourlet变换[19]和Shearlet变换[20]等,以及基于块匹配的三维协同小波去噪方法[21]。值得一提的是,近年来随着深度学习在图像处理领域的广泛应用,也出现了一些基于深度学习的相干斑滤波方法[22-26]。

视频SAR相干斑滤波目前主要有两种方式:第1种是基于单幅SAR图像空域信息的相干斑滤波方法,这类方法将图像序列分成单幅SAR图像,分别进行滤波;第2种是迁移应用针对光学视频的滤波方法,如视频块匹配三维滤波方法(video-block match and 3D filtering, V-BM3D)[27-28],该方法在一定程度上也能滤除相干斑,然而这些方法都没能充分结合SAR图像特性进行改进,在噪声标准差过大时,块匹配搜索的准确性严重降低,去噪性能会出现急剧下降[29]。

视频SAR能够获取成像场景的方位-距离-时间三维动态散射信息,显著提升了信息量[30],提供了从时间维进行相干斑滤波的新途径。基于此,本文首先构造了空时上下文散射矢量和空时上下文协方差矩阵,用于表征以视频SAR为代表的序列SAR图像全维度信息。在此基础上,提出了结合空时上下文协方差矩阵的视频SAR图像相干斑滤波方法。最后,结合视频SAR数据开展实验,对所提方法的滤波性能进行评估分析。

1 视频SAR空时上下文协方差矩阵

1.1 视频SAR数据表征

对于单幅SAR图像,其像素分布在方位向和距离向两个空间维度上,记作p(x,y),其中下标(x,y)表示像素的距离向和方位向坐标。

视频SAR图像序列如图1所示,其像素分布在距离向、方位向两个空间维度和帧序这一时间维度上,记作p(x,y,t),下标(x,y,t)表示像素的距离向坐标、方位向坐标和时间向坐标。

图1 视频SAR数据表征三维散射信息Fig.1 3D scattering information represented by VideoSAR data

1.2 空时上下文协方差矩阵

视频SAR数据能够提供三维散射信息,为利用视频SAR图像序列的空时上下文信息,在前期SAR图像上下文协方差矩阵表征方法[5]基础上,本文提出了序列SAR图像空时上下文协方差矩阵新表征形式。对任一像素p(x,y,t),可构造空时上下文协方差矩阵M,用以代替复数p表征该处的电磁散射特性。空时上下文协方差矩阵的构造过程主要分为两步:构造空时上下文矢量和计算空时上下文协方差矩阵。

1.2.1 构造空时上下文散射矢量

以3×3×3邻域为例,像素p(x,y,t)的空时上下文散射矢量v的构造原理如图2所示。

图2 空时上下文矢量示意图Fig.2 Schematic diagram of space-time context vector

构造的13个空时上下文矢量分别为

(1)

式中:上标“T”表示转置。

值得说明的是,对其他n1×n2×n3邻域也可类似构造空时上下文散射矢量。考虑到维度越大,空时上下文散射矢量的数量将变得越多,通常可取3×3×3、5×5×5的邻域进行处理。本文采用3×3×3邻域。

1.2.2 构造空时上下文协方差矩阵

对视频SAR,空时上下文协方差矩阵M(x,y,t)的构造方法为

(2)

式中:上标“H”表示共轭转置;可知M(x,y,t)是一个3×3的Hermitian矩阵,其各个元素含义如下:

(1)M(1,1)表示前一帧9个邻域像素(见图2中的绿色像素)和当前帧左上角4个邻域像素(像素p(x-1,y-1,t)、p(x-1,y,t)、p(x-1,y+1,t)、p(x,y-1,t))的强度均值;

(2)M(1,2)、M(2,1)表示前一帧9个邻域像素(见图2中的绿色像素)和当前帧左上角4个邻域像素(像素p(x-1,y-1,t)、p(x-1,y,t)、p(x-1,y+1,t)、p(x,y-1,t))与中心像素p(x,y,t)(见图2中的红色像素)相关值的均值;

(3)M(1,3)、M(3,1)表示13对关于中心像素p(x,y,t)中心对称的像素(见图2中13个空时上下文矢量两端像素)相关值的均值;

(4)M(2,2)表示中心像素p(x,y,t)的强度;

(5)M(2,3)、M(3,2)表示后一帧9个邻域像素(见图2中的黄色像素)和当前帧右下角4个邻域像素(像素p(x,y+1,t)、p(x+1,y-1,t)、p(x+1,y,t)、p(x+1,y+1,t))与中心像素p(x,y,t)(见图2中的红色像素)相关值的均值;

(6)M(3,3)表示后一帧9个邻域像素(见图2中的黄色像素)和当前帧右下角4个邻域像素(像素p(x,y+1,t)、p(x+1,y-1,t)、p(x+1,y,t)、p(x+1,y+1,t))的强度均值。

对于实际视频SAR图像数据,为使得上下文协方差矩阵M(x,y,t)满秩,通常还需计算出像素p(x+i,y+j,t)对应的上下文协方差矩阵M(x+i,y+j,t),以它们的均值作为最终的M(x,y,t),一般取-1≤i≤1,-1≤j≤1,即

(3)

1.3 空时上下文协方差矩阵的统计分布

一般可用复正态分布描述SAR数据的统计特性。复正态分布被广泛用于描述从低分辨率到高分辨率的同质区SAR数据[31-32]。同理可认为所构造的空时上下文矢量v服从均值为0、协方差矩阵为M的多元复高斯分布,记作v∈NC(0,M),分布函数为

(4)

式中:协方差矩阵M=E[vvH];M-1表示M的逆矩阵;|·|表示计算行列式。

设X=nM,则方阵X服从复Wishart分布,记作X∈WC(p,n,M)。分布函数为

(5)

式中:Tr(·)表示求矩阵的迹,则

(6)

式中:p为所构造的空时上下文矢量的维度,有p=3;Γ(·)为伽马函数。

2 视频SAR图像相干斑滤波方法

2.1 空时上下文协方差矩阵相似性检验

根据Conradsen等[33]的工作,可用似然比Q及其自然对数形式lnQ来检验两个相互独立的、服从复Wishart分布的p×p维的正定矩阵X、Y的相似性。

对于两个相互独立的服从复Wishart分布的p×p维的矩阵X、Y,记作

X∈WC(p,n,Mx),Y∈WC(p,m,My)

(7)

在假设Mx=My下的似然比Q为

(8)

如果n=m,并对似然比Q取自然对数形式,则有

lnQ=n(2pln 2+ln|X|+ln|Y|-2ln|X+Y|)

(9)

特别地,如果X=Y,则lnQ=0,而其他情况则有lnQ<0。

2.2 像素相似样本选取

空时上下文协方差矩阵服从复Wishart分布,利用对数似然比lnQ可检验其两两之间的相似性,进而选取相似样本像素用于相干斑滤波。

为检验上述思路的有效性,本文采用2组公开视频SAR数据[34]进行验证。视频SAR图像及其对应光学图如图3所示,视频1成像区域为基地大门,场景中主要有运动车辆、道路、停车场等。视频2成像区域为太阳能塔,观测时间为晚上9点左右,场景中存在一座高塔和大量反射镜。

图3 视频SAR图像及对应光学图Fig.3 VideoSAR screenshots and optical images of imaging scene

采用数据第247帧开展相似样本选取实验,光学图像和视频SAR图像如图4所示,选取的中心像素点在图4中标记为红色。

图4 光学图像和视频SAR图像Fig.4 Optical image and VideoSAR image

目标像素的相似性检验如图5所示。图5(a)展示的是大门区域的3个感兴趣点,代表强散射区域。像素1为大门像素点,像素2为大门与道路交界处像素点,像素3为道路像素点。lnQ子图越接近白色,表示lnQ数值越接近0,即与中心像素点越相似;掩膜中的白色点和视频SAR图像中的蓝色点表示与中心像素相似的样本。对于大门像素,其邻域内的大门的其他像素基本被选为相似样本,而没有交界处或道路上的像素点被选为相似样本;对边缘像素,其邻域内的同为交界处的像素点基本被选为相似样本,相似样本整体沿着大门两侧分布;对道路像素,其邻域内的同为道路上的像素点大部分被选为相似样本,而没有边界处或大门处的像素被误选。

图5 利用空时上下文协方差矩阵相似性检验区分异质像素Fig.5 Space-time context covariance matrix similarity test to distinguish the heterogeneous pixels

图5(b)展示的是停车场区域的3个感兴趣点,代表弱散射区。红点表示感兴趣点,蓝点表示相似样本点。像素点1为空地像素点,像素点2为停车场与空地交界处像素点,像素点3为停车场处像素点。

对空地像素和停车场像素,其邻域内大部分像素被选为相似样本;对边缘像素,其邻域内同为交界处的像素点基本被选为相似样本,仅少量其他像素被误选,相似样本整体沿着交界纹理呈曲线状分布。

综合来看,无论在强散射区还是弱散射区,空时上下文协方差矩阵相似性度量指标都能有效表征像素间的相似性,可用于相似样本像素选取。

2.3 基于空时上下文协方差矩阵的视频SAR图像相干斑滤波方法

基于空时上下文协方差矩阵和矩阵相似性度量方法,本节提出一种视频SAR图像相干斑滤波方法。考虑到空时上下文协方差矩阵表征了视频SAR三维数据块信息,因此该相干斑滤波方法也被称为3D-CCM(3-dimensional context covariance matrix)方法。该方法的主要步骤如下:

算法 1 3D-CCM滤波算法输入 视频SAR图像步骤 1 构造空时上下文协方差矩阵M(x,y,t)。步骤 2 根据相似性指标ln Q(x,y,t)-(x',y',t')选择待滤波像素p(x,y,t)的相似样本。根据公式计算滤波像素p(x,y,t)和候选样本像素p(x',y',t')的相似性指标ln Q(x,y,t)-(x',y',t')。若像素p(x,y,t)与候选样本像素p(x',y',t')的相似性指标ln Q(x,y,t)-(x',y',t')大于等于相似性判定阈值thsimitest,则将像素p(x',y',t')归为像素p(x,y,t)的相似样本,得到滤波像素p(x,y,t)在滤波窗内的所有相似样本集合U。步骤 3 构建均值滤波器得到滤波像素p(x,y,t)的估计值p^(x,y,t),记作p^(x,y,t)=1N∑Ni=1|p(x',y',t')|,其中N是相似样本集合U中像素点的个数,|p(x',y',t')|表示像素p(x',y',t')的模值。输出 相干斑滤波后的视频SAR图像

在算法实现过程中,为保持滤波后的视频SAR数据大小不变,采用了镜像拓展方法对边缘进行镜像拓展;为保护目标的强散射细节,本文设置了强散射点保护机制[11]。

3 对比实验

为验证所提滤波方法的有效性,本节以官网公开的视频SAR数据开展相干斑滤波对比实验。对比方法选用均值滤波、Improved Sigma[11]、SAR-BM3D[21]、V-BM3D[28]4种方法。参照相应参考文献,各对比方法主要参数设置为:均值滤波采用5×5的滤波窗;Improved Sigma采用9×9的滤波窗,Sigma取0.9;V-BM3D方法Sigma取25。所提方法(3D-CCM)主要参数有2个:滤波窗大小和相似性阈值。通常滤波窗越大,候选样本越多,滤波所需时间也越长;相似性阈值越小,选中的相似样本越多,则误选的相似样本会增多,可能出现过度平滑的现象。综合考虑运行时间和滤波效果,将所提方法滤波窗大小设置为25×25×3(距离×方位×时间),将相似性阈值设置为-1。

采用相干斑滤波常用评价指标:等效视数(equivalent number of look, ENL)、目标杂波比(target to clutter ratio, TCR)和边缘保持品质因数(figure of merit, FOM)来定量评价滤波效果[3,5]。

ENL的计算公式为

(10)

式中:|p|表示感兴趣区域(region of interest, ROI)的像素的幅度值。通常选择一个包含同质像素的ROI进行计算。

TCR的计算公式为

(11)

式中:|pTarget|表示目标像素的幅度值;|pClutter|表示杂波像素的幅度值。计算时选择一个包含强散射人造目标的ROI并标注目标真值,考虑到阴影效应,一般选择真值中前30%的像素作为目标像素pTarget参与计算,将ROI中真值以外的像素作为杂波像素pClutter。

FOM的计算公式为

(12)

3.1 基地大门数据(数据1)

本节以公开的基地大门视频SAR数据的第321～325开展对比实验,每一帧数据大小为720×640。

整个成像场景的原图和用5种滤波方法滤波后的图像如图6所示。图6(a)标记了6个大小为70×70的ROI,用于进一步评估滤波效果,其中包含2个匀质区ROI(红色框标记,分别记作ROI1、ROI2)、2个强散射区ROI(黄色框标记、分别记作ROI3、ROI4)和2个含边缘弱散射区的ROI(绿色框标记,分别记作ROI5、ROI6)。

图6 基地大门第323帧相干斑滤波结果对比Fig.6 Comparison of speckle filtering results of the 323rd frame at the gate

从全图来看,5种滤波方法都能抑制相干斑。其中均值滤波方法抑制效果有限且在强散射区域有明显模糊;improved Sigma方法相较均值滤波方法能有效抑制相干斑但存在较明显的模糊边缘现象;所提方法滤波、SAR-BM3D和V-BM3D则表现出相近的滤波性能,既能有效抑制相干斑噪声,又能很好地保持图像的细节;特别地,所提方法在匀质区的平滑效果更加明显。

图7展示的是放大的匀质区和强散射区ROI,用于评估所提方法在这两类区域的滤波性能。在匀质区ROI1,所提方法的滤波效果较为平滑;在匀质区ROI2,除均值滤波方法外,各方法的滤波效果较为接近。从图7(n)、图7(t)、图7(o)、图7(p)、图7(v)、图7(q)和图7(w)可以看到对比方法则存在不同程度的模糊强散射点的现象,而从图7(r)、图7(x)可以看到所提方法在包含强散射的目标区,既能有效滤除相干斑,又能保护强散射目标。

图7 数据1 ROI1-4相干斑滤波结果对比Fig.7 Comparison of speckle filtering results of ROI1-4 of data 1

为进一步定量评价所提方法的滤波效果,利用所选的6个ROI计算ENL、FOM、TCR指标。图8(a)、图8(b)为计算TCR指标所用掩膜,黑色区域像素表示杂波像素,白色区域表示目标真值中前30%的目标像素。图8(c)、图8(d)为用于计算FOM指标的边缘真值。图9、图10为包含边缘的ROI、边缘检测结果和二值化边缘图,其中边缘检测结果由ROA算法得到,二值化边缘图由检测阈值取0.8得到。可以看到,所提方法能在有效平滑弱散射区域的同时较好地保持边缘纹理,且滤波后能有效改善边缘检测结果。

图8 数据1真值图Fig.8 Ground truth of data 1

图9 曲线边缘区域检测对比Fig.9 Edge detection comparison for the curved edge region

图10 直线边缘区域检测对比Fig.10 Edge detection comparison for the linear edge region

定量计算结果如表1所示。从表1可以看出,各方法滤波后匀质区域的ENL指标都有提升。未滤波前,匀质区ROI1和ROI2的ENL指标为99.66和101.8;采用所提方法滤波后分别提升到372.2、833.77,优于其他4种对比方法。在强散射区域,各方法滤波后TCR指标都有不同程度下降,说明各方法对强散射目标存在不同程度的模糊作用。在强散射区ROI3,所提方法滤波后TCR指标由原图的15.95下降到14.86,比SAR-BM3D方法、V-BM3D方法的15.62、15.57低一些,而优于improved Sigma方法的14.38及均值滤波方法的13.87;在强散射区ROI4,所提方法滤波后TCR指标由原图的23.4下降到23.32,优于improved Sigma方法的23.05、SAR-BM3D方法和V-BM3D方法的23.27,以及均值滤波方法的21.06。在含边缘弱散射区ROI5,所提方法滤波后FOM指标由原图的0.25提升到0.63,略低于V-BM3D方法的0.70但优于SAR-BM3D方法的0.57、Improved Sigma方法的0.18;在含边缘弱散射区ROI6,所提方法滤波后FOM指标由原图的0.15提升到0.60,低于均值滤波方法的0.63和V-BM3D方法的0.68而优于Improved Sigma方法的0.34、SAR-BM3D方法的0.51。

表1 数据1相干斑滤波定量评价指标

综合ENL、TCR和FOM定量指标和目视定性评价,所提方法能在有效提高ENL的同时保持较高的TCR和FOM,有更好的相干斑滤波性能。

3.2 太阳能塔数据(数据2)

为进一步评估所提滤波方法的滤波性能,以太阳能塔视频SAR数据第61～65帧开展对比实验,每一帧数据大小为700×440。

整个成像场景的原图和用5种滤波方法滤波后的图像如图11所示。

图11 太阳能塔第63帧相干斑滤波结果对比Fig.11 Comparison of speckle filtering results of the 63rd frame at solar tower

选取6个大小为50×50的ROI用于进一步评估滤波效果,其中包含2个匀质区ROI(红色框标记)、2个包含反射镜的强散射区ROI(黄色框标记)和2个含阴影的弱散射区ROI(绿色框标记)。5种滤波方法能够有效抑制原图存在的相干斑噪声,但均值滤波、improved Sigma方法存在较明显的模糊现象。

为进一步展示滤波效果,放大的匀质区和强散射区ROI结果如图12所示。在匀质区,由于整体散射强度较低,图像整体偏暗,难以从目视的角度区分各滤波方法的性能差异。在强散射区,可以看到均值滤波、improved Sigma存在较明显的过模糊问题;所提方法和SAR-BM3D、 V-BM3D方法能较好地保持边缘。需要指出的是,所提方法将反射镜区域进行了强散射点保护后的统一平滑,滤掉了反射镜区域的一些细小纹理。

图12 数据2 ROI1-4相干斑滤波结果对比Fig.12 Comparison of speckle filtering results of ROI1-4 of data 2

利用所选的匀质区ROI1和ROI2计算ENL指标,利用强散射区ROI3和ROI4计算TCR指标,利用弱散射区ROI5和ROI6计算FOM指标。图13(a)、图13(b)为计算TCR指标所用掩膜,图13(c)、图13(d)为用于计算FOM指标的边缘真值。其中,图13(c)为ROI5的边缘真值,由2条太阳塔阴影边缘、4个点状目标边缘构成;图13(d)为ROI6的边缘真值,由2条目标边缘、3个反射镜阴影边缘构成。图14、图15分别为含太阳能塔和反射镜阴影的弱散射区的边缘检测结果,其中边缘检测结果由ROA算法得到,二值化边缘图由检测阈值取0.75得到。可以看到,在两个弱散射区,所提方法较好地保持了边缘纹理信息。

图13 数据2真值图Fig.13 Ground truth of data 2

图14 ROI3边缘检测对比Fig.14 Edge detection comparison for ROI3

图15 ROI4边缘检测对比Fig.15 Edge detection comparison for ROI4

定量计算结果如表2所示。从表2可以看出,各方法滤波后匀质区域的ENL指标都有提升。未滤波前,匀质区ROI1和ROI2的ENL指标为162.38和131.35;所提方法滤波后分别提升到494.80和361.49,优于其他4种对比方法,说明所提方法有着较优的相干斑抑制效果。在强散射区域,各方法滤波后TCR指标都有不同程度下降。由于强点保护机制,反射镜中有较多强散射点被保护起来,在强散射区所提方法取得较优的TCR指标。

表2 数据2相干斑滤波定量评价指标

在强散射区ROI3,所提方法滤波后TCR指标由原图的30.04下降到28.74,优于V-BM3D方法、SAR-BM3D方法及improved Sigma方法的27.91、27.86和26.38;在强散射区ROI4,所提方法滤波后TCR指标由原图的27.87下降到26.00,优于V-BM3D方法、SAR-BM3D方法及improved Sigma对比方法的25.35、25.28和22.07。在边缘保持方面,所提方法滤波后FOM指标有略微提升。在ROI5,FOM由原图的0.31提升到0.45,低于V-BM3D方法的0.57而高于SAR-BM3D方法的0.43,以及improved Sigma方法的0.23;在ROI6中,所提方法滤波后FOM指标由原图的0.43提升到0.63,低于V-BM3D方法的0.67而高于improved Sigma方法的0.12和均值滤波方法的0.48。从边缘检测结果图14(f)、图14(i)和图14(r),以及图15(f)、图15(i)和图15(r)可以看到,所提方法较好地保持了弱散射区的边缘纹理信息。

从两组视频SAR数据的滤波结果来看,所提方法在ENL、TCR、FOM等指标上综合性能较高,相干斑滤波效果优于对比方法。

4 结 论

本文提出了一种基于空时上下文协方差矩阵的视频SAR图像相干斑滤波方法。该方法利用视频SAR图像的空时上下文信息,通过构造空时上下文协方差矩阵来表征视频SAR图像距离-方位-时间三维动态散射信息,进而利用矩阵相似性检验方法和度量指标进行相似样本选取,最后根据相似样本实现相干斑滤波。使用公开的视频SAR数据开展对比实验。实验结果表明,所提方法能在较好地保持强散射点和边缘信息的同时有效抑制视频SAR中的相干斑,总体上得到较为可观的滤波结果。