小基线技术在玉溪市的地面沉降研究

张航，左小清，李勇发，熊鹏，王志红，管庆丹，游洪

(1.昆明理工大学 国土资源工程学院，云南 昆明 650093； 2.贵州工程应用技术学院，贵州 毕节 551700)

1 引 言

地面沉降是自然和人为因素引发的，松散地层压缩导致的地面高程降低的地质现象[1]。近年来，由于我国经济迅速发展，城镇化进程加快，城市地面沉降速度明显加快，地面沉降给工农业、交通物流、城市建设造成严重危害，并使居民生活质量下降。目前，中国在21个省份中超过96个城市发生了不同程度的地面沉降，累计沉降量超过 200 mm的总面积超过7.9万平方公里，地面沉降问题愈发严峻，对城市进行地面沉降监测刻不容缓[2，3]。然而传统的城市地面监测方法如水准测量、三角高程测量、GPS(global positioning system)测量等，不仅存在工作成本高、观测周期长、观测过程中测点缺失等缺点，更不能进行大范围沉降监测[4，5]。合成孔径雷达干涉测量(interferometric synthetic aperture radar，InSAR)技术的出现，以其全天候、低成本、范围广、高精度等优势在城市地表形变监测中得到广泛应用，使得对于城市的大范围变形监测更加快捷、方便[6]。2002年SBAS-InSAR技术的出现，对D-InSAR技术中的大气影响、失相干等问题，可以通过多景雷达影像和小基线技术，有效降低时空失相干和大气因素的影响，以得到更加精确的形变结果[7，8]。

近年来，国内外学者运用SBAS-InSAR技术进行城市地面沉降监测，已取得一系列研究成果。2020年莫莹等利用小基线集技术对绍兴市上虞区进行了监测，获取了该地区的沉降信息，并利用同时期的水准数据验证了SBAS-InSAR监测结果的可靠性[9]。2018年谢文斌等利用SBAS-InSAR对昆明市进行了沉降监测，得到相关形变结果，昆明最大沉降速率为 -37.6 mm/a[10]。2020年熊鹏等对昆明市进行SBAS-InSAR双极化实验，其沉降区域与谢文斌所得到的区域沉降速率基本吻合[11]。2020年Shi Wei等利用了Sentinel-1A与ALOS两种数据对西安市进行监测，结果表明，西安市在2015年～2019年由于地下水位的上升，出现显著的局部回弹变形，最大抬升速率为 22 mm/a[12]。2021年冉培廉等采用SBAS-InSAR技术对西安市进行了沉降监测，结果显示，西安市的地面沉降趋势正在减缓，但其总体沉降依旧十分严重[13]。2017年JimingGuo等将SBAS-InSAR得到的天津市总体沉降结果与2015年地面沉降年报进行对比分析，结果表明，两者的沉降趋势和分布基本一致，并运用46个参考站数据验证了沉降结果的可靠性[14]。2020年张文哲等研究天津市SBAS时序沉降结果，发现与GNSS点观测结果相吻合，并依据《地面沉降干涉雷达数据测量技术规程》计算出两者中误差为 12 mm，满足技术规程要求[15]。因此，运用SBAS-InSAR技术进行地面沉降监测，可以得到较好的城市地面沉降效果。

目前，还未有学者利用SBAS-InSAR技术对玉溪市进行地面沉降监测，本文利用SBAS-InSAR技术适宜大范围地面沉降监测的优势，对玉溪市进行大面积沉降监测，可以克服传统监测方法的缺陷。本文基于SBAS-InSAR技术选取覆盖玉溪市中心区域的Sentinel-1A数据影像进行实验，以得到玉溪市地面沉降结果。并且结合水文地质因素分析玉溪市地面沉降产生的原因，为玉溪市的地面沉降防治提供相关依据。

2 研究区概况与实验数据

2.1 研究区概况

