基于无人机的时序InSAR地表形变数据空洞处理方法

基于无人机的时序insar地表形变数据空洞处理方法
技术领域
1.本发明涉及insar地表形变监测技术领域，特别是一种基于无人机的时序insar地表形变数据空洞处理方法。


背景技术：

2.地表形变监测长期以来受到各研究领域的密切关注，尤其涉及自然灾害、地球物理、岩土工程、勘测工程等。地表形变对社会稳定性评估与经济财产风险评估、地形地貌变化研究与工程地址勘探研究等人类活动有着深远影响。长期处于高形变幅度与形变速率的地表结构会对其上的基础构筑物产生大量的应力与能量堆积，导致地面崩塌、滑坡、泥石流等灾害发生，进而破坏构筑物的稳定性，威胁人类活动安全。因此，对地表形变制定科学、规律、有效的监测方法对人类社会发展有重要意义。
3.近年来，合成孔径雷达干涉测量insar(interferometric synthetic aperture radar)技术多被用于时序地表形变监测工作。insar技术来源于合成孔径雷达sar(synthetic aperture radar)工作原理。sar数据的构成形式为复数对，一组复数对包含一对振幅和相位，用于记录sar天线测量地表的后向散射强度和微波传递周期。由于单幅sar数据格式中用以记录周期的相位值域限制，对地表三维结构的解析能力弱，难以成体系运用。在经历了双幅、多幅sar影像差分干涉技术发展后，基于1801年thomas young的“杨氏双狭缝光干涉实验”原理，insar技术被提出。当以固定时序对同一地表区域进行测量时，由于卫星轨道差异，造成sar天线具有邻点相位差。获得多幅sar数据后，以像对形式提取雷达影像数据像素中的相位差，获得干涉相位图。对干涉相位图进行进一步特征处理，可以得到包含地形起伏、地表粗糙度、地表纹理信息、地表形变在内的多种数据。
4.目前，主流时序insar技术分为永久散射体干涉测量方法ps-insar(persistent scatterers)和小基线集干涉测量方法sbas-insar(small baseline subset)。ps-insar以高相干性的永久散射体为分析目标，通过对分辨单元内的强烈散射体进行ps点定位，获得地表形变的点集合信息；sbas-insar以时间、视角基线控制手段减少去相关和高程误差影像，在对数据进行滤波、多视处理、相位解缠后，获得地表形变的面集合信息。
5.主流时序insar目前的精度范畴都可达毫米级，但存在的共同问题为：获得的地表形变数据均为不连续集合，存在数据空洞。原因在于相位跨度存在小于一个像元的情况，或相位分布由于大气干扰和dem混淆过密集。尽管目前有相关算法尝试从地表形变结果信息中对数据空洞部分进行预测性填补，但数据结果缺乏检验，难以应用于地表形变监测的真实数据采集工作中。
6.因此，目前需要一种可用于填补时序insar地表形变数据空洞的技术方法，以用于真实数据采集工作，提高地表形变监测数据的完整性与精确性。


技术实现要素：

