一种面向黄土填土沉降区的边界识别与提取方法

摘要

本发明公开了一种面向黄土填土沉降区的边界识别与提取方法，包括步骤：获取影像覆盖范围内黄土填土区的地表不均匀沉降数据；建立地面沉降点数据库；统计出每个沉降区域内沉降速率相同或相近的沉降观测点，将不同区域内所有等值点通过连接形成等值线；设定沉降速率和累计沉降量阈值，将大于一定阈值的沉降点视为位于填土区内的沉降点，相应的沉降点构成的等值线即可视为该区域内的填土区边界。本申请通过建立沉降区域的沉降观测点数据库，对沉降漏斗进行分类管理，形成了一个完整的黄土填土沉降区边界的识别提取技术体系，显著提高了黄土填土沉降区边界的识别提取能力，可以对填土沉降区进行科学有效的防治，及时发现问题，减少居民和建筑损失。

权利要求

1.一种面向黄土填土沉降区的边界识别与提取方法，其特征在于，包括以下步骤：S1、获取原始SAR影像数据并进行预处理，使用时间序列InSAR技术获取影像覆盖范围内黄土填土区的地表不均匀沉降数据；

S2、根据获得的地面沉降空间分布结果将沉降区划分为m个区域，其中，m为沉降漏斗个数；再按照划分的区域导出沉降区的沉降点的基本属性表，并以此建立地面沉降点数据库，将沉降区域进行分类管理；基本属性表包括形变速率表和累计形变量表；形变速率表中包含沉降观测点的坐标、对应的形变速率信息；累计形变量表中包含沉降观测点的坐标、SAR影像时间序列中各个时间段的累计形变量信息；

S3、在填土沉降区的基础地表形变数据库基础之上，统计出每个沉降区域内沉降速率相同或相近的沉降观测点，将不同区域内所有等值点通过连接形成等值线；

S4、设定沉降速率和累计沉降量阈值，将大于一定阈值的沉降点视为位于填土区内的沉降点，相应的沉降点构成的等值线即可视为该区域内的填土区边界；

步骤S1包括：

在影像预处理阶段，首先利用外部精密轨道数据完成轨道校正；其次选定主影像后利用外部数字高程模型数据完成地形校正；再将副影像配准到主影像几何，通过干涉生成干涉图；

步骤S1还包括：

获得影像覆盖范围内的黄土填土区地面沉降分布情况：首先设定时空基线阈值将预处理获得的干涉图构成干涉对基线网络，然后设定时空相干性阈值筛选相干点目标，以去除相干性较差的点；最后通过相位解缠得到干涉相位时间序列结果，通过网络反演得到研究区内的相应沉降结果信息。

2.根据权利要求1所述的一种面向黄土填土沉降区的边界识别与提取方法，其特征在于，每个时间段的沉降观测点都被赋予唯一的ID进行存储。

说明书

一种面向黄土填土沉降区的边界识别与提取方法

技术领域

[0001] 本发明涉及高原地表形变监测和黄土填土区边界识别技术领域，尤其涉及一种面向黄土填土沉降区的边界识别与提取方法。

背景技术

[0002] 随着现代城市的快速扩张，许多位于黄土高原上的城市通过削山填沟的方式来解决城市用地短缺的问题。近年来，为了有效的解决城市人口快速膨胀但城市用地十分短缺的问题，中国黄土高原上的许多城市都采取了削山填沟的方式来增加城市用地，如延安、兰州等。由于新建城市基本位于填方区域上方，黄土地基薄弱，会不可避免发生黄土湿陷等情况，进而出现严重的地面沉降问题，尤其是不均匀地面沉降存在重大的安全隐患。因此，为了保障新修建筑和居民的安全，确保地面的稳定性，十分有必要对黄土填土区的边界信息进行及时识别提取和科学管理规划。

[0003] 目前，填土沉降区的边界识别方法主要包括基于点目标的全球定位系统(GPS)方法和传统的人口踏勘方法两种。传统方法需要投入大量的人力物力，耗时耗力，特别是对于面积大，地形复杂的区域来说缺点尤为明显。

[0004] 因此需要研发出一种面向黄土填土沉降区的边界识别与提取方法来解决上述问题。

发明内容

[0005] 本发明的目的就在于为了解决上述问题设计了一种面向黄土填土沉降区的边界识别与提取方法。

[0006] 本发明通过以下技术方案来实现上述目的：

[0007] 一种面向黄土填土沉降区的边界识别与提取方法，包括以下步骤：

[0008] S1、获取原始SAR影像数据并进行预处理，使用时间序列InSAR技术获取影像覆盖范围内黄土填土区的地表不均匀沉降数据；

[0009] S2、根据获得的地面沉降空间分布结果将沉降区划分为m个区域，其中，m为沉降漏斗个数；再按照划分的区域导出沉降区的沉降点的基本属性表，并以此建立地面沉降点数据库，将沉降区域进行分类管理；

