地面沉降是由自然或人为因素引发地下松散岩层固结压缩造成的，属于缓变性地质灾害，具有不均匀、不可逆等特点，致使一定区域范围内地面高程降低的地质现象，会对国家基础设施及大型工程等造成不同程度的损坏，特别是对路基结构稳定性、连续桥梁及轨道平整性要求极其严格的高速铁路(葛大庆，2013；贾洪果，2012)。截止2009年，全国累计地面沉降量超过200 mm的地区达到8万km²，发生地面沉降的城市超过50个，其中京津冀地区是中国地表沉降较严重的区域，包含多条运营中的高铁(高文峰 等，2015)，累计经济损失达3328.28亿元(于冰，2015)。因此，京津冀及交通设施沿线区域急需范围大、精度高与自动化程度高等特点的时序雷达差分干涉TS-DInSAR沉降监测技术。

TS-DInSAR已被广泛用于地表沉降探测，除了经典的PSI (Ferretti 等，2000，2001)、SBAS方法(Berardino 等，2002)，还有干涉点目标分析IPTA算法(Werner 等，2003a，2003b)、StaMPS (Hooper 等，2004，2006)、时域相干点目标分析TCP-InSAR (Zhang 等，2011)等改进算法。TS-DInSAR算法的主要思想是根据多幅覆盖同一区域的SAR影像序列，选择PS点(Persistent Scatterer)，建立干涉相位模型，通过对模型解译获取沉降结果(Liu 等，2011)。目前，Werner等人(2003a，2003b)采用IPTA算法分别分析了C波段的ERS数据和L波段的JERS-1数据的时序应用，但均是应用于PS点密集的城镇小区域；孟秀军等人(2014)采用IPTA算法处理C波段的ASAR数据来探测秦岭断裂带地壳形变，但研究区域较小、地表覆盖单一。由上可知，在地面沉降调查监测中，TS-DInSAR受制于传统商业SAR卫星影像昂贵和卫星大范围成像能力较低，如何针对高现势性的影像数据准确提取大范围地表形变仍需要进一步深入。因此，本文研究成像能力较强且免费开放的新型Sentinel-1A数据在大区域沉降监测中的干涉点目标时序分析方法应用。

Sentinel-1A是欧洲太空局于2014年4月发射升空，可持续提供重访周期为12天的覆盖全球的开放数据。其中，TOPS成像技术获得的IW (Interferometric Wide Swath)模式影像覆盖范围为250 km×170 km，距离向分辨率为2.3 m，方位向分辨率为14.0 m，空前的成图能力足以满足大区域多时相InSAR处理和应用分析(Yague-Martinez 等，2015)。但是，在TOPS数据方位向上的9个脉冲序列(bursts)易发生多普勒中心不一致问题(Nannini 等，2016)，表现为方位向干涉相位跳跃，致使影像难以精确配准。

针对Sentinel-1A数据配准问题以及TS-DInSAR应用，吴文豪等人(2016)采用联合估计方法确定方位向burst之间的多普勒中心频率，并参考零多普勒时间对burst干涉图进行拼接；Yague-Martine等人(2016)利用轨道和外部DEM (Digital Elevation Model)进行几何粗配准的基础上，采用频移滤波、ESD (Enhanced Spectral Diversity)等方法提高配准效果；Dai等人(2016)采用强度互相关配准、ESD等方法进行迭代配准，并使用TSInSAR-AEM时序方法监测滑坡；Prats-Iraola等人(2015a)采用PSI方法分析了意大利那不勒斯地区和墨西哥城的时序应用；Sowter等人(2016)采用SBAS技术提取了墨西哥城市的地表沉降速率。

