基于SBAS-InSAR 技术的安徽亳州市地面沉降时空分布特征与影响因素分析

2023-11-11何清魏路肖永红

中国地质灾害与防治学报 2023年5期

何清，魏路，肖永红

（安徽省地质环境监测总站，安徽合肥 230001）

0 引言

地面沉降又称地面下降或地陷，是指在自然环境和人类建设活动影响下，由于地下松散土层及岩层压缩固结，导致地表标高损失的一种地质现象[1]。造成该现象的自然因素和社会经济因素被认为是地面沉降的驱动力，主要包括水文地质条件、矿产资源开发及地下水开采状况等[2-3]。地面沉降形成原因复杂，防护与治理难度较高[4-6]，对人民的生命财产安全造成较大的损害。为深入探索沉降机制及其变化规律，国内外学者对沉降监测方法和驱动力因素进行研究，并不断发展新的理论技术[7-10]。

传统的地面沉降监测手段主要为水准测量和GNSS，这种局部单点测量的技术不仅成本高而且空间分辨率低，难以识别和监测大面积的地表形变[11]。时序合成孔径雷达干涉测量（time-series interferometric synthetic aperture radar，TS-InSAR）技术不仅周期短、精度高，还能够全天时、全天候地对大范围的地表形变进行监测，极大地弥补了传统监测手段的不足[12]。其中永久散射体雷达干涉[13]（permanent scatterer interferometric synthetic aperture radar，PS-InSAR）和小基线集雷达干涉[14]

（small baseline subset interferometric synthetic aperture radar，SBAS-InSAR）最具代表性。PS-InSAR 方法能够有效获取高相干目标（如建筑、桥梁及裸岩等）的时序形变信息，但在植被茂盛、稳定散射体稀少的区域，无法取得足量的稳定目标点，易使得形变解算结果产生偏差[15]。而SBAS-InSAR 方法是利用慢失相关滤波相位像素点获取地表形变信息，该类点能够在短时段内保持较强的相干性，且普遍存在于自然界中（如草地、裸土等）。因此SBAS-InSAR 方法比PS-InSAR方法更适用于大范围区域的形变监测[16-17]。

皖北地区的地面沉降问题历来较为突出，亳州市作为安徽省重要的新兴产业基地，其城市地质灾害监测一直受到政府及相关管理部门广泛关注。以往对该区域地面沉降的研究更侧重于观测数据的处理方法以及成因的简单分析[18-19]，对于其驱动力的量化研究尚且不足。探求地面沉降的主要驱动因素，能够为地质灾害防治和城市建设提供科学指导。本文以亳州市为研究区，选取2021 年10月至2022 年10月共62 景Sentinel-1数据，利用SBAS-InSAR 技术对亳州市地面沉降进行监测，分析亳州市地面沉降的时空分布特征，并基于地理加权回归（geographically weighted regression，GWR）模型，从地质环境、水文地质条件、人类工程活动和经济发展状况等方面对亳州市地面沉降的空间分异进行分析，探究亳州市地面沉降的主要驱动因素。

1 研究区与数据

1.1 研究区概况

亳州市位于黄淮海平原南端，皖、豫两省交界，全市下辖一区三县，总面积约8 522.58 km2（图1）。亳州市地处中朝准地台的淮河台坳二级构造单元，主要发育有褶皱、断裂构造。该地区地形起伏较小，地势西北高、东南低，辖境与黄河决口扇形地相连，总体呈典型的黄淮堆积型地貌[20]。地层属华北地层大区徐淮地层分区，第四系覆盖区内大部分基岩，第四系及新近系松散地层厚度在800～1 000 m。亳州市煤炭资源丰富，根据《亳州市矿产资源总体规划（2021—2025 年）》，区内现有煤矿产地17 处，主要分布在涡阳、蒙城等地，保有资源储量43.50 亿吨，占全省煤炭资源储量17.17%。主要含煤层为石炭、二叠系地层，煤层埋深600～1 000 m。

图1 研究区范围Fig.1 Study area scope

