基于InSAR相位梯度堆叠的广域地表形变探测方法及系统

基于insar相位梯度堆叠的广域地表形变探测方法及系统
技术领域
1.本发明涉及雷达干涉测量技术领域与地质灾害监测技术领域，提出实现了一种基于insar相位梯度堆叠的广域地表形变探测方案。


背景技术：

2.合成孔径雷达干涉测量(interferometric synthetic aperture radar,insar)技术是上世纪90年代发展起来的空间大地测量技术。insar技术的优势在于全天时、全天候、高精度、高效率、低成本、可重复观测，特别对于自然条件恶劣以及大区域的数据获取更具优势。因此，insar技术为获取大范围的高精度地表形变信息提供了一种全新的技术手段。insar技术已被成功应用于火山运动、地震活动、山体滑坡、冰川漂移、区域地面沉降、基础设施安全监测等领域。
3.然而，常规的差分干涉测量(dinsar)会受到很多因素的影响，如大气延迟、dem误差，时间失相干等，使得最终形变测量结果的准确性受到影响。为克服这些不利因素的影响，时序insar技术应运而生并在监测各种类型的地表形变中取得了巨大成功。时序insar技术是通过分析大量sar数据，关注高相干目标，有效抑制了常规dinsar技术中的各种误差的影响，提高了形变测量精度。然而，在一些地形复杂或失相干严重的区域，可能存在insar相位解缠困难的问题；并且面对大范围研究区域时，时序insar数据处理效率很低；虽然时序insar技术能够在一定程度上削弱大气延迟扰动的影响，但普遍存在算法复杂、不具有普适性等特点。所以，时序insar在大范围地表形变测量中仍然存在较多问题。
4.传统的insar数据处理中相位解缠是其中一个重要的步骤，例如已有专利文献cn111998766a、cn112051571a、cn113340191a，但是在一些低相干或地表复杂区域往往造成相位数据不连续，给相位解缠工作带来极大的困难，解缠的结果往往带有很大的误差。


技术实现要素：

