CN111059998B - 一种基于高分辨率的时序InSAR形变监测方法及系统 - Google Patents

一种基于高分辨率的时序InSAR形变监测方法及系统 Download PDF

Info

Publication number
CN111059998B
CN111059998B CN201911416936.1A CN201911416936A CN111059998B CN 111059998 B CN111059998 B CN 111059998B CN 201911416936 A CN201911416936 A CN 201911416936A CN 111059998 B CN111059998 B CN 111059998B
Authority
CN
China
Prior art keywords
image
points
phase
point
processing
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201911416936.1A
Other languages
English (en)
Other versions
CN111059998A (zh
Inventor
杨红磊
刘杰
彭军还
余洋
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
CHINA GEOLOGICAL ENVIRONMENTAL MONITORING INSTITUTE
China University of Geosciences Beijing
Original Assignee
CHINA GEOLOGICAL ENVIRONMENTAL MONITORING INSTITUTE
China University of Geosciences Beijing
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by CHINA GEOLOGICAL ENVIRONMENTAL MONITORING INSTITUTE, China University of Geosciences Beijing filed Critical CHINA GEOLOGICAL ENVIRONMENTAL MONITORING INSTITUTE
Priority to CN201911416936.1A priority Critical patent/CN111059998B/zh
Publication of CN111059998A publication Critical patent/CN111059998A/zh
Application granted granted Critical
Publication of CN111059998B publication Critical patent/CN111059998B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01BMEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B7/00Measuring arrangements characterised by the use of electric or magnetic techniques
    • G01B7/16Measuring arrangements characterised by the use of electric or magnetic techniques for measuring the deformation in a solid, e.g. by resistance strain gauge
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section

Abstract

本发明公开了一种基于高分辨率的时序InSAR形变监测方法及系统。该方法包括:对地表SAR图像进行配准、差分干涉处理和去高程处理;识别图像中的PS点;并对点密度超过设定阈值的PS点进行抽稀处理;对抽稀后的PS点构建三角网;对三角网中相邻的PS点的相位进行差分,通过解空间搜索的方法得到相邻的PS点之间的形变速率差和高程改正量差,并据此计算各PS点的残差相位;从残差相位中分离出大气延迟相位,并对大气延迟相位进行插值处理;从差分干涉图中除去大气延迟相位,重复上述步骤进行迭代,得到多次去除大气延迟相位后的图像,再对图像进行相位解缠以及对所有PS点进行分析,确定地表的形变速率。本发明提高了图像的处理速率。

Description

一种基于高分辨率的时序InSAR形变监测方法及系统
技术领域
本发明涉及微波遥感技术领域,特别是涉及一种基于高分辨率的时序InSAR形变监测方法及系统。
背景技术
IInSAR是一种主动式微波雷达探测技术。它具有大规模,低成本和高精度等技术优势。常规的差分干涉测量(Differential SAR Interferomtry,DInSAR)技术通过对两幅影像作差分来计算形变,但是由于DInSAR技术受失相干现象的影响,目前国际上对InSAR技术的研究主要集中在长时间序列DInSAR技术的研究上,如:PSInSAR技术(Ferretti等,2002),相干目标探测技术(Vander Kooij,2003),干涉点目标分析(Wegmuller,2003;Werner et al.,2003),稳定点网络分析(Arnaud et al.,2003),小基线集技术(Berardinoet al 2003)等。这些技术通过对时间序列中具有稳定、可靠相位的点进行分析,采用迭代的方法计算高精度的形变值和高程改正量,已经成功应用于地表形变、地震形变、冰川移动、火山运动、山地滑坡和大型构建物形变监测中。
而时间序列InSAR技术的难点在于准确分离高程改正和大气延迟量,由于先验的大气改正模型和准确的地表高程值未知,需要从时间序列差分干涉相位通过回归分析多次迭代提取;但随着高分辨率SAR传感器的发展,如TerraSAR-X、COSMO-Skymed和RADARSAT-2等,数据质量、重访周期和分辨率都有明显提高,分辨率达到1m,为精细化形变监测提供了数据保障。分辨率提高的同时,带来的是巨大的数据量;大数据量的迭代严重影响了数据处理的效率,采用普通PC机,甚至图像服务器无法实现对数据快速处理。成为了阻碍InSAR技术向工程化推广的一大问题。
发明内容
本发明的目的是提供一种基于高分辨率的时序InSAR形变监测方法及系统,能够高效的去除图像中的高程改正量和大气延迟量,得到地表的变形速率。
为实现上述目的,本发明提供了如下方案:
一种基于高分辨率的时序InSAR形变监测方法,包括:
获取被监测区域时间序列上的多个地表SAR图像;
对所述地表SAR图像进行配准,并对各地表SAR图像进行差分干涉处理,得到多个差分干涉图;
采用外部DEM数据对所述差分干涉图进行去高程处理,得到第一图像;
对所述第一图像中的PS点进行识别;
对所述第一图像中点密度超过设定阈值的PS点进行抽稀处理;
对抽稀后的PS点构建三角网;
对所述三角网中相邻的PS点的相位进行差分,通过解空间搜索的方法得到相邻的PS点之间的形变速率差和高程改正量差,并根据所述形变速率差和高程改正量差计算得到各PS点的残差相位;
从所述残差相位中分离出大气延迟相位,并对所述大气延迟相位进行插值处理;
从所述第一图像中除去所述大气延迟相位,得到第二图像;以所述第二图像替换所述第一图像,跳转至“对抽稀后的PS点构建三角网”步骤,直至满足设定条件后,停止跳转,得到多次去除大气延迟相位后的图像,记为第三图像;
对各差分干涉图对应的第三图像进行相位解缠;
根据相位解缠后的各第三图像中所有PS点的形变,确定地表的形变速率。
可选的,所述对所述第一图像中的PS点进行识别,具体包括:
采用振幅离差法或相干系数法对所述第一图像中的PS点进行识别。
可选的,在采用所述振幅离差法对所述第一图像中的PS点进行识别时,选取振幅离差小于0.6的点作为PS点。
可选的,所述对所述第一图像中点密度超过设定阈值的PS点进行抽稀处理,具体包括:
确定各所述PS点的点密度;
抽取设定区域内点密度超过设定阈值的PS点中时间序列相干均值相对高的PS点。
可选的,所述设定阈值为50。
可选的,所述从所述残差相位中分离出大气延迟相位,具体包括:
对所述残差相位进行多视处理;
采用空间滤波法对多视处理后的残差相位进行滤波;
对滤波后的残差相位进行插值处理。
本发明还提供了一种基于高分辨率的时序InSAR形变监测系统,包括:
图像获取模块,用于获取被监测区域时间序列上的多个地表SAR图像;
图像预处理模块,用于对所述地表SAR图像进行配准,并对各地表SAR图像进行差分干涉处理,得到多个差分干涉图;
去高程处理模块,用于采用外部DEM数据对所述差分干涉图进行去高程处理,得到第一图像;
PS点识别模块,用于对所述第一图像中的PS点进行识别;
PS点抽稀模块,用于对所述第一图像中点密度超过设定阈值的PS点进行抽稀处理;
三角网构建模块,用于对抽稀后的PS点构建三角网;
相邻点相位差分模块,用于对所述三角网中相邻的PS点的相位进行差分,通过解空间搜索的方法得到相邻的PS点之间的形变速率差和高程改正量差,并根据所述形变速率差和高程改正量差计算得到各PS点的残差相位;
大气延迟相位分离模块,用于从所述残差相位中分离出大气延迟相位,并对所述大气延迟相位进行插值处理;
大气延迟相位去除模块,用于从所述第一图像中除去所述大气延迟相位,得到第二图像;以所述第二图像替换所述第一图像,跳转至所述三角网构建模块,直至满足设定条件后,停止跳转,得到多次去除大气延迟相位后的图像,记为第三图像;
相位解缠模块,用于对各差分干涉图对应的第三图像进行相位解缠;
形变速率确定模块,用于根据相位解缠后的各第三图像中所有PS点的形变,确定地表的形变速率。
可选的,所述PS点识别模块,具体包括:
PS点识别单元,用于采用振幅离差法或相干系数法对所述第一图像中的PS点进行识别。
可选的,所述PS点抽稀模块,具体包括:
点密度计算单元,用于确定各所述PS点的点密度;
PS点抽取单元,用于抽取设定区域内点密度超过设定阈值的PS点中时间序列相干均值相对高的PS点。
可选的,所述大气延迟相位分离模块,具体包括:
多视处理单元,用于对所述残差相位进行多视处理;
空间滤波单元,用于采用空间滤波法对多视处理后的残差相位进行滤波;
插值处理单元,用于对滤波后的残差相位进行插值处理。
根据本发明提供的具体实施例,本发明公开了以下技术效果:本发明提供的基于高分辨率的时序InSAR形变监测方法及系统,对PS点进行抽稀处理,在计算出抽稀后的PS点对应的大气延迟相位后,通过插值的方式获取所有PS点对应的大气延迟相位,降低了数据量,有效的降低了数据迭代花费的时间,便于快速获得大气延迟相位,实现了对高分辨率InSAR数据的快速处理。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例中所需要使用的图作简单地介绍,显而易见地,下面描述中的图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些图获得其他的图。
图1为本发明实施例中基于高分辨率的时序InSAR形变监测方法流程图;
图2a和图2b分别为本发明实施例中PS点位分布图和密度图;
图3为本发明实施例中基于高分辨率的时序InSAR形变监测系统结构图;
图4为本发明示例中的差分干涉图;
图5为本发明示例中原始的PS点图;
图6为本发明示例中抽稀后的PS点图;
图7为本发明示例中去除大气后的差分干涉图;
图8为本发明示例中的地表形变速率图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
本发明的目的是提供一种基于高分辨率的时序InSAR形变监测方法及系统,能够高效的去除图像中的高程改正量和大气延迟量,得到地表的变形速率。
为使本发明的上述目的、特征和优点能够更加明显易懂,下面结合图和具体实施方式对本发明作进一步详细的说明。
如图1所示,本发明的第一方面提供了一种基于高分辨率的时序InSAR形变监测方法包括以下步骤:
步骤101:获取被监测区域时间序列上的多个地表SAR图像;
步骤102:对所述地表SAR图像进行配准,并对各地表SAR图像进行差分干涉处理,得到多个差分干涉图;
步骤103:采用外部DEM数据对所述差分干涉图进行去高程处理,得到第一图像;
步骤104:对所述第一图像中的PS点进行识别;
步骤105:对所述第一图像中点密度超过设定阈值的PS点进行抽稀处理,保留PS点密度没有超过设定阈值的区域中的PS点;
步骤106:对抽稀后的PS点构建三角网;
步骤107:对所述三角网中相邻的PS点的相位进行差分,通过解空间搜索的方法得到相邻的PS点之间的形变速率差和高程改正量差,并根据所述形变速率差和高程改正量差计算得到各PS点的残差相位;
步骤108:从所述残差相位中分离出大气延迟相位,并对所述大气延迟相位进行插值处理;
步骤109:从所述第一图像中除去所述大气延迟相位,得到第二图像;
以所述第二图像替换所述第一图像,跳转至步骤106,直至满足设定条件后,停止跳转,得到多次去除大气延迟相位后的图像,记为第三图像;其中,设定条件可以是满足设定的迭代次数,也可以是直到图像不再受到大气相位影响,即稳定区域相位值接近于0。
步骤110:对各差分干涉图对应的第三图像进行相位解缠;
步骤111:根据相位解缠后的各第三图像中所有PS点的形变,确定地表的形变速率。
在上述实施例中,步骤102中所述的配准可以包括粗配准和精配准,首先,可以以最大化所有干涉像对间的相干性的原则选择参考影像,Zebker和Villasenor(1992)指出,干涉像对间的相干值依赖像对间的时间基线(T)、空间垂直基线(B)、多普勒中心基线(FDC)和热噪声四部分,是四部分的乘积形式,因此从多时相SAR影像中选择主影像时,可以综合考虑这四部分的分布情况,从中选择使得整体最优的一个影像作为主影像。然后通过卫星轨道参数法进行粗配准,再通过相干系数法进行精配准,在精配准的过程中,需要对影像进行重采样,使每幅图的像元精确配准。进而,获取两两幅图间的差分干涉图。
在上述实施例中,在步骤104中,可以采用振幅离差法或相干系数法对所述第一图像中的PS点进行识别。
振幅离差法的核心思想是利用振幅离差和相位标准差统计特性上的关系,当σφ<0.25rad,即在高信噪比(g/σn>4)的条件下,相位标准偏差σφ近似等于振幅离差DA。振幅离差和相位标准差的数学表达式为:
Figure GDA0002413992490000061
式中,σnR和σnl分别为噪声实部和虚部的标准差;μA、σA分别为振幅的均值和标准偏差;σφ和DA分别表示相位标准差和振幅离差。上式表明,即在高信噪比情况下,可直接用振幅离差来表示目标点的稳定性。在本发明中,如果采用振幅离差法识别PS点,识别阈值的优选值可以选取0.6,即可以选取振幅离差大于0.6的点作为PS点。
相干系数法为:设时间序列内有N景SAR影像,配准重采样后形成N-1个干涉对,可以求得N-1个相干值图像,每个象元的相干系数γ可以采用下式计算:
Figure GDA0002413992490000071
式2中,M和S分别表示主从SAR影像的局部象元值,*表示共轭相乘。这样可以得到γ1,γ2,…,γN-1个相干值时间序列。采用式3逐象元计算时间序列相干系数平均值来评价PS点质量的好坏,根据质量的好坏来识别PS点。
Figure GDA0002413992490000072
在上述实施例中,步骤104中对PS点进行抽稀处理可以基于PS点的点密度和质量进行抽稀:
1)计算PS点的点密度:如图2所示,以点(i,j)为中心,设置计算半径r,半径内其它点位坐标为(ii.jj),其中(-r≤ii≤r,-r≤jj≤r),半径内黑点为PS点,标示为1,绿点为非PS点,标示为0,计算该点的点位密度m,点密度等于:以该PS点为中心,半径为r的区域内所有PS点的标示之和,也可以理解为以该PS点为中心,半径为r的区域内所有PS点的数量之和,计算公式如下式4所示:
Figure GDA0002413992490000073
2)判断所述PS点的点密度是否超过所述设定阈值,设定阈值可以优选为50,对于点密度超过所述设定阈值的PS点进行抽稀处理,对其他区域的PS点进行保留。
3)计算PS点的时间序列相干均值。
统计在一定半径内PS点的时间序列相干均值,抽稀时,以此来估计每个点的质量,剔除低质量点,保留高质量点。设置动态抽稀步长,对点密度值和点质量高的区域,采用小步长,反之,采用较大步长。
在上述实施例中,步骤105和步骤106中,对抽稀后的PS点构建delaunay三角网,设x和y点为同一三角网上的两个点,对x和y点做差分,得到差分干涉相位:
Figure GDA0002413992490000081
式中,
Figure GDA0002413992490000082
为相对参考点高程差引起的相位值,
Figure GDA0002413992490000083
是相对参考点视线方向位移差对应的相位值,
Figure GDA0002413992490000084
为相对参考点的大气延迟相位差,
Figure GDA0002413992490000085
为失相关噪声相位。
由于相邻点之间的高程差是垂直基线的函数,相邻点之间的形变速率与时间间隔相关,因此式5也可以表示为:
Figure GDA0002413992490000086
式中,βk为垂直基线,Δh(x,y)为相对高程改正量,Tk为时间间隔,Δv(x,y)为相对形变速度。采用解空间搜索法求解,求解过程中,用到一个重要的参数——整体相位相关系数(ensemble phase coherence),如下式表示:
Figure GDA0002413992490000087
式中,j为虚数单位,
Figure GDA0002413992490000088
为第k幅干涉图中点x和点y之间模型相位和真实相位的差值,如下式表示:
Figure GDA0002413992490000089
Figure GDA00024139924900000810
中的∧用来表示对未知相干值的一个估计值。实际计算中,在二维解空间中,以定长的采样间隔计算
Figure GDA00024139924900000811
使
Figure GDA00024139924900000812
最大的那组解就是两点之间的高程改正量差和形变速度差。
求得相邻点的高程改正量差和形变速度差后,对其积分,可以得到每个永久散射体的绝对速度和高程差,同时完成了每个点的相位解缠,得到每个点的真实差分相位。采用式6求出形变速率差和高程改正量差后,利用式6并通过对弧段积分的形式可计算出残差相位,残差相位中包含了大气延迟和噪声相位。
在上述实施例中,步骤108从所述残差相位中分离出大气延迟相位,并对所述大气延迟相位进行Kriging插值处理,由于大气延迟相位是通过抽稀后的PS点所求取得到的,所以要对大气延迟相位进行Kriging插值处理,以得到所有PS点的大气延迟相位。具体的插值处理如下:
1)计算离散点之间的距离sij和大气延迟相位方差
Figure GDA0002413992490000091
Figure GDA0002413992490000092
Figure GDA0002413992490000093
其中i,j=1,2,…NPS,NPS为PS点总数,xi,yi为PS点坐标,
Figure GDA0002413992490000094
为i点的大气延迟相位;
2)将数据分成距离组(用{h'm}表示),计算各距离组对应的变异函数的估计值γ*(h):
Figure GDA0002413992490000095
Figure GDA0002413992490000096
式中,NH表示距离组的个数。
{h'm}的选择要注意两点:①保证变异函数中有意义的参数,至少要划分3-4组来计算变异函数γ*(h'm),也即是NH≥4;②保证每个距离组包括足够多的数据,以便使γ*(h'm)值更可靠,也即是N(h'm)足够大;
3)根据得到的γ*(h'm)的点的分布形状,选择合适的变异函数进行拟合;
4)求得变异函数估计值γ*(h)后,已知了Kriging方程的系数矩阵K,可以得到未知点的估计值。
在步骤108中,从所述残差相位中分离出大气延迟相位,具体可以采用以下方法:
对所述残差相位进行多视处理;
采用空间滤波法对多视处理后的残差相位进行滤波;
对滤波后的残差相位进行插值处理。
从差分干涉图中去除高程改正量和线性形变后,残差相位中剩余大气延迟相位、非线性形变和噪声。通常情况下,可以通过空间域滤波的方式分离出大气,但空间域滤波效率与滤波窗口成反比,窗口越大,效率越低,不适合对高分辨InSAR数据处理。因此,本发明提出了一种新的大气延迟相位滤波方法:首先对残差相位数据作多视平均,对数据进行多视平均后,虽然降低了空间分辨率,但是去除了噪声的影响。当采用视数较大的窗口多视处理后,非线性形变表现为噪声点,因此需要采用空间滤波方法进行去除。由于空间滤波后的大气延迟相位分辨率和原始图像不一致,因此再采用Kriging插值方法恢复到原始分辨率即可。
在上述实施例中,步骤110和步骤111分别对各差分干涉图对应的第三图像进行相位解缠以及根据相位解缠后的各第三图像中所有PS点的形变,确定地表的形变速率。具体方法如下:
对所有的PS点进行分析
可以表示为:
A·φ=δφ (10)
式中A为M×N矩阵,假设δφ1=φ42,δφ2=φ31,δφM=φNN-2,则上式的矩阵形式如下:
Figure GDA0002413992490000101
A表示干涉图的组合形式。当采用单主影像时,矩阵A的秩为N,采用最小二乘可以求解式11
Figure GDA0002413992490000102
由于本文采用的多主影像组合方式,方程秩亏,即ATA为奇异矩阵,方程存在无数多个解。因此对A进行奇异值分解,获得最小范数意义上的最小二乘解。对A进行SVD分解得到下式:
A=USVT (13)
式13中,S为M×M的对角矩阵,对角线元素为AAt的特征值,U是M×M的正交矩阵,U的列为AAt的正交特征向量,V为N×N矩阵,V的列为AtA的正交特征向量。通常M>N,有M-N个特征值为0,同时由于A秩亏,还有另外L-1个0特征值。总的来讲,S的结构如下式:
S=diag(σ1,…,σN-L+1,0,…,0) (14)
式中,σi为奇异值。
采用奇异值分解可以得到式12最小范数意义上的最小二乘解:
Figure GDA0002413992490000111
式中,S+=diag(1/σ1,…,1/σN-L+1,0,…,0),最后可得:
Figure GDA0002413992490000112
式中,ui和vi分别为U和V的列向量。
采用式10求得相位在时间维上表现不连续,呈现上下跳跃,是一个没有物理意义的解。为了得到一个具有物理意义的解,修改了方程10,采用相邻时间的平均相位速度代替方程中的未知数,新的未知数如下式所示:
Figure GDA0002413992490000113
代入式10。可得:
Figure GDA0002413992490000114
矩阵形式为:
Bv=δφ (19)
式中,B仍是一个M×N矩阵,对第j行,位于主副影像获取时间之间的列,B(i,j)=tk-tk-1,其他情况,B(i,j)=0。对B进行SVD分解可以得到形变速度最小范数意义的最小二乘解,进而获得形变时间序列值。
本发明解决了高分辨率InSAR数据处理效率低的问题,采用基于点密度的数据抽稀方法,对原始PS点集抽稀,依据抽稀的PS点获得大气延迟量和高程改正量,为整个PS点集提供初始值,缩短了用于回归迭代的时间,能有效提高高分辨率InSAR数据处理效率。
本发明的第二方面提供了基于高分辨率的时序InSAR形变监测系统,如图3所示,该系统包括:
图像获取模块301,用于获取被监测区域的地表SAR图像;
图像预处理模块302,用于对所述地表SAR图像进行配准,并对各地表SAR图像进行差分干涉处理,得到多个差分干涉图;
去高程处理模块303,用于采用外部DEM数据对所述差分干涉图进行去高程处理,得到第一图像;
PS点识别模块304,用于对所述第一图像中的PS点进行识别;
PS点抽稀模块305,用于对所述第一图像中点密度超过设定阈值的PS点进行抽稀处理,保留PS点密度没有超过设定阈值的区域中的PS点;
三角网构建模块306,用于对抽稀后的PS点构建三角网;
相邻点相位差分模块307,用于对所述三角网中相邻的PS点的相位进行差分,通过解空间搜索的方法得到相邻的PS点之间的形变速率差和高程改正量差,并根据所述形变速率差和高程改正量差计算得到各PS点的残差相位;
大气延迟相位分离模块308,用于从所述残差相位中分离出大气延迟相位,并对所述大气延迟相位进行插值处理;
大气延迟相位去除模块309,用于从所述第一图像中除去所述大气延迟相位,得到第二图像;以所述第二图像替换所述第一图像,跳转至PS点抽稀模块,直至满足设定条件后,停止跳转,得到多次去除大气延迟相位后的图像,记为第三图像;其中,设定条件可以使满足设定的迭代次数。
相位解缠模块310,用于对各差分干涉图对应的第三图像进行相位解缠;
形变速率确定模块311,用于根据相位解缠后的各第三图像中所有PS点的形变,确定地表的形变速率。
在上述实施例中,所述PS点识别模块304,具体包括:
PS点识别单元,用于采用振幅离差法或相干系数法对所述第一图像中的PS点进行识别。
在上述实施例中,所述大气延迟相位分离模块308,具体包括:
多视处理单元,用于对所述残差相位进行多视处理;
空间滤波单元,用于采用空间滤波法对多视处理后的残差相位进行滤波;
插值处理单元,用于对滤波后的残差相位进行插值处理。
在上述实施例中,PS点抽稀模块305,具体包括:
PS点密度计算单元,用于确定抽稀区域内的PS点密度;
判断单元,用于判断所述PS点密度是否超过所述设定阈值;
PS点抽取单元,用于在所述PS点密度超过所述设定阈值时,根据PS点密度确定PS点的抽取数量,并确定抽稀区域内各PS点的时间序列相干均值,根据所述抽取数量抽取相干均值相对高的PS点。
下面以示例的方式对本发明进行解释说明
如图4~8所示,具体实施步骤如下:
步骤一:数据采用2012-2013年期间郑州市的17景TerraSAR-X数据,采用2013年01月11日的数据作为主影像,进行影像配准,其中一副差分干涉图见图4。
步骤二:采用振幅离差指数阈值方法识别初始PS点,阈值设为0.6,共确定9081850个PS点,见图5。
步骤三:对PS点进行抽稀时,点密度阈值设为50,低于50的点保留,不做抽稀,高于50的点,依照时间序列相干值进行抽稀,抽稀后的PS点个数为521716,占PS点总数的5.74%。抽稀后的PS点分布图如图6所示。
步骤四:对抽稀后的PS点采用解空间搜索法计算出抽稀PS点的线性形变和高程差,并计算得到残差相位。
步骤五:对残差相位进行多视和空间域滤波得到大气相位。
步骤六:对上一步得到的大气相位进行Kriging插值以得到原始分辨率的大气相位。
步骤七:将原始干涉图与大气相位作差分以除去大气,见图7,再次使用抽稀的PS点进行迭代计算以求解大气相位,最终的大气相位为多次迭代得到的大气相位的总和;将多次迭代后得到的最终的大气相位和高程差相位从原始干涉图中减去,并对所有PS点进行分析,得到最终的形变速率,见图8。
从图5中可以看出,PS点主要集中在郑州市城区,在郊区农田也有稀疏的点分布。
从图6中可以看出,抽稀PS点明显降低了高密度区域的PS点数量,同时保留了PS点稀疏区域的点。
从图8中可以看出,最大形变速率为0.121m/a,主要集中在郑州市北部和西部区域,其中北部区域最为严重,具有形变速率大和形变漏斗面积大的特点,而且漏斗显现出连通趋势。
表一统计分析了两种方法各环节的处理时间,时间主要集中在差分干涉处理、时间序列分析和时空域滤波。由表可知,采用本专利提出的方法效率明显优于传统方法,效提高将近一倍。
表一抽稀PS点和不抽稀PS点时间对比
Figure GDA0002413992490000141
表二验证InSAR监测成果。在监测区部署了水准监测,采用InSAR、水准同步观测的方式,利用水准观测值对InSAR精度进行分析与评价,从表中可看出本发明提出方法是可靠的,在提高数据处理效率的同时能够实现mm级的形变监测精度。
表二精度评定表
观测时期 平均误差(mm) 中误差(mm)
2012.11-2013.6 ±3.8 ±4.6
2012.11-2013.8 ±3.7 ±4.5
2013.6-2013.8 ±1.5 ±1.9
年平均 ±3.5 ±4.4
本发明提供的基于高分辨率的时序InSAR形变监测方法及系统,对高点位密度区域进行抽稀,同时保留低点位密度区域的PS点,然后通过插值的方式获取所有PS点的大气延迟和高程改正,降低了数据量,有效的降低了数据迭代花费的时间,便于快速获得大气延迟相位和高程改正量,实现了对高分辨率InSAR数据的快速处理。
本说明书中各个实施例采用递进的方式描述,每个实施例重点说明的都是与其他实施例的不同之处,各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言,由于其与实施例公开的方法相对应,所以描述的比较简单,相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述,以上实施例的说明只是用于帮助理解本发明的方法及其核心思想;同时,对于本领域的一般技术人员,依据本发明的思想,在具体实施方式及应用范围上均会有改变之处。综上所述,本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种基于高分辨率的时序InSAR形变监测方法,其特征在于,包括:
获取被监测区域时间序列上的多个地表SAR图像;
对所述地表SAR图像进行配准,并对各地表SAR图像进行差分干涉处理,得到多个差分干涉图;
采用外部DEM数据对所述差分干涉图进行去高程处理,得到多个第一图像;
对所述第一图像中的PS点进行识别;
对所述第一图像中点密度超过设定阈值的PS点进行抽稀处理;
对抽稀后的PS点构建三角网;
对所述三角网中相邻的PS点的相位进行差分,通过解空间搜索的方法得到相邻的PS点之间的形变速率差和高程改正量差,并根据所述形变速率差和高程改正量差计算得到各PS点的残差相位;
从所述残差相位中分离出大气延迟相位,并对所述大气延迟相位进行插值处理;
从所述第一图像中除去所述大气延迟相位,得到第二图像;以所述第二图像替换所述第一图像,跳转至“对抽稀后的PS点构建三角网”步骤,直至满足设定条件后,停止跳转,得到多次去除大气延迟相位后的图像,记为第三图像;
对各差分干涉图对应的第三图像进行相位解缠;
根据相位解缠后的各第三图像中所有PS点的形变,确定地表的形变速率。
2.根据权利要求1所述的基于高分辨率的时序InSAR形变监测方法,其特征在于,所述对所述第一图像中的PS点进行识别,具体包括:
采用振幅离差法或相干系数法对所述第一图像中的PS点进行识别。
3.根据权利要求2所述的基于高分辨率的时序InSAR形变监测方法,其特征在于,在采用所述振幅离差法对所述第一图像中的PS点进行识别时,选取振幅离差小于0.6的点作为PS点。
4.根据权利要求1所述的基于高分辨率的时序InSAR形变监测方法,其特征在于,所述对所述第一图像中点密度超过设定阈值的PS点进行抽稀处理,具体包括:
确定各所述PS点的点密度;
抽取设定区域内点密度超过设定阈值的PS点中时间序列相干均值相对高的PS点。
5.根据权利要求1或4所述的基于高分辨率的时序InSAR形变监测方法,其特征在于,所述设定阈值为50。
6.根据权利要求1所述的基于高分辨率的时序InSAR形变监测方法,其特征在于,所述从所述残差相位中分离出大气延迟相位,具体包括:
对所述残差相位进行多视处理;
采用空间滤波法对多视处理后的残差相位进行滤波;
对滤波后的残差相位进行插值处理。
7.一种基于高分辨率的时序InSAR形变监测系统,其特征在于,包括:
图像获取模块,用于获取被监测区域时间序列上的多个地表SAR图像;
图像预处理模块,用于对所述地表SAR图像进行配准,并对各地表SAR图像进行差分干涉处理,得到多个差分干涉图;
去高程处理模块,用于采用外部DEM数据对所述差分干涉图进行去高程处理,得到多个第一图像;
PS点识别模块,用于对所述第一图像中的PS点进行识别;
PS点抽稀模块,用于对所述第一图像中点密度超过设定阈值的PS点进行抽稀处理;
三角网构建模块,用于对抽稀后的PS点构建三角网;
相邻点相位差分模块,用于对所述三角网中相邻的PS点的相位进行差分,通过解空间搜索的方法得到相邻的PS点之间的形变速率差和高程改正量差,并根据所述形变速率差和高程改正量差计算得到各PS点的残差相位;
大气延迟相位分离模块,用于从所述残差相位中分离出大气延迟相位,并对所述大气延迟相位进行插值处理;
大气延迟相位去除模块,用于从所述第一图像中除去所述大气延迟相位,得到第二图像;以所述第二图像替换所述第一图像,跳转至所述三角网构建模块,直至满足设定条件后,停止跳转,得到多次去除大气延迟相位后的图像,记为第三图像;
相位解缠模块,用于对各差分干涉图对应的第三图像进行相位解缠;
形变速率确定模块,用于根据相位解缠后的各第三图像中所有PS点的形变,确定地表的形变速率。
8.根据权利要求7所述的基于高分辨率的时序InSAR形变监测系统,其特征在于,所述PS点识别模块,具体包括:
PS点识别单元,用于采用振幅离差法或相干系数法对所述第一图像中的PS点进行识别。
9.根据权利要求7所述的基于高分辨率的时序InSAR形变监测系统,其特征在于,所述PS点抽稀模块,具体包括:
点密度计算单元,用于确定各所述PS点的点密度;
PS点抽取单元,用于抽取设定区域内点密度超过设定阈值的PS点中时间序列相干均值相对高的PS点。
10.根据权利要求7所述的基于高分辨率的时序InSAR形变监测系统,其特征在于,所述大气延迟相位分离模块,具体包括:
多视处理单元,用于对所述残差相位进行多视处理;
空间滤波单元,用于采用空间滤波法对多视处理后的残差相位进行滤波;
插值处理单元,用于对滤波后的残差相位进行插值处理。
CN201911416936.1A 2019-12-31 2019-12-31 一种基于高分辨率的时序InSAR形变监测方法及系统 Active CN111059998B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911416936.1A CN111059998B (zh) 2019-12-31 2019-12-31 一种基于高分辨率的时序InSAR形变监测方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911416936.1A CN111059998B (zh) 2019-12-31 2019-12-31 一种基于高分辨率的时序InSAR形变监测方法及系统

Publications (2)

Publication Number Publication Date
CN111059998A CN111059998A (zh) 2020-04-24
CN111059998B true CN111059998B (zh) 2020-11-13

Family

ID=70305777

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911416936.1A Active CN111059998B (zh) 2019-12-31 2019-12-31 一种基于高分辨率的时序InSAR形变监测方法及系统

Country Status (1)

Country Link
CN (1) CN111059998B (zh)

Families Citing this family (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111562575B (zh) * 2020-06-01 2022-12-09 江苏中煤地质工程研究院有限公司 一种用于地面沉降的监测方法
CN111895903B (zh) * 2020-07-21 2021-06-01 太原理工大学 一种检测区域积雪深度的遥感估算方法
CN111999704B (zh) * 2020-08-18 2023-05-26 中国电子科技集团公司第三十八研究所 一种基于vpx总线的车载雷达时序产生系统及方法
CN112596055B (zh) * 2020-12-08 2023-04-25 中国地质大学(武汉) 一种改正InSAR DEM残余系统误差的方法
CN112540370A (zh) * 2020-12-08 2021-03-23 鞍钢集团矿业有限公司 InSAR和GNSS融合的露天矿边坡形变测量方法
CN112946647A (zh) * 2021-02-02 2021-06-11 河海大学 大气误差改正InSAR干涉图堆叠地质灾害普查方法和装置
CN112986990B (zh) * 2021-02-04 2023-02-17 中国地质大学(北京) 一种大气相位改正方法及系统
BR102021005841A2 (pt) * 2021-03-25 2022-09-27 Radaz Industria E Comercio De Produtos Eletronicos Ltda Método de localização de alterações no subsolo
CN113239136B (zh) * 2021-05-17 2023-12-19 北京车和家信息技术有限公司 一种数据处理方法、装置、设备及介质
CN113513978B (zh) * 2021-06-02 2023-04-14 北京卫星制造厂有限公司 高低温环境下的端面位姿相对变化高精度测量方法和系统
CN113253234B (zh) * 2021-06-17 2021-10-29 中国人民解放军国防科技大学 目标微小形变观测雷达系统的信号处理方法及雷达系统
CN114594479B (zh) * 2022-05-07 2022-07-26 中国测绘科学研究院 一种全散射体FS-InSAR方法及系统
CN115201825B (zh) * 2022-09-16 2023-01-17 眉山环天智慧科技有限公司 一种InSAR震间形变监测中的大气延迟校正方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104006791A (zh) * 2014-05-06 2014-08-27 国家基础地理信息中心 基于多源遥感影像的城市地区高程维度变化信息提取方法
CN104091064A (zh) * 2014-07-02 2014-10-08 北京航空航天大学 基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法
CN104111456A (zh) * 2014-07-23 2014-10-22 中国国土资源航空物探遥感中心 一种高速铁路沿线地表形变高分辨率InSAR监测方法

Family Cites Families (19)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102608584B (zh) * 2012-03-19 2014-04-16 中国测绘科学研究院 基于多项式反演模型的时间序列InSAR形变监测方法及装置
CN103675790B (zh) * 2013-12-23 2016-01-20 中国国土资源航空物探遥感中心 一种基于高精度DEM提高InSAR技术监测地表形变精度的方法
CN104111457B (zh) * 2014-07-23 2016-10-12 中国国土资源航空物探遥感中心 一种升降轨PSInSAR地面沉降监测结果的互检验与时序融合方法
CN104407332B (zh) * 2014-11-25 2017-08-11 沈阳建筑大学 一种地基sar更新dem的校正方法
CN104459696B (zh) * 2014-12-24 2017-02-22 中南大学 一种基于平地相位的sar干涉基线精确估计方法
CN105954747A (zh) * 2016-06-20 2016-09-21 中国电力工程顾问集团中南电力设计院有限公司 基于电网不良地质体三维形变监测的塔基稳定性分析方法
CN106226764A (zh) * 2016-07-29 2016-12-14 安徽理工大学 一种基于D‑InSAR的煤矿开采地沉陷区域的测定方法
CN106556834B (zh) * 2016-11-24 2018-11-16 首都师范大学 一种从两平行轨道sar数据集中精确提取地面垂直形变方法
CN106772377A (zh) * 2017-01-18 2017-05-31 深圳市路桥建设集团有限公司 一种基于InSAR的建筑物变形监测方法
CN107037428B (zh) * 2017-03-27 2019-11-12 中国科学院遥感与数字地球研究所 一种提高星载双站差分InSAR提取形变精度的方法
CN108007401A (zh) * 2017-11-20 2018-05-08 贵州省水利水电勘测设计研究院 一种基于船载InSAR平台的河湖库沿岸形变检测装置及方法
WO2019126972A1 (zh) * 2017-12-26 2019-07-04 深圳市城市公共安全技术研究院有限公司 基于InSAR的形变信息提取方法、终端及存储介质
CN108710984B (zh) * 2018-04-04 2021-06-08 中国地质环境监测院 一种矿山地质环境综合评价方法及系统
CN108387899B (zh) * 2018-04-17 2020-07-31 南京师范大学 合成孔径雷达干涉测量中地面控制点自动选取方法
CN108663017A (zh) * 2018-08-13 2018-10-16 伟志股份公司 一种监测城市地铁沿线地表沉降的方法
CN109031301A (zh) * 2018-09-26 2018-12-18 云南电网有限责任公司电力科学研究院 基于PSInSAR技术的山区地形形变提取方法
CN109886910B (zh) * 2019-02-25 2021-07-13 榆林市国土资源局 外部数字高程模型dem修正方法及装置
CN110174044B (zh) * 2019-04-16 2021-08-03 中国矿业大学 一种基于psi技术的桥梁纵向位移形变监测的方法
CN110412574B (zh) * 2019-09-05 2021-05-25 河海大学 一种时空相干性增强的分布式目标InSAR时序处理方法和装置

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104006791A (zh) * 2014-05-06 2014-08-27 国家基础地理信息中心 基于多源遥感影像的城市地区高程维度变化信息提取方法
CN104091064A (zh) * 2014-07-02 2014-10-08 北京航空航天大学 基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法
CN104111456A (zh) * 2014-07-23 2014-10-22 中国国土资源航空物探遥感中心 一种高速铁路沿线地表形变高分辨率InSAR监测方法

Also Published As

Publication number Publication date
CN111059998A (zh) 2020-04-24

Similar Documents

Publication Publication Date Title
CN111059998B (zh) 一种基于高分辨率的时序InSAR形变监测方法及系统
Gardner et al. Increased West Antarctic and unchanged East Antarctic ice discharge over the last 7 years
CN108010320B (zh) 一种基于自适应时空约束低秩算法的路网交通数据的补全方法
CN110673145B (zh) 一种基于间断相干的InSAR地表形变监测方法及系统
Liu et al. Estimating Spatiotemporal Ground Deformation With Improved Persistent-Scatterer Radar Interferometry $^\ast$
CN110763187B (zh) 一种稳健的基于雷达分布式目标的地面沉降监测方法
CN111856459B (zh) 一种改进的DEM最大似然约束多基线InSAR相位解缠方法
CN111398959B (zh) 基于地表应力应变模型的InSAR时序地表形变监测方法
CN110532953B (zh) 基于纹理特征辅助的sar影像冰川识别方法
CN111273293A (zh) 一种顾及地形起伏的InSAR残余运动误差估计方法及装置
CN116047519B (zh) 一种基于合成孔径雷达干涉测量技术的选点方法
CN112051572A (zh) 一种融合多源sar数据三维地表形变监测方法
CN108983231B (zh) 一种基于视频合成孔径雷达的干涉视频测量方法
CN116381680A (zh) 一种基于时序雷达干涉测量技术的城市地表形变监测方法
CN115079172A (zh) 一种MTInSAR滑坡监测方法、设备及存储介质
CN113238228B (zh) 基于水准约束的三维地表形变获取方法、系统及装置
Warner et al. Pine Island Glacier (Antarctica) velocities from Landsat7 images between 2001 and 2011: FFT-based image correlation for images with data gaps
CN113552565B (zh) 针对sar数据高噪声及大梯度变化区域的相位解缠方法
Chen et al. Surface velocity estimations of ice shelves in the northern Antarctic Peninsula derived from MODIS data
CN114812491A (zh) 基于长时间序列分析的输电线路地表形变预警方法及装置
CN110297243B (zh) 合成孔径雷达层析三维成像中的相位误差补偿方法及装置
CN113204023A (zh) 联合ps目标与ds目标的双极化相位优化地表形变监测方法
CN111915570A (zh) 基于反向传播神经网络的大气延迟估计方法
CN113866765B (zh) 基于多成分时间相干模型的PS-InSAR测量方法
CN115546264A (zh) 一种星载InSAR图像精细配准与立体测量方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant