使用地面控制点检核ALS点云精度

任务描述

​ 使用地面布设的外业控制点GCP来检核ALS点云条带平差之后的绝对位置精度。如下图所示，红圈中的是GCP点，紫色和灰色是两个条带的ALS点云。图中显示了在\(z\)方向上的偏差。

​ 通过观察，有一部分控制点是在平面区域采集的，有一些是在边缘比较锋利的区域采集的。

目前存在的问题

​ 本身ALS点云就比较稀疏，再加上GCP也是离散的布设的，导致很难找到对应同名关系，想要找到对应位置的值就只能通过内插来实现。

​ 内插方法有很多，但是直接在ALS里面内插不是很现实，可以考虑转成格网，mesh的形式，然后基于新的数据表达方式再内插。例如可以用arcgis生成DTM，DEM，内插到高分辨率，然后再采样点和GCP做计算。

解决方案

最终论证，通过地面控制点应该是无法做到水平精度检核的，只能做到高程方向的检核。

相关方法

A Quantitative Assessment of LIDAR Data Accuracy¹

​ 文章里面的高程和平面精度是分开计算的。

• 对于高程精度 首先只能选择平面上的控制点。GCP是有规律线性分布的方式布设的，因此作者就选择1m范围内的点作为领域点，然后拟合平面，计算高程方向上的偏差。

• 对于平面精度

  由于对应关系不好确定，采用两个曲线先平移后相交的方式确定对应关系。首先找到地面控制点一定半径内的点作为候选对应同名点，然后绘制曲线，通过在东(北)方向以0.01\(m\)平移曲线，找到和GCP最吻合的位置（好像是逐位置计算误差，找到误差最小的就以对应时刻的步长作为偏移量），就说明找到了ALS和GCP之间的同名关系。曲线绘制方式是，在东(北)方向和高程方向构成的平面中分别投影GCP和ALS点云，消除高度偏差之后，就可以得到两条大致重叠的曲线。

• 问题：

  依然确实存在平面点不对应的问题，导致精度计算出现问题。

  此外，要求GCP是线性分布的，不然没法计算曲线。对检核点的布设要求很高。

Accuracy of ground surface interpolation from airborne laser scanning (ALS) data in dense forest cover²

简介

​ DTM的精度取决于激光系统测量值的精度和准确度、点云滤波的精度（因为需要滤除地面点才能做DTM），以及需要内插来解决空洞问题。其中内插精度为影响DTM精度第二主要因素，因此文章主要是想在森林密集的区域，通过比较生成的DTM的精度来进一步验证9种内插方式的精度。涉及到的内插方式有：==反距离加权、最近邻、自然内插、Delauney三角网、多层-B-Spline、三次样条、Thin-Plate Spline、Thin-Plate Spline by TIN、克里金插值法==。对比方法为使用cross-validation的方式，将预测值与真值GT做差，计算偏差和RMSE。

结论

​ 内插方法的精度和想要生成的DTM的分辨率有关，DTM分辨率越低，内插精度越差，但是不同内插方法之间的RMSE在点云密度、熵密度、地面坡度等因素下变化不明显。NN和Spline类的方法总体精度相近，但是NN的方法不需要额外的输入，而且效率要比样条线的高一些。其中：

• 反距离权重的精度最低。同时NeN、克里金插值方法受到点密度和陡峭坡度地形影响大，精度会降低。

• 样条线类的方法得到的表面更加平滑，也接近真值，不同函数之间的RMSE均比较小，适用于崎岖的地形=坡度大、陡峭的地方。不同的样条基函数之间没有明显差异。反倒是越复杂计算时间越长。

• 只有克里金插值方法受熵密度的影响大，其他方法由于测试数据的问题，并不能得出很充分的相关关系。

• NN，DT，Spline-based 方法受到地面点密度、地面坡度变化的影响小。

Scattered data interpolation with multilevel B-splines⁴

Motivation

​ 目的是从离散点云中实现点云表面拟合B-spline approximation (BA)。现有方法也有用B-spline的，但是依然受限于平滑性，时间复杂性，限定范围内的数据分布等的限制。因此本文提出使用层级的bicubic B-spline修正方法来逼近直接拟合B-spline，并且实现了最小的内存占用。本文是从Image morphing中借鉴过来的。

Basic Idea

