一种多维层析SAR古迹遗址形变监测方法

摘要

本发明公开了一种多维层析SAR古迹遗址形变监测方法，包括：收集待监测古迹遗址的SAR影像；提取SAR影像中的第一PS候选点，对第一PS候选点构建Delaunay三角网，进行弧段筛选，得到初选后的三角网及其弧段所连接的第二PS候选点；对所述进行预处理后的SAR影像进行干涉对组合，得到N幅差分干涉图；将所述初选后的三角网的每个弧段所连接的两个第二PS候选点作为邻近点，对所述差分干涉图建立邻近点之间相位差分模型；建立所述第二PS候选点沿多维空间坐标轴的离散化表达模型，采用非参数谱估计法中的波束形成法重构所述层析信号，进行弧段筛选，提取第一高置信度PS点和第一高置信度MPS点。通过本发明，为古迹遗址的监测和管理提供了科学依据。

说明书

技术领域

本发明涉及遥感影像的数字摄影测量技术领域，更具体地，涉及一种多维层析SAR古迹遗址形变监测方法。

背景技术

古迹遗址本体及其赋存环境的形变监测一直以来都是各国科学家关注的重点。传统的遗产本体结构不稳定性监测主要采用两种方式，一种是在遗产目标本体上安装传感器自动获取并记录数据，另一种是通过架设便携式设备，人工周期性获取观测数据。第一类方法特点是基于点目标的观测，其不足之处在于成本高、效率低，耗费大量人力物力，且观测区域有限，尤其是地势险峻、人力不可达地区。因受限于传感器布设位置，该方法空间维采样密度低，无法获取大范围地表形变信息。此外，该方法因需要在文物本体安装设备，属于有损(微损)侵入式监测。第二类方法不足之处在于其时间维采样数低，难以持续开展长时序观测，且观测的空间范围受仪器测距限制。

差分合成孔径雷达干涉(Differential Synthetic Aperture RadarInterferometry，DInSAR)技术通过多次重复观测进行地表形变监测，获取雷达视线方向(line of sight,LOS)上的形变量，在高精度捕获各种地球物理现象引起的地表位移表现出特有能力。传统DInSAR侧重于单次形变或两时刻的累积形变，其因时间、空间去相干和大气扰动的影响，较难获得高精度、时间序列形变信息。长时间序列雷达干涉(Multi-Temporal Interferometric SAR，MT-InSAR)克服了上述局限性，使用观测时间段内所获取的SAR影像集，通过将影像集中的影像按照一定的组合模式进行配对形成干涉图，对干涉相位进行信号分离和提取，从而获得地表形变信息。该方法现已成功应用于多种地表形变监测，包括火山活动、地震灾害、城市地表沉降等所引起的地表形变，但对于古迹遗址形变监测，仍面临如下挑战：1)现存的较多文化遗产地位于植被覆盖较为密集区域，如长城，因植被的体散射失相干严重而使所获取的时序SAR影像相干性低，难以提取足够数量的散射体作为量测目标；2)此外，因长期受自然环境风化侵蚀和人类活动的影响，古迹遗址均存在不同程度的破损，且其材质多为土、石、木类，有别于以混凝土为主的现代建筑，导致能获取的强散射体(persistent scatters，PS)数量有限且空间分布不均匀；3)部分古迹遗址位于大气效应显著区域，如吴哥遗产地所在区域属热带季风气候，采用传统的时空滤波法分离大气相位易引入误差，从而导致形变参数估算错误；4)无法分离大型复杂古建筑高度向混叠信号(叠掩现象)，从而无法提取古建筑在高度向的非均匀形变信息。

由此可见，现有的干涉SAR技术在针对文化遗产地监测方面仍存在技术瓶颈，如何创设一种适用于古迹遗址的精细形变监测方法，成为业界研究的目标。

发明内容

有鉴于此，本发明提供了一种多维层析SAR古迹遗址形变监测方法，以获得古迹遗址高精度的精细形变信息，从而为古迹遗址的预警、修复、保护、监测和管理提供科学依据。

一种多维层析SAR古迹遗址形变监测方法，所述方法包括：

收集待监测古迹遗址的N+1幅SAR影像，对所述SAR影像进行预处理；

采用振幅离散度阈值法和散射强度阈值法提取所述进行预处理后的N+1幅SAR影像中的第一PS候选点，对所述第一PS候选点构建Delaunay三角网，设定大气距离阈值进行弧段筛选，得到初选后的三角网及其弧段所连接的第二PS候选点；

