InSAR变形监测方法与研究进展

朱建军, 李志伟, 胡俊     

中南大学地球科学与信息物理学院, 湖南 长沙 410083

基金项目：国家重点基础研究发展计划（973计划）（2013CB733303）；国家自然科学基金（41371335；41404011）；高分辨率对地观测系统重大专项（民用部分）（03-Y20A11-9001-15/16）

第一作者简介：朱建军(1962-), 男, 教授, 博士生导师, 研究方向为测量平差与InSAR数据处理。E-mail:zjj@mail.csu.edu.cn

星载合成孔径雷达干涉测量(interferometric synthetic aperture radar, InSAR)因其高精度、高分辨率、全天候等优点已迅速成为常用的大地测量技术之一，旨在通过计算两次过境时SAR影像的相位差来获取数字高程模型[1-2]。随之而来的差分雷达干涉技术(differential interferometric synthetic aperture radar, D-InSAR)则是通过引入外部DEM或三轨/四轨差分实现了地表变形监测。除了地形因素外，时空失相干、大气、轨道等也是影响DInSAR技术的形变监测精度的主要因素[3-4]。与此同时，卫星硬件的不断发展，重返周期越来越短会导致微小形变与噪声之间的混叠，空间分辨率越来越高也会导致噪声更加复杂，这些都对常规D-InSAR的形变监测提出更高的挑战。此外，D-InSAR技术获取的形变是地表真实三维形变在雷达视线方向(line-of-sight, LOS)上的投影，目前多通过除以入射角的余弦值得到垂直向的形变。但是，这需要假设水平方向上无形变[5]，很多情况下并不满足，从而导致形变信号的错误解译。

为了突破DInSAR技术的这些限制，学者们提出了多时相InSAR技术(multi-temporal InSAR, MT-InSAR)来进行高精度的形变监测，如永久散射体(persistent scatterer InSAR, PS-InSAR)、小基线集(small baseline subsets InSAR, SBAS-InSAR)和分布式散射体(distributed scatterer InSAR, DS-InSAR)等。此外，为了弥补D-InSAR或MT-InSAR技术只能获取LOS形变的缺陷，学者又提出多孔径InSAR(multi-aperture InSAR, MAI)技术来获取方位向(即卫星飞行方向)上的形变信息。以上这些InSAR技术已广泛应用于各类形变监测中，如城市地面沉降、矿区沉降、地震及板块运动、火山喷发、基础设施变形、冰川漂移、冻土形变、滑坡等。

本文旨在通过分析以上这些InSAR技术及其在各应用领域的研究现状，使读者能够对InSAR变形监测有更为系统性的理解。本文章节安排如下：第1节介绍InSAR技术与现状，包括雷达干涉测量原理和SAR卫星现状；第2节介绍目前主流的几种InSAR变形监测方法及进展；第3节介绍InSAR变形监测在重点应用领域的研究现状；第4节则对InSAR变形监测的挑战进行分析。

1 InSAR技术与现状1.1 雷达干涉测量原理

合成孔径雷达干涉测量技术(InSAR)集合成孔径雷达技术与干涉测量技术于一体。合成孔径雷达(SAR)是一种主动式微波遥感，用来记录地物的散射强度信息及相位信息，前者反映了地表属性(含水量、粗糙度、地物类型等)，后者则蕴含了传感器与目标物之间的距离信息。干涉的基本原理是同一区域两次或多次过境的SAR影像的复共轭相乘，来提取地物目标的地形或者形变信息。雷达干涉的模式分为：沿轨道干涉法、交叉轨道干涉法、重复轨道干涉法。其中利用重轨干涉最为常用的方式，以此为例，所得到的干涉相位表达式为[4]  (1)

式中，φflat为平地相位；φtopo为地形相位，该相位可以用来恢复地形信息；φdef为地表形变引起的相位；φatmo为大气延迟相位；φnoise为观测噪声引起的相位。

 将平地相位、地形相位、噪声相位、大气相位去除，即可得到地表形变相位，这样的雷达干涉技术被称为雷达差分干涉技术(D-InSAR)，文献[2]首次论证该项技术可用于探测厘米级的地表形变。其几何关系如图 1所示。

图 1 差分干涉测量几何关系Fig. 1 Geometric relation of differential interferometry

