一种基于相位补偿的干涉图自适应滤波算法

谭宏伟，周芳，贺丰收

(1.合肥工业大学 计算机与信息学院，合肥 230009；2.工业安全与应急技术安徽省重点实验室，合肥 230601；3.中国航空工业集团公司雷华电子技术研究所，江苏 无锡 214063)

0 引言

合成孔径雷达干涉测量技术[1-2](interferometric synthetic aperture radar，InSAR)利用覆盖同一地面区域的两幅合成孔径雷达图像生成干涉图，进而反演获得地面的高程或形变。InSAR具有观测范围大、测量精度高、可以全天时全天候对地面进行观测的优点[3-4]，近年来受到人们的广泛关注。然而，受各种去相干因素的影响，InSAR生成的干涉图相位中往往含有大量噪声，严重影响相位展开、高程重建等后续处理的性能[5-6]，因此，在干涉测量过程中需要对干涉图进行滤波。

目前，研究者已经提出多种干涉图滤波算法，这些算法大致可以分为两类：空间域滤波和变换域滤波。空间域滤波中多视滤波[7]是应用最早的一种算法，它通过将均值滤波模板在干涉图上逐像素移动卷积，从而获得每一像素的滤波结果。多视滤波的噪声抑制效果是以牺牲图像分辨率为代价的，不利于相位细节的保持。在多视滤波基础上，Trouve等[8]提出坡度相位补偿滤波算法，通过补偿滤波窗口内的线性坡度相位，对残余相位进行窗口大小自适应的多视滤波，从而能更好地去除噪声并降低相位信息损失。然而受线性相位模型的限制，当局部区域地形复杂或相位噪声严重时，坡度相位补偿滤波的滤波效果不佳。因此，Wang等[9]提出主相位成分补偿滤波算法，它通过对局部干涉图的频谱进行加权截断处理来获得反映地形轮廓的主要相位成分，然后对补偿了主要相位成分的残余相位进行自适应方向窗口滤波。主相位成分补偿滤波不受线性相位模型的约束，在一定程度上提高了滤波的效果，但该算法提取的主要相位分量存在过量噪声残留的情况，并且残余相位滤波时滤波像元数目固定，因此存在过滤波或欠滤波的情况。比较常用的空间域滤波算法还包括Lee滤波[10-12]、非局部均值滤波[13-14]等，这些算法在滤波处理时均存在一些问题，如Lee滤波中方向模板窗口不能适应所有的相位条纹模式，在干涉图相位变化复杂的区域会破坏条纹结构；非局部均值滤波计算复杂度高，且在相位噪声严重的区域难以有效估计像素间的相似性。

Goldstein滤波[15]是变换域滤波算法中最为经典的一种算法，它通过对局部干涉图的幅度谱进行平滑加权处理从而得到滤波器的响应。Goldstein滤波能有效保持相位的条纹结构，然而其噪声去除效果受滤波参数的影响较大，当参数设置不当时，干涉图会出现过度平滑或者大量噪声残留的情况。此后，学者们针对Goldstein滤波的不足提出一系列改进算法，如Baran等[16]基于干涉图局部相干性均值设置频谱加权参数的大小，Suo等[17]结合坡度相位补偿以及自适应频谱平滑窗口大小和频谱加权参数大小进行滤波。这些算法的滤波效果虽然较Goldstein滤波有所提升，然而均难以说明参数设置的最优性。干涉图小波滤波[18-21]能对干涉图像进行多分辨率分析，因此能有效避免分辨率损失。然而，小波变换不能精确地表示图像中的边缘、条纹信息[22]，因此不利于条纹结构的保持。

针对现有干涉图滤波算法存在噪声去除不充分或相位细节保持困难的问题，本文提出一种更加有效的自适应滤波算法。该算法通过对干涉图预滤波并对预滤波干涉图进行频谱截断的Goldstein平滑处理，从而提取更加精确的地形轮廓相位。在对补偿了地形轮廓相位的残余相位滤波时，基于像元相干系数设计了一种自适应像元数目的多视滤波方法，在有效去除噪声的同时充分保持相位细节。

1 基础理论分析

1.1 相位噪声统计特性

干涉图相位噪声的严重程度可以由生成它的两幅SAR图像的相干系数来衡量，相干系数越大，相位噪声越小。Franceschetti等[23]给出了相位噪声标准差与相干系数以及多视视数的关系，如式(1)所示。

(1)

式中：σφ为相位噪声标准差；γ为像元相干系数；L为多视视数；φ为干涉相位；φ0为相位的数学期望；ρ(φ，γ，L，φ0)代表相位的概率密度函数，其解析表达如式(2)所示。

(2)

式中：Γ(g)表示伽玛函数；β=γcos(φ-φ0)；F(g)表示高斯超几何函数。

