CN113687353A - 基于同质像元时序相位矩阵分解的ds目标相位优化方法 - Google Patents

Abstract

本发明公开了基于同质像元时序相位矩阵分解的DS目标相位优化方法，具体方法如下：获取监测区域的时序SAR影像，形成差分干涉图，识别同质像元；针对分布式散射体进行相位优化；基于优化后的时序差分干涉图，选取便于观察的时序监测点，去除轨道误差、DEM残差、大气相位改正，然后进行三维相位解缠以得到监测点形变结果，完成监测区域的形变监测，本方法所需数据存储和计算量相对较小；具有更好表现；在地表形变监测方面具有较高的可靠性；在监测点数量上较原始干涉图的PS方法提升了257％，较加权平均方法提升了194％，能够更好地反演地表形变。

Description

基于同质像元时序相位矩阵分解的DS目标相位优化方法

技术领域

本发明涉及DS相位优化方法技术领域，具体涉及基于同质像元时序相位矩阵分解的DS目标相位优化方法。

背景技术

断层活动、地震等地壳构造活动，地下水开采导致的地面沉降、塌陷，矿山开采引起的地面开采沉陷，以及由于水土流失造成的滑坡、泥石流等地质灾害对人民生命财产安全和社会发展产生了重大威胁。快速准确地监测地表形变，对于地质灾害评估与防治至关重要。

目前地表形变监测技术主要包括以水准测量、导线测量以及GNSS(GlobalNavigationSatelliteSystem)测量为代表的传统技术以及近几十年发展起来的合成孔径雷达干涉测量(InterferometrySyntheticApertureRadar，InSAR)技术。其中，传统监测技术监测精度高、时间分辨率可控，但就大范围形变监测而言，工程实施复杂，耗费大量人力物力，不适宜大区域的长期重复监测，且均是基于离散点的观测技术。而InSAR技术是基于面域的对地观测技术，具有大范围、高精度、高分辨率、全天时全天候等优势，已经被广泛地用于城市地表及建构筑物变形、矿区开采沉陷、冰川移动、火山喷发、滑坡、泥石流等监测中。

1993年，Massonnet和Rossi首次利用D-InSAR技术获取由地震导致的地表形变。自此，如何利用InSAR技术实现地表形变监测成为研究热点。由于D-InSAR技术监测地表形变易受到大气和失相干等因素影响导致其监测精度较低。为克服上述问题，研究者们提出了基于永久散射体的PS-InSAR(PermanentScattersInSAR,PS-InSAR)技术，该技术利用受失相干影响较小的永久散射体(PS目标)进行形变监测，较好地解决时空失相干问题。然而，在非城市地区，如农田、荒漠、戈壁等，存在PS目标少，导致形变监测点密度低的缺点。而在这些区域往往存在大量的分布式散射体(DistributedScatters,DS)目标，因此，结合DS目标信息的InSAR监测技术，即DS-InSAR技术，应运而生。然而，DS目标在长时、空基线或者植被覆盖区域等情况下，易受到噪声、失相干的影响，进而导致其形变监测效果不理想。因此，如何优化DS目标相位，提高其相位质量，是DS-InSAR技术关键。目前，针对DS目标的相位优化方法大多是通过像元协方差矩阵分解实现的，这些方法需要构建像元协方差矩阵，像元协方差矩阵的数据存储和计算量较大，此外，部分方法通过对同质像元加权平均来优化DS目标相位，这些方法DS目标相位优化反演地表形变的效果较差。

发明内容

为全面解决上述问题，尤其是针对现有技术所存在的不足，本发明提供了基于同质像元时序相位矩阵分解的DS目标相位优化方法能够全面解决上述问题。

为实现上述目的，本发明采用以下技术手段：

基于同质像元时序相位矩阵分解的DS目标相位优化方法，具体方法如下：

S1、获取监测区域的时序SAR影像，形成差分干涉图，识别同质像元；

S2、针对分布式散射体(DS目标)进行相位优化，以步骤S1识别的DS目标为基础，构建DS目标的同质像元时序相位矩阵，采用奇异值分解、SHP(statisticallyhomogeneouspixels)加权平均进而完成相位优化；

S3、基于优化后的时序差分干涉图，选取便于观察的时序监测点，去除轨道误差、DEM残差、大气相位改正，然后进行三维相位解缠以得到监测点形变结果，完成监测区域的形变监测。

优选的，步骤S1的具体处理步骤为：

步骤一：时序SAR影像获取；

步骤二：选取主影像；

步骤三：SAR数据的定标和配准；

步骤四：差分干涉SAR数据处理；

步骤五：同质像元识别。

优选的，所述差分干涉影像数据集后进行滤波处理或不进行滤波处理。

优选的，所述同质像元识别的结果通过设置阈值能够提取永久散射体(PS)与分布式散射体(DS目标)。

优选的，所述同质像元识别的算法包括BWS假设检验、双样本t假设检验、广义似然比假设检验、FaSHPS检验算法和HTCI算法。

优选的，步骤S2的具体方法为：

步骤一：构建同质像元时序相位矩阵

在同质像元时序相位矩阵中，每一行代表同一像元的时序相位，每一列代表同一干涉图同质像元相位，即其行数对应待优化像元的同质像元个数，列数对应时序干涉图幅数，构建同质像元时序相位矩阵，只需将窗口内的同质像元时序相位按行依次排序，构建方式如下：

式中，α_n，t代表窗口内第n个同质像元在第t幅干涉图中的相位值，α_ref，t为参考像元在第n幅影像上的相位值；

步骤二：针对同质像元时序相位矩阵的奇异值分解

在构建好同质像元时序相位矩阵D后，对该矩阵进行奇异值分解，公式如下：

式中，r为矩阵D的非零奇异值个数，得到r≤min(n，t)；S_i、U_i和V_i分别为奇异值矩阵、左右奇异向量，基于主成分分析方法和同一窗口内同质像元差分干涉相位基本一致原则，认为第一主成分，即最大的奇异值S₁及其对应的左右奇异向量U₁和V₁，可重构去噪后的干涉相位。至此，获得该窗口降噪后的矩阵D′，如下式：

α′_ref＝[α′_ref，1 α′_ref，2 … α′_ref，t]

对所有同质像元时序相位矩阵D进行上述过程，并提取每个D′中的参考像元相位α′_ref，进而获得分解后的干涉图；

步骤三：对分解后的干涉图进行加权平均获得优化干涉相位

对分解后的干涉图进行相干性估计，并利用相干性对DS目标同质像元进行加权平均得到最终的DS目标优化干涉相位，如下式：

式中，

表示参考像元的相位值，

代表周围同质像元点的相位值，Coh_i为估计得到的相干性。

本发明的有益效果：

1、本方法所需数据存储和计算量相对较小；

2、无论在模拟数据中还是在真实数据中，本发明较对比DS优化方法(加权平均方法)都具有更好表现；

3、本发明与加权平均方法、原始PS方法对地表形变监测的结果整体一致，因而本发明在地表形变监测方面具有较高的可靠性；

4、在监测点数量上本发明较基于原始干涉图的PS方法提升了257％，较加权平均方法提升了194％，能够更好地反演地表形变。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明基于同质像元时序相位矩阵分解的DS-InSAR形变监测流程图；

图2为本发明DS目标优化过程示意图；

图3为本发明模拟数据原始信息图；

图4为本发明模拟数据优化相位对比图；

图5为本发明模拟数据优化效果对比图；

图6为本发明真实数据优化相位对比图；

图7为本发明真实数据相位相干性对比图；

图8为本发明真实数据形变监测结果对比图。

具体实施方式

下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

在本发明的描述中，需要说明的是，术语“中心”、“上”、“下”、“左”、“右”、“竖直”、“水平”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性。

在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。

如图1所示，本发明实施例的基于同质像元时序相位矩阵分解的DS-InSAR地表形变监测方法，包括以下步骤：

S1，获取监测区域的时序SAR影像，形成差分干涉图，识别同质像元。

该步骤可分为以下几个主要处理步骤：时序SAR影像获取、选取主影像、SAR数据的定标和配准、差分干涉SAR数据处理和同质像元识别。在获取差分干涉影像数据集后可根据需求决定是否对其进行滤波处理；在同质像元识别结果基础上，通过设置相关阈值提取永久散射体(PS)与分布式散射体(DS目标)。所述的同质像元识别算法主要包括BWS假设检验、双样本t假设检验、广义似然比假设检验、FaSHPS检验算法和HTCI算法等。