采用单主影像模式对所述进行预处理后的SAR影像进行干涉对组合，进行干涉、去除地形相位和平地相位后，得到N幅差分干涉图；

将所述初选后的三角网的每个弧段所连接的两个第二PS候选点作为邻近点，对所述差分干涉图建立邻近点之间相位差分模型，所述相位差分模型包括：高程残余误差相位模型、线性形变相位模型以及非线性形变相位模型，所述非线性形变相位模型包括基于不同关键驱动因子建立的多元非线性形变相位模型，所述高程残余误差相位模型中的待求解参数为：高程残余误差增量，所述线性形变相位模型中的待求解参数为：年平均线性形变速率增量，所述多元非线性形变相位模型中的待求解参数为：多元非线性形变参数增量；

建立所述第二PS候选点沿多维空间坐标轴的离散化表达模型，所述离散化表达模型包括：采样矩阵与层析信号；所述采样矩阵包括：高程残余误差相位模型、线性形变相位模型以及多元非线性形变相位模型；

采用非参数谱估计法中的波束形成法重构所述层析信号，进行弧段筛选，提取第一高置信度PS点，对所述第一高置信度PS点的高程残余误差增量、年平均线性形变速率增量以及多元非线性形变参数增量进行参数增量集成，得出所述第一高置信度PS点上的高程残余误差、年平均线性形变速率以及多元非线性形变参数；

将所述第一高置信度PS点作为第一参考点；计算所述预处理后的N+1幅SAR图像的平均幅度值，将所述平均幅度值与设定幅度阈值相比较，并将所述平均幅度值大于设定幅度阈值的像元作为PS-MPS候选点；将所述PS-MPS候选点分别与最邻近第一参考点相连，得到星状局域子网，使用波束形成法进行多维层析成像，对所述星状局域子网进行弧段筛选，提取第二高置信度PS点和第一高置信度MPS点；

对所述第二高置信度PS点和所述第一高置信度MPS点进行参数增量集成，得出所述第二高置信度PS点和所述第一高置信度MPS点上的高程残余误差、年平均线性形变速率以及多元非线性形变参数；

对所述第二高置信度PS点和所述第一高置信度MPS点上的高程残余误差、所述年平均线性形变速率以及所述多元非线性形变参数进行地理编码得到所述第二高置信度PS点和所述第一高置信度MPS点上的实际的高程残余误差、实际的年平均线性形变速率以及实际的多元非线性形变参数。

可选地，所述方法还包括：

对所述实际的高程残余误差进行精度验证，对所述实际的线性形变速率以及实际的非线性形变参数进行考察以验证古迹遗址不稳定现象。

可选地，所述方法还包括：

在进行弧段筛选的过程中，采用加权最小二乘迭代法进行参数优化。

可选地，所述采用非参数谱估计法中的波束形成法重构所述层析信号，进行弧段筛选，提取第一高置信度PS点具体包括：

搜索层析信号的最大值；

检测所述层析信号的最大值是否大于第一设定阈值；

如果是，保留与当前层析信号相对应的弧段，并且将所述弧段上的两个第二PS候选点作为第一高置信度PS点，提取第一高置信度PS点。

可选地，所述方法还包括：

对所述第一高置信度PS点进行参数增量集成之前，对所有提取出的所述第一高置信度PS点再次构建Delaunay三角网，使用波束形成法进行多维层析成像，对所述再次构建Delaunay三角网进行弧段筛选，提取第三高置信度PS点，采用加权最小二乘迭代法，进行参数在连续空间的优化；

搜寻最大连通子网，对最大连通子网中的所述第三高置信度PS点进行参数增量集成，得出所述第三高置信度PS点上的高程残余误差、年平均线性形变速率以及多元非线性形变参数；

将所述第三高置信度PS点作为第二参考点；计算所述预处理后的N+1幅SAR图像平均幅度值，将所述平均幅度值大于设定幅度阈值的像元作为PS-MPS候选点；将所述PS-MPS候选点分别与最邻近第二参考点相连，得到星状局域子网，使用波束形成法进行多维层析成像，对所述星状局域子网进行弧段筛选，提取第四高置信度PS点和第二高置信度MPS点；

