CN109061641A - 一种基于序贯平差的InSAR时序地表形变监测方法 - Google Patents

Abstract

本发明公开了一种基于序贯平差的InSAR时序地表形变监测方法，首先获取研究区域的历史形变结果，包括时序形变，形变速率和地形残差；当有新影像观测数据加入时，利用序贯平差结合新InSAR观测数据及历史数据结果对形变成果进行更新以达到整体解算的目的。该方法突破了常规计算过程，巧妙的利用已经求解结果的历史观测数据，当有新的观测数据加入时，无需融合所有数据进行整体解算，利用序贯平差的计算思想，仅利用历史解算结果作为解算基础，结合新观测数据进行解算，即可达到整体解算的目的，其计算效率大大提高。

Description

一种基于序贯平差的InSAR时序地表形变监测方法

技术领域

本发明属于基于遥感影像的大地测量领域，特别涉及一种基于序贯平差的InSAR时序地表形变监测方法。

背景技术

合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar,SAR，InSAR)在监测地表形变方面具有全天候，全天时，大范围，高精度与高空间分辨率等特点，已经被广泛应用于监测地表形变领域。然而InSAR技术易受时空失相关，大气误差及高程残差等影响，使得InSAR技术在应用过程中受到一定限制。多时相InSAR(Multi-temporalInSAR，MT-InSAR)技术对同一研究区域时域上的一系列差分干涉影像进行整体处理，以达到消除或减弱相关误差对地表形变监测的影响。

随着技术的发展，SAR卫星重访周期不断缩短，分辨率不断提高。海量SAR数据为MT-InSAR提供数据基础的同时，也给快速、高效的形变结果解算带来了挑战。对于传统MT-InSAR技术，每获取一景新影像时，都需将新影像和历史影像进行整体解算，计算效率较低，冗余度大。

发明内容

本发明提出了一种基于序贯平差的InSAR时序地表形变监测方法，该方法在动态InSAR时序地表形变监测中可实现快速、高效的形变结果解算。

一种基于序贯平差的InSAR时序地表形变监测方法，包括以下步骤：

步骤1：获取待监测的地表区域N+1幅时序SAR影像，并获得对应的M幅历史解缠差分干涉相位图；

步骤2：基于历史解缠差分干涉相位结合和W⁰，求解待监测区域的历史时序形变速率v^p,0、历史地形残差

其中，表示地形残差求解系数矩阵，Δt⁰为历史干涉图时间间隔向量，H⁰为历史干涉图对应的高程转换系数向量，W⁰为M×M的单位矩阵；

步骤3：利用去除地形残差相位的M幅历史解缠差分干涉相位结合和W⁰求解历史时序形变相位从而获得历史时序形变

其中，λ表示雷达波长，表示时序形变相位求解系数矩阵，大小为M*N，M是历史解缠差分干涉相位图的数量，矩阵的第k行表示第k幅干涉图，第k行中的第i个元素值为-1，第j个元素值为1，其余元素值均为0，第k幅干涉图是由第i，j两幅SAR影像干涉生成；

先假定时序形变模型为线性模型，建立时序形变速率、地形残差与InSAR数据的观测方程，此时即可求得时序形变速率与地形残差；

在InSAR观测值中将高程残差相位去除，结合最小二乘准则，利用一系列去除地形残差的InSAR观测值求解时序形变。

步骤4：获取新增SAR影像，并从新增的SAR影像中选取满足时空基线阈值的SAR影像对，获取新增的M₁幅解缠差分干涉相位图；

步骤5：利用新增的M₁幅解缠差分干涉相位和待监测区域的W⁰，基于序贯平差求得历史地形残差相关未知参数向量的改正数向量进而得到更新后的地形残差相关未知参数向量即获得了待监测区域更新后的时序形变速率v^p,a和地形残差

其中，表示历史地形残差相关未知参数向量， 表示更新后的地形残差相关未知参数向量，

步骤6：利用W⁰和对历史时序形变相位进行初步更新，得到时序形变相位初步更新值

