激光雷达技术与地学应用（双语）

本文链接：https://blog.youkuaiyun.com/uwpk123/article/details/139922464

Ch 1 Introduction

LIDAR = Light Detection And Ranging
实现空间三维坐标的同步、快速和精确获取

LiDAR technology can realize the synchronization, fast and accurate acquisition of spatial three-dimensional coordinates, and reproduce the real-time and real morphological characteristics of targeting objects, providing a simple and effective means for quick acquisition of spatial information.
激光雷达技术可以实现空间三维坐标的同步、快速和精确获取，并重现目标对象的实时和真实形态特征，为快速获取空间信息提供了简单有效的手段。

LIDAR 的优点
- 非接触 Non-contact
- 高精度 High Precision，毫米/亚毫米级
- 快速 Fast Speed，高达120万个点每秒
- 高密度 High Density，点间距高达毫米级
- 灵活数据采集 Flexible data acquisition，需要低光照和低温环境

1.1 测量原理

测量原理：光线传输时间（Light transit time）
脉冲测距（Pulse ranging）：
- 直接测量时间
- 激光雷达发射的光脉冲通过物体表面反射回激光雷达，测量接收到的光脉冲与发射的脉冲之间的时间差 $\frac{1}{2}\Delta t \cdot{c}$ ，D 代表水平距离
- 时间测量的精度决定测距精度（D 越大精度越差）
- 测距长度远（高达 6 km）；精度低
相位测距（Phase ranging）：
- 间接测量时间
- 通过确定“调制光波”的往返相位差来测量距离， $\Delta N$ 是相位角的比值，f 是频率
  $D=\frac{1}{2}c\frac{\Delta N}{f}\\ \Delta N=\frac{\Delta\varphi}{2\pi}$
- 局限性：测量范围接近100m；不能确定整周数；不能测超过 $2\pi$ 以上的相位差，测距范围短；相位差的分辨率限制测距的精度
- 解决：运用多种频率的波：长尺（386.5m）+短尺（精度尺）（6.57m）
- 精度高，采样率高
三角测量
$Z\tan\alpha\quad Z=\frac{Bf_{0}}{p+f_{0}\tan\alpha}$

1.2 扫描模式

窗口扫描 Window Scanning：安全性高
全景扫描 Panoramic Scanning：320-330 度，安全性相对差
移动扫描 Mobile Scanning：一个自由度变化，随平台动；装载 gnss 以及惯导定位和定姿设备，将局部坐标系归算到统一坐标系
结构光扫描 Structured Light Scanning：一次扫描，每条扫描线包括 640 激光点，一次发射 480 条激光线；一次投影 $640\times480$ 激光点阵
机载 LiDAR 系统
1. 摆动镜扫描 Oscillating mirror scan
2. 旋转棱镜扫描 Rotating polygonal mirror scan
3. 光纤维电扫描 Fibre scan：电距最规则
4. 圆锥镜扫描 Palmer mirror scan：有大量重复点，覆盖无空白，质量更高

1.3 LiDAR 分类

星载激光雷达 Satellite-based Laser Altimeter(激光测高仪)
1. Distance H from the center of the instrument to the surface，从仪器中心到表面的距离H（基于飞行时间）
2. The location of the instrument center o，仪表中心o的位置（基于精确轨道测定(POD)）
3. Laser Beam Attitude Information u，激光束的姿态信息u（基于精确姿态测定(PAD)）
4. Determine the plane coordinates and geodetic height of the laser point，确定激光点的平面坐标和大地测量高度（基于地球中心）
5. 星载系统：ICESat(-2)（2003.1.12 发射第一个），GF-7
6. 优点：a high orbit高轨道、wide field of view宽视场，为海外地区三维控制点和 DEM 的获取提供新途径；具有观测整个天体的能力；在植被垂直分布测量、海面高度测量、云层和气溶胶垂直分布测量以及气候现象检测方面发挥作用
机载 Airborne-LiDAR：laser ranging,Global Positioning System (GPS) and Inertial Navigation System (INS) 激光测距，GPS，惯导系统
- 基于 GNSS 定位大地坐标，基于 INS 确定姿态
- 脉冲柱：返回的点不止一个，在范围内都返回，不同的峰值对应不同的高度
- 得到树的垂直结构信息：反演生物量
- 特点
  - 直接主动测量，受天气影响小（雨点大噪音大）
  - 减少地面工作
  - 运行安全，可在危险区域进行测绘工作
  - 处理速度快，easy to update and high currency（更新快，时效性高）
  - 可以穿过植被
  - 应用于检测输电线
  - 不局限于地表，可以测量地面和地下
