面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法

摘要

本发明提供了一种面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法，包括：选择合适的时序DS‑InSAR技术分别对每个平台的多时域数据根据时空基线进行干涉对组合和单视干涉处理，引入外部DEM数据来削弱区域地形引起的密集条纹，对得到的差分干涉图进行nonLocal滤波和相干性重估计，得到最佳的参数估计值；利用地形残差估计值纠正每个分布式散射体的编码高程值，获取所有DS的三维大地坐标，利用ICP技术对多个平台DS点云进行精确配准融合。本发明通过对多源SAR数据进行精确融合，减少了参数估计过程的观测值误匹配率，提高了监测的精度和准确性，结合少量地面控制点，可以获得研究区三维立体增强信息。

权利要求

1.一种面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法，其特征在于，包括：步骤1，选择合适的时序DS-InSAR技术分别对每个平台的多时域数据根据时空基线进行干涉对组合和单视干涉处理，引入外部DEM数据来削弱区域地形引起的密集条纹，对得到的差分干涉图进行nonLocal滤波和相干性重估计，得到最佳的参数估计值；

步骤2，利用地形残差估计值纠正每个分布式散射体(DS)的编码高程值，获取所有DS的三维大地坐标，利用迭代最近点(ICP)算法对多个平台DS点云进行精确配准融合；

步骤3，对融合后的DS点云，在局部球体空间内查找每个点的观测值，选择合适的加权应力应变模型(SM)去求解形变参数，得到融合多平台、多时域SAR影像的非城市研究区的垂直、东西、南北方向的精确三维形变。

2.根据权利要求1所述的面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法，其特征在于，所述步骤1包括：根据获取的覆盖研究区域的多时域的N幅SAR影像，在充分考虑N幅SAR影像之间的时间、空间基线和多普勒质心频率差最优的基础上，根据建立的综合相关函数选取相关系数最大时所对应的影像作为公共主影像，将其余的影像作为从影像分别与主影像进行配准；

根据设定的时空基线阈值对已经配准后的多时域SAR影像自由组合出多个多主影像的短基线干涉对，并对所有干涉对进行常规的干涉处理，接着引入外部DEM数据来削弱区域地形引起的密集条纹，最后获得M幅差分干涉图；

在M幅干涉图中选择干涉质量高的影像，并选取nonLocal滤波方法对干涉结果进行预处理。

3.根据权利要求1所述的面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法，其特征在于，所述步骤2具体包括：利用DS-InSAR获取地形残差以及编码坐标，将二维SAR影像坐标扩展到三维大地坐标；

对每个平台的每个DS点，选择一个固定的球半径r，利用kd-tree快速计算该平台所有点与当前DS点的距离，并统计距离小于半径r的近邻点数，去掉数量少于阈值的噪声点；

将一个指定点与其k邻域内的几何特征信息以统计分布的形式体现在直方图中，比较两个点云数据中点的点特征直方图(PFH)的相似度来确立点对关系，并利用采样一致性原理(RANSAC)对贪婪的初始配准方法进行优化以达到两点云粗配准；

基于ICP算法，采用迭代的方式快速收敛逼近最优的结果，通过估计出某源点云到另一点云的转换矩阵，根据变换矩阵迭代求解使目标函数最小的变换关系，获取最佳配准效果的点云数据。

4.根据权利要求1所述的面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法，其特征在于，所述步骤3包括：对精确融合后的多源DS点云，筛选抛弃非重叠区域点；

为每个DS点选择一个半径为r的局部球体，球体半径根据地面同质区域面积而定，兼顾球体内DS点的同质性和数量，为参数解算提供观测值；

选择合适的先验定权方法，利用点间距离加权，为不同平台观测值提供正确的先验权值；

利用SM模型求解每个DS目标的三维形变参数，充分考虑选取的观测值的邻域相关特性，保证参数的准确性和精度。

说明书

面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法

技术领域

[0001] 本发明涉及基于DS-InSAR技术的精细三维形变监测技术领域，特别涉及一种面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法。

背景技术

[0002] 合成孔径雷达干涉测量技术(Interferometric Synthetic Aperture Radar,InSAR)是近几十年发展起来的一种新型空间大地测量方法。凭借着微波遥感全天时、全天候、大范围(几十公里到几百公里)监测的特点以及SAR卫星影像高精度(厘米到毫米级)和高空间分辨率(几十米到亚米级)的优势，InSAR技术已经越来越得到专家学者的认可，并被广泛用于监测火山、地震运动、山体滑坡、板块运动、冰川漂移、以及矿山开采、地下水抽取和填海等引起的各种地表形变。

