基于DEM的低相干区SAR干涉图卡尔曼滤波相位解缠算法

郝华东，刘国林，陈贤雷，曹振坦

(1.舟山市质量技术监督检测研究院计量检定测试中心，舟山 316021;2.山东科技大学测绘科学与工程学院，青岛 266510)

0 引言

InSAR是一种极具潜力的微波遥感技术，在地形测绘、地震形变、火山运动、地面沉降以及滑坡监测等方面都具有广泛应用。然而，在星载或机载InSAR遥感图像上，由于雷达热噪声、阴影、几何重叠、时间去相关和相位失真等因素引起的噪声和伪信号往往导致干涉相位的不连续，进而引起局部相位误差[1]，这使得相位解缠成为国内外相关学者的研究热点。相位解缠算法主要有路径跟踪算法[2-4]、最小范数算法[5-7]、基于最优估计方法[8-13]以及基于区域分割和编码的相位解缠算法[14]等。卡尔曼滤波相位解缠算法作为一种基于最优估计的方法，通过将相位的解缠问题转化为相位的状态估计问题，运用卡尔曼滤波技术，最终实现相位解缠最优解的求取。其优点是直接对SAR复干涉图像进行卡尔曼滤波处理，克服了传统相位解缠算法在解缠之前需要进行预滤波处理的弊端，能够同时实现滤波和解缠，简化了数据处理过程。

相位解缠的卡尔曼滤波状态模型一般假设相位的状态变化是一个缓慢的过程。在地形平坦或起伏不大时，常规卡尔曼滤波相位解缠算法通常能够得到可靠的解缠结果，但在地形陡峭或起伏较大时，则因相位的状态变化较快，导致解缠结果存在误差传递[12]。因此，地形起伏地区的 SAR干涉图解缠需要充分考虑地形因素的影响，通过引入相关的地形信息来提高解缠的精度。文献[12]提出了一种顾及地形因素的卡尔曼滤波相位解缠算法，该算法采用线性调频Z变换(chirp z transform，CZT)进行局部频率估计来实现对地形的自适应模拟，但该方法的精度不高且计算效率偏低。为此，本文提出了一种基于SRTM DEM的卡尔曼滤波相位解缠方法，并通过采用InSAR实际图像数据对该算法的可行性和有效性进行了验证。

1 基于地形的卡尔曼滤波状态空间模型

卡尔曼滤波的观测方程以及相关参数的定义详见文献[12]。这里只介绍基于地形的卡尔曼滤波状态空间模型。

当SAR干涉相位在离散时间情况下，考虑地形影响的卡尔曼滤波相位解缠状态空间模型为

式中:A为给定状态空间模型的系统矩阵;B为输入控制矩阵;η为与地形有关的输入控制变量;w(k)为状态噪声;E{}表示数学期望;Q(k)是状态噪声方差阵;x(k)为k点处的真实相位。这里，将SRTM DEM引入与地形有关的输入控制变量η。在引入之前，需要对SRTM DEM进行一系列处理，首先填充SRTM DEM上的空白区，然后与SAR主图像配准，再模拟相位，最终将其结果应用于状态空间模型中。

2 基于SRTM DEM的卡尔曼滤波相位解缠算法实施过程

基于SRTM DEM的二维卡尔曼滤波相位解缠算法的基本流程如图1所示。

图1 基于SRTM DEM的卡尔曼滤波相位解缠算法流程Fig.1 Flow chart of Kalman filter phase unwrapping based on SRTM DEM

具体计算步骤[13]如下:

式中:矩阵A和B已在式(1)中定义;η(k，k)为局部频率估计;Q(k，k)为状态噪声方差阵。这里，矩阵A和B均取单位矩阵。设和 P(1，1)分别为初始相位以及对应方差阵(根据经验设定)。

式中:I为单位矩阵;J(k+1)为卡尔曼滤波的增益矩阵;r(k+1，k+1)为残差;C(k+1，k)为线性化的观测矩阵。且

这里，R(k+1，k+1)为观测噪声方差阵，具体计算详见文献[12]。

若式(8)中计算出的|r(k+1，k+1)|＞π，则

在对SAR干涉图进行相位解缠时，由于任一估计相位都是通过2个相邻区域计算得到的，所以在卡尔曼滤波相位解缠算法中采用式(11)(12)进行估计，即

