基于小基线集干涉测量的天津地区地表沉降时空分析

朱子林, 任 超,2, 周 吕,3,4*, 施显健, 李现广

(1.桂林理工大学测绘地理信息学院，桂林 541004；2.广西空间信息与测绘重点实验室，桂林 541004；3.武汉大学地球空间环境与大地测量教育部重点实验室，武汉 430079；4.城市空间信息工程北京市重点实验室，北京 100038)

城市地表沉降是当今世界各国城市化发展进程中广泛存在的一种缓慢的区域性地质灾害，该现象是由自然地质灾害(如地震、火山等)或人类活动(如地下水开采、建筑荷载等)引起的，城市建设、农业发展、防洪排涝等方面带来严重影响[1-3]，严重影响着人民的生命财产安全以及社会的可持续发展。天津市大部分地区由于大量的地下水开采以及工业和人口的迅速增长，从而导致地面发生严重沉降，监测整个天津地区沉降过程，保证人民群众的人身和财产安全是亟待完成的工作[4]。目前，对于地表形变监测和分析，传统的方法主要包括水准测量、GNSS(global navigation satellite system)测量等可以得到研究区域内精度较高以及高时间分辨率的形变量[5],但是这些方法的监测成本较高，监测点在整个研究范围内密度有限，而且其空间分辨率较低，很难有效地分析整个研究区域范围内的形变。因此，对城市地表形变进行大规模的观测，分析更为复杂的地表沉降现象是当前研究的热点问题。常规的合成孔径雷达差分干涉测量技术(differential interferometric synthetic aperture radar，D-InSAR)因其在地表沉降监测方面具有较高的空间分辨率，逐渐成为国际上公认的开展地表沉降监测的有效方法[6-7]。但是因为空间失相干和大气延迟的影响，传统的D-InSAR技术在应用过程中精度和可靠性会降低，为了减小上述影响，小基线集合成孔径雷达干涉测量(small baseline subset interferometric synthetic aperture radar,SBAS-InSAR)技术就是在常规的D-InSAR基础上新兴的一种技术，该技术不仅限制了时空失相干，提高影像的利用率，同时也丰富了影像的观测信息[8]。姜德才等[9]已经利用MCTSB-InSAR(multiple images coherent targets small baselines interferometric synthetic aperture radar)及其相关技术监测天津市地铁沉降的时序变化，并验证天津市地铁沿线沉降比较严重的区域；白泽朝等[10]基于Sentinel-1A数据利用PS-InSAR(permanent scatters interferometric synthetic aperture radar)技术监测天津市地表形变并确定天津市地表形变是地下水过度开采造成的；雷坤超等[11]利用PS-InSAR技术研究天津市地表沉降并确定其沉降区域有向周边延伸的趋势，但同时利用SBAS-InSAR技术和Sentinel-1A数据研究天津地区地表沉降的工作尚不多见。因此，利用2017—2019年的30景Sentinel-1A 数据，采用SBAS-InSAR技术监测天津市地表沉降，获取了天津地区的沉降分布情况以及沉降形变速率，通过与相同时间段内天津市二等水准测量数据对比分析，验证了SBAS-InSAR技术监测城市地表沉降的可行性，主要讨论分析了天津地区的地表沉降时空变化与地下水和城市化发展的关系，以期为掌握天津地区地表形变提供理论支撑。

1 SBAS-InSAR技术原理

SBAS-InSAR是由Bernardino等[12]和Lanari等[13]于2002年正式提出的一种时序InSAR的分析方法。此方法主要是将多景SAR影像进行一系列的排列组合，得到相关空间基线差分干涉图，同时可以连接被较大空间基线分开的SAR数据集，确保生成的像对中没有孤立的集群，以增加观测数据的时间采样率，可以更好地克服空间失相干现象。SBAS-InSAR技术主要使用奇异值分解法(singular value decomposition,SVD)来求解研究范围内的形变速率，进而获得高相干点的沉降形变速率和时间序列。其主要步骤如下[14-15]。

(1)在(t0,t1,…,tn)时间内，对获取的N+1景同一研究区域的Sentinel-1A影像数据，选择其中任意一幅影像作为超主影像，并将其他辅影像配准到超主影像上。假设每景影像都至少有一景影像与之干涉，则N+1景影像生成M景差分干涉图，M满足条件：

(1)