[0003] 近几年来，高分辨率(米级至亚米级)SAR卫星(TerraSAR-X/TanDEM-X、COSMO_SkyMed、ALOS-2等)的发射，使得InSAR技术应用范围愈加的广泛，像针对单个建筑物的形变和应力状况监测，制作高分辨率DEM，精细三维建模等等，这些都是中低分辨率SAR卫星数据实现不了的。SAR卫星的升级为我们带来了巨大的便利，但是也对传统的SAR数据处理技术提出了挑战。众所周知，地表形变通常是发生在三维空间框架下的，但是，由于SAR侧视成像特点，单个平台SAR数据只能获取雷达视线方向(LOS)上的地表形变量。因此，在InSAR形变监测中，通常需要融合多源观测值来获取目标三维形变。常规的多源SAR数据融合是将不同平台影像结果地理编码后，重采样到统一坐标系下获取同名观测量。一方面，这种融合方法并不能用于点目标InSAR技术(如PSI、SBAS等)，另一方面，对于高分辨率数据而言，不同成像轨道的差异，参考点参数误差，传播延迟等都会对编码结果产生影响，最终降低多维参数反演的精度，有必要研究更精确的多源数据融合方法。因此，本研究主要目标之一是实现多源高分辨率SAR影像的精确融合。

[0004] 分布式雷达目标InSAR技术是目前InSAR形变监测领域的主流方向，常用的基于分布式散射体(DS)的时序InSAR技术主要有小基线集算法(SBAS)、时域相干目标算法(TCP-InSAR)、相干点目标分析(IPTA)以及SqueeSAR等。相较于永久散射体(PS)点而言，DS点数量更多，应用范围更广泛，在非城市区域同样可以获得丰富的点密度信息。因此，利用基于多源DS-InSAR技术可以监测许多复杂地形区的特定不良地质体的精细多维地表形变。结合以上各种方法的优点，本研究拟将常规二维SAR影像融合推广至三维，发展基于DS-InSAR的多源3D点云数据融合技术，并结合一种新颖的多维参数反演方法，来获取研究区精确三维形变。

发明内容

[0005] 本发明提供了一种面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法，其目的是为了减少参数估计过程的观测值误匹配率，提高三维形变监测的精度和准确性。

[0006] 为了达到上述目的，本发明的实施例提供了一种面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法，包括：

[0007] 步骤1，选择合适的时序DS-InSAR技术分别对每个平台的多时域数据根据时空基线进行干涉对组合和单视干涉处理，引入外部DEM数据来削弱区域地形引起的密集条纹，对得到的差分干涉图进行nonLocal滤波和相干性重估计，得到最佳的参数估计值；

[0008] 步骤2，利用地形残差估计值纠正每个分布式散射体(DS)的编码高程值，获取所有DS的三维大地坐标，利用ICP技术对多个平台DS点云进行精确配准融合；

[0009] 步骤3，对融合后的DS点云，在局部球体空间内查找每个点的观测值，选择合适的加权SM模型去求解形变参数，得到融合多平台、多时域SAR影像的非城市研究区的垂直、东西、南北方向的精确三维形变。

[0010] 其中，所述步骤1包括：

[0011] 根据获取的覆盖研究区域的多时域的N幅SAR影像，在充分考虑N幅SAR影像之间的时间、空间基线和多普勒质心频率差最优的基础上，根据建立的综合相关函数选取相关系数最大时所对应的影像作为公共主影像，将其余的影像作为从影像分别与主影像进行配准；

[0012] 根据设定的时空基线阈值对已经配准后的多时域SAR影像自由组合出多个多主影像的短基线干涉对，并对所有干涉对进行常规的干涉处理，接着引入外部DEM数据来削弱区域地形引起的密集条纹，最后获得M幅差分干涉图；

[0013] 在M幅干涉图中选择干涉质量高的影像，并选取合适的滤波方法对干涉结果进行预处理。

[0014] 其中，所述步骤2具体包括：

[0015] 利用DS-InSAR获取地形残差以及编码坐标，将二维SAR影像坐标扩展到三维大地坐标；

