王莉，王志明，杨莹

1.西藏自治区重点公路建设项目管理中心, 拉萨 850001；2.交通运输部天津水运工程科学研究所, 天津 300456

0 引言

随着中国国民经济的快速发展，道路交通网络迅猛扩张，但是区域性的地表形变给道路交通的安全运行带来了许多潜在危害。因此，对道路交通区域地表形变的监测日益受到关注[1--2]。

传统地表形变监测主要通过精密水准测量、GPS测量等，采用人工方式展开，精度很高，但存在效率低、时效性差和覆盖范围有限等问题，很难实现大范围的动态连续监测。而合成孔径雷达干涉测量(InSAR)技术，因具有全天时、范围广、受天气影响小等优点，可以获取大范围、高精度(mm级)的地表形变信息，在地表形变监测方面得到了越来越广泛的应用。Heleno et al.利用永久散射体干涉测量(PS--InSAR)技术对Lisbon地区的主要公路和铁路线路进行沉降监测，在道路沿线发现了两个新形成的沉降区[3]。李英会基于稳定点目标分析方法，利用PS--InSAR和小基线集干涉测量(SBAS--InSAR)方法开发了一种新的InSAR技术，对京津高铁天津段2009—2010年的地面沉降进行监测，发现区域最大沉降量达400 mm[4]。Perissin et al.使用InSAR技术监测上海的地表形变情况，发现地铁隧道和公路是地形沉降的主要发生区域[5]。杨艳等分析了京津城际高铁线路运营五年多来的差异沉降与坡度变化特征，认为区域沉降导致的差异沉降量较大，但沉降变化仍符合高铁轨道的设计要求[6]。张学东等融合PS--InSAR技术和SABS--InSAR技术，分析了京沪高速公路和上海市地铁沿线的地表形变，得到了较好的结果[7]。刘欢欢等监测了京津高铁北京段地面沉降情况，并对其成因及影响因素进行了分析[1]。与传统方法相比，PS--InSAR技术在道路区地表形变监测方面，表现出精度较高、连续性好的优势。

但是，当前尚缺乏对中国东北高铁运行区地表形变的监测。哈尔滨到大连高速铁路(哈大高铁)是中国“四纵四横”快速铁路网京哈高铁的重要组成部分，对东北地区的经济社会发展具有重要战略价值；也是国内开通的首条高寒地带客运专线，运行路线经过松花江、拉林河、饮马河及辽河等多条大河，恶劣的气候、不稳定的地下特征和复杂的地质结构都会带来不可避免的地表形变，从而给高铁的安全运行带来潜在风险[8--9]。因此，笔者以跨越松花江、饮马河的哈大高铁扶余至德惠段为例，拟采用Sentinel--1B卫星获取的合成孔径雷达(SAR)数据，尝试基于PS--InSAR方法，进行哈大高铁运行以来的地表形变分析。

1 研究区域和数据

1.1 研究区概况

哈大高铁北起哈尔滨西站、南至大连北站，线路全长921 km。高铁的运行极大缩短了中国东三省主要城市间的时空距离，为东北区域的经济一体化创造条件，是中国东北地区的干线铁路之一。哈大高铁也是世界上首条在高寒地区运行的高铁，路基冻胀等地表形变也是高铁安全运行面临的主要问题之一。

为了验证基于PS--InSAR技术监测高铁行进区间地表形变的可行性，选择哈大高铁扶余至德惠段作为研究区(图1)。研究区内高铁全长约60 km，大部分位于中国东北季节性冻土区；高铁运行区间经过了松花江及其重要支流饮马河，高铁沿线的土地利用类型以旱田、水田为主。耕地的季节性特征、城市的扩张，都增加了铁路沿线地表形变监测的难度。

图1 研究区域示意图(底图为Sentinel--2B影像)Fig.1 Sketch map of study area

1.2 Sentinel--1B SAR数据

Sentinel--1B卫星是欧空局(ESA)实施的涉及环境和安全信息服务的哥白尼计划的一部分，搭载有C波段(3.8～7.5 cm)合成孔径雷达(SAR)，主要用于陆地、海洋和海冰监测服务。卫星于2016年发射，轨道高度为693 km，重访周期为12 d；中心入射角为38.89°，距离向分辨率为 5 m，方位向分辨率为20 m。合成孔径雷达采用了4种工作模型，本研究选择了宽幅干涉模式成像，VV极化方式。所有影像均为单视复数类型的数据产品。

研究中，选择了2016年12月至2019年11月的32景SAR数据，覆盖研究区近3年时间，采样的平均时间间隔约为34 d。

