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

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

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

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

x(k+1)=Ax(k)+Bη +w(k), (1)

E{w(k)}=0; E{w(k)w(j)}=Q(k)δ (k, j), (2)

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

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

图1

Fig.1

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

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

1)由式(3)(4)分别计算相位预测值 x^(k+1, k)和对应的方差阵P_{(k+1, k)}, 即

x^(k+1, k)= Ax^(k, k)+Bη _{(k, k)}(3)

P_{(k+1, k)}=AP_{(k, k)}A^T+Q_{(k, k)}(4)

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

2)由式(5)(6)分别计算相位状态估计值 x^(k+1, k+1)和对应的方差阵P_{(k+1, k+1)}, 即

x^(k+1, k+1)= x^(k+1, k)+J_(k+1)r_{(k+1, k+1)}(5)

P_{(k+1, k+1)}=(I-J_(k+1)C_{(k+1, k)})P_{(k+1, k)}(6)

式中: I为单位矩阵; J_(k+1)为卡尔曼滤波的增益矩阵; r_{(k+1, k+1)}为残差; C_{(k+1, k)}为线性化的观测矩阵。且

J_(k+1)=P_{(k+1, k)} CT(k+1, k)[C_{(k+1, k)}P_{(k+1, k)} CT(k+1, k)+(7)

r_{(k+1, k+1)}=y_{(k+1, k+1)}-C_{(k+1, k)} x^(k+1, k); (8)

C_{(k+1, k)}= ddxh(x) |x^(k+1, k)=[-sin x^(k+1, k), cos x^(k+1, k)]^T。 (9)

这里, R_{(k+1, k+1)}为观测噪声方差阵, 具体计算详见文献^[12]。

若式(8)中计算出的|r_{(k+1, k+1)}|> π , 则

r¯(k+1, k+1)=r_{(k+1, k+1)}-2π sign(r_{(k+1, k+1)})。 (10)

将 r¯(k+1, k+1)代入式(5), 以实现对残差的有效控制, 阻止误差的扩散。

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

x^(k+1, k)(m, n)= 12[ x^(k, k)(m-1, n)+ x^(m, n-1)]+ 12[η _m(m-1, n)+η _n(m, n-1)]; (11)

P_{(k+1, k)}(m, n)= 14[P_{(k, k)}(m-1, n)+P_{(k, k)}(m, n-1)]+M_wmQ(m-1, n) MTwm+M_wnQ(m, n-1) MwnT。 (12)

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

通过由步骤1)计算求得的相位预测值 x^(k+1, k)及对应的方差阵P_{(k+1, k)}, 进而完成步骤2)的计算, 最终求得SAR干涉图的解缠相位 x^(k+1, k+1), 实现相位解缠和噪声消除的同时处理。

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

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 真实数据相位图像Fig.2 Real data phase images

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

图3

Fig.3

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

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

图4

Fig.4

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

图5

Fig.5

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

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

Fig.6

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

3.2.2 定量分析

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

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

ε = 1MN∑i=0M-2 ∑j=0N-1ωi, jxϕi+1, j-ϕi, j-Δi, jxp+ 1MN∑i=0M-1 ∑j=0N-2ωi, jyϕi, j+1-ϕi, j-Δi, jxp, (13)

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

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

表1

Tab.1

表1(Tab.1)

表1 4种相位解缠算法的运算时间、不连续点数目和ε 值 Tab.1 Computation time、number of discontinuous point and ε value of 4 phase unwrapping methods 算法① 运算时间/s 不连续点/个 ε /rad 算法 运算时间/s 不连续点/个 ε /rad 原方法 19.0 280 0.386 04 本文方法 15.187 2 0.353 12 CZT方法 95.359 179 0.379 31 Snaphu方法 / 339 0.353 08 ①原方法、CZT方法和本文方法均在Windows XP操作系统下的Matlab7.0环境(计算机配置: CPU为Intel(R) Core (M) i3 CPU M370 @2.40GHz; 内存为2G)中运行; Snaphu方法在Ubuntu8.04下运行。

表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

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

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

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