2. 电子科技大学 电子工程学院,四川 成都 611731
2. College of Electronic Engineering, University of Electronic Science and Technology of China, Chengdu 611731, China
1 引 言
干涉相位展开是干涉合成孔径雷达三维成像处理中的关键步骤,一直是干涉测量技术应用研究的热点和难点[1, 2]。利用多基线数据融合提高干涉相位精度是近些年来出现的一种新思想,它能克服单基线干涉SAR处理的缺陷,获得更为精确的高程地图(DEM),受到了越来越多的关注。在传统的干涉SAR三维成像中,干涉基线长度与高度模糊数成反比,长基线对应的干涉图条纹密集,能反映地形的精细结构,但增加了基线去相关,且条纹太密容易被噪声污染,将增加相位展开误差。另一方面,短基线对应的干涉相位图条纹稀疏,有利于相位展开,但却失去了地形的精细信息。多基线干涉SAR利用长短基线各自优点,根据长短基线之间的关系更加准确地估计出长基线干涉相位,既有利于干涉相位展开,又能减小高度模糊数,保留地形的精细结构[3, 4, 5, 6, 7]。过去十几年以来,有关多基线SAR干涉相位估计方法研究的文献[8, 9, 10, 11, 12, 13, 14]已陆续发表。
目前,较为有效的多基线干涉相位估计方法主要有最小二乘(LS)和最大似然估计(ML)等方法[6, 11]。其中LS算法能直接展开缠绕相位图,但如果权值选择不当或离散相位梯度估计不能反映真实的相位梯度,则会引入较大的误差。尽管ML算法没有LS算法中相邻像素点的干涉相位差小于π的假设条件[6],但它通常必须与其他方法结合才能获得最终的展开相位。如文献[6]中,ML算法只能获得折叠的相位,还需用其他相位展开算法对折叠后的相位进行展开。另外,文献[5]提出一种多基线迭代的相位展开方法,首先展开短基线干涉图以获得大致的DEM数据,再利用获得的粗略DEM数据估计较长基线对应的干涉相位,以此类推,最终获得最长基线对应的干涉相位估计值,迭代过程较为繁琐。文献[7]首先展开短基线对应的干涉图,其次根据短基线和长基线的关系,估计出长基线对应的干涉相位被2π模糊的倍数,再根据长基线对应的干涉相位的主值和估计出来的模糊倍数恢复出长基线对应的干涉相位。但该方法会把前一级误差扩大,如果扩大的误差影响到后一级干涉相位模糊倍数估计时,则会导致误差严重扩大而使得该相位展开后成为残差点。文献[8]把卡尔曼滤波(UKF)应用到多基线干涉相位估计之中,利用UKF强大的数据融合能力进行长基线干涉相位估计,具有较强的稳健性和较高的估计精度。但该方法需对每一干涉图进行滤波和展开处理,计算量较大,实时性不强。
本文提出一种多基线干涉相位估计方法,该方法选取适当的基线组合,用最大似然估计器[6]提取每一复像元随基线变化的频率,在展开短基线干涉相位基础上,最终获得长基线干涉相位估计值。该方法简单有效,误差传递较小,且不需要进行大量的迭代运算,甚至当短基线较短时,可不需要进行干涉相位展开。
2 多基线相位展开方法原理设多基线干涉中垂直轨迹基线长度为Bi(i=1,2,…,S),对应的基线倾角为ai(i=1,2,…,S)。则归一化的复干涉测量值可以表示为
式中,zi(k)为基线Bi所对应的干涉图k像元复干涉测量值;Φi(k)和Φ′i(k)分别为基线Bi所对应的干涉图k像元真实的以及含噪声的干涉相位值;S指干涉基线的数目;ϑi(k)和δi(k)分别为雷达视线角以及几何去相干等导致的干涉相位噪声;λ为雷达波长。在多基线干涉中,通常有ϑi(k)=ϑ1(k),ai=a1,因此,有 从而有 如果垂直轨迹基线长度Bi满足 则有 式中,d(k)=exp{-jB1sin[ϑ1(k)-a1]}和fk=-sin[ϑ1(k)-a1]ΔB是与i无关的;vk(i)=exp[jδi(k)]可以看做为复高斯噪声。于是,可以直接用最大似然估计器从序列y(k)=[y1(k) y2(k)…yS-1(k) yS(k)]中获得fk的估计值,于是Φi(k)可以通过下式求得只需要展开短基线B1对应的干涉图Φ1,就可以通过式(6)获得最长基线BS所对应的干涉相位ΦS(k)。
而基线间隔ΔB必须满足采样定理,即
需要注意的是满足式(7)要求的基线间隔非常小,因此很难满足实用要求。实际的干涉图由于平地效应的影响,干涉条纹密集不利于相位展开,通常需要先去除平地效应,再进行展开,且去平地效应之后的干涉图能大致反映场景的相对高度。InSAR几何关系如图 1所示,B为基线的长度,θ1为主天线入射角,雷达高度为H,α为基线倾角,R1和R2为主副天线到目标点P的斜距。
根据InSAR成像的几何关系,垂直雷达视线方向基线长度为
去平地效应后的干涉相位为
因此,基线Bi(i=1,2,…,S)去平地效应后的干涉相位可表示为
式中,h(k)和R1(k)分别是干涉图k像元所对应的场景目标点高度以及到主天线的斜距,它们大致范围通常是已知的;,同样可以用最大似然估计器获取fk的估计值。这时基线间隔ΔB应该满足 3 基线间隔均匀时的频率估计以去平地效应后干涉图为例来估计每一像素随基线变化的频率。去平地效应的复干涉图k像元可以表示为
式中对复干涉图k像元信号作关于变量i的傅里叶变换,并作取模运算,即可得
式中,F[x]表示取x的傅里叶变换。图 2即为图 3所示的干涉图对应的复干涉图像在(23,25)像素的归一化功率密度谱。本文使用最大似然频率估计器来估计每一像素随基线变化的频率fk。具体实现时,最大似然频率估计值通过搜索窗口内一维信号频谱最大值处的频率而得到。为了进一步提高估精度,可以通过给窗口内的信号末尾补零来增加频率域的采样数目。4 基线间隔非均匀时的频率估计
前述方法实质上相当于对干涉基线B以满足采样定理的方式进行均匀抽样,要求基线以相同的长度ΔB增加,条件比较严格,实现起来可能比较困难。因此,笔者也研究了基线间隔非均匀时的相位估计方法。
由前面的推导知道,去平地效应后复干涉图k像元可以表示为
式中,A(k)=exp[jΦ1(k)];·ΔB。当基线间隔均匀时,i取间隔一致的整数,可通过对式(14)作关于变量i的傅里叶变换来估计随基线变化的频率fk,从而获得长基线干涉相位估计值。然而当基线间隔不等时,i将取间隔不一致的数(部分是整数,部分是非整数),不能直接对其作傅里叶变换。针对式(14)中所示的这种单频信号,下面介绍一种用非均匀采样求其频率的方法。一般的,一维连续单频信号x(t)的频谱X(f)被定义为
由于不可能获得时间无限长的信号测量值,因此,在实际的处理中通常把时间限定在一定范围内,这里假设为信号时间宽度为T0。为了便于实际的计算,需要对连续的信号x(t)进行离散化处理,用离散的N0点采样值代替连续信号x(t)。对于非均匀N0点采样,用集合{t1,t2,…,tN-1,tN}代表采样时刻,则连续信号x(t)的频谱X(f)近似表示为[15]
式中,δ(t-tn)为脉冲函数,频谱X(f)是连续的,需要进一步离散化为通常采样点数目是非常有限的,这导致频率域采样间隔较大,不利于高精度的估计其主频。因此,为了进一步提高估精度,可以通过给窗口内的信号末尾补零来减小频率域采样间隔,假设补零后的信号采样点数目为Nd,则此离散化频率谱可进一步表示为
同样的,通过索窗口内一维信号频谱最大值的频率,从而得到其频率估计值。
5 仿真结果与分析 5.1 基线间隔均匀条件下仿真结果为了验证本文方法的有效性,对仿真多基线干涉SAR数据进行了处理。仿真基本参数如下:轨道高度为590 km,下视角为45°,波长为0.136 6 m,基线倾向角为10°,仿真场景为锥形场景(场景高度为390 m),地面分辨率为8 m×8 m。本文仿真中,最短基线为200 m,最长基线为550 m,基线间隔为50 m,相应的基线条数为8。
图 3(a)到图 3(h)为基线间隔为50 m时的缠绕相位图(图 3(e)、(f)省略),信噪比(SNR)为1.41 dB。图 4(a)为干涉基线长度为550 m时的真实干涉相位曲面。图 4(b)是枝切法[16]直接展开图 3(h)(即长基线干涉图)的结果,其展开相位误差及误差统计直方图见图 4(c)和图 4(d)。可以看出尽管枝切法展开相位误差范围较大,且主要分布在(0.2,0.5)区间,其平均展开相位误差为0.35 rad左右。
本文方法首先利用枝切法展开了图 3(a)所示的干涉图(即最短基线干涉图),其展开结果见图 5(a)。在展开最短基线干涉图的基础之上,利用最大似然估计器提取每一复像元随基线变化的频率,再利用公式(10)即可获得图 3(h)所示的长基线干涉图的展开结果,见图 5(b)。本文方法展开相位误差及误差统计直方图见图 5(c)和图 5(d),可以看出本文方法展开相位误差范围很小,且误差分布是以0为中心分布,其平均展开相位误差相当小。
另外,也用LS算法进行多基线相位展开,其展开相位图见图 6(a)。LS算法展开相位误差及误差统计直方图见图 6(b)和图 6(c),可以看出多基线LS算法展开相位误差范围比本文方法展开误差范围大,且相位误差主要分布在(0.1,0.55)区间,其平均展开相位误差为0.31 rad左右。为了定量分析相位展开算法精度,表 1列出了在不同信噪比条件下多基线LS算法、枝切法(直接展开长基线干涉图)以及本文方法的均方根误差。另外,表 1也列出了多基线ML算法展开误差,需要注意的是ML算法通常通过2基线或者3基线联合处理来估计较长基线的干涉相位,更多基线的联合处理是非常困难的(且暂时尚未发现有文献进行3基线以上的联合处理)。表 1列出的ML算法误差是ML算法联合处理B1基线干涉图、B2基线干涉图和B3干涉图的结果。
均方根误差/rad | ||||
SNR/dB | 枝切法 单基线 |
本文 方法 |
最小二乘法 多基线 |
最大似然 多基线 |
6 | 0.202 5 | 0.069 7 | 0.232 1 | 0.193 7 |
1.41 | 0.345 8 | 0.088 6 | 0.305 7 | 0.361 1 |
0 | 0.433 9 | 0.098 1 | 0.522 3 | 0.817 2 |
基线间隔均匀条件比较严格,现实起来可能比较困难。于是,对基线间隔非均匀时的相位估计也进行了研究,用本文方法实现了非均匀基线间隔条件下长基线干涉相位估计。仿真场景与主要参数与上一节相同,仿真基线长度分别为:B1=200 m、B2=260 m、B3=290 m、B4=350 m、B5=400 m、B6=440 m、B7=490 m、B8=550 m。图 7(a)为本文方法展开信噪比为1.41 dB的干涉图的结果,图 7(b)是展开误差图,图 7(c)为展开误差直方图。表 2列出了在基线间隔非均匀与均匀条件下本文方法展开精度,可以看出非均匀间隔条件下本文方法展开精度要略微高一些。这是因为非均匀频率估计具有更好的频率检测能力。
本文提出一种多基线干涉相位估计方法,并在多基线仿真数据处理中获得了较好的效果。与最小二乘多基线展开方法等相比,本文方法具有简单有效、计算量较小和误差传递较小的特点。但由于该方法需要用最大似然估计器提取每一像元随基线变化的频率,基线间隔不能太大,否则提取出来的频率是折叠的,这将引入严重的误差。因此用本文方法来估计较长基线干涉相位时通常需要的干涉基线条数较多,代价较高,但这对于多航迹重复轨道高精度三维测高系统来说是可以接受的。另外,笔者也正在研究增大基线间隔(亦即基线间隔不必满足采样定理)的方法,这也是下一步工作的研究重点。
[1] | XIANG Maosheng, LI Shukai. InSAR Phase Unwrapping Using Poisson’s Equation with Dirichlet Boundary Conditions[J].Acta Geodaetica et Cartographica Sinica, 1998, 27(3): 204-210. (向茂生,李树楷.用狄氏边界的泊松方程进行InSAR相位解缠的研究[J].测绘学报, 1998, 27(3): 204-210.) |
[2] | LIU Zilong, CAI Bin, DONG Zhen. A New InSAR Phase Unwrapping Method Based on Local Frequency Estimation and Multigrid Technique[J].Acta Geodaetica et Cartographica Sinica, 2010, 39(1): 82-87. (刘子龙,蔡斌,董臻.基于局部频率和多网格技术的InSAR相位解缠算法[J].测绘学报, 2010, 39(1): 82-87.) |
[3] | LOMBARDINI F, LOMBARDO P. Maximum Likelihood Array SAR Interferometry[C]//IEEE Digital Signal Processing Workshop Proceedings. Loen: IEEE, 1996: 358-361. |
[4] | LOMBARDO P, LOMBARDINI F. Multi-baseline SAR Interferometry for Terrain Slope Adaptivity[C]//Procee-dings of the 1997 IEEE National Radar Conference. Syracuse: IEEE, 1997: 196-201. |
[5] | THOMPSON D G, ROBERTSON A E, ARNOLD D V, et al. Multi-baseline Interferometric SAR for Iterative Height Estimation[C]//IGARSS’99 Proceedings: IEEE 1999 International Geoscience and Remote Sensing Symposium: 1. Hamburg: IEEE, 1999: 251-253. |
[6] | ZHANG Qiuling, WANG Yanfei. Improving the Interfero-metric Phase Accuracy of Distributed Satellites InSAR System with Multibaseline Data Fusion[J]. Journal of Electronics and Information Technology, 28(11): 2011-2014. (张秋玲,王岩飞.利用多基线数据融合提高分布式卫星InSAR系统的干涉相位精度[J].电子与信息学报, 2006, 28(11): 2011-2013.) |
[7] | LIU Hui, ZHOU Yinqing, XU Huaping. An Estimate Method for Multi-baseline INSAR Phase[J]. Journal of Astronautics, 2008, 29(6): 1992-1994. (刘慧,周荫清,徐华平.多基线干涉SAR的相位估计[J].宇航学报, 2008, 29(6): 1992-1994.) |
[8] | XIE Xianming, PI Yiming. Multi-baseline Phase Unwrapping Algorithm Based on the Unscented Kalman Filter[J]. IET Radar, Sonar & Navigation, 2011, 5(3): 296-304. |
[9] | SHI X J, ZHAN Y H. A Multi-baseline Data Fusion Algorithm for Distributed Satellites SAR Interferometry by Combining Iterative and Maximum-likelihood Methods[C]//Proceedings of 2009 Asia-Pacific Conference on Synthetic Aperture Radar(APSAR’09). Xi’an: [s.n.], 2009: 165-168. |
[10] | XU W, CHANG E C, KWOH L K, et al. Phase Unwrapping of SAR Interferogram with Multi-frequency or Multi-baseline[C]//1994 International Geoscience and Remote Sensing Symposium (IGARSS’94). Piscataway:[s.n.], 1994: 730-732. |
[11] | GHIGLIA D C, WAHL D E. Interferometric Synthetic Aperture Radar Terrain Elevation Mapping from Multiple Observations[C]//Proceedings of 6th IEEE Digital Signal Processing Workshop. Yosemite National Park: [s.n.], 1994: 33-36. |
[12] | KIM M G, GRIFFITHS H D. Phase Unwrapping of Multi-baseline Interferometry Using Kalman Filtering[J]. Seventh International Conference on Image Processing and Its Applications. Manchester: IEEE, 1999: 813-817. |
[13] | LI Z F, BAO Z. A Joint Image Coregistration, Phase Noise Suppression, and Phase Unwrapping Method Based on Subspace Projection for Multi-baseline InSAR Systems[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(3): 584-591. |
[14] | LACHAISE M, EINEDER M, FRITZ T. Multi Baseline SAR Acquisition Concepts and Phase Unwrapping Algorithms for the TanDEM-X Mission[C]//2007 International Geoscience and Remote Sensing Symposium (IGARSS’07), Barcelona: [s.n.], 2007: 5272-5276. |
[15] | GAO Yukai, ZHANG Wei. Spectrum Estimation from Non-uniformly Sampled Signals[J]. Electrical Measurement & Instrumentation, 2009, 46(2): 8-11. (高玉凯,张维.非均匀采样信号的频谱分析方法[J].电测与仪表, 2009, 46(2): 8-11.) |
[16] | GOLDSTEIN R M, ZERBER H A, WERNER C L. Satellite Radar Interferometry: Two-dimensional Phase Unwrapping[J]. Radio Science, 1988, 23 (4):713-720. |