一种基于机载激光雷达的水库库容测算方法与流程

1.本技术涉及水库库容测算技术领域，特别是涉及一种基于机载激光雷达的水库库容测算方法。


背景技术：

2.目前，无人机激光雷达具有结构小巧美观、集成化以及自动化程度高等优点，在水利工程测绘等领域取得了较为理想的效果。水库库容是水利工程建设和运行管理的重要参数，是确定防洪调度、大坝安全监测、流域下游安全以及水利功能是否正常的重要依据。水库库容主要通过获取地下和地面地形数据计算得到，过去地形数据主要通过人工实地测量获取，该方法不仅费时费力，且获取周期较长，无法满足实际数据更新的需求。随着机载激光雷达技术发展，越来越多研究者将这项技术应用于水利领域中。机载激光雷达具有可操作性高、成果精度高、效率高、控制点少、工作周期短且不受作业时间和作业范围限制等特点，该技术在水库库容监测中有不可比拟的优势。
3.水库库容是水库调度的重要参数，其精度直接关系到水库的防洪安全与蓄水兴利。传统的水库库容主要通过人工水深测量结合地形资料或者卫星遥感和测高数据计算所得，获取周期较长，并且随着时间推移，水库本身也会引起局部变化，已有的库容数据无法满足实际数据更新的需求。库容大小决定着水库调节径流的能力和它所能提供的供水、防洪、调蓄等效益。因此，确定水库特征水位及相应库容是水利水电设计的重要任务之一。


技术实现要素：

