多时相Sentinel-1A InSAR的连盐高铁沉降监测分析

何秀凤，高 壮，肖儒雅，罗海滨，冯 灿

1. 河海大学地球科学与工程学院,江苏 南京 211100； 2. 南京信息工程大学遥感与测绘工程学院，江苏 南京 210044； 3. 北方信息控制研究院集团有限公司，江苏 南京 211153

近年来，随着我国经济发展和城市化进程的不断加快，我国高铁事业得到了快速发展。2013—2019年是我国铁路投资最集中、强度最大的时期，截至2017年年末，我国高铁运营里程为25 000 km，占世界高铁里程总量的66%[1]。虽然高铁为人口产业集聚、要素流动、时空压缩和地方高铁新城的涌现提供了一系列有利条件，但是在建造长距离轨道线路的过程中，脆弱的工程地质环境、纵横交错的地表构筑物载荷和地下水的过度开采等自然和人为因素，极易引发区域地表沉降和不均匀的纵向变形，给高铁安全运营带来了极大的隐患，甚至造成严重的经济损失和不良社会影响。高铁线路通常呈线性走向分布，具有距离长、跨度大、分布不均等特点[2]。在设计、施工和运营监测过程中，线路需要进行形变监测。而常规单一的形变测量手段，如水准测量、GNSS测量[3]和三维激光扫描，主要依赖监测站点的重复观测，对观测环境条件要求高，需要耗费大量的人力物力，难以监测高铁沿线大跨度、长时间的沉降趋势情况，从而给高铁的运营和维护带来极大的困难。

合成孔径雷达干涉测量(interferometry synthetic aperture radar，InSAR)作为近些年发展迅猛的大地测量手段，具有监测范围大、时空分辨率高，可利用数据源多，可以全天时全天候观测的优点，在形变监测方面具有独特的优势。文献[4—5]利用D-InSAR(differential InSAR)手段进行地表沉降监测，并融合GPS技术获取三维形变场信息。但是常规D-InSAR方法易受时空去相关和大气延迟误差的影响，形变监测的精度和可靠性受到限制。为此，多时相InSAR(multiple temporal InSAR，MT-InSAR)技术在D-InSAR的基础上逐渐发展起来，如永久散射体干涉测量(permanent scatterer interferometry，PSI)技术[6-7]和小基线集(small baseline subset，SBAS)技术[8-9]。它克服了常规D-InSAR中的时空去相关和大气效应问题，在监测人工线状地物形变上具有独特的优势，展现出极具潜力的应用前景。文献[10]利用26景TerraSAR-X数据，通过PS-InSAR技术获得了上海轨道交通的形变结果；文献[11]利用相同的数据源和处理方法对上海宝山区1、3号地铁线进行监测，揭示了地铁沿线沉降的时空格局；文献[12]利用天津地区某条铁路沿线的水准数据对时序InSAR监测结果进行验证，二者结果取得很好的一致性；文献[13]利用相同的方法对京津线进行监测，获取了京津线在运营过程中的沉降情况；文献[14]对青海省内G214高速公路进行PS-InSAR监测，并探讨了冻土季节性解冻现象对公路沉降的影响；文献[15]利用X波段的COSMO-SkyMed数据，联合PS-InSAR与DS-InSAR技术对北京至石家庄的高铁展开监测，验证了COSMO数据能够获取亚毫米级线性形变速率。

应用多时相InSAR技术监测人工线状地物(如铁路、高速公路)虽取得了一些研究成果，但大多都是宏观的对监测结果进行解释，加以定性的分析，对InSAR监测结果进行定量评价的较少，同时就应用的SAR卫星数据源来说大都是X波段的数据，应用C波段数据的案例鲜少。因此，本文以刚开始运营的连云港至盐城的高铁作为研究对象，获取研究区2018年1月—2019年9月期间共47景C波段Sentinel-1A数据，利用MT-InSAR技术进行处理分析。此外，采用同时期内的BDS观测数据进行对比验证，探究应用C波段SAR数据的多时相InSAR技术在高铁形变监测的应用能力。由结果可知，利用C波段的Sentinel-1A数据能够获取毫米级的线性形变速率和时序位移序列，与BDS监测结果相比取得很好的一致性，且连盐高铁线路整体表现稳定。但在与BDS监测站验证过程中也发现，利用C波段SAR数据探测出的高相干性点在线路空间分布上具有不均匀性，在一些横向形变明显区域很难获得高相干性点。

