基于局部条纹特征的频率估计方法

基于局部条纹特征的频率估计方法

合成孔径干扰测量技术是利用合成孔内数据提取的相位信息为信息源获得三维信息和变化信息的技术。干涉合成孔径雷达/声纳通过两幅接收阵获取地表同一场景的复图像对,经共轭复乘后得到干涉图。干涉图中包含了斜距向上的点到两接收阵之间精确的距离差信息。但是,由于受系统噪声、去相干、配准误差和地表变化等因素的影响,实际获得的干涉图会存在不同程度的噪声,掩盖了真实值的同时,也影响了后续相位解缠的速度和精度。为了降低噪声的影响,需要对干涉图进行滤波。目前干涉图滤波方法有很多,滤波可以在空域或频域进行。空域滤波主要有圆周期均值滤波,圆周期中值滤波和Lee滤波等。空域滤波器实现简单,运算速度快,但容易破坏陡峭地形产生的高频干涉条纹。频域滤波的代表为Goldstein滤波,该算法利用干涉图自身的频谱幂指数进行频域加权实现滤波,对相位噪声较小的干涉图具有良好的降噪性能,但在相位噪声大的区域滤波效果下降。无论是空域还是频域滤波,都很少具有根据局部条纹特征调整滤波器的能力。近些年,一些文献提出利用局部频率估计的空频域自适应滤波算法,相比传统方法,这类算法对局部区域的滤波更有针对性,能够在保持干涉条纹的前提下更好地去除相位噪声。但存在的问题是:滤波性能依赖局部信噪比的大小。当信噪比较低时,可能产生错误的局部频率估计结果,进而致使相位滤波出现错误。为了解决局部频率估计类算法存在的问题,本文提出了一种改进的干涉图滤波方法。其主要思想是:在提取局部条纹方向和频率特征的基础上,利用方位和频率平稳的特征,对二维频率估计结果进行分析和判定,对认定为频率估计出现错误的点,通过适当调整估计窗口改善信噪比,实现对错误频率的修正,然后利用修正后的频率作为滤波器的参数完成滤波。1反映地形轮廓相位t-c假设单视主辅复图像s1和s2为s1=C+n1s2=Ce-jϕΤ+n2(1)s1=C+n1s2=Ce−jϕT+n2(1)式中:C为主辅图像的公共部分,n1和n2为热噪声,ϕT为真实的干涉相位。s1和s2经共轭相乘后得到的干涉信号可以表示为I=s1s*2=|C|2ejϕT+Cn*2+C*n1ejϕT+n1n*2(2)式(2)可以进一步写为I=ejϕc{|C|2ej(ϕT-ϕc)+e-jϕcEJϕT+N1N*2}(3)式中:ϕc为反映地形轮廓的相位,ϕT-ϕc为反映地形细节的相位,方括号内的项可视为噪声。根据式(3)滤波过程可分两步进行,首先提取反映地形轮廓的相位,然后对剩余地形细节相位进行提取。ϕc在局部区域内(窗口)可以近似认为线性不变,利用平面近似,即ejϕc≈ej(4)式中:m和n表示相对窗口中心点(x,y)的位置偏移,ϕ0为初始相位,可根据文献中的方法进行估计;fx和fy为点(x,y)对应的二维空间频率。在提取了地形轮廓对应的干涉相位后,再对细节相位ϕT-ϕc进行提取,可以通过计算相位均值的方法获得。反映地形轮廓相位ϕc的提取对于该滤波算法非常重要,因此二维频率fx和fy的参数估计就成为了滤波的关键。最常用的二维频率估计方法是最大似然法,它通过二维信号周期图的峰值位置来估计二维频率fx和fy,可以采用二维离散傅立叶变换(DiscreteFouriertransform,DFT)实现[ˆfx,ˆfy]=maxfx,fy{|Μ∑m=1Ν∑n=1Ι(x,y)e-j(mfx+nfy)|}(5)[fˆx,fˆy]=maxfx,fy{∣∣∣∣∑m=1M∑n=1NI(x,y)e−j(mfx+nfy)∣∣∣∣}(5)式中:M×N为选取的样本窗口尺寸,ˆffˆx和ˆffˆy分别为fx和fy的估计值。为了提高估计精度,需要对二维DFT进行插值。插值可以采用补零实现,但由于图像尺寸通常很大,直接补零插值计算量巨大。文献提出先利用DFT对局部条纹频率粗估计,然后再采用线性调频Z变换(Chirp-Ztransform,CZT)在粗估计值附近细化频谱的方法,大大减少了二维频率估计所需的计算量。2干涉图信噪比较低时频率估计错误的发生局部频率估计干涉图滤波利用局部相位特征进行滤波,滤波时很有针对性,滤波能力强。但不足的是,由于该滤波方法过于依赖频率估计结果,一旦局部频率估计出现较大误差或错误,就会破坏原有的干涉图,或在干涉图中引入额外的奇异点。局部频率估计的错误主要是由噪声引起的,当局部区域干涉图信噪比较低时,可能出现的错误情况有两种:一是频谱中相位噪声的峰值超过真实条纹所对应的频率峰值,如图1(a)所示。在采用频谱峰值位置估计条纹频率时,错误地选择了噪声峰值。另一种频率估计错误是由DFT频率粗估计产生的,如图1(b)所示。虽然干涉条纹对应的峰值比相位噪声的峰值高,但超过幅度并不大,当DFT的点数较少时,就可能出现因频率采样点少,导致频率采样后的噪声频率峰值相比条纹频率峰值高,产生频率粗估计错误,从而导致后续的CZT频谱细化的是噪声频率。虽然干涉图信噪比较低时容易发生频率估计错误,但并非完全无法克服。条纹信号的频率峰值仍然存在,只是被噪声峰值掩盖了。只要能够对局部频率估计错误值进行判定及识别,对频谱估计值进行修正,仍然可以得到精确的滤波结果。3条纹方向和条纹频率干涉条纹可以看成一种方向流场,并且其局部可认为是平稳的,这一特征与激光干涉条纹和指纹条纹的局部特征非常相似,可采用类似方法对局部干涉条纹进行建模。局部区域内的干涉条纹可以用个两个参数来描述:条纹方向和条纹频率。其中方向性是局部干涉条纹最重要的特征。相对于干涉条纹的频率,滤波对干涉条纹的方向更加敏感。若条纹方向错误,即使条纹频率估计正确,得到的滤波结果也将完全偏离真实值。条纹频率是局部干涉条纹重要的信息,决定滤波后相位误差大小。下文给出条纹方向和条纹频率的定义。假设图像中任意点(x,y)的相位为ϕ(x,y),则干涉条纹方向定义为相位场方向导数为零的方向,即∂ϕ∂l=∂ϕ∂xcosθ′+∂ϕ∂ysinθ′=Gxcosθ′+Gysinθ′=0(6)∂ϕ∂l=∂ϕ∂xcosθ′+∂ϕ∂ysinθ′=Gxcosθ′+Gysinθ′=0(6)式中:θ′为l到x轴的转角,称为条纹的方向角,求解后可写为θ′=arctan(Gy/Gx)-π/2(7)需要注意的是θ′取值范围为[0,π),而实际条纹方向取值范围应为),经修正后的方向角θ可表示为【MATH244Z】由于二维频率表示为【MATH245Z】因此方向图估计也可以用二维频率估计结果表示,即【MATH246Z】条纹的频率定义为沿条纹方向的相位变化率,可以直接计算二维频率矢量的模值得到频率估计【MATH247Z】4维频率估计的判别进行频率估计修正,首先需要对频率估计错误的像素点进行判定。由于局部区域干涉条纹是平稳的,条纹频率的变化比较平缓,局部区域干涉条纹方向和频率不会存在大的跳变。因此,可以通过考察相邻像素点的方向或频率的连续性,对当前像素点二维频率估计的可信度进行判断。可信度可以分为方向估计可信度和频率估计可信度,分别从方向图和频率图中获得。图2给出了一组仿真数据的干涉相位图、方向图和频率图。其中,图2(a)为真实相位为圆锥体的干涉相位图,信噪比为2dB,图2(b)和(c)分别为利用式(8)和(11)计算得到的方向图和频率图,部分像素点方向和频率估计存在明显的不连续,需要进行修正。本文以条纹方向和频率不连续程度作为测度提出一种可信度评价标准。假设像素点(x,y)的条纹方向和频率分别为θ(x,y)和f(x,y),则其方向和频率估计可信度Cd(x,y)和Cf(x,y)可表示为Cd(x,y)=L∑u=1Κ∑v=1∑u=1L∑v=1K|arg{exp[jˆθ(x,y)-jˆθ(x-u,y-v)]}|(12)Cf(x,y)=L∑u=1Κ∑v=1|ˆf(x,y)-ˆf(x-u,y-v)|(13)式中:u和v表示相对中心点(x,y)的位置偏移,L×K为相邻样本点数。当干涉图信噪比较高时,局部区域的干涉信号近似平稳,相邻点的条纹方向和频率基本保持不变,Cd和Cf值小,方向和频率估计可信度高。当干涉图信噪比较低时,局部区域的非平稳性增强,Cd和Cf值增大,方向和频率估计可信度降低。由于有可能出现频率相同方向不同或方向相同频率不同的情况,为了对频率估计错误点判定更加准确,可信度CA和Cf之间需要进行互补,这时干涉图中需要修正的区域R可以表示为R=Rd∪Rf(14)Rd={(x,y)|Cd(x,y)>DTd}(15)Rf={(x,y)|Cf(x,y)>DTf}(16)需修正区域R和非需修正区域ˉR可以视为二模态的,因此,最佳阈值DTd和DTf可利用最大类别方差法确定。下文从更直观的角度,以当前像素点及其四邻域点为例说明二维频率估计的判定过程。如图3(a)所示,中心灰色A0为待修正点,其邻点分别为A1、A2、A3和A4。当信噪比较高时,各点的二维频率估计是准确的,方向和频率变化都很小,A0频率估计的可信度小于阈值,无需修正。当信噪比较低时,由于局部噪声是非平稳的,空间位置不同点的噪声频率峰值也是不确定的,因此必然会导致A0与4个相邻点的频率估计值相关性变差,可信度大于阈值,需要进行频率修正。频率修正采用增大空间估计窗口长度的办法,通过将估计窗口由W扩大为W′,提高窗口内干涉条纹信噪比。需要注意的是,W′尺寸选择不能过大,以免破坏干涉图局部区域平稳的假设,同时也要确保窗口内信噪比能明显提高。另外,为减小由于窗口增大可能引入的频率估计误差的影响,可采取以下步骤:将W′窗口的频率估计结果(ˆfx′,ˆfy′)作为初步频率估计。在W窗口内的信号二维频谱中,搜索(ˆfx′,ˆfy′)附近的频谱峰值,如图3(b)所示,中间虚线框内为峰值搜索窗口,将最接近(ˆfx′,ˆfy′)的峰值所对应的频率作为最终估计频率(ˆfx,ˆfy)。若W窗口内所有的频谱峰值离(ˆfx′,ˆfy′)都比较远,则采用(ˆfx′,ˆfy′)作为A0区域的频率估计值。5研究结果与分析为了验证方法的滤波性能,本文采用几种经典的滤波算法进行比较。分别为圆周期均值滤波、Goldstein滤波和基于局部频率估计的地形自适应滤波。采用仿真数据和实际数据分别进行验证。仿真的理想干涉图如图4(a)所示。为了比较各算法在不同信噪比下的滤波性能,将干涉图分为上下两部分加入相位噪声,上半部分信噪比为9dB,下半部分信噪比为0.5dB,加入噪声后的干涉图如图4(b)所示。各滤波算法的样本窗口大小均为11×11,Goldstein滤波算法中的幂指数α取值0.5。图4(c)～(f)依次给出了采用圆周期均值滤波、Goldstein滤波、基于局部频率估计的地形自适应滤波以及本文方法的处理结果。由图4可见,在两种不同信噪比下,圆周期均值滤波对四周比较稀疏的条纹滤波可以取得比较满意的结果,但是中间条纹密集的区域效果不好,即使噪声较小的区域也影响了干涉条纹的致密性。Goldstein算法只在高信噪比区域有较好的滤波效果。基于局部频率估计的地形自适应滤波方法在低信噪比条件下,滤波性能较前两种方法有较大改善,但仍然存在明显的小块相位噪声,这是由于频率估计错误所致,采用本文方法可获得更优的滤波效果。为进一步检验本文方法的可靠性和有效性,采用航天飞机合成孔径雷达X-SAR对意大利Enta火山的实测数据进行验证,处理结果如图5所示。各滤波算法的样本窗口大小均为15×15,Goldstein滤波算法中的幂指数α取值0.9。从图5可以看出,圆周期均值滤波虽然可以显著去除相位噪声,但是在干涉条纹密集区域其致密性受到严重影响,Goldstein算法滤波后仍存在大量相位噪声,局部频率估计方法相比前两种算法,处理性能有明显改善,但在中间区域还是出现部分明显错误。而采用本文方法处理得到的干