[0016] 对每个平台的每个DS点，选择一个固定的球半径r，利用kd-tree快速计算该平台所有点与当前DS点的距离，并统计距离小于半径r的近邻点数，去掉数量少于阈值的噪声点；

[0017] 将一个指定点与其k邻域内的几何特征信息以统计分布的形式体现在直方图中，比较两个点云数据中点的PFH的相似度来确立点对关系，并利用采样一致性原理(RANSAC)对贪婪的初始配准方法进行优化以达到两点云粗配准；

[0018] 基于ICP算法，采用迭代的方式快速收敛逼近最优的结果，通过估计出某源点云到另一点云的转换矩阵，根据变换矩阵迭代求解使目标函数最小的变换关系，获取最佳配准效果的点云数据。

[0019] 其中，所述步骤3包括：

[0020] 对精确融合后的多源DS点云，筛选抛弃非重叠区域点；

[0021] 为每个DS点选择一个半径为r的局部球体，球体半径根据地面同质区域面积而定，兼顾球体内DS点的同质性和数量，为参数解算提供观测值；

[0022] 选择合适的先验定权方法，利用点间距离加权，为不同平台观测值提供正确的先验权值；

[0023] 利用SM模型求解每个DS目标的三维形变参数，充分考虑选取的观测值的领域相关特性，保证参数的准确性和精度。

[0024] 本发明的上述方案有如下的有益效果：

[0025] 本发明所述的面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法利用InSAR像素分析和3D点云建模的相似性，将传统二维SAR影像融合扩展到了三维，并利用成熟DS-InSAR技术扩展SAR点云融合应用范围。通过对多源SAR数据进行精确融合，减少了参数估计过程的观测值误匹配率，提高了监测的精度和准确性，结合少量地面控制点，还可以获得研究区三维立体增强信息。在提高融合精度的前提下，结合弹性理论中的应力应变模型，利用局部空间里的多个观测值可以获取地面目标的精确三维形变。

附图说明

[0026] 图1是单点的点特征直方图影响区域。

[0027] 图2是两邻域点的UVW局部坐标系示意图。

[0028] 图3是基于多平台SAR影像的三维形变求解关系。

[0029] 图4是基于同质区域点的三维形变反演图。

[0030] 图5是基于DS点云精确配准的多维形变监测流程。

具体实施方式

[0031] 为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述。

[0032] 本发明针对现有的问题，提供了一种面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法。

[0033] 如图5所示，本发明的实施例提供了一种面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法，包括：

[0034] 步骤1，选择合适的时序DS-InSAR技术分别对每个平台的多时域数据根据时空基线进行干涉对组合和单视干涉处理，引入外部DEM数据来削弱区域地形引起的密集条纹，对得到的差分干涉图进行nonLocal滤波和相干性重估计，得到最佳的参数估计值；

[0035] 步骤2，利用地形残差估计值纠正每个分布式散射体(DS)的编码高程值，获取所有DS的三维大地坐标，利用ICP技术对多个平台DS点云进行精确配准融合；

[0036] 步骤3，对融合后的DS点云，在局部球体空间内查找每个点的观测值，选择合适的加权SM模型去求解形变参数，得到融合多平台、多时域SAR影像的非城市研究区的垂直、东西、南北方向的精确三维形变。

[0037] 本发明的上述实施例所述的面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法利用InSAR像素分析和3D点云建模的相似性，将传统二维SAR影像融合扩展到了三维，并利用成熟DS-InSAR技术扩展SAR点云融合应用范围。通过对多源SAR数据进行精确融合，减少了参数估计过程的观测值误匹配率，提高了监测的精度和准确性，结合少量地面控制点，还可以获得研究区三维立体增强信息。在提高融合精度的前提下，结合弹性理论中的应力应变模型，利用局部空间里的多个观测值可以获取地面目标的精确三维形变。

[0038] 其基本思路可以概括为：

[0039] (1)对获取的多平台的数据分别进行单视DS-InSAR时序处理，如干涉、去除平地效应以及去除轨道残留趋势相位、解缠、强散射体选择，得到高相干点的解缠干涉相位。但在相位解缠之前，相位噪声不可避免地会增加干涉图中相位残差的数量(尤其对于单视影像)，最终影响解缠质量。残差是由时间去相关、基线去相关、欠采样、配准误差等多种去相关因素引起的。这些模糊去相关因素因此使得相位噪声滤波成为一项具有挑战性的任务，不采取措施将会影响形变估计准确性以及数字高程模型的生成精度。本研究选择对干涉图进行非局部滤波处理，去掉因大气，噪声，以及各种去相关因子引起的噪声相位，改善干涉图质量，并用于精确估计相干，提高干涉相位相干性，保证DS点密度。

