基于SBAS-InSAR的北京地区地表沉降监测与分析

周　吕　郭际明　李　昕　胡纪元

1　武汉大学测绘学院，武汉市珞喻路129号，430079 2　北京市测绘设计研究院城市空间信息工程北京市重点实验室，北京市羊坊店路15号，100038 3　桂林理工大学广西空间信息与测绘重点实验室，桂林市雁山街319号，541004 4　武汉大学精密工程与工业测量国家测绘地理信息局重点实验室，武汉市珞喻路129号，430079



基于SBAS-InSAR的北京地区地表沉降监测与分析

周吕1,2,3郭际明1,4李昕1胡纪元1

1武汉大学测绘学院，武汉市珞喻路129号，430079 2北京市测绘设计研究院城市空间信息工程北京市重点实验室，北京市羊坊店路15号，100038 3桂林理工大学广西空间信息与测绘重点实验室，桂林市雁山街319号，541004 4武汉大学精密工程与工业测量国家测绘地理信息局重点实验室，武汉市珞喻路129号，430079

运用SBAS-InSAR获取北京地区的地表沉降信息，采用18景ENVISAT ASAR影像完成北京地区2007～2010年地表沉降的时空分析。结果表明，北京地区沉降不均匀较为严重，在昌平区、顺义区、通州区等区域出现多处沉降漏斗，且有连成一片并向东扩张的趋势；大部分地区的平均沉降速率在-150 ～10 mm/a，沉降中心的最大沉降量超过400 mm；地表沉降受地下水开采与城市化影响明显。

地表沉降；SBAS-InSAR；北京地区；时空分析

对于地表形变监测与分析，传统的监测方法(如水准测量、GNSS测量、全站仪测量、分层标测量等[1])可以获取监测点较高时间分辨率与测量精度的形变量，但监测结果空间分辨率低，难以有效地监测与分析整个城市的区域性形变。合成孔径雷达差分干涉测量(differential interferometric synthetic aperture radar, D-InSAR)技术可以探测亚cm级的地表沉降，但其易受时间、空间失相干与大气延迟的影响，较难完成高精度的长时间间隔地表监测[2]。永久散射体合成孔径雷达干涉测量(persistent scatterer interferometric synthetic aperture radar, PS-InSAR)技术通过对高相干目标点进行相位分析，较好地克服了失相干与大气延迟影响，但需要较多的SAR影像数据[3-4]。小基线集合成孔径雷达干涉测量(small baseline subset interferometric synthetic aperture radar, SBAS-InSAR)技术将SAR影像数据组成若干个子集，采用最小二乘法求解子集的形变时间序列，同时利用奇异值分解法(singular value decomposition, SVD)将多个小基线集联合求解，获取整个时间跨度的形变序列，相对于PS-InSAR技术，SBAS-InSAR需要的SAR影像数目较少且获取非线性形变信息的能力较强[5-12]。

本文选取覆盖北京地区的18景C波段的ENVISAT ASAR 影像数据，采用SBAS-InSAR技术对地表沉降进行形变监测与分析，获取了研究区域2007～2010年的沉降分布情况与沉降速率场，并验证了SBAS-InSAR技术监测城市地表沉降形变的可行性。

1　SBAS-InSAR技术原理

假设获取了N+1景SAR影像，且影像获取时间t按时间顺序(t0,…,tN)排列,依据干涉基线组合可以生成M幅干涉图，且M满足：

(1)

假设干涉图j由tA时间获取的影像与tB时间获取的影像进行干涉生成(tB>tA),在去除平地效应与地形相位影响后，干涉图j中距离向为r与方位向为x的某一像素的干涉相位可表示为：

(2)

式中，φ(tB,x,r)和φ(tA,x,r)分别为tB与tA时刻SAR影像上的相位值，φdef,j(x,r)为tA时刻至tB时刻间卫星视线向的形变相位，φtopo,j(x,r)为参考DEM不精确引起的地形相位误差，φatm,j(x,r)为大气相位误差，φnoise,j(x,r)为噪声相位。其中，φdef,j(x,r)、φtopo,j(x,r)和φatm,j(x,r)可表示为：

(3)

式中，λ为雷达传播信号的波长，d(tB,x,r)和d(tA,x,r)分别是tB与tA时刻相对于参考时刻t0的雷达视线向的累积形变量，B⊥j为干涉图j的垂直基线，Δz为DEM误差，R为雷达与目标点之间的距离，θ为入射角，φatm,j(tB,x,r)和φatm,j(tA,x,r)分别为tB与tA时刻SAR影像中的大气相位分量。

为了获取研究区域的形变序列，需要精确估计出地形相位误差分量、大气延迟相位误差分量以及噪声相位分量，并将这3个分量从干涉相位δφj(x,r)中去除。