- 缺点
  1. 昂贵
  2. 水域边界不好检测
  3. 数据分布不规则，坐标分布不连续
- DEM 的生产
- 主要应用：生成 DTM，三维城市建模，等
移动 Mobile-LiDAR（车载）
- 组成部分：GPS 和INS，数字成像传感器，激光扫描器
- 三维激光扫描技术的应用：
  - 生成输出地形剖面图和轮廓图；
  - 断裂调查和辅助地质记录；
  - 隧道施工验收和变形检测；
  - 斜坡变形监测的应用；
  - 确定危岩体的空间分布位置和边界范围（determination of spatial distribution location and boundary range of dangerous rock mass）；
  - 自动驾驶（self-dricing automobile）
  - 天体无人软着陆
三维 Stationary Laser Scanner（静态）
- 非接触式高速激光测量方法，对被测物体快速扫描和测量
- 直接获取激光点所接触物体表面的水平方向、天顶距离、斜距和反射强度，并自动存储、计算和获取点云数据

Ch 2 Visualization and Structuring of Point Clouds

机载点云 Airbone Point Cloud
车载点云 Vehicle-based Point Cloud
地面点云 Tessestrial Point Cloud

2.1 Data Visualization

二维图像：
1. 垂直投影：定义网格；确定每个像素的属性(双线性差值等)；将属性转换为灰度值或颜色
2. 全景排列 panoramic arrangement：球坐标系转换

2.2 Data Structure

2D range images（深度信息）
Delaunay triangulation（德劳内三角网）
- 特征：
  - 凸边
  - 没有交叉边
  - 没有剩余点
- 满足条件：
  - Empty Circle Property 空圆特性
  - Largest minimum angle 最小角：
    
    the minimum angle of each triangle formed by the Delaunay triangulation should be the largest one
    - 德劳内三角网形成的每个三角形的最小角度应为最大角度
- Lawson algorithm to insert point by point：逐点插入的劳森算法
Octrees 八叉树：如果不是空树，则树的每个节点只能包含 0 或 8 个子节点。重建原则：
1. 设置最大递归深度
2. 获取场景的最大尺寸，并以该尺寸构建第一个立方体
3. 将包含的激光点按顺序放入没有子节点的立方体中
4. 如果没有达到最大递归深度，立方体会被细分成八个相等的部分，然后立方体中的所有激光点都会被分配到八个子立方体中
5. 如果分配给子立方体的激光点数量不为零，而且与其母立方体的激光点数量相等，则应结束细分，因为根据空间分割理论，细分后的空间体积应该减小。如果数量保持不变，细分可能会导致无限切割
6. 重复第三步直到达到最大递归深度
K-d tree：k 维二叉树 K-d
- 按照点的离散程度来划分，大的先分（方差），选中间点分
- 时间复杂度 $O(n\log_{2}(n))$

Ch 3 Point Cloud Registration（点云配准）

定义：两个或多个坐标系中的大型三维空间数据点集转换为统一坐标系的数学计算过程。

Point cloud registration is a mathematical calculation process of transforming large 3-D spatial data point sets in two or more coordinate systems into a unified coordinate system.
- 数学配准模型：七参数转换（欧拉角，角坐标系中四元数表示）
- Gimbal Lock 万向节锁：当绕Y轴旋转90度后，X轴和Z轴会重合，指向相同的方向。在这种情况下，物体的两个旋转轴（X轴和Z轴）变得等效于同一个轴。这意味着，系统失去了一个自由度，因为现在只能绕两个独立的轴旋转，而不能绕三个独立的轴旋转。
Corresponding feature recognition 相应的特征识别
- Artificial target: 人工标靶
- Image Feature-based registration：基于图像特征的配准
- Feature-based registration
- Iterative Closest Point （ICP）：精准拼接
Artificial target 人工标靶：
- 手动标靶：Sphere target and Planar target球形目标和平面目标
- 自动标靶
ICP：迭代最近点逼近问题
- 将一个点与相邻点集中的最近点进行匹配
- 根据匹配的对应点计算两个相邻点集之间的变换
- 要求相邻点云之间有很大程度的重叠，对初始值敏感
- 如何确定相应的点：假设最近的点就是对应点，当初始值足够接近时收敛；
  - 2D影像时间复杂度最低，为 $O (n)$
