沿海城市时序InSAR沉降监测分析方法、装置、设备及介质

摘要

本发明提供了一种沿海城市时序I nSAR沉降监测分析方法、装置、设备及介质，包括获取待处理数据集，并对待处理数据集进行预处理，生成输入参数集；对输入参数集进行PS‑I nSAR形变处理，生成PS‑I nSAR沉降结果；对输入参数集进行SBAS形变处理，生成SBAS沉降结果；对PS‑I nSAR和SBAS沉降结果进行比较分析处理，生成形变结果，从形变空间分布和形变量方面对沉降结果进行对比验证，获得同时间段、不同解算方法的形变信息；结合待处理数据集的沉降导致因素和形变结果进行沉降相关性分析，生成沉降形变速率，根据沉降形变速率确定城市中的沉降速率加速地区，以使研究人员对沉降速率加速地区进行监控。此外，传统差分干涉技术在持续缓慢形变的监测中，存在着精度低，结果可靠性差的缺陷。

说明书

技术领域

本发明涉及合成孔径雷达干涉测量技术领域，具体涉及沿海城市时序 InSAR沉降监测分析方法、装置、设备及介质。

背景技术

地面沉降是在自然因素和人为因素的共同作用下，所产生的一种地表高程缓慢降低的地质现象。城市地表的监测与预防直接关系到城市环境和社会经济的可持续发展，然而近年来，由于地下水抽取、地下空间开发、冲填土的压实、地表负载过重等引起的地面沉降和基础设施形变问题日益严重，造成建筑物的结构变形和损毁，其所带来的安全隐患和经济损失越来越大。截至2011年12月，中国已有50余个城市出现地面沉降，地面沉降这一问题已经遍布全国，例如上海市从1921年至1965年最大累计沉降量已达2.63m，沉降区域达400km2；天津市目前最大累计沉降量已达2.5m，沉降100mm以上的受灾区达900km2；西安市的地面沉降最为严重，高陵县、三原县和礼泉县地表比较稳定，其余7个县均都存在不同程度的地面沉降。因此，对城市地表沉降进行精确、快速和低成本的监测已成为重要方向，对监测结果进行有效的分析，在此分析的基础上，有关部门可以做出正确的防范措施。

目前，在经济快速发展与土地资源有限的背景下，围填海活动成为沿海地区缓解土地供求矛盾、拓展社会生存空间、扩大经济发展空间的有效手段。同时，随着我国城市建设的不断发展，一大批城市基础设施陆续建成，如轨道交通、地下管网，快速路，桥梁等。城市基础设施是人民生产生活的重要保证，同时也面临地面沉降及基础设施形变等严重问题。因此，针对沿海地区，在围填海造地和城市基础设施大量建设的前提下，如何建立面向沿海城市普遍存在的由填海造地、地铁、高速道路兴建导致的沉降监测方案及沉降成因分析是具有非常重要的应用价值。

由于传统的地面沉降测量手段存在成本高，无法持续监测等问题，星载合成孔径雷达干涉(Interferometric Synthetic Aperture Radar，InSAR) 作为一种新兴的测量技术，以其范围广，精度高，全天时等特点，可实现对城市基础设施持续性观测。InSAR相位包含着高程与形变相关信息，因此 InSAR技术起初用来提取大范围内的高程信息，后也广泛应用于地表沉降监测——利用星载雷达的系统参数和成像几何关系反演出某一点的三维空间位置及其微小变化。目前，该技术被广泛应用于城市沉降监测活动过程中，但因受制于时空失相干、大气延时相位、数字高程模型(Digital Elevation Model,DEM)误差等因素，现有的传统差分干涉(Differential InSAR,D-InSAR)技术在持续缓慢形变的监测中，存在着监测精度低，监测结果可靠性差的缺陷。

有鉴于此，提出本申请。

发明内容

本发明公开了一种沿海城市时序InSAR沉降监测分析方法、装置、设备及介质，能够有效解决传统差分干涉技术在持续缓慢形变的监测中，存在着监测精度低，监测结果可靠性差的问题。

本发明公开了一种沿海城市时序InSAR沉降监测分析方法,包括：

获取待处理数据集，并对所述待处理数据集进行预处理，生成输入参数集；

对所述输入参数集进行PS-InSAR形变处理，生成PS-InSAR沉降结果；

对所述输入参数集进行SBAS形变处理，生成SBAS沉降结果；

对所述PS-InSAR沉降结果和所述SBAS沉降结果进行比较分析处理，生成形变结果，其中，从形变空间分布和形变量两个方面对所述PS-InSAR 沉降结果和所述SBAS沉降结果进行对比验证，以获得同时间段、不同解算方法的形变信息；

结合所述待处理数据集的沉降导致因素和所述形变结果进行沉降相关性分析，生成沉降形变速率，并根据所述沉降形变速率确定城市中的沉降速率加速地区，以使得研究人员对所述沉降速率加速地区进行监控。

本发明还公开了一种沿海城市时序InSAR沉降监测分析装置，包括：

数据获取单元，用于获取待处理数据集，并对所述待处理数据集进行预处理，生成输入参数集；

PS-InSAR沉降结果获取单元，用于对所述输入参数集进行PS-InSAR形变处理，生成PS-InSAR沉降结果；