对于整个集合中的所有干涉图，在去除各项误差分量之后，由式(2)可以得到一个方程组系统，其矩阵形式为：

(4)

式中，A为M×N的系数矩阵，且∀j=1,…,M,M对应于干涉图数量，N对应于SAR影像数量，φT=[φ(t1),…,φ(tN)]为每一景SAR影像中高相干点对应的相位值所组成的向量,δφT=[δφ1,…,δφM]为各差分干涉图对应的解缠相位值所组成的向量。

为求解研究区域各高相干点的形变速率，可用两景影像间的平均相位速率代替相位值，则式(4)变为：

(5)

式中，B为M×N的系数矩阵，vT可以表示为：

(6)

当系数矩阵B为满秩(即M≥N)时，可用最小二乘法则求解出形变速率；而当M

2　研究区概况

北京市位于华北平原，地势较为平坦，平均海拔为43.5 m。北京市的平原地区主要是由永定河、潮白河、温榆河、泃河等河流联合作用而形成的山前洪积、冲积平原。该地区多处沉降区位于平原地区。研究区域如图1中的黑色方框所示，该区域西临北京西山，东临三河市燕郊镇，南临廊坊市，北临顺义区高丽营镇，区域内多为平原，地形起伏较小。研究区域的中心经纬度为39°51′N、116°28′E，面积约为4 013 km2。

图1　北京地区地理位置Fig.1　The geographical location of Beijing area

3　数据处理与结果分析

3.1数据处理

选取2007-08-01～2010-09-29的18景覆盖研究区域的ENVISAT ASAR 0级影像数据作为实验数据，其中影像数据采用波长为5.6 cm 的C波段，轨道方向为降轨，方位向与距离向的分辨率分别为4.570 m 与7.804 m，影像中心的入射角约为20.8°，极化方式为VV。采用美国宇航局提供的分辨率为90 m×90 m的SRTM DEM数据去除地形相位影响，同时利用欧洲空间局发布的DORIS轨道数据提高影像的轨道数据精度。

数据处理过程如下：1)将影像裁剪为本文的研究区域，对影像进行1×5(距离向×方位向)多视处理，选取2009-10-14获取的影像为公共主影像，将所有辅影像进行配准并重采样至公共主影像；2)选取时间基线阈值为400 d和空间基线阈值为900 m进行差分干涉图处理，共生成64对小基线差分干涉图集(图2，圆点代表SAR影像，线段代表干涉图，方框代表公共主影像)；3)采用DORIS轨道数据去除平地效应，同时利用SRTM3 DEM 数据消除地形相位；4)滤波处理，消除相关噪声；5)采用最小费用流法对小基线干涉图集进行相位解缠，结果见图3；6)依据振幅与相位的稳定性选择研究区域内的PS点，去除由DEM误差引起的相位分量，通过奇异值分解(SVD)法求解高相干目标点的沉降速率。

图2　时空基线结构Fig.2　The spatial-temporal structure of baselines

求解研究区的形变速率场时，选取区域内一个稳定点作为参考点R(依据已有的水准测量数据选取)。假定该点的形变速率为0，则各PS点的沉降速率均是相对于该参考点而言。获取研究区各PS点的沉降速率后，采用克里金插值算法计算出整个研究区域内的沉降形变速率场。最后，依据监测时段内速度的积分便可得到该时段的形变量。

3.2结果分析

图3为对各小基线差分干涉图进行解缠后的相位形变图，正值(蓝色)表示在视线方向上辅影像相对于主影像沉降，负值(红色)表示辅影像相对于主影像上升。从各时间段内的相位形变图可以看出，研究区内形变较明显，存在沉降漏斗。研究区内共识别出202 742个高相干目标点，平均每1 km2包含约51个高相干目标点。

图3　解缠后的小基线差分干涉图Fig.3　Small baseline interferograms stack after phase unwrapping

图4为2007-08-01～2010-09-29研究区的沉降形变平均速率图。可以看出，北京地区从北至南的主要沉降区分别为昌平区、顺义区、朝阳区、通州区和大兴区，与Ng等[13]采用PS-InSAR技术利用ALOS PALSAR 数据获取的形变结果一致，验证了本文采用SBAS-InSAR技术监测北京地区地表沉降形变的可行性与可靠性。

图4　研究区LOS向平均速率Fig.4　Mean LOS velocity of the studied area

为进一步验证实验结果的可靠性，收集研究区域内2010～2011年的12个水准点测量数据(水准点L1至L12的分布见图4)，将水准测量与SBAS-InSAR重叠时段的处理结果进行对比，见图5。可以看出，水准测量结果与SBAS-InSAR获取的结果趋势符合较好。对比分析两种方法的数据结果可以发现，两种方法的最大相对误差为 7.1 mm/a, 最小相对误差为 -0.4 mm/a，说明两种方法获取的结果具有较好的一致性。

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

