京津冀地区高速铁路沿线区域地表沉降监测及时序演化态势分析

杨 斌，刘国祥，姚京川，李 政，冯海龙，毛文飞，张 瑞，吴明杰

(1.中国国家铁路集团有限公司 工程管理中心，北京 100844； 2.西南交通大学 地球科学与环境工程学院，四川 成都 610031； 3.中国铁道科学研究院集团有限公司 铁道建筑研究所，北京 100081； 4.雄安高速铁路有限公司，河北 雄安 071800)

地表沉降是自然或人类活动引起的一种具有缓变性、不均匀及不可逆性的地质灾害[1-2]。长期以来，因我国北方干旱少雨持续超采地下水，致使多地存在显著的不均匀地表沉降。随着时间的累积，已对铁路基础设施结构稳定性及轨道平整性等造成严重影响[1-3]。据统计，截至2009年，我国有8万km2的地区累计沉降量超过200 mm，包括多个知名城市及周边区域。因此，及时把握区域性沉降的空间分布、及时探明其发展趋势与演化态势，对于保障重大交通基础设施安全运营和制定相关管控政策，均具有重要的参考价值和意义。

在地表沉降的监测领域，基于时域建模分析的合成孔径雷达时序干涉(TS-InSAR)技术克服了时空失相干、轨道和地形累积误差、相位解缠误差以及大气延迟误差等负面因素[4-11]，能够准确地提取长时间尺度上的地表形变信息，现已广泛运用于各种地灾调查、冻土监测[12]及地表沉降监测等。其中，在算法上最有代表性的技术方法包括Ferretti等于2000年前后提出的永久散射体雷达干涉(Permanent Scatterers InSAR, 即PSInSAR)[13-16]，Berardino等[16-17]提出的小基线集(Small Baseline Subsets, SBAS)差分雷达干涉方法，以及Hooper等[18]充分考虑地物的稳定性并利用三维解缠法有效消除干涉相位模糊度提出的基于后验估计的PSI方法[19-21]。而侧重于干涉点目标分析的技术方法，则主要包括Werner等[22-23]提出的IPTA算法、Zhang等[24]提出的时域相干点目标分析(TCP-InSAR)等，这些方法均可统称为PSI方法[25]。需要注意的是，PSI方法在进行PS点探测时，振幅离差指数(Amplitude Dispersion Index, ADI)小于0.25的像元方可认为对应于PS点目标，且该方法至少需要30景以上 SAR 影像进行PS探测以及形变建模和提取。

为了探究京津冀地区的不均匀沉降对高速铁路的影响，本文选取2015年11月至2018年3月间覆盖研究区域的51景C波段SAR影像，对影像进行精密配准后[1]，利用干涉点目标时序分析方法提取该区域的沉降信息，并结合该区域人类活动信息和水文地质资料等对高铁沿线的沉降漏斗进行演化态势分析，为该区域铁路安全运营及维护等提供数据支撑。

1 地表形变信息提取方法

1.1 干涉点目标时序建模原理

PSI技术以某一影像为主影像，将其他多幅影像配准到该主影像上，并经重采样处理，运用统计分析的方法获取时序的影像强度和相位信息，并识别出不受时空基线和大气延迟影响的点目标，其主要提取后向散射特性稳定的相干目标，对这些相干目标进行干涉处理，在处理过程中将差分干涉相位φdiff分解为[1]

φdiff=φl+φtopo+φres

(1)

其中，φres=φbaseline+φatm+φnoise+φnl

式中：φl为线性形变差分干涉相位；φtopo为外部DEM引入的地形误差相位；φres为残余相位；φbaseline为基线误差造成的残余相位；φatm为大气引起的相位；φnoise为噪声引起的相位；φnl为非线性引起的相位。

干涉点目标时序分析方法以差分干涉相位φdiff为因变量，高程改正Δhe和线性形变速率Δv为自变量建立相位函数模型，残余相位为常量建立的相位回归分析模型[1]为

φdiff=aΔvT+bΔhe+φres

(2)

式中：T为时间基线；B为空间垂直基线；λ，R和θ分别为雷达波长、雷达至地面的斜距和入射角。

利用式(2)在时空域进行回归分析，首先解算出Δhe和Δv，然后更新差分干涉相位，去除残余相位，接着进行迭代回归分析，最终得到与点目标观测值最佳拟合的形变速率和沉降值。

1.2 干涉点目标时序处理方法

干涉点目标时序分析方法主要步骤包括：

(1) 获取N(N≥30)幅SAR影像数据[9]，裁剪出研究区域，选取其中1幅为主影像，将其他影像与主影像进行精密配准；