SBAS沉降结果获取单元，用于对所述输入参数集进行SBAS形变处理，生成SBAS沉降结果；

沉降结果比较单元，用于对所述PS-InSAR沉降结果和所述SBAS沉降结果进行比较分析处理，生成形变结果，其中，从形变空间分布和形变量两个方面对所述PS-InSAR沉降结果和所述SBAS沉降结果进行对比验证，以获得同时间段、不同解算方法的形变信息；

沉降因素分析单元，用于结合所述待处理数据集的针对沿海城市易导致沉降的因素和所述形变结果进行沉降相关性分析，并根据沉降形变速率确定城市中的沉降速率加速地区，以使得研究人员对所述沉降速率加速地区进行监控。

本发明还公开了一种沿海城市时序InSAR沉降监测分析设备，包括处理器、存储器以及存储在所述存储器中且被配置由所述处理器执行的计算机程序，所述处理器执行所述计算机程序时实现如上任意一项所述的沿海城市时序InSAR沉降监测分析方法。

本发明还公开了一种可读存储介质，存储有计算机程序，所述计算机程序能够被该存储介质所在设备的处理器执行，以实现如上任意一项所述的沿海城市时序InSAR沉降监测分析方法。

综上所述，本实施例提供的一种沿海城市时序InSAR沉降监测分析方法、装置、设备及介质，所述时序InSAR技术沉降监测分析方法通过选用 SAR卫星数据，结合时序InSAR技术对快速城市化地区进行沉降反演，分别利用PS-InSAR和SBAS两种时序InSAR技术对沉降结果进行对比验证，减少了传统人工水准测量的工作量；并对沿海城市沉降的因素和根据时序 InSAR反演获得的沉降量进行分析，得到沿海城市沉降变化及成因特点。从而解决现有的传统差分干涉技术在持续缓慢形变的监测中，存在着监测精度低，监测结果可靠性差的问题。

附图说明

图1是本发明提供的沿海城市时序InSAR沉降监测分析方法的流程示意图。

图2是本发明实施例提供的研究区范围示意图。

图3是本发明实施例提供的POD精轨文件名示例示意图。

图4是本发明实施例提供的覆盖厦门市的DEM--SRTM v3.0示意图。

图5是本发明实施例提供的厦门本岛建筑物等级情况示意图。

图6是本发明实施例提供的厦门本岛地质简图(含填海区)。

图7是本发明实施例提供的地层类型与沉降率的箱型示意图。

图8是本发明实施例提供的厦门岛典型沉降区(SBAS结果)示意图。

图9是本发明实施例提供的A区沉降情况示意图。

图10是本发明实施例提供的B区沉降情况示意图。

图11是本发明实施例提供的C沉降区情况示意图。

图12是本发明实施例提供的厦门岛地铁300米缓冲区沉降示意图。

图13是本发明实施例提供的厦门岛部分地铁站的时序沉降量示意图。

图14是本发明提供的沿海城市时序InSAR沉降监测分析装置的模块示意图。

图15是本发明实施例提供的Sentinel-1A卫星的SAR数据列表示意图。

图16是本发明实施例提供的PSInSAR方法情况示意图。

图17是本发明实施例提供的SBAS方法情况示意图。

图18是本发明实施例提供的厦门岛形变结果示意图。

具体实施方式

为使本发明实施方式的目的、技术方案和优点更加清楚，下面将结合本发明实施方式中的附图，对本发明实施方式中的技术方案进行清楚、完整地描述，显然，所描述的实施方式是本发明一部分实施方式，而不是全部的实施方式。基于本发明中的实施方式，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施方式，都属于本发明保护的范围。因此，以下对在附图中提供的本发明的实施方式的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施方式。基于本发明中的实施方式，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施方式，都属于本发明保护的范围。

以下结合附图对本发明的具体实施例做详细说明。

请参阅图1，本发明公开了的第一实施例提供了一种实施例提供的,包括：

S101，获取待处理数据集，并对所述待处理数据集进行预处理，生成输入参数集。

具体地，步骤S101包括：对所述待处理数据集进行数据导入处理，生成预设标准数据格式的导入数据，其中，所述待处理数据集包括单视复数影像数据、DEM数据和轨道数据，在导入过程中，采用精轨数据对所述轨道数据进行矫正，并根据卫星入射角计算多视参数；对所述导入数据进行图像裁剪处理，生成输入参数集，其中，基于SAR坐标系按预设研究区域范围对所述导入数据进行裁剪，以生成输入参数集。

在本实施例中，可以采集星载SAR遥感影像Sentinel-1 A/B IW数据进行沉降监测分析；首先根据监测范围坐标，选择并获取20景以上的星载 SAR的单视复数影像数据，以及对应监测区域的DEM数据和对应SAR数据的轨道数据；选择合成孔径雷达处理软件，进行数据导入、多视、滤波、裁剪、配准等数据预处理工作，便于后续的差分干涉相位图生成。以Envi软件的SARscape插件为例，进行sentinel-1A SLC单视复数数据和对应DEM 数据导入，获得Sarscape软件能够处理的数据格式；同时，由于 Sentinel-1A/B数据幅宽较大，为提高运行效率，利用SARscape软件利用行列号进行裁剪，获得研究区范围。