沉降区主要分布于潮白河、温榆河和泃河流域的冲积、洪积扇平原上，昌平沉降区、顺义沉降区、朝阳沉降区与通州沉降区逐渐连成一片并有向东扩张的趋势，这与近几年北京地区的城市化扩张紧密相关。同时，在东部出现了新的沉降区即燕郊镇沉降区。

由图4可知，研究区内大部分地区的平均沉降速率在-150～10 mm/a之间，不均匀沉降情况明显。北京中心城区、海淀区、房山区、丰台区与石景山区的地表沉降量小且相对稳定，大部分地区的沉降速率小于10 mm/a，这与市区严格控制地下水的开采、北京地区西部的地表沉积物多为压缩性较小的粗颗粒沙卵砾石有关。昌平区出现多处沉降漏斗，其中沙河镇与上庄镇沉降中心较为严重，年沉降速率超过70 mm/a，且两处沉降中心的最大累计沉降量均超过180 mm。朝阳区的来广营乡、孙河乡和金盏乡一带沉降较为明显，其中金盏乡沉降区最为严重，该区域最大平均沉降速率超过130 mm/a，累计最大沉降量达到300 mm。

通州区的管庄乡、三间房乡和豆各庄乡一带(图4中的区域S)沉降较为严重。图6为通州区的沉降带区域S的沉降形变速率情况。图7为区域S浅地表空间的利用情况，该区域内地铁八通线、地铁6号线、通惠河沿东西向通过，并有多条高速公路与铁路穿过，同时其地下水开采较为严重，因此其沉降可能同时受动载荷、地下空间应用、地质构造以及地下水开采等影响。对比图6与图7可以发现，沉降中心主要位于同时有地铁、公路、铁路以及河流穿过的区域，且最大年平均沉降速率超过150 mm/a，最严重沉降中心的累计沉降量超过400 mm，不均匀沉降较为明显，导致该地区出现部分建筑物墙面开裂现象。说明复杂的浅地表空间利用对地表沉降具有一定的贡献率，同时地表沉降严重时会威胁建筑物的安全。

图6　S区域的沉降形变速率Fig.6　Subsidence deformation velocity in S area

图7　S区域的浅地表空间利用情况Fig.7　Shallow surface space utilization in S area

相对于通州沉降区、朝阳沉降区与昌平沉降区，顺义沉降区与大兴沉降区的沉降量较小。

2007～2010年，北京地区地下水开采量逐年增长，使得地下水位基本处于持续下降状态。同时，受经济发展、人口增加以及政策影响，北京的城市化面积逐渐增长，使得北京地区的地表沉降日趋严重，通州沉降区、朝阳沉降区、昌平沉降区等较为严重的沉降漏斗的地表沉降有逐渐连成一片并向东扩张的趋势。

4　结　语

本文采用SBAS-InSAR技术，利用18景ENVISAT ASAR 数据对北京地区的地表沉降形变进行时空分析，获取了该地区2007-08-01～2010-09-29的平均沉降速率以及各时段内的相位形变，并与前人的研究结果和实测水准数据进行对比，证明了本文方法的有效性。通过对研究区域的分析得出，2007～2010年北京地区沉降不均匀较为明显，并出现多处沉降漏斗，且各沉降漏斗逐渐连成一片并有向东发展的趋势，多处沉降中心的年平均沉降速率超过100 mm/a；较严重的沉降区主要分布在潮白河、温榆河和泃河等流域的冲积、洪积扇平原上，且部分地区的累计沉降量达到400 mm；地下水的过度开采与北京城市化的扩张严重影响地表形变的稳定性，使得地表沉降的幅度与范围逐渐增大。

[1]Poland M, Bürgmann R, Dzurisin D, et al. Constraints on the Mechanism of Long-Term, Steady Subsidence at Medicine Lake Volcano, Northern California, from GPS, Leveling, and InSAR[J]. Journal of Volcanology and Geothermal Research, 2006, 150(1):55-78

[2]Gabriel A K, Goldstein R M,Zebker H A. Mapping Small Elevation Changes over Large Areas: Differential Radar Interferometry [J].Journal of Geophysical Research, 1989, 94(B7):9 183-9 191

[3]Hooper A. A Multi-Temporal InSAR Method Incorporating Both Persistent Scatterer and Small Baseline Approaches[J]. Geophysical Research Letters, 2008, 35(16):96-106[4]刘媛媛, 张勤, 赵超英,等. PS-InSAR技术用于太原市地面沉降形变分析[J]. 大地测量与地球动力学, 2014, 34(2):6-9(Liu Yuanyuan, Zhang Qin, Zhao Chaoying, et al. Analysis of Ground Subsidence Deformation in Taiyuan City with PS-InSAR Technique[J].Journal of Geodesy and Geodynamics, 2014, 34(2):6-9)

