基于PS-InSAR数据的北京市地面沉降影响因素及其交互作用探测

0　引言

地面沉降是在自然与人类活动交互影响下，由于地下松散地层固结压缩导致区域性地面高程降低的一种不可逆的地质灾害现象，具有发展缓慢、不易引起人们察觉、波及范围广、持续时间长等特点。近几十年来，随着工业化进程和大规模城市化建设的快速推进，沉降现象越来越普遍，严重的沉降会导致地面塌陷、建筑物开裂、市政设施被破坏等一系列灾害问题，影响人类社会的可持续发展。截至2011年底，中国50多个城市出现地面沉降，其中，华北平原、长三角地区和汾渭盆地已发展成了沉降重灾区。地面沉降通常呈现区域性特征，即沉降量及成因随区域的不同而发生改变，因此，探测区域地面沉降的影响机制能够为地面沉降有效防治防控提供科学的决策依据。

目前，关于地面沉降影响机制的研究主要集中于以下两方面：①地面沉降与影响因素的定性研究，如Teije等人利用箱线图等分析缅甸仰光地区的地面沉降，发现地下水抽取是引起该地区地面沉降的决定性因素；薛廉等人在分析成都市地面变形及发生机制的过程中，认为快速的城市化进程和轨道交通建设的迅猛发展是导致市中心区域地面沉降的主因；荆创利等人揭示了地壳运动、地下水开采、地表建设等自然因素和人为因素所导致的地面沉降规律。②采用基于机器学习或空间统计方法分析各影响因素对地面沉降的作用，如Zhou等人基于梯度提升决策树模型分别定量分析了地下水位变化、可压缩沉积层厚度、地面荷载对沉降的影响程度；Omid等人采用最大熵值法和规则集的遗传算法研究了地面沉降特征与地质环境影响因素的关系。虽然学者们对地面沉降的影响机制进行了一系列研究并取得了较为丰硕的成果，但大多数研究侧重于分析各影响因素对地面沉降的独立影响作用，对因素间交互作用鲜有探讨。而实际上，多种因子之间的复杂相互作用共同决定了地面沉降的空间格局。探测地面沉降分布的主要影响因子以及影响因子间交互作用有助于更全面地掌握区域地面沉降分布机制。

传统的地面沉降监测手段（如地面监测站、高精度水准测量网、全球导航卫星系统（ G l o b a l Navigation Satellite System，GNSS）测量网）能提供高精度的变形信息，但具有监测成本高、监测范围小等缺点。永久散射体干涉测量（Persistent Scatterer Interferometric Synthetic Aperture Radar，PSInSAR）技术利用稳定且小于像元尺寸的永久散射体，能有效地去除大气效应贡献值，获得高精度的地表形变值，具有监测范围广、成本低、精度高等特点，被广泛地用于大范围的地表变形监测。因此，本文利用PS-InSAR技术获得北京市地面沉降点（Persistent Scatterer，PS）数据，借助于地理探测器模型分析自然与人为要素对地面沉降分异的影响，识别出地面沉降影响因素及其交互作用机制，为理解北京市地面沉降的影响机制以及科学防治提供重要的依据。

1　研究区概况及方法

1.1 研究区与数据来源

北京，是我国的首都，也是世界上著名的国际化大都市，位于115°25′～117°30′E，39°26′～41°30′N，是国内地面沉降严重的城市之一。选取北京市主城区（海淀区、朝阳区、东城区、西城区、石景山区、丰台区、顺义区、通州区）为研究区，如图1所示，该区域位于北京市东南部，是人口较为集中、社会经济发展水平较高的区域，也是沉降严重区域。据《2018年北京统计年鉴》记载，截至2018年底，该区面积是3 303.25 km²，常住人口达1 440.6万，2018年实现地区生产总值20 267亿元。研究区地形以平原为主，西部与东北部有少量山地丘陵，平原面积约占总面积的98%，属于暖温带半湿润大陆性季风气候，夏季高温多雨，冬季寒冷干燥，降水季节分布极不均匀，全年约80%的降水集中在夏季6—8月。