鉴于此，基于GAMMA软件平台，本文借助外部高精度的POD轨道和ALOS World 3D(分辨率30 m)，在强度互相关迭代配准的基础上，采用频移滤波和ESD结合迭代精化方位向偏移，从而提高Sentinel-1A空间几何配准的精度，并针对Sentinel-1A数据特征集成优化了基于点目标时序分析的大区域地表形变监测方法。根据过往资料，本文选取京津冀沉降较为严重的北京东、廊坊及天津西等区域作为研究区，覆盖面积约为1.02万km²，使用2015年6月—2016年8月29景C波段的Sentinel-1A数据，获取了大范围研究区域内最新的TS-DInSAR监测成果，并结合区域内人口分布、产业分布、地表覆盖以及线路剖面图等，深入分析了沉降成因和交通基础设施沿线的沉降时空分布，验证了本文针对Sentinel数据集成优化的一套影像配准流程和干涉点目标时序分析方法在获取大范围地表沉降结果方面具有可靠的应用精度。

TOPS模式Sentinel-1A数据的配准是在轨道信息和外部DEM的基础上，利用相邻burst之间的重叠区域信息，并以burst为基础单元进行的几何空间配准。为充分利用Sentinel-1A数据burst重叠区域信息，建议使用全幅数据配准后，再裁剪感兴趣区域进行下一步处理。Sentinel-1A SLC数据配准的具体流程如图1所示，其具体步骤如下：

步骤1 数据输入，生成S1 TOPS模式的SLC (Single Looks Complex)数据，利用burst之间和子条带之间的重叠区域信息生成mosaic SLC和mosaic MLI(Multi Looks Intensity)，由于S1 TOPS SLC的距离向分辨率与方位向分辨率悬殊较大，因此，和其他SAR数据不同，距离向和方位向的多视比较大(多视系数分别为10和2)。

步骤2 利用ESA (European Space Agency)提供的精轨POD数据，精化参数文件中的轨道信息，这有利于提高查询列表精度、初步减小卫星定时误差以及后续轨道误差对干涉相位的影响。

步骤3 借助外部DSM地理编码生成DSM模拟相位和mosaic MLI参数文件，生成主副影像之间的几何对应关系(查询列表)，并利用主副影像之间的偏移值改正来进一步精化查询列表。DSM相对于DEM来说，包含了建筑物等的高度，更有利于提高城市地区查询列表精度和地形模拟相位精度，从而提高几何配准的精度。

步骤4 以burst为基础单元，在主影像几何空间内，采用强度互相关方法迭代进行配准，拟合主副影像配准多项式，直至方位向偏移值小于0.02像素。

步骤5 为抑制定时误差造成频率失相干而引起的方位向相位跳跃和进一步精确burst在方位向的偏移值，本文提出结合频移滤波和ESD方法迭代直至消除方位向干涉条纹跳跃，能使得影像配准精度高于0.001像素。频移滤波是通过将主副影像的burst滤波到同一带宽，进而修正轨道定时误差引起的频谱失相关；ESD方法则是利用burst重叠区域之间的相位差分信息在强度互相关配准后进一步精确方位向偏移(Yague-Martinez 等，2016)。

图 1 Sentinel-1A数据配准流程图 Figure 1 Flowchart of Sentinel-1A images registration

干涉点目标时序分析方法属于TS-DInSAR方法。TS-DInSAR方法主要是基于已选PS点的差分干涉相位 $\varphi _{\rm {diff}} $ ，通过建立相位函数模型和采用合理算法，分离残余相位并获取地表形变信息。

式中， $\varphi _{\rm {def}} $ 是线性形变对应的差分干涉相位， $\varphi _{\rm {topo}} $ 是外部DEM引入的地形误差相位，残余相位 $\varphi _{\rm {res}} $ 包括基线误差造成的残余相位 $\varphi _{\rm {baseline}} $ ，大气相位 $\varphi _{\rm {atm}} $ ，噪声相位 $\varphi _{\rm {noise}} $ 以及非线性相位 $\varphi _{\rm {nldef}} $ 。

干涉点目标时序分析方法所采用的相位函数模型则是以差分干涉相位 $\varphi _{\rm {diff}} $ 为因变量，高程改正 $ \Delta h_{\rm{e}}$ 和线性形变速率 $ \Delta v$ 为自变量，残余相位为常量建立的相位回归分析模型，如式(3)：

式中，T为时间基线，B_p为空间垂直基线， $\lambda $ 、R、 $\theta $ 分别为雷达波长、雷达至地面目标的斜距及入射角，a和b为模型系数。

