城市地表形变问题是当今世界的一个常见问题，严重威胁着城市的基础建筑和人类的生命财产安全。引发城市地表形变的原因包括人为因素和自然地质因素。人为因素主要是由城市化的建设需要引起，自然因素则包括降雨、温度、地质等方面的影响。人类早在19世纪就开始桥梁的建设，自20世纪起，高层建筑设计逐渐兴起。到21世纪，地下空间的开发越发加快了城市化的进展^[1]。城市地表形变的现象也因此变得更加严峻，这严重影响了城市公共基础设施的稳定性，直接损害了当地的经济、环境的可持续发展^[2-3]。

监测地表形变是防灾减灾的关键。针对于大范围的研究区域，主要的技术手段有全球定位系统(GPS) 、合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar, InSAR)等方法^[4]。城市地表形变在空间上具有大范围的不确定性，传统的大地测量手段和GPS技术无法满足监测的需要^[5]。InSAR技术是具有大范围、全天时全天候监测能力的空间大地测量手段^[6]。InSAR技术可以利用雷达卫星的参数和空间几何位置关系，精确测量地面的空间位置^[7-8]。合成孔径雷达干涉测量差分技术(D-InSAR)通过差分手段可以获得某一位置上不同时间的微小形变^[9-10]。为了获得时间序列上的累计形变值，提出了时序合成空间雷达干涉测量技术(TS-InSAR)。其中包括Berardino在D-InSAR的基础上提出了小基线集合差分干涉测量技术(SBAS-InSAR)，以及Ferretti提出的永久散射体干涉测量技术(PS-InSAR)^[11-12]。TS-InSAR克服了常规InSAR中大气、失相干、DEM等误差，提高了时间空间分辨率，因此TS-InSAR技术已经广泛应用于城市地表形变的监测。

南昌市位于中国的江西省，处于赣江、抚河交汇处，是典型的冲击平原地貌。南昌市大致可以划分为3个区域，分别是中心城区、生态农业区、自然环境地貌区。中心城区是近年来城市化发展的主要区域，其中由于地铁施工运营导致的地表坍塌现象也出现过不少。另外，南昌市近年来年平均降水量逐渐增加，降水带来的地下水变化同样影响着地表的形变^[13]。目前使用InSAR技术对南昌市进行地表形变的研究有很多，大多数都是基于前人对当地地质资料的研究，结合研究处理结果的时空分布，判断形变结果的影响因素与程度。单从时空分布确认的影响因素存在不确定性，本文从周期方面考虑，进一步确定地表形变的影响因素。

本文基于SBAS-InSAR和PS-InSAR技术，利用Senitinel-1A卫星图像获取南昌市时间序列地表形变，揭示了自然地势、降水量、人类活动等因素对南昌市地表形变的影响，并结合外部资料对危险沉降区域的演化进行解析。研究分4个阶段进行：(1) 使用SBAS-InSAR和PS-InSAR技术获取长时间序列的地表形变值；(2) 对形变结果进行多角度相互验证，确保实验结果的可靠性；(3) 利用小波分析和灰色关联度分析评估降水量与地表形变的相关程度；(4) 结合外部资料，基于形变幅度较大的区域，解析该区域地铁在时间空间上的形变演化。

1 材料和方法 1.1 研究区域

南昌市位于江西省中部偏北，是中国长江中游地区重要的中心城市，其地理位置如图1所示。南昌市地貌属于赣江冲积平原，其地势西高东低。中部主城区多为软质的粉砂岩，土层性质较为柔软脆弱，东北地区相对平坦，西北地区丘陵起伏。作为连接海峡西岸经济区、珠江三角洲和长江三角洲的重要交通枢纽，南昌市近年来的城市化交通发展极为迅速。

1.2 数据来源

本文使用的数据有：欧洲航空局对地观测卫星系列Senitinel-1以IW扫描模式获取的单视角 SAR数据集(所使用的影像时间如表1所示)，如图1所示，覆盖南昌市地区需要两相邻的SAR影像拼接处理，轨道号分别是42287和42112；日本宇宙航空开发机构提供的30 m空间分辨率的ALOS World 3D 30 m DEM；POD精密定轨星历数据；天地图的光学影像，以及美国物理研究所提供的全球区域降雨数据。

