利用SBAS-InSAR技术监测北京地面沉降

祝 杰 刘洋洋 邵银星 陶志刚

（中国北京100045中国地震台网中心）

0 引言

地面沉降是在自然或人为因素作用下，由于地下松散地层固结压缩导致地壳表面标高缓慢下降的一种局部地质现象，具有形成时间长、影响范围广、防治难度大和过程不可逆等特点（何庆成等，2006）。研究表明，我国已形成长江三角洲、华北平原、汾渭盆地呈三足鼎立的地面严重沉降区，近百个城市正在遭受不同程度的地面沉降。北京地区的地面沉降经历了形成阶段（20世纪50—70年代）、发展阶段（20世纪70—80年代）、扩展阶段（20世纪80—90年代）和快速发展阶段（2000年至今）（周朝栋等，2017），沉降演化过程呈现范围不断扩大、变形逐渐加速的趋势。

对地面沉降进行大范围、长时间、高精度监测是开展沉降控制工作的基础。基岩标、分层标、GNSS（Global Navigation Satellite System）、水准点等监测方式存在空间密度低、无法提供广域形变信息等问题，难以为防灾减灾提供可靠的决策依据（张永红等，2016）。合成孔径雷达干涉测量（Interferometric Synthetic Aperture Radar，InSAR）技术具有明显的大范围、高精度、高时空分辨率、全天候监测等特点，弥补了现有常规沉陷监测手段的不足，已在地震同震、震间和震后地壳形变监测、地质灾害隐患早期识别、区域地面沉降监测等领域（王超等，2002；常占强等，2008；张景发等，2008；许强等，2019）成功应用。在北京沉降研究中，李永生等（2013）利用2个相邻轨道的ASAR数据提取北京地区沉降信息，分析了沉降空间分布特征，发现2003—2010年北京主要沉降区域快速下沉。张永红等（2016）对京津冀地区1992—2014年3阶段地面沉降进行了InSAR时序监测，也发现北京地区22年间地面沉降呈不断加重的趋势。Chen等（2016）分别利用ASAR和TerraSAR-X数据提取了北京2003—2011年的地表形变，得出最大形变超过100 mm/a的结论。上述研究均分析了北京地面沉降演化趋势，为沉降防控提供了参考。2014年12月南水北调中线工程正式运营，每年向京输水超1×108m3，改变了供水格局，也为沉降防控创造了有利条件（于海若等，2020；雷坤超，2023）。在北京水情发生显著变化的背景下，地下水深度的增加和新的集中开采区的出现，肯定会造成沉降分布出现新的特征，再加之北京近年地面沉降危险性评估工作鲜有研究，因此对南水北调工程运营后北京地区近年的沉降开展时序监测尤为必要。鉴于此，本文利用SBAS（Small Baseline Subset）时间序列InSAR技术对覆盖北京地区的2019年12月24日至2021年8月27日的52景Sentinel-1数据进行差分干涉处理，提取了北京近期的地面沉降信息，分析了沉降时空分布和演化特征，并开展了地面沉降危险性评估，为北京地区地面沉降防治等工作提供了参考依据。

1 研究区概况及数据来源

1.1 研究区概况

北京地处华北平原西北缘，地理范围为（115°25′—117°35′E，39°28′—41°05′N）（图1），总体地势西北高、东南低。西部山区属太行山山脉，北部山区属燕山山脉，东南部平原地区是由拒马河、永定河、潮白河、北运河、蓟运河五大水系共同作用形成的冲积扇群构成。目前探明的北京地区发育的活动断裂带主要有南口—孙河断裂、黄庄—高丽营断裂、顺义—良乡断裂、南苑—通县断裂、南口山前断裂、紫荆关断裂、永定河断裂、大兴凸起东源断裂、夏垫断裂、二十里长山断裂等。

图1 北京地区地理位置及SAR影像覆盖范围Fig.1 Geographical location of Beijing and SAR image coverage

1.2 数据来源

选取覆盖北京地区的52景142轨道的Sentinel-1A TOPS SAR数据进行干涉处理，时间跨度为2019年12月24日至2021年8月27日。Sentinel-1A数据的空间分辨率是5 m×20 m，C波段降轨VV极化，主要参数见表1。