(2)对于任意差分干涉图j在辅影像tA时刻和主影像tB(tB>tA)时间内获得的影像进行干涉处理，方位向坐标x和距离向坐标r像素的干涉相位(忽略大气延迟相位、残余地形相位和噪声相位的影响)可以表示为

δφj=φB(x,r)-φA(x,r)≈

(2)

式(2)中：φ为干涉相位；j为差分干涉图的景号，j∈(1,2,…,M)；λ为雷达信号波长；φA(x,r)和φB(x,r)分别为tA和tB时刻SAR影像的相位值；d(tA,x,r)和d(tB,x,r)分别为tA时刻和tB时刻相对于d(t0,x,r)=0的视线方向上的累积形变量。

(3)为了更加清晰地表示地表沉降的时间序列，可以将干涉相位值用两个采集时间内的平均相位速率vj和该时间段之间的乘积来表示，即

(3)

则第j幅干涉图的相位值可以表示为

(4)

即在主影像时间间隔和辅影像时间间隔上每个时间周期上速度的积分，表达成矩阵的形式,即

Bv=δφ

(5)

在计算系数矩阵B(m×n)的过程中，因为SBAS-InSAR技术是采用多主影像时空基线的方法获取差分干涉图，所以矩阵B可能会出现秩亏的情况。矩阵B的广义逆矩阵可以用奇异值分解法(SVD)求解，从而得到速度矢量的最小范数解，最终通过对每个时间段的速度进行积分得到每个时间段内的形变量。

2 研究区概况及数据源

天津市是中国北方的一个沿海城市，部分地区被渤海环绕，该城市中心坐标约为117°20′E和39°12′N。由于大量的地下水开采、工业和人口的快速增长，天津大部分城市地区遭受了严重的沉降。水准历史资料显示，自20世纪80年代以来最大累积沉降已经超过1 m，自天津市开始实施限水政策控制后，沉降速率逐渐变慢，但是天津市中心城区以及周边城区仍然缺乏控制和沉降监测[4,16]。近年来，天津地区城市与轨道交通建设飞速发展，现在已基本形成完整的交通网络，新的沉降中心也在不断产生。为了进一步研究天津地区地表沉降情况，选择天津地区升轨模式下Sentinel-1A数据，由于原始影像数据量较大，后期工作量较大，处理起来比较困难，因此提前对原始数据进行裁剪，裁剪后的研究区范围如图1所示。

图1 研究区地理范围Fig.1 Geographical scope of the study area

研究选择了覆盖天津地区的30景升轨模式下Sentinel-1A影像数据[均为垂直发射水平接收(VH)极化方式]，时间跨度为2017年1月15日—2019年7月16日由欧洲航天局(European Space Agency，ESA)提供。精密轨道数据使用的是欧洲航天局提供的卫星精密轨道数据POD(precise orbit ephemerides)，用来提高配准和基线估算的精度[17]。此外，数据处理过程中所采用的DEM数据来自美国国家航空局(national aeronautics and space administration,NASA)提供的SRTM V4 DEM，分辨率为90 m。

3 数据处理与结果分析

3.1 数据处理

使用SARscape数据处理软件和30景Sentinel-1A影像数据监测天津市2017—2019年的地表沉降情况，主要处理流程有：①数据预处理：由于原始数据量较大，处理起来比较困难，因此将30景原始Sentinel-1A影像数据进行预处理，并通过裁剪得到本文研究区域；②生成连接像对：选取2018年5月10日获取的影像为超主影像，将其他影像配准到超主影像上，选取时间基线阈值为400 d和空间基线阈值为400 m，共生成139对小基线差分干涉图集，其相对空间基线连接图如图2所示；③干涉工作及相位解缠：将139对干涉像对进行干涉处理，主要包括生成干涉图，干涉图的去平、滤波和相位解缠以及相干性计算，将所有的数据对都配准到超主影像上；④轨道精炼和重去平：该步骤主要用于估算和去除不符合要求的残余恒定相位以及解缠后还存在的坡道相位；⑤生成时间序列结果：利用奇异值分解法对解缠后的相位进行解算，得到每个时段内的形变速率，最后通过地理编码得到每个时间段在视线方向(LOS)的形变量(负值表示沉降量，正值表示抬升量)。

绿点代表SAR影像；线段代表干涉图；黄点代表超主影像图2 SAR干涉图时空基线分布Fig.2 Temporal and spatial baseline distribution of SAR interferogram

3.2 结果分析

3.2.1 研究区域沉降速率分析