本文数据主要包括地面沉降PS点集数据和影响因子数据。PS点集数据是利用GAMMA软件处理24景空间分辨率为30 m的SAR影像，获取了2003年6月至2010年9月的23个时间段总计230 028个PS点数据，包括经纬度坐标、时间、年均沉降速率、累积沉降量等属性信息，获取时间间隔不定，由2～10个月不等，PS点集密度如图1所示。

图1　研究区及PS点密度

Fig.1 Study area and PS point density

地面沉降发生受自然因素和人类活动影响，两者的作用方向和强度共同决定了地面沉降发生的状况，因此在选取沉降影响因子时需要综合考虑这两方面的因素，考虑到研究区的实际情况和数据的可获取性，选定下列7类共14种影响因子，包括：①粘性土可压缩层厚度。可压缩层粘性土受地面荷载的等外力作用容易压缩，进而引起地面沉降，研究区内可压缩层厚度取值范围是51.2～236.8 m，其中大于200 m部分占很小比例，部分区域数据缺失，为此，将其分为5个等级（50～100 m、100～150 m、150～200 m、大于200 m以及缺失类）。②土壤类型。通州区、顺义区主要以潮土为主，西城区、东城区、朝阳区、海淀区、丰台区等建成区以城市混凝土为主。③地形因子。包括高程、坡度和坡向。④地下水位高度。以统一高程作为基准面，地下水水位高度数值越小，表示地下水水位越低，地下水含量越少，地下水被超采。考虑到研究区地下水水位高度范围，以5 m和10 m为等高距，将其分为12个等级（小于-10 m、-10～-5 m、-5～0 m、0～5 m、5～10 m、10～15 m、15～20 m、20～25 m、25～30 m、30～40 m、40～50 m、50～60 m）。⑤地质类影响因子。包括距断裂带距离、距5级以下地震距离(简称小地震源)和距5级以上地震距离(简称大地震源)。⑥土地利用类型。分别为草地、林地、地表水、农用地、建设用地和未利用地。⑦交通可达类影响因子。包括距国道距离、距省道距离、距铁路距离、距河流距离。影响因子数据来源于地理空间数据云平台、国家地球系统科学数据中心等，详细数据来源见表1。

表1　影响因子数据来源

Tab.1 Data source of associated factors

利用ArcGIS软件将点、线类型的因子数据转换成到格网的距离值，在此基础上，通过Python语言对各因子属性数据表格进行处理，基于面积占优原理对各影响因子数据实施格网化，并与PS点网格一一对应，实现属性关联。

1.2 研究方法

1.2.1 基于空间自相关的PS点格网化

由于地面沉降数据是离散PS点，空间上呈离散分布，为了便于后续的空间关联分析，需对离散PS点进行不同尺寸的网格化。空间自相关性是空间数据所特有的性质之一，网格化后的数据应以最大限度保留原始数据的空间自相关性为前提，常用的测度方法有莫兰指数（Moran's I），故本文采用空间自相关的Moran's I系数差法来确定最优空间格网，其计算公式如下：

式中，Ig和Ib分别为格网化后数据和原始沉降均值数据的Moran's I系数；S为格网化数据与原始沉降均值数据的Moran's I系数差值；对莫兰指数差值取绝对值，S值越小，表明该尺寸下格网化后数据与原始数据的自相关系数越相近，网格采样后的数据较好地保留了原始数据的自相关结构，则可将该尺寸作为最优的空间格网化参数。

1.2.2 地理探测器

地理探测器是探测空间分异性并揭示其背后驱动力的一组统计学方法，广泛应用于地理要素空间变异的驱动力分析、精准扶贫、共享单车骑行时空分布特征等诸多领域。地理探测器模型包括因子探测、交互探测、风险探测和生态探测4个模块，这些模块既可以分析单影响因子的作用，也可以分析两影响因子交互的作用。本文分别采用因子探测与交互探测两个模块解释单因子和不同影响因子之间的交互作用对地面沉降影响力的作用。