对相位模型式(3)在时空域进行回归分析，利用最小二乘方法求解相位模型中的高程改正值和线性形变速率改正值，更新对应的高程相位和线性形变相位，并通过逐步去除残余相位，对更新的相位模型进行迭代回归分析，从而得到与点目标观测值最佳拟合的形变速率和沉降值。

Sentinel-1A时序处理流程如图2，具体过程如下。

图 2 Sentinel-1A时序处理流程图 Figure 2 Flowchart of Sentinel-1A time series processing

(1) 预处理。经过全幅Sentinel-1A数据配准后，裁剪出感兴趣区域。根据SLC数据集和参数文件，生成基线文件和记录干涉对信息的文本文件。并使用Sentinel-1A数据配准时的地理编码流程，生成感兴趣区域的DEM模拟相位。

(2) 选取PS点。考虑到大区域内城镇和非城镇区域的同时存在，本文干涉点目标时序分析方法的初始选点密度不要低于500个/km²。为充分保证大区域内选点密度以及点较均匀分布，本文利用两种选点方法，一是基于时间域像素强度稳定性，二是基于影像点目标的光谱特性信息，通过合并处理两种方法获取的PS点，确保后续分块回归分析在参考点之间的传递。基于时间域像素强度稳定性是以均值与标准差的比值和时间域上的强度稳定性双重准则作为阈值进行选点，该阈值在影像数量超过20景时才有较可靠的统计；根据影像点目标的光谱特性信息是基于单景SLC影像计算的，可知该方法选出的PS点很有可能是移动的目标，后续可通过回归分析的离差阈值进行筛选(Werner 等，2003a，2003b)。

(3) 生成差分干涉相位集和相位初始回归分析。由于PS点选取精度、地形误差、基线误差和大气延迟等影响，生成的初始差分干涉图集中包含一定的残余相位，需在后续的回归分析迭代过程中逐步去除。初始回归分析成功的关键在于回归分析参考点的选取以及回归模型的参数设置。参考点应位于稳定、相干性好和PS点密集的区域。回归分析参数的设置是保证整个区域全部参与运算和解缠结果可靠的重要因素，遵循先大后小、单一变量原则，合适的参数将会增大回归模型的可靠性以及避免引入较大的传递误差和解缠误差。

(4) 残余相位解译和迭代回归分析。残余相位解译的目的在于改善差分干涉相位，优化回归分析模型。对差分干涉相位影响较大的残余误差包含地形误差、基线误差和大气延迟等。对于Sentinel-1A数据而言，TOS模式成像的非连续获取，导致其对基线误差更加敏感(Nannini 等，2016；Prats-Iraola 等，2015b)；同时POD轨道数据的厘米级精度仍存在一定的基线误差(Potin 等，2016)；在平常数据处理过程中，Sentinel-1A数据在大区域范围内普遍存在较明显大气延迟相位。因此有必要针对Sentinel-1A数据在大范围时序处理过程中进行基线精化和大气等误差相位分离。

基线精化是基于解缠的干涉相位和高程相位之间的最小二乘完成的，为尽可能获取足够高的精度，应遵循3个准则：(1)使用较稳定的PS点的解缠干涉相位，防止形变相位影响基线精化；(2)不要进行任何大气相位改正；(3)使用初始DEM模拟相位。尽管使用经过大气和DEM模拟相位改正的干涉相位会减小模型相位与DEM相位之间的标准差，但是也会引入基线精化之前基线误差造成的大趋势相位。在基线精化后，利用拟合的均方根进行检核，若值较小则代表精化成功。

残差相位中仍包含大气延迟、非线性形变和噪声等信息，这3种相位信息可根据在时空域的不同表现特征分离(Ferretti 等，2001)。大气延迟相位在大区域地表形变时序分析中具有较大影响，在空间域表现为低通且高相关，空间域差分的回归分析可以削弱大气对形变速率解算的影响；非线性形变相位在时空域一般是低通的；噪声信息在时空域均是随机高频的，因此可以通过滤波窗口为25的空间域最小二乘滤波器滤除噪声，在残余相位中分离出大气延迟和非线性形变相位，更新差分干涉相位并迭代回归分析。再通过时间域最小二乘滤波器从残余相位中分离出非线性形变相位(Wegmuller 等，2010)，将其和线性形变相位叠加获得最终沿视线向的形变时序结果。

