一种广域PS-InSAR地表形变监测方法与流程

本发明涉及地表变形监测，尤其涉及一种广域ps-insar地表形变监测方法。

背景技术：

1、合成孔径雷达干涉测量(insar-interferometric synthetic aperture radar)是一种新型对地观测工具，具有全天时、全天候，不受云雨影响的特点，广泛用于大面积的地面高程及变形测量。永久散射体技术（ps-insar）是在insar技术基础上发展的一种长时间序列监测方法，其监测精度可达毫米级别，广泛应用于城市地表沉降、基础设施、滑坡的等应用场景。隨着近年来sar数据源的丰富、空间分辨率的提高，ps-insar监测能力得到极大提升。

2、然而，面对海量sar数据及大范围、复杂地域环境的精细化的时序变形监测需求，ps-insar面临新的挑战。ps-insar技术用于大范围区域、复杂地域形变监测目前主要有以下缺陷：

3、１）对于复杂地域的ps-insar变形监测，由于受大气相关距离限制以及模型相干系数阈值影响，连接点目标的delaunay三角网常常不能完整地连接所有点目标，出现若干个不连通的子网，导致无法获取部分点目标的形变信息，这种情况在非城镇地区影响尤为明显。

4、２）大范围海量sar数据对计算机硬件提出了更高的要求。随着sar影像时空分辨率的提升，同一区域选取的永久散射体及弧段成指数级增长，网平差运算产生大量的内存消耗且解算过程漫长。

5、３）针对大范围ps-insar的变形解算，传统方法是通过手动对数据分块的方法，但是由于不同分块内部永久散射体数量不固定，设定分块尺寸大小的方法很难达到自动化的监测要求。此外，简单的分块缺少整体的误差约束，导致解算误差、大气延迟等误差在分块间进行误差传递，出现明显的分块间形变速率系统性偏差，影响最终解算结果的精度。

6、因此，现有技术中将ps-insar应用于大范围和复杂区域的形变监测，受到算法和计算机硬件的局限，在现有条件和需求下难以实施；此外在大范围变形监测中难以实现自动化的分块变形解算及分块间的误差约束，从而影响了该技术的应用和推广。

技术实现思路

1、针对现有技术存在的缺陷，本发明提出一种广域ps-insar地表形变监测方法，通过对sar影像进行分块处理，大幅度地降低了解算所需的时间；解决了ps稀疏分布形成的多个子网的连通性差的问题，并充分利用块间重叠区域的信息，消除了自动化分块解算及分块结果间的系统性存在偏差的问题，用于大范围、复杂区域的地表形变监测中，降低了对计算机硬件设备的要求，同时还具有高效率和高精度的优势。

2、本发明提供了一种广域ps-insar地表形变监测方法，包括：

3、s1选取目标区域和时间段内的多幅sar影像，选取其中一幅sar影像作为主影像，将其他sar影像配准并采样到主影像网格内，生成辅影像集；

4、s2基于主影像和辅影像集生成多幅干涉图、多幅差分干涉图和多幅相干图；

5、s3利用多幅sar影像中各同名像元的振幅离差值和相干系数选取ps点；根据ps点的数量分布情况，依次对每幅sar影像进行分块处理，得到多个子块区域；其中相邻子块区域间设置预设面积的重叠区域；

6、s4依次将各子块区域内的ps点采用delaunay三角网构建第一高级弧段网络，所述第一高级弧段网络包括多个独立子网；采用弧段连通性分法计算独立子网的数量及独立子网上ps点的数量；在多个独立子网中筛选出预处理子网，依次将预处理子网边界上的ps点与其余的独立子网上最邻近的ps连接构建第二高级弧段网络；

7、s5依次将各子块区域内的第一和第二高级弧段网络合并得到第三高级弧段网络，采用弧段连通性分析法在第三高级弧段网络中选取连通域最大的一组独立子网作为ps形变解算网络，得到ps形变解算网络集合；

8、s6依次在各子块区域内的ps形变解算网络中选取任一ps点作为解算起点，以解算起点对应高级弧段的高程残差和线性形变速率为第一观测值，建立相应的ps形变解算网络的误差方程，通过加权最小二乘依次求解所有ps点的高程残差和线性形变速率，得到所有ps点的解算结果；

9、s7采用分块连通性分析法对各子块区域间的连通性进行分析，选取连通面积最大的子块区域集合，将该子块区域集合中重叠区域ps点的解算结果的差异作为第二观测值，以真实起算点所在子块区域中重叠区域ps点的解算结果作为起算数据，建立覆盖连通面积最大子块区域的误差方程，通过网平差校准所有ps点的高程残差和线性形变速率，得到地表形变监测结果。

10、具体地，步骤s3还包括：

11、s31对多幅sar影像进行相对辐射校正，计算相对辐射校正后的多幅sar影像中各同名像元的振幅离差值，将振幅离差值小于预设的离差阈值的同名像元作为候选点；