结合式(1)和式(2)，可以计算不同相干系数以及多视视数情况下的相位噪声标准差。当视数为1、2、5、10、15、20时，噪声标准差随相干系数的变化曲线如图1所示。从图1可以看出，相位噪声标准差随着相干系数的增大而逐渐减小，并且增大多视视数可以有效降低噪声。

图1 相位噪声标准差与相干系数以及多视视数的关系

研究表明，当多视视数大于4且相干系数大于0.2时，相位噪声标准差的近似表示如式(3)所示[24]。

(3)

1.2 相位补偿滤波分析

相位补偿滤波中的干涉图相位模型如式(4)所示。

φz=φm+φr+v

(4)

式中：φz为测量相位；φm为地形变化引起的相位或反应地形轮廓的相位；φr为残余信号相位；v为噪声相位，v均值为零且与φm和φr相互独立。相位补偿滤波算法通过在局部窗口内估计并补偿地形变化相位或地形轮廓相位，然后对残余相位进行更加有效的滤波处理，从而能有效增强滤波的去噪效果并减少相位信息损失。

坡度相位补偿滤波假设局部窗口内的相位近似满足线性模型，即式(4)中地形变化相位φm可以表示为式(5)。

φm(x，y)=2π(xfx+yfy)

(5)

式中：(x，y)为像素在窗口内的位置坐标；fx和fy分别为方位向和距离向的相位频率。通过估计频率fx和fy可以得到线性坡度相位，进而对残余相位进行窗口大小自适应的多视滤波。

然而在坡度相位补偿滤波中当局部区域地形复杂时，干涉相位在滤波窗口内难以用线性模型近似，滤波会存在细节损失严重的问题。此外，当相位噪声较大时，相位频率fx和fy的估计精度会降低，滤波可能引入虚假相位条纹。残余相位滤波时，由于不同大小的窗口内像元数目不连续变化，因此坡度相位补偿滤波不能实现滤波强度自适应。在此基础上，文献[9]提出了主相位成分补偿滤波算法，算法在更大的局部窗口内通过对干涉图的频谱进行加权截断处理，获得反映地形轮廓的主要相位成分，避免了线性模型近似，并对残余相位进行自适应方向窗口滤波。主相位成分补偿滤波中主要相位分量的计算方法如式(6)所示。

φm(x，y)=arg{FFT2-1{S(u，v)Appc(u，v)}}

(6)

式中：arg{g}为求复数相位；FFT2-1{g}为二维逆傅里叶变换；(u，v)表示频率坐标；S(u，v)代表局部窗口内的干涉图频谱；Appc(u，v)代表频谱截断后的干涉图幅度谱。在文献[9]中，频谱截断的阈值设为半功率点。

然而，由于噪声分布在干涉图频谱的所有频带上[25]，因此主相位成分补偿滤波算法提取的主要相位分量存在过量噪声残留的情况。并且在实际数据处理时，频谱截断阈值设为半功率点时提取的主要相位分量往往无法反映干涉相位的轮廓结构，而且残余相位滤波时没有进行自适应处理，因此滤波存在噪声去除不足或相位细节损失严重的问题。

2 改进的相位补偿滤波算法

针对干涉图相位补偿滤波中存在的上述问题，本文提出一种改进算法，以进一步提高滤波的性能。提出算法将滤波过程分为以下四步。

1)对干涉图预滤波，降低相位噪声对干涉图频谱的影响。

2)将预滤波干涉图划分为一系列重叠的图像块，并采用频谱截断的Goldstein滤波方法提取每一图像块的地形轮廓相位。

3)从相应位置的噪声干涉图中去除图像块的地形轮廓相位，然后基于像元相干系数对残余相位进行自适应滤波。

4)将步骤3)中的滤波后相位与步骤2)中提取的地形轮廓相位相加，得到每一图像块的滤波结果，并对图像块的边缘重叠部分进行加权平滑处理。

算法流程如图2所示。

图2 滤波算法流程图

2.1 干涉图预滤波

噪声会影响地形轮廓相位的提取精度，因此有必要对干涉图预滤波。本文采用坡度相位补偿滤波与多视滤波相结合的方法进行预滤波处理。

坡度相位补偿滤波可以表示为式(7)。

(7)

(8)

式中：max{g}表示求函数最大值；P、Q为估计窗口的方位向和距离向像素半径；j为虚部；|g|表示求复数的模。

由1.2节分析可知，坡度相位补偿滤波在干涉图噪声严重的区域滤波效果较差，严重时可能因为相位频率的错误估计而引入虚假相位条纹。因此本文在低相干区域采用多视滤波方法进行预滤波处理，多视滤波计算方法如式(9)所示。

(9)