1）因子探测。其核心思想是比较某一地理环境因子（自变量X）与地理事物（因变量Y）的变化在空间上是否具有显著的一致性，如地理环境因子与地理事物的变化具有一致性，则说明该地理环境因子对地理事物的发生和发展具有决定意义，用q值来度量，其计算公式为：

式中，h=1,2,…,L为地面沉降发生风险各指标的分层；Nh和N分别为h层和全区的单元数；σh²和σ²分别为h层和全区因变量Y值的方差；SSW和SST分别为层内方差之和和全区方差；q值的范围为[0,1]，q值越大说明因子X对地面沉降分布的空间分布影响力越强。进一步，通过非中心的F分布的显著性检验结果（用p值表示）对q统计量进行显著性检验，若p值小于预设的显著性水平（0.05），则拒绝原假设（分层条件下Y变量呈现随机分布），而接受备选假设（分层条件下Y变量呈现空间聚集分布），否则，q值不显著，没有充分的证据拒绝原假设。

2）交互探测。主要用于研究两个不同影响因子之间的交互对地面沉降格局的作用关系，如对于影响地面沉降格局影响因子X1和X2，通过空间叠加X1和X2形成新图层X3，X3的属性由X1和X2共同决定。通过比较X1和X2层的影响力和X3的影响力，可以判断两个因子交互作用对地面沉降空间格局与单个影响因子对地面沉降空间格局的影响力是强化还是弱化。

交互探测由下列表达式组成：若q (X1∩X2 )＜min(q (X1 ) ,q (X2 ))，则表示因子X1和X2交互后非线性减弱；若min(q (X1 ) ,q (X2 ))＜q (X1∩X2 )＜max(q(X1),q(X2))，则表示因子X1和X2交互后单因子非线性减弱；若q(X1∩X2)＞max(q(X1),q(X2))且q(X1∩X2) ＜q(X1)+q(X2)，则表示因子X1和X2交互后双因子增强；若q(X1∩X2)＞q(X1)+q(X2)，则表示因子X1和X2交互后非线性增强；若q(X1∩X2)=q(X1)+q(X2)，则说明X1和X2相互独立。本文计算结果均使用软件GeoDetector计算获得。

2　结果分析

2.1 最优网格尺寸

采用ArcGIS软件分别计算PS原始点集与不同格网尺寸下的Moran's I系数，并求得其系数差。由图2可知，当格网尺寸为1 700 m时，Moran's I系数差值最小为0.062，因此，本文基于1 700 m×1 700 m网格对原始沉降PS点实施重采样，得到地面沉降栅格数据，作为后续关联分析和地理探测的数据源。

图2　格网数据Moran's I系数差

Fig.2 Moran's I difffference of grid data

2.2 影响因子离散数据分级

以时序沉降数据获得的多年年平均沉降速率定义为地理探测器中因变量Y，其分布格局如图3所示，上文提到的14种影响因子作为地理探测器中的自变量X1， X2，…，X14，其中，因变量多年年平均速率和自变量距断裂带距离、距河流距离等10个影响因子均为连续型变量。由于地理探测器中的影响因子变量应为离散型，因此，首先需对连续型的自变量与因变量进行基于自然断点法的离散化。自然断点法是一种通过减少类内距离而增大类间距离的自然聚类方法，即通过迭代比较每个分组和分组中元素的均值和观测值之间的平方差之和来确定值分组中的最佳排列，分类后数据具有强度大、表达效果良好的特点，被广泛应用于离散化数值型数据。对因变量多年年平均速率和自变量距断裂带距离、距河流距离等10个连续型影响因子进行分级，当类间距离与类内距离之差最大，则表明该分类级别是最优的。其中，多年平均沉降速率的分级如图4所示，总的分类结果见表2（表中：Y为多年年平均速率，X1为距断裂带距离，X2为距河流距离，X3为距国道距离，X4为距省道距离，X5为距铁路距离，X6为高程，X7为坡度，X8为坡向，X9为距5级以上地震源距离，X10为距5级以下地震源距离）。如离散指标Y，当分级数量为11时，类间距离与类内距离之差为24.81达到最大，故对指标Y分类的最优分类数目为11。类似地，计算其他指标因子的类间距离和类内距离，以两者差值最大为原则来确定其相应的最优分类数目。

