基于小基线集时序干涉雷达分析武汉中心城区地表沉降时空特征

周 吕,鲁金金,文学霖,章 迪,王 成,杨 飞,

(1.桂林理工大学测绘地理信息学院，桂林 541004；2.桂林理工大学博文管理学院，桂林 541004；3.武汉大学测绘学院，武汉 430079；4.广西壮族自治区自然资源调查监测院，南宁 530023)

由于非人为(如地震、土层自然固结等)或人为因素(如地下水开采、地下建筑工程等)会引起地下松散地层的固结与压缩，进而导致地表沉降[1-2]。地表沉降是主要的区域性地质灾害之一，会对城市中的建筑物、基础设施、道路及桥梁等造成严重破坏，并影响人类的人身安全[3-4]。近年来，干涉雷达(interferometric synthetic aperture radar，InSAR)技术已经被广泛应用于研究地表沉降，相比于传统空间大地测量技术(如全球导航卫星系统、水准测量等)，该技术可以实现厘米至毫米级区域性尺度的地表沉降监测[5-6]。然而，散射体随时间产生的变化使采用InSAR技术时会存在一些问题，如时间、空间失相干，这些问题会导致信号失相干，降低形变监测精度[7]。为了克服InSAR技术存在的上述问题，提高监测精度，时序InSAR(time-series InSAR,TSInSAR)技术便应运而生了[8-9]。TS-InSAR通过同时处理不同时间获取的多景合成孔径雷达(synthetic aperture radar，SAR)影像来提取形变信息。目前，常用的TS-InSAR方法主要有PSInSAR(persistent scatterer InSAR)、SBASInSAR(small baseline subset interferometric synthetic aperture radar)、ISBAS(intermittent SBAS)及MTInSAR(multi-temporal interferometric synthetic aperture radar)等[10]。这些TS-InSAR方法已被广泛应用于不同城市地表形变监测，如北京[11-12]、墨西哥城[13]、河内[14]以及上海[15]等。

武汉作为中国中部地区的中心城市，由于快速的城市发展，近几十年来地面沉降严重[16]。武汉地区首次记录地面塌陷是1931年，这次地面塌陷导致了长江堤坝溃口，淹没了白沙洲地区[17]。自1978年以来，武汉市不同地区(如汉阳轧钢厂、烽火村以及青菱乡等)间歇性地发生了不同程度的地面塌陷[18]。近年来，由于地下水过度开采、地铁建设和岩溶塌陷等原因，在武汉形成了多个沉降区(如后湖与江汉沉降区等)，且这些沉降区的沉降范围与量级在逐渐扩大[19]。目前，武汉市测绘研究院在武汉市范围内布设了大约300个水准点监测武汉市地表沉降。然而，由于水准测量空间分辨率低、成本高，这使武汉市的区域性沉降信息和沉降漏斗的分布难以获得。因此，为了城市地质灾害的防治和可持续发展，需要有效的数据(如SAR数据)和方法来监测武汉的地表沉降状态和空间分布特征，掌握城市地面的形变规律。已有多项研究利用SAR数据，采用TSInSAR方法进行了武汉的地表沉降监测与分析，如Costantini等[19]利用2013年6月—2014年6月的45景COSMO-SkyMed SAR影像研究了汉口地区的地表沉降时空特征；白林等[20]利用27景TerraSAR-X影像(2013—2015年)提取了武汉的地表沉降信息，且研究时段内研究区域的沉降速率范围为-78.1～16.1 mm/a。然而，上述研究采用高成本的SAR数据提取了武汉地区的变形信息。采用Sentinel-1A IW SLC(single look complex)产品获取武汉中心城区沉降信息，同时依据水准数据、长江日水位变化数据、地铁建设分布数据以及降雨数据等数据研究区域地表沉降与各沉降影响因子(如城市建设、降雨量以及水位变化等)之间的相互关系。

首先基于水准数据与15景Sentinel-1A TOPS影像利用SBASInSAR时序方法获取武汉中心城区2015—2016年的地表沉降速率与沉降时间序列，并分析研究区域内的地表沉降时空特征；然后采用32个水准点数据验证SBASInSAR监测该地区地表沉降的精度与可行性；最后将长江水位变化、降雨量、地下水以及城市建设考虑为影响因子，研究其对地表沉降的影响机理，并分析地表沉降与影响因子间的相关性。

1 研究区域与数据源

1.1 研究区域概况