对所述第四高置信度PS点和所述第二高置信度MPS点进行参数增量集成，得出所述第四高置信度PS点和所述第二高置信度MPS点上的高程残余误差、年平均线性形变速率以及多元非线性形变参数；

对所述第四高置信度PS点、所述第三高置信度PS点和所述第二高置信度MPS点上的高程残余误差、所述年平均线性形变速率以及所述多元非线性形变参数进行地理编码得到所述第四高置信度PS点、所述第三高置信度PS点和所述第二高置信度MPS点上的实际的高程残余误差、实际的年平均线性形变速率以及实际的多元非线性形变参数。可选地，所述相位差分模型为：

2πε_iΔs为高程残余误差相位模型，2πη_iΔv为线性形变相位模型，其中，η_i＝2t_i/λ为时间频率，Δv为年平均线性形变速率增量，Δm_k(k＝1,...,n)为多元非线性形变参数增量，τ_k,i为非线性形变基函数，多元非线性形变参数增量与非线性形变基函数的乘积为多元非线性形变相位模型；Δω_i为残余相位增量。

可选地，所述关键驱动因子包括：因降雨引起的地下水位季相性波动；

所述多元非线性形变相位模型包括：

λ为入射波波长，m₁为地下水位波动引起的地表季相性形变幅度，t_i为第i幅SAR影像时间基线，t₀为地表季相性形变相对地下水位波动的延迟时间，为通过先验知识获得的常数。

可选地，所述关键驱动因子包括：因降雨引起的地下水位季相性波动；

所述多元非线性形变相位模型包括：

λ为入射波波长，m₁、m₂为地下水位波动引起的地表季相性形变幅度，t_i为第i幅SAR影像时间基线，t_a0、t_b0为地表季相性形变相对地下水位波动的延迟时间，为通过先验知识获得的常数。

可选地，所述关键驱动因子还包括：因温度产生的热膨胀形变；

所述多元非线性形变相位模型还包括：

其中，T_i为时间基线，m₃为热膨胀形变幅度。

与现有技术相比，本发明提供的多维层析SAR古迹遗址形变监测方法，至少实现了如下的有益效果：

1)本发明提供的多维层析SAR古迹遗址形变监测方法，将PS-InSAR(Persistentscatterer SAR Interferometry，永久散射体雷达干涉)中相邻的永久散射体通过相位差分模型的参数求解方法引入到多维层析SAR中，通过多维层析成像求解相邻永久散射体的参数增量，从而无需提前估计大气相位，避免了大气相位估计不准确而在形变参数解算中引入的误差，提高了参数估计的可靠性。

2)本发明多维层析SAR古迹遗址形变监测方法，解决植被覆盖密集区域的PS点(persistent scatters，永久散射体)和MPS点(mixed persistent scatters，混叠永久散射体)构网和提取问题，削弱失相干严重和永久散射体空间分布不均匀的影响，扩大形变监测范围，提高自然地表环境下形变监测的能力。

3)本发明多维层析SAR古迹遗址形变监测方法，通过形变模型的建立和验证，来分离文化遗产地的形变驱动因子，从而定量阐明形变成因和时空演变过程。

4)本发明提供的多维层析SAR古迹遗址形变监测方法，通过二次构建Delaunay三角网的方法将大大提高永久散射体点最大连通子网的空间覆盖度和网络连接的稳定性。

当然，实施本发明的任一产品必不特定需要同时达到以上所述的所有技术效果。

通过以下参照附图对本发明的示例性实施例的详细描述，本发明的其它特征及其优点将会变得清楚。

附图说明

被结合在说明书中并构成说明书的一部分的附图示出了本发明的实施例，并且连同其说明一起用于解释本发明的原理。

图1是本发明实施例多维层析SAR古迹遗址形变监测方法的一种流程图；

图2是本发明实施例多维层析SAR古迹遗址形变监测方法的另一种流程图。

具体实施方式

现在将参照附图来详细描述本发明的各种示例性实施例。应注意到：除非另外具体说明，否则在这些实施例中阐述的部件和步骤的相对布置、数字表达式和数值不限制本发明的范围。

以下对至少一个示例性实施例的描述实际上仅仅是说明性的，决不作为对本发明及其应用或使用的任何限制。

对于相关领域普通技术人员已知的技术、方法和设备可能不作详细讨论，但在适当情况下，所述技术、方法和设备应当被视为说明书的一部分。