预滤波在初步去除噪声的同时应该避免过度模糊相位，因此本文设置较小的滤波窗口，如坡度相位补偿滤波的滤波窗口大小设为5像素×5像素，多视滤波的滤波窗口大小设为3像素×3像素，以充分保留相位信息。

2.2 局部区域地形轮廓相位提取

局部干涉图频谱中地形相位分量主要分布在较窄的频带范围内，而噪声分量则分布在所有频带上[25]，本文利用干涉图频谱的这一特点从预滤波干涉图中提取地形轮廓相位。

首先，将预滤波干涉图划分为一系列一定大小的图像块，并且为了避免滤波后相邻图像块间相位条纹不连续，本文采用与Goldstein滤波方法相同的策略，设置图像块间有一定的重叠，在滤波完成后，对重叠部分的滤波干涉图加权平滑处理。

对每一块预滤波干涉图进行二维傅里叶变换，得到其频谱S(u，v)。由于预滤波初步降低了相位噪声，因此信号频谱中混叠的噪声分量降低，更有利于地形轮廓相位的提取。为了进一步避免提取的相位中包含噪声，本文对预滤波干涉图进行频谱截断的Goldstein滤波处理，即只保留频谱中幅度超过阈值的窄带信号分量并对其加权平滑以进一步降低噪声的影响。提取的轮廓相位如式(10)所示。

φ′m(k，l)=arg{FFT2-1{SM(|S(u，v)|)T(S(u，v))}}

(10)

式中：SM(g)表示幅度谱平滑，通常采用3像素×3像素大小的均值窗口进行平滑处理；T(g)表示频谱阈值截断，如式(11)所示。

(11)

式中：b为频谱截断阈值。经过预滤波处理，地形相位分量的窄带特性进一步增强。该频谱会在较少频率上形成明显的峰值，因此本文将b设为幅度谱均值的3倍，即b=mean(|S(u，v)|)×3，经过测试提取的地形轮廓相位较为平滑且能精确地反映干涉相位的结构。

2.3 残余相位自适应滤波

将地形轮廓相位从噪声干涉图中去除，得到残余相位，如式(12)所示。

φR(k，l)=φr(k，l)+v(k，l)=arg{φ(k，l)e-jφ′m(k，l)}

(12)

由于残余信号相位φr主要是干涉相位的细节分量，因此残余相位滤波应能在有效去除噪声的同时充分避免细节损失。本文基于像元相干系数设计了一种自适应像元数目的多视滤波方法，避免滤波在低相干区域噪声去除不足或者在高相干区域相位细节损失严重。

由1.1节的分析可知，当相干系数大于0.2且多视视数大于4时，相位噪声标准差与相干系数以及多视视数间的关系可以表示为式(3)。因此，可以综合考虑噪声去除以及细节保持两方面因素，设定滤波后的噪声标准差值，根据此设定值就可以得到不同相干系数的像元应该采用的视数。另一方面，当像元相干系数很低时，需要较大视数才能将噪声标准差降至期望值，然而视数太大可能使相位损失严重，因此本文定义最大视数值，当视数达到最大视数时不再增加。当像元相干系数较高时，得到的视数可能较小，为了充分降低噪声，本文定义最小视数值，当视数降至最小视数时也不再减小。由此，可以得到多视视数，如式(13)所示。

(13)

式中：Nmax、Nmin、σ分别为设定的最大视数、最小视数以及滤波后的相位噪声标准差；round(g)代表四舍五入。

由于已经将轮廓相位从干涉图中去除，残余相位基本不含有相位趋势，因此可以从滤波像元邻域内选择相应视数的滤波样本进行多视处理。残余相位滤波的结果如式(14)所示。

(14)

式中：x为选择的滤波样本位置。

2.4 滤波结果

将提取的地形轮廓相位与滤波后的残余信号相位相加，即可得到滤波结果，如式(15)所示。

φ(k，l)=wrap{φ′m(k，l)+φ′r(k，l)}

(15)

式中：wrap{g}表示相位缠绕操作。

3 实验结果与分析

为验证提出算法的有效性，分别采用仿真干涉图和实测干涉图进行滤波实验，并将提出算法的滤波结果与Lee滤波、坡度相位补偿滤波、Goldstein滤波、文献[9]中的主相位成分补偿滤波以及文献[14]中的非局部均值滤波的滤波结果进行对比分析。

3.1 仿真干涉图滤波

在MATLAB中仿真两种具有不同相位变化特点的无噪干涉图，图像大小均为300像素×300像素，干涉相位如图3(a)和图3(c)所示。其中图3(a)的相位频率呈线性变化，频率从上到下由1/28逐渐增大到1/8，而图3(c)的条纹结构较为复杂，相位非线性变化严重。为验证算法在不同相干质量条件下的滤波性能，在图像上加入不同强度的相位噪声，使干涉图相干性从左到右由0.1增大到0.9，噪声干涉相位分别如图3(b)和图3(d)所示。