本文研究区域如图3黄框所示，覆盖面积约为1.02万km²，包含了北京东部顺义区、朝阳区、通州区；廊坊市广阳区、安次区、经济技术开发区(国家级)、三河市、霸州市以及下属县镇；天津市市区大部分、静海区、西青区、北辰区、武清区等。该研究区域是京津冀一体化的核心区域，包含众多的公路、铁路等国家重要基础设施，具有人口密度大，经济开发区众多，道路网复杂且密集等特点。区域内的水文地质条件复杂，地表水资源匮乏，降雨量少，存在严重的过度地下水开采现象；地基主要为软土、松软土层，天津段基岩位于1000 m以下，这类土质具有含水量高、透水性差、压缩性高和强度低等特点。

针对本文研究区域，使用2015-06-24—2016-08-29(历时14个月)的29景Sentinel-1A数据作为实验的数据源，其距离向和方位向分辨率分别约为2.3 m和13.9 m，如表1干涉对和时间基线所示；选取2016年2月7日的影像作为主影像，其余副影像相对于主影像的空间基线最长为176.6291 m，最短为12.2570 m。根据本文提出的Sentinel-1A配准流程(图1)，所有副影像与主影像20160207的配准精度均高于0.001像素，且未出现相位跳跃现象。

图 3 研究区域 Figure 3 Interesting area

表 1(Table 1)

表 1 Sentinel-1A干涉对及其时空基线 Table 1 Temporal and spatial (perpendicular) baselines of Sentinel-1A interferograms 影像编号 成像时间 T/d Bp/m 影像编号 成像时间 T/d Bp/m 1 2015–06–24 –228 –106.1099 16 2016–01–26 –12 –145.3619 2 2015–07–06 –216 –176.6291 17 2016–02–07 0 0.0000 3 2015–07–18 –204 –70.8022 18 2016–02–19 12 –50.7977 4 2015–07–30 –192 –28.9481 19 2016–03–02 24 –130.9560 5 2015–08–11 –180 –114.3794 20 2016–03–14 36 –169.4663 6 2015–08–23 –168 –70.9494 21 2016–03–26 48 –72.0323 7 2015–09–16 –144 –162.8341 22 2016–04–07 60 –47.8485 8 2015–10–10 –120 –81.5029 23 2016–04–19 72 –76.9791 9 2015–10–22 –108 –161.5893 24 2016–05–13 96 –139.5867 10 2015–11–03 –96 –91.9303 25 2016–05–25 108 –122.7468 11 2015–11–27 –72 –12.2570 26 2016–06–06 120 –72.3288 12 2015–12–09 –60 30.1009 27 2016–06–30 144 –64.5355 13 2015–12–21 –48 –81.5897 28 2016–08–17 192 –67.4224 14 2016–01–02 –36 –91.3590 29 2016–08–29 204 –149.6802 15 2016–01–14 –24 –124.3575

表 1 Sentinel-1A干涉对及其时空基线 Table 1 Temporal and spatial (perpendicular) baselines of Sentinel-1A interferograms

根据2.3节所述的时序处理思想和流程基础上设置初始回归分析参数，依据强度阈值1.4和光谱阈值0.35选择的PS点分布如图4所示，数目为5557095个，初始点密度约为544个/km²；外部地形数据使用的是日本宇宙航空局研发的30 m空间分辨率的AW3D30 DSM；结合官方网站公布DSM的标准差为5 m，设置高程改正阈值为40 m；在保证绝大数点参与回归分析的基础上，初始拟合的沉降速率范围设置是±0.03 m/a，块大小设置为200像素，块间标准离差设置为0.6；结合Yan等人(2012)所做模拟实验可知，当块内标准离差低于1.1时，参与回归分析的PS点总是正确的，因此将块内标准离差设置为1.1。在逐步去除残余相位精化差分干涉相位后，由回归分析对形变相位进行改正的结果如图5所示，图5(a)是去除地形误差后线性形变相位的改正结果，图5(b)是去除基线误差后线性形变相位的改正结果，图5(c)是去除残余相位(大气、非线性和噪声)的改正结果，从结果可以看出，回归分析对去除地形和基线误差后的线性形变相位改正明显，由于空间域差分已经削弱大气影响，故去除大气等残余相位后的回归分析对形变相位已无明显改正，这也证明回归分析的解算精度达到最优。由于华北平原水平位移较小(姚宜斌，2008；张瑞，2012)，且较为关注的是垂直向沉降，因此可假定水平位移可以忽略，并根据式(4)将视线向(LOS)位移转化为垂直向。

