CN113159174B - 无控制点的dem匹配与形变量探测方法及装置 - Google Patents
无控制点的dem匹配与形变量探测方法及装置 Download PDFInfo
- Publication number
- CN113159174B CN113159174B CN202110428378.1A CN202110428378A CN113159174B CN 113159174 B CN113159174 B CN 113159174B CN 202110428378 A CN202110428378 A CN 202110428378A CN 113159174 B CN113159174 B CN 113159174B
- Authority
- CN
- China
- Prior art keywords
- dem data
- points
- dem
- deformation
- period
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 86
- 238000012544 monitoring process Methods 0.000 claims abstract description 78
- 238000007619 statistical method Methods 0.000 claims abstract description 35
- 230000000694 effects Effects 0.000 claims abstract description 28
- 238000012216 screening Methods 0.000 claims description 58
- 238000001514 detection method Methods 0.000 claims description 36
- 238000004422 calculation algorithm Methods 0.000 claims description 32
- 238000009616 inductively coupled plasma Methods 0.000 claims description 30
- 238000006243 chemical reaction Methods 0.000 claims description 27
- 239000011159 matrix material Substances 0.000 claims description 25
- 238000012545 processing Methods 0.000 claims description 21
- 230000009466 transformation Effects 0.000 claims description 19
- 238000004590 computer program Methods 0.000 claims description 16
- 230000008030 elimination Effects 0.000 claims description 16
- 238000003379 elimination reaction Methods 0.000 claims description 16
- 230000006870 function Effects 0.000 claims description 15
- 238000004364 calculation method Methods 0.000 claims description 9
- 238000013507 mapping Methods 0.000 claims description 9
- 238000003860 storage Methods 0.000 claims description 7
- 238000000354 decomposition reaction Methods 0.000 claims description 6
- 230000009897 systematic effect Effects 0.000 claims description 6
- 238000009826 distribution Methods 0.000 claims description 5
- 238000010276 construction Methods 0.000 claims description 4
- 238000001914 filtration Methods 0.000 claims description 4
- 238000002156 mixing Methods 0.000 claims description 4
- 238000012360 testing method Methods 0.000 claims description 4
- 238000010586 diagram Methods 0.000 description 20
- 230000008569 process Effects 0.000 description 8
- 230000008859 change Effects 0.000 description 6
- 238000000605 extraction Methods 0.000 description 6
- 238000013519 translation Methods 0.000 description 6
- 230000004044 response Effects 0.000 description 5
- 230000005484 gravity Effects 0.000 description 4
- 230000009286 beneficial effect Effects 0.000 description 3
- 238000011439 discrete element method Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 239000000284 extract Substances 0.000 description 3
- 238000007476 Maximum Likelihood Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 239000013598 vector Substances 0.000 description 2
- 238000009395 breeding Methods 0.000 description 1
- 230000001488 breeding effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000002427 irreversible effect Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000004660 morphological change Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 238000011155 quantitative monitoring Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000001629 suppression Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/74—Image or video pattern matching; Proximity measures in feature spaces
- G06V10/75—Organisation of the matching processes, e.g. simultaneous or sequential comparisons of image or video features; Coarse-fine approaches, e.g. multi-scale approaches; using context analysis; Selection of dictionaries
- G06V10/757—Matching configurations of points or features
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
- G06T3/4023—Scaling of whole images or parts thereof, e.g. expanding or contracting based on decimating pixels or lines of pixels; based on inserting pixels or lines of pixels
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Physics & Mathematics (AREA)
- Physics & Mathematics (AREA)
- Evolutionary Computation (AREA)
- Databases & Information Systems (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Software Systems (AREA)
- Computing Systems (AREA)
- Artificial Intelligence (AREA)
- Health & Medical Sciences (AREA)
- Multimedia (AREA)
- Image Analysis (AREA)
Abstract
本发明公开了一种无控制点的DEM匹配与形变量探测方法及装置，其中该方法包括：获取监测区域的多期无人机航摄影像；根据多期无人机航摄影像，得到监测区域多期DEM数据及对应的RGB图像；通过特征识别方法提取出多期RGB图像间的同名特征点，基于此获得监测区域多期DEM数据间优选的初始位置关系；在获得优选的初始位置关系后，首先利用统计学方法监测出最大概率位于形变区域的点，将其剔除；其次利用距离衰减效应，探测出可能位于形变区域的微小形变点，将其剔除；在剔除干扰形变点后，进行精确配准并获取监测区域精确形变量，本发明可以在难以定向布设控制点或未布设控制点的区域进行形变量精确探测，对于地质灾害监测具有重要意义。
Description
技术领域
本发明涉及地质灾害分析技术领域，尤其涉及一种无控制点的DEM匹配与形变量探测方法及装置。
背景技术
本部分旨在为权利要求书中陈述的本发明实施例提供背景或上下文。此处的描述不因为包括在本部分中就承认是现有技术。
地表形变是一种缓慢变化和不可逆的地质现象，地表形变场的获取与精化是认识地壳运动与自然灾害孕育机理、探求其对环境的影响的重要基础性工作。地表形变场的变化通常具有连续性、动态性、区域性和复杂性等特点，受到地球内部构造运动和人类活动的共同影响，其形变场能直接反映地球表面的动态变化，其变化结果与人类生存环境密切相关。在地表变化过程中，形变量是反映地表当前稳定性及运动状态最直接的物理量，对了解掌握地表变化以及准确预测其发展趋势具有重要意义。
数字高程模型(Digital Elevation Model,DEM)蕴含大量地质地貌、地表特征等信息，高精度、高分辨率的DEM使得对地貌参数精确提取及动态变化监测得以实现，目前已成为地表形态变化定量化监测与分析的关键性数据。目前基于DEM的地表形变监测技术其基本方法都是利用多时相地面监测数据，基于地面控制点建立多期高精度DEM数据，通过多期DEM构建差分模型来监测研究区域形变。
然而，各种监测技术都需要地面控制点对DEM数据进行绝对定向，但是在地表监测过程中，难以保证地面控制点的稳定性与不变形，因此现有对监测区域进行准确形变探测存在一定困难。
发明内容
本发明实施例提供一种无控制点的DEM匹配与形变量探测方法，用以在控制点难以定向布设或未布设控制点的区域进行形变量精确探测，该方法包括：
获取监测区域的多期无人机航摄影像；
根据多期无人机航摄影像，得到监测区域多期DEM数据及对应的RGB图像；
通过特征识别方法提取出多期RGB图像间的同名特征点，在DEM数据和RGB 图像间建立模型点与像素的对应关系，根据对应关系将同名特征点映射至多期DEM 数据中，确定多期DEM数据中的同名特征点间的转换关系，根据转换关系对DEM 数据进行坐标变换，获得监测区域多期DEM数据间优选的初始位置关系；
对获得优选的初始位置关系的DEM数据，通过插值方法建立多期DEM数据的高程对应关系，根据所述高程对应关系，利用统计方法检测出多期DEM数据中位于形变区域的点，将位于形变区域的点剔除，得到初步筛选后的多期DEM数据；
利用距离衰减效应，检测出初步筛选后的多期DEM数据中未被统计方法检测出的位于形变区域的点，将未被统计方法检测出的位于形变区域的点剔除，得到二次筛选后的多期DEM数据；
对二次筛选后的多期DEM数据，采用ICP算法进行整体精匹配，得到精匹配后的多期DEM数据；根据精匹配后的多期DEM数据构建差分模型，得到监测区域形变量；所述形变量用于分析监测区域的滑坡灾害。
本发明实施例还提供一种无控制点的DEM匹配与形变量探测装置，用以在控制点难以定向布设或未布设控制点的区域进行形变量精确探测，该装置包括：
获取单元，用于获取监测区域的多期无人机航摄影像；
DEM数据及RGB图像确定单元，用于根据多期无人机航摄影像，得到监测区域多期DEM数据及对应的RGB图像；
初始位置关系确定单元，用于通过特征识别方法提取出多期RGB图像间的同名特征点，在DEM数据和RGB图像间建立模型点与像素的对应关系，根据对应关系将同名特征点映射至多期DEM数据中，确定多期DEM数据中的同名特征点间的转换关系，根据转换关系对DEM数据进行坐标变换，获得监测区域多期DEM数据间优选的初始位置关系；
初步筛选单元，用于对获得优选的初始位置关系的DEM数据，通过插值方法建立多期DEM数据的高程对应关系，根据所述高程对应关系，利用统计方法检测出多期DEM数据中位于形变区域的点，将位于形变区域的点剔除，得到初步筛选后的多期DEM数据；
二次筛选单元，用于利用距离衰减效应，检测出初步筛选后的多期DEM数据中未被统计方法检测出的位于形变区域的点，将未被统计方法检测出的位于形变区域的点剔除，得到二次筛选后的多期DEM数据；
形变量确定单元，用于对二次筛选后的多期DEM数据，采用ICP算法进行整体精匹配，得到精匹配后的多期DEM数据；根据精匹配后的多期DEM数据构建差分模型，得到监测区域形变量；所述形变量用于分析监测区域的滑坡灾害。
本发明实施例还提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述无控制点的DEM匹配与形变量探测方法。
本发明实施例还提供一种计算机可读存储介质，所述计算机可读存储介质存储有执行上述无控制点的DEM匹配与形变量探测方法的计算机程序。
本发明的有益效果是：本发明实施例提供的无控制点DEM匹配与形变量探测方案通过：采用无人机航摄影像，基于鲁棒的特征提取方法提取出多期DEM数据间的同名特征点，以获得优选的初始位置姿态；在获取优选的初始位置姿态后，通过空间插值方法建立对应点关系，利用统计学方法监测出最大概率位于形变区域的点，将其剔除；同时，充分考虑DEM数据间的空间效应，利用距离衰减效应，探测出可能位于形变区域的微小形变点，将其剔除；在剔除会对精匹配过程造成干扰的形变点后，进行精确配准并获取精确形变量，可分析滑坡灾害的时空形变特征，为灾害早期预报预警提供决策支撑，可在难以定向布设控制点或未布设控制点的区域进行形变量精确探测，对于地质灾害监测具有重要意义。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。在附图中：
图1为本发明实施例中无控制点的DEM匹配与形变量探测方法的整体流程图；
图2为本发明实施例中一幅无人机航摄影像示意图，所述无人机航摄影像为RGB图像；
图3(a)、图3(b)为本发明实施例中DEM数据示意图，其中图3(a)为由第一期航摄影像构建的DEM数据，图3(b)为第二期航摄影像构成的DEM数据；
图4为本发明实施例中具有正确对应关系的特征点对示意图；
图5为本发明实施例中具有初始匹配关系的DEM数据示意图；
图6为本发明实施例中某监测区域形变探测结果图；
图7为本发明实施例中无控制点的DEM匹配与形变量探测方法的流程示意图；
图8为本发明实施例中无控制点的DEM匹配与形变量探测装置的结构示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚明白，下面结合附图对本发明实施例做进一步详细说明。在此，本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。
目前基于DEM的地表形变监测技术其基本方法都是利用多时相地面监测数据，基于地面控制点建立多期高精度DEM数据，通过多期DEM构建差分模型来监测研究区域形变。然而，在实际情况中，尤其是在地形特征较复杂的地区，往往难以布设控制点，故无法进行准确绝对定向。因此，无控制点匹配方法具有重大作用。传统无控制点匹配方法基于最小二乘策略，要求数据具有较好的初始位置，否则匹配成功率大大较低；并且，由于发生形变，导致数据中仅存在部分重叠区域，具有较大的干扰性，使传统无控制点匹配方法求解时容易陷入局部最优，甚至无法求解。
针对现有无控制DEM匹配要求条件高、精度低、成功率低等问题，本发明实施例提供一种无控制点的DEM匹配与形变量探测方案，该方案解决了多期监测数据需要控制点定向的问题，该方案基于无人机高分辨率航摄影像，首先基于DEM数据对应的RGB影像提取同名特征点，获取初始转换参数，为多期DEM数据提供较好的初始姿态；其次，将统计学方法与地理学方法引入粗差剔除中，有效剔除会对精匹配过程造成干扰的形变点；最后，采用ICP算法进行精匹配并通过构建差分模型获取监测区域形变量。下面对该无控制点的DEM匹配与形变量探测方案进行详细介绍。
图7为本发明实施例中无控制点的DEM匹配与形变量探测方法的流程示意图，如图7所示，该方法包括如下步骤：
步骤101：获取监测区域的多期无人机航摄影像；
步骤102：根据多期无人机航摄影像，得到监测区域多期DEM数据及对应的RGB 图像；
步骤103：通过特征识别方法提取出多期RGB图像间的同名特征点，在DEM数据和RGB图像间建立模型点与像素的对应关系，根据对应关系将同名特征点映射至多期DEM数据中，确定多期DEM数据中的同名特征点间的转换关系，根据转换关系对DEM数据进行坐标变换，获得监测区域多期DEM数据间优选的初始位置关系；
步骤104：对获得优选的初始位置关系的DEM数据，通过插值方法建立多期DEM 数据的高程对应关系，根据所述高程对应关系，利用统计方法检测出多期DEM数据中位于形变区域的点，将位于形变区域的点剔除，得到初步筛选后的多期DEM数据；
步骤105：利用距离衰减效应，检测出初步筛选后的多期DEM数据中未被统计方法检测出的位于形变区域的点，将未被统计方法检测出的位于形变区域的点剔除，得到二次筛选后的多期DEM数据；
步骤106：对二次筛选后的多期DEM数据，采用ICP算法进行整体精匹配，得到精匹配后的多期DEM数据；根据精匹配后的多期DEM数据构建差分模型，得到监测区域形变量；所述形变量用于分析监测区域的滑坡灾害。
针对现有技术中的缺陷和不足，本发明提供一种无控制点的DEM匹配与形变量探测方法以解决多期监测数据需要控制点定向的问题。为实现上述目的，本发明实施例提供的无控制点的DEM匹配与形变量探测方法工作时：获取监测目标发生形变前与发生形变后过程中的多期无人机航摄影像；将所述多期无人机航摄影像建立出监测区域DEM数据，并将多期DEM数据导出为RGB图像；对上述DEM数据导出的 RGB图像，通过特征识别方法提取出多幅影像间的同名特征点，在上述DEM数据和 RGB图像间建立模型点与像素的对应关系，根据对应关系将上述同名点映射至DEM 数据中，计算出上述DEM数据中的同名点间的转换关系，并以此转换关系对DEM 数据进行坐标变换，获得上述多期DEM间较好的“初始位置关系”(优选的初始位置关系)；对上述已具有较好初始转换的DEM数据，通过插值方法建立多期DEM数据的高程对应关系，利用统计方法检测出多期DEM数据中位于形变区域的点，将其剔除；对上述DEM数据，利用距离衰减效应，检测出多期DEM数据中未被上述统计方法检测出的位于形变区域的点，将其剔除；对上述剔除形变点的DEM数据，采用ICP算法进行整体精匹配，并将上述精匹配后的多期DEM数据构建差分模型，得到监测区域形变量。
本发明实施例的有益效果是：本发明的无控制点DEM匹配与形变量探测方法采用无人机航摄影像，基于鲁棒的特征提取方法提取出多期DEM数据间的同名特征点，以获得较好的初始位置姿态；在获取初始姿态后，通过空间插值方法建立对应点关系，利用统计学方法监测出最大概率位于形变区域的点，将其剔除；同时，充分考虑DEM 数据间的空间效应，利用距离衰减效应，探测出可能位于形变区域的微小形变点，将其剔除；在剔除会对精匹配过程造成干扰的形变点后，进行精确配准并获取精确形变量，可分析滑坡灾害的时空形变特征，为灾害早期预报预警提供决策支撑，可在难以布设控制点或未布设控制点的区域进行差异探测，对于地质灾害监测具有重要意义。
下面结合图1至图6对本发明实施例涉及的各个步骤进行详细介绍。
如图1所示，本发明提供的无控制DEM匹配与形变量探测方法包括以下步骤：
步骤S1：获取监测目标发生形变前与发生形变后过程中的多期无人机航摄影像，即上述步骤101的详细步骤，以下实例设定以两期无人机航摄影像为例进行介绍。
图2为本发明提供的一幅无人机航摄影像示意图，所述无人机航摄影像为RGB 图像。
步骤S2：将所述多期无人机航摄影像(至少两期无人机航摄影像)建立出监测区域DEM数据，并将多期DEM数据导出为RGB图像，即上述步骤102的详细步骤。
步骤S2.1：将获取的多期监测区域航摄影像进行空三计算(空中三角测量)、密集匹配、点云滤波以及规则格网构建处理，生成监测区域DEM数据。图3(a)、图3(b)为本发明提供的DEM数据示意图，其中图3(a)为由第一期航摄影像构建的 DEM数据，图3(b)为第二期航摄影像构成的DEM数据。
步骤S2.2：将上述DEM数据按原分辨率导出对应的RGB图像，使DEM格网数与RGB图像的像素数相对应。
通过上述可知，在一个实施例中，根据多期无人机航摄影像，得到监测区域多期DEM数据及对应的RGB图像，可以包括：
将监测区域的多期无人机航摄影像进行空三计算、密集匹配、点云滤波及规则格网构建处理，生成监测区域多期DEM数据；
将监测区域多期DEM数据按原分辨率导出对应的RGB图像，使DEM格网数与 RGB图像的像素数相对应。
步骤S3：对上述DEM数据导出的RGB图像，通过特征识别方法提取出多幅影像间的同名地物点，在上述DEM数据和RGB图像间建立模型点与像素的对应关系，根据对应关系将上述同名点映射至DEM数据中，利用奇异值分解法解算出上述DEM 数据中的同名点间的转换关系，并以此转换关系对DEM数据进行坐标变换，获得上述多期DEM间较好的“初始位置关系”，即上述步骤103的详细步骤。
在一个实施例中，在上述步骤103中，通过特征识别方法提取出多期RGB图像间的同名特征点，在DEM数据和RGB图像间建立模型点与像素的对应关系，根据对应关系将同名特征点映射至多期DEM数据中，确定多期DEM数据中的同名特征点间的转换关系，根据转换关系对DEM数据进行坐标变换，获得监测区域多期DEM 数据间优选的初始位置关系，可以包括：
使用SURF算法对多期RGB图像提取同名特征点；
使用RANSAC算法对所述同名特征点进行一致性检测，剔除错误配对的同名特征点，得到剔除处理后的RGB同名特征点；
以DEM数据左上角为坐标原点建立DEM平面直角坐标系，以RGB图像的左上角为坐标原点建立RGB平面直角坐标系，记录剔除处理后的RGB同名特征点在RGB 图像平面坐标系中的位置，根据DEM数据和RGB图像间的矩阵对应关系在DEM平面直角坐标系中选择与剔除处理后的RGB同名特征点相同位置的点作为DEM同名特征点；
使用奇异值分解法求解多期DEM同名特征点间的初始转换关系；
根据初始转换关系对DEM数据进行坐标变换，获得监测区域多期DEM数据间优选的初始位置关系。
下面对该步骤进行详细介绍。
步骤S3.1：使用SURF算法对两期RGB图像提取同名特征点。
步骤S3.1.1提取特征点
采用Hessian矩阵检测上述导出的RGB图像特征点，图像中某点(x,y)在某一尺度σ下的Hessian矩阵定义为：
上式中，Lxx(x,y,σ)代表高斯二阶偏导数
与图像中点(x,y)的像素作卷积，Lxy(x,y,σ)代表
与图像中点(x,y)的像素作卷积，Lyy(x,y,σ)代表
与图像中点(x,y)的像素作卷积。
为了提高运算速度，采用盒式滤波器来近似高斯滤波器，与原始图像进行卷积后以Dxx、Dxy、Dyy近似代替上述Lxx、Lxy、Lyy，形成快速Hessian矩阵，即：
快速Hessian矩阵描述了函数的局部曲率，通过求取矩阵的行列式来判断当前检测点是否为极值点，即：
det(Hfast)＝DxxDyy-(ωDxy)2；
上式中，ω为权值，为了减小使用因使用盒式滤波器近似代替高斯滤波器所带来的误差，此处取值为0.9。当快速Hessian矩阵的行列式值为正值时，表示当前检测点是极大值或极小值，表明该点周围像素变化剧烈，可将其判断为特征点。
步骤S3.1.2特征点的尺度空间定位
构建图像的尺度空间，尺度空间由4组4层组成，同一组间不同层使用相同尺寸的滤波器，但是滤波器的模糊系数逐渐增大；
在不同尺度特征点的响应图像上采用领域非极大值抑制，将每个像素点与二维图像空间和尺度空间领域内的26个点进行比较，选出特征点候选点；
利用三维线性插值法对候选点进行定位，获得亚像素级别的特征点，完成特征点的提取。
步骤S3.1.3确定特征点主方向
以特征点为圆心，6S为半径形成一个圆形区域，计算区域内所有特征点在水平和垂直方向上的Haar小波响应值；
以特征点为圆心，以60°的扇形窗口扫描整个圆形区域，将扇形区域内所有特征点的Haar小波响应值进行累加，以累加值最大的那个方向作为该特征点的主方向。
步骤S3.1.4特征描述符
以特征点为坐标原点，坐标轴旋转至特征点主方向，建立局部坐标系；
在该局部坐标系上构建一个尺寸为20×20的正方形窗口，再对该窗口进行分割，分成4×4个子窗口，即每个子窗口的尺寸为5×5；
计算每个子窗口内的特征点在水平方向Haar小波响应值dx之和、垂直方向Haar小波响应值dy之和以及dx与dy的绝对值之和|dx|、|dy|，因此，每个子区域各具有一个四维向量V＝[∑dx ∑dy ∑|dx| ∑|dy|]；
统计该特征点区域内的所有特征向量V，即得到一个4×4×4＝64维的特征描述符。
步骤S3.1.5对应特征点计算
根据每个特征点的64维特征描述符，两两计算特征点描述符间的欧氏距离，以确定匹配度，欧氏距离最短的特征点，作为多期RGB数据间的对应同名特征点。
步骤S3.2：使用RANSAC检查特征点一致性，剔除错误匹配点。
上述SURF算法可有效提取出多期RGB图像中的对应特征点，然而在实际情况中，往往存在一些错误对应点，使用RANSAC算法对上述同名特征点进行一致性检测，剔除错误配对的特征点，RANSAC算法可有效剔除错误对应点，提高匹配精度 RANSAC通过反复选择数据中的一组随机子集来达成目标。假设上述SURF算法提取出的对应点集中正确的对应点对为“内点”，错误的对应点对为“外点”。
设SURF在第一期RGB图像中找到的特征点集为{Qfeature|(xi,yi),i＝1,2...,n}，其在第二期RGB图像中对应的特征点为{Pfeature|(x'i,y'i),i＝1,2...,n}，n为特征点对的数量。计算步骤：
步骤S3.2.1随机从Qfeature中抽取4个点，并保证4个点不共线，由式
计算出单应性矩阵h，记为模型M，即从数据集中随机选取一组局内点，并计算模型参数M；
上式中，(x,y)为Qfeature点集中的点，(x’,y’)为Pfeature点集中对应的点，h为
步骤S3.2.2用上述计算出的模型参数M去测试其它所有数据点，若各点与模型的误差小于设定的阈值，则将其视为“内点”，否则视为“外点”，即计算点集Qfeature与点集Pfeature中满足模型M的点对数量，若小于设定的误差阈值则将其视为“内点”，计算适应度函数，若此模型越优越，则对应的适应度函数越小，适应度函数为：
步骤S3.2.3将当前模型的适应度函数与最佳模型进行比较，若当前模型更优越，则将当前模型设定为最佳模型；即，如果有足够多的点被归类为“内点”，那么该模型就足够合理，保留目前为止局内点数目最多的模型，将其记录为最佳模型；反之，则返回第一步；
步骤S3.2.4重复上述步骤直至达到设定的迭代次数，选用最佳模型(确定为最佳模型的最终模型)对应的“内点”集作为多期数据的对应特征点集，即重复上述步骤直至达到设定的迭代次数，并选用最佳模型作为最终模型。图4为本发明提供的经筛选后具有正确对应关系的特征点对示意图。
通过上述可知，在一个实施例中，使用RANSAC算法对所述同名特征点进行一致性检测，剔除错误配对的同名特征点，得到剔除处理后的RGB同名特征点，可以包括：
从第一期RGB图像的特征点集中随机选取一组内点，计算模型参数M；
用计算出的模型参数M去测试所有数据点，若各数据点与模型的误差小于设定的阈值，则将其视为内点，否则视为外点；所述内点为正确配对的特征点，外点为错误配对的特征点；
如果有大于预设数目的数据点被归类为内点，确定对应模型为合理模型，在多个合理模型中选出内点数目最多的模型作为最佳模型；
重复上述步骤直至达到设定的迭代次数，将达到设定的迭代次数时选出的最佳模型作为最终模型，选用最终模型对应的内点集作为剔除处理后的RGB同名特征点。
步骤S3.3：对DEM数据与RGB图像，以数据(图像)左上角为坐标原点，设横轴为x轴，纵轴为y轴，分别建立平面直角坐标系。记录上述特征点在RGB图像平面坐标系中的位置，通过矩阵对应关系在DEM数据中选择相同位置的点作为DEM 数据特征点，即根据DEM数据和RGB影像的矩阵对应关系，将上述正确的特征点由RGB影像坐标系映射至DEM数据的坐标系中，提取多期DEM数据间的同名特征点。
步骤S3.4：使用奇异值分解法求解两期DEM数据同名特征点间的初始转换关系。
步骤S3.4.1现有上述特征提取得到的两组点，第一组来组第一期影像，假设其为点集{q|(xi,yi,zi),i＝1,2...,n}，第二组来自第二期影像，并与点集q具有对应关系，设其为{p|(x'i,y'i,z'i),i＝1,2...,n}；
步骤S3.4.2求重心：求取点集p的重心p'以及点集q的重心q'；
步骤S3.43去重心化：计算每个点的去重心坐标，即pi＝pi-p'，qi＝qi-q'；
步骤S3.4.5计算H矩阵：
式中，n为特征点的个数；
步骤S3.4.6求解旋转矩阵：对H矩阵进行SVD分解，即H＝UΣVT，旋转矩阵 R0＝VUT；
步骤S3.4.7求解平移矩阵：平移矩阵T0＝p-R0p'。
步骤S3.5：粗匹配，利用上述求解出的旋转矩阵R0和平移矩阵T0，采用公式 P＝R×p+T对待匹配DEM数据进行坐标变换，此时多期DEM已具有较正确的位置关系。图5为本发明提供的具有初始匹配关系的DEM数据示意图。
上式中，p为待匹配DEM数据，P为待匹配DEM数据初始变换后的DEM数据， R0为旋转矩阵，T0为平移矩阵。
步骤S4：对上述已具有较好初始转换的DEM数据，通过插值方法建立多期DEM 数据的高程对应关系，利用统计方法检测出多期DEM数据中位于形变区域的点，将其剔除，即上述步骤104的详细步骤。
步骤S4.1设基准DEM为Q，其格网数为n，高程值为{Zqi|i＝1,2...,n}，待匹配 DEM为P，高程值为{Zpj|j＝1,2...,m},其格网数为m。
采用空间插值方法，建立多期DEM之间点的对应关系，即插值出基准DEM数据与待匹配DEM数据在相同平面坐标位置处的对应高程值，即插值出Q在P中相同平面坐标处的高程值Z，其数量为m。
将Zp与Z作差，获取高差关系△h。
假设上述高差关系△h由两个高斯模型混合而成，即系统误差模型N1与形变误差模型N2，根据实际情况设定N1与N2的初始值，包括N1的μ1、σ1、w1及N2的μ2、σ2、w2。通过E-step和M-step迭代计算出μ1、σ1、w1及μ2、σ2、w2的真值。
上式中，其中μ为高斯模型中的均值，σ为高斯模型中的方差，w为各个模型所占的权重，即μ1为系统误差模型N1的均值，μ2为形变误差模型N2的均值，σ1为系统误差模型N1的方差，σ2为形变误差模型N2的方差，w1为系统误差模型N1的权重， w2为形变误差模型N2的权重。
E-step：根据上一次迭代的模型参数(第一次迭代时为初始值)，计算出隐性变量的后验概率，作为隐藏变量的现估计值，即
M-step：根据上一步计算出的隐藏变量的估计值，将似然函数最大化以获得新的参数值，即最大似然估计，将似然函数最大化以获得新的参数值，即
重复上述E-step和M-step直至μ1、σ1、w1及μ2、σ2、w2趋于稳定，将最后一次计算出的值作为模型N1和模型N2的真值，即将上述步骤迭代计算至值收敛，并以最后一次的值作为模型的真值。
步骤S4.2根据正态分布中“3σ”准则，即模型上任一点落在区间[μ-3σ,μ+3σ]的概率是99.7％，计算μ1-3σ1与μ2-3σ2，将高差位于该区间内对应的P中的点，作为P的有效点，并剔除其余的点，即根据正态分布中“3σ”准则，即模型上任一点落在区间 [μ-3σ,μ+3σ]的概率是99.7％，据此提取出属于系统误差的点，并剔除掉属于形变模型的点。
通过上述可知，在一个实施例中，对获得优选的初始位置关系的DEM数据，通过插值方法建立多期DEM数据的高程对应关系，根据所述高程对应关系，利用统计方法检测出多期DEM数据中位于形变区域的点，将位于形变区域的点剔除，得到初步筛选后的多期DEM数据，可以包括：
设获得优选的初始位置关系的DEM数据中任一期DEM数据为基准DEM数据，利用空间插值方法，插值出基准DEM数据与其它期待匹配DEM数据在相同平面坐标位置处的对应高程值；
将多期DEM数据的对应高程值作差，获得高差关系；
假设所述高差关系由如下两个高斯模型混合而成：系统误差模型N1与形变误差模型N2；根据所述高差关系，按照如下步骤将位于形变区域的点剔除，得到初步筛选后的多期DEM数据：
E-step：根据上一次迭代的模型参数，计算出隐性变量的后验概率，作为隐藏变量的估计值；
M-step：根据计算出的隐藏变量的估计值，将似然函数最大化以获得新的参数值；
重复所述E-step和M-step迭代计算直至值收敛，将最后一次的值作为模型N1和模型N2的真值；
基于正态分布中3σ准则，根据模型N1和模型N2的真值，提取出属于系统误差模型的点，并剔除掉属于形变误差模型的点，得到初步筛选后的多期DEM数据，并假定初步筛选后的第二期DEM为Pem。
具体实施时，基准DEM数据可以是任一期DEM数据，但是通常以第一期的为基准DEM数据，因为其未发生形变，并且便于最后与形变数据进行差值分析。
步骤S5：对上述DEM数据，利用距离衰减效应，检测出多期DEM数据中未被上述统计方法检测出的位于形变区域的点，将其剔除，即上述步骤105的详细步骤。
在某些特殊情况下，形变量较小或监测区域中存在小部分较小形变的点，由于其较小的高差从而无法被上述统计方法探测出，但仍会对终结精匹配造成干扰。
根据地理学定律，距离形变区域越近的点，其本身是形变点的概率也越大，用 Gij＝KWiWjf(dij)来量化距离与空间交互的影响，根据上述步骤S4的计算结果，将形变误差模型N2中的μ2设定为阈值，即高差大于该阈值的点对，确定为形变点，将该点集记为Cp。
上式中，K为常量系数、Wi和Wj分别为空间实体i和j的规模，f(d)用以刻画距离因素的影响。
本发明针对实际数据对其进行约束和化简，空间实体规模用点数衡量，Wi和Wj分别代表点云P中的每一个点和点集Cp，因此Wi与Wj恒不变，令其为1，令常量系数K＝1，同时，为使探测结果更准确，添加高差因素，即对于Pem(初步筛选后的DEM数据)中的每个点，令hi＝|hi|/Hmax，即由式hi＝|hi|/Hmax计算每个点高程占比值，其中hi为第i个点对关系的高差，Hmax为所有点对中最大高差。最终，采用给定点集Cp对Pem中每个点的空间影响，令vi＝Gi/∑G为每个点的权值。具体做法如下：
对Pem中每个点，计算其到点集Cp中每个点的距离d；
由
计算距离衰减函数，每个点的衰减距离，其中距离d的个数为n 个(与Cp中的点数相同)，
为d的平均值，计算出f(d)的个数为n个；
计算hi＝|hi|/Hmax；
计算
其中
为f(d)的平均值；
用vi＝Gi/∑G计算每个点的权值，并归一化至[0,1]区间；
采用二分法剔除形变点，即：设定阈值m，若vi大于阈值m，则认定其对应的 Pem中的点为形变点，与予剔除，即由每个点的衰减距离和高程占比值计算每个点的空间衰减效应值，并作为形变点与非形变点判断依据，采用二分法分割形变点与非形变点。
通过上述可知，在一个实施例中，利用距离衰减效应，检测出初步筛选后的多期DEM数据中未被统计方法检测出的位于形变区域的点，将未被统计方法检测出的位于形变区域的点剔除，得到二次筛选后的多期DEM数据，可以包括：
计算初步筛选后的多期DEM数据中每个点到位于形变区域的每个点的衰减距离；
计算初步筛选后的多期DEM数据中每个点的高程占比值；
根据所述衰减距离和高程占比值，确定初步筛选后的多期DEM数据中每个点的空间衰减效应值；
根据所述空间衰减效应值，采用二分法剔除形变点，得到二次筛选后的多期DEM数据。
步骤S6：对上述剔除形变点的DEM数据，采用ICP算法进行整体精匹配，并将上述精匹配后的多期DEM数据构建差分模型，得到监测区域形变量，即上述步骤106 的详细步骤。
步骤S6.1：ICP算法方法，即采用ICP算法对上述基准DEM数据与剔除形变点后的待匹配DEM数据进行精确整体配准。
i.对Pse(经过二次筛选后的DEM数据)中的每一个点Psei，找到其在Q中最近的点Qi；
ii计算旋转矩阵R和平移矩阵T，使得误差函数
最小；
iii.对Pse使用上一步获取的平移和旋转参数进行刚体变换，得到新的点集P’；
iv.重复上述i-iv步直至误差小于预先设定的阈值，完成精匹配。
步骤S6.2：差异探测方法，即将上述配准后的多期DEM数据作差，获取监测区域形变量。
完成精匹配之后，在多期DEM之间构建差分模型，即△Z＝f(x,y)-g(x,y)。
上式中，△Z即为形变量，f(x,y)为形变前DEM在平面位置(x,y)处的高程值，g(x,y) 为形变后DEM在平面位置(x,y)处的高程值。图6为本发明提供的某区域形变探测结果图，根据该图以及发明人的大量实验表明，本发明实施例提供的无控制点的DEM 匹配与形变量探测方法提高了形变量探测的效率和精度，进而对地质灾害监测起到了重大的指导作用。
通过上述可知，在一个实施例中，对二次筛选后的多期DEM数据，采用ICP算法进行整体精匹配，得到精匹配后的多期DEM数据；根据精匹配后的多期DEM数据构建差分模型，得到监测区域形变量，可以包括：
采用ICP算法对基准DEM数据与二次筛选后的待匹配DEM数据按预设精确度进行整体配准处理，得到配准后的多期DEM数据；
将配准后的多期DEM数据作差，得到监测区域形变量。
综上，本发明实施例提供的无控制点的DEM匹配与应变量探测方法：基于无人机高分辨率航摄影像，首先基于DEM数据对应的RGB影像提取同名特征点，获取初始转换参数，为多期DEM数据提供较好的初始姿态；其次，将统计学方法与地理学方法引入粗差剔除中，有效剔除会对精匹配过程造成干扰的形变点；最后，采用ICP 算法进行精匹配并通过构建差分模型获取精确的监测区域形变量，可分析滑坡灾害的时空形变特征，为灾害早期预报预警提供决策支撑，可在难以布设控制点或未布设控制点的区域进行差异探测，对于地质灾害监测具有重要意义。
本发明实施例中还提供了一种无控制点的DEM匹配与形变量探测装置，如下面的实施例所述。由于该装置解决问题的原理与无控制点的DEM匹配与形变量探测方法相似，因此该装置的实施可以参见无控制点的DEM匹配与形变量探测方法的实施，重复之处不再赘述。
图8为本发明实施例中无控制点的DEM匹配与形变量探测装置的结构示意图，如图8所示，该装置包括：
获取单元01，用于获取监测区域的多期无人机航摄影像；
DEM数据及RGB图像确定单元02，用于根据多期无人机航摄影像，得到监测区域多期DEM数据及对应的RGB图像；
初始位置关系确定单元03，用于通过特征识别方法提取出多期RGB图像间的同名特征点，在DEM数据和RGB图像间建立模型点与像素的对应关系，根据对应关系将同名特征点映射至多期DEM数据中，确定多期DEM数据中的同名特征点间的转换关系，根据转换关系对DEM数据进行坐标变换，获得监测区域多期DEM数据间优选的初始位置关系；
初步筛选单元04，用于对获得优选的初始位置关系的DEM数据，通过插值方法建立多期DEM数据的高程对应关系，根据所述高程对应关系，利用统计方法检测出多期DEM数据中位于形变区域的点，将位于形变区域的点剔除，得到初步筛选后的多期DEM数据；
二次筛选单元05，用于利用距离衰减效应，检测出初步筛选后的多期DEM数据中未被统计方法检测出的位于形变区域的点，将未被统计方法检测出的位于形变区域的点剔除，得到二次筛选后的多期DEM数据；
形变量确定单元06，用于对二次筛选后的多期DEM数据，采用ICP算法进行整体精匹配，得到精匹配后的多期DEM数据；根据精匹配后的多期DEM数据构建差分模型，得到监测区域形变量；所述形变量用于分析监测区域的滑坡灾害。
在一个实施例中，DEM数据及RGB图像确定单元具体可以用于：
将监测区域的多期无人机航摄影像进行空三计算、密集匹配、点云滤波及规则格网构建处理，生成监测区域多期DEM数据；
将监测区域多期DEM数据按原分辨率导出对应的RGB图像，使DEM格网数与 RGB图像的像素数相对应。
在一个实施例中，初始位置关系确定单元具体可以用于：
使用SURF算法对多期RGB图像提取同名特征点；
使用RANSAC算法对所述同名特征点进行一致性检测，剔除错误配对的同名特征点，得到剔除处理后的RGB同名特征点；
以DEM数据左上角为坐标原点建立DEM平面直角坐标系，以RGB图像的左上角为坐标原点建立RGB平面直角坐标系，记录剔除处理后的RGB同名特征点在RGB 图像平面坐标系中的位置，根据DEM数据和RGB图像间的矩阵对应关系在DEM平面直角坐标系中选择与剔除处理后的RGB同名特征点相同位置的点作为DEM同名特征点；
使用奇异值分解法求解多期DEM同名特征点间的初始转换关系；
根据初始转换关系对DEM数据进行坐标变换，获得监测区域多期DEM数据间优选的初始位置关系。
在一个实施例中，使用RANSAC算法对所述同名特征点进行一致性检测，剔除错误配对的同名特征点，得到剔除处理后的RGB同名特征点，可以包括：
从第一期RGB图像的特征点集中随机选取一组内点，计算模型参数M；
用计算出的模型参数M去测试所有数据点，若各数据点与模型的误差小于设定的阈值，则将其视为内点，否则视为外点；所述内点为正确配对的特征点，外点为错误配对的特征点；
如果有大于预设数目的数据点被归类为内点，确定对应模型为合理模型，在多个合理模型中选出内点数目最多的模型作为最佳模型；
重复上述步骤直至达到设定的迭代次数，将达到设定的迭代次数时选出的最佳模型作为最终模型，选用最终模型对应的内点集作为剔除处理后的RGB同名特征点。
在一个实施例中，初步筛选单元具体可以用于：
设获得优选的初始位置关系的DEM数据中任一期DEM数据为基准DEM数据，利用空间插值方法，插值出基准DEM数据与其它期待匹配DEM数据在相同平面坐标位置处的对应高程值；
将多期DEM数据的对应高程值作差，获得高差关系；
假设所述高差关系由如下两个高斯模型混合而成：系统误差模型N1与形变误差模型N2；根据所述高差关系，按照如下步骤将位于形变区域的点剔除，得到初步筛选后的多期DEM数据：
E-step：根据上一次迭代的模型参数，计算出隐性变量的后验概率，作为隐藏变量的估计值；
M-step：根据计算出的隐藏变量的估计值，将似然函数最大化以获得新的参数值；
重复所述E-step和M-step迭代计算直至值收敛，将最后一次的值作为模型N1和模型N2的真值；
基于正态分布中3σ准则，根据模型N1和模型N2的真值，提取出属于系统误差模型的点，并剔除掉属于形变误差模型的点，得到初步筛选后的多期DEM数据。
在一个实施例中，二次筛选单元具体可以用于：
计算初步筛选后的多期DEM数据中每个点到位于形变区域的每个点的衰减距离；
计算初步筛选后的多期DEM数据中每个点的高程占比值；
根据所述衰减距离和高程占比值，确定初步筛选后的多期DEM数据中每个点的空间衰减效应值；
根据所述空间衰减效应值，采用二分法剔除形变点，得到二次筛选后的多期DEM数据。
在一个实施例中，形变量确定单元具体可以用于：
采用ICP算法对基准DEM数据与二次筛选后的待匹配DEM数据按预设精确度进行整体配准处理，得到配准后的多期DEM数据；
将配准后的多期DEM数据作差，得到监测区域形变量。
本发明实施例还提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述无控制点的DEM匹配与形变量探测方法。
本发明实施例还提供一种计算机可读存储介质，所述计算机可读存储介质存储有执行上述无控制点的DEM匹配与形变量探测方法的计算机程序。
本发明实施例的有益效果是：本发明实施例提供的无控制点DEM匹配与形变量探测方案通过：采用无人机航摄影像，基于鲁棒的特征提取方法提取出多期DEM数据间的同名特征点，以获得优选的初始位置姿态；在获取优选的初始位置姿态后，通过空间插值方法建立对应点关系，利用统计学方法监测出最大概率位于形变区域的点，将其剔除；同时，充分考虑DEM数据间的空间效应，利用距离衰减效应，探测出可能位于形变区域的微小形变点，将其剔除；在剔除会对精匹配过程造成干扰的形变点后，进行精确配准并获取精确形变量，可分析滑坡灾害的时空形变特征，为灾害早期预报预警提供决策支撑，可在难以定向布设控制点或未布设控制点的区域进行形变量精确探测，对于地质灾害监测具有重要意义。
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等) 上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (10)
1.一种无控制点的DEM匹配与形变量探测方法，其特征在于，包括：
获取监测区域的多期无人机航摄影像；
根据多期无人机航摄影像，得到监测区域多期DEM数据及对应的RGB图像；
通过特征识别方法提取出多期RGB图像间的同名特征点，在DEM数据和RGB图像间建立模型点与像素的对应关系，根据对应关系将同名特征点映射至多期DEM数据中，确定多期DEM数据中的同名特征点间的转换关系，根据转换关系对DEM数据进行坐标变换，获得监测区域多期DEM数据间优选的初始位置关系；
对获得优选的初始位置关系的DEM数据，通过插值方法建立多期DEM数据的高程对应关系，根据所述高程对应关系，利用统计方法检测出多期DEM数据中位于形变区域的点，将位于形变区域的点剔除，得到初步筛选后的多期DEM数据；
利用距离衰减效应，检测出初步筛选后的多期DEM数据中未被统计方法检测出的位于形变区域的点，将未被统计方法检测出的位于形变区域的点剔除，得到二次筛选后的多期DEM数据；
对二次筛选后的多期DEM数据，采用ICP算法进行整体精匹配，得到精匹配后的多期DEM数据；根据精匹配后的多期DEM数据构建差分模型，得到监测区域形变量；所述形变量用于分析监测区域的滑坡灾害。
2.如权利要求1所述的无控制点的DEM匹配与形变量探测方法，其特征在于，根据多期无人机航摄影像，得到监测区域多期DEM数据及对应的RGB图像，包括：
将监测区域的多期无人机航摄影像进行空三计算、密集匹配、点云滤波及规则格网构建处理，生成监测区域多期DEM数据；
将监测区域多期DEM数据按原分辨率导出对应的RGB图像，使DEM格网数与RGB图像的像素数相对应。
3.如权利要求1所述的无控制点的DEM匹配与形变量探测方法，其特征在于，通过特征识别方法提取出多期RGB图像间的同名特征点，在DEM数据和RGB图像间建立模型点与像素的对应关系，根据对应关系将同名特征点映射至多期DEM数据中，确定多期DEM数据中的同名特征点间的转换关系，根据转换关系对DEM数据进行坐标变换，获得监测区域多期DEM数据间优选的初始位置关系，包括：
使用SURF算法对多期RGB图像提取同名特征点；
使用RANSAC算法对所述同名特征点进行一致性检测，剔除错误配对的同名特征点，得到剔除处理后的RGB同名特征点；
以DEM数据左上角为坐标原点建立DEM平面直角坐标系，以RGB图像的左上角为坐标原点建立RGB平面直角坐标系，记录剔除处理后的RGB同名特征点在RGB图像平面坐标系中的位置，根据DEM数据和RGB图像间的矩阵对应关系在DEM平面直角坐标系中选择与剔除处理后的RGB同名特征点相同位置的点作为DEM同名特征点；
使用奇异值分解法求解多期DEM同名特征点间的初始转换关系；
根据初始转换关系对DEM数据进行坐标变换，获得监测区域多期DEM数据间优选的初始位置关系。
4.如权利要求3所述的无控制点的DEM匹配与形变量探测方法，其特征在于，使用RANSAC算法对所述同名特征点进行一致性检测，剔除错误配对的同名特征点，得到剔除处理后的RGB同名特征点，包括：
从第一期RGB图像的特征点集中随机选取一组内点，计算模型参数M；
用计算出的模型参数M去测试所有数据点，若各数据点与模型的误差小于设定的阈值，则将其视为内点，否则视为外点；所述内点为正确配对的特征点，外点为错误配对的特征点；
如果有大于预设数目的数据点被归类为内点，确定对应模型为合理模型，在多个合理模型中选出内点数目最多的模型作为最佳模型；
重复上述步骤直至达到设定的迭代次数，将达到设定的迭代次数时选出的最佳模型作为最终模型，选用最终模型对应的内点集作为剔除处理后的RGB同名特征点。
5.如权利要求1所述的无控制点的DEM匹配与形变量探测方法，其特征在于，对获得优选的初始位置关系的DEM数据，通过插值方法建立多期DEM数据的高程对应关系，根据所述高程对应关系，利用统计方法检测出多期DEM数据中位于形变区域的点，将位于形变区域的点剔除，得到初步筛选后的多期DEM数据，包括：
设获得优选的初始位置关系的DEM数据中任一期DEM数据为基准DEM数据，利用空间插值方法，插值出基准DEM数据与其它期待匹配DEM数据在相同平面坐标位置处的对应高程值；
将多期DEM数据的对应高程值作差，获得高差关系；
假设所述高差关系由如下两个高斯模型混合而成：系统误差模型N1与形变误差模型N2；根据所述高差关系，按照如下步骤将位于形变区域的点剔除，得到初步筛选后的多期DEM数据：
E-step：根据上一次迭代的模型参数，计算出隐性变量的后验概率，作为隐藏变量的估计值；
M-step：根据计算出的隐藏变量的估计值，将似然函数最大化以获得新的参数值；
重复所述E-step和M-step迭代计算直至值收敛，将最后一次的值作为模型N1和模型N2的真值；
基于正态分布中3σ准则，根据模型N1和模型N2的真值，提取出属于系统误差模型的点，并剔除掉属于形变误差模型的点，得到初步筛选后的多期DEM数据。
6.如权利要求1所述的无控制点的DEM匹配与形变量探测方法，其特征在于，利用距离衰减效应，检测出初步筛选后的多期DEM数据中未被统计方法检测出的位于形变区域的点，将未被统计方法检测出的位于形变区域的点剔除，得到二次筛选后的多期DEM数据，包括：
计算初步筛选后的多期DEM数据中每个点到位于形变区域的每个点的衰减距离；
计算初步筛选后的多期DEM数据中每个点的高程占比值；
根据所述衰减距离和高程占比值，确定初步筛选后的多期DEM数据中每个点的空间衰减效应值；
根据所述空间衰减效应值，采用二分法剔除形变点，得到二次筛选后的多期DEM数据。
7.如权利要求1所述的无控制点的DEM匹配与形变量探测方法，其特征在于，对二次筛选后的多期DEM数据，采用ICP算法进行整体精匹配，得到精匹配后的多期DEM数据；根据精匹配后的多期DEM数据构建差分模型，得到监测区域形变量，包括：
采用ICP算法对基准DEM数据与二次筛选后的待匹配DEM数据按预设精确度进行整体配准处理，得到配准后的多期DEM数据；
将配准后的多期DEM数据作差，得到监测区域形变量。
8.一种无控制点的DEM匹配与形变量探测装置，其特征在于，包括：
获取单元，用于获取监测区域的多期无人机航摄影像；
DEM数据及RGB图像确定单元，用于根据多期无人机航摄影像，得到监测区域多期DEM数据及对应的RGB图像；
初始位置关系确定单元，用于通过特征识别方法提取出多期RGB图像间的同名特征点，在DEM数据和RGB图像间建立模型点与像素的对应关系，根据对应关系将同名特征点映射至多期DEM数据中，确定多期DEM数据中的同名特征点间的转换关系，根据转换关系对DEM数据进行坐标变换，获得监测区域多期DEM数据间优选的初始位置关系；
初步筛选单元，用于对获得优选的初始位置关系的DEM数据，通过插值方法建立多期DEM数据的高程对应关系，根据所述高程对应关系，利用统计方法检测出多期DEM数据中位于形变区域的点，将位于形变区域的点剔除，得到初步筛选后的多期DEM数据；
二次筛选单元，用于利用距离衰减效应，检测出初步筛选后的多期DEM数据中未被统计方法检测出的位于形变区域的点，将未被统计方法检测出的位于形变区域的点剔除，得到二次筛选后的多期DEM数据；
形变量确定单元，用于对二次筛选后的多期DEM数据，采用ICP算法进行整体精匹配，得到精匹配后的多期DEM数据；根据精匹配后的多期DEM数据构建差分模型，得到监测区域形变量；所述形变量用于分析监测区域的滑坡灾害。
9.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至7任一所述方法。
10.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有执行权利要求1至7任一所述方法的计算机程序。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110428378.1A CN113159174B (zh) | 2021-04-21 | 2021-04-21 | 无控制点的dem匹配与形变量探测方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110428378.1A CN113159174B (zh) | 2021-04-21 | 2021-04-21 | 无控制点的dem匹配与形变量探测方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113159174A CN113159174A (zh) | 2021-07-23 |
| CN113159174B true CN113159174B (zh) | 2022-05-06 |
Family
ID=76867530
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110428378.1A Active CN113159174B (zh) | 2021-04-21 | 2021-04-21 | 无控制点的dem匹配与形变量探测方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113159174B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115908954B (zh) * | 2023-03-01 | 2023-07-28 | 四川省公路规划勘察设计研究院有限公司 | 基于人工智能的地质灾害隐患识别系统、方法及电子设备 |
Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105893972A (zh) * | 2016-04-08 | 2016-08-24 | 深圳市智绘科技有限公司 | 基于影像进行的违章建筑物自动监测方法及其实现系统 |
| CN107132539A (zh) * | 2017-05-03 | 2017-09-05 | 中国地质科学院探矿工艺研究所 | 一种基于小基线集的时间序列InSAR的滑坡早期识别方法 |
| CN108595373A (zh) * | 2018-01-31 | 2018-09-28 | 中南林业科技大学 | 一种无控制dem配准方法 |
| CN110084817A (zh) * | 2019-03-21 | 2019-08-02 | 西安电子科技大学 | 基于深度学习的数字高程模型生产方法 |
| CN111207728A (zh) * | 2020-03-03 | 2020-05-29 | 王刚强 | 一种无人机测量边坡及基坑开挖量的设备 |
| CN111340950A (zh) * | 2020-02-18 | 2020-06-26 | 云南师范大学 | 数字高程模型重建方法及装置 |
| CN111683221A (zh) * | 2020-05-21 | 2020-09-18 | 武汉大学 | 嵌入矢量红线数据的自然资源实时视频监测方法及系统 |
| CN112434709A (zh) * | 2020-11-20 | 2021-03-02 | 西安视野慧图智能科技有限公司 | 基于无人机实时稠密三维点云和dsm的航测方法及系统 |
| CN112561981A (zh) * | 2020-12-16 | 2021-03-26 | 王静 | 一种融合影像信息的摄影测量点云滤波方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9488492B2 (en) * | 2014-03-18 | 2016-11-08 | Sri International | Real-time system for multi-modal 3D geospatial mapping, object recognition, scene annotation and analytics |
-
2021
- 2021-04-21 CN CN202110428378.1A patent/CN113159174B/zh active Active
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105893972A (zh) * | 2016-04-08 | 2016-08-24 | 深圳市智绘科技有限公司 | 基于影像进行的违章建筑物自动监测方法及其实现系统 |
| CN107132539A (zh) * | 2017-05-03 | 2017-09-05 | 中国地质科学院探矿工艺研究所 | 一种基于小基线集的时间序列InSAR的滑坡早期识别方法 |
| CN108595373A (zh) * | 2018-01-31 | 2018-09-28 | 中南林业科技大学 | 一种无控制dem配准方法 |
| CN110084817A (zh) * | 2019-03-21 | 2019-08-02 | 西安电子科技大学 | 基于深度学习的数字高程模型生产方法 |
| CN111340950A (zh) * | 2020-02-18 | 2020-06-26 | 云南师范大学 | 数字高程模型重建方法及装置 |
| CN111207728A (zh) * | 2020-03-03 | 2020-05-29 | 王刚强 | 一种无人机测量边坡及基坑开挖量的设备 |
| CN111683221A (zh) * | 2020-05-21 | 2020-09-18 | 武汉大学 | 嵌入矢量红线数据的自然资源实时视频监测方法及系统 |
| CN112434709A (zh) * | 2020-11-20 | 2021-03-02 | 西安视野慧图智能科技有限公司 | 基于无人机实时稠密三维点云和dsm的航测方法及系统 |
| CN112561981A (zh) * | 2020-12-16 | 2021-03-26 | 王静 | 一种融合影像信息的摄影测量点云滤波方法 |
Non-Patent Citations (3)
| Title |
|---|
| Point Cloud Registration Based on One-Point RANSAC and Scale-Annealing Biweight Estimation;Jiayuan Li等;《IEEE》;20210108;第9716-9729页 * |
| 无人机RGB影像中人工林单木位置的提取;李晨 等;《林业工程学报》;20210106;第141-147页 * |
| 高程异常模型的已知点框架约束;徐新强 等;《海洋测绘》;20061130;第59-61页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113159174A (zh) | 2021-07-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107038717B (zh) | 一种基于立体栅格自动分析3d点云配准误差的方法 | |
| CN108304883B (zh) | 基于改进sift的sar图像匹配方法 | |
| CN111524168B (zh) | 点云数据的配准方法、系统、装置及计算机存储介质 | |
| Theiler et al. | Markerless point cloud registration with keypoint-based 4-points congruent sets | |
| CN106651942A (zh) | 基于特征点的三维旋转运动检测与旋转轴定位方法 | |
| WO2008155259A1 (en) | Method of segmenting anatomic entities in 3d digital medical images | |
| CN106558072A (zh) | 一种基于改进sift特征在遥感图像上配准的方法 | |
| CN115601574A (zh) | 一种改进akaze特征的无人机影像匹配方法 | |
| CN114663373A (zh) | 一种用于零件表面质量检测的点云配准方法及装置 | |
| Sveier et al. | Object detection in point clouds using conformal geometric algebra | |
| CN113159174B (zh) | 无控制点的dem匹配与形变量探测方法及装置 | |
| CN110793441B (zh) | 一种高精度物体几何尺寸测量方法及装置 | |
| CN113450269A (zh) | 一种基于3d视觉的点云关键点提取方法 | |
| CN109241981B (zh) | 一种基于稀疏编码的特征检测方法 | |
| WO2008155255A1 (en) | Method of constructing a gray value model and/or a geometric model of an anatomic entity in a 3d digital medical image | |
| CN114663880A (zh) | 基于多层级跨模态自注意力机制的三维目标检测方法 | |
| CN114463503A (zh) | 三维模型和地理信息系统的融合方法及装置 | |
| CN110223356A (zh) | 一种基于能量生长的单目相机全自动标定方法 | |
| CN115713548A (zh) | 一种多期实景三维模型自动配准方法 | |
| CN116645500A (zh) | 边界框的确定方法、装置、设备、存储介质和程序产品 | |
| Morell et al. | 3d maps representation using gng | |
| Wan et al. | A performance comparison of feature detectors for planetary rover mapping and localization | |
| CN116091562A (zh) | 基于二维投影线段的建筑物点云自动配准的方法 | |
| Dierenbach et al. | Next-Best-View method based on consecutive evaluation of topological relations | |
| WO2020197495A1 (en) | Method and system for feature matching |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |