基于张量分解的时序InSAR影像时间-空间维降维压缩方法

基于张量分解的时序insar影像时间-空间维降维压缩方法
技术领域
1.本发明涉及长时序、大范围时序insar地表形变监测领域，尤其涉及一种基于张量分解的时序insar影像时间-空间维降维压缩方法。


背景技术：

2.对于较短时间序列和小范围的形变监测，现有时序insar技术较为有效。为获取长时序、大范围连续的地表形变信息，卫星将朝着短重访周期、大宽度以及长运行周期的方向发展，每周进行一次全球系统形变监测将成为可能，意味着insar将面临海量数据的处理难题。在海量sar影像数据条件下，时序insar技术可以获取更大范围、更长时间、更精确的地表形变信息，甚至孕育出新的应用前景和机遇。一方面，海量sar数据可以提供更加丰富的信息，从而提高insar形变估计的精度，众多应用领域将受益于这种前所未有的数据连续性，例如，地质灾害的预警和预报；另一方面，在新卫星不断发射，卫星编队数量不断增加，海量sar影像数据不断累积更新的条件下，对数据处理效率提出了新的挑战。滑坡蠕变、地壳形变等地质灾害具有形变周期长、形变缓慢且隐蔽性强的特点，需进行长时序insar监测才能有效提取其形变特征。但在长时序监测分析中，insar将面临失相干效应加重、计算效率偏低的问题。因此，如何充分利用存档的海量sar影像数据，研究高效率、高精度的数据处理方法，满足长时序地表形变监测需求，成为行业用户关心的重点问题，也是目前insar领域所面临的关键问题和难题之一。
3.随着sar数据量的不断积累，如sentinel-1每天提供大约10tb的sar产品，并且每个场景的时间序列sar影像的大小都超过1tb。此外，全国范围的sar存档数据量已达到pb级。insar数据量的剧增不仅需要的巨大存储设备和计算资源，而且极大地增加了数据处理负担和复杂性。为此，近几年insar界开始着手解决数据处理中的计算效率问题。在这方面，一种态度是保持传统的处理技术并依赖于并行计算以减少计算等待时间；应对海量sar数据处理的第二种态度是将传统的先进算法迁移到insar技术方案上(如降维压缩方法)，以便有效的处理和利用海量sar数据。然而，现有的时序insar影像降维压缩只考虑时间维的压缩处理，未考虑空间维的信息冗余。
4.在高光谱领域对于海量数据的处理主要是通过数据降维压缩来减少数据处理量，且具有较为成熟的理论方法，如主成分分析方法、因子分析方法、分形理论、经验模式分解、小波分析、张量分解等，而在时序insar领域的应用较为稀缺。
5.综上，对于长时序、大范围的时序insar地表形变监测存在数据量大、计算效率低等特点。


技术实现要素：