图 4 PS点分布 Figure 4 PS distribution

图 5 线性形变相位改正 Figure 5 Linear deformation phase correction

最终垂直向沉降速率如图6所示，回归分析选择的参考点RP位于三河市，研究区域内整体沉降较为显著，其中较严重的沉降区有北京朝阳区及通州区、燕郊镇、廊坊市区及东冯务村、天津武清区、天津北辰区、天津静海区等，在这些沉降区的沉降中心位置选出以下沉降点金盏地区(a)、通州区黑庄户(b)、燕郊镇(c)、东冯务村(d)、廊坊市(e)、武清区(f)、北辰区(g)、王庆坨镇(h)、胜芳镇(i)、左各庄镇(j)和静海县区(k)，并绘制出上述沉降点的时间序列图(图7)。根据图7可知，a点沉降速率为138 mm/a，累积沉降量达到163 mm；b点沉降速率为159 mm/a，累积沉降量达到187 mm；c点沉降速率为104 mm/a，累积沉降量达到121 mm；d点沉降速率为161 mm/a，累积沉降量达到191 mm；e点沉降速率为108 mm/a，累积沉降量达到120 mm；f点沉降速率为98 mm/a，累积沉降量达到120 mm；g点沉降速率为106 mm/a，累积沉降量达到125 mm；h点沉降速率为224 mm/a，累积沉降量达到265 mm；i点沉降速率为197 mm/a，累积沉降量达到234 mm；j点沉降速率为163 mm/a，累积沉降量达到191 mm；k点沉降速率为111 mm/a，累积沉降量达到265 mm。

图 6 京津冀年平均沉降速率图 Figure 6 Subsidence rate per year of Beijing-Tianjin-Hebei area

图 7 沉降区内点对应时间序列图 Figure 7 Subsidence history of ground points in subsidence funnel

结合张永红等人(2016)、罗文林等人(2015)已有监测成果，根据本文最新的监测成果(如图8所示)可知，北京五环外的顺义–朝阳–通州沉降带的沉降速率正在逐步加剧，且通州区在2015年至2016年的最大沉降速率达到159 mm/a，其次是朝阳金盏地区最大沉降速率达到138 mm/a，同时该沉降带已扩展到三河市燕郊镇。北京地区沉降与人口密度有极大的关系，2015年北京人口普查显示，五环之外的人口占全市总人口的51%，达到1098万常住人口，北京功能扩展区朝阳区的常住外来人口密度占比位居第一；同时，因为北京市房价走高和京津冀一体化的发展，导致北京城市发展新区(顺义区、通州区、大兴区)和河北周边区域(如较发达的三河市燕郊镇)人口密度比较大(李国平和席强敏，2015)，从图8也可以看出，这些地方相对于其他区域的沉降较为显著，并结合杨艳等人(2012)和姜媛等人(2014)年研究可知，北京功能和人口疏解是解决地面沉降加剧的方法之一。

图 8 北京市年平均沉降速率图 Figure 8 Subsidence rate per year of Beijing

