CN109884635A - 大范围高精度的InSAR形变监测数据处理方法 - Google Patents

Abstract

本发明提供了一种大范围高精度的InSAR形变监测数据处理方法，包括：对多幅SAR影像进行辐射定标，通过选取不同程度的相干系数阈值和振幅离差阈值，选出不同质量等级的PS点；将所选PS点的像素位置和坐标信息归入不同的规则格网内，相邻格网间设置第一预设范围的重叠区域；每个格网内的解算高等级PS点构成该格网内PS点的候选参考点集，通过某一候选点与其余候选点构网解算求解判断残差大小的形式，选出该格网内的最优参考点。本发明能够有效解决海量数据下PS‑InSAR用于大范围形变监测中的误差传递导致解算结果精度较低的问题，以及数据处理和存储过程中对计算机内存需求过大易导致的内存溢出解算失败问题。

Description

大范围高精度的InSAR形变监测数据处理方法

技术领域

本发明涉及基于遥感影像的大地测量领域，特别涉及一种大范围高精度的InSAR形变监测数据处理方法。

背景技术

合成孔径干涉雷达测量(InSAR-Interferometric Synthetic Aperture radar)具有高精度、大范围、全天时全天候的优点，在地表形变监测中得到极大的应用和发展。基于其发展出来的永久散射体技术(PS-InSAR)，已经被成功应用于各类形变监测中，包括滑坡、城市沉降、基础设施沉降、矿区沉陷等，并且得到了精度达到毫米级的形变监测结果。得益于近些年各类SAR卫星的发射尤其是TerraSAR-X、COSMO-SkyMed、Sentinel等，我们可获取使用的影像数据变得越来越丰富，空间分辨率也越来越高，这极大地提高了PSInSAR技术监测地表形变的能力。

然而，面对海量SAR数据以及大范围精细形变监测的需求，PSInSAR技术面临着新的挑战。PS-InSAR技术用于大范围区域形变监测目前主要有以下缺陷：1)大范围形变监测过程中的误差传递。PSInSAR通过筛选符合要求的PS点，通过一定的方式在相邻点间构建弧段组成一定三角网，通过积分解算。在解算范围较小的情况下(小于5km*5km)，能够保证形变结果的可靠性。但是，面对大范围的形变监测(上百公里)，三角网积分解算在缺少丰富地面控制点的情况下，会出现误差的传递，最终影响形变解算结果的精度和可靠性。

2)大范围海量SAR数据解算对电脑硬件提出的过大需求。随着SAR影像时间空间分辨率的不断提高，对地物的辨析能力越来越强，因此，大范围精细监测意味着PS点和弧段成倍的增长。这就导致PS-InSAR解算过程中，积分运算需要占据异常巨大的计算机内存，同时，结算过程十分耗时，甚至以天计算。例如，通过对834056个PS点进行构网解算，会构建2501630个弧段，积分解算时，单设计矩阵就需要消耗内存15545.6G，这在较高配置的计算机上也难以满足。即使有计算机能够满足数据运算和存储的要求，解算过程需要的时间也难以想象。这些都与现如今简单高效的应用需求相违背。

3)针对PS-InSAR应用于大范围形变监测的问题，有人尝试通过分块解算的方法来求解。但是，简单的分块解算缺少整体的误差约束，导致解算误差、大气延迟等误差在分块间的传递，出现部分分块出现整体性的偏差，即所谓的分块现象，使得解算结果缺乏连续性，影响最终结果的精度。

根据以上对现有的PS-InSAR大范围形变监测方法的分析，可以看出现有方法受到算法和计算机硬件的局限，在现有条件和需求下难以实施，同时受误差传递的影响，PS-InSAR干涉测量技术无法达到其应有的理论精度，从而影响了该技术的应用和推广。

发明内容

本发明提供了一种大范围高精度的InSAR形变监测数据处理方法，其目的是为了解决海量数据下PS-InSAR用于大范围形变监测中的误差传递导致解算结果精度较低的问题，以及数据处理和存储过程中对计算机内存需求过大易导致的内存溢出解算失败的问题。

