CN117523124B - 一种基于无人机激光雷达的大流域地形变化切片方法 - Google Patents

一种基于无人机激光雷达的大流域地形变化切片方法 Download PDF

Info

Publication number
CN117523124B
CN117523124B CN202311507747.1A CN202311507747A CN117523124B CN 117523124 B CN117523124 B CN 117523124B CN 202311507747 A CN202311507747 A CN 202311507747A CN 117523124 B CN117523124 B CN 117523124B
Authority
CN
China
Prior art keywords
point
point cloud
computer
slice
obtaining
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
Application number
CN202311507747.1A
Other languages
English (en)
Other versions
CN117523124A (zh
Inventor
李朋飞
杨鑫
高晨迪
胡晋飞
白晓
严露
李豆
丹杨
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian University of Science and Technology
Original Assignee
Xian University of Science and Technology
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Xian University of Science and Technology filed Critical Xian University of Science and Technology
Priority to CN202311507747.1A priority Critical patent/CN117523124B/zh
Publication of CN117523124A publication Critical patent/CN117523124A/zh
Application granted granted Critical
Publication of CN117523124B publication Critical patent/CN117523124B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05—Geographic models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/60—Analysis of geometric attributes
    • G06T7/62—Analysis of geometric attributes of area, perimeter, diameter or volume
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2200/00—Indexing scheme for image data processing or generation, in general
    • G06T2200/08—Indexing scheme for image data processing or generation, in general involving all processing steps from image acquisition to 3D model generation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10028—Range image; Depth image; 3D point clouds

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Software Systems (AREA)
  • Remote Sensing (AREA)
  • Electromagnetism (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Computer Graphics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Optical Radar Systems And Details Thereof (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种基于无人机激光雷达的大流域地形变化切片方法，该方法包括以下步骤：一、待测流域的点云获取；二、待测流域的点云的切片；三、待测流域点云切片的收缩；四、获取后一期点云切片相对前一期点云切片的切片变化面积；五、待测流域的体积变化量获取。本发明设计合理，获取沿待测流域的主沟道延伸方向各个点云切片，对各个点云切片进行拉普拉斯收缩，并剔除空间分布式LoD后获取切片面积，以提高获取切片体积变化的精度，有效地适应流域地形变化监测。

Description

一种基于无人机激光雷达的大流域地形变化切片方法
技术领域
本发明属于流域体积变化测量技术领域，尤其是涉及一种基于无人机激光雷达的大流域地形变化切片方法。
背景技术
流域尺度是指地理学和水文学中研究流域特征和过程的尺度范围，流域是指有一片地表的区域组成，其内部的水流汇集到一个共同的出口点，流域尺度涵盖了从小到大不同层次的地理范围，例如小流域、中等流域和大流域。
高精度地形变化监测为流域区域变化、地质灾害、地表沉降等过程研究提供了有效手段。目前切片体积方法对地形点云进行切片，利用模糊C均值聚类方法获取切片点云轮廓，并利用点云模型的LoD均值剔除地形变化点云的不确定性误差，该方法利用聚类方法得到的点云切片点较少，但由于沟道侵蚀地形变化体积较大，导致切片点云轮廓误差以及LoD计算误差对整体精度影响较小，因此该方法在小尺度沟道侵蚀地形变化监测中能够达到较高的精度，但在大尺度地形变化监测中，地形变化量较小，聚类方法难以描述微变化地形的细节信息，并且均值LoD计算方法难以捕获到细微地形变化且易放大重大地形变化，精度较低，因此传统切片体积方法在大流域地形变化监测中并不适用。
发明内容
本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种基于无人机激光雷达的大流域地形变化切片方法，其方法步骤简单，设计合理，获取沿待测流域的主沟道延伸方向各个点云切片，对各个点云切片进行拉普拉斯收缩，并剔除空间分布式LoD后获取切片面积，以提高获取切片体积变化的精度，有效地适应流域地形变化监测。
为解决上述技术问题，本发明采用的技术方案是：一种基于无人机激光雷达的大流域地形变化切片方法，其特征在于，该方法包括以下步骤：
步骤一、待测流域的点云获取：
采用无人机搭载三维激光雷达系统对待测流域进行扫描，得到第q期初始激光点云和第q+1期初始激光点云；其中，q为正整数，且1≤q；第q期初始激光点云和第q+1期初始激光点云的点云坐标为CGCS2000坐标系；
步骤二、待测流域的点云的切片：
步骤201、采用计算机建立OXYZ坐标系；其中，OXYZ坐标系中X轴的正方向沿待测流域的主沟道延伸方向，OXYZ坐标系中Y轴垂直X轴，OXYZ坐标系中Z轴垂直X轴和Y轴垂直向上；
步骤202、采用计算机将第q期初始激光点云和第q+1期初始激光点云转换为OXYZ坐标系下的点云坐标，得到第q期激光点云和第q+1期激光点云；
步骤203、当X轴的正方向沿待测流域的主沟道延伸方向，在OXYZ坐标系下，采用计算机按照相邻两期的点云切片厚度δ(q,q+1)，将第q期激光点云沿X轴方向做多个切平面并将多个切平面沿X轴正方向排序依次记作第1个切平面，...，第j个切平面，...，第J个切平面；其中，j和J均为正整数，且1≤j≤J；任一个切平面垂直X轴，第1个切平面和第J个切平面之间相邻两个切平面的间距为δ(q,q+1)，第1个切平面距离第q期激光点云与第q+1期激光点云中最小X轴坐标处的间距为δ(q,q+1)/2，第J个切平面距离第q期激光点云与第q+1期激光点云中最大X轴坐标处的间距记作L′(J,J-1),且L′(J,J-1)不大于δ(q,q+1)；
步骤三、待测流域点云切片的收缩：
步骤301、当X轴的正方向沿待测流域的主沟道延伸方向，当j处于1～J-1时，采用计算机将第j个切平面左侧及右侧X轴方向δ(q,q+1)/2范围内的点云投影至第j个切平面，得到第q期第j个点云切片；
当j取J时，采用计算机将第J个切平面右侧X轴方向的剩余点云及左侧X轴方向δ(q,q+1)/2范围内的点云投影至第J个切平面，得到第q期第J个点云切片；
步骤302、采用计算机获取第q期第j个点云切片上投影点的总数并记作P；
步骤303、采用计算机利用K近邻方法对第q期第j个点云切片上的点云进行处理，得到第q期第j个点云切片上第p个投影点的K邻域Np，再对Np进行狄洛尼三角剖分，得到第q期第j个点云切片上第p个投影点的单环邻域点集合；其中，第q期第j个点云切片上第p个投影点的单环邻域点集合记作第p个单环邻域点集合；p为正整数，且1≤p≤P；
步骤304、采用计算机将第p个投影点和第p个投影点的单环邻域点集合记作第p个投影点集合，并对第q期第j个点云切片上第p个投影点集合进行处理，得到余切权拉普拉斯矩阵L；
步骤305、采用计算机根据余切权拉普拉斯矩阵L，对第p个投影点集合进行收缩，得到第p个投影点集合中各个收缩后的坐标，提取得到第p个投影点收缩后的坐标；
步骤306、多次重复步骤303至步骤305的方法，完成第q期第j个点云切片上各个投影点的收缩，将收缩后的各个投影点组合为第q期第j个理想点云切片；
步骤307、按照步骤301至步骤306的方法，完成第q+1期第j个点云切片上各个投影点的收缩，得到第q+1期第j个理想点云切片；
步骤四、获取后一期点云切片相对前一期点云切片的切片变化面积：
步骤401、采用计算机将第q期第j个理想点云切片中距离最远的两个点中靠近原点的点作为第一个点，另一个作为终点，利用贪心算法以两点距离最近对第q期第j个理想点云切片中各个点进行排序，得到第q期第j个有序点云；
步骤402、按照步骤401的方法，得到第q+1期第j个有序点云；
步骤403、当X轴的正方向沿待测流域的主沟道延伸方向，采用计算机将第q期第j个有序点云和第q+1期第j个有序点云/>同步绘制在平面OYZ上，且/>和/>中相交的相邻两个交点将/>和/>之间围设的区域划分为第1个切平面区块，..，第e个切平面区块，..，第Ej个切平面区块；其中，e和Ej均为正整数，且1≤e≤Ej，Ej表示第j个点云切片上切平面区块的总数；
步骤404、获取剔除空间分布式不确定性值LoD后第j个有序点云上第e个切平面区块的面积记作
步骤五、待测流域的体积变化量获取：
步骤501、根据第j个有序点云上第e个切平面区块的面积，获取第j个有序点云上第e个切平面区块的变化体积，进而得到第j个点云切片的变化体积；
步骤502、将各个点云切片的变化体积累加得到待测流域区域第q+1期相对第q期的体积变化量V(q,q+1)。
上述的一种基于无人机激光雷达的大流域地形变化切片方法，其特征在于：步骤一中采用无人机搭载三维激光雷达系统对待测流域进行扫描，获取第q期第1架次激光点云，...，第q期第f架次的激光点云，...，第q期第F架次的激光点云；其中，f和F均为正整数，1≤f≤F，F≥5；
步骤102、采用计算机利用CloudCompare软件中“Merge”工具，对第q期第1架次激光点云至第q期第F架次的激光点云进行合并，得到合并后的第q期激光点云；
步骤103、采用计算机利用CloudCompare软件中“Segment”工具，按照待测流域实地野外踏勘结果对合并后的第q期激光点云进行裁剪，得到第q期待测流域点云；
步骤104、采用计算机利用MCC点云滤波方法对第q期待测流域点云进行一次滤波，得到一次滤波后第q期激光点云；
步骤105、采用计算机利用TerraSolid软件对一次滤波后第q期激光点云进行二次滤波，得到二次滤波后第q期激光点云，并记作第q期初始激光点云；
步骤106、按照步骤101至步骤105所述的方法，得到第q+1期初始激光点云。
上述的一种基于无人机激光雷达的大流域地形变化切片方法，其特征在于：步骤二中相邻两期的点云切片厚度δ(q,q+1)的获取，具体过程如下：
步骤201、采用计算机从第q期激光点云中随机选择点云组成第q个点云集合；其中，第q个点云集合中第i个点记作1≤i≤I，I表示第q个点云集合的总数；
步骤202、采用计算机获取第q个点云集合中第i个点与第q期激光点云中各个点的欧式距离，并将各个欧式距离按照从小到大的顺序排序，获取前M个欧式距离；其中，与/>对应的前个欧式距离中第m个欧式距离记作/>m和M均为正整数，且1≤m≤M；
步骤203、根据公式得到第q期激光点云的平均点间距ρ(q)；
步骤204、按照步骤201至步骤203所述的方法，得到第q+1期的激光点云的平均点间距ρ(q+1)；并对ρ(q)和ρ(q+1)进行平均值处理，得到相邻两期的平均点间距ρ(q,q+1)；
步骤205、根据公式δ(q,q+1)＝B×ρ(q,q+1)，得到相邻两期的点云切片厚度δ(q,q+1)；其中，B为常数且B取值为0.5～1.5。
上述的一种基于无人机激光雷达的大流域地形变化切片方法，其特征在于：步骤304对第q期第j个点云切片上第p个投影点集合进行处理，得到余切权拉普拉斯矩阵L，具体过程如下：
步骤3041、采用计算机将第p个投影点记作第1个点，将第p个单环邻域点集合中点云按照逆时针顺序标记为第2个点，...，第g个点；其中，将第1个点至第g个点中任一个点记作第a个点或者第b个点，且a和b均为正整数且取值为1～g；
步骤3042、当a和b取值位于2～g，且b大于a，a和b为相邻两个点时，采用计算机根据公式ωab＝cotαab，得到余切权拉普拉斯矩阵L的第a行第b列元素值ωab；其中，αab为第1个点和第a个点连线与第1个点和第b个点连线之间的夹角，g和2为相邻两个点；并根据ωba＝ωab，得到余切权拉普拉斯矩阵L的第b行第a列元素值ωba；
当a取值为1，当b取值为3～g-2时，采用计算机根据公式ω1b＝cotα1b+cotβ1b，得到余切权拉普拉斯矩阵L的第1行第b列元素值ω1b；其中，α1b为第1个点和第b-1个点连线与第b-1个点和第b个点连线之间的夹角；βab为第1个点和第b+1个点连线与第b+1个点和第b个点连线之间的夹角；并根据ω1b＝ωb1，得到余切权拉普拉斯矩阵L的第b行第1列元素值ωb1；
当a取值为1，当b取值为2时，采用计算机根据公式ω12＝cotα12+cotβ12，得到余切权拉普拉斯矩阵L的第1行第2列元素值ω12；其中，α12为第1个点和第g个点连线与第g个点和第2个点连线之间的夹角；β12为第1个点和第3个点连线与第3个点和第2个点连线之间的夹角；并根据ω12＝ω21，得到余切权拉普拉斯矩阵L的第2行第1列元素值ω21；
当a取值为1，当b取值为g-1时，采用计算机根据公式ω1(g-1)＝cotα1(g-1)+cotβ1(g-1)，得到余切权拉普拉斯矩阵L的第1行第g-1列元素值ω1(g-1)；其中，α1(g-1)为第1个点和第g-2个点连线与第g-2个点和第g-1个点连线之间的夹角；β1(g-1)为第1个点和第g个点连线与第g个点和第g-1个点连线之间的夹角；并根据ω1(g-1)＝ω(g-1)1，得到余切权拉普拉斯矩阵L的第g-1行第1列元素值ω(g-1)1；
当a取值为1，当b取值为g时，采用计算机根据公式ω1g＝cotα1g+cotβ1g，得到余切权拉普拉斯矩阵L的第1行第g列元素值ω1g；其中，α1g为第1个点和第g-1个点连线与第g-1个点和第g个点连线之间的夹角；β1g为第1个点和第2个点连线与第2个点和第g个点连线之间的夹角；并根据ω1g＝ωg1，得到余切权拉普拉斯矩阵L的第g行第1列元素值ωg1；
当a和b取值为1，采用计算机根据公式
当a和b取值相同且均为2时，采用计算机根据公式ω22＝-(ω21+ω23+ω2g)，得到余切权拉普拉斯矩阵L的第2行第2列元素值ω22；其中，ω21＝ω12，ω23＝ω32；
当a和b取值相同且均为3～g-2时，采用计算机根据公式ωaa＝-(ωa1+ωa(b-1)+ωa(b+1))，得到余切权拉普拉斯矩阵L的第a行第a列元素值ωaa；其中，ωa1＝ω1a；
当a和b取值相同且均为g-1时，采用计算机根据公式ω(g-1)(g-1)＝-(ω(g-1)1+ω(g-1)(b-1)+ω(g-1)(b+1))，得到余切权拉普拉斯矩阵L的第g-1行第g-1列元素值ω(g-1)(g-1)；其中，ω(g-1)1＝ω1(g-1)；
当a和b取值相同且均为g时，采用计算机根据公式ωgg＝-(ωg1+ωg(g-1)+ωg2)；其中，ωg1＝ω1g，ωg2＝ω2g；
步骤3043、将其余元素值设置为零，则得到余切权拉普拉斯矩阵L；其中，L为对称矩阵。
上述的一种基于无人机激光雷达的大流域地形变化切片方法，其特征在于：步骤305，具体过程如下：
步骤3051、采用计算机根据公式得到第t+1次第p个投影点集合收缩后的坐标矩阵Pt+1；其中，t为大于等于0的自然数，Wl t为第t次更新的第一对角矩阵，/>为第t次更新的第二对角矩阵，Lt为第t次第p个投影点集合收缩后的余切权拉普拉斯矩阵，Pt为第t次第p个投影点集合收缩后的坐标矩阵，O表示零矩阵；当t取值为0时，Wl 0为初始的第一对角矩阵，且对角线元素值为1；/>为初始的第二对角矩阵，且对角线元素值为S0表示初始的第p个单环邻域点集合围成的邻域面积，P0为初始的第p个投影点集合的坐标矩阵；L0为初始的余切权拉普拉斯矩阵即L；
步骤3052、采用计算机根据公式Wl t+1＝Sl×Wl t，得到第t+1次第一对角矩阵的更新矩阵Wl t+1；其中，Sl取值为3.0；
采用计算机根据公式得到第t+1次第二对角矩阵的更新矩阵其中，St表示第t次收缩后第p个单环邻域点集合围成的邻域面积；
步骤3053、多次重复步骤3051和步骤3052，直至到达设定的迭代次数T，则得到第T次第p个投影点结合收缩后的坐标矩阵PT，并通过PT得到第p个投影点集合中各个收缩后的坐标。
上述的一种基于无人机激光雷达的大流域地形变化切片方法，其特征在于：步骤404中获取剔除空间分布式不确定性值LoD后第j个有序点云上第e个切平面区块的面积记作，具体过程如下：
步骤4041、采用计算机第j个有序点云上第e个切平面区块的Y轴方向最大值Ymax和Y轴方向最小值Ymin，从Y轴方向最小值开始以步长Δx，沿Y轴正方向依次设置多条直线R；
步骤4042、将第i＇条直线记作直线Ri＇，且Ri＇＝Ymin+n×Δx；其中，n表示直线数量，i＇为正整数；且n＝Ceiling[(Ymax-Ymin)/Δx]；1≤i＇≤n；
步骤4043、采用计算机将第q期第j个有序点云上距离直线Ri＇最近且分布在直线Ri＇两侧的两个点的连线和直线Ri＇的交点记作第一个交点，将第q+1期第j个有序点云上距离直线Ri＇最近且分布在直线Ri＇两侧的两个点的连线和直线Ri＇的交点记作第二个交点，采用计算机将第一个交点和第二个交点之间的距离记作距离Hi＇；
步骤4044、采用计算机利用CloudCompare软件获取距离Hi＇的不确定性值LoD(i＇)；
步骤4045、采用计算机将距离Hi＇和LoD(i＇)进行判断修正，如果距离Hi＇的绝对值小于等于LoD(i＇)，则距离Hi＇的修正值为Hi″＝0；如果距离Hi＇为负值且绝对值大于LoD(i＇)，则距离Hi＇的修正值为Hi″＝Hi＇+LoD(i＇)；如果距离Hi＇为正值且绝对值大于LoD(i＇)，则距离Hi＇的修正值为Hi″＝Hi＇-LoD(i＇)；
步骤4046、采用计算机根据公式S(i＇)＝Hi″×Δx，得到第i＇个直线处面积；
步骤4047、多次重复步骤4042至步骤4046，得到各个直线处面积，并将各个直线处面积求和，得到第j个有序点云上第e个切平面区块的面积
上述的一种基于无人机激光雷达的大流域地形变化切片方法，其特征在于：步骤4034中采用计算机利用CloudCompare软件获取距离Hi＇的不确定性值LoD(i＇)，具体过程如下：
步骤A、采用计算机利用CloudCompare软件“M3C2 distance”工具，通过输入第q期激光点云和第q+1期激光点云进行处理，得到带有不确定性误差属性值的M3C2点云；其中，所述M3C2点云中每个点具有不确定性值；
步骤B、采用计算机获取第一个交点和M3C2点云中所有点的欧式距离，并筛选出距离第一个交点小于0.5m的点集D；
步骤C、采用计算机根据公式得到第一个交点的不确定性值LoD，并记作距离Hi＇的不确定性值LoD(i＇)；其中，zk表示点集D中第k个点的不确定性值，dk表示点集D中第k个点距离第一个交点的欧式距离，k为正整数，且1≤k≤K＇，K＇表示点集D中点的总数。
上述的一种基于无人机激光雷达的大流域地形变化切片方法，其特征在于：步骤501中根据第j个有序点云上第e个切平面区块的面积，获取第j个有序点云上第e个切平面区块的变化体积，进而得到第j个点云切片的变化体积；
步骤5011、当j处于1～J-1时，则根据公式得到第j个点云切片上第e个切平面区块的变化体积/>
当j取J时，则根据公式得到第J个点云切片上第e个切平面区块的变化体积/>
步骤5012、当j处于1～J-1时，采用计算机根据得到第j个点云切片的变化体积Vj(q,q+1)；
当j取J时，采用计算机根据得到第J个点云切片的变化体积VJ(q,q+1)；
步骤502中将各个点云切片的变化体积累加得到待测流域区域第q+1期相对第q期的体积变化量V(q,q+1)，具体过程如下：
采用计算机根据得到待测流域区域第q+1期相对第q期的体积变化量V(q,q+1)。
本发明与现有技术相比具有以下优点：
1、本发明方法步骤简单、设计合理且实现方便，精度高。
2、本发明采用无人机搭载三维激光雷达系统对待测流域进行扫描，获取多期激光点云，便于后续基于激光点云进行相邻两期流域体积变化量计算。
3、本发明对任一期的激光点云沿X轴主沟道延伸方向做切平面，并将切平面沿待测流域的主沟道延伸方向点云切片厚度范围内的点云投影至切平面上，形成点云切片，从而将流域变形点云的三维问题简化为二维问题，进而有效获取和区分每个点云切片的正负地形变化量，实现复杂流域地形中各个点云切片体积变化的准确量化。
4、本发明采用余切权拉普拉斯矩阵对点云切片上的点进行收缩，是为了得到理想点云切片，进而便于后续利用贪心算法对理想点云切片中各个点进行排序，得到有序点云，从而提供更加准确的切片点云，提高了切片体积精度。
5、本发明在获取各个切平面区块的面积时，采用空间分布式LoD计算方法克服了传统切片体积方法使用统一均值LoD计算方法导致部分地形LoD计算误差远高于真实地形变化的问题，提高了切片体积方法的精度。
综上所述，本发明方法方法步骤简单，设计合理，获取沿待测流域的主沟道延伸方向各个点云切片，对各个点云切片进行拉普拉斯收缩，并剔除空间分布式LoD后获取切片面积，以提高获取切片体积变化的精度，有效地适应流域地形变化监测。
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
图1为本发明的方法流程图。
图2为本发明第p个投影点的单环邻域点集合的示意图。
图3为本发明距离Hi＇的示意图。
具体实施方式
如图1至图3所示的一种基于无人机激光雷达的大流域地形变化切片方法，包括以下步骤：
步骤一、待测流域的点云获取：
采用无人机搭载三维激光雷达系统对待测流域进行扫描，得到第q期初始激光点云和第q+1期初始激光点云；其中，q为正整数，且1≤q；第q期初始激光点云和第q+1期初始激光点云的点云坐标为CGCS2000坐标系；
步骤二、待测流域的点云的切片：
步骤201、采用计算机建立OXYZ坐标系；其中，OXYZ坐标系中X轴的正方向沿待测流域的主沟道延伸方向，OXYZ坐标系中Y轴垂直X轴，OXYZ坐标系中Z轴垂直X轴和Y轴垂直向上；
步骤202、采用计算机将第q期初始激光点云和第q+1期初始激光点云转换为OXYZ坐标系下的点云坐标，得到第q期激光点云和第q+1期激光点云；
步骤203、当X轴的正方向沿待测流域的主沟道延伸方向，在OXYZ坐标系下，采用计算机按照相邻两期的点云切片厚度δ(q,q+1)，将第q期激光点云沿X轴方向做多个切平面并将多个切平面沿X轴正方向排序依次记作第1个切平面，...，第j个切平面，...，第J个切平面；其中，j和J均为正整数，且1≤j≤J；任一个切平面垂直X轴，第1个切平面和第J个切平面之间相邻两个切平面的间距为δ(q,q+1)，第1个切平面距离第q期激光点云与第q+1期激光点云中最小X轴坐标处的间距为δ(q,q+1)/2，第J个切平面距离第q期激光点云与第q+1期激光点云中最大X轴坐标处的间距记作L′(J,J-1),且L′(J,J-1)不大于δ(q,q+1)；
步骤三、待测流域点云切片的收缩：
步骤301、当X轴的正方向沿待测流域的主沟道延伸方向，当j处于1～J-1时，采用计算机将第j个切平面左侧及右侧X轴方向δ(q,q+1)/2范围内的点云投影至第j个切平面，得到第q期第j个点云切片；
当j取J时，采用计算机将第J个切平面右侧X轴方向的剩余点云及左侧X轴方向δ(q,q+1)/2范围内的点云投影至第J个切平面，得到第q期第J个点云切片；
步骤302、采用计算机获取第q期第j个点云切片上投影点的总数并记作P；
步骤303、采用计算机利用K近邻方法对第q期第j个点云切片上的点云进行处理，得到第q期第j个点云切片上第p个投影点的K邻域Np，再对Np进行狄洛尼三角剖分，得到第q期第j个点云切片上第p个投影点的单环邻域点集合；其中，第q期第j个点云切片上第p个投影点的单环邻域点集合记作第p个单环邻域点集合；p为正整数，且1≤p≤P；
步骤304、采用计算机将第p个投影点和第p个投影点的单环邻域点集合记作第p个投影点集合，并对第q期第j个点云切片上第p个投影点集合进行处理，得到余切权拉普拉斯矩阵L；
步骤305、采用计算机根据余切权拉普拉斯矩阵L，对第p个投影点集合进行收缩，得到第p个投影点集合中各个收缩后的坐标，提取得到第p个投影点收缩后的坐标；
步骤306、多次重复步骤303至步骤305的方法，完成第q期第j个点云切片上各个投影点的收缩，将收缩后的各个投影点组合为第q期第j个理想点云切片；
步骤307、按照步骤301至步骤306的方法，完成第q+1期第j个点云切片上各个投影点的收缩，得到第q+1期第j个理想点云切片；
步骤四、获取后一期点云切片相对前一期点云切片的切片变化面积：
步骤401、采用计算机将第q期第j个理想点云切片中距离最远的两个点中靠近原点的点作为第一个点，另一个作为终点，利用贪心算法以两点距离最近对第q期第j个理想点云切片中各个点进行排序，得到第q期第j个有序点云；
步骤402、按照步骤401的方法，得到第q+1期第j个有序点云；
步骤403、当X轴的正方向沿待测流域的主沟道延伸方向，采用计算机将第q期第j个有序点云和第q+1期第j个有序点云/>同步绘制在平面OYZ上，且/>和/>中相交的相邻两个交点将/>和/>之间围设的区域划分为第1个切平面区块，..，第e个切平面区块，..，第Ej个切平面区块；其中，e和Ej均为正整数，且1≤e≤Ej，Ej表示第j个点云切片上切平面区块的总数；
步骤404、获取剔除空间分布式不确定性值LoD后第j个有序点云上第e个切平面区块的面积记作
步骤五、待测流域的体积变化量获取：
步骤501、根据第j个有序点云上第e个切平面区块的面积，获取第j个有序点云上第e个切平面区块的变化体积，进而得到第j个点云切片的变化体积；
步骤502、将各个点云切片的变化体积累加得到待测流域区域第q+1期相对第q期的体积变化量V(q,q+1)。
本实施例中，步骤一中采用无人机搭载三维激光雷达系统对待测流域进行扫描，获取第q期第1架次激光点云，...，第q期第f架次的激光点云，...，第q期第F架次的激光点云；其中，f和F均为正整数，1≤f≤F，F≥5；
步骤102、采用计算机利用CloudCompare软件中“Merge”工具，对第q期第1架次激光点云至第q期第F架次的激光点云进行合并，得到合并后的第q期激光点云；
步骤103、采用计算机利用CloudCompare软件中“Segment”工具，按照待测流域实地野外踏勘结果对合并后的第q期激光点云进行裁剪，得到第q期待测流域点云；
步骤104、采用计算机利用MCC点云滤波方法对第q期待测流域点云进行一次滤波，得到一次滤波后第q期激光点云；
步骤105、采用计算机利用TerraSolid软件对一次滤波后第q期激光点云进行二次滤波，得到二次滤波后第q期激光点云，并记作第q期初始激光点云；
步骤106、按照步骤101至步骤105所述的方法，得到第q+1期初始激光点云。
本实施例中，步骤二中相邻两期的点云切片厚度δ(q,q+1)的获取，具体过程如下：
步骤201、采用计算机从第q期激光点云中随机选择点云组成第q个点云集合；其中，第q个点云集合中第i个点记作1≤i≤I，I表示第q个点云集合的总数；
步骤202、采用计算机获取第q个点云集合中第i个点与第q期激光点云中各个点的欧式距离，并将各个欧式距离按照从小到大的顺序排序，获取前M个欧式距离；其中，与/>对应的前个欧式距离中第m个欧式距离记作/>m和M均为正整数，且1≤m≤M；
步骤203、根据公式得到第q期激光点云的平均点间距ρ(q)；
步骤204、按照步骤201至步骤203所述的方法，得到第q+1期的激光点云的平均点间距ρ(q+1)；并对ρ(q)和ρ(q+1)进行平均值处理，得到相邻两期的平均点间距ρ(q,q+1)；
步骤205、根据公式δ(q,q+1)＝B×ρ(q,q+1)，得到相邻两期的点云切片厚度δ(q,q+1)；其中，B为常数且B取值为0.5～1.5。
本实施例中，步骤304对第q期第j个点云切片上第p个投影点集合进行处理，得到余切权拉普拉斯矩阵L，具体过程如下：
步骤3041、采用计算机将第p个投影点记作第1个点，将第p个单环邻域点集合中点云按照逆时针顺序标记为第2个点，...，第g个点；其中，将第1个点至第g个点中任一个点记作第a个点或者第b个点，且a和b均为正整数且取值为1～g；
步骤3042、当a和b取值位于2～g，且b大于a，a和b为相邻两个点时，采用计算机根据公式ωab＝cotαab，得到余切权拉普拉斯矩阵L的第a行第b列元素值ωab；其中，αab为第1个点和第a个点连线与第1个点和第b个点连线之间的夹角，g和2为相邻两个点；并根据ωba＝ωab，得到余切权拉普拉斯矩阵L的第b行第a列元素值ωba；
当a取值为1，当b取值为3～g-2时，采用计算机根据公式ω1b＝cotα1b+cotβ1b，得到余切权拉普拉斯矩阵L的第1行第b列元素值ω1b；其中，α1b为第1个点和第b-1个点连线与第b-1个点和第b个点连线之间的夹角；βab为第1个点和第b+1个点连线与第b+1个点和第b个点连线之间的夹角；并根据ω1b＝ωb1，得到余切权拉普拉斯矩阵L的第b行第1列元素值ωb1；
当a取值为1，当b取值为2时，采用计算机根据公式ω12＝cotα12+cotβ12，得到余切权拉普拉斯矩阵L的第1行第2列元素值ω12；其中，α12为第1个点和第g个点连线与第g个点和第2个点连线之间的夹角；β12为第1个点和第3个点连线与第3个点和第2个点连线之间的夹角；并根据ω12＝ω21，得到余切权拉普拉斯矩阵L的第2行第1列元素值ω21；
当a取值为1，当b取值为g-1时，采用计算机根据公式ω1(g-1)＝cotα1(g-1)+cotβ1(g-1)，得到余切权拉普拉斯矩阵L的第1行第g-1列元素值ω1(g-1)；其中，α1(g-1)为第1个点和第g-2个点连线与第g-2个点和第g-1个点连线之间的夹角；β1(g-1)为第1个点和第g个点连线与第g个点和第g-1个点连线之间的夹角；并根据ω1(g-1)＝ω(g-1)1，得到余切权拉普拉斯矩阵L的第g-1行第1列元素值ω(g-1)1；
当a取值为1，当b取值为g时，采用计算机根据公式ω1g＝cotα1g+cotβ1g，得到余切权拉普拉斯矩阵L的第1行第g列元素值ω1g；其中，α1g为第1个点和第g-1个点连线与第g-1个点和第g个点连线之间的夹角；β1g为第1个点和第2个点连线与第2个点和第g个点连线之间的夹角；并根据ω1g＝ωg1，得到余切权拉普拉斯矩阵L的第g行第1列元素值ωg1；
当a和b取值为1，采用计算机根据公式
当a和b取值相同且均为2时，采用计算机根据公式ω22＝-(ω21+ω23+ω2g)，得到余切权拉普拉斯矩阵L的第2行第2列元素值ω22；其中，ω21＝ω12，ω23＝ω32；
当a和b取值相同且均为3～g-2时，采用计算机根据公式ωaa＝-(ωa1+ωa(b-1)+ωa(b+1))，得到余切权拉普拉斯矩阵L的第a行第a列元素值ωaa；其中，ωa1＝ω1a；
当a和b取值相同且均为g-1时，采用计算机根据公式ω(g-1)(g-1)＝-(ω(g-1)1+ω(g-1)(b-1)+ω(g-1)(b+1))，得到余切权拉普拉斯矩阵L的第g-1行第g-1列元素值ω(g-1)(g-1)；其中，ω(g-1)1＝ω1(g-1)；
当a和b取值相同且均为g时，采用计算机根据公式ωgg＝-(ωg1+ωg(g-1)+ωg2)；其中，ωg1＝ω1g，ωg2＝ω2g；
步骤3043、将其余元素值设置为零，则得到余切权拉普拉斯矩阵L；其中，L为对称矩阵。
本实施例中，步骤305，具体过程如下：
步骤3051、采用计算机根据公式得到第t+1次第p个投影点集合收缩后的坐标矩阵Pt+1；其中，t为大于等于0的自然数，Wl t为第t次更新的第一对角矩阵，/>为第t次更新的第二对角矩阵，Lt为第t次第p个投影点集合收缩后的余切权拉普拉斯矩阵，Pt为第t次第p个投影点集合收缩后的坐标矩阵，O表示零矩阵；当t取值为0时，Wl 0为初始的第一对角矩阵，且对角线元素值为1；/>为初始的第二对角矩阵，且对角线元素值为S0表示初始的第p个单环邻域点集合围成的邻域面积，P0为初始的第p个投影点集合的坐标矩阵；L0为初始的余切权拉普拉斯矩阵即L；
步骤3052、采用计算机根据公式得到第t+1次第一对角矩阵的更新矩阵Wl t+1；其中，Sl取值为3.0；/>
采用计算机根据公式得到第t+1次第二对角矩阵的更新矩阵/>其中，St表示第t次收缩后第p个单环邻域点集合围成的邻域面积；
步骤3053、多次重复步骤3051和步骤3052，直至到达设定的迭代次数T，则得到第T次第p个投影点结合收缩后的坐标矩阵PT，并通过PT得到第p个投影点集合中各个收缩后的坐标。
本实施例中，步骤404中获取剔除空间分布式不确定性值LoD后第j个有序点云上第e个切平面区块的面积记作，具体过程如下：
步骤4041、采用计算机第j个有序点云上第e个切平面区块的Y轴方向最大值Ymax和Y轴方向最小值Ymin，从Y轴方向最小值开始以步长Δx，沿Y轴正方向依次设置多条直线R；
步骤4042、将第i＇条直线记作直线Ri＇，且Ri＇＝Ymin+n×Δx；其中，n表示直线数量，i＇为正整数；且n＝Ceiling[(Ymax-Ymin)/Δx]；1≤i＇≤n；
步骤4043、采用计算机将第q期第j个有序点云上距离直线Ri＇最近且分布在直线Ri＇两侧的两个点的连线和直线Ri＇的交点记作第一个交点，将第q+1期第j个有序点云上距离直线Ri＇最近且分布在直线Ri＇两侧的两个点的连线和直线Ri＇的交点记作第二个交点，采用计算机将第一个交点和第二个交点之间的距离记作距离Hi＇；
步骤4044、采用计算机利用CloudCompare软件获取距离Hi＇的不确定性值LoD(i＇)；
步骤4045、采用计算机将距离Hi＇和LoD(i＇)进行判断修正，如果距离Hi＇的绝对值小于等于LoD(i＇)，则距离Hi＇的修正值为Hi″＝0；如果距离Hi＇为负值且绝对值大于LoD(i＇)，则距离Hi＇的修正值为Hi″＝Hi＇+LoD(i＇)；如果距离Hi＇为正值且绝对值大于LoD(i＇)，则距离Hi＇的修正值为Hi″＝Hi＇-LoD(i＇)；
步骤4046、采用计算机根据公式S(i＇)＝Hi″×Δx，得到第i＇个直线处面积；
步骤4047、多次重复步骤4042至步骤4046，得到各个直线处面积，并将各个直线处面积求和，得到第j个有序点云上第e个切平面区块的面积
本实施例中，步骤4034中采用计算机利用CloudCompare软件获取距离Hi＇的不确定性值LoD(i＇)，具体过程如下：
步骤A、采用计算机利用CloudCompare软件“M3C2 distance”工具，通过输入第q期激光点云和第q+1期激光点云进行处理，得到带有不确定性误差属性值的M3C2点云；其中，所述M3C2点云中每个点具有不确定性值；
步骤B、采用计算机获取第一个交点和M3C2点云中所有点的欧式距离，并筛选出距离第一个交点小于0.5m的点集D；
步骤C、采用计算机根据公式得到第一个交点的不确定性值LoD，并记作距离Hi＇的不确定性值LoD(i＇)；其中，zk表示点集D中第k个点的不确定性值，dk表示点集D中第k个点距离第一个交点的欧式距离，k为正整数，且1≤k≤K＇，K＇表示点集D中点的总数。
本实施例中，步骤501中根据第j个有序点云上第e个切平面区块的面积，获取第j个有序点云上第e个切平面区块的变化体积，进而得到第j个点云切片的变化体积；
步骤5011、当j处于1～J-1时，则根据公式，得到第j个点云切片上第e个切平面区块的变化体积/>
当j取J时，则根据公式得到第J个点云切片上第e个切平面区块的变化体积/>
步骤5012、当j处于1～J-1时，采用计算机根据得到第j个点云切片的变化体积Vj(q,q+1)；
当j取J时，采用计算机根据得到第J个点云切片的变化体积VJ(q,q+1)；
步骤502中将各个点云切片的变化体积累加得到待测流域区域第q+1期相对第q期的体积变化量V(q,q+1)，具体过程如下：
采用计算机根据得到待测流域区域第q+1期相对第q期的体积变化量V(q,q+1)。
本实施例中，在OXYZ坐标系中也可设置Y轴正方向沿待测流域的主沟道延伸方向，则步骤203和步骤301中X轴更换为Y轴；步骤403中平面OYZ更换为平面OXZ；步骤4041中Y轴更换为X轴；
本实施例中，采用DJI M600无人机搭载SZT-R250三维激光雷达系统沿设定航线进行扫描数据采集；其中，飞行高度高于地形海拔70m，飞行速度为7km/h，扫描条带宽度为60m，扫描角度为90°-270°，扫描脉冲速率为100kHz。
本实施例中，设定航线是利用pix4D软件提供的“DOUBLE GRi＇D For3D Models”功能对航线进行规划，航线整体呈“井”字型。
本实施例中，设定航线之前需要进行现场踏勘。
本实施例中，步骤303的K近邻方法中K的取值为int(0.012P)，且如果int(0.012P)小于8，则K取8，如果int(0.012P)大于30，则K取30。
本实施例中，相邻两期的时间间隔为10days～20days，可以根据实际需要进行调整。
本实施例中，需要说明的是，第q期的激光点云和第q+1期的激光点的切平面在X轴方向的位置相同，以及总的切平面个数J相同。
本实施例中，当j取1时，E1表示第1个点云切片上切平面区块的总数；当j取2时，E2表示第2个点云切片上切平面区块的总数；...；当j取J时，EJ表示第J个点云切片上切平面区块的总数，且E1、E2和EJ均为正整数。
本实施例中，Wl t、、Wl 0和/>的矩阵均为g×g；O的大小为g×2；采用步骤304的方法得到第t次第p个投影点集合收缩后的余切权拉普拉斯矩阵Lt；
本实施例中，设定的迭代次数T取值为5～10，且t小于等于T。
本实施例中，Wl t为控制收缩力度的对角矩阵，为控制保持原有形状力度的对角矩阵。
本实施例中，Ceiling[]表示取整函数。
本实施例中，步长Δx取值为相邻两期的点云切片厚度δ(q,q+1)的1/30～1/10。
本实施例中，步骤201中第q个点云集合的总数I的取值为第q期激光点云的总数的2％～5％；M为正整数，M小于I，且M的取值范围为3～6。
本实施例中，需要说明的是，同步绘制在平面OYZ上的第q期第j个有序点云和第q+1期第j个有序点云/>的最小Y轴点处两点取均值连接闭合，同步绘制在平面OYZ上的第q期第j个有序点云/>和第q+1期第j个有序点云/>的最大Y轴点处两点取均值连接闭合。
本实施例中，步骤401中原点为OXYZ坐标系的原点，该原点为待测流域的设计点，可根据监测要求设置。
本实施例中，步骤C中点集D中第k个点距离第一个交点的欧式距离dk为三维，且第一个交点的X轴坐标为相应切片处的X轴坐标。
综上所述，本发明方法步骤简单，设计合理，获取沿待测流域的主沟道延伸方向各个点云切片，对各个点云切片进行拉普拉斯收缩，并剔除空间分布式LoD后获取切片面积，以提高获取切片体积变化的精度，有效地适应流域地形变化监测。
以上所述，仅是本发明的较佳实施例，并非对本发明作任何限制，凡是根据本发明技术实质对以上实施例所作的任何简单修改、变更以及等效结构变化，均仍属于本发明技术方案的保护范围内。

Claims (7)

1.一种基于无人机激光雷达的大流域地形变化切片方法，其特征在于，该方法包括以下步骤：
步骤一、待测流域的点云获取：
采用无人机搭载三维激光雷达系统对待测流域进行扫描，得到第q期初始激光点云和第q+1期初始激光点云；其中，q为正整数，且1≤q；第q期初始激光点云和第q+1期初始激光点云的点云坐标为CGCS2000坐标系；
步骤二、待测流域的点云的切片：
步骤201、采用计算机建立OXYZ坐标系；其中，OXYZ坐标系中X轴的正方向沿待测流域的主沟道延伸方向，OXYZ坐标系中Y轴垂直X轴，OXYZ坐标系中Z轴垂直X轴和Y轴垂直向上；
步骤202、采用计算机将第q期初始激光点云和第q+1期初始激光点云转换为OXYZ坐标系下的点云坐标，得到第q期激光点云和第q+1期激光点云；
步骤203、当X轴的正方向沿待测流域的主沟道延伸方向，在OXYZ坐标系下，采用计算机按照相邻两期的点云切片厚度δ(q,q+1)，将第q期激光点云沿X轴方向做多个切平面并将多个切平面沿X轴正方向排序依次记作第1个切平面，...，第j个切平面，...，第J个切平面；其中，j和J均为正整数，且1≤j≤J；任一个切平面垂直X轴，第1个切平面和第J个切平面之间相邻两个切平面的间距为δ(q,q+1)，第1个切平面距离第q期激光点云与第q+1期激光点云中最小X轴坐标处的间距为δ(q,q+1)/2，第J个切平面距离第q期激光点云与第q+1期激光点云中最大X轴坐标处的间距记作L′(J,J-1),且L′(J,J-1)不大于δ(q,q+1)；
步骤三、待测流域点云切片的收缩：
步骤301、当X轴的正方向沿待测流域的主沟道延伸方向，当j处于1～J-1时，采用计算机将第j个切平面左侧及右侧X轴方向δ(q,q+1)/2范围内的点云投影至第j个切平面，得到第q期第j个点云切片；
当j取J时，采用计算机将第J个切平面右侧X轴方向的剩余点云及左侧X轴方向δ(q,q+1)/2范围内的点云投影至第J个切平面，得到第q期第J个点云切片；
步骤302、采用计算机获取第q期第j个点云切片上投影点的总数并记作P；
步骤303、采用计算机利用K近邻方法对第q期第j个点云切片上的点云进行处理，得到第q期第j个点云切片上第p个投影点的K邻域Np，再对Np进行狄洛尼三角剖分，得到第q期第j个点云切片上第p个投影点的单环邻域点集合；其中，第q期第j个点云切片上第p个投影点的单环邻域点集合记作第p个单环邻域点集合；p为正整数，且1≤p≤P；
步骤304、采用计算机将第p个投影点和第p个投影点的单环邻域点集合记作第p个投影点集合，并对第q期第j个点云切片上第p个投影点集合进行处理，得到余切权拉普拉斯矩阵L；
步骤305、采用计算机根据余切权拉普拉斯矩阵L，对第p个投影点集合进行收缩，得到第p个投影点集合中各个收缩后的坐标，提取得到第p个投影点收缩后的坐标；
步骤306、多次重复步骤303至步骤305的方法，完成第q期第j个点云切片上各个投影点的收缩，将收缩后的各个投影点组合为第q期第j个理想点云切片；
步骤307、按照步骤301至步骤306的方法，完成第q+1期第j个点云切片上各个投影点的收缩，得到第q+1期第j个理想点云切片；
步骤四、获取后一期点云切片相对前一期点云切片的切片变化面积：
步骤401、采用计算机将第q期第j个理想点云切片中距离最远的两个点中靠近原点的点作为第一个点，另一个作为终点，利用贪心算法以两点距离最近对第q期第j个理想点云切片中各个点进行排序，得到第q期第j个有序点云；
步骤402、按照步骤401的方法，得到第q+1期第j个有序点云；
步骤403、当X轴的正方向沿待测流域的主沟道延伸方向，采用计算机将第q期第j个有序点云和第q+1期第j个有序点云/>同步绘制在平面OYZ上，且/>和/>中相交的相邻两个交点将/>和/>之间围设的区域划分为第1个切平面区块，..，第e个切平面区块，..，第Ej个切平面区块；其中，e和Ej均为正整数，且1≤e≤Ej，Ej表示第j个点云切片上切平面区块的总数；
步骤404、获取剔除空间分布式不确定性值LoD后第j个有序点云上第e个切平面区块的面积记作
步骤五、待测流域的体积变化量获取：
步骤501、根据第j个有序点云上第e个切平面区块的面积，获取第j个有序点云上第e个切平面区块的变化体积，进而得到第j个点云切片的变化体积；
步骤502、将各个点云切片的变化体积累加得到待测流域区域第q+1期相对第q期的体积变化量V(q,q+1)；
步骤304对第q期第j个点云切片上第p个投影点集合进行处理，得到余切权拉普拉斯矩阵L，具体过程如下：
步骤3041、采用计算机将第p个投影点记作第1个点，将第p个单环邻域点集合中点云按照逆时针顺序标记为第2个点，...，第g个点；其中，将第1个点至第g个点中任一个点记作第a个点或者第b个点，且a和b均为正整数且取值为1～g；
步骤3042、当a和b取值位于2～g，且b大于a，a和b为相邻两个点时，采用计算机根据公式ωab＝cotαab，得到余切权拉普拉斯矩阵L的第a行第b列元素值ωab；其中，αab为第1个点和第a个点连线与第1个点和第b个点连线之间的夹角，g和2为相邻两个点；并根据ωba＝ωab，得到余切权拉普拉斯矩阵L的第b行第a列元素值ωba；
当a取值为1，当b取值为3～g-2时，采用计算机根据公式ω1b＝cotα1b+cotβ1b，得到余切权拉普拉斯矩阵L的第1行第b列元素值ω1b；其中，α1b为第1个点和第b-1个点连线与第b-1个点和第b个点连线之间的夹角；βab为第1个点和第b+1个点连线与第b+1个点和第b个点连线之间的夹角；并根据ω1b＝ωb1，得到余切权拉普拉斯矩阵L的第b行第1列元素值ωb1；
当a取值为1，当b取值为2时，采用计算机根据公式ω12＝cotα12+cotβ12，得到余切权拉普拉斯矩阵L的第1行第2列元素值ω12；其中，α12为第1个点和第g个点连线与第g个点和第2个点连线之间的夹角；β12为第1个点和第3个点连线与第3个点和第2个点连线之间的夹角；并根据ω12＝ω21，得到余切权拉普拉斯矩阵L的第2行第1列元素值ω21；
当a取值为1，当b取值为g-1时，采用计算机根据公式ω1(g-1)＝cotα1(g-1)+cotβ1(g-1)，得到余切权拉普拉斯矩阵L的第1行第g-1列元素值ω1(g-1)；其中，α1(g-1)为第1个点和第g-2个点连线与第g-2个点和第g-1个点连线之间的夹角；β1(g-1)为第1个点和第g个点连线与第g个点和第g-1个点连线之间的夹角；并根据ω1(g-1)＝ω(g-1)1，得到余切权拉普拉斯矩阵L的第g-1行第1列元素值ω(g-1)1；
当a取值为1，当b取值为g时，采用计算机根据公式ω1g＝cotα1g+cotβ1g，得到余切权拉普拉斯矩阵L的第1行第g列元素值ω1g；其中，α1g为第1个点和第g-1个点连线与第g-1个点和第g个点连线之间的夹角；β1g为第1个点和第2个点连线与第2个点和第g个点连线之间的夹角；并根据ω1g＝ωg1，得到余切权拉普拉斯矩阵L的第g行第1列元素值ωg1；
当a和b取值为1，采用计算机根据公式
当a和b取值相同且均为2时，采用计算机根据公式ω22＝-(ω21+ω23+ω2g)，得到余切权拉普拉斯矩阵L的第2行第2列元素值ω22；其中，ω21＝ω12，ω23＝ω32；
当a和b取值相同且均为3～g-2时，采用计算机根据公式ωaa＝-(ωa1+ωa(b-1)+ωa(b+1))，得到余切权拉普拉斯矩阵L的第a行第a列元素值ωaa；其中，ωa1＝ω1a；
当a和b取值相同且均为g-1时，采用计算机根据公式ω(g-1)(g-1)＝-(ω(g-1)1+ω(g-1)(b-1)+ω(g-1)(b+1))，得到余切权拉普拉斯矩阵L的第g-1行第g-1列元素值ω(g-1)(g-1)；其中，ω(g-1)1＝ω1(g-1)；
当a和b取值相同且均为g时，采用计算机根据公式ωgg＝-(ωg1+ωg(g-1)+ωg2)；其中，ωg1＝ω1g，ωg2＝ω2g；
步骤3043、将其余元素值设置为零，则得到余切权拉普拉斯矩阵L；其中，L为对称矩阵。
2.按照权利要求1所述的一种基于无人机激光雷达的大流域地形变化切片方法，其特征在于：步骤一中采用无人机搭载三维激光雷达系统对待测流域进行扫描，获取第q期第1架次激光点云，...，第q期第f架次的激光点云，...，第q期第F架次的激光点云；其中，f和F均为正整数，1≤f≤F，F≥5；
步骤102、采用计算机利用CloudCompare软件中“Merge”工具，对第q期第1架次激光点云至第q期第F架次的激光点云进行合并，得到合并后的第q期激光点云；
步骤103、采用计算机利用CloudCompare软件中“Segment”工具，按照待测流域实地野外踏勘结果对合并后的第q期激光点云进行裁剪，得到第q期待测流域点云；
步骤104、采用计算机利用MCC点云滤波方法对第q期待测流域点云进行一次滤波，得到一次滤波后第q期激光点云；
步骤105、采用计算机利用TerraSolid软件对一次滤波后第q期激光点云进行二次滤波，得到二次滤波后第q期激光点云，并记作第q期初始激光点云；
步骤106、按照步骤101至步骤105所述的方法，得到第q+1期初始激光点云。
3.按照权利要求1所述的一种基于无人机激光雷达的大流域地形变化切片方法，其特征在于：步骤二中相邻两期的点云切片厚度δ(q,q+1)的获取，具体过程如下：
步骤201、采用计算机从第q期激光点云中随机选择点云组成第q个点云集合；其中，第q个点云集合中第i个点记作I表示第q个点云集合的总数；
步骤202、采用计算机获取第q个点云集合中第i个点与第q期激光点云中各个点的欧式距离，并将各个欧式距离按照从小到大的顺序排序，获取前M个欧式距离；其中，与/>对应的前个欧式距离中第m个欧式距离记作/>m和M均为正整数，且1≤m≤M；
步骤203、根据公式得到第q期激光点云的平均点间距ρ(q)；
步骤204、按照步骤201至步骤203所述的方法，得到第q+1期的激光点云的平均点间距ρ(q+1)；并对ρ(q)和ρ(q+1)进行平均值处理，得到相邻两期的平均点间距ρ(q,q+1)；
步骤205、根据公式δ(q,q+1)＝B×ρ(q,q+1)，得到相邻两期的点云切片厚度δ(q,q+1)；其中，B为常数且B取值为0.5～1.5。
4.按照权利要求1所述的一种基于无人机激光雷达的大流域地形变化切片方法，其特征在于：步骤305，具体过程如下：
步骤3051、采用计算机根据公式得到第t+1次第p个投影点集合收缩后的坐标矩阵Pt+1；其中，t为大于等于0的自然数，/>为第t次更新的第一对角矩阵，为第t次更新的第二对角矩阵，Lt为第t次第p个投影点集合收缩后的余切权拉普拉斯矩阵，Pt为第t次第p个投影点集合收缩后的坐标矩阵，O表示零矩阵；当t取值为0时，Wl 0为初始的第一对角矩阵，且对角线元素值为1；/>为初始的第二对角矩阵，且对角线元素值为S0表示初始的第p个单环邻域点集合围成的邻域面积，P0为初始的第p个投影点集合的坐标矩阵；L0为初始的余切权拉普拉斯矩阵即L；
步骤3052、采用计算机根据公式得到第t+1次第一对角矩阵的更新矩阵其中，Sl取值为3.0；
采用计算机根据公式得到第t+1次第二对角矩阵的更新矩阵/>其中，St表示第t次收缩后第p个单环邻域点集合围成的邻域面积；
步骤3053、多次重复步骤3051和步骤3052，直至到达设定的迭代次数T，则得到第T次第p个投影点结合收缩后的坐标矩阵PT，并通过PT得到第p个投影点集合中各个收缩后的坐标。
5.按照权利要求1所述的一种基于无人机激光雷达的大流域地形变化切片方法，其特征在于：步骤404中获取剔除空间分布式不确定性值LoD后第j个有序点云上第e个切平面区块的面积记作具体过程如下：
步骤4041、采用计算机第j个有序点云上第e个切平面区块的Y轴方向最大值Ymax和Y轴方向最小值Ymin，从Y轴方向最小值开始以步长Δx，沿Y轴正方向依次设置多条直线R；
步骤4042、将第i＇条直线记作直线Ri＇，且Ri＇＝Ymin+n×Δx；其中，n表示直线数量，i＇为正整数；且n＝Ceiling[(Ymax-Ymin)/Δx]；1≤i＇≤n；
步骤4043、采用计算机将第q期第j个有序点云上距离直线Ri＇最近且分布在直线Ri＇两侧的两个点的连线和直线Ri＇的交点记作第一个交点，将第q+1期第j个有序点云上距离直线Ri＇最近且分布在直线Ri＇两侧的两个点的连线和直线Ri＇的交点记作第二个交点，采用计算机将第一个交点和第二个交点之间的距离记作距离Hi＇；
步骤4044、采用计算机利用CloudCompare软件获取距离Hi＇的不确定性值LoD(i＇)；
步骤4045、采用计算机将距离Hi＇和LoD(i＇)进行判断修正，如果距离Hi＇的绝对值小于等于LoD(i＇)，则距离Hi＇的修正值为Hi＇＇＝0；如果距离Hi＇为负值且绝对值大于LoD(i＇)，则距离Hi＇的修正值为Hi＇＇＝Hi＇+LoD(i＇)；如果距离Hi＇为正值且绝对值大于LoD(i＇)，则距离Hi＇的修正值为Hi＇＇＝Hi＇-LoD(i＇)；
步骤4046、采用计算机根据公式S(i＇)＝Hi＇＇×Δx，得到第i＇个直线处面积；
步骤4047、多次重复步骤4042至步骤4046，得到各个直线处面积，并将各个直线处面积求和，得到第j个有序点云上第e个切平面区块的面积
6.按照权利要求5所述的一种基于无人机激光雷达的大流域地形变化切片方法，其特征在于：步骤4034中采用计算机利用CloudCompare软件获取距离Hi＇的不确定性值LoD(i＇)，具体过程如下：
步骤A、采用计算机利用CloudCompare软件“M3C2 distance”工具，通过输入第q期激光点云和第q+1期激光点云进行处理，得到带有不确定性误差属性值的M3C2点云；其中，所述M3C2点云中每个点具有不确定性值；
步骤B、采用计算机获取第一个交点和M3C2点云中所有点的欧式距离，并筛选出距离第一个交点小于0.5m的点集D；
步骤C、采用计算机根据公式得到第一个交点的不确定性值LoD，并记作距离Hi＇的不确定性值LoD(i＇)；其中，zk表示点集D中第k个点的不确定性值，dk表示点集D中第k个点距离第一个交点的欧式距离，k为正整数，且1≤k≤K＇，K＇表示点集D中点的总数。
7.按照权利要求5所述的一种基于无人机激光雷达的大流域地形变化切片方法，其特征在于：步骤501中根据第j个有序点云上第e个切平面区块的面积，获取第j个有序点云上第e个切平面区块的变化体积，进而得到第j个点云切片的变化体积；
步骤5011、当j处于1～J-1时，则根据公式得到第j个点云切片上第e个切平面区块的变化体积/>
当j取J时，则根据公式得到第J个点云切片上第e个切平面区块的变化体积/>
步骤5012、当j处于1～J-1时，采用计算机根据得到第j个点云切片的变化体积Vj(q,q+1)；
当j取J时，采用计算机根据得到第J个点云切片的变化体积VJ(q,q+1)；
步骤502中将各个点云切片的变化体积累加得到待测流域区域第q+1期相对第q期的体积变化量V(q,q+1)，具体过程如下：
采用计算机根据得到待测流域区域第q+1期相对第q期的体积变化量V(q,q+1)。
CN202311507747.1A 2023-11-13 2023-11-13 一种基于无人机激光雷达的大流域地形变化切片方法 Active CN117523124B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202311507747.1A CN117523124B (zh) 2023-11-13 2023-11-13 一种基于无人机激光雷达的大流域地形变化切片方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202311507747.1A CN117523124B (zh) 2023-11-13 2023-11-13 一种基于无人机激光雷达的大流域地形变化切片方法

Publications (2)

Publication Number Publication Date
CN117523124A CN117523124A (zh) 2024-02-06
CN117523124B true CN117523124B (zh) 2024-04-26

Family

ID=89750818

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202311507747.1A Active CN117523124B (zh) 2023-11-13 2023-11-13 一种基于无人机激光雷达的大流域地形变化切片方法

Country Status (1)

Country Link
CN (1) CN117523124B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113066162A (zh) * 2021-03-12 2021-07-02 武汉大学 一种用于电磁计算的城市环境快速建模方法
CN114998419A (zh) * 2022-08-02 2022-09-02 西安科技大学 一种基于地形点云的沟壑体积变化切片方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012141235A1 (ja) * 2011-04-13 2012-10-18 株式会社トプコン 三次元点群位置データ処理装置、三次元点群位置データ処理システム、三次元点群位置データ処理方法およびプログラム

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113066162A (zh) * 2021-03-12 2021-07-02 武汉大学 一种用于电磁计算的城市环境快速建模方法
CN114998419A (zh) * 2022-08-02 2022-09-02 西安科技大学 一种基于地形点云的沟壑体积变化切片方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Topology authentication for CAPD models based on Laplacian coordinates;Su, ZY 等;COMPUTERS & GRAPHICS-UK;20130703;全文 *
基于切片点云中心的形变监测;苏芬;余锐;;测绘工程;20160825(第08期);全文 *
局部信息约束的三维重建方法研究;张卫龙;中国博士学位论文全文数据库信息科技辑;20200615;全文 *
机载LiDAR监测黄土高原土壤侵蚀的能力评估;李朋飞 等;测绘学报;20230815;全文 *

Also Published As

Publication number Publication date
CN117523124A (zh) 2024-02-06

Similar Documents

Publication Publication Date Title
Kraus et al. Advanced DTM generation from LIDAR data
CN106597416B (zh) 一种地面GPS辅助的LiDAR数据高程差的误差修正方法
CN107038717A (zh) 一种基于立体栅格自动分析3d点云配准误差的方法
CN107167786A (zh) 卫星激光测高数据辅助提取高程控制点方法
CN110780307B (zh) 基于电瓶车车载式激光点云移动测量系统获取道路横断面的方法
CN103020342A (zh) 一种从地面LiDAR数据中提取建筑物轮廓和角点的方法
CN111598780B (zh) 一种适用于机载LiDAR点云的地形自适应插值滤波方法
CN113920262B (zh) 一种增强边缘取样与改进Unet模型的矿区FVC计算方法及系统
CN109100719B (zh) 基于星载sar影像与光学影像的地形图联合测图方法
CN107792115A (zh) 一种利用三维激光点云自动提取既有线轨顶高程方法
CN110726998B (zh) 一种激光雷达扫描测定矿区采煤塌陷盆地的方法
CN113470090A (zh) 基于sift-shot特征的多固态激光雷达外参标定方法
CN107169301B (zh) 一种分而治之航迹关联方法
CN114998419B (zh) 一种基于地形点云的沟壑体积变化切片方法
CN114877838B (zh) 一种基于车载激光扫描系统的道路几何特征检测方法
CN114689015B (zh) 一种提高光学卫星立体影像dsm高程精度的方法
CN114581619A (zh) 一种基于三维定位和二维建图的煤仓建模方法
CN107393004A (zh) 一种获取输电线走廊中建筑物拆迁量的方法及装置
Baily et al. Comparative assessment of analytical and digital photogrammetric methods in the construction of DEMs of geomorphological forms
CN116758234A (zh) 一种基于多点云数据融合的山地地形建模方法
CN117523124B (zh) 一种基于无人机激光雷达的大流域地形变化切片方法
CN113960532A (zh) 一种基于假想源的二次定位计算的微地震定位方法
CN106018198A (zh) 一种气泡粒径的反演计算方法
CN113763570B (zh) 隧道点云高精度快速自动拼接方法
CN110927765B (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