在这里示出和讨论的所有例子中，任何具体值应被解释为仅仅是示例性的，而不是作为限制。因此，示例性实施例的其它例子可以具有不同的值。

应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步讨论。

实施例一

如图1所示是本发明实施例多维层析SAR古迹遗址形变监测方法的一种流程图，包括以下步骤：

步骤100：收集待监测古迹遗址的N+1幅SAR影像，对所述SAR影像进行预处理。

需要说明的是，对所述SAR影像进行预处理主要是指对SAR影像进行主影像选取、配准、辐射校正处理。

步骤101：采用振幅离散度阈值法和散射强度阈值法提取所述进行预处理后的N+1幅SAR影像中的第一PS候选点。

需要说明的是，PS(persistent scatters，永久散射体)，第一PS候选点是初筛出的永久散射体目标。

步骤102：对所述第一PS候选点构建Delaunay三角网，设定大气距离阈值进行弧段筛选，得到初选后的三角网及其弧段所连接的第二PS候选点。

需要说明的是，第二PS候选点是在所述第一PS候选点的基础上构建Delaunay三角网进行弧段筛选得到的永久散射体目标。进一步，大气距离阈值根据SAR数据以及研究区域的地形地貌、气象特性等因素确定，比如，针对X波段热带雨林地区的SAR数据，大气距离阈值可设为800m。

具体地，Delaunay三角网，又名德洛内三角网，Delaunay三角网是一系列相连的但不重叠的三角形的集合,而且这些三角形的外接圆不包含这个面域的其他任何点。Delaunay三角网具有两个特有的性质:

(1)每个Delaunay三角形的外接圆不包含面内的其他任何点,称之为Delaunay三角网的空外接圆性质,这个特征已经作为创建Delaunay三角网的一项判别标准。

(2)它的另一个性质最大最小角性质:每两个相邻的三角形构成的凸四边形的对角线，在相互交换后，六个内角的最小角不再增大。

进一步，本发明的另一个实施例中，在进行弧段筛选的过程中，可以采用加权最小二乘迭代法进行参数优化。

具体地，采用波束形成法进行层析成像，其参数估算具有正负1/2采样间隔的固有误差。虽然可以通过减小高程残余误差增量、年平均线性形变速率增量以及多元非线性形变参数增量的搜索步长来提高参数估算精度，但代价是将增大计算负荷，严重影响计算效率。然而在连续空间采用加权迭代最小二乘法进行参数优化，可以通过每次迭代计算出的N个干涉对的残余相位来建立约束并修改权值矩阵，修改准则为第i个干涉对的残余相位值越大，该干涉对对应的权值矩阵中的权重则越小。通过多次迭代逐步削弱相位解缠值异常点对整体参数解算结果的影响，从增强了参数解算结果的可靠性和稳健型，以获得参数在连续空间的最优解，从而实现在不显著增加计算负荷的情况下突破固有误差限制。

步骤103：采用单主影像模式对所述进行预处理后的SAR影像进行干涉对组合，进行干涉、去除地形相位和平地相位后，得到N幅差分干涉图。

步骤104：将所述初选后的三角网的每个弧段所连接的两个第二PS候选点作为邻近点，对所述差分干涉图建立邻近点之间相位差分模型，所述相位差分模型包括：高程残余误差相位模型、线性形变相位模型以及非线性形变相位模型，所述非线性形变相位模型包括基于不同关键驱动因子建立的多元非线性形变相位模型，所述高程残余误差相位模型中的待求解参数为：高程残余误差增量，所述线性形变相位模型中的待求解参数为：年平均线性形变速率增量，所述多元非线性形变相位模型中的待求解参数为：多元非线性形变参数增量。

具体地，所述相位差分模型为：

2πε_iΔs为高程残余误差相位模型，2πη_iΔv为线性形变相位模型，其中，η_i＝2t_i/λ为时间频率，Δv为年平均线性形变速率增量，Δm_k(k＝1,...,n)为多元非线性形变参数增量，τ_k,i为非线性形变基函数，多元非线性形变参数增量与非线性形变基函数的乘积为非线性形变相位模型；Δω_i为残余相位增量。需要说明的是，模型中涉及到的n大于等于1，i大于等于1。

具体地，本发明的一个实施例中，所述关键驱动因子包括：因降雨引起的地下水位季相性波动。

所述多元非线性形变相位模型包括：