12、s32计算多幅干涉图的平均相干性，得到各候选点的相干系数，根据预设的相干系数阈值，将相干系数大于预设的相干系数阈值的候选点作为ps点；

13、s33根据ps点的数量分布情况设定分块尺寸；按照分块尺寸对sar影像进行分块处理，得到多个子块区域，使得各子块区域内ps点数量均小于预设的第一数量阈值，并在相邻子块区域间设置预设面积的重叠区域。

14、具体地，步骤s4还包括：

15、s41依次将各子块区域内的ps点采用delaunay三角网构建初始弧段网络，并计算初始弧段网络中各初始弧段的长度，去除弧段长度大于距离阈值的初始弧段，得到第一初级弧段网络；

16、s42利用解空间搜索算法依次解算第一初级弧段网络中各第一初级弧段的高程残余差、线性速率差和弧段模型相关系数，根据弧段模型相关系数剔除无效的第一初级弧段，得到第一高级弧段网络；所述第一高级弧段网络包括多个立子网；

17、s43采用弧段连通性分析法计算独立子网的数量及独立子网上ps点的数量；筛选出独立子网上ps点的数量大于第二数量阈值的独立子网作为预处理子网；

18、s44通过光栅扫描法依次识别各个预处理子网的边界，将各预处理子网的边界上的ps点与其余的独立子网上最邻近的ps点连接构建第二初级弧段网络，去除弧段长度大于距离阈值和无效的第二初级弧段，得到第二高级弧段网络。

19、具体地，步骤s7还包括：

20、s71采用分块连通性分析法判断子块区域间的连通性，选取连通面积最大的子块区域集合；

21、s72将该子块区域集合中重叠区域ps点的解算结果的差异作为第二观测值，以真实起算点所在子块区域中重叠区域ps点的解算结果作为起算数据，建立覆盖连通面积最大子块区域的误差方程；

22、s73通过网平差快速校准子块区域集合中的子块区域间的误差，实现高程残差和线性形变速率的拼接；

23、s74对所有各子块区域内的ps点的高程残差和线性形变速率进行整体平移校正，用整体平移校正后高程残差和线性形变速率的均值分别作为重叠区域的ps点的高程残差和线性形变速率，得到地表形变监测结果。

24、具体地，步骤s1还包括：

25、s11选取目标区域和时间段内的多幅sar影像；

26、s12根据每幅sar影像的时间基线、空间基线和多普勒质心频率差，计算每幅sar影像的整体相关系数，选择整体相关系数最大的影像作为主影像；

27、s13以主影像为基准影像，将其他影像都配准到基准sar影像的网格内，生成辅影像。

28、具体地，步骤s2还包括：

29、s21获取初始dem数据，将初始dem数据采样成与主影像一致的分辨率，并与主影像进行配准，并将初始dem范围裁剪成与主影像范围一致，生成覆盖目标区域的目标dem数据；

30、s22依次将主影像和辅影像集中的辅影像两两组队得到多个影像对，并对每对影像对进行前置滤波处理，对前置滤波处理后的影像对中的像元对进行共轭相乘，得多幅干涉图；

31、s23通过目标dem数据去除平地和地形相位获取差分干涉相位，依次在每幅干涉图中逐像元对去除差分干涉相位，得到多幅差分干涉图；

32、s24依次对每幅差分干涉图进行相位滤波处理，对相位滤波后的每幅差分干涉图对应的差分干涉相位像元对选择窗体大小，逐像元对计算其相干系数，得到多幅相干图。

33、具体地，步骤s5还包括：

34、s51将第一高级网络弧段网络和第二高级弧段网络合并得到第三初级弧段网络；

35、s52利用解空间搜索算法依次解算第三初级弧段网络中各第三初级弧段的弧段长度、高程残余差、线性速率差和弧段模型相关系数，去除弧段长度大于距离阈值的第三初级弧段，根据弧段模型相关系数剔除无效的第三初级弧段，得到第三高级弧段网络。

36、本发明的有益效果包括：

37、（1）本发明通过将目标区域的sar影像进行分块处理，得到多个子块区域，在多个子块区域中通过并行运算设计，大幅度降低了解算所需的时间，降低了对计算机硬件设备的要求；

38、（2）本发明充分利用相邻子块区域间重叠区域的信息，避免了分块解算过程中出现监测结果不连续的情况，消除了分块自动化分块解算及分块结果间的系统性偏差；

39、（3）本发明将通过不同步长扩展的最近邻点目标子网连接方法逐块连网的方式，对不连通的子网进行快速连接，减少子网数目，从而有利于大范围特别是ps点目标稀疏的非城镇地区地表形变监测，同时该方法易于实施，形变监测效率和精度均较高，扩展了insar干涉测量技术的应用范围。