从图3中可以导出以下关系

  (2)

式中，Δr为所要求解的地面目标点的LOS向形变信息；λ为雷达波长。

 1.2SAR卫星的发展现状

 1978年，美国发射了L波段SEASAT卫星，在1981、1984年又相继发射了SIR-A和SIR-B卫星，为最初的研究提供了少量的星载数据。随后其他国家也相继研究并研发了不同波段的中低分辨率的SAR卫星。1991年7月、1995年4月欧洲空间局(ESA)相继发射了C波段ERS-1和ERS-2卫星，采用相同的轨道参数，组成串行飞行模式，为对地观测提供了大量的数据，极大地推动了InSAR技术的发展。1992年日本空间局发射了L波段JERS-1卫星，长波段的穿透力强，使其在InSAR监测中具有很强的优势，但JERS-1卫星并不是专门针对干涉测量而设计的，轨道误差较为严重，影响了它的适用性。1995年，加拿大发射了C波段RADARSAR-1卫星，获取了覆盖全球的大量数据。21世纪以后，SAR卫星系统的硬件得到更大的改进，专门针对干涉测量的卫星陆续发射成功。2000年美国实施的SRTM计划，利用航天飞机获取了全球80%的30m和90m分辨率的DEM数据，对差分干涉测量提供了可靠的外部DEM数据。欧空局2002年3月发射的ERS系列后继卫星ENVISAT不仅提供了SAR数据，还提供了MERIS水汽数据，对InSAR数据中大气延迟误差的改正提供了极大的帮助。2006年1月日本空间局发射了L波段的ALOS-1卫星，因其长波长穿透性强、抗失相干的特性在地震、矿山等方面的监测得到了广泛的应用。ERS-1/2、ENVISAT、ALOS-1、RADARSAT-1获取丰富的中低分辨率SAR数据成为21世纪早期InSAR技术发展的主要数据源，大部分的研究成果也是基于这些数据得到。

在上述的中低分辨率的SAR卫星陆续停止服务后，新一代高分辨率卫星相继发射，而且这些SAR卫星系统也从单一极化、模式、波段、固定入射角向多极化、多模式、多波段、可变入射角发展。由此，SAR卫星系统迎来了一个全新的时代。2007年德国航空航天局成功发射了X波段TerraSAR-X卫星，随后在2010年发射了其姊妹星TanDEM-X，构成分布式协同工作模式，可以提供全球高精度的数据及数字高程模型，分辨率最高可达0.25m。意大利空间局设计的X波段COSMO-SkyMed星座包括4颗星，分别于2007年6月、2007年12月、2008年10月、2010年11月成功发射，分辨率高达1m，重访周期短，在风险预警、灾害管理等方面发挥着重要的作用。2007年12月加拿大发射的C波段RADARSAT-2卫星继承了RADARSAT-1卫星的工作模式，还增加了多极化成像等，增强了其在地形监测领域的观测能力。2014年4月欧空局成功发射了首颗环境监测卫星C波段Sentinel-1A卫星，并于2016年4月发射了Sentinel-1B卫星，对全球大范围的地质、环境灾害的监测提供了丰富的数据。2014年5月日本成功发射的L波段ALOS-2卫星，空间分辨率最高达到1m，其在大型地质灾害(如地震)监测领域具有独一无二的优势。与此同时，我国也在积极的研发SAR卫星，并于2016年8月10日在太原成功发射了高分三号(GF-3)卫星，实现了我国卫星SAR影像干涉测量零的突破。

现有的SAR卫星系统具体参数如表1所示。在未来，SAR卫星系统的设计将会针对不同的研究领域而进行规划，例如，欧空局规划在2020年发射P波段BIOMASS卫星是为森林树高测绘和生物量反演而设计；德国航天局预计在2022年发射的L波段TanDEM-L卫星，主要针对全球陆表动态变化监测而设计。

表1现有SAR卫星系统具体参数Tab.1SpecificparametersofcurrentSARsatellites

2InSAR变形监测方法与进展

2.1D-InSAR方法