武汉市是华中地区的中心城市，位于江汉平原东部与大别山南坡。世界第三大河长江与其最大支流汉江流经武汉中心城区，将武汉分为三大主要区域，即汉口、武昌和汉阳。许多湖泊(如东湖、汤逊湖与梁子湖等)与河流分布于武汉市，这使得该市的水域面积达到2 205.06 km2，占武汉市总面积的25.79%。武汉地势为南低北高，北部为低山丘陵地带，中部主要是相对平坦的长江中游平原，南部被丘陵包围，其平均海拔约55 m。自1931年以来，武汉地区至少发生过15次岩溶地表塌陷灾害。由于地理位置的原因，武汉是中国重要的工业基地、科教基地和综合交通枢纽。近年来，随着经济的发展，大量的工程建筑物和地铁都在武汉地区施工建设，这导致武汉地区逐渐形成了多个地表不均匀沉降区(如后湖沉降区)。研究区域的地理位置如图1所示。

图1 研究区域地理位置Fig.1 Geographical location of the study area

1.2 数据源概述

现选取2015年4月—2016年4月期间覆盖武汉地区的15景升轨Sentinel-1A TOPS SAR影像(C波段)来估计垂直向平均沉降速率和沉降时间序列，SAR影像由欧洲空间局提供，在其官网下载。研究区域内的SAR影像的轨道方向为升轨，所有SAR影像的极化方式均为VV，且斜距分辨率为2.3 m，方位向分辨率为13.9 m。利用美国国家航空航天局(National Aeronautics and Space Administration,NASA)提供的SRTM3 DEM去除地形相位，采用ESA发布的精密轨道(precise orbit determination,POD)数据进行轨道精化与相位重去平。采用武汉市测绘研究院提供的32个水准点数据检校基于SBASInSAR获取的沉降结果，水准测量数据获取的时间跨度为2013—2016年，监测周期为1年。长江2015年4月—2016年4月期间的日水位变化数据由长江武汉航道局提供，这些数据和武汉地区的降雨量数据(数据在国家气象科学数据中心官网下载)被用来研究地表沉降与水位变化之间的相关性。

2 小基线集时序InSAR数据处理方法

利用SBASInSAR技术处理覆盖研究区的15景Sentinel-1A TOPS SAR影像获取武汉地区的沉降信息。首先设置合适的基线阈值将上述SAR数据生成合适的差分干涉图组合；然后估计出每一幅差分干涉图的形变信息，并把它们作为观测值；最后采用最小二乘算法或奇异值分解法(singular value decomposition，SVD)法求解形变速率，进而获得研究区域的沉降时间序列信息。主要的数据处理方法步骤如下[21-22]：

(1)差分干涉图组合生成：在生成差分干涉图之前，选择获取时间为2015年10月20日的SAR影像作为超主影像，并将所有辅影像配准和重采样至超主影像。依据空间基线小于300 m与时间基线小于200 d的原则选择干涉对，在此次分析中干涉对时间基线最长的为155 d。此外，对相干性低、解缠效果差的干涉对进行剔除处理，最后此次研究的差分干涉图组合中包含92对差分干涉图(如图2所示)。

图2 差分干涉基线图Fig.2 Differential interference baseline diagram

(2)平地相位与地形相位去除：采用ESA发布的POD数据去除平地相位，而地形相位利用NASA提供的90空间分辨率的SRTM3 DEM去除。

(3)轨道精化与相位重去平：在完成对干涉相位的自适应滤波与相位解缠(采用最小费用流算法)之后，从由武汉市测绘研究院提供的水准数据中选取稳定点对所有干涉对进行轨道精化与相位重去平处理。

(4)沉降速率与沉降时间序列提取：由图2可以看出，此次差分干涉图组合中没有孤立的干涉图集群，故采用最小二乘算法估计沉降速率。在计算沉降时间序列之前，首先将估计的线性沉降从原始沉降时间序列中减去，则余下的分量由大气相位、随机噪声相位及非线性沉降分量构成。由于大气相位和噪声相位分量具有高空间相关性，但时间相关性很低，因此可通过时空带通滤波将非线性沉降与它们分离。最后，将线性沉降分量与非线性沉降分量相加即可获取最终的沉降时间序列。

3 地表沉降结果分析及相互关系讨论

3.1 地表沉降速率结果与分析

图3显示了采用SBASInSAR时序分析方法获得的武汉中心城区的沉降速率图，该速率图叠加至2016年获取的武汉地区的Google Earth影像上。依据Sentinel-1A TOPS数据集通过SBASInSAR时序分析方法在研究区域内识别出的高相干目标点密度约为1 300个/km2。由图3可以看出，在研究时段内研究区域的年平均沉降速率范围为-81.5～18.2 mm/a，最大沉降速率位于汉口区的后湖地区。

