在InSAR数据处理中，干涉图的质量直接关系到数字高程模型的高程精度和地表形变监测^[1]的精度，而干涉图往往受到热噪声、时空去相关、多普勒中心频率失相关^[2]等因素的影响，含有大量噪声，因此有必要对干涉图相位滤波，为后续相位解缠以及保证地表形变监测精度奠定基础。目前，用于干涉图滤波的方法有很多种，现有的干涉图像滤波方法主要分为频率域滤波、空间域滤波和时频域滤波。经典频率域滤波为Goldstein滤波^[3]，算法简单，得到了广泛应用。Goldstein算法虽然能滤除大部分噪声，但也会丢失大量边缘信息，给数字高程模型(DEM)和测量形变带来误差。文献[4-6]在Goldstein算法基础上又提出了改进算法，这些算法主要是调整了滤波因子，使干涉图滤波效果有所改善。空间域滤波有代表性的是Lee滤波。Lee滤波^[7]通过局部解缠，估计中心像素邻域内的坡度并改善中心像素，但局部解缠计算复杂、费时且目视效果差。此外，空间域滤波还发展了复数空间自适应滤波^[8]、最优方向融合滤波^[9]、数学形态法^[10]、EMD自适应滤波^[11]等方法，文献[8]在梯度计算时只考虑了水平和垂直方向，文献[12]进行了改进，有更强的自适应性。文献[9]的最优方向融合滤波结合相干性选择线性方向窗口滤波，增强了Lee算法的稳健性。数学形态法^[10]在Lee滤波中进行膨胀腐蚀得到边缘，增强了边缘信息，但强噪区滤波效果有待提高。EMD自适应滤波^[11]根据信号和噪声经验模态分解后表现的不同特征再进行自适应滤波，但尺度函数和IMF滤波个数的确定计算复杂。干涉图空间域和频率域滤波是以干涉相位的频率特性时不变或统计特性平稳为前提的，但实际干涉信号经常出现非平稳性^[9]，因此对其进行空间域或频率域滤波时得不到理想的滤波结果，可以把空间域与频率域分析结合起来，采用时(空)频分析方法对干涉图进行滤波处理，以便于提高精度，增强滤波效果。目前时频分析滤波方法主要小波滤波方法等。文献[13]提出矢量分离式小波软阈值滤波方法，能有效地保持干涉图中的相位信息。文献[14]提出对各个方向高频系数对应方向的边缘位置先平滑后再作中值滤波，很好地保持了图像的边缘特征。文献[15]提出采用小波维纳滤波对相位滤波，能更有效地滤除噪声，然其多次进行的小波维纳滤波，使滤波后的干涉图产生了块状效应。文献[16]将双树复小波用于干涉图去噪，去除了噪声还保持了干涉图的细节纹理信息，但在滤波函数模型中没有考虑复小波系数为复数的情况，丢失了虚部信息，导致信号相位噪声的增加影响滤波效果。

本文用复小波进行干涉图滤波，有其特殊的优越性。首先文献[17]已证明干涉相位噪声是加性的，而含加性噪声的干涉图经小波变换后，小波系数服从高斯分布，这决定了基于小波变换的干涉图滤波具有一定的优越性。其次是双树复小波(DC-CWT)除具有一般小波的优点外，还具有近似平移不变性、有限冗余性、良好的方向选择性和完全重构性等特点，具有提高精度和保留图像细节信息等优点，比实数小波滤波更有优势。

一种好的滤波方法除了变换的方法具备良好的性能外，还应该对变换系数有一个好的统计以及好的系数修改算法，而双变量贝叶斯估计算法正满足此要求。本文将复小波变换与双变量收缩函数相结合，给出了表达复噪声系数实部与虚部的相关性的概率密度函数，推导了复噪声系数概率密度函数、分解系数与其父节点系数的联合概率密度函数、复数域双变量收缩函数过程，给出了含噪信号小波复系数方差、噪声方差等，把双变量模型干涉图滤波函数模型从实数域推广到了复数域，得到了不错的效果。

1 复小波变换和复数域双变量模型干涉图滤波算法 1.1 复小波变换

复小波变换是通过两个或多个实小波变换平行对信号处理，一维双树复小波变换表示为^[18]

对于二维双树复小波变换

双树复小波在分解二维信号时，分别对二维数据的行和列进行变换，每一级分解出两个低频部分和6个高频细节部分，比二维实数小波变换多了3个方向，且由于干涉图中的噪声为加性噪声，含加性噪声的干涉图经小波变换后，其信号小波系数能量主要集中在低频子带，噪声小波系数主要反映在高频子带上，并且高频小波系数服从高斯正态分布，而双变量贝叶斯估计算法针对高频小波系数服从于高斯正态分布，给出了一个好的系数修改方法^[19]。因此，将双树复小波变换与双变量贝叶斯估计算法相结合是一种不错的滤波方法，该方法更适用于干涉图滤波，更有利于干涉图干涉条纹边缘等细节信息的保留与检测。

