基于SBAS-InSAR苏通GIL综合管廊周边地表形变监测

周一鸣， 刘国林， 牛冲， 张亚凤， 郭在洁

1. 山东科技大学 测绘与空间信息学院， 山东 青岛 266590； 2. 山东省地质测绘院， 济南 250014

随着国民经济快速发展， 城市规模不断扩大， 城市地表空间资源变得格外紧张． 对地下空间的开发利用可以有效地缓解地表空间不足的问题[1]． 综合管廊是建于城市地下用于容纳两类及以上城市工程管线的构筑物及附属设施[2]． 综合管廊的施工建设会对其周边地表产生影响． 为保障综合管廊的健康长效运行， 对其进行实时、 准确的形变监测有着非常重要的意义[3]．

地表形变监测手段主要有大地水准测量、 GPS测量等传统方法， 这类测量方法工作量大， 速度慢， 易受各种自然因素的影响， 无法对大区域地表进行全面有效的监测， 很难满足时间和空间尺度的需要[4]． 近年来， 合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar， InSAR)技术广泛地应用于地表形变监测、 自然灾害监测以及地下工程监测， 其具有全天候、 覆盖范围广、 效率高、 精度高等优点[5-6]． 小基线集技术(Small Baseline Subset， SBAS)作为一种在InSAR技术基础上发展起来的重要方法， 已经成为地表形变监测的主要手段， 其可以有效减弱时空基线和大气相位的影响， 获取时间序列上的形变发展规律[7-8]． 杨振等[9]通过SBAS-InSAR技术获取到南京河西区的地表沉降监测结果， 该结果与实测水准数据具有较好的一致性， 监测结果可以为城市发展提供参考． 冉培廉等[10]利用SBAS-InSAR技术对西安市进行地表沉降监测， 全面细致地分析了西安市的地表沉降状况以及沉降原因． 范雪婷等[11]利用52景Sentinel-1A影像对南通市进行时序地表沉降监测， 利用水准数据进行精度评定， 得到了较好的结果， 为南通规划工作提供参考． ZHANG等[12]利用SBAS-InSAR技术对武汉市进行形变监测， 研究了地面各种因素与地表形变的关系， 并将监测结果与水准数据作对比， 验证了监测的可靠性． 在隧道监测方面， 杨帆等[13]以山东青岛胶州机场隧道为研究区， 采用SBAS-InSAR技术对其进行监测， 结合水准数据对比分析， 取得了较好的试验结果． 葛大庆等[14]通过高分辨率InSAR时序分析方法揭示上海10号线地铁地表形变场在建设和运营初期的分布特征、 影响范围及其随时间变化的规律， 取得较好的试验结果． HU等[15]利用SBAS-InSAR技术得到福州市区地铁建设项目后的地面沉降信息， 为城市建设和维护提供参考． 叶思远[16]利用5景COSMO-SkyMed影像对广州地铁进行时序InSAR监测， 并结合GPS监测数据， 对广州地铁沿线进行了细致、 全面的沉降分析． 李红梅等[17]利用SBAS-InSAR技术得到青岛市地铁1号线西海岸新区路段沉降信息， 结合Logistic函数模型对部分区域进行预测， 为地铁的施工及运行提供数据参考． 张严等[18]利用56景Sentinel-1A数据， 基于SBAS-InSAR技术获取佛山市禅城区地铁2号线的时空形变信息， 通过分析形变特征和地质资料， 对盾构施工中发生的塌陷事故进行了机理推测， 为今后盾构施工提供了依据．

本研究以苏通GIL综合管廊周边地表部分为研究对象， 结合其建设情况， 选取了2017年6月-2018年10月的38景Sentinel-1A数据， 利用SBAS-InSAR技术获取了苏通GIL综合管廊周边地表形变信息， 详细分析了苏通GIL综合管廊周边地表形变时空分布特征和主要成因， 并与实际水准监测数据进行对比分析和验证．

1 研究区与数据源

1.1 研究区概况

苏通GIL综合管廊工程位于江苏省南通市和苏州市交界处， 紧靠G15沈海高速苏通长江大桥(图1a)． 江南起始井(图1b)距离南岸长江大堤约250 m． 江北接收井(图1b)距离北岸长江大堤约230 m． 管廊隧道由盾构机盾构施工， 自南端起始井向北敷设， 在苏通大桥展览馆东侧避让展览馆， 下穿南岸长江大堤进入长江河道， 再下穿北岸大堤抵达江北接收井， 包括两岸工作井在内， 隧道总长度约为5 531 m， 其中盾构段约为5 469 m． 研究区地处长江河漫滩地貌单元和人工吹填地基， 地层主要有淤泥质粉质黏土、 粉砂以及粉质黏土与粉砂互层．