1.3 技术方法 1.3.1 技术流程

研究主要分为3个部分：(1) 使用InSAR技术获取研究区域长时间序列形变结果，进行交叉验证证明方法的可靠性；(2) 利用灰色关联度和小波分析方法，评估降水量和沉降区域的关联度；(3) 基于以上数据，结合外部资料可视化危险沉降区域，解析其时空演化过程。本文从5个角度验证了InSAR结果的可靠性，评估了两种InSAR技术的适用性。此外，还探索了研究区域已有的危险沉降点，对大范围的城市地表形变进行初步安全排查。图2描述了研究的步骤。

1.3.2 灰色关联度分析方法

灰色关联度是可以在灰色系统中评估某个项目受其他因素影响的相对强弱的一种分析方法。灰色关联度分析前，对评估量进行无量纲化处理，减少数据绝对数值的差异。将各项目代入式(1)可求得灰色关联系数。将各序列的关联系数取平均值，即为该序列与参考序列的灰色关联度。

$ \begin{split} &{{\boldsymbol{\zeta }}}_{i}\left(k\right) = \\ & \frac{{{\rm{min}}}_{i=1}^{n}{{\rm{min}}}_{k=1}^{m}\left|{{\boldsymbol{X}}}_{0}\left(k\right) - {{\boldsymbol{X}}}_{i}\left(k\right) \right| + \rho {{\rm{max}}}_{i=1}^{n}{{\rm{max}}}_{k=1}^{m}\left|{{\boldsymbol{X}}}_{0}\left(k\right) - {{\boldsymbol{X}}}_{i}\left(k\right) \right|}{\left|{{\boldsymbol{X}}}_{0}\left(k\right) - {{\boldsymbol{X}}}_{i}\left(k\right) \right| + \rho {{\rm{max}}}_{i=1}^{n}{{\rm{max}}}_{k=1}^{m}\left|{{\boldsymbol{X}}}_{0}\left(k\right) - {{\boldsymbol{X}}}_{i}\left(k\right) \right|} \end{split} $

(1)

式中： $ {\boldsymbol{X}}_{0} $ 为母序列， $ {\boldsymbol{X}}_{{i}} $ 为各项目， $ {\boldsymbol{X}}_{{i}}\left({k}\right) $ 为各项目的序列，ρ为分辨系数，其系数越小，关联系数间差异越大，区分能力越强。

1.3.3 小波分析方法

小波分析基于合适的小波函数，通过小波变换计算小波系数，对信号逐步细化，得到信号的任意细节，进而分析信号在时间序列上的时频变化特征。本文利用小波方差，分析研究区域地表形变和降水量变化的周期成分。

本文使用的Cmor小波是复数形式的morlet小波，其表达式如式(2)所示。

$ {\psi }\left({x}\right) =\sqrt{\text{π}{{f}}_{\mathrm{b}}}{\mathrm{e}}^{2{i}\text{π}{{f}}_{\mathrm{c}}{x}}{\mathrm{e}}^{\frac{{x}}{{{f}}_{\mathrm{b}}}} $

(2)

式中： $ {{f}}_{\mathrm{b}} $ 为带宽参数， $ {{f}}_{\mathrm{c}} $ 为小波中心频率。

将母函数 $ {\psi }\left({t}\right) $ 进行伸缩和平移后，得到的小波序列如式(3)所示。

$ {{\psi }}_{{a},{b}}\left({t}\right) =\frac{1}{\sqrt{\left|{a}\right|}}{\psi }\left(\frac{{t}-{a}}{{a}}\right),\; {a},{b}\in {\bf{R}};{a}\ne 0 $

(3)

式中： $ {{\psi }}_{{a},{b}}\left({t}\right) $ 为子小波，a为伸缩因子；b为平移因子。

本文的时间序列是离散的数据，将其函数看作f (k $ \Delta $ t) ，则其连续小波变换系数如式(4)所示。

$ {\boldsymbol{W}}_{\mathrm{f}}({a},{b}) ={\left|{a}\right|}^{-\frac{1}{2}}\Delta {t}{\sum }_{{k}=1}^{{N}}\boldsymbol{f}({k}\Delta {t}) \overline {\boldsymbol{\varPsi }}\left(\frac{{k}\Delta {t}-{b}}{{a}}\right) $