表1 SAR数据基本参数Table 1 Basic parameters of SAR data

还获取了欧空局发布的POD精密定轨星历数据，该数据通过GNSS下行后处理获得，数据定位精度为5 cm，主要用于消除轨道误差。选取的数字高程模型（Digital Elevation Model，DEM）是由美国航空航天局提供的30 m空间分辨率的SRTM3 DEM数据，用于进行地理编码和去除地形相位误差。Sentinel-1A、POD精密定轨星历数据下载网址：https://scihub.copernicus.eu/dhus/#/home；DEM下载网址：https://dwtkns.com/srtm30m/。

2 SBAS-InSAR处理

小基线集方法（SBAS-InSAR）是Berardino等（2002）提出的一种用于监测地表形变时间演化的差分干涉测量算法，基本思想是将SAR数据按照一定的规则进行组合，使同一子集内干涉像对的时间和空间基线较小，而各子集之间的时间基线和空间基线都较大，得到一系列短空间基线差分干涉图。通过集合内的最小二乘和集合间的奇异值分解（SVD）法，得到整个时间序列地表形变信息的联合求解结果。假设有N+1幅覆盖同一地区的SAR影像，选择一幅作为主影像，将其他影像配准，设定合适的时间基线阈值和空间基线阈值，将时间与空间基线都小于该阈值的影像组成若干集合，其中每幅影像至少可与其他N幅影像中的一幅形成干涉像对。

设定空间基线阈值为150 m，时间基线阈值为36 d进行短基线干涉影像对组合，生成的111个干涉像对连接如图2所示，其中空间基线最小1.2 m，最大149.7 m，时间基线最小12 d，最大36 d。

图2 Sentinel-1数据干涉对时空基线分布Fig.2 Spatiotemporal baseline distribution of Sentinel-1 data interferograms

对所有干涉像对进行干涉处理，距离向方位向多视比设置为10:2，以保证数据解算中输出的相位图保持高分辨率。利用欧空局提供的POD精密定轨星历数据进行SLC（Single-Look Complex)数据轨道校正，去除轨道误差。用SRTM3 DEM数据去除地形相位，得到111个差分干涉图。采用Goldstein自适应滤波方法对差分干涉图进行滤波，减弱大气和各类噪声误差，采用最小费用流算法（Eineder et al，1998）进行相位解缠，相位解缠中相干性阈值设置为0.4，对于相干性阈值以下的区域不进行解缠计算，采用高低通滤波方法进一步消除大气延迟相位影响，得到差分解缠相位（Gabriel et al，1989；Massonnet et al，1993；Wright et al，2004）。检查所有相位解缠图质量，剔除相位解缠错误和质量较差的干涉像对。然后基于解缠的相位建立观测方程，采用奇异值分解方法（李永生等，2013）求解形变速率的最小二乘范数解，对求取的形变速率在时域上进行积分，得到形变时间序列结果。

3 北京地区InSAR时间序列结果及空间演化特征分析

3.1 北京地区InSAR地表形变监测

运用SBAS-InSAR技术提取北京地区2019年12月24日至2021年8月27日的地表形变信息。52景Sentinel-1A SAR数据提取的地表沿雷达视线方向的平均形变速率见图3，可以看出，提取的高相干点位空间分布密度高，典型形变区域的高相干点位呈现清晰、连续的分布状态，未发生明显空洞现象。需说明，InSAR技术提取到怀柔和密云部分地区存在微抬升现象，初步分析可能是数据处理过程中边缘效应导致的，实际并未发生。

3.2 空间分布特征分析

从2019年12月24日至2021年8月27日的地表平均形变速率可以看出，近垂向形变呈现明显的不均匀、不规则、漏斗式运动特征，沉降现象在朝阳、通州、大兴和海淀4个区比较显著，出现4个明显的片状沉降区域：通州沉降区[图4（a）]、朝阳金盏沉降漏区[图4（b）]、海淀稻香湖公园沉降区[图4（c）]和大兴机场周边沉降区[图4（d）]。监测到最大地表形变速率位于北京大兴安定垃圾清埋场，约-156 mm/a，沉降面积约88×104m2，该区域沉降暂未展现进一步向外侵蚀的趋势。下面具体分析4个沉降区的地表形变。

