一种适用于大梯度地表沉降监测的时序差分雷达干涉方法技术

本申请公开了一种适用于大梯度地表沉降监测的时序差分雷达干涉方法，该方法通过短时间基线差分干涉图筛选、离散点相位解缠、基于短时间基线差分干涉图的形变分量建模和解算、形变分量可靠性检验、差分干涉图相位梯度修正、对上述过程进行迭代以确保形变分量解算正确，以及基于修正后差分干涉图的形变时间序列建模和解算等过程，并形成整体的技术方案，解决原有时序差分雷达干涉技术中由于形变和相位梯度较大导致形变相位模糊度解算困难或失败的问题，最终达到正确提取大梯度地表形变速率和形变时间序列的目的，并起到降低大梯度形变建模和解算所需合成孔径雷达影像数量的效果，节约时序差分雷达干涉的应用经济成本。

A method of time difference radar interferometry for large scale land surface subsidence monitoring

The invention discloses a suitable for the large gradient of surface subsidence monitoring of temporal differential radar interferometry method, this method through the short baseline differential interferogram filtering, discrete phase unwrapping, short baseline differential interferometry diagram, component modeling and deformation component reliability test, differential interferogram phase gradient the correction and iterative process to ensure that the solution is correct based on the deformation of components, and based on the modified differential interferometric deformation time series modeling and calculation process diagram, and the formation of the overall technical scheme, the solution of the original time series differential Radar Interferometry Technology due to the deformation and phase gradient due to the large deformation phase ambiguity difficulties or failure problems, finally achieve the correct extraction of large gradient surface strain rate and deformation time series, and to reduce the high gradient deformation modeling The computational cost of the SAR image is calculated, and the economical cost of the application of the time series differential radar interference is saved.