S2，针对分布式散射体(DS目标)进行相位优化。

以步骤S1识别的DS目标为基础，构建DS目标的同质像元时序相位矩阵，采用奇异值分解、SHP(statisticallyhomogeneouspixels)加权平均进而完成相位优化，如图2所示。该过程进一步描述如下：

步骤一：构建同质像元时序相位矩阵

在同质像元时序相位矩阵(图2中D矩阵)中，每一行代表同一像元的时序相位(不同干涉图中的相位)，每一列代表同一干涉图同质像元相位，即其行数对应待优化像元的同质像元个数，列数对应时序干涉图幅数。构建同质像元时序相位矩阵，只需将窗口内的同质像元时序相位按行依次排序，构建方式如下：

式中，α_n，t代表窗口内第n个同质像元在第t幅干涉图中的相位值，α_ref，t为参考像元在第n幅影像上的相位值。

步骤二：针对同质像元时序相位矩阵的奇异值分解

在构建好同质像元时序相位矩阵D后，对该矩阵进行奇异值分解，公式如下：

式中，r为矩阵D的非零奇异值个数，得到r≤min(n，t)；S_i、U_i和V_i分别为奇异值矩阵、左右奇异向量。基于主成分分析方法和同一窗口内同质像元差分干涉相位基本一致原则，认为第一主成分，即最大的奇异值S₁及其对应的左右奇异向量U₁和V₁，可重构去噪后的干涉相位。至此，获得该窗口降噪后的矩阵D′，如下式：

α′_ref＝[α′_ref，1 α′_ref，2 … α′_ref，t]

对所有同质像元时序相位矩阵D进行上述过程，并提取每个D′中的参考像元相位α′_ref，进而获得分解后的干涉图_。

步骤三：对分解后的干涉图进行加权平均获得优化干涉相位

对分解后的干涉图进行相干性估计，并利用相干性对DS目标同质像元进行加权平均得到最终的DS目标优化干涉相位，如下式：

式中，

表示参考像元的相位值，

代表周围同质像元点的相位值，Coh_i为估计得到的相干性。

S3，基于优化后的时序差分干涉图，选取便于观察的时序监测点，去除轨道误差、DEM残差、大气相位改正，然后进行三维相位解缠以得到监测点形变结果，完成监测区域的形变监测。

模拟数据实施例：

本发明提出方法针对分布式散射体，故模拟数据只模拟DS点目标，使用matlab自带peaks函数制作无噪声情况下的相位信息β，使用图3中原始图像模拟无噪声情况下的相干性γ和强度信息R，为尽量接近真实情况，对模拟数据添加噪声信息，本实施例共模拟27景SLC数据和27对干涉图(包含自身干涉图)。详细流程如下：

步骤一：计算无噪声情况下的相位信息β，相干性γ和强度信息R

使用matlab自带peaks函数制作无噪声情况下的相位信息β，并进行相位缠绕。在本实施例中所使用像元个数为463*463，故使用peaks(463)生成相位信息，利用取余函数进行相位缠绕，得到模拟相位图如图3(b)所示。

根据输入图像的灰度信息计算模拟相干性γ。本实施例中，首先，根据图像灰度值计算参考相干性γ_ref，将灰度值大于100的相干性设为0.8，将灰度值小于100的相干性设为0.4。为了更接近真实情况，模拟的干涉图相干性随着影像时间间隔的增加而减少，减少模式如下：

式中，γ_ref为参考相干性，γ_n为第n幅SLC图像相干性，N代表模拟影像的总数。

根据输入图像的灰度信息计算模拟强度图R。在本实施例中，将图像灰度值按

进行归一化，再将归一化值进行平方，得到参考强度R_ref。为了更接近真实情况，本实施例模拟强度随时间变化，不失一般性，强度图模拟为周期性变化，模式如下：

式中，R_ref为参考强度图，R_n为第n幅SLC图像强度图，N代表模拟影像的总数。

步骤二：计算模拟干涉图协方差矩阵CV

步骤三：对CV矩阵进行Cholesky分解得到L矩阵