图4 主要沉降漏斗区空间分布特征Fig.4 Spatial distribution characteristics of the main subsidence funnel zones

3.2.1 通州沉降区。通州沉降区主要由通州次渠镇、亦庄火车站、北京环球影城、北马头村4个离散的沉降漏斗组成，呈现近NE向展布，沉降总面积达697×104m2，最大形变速率约-120 mm/a。其中，次渠镇漏斗约186×104m2，平均形变速率约-25 mm/a。亦庄火车站站前街道的京津高速桂家坟桥—麦庄桥周边区域沉降面积约60×104m2，最大形变速率达-120 mm/a。通州沉降区域中北京环球影城周边地区的形变面积最大，约371×104m2，最大形变速率达-30 mm/a。北马头村周边沉降区域面积约80×104m2，最大形变速率达-37 mm/a。

3.2.2 朝阳金盏沉降区。该沉降区由朝阳区金盏乡、黎各庄乡、马各庄环绕的不连续沉降漏斗组成，呈扇形分布，总沉降面积约4500×104m2，最大地表形变速率达-110 mm/a。其中，金盏乡沉降区主要由金盏乡、雷庄村和黑桥村3个子沉降漏斗组成，总沉降面积约1578×104m2，金盏乡、雷庄村和黑桥村3个子沉降区的最大形变速率分别约为-27 mm/a、-40 mm/a、-34 mm/a，雷庄村演播大楼附近地表形变速率最大，约-50 mm/a。黎各庄乡附近沉降范围涉及整个宋庄水井公园，面积约210×104m2，最大形变速率达-110 mm/a。马各庄周边地面沉降呈离散点状分布，面积约144×104m2，位于马各庄高安屯附近地表形变速率最大，达-82 mm/a，将形变速率结果叠加到谷歌地图卫星影像上进行目视解译推断，结果显示，马各庄周边地面沉降现象大概率由区域工程活动引起。

3.2.3 海淀稻香湖公园沉降区。海淀稻香湖公园沉降区位于海淀西北部，由苏家坨镇温阳路、稻香湖公园和南沙河3个子沉降区组成，呈三角形分布。此区域沉降面积较小，约114×104m2，最大形变速率约-30 mm/a，根据形变分布可以看出，稻香湖公园和南沙河子沉降区有进一步相向扩张的运动趋势。海淀区沉降速率最大的区域位于六里屯垃圾清埋场，达-152 mm/a，但该区域沉降现象仅发生在内部，暂未呈现进一步向外扩张的趋势。

3.2.4 大兴机场周边沉降区。大兴机场周边沉降区由北京大兴机场周边的20余个小范围沉降区组成，离散分布在大兴机场北侧、西侧和南侧，在机场东侧未见明显沉降现象。总沉降面积约420×104m2，最大形变速率约-80 mm/a。最近的义和场村沉降漏斗距机场东侧飞机跑道仅约300 m，因此北京大兴机场周边区域的地面沉降态势有必要重点关注。

4 北京地区地面沉降时间演化特征分析

选取6个特征点（图5）的形变量时间序列来表征北京典型区域地面沉降随时间的演化过程。特征点选取依据为：优先选择位于局部沉降速率较大的漏斗区且演化过程中未出现明显突跳的点。这种依据能使选点基本符合局部区域沉降的共性特征。下面具体分析6个特征点形变量时间序列。

图5 6个特征点形变量时间序列（a）首都机场跑道特征点1；（b）首都机场跑道特征点2；（c）金盏地区特征点；（d）故宫特征点；（e）高安屯特征点；（f）大兴机场特征点Fig.5 Time-series of deformation quantity for six feature points

（1）从首都机场跑道上2个特征点的形变量变化过程可以看出，首都机场部分飞机跑道正在遭受着地面沉降现象，2点位置在研究期内的平均形变量约-31 mm，形变量基本呈线性发展趋势，沉降速率分别是-19 mm/a和-22 mm/a，因此可知首都机场虽正在遭受地面沉降，但整体形变速率较小且较稳定，影响不大。