如图9所示的廊坊市区沉降速率图，广阳区沉降较安次区严重，这与两区域的人口密度、产业分布和水资源分布有关(易立新，2005)。广阳区人均密度为1133人/km²，而安阳区人口密度为623人/km²；广阳区产业主要以机械制造、建材、电子电器、食品、调料、服装、玻璃制品、印刷、化工等为主，且廊坊市经济开发区位于区内，安阳区则以文体乐器、印刷包装、精密电子配件为主导产业；广阳区较安阳区具有丰富的地热、温泉和地下水分布，这些产业和地下水分布的不同可能是导致广阳区地下水开采较安阳区更甚的原因，从而致使广阳区地面沉降较安次区更加严重。廊坊市东冯务村漏斗沉降较为严重，最大沉降速率达到161 mm/a，最大沉降量达到191 mm，并逐渐向周围区域发展。由图7廊坊市东冯务村e点对应的时间序列图和遥感影像图10(a)可知，e点及周围区域包含了大面积的竹柳、金花葵种植基地和食品饮料加工厂等，在2013年后，竹柳和金花葵种植基地实现机械化种植，竹柳属速生乔木，金花葵喜湿，这些可能致使地下水过度开采，导致地面持续沉降。

图 9 廊坊市年平均沉降速率图 Figure 9 Subsidence rate per year of Langfang

图 10 遥感影像 Figure 10 Remote sensing images

图 11 天津市年平均沉降速率图 Figure 11 Subsidence rate per year of Tianjin

如图11所示天津市和相邻河北区域的沉降速率图，并结合唐嘉等人(2015)计算2007年—2010年和张永红等人(2016)计算2012年—2014年的左各庄–胜芳镇–王庆坨镇沉降带(左胜王沉降带)结果分析，目前该沉降带正往东北方向发展，其总体沉降速率在125—244 mm/a的范围内，且武清区王庆坨镇沉降正逐步加剧。如图10(b)所示为2016年6月份的武清区王庆坨镇遥感影像，右下角为对应的沉降速率图，该地区在工业持续发展的同时，新建了尤张堡玉兰园项目、花卉苗木观光园等农业项目，包含辣木、油牡丹、树莓和黑森等需水量较高的高贵植物，该沉降区域毗邻居民区，并种植了大量需水较高的植物，这些因素可能导致该区域沉降在研究区域内最大且呈现加剧趋势。京沪和津保高速周围区域内沉降速率在120—220 mm/a，严重影响了两条高速的安全运营(张学东，2014)。同时，因为居民聚集区内人口和工业园区生产需水量的只增不减，左胜王沉降带正向东发展，并延伸至北辰区双口镇和青光镇，与武清区—北辰区沉降带(武北沉降带)连成一片。

综上，整体区域的沉降分布与以往资料相吻合，但根据本文结果显示，京津沉降区在逐渐连成片，且正在加剧。其中由于近期经济发展重心转移和市区房价走高等可能因素，引发城市郊区人口密度增加和工业园区增多，从而造成该区域沉降逐步加剧，天津武清区和廊坊市东冯务村的农业定位发展也成为了区域沉降加速的可能因素之一，尤其以王庆坨镇沉降最为显著。地区功能疏解和经济转型将可能是缓解研究区域内沉降危害的方法。

图 12 北京—天津交通基础设施沿线10 km沉降速率缓冲区 Figure 12 Ten kilometers buffer area of the Beijing-Tianjin line

图 13 北京—天津交通基础设施沿线剖面 Figure 13 Subsidence rate profile of the Beijing-Tianjin line

如图12和图13所示，分别是研究区域内北京—天津交通基础设施沿线区域及线路剖面沉降监测结果，结合两图可以看出区域内的不均匀沉降较为显著且主要发生在城市发展的边缘部分，研究区域内的北京沉降主要位于通州区，沉降速率在45—94 mm/a，该区边缘部分的沉降速率在20 mm/a上下；线路绕过廊坊市区且到该沉降区边缘的距离约为2 km左右；线路穿过天津武清区时，交通基础设施沿线周围的不均匀沉降较为严重，沉降速率在70—107 mm/a。结合图13交通设施沿线剖面结果中的红线框内局部放大图可知，铁路沿线沉降主要位于城市郊区，在北京东部的通州区的最大沉降量为94 mm/a，在天津郊区的最大沉降量为107 mm/a，而且郊区的边缘区域沉降梯度较为明显，这是交通设施安全运行的隐患区(Ge 等，2008)。郊区沉降较为严重的原因在于人口密度较大和工业园区较发达引发地下水超采；但相比于北京东部郊区，天津武清区域内沿线沉降较为严重的原因可能是武清区除工业园区外还有较为发达的种植业；相比于天津郊区，天津市区的沉降较为稳定，这是因为天津市区采取了地下水回灌措施(王威和陆阳，2016)，缓解了沉降状况。