亳州市地下水类型可划分为松散岩类孔隙水、碳酸盐岩类裂隙溶洞水和基岩裂隙水三种类型。按照含水层的埋藏条件，可进一步划分为浅层孔隙含水层组（50 m 以浅）、中深层隙含水层组（50～165 m）、深层孔隙含水层组（165～660 m）、超深层孔隙含水层组（660～900 m）。根据《2021 年亳州市水资源公报》，亳州市地下水资源总量约为15.67×108m3，浅层地下水供水量4.59×108m3，中深层地下水供水量为1.47×108m3。全市域内已形成9 个超采区，总面积约980.9 km2，其中浅层地下水超采区开采量约0.4×108m3/a，中深层地下水超采区开采量约1.38×108m3/a。

1.2 数据

1.2.1 Sentinel-1

Sentinel-1 卫星是欧洲航天局发射的地球观测卫星，重访周期为12 天，具有干涉宽幅（IW）、超宽幅（EW）、波（WV）和带状图（SM）四种工作模式。本文选取2021年10 月至2022 年10 月间，共62 景升轨Sentinel-1 干涉宽幅（IW）模式的SLC 影像用于形变监测，数据的基本参数见表1。

表1 Sentinel-1 卫星数据参数表Table 1 Parameters of Sentinel-1 satellite data

1.2.2 SRTM DEM

SRTM 数据由美国国家航空航天局（NASA）和美国国家地理空间情报局（NGA）生产并面向全球用户免费发布，该数据覆盖了全球约五分之四的陆地表面，分辨率为30 m，高程精度为±16 m。本文中用于去除干涉测量过程中由地形起伏因素导致的地形相位。

1.2.3 其它数据

驱动力因子的选择主要依据研究区地质环境、水文地质条件、人类工程活动和经济发展状况四方面的综合影响，共选取8 个指标，分别为松散层厚度、中深层地下水埋深、深层地下水埋深、中深层水位变幅、深层水位变幅、道路密度、人口密度和单位面积GDP。

2 方法原理与数据处理

2.1 地面沉降监测

2.1.1 SBAS-InSAR 原理

将覆盖研究区的N+1幅影像进行配准后，参照一定的阈值组成M个干涉对[21]，则有：

假设第i（i∈1,2,···,M）个干涉对的主辅影像获取时间为ta和tb（其中tb在ta之后），并且其干涉相位中除去形变相位的部分已被剔除，则该干涉对的相位可以表示为：

那么M个干涉对的形变相位可以表示为如下矩阵形式：

每个干涉都可以产生一个观测方程，结合式（2）（3），可以组成具有M个观测方程的方程组，其中有N个待求未知数，矩阵形式方程组如下：

式中：A——M×N的系数矩阵。

若矩阵A的秩r（A）大于N，则可以通过最小二乘法求解式（4），公式如下：

式中：U——M×M阶正交矩阵，由AAT的特征向量组成；

S——M阶对角矩阵；

V——由ATA的特征向量组成的N×M阶正交矩阵。

式中：UT——U的转置矩阵；

A+、S+——矩阵A、矩阵S的广义逆矩阵。

将求解出时序形变量φ，除以形变所对应的时间间隔，即可求解出对应的形变速率。

2.1.2 数据处理流程

SBAS-InSAR 技术路线如图2 所示，获取地面形变信息的流程主要包括两部分：数据预处理和SBAS-InSAR工作流。本研究使用ENVI 平台的SARscape 对Sentinel-1数据进行处理，SARscape 是由sarmap 公司开发的一款专业的雷达影像处理软件，已被广泛应用于处理ERS-1/2、RADARSAT-1/2、ENVISAT ASAR、ALOS PALSAR以及 Sentinel-1（哨兵）等一系列星载雷达数据[22-25]。

3.3.2 其他个体因素除性别外，有研究发现年龄、职业和文化程度等对配偶间HIV传播也会产生影响［20，27］。

图2 SBAS-InSAR 技术路线Fig.2 SBAS-InSAR technical workflow