具体地，先将数据导入成SARscape标准数据格式，在数据导入过程中加入精轨数据进行轨道精矫正，使用卫星精密轨道数据对轨道信息进行修正，可有效去除因轨道误差引起的系统性误差。针对Sentinel1 A/B卫星，使用欧空局会提供的POD精密定轨星历数据(POD,Precise Orbit Ephemerides)对轨道信息进行修正，去除因轨道误差引起的系统性误差；同时，在数据导入的过程中，设置地面采样间隔，根据卫星入射角计算多视参数，进行多视处理，并基于SAR坐标系按研究区域范围对SLC文件进行裁剪。感兴趣区域的框选则可以多选一些像素，如缓冲区2～3公里，至少200个像素，以避免处理过程中出现明显的边界效应。

厦门市地处福建省东南角，是一个高人口密度的海湾型城市，全市陆地总面积约为1579km，如图2(a)所示。面对日益严重的“土地赤字”问题，向海要地成为厦门地区发展的重要手段，如图2(b)所示。由于厦门市属亚热带海洋性季风气候，强降水天气多发，经常突发路面塌陷的情况，如图2(c)所示。

具体地，在本实施例中，以厦门岛为例，对时序InSAR技术沉降监测分析方法进行说明。获取厦门岛的Sentinel-1时序数据，其中，entinel 系列卫星是欧洲航天局哥白尼计划中的地球观测卫星，Sentinel-1卫星搭载着欧洲极地轨道C波段雷达成像系统，具有四种传感器模式，单颗卫星重返周期为12天，双星重访周期缩短至6天。如下表1所示：

表1：

选用Sentinel-1A IW升轨数据作为SAR数据源，下载了时间跨度从2017 年3月14日至2018年3月9日共计29景影像，如图15所示，在InSAR 处理过程中，若轨道信息存在误差，则会造成基线误差，基线误差会以残差条纹的方式存在于干涉图中，进而造成系统性误差。因此，必须使用卫星精密轨道数据对轨道信息进行修正，以去除可能出现的系统性误差。 Sentinel-1所使用的是POD精密轨道数据，POD文件是在数据获取后21天内即可下载，其精度在5cm以内。POD轨道数据文件名中带有3个日期，从左至右依次是：精密轨道数据发布日期；数据成像时间的前一天；数据成像时间的后一天。因此，POD轨道数据文件名应为2016年9月15日哨兵数据的精轨文件，如图3所示。

通过获取DEM数据、道路数据、建筑物数据和厦门岛的地质图作为用于最后分析沉降结果的辅助数据。首先，获取DEM数据；由于在利用InSAR 技术监测地表沉降情况时，需要引入外部的DEM，用来消除地形信息而获取形变信息。在本实施例中，采用美国国家航空航天局航天飞机雷达地形测绘任务在2015年所发布的DEM，其分辨率为30米。其中，厦门市高程介于 -33至1165米，请参阅图4，此为厦门市的大地高程，即地面点距离WGS84 椭球体表面的高程数值。

其次，获取道路数据；路网承载着物质世界中的能量与物质交换，与人们的生活紧密相关，因此探究路网与地表沉降的内在关系，能够对道路健康状况进行实时监控并及时预警显得尤为重要。在本实施例中，将获取到的线状矢量数据作为道路数据，其包含属性字段——道路等级，我国将公路分为以下5个等级，将省道、国道归为一级公路，其中，根据统计结果表明，厦门市道路总长达14727468.25米。如下表2所示：

表2：

再次，获取建筑物数据；现如今，钢筋混凝土成为了建筑物的主要材料，并且建筑物越建越高，导致地表载荷越来越大。在本实施例中，获取到的建筑数据为面状矢量数据，其包含了楼高属性字段。将建筑物数据分为四种等级，分别为多层(1～7层)、小高层(8～18层)、高层(19～32层)与超高层(大于32层)，如图5所示。根据统计结果表明，厦门岛总建筑面积达24.89km2，占厦门市建筑物总面积的61.98％，建筑物密度高达15.78％，平均楼高4层。

最后，获取地质数据；使用ArcGIS软件加载此Web地图切片服务中的1:25万建造构造图G50C004003厦门市幅，以研究区为范围进行矢量化，最终获得厦门岛的地质数据，如图6所示。再结合研究区界线，叠置出厦门岛近年的填海区域，以便于后期分析。其中，厦门岛地质所涉及地层的按从古至今的时间顺序为：三叠系地层、侏罗纪地层、燕山期花岗岩地层、第四纪地层、填海区地层(属于第四纪期)。

S102，对所述输入参数集进行PS-InSAR形变处理，生成PS-InSAR沉降结果。

具体地，步骤S102包括：对所述输入参数集进行处理，生成连接图，其中，在输入的N副单视复数影像数据中定义一副主影像，并根据所述主影像生成N-1个主从相对，根据所述主影像和所述主从相对生成连接图；

对所述主从相对进行干涉处理，生成去已知地形后的干涉图和SAR影像的强度数据，其中，对每一所述主从相对和进行配准处理，生成干涉图，并利用所述DEM数据和获取到的大地控制点对所述干涉图进行处理，生成去已知地形后的干涉图和SAR影像的强度数据；

