地面沉降是较为普遍的地质灾害，具有形成缓慢、持续时间长、影响范围广、成因机制复杂和防治难度大等特点。地面沉降已在我国20多个省(区、市)发生，每年造成数百亿元的直接经济损失，是影响经济社会发展的最重要的灾害因素之一。京津冀地区所在的华北平原是世界上最大的沉降漏斗^[1]。对地面沉降进行精确、大范围、长时间监测是开展沉降控制工作的基础。以基岩标、分层标、 GNSS、水准点等为主要监测手段的地面沉降监测网络存在空间密度低、无法提供整个区域的沉降信息等问题，难以为防灾减灾提供可靠的决策依据。

始于1989年的差分合成孔径雷达干涉测量(differential interferometric SAR, D-InSAR)技术^[2]曾在地表形变监测中展现了巨大潜力，因而受到广泛关注。但是，传统的两轨或三轨D-InSAR技术受到空间失相干、时间失相干和大气非均匀等因素^[3-4]的限制，难以真正实用化。近10多年来得到迅速发展的时间序列InSAR技术是对传统D-InSAR技术的革新。它对时间序列干涉图进行分析，处理的对象不是整幅影像的全部像元, 而是其中具有稳定散射特性的高相干像元子集(下称“高相干点”)。时间序列InSAR技术总体上包括两类：以永久散射体干涉(permanent scatterer interferometry，PS-InSAR)^[5-6]或者Persistent Scatterer Interferometry^[7]为代表的单一主影像时间序列InSAR技术和以小基线集技术(small baseline subset interferometry, SBAS InSAR)^[8]为代表的多主影像时间序列InSAR技术。时间序列InSAR技术已在火山、地震、滑坡、地面沉降等领域得到大量应用。

鉴于京津冀地区地面沉降的严峻局面，2014年1月国家测绘地理信息局在开展第一次全国地理国情普查工作的同时，启动了京津冀地区重要地理国情监测项目，地面沉降是其中的4项监测内容之一。目前，一些研究者已经开展过北京^[9-10]、天津^[11-13]、廊坊^[14]、沧州^[15-16]等地的地面沉降InSAR监测，但都是针对局部地区、较短时段的试验性研究。本研究利用1992年以来多颗卫星获取的SAR影像对京津冀大区域长时间跨度的地面沉降进行监测，面临着时间序列SAR数据集小、大区域点目标密度不均一、不同轨道沉降监测结果镶嵌等难题。因此，数据处理方法不能照搬现有方法。本文提出一种改进的时间序列InSAR技术，即多主影像相干目标小基线技术，对京津冀地面沉降进行遥感监测，并对地面沉降监测结果进行了精度评价，最后分析了京津两地22年间地面沉降时空变化的主要特点。

1 SAR影像数据源

根据不同合成孔径雷达卫星在京津冀地区的数据存档情况，分别选用ERS-1/2、ENVISAT、RADARSAT-2 SAR影像进行1992-2000年、2003-2010年、2012-2014年的地面沉降监测。所有影像均为单视复数据(SLC)格式。不同传感器影像的覆盖范围如图 1所示。

ERS/ENVISAT SAR影像为条带扫描模式，距离向和方位向像元大小为7.90 m和3.98 m，单景影像的覆盖区域约为100 km×100 km。由于条件所限，只获取了7个Frame的ERS/ENVISAT SAR影像。轨道(Track)号分别为490、218、447、175。这7景ERS/ENVISAT SAR影像的覆盖面积为5.03万km²，完整覆盖了北京市和天津市的平原地区，同时覆盖了河北省唐山市、廊坊市、保定市的大部分平原地区。

RADARSAT-2 SAR影像为宽模式，距离向和方位向像元大小分别为11.83 m和5.21 m，单景影像的覆盖面积约为150 km×150 km。本次监测使用了8景RADARSAT-2 SAR宽模式影像，其覆盖面积为9.36万km²，包含北京市和天津市的全部平原地区以及北京以南的河北省平原地区。为节省数据成本，每个RADARSAT-2单景覆盖区只获取了19期时间序列影像。京津冀全区域共获取了152景RADARSAT-2影像。