【技术实现步骤摘要】
一种适用于大梯度地表沉降监测的时序差分雷达干涉方法
本专利技术属于地表沉降监测
，具体地说，涉及一种适用于大梯度地表沉降监测的时序差分雷达干涉方法。
技术介绍
地表沉降(垂直向地表形变，为便于表述，后续采用“形变”代表“沉降”)是发生范围最广的地质灾害之一，具有持续时间长的特点，且多发于城市及其近郊等经济发达和人口聚集区，对经济发展、城市建设和人民生活均会构成持久危害。我国是世界上地表沉降灾害最为严重的国家之一，累积沉降大于200毫米的面积超过15万平方公里，主要集中在华北平原、长江三角洲和汾渭盆地等经济发达地区，且出现了严重的沉降漏斗，造成了严重的经济损失。对地表沉降(尤其是沉降漏斗)开展大范围精确监测，对沉降防控及避免相应的危害具有十分重要的现实意义。目前，时序差分合成孔径雷达(SyntheticApertureRadar,SAR)干涉(DifferentialSARInterferometry,DInSAR)技术已在地表形变监测中得到广泛应用，且是微波遥感、卫星大地测量以及地球物理学领域研究和应用的热点之一。时序DInSAR(TimeSeriesDInSAR,TS-DInSAR)具有覆盖范围广、空间分辨率高、效率高、精度高且不易受云雾和阴雨天气影响的技术优势，非常适用于开展大范围地表形变监测。TS-DInSAR对缓慢累积性地表形变具有较好的监测效果和精度，但当形变较快和形变梯度较大时易造成解算精度较低或解算失败(如形变欠估计和形变漏斗区的“空洞”现象，即结果缺失)。
技术实现思路
有鉴于此，本申请针对TS-DInSAR难以满足快速或大梯度地表形变监测需求的问题，提供了一种适用于大梯度地表沉降监测的时序差分雷达干涉方法。为了解决上述技术问题，本申请公开了一种适用于大梯度地表沉降监测的时序差分雷达干涉方法，具体按照以下步骤实施：步骤1，对所有筛选后的SAR影像进行任意干涉组合并计算时间基线(形成一个干涉对的两幅SAR影像的获取时间差)和空间基线(两次成像时刻SAR传感器之间的空间距离，一般取该距离垂直于SAR视线方向的分量)；步骤2，设置时空基线阈值进行干涉(干涉组合)对初始筛选，在保证干涉对数量的前提下限制时间失相干和空间失相干，按照干涉组合进行差分干涉处理得到差分干涉图集；步骤3，永久散射体(PersistentScatterer,PS)和相干目标(CoherentTarget,CT)探测及点集合并，得到相干散射体(CoherentScatterer,CS)点集，提取CS点集上的差分干涉相位，并以CS为节点构建不规则三角网(TriangularIrregularNetwork,TIN)，本专利技术技术方案中采用Delaunay三角网；步骤4，基于Delaunay三角网和最小费用流(MinimumCostFlow,MCF)方法的CS相位解缠；步骤5，线性形变速率和高程误差建模及解算；步骤6，基于数据模拟或地面测量数据的线性形变解算结果检验，若结果不可靠则执行步骤7，否则转向步骤8；步骤7，重新筛选参与计算的差分干涉图，并重新执行步骤5和步骤6，当达到计算终止条件时，执行步骤8；步骤8，从所有原始差分干涉图中减去线性形变相位分量，对残差相位(从差分干涉图中减去线性形变相位后剩余的相位)进行重新解缠，然后再将线性形变相位分量加回重新解缠后的相位中，重新执行步骤5，得到更新后的线性形变速率和高程误差(这一步主要是让更多干涉对参与计算，提高线性形变速率和高程误差的解算精度)；步骤9，记录步骤8中所得到的线性形变速率，并从所有的原始差分干涉图中减去更新后的线性形变和高程误差相位分量，重新执行基于离散点的相位解缠，然后将步骤8中的线性形变分量重新加回新的解缠相位中，执行形变时间序列建模和解算过程，最终得到形变时间序列。与现有技术相比，本申请可以获得包括以下技术效果：本专利技术技术方案依据TS-DInSAR时序差分干涉图中形变相位大小以及形变相位梯度大小与时间间隔(差分干涉图的时间基线，即形成差分干涉图的两幅SAR影像获取时间差)成正相关关系这一特点，提出短时间基线差分干涉图筛选、离散点相位解缠、基于短时间基线差分干涉图的形变分量建模和解算、形变分量可靠性检验、差分干涉图相位梯度修正、对上述过程进行迭代以确保形变分量解算正确，以及基于修正后差分干涉图的形变时间序列建模和解算这一整体的技术方案，解决原有TS-DInSAR技术中由于形变和相位梯度较大导致形变相位模糊度解算困难或失败的问题，最终达到正确提取大梯度地表形变速率和形变时间序列的目的。采用本专利技术技术方案，只要差分干涉图序列中存在至少6个短时间基线干涉对(由4幅SAR影像构成)即可实现线性形变分量的解算，然后开展长时间基线差分干涉图的形变相位梯度修正，促使更多的差分干涉图得到正确解缠，并参与形变分量的重新估算以及形变时间序列的建模和解算。因此，无需很高的SAR影像使用量便可实现大梯度形变的有效提取。当然，实施本申请的任一产品并不一定需要同时达到以上所述的所有技术效果。附图说明此处所说明的附图用来提供对本申请的进一步理解，构成本申请的一部分，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。在附图中：图1是一种适用于大梯度地表沉降监测的时序差分雷达干涉方法的实施流程图；图2是35幅SAR影像任意组合所形成的干涉对的时间基线和空间基线分布图；图3是空间基线阈值为30米时干涉对的时空基线分布图；图4是第一次解算所得CS点上线性形变速率结果图；图5是使用第一次解算所得的线性形变模拟的差分干涉图与原始差分干涉图的对比；其中，a为原始差分干涉图，b为使用第一次解算所得的线性形变模拟的差分干涉图；图6是经3次迭代后解算所得CS点上的线性形变速率结果图；图7是使用第三次解算所得的线性形变模拟的差分干涉图与原始差分干涉图的对比；其中，a为原始差分干涉图，b为使用第三次解算所得的线性形变模拟的差分干涉图；图8是2009年6月23日至2010年7月2日期间的累积形变量；图9是图8中所标示三个特征点(A、B和C)的形变时间序列。具体实施方式以下将配合附图及实施例来详细说明本申请的实施方式，藉此对本申请如何应用技术手段来解决技术问题并达成技术功效的实现过程能充分理解并据以实施。一种适用于大梯度地表沉降监测的时序差分雷达干涉方法，如附图1所示，具体按照以下步骤实施：步骤1，对所有筛选后的SAR影像进行任意干涉组合并计算时间基线和空间基线；在筛选出合适的SAR影像(排除受雨雪等天气以及积雪影响的SAR影像)后，进行任意干涉组合配对，假设有N+1幅SAR影像，通过任意干涉组合可形成N(N+1)/2个干涉对。每个干涉对由主、从两幅SAR影像构成。然后，根据每个干涉对主、从SAR影像的获取时间计算该干涉对的时间基线(即SAR影像获取的时间差)，根据主、从SAR影像的参数文件所记录的SAR传感器运行位置及相应的时间参数计算该干涉对的空间基线(主、从影像成像时SAR传感器的空间距离，实际中一般取该空间距离在垂直于SAR传感器视线方向上的分量为空间基线，也即垂直基线)。本实施例中采用覆盖天津市精武镇的35幅SAR影像为数据进行展示。附图2所示为35幅SAR影像进行任意组合所本文档来自技高网...