1 研究区与数据集

1.1 SAR数据集

图1 研究区地理位置Fig.1 Geographic location of the studied HSR in Jiangsu province

连盐高铁位于我国苏北平原，沿线分布众多城镇，人口分布密集，平均每平方千米600～700人。在过去的几十年中，苏北平原地下水面临着过量开采的问题，地面沉降现象也因此变得严重起来。地下水开采和重大工程建设一直是引起该区域地表沉降的主要因素，特别是在过去15年，沿海城镇化的加速与扩张加剧着部分地区地表沉降，对人工线状构筑物带来了极大的影响。而高铁对轨道的稳定程度又要求极高，地表轻微的形变就会给高铁运营带来较高的潜在风险。为此，本文收集了高铁建成后检测阶段共47景TOPS成像模式的S1数据，具体时间跨度为2018年1月10日至2019年9月2日。S1数据可以免费从欧空局下载，入射角为33.6°，影像覆盖区域约为110×250 km2，像元分辨率在距离向是5 m，方位向是20 m。由于研究区在影像中相对较小，因而本文只对裁剪后的影像进行处理，裁剪区域长约30 km，宽约36 km，面积约为1080 km2。表1列出了影像的获取日期以及时空基线分布情况。选取2018年8月26日获取的影像作为主影像，其他从影像均配准并重采样到主影像上。

表1 Sentinel-1A数据获取时间以及干涉对参数

* 为参考主影像。

1.2 BDS实时形变监测系统

中国兵器工业集团下的北方信息控制研究院作为国内首家将北斗高精度卫星定位技术应用于监测铁路基准点坐标准确性及路基沉降的单位，负责连盐线北斗基础测量基准控制网的布设。布设的控制网分为基础控制网和加密控制网，其中基础控制网由15个基准点组成，间距2 km，全长29.60 km；加密基准点布设在灌河特大桥的CPⅡ控制点处，共布设12个基准点，间距500 m，全长10.396 km，形成北斗铁路加密测量基准控制网，为目标路段范围内的变形监测及轨道控制提供测量基准。此外，在路基段四标段一工区选择3个监测断面，共布设9个路基沉降监测点，各断面之间的距离约为80 m，进行基于北斗的铁路路基沉降监测应用研究，并基于此研究确定路基沉降监测的方式及预警机制。各个BDS基准点、加密点和断面位置如图2所示。由于各个监测站与基准站距离相对较近，故BDS监测数据处理主要是采用静态相对定位[16]。对于短基线的相对定位，大气延迟(对流层延迟和电离层延迟)能较好地被削弱或消除。其中采用LAMBDA方法固定整周模糊度[17]，采用高次差法和Turbo Edit组合法进行周跳探测与修复。

图2 连盐线BDS形变监测系统装置Fig.2 Photos of BDS deformation monitoring system alone Lian-Yan HSR

2 算法与结果分析

相比于城市地面沉降监测等应用领域，高铁等这种线状人工地物具有距离长、跨度大、分布不均等特点，给时序InSAR数据处理造成了一定困难，主要表现在：①高铁在影像上分布不均，如走向上占很多像元，而铁路的宽度可能只占几个像元；②大部分线路位于非城市区域，成像范围内后向散射强度低，经典PS选取方法得到的空间点分布较为稀疏，无法满足沉降监测的要求，而简单放宽选取条件可能引起大量错选；③PSI算法无法像SBAS算法一样可以使用相位环闭合差进行解缠相位误差检验，当解缠相位含有明显误差时，一般最小二乘估计的结果会发生偏差。针对上述现实问题，本文对常规时序InSAR数据处理流程进行改进，通过应用适当的策略，具体包括采用全分辨率干涉、相位稳定性分析探测PS点和基于L1范数的抗差M估计解算时序形变等。