Feature-based registration:registration method based on feature plane matching 基于特征平面的拼接
- 要求至少提取一组相互垂直的功能平面
  - 对点云 P（扫描 1）和 Q（扫描 2）进行分割，以确定平面之间的一一对应关系。选择三对不平行的平面，提取它们的法向量。
  - 求点云 P 中三个平面交点 M 的坐标，点云 Q 中三个相应平面交点 N 的坐标，并分别将点云 P 和 Q 的原点移动到交点 M 和 N 处，使三个平面分别交于原点。
- P 中平面的法向量归一化为 $n_1$ ，Q 中相应平面的法向量归一化为 $n_2$ 。已知两个平面之间的夹角为 $\theta = \cos^{-1}<n_1,n_2>$ ，而 $n_1$ 和 $n_2$ 的垂直轴是旋转轴。根据旋转轴 $A$ 和角度 $\theta$ 可以建立第一个旋转关系。
- 用同样的步骤确定第二和第三个选定平面对之间的旋转关系。将三个旋转关系合并为 P 和 Q 之间的旋转矩阵 R。
- 旋转后，三对平面分别相互平行，然后确定每三个平面的交点，由此计算出平移矢量 T。
- 局限性
  - 3个面中2个面互相平行，只有两个方向，无法解算平移量
- 可以利用成熟的图像匹配方法间接获得点云之间的三维对应点，从而计算出它们之间的配准参数
图像特征匹配
颜色图像
- 先决条件：照相机校准，照相机和扫描仪之间的集成
全景反射图像
- 优点：无需校准摄像机；固有的像素点映射关系
Global registration
- 改正闭合差：对每两站之间的参数进行优化
- 多作一次旋转平移进行平差

Ch 4 Filtering and Classification of Point Clouds

条件
- Filtering：尽可能完整地保留重要的地形特征
- Classification：尽量减少误差，包括第一类误差：剔除属于地形面的激光脚点；第二类误差：接受不属于地形面的激光脚点。
原则
- 高程异常：建筑物引起的高程异常则会在各个方向上形成跳跃边界；高程不连续的轮廓反映了建筑物的形状。
- 判断一个点是否在地表要考虑该点到参考点的距离：随着两点之间距离的增大，阈值也应增大，这主要是为了同时考虑地形起伏引起的高程变化。
方法
- Mathematical Morphology Method 数学形态学方法：移动窗口中最低点视为地形点，与确认点相比，高差在一定范围内的也被视为地形点
- Moving-window filtering method 移动窗口拟合法：使用大尺寸移动窗口找到最低点，计算粗略的地形模型，然后过滤掉所有高差（参照第一步计算的地形模型）超过给定阈值的点，计算出更精确的 DEM。
  - 选取窗口尺寸和预设判断阈值受到场景影响；
  - 滤波参数的设置取决于测量区域的实际地形条件
- Iterative Linear Least Squares Interpolation 迭代线性最小二乘法插值
- 利用残差确定权重，负残差值越大权重越大，残差大于一定范围的视为非地形点
- 进行多次筛选和验证最小值的可靠性
点云分类
- Height texture 高程纹理：The elevation fluctuation formed by the change of local elevation of different objects or different parts of the same object（不同物体或同一物体不同部分的局部高程变化所形成的高程波动）
  - 局部范围内激光脚点之间的高程变化
  - 原始高程数据：分割高处物体
  - 高度差：每个像素周围一定范围内最大高度和最小高度的差值
  - 高度变化：描述一定窗口范围内高度值的变化
  - 地形坡度：根据 X 和 Y 方向的坡度确定每个像素邻域的最大坡度
  - 坡度图像可作为区分倾斜屋顶和水平屋顶、街道和树木的重要标准
  - 缺点：对规则网格进行插值导致插值误差；需要有足够的点密度
- Fusion of laser echo intensity and laser foot points height 融合激光回波强度和激光脚点高度
  - 相邻物体的介质属性差异明显时，通过分析激光回波强度的信息，可以很容易地区分出不同的物体
- Height difference between two echoes of a laser pulse 两次回波的高度差：机载LiDAR
多源数据融合
1. 对激光点云的最后一个回波进行过滤，并将其分为地形点和地面物体点。然后计算第一个和最后一个高程差，作为判断激光点云是否为树木的标准
2. 分析假彩色近红外图像中的植被信息（主要提取图像的红波段信息），以此作为判断是否为植被的线索
3. 从彩色航空图像中提取线条，获得线条特征作为判断建筑物的线索
4. 对彩色正射影像进行分割，将图像划分为许多相对均匀的区域。每个同质区域块都被视为一个对象。使用传统的统计或模糊分类方法来处理这些同质对象

Ch5 Segmentation and Fitting of Point Clouds

5.1 Segmentation of point cloud data

