基于时序InSAR和GPS技术的北京平原区地表三维形变场特征*

雷坤超 马凤山 陈蓓蓓 罗 勇 崔文君 刘 贺 田 芳

(①中国科学院地质与地球物理研究所，中国科学院页岩气与地质工程重点实验室，北京 100029，中国)(②北京市水文地质工程地质大队(北京市地质环境监测总站)，北京 100195，中国)(③首都师范大学地面沉降机理与防控教育部重点实验室，北京 100048，中国)

0 引 言

许多地球物理现象(如构造运动、活动断裂、地震、地面沉降等)所引发的地表变形，在一定程度上均表现出三维形变的特征。准确获取这些地球物理现象三维形变场信息，对于查明其内在发生机制、运动过程和发展趋势，具有重要的地学数据支撑作用(陈强等，2012)。北京是国际上为数不多的以地下水作为主要供水水源的特大型城市之一，地下水开采量占到全市供水总量的50%～70%，地下水资源长期处于超量开采状态(1961～2018年，地下水资源累计亏损96亿立方米)，由此导致平原地区第四纪地层出现了较为严重的地表变形问题(雷坤超等，2019)。以往研究中，大多数学者针对超量开采地下水引发的多孔介质垂向固结压密变形，即地面沉降问题，利用大量观测资料，围绕地面沉降监测信息获取(张勤等，2009；雷坤超等，2014；周毅等，2016)、地面沉降演化规律(李曼等，2016；张永红等，2016；周超凡等，2017；张静等，2018)和成因机制(雷坤超等，2016；郭海朋等，2017；Gong et al.，2018；Chen et al.，2019)等方面开展了较为深入的研究，取得了丰富的研究成果。但对于抽取地下水引起的多孔介质固体骨架水平运动或变形特征问题研究较少。主要原因是由于多孔介质水平运动不如垂直固结压密变形或地面沉降那样直观、容易监测(王庆良等，2002)。但大量证据表明，抽水引起的含水层水平运动和变形确实存在(Helm，1994；Burbey，1999; Burbey et al.，2006；罗跃等，2018)。

目前，对于大尺度区域地表形变的监测，通常采用全球定位系统(Global Positioning System，GPS)和卫星雷达干涉测量(Interferometric Synthetic Aperture Radar，InSAR)方法。其中：GPS技术的优势在于其水平方向(E、N向)的地表位移具有较高的监测精度和灵敏度，采用较长的观测时段，能够探测出毫米级地表水平位移(Hu et al.，2006)。但在垂直方向的形变监测精度稍低，比水平精度低2～3倍(杨建图等，2006)。InSAR技术的优势在于其基于影像面的大范围测量，具有覆盖范围广、空间分辨率高、无接触遥感等特点，特别对于大尺度地表形变，可以有效获取空间连续的地表形变场信息(Galloway et al.，1998；Ferretti et al.，2001；Hooper，2008；王晨兴等，2018)。但由于SAR传感器侧视成像的特点，InSAR对于地表垂向形变较为敏感，而对于南北向的水平形变探测较为困难，由此使得该技术在三维形变信息获取及相应地球物理机制解释方面受到制约(陈强等，2012)。因此，综合GPS与InSAR两种技术各自的优势，将高精度的GPS水平位移测量与高密度的InSAR垂直形变观测相结合，对于获取完整的地表三维形变场信息，进一步认识地表变形的形成机制具有重要意义。

本文以北京平原区为研究对象，采用GPS与InSAR技术相结合，获取平原区时序地表三维形变场信息，分析其空间分布特征与演化规律。研究内容主要包括以下3个方面：(1)基于永久散射体干涉测量(PS-InSAR)和GPS技术获取平原区2013～2018年时序地表垂向形变场特征，分析其演化规律。(2)利用2013～2018年GPS同步观测数据进行各期联网解算和平差处理，获取各期GPS点在ITRF2005(International Terrestrial Reference Frame 2005)框架下和欧亚参考框架下精确三维坐标，分析不同参考框架下GPS点水平运动特征。(3)将欧亚参考框架下GPS点水平形变速率和InSAR垂向形变速率与不同含水层系统等水位线进行耦合分析，查明地下水开采对地表水平形变产生的影响。本文结果可为后续深入研究北京平原区地表三维形变成因机制及水-土三维耦合模型提供重要的数据基础与参考。

1 研究区概况

北京平原区第四纪地层是由5大水系(永定河水系、潮白河水系、北运河水系、大清河水系、蓟运河水系)联合作用形成的冲洪积扇群构成。在山前和冲洪积扇顶部，第四系厚度在20～40 m左右，为单一的砂、砂砾石层或顶部覆盖较薄的黏性土层。在冲洪积扇中下部和冲积平原地区，沉积物厚度逐渐增大，层次增多，颗粒逐渐由粗变细，岩性过渡为砂、砂砾石与黏性土层交错出现，并以黏性土为主。在部分沉积凹陷中心，第四系厚度达1000余米(雷坤超等，2016)(图1)。

