地铁施工塌陷的SBAS-InSAR长时序监测与早期识别

刘金昌，王慧妮，徐 吟，李康伦，张华睿，刘天鹤

1. 武汉工程大学土木工程与建筑学院，湖北 武汉 430074；

2. 中铁十一局集团有限公司，湖北 武汉 430061；

3. 武汉综合交通研究院有限公司，湖北 武汉 430015

在大规模的地铁施工建设过程中，由于受到地质环境、周边环境、施工方法和方案等因素的影响，由城市地铁施工引起的地面塌陷、突涌水［1］等灾害屡见不鲜。地面塌陷是指在人为因素或者自然因素的作用下，地表岩、土体向下陷落，并形成塌陷坑（洞）的一种动力地质现象［2］。地面塌陷是突发性的，尤其是在城市建筑物密集的地方，一旦发生塌陷事故将会造成特别严重的经济损失和生命安全问题。以青岛地铁施工为例，2019 年5 月27 日，青岛地铁4 号线沙子口静沙区间施工段发生坍塌，同年7 月4 日，青岛地铁1 号线胜利桥站施工围挡处再次发生塌陷。

在城市地铁施工区域开展地面塌陷的监测和早期预警研究具有重要意义。在地面塌陷灾害形成过程中，地表早期会出现缓慢下沉、开裂等异常现象［3］，可以利用塌陷灾变演化过程特性基于长时序大范围地表形变开展监测，以实现潜在的地面塌陷灾害隐患点的早期识别及预警。合成孔径雷达干涉测量（interferometric synthetic aperture radar，InSAR）技术是利用合成孔径雷达（synthetic aperture radar，SAR）对同一地区观测到的2 幅复数值影像（既有幅值又有相位的影像）数据进行相干处理，以获取地表高程信息的新技术，因其具有全天候、广覆盖、高精度、低成本、低风险获取的特点［4-7］，现如今已经广泛应用于冰川融化、地质灾害监测、城市地下水变化和地面塌陷等沉降的监测中［8-11］。InSAR 技术的测量范围大，形变测量精度可达到mm 级别，非常适合对城市地表的形变监测研究。目前，国内外学者利用InSAR 技术开展的地面沉降方面研究工作主要集中于地铁沿线的沉降监测［12-16］，但基于长时序地表形变监测数据结合历史塌陷事故，开展地面塌陷预警的研究工作相对较少。

本文以青岛地铁工程施工过程中的2 次地铁塌陷事故为研究对象，收集2018 年8 月至2020 年8 月期间青岛市的60 景哨兵1 号（Sentinel-1A）卫星数据，基于差分干涉测量短基线集时序分析（small baseline subset-InSAR，SBAS-InSAR）技术获取塌陷事故前后期地表的时序形变过程数据，结合历史塌陷事故开展地表时序形变与塌陷灾害演化过程相关性分析，建立利用SAR 卫星遥感数据开展城市地铁施工地面塌陷早期识别与监测预警方法，为城市地铁施工监测预警提供技术支撑。

1 基本情况

青岛地处山东半岛东南部，位于东经119°30′~121°00′、北纬35°35′~37°09′。青岛所处大地构造位置为新华夏隆起带次级构造单元——胶南隆起区东北缘和胶莱凹陷区中南部。区内缺失整个古生界地层及部分中生界地层，但白垩系青山组火山岩层发育充分，在青岛市出露十分广泛。岩浆岩以元古代胶南期月季山式片麻状花岗岩及中生代燕山晚期的艾山式花岗闪长岩和崂山式花岗岩为主。市区全部坐落于该类花岗岩之上，建筑地基条件优良。构造以断裂构造为主。

沙子口事故发生于2019 年5 月27 日，青岛地铁4 号线沙子口静沙区间施工段发生坍塌。胜利桥事故发生于2019 年7 月4 日，青岛地铁1 号线胜利桥站施工围挡处发生直径约10 m 的塌陷。

2 实验数据和方法

2.1 SAR 数据

本研究选择了地球观测卫星Sentinel-1 的遥感数据。该卫星配备了C 波段SAR，由2 颗卫星A，B 组成，卫星重访期为12 d。本实验搜集了覆盖研究区2018-08－2020-08 期间的60 景Sentinel-1A 数据，所选数据极化方式均为VV（垂直极化发射，垂直极化接收）极化，所处波段均为C 波段，轨道方向为升轨，方位向和距离向分辨率分别为20 m 和5 m。研究区域使用的外部数字高程模型（digital elevation model，DEM）数据为德国航空航天中心使用TerraSAR-X/TanDEM-X 卫星获取的全球陆地高程数据，空间分辨率为30 m。为了提高Sentinel-1 影像的轨道精度，采用了欧空局发布的精确轨道数据，将研究区域数据按照小基线原则（只保留平均相干性大于0.4 的干涉对，时间基线设置为12 d/24 d，空间基线由SBAS-InSAR 方法进行计算）组合而成干涉对的时空基线图，如图1 所示。

图1 时空基线图Fig.1 Space-time baseline map

2.2 SBAS-InSAR 数据处理方法