ERS-1/2、ENVISAT及RADARSAT-2时间序列SAR影像集的大小及获取的起始时间列于表 1和表 2中。可以看出，ERS-1/2 SAR影像数据集普遍较小，最小的只有11期影像。同时，ERS-1/2卫星在其生命周期有不同的任务阶段(mission phases)，导致ERS-1/2 SAR序列影像的时间和空间基线很不均匀，造成ERS-1/2 SAR干涉失相干较严重，因此形变提取的精度不如ENVISAT ASAR影像的精度高。后面的监测结果也将印证这一点。

2 时间序列InSAR处理方法

本文提出了一种改进的时间序列InSAR技术，即多主影像相干目标小基线干涉技术(MCTSB-InSAR)。先利用MCTSB-InSAR技术进行单景覆盖区数据处理，获取单景影像覆盖区的地面沉降(形变)信息，包括平均沉降速率和时间序列沉降量。大区域通常需要多景SAR影像才能覆盖，因此需要将不同景，特别是不同轨道的沉降提取结果进行拼接，以得到大区域覆盖完整且一致的沉降结果。MCTSB-InSAR技术能实现多景形变结果的平差及镶嵌处理。这是MCTSB-InSAR方法区别于其他方法的一个优点。除文献[17]外，多景形变结果的连接问题很少被讨论。

2.1 单景覆盖区的时间序列InSAR处理

MCTSB-InSAR技术中的单景覆盖区时间序列处理综合了PS-INSAR和SBAS InSAR方法的优点，克服了相应的缺陷。它不以单一SAR作为全部干涉像对的固定主影像，而是考虑满足小基线要求的时间序列SAR数据集任意组合的干涉对，因此不像PS-INSAR方法需要25景以上的SAR影像才能进行有效的分析。MCTSB-InSAR方法也不像SBAS方法需要先完成每个差分干涉图的相位解缠，然后再进行时间序列分析，而是通过考察空间邻近的两个高相干点的二次相位差分，同时完成形变参数解算和相位解缠。一般情况下，MCTSB-InSAR将形变分解为线性形变(主形变)加上非线性形变(次形变)。如果形变呈现强烈的非线性，则可将形变分解为高阶多项式形变加上残余形变^[18]。MCTSB-InSAR方法中的单景覆盖区时间序列处理主要包括以下5个方面。

2.1.1 小基线干涉组合

对SLC格式的SAR影像完成精确配准后，根据预先设定的关于时间基线和空间基线的阈值条件自由组合生成干涉图。这种方式在处理小数据集时尤为重要，因为可以使用较小数目的SAR影像得到较多的干涉图。MCTSB-InSAR方法要求这些干涉图组成一个连通网络。如果这些干涉图形成了两个或多个小基线干涉图集，则需要手动添加几个干涉组合将这些分离的小基线集连通起来。这也是与SBAS方法的一个不同之处。因为SBAS方法要求先对干涉图进行相位解缠，然后通过奇异值分解技术解决孤立的干涉图集之间相位的一致性问题。MCTSB-InSAR方法不要求事先对干涉图进行相位解缠，孤立的干涉图容易造成形变模型参数求解不唯一。

2.1.2 高相干点提取

高相干点指干涉相干性较高的像素，既包括具有优势散射体的点目标(permanent scatterer or persistent scatterer，PS)，也包括能在一段时间内保持较高相干性的分布式目标。对点目标和高相干分布式目标可以设置不同的参数及阈值条件分别提取，然后求其点集的并集作为最终的高相干点。提取PS点的参数有幅度离差^[6]和平均振幅。提取高相干分布式目标的参数有平均相干系数^[19]等。大区域监测中将点目标和高相干分布式目标都包括进来非常重要，因为城市区域点目标密集，但是非城区点目标很少。如果不考虑具有较高相干性的分布式目标，在非城区提取的高相干点密度不足，将难以建立高相干点有效的网络连接。

2.1.3 高相干点网络连接

