基于SBAS-InSAR的滨州市地表沉降及其对影响因素分析

杜自豪，孙哲，刘智敏，武鹏达，戴昭鑫

(1.山东科技大学 测绘与空间信息学院，山东 青岛 266590； 2.中国测绘科学研究院，北京 100039)

1 引 言

地表沉降如何解决是在全球范围内的主要工程地质问题之一，它是由非人类或人类相关导致的地表以下松散层的巩固和压实[1～3]。近数十年来，由于城市工业化进程和大型城镇化工程的迅速实施，城市沉降现象越来越频繁，其会造成地表塌陷、城市市政基础设施被严重破坏，沿海城市还面临着海水倒灌的威胁，严重威胁到社会的可持续发展[4～6]。地表沉降严重程度和因素随着地域的差异而产生巨大变化，所以，通过探究地表沉降的危害机理就可以为地表沉降有效防止提出科学合理的决策根据。

目前，关于地表沉降时空演化规律和其形成机制的研究成果，主要集中在两方面：①对地表下沉及其影响因子的定性研究，如张金芝[7]等人研究发现，黄河三角洲地表沉降与油层开采和地表荷载有明显的相关性；Le[8]等人深入研究了岩溶坍塌、建筑荷载、土地开垦和利用，以及人类活动对地表下沉的影响；Tei je[9]等通过箱线图解析缅甸仰光地区的地表沉降发现，地下水抽取是导致该区域地表下沉的决定性原因；在剖析成都市地面变形和发生机制的过程中，薛廉[10]等人提出正在加快的城市城镇化建设和轨道交通快速发展，是造成城市地面沉降的主因。肖巍峰[11]在分析北京市地表沉降时，发现土地利用变化对地表沉降是有影响的；程霞[12]等研究黄河三角洲地表形变时发现盐化工厂开采地下卤水对地表沉降带来的影响是很大。②各因素之间的“横向”关系定量研究，如Merad[13]等人利用多因子评判模型，对矿井塌陷的影响开展了深入研究；Zhou[14]等人使用梯度提升决策树模型，定量分析了地下水位、可压缩沉积物厚度和地表负荷变化对沉降的影响；Omid[15]等人通过最大熵值法和规则集的遗传算法，研究了地表下沉特征与地质环境影响因子的相互关系。总的来说，地表沉降的原因是复杂的。尽管已有许多研究讨论过地面沉降的诱因，但目前仍面临着若干问题：对传统的研究中沉降影响各种因素(如地壳运动、地下水开采、地下水水位、断裂带等)都无法掌握，且资料收集时间较长，在适用程度上不能广泛的应用分析，作出合理有效的防治措施；关于影响地表沉降的因素，多是各种因素的“横向”对比分析，并没有对影响因素进行“纵向”分析，而这些研究偏重分析各影响因素对地表下沉的独立影响与意义，对影响因素相互影响则有相对较少的研究。

关于地表沉降监测，常规的监测手段有水准测量、全球导航卫星系统(GNSS)等科学技术虽可提升较大的准确度，但具有监测范围较小、工作量大和监测成本高以及站点较少等缺点。SBAS-InSAR技术采用的是多主影像的沉降观测，降低了主从影像发生时间和空间上的失相干现象的可能，且节约人力物力，大范围的沉降监测更加容易，同时结果展示也更加直观。为此，本文通过运用SBAS-InSAR技术获得滨州市地表沉降点数据，并利用地理探测器模型分析自然与人类活动影响因素对地表沉降空间分异性的影响，从而确定出地表沉降影响原因以及交互作用机理，为滨州市地面沉降科学治理提出了重要的理论依据。

2 研究区概况

滨州市(36°41′～38°16′N，117°15′～118°37′E)位于黄河三角洲西北部边缘，渤海湾西南岸，沿海湿地十分广阔，能源矿产和非金属矿产资源丰富，滨州地势由南向东北略有坡度。地表坡度下降1/6000～1/8000，整体海拔约为 -10 m～4 m[16]。2010年以来，滨州市市政府批复了《滨州北海经济开发区城市总体规划(2011-2030)》，使得土地扩张，城镇化建设加快，加快了地面沉降。选取了滨州市部分区县(惠民县、阳信县、滨城区、无棣县、沾化区)作为研究区，位置如图1所示。

图1 研究区地理位置示意图

3 数据与方法

