一种基于干涉图闭合环的InSAR解缠误差探测方法与流程

本发明涉及一种insar解缠误差的探测方法，具体是基于重访周期较短的多个复数影像，使用小基线技术组成的解缠对，适用于任何波段、平台、频率获取的大量干涉对解缠误差探测的情况。

背景技术：

为了更准确和高效的获取地球表面的变化信息及地球内部的变化规律，已有一些较为成熟的大地测量观测手段，如：水准、经纬仪、全球观测系统等广泛应用于地球环境监测。然而这些技术虽灵活便利、具有较高的时间分辨率，但是对于大区域长时间的形变监测仍存在许多不足：1)网点布设需要较大的劳动强度及昂贵的成本花费；2)点状观测特征难以满足重点区域高分辨率的需求；3)在危险性较高、人无法到的区域、现代高楼林立的城市特征，这些技术无法获取数据或者受误差影响较大，无法满足形变观测要求。

合成孔径雷达干涉测量(inasr)自上世纪九十年代成功应用于地震形变监测以来，因其高分辨率、全天候获取、无需地面辅助的优点得到广泛发展，而且已经成功应用在火山喷发、冰川漂移、滑坡、泥石流、地震、矿山沉陷等地质灾害的地表形变监测^[1-3]。insar采用的主动式微波成像，相位值为地表与传感器之间距离的变化量，强度值反映了地表的散射特征。同时，微波相位易受时空失相干和大气误差的影响，针对上述两项误差，目前已有永久散射体技术和小基线集技术和两者的结合来最大限度的消除上述误差。然而，除了轨道误差、大气效应、失相干等误差影响，一个周期的解缠误差对形变的贡献量为λ/2，对c波段而言即是2.8cm，因此，解缠误差，作为insar技术中一个关键性问题，已不容忽略，否则会造成污染到最终参数的求解。

如何正确探测干涉图中的解缠误差并最大化的保留干涉图的相位信息，目前国际上主要的方法主要有以下几种：1)将解缠误差和形变量及dem误差参数同时解算，将含有解缠误差的相干点作为粗差滤除。永久散射体技术中，选择裸露岩石等长时间内保持相干性的相干点，组成三角网，弧段间因距离较短，大气误差可忽略，相差值如果大于设定的阈值，则认为含有解缠，并将含有解缠误差的弧段作为粗差去除，同时解算弧段间的形变及dem误差，通过网平差来最终获取全局相干点的形变量和dem误差，具体详见参考文献^[4]。2)将干涉图预处理，探测出解缠误差，小基线技术中通过干涉图闭合环中的粗差值探测解缠误差，将闭合环残差较大的像素在组成闭合环的三个干涉图中掩模掉，详细步骤在干涉图^[5]。将解缠误差和形变参数即dem误差一起解算中，没有考虑大气误差的影响，容易造成解缠误差的不准确判断；仅用单个闭合环检查干涉图的预处理方法中，并不能确定误差所在哪个干涉图上，而且闭合环会有呈现有规律的轨道误差，收部分低相干区域的影响，解缠误差并不总是呈现明显的块状分布，逐个手动检查来定位解缠误差的位置和范围并不现实。因此，上述方法无法有效确定干涉图上的解缠误差。

通过以上分析可总结，在预处理阶段有效定位解缠误差所在像素并最大化保留干涉图的干净的相位信息显得尤为重要，现有的方法并不能满足此要求。实际上，干涉图中会同时存在形变量、轨道误差、dem误差、大气噪声、解缠误差、随机噪声，如何有效的分离出解缠误差并没有得到有效的解决。

[1]massonnet,d.,etal.(1993)."thedisplacementfieldofthelandersearthquakemappedbyradarinterferometry."natureinternationalweeklyjournalofscience364(6433):138‐142.

[2]lu,z.,etal.(1997)."deformationofnewtridentvolcanomeasuredbyers‐1sarinterferometry,katmainationalpark,alaska."geophysicalresearchletters24(6):695‐698.

[3]massonnet,d.andk.l.feigl(1998)."radarinterferometryanditsapplicationtochangesintheearth'ssurface."reviewsofgeophysics36(4):441–500.

[4]ferretti,a.,etal.(2011)."anewalgorithmforprocessinginterferometricdata‐stacks:squeesar."ieeetransactionsongeoscienceandremotesensing49(9):3460‐3470.