图3为天津市2017—2019年地表沉降形变速率。从图3可以看出：研究时间段内80%的研究区域沉降形变速率在-15.8～4.1 mm/a(正值表示地表上升，负值表示地表下降)。天津中心城区的沉降范围为-5.9～4.1 mm/a，表明天津市中心城区地表相对稳定，这与天津市在20世纪80年代以来施行市内六区严格控制地下水的开采政策有关，因此天津地区市内六区的地表小沉降相对较小且比较稳定；西青区以及胜芳镇(河北省霸州市所属)出现明显的沉降漏斗，其中王庆坨镇沉降中心较为严重，沉降中心的工业区沉降最为明显，年平均沉降速率达到-134.7 mm/a，最大累积沉降达到 324.6 mm，工业区以及居民小区都集中在该地区，工业用水以及生活用水都需要地下水，地下水开采严重以及地面负荷的增加是王庆坨镇沉降漏斗形成的主要原因；胜芳镇沉降漏斗位于河北省霸州市与西青区临近，该区域的郝家堡村沉降较为严重，年平均沉降速率达到-114.9 mm/a，最大累积沉降量达到230.8 mm，工业化的发展以及城市居民增加导致该地区地下水的抽取较为严重，过度的地下水开采是导致该地区出现沉降漏斗的主要原因；研究区域内静海县沉降区中团泊新城沉降相对比较明显，该区域年平均沉降速率达到-85.2 mm/a，最大累积沉降量达到167.4 mm，团泊镇被列为第三批国家新型城镇化综合试点地区，近年来城镇化发展较为迅速，新城突出“水和生态”逐渐成为多功能于一体的中国北方水上旅游城市，城镇化的迅速发展，城市扩张以及城镇化建设是该地区地表沉降的主要原因。

图3 天津市2017—2019年地表沉降形变速率Fig.3 Deformation rate of surface subsidence in Tianjin from 2017 to 2019

3.2.2 SBAS-InSAR结果与水准数据对比分析

研究区域内主要沉降中心位于王庆坨镇、胜芳镇以及静海县，与Li等[16]利用SBAS-InSAR技术获取的结果时空分布以及最大沉降区域(王庆坨镇和胜芳镇)基本吻合，验证了本文方法获取天津地区地表沉降形变的可靠性。为了进一步验证基于Sentinel-1A数据采用SBAS-InSAR准确性，选取研究区域内相同时间段内15个二等水准点数据(水准点 1～15分布如图4所示)，将二等水准测量数据与SBAS-InSAR相同时间段内计算结果进行比对，结果如图5所示。从图5可以看出，相同时间段内SBAS-InSAR计算的结果与二等水准数据结果吻合较好。对比分析两种数据结果可知，两种方法的最大相对误差为8.1 mm，最小相对误差为-1.54 mm，中误差为2.9 6 mm，结果表明两种方法得到的结果具有较好的一致性。

1～15为水准点位置图4 天津地区二等水准分布情况Fig.4 Second-class level distribution in Tianjin

图5 水准测量与SBAS-InSAR处理结果对比Fig.5 Comparison between leveling and SBAS-InSAR processing results

3.2.3 研究区地表沉降时空变化分析

利用时序特征点(特征点分布如图6所示)分析天津地区地表沉降能够较为直观地看出各特征点以及重点区域的沉降时空变化过程，选择四处不同沉降区域进行时序特征分析，结果如下：总体来看天津地区地表沉降不均匀较为明显(图7、图8)，天津市中心城区沉降相对稳定，王庆坨镇和胜芳镇出现明显的沉降漏斗，两地的最大累积沉降量分别达到324.6和230.8 mm，与这两个明显的沉降漏斗相比，静海县的沉降相对平缓，最大累积沉降为167.4 mm。由图7、图8(a)可以看出，王庆坨镇工业区沉降较为明显，出现明显的沉降漏斗，特征点沉降量成逐年增加趋势，研究时段内该区域的最大累积沉降达到324.6 mm，工业化发展以及地下水过度开采是导致沉降的主要原因；研究时段内胜芳镇沉降区的沉降随时间变化逐渐增加，其中郝家堡村沉降最为明显，最大累积沉降达到230.8 mm，该地区工业化发展较为迅速，工业和生活用水基本来自地下水的开采，过度开采地下水和工业的快速发展是该地区沉降的主要原因；静海县的团泊镇相较于王庆坨镇和郝家堡沉降相对较为平缓，最大累积沉降量为167.4 mm，团泊新城的发展迅速，城市化扩张是该地区沉降的主要原因；相较于三个沉降区明显的沉降漏斗，天津市中心城区沉降相对稳定，由图8(d)可以看出研究时段内沉降量随时间的变化有明显的起伏现象，中心城区发展相对较为稳定，同时中心城区施行限水政策严格控制地下水的开采，因此中心城区的沉降相对较为稳定。