通常的时间序列InSAR方法^[19-20]都是采用整体的Delaunay三角网将高相干点连接在一起。如果研究区仅限于城市区域，由于高相干点密度很高，整体Delaunay三角网连接比较稳健。但是，如果研究区包含很大的非城区范围(例如整景影像)，由于非城区高相干点密度不高，或者高相干点质量不高，有些连接将会在后续处理中被剔除。这样，整体Delaunay三角网连接由于冗余度不够，在剔除不满足要求的连接边后极有可能分裂成许多不连通的子网络，给后面的积分操作带来困难。为了增加网络连接的冗余度，MCTSB-InSAR技术采用了局部Delaunay三角网连接方法^[21]，即将整景影像分割互相重叠的网格，在每个网格中采用Delaunay三角网连接。最终连接的冗余度依赖于网格大小和相邻网格的重叠度，网格越小，重叠度越大，整体冗余度越大。为了减弱连接边上大气相位的影响，要求连接边的长度不大于大气相关距离(根据SAR影像分辨率不同，通常设为1000~2000 m)。

2.1.4 线性形变求解

假设一共有M幅干涉图，去掉平地和去地形相位后，第k(1≤k≤M)个差分干涉图上的第i(1≤i≤N)个高相干点的相位可以用式(1)表示

式中，wrap{}表示相位缠绕算子；φ_i^k表示差分干涉相位；φ_{i, atmo}^k表示大气延迟所导致的相位贡献；φ_{i, non-linear}^k是非线性形变相位；φ_{i, noise}^k是噪音相位；Δh_i是DEM误差；t_i^k是干涉图的时间基线；v_i是线性形变速率；B_{i, ⊥}^k是垂直基线；系数a₁由该点的斜距、入射角以及雷达波长λ决定，a₂=4π/λ。

在完成高相干点的网络连接后，第k个干涉图上由相邻高相干点i, j构成的连接边的干涉相位差可以表示为

式中，Δh_{i, j}表示两高相干点间的高程误差之差；Δv_{i, j}表示两点间的线性形变速率之差；Δφ_{i, j, res}^k=Δφ_{i, j, non-linear}^k+Δφ_{i, j, atmo}^k+Δφ_{i, j, noise}^k，是两点上的大气贡献相位、非线性形变相位和噪声贡献相位的差分之和。由于大气延迟相位具有较强的空间相关性，非线性形变相位尽管表现为空间非相关性，但是相位贡献相对较小，因为是高相干点，单点的噪声相位很小，因此|Δφ_{i, j, res}^k＜π成立。如此，按式(3)定义的模型相干因子最大化即可求解未知数Δh_{i, j}、Δv_{i, j}

式中，J表示单位复数。求解未知数Δh_{i, j}、Δv_{i, j}的方法有二维周期图法^[5]或者空间搜索法^{[18, 21]}，二者是等价的。当对所有的边完成最大化求解后，选择γ_{i, j}≥0.7的边作为可靠的连接，那些没有被可靠边连接的高相干点将被剔除。最后，从某一参考点开始对连接边上解算出的速度和高程误差增量进行积分，将得到各相干点上形变速率和高程误差的绝对量。

2.1.5 非线性形变估计

从差分干涉相位中去除高相干点的线性形变相位和高程误差相位，得到高相干点的残余相位，包括大气影响相位、非线性形变相位以及噪声相位。大气影响相位在时间序列上是不相关的，为高频信号；在空间分布是相关的，为低频信号；非线性形变相位在时间序列上是低频信号，在空间上是不相关的，为高频信号；噪声相位则是时间和空间都不相关的随机高频信号。利用残余相位3个分量的不同的时空频谱特征，可以将三者分离出来^[5]。将线性形变和非线性形变相位联合起来就可得到每个影像获取时刻的累积地表形变。假定京津冀地表形变主要体现为地面沉降，基于雷达波入射角关系，将获取的雷达视线向形变量转换为垂直方向上的地面沉降值^[22-23]。

2.2 多景覆盖区的地面沉降结果平差与镶嵌

单景影像覆盖区的时间序列SAR影像的起始时间通常不一致，累积沉降量自然不同。因此累积沉降量难以镶嵌。将多景覆盖区的沉降速率结果进行平差，可以获得大区域的平均沉降速率镶嵌图。