另外，在PS--InSAR数据处理方面，卫星星历轨道信息在不同时相SAR的影像配准和形变信息生成等环节都起着重要作用。原始的卫星轨道信息存在诸多误差，会造成干涉图中存在因基线误差而出现的残余条纹误差。因此，基于精密的卫星星历轨道数据对卫星原始轨道信息进行校正，可有效去除因轨道误差引起的残余条纹误差[1,5]。为了确保Sentinel--1B数据的轨道精度，本实验采用的是欧空局公布的POD精密定轨星历数据，定位精度优于5 cm。

从美国太空总署官方网站下载了由美国太空总署和国防部国家测绘局联合测量的航天飞机雷达地形测绘使命(SRTM--DEM)数据，空间分辨率为30 m，为Sentinel--1B SAR数据的处理提供参考DEM，用于SAR影像集差分干涉处理时的地形相位去除[1,5]。

实地调查发现，铁路附近90%以上的区域为耕地，平均气温、地表特征随季节的变化很大；而PS--InSAR技术可以较好地抑制由于季节性温度变化和植被覆盖等因素引起的失相干现象。因此，本研究采用PS--InSAR技术进行研究区的地表形变分析。

2 数据处理

2.1 PS--InSAR方法

PS--InSAR技术是Ferretti et al.提出的，利用同一区域的多时相SAR影像，通过适当的数据处理和分析程序，将形变相位分量与其他相位分量分离[10]。PS--InSAR技术的关键是确定PS点，用以估计形变时间序列和形变速率。PS--InSAR方法的另一个结果是地形残差，即给定PS点的散射相位中心的真实高度和该点的DEM高度之差，这也是获得精确PS地理编码的关键参数之一[10]。PS--InSAR技术处理流程包括以下几个步骤：

SAR影像定标和配准进行SAR影像配准和辐射定标，将时间序列数据的几何位置和辐射强度统一起来，用以消除成像过程中受卫星视角、轨道高度、入射角、大气效应以及SAR系统热噪声等造成的振幅信息变化。

生成差分干涉数据集选取其中一个时相的SAR数据作为主影像，其他时相数据作为副影像，分别进行干涉处理，生成干涉图集；利用DEM数据去除地形相位，生成时间序列差分干涉图集。

PS候选点的识别利用振幅离差指数法与相关系数法双阈值法进行PS候选点的识别，其中，相关系数双阈值法能够快速缩小PS点探测范围，振幅离差指数法能够逐一进行PS点的精确筛选。

基于多图像稀疏网格的相位解缠利用Delaunay不规则三角网建立PS候选点之间的连接关系，采用加权最小二乘法进行相位解缠。

大气相位的估计和去除在估计每个PS点的线性畸变和DEM误差后，从初始差分干涉相位减去残余相位，包括非线性畸变相位和大气相位，采用时域滤波和空间滤波进行分离。在估计时间序列SAR影像上PS候选点的大气相位后，通过Kriging插值方法，得到SAR影像任意像素的大气相位。

PS点的时序分析与形变量估计利用相干系数法，确定最终的PS点，重复第(4)步，重新构建Delaunay不规则三角网，进行多影像稀疏格网相位解缠，获得所有PS点的时序形变值与形变速率。

2.2 Sentinel--1B SAR数据处理

本实验使用Esri公司SARscape软件进行PS--InSAR处理，选取高铁沿线左右各延伸12 km作为研究区域。根据所选择的32景Sentinel--1B SAR数据，时空基线的最优值(图2)，选择2018年8月10日的SAR数据作为主影像，其余31景数据为从影像，分别与主影像进行配准和干涉处理。

经过对比分析，选取时成像时间为2018年8月10日的SAR数据作为主影像，共生成31个干涉对，其中最大时间基线600 d，最小时间基线为24 d(图2)。

图2 Sentinel--1B SAR数据的时间位置(a)和Sentinel--1B SAR干涉对的时间基线(b)Fig.2 Time position of Sentinel--1B SAR data (a) and time baseline of Sentinel--1B SAR interference pairs (b)

首先进行Sentinel--1B SAR影像集预处理，包括SAR数据导入、研究区裁剪和SAR影像配准、辐射定标处理。然后，根据所选择的32景SAR数据，制作了空间基线和时间基线分布图(图2)；根据空间基线和时间基线的最优值和多普勒质心差阈值，选择2018年8月10日影像为主影像，由此将数据集划分为不同的子集，生成31个干涉测量对；这样，对研究区地表形变的平均采样间隔约为34 d，最大时间基线600 d，最小时间基线为24 d，平均绝对基线和最大绝对基线分别为38.01 m和106.2 m。