Sentinel-1 数据预处理步骤如下：①将数据导入为SARscape 的标准格式；②对同一时期两景SAR 影像进行镶嵌；③按照研究区范围对数据进行裁剪。对预处理后的影像进行SBAS-InSAR 处理，主要流程包括：①对输入数据以最优的组合方式配对；②配对后的像对进行干涉处理；③利用控制点对所有数据重去平；④去除大气相位并估算形变速率；⑤地理编码，将形变结果投影到地理坐标系上。

2.2 地理加权回归模型

空间关系具有异质性和非平稳性规律，为了对空间数据进行精确局部描述，Fotheringham 基于局部光滑的思想提出了地理加权回归模型（geographical weighted regression,GWR）[26]。GWR 实质上是一种空间变系数回归模型[27]，可以根据空间数据的位置信息生成对应的局部回归系数，从而对变量的局部空间关系与空间异质性进行合理的解释[28]。

2.2.1 模型构建

运用GWR 模型进行回归分析时，考虑到因子间的多重共线性问题会影响模型的可靠性，因此本文首先计算各因子的方差膨胀因子（VIF）。结果显示（表2），各因子的VIF 均处于0 到10 之间[29]，表明因子间不存在多重共线性。

表2 模型多重共线性检验Table 2 Model multicollinearity test

在此基础上，对亳州市地面沉降建立GWR 模型如下：

式中：yi——响应变量；

βj(ui,vi)——第i个样本点在(ui,vi)处的第j个回归参数；

xij——影响因素；

εi——随机误差项。

采用赤池信息准则最优带宽策略，构建不同运行模式下的活动强度GWR 模型[27]，结果如表3 所示。地面沉降GWR 模型的可决系数（R2）为0.394 583，表明自变量与因变量之间具有相关性。校正可决系数（adjustedR2）是0.373 125，说明可解释因变量在模型中具有较高比例。

表3 2022 年地面沉降GWR 回归模型参数Table 3 Ground subsidence GWR regression model parameters for 2022

3 结果

3.1 精度验证

为验证本文中研究区地面沉降数据的可靠性，选取谯城区周围3 个同期水准点测量值与SBAS-InSAR 监测结果进行对比。结果（表4）显示，SBAS-InSAR 结果与水准测量值的误差在1 mm 以内，说明监测结果具有较高的可信度。

表4 SBAS-InSAR 监测结果与水准数据对比Table 4 Comparison between SBAS-InSAR monitoring results and leveling data

监测结果与实测数据间存在一定误差，主要是因为SAR 影像在干涉过程中受到大气延迟、地形起伏和失相干等多种因素影响产生的误差。此外，水准测量获取的是单个监测点的高程变化，而SBAS-InSAR 结果则是一个单元格网（面状）的平均形变量，二者不一定完全对应。

3.2 地面沉降时空分布特征

通过SBAS-InSAR 处理，得到2021 年10 月至2022年10 月内亳州市地表形变速率（图3）。亳州市整体沉降速率为5～30 mm/a，平均沉降速率为5.7 mm/a。地面沉降主要分布于谯城区东北部、涡阳县城、利辛县城以及蒙城县的部分地区；沉降最严重区域位于涡阳县公吉寺镇以北，受煤矿开采影响，沉降速率幅值达到84.3 mm/a；谯城区东北侧的地面沉降幅值为25.8 mm/a；在利辛县城及蒙城县的部分地区内，大多数区域地面沉降速率幅度小于10 mm/a 水平，局部区域地面沉降幅度达到20 mm/a 水平。

图3 亳州市2021 年10 月至2022 年10 月形变速率分布图Fig.3 Distribution map of the subsiding rate of Bozhou from October 2021 to October 2022

监测时段内沿雷达视线向（Line of Sight,LOS）的时序累计形变量如图4 所示。在涡阳县中部、谯城区东北部、利辛县西部以及蒙城县中部，均监测到明显形变，形变量随时间推移逐渐增大。截至观测结束，涡阳县受煤矿开采影响区域，累积沉降量幅值达到83.4 mm，其余地区最大累计沉降量为27.3 mm；亳州市大部分区域地表累计沉降量处于5～30 mm 水平，平均累计沉降量为7.3 mm 左右。

