基于地表应力应变模型的InSAR时序地表形变监测方法技术

本发明专利技术提供了一种基于地表应力应变模型的InSAR时序地表形变监测方法，包括：步骤1，收集覆盖待监测区域的时序SAR影像，基于现有软件实现时序SAR影像的配准、差分干涉，生成满足预设时空基线阈值的干涉图；步骤2，基于干涉图中的所有像素点构建局部狄洛尼三角网。本发明专利技术利用地表应力应变模型来描述地表临近点形变之间的物理力学关系，可显著提高InSAR时序地表形变监测的精度，基于邻近点组成的弧段进行建模和参数求解，因此无需对干涉图进行解缠，可有效提高InSAR数据处理效率，考虑了观测数据在空间上的物理力学关系，因此无需对干涉图中的低相干区域进行掩膜，即可获取空间覆盖密度大的高精度InSAR形变监测结果。

【技术实现步骤摘要】
基于地表应力应变模型的InSAR时序地表形变监测方法
本专利技术涉及基于卫星干涉合成孔径雷达的大地测量领域，特别涉及一种基于地表应力应变模型的InSAR时序地表形变监测方法。
技术介绍
合成孔径雷达干涉(InterferometricSyntheticApertureRadar,SAR，InSAR)技术利用两景SAR影像可进行地表形变测量。时序InSAR技术(MultipleTemporalInSAR，MTInSAR)对时序上的一系列SAR影像进行分析，可获取研究区域的时序地表形变结果，同时可有效削弱大气延迟，地形残差，失相干误差等因素对地表形变结果的影响。MTInSAR技术的监测目标主要分为永久散射体(PermanentScatterer，PS)和分布式散射体(DistributedScatterer，DS)。PS在整个检测周期内均保持较高的信噪比，但是在自然环境下，PS的分布密度极低，难以满足监测需求。DS在某些短时空基线干涉图上保持了较好的信噪比，并且相邻像元的时空变化特征较为相似，在自然环境下，DS分布密度较高(如稀疏植被区，沙漠等)，应用较为广泛。在获取DS目标处的时序地表形变时，现有方法往往需要基于所有可能干涉对进行复杂繁琐的数学最优化估计(即相位重构)，或者利用解缠后的短时空基线干涉图基于最小二乘准则逐点进行估计。实践证明，前者在估计地表形变时效率较低，而后者由于需要进行事先解缠，其监测结果易受解缠误差等因素影响。并且，现有方法在估计时序地表形变时，均未考虑邻近DS间的物理力学关系，且在低相干区域无法获得较为精确的形变结果。
技术实现思路
本专利技术提供了一种基于地表应力应变模型的InSAR时序地表形变监测方法，其目的是为了解决传统方法仅针对高相干点进行数学上的空间滤波，或者先对干涉图解缠、掩膜，再进行逐点的时序求解，均未考虑临近高相干点间的物理力学关系，且易导致监测结果的空间密度较低的问题。为了达到上述目的，本专利技术的实施例提供了一种基于地表应力应变模型的InSAR时序地表形变监测方法，包括：步骤1，收集覆盖待监测区域的时序SAR影像，基于现有软件实现时序SAR影像的配准、差分干涉，生成满足预设时空基线阈值的干涉图；步骤2，基于干涉图中的所有像素点构建局部狄洛尼三角网；步骤3，基于地表应力应变模型建立目标弧段的时序形变相位梯度与周围1km×1km范围内的所有弧段的干涉相位差之间的函数关系；步骤4，引入抗差估计剔除包含2π模糊度的弧段，最终基于加权最小二乘准则实现目标弧段时序形变相位差的解算；步骤5，对局部狄洛尼三角网的每一条边重复所述步骤3和所述步骤4，直至所有弧段上的时序形变相位差完成解算，以三角网中的一个地表稳定点或已知形变点为参考，对所有弧段进行空间积分，即可获得时序形变结果。其中，所述步骤1具体包括：收集覆盖待监测区域T景时序SAR影像，基于现有软件实现时序SAR影像的配准、差分干涉，最终生成满足预设时空基线阈值的M幅干涉图。其中，所述步骤2具体包括：基于M幅干涉图中的所有像素点，构建了局部狄洛尼三角网，局部狄洛尼三角网中共包含条弧段，最长的弧段应小于弧段长度预设阈值，弧段长度预设阈值的设定是根据假定同质区域的范围来确定的。其中，所述步骤2还包括：假设所构建的局部狄洛尼三角网中的某一弧段A0是以点Pi和点Pj为起止点，点Pi和点Pj所对应的三维坐标分别为和点Pi和点Pj处在第m幅干涉图所对应的SAR影像获取时刻之间发生的三维地表形变分别为和根据地表应力应变模型得到：其中，e代表东西向，n代表南北向，u代表垂直向，代表地表应力应变模型未知参数矩阵，可表示为：其中，代表求偏导数，dm＝[dm，edm，ndm，u]T代表第m幅干涉图所对应的主辅SAR影像获取时刻之间发生的三维地表形变场，x＝[xexnxu]T代表东西向、南北向和垂直向的三维方向。其中，所述步骤2还包括：由于点Pi和点Pj距离较近，SAR卫星对这两点进行观测时的几何视角差异可以忽略不计，即在点Pi和点Pj处，将三维形变投影到InSAR视线向形变的系数矩阵Bgeo是相同的，系数矩阵Bgeo如下所示：Bgeo＝[abc](3)其中，a代表东西向形变在InSAR视线方向形变的投影系数，b代表南北向形变在InSAR视线方向形变的投影系数，c代表垂直向形变在InSAR视线方向形变的投影系数，可根据SAR卫星的成像几何确定；在式(1)等号左右同乘以Bgeo，得到以下公式：其中，分别代表点Pi和点Pj处三维地表形变沿InSAR视线向的投影形变，即在第m幅干涉图中两点处的观测值，另外，可认为是InSAR一维形变观测值在弧段A0处沿东西向、南北向和垂直向的形变梯度参数。其中，所述步骤3具体包括：基于地表应力应变模型可知，在第m幅干涉图中，弧段A0一定空间范围内的地表形变梯度可以假设为恒定的，因此，基于式(4)，即可建立第m幅干涉图中目标弧段A0上地表形变梯度与周围一定范围内的所有K个弧段Ak，包含弧段A0的形变相位差之间的函数关系，其中，(k＝1，2，...，K)，得到：其中，Bsm＝[Δ1，Δ2，...，Δk，…，ΔK]T，Δk代表弧段Ak两个端点间的坐标增量，代表弧段Ak两个端点间InSAR形变观测值的差值。其中，所述步骤3还包括：以M幅干涉图的第一景SAR影像为参考，易得InSAR干涉图对应的形变观测值与InSAR视线向时序地表形变之间的转换矩阵Btrs，其中，Btrs的矩阵大小为M×(T-1)，每一行对应一个干涉图，每一列对应一个SAR影像，除去参考影像，每一行中，M幅干涉图主影像对应的列为-1，辅影像对应的列为+1，其他元素为0；基于式(5)和矩阵Btrs，即可建立M幅干涉图中所有弧段上的形变相位差ΔLdefo与弧段A0上的时序地表形变相位梯度之间的函数关系：其中，ΔLdefo＝[(ΔL1)T，(ΔL2)T，…，(ΔLm)T，…，(ΔLM)T]T，其中，t＝2，3，...，T，代表克罗内克积运算符号，知分别代表第t个SAR影像获取时InSAR视线向形变相位在弧段A0处沿东西向、南北向和垂直向的梯度；目标弧段A0上的实际InSAR观测值包含地表形变相位差包含地形残差相关相位及噪声相位，其中，目标弧段A0上的地形残差相位表示为：其中，其中，表示第m幅干涉图中弧段A0上的地形残差相位，代表第m幅干涉图的垂直基线，dz0表示弧段A0上的地形残差，λ，ρ，θ分别代表卫星的波长，卫星到目标点距离，目标点处的卫星入射角。其中，所述步骤3还包括：假定地形残差相关相位在空间上是不相关的，则弧段Ak上的地形残差相位可认为是噪声，结合式(6)和式(7)，基于地表应力应变模型，即可建立目标弧段的时序形变相位梯度与周围一定范围内本文档来自技高网...