[0010] S3、在填土沉降区的基础地表形变数据库基础之上，统计出每个沉降区域内沉降速率相同或相近的沉降观测点，将不同区域内所有等值点通过连接形成等值线；

[0011] S4、设定沉降速率和累计沉降量阈值，将大于一定阈值的沉降点视为位于填土区内的沉降点，相应的沉降点构成的等值线即可视为该区域内的填土区边界。

[0012] 具体地，步骤S1包括：

[0013] 在影像预处理阶段，首先利用外部精密轨道数据完成轨道校正；其次选定主影像后利用外部数字高程模型数据完成地形校正；再将副影像配准到主影像几何，通过干涉生成干涉图。

[0014] 进一步地，步骤S1还包括：

[0015] 获得影像覆盖范围内的黄土填土区地面沉降分布情况：首先设定时空基线阈值将预处理获得的干涉图构成干涉对基线网络，然后设定时空相干性阈值筛选相干点目标，以去除相干性较差的点；最后通过相位解缠得到干涉相位时间序列结果，通过网络反演得到研究区内的相应沉降结果信息。

[0016] 优选地，基本属性表包括形变速率表和累计形变量表；形变速率表中包含沉降观测点的坐标、对应的形变速率信息；累计形变量表中包含沉降观测点的坐标、SAR影像时间序列中各个时间段的累计形变量信息。

[0017] 具体地，每个时间段的沉降观测点都被赋予唯一的ID进行存储。

[0018] 本发明的有益效果在于：

[0019] 本申请通过建立沉降区域的沉降观测点数据库，对沉降漏斗进行分类管理，形成了一个完整的黄土填土沉降区边界的识别提取技术体系；本申请投入的人力物力、时间显著少于现有技术，实际可操作性强，显著提高了黄土填土沉降区边界的识别提取能力，可以对填土沉降区进行科学有效的防治，及时发现问题，减少居民和建筑损失。

附图说明

[0020] 图1为本申请的流程框图；

[0021] 图2是本申请实施例中研究时间内的平均形变速率；

[0022] 图3是本申请实施例中基础地面沉降数据库图；

[0023] 图4是本申请实施例中填土沉降区边界识别提取结果示意图；

[0024] 图5是光学历史影像判别对比验证图。

具体实施方式

[0025] 为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。

[0026] 因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

[0027] 应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。

[0028] 在本发明的描述中，需要理解的是，术语“上”、“下”、“内”、“外”、“左”、“右”等指示的方位或位置关系为基于附图所示的方位或位置关系，或者是该发明产品使用时惯常摆放的方位或位置关系，或者是本领域技术人员惯常理解的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的设备或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

[0029] 此外，术语“第一”、“第二”等仅用于区分描述，而不能理解为指示或暗示相对重要性。

[0030] 在本发明的描述中，还需要说明的是，除非另有明确的规定和限定，“设置”、“连接”等术语应做广义理解，例如，“连接”可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接连接，也可以通过中间媒介间接连接，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。

[0031] 下面结合附图，对本发明的具体实施方式进行详细说明。

[0032] 如图1所示，一种面向黄土填土沉降区的边界识别与提取方法，包括以下步骤：

[0033] S1、获取原始SAR影像数据并进行预处理，使用时间序列InSAR技术(TS‑InSAR)获取影像覆盖范围黄土填土区的地表不均匀沉降数据。

[0034] 在影像预处理阶段，利用外部精密轨道数据完成轨道校正，其次选定主影像后利用外部数字高程模型(DEM)数据完成地形校正，再将副影像配准到主影像几何，通过干涉生成干涉图。

[0035] 为获得影像覆盖范围内的研究区地面沉降分布情况，首先设定时空基线阈值将预处理获得的干涉图构成干涉对基线网络，然后设定时空相干性阈值筛选相干点目标，以去除相干性较差的点。最后通过相位解缠得到干涉相位时间序列结果，通过网络反演得到研究区内的相应沉降结果信息。

[0036] S2、根据获得的地面沉降空间分布结果将沉降区划分为m个区域，其中，m为沉降漏斗个数，再按照划分的区域导出沉降区的沉降点的基本属性表，并以此建立地面沉降点数据库，将沉降区域进行分类管理。

[0037] 填土区的地面沉降点信息都存储在数据库中，该数据库由多个数据表单组成，各个属性表按照划分的沉降区域分别排列，数据表单主要由形变速率表和累计形变量表两个属性表组成。形变速率表中主要包含沉降观测点的坐标、对应的形变速率等信息；累计形变量表中主要包含沉降观测点的坐标、SAR影像时间序列中各个时间段的累计形变量等信息，此外，每个时间段的沉降观测点都被赋予唯一的ID进行存储。

[0038] S3、在填土沉降区的基础地表形变数据库基础之上，统计出每个沉降区域内沉降速率相同或相近的沉降观测点，将不同区域内所有等值点通过连接形成等值线。

[0039] S4、设定沉降速率和累计沉降量阈值，将大于一定阈值的沉降点视为位于填土区内的沉降点，相应的沉降点构成的等值线即可视为该区域内的填土区边界。实施例