(4)

式中： $ {\boldsymbol{W}}_{\mathrm{f}}({a},{b}) $ 为小波变换系数，f(k $ \Delta $ t) 为一个信号， $ \overline {\boldsymbol{\varPsi }}\left(\dfrac{{k}\Delta {t}-{b}}{{a}}\right) $ 为 $ {\varPsi }\left(\dfrac{{k}\Delta {t}-{b}}{{a}}\right) $ 的共轭复数。

将小波变换系数的平方值在b范围内积分，得到小波方差如式(5)所示，进而通过峰值查找信号的主要时间尺度。

$ {\boldsymbol{D}}\left({a}\right) ={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\left|{\boldsymbol{W}}_{\mathrm{f}}({a},{b}) \right|}^{2}\mathrm{d}{b} $

(5)

1.4 数据处理

SBAS数据处理所使用的时间基线阈值为120 d，常规基线阈值为45%，得到86对干涉对，其连接图如图3(a)所示。选择的相位解缠方法为Delaunay MCF。

PS-InSAR技术根据时间基线、空间基线、多普勒中心的因素影响自动选择2019年12月5日的影像为主影像，最长空间基线为116.612 m，其连接图如图3(b)所示。

在高程方面，由于SBAS-InSAR是基于面的处理，在相位解缠的时候会对高程进行平滑处理；而PS的分析对象是稳定的散射体像元，其保存了真实的高程信息，得到的高程精度较高^[14]。但是PS-InSAR 对数据集数量要求较高，至少需要20至30景以上的SAR影像，否则容易高估相干性，识别出虚假的PS点。南昌市新建区的大量翻新更对PS-InSAR技术的PS点识别造成影响。所以SBAS技术在本文的适用性会比较高一些，后续的时间序列分析使用的形变数据就是SBAS-InSAR的结果。

2 结果 2.1 InSAR结果

InSAR技术反演的南昌市雷达视线向的形变速率分别如图4和图5所示，其中速率的正值表示地表形变趋势为沿着雷达视线向上抬升，负值表示地表形变趋势为沿着雷达视线向下沉降。利用SBAS技术得到的形变结果的范围在−56.5~51.9 mm/a之间，其中形变范围在−16.5~9 mm/a区间内的点占99.4%，时间相干性在0.5以上的点的形变范围在−24.9~13.7 mm/a之间。而利用PS技术得到的形变结果的范围则在−24.3~24.4 mm/a之间。

2.2 准确性评估 2.2.1 时间空间一致性

图4和图5展示了两种方法所计算的南昌市地表形变具有空间上的相似性。值得注意的是，图5输出分辨率比图4低5倍，即图5若以图4的输出分辨率输出，则图5的PS点密度将更稀疏。其主要原因是在监测时间范围内，南昌市对新建区做了大量的旧城改造，PS-InSAR技术识别PS点时忽略了这些旧城改造的大片区域，而这些区域在SBAS技术上正是沉降量级较大的区域，因此图4和5会出现两种技术形变量级差别大的视觉误判。

为检验两种技术在量级上的差异，对比分析了两种技术的同名点位的时间序列形变，如图6所示，其中1号点位于28°38′34.53″N，116°11′59.27″E；2号点位于28°46′22.51″N，116°5′7.09″E。SBAS-InSAR和PS-InSAR得到的结果具有相近的累计形变量，但SBAS-InSAR的时序波动较大，具有更明显的周期信号。

2.2.2 误差分析

为进一步验证InSAR结果的可靠性。基于PS-InSAR与SBAS年平均形变结果的同名点，利用式(6)和(7)的中误差 $\delta $ 和平均误差 $\theta $ 两个指标进行精度评估。式中： $d_{{\rm{L}}_i} $ 为SBAS获取的形变结果， $d_{{\rm{I}}_i} $ 为PS-InSAR获取的形变结果，N为同名点的数量。得到中误差为±4.853 mm/a，平均误差为±3.847 mm/a。可见InSAR技术在相对精度上具有较好的可靠性。