图1 研究区位置和地质条件概况

与第四系沉积特征相对应，北京平原区含水层系统在垂向上主要分为3个含水层组：第1含水层组(潜水和浅层承压水)为第四系全新统(Q4)和上更新统(Q3)冲(洪)积物，含水层底板埋深在25 m和80～120 m左右。第2含水层组(中深层承压水)为第四系中更新统(Q2)地层，多层结构，岩性以中粗砂为主，部分地区含砾石，底板埋深在180 m左右。第3含水层组(深层承压水)为第四系下更新统(Q1)地层，多层结构，以中粗砂和砾石为主，底界为第四系基底，部分地区底界在260～300 m左右(图2)。

与第四系含水层组划分相对应，北京平原区可划分3个主要的压缩层组：第1压缩层组(Q4+Q3)，底板埋深小于100 m，广泛分布于平原区。岩性主要为全新统和上更新统冲积、冲湖积的粉土和黏性土。第2压缩层组(Q2)，底板埋深150～180 m，主要分布在冲洪积扇中下部、冲积平原地区。岩性主要为中更新统冲洪积、冲湖积粉土和黏性土层。第3压缩层组(Q1)，底板为第四系基底，部分地区在260～300 m左右，主要分布在几大沉积凹陷地区。岩性主要为下更新统河湖相沉积的黏性土层(贾三满等，2007)(图2)。

图2 A—A1处水文地质剖面

2 地表形变监测方法

2.1 时序PS-InSAR监测方法

永久散射体干涉测量(Permanent Scatterer for SAR Interferometry，PS-InSAR)方法主要是从长时间序列的SAR影像中选取那些保持高相干的像素点作为研究对象，利用它们的散射特性具有很好稳定性的特点，获得可靠的相位信息。通过对相位进行时空分析，分解各相干目标点上相位的组成，主要包括高程、地表形变以及由于大气引起的相位变化，最后得到地表的位移。该方法可以有效克服DinSAR在空间和时间方面失相干以及大气延迟的影响，从而提高地表形变的监测精度，达到毫米级(Ferretti et al.，2001)。主要步骤包括：SAR影像配准、差分干涉处理、PS点选取、地表形变速率提取、投影变换与精度验证(雷坤超等，2016)。文中采用GAMMA软件对覆盖北京平原区54景RadarSAT-2 雷达卫星 SAR影像(Wide strip模式，C波段，分辨率30 m×30 m，入射角为27.8°，2013/01/16～2018/12/22)，进行时序差分干涉处理及PS点提取，剔除了基线距大于300 m的干涉像对，利用30 m分辨率SRTM数据为参考DEM，进行地形相位去除，经地理编码及视线向(LOS)到垂直向投影变换后，获取了北京平原区2013～2018年时序地表垂向形变速率。

2.2 GPS监测方法

北京市自2008年至今，每年8～10月均对平原区地表形变GPS观测点进行联网观测。结合北京房山(BJFS)、北京十三陵(BJSH)和天津蓟县(JIXN)3个基岩基准点和平原区14个CORS(Continuously Operating Reference Station)站，采用美国麻省理工学院研制的GAMIT/GLOBK软件包进行联网解算和平差处理。每期(年)GPS数据处理均包括61个GPS与水准测量一体点(GPS观测墩底部同时建设水准测量点)和14个CORS站，共计75个GPS同步观测数据。文中利用上述GPS处理方法，获取了平原区2013～2018年共6期(年)GPS点在ITRF2005框架下精确三维坐标。同时，为了更好地研究北京区域内部GPS点水平运动特征，以3个基岩基准点为稳定参考基准，采用NNR-NUVEL1A模型，将ITRF框架下的水平速度场扣除欧亚板块速度，得到所有GPS点在欧亚参考框架下的区域内部水平速度场。

3 平原区地表垂向形变特征

3.1 区域地面沉降特征

图3 利用InSAR和GPS技术获取北京平原区2013～2018年地表垂向形变特征

图4 北京平原区2013～2018年平均沉降速率、沉降区面积和体积统计

通过对北京平原区时间序列地面沉降分布特征及演化规律分析发现，2013～2018年地面沉降总体呈减缓趋势，但各地沉降减缓时间存在差异性。平原区北部在2016年开始减缓，至2017年和2018年沉降减缓现象越发明显；平原区东部在2018年出现明显减缓。出现这种变化趋势的原因主要包括3个方面：(1)南水北调工程产生的环境正效益。2014年底，南水正式进京，每年向北京输水超过10亿立方米，有效缓解了北京市水资源紧缺形势，为北京市生产、生活供水及地表水-地下水联合调蓄提供了重要的水源基础。(2)一系列地面沉降防控措施初见成效。南水进京后，北京市先后对密怀顺应急水源地进行地下水压采，利用南水在冲洪积扇顶部开展地下水回补，同步实施城区自备井置换等，致使地下水开采量逐年下降，由2013年的21.10亿立方米下降到2018年的16.20亿立方米。2016～2018年3年间，部分地区含水层组水位出现不同程度回升，特别是冲洪积扇顶部(平原区北部)最为明显。(3)平原区北部与东部地区在地层岩性及结构特征方面差异明显。北部地区主要是单一砂砾石层，对地下水位变化响应较为敏感，水位抬升后，地层几乎同步回弹或沉降明显减缓。而东部地区主要由砂、砂砾石和黏性土层组成的多层结构区，并以黏性土层为主，层次较多，颗粒较细，这种地层结构对地下水位变化的响应存在较为明显的滞后性。