图3　多年平均沉降速率分布格局

Fig.3 Distribution pattern of annual average settlement rate

图4　多年平均沉降速率分级

Fig.4 Classifification of annual average settlement rate

表2　离散数据自然断点法分级

Tab.2 Classifification of discrete data based on natural break method

2.3 地面沉降影响单因子探测

从图3可以发现，研究区内的沉降速率具有明显的空间分异性，因此有必要采用地理探测器模型进一步定量分析地面沉降空间分异的影响机制。通过地理探测器计算各影响因子对地面沉降速率分布的影响力（q值），结果表明只有8种地理环境要素影响了北京市地面沉降速率分布的空间格局，其作用强度由大到小分别为：地下水位(q=0.41,p=0.00)、可压缩层厚度(q=0.33,p=0.00)、坡度(q=0.08,p=0.00)、土地利用类型(q=0.07,p=0.00)、土壤类型(q=0.04,p=0.00)、高程(q =0.03,p =0.00)、距5级以上地震源距离(q=0.02,p=0.01)和距铁路距离(q=0.01,p=0.01)，其他6种因子没明显的影响力。

探测结果表明人类活动引起的地下水位下降对研究区的地面沉降作用强度远远大于其他因子。研究区北京是一个严重缺水的特大型城市，城区2/3以上的供水来源于地下水，地下水开采引起的地下水位下降是北京地面沉降形成与发展的关键外在因素。在满足显著性水平0.05的前提下，可压缩层厚度、坡度、土地利用类型、土壤类型、高程、距5级以上地震源距离这6种因子的作用强度存在较大差异，其中，可压缩层厚度因子的影响力远远大于其他因子，与地下含水层所对应的可压缩层组的分布、结构、厚度是形成地面沉降的主要内在因素。此外，距铁路距离对整个研究区沉降的作用强度不大，但对铁路沿线地区的沉降还是有一定作用强度。

2.4 地面沉降影响因子交互探测

交互作用探测结果显示见表3，大部分不同因素交互作用对地面沉降的影响力均大于单独作用的影响力，各影响因素间的交互作用类型有非线性增强型和双因子增强型两种（表中A↑B表示A和B非线性增强；A↑↑B表 示A和B双因子增强）。其中，地下水分别与距大地震源距离、距小地震源距离、土壤类型、距省道距离、距河流距离、距国道距离、距断裂带距离、距铁路距离、坡向在影响地面沉降方面是非线性增强的；而地下水分别与可压缩层厚度、土地利用、坡度具有交互协同的作用，即两因子在影响地面沉降方面是双因子增强的，效果没有非线性增强型显著。可压缩层厚度分别与土壤类型、距大地震源距离、距省道距离、距小地震源距离在影响地面沉降方面是非线性增强的。

表3　交互因子对地面沉降的影响(q>0.4)

Tab.3 Inflfluence of interaction factors on land subsidence (q>0.4) 

各类因素交互探测的结果显示，地下水位和可压缩层厚度作为对地面沉降最大影响的单因子，在与其它较小影响力的因子叠加时，对地面沉降的综合作用是增强的。根据表3的结果，其可能解释如下：

1）单因子影响力最大的地下水位和可压缩层厚度指标，交互作用后的影响力并不是最大，说明两种因子相互叠加后在解释地面沉降过程中可能存在作用重复，地下水水位与可压缩层厚度相互影响，地下水大量抽取导致地下水水位下降，可压缩层内缺少填充物而发生压缩，外来补充水源无法充分吸收补充地下水，两者相互促进，地面沉降现状不断恶化。