玉溪市是云南省第五大城市，地处云南省中部，玉溪市东南靠红河哈尼族彝族自治州，西南连普洱市，东北接省会昆明市，西北邻楚雄彝族自治州，总面积 15 285 km2，玉溪市地质构造属于扬子准地台，位于其西南边缘。若按二级构造划分，玉溪市东部为滇东台褶带，是元古代至三叠纪的盖层沉积，主要岩石为震旦系的砂页岩、碳酸盐岩，厚度达到了 800 m～2 600 m。并且由于断层陷落的影响东部存在抚仙湖、星云湖、杞麓湖三个大面积湖泊，其中抚仙湖是我国最大的深水型淡水湖，平均水深达到 95 m。玉溪市西部为川滇台背斜，地质年代为中晚元古界晋宁亚构造层，地面主要是昆阳群砂泥质板岩以及碳酸盐岩，岩层厚度在 10 km左右[17]。本文研究区范围如图1所示，主要选取红塔区、江川区全区，以及通海县、华宁县、峨山彝族自治县部分地区，其总面积 2 700余平方公里。

图1 研究区范围示意图

2.2 实验数据

本文实验数据，包括星载SAR数据，精密轨道数据，数字高程模型(digital elevation model，DEM)。欧洲航天局(European Space Agency)于2014年发射了Sentinel-1A卫星，该卫星包含一个用于数据采集的C波段合成孔径雷达，本文实验主要使用该卫星提供的SLC数据。Sentinel-1A数据为2018年4月～2020年11月范围内的69景降轨数据，影像主入射角为39.3°，距离向分辨率为 2.3 m，方位向分辨率为 13.9 m，极化方式采用垂直极化(VV)，详细信息如表1所示。精密轨道数据为哨兵数据对应轨道数据。DEM数据为美国航空航天局(NASA)的SRTM-1 DEM数据，平面分辨率为 30 m。

实验数据基本信息 表1

3 SBAS-InSAR数据处理

3.1 SBAS-InSAR原理

小基线干涉测量(synthetic baseline sub-set interferometry，SBAS-InSAR)的基本原理是选取多主影像，将多主影像和超级主影像组合出短时间基线、空间基线的干涉图，然后用奇异值分解(Singular Value Decomposition，SVD)法进行反演，以得到最终的沉降速率结果[16]，SBAS-InSAR技术路线如图2所示。大致原理如下：

假定某一研究区的影像有N+1幅，选取其中一幅作为主影像，将其余的影像配准到主影像的成像空间上，根据空间基线和时间基线的阈值组合干涉对，得到M幅差分干涉图，M数量与影像数量有以下关系(图2)：

图2 SBAS-InSAR技术路线图

(1)

对经过两个不同时刻tA和tB(tA>tB)获取的影像生成第j景差分干涉图，该差分干涉图中的干涉相位可以表示为：

(2)

上式(1)中的j代表任意一景差分干涉图，dtA和dtB分别为tA和tB时刻相对于t0的累积形变量。当忽略高程和大气等因素的影响时，可以简化为：

(3)

为了简化公式(3)，假定形变速率在相邻时间间隔以内是线性变化的，在整个的时间段内都是分段线性的，上式(3)的形变相位值可以写成：

(4)

将处理后的相位写成矩阵形式为：

Bv=δφ

(5)

式中：B为一个M×N的矩阵。在理想情况下同一小基线合内包含了所有影像，可以用最小二乘法直接进行求解以得到形变相位。但实际上这种可能性很小，因为在进行差分干涉图组合时，由于设置不同的时间及空间基线阈值，会出现在单个集合内采样不够，此时会使得矩阵B出现秩亏。为了得到最小范数意义上的最小二乘解，奇异值分解(SVD)方法可以解决秩亏问题，以求得最终的形变量信息。

3.2 数据处理