同名高相干点的识别是平差及镶嵌的基础。虽然利用影像匹配技术能通过相邻SAR振幅影像的匹配完成同名高相干点的自动识别，但耗时太长且易出错。文献[18]采用将重叠区的不同轨SAR影像转化到一个统一的坐标系统中，并将其中某轨SAR影像坐标作为统一的坐标系，在重叠区高差不大的情况下，可以通过多项式将其他轨SAR影像坐标变换到这个统一坐标系。显然，如果重叠区高差大，这种做法难以确认同名高相干点。事实上，最好的统一坐标系是地理编码后得到的大地坐标系(此时，高程导致的不同轨像元定位差异已经得到改正)，利用地理编码后同名高相干点有相同大地坐标这一原理即可有效解决这一问题。必须注意的是所有SAR影像的地理编码应该使用同种类型的DEM，如SRTM DEM，这样能保证相邻景SAR影像地理编码后有相同的平面精度。

理论上，用一个偏移值可足够说明不同轨所得形变结果的差异^[18]。文献[18]采用了距离和方位向的一次多项式来说明这个差异。这种做法极易将形变场中可能存在的线性趋势分量吸收，当成是不同轨上由于轨道误差、大气误差等造成的InSAR形变差异，反而造成平差后更大的不一致。另外，文献[17]讨论的是雷达视线向形变结果的连接，真实的形变投影到不同轨的雷达视线向后可能存在沿距离和方位向线性变化的趋势性差异。本文已经将视线向形变转化为垂向沉降，消除了不同视角投影的影响，因此在进行多轨沉降结果的平差时只考虑了偏移(即零次多项式)。

平差方程由两类构成。第1类方程是根据同名高相干点的形变速率一致的准则建立的，形如式(4)

式中，Vi(x, y)、Vj(x, y)分别代表相邻单景覆盖区(第i景和第j景)同名高相干点(x, y)上的形变速率值；Δi、Δj为待求的第i景和第j景的形变速率偏移值。第1类方程进一步可写成

第2类方程是将水准(或GPS)测量结果当成控制点构成的，如式(6)

式中，GCP_V(x, y)表示在点(x, y)处经由水准测量得到的形变速率值。将如式(5)和式(6)的所有方程联立在一起，即可解得所有单景覆盖区InSAR形变速率的偏移值。在有控制点的情况下，平差求得的偏移值同时也是每个单景覆盖区形变结果的绝对定标值，这样同时完成了每景覆盖区形变速率的平差与绝对定标。如果没有控制点，则没有第2类方程。通过平差将会得到以某一个单景覆盖区形变结果为参考的相对偏移值。最后，将经偏移值改正后的同名点的形变速率取算术平均，即可实现整个大区域的形变速率的镶嵌。

3 京津冀地区22年地面沉降监测结果分析

根据上述方法，提取了京津冀地区1992-2000年、2003-2010年、2012-2014年的地面沉降速率和累积沉降量等信息，并得到了3个时段地面沉降速率的大区域镶嵌结果，如图 2-图 4所示。在大区域沉降结果平差之前，已经对沉降结果进行了地理编码。ERS-1/2和ENVISAT沉降结果地理编码的输出像元大小为40 m，RADARSAT-2沉降结果地理编码的输出像元大小为30 m。图 2-图 4中负值表示地面下沉，正值表示地面抬升，负值越小表示沉降越严重。从图 2可以看到ERS-1/2影像得到的大区域镶嵌图上存在一些不一致，特别是T175/F2799与T175/F2817之间，这主要是由上述的ERS-1/2 SAR数据集小、时空基线差异大等因素造成的。相比之下，图 3和图 4的大区域沉降速率镶嵌结果则有很好的一致性。这说明本文的大区域形变结果的平差与镶嵌方法是有效的。

