一种全散射体FS-InSAR方法及系统与流程

一种全散射体fs-insar方法及系统
技术领域
1.本发明涉及星载合成孔径雷达时间序列干涉技术领域，具体涉及一种全散射体fs-insar方法及系统。


背景技术：

2.时间序列insar技术包括ps类 insar技术和 ds类 insar技术，在地面沉降、采空区塌陷、地震、火山、滑坡以及重要基础设施安全等形变监测领域得到了广泛应用。然而，时间序列insar高精度地表形变监测仍然面临两方面挑战：一是时间高度不相关的对流层大气效应干扰，特别是小尺度大气湍流效应，其延迟相位甚至超过形变相位。二是高植被覆盖地区地表严重失相干，导致监测密度稀疏，无法获取足够的形变信息。
3.目前大气校正方法可分为两大类：一是利用外部大气数据建模估计大气延迟相位，如气象数值模型、gps数据、成像频谱仪数据等。这些方法对大中尺度的大气湍流效应具有一定的改正效果，但受天气条件或时空分辨率等因素干扰，难以改正小尺度大气湍流效应。
4.二是基于insar数据本身的大气校正，如stacking insar、干涉图叠加大气估计方法以及ps或ds方法。stacking insar方法和干涉图叠加大气估计方法均假设地表形变为线性，所以不适用于非线性形变地区。ps或ds方法假定大气相位在一定距离内空间相关（通常小于2 km），通过相邻点相位差分抑制大气效应影响。
5.然而当出现小尺度大气效应时，相隔几百米甚至几十米范围的大气都可能出现较大差异，通过相位差分方式难以去除大气影响。
6.地表失相干问题目前主要通过分布式目标（ds）insar技术解决。ds提取方法有两种，一是非参数统计方法，如ks（kolmogorov-smirnov）检验、ad（anderson-darling）检验、bws（baumgartner-weiβ-schindler）检验等，其场景适应性强，但是计算效率低；二是参数统计方法，如似然比检验、fashps（fast statistically homogeneous pixel selection）、htci（hypothesis test of confidence interval）等方法，虽然计算效率较高，但要求地表服从高斯分布，对于植被覆盖地区难以满足该假设条件。因此现有ds方法主要应用于裸地、草地、农田等中度相干目标，而对于高密度植被覆盖的低相干地区应用较少。
7.所以，现有技术存在以下问题：1、利用外部大气数据建模估计大气延迟相位时，容易受天气条件或时空分辨率等因素干扰，难以改正小尺度大气湍流效应。
8.2、基于insar数据本身的大气校正，采用stacking insar、干涉图叠加大气估计方法时，不适用于非线性形变地区；ps或ds方法在出现小尺度大气效应时，通过相位差分的方式难以抑制大气影响。
9.3、针对地表失相干问题，通过分布式目标insar技术解决时容易出现计算效率低以及对于植被覆盖地区难以满足假设条件的问题。


技术实现要素：