本文使用基于Envi5.3平台的SARscape 5.2.1软件进行处理。在数据处理之前，首先对哨兵数据进行预处理，将哨兵数据进行数据导入，之后根据本文研究区范围，进行裁剪。哨兵数据在进行导入、裁剪等预处理操作后，将进行以下五步处理操作，以得到研究区的沉降数据。详细处理步骤如下：

(1)生成基线连接图对输入的69景数据进行配对，生成的连接图会自动选择最优的组合方式进行配对，结果会出现时间、空间基线连接图。生成的像对会进行干涉工作流处理，用于SBAS反演。本文处理数据时自动选择2018年11月5日的影像作为超级主影像，所有的像对与超级主影像进行配准。由于本文时间跨度大，故时间基线阈值设置为100天，生成的空间连接如图3所示。

图3 基线连接图

(2)干涉处理软件生成干涉工作流，主要是对通过配对的干涉对进行干涉处理，将所有的数据对都配准到超级主影像上。在多视视数设置时，通过查询影像的多视视数，设置为5和1。滤波采用Goldstein方法减少噪声，Minimum Cost Flow方法进行解缠。这一步完成后，查看去平和干涉后的结果，发现干涉图总体干涉效果达到了预期。

(3)轨道精炼和重去平这一步主要是为了估算和去除残余的恒定相位以及解缠后还存在的相位坡道。在卫星不准确或者DEM定位不准确的情况下，用地面控制点(GCP)来修正SAR数据。为达到以上目的，主要通过选择一个去平干涉效果好的干涉图选取GCP，利用GCP进行重去平。选择GCP是在干涉图中无条纹，且平滑之处。

(4)SBAS反演分两次进行。第一次反演是SBAS反演的核心步骤，第一次估算形变速率地形，并且进行第二次空间解缠。处理完成后加载结果，符合预期结果。第二次反演是在第一次反演结果的基础上，通过时间域上进行高通滤波，空间域上进行低通滤波，还利用GCP去除残余的相位，得到更加准确的相位形变信息。

(5)地理编码该步骤是将处理得到的相位信息，转化为形变信息，并且将其转换到地理坐标系上，得到相关的矢量文件和栅格文件。

4 结果与分析

4.1 沉降结果

本文基于SBAS时序分析方法，通过对Sentinel-1SAR影像进行处理，得到玉溪市2018年4月～2020年11月的沉降结果，如图4所示，有七处明显的沉降区域(A～G)。其中红塔区沉降最为严重，出现四个明显的沉降区域(ABCD)，江川区出现两个沉降漏斗(EF)，通海县出现一处沉降漏斗(G)。七处沉降漏斗中A、C、E沉降漏斗呈片状分布，B、D、F、G沉降漏斗面积小，并且较为分散。本文统计沉降点的个数如表2所示，发现研究区内共 3 208 507个沉降点，其中超过85%的沉降点速率在 -13.58 mm/a～6.04 mm/a之间，表明研究区内总体保持稳定状态，沉降速率在 -13.59 mm/a以下的占比11.91%，抬升速率超过 6.04 mm/a占比仅2.07%，抬升和沉降的区域只占一小部分。

A沉降区主要以北城街道为中心，形成了一个巨大的沉降漏斗，并且其面积达到 30.1 km2，该区域总体的沉降速率在 -20 mm/a左右，累积沉降量最大为 -139.8 mm，最大沉降速率达到 -59.23 mm/a。该区域不仅包含北城街道的中心，还有其周围的梅园村、莲池村、刘张营等村落，不仅有建筑用地，还包括一大部分农业用地。分析其沉降原因，发现该研究区不仅城市扩张迅速，北城西侧为莲花池洪积扇，东北部是洪积裙及冲积台地，向内挤压造成了大范围的沉降。并有一定数量的深井开采地下水，导致北城地区地下水水位下降，引起该地区的沉降漏斗。