图3 仿真干涉相位图

分别采用本文算法以及几种对比算法对仿真的噪声干涉图进行处理，结果如图4和图5所示。在滤波参数设置时，本文将对比算法的参数大小设置为该类算法所采用的一般性大小，如Goldstein滤波、主相位成分补偿滤波的图像块划分尺寸设为32像素×32像素，重叠率设为50%，Goldstein滤波采用3像素×3像素大小的频谱平滑窗口，频谱加权参数α设为0.6。将本文算法与对比算法中具有相同含义的参数设置为一致，以观察相同参数下的滤波效果。对于本文算法独有的参数，本文则是参考相关文献并结合滤波实验进行的设置，如滤波后噪声标准差是参考文献[26]中的相关说明将其设为0.2 rad，最大视数、最小视数则是参考坡度相位补偿滤波中的最大窗口、最小窗口设置，将其分别设为81和9。

图4 仿真干涉图1滤波结果

图5 仿真干涉图2滤波结果

由图4和图5可以看出，各种滤波算法在仿真干涉图的高相干区域滤波结果相近，然而随着相干性的降低，滤波算法的性能出现较大差异。Lee滤波在低相干且相位变化较快的区域造成严重的条纹结构损失，这主要是由于Lee滤波中方向窗口不能适应各种条纹模式，并且噪声会影响方向窗口检测的准确性。坡度相位补偿滤波在中高相干性的图像区域滤波性能较好，然而当相干性较低时去噪效果较差。Goldstein滤波以及非局部均值滤波均能有效保持相位结构，然而其去噪性能较差，在图像左侧残留有大量的相位噪声。主相位成分补偿滤波在仿真干涉图1上滤波效果较好，然而在仿真干涉图2上存在低相干区域去噪不充分，干涉条纹复杂区域的相位结构损失严重的情况。综合来看，本文算法在不同相干性条件下均能有效去除相位噪声并保持相位的结构信息不被破坏，取得最好的滤波效果。

采用残差点减少率以及滤波相位相对于真实相位的均方根误差来客观评价滤波算法的性能，残差点减少率越大说明算法的噪声抑制效果越显著，均方根误差越小说明算法的相位保持效果越好。各种滤波算法的滤波性能分别如表1和表2所示，可以看出，本文算法在两幅干涉图上均具有最优的滤波性能。

表1 不同算法对仿真干涉图1滤波的效果评价

表2 不同算法对仿真干涉图2滤波的效果评价

3.2 实测干涉图滤波

实测干涉图由覆盖重庆地区的两幅ERS-1/ERS-2卫星数据生成，主辅图像成像日期分别为1996年4月4日和1996年4月5日，对应成像区域地形起伏剧烈，干涉相位具有条纹密集且噪声分布不平稳的特点。滤波干涉图大小为256像素×256像素，干涉相位如图6(a)所示，其对应的相干系数如图6(b)所示。分别采用本文算法以及五种对比算法对干涉图进行滤波处理，参数设置情况与前文一致，结果如图6(c)～图6(h)所示。

图6 实测干涉图滤波结果

从图6可以看出，几种对比算法在实测干涉图上的滤波表现与模拟干涉图滤波较为相似，Lee滤波造成严重的条纹结构损失；坡度相位补偿滤波、Goldstein滤波以及非局部均值滤波均有效保持了相位的条纹结构，然而它们在图像低相干区域滤波不足；主相位成分补偿滤波在一些区域改变了相位的结构，并且在低相干区域噪声去除不够充分。从目视效果来看，本文算法的滤波结果平滑度最好，恢复出的条纹结构也最清晰。

由于没有无噪干涉图作为参考，无法计算滤波相位的均方根误差。各滤波算法的残差点减少率分析如表3所示。从表3可知，本文算法的残差点减少率最大，噪声去除效果最好。

表3 不同算法对实测干涉图滤波的效果评价

4 结束语

本文提出一种基于相位补偿的干涉图自适应滤波算法。算法通过对干涉图预滤波以及频谱截断的Goldstein平滑滤波，解决了相位补偿滤波中地形轮廓相位提取不精确的问题，并且在残余相位滤波时基于像元相干系数进行滤波强度的自适应设置，从而使算法在噪声去除以及细节保持两方面取得平衡。实验结果表明，相比于五种对比算法，本文算法的残差点减少率最大，相位信息损失最小，具有最优的滤波性能。然而本文算法存在最优参数选择困难的问题，在实验部分参数的设置是根据经验进行的，如何从理论上得到自适应的最优参数仍然有待研究。