式中:Mwm=[0.5，0];Mwn=[0，0.5];ηm(m-1，n)和ηn(m，n-1)分别表示方位向和距离向上的局部频率估计，这里采用的是SRTM DEM模拟相位在方位向和距离向上的差分代替局部频率估计。

3 实验结果分析

以2003年6月和2004年1月伊朗Bam地区2景ENVISAT SAR图像为数据源，从图像上选取地形起伏较大、细节较丰富的150像元×150像元大小的区域(图2(a))，其相干图如图2(b)所示(图中矩形框区域中的黑色部分噪声污染严重，相干系数在0.1左右)。经过与SAR主图像配准的SRTM DEM模拟得到的相位图如图2(c)所示。

图2 真实数据相位图像Fig.2 Real data phase images

运用原卡尔曼滤波相位解缠方法(简称原方法)[11]、CZT顾及地形因素卡尔曼滤波相位解缠方法(简称CZT方法)[12]、本文基于SRTM DEM 的卡尔曼滤波相位解缠方法(简称本文方法)和目前较为流行的Snaphu网络流相位解缠算法[8]分别进行解缠，解缠结果如图3(a)—(h)所示。

图3 真实数据解缠相位图像(rad)Fig.3 Real data unwrapped phase images(rad)

图4 和图5分别表示上述前3种方法在方位向 和距离向的频率估计结果。

图4 真实数据在方位向的频率估计结果Fig.4 Results of real data frequency estimation in azimuth direction

图5 真实数据在距离向的频率估计结果Fig.5 Results of real data frequency estimation in range direction

3.1 频率估计结果分析