【技术保护点】
1.一种基于地表应力应变模型的InSAR时序地表形变监测方法，其特征在于，包括：/n步骤1，收集覆盖待监测区域的时序SAR影像，基于现有软件实现时序SAR影像的配准、差分干涉，生成满足预设时空基线阈值的干涉图；/n步骤2，基于干涉图中的所有像素点构建局部狄洛尼三角网；/n步骤3，基于地表应力应变模型建立目标弧段的时序形变相位梯度与周围1km×1km范围内的所有弧段的干涉相位差之间的函数关系；/n步骤4，引入抗差估计剔除包含2π模糊度的弧段，最终基于加权最小二乘准则实现目标弧段时序形变相位差的解算；/n步骤5，对局部狄洛尼三角网的每一条边重复所述步骤3和所述步骤4，直至所有弧段上的时序形变相位差完成解算，以三角网中的一个地表稳定点或已知形变点为参考，对所有弧段进行空间积分，即可获得时序形变结果。/n

【技术特征摘要】
1.一种基于地表应力应变模型的InSAR时序地表形变监测方法，其特征在于，包括：
步骤1，收集覆盖待监测区域的时序SAR影像，基于现有软件实现时序SAR影像的配准、差分干涉，生成满足预设时空基线阈值的干涉图；
步骤2，基于干涉图中的所有像素点构建局部狄洛尼三角网；
步骤3，基于地表应力应变模型建立目标弧段的时序形变相位梯度与周围1km×1km范围内的所有弧段的干涉相位差之间的函数关系；
步骤4，引入抗差估计剔除包含2π模糊度的弧段，最终基于加权最小二乘准则实现目标弧段时序形变相位差的解算；
步骤5，对局部狄洛尼三角网的每一条边重复所述步骤3和所述步骤4，直至所有弧段上的时序形变相位差完成解算，以三角网中的一个地表稳定点或已知形变点为参考，对所有弧段进行空间积分，即可获得时序形变结果。


