大气误差改正InSAR干涉图堆叠地质灾害普查方法和装置

摘要

本发明公开了大气误差改正InSAR干涉图堆叠地质灾害普查方法和装置，其中方法包括获取覆盖研究区的统一影像空间的差分干涉图；获取SAR成像时刻监测区域的大气延迟图，并计算得到与差分干涉图所对应的差分大气延迟图；利用各差分干涉图及与其对应的差分大气延迟图进行大气延迟误差改正获得改正后的InSAR差分干涉图；利用大气延迟误差改正后的InSAR差分干涉图利用堆叠法获取监测区平均形变速率；根据监测区平均变形速率进行地质灾害判读。本发明利用大气改正对干涉图进行误差改正后，再进行干涉图堆叠获取形变速率这一组合流程，有效消除大气系统性误差影响的方法的影响，降低大气噪声对InSAR测量的影响，提高观测精度。

说明书

技术领域

本发明涉及合成孔径雷达干涉测量与地质灾害监测技术领域，具体地，涉及一种利用基于大气延迟误差改正的InSAR干涉图堆叠进行地质灾害快速普查的方法。

背景技术

利用合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar,InSAR)技术对地面沉降、滑坡、地裂缝等地质灾害进行调查、监测，具有全天时、全天候、精度高等优点，得到了广泛关注和应用。搭载于卫星上的合成孔径雷达(Synthetic ApertureRadar,SAR)传感器发出的电磁波信号两次穿过大气层，在传播过程中受到大气层，尤其是对流层中大气水汽的影响产生相位延迟，被称为大气效应。大气效应是InSAR测量的主要误差源之一，严重制约着其应用于地质灾害调查和监测时的精度。此外，现有的一些诸如基于永久散射体、分布式散射体或基于小基线集技术的InSAR时间序列分析方法，普遍存在算法复杂、运算效率不高的缺点，且在对InSAR干涉图的大气噪声进行处理时，所利用的滤波方法多基于一定的假设条件。因此在针对大范围地区开展地质灾害调查和监测时，现有方法无法有效消除大气效应影响，观测精度有待提高，且处理过程复杂，耗时较长，难以实现快速普查等目的。

发明内容

本发明旨在有效消除大气效应的影响，提出一种基于大气延迟误差改正的InSAR干涉图堆叠地质灾害普查方法，降低大气噪声对InSAR测量的影响，提高观测精度。

为实现上述技术目的，本发明采用以下技术方案。

一方面，本发明提供大气误差改正InSAR干涉图堆叠地质灾害普查方法，包括以下步骤：

根据获取的拟监测区域的SAR卫星数据获得覆盖研究区的统一影像空间的差分干涉图；

获取SAR卫星影像成像时刻对应的研究区的大气延迟图，并计算得到与差分InSAR数据集所对应的差分大气延迟图；

利用差分干涉图及与其对应的差分大气延迟图进行大气延迟误差改正获得改正后的InSAR差分干涉图；

对利用大气延迟误差改正后的InSAR差分干涉图利用堆叠法获取监测区平均形变速率；

根据监测区平均变形速率进行地质灾害判读。

可选地，获取SAR成像时刻监测区域的大气延迟图并计算得到与差分干涉图所对应的差分大气延迟图包括：获取N幅与SAR卫星影像成像时刻对应的研究区大气延迟图DM，计算得到差分干涉图所对应的差分大气延迟图dDM；

其中确定差分干涉图INT

dDM

可选地，利用差分干涉图及与其对应的差分大气延迟图进行大气延迟误差改正获得改正后的InSAR差分干涉图包括：对应像素的差分干涉图减去差分大气延迟得到改正后的差分干涉图

其中INT

可选地，对改正后的InSAR差分干涉图利用堆叠法获取监测区平均形变速率包括以下步骤：

对于大气延迟误差改正后的差分干涉图

其中v为各像素点在监测时段内的平均形变速率，ΔT

求解观测方程获得平均形变速率。

可选地，求解观测方程获得平均形变速率的方法包括：

利用最小二乘法计算其平均形变速率v，公式如下：

可选地，若无效干涉图数量超过InSAR干涉图总数的设定比例，在最终平均形变速率图上舍弃该点，无效干涉图指的是InSAR差分干涉图上存在点为空值或无效值的情况。