λ为入射波波长，m1为地下水位波动引起的地表季相性形变幅度，t_i为第i幅SAR影像时间基线，t₀为地表季相性形变相对地下水位波动的延迟时间，为通过先验知识获得的常数。

进一步，在当地水文环境的影响因素较为复杂时，单一的基于正弦或余弦表达的季相形变模型并不能完全表征文化遗产地的非线性形变特征，则可考虑建立双重余弦函数表征文化遗产地的季相形变模型，本发明的另一个实施例中，所述关键驱动因子包括：因降雨引起的地下水位季相性波动。

所述多元非线性形变相位模型包括：

λ为入射波波长，m1、m₂为地下水位波动引起的地表季相性形变幅度，t_i为第i幅SAR影像时间基线，t_a0、t_b0为地表季相性形变相对地下水位波动的延迟时间，为通过先验知识获得的常数。

进一步，通常古建筑多为土、石、木材质，长期暴露于大自然，其结构和力学特性对温度和湿度具有一定的响应，易受气象条件影响产生非均匀形变，针对上述关键驱动因子以及多元非线性形变相位模型，所述关键驱动因子还可以包括：因温度产生的热膨胀形变。所述多元非线性形变相位模型还包括：

其中，T_i为时间基线，m₃为热膨胀形变幅度。

步骤105：建立所述第二PS候选点沿多维空间坐标轴的离散化表达模型，所述离散化表达模型包括：采样矩阵与层析信号；所述采样矩阵包括：高程残余误差相位模型、线性形变相位模型以及多元非线性形变相位模型。

具体地，基于D-TomoSAR(Differential SAR Tomography，差分SAR层析成像)理论，在由N景SAR复数影像组成的数据集中，SAR复数信号g可以看作反射信号沿多维空间坐标轴(斜距高程-线性形变速率-非线性形变幅度)(s-v-m₁-…-m_n)的多重积分，经离散化表达即为g＝Lγ，其中，L为N×M维矩阵，M为在多维空间中沿不同坐标轴方向离散采样的个数；γ为M×1维向量，即待重建的反射信号(层析信号)。假设在(s-v-m₁-…-m_n)多维空间各有个离散采样，那么，L共有个采样值：

因此，参数估计的过程即为重构层析基γ的逆问题

由于D-TomoSAR可以分离混叠信号中的多个散射体，而包含混叠信号的像元相位稳定性已无法仅依赖幅度离散度来判断，因此需考虑同时利用幅度和相位信息来重构层析基。鉴于非参数法谱估计法的鲁棒性好，能稳健提取PS点和MPS点，本课题拟采用非参数谱估计法中的波束形成法进行参数估计。对于邻近PS候选点的相位差分模型(式5)，对两个PS候选点的复数信号进行共轭相乘，则层析信号表示如下：

其中，(·)^H为共轭转置运算符，||·||₂为2-范数，且层析信号已归一化处理，即为最优解。

步骤106：采用非参数谱估计法中的波束形成法重构所述层析信号，进行弧段筛选，提取第一高置信度PS点。

具体地，所述采用非参数谱估计法中的波束形成法重构所述层析信号，进行弧段筛选，提取第一高置信度PS点具体包括A～C三个步骤：

步骤A.搜索层析信号的最大值；

步骤B.检测所述层析信号的最大值是否大于第一设定阈值；如果是，执行步骤C。

需要说明的是，第一设定阈值T_γ1，一般要求大于0.7，即0.7<T_γ1<1。

步骤C.保留与当前层析信号相对应的弧段，并且将所述弧段上的两个第二PS候选点作为第一高置信度PS点，提取第一高置信度PS点。

步骤107：对所述第一高置信度PS点的高程残余误差增量、年平均线性形变速率增量以及多元非线性形变参数增量进行参数增量集成，得出所述第一高置信度PS点上的高程残余误差、年平均线性形变速率以及多元非线性形变参数，将所述第一高置信度PS点作为第一参考点。

步骤108：计算所述预处理后的N+1幅SAR图像的平均幅度值，将所述平均幅度值与设定幅度阈值相比较，并将所述平均幅度值大于设定幅度阈值的像元作为PS-MPS候选点。

需要说明的是，平均幅度值是N+1幅SAR图像的幅度求平均得到的值，而设定幅度阈值由N+1幅SAR图像数据决定，比如，设定幅度阈值为0.033。