为了达到上述目的，本发明的实施例提供了一种大范围高精度的InSAR形变监测数据处理方法，包括：

步骤1，对多幅SAR影像进行辐射定标，通过选取不同程度的相干系数阈值和振幅离差阈值，选出不同质量等级的PS点；

步骤2，将所选PS点的像素位置和坐标信息归入不同的规则格网内，相邻格网间设置第一预设范围的重叠区域；

步骤3，每个格网内的解算高等级PS点构成该格网内PS点的候选参考点集，通过某一候选点与其余候选点构网解算求解判断残差大小的形式，选出该格网内的最优参考点，与该格网内的其余PS点连接构建星型网络，形成PS-InSAR解算的二级网络，每个格网内的最优参考点之间采用Delaunay三角网形式构建PS-InSAR解算的一级网络；

步骤4，在每个格网内，假定参考点的形变速率和高程残差为0，在格网内星型网间，采用解空间搜索法估算出每个PS点的高程残差和线性形变速率；

步骤5，通过利用相邻格网重叠区域PS点解算结果的差异作为观测值，以真实起算点所在格网的重叠区域PS解算结果作为起算数据，建立覆盖所有格网的误差方程和法方程，通过网平差，校准格网PS点相对高程残差和线性形变速率。

其中，所述步骤2具体包括：

将所选PS点的像素位置和坐标信息将其归入不同的规则格网时，大气延迟在不足1km的两点大气相位差的方差不足0.1rad²，将研究区域划分成1km*1km的格网，并且相邻格网间设定20％的重叠，确保足够数量的PS点落在重叠区域。

其中，所述步骤3具体包括：

根据阈值选取格网内的高等级PS点，剔除掉落在重叠区域的点，组成该格网内的候选PS参考点集；

分别以每一个候选参考点为中心，在候选参考点集中构建星型解算网，利用解空间搜索法估算出每一个弧段的高程残差之差Δε_h和线性形变速率之差Δν，并计算该弧段参数估计后平均残余相位

其中，i表示第i幅干涉图，t_i和分别表示该幅干涉图的时间基线和空间垂直基线；λ、R、θ分别表示波长、斜距以及入射角；N表示干涉对的数量，ΔΦ(x,y,i)表示弧段两端PS点差分干涉相位之差；

计算每一个星型网中所有弧段的残差均值，残差均值最小的星型网对应的中心候选参考点作为该格网的参考点，真实起算数据点所在的格网，以真实起算数据点作为该格网的参考点。

其中，所述步骤5具体包括：

计算每一个重叠区域PS点解算结果之差，统计其分布情况，取其主要分布聚集区间的中值作为该重叠区域的差异值，即相邻格网的计算结果差异值；

以计算结果差异值作为观测值，以每个格网PS点解算结果整体改正量为代求参数，建立覆盖所有格网的平差网络，设定真实起算点所在格网的改正量为零，并以此为起算数据，进行最小二乘平差求解每个格网的校准改正量；

根据校准改正量对每一个格网估算的高程残差和香型形变速率进行整体平移校正，便得到整个研究区域最终连续的解算结果。

其中，所述求解每个格网的校准改正量的步骤具体包括：

假定以部分重叠的方式将研究区域划分为M个格网，形成N个两两重叠的区域(N>M)，并附加真实起算点所在格网的改正量为零的限定，则有条件方程：

B·x-L＝0

其中B为参数x的系数矩阵，且有：

L＝[L₁,L₂,…,L_N,0]^T，x＝[x₁,x₂,…,x_M]T

根据最小二乘平差原理，可以得到以下法方程：

B^TPBx＝B^TPL

其中，P为观测向量L的权矩阵，在此取单位权；

通过

计算出每个格网结算值的校准改正量。

本发明的上述方案有如下的有益效果：

本发明的上述实施例所述的大范围高精度的InSAR形变监测数据处理方法通过充分利用格网重叠区域的信息，避免了分块解算过程中出现的监测结果不连续的现象，并且在格网解算过程中通过并行运算设计，大幅度地降低了解算所需的时间；可以实现PS-InSAR技术大范围、高效率、高精度地表形变监测，方法易于实施，应用广泛，有助于扩展InSAR干涉测量技术的应用范围。