(2) 生成N-1幅干涉图，经过去除地形相位和参考椭球相位等处理，生成差分干涉图；

(3) 选取振幅离差指数(ADI)阈值，选出高振幅的像素作为PS候选点，但考虑到大气相位(Atmospheric Phase Screen, APS)对形变速率和高程误差估算模型的影响，选点后，需要估算形变速率、高程误差以及APS[10]；

(4) 根据Ferretti等提出的二维周期图(Two Dimensions Periodogram, TDP)[11]解算方法，获取每个点的线性形变速率和高程误差；

(5) 对残余相位时空滤波，分离出大气相位、非线性形变和噪声相位(忽略不计)，最终得到研究区域内PS点的形变速率和高程信息。

1.2.1 PS点探测

振幅离差指数PS点选取是Ferretti等[10]提出，其主要是利用PS点的强反射特性，挑选出高振幅值的像素作为PS的候选点，而刘国祥等[26]通过对PS散射特性的稳定性进行分析，提出利用振幅离差指数阈值从PS候选点中进一步精选出PS。本文利用后者提出的方法进行PS点选取，其选取条件为

(3)

式(3)中第1个公式确定相干性，第2个公式确定的是稳定性。

1.2.2 时序建模及线性参数求解

永久散射体时序差分雷达干涉(PSI)是一种单一主影像的时序InSAR技术，利用1.2.1节中的方法进行PS点选取，借助外部DEM生成PS点在时间维度上的差分干涉相位。差分相位包括形变相位、大气延迟相位、高程误差及随机噪声，具体可表达为

φdiff=φref+φtopo+φdef+φatm+φnoise

(4)

式中：φref为参考椭球相位；φtopo为地形相位；φdef为真实形变相位。

为了克服大气延迟等影响，PSI技术采用邻域差分建模方法。相邻2个PS点在第i幅干涉图中的相位差Δφdiff,i可表示为

(5)

式中：Bi和Ti分别为第i幅干涉对的空间垂直基线和时间基线；(xl,yl)和(xp,yp)分别为第i幅干涉图中点l和点p的位置坐标； Δε(xl,yl;xp,yp)， Δφres，i(xl,yl;xp,yp;Ti)和Δv(xl,yl;xp,yp)分别为第i幅干涉图中点l和点p的高程误差、残留相位增量和LOS向的形变速率增量。

利用Ferretti等[11]提出的二位周期图解算方法，在|Δφres，i|<π的条件下，可基于多个缠绕的基线增量，通过使目标函数最大化，即式(6)，并搜索估计出二维解空间最优的Δε和Δv。

(6)

式中：γ为基线模型相干系数；Δωi为观测值与拟合值之差。

1.2.3 非线性形变及大气延迟相位

在差分干涉相位中，除了形变速率和高程误差外，还有非线性形变相位，大气延迟相位和噪声相位，即残留相位φres，i为

φres，i=φnl，i+φatm，i+φnoise，i

(7)

式中：φnl，i为第i幅干涉图的非线性形变相位；φatm，i为第i幅干涉图的大气延迟相位；φnoise，i为第i幅干涉图的噪声。

由于各幅干涉图相位的残留相位在时空域的特性不同，因此可对以上几个相位进行分离。利用空间低通滤波消除干涉失相关和其他随机相位分量，然后在时间域滤波，便可以通过时序变化趋势识别对应的非线性形变相位。得到的真实形变相位φdef为

φdef=φl+φnl

(8)

2 研究区域及试验数据

2.1 研究区域概况

研究区域为京津冀地区，主要包含雄安县、霸州市、廊坊市、北京以及部分县镇等，面积约为1.5万(km)2。该研究区域是京津冀一体化的核心区域，人口密度较大，经济开发区众多，交通基础设施较多且路网复杂。区域内大部分地区地下水开采处于超采状态；地质条件复杂，存在地裂缝、崩滑流和地面塌陷、地表沉降以及地下水污染等环境地质问题[27]；地基主要为软土、松软土层，压缩性较高[1]。研究区域如图1所示。

图1 研究区域

2.2 试验数据选取

选取研究区域2015年11月至2018年3月51景C波段的SAR影像，具体干涉对及其时空基线见表1，分辨率(距离向×方位向)约为2.3 m×13.9 m，宽幅约为250 km，2016年12月14日的影像被选为主影像，从影像与主影像的空间基线最长为257.104 6 m，最短为3.291 6 m，而时间基线最长为468 d，最短为12 d。此外，为了去除参考椭球和地形相位的影响，选取30 m空间分辨率的AW3D30数字地表模型(DSM)数据作为地形数据，该数据官方公布的高程精度为5 m标准差。