1～15为时序特征点位置图6 时序特征点分布情况Fig.6 Time series feature point distribution

图7 沉降漏斗研究区Fig.7 Subsidence funnel study area

图8 各特征点累积沉降Fig.8 Accumulated settlement of characteristic each points

3.2.4 工业区沉降分析

王庆坨镇和胜芳镇这一带区域沉降较为严重，出现明显的沉降漏斗，并且沉降区域有向西延伸的趋势，针对两地沉降的情况，单独分析两地沉降形成的主要原因。图9(a)为王庆坨镇沉降区域形变速率情况，图9(b)王庆坨镇沉降区地表利用情况，该地区工业化发展较快，素有“自行车之乡”之称的王庆坨镇是中国北方最大的自行车生产基地，有260多家自行车厂家坐落于此，在工业化发展的同时该地区的城镇化建设达到了新水平，该地区的工业和生活用水均来自地下水，地下水的过度开采以及地面负荷的增加等因素是造成该地区沉降的主要原因。由图9可以发现，沉降最为严重的区域主要集中在工厂园区，年平均沉降速率最大超过-134 mm/a，最大累积沉降达到了324.6 mm，工业的快速发展带来了地下水的过度开采以及较大的地表空间利用率对地表沉降起着至关重要的影响。

图9 王庆坨镇沉降形变速率及地表空间利用Fig.9 Settlement deformation rate and utilization of surface space in Wangqingtuo town

近年来，胜芳镇经济持续发展，随着中国科学院中试基地和中国金属玻璃家具研发中心等国家级科研基地落户于此，胜芳镇被称为“中国金属玻璃家具产业基地”，胜芳镇工业发展相对较为迅速，主要有金属玻璃家具产业、钢铁产业、板材产业以及钢铁产业。由图10可以看出，胜芳镇沉降最为严重的区域位于郝家堡村，年平均沉降速率达到-134 mm/a，最大累积沉降量达到230 mm，钢铁产业和板材产业聚集于此，工业快速发展的同时带来了地下水的过度开采和地面荷载的增加，城镇化建设使地表空间利用率增大，地下水的过度开采和地面荷载的增加等因素是胜芳镇沉降的主要原因。

图10 胜芳镇沉降形变速率地表空间利用Fig.10 Settlement deformation rate and utilization of surface space in Shengfang town

2017—2019年，西青区和胜芳镇两个地区地下水过度开采导致地下水位基本处于下降状态。同时，由于工业化和人口增长的影响，两地的城镇化区域逐渐增加，因此天津地区地表沉降日益严重，且沉降区域有向西延伸的趋势，需要指出的是，胜芳镇和西青区之间的一大片区域受耕地和农作物的影响，沉降的连贯性缺失。

4 结论

研究采用SBAS-InSAR技术对天津地区30幅Sentinel-1A影像数据进行处理，揭示了天津地区2017—2019年的地表沉降时空变化特征。得到如下结论。

2017—2019年研究区域内出现明显的不均匀沉降，天津地区中心城区地表沉降相对平稳，平均沉降速率-5.9 mm/a。王庆坨镇和胜芳镇是天津地区两个最大沉降漏斗，年平均沉降速率分别达到了-134、-114 mm/a。此外研究区内还发现三个沉降较为明显的区域，即王庆坨镇和胜芳镇周边区域(-95.1～-65.3 mm/a)、这两个城镇之间的大片区域(-65.3～-55.4 mm/a)和静海县团泊镇(-104.9～-75.2 mm/a)；将SBAS-InSAR处理结果与同时期天津市二等水准数据进行对比分析，得到同一沉降趋势和分布。最大相对误差为8.1 mm，最小相对误差为-1.54 mm，中误差为2.96 mm，两种方法得到的数据结果具有较好的一致性，验证了利用SBAS-InSAR技术和Sentinel-1A数据用于城市地表沉降监测的可靠性；地下水的过度开采、工业化的迅速发展以及城镇化建设是天津地区地表沉降的主要因素，同时也是导致天津市地表沉降范围和幅度不断增大的主要原因。