一种利用LiDAR点云数据提取水体边缘线的方法

摘要

本发明公开了一种利用LiDAR点云数据提取水体边缘线的方法，只利用机载LiDAR散乱点云数据即可提取水体轮廓线的方法，可以应用于机载LiADR点云数据水体边缘检测。本发明方法简单易行，执行效率高,不利用其他辅助数据就能够直接从机载Lidar点云中获得水体轮廓线，能够较为精确和全面的提取出水体。

说明书

技术领域

本发明涉及一种提取水体边缘线的方法，特别涉及一种利用LiDAR点云数据提取水体边 缘线的方法。

背景技术

水是决定湿地存亡的必要条件，利用遥感技术获取湿地水体面积对湿地具有重要指示意 义，遥感技术作为一种新型的空间对地观测技术，一直被广泛用于滩涂湿地的动态监测和资 源调查。但传统光学遥感技术在滩涂水体调查方面存在数据质量受天气状况影响严重、夜间 无法工作以及植被覆盖较浓郁区域光学遥感不能穿透植被等影响，因而光学遥感技术在湿地 应用中受到一定限制。

机载激光扫描探测系统是一种主动式对地观测系统，具有高可靠性、高分辨率和高精度 等特征，已在诸多领域得到广泛应用，将机载LiDAR技术应用于湿地潜在储水能力探测方面 具有较大前景，基于散乱点云数据的水体轮廓线精确提取则是获得储水能力的关键步骤。目 前，利用散乱点云数据直接提取出边界点是边缘提取算法的主要手段之一。

发明内容

本发明所要解决的技术问题是提供一种利用LiDAR点云数据提取水体边缘线的方法，只 利用机载LiDAR散乱点云数据即可提取水体轮廓线的方法，可以应用于机载LiADR点云数 据水体边缘检测。

本发明为解决上述技术问题采用以下技术方案：

本发明提供一种利用LiDAR点云数据提取水体边缘线的方法，包括以下具体步骤：

步骤1，读取目标区域的原始机载LiDAR点云数据，得到原始LiDAR点云数据的总个 数n；

步骤2，确定格网划分的步长并确定x方向和y方向格网的个数n_x和n_y；

步骤3，遍历步骤2中得到的n_x×n_y个格网，根据格网内点云的个数来判断格网是否为水 体内网格，并记录每个水体内网格的位置，从而得到所有水体内网格及其位置，具体步骤如 下：

301：令w＝1，v＝1；

302：对第w行第v列格网的四个顶点进行编号，其中该网格的右下角顶点编号为1，其 他顶点按逆时针顺序依次编号为2、3和4，利用公式4计算该格网1～4顶点的坐标

(x_right,y_down)、(x_right,y_up)、(x_left,y_up)和(x_left,y_down)：

$\begin{array}{c}{x}_{left}=\min x+(w-1)\times step\\ x_{right}=\min x+w\times step\\ y_{down}=\min y+(v-1)\times step\\ y_{up}=\min y+v\times step\end{array}---\left(4\right)$

303：计算当前格网内点云的个数z，如果z≥5，则执行步骤304；若z＜5，则将当前格 网标记为水体内格网，转入步骤305；

304：令w＝w+1，如果w≤n_x，则v保持不变，返回步骤302；如果w＞n_x，则令w＝1， v＝v+1，返回步骤302；

305：记录每个水体内网格的位置，从而得到所有水体内网格及其位置；

步骤4，逐个搜索步骤3得到的l个水体内格网的八连通格网，获得目标区域内所有水陆 交界格网，并确定所有水陆交界格网位于水体内的顶点，具体步骤如下：

401：令I＝1；

402；将第I个水体内格网作为种子格网，对种子格网的八连通格网进行编号，其中右下 角格网编号为1，其他格网按逆时针方向依次编号为2～8；

403：分别计算1～8号格网内的点云个数h，若h≥5则将该格网标记为水陆交界格网；

404：令I＝I+1，若I≤l则返回步骤402，否则终止迭代，从而得到目标区域内所有水 陆交界格网，执行步骤405；

405：确定步骤403中标记的水陆交界格网位于水体内的顶点；

步骤5，逐个搜索步骤4得到的K个水陆交界格网，得到每个水陆交界格网的水体边缘点， 具体步骤如下：

501：令t＝1；