[5]berardino,p.,etal.(2002)."anewalgorithmforsurfacedeformationmonitoringbasedonsmallbaselinedifferentialsarinterferograms."ieeetransactionsongeoscienceandremotesensing40(11):2375‐2383.

技术实现要素：

针对上述现有技术中存在的问题，本发明的目的在于，克服出现有的insar解缠误差的探测技术局限性，基于多个干涉对闭合环残差的探测，提供一种不受其他误差的影响，探测出insar干涉对中解缠误差的方法。

为了实现上述任务，本发明采用以下技术方案：

一种基于干涉图闭合环的insar解缠误差的探测方法，包括以下步骤：

步骤一，获取研究区域的多景单视复数影像，进行多视、滤波、干涉去平、去除地形相位、解缠操作，获取解缠的insar差分干涉图；

步骤二，选择要探测解缠误差的目标干涉对，探测能与目标干涉对组成闭合环的其他两个干涉对，得到多组干涉图闭合环，解算每组中三个干涉对的残差；

步骤三，对每组干涉对闭合环残差进行多项式拟合，并用最小二乘法求解多项式参数，分别求解闭合环残差与多项式拟合结果的差值；

步骤四，取各组闭合差结果中绝对值小于0.5弧段的像素，其他像素掩模掉，将保留的各个差值干涉对取并集；

步骤五，将并集结果与目标干涉对像素取交集，目标干涉对掩模掉的像素便是解缠误差存在的像素。

进一步地，步骤二中，解算每组中三个干涉对的残差的公式为：

上式中，为解算后得到的残差，为目标干涉对，为与目标干涉对组成闭合环的其他两个干涉对，a1，a2，…，ak为与干涉对组成闭合环的所有slc编号。

进一步地，步骤三的具体过程包括：

检查闭合环残差结果，对其闭合环残差通过多项式估计去除轨道误差；如不含有轨道误差，则令改正后的闭合环残差

如果含有轨道误差，则对进行二项式拟合，采用的多项式为：

pi＝b0+b1x+b2y+b3xy

上式中，b0，b1，b2和b3为多项式的系数，x，y为多项式的参数；通过最小二乘法求解多项式的系数，并令改正后的闭合环残差：

上式中，和为通过最小二乘法估计的多项式系数。

进一步地，所述的步骤四的过程包括：

针对于步骤三得到的改正后的闭合环残差如果说明满足条件的像素在干涉对中均不存在解缠误差，仅含有随机性误差，则将不满足要求的像素掩模掉；其中，ε为设定的阈值，ε≤0.5·2π；

将保留的各个差值干涉对取并集：

上式中，表示所有含有干涉对的闭合环并集，a1，a2，…，ak为与干涉对组成闭合环的所有slc编号。

本发明与现有技术相比具有以下技术特点：

本发明方法较现有方法相比，更容易实现，充分去除了inasr观测量中的大气误差、残余轨道误差和形变信息，更有利于研究解缠误差位置，有效的分离出解缠误差，并避免对其他参数的求解带来误差传播，适合用于有大量干涉对的情况，进一步推动了insar大地测量技术的应用进程。

附图说明

图1为本发明方法的流程图；

图2为实施例中三个干涉对：20150413_20150718、20150413_20150730、20150718_20150730组成的闭合环及多项式拟合结果，其中：a为三个干涉对的原始闭合环残差，b为a结果的直方图统计图，c为经轨道改正后的结果，d为c结果的直方图估计，e为估计的轨道误差，f为设置阈值去除粗大残差后的结果，单位为弧段/2pi；

图3为改正前后的20150413_20150718干涉图对比，其中a为未去除解缠误差的差分干涉图，b为经处理后的干涉图，c、d分别为a和b红色框选的放大部分，单位为弧度。

具体实施方式

本发明的理论基础阐述如下：

地球表面的各种地质灾害均以地表形变的形式体现出来，而利用insar观测技术探测地表形变时，解缠误差是一项不可忽略的因素。

不存在解缠误差的干涉对组合中，系列公式恒等：

为与干涉对组成闭合环的其他两个干涉对，指在时间点i和时间点j获取的slc交叉相乘形成的干涉对(相位)，即在时间点i和时间点j之间，测得卫星与地面目标距离的变化量；其贡献量由以下各项构成：