附图说明

图1为本发明的流程示意图；

图2为本发明的整体流程示意图。

具体实施方式

为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述。

本发明针对现有的海量数据下PS-InSAR用于大范围形变监测中的误差传递导致解算结果精度较低的问题，以及数据处理和存储过程中对计算机内存需求过大易导致的内存溢出解算失败的问题，提供了一种大范围高精度的InSAR形变监测数据处理方法。

如图1所示，本发明的实施例提供了一种大范围高精度的InSAR形变监测数据处理方法，具体方案如下：

1)前期数据处理：

根据研究区域的N景(一般>25景)SAR影像，充分考虑各影像件的有效空间基线、时间基线和多普勒质心频率差的基础上，通过计算整体相关系数，选取最大系数对应的影像为主影像，对其余影像进行配准和重采样。随后通过主影像和其余从影像两两组队进行干涉处理，并引入外部DEM来去除干涉相位中的地形相位，得到一系列覆盖研究区域的差分干涉图。第i幅干涉图中的某点x的差分干涉相位可以表示为：

Φ(x,t_i)＝Φ_defo+Φ_{topo_res}+Φ_atm+Φ_noise (1)

其中，Φ(x,t_i)表示第i幅干涉图中x点处的差分干涉相位，该相位由形变相位Φ_defo、地形残余相位Φ_{topo_res}、大气延迟相位Φ_atm以及噪声相位Φ_noise组成。

由于PS-InSAR技术主要用于监测目标的长期缓慢形变，因此我们可以将形变相位分解为一个随时间均匀变化的线性分量，以及实际形变量与线性形变量的差异部分，该差异部分也即非线性分量。则有

Φ_defo＝Φ_{defo_linear}+Φ_{defo_unlinear} (2)

同时，高程残余相位被证实和干涉对的空间基线存在一定的关系。因此，式(1)可以改写为：

其中，t_i和分别表示该幅干涉图的时间基线和空间垂直基线；λ、R、θ分别表示波长、雷达到点x的斜距以及改点处的入射角；ν(x)、分别表示PS点x处的形变速率以及高程改正值；表示残余相位，主要包含非线性形变相位大气延迟相位和噪声相位

另外，该部分也可以获得干涉对的相干性信息和SAR影像的强度信息，可以用于后续PS点的筛选工作。

2)PS点的选取与分块格网的划分：

首先，对获取的N景SAR影像强度信息进行辐射校正，步骤如下：

步骤一：获取SAR影像的幅度均值I_i；

其中，为第i幅SAR影像的幅度值，m，n为该幅影像方位向和距离向的采样数。

步骤二：计算SAR影像序列的幅度均值

步骤三：对每幅SAR影像进行辐射校正，获得校正后的幅度值

完成辐射校正后，根据上述获取的干涉相干性序列和校正后的SAR幅度值序列，通过振幅离差阈值法和平均相干性阈值法，设定不同程度的阈值，筛选出不同质量等级的永久散射体点。

选取PS点后，根据所选PS点的分布和像素坐标，将其划分入不同的格网内。格网均等划分，划分格网大小时应充分考虑大气延迟在空间域的特性，即大气延迟在一定空间范围内具有较强的相关性，不足1km的两点大气相位差的方差不足0.1rad²。因此，该方法将研究区域划分成1km*1km的格网，并且相邻格网间设定20％的重叠，确保足够数量的PS点落在重叠区域内。

3)格网内参考点的选取与格网解算

为实现分块间的独立并行计算，需要为每个分块选出起算参考点。由于各格网内采用星型方式构网，因此参考点需要保证较高的相位质量。可以通过以下步骤实现格网参考点的选取与解算：

步骤一：设定较高的相干性阈值和较低的振幅离差阈值，对每个格网内的PS点进行筛选，剔除掉落在重叠区域的点后，获得在相干性和振幅离差两方面均表现优异的点序列，作为该格网内的候选参考点；