5.本发明目的在于提出一种基于insar差分干涉相位梯度堆叠进行大范围地表形变探测的方案，基于大气延迟噪声在局部空间尺度上相关和时间上不相关的特点，对差分干涉相位进行空间上差分，在时间上进行堆叠，可以有效降低大气延迟噪声对形变探测的影响。干涉相位梯度堆叠对局部形变信号比较敏感，无需进行相位解缠，可以快速识别大范围区域内的局部形变信号，可用于大范围地质灾害普查。
6.为实现上述目的，本发明采用的技术方案为一种基于insar相位梯度堆叠的广域地表形变探测方法，基于insar差分干涉相位梯度堆叠支持大范围地表形变探测，实现过程包括以下步骤，
7.步骤一、通过获取覆盖研究区域的多幅sar影像，选取主影像并将其他影像与主影像配准，通过设置时空基线阈值，进一步处理获得干涉图；
8.步骤二、根据sar影像之间的时间基线和空间基线挑选干涉对，采用外部dem模拟地形相位，并从干涉图中减去地形相位得到差分干涉图，并进行空间自适应滤波处理获取
滤波后的差分干涉图；
9.步骤三、对步骤二挑选并处理后的差分干涉图，在像元级上按照预设的步长分别进行像元周围多个方向的梯度计算，获得差分干涉图对应的各方向梯度图；
10.设差分干涉相位对应的某一方向d梯度为表示如下：
[0011][0012]
其中，为sar影像经过处理得到的第i幅差分干涉图，m和n分别代表干涉图中像元所在的行列号，为差分干涉图对应的某一方向d的梯度图，s为步长；
[0013]
步骤四、对获取的各个方向梯度图分别进行堆叠以增强形变信号并再次进行滤波；分别将滤波后的各方向计算的相位梯度堆叠图进行合并，获得研究区域梯度堆叠结果；
[0014]
步骤五、对获取的研究区域梯度结果映射到预设的范围得到最终的研究区域梯度堆叠结果。
[0015]
而且，设求取东、西、南、北、东南、东北、西南、西北方向的相位梯度，东、西、南、北、东南、东北、西南、西北方向的相位梯度分别用东南、东北、西南、西北方向的相位梯度分别用表示。
[0016]
而且，将各方向梯度图进行堆叠以增强形变信号时，进行包括以下步骤，
[0017]
对于不同时间基线的梯度图上的沿某一方向d，间隔预设步长的两个像元，其干涉相位梯度堆叠公式如下：
[0018][0019]
上式中，φd为获取的n个不同时间基线差分干涉图对应的沿某一方向间隔预设步长的梯度图叠加的结果，为第i幅差分干涉图对应的沿某一方向d获得的梯度，i＝1，2，3,
…
，n，n为干涉对总数；
[0020]
然后，对获得的梯度堆叠结果φd采用空间滤波方法进行滤波，获得滤波后的结果用ψd表示，东、西、南、北、东南、东北、西南、西北方向的相位梯度堆叠的滤波结果分别用ψe、ψw、ψs、ψn、ψ
en
、ψ
es
、ψ
wn
、ψ
ws
表示。
[0021]
而且，将获得的空间滤波后的各方向梯度堆叠图进行合并以获取最终的差分干涉相位梯度堆叠信息，其相应公式表示如下：
[0022][0023]
其中，ψ为合并各个方向干涉相位梯度图的结果，总的方向个数为8。
[0024]
而且，对获取的研究区域干涉相位梯度的结果值映射到预设的范围后，在相位梯度堆叠图上出现聚集的与周边对比明显的大梯度值即为探测所得的地表形变信息。
[0025]
另一方面，本发明提供一种基于insar相位梯度堆叠的广域地表形变探测系统，用于实现如上所述的一种基于insar相位梯度堆叠的广域地表形变探测方法。
[0026]
而且，用于基于insar差分干涉相位梯度堆叠支持大范围地表形变探测，包括以下模块，
[0027]
第一模块，用于通过获取覆盖研究区域的多幅sar影像，选取主影像并将其他影像与主影像配准，通过设置时空基线阈值，进一步处理获得干涉图；
[0028]
第二模块，用于根据sar影像之间的时间基线和空间基线挑选干涉对，采用外部dem模拟地形相位，并从干涉图中减去地形相位得到差分干涉图，并进行空间自适应滤波处理获取滤波后的差分干涉图；
[0029]
第三模块，用于对第二模块挑选并处理后的差分干涉图，在像元级上按照预设的步长分别进行像元周围多个方向的梯度计算，获得差分干涉图对应的各方向梯度图；
[0030]
设差分干涉相位对应的某一方向d梯度为表示如下：
[0031][0032]
其中，为sar影像经过处理得到的第i幅差分干涉图，m和n分别代表干涉图中像元所在的行列号，为差分干涉图对应的某一方向d的梯度图，s为步长；
[0033]
第四模块，用于对获取的各个方向梯度图分别进行堆叠以增强形变信号并再次进行滤波；分别将滤波后的各方向计算的相位梯度堆叠图进行合并，获得研究区域梯度堆叠结果；
[0034]
第五模块，对获取的研究区域梯度结果映射到预设的范围得到最终的研究区域梯度堆叠结果。
[0035]
或者，包括处理器和存储器，存储器用于存储程序指令，处理器用于调用存储器中的存储指令执行如上所述的一种基于insar相位梯度堆叠的广域地表形变探测方法。
[0036]
或者，包括可读存储介质，所述可读存储介质上存储有计算机程序，所述计算机程序执行时，实现如上所述的一种基于insar相位梯度堆叠的广域地表形变探测方法。
[0037]
本发明的创新之处在于无需相位解缠，采用insar差分干涉相位梯度堆叠的方法识别研究区的形变信号，不仅有效抑制大气延迟噪声的影响，而且可以在大范围的研究区域快速的识别形变信号，从而可以及时的进行地质灾害的普查。
[0038]
本发明提出的方案没有直接对差分干涉图进行相位解缠及其后续相位叠加或时间序列insar分析等过程，而是直接对缠绕的差分干涉图在时域上进行叠加以增强形变信号，有效削弱了大气噪声的影响和避免了在复杂区域进行相位解缠可能带来的误差影响。在大范围的地表形变探测中可以实现对地表形变的位置进行快速的识别，有效的降低了计算的复杂度。所以，可以高效可靠地探测大范围区域内局部形变信号，在滑坡等地质灾害领域具有极大的应用价值。
附图说明
[0039]
图1为本发明实施例对insar差分干涉图进行相位梯度堆叠的技术流程图。
具体实施方式
[0040]
下面结合实施例，对本发明作进一步描述。以下实施例仅用于更清楚的说明本发明，不能因此限制本发明的保护范围。
[0041]
参见图1，本发明实施例提供一种基于insar相位梯度堆叠的广域地表形变探测方法，包含以下步骤：
[0042]
步骤一、获取覆盖研究区域的n景sar影像，选择合适的sar影像作为主影像，然后与其它sar影像进行配准并进行重采样，通过设置时空基线阈值，生成若干insar干涉图。
[0043]
本例中，以贵州省西部区域为研究区，范围为东经104度至106度，北纬24度至27度，范围约为250km*250km。研究区属于中国西南部高原山地，境内地势西高东低，大部分属于山地，地势起伏明显，植被茂密，大气噪声较为严重。数据源为覆盖贵州省西部区域的欧洲空间局哨兵1号数据，共获取2017年1月至2020年12月间的100幅左右sar卫星影像数据。
[0044]
步骤二、根据sar影像之间的时间基线和空间基线挑选干涉对，采用外部dem模拟地形相位，并从干涉图中减去地形相位得到差分干涉图，进行空间自适应滤波处理获取滤波后的差分干涉图；
[0045]
根据时间基线越长，相干性越差的原则出发，优选挑选相干性较好的insar差分干涉图作为下一步处理的数据源。针对哨兵数据，采用的干涉组合策略为时间基线小于60天和空间基线小于100米，实施例在数据处理过程中采用的多视比为10:2。利用sar卫星轨道信息从insar干涉图中去除平地相位和利用外部数字高程模型(如srtm)模拟地形相位以去除insar干涉相位中的地形相位，最终获得insar差分干涉图。为了平滑差分干涉图中的噪声信息，实施例采用自适应滤波的方法对差分干涉图进行滤波，在平滑噪声的同时也较好的保留了干涉相位的边缘信息。
[0046]
步骤三、对步骤二挑选并处理后的差分干涉图，在像元级上按照预设的步长分别进行东、西、南、北、东南、东北、西南、西北向的梯度计算，获得差分干涉图对应的各方向梯度图；
[0047]
一个像元周围有8个像元，每个像元代表一个方向，通过大量实验证明，8个方向结果合并在一起效果最好，因此实施例采用这种优选建议方案。
[0048]
利用步骤二处理得到的差分干涉图通过梯度计算获得差分干涉图所对应的各方向梯度图，
[0049]
设差分干涉相位对应的某一方向d梯度为表示如下：
[0050][0051]
其中，为sar卫星影像经过配准、差分等一系列处理得到的第i幅差分干涉图，m和n分别代表干涉图中像元所在的行列号，为差分干涉图对应的某一方向d的梯度图，s为步长；
[0052]
例如，其中差分干涉相位对应的东向梯度表示如下：
[0053][0054]
其中，为sar卫星影像经过配准、差分等一系列处理得到的第i幅差分干涉图，m和n分别代表干涉图中像元所在的行列号，为差分干涉图对应的东向的梯度图，s为步长。
[0055]
同理西、南、北、东南、东北、西南、西北七个方向的相位梯度分别用同理西、南、北、东南、东北、西南、西北七个方向的相位梯度分别用表示。具体实施时，可以根据经验值设置步长取值。
[0056]
步骤四、对获取的各个方向梯度图分别进行堆叠以增强形变信号并再次进行滤
波；分别将滤波后的各方向计算的相位梯度堆叠图进行合并，获得研究区域梯度堆叠结果。
[0057]
实施例采用的优选方案是，对获取的差分干涉图对应的各个方向梯度图分别进行时间上的堆叠以增强形变信号；将梯度堆叠结果采用空间滤波的方法进行滤波；将滤波后的各方向计算的梯度堆叠图的结果进行合并，获得研究区域差分干涉相位梯度堆叠图。
[0058]
对获取的差分干涉图对应的各个方向梯度图分别进行时间上的堆叠以增强形变信号时，实现方式为，
[0059]
对于不同时间基线的梯度图上的沿某一方向d，间隔预设步长的两个像元，其干涉相位梯度堆叠公式如下：
[0060][0061]
上式中，φd为获取的n个不同时间基线差分干涉图对应的沿某一方向间隔预设步长的梯度图叠加的结果，为第i幅差分干涉图对应的沿某一方向d获得的梯度，i＝1，2，3,
…
，n，n为干涉对总数；
[0062]
然后，对获得的梯度堆叠结果φd采用空间滤波的方法进行滤波，获得滤波后的结果用ψd表示。
[0063]
例如，对于不同时间基线的梯度图上的沿东方向，间隔一定预设步长的两个像元，其干涉相位梯度堆叠公式如下：
[0064][0065]
上式中，φe为获取的n个不同时间基线差分干涉图对应的沿东方向间隔预设步长的梯度图叠加的结果，为第i幅差分干涉图对应的沿东方向获得的梯度，i＝1，2，3,
…
，n，n为干涉对总数。
[0066]
然后，对获得的梯度堆叠结果φe采用空间滤波的方法进行滤波，获得滤波后的结果用ψe表示，同理西、南、北、东南、东北、西南、西北七个方向的相位梯度堆叠的滤波结果分别用ψw、ψs、ψn、ψ
en
、ψ
es
、ψ
wn
、ψ
ws
表示。
[0067]
最后，对进行空间滤波的各方向梯度堆叠图进行合并以获取完整的研究区域差分干涉相位梯度堆叠信息，其公式如下：
[0068][0069]
其中，ψ为合并各个方向干涉相位梯度图的结果，8为总的方向个数。
[0070]
需要说明的是，本发明将各个方向的梯度堆叠结果采用均方根的方法来进行各个方向梯度堆叠结果合并。
[0071]
步骤五、对获取的研究区域梯度结果ψ映射到预设范围得到最终的研究区域梯度堆叠结果。
[0072]
考虑到实际实施时，可能涉及多个研究区域，对每个研究区域分别执行以下流程，即分别进行梯度堆叠处理最终获取的结果值范围可能存在差异，因此本发明提出可把区间结果采用映射的方式统一到同一个区间，比如[-5 5]区间上，便于后续的输出、展示等应用。映射采用的方法是线性变换，之前得到的梯度无量纲的，所以可以进行区间映射。
[0073]
在进行了相位梯度堆叠和滤波处理后，大气信号和随机噪声被有效的抑制。由于采用大量的干涉图进行相位梯度堆叠，使得局部地表形变信号明显的增强，因此在获得的差分干涉相位梯度堆叠图上，明显的局部信号即为采用上述方法获取的地表形变探测信号，得到insar差分干涉相位梯度堆叠图后，通过局部放大图可以明显的看出形变信号。因此，对获取的研究区域差分干涉相位梯度的结果值映射到一定的范围得到最终的研究区域差分干涉相位梯度堆叠结果，在相位梯度堆叠图上出现聚集的与周边对比明显的大梯度值即为探测所得的地表形变信息。
[0074]
至此，本发明基于insar差分干涉相位梯度堆叠进行大范围地表形变探测的方法被完成和实现。通过以上流程，无需现有技术中的常规后期处理，即相位解缠(一是计算耗时；二是可能由此引起一些误差)，就可以快速准确的探测地表大范围的微小形变区域的位置，对于支持地质灾害的实时自动普查有一定的作用。
[0075]
具体实施时，本发明技术方案提出的方法可由本领域技术人员采用计算机软件技术实现自动运行流程，实现方法的系统装置例如存储本发明技术方案相应计算机程序的计算机可读存储介质以及包括运行相应计算机程序的计算机设备，也应当在本发明的保护范围内。
[0076]
在一些可能的实施例中，提供一种基于insar相位梯度堆叠的广域地表形变探测系统，包括以下模块，
[0077]
第一模块，用于通过获取覆盖研究区域的多幅sar影像，选取主影像并将其他影像与主影像配准，通过设置时空基线阈值，进一步处理获得干涉图；
[0078]
第二模块，用于根据sar影像之间的时间基线和空间基线挑选干涉对，采用外部dem模拟地形相位，并从干涉图中减去地形相位得到差分干涉图，并进行空间自适应滤波处理获取滤波后的差分干涉图；
[0079]
第三模块，用于对第二模块挑选并处理后的差分干涉图，在像元级上按照预设的步长分别进行像元周围多个方向的梯度计算，获得差分干涉图对应的各方向梯度图；
[0080]
设差分干涉相位对应的某一方向d梯度为表示如下：
[0081][0082]
其中，为sar影像经过处理得到的第i幅差分干涉图，m和n分别代表干涉图中像元所在的行列号，为差分干涉图对应的某一方向d的梯度图，s为步长；
[0083]
第四模块，用于对获取的各个方向梯度图分别进行堆叠以增强形变信号并再次进行滤波；分别将滤波后的各方向计算的相位梯度堆叠图进行合并，获得研究区域梯度堆叠结果；
[0084]
第五模块，对获取的研究区域梯度结果映射到预设的范围得到最终的研究区域梯度堆叠结果。
[0085]
在一些可能的实施例中，提供一种基于insar相位梯度堆叠的广域地表形变探测系统，包括处理器和存储器，存储器用于存储程序指令，处理器用于调用存储器中的存储指令执行如上所述的一种基于insar相位梯度堆叠的广域地表形变探测方法。
[0086]
在一些可能的实施例中，提供一种基于insar相位梯度堆叠的广域地表形变探测系统，包括可读存储介质，所述可读存储介质上存储有计算机程序，所述计算机程序执行
时，实现如上所述的一种基于insar相位梯度堆叠的广域地表形变探测方法。
[0087]
本文中所描述的具体实施例仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。