步骤109：将PS-MPS候选点分别与最邻近第一参考点相连，得到星状局域子网，使用波束形成法进行多维层析成像，对所述星状局域子网进行弧段筛选，提取第二高置信度PS点和第一高置信度MPS点

步骤110：对所述第二高置信度PS点和第一高置信度MPS点进行参数增量集成，得出所述第二高置信度PS点和所述第一高置信度MPS点上的高程残余误差、年平均线性形变速率以及多元非线性形变参数。

步骤111：对所述第二高置信度PS点和所述第一高置信度MPS点上的高程残余误差、所述年平均线性形变速率以及所述多元非线性形变参数进行地理编码得到所述第二高置信度PS点和所述第一高置信度MPS点上的实际的高程残余误差、实际的年平均线性形变速率以及实际的多元非线性形变参数。

本发明实施例提供的多维层析SAR古迹遗址形变监测方法，将PS-InSAR(Persistent scatterer SAR Interferometry，永久散射体雷达干涉)中相邻的永久散射体通过相位差分模型的参数求解方法引入到多维层析SAR中，通过多维层析成像求解相邻永久散射体的参数增量，从而无需提前估计大气相位，避免了大气相位估计不准确而在形变参数解算中引入的误差，提高了参数估计的可靠性。进一步，本发明解决了植被覆盖密集区域的PSs和MPSs构网和提取问题，削弱失相干严重和永久散射体空间分布不均匀的影响，扩大形变监测范围，提高自然地表环境下形变监测的能力。

实施例二

进一步，为了更好地对本实施例的监测方法进行验证，本发明实施例提供的多维层析SAR古迹遗址形变监测方法，具体包括以下步骤：

步骤200：收集待监测古迹遗址的N+1幅SAR影像，对所述SAR影像进行预处理。

步骤201：采用振幅离散度阈值法和散射强度阈值法提取所述进行预处理后的N+1幅SAR影像中的第一PS候选点。

步骤202：对所述第一PS候选点构建Delaunay三角网，设定大气距离阈值进行弧段筛选，得到初选后的三角网及其弧段所连接的第二PS候选点。

步骤203：采用单主影像模式对所述进行预处理后的SAR影像进行干涉对组合，进行干涉、去除地形相位和平地相位后，得到N幅差分干涉图。

步骤204：将所述初选后的三角网的每个弧段所连接的两个第二PS候选点作为邻近点，对所述差分干涉图建立邻近点之间相位差分模型，所述相位差分模型包括：高程残余误差相位模型、线性形变相位模型以及非线性形变相位模型，所述非线性形变相位模型包括基于不同关键驱动因子建立的多元非线性形变相位模型，所述高程残余误差相位模型中的待求解参数为：高程残余误差增量，所述线性形变相位模型中的待求解参数为：年平均线性形变速率增量，所述多元非线性形变相位模型中的待求解参数为：多元非线性形变参数增量。

步骤205：建立所述第二PS候选点沿多维空间坐标轴的离散化表达模型，所述离散化表达模型包括：采样矩阵与层析信号；所述采样矩阵包括：高程残余误差相位模型、线性形变相位模型以及多元非线性形变相位模型。

步骤206：采用非参数谱估计法中的波束形成法重构所述层析信号，进行弧段筛选，提取第一高置信度PS点。

步骤207：对所述第一高置信度PS点的高程残余误差增量、年平均线性形变速率增量以及多元非线性形变参数增量进行参数增量集成，得出所述第一高置信度PS点上的高程残余误差、年平均线性形变速率以及多元非线性形变参数，将第一高置信度PS点作为第一参考点。

步骤208：计算所述预处理后的N+1幅SAR图像的平均幅度值，将平均幅度值与设定幅度阈值相比较，并将所述平均幅度值大于设定幅度阈值的像元作为PS-MPS候选点。

步骤209：将PS-MPS候选点分别与最邻近第一参考点相连，得到星状局域子网，使用波束形成法进行多维层析成像，对所述星状局域子网进行弧段筛选，提取第二高置信度PS点和第一高置信度MPS点

步骤210：对所述第二高置信度PS点和第一高置信度MPS点进行参数增量集成，得出所述第二高置信度PS点和所述第一高置信度MPS点上的高程残余误差、年平均线性形变速率以及多元非线性形变参数。