[0040] (2)基于点特征直方图(PFH)的SAR三维点云粗配准。本研究拟采用PFH方法进行SAR点云粗配准，PFH特征描述算子是一种新的高维特征描述算子，它通过对某目标点云特征点及其邻域点之间的空间几何特征进行参数化描述，将一个指定点与其k邻域内的其他点形成的局部空间用统计分布的形式体现在直方图中，然后比较两个点云数据的PFH相似度来确立点对应关系，这种方法无需手动选同名点，相较于最近点选点方法，具有更高的准确度。找出点对应关系后进行粗配准，我们选择基于采样一致性原理(RANSAC)的点云初配准方法。RANSAC算法是一种随机参数估计算法，它与最小二乘这类普通参数估计方法的区别之处在于，它不会直接利用模型拟合全体数据，而是先随机采集能代表点集的部分数据去模拟出一个模型，再利用其余部分数据来完成该模型的检验，该方法的优势在于极大程度地减轻了异常干扰数据对参数模拟的影响。贪婪的初始配准方法(SAC-IA，RANSAC-InitialAlignment)是对点云初始变换矩阵进行估计的一种常用方法，它利用了点云内部数据的刚体旋转不变特性，是一个通过采样并得出最优解的过程。

[0041] (3)基于迭代最近点(ICP)算法的SAR三维点云精配准。基于PFH的三维点云初始配准算法，是把两个不同视角所采集到的三维点云数据通过刚体变换矩阵统一到一个坐标系下，这种初始配准只能使两个点云大致对齐,一般不能达到配准精度的要求，只能作为三维点云配准过程的第一步。因此需要一个精度更高的配准算法，对两个经过初始配准的点云数据进行微调，以进一步提高点云配准的精度。ICP算法以其稳定性和高效性，在三维点云精配准中被广泛的应用和改进，具有很强的鲁棒性。ICP算法通过利用两点云中同名点估计出源点云到另一点云的转换矩阵，然后根据转换矩阵迭代求解使目标函数最小的变换关系，利用该变换关系进行点云精配准。在利用ICP配准SAR点云过程中，我们任然可以利用PFH方法选取同名点，并利用RANSAC随机采样数据，一旦配准达到了预先设置的标准，就认为这两个点云是精确配准的，这样便极大地减弱了地理编码误差对三维形变参数反演的影响。

[0042] (4)基于应力应变模型(SM)的多维形变分解。使用多平台、多时域的InSAR观测值进行物体多维形变监测，其实质是利用测量中的交会原理来获取地表的三维变化，如图3所示。然而，我们都知道从不同平台上不可能获取完全意义上的同名观测值，这和不同卫星空间分辨率，成像几何不同有关，像地面上一个细小目标，在升降轨上成像的可能是该目标的不同部位，但我们同样将其看成同一地物观测值来计算其多维形变。尤其在三维点云空间中，还可能会受到高度的影响，以及配准精度的限制。因此，常规的以像素为基础，用最小二乘法来分解形变的方法在这里就不太适用。我们采用空间局部球体来寻找我们需要的同名点观测信息，这个空间局部球体的半径可以随地面同质区域面积的改变而改变。球体在每个点上滑动，直到覆盖整个感兴趣区域，而在球体内部，如图4所示，则包含了用于后续形变估计的各类观测值信息。在以往研究中，三维变形通常在逐像素的基础上解决，因此忽略了相邻点之间的地面变形的空间相关性。然而，根据弹性变形理论，这也是不合理的。为了利用相邻点变形的空间相关性，有学者提出了基于SM的InSAR和GPS测量相结合的综合应变张量估计方法，该方法可以反演三维形变和应变参数，不需要对稀疏的测量数据进行插值。因此，在获取了局部空间里大量观测量的基础上，我们利用SM模型去解算每个点目标的三维形变，既避免了二维影像融合过程的插值和重采样过程，也考虑了各类观测值的空间相关性。