3.2 地面沉降中心特征

图5 2013～2018年主要沉降中心地面沉降速率变化特征

图6 2013年东部4个沉降中心剖面和2013～2018年平原区沉降中心速率统计

3.3 累计地面沉降特征

在累计地面沉降变化方面，文中利用InSAR、GPS和水准测量3种手段，获取了各主要沉降中心2013～2018年时间序列累计沉降量曲线(图7)。由于选取的PS点和GPS-水准一体点不能完全重合，故以GPS-水准一体点为中心，以100 m为半径，选取邻近PS点进行统计分析。结果显示，随着时间的推移，各沉降中心累计沉降量逐渐增大。2013～2018年朝阳金盏(P1)最大累计沉降量达816.77 mm，是北京市累计地面沉降量最大的地区。

图7 2013～2018年InSAR、GPS和水准测量获取的累计沉降量对比图

其次是黑庄户地区(P2)最大累计沉降量792.10 mm；朝阳三间房(P3)、通州城区(P4)、昌平八仙庄(P5)和海淀上庄(P6)最大累计沉降量分别为749.00 mm、627.00 mm、456.00 mm和437.00 mm。通过InSAR、GPS分别与水准测量结果对比发现，2013～2018年InSAR累计沉降量与水准结果偏差在4～30 mm之间；GPS累计沉降量与水准结果偏差在5～50 mm之间，但3种手段获得的累计沉降量在发展趋势上具有很好的一致性。

3.4 InSAR和GPS垂向形变精度分析

为定量评定InSAR和GPS两种手段垂向形变监测精度，在InSAR和GPS获得的2013～2018年沉降速率图上选取40个相同点位同期水准测量数据进行单点精度检验(表1、图8)。需要说明的是，GPS测量高程值属于大地高系统，水准测量属于正常高系统，两者之间存在高程异常。但考虑到研究区较小，比较的是垂向形变量(即高程差值)，因此两者之间高程异常可以忽略不计。同时，InSAR获取的是视线向(LOS)形变量，在与水准测量结果比较时，需将视线向(LOS)形变量投影变换到垂直方向。经对比分析发现：InSAR监测结果与水准测量结果吻合程度较高，表现出相同的形变特征，两者之间的互差均在2～8 mm范围之间，互差均方差为4.33 mm，并且具有显著的线性相关性，相关系数(R2)达到0.983。同时，将GPS测量结果与水准结果对比分析发现，80%以上的GPS点与水准点互差不大于10 mm，个别点位上两者之间互差较大，最大为15.80 mm，互差均方差为7.56 mm，相关系数为0.946。由此说明，本次InSAR和GPS测量成果均与水准测量具有很好的一致性，单点监测精度较高，可以较为详细地揭示出北京平原区地表垂向形变特征。

表1 InSAR、GPS垂向形变量与水准测量结果对比

图8 利用水准测量结果对InSAR和GPS垂向形变量进行精度检验

4 平原区地表水平形变特征

4.1 ITRF2005参考框架下GPS水平速度场

图9 ITRF2005参考框架下北京平原区GPS点水平速度场

4.2 GPS水平速度场精度分析

表2 ITRF2005框架下GPS观测点水平运动速率及精度

4.3 欧亚参考框架下GPS水平速度场

图10 欧亚参考框架下北京平原区内部GPS点水平速度场

图11 地下水位漏斗与地表三维形变响应关系

5 结 论

本文将高精度的GPS水平位移测量与高密度的InSAR垂直形变观测相结合，获取了北京平原区时序地表三维形变场信息，分析了其分布特征与演化规律。研究结果表明：

(1)北京平原区在抽水引发的第四系附加应力场作用下，地表呈现出显著的三维变形特征，以垂向变形为主，并辅以水平向位移。

(4)在欧亚参考框架下，GPS点水平运动速率明显减小，各点之间非一致性变化较为明显，各点运动速率的大小和方向不具备整体趋势性活动特征。特别是几大活动断裂交接部位的地面沉降严重区，往往也是GPS点水平运动速率较大的地区。GPS点水平运动方向总体指向地面沉降或地下水位降落漏斗中心，或由地下水位较高的地区指向地下水位较低的地区。这主要是抽取地下水导致第四系含水层系统在水平向产生的变形分量引起的。