利用SRTM DEM数据去除地形相位，利用POD精密轨道数据消除因轨道误差导致的基线误差和卫星系统误差导致的平地效应相位，生成31幅时间序列差分干涉图集。将振幅离差指数法与相干系数双阈值法相结合用以识别和确定候选PS点；基于Delaunay三角网法，确定相邻PS点间最优解关系，并进行空间上的相位解缠处理，利用已知的地面控制点解算任一点的线性形变速率与高程残差；基于时空滤波处理，去除大气相位误差，进而计算各PS点整体相干性；进行二次解缠处理，解算最终PS点的稀疏格网，得到雷达视线(line of sight, LOS)方向年平均形变速率分布图(图3)。

3 结果与讨论

3.1 研究区年平均地表形变特征

图3是高铁扶余至德惠段LOS向的年平均形变速率图。图3表明，从2016年12月到2019年11月，高铁沿线的年平均形变速率为-16.8～9.1 mm/a，形变速率集中在-2.4～2.0 mm/a的变化区间，表明研究区总体的地表形变特征是比较稳定的。

图3 雷达LOS向的年平均地表形变速率分布图(底图为Sentinel--2B影像)Fig.3 Annual average surface deformation rate map in radar LOS direction

研究中共得到51 169个PS点，其中发生上升形变的PS点有19 626个，占总数的38.4%；发生沉降形变的PS点31 543个，占总数的61.6%；表明研究区的地表形变还是以沉降为主。发生上升形变的PS点集中在德惠市区的南部，结合实地调查和Google Earth影像，表明上升形变与城市区的扩张有关。沉降形变集中在扶余市区北部和南部、德惠市区西北部和松花江两岸。其中，松花江附近的地表形变均为负值，最大的年平均沉降速率可达到-12.2 mm/a。而在扶余市区，地表形变特征表现复杂，西部和南部以沉降为主，最大沉降量达16 mm/a；而中部和东部则是上升形变和沉降形变同时出现。德惠市区以上升形变为主，仅在中部沉降形变以线性形式呈北西向展布。

在扶余至德惠段的高铁沿线，地表形变特征表现特殊，以上升形变为主。统计直方图(图4)表明，高铁沿线的年平均形变速率集中在-1.8～2.0 mm/a，整体的年平均形变速率均值为0.44 mm/a。共采集PS点1 295个，发生上升形变的点883个，占总数的68.2%。发生沉降的PS点集中在松花江附近和扶余市区附近。

图4 高铁沿线PS点的年平均地表形变速率直方图Fig.4 Annual average surface deformation histogram of PS points along high-speed railway

3.2 典型位置累积地表形变特征

为了更清楚地掌握高铁沿线的地表形变特征，选择4个典型PS点位置(图3)，分别是扶余北站(A点)、松花江大桥(B点)、耕地(C点)和德惠西站(D点)，根据所生成的31个干涉影像对，制作了典型位置累积地表形变量的时序变化曲线(图5)。图5表明，A点和D点表现出了总体的沉降特征，B点的沉降量保持不变，而C点则出现了缓慢上升的地表形变特征。

图5 典型PS点(图3)地表累积形变量的时序变化曲线Fig.5 Variation of accumulated surface deformation with time in typical PS points as shown in Fig.3

但是，累积地表形变量随时间的变化幅度非常大。A点在扶余北站附近，从2016年12月至2017年3月，累积地表形变呈明显增加的趋势；而在2017年3至4月，出现大幅度下降；从2017年4月至8月，累积地表沉降呈上升趋势；之后，至2017年12月，累积地表形变呈下降趋势，仅在11月份有小幅波动；从2017年12月至2018年3月，累积地表形变呈快速增大趋势；至2018年5月，再次呈直线下降趋势；随后，至2018年9月，累积地表形变表现为缓慢上升趋势，之后，至2018年12月，再次呈现快速沉降特征；2018年12月至2019年4月，总体表现为上升趋势；随后至2019年7月，再次表现为直线下降趋势；2019年7月至11月，总体上表现为上升趋势。

B点在松花江大桥附近，累积地表形变随时间的变化表现复杂。从2016年12月至2017年4月，表现为增大的沉降形变；2017年4月至10月，累积地表形变波动很大，而2017年11月至2018年2月，表现为上升形变；2018年2月至11月，累积地表形变波动再次增大；2018年11月至2019年3月，再次表现为上升形变；之后，累积地表形变再次呈大幅波动形式。