​ Let \(\Phi_{ij}\) be the value of the \(ij\)-th control point on lattice \(\Phi \in \Omega ^{(m+3) \times (n+3)}, \Omega = \{ (x,y)| 0 \leq x <m, 0 \leq y < n \}\), located at \((i,j)\) for $ i =-1,0,...,m+1$ and \(j =-1,0,...,m+1\). The ==approximation function== \(f\) is defined in terms of these control points by: \[ f(x,y)=\Sigma_{k=0}^3\Sigma_{l=0}^3B_k(s)B_l(t)\sigma_{(i+k)(j+1)} \] where \((x,y)\) is a scattered point in \(\Omega\), \(i=\lfloor x \rfloor -1,i=\lfloor y \rfloor-1,s= x-\lfloor x \rfloor\), and \(t=y-\lfloor y \rfloor,t \in [0,1)\). \(B_k,B_l\)是均匀三次B-spline函数. 即 \(f: \Omega \mapsto \mathbb{R}\)表明了函数\(f\)将平面domain的点映射为一个表面。

With this formulation, the problem of deriving function \(f\) is reduced to solving for the control points in \(\Phi\) that best approximate the scattered data in \(P=\{(x_c,y_c,z_c)\}\).

Methodology

​ 本文通过多层B-spline来实现表面近似任务中平滑性和准确性之间的平衡。使用控制顶点的层级性得到一系列的函数，通过求和实现近似。coarse的顶点提供粗略的近似，使用finer的顶点构成的函数来优化精度。

​ 那么使用B-spline来做数据拟合最终调整的是什么呢？根据B-pline的定义和另一篇博客可知。近似的目的就是要不断调整样条曲线中控制点的权值，使得能够最大限度的拟合所给的离散数据。理解这个很重要。控制点的位置是固定分布在domain的，也就是随着knot的。表面拟合的结果和控制点的分辨率是有关系的。

• 分辨率越低，越能覆盖到更多的数据点，大量的数据点会得到比较平滑的拟合结果。

• 分别率越低，因为指用到了局部的信息进行拟合，因此局部拟合的更好，但是全局来看不同区域之间的过渡不是也很好，有很多零的区域。

​ 作者提出首先在\(\Omega\)域内构建层级的控制点\(\Phi_0,\Phi_1,...,\Phi_h\)。相邻层间控制点数量比为1:2，首先在coarest的控制点\(\Phi_0\)上运行一次BA算法作为初始平滑近似\(f_0\)。这样会在一些散点上有较大的偏差，例如在点\(P_1\{(x_1,y_1,z_1)\}\)上存在较大的偏差\(\Delta^1 z_c\)(上标1表示是第一次近似得到的偏差)，用\(\Delta^1 z_c\)替换\(z_1\)，变成\(P_1\{(x_1,y_1,\Delta^1 z_c)\}\)。依次类推，用\(f_0+f_1\)可以得到一个更小的偏差\(\Delta^2 z_c=Z_c-f_0(x_c,y_c)-f_1(x_c,y_c)\)。那么从coarest的\(\Phi_0\)开始直到finer的控制点\(\Phi_h\)，最终的近似结果可以表示为\(f=\Sigma_{k=0}^hf_k\)，即一系列层级近似结果的叠加。只需要在全部的点上做一次分辨率低的近似就可以得到全局的形状。按照上述方案，就实现渐进式的平滑近似。时间复杂度是\(\mathcal{O}=hp+\frac{4}{3}mn\)，空间复杂度是\(\mathcal{O}=p+\frac{4}{3}mn\)，并且仍然可以实现\(C^2\)的平滑性。\(mn\)是domain的范围，是正整数。\(p\)是3D点。

​ 由上可知，需要不断的计算B-Spline函数，导致计算效率下降，作者因此提出一种局部B-Spline精化方法实现计算效率的提升。此外作者还根据数据分布提出一种自适应的分层方法。这里就先略了。重点关注前面。

相关实现

• SINTEF，实现了文献 ⁴ 中的多尺度样条曲线表面拟合，但是代码比较老了，需要额外编译，除了代码框架复杂一些，编译不是问题。
• ddemidov: 基于C++11做了一个文献 ⁴ 新的版本，实现了adaptive 的MBA，开发了python的API，并且是header-only的方式，方便使用。根据官方提供的jupyter notebook的测试显示，两个库对于相同数据的内插结果之间并没有很大的误差。

Smooth Approximation and Rendering of Large Scattered Data Sets⁵

Motivation

