基于InSAR技术的深圳市地表沉降空间格局分析

文｜深圳市城市公共安全技术研究院有限公司、深圳城安软通科技集团有限公司 金典琦/深圳市城市公共安全技术研究院有限公司 陈天东/深圳城安软通科技集团有限公司 况达/通讯作者：深圳市城市公共安全技术研究院有限公司 张会

1 引言

地表沉降是在自然与人为因素的作用下，由于地壳表层土体压缩等原因造成的地表标高变化的现象，对人类生活特别是人力城市发展的影响较大。地表沉降在城市中造成的危害，由于使用水准仪进行的水准测量、利用GNSS终端的测绘设备进行的测量给出的变化都是点状形变量，观测点的密度、成本和周期都难以覆盖全域。近年来，由于具备高精度沉降监测能力，合成孔径雷达干涉（InSAR）技术在地表形变监测方面得到应用和发展。与其他测量技术的对比，InSAR技术具有独特优势特点：（1）全天时全天候，受天气影响小；（2）测量范围广；（3）历史沉降信息回溯；（4）实施方便，建设及运营成本相对较低。为城市区域沉降变形的及时定量监测提供了可能。形成了差分干涉测量（DInSAR）、小基线干涉测量（SBAS-InSAR）和永久散射体干涉测量（PS-InSAR）等地表沉降监测方法。

沉降灾害的诱发因素包括自然因素（气候变化、地下水位随季节的涨跌以及地质运动等）和人为因素（采矿、高层建筑、深基坑施工、地下水枯竭和植被破坏等）。自然和人类因素本身具有的空间异质性特征，决定了区域地表沉降也应表现出空间异质性，本文利用永久散射体干涉测量（PSInSAR）技术获取城市大范围的地表沉降量，以街道-社区为基本空间单位，利用全局莫兰指数（Global Moran's I）和聚类和异常值分析 （Anselin Local Moran's I）以及热点分析（Getis-Ord Gi*）分析地表沉降现象的空间格局特征，揭示深圳市地表沉降空间格局和空间分异规律，为研究地表沉降诱发的各类灾害触发机理、为沉降风险预警、灾害预防和治理提供支持。

2 研究区概况与数据

2.1 研究区概况

深圳市位于广东省南部,陆域位置是东经113°46′～114°37′,北纬22°27′～22°52′。深圳地貌主要包括三类地貌带，即山地地貌带、台地地貌带、海岸地貌带。其中，山地地貌带主要是低山、高丘陵、低丘陵三种类型；台地地貌带主要是高台地、中台地、低台地；海岸地貌带主要是海河混合堆积平原、生物成堆积滩地、阶地滩地、海成堆积、河成堆积阶地平原等。西部地区山地地貌带基本以羊台山为中心,南至塘朗山、西至凤凰岩、北至吊神山、东到打石坑鸡公头,由高丘陵和低丘陵构成一个近于东西轴向的椭圆环状分布区；环绕着这个内环分布区外圈大沙河下游、铁岗水库北侧、石岩、光明横岗水库、平湖、布吉河、梅林等,由高台地、中台地和低地构成了不连续台地地貌带外环分布区。东部地区大致以葵涌盆地北北东向为界划分为东西两区：以梧桐山、梅沙尖、马峦山、田心山为代表的西区山地地貌带，主要由低山、高丘陵组成，呈北东方向展布，延续至北部依次为高丘陵或低丘陵中间分别被坪山河、坑梓河、龙岗河等谷地间隔，这些谷地的两侧往往形成台地地貌带。在龙岗、坪山、坪地、坑梓、横岗、西坑、葵涌等谷地平原之下，常形成隐伏的岩溶地貌；东区的大鹏半岛被夹于大亚湾与大鹏湾之间，从南面的七娘山、排牙山、笔架山等由低山、高丘陵、低丘陵组成的山地地貌带，呈北西走向。深圳海岸线总长257km，东海岸线长13.22km，以基岩一砂砾质海岸类型为主；西海岸线长123.78km属淤泥质海岸，部分属三角洲河口岸及红树林海岸。

2.2 数据

研究的地面沉降监测，共获取了2016年月-2020年8月时间范围内覆盖研究区的期高分辨率、长时间序列COSMO-SkyMed重复轨道SAR影像，数据参数及监测时间如表1所示。

表1 监测区域SAR数据的基本参数

图1 研究区概况

3 数据处理方法及结果

3.1 PS-InSAR沉降测量技术

PS-InSAR的处理过程主要包括以下步骤，对于K+1幅SAR单视复数影像，主要经过图像的几何校正、定标、永久散射体探测、干涉和差分干涉等步骤，利用监测区域的激光点云数据生成的数字高程模型进行差分干涉，得到K幅干涉、差分干涉图和差分干涉相位集。每个PS点在差分干涉相位都受到了监测区的地表沉降、高程误差、大气影响及失相关的综合影响，PS点的差分干涉相位组成如下：