D-InSAR是在传统的InSAR技术上发展起来的方法，主要通过引入外部DEM去除InSAR获取的干涉图中的地形相位，进而得到差分干涉图[2]。而这差分干涉相位中除了地形相位外，还存在平地效应、形变信号、大气以及噪声成分。由于该差分干涉相位仍然是缠绕的，因此需对剩余的相位成分进行整周模糊度求解。具体地，通过精确的轨道信息来去除平地效应，利用一定的滤波来降噪，通过外部数据或模型来去除大气相位[6]，通过相应的相位解缠方法来获取整周模糊度[7]。自Gabriel在1989年首次利用D-InSAR成功获取了农田区的形变信息后，学者们相继对DInSAR技术做了一系列的改进，如滤波算法的提出和改进[8-9]、轨道相位的去除[10]、相干性的稳健性估计、电离层的去除[12]等。同时，D-InSAR也从最初的地震、火山等大范围地表形变监测慢慢发展到局部形变监测，如结合角反射器的矿区、滑坡、油气田以及城市地表形变。虽然D-InSAR已取得了一定的成果，但因监测精度仍受时空失相干、大气等因素的影响，因此基于D-InSAR发展起来的MT-InSAR是现今InSAR技术监测变形的研究热点。

2.2PS-InSAR方法

文献[13—14]首次提出了PS-InSAR方法，其基本思想是：第一，利用覆盖同一研究区的多景单视SAR影像，选取其中一景SAR影像作为主影像，其余SAR影像分别与主影像配准，依据时间序列上的幅度和(或)相位信息的稳定性选取永久散射体(PersistentScatterer，PS)目标；第二，经过干涉和去地形处理，得到基于永久散射体目标的差分干涉相位，并对相邻的永久散射体目标的差分干涉相位进行再次差分；第三，根据两次差分后的干涉相位中各个相位成分的不同特性，采用构建形变相位模型和时空滤波或方式估计形变和地形残余信息。

PS-InSAR技术不是针对SAR影像中的所有像元进行数据处理，而是选取在时间上散射特性相对稳定、回波信号较强的PS点作为观测对象。这些PS点通常包括人工建筑物、灯塔、裸露的岩石以及人工布设的角反射器等[14]。PS点的准确选取可以确保即便在干涉对的时间或空间基线很长的条件下(甚至达到临界基线)，PS点依然呈现出较好的相干性和稳定性。常用PS点选取方法包括振幅离差阈值法[14]、相干系数法[16]、相位分析法[17]，以及这些方法的组合等。而形变和地形残差解算通常采用解空间搜索法[14]、LAMBDA方法[18]和StamPS中的三维解缠法[19]等。

目前，PS-InSAR技术已在多个城市的高分辨率地面沉降监测中得到广泛应用，特别是城市重点基础设施的高分辨率形变监测。通过对比同期的水准和GPS测量数据，证实了PS-InSAR技术具有较高的可靠性(mm级的精度)[20]，因此对掌握城市地面沉降动态规律以及分析地面沉降的成因等具有重要意义。

然而PS-InSAR技术也存在自身的缺陷。首先，其通常要求覆盖同一区域的SAR影像数目较多(通常要求大于25景)，便于保证模型解算的可靠性。其次，PS-InSAR技术由于是基于大量PS点的迭代回归或网平差计算，运算效率不高，因此不适合大范围地区高分辨率的地面沉降监测。

2.3SBAS-InSAR方法

SBAS-InSAR是一种基于多主影像的InSAR时间序列方法，只利用时空基线较短的干涉对来提取地表形变信息，文献[21]提出并利用该技术进行意大利南部CampiFlegrei火山口和Naples市的时间序列形变估计，证明了该方法具有能够准确估计形变的能力。SBAS-InSAR技术克服了PS-InSAR因选取一幅影像作为公共主影像而引起的部分干涉图相干性较差的不足，同时降低了对SAR数据的需求量，运算效率较高。

SBAS-InSAR技术是以多主影像的干涉对为基础，基于高相干点恢复研究区域的时间序列形变信息，其原理如下：首先对覆盖某个地区的不同时间段的多景SAR影像计算时间空间基线，选择恰当的时空基线阈值选取干涉对；然后对选取的干涉对进行差分干涉处理并进行相位解缠；最后根据自由组合的干涉图形成子集的情况，对所有干涉图组成相位方程采用最小二乘法或者SVD方法进行形变参数的估计。在实际处理中会采用时空滤波的方法去除大气延迟的影像分离出非线性形变，估算的低频形变与此非线性形变的总和即为整个研究区的形变信息。