上式中，表示地表在视线方向形变量，视线方向为卫星和地面目标连线的方向；表示因两个slc入射角不完全重合，残余的轨道误差；表示用于地形改正的外部dem误差；表示干涉图上解缠误差；表示干涉图上的大气误差；表示时间失相干和传感器热噪声的随机误差的综合影响。

在干涉对组成闭合环残差中：

上式中：

为闭合环残差中的形变相位分量，

为大气相位残差，

为残余轨道误差，

为dem相位分量，

解缠误差分量，

随机噪声，

干涉图中的形变和大气影响仅与获取时间有关并与时间呈线性关系，闭合环残差后恒等于零，残差结果仅含有后四项：

上式中，与均表现出随机性，呈趋势面分布，在干涉图闭合差表现为跳变，因为通过多项式拟合去除趋势面则可得到所在区域：

为将去除轨道误差的结果，n是随机误差，含有和

将含有轨道误差的闭合环残差通过多项式估计去除轨道误差；采用的多项式pi＝b0+b1x+b2y+b3xy，通过最小二乘求解多项式系数，并改正闭合环残差式中为最小二乘估计的多项式系数。

假设闭合环不存在解缠误差的像素，在组成闭合环的三个中也不存在解缠误差。

如果满足条件的像素在干涉对中均不存在解缠误差，仅含有随机性误差，将不满足要求的像素掩模掉，ε为事先设定的阈值，一般取ε小于或者等于0.5·2。

但，存在解缠误差的情况下，三个干涉图的相位值全部被掩模掉，即，如果或某像素存在有解缠误差，的会被掩模，造成有效信息丢失。

通过求解多个闭合环中，例如干涉图组成的闭合环像素并不存在解缠误差，则通过两个干净闭合环残差的并集，的像素的相位又被包括进来，不至于受其他两个干涉对的影响太大。

对于干涉图ij，如果有k个闭合环，则取闭合环之间的并集，保存的像素中均不含有解缠误差。

上式中，为不含有解缠误差的闭合环，即改正后的闭合环残差；为原始的干涉图；为以为模板，保留不含有解缠误差的相位，含有误差的相位已通过模板掩模掉。获取的既去除了解缠误差，又保留了有效消息。

本发明的具体步骤如下：

一种基于干涉图闭合环的insar解缠误差的探测方法，包括以下步骤：

步骤一，获取研究区域的多景单视复数(slc)影像，进行多视、滤波、干涉去平、去除地形相位、解缠操作，获取解缠的insar差分干涉图；

本方案中，研究数据获取自欧空局网站，查询研究区域对应数据并下载单数复数影像，进行方位向6距离向30的多视操作，减少斑点噪声，然后进行滤波、干涉去平处理，下载对应区域的dem，并地理编码到sar坐标系下，减去高程相位，再利用最小费用流方法解缠，获取解缠后的差分干涉图。上述步骤的原理是通过两个单视复数相位交叉相乘，减去模拟的高程相位，获取常规的差分干涉图。

步骤二，选择要探测解缠误差的目标干涉对，探测能与目标干涉对组成闭合环的其他两个干涉对，得到多组干涉图闭合环，解算每组中三个干涉对的残差；具体过程为：

选择要检验的差分干涉图，并搜索与此干涉图组成的闭合环的个数；

解算多个闭合环残差，其公式表示为：

上式中，为解算后得到的残差，为目标干涉对，为与目标干涉对组成闭合环的其他两个干涉对，其中的闭合环，指三个干涉对构成的相位残差；a1，a2，…，ak为与干涉对组成闭合环的所有slc编号。

步骤三，对每组干涉对闭合环残差进行多项式拟合，并用最小二乘法求解多项式参数，分别求解闭合环残差与多项式拟合结果的差值；具体过程包括：

检查闭合环残差结果，对其闭合环残差通过多项式估计去除轨道误差；如不含有明显的轨道误差，则令改正后的闭合环残差所述的轨道误差，是指中如果含有趋势项，即表明有明显的轨道误差。

如果含有轨道误差，则对进行二项式拟合，采用的多项式为：

pi＝b0+b1x+b2y+b3xy