如图14所示，廊坊—天津交通基础设施沿线周围区域的不均匀沉降较为显著，经过了廊坊市中心沉降区、天津武北沉降带和静海沉降区，并将其中的3个较为严重的沿线沉降区域进行局部放大。廊坊市内的交通设施沿线周围区域沉降主要以广阳区和安次区为中心向四周扩展，较严重的广阳区东冯务村(d)沉降漏斗边缘距离铁路线约3.3 km；天津界内区域沉降主要发生在北辰区和静海区，其中北辰区双口镇(S)和青光镇(Q)沉降漏斗已经与左胜王沉降带连接，最大沉降速率已达到200 mm/a，且该漏斗边缘距线路不到2 km，这种不均匀沉降会对线路工程安全造成很大威胁，同时由于线路经过两镇居民聚集地边缘，且两镇均大力发展蔬菜大棚种植和工业园区，导致周围区域沉降持续加剧；线路在经过静海区时，西面紧邻静海工业园区(A)，区域沉降速率在110—150 mm/a，东面距离最近的钢铁厂(B)等工业园区不到2 km，且该区域的沉降速率已达到150 mm/a。图15是交通基础设施沿线沉降速率剖面结果，廊坊市区内沿线的最大沉降速率为106 mm/a，天津郊区内沿线最大沉降速率为104 mm/a，影响安全运行的隐患区位于廊坊市区边界处和天津郊区边界及与市中心的交界处，其中廊坊市区的沉降防治可以借鉴天津市区沉降缓解的经验(熊隽，2013)。廊坊—天津交通基础设施沿线及周围区域沉降除受人口密度和工业园区影响外，大量需水作物的种植也呈现出重要影响，需要提前预防。

图 14 廊坊—天津交通基础设施沿线10 km沉降速率缓冲区 Figure 14 Ten kilometers buffer area of the Langfang-Tianjin line

图 15 廊坊—天津交通基础设施沿线剖面 Figure 15 Subsidence rate profile of the Langfang-Tianjin line

相对于传统的监测手段，TS-DInSAR方法在勘测区域地表形变方面具有极大优势。根据本文的监测结果显示：研究区域内分布着较多沉降区，不同的沉降带有连成片的趋势，城市和工业园区沉降正在向外逐步扩展；研究区内最大沉降位于王庆坨镇，主要原因包括该区域发展种植需水量高的经济植物以及毗邻用水较多的居民区和工业园区等；区域内沉降主要是由较密集的工业园区和人口致使地下水超采造成的，除此之外，城市郊区的沉降一定程度上受到大需水量作物种植的影响；交通基础设施沿线周围存在较多不均匀沉降区域，这些安全隐患区主要位于城市郊区和市区的边界区域。对比分析不同区域沉降的时空分布和成因可知，功能疏解、经济转型和地下水回灌可能是缓解沉降的方法。

在TS-DInSAR监测大区域地表形变中，面临着传统SAR影像覆盖范围有限和TOPS成像多普勒频率不一致的问题，本文针对新型Sentinel-1A影像序列开展TS-DInSAR分析，集成优化了一套配准流程和干涉点目标时序分析方法。重点针对多时相Sentinel-1A影像的配准问题，在强度互相关迭代配准的基础上，采用频移滤波和ESD结合迭代精化方位向偏移，有效抑制了burst方位向上的相位跳变并充分保证了配准精度；在针对大区域地表形变探测的时序方法应用中，探讨了回归分析参数的设置和处理流程的设计。

研究结果证明，Sentinel-1A数据采用本文所集成优化的数据处理方法在大范围地表形变监测领域中具有可靠的应用精度，可以用于监测基础设施运营的安全隐患区域。下一步研究工作将着重于探讨不同时序分析方法在监测大区域地表形变的应用，与本文监测结果进行交叉验证。

志 谢 欧洲太空局为本研究提供的Sentinel-1A数据和POD轨道数据，以及日本宇宙航空研究开发机构提供的DSM(AW3D)，在此一并表示感谢。