图4 亳州市2021 年10 月至2022 年10 月时序累计形变量图Fig.4 Time-series accumulated deformation map in Bozhou City from October 2021 to October 2022

3.3 地理加权回归模型构建结果

为了能够有效地掌握建模数据的分布情况，本文使用最小值、中值、最大值及平均值对模型运算结果进行叙述性统计。各建模变量拟合系数如表5 所示，当系数为正时，自变量与因变量呈正相关关系；当系数为负时，自变量与因变量呈负相关关系，且拟合系数的绝对值越大，相关性越强。因此，各因素对地面沉降的贡献度排序依次为深层水位变幅、中深层水位变幅、中深层地下水埋深、深层地下水埋深、单位面积GDP、松散层厚度、道路密度、人口密度。

表5 模型运算结果叙述性统计Table 5 Descriptive statistics of model calculation results

4 讨论

4.1 采煤区地面沉降影响因素

亳州市煤矿资源主要分布在涡阳、蒙城两县，其中涡阳县现有矿产地13 处，是区内主要的采煤区。结合驱动力因子回归结果与煤炭实际开采情况可知，采煤区沉降受地下水抽取与煤矿开采共同影响，而煤矿开采是沉降严重区域形变的主导因素。驱动力因子回归系数显示，采煤区地面沉降与中深层水位变幅，见图5（c）、深层水位变幅，见图5（d）呈显著正相关，说明地下水水位变化对地面沉降具有一定贡献。而该地区沉降最严重区域位于涡阳县公吉寺镇以北的信湖煤矿，最大沉降达83.4 mm。信湖煤矿于2021 年9 月16 日正式投产，随着采矿活动的进行，矿区沉降速率持续加快。煤炭被采出后，形成采空区，随着采空区范围不断扩张，采空区上部覆岩和周围岩体的应力平衡遭到破坏，覆岩受到的重力作用逐渐增加，当压力超过临界值后，煤层顶板及周围岩体发生弯曲、断裂和垮落，导致整个上覆岩层的变形和移动，最终在地表形成大范围塌陷坑[30]。

图5 各因子对地面沉降影响的回归系数图Fig.5 Regression coefficients of different factors influencing ground subsidence

4.2 非采煤区地面沉降影响因素

非采煤沉降区主要位于谯城区东北部与利辛县西部，累计沉降量幅值为27.3 mm。驱动力因子回归系数显示，地下水状况与非采煤区地面沉降相关性较强，中深层水位变幅，见图5（c）和深层水位变幅，见图5（d）与地面沉降呈现正相关。自20 世纪80 年代以来，亳州市城市经济持续发展，人口迅速增长，对地下水资源的需求量也逐年增加。据有关资料估测[20]，亳州市目前地下水日开采量在3.5×105m3左右，深层地下水开采量在1.9×105m3。对深层地下水的过度开采，诱使承压水头持续降低，降落漏斗面积不断扩大。当水头压力差作用于下伏黏性土层时，黏性土层的中低压缩性，会导致其越流或者压密释水，引起自身测压水头下降，使土体被纵向压缩；而砂性含水层受水头减小的影响，会释放出一部分储存的水，使得含水层内部的应力状态发生变化。原本承压水头支撑的上覆载荷被转移至含水层砂砾间，致使砂砾间的有效压力增加，含水层被垂直压缩。在地表上监测到的沉降量，即为降落漏斗范围内黏性土层与含水砂性土层的压缩量之和[31]。

4.3 地质因素对地面沉降的影响

松散层厚度与地面沉降的关系如图5（e）所示，在亳州市中部和西南部回归系数为负，对地面沉降起抑制作用；在北部和南部回归系数为正，对地面沉降起促进作用。亳州市地处淮北平原，地下发育有第四系和新近系松散沉积物[20]，构成了地面沉降的物质基础。城市公共设施建设快速发展，城市建筑物的荷载不断增加，导致松散层被压实，进而引发地面沉降[22]；另一方面，松散层中不同土层持水性具有的很大差异，过度开采地下水，使得土层颗粒间的有效应力增大，孔隙体积被压缩，导致地面沉降。