根据式（2）差分干涉相位模型可采用数字高程模型的高程误差和地表形变带来的线性沉降速率函数模型进行表达，即：

模型式（3）较难直接求解，因此，采用PS点的相位各分量的空间相关性进行间接求解，即通过建立PS邻域差分相位模型进行求解，

则：

如有K幅差分干涉图，对于每对PS点，可计算出K个时序邻域差分相位，并建立K个（4）式组成的方程组：

通过对此积分，可以获得相对于主参考点的解缠线性相位残差，即：

对时序残留相位取均值，并作空间平滑即可消除大气不均质分布对于相位的贡献量，即主影像大气相位：

利用上述计算得到PS主影像和各干涉图中的大气相位后，主影像的大气相位与各干涉图大气相位之差，即可得到各干涉图相对应的从影像的大气相位：

其在干涉图中的非线性沉降相位为：

3.2 数据处理结果

采用PS-InSAR算法计算沉降数据共5104827个PS点，计算结果及数据分布情况如图3所示。

图3 地面沉降空间分布图

4 地表沉降的空间格局

图2 PS-InSAR方法基本流程图

深圳的地表沉降多由人类活动扰动引起，因而其地表沉降具有较强的空间关联特征，采用Moran'sI指数定量的分析其全局空间格局，采用local Moran'sI定量分析人类活动带来的地表沉降现象的空间异质性带来的局部空间格局特征，根据Anselin Local Moran'I算法，分析沉降现象的空间格局的聚类/异常值特征类型。将全市的675个社区作为基本分析单元，利用各分析单元内PS点的地表累积沉降平均值、中位数，进行空间权重矩阵、Moran'sI指数、local Moran'sI指数计算。

4.1 全局空间格局

Moran'sI指数的取值范围在-1～1之间，-1表明存在强烈的空间负相关关系，0表示属于完全随机模式，1表示存在显著的空间正相关。Moran'sI指数的计算如公式（11）所示。

式中：r2为地表沉降速率的方差值；xi和xj则分别为i社区和j社区的地表沉降速率；wij则为空间矩阵中的元素。

图4 PS点累积沉降值分布图

采用地表沉降均值和中位数计算的675个社区的地表沉降的Moran'sI指数均达到了0.45以上，且p值小于0.01，这说明研究区内675个社区的地表沉降统计值在空间分布上有较强的集聚效应。

4.2 聚类分析

Moran'sI指数从全局的角度对研究区域内社区单元的地表沉降统计值的空间分布规律进行评估，评估结果很大程度上掩盖了局部反常和规律的不确定性问题。因此本文采用AnselinLocal Moran'sI进行空间格局分析，AnselinLocal Moran'sI又称为聚类与异常值分析，该算法是Anselin基于Global Moran'sI算法提出的一个用于识别局部聚类和空间异常值的算法。AnselinLocal Moran'sI队点集内的每个点都计算其对应的Moran'sI指数和z值，并对计算结果进行空间分布特征的分类，分为四种情况：HH型高值聚类，又称为“热点”，即改点周围都是和它一类的高值点；LL型低值聚类，又称为“冷点”，即改点周围都是和它一类的低值点；HL型异常值,高值点被值点包围；LH型异常值，低值点被高值点包围。

将全市的675个社区内的PS监测点数据的地表累积沉降平均值、中位数分别计算。结果显示置信度为 95% 的统计显著性聚类和异常值如图5、图6所示。热点分析结果表明全市总体上沉降热点有三个片区：宝安-南山-福田片区、大鹏片区、龙岗片区。

图5 累积沉降中位数局部空间自相关分布图

图6 累积沉降中位数局部热点分布图

聚类分析结果表明全市共有170 个社区沉降指标的分析结果为HH，主要分布在南山区、福田区、罗湖区等，这些区域是深圳市城市建设、填海工程等密集的区域，说明了深圳的城市建设开发强度和市域范围的地面沉降有较强的相关性。

5 结束语

本文利用PS干涉测量技术（PS-InSAR）获取了2016年5月至2020年8月的深圳的地表沉降数据，利用莫兰指数（Moran'sI）热点分析 （Getis-Ord Gi*）对研究区地表沉降格局特征进行分析。得出以下几个结论：

（1）采用区域统计计算的各社区地表沉降均值和中位数进行计算的研究区地表沉降的Moran'sI指数均达到了0.45以上，且p值小于0.01，这说明区域分布上675个社区的地表沉降分布呈现出较强的空间集聚现象。

（2）将全市的675个社区内的PS监测点数据的地表累积沉降平均值、中位数分别计算。结果显示置信度为 95%热点分析结果表明全市共有170个社区沉降指标的分析结果为HH，说明深圳的地表沉降空间上集中分布在城市的高强度开发区域。