502，任选一个第t个水陆交界格网位于水体内的顶点作为种子顶点w₁，采用K-D树组织 种子顶点和该水陆交界格网内的点云，利用最近邻算法获得种子顶点的最近邻点，将最近邻 点作为水体边边缘点的初始种子点o₁；

503：以初始种子点o₁为起点，确定初始种子点o₁逆时针一侧的水体边缘点，具体步骤如 下：

50301：采用k近邻算法搜索初始种子点o₁在该水陆交界格网中的最近邻点云点 p₁,p₂,...,p_q，其中q为搜索到的最近邻点云点的个数；

50302：利用公式5计算初始角α₀：

式中，a为o₁和w₁之间的距离，b为w₁和w₂之间的距离，c为o₁和w₂之间的距离，w₂为 该水陆交界格网中w₁逆时针方向的下一个顶点；

50303：利用公式6计算p_i、w₁和w₂构成的∠p_iw₁w₂的度数，并记为α_i：

式中，i＝1,2,...,q；a_i为p_i和w₁之间的距离，b为w₁和w₂之间的距离，c_i为p_i和w₂之 间的距离；

50304：从α₁,α₂,…,α_q中选择大于初始角α₀的角，并将其对应的最近邻点云点重新记为 N₁,…,N_s，其中s为α₁,α₂,.….,α_q中大于初始角α₀的角的个数；

50305：利用公式7计算∠N_jo₁w₁的度数，并记为β_j：

式中，j＝1,2,…,s；A_j为N_j和o₁之间的距离，a为o₁和w₁之间的距离，C_j为N_j和w₁之 间的距离；

50306：从β₁,β₂,…,β_s中选择度数最小的角，该角对应的最近邻点云点即为水体边缘点；

50307：将步骤50306得到的水体边缘点作为新的初始种子点，重复步骤50301～50306， 继续寻找下一个水体边缘点，直至无新的水体边缘点出现；

504：以初始种子点o₁为起点，确定初始种子点o₁顺时针一侧的水体边缘点，具体步骤如 下：

50401：利用公式8计算初始角θ₀：

式中，a为o₁和w₁之间的距离，e为w₁和w₃之间的距离，f为o₁和w₃之间的距离，w₃ 为该水陆交界格网中w₁顺时针方向的下一个顶点；

50402：利用公式9计算p_i、w₁和w₃构成的∠p_iw₁w₃的度数，并记为θ_i：

式中，a_i为p_i和w₁之间的距离，e为w₁和w₃之间的距离，f_i为p_i和w₃之间的距离；

50403：从θ₁,θ₂,...,θ_q中选择大于初始角θ₀的角，并将其对应的最近邻点云点重新记为 M₁,M₂,...,M_d，其中d为θ₁,θ₂,...,θ_q大于初始角θ₀的角的个数；

50404：利用公式10计算∠M_uo₁w₁的度数，并记为ω_u：

式中，u＝1,2,…,d；D_u为M_u和o₁之间的距离，a为o₁和w₁之间的距离，F_u为M_u和w₁之 间的距离；

50405：从ω₁,ω₂,...,ω_d中选择度数最小的角，该角对应的最近邻点云点即为水体边缘点；

50406：将步骤50305得到的水体边缘点作为新的初始种子点，重复步骤50401～50405， 继续寻找下一个水体边缘点，直至无新的水体边缘点出现；

505：整合步骤503和步骤504获得的水体边缘点，得到当前水陆交界格网的所有水体边 缘点；

506，令t＝t+1，若t≤K，则返回步骤502，搜索下一个水陆交界格网的水体边缘点；若 t＞K，则终止迭代，执行步骤6；从而得到所有水陆交界格网的水体边缘点；

步骤6，整合步骤5中得到的所有水陆交界格网的水体边缘点，从而得到目标区域水体的 边缘图像。

作为本发明的进一步优化方案，步骤2中确定格网划分的步长并确定x方向和y方向格 网的个数，具体步骤如下：

201：预设每个格网期望包含的点云个数m；

202：计算点云x方向的坐标最小值minx和最大值maxx，以及y方向的坐标最小值miny 和最大值maxy；

203：计算目标区域点云的平均密度；

204：计算格网划分的步长；

205：计算x方向和y方向格网的个数，从而得到n_y行n_x列的格网。

作为本发明的进一步优化方案，利用公式1计算步骤203中的目标区域点云的平均密度， 具体为：

$\sigma =\frac{n}{(\max x-\min x)\times (\max y-\min y)}---\left(1\right)$