步骤211：对所述第二高置信度PS点和所述第一高置信度MPS点上的高程残余误差、所述年平均线性形变速率以及所述多元非线性形变参数进行地理编码得到所述第二高置信度PS点和所述第一高置信度MPS点上的实际的高程残余误差、实际的年平均线性形变速率以及实际的多元非线性形变参数。

步骤212：对所述第二高置信度PS点和所述第一高置信度MPS点上的实际的高程残余误差进行精度验证。

步骤213：对所述第二高置信度PS点和所述第一高置信度MPS点上的实际的线性形变速率以及实际的非线性形变参数进行考察以验证古迹遗址不稳定现象。

具体地，在离散化表达模型中，高程误差与形变参数同时进行解算，因此，修正后的相对高程精度亦可间接反映形变反演的精度。将高程误差修正后的相对高程值与通过实地测量或其他手段获取的建筑物高程值进行对比，分析参数估算精度。此外，通过古迹遗址野外实地调研，对存在明显趋势性形变的散射体对应位置进行现场考察，验证是否存在不稳定现象或趋势。

本发明实施例提供的多维层析SAR古迹遗址形变监测方法，高程残余误差与形变参数(线性形变速率与非线性形变参数)同时进行解算，因此，修正后的相对高程精度亦可间接反映形变反演的精度；将高程误差修正后的相对高程值与通过实地测量或其他手段(如机载LiDAR数据)获取的建筑物高程值进行对比，便于分析参数估算精度。

实施例三

进一步，为了进一步削弱永久散射体空间分布不均匀的影响，如图2所示是本发明实施例多维层析SAR古迹遗址形变监测方法的另一种流程图，包括以下步骤：

步骤300：收集待监测古迹遗址的N+1幅SAR影像，对所述SAR影像进行预处理。

步骤301：采用振幅离散度阈值法和散射强度阈值法提取所述进行预处理后的N+1幅SAR影像中的第一PS候选点。

步骤302：对所述第一PS候选点构建Delaunay三角网，设定大气距离阈值进行弧段筛选，得到初选后的三角网及其弧段所连接的第二PS候选点。

步骤303：采用单主影像模式对所述进行预处理后的SAR影像进行干涉对组合，进行干涉、去除地形相位和平地相位后，得到N幅差分干涉图。

步骤304：将所述初选后的三角网的每个弧段所连接的两个第二PS候选点作为邻近点，对所述差分干涉图建立邻近点之间相位差分模型，所述相位差分模型包括：高程残余误差相位模型、线性形变相位模型以及非线性形变相位模型，所述非线性形变相位模型包括基于不同关键驱动因子建立的多元非线性形变相位模型，所述高程残余误差相位模型中的待求解参数为：高程残余误差增量，所述线性形变相位模型中的待求解参数为：年平均线性形变速率增量，所述多元非线性形变相位模型中的待求解参数为：多元非线性形变参数增量。

步骤305：建立所述第二PS候选点沿多维空间坐标轴的离散化表达模型，所述离散化表达模型包括：采样矩阵与层析信号；所述采样矩阵包括：高程残余误差相位模型、线性形变相位模型以及多元非线性形变相位模型

步骤306：采用非参数谱估计法中的波束形成法重构所述层析信号，进行弧段筛选，提取第一高置信度PS点。

步骤307：对所有提取出的所述第一高置信度PS点再次构建Delaunay三角网。

步骤308：使用波束形成法进行多维层析成像，对所述再次构建Delaunay三角网进行弧段筛选，提取第三高置信度PS点，采用加权最小二乘迭代法，进行参数在连续空间的优化。

具体地，利用幅度阈值法建立规则进一步提取PS候选点，并与上一步中提取出的第一高置信度PS点采用最邻近原则相连，构建局域子网；对于每个局域网的各连接弧段，同样采用波束形成法进行层析成像设置第二设定阈值并且一般大于0.7，若则保留该弧段，且认为该弧段连接的候选点为第三高置信度PS点。

步骤309：搜寻最大连通子网，对最大连通子网中的所述第三高置信度PS点进行参数增量集成，得出所述第三高置信度PS点上的高程残余误差、年平均线性形变速率以及多元非线性形变参数，将所述第三高置信度PS点作为第二参考点。