​ 表面拟合任务要求：近似程度要高、可视质量要好、适用于大数据量、高效、数值计算稳定、自适应性好以及算法简单易实现等。基于spline的方法有很多，e.g. NURBS要求数据分布均匀，存在依赖矩形的分布空间，近似误差大等问题。其他spline的方法，依赖三角网，而且对数据分布要求较高，还需要全局最小二乘近似等。Related Work写的不错。

Methodology

​ 本文提出一种两阶段方法实现表面美观、属性好的表面拟合方法。Step1：仅使用==局部==少量的数据点计算离散的最小二乘多边形段discrete least squares polynomial pieces，基于==SVD==构建了\(Bernstein-B \'{e} zier\)形式以根据实际点云分布空间初始多边形的次数。Step2：基于\(C^1\)连续条件结合相邻的三角化\(Bezier\)块直接得到剩余多边形段。算法不需要全局计算，不依赖三角网，不需要计算导数。

​ 基于\(C^1\)连续的bivariate cubic spline实现了大场景下散点的最优近似拟合。使用局部多边形最小二乘近似的方式来减少数据的局部差异和非均匀分布的问题。但是考虑到实际数据分布，算法使用的cubic spline的次数是自适应变化的。算法的复杂度是和点数成线性关系的。这里去掉了从数据生成三角网的过程，而算法在介绍的时候依赖三角网的生成。中间使用了LOD加速计算。

​ 最后使用了\(Bernstein-B \'{e} zier\)来评估和渲染spline。通过计算spline点和散点之间的单向\(Housdorff\)距离来评估质量。

结论

• 点云越密，构建的多边形阶数越高。拟合误差最大为\(m\)级。
• 在一块a Vpro/8 graphics board可以做到10FPS。

相关实现

• [SAGA][]：SAGA - System for Automated Geoscientific Analyses - is a Geographic Information System (GIS) software with immense capabilities for geodata processing and analysis. SAGA is programmed in the object oriented C++ language and supports the implementation of new functions with a very effective Application Programming Interface (API). Functions are organised as modules in framework independent Module Libraries and can be accessed via SAGA’s Graphical User Interface (GUI) or various scripting environments (shell scripts, Python, R, ...). 提供了编译好的可执行程序和源代码。针对内插方法，SAGA已经实现了B-spline Approximation、Cubic Spline Approximation、Multilevel B-Spline、Multilevel B-Spline (3D)、Multilevel B-Spline for Categories、Multilevel B-Spline form Grid Points、Thin Plate Spline、Thin Plate Spline (TIN)。包含了文章[^ 4]中提到的大部分算法。每种内插方法配有注释。简单测试之后发现确实和文章中提到的结论一致，起码效率方面表现相同。

Approximate thin plate spline mappings⁶

​ Thin plate spline是一种常见的使用基函数实现从\(\mathbb{R}^2\)到\(\mathbb{R}^2\)坐标映射的方法。本身TPS是cubic spline的2D泛化形式，TPS模型包含一个仿射模型作为一个特例。TPS模型需要对一个超大的密集矩阵\(p \times p\)进行求逆，\(p\)是点数。为了减少求逆的运算量，提出了很多的改进方法，包括使用局部对应点集来估计全部对应关系的方法，本文结合高斯径向基网络领域中的相关方法来讨论基于这个思路的近似方法的优缺点。指出借鉴径向基函数的方法不满足principal warp分析⁷的要求，因此借助经典矩阵补全技术进行完善。最后，作者讨论了一种使用全部近似基函数来近似目标集数据的方法，即\(Nystor\ddot{o}m\)近似法实现TPS映射。最后作者并总结了三种近似方法之间的优缺点：基于局部点采样并考虑全部目标集的方法的效果都不错，基于\(Nystor\ddot{o}m\)近似的方法不但可以实现principal warp 分析，同时在子集采样结果不好时表现出最优的性能。并用实验突出了在图像编辑方面的优势。

• 数学有点啃不动。详见文章6和8。

LR B-Spline

​ 最新版书籍⁸中全面细致的总结了不同地形表面表达方式之间的区别，并肯定了local refine B-Spline表达形式（其实是一类方法）的好处。

相关实现

SINTEF下属的应用数学研究所开发的GoTools，The core module contains generic tools and spline functionality. The additional modules contain functionality for intersections, approximative implicitization, parametrization, topology, and more. 包含NURBS 和TIN模板库。

References