可选地，差分大气延迟图经过空间插值或重采样处理，具有与InSAR干涉图相同的空间分辨率。

可选地，若改正后的InSAR差分干涉图为缠绕相位，则在改正后对InSAR差分干涉图进行相位解缠处理。

第二方面，本发明提供基于大气延迟误差改正InSAR干涉图堆叠的地质灾害普查装置，包括差分干涉图获取模块、差分大气延迟图获取模块、差分大气延迟图获取模块、大气延迟误差改正模块、平均形变速率求解模块以及地质灾害判读模块；

所述差分干涉图获取模块，用于根据获取的拟监测区域的SAR卫星数据获得覆盖研究区的统一影像空间的差分干涉图；

差分大气延迟图获取模块，用于获取SAR卫星影像成像时刻对应的研究区的大气延迟图，并计算得到与差分InSAR数据集所对应的差分大气延迟图；

大气延迟误差改正模块，用于利用差分干涉图及与其对应的差分大气延迟图进行大气延迟误差改正获得改正后的InSAR差分干涉图；

平均形变速率求解模块，用于对利用大气延迟误差改正后的InSAR差分干涉图利用堆叠法获取监测区平均形变速率；

地质灾害判读模块，用于根据监测区平均变形速率进行地质灾害判读。

本发明提供一种计算机可读存储介质,其特征在于,所述计算机存储介质存储有计算机程序,所述计算机程序包括程序指令,所述程序指令当被处理器执行时,使所述处理器执行如权利要求1～8任意一条权利要求所述的基于大气改正InSAR干涉图堆叠的地质灾害普查方法。

本发明所取得的有益技术效果：本发明利用大气延迟误差改正对干涉图进行误差改正后，再进行干涉图堆叠获取形变速率这一组合流程，有效消除大气延迟误差中非随机性误差的影响，降低大气噪声对InSAR测量的影响，提高观测精度。

传统的干涉图堆叠无法克服InSAR大气延迟误差中的非随机性误差，本发明利用大气延迟误差改正先对其进行去除，再利用干涉图堆叠这一简单策略消除残余的随机性大气误差，对获取得到的形变速率结果中的大气延迟误差改正效果非常显著。

附图说明

图1本发明实施例提供的方法流程图；

图2实施例中所形成的InSAR数据集时间、空间基线图；

图3实施例中InSAR数据集的一幅未经大气改正的InSAR差分干涉图示例；

图4实施例中所采用的两幅大气延迟相位图示例；

图5实施例中经计算得到的一幅差分大气延迟相位图示例；

图6实施例中InSAR数据集的一幅经大气改正的InSAR差分干涉图示例；

图7实施例中大气改正InSAR干涉图堆叠获得的地表形变平均速率图；

图8干涉图堆叠法求解平均形变速率时空值/无效点取舍示意图；

图9实施例中未作大气改正InSAR干涉图堆叠获得的地表形变平均速率图。

具体实施方式

下面结合附图对本发明作进一步描述。以下实施例仅用于更清楚地说明本发明的技术方案和实施形式，而不能以此来限制本发明的保护范围。

实施例一、图1是本发明实施例提供的基于大气改正InSAR干涉图堆叠的地质灾害普查方法流程图。具体实施包括五个步骤：

步骤一：根据获取的拟监测区域的SAR卫星数据获得覆盖监测区域的统一影像空间的差分InSAR数据集，所述差分InSAR数据集包括各差分干涉图。

获取拟监测区域的SAR卫星数据，经过影像配准、重采样、平地相位去除、地形相位去除、相位解缠、地理编码等常规InSAR差分干涉数据处理，获取覆盖监测区域的统一影像空间的差分InSAR数据集。InSAR数据集包括，InSAR干涉图数据，为方便阐述，InSAR干涉图记作INT

对于N幅SAR卫星影像，最多可生成N(N-1)/2幅InSAR干涉图。

可选地，经干涉处理所得到的InSAR数据集中干涉图数量不益小于最大干涉组合数N(N-1)/2的30％-40％。所使用的InSAR数据集应尽量均匀地遍布监测时段。InSAR数据集中干涉图的时间基线应有长短搭配。InSAR数据集的影像空间可采用地理坐标系。

本例中，以欧洲空间局哨兵1号(Sentinel-1)卫星获取的SAR卫星影像为数据源，共获取2016年1月至2019年1月间的75幅SAR卫星影像数据，影像大致覆盖范围为东经117度至120度，北纬35度至37度，约270km*210km，成像时刻约为北京时间18:04。

