Sentinel-1A京津冀平原区2016—2018年地面沉降InSAR监测

史 珉， 宫辉力， 陈蓓蓓， 高明亮， 张舜康

(1.首都师范大学地面沉降机理与防控教育部重点实验室，北京 100048； 2.首都师范大学水资源安全北京实验室，北京 100048； 3.首都师范大学北京市城市环境过程与数字模拟国家重点实验室培育基地，北京 100048； 4.首都师范大学三维信息获取与应用教育部重点实验室，北京 100048)

0 引言

地面沉降是一种地表高程下降的环境地质现象，其发生发展增加了城市内涝、海水倒灌等风险，并可以诱发地面塌陷、地裂缝等系列环境灾害，形成灾害链。区域地面沉降不仅对城市基础设施、高速铁路等重大工程产生潜在威胁，还影响着城市地下空间资源利用效率，成为制约区域可持续发展的重大问题。

京津冀地区所处的华北平原是我国沉降速率最大、影响面积最大的地区[1]。自1923年首次在天津发现地面沉降以来[2]，华北平原地面沉降按期发展过程可分为4个阶段： 局部孕育阶段(1954年之前)、加速发展阶段(1955—1974年)、急速扩张阶段(1975—1984年)和控制减缓阶段(1985年至今)[3]。截至2012年，华北平原沉降速率大于20 mm/a的区域面积约占其总面积的12.5%，其中90%的沉降区位于北京、天津、河北三地，局部地区累计沉降量超过3.4 m[4]。京津冀协同发展地质调查中将地面沉降列为京津冀一体化进程中需要关注的重大地质问题之一。对其地面沉降进行监测，认识其分布规律、演化特征，对京津冀城市和重大工程规划建设，为实现京津冀协同发展的重大国家战略提供地质环境保障。

合成孔径雷达干涉测量技术(interferometric synthetic aperture Radar，InSAR)已被广泛应用于地表形变监测中。与水准测量、分层标、全球定位系统(global positioning system，GPS)测量技术相比，InSAR技术具备大范围、全天时全天候获取高时空分辨率地表形变的能力。多时相InSAR(multi-temporal InSAR，MT-InSAR)技术的发展(如永久散射体干涉(persistent scatterer InSAR，PS-InSAR)、小基线 (small baseline subsets，SBAS)等)更是克服了大气延迟、时空失相干对传统InSAR方法的影响[5-6]。目前，MT-InSAR技术已被广泛应用于京津冀地面沉降研究中，但这类研究多关注北京[7-11]、天津[12-15]、沧州[16-18]和廊坊[19-20]等城市区域，通过多源SAR数据获取其地面沉降信息，而对京津冀平原地面沉降的分布、范围及沉降程度缺乏整体性认识。针对单轨SAR数据源幅宽对区域地面沉降监测范围的限制问题，本研究基于多轨SAR数据融合技术，获取了2016—2018年京津冀平原地面沉降信息； 为保证InSAR监测结果的可靠性，选用同时期水准测量数据对InSAR结果进行精度验证，并对相邻轨道结果进行一致性检验； 并结合GIS空间分析、剖面分析等方法，获取InSAR监测时段内京津冀平原区地面沉降演化特征。

1 研究区概况及其SAR数据源