步骤7：基于去除地形残差相位的M₁幅新增解缠差分干涉相位和初步更新的时序形变相位利用序贯平差对地表时序形变相位进行更新，得到待监测区域的近实时时序形变相位从而获得对应的时序形变

利用地形残差乘上高程转换系数向量，可以计算出高程残差相关的相位，只要得到了历史地形残差之后，就简单快速的算出去除地形残差的观测值，故得到了新的地形残差就可以得到去除地形残差的新观测值；

进一步地，所述历史地形残差相关未知参数向量的改正数向量按以下公式计算获得：

其中，Δv^p,a和分别表示待监测区域历史时序形变速率和地形残差的改正数向量；J₂、表示中间变量；

表示新增M₁幅解缠差分干涉相位图对应的地形残差求解系数矩阵， Δt^a为新增M₁干涉图对应的时间间隔向量，H^a为新增M₁干涉图对应的高程转换系数向量；

W^a表示新增M₁幅解缠差分干涉相位图的权重，取值为M₁×M₁的单位矩阵。

进一步地，所述时序形变相位初步更新值的计算过程如下：

其中，表示表示p点处第k幅干涉图的垂直基线，R^p表示p点到卫星的斜距，θ^p表示p点处的雷达入射角。

在已有的时序形变速率和高程残差的基础上，建立时序形变速率、地形残差与新InSAR数据的观测方程，此时利用序贯平差即可求得时序形变速率与地形残差改正数，利用该改正数对时序形变速率，地形残差进行更新，同时对历史时序形变相位进行初步更新；在新InSAR观测值中将更新后的地形残差相位去除，基于已初步更新的时序形变相位和新的去除地形残差的InSAR观测值，利用序贯平差即可求解得到最终的时序形变结果。当有新影像观测数据加入时，利用序贯平差可大大提高计算效率。

有益效果

本发明提供了一种基于序贯平差的InSAR时序地表形变监测方法，首先获取研究区域的历史形变结果，包括时序形变，形变速率和地形残差；当有新影像观测数据加入时，利用序贯平差结合新InSAR观测数据及历史数据结果对形变成果进行更新以达到整体解算的目的。该方法突破了常规计算过程，巧妙的利用已经求解结果的历史观测数据，当有新的观测数据加入时，无需融合所有数据进行整体解算，利用序贯平差的计算思想，仅利用历史解算结果作为解算基础，结合新观测数据进行解算，即可达到整体解算的目的，其计算效率大大提高。随着InSAR大数据时代的来临，本发明可充分发挥InSAR监测地表形变的时效性。并且，基于本发明方法也可以大大降低InSAR大数据量下对计算机性能的要求，为InSAR技术的推广及应用奠定了良好基础。

附图说明

图1是基于序贯平差的InSAR时序地表形变监测方法的流程图；

图2是模拟时序InSAR数据的时空基线图，三角形和圆圈代表SAR影像获取的相对时空位置，连线代表对应SAR影像组成干涉对，三角形及实线代表历史已有数据，圆圈及虚线代表新增观测数据；

图3是利用序贯平差和传统整体平差结果对比图，直线代表模拟的时序形变，“*”代表传统方法得到的时序形变结果，“o”代表序贯平差得到的时序形变结果；斜杠前数据代表对应方法求解过程的时间消耗，斜杠后数据代表对应方法所得结果相比于模拟数据的均方根误差；

图4是利用序贯平差和传统整体平差过程中运算效率对比图，圆圈代表50次模拟实验中传统方法求解形变过程所用的时间，三角形代表对应的序贯平差所用的时间。

具体实施方式

为了使本技术相关领域的人员能够更好地理解本发明，下面将对本发明的实施方案进行清楚、详细的描述。

一种基于序贯平差的InSAR时序地表形变监测方法，如图1所示，为本发明的流程图，包括以下步骤：

步骤1：历史解缠差分干涉影像集数据获取。选取关于研究区域同一轨道的N+1幅时序SAR影像，对应的时刻分别为t₀,t₁,t₂…t_N，并将其配准重采样至同一影像坐标系下；对符合时空基线阈值的影像对分别进行干涉、差分、解缠过程得到M幅解缠差分干涉相位图；结合解缠差分干涉图及相应的相干性图挑选出高相干点作为解算对象。