2.根据权利要求1所述的基于地表应力应变模型的InSAR时序地表形变监测方法，其特征在于，所述步骤1具体包括：
收集覆盖待监测区域T景时序SAR影像，基于现有软件实现时序SAR影像的配准、差分干涉，最终生成满足预设时空基线阈值的M幅干涉图。


3.根据权利要求2所述的基于地表应力应变模型的InSAR时序地表形变监测方法，其特征在于，所述步骤2具体包括：
基于M幅干涉图中的所有像素点，构建了局部狄洛尼三角网，局部狄洛尼三角网中共包含条弧段，最长的弧段应小于弧段长度预设阈值，弧段长度预设阈值的设定是根据假定同质区域的范围来确定的。


4.根据权利要求3所述的基于地表应力应变模型的InSAR时序地表形变监测方法，其特征在于，所述步骤2还包括：
假设所构建的局部狄洛尼三角网中的某一弧段A0是以点Pi和点Pj为起止点，点Pi和点Pj所对应的三维坐标分别为和点Pi和点Pj处在第m幅干涉图所对应的SAR影像获取时刻之间发生的三维地表形变分别为和根据地表应力应变模型得到：



其中，e代表东西向，n代表南北向，u代表垂直向，代表地表应力应变模型未知参数矩阵，可表示为：



其中，代表求偏导数，dm＝[dm，edm，ndm，u]T代表第m幅干涉图所对应的主辅SAR影像获取时刻之间发生的三维地表形变场，x＝[xexnxu]T代表东西向、南北向和垂直向的三维方向。


5.根据权利要求4所述的基于地表应力应变模型的InSAR时序地表形变监测方法，其特征在于，所述步骤2还包括：
由于点Pi和点Pj距离较近，SAR卫星对这两点进行观测时的几何视角差异可以忽略不计，即在点Pi和点Pj处，将三维形变投影到InSAR视线向形变的系数矩阵Bgeo是相同的，系数矩阵Bgeo如下所示：
Bgeo＝[abc](3)
其中，a代表东西向形变在InSAR视线方向形变的投影系数，b代表南北向形变在InSAR视线方向形变的投影系数，c代表垂直向形变在InSAR视线方向形变的投影系数，可根据SAR卫星的成像几何确定；
在式(1)等号左右同乘以Bgeo，得到以下公式：



其中，分别代表点Pi和点Pj处三维地表形变沿InSAR视线向的投影形变，即在第m幅干涉图中两点处的观测值，另外，可认为是InSAR一维形变观测值在弧段A0处沿东西向、南北向和垂直向的形变梯度参数。


6.根据权利要求5所述的基于地表应力应变模型的InSAR时序地表形变监测方法，其特征在于，所述步骤3具体包括：
基于地表应力应变模型可知，在第m幅干涉图中，弧段A0一定空间范围内的地表形变梯度可以假设为恒定的，因此，基于式(4)，即可建立第m幅干涉图中目标弧段A0上地表形变梯度与周围一定范围内的所有K个弧段Ak，包含弧段A0的形变相位差之间的函数关系，其中，(k＝1，2，...，K)，得到：



其中，Bsm＝[Δ1，Δ2，...，Δk，...，ΔK]T，Δk代表弧段Ak两个端点间的坐标增量，代表弧段Ak两个端点间InSAR形变观测值的差值。


7.根据权利要求6所述的基于地表应力应变模型的InSAR时序地表形...

【专利技术属性】
技术研发人员：胡俊，刘计洪，孙倩，李志伟，朱建军，
申请(专利权)人：中南大学，湖南师范大学，
类型：发明
国别省市：湖南;43