【技术保护点】
一种适用于大梯度地表沉降监测的时序差分雷达干涉方法，其特征在于，具体按照以下步骤实施：步骤1，对所有筛选后的SAR影像进行任意干涉组合并计算时间基线和空间基线；步骤2，设置时空基线阈值进行干涉对初始筛选，在保证干涉对数量的前提下限制时间失相干和空间失相干，按照干涉组合进行差分干涉处理得到差分干涉图集；步骤3，永久散射体PS和相干目标CT探测及点集合并，得到相干散射体CS点集，提取CS点集上的差分干涉相位，并以CS为节点构建不规则Delaunay三角网TIN；步骤4，基于Delaunay三角网和MCF方法的CS相位解缠；步骤5，线性形变速率和高程误差建模及解算；步骤6，基于数据模拟或地面测量数据的线性形变解算结果检验，若结果不可靠则执行步骤7，否则转向步骤8；步骤7，重新筛选参与计算的差分干涉图，并重新执行步骤5和步骤6，当达到计算终止条件时，执行步骤8；步骤8，从所有原始差分干涉图中减去线性形变相位分量，对残差相位进行重新解缠，然后再将线性形变相位分量加回重新解缠后的相位中，重新执行步骤5，得到更新后的线性形变速率和高程误差；步骤9，记录步骤8中所得到的线性形变速率，并从所有的原始差分干涉图中减去更新后的线性形变和高程误差相位分量，重新执行基于离散点的相位解缠，然后将步骤8中的线性形变分量重新加回新的解缠相位中，执行形变时间序列建模和解算过程，最终得到形变时间序列。...