3 结果及分析

3.1 沉降

根据上面提出的技术方法，得到研究区域的沉降速率场，如图2所示。不难发现，研究区域局部沉降较为明显，沉降漏斗主要分布在雄县、霸州市以东的胜芳镇和左各庄镇、固安县、廊坊市以及北京通州区，全局沉降速率分布在20～206 mm·a-1之间。其中沉降最为严重的是胜芳镇，最高沉降速率达206 mm·a-1，漏斗中心累积沉降量为248 mm；左各庄镇、通州地区、雄县和廊坊市的沉降漏斗中心沉降速率分别为159，152，185和110 mm·a-1，累计沉降量分为186，181，200和124 mm；而沉降最为缓慢的是固安县，但沉降速率也达79 mm·a-1，漏斗中心累积沉降量为116 mm。李广宇等[1]利用2015年6月至2016年8月间的C波段SAR数据，监测到左各庄镇、通州区黑庄户、胜芳镇和廊坊市的沉降速率分别为163，159，197和108 mm·a-1，累积沉降量分别为191，187，234和120 mm，与本文监测结果基本相当。与李广宇等[1]文中采用的数据集相比，本文观测的时间跨度为2016年11月至2018年3月，新增了2016年9月至2018年3月的数据，一方面，验证了前人在研究区域监测结果的可靠性；另一方面，在一定程度上表明2016年9月至2018年3月间，研究区域的沉降变化趋势基本稳定。由图2可知，高铁线1沿线经过雄县、霸州及固安县等区域，高铁线2沿线河北段经过了廊坊市，高铁线3沿线北京段经过通州区，而以上分析表明，这几个区域均存在年沉降速率超过100 mm·a-1的大范围沉降漏斗，其中漏斗中心累积沉降量最大的是雄县，达200 mm。这些不均沉降可能会对高速铁路的安全运营及相关基础设施的稳定造成潜在的危害，应引起关注。

3.2 演化态势

为了进一步分析高速铁路沿线区域的不均匀沉降演变过程，以高铁线1为例，对铁路左右两边10 km范围内的沉降进行分析。图3给出了高铁线1沿线区域的年沉降速率。从图3可以看出，线路沿线的沉降漏斗主要集中在雄县至黄村段，沉降明显，范围分布较大，其中，雄县至霸州段不均匀沉降最为显著，漏斗中心相对的最大沉降速率达185 mm·a-1。为了分析铁路沿线的不均匀沉降变化，在铁路沿线选取9个PS点(图3中点a—k) ，对其进行时序分析。

图3 高铁线1沿线10 km沉降速率缓冲图

图4给出了9个PS点累积形变量的变化情况。由图4可以看出：从2015年11月8日至2018年3月27日，雄县至新机场段总体呈下沉趋势，沉降速率相对稳定；点a—f的累积沉降量均超过130 mm，其中点b的累积沉降量达200 mm，点g和h处的沉降速率缓慢，累积沉降量分别为38和68 mm；点a—k的累积形变量曲线中有4个明显的加速段 (如图4中粉红色框所示)，分别是2016年4月12日至2016年5月30日(图4中粉红框1所示) 、2017年5月31日至2017年6月12日(图4中粉红框2所示) 、2017年10月22日至2017年11月03日(图4中粉红框3所示)以及2018年1月2日至2018年2月7日(图4中粉红框4所示)。而这9个PS点中，除点h外，其余各点基本上都选在铁路附近，因此，这些不均匀沉降可能会对铁路安全运营及其基础配套设施的稳定造成潜在的危害，必要时需提前预防。

图4 沉降区域内所选PS点的累积沉降曲线

为了探究高铁线1沿线的沉降速率分布，图5给出了铁路沿线沉降速率。由图5可以看出：整体上，铁路沿线呈下降趋势；里程在0～15 km之间(李营至黄村附近)的沉降速率基本在6 mm·a-1左右波动；里程在15～55 km之间(黄村至固安县附近)的沉降速率出现3次陡崖式下降，位置分别在约DK23，DK30以及DK44处；里程55～85 km的沉降速率稳定在12～38 mm·a-1之间，里程85～95 km之间，沉降速率急剧下降，在DK98前后，出现陡崖式下降，沉降速率约达88 mm·a-1，而相应的地理位置是霸州至雄县段。通过距离量算发现，该结果与图3和图4的结果及地理位置基本对应。

图5 高铁线1沿线剖面沉降速率