SBAS-InSAR 是一种新开发的星载合成孔径雷达差分干涉测量技术（differential intereferometric synthetic aperture radar，D-InSAR）时间序列分析方法，它将采集到的所有SAR 数据组合成干涉对，然后设置空间基线和时间基线阈值，并获得几个时空基线小于阈值的小基线集。每个集合的表面变形时间序列可以通过奇异值分解（singular value decomposition，SVD）求解集合来解决。这有效地解决了由于基线时间长不连续性概率引起的基线组间差异小基线，提高了监测时间分辨率，有效解决了一般方程秩缺陷问题，最终得到覆盖整个地面沉降观测时间的序列，SBAS-InSAR 技术的数据处理流程如图2 所示。

图2 SBAS-InSAR 技术流程图Fig.2 SBAS-InSAR technical flowchart

（1）将监测区域的M+1 景SAR 影像数据集根据给定的干涉组合条件，生成N景小基线距的干涉图，如图3（a）所示，将干涉图中的地形相位去除生成去地形干涉图。

图3 SBAS-InSAR 干涉及相位解缠过程：（a）干涉，（b）掩膜，（c）相位解缠，（d）去除轨道倾斜Fig.3 SBAS-InSAR interferogram and phase unwrapping process：（a）interference，（b）mask，（c）phase unwrapping，（d）removal of orbital tilt

（2）以相干系数为基础生成掩膜图，如图3（b）所示，在去地形干涉图中选择高相干点进行相位解缠处理，如图3（c）所示，并对相干点进行定标处理。

数据处理流程中生成的干涉图、去除轨道倾斜改正、相位解缠、掩膜如图3 所示。

（3）根据不同时段的平均形变速率、高程误差、差分相位进行模型方程组的建立，利用SVD 算法获得最小范数条件下的最小二乘解。在开展SAR 数据时序分析前，使用斜率格网法估计轨道效应，进行轨道误差校正，如图3（d）去除轨道倾斜所示。使用最小二乘法在每个干涉图上独立估计轨道误差。默认情况下，对于由指数为i和j的2 个SAR 采集组成的干涉图，采用最小范数反演算法解决中dij(x,y)的值是方位x，y处的像素值。

其中aij、bij、cij和eij称为干涉图ij的轨道参数。

当像素格网为1 时，

当像素格网为3 时，

当像素格网为4 时，

为了确保干涉测量网络的一致性，在网络感知中重新估计轨道参数。通过i反转线性系统来估计每个SAR 采集的轨道参数。

结合重新估计的轨道参数，生成与干涉测量网络一致的轨道校正图，并校正每个干涉图。

假设使用覆盖同一研究区域的N+1 幅SAR 影像图，则按时间顺序对它们进行排序：

N个干涉对是通过最大空间基线阈值获得的M对多视差分干涉对生成的，M满足：

对生成的多视差分干涉图进行滤波并解开相位，使用相干性在未包覆的差分干涉图中选择适当的参考点。对于任意2 个微分干涉图（ta

式（11）中：ϕx,i是i图像和x像素的干涉相位，ϕx,i(tb)和ϕx,i(ta)是x像素2 个不同时间的干涉相位。ϕflat,x,i是平坦相位的相位差，ϕele,x,i是地形相位的相位差，ϕdisp,x,i是地面变形的相位差，ϕatm,x,i是大气相位的误差，ϕnoise是噪声的误差。

公式中的平坦相位、地形相位、变形相位、大气相位可以进一步表示为：

式（12）~式（13）中：λ是雷达波长，R是斜率距离，Bn是垂直基线，Δq是DEM 高程差，θ是入射角。

假设2 个不同干涉图之间的变形速率为Vi(i= 1,2,3,…,n)，则ta和tb之间的累积变形可以表示为：

变形速率的积分在干涉图的主图像和从图像之间的时间间隔内。

直接在主图像和从图像间隔之间的每个时间间隔处进行速率积分。以矩阵形式书写：

通过使用SVD方法，可以获得系数矩阵B的广义逆矩阵，从而得到数据向量的最小范数解。最后，每个时间段的速度决定了每个时间段的变形量。

（4）开展估算大气相位与非线性形变。对于给定数据集，选择与SAR 场景的空间覆盖范围重叠的格网点。大气变量以精确的压力水平提供。将这些值垂直插值到表面和参考高程之间的规则网格中，计算每个选定网格点上的延迟函数作为高程的函数，高程处的视线方向（line of sight，LOS）单路径延迟

式（18）中：θ为入射角；Rd＝287.05 J·kg-1·K-1，Rν＝461.495 J·kg-1·K-1，分别是干燥和湿润环境下的气体常数；gm是z和zref之间重力加速度的加权平均值；P是干燥空气分压（Pa）；e是水蒸气分压（Pa）；T是温度（K）。式（18）中：k1＝0.776 K·Pa-1，k2＝0.716 K·Pa-1，k3＝3.75×103K2·Pa-1。

3 结果与讨论

根据上述数据处理流程，使用改进的SBASInSAR 时序分析方法获取了青岛市范围内2018 年8 月9 日至2020 年8 月22 日期间的地表累计形变计算结果，将结果叠加在天地图影像上，如图4所示。

