地面沉降是当今世界各国城市化建设进程中普遍存在且不容忽视的环境地质问题之一。这种连续的、渐进的、累积式地质灾害过程缓慢且不可逆，一旦形成便难以恢复^[1]。引发地面沉降的主要原因是人类活动和地质作用，而过度抽取地下水及大规模市政工程建设是其中最为重要的诱发因素^[2]。地面沉降对地下和地面建筑物、构筑物产生破坏，威胁地铁、高铁等重大基础设施运营安全^[3]，给人们的生产、生活带来巨大危害，极大地影响和制约着当地经济建设可持续发展和社会安定。

武汉市地处江汉平原东部，地貌上属于残丘性的河湖冲积平原，软土分布广泛^[4]。近年来，作为长江中游城市群中心，武汉的城市化建设高速发展，大量大型基础设施工程及高层建筑的建设引起大范围的地面快速沉降，如汉口后湖地区、中山公园、武汉火车站周边等区域均有不同程度的地面沉降现象^[3]。然而，武汉目前缺乏长期、有效的地面沉降监测网，仅有少量工程施工期间的水准测量结果^[5-6]，难以进行大范围的地面沉降监测。

相对于水准测量等常规地面沉降测量方法，时间序列InSAR(time series InSAR)方法具有全天时、全天候、覆盖范围大、时间/空间分辨率高等优点^[7]，已广泛应用于地面沉降监测^[8-9]。目前已有相关研究利用该方法监测武汉地区的地面沉降：Bai等^[3]在2016年首次利用高分辨率TerraSAR-X数据监测武汉城区2009~2010年地面沉降，分析表明，城市建设和岩溶地质是沉降的主要诱因；Zhou等^[10]在2017年使用Sentinel-1A数据监测武汉地区2015~2016年地面沉降空间分布。2011~2015年武汉市有近万个工地同时开工，导致了严重的地面沉降及塌陷，如2015-02-10下午17时许，轨道交通3号线惠济二路至香港路站区间在施工过程中发生涌水涌砂险情，引起建设大道与黄孝河路交会处路面沉降达2 m，影响面积约400 m^2[11]。本文共收集到2013~2015年27景3 m分辨率的TerrsSAR-X数据(覆盖范围如图 1方框所示)，用于监测武汉市地面沉降情况。

1 时序InSAR数据处理

本文使用StaMPS-SB方法^[12-13]对武汉城区进行地面沉降监测。该方法使用基于干涉相位空间相关性的算法，识别滤波相位失相干缓慢(slowly-decorrelating filtered phase，SDFP)的目标点，能够识别相位稳定但幅度值较低的SDFP点，且不需要地面形变的先验知识^[14]。

1.1 InSAR数据预处理

首先根据时间、空间基线选择2014-02-07获取的影像为主影像，并将其余SAR影像与主影像进行配准重采样。为了抑制由于时间/空间基线过长引起的失相干，设置时间基线、空间基线阈值分别为110 d、180 m，共获取78个干涉对，时间-空间基线分布如图 2所示。为进一步提高相干性，分别在方位向去除不重叠的多普勒频率，在距离向进行频域滤波^[14]，利用30 m分辨率SRTM(shuttle radar topography mission)DEM^[15]去除地形相位，获取差分干涉图，并利用轨道数据和DEM对SAR影像进行地理编码。

1.2 SDFP点选取

为提高计算效率，首先使用幅度差离差法选取SDFP候选点。幅度差离差定义为：

${D_{\Delta A}} = {\delta _{\Delta A}}/{\mu _A} $

(1)

式中，δ_ΔA为主副影像幅度差的标准差，μ_A为平均幅度值。设置D_ΔA阈值为0.6，选取SDFP候选点。其次，对SDFP候选点干涉相位进行空间相关性分析，选取SDFP点。

在使用DEM和卫星轨道数据去除平地效应和地形相位后，第i个差分干涉对的第x个像素干涉相位为：