3.1 研究数据

本文数据中主要包含地面沉降点数据和影响因子数据。地表沉降点数据是基于SBAS-InSAR方法利用2016-1-9—2020-9-20获取的覆盖滨州市部分地区的18幅上升Sentinel-TOPS合成孔径雷达图像(C波段)和辅助数据SRTM-1(shuttle radar topography mission)DEM数据及ESA提供的相应精密轨道数据进行地表沉降速率估算得到的。为了验证SBAS干涉测量结果，使用滨州市2016-2019年水准点数据进行校验。

自然因素和人类活动是影响地表沉降的两个因素，二者共同决定了地表沉降发生的情况。所以在选择沉降影响因子时必须兼顾以上这两个方面的因素外，还需要兼顾可获取性与适用性，这里选定3类共16种影响因子，包括：土地利用类型。这里采用的是第三次国土调查数据，按照 1 km×1 km分别统计格网内一级类土地利用类型。土壤质地，其根据 1∶100万土壤类型图和第二次土壤普查数据编制而成，包括砂土、粉砂土、黏土。POI数据从高德开放平台抓取获得，这里筛选出盐化工产业和化学化工类企业，对其进行核密度分析和重分类处理。

3.2 研究方法

3.2.1 SBAS-InSAR技术原理

SBAS-InSAR技术旨在提高变形观测的质量，通过使用多景SAR数据来准确观测长时间的小变形现象，并使用新形成的具有相对较小时间和空间基线特征的干涉图序列来最好地补偿过去相关的影响。

假设测区N+1景雷达影像所在(t0，…tN)的时间段里，根据干涉关系的自由组合，可以得到M幅干涉图，并满足：

(1)

首先，从tA，tB两景的图像中去除不同时间的地形相位，根据假设tA

δφi(x，y)=φ(tB，x，y)-φ(tA，x，y)

(2)

其中φ(tA，x，y)和φ(tB，x，y)分别表示与tA，tB时间相对应的相位，λ表示雷达波长，d(tA，x，y)和d(tB，x，y)分别代表在tA和tB各自时间每个像元(x，y)在LOS(line of sight)方向上的地表形变，相对于参考时刻t0。

方程(2)中的相位φ表示为平均速度v和任何两个测量时间t之间的时间的乘积，被定为具有物理含义的沉降序列，表达式为：

(3)

故，第i幅干涉图的相位可表示为：

(4)

将其改写为矩阵形式：

Bv=δφ

(5)

式中：B—M×N矩阵。最后，用SVD方法可以得到速率矢量v的广义逆解，然后通过对每个时间段的速率在一个时间区间内进行微分，可以找到每个时间段的形变量。

3.2.2 精度验证

为检验使用SBAS-InSAR技术分析成果的真实性，根据《CHT6006-2018时间序列InSAR地表形变监测数据处理规范》[17]，在这里使用二等水准实测数值对InSAR地表变形监测结果进行精度分析。为确保精度分析的客观与精确，采用邻近点原则，即对水准点数据设置 80 m缓冲区，并比较水准点数据与InSAR监测结果的差值，精度评价指标为：

①误差绝对值，宜优于15 mm/a；

②误差标准差，宜≤8 mm/a。

3.2.3 地理探测器

地理探测器是王劲峰[18]发明的一种统计学方法，用于检测空间结构异质性现象并揭示其背后的驱动力，它被广泛用于分析自然环境和人类经济机制的空间变化的驱动力。地理探测器模型由四个模块组成：因子探测、交互探测、风险探测和生态探测。地理探测器不仅可以分析每个影响因素之间的相互作用，还可以分析两个影响因素中每个因素之间的相互作用。本文将分别使用单因素探测和交互探测模块，解析单因素与不同影响因素之间的交互作用及其对地面沉降的影响。

(1)因子探测器

因子探测器用来分析自变量(不同驱动因素X)对因变量(Y)空间分布形成的影响程度，识别不同驱动因素对因变量空间的驱动力大小(q值度量)，表达式为：

(6)

(2)交互探测

交互探测主要是为了确定各种影响因子间的交互多对地面沉降格局的影响关系，即各种影响因子在联合相互作用时，是否会增加或降低对因变量Y空间分布的影响程度。通过分别计算值q(X1)、q(X2)值及两因子图层叠加后的q(X1∩X2)数值，并对比前三者数值大小，从而确定了两因子之间是否具有相互作用，以及对地面沉降空间格局的影响力是加强还是削弱，并由此推广到三因子的相互影响。

