在城市地铁施工区域开展地面塌陷的监测和早期预警研究具有重要意义。在地面塌陷灾害形成过程中，地表早期会出现缓慢下沉、开裂等异常现象［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所示。
2.2 SBAS-InSAR数据处理方法
SBAS-InSAR是一种新开发的星载合成孔径雷达差分干涉测量技术（differential intereferometric synthetic aperture radar，D-InSAR）时间序列分析方法，它将采集到的所有SAR数据组合成干涉对，然后设置空间基线和时间基线阈值，并获得几个时空基线小于阈值的小基线集。每个集合的表面变形时间序列可以通过奇异值分解（singular value decomposition，SVD）求解集合来解决。这有效地解决了由于基线时间长不连续性概率引起的基线组间差异小基线，提高了监测时间分辨率，有效解决了一般方程秩缺陷问题，最终得到覆盖整个地面沉降观测时间的序列，SBAS-InSAR技术的数据处理流程如图2所示。
<G:\武汉工程大学\2023\第6期\刘金昌-2.tif>
图2 SBAS-InSAR技术流程图
Fig. 2 SBAS-InSAR technical flowchart
（1）将监测区域的M+1景SAR影像数据集根据给定的干涉组合条件，生成N景小基线距的干涉图，如图3（a）所示，将干涉图中的地形相位去除生成去地形干涉图。
（2）以相干系数为基础生成掩膜图，如图3（b）所示，在去地形干涉图中选择高相干点进行相位解缠处理，如图3（c）所示，并对相干点进行定标处理。
数据处理流程中生成的干涉图、去除轨道倾斜改正、相位解缠、掩膜如图3所示。
<G:\武汉工程大学\2023\第6期\刘金昌-3.tif>[0 876 1 752 2 628
像素值 / mm][0
528
1 057
1 585][像素值 / mm][累计形变值 / mm][0 876 1 752 2 628
像素值 / mm][45
30
15
0
-15
-30
-45
-60][累计形变值 / mm][1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0][累计形变值 / mm][45
30
15
0
-15
-30
-45
-60][0
528
1 057
1 585][像素值 / mm][0 876 1 752 2 628
像素值 / mm][0 876 1 752 2 628
像素值 / mm][45
30
15
0
-15
-30
-45
-60][（b）][（a）][（d）][（c）][累计形变值 / mm]
图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
（3）根据不同时段的平均形变速率、高程误差、差分相位进行模型方程组的建立，利用SVD算法获得最小范数条件下的最小二乘解。在开展SAR数据时序分析前，使用斜率格网法估计轨道效应，进行轨道误差校正，如图3（d）去除轨道倾斜所示。使用最小二乘法在每个干涉图上独立估计轨道误差。默认情况下，对于由指数为[i]和[j]的2个SAR采集组成的干涉图，采用最小范数反演算法解决，[dij(x,y)-Rij(x,y)2]中[dij(x,y)]的值是方位[x]，[y]处的像素值。
[Rijx,y=eij?xy+aij?x+bij?y+cij] （1）
其中[aij]、[bij]、[cij]和[eij]称为干涉图[ij]的轨道参数。
当像素格网为1时，
[Rijx,y=cij] （2）
当像素格网为3时，
[Rijx,y=aij?x+bij?y+cij] （3）
当像素格网为4时，
[Rijx,y=eij?xy+aij?x+bij?y+cij] （4）
为了确保干涉测量网络的一致性，在网络感知中重新估计轨道参数。通过i反转线性系统来估计每个SAR采集的轨道参数。
[aij=ai-aj] （5）
[bij=bi-bj] （6）
[cij=ci-cj] （7）
[dij=di-dj] （8）
结合重新估计的轨道参数，生成与干涉测量网络一致的轨道校正图，并校正每个干涉图。
假设使用覆盖同一研究区域的N+1幅SAR影像图，则按时间顺序对它们进行排序：
[t=t0, t1, t2, ?, tM-1T] （9）
N个干涉对是通过最大空间基线阈值获得的M对多视差分干涉对生成的，M满足：
[M2≤n≤MM-12] （10）
对生成的多视差分干涉图进行滤波并解开相位，使用相干性在未包覆的差分干涉图中选择适当的参考点。对于任意2个微分干涉图（[ta<tb]），未包覆的相位可以表示为：
[?x,i=?x,itb-?x,ita =?flax,x,i+?ele,x,i+]
[?disp,x,i+?atm,x,i+?noise] （11）
式（11）中：[?x,i]是i图像和x像素的干涉相位，[?x,itb]和[?x,ita]是x像素2个不同时间的干涉相位。[?flat,x,i]是平坦相位的相位差，[?ele,x,i]是地形相位的相位差，[?disp,x,i]是地面变形的相位差，[?atm,x,i]是大气相位的误差，[?noise]是噪声的误差。
公式中的平坦相位、地形相位、变形相位、大气相位可以进一步表示为：
[?flat,x,i=-4πλBnRtanθ] （12）
[?ele,x,i=-ΔqsinθBnR4πλ] （13）
[?disp,x,i=4πλdtb,x,i-dta,x,i] （14）
[?atm,x,i=?atmtb,x,i-?atmta,x,i] （15）
式（12）~式（13）中：λ是雷达波长，R是斜率距离，Bn是垂直基线，Dq是DEM高程差，q是入射角。
假设2个不同干涉图之间的变形速率为[Vii=1, 2, 3,?, n]，则ta和tb之间的累积变形可以表示为：
[?disp,x,i=4πλi=tatbVi, j+1ti+1-ti] （16）
变形速率的积分在干涉图的主图像和从图像之间的时间间隔内。
直接在主图像和从图像间隔之间的每个时间间隔处进行速率积分。以矩阵形式书写：
[BV=δ?] （17）
通过使用SVD方法，可以获得系数矩阵B的广义逆矩阵，从而得到数据向量的最小范数解。最后，每个时间段的速度决定了每个时间段的变形量。
（4）开展估算大气相位与非线性形变。对于给定数据集，选择与SAR场景的空间覆盖范围重叠的格网点。大气变量以精确的压力水平提供。将这些值垂直插值到表面和参考高程之间的规则网格中，计算每个选定网格点上的延迟函数作为高程的函数，高程处的视线方向（line of sight，LOS）单路径延迟[δsLOSz]［17］。
[δsLOSz=10-6cosθk1RdgmP?-Pzref+]
[?Zrefk2-RdRvk1eT+k3eT2dz] （18）
式（18）中：q为入射角；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×103 K2·Pa-1。
3 结果与讨论
根据上述数据处理流程，使用改进的SBAS-InSAR时序分析方法获取了青岛市范围内2018年8月9日至2020年8月22日期间的地表累计形变计算结果，将结果叠加在天地图影像上，如图4所示。
以事故点为中心，选取方圆300 m范围内的时序形变点开展数理统计分析。分别在2个事故点提取分析多个点位的时序形变信息，如图5所示。
将InSAR数据获取的地表变形量数据转化为日均形变速率，处理后得到该点的日均形变速率，单位为mm/d，如图6所示。形变速率计算与分析过程如下：
（1）将时间信息和沉降信息一一对应：D0-S0，D1-S1，D2-S2，D3-S3，…，Dn-Sn；
（2）计算日均地表形变速率
[Vn=Sn-Sn-1Dn-1-Dn] （19）
式（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数据，基于改进的SBAS-InSAR方法，计算长时序、高精度地表形变数据，通过与青岛市2次事故对比分析，提出了地铁施工塌陷监测早期识别方法和监测预警指标。主要结论如下：
（1）本文中通过改进的SBAS-InSAR方法，利用SVD算法获得最小范数条件下的最小二乘解，使用斜率格网法估计轨道效应，并进行轨道误差校正，同时结合高程、气温、气压及湿度估计值构建估算函数，使用大气变量的3D分布来确定每个干涉图的每个像素上的大气相位延迟，以提高SBAS-InSAR方法处理的地表形变数据精度。实验表明该方法获取的地表监测数据能够满足对城市地铁施工监测的需求。
（2）通过统计分析方法对沙子口和胜利桥两处事故点长时序地表形变监测数据进行对比分析，将长时序累积形变量转换为日均沉降速率，提出地表沉降监测早期识别指标为单周期日均沉陷速率大于0.02 mm/d，施工预警阈值为连续3个观测周期（36 d）大于0.02 mm/d。
因此，基于改进的SBAS-InSAR获取的高精度、长时序地表形变监测信息能够为地铁施工建设中的塌陷早期识别和监测预警提供辅助决策依据，本研究结果可为在建地铁工程及未来规划线路的建设提供大范围监测的理论技术支撑。