[0040] 本实施例中，我们使用2017年12月6日至2020年12月26日的89景Sentinel‑1A SAR影像数据，选择陕西延安新区作为研究区域。在2012年，延安新区（北区）一期项目正式开工，建设期为4年，北区一期工程位于宝塔区桥沟镇；后续又扩建了高家沟等地，逐渐形成了今天的规模。由于延安新区的大部分房屋和基础设施都是建立在不同时期的地基上，导致地表很容易出现不均匀形变现象，所以开展延安新区的的地表形变研究并对填土区进行识别提取对评估整个延安新区的整体稳定性和安全情况十分重要。

[0041] 具体地，一种面向黄土填土沉降区边界识别与提取方法，包括：

[0042] 获取覆盖研究区域时间序列SAR影像，利用外部精密轨道数据完成轨道校正以去除参考椭球面相位，选定主影像后利用外部数字高程模型(DEM)数据完成地形校正以去除地形相位，再将副影像配准到主影像几何。在时间序列InSAR中需要将所有SAR影像配对组合进行干涉构网。以特定时间顺序（t0、t1、…、tN）排列的覆盖相同区域的N+1景SAR影像，将所有影像一起配准到主影像山应满足以下条件

[0043]                 （1）

[0044] 其中M为干涉对数量；

[0045] 对所有配对的影像进行干涉处理，生成相应的干涉图，再经过相位解缠等步骤获得各沉降观测点的形变速率（见图2）。

[0046] 假设干涉图j由在时间tA和t（B tA

[0047]                 （2）

[0048] 上述方程式也可以表示为：

[0049] （3）

[0050] 式中 和 分别表示像素点(l,r)处tB和tA时刻的相位值，为两幅影像期间的形变相位， 为地形相位， 为大气效应引起
的相位延迟误差， 为随机噪声产生的相位信息。

[0051] 式(3)中各分量分别可表示为：

[0052]     （4）

[0053]  （5）

[0054]  （6）

[0055] λ为雷达波长， 和 表示tA、tB时刻相对于初始时间的LOS向累计形变信息； 表示垂直基线值；为DEM高程差；γ为雷达视线到地面观测物的距离，θ则表示雷达入射角。在去除各种干扰相位后，每个干涉图的相位信息可用以下模型表示：

[0056]    （7）

[0057] 其中， 代表每个干涉图的干涉相位； 是其他影像相对于参考影像的时间干涉相位（假设参考影像的相位为 ），
是残留的干涉相位误差，A是一个M×N矩阵，代表组合的干涉图，由1、0和‑1组成，其中‑1代表主影像，1代表从影像。时间序列的干涉相位的最佳估计是用最小二乘法计算的。

[0058]  （8）

[0059] 其中， 是时间序列干涉相位的最佳估计，W是Fisher Information Matrix（FIM）对角线权重矩阵。方程9表示计算时间相干性（TC）的模型，它可以用来选择具有高相干性的相干点。随后，通过对小型基线网络的加权最小二乘法（WLS）估计，可以得到研究区域内相干点的高精度时间序列变形结果。

[0060] （9）

[0061] 其中，k是虚数单位。

[0062] 通过InSAR技术进行全天时全天候监测可以大大提高监测效率，且可以替代人工进行广域监测，具有极强的适应能力。

[0063] 步骤3、根据沉降漏斗分布将沉降区域划分为SHQ、YZD、JSQ三个区域，然后计算各区域的平均形变速率及累计沉降量。不同用地情况会导致不同沉降类型，有助于区分不同时期填方区并进行不同沉降因素分析。

[0064] 结合图2，从2017年12月到2020年12月，延安新区大部分地区保持相对稳定，变形率主要集中在−10至15 mm/yr。然而，图中有三个区域的地面变形显著，即SHQ、YZD和JSQ区域。所有三个沉降区域均呈带状分布，与延安新区的规划填方区域重叠，最大沉降速率达到−50 mm/yr。SHQ区域有许多建筑工程。YZD区域为交通要道沿线区域。JSQ区域目前正在建造土地。时间序列变形空间分布图还显示，2017年至2020年，延安新区的地表变形范围逐渐减小且速度减缓，这是由于2017年平山造城项目基本完成，沉降恢复工作逐步完成。

[0065] 步骤4、基于划分的三个区域，建立基础地表形变数据库，这样有助于分类统一管理数据。

[0066] 结合图3，基础地面沉降数据库是由三个数据表单组成，数据表单的个数是由划分的区域个数决定。每个数据表单包括两个属性表，分别为沉降速率表和累计沉降量表。

[0067] 步骤5、分析三个沉降区的沉降趋势，并根据各区填方的情况做出不同沉降速率与累计沉降量阈值方案，可以对不同沉降情况区域作出不同预警方案。

[0068] 结合图4，将数值相同或相近的沉降观测点连接为等值线，最后根据设定的阈值将填土沉降区边界识别并提取出来，并结合图5光学历史影像对提取结果进行判别验证。