以图4、图5中真实数据的方位向和距离向频率估计结果对比来看，原方法中固定分析窗口获得的结果虽然平滑，但是在噪声污染严重的局部区域出现了错误的估计结果。无论是方位向(图4(a))还是距离向((图5(a))的结果图上都存在大量的残差点;CZT方法(图4(b)和图5(b))在平滑度上有了一定的提高，但是仍然存在一些残差点;而本文方法(图4(c)和图5(c))无论在平滑度和估计精度上都有明显的优势，在噪声污染严重的区域(图2(b)中矩形框区域)能够获得相对平滑的估计结果。

3.2 相位解缠结果分析

3.2.1 定性分析

从图3可以看出，原方法的解缠结果(图3(a)(e))有明显的误差传递，且丢失部分边缘信息，“拉线”现象较严重;CZT方法结果(图3(b)(f))虽然对地形进行了有效估计，但是仍然有误差传递;而本文方法结果(图3(c)(g))更加平滑，说明图2(b)中所框的噪声污染严重区域已经成功解缠，几乎没有误差传递，与Snaphu方法的解缠结果(图3(d)(h))比较接近。从解缠重缠绕结果(图6)上看，Snaphu方法(图6(d))中噪声仍然存在，故Snaphu方法不能实现对SAR干涉图的滤波;但是卡尔曼滤波相位解缠算法能够在对SAR干涉图解缠的同时实施滤波，而本文方法(图6(c))滤波效果也好于原方法(图6(a))和CZT方法(图6(b))，噪声几乎完全被去除。这充分说明本文方法在解缠和滤波方面的优势。

图6 解缠重缠绕结果(rad)Fig.6 Rewrapped phase images of unwrapping result(rad )

3.2.2 定量分析

在对解缠结果进行定性分析的基础上，运用定量分析的方法，分别从运算时间、不连续点数目和解缠质量ε值[15]3个方面评价本文基于SRTM DEM方法的性能:运算时间越短，表明算法计算效率越高;不连续点数目越少，说明抗相位畸变的性能越好;ε值越小，则解缠质量越高。

这里给出评价解缠结果质量的ε值表达式[15]为

式中:M和N分别为图像的列数和行数;φ(i，j)是点(i，j)处的解缠相位;分别是与缠绕相位梯度相对应的权重，此权重一般由相位质量图给出，其值一般在0～1之间;p的取值通常有0，1，2，这里p取1。本文实验中以相干图作为相位质量图。

表1给出了原方法、CZT方法、本文方法和Snaphu方法的运算时间、不连续点数目和ε值。

表1 4种相位解缠算法的运算时间、不连续点数目和ε值Tab.1 Computation time、number of discontinuous point and ε value of 4 phase unwrapping methods

图7给出了原始缠绕相位及原方法、CZT方法、本文方法和Snaphu方法的不连续点分布图，其中原始缠绕相位的不连续点数目约为2 239个。

图7 不连续点分布图(白色为不连续点)Fig.7 Discontinuity maps (white is discontinuous point)

通过表1和图7可以看出，本文方法不仅运算时间和ε值均小于原方法和CZT方法的，而且不连续点数目最少(仅剩2点);在与Snaphu方法计算出的ε值相当的情况下，本文方法得到的不连续点数目相对更少，表明本文方法抗相位畸变性能强，解缠质量相对较高。

4 结论

本文针对现有卡尔曼滤波相位解缠算法在地形起伏较大或噪声高时无法准确反演地表形变信息这一情况，通过在卡尔曼滤波状态空间模型中引入SRTM DEM模拟相位来估计地形因素对解缠的影响，提出了一种基于SRTM DEM的卡尔曼滤波相位解缠算法。该算法通过利用SRTM DEM地形信息指导解缠，提高了运算的速度和精度。实际数据解缠结果表明，该算法的运算效率和解缠质量有了明显的提高，尤其能够提高低相干区域的解缠精度。

[1]张 红，王 超，刘 智.星载合成孔径雷达干涉测量[M].北京:科学出版社，2002:100-107.Zhang H，Wang C，Liu Z.Spaceborne synthetic aperture Radar interferometry[M].Beijing:Science Press，2002:100-107.

[2]Goldstein R M，Zerber H A，Werner C L.Satellite Radar interferometry:Two-dimensional phase unwrapping[J].Radio Science，1988，23(4):713-720.

[3]Xu W，Cumming I.A region-growing algorithm for InSAR phase unwrapping[J].IEEE Trans Geosci Remote Sens，1999，37(1):124-134.

[4]Flynn T J.Two-dimensional phase unwrapping with minimum weighted discontinuity[J].J Opt Soc Am A，1997，14(10):2692-2701.

[5]Ghiglia D C.Romero L A.Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods[J].J Opt Soc Am A，1994，11(1):107-117.

[6]Pritt M D，Shipman J S.Least-squares two-dimensional phase unwrapping using FFTs[J].IEEE Trans Geosci Remote Sens，1994，32(3):706-708.

[7]Pritt M D.Phase unwrapping by means of multigrid techniques for interferometric SAR[J].IEEE Trans Geosci Remote Sens，1996，34(3):728-738

[8]Chen C W.Statistical-cost network-flow approaches to two-dimensional phase unwrapping for radar interferometry[D].California，USA，Stanford:Stanford University，2001.

[9]Krämer R.Auf kalman-filtern basierende phasen-und parameter estimation zur lösung der phasen vieldeutigkeits problematik bei der Höhenmodellerstellung aus SAR-interferogrammen[D].Germany，Siegen:Univer-sität-GH Siegen，1989.

[10]Loffeld O，Nies H，Knedlik S，et al.Phase unwrapping for SAR interferometry—a data fusion approach by Kalman filtering[J].IEEE Trans Geosci Remote Sens，2008，46(1):47-58.

[11]刘国林，郝华东，陶秋香.卡尔曼滤波相位解缠及其与其他方法的对比分析[J].武汉大学学报:信息科学版，2010，35(10):1174-1178.Liu G L，Hao H D，Tao Q X.Kalman filter phase unwrapping algorithm and comparison and analysis with other methods[J].Geomatics and Information Science of Wuhan University，2010，35(10):1174-1178.

[12]刘国林，郝华东，闫 满，等.顾及地形因素的卡尔曼滤波相位解缠算法[J].测绘学报，2011，40(3):283-288.Liu G L，Hao H D，Yan M，et al.Phase unwrapping algorithm by using Kalman filter based on topographic factors[J].Acta Geodaetica et Cartographica Sinica，2011，40(3):283-288.

[13]郝华东.卡尔曼滤波在InSAR相位解缠中的应用研究[D].青岛:山东科技大学，2010.Hao H D.Study of Kalman filter applied in phase unwrapping of InSAR[D].Qingdao:Shandong University of Science and Technology，2010.

[14]詹总谦，钱 俊，舒 宁.基于区域分割和编码的相位解缠方法[J].武汉大学学报:信息科学版，2002，27(3):316-320.Zhan Z Q，Qian J，Shu N.A method of phase unwrapping based on region-segmentation and encoding[J].Geomatics and Information Science of Wuhan University，2002，27(3):316-320.

[15]Ghiglia D C，Pritt M D.Two-dimensional phase unwrapping:Theory，algorithms and software[M].New York:John Wiley ＆ Sons Inc，1998.