2.1 MT-InSAR技术及数据预处理

考虑到本文监测对象是线状人工地物，影像时间间隔小且影像数量也足够，故采用最具有代表性的PS-InSAR处理流程，并采用全分辨率法提取小尺度形变。PS-InSAR处理流程具体包括主影像的选择、影像配准、差分干涉图的生成、PS点的探测、PS点网络构建、相位解缠、大气效应去除、形变速率的估算。本文数据预处理利用POD(precise orbit ephemerides)精密定轨星历数据去除平地相位，定位精度可达5 cm。利用30 m分辨率的SRTM DEM数据模拟地形相位。由于雷达是采用侧视成像方式，波束斜向照射地表时会导致雷达图像中的叠掩和阴影现象[18]，借助DEM(digital elevation model)的局部入射角影像和叠掩、阴影影像补偿由倾斜地形引起的SAR后向散射。在将DEM配准、重采样到SAR坐标系过程中，需要利用基线信息模拟干涉条纹。经验显示，即使是利用精密轨道数据计算得到的基线信息也不能保证完全正确，生成的差分干涉图中会存在以斜坡形式、与地形或系统有关的残余相位。因此，需要对生成的差分干涉图进行快速傅里叶变换(FFT)确定残余的基线分量，对基线模型进行改正，以达到最佳的干涉效果。

2.2 主图像的选取

PS-InSAR算法首先需要根据总体相干性最大的原则选择主影像，然后将其他从影像与之配准，生成差分干涉图集。公共主影像的选取需要顾及时空基线和多普勒质心频率差，其函数模型[19]为

（3）完善国家创新体制机制，提高政府管理效率。从美国创新战略可以看出，美国政府将民众、学生和企业加入创新体制中，共同关注世界热点问题，解决现实难题。在全球科技创新主体多元化、创新形式多样化的大背景下，为满足我国科技创新发展需求，提升科技成果转化效率，我国的科技创新体制机制也需要进一步的完善和更新。政府在科技创新中的地位和作用有待进一步明确，在坚持以市场为导向、企业为主体、政策为引导、产学研协同创新的体制机制下，确立企业在产业导向的科技计划中决策者、组织者、投资者的地位。与此同时，无论是从法律、法规层面的约束，还是普惠性政策的实施，都要确保实施过程得到了有效的监督、评估，提高政府的管理效率。

(1)

式中，γtotal为综合选取系数；T为时间基线；B为垂直基线；FDC为多普勒质心频率差；C表示各自的临界参数值。根据式(1)计算每一幅影像的选取系数，系数最大的作为主影像。最大综合选取系数的影像为2018年8月26日获取的影像，将其作为主影像，共生成46个干涉对，平均相干性为0.82。连盐高铁区域生成的时空基线图如图3所示。

图3 46对差分干涉对时空基线分布Fig.3 Baseline-time plots for the 46 differential interferograms

2.3 高相干点目标的选取

像素子视相干系数阈值法、相位离差法和振幅离差指数法是进行PS点探测最常用的几种方法[20]。这几种方法均是利用高相干点的稳定相位特征或者低相位离散特性，但是各有优缺点，实际数据处理中需要结合其中两种或多种方法进行高相干点的选取。由于连盐铁路沿线大部分位于非城市区域，而振幅离差方法适合城区高相干点的探测，故采用常用的振幅离差方法探测出的PS点密度可能不高，不利于分析线路的沉降细节特征。为了增加沿线区域PS点的密度，本文采用相位稳定性分析方法。该方法判定点在时空范围内的相位特征和基于点的空间距离越小相位相关性越高的假设。因此，首先进行空间低通滤波，然后计算相位离散度，进而选取具有稳定相位的点。这种选取方案，既保证了点位密度，也为三维相位解缠提供足够多可靠的点。首先利用振幅离差指数阈值(阈值设为0.4)方法，选择在长时间间隔内保持稳定的候选点集(PSCs)；然后运用相位稳定性分析估计每个候选点的噪声相位，将PS点概率函数值作为权重迭代进行精化。基于计算出每个候选点的时序相位相干因子γx值并结合振幅离差选取PS像元。γx计算公式为