定义：the process of aggregating points according to specific similarity principle. 根据特定的相似性原则聚合点的过程
聚集原则
- 平面化：点到平面的范围在一定的范围
- 特征面平滑度：二次曲面拟合(Quaric surtace)
- 邻域法向量：局部范围内点构成平面计算的法向量平行
Surface Growing 区域增长法：考虑激光点与其空间邻域点之间的关系，并将具有相似属性的激光点集合起来形成一个区域。
- 实施：找到一个种子点作为增长的起点，然后合并具有相同或相似属性的其他点。从未分离的点集中拾取一个新的种子点，重复上述过程，直到每个点都属于一个同质性区域，且相邻区域的同质性不同。
- Seed point selection 选取种子点：
  - 随机选择一个点及其邻近点，形成候选种子点集
  - 平面/曲线表面检测
  - 计算符合特征模型的候选点比例
  - 如果比率大于给定的临界值，则确定种子表面；否则重复
- Growing criterion 生长标准：聚集原则
- Growing ending condition 生长结束条件：没有满足生长条件的点
- 算法步骤：八个方向
  - 建立一个空栈；从点云数据中随机选择未分配给曲面的点
  - 搜索点 a 邻域，构建候选种子点集，并进行有效性判断
  - 将种子表面放入堆栈，并将每个点对应的标记修改为非零值
  - 判断堆栈是否为空。如果堆栈中有元素，则继续下面的步骤。否则，完成生长过程
  - 如果邻域点满足生长条件，它们就会与种子点合并，更新生长模型，然后压入堆栈。如此循环往复，直到邻域点全部被遍历。
  - 在完成点 a 的邻域点生长后，从堆栈中提取另一个点作为种子点，并返回步骤 4
  - 栈空后重复1-6步，直到每个点都被分配
- 分割阈值的选择：过大欠分割(Under Segmentation)，过小过分割(Over Segmentation)，宁愿欠分割
- 邻域搜索：二维图像和 k-d 树
- 存在问题：选择种子点随机性，导致结果不稳定；验证种子点表面运算成本高
- 基于反射影像生长：对原始点云生成的反射图像进行分割，并将拟合的基元面作为种子面进行生长。实现分割和拟合的一体化。

5.2 Fitting of point cloud data

定义：The process of calculating the parameters of the feature model from the coordinates of discrete laser points. 根据离散激光点坐标计算特征模型参数的过程
方法：最小二乘法；三维霍夫变换；RANdom SAmple Consensus (RANSAC) 随机抽样一致性；BAYes SAmple Consensus (BAYSAC)贝叶斯抽样一致性
三维霍夫变换：影像空间到参数空间；可用于检测图像中可参数化的任意形状物体（线、圆等）；根据变换方程，图像中的每个点都可以代表方程中的参数值，参数空间（或累加器数组）的峰值则代表形状对象的参数。
- 如果激光脚点数据集包含平面上的点，则这些点对应的参数空间平面相交于某一点。一般来说，根据投票原则，得票最多的参数集被视为最终的平面参数。
- 只能发现平面，不能准确计算平面的参数
RANdom SAmple Consensus：通过重复选择随机数据子集来实现其目标
- 一个模型适用于假设的正确点，所有未知参数都可以根据假设的正确点计算出来（根据随机的点构建一个参数模型）
- 用 1 中得到的模型测试所有其他数据，如果某个点符合估计模型，则认为该点正确
- 如果有足够多的点被归类为假设正确点，那么估计模型就足够合理
- 用所有假设的正确点对模型进行重新估计，因为只有最初假设的正确点才能对模型进行估计
- 通过估算模型的正确点和错误率对模型进行评估
- 参数：
  - t：用于确定数据是否符合模型的阈值
  - k：算法迭代次数
    $k=\frac{\log (1-p)}{\log (1-w^n)}$
    1-p代表至少有一个点是错的概率，w=正确点个数/数据集的数目
  - d：正确点达到多少算好
- 限制性：当原始点云中的异常值过多时，RANSAC 的计算成本就会大大增加。
- - maybe_inliers：可能的局内点
BAYes SAmple Consensus (BAYSAC)贝叶斯抽样一致性：在RANSAC基础上结合贝叶斯估计原理，并在采样过程中引入内点先验概率的概念。
- 确定内点的先验概率
- 选取数据集中内点概率最高的 n 个点作为样本集
- 初始模型参数取自 n 个点
- 检查所有剩余样本数据，确定内点和异常值
- 每次迭代结束前，根据贝叶斯估计公式更新每个点的内部概率，并开始下一次迭代