由图3可知，研究区域内不均匀沉降明显，但郊区的地表沉降较小且稳定。研究区域内存在4个明显的沉降区域，即后湖沉降区、江汉沉降区、中山大道沉降区和武汉中央商务区(central business district，CBD)沉降区，且这些地表沉降分布与Costantini等[19]获取的结果相似。沉降最为严重的是后湖沉降区，该地区的最大沉降速率达到了-81.5 mm/a。江汉沉降区与中山大道沉降区的最大沉降速率分别达到-51 mm/a和-65 mm/a。此外，在研究时段内武汉CBD存在许多大型建筑物(如438 m高的在建武汉中心大厦)，这导致CBD探测到的PS点密度较其他区域稀疏(如图3所示)。CBD沉降区探测到的最大沉降速率为-59 mm/a。

沉降速率为正值表示垂直方向的地表抬升，沉降速率为负值则表示垂直方向的地表沉降图3 研究区域2015—2016年垂直沉降速率图Fig.3 Vertical deformation rates derived for the study area during the period from 2015 to 2016

3.2 时序InSAR结果精度评定

为了评估采用SBASInSAR时序分析方法处理Sentinel-1A TOPS数据来提取研究区域地表沉降速率的内符合精度，对沉降速率的标准差进行统计分析。图4显示了研究区域内所有监测目标点的沉降速率标准差的概率密度函数，目标点的最大标准差达11 mm/a，研究区域内86%目标监测点的标准差小于6 mm/a。此外，基于Sentinel-1A TOPS数据反演的沉降结果的标准差为3 mm/a。综合上述分析表明，依据Sentinel-1A TOPS数据采用SBASInSAR反演的地表沉降结果具有较高的可靠性和精度。

图4 平均速率标准差的概率密度函数Fig.4 The probability density function of the standard deviations of subsidence rates

为了进一步定量验证依据Sentinel-1A TOPS数据利用SBASInSAR时序分析方法进行地表沉降监测的精度与可靠性，故对SBASInSAR和水准测量结果进行了对比分析，此次对比分析采用了32个监测研究区域地表沉降的水准点(即BM1-BM32)。图3显示了上述32个水准监测点的位置分布。在对比分析之前，先选择时间上与Sentinel-1A TOPS数据相对应的水准测量数据。由于SBASInSAR识别的监测目标点与水准点在地理位置上不一定完全重合，因此针对上述现象，本节依据水准点附近的SBASInSAR识别的监测目标点采用Kriging方法内插出与水准点地理位置一致的监测目标点的沉降形变信息，两种方法的对比结果如图5所示。

由图5(a)可以看出，水准测量与SBASInSAR结果呈现出较高的一致性。假定水准测量结果为真值，图5(b)显示了SBASInSAR结果相对于水准测量结果的残差柱状分布，大部分差值分布在±10 mm/a以内。对两种方法的差值进行了统计分析，结果如表1所示。平均误差、最大残差及最小残差分别为-2.47、7.99、-11.09 mm/a，SBASInSAR时序分方法计算的沉降结果的均方根误差为5.82 mm/a。上述对比分析结果表明，采用SBASInSAR时序分析方法获取的结果与水准测量结果吻合得较好；依据Sentinel-1A TOPS影像数据利用SBASInSAR可以成功提取研究区域的沉降信息，且监测精度达5.82 mm/a。

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

表1 SBASInSAR与水准测量结果对比Table 1 Comparison between SBASInSAR and leveling-derived result

3.3 地表沉降时间序列与降雨量的相关性分析

为了分析地表沉降随时间的演变以及地表沉降与降雨量的关系，选取了分布在不同沉降区的6个PS点(即图3所示的PS1、PS2、PS3、PS4、PS5和PS6)，对研究时段内这6个PS点沉降时间序列和月平均降水量进行了对比分析，结果如图6所示。图6中6个PS点的沉降时间序列呈非线性下降，并伴随季节性变化。由于6个PS点分布的空间位置不一，它们的累计沉降量也不尽相同。位于后湖地区的PS1点处的累计沉降量明显大于其他点，这与后湖地区淤泥层的土壤承载能力较差有关。由图6可知，各PS点在2015年5月—7月间的沉降量非常小，而相应时段内武汉地区的降雨量却是全年中最大的时期。这些降水有效地补充了地下水，减缓了地表沉降速率。7月以后，降雨量急剧下降，同时随着武汉地区的气温与生活用水明显增加，导致地下水的抽取不断增加，进而引起地表沉降速率的增加。上述对比分析表明，研究区域的地表沉降伴随有季节性变化，地表沉降不仅与降雨量有关，还有与其他影响因素(如土层压缩性、建筑施工等)有关；降雨量是影响地表沉降季节性分量的主要因素。

图6 沉降时间序列与研究区域月平均降雨量对比Fig.6 Subsidence time series versus average monthly precipitation of study area

3.4 城市建设对地表沉降的影响分析

