1、四川测绘第 28 卷第 1 期 2005 年 3 月 45 InSAR 系列讲座 4 I nSAR 数据处理及关键算法 * 刘国祥 ( 西南交通大学测量工程系， 四川 成都 610031 ) 摘要 作为 I nSAR 系列讲座的第四篇 , 本文介绍并讨论 InSAR 数据处理方法和关键算法 , 并对主要的数据处理流程以数学模型的形式进行了概括。 关键词 合成孔径雷达干涉 ; 数据处理 ; 算法 中图分类号 P237 文献标识码 A 文章编号 1001-8379(2005)01-0045-03 DATA REDUCTION AND KEY ALGORITHMS IN INSAR LIU Guo-
2、xiang (Dept. of Surveying Engineering, Southwest Jiaotong University, Chengdu 610031, China ) Abstract: As the fourth part of the tutorial, the paper will introduce and discuss the data reduction procedures and the key algorithms in InSAR. The primary data processing procedures are also represented
3、in the form of mathematical models. Key words: InSAR; data reduction; algorithms 1 引言 基于数字信号处理技术， InSAR的数据处理过程可以被高度自动化，以提取地表三维信息和地表形变结果。 前一讲座已经介绍了干 涉相位信号的构成、InSAR测高与探测形变的基本原理和方法，本讲将主要讨论 InSAR的数据处理流程及关键算法。 这里所涉及的干涉数据处理一般是从单视复数 SAR图像开始的，而有关 SAR数据处理的方法可参见文献 1。 在干涉数据处理实施之前，必须选择合适的干涉像对和其它辅助数据（如外部 DEM，用于地
4、形相位的去除）。 由上一讲可知，干涉像对的选择准则是：对 DEM生成来说，干涉基线既不能太长也不能太短；对于形变探测来说，干涉基线越短越好。 在得到有效的干涉数据集后，要对它们进行必要的处理，这些处理步骤包括 SAR图像配准、干 涉图生成、参考面/地形影响去除、几何变换、相位解缠等 2。 下面在描述这些关键处理步骤与算法的同时，也将对中间处理结果作一定的解释。 2 图像配准 从多时相的 SAR复数图像来提取地形起伏或地表形变信息，首要面临的问题便是将沿重复轨道（存在轻微偏移）获取的覆盖同一地区的图像进行精确配准。 SAR影像的配准就是计算参考影像（主影像）与待配准影像（从影像）之间的影像坐标映
5、射关系，再利用这个关系对待配准影像实行坐标变换和重采样。 因为轨道偏移量较小（一般在 1 km左右），而轨道高度为数百公里，因此，在重复轨道影像重叠区域内，同名像点对间的坐标偏移量具有一定的变化规律，一般可使用一个高阶多项式来拟合 2。 干涉测量要求影像配准精度必须达到子像元级。一般分两个阶段来实施，即粗配准和精配准。粗配准可利用卫星轨道数据或选取少量的特征点计算待配准影像相对于参考影像在方位向（影像列方向，以 t 坐标表示）和斜距向（影像行方向，以 r 坐标表示）的粗略偏移量，目的是为影像精确配准中的同名像素搜索提供初值。 而精配准首先是基于粗略影像偏移量和影像匹配算法，从主从影像上搜索出足
6、够数量的且均匀分布在重叠区域内的同名像点对，然后使用多项式模型来描述两 影像像素坐标偏差，即主从影像同名像点对的坐标差 （ ? r, ? t）可表示为主影像坐标（ r, t） 的函数表达式（ 多项式 ，如二次） ， rtbtbrbtbrbbtrtatarataraar5242321052423210+=+=DD（ 1） 基于所得到的同名像点坐标偏移观测量和最小二乘算法，多项式模型参数 ai 和 bi （ i=0,1, , 5）可以被求解出来，这样便完成了影像对坐标变换关系的建立。 最后利用这一模型对待配准影像进行重采样处理，使从影像取样到主影像的空间。重采样的可选方法较多，包括最邻近点法、双线
7、性内插法、双三次样条内插法等。 实验证明，双三次样条内 插法能46 四川测绘第 28 卷第 1 期 2005 年 3 月 满足干涉处理高精度的要求 2。 3 干涉图生成 将主影像与重取样后的从影像对应像素的相位相减，便可很容易地得到相位差图。 实际计算处理中，是先将主从影像作复数共轭相乘，其数学表达式为 2， = ),(),(),( trStrMtrI （ 2） 式中， M(r, t)和 S(r, t)分别代表主从图像对应像素的复数值， *表示复数共轭，而 I(r, t)表示所生成的干涉信息，也是复数值。 由此所产生的结果称为复数形式的干涉图（ Interferogram）。然后，从此干涉图中
8、提取相位主值分量图，即 可得到一次相位差图，注意，干涉相位在 - pi到 +pi内变化，一个完整的变化呈现为一个干涉条纹，但每一像素上存在相位整周模糊度问题。 为了反映干涉相位质量的高低，一般采用相关性测度 g 来衡量，实际计算中，利用主影像和取样后的从影像一定大小窗口内对应的多个（设为 N，如 N 5 5）像素值来估计，即， =Ni iiNi iiNi iiSSMMSM111g （ 3） 式中， Mi 与 Si 为复数值。 相关系数 g )1,0( 越大，干涉相位的噪声越低。 这样，我们可以得到一个与干涉相位图相对应的相关图（ Coherence map），以此来描述干涉相位的质量。 从计算
9、机显示效果来看，越亮的区域，相关系数越高，干涉相位的条纹越清晰，干涉相位观测量越可靠，反之，相关系数越低，干涉条纹越模糊，干涉相位噪声越大。 为了提高干涉相位的信噪比，我们可以采用空间滤波方法，也就是所谓的多视（ Multilooking）处理 2，即， = Liii SMv1（ 4） 上式是复数信号平滑的方法， 与实数信号的平滑方法（如加权平均）是有一定区别的。 4 参考面 / 地形影响去除 前已指出，一次差分干涉相位图是多种因素如参考趋势面、地形起伏、地表位移和噪声等方面的调和反映。对于地形测量来说，一般事先根据先验信息，选择不包含形变信息的干涉对来进行处理，以避免不必要的麻烦，因此，直接
10、相位差分值主要包含参考面（一般选择为参考椭球面）和地形起伏的贡献，为了使后续相位解缠变得容易，一般先将椭球参考面的相位分量从直接差分相位中去除。值得注意的是，相对于地形贡献来说，参考椭球面的贡献是占主导地位的，这就是为什么一次差分干涉相位图看起来呈现为大致与轨道相平行的条纹，有效干涉基线越长，干涉条纹越密集，地形坡度越大，干涉条纹越密集，地形越复杂，条纹曲率变化越明显。 当我们去除掉参考面的贡献后，地形相位条纹便清晰地显现出来，其表现形状与地形等高线的形状一致 2。 对于形变测量来说，尽管可以选择短基线的干涉对来处理，但参考椭球面和地形贡献的影响仍然是不可忽视的。为了提取形变信息，参考趋势面和
11、地形起伏的相位贡献必须从一次差分相位中去除，也就是所谓的二次差分。 所剩余的干涉条纹越密集，形变越显著剃度越大，反之，形变越弱且剃度越小 3。 4 .1 参考面干涉相位去除 在所选定的坐标参考框架系（如 WGS84）下，我们可以利用轨道（基线）数据来计算同名像点的参考斜距差 refRD 及参考面干涉相位，具体方法如文献 4 中所述。 当求出数量足够且均匀分布在重叠区域内的同名像点的多个参考斜距差之后，可使用 P阶（如 3 阶）多项式模型来描述这种参考斜距差大小在重叠区域内的变化规律 2，即， = =piijjjij, jiref tra(r, t)?R0 0（ 5） 这里的多项式系数可采用最小
12、二乘方法来求解。 一旦上述模型被建立，每一分辨元的参考椭球面干涉信 息 RI（复数值）便可采用如下公式来计算 2， += ),(4sin),(4cos)( trRitrRr, tRI refref DlpDlp （ 6） 最后依据下式，参考面相位可从初始干涉图中去除，以形成新的干涉图 I12。 ),(),(),(),(),(),(1trRItrStrMtrRItrItrI= （ 7） 4.2 地形干涉相位去除 为了提取形变信息，我们不仅需要 从干涉图 I1中去除参考面相位的贡献，而且还必须使用仅包含地形起伏信息的其它干涉图 2IT 进一步扣除地形的相位贡献。 对于“三轨”或“四轨”差分干涉方法
13、来说， 2IT 使用 一个 SAR 图像对来生成，且基于上讲式（ 12b）作比例化处理，而对于“两轨”方法来四川测绘第 28 卷第 1 期 2005 年 3 月 47 说， 2IT 使用外部 DEM（可采用数字地图来生成）结合轨道数据来模拟，具体的数学模型如上讲所述。设绝对意义下的地形干涉相位为 ),( trtopf ，则相应的复数值可表达为 2， ( ) ( ) ,(sin) ,(cos)(2 tritrr, tIT toptop ff += （ 8） 最后，形变干涉图 ID 可表达为 ),(),(),( 21 trITtrItrID = （ 9） 5 相位解缠 为了获得地表高程或 沿雷达斜
14、距方向上的地表位移量，我们必须确定干涉相位图（地形 I1、形变ID1）中每一像素的相位差整周数，这类似于 GPS 中的整周模糊度确定问题，在 InSAR 中称为相位解缠（ Phase unwrapping），是干涉数据处理中的难点和重点，也是干涉数据处理中的关键算法 5。 目前，相位解缠算法较多，但主要归为两类： (1)基于路径控制的积分法；和 (2)基于最小二乘的整体求解算法。 积分法的思路是：对缠绕相位图的每一像素，首先求其沿行向和列向的一阶差分，然后对一阶差分连续积分即可求得解缠相位。由于干涉相位图存在奇异点 （在复变函数里称为留数点），积分路径应受到约束以免局部干涉相位的误差传播，故这
15、种算法的关键是按一定的原则对奇异点定位并连接它们作为积分路径的“防火墙”，即积分时不能穿越这些路径。最小二乘算法的思想是：在解缠后的相位梯度与缠绕相位梯度差异平方和为最小的意义下整体求解，使用带权估计方法可削弱奇异相位对解缠结果的影响。 这两类算法的详细讨论见文献 5，然而，大量的实验已经证明积分法比最小二乘法更加稳健、可靠。 必须强调，相位解缠是一个十分棘手的问题。主要原因是三方面的 2，首先，尽管目前 SAR 成像系统的空间分辨率较 高，但是其均匀的取样率对一些特别的区域（呈现高频信号）来说是不够的，典型的例子如快速起伏的地形（悬崖、陡坡等）和过渡的地表位移区域（地震引起的断层破裂带），不
16、足的取样率导致相位信号失真（ Phase aliasing），这严重阻碍了这些特别区域内的连续相位重建。 其次，相位噪声会导致相位解缠失败，当相位噪声达到 p时，很难将它们与真实相位信号区分开来。 最后，干涉图可能包含一些“岛屿”，在这些岛屿内，因相位质量高，相位解缠能够顺利进行，然而，其它区域噪声太严重，致使相位解缠无法进行，在这种情况下，要将这些孤立的岛屿 内的解缠结果统一到同一系统中（即确定相对相位），没有先验信息的支持是无法完成的。典型实例如由河道隔开的区域、海岛区域、多植被覆盖地区等。 目前，相位解缠算法仍然是研究热点问题之一。 6 小结 我们已对 InSAR 中一些关键的数据处理步
17、骤和相关技巧逐一作了介绍，并给出了相应的数学模型。指出了相位解缠是数据处理难点的具体原因，并概括地总结和分析了已有的处理算法。 利用一次差分干涉图并结合雷达平台的姿态数据可大规模地重建 DEM。 对一次差分干涉图进一步扣除地形起伏的贡献，即所谓的二次差分干涉处理，则残余干涉相位图可用于监 测厘米级或更小的地球表面形变，以揭示许多地球物理现象并监测自然灾害，如地震形变、火山运动、冰川漂移、地面下沉以及山体滑坡等。 因此， InSAR 数据处理的焦点是干涉相位图， SAR 图像配准、干涉图生成与滤波、相位解缠、参考面相位去除、地形相位去除等都是关键的数据处理步骤和算法。 此外，由于雷达侧视成像几何
18、和因地形起伏引起的像点位移，我们还必须将具有几何变形的 SAR 图像及相关干涉结果图像作投影转换，以便得到正射投影产品并且纳入到地面坐标系统，也就是所谓的纠正处理。 参考文献 1Curlander, J.C. and R.N. McDonough, Synthetic Aperture Radar: Systems and Signal Processing. New York: John Wiley & Sons, Inc., pp.1-451, 1991. 2Liu, Guoxiang. Mapping of Earth Deformations with Satellite Radar
19、Interferometry: A Study of Its Accuracy and Reliability Performances. Ph.D. thesis: Dept. of Land Surveying & Geo-Informatics, The Hong Kong Polytechnic University, pp.250, 2003. 3Liu, G.X., X.L. Ding, Z.L. Li, Z.W. Li, Y.Q. Chen, SB Yu, Pre- and co-seismic ground deformations of the 1999 Chi-Chi, T
20、aiwan earthquake, measured with SAR interferometry. Computers and Geoscience, Vol.30, pp.333-343, 2004. 4刘国祥 , 丁晓利 , 李志林 , 陈永奇 , 章国宝 , 卫星 SAR 复数图像的空间配准 . 测绘学报 ， Vol.30, pp.60-66, 2001. 5Ghiglia, D. C. and M. D. Pritt, Two-dimensional Phase Unwrapping: Theory, Algorithms, and Software. New York: John Wiley & Sons, Inc, pp.1-200, 1998. 收稿日期 200 5 - 01 - 18 作者简介 刘国祥 ( 19 6 8 ) , 男 , 博士 , 现任西南交通大学测量工程系教授 , 主要从事遥感及 GIS 研究。 本讲座得到了国家自然科学基金项目的资助 (批准号：40374003)