步骤2：求解地形残差解缠差分干涉相位图均是基于已有的数字高程模型进行地形相位改正，但是已有的数字高程模型数据往往与真实高程不符，这在一定程度上降低了利用InSAR技术监测地表形变的精度及可靠性。一系列的解缠差分干涉相位中包含影像获取期间的累计地表形变和地形残差信息，若假定一定的时序形变模型(如线性形变)，利用一系列的解缠差分干涉相位，即可求解出相应的时序形变速率及地形残差。

假设第p个高相干点在第k幅干涉图中的解缠相位为(k＝1,2…M，上标0代表该变量为历史解算过程变量)，其中包含时间t_i到t_j的累计形变相位(0≤i≤N-1,1≤j≤N)及地形残差相位即可写成：

假设地表发生时域上的线性形变，即t_i时刻相对于t₀时刻的累计形变相位t_i,v^p,0代表时序形变速率，那么根据SAR卫星成像几何可得其中代表第k幅干涉图的垂直空间基线，代表第p个高相干点处的DEM残差，λ代表SAR卫星的波长，R^p代表卫星到地面点的斜距，θ^p代表第p个点上卫星的入射角。即：

考虑共有M幅干涉图可得：

其中

令并将M幅干涉图做等权处理，即权重W⁰＝I，I为M×M的单位矩阵，基于加权最小二乘准则可得未知参数初始解：

此时即可得到第p个高相干点的时序形变速率v^p,0和地形残差

步骤3：利用去除地形残差相位的M幅解缠差分干涉相位图进行时序形变求解。

在解缠差分干涉相位中去除地形残差相关相位可得到仅包含形变信息和随机误差的解缠差分干涉相位(不考虑大气等其他误差项)，利用结合加权最小二乘准则即可求得时序形变相位。

其中：

设为大小为M×N的矩阵，其第k行代表的第k幅干涉图由第i,j(0≤i<j≤N)幅SAR影像干涉生成，则对应矩阵的第k行的第i个元素(i>0)为-1，第j个元素为1，其余元素均为0。

进而基于加权最小二乘准则即可求解时序形变相位：

其中

代表时序形变相位，其中序列中的每一个元素代表对应时间序列中每一个时刻相对于t₀的累计形变相位，即时序形变

通过步骤1-3即可基于历史已有数据进行InSAR时序地表形变监测，解算出历史数据结果(即时序形变时序形变速率v^p,0和地形残差)。

对应流程图的d0，v^p,0对应v0和对应h0；

步骤4：新增影像数据。利用新增N₃幅SAR影像(对应时刻分别为)以及与新增影像满足时空基线阈值的历史N₂幅SAR影像进行类似步骤1的处理方式得到新增的M₁幅解缠差分干涉相位图。

步骤5：更新地形残差值。此步骤将步骤2中的W⁰和新增解缠差分干涉相位(带有上标a的变量表示新增(added)解缠差分干涉相位之后序贯平差的中间变量)及其权阵W^a作为输入参数进行时序形变速率和地形残差的更新

同样假设第p个高相干点在第r幅干涉图中的解缠相位为(r＝1,2…M₁),其中包含时间t_u到t_v的累计形变相位(0≤u≤N+N₃-1,N+1≤v≤N+N₃)及地形残差相位即可写成：

同步骤2可得

即：

考虑共有M₁幅干涉图可得：

其中

Δt_r＝t_v-t_u

令考虑已求解的未知参数可得误差方程为：

其中表示观测值改正数向量，代表未知参数的改正数向量。

令

其中W^a为新增M₁幅干涉图的权阵，一般设置为M₁×M₁的单位矩阵，则可得到时序形变速率和地形残差改正数向量

即经过序贯平差更新后的未知参数解为：

此时即可求得第p个高相干点的更新后的时序形变速率v^p,a和地形残差

令

其中，A₀为M₁×M大小的零矩阵，即可将时序形变速率与地形残差向量对应的系数矩阵以及对应观测值的权阵W作为下次确定时序形变速率及地形残差的序贯平差输入参数。