根据所述干涉图的振幅离差指数和相干性选定形变反演结果，即PS点，并利用线性模型从所有干涉图中估算出PS点的残余地形和形变速率；

通过低通、高通滤波将所述残余地形和所述形变速率中的大气信息提出后剔除，提高PS方法的精度；

将所述PS点从SAR坐标系中转换成地理坐标系或投影坐标系中，生成PS-InSAR沉降结果。

具体地，在本实施例中，通过PS-InSAR技术获取形变速率；首先，输入N幅SLC数据，定义一幅主影像，并根据主影像生成N-1个主从相对，最后生成了连接图，其中，主影像是指以该影像为标准，其他影像都配准到该影像上。其次，在干涉过程中为每一对主从影像进行配准、干涉图生成、去平地效应、振幅差指数计算，生成去已知地形后的干涉图与SAR影像的强度数据。其中，此干涉过程需要DEM与GCP文件这两个辅助数据，并利用DEM数据消除干涉处理过程中的平地效应；大地控制点的选点对轨道进行校正、主影像和DEM进行几何校正，以增加配准精度。再次，根据振幅离差指数、相干性初步选定最佳PS点，即处于能被SAR天线探测到的方向、波动小于1mm/y，并用线性模型从所有差分干涉图中估算出PS点的残余地形和形变速率，以便在接下来推算出其他PS点的结果，其中，所述结果包含大气信息的相干性图、差分干涉图、沉降速率与高程信息。从次，大气因素是传统D-InSAR技术停留在分米级到毫米级精度的主要原因之一，在本实施例中，通过低通、高通滤波将大气信息提出后剔除，提高PS方法的精度，以剔除大气信息后，基于最佳PS点进一步解算其他PS点的信息，反演结果为PS散射点，每个永久散射点具有时序变化的形变速率。最后，将形变反演结果即PS点从SAR坐标系转换成地理坐标系或投影坐标系，方便在后续叠加地图底图，进行监测区域时序形变特征分析。

在本实施例中，所述PS-InSAR技术的核心思想是利用永久性散射体相位特征，如高楼、桥梁、裸露的岩石等，对这些PS目标点的相位特性进行建模和解算，能以mm级监测精度获取地表在相应时段的缓慢形变过程和特征，已在国内外多个城市及地区地表形变中得到广泛的应用。其中，所述 PS-InSAR技术较传统的DInSAR技术可以减少时间失相干和大气延迟相位的影响，尤其适用于沿海城市地区的形变监测。

具体地，在本实施例中，以上述厦门岛为例：利用PS-InSAR技术对 sentinel-1数据进行沉降反演，其生成的连接图，如图16所示；在连通图生成的过程中，将基线阈值设置为400％，并未出现超出阈值的单幅影像；由于sentinel-1 IW成像方式中方位向和距离向分辨率为5m*20m，故在进行多视处理时，距离向设置为1倍，方位向设置为4倍。在干涉图生成的过程中，距离向上会过采样4倍，如此可避免由于基线太长而产生快速变化的干涉条纹；楼层数按3.3米计算，厦门岛平均楼高约15米，故设定高程残差值为15米；根据多次反复实验，研究区形变范围大概介于-30mm/y 至30mm/y，因此以30为阈值来预估研究区的形变范围；大气高通滤波阈值设置为365天；大气低通滤波设置为1200米，若研究区内有较多植被覆盖，可减少低通滤波的设置。

S103，对所述输入参数集进行SBAS形变处理，生成SBAS沉降结果。

具体地，步骤S103包括：根据预设的空间基线阈值与时间基线阈值，自动从所述输入参数集中生成多幅主从影像对；

将所有影像对配准到超级主影像上，并对每一所述主从影像对进行干涉处理，生成多对干涉结果；

对每一所述干涉结果进行检查，剔除相干性结果较差的所述主从影像对；

获取输入的GCP数据，并设置为轨道精炼的控制点；

对所述主从影像对进行第一次反演处理，估算形变速率和残余地形；

利用大气高通、大气低通对所述主从影像对进行第二次反演处理，生成大气因素，并将每一所述干涉结果都减去所述大气因素；

将不同时序相同区域内两幅所述主从影像对进行共轭相乘后相减，生成地理坐标系下的空间连续的形变结果。