C点在接近德惠的典型耕地区，2016年12月至2017年4月，累积地表形变量呈明显增大趋势；随后，至2017年8月，呈沉降趋势；2017年8月至2018年3月，累积地表形变一直呈缓慢增大趋势，仅在2017年11月份有小幅波动；2018年3月至10月，地表形变波动很大；2018年10月至2018年12月、2019年2月至3月，累积地表形变均表现为增大特征，中间出现小幅减小特征；随后，至2019年7月，呈现快速下降趋势；7月至11月，累积地表形变特征表现复杂。

D点在德惠西站附近，2016年12月至2017年1月，地表形变有小幅增大；随后，至2017年6月，呈缓慢下降趋势；8月有小幅上升，之后再次沉降；至9月到达最大累积沉降量；2017年9月至2018年2月，累积地表形变再次表现为上升趋势；之后，至2018年4月，再次呈现沉降特征；2018年4月至10月，累积地表形变量随时间的变化波动很大；2018年10月至2019年2月，累积地表形变表现为快速上升趋势；2019年2月至6月，表现为沉降特征；随后至2019年11月，再次出现较大波动的累积地表形变量。

总体上，4个典型位置的累积地表形变都表现出了随时间的明显变化特征。

3.3 高铁安全运行的潜在风险分析

综合研究区总体的地表形变特征，认为高铁扶余至德惠段安全运行的潜在风险有以下两个：

(1)地表的季节性形变特征

对4个典型区的3年累积地表形变分析结果表明，无论是在高铁站点(图5，A点、D点)，还是河流(图5，B点)、耕地区域(图5，C点)，地表在每年12月至第二年2月(或3月)表现为上升形变，3月至6月表现为沉降形变，6月至10月的地表形变量变化很大。结合东北区的气候特征，每年12月至第二年2月(或3月)的上升地表形变，明显是寒冷气候条件下由于土壤的冻胀生成的；而3月至6月的沉降形变应该与冻土层的融化有关；出现较大波动地表形变的6月到10月，正好处于东北地区的雨季。

因此，东北地区高铁行进路线的地表形变，主要是季节性的路基冻胀和降雨造成的。

(2)高铁线路附近的潜在风险

从年平均沉降速率图(图3)可以看出，在德惠西站和扶余北站附近都出现了大面积的地表沉降；从累积地表形变曲线(图5)来看，这两个点的时序变化曲线与其他两个点的区别，在于每年8月或9月有一个明显的上升形变。

结合东北区的气候特征，初步分析很可能是地下水开采造成的地下漏斗，总体的沉降与沉降漏斗的加剧有关；而在8月或9月的上升形变，则是与该时期较丰富的地表降水有关，减少了地下水的开采量，并对地下水产生一定的补充。目前，尚缺乏扶余地区的地下水勘探资料。徐世明和李虎等分别在2005和2002年的地下水资源调查中，已明确指出，扶余地区地下水位呈连年下降趋势[11--12]。另外，图3表明，扶余市区的地表形变特征表现复杂，市区中部和东部的上升形变和沉降形变同时出现，进一步证实了该结论的可靠性；这里，上升形变很可能与城市东扩有关，而沉降形变则很可能来自新居住区的地下水开采。

总体地形沉降量在德惠市<4 mm，对高铁的安全运行影响不大；扶余地区的累积地表形变量在2019年7月之前较为稳定，之后变化幅度加剧，变化量接近或超过5 mm，这已经对高铁的安全运行具有潜在的风险。

4 结论

(1)总体上，高铁行进区间的地表形变表现稳定，LOS向年平均形变速率在-2.4 ～2.0 mm/a之间；尤其是在松花江大桥附近，累积地表形变量随时间的变化表现稳定,表明研究区的地表形变特征符合高铁安全运行的要求。

(2)季节性的路基冻胀和降雨是东北高铁区的地表形变的主要影响因素。

(3)在德惠西站和扶余北站附近，都发现了可能与地下水开采有关的明显地表沉降；其中扶余北站的地表沉降有加剧趋势，已对高铁的安全运行构成潜在风险。

由于缺少高铁线路和地下水位监测数据，该研究未进行地表形变监测精度的评价。但从雷达LOS向的年平均地表形变速率分布图(图3)和高铁沿线PS点的年平均地表形变速率直方图(图4)来看，研究区的地表形变速率在合理的变化范围内，表明了该研究成果的可靠性。后续将结合精确的地表形变测量结果，进行哈大高铁沿线的长时间序列地表形变监测，为高铁的安全运行提供科学的数据支撑。

致谢感谢美国地质调查局STRM DEM数据和欧空局Sentinel--1 B SAR数据的免费使用。