步骤6：利用地形残差改正数对步骤3得到的时序形变相位进行初步更新；

此步骤将步骤3中的W⁰和步骤5中的为输入参数更新得到时序形变相位若对历史数据和新增数据进行传统整体平差，则用于求解地表时序形变的观测值为去除地形残差更新解之后的解缠差分干涉相位，而单对历史数据进行平差得到时序形变所用的观测值是去除地形残差初始解之后的解缠差分干涉相位。为了使序贯平差解算结果达到整体平差的目的，需对由历史数据解算得到的时序形变相位进行改正。

改正后的时序形变相位

其中为步骤5中所得的地形残差改正数相关的相位

步骤7：基于去除地形残差相位的M₁幅新增解缠差分干涉相位和步骤6更新的时序形变相位利用序贯平差对地表时序形变相位进行更新

此步骤将步骤3中的W⁰，步骤5中的步骤6中的和新增解缠差分干涉相位及其权阵W^a作为输入参数得到序贯平差时序形变相位结果在解缠差分干涉相位中去除地形残差相关相位可得到仅包含形变信息和随机误差的解缠差分干涉相位利用和历史时序形变相位结果结合序贯平差即可求得最终时序形变相位解。

其中：

假设新增M₁幅干涉图是由N₂景历史SAR影像和N₃幅新增SAR影像干涉生成，那么则有N₁(N₁＝N-N₂)景历史SAR影像未参与新增M₁幅干涉图的组成。同时令N₁景未参与新增M₁幅干涉图组成的历史SAR影像对应的时刻集合为t_Ⅰ，N₂景参与新增M₁幅干涉图组成的历史SAR影像对应的时刻集合为t_Ⅱ，N₃幅新增SAR影像对应的时刻集合为t_Ⅲ。如步骤3所述，矩阵中的每一列对应每一个时刻的参数，这里将矩阵中对应t_Ⅰ时刻的列依次提取组成矩阵对应t_Ⅱ时刻的列依次提取组成矩阵同样历史时序形变相位解算结果向量中对应t_Ⅰ时刻的元素依次提取组成向量对应t_Ⅱ时刻的元素依次提取组成向量

同步骤3类似，创建系数矩阵其大小为M₁×(N₂+N₃)的矩阵，其第r行代表的新增第r幅干涉图由第t_u和t_v(0≤u≤N+N₃-1,N+1≤v≤N₁+N₂+N₃)时刻的SAR影像干涉生成，则对应矩阵的第r行的第u个元素(u>0)为-1，第v个元素为1，其余元素均为0。此时将中对应t_Ⅱ时刻的列依次提取组成矩阵对应t_Ⅲ时刻的列依次提取组成矩阵

此时令

即可得时刻t_Ⅱ和t_Ⅲ对应的时序形变相位的改正数为

中包含的t_Ⅱ和t_Ⅲ时刻对应的时序形变相位改正数分别为和则此时t_Ⅰ对应的时序形变相位的改正数为

此时，可得最终由序贯平差更新得到的时序形变相位向量为

即卫星视线方向的地表时序形变向量

令A₁为M×N₃大小的零矩阵，

基于历史SAR数据得到的地表形变结果，当获取新的SAR影像之后，针对每一个高相干点p进行上述步骤4-7，即可得到每个高相干点最终的时序形变，时序形变速率及地形残差。

每当有新观测数据加入时，将新的观测数据以及上次平差得到的W作为历史解算结果(输入参数)重复上述步骤4-7，以此类推，即可实现快速、高效地动态地表形变解算(包括地表时序形变，地表形变速率及地形残差)，为由地表形变引起的地质灾害预测和解译提供有力保障。

本发明效果可以通过以下模拟实验进一步说明。