具体地，在本实施例中，通过SBAS获取形变速率；首先，根据设置的空间基线阈值与时间基线阈值，自动生成多幅主从影像对，并反复设置连接阈值，避免生成的像对过多而导致数据冗余。其次，将所有影像配准到超级主影像上，对每一个主从相对进行干涉处理、相干性生成、去平、滤波和相位解缠，对干涉结果进行检查，将相干性结果较差的主从相对，从连接图中去除。再次，将GCP点作为轨道精炼的控制点，并且在接下来的解算中其形变值被定义为0mm，因此GCP点将直接影响着SBAS技术结果的准确性；其中，主要选用相对稳定的、没有残余地形条纹、没有形变条纹、远离形变区域、没有相位跃变的GCP点。从次，进行第一次反演估算形变速率和残余地形；进行第二次反演，并利用大气高通、大气低通对大气影响进行估计，将每一幅的形变结果都减去这些大气因素。最后，将不同时序相同区域内两幅影像进行共轭相乘后相减，但考虑了时序上的多幅影像，并且主从影像对极其丰富，提取的形变结果为栅格数据，能够获得较为连续的沉降信息，结合DEM数据进行地理编码，以获得地理坐标系下的空间连续的形变结果。其中，所述SBAS技术的核心思想是通过设置时空基线阈值将满足时空基线阈值的影像两两进行差分干涉处理，利用这些具有较短时空基线的影像产生具有较高相干性的差分干涉图，通过对差分干涉图的多视处理、滤波处理和相位解缠处理，提取高相干点，利用奇异值分解得到地表形变速率的最小范式最小二乘解，最终得到长时间形变时间序列。其中，所述SBAS技术较传统的DInSAR技术可以减少时间失相干和大气延迟相位的影响，尤其适用于沿海城市地区的形变监测。

具体地，在本实施例中，以上述厦门岛为例：利用SBAS技术获取厦门岛沉降信息，分别进行沉降信息提取，过程数据的示例，如图17所示，由于Sentinel数据的基线较短，若两两连接的空间基线阈值设太宽，会造成较大的数据冗余，经过反复试验，将临界基线最大百分比设置为1％； Sentinel-1的重返周期为12天，故将最大时间基线设置为120天；在干涉过程中，选择3D解缠，即满足3D解缠的像对将进行3D解缠，以增加精度；解缠方法选用最小费用流法；滤波方法选用Goldstein；线性模型稳定性较强、适用性较广，故而选用线性模型；大气高通滤波阈值设置为365天；大气低通滤波设置为1200m。在轨道精炼的过程中，选取一个能代表绝大多数像对情况的像对，用其解缠结果、滤波后的差分干涉结果作为参考图，选取控制点，这些控制点的形变值将定义为零，这些控制点直接影响着结果的准确性。

S104，对所述PS-InSAR和所述SBAS沉降结果进行比较分析处理，生成形变结果，其中，从形变空间分布和形变量两个方面对沉降结果进行对比验证，以获得同时间段、不同解算方法的形变信息。

具体地，步骤S104包括：对所述PS-InSAR沉降结果和所述SBAS沉降结果进行反演结果评价，生成反演异同结果，其中，对所述PS-InSAR沉降结果和所述SBAS沉降结果的沉降数值进行直方图统计，获得两个沉降结果的标准差、平均值、最大值、最小值以及直方图形态；

对所述PS-InSAR沉降结果和所述SBAS沉降结果的精度图进行相同量级程度的颜色渲染，比较二者在空间分布上的差异，生成空间差异结果；对所述PS-InSAR沉降结果和所述SBAS沉降结果中沉降明显区域进行随机抽样处理，建立两种反演结果的变化趋势以及变化精度差异，其中，选择若干个控制点，分别获取在两种方法反演结果中，选择的控制点的沉降形变信息；结合所述反演异同结果、所述空间差异结果、所述变化趋势以及所述变化精度差异，生成形变结果。

现市面上的基于地面水准测量数据的对比验证方法是对沉降结果进行精度验证的常用方法，在时间序列的角度上，水准测量结果往往间隔数月， InSAR的沉降结果取决于卫星的重返周期，因此后者的时间频率高于前者；同时，水准测量以点形式表现，通常无法获取大范围的水准测量结果，水准测量的控制点的密度远远小于InSAR反演结果，并且两者的点位不能够完全对应，若使用空间插值的方法解决此问题，会造成精度的丢失与准确性下降，若使用最近点匹配的方法解决此问题，会影响精度验证的结果。

具体地，在本实施例中，在对比验证的过程中，分别从定性和定量两个方面进行反演结果精度对比。可以利用SAR处理软件自带的精度评价模块，进行反演结果评价，其能够导出精度评价模型结果，从精度量级上对两种反演结果的精度进行对比，对沉降数值进行直方图统计，分别统计两种方法获得沉降结果的标准差、平均值、最大值、最小值以及直方图形态，对比反演精度结果的异同。从空间分布格局上，可以利用ArcGIS软件对获得沉降反演结果的精度图进行相同量级程度的颜色渲染，比较二者在空间分布上的差异。对反演结果沉降明显区域随机抽样，选择若干个控制点，分别获取在两种方法反演结果中，选择的控制点的沉降形变信息，建立两种反演结果的变化趋势以及变化精度差异；对比分析结果可以验证利用两种时序InSAR技术获取沉降信息的可靠性，若两种方法获取的沉降信息在空间分布、沉降量级上相似度很高，则可确定为该区域的沉降结果是可靠的，以供后续对沿海区域沉降因素进行分析，寻找易沉降典型地区，在城市管理中加以监控。因此，将利用PSInSAR、SBAS方法对研究区域分别获得沉降反演结果进行对比，若两种方法获得的反演结果量级相似，且空间分布范围相同，进而可以证明获取的研究区沉降信息相对准确。