B沉降区官村-黄泥田一带，沉降点分散，累积沉降量最大达 -108.59 mm，最大沉降速率 -49.72 mm/a，沉降面积达 3.92 km2。有两处明显的形变区，沉降漏斗出现在官村-小犁铧营以及橄榄湾-黄泥田之间。不难看出，该区域西侧为莲花池洪积扇，受地质因素影响较大，判断主要因素为二元结构的以粗砾沉积为主的冲积物的影响。

C沉降区大营街-高仓呈现两片相连的趋势，尤其是杯湖村到上小屯之间，沉降现象非常严重，累积沉降量为 -117 mm，最大沉降速率达到 -72.83 mm/a。沉降区在中间的冲洪积层构成冲积台的西侧，而台面向西倾斜，造成表层黏土层发生沉降变化。随着该地区洗浴休闲为主的地热开发的蓬勃发展，已经造成了地热水资源的严重超采，地热水属于一种深层地下水，导致该地区地下水下降区域广，间接导致了地面沉降的发生。

D沉降区研和街道整体沉降区域分为两部分即潘井-和乐村与向家庄-南厂村，造成的沉降区域范围较小，沉降面积为 3.41 km2，但累积沉降量达到了 -146.04 mm，最大值沉降速率达 -67.27 mm/a。分析其沉降原因，研和街道地表岩层由昆阳群浅变质砂板岩组成，即第四系松散岩其厚度在 10 m以上，易造成沉降现象的发生。

E沉降区为江城镇中心城区，是一个较大的沉降漏斗，该地区位于地处星云湖北岸，抚仙湖西南，江城镇镇中心沉降尤为突出，沉降面积达到 7.08 km2，累积沉降量为 -142.38 mm，最大值沉降速率达 -68.67 mm/a。沉降漏斗主要为镇中心的城市建筑范围内，该处居民生活用水为抚仙湖水，地下水因素影响较小，该区域城镇化进程明显，分析为该地处地面载荷较大引起的地面沉降。

F沉降区位于星云湖南岸，该沉降区沉降范围分散，存在两个主要的沉降漏斗，即大寨河入湖区与烟草小区—鱼文化广场之间区域，总沉降面积达到 6.55 km2，累积沉降量为 -138.62 mm，最大值沉降速率达 -80.20 mm/a。分析其原因，大寨河入湖区初步判断为该区域常年被星云湖湖水浸泡，河流丰水期冲刷河道冲卓导致岸边土质松软，易导致地面沉降的发生。烟草小区—鱼文化广场之间，则是一处大型工地，建设过程中造成该地区的地面发生沉降。

G沉降区位于杞麓湖南岸。两处明显的沉降区域，即通海县杞麓湖一号别墅区、独房子村以及其周边区域。该区域沉降面积达到 3.83km2，累积沉降量为 -138.25 mm，最大值沉降速率达 -72.70 mm/a。杞麓湖一号别墅区，北面毗邻杞麓湖，东面有人工湖，地基受湖水影响，造成该区域的沉降较大。独房子村以及其周边地区，距离杞麓湖有一定区域，但由于湖水的浸泡，湖浪和沿岸流的冲刷和搬运作用形成各种侵蚀地形和沉积砂体，在农业生产时主要利用地下水进行灌溉，造成该处沉降的发生。

时序沉降点个数统计 表2

地面沉降是一个长期缓慢的过程，不易被人察觉，为了分析玉溪市地面沉降随时间变化的规律，提取各研究区沉降中心的沉降点(图4中1-7点)，形成如图5的时序累积沉降图。从图中可以明显看出，七处沉降区域沉降随时间的增加而逐渐加重。1(A)在2019年3月～2019年11月沉降加剧，沉降量达到 -40 mm，在2019年11月后，沉降趋势减缓。2(B)在2018年9月～2019年3月和2019年9月～2020年3月，整体保持稳定，分别在 -20 mm和 -60 mm上下波动，在2019年4月～2019年8月，急剧沉降 -30 mm。3(C)整体呈现线性分布，只在2020年5月～2020年7月沉降加速。4(D)在2019年2月～2019年7月，其沉降量达到了 -50 mm，其总体沉降量达到了 -132.55 mm。5(E)和7(G)两者沉降趋势大概一致，在2018年6月～2019年6月，两者差距在 20 mm，之后趋于重合，7(G)在2020年9月，急剧沉降，最后达到了 -129.92 mm。6(F)在2018年12月～2019年7月沉降速度加快，之后趋于缓慢，但在2019年7～10月与2020年5月数据出现异常，分析为两期数据去噪效果差。