$ \begin{array}{l} {\varphi _{{\rm{int}}, x, i}} = W\{ {\varphi _{{\rm{def}}, x, i}} + {\varphi _{{\rm{atm}}, x, i}} + \\ \;\;\;\;\;{\varphi _{{\rm{orb}}, x, i}} + \Delta {\varphi _{\theta , x, i}} + {\varphi _{{\rm{n}}, x, i}}\} \end{array} $

(2)

式中，φ_{def, x, i}为视线方向形变引起的相位，φ_{atm, x, i}为大气延迟引起的相位，φ_{orb, x, i}为轨道误差引起的相位，Δφ_{θ, x, i}为入射角误差引起的相位，φ_{n, x, i}为噪声(像素内散射体变化、热噪声配准误差和方位向相位中心不确定性等)，W{}表示缠绕。

SDFP点相位应在长时间内保持稳定，其噪声需要足够小，为了识别SDFP点，需准确估算φ_n。

式(2)中，φ_def、φ_atm、φ_orb是空间相关的，Δφ_θ是部分空间相关的(即在空间域上是低频的)，因此空间域的低通滤波即可得到干涉相位中的空间相关部分。大气效应和形变速度会影响相位空间相关的距离，为了适应不同的相位梯度，在频率域进行自适应滤波，获取相位的空间相关部分$ {{\tilde \varphi }_{{\rm{int}}}}$。

干涉相位φ_int减去空间相关部分$ {{\tilde \varphi }_{{\rm{int}}}}$得：

$ \begin{array}{l} \;\;\;W\{ {\varphi _{{\rm{int}}, x, i}} - {{\tilde \varphi }_{{\rm{int}}, x, i}}\} = \\ W\{ \Delta \varphi _{\theta , x, i}^{{\rm{nc}}} + \varphi _{{\rm{n}}, x, i}^{{\rm{nc}}} + {\delta _{x, i}}\} \end{array} $

(3)

式中，Δφ_{θ, x, i}^nc为Δφ_{θ, x, i}的空间不相关部分，φ_{n, x, i}^nc为φ_{n, x, i}的空间不相关部分，δ_{x, i}为φ_{def, x, i}、φ_{atm, x, i}和φ_{orb, x, i}的非空间相关部分，假设其非常小、可忽略。

利用参数空间搜索和最小二乘反演联合估算Δφ_{θ, x, i}^nc，并将其从式(3)中减去：

$\begin{array}{l} W\{ {\varphi _{{\rm{int}}, x, i}} - {{\tilde \varphi }_{{\rm{int}}, x, i}} - \Delta \varphi _{\theta , x, i}^{{\rm{nc}}}\} = \\ \;\;\;\;\;\;W\{ \varphi _{{\rm{n}}, x, i}^{{\rm{nc}}} + \delta {\prime _{x, i}}\} \end{array} $

(4)

式中，δ′_{x, i}=δ_{x, i}+Δφ_{θ, x, i}^nc-Δ_{θ, x, i}^nc。

式(4)中相位主要为噪声，可定义像素噪声水平的测度参数γ_x作为SDFP点的判定标准：

$ {\gamma _x} = \frac{1}{N}\left| {\sum\limits_{i = 1}^N {{\rm{exp\{ j(}}{\varphi _{{\rm{int}}, x, i}} - {{\tilde \varphi }_{{\rm{int}}, x, i}} - \Delta \varphi _{\theta , x, i}^{{\rm{nc}}})\} } } \right| $

(5)

其中，γ_x越大，该像素为SDFP点的可能性越大。剔除γ_x较小的SDFP候选点，并将SDFP候选点的γ_x值作为下次计算的权重进行迭代计算，当2次迭代γ_x之差小于0.005时迭代结束，获取SDFP点集。

1.3 相位解缠与形变量估算

SDFP点的干涉相位是其真实相位的主值，取值范围在(-π, +π]，要得到真实相位需恢复其相位整周数。这种将相位主值恢复为真实值的过程称为相位解缠^[16]，是时序InSAR分析的关键步骤之一。在完成SDFP点识别与选取后，需对其干涉相位进行解缠，获取其真实值。为保证解缠正确，需要在解缠前估算并去除干涉相位中空间不相关部分，然后利用三维解缠算法进行相位解缠，获取SDFP点解缠相位^[17]。而形变相位则需利用时空带通滤波消除解缠相位中大气延迟、轨道误差、入射角误差及噪声等干扰项的影响，最后根据形变相位与形变量的线性关系，即可反演获取每个SDFP点的形变量时间序列及平均形变速率。