具体地，在本实施例中，以上述厦门岛为例：利用PS-InSAR、SBAS技术对厦门岛进行沉降反演，最终获得了同时间段、两种解算方法的毫米级的形变结果，两种结果的直方图都趋于正态曲线，如图18所示，其峰值均接近0，说明两个提取结果较为理想。根据形变率数值，将厦门岛的沉降状况进行分级，以箱型图的下四分位数-3.55mm/y、箱型图的最小值 -11.75mm/y为断点，将厦门岛的沉降区域分为三个等级，等级越高沉降越严重，根据表7可以看出，厦门岛共有40.74％的地区在不同程度上存在地表沉降，并主要集中在西部区域。如表3所示：

表3：

在本实施例中，由表3可以看出根据箱型图，厦门岛的形变量主要聚集在-10mm/y至10mm/y之间，总体沉降量较小；从整体上来看，除无法检测到的地区外，PS-InSAR结果中沉降的PS点个数占总数的52.34％，SBAS结果中沉降面积占总面积的60.25％，说明除无法检测到的地区外，厦门沉降地区面积占比较多，城市在整体上存在一定的沉降风险；从空间分布来看，东部沿海区域总体呈抬升趋势，形变值最大达27.44mm/y，这可能与厦门地区的板块运动有关；西部区域沉降严重，最大沉降率达-38.06mm/y；云顶隧道北出口及其附近区域沉降严重，最大沉降率达-19.17mm/y，其他地区也出现不同程度的沉降；从形变量来看，形变信息密集分布在建筑物、道路等区域，对PS-InSAR沉降结果进行统计，建筑物、道路等区域PS点占比高达82％；而在植被与水体覆盖区域较少，这与该区域受植被和水体的影响，干涉相对的相干性较低的缘故有关；零星分布在植被与水体区域间的 PS点位于人工建筑与裸地等处，该形变PS点的分布结果与永久散射体方法的特点有关。对比两种方法，PS方法能够获得大量的形变信息，共有6× 105个相干性高的PS点被提取出来，SBAS结果的像素个数有5×105个，二者沉降反演结果在数值和空间上均具备一致性。在空间表现上，PS获取的沉降信息以PS点位展示，结果是点矢量数据；SBAS反演的形变结果是面状结果，空间表现上是连续的，因此，利用SBAS反演结果进行沉降空间格局分析。

S105，结合所述待处理数据集的针对沿海城市易导致沉降的因素和所述形变结果进行沉降相关性分析，并根据沉降形变速率确定城市中的沉降速率加速地区，以使得研究人员对所述沉降速率加速地区进行监控。

具体地，步骤S105包括：梳理所述待处理数据集的沉降导致因素和沉降量，确定主要因素；

统计所述主要因素与所述沉降量，并与所述形变结果进行相关性分析，确定沉降量与主要因素之间的关系；

对所述形变结果进行制图处理，生成形变结果图；

确实所述形变结果图中沉降量较大的区域，对所述沉降量较大的区域进行监测处理，生成沉降形变速率；

根据所述沉降形变速率确定城市中的沉降速率加速地区。具体为：对所述沉降形变速率进行时序数据变化处理，并监测所述沉降形变速率，确定沉降速率加速区域。

通常情况下，地表形变或沉降的监测结果一般用沉降量或沉降速率来描述，即形变过程在时间上累积结果或变化快慢来描述沉降过程；而在工程实际和灾害预防中，空间上均匀的沉降所造成的危害远小于不均匀沉降。

通常情况下，由于地质结构、地理位置等的差异，不同的城市沉降因素不同；为缓解土地资源紧缺压力，海岸带城市有不少填海造地工程。基于此，在海岸带城市的沉降因素分析中，将填海造地、地质条件作为影响城市沉降的首要考虑因素；同时，快速城市化过程中，存在大量的地下基础设施，以及大型建筑物导致地面载荷过量等也对城市沉降产生极大的影响。因此，将填海造地、地质底图、地铁线路图及其缓冲图、土地利用图及楼高数据作为影响沉降的参数，用于对沉降结果进行影响因素分析和沉降易发生区沉降结果分析。

具体地，在本实施例中，首先需要梳理易导致沿海城市沉降主要因素的基础底图，如填海造陆部分、地质图、大型建筑物即人类活动密集、地下基础设施的建设区用图等，并与沉降反演结果进行配准；在统一的空间框架下，以PS点为单位，统计易导致沉降的主要因素数值和沉降量，并对二者进行相关性分析，寻找沉降量与主要因素之间的关系；其次，对步骤 S104中的地面沉降结果进行制图，寻找沉降量较大的区域，提取该区域形变速率，在时序数据变化过程中，找到沉降速率加速区域；最后，对沉降加速区域，重点考察填海造陆形成的新成陆区、地下基础设施建设的缓冲单元，结合相关性分析，划分沉降重点关注区域，建立城市沉降风险等级，为重点区域沉降监控提供参考。

具体地，在本实施例中，以上述厦门岛为例：利用楼高数据对PS结果进行分区统计，分析楼高与沉降之间的相关关系，其结果如表4所示：

表4：