图4 地表累计形变图：（a）沙子口事故区域，（b）胜利桥事故区域Fig.4 Cumulative surface deformation maps：（a）Shazikou accident area，（b）Shengliqiao accident area

以事故点为中心，选取方圆300 m 范围内的时序形变点开展数理统计分析。分别在2 个事故点提取分析多个点位的时序形变信息，如图5 所示。

图5 事故区域部分点位累计沉降：（a）沙子口，（b）胜利桥Fig.5 Cumulative subsidence of some points in accident areas：（a）Shazikou，（b）Shengliqiao

将InSAR 数据获取的地表变形量数据转化为日均形变速率，处理后得到该点的日均形变速率，单位为mm/d，如图6 所示。形变速率计算与分析过程如下：

图6 沙子口和胜利桥研究区域日均沉降速率Fig.6 Average daily subsidence rates of Shazikou and Shengliqiao study areas

（1）将时间信息和沉降信息一一对应：D0-S0，D1-S1，D2-S2，D3-S3，…，Dn-Sn；

（2）计算日均地表形变速率

式（19）中：Di（i=1，2，3，…，n）是日期，Si（i=1，2，3…，n）是沉降量。因为用的是前1 个周期的日期减去后1 个周期的日期，所以地面发生沉降时，沉降速率为正值。

基于SBAS-InSAR 时序形变观测信息对沙子口研究区开展时序分析，如图6 沙子口事故所示，2018 年8 月21 日，该区域便开始呈现下沉趋势，下沉速率在0.011 4 mm/d 左右，有地表沉降形变但速率相对较为缓慢。2018 年8 月21 日至2019 年1 月12 日该研究区沉降速率呈现加速增长趋势，于2019 年1 月12 日沉降速率突破0.02 mm/d，在这一段时间内该研究区累计沉降量为2.37 mm。从2019 年1 月12 日至2019 年5 月27 日发生地面塌陷事故，在连续132 d 内的11 次观测（Sentinel-1 卫星的重访周期为12 d）中，研究区日均沉降速率均大于0.02 mm/d，期间地面累计沉降量为3.29 mm。在事故发生后，地表形变沉降速率显著下降，地表形变逐渐趋于稳定。

在胜利桥区域，如图6 胜利桥事故所示，从2018 年8 月21 日开始，该区域沉降速率处于-0.010 5 mm/d 左右，这是由于青岛市所处板块的地壳运动，呈现缓慢的上升趋势，但是到2018 年10 月20 日开始出现沉降，沉降速率开始快速增大。出现沉降后，该研究点的沉降速率越来越快，于2019 年5 月12 日沉降速率突破0.02 mm/d，此时累计沉降量为1.41 mm。从2019 年5 月12 日到2019 年7 月4 日发生地面塌陷事故，连续4 次卫星观测（53 d）的日均沉降速率均大于0.02 mm/d，期间地面累计沉降量为1.44 mm。胜利桥区域在日均沉降速率最大（0.024 5 mm/d）处发生塌陷事故，事故发生后，该区域沉降速率逐渐降低，于2020 年1 月7 日沉降速率开始小于0 mm/d，期间累计沉降量为4.88 mm。最后于2020 年3 月开始趋于稳定的上升状态。

实验数据说明：日均沉降速率大于0.02 mm/d可作为监测预警依据，若在施工期间，连续在3 个监测周期（36 d）出现日均沉降速率大于0.02 mm/d时，应对施工单位发出塌陷事故预警，提醒施工单位此地存在塌陷灾害隐患，及时采取更加精密的监测手段来监测地面沉降情况并作出事故预防措施。

4 结 论

通过获取2018 年8 月—2020 年8 月期间青岛市的60 景Sentinel-1A 数据，基于改进的SBASInSAR 方法，计算长时序、高精度地表形变数据，通过与青岛市2 次事故对比分析，提出了地铁施工塌陷监测早期识别方法和监测预警指标。主要结论如下：

（1）本文中通过改进的SBAS-InSAR 方法，利用SVD 算法获得最小范数条件下的最小二乘解，使用斜率格网法估计轨道效应，并进行轨道误差校正，同时结合高程、气温、气压及湿度估计值构建估算函数，使用大气变量的3D 分布来确定每个干涉图的每个像素上的大气相位延迟，以提高SBAS-InSAR 方法处理的地表形变数据精度。实验表明该方法获取的地表监测数据能够满足对城市地铁施工监测的需求。

（2）通过统计分析方法对沙子口和胜利桥两处事故点长时序地表形变监测数据进行对比分析，将长时序累积形变量转换为日均沉降速率，提出地表沉降监测早期识别指标为单周期日均沉陷速率大于0.02 mm/d，施工预警阈值为连续3 个观测周期（36 d）大于0.02 mm/d。

因此，基于改进的SBAS-InSAR 获取的高精度、长时序地表形变监测信息能够为地铁施工建设中的塌陷早期识别和监测预警提供辅助决策依据，本研究结果可为在建地铁工程及未来规划线路的建设提供大范围监测的理论技术支撑。