时序InSAR解译2017～2020年北京地面沉降时空变化

张双成 许 强 罗 勇 雷坤超,4 牛玉芬 庞校光

1 长安大学地质工程与测绘学院，西安市雁塔路126号，7100542 地理信息工程国家重点实验室，西安市雁塔中路1号，7100543 北京市水文地质工程地质大队，北京市西四环北路123号，1001954 中国科学院地质与地球物理研究所，北京市北土城西路19号，1000295 河北工程大学矿业与测绘工程学院，河北省邯郸市太极路19号，056038

地面沉降是我国平原地区主要的区域性环境地质灾害，可能对人民的生命和财产安全造成危害[1]。北京是中国地面沉降最严重的地区之一[2]，北京的沉降现象最早出现在20世纪50年代，并在近几十年来迅速扩散[3]。截至2009年底，北京市最大年沉降量已达137.51 mm，最大累积沉降量为1 163 mm[4-5]。因此，监测和分析北京市地面沉降活动的时空变化，对于防灾减灾和城市化的可持续发展具有重要意义。

InSAR具有大范围、高精度、全天候、全天时的特点，已逐渐成为监测城市地面沉降的主要技术手段[6]。Hu等[7]利用短基线集(small baseline subset，SBAS)方法获取2003～2010年北京地面沉降时空演变特征，结果表明，研究区内存在3个大型沉降漏斗；周吕等[8]运用SBAS方法获取2007～2010年北京地区的地表沉降分布，结果表明，不均匀沉降较为明显，各个沉降漏斗逐渐连成一片，且有向东发展的趋势；Zhang等[9]基于一种改进的多时相InSAR方法获取1992～2014年京津冀地区3个时间跨度的沉降速率，结果表明，北京市1992～2000年间只有部分小规模沉降点，2003～2010年地面沉降范围急剧扩大，2012～2014年间地面沉降持续扩散。南水北调工程于2014-12开始向北京输送水源，截至2020-09，平原区地下水埋深22.49 m，较2015年底的25.75 m上升了3.26 m，地下水水位连续5 a持续回升。南水北调工程为北京市提供新的水源，改变了2015年以来北京市地面沉降演变的格局。因此，不断更新北京地面沉降速率以预测危急情况，并在必要时采取缓解措施是很有必要的。

选取覆盖北京地区的85景C波段的Sentinel-1数据集，采用SBAS技术获取研究区域2017～2020年的沉降分布情况和沉降速率场，并通过对比分析2018年和2019年的年形变速率来探究北京地面沉降的时空变化规律，为城区的沉降监测和预警提供理论依据。

1 研究区域及实验数据

1.1 研究区域

北京市位于华北平原的西北边缘，中心位置为39°54′N、116°23′E。其西部、北部均为山区，东南部为向渤海缓慢倾斜的平原。山区面积为10 072 km2(占总面积的61.4%)，平原面积为6 388 km2(占38.6%)；山区平均海拔为1 000～1 500 m，平原地区平均海拔为20～60 m。北京属暖温带半湿润季风型大陆性气候，夏季炎热多雨，冬季寒冷干燥；年平均气温为10～12 ℃，年平均降雨量约为600 mm。

1.2 实验数据

选取2017-06～2020-06覆盖实验区的Sentinel-1A序列数据集，共包含85景升轨影像。采用覆盖研究区的30 m分辨率SRTM(shuttle radar topography mission)DEM作为参考高程数据，用于辅助去除地形相位。同时获取研究区域内8个GPS站点数据，用于InSAR结果的对比验证。

2 数据处理方法

SBAS方法的原理是通过设定一定的空间和时间阈值，将所有的SAR数据组成若干个小集合，集合内基线较小，集合间基线较大。最后通过集合内的最小二乘求解和集合间的奇异值分解方法，得到整个时间序列地表形变信息的联合求解结果[10-11]。

假设获取N+1幅SAR影像，其获取的时间序列为(t0,…,tN)，依据干涉基线组合可生成M幅干涉图，则在像元(x,r)处由tA、tB两个时刻SAR影像生成的差分干涉相位可表示为：

δφj(x,r)=φ(tB,x,r)-φ(tA,x,r)≈

φdef,j(x,r)+φtopo,j(x,r)+φatm,j(x,r)+

φnoise,j(x,r)，j=1,…,M

(1)

式中，(x,r)为方位向和距离向坐标，φ(tB,x,r)和φ(tA,x,r)分别为生成干涉相位的影像相位，φdef,j(x,r)为tA时刻至tB时刻间卫星视线向的形变相位，φtopo,j(x,r)为参考DEM不精确引起的地形相位误差，φatm,j(x,r)为大气相位误差，φnoise,j(x,r)为噪声相位。