7.本发明要解决的技术问题是针对上述现有技术的不足，而提供一种基于无人机的
时序insar地表形变数据空洞处理方法，该基于无人机的时序insar地表形变数据空洞处理方法能实现以高精度真实地表形变数据耦合时序insar地表形变计算结果，并依靠算法完成数据的拟合优度检验。
8.为解决上述技术问题，本发明采用的技术方案是：
9.一种基于无人机的时序insar地表形变数据空洞处理方法，包括如下步骤。
10.步骤1、获取前置时序insar地表形变数据集：设置insar地表形变数据孔洞处理的时序变更点，将位于时序变更点之前的insar地表形变数据称为前置时序insar地表形变数据，将位于时序变更点之后的insar地表形变数据称为后置时序insar地表形变数据；按照确定的insar数据范围、重访周期和时序变更点，在每个重访周期内均采集一次前置时序insar地表形变数据，从而得到包含空洞位置信息的前置时序insar地表形变数据集。
11.步骤2、获取无人机平台后置时序空洞填补数据集：无人机根据步骤1中的空洞位置信息设置航迹，并从步骤1的时序变更点开始，参照步骤1的重访周期和时序变更点，在每个重访周期的时序变更点开始时，均根据设置的航迹采集一次台后置时序空洞填补数据，从而得到无人机平台后置时序空洞填补数据集。
12.步骤3、获取后置时序insar地表形变数据集：根据步骤2中无人机平台后置时序空洞填补数据集的数据范围，从步骤1的时序变更点开始，参照步骤1的重访周期和时序变更点，在每个重访周期的时序变更点开始时，均采集一次后置时序insar地表形变数据，从而得到包含空洞位置信息的后置时序insar地表形变数据集。
13.步骤4、获取前置时序空洞填补预测集：对步骤1获取的前置时序insar地表形变数据集中的所有空洞均进行预测填补，从而形成前置时序空洞填补预测集。
14.步骤5、获取后置时序空洞填补数据集：使用步骤2得到的无人机平台后置时序空洞填补数据集，对步骤3得到的后置时序insar地表形变数据集中空洞进行融合与填补，将空洞填补后的后置时序insar地表形变数据集称为后置时序空洞填补数据集。
15.步骤6、拟合优度检验：将步骤5获取的后置时序空洞填补数据集，采用步骤4获取的前置时序空洞填补预测集进行拟合优度检验，当拟合优度检验指标满足设定阈值时，判断空洞填补完成；否则，重复步骤5至步骤6，直至拟合优度检验指标满足设定阈值要求。
16.步骤5中，后置时序空洞填补数据集fixed set的获取方法，包括如下步骤。
17.步骤51、坐标系统一：设步骤2获取的无人机平台后置时序空洞填补数据集为patch set，步骤3获取的后置时序insar地表形变数据集为raw set；接着，将patch set和raw set分别进行投影并转换至同一坐标系下。
18.步骤52、raw set边界提取：对raw set中的每个空洞均进行边界提取。
19.步骤53、patch set像元裁剪：对patch set进行像元裁剪处理。
20.步骤54、生成fixed set：将步骤52中边界提取后的raw set和步骤53中像元裁剪后的patch set，进行格栅数据融合，从而形成fixed set。
21.步骤54中，raw set和patch set进行格栅数据融合的方法，包括如下步骤。
22.步骤54-1、像元灰度匹配：根据像元灰度匹配的方法，对patch set中的像元执行嵌入及覆盖，对象为raw set中的经过边界提取的对应空洞。
23.步骤54-2、边界平滑：对嵌有像元的每个空洞边界进行平滑处理，从而形成fixed set。
24.步骤54-2中，每个空洞边缘平滑处理的方法，包括如下步骤。
25.步骤54-2a、计算点位权值：将位于每个空洞边界领域的每个像元分别计算点位权值；其中，点位权值m的计算公式为。
[0026][0027]
式中，w1和w2均为权重系数，设定值；a为当前像元对应的地表形变数据；为当前空洞边界领域中所有像元对应的地表形变数据的均值；b为当前像元至对应空洞边界的距离；表示对应空洞中心点至当前空洞边界的距离均值。
[0028]
步骤54-2b、平滑：将点位权值高于设定权值阈值的像元，采用高斯核函数进行当前像元的平滑计算，并将计算结果赋值给当前像元，实现当前像元的平滑。
[0029]
步骤1中，前置时序insar地表形变数据集获取的方法，包括如下步骤。
[0030]
步骤11、确定目标区域范围：设定需要收集的前置时序insar地表形变数据范围与地理位置。
[0031]
步骤12、确定重访周期：根据待填补地表所处地理位置，确定insar数据系统类型；根据insar数据系统所对应的卫星型号，确定重访周期。
[0032]
步骤13、时序变更点定位：根据连续测量的前置时序insar地表形变数据，定位时序变更点。
[0033]
步骤14、生产前置时序insar地表形变数据集：在每个重访周期，均使用对应的insar数据系统，采集一次前置时序insar地表形变数据，从而得到包含空洞位置信息的前置时序insar地表形变数据集。
[0034]
步骤14中，前置时序insar地表形变数据集中包括至少15幅时序间距均匀的前置时序insar地表形变数据；步骤12中，insar数据系统类型的确定方法为。
[0035]
a、当待填补地表所处地理位置中人工构筑物占地面积超过设定值或地势起伏度平缓时，insar数据系统使用ps-insar系统。
[0036]
b、当待填补地表所处地理位置中人工构筑物占地面积不超过设定值或地势起伏度陡峭时，insar数据系统使用sbas-insar。
[0037]
步骤2中，无人机中内置有雷达信号源和波段切换器。
[0038]
雷达信号源能产生指定波段的雷达信号源，雷达波段的覆盖范围包含波长范围0.1米至1米的空间波雷达。
[0039]
波段切换器能切换当前工作雷达信号源的波长范围。
[0040]
步骤4中，通过生成roc曲线和采用roc曲线进行插值预测的方法，对步骤1获取的前置时序insar地表形变数据集中的所有空洞均进行预测填补。
[0041]
步骤6中，拟合优度检验指标为可决系数r2。
[0042]
步骤6中，设步骤4获取的前置时序空洞填补预测集为validation set，则r2的计算公式为。
[0043][0044]
式中，rss为validation set与fixed set的残差平方和；tss为validation set与fixed set的总离差。
[0045]
本发明具有如下有益效果：
[0046]
1、通过特高频雷达包含的多波段雷达波长发射范围，对各类地物反射信息有着良好的勘察效果，并与多种insar数据采集平台的元数据格式的匹配性良好。
[0047]
2、通过无人机工作平台与地面信息处理平台的协同工作，解决雷达原始数据难以解译和判别可用性的问题。
[0048]
3、通过无人机平台的高机动性满足数据空洞填补位置随机性的需求，解决时序insar地表形变数据可用性不稳定的问题。
[0049]
4、通过以无人机采样获得的高精度真实地表形变数据耦合时序insar地表形变计算结果，解决地表形变数据空洞难以用insar数据算法训练生成可用性、可靠性较高的填补数据问题。
[0050]
5、通过前、后置insar时序数据配合无人机采样数据进行填补数据的拟合优度检验，解决地表形变数据空洞填补置信度控制的问题。
附图说明
[0051]
图1显示了本发明基于无人机的时序insar地表形变数据空洞处理方法的流程图。
[0052]
图2显示了本发明基于无人机的时序insar地表形变数据空洞处理工作系统的原理图。
[0053]
图3显示了用于实施本发明处理方法的insar原始数据截取。
[0054]
图4显示了本发明处理方法的实施对象截取，时序insar数据空洞。
[0055]
图5显示了实施本发明处理方法的成果截取，时序insar数据空洞填补效果。
具体实施方式
[0056]
下面结合附图和具体较佳实施方式对本发明作进一步详细的说明。
[0057]
本发明的描述中，需要理解的是，术语“左侧”、“右侧”、“上部”、“下部”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，“第一”、“第二”等并不表示零部件的重要程度，因此不能理解为对本发明的限制。本实施例中采用的具体尺寸只是为了举例说明技术方案，并不限制本发明的保护范围。
[0058]
本实施例以某滑坡、泥石流易发山区为研究对象，通过本发明实验，验证研究区的地表形变数据空洞填补工作的有效性。
[0059]
如图2所示，基于无人机的时序insar地表形变数据空洞处理工作系统，包括：
[0060]
多旋翼的无人机，无人机搭载有地表三维信息采集模块、数据存储模块、数据传输模块和地面信息处理平台。
[0061]
无人机平台搭载的地表三维信息采集模块包括：
[0062]
特高频雷达信号源，用以产生指定波段的雷达信号源，雷达波段覆盖范围包含波长范围0.1米至1米的空间波雷达，与目前主流的c波段、l波段insar数据采集信号有良好适配性。
[0063]
波段切换器，用以切换当前工作的雷达信号源波长范围。
[0064]
发射机，用以发射制定波段的雷达信号。
[0065]
接收机，用以接受经地表吸收和后向散射之后的发射雷达信号。
[0066]
信号转换机，用以将接收机接受到的雷达信号转换为计算机可识别的特殊格式数据。
[0067]
无人机平台搭载的数据存储模块包括：
[0068]
高性能固态硬盘，用以存储经过初步信号处理的无人机平台地表信息采集的原始数据。
[0069]
便携电源，用以微型中央处理计算单元和高性能固态硬盘供电，防止由于供电不稳造成的数据损坏。
[0070]
微型中央处理计算单元，用以计算处理雷达信号原始数据处理，数据压缩存储处理，及原始数据质量控制信息到地面信息处理平台的数据传输处理。
[0071]
地面信息处理平台包括：
[0072]
高性能计算工作站，用以进行原始数据的数据转移与解压缩，后续时序insar地表形变数据处理步骤，及无人机数据采集过程中的质量控制信息监测。
[0073]
显示单元，用以无人机信息采集过程中的质量控制信息展示，及后续时序insar地表形变数据处理成果展示。
[0074]
上述工作系统的有益效果为：利用无人机平台的高机动性与工作模式高适应性，自行根据需要切换基于不同雷达信号源，对各类地形环境进行地表信息采集；通过信号切换，产生多波段雷达信号源，可以与目前主流insar数据采集平台的采用的雷达信号波段兼容，提高数据空洞填补工作的准确性。
[0075]
基于无人机的时序insar地表形变数据空洞处理方法，包括如下步骤。
[0076]
步骤1、获取前置时序insar地表形变数据集：设置insar地表形变数据孔洞处理的时序变更点，将位于时序变更点之前的insar地表形变数据称为前置时序insar地表形变数据，将位于时序变更点之后的insar地表形变数据称为后置时序insar地表形变数据；按照确定的insar数据范围、重访周期和时序变更点，在每个重访周期均采集一次前置时序insar地表形变数据，从而得到包含空洞位置信息的前置时序insar地表形变数据集。
[0077]
上述前置时序insar地表形变数据集获取的方法，优选包括如下步骤。
[0078]
步骤11、确定目标区域范围：设定需要收集的前置时序insar地表形变数据范围与地理位置。
[0079]
步骤12、确定重访周期：根据待填补地表所处地理位置，确定insar数据系统类型；根据insar数据系统所对应的卫星型号，确定重访周期。
[0080]
步骤13、时序变更点定位：根据连续测量的前置时序insar地表形变数据，定位时序变更点。其中，时序变更点是用于区分前置时序insar与后置时序insar的时间节点，其定位意义为不同的卫星数据源具有不同的重访周期，在对同一区域的进行不同时段的insar数据采集时，使用时序变更点固定每一时段的前、后置时序切换。
[0081]
步骤14、生产前置时序insar地表形变数据集：在每个重访周期，均使用对应的insar数据系统，采集一次前置时序insar地表形变数据，从而得到包含空洞位置信息的前置时序insar地表形变数据集。
[0082]
本步骤14中，前置时序insar地表形变数据集中包括至少15幅时序间距均匀的前置时序insar地表形变数据；步骤12中，insar数据系统类型的确定方法为：
[0083]
a、当待填补地表所处地理位置中人工构筑物占地面积超过设定值或地势起伏度平缓时，insar数据系统使用ps-insar系统。
[0084]
b、当待填补地表所处地理位置中人工构筑物占地面积不超过设定值或地势起伏度陡峭时，insar数据系统使用sbas-insar。
[0085]
步骤2、获取无人机平台后置时序空洞填补数据集：无人机根据步骤1中的空洞位置信息设置航迹，并从步骤1的时序变更点开始，参照步骤1的重访周期和时序变更点，在每个重访周期的时序变更点开始时，均根据设置的航迹采集一次台后置时序空洞填补数据，从而得到无人机平台后置时序空洞填补数据集。
[0086]
步骤21、通过gnss和无人机平台的航迹规划模块与步骤14生成的存在空洞的地表形变前置时序数据集中的空洞位置信息，确定无人机数据采集轨迹，制定目标区域无人机数据采集航迹规划方案。
[0087]
步骤22、通过无人机平台的时序解算模块，拟合步骤1所选用的insar数据采集平台的重访周期和定位的时序变更点，制定无人机数据采集时序方案。
[0088]
步骤23、通过无人机平台的数据采集模块，根据步骤21、步骤22步骤制定的航迹规划方案和时序方案，对目标区域进行后置时序数据采集；数据存放在无人机平台的数据存储模块中。
[0089]
步骤24、将步骤23所获数据由无人机平台的数据存储模块传输至信息处理平台，并经由数据处理步骤后，生成空洞填补数据集。
[0090]
上述步骤24生成空洞填补数据集的方法，优选包括以下分步骤。
[0091]
步骤24-1、将无人机平台存储的数据传输至信息处理平台，包含所采集的后置时序空洞填补数据集、无人机位置信息、数据采集时序方案。
[0092]
步骤24-2、对数据格式进行标准化处理，并统一地理坐标系。
[0093]
步骤24-3、根据时序方案，对数据进行后置时序梳理，以精确对齐in步骤ar数据采集平台。
[0094]
步骤24-4、对冗余数据进行精简化处理，以剔除数据瑕疵和噪声。
[0095]
步骤24-5、生成以处理后的时序为轴的地面高程纹理原始模型。
[0096]
步骤24-6、根据无人机位置信息，对地面高程纹理原始模型进行修正。
[0097]
步骤24-7、根据时序方案，生成后置时序的地表形变数据集，作为空洞填补数据集。
[0098]
上述进一步方案的有益效果为：时序解算是以insar数据采集平台的重访周期为源，让无人机平台的数据采集工作与insar数据采集平台工作产生时空交集，使数据空洞填补处理中的时空序列对齐，可以让最终数据成品更加精确。
[0099]
步骤3、获取后置时序insar地表形变数据集：根据步骤2中无人机平台后置时序空洞填补数据集的数据范围，从步骤1的时序变更点开始，参照步骤1的重访周期和时序变更点，在每个重访周期的时序变更点开始时，均采集一次后置时序insar地表形变数据，从而得到包含空洞位置信息的后置时序insar地表形变数据集。
[0100]
上述后置时序insar地表形变数据集的获取方法，具体参照步骤1，这里将不再赘述。
[0101]
步骤4、获取前置时序空洞填补预测集：对步骤1获取的前置时序insar地表形变数
set。
[0123]
步骤54-2中，每个空洞边缘平滑处理的方法，包括如下步骤。
[0124]
步骤54-2a、计算点位权值：将位于每个空洞边界领域的每个像元分别计算点位权值；其中，点位权值m的计算公式为。
[0125][0126]
式中，w1和w2均为权重系数，设定值；a为当前像元对应的地表形变数据；为当前空洞边界领域中所有像元对应的地表形变数据的均值；b为当前像元至对应空洞边界的距离；表示对应空洞中心点至当前空洞边界的距离均值。
[0127]
步骤54-2b、平滑：将点位权值高于设定权值阈值的像元，采用高斯核函数进行当前像元的平滑计算，并将计算结果赋值给当前像元，实现当前像元的平滑；若当前像元平滑处理效果不能较好表示边界领域的过渡效果，尝试不同步长，重复当前步骤，直至对应像元平滑全部平滑完成。
[0128]
上述进一步方案的有益效果为：通过后置时序insar地表形变空洞边缘监测与边界数据平滑，对空洞边缘的特殊地形，如沟堑、山谷、山脊等存在地形突变的地表情况反映更加贴近真实地貌。根据数据计算估计，使用本处理方法的拟合地表曲率值大于20％的空洞边界邻域数据处理成果对比实测空洞边界邻域地表形变数据的拟合度优于未经本方法处理的数据结果。
[0129]
步骤6、拟合优度检验：将步骤5获取的后置时序空洞填补数据集，采用步骤4获取的前置时序空洞填补预测集进行拟合优度检验，当拟合优度检验指标满足设定阈值时，判断空洞填补完成；否则，重复步骤5至步骤6，直至拟合优度检验指标满足设定阈值要求。
[0130]
上述拟合优度检验指标为可决系数r2，设步骤4获取的前置时序空洞填补预测集为validation set，则r2的计算公式为：
[0131][0132]
式中，rss为validation set与fixed set的残差平方和；tss为validation set与fixed set的总离差。
[0133]
式中，r2的值域为[0,1]，当越接近于1，证明实际观测点离样本线越近，拟合优度越高。
[0134]
如图3所示，通过本发明处理方法的步骤，记录时序变更点，收集、处理获得前置时序与后置时序的时序insar原始数据。时序insar处理模式为sbas-insar。建立时间基线与空间基线连接。对干涉相位进行差分处理，并经过处理运算后获得干涉系数和干涉图。通过相位解缠，获取干涉效果较好的像元位置上被隐藏的主相位。通过对解缠结果解译确定是否继续对时间基线与空间基线进行修改，以获得更好的相位干涉效果和相位解缠结果。处理后获得前置时序与后置时序的时序sbas-insar地表形变数据集，展示成果选取地表平均位移速率结果数据栅格。其中，前置时序sbas-insar地表形变数据集在经s4步骤处理后，在本实施例中的有益效果为仅根据自然条件下的地表形变发生趋势作为数据集参考标准，在毫米级尺度上的山区地表形变评价中，与未经过s4步骤处理而使用常规插值的数据，对后续步骤的拟合优度判断具备更好的参考性。
[0135]
如图4所示，时序sbas-insar地表形变速率结果数据栅格，可见明显数据空洞。使用本发明的处理方法，对前置时序insar地表形变数据集进行时序预测性模型训练，以生成空洞填补预测集，用于空洞填补数据集的拟合优度验证；使用本发明的处理方法，导入无人机平台采集所得的用于地表形变空洞填补的数据集，对该数据集进行投影坐标系转换、重采样、边界提取等操作，获得地表形变空洞填补数据集，并经由边界灰度匹配检验，验证空洞填补有效性，并根据检验结果对边界邻域使用边缘平滑计算。计算结果表明，本研究区的拟合地表曲率值大于20％的部分，使用本方法提出的边界灰度匹配检验与边缘平滑处理结果数据对比未使用本方法的数据，对地表形变的反映更加接近真实山区地貌情况，对空洞边缘的特殊地形，如沟堑、山谷、山脊等存在地形突变的地表情况的反映细节更加明显。
[0136]
如图5所示，实施本发明处理方法的时序insar数据空洞填补，经过拟合优度检验后，完成时序insar地表形变数据空洞处理。空洞填补数据中地表形变年平均速率单位为mm/年，数据范围为-对比未经处理的数据。
[0137]
本发明的实施过程为：采集目标区域的前置时序insar地表形变数据栅格；无人机工作平台后置时序insar地表形变数据空洞填补数据采集与处理；采集目标区域的后置时序insar地表形变数据栅格；前置时序insar地表形变数据空洞填补预测训练；后置时序insar地表形变数据空洞填补处理；根据预测训练模型，验证空洞填补数据的拟合优度，完成时序insar地表形变数据空洞处理。
[0138]
以上通过实施例对本发明作了具体描述。需要理解的是，术语“patch set”、“validation set”、“raw set”、“fixed set”等指代的数据集名称，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的数据集合必须命名为特定的称呼，因此不能理解为本发明的限制。因此，限定由“patch set”、“validation set”、“raw set”、“fixed set”的数据集可以明示或隐含地包括一个或多个本发明提供的数据处理集合。
[0139]
以上详细描述了本发明的优选实施方式，但是，本发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种等同变换，这些等同变换均属于本发明的保护范围。