(2)

图4 研究区PS点分布Fig.4 Distribution maps of PS in the study area

2.4 沉降速率和沉降量分析

对于监测高铁这种线状人工地物，C波段S1数据的分辨率略低于TerraSAR-X和COSMO等分辨率达1～3 m的卫星影像数据。为了尽可能地提高可监测形变梯度，文中采用全分辨率法提取小尺度形变。所谓的全分辨率干涉就是不对影像进行多视操作，直接使用SLC数据进行干涉处理。这样既可以保证对强相干目标如铁路、桥等构筑物准确识别，也提高像元中强散射体占主导地位的概率，获得更多的高相干点目标，不会造成孤立强散射体。

基于选取的高相干点目标，采用3D相位解缠方法[22]：①在时间域内，基于这些PS点建立Delaunay三角网，建立相邻两点之间的相位模型，经过去除大气效应和残余的轨道误差后，对时间序列相位差进行低通滤波和相位解缠；②运用解缠的差分相位时间序列，对每一个解缠相位值构建一个费用函数，最终得到解缠后的相位如下

Δφx,i=W{Δφx,i}+2kx,iπ=Δφh,x,i+φdef,x,i+φorb,x,i+φatm,x,i+φnoise,x,i

(3)

式中，Δφx,i为第i对干涉图像元x的解缠相位；kx,i为相位整周数，若解缠精确，则kx,i对干涉图上大部分像元是相同的；W{·}代表缠绕；Δφh,x,i为DEM不精确导致的地形残差相位；φdef,x,i为视线向上形变相位；φatm,x,i为传播过程中大气引起的路径延迟相位；φorb,x,i为残余轨道误差相位；φnoise,x,i为热噪声、配准误差等引起的噪声相位。

为了获取最终的形变相位，还需要扣除各噪声分量。对于轨道误差相位，文中采用二阶多项式对方位向和距离向进行拟合去除。对于残余高程和形变相位的估计，一般采用最小二乘(LS)估计。然而LS估计没有考虑到残余轨道误差和大气相位的影响，它只是一个理想的全局最优解。只有当观测相位值是相互独立并且服从正态分布，解算的结果才是无偏和有效的，当解缠相位观测值含有明显残余误差时，LS估计的结果会发生明显歪曲。因此，本文在解算参数时，采用抗差M估计[23]为

(4)

由于试验区地形平坦，大气中垂直分层效应不明显，湍流效应占据主导地位。而湍流大气一般被认为在时间和空间上是随机分布的，为了避免引入大气误差，数据处理中并没有直接利用外部数据进行大气相位的去除，而是采用时空滤波进行大气相位估计，最终得到平均沉降速率图和监测周期内时序位移序列。图5为分别利用抗差M估计和最小二乘方法解算的形变速率标准差分布情况，可以看出利用抗差M估计解算的形变速率标准差比LS估计的结果普遍要小，形变速率标准差降低的PS点约占89.2%。图6为利用抗差M估计得到的线性沉降速率结果。

图5 形变速率标准差分布图Fig.5 Distribution map of the standard deviation of deformation rate

图6显示了利用时序InSAR获取的整个区域在2018年1月至2019年9月期间的平均线性形变速率，整个区域最大形变速率达到了15 mm/a，而高铁沿线区域形变更是微小，平均形变速率约为5 mm/a。在将近两年的监测时间内，人工构筑物保持了很好的相干性。就高铁沿线区域而言(图6(b)所示)，尽管相较X波段，C波段的S1数据地面分辨率中等，但本文采用的是全分辨率干涉，保证了最高分辨率，共获取10 199个PS点，密度远远超过BDS监测点。结果表明利用C波段S1数据监测像高铁这类人造线状地物目标，具有很大的应用潜力。

图6 线性形变速率Fig.6 Linear deformation velocity map

3 沉降监测结果对比验证及分析