由式(1)可以得到一个方程组，包含N个未知数的M个方程：

δφ=Aφ

(2)

式中，δφ为M幅干涉相位构成的矩阵，φ为N幅SAR影像上的待求形变相位构成的矩阵；A为M×N矩阵，其每一行对应一个干涉对。当小基线子集个数L=1时，A为列满秩矩阵，其最小二乘解为：

(3)

当小基线子集个数L>1时，式(2)是秩亏的，秩亏数为N-L+1，可以对A进行奇异值分解，求出形变相位φ最小范数意义下的最小二乘解[12]。

3 时序InSAR解译地面沉降及精度评定

3.1 地面沉降空间分布

利用SBAS方法对覆盖北京市的Sentinel-1影像数据进行地表形变特征信息的提取，获取北京平原地区沿卫星雷达视线向上的地表年平均形变速率，如图1所示(负值表示地表正在远离卫星，正值表示地表正在靠近卫星)。在解缠相位时，需要选择一个像素作为解缠的起始点(即参考点)，该点位于高相干区域的局部最高处，变形速度场中各监测点的平均沉降速度均与该点有关，该点的位置如图1所示。

图1 北京2017-06～2020-06年平均形变速率Fig.1 Annual mean deformation velocity of Beijing from June 2017 to June 2020

由图1可见，研究区域内不均匀沉降较为突出。2017-06～2020-06沉降主要发生在昌平(以下简称CP)、顺义(SY)、通州(TZ)、大兴(DX)以及北京和天津廊坊交界处(LF)，最大年形变速率在SY区域(为-111.3 mm/a)。CP沉降区主要的形变速率范围为-20～-60 mm/a，平均值为-37.8 mm/a，最大形变速率达-84 mm/a；SY沉降区主要的形变速率范围为-20～-70 mm/a，平均值为-43.1 mm/a，最大形变速率达-111.3 mm/a；TZ沉降区主要的形变速率范围为-25～-65 mm/a，平均值为-45.3 mm/a，最大形变速率达-89.7 mm/a；DX沉降区主要的形变速率范围为-20～-45 mm/a，平均值为-31.9 mm/a，最大形变速率达-63.6 mm/a；LF沉降区主要的形变速率范围为-20～-40 mm/a，平均值为-30.2 mm/a，最大形变速率达-59.5 mm/a。

研究区域2017-06～2020-06期间根据C波段的Sentinel数据获取的累积形变如图2所示。2017～2020年期间，TZ沉降区最大累积量约为255 mm，SY沉降区最大累积量约为287 mm，CP沉降区最大累积量约为241 mm，DX和LF沉降区相对较小，最大累积量分别约为157 mm和175 mm。

图2 北京2017-06～2020-06累积形变Fig.2 Cumulative deformation of Beijing from June 2017 to June 2020

3.2 地面沉降精度评定

对2017-06～2020-06视线向上的平均形变速率进行误差评估，通过计算速度线性拟合的偏差获取其对应的均方根误差RMSE(图3)。由图3可见，总体上RMSE不超过4 mm/a，平均值为1.4 mm/a，这对于InSAR时间序列结果来说精度较高。北京市中心的RMSE较小，小于1 mm/a，形变区和西北区域的RMSE为2～4 mm/a。推测造成此误差分布格局的原因是，北京市中心相干性较好；由于形变区沉降明显不均匀及西北区处于山区，相干性较差。

图3 平均形变速率的RMSEFig.3 RMSE of the mean velocity results

获取2017～2019年期间8个GPS站点的时序形变监测值，总共3期数据，各期起始监测时间分别为2017-08-24、2018-08-25和2019-07-23，各期结束时间分别为2017-10-23、2018-10-31和2019-09-03。为了将InSAR结果与GPS观测值在时间段上进行统一，采用2017-10-23～2019-09-03期间内的影像数据生成InSAR视线方向上的平均速率图。

提取各个GPS点对应的InSAR结果与视线方向上的GPS形变量进行对比，结果如表1和图4所示，最大绝对误差为4.45 mm/a，最小绝对误差为-5.48 mm/a，各站测量值与InSAR结果之间的RMSE为3.04 mm/a，表明两者之间具有良好的一致性。GPS与InSAR结果的相关系数为0.96，说明InSAR获取的视线方向的平均速率的精度较高。

表1 InSAR结果与GPS观测值对比

图4 InSAR结果与GPS观测值对比Fig.4 Comparison between time series ofInSAR and GPS measurements

3.3 地面沉降时空变化