[0043] 根据这一思路，可以得到如图5所示的基于DS点云精确配准的目标多维形变动态监测方法的实现框图。从图中可以看出，整个方法充分应用了高分辨率SAR数据和3D点目标特点，旨在解决高分辨率SAR数据融合中出现的空间几何不一致问题，利用多平台、多时域的LOS向观测值，加入先进的SM模型还可以获取研究区多维形变趋势以及三维增强地形。整个流程结构清晰，具有实现简单、监测精度高、监测范围大、自动化程度高等优点。

[0044] 其中，所述步骤1包括：根据获取的覆盖研究区域的多时域的N幅SAR影像，在充分考虑N幅SAR影像之间的时间、空间基线和多普勒质心频率差最优的基础上，根据建立的综合相关函数选取相关系数最大时所对应的影像作为公共主影像，将其余的影像作为从影像分别与主影像进行配准；根据设定的时空基线阈值对已经配准后的多时域SAR影像自由组合出多个多主影像的短基线干涉对，并对所有干涉对进行常规的干涉处理，接着引入外部DEM数据来削弱区域地形引起的密集条纹，最后获得M幅差分干涉图；在M幅干涉图中选择干涉质量高的影像，并选取合适的滤波方法对干涉结果进行预处理。

[0045] 如图3所示，所述步骤2具体包括：利用DS-InSAR获取地形残差以及编码坐标，将二维SAR影像坐标扩展到三维大地坐标；对每个平台的每个DS点，选择一个固定的球半径r，利用kd-tree快速计算该平台所有点与当前DS点的距离，并统计距离小于半径r的近邻点数，去掉数量少于阈值的噪声点；将一个指定点与其k邻域内的几何特征信息以统计分布的形式体现在直方图中，比较两个点云数据中点的PFH的相似度来确立点对关系，并利用采样一致性原理(RANSAC)对贪婪的初始配准方法进行优化以达到两点云粗配准；基于ICP算法，采用迭代的方式快速收敛逼近最优的结果，通过估计出某源点云到另一点云的转换矩阵，根据变换矩阵迭代求解使目标函数最小的变换关系，获取最佳配准效果的点云数据。

[0046] 其中，所述步骤3包括：对精确融合后的多源DS点云，筛选抛弃非重叠区域点；为每个DS点选择一个半径为r的局部球体，球体半径根据地面同质区域面积而定，兼顾球体内DS点的同质性和数量，为参数解算提供观测值；选择合适的先验定权方法，利用点间距离加权，为不同平台观测值提供正确的先验权值；利用SM模型求解每个DS目标的三维形变参数，充分考虑选取的观测值的邻域相关特性，保证参数的准确性和精度。

[0047] 本发明的具体实施方案包含以下几个步骤：

[0048] (1)SAR数据选取

[0049] 选取覆盖研究区域的TerraSAR-X/TanDEM-X、ALOS-2、COSMO_SkyMed等多平台、多时域的高分辨率SAR数据。TerraSAR-X/TanDEM-X系统是由EADS Astrium公司与德国宇航中心(DLR)联合开发，TerraSAR-X是2007年升空的，TanDEM-X于2010年6月21日成功送入预定轨道。在运行轨道上，两颗卫星将以不到200米的距离同步飞行，采用X波段精确扫描地球表面来获取影像。ALOS-2卫星是日本宇航局于2014年成功送入预定轨道，采用穿透力强的L波段来监测更为广泛的地表变化。高分辨率雷达卫星COSMO-SkyMed是意大利航天局和意大利国防部共同研发的高分辨率雷达卫星星座，该卫星星座由4颗X波段合成孔径雷达(SAR)卫星组成，整个卫星星座的发射任务于2008年底前完成。

[0050] (2)单平台DS-InSAR时序参数反演

[0051] 得到不同平台下的系列SAR数据后，首先要对每个单平台时间序列的SAR数据进行DS-InSAR处理，提取强散射目标的观测序列值。主要包括以下步骤：

[0052] 1)首先根据获取的覆盖研究区域的多时域的N幅SAR影像，在充分考虑N幅SAR影像之间的时间、空间基线和多普勒质心频率差最优的基础上，根据建立的综合相关函数选取相关系数最大时所对应的影像作为公共主影像。将其余的影像作为从影像分别与主影像进行配准。根据设定的时空基线阈值对已经配准后的多时域SAR影像自由组合出多个多主影像的短基线干涉对，并对所有干涉对进行常规的干涉处理，接着引入外部DEM数据来削弱区域地形引起的密集条纹，最后获得M幅差分干涉图。