式中，σ为目标区域点云的平均密度。

作为本发明的进一步优化方案，利用公式2计算步骤204中的格网划分的步长，具体为：

$step=\sqrt{\frac{m}{\sigma }}---\left(2\right)$

式中，step为格网划分的步长。

作为本发明的进一步优化方案，利用公式3计算步骤205中x方向和y方向格网的个数， 具体为：

$\begin{array}{c}{n}_x=fix\left(\frac{\max x-\min x}{step}\right)\\ n_y=fix\left(\frac{\max y-\min y}{step}\right)\end{array}---\left(3\right)$

式中，n_x和n_y分别为x方向和y方向格网的个数，fix表示截取尾数函数。

作为本发明的进一步优化方案，步骤405中水陆交界格网位于水体内的顶点即为：水陆 交界格网与种子格网的共同顶点。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：本发明方法简单易行， 执行效率高,不利用其他辅助数据就能够直接从机载LiDAB点云中获得水体轮廓线，能够较为 精确和全面的提取出水体。

附图说明

图1是本发明的方法流程图。

图2是格网顶点编号方法。

图3是种子格网的八连通格网编号方法。

图4是搜索水体边缘点的原理图。

图5是原始机载LiDAR点云图像。

图6是提取的水体边缘图像。

具体实施方式

下面结合附图以及具体实施例对本发明的技术方案做进一步的详细说明：

如图1所示，一种利用LiDAR点云数据提取水体边缘线的方法，包括以下步骤：

1、某湖泊区域的原始机载LiDAR点云图像如图5所示，读取原始机载LiDAR点云数据， 点云数据的总个数n＝1244620。

2、确定格网划分的步长并确定x方向和y方向格网的个数，具体步骤如下：

201：预设每个格网期望包含的点云个数k＝100；

202：计算X方向最小值minx和最大值maxx，以及Y方向最小值miny和最大值maxy；

204：利用公式1计算目标区域的点云的平均密度σ：

$\sigma =\frac{n}{(\max x-\min x)\times (\max y-\min y)}---\left(1\right)$

205：利用公式2计算格网划分的步长step：

$step=\sqrt{\frac{m}{\sigma }}---\left(2\right)$

206：利用公式3计算x方向和y方向格网的个数，共得到了n_y行n_x列个格网。

$\begin{array}{c}{n}_x=fix\left(\frac{\max x-\min x}{step}\right)=113\\ n_y=fix\left(\frac{\max y-\min y}{step}\right)=109\end{array}---\left(3\right)$

步骤3，遍历109×113个格网，根据格网内点云的个数来判断格网是否为水体内网格， 并记录每个水体内网格的位置，从而得到所有水体内网格及其位置，具体步骤如下：

301：令w＝1，v＝1；

302：对第w行第v列格网的四个顶点进行编号，如图2所示，其中该网格的右下角顶点 编号为1，其他顶点按逆时针顺序依次编号为2、3和4，利用公式4计算该格网1～4顶点的 坐标(x_right,y_down)、(x_right,y_up)、(x_left,y_up)和(x_left,y_down)：

$\begin{array}{c}{x}_{left}=\min x+(w-1)\times step\\ x_{right}=\min x+w\times step\\ y_{down}=\min y+(v-1)\times step\\ y_{up}=\min y+v\times step\end{array}---\left(4\right)$

303：计算当前格网内点云的个数z，如果z≥5，则执行步骤304；若z＜5，则将当前格 网标记为水体内格网，转入步骤305；

304：令w＝w+1，如果w≤109，则v保持不变，返回步骤302；如果w＞109，则令w＝1， v＝v+1，返回步骤302；

305：记录每个水体内网格的位置，从而得到所有水体内网格及其位置。

步骤4，逐个搜索步骤3得到的l个水体内格网的八连通格网，获得目标区域内所有水陆 交界格网，并确定所有水陆交界格网位于水体内的顶点，具体步骤如下：

401：令I＝1；

402；将第I个水体内格网作为种子格网，对种子格网的八连通格网进行编号，如图3所 示，其中右下角格网编号为1，其他格网按逆时针方向依次编号为2～8；

403：分别计算1～8号格网内的点云个数h，若h≥5则将该格网标记为水陆交界格网；

404：令I＝I+1，若I≤l则返回步骤402，否则迭代结束，从而得到目标区域内所有水陆 交界格网，执行步骤405；