审图号： GS(2019)1822号

1.2 数据源

Sentinel-1卫星载有C波段合成孔径雷达， 具备多种工作模式， 是欧洲航天局哥白尼计划(Global Monitoring for Environment and Security， GMES)中的地球观测卫星． 本研究选取了覆盖研究区域2017年6月27日-2018年10月20日期间的38景Sentinel-1A的SAR影像， 对苏通GIL综合管廊周边地表形变进行监测， 数据具体参数信息如表1所示．

表1 SAR影像的具体参数

本研究采用由美国国家航空航天局(National Aeronautics and Space Administration， NASA)提供的90 m空间分辨率的数字高程模型(Digital Elevation Model， DEM)作为参考DEM数据， 去除地形相位的影响．

水准监测数据采用的是山东省鲁南地质工程勘察院在2017年7月-2018年9月的水准观测结果， 数据信息如表2所示， 其严格按照《城市轨道交通工程测量规范》(GB/T 50308-2017)技术要求进行施测， 监测结果具有较好的可靠性和真实性． 其中A，B分布于南岸地表， C，D分布于北岸地表， 每个断面共有11个水准点．

表2 水准测量数据信息

2 技术原理与数据处理流程

2.1 SBAS-InSAR技术的原理

假设有覆盖同一研究区的N+1景SAR影像， 影像的获取时间序列为t0，t1， …，tN． 选取其中一景为公共主影像， 其他N景影像作为辅影像配准到主影像上， 通过设置时空基线阈值， 生成M幅差分干涉图， 于是得到M与N之间的关系[19-20]：

(1)

如果利用tA和tB时刻获取的两景影像生成第k幅差分干涉图，tA时刻早于tB时刻， 那么生成的差分干涉图上像元(x，r)的解缠后差分干涉相位δφk(x，r)可以表示为[19-20]：

(2)

式中：λ为雷达波长；d(tB，x，r)和d(tA，x，r)分别为tB和tA两个时刻相对于t0时刻的沿雷达视线向累积形变量．

在生成干涉图的过程中， 主、 辅影像的时间序列为IE=[IE1，IE2， …，IEM]，IS=[IS1，IS2， …，ISM]， 所有的干涉对都满足IEk>ISk． 那么， 干涉相位可由下式表示：

δφk=φ(tIEk)-φ(tISk)

(3)

其矩阵形式可表示为[19-20]：

Aφ=δφ

(4)

系数矩阵A为一个M×N的矩阵， 当处理所用的影像数据属于一个共同的小基线子集时，M≥N， 此时矩阵A的秩为N． 由最小二乘法可得其唯一解φ的估值为[19-20]：

(5)

若所处理的影像数据属于L个不同的基线子集， 那么方程(5)是秩亏的， 其秩亏数为N-L+1， 需要采用奇异值分解(Singular Value Decomposition， SVD)方法求解， 然后再得出φ的最小二乘最小范数解．

2.2 数据处理流程

通过对38景数据进行基线估算， 确定2018年3月6日的SAR影像为公共主影像， 设置时间基线阈值为55 d， 空间基线阈值设为200 m， 共生成了89组干涉对． 利用DEM数据对生成的干涉对进行差分干涉处理， 并采用Goldstein滤波方法对差分干涉图进行滤波处理， 采用最小费用流(Minimum Cost Flow， MCF)方法对滤波后干涉图进行解缠处理， 设置解缠的相干系数阈值为0.25． 然后进行轨道精炼与重去平， 在远离形变区域的地方选取了25个地面控制点(Ground Control Points， GCP)， 用于估算和去除残余的误差相位的影响． 为得到较好的处理结果， 进行估算形变速率和残余地形处理及二次解缠和优化工作， 在此基础上进行大气滤波处理， 进一步去除大气延迟相位， 获取时间序列上的位移结果． 最后对其进行地理编码， 得到该研究区地表的最终形变结果． 对所得结果进行地表形变分析和水准精度验证．

数据处理的流程如图2所示．

图2 数据处理流程

3 结果与分析

3.1 研究区地表形变监测结果

本次实验运用SBAS-InSAR技术对苏通GIL综合管廊周边地表进行形变监测， 获得了研究区在2017年6月-2018年10月期间内的地表形变信息， 图3给出了基于SBAS-InSAR监测得到的苏通GIL综合管廊周边2017年6月-2018年10月的地表形变速率， 其中图3a显示了苏通GIL综合管廊南岸地表形变分布及其形变速率； 图3b显示了苏通GIL综合管廊周边地表南岸和北岸的地表形变分布及其形变速率； 图3c显示了苏通GIL综合管廊北岸地表形变分布及其形变速率． 图4是以2017年6月-2018年10月期间， 4个时间段内南、 北岸累计地表形变情况．