$ {\delta }=\pm {\Bigg(\frac{{\displaystyle\sum }_{{i}=1}^{{N}}{\left({{d}}_{{{{\rm{L}}}}_{{i}}}-{{d}}_{{{{\rm{I}}}}_{{i}}}\right) }^{2}}{{N}-1}\Bigg) ^{\tfrac{1}{2}}} $

(6)

$ {\theta }=\pm \frac{{\displaystyle\sum }_{{i}=1}^{{N}}\left|{{d}}_{{{{\rm{L}}}}_{{i}}}-{{d}}_{{{{\rm{I}}}}_{{i}}}\right|}{{N}} $

(7)

3 讨论

SBAS和PS的结果表明，在研究时间范围内，南昌市地表形变出现不均匀的区域主要分布在如图7所示的A、B、C 3个区域，其中图7中的形变量是SBAS的结果。

3.1 季节性的水位变化

南昌市地处赣江抚河交汇处，是典型的河流侵蚀堆积平原。赣江的水位变化受当地降水量的影响。降水的增加使得地下水水位上升，进而降低了地下水与地面之间的应力，抬高了地表。降水量减少时，地下水补给不足，农业灌溉的需要使得地下水水位下降，增加了地下水与地面之间的应力，导致地面沉降。

降水量形变值数据是由美国物理研究所提供的，2019年1月至2021年1月的南昌站月平均降水量形变值如图8中黑色折线所示。图中展现了南昌市降雨情况呈现出季节性的周期变化，且2020年的年降水量形变值比2019年的大。

研究区域的A、B、C 3个范围内，分别选取了#1、#2、#3采样点，其地理位置如图7黑色三角形所示。将各采样点相邻两时间序列的形变值作差，得到每个时间点相对于上一时间点的变化量，如图8中橙、绿、黄3条折线所示。图8表明地表形变随降水量的变化而变化。2019年7月至10月，降水量大幅度下降，地表形变表现出下沉的信号。2020年1月至8月降水量大幅度增加，地表形变表现出抬升趋势。3个采样点的地表形变趋势都较为接近，说明区域A、B、C的地表时间序列形变趋势仍然具有很强的相似性，且与降水量的变化相关联。

为更进一步阐述地表形变与降雨量的关系，对#1、#2、#3处的时序地表形变进行无量纲处理，分别与降雨量进行灰色关联分析。灰色关联度越接近1，表明两时间序列相关性越好，分析结果如表2所示，其灰色关联度高达0.9以上。在去除地表形变的线性趋势的情况下，其与降水量的灰色关联度有所增加，说明区域A、B、C时间序列趋势的相似性与降水量有很强的相关性。

为了佐证灰色关联分析方法得到的灰色关联度的准确性，采用小波分析的方法，挖掘地表形变和降水量的周期信息。如图9(a)降水量的小波系数实部等值线图和图9(b)采样点#2的小波系数实部等值线图所示，降水量和采样点#2在第24个时间尺度附近具有最大时间跨度的周期性，而在第24个时间尺度上能量强度变化最大，所以第一主周期在第24个时间尺度上。同样的道理，在第9个时间尺度附近具有第二大时间跨度的周期性，而在第9个时间尺度上能量强度变化最大，所以第二主周期在第9个时间尺度上。在图9(c)的方差图中，降水量和采样点#2的第一极大值和第二极大值分别在第24和第9个时间尺度上，印证了主周期所在时间尺度的准确性。在图9(d)的第24个时间尺度的能量折线图中，采样点#2和降水量的周期都是16个时间序列，约为一年。在图9(e)的第9个时间尺度的能量折线图中，采样点#2和降水量的周期都是6个时间序列，约为5个月。

基于以上小波周期分析，采样点#2的地表形变和降水量具有相同的周期信号，佐证了灰色关联度分析的可靠性，证明了区域A、B、C时间序列地表形变受降水量影响。

3.2 地下水的影响

区域C是南昌市主要的生态农业区域，城市建筑群较稀疏，其地表形变由城市建设引发的可能性较小。该地生态农业区的灌溉主要通过抽取地下水来进行，因此地表形变可能是由地下水的水位变化引起的。