405：确定步骤403中标记的水陆交界格网位于水体内的顶点，即为水陆交界格网与种子 格网的共同顶点。

步骤5，逐个搜索步骤4得到K个水陆交界格网，得到每个水陆交界格网的水体边缘点， 如图4所示，具体步骤如下：

501：令t＝1；

502，任选一个第t个水陆交界格网位于水体内的顶点作为种子顶点w₁，采用K-D树组织 种子顶点和该水陆交界格网内的点云，利用最近邻算法获得种子顶点的最近邻点，将最近邻 点作为水体边缘点的初始种子点o₁；

503：以初始种子点o₁为起点，确定初始种子点o₁逆时针一侧的水体边缘点，具体步骤如 下：

50301：采用k近邻算法搜索初始种子点o₁在该水陆交界格网中的最近邻点云点 p₁,p₂,...,p_q，其中q为搜索到的最近邻点云点的个数；

50302：利用公式5计算初始角α₀：

式中，a为o₁和w₁之间的距离，b为w₁和w₂之间的距离，c为o₁和w₂之间的距离，w₂为 该水陆交界格网中w₁逆时针方向的下一个顶点；

50303：利用公式6计算p_i、w₁和w₂构成的∠p_iw₁w₂的度数，并记为α_i：

式中，i＝1,2,…,q；a_i为p_i和w₁之间的距离，b为w₁和w₂之间的距离，c_i为p_i和w₂之 间的距离；

50304：从α₁,α₂,...,α_q中选择大于初始角α₀的角，并将其对应的最近邻点云点重新记为 N₁,...,N_s，其中s为α₁,α₂,…,α_q大于初始角α₀的角的个数；

50305：利用公式7计算∠N_jo₁w₁的度数，并记为β_j：

式中，j＝1,2,…,s；A_j为N_j和o₁之间的距离，a为o₁和w₁之间的距离，C_j为N_j和w₁之 间的距离；

50306：从β₁,β₂,...,β_s中选择度数最小的角，该角对应的最近邻点云点即为水体边缘点；

50307：将步骤50306得到的水体边缘点作为新的初始种子点，重复步骤50301～50306， 继续寻找下一个水体边缘点，直至无新的水体边缘点出现；

504：以初始种子点o₁为起点，确定初始种子点o₁顺时针一侧的水体边缘点，具体步骤如 下：

50401：利用公式5计算初始角θ₀：

式中，a为o₁和w₁之间的距离，e为w₁和w₃之间的距离，f为o₁和w₃之间的距离，w₃为该水陆交界格网中w₁顺时针方向的下一个顶点；

50402：利用公式6计算p_i、w₁和w₃构成的∠p_iw₁w₃的度数，并记为θ_i：

式中，a_i为p_i和w₁之间的距离，e为w₁和w₃之间的距离，f_i为p_i和w₃之间的距离；

50403：从θ₁,θ₂,…,θ_q中选择大于初始角θ₀的角，并将其对应的最近邻点云点重新记为 M₁,M₂,…,M_d，其中d为θ₁,θ₂,…,θ_q大于初始角θ₀的角的个数；

50404：利用公式7计算∠M_uo₁w₁的度数，并记为ω_u：

式中，u＝1,2,…,d；D_u为M_u和o₁之间的距离，a为o₁和w₁之间的距离，F_u为M_u和w₁之 间的距离；

50405：从ω₁,ω₂,...,ω_d中选择度数最小的角，该角对应的最近邻点云点即为水体边缘点；

50406：将步骤50305得到的水体边缘点作为新的初始种子点，重复步骤50401～50405， 继续寻找下一个水体边缘点，直至无新的水体边缘点出现；

505：整合步骤503和步骤504获得的水体边缘点，得到当前水陆交界格网的所有水体边 缘点；

506，令t＝t+1，若t≤K，则返回步骤502，搜索下一个水陆交界格网的水体边缘点；若 t＞K，则终止执行，执行步骤6；从而得到所有水陆交界格网的水体边缘点。

步骤6，整合步骤5中得到的所有水陆交界格网的水体边缘点，从而得到目标区域水体的 边缘图像，如图6所示。

从图5和图6可以看出，在不利用其他辅助数据的情况下，本发明就能够直接从原始机 载LiDAR点云图像中获取到水体的轮廓点云数据，而且提取的点云图像能够较为精确和全面 地体显示水体的轮廓。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟 悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明 的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。