MT-InSAR技术监测的是视线向形变，而BDS监测的是垂直方向和水平方向形变，为了验证MT-InSAR的结果，一般假设没有水平变形，将InSAR视线向形变转换到垂直方向上。根据雷达成像几何，将InSAR监测结果除以cosθ(θ为入射角，每个PS点入射角不同)即可完成转换。BDS监测数据后处理采用短基线相对定位，精度在亚太地区优于GPS，在垂直向优于2 cm，水平向上优于1 cm[24-25]。关于InSAR结果的验证，本文将从两个方面展开，即依据BDS监测结果验证比较InSAR的平均形变速率和时序位移量结果。

3.1 平均形变速率结果验证

根据BDS测站的位置选取相应PS点，将计算得到的BDS形变速率与PS结果进行比较。值得注意的是，BDS监测时间段与InSAR数据时间范围并不完全一致，仅是其一部分，因此本文仅比较两者在重叠时间内的形变速率，并引入均方根误差(RMSE)评价精度。RMSE计算公式为

(5)

式中，Δvk为第k个BDS监测站形变速率与InSAR结果形变速率之差；n为测站点的个数。

其比较结果见表2和图7所示。

表2 BDS监测点与PS点形变速率结果比较

由表2可以看出，BDS监测点与PS点形变速率均方根误差为6.2 mm/a，最大偏差为13.9 mm/a，且偏差落在5 mm/a内占70.6%。偏差大于1 cm/a仅有2个测站(图7(b))，对应图7(a)中的序号2和17，回归分析结果负相关系数R为0.79，表现出很好的一致性。关于图7(a)中序号2和17对应的测站出现较大偏差的原因，4.2节将详细分析。图8显示了MT-InSAR在时间跨度为21个月内，高铁沿线各个BDS监测点位置的平均形变速率变化情况，从图8中可以发现整条线路平均形变速率在1 cm/a内，线路比较稳定。为了对比验证MT-InSAR监测的精度，又分别获取了二者在重叠时间内(2018年11月～2019年9月)的沿线平均形变速率变化情况，可以发现InSAR结果与BDS监测的结果较为吻合，趋势相近。

图7 BDS监测与MT-InSAR监测结果回归分析与误差直方图Fig.7 The regression analysis and error histogram of velocity and displacement comparison between PS results and BDS measurements along a railway

图8 高铁沿线各个BDS监测站与PS点平均形变速率结果Fig.8 Average subsidence velocity comparison of BDS monitoring stations along the railway

3.2 时序位移量结果对比验证

以BDS监测站天线位置为中心，100 m作为半径搜选所有PS点，并求取该范围内所有PS点的平均时序位移序列。为了匹配两种技术获取的位移序列，需要将位移时间序列转化为重叠时间内共同的时间参考。由于篇幅所限，文中只显示了其中10个监测站与MT-InSAR获取的时序位移量的比较结果。

为了比较MT-InSAR监测结果，本文将2018年11月6日设为时间参考基准点，求取BDS与MT-InSAR相对基准点在垂直方向上的位移变化量。从图9中可以看出，除了JM0006Q和JC0006Q测站，其他测站的结果十分吻合。JM0006Q和JC0006Q测站偏差相对较大是因其周围的PS点很少，且距离周围最近的PS点相对较远。选取的PS点并不能完全反映这两个测站的真实形变情况，而地表形变具有空间相关性，故虽然整体略有偏差，但整体形变趋势较为相近。从图9中还可以发现，2018年11月之前，形变位移序列起伏较大，例如JC0007Q、JC0015Q站，除了有一部分原因是残余大气的影响，还可能是因为前期施工过程中地基需要一个稳定的过程，即下沉-抬升过程。

为了揭示前期沉降过程，利用同样的算法，但只处理了前期即2018年间获取的影像，结果如图10所示。对比图6形变速率，可以发现前期地面沉降更为明显。尤其是第2段线路(图10(c))，形变速率最大达到2 cm/a。从早期整个线路的形变来看，可以看出前期施工后一段时间内线路出现地基回升与缓慢下沉，这在铁路建造过程中是常出现的现象。再结合图9，可以发现线路在进入2019年后，地基逐渐变得稳定，无论抬升还是下沉均在数毫米范围内波动。