针对SBAS-InSAR技术，各国学者展开了许多研究，在高相干点选取方面主要有相位稳定性选点法[22]、振幅离差指数选点法[13]、空间相干性[23]等方法；在形变观测的数学模型选取方面，主要有线性模型[23]、二阶多项式形变模型[21]、周期性形变项[24]等；在参数估算方面主要有最小二乘法[23]、SVD法[21]、基于L1范数的解算方法[26]等。

为了提高高质量相干点的密度、减弱误差的影像、改善形变参数估计结果的精度，国内外学者提出了一些融合PS-InSAR和SBAS-InSAR优点的算法，如文献[23]提出了基于最小二乘网络的平差法，利用少量的SAR数据获取干涉图中高相干点的地表形变信息。文献[27]应用PS-InSAR和SBAS-InSAR融合技术于Eyjafjallajokull火山的形变监测中，获取了该区域的时间序列形变结果。文献[28]提出了时域相干点法(TCP-InSAR)，对形变参数、轨道误差同时进行估计，提高了监测精度。

2.4DS-InSAR方法

自2011年文献[29]发布第二代永久散射体技术SqueeSAR以来，时序InSAR领域的研究热点逐渐转向对分布式目标(distributedscatterer,DS)的探索。与PS目标的物理属性不同，DS是指在雷达分辨率单元内没有任何散射体的后向散射占据统治地位的点目标[30]。事实上，利用DS-InSAR进行形变监测的概念早在SBAS技术和StamPS技术中已经建立[21,27]，当时学界更多的是以相干、非相干目标加以区分，而弱化了PS与DS之间的物理界限。因此，从数据处理的角度来说，以SBAS为代表的小基线技术与SqueeSAR技术及其变种有本质不同，但它们又同时强调DS目标的信噪比，因而均属于相位滤波。

SqueeSAR的技术要点是：①通过同质点选取算法增强时序InSAR协方差矩阵的估计精度，并同时辅助PS与DS目标的分离；②通过相位优化算法从协方差矩阵中恢复时序SAR影像的相位。在第1个步骤中，其前提条件是相同SAR影像质地的像素具有相同相位中心，因此在时序统计推断的框架下，选取具有相同SAR统计分布的像素参与平均不仅可以提升相位信噪比，还能保留图像的空间分辨率[32]。相比之下，SBAS多采用常规多视处理或空间(自适应)滤波，是一种以牺牲空间分辨率为代价换取相位质量提升的方法[29]。在复杂形变特征环境下，这类方法极易损失细节，并伴随形变非形变区域的误读。在第2个步骤中，SqueeSAR在样本协方差矩阵服从复Wishart分布的基础之上，采用极大似然估计方法得到优化后的时序相位。由于似然估计量无解析表达式，需采用非线性优化“挤压”待估参数，便成为“Squeeze”一词命名的由来。值得注意的是，SqueeSAR一方面运用了所有干涉对的信息，另一方面却并未从滤波后的干涉相位中直接提取时序相位，这与大多数时序InSAR技术采用相位三角关系从滤波后的干涉相位中直接解算时序相位有本质区别[29]。在优化DS之后，与PS目标一起融入传统PS-InSAR数据处理框架就可以获得精度更高、空间分辨率增强的时序形变产品。

SqueeSAR技术也有其缺陷。在同质选点算法中，使用的KS假设检验功效低，特别在小样本条件下易使选取的同质集合中包含许多异质点，而且计算效率低下[32]。因此，近年来国内外学者纷纷提出了计算更为高效或估计精度更高的统计推断方法，具有代表性的是AD检验[33]、似然比检验[32]、自适应检验和置信区间估计等[34]。在相位优化算法中，计算效率是该方法的主要诟病之处，使得大范围的精细监测几乎无法施展，备择的方案是文献[35]提出的协方差矩阵奇异值分解方法。在优化可靠性方面，可采用抗差性更强的估计量，如M估计量等[36]。

上述改进算法与SqueeSAR技术一起构成了当今DS-InSAR的雏形，旨在提高计算效率、增加观测密度和改善估计精度。广义的，除基于PS目标之外的时序InSAR技术都或多或少地应用了DS的概念，而联合PS-InSAR和DS-InSAR或许是提高最终产品精度的最佳方案[27,29]。