[0053] 2)在M幅干涉图中选择干涉质量高的影像，并选取合适的滤波方法对干涉结果进行预处理。通常，影像滤波可以极大地降低噪声数量，而代价是很强的分辨率损失。因此，为了保存点状结构和精细结构和纹理我们会选择采用自适应方法去调整局部估计窗口。非局部滤波器(Nonlocal)通过从小图像块之间的相似性导出数据驱动的权重来成功地对影像进行平滑处理，适合应用在高分辨率SAR影像中，从而为后续处理提供最佳参数配置。

[0054] (3)基于多平台、多时域InSAR三维点云精确配准

[0055] 经过单平台DS-InSAR粗略估计出了DS点在各自成像几何中的地形残差，这里将其用于补偿DEM误差以获取每个DS目标点高程值，再通过获取这些点经纬度坐标可以将二维SAR像素坐标扩展到三维大地坐标。对于高分辨率SAR数据点，由于受到不同轨道成像几何，参考点初始参数选择误差，DEM精度误差等因素的影响，不同平台之间的点云会产生一定的偏移，若不对其进行校正，对多平台同名观测值的选择会产生影响，最终影响形变监测精度，因此，这里拟采用ICP技术对InSAR三维空间点云进行精确配准，从而提高融合质量。点云是指在统一空间坐标系下用来表示目标表面特征及其空间分布的数据集，它以离散点的形式记录了物体表面上各类物理参数。利用多平台、多角度融合点云数据可以获取目标三维变化信息，主要融合步骤如下：

[0056] ①多源数据同名点查找

[0057] 1)基于kd-tree的点邻域搜索

[0058] 设三维数据点集P＝{p1,p2,…,pm}，P中与所选点pi,i＝1,2…m欧式距离最近的k个点即为点pi的k近邻。有两种方法可以确定点的k邻域，一种是直接找到点集中k个与待测点最近的邻域元素，另一种方法则是预先设定一个搜索半径r，找到所有位于以r为半径的球体内的点作为待测点的邻域元素。kd-tree是一种二叉索引树结构，主要基于k维空间，它的每一个节点都表示k维空间的点，这是它与常规二叉树不同之处，可以用来检索多属性数据或多维点数据。利用kd-tree来检索点的k邻域可以极大的提高计算效率，因此，该方法在邻域特征提取、对应点的匹配、特征描述子计算等方面都发挥着重要的作用，成为处理大规模散乱点云数据的重要步骤之一。

[0059] 2)点云特征点提取及表面法向量估计

[0060] 特征点就是指通过对点云模型或2D图像、3D曲面上的点按照一定标准所获取的点集，一般来说，特征点的数量相较于原始点云或图像像素来说是大大减小的，且这些点与其他点相比具有稳定性和代表性。一般特征点会位于山脊，山顶，或者边缘等区域，它可以最真实的表示原始点云数据同时不失描述性，在计算中只针对特征点可以极大地减少后续的点云数据处理的时间。当邻域搜索结束后，各平台点云的局部表面特征就可以用特征点的邻域点集的空间几何分布来描述，而最能体现一个曲面的空间几何特征就是曲面的法线，常用计算法线的方法为利用协方差矩阵分析求解单点的表面法向量，即主成分分析法(PCA)。对数据点集P＝{p1,p2,…,pm}中某点pi,其邻域重心px可表示如下：

[0061]

[0062] 式中，k为点pi在点集P中的邻域点个数。

[0063] 法向量可由式(2)所示的协方差矩阵C来求解。点pi与其邻域点集构成的协方差矩阵可以表达为：

[0064]

[0065] 式中，ξ代表每个邻域点的权重值，最简单的可由邻域点与其重心点距离的倒数来表示。矩阵C是一个对称半正定矩阵，反映了点pi的邻域点相对于重心px的离散程度，λ为矩阵C特征值， 为特征值相对应的特征向量，如果0＜λ0＜λ1＜λ2,则最小特征值λ0对应的特征向量 就是对应点pi的法向 (或者 )。由上述方法算出来的是单位法向矢量，该法向量的方向还未确定，需要规定其正负值已达到所有点的法向一致化，这里采用欧式最小生成树(EMST)原理，即常用的几何一致化方法，随机在树节点中挑选某一节点作为初始点，逐渐向周围扩散，使得该节点 与前后相邻节点的法向 满足下述表达式：