步骤二：分别以每一个候选参考点为中心，在候选参考点集中构建星型解算网，利用解空间搜索法估算出每一个弧段的高程残差之差Δε_h和线性形变速率之差Δν，并计算该弧段参数估计后平均残余相位

其中，i表示第i幅干涉图，t_i和分别表示该幅干涉图的时间基线和空间垂直基线；λ、R、θ分别表示波长、斜距以及入射角；N表示干涉对的数量，ΔΦ(x,y,i)表示弧段两端PS点差分干涉相位之差；

步骤三：计算每一个星型网所有弧段的残差均值，残差均值最小的星型网对应的中心候选参考点作为该格网的参考点。真实起算数据点所在的格网，以真实起算数据点作为该格网的参考点；

步骤四：在筛选出每个格网的参考点后，针对每个格网，分别连接参考点和格网内其余的每一个PS点组成弧段，建立以参考点为中心的星型网。对弧段两边的x和y两个PS点的差分干涉相位Φ(x,t_i)和Φ(y,t_i)再次进行差分，其差分相位如下：

其中，参数Δν、Δε_h分别代表x、y两点间的线性形变速率之差和高程改正之差，为x、y两点的残余相位之差；利用解空间搜索法，在一定的解空间合理范围内，通过迭代搜索出参数的最优解。

步骤五：假定格网参考点的线性形变速率和高程残差为零，并以此作为起算数据，利用步骤四获得的各弧段参数，求解格网内各PS点的线性形变速率值和高程改正值。值得注意的是，由于参考点的线性形变速率和高程残差起算值是假定值，因此这一步格网PS点解算结果为相对值，在后续步骤中会对此进行校正。

4)参考点解算

通过以上，实现了PS-InSAR的分块解算，并且各分块内的解算过程相互独立，可以实现并行运算。但是上述结果为相对各分块内参考点的相对量。为实现相对结果的绝对纠正，需要对各参考点的零值假设进行绝对校正。

设定每个格网间有一定程度的重叠，重叠区域的PS点在不同格网解算中会有不同的解算结果。本方法主要利用相邻格网重叠区域PS点解算结果的差异校准格网PS点相对高程残差和相对线性形变速率。具体步骤如下：

步骤一：计算每一个重叠区域PS点解算结果之差，统计其分布情况，取其主要分布聚集区间的中值作为该重叠区域的差异值，也即相邻格网计算结果的差异值；

步骤二：以计算结果差异值作为观测值，以每个格网PS点解算结果整体改正量为代求参数，建立覆盖所有格网的平差网络。设定真实起算点所在格网的改正量为零，并以此为起算数据，可以进行最小二乘平差求解每个格网的校准改正量。具体地：

假定以部分重叠的方式将研究区域划分为M个格网，形成N个两两重叠的区域(N>M)，并附加真实起算点所在格网的改正量为零的限定，则有条件方程：

B·x-L＝0

其中B为参数x的系数矩阵，且有：

L＝[L₁,L₂,…,L_N,0]^T，x＝[x₁,x₂,…,x_M]T

根据最小二乘平差原理，可以得到以下法方程：

B^TPBx＝B^TPL

其中，P为观测向量L的权矩阵，在此取单位权。

从而可求得基线参数向量x最优线性无偏估计

自此求出各格网结算值的改正量。

步骤三：利用步骤二求解出的改正量对每一个格网估算的高程残差和香型形变速率进行整体平移校正，便得到整个研究区域最终连续的解算结果。

本发明的上述实施例所述的大范围高精度的InSAR形变监测数据处理方法通过充分利用格网重叠区域的信息，避免了分块解算过程中出现的监测结果不连续的现象，并且在格网解算过程中通过并行运算设计，大幅度地降低了解算所需的时间；可以实现PS-InSAR技术大范围、高效率、高精度地表形变监测，方法易于实施，应用广泛，有助于扩展InSAR干涉测量技术的应用范围。

以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明所述原理的前提下，还可以作出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (5)