1.2 复数域双变量贝叶斯估计算法

经研究发现，图像经小波分解后，不同尺度间有很大的相关性，而双变量模型则在充分挖掘图像小波系数的尺度内、尺度间相关性的基础上，采用非高斯密度函数对“父子”系数分布情况进行建模，通过贝叶斯统计理论可求得原始影像的最大后验估计^[19]，能很好地刻画小波系数边缘分布形状和局部邻域系数幅值间相关性。该算法首先将含噪InSAR图像经多尺度双树复小波变换后分解成两个低频复数近似子图和若干高频复数方向细节子图，然后对各尺度高频复数方向子图利用双变量贝叶斯模型建模，建立复数双变量贝叶斯模型，并利用该模型对复系数进行贝叶斯最小均方估计，从而达到去噪目的。

对受噪声污染的影像经多尺度双树复小波变换后，得到变换后的复系数

设w₁是当前小波系数，w₂是上一尺度空间位置与w₁相同的系数，即w₁的“父”系数。用矢量表示为y=w+n，其中

式中，y、w、n分别代表含噪影像、原始影像和独立同分布高斯噪声的小波复系数矢量。

由于双树复小波的分解系数是复数，而双变量贝叶斯模型往往是在实数域的范围内使用，以下将实数域双变量滤波模型推广到复数域。

设p_w(w)、p_n(y－w)分别是真实系数w、噪声系数n的概率密度函数。在贝叶斯统计理论中，最大后验MAP估计即在给定y条件下，使后验概率密度p_w|y(w|y)取最大的w

利用贝叶斯准则，得到

由于P(Y≤y)=P(w+n≤y)=P(n≤y－w)，Y为随机变量，则p_y|w(y|w)=p_n(y－w)，进一步可得，。

对式(7)取对数可得

从式(8)可知，欲求原始影像复小波系数矢量w的估计值，必须知道噪声复系数的概率密度函数P_n(y－w)以及原始影像复小波系数的概率密度函数P_w(w)。

对于原始影像小波系数的概率密度函数P_w(w)，Sendur等通过对大量小波系数直方图中“父子”系数间概率分布的统计，提出了一种用非高斯的双变量概率分布函数来刻画小波分解系数与其父节点系数之间的相关性。

式中，w₂是w₁上一个尺度同位置处父节点系数。

由于w₁、w₂为复系数，考虑复数虚部的影响，w₁、w₂的概率密度函数为

干涉图经复小波分解后，复噪声系数n₁和n₂符合复高斯正态分布，设n₁=x₁+iy₁、n₂=x₂+iy₂，概率密度函数p_n1、p_n2给出了复噪声系数实部与虚部的相关性。n₁、n₂的概率密度函数分别为

利用噪声独立同分布假设，则n₁和n₂的联合概率密度函数

式中，n₂是n₁上尺度噪声复系数，σ_n²复噪声系数方差。

将矢量y、w、n以及式(13)代入式(8)，定义f(w)=lg[p_w(w)]，则

对于w(y)的极值，这里是复数的模的极值，即将复数函数取模之后，变成实函数，再用实函数求极值的方法，求出复值函数的极值。

式(14)的等价方程为

式中，分别代表 对复系数w₁、w₂的模的一阶导数

将式(17)、(18)代入式(15)、(16)，经简单推导

式(19)为复数域双变量收缩函数，考虑了复数的虚部对收缩函数的影响，增强了滤波效果。y₁、y₂为含噪影像的复小波系数，y₂是y₁下一级同位置处的父节点复系数。

此处设阈值 ，由函数模型y=w+n可知，σ_y²=σ²+σ_n²，σ_y²为含噪信号小波复系数方差，σ²为原始影像小波复系数方差，σ_n²为噪声方差。

由于含噪信号小波复系数均值为零^[20]，文献[17]证明符合高斯分布^[17]

M×M为正方形窗口w(k)大小。

噪声标准差在Donoho的鲁棒性中值估计方法的基础上

|y_i|是各层高频子带复系数的模。

则

1.3 算法实现

本文提出的复小波域复数双变量模型干涉图滤波的步骤如下：

(1) 取干涉相位的实部和虚部进行双树复小波分解，每层分解为2个低频子带和6个高频子带，6个高频子带分别对应图像中的6个方向(±15°，±45°，±75°)的复系数。

(2) 对每一个尺度下的6个高频子带的复系数按照式(19)—式(22)进行如下步骤的处理：①计算各尺度下复系数子带的噪声标准差，其算法按照式(21)计算，然后计算噪声方差;②按式(20)计算出各尺度下含噪信号小波复系数的方差，窗口大小为7×7;③按式(22)计算出各尺度下原始影像小波复系数方差;④按式(19)进行复小波系数的收缩，得到高频复信号系数的估计。