图9 InSAR与BDS在垂直方向上时序位移对比结果Fig.9 The displacement comparison along the vertical direction between InSAR and BDS results along Lian-Yan HSR

表3为各个BDS监测站与其相应的PS时序位移统计结果，其中最大均方根误差为6.1 mm，最小为1.6 mm，平均均方根误差为3.8 mm。选取的所有样本点中最大偏差为14.1 mm，偏差位于2 mm内的样本占样本总数的47%，最小偏差为0，二者偏差大于5 mm仅占19%，以上可以看出二者结果是十分吻合的。

表3 BDS监测点与PS点均方根误差

3.3 连盐高铁沿线地面沉降分析

为了更加全面的分析连盐线的形变特征，图11可视化了3个区域(A、B和C区域)中的一些探测点的形变序列。在A和C区域中，探测到的PS点主要分布在线路上以及相邻的人工建筑物。其中A区域位于第1段线路(图10)，选取的两个PS点，形变速率分别为10.7 mm/a和12.3 mm/a。从其形变序列上可以看出，二者有着相近的形变趋势，形变主要发生在2018年5月至8月，具体表现为抬升，幅度大约在1 cm左右，8月之后，地基逐渐变得稳定。C区域主要位于第2段线路灌河特大桥北侧，可以发现区域内的形变特征虽然和A区域相近，但是在2018年间形变更为剧烈，这也验证了图10得到的第2段线路较第1段线路形变速率更为明显的结论。推测这可能是由于该区间建成时间比第1段线路晚，且位于软土分布区，主要表现为高含水量、高压缩性、低承载力，因此地表固结稳定的时间相对延后。相比A区域，C区域在2018年5月至2018年8月地基回弹更为剧烈，时间更长，直到12月中旬线路开通后才逐渐稳定。无论A区域还是C区域在线路开通运营后，线路均变得较为稳定。而B区域是连盐线路基段四标段一工区，区域内有3个监测断面，共布设9个路基沉降监测点，各断面之间的距离约为80 m。可以发现在该区域并没有探测到相干点，显现了高相干性点在线路空间分布上具有不均匀性的特征。但是从9个监测站的监测结果来看，3个断面在2018年11月至次年4月的监测周期内，没有明显的形变发生，只在1～2 mm之间变化。以上可以得出结论，连盐线总体表现稳定，尤其是12月中旬线路开通运营后。

图11 连盐线重点监测区域时序形变Fig.11 Time series deformation of the points in the areas A, B, C along the Lian-Yan Line

4 结 语

本文将C波段Sentinel-1A数据应用到连盐高铁沉降监测，采用多时相InSAR技术进行处理并与同时期获取的BDS监测结果进行对比验证分析。结果发现，在连盐高铁建造完成至试运营阶段，平均形变速率在1 cm/a之内，整体线路路基在2018年11月之前经历一段时间的沉降与回弹后，11月以后，尤其在进入2019年后，沿线逐渐平稳。与BDS监测数据验证结果来看，二者线性平均速率均方根误差为6.2 mm/a，最大偏差为13.9 mm/a，最小偏差为1.1 mm/a，相关系数达到了0.79；10个BDS监测站平均时序位移均方根误差为3.8 mm，且与BDS监测结果呈现微小的系统误差影响，偏差在2 mm内约占47%。无论是获取的平均形变速率还是时序形变序列，均取得很好的一致性，显现出利用C波段S1数据对人工线状地物进行时序InSAR监测的巨大潜力。同时在监测过程中，虽然利用多时相InSAR技术获得了相当高密度的PS点，但是试验发现在3个重点监测的断面处，并没有获取足够的PS点，导致无法与布设的BDS监测点进行比较，这也暴露了在利用C波段SAR数据进行时序InSAR监测过程中，虽可以获取可观密度的PS点，但是在空间分布上具有不对称、不均匀的特点。因此，综合考虑高铁线路的走向、坡度以及雷达成像视角，探测沿线区域分布更加均匀、稳定可靠的相干性点，并获取高铁线路的三维形变场，有待进一步深入研究。

致谢：感谢欧空局提供的2018年1月至2019年9月的Sentinel-1A数据。