1.一种大范围高精度的InSAR形变监测数据处理方法，其特征在于，包括：

步骤1，对多幅SAR影像进行辐射定标，通过选取不同程度的相干系数阈值和振幅离差阈值，选出不同质量等级的PS点；

步骤2，将所选PS点的像素位置和坐标信息归入不同的规则格网内，相邻格网间设置第一预设范围的重叠区域；

步骤3，每个格网内的解算高等级PS点构成该格网内PS点的候选参考点集，通过某一候选点与其余候选点构网解算求解判断残差大小的形式，选出该格网内的最优参考点，与该格网内的其余PS点连接构建星型网络，形成PS-InSAR解算的二级网络，每个格网内的最优参考点之间采用Delaunay三角网形式构建PS-InSAR解算的一级网络；

步骤4，在每个格网内，假定参考点的形变速率和高程残差为0，在格网内星型网间，采用解空间搜索法估算出每个PS点的高程残差和线性形变速率；

步骤5，通过利用相邻格网重叠区域PS点解算结果的差异作为观测值，以真实起算点所在格网的重叠区域PS解算结果作为起算数据，建立覆盖所有格网的误差方程和法方程，通过网平差，校准格网PS点相对高程残差和线性形变速率。

2.根据权利要求1所述的大范围高精度的InSAR形变监测数据处理方法，其特征在于，所述步骤2具体包括：

将所选PS点的像素位置和坐标信息将其归入不同的规则格网时，大气延迟在不足1km的两点大气相位差的方差不足0.1rad²，将研究区域划分成1km*1km的格网，并且相邻格网间设定20％的重叠，确保足够数量的PS点落在重叠区域。

3.根据权利要求1所述的大范围高精度的InSAR形变监测数据处理方法，其特征在于，所述步骤3具体包括：

根据阈值选取格网内的高等级PS点，剔除掉落在重叠区域的点，组成该格网内的候选PS参考点集；

分别以每一个候选参考点为中心，在候选参考点集中构建星型解算网，利用解空间搜索法估算出每一个弧段的高程残差之差Δε_h和线性形变速率之差Δν，并计算该弧段参数估计后平均残余相位

其中，i表示第i幅干涉图，t_i和分别表示该幅干涉图的时间基线和空间垂直基线；λ、R、θ分别表示波长、斜距以及入射角；N表示干涉对的数量，ΔΦ(x,y,i)表示弧段两端PS点差分干涉相位之差；

计算每一个星型网中所有弧段的残差均值，残差均值最小的星型网对应的中心候选参考点作为该格网的参考点，真实起算数据点所在的格网，以真实起算数据点作为该格网的参考点。

4.根据权利要求1所述的大范围高精度的InSAR形变监测数据处理方法，其特征在于，所述步骤5具体包括：

计算每一个重叠区域PS点解算结果之差，统计其分布情况，取其主要分布聚集区间的中值作为该重叠区域的差异值，即相邻格网的计算结果差异值；

以计算结果差异值作为观测值，以每个格网PS点解算结果整体改正量为代求参数，建立覆盖所有格网的平差网络，设定真实起算点所在格网的改正量为零，并以此为起算数据，进行最小二乘平差求解每个格网的校准改正量；

根据校准改正量对每一个格网估算的高程残差和香型形变速率进行整体平移校正，便得到整个研究区域最终连续的解算结果。

5.根据权利要求4所述的大范围高精度的InSAR形变监测数据处理方法，其特征在于，所述求解每个格网的校准改正量的步骤具体包括：

假定以部分重叠的方式将研究区域划分为M个格网，形成N个两两重叠的区域(N>M)，并附加真实起算点所在格网的改正量为零的限定，则有条件方程：

B·x-L＝0

其中B为参数x的系数矩阵，且有：

L＝[L₁,L₂,…,L_N,0]^T，x＝[x₁,x₂,…,x_M]^T

根据最小二乘平差原理，可以得到以下法方程：

B^TPBx＝B^TPL

其中，P为观测向量L的权矩阵，在此取单位权；

通过

计算出每个格网结算值的校准改正量。