具体地，对最大连通网络中的第三高置信度PS点进行参数增量集成，即在得到邻近第三高置信度PS点间的形变参数增量(Δv，Δm_k(k＝1,...,n))和高程残余误差增量Δs后，则需要通过增量集成算法，求解对于选定参考点，其余各PS点上形变参数(v，m_k(k＝1,...,n))和s的绝对值。假设经二次构网筛选后的最大连通子网中，弧段数为P，第三高置信度PS点个数为Q，则可将增量集成表达为Y＝A·X，其中，Y表征参数增量矩阵，X表征待求解的参数绝对值，定义A为P×Q维调整矩阵，将Q个第三高置信度PS点排序后，每个弧段对应矩阵A中的一行元素，将该弧段起始第三高置信度PS点在该行中对应元素的值记为-1，将弧段结束第三高置信度PS点在该行中的对应元素的值记为1，其余元素值为0，则可通过矩阵运算方程式求解X。

步骤310：计算所述预处理后的N+1幅SAR图像的平均幅度值，将平均幅度值与设定幅度阈值相比较，并将所述平均幅度值大于设定幅度阈值的像元作为PS-MPS候选点。

步骤311：将所述PS-MPS候选点分别与最邻近第二参考点相连，得到星状局域子网。

步骤312：使用波束形成法进行多维层析成像，对所述星状局域子网进行弧段筛选，提取第四高置信度PS点和第二高置信度MPS点。

本发明的实施例，设置第三设定阈值并且一般大于0.7，将满足的弧段所连接的候选点作为PS-MPS候选点，根据MPS点层析信号幅度应具有的多峰效应特点判定该候选点是否为MPS点。

步骤313：对所述第四高置信度PS点和第二高置信度MPS点进行参数增量集成，得出第四高置信度PS点和第二高置信度MPS点上的高程残余误差、年平均线性形变速率以及多元非线性形变参数。

步骤314：对所述第四高置信度PS点、第三高置信度PS点和第二高置信度MPS点上的高程残余误差、所述年平均线性形变速率以及所述多元非线性形变参数进行地理编码得到所述第四高置信度PS点、第三高置信度PS点和第二高置信度MPS点上的实际的高程残余误差、实际的年平均线性形变速率以及实际的多元非线性形变参数。

步骤315：对所述第四高置信度PS点、第三高置信度PS点和第二高置信度MPS点上的实际的高程残余误差进行精度验证。

步骤316：对所述第四高置信度PS点、第三高置信度PS点和第二高置信度MPS点上的实际的线性形变速率以及实际的非线性形变参数进行考察以验证古迹遗址不稳定现象。

本发明实施例提供的多维层析SAR古迹遗址形变监测方法，通过二次构建Delaunay三角网的方法将大大提高永久散射体点最大连通子网的空间覆盖度和网络连接的稳定性。解决了植被覆盖密集区域的PSs和MPSs构网和提取问题，削弱失相干严重和永久散射体空间分布不均匀的影响，扩大形变监测范围，提高自然地表环境下形变监测的能力。

通过上述实施例可知，本发明提供的多维层析SAR古迹遗址形变监测方法，至少实现了如下的有益效果：

1)本发明提供的多维层析SAR古迹遗址形变监测方法，将PS-InSAR(Persistentscatterer SAR Interferometry，永久散射体雷达干涉)中相邻的永久散射体通过相位差分模型的参数求解方法引入到多维层析SAR中，通过多维层析成像求解相邻永久散射体的参数增量，从而无需提前估计大气相位，避免了大气相位估计不准确而在形变参数解算中引入的误差，提高了参数估计的可靠性。

2)本发明多维层析SAR古迹遗址形变监测方法，解决植被覆盖密集区域的PSs和MPSs构网和提取问题，削弱失相干严重和永久散射体空间分布不均匀的影响，扩大形变监测范围，提高自然地表环境下形变监测的能力。

3)本发明多维层析SAR古迹遗址形变监测方法，通过形变模型的建立和验证，来分离文化遗产地的形变驱动因子，从而定量阐明形变成因和时空演变过程。

4)本发明提供的多维层析SAR古迹遗址形变监测方法，通过二次构建Delaunay三角网的方法将大大提高永久散射体点最大连通子网的空间覆盖度和网络连接的稳定性。

虽然已经通过例子对本发明的一些特定实施例进行了详细说明，但是本领域的技术人员应该理解，以上例子仅是为了进行说明，而不是为了限制本发明的范围。本领域的技术人员应该理解，可在不脱离本发明的范围和精神的情况下，对以上实施例进行修改。本发明的范围由所附权利要求来限定。