10.针对上述现有技术的不足，本技术提供一种全散射体fs-insar方法及系统。
11.第一方面本技术提出了一种全散射体fs-insar方法，包括以下步骤：选定监测区，基于第一干涉像对组合并结合所述监测区的时序sar影像生成差分干涉图；对所述差分干涉图进行相位解缠，得到解缠相位；根据所述解缠相位得到每个sar成像时刻对应的时序相位；根据双尺度时域低通滤波，抑制所述时序相位中的大气相位和噪声相位，得到时序低频相位；重新计算干涉像对组合，得到第二干涉像对组合，根据所述第二干涉像对组合对所述时序低频相位进行差分，得到差分低频相位；根据所述时序sar影像得到平均幅度，根据平均幅度提取所述监测区中非水体的全散射体；基于所述全散射体和所述差分低频相位得到地表形变信息。
12.在一些实施例中，上述一种全散射体fs-insar方法中，所述选定监测区，基于第一干涉像对组合并结合所述监测区的时序sar影像生成差分干涉图，包括：设置第一时间基线阈值和第一垂直基线阈值得到所述第一干涉像对组合，所述第一干涉像对组合为连通网络，所述第一时间基线阈值小于等于两个卫星重访周期。
13.在一些实施例中，上述一种全散射体fs-insar方法中，所述对所述差分干涉图进行相位解缠，得到解缠相位，包括：采用相位解缠方法对所述差分干涉图进行相位解缠，得到解缠相位；所述相位解缠方法包括最小费用流法、枝切法和区域生长法。
14.在一些实施例中，上述一种全散射体fs-insar方法中，所述根据所述解缠相位得到每个sar成像时刻对应的时序相位，包括：利用最小二乘方法计算每个sar成像时刻的所述解缠相位，得到对应的时序相位，所述时序相位可表示为：（1）其中，为形变相位，为dem高程误差相位，为轨道误差相位，为大气延迟相位，为噪声相位。
15.在一些实施例中，上述一种全散射体fs-insar方法中，在所述根据双尺度时域低通滤波，抑制所述时序相位中的大气相位和噪声相位，得到时序低频相位前，包括：确定双尺度时间窗口；所述双尺度时间窗口包括第一尺度时间窗口和第二尺度时间窗口，所述第一尺度时间窗口小于第二尺度时间窗口。
16.在一些实施例中，上述一种全散射体fs-insar方法中，所述根据双尺度时域低通滤波，抑制所述时序相位中的大气相位和噪声相位，得到时序低频相位，包括：采用第一尺度时间窗口并通过低通滤波器对所述时序相位进行低通滤波，得到第
一低频相位：（2）（3）其中，为时序相位，为高斯函数，为第一尺度时间窗口下输出的时刻的低通滤波相位，为n幅sar影像的成像日期，为第一尺度时间窗口；利用所述时序相位减去所述第一低频相位，得到残余相位，所述残余相位为：（4）采用第二尺度时间窗口并通过低通滤波器对所述残余相位再次进行低通滤波，得到第二低频相位：（5）（6）其中，为第二尺度时间窗口，为第二尺度时间窗口下输出的时刻的低通滤波相位；将第一低频相位和第二低频相位叠加，得到时序低频相位，时序低频相位为：（7）其中，为所述时序低频相位。
17.在一些实施例中，上述一种全散射体fs-insar方法中，所述重新计算干涉像对组合，得到第二干涉像对组合，根据所述第二干涉像对组合对所述时序低频相位进行差分，得到差分低频相位，包括：设置第二时间基线阈值和第二垂直基线阈值得到所述第二干涉像对组合，所述第二干涉像对组合的数目大于所述第一干涉像对组合的数目；根据所述第二干涉像对组合，对所述时序低频相位进行差分，得到差分低频相位：（8）其中，为主影像对应的时序低频相位，为辅影像对应的时序低频相位，,。
18.在一些实施例中，上述一种全散射体fs-insar方法中，所述根据所述时序sar影像
得到平均幅度，根据平均幅度提取所述监测区中非水体的全散射体，包括对所述时序sar影像进行相对归一化处理得到平均幅度；设置筛选阈值，若所述时序sar影像的平均幅度大于所述筛选阈值，则提取对应的像元为全散射体，所述全散射体表示为：（9）其中，为时序sar影像的平均幅度，为筛选阈值。
19.在一些实施例中，上述一种全散射体fs-insar方法中，所述基于所述全散射体和所述差分低频相位得到地表形变信息，包括：根据所述全散射体和所述差分低频相位，提取所述全散射体对应的差分低频相位；利用常规时序insar方法，处理所述全散射体的差分低频相位，得到地表形变信息。
20.第二方面本技术提出了一种全散射体fs-insar系统，包括差分干涉图生成模块、相位解缠模块、时序相位生成模块、双尺度时域低通滤波模块、差分低频相位生成模块、全散射体提取模块和地表形变获取模块；所述差分干涉图生成模块，用于选定监测区，基于第一干涉像对组合并结合所述监测区的时序sar影像生成差分干涉图；所述相位解缠模块，用于对所述差分干涉图进行相位解缠，得到解缠相位；所述时序相位生成模块，用于根据所述解缠相位得到每个sar成像时刻对应的时序相位；所述双尺度时域低通滤波模块，用于根据双尺度时域低通滤波，抑制所述时序相位中的大气相位和噪声相位，得到时序低频相位；所述差分低频相位生成模块，用于重新计算干涉像对组合，得到第二干涉像对组合，根据所述第二干涉像对组合对所述时序低频相位进行差分，得到差分低频相位；所述全散射体提取模块，用于根据所述时序sar影像得到平均幅度，根据平均幅度提取所述监测区中非水体的全散射体；所述地表形变获取模块，用于基于所述全散射体和所述差分低频相位得到地表形变信息。
21.本发明的有益效果：对于本领域技术人员而言，目前大家公认的是，利用低通滤波器对点目标的残余形变（即从原始差分干涉相位中去除主形变和高程误差相位之后的相位）进行处理，不但需要时间低通滤波，而且还需要空间低通滤波处理，从而达到抑制大气相位和抑制除噪声相位得效果；本方案仅利用时间低通滤波处理就能实现从差分干涉相位中直接抑制高程误差相位、轨道误差相位、大气相位和噪声相位，且不需要空间低通滤波处理，滤波后的相位质量大幅提高，而且处理对象是非水体的全部像元，无需选择点目标或分布式目标，并且解决了现有技术中难以改正小尺度大气湍流效应；ps或ds方法在出现小尺度大气效应时，通过相位差分的方式难以去除大气影响；通过分布式目标insar技术解决时容易出现计算效率
低以及对于植被覆盖地区难以满足假设条件的问题。
附图说明
22.图1为本发明的总体流程图。
23.图2为第一干涉像对组合的时空基线分布图。
24.图3为第二干涉像对组合的时空基线分布图。
25.图4为差分干涉相位图。
26.图5为双尺度时域低通滤波后的差分高频相位图。
27.图6为双尺度时域低通滤波后的差分低频相位图。
28.图7为平均幅度图。
29.图8为点目标沉降速率图。
30.图9为全散射体沉降速率图。
31.图10为本发明的系统原理框图。
具体实施方式
32.本技术提出一种全散射体fs-insar方法及系统，解决了现有技术中难以改正小尺度大气湍流效应；ps或ds方法在出现小尺度大气效应时，通过相位差分的方式难以抑制大气影响；通过分布式目标insar技术解决时容易出现计算效率低以及对于植被覆盖地区难以满足假设条件的问题，仅通过双尺度时域低通滤波便能实现大幅提高相位质量以及提升地表形变监测密度的技术效果。
33.下面结合附图对本技术作进一步描述，以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本技术的保护范围。
34.第一方面本技术提出了一种全散射体fs-insar方法，如图1所示，包括以下步骤：s100选定监测区，基于第一干涉像对组合并结合所述监测区的时序sar影像生成差分干涉图；以天津市20201002
ꢀ‑
20211009获取的sentinel-1影像为例，共有30景影像，其时间基线如下表所示：
为了便于后续相位解缠，在保证差分干涉图时空基线连通的前提下，对于第一干涉像对组合，应尽可能减小时间基线和垂直基线，所以设置第一时间基线阈值和第一垂直基线阈值得到所述第一干涉像对组合，其中，设置的第一时间基线阈值小于等于两个卫星重访周期，本实施例中设置≤24天，≤200米，从而得到53个第一干涉像对组合，即小基线干涉像对组合，如图2所示；特殊情况下，当相邻两景sar影像时间间隔较大时，第一时间基线阈值也可以大于两个重访周期，本实施例中当相邻两景sar影像时间间隔为36天时，此时就需要大于两个重访周期；值得一提的是，对于不同类型的卫星，其重访周期不同，所以需要根据卫星类型具体确定对应的重访周期。
35.s200：对所述差分干涉图进行相位解缠，得到解缠相位；采用相位解缠方法对所述差分干涉图进行相位解缠，得到解缠相位；值得说明的是，相位解缠方法包括但不限于最小费用流法、枝切法和区域生长法，技术人员可以根据实际应用适应性选择相位解缠方法。
36.其中，本实施例中采用最小费用流相位解缠方法对所述差分干涉图进行相位解缠，得到解缠相位。
37.s300：根据所述解缠相位得到每个sar成像时刻对应的时序相位；对所述解缠相位利用最小二乘方法计算每个sar成像时刻的相位，得到对应的时序相位，所述时序相位可表示为：
（1）其中，为形变相位，为dem高程误差相位，为轨道误差相位，为大气延迟相位，为噪声相位。
38.在所述根据双尺度时域低通滤波，得到时序低频相位前，需要确定双尺度时间窗口，确定双尺度时间窗口的方法如下：通过双尺度时域低通滤波模拟实验，确定双尺度时间窗口；其中，双尺度时域低通滤波模拟实验为：采用模拟数据验证双尺度时域低通滤波分离形变相位的有效性：设时间采样间隔为12天，形变发生时段为730天，模拟振幅为10rad、周期为365天的正弦形变相位。参考多个无形变区域的真实sar时序差分干涉相位，以均值为0rad、标准差为8rad的高斯噪声作为高频信号，与形变相位叠加，得到模拟的时序差分干涉相位。
39.首先采用第一尺度时间窗口进行低通滤波。设时间窗口大小为1~365天，分别对时序相位滤波，并计算每次滤波后的低频相位与模拟形变相位的差值标准差和均值；结果显示，当时间窗口设为34天时，低频相位与模拟形变相位误差的标准差最小为1.85rad，均值为-1.34rad。
40.以34天的滤波结果为基础，将模拟形变相位减去一次低频相位，得到残余形变相位；将模拟时序差分相位减去一次低频相位，得到残余差分相位，仍设时间窗口大小为1~365天，分别对该残余差分相位滤波，并计算二次低频相位与残余形变相位的差值标准差和均值；结果显示，当时间窗口为77天时，二次滤波的误差标准差最小为3.85rad，均值为-1.35rad。
41.将两次滤波的时序低频相位叠加，得到总的时序低频相位，总的时序低频相位与模拟形变相位的差值标准差为1.74rad，均值为-1.32rad。与第一次滤波相比，差值标准差降低了5.43%，表明双尺度时域低通滤波比单次滤波能更好的拟合真实形变。
42.上述模拟实验表明，在双尺度时域低通滤波中，第一次滤波的时间窗口应较小，可设为30-40天；第二次滤波的时间窗口应较大，可设为60-80天。
43.所以得到第一尺度时间窗口的范围为30-40天，第二尺度时间窗口的范围为60-80天。
44.s400：根据双尺度时域低通滤波，抑制所述时序相位中的大气相位和噪声相位，得到时序低频相位；在时间序列上，形变相位表现为低频信号；轨道误差相位、大气延迟相位、噪声相位均为高频信号；高程误差相位与垂直基线有关，表现为高频信号；因此，根据不同相位成分的时域频谱特征，利用低通滤波器，可将形变相位与其他相位分离，从而有效抑制大气延迟相位和噪声相位的影响；其中低通滤波器可以是高斯滤波器，也可以是其他低通滤波器，本实施例以高斯滤波器为例进行详细说明：采用第一尺度时间窗口并通过低通滤波器对所述时序相位进行低通滤波，得到第
一低频相位：（2）（3）其中，为时序相位，为高斯函数，为第一尺度时间窗口下输出的时刻的低通滤波相位，为n幅sar影像的成像日期，为第一尺度时间窗口；由于高斯滤波采用加权平均方式实现，距离中心点越近的相位权重越大，距离中心点越远的相位权重越小。因此滤波后的一次低频相位虽然以形变相位为主，还含有少量的高程误差相位、轨道误差相位、大气相位以及噪声相位；所以，利用所述时序相位减去所述第一低频相位，得到残余相位，所述残余相位为：（4）残余相位中不仅包含高程误差相位、轨道误差相位、大气相位以及噪声相位等高频相位，而且还有少量的形变相位。为了提取该部分形变相位，所以采用第二尺度时间窗口并通过低通滤波器对所述残余相位再次进行低通滤波，得到第二低频相位：（5）（6）其中，为第二尺度时间窗口，为第二尺度时间窗口下输出的时刻的低通滤波相位；将第一低频相位和第二低频相位叠加，得到时序低频相位，时序低频相位为：（7）其中，为所述时序低频相位。
45.s500：重新计算干涉像对组合，得到第二干涉像对组合，根据所述第二干涉像对组合对所述时序低频相位进行差分，得到差分低频相位；以天津市20201002
ꢀ‑
20211009获取的sentinel-1影像为例，共有30景影像，其时间基线表不再赘述，设置第二时间基线阈值和第二垂直基线阈值得到所述第二干涉像对组合，所述第二干涉像对组合的数目大于所述第一干涉像对组合的数目；其中，≤100天，≤200米，从而得到194个第二干涉像对组合，如图3所示；
根据所述第二干涉像对组合，对所述时序低频相位进行差分，得到差分低频相位：（8）其中，为主影像对应的时序低频相位，为辅影像对应的时序低频相位，,。
46.s600：根据所述时序sar影像得到平均幅度，根据平均幅度提取所述监测区中非水体的全散射体；对所述时序sar影像进行相对归一化处理（即每幅sar影像幅度除以其幅度均值）得到平均幅度；设置筛选阈值，若所述时序sar影像的平均幅度大于所述筛选阈值，则提取对应的像元为全散射体，设置筛选阈值是为了全散射体中不含有水体像元，所述全散射体表示为：（9）其中，为时序sar影像的平均幅度，为筛选阈值，筛选阈值可设为0.5，由于不同计算方法会得到不同的平均幅度，因此筛选阈值可根据具体情况而定。
47.需要说明的是，由于时序低频相位以形变相位为主，其他高频相位部分已基本去除，因此我们可以选择非水体像元作为全散射体（fullscatterers，fs）进行地表形变监测，与ps-insar或ds-insar方法相比，这将极大地提高监测密度，获取全域地表形变信息。
48.可选的，由于全散射体像元数量巨大，当选定的监测区范围较大时运算耗时较长，可以采用常规稳定点目标选择方法，如幅度离差法、幅度差分离差法、平均相干系数法等提取需要的监测对象。
49.s700：基于所述全散射体和所述差分低频相位得到地表形变信息。
50.根据所述全散射体和所述差分低频相位，提取所述全散射体对应的差分低频相位；利用常规时序insar方法，处理所述全散射体的差分低频相位，得到地表形变信息。
51.其中，具体处理过程包括全散射体目标三角网连接、求解相邻目标的相对线性形变速率和相对高差误差、求解每个目标的线性形变速率和高程误差、残余相位计算、非线性形变估计、累计形变信息获取等。
52.下面结合图4-图9说明本方案的核心效果：其中，图4为天津市sentinel-1影像的差分干涉相位，图5为天津市sentinel-1影像的双尺度时域低通滤波后的差分高频相位，图6为天津市sentinel-1影像的双尺度时域低通滤波后的差分低频相位；从中可以看出仅通过时间低通滤波处理便有效抑制了大气相位和噪声相位，相位质量得到大幅提高。
53.其中，图7为400*400范围内平均幅度，图8为点目标沉降速率，图9为全散射体沉降速率，结果对比：点目标数量为23782个，监测密度为：14.86%；全散射体数量为157685个，监测密度为：98.55%，是点目标的6.63倍，对比可知监测对象是非水体的全散射体（full scatterers，fs），与现有的点目标或分布式目标方法相比，监测密度得到了突破性提升。
54.第二方面本技术提出了一种全散射体fs-insar系统，如图10所示，包括差分干涉图生成模块、相位解缠模块、时序相位生成模块、双尺度时域低通滤波模块、差分低频相位生成模块、全散射体提取模块和地表形变获取模块；所述差分干涉图生成模块，用于选定监测区，基于第一干涉像对组合并结合所述监测区的时序sar影像生成差分干涉图；所述相位解缠模块，用于对所述差分干涉图进行相位解缠，得到解缠相位；所述时序相位生成模块，用于根据所述解缠相位得到每个sar成像时刻对应的时序相位；所述双尺度时域低通滤波模块，用于根据双尺度时域低通滤波，抑制所述时序相位中的大气相位和噪声相位，得到时序低频相位；所述差分低频相位生成模块，用于重新计算干涉像对组合，得到第二干涉像对组合，根据所述第二干涉像对组合对所述时序低频相位进行差分，得到差分低频相位；所述全散射体提取模块，用于根据所述时序sar影像得到平均幅度，根据平均幅度提取所述监测区中非水体的全散射体；所述地表形变获取模块，用于基于所述全散射体和所述差分低频相位得到地表形变信息。
55.以上仅是本发明优选的实施方式，需指出的是，对于本领域技术人员在不脱离本技术方案的前提下，作出的若干变形和改进的技术方案应同样视为落入本权利要求书要求保护的范围。