在研究时段内研究区域存在大量的在建的城市建设工程，如地铁6号线、武汉世贸中心大楼(位于图3的CBD沉降区)等。在深基坑和地下结构工程的开挖过程中，地下水位往往会出现高于施工作业面的情况，为了确保施工顺利进行，避免水下作业，需要对基坑进行排水。在研究区域内，整条地铁6号线与CBD地区的许多建筑物均在同时施工，因此在施工过程中抽取了大量的地下水，地下水位逐渐下降。这个现象降低了土层中地下水对地面上建(构)筑物的浮托力，导致土层压缩而引起地表沉降。研究区域主要位于由长江、汉江及湖泊共同作用而形成的冲积平原，故该地区土层的压缩性较高，易诱发地表沉降。此外，为了进一步验证计算结果，对后湖沉降区、中山大道沉降区和武汉CBD沉降区进行了野外实地调查，图7是现场拍摄的由不均匀沉降导致的某两栋楼的裂缝。通过野外实地调查进一步说明了研究区域存在明显的不均匀沉降现象。上述分析表明，大量的城市建设在研究区域的地表沉降中起主导作用；同时，地表沉降也受土层压缩性的影响。

图7 研究区域内某两栋楼的裂缝实景图Fig.7 Real view of the cracks of two buildings in study area

3.5 地表沉降与长江水位变化的比较分析

为了分析地表沉降与长江水位变化的相关性，选择沿长江边上分布的4个特征点(图3中的T1、T2、T3和T4)的沉降数据和由汉口水位站(图3中的红色方形)观测的长江日水位变化数据对比分析它们的变化时间序列。由于重点讨论分析地表沉降季节性变化与水位变化的关系，在对比分析地表沉降与长江水位变化的相关性之前，首先采用最小二乘方法去除地表沉降变化与水位变化中的线性趋势，去除线性趋势之后两者的对比结果如图8所示。由图8可以看出，这些特征点(T1、T2、T3和T4)的沉降变化呈现出一定的周期性特征；去趋势后的水位变化与沉降变化中呈现的峰值具有一定的相似性，但是沉降变化时间序列中的峰值比水位变化出现得要晚一些；形变曲线的波动幅度为3～4 mm。为了定量分析地表沉降与水位变化之间的相关性，采用了灰色关联分析方法来分析沉降与水位变化之间的接近度。接近度由灰关联度来描述，灰关联度是衡量不同离散时间序列数据间相似度的度量，灰关联度的数值越接近1，表示两离散时间序列之间的相关性越好。上述4个特征点去趋势后的沉降与水位变化时间序列的灰关联度均大于0.9，这说明图8中特征点的沉降与水位变化存在较大的关系。

图8 去趋势后的T1～T4地表垂直形变时序与长江水位变化对比Fig.8 Comparison between detrended surface vertical deformation time series relevant to the points (T1～T4) and detrended water level changes in Yangtze River

总体而言，长江两岸的地表沉降与水位变化相关，沉降时间序列的季节性分量受到水位变化影响的可能性较大。由于数据资料不充足，难以更为详细地研究沉降时间序列的季节变化。另外，较长的时间序列有助于进一步确定地表沉降与水位变化之间的相关性。因此，在未来的研究中，将利用更多的数据来研究信号的周期性特征。

4 结论

依据Sentinel-1A TOPS影像数据，采用小基线SBASInSAR分析方法反演武汉中心城区2015年4月—2016年4月的地表沉降时空分布特征，并采用水准数据验证。同时结合水准测量资料、长江水位日变化资料以及降雨量资料等，详细分析了地表沉降模式与人为活动(如地铁建设、大型建筑施工)和自然因素(如降水和水位变化)之间的相互关系。主要结论如下。

(1)研究区域内地表沉降不均匀明显，且存在4个明显沉降区。这些地表沉降区主要分布在后湖沉降区、江汉沉降区、中山大道沉降区以及武汉CBD沉降区。研究时段内武汉中心城区的沉降速率范围为-81.5～18.2 mm/a，最严重的沉降漏斗位于后湖地区。此外，武汉地区的地表沉降时间序列呈伴随明显季节性变化的非线性沉降。

(2)内符合精度检验结果显示，基于Sentinel-1A数据采用SBASInSAR获取研究区域沉降结果的标准差为3 mm/a，86%监测目标点的标准差小于6 mm/a。此外，SBASInSAR结果与水准处理结果吻合得较好，均方根误差和平均误差分别为5.82 mm/a和-2.47 mm/a。上述检验表明小基线集时序InSAR分析方法可有效地监测研究区域地表沉降，且具有较高的可靠性与精度。

(3)研究区域的地表沉降受城市建设影响严重；地表沉降的季节性变化与水位变化和降水有关。去趋势后的沉降与水位变化的时间序列的灰关联度均大于0.9，表明沉降变化与水位变化密切相关。但人为活动对武汉地表沉降的影响比自然因素更为显著。