CN116125490B - 面向滑坡体形变场时序监测的tls多目标优化选址方法 - Google Patents
面向滑坡体形变场时序监测的tls多目标优化选址方法 Download PDFInfo
- Publication number
- CN116125490B CN116125490B CN202310053117.5A CN202310053117A CN116125490B CN 116125490 B CN116125490 B CN 116125490B CN 202310053117 A CN202310053117 A CN 202310053117A CN 116125490 B CN116125490 B CN 116125490B
- Authority
- CN
- China
- Prior art keywords
- landslide
- individual
- voxels
- candidate
- site
- 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
- 238000012544 monitoring process Methods 0.000 title claims abstract description 26
- 238000005457 optimization Methods 0.000 title claims abstract description 21
- 238000010187 selection method Methods 0.000 title claims abstract description 10
- 238000000034 method Methods 0.000 claims abstract description 34
- 238000004088 simulation Methods 0.000 claims abstract description 16
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 15
- 238000012216 screening Methods 0.000 claims abstract description 14
- 239000013598 vector Substances 0.000 claims description 33
- 238000006073 displacement reaction Methods 0.000 claims description 17
- 101100397224 Bacillus subtilis (strain 168) isp gene Proteins 0.000 claims description 15
- 101100052502 Shigella flexneri yciB gene Proteins 0.000 claims description 15
- 101150064873 ispA gene Proteins 0.000 claims description 15
- 238000004364 calculation method Methods 0.000 claims description 14
- 230000002776 aggregation Effects 0.000 claims description 3
- 238000004220 aggregation Methods 0.000 claims description 3
- 238000013507 mapping Methods 0.000 claims description 3
- 101100533306 Mus musculus Setx gene Proteins 0.000 claims 1
- 230000008569 process Effects 0.000 description 11
- 238000011160 research Methods 0.000 description 7
- 238000012360 testing method Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 2
- 238000002310 reflectometry Methods 0.000 description 2
- 239000000758 substrate Substances 0.000 description 2
- 238000007792 addition Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000011835 investigation 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
- 230000002265 prevention Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
- G01S17/88—Lidar systems specially adapted for specific applications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B11/00—Measuring arrangements characterised by the use of optical techniques
- G01B11/16—Measuring arrangements characterised by the use of optical techniques for measuring the deformation in a solid, e.g. optical strain gauge
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T10/00—Road transport of goods or passengers
- Y02T10/10—Internal combustion engine [ICE] based vehicles
- Y02T10/40—Engine management systems
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Computer Networks & Wireless Communication (AREA)
- Electromagnetism (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Geophysics And Detection Of Objects (AREA)
- Testing Or Calibration Of Command Recording Devices (AREA)
Abstract
本发明公开了面向滑坡体形变场时序监测的TLS多目标优化选址方法,生成观测区域的三维点云,将三维点云转换为体素,建立观测区域的初始三维模型;在初始三维模型中确定候选站点;确定TLS扫描仪的参数,基于扫描仪参数以及初始三维模型为约束进行扫描模拟,筛选视线不被遮挡的滑坡体素为已观测体素,筛选视线被遮挡的滑坡体素为未观测体素;计算观测区域的站点基数;采用约束NSGA‑II算法生成Pareto最优解集;从Pareto最优集中决策出最终个体。本发明更全面地考虑了影响滑坡隐患区域TLS站点采集质量的因素,并进行多目标优化,进而提高了TLS站点的可靠性。
Description
技术领域
本发明属于滑坡隐患区域时序监测中的地面激光雷达点云采集领域,具体涉及面向滑坡体形变场时序监测的TLS多目标优化选址方法。
背景技术
滑坡监测预警,是当前地学领域研究的热点和难点。目前,主要通过探测滑坡区域地表厘米级微小形变作为滑坡即将发生的预警信息。地面三维激光扫描(TLS),是一种地表高密度、高精度三维坐标快速获取手段,可对滑坡隐患区域进行不同时相观测,并经过一系列后处理,以此获取整个滑坡隐患区域高分辨形变场。在滑坡监测预警过程中,缩小激光点云观测误差使TLS探测到滑坡区域毫米级形变至关重要。为此,大量学者从仪器误差、气候条件等方面提高激光点云观测精度。然而,观测几何对点云精度也产生了影响,却鲜有关注和研究,特别是大范围的自然场景区域无相关研究。目前,滑坡区域TLS站点构建主要依赖于用户的主观经验,几乎不考虑站点对精度的影响,严重影响成本和数据质量。尤其是在需要设置固定观测墩的大范围滑坡隐患区域长时序监测中,此类问题将变得愈加严重。因此,如何科学布设滑坡区域TLS观测站点,是一个重要问题。
TLS站点规划,主要采用模型模拟方法,在现有的观测场景模型基础上,模拟候选测站TLS扫描,建立观测完整、成本等约束方程,并采用贪婪等优化方法确定TLS站点位置。目前,国内外相关研究较少,现有研究主要以基于二维地形图的小范围区域模拟规划为主。但是,滑坡隐患区域通常具有较大的高差,这导致使用基于二维模型的方法必然导致模拟失真,进而影响规划的TLS站点的可靠性。因此,少数研究人员开始对TLS站点三维模拟规划方法进行了探索,并在建筑物、雕塑等观测方面进行了尝试。然而,在大范围滑坡隐患区域TLS站点三维模拟规划方面,却还未有相关研究,其中面临如下问题和挑战:1)现有方法无法从观测几何的角度制定出精度最高的滑坡隐患区域TLS站点,需考虑如何量化TLS观测精度,并使其作为站点优化的目标;2)现有方法仅考虑单目标,影响滑坡隐患区域最佳TLS站点的条件是多元的、复杂的,需要引入多目标优化算法;3)现有方法规划结果在采集阶段获取的多站数据不一定能在处理阶段进行自动配准,需考虑如何保证优化后TLS站点中站点之间的可配准性。
发明内容
本发明的目的在于针对现有技术中缺乏对滑坡隐患区域TLS站点规划方法,导致的观测精度不足,进而无法探测微小形变的问题,提供面向滑坡体形变场时序监测的TLS多目标优化选址方法。利用本发明可以设计出滑坡隐患区域的高精度TLS站点,并可以保障站与站之间的可配准性,从而优化激光点云观测几何,降低激光点云的点位精度,以此服务于大范围滑坡隐患区域长时序监测,提高TLS形变监测能力,为灾害机理研究及防灾减灾提供重要的技术支撑。
为了解决上述技术问题,本发明采用如下的技术方案:
面向滑坡体形变场时序监测的TLS多目标优化选址方法,包括以下步骤:
步骤1:生成观测区域的三维点云,将三维点云转换为体素,根据体素建立观测区域的初始三维模型;
步骤2:在初始三维模型中确定候选站点,构成候选站点集合J,候选站点总数为N;
步骤3:确定TLS扫描仪的参数,基于扫描仪参数以及初始三维模型为约束进行扫描模拟,筛选视线不被遮挡的滑坡体素为已观测体素,筛选视线被遮挡的滑坡体素为未观测体素;
步骤4:计算观测区域的站点基数B sp ;
步骤5:随机生成长度为候选站点总数N的二进制编码,每一条二进制编码作为一个个体,共生成数量为n的个体组成为初始父代P t ,采用约束NSGA-II算法生成Pareto最优解集;
步骤6:基于PEG-MCDM算法从Pareto最优集中决策出最终个体。
如上所述步骤1包括以下步骤:
步骤1.1:用无人机获取观测区域垂直影像和倾斜影像;
步骤1.2:将无人机采集的垂直影像和倾斜影像重建为观测区域的三维点云;
步骤1.3:对观测区域的三维点云进行点云分类,将位于滑坡隐患区域内的点云标记为滑坡点,将滑坡隐患区域外的点云依照其属性划分为地面点与其他点;
步骤1.4:将空间分割为具有定义大小的三维体素,并将点云与点云类型分配给相应的体素,如果一个体素包含点云,则将保留该体素;否则,该体素将被删除,然后将保留的体素作为观测区域的初始三维模型,其中滑坡点云转换为滑坡体素的数量为N v 。
如上所述步骤2包括以下步骤:
步骤2.1:在初始三维模型构建网格网络,并将网格的中心点作为拟候选站点位置;
步骤2.2:遍历所有的拟候选站点位置,如果与该拟候选站点位置距离最近体素的点云类型为其他点,则将该拟候选站点删除;
步骤2.3:计算剩余拟候选站点所处位置的坡度,并将所有坡度值大于45的拟候选站点位置删除;
步骤2.4:剩余的所有拟候选站点将被确定为候选站点,组成候选站点集合J,候选站点总数为N。
如上所述步骤3包括以下步骤:
步骤3.1:确定TLS扫描仪的参数,依照遍历顺序选出第一个候选站点,确定候选站点的扫描范围;
步骤3.2:遍历所有滑坡体素,计算候选站点到滑坡体素的方向向量与观测目标点所在曲面法向量的夹角关系,筛选满足预设夹角条件的滑坡体素;
步骤3.3:计算候选站点与滑坡体素的水平距离,筛选满足预设水平距离条件的滑坡体素;
步骤3.4:筛选视线不被遮挡的滑坡体素为已观测体素;
步骤3.5:遍历所有候选站点,已观测体素为满足步骤3.1 - 3.4的滑坡体素,不满足的为未观测体素。
如上所述步骤4包括以下步骤:
将所有滑坡体素设置为未覆盖,每次都选出覆盖最多未覆盖滑坡体素的候选站点,并将所选候选站点覆盖的未覆盖滑坡体素标记为已覆盖滑坡体素,当所选候选站点需要标记的未覆盖滑坡体素数量低于所设定的阈值时,停止迭代,此时已选出的候选站点数量为站点基数B sp 。
如上所述步骤5包括以下步骤:
步骤5.1:设定种群大小n,随机生成长度为候选站点总数N的二进制编码,其中仅有1与0组成,1代表该位置序号的候选站点处设站,0代表该位置序号的候选站点不设站,同时二进制编码中1对应的位置序号,组成该个体的所选站点序号集A,记个体的所选站点序号集A中序号的个数为u,每一条二进制编码作为一个个体,共生成数量为n的个体组成为初始父代P t ;
步骤5.2:对父代种群P t 中的个体进行交叉和变异,生成一个种群大小同样为n的子代种群Q t ,合并父代种群P t 与子代种群Q t ,得到种群大小为2n的新种群R t ;
步骤5.3:计算新种群R t 中各个体的控制规划结果成本的适应度函数F 1和控制规划结果质量的适应度函数F 2;
步骤5.4:基于各个体的控制规划结果成本的适应度函数F 1和控制规划结果质量的适应度函数F 2对新种群R t 进行非支配排序,判断新种群R t 中所有个体之间的支配关系;
步骤5.5:通过拥挤度对新种群R t 中同一Pareto等级中的个体进行排序,同一Pareto等级中的个体拥挤度越高,排名越靠前;
步骤5.6:从排序后的新种群R t 中选择拥挤度靠前的n个个体组成下一代父代种群P t+1 ;
步骤5.7:重复步骤5.2 - 5.6直到迭代次数大于设定阈值,输出最后一代种群的Pareto等级最高的一组个体作为Pareto最优集。
如上所述步骤5.3包括以下步骤:
控制规划结果成本的适应度函数F 1基于以下公式:
其中,j为滑坡体素的序号,N v 为滑坡体素的总数,a为候选站点的序号,N为候选站点的总数,lv j 为第j个滑坡体素的标签,已测量滑坡体素表示该滑坡体素在个体中的多个候选站点中有被作为已观测体素,未测量滑坡体素表示该滑坡体素在个体的多个候选站点中没有被作为已观测体素,p 0是惩罚因子,cp a 是第a个候选站点的标签,如果cp a 为1表示已选择该候选站点,而cp a 为0则表示未选择该候选站点,B sp 为站点基数;
控制规划结果质量的适应度函数F 2基于以下公式:
其中,k为已测量滑坡体素的索引;A为个体所选候选站点序号集;为第a个候选站点与第k个已测量滑坡体素之间的入射角;为第a个候选站点与第k个已测量滑坡体素之间的距离;r max为所使用的扫描仪最大射程;r min为所使用的扫描仪最大射程。
如上所述步骤5.4包括以下步骤:
假定个体p与个体q,
若个体p的F 1和F 2均分别大于个体q的F 1和F 2时,则视为p支配q,
若个体p的F 1等于个体q的F 1,个体p的F 2大于个体q的F 2,则视为p支配q,
若个体p的F 2等于个体q的F 2,个体p的F 1大于个体q的F 1,则视为p支配q,
若个体p的F 1大于个体q的F 1,而个体p的F 2小于个体q的F 2,则视为p与q为非支配关系,
若个体p的F 2大于个体q的F 2,而个体p的F 1小于个体q的F 1,则视为p与q为非支配关系,
依照个体被支配的次数,划分Pareto等级,被支配的次数越少,个体越优,等级越高,个体不满足约束条件,则直接视为Pareto最低等级。
如上所述步骤5.4中的约束条件为:
重叠度O m 计算包括以下步骤:
若个体的候选站点数为1,则不考虑O lim约束;
若个体的候选站点数为2,则直接计算两候选站点的重叠度,相邻两候选站点同时覆盖的滑坡体素数量与两候选站点中覆盖滑坡体素数量较低的候选站点所覆盖的滑坡体素数量的比值为两站重叠度;
若个体的候选站点数大于等于3,则构建候选站点的Delaunay三角网,计算Delaunay三角网的各边的重叠度,并以重叠度倒数为权重构建最小生成树,该最小生成树为个体的配准路径,其中路径的序号为m。
如上所述步骤6包括以下步骤:
步骤6.2:将Pareto最优集中个体的归一化适应度函数值从最小值到最大值顺序重新排序形成主准则向量x 1,将Pareto最优集中个体的归一化适应度函数值从最小值到最大值顺序重新排序形成主准则向量x 2,将主准则向量x 1与主准则向量x 2反向排列,并相互交换获得主聚合向量y 1与主聚合向量y 2;
步骤6.5:根据PEG理论估计两个适应度函数的最优值F 1 0与F 2 0;
与现有技术相比,本发明具有以下有益效果:
多目标:本发明提供了一种建立新的多目标TLS选址方法,更全面地考虑了影响滑坡隐患区域TLS站点采集质量的因素,并进行多目标优化,进而提高了TLS站点的可靠性;
高精度:本发明通过量化观测几何的精度指标,并使其作为优化目标,有效提高了滑坡隐患区域TLS站点采集数据的精度,进而使本发明具备更高的TLS形变监测能力;
可配准性:本发明在约束NSGA-II优化过程中进行非支配排序时,融入了对重叠度约束。通过对站点的重叠度进行了检验,降低了存在低于20%相邻站点的方案的Pareto等级,使得Pareto最优解集中所有方案均满足配准所需的重叠度需求,解决了传统方法所得TLS站点无法确保自动配准的问题,具有更好的普适性和实用性。
附图说明
图1为本发明的流程示意图。
图2武汉市江夏区滑坡测试区1模拟规划结果图,其中 (a) Pareto最优前沿图;(b)F 1的收敛性图;(c)F 2的收敛性图;(d) 最终TLS站点位置图。
图3武汉市江夏区滑坡测试区2模拟规划结果图 ,其中(a) Pareto最优前沿图;(b)F 1的收敛性图;(c)F 2的收敛性图;(d) 最终TLS站点位置图。
具体实施方式
为了便于本领域普通技术人员理解和实施本发明,下面结合实施例对本发明作进一步的详细描述,应当理解,此处所描述的实施示例仅用于说明和解释本发明,并不用于限定本发明。
实施例1:
面向滑坡体形变场时序监测的TLS多目标优化选址方法,具体包括以下步骤:
步骤1:生成观测区域的三维点云,将三维点云转换为体素,根据体素建立观测区域的初始三维模型。
具体的,生成观测区域的初始三维模型的过程如下:
步骤1.1:用无人机快速获取观测区域垂直影像和倾斜影像。
步骤1.2:使用Bentley公司旗下ContextCapture软件,将无人机采集的垂直影像和倾斜影像重建为观测区域的三维点云。
步骤1.3:对观测区域的三维点云进行点云分类。将位于滑坡隐患区域内的点云标记为滑坡点,将滑坡隐患区域外的点云依照其属性划分为地面点与其他点。
步骤1.4:将观测区域带有语义信息的三维点云转换为体素。在体素化的过程中,将空间分割为具有定义大小(如1 m)的三维体素,并将点云与点云类型分配给相应的体素。如果一个体素包含点云,则将保留该体素;否则,该体素将被删除。然后将保留的体素作为观测区域的初始三维模型,其中滑坡点云转换为的滑坡体素数量为N v 。
步骤2:在初始三维模型中确定候选站点,构成候选站点集合J,候选站点总数为N;。
具体的,基于设站可行性原则确定候选站点的过程如下:
步骤2.1:在初始三维模型构建一个均匀网格网络,并将网格的中心点作为拟候选站点位置。
步骤2.2:遍历所有的拟候选站点位置,检查与该拟候选站点位置距离最近体素的点云类型。如果点云类型为其他点,则将该拟候选站点删除。
步骤2.3:计算剩余拟候选站点所处位置的坡度,并将所有坡度值大于45的拟候选站点位置删除。
步骤2.4:剩余的所有拟候选站点将被确定为候选站点,且组成候选站点集合J,候选站点总数为N。
步骤3:输入所选TLS扫描仪各项参数(包含:扫描范围、垂直视场、入射角阈值等)。并在所有候选站点上,基于扫描仪参数以及初始三维模型为约束进行扫描模拟,筛选视线不被遮挡的滑坡体素为已观测体素,筛选视线被遮挡的滑坡体素为未观测体素。获取所有候选站点对滑坡体素的模拟覆盖情况。
具体地,在以扫描仪参数与初始三维模型为约束下进行观测模拟的过程为:
步骤3.1:依照遍历顺序选出第一个候选站点,通过限制候选站点的最大搜索范围和最小搜索范围,从而达成对扫描范围的模拟,进而确定候选站点的扫描范围;
步骤3.2:遍历所有滑坡体素,计算候选站点到滑坡体素的方向向量与观测目标点所在曲面法向量的夹角关系,筛选满足预设夹角条件的滑坡体素;
步骤3.3:计算候选站点与滑坡体素的水平距离(在扫描仪仰角已经确定且只考虑仰角的情况下,滑坡体素与候选站点的高差将直接决定满足两者可见性的最短水平距离),筛选满足预设水平距离条件的滑坡体素;
步骤3.4:将候选站点作为起点,滑坡体素作为终点,以固定步长在视线上内插点,并计算内插点的高程。将视线上内插点的高程与三维模型中同一水平坐标的高程进行比较。当出现内插点高程小于三维模型中同一水平坐标的高程时,该视线视为被遮挡并停止计算,筛选视线不被遮挡的滑坡体素为已观测体素。此方法的计算量不会随着扫描范围的增加而成指数型增长。
步骤3.5:遍历所有候选站点,并使用列表记录所有候选站点的已观测体素,已观测体素为满足步骤3.1 - 3.4的滑坡体素,不满足的为未观测体素。
本步骤中提及的约束主要包含模拟观测条件的TLS扫描仪各项参数和与初始三维模型的场景约束。本发明主要考虑的TLS扫描仪各项参数有:扫描范围(步骤3.1),由于沙石、悬崖等物体的目标反射率为60%,则扫描范围为在此目标反射率下的最大测量距离;入射角(即步骤3.2中的夹角关系间接体现入射角),即扫描仪的激光束与被观测物体表面法线的夹角;仰角(即步骤3.3中的水平距离间接体现仰角),即扫描仪的最大垂直角。扫描范围与仰角在仪器出厂时便以确定,可直接查询参数确定。而入射角影响因素颇多,本发明从回波角度考虑将最大入射角设置为75°。环境约束(步骤3.4)包含初始三维模型中所有会对扫描仪产生遮挡的物体,这其中包含观测目标自遮挡、非观测目标地物遮挡与地形遮挡。
步骤4:计算观测区域的站点基数B sp 。
通常情况下,用户依照观测区域大小主观确定需要布设的站点数量,该数量只要不超过项目承受范围都是可以接受的。但从经济角度来看,更少的站点意味着更低的观测墩布设成本和每次监测更少的扫描时间。因此在不影响扫描质量的情况下,用户更青睐于站点数量更低的规划方案。鉴于上述思路,本方法提出了“站点基数”B sp 作为观测区域站点数量的最低值,即布设数量为B sp 的候选站点时,就可以覆盖绝大多数的观测区域。同时,使用“额外站点率”E sp 来确定观测区域中允许布设的额外站点数量,进而得到观测区域中允许布设的站点区间,如公式1所示。
其中,a为候选站点的序号,N为候选站点的总数。因此,cp a 是第a个候选站点的标签,如果cp a 为1表示已选择该候选站点,而cp a 为0则表示未选择该候选站点。
在优化过程中快速得到“站点基数”B sp 非常重要。因此,本发明使用了操作性强、效率极高的贪婪算法在扫描模拟结束后立刻计算出“站点基数”。基于贪婪算法计算“站点基数”的流程可以被描述为:将所有滑坡体素设置为未覆盖,每次都选出覆盖最多未覆盖滑坡体素的候选站点,并将所选候选站点覆盖的未覆盖滑坡体素标记为已覆盖滑坡体素。当所选候选站点需要标记的未覆盖滑坡体素数量低于所设定的阈值(如:5%)时,停止迭代。此时已选出的候选站点数量为“站点基数”B sp 。
步骤5:随机生成长度为候选站点总数N的二进制编码,每一条二进制编码作为一个个体,共生成数量为n的个体组成为初始父代P t ,采用约束NSGA-II算法生成Pareto最优解集。
具体地,约束NSGA-II算法选址优化过程为:
步骤5.1:初始化种群。设定种群大小n,随机生成长度为候选站点总数N的二进制编码,其中仅有1与0组成,1代表该位置序号的候选站点处设站,0代表该位置序号的候选站点不设站,同时二进制编码中1对应的位置序号,组成该个体的所选站点序号集A,记个体的所选站点序号集A中序号的个数为u。每一条二进制编码作为一个个体,共生成数量为n的个体组成为初始父代P t 。
步骤5.2:对父代种群P t 中的个体进行交叉和变异,生成一个种群大小同样为n的子代种群Q t 。合并父代种群P t 与子代种群Q t ,得到一个为种群大小为2n的新种群R t 。
步骤5.3:计算新种群R t 中各个体(i=1,2,...,2n)的各适应度函数值,主要包含以下两个适应度函数:控制规划结果成本的适应度函数F 1,以及控制规划结果质量的适应度函数F 2。具体计算公式如式2-3所示。F 1和F 2的值越大,个体的经济性与扫描质量越好。
其中,j为滑坡体素的序号,其包含步骤3中所有已观测体素和未观测体素。N v 为滑坡体素的总数,a为候选站点的序号,N为候选站点的总数。因此,lv j 为第j个滑坡体素的标签,其中1表示已测量滑坡体素,0表示未测量滑坡体素。已测量滑坡体素表示该滑坡体素在个体中的多个候选站点中有被作为已观测体素,未测量滑坡体素表示该滑坡体素在个体的多个候选站点中没有被作为已观测体素。p 0是定义的惩罚因子,该值取值范围为[0, 1],且值越大惩罚力度越大,也意味着优化过程中越不能接受冗余站点,cp a 是第a个候选站点的标签,如果cp a 为1表示已选择该候选站点,而cp a 为0则表示未选择该候选站点,B sp 为站点基数。
其中,k为已测量滑坡体素的索引;A为个体所选候选站点序号集;θ a.k 为第a个候选站点与第k个已测量滑坡体素之间的入射角;r a.k 为第a个候选站点与第k个已测量滑坡体素之间的距离;r max为所使用的扫描仪最大射程;r min为所使用的扫描仪最大射程。
步骤5.4:基于各个体的各个适应度函数值对新种群R t 进行带约束条件的非支配排序。判断新种群R t 中所有个体之间的支配关系。
支配关系的判断方法如下具体实例,假定存在个体p与个体q,由于F 1和F 2皆为最大最优,
若个体p的F 1和F 2均分别大于个体q的F 1和F 2时,则视为p支配q,
若个体p的F 1等于个体q的F 1,个体p的F 2大于个体q的F 2,则视为p支配q,
若个体p的F 2等于个体q的F 2,个体p的F 1大于个体q的F 1,则视为p支配q,
若个体p的F 1大于个体q的F 1,而个体p的F 2小于个体q的F 2,则视为p与q为非支配关系,
若个体p的F 2大于个体q的F 2,而个体p的F 1小于个体q的F 1,则视为p与q为非支配关系。
依照个体被支配的次数,划分Pareto等级,被支配的次数越少,代表个体越优,等级越高。如果个体不满足基于方案可行性的约束条件,则跳过计算支配关系的步骤,直接视为Pareto最低等级。本发明中考虑的基于方案可行性的约束条件,如式4。
其中,O lim是用户设定的满足可配准需求的最低重叠度约束(通常为20%)。O m 是第m条配准路径的重叠度,其中m属于集合 。Cov min是用户设定的对滑坡区域最低的覆盖范围要求,E sp 为额外站点率。
其中重叠度O m 计算大致流程为:若个体的候选站点数为1,则不考虑公式4的最后一个约束。若个体的候选站点数为2,则直接计算两候选站点的重叠度。相邻两候选站点同时覆盖的滑坡体素数量与两候选站点中覆盖滑坡体素数量较低的候选站点所覆盖的滑坡体素数量的比值为两站重叠度。若个体的候选站点数大于等于3,则构建候选站点的Delaunay三角网,计算Delaunay三角网的各边的重叠度,并以重叠度倒数为权重构建最小生成树。该最小生成树为个体的配准路径,其中路径的序号为m。
步骤5.5:通过拥挤度计算对新种群R t 中同一Pareto等级中的个体进行排序。第i个体的拥挤度为第i+1个体与第i-1个体所有适应度函数值之差归一化后的和。其中第一个和最后一个个体的拥挤距离设置为无穷大。同一Pareto等级中的个体拥挤度越高,排名越靠前。第i个体的拥挤度C i 的具体计算公式如5所示。
其中:
步骤5.6:从排序后的新种群R t 中选择拥挤度靠前的n个个体组成下一代父代种群P t+1 ;
步骤5.7:重复步骤5.2 - 5.6直到迭代次数大于设定阈值。输出最后一代种群的Pareto等级最高的一组个体作为Pareto最优集。
步骤6:基于PEG-MCDM算法从Pareto最优集中决策出最终个体。
Pareto最优集中任何一个解在实际滑坡监测中都是可行的,但用户往往仅需要一个唯一解。PEG-MCDM算法是一种不需要用户提供任何权重,仅通过输入的Pareto最优集中各适应度函数值的分布情况,输出一个权衡各项目标的折中解。基于PEG-MCDM算法从Pareto最优集中决策出最终方案的具体步骤如下。
步骤6.1:将Pareto最优集中个体的适应度函数值都映射到[0,1]的范围,其中1对应于设计者的愿望。Pareto最优集中第i个体两个适应度函数值归一化后,得到Pareto最优集中归一化适应度函数值与归一化适应度函数值,具体计算公式如6所示。
其中,是Pareto最优集中第i个个体的适应度函数F 1的值;是Pareto最优集中第i个个体的适应度函数F 2的值;和为Pareto最优集中所有适应度函数值F 1的最大值与最小值;和为Pareto最优集中所有适应度函数值F 2的最大值与最小值。
步骤6.2:将Pareto最优集中个体的归一化适应度函数值从最小值到最大值顺序重新排序形成主准则向量x 1,将Pareto最优集中个体的归一化适应度函数值从最小值到最大值顺序重新排序形成主准则向量x 2,如公式7。将主准则向量x 1与主准则向量x 2反向排列,并相互交换获得主聚合向量y 1与主聚合向量y 2,如公式8。
步骤6.5:根据PEG理论估计两个适应度函数的最优值F 1 0与F 2 0,具体计算公式如11所示。
实施例2:
为说明本发明的有效性,利用实施例1所述面向滑坡体形变场时序监测的TLS多目标优化选址方法对中国湖北省武汉市江夏区的两处滑坡场景进行模拟。其中测区1滑坡隐患区域占地面积约为30850m2,研究区海拔高度为34.6m~146.7m;测区2滑坡隐患区域占地面积约为99238m2,研究区海拔高度为13.4m~112.9m。两处场景均发生过滑坡,对山地下方居住村民的生命与财产构成了威胁。因此,有必要对这两块研究区域进行监测,同时该区域也非常适合测试本发明的效果。模拟使用的扫描仪为Riegl VZ400,其扫描范围为500m、仰角为60°,设置最大水平入射角为75°。优化过程中,种群大小设为50,最大迭代次数设置为400,E sp 设置为50%。
在测区1放置候选站点的空间中划分5m的候选网格,得到5429个TLS候选站点,并进行扫描模拟。通过贪婪算法得到该区域的站点基数为1,这表明该研究区域在E sp 为50%时,允许布设的站点数量区间为[0,1]。紧接着,通过约束NSGA-II生成该测区的Pareto最优集,如图2 (a)所示,并通过PEG-MCDM输出TLS最终解,如图2(a)中的红点。图2(b)和图2(c)分别展示了成本适应度函数值与质量适应度函数值的收敛情况。其中,成本适应度函数值约在70代时收敛,质量适应度函数值约在350代时收敛。最终结果包含2个站点,其覆盖率为95.02%,精度值为0.6163,所选站点在研究区域中的分布情况如图2(d)所示。
在测区2放置候选站点的空间中划分5m的候选网格,得到10216个TLS候选站点,并进行扫描模拟。通过贪婪算法得到该区域的站点基数为3,这表明该研究区域在E sp 为50%时,允许布设的站点数量区间为[3,5]。紧接着,通过NSGA-II生成该研究区域的Pareto最优集,如图3 (a)所示,并通过PEG-MCDM输出TLS最终解,如图3 (a)中的红点。图3 (b)和图3(c)分别展示了成本适应度函数值与质量适应度函数值的收敛情况。其中,成本适应度函数值约在50代时收敛,质量适应度函数值约在200代时收敛。最终结果包含4个站点,其覆盖率为87.99%,精度值为0.4021,所选站点在研究区域中的分布情况如图3 (d)所示。
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代,但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。
Claims (9)
1.面向滑坡体形变场时序监测的TLS多目标优化选址方法,其特征在于,包括以下步骤:
步骤1:生成观测区域的三维点云,将三维点云转换为体素,根据体素建立观测区域的初始三维模型;
步骤2:在初始三维模型中确定候选站点,构成候选站点集合J,候选站点总数为N;
步骤3:确定TLS扫描仪的参数,基于扫描仪参数以及初始三维模型为约束进行扫描模拟,筛选视线不被遮挡的滑坡体素为已观测体素,筛选视线被遮挡的滑坡体素为未观测体素;
步骤4:计算观测区域的站点基数B sp ;
步骤5:随机生成长度为候选站点总数N的二进制编码,每一条二进制编码作为一个个体,共生成数量为n的个体组成为初始父代P t ,采用约束NSGA-II算法生成Pareto最优解集;
步骤6:基于PEG-MCDM算法从Pareto最优集中决策出最终个体,
所述步骤3包括以下步骤:
步骤3.1:确定TLS扫描仪的参数,依照遍历顺序选出第一个候选站点,确定候选站点的扫描范围;
步骤3.2:遍历所有滑坡体素,计算候选站点到滑坡体素的方向向量与观测目标点所在曲面法向量的夹角关系,筛选满足预设夹角条件的滑坡体素;
步骤3.3:计算候选站点与滑坡体素的水平距离,筛选满足预设水平距离条件的滑坡体素;
步骤3.4:筛选视线不被遮挡的滑坡体素为已观测体素;
步骤3.5:遍历所有候选站点,已观测体素为满足步骤3.1 - 3.4的滑坡体素,不满足的为未观测体素。
2.根据权利要求1所述面向滑坡体形变场时序监测的TLS多目标优化选址方法,其特征在于,所述步骤1包括以下步骤:
步骤1.1:用无人机获取观测区域垂直影像和倾斜影像;
步骤1.2:将无人机采集的垂直影像和倾斜影像重建为观测区域的三维点云;
步骤1.3:对观测区域的三维点云进行点云分类,将位于滑坡隐患区域内的点云标记为滑坡点,将滑坡隐患区域外的点云依照其属性划分为地面点与其他点;
步骤1.4:将空间分割为具有定义大小的三维体素,并将点云与点云类型分配给相应的体素,如果一个体素包含点云,则将保留该体素;否则,该体素将被删除,然后将保留的体素作为观测区域的初始三维模型,其中滑坡点云转换为滑坡体素的数量为N v 。
3.根据权利要求2所述面向滑坡体形变场时序监测的TLS多目标优化选址方法,其特征在于,所述步骤2包括以下步骤:
步骤2.1:在初始三维模型构建网格网络,并将网格的中心点作为拟候选站点位置;
步骤2.2:遍历所有的拟候选站点位置,如果与该拟候选站点位置距离最近体素的点云类型为其他点,则将该拟候选站点删除;
步骤2.3:计算剩余拟候选站点所处位置的坡度,并将所有坡度值大于45的拟候选站点位置删除;
步骤2.4:剩余的所有拟候选站点将被确定为候选站点,组成候选站点集合J,候选站点总数为N。
4.根据权利要求3所述面向滑坡体形变场时序监测的TLS多目标优化选址方法,其特征在于,所述步骤4包括以下步骤:
将所有滑坡体素设置为未覆盖,每次都选出覆盖最多未覆盖滑坡体素的候选站点,并将所选候选站点覆盖的未覆盖滑坡体素标记为已覆盖滑坡体素,当所选候选站点需要标记的未覆盖滑坡体素数量低于所设定的阈值时,停止迭代,此时已选出的候选站点数量为站点基数B sp 。
5.根据权利要求4所述面向滑坡体形变场时序监测的TLS多目标优化选址方法,其特征在于,所述步骤5包括以下步骤:
步骤5.1:设定种群大小n,随机生成长度为候选站点总数N的二进制编码,其中仅有1与0组成,1代表该位置序号的候选站点处设站,0代表该位置序号的候选站点不设站,同时二进制编码中1对应的位置序号,组成该个体的所选站点序号集A,记个体的所选站点序号集A中序号的个数为u,每一条二进制编码作为一个个体,共生成数量为n的个体组成为初始父代P t ;
步骤5.2:对父代种群P t 中的个体进行交叉和变异,生成一个种群大小同样为n的子代种群Q t ,合并父代种群P t 与子代种群Q t ,得到种群大小为2n的新种群R t ;
步骤5.3:计算新种群R t 中各个体的控制规划结果成本的适应度函数F 1和控制规划结果质量的适应度函数F 2;
步骤5.4:基于各个体的控制规划结果成本的适应度函数F 1和控制规划结果质量的适应度函数F 2对新种群R t 进行非支配排序,判断新种群R t 中所有个体之间的支配关系;
步骤5.5:通过拥挤度对新种群R t 中同一Pareto等级中的个体进行排序,同一Pareto等级中的个体拥挤度越高,排名越靠前;
步骤5.6:从排序后的新种群R t 中选择拥挤度靠前的n个个体组成下一代父代种群P t+1 ;
步骤5.7:重复步骤5.2 - 5.6直到迭代次数大于设定阈值,输出最后一代种群的Pareto等级最高的一组个体作为Pareto最优集。
6.根据权利要求5所述面向滑坡体形变场时序监测的TLS多目标优化选址方法,其特征在于,所述步骤5.3包括以下步骤:
其中,j为滑坡体素的序号,N v 为滑坡体素的总数,a为候选站点的序号,N为候选站点的总数,lv j 为第j个滑坡体素的标签,已测量滑坡体素表示该滑坡体素在个体中的多个候选站点中有被作为已观测体素,未测量滑坡体素表示该滑坡体素在个体的多个候选站点中没有被作为已观测体素,p 0是惩罚因子,cp a 是第a个候选站点的标签,如果cp a 为1表示已选择该候选站点,而cp a 为0则表示未选择该候选站点,B sp 为站点基数;
其中,k为已测量滑坡体素的索引;A为个体所选候选站点序号集;θ a.k 为第a个候选站点与第k个已测量滑坡体素之间的入射角;r a.k 为第a个候选站点与第k个已测量滑坡体素之间的距离;r max为所使用的扫描仪最大射程;r min为所使用的扫描仪最小射程。
7.根据权利要求6所述面向滑坡体形变场时序监测的TLS多目标优化选址方法,其特征在于,所述步骤5.4包括以下步骤:
假定个体 p与个体q,
若个体p的F 1和F 2均分别大于个体q的F 1和F 2时,则视为p支配q,
若个体p的F 1等于个体q的F 1,个体p的F 2大于个体q的F 2,则视为p支配q,
若个体p的F 2等于个体q的F 2,个体p的F 1大于个体q的F 1,则视为p支配q,
若个体p的F 1大于个体q的F 1,而个体p的F 2小于个体q的F 2,则视为p与q为非支配关系,
若个体p的F 2大于个体q的F 2,而个体p的F 1小于个体q的F 1,则视为p与q为非支配关系,
依照个体被支配的次数,划分Pareto等级,被支配的次数越少,个体越优,等级越高,个体不满足约束条件,则直接视为Pareto最低等级。
重叠度O m 计算包括以下步骤:
若个体的候选站点数为1,则不考虑O lim约束;
若个体的候选站点数为2,则直接计算两候选站点的重叠度,相邻两候选站点同时覆盖的滑坡体素数量与两候选站点中覆盖滑坡体素数量较低的候选站点所覆盖的滑坡体素数量的比值为两站重叠度;
若个体的候选站点数大于等于3,则构建候选站点的Delaunay三角网,计算Delaunay三角网的各边的重叠度,并以重叠度倒数为权重构建最小生成树,该最小生成树为个体的配准路径,其中路径的序号为m。
9.根据权利要求8所述面向滑坡体形变场时序监测的TLS多目标优化选址方法,其特征在于,所述步骤6包括以下步骤:
步骤6.1:将Pareto最优集中个体的适应度函数值都映射到[0,1]的范围,其中1对应于设计者的愿望,计算Pareto最优集中归一化适应度函数值X i,1与归一化适应度函数值X i,2;
步骤6.2:将Pareto最优集中个体的归一化适应度函数值X i,1从最小值到最大值顺序重新排序形成主准则向量x 1,将Pareto最优集中个体的归一化适应度函数值X i,2从最小值到最大值顺序重新排序X i,2形成主准则向量x 2,将主准则向量x 1与主准则向量x 2反向排列,并相互交换获得主聚合向量y 1与主聚合向量y 2;
步骤6.5:根据PEG理论估计两个适应度函数的最优值F 1 0与F 2 0;
步骤6.6:计算Pareto最优集中每个个体I i 的均方误差MSE i ,MSE i 最低的个体视为最终个体进行输出。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202310053117.5A CN116125490B (zh) | 2023-02-03 | 2023-02-03 | 面向滑坡体形变场时序监测的tls多目标优化选址方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202310053117.5A CN116125490B (zh) | 2023-02-03 | 2023-02-03 | 面向滑坡体形变场时序监测的tls多目标优化选址方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN116125490A CN116125490A (zh) | 2023-05-16 |
CN116125490B true CN116125490B (zh) | 2023-07-04 |
Family
ID=86306049
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202310053117.5A Active CN116125490B (zh) | 2023-02-03 | 2023-02-03 | 面向滑坡体形变场时序监测的tls多目标优化选址方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN116125490B (zh) |
Family Cites Families (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2007069724A1 (ja) * | 2005-12-16 | 2007-06-21 | Ihi Corporation | 三次元形状データの位置合わせ方法と装置 |
JP5631688B2 (ja) * | 2010-10-13 | 2014-11-26 | 国土防災技術株式会社 | 地上3dレーザースキャナーを用いた、ターゲット無設置による斜面の変動監視方法。 |
JP2013096745A (ja) * | 2011-10-28 | 2013-05-20 | Hokuriku Kosoku Co Ltd | 三次元モデルの作成方法 |
KR101507333B1 (ko) * | 2014-05-14 | 2015-03-31 | 숭실대학교산학협력단 | 경사면 상황 인식 방법, 이를 이용한 시스템 및 이를 수행하기 위한 기록매체 |
JP7076372B2 (ja) * | 2015-10-12 | 2022-05-27 | グラウンドプルーブ・ピーティーワイ・リミテッド | 斜面安定性ライダー |
CN106951994A (zh) * | 2017-03-21 | 2017-07-14 | 武汉理工大学 | 一种海上应急救援站点的选址方法 |
CN110443002A (zh) * | 2019-08-16 | 2019-11-12 | 中国水利水电科学研究院 | 一种高边坡形变预测方法及系统 |
CN114549601B (zh) * | 2022-02-11 | 2023-03-28 | 中国科学院精密测量科学与技术创新研究院 | 一种顾及点对可靠性的滑坡多时相tls点云精配准方法 |
CN115166800A (zh) * | 2022-06-01 | 2022-10-11 | 河北工业大学 | 联合gnss和三维激光扫描的滑坡监测数据的处理方法 |
-
2023
- 2023-02-03 CN CN202310053117.5A patent/CN116125490B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN116125490A (zh) | 2023-05-16 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN105160702B (zh) | 基于LiDAR点云辅助的立体影像密集匹配方法及系统 | |
CN101105396B (zh) | 3d扫描数据自动对齐的系统和方法 | |
CN106846308B (zh) | 基于点云的地形图精度的检测方法和装置 | |
Király et al. | Tree height estimation methods for terrestrial laser scanning in a forest reserve | |
CN100485662C (zh) | 基于数据动态存取模型的产品点云型面特征分析方法 | |
CN102096072B (zh) | 一种城市部件自动化测量方法 | |
CN112347550A (zh) | 耦合式室内三维语义建图及建模方法 | |
CN113916130B (zh) | 一种基于最小二乘法的建筑物位置测量方法 | |
CN109348416B (zh) | 基于二分k均值的指纹室内定位方法 | |
CN102054166A (zh) | 一种新的用于户外增强现实系统的场景识别技术 | |
CN110047133A (zh) | 一种面向点云数据的列车边界提取方法 | |
CN115761303A (zh) | 基于机载激光雷达点云和遥感影像数据的地物分类方法 | |
CN112258624A (zh) | 一种三维实景融合建模的方法 | |
CN110114802A (zh) | 地形信息处理装置、地形信息处理方法及程序 | |
CN115187648A (zh) | 输电线路本体逆向建模方法、装置、电子设备及存储介质 | |
CN117609774A (zh) | 一种基于配电网的三维场景点云数据提取与识别系统 | |
CN116125490B (zh) | 面向滑坡体形变场时序监测的tls多目标优化选址方法 | |
CN117876465A (zh) | 一种基于落水洞内部点云数据的体积计算方法、系统、设备和介质 | |
Metawie et al. | Optimizing laser scanning positions in buildings exteriors: heritage building application | |
Ma et al. | Hybrid model for realistic and efficient estimation of highway sight distance using airborne LiDAR data | |
CN117218063A (zh) | 基于施工阶段吊篮空间位置变化的桥体施工进度监测方法 | |
CN114330526B (zh) | 基于tgo-ssa-fcm算法的交通流量缺失数据修复方法 | |
CN115018973B (zh) | 一种低空无人机点云建模精度的无靶标评估方法 | |
CN108986217B (zh) | 基于图型矢量距离的多点地质统计学建模方法 | |
CN115761162A (zh) | 一种气象站点等值面构造方法、装置及存储介质 |
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 |