4.基于此，针对上述技术问题，提供一种能够实现对水库库容的实时高精度快速测算的基于机载激光雷达的水库库容测算方法。
5.一种基于机载激光雷达的水库库容测算方法，所述方法包括：
6.获取待测水库区域，并根据待测水库区域生成多个kml文件；
7.选取无人机机载激光雷达，将所述kml文件导入无人机机载激光雷达遥控器中；
8.根据当天飞行条件和周边地形特征依次在所述无人机机载激光雷达遥控器预设飞行高度、飞行速度、扫描方式、回波强度、回波次数、相机参数、旁向重叠度、航向重叠度，对待测水库区域进行数据采集；
9.对所述待测水库区域进行实地激光点云数据采集；
10.获取无人机采集到的lidar点云数据，生成多个无人机飞行区域点云回波强度图和高程分布图；
11.采用三维正态分布变换算法对无人机采集到的lidar点云数据进行拼接；
12.对拼接后的lidar点云数据进行去噪，并对拼接后的lidar点云数据进行滤波处理，得到真实地面点；
13.提取真实地面点坐标，生成多个arcgis文件，得到水面和水库边界；
14.基于反距离加权插值法对滤波后点云数据缺失部分进行填充或修复，生成多个数
字高程模型；
15.根据得到的数字高程模型，计算该水库坡度，并依次设置不同的参考水位平面和水库坡度，计算水库库容。
16.上述方案中，可选的，所述选取无人机机载激光雷达具体为：选用大疆经纬m300无人机搭载rieglminivux-1uav激光雷达，在飞行前接通飞行器和控制器电源，完成惯导调试、激光雷达预热工作。
17.上述方案中，进一步可选的，所述对待测水库区域进行数据采集包括：所述无人机按照规划航线飞行，对无人机在起飞前和降落后分别进行静态观测，在航飞任务结束之后，统计检查采集数据范围和点云密度。
18.上述方案中，进一步可选的，所述对所述待测水库区域进行实地激光点云数据采集具体为：利用gps测量设备选取水库边界区域进行实地激光点云数据采集；
19.根据坝体、道路、草地不同地形条件，在水库周边按不同地物类型选取区域采集高程。
20.上述方案中，进一步可选的，所述三维正态分布变换算法是利用一系列正态分布概率描述表面信息，将点云分割成格则均匀的体素，并计算体素内点云的概率分布，作为模型的匹配基础；
21.求取参数变换后的点云在对应体素内的概率；
22.利用牛顿迭代法将概率最大时变换参数的最优解作为最终解。
23.上述方案中，进一步可选的，所述对无人机采集到的lidar点云数据进行拼接包括：初始化参数，对于需要拼接的无人机a和b采集到的航带a和航带b，设定航带b看做基准点云；
24.设定6个变换参数模型p6，建立体素单元b，将参考点云y＝{y1,
…
,yn}中所有点yk∈y分配到对应的体素单元bi∈b；其中，所述6个变换参数包括3个空间平移参数和3个空间旋转参数，分别用t
x
、ty、tz和φ
x
、φy、φz表示，对于空间中任意一点x经过“六参数”模型p6＝[t
x t
y t
z φ
x φ
y φz]
t
变换得到te(p6，x)；
[0025][0026][0027][0028][0029]
式中ci＝cosφi,si＝sin φi；
[0030]
在所有体素单元bi∈b利用以下公式式求取其中的正态分布函数：
[0031][0032][0033][0034]
式中，y
k＝1,
…
,n
为体素内点所对应的维度值，μ为所有体素内点的向量均值，c为所有体素内点云的协方差矩阵，是调整概率分布的和为1，用一个常量c0来代替；
[0035]
将航带a中点云x＝{x1,
…
,xn}中所有的点经过模型p6变换，得到te(p6,x)；
[0036]
将变换后航带a的点云与航带b的点云y＝{y1,
…
,yn}对应，即求te(p6,x)与y＝{y1,
…
,yn}最大相似时p6，即
[0037][0038]
利用牛顿迭代求取最优解：hδp＝-g式中，h和g分别是黑塞矩阵和梯度向量，δp是变换模型每次迭代的增量；
[0039][0040][0041]
根据得到的δp6，并计算点云的p6，其中p6＝p6+δp6，完成所有点云的拼接。
[0042]
上述方案中，进一步可选的，所述对拼接后的lidar点云数据进行去噪具体为：采用k-d树自上而下的递归方式作空间划分及近邻搜索；设标准差倍数，对每一个点搜索指定邻域点个数的相邻点，计算点到相邻点的平均值距离，若距离平均值大于最大距离，则判定为噪点并去除。
[0043]
上述方案中，进一步可选的，所述水库坡度计算之前，还包括：对lidar点云赋色；具体为：采用正射赋色法对lidar点云与高分辨率影像匹配，将影像中颜色信息赋予lidar点云，获得真彩色lidar点云。
[0044]
上述方案中，进一步可选的，所述计算水库库容，具体方法为：
[0045][0046]
其中，v为库容变化量，ka为面积系数，α为系数，bc为特征宽度，i为坡度，h0为坝前水深。
[0047]
本发明至少具有以下有益效果：
[0048]
本发明基于对现有技术问题的进一步分析和研究，通过获取待测水库区域，并根据待测水库区域生成多个kml文件，选取无人机机载激光雷达，将所述kml文件导入无人机
minivux-1uav激光雷达，飞行前接通飞行器和控制器电源，完成惯导调试、激光雷达预热等工作。
[0060]
根据当天飞行条件和周边地形特征依次在所述无人机机载激光雷达遥控器预设飞行高度、飞行速度、扫描方式、回波强度、回波次数、相机参数、旁向重叠度、航向重叠度，对待测水库区域进行数据采集。
[0061]
对所述待测水库区域进行实地激光点云数据采集，其中，对所述待测水库区域进行实地激光点云数据采集可以为，利用gps测量设备选取水库边界区域进行实地激光点云数据采集。根据坝体、道路、草地等不同地形条件，在水库周边按不同地物类型选取区域采集高程。
[0062]
获取无人机采集到的lidar点云数据，生成多个无人机飞行区域点云回波强度图和高程分布图。
[0063]
采用三维正态分布变换算法对无人机采集到的lidar点云数据进行拼接，对拼接后的lidar点云数据进行去噪，并对拼接后的lidar点云数据进行滤波处理，得到真实地面点。
[0064]
其中，采用三维正态分布变换(3d-ndt)对多架次点云数据集进行拼接。3d-ndt算法是利用一系列正态分布概率描述表面信息，首先将点云分割成格则均匀的体素，并计算体素内点云的概率分布，作为模型的匹配基础；然后求取参数变换后的点云在对应体素内的概率；最后利用牛顿迭代将概率最大时变换参数的最优解作为最终解。由于机载lidar点云存在大量离散的粗差点，且分布无规律，本研究采用k-d树自上而下的递归方式作空间划分及近邻搜索。设标准差倍数，对每一个点搜索指定邻域点个数的相邻点，计算点到相邻点的平均值距离，如果距离平均值大于最大距离，则判定为噪点并去除。
[0065]
三维正态分布变换算法是利用一系列正态分布概率描述表面信息，将点云分割成格则均匀的体素，并计算体素内点云的概率分布，作为模型的匹配基础；
[0066]
求取参数变换后的点云在对应体素内的概率；
[0067]
利用牛顿迭代法将概率最大时变换参数的最优解作为最终解。
[0068]
上述方案中，进一步可选的，所述对无人机采集到的lidar点云数据进行拼接包括：初始化参数，对于需要拼接的无人机a和b采集到的航带a和航带b，设定航带b看做基准点云；
[0069]
设定6个变换参数模型p6，建立体素单元b，将参考点云y＝{y1,
…
,yn}中所有点yk∈y分配到对应的体素单元bi∈b；其中，所述6个变换参数包括3个空间平移参数和3个空间旋转参数，分别用t
x
、ty、tz和φ
x
、φy、φz表示，对于空间中任意一点x经过“六参数”模型p6＝[t
x t
y t
z φ
x φ
y φz]
t
变换得到te(p6，x)；
[0070][0071]
[0072][0073][0074]
式中ci＝cos φi,si＝sin φi；
[0075]
在所有体素单元bi∈b利用以下公式式求取其中的正态分布函数：
[0076][0077][0078][0079]
式中，y
k＝1,
…
,n
为体素内点所对应的维度值，μ为所有体素内点的向量均值，c为所有体素内点云的协方差矩阵，是调整概率分布的和为1，用一个常量c0来代替；
[0080]
将航带a中点云x＝{x1,
…
,xn}中所有的点经过模型p6变换，得到te(p6,x)；
[0081]
将变换后航带a的点云与航带b的点云y＝{y1,
…
,yn}对应，即求te(p6,x)与y＝{y1,
…
,yn}最大相似时p6，即
[0082][0083]
利用牛顿迭代求取最优解：hδp＝-g式中，h和g分别是黑塞矩阵和梯度向量，δp是变换模型每次迭代的增量；
[0084][0085][0086]
根据得到的δp6，并计算点云的p6，其中p6＝p6+δp6，完成所有点云的拼接。
[0087]
提取真实地面点坐标，生成多个arcgis文件，得到水面和水库边界；基于反距离加权插值法对滤波后点云数据缺失部分进行填充或修复，生成多个数字高程模型。
[0088]
其中，滤波后地面点云会有部分空白，在生成多个数字高程模型之前，需要根据临近点云数据的高程进行插值。本研究基于反距离加权插值法对滤波后点云数据缺失部分进行填充和修复，从而生成多个高精度的数字高程模型。本研究选择渐进加密不规则三角网(ptind)滤波算法对去噪后数据进行滤波处理，滤波处理后的结果为真实地面点。具体的滤波过程为，读取剔除噪声点的点云数据，获取地面种子点，然后对这些地面种子点建立不规
则三角网(tin)并进行迭代加密；在每次迭代过程中，都对其余各点到所在三角形的反复角和反复距离进行阈值判断，将满足条件的点加入tin。迭代过程一直持续到没有点可以再加入tin为止；分离出分出地面点云和非地面点云数据，同时输出地面点云数据。
[0089]
得到滤波后的地面点云后，对点云赋色，lidar点云数据具备高精度坐标信息，但缺乏颜色信息，难以用真实的色彩对现实世界进行直观展示。本研究采用正射赋色法对lidar点云与高分辨率影像匹配，将影像中颜色信息赋予lidar点云，从而获得真彩色lidar点云。具体的，将lidar点云进行二维投影转换，获得二维点云数据；将高分辨率遥感影像中像素与其对应的lidar点云进行匹配，把像素中的颜色信息添加到对应的lidar点云中，获得真彩色lidar点云数据。
[0090]
根据得到的数字高程模型，计算该水库坡度，并依次设置不同的参考水位平面和水库坡度，计算水库库容。其中，水库坡度计算是利用arcgis软件中3d analyst tools-raster surface-slope功能，计算水库坡度。
[0091]
其中，计算水库库容，具体方法为：
[0092][0093]
其中，v为库容变化量，ka为面积系数，α为系数，bc为特征宽度，i为坡度，h0为坝前水深。
[0094]
点云滤波是指在激光扫描的点云数据中，提取出精确的地面点，进而得到dem的过程。经过滤波后的能将lidar点云数据中的地面点和非地面点进行分离。经过滤波处理后的点云可以将飞行时生物和大气折射的影响去除，获取地面种子点，然后对这些地面种子点建立不规则三角网(tin)并进行迭代加密。在每次迭代过程中，都对其余各点到所在三角形的反复角和反复距离进行阈值判断，将满足条件的点加入tin，滤波效果图1所示。
[0095]
在机载lidar点云数据处理中，由于地面植被以及建筑的影响，激光难以穿透障碍物到达地表，导致有的点位于真实地面，有的点则位于植被或者建筑物上。由于机载lidar点云存在大量离散的粗差点，且分布无规律，本研究采用k-d树自上而下的递归方式作空间划分及近邻搜索。设标准差倍数，对每一个点搜索指定邻域点个数的相邻点，计算点到相邻点的平均值距离，如果距离平均值大于最大距离，则判定为噪点并去除，去噪效果如图2所示。经过去噪处理可将水面反射的大部分异常点清除，并保留完整清晰的水边线。
[0096]
对比lidar测量得出的水库面积与前人实测结果(图3)，lidar测量的水域面积范围为2.12-3.9km2，前人测量的水域面积范围为2-4.1km2。在52米以后，基于lidar测量的水域面积小于1950s水库设计值，差值范围在0.01-0.46km2，主要原因为长期的淤泥堆积和人类活动。图形结果趋势一致，在浮动误差允许范围内，由此可得本实验的研究方法可满足适用性分析需求，可作为水库监测水体库容变化量的参考。
[0097]
本实施例中，通过获取待测水库区域，并根据待测水库区域生成多个kml文件，选取无人机机载激光雷达，将所述kml文件导入无人机机载激光雷达遥控器中，根据当天飞行条件和周边地形特征依次在所述无人机机载激光雷达遥控器预设飞行高度、飞行速度、扫描方式、回波强度、回波次数、相机参数、旁向重叠度、航向重叠度，对待测水库区域进行数据采集，对所述待测水库区域进行实地激光点云数据采集，获取无人机采集到的lidar点云数据，生成多个无人机飞行区域点云回波强度图和高程分布图，采用三维正态分布变换算
法对无人机采集到的lidar点云数据进行拼接，对拼接后的lidar点云数据进行去噪，并对拼接后的lidar点云数据进行滤波处理，得到真实地面点，提取真实地面点坐标，生成多个arcgis文件，得到水面和水库边界，基于反距离加权插值法对滤波后点云数据缺失部分进行填充或修复，生成多个数字高程模型，根据得到的数字高程模型，计算该水库坡度，并依次设置不同的参考水位平面和水库坡度，计算水库库容，本发明采用无人机机载激光雷达该技术具有可操作性高、成果精度高、效率高、控制点少、工作周期短且不受作业时间和作业范围限制等特点。并利用地理信息技术对采集的点云进行后处理，快速更新水库现状库容，并拟合库容曲线，实现对各特征水位下库容的快速估算。
[0098]
应该理解的是，虽然图1的流程图中的步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图1中的至少一部分步骤可以包括步骤或者阶段，这些步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤中的步骤或者阶段的至少一部分轮流或者交替地执行。
[0099]
以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
[0100]
以上所述实施例仅表达了本技术的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本技术构思的前提下，还可以做出若干变形和改进，这些都属于本技术的保护范围。因此，本技术专利的保护范围应以所附权利要求为准。