因子之间交互作用分别表现为增强、双因子增强、非线性增强、减弱、非线性减弱、独立，判断方式如表1所示：

交互作用类型 表1

4 结果与讨论

4.1 基于SABS-InSAR的地表沉降监测结果及精度验证

通过SARscape软件对2016-1-9—2020-9-2018景Sentinel-1A数据进行SBAS-InSAR处理，获得地表沉降速率矢量点，经过克里金插值得到研究区的沉降速率图，如图2所示，正值代表地表抬升，负值代表地表下降。其中，研究时段内地表沉降平均沉降速率范围为 -62.8 mm/a～34.6 mm/a，沉降较严重区域集中分布在滨州市无棣县和沾化区的北部沿海地区。

图2 基于SBAS-InSAR技术监测的2016-1—2020-9研究区沉降速率分布结果

为了验证SBAS-InSAR技术对地表监测的可靠性，现用21个二等水准点数据建立80米缓冲区与缓冲区内的沉降数据进行比较，得到结果如表2，可以得到缓冲区内的沉降点有259个，与水准点进行比较，可以看出，误差与标准差都符合《规范》[17]，其中差值绝对值的 ≤6 mm/a占134个，差值绝对值分布如图3所示，证明InSAR监测结果具有较高的可靠性，可用于研究和分析。

水准数据与InSAR监测结果比较 表2

图3 差值绝对值分布

按照《地质灾害危险性评估规范》对地表沉降速率划分等级，探究地表沉降空间分布情况(图2)和沉降面积统计(图4)。由图4可知：研究区以较低沉降为主，占有 3 351.17 km2，其主要分布在盐田区域的埕镇和马山子镇北部，以及滨城镇；沉降区域最严重仅有 4.66 km2，分布在滨州市北部与渤海相接区域。研究区总体以较低沉降为主，呈北部沉降，南部沉降较小分布状态。

图4 研究区沉降面积统计

4.2 基于地理探测器的地表沉降影响因素分析

4.2.1 地表沉降影响因子单因子探测

图2显示，研究区域内的沉降速率存在明显空间异质性现象，因此，可以利用地理探测器对其进行定量分析，研究地面沉降-空间分布的影响机理。在本研究中，以2016-1-9—2020-9-20地表沉降速率为因变量Y，以土地利用类型分类面积为自变量X，进行地理探测器探测，计算其类内影响因子对地表沉降速率分布的驱动力值(q值)，其强度大于0.01的驱动力如表3所示。

土地利用类型对地表沉降速率的影响力 表3

从表3中可以看出，工矿用地对地表沉降速率的作用强度最大，工矿用地包含盐田、工业用地、采矿用地。盐田区的晒盐及地下卤水的开采，导致地下水位下降，从而引起地表沉降。

针对土地利用类型对地表沉降速率的影响，得到工矿用地是影响地表沉降的主要因素，在研究区内，工矿用地占比最大的是盐田，如图5(c)所示，因此，这里对研究区内的POI类型中的盐化工产业及化学化工类企业进行核密度分析，对其按照自然间断法重分类，作为自变量X输入地理探测器，得到其驱动力如表4所示。从结果表中可以得到盐化工产业的q值最大，盐化工产业分布如图5(a)所示，从图中可以看到，沉降区内的盐化工产业分布相对较多，而盐化厂是生产盐和化工原料的基地，如卤水开采，溴素提取。如果是长期或长期不合理利用地下卤水将会导致地下卤水水位的降低，从而产生地面沉降。

图5 实验数据

POI类型对地表沉降速率的驱动力 表4

土壤质地对地表沉降速率的驱动力 表5

从图5(a)盐化工产业POI的分布可以看出，除了沉降区分布有盐化工厂外，在滨州的东部和西南部均有分布，但对比发现，这些地方未发生严重沉降，初步猜想，可能是工矿用地区域与其他因子结合，对沉降区域的影响有所增强，因此，引入土壤质地类型数据，对地表沉降进行地理探测器分析，进一步探究沉降分布机制，得到土壤质地对地表沉降的驱动力如表5所示。由此可以得出，砂土的驱动力远大于粉砂土与黏土对地表沉降的影响，砂土液化导致地表沉降。