[0066]

[0067] 一旦上述表达式成立，则表示前后两点的法向朝向一致；如果不成立，则令即可达到一致化效果。

[0068] 3)计算点云PFH描述算子

[0069] PFH是一种高维特征描述算子，其充分考虑了点云的邻域空间几何属性，在多维直方图所在的高维信息空间中是可度量的，相比较与点的三维坐标相似性度量，具有更高的可信度。PFH所描述的空间点云曲面特征具有不变性，同时适用于采样密度不同或者噪声等级不同的点云数据，具有很好的鲁棒性。

[0070] 附图1表示单平台点云里某目标点pi的影响区域，虚线圆圈表示以点pi为中心，搜索半径为R的领域范围，所有位于圆圈内的点均属于该点的领域点。该目标点pi的PFH就是利用该圆圈里的所有领域点两两之间计算的法向关系来表示的。一般来说，领域内所有点之间的法向之间相互作用关系主要用法向之间的三个夹角偏差参数α、φ、θ来描述，如附图2，以点pi为原点，建立一个坐标轴为UVW的局部坐标系：

[0071]

[0072] 式中，pt和ps分别为点pi领域点集内两点，nt和ns分别为其法向量，||pt-ps||2为pt和ps两点间距离。则参数α、φ、θ可由式(5)-(7)得到：

[0073] α＝cos-1(V·nt)   (5)

[0074]

[0075]

[0076] 式中，d是两点pt和ps之间的几何距离。将每两个领域点之间计算的三参数α、φ、θ作为描述特征，依次放入该点pi的直方图内。该过程首先对各特征参数值范围进行划分，可以将三个参数分成b个子区间，可以得到一个3*b维空间，在该空间内对点pi的领域内所有两点之间的参数进行统计，得到最终直方图作为该点pi的点特征直方图，用以和其他平台点云特征直方图进行比较，找到直方图最相近的两点，作为同名点进行粗配准。

[0077] ②多源点云数据粗配准

[0078] 点云初始配准主要利用多源点云数据内部的刚体旋转不变特性，主要过程为首先从某源点集P中选取n个样本点，对于这n个样本点，在另一源点云Q中找出能够满足相似度条件的点(本文是满足样本点直方图相同或相似的点)，并存入一个列表中，利用随机采样原则从这些点中选择出一些能代表所有匹配点的采样点，根据选择出的对应点来计算两点云数据集的刚体变换参数：

[0079]

[0080] 式中

[0081]

[0082] 如公式(8)所示，点云配准的空间变换就是首先利用旋转矩阵分别对坐标系o-xyz中的三个坐标轴x,y,z旋转一定的角度，使得两平台点云坐标轴指向一致，再利用平移矢量将坐标系原点o平移到源点云坐标系O-XYZ上。坐标转换一般可以利用七参数模型，但点云数据之间不需要进行缩放，可根据旋转和平移两种变换建立点云数据坐标之间的6独立参数模型，该模型包括3个角度旋转参数α、β、γ和3个平移参数x0,y0,z0，利用该模型可建立坐标系o-xyz和基准坐标系O-XYZ的精确转换关系。一般来说，至少要找到3对控制点，才能确定七参数模型的7个变换参数，且3对控制点不能共线才能完成欧式变换矩阵的参数估计，最终完成数据配准。同样，用于配准两组不同角度获取的SAR点云数据的控制点我们要求其不能共线，且尽量均匀分布。然而，通常在两组数据集中很难精准的完成正确对应点的选取，尤其当研究区域为山区或地形复杂区域时，只选用最少的点参与线性解算，解算出来的刚体变换矩阵会存在比较大的误差，因此，为了提髙其计算的精度，最好是利用尽可能多的对应点对参与到计算中来，以此来尽可能地约束变换方程，然后用相关的数学方法进行参数估算。

[0083] 现有的计算刚体矩阵的方法很多，相比较而言，单位四元数法抗噪性更强且更实用，算法具体如下：

[0084] 对于随机采样的两对应点集P＝{p1,p2,…,pn}和Q＝{q1,q2,…,qn}，分别计算两点集质心：

[0085]

[0086] 计算两数据集的互协方差矩阵：

[0087]

[0088] 并求协方差矩阵的循环列向量Δ＝(m23,m31,m12)，其中mij为mij＝(Cpq-CpqT)ij，再根据协方差矩阵构造对称阵：