上式中，b0，b1，b2和b3为多项式的系数，x，y为多项式的参数；通过最小二乘法求解多项式的系数，并令改正后的闭合环残差：

上式中，和为通过最小二乘法估计的多项式系数。

检查残差值轨道误差在干涉图上表现为趋势项，同样，可传递到闭合环残差，去除趋势项后更好的探测解缠误差。

步骤四，取各组闭合差结果中绝对值小于0.5弧段的像素，其他像素掩模掉，将保留的各个差值干涉对取并集；具体过程包括：

针对于步骤三得到的改正后的闭合环残差如果说明满足条件的像素在干涉对中均不存在解缠误差，仅含有随机性误差，则将不满足要求的像素掩模掉；其中，ε为设定的阈值，ε≤0.5·2π；

将保留的各个差值干涉对取并集：

上式中，表示所有含有干涉对的闭合环并集，并且已经将含有解缠误差的像素掩模掉，仅含有热噪声失相关等随机性误差；a1，a2，…，ak为与干涉对组成闭合环的所有slc编号。不含有解缠误差的干涉图其闭合环残差接近零值，如果某一闭合环残差值较小，则组成闭合环残差的三个干涉对均不含有解缠误差。中所有像素不含解缠误差。

步骤五，将并集结果与目标干涉对像素取交集，目标干涉对掩模掉的像素便是解缠误差存在的像素，具体过程包括：

其中，表示时间点i和时间点j组成的干涉对，分别是含有解缠误差和改正后的干涉对。通过多个含有干涉图的闭合环，最大可能的探测每一个闭合环中干涉图表现的非解缠误差区域，取每个闭合环结果的并集，最大限度的保留相位的相位信息，更为真实的反映解缠误差像素。

实施例：

步骤一，获取研究区域的多景单视复数影像，进行多视、滤波、干涉、去平、去除地形相位、解缠操作，获取解缠的insar差分干涉图；通过实测墨西哥城区的sentinel‐1a数据(欧空局获取)，34景slc影像组成560个干涉图。基于gamma平台，干涉对距离上31和方位向6的多视处理获得分辨率150x150的像素，经自适应滤波、由外部dem模拟地形相位、去除地形相位的贡献量、最小费用流解缠，获得干涉对。

步骤二，选择要探测解缠误差的目标干涉对，探测能与目标干涉对组成闭合环的其他两个干涉对，得到多组干涉图闭合环，解算每组中三个干涉对的残差；选择20150413_20150718干涉对(表1所示)作为检验的目标干涉对，目标干涉与其他干涉对组成闭合环的其他基线情况如表1所示，共组成31个闭合环。此处仅随机选择一个闭合环的结果呈现，在表1第5行加重字体：20150413_20150718、20150413_20150730、20150718_20150730三个干涉对。为了更方便的可视化解缠误差，图2将残余的相位除以2π，不在以相位弧度为单位，闭合环原始结果如a所示，b中统计了闭合环结果的直方图，从直方图中可看到原始干涉图中存在明显的解缠误差，及多项式拟合结果。从而由最小二乘模拟出多项式系数，如图e所示，将原始闭合差结果中去除模拟的多项式结果如c所示，图d为c结果的直方图估计；去除一些粗大的残差，结果如f所示。

表1.与20150413_20150718组成闭合环的情况

步骤三，对每组干涉对闭合环残差进行多项式拟合，并用最小二乘法求解多项式参数，分别求解闭合环残差与多项式拟合结果的差值；

步骤四，取各组闭合差结果中绝对值小于0.5弧段的像素，其他像素掩模掉，将保留的各个差值干涉对取并集；

步骤五，将并集结果与目标干涉对像素取交集，目标干涉对掩模掉的像素便是解缠误差存在的像素。取另外31个闭合环非被解缠误差污染像素的并集，再与所求干涉对解算交集，结果如图3所示，对a和b可发现并没有很多的信息被掩模掉，干涉图左部分因失相干较为严重，表现了较大误差，但此部分不是我们的研究范围。将试验区域重点方法对比图c和d可发现，解缠误差被掩模掉，并保留了不含解缠误差的部分，通过对比探测出的解缠误差的位置与相邻像素的对比，基本上均存在2π的整数倍跳变，进一步验证了此方法的正确性。