6.本发明要解决的技术问题在于克服现有技术的不足，提供基于张量分解的时序insar影像时间-空间维降维压缩方法，本发明借鉴高光谱领域影像降维压缩方法，引入张量分解理论，同时实现时序insar影像在时间维和空间维的降维压缩处理，旨在减少时序
insar数据处理量，提高其处理效率。为长时序、大范围地表形变监测(如大范围地质灾害隐患识别)提供了新的数据处理方案，以解决长时序、大范围地表形变监测数据处理量大、计算效率低等缺点。
7.为解决上述技术问题，本发明采用技术方案的基本构思是：
8.基于张量分解的时序insar影像时间-空间维降维压缩方法，用于长时序、大范围时序insar地表形变监测分析，所述降维压缩方法包括以下步骤：
9.1)将原始时序insar影像进行分组。具体为：采用人工干预手段将精配准后的原始时序insar影像进行平均分组，在每组中进行差分干涉处理和phase-linking。
10.2)子空间划分：在每组中将时序insar影像划分为多个小的子空间，并获取每个小的子空间范围内的像素协方差矩阵。
11.3)建立三阶张量模型：根据较小子空间内各像素之间具有相似的统计特性，将每个小子空间范围内的像素协方差矩阵进行堆叠，构建协方差矩阵三阶张量模型。
12.4)时序insar影像降维压缩张量分解模型：采用张量分解中的tucker分解算法对上述构建的协方差矩阵三阶张量模型进行分解。将原始张量分解为一个低秩核心张量与对应模式下的因子矩阵的乘积，得到时序insar影像的低秩表示，通过因子矩阵对核张量进行多线性反向投影，即可得到近似的低秩时序insar影像(称为虚拟影像)，具体如下：
13.根据步骤3)将给定子空间大小内的像素协方差矩阵表示成三阶张量其中l1和l2表示时间维(影像数量)，l3表示空间维(子空间像素个数)。根据tucker分解原理，三阶张量可用核张量和对应模式的因子矩阵表示，公式如下：
[0014][0015]
其中，以及为因子矩阵，且ji≤li,(i＝1,2,3)。对公式(1)进行线性投影变换，可获得核张量的表达式如下：
[0016][0017]
其中，以及为投影矩阵，公式(2)中的核张量也称为压缩张量，由于压缩张量的大小远小于原始张量因此，通过该方式即可实现时序insar影像的降维压缩。在实际计算过程中，采用了张量分解中的tucker分解，也称高阶奇异值分解(hosvd)，具体算法流程见算法1如下。
[0018]
hosvd算法：
[0019]
输入：输入阶数为l1×
l2×…×in
的张量数据
[0020]
j1×
j2×…×jn
：每个因子矩阵包含的列数
[0021]
[0022][0023]
5)生成虚拟干涉图
[0024]
利用步骤4)得到的虚拟影像进行干涉处理，生成虚拟干涉图。
[0025]
6)利用步骤5)生成的虚拟干涉图进行时序处理，获取地表形变信息。
[0026]
采用上述技术方案后，本发明与现有技术相比具有以下有益效果。
[0027]
本发明与现有sar影像降维压缩相比，其优势在于：
[0028]
1)本发明将张量理论引入时序insar影像降维压缩中，根据较小子空间范围内像素之间具有相似的统计特性，可将子空间内的像素协方差矩阵堆叠成三阶张量形式，然后构建时序insar影像降维压缩张量分解模型，同时实现了时间维和空间维的降维压缩处理，为长时序、大范围地表形变监测(如大范围地质灾害隐患识别)提供了新的数据处理方案；
[0029]
2)本发明采用tucker分解算法实现时序insar影像在时间维和空间维的降维压缩处理，获得数据量较少的虚拟影像，再利用虚拟影像进行干涉处理生成虚拟干涉图，进而提取时序形变信息。经降维压缩处理后的时序insar数据量大大减少，时序insar数据处理效率得到极大提升，对于快速获取长时序、大范围区域地表形变信息具有重要意义。
[0030]
下面结合附图对本发明的具体实施方式作进一步详细的描述。
附图说明
[0031]
附图作为本技术的一部分，用来提供对本发明的进一步的理解，本发明的示意性实施例及其说明用于解释本发明，但不构成对本发明的不当限定。显然，下面描述中的附图仅仅是一些实施例，对于本领域普通技术人员来说，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。在附图中：
[0032]
图1为时序insar影像时间-空间维降维压缩方法流程图；
[0033]
图2为时序insar影像降维压缩前后时空基线分布图；
[0034]
图3为时序insar影像时间-空间维降维压缩原理图；
[0035]
图4为时序insar影像降维压缩前后干涉图对比；
[0036]
图5为时序insar影像降维压缩前地表形变速率图；
[0037]
图6为压缩子空间2
×
3下获取的地表形变速率图；
[0038]
图7为压缩子空间2
×
5下获取的地表形变速率图；
[0039]
图8为压缩子空间2
×
10下获取的地表形变速率图；
[0040]
图9为压缩子空间2
×
15下获取的地表形变速率图；
[0041]
图10为时序insar影像降维压缩前后形变速率相关性图。
[0042]
需要说明的是，这些附图和文字描述并不旨在以任何方式限制本发明的构思范围，而是通过参考特定实施例为本领域技术人员说明本发明的概念。
具体实施方式
[0043]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对实施例中的技术方案进行清楚、完整地描述，以下实施例用于说明本发明，但不用来限制本发明的范围。
[0044]
实施例
[0045]
为验证本发明的可行性和可靠性，选取2018年1月至2021年5月期间的共101景升轨sentinel-1a数据(轨道号为26)进行实验分析，其波长为5.6cm，分辨率为5m
×
20m，重访周期为12天。每个影像由3个子带组成，每个子带由多个burst组成。在时序insar数据处理过程中，利用外部高程数字模型(dem)数据作为参考，dem数据可提供参考地形和地理坐标系，以便去除地形相位的影响，本文选用由美国航空航天局(nasa)和国防部国家测绘局(nima)联合测量得到的srtm-1高程数据作为参考，其分辨率可达30m。另外，在处理sentinel-1a数据之前，需要获取与之对应的精密轨道数据(pod)，以便去除因轨道误差引起的系统性误差。
[0046]
基于张量分解的时序insar影像时间-空间维降维压缩方法，技术方案实施流程如附图1所示，其步骤如下：
[0047]
1)对全堆栈时序影像进行分组处理，每组3景影像压缩为1景。按照影像获取时间将101景升轨影像分为33组，其中前32组中每组包含3景影像，最后一组包含5景影像，选取每组中的第2景作为参考影像。图2展示了全堆栈sar影像数据和压缩后的影像时空基线分布情况。
[0048]
2)将原始时序insar影像的整个空间划分为多个小的子空间，然后在每个子空间中构建三阶张量模型，利用张量分解方法将原始张量分解为一个低秩核心张量与对应模式下的因子矩阵的乘积，得到时序insar影像的低秩表示，通过因子矩阵对核张量进行多线性反向投影，即可得到近似的低秩时序insar影像(称为虚拟影像)，如图3所示。
[0049]
3)图4为时序insar影像压缩前和压缩后的干涉图，其中(a)和(b)分别为干涉对20190911-20180403影像压缩前和影像压缩后的干涉图，(c)和(d)分别为干涉对20190911-20201001影像压缩前和影像压缩后的干涉图。从图可以看出，当同时进行时间维和空间维压缩时，压缩后的干涉相位信息比压缩前的干涉相位信息少，主要原因是在时间-空间维压缩过程中，降低了空间分辨率，导致干涉相位信息减少。
[0050]
4)图5-图9展示了时间-空间维影像压缩前和影像压缩后得到的los方向上地表形变速率，其中图5为影像压缩前的地表形变速率，图6-图9分别为选取子空间大小为2
×
3、2
×
5、2
×
10以及2
×
15时进行空间维压缩得到的地表形变速率。从图中可以看出在各种压缩子空间下均能有效识别出地表形变位置，且其形变空间位置分布与压缩前保持高度一致，表明本文的影像压缩方案具有可行性，但获取的相干目标点空间密度具有一定差异，压缩后的相干目标点密度明显降低，该特征在山区较为明显，特别是压缩子空间为2
×
10和2
×
15时，在植被茂密的山区几乎没有获取到相干目标点。主要原因是山区每个压缩子空间内的地物特征相似性不高，相位信息冗余少，导致空间为压缩后的像素相干性降低，从而难以提取到相干目标点。图5-图9中的(b)展示了不同压缩子空间下得到的形变速率统计直方图，从图中可以看出各种压缩情况下得到的形变速率值具有较高的相似性，绝大部分分布在
±
30mm/yr以内，少数点形变速率绝对值大于30mm/yr。
[0051]
5)选取型形变区进行可靠性验证分析。分别以不同压缩子空间下获取的形变速率为参考，采用距离最近原则，在一定距离范围内(50m)，从影像压缩前获取的形变速率点中搜索与之对应的同名相干目标点，并绘制出各同名相干目标点影像前后的相关性图，如图10所示，图中的(a)、(b)、(c)、(d)分别表示压缩子空间为2
×
3、2
×
5、2
×
10、2
×
15下与影像压缩前获取的形变速率相关性图。从图中可以看出，在压缩子空间为2
×
3和2
×
5时相关性
较高，均达到了0.9以上，均方差分别为5.25mm/yr和5.43mm/yr，如图10(a)和(b)。在压缩子空间为2
×
10和2
×
15时相关性较低，分别为0.54和0.59，均方差分别达到9.86mm/yr和10.93mm/yr。因此，从上述结果看，在压缩子空间为2
×
3和2
×
5时获得的形变速率与影像压缩前的形变速率值较为接近，具有较高可靠性。而在压缩子空间为2
×
10和2
×
15时获得的形变速率与影像压缩前获得的形变速率值存在一定差距，可靠性相对较差。
[0052]
以上所述仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专利的技术人员在不脱离本发明技术方案范围内，当可利用上述提示的技术内容作出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明方案的范围内。