(3) 进行复小波逆变换，得到去噪后的复系数并生成干涉图。

2 干涉图试验结果及分析

为了验证算法的可靠性和有效性，分别采用模拟干涉图和真实干涉图将本文算法与Goldstein滤波、小波滤波、实数域复小波双变量滤波、最优化融合滤波进行了分析与比较。

2.1 模拟干涉图试验结果

采用模拟干涉图比较，相当于知道干涉图的真值，可以更有效的判断本文算法的性能。图 1为模拟干涉图各种方法滤波结果，其大小为592×592像素。Goldstein滤波的滤波参数根据经验选取为0.5，Sym4小波滤波、实数域复小波双变量滤波以及本文算法小波分解尺度皆为3层，实数域复小波双变量滤波(以下简称复小波双变量滤波)只考虑了小波复系数的实部，未考虑虚部，以便于与本文算法比较。其中图 1(d)采用的噪声方差、信号方差估计、双变量收缩函数只考虑复数的实部;图 1(e)该方法考虑了干涉图的相干性。从目视效果看，Goldstein滤波、复小波双变量滤波去噪效果不理想，去噪不均匀，部分区域存在大量噪声，局部地区的细节信息淹没在噪声中;最优化融合方向滤波和Sym4小波滤波去噪效果较好，去除了大部分噪声，细节信息保持较好，但在局部地区还存在明显的噪声斑点;复小波双变量滤波与本文算法相比，复小波双变量滤波由于忽略复小波系数的虚部，未考虑实部与虚部的相关性，导致信号相位噪声的增加影响滤波效果。本文算法考虑了复小波分解后小波系数虚部对干涉图的影响，整体去噪平滑，不存在明显噪声斑点，且细节信息保持良好，与未含噪相位图即真实值相接近。

在模拟的干涉图滤波结果中，选择具有代表性的第510行的相位值进行比较(干涉条纹既有密集的区域，又有疏松的区域)。从图 2中可以看出，Goldstein滤波剖面图和复小波双变量滤波剖面图含有很多毛刺，噪声未滤除干净;效果较好的最优方向融合滤波、Sym4滤波，与未含噪图像的剖面图相比，最优方向融合滤波在图中圆圈处存在剧烈震荡，Sym4滤波由于缺乏平移不变性在图中圆圈处即信号突变点附近存在较大的震荡，即伪吉布斯现象，而本文算法由于复小波具备的优点得到了很好的体现，既保留了细节信息又无离散小波的伪吉布斯现象，与未含噪图像剖面图最接近。

2.2 真实干涉图试验结果及分析

采用意大利Etna火山地区的真实干涉图进行研究，该干涉图利用ERS卫星分别于2000年9月6日和10月11日获取，选取范围为400×400像素的区域进行对比分析。从目视效果看，滤波使噪声大量减少，但Goldstein滤波、复小波双变量滤波在局部地区还存在大量明显噪声;效果较好的最优方向融合滤波和Sym4小波滤波去除了大部分噪声，在局部地区出现了纹理断裂现象;效果最好的本文算法，最大程度的消除了图像中的噪声，保持了图像的主要特征，图像纹理更清楚，噪声的抑制最明显，效果最优(图 3)。

2.3 干涉图滤波效果定量比较

以上是通过视觉效果定性评定干涉图滤波结果，下面通过残差点数、相位差和值(SPD)、相位标准差(PSD)指标定量评定干涉图滤波效果，在模拟干涉图中，由于知道干涉图真实值，还采用了均方误差(RMS)指标评定。

从表 1模拟干涉图滤波结果中看出，各种方法滤波后各评价指标均有所改善。均方误差(RMS)是利用未加噪的相位图即真值作为参考值，来计算滤波后的相位图和真值之间的偏差，其值越小代表滤波的保真性越好，越接近真值。本文算法RMS最小，滤波效果最好;从残差点的较少数、相位差和值(SPD)以及相位标准差(PSD)等指标来看，本文方法要优于其他算法。从表 2真实干涉图滤波结果看出，Sym4小波滤波和本文算法在各评价指标中表现最优，Sym4小波滤波残差点减少率达到94.31%，本文算法残差点减少率达到95.63%。

3 结 论

本文提出了复数域双变量模型干涉图滤波算法，该算法具有以下特点：①将复小波的双变量去噪模型从实数域推广到复数域，考虑了复数的实部与虚部对信号方差、噪声方差的影响，建立了复小波复数域双变量模型;②从干涉图滤波结果可以看出，二维双树复小波具有的多方向选择性，可提高图像滤波的精度和保留图像细节信息如边缘等，单小波由于缺乏平移不变性在信号突变点附近存在较大的震荡，而复小波由于具有的平移不变性使信号更为平滑;③对各种方法的滤波结果研究对比表明，无论是从目视效果还是定量评价指标来看，本文方法优于其他方法，残差点减少率高达90%以上，均方误差更小，更接近于真实值，达到了滤波的目的。