京津冀地区(N36°05′～42°37′，E113°11′～119°45′)位于我国华北地区，包含北京、天津和河北省的11个地级市，属暖温带大陆性季风型气候，降水多集中在7—8月。该地区地形条件复杂，总体地势为西北高、东南低，自西北向东南依次分布坝上高原、燕山—太行山山地和东南部平原[21]。京津冀地区存在多种环境地质问题，其中以地面沉降等为例的缓变性地质主要发生在平原区。京津冀平原区是华北平原的重要组成部分，自山麓到海岸，按成因和形态特征可分为第四纪山前冲积、洪积倾斜平原(山前平原)，中部冲积、湖积多层叠加平原(中部平原)和东部冲积、湖积夹海积的滨海平原(滨海平原)[22]。截至2015年，京津冀平原区累积沉降量大于200 mm的面积约为6.4万km2[21]，约占平原区面积的73%，沉降速率大于50 mm/a的严重沉降区面积占全国严重沉降区面积的92%(http: //www.tianjin.cgs.gov.cn)。

Sentinel-1A卫星是欧洲空间局于2014年发射的重访周期为12 d的C波段对地观测卫星。该卫星使用近极地太阳同步轨道，轨道高度、倾角和周期分别为693 km，98.18°和99 min，具多种成像模式和极化方式(https: //sentinel.esa.int/)。单轨SAR影像受其幅宽限制，监测范围有限，为获取京津冀平原区的地面沉降信息，本次研究共选用了3轨时间跨度为2016—2018年的Sentinel-1A数据(Track 40，Track 142和Track 69)。所用3轨Sentinel-1A数据均为宽幅干涉模式(interferometric wide，IW)获取的单视复数(single look complex，SLC)升轨存档影像，其极化方式均为垂直极化(vertical-vertical，VV)，空间分辨率为5 m×20 m，单幅幅宽达250 km。影像覆盖范围及主要参数分别如图1和表1所示。

图1 研究区概况示意图Fig.1 Location of study area

表1 Sentinel-1A影像信息Tab.1 Satellite information for the data used in this study

2 研究方法

2.1 MT-InSAR处理过程

通过SARPROZ软件(https: //www.sarproz.com/)对3轨Sentinel-1A影像分别进行PS-InSAR处理，以获取每轨SAR影像范围内时序地面沉降信息。处理过程包括以下几个主要步骤： 主影像选取，数据集配准及重采样，基线建立，数字高程模型(digital elevation model,DEM)模拟，生成差分干涉图，选取PS候选点，相位解缠，大气相位估计和去除，重新选取PS点和估计PS点形变速率[23-24]。由前人研究成果已知，京津冀地区水平向形变值约为1～3 mm/a[25-28]，因此在本次研究中忽略水平向形变结果，通过式(1)将视线向(line-of-sight，LOS)InSAR形变结果转为垂向，并将经过处理后的垂向结果用于与水准结果的精度验证和后续的SAR影像结果拼接中，即

dv=dLOS/cosθ

(1)

式中：dv为垂向形变结果，mm；dLOS为LOS向形变结果，mm；θ为雷达LOS向与垂向的夹角。

2.2 多轨SAR数据拼接

为了获取大区域的地面形变信息，需对上文所得单轨SAR影像结果进行拼接。首先对各轨SAR影像垂向形变结果进行了一致性检验和精度验证，在保证数据精度可靠的情况下，对相邻轨道数据进行融合。在相邻轨道沉降速率融合中，Ge等[29]通过引入与参考点距离权重的关系进行整体平差; 孙赫[30]通过影像重叠区内PS点形变量差值众数对相邻轨道数据结果进行校正; 熊思婷等[31]比较了区块法和插值法在求解异轨重叠区的形变差中的应用; 张永红等[32]通过大量的水准观测数据对SAR数据结果进行了平差控制。本次实验中，依据重叠区内同名高相干点的形变速率一致准则[32-33]，计算平差权重。

由于受不同轨道卫星成像几何差异，不同轨道SAR数据所得LOS向地表形变结果间存在系统性偏差。为消除不同视角投影的影响，首先依据式(1)将InSAR监测所得各轨LOS向形变值转为垂向形变值。如图1所示，Track 142与Track 69和40间都具有重叠区，因此在多轨SAR数据融合时，将Track 142形变结果作为基准。Track 142与Track 40影像重叠区域内，雷达LOS向与垂向的夹角分别为33.72° 和43.75°； Track 142与Track 69影像重叠区内LOS向与垂向夹角分别为43.77°和33.67°。根据上述θ值，将重叠区域内每轨影像的各PS点LOS向结果依据式(1)分别投影到垂直向上，选取的PS点的相干性均在0.9以上。之后，选用最邻近法对重叠区内各轨PS点进行同名点匹配。最后，采用最小二乘法(least square method, LSM)构建相邻轨道重叠区域的平差方程，并分别对Track 40和Track 69垂向形变结果进行偏移量补偿。该方法是通过实际观测值与模型计算值的误差的平方和最小来寻找最优的匹配函数。假设相邻轨道分别为x和y，其中x为参考基准，y为需要进行平差的数据，相邻轨道之间存在n个重合点，基于LSM的原理，则

(2)

此时，通过对δ求极限，即可求得a和b的值。最终，获得多轨SAR数据融合结果。

3 结果及精度验证

3.1 京津冀2016—2018年地面沉降监测结果分析

为获取京津冀平原区地面沉降信息，基于MT-InSAR技术分别获取了3轨Sentinel-1A数据源地面沉降信息，后对相邻轨道数据进行融合获取了京津冀平原区内2016—2018年地面沉降速率(图2)。为了显示清晰，对结果进行了栅格化处理，像元大小为1 km。InSAR监测结果显示， 2016—2018年间京津冀平原区内最大沉降速率达164 mm/a，分布在北京朝阳—通州沉降区和廊坊—天津沉降区内。研究区内沉降分布范围广，且空间分布不均，在北京、天津及河北省各地级市均有不同程度的地面沉降问题，且三地沉降区有连片趋势。其中北京平原区沉降较为严重的地区主要分布在北京东南部朝阳、通州等地，且北京东南部沉降区已扩张到廊坊三河市界内。天津、廊坊两地沉降较为严重的区域分布在两地交界处。河北平原沉降较为严重的区域主要分布在其中东部地区，且该地区沉降区呈南北贯通趋势。研究区内共获得了1 426 967个PS点，其中2016—2018年沉降速率大于20, 50, 100 mm/a的点数分别占总点数的53.9%, 23.9%和3.6%。沉降速率大于20, 50, 100 mm/a的面积分别为3.3万, 1.4万和0.14万km2。据统计，2016年沉降速率大于20, 50, 100 mm/a的点数分别占总点数的53.8%, 23.7%和3.6%； 到2017年沉降速率大于20和50 mm/a的点数略有增长，分别占总点数的54.4%和23.8%，沉降速率大于100 mm/a的点数占比则降至3.5%。2018年沉降速率大于20和50 mm/a的PS点数与2017年比有所增长，占比达57.2%和24.0%，沉降速率大于100 mm/a的点数与2017年相比保持稳定。

图2 2016—2018年京津冀平原区PS-InSAR地表形变速率Fig.2 Mean displacement velocities throughout the Beijing-Tianjin-Hebei derived from the Sentinel-1A data by using PS-InSAR from 2016 to 2018

研究区内地面沉降大于50 mm的区域主要分布在北京东南部，廊坊北部，廊坊东部到天津西部，天津南部，保定东部到廊坊西部，保定西部，唐山—秦皇岛沿海区域，以及河北平原东南部区域。对监测时段内研究区各子平原内沉降速率大于50 mm/a的严重沉降区面积进行统计发现，山前冲积洪积平原沉降量大于50 mm的面积有轻微增长(2016—2018年)，由2 937 km2变为2 982 km2，涨幅为1.5%。2016—2018年中部平原沉降量大于50 mm的面积的涨幅为0.3%，基本维持稳定。2016—2018年，严重沉降区的面积分别为10 547 km2，10 574 km2和10 577 km2。与2016年相比，滨海平原严重沉降区的面积在2018年增长较为明显，由382 km2变为577 km2，涨幅达33.8%。其中唐山—秦皇岛沿海沉降区范围明显扩大(如图3所示)。为进一步分析研究时段内地面沉降发展过程，对沉降较严重区

(a) 2016年 (b) 2017年 (c) 2018年

图3 2016—2018年各年度InSAR平均形变速率Fig.3 Mean displacement velocities derived by InSAR between 2016 and 2018

域沿东西、南北方向分别进行剖面分析，剖面线位置如图3(c)中虚线所示。2016—2018年，不同地面沉降区的沉降空间演化存在显著差异(如图4所示)。其中位于唐山沉降区的TS剖面线处(图4(e)和(f)所示)的沉降速率在监测时段呈增长趋势。而位于北京、廊坊—天津、保定—沧州—衡水、衡水—邢台、邯郸沉降区的BJ，LT，BCH，HX和HD剖面处2016—2018年沉降速率基本保持一致，沉降中心发育较为稳定。综上，InSAR监测时段内山前平原和中部平原地面沉降发展相对稳定，滨海平原地面沉降范围及沉降速率在监测时段内增长明显，尤其是唐山沿海地区，沉降速率涨幅最为明显。

(a) BJ西—东 (b) BJ北—南

(c) LT西—东 (d) LT北—南

(e) TS西—东 (f) TS北—南

(g) BCH西—东 (h) BCH 北—南

(i) HX西—东 (j) HX北—南

(k) HD西—东 (l) HD北—南

图4 剖面时间序列InSAR形变速率结果Fig.4 Profiles of the InSAR displacement velocities

为了获取2016—2018年研究区内地面沉降时序特征，选取了8个不同位置上的沉降特征点，各沉降点的位置如图3(b)所示，其对应时序特征如图5所示。北京朝阳金盏、天津王庆坨、保定大崔营、保定高阳、邢台南宫和邯郸临漳等沉降区，形变特征以线性下沉为主，最大累积沉降量超过300 mm。 在唐山丰南、乐亭沉降区，区域时序沉降具有明显的季节性特征，在2016—2018年沉降呈下降趋势，如图5(g)和(h)所示，但在2016年9月—2017年1月和2017年9月—2018年1月沉降减缓甚至出现回弹现象，之后地表继续下沉，至2018年10月地面累计沉降量达180～254 mm。

(a) 北京朝阳 (b) 天津王庆坨 (c) 保定大崔营 (d) 保定高阳

(e) 邢台南宫 (f) 邯郸临漳 (g) 唐山丰南 (h) 唐山乐亭图5 研究区内沉降特征点对应时序特征Fig.5 Accumulative time-series deformation revealed by the InSAR technique

3.2 精度验证

为保证InSAR数据结果的准确性，选用了26个2017年9月—2018年9月的水准观测值对Sentinel-1A数据所得垂向形变速率进行精度验证。水准数据位置如图1所示，位于Track 142和Track 69影像覆盖范围内，因此通过水准测量数据对这2轨影像所得形变结果进行验证。在验证过程中，以水准点为圆心，100 m为半径建立缓冲区，将缓冲区内所有PS点的平均值视为该水准点对应的InSAR结果。验证结果如图6和表2所示。

(a) Track 142 (b) Track 69图6 InSAR结果与水准测量结果精度验证Fig.6 Comparison between the InSAR measurements and levelling data

表2 InSAR所得形变信息与水准测量结果比对Tab.2 Comparison of the mean subsidence rate between the Sentinel-1A PS and leveling data

其中Track 142所得研究区内垂向形变结果与水准测量结果相比，误差最大值为11.6 mm/a，最小值为0.2 mm/a，均方根误差为5.5 mm/a。对2组数据进行线性回归分析，R达0.97。Track 69所得研究区内垂向形变结果与水准测量结果相比，误差最大值为13.5 mm/a，最小值为1.0 mm/a，均方根误差为7.2 mm/a，R达0.98。验证结果表明InSAR监测结果与水准测量结果在同时间内具有较好的吻合度，表明InSAR结果的可靠性。

为进一步验证不同轨SAR影像结果的准确性，对其进行了交叉验证。由于Track 142的影像与其他2轨影像之间都具有重叠区，将其他2轨数据所得地表形变结果分别于Track 142所得结果进行比对。结果如图7所示，由Track 142与Track 69所得形变结果之间的相关性达0.96，由Track 142与Track 40所得形变结果之间的相关性为0.97。多轨SAR自相关验证结果表明了本次研究所用多轨InSAR数据结果的可靠性。

(a) Track 69与Track 142 (b) Track 40与Track 142图7 多轨SAR影像形变速率交叉验证结果Fig.7 Consistency between the vertical displacement rates derived from different datasets

4 结论

本研究基于多轨SAR数据融合技术，获取了2016—2018年京津冀平原区地表形变信息。不同轨SAR数据集交叉验证和水准验证结果均表明本次研究InSAR测量所得地面沉降信息的可靠性，满足地面沉降监测及其演化规律的研究。

从沉降场演变过程来看，2016—2018年京津冀平原区内沉降速率大于20 mm/a和50 mm/a的PS点数占总点数的比例均呈增长趋势，沉降速率大于100 mm/a的点数占总点数的比例基本维持不变(3.6%，3.5%，3.5%)。InSAR监测时段内，滨海平原范围内地面沉降发展较为迅速，其中2016—2018年，该区域地面沉降严重区(年沉降量>50 mm)面积增长195 km2，涨幅达33.8%。沉降区时序剖面结果显示，唐山沉降区呈加速发展趋势，其余沉降区发展较为稳定。InSAR所得地面沉降时序结果表明，研究时段内，京津冀平原区内地面沉降时序特征以线性下降为主，且唐山等沿海地区沉降区内地面沉降时序具有明显的季节性波动特征。

研究成果一定程度上填补了京津冀平原区地面沉降整体演化特征研究较少的空缺，为京津冀地区地面沉降防治、调控，城市地下空间安全等提供科学依据。