步骤四：对主副影像添加噪声。首先生成两个服从零均值复圆高斯分布的独立随机复数变量x₁和x₂，利用L矩阵计算模拟干涉图对应的单视复数影像z₁和z₂：

所得到干涉图的噪声水平与给定相干性相关。

步骤五：重复步骤三和步骤四，生成时序SLC影像和对应时序干涉图集。在本实施例中共生成27景SLC影像和对应时序干涉图集，如图3所示(b)、(c)和(d)从左至右分别为第一幅模拟干涉图所对应的真实相位、相干性及强度。

步骤六：模拟数据进行同质像元识别以及相位优化。本实施例中，同质像元识别采用的是HTCI算法，设置同质像元检测窗口为15*15，并设置同质像元阈值为20，同质像元个数小于阈值的为永久散射体目标，否则为分布式散射体目标。相位优化操作流程如S2所述，得到最终相位优化结果。

为进一步地说明本发明方法的优化提升效果，本实施例中，又采用加权平均方法(S2中的步骤三)对模拟数据进行相位优化。如图4所示，依次为第13幅的真实相位，模拟相位，加权平均优化结果和本发明方法优化结果。可以看出本发明方法优化结果更接近真实相位，能够更好的抑制噪声。更进一步地计算了优化相位平均中误差，并绘制了第460行第180-300列的相位剖面图以及(250，200)点的时序相位图。如图5所示，其中(a)、(b)分别为本发明方法和加权平均优化相位平均中误差，(c)、(d)分别为相位剖面图和点时序相位图。从图中可以看出，无论在高相干区域还是在低相干区域，本发明方法优化结果均更加接近真值，同时计算了总优化相位平均中误差分别为0.7622和0.3639，本发明方法较加权平均方法提升了52.26％。

真实数据实施例：

为验证发明方法在真实数据中的有效性，以徐州大泉街道周边地区为实例，采用74景Sentinel-1A数据(2017-11-04到2020-05-16)进行实验验证。具体实施方案如下：

步骤一：利用时序SAR影像形成差分干涉图，并识别同质像元。

该步骤可分为以下几个主要处理流程：时序SAR影像获取、选取主影像、SAR数据的定标和配准、差分干涉SAR数据处理和同质像元识别。在获取差分干涉影像数据集后可根据需求决定是否对其进行滤波处理，在本实施例中对干涉影像进行了Goldstein滤波处理；采用FaSHP算法对时序SAR影像进行识别同质像元，其中监测窗口设置为15*15，阈值设置为20，，同质像元个数小于阈值的为永久散射体目标(PS目标)，否则为分布式散射体目标(DS目标)。

步骤二：DS目标相位优化。

以步骤一识别的DS目标为基础，构建DS目标的同质像元时序相位矩阵，采用奇异值分解、SHP加权平均进而完成相位优化。详细流程见步骤S2.

步骤三：基于优化干涉图进行形变监测

基于优化后的时序差分干涉图，选取便于观察的时序监测点，去除轨道误差、DEM残差、大气相位改正，然后进行三维相位解缠以得到监测点形变结果，完成监测区域的形变监测。

为对比本发明方法与原始PS方法以及加权平均方法的监测效果，本实施例中，分别采用原始PS方法、加权平均方法和本发明提出的优化方法对SAR时序影像进行处理。

如图6所示，(a)、(b)、(c)图分别为第34幅原始差分干涉相位图、加权平均优化相位图以及本发明方法优化相位图，可以看出在高相干地区，两种方法均能取得不错的优化效果，本发明方法优于加权平均算法；而在低相干地区，加权平均表现不佳，而本发明方法仍能取得理想效果。进一步，计算了各个干涉图的相位相干性γ_ph，具体计算公式如下：

其中，

为像元相位值；M为窗口大小，这里为7*7。如图7中(a)、(b)、(c)图分别为原始相位相干性、加权平均优化相位相干性以及本发明方法优化相位相干性。从图中可以明显看出，经过本发明方法优化后的相位相干性明显优于加权平均方法和原始干涉图，同时，计算原始干涉图、加权平均干涉图和本发明方法干涉图的总平均相位相干性分别为0.2541，0.3305和0.6416，本发明方法较原始干涉结果和加权平均结果分别提升了153％和94％，说明本发明的优化效果提升明显。