[5]Beradino P, Fornaro G, Lanari R, et al. A New Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms[J]. IEEE Transactions on Geoscience & Remote Sensing, 2002, 40(11):2 375-2 383

[6]Casu F, Manzo M, Lanari R. A Quantitative Assessment of the SBAS Algorithm Performance for Surface Deformation Retrieval from DInSAR Data[J]. Remote Sensing of Environment, 2006, 102(3-4):195-210

[7]Chaussard E, Wdowinski S, Cabral-Cano E, et al. Land Subsidence in Central Mexico Detected by ALOS InSAR Time-Series[J]. Remote Sensing of Environment, 2014, 140(1):94-106

[8]Dong S C, Samsonov S, Yin H W, et al. Time-Series Analysis of Subsidence Associated with Rapid Urbanization in Shanghai, China Measured with SBAS InSAR Method[J]. Environmental Earth Sciences, 2014, 72(3):677-691

[9]Lanari R, Mora O, Manunta M, et al. A Small-Baseline Approach for Investigating Deformations on Full-Resolution Differential SAR Interferograms[J]. IEEE Transactions on Geoscience & Remote Sensing, 2004, 42(7):1 377-1 386

[10] Lanari R, Casu F, Manzo M, et al. An Overview of the Small Baseline Subset Algorithm: A DInSAR Technique for Surface Deformation Analysis[J].Pure & Applied Geophysics, 2007,164(4):637-661

[11] Usai S. A Least Squares Database Approach for SAR Interferometric Data[J]. IEEE Transactions on Geoscience & Remote Sensing, 2003, 41(4):753-760

[12] Usai S, Klees R. SAR Interferometry on a Very Long Time Scale: A Study of the Interferometric Characteristics of Man-Made Features[J]. IEEE Transactions on Geoscience & Remote Sensing, 1999, 37(4):2 118-2 123

[13] Ng A H M, Ge L, Li X, et al. Monitoring Ground Deformation in Beijing, China with Persistent Scatterer SAR Interferometry[J]. Journal of Geodesy, 2012, 86(6):375-392

Foundation support:National Natural Science Foundation of China,No.41474004,41461089; Fund of Beijing Key Laboratory of Urban Spatial Information Engineering,No. 2016204;Fund of Guangxi Key Laboratory of Spatial Information and Geomatics,No. 15-140-07-32; Open Fund of Key Laboratory of Precise Engineering and Industry Surveying, NASMG,No. PF2013-10.

About the first author:ZHOU Lü, PhD candidate, majors in InSAR data processing, E-mail: zhoulv_whu@163.com.

Monitoring and Analyzing on Ground Settlement in Beijing Area Based on SBAS-InSAR

ZHOULü1,2,3GUOJiming1,4LIXin1HUJiyuan1

1School of Geodesy and Geomatics, Wuhan University, 129 Luoyu Road, Wuhan 430079,China 2Beijing Key Laboratory of Urban Spatial Information Engineering, Beijing Institute of Surveying and Mapping,15 Yangfangdian Road,Beijing 100038,China 3Guangxi Key Laboratory for Spatial Information and Geomatics, Guilin University of Technology,319 Yanshan Street,Guilin 541004,China 4Key Laboratory of Precise Engineering and Industry Surveying, NASMG, Wuhan University,129 Luoyu Road, Wuhan 430079,China

In this paper, SBAS-InSAR is used to obtain high resolution ground subsidence information for the Beijing region. A spatial-temporal analysis of the ground subsidence in the region during the years of 2007 to 2010 is performed utilizing eighteen ENVISAT ASAR images. The results show that subsidence in the Beijing region is severely uneven; that multiple subsidence funnels formed in Changping district, Shunyi district, Tongzhou district, etc. are interconnected and have an eastward expansion trend; that the subsidence velocities in most areas are in the range of -150 mm/a to 10 mm/a and the maximum subsidence is over 400 mm; and that ground subsidence is influenced by groundwater exploitation and urbanization significantly.

ground subsidence; SBAS-InSAR; Beijing area; spatial-temporal analysis

GUO Jiming, professor, PhD supervisor, majors in precise engineering surveying and deformation monitoring, E-mail: jmguo@sgg.whu.edu.cn.

2015-09-18

周吕，博士生，主要从事InSAR数据处理研究，E-mail: zhoulv_whu@163.com。

郭际明，教授，博士生导师，主要从事精密工程测量与形变监测研究，E-mail: jmguo@sgg.whu.edu.cn。

10.14075/j.jgg.2016.09.009

1671-5942(2016)09-0793-05

P237

A

项目来源：国家自然科学基金(41474004，41461089)；城市空间信息工程北京市重点实验室基金(2016204)；广西空间信息与测绘重点实验室基金(15-140-07-32)；精密工程与工业测量国家测绘地理信息局重点实验室开放基金(PF2013-10)。