结果表明，小高层的沉降率最低，达3.85mm/y；其次是多层的建筑物，沉降率达4.12mm/y；随后，高层的建筑物，沉降率达5.59mm/y；超高层的沉降率最高，达5.70mm/y，可见四种等级建筑物的沉降率相差较小，而且与楼高的关系并不明确。但随着楼高的增大，单位面积内的年平均沉降率呈指数级别的增长，表明建筑物楼高与沉降成正相关关系。统计单位面积内形变率，多层建筑物最小，达0.19；其次是小高层建筑物，达1.57；接着是高层建筑物，达9.02；最后是超高层建筑物，单位面积内的形变量最高，高达95，说明对于在城市发展过程中，高层建筑导致的地面过载，会导致一定的沉降隐患。根据地层类型，对PS点沉降率进行统计，分析地层类型与研究区沉降的关系。各个地层类型的沉降率分散情况，如图7所示；其中，填海区的沉降最为严重。从箱型图中位数的角度分析，侏罗纪地层覆盖区域，侏罗纪地层距今年代较为久远，地层压密极好，且无构造带穿过，因此沉降率最小；三叠系地层覆盖区域，三叠系地层为该区域年代最为久远的地层，该时期岩层具有明显的三分性，而且其周边均存在断裂带，对该地层有极大影响，严重破坏其结构，所以该区域沉降率比较大；第四纪地层覆盖区域，第四纪地层沉积距今时间较近，地层密实程度不高，因此沉降率较大；燕山期花岗岩覆盖区域，燕山期花岗岩是侏罗纪至白垩纪之间侵入的岩浆岩，并裸露出地表，受风化作用较大，并且地区花岗岩节理裂隙发育，因此中位数比第四纪大，但从总体数据的分散情况而言，第四纪的沉降更大；填海区，填海区为人工填海所得，属于第四纪人工填积物，无论是中位数，还是整体数据分散情况，其沉降都最为严重，因此沿海城市应密切关注填海区的沉降情况。

在本实施例中，在2017至2018年期间，厦门岛出现了许多不同程度的沉降漏斗，如图8所示；结合建筑物数据、道路数据、地层类型数据、 POI信息及谷歌地球等资料获得各个区域的基本情况，如表5所示：

表5：

根据表5的结果显示，厦门岛整体沉降虽较为稳定，但也存在许多沉降漏斗，并且类型较多，主要发生在填海区、交通要道、高楼区、活动密集区等，同时也验证了沿海地区沉降的主要因素为填海成陆区、地面高负载区及地下基础设施建筑区。利用SBAS结果对A、B、C沉降区及厦门岛的两条地铁区域的沉降进行详细的空间分析与成因分析。

A沉降区位于厦门本岛西北部的一个小港口，均由填海所得，区域功能以商业服务用地为主，包含生鲜市场与海鲜城等。此区域在2017年至2018 年沉降严重，出现了一个明显的沉降漏斗，沉降率最大达38.06mm/y，最小为2.59mm/y，平均沉降率达27.17mm/y，该区域在长半轴上为一个U型的沉降漏斗，如图9(e)所示；剖面中点a1在2017至2018年内累积沉降量达41.9毫米。由沉降折线可知，如图9(f)所示，时序点a1、a2、a3在该时间段内发生持续沉降，三者在2017年7、8、10月有少量抬升，但总体而言三者的沉降趋势较为一致。其中，a1为建筑物地区，沉降速率为区域内的最大值38.06mm/y；a2亦为建筑物地区，沉降率为22.98mm/y；a3 为该区域的一个道路交叉路口，沉降率达22.35mm/y。由此可见，该区域建筑物区域沉降比道路区域沉降更为严重，如图9(c)、图9(d)所示，该区域在2015年10月为施工区，2016年5月基本完成，2017年1月道路已全部建设完成。因此，此区域的沉降与该地区压密活动密切相关，尤其为建筑物区域，地层压密活动更为剧烈，沉降率达到区域最大值。

B沉降区位于海峡博览中心区域，如图10(a)所示，该地为早期的填海区域。在观测的时间段内，此区域处于建设期，形变率介于-19.96mm/y 至11.99mm/y，如图10(b)、图10(c)、图10(d)所示。b1所在建筑物为厦门佰翔会展中心，于2016年左右开始破土动工，2017年3月建筑物主体结构已经完工。根据其沉降折线可知，如图10(d)所示，从2017年 7月开始，该建筑物地表开始发生沉降。结合A填海区中的分析，b1建筑物的工程规模与质量高于a1的建筑物，并且相比于a1时序点而言，b1处的沉降速率缓慢，沉降量较小，至2018年3月该地累积沉降量达10.1毫米，说明建筑物在建成之后确实存在不同程度的地层压密活动，压密活动的时间长短和程度可能与建筑物地基工程有关。虽然b1地区以沉降为主， b2地区以抬升为主，但两者在2017年的抬升趋势基本相同，在20170407、 20170501、20170606、20170712等四个时间点中b1、b2都有所抬升，这可能与该区域的风速、温度有关，存在周期性的摆动。b3在2017年3月至 2018年3月期间还处于施工状态，但其沉降值也较为明显，沉降率达 14.03mm/y，并出现了小范围内沉降漏斗，由此，毫米级InSAR技术可对地表施工情况进行监测，可利用该技术对非法建设等活动进行大范围的监控。