为了标定和验证InSAR监测结果，从北京市和天津市测绘部门获取了相近时段的水准数据，包括：天津市1993-2000年184个观测数据，用于ERS-1/2结果定标(即平差控制)和精度验证；北京市2005-2010年74个观测数据、天津市2003-2010年175个观测数据，用于ENVISAT ASAR结果定标和精度验证；北京市2012-2014年26个观测数据、天津市2011-2013年185个观测数据，用于RADARSAT-2结果定标和精度验证。其中1/3数据用于平差中的控制，2/3数据用于精度验证。对RADARSAT-2结果进行定标和精度验证的水准点分布如图 5所示，其中蓝色点为平差控制点，共70个，红色点为精度验证点，共141个。表 3为各时段InSAR沉降速率的精度验证结果。可以看出，除1992-2000年时段的沉降速率精度略差外，ENVISAT和RADARSAT-2获取的沉降速率精度约在5 mm/a。

从图 4中可以清晰地看到整个京津冀三地近3年的沉降中心分布。京津冀地区地面沉降非常严重，三地沉降漏斗有连接成片的趋势。

由于3个时段的监测区基本完整地包括了京津两市的平原地区，笔者对京津两市1992-2014年的地面沉降(图 6和图 7)的时空变化特点进行了进一步分析。在分析沉降区面积时，首先采用克里金(Kriging)插值方法将高相干点上的沉降速率值内插至每个像元，然后进行面积统计。

图 6所示，北京平原地区地面沉降1992-2000年处于初步发展阶段，零星分布一些较小的沉降中心，最大沉降速率为48 mm/a，位于门头沟区龙泉镇。经查实，该沉降主要由龙泉镇煤矿地下采煤引起(该煤矿2000年被关闭)。2003-2010年北京地区地面沉降快速发展，在东部朝阳-通州形成了沉降带，在北部海淀、昌平、顺义等地形成了5个散落分布的沉降中心，沉降速率最大达143 mm/a，位于朝阳区东坝乡(门头沟煤矿的沉降已经停止)；沉降速率超过50 mm/a的面积为265.41 km²。2012-2014年东部沉降带继续发展，沉降范围和速率均继续增大；北部散落分布的沉降中心继续扩展，连接成片，形成海淀-昌平-顺义沉降带；最大沉降速率达152 mm/a，仍然位于朝阳区东坝乡；沉降速率超过50 mm/a的面积达到433.25 km²。整体来看，在过去3个时段，北京市地面沉降呈不断加重的趋势。

从图 7可以看出，天津市1992-2000年地面沉降已有一定规模，最大沉降速率为81 mm/a，位于汉沽区，沉降速率超过50 mm/a的面积达到85.73 km²。2003-2010年，天津市地面沉降的速率和范围都在迅速发展，在西北部、南部、东南部形成了3个大的沉降带；最大沉降速率为135 mm/a，位于静海区；沉降速率超过50 mm/a的面积增长到2749.1 km²。2012-2014年，整体上呈现“北增南减”的格局，北部武清、北辰和汉沽等地区沉降速率增加，南部地区沉降均有减缓；最大沉降速率达153 mm/a，位于武清区南部的王庆坨镇；沉降速率超过50 mm/a的面积缩小到1 117.55 km²。总体上，天津市地面沉降在2010年以后有所减缓。

从笔者收集的北京地区地下水流场监测数据、基础水文地质数据等资料来看，北京地区地面沉降分布范围与地下水降落漏斗分布在空间位置上基本一致，说明长期超量开采地下水是导致北京地面沉降发生的主要原因。

4 结论

本文重点介绍了京津冀大区域1992-2014年地面沉降InSAR监测的技术方法和结果。该监测结果首次提供了京津冀地区过去22年间地面沉降的完整、可靠信息，对于京津冀地区地面沉降灾害防治、规划制定、水资源合理使用、经济布局调整等具有重要的参考价值。研究成果同时表明，本文提出的大区域时间序列InSAR形变监测技术MCTSB-InSAR是有效的，可在未来全国大范围的地面沉降监测中推广应用。

致谢: 特别感谢欧空局提供的ERS-1/2 SAR及ENVISAT ASAR数据，加拿大MDA公司提供的RADARSAT-2数据，北京市测绘设计研究院、北京市水文地质工程地质大队、天津市测绘院提供的水准测量数据。