本例中所采用的干涉组合策略为：时间基线小于200天，按照成像时间顺序从前往后进行差分干涉处理，每幅影像最多形成9个干涉图，利用具有30米空间分辨率的SRTM数字高程模型数据进行地形相位去除，数据处理过程中采用30:6多视处理系数，采用最小费用流法进行二维相位解缠，解缠时干涉相干性阈值设置为0.4。本例中共形成616个干涉图构成InSAR数据集，其时间、空间基线情况如附图2所示。

图3举例显示了在进行大气改正前的一幅经过相位解缠，并地理编码至地理坐标系下的差分InSAR干涉图，该干涉图由2017年9月12日(图中记作20170912，下同)和2017年10月6日获取的SAR影像经干涉处理得到。

步骤二：获取SAR成像时刻监测区域的大气延迟图，并计算得到与差分InSAR数据集所对应的差分大气延迟图。可选地，差分大气延迟图应投影变换到对应的卫星视线方向。

本实施例获取75幅SAR数据成像时刻研究区的大气延迟相位图。

图4举例显示了两幅大气延迟相位图，时间为北京时间2017年9月12日18:04和2017年10月6日18:04。

本实施例中，大气相位延迟图由www.gacos.net提供。也可利用地面气象观测资料、全球导航卫星系统(GNSS)观测资料或其他方法得到。

图5显示了由图4中两幅大气延迟图计算得到的差分大气延迟图，并截取得到与图3所示的差分InSAR干涉图同样的空间覆盖范围。

步骤三：利用各差分干涉图及与其对应的差分大气延迟图进行大气延迟误差改正获得改正后的InSAR差分干涉图。

本实施例采用一种通用型InSAR大气延迟误差改正系统数据(GACOS)进行大气延迟误差改正。具体实施例中对每一幅干涉图进行大气改正的方法有很多种，如基于外部气象数据的、基于干涉图本身相位特征的、基于GNSS观测资料的、基于气象模型的等等。因此不限于本实施例采用的方法，可以利用现有技术中的方法进行大气延迟误差改正，本发明不做详细说明。

本实施例中利用各差分干涉图及与其对应的差分大气延迟图进行大气延迟误差改正获得改正后的InSAR差分干涉图包括：对应像素的差分干涉图减去差分大气延迟得到改正后的差分干涉图

其中INT

若进行大气延迟误差改正时的InSAR数据集为缠绕相位，则应在改正后对干涉图进行相位解缠处理。

利用图5所示的差分大气延迟图，对图3的InSAR差分干涉图进行大气延迟误差改正，得到经大气改正的差分干涉图20170912-20171006，如附图6所示。

至此，完成对InSAR差分干涉图20170912-20171006的大气延迟误差改正。对其他615个干涉图进行同样处理，得到经大气改正的InSAR数据集。

步骤四：利用大气延迟误差改正后的InSAR差分干涉图利用堆叠法获取监测区平均形变速率。

利用改正后的InSAR差分干涉图利用堆叠法获取平均形变速率。

对于大气改正后的非缠绕的差分干涉图

其中v为各像素点在监测时段内的平均形变速率，ΔT

可选地，上述公式中平均形变速率v的单位为毫米每年(mm/a)，时间基线ΔT

需要说明的是最小二乘法(Least Square)是求解方程组的通用数学方法，并非堆叠法所特有或唯一方法。凡是数学上能用于求解线性方程组的方法，都可以用于此处。比如最小二乘方法的相关改进方法，加权最小二乘法等，又如最大似然估计方法等。在一定情况下，这些数学方法是相同或可相互转化的。本专利保护的重点不在于某一特定求解方法，最小二乘只是求解的一种路径。

可选地，针对某像素点在若干干涉图上为空值或无效值的情况，若无效干涉图数量超过数据集干涉图总数的设定比例(本实施例设定为1/3)，应在最终平均形变速率图上舍弃该点。

本实施例利用大气改正的InSAR数据集，经干涉图堆叠法获得地表形变平均速率，如图7所示。