C区位于殿前一路与兴湖路形成的锐角区域内，四周的道路为沥青路面的三级道路，区域内的路网多为混泥土路面的四级道路，如图11(a)所示。该区域为商业服务用地，区域内布满餐厅、酒店、娱乐设施与停车场，该区在2017至2018年沉降严重，整个区域都发生了不同程度的沉降，沉降率介于-11.09mm/y至-1.57mm/y之间。L2线上跨多种地物类型，导致了曲线的曲折性，剖面线的低谷是该区域的主要道路，且其为四级道路，如图 11(c)所示。c1点位于厦门翔鹭大酒店中心，c2位于该区域的交通要道， c1在研究时间段内都发生沉降，而c2在2017年8月之前地表呈抬升状态，至2017年12月3日两者的沉降量才逐渐趋于一致。说明道路载荷过多会导致路面沉降速率变快，其沉降的数量级能够很快赶超建筑物的沉降量。此区域的沉降主要是由于该区域的商业服务活动活跃，车流量、人流量大，导致道路、建筑物的载荷过多，引发了区域性的沉降。由此可见，车流量可能是道路影响沉降的主要贡献因子。

利用SBAS结果对厦门岛内地铁沉降情况进行分析，在2017年3月至 2018年3月期间，厦门岛地铁300米缓冲区域内形变率，如图12所示，介于-24.71mm/y至10.09mm/y之间。地铁1号线于2017年12月31日开始试运营，线路上许多地方已发生不同程度的沉降，其中火炬园、莲坂、莲花路口和吕厝等地铁站沉降较为明显。请参阅图13，吕厝地铁站从2017年 12月份开始缓慢下沉，至2018年3月沉降量达5.73mm。在湖滨南路与湖滨中路交叉路口出现一个沉降漏斗，此为典型的交通要道沉降区，沉降率最大高达-21.88mm。地铁3号线一期于2015年年底开工，计划2020年底开通试运营，但在厦门火车站、体育中心、湖里公园、华荣路等地铁站上已发生不同程度的沉降，其可能是由地下施工所导致的。人才中心地铁站至火炬园地铁站路段发生了连续沉降，沉降从仙岳山出口处开始，至华荣路路段沉降较为严重，从华荣路之后沉降情况有所好转。华荣路地铁站沉降较为平稳，从2017年3月开始下沉，至2018年1月趋于稳定，停留在 8mm上下；但体育中心地铁站沉降较严重，从2017年12月开始，沉降速度加快，至2018年3月已下沉12.1mm；厦门火车站受温度、载荷度、客流量等多因子影响，因此沉降曲线较为复杂，但其在2017年12月份也开始下沉，至2018年3月沉降量达5.92mm。

综上所述，所述时序InSAR技术沉降监测分析方法，具有较高的时间和空间观测分辨率，随着技术的发展，将其应用于城市大区域范围的沉降监测具有极强的实际应用价值。其使用的PS-InSAR和SBAS技术相较于传统的DInSAR技术可以减少时空失相干和大气延时相位的影响，对两种技术获得的形变结果进行综合对比验证，以此确定研究区域沉降结果，克服了无法及时获取大面积实地水准观测的缺陷，提高了城市沉降实时监测能力；通过找出沉降发生的典型区域，并对该区域的沉降特点及诱发因素进行重点分析，为关键地区沉降监测提供了参考，对城市沉降预警提供了依据，具有很强的实用价值。

请参阅图14，本发明的第二实施例提供了一种沿海城市时序InSAR沉降监测分析装置，包括：

数据获取单元201，用于获取待处理数据集，并对所述待处理数据集进行预处理，生成输入参数集；

PS-InSAR沉降结果获取单元202，用于对所述输入参数集进行PS-InSAR 形变处理，生成PS-InSAR沉降结果；

SBAS沉降结果获取单元203，用于对所述输入参数集进行SBAS形变处理，生成SBAS沉降结果；

沉降结果比较单元204，用于对所述PS-InSAR沉降结果和所述SBAS沉降结果进行比较分析处理，生成形变结果，其中，从形变空间分布和形变量两个方面对所述PS-InSAR沉降结果和所述SBAS沉降结果进行对比验证，以获得同时间段、不同解算方法的形变信息；

沉降因素分析单元205，用于结合所述待处理数据集的针对沿海城市易导致沉降的因素和所述形变结果进行沉降相关性分析，并根据沉降形变速率确定城市中的沉降速率加速地区，以使得研究人员对所述沉降速率加速地区进行监控。

本发明的第三实施例提供了一种沿海城市时序InSAR沉降监测分析设备，包括处理器、存储器以及存储在所述存储器中且被配置由所述处理器执行的计算机程序，所述处理器执行所述计算机程序时实现如上任意一项所述的沿海城市时序InSAR沉降监测分析方法。

本发明的第四实施例提供了一种可读存储介质，存储有计算机程序，所述计算机程序能够被该存储介质所在设备的处理器执行，以实现如上任意一项所述的沿海城市时序InSAR沉降监测分析方法。

需说明的是，以上所描述的装置实施例仅仅是示意性的，其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本实施例方案的目的。另外，本发明提供的装置实施例附图中，模块之间的连接关系表示它们之间具有通信连接，具体可以实现为一条或多条通信总线或信号线。本领域普通技术人员在不付出创造性劳动的情况下，即可以理解并实施。

以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。