图3 苏通GIL综合管廊周边地表形变速率

图4 苏通GIL综合管廊周边4个时间段内南、 北岸累计地表形变

3.2 研究区时空演化特征分析

由图3、 图4可以看出：

在苏通GIL综合管廊的南岸地表区域， 如图3a所示， 管廊轴线地表形变呈带状分布， 主要的形变速率范围为0～5 mm/a， 在江南起始井的南部和苏通大桥展览馆附近等局部区域形变速率范围为5～15 mm/a． 在南岸管廊轴线区域这一位置， 盾构机处于工作阶段， 盾构机掘进姿势的调整和盾构施工中注浆压力的试验调整等因素可能是该区域地表形变的主要原因． 在管廊的轴线外围西侧， 地表形变速率范围为-15～0 mm/a， 在管廊的轴线外围东侧， 地表形变呈现沉降漏斗状分布， 地表形变速率范围为-40～0 mm/a， 由于该地区为长江河漫滩地貌， 经过地表人类活动， 使得地表土层不稳定产生地表形变．

在苏通GIL综合管廊的北岸地表区域， 如图3c所示， 江北接收井以及管廊轴线地表形变以沉降漏斗状分布， 最大形变速率在D断面达到-49 mm/a． 北岸地质为人工吹填地， 出现沉降漏斗与不稳定的软土地质有关， 同时也与江北接收井排水以及盾构机扰动土层下沉有关． 北岸东南角靠江沿岸出现较大形变， 形变速率范围为-40～-15 mm/a， 这与该处靠江建筑有关．

2017年7月9日， 如图4a所示， 北岸地面相对比较稳定， 沉降量整体比较小， 未出现比较严重的沉降现象， 沉降量范围主要为-5～0 mm， 靠江边地区沉降量范围为-20～-10 mm． 2017年11月6日， 如图4b所示， 北岸接收井以及北部附近区域出现不均匀沉降， 其沉降量范围为-30～-10 mm， 这可能与北岸接收井地面施工以及人类活动影响有关． 在图4b东南角也出现了较为严重的沉降， 其沉降量为-50～-20 mm， 这与江边一建筑有关． 2018年5月5日， 如图4c所示， 盾构机即将跨过长江开始在北岸进行盾构工作， 此时地表的施工建设重心转移至北岸， 北岸地表沉降现象加重， 沉降加重区集中在北岸接收井附近， 沉降量范围为-30～0 mm． 在北岸管廊轴线上出现严重沉降， 最大沉降量达到-38 mm， 且东南角沉降仍在继续， 沉降加剧． 2018年10月20日， 如图4d所示， 盾构隧道已经贯通， 北岸沉降量进一步加重并向着北岸管廊轴线扩展， 沉降以沉降漏斗形式分布， 沉降量范围为-91～-30 mm， 从整个北岸时空演化过程来看， 北岸一直处在地表沉降的状态， 地面沉降主要分布在北岸接收井附近区域， 且随着工程的发展， 沉降加重并向江北接收井四周发展．

2017年7月9日， 如图4e所示， 南岸开始管廊的施工建设， 南岸整体出现沉降， 大部分地区沉降量为-5～0 mm， 个别区域沉降量为-10～-5mm． 2017年11月6日， 如图4f所示， 南岸盾构工作完成， 盾构机进入江底工作， 管廊轴线区域地表整体继续出现沉降现象， 沉降量范围为-5～0 mm． 在管廊轴线西侧靠江区域出现沉降， 最大沉降量达到-66 mm， 东侧规划用地沉降量范围主要为-20～-5 mm， 东侧规划用地局部形变量范围为0～5 mm， 该地表形变由人类活动和盾构后软土地质不稳定造成． 2018年5月5日， 在图4g中， 南岸管廊施工已经结束数月， 南岸整体沉降量范围为-20～-5 mm． 西北部靠江沿岸沉降进一步加重， 沉降量范围为-75～-30 mm． 2018年10月20日， 如图4h所示， 南岸起始井南部地表抬升， 形变量为0～5 mm． 苏通大桥展览馆区域出现抬升， 形变量范围为0～5 mm． 管廊轴线其余地区沉降量范围为-5～0 mm， 管廊西侧沉降量范围为-30～-10 mm， 东侧规划用地沉降量范围为-50～-10 mm， 出现沉降漏斗． 从整个南岸地表形变的时空演化过程来看， 南岸整体地表形变经历了先沉降后抬升的过程， 南岸管廊沿线区域未发生较为严重的不均匀地表形变， 地表形变最终相对稳定．

