雄安新区2017—2019年地面沉降SBAS-InSAR监测与分析

冉培廉，李少达，戴可人,2，杨晓霞，曾森,3，冉中鑫,4

(1.成都理工大学 地球科学学院，四川 成都 610059；2.成都理工大学 地质灾害防治与地质环境保护国家重点试验室，四川 成都 610059；3.四川省自然资源厅 信息中心，四川 成都 610059；4.长江水利委员会水文局 汉江水文水资源勘测局，湖北 襄阳 441100)

0 引 言

雄安新区于2017年4月1日在中共中央和国务院指导下成立。该区由雄县、容城县、安新县及周边部分地区合并而成，地处北京市、天津市和保定市腹地，是河北省的国家级新区，地热资源丰富，开采条件良好，主要集中在容城凸起、牛驼镇凸起和高阳低凸起，且已有多年的开采历史[1-3]。地处暖温带大陆性季风气候，年平均降水量为516 mm，水资源总量1.73×108m3，其中90%以上为地下水[4]。因地下水长期过度开采，华北平原遭受了较严重的地面沉降灾害，如北京市[5-6]和天津市[7-8]。随着城市快速发展，城市人口增加，地热资源和水资源需求急剧上升，将导致地面沉降发展。及时掌握雄安新区地面沉降状况，可为该区的健康发展和长期规划提供重要参考。

合成孔径雷达干涉(synthetic aperture radar interferometry,InSAR)测量技术是一种有效监测地面微小形变的新型遥感技术，与传统监测技术相比，InSAR技术具有覆盖范围广、全天候、低成本等优点，已广泛应用于城市地面沉降监测[9-11]。2018年前仅有覆盖该区的地面沉降研究结果，没有对雄安新区地面沉降的原因进行详细分析，如：葛大庆等[12]基于SBAS-InSAR技术，利用RADARSAT-2数据获取华北平原2012—2015年的地面沉降结果，但缺少对雄县、容城县、安新县及高阳县区域内地面沉降情况的专门描述；张永红等[13]基于SBAS-InSAR技术，利用ERS，ENVISAT，RADARSAT-2数据获取京津冀地区1992—2014年的地面沉降结果，不仅利用水准数据验证SBAS-InSAR监测结果的精度，还对北京市和天津市22年间地面沉降的时空演变特征进行分析，但是未对雄县、容城县、安新县及高阳县区域内的地面沉降进行专门分析；张永红等[14]基于SBAS-InSAR技术，利用RADARSAT-2数据监测雄安新区2012—2016年的地面沉降结果，利用相近时段的水准测量数据验证SBAS-InSAR监测结果的精度，不仅详细分析了该区地面沉降的原因，而且还对地面沉降危险性进行了评价；DAI K R等[15]对雄安新区2017—2018年地面沉降进行监测，指出该区地面沉降分布与区域内地热井开采的空间分布有关。

综上所述，根据已有的地面沉降研究结果，雄安新区存在较大范围地面沉降，同时其与地热能及地下水开采分布的关系值得进一步探究。因此，本文基于SBAS-InSAR技术，利用2017—2019年覆盖雄安新区的32景Sentinel-1B影像对该区地面沉降进行监测，获取该区地面沉降速率和累积沉降量，分析地面沉降与地热资源和地下水之间的关系，进而揭示雄安新区地面沉降状况。

1 研究区概况及试验数据

雄安新区位于华北平原，毗邻北京市、天津市和保定市，是河北省的国家级新区。该区总面积1 576.6 km2，境内有白洋淀、孝义河、府河、漕河等水系，但依然缺乏淡水资源。雄安新区成立后将会修建大量基础设施，且由于地下水的长期过度开采，华北平原存在较严重的地面沉降，为了及时掌握该区地面沉降状况，基于SBAS-InSAR技术，利用Sentinel-1B(1w)数据监测雄安新区2017—2019年地面沉降状况。由于Sentinel-1B干涉宽幅影像覆盖面积为(250×250)km2，远大于本文研究区域，为了提高运算速度和监测结果精度，对合成孔径雷达(synthetic aperture radar,SAR)影像进行裁剪得到本文研究区域，图1中蓝色矩形为研究区影像范围。

图1 研究区位置Fig.1 Location of study area