（2）金盏地区的地面沉降也基本呈线性发展趋势，累计形变量约为-166 mm，形变速率约-127 mm/a，形变速率较大，地面沉降形势较严重。

（3）故宫在研究期内垂向形变基本保持为零，因此故宫未发生地面沉降现象。

（4）高安屯地区也是一个明显的沉降漏斗，累计形变量持续加大，约-167 mm，平均形变速率达-121 mm/a。

（5）在研究期内大兴机场周边出现20余个离散且面积很小的沉降漏斗，选择其中一个漏斗分析其沉降量的变化，可以看出该特征点形变也呈线性趋势，累计形变量约-122 mm，形变速率约为-71 mm/a。

从6个特征点的形变量时间序列及地表形变空间分布可以看出，北京朝阳金盏地区和高安屯地区的地表形变速率较大，均超过-120 mm/a，值得进一步关注。

5 北京地区地面沉降危险性浅析

在地表平均形变速率基础上，根据北京地面沉降灾害危险性分区量化指标（表2）（杨艳等，2010；杨艳，2015）对北京地区地面沉降危险性进行初步分析。北京地面沉降危险性评价结果如图6所示。

表2 北京地面沉降危险性评价形变速率分区Table 2 Risk assessment of land subsidence in Beijing based on deformation rate zoning

图6 北京地表形变速率危险性分区Fig.6 Risk zoning of surface deformation rate in Beijing

从北京地区地面沉降危险性分区结果（表3）可知，在2019年12月24日至2021年8月27日期间，北京地区严重沉降区占比较小，仅为0.04%，主要位于海淀区六里屯垃圾清埋场、大兴区安定垃圾清理场和宋庄水井公园；较严重沉降区和一般沉降区占据的沉降面积比例约为6%，主要位于前述的通州沉降区、朝阳金盏沉降区、海淀稻香湖公园沉降区和大兴机场周边沉降区。这些区域的沉降速率均超过-30 mm/a，属于地面沉降的高危险区，进一步发育可能会对城市基础设施、重要建筑物造成不可挽回的损害。轻微沉降区和较稳定区占比分别为46.65%和47.15%，二者之和达93.80%，可认为在研究期内北京绝大部分地区未发生明显的地面沉降现象。

6 结论

在南水北调新水情下，北京地区地面沉降发生明显变化，为揭示地面沉降分布及其演化特征，利用SBAS-InSAR技术成功提取了北京地区2019年12月24日至2021年8月27日地面沉降信息，进而分析地面沉降时空分布特征和演化规律，开展地面沉降危险性评估分析。主要结论如下：

（1）2019年12月至2021年8月北京地区出现地面沉降现象，沉降分布在空间上呈现明显的不均匀、不规则、漏斗式运动特征。但总体上沉降面积占比较小，最大形变速率约-156 mm/a。

（2）地面沉降现象在朝阳、通州、大兴、海淀4个区比较显著，出现了4个明显的片状沉降区域，分别是通州沉降区、朝阳金盏沉降区、海淀稻香湖公园沉降区和大兴机场周边沉降区。通州沉降区4个沉降漏斗的地面沉降面积达697×104m2，最大平均沉降速率约-120 mm/a。金盏沉降区总沉降面积约4500×104m2，最大地表形变速率达-110 mm/a。海淀稻香湖公园沉降区沉降面积较小为114×104m2，最大形变速率约-30 mm/a。大兴机场周边沉降区沉降面积约420×104m2，最大形变速率约-80 mm/a。

（3）从典型形变点的形变量时间序列可以看出，北京朝阳金盏和高安屯地区形变量基本呈线性发展趋势，地表形变速率较大，均超过-120 mm/a，值得进一步关注。

（4）地面沉降速率分类结果表明，严重沉降区、较严重和一般沉降区占北京总面积的6.2%，主要位于上述4个地面沉降区，危险程度较高，值得重点关注，其他区域地面沉降危险性较低。