[0089]

[0090] 式中，trace(Cpq)为协方差矩阵Cpq的迹，I3为3×3的单位阵。计算协方差矩阵C的特征值、特征向量，找出最大特征值对应的特征向量v＝(v1,v2,v3,v4)；最后求解刚体旋转矩阵R和T：

[0091]

[0092] T＝p-Rq   (13)

[0093] 利用计算出的刚体变换矩阵对点云数据进行旋转和平移，以达到两点云粗配准。

[0094] ③基于ICP点云精确配准方法

[0095] 在完成初始配准的条件下，ICP算法在点云数据精配准中有着广泛的应用。但在利用ICP算法进行精细配准之前，同样需要使用随机采样方法对参与配准的初配准点云数据进行精简，从而减少参与运算的数据量，可以加快配准的速度。由于利用最近点搜索同名点往往会造成错误匹配点对的产生，从而对ICP配准产生负面的影响，本研究拟结合SAR影像像素相关性以及PFH特征来选取同名点，以减少点错误配对的情况。ICP算法主要是利用选取的同名点去迭代估计出某源点云到另一点云的最优转换矩阵，可以用一个目标函数表示，如公式(14)，该方法通过迭代的方法对旋转矩阵和平移矩阵进行优化，对于每一次迭代，都在参考点集中寻找目标点集中的对应点，并由此对坐标变换量进行改进，直到误差收敛，进而利用得到的最优变换矩阵实现两点集的精确配准。

[0096]

[0097] 其中，f(R,T)为相应的均方误差，通常会预先设置一定的阈值，作为判断配准好坏的依据，R为旋转变换矩阵，T为平移变换矩阵。pi,qi分别为两粗匹配点云。

[0098] (4)基于SM模型的多维形变监测方法

[0099] 传统的InSAR一般采用逐像素的方法来求解三维形变，从而忽略了相邻点之间的空间相关性，因此，本研究拟考虑利用弹性变形理论中的SM模型来解决这个问题。通过考虑空间相邻元素的观测值，SM模型可以检索更精确的3D形变量以及地表应变参数,并且不需要对稀疏SAR观测值进行内插。利用球体内包含的所有配准后的点(包括不同平台的点)来计算当前点的三维形变，增加了方程观测量的同时，可以获得精确的地面三维形变量。

[0100] 地球表面的一部分由于地球动力学过程(如岩浆侵入或断层活动)而变形的区域可视为均质应变场，假设场中一个以地面位置分量x0＝[xe0,xn0,xu0]为特征的感兴趣点P0，为了求其三维形变矢量d0＝[de0,dn0,du0]，我们在局部空间里选择了其周围k个观测目标i＝1,2…k，获得相应的三维形变矢量di＝[dei,dni,dui]和位置矢量xi＝[xei,xni,xui]。在这个局部空间里，我们可以为当前感兴趣点P0建立如下关系：

[0101] di＝H·Δi+d0   (15)

[0102] 其中，Δi＝[Δxei Δxni Δxui]T表示点P0和Pi之间的坐标增量向量，H是应变参数矩阵，可分为对称E和非对称部分R：

[0103]

[0104] 其中，E和R分别表示应变张量和刚体旋转张量，ξ，ω为SM模型的未知参数。我们利用常规三维分解模型可以建立起观测量L和未知参数的关系：

[0105]

[0106] 式中，l＝[de0,dn0,du0,ξ11,ξ12,ξ13,ξ22,ξ23,ξ33,ω1,ω2,ω3]，d0＝[de0,dn0,du0]需要形变分解的点的三维形变值，ξ、ω代表SM模型的未知参数，Li为第i个点的LOS向形变量，Bsm为SM模型的设计矩阵，Bgeo为常规三维分解关系矩阵。通过计算所有点(不同平台的点)的多维形变量，可以获取更丰富的区域形变信息。对于重叠区域内所有点，都可用上述方法来求解多维形变。

[0107] 如图5所示，首先对每个平台获取的高分辨率SAR数据集分别进行DS-InSAR处理，获取最优的参数(形变、地形残差值)估值。再选择一个合适的参考SAR点云基准，采用迭代的方法将其余目标点云分别和它配准，达到相对位置最优。最后在融合后的点云局部空间里选择观测值，并利用SM模型提取出各个点目标的三维形变矢量。

[0108] 以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明所述原理的前提下，还可以作出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。