4 验证分析

为了进一步验证SBAS-InSAR技术的精确性和可靠程度， 在图3a和图3c中共选取了苏通大桥展览馆附近的10个水准点监测数据和垂直于隧道轴线方向的A，B，C，D 4个横向监测断面上44个水准点监测数据． 其中， A，B分布于南岸地表， C，D分布于北岸地表， 每个断面共有11个水准点， 水准点以隧道轴线点为中心， 两侧各间隔3，5，10，10，10 m分布． A，B断面点号由东到西为A1-A11和B1-B11， C，D断面点号由西到东为C1-C11和D1-D11， 苏通大桥展览馆点号从西南角点开始逆时针依次为J1-J10． 将水准点的水准监测数据与SBAS-InSAR监测最临近时间段内监测数据分别作为x和y绘制在直角坐标系中， 具体时间段如表3所示， 根据其与蓝色虚线y=x的离散程度就可以对其进行精度验证， 如图5所示．

表3 监测时段信息

图5 SBAS-InSAR监测结果与水准监测结果对比

由图5可以看出， 绝大多数点的坐标都分布于直线y=x附近， 只有少数点距离直线y=x较远， 说明x和y的数值较为接近， SBAS-InSAR监测的形变量与水准监测的形变量大小较为一致． 两种方法监测结果的拟合优度R2为0.819， 均方根误差(Root Mean Square Error， RMSE)约为3.597， 最大偏差绝对值(Maximum Difference， MaxD)为12 mm， 最小偏差绝对值(Minimum Difference， MinD)为1 mm． 两种监测结果基本相同， 说明了SBAS-InSAR对苏通GIL综合管廊周边地表以及建筑物的监测精度结果较好．

为了更加全面地验证SBAS-InSAR监测的准确性， 分别选取了南岸B断面、 北岸D断面和苏通大桥展览馆各3个点， 将水准监测时段内所得的数据与SBAS-InSAR在该时段内监测的数据进行对比， 如图6所示． 同时， 选取断面B，D和苏通大桥展览馆在水准监测最后一天的结果进行空间上的对比， 如图7所示． 由于SBAS-InSAR监测数据时间间隔要比水准监测时间间隔大， 因此对SBAS-InSAR监测数据进行了线性插值处理．

图6 9个特征点的时间序列对比

图7 断面B，D和苏通大桥展览馆空间对比

绝大多数点的SBAS-InSAR监测数据值与水准监测数据值非常接近， 如图6a，6c，6f和6g，6h，6i， 但在如图6b，6d和6e中， B6，D1和D6这种管廊正上方位置和北岸盾构后形变大的地方， 两种监测数据值出现差别． 在时间尺度上， 大部分点的SBAS-InSAR监测值与水准监测值在时间上形变趋势一致， 但在短时间内水准变化大的时间段， SBAS-InSAR监测只能反映大致的形变变化． 由图7可知， SBAS-InSAR监测的数据值与水准监测的数据值比较吻合， 在空间尺度上， 形变趋势保持一致．

5 结论

本研究利用2017年6月29日-2018年10月20日覆盖苏通GIL综合管廊周边地表的38景Sentinel-1A哨兵数据和54个水准点观测数据， 系统研究和分析了基于SBAS-InSAR苏通GIL综合管廊周边地表形变监测． 主要结论有：

SBAS-InSAR监测得到苏通GIL综合管廊周边地表出现了不同程度的地表形变． 南岸地表整体比较稳定， 只有在轴线外围东侧出现了明显形变； 北岸地表发生了较为明显的地表形变， 主要的形变区域集中在管廊轴线位置并且出现了沉降漏斗． 在整个监测期间内， 南岸整体地表形变经历了先沉降后抬升的过程， 北岸一直处在地表沉降的状态， 并且随着盾构施工， 沉降加重并且沉降范围进一步扩大．

SBAS-InSAR监测值与水准监测值具有良好的一致性， 在空间尺度上， 二者监测到的地表形变趋势一致． 在时间尺度上， 在形变较小的水准点上， SBAS-InSAR监测值与水准监测值在时间上形变趋势一致， 但形变较严重的水准点上， SBAS-InSAR监测只能反映大致的形变变化， 两种方法监测得到的形变量值出现较大差别．

总体而言， 利用SBAS-InSAR技术能够有效地监测和分析地下综合管廊地表及周边区域大范围、 长时间序列的形变情况． 本研究有助于推动InSAR在地下综合管廊、 隧道等地下工程形变监测中的应用， 为地下工程的施工和运营提供保障．