进一步，本实施例中，分别采用基于原始干涉图、加权平均和本发明优化方法的干涉图，统一采用TPC筛选高质量点目标。阈值设定为0.9，对应相位标准偏差约为25°，分别基于原始和优化后的干涉图及选取的高质量像元，使用StaMPS做时序分析，得到结果如图8所示。图8中(a)、(b)、(c)图分别为原始相位监测结果图、加权平均监测结果图以及本发明方法监测结果图。从图中可以看出，三种方法监测的整体沉降趋势一致，侧面说明本发明方面监测的可靠性。三种监测方法的监测点分别为284186、344856和1014583，本发明方法较原始干涉结果和加权平均结果分别提升了257％和194％，监测点密度有了较大的提升，能够更好的反应实例区的地表形变细节_。

无论是在模拟数据还是在真实SAR影像数据中，本发明方法均表现出来较大的优势，相位优化结果更加接近真值，监测结果可靠，监测点密度高，能够更加有效、描述地反演地表形变。

本发明所作的举例，而并非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无需也无法对所有的实施方式予以穷举，而由此所引申出的显而易见的变化或变动仍处于本发明的保护范围中。

Claims (6)

1.基于同质像元时序相位矩阵分解的DS目标相位优化方法，其特征在于，具体方法如下：

S1、获取监测区域的时序SAR影像，形成差分干涉图，识别同质像元；

S2、针对分布式散射体即DS目标进行相位优化，以步骤S1识别的DS目标为基础，构建DS目标的同质像元时序相位矩阵，采用奇异值分解、SHP加权平均进而完成相位优化；

S3、基于优化后的时序差分干涉图，选取便于观察的时序监测点，去除轨道误差、DEM残差、大气相位改正，然后进行三维相位解缠以得到监测点形变结果，完成监测区域的形变监测。

2.如权利要求1所述的基于同质像元时序相位矩阵分解的DS目标相位优化方法，其特征在于，所述步骤S1的具体处理步骤为：

步骤一：时序SAR影像获取；

步骤二：选取主影像；

步骤三：SAR数据的定标和配准；

步骤四：差分干涉SAR数据处理；

步骤五：同质像元识别。

3.如权利要求2所述的基于同质像元时序相位矩阵分解的DS目标相位优化方法，其特征在于，步骤S1中，获取差分干涉影像数据集并进行滤波处理。

4.如权利要求2所述的基于同质像元时序相位矩阵分解的DS目标相位优化方法，其特征在于，步骤S1中，同质像元识别的结果通过设置阈值提取永久散射体与分布式散射体。

5.如权利要求2所述的基于同质像元时序相位矩阵分解的DS目标相位优化方法，其特征在于，步骤S1中，同质像元识别的算法包括BWS假设检验、双样本t假设检验、广义似然比假设检验、FaSHPS检验算法和HTCI算法。

6.如权利要求1所述的基于同质像元时序相位矩阵分解的DS目标相位优化方法，其特征在于，步骤S2的具体方法为：

步骤一：构建同质像元时序相位矩阵D，构建方式如下：

式中，α_n，t代表窗口内第n个同质像元在第t幅干涉图中的相位值，α_ref，t为参考像元在第n幅影像上的相位值；

步骤二：针对同质像元时序相位矩阵的奇异值分解，公式如下：

式中，r为矩阵D的非零奇异值个数，得到r≤min(n，t)；S_i为奇异值矩阵，U_i和V_i分别为左右奇异向量，至此，获得该窗口降噪后的矩阵D′，如下式：

α′_ref＝[α′_ref，1 α′_ref，2 … α′_ref，t]

对所有同质像元时序相位矩阵D进行上述过程，并提取每个D′中的参考像元相位α′_ref，进而获得分解后的干涉图；

步骤三：对分解后的干涉图进行加权平均获得优化干涉相位；

对分解后的干涉图进行相干性估计，并利用相干性对DS目标同质像元进行加权平均得到最终的DS目标优化干涉相位，如下式：

式中，

表示参考像元的相位值，

代表周围同质像元点的相位值，Coh_i为估计得到的相干性。

Images