模拟数据描述：①结合哨兵-1A的数据，模拟出时间轴上间隔为12天的91个时刻(2017年4月3日至2020年3月18日)，及相应的空间位置，即可认为模拟了91景SAR影像；②设定时空基线阈值分别为40天和100米，得到了236幅干涉对(其时空基线如图2)。③假定该序列SAR影像为升轨右视影像，方位角为-12°，入射角为33.8°，东西向、南北向和垂直向形变速率均为1mm/天，高程残差为10米，卫星到地面点斜距为800千米，电磁波波长为5.6厘米，根据卫星成像几何等信息即可得到236个形变数据；④在得到的236个形变数据中加入标准差为3毫米的高斯噪声作为模拟实验的观测数据。其中以2020年1月18日之前组成的干涉对(即前224个数据)作为历史观测数据，之后得到的干涉对数据作为新增观测数据(共12个数据)。

当新增观测数据时，传统方法是融合历史及新增数据进行整体解算，效率较低、冗余度大。本次模拟实验分别利用传统解法与序贯平差方法对模拟数据进行时序地表时序形变的解算，为了更好地说明问题，本次重复进行了50次模拟实验。在每次模拟实验中，两种方法得到的地表形变及高程残差结果均相同，如图3，是其中某一次模拟实验得到的结果图。

另外，在50次模拟实验中利用Matlab中的计时函数统计了每次实验运行的时间，如图4。

由模拟实验可知，相比于传统算法，本发明提及的算法在保证计算精度的同时，大大提高了InSAR时序地表形变计算的效率，为由地表形变引起的地质灾害预测和解译提供有力保障。

Claims (3)

1.一种基于序贯平差的InSAR时序地表形变监测方法，其特征在于，包括以下步骤：

步骤1：获取待监测的地表区域N+1幅时序SAR影像，并获得对应的M幅历史解缠差分干涉相位图；

步骤2：基于历史解缠差分干涉相位结合和W⁰，求解待监测区域的历史时序形变速率v^p,0、历史地形残差

其中，表示地形残差求解系数矩阵，Δt⁰为历史干涉图时间间隔向量，H⁰为历史干涉图对应的高程转换系数向量，W⁰为M×M的单位矩阵；

步骤3：利用去除地形残差相位的M幅历史解缠差分干涉相位结合和W⁰求解历史时序形变相位从而获得历史时序形变

其中，λ表示雷达波长，表示时序形变相位求解系数矩阵，大小为M*N，M是历史解缠差分干涉相位图的数量，矩阵的第k行表示第k幅干涉图，第k行中的第i个元素值为-1，第j个元素值为1，其余元素值均为0，第k幅干涉图是由第i，j两幅SAR影像干涉生成；

步骤4：获取新增SAR影像，并从新增的SAR影像中选取满足时空基线阈值的SAR影像对，得到新增的M₁幅解缠差分干涉相位图；

步骤5：利用新增的M₁幅解缠差分干涉相位和待监测区域的W⁰，基于序贯平差求得历史地形残差相关未知参数向量的改正数向量进而得到更新后的地形残差相关未知参数向量即获得了待监测区域更新后的时序形变速率v^p,a和地形残差

其中，表示历史地形残差相关未知参数向量， 表示更新后的地形残差相关未知参数向量，

步骤6：利用W⁰和对历史时序形变相位进行初步更新，得到时序形变相位初步更新值

步骤7：基于去除地形残差相位的M₁幅新增解缠差分干涉相位和初步更新的时序形变相位利用序贯平差对地表时序形变相位进行更新，得到待监测区域的实时时序形变相位从而获得对应的时序形变

2.根据权利要求1所述的方法，其特征在于，所述历史地形残差相关未知参数向量的改正数向量按以下公式计算获得：

其中，Δv^p,a和分别表示待监测区域历史时序形变速率和地形残差的改正数向量；J₂、表示中间变量；

表示新增M₁幅解缠差分干涉相位图对应的地形残差求解系数矩阵， Δt^a为新增M₁干涉图对应的时间间隔向量，H^a为新增M₁干涉图对应的高程转换系数向量；

W^a表示新增M₁幅解缠差分干涉相位图的权重，取值为M₁×M₁的单位矩阵。

3.根据权利要求2所述的方法，其特征在于，所述时序形变相位初步更新值的计算过程如下：

其中，表示表示p点处第k幅干涉图的垂直基线，R^p表示p点到卫星的斜距，θ^p表示p点处的雷达入射角。