本文使用欧洲航天局网站上2017—2019年覆盖雄安新区的32景(按时间顺序t0,…，t31排列)Sentinel-1B干涉宽幅影像，采用C波段(波长5.6 cm)，分辨率5 m × 20 m，并采用单极化(VV)观测方式，Sentinel-1B降轨影像数据基本参数信息如表1所示。采用Sentinel-1B卫星的精密轨道数据(precise orbit ephemerides，POD)修正轨道误差；采用日本宇航局网站上ALOS World 3D(AW3D30)数字表面模型(digital surface model,DSM)作为外部数字高程模型(digital elevation model,DEM)，以去除地形效应。

表1 Sentinel-1B IW模式影像数据基本参数信息Tab.1 Basic parameter informations of Sentinel-1B IW mode image data

2 数据处理

短基线干涉测量技术(small baseline subset interferometry,SBAS-InSAR)是由P.Berardino等[16]于2002年提出的一种用于监测地面微小形变的新型时间序列分析方法。SBAS-InSAR技术能有效克服合成孔径雷达差分干涉(differential interferometric synthetic aperture radar,D-InSAR)测量技术在数据处理过程中的限制因素(时空失相干、大气效应等)，可连续监测地面微小形变[17-18]，监测精度可以达到毫米级。SBAS-InSAR技术流程如图2所示。

图2 SBAS-InSAR技术流程Fig.2 Flow chart of SBAS-InSAR technology

选择32景Sentinel-1B干涉宽幅影像中任意一景影像作为超级主影像，其余影像与超级主影像进行配准，设置时间基线阈值(temporal baseline threshold,TBT)和空间基线阈值(perpendicular baseline threshold,PBT)，N+1幅SAR影像最多可以生成M幅差分干涉图，并满足关系式(1)[16]：

(1)

设置时间基线阈值为120 d，空间基线阈值为154 m，组合生成132幅干涉对。干涉对时空基线分布如图3所示。采用最小费流法(minimum cost flow,MCF)对132幅干涉对进行相位解缠，假设从影像tA和主影像tB时刻获取SAR影像生成第j幅差分干涉图，像素的任意干涉相位可以用SAR坐标表示为

图3 时空基线Fig.3 Time-space baseline

(2)

相位解缠完成后，线性形变和DEM误差相位之间会建立一个新的SBAS线性方程，即基于线性模型估算残余高程和形变速率，未解缠的相位将进行二次解缠，其表达式为

[B,C]p′=δφ,

(3)

式中：B为M×N的矩阵；C为与空间基线距离相关的系数矩阵。p′满足关系式(4)：

p′=[vΔZ]T,

(4)

其中，v为平均形变速率的相位值矩阵。

首先，从生成差分干涉图到产生残余相位的过程中，去除线性形变相位和DEM误差，包括时空失相干、非线性形变相位和大气延迟相位；其次，空间上利用高通滤波，时间上利用低通滤波，将大气延迟相位和非线性形变相位分离；然后，运用最小二乘法(least square method,LSM)和奇异值分解法(singular value decomposition,SVD)对所有的干涉对进行解缠，从而得到时间序列的形变结果；最后，基于地理编码将SAR坐标系下的累积沉降量和形变速率转换到地理坐标系,得到地理坐标系下的累积沉降量和形变速率。

3 结果与分析

3.1 地面沉降现状