为了探究研究区地面沉降的时空变化规律，分别利用2018年和2019年期间影像数据组成干涉对获取研究区2018年和2019年视线向上的形变速率，如图5所示(负值表示地面位移方向为远离卫星，正值表示地面位移方向为朝向卫星)。由图可见，北京地区2018年与2019年的沉降区空间分布具有较高的一致性，但从沉降速率及沉降范围来看，北京地区出现了减速沉降的现象。例如，2018年北京地区最大沉降速率为137 mm/a，而2019年最大沉降速率为110.1 mm/a，下降了26.9 mm/a。2018年和2019年形变速率图栅格数量分别为9 090 995和9 177 465个，之间的差量相比于总数可以忽略不计。对整幅速率图进行统计分析可知，2018年平均沉降速率为10.2 mm/a，而2019年平均沉降速率为5.6 mm/a，下降了4.6 mm/a。由此可知，北京市整体沉降减缓。

图5 北京视线向平均形变速率Fig.5 Mean deformation velocity in LOS direction of Beijing

2018年和2019年沉降速率变化如图6所示(栅格大小为46.60 m×69.64 m)，考虑到SBAS-InSAR结果的不确定性，认定沉降速率大于10 mm/a的区域为沉降区。由图可见，在各个沉降速率范围内，2019年形变速率图包含的栅格数量都比2018年的少。在沉降速率为10～20 mm/a的范围内，减少的栅格数量最多，达到18万个，减少幅度为13.62%；其次是沉降速率为20～30 mm/a的范围，减少16万个；在沉降速率为60～70 mm/a、70～80 mm/a、80～90 mm/a和90～140 mm/a的范围内，共计减少15万个栅格，减少幅度较大，分别为74.72%、75.41%、81.96%和75.54%。由此可知，在各个沉降范围内沉降的面积都在减小。

图6 2018年和2019年沉降速率变化Fig.6 The deformation velocity change in 2018 and 2019

将2018年和2019年形变速率图进行差异化分析，获取这2 a速度变化幅度的统计数据，如图7所示(负值表示2018～2019年沉降减缓，正值表示2018～2019年沉降加速)。考虑到SBAS-InSAR结果的不确定性，将速度变化为-10～10 mm/a的栅格视为稳定点。图7中处于沉降减缓的区域占总面积的10.3%，而加速区域仅占0.6%，减缓的面积是加速面积的16倍。由此可知，沉降减缓的面积远大于沉降加快的面积。

图7 2018年和2019年形变速率的差异Fig.7 Difference indeformation velocity between 2018 and 2019

分析CP、SY、TZ、DX、LF五处沉降区(图5中沉降区对应的蓝色长方形区域)的沉降状况，并对局部沉降速度变化进行调查(图8)。由图可见， SY、TZ、DX、LF 四处沉降区2018～2019年的最大形变速率和平均形变速率都在下降，其中TZ区下降最为明显，最大沉降速率由113.0 mm/a下降到77.8 mm/a，平均沉降速率由50.7 mm/a下降到23.7 mm/a；CP区沉降仍在加速，最大沉降速率由88.4 mm/a提高到110.0 mm/a，平均沉降速率由37.5 mm/a提高到49.3 mm/a。由此可知，SY、TZ、DX、LF四处沉降区沉降均在减缓，而CP区沉降仍在加快，其中TZ区沉降减缓的幅度大于CP区沉降加快的幅度。

图8 5处沉降区形变速率统计Fig.8 Statistics of deformation velocity in five subsidence areas

4 结 语

基于2017-06～2020-06共85景Sentinel-1影像数据，利用SBAS-InSAR技术对北京地表形变进行监测，分析2017～2020年北京地面沉降时空变化。得到的SBAS-InSAR结果与GPS观测值吻合良好，R2为0.96，RMSE为3.04 mm/a。结果表明：

1)北京不均匀沉降较为突出，2017-06～2020-06地表形变呈现5处沉降区，最大年沉降速率为-111.3 mm/a，最大累积量为287 mm。

2)综合2018年和2019年形变速率图可以看出，北京整体地面沉降情况在减缓，在各个沉降范围内的沉降速率都在减缓，且沉降减缓的面积远大于沉降加速的面积。

3)对局部沉降调查发现，北京市5处沉降区中除1处仍然在加速外，其他4处沉降速度均在减缓。

由于本次实验获取的SAR数据集时间较短，可能会对实验结果造成一定的影响。随着Sentinel-1卫星进入业务化运行，SAR数据的不断累积可为更好地揭示研究区域地面沉降时空演变规律提供丰富的研究资料。地下水是维持土壤应力平衡的重要因素，研究地表沉降与地下水变化的关系，可进一步全面解译南水进京后北京地面沉降的发展趋势。

致谢：感谢欧空局提供Sentinel-1雷达影像数据、美国航空航天局提供SRTM DEM地形数据以及北京市水文地质工程地质大队提供部分GPS监测数据。