4.2.2 地表沉降影响因子交互探测

土地利用类型、设施类型POI数据、土壤质地中对地表沉降影响的最强的是工矿用地(q=0.294，p=0.000)、盐化工产业(q=0.063，p=0.000)、砂土(q=0.373，p=0.000)，对研究区进行双因子交互和三因子交互，其作用探测显示如表6、表7，在表中可发现，不同因素相互影响对地表沉降的影响力均超过单一作用的影响力，并体现为双因子相互增强，其中，砂土分别和盐化工产业、工矿用地进行交互作用，对地表沉降的影响力均达到了 0.415 9、0.402 9，三因子交互作用影响力达到了0.624。

双因子交互对地表沉降的影响 表6

三因子交互对地表沉降的影响 表7

根据表6，表7的结果，其可能解释如下：

(1)滨州市位于断裂带之上[19]又位于渤海湾南岸，历史上曾出现多次地震[20]，使得滨州无棣县和沾化区北部沿海地区沙土地质受到外力后，使得大面积区域的砂土发生液化，其位置区域如图6所示。当沙土变成液体时，孔隙水在超孔压力下上下移动，当沙土上不再有不透水的覆盖物时，地下水就会自动排出土壤；如果沙土层上有透水性低的黏土，超孔水的压力超过了覆盖物的强度，地下水就会通过覆盖物输送沙粒，沙土的保水能力很差。因此导致地表不均匀沉降，再加上盐化工产业主要是以开采地下卤水为主要原料制取溴化合物和盐，而过量利用地下卤水都会导致地下卤水水位下降，两者相互促进，从而导致地面沉降状况越来越严峻。

图6 滨州市无棣县、沾化区砂土液化程度分区图[16]

(2)从图5(c)可以看出，无棣、沾化和北海新区工矿区域分布最多的是盐田，宜盐面积可达9.6万公顷，现已有盐田100万亩，年产均200万吨，根据2016-2019年滨州《统计年鉴》可得原盐产量如表8所示，“井滩晒盐”的主要采用的是地下卤水，已探明静态储

量30亿立方米。因此，地下卤水的长时间开采，再加上砂土质地的液化，导致地下水位下降，加剧地表沉降。

2016-2019滨州原盐产量(万吨) 表8

(3)土地利用类型因子与设施类型POI因子中影响最大的工矿用地和盐化工产业，虽然这种交互作用后的影响力是双因子增强，但并不像砂土分别与盐化工产业和工矿用地交互后，对地表沉降作用强度大，表明这两个因素在影响地表沉降过程中可能有重叠，盐田制盐和盐化工产业的生产都需要对地下卤水进行开采，导致地表沉降。

5 结 语

以滨州市5个区县为主要研究区，利用SBAS-InSAR技术获取了2016-1-9—2020-9-20的18个时间段地表沉降数据，并与土地利用类型数据、土壤质地数据和设施类型POI等多源数据相结合，基于地理探测器模型，定量监测了影响地表沉降的各种因素，并深入分析了影响因素之间交互作用对地表沉降的影响力。该研究表明：

(1)研究区地表沉降以中等沉降为主，沉降区域主要分布于滨州市无棣县和沾化区北部，南部沉降区域较少。

(2)单因子探测结果。对土地利用数据、设施类型POI数据、土壤质地数据中对地面沉降影响力较大的因素依次是工矿用地、盐化工产业和砂土土质。

(3)因子交互探测结果。单因子影响力最强的三个因子双交互作用后都表现为双因子增强，其中砂土分别与工矿用地、盐化工产业相互叠加后，综合影响强度体现得更为突出；三因子作用后对地表沉降的影响则尤为突出。

由于滨州市北部盐田产量的逐渐增加和《滨州北海经济开发区城市总体规划(2011-2030)》，随着地下盐卤过量开发以及城市化飞跃型发展趋势，目前滨州市及其以北区域已发生了强烈的地面沉降现象，建议在以后的开发建设过程中，合理利用资源，减缓地表沉降。本研究在考虑研究区位置和人文情况下，仅选用了土地利用和土壤质地因子，并引入大数据方面的POI因子，并将空间统计学方法应用于地表沉降，定量分析影响因子对地表沉降的作用强度，今后可收集地下卤水或地下水位数据作为补充数据研究地表沉降，并探测其与他影响因子交互对地表沉降的作用。