图5 玉溪市典型沉降区累积形变序列

4.2 形变典型区域分析

为了进一步揭示沉降严重地区的垂直空间特征，本文选取A、C、E三个典型沉降漏斗区域，提取剖面线a-a′、b-b′、c-c′，位置如图6所示，得到沉降剖面线如图7所示。

图6 玉溪市三大沉降区域示意图

图7 沉降剖面线图

A区中北城横断面以大米饭馆为起点至春和居委会北部为终点，跨度约为 6 km。该区域出现了两个典型的“V”字形沉降漏斗(1)与(2)，分别位于陈大场村南部和高桥东南部，最大沉降速率分别为 -44.83 mm/a、-40.73 mm/a，该区域总体平均沉降速率在 -25 mm/a左右，并且两个沉降漏斗有进一步发展的趋势。C区中为大营街、高仓一带横断面从蔡官营开始至棠梨树东部结束，该区域有两个“U”字形沉降漏斗即(3)和(4)，(3)沉降漏斗的宽度约为 0.6 km，沉降漏斗沉降速率在 -35 mm/a以下，(4)沉降漏斗的宽度约为 1.4 km，该沉降区域沉降速率在 -30 mm/a以下。从两个沉降漏斗的发展来看，两个沉降漏斗的下沉趋势正在减缓。E区江城镇横截面起点自伊旗村东北与张官营交接处，终点在云岩村东北部，沉降区的跨度约为 4.5 km，该区域整体沉降速率在 -30 mm/a左右，在约 1.5 km处有一个明显的沉降漏斗(5)，江城镇财政所最大沉降速率达到 -55.99 mm/a，该区域存在两处沉降脊，沉降脊两侧沉降漏斗有进一步发育的趋势。

4.3 形变原因分析

现代社会每年用水量巨大，大致用途可分为农业、工业、服务业用水。而玉溪市内有珠江、红河等水系以及星云湖、杞麓湖、抚仙湖等湖泊，主要居民点的服务业用水以地表水为主，少部分以地下水供水。玉溪市的工农业用地占土地总面积的15%左右，工农业用水量超过总用水量的80%，地表用水已不足以供应其用水，绝大部分抽取地下水用以生产，造成地下水水位一直处于下降状态。2000年以来，而伴随着休闲服务行业的兴起，玉溪市内出现了多家以地热水为主的洗浴休闲中心，仅红塔区内处于开采中的地热井估计在10眼以上，年采热水超过122.6×104m3，将红塔区的沉降区域与张苗红[19]研究得到的地下水沉降区域进行对比发现，两者有较高的吻合度，表明地下水水位的下降是导致地面沉降的一个重要因素。

5 结 论

本文基于2018年7月～2020年11月Sentinel-1A的69景SAR数据，提取出玉溪市的沉降速率图。结果表明，有7处明显的形变区，最大累计沉降量为 -175.27 mm，最大年沉降速率达 -85.25 mm/a。北城沉降区沉降面积最大，达到了 30.1 km2。分析7处典型区域中心点的累积形变序列，发现在2019年3月～2019年9月沉降速度加快，之后沉降速度减缓。并对3个典型区域进行垂直空间特征分析，大部分漏斗的沉降速率在 -30 mm/a～-40 mm/a范围内，江城镇沉降漏斗最大沉降速率达到 -55.99 mm/a。通过进一步分析玉溪市沉降原因，发现玉溪市的经济快速发展导致的地下水位下降是主要原因。