【技术特征摘要】
1.一种适用于大梯度地表沉降监测的时序差分雷达干涉方法，其特征在于，具体按照以下步骤实施：步骤1，对所有筛选后的SAR影像进行任意干涉组合并计算时间基线和空间基线；步骤2，设置时空基线阈值进行干涉对初始筛选，在保证干涉对数量的前提下限制时间失相干和空间失相干，按照干涉组合进行差分干涉处理得到差分干涉图集；步骤3，永久散射体PS和相干目标CT探测及点集合并，得到相干散射体CS点集，提取CS点集上的差分干涉相位，并以CS为节点构建不规则Delaunay三角网TIN；步骤4，基于Delaunay三角网和MCF方法的CS相位解缠；步骤5，线性形变速率和高程误差建模及解算；步骤6，基于数据模拟或地面测量数据的线性形变解算结果检验，若结果不可靠则执行步骤7，否则转向步骤8；步骤7，重新筛选参与计算的差分干涉图，并重新执行步骤5和步骤6，当达到计算终止条件时，执行步骤8；步骤8，从所有原始差分干涉图中减去线性形变相位分量，对残差相位进行重新解缠，然后再将线性形变相位分量加回重新解缠后的相位中，重新执行步骤5，得到更新后的线性形变速率和高程误差；步骤9，记录步骤8中所得到的线性形变速率，并从所有的原始差分干涉图中减去更新后的线性形变和高程误差相位分量，重新执行基于离散点的相位解缠，然后将步骤8中的线性形变分量重新加回新的解缠相位中，执行形变时间序列建模和解算过程，最终得到形变时间序列。2.如权利要求1所述的一种适用于大梯度地表沉降监测的时序差分雷达干涉方法，其特征在于，步骤1具体实施方式为：排除受雨雪天气以及积雪影响的SAR影像后进行任意干涉组合配对，假设有N+1幅SAR影像，通过任意干涉组合形成N(N+1)/2个干涉对，每个干涉对由主、从两幅SAR影像构成，然后，根据每个干涉对主、从SAR影像的获取时间计算该干涉对的时间基线即SAR影像获取的时间差；根据主、从SAR影像的参数文件所记录的SAR传感器运行位置及相应的时间参数计算该干涉对的空间基线，主、从影像成像时SAR传感器的空间距离，取该空间距离在垂直于SAR传感器视线方向上的分量为空间基线，也即垂直基线。3.如权利要求1所述的一种适用于大梯度地表沉降监测的时序差分雷达干涉方法，其特征在于，步骤2中在进行时序差分干涉处理时，为降低时间失相干和空间失相干的影响，采用时间基线阈值和空间基线阈值方法排除时间基线和空间基线大于某给定的阈值的干涉对，当SAR影像获取时间整体跨度小于2年时不考虑对时间基线进行限制。4.如权利要求1所述的一种适用于大梯度地表沉降监测的时序差分雷达干涉方法，其特征在于，步骤3具体实施方式为：PS探测采用振幅阈值和振幅离差指数ADI阈值双阈值方法，对于SAR影像中某一特定的像元，其在N+1幅SAR影像中的振幅值直接从SAR振幅影像中提取，则其时序振幅均值和ADI值为：其中，i为SAR影像索引号(序号)；ai和分别为该像元在第i幅SAR影像中的振幅值及在所有影像中的振幅值的均值；Da和σa分别为该像元的ADI值及时序振幅标准差。当该像元的振幅信息满足如下条件时，认为其为PS：其中，和σA分别为N+1幅SAR影像所有像元时序振幅值的均值和标准差；c和l分别为像元的列号和行号，C和L分别为影像列数和行数；Ai和acl分别为第i幅SAR影像所有像元振幅值的均值和行、列号分别为c和l的某像元的振幅值，公式(3)中0.25即为ADI阈值，即为振幅阈值；CT探测采用相干系数阈值方法，假设由N+1幅SAR影像形成了L个干涉对，并通过差分干涉数据处理得到L幅差分干涉图，通过下式计算每个像元在所有干涉图中相应的相干系数：其中，为某像元在第l幅干涉图中的相干系数值；IMl和ISl分别为第l个干涉对的主影像和从影像；Z为相干系数估计窗口内像元数目，z为窗口内像元索引；l为干涉图索引；当某个像元满足以下条件时被识别为CT：其中，min(·)表示取变量的最小值；γcrit为判别某个像元是否为CT的相干系数阈值，取0.25至0.3；当探测出所有的PS和CT后，将二者进行合并，并去除重复点，得到CS点集，最后从所有差分干涉图中提取CS点上的相位值。5.如权利要求1所述的一种适用于大梯度地表沉降监测的时序差分雷达干涉方法，其特征在于，步骤4中相位解缠采用基于离散点的相位解缠方法对CS上的差分干涉相位进行解缠处理；对于某个CS点P，其在差分干涉图中的相位和该点上的真实相位的关系为：其中，为P点上的真实相位，即相位解缠要得到的相位值；φP为P点上的缠绕相位值即P在差分干涉图中对应的相位值，位于[-π,π)区间内，只记录了不足整周(2π)的小数部分，存在整周模糊度；nP为整周模糊数；首先，根据Delaunay剖分法则，以所有的CS为顶点构造Delaunay三角网，然后以网络边端点对应的两个CS之间的缠绕相位差为观测量，估算两点间的绝对相位差，通过最小费用流MCF方法对与相位不连续性相关的网络费用流进行估算，并寻找最小费用流对应的积分路径，进行相位积分，完成相位解缠，此过程也即求解公式(6)中nP的过程。6.如权利要求1所述的一种适用于大梯度地表沉降监测的时序差分雷达干涉方法，其特征在于，步骤5中在对CS上的差分干涉相位进行相位解缠后，根据步骤4中构建的Delaunay三角网重新计算相邻CS点间的解缠相位差，此处，相邻的CS点是指Delaunay三角网中三角形边的端点，以任意一边的两端点P和Q为例，二者在第l个差分干涉图中对应于的解缠相位为：其中，和分别为第l幅差分干涉图的时间基线和垂直基线；λ为雷达波波长；θP和θQ分别为P和Q点处的雷达波入射角；RP和RQ分别为雷达天线到P和Q之间的斜距；和分别为P和Q点上解缠后的差分干涉相位；vP和vQ分别为P和Q点的线性形变速率；δhP和δhQ分别为P和Q点的高程误差；和分别为P和Q点在第l幅差分干涉图中的非线性形变相位；和分别为P和Q点在第l幅差分干涉图中的轨道误差相位；和分别为P和Q点在第l幅差分干涉图中的大气延迟相位；和分别为P和Q点在第l幅差分干涉图中的噪声相位；线性形变和高程误差建模及解算的目的是对vP和vQ及δhP和δhQ进行估算，以P和Q上的解缠相位为观测量进行网络邻域差分建模，对于P和Q，二者之间的网络邻域差分相位增量为：其中，和分别为P和Q点处斜距的平均值和入射角的平均值；为P和Q之间的邻域差分相位增量；ΔvPQ为P和Q之间的线性形变速率增量；ΔδhPQ为P和Q之间的高程误差增量；为第l幅差分干涉图中P和Q之间的相位残差增量，即P和Q点上非线性形变、轨道误差、大气延迟和噪声相位之间的差值之和；以L幅差分干涉图为例，对于P和Q所对应的网络边而言，列出L个与式(9)相同的方程，组成相应的线性方程组，将其表达为矩阵的形式有：ΔΨ＝AX+W(10)其中，A＝[κ,η](1...

【专利技术属性】
技术研发人员：于冰，马德英，肖东升，杨莹辉，刘福臻，贾宏亮，苏勇，王继燕，
申请(专利权)人：西南石油大学，
类型：发明
国别省市：四川,51