2 结果及分析 2.1 地面沉降整体趋势分析

2013~2015年武汉市地面沉降监测结果如图 3所示。在实验区域内分别识别并提取了3 255 386个高质量SDFP点，点密度达6 900个/km²以上，高密度的SDFP点能够对武汉市地面沉降的空间分布进行详细的描述。选择地理坐标为114.357 259°E、30.531 654°N(位于武汉大学IGS站附近)的SDFP点作为参考点，图 3中趋于蓝色标记的SDFP点表示沉降不明显或轻微抬升，而趋于红色标记的SDFP点则表明可能出现了明显下沉。从图中可以看出，2013~2015年武汉市整体比较稳定，地面沉降主要分布在汉口后湖、青山区工业路及光谷广场附近区域，沉降速率为5~78 mm/a。其中，后湖地区沉降范围最广，沉降速率最大，本文将详细分析其空间分布及成因。

2.2 地面沉降与软土分布的关系

从图 3可以看出，武汉后湖地区沉降最为严重，已发展成以竹叶山-解放公园为中心的沉降漏斗，最大沉降速度达78.1 mm/a。该地区历史上均为水田、湖滩，土质多为淤泥质软粉土或淤泥(图 4)，承载能力差，在人为抽取地下水及震动影响下，极易发生不均匀沉降。近年来，后湖地区大规模围湖造田，并修建了大量高层建筑与地铁，由于武汉地区地下水位较高，此举会影响地基承载力和边坡稳定性，在开挖基坑时，为确保施工安全，需采取排水措施，降低地下水位。从图 4(a)可以看出，时序InSAR监测结果的稳定区域与沉降区域分界线和1990年建筑用地与水田、湖滩的分界线(图 4(a)中红色虚线)十分一致，表明该区域地面沉降最可能的原因是软土地基上大规模城市建设引起的地下水下降和荷载增加。图 4(b)为该区域沉降中心3个SDFP点的沉降时间序列，可以看出，该地区沉降趋势明显，最大累积沉降量达86 mm。图 4(c)为后湖地区地面沉降引起建筑物损坏的实地照片，表明该地区地面沉降严重，验证了时序InSAR的监测结果。

2.3 地表形变与长江水位变化的关系

长江水位季节性变化剧烈，为研究长江水位变化对沿岸的影响，将沿岸SDFP点形变时间序列与长江水位进行对比分析(图 5)，结果表明，二者变化趋势较为一致。初步推断原因为沿岸地区地下水位与长江水位关系密切，丰水期长江水位上升，补给沿岸地下水，地下水位上升引起地面抬升；枯水期长江水位下降，沿岸地下水补给长江水，地下水位下降引起地面沉降。本文所用SAR数据时间跨度为2013-04~2015-03，长江水位完成2个周期的变化，地面形变表现为随长江水位上升或下降的同步变化，长江沿岸SDFD点形变时序表现为稳定信号(图 5(b))。

3 结语

利用时间序列InSAR方法对27景3 m分辨率的TerraSAR-X数据进行时序分析，获取2013~2015年武汉城区地面沉降分布。结果表明，武汉城区整体比较稳定，地面沉降主要分布在汉口后湖、青山区工业路及光谷广场附近区域，最大沉降速率达78.1 mm/a。结合地质、水文及实地踏勘资料分析表明，软土地基上大规模的城市建设是造成沉降的主要原因。另外，长江沿岸地表形变与长江水位变化密切相关，初步推测与长江水位对两岸承压水的动态影响有关。

本文填补了武汉城区2013~2015年地面沉降监测的空白，丰富了已有的地面沉降监测结果，为分析武汉城区地面沉降的长期时空演变特征和地面沉降的预防及治理提供了基础数据。随着Sentinel-1、ALOS-2卫星的相继发射，武汉地区的SAR数据愈加丰富，使得该地区定期开展时序InSAR地面沉降监测成为可能。