2）土壤类型与地下水位存在密切关系，降水是地下水补给的主要来源，降水形成的土壤水分入渗将穿过非饱和地带形成对地下水的补给。由于不同土壤类型的土壤颗粒间缝隙不同，其渗水性能不一致，同时土壤蒸发锋面也不同，导致降水形成地下水补给的时间也不同。降水对地下水补给与地下水位埋深具有密切联系，如地下水位埋深较浅时，降雨形成的土壤水分入渗将在较短时间内穿过非饱和地带，形成对地下水的补给，雨停后发生的蒸发，其土壤蒸发锋面将缓慢下移最终抵达地下水面，形成潜水蒸发。此外，土壤类型作为地面沉降的下垫面，在影响地下水分布的同时，导致土壤空隙减小，有可能在一定程度上影响可压缩层状况。

3）土地利用类型变化主要是由于人类活动导致的，如成片的高大建筑物和铁路工程建设，大大增加了地面荷载，导致地面高程损失而直接形成地面沉降。同时，由于工程建设导致的土地利用类型变化而形成的不透水地面，严重阻碍了降水的地面渗水率，不能快速有效地形成地下水补给，间接影响了地下水位的提升，加剧了地面沉降的发生。不同类型的地面植被，需要不同的灌溉条件，也会对地下水位的变化产生影响。

4）在北京区域地下水超采不严重的地区引起的水位趋势性变化中，包含了地震前兆信息，也就意味着地下水位降低有可能会诱发地震的发生。同时，离震源越近的地区地下水位变化的可能性越大，可能的情况是地震的发生与地下水位变化具有密切联系，而地下水位变化又与地面沉降的发生密切相关。

3　结束语

以北京市城中心的8个区为例，基于PS-InSAR技术获取了时间跨度为2003年6月至2010年9月的23个时间段PS点地面沉降数据，以及包括地下水、可压缩层厚度等多源数据，建立了以格网为单元的地理探测器模型，精细地探测了地面沉降分布的影响因素，并深入分析了各影响因素间交互作用产生的影响力。研究发现：

1）单因子探测结果。对地面沉降影响力较大的因子依次是地下水位、可压缩层厚度、坡度、土地利用类型、土壤类型、高程、距大地震源距离、距铁路距离，其中，除了地下水位和可压缩层厚度两个因子对地面沉降产生显著的作用外，其余要素产生的作用强度显著减弱。

2）因子交互探测结果。作用强度较弱的单因子与地下水位或可压缩层厚度分别相互叠加后，其综合作用强度均表现为增强作用，大部分表现为非线性增强，极少部分表现为双因子增强。

随着地下水的过度开采与城市化进程的迅猛发展，中国的许多城市和地区出现了严重的地面沉降现象，带来了一系列的地质环境问题。本文主要聚焦地面沉降分布空间分布特征与其影响因素差异，其研究结果可为城市地面沉降防控提高科学的参考依据。本研究是将空间统计学方法应用于PS-InSAR大数据的一次有益的尝试，但目前还处于探索阶段，在后续研究中将从地面沉降时空变化特征来探测其影响因子的交互作用变化特征。

致谢：感谢中国测绘科学研究院张永红研究员提供的研究区PS-InSAR数据以及深圳大学黄金彩博士对研究区PS-InSAR数据的整理与预处理。

作者简介：

肖巍峰（1977-），男，湖南湘潭人，讲师，博士，主要从事三维点云数据处理及表面重构、地理空间数据分析等研究工作。

E-mail:472364382@qq.com

通讯作者：杨文涛（1988-），男，湖北孝感人，讲师，博士，硕士生导师，主要从事地理空间数据分析与挖掘理论、方法与技术等研究工作。

E-mail:yangwentao8868@126.com