从图4可以看出，研究区内存在2个沉降区(负值代表地面沉降，正值代表地面抬升)，分别位位于雄安新区北部(图4(b))和南部(图4(a))。北部沉降区位于雄县，沉降面积约98 km2，沉降速率超过80 mm/a,南部沉降区位于高阳县，沉降面积约152 km2，沉降速率超过90 mm/a，南部地面沉降比北部严重。假设以2017年6月12日的SAR影像为参考影像，即形变量为0 mm，余下SAR影像以该影像为基准，在时间序列中会产生累积形变。由图5可以看出，从2017年6月12日至2019年8月25日，随着时间推移，地面沉降现象逐渐明显，沉降区范围也不断扩大，呈现出不断加剧沉降趋势。为进一步探讨雄安新区2017—2019年的地表形变时间序列情况，在南部沉降区选取6个采样点，依次编号为1～6，在北部沉降区选取4个采样点，依次编号为7～10，在中部选取5个采样点，依次编号为11～15，采样点具体分布如图4所示。由图6可知，南部沉降区(采样点1～6)的累积沉降量分别为130，132，184，180，150，125 mm，北部沉降区(采样点7～10)的累积沉降量分别为170，180，140，129 mm，中部区域采样点11，13～15的累积沉降量分别为12，5，8，20 mm，采样点12的累积抬升量为2 mm，其中采样点3～4，7～8，10起初表现为小幅抬升，此后虽偶有起伏，但不影响整体下沉趋势，采样点1～2，5～6，9均表现为整体下沉趋势，采样点12表现为向上缓慢抬升，采样点11，13～15仅发生微小沉降变化，采样点11～15均表现为稳定发展趋势。比较年平均形变速率图与以往研究结果[14-15]，沉降区位置和沉降区边界表现出良好的一致性，说明本文InSAR监测结果可靠。

图4 年平均形变速率图Fig.4 Annual average deformation rate image

图5 雄安新区2017年6月至2019年8月累积沉降量序列Fig.5 Accumulative subsidence sequence of Xiongan New Area from June 2017 to August 2019

图6 采样点时间序列形变Fig.6 Time sequence deformation of sampling points

3.2 地面沉降原因分析

雄县是华北平原地热资源最丰富的地区，地热资源面积约320 km2，占雄县面积的60%[20]。自20世纪70年代首次发现地热资源以来[21]，地热资源一直被不断地开发，至2017年雄县的地热井已达到68眼[14]，且随着经济不断发展，地热井数量仍在迅速增加，地热水资源的消耗将加剧雄安新区的地面沉降。

雄安新区的地面沉降与地热资源的开发密切相关。由图7可以看出，雄县地热井分布在雄县东部和北部，2017年6月至2019年8月雄县东部未探测到明显地面沉降，雄县北部(图4(b))地面沉降严重，沉降面积约为98 km2，累积沉降量高达180 mm，地热井较多。由此可见，雄安新区北部狭长的沉陷区(图4(b))与地热井分布(图7)具有较高的相关性。此外，雄县东部地区未开发地热[22]，未发现地面沉降。雄安新区南部沉降区(图4(a))位于高阳县北部，地面沉降面积约152 km2，累积沉降量高达184 mm。雄安新区南部沉降区处于地热田开采范围(图8)，地热田内存在地热井(图7)。由此可见，地热资源开采会对地面沉降产生一定影响，以上研究结果可为雄安新区地热资源和地下水开发的长期规划提供重要参考。

图7 2019年雄安新区地热井位置(修改自文献[21])Fig.7 Geothermal well locations of Xiongan New Area in 2019(modified from reference[21])

图8 2017年雄安新区地热田分布图(修改自文献[23])Fig.8 Geothermal field distribution map of Xiongan New Area in 2017(modified from reference[23])

其次，塑料包装业和纺织业消耗大量的地下水，地下水过度开采是造成雄安新区地面沉降的又一重要因素。雄县是华北平原最大的塑料包装业基地[24-25]。由于华北平原水资源较匮乏[26]，且塑料包装业在生产过程中耗费大量的水资源，导致华北平原的地下水过度开采。高阳县纺织业是该县的支柱产业[27-28]，纺织业亦属于高耗水行业。因此可以推测地下水过度开采与雄安新区地面沉降有关。

4 结 论

(1)雄安新区存在2个地面沉降区：北部沉降区和南部沉降区。北部沉降区沉降速率超过80 mm/a,沉降面积约为98 km2；南部沉降区沉降速率超过90 mm/a，沉降面积约为152 km2。

(2)雄安新区地面沉降与地下水和地热资源的过度开采有关。北部沉降区可能由地热资源过度开采造成，南部沉降区可能由地下水过度开发造成。

(3)由于未获取雄安新区2017—2019年的水准数据，无法直接验证SBAS-InSAR技术获取的雄安新区2017—2019年地面沉降结果精度，但通过与现有研究成果比较分析,可间接验证。本文获取的沉降区位置和沉降区边界表现出良好的一致性，表明利用SBAS-InSAR技术获取雄安新区2017—2019年的地面沉降结果具有一定可靠性。