图8显示了在利用大气改正InSAR数据集进行干涉图堆叠法计算地表形变速率过程中，某些点为空值或无效值的情况。图中方块所代表的像素在所有干涉图中都存在有效值，三角形代表在多数干涉图上存在有效值，圆形代表该像素仅在少数干涉图上存在有效值。在本例实施过程中，仅保留在多于400个干涉图(占616个干涉图的约2/3)上存在有效值的像素。

作为对比，图9给出未经大气改正的InSAR干涉图堆叠法得到的地表形变平均速率图，可以看出，该结果受大气延迟误差影响严重，易导致对地表变形情况产生误判。

附图7和图9中五角星标志位置为参考点。

步骤五：根据监测区平均变形速率进行地质灾害判读。

根据步骤四得到的监测区平均速率专题图，对地表变形情况进行计算机或专家人工判读，圈定出地质灾害危险区，从而实现地质灾害普查，为进一步的人工详查提供基础资料和指引。

可选地，若前四步处理时InSAR数据集均处于雷达坐标系下，应将干涉图经地理编码处理转换至地理坐标系下。

本实施例根据得到的地表形变平均速率图，进行人工判读，圈定地质灾害危险区，如图7中椭圆圈定范围所示。本实施例中，该区域的地质灾害类型为地表沉降，在进一步的详查中应重点关注该区域情况。

本发明对InSAR差分干涉图首先进行大气改正，有效去除干涉图中的大气延迟误差，为后续处理提供高质量观测数据集；干涉图堆叠算法原理和技术流程相对简单，数据处理速度快，适合地质灾害快速普查推广应用；利用遥感InSAR技术对地质灾害进行大面积普查，首先圈定危险区，可以快速有效缩小地面排查范围，是一种高效手段。

实施例二、基于相同的技术构思，本实施例提供了基于大气改正InSAR干涉图堆叠的地质灾害普查装置，包括：差分干涉图获取模块、差分大气延迟图获取模块、差分大气延迟图获取模块、大气延迟误差改正模块、平均形变速率求解模块以及地质灾害判读模块；

所述差分干涉图获取模块，用于根据获取的拟监测区域的SAR卫星数据获得覆盖的研究区的统一影像空间的差分干涉图；

所述差分大气延迟图获取模块，用于获取SAR成像时刻对应的研究区的的大气延迟图，并计算得到与差分干涉图所对应的差分大气延迟图；

大气延迟误差改正模块，用于利用差分干涉图对及与其对应的差分大气延迟图进行大气延迟误差改正获得改正后的InSAR差分干涉图；

平均形变速率求解模块，用于对利用大气延迟误差改正后的InSAR差分干涉图利用堆叠法获取监测区平均形变速率；

地质灾害判读模块，用于根据监测区平均变形速率进行地质灾害判读。

本实施例中，差分大气延迟图获取模块具体执行以下步骤：获取N幅与SAR卫星影像成像时刻对应的研究区大气延迟图DM，计算得到差分干涉图所对应的差分大气延迟图dDM；

其中确定差分干涉图INT

dDM

上式中DM

本实施例中，大气延迟误差改正模块具体执行以下步骤：对应像素的差分干涉图减去差分大气延迟得到改正后的差分干涉图

其中INT

本实施例中，平均形变速率求解模块具体执行以下步骤：

对于大气延迟误差改正后的差分干涉图

其中v为各像素点在监测时段内的平均形变速率，ΔT

求解观测方程获得平均形变速率。

进一步地，利用最小二乘法计算其平均形变速率v，公式如下：

可选的，若无效干涉图数量超过InSAR干涉图总数的设定比例，在最终平均形变速率图上舍弃该点，无效干涉图指的是InSAR差分干涉图上存在点为空值或无效值的情况。

可选的，差分大气延迟图经过空间插值或重采样处理，具有与InSAR干涉图相同的空间分辨率。

可选的，改正后的InSAR差分干涉图为缠绕相位，则在改正后对InSAR差分干涉图进行相位解缠处理。

所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统，装置、单元或者模块的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述

在以往的研究中，由于卫星影响覆盖范围等因素，大气延迟误差中可能并不显著，但是现在由于卫星影响覆盖范围更广，需要有效消除大气的系统性误差。由于干堆叠虽然是该领域的一种手段，但其本身并不能有效消除InSAR大气的非随机性误差。本发明的核心点在于，先利用大气改正消除系统性误差，再进行干涉图堆叠以消除偶然(随机)大气误差，这种方法可以有效消除大气系统性误差影响的方法。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。