莫莹在2020年的研究中，利用2016年1月至2018年7月的Sentinel-1A，基于SBAS和PS技术，得到南昌市东南部生态农业区的形变表现为下沉信号^[15]。结合图10所示的2016年至2020年的降水量数据，可见2016年至2018年，南昌市降水量逐年降低。地表水的缺失以及农业灌溉的需要，导致在2016年至2018年期间，该生态农业区的地下水水位下降，从而引起地表沉降。

根据图11可知，2019年至2020年间，降水量的增加带来充足的地表水，可以持续补充到地下水中，缓和了之前地下水不足导致的沉降现象，所以区域C具有大范围的抬升信号。莫莹等^[15]计算的实验结果，符合本研究的结果和分析。综上所述，区域C的地表形变主要受到地下水位的影响，且2019年至2020年间，灌溉用水小于地下水的补充量，引起地表形变的回弹，表现为抬升趋势。

3.3 城市建设的作用

在研究时间范围内，南昌市的城市建设类型属地铁建设最为频繁，为保证地铁修建的安全性，本文作了南昌市地铁分布图，如图12(a)所示。从图12(a)中可以看出区域A、B是南昌市城市建设的主要区域，南昌市地铁运营线路都建设在这区域中。图12中的(b)(c)(d)(e)分别代表了地铁一号线、二号线、三号线、四号线周围1 km缓冲区范围内的形变速率图。

图12(a)展示了地铁沿线的沉降主要集中在赣江以北。结合中国地质科学院水文地质环境地质研究所的南昌市调查区第四纪地质与地貌图，可见赣江以北的黏土层较赣江以南厚，且近两年来南昌市的城市建设也主要集中在赣江以北的新建区。城市建设越频繁，地表的沉将越大。因此在地铁线跨河部分需要进行一系列的维护，防止不均匀的形变导致地铁设施断裂。

在赣江以北区域中，地铁四号线作为沉降量级较大的线路，其周边形变速率如图13(a)所示。图13(b)代表的是地铁沿线站点的累计形变值。从图13(b)可以看出，近两年的时间里白马山站、裕丰街站、璜溪站、中堡站以及礼庄山站的累计沉降量在−11.5~−15 mm之间，而在礼庄山站之后的西站南广场站、怀玉山大道站以及安丰站的累计沉降量在−5.5~−6.5 mm之间。礼庄山站与西站南广场站之间具有明显的不均匀沉降，两个站点及其周边的年平均形变速率图如图13(c)所示。

据报道，南昌市地铁四号线一期工程起点为白马山站，终点为鱼尾洲站，线路全长39.60 km，其中白马山站到礼庄山西段区间为高架线，长约5.49 km，其余为地下线，约33.91 km。地铁四号线高架段各区间陆续在2020年7至2020年10月完成箱梁工作。施工时间在本研究监测范围内，所以白马站到礼庄山站的大幅度沉降量可能与高架桥的施工有关。为了防止接下来地铁和高架桥的施工对地面沉降造成大幅度的影响，需要在礼庄山站到西站南广场地铁站直接增加维护结构，并且继续对这一区域进行地表形变监测。

4 总结

本文利用SBAS-InSAR技术和PS-InSAR技术，对南昌市2019年1月至2021年1月的25景Sentinel-1A雷达影像进行干涉，得到时间跨度内的地表形变速率图。从多种角度将两种技术得到的结果进行交叉验证，证明了结果的可靠性。满足相干性高于0.5的SBAS结果表明，南昌市地表形变的范围在−24.9~13.7 mm/a。

基于形变量级划分区域A、B、C，采用灰色关联度分析、小波周期分析方法，精细化评估了形变与降水量的关系。研究表明，A、B、C 3个区域其变化周期皆受降水量周期性变化影响。其中区域B形变量级较小，形变因素主要是降水量的周期性变化。

A区域赣江以北黏土层较厚，频繁的城市建设压缩了黏土层，引发沉降。区域C作为南昌市东南部的生态农业区，地下水的补充增加了其与土层的应力，引起了地表的抬升。

本文还分析了地铁沿线在沉降幅度较大的赣江以北的安全性，发现了若干不均匀沉降区域，可能会存在塌陷的危险。