为了更好地解释高铁线1沿线的沉降演变过程，还获取了该区域的水文地质资料。通过水文地质资料发现，固安及其附近区域的浅层水位埋深大致分布在15～25 m之间，霸州地区大致在28～33 m之间，而雄县地区大致在13～35 m之间。整个研究区域地下水埋深变化较大，水位季节性变化幅度在3～5 m之间，局部区域达7～9 m，而这些水位变化与不均匀沉降息息相关。表2统计了研究区域部分的集中供水井。由表2可以看出，供水井的供水面积普遍较广，深水井多，年供水量均在10万m3以上。据调查，李营至北京新机场路段两侧200 m范围内分布12 个集中供水井，59 个灌溉井，新机场至雄安路段两侧200 m范围内分布36个集中供水井，125 个灌溉井，由于浅部地层更松散，工程地质性质较深部地层差，多属中等压缩性地层，地层结构是黏性土层与砂层交互成层，为双面排水的地层结构，在相同的水位下降条件下，同样厚度的地层，浅部地层引起的局部不均匀沉降比深部地层更严重一些。通过以上资料分析，一定程度上解释了雄县至黄村段不匀沉降严重的原因。

图6给出了高铁线1—3沿线及周边的几个典型的沉降漏斗，并结合沉降漏斗周边的人口、交通、工业以及自然资源等对沉降漏斗演化进行分析。固安县位于北京市、廊坊市、天津市以及涿州市等中间，京九铁路纵穿全县，还有106国道、大广高速等横贯东西，北面有大兴机场，工业园区众多，流动人口较大，常住人口达52万左右，且该地区长期的天然气、地热以及石油等自然资源开采，是该地区沉降的主要原因之一。廊坊市的沉降主要集中在广阳和安次2个区。已有研究表明，2个区人口密度大，产业和水资源丰富[1，28]，地热和地下水的过度开采，导致该区域持续的地面沉降；胜芳和左各庄2镇工业园区众多，区域人口密度大，人类相关活动频繁，结合唐嘉等[26]研究，初步认为该因素为霸州市发生沉降的主要原因之一；雄县城镇人口达17.9万，城镇化率达47.82%，降水量为600 mm左右，天然气、石油、地热及地热水等储量丰富，该县原油年产70万t，天然气1 800万m3，人口的集中，自然资源的大量开采，是该区域沉降的主要原因之一。前已述及，9各PS点的累积沉降量均存在4个明显的加速段，而从这4个加速段的季节分布看，均在夏季或冬季，而夏季用水量巨大，冬季对地热和天然气等资源消耗较大，这在一定程度上呼应了以上分析。

表2 研究区域部分集中供水井统计表

图6 典型的沉降漏斗

4 结 语

通过对研究区域全局的沉降速率进行分析，发现严重的沉降漏斗主要分布在雄县、霸州市、固安县、廊坊市及北京通州区，最大沉降速率达206 mm·a-1。将本文的监测结果与前人相同区域的监测结果进行对比，发现沉降速率基本相当，表明研究区域的沉降态势趋于稳定。分析研究区域内3条重要的高速铁路沿线的沉降信息，发现它们经过了雄县、霸州市、廊坊市及北京通州区等这几个沉降严重的区域，漏斗中心的相对沉降速率均在100 mm·a-1以上，且累积沉降量最高达200 mm。利用高铁线1的沉降信息，发现雄县至霸州段不均匀沉降最为显著，漏斗中心相对的最大沉降速率达185 mm·a-1；分析雄县至新机场段的9个PS点的累积沉降速率发现，该段整体呈下降势态，且沉降速率相对稳定，但雄县至霸州段，累积沉降均超过130 mm，且各点在相同的时段出现4明显的加速沉降段，而4个加速段均在夏季或者冬季；从铁路沿线的剖面沉降速率来看，其结果基本与沿线的沉降信息及地理位置相对应。为了更好地解释沉降漏斗的演化态势，利用水文地质和相关的人类活动资料以及前人的相关解释对研究区域进行分析，结果表明各地区的不均匀沉降主要与区域本身的水文地质条件、自然资源的过度开采、区域人口高度集中、工农业活动频繁及交通情况相关。

以上研究和分析表明，研究区域内的不均匀沉降可能会对高速铁路的安全运营及相关附属建筑的稳定造成潜在的危害，需要密切关注。时序InSAR技术克服了传统InSAR技术的不足，能够精确，并全覆盖地探测出高速铁路沿线的沉降漏斗，这是传统监测技术无法相比的，但由于该技术不能进行定位，需借助外部DEM进行地理编码，完成雷达坐标到地理坐标的转换，其精度与DEM的分辨率密切相关；另外，星载SAR卫星是侧视成像，这势必会产